text:
"""
    LiftingLine{NA,TA} <: AbstractModel

Lifting line model with `NA` cross sections, using the aerodynamic models in `TA`.
State variables, inputs, and parameters correspond to the state variables, inputs,
and parameters of each of the cross sections concatenated
"""
struct LiftingLine{NA,TA<:NTuple{NA,Any}} <: AbstractModel
    models::TA
end

# --- Constructors --- #

"""
    LiftingLine(models)

Construct a lifting line aerodynamic model given a tuple of 2D aerodynamic models.
"""
LiftingLine(models)

LiftingLine{NA}(models::TA) where TA<:NTuple{NA,Any} where NA = LiftingLine{NA,TA}(models)

"""
    LiftingLine{NA}(model)

Construct a lifting line aerodynamic model using `NA` instances of `model`.
"""
LiftingLine{NA}(model) where NA = LiftingLine(ntuple(i -> model, NA))

# --- Traits --- #

number_of_states(model::LiftingLine) = sum(number_of_states.(model.models))
number_of_inputs(model::LiftingLine) = sum(number_of_inputs.(model.models))
number_of_parameters(model::LiftingLine) = sum(number_of_parameters.(model.models))

inplaceness(::Type{<:LiftingLine}) = InPlace()

function rate_jacobian_type(::Type{LiftingLine{NA,TA}}) where {NA,TA}
    model_types = (TA.parameters...,)
    if all(isempty.(rate_jacobian_type.(model_types)))
        return Empty()
    elseif all(iszero.(rate_jacobian_type.(model_types)))
        return Zeros()
    elseif all(isidentity.(rate_jacobian_type.(model_types)))
        return Identity()
    elseif all(isinvariant.(rate_jacobian_type.(model_types)))
        return Invariant()
    elseif all(isconstant.(rate_jacobian_type.(model_types)))
        return Constant()
    elseif all(islinear.(rate_jacobian_type.(model_types)))
        return Linear()
    else
        return Nonlinear()
    end
end

function state_jacobian_type(::Type{LiftingLine{NA,TA}}) where {NA,TA}
    model_types = (TA.parameters...,)
    if all(isempty.(state_jacobian_type.(model_types)))
        return Empty()
    elseif all(iszero.(state_jacobian_type.(model_types)))
        return Zeros()
    elseif all(isidentity.(state_jacobian_type.(model_types)))
        return Identity()
    elseif all(isinvariant.(state_jacobian_type.(model_types)))
        return Invariant()
    elseif all(isconstant.(state_jacobian_type.(model_types)))
        return Constant()
    elseif all(islinear.(state_jacobian_type.(model_types)))
        return Linear()
    else
        return Nonlinear()
    end
end

function input_jacobian_type(::Type{LiftingLine{NA,TA}}) where {NA,TA}
    model_types = (TA.parameters...,)
    if all(isempty.(input_jacobian_type.(model_types)))
        return Empty()
    elseif all(iszero.(input_jacobian_type.(model_types)))
        return Zeros()
    elseif all(isidentity.(input_jacobian_type.(model_types)))
        return Identity()
    elseif all(isinvariant.(input_jacobian_type.(model_types)))
        return Invariant()
    elseif all(isconstant.(input_jacobian_type.(model_types)))
        return Constant()
    elseif all(islinear.(input_jacobian_type.(model_types)))
        return Linear()
    else
        return Nonlinear()
    end
end

function parameter_jacobian_type(::Type{LiftingLine{NA,TA}}) where {NA,TA}
    model_types = (TA.parameters...,)
    if all(isempty.(parameter_jacobian_type.(model_types)))
        return Empty()
    elseif all(iszero.(parameter_jacobian_type.(model_types)))
        return Zeros()
    elseif all(isidentity.(parameter_jacobian_type.(model_types)))
        return Identity()
    elseif all(isinvariant.(parameter_jacobian_type.(model_types)))
        return Invariant()
    elseif all(isconstant.(parameter_jacobian_type.(model_types)))
        return Constant()
    elseif all(islinear.(parameter_jacobian_type.(model_types)))
        return Linear()
    else
        return Nonlinear()
    end
end

function time_gradient_type(::Type{LiftingLine{NA,TA}}) where {NA,TA}
    model_types = (TA.parameters...,)
    if all(isempty.(time_gradient_type.(model_types)))
        return Empty()
    elseif all(iszero.(time_gradient_type.(model_types)))
        return Zeros()
    elseif all(isinvariant.(time_gradient_type.(model_types)))
        return Invariant()
    elseif all(isconstant.(time_gradient_type.(model_types)))
        return Constant()
    elseif all(islinear.(time_gradient_type.(model_types)))
        return Linear()
    else
        return Nonlinear()
    end
end

# --- Methods --- #

function get_residual!(resid, model::LiftingLine, dx, x, y, p, t)

    models = model.models

    ix = state_indices(models)
    iy = input_indices(models)
    ip = parameter_indices(models)

    resids = view.(Ref(resid), ix)
    dxs = view.(Ref(dx), ix)
    xs = view.(Ref(x), ix)
    ys = view.(Ref(y), iy)
    ps = view.(Ref(p), ip)

    get_residual!.(resids, models, dxs, xs, ys, ps, t)

    return resid
end

# --- Performance 0verloads --- #

function get_rate_jacobian!(J, model::LiftingLine)

    J .= 0

    models = model.models

    ix = state_indices(model.models)

    Js = view.(Ref(J), ix, ix)

    get_rate_jacobian!.(Js, models)

    return J
end

function get_rate_jacobian!(J, model::LiftingLine, p)

    J .= 0

    models = model.models

    ix = state_indices(models)
    ip = parameter_indices(models)

    Js = view.(Ref(J), ix, ix)
    ps = view.(Ref(p), ip)

    get_rate_jacobian!.(Js, models, ps)

    return J
end

function get_rate_jacobian!(J, model::LiftingLine, dx, x, y, p, t)

    J .= 0

    models = model.models

    ix = state_indices(models)
    iy = input_indices(models)
    ip = parameter_indices(models)

    Js = view.(Ref(J), ix, ix)
    dxs = view.(Ref(dx), ix)
    xs = view.(Ref(x), ix)
    ys = view.(Ref(y), iy)
    ps = view.(Ref(p), ip)

    get_rate_jacobian!.(Js, models, dxs, xs, ys, ps, t)

    return J
end

function get_state_jacobian!(J, model::LiftingLine)

    J .= 0

    models = model.models

    ix = state_indices(model.models)

    Js = view.(Ref(J), ix, ix)

    get_state_jacobian!.(Js, models)

    return J
end

function get_state_jacobian!(J, model::LiftingLine, p)

    J .= 0

    models = model.models

    ix = state_indices(models)
    ip = parameter_indices(models)

    Js = view.(Ref(J), ix, ix)
    ps = view.(Ref(p), ip)

    get_state_jacobian!.(Js, models, ps)

    return J
end

function get_state_jacobian!(J, model::LiftingLine, dx, x, y, p, t)

    J .= 0

    models = model.models

    ix = state_indices(models)
    iy = input_indices(models)
    ip = parameter_indices(models)

    Js = view.(Ref(J), ix, ix)
    dxs = view.(Ref(dx), ix)
    xs = view.(Ref(x), ix)
    ys = view.(Ref(y), iy)
    ps = view.(Ref(p), ip)

    get_state_jacobian!.(Js, models, dxs, xs, ys, ps, t)

    return J
end

function get_input_jacobian!(J, model::LiftingLine)

    J .= 0

    models = model.models

    ix = state_indices(models)
    iy = input_indices(models)

    Js = view.(Ref(J), ix, iy)

    get_input_jacobian!.(Js, models)

    return J
end

function get_input_jacobian!(J, model::LiftingLine, p)

    J .= 0

    models = model.models

    ix = state_indices(models)
    iy = input_indices(models)
    ip = parameter_indices(models)

    Js = view.(Ref(J), ix, iy)
    ps = view.(Ref(p), ip)

    get_input_jacobian!.(Js, models, ps)

    return J
end

function get_input_jacobian!(J, model::LiftingLine, dx, x, y, p, t)

    J .= 0

    models = model.models

    ix = state_indices(models)
    iy = input_indices(models)
    ip = parameter_indices(models)

    Js = view.(Ref(J), ix, iy)
    dxs = view.(Ref(dx), ix)
    xs = view.(Ref(x), ix)
    ys = view.(Ref(y), iy)
    ps = view.(Ref(p), ip)

    get_input_jacobian!.(Js, models, dxs, xs, ys, ps, t)

    return J
end

function get_parameter_jacobian!(J, model::LiftingLine)

    J .= 0

    models = model.models

    ix = state_indices(model.models)
    ip = parameter_indices(models)

    Js = view.(Ref(J), ix, ip)

    get_parameter_jacobian!.(Js, models)

    return J
end

function get_parameter_jacobian!(J, model::LiftingLine, p)

    J .= 0

    models = model.models

    ix = state_indices(models)
    ip = parameter_indices(models)

    Js = view.(Ref(J), ix, ip)
    ps = view.(Ref(p), ip)

    get_parameter_jacobian!.(Js, models, ps)

    return J
end

function get_parameter_jacobian!(J, model::LiftingLine, dx, x, y, p, t)

    J .= 0

    models = model.models

    ix = state_indices(models)
    iy = input_indices(models)
    ip = parameter_indices(models)

    Js = view.(Ref(J), ix, ip)
    dxs = view.(Ref(dx), ix)
    xs = view.(Ref(x), ix)
    ys = view.(Ref(y), iy)
    ps = view.(Ref(p), ip)

    get_parameter_jacobian!.(Js, models, dxs, xs, ys, ps, t)

    return J
end

function get_time_gradient!(dT, model::LiftingLine)

    dT .= 0

    models = model.models

    ix = state_indices(model.models)

    dTs = view.(Ref(dT), ix)

    get_time_gradient!.(dTs, models)

    return dT
end

function get_time_gradient!(dT, model::LiftingLine, p)

    dT .= 0

    models = model.models

    ix = state_indices(models)
    ip = parameter_indices(models)

    dTs = view.(Ref(dT), ix)
    ps = view.(Ref(p), ip)

    get_time_gradient!.(dTs, models, ps)

    return dT
end

function get_time_gradient!(dT, model::LiftingLine, dx, x, y, p, t)

    dT .= 0

    models = model.models

    ix = state_indices(models)
    iy = input_indices(models)
    ip = parameter_indices(models)

    dTs = view.(Ref(dT), ix)
    dxs = view.(Ref(dx), ix)
    xs = view.(Ref(x), ix)
    ys = view.(Ref(y), iy)
    ps = view.(Ref(p), ip)

    get_time_gradient!.(dTs, models, dxs, xs, ys, ps, t)

    return dT
end

# --- Convenience Methods --- #

function set_states!(x, model::LiftingLine; section_states)

    section_models = model.models

    section_indices = state_indices(section_models)

    vxs = view.(Ref(x), section_indices)

    bfn! = (x, model, states) -> set_states!(x, model; states...)

    bfn!.(vxs, section_models, section_states)

    return x
end

function set_inputs!(y, model::LiftingLine; section_inputs)

    section_models = model.models

    section_indices = input_indices(section_models)

    vys = view.(Ref(y), section_indices)

    bfn! = (y, model, inputs) -> set_inputs!(y, model; inputs...)

    bfn!.(vys, section_models, section_inputs)

    return y
end

function set_parameters!(p, model::LiftingLine; section_parameters)

    section_models = model.models

    section_indices = parameter_indices(section_models)

    vps = view.(Ref(p), section_indices)

    bfn! = (p, model, parameters) -> set_parameters!(p, model; parameters...)

    bfn!.(vps, section_models, section_parameters)

    return p
end

function separate_states(model::LiftingLine, x)

    section_models = model.models

    section_indices = state_indices(section_models)

    vxs = view.(Ref(x), section_indices)

    return (section_states = separate_states.(section_models, vxs),)
end

function separate_inputs(model::LiftingLine, y)

    section_models = model.models

    section_indices = input_indices(section_models)

    vys = view.(Ref(y), section_indices)

    return (section_inputs = separate_inputs.(section_models, vys),)
end

function separate_parameters(model::LiftingLine, p)

    section_models = model.models

    section_indices = parameter_indices(section_models)

    vps = view.(Ref(p), section_indices)

    return (section_parameters = separate_parameters.(section_models, vps),)
end

# --- Internal Methods --- #


text:
using Documenter, GraffSDK

makedocs(
    modules = [GraffSDK],
    format = :html,
    sitename = "GraffSDK.jl",
    pages = Any[
        "Home" => "index.md",
        "Getting Started" => [
            "Introduction" => "getting_started.md",
            "Making Robust Calls" => "handling_errors.md",
            "Common Structures and Functions" => "ref_common.md",
            "Working with Data" => "working_with_data.md",
            "Caching Data Locally" => "local_caching.md"
        ],
        "Examples" => [
            "Introduction" => "examples/examples.md",
            "Basic Initialization" => "examples/basics_initialization.md",
            "Basic Robot" => "examples/basics_robot.md",
            "Basic Session" => "examples/basics_session.md",
            "Building Graphs" => "examples/basics_variablesandfactors.md"
        ],
        "End-to-End Examples" => [
            "Hexagonal Robot" => "examples/hexagonal.md",
            "Brookstone Rover" => "examples/brookstone.md"
        ],
        "Service Reference" => [
            "User Service" => "ref_user.md",
            "Robot Service" => "ref_robot.md",
            "Session Service" => "ref_session.md"
        ],
        "Reference" => "reference.md"
    ]
    # html_prettyurls = !("local" in ARGS),
    )

deploydocs(
    repo   = "github.com/GearsAD/GraffSDK.jl.git",
    target = "build",
    deps   = Deps.pip("mkdocs", "python-markdown-math"),
    make   = nothing,
    julia  = "0.6",
    osname = "linux"
)

ext: jl


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: jl:light,ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.1
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# + [markdown] run_control={"frozen": false, "read_only": false}
# # Wavefunctions in  Psi4
#
# One very advantageous feature of  Psi4  is the
# ability to generate, return, and manipulate wavefunctions both from computations and as independent entities.
# This is particularly useful because of the depth of information carried by a wavefunction -- which is formally
# an instance of the [`psi4.core.Wavefunction`](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Wavefunction
# "Go to API") class.  This tutorial introduces these objects
# and provides an overview of their capabilities that will be leveraged in future tutorials.
#
# Let's begin our discussion by importing  Psi4  and NumPy, and setting
# some basic options for  Psi4, like the memory, to direct output to a file
# named `output.dat`, and options to be used when performing a computation.

# + run_control={"frozen": false, "read_only": false}
# ==> Basic Options <==
using PyCall: pyimport
psi4 = pyimport("psi4")
np   = pyimport("numpy")

# Memory & Output File
psi4.set_memory(Int(2e9))
psi4.core.set_output_file("output.dat", false)

# Computation options
psi4.set_options(Dict("basis" => "aug-cc-pvdz",
                      "scf_type" => "df",
                      "e_convergence" => 1e-8,
                      "d_convergence" => 1e-8))

# + [markdown] run_control={"frozen": false, "read_only": false}
# Now that we've set the basics, let's use what we learned in the Molecule tutorial to define a water molecule, in Z-matrix
# format, specifying that we want $C_1$ symmetry (instead of letting  Psi4
#  detect the real symmetry $C_{\text{2v}}$):

# + run_control={"frozen": false, "read_only": false}
# ==> Define C_1 Water Molecule <==
h2o = psi4.geometry("""
O
H 1 0.96
H 1 0.96 2 104.5
symmetry c1
""")

# + [markdown] run_control={"frozen": false, "read_only": false}
# Other than molecules, of course, the quintessential object within quantum chemistry (and arguably **the**
# quintessential object in all of quantum mechanics) is the _wavefunction_.  Every method in quantum mechanics seeks to
# find the wavefunction which describes the state of the system of interest.  So, how can we create these objects with
#  Psi4?  If we simply wish to perform a computation (e.g., Hartree–Fock
# or MP2), all we need to do is to define the molecule, call [``psi4.energy()``](http://psicode.org/psi4manual/master/api/psi4.driver.energy.html#psi4.driver.energy "Go to API"), and  Psi4  will do the rest.  What about if we need a wavefunction _before_ performing a computation, or in
# order to implement a method?  Fortunately, the class method [`Wavefunction.build()`](http://psicode.org/psi4manual
# /master/psi4api.html#psi4.core.Wavefunction.build "Go to Documentation") allows us to build one from scratch, given
# a molecule and a basis set.  In the cell below, we've illustrated how to invoke this function:

# + run_control={"frozen": false, "read_only": false}
# ==> Build wavefunction for H2O from scratch with Wavefunction.build() <==
h2o_wfn = psi4.core.Wavefunction.build(h2o, psi4.core.get_global_option("basis"))

# + [markdown] run_control={"frozen": false, "read_only": false}
# Notice that we have passed the variable `h2o` (itself an instance of the `psi4.core.Molecule` class) and the AO basis
# set we wish to use to construct the wavefunction for this molecule.  We could have just as easily passed the string
# `'aug-cc-pvdz'` as an argument, but then we would have to remember to change the argument if we ever changed the
#  Psi4  option in the `psi4.set_options()` block above. Generally, when
# creating something like a wavefunction or a basis  Psi4-side, the class
# instances themselves are what is used to do so.  (Don't worry too much about creating basis sets yet, we'll cover
# these in more detail later.)
#
# Now that we have built an instance of the `Wavefunction` class, we can access our wavefunction's information by
# calling any of the member functions of the `Wavefunction` class on our object.  For instance, the number of spin-up
# ($\alpha$) electrons in our wavefuntion can be found using the [`Wavefunction.nalpha()`](http://psicode.org/psi4manual
# /master/psi4api.html#psi4.core.Wavefunction.nalpha "Go to Documentation") function:
# ~~~julia
# h2o_wfn.nalpha()
# ~~~
# Since the water molecule
# above was defined to be a neutral singlet, we expect that the total number of electrons in our wavefunction should
# be $\alpha + \beta = 2\alpha$.  Let's check:
# -

# Compute the number of electrons in water
println("Water has $(Int(2h2o_wfn.nalpha())) electrons, according to our wavefunction.")

# + [markdown] run_control={"frozen": false, "read_only": false}
# Good,  Psi4  and every General Chemistry textbook on the planet agree.
# What other information can be gleaned from our wavefunction object? For now, not a whole lot.  Since we have built our wavefunction from our molecule and choice of basis set but haven't yet computed anything, the wavefunction doesn't have the orbitals, electron density, energy, or Psi variables attributes set.  Once a computation has been performed,
# however, all this information may be accessed.  This may be accomplised by _returning the wavefunction_ from a
# successful computation, a concept referred to in  Psi4-lingo as
# _wavefunction passing_.  To run a computation with  Psi4, the function
# [`psi4.energy()`](http://psicode.org/psi4manual/master/api/psi4.driver.energy.html#psi4.driver.energy "Go to
# Documentation") is invoked with a particular quantum chemical method, like `'scf'`, `'mp2'`, or `'ccsd(t)`'.  To
# return the wavefunction from that method, the additional argument `return_wfn=True` can be specified:
# ~~~julia
# # Returning a CCSD(T) wavefunction
# energy, wfn = psi4.energy("ccsd(t)", return_wfn=true)
# ~~~
# Then, both the energy and the wavefunction are returned.  Give it a try yourself, for the Hartree–Fock computation
# in the cell below:

# + run_control={"frozen": false, "read_only": false}
# Get the SCF wavefunction & energies for H2O

scf_e, scf_wfn = psi4.energy("scf", return_wfn=true)
println("A float and a Wavefunction object returned: ", scf_e, "\n", scf_wfn)

# + [markdown] run_control={"frozen": false, "read_only": false}
# Now, we can access information you would expect a wavefunction to carry — basically everything we couldn't before.
# Below is summarized several quantities which will be used throughout the modules and tutorials to come.  All these  wavefunction attributes are available after a Hartree–Fock computation; make sure to try them out on our `scf_wfn`!
#
# | Quantity | Function(s) | Description |
# |----------|-------------|-------------|
# | Orbital Coefficients, **C** | [wfn.Ca()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Wavefunction.Ca "Go to Documentation"), [wfn.Cb()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Wavefunction.Cb "Go to Documentation") | Returns orbital coefficient matrix for $\alpha$ (Ca) or $\beta$ (Cb) orbitals. (Identical for restricted orbitals) |
# | Electron Density, **D** | [wfn.Da()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Wavefunction.Da "Go to Documentation"), [wfn.Db()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Wavefunction.Db "Go to Documentation") | One-particle density matrices for $\alpha$ (Da) and $\beta$ (Db) electrons. (Identical for restricted orbitals) |
# | Fock Matrix, **F** | [wfn.Fa()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Wavefunction.Fa "Go to Documentation"), [wfn.Fb()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Wavefunction.Fb "Go to Documentation") | Returns the Fock matrix. For wavefunction with unrestricted orbitals, distinct Fock matrices $F^{\alpha}$ and $F^{\beta}$ for $\alpha$ and $\beta$ orbitals, respectively, are created.|
# | Basis Set | [wfn.basisset()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Wavefunction.basisset "Go to Documentation") | Returns basis set associated with the wavefunction. |
# | $\alpha$ ($\beta$) electrons | [wfn.nalpha()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Wavefunction.nalpha "Go to Documentation"), [wfn.nbeta()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Wavefunction.nbeta "Go to Documentation") | Returns number of $\alpha$ ($\beta$) electrons of the wavefunction. |
# | Irreducible Representations | [wfn.nirrep()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Wavefunction.nirrep "Go to Documentation") | Returns number of irreducible representations (number of symmetry elements). Several objects can utilize molecular symmetry in the wavefunction. |
# | Occupied Orbitals | [wfn.doccpi()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Wavefunction.doccpi "Go to Documentation") | Returns number of doubly occupied orbitals per irrep in the wavefunction. |
# | Psi Variables | [wfn.variables()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Wavefunction.variables "Go to Documentation") | Returns all Psi variables associated with the method which computed the wavefunction. |
# | Energy   | [wfn.energy()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Wavefunction.energy "Go to Documentation") | Returns current energy of the wavefunction. |
# | Orbital Energies, $\boldsymbol{\epsilon}$ | [wfn.epsilon_a()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Wavefunction.epsilon_a "Go to Documentation"), [wfn.epsilon_b()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Wavefunction.epsilon_b "Go to Documentation") | Returns $\alpha$ (a) and $\beta$ (b) orbital energies. (Identical for restricted orbitals) |
#
# Note: The functions returning any of the matrices mentioned above (**C**, **D**, $\boldsymbol{\epsilon}$), actually
# return instances of the `psi4.core.Matrix` class (noticing a pattern here?) and not viewable arrays.  Fortunately,
# the previous tutorial introduced how to modify these arrays Julia-side, using Julia views created through our `psi4view()`.
#
# The full list is quite extensive; however, this likely comprises the most utilized functions. It should be noted that the "a" stands for alpha and conversely the beta quantities can be accessed with the letter "b". For now let's ensure that all computations have C1 symmetry; molecular symmetry can be utilized in Psi4NumPy computations but adds significant complexity to the code.

# + run_control={"frozen": false, "read_only": false}
# Try out the wavefunction class member functions!


text:
#
# Tests
#

using ELM

tests = ["xor.jl"]

println("Running test:")

for tst in tests
	println(" * $(tst)")
	include(tst)
end


text:
# functions to deal with MATLAB arrays

mutable struct MxArray
    ptr::Ptr{Cvoid}
    own::Bool

    function MxArray(p::Ptr{Cvoid}, own::Bool)
        p == C_NULL && error("NULL pointer for MxArray.")
        self = new(p, own)
        if own
            finalizer(release, self)
        end
        return self
    end
end
MxArray(p::Ptr{Cvoid}) = MxArray(p, true)

mxarray(mx::MxArray) = mx

function release(mx::MxArray)
    if mx.own && mx.ptr != C_NULL
        ccall(mx_destroy_array[], Cvoid, (Ptr{Cvoid},), mx.ptr)
    end
    mx.ptr = C_NULL
    return nothing
end

# delete & copy

function delete(mx::MxArray)
    if mx.own
        ccall(mx_destroy_array[], Cvoid, (Ptr{Cvoid},), mx)
    end
    mx.ptr = C_NULL
    return nothing
end

function copy(mx::MxArray)
    pm = ccall(mx_duplicate_array[], Ptr{Cvoid}, (Ptr{Cvoid},), mx)
    return MxArray(pm)
end

function unsafe_convert(::Type{Ptr{Cvoid}}, mx::MxArray)
    ptr = mx.ptr
    ptr == C_NULL && throw(UndefRefError())
    return ptr
end
# functions to create mxArray from Julia values/arrays

const MxRealNum = Union{Float64, Float32, Int32, UInt32, Int64, UInt64, Int16, UInt16, Int8, UInt8, Bool}
const MxComplexNum = Union{ComplexF32, ComplexF64}
const MxNum = Union{MxRealNum, MxComplexNum}

###########################################################
#
#  MATLAB types
#
###########################################################

const mwSize = UInt
const mwIndex = Int

@enum mxClassID::Cint begin
    mxUNKNOWN_CLASS
    mxCELL_CLASS
    mxSTRUCT_CLASS
    mxLOGICAL_CLASS
    mxCHAR_CLASS
    mxVOID_CLASS
    mxDOUBLE_CLASS
    mxSINGLE_CLASS
    mxINT8_CLASS
    mxUINT8_CLASS
    mxINT16_CLASS
    mxUINT16_CLASS
    mxINT32_CLASS
    mxUINT32_CLASS
    mxINT64_CLASS
    mxUINT64_CLASS
    mxFUNCTION_CLASS
    mxOPAQUE_CLASS
    mxOBJECT_CLASS
end

@enum mxComplexity::Cint begin
    mxREAL
    mxCOMPLEX
end

mxclassid(::Type{Bool})    = mxLOGICAL_CLASS
mxclassid(::Union{Type{Float64}, Type{ComplexF64}}) = mxDOUBLE_CLASS
mxclassid(::Union{Type{Float32}, Type{ComplexF32}}) = mxSINGLE_CLASS
mxclassid(::Type{Int8})    = mxINT8_CLASS
mxclassid(::Type{UInt8})   = mxUINT8_CLASS
mxclassid(::Type{Int16})   = mxINT16_CLASS
mxclassid(::Type{UInt16})  = mxUINT16_CLASS
mxclassid(::Type{Int32})   = mxINT32_CLASS
mxclassid(::Type{UInt32})  = mxUINT32_CLASS
mxclassid(::Type{Int64})   = mxINT64_CLASS
mxclassid(::Type{UInt64})  = mxUINT64_CLASS

mxcomplexflag(::Type{T}) where {T<:MxRealNum}    = mxREAL
mxcomplexflag(::Type{T}) where {T<:MxComplexNum} = mxCOMPLEX

const classid_type_map = Dict{mxClassID,Type}(
    mxLOGICAL_CLASS => Bool,
    mxCHAR_CLASS    => Char,
    mxDOUBLE_CLASS  => Float64,
    mxSINGLE_CLASS  => Float32,
    mxINT8_CLASS    => Int8,
    mxUINT8_CLASS   => UInt8,
    mxINT16_CLASS   => Int16,
    mxUINT16_CLASS  => UInt16,
    mxINT32_CLASS   => Int32,
    mxUINT32_CLASS  => UInt32,
    mxINT64_CLASS   => Int64,
    mxUINT64_CLASS  => UInt64
)

function mxclassid_to_type(cid::mxClassID)
    ty = get(classid_type_map, cid, nothing)
    ty === nothing && throw(ArgumentError("The input class id is not a primitive type id."))
    return ty
end

###########################################################
#
#  Functions to access mxArray
#
#  Part of the functions (e.g. mxGetNumberOfDimensions)
#  are actually a macro replacement of an internal
#  function name as (xxxx_730)
#
###########################################################

macro mxget_attr(fun, ret, mx)
    :(ccall($(esc(fun)), $(esc(ret)), (Ptr{Cvoid},), $(esc(mx))))
end

classid(mx::MxArray) = @mxget_attr(mx_get_classid[], mxClassID, mx)
nrows(mx::MxArray)   = convert(Int, @mxget_attr(mx_get_m[], UInt, mx))
ncols(mx::MxArray)   = convert(Int, @mxget_attr(mx_get_n[], UInt, mx))
nelems(mx::MxArray)  = convert(Int, @mxget_attr(mx_get_nelems[], UInt, mx))
ndims(mx::MxArray)   = convert(Int, @mxget_attr(mx_get_ndims[], mwSize, mx))

eltype(mx::MxArray)  = mxclassid_to_type(classid(mx))
elsize(mx::MxArray)  = convert(Int, @mxget_attr(mx_get_elemsize[], UInt, mx))
data_ptr(mx::MxArray) = convert(Ptr{eltype(mx)}, @mxget_attr(mx_get_data[], Ptr{Cvoid}, mx))
real_ptr(mx::MxArray) = convert(Ptr{eltype(mx)}, @mxget_attr(mx_get_pr[], Ptr{Cvoid}, mx))
imag_ptr(mx::MxArray) = convert(Ptr{eltype(mx)}, @mxget_attr(mx_get_pi[], Ptr{Cvoid}, mx))

mxnfields(mx::MxArray) = convert(Int, @mxget_attr(mx_get_nfields[], Cint, mx))

# validation functions

macro mx_test_is(fun, mx)
    :((ccall($(esc(fun)), Bool, (Ptr{Cvoid},), $(esc(mx)))))
end

is_double(mx::MxArray) = @mx_test_is(mx_is_double[], mx)
is_single(mx::MxArray) = @mx_test_is(mx_is_single[], mx)
is_int64(mx::MxArray)  = @mx_test_is(mx_is_int64[], mx)
is_uint64(mx::MxArray) = @mx_test_is(mx_is_uint64[], mx)
is_int32(mx::MxArray)  = @mx_test_is(mx_is_int32[], mx)
is_uint32(mx::MxArray) = @mx_test_is(mx_is_uint32[], mx)
is_int16(mx::MxArray)  = @mx_test_is(mx_is_int16[], mx)
is_uint16(mx::MxArray) = @mx_test_is(mx_is_uint16[], mx)
is_int8(mx::MxArray)   = @mx_test_is(mx_is_int8[], mx)
is_uint8(mx::MxArray)  = @mx_test_is(mx_is_uint8[], mx)

is_numeric(mx::MxArray) = @mx_test_is(mx_is_numeric[], mx)
is_logical(mx::MxArray) = @mx_test_is(mx_is_logical[], mx)
is_complex(mx::MxArray) = @mx_test_is(mx_is_complex[], mx)
is_sparse(mx::MxArray)  = @mx_test_is(mx_is_sparse[], mx)
is_struct(mx::MxArray)  = @mx_test_is(mx_is_struct[], mx)
is_cell(mx::MxArray)    = @mx_test_is(mx_is_cell[], mx)
is_char(mx::MxArray)    = @mx_test_is(mx_is_char[], mx)
is_empty(mx::MxArray)   = @mx_test_is(mx_is_empty[], mx)

# size function

function size(mx::MxArray)
    nd = ndims(mx)
    pdims::Ptr{mwSize} = @mxget_attr(mx_get_dims[], Ptr{mwSize}, mx)
    _dims = unsafe_wrap(Array, pdims, (nd,))
    dims = Vector{Int}(undef, nd)
    for i = 1:nd
        dims[i] = Int(_dims[i])
    end
    tuple(dims...)
end

function size(mx::MxArray, d::Integer)
    d <= 0 && throw(ArgumentError("The dimension must be a positive integer."))

    nd = ndims(mx)
    if nd == 2
        d == 1 ? nrows(mx) :
        d == 2 ? ncols(mx) : 1
    else
        pdims::Ptr{mwSize} = @mxget_attr(mx_get_dims[], Ptr{mwSize}, mx)
        _dims = unsafe_wrap(Array, pdims, (nd,))
        d <= nd ? Int(_dims[d]) : 1
    end
end

###########################################################
#
#  functions to create & delete MATLAB arrays
#
###########################################################

function _dims_to_mwSize(dims::Tuple{Vararg{Integer,N}}) where {N}
    _dims = Vector{mwSize}(undef,N)
    for i = 1:N
        _dims[i] = mwSize(dims[i])
    end
    _dims
end

function mxarray(::Type{T}, dims::Tuple{Vararg{Integer,N}}) where {T<:MxNum,N}
    pm = ccall(mx_create_numeric_array[], Ptr{Cvoid},
        (mwSize, Ptr{mwSize}, mxClassID, mxComplexity),
        N, _dims_to_mwSize(dims), mxclassid(T), mxcomplexflag(T))
    MxArray(pm)
end
mxarray(::Type{T}, dims::Integer...) where {T<:MxNum} = mxarray(T, dims)

# create scalars

function mxarray(x::Float64)
    pm = ccall(mx_create_double_scalar[], Ptr{Cvoid}, (Cdouble,), x)
    MxArray(pm)
end

function mxarray(x::Bool)
    pm = ccall(mx_create_logical_scalar[], Ptr{Cvoid}, (Bool,), x)
    MxArray(pm)
end

function mxarray(x::T) where T<:MxRealNum
    pm = ccall(mx_create_numeric_matrix[], Ptr{Cvoid},
        (mwSize, mwSize, mxClassID, mxComplexity),
        1, 1, mxclassid(T), mxcomplexflag(T))

    pdat = ccall(mx_get_data[], Ptr{T}, (Ptr{Cvoid},), pm)

    unsafe_wrap(Array, pdat, (1,))[1] = x
    MxArray(pm)
end
mxarray(x::T) where {T<:MxComplexNum} = mxarray([x])

# conversion from Julia variables to MATLAB
# Note: the conversion is deep-copy, as there is no way to let
# mxArray use Julia array's memory

function mxarray(a::Array{T}) where T<:MxRealNum
    mx = mxarray(T, size(a))
    ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt), data_ptr(mx), a, length(a)*sizeof(T))
    return mx
end

function mxarray(a::Array{T}) where T<:MxComplexNum
    mx = mxarray(T, size(a))
    na = length(a)
    rdat = unsafe_wrap(Array, real_ptr(mx), na)
    idat = unsafe_wrap(Array, imag_ptr(mx), na)
    for i = 1:na
        rdat[i] = real(a[i])
        idat[i] = imag(a[i])
    end
    mx
end

mxarray(a::BitArray) = mxarray(convert(Array{Bool}, a))
mxarray(a::AbstractRange) = mxarray([a;])

# sparse matrix

function mxsparse(ty::Type{Float64}, m::Integer, n::Integer, nzmax::Integer)
    pm = ccall(mx_create_sparse[], Ptr{Cvoid},
        (mwSize, mwSize, mwSize, mxComplexity), m, n, nzmax, mxREAL)
    MxArray(pm)
end

function mxsparse(ty::Type{ComplexF64}, m::Integer, n::Integer, nzmax::Integer)
    pm = ccall(mx_create_sparse[], Ptr{Cvoid},
        (mwSize, mwSize, mwSize, mxComplexity), m, n, nzmax, mxCOMPLEX)
    MxArray(pm)
end

function mxsparse(ty::Type{Bool}, m::Integer, n::Integer, nzmax::Integer)
    pm = ccall(mx_create_sparse_logical[], Ptr{Cvoid},
        (mwSize, mwSize, mwSize), m, n, nzmax)
    MxArray(pm)
end

function _copy_sparse_mat(a::SparseMatrixCSC{V,I}, ir_p::Ptr{mwIndex}, jc_p::Ptr{mwIndex}, pr_p::Ptr{Float64}, pi_p::Ptr{Float64}) where {V<:ComplexF64,I}
    colptr::Vector{I} = a.colptr
    rinds::Vector{I} = a.rowval
    vr::Vector{Float64} = real(a.nzval)
    vi::Vector{Float64} = imag(a.nzval)
    n::Int = a.n
    nnz::Int = length(vr)

    # Note: ir and jc contain zero-based indices

    ir = unsafe_wrap(Array, ir_p, (nnz,))
    for i = 1:nnz
        ir[i] = rinds[i] - 1
    end

    jc = unsafe_wrap(Array, jc_p, (n+1,))
    for i = 1:n+1
        jc[i] = colptr[i] - 1
    end

    copyto!(unsafe_wrap(Array, pr_p, (nnz,)), vr)
    copyto!(unsafe_wrap(Array, pi_p, (nnz,)), vi)
end

function _copy_sparse_mat(a::SparseMatrixCSC{V,I}, ir_p::Ptr{mwIndex}, jc_p::Ptr{mwIndex}, pr_p::Ptr{V}) where {V,I}
    colptr::Vector{I} = a.colptr
    rinds::Vector{I} = a.rowval
    v::Vector{V} = a.nzval
    n::Int = a.n
    nnz::Int = length(v)

    # Note: ir and jc contain zero-based indices

    ir = unsafe_wrap(Array, ir_p, (nnz,))
    for i = 1:nnz
        ir[i] = rinds[i] - 1
    end

    jc = unsafe_wrap(Array, jc_p, (n+1,))
    for i = 1:n+1
        jc[i] = colptr[i] - 1
    end

    copyto!(unsafe_wrap(Array, pr_p, (nnz,)), v)
end

function mxarray(a::SparseMatrixCSC{V,I}) where {V<:Union{Float64,ComplexF64,Bool},I}
    m::Int = a.m
    n::Int = a.n
    nnz = length(a.nzval)
    @assert nnz == a.colptr[n+1]-1

    mx = mxsparse(V, m, n, nnz)
    ir_p = ccall(mx_get_ir[], Ptr{mwIndex}, (Ptr{Cvoid},), mx)
    jc_p = ccall(mx_get_jc[], Ptr{mwIndex}, (Ptr{Cvoid},), mx)

    if V <: ComplexF64
        pr_p = ccall(mx_get_pr[], Ptr{Float64}, (Ptr{Cvoid},), mx)
        pi_p = ccall(mx_get_pi[], Ptr{Float64}, (Ptr{Cvoid},), mx)
        _copy_sparse_mat(a, ir_p, jc_p, pr_p, pi_p)
    else
        pr_p = ccall(mx_get_pr[], Ptr{V}, (Ptr{Cvoid},), mx)
        _copy_sparse_mat(a, ir_p, jc_p, pr_p)
    end
    return mx
end

# char arrays and string

function mxarray(s::String)
    utf16string = transcode(UInt16, s)
    pm = ccall(mx_create_char_array[], Ptr{Cvoid}, (mwSize, Ptr{mwSize},), 2,
               _dims_to_mwSize((1, length(utf16string))))
    mx = MxArray(pm)
    ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt), data_ptr(mx), utf16string,
          length(utf16string)*sizeof(UInt16))
    return mx
end

# cell arrays

function mxcellarray(dims::Tuple{Vararg{Integer,N}}) where {N}
    pm = ccall(mx_create_cell_array[], Ptr{Cvoid}, (mwSize, Ptr{mwSize}),
        N, _dims_to_mwSize(dims))
    MxArray(pm)
end
mxcellarray(dims::Integer...) = mxcellarray(dims)

function get_cell(mx::MxArray, i::Integer)
    pm = ccall(mx_get_cell[], Ptr{Cvoid}, (Ptr{Cvoid}, mwIndex), mx, i-1)
    MxArray(pm, false)
end

function set_cell(mx::MxArray, i::Integer, v::MxArray)
    v.own = false
    ccall(mx_set_cell[], Cvoid, (Ptr{Cvoid}, mwIndex, Ptr{Cvoid}), mx, i - 1, v)
    return nothing
end

function mxcellarray(a::Array)
    pm = mxcellarray(size(a))
    for i = 1:length(a)
        set_cell(pm, i, mxarray(a[i]))
    end
    return pm
end

mxarray(a::Array) = mxcellarray(a)

# struct arrays

function _fieldname_array(fieldnames::String...)
    n = length(fieldnames)
    a = Vector{Ptr{UInt8}}(undef, n)
    for i = 1:n
        a[i] = unsafe_convert(Ptr{UInt8}, fieldnames[i])
    end
    return a
end

function mxstruct(fns::Vector{String})
    a = _fieldname_array(fns...)
    pm = ccall(mx_create_struct_matrix[], Ptr{Cvoid},
        (mwSize, mwSize, Cint, Ptr{Ptr{UInt8}}), 1, 1, length(a), a)
    MxArray(pm)
end

function mxstruct(fn1::String, fnr::String...)
    a = _fieldname_array(fn1, fnr...)
    pm = ccall(mx_create_struct_matrix[], Ptr{Cvoid},
        (mwSize, mwSize, Cint, Ptr{Ptr{UInt8}}), 1, 1, length(a), a)
    MxArray(pm)
end

function set_field(mx::MxArray, i::Integer, f::String, v::MxArray)
    v.own = false
    ccall(mx_set_field[], Cvoid, (Ptr{Cvoid}, mwIndex, Ptr{UInt8}, Ptr{Cvoid}), mx, i-1, f, v)
    return nothing
end

set_field(mx::MxArray, f::String, v::MxArray) = set_field(mx, 1, f, v)

function get_field(mx::MxArray, i::Integer, f::String)
    pm = ccall(mx_get_field[], Ptr{Cvoid}, (Ptr{Cvoid}, mwIndex, Ptr{UInt8}), mx, i-1, f)
    pm == C_NULL && throw(ArgumentError("Failed to get field."))
    MxArray(pm, false)
end

get_field(mx::MxArray, f::String) = get_field(mx, 1, f)

function get_field(mx::MxArray, i::Integer, fn::Integer)
    pm = ccall(mx_get_field_bynum[], Ptr{Cvoid}, (Ptr{Cvoid}, mwIndex, Cint), mx, i-1, fn-1)
    pm == C_NULL && throw(ArgumentError("Failed to get field."))
    MxArray(pm, false)
end

get_field(mx::MxArray, fn::Integer) = get_field(mx, 1, fn)

function get_fieldname(mx::MxArray, i::Integer)
    p = ccall(mx_get_fieldname[], Ptr{UInt8}, (Ptr{Cvoid}, Cint), mx, i-1)
    unsafe_string(p)
end

const Pairs = Union{Pair, NTuple{2}}

function mxstruct(pairs::Pairs...)
    nf = length(pairs)
    fieldnames = Vector{String}(undef, nf)
    for i = 1:nf
        fn = pairs[i][1]
        fieldnames[i] = string(fn)
    end
    mx = mxstruct(fieldnames)
    for i = 1:nf
        set_field(mx, fieldnames[i], mxarray(pairs[i][2]))
    end
    return mx
end

function mxstruct(d::T) where T
    names = fieldnames(T)
    names_str = map(string, names)
    mx = mxstruct(names_str...)
    for i = 1:length(names)
        set_field(mx, names_str[i], mxarray(getfield(d, names[i])))
    end
    return mx
end

function mxstructarray(d::Array{T}) where T
    names = fieldnames(T)
    names_str = map(string, names)
    a = _fieldname_array(names_str...)

    pm = ccall(mx_create_struct_array[], Ptr{Cvoid}, (mwSize, Ptr{mwSize}, Cint,
        Ptr{Ptr{UInt8}}), ndims(d), _dims_to_mwSize(size(d)), length(a), a)
    mx = MxArray(pm)

    for i = 1:length(d), j = 1:length(names)
        set_field(mx, i, names_str[j], mxarray(getfield(d[i], names[j])))
    end
    return mx
end

mxstruct(d::AbstractDict) = mxstruct(collect(d)...)
mxarray(d) = mxstruct(d)

###########################################################
#
#  convert from MATLAB to Julia
#
###########################################################

# use deep-copy from MATLAB variable to Julia array
# in practice, MATLAB variable often has shorter life-cycle

function _jarrayx(fun::String, mx::MxArray, siz::Tuple)
    if is_numeric(mx) || is_logical(mx)
        @assert !is_sparse(mx)
        T = eltype(mx)
        if is_complex(mx)
            rdat = unsafe_wrap(Array, real_ptr(mx), siz)
            idat = unsafe_wrap(Array, imag_ptr(mx), siz)
            a = complex.(rdat, idat)
        else
            a = Array{T}(undef, siz)
            if !isempty(a)
                ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt), a, data_ptr(mx), length(a)*sizeof(T))
            end
        end
        return a
        #unsafe_wrap(Array, data_ptr(mx), siz)
    elseif is_cell(mx)
        a = Array{Any}(undef, siz)
        for i = 1:length(a)
            a[i] = jvalue(get_cell(mx, i))
        end
        return a
    else
        throw(ArgumentError("$(fun) only applies to numeric, logical or cell arrays."))
    end
end

jarray(mx::MxArray) = _jarrayx("jarray", mx, size(mx))
jvector(mx::MxArray) = _jarrayx("jvector", mx, (nelems(mx),))

function jmatrix(mx::MxArray)
    if ndims(mx) != 2
        throw(ArgumentError("jmatrix only applies to MATLAB arrays with ndims == 2."))
    end
    return _jarrayx("jmatrix", mx, (nrows(mx), ncols(mx)))
end

function jscalar(mx::MxArray)
    if !(nelems(mx) == 1 && (is_logical(mx) || is_numeric(mx)))
        throw(ArgumentError("jscalar only applies to numeric or logical arrays with exactly one element."))
    end
    @assert !is_sparse(mx)
    if is_complex(mx)
        return unsafe_wrap(Array, real_ptr(mx), (1,))[1] + im*unsafe_wrap(Array, imag_ptr(mx), (1,))[1]
    else
        return unsafe_wrap(Array, data_ptr(mx), (1,))[1]
    end
end

function _jsparse(ty::Type{T}, mx::MxArray) where T<:MxRealNum
    m = nrows(mx)
    n = ncols(mx)
    ir_ptr = ccall(mx_get_ir[], Ptr{mwIndex}, (Ptr{Cvoid},), mx)
    jc_ptr = ccall(mx_get_jc[], Ptr{mwIndex}, (Ptr{Cvoid},), mx)

    jc_a::Vector{mwIndex} = unsafe_wrap(Array, jc_ptr, (n+1,))
    nnz = jc_a[n+1]

    ir = Vector{Int}(undef, nnz)
    jc = Vector{Int}(undef, n+1)

    ir_x = unsafe_wrap(Array, ir_ptr, (nnz,))
    for i = 1:nnz
        ir[i] = ir_x[i] + 1
    end

    jc_x = unsafe_wrap(Array, jc_ptr, (n+1,))
    for i = 1:n+1
        jc[i] = jc_x[i] + 1
    end

    pr_ptr = ccall(mx_get_pr[], Ptr{T}, (Ptr{Cvoid},), mx)
    pr::Vector{T} = copy(unsafe_wrap(Array, pr_ptr, (nnz,)))
    if is_complex(mx)
        pi_ptr = ccall(mx_get_pi[], Ptr{T}, (Ptr{Cvoid},), mx)
        pi::Vector{T} = copy(unsafe_wrap(Array, pi_ptr, (nnz,)))
        return SparseMatrixCSC(m, n, jc, ir, pr + im.*pi)
    else
        return SparseMatrixCSC(m, n, jc, ir, pr)
    end
end

function jsparse(mx::MxArray)
    if !is_sparse(mx)
        throw(ArgumentError("jsparse only applies to sparse matrices."))
    end
    return _jsparse(eltype(mx), mx)
end

function String(mx::MxArray)
    if !(classid(mx) == mxCHAR_CLASS && ((ndims(mx) == 2 && nrows(mx) == 1) || is_empty(mx)))
        throw(ArgumentError("String(mx::MxArray) only applies to strings (i.e. char vectors)"))
    end
    return transcode(String, unsafe_wrap(Array, Ptr{UInt16}(data_ptr(mx)), ncols(mx)))
end

function Dict(mx::MxArray)
    if !(is_struct(mx) && nelems(mx) == 1)
        throw(ArgumentError("Dict(mx::MxArray) only applies to a single struct"))
    end
    nf = mxnfields(mx)
    fnames = Vector{String}(undef, nf)
    fvals = Vector{Any}(undef, nf)
    for i = 1:nf
        fnames[i] = get_fieldname(mx, i)
        pv = ccall(mx_get_field_bynum[], Ptr{Cvoid}, (Ptr{Cvoid}, mwIndex, Cint), mx, 0, i-1)
        fx = MxArray(pv, false)
        fvals[i] = jvalue(fx)
    end
    Dict(zip(fnames, fvals))
end

function jvalue(mx::MxArray)
    if is_numeric(mx) || is_logical(mx)
        if !is_sparse(mx)
            nelems(mx) == 1 ? jscalar(mx) :
            ndims(mx) == 2 ? (ncols(mx) == 1 ? jvector(mx) : jmatrix(mx)) :
            jarray(mx)
        else
            jsparse(mx)
        end
    elseif is_char(mx) && (nrows(mx) == 1 || is_empty(mx))
        String(mx)
    elseif is_cell(mx)
        ndims(mx) == 2 ? (ncols(mx) == 1 ? jvector(mx) : jmatrix(mx)) :
        jarray(mx)
    elseif is_struct(mx) && nelems(mx) == 1
        Dict(mx)
    else
        throw(ArgumentError("Unsupported kind of variable."))
    end
end

# deep conversion from MATLAB variable to Julia array

jvalue(mx::MxArray, ::Type{Array})  = jarray(mx)
jvalue(mx::MxArray, ::Type{Vector}) = jvector(mx)
jvalue(mx::MxArray, ::Type{Matrix}) = jmatrix(mx)
jvalue(mx::MxArray, ::Type{Number}) = jscalar(mx)::Number
jvalue(mx::MxArray, ::Type{String}) = String(mx)
jvalue(mx::MxArray, ::Type{Dict}) = Dict(mx)
jvalue(mx::MxArray, ::Type{SparseMatrixCSC}) = jsparse(mx)

# legacy support (eventually drop, when all constructors added)
jdict(mx::MxArray) = Dict(mx)
jstring(mx::MxArray) = String(mx)


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using AdFem
using Test
Random.seed!(233)

function tpfa_op(kvalue,bc,pval,m,n,h)
    tpfa_op_ = load_op_and_grad("./build/libTpfaOp","tpfa_op", multiple=true)
    kvalue,bc,pval,m_,n_,h = convert_to_tensor(Any[kvalue,bc,pval,m,n,h], [Float64,Int64,Float64,Int64,Int64,Float64])
    ii, jj, vv, rhs = tpfa_op_(kvalue,bc,pval,m_,n_,h)
    SparseTensor(ii + 1, jj + 1, vv, m*n, m*n), set_shape(rhs, m*n)
end

# TODO: specify your input parameters
m = 20
n = 10
h = 0.1
kvalue = rand(m*n)
bc = bcedge("all", m, n, h)
pval = rand(size(bc, 1))
K, rhs = compute_fvm_tpfa_matrix(kvalue,bc,pval,m,n,h)
K1, rhs1 = tpfa_op(kvalue,bc,pval,m,n,h)
sess = Session(); init(sess)
K2 = run(sess, K1)
rhs2 = run(sess, rhs1)

@test rhs≈rhs2
@test K ≈ K2

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(x)
    # A, rhs = tpfa_op(kvalue,bc,x,m,n,h)
    # return sum(values(A)^2) + sum(rhs^2)
    A, rhs = tpfa_op(x,bc,pval,m,n,h)
    return sum(values(A)^2) + sum(rhs^2)
end

# TODO: change `m_` and `v_` to appropriate values
# m_ = constant(rand(size(bc,1)))
# v_ = rand(size(bc,1))

m_ = constant(rand(m*n))
v_ = rand(m*n)

y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
using SimpleDiffEq, StaticArrays, OrdinaryDiffEq, Test

function loop(u, p, t)
    @inbounds begin
        σ = p[1]; ρ = p[2]; β = p[3]
        du1 = σ*(u[2]-u[1])
        du2 = u[1]*(ρ-u[3]) - u[2]
        du3 = u[1]*u[2] - β*u[3]
        return SVector{3}(du1, du2, du3)
    end
end
function liip(du, u, p, t)
    σ = p[1]; ρ = p[2]; β = p[3]
    du[1] = σ*(u[2]-u[1])
    du[2] = u[1]*(ρ-u[3]) - u[2]
    du[3] = u[1]*u[2] - β*u[3]
    return nothing
end

u0 = 10ones(3)
dt = 1e-2

odeoop = ODEProblem{false}(loop, SVector{3}(u0), (0.0, 100.0),  [10, 28, 8/3])
sol = solve(odeoop,SimpleATsit5(),dt=dt)
sol2 = solve(odeoop,GPUSimpleATsit5(),dt=dt)

@test sol.u == sol2.u
@test sol.t == sol2.t


text:
using SimpleSVMs
using Test
import Clp
using JuMP
import Random
import Ipopt

@testset "Improving accuracy with relaxed w" begin
    Random.seed!(33)
    X = vcat(rand(20, 2), rand(30,2) .+ [1,0.5]')
    y = append!(ones(20), -ones(30))
    prevloss = 500.0
    for rhs in [0.0, 0.1, 0.5, 2.0, 10.0, 40.0, 100.0]
        (m, w, b) = SimpleSVMs.build_svm(SimpleSVMs.L1Penalty(rhs), X, y, with_optimizer(Clp.Optimizer, LogLevel = 0))
        optimize!(m)
        res = JuMP.objective_value(m)
        @test res < prevloss
        prevloss = res
        if res ≈ 0.0
            break
        end
    end
    @test prevloss ≈ 0
end

@testset "Non-separable domains" begin
    Random.seed!(33)
    X = vcat(randn(20, 2), randn(30,2) .+ [3.0,1.5]')
    y = append!(ones(20), -ones(30))
    prevloss = 500.0
    for rhs in [0.0, 0.1, 0.5, 1.5, 100.0]
        (m, w, b) = SimpleSVMs.build_svm(SimpleSVMs.L1Penalty(rhs), X, y, with_optimizer(Clp.Optimizer, LogLevel = 0))
        optimize!(m)
        res = JuMP.objective_value(m)
        @test res ≤ prevloss
        prevloss = res
    end
    @test prevloss > 0
end

@testset "L2-restriction" begin
    Random.seed!(33)
    X = vcat(randn(20, 2), randn(30,2) .+ [3.0,1.5]')
    y = append!(ones(20), -ones(30))
    prevloss = 500.0
    for rhs in [0.0, 0.1, 0.5, 1.5, 100.0]
        (m, w, b) = SimpleSVMs.build_svm(SimpleSVMs.L2Penalty(rhs), X, y, with_optimizer(Ipopt.Optimizer, print_level = 0))
        optimize!(m)
        res = JuMP.objective_value(m)
        @test res ≤ prevloss
        prevloss = res
    end
    @test prevloss > 0
end


text:
module Commands

using Compat: addenv
@static if VERSION >= v"1.6"
    using Preferences: @load_preference, @set_preferences!
end

export mpiexec

@static if VERSION >= v"1.6"
    get_path() = @load_preference("mpiexec path", "mpiexec")
    function set_path(path::String)
        @assert ispath(path)
        @set_preferences!("mpiexec path" => path)
    end
else
    get_path() = "mpiexec"
end

"Represent the configurations of a command."
abstract type CommandConfig end

# See https://www.open-mpi.org/doc/v3.0/man1/mpiexec.1.php & https://www.mpich.org/static/docs/v3.1/www1/mpiexec.html
const SHORT_OPTIONS = (
    # Open MPI options
    "h",
    "q",
    "v",
    "v",
    "n",
    "h",
    "c",
    "n",
    "np",
    "rf",
    "s",
    "wd",
    "wdir",
    "x",
    "am",
    "cf",
    "d",
    # MPICH options
    "f",
    "l",
)
const LONG_OPTIONS = (
    # Open MPI options
    "help",
    "quiet",
    "verbose",
    "version",
    "display-map",
    "display-allocation",
    "output-proctable",
    "dvm",
    "max-vm-size",
    "novm",
    "hnp",
    "host",
    "hostfile",
    "default-hostfile",
    "machinefile",
    "cpu-set",
    "npersocket",
    "npernode",
    "pernode",
    "map-by",
    "bycore",
    "byslot",
    "nolocal",
    "nooversubscribe",
    "oversubscribe",
    "bynode",
    "cpu-list",
    "rank-by",
    "bind-to",
    "cpus-per-proc",
    "cpus-per-rank",
    "bind-to-core",
    "bind-to-socket",
    "report-bindings",
    "rankfile",
    "output-filename",
    "stdin",
    "merge-stderr-to-stdout",
    "tag-output",
    "timestamp-output",
    "xml",
    "xml-file",
    "xterm",
    "path",
    "prefix",
    "noprefix",
    "preload-binary",
    "preload-files",
    "set-cwd-to-session-dir",
    "gmca",
    "mca",
    "tune",
    "debug",
    "get-stack-traces",
    "debugger",
    "timeout",
    "tv",
    "allow-run-as-root",
    "app",
    "cartofile",
    "continuous",
    "disable-recovery",
    "do-not-launch",
    "do-not-resolve",
    "enable-recovery",
    "index-argv-by-rank",
    "leave-session-attached",
    "max-restarts",
    "ompi-server",
    "personality",
    "ppr",
    "report-child-jobs-separately",
    "report-events",
    "report-pid",
    "report-uri",
    "show-progress",
    "terminate",
    "use-hwthread-cpus",
    "use-regexp",
    "debug-devel",
    "debug-daemons",
    "debug-daemons-file",
    "display-devel-allocation",
    "display-devel-map",
    "display-diffable-map",
    "display-topo",
    "launch-agent",
    "report-state-on-timeout",
    # MPICH options
    "arch",
    "file",
    "soft",
    "configfile",
)

"""
    mpiexec(env = Pair{String,String}[]; kwargs...)

Generate a function from `kwargs` and `env`.
"""
function mpiexec(env = Pair{String,String}[]; kwargs...)
    args = [get_path()]
    for (arg, val) in kwargs
        if arg in (:env, :genv, :envlist, :genvlist)
            throw(ArgumentError("Please treat `$arg` as a positional argument `env`."))
        end
        _pusharg!(args, string(arg), val)
    end
    return function (exec)
        append!(args, exec)
        cmd = Cmd(args)
        return addenv(cmd, env...)
    end
end

function _pusharg!(args, arg, val)
    arg = replace(arg, '_' => '-')
    option = (arg in LONG_OPTIONS ? "--" : '-') * arg
    if val isa AbstractVector{<:AbstractString}
        join(val, ',')
        return push!(args, option, val)
    elseif val isa Bool  # flag
        return push!(args, option)
    elseif val isa Pair
        return push!(args, option, string(val.first), string(val.second))
    elseif val isa AbstractVector{<:Pair} || val isa AbstractDict
        for v in val
            push!(args, option, string(v.first), string(v.second))
        end
    else
        return push!(args, option, string(val))
    end
end

end


text:
"""
    PartialRange(start, stop)
"""
struct PartialRange
    start::Number
    stop::Number
end

(r::PartialRange)(length) = Base.range(r.start, r.stop; length)
(r::PartialRange)(; length=nothing, step=nothing) = Base.range(r.start, r.stop; length, step)

Base.step(r::PartialRange, s) = Base.range(r.start, r.stop; step=s)
Base.length(r::PartialRange, l) = Base.range(r.start, r.stop; length=l)


text:
using CompileBot

snoop_bench(
    BotConfig(
        "Plots",
    ),
    joinpath(@__DIR__, "precompile_script.jl"),
)


text:
module Tasks
export TaskParams, PredictorTask

using FactorizationMachines: FMMatrix, FMFloat
using FactorizationMachines.Common: nlogsig, nlogsig_deriv, sqerr, sqerr_deriv

abstract TaskParams
abstract PredictorTask

"""Represents a classification task"""
immutable ClassificationTaskParams <: TaskParams
end
const classification = ClassificationTaskParams

"""Classification parameters derived from data"""
immutable ClassificationTask <: PredictorTask
end

"""
Given data `X` and `y`, initializes a `ClassificationTask`
"""
function init(::ClassificationTaskParams, X::FMMatrix, y::Vector{FMFloat})
  ClassificationTask()
end

loss(::ClassificationTask, p::Number, y::Number)       = nlogsig(p, y)
loss_deriv(::ClassificationTask, p::Number, y::Number) = nlogsig_deriv(p, y)

"""Represents a regression task"""
immutable RegressionTaskParams <: TaskParams
end
const regression = RegressionTaskParams

"""Regression parameters derived from data"""
immutable RegressionTask <: PredictorTask
  target_min::FMFloat
  target_max::FMFloat

  RegressionTask(; target_min::FMFloat = typemin(FMFloat),
                   target_max::FMFloat = typemax(FMFloat)) = new(target_min, target_max)
end

"""
Given data `X` and `y`, initializes a `RegressionTask`
"""
function init(::RegressionTaskParams, X::FMMatrix, y::Vector{FMFloat})
  RegressionTask(target_min = minimum(y), target_max = maximum(y))
end

loss(::RegressionTask, p::Number, y::Number)       = sqerr(p, y)
loss_deriv(::RegressionTask, p::Number, y::Number) = sqerr_deriv(p, y)

end


text:
using Statistics
using Test

using MarketData

using TimeSeries

@testset "basemisc" begin

@testset "cumulative functions" begin
    let ta = cumsum(cl)
        @test values(ta) == cumsum(values(cl))
        @test meta(ta) == meta(cl)
    end

    let ta = cumsum(ohlc, dims = 1)
        @test values(ta) == cumsum(values(ohlc), dims = 1)
        @test meta(ta) == meta(ohlc)
    end

    let ta = cumsum(cl, dims = 2)
        @test values(ta) == cumsum(values(cl), dims = 2)
        @test meta(ta) == meta(cl)
    end

    let ta = cumsum(ohlc, dims = 2)
        @test values(ta) == cumsum(values(ohlc), dims = 2)
        @test meta(ta) == meta(ohlc)
    end

    @test_throws DimensionMismatch cumsum(cl, dims = 3)
    @test_throws DimensionMismatch cumsum(ohlc, dims = 3)

    let ta = cumprod(cl[1:5])
        @test values(ta) == cumprod(values(cl[1:5]))
        @test meta(ta) == meta(cl[1:5])
    end

    let ta = cumprod(ohlc[1:5], dims = 1)
        @test values(ta) == cumprod(values(ohlc[1:5]), dims = 1)
        @test meta(ta) == meta(ohlc[1:5])
    end

    let ta = cumprod(cl[1:5], dims = 2)
        @test values(ta) == cumprod(values(cl[1:5]), dims = 2)
        @test meta(ta) == meta(cl[1:5])
    end

    let ta = cumprod(ohlc[1:5], dims = 2)
        @test values(ta) == cumprod(values(ohlc[1:5]), dims = 2)
        @test meta(ta) == meta(ohlc[1:5])
    end

    @test_throws DimensionMismatch cumprod(cl, dims = 3)
    @test_throws DimensionMismatch cumprod(ohlc, dims = 3)
end

@testset "reduction functions" begin
    for (fname, f) ∈ ([(:sum, sum), (:mean, mean)])
        for (name, src) ∈ [(:cl, cl), (:ohlc, ohlc)]
            @testset "$fname::$name" begin
                let ta = f(src)
                    @test meta(ta)   == meta(src)
                    @test length(ta) == 1
                    @test values(ta) == f(values(src), dims = 1)
                end

                let ta = f(src, dims = 2)
                    @test meta(ta)     == meta(src)
                    @test length(ta)   == length(timestamp(src))
                    @test values(ta)   == f(values(src), dims = 2)
                    @test colnames(ta) == [fname]
                end

                @test_throws DimensionMismatch f(src, dims = 3)
            end  # @testset
        end
    end  # for func

    for (fname, f) ∈ ([(:std, std), (:var, var)])
        for (name, src) ∈ [(:cl, cl), (:ohlc, ohlc)]
            @testset "$fname::$name" begin
                let ta = f(src)
                    @test meta(ta)   == meta(src)
                    @test length(ta) == 1
                    @test values(ta) == f(values(src), dims = 1)
                end

                let ta = f(src, dims = 2, corrected = false)
                    @test meta(ta)     == meta(src)
                    @test length(ta)   == length(timestamp(src))
                    @test values(ta)   == f(values(src), dims = 2, corrected = false)
                    @test colnames(ta) == [fname]
                end

                @test_throws DimensionMismatch f(src, dims = 3)
            end  # @testset
        end
    end  # for func
end

@testset "Base.any function" begin
    let  # single column
        """
        julia> any(x .== lag(x))
        1x1 TimeSeries.TimeArray{Bool,1,Date,BitArray{1}} 2000-01-07 to 2000-01-07
        │            │ _     │
        ├────────────┼───────┤
        │ 2000-01-07 │ false │
        """
        let
            ts = timestamp(cl)[1:5]
            ta = TimeArray(ts, 1:5)
            xs = any(ta .== lag(ta))
            @test timestamp(xs)[1] == ts[end]
            @test values(xs)[1]    == false
        end

        """
        julia> any(x .== lag(x), 2)
        4x1 TimeSeries.TimeArray{Bool,1,Date,BitArray{1}} 2000-01-04 to 2000-01-07
        │            │ any   │
        ├────────────┼───────┤
        │ 2000-01-04 │ false │
        │ 2000-01-05 │ true  │
        │ 2000-01-06 │ false │
        │ 2000-01-07 │ false │
        """
        let
            ts = timestamp(cl)[1:5]
            ta = TimeArray(ts, [1, 2, 2, 3, 4])
            xs = any(ta .== lag(ta), dims = 2)

            @test timestamp(xs)   == ts[2:end]
            @test any(values(xs) .== [false, true, false, false])
            @test colnames(xs)    == [:any]
        end
    end  # single column

    let  # multi column
        """
        julia> x .> 3
        3x2 TimeSeries.TimeArray{Bool,2,Date,BitArray{2}} 2000-01-03 to 2000-01-05
        │            │ A     │ B     │
        ├────────────┼───────┼───────┤
        │ 2000-01-03 │ false │ false │
        │ 2000-01-04 │ false │ true  │
        │ 2000-01-05 │ true  │ true  │

        julia> any(x .> 3, 2)
        3x1 TimeSeries.TimeArray{Bool,2,Date,BitArray{2}} 2000-01-03 to 2000-01-05
        │            │ any   │
        ├────────────┼───────┤
        │ 2000-01-03 │ false │
        │ 2000-01-04 │ true  │
        │ 2000-01-05 │ true  │
        """
        ts = timestamp(cl)[1:3]
        ta = TimeArray(ts, [1 2; 3 4; 5 6])
        xs = any(ta .> 3, dims = 2)

        @test timestamp(xs) == ts
        @test any(values(xs) .== [false, true, true])
        @test colnames(xs) == [:any]
    end  # multi column
end  # @testset "Base.any" function

@testset "Base.all function" begin
    let  # single column
        """
        julia> all(x .== lag(x))
        1x1 TimeSeries.TimeArray{Bool,1,Date,BitArray{1}} 2000-01-07 to 2000-01-07
        │            │ A     │
        ├────────────┼───────┤
        │ 2000-01-07 │ false │
        """
        let
            ts = timestamp(cl)[1:5]
            ta = TimeArray(ts, 1:5)
            xs = all(ta .== lag(ta))
            @test timestamp(xs)[] == ts[end]
            @test values(xs)[]    == false
        end

        """
        julia> all(x .== lag(x), 2)
        4x1 TimeSeries.TimeArray{Bool,1,Date,BitArray{1}} 2000-01-04 to 2000-01-07
        │            │ any   │
        ├────────────┼───────┤
        │ 2000-01-04 │ false │
        │ 2000-01-05 │ true  │
        │ 2000-01-06 │ false │
        │ 2000-01-07 │ false │
        """
        let
            ts = timestamp(cl)[1:5]
            ta = TimeArray(ts, [1, 2, 2, 3, 4])
            xs = all(ta .== lag(ta), dims = 2)

            @test timestamp(xs) == ts[2:end]
            @test all(values(xs) .== [false, true, false, false])
            @test colnames(xs) == [:all]
        end
    end  # single column

    let  # multi column
        """
        julia> x .> 3
        3x2 TimeSeries.TimeArray{Bool,2,Date,BitArray{2}} 2000-01-03 to 2000-01-05
        │            │ A     │ B     │
        ├────────────┼───────┼───────┤
        │ 2000-01-03 │ false │ false │
        │ 2000-01-04 │ false │ true  │
        │ 2000-01-05 │ true  │ true  │

        julia> all(x .> 3, 2)
        3x1 TimeSeries.TimeArray{Bool,2,Date,BitArray{2}} 2000-01-03 to 2000-01-05
        │            │ all   │
        ├────────────┼───────┤
        │ 2000-01-03 │ false │
        │ 2000-01-04 │ false │
        │ 2000-01-05 │ true  │
        """
        ts = timestamp(cl)[1:3]
        ta = TimeArray(ts, [1 2; 3 4; 5 6])
        xs = all(ta .> 3, dims = 2)

        @test timestamp(xs) == ts
        @test all(values(xs) .== [false, false, true])
        @test colnames(xs) == [:all]
    end  # multi column
end  # @testset "Base.all" function

end  # @testset "basemisc"


text:
import Pkg
Pkg.add("Plots")
Pkg.add("PyPlot")
using Plots

using Random
Random.seed!(2018)

x = randn(1000)
y = randn(1000)
z = randn(1000)

histogram(x, bins=20, alpha=0.4, label="A")
histogram!(y, bins=20, alpha=0.6, label="B")
histogram!(z, bins=20, alpha=0.8, label="C")

#Makie
using CairoMakie

data = randn(1000)

hist(data, normalization = :pdf, bar_labels = :values,
     label_formatter=x-> round(x, digits=2), label_size = 15,
     strokewidth = 0.5, strokecolor = (:black, 0.5), color = :values)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
#=
Copyright (c) 2016 Saurav Sachidanand

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
=#

module CBOR

include("constants.jl")
include("encoding-common.jl")
include("decoding-common.jl")

export encode
export decode, decode_with_iana
export Simple, Null, Undefined

function decode(cbor_bytes::Array{UInt8, 1})
    data, _ = decode_next(1, cbor_bytes, false)
    return data
end

function decode_with_iana(cbor_bytes::Array{UInt8, 1})
    warn("Results from decode_with_iana may change in the future.")
    data, _ = decode_next(1, cbor_bytes, true)
    return data
end

function encode(data)
    bytes = UInt8[]
    encode(data, bytes)
    return bytes
end

end


text:
module ParallelBarnesHut

# Imports

using Plots

using Random

# Exports

export Particle

export simulate!, grav_acc

export show_particles, show_quadview, show_boxes, animate_frames

export rand_particles

# Structs

"""
Represents a particle in a Barnes Hut simulation. For a galaxy simulation, a
particle is a star or a group of stars.
"""
mutable struct Particle
    "Position coordinates of a particle (in meters)."
    pos::Array{Float64,1}

    "Velocity vector of a particle (in meters per second)."
    vel::Array{Float64,1}

    "Mass of a particle (in kilograms)."
    mass::Float64
end

"""
Used in the octree for the Barnes Hut algorithm. Represents a cubical octant of
space and stores the mass and center of mass of the contained particles.
"""
struct Node
    "Side length of the cubical region of space represented by a node."
    size::Float64

    "Total mass of the particles contained in the region represented by a node."
    total_mass::Float64

    "Center of mass of the particles contained in the region represented by a node."
    center_of_mass::Array{Float64,1}

    "Children of a node in the octree."
    children::Array{Union{Node, Nothing},1}
end

# Barnes Hut Tree Generation

"""
    generate_tree(particles; return_boxes = false)

Generate the Barnes Hut octree for a given array of particles. Optionally
return coordinates of the boxes represented by the nodes in the octree.
"""
function generate_tree(particles::Array{Particle,1}; return_boxes::Bool = false)
    # Calculate bounds of a box enclosing all particles
    minx, maxx = minimum(p -> p.pos[1], particles), maximum(p -> p.pos[1], particles)
    miny, maxy = minimum(p -> p.pos[2], particles), maximum(p -> p.pos[2], particles)
    minz, maxz = minimum(p -> p.pos[3], particles), maximum(p -> p.pos[3], particles)

    # Calculate corner and size of a cube enclosing all particles
    size = max(maxx - minx, maxy - miny, maxz - minz)
    corner = [minx, miny, minz]

    if return_boxes
        # Recursively compute tree of particles and return coordinates of boxes
        boxes = Tuple{Float64,Array{Float64,1}}[]
        boxes_threaded = [Tuple{Float64,Array{Float64,1}}[] for i = 1:Threads.nthreads()]
        generate_tree_helper(particles, size, corner, return_boxes, boxes_threaded)
        for thread = eachindex(boxes_threaded)
           append!(boxes, boxes_threaded[thread])
        end
        boxes
    else
        # Recursively compute tree of particles and return that tree
        generate_tree_helper(particles, size, corner, false, Array{Tuple{Float64,Array{Float64,1}},1}[])
    end
end

"""
    generate_tree_helper(particles, size, corner, save_boxes, boxes_threaded)

Generate subtree of Barnes Hut octree given a list of particles and the corner
size of a cubical region of space. Optionally save coordinates of this region in
the array `boxes_threaded`. Helper function for `generate_tree`.
"""
function generate_tree_helper(particles::Array{Particle,1}, size::Float64, corner::Array{Float64,1},
                              save_boxes::Bool, boxes_threaded::Array{Array{Tuple{Float64,Array{Float64,1}},1},1})
    # If saving boxes, save current box coordinates
    if save_boxes
        push!(boxes_threaded[Threads.threadid()], (size, corner))
    end

    if length(particles) == 0
        # Base case 1: if no particles in octant, return empty leaf node
        Node(size/2, 0, corner .+ size/2, [nothing for i = 1:8])
    elseif length(particles) == 1
        # Base case 2: if one particle in octant, return leaf node with that particle
        Node(size/2, particles[1].mass, particles[1].pos, [nothing for i = 1:8])
    else
        # Recursive case

        # Calculate total mass and center of mass
        total_mass = sum(particle.mass for particle in particles)
        center_of_mass = sum(particle.mass .* particle.pos for particle in particles) / total_mass

        # Sort particles into 8 octants (using threading)
        octants = [Particle[] for oct = 1:8]
        octants_threaded = [[Particle[] for oct = 1:8] for thread = 1:Threads.nthreads()]
        new_corners = vec([corner + [i, j, k] for i = (0, size/2), j = (0, size/2), k = (0, size/2)])
        @sync for particle in particles
            Threads.@spawn begin
                oct_num = 0
                for oct = 1:8
                    if all(new_corners[oct] .<= particle.pos .<= new_corners[oct] .+ size/2)
                        oct_num = oct
                    end
                end
                if oct_num != 0
                    push!(octants_threaded[Threads.threadid()][oct_num], particle)
                end
            end
        end
        for thread = eachindex(octants_threaded)
            for oct = 1:8
               append!(octants[oct], octants_threaded[thread][oct])
            end
        end

        # Calculate child nodes recursively (using threading)
        children::Array{Union{Node, Nothing},1} = fill(nothing, 8)
        @sync for oct = 1:8
            Threads.@spawn children[oct] = generate_tree_helper(octants[oct], size/2, new_corners[oct],
                                                                save_boxes, boxes_threaded)
        end

        # Return node
        Node(size, total_mass, center_of_mass, children)
    end
end

# Barnes Hut Simulation

"""
    simulate!(particles, steps; Δt = 0.1, θ = 0.5, acc_func = grav_acc)

Approximately solve N-body problem for an array of particles for number of time
steps into the future using the Barnes Hut method. Return an array of particle
arrays, where each particle array is the state of all particles at each time
step.

# Arguments
- `particles::Array{Particle,1}`: an array of particles.
- `steps::Int64`: the number of time steps to run the simulation.
- `Δt::Float64 = 0.1`: length of time step in seconds (keyword argument).
- `θ::Float64 = 0.5`: node size to particle distance threshold for Barnes Hut
  simulation (keyword argument).
- `acc_func(mass::Float64, r::Array{Float64,1}) = grav_acc`: calculates particle
  acceleration from mass and distance vector.
"""
function simulate!(particles::Array{Particle,1}, steps::Int64;
                   Δt::Float64 = 0.1, θ::Float64 = 0.5, acc_func = grav_acc)
    # Initialize array of snapshots of particles at each time step
    frames = Array{Particle,1}[]

    # Store initial conditions as first frame
    push!(frames, deepcopy(particles))

    # Advance velocity by half time step to do leapfrog integration method
    vel_step!(particles, generate_tree(particles), Δt/2, θ, acc_func)

    # Simulate steps and save each frame
    for i = 1:steps
        step!(particles, generate_tree(particles), Δt, θ, acc_func)
        push!(frames, deepcopy(particles))
    end

    # Return frames
    frames
end

"""
    vel_step!(particles, tree, Δt, θ, acc_func)

Calculate the acceleration of each particle using `net_acc` and approximate new
velocity using `Δt`.
"""
function vel_step!(particles::Array{Particle,1}, tree::Node, Δt::Float64, θ::Float64, acc_func)
    # Calculate acceleration and approximately advance velocity (using threading)
    @sync for particle in particles
        Threads.@spawn begin
            particle.vel += net_acc(particle, tree, θ, acc_func) * Δt
        end
    end
end

"""
    step!(particles, tree, Δt, θ, acc_func)

Calculate the acceleration of each particle using `net_acc` and approximate new
velocity and position using `Δt`.
"""
function step!(particles::Array{Particle,1}, tree::Node, Δt::Float64, θ::Float64, acc_func)
    # Calculate acceleration, and approximately advance velocity and position (using threading)
    @sync for particle in particles
        Threads.@spawn begin
            particle.vel += net_acc(particle, tree, θ, acc_func) * Δt
            particle.pos += particle.vel * Δt
        end
    end
end

"""
    net_acc(particle, node, θ, acc_func)

Recursively calculate the net acceleration on a particle given the Barnes Hut
octree, a threshold `θ`, and the acceleration function `acc_func`.
"""
function net_acc(particle::Particle, node::Union{Node, Nothing}, θ::Float64, acc_func)
    # Recusively calculates net acceleration on a particle given a particle tree
    if node == nothing
        # Base case 1: if tree is empty, return 0 acceleration
        zeros(Float64, 3)
    else
        # Recursive case, Barnes Hut method
        s = node.size                            # Width of the cube that this node represents
        r = node.center_of_mass .- particle.pos  # Vector from particle to node center of mass
        if s/√sum(r.^2) < θ
            # Base case 2: if node size divided by distance is less than threshold θ,
            #              calculate gravitational acceleration
            acc_func(node.total_mass, r)
        else
            # Recursive case
            sum(net_acc(particle, child, θ, acc_func) for child in node.children)
        end
    end
end

"""
    grav_acc(mass, r; ϵ = 0.02)

Calculate the gravitation acceleration on a particle given the mass of another
object and the distance vector from particle to object. The usual gravitational
force is softened by a factor `ϵ` to prevent it from blowing up when the
distance is very small.
"""
function grav_acc(mass::Float64, r::Array{Float64,1}; ϵ::Float64 = 0.02)
    # Calculate gravitational acceleration, using node center of mass
    # and softening factor ϵ (to prevent blow up at d2 = 0)
    G::Float64 = 6.67430e-11
    ((G * mass) / ((sum(r.^2)+ϵ^2)^(3/2))) .* r
end

# Display and Animate Particles

"""
    show_particles(particles, θ = 30.0, ϕ = 30.0)

Using the `Plots` package, display an array of particles in 3D with azimuthal
camera angle θ and elevation angle ϕ.
"""
function show_particles(particles::Array{Particle,1}, θ::Float64 = 30.0, ϕ::Float64 = 30.0)
    scatter([Tuple(p.pos) for p = particles],
        lims = (-1, 1),
        camera = (θ, ϕ),
        size = (500, 500),
        label = "",
        background_color = :black,
        marker = (:circle, 3, 0.8, :white, stroke(0)))
end

"""
    show_quadview(particles)

Display four views of the particles in 3D at different camera angles.
"""
function show_quadview(particles::Array{Particle,1})
    scatter([Tuple(p.pos) for p = particles],
        lims = (-1, 1),
        camera = (30, 30),
        size = (800, 800),
        label = "",
        background_color = :black,
        marker = (:circle, 3, 0.8, :white, stroke(0)),
        layout = 4,
        subplot = 1)
    scatter!([Tuple(p.pos) for p = particles],
        lims = (-1, 1),
        camera = (0, 90),
        label = "",
        marker = (:circle, 3, 0.8, :white, stroke(0)),
        subplot = 2)
    scatter!([Tuple(p.pos) for p = particles],
        lims = (-1, 1),
        camera = (90, 0),
        label = "",
        marker = (:circle, 3, 0.8, :white, stroke(0)),
        subplot = 3)
    scatter!([Tuple(p.pos) for p = particles],
        lims = (-1, 1),
        camera = (60, 30),
        label = "",
        marker = (:circle, 3, 0.8, :white, stroke(0)),
        subplot = 4)
end

"""
    show_boxes(particles)

Calculate and display in 3D the cubes represented by the Barnes Hut octree.
"""
function show_boxes(particles::Array{Particle,1})
    p = scatter([Tuple(p.pos) for p = particles],
        camera = (30.0, 30.0),
        size = (500, 500),
        label = "",
        background_color = :black,
        marker = (:circle, 3, 0.8, :white, stroke(0)))
    boxes = generate_tree(particles, return_boxes = true)
    for (s, (x, y, z)) in boxes
        plot!([(x,y,z),(x+s,y,z),(x+s,y+s,z),(x,y+s,z),(x,y,z),(x,y,z+s),
               (x+s,y,z+s),(x+s,y+s,z+s),(x,y+s,z+s),(x,y,z+s),(x+s,y,z+s),
               (x+s,y,z),(x+s,y+s,z),(x+s,y+s,z+s),(x,y+s,z+s),(x,y+s,z)],
            linecolor = :red,
            label = "")
    end
    p
end

"""
    animate_frames(frames, frame_rate = 1, file_name = "animation"; θ = 30.0, ϕ = 30.0, quadview = false)

Plot particles for every `frame_rate` frames resulting from a Barnes Hut
simulation, and generate a gif file of this animation with a custom name.
Optionally, specify azimuthal camera angle θ and elevation angle ϕ, or set
quadview to `true` in order to animate at four different angles.
"""
function animate_frames(frames::Array{Array{Particle,1},1}, frame_rate::Int64 = 1, file_name = "animation";
                        θ::Float64 = 30.0, ϕ::Float64 = 30.0, quadview::Bool = false)
    anim = Animation()
    for i = 1:frame_rate:length(frames)
        if quadview
            show_quadview(frames[i])
        else
            show_particles(frames[i], θ, ϕ)
        end
        frame(anim)
    end
    gif(anim, "$file_name.gif")
end

# Generate Random Particles

"""
    rand_particles(num_particles::Int64)

Generate `num_particles` random particles using a distribution similar to a
disk-shaped galaxy.
"""
function rand_particles(num_particles::Int64)
    particles = Particle[]
    for i = 1:num_particles
        θ = 2π*rand()
        R = randexp()/ℯ
        z = (rand() > 0.5 ? 1 : -1) * (randexp()/(10ℯ))
        v = √(R*num_particles/1e6)
        mass = 1e4
        push!(particles, Particle([R*cos(θ), R*sin(θ), z], [v*sin(θ), -v*cos(θ), 0.0], mass))
    end
    particles
end

end # module


text:
println("Hello world!");

@everywhere using Printf;

@show workdir = if length(ARGS) > 0
  ARGS[1];
else
  ".";
end

@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x));
@everywhere fmt_int(x) = @sprintf("%03d", x);

@everywhere function prefix(case)
  "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))_run-$(fmt_int(case[:run]))";
end

cases = Any[];
#E0s = [1.0; 5.0];
#K1s = [1e-1; 1.0];
#kTs = [1.0; 0.1];
#Fzs = [0.0; 0.5; 1.0];
#Fxs = [0.0; 0.5; 1.0];
#ns = Int[100];
#bs = [1.0];
runs = 1:25;
for run in runs
  push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0,
                    :kT => 1, :Fz => 0,
                    :Fx => 0, :n => 100, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 5, :K1 => 1, :K2 => 0.0,
                    :kT => 1, :Fz => 0,
                    :Fx => 0, :n => 100, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 5, :K1 => 0.1, :K2 => 0.0,
                    :kT => 1, :Fz => 0,
                    :Fx => 0, :n => 100, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 0.1, :K2 => 0.0,
                    :kT => 1, :Fz => 0,
                    :Fx => 0, :n => 100, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0,
                    :kT => 1, :Fz => 1,
                    :Fx => 0, :n => 100, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0,
                    :kT => 1, :Fz => 2,
                    :Fx => 0, :n => 100, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0,
                    :kT => 1, :Fz => 0,
                    :Fx => 1, :n => 100, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0,
                    :kT => 0.1, :Fz => 0,
                    :Fx => 1, :n => 100, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0,
                    :kT => 0.1, :Fz => 0,
                    :Fx => 0, :n => 100, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0,
                    :kT => 0.1, :Fz => 1,
                    :Fx => 0, :n => 100, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 0.1, :K2 => 0.0,
                    :kT => 0.1, :Fz => 1,
                    :Fx => 0, :n => 100, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0,
                    :kT => 1, :Fz => 0,
                    :Fx => 0, :n => 200, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0,
                    :kT => 1, :Fz => 0,
                    :Fx => 0, :n => 100, :b => 0.5, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0,
                    :kT => 1, :Fz => 0,
                    :Fx => 0, :n => 100, :b => 2.0, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0,
                    :kT => 1, :Fz => 1,
                    :Fx => 0, :n => 100, :b => 0.5, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0,
                    :kT => 1, :Fz => 0,
                    :Fx => 1, :n => 100, :b => 0.5, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0,
                    :kT => 1, :Fz => 1,
                    :Fx => 1, :n => 100, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 3, :K1 => 1, :K2 => 0.0,
                    :kT => 1, :Fz => 3,
                    :Fx => 3, :n => 100, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0,
                    :kT => 1, :Fz => 3,
                    :Fx => 0, :n => 100, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0,
                    :kT => 1, :Fz => 5,
                    :Fx => 0, :n => 100, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0,
                    :kT => 1, :Fz => 0,
                    :Fx => 3, :n => 100, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0,
                    :kT => 1, :Fz => 0,
                    :Fx => 5, :n => 100, :b => 1, :run => run));
  push!(cases, Dict(:E0 => 5, :K1 => 1, :K2 => 0.0,
                    :kT => 1, :Fz => 5,
                    :Fx => 5, :n => 100, :b => 1, :run => run));
end

@info "total number of cases to run: $(length(cases))";

mkpath(workdir);
mkpath(joinpath(workdir, "standard"));
mkpath(joinpath(workdir, "umbrella"));
mkpath(joinpath(workdir, "clustering"));

pmap(case -> begin;

  outfile = joinpath(workdir, "standard", "$(prefix(case)).out");
  if !isfile(outfile)
    println("Running case: $case.");
    command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "standard", prefix(case)))`;
    output = read(command, String);
    write(outfile, output);
  else
    println("Case: $case has already been run.");
  end

  outfile = joinpath(workdir, "umbrella", "$(prefix(case)).out");
  if !isfile(outfile)
    println("Running case: $case.");
    command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type Ising --umbrella-sampling -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "umbrella", prefix(case)))`;
    output = read(command, String);
    write(outfile, output);
  else
    println("Case: $case has already been run.");
  end

  outfile = joinpath(workdir, "clustering", "$(prefix(case)).out");
  if !isfile(outfile)
    println("Running case: $case.");
    command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "clustering", prefix(case)))`;
    output = read(command, String);
    write(outfile, output);
  else
    println("Case: $case has already been run.");
  end

end, cases);

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats