text:
# VA Algorithm 10.2, pg 362
"""
    maxdet_completion(A::SparseMatrixCSC{Tv, Ti}) where {Tv <: AbstractFloat, Ti <: Integer}

Returns the maximum determinant completion of chordal sparse matrix `A` using
Algorithm 10.2 in [VA15].

# Reference
[VA15] [Chordal Graphs and Semidefinite Optimization](https://www.seas.ucla.edu/~vandenbe/publications/chordalsdp.pdf)
by Lieven Vandenberghe and Martin Andersen
"""
function maxdet_completion(A::SparseMatrixCSC{Tv, Ti}) where {Tv <: AbstractFloat, Ti <: Integer}
    !issymmetric(A) && error(ArgumentError("A must be symmetric"))
    n = size(A, 1)

    sp = sparsity_pattern(A)
    L = get_chordal_extension(sp; perm=nothing, verbose=false)[3]
    ct = CliqueTree(L)
    order_snds!(ct)

    W = Matrix(A[ct.perm, ct.perm])
    # cache = zeros(maximum(length.(ct.seps)), maximum(length.(ct.snds)))

    for j in (length(ct.snds)-1):-1:1
        vrep_ind = ct.postordering[j]
        vrep = ct.vreps[vrep_ind]

        ν = ct.snds[vrep_ind]
        α = ct.seps[vrep_ind]
        η = filter(x->(!(x in ν) && !(x in α)), vrep+1:n)

        # cache = W_αα^† * W_αν
        cache = zeros(length(α), length(ν))
        try
            cache = W[α, α] \ W[α, ν]
        catch e
            # Occasionally get a singular exception. Ignoring for now
            cache = pinv(W[α, α]) * W[α, ν]
        end

        # res = W_ηα * cache = W_ηα * W_αα^† * W_αν
        @views mul!(W[η, ν], W[η, α], cache)
        W[ν, η] .= W[η, ν]'
    end
    iperm = invperm(ct.perm)
    return W[iperm, iperm]
end

# Logarithmic Barriers for Sparse Matrix Cones
# Andersen, Dahl, Vandenberghe
# Algorithm 4.2
# NOTE: This is very inefficient -- needs some significant work
# TODO: Clean up
"""
    maxdet_completion_etree(A::SparseMatrixCSC{Tv, Ti}) where {Tv <: AbstractFloat, Ti <: Integer}

Returns the cholesky factors of the inverse of the maximum determinant completion
of chordal sparse matrix `A` using Algorithm 4.2 in [ADV14]. This algorithm uses
the elimination tree of `A` and, therefore, BLAS level 2 operations.

# Reference
[ADV12] [Logarithmic barriers for sparse matrix cones](https://arxiv.org/abs/1203.2742)
by Martin S. Andersen, Joachim Dahl, Lieven Vandenberghe
"""
function maxdet_completion_etree(A::SparseMatrixCSC{Tv, Ti}) where {Tv <: AbstractFloat, Ti <: Integer}
    !issymmetric(A) && error(ArgumentError("A must be symmetric"))
    n = size(A, 1)

    sp = sparsity_pattern(A)
    L = get_chordal_extension(sp; perm=nothing, verbose=false)[3]
    et = EliminationTree(L; peo=collect(1:n))
    postordering = get_postordering(et.par, et.child)

    L_chol = spzeros(n,n)
    L_chol[diagind(L_chol)] .= 1
    D_chol = Diagonal(zeros(n))
    V = Vector{Matrix{Float64}}(undef, n)

    for i in n:-1:1
        j = postordering[i]
        Ij = rowvals(L)[nzrange(L, j)]

        Vj = (i == n) ? [A[end, end]] : V[j]
        if i != n
            try
                L_chol[Ij, j] = (-Vj) \ Vector{Float64}(A[Ij, j])
            catch e
                L_chol[Ij, j] = -pinv(Vj) * Vector{Float64}(A[Ij, j])
            end
            D_chol[j,j] = 1/(A[j,j] + dot(A[Ij, j], L_chol[Ij, j]))
        else
            D_chol[j,j] = 1 / A[j,j]
        end

        for ch in et.child[j]
            Ich = rowvals(L)[nzrange(L, ch)]

            nv = length(Ij) + 1
            tmp = zeros(nv, nv)
            tmp[1,1] = A[j,j]
            tmp[2:end, 1] .= A[Ij, j]
            tmp[1, 2:end] .= A[Ij, j]
            tmp[2:end, 2:end] .= Vj

            E_Jj_Ich = spzeros(length(Ij) + 1, length(Ich))
            Jj = vcat([j], Ij)
            for ii in 1:length(Jj), jj in 1:length(Ich)
                if Jj[ii] == Ich[jj]
                    E_Jj_Ich[ii, jj] = 1.0
                end
            end

            V[ch] = E_Jj_Ich' * tmp * E_Jj_Ich
        end
    end

    return L_chol, D_chol
end

# Logarithmic Barriers for Sparse Matrix Cones
# Andersen, Dahl, Vandenberghe
# Algorithm 7.3
# NOTE: This is very inefficient -- needs some significant work
# TODO: cleanup
"""
    maxdet_completion_etree(A::SparseMatrixCSC{Tv, Ti}) where {Tv <: AbstractFloat, Ti <: Integer}

Returns the cholesky factors of the inverse of the maximum determinant completion
of chordal sparse matrix `A` using Algorithm 7.3 in [ADV14]. This algorithm uses
the supernodal elimination tree of `A` and, therefore, BLAS level 3 operations.

# Reference
[ADV12] [Logarithmic barriers for sparse matrix cones](https://arxiv.org/abs/1203.2742)
by Martin S. Andersen, Joachim Dahl, Lieven Vandenberghe
"""
function maxdet_completion_factors(A::SparseMatrixCSC{Tv, Ti}) where {Tv <: AbstractFloat, Ti <: Integer}
    !issymmetric(A) && error(ArgumentError("A must be symmetric"))
    n = size(A, 1)

    sp = sparsity_pattern(A)
    L = get_chordal_extension(sp; perm=nothing, verbose=false)[3]
    etree_par = get_etree(L)
    vreps, snd_par, snd_membership = max_supernode_etree(L, etree_par)
    n_snds = length(vreps)
    snds = [Vector{Int}(undef, 0) for _ in 1:n_snds]
    for v in 1:n
        push!(snds[snd_membership[v]], v)
    end
    snd_children = get_children_from_par(snd_par)
    post_ord = get_postordering(snd_par, snd_children)

    L_chol = spzeros(n,n)
    L_chol[diagind(L_chol)] .= 1
    D_chol = spzeros(n,n)
    V = Vector{Matrix{Float64}}(undef, n_snds)

    for i in n_snds:-1:1
        vrep_ind = post_ord[i]
        j = vreps[vrep_ind]
        Jj = vcat([j], rowvals(L)[nzrange(L, j)])
        Nj = snds[vrep_ind]
        Aj = filter(x->!(x in Nj), Jj)

        Vj = (i == n_snds) ? [A[end, end]] : V[vrep_ind]
        if i != n_snds
            try
                L_chol[Aj, Nj] = -Vj \ A[Aj, Nj]
            catch e
                L_chol[Aj, Nj] = pinv(-Vj) * A[Aj, Nj]
            end
            D_chol[Nj,Nj] = inv(Matrix(A[Nj,Nj] + A[Aj, Nj]'*L_chol[Aj, Nj]))
        else
            D_chol[Nj,Nj] = inv(Matrix(A[Nj,Nj]))
        end

        for ch_ind in snd_children[vrep_ind]
            Nch = snds[ch_ind]
            ch = vreps[ch_ind]
            Jch = vcat([ch], rowvals(L)[nzrange(L, ch)])
            Ach = filter(x->!(x in Nch), Jch)

            nv = length(Jj)
            tmp = zeros(nv, nv)
            len_Nj = length(Nj)
            len_Aj = length(Aj)
            tmp[1:len_Nj,1:len_Nj] = A[Nj,Nj]
            if len_Aj > 0
                tmp[len_Nj+1:end, 1:len_Nj] .= A[Aj, Nj]
                tmp[1:len_Nj, len_Nj+1:end] .= A[Aj, Nj]'
                tmp[len_Nj+1:end, len_Nj+1:end] .= Vj
            end
            E_Jj_Ach = spzeros(length(Jj), length(Ach))
            for ii in 1:length(Jj), jj in 1:length(Ach)
                if Jj[ii] == Ach[jj]
                    E_Jj_Ach[ii, jj] = 1.0
                end
            end
            V[ch_ind] = E_Jj_Ach' * tmp * E_Jj_Ach
        end
    end

    return L_chol, D_chol
end

"""
    minrank_completion(A::SparseMatrixCSC{Tv, Ti}) where {Tv <: AbstractFloat, Ti <: Integer}

Returns the cholesky factor of the minimum rank completion (`A_complete = Y*Yᵀ`)
of chordal sparse matrix `A` using Algorithm 2 in [Sun15]. This algorithm uses
the supernodal elimination tree of `A` and, therefore, BLAS level 3 operations.

# Reference
[Sun15] [Decomposition Methods for Semidefinite Optimization (PhD Thesis)](https://escholarship.org/content/qt1cv6981p/qt1cv6981p.pdf)
by Yifan Sun
"""
function minrank_completion(A::SparseMatrixCSC{Tv, Ti}) where {Tv <: AbstractFloat, Ti <: Integer}
    !issymmetric(A) && error(ArgumentError("A must be symmetric"))
    n = size(A, 1)

    sp = sparsity_pattern(A)
    L = get_chordal_extension(sp; perm=nothing, verbose=false)[3]
    ct = CliqueTree(L)
    n_snds = length(ct.vreps)

    # Determine rank
    r = maximum([rank(Matrix(A[c,c]), rtol=1e-10) for c in get_cliques(L)])

    Y = zeros(n, r)
    for j in n_snds:-1:1
        vrep_ind = ct.postordering[j]
        vrep = ct.vreps[vrep_ind]

        # col_j = vcat([vrep], rowvals(L)[nzrange(L, vrep)])
        ν = ct.snds[vrep_ind]
        α = filter(x->!(x in ν), rowvals(L)[nzrange(L, vrep)])
        col_j = vcat(ν, α)

        dd, VV = eigen(Symmetric(Matrix(@view(A[col_j,col_j])), :L), sortby=x->-real(x))
        r_ = min(length(dd), r)
        Z = VV[:,1:r_]*Diagonal(sqrt.(max.(real.(dd[1:r_]), 0.0)))

        # First block to be filled
        if j == n_snds
            Y[ν, 1:r_] .= Z[1:length(ν), :]
            continue
        end

        # Zero padding to make dimensions line up appropriately
        U = zeros(length(ν), r)
        V = zeros(size(Z,1) - length(ν), r)
        U[:, 1:r_] .= Z[1:length(ν), 1:r_]
        V[:, 1:r_] .= Z[length(ν)+1:end, 1:r_]

        # Find orthogonal transform s.t. Y[α, :] == V*Q
        Q2, Σ, Q1 = svd(V'*Y[α, :], full=true)
        Q = Q2*Q1'

        if !≈(Y[α, :], V*Q)
            # Take generalized SVD
            F = svd(Matrix(V'), Matrix(Y[α, :]'))
            Q = F.U*F.V'
        end

        Y[ν, :] .= U*Q
    end
    return Y
end


text:
"""
NCDatasets is a module to read and write NetCDF files.
This is minimal example how to read a NetCDF file.

```julia
# open the file and show its metadata (if called in the REPL without ending semicolon)
ds = Dataset("filename.nc","r")
# load all data of the variable temperature
v = ds["temperature"][:,:]
# load the attribute units
unit = v.attrib["units"]
# close the file
close(ds)
```

More information is available at https://github.com/Alexander-Barth/NCDatasets.jl .
"""
module NCDatasets
using Dates
using Printf

using Base
using Compat
using DataStructures: OrderedDict
import Base.convert
import Compat: @debug, findfirst

import Base: close
import Base: Array

using CFTime

include("CatArrays.jl")
export CatArrays

# NetCDFError, error check and netcdf_c.jl from NetCDF.jl (https://github.com/JuliaGeo/NetCDF.jl)
# Copyright (c) 2012-2013: Fabian Gans, Max-Planck-Institut fuer Biogeochemie, Jena, Germany
# MIT

# Exception type for error thrown by the NetCDF library
mutable struct NetCDFError <: Exception
    code::Cint
    msg::String
end

"""
    NetCDFError(code::Cint)

Construct a NetCDFError from the error code.
"""
NetCDFError(code::Cint) = NetCDFError(code, nc_strerror(code))

#function Base.showerror(io::IO, err::NetCDFError)
#    print(io, "NetCDF error code $(err.code):\n\t$(err.msg)")
#end

"Check the NetCDF error code, raising an error if nonzero"
function check(code::Cint)
    # zero means success, return
    if code == Cint(0)
        return
        # otherwise throw an error message
    else
        throw(NetCDFError(code))
    end
end

include("netcdf_c.jl")

# end of code from NetCDF.jl

const default_timeunits = "days since 1900-00-00 00:00:00"

### type definition

# -----------------------------------------------------
# base type of attribytes list
# concrete types are Attributes (single NetCDF file) and MFAttributes (multiple NetCDF files)

abstract type BaseAttributes
end

abstract type AbstractDataset
end

abstract type AbstractVariable{T,N} <: AbstractArray{T,N}
end

abstract type AbstractDimensions
end

abstract type AbstractGroups
end

function Base.get(a::BaseAttributes,name::AbstractString,default)
    if haskey(a,name)
        return a[name]
    else
        return default
    end
end

# -----------------------------------------------------
# List of attributes (for a single NetCDF file)
# all ids should be Cint

mutable struct Attributes <: BaseAttributes
    ncid::Cint
    varid::Cint
    isdefmode::Vector{Bool}
end

mutable struct Groups <: AbstractGroups
    ncid::Cint
    isdefmode::Vector{Bool}
end

mutable struct Dimensions <: AbstractDimensions
    ncid::Cint
    isdefmode::Vector{Bool}
end

mutable struct Dataset <: AbstractDataset
    ncid::Cint
    # true of the NetCDF is in define mode (i.e. metadata can be added, but not data)
    # need to be an array, so that it is copied by reference
    isdefmode::Vector{Bool}
    attrib::Attributes
    dim::Dimensions
    group::Groups
end

# Mapping between NetCDF types and Julia types
const jlType = Dict(
    NC_BYTE   => Int8,
    NC_UBYTE  => UInt8,
    NC_SHORT  => Int16,
    NC_USHORT => UInt16,
    NC_INT    => Int32,
    NC_UINT   => UInt32,
    NC_INT64  => Int64,
    NC_UINT64 => UInt64,
    NC_FLOAT  => Float32,
    NC_DOUBLE => Float64,
    NC_CHAR   => Char,
    NC_STRING => String
)

# Inverse mapping
const ncType = Dict(value => key for (key, value) in jlType)

# default fill value per types
const ncFillValue = Dict(
    Int8    => NC_FILL_BYTE,
    UInt8   => NC_FILL_UBYTE,
    Int16   => NC_FILL_SHORT,
    UInt16  => NC_FILL_USHORT,
    Int32   => NC_FILL_INT,
    UInt32  => NC_FILL_UINT,
    Int64   => NC_FILL_INT64,
    UInt64  => NC_FILL_UINT64,
    Float32 => NC_FILL_FLOAT,
    Float64 => NC_FILL_DOUBLE,
    Char    => NC_FILL_CHAR,
    String  => NC_FILL_STRING
)

"Return all variable names"
listVar(ncid) = String[nc_inq_varname(ncid,varid)
                       for varid in nc_inq_varids(ncid)]

"Return all attribute names"
function listAtt(ncid,varid)
    natts = nc_inq_varnatts(ncid,varid)
    names = Vector{String}(undef,natts)

    for attnum = 0:natts-1
        names[attnum+1] = nc_inq_attname(ncid,varid,attnum)
    end

    return names
end

"Make sure that a dataset is in data mode"
function datamode(ncid,isdefmode::Vector{Bool})
    if isdefmode[1]
        nc_enddef(ncid)
        isdefmode[1] = false
    end
end

"Make sure that a dataset is in define mode"
function defmode(ncid,isdefmode::Vector{Bool})
    if !isdefmode[1]
        nc_redef(ncid)
        isdefmode[1] = true
    end
end

function Base.show(io::IO,a::BaseAttributes; indent = "  ")
    try
        # use the same order of attributes than in the NetCDF file
        for (attname,attval) in a
            print(io,indent,@sprintf("%-20s = ",attname))
            printstyled(io, @sprintf("%s",attval),color=:blue)
            print(io,"\n")
        end
    catch err
        if isa(err,NetCDFError)
            if err.code == NC_EBADID
                print(io,"NetCDF attributes (file closed)")
                return
            end
        end
        rethrow
    end
end

### Dimensions

function Base.keys(d::Dimensions)
    return String[nc_inq_dimname(d.ncid,dimid)
                  for dimid in nc_inq_dimids(d.ncid,false)]
end

function Base.getindex(a::Dimensions,name::AbstractString)
    dimid = nc_inq_dimid(a.ncid,name)
    return nc_inq_dimlen(a.ncid,dimid)
end

"""
    unlimited(d::Dimensions)

Return the names of all unlimited dimensions.
"""
function unlimited(d::Dimensions)
    return String[nc_inq_dimname(d.ncid,dimid)
                  for dimid in nc_inq_unlimdims(d.ncid)]
end

export unlimited

"""
    Base.setindex!(d::Dimensions,len,name::AbstractString)

Defines the dimension called `name` to the length `len`.
Generally dimension are defined by indexing, for example:

```julia
ds = Dataset("file.nc","c")
ds.dim["longitude"] = 100
```

If `len` is the special value `Inf`, then the dimension is considered as
`unlimited`, i.e. it will grow as data is added to the NetCDF file.
"""
function Base.setindex!(d::Dimensions,len,name::AbstractString)
    defmode(d.ncid,d.isdefmode) # make sure that the file is in define mode
    dimid = nc_def_dim(d.ncid,name,(isinf(len) ? NC_UNLIMITED : len))
    return len
end

# Attributes

"""
    getindex(a::Attributes,name::AbstractString)

Return the value of the attribute called `name` from the
attribute list `a`. Generally the attributes are loaded by
indexing, for example:

```julia
ds = Dataset("file.nc")
title = ds.attrib["title"]
```
"""
function Base.getindex(a::Attributes,name::AbstractString)
    return nc_get_att(a.ncid,a.varid,name)
end

"""
    Base.setindex!(a::Attributes,data,name::AbstractString)

Set the attribute called `name` to the value `data` in the
attribute list `a`. Generally the attributes are defined by
indexing, for example:

```julia
ds = Dataset("file.nc","c")
ds.attrib["title"] = "my title"
```
"""
function Base.setindex!(a::Attributes,data,name::AbstractString)
    defmode(a.ncid,a.isdefmode) # make sure that the file is in define mode
    return nc_put_att(a.ncid,a.varid,name,data)
end

"""
    Base.keys(a::Attributes)

Return a list of the names of all attributes.
"""
Base.keys(a::Attributes) = listAtt(a.ncid,a.varid)

# -----------------------------------------------------

mutable struct MFAttributes{T} <: BaseAttributes where T <: BaseAttributes
    as::Vector{T}
end

function Base.getindex(a::MFAttributes,name::AbstractString)
    return a.as[1][name]
end

function Base.setindex!(a::MFAttributes,data,name::AbstractString)
    for a in a.as
        a[name] = data
    end
    return data
end

Base.keys(a::MFAttributes) = keys(a.as[1])

# -----------------------------------------------------

struct Resource
    filename::String
    mode::String
    metadata::OrderedDict
end

mutable struct DeferAttributes <: BaseAttributes
    r::Resource
    varname::String # "/" for global attributes
    data::OrderedDict
end

mutable struct DeferDimensions <: AbstractDimensions
    r::Resource
    data::OrderedDict
end

mutable struct DeferGroups <: AbstractGroups
    r::Resource
    data::OrderedDict
end

# -----------------------------------------------------

mutable struct MFDimensions{T} <: AbstractDimensions where T <: AbstractDimensions
    as::Vector{T}
    aggdim::String
end

mutable struct MFGroups{T} <: AbstractGroups where T <: AbstractGroups
    as::Vector{T}
    aggdim::String
end

function Base.getindex(a::MFDimensions,name::AbstractString)
    if name == a.aggdim
        return sum(d[name] for d in a.as)
    else
        return a.as[1][name]
    end
end

function Base.setindex!(a::MFDimensions,data,name::AbstractString)
    for a in a.as
        a[name] = data
    end
    return data
end

Base.keys(a::Union{MFDimensions,MFGroups}) = keys(a.as[1])

function Base.getindex(a::MFGroups,name::AbstractString)
    ds = getindex.(a.as,name)
    attrib = MFAttributes([d.attrib for d in ds])
    dim = MFDimensions([d.dim for d in ds],a.aggdim)
    group = MFGroups([d.group for d in ds],a.aggdim)

    return MFDataset(ds,a.aggdim,attrib,dim,group)
end

#---
mutable struct MFDataset{T,N,TA,TD,TG} <: AbstractDataset where T <: AbstractDataset
    ds::Array{T,N}
    aggdim::AbstractString
    attrib::MFAttributes{TA}
    dim::MFDimensions{TD}
    group::MFGroups{TG}
end

mutable struct DeferDataset <: AbstractDataset
    r::Resource
    groupname::String
    attrib::DeferAttributes
    dim::DeferDimensions
    group::DeferGroups
    data::OrderedDict
end

"""
    Base.keys(g::NCDatasets.Groups)

Return the names of all subgroubs of the group `g`.
"""
function Base.keys(g::Groups)
    return String[nc_inq_grpname(ncid)
                  for ncid in nc_inq_grps(g.ncid)]
end

"""
    group = getindex(g::NCDatasets.Groups,groupname::AbstractString)

Return the NetCDF `group` with the name `groupname`.
For example:

```julia-repl
julia> ds = Dataset("results.nc", "r");
julia> forecast_group = ds.group["forecast"]
julia> forecast_temp = forecast_group["temperature"]
```

"""
function Base.getindex(g::Groups,groupname::AbstractString)
    grp_ncid = nc_inq_grp_ncid(g.ncid,groupname)
    return Dataset(grp_ncid,g.isdefmode)
end

"""
    defGroup(ds::Dataset,groupname, attrib = []))

Create the group with the name `groupname` in the dataset `ds`.
`attrib` is a list of attribute name and attribute value pairs (see `Dataset`).
"""
function defGroup(ds::Dataset,groupname; attrib = [])
    grp_ncid = nc_def_grp(ds.ncid,groupname)
    ds = Dataset(grp_ncid,ds.isdefmode)

    # set global attributes for group
    for (attname,attval) in attrib
        ds.attrib[attname] = attval
    end

    return ds
end

group(ds::AbstractDataset,groupname) = ds.group[groupname]

# -----------------------------------------------------
# Dataset

"""
    Dataset(filename::AbstractString,mode::AbstractString = "r";
                     format::Symbol = :netcdf4, attrib = [])

Create a new NetCDF file if the `mode` is `"c"`. An existing file with the same
name will be overwritten. If `mode` is `"a"`, then an existing file is open into
append mode (i.e. existing data in the netCDF file is not overwritten and
a variable can be added). With the mode set to `"r"`, an existing netCDF file or
OPeNDAP URL can be open in read-only mode.  The default mode is `"r"`.
The optional parameter `attrib` is an iterable of attribute name and attribute
value pairs, for example a `Dict`, `DataStructures.OrderedDict` or simply a
vector of pairs (see example below).

# Supported formats:

* `:netcdf4` (default): HDF5-based NetCDF format.
* `:netcdf4_classic`: Only netCDF 3 compatible API features will be used.
* `:netcdf3_classic`: classic netCDF format supporting only files smaller than 2GB.
* `:netcdf3_64bit_offset`: improved netCDF format supporting files larger than 2GB.

Files can also be open and automatically closed with a `do` block.

```julia
Dataset("file.nc") do ds
    data = ds["temperature"][:,:]
end
```

```julia
Dataset("file.nc", "c", attrib = ["title" => "my first netCDF file"]) do ds
   defVar(ds,"temp",[10.,20.,30.],("time",))
end;
```

"""
function Dataset(filename::AbstractString,
                 mode::AbstractString = "r";
                 format::Symbol = :netcdf4,
                 diskless = false,
                 persist = false,
                 attrib = [])

    ncid = -1
    isdefmode = [false]

    if (mode == "r") || (mode == "a")
        ncmode =
            if (mode == "r")
                NC_NOWRITE
            else
                NC_WRITE
            end

        if diskless
            ncmode = ncmode | NC_DISKLESS
        end

        ncid = nc_open(filename,ncmode)
    elseif mode == "c"
        ncmode  = NC_CLOBBER

        if format == :netcdf3_64bit_offset
            ncmode = ncmode | NC_64BIT_OFFSET
        elseif format == :netcdf4_classic
            ncmode = ncmode | NC_NETCDF4 | NC_CLASSIC_MODEL
        elseif format == :netcdf4
            ncmode = ncmode | NC_NETCDF4
        elseif format == :netcdf3_classic
            # do nothing
        else
            throw(NetCDFError(-1, "Unkown format '$(format)' for filename '$(filename)'"))
        end

        if diskless
            ncmode = ncmode | NC_DISKLESS

            if persist
                ncmode = ncmode | NC_PERSIST
            end
        end

        ncid = nc_create(filename,ncmode)
        isdefmode[1] = true
    else
        throw(NetCDFError(-1, "Unsupported mode '$(mode)' for filename '$(filename)'"))
    end

    ds = Dataset(ncid,isdefmode)

    # set global attributes
    for (attname,attval) in attrib
        ds.attrib[attname] = attval
    end

    return ds
end

function Dataset(ncid::Integer,
                 isdefmode::Vector{Bool})
    attrib = Attributes(ncid,NC_GLOBAL,isdefmode)
    dim = Dimensions(ncid,isdefmode)
    group = Groups(ncid,isdefmode)
    return Dataset(ncid,isdefmode,attrib,dim,group)
end

function Dataset(f::Function,args...; kwargs...)
    ds = Dataset(args...; kwargs...)
    try
        f(ds)
    finally
        #@debug "closing netCDF Dataset $(NCDatasets.path(ds))"
        close(ds)
    end
end

"""
    defDim(ds::Dataset,name,len)

Define a dimension in the data set `ds` with the given `name` and length `len`.
If `len` is the special value `Inf`, then the dimension is considered as
`unlimited`, i.e. it will grow as data is added to the NetCDF file.

For example:

```julia
ds = Dataset("/tmp/test.nc","c")
defDim(ds,"lon",100)
```

This defines the dimension `lon` with the size 100.
"""
function defDim(ds::Dataset,name,len)
    defmode(ds.ncid,ds.isdefmode) # make sure that the file is in define mode
    dimid = nc_def_dim(ds.ncid,name,(isinf(len) ? NC_UNLIMITED : len))
    return nothing
end

function renameDim(ds::Dataset,oldname,newname)
    defmode(ds.ncid,ds.isdefmode) # make sure that the file is in define mode
    dimid = nc_inq_dimid(ds.ncid,oldname)
    nc_rename_dim(ds.ncid,dimid,newname)
    return nothing
end
export renameDim

"""
    defVar(ds::Dataset,name,vtype,dimnames; kwargs...)
    defVar(ds::Dataset,name,data,dimnames; kwargs...)

Define a variable with the name `name` in the dataset `ds`.  `vtype` can be
Julia types in the table below (with the corresponding NetCDF type). The
parameter `dimnames` is a tuple with the names of the dimension.  For scalar
this parameter is the empty tuple `()`.
The variable is returned (of the type CFVariable).

Instead of
providing the variable type one can directly give also the data `data` which
will be used to fill the NetCDF variable. In this case, the dimensions with
the appropriate size will be created as required using the names in `dimnames`.

!!! note

    Note if `data` is a vector or array of `DateTime` objects, then the dates
    are saved as double-precision floats and units
    "$(CFTime.DEFAULT_TIME_UNITS)" (unless a time unit
    is specifed with the `attrib` keyword as described below). Dates are
    converted to the default calendar in the CF conversion which is the
    mixed Julian/Gregorian calendar.

## Keyword arguments

* `fillvalue`: A value filled in the NetCDF file to indicate missing data.
   It will be stored in the _FillValue attribute.
* `chunksizes`: Vector integers setting the chunk size. The total size of a chunk must be less than 4 GiB.
* `deflatelevel`: Compression level: 0 (default) means no compression and 9 means maximum compression. Each chunk will be compressed individually.
* `shuffle`: If true, the shuffle filter is activated which can improve the compression ratio.
* `checksum`: The checksum method can be `:fletcher32` or `:nochecksum` (checksumming is disabled, which is the default)
* `attrib`: An iterable of attribute name and attribute value pairs, for example a `Dict`, `DataStructures.OrderedDict` or simply a vector of pairs (see example below)
* `typename` (string): The name of the NetCDF type required for [vlen arrays](https://web.archive.org/save/https://www.unidata.ucar.edu/software/netcdf/netcdf-4/newdocs/netcdf-c/nc_005fdef_005fvlen.html)

`chunksizes`, `deflatelevel`, `shuffle` and `checksum` can only be
set on NetCDF 4 files.

## NetCDF data types

| NetCDF Type | Julia Type |
|-------------|------------|
| NC_BYTE     | Int8 |
| NC_UBYTE    | UInt8 |
| NC_SHORT    | Int16 |
| NC_INT      | Int32 |
| NC_INT64    | Int64 |
| NC_FLOAT    | Float32 |
| NC_DOUBLE   | Float64 |
| NC_CHAR     | Char |
| NC_STRING   | String |

## Example:

```julia-repl
julia> data = randn(3,5)
julia> Dataset("test_file.nc","c") do ds
          defVar(ds,"temp",data,("lon","lat"), attrib = [
             "units" => "degree_Celsius",
             "long_name" => "Temperature"
          ])
       end;

```
"""
function defVar(ds::Dataset,name,vtype::DataType,dimnames; kwargs...)
    # all keyword arguments as dictionary
    kw = Dict(k => v for (k,v) in kwargs)

    defmode(ds.ncid,ds.isdefmode) # make sure that the file is in define mode
    dimids = Cint[nc_inq_dimid(ds.ncid,dimname) for dimname in dimnames[end:-1:1]]

    typeid =
        if vtype <: Vector
            # variable-length type
            typeid = nc_def_vlen(ds.ncid, kw[:typename], ncType[eltype(vtype)])
        else
            # base-type
            ncType[vtype]
        end

    varid = nc_def_var(ds.ncid,name,typeid,dimids)

    if haskey(kw,:chunksizes)
        storage = :chunked
        chunksizes = kw[:chunksizes]

        # this will fail on NetCDF-3 files
        nc_def_var_chunking(ds.ncid,varid,storage,reverse(chunksizes))
    end

    if haskey(kw,:shuffle) || haskey(kw,:deflatelevel)
        shuffle = get(kw,:shuffle,false)
        deflate = haskey(kw,:deflatelevel)
        deflate_level = get(kw,:deflatelevel,0)

        # this will fail on NetCDF-3 files
        nc_def_var_deflate(ds.ncid,varid,shuffle,deflate,deflate_level)
    end

    if haskey(kw,:checksum)
        checksum = kw[:checksum]
        nc_def_var_fletcher32(ds.ncid,varid,checksum)
    end

    if haskey(kw,:fillvalue)
        fillvalue = kw[:fillvalue]
        nofill = get(kw,:nofill,false)
        nc_def_var_fill(ds.ncid, varid, nofill, vtype(fillvalue))
    end

    if haskey(kw,:attrib)
        v = ds[name]
        for (attname,attval) in kw[:attrib]
            v.attrib[attname] = attval
        end
    end

    return ds[name]
end

# data has the type e.g. Array{Union{Missing,Float64},3}
function defVar(ds::Dataset,
                name,
                data::AbstractArray{Union{Missing,nctype},N},
                dimnames;
                kwargs...) where nctype <: Union{Int8,UInt8,Int16,Int32,Int64,Float32,Float64} where N
    _defVar(ds::Dataset,name,data,nctype,dimnames; kwargs...)
end

# data has the type e.g. Vector{DateTime}, Array{Union{Missing,DateTime},3} or
# Vector{DateTime360Day}
# Data is always stored as Float64 in the NetCDF file
function defVar(ds::Dataset,
                name,
                data::AbstractArray{<:Union{Missing,nctype},N},
                dimnames;
                kwargs...) where nctype <: Union{DateTime,AbstractCFDateTime} where N
    _defVar(ds::Dataset,name,data,Float64,dimnames; kwargs...)
end

function defVar(ds::Dataset,name,data,dimnames; kwargs...)
    nctype = eltype(data)
    _defVar(ds::Dataset,name,data,nctype,dimnames; kwargs...)
end

function _defVar(ds::Dataset,name,data,nctype,dimnames; kwargs...)
    # define the dimensions if necessary
    for i = 1:length(dimnames)
        if !(dimnames[i] in ds.dim)
            ds.dim[dimnames[i]] = size(data,i)
        end
    end

    v =
        if Missing <: eltype(data)
            # make sure a fill value is set (it might be overwritten by kwargs...)
            defVar(ds,name,nctype,dimnames;
                   fillvalue = ncFillValue[nctype],
                   kwargs...)
        else
            defVar(ds,name,nctype,dimnames; kwargs...)
        end
    v[:] = data
    return v
end

function defVar(ds::Dataset,name,data::T; kwargs...) where T <: Number
    v = defVar(ds,name,T,(); kwargs...)
    v[:] = data
    return v
end

function renameVar(ds::Dataset,oldname,newname)
    defmode(ds.ncid,ds.isdefmode) # make sure that the file is in define mode
    varid = nc_inq_varid(ds.ncid,oldname)
    nc_rename_var(ds.ncid,varid,newname)
    return nothing
end
export renameVar

"""
    keys(ds::Dataset)

Return a list of all variables names in Dataset `ds`.
"""
Base.keys(ds::Dataset) = listVar(ds.ncid)

"""
    path(ds::Dataset)
Return the file path (or the opendap URL) of the Dataset `ds`
"""
path(ds::Dataset) = nc_inq_path(ds.ncid)

"""
    groupname(ds::Dataset)
Return the group name of the Dataset `ds`
"""
groupname(ds::Dataset) = nc_inq_grpname(ds.ncid)

"""
    sync(ds::Dataset)

Write all changes in Dataset `ds` to the disk.
"""
function sync(ds::Dataset)
    datamode(ds.ncid,ds.isdefmode)
    nc_sync(ds.ncid)
end

"""
    close(ds::Dataset)

Close the Dataset `ds`. All pending changes will be written
to the disk.
"""
Base.close(ds::Dataset) = nc_close(ds.ncid)

"""
    v = variable(ds::Dataset,varname::String)

Return the NetCDF variable `varname` in the dataset `ds` as a
`NCDataset.Variable`. No scaling or other transformations are applied when the
variable `v` is indexed.
"""
function variable(ds::Dataset,varname::AbstractString)
    varid = nc_inq_varid(ds.ncid,varname)
    name,nctype,dimids,nattr = nc_inq_var(ds.ncid,varid)
    ndims = length(dimids)
    #@show ndims
    shape = zeros(Int,ndims)
    #@show typeof(shape),typeof(Int(1))

    for i = 1:ndims
        shape[ndims-i+1] = nc_inq_dimlen(ds.ncid,dimids[i])
    end
    #@show shape
    #@show typeof(shape)

    attrib = Attributes(ds.ncid,varid,ds.isdefmode)

    # reverse dimids to have the dimension order in Fortran style
    return Variable{nctype,ndims}(ds.ncid,varid,
                                  #(shape...),
                                  (reverse(dimids)...,),
                                  attrib,ds.isdefmode)
end

function Base.show(io::IO,ds::AbstractDataset; indent="")
    try
        dspath = path(ds)
        printstyled(io, indent, "Dataset: ",dspath,"\n", color=:red)
    catch err
        if isa(err,NetCDFError)
            if err.code == NC_EBADID
                print(io,"closed NetCDF Dataset")
                return
            end
        end
        rethrow
    end

    print(io,indent,"Group: ",groupname(ds),"\n")
    print(io,"\n")

    dims = collect(ds.dim)

    if length(dims) > 0
        printstyled(io, indent, "Dimensions\n",color=:red)

        for (dimname,dimlen) in dims
            print(io,indent,"   $(dimname) = $(dimlen)\n")
        end
        print(io,"\n")
    end

    varnames = keys(ds)

    if length(varnames) > 0

        printstyled(io, indent, "Variables\n",color=:red)

        for name in varnames
            show(io,variable(ds,name); indent = "$(indent)  ")
            print(io,"\n")
        end
    end

    # global attribues
    if length(ds.attrib) > 0
        printstyled(io, indent, "Global attributes\n",color=:red)
        show(io,ds.attrib; indent = "$(indent)  ");
    end

    # groups
    groupnames = keys(ds.group)

    if length(groupnames) > 0
        printstyled(io, indent, "Groups\n",color = :red)
        for groupname in groupnames
            show(io,group(ds,groupname); indent = "  ")
        end
    end

end

"""
    v = getindex(ds::Dataset,varname::AbstractString)

Return the NetCDF variable `varname` in the dataset `ds` as a
`NCDataset.CFVariable`. The CF convention are honored when the
variable is indexed:
* `_FillValue` will be returned as `missing`
* `scale_factor` and `add_offset` are applied
* time variables (recognized by the units attribute) are returned usually as
`DateTime` object. Note that `DateTimeAllLeap`, `DateTimeNoLeap` and
`DateTime360Day` cannot be converted to the proleptic gregorian calendar used in
julia and are returned as such.

A call `getindex(ds,varname)` is usually written as `ds[varname]`.
"""
function Base.getindex(ds::AbstractDataset,varname::AbstractString)
    v = variable(ds,varname)

    # return element type of any index operation
    if eltype(v) <: Number
        rettype = Union{Missing,Number,DateTime,AbstractCFDateTime}
    else
        rettype = Union{Missing,eltype(v)}
    end

    return CFVariable{rettype,ndims(v),typeof(v),typeof(v.attrib)}(v,v.attrib)
end

# -----------------------------------------------------
# Variable (as stored in NetCDF file, without using
# add_offset, scale_factor and _FillValue)

mutable struct Variable{NetCDFType,N} <: AbstractVariable{NetCDFType, N}
    ncid::Cint
    varid::Cint
    dimids::NTuple{N,Cint}
    attrib::Attributes
    isdefmode::Vector{Bool}
end

mutable struct MFVariable{T,N,M,TA} <: AbstractVariable{T,N}
    var::CatArrays.CatArray{T,N,M,TA}
    attrib::MFAttributes
    dimnames::NTuple{N,String}
    varname::String
end

mutable struct DeferVariable{T,N} <: AbstractVariable{T,N}
    r::Resource
    varname::String
    attrib::DeferAttributes
    data::OrderedDict
end

Dataset(var::Variable) = Dataset(var.ncid,var.isdefmode)

# the size of a variable can change, i.e. for a variable with an unlimited
# dimension
Base.size(v::Variable) = (Int[nc_inq_dimlen(v.ncid,dimid) for dimid in v.dimids]...,)

"""
    dimnames(v::Variable)

Return a tuple of the dimension names of the variable `v`.
"""
function dimnames(v::Variable)
    return (String[nc_inq_dimname(v.ncid,dimid) for dimid in v.dimids]...,)
end

"""
    name(v::Variable)

Return the name of the NetCDF variable `v`.
"""
name(v::Variable) = nc_inq_varname(v.ncid,v.varid)

chunking(v::Variable,storage,chunksizes) = nc_def_var_chunking(v.ncid,v.varid,storage,reverse(chunksizes))

"""
    storage,chunksizes = chunking(v::Variable)

Return the storage type (`:contiguous` or `:chunked`) and the chunk sizes
of the varable `v`.
"""
function chunking(v::Variable)
    storage,chunksizes = nc_inq_var_chunking(v.ncid,v.varid)
    return storage,reverse(chunksizes)
end

"""
    isshuffled,isdeflated,deflate_level = deflate(v::Variable)

Return compression information of the variable `v`. If shuffle
is `true`, then shuffling (byte interlacing) is activaded. If
deflate is `true`, then the data chunks (see `chunking`) are
compressed using the compression level `deflate_level`
(0 means no compression and 9 means maximum compression).
"""
deflate(v::Variable,shuffle,deflate,deflate_level) = nc_def_var_deflate(v.ncid,v.varid,shuffle,deflate,deflate_level)
deflate(v::Variable) = nc_inq_var_deflate(v.ncid,v.varid)

checksum(v::Variable,checksummethod) = nc_def_var_fletcher32(v.ncid,v.varid,checksummethod)

"""
   checksummethod = checksum(v::Variable)

Return the checksum method of the variable `v` which can be either
be `:fletcher32` or `:nochecksum`.
"""
checksum(v::Variable) = nc_inq_var_fletcher32(v.ncid,v.varid)

function fillmode(v::Variable)
    no_fill,fv = nc_inq_var_fill(v.ncid, v.varid)
    return no_fill,fv
end

function fillvalue(v::Variable)
    no_fill,fv = nc_inq_var_fill(v.ncid, v.varid)
    return fv
end

function Base.getindex(v::Variable,indexes::Int...)
    #    @show "ind",indexes
    return nc_get_var1(eltype(v),v.ncid,v.varid,[i-1 for i in indexes[ndims(v):-1:1]])
end

function Base.setindex!(v::Variable{T,N},data,indexes::Int...) where N where T
    @debug "$(@__LINE__)"
    datamode(v.ncid,v.isdefmode)
    # use zero-based indexes and reversed order
    nc_put_var1(v.ncid,v.varid,[i-1 for i in indexes[ndims(v):-1:1]],T(data))
    return data
end

function Base.getindex(v::Variable{T,N},indexes::Colon...) where {T,N}
    # special case for scalar NetCDF variable
    if N == 0
        data = Vector{T}(undef,1)
        nc_get_var!(v.ncid,v.varid,data)
        return data[1]
    else
        #@show v.shape,typeof(v.shape),T,N
        #@show v.ncid,v.varid
        data = Array{T,N}(undef,size(v))
        nc_get_var!(v.ncid,v.varid,data)
        return data
    end
end

function Base.setindex!(v::Variable{T,N},data::T,indexes::Colon...) where {T,N}
    @debug "setindex! colon $data"
    datamode(v.ncid,v.isdefmode) # make sure that the file is in data mode
    tmp = fill(data,size(v))
    #@show "here number",indexes,size(v),fill(data,size(v))
    nc_put_var(v.ncid,v.varid,tmp)
    return data
end

# call to v .= 123
function Base.setindex!(v::Variable{T,N},data::Number) where {T,N}
    @debug "setindex! $data"
    datamode(v.ncid,v.isdefmode) # make sure that the file is in data mode
    tmp = fill(convert(T,data),size(v))
    #@show "here number",indexes,size(v),tmp
    nc_put_var(v.ncid,v.varid,tmp)
    return data
end

Base.setindex!(v::Variable,data::Number,indexes::Colon...) = setindex!(v::Variable,data)

function Base.setindex!(v::Variable{T,N},data::AbstractArray{T,N},indexes::Colon...) where {T,N}
    datamode(v.ncid,v.isdefmode) # make sure that the file is in data mode
    #@show @__LINE__,@__FILE__
    nc_put_var(v.ncid,v.varid,data)
    return data
end

function Base.setindex!(v::Variable{T,N},data::AbstractArray{T2,N},indexes::Colon...) where {T,T2,N}
    datamode(v.ncid,v.isdefmode) # make sure that the file is in data mode
    tmp =
        if T <: Integer
            round.(T,data)
        else
            convert(Array{T,N},data)
        end

    nc_put_var(v.ncid,v.varid,tmp)
    return data
end

function ncsub(indexes)
    count = [length(i) for i in indexes[end:-1:1]]
    start = [first(i)-1 for i in indexes[end:-1:1]]     # use zero-based indexes
    stride = [step(i) for i in indexes[end:-1:1]]
    jlshape = (count[end:-1:1]...,)
    return start,count,stride,jlshape
end

function Base.getindex(v::Variable{T,N},indexes::TR...) where {T,N} where TR <: Union{StepRange{Int,Int},UnitRange{Int}}
    start,count,stride,jlshape = ncsub(indexes[1:ndims(v)])
    data = Array{T,N}(undef,jlshape)
    nc_get_vars!(v.ncid,v.varid,start,count,stride,data)
    return data
end

function Base.setindex!(v::Variable{T,N},data::T,indexes::StepRange{Int,Int}...) where {T,N}
    #@show @__FILE__,@__LINE__,indexes
    datamode(v.ncid,v.isdefmode) # make sure that the file is in data mode
    start,count,stride,jlshape = ncsub(indexes[1:ndims(v)])
    tmp = fill(data,jlshape)
    nc_put_vars(v.ncid,v.varid,start,count,stride,tmp)
    return data
end

function Base.setindex!(v::Variable{T,N},data::Number,indexes::StepRange{Int,Int}...) where {T,N}
    #@show @__FILE__,@__LINE__,indexes
    datamode(v.ncid,v.isdefmode) # make sure that the file is in data mode
    start,count,stride,jlshape = ncsub(indexes[1:ndims(v)])
    tmp = fill(convert(T,data),jlshape)
    nc_put_vars(v.ncid,v.varid,start,count,stride,tmp)
    return data
end

function Base.setindex!(v::Variable{T,N},data::Array{T,N},indexes::StepRange{Int,Int}...) where {T,N}
    #@show "sr",indexes
    datamode(v.ncid,v.isdefmode) # make sure that the file is in data mode
    start,count,stride,jlshape = ncsub(indexes[1:ndims(v)])
    nc_put_vars(v.ncid,v.varid,start,count,stride,data)
    return data
end

# data can be Array{T2,N} or BitArray{N}
function Base.setindex!(v::Variable{T,N},data::AbstractArray,indexes::StepRange{Int,Int}...) where {T,N}
    #@show "sr2",indexes
    datamode(v.ncid,v.isdefmode) # make sure that the file is in data mode
    start,count,stride,jlshape = ncsub(indexes[1:ndims(v)])

    tmp = convert(Array{T,ndims(data)},data)
    nc_put_vars(v.ncid,v.varid,start,count,stride,tmp)

    return data
end

function normalizeindexes(sz,indexes)
    ndims = length(sz)
    ind = Vector{StepRange}(undef,ndims)
    squeezedim = falses(ndims)

    # normalize indexes
    for i = 1:ndims
        indT = typeof(indexes[i])
        # :
        if indT == Colon
            ind[i] = 1:1:sz[i]
            # just a number
        elseif indT == Int
            ind[i] = indexes[i]:1:indexes[i]
            squeezedim[i] = true
            # range with a step equal to 1
        elseif indT == UnitRange{Int}
            ind[i] = first(indexes[i]):1:last(indexes[i])
        elseif indT == StepRange{Int,Int}
            ind[i] = indexes[i]
        else
            #@show indT
            error("unsupported index")
        end
    end

    return ind,squeezedim
end

function Base.getindex(v::Variable,indexes::Union{Int,Colon,UnitRange{Int},StepRange{Int,Int}}...)
    #    @show "any",indexes
    ind,squeezedim = normalizeindexes(size(v),indexes)

    data = v[ind...]
    # drop any dimension which was indexed with a scalar
    if any(squeezedim)
        return dropdims(data,dims=(findall(squeezedim)...,))
    else
        return data
    end
end

function Base.setindex!(v::Variable,data,indexes::Union{Int,Colon,UnitRange{Int},StepRange{Int,Int}}...)
    #@show "any",indexes
    ind,squeezedim = normalizeindexes(size(v),indexes)

    # make arrays out of scalars
    if ndims(data) == 0
        data = fill(data,([length(i) for i in ind]...,))
    end

    if ndims(data) == 1 && size(data,1) == 1
        data = fill(data[1],([length(i) for i in ind]...,))
    end

    # return data
    return v[ind...] = data
end

"""
    load!(ncvar::Variable, data, indices)

Loads a NetCDF variables `ncvar` and puts the result in `data` along the
specified `indices`.

```julia
data = zeros(5,6); # must have the right shape and type
load!(ds["temp"].var,data,:,:) # loads all data

data = zeros(5); # must have the right shape and type
load!(ds["temp"].var,data,:,1) # loads the 1st column
```
"""
@inline function load!(ncvar::NCDatasets.Variable{T,N}, data, indices::Union{Integer, UnitRange, StepRange, Colon}...) where {T,N}
    ind = to_indices(ncvar,indices)
    start,count,stride,jlshape = ncsub(ind)
    nc_get_vars!(ncvar.ncid,ncvar.varid,start,count,stride,data)
end

@inline function load!(ncvar::NCDatasets.Variable{T,2}, data, i::Colon,j::UnitRange) where T
    # reversed and 0-based
    start = [first(j)-1,0]
    count = [length(j),size(ncvar,1)]
    nc_get_vara!(ncvar.ncid,ncvar.varid,start,count,data)
end

"""
     data = loadragged(ncvar,index::Colon)

Load data from `ncvar` in the [contiguous ragged array representation](https://web.archive.org/web/20190111092546/http://cfconventions.org/cf-conventions/v1.6.0/cf-conventions.html#_contiguous_ragged_array_representation) as a
vector of vectors. It is typically used to load a list of profiles
or time series of different length each.

The [indexed ragged array representation](https://web.archive.org/web/20190111092546/http://cfconventions.org/cf-conventions/v1.6.0/cf-conventions.html#_indexed_ragged_array_representation) is currently not supported.
"""
function loadragged(ncvar,index::Colon)
    ds = Dataset(ncvar.var.ncid,ncvar.var.isdefmode)

    dimensionnames = dimnames(ncvar)
    if length(dimensionnames) !== 1
        throw(NetCDFError(-1, "NetCDF variable $(name(ncvar)) should have only one dimensions"))
    end
    dimname = dimensionnames[1]

    ncvarsizes = varbyattrib(ds,sample_dimension = dimname)
    if length(ncvarsizes) !== 1
        throw(NetCDFError(-1, "There should be exactly one NetCDF variable with the attribiute sample_dimension equal to '$(dimname)'"))
    end

    ncvarsize = ncvarsizes[1]
    varsize = ncvarsize.var[:]

    istart = 0;
    tmp = ncvar[:]

    T = typeof(view(tmp,1:varsize[1]))
    data = Vector{T}(undef,length(varsize))

    for i = 1:length(varsize)
        data[i] = view(tmp,istart+1:istart+varsize[i]);
        istart += varsize[i]
    end
    return data
end

# -----------------------------------------------------
# Variable (with applied transformations following the CF convention)

mutable struct CFVariable{T,N,TV,TA}  <: AbstractArray{T, N}
    var::TV
    attrib::TA
end

Dataset(var::CFVariable) = Dataset(var.var)

Base.size(v::CFVariable) = size(v.var)
dimnames(v::CFVariable)  = dimnames(v.var)
name(v::CFVariable)  = name(v.var)
chunking(v::CFVariable,storage,chunksize) = chunking(v.var,storage,chunksize)
chunking(v::CFVariable) = chunking(v.var)

deflate(v::CFVariable,shuffle,dodeflate,deflate_level) = deflate(v.var,shuffle,dodeflate,deflate_level)
deflate(v::CFVariable) = deflate(v.var)

checksum(v::CFVariable,checksummethod) = checksum(v.var,checksummethod)
checksum(v::CFVariable) = checksum(v.var)

fillmode(v::CFVariable) = fillmode(v.var)
fillvalue(v::CFVariable) = fillvalue(v.var)

function Base.getindex(v::CFVariable,indexes::Union{Int,Colon,UnitRange{Int},StepRange{Int,Int}}...)
    attnames = keys(v.attrib)

    data = v.var[indexes...]
    isscalar =
        if typeof(data) == String || typeof(data) == DateTime
            true
        else
            ndims(data) == 0
        end

    if isscalar
        data = [data]
    end

    if "_FillValue" in attnames
        fillvalue = v.attrib["_FillValue"]
        mask = isfillvalue(data,fillvalue)
    else
        mask = falses(size(data))
    end

    # do not scale characters and strings
    if eltype(v.var) != Char
        if "scale_factor" in attnames
            data = v.attrib["scale_factor"] * data
        end

        if "add_offset" in attnames
            data = data .+ v.attrib["add_offset"]
        end
    end

    if "units" in attnames
        units = v.attrib["units"]
        if occursin(" since ",units)
            # type of data changes
            calendar = lowercase(get(v.attrib,"calendar","standard"))
            # time decode only valid dates
            tmp = timedecode(data[.!mask],units,calendar)
            data = similar(tmp,size(data))
            data[.!mask] = tmp
        end
    end

    if isscalar
        if mask[1]
            missing
        else
            data[1]
        end
    else
        datam = Array{Union{eltype(data),Missing}}(data)
        datam[mask] .= missing
        return datam
    end
end

function Base.setindex!(v::CFVariable,data::Array{Missing,N},indexes::Union{Int,Colon,UnitRange{Int},StepRange{Int,Int}}...) where N
    v.var[indexes...] = fill(v.attrib["_FillValue"],size(data))
end

function Base.setindex!(v::CFVariable,data::Missing,indexes::Union{Int,Colon,UnitRange{Int},StepRange{Int,Int}}...)
    v.var[indexes...] = v.attrib["_FillValue"]
end

isfillvalue(data,fillvalue) = data .== fillvalue
isfillvalue(data,fillvalue::Number) = (isnan(fillvalue) ? isnan.(data) : data .== fillvalue)

function Base.setindex!(v::CFVariable,data::Union{T,Array{T,N}},indexes::Union{Int,Colon,UnitRange{Int},StepRange{Int,Int}}...) where N where T <: Union{AbstractCFDateTime,DateTime,Union{Missing,DateTime,AbstractCFDateTime}}
    attnames = keys(v.attrib)
    units =
        if "units" in attnames
            v.attrib["units"]
        else
            @debug "set time units to $CFTime.DEFAULT_TIME_UNITS"
            v.attrib["units"] = CFTime.DEFAULT_TIME_UNITS
            CFTime.DEFAULT_TIME_UNITS
        end

    if occursin(" since ",units)
        if !("calendar" in attnames)
            # these dates cannot be converted to the standard calendar
            if T <: Union{DateTime360Day,Missing}
                v.attrib["calendar"] = "360_day"
            elseif T <: Union{DateTimeNoLeap,Missing}
                v.attrib["calendar"] = "365_day"
            elseif T <: Union{DateTimeAllLeap,Missing}
                v.attrib["calendar"] = "366_day"
            end
        end
        calendar = lowercase(get(v.attrib,"calendar","standard"))
        # can throw an convertion error if calendar attribute already exists and
        # is incompatible with the provided data
        v[indexes...] = timeencode(data,units,calendar)
        return data
    end

    throw(NetCDFError(-1, "time unit ('$units') of the variable $(name(v)) does not include the word ' since '"))
end

function transform(v,offset,scale)
    if offset !== nothing
        if scale !== nothing
            return (v - offset) / scale
        else
            return v - offset
        end
    else
        if scale !== nothing
            return v / scale
        else
            return v
        end
    end
end

# the transformv function is necessary to avoid "iterating" over a single character
# https://discourse.julialang.org/t/broadcasting-and-single-characters/16836
transformv(data,offset,scale) = transform.(data,offset,scale)
transformv(data::Char,offset,scale) = data

function Base.setindex!(v::CFVariable,data,indexes::Union{Int,Colon,UnitRange{Int},StepRange{Int,Int}}...)
    offset,scale =
        if eltype(v.var) == Char
            (nothing,nothing)
        else
            (get(v.attrib,"add_offset",nothing),
             get(v.attrib,"scale_factor",nothing))
        end

    fillvalue = get(v.attrib,"_FillValue",nothing)

    v.var[indexes...] =
        if fillvalue == nothing
            transformv(data,offset,scale)
        else
            coalesce.(transformv(data,offset,scale),fillvalue)
        end

    return data
end

function Base.show(io::IO,v::AbstractVariable; indent="")
    delim = " × "
    dims =
        try
            dimnames(v)
        catch err
            if isa(err,NetCDFError)
                if err.code == NC_EBADID
                    print(io,"NetCDF variable (file closed)")
                    return
                end
            end
            rethrow
        end
    sz = size(v)

    printstyled(io, indent, name(v),color=:green)
    if length(sz) > 0
        print(io,indent," (",join(sz,delim),")\n")
        print(io,indent,"  Datatype:    ",eltype(v),"\n")
        print(io,indent,"  Dimensions:  ",join(dims,delim),"\n")
    else
        print(io,indent,"\n")
    end

    if length(v.attrib) > 0
        print(io,indent,"  Attributes:\n")
        show(io,v.attrib; indent = "$(indent)   ")
    end
end

Base.show(io::IO,v::CFVariable; indent="") = Base.show(io::IO,v.var; indent=indent)

Base.display(v::Union{Variable,CFVariable}) = show(Compat.stdout,v)

Base.Array(v::Union{CFVariable,Variable}) = v[:]

# Common methods

const NCIterable = Union{BaseAttributes,AbstractDimensions,AbstractDataset,AbstractGroups}

Base.length(a::NCIterable) = length(keys(a))

"""
    haskey(ds::Dataset,varname)

Return true if the Dataset `ds` has a variable with the name `varname`.
For example:

```julia
ds = Dataset("/tmp/test.nc","r")
if haskey(ds,"temperature")
    println("The file has a variable 'temperature'")
end
```

This example checks if the file `/tmp/test.nc` has a variable with the
name `temperature`.
"""
Base.haskey(a::NCIterable,name::AbstractString) = name in keys(a)
Base.in(name::AbstractString,a::NCIterable) = name in keys(a)

function Base.iterate(a::NCIterable, state = keys(a))
    if length(state) == 0
        return nothing
    end

    return (state[1] => a[popfirst!(state)], state)
end

"""
    escape(val)

Escape backslash, dollar and quote from string `val`.
"""
function escape(val)
     valescaped = val
     # backslash must come first
     for c in ['\\','$','"']
        valescaped = replace(valescaped,c => "\\$c")
    end
	return valescaped
end

function ncgen(io::IO,fname; newfname = "filename.nc")
    ds = Dataset(fname)
    unlimited_dims = unlimited(ds.dim)

    print(io,"ds = Dataset(\"$(escape(newfname))\",\"c\")\n")

    print(io,"# Dimensions\n\n")
    for (d,v) in ds.dim
        if d in unlimited_dims
            print(io,"ds.dim[\"$d\"] = Inf # unlimited dimension\n")
        else
            print(io,"ds.dim[\"$d\"] = $v\n")
        end
    end

    print(io,"\n# Declare variables\n\n")

    for (d,v) in ds
        print(io,"nc$d = defVar(ds,\"$d\", $(eltype(v.var)), $(dimnames(v)))\n")
        ncgen_setattrib(io,"nc$d",v.attrib)
        print(io,"\n")
    end

    print(io,"# Global attributes\n\n")

    ncgen_setattrib(io,"ds",ds.attrib)

    print(io,"\n# Define variables\n\n")

    for d in keys(ds)
        print(io,"# nc$d[:] = ...\n")
    end

    print(io,"\nclose(ds)\n")
    close(ds)
end

"""
    ncgen(fname; ...)
    ncgen(fname,jlname; ...)

Generate the Julia code that would produce a NetCDF file with the same metadata
as the NetCDF file `fname`. The code is placed in the file `jlname` or printed
to the standard output. By default the new NetCDF file is called `filename.nc`.
This can be changed with the optional parameter `newfname`.
"""
ncgen(fname; kwargs...)  = ncgen(Compat.stdout,fname; kwargs...)

function ncgen(fname,jlname; kwargs...)
    open(jlname,"w") do io
        ncgen(io, fname; kwargs...)
    end
end

function ncgen_setattrib(io,v,attrib)
    for (d,val) in attrib
        litval = if typeof(val) == String
            "\"$(escape(val))\""
        elseif typeof(val) == Float64
            val
        elseif typeof(val) == Float32
            "$(eltype(val))($(val))"
        else
            val
        end

        print(io,"$(v).attrib[\"$d\"] = $litval\n");
    end
end

"""
    varbyattrib(ds, attname = attval)

Returns a list of variable(s) which has the attribute `attname` matching the value `attval`
in the dataset `ds`.
The list is empty if the none of the variables has the match.
The output is a list of `CFVariable`s.

# Examples

Load all the data of the first variable with standard name "longitude" from the
NetCDF file `results.nc`.

```julia-repl
julia> ds = Dataset("results.nc", "r");
julia> data = varbyattrib(ds, standard_name = "longitude")[1][:]
```

"""
function varbyattrib(ds::Dataset; kwargs...)
    # Start with an empty list of variables
    varlist = []

    # Loop on the variables
    for v in keys(ds)
        var = ds[v]

        matchall = true

        for (attsym,attval) in kwargs
            attname = String(attsym)

            # Check if the variable has the desired attribute
            if haskey(var.attrib, attname)
                # Check if the attribute value is the selected one
                if var.attrib[attname] != attval
                    matchall = false
                    break
                end
            else
                matchall = false
                break
            end
        end

        if matchall
            push!(varlist, var)
        end
    end

    return varlist
end

"""
    a = nomissing(da)

Retun the values of the array `da` of type `Array{Union{T,Missing},N}`
(potentially containing missing values) as a regular Julia array `a` of the same
element type and checks that no missing values are present.

"""
function nomissing(da::Array{Union{T,Missing},N}) where {T,N}
    if any(ismissing.(da))
        error("arrays contains missing values (values equal to the fill values attribute in the NetCDF file)")
    end
    if VERSION >= v"1.2"
        # Illegal instruction (core dumped) in Julia 1.0.5
        # but works on Julia 1.2
        return Array{T,N}(da)
    else
        # old
        if isempty(da)
            return Array{T,N}([])
        else
            return replace(da, missing => da[1])
        end
    end
end

"""
    a = nomissing(da,value)

Retun the values of the array `da` of type `Array{Union{T,Missing},N}`
as a regular Julia array `a` by replacing all missing value by `value`.
"""
function nomissing(da::Array{Union{T,Missing},N},value) where {T,N}
    return replace(da, missing => T(value))
end

export defVar, defDim, Dataset, close, sync, variable, dimnames, name,
    deflate, chunking, checksum, fillvalue, fillmode, ncgen
export nomissing
export varbyattrib
export path
export defGroup
export loadragged

include("multifile.jl")
export MFDataset, close

include("defer.jl")
export DeferDataset

include("cfconventions.jl")

# it is good practise to use the default fill-values, thus we export them
export NC_FILL_BYTE, NC_FILL_CHAR, NC_FILL_SHORT, NC_FILL_INT, NC_FILL_FLOAT,
    NC_FILL_DOUBLE, NC_FILL_UBYTE, NC_FILL_USHORT, NC_FILL_UINT, NC_FILL_INT64,
    NC_FILL_UINT64, NC_FILL_STRING

export CFTime
export daysinmonth, daysinyear, yearmonthday, yearmonth, monthday
export dayofyear, firstdayofyear

export DateTimeStandard, DateTimeJulian, DateTimeProlepticGregorian,
    DateTimeAllLeap, DateTimeNoLeap, DateTime360Day, AbstractCFDateTime

end # module

ext: jl


text:
name = "Chemfiles"
version = v"0.10.2"
chemfiles_sha256 = "2e3b58167f25d561ab19ae06acdc02f26b5640bd6c85e0a5b10fedfec59f5285"

max_stars_repo_name: chemfiles/chemfiles-prebuilt

max_issues_repo_name: chemfiles/chemfiles-prebuilt

max_forks_repo_name: chemfiles/chemfiles-prebuilt


text:
using GPMLj, Test, Random, PyCall

Random.seed!(123)

@testset "Models" begin
    @testset "GPR" begin
        kern= gpflow.kernels.Matern52(2)
        gpr = gpflow.models.GPR(randn(100,2), randn(100,1), kern;)
        @test typeof(gpr)<:gpflow.models.GPR
        @test typeof(gpr)<:gpflow.Model
        @test typeof(gpr.o)<:PyObject
        @test typeof(gpr.kern.o)<:PyObject

    end

    # FIXME
    #= @testset "SGPR" begin
        kern= gpflow.kernels.Matern52(2)
        sgpr = gpflow.models.SGPR(randn(100,2), randn(100,1), kern;)
        @test typeof(sgpr)<:gpflow.models.SGPR
        @test typeof(sgpr)<:gpflow.Model
        @test typeof(sgpr.o)<:PyObject
        @test typeof(sgpr.kern.o)<:PyObject

    end =#

    @testset "VGP" begin
        kern = gpflow.kernels.Matern52(2)
        like = gpflow.likelihoods.Gaussian()
        vgp = gpflow.models.VGP(randn(100,2), randn(100,1), kern, like)
        @test typeof(vgp)<:gpflow.models.VGP
        @test typeof(vgp)<:gpflow.Model
        @test typeof(vgp.o)<:PyObject
        @test typeof(vgp.kern.o)<:PyObject
        @test typeof(vgp.likelihood.o)<:PyObject
    end

    @testset "GPMC" begin
        kern= gpflow.kernels.Matern52(2)
        like = gpflow.likelihoods.Gaussian()
        gpmc = gpflow.models.GPMC(randn(100,2), randn(100,1), kern, like)
        @test typeof(gpmc)<:gpflow.models.GPMC
        @test typeof(gpmc)<:gpflow.Model
        @test typeof(gpmc.o)<:PyObject
        @test typeof(gpmc.kern.o)<:PyObject
        @test typeof(gpmc.likelihood.o)<:PyObject
    end

end #module


text:
module LambdaMaker

using Mustache
using UUIDs

export create_lambda_package, DirectoryExists

struct DirectoryExists <: Exception
    msg::AbstractString
end
Base.show(io::IO, e::DirectoryExists) = println(io, e.msg)

_src_path(filename::AbstractString) = joinpath(@__DIR__, "..", "template_files", filename)

function _copy_file(dest_dir::AbstractString, filename::AbstractString; force::Bool=false)
    src_path = _src_path(filename)
    dest_path = joinpath(dest_dir, filename)

    cp(src_path, dest_path; force=force)
end

function _render_file(filename::AbstractString, dest_path::AbstractString, render_args::AbstractDict)
    file_path = _src_path(filename)
    rendered_contents = render(read(file_path, String), render_args)

    dest_file_path = joinpath(dest_path, filename)

    open(dest_file_path, "w") do f
        print(f, rendered_contents)
    end
end

function create_lambda_package(package_name::AbstractString; force::Bool=false)
    render_args = Dict(
        "PKG" => package_name,
        "UUID" => string(UUIDs.uuid4())
    )

    # Create a project named directory at the current location
    if isdir(package_name) && !force
        throw(DirectoryExists("Directory $(package_name) already exists. Pass force=true to overwrite."))
    end

    dest_src_path = mkpath(joinpath(package_name, "src"))
    dest_path = dirname(dest_src_path)

    # Copy over bootstrap, Dockerfile, Manifest.toml, handle_requests.jl
    files_to_copy = ["bootstrap", "Manifest.toml", "handle_requests.jl", "template.yml"]
    _copy_file.(dest_path, files_to_copy; force=force)

    # Use Mustache to render the Dockerfile and Project.toml files
    _render_file("Project.toml", dest_path, render_args)
    _render_file("Dockerfile", dest_path, render_args)

    # Special case writing the rendered module.jl into src/package_name.jl
    module_path = joinpath(@__DIR__, "..", "template_files", "src", "module.jl")
    module_rendered = render(read(module_path, String), render_args)

    open(joinpath(dest_src_path, package_name * ".jl"), "w") do f
        print(f, module_rendered)
    end
end

end


text:
for dir in (:x, :y, :z)
    func! = Symbol(dir, :ticks!)
    attr1 = Symbol(dir, :axis)
    attr2 = Symbol(:set_, dir, :ticklabels)
    types = (AbstractRange, NTuple{N} where N)
    for X in types
        @eval begin
            function $(func!)(ax, locs::$X)
                ax.$(attr1).set_major_locator(plt.matplotlib.ticker.MaxNLocator(length(locs)))
                ax.$(attr1).set_major_locator(plt.matplotlib.ticker.FixedLocator(locs))
                return nothing
            end
            function $(func!)(ax, ticks::NTuple{N,S} where {N,S<:AbstractString}, locs::$X)
                $(func!)(ax, locs); ax.$(attr2)(ticks)
                return nothing
            end
        end
        for Y in types
            @eval begin
                function $(func!)(ax, ticks::$Y, locs::$X)
                    $(func!)(ax, locs); ax.$(attr2)(tick2label(ticks))
                    return nothing
                end
            end
        end
    end
end

function tick2label(ticks::AbstractRange{Int})
    labels = Vector{String}(undef, length(ticks))

    @inbounds @simd for i in eachindex(ticks)
        labels[i] = string(ticks[i])
    end

    return labels
end

function tick2label(ticks::AbstractRange{T}) where T<:Real
    labels = Vector{String}(undef, length(ticks))

    @inbounds @simd for i in eachindex(ticks)
        labels[i] = @sprintf("%.1f", ticks[i])
    end

    return labels
end

function tick2label(ticks::NTuple{N,Int}) where N
    if @generated
        e = Expr(:vect)
        a = Vector{Any}(undef, N)

        @inbounds for i in eachindex(a)
            a[i] = :(string(ticks[$i]))
        end

        e.args = a

        return quote
            $(Expr(:meta, :inline))
            @inbounds return $e
        end
    else
        labels = Vector{String}(undef, N)

        @inbounds @simd for i in eachindex(ticks)
            labels[i] = string(ticks[i])
        end

        return labels
    end
end

function tick2label(ticks::NTuple{N,T}) where {N,T<:Real}
    if @generated
        e = Expr(:vect)
        a = Vector{Any}(undef, N)

        @inbounds for i in eachindex(a)
            a[i] = :(@sprintf("%.1f", ticks[$i]))
        end

        e.args = a

        return quote
            $(Expr(:meta, :inline))
            @inbounds return $e
        end
    else
        labels = Vector{String}(undef, N)

        @inbounds @simd for i in eachindex(ticks)
            labels[i] = @sprintf("%.1f", ticks[i])
        end

        return labels
    end
end


text:
using Munkres
import Base: eye, rand, *, inv

# Type definitions

immutable SymmetricGroup{N} <: Group
    σ::Vector{Int}
end

Sn{N} = SymmetricGroup{N}
SymmetricGroup(σ::Vector{Int}) = SymmetricGroup{length(σ)}(σ)

# Identity element

eye{N}(::Type{Sn{N}}) = Sn{N}(1:N)

# Random group element

function rand{N}(::Type{Sn{N}})
    σ = shuffle(1:N)
    Sn(σ)
end

# Defining representation

@inbounds function defrep!{N, T}(M::Matrix{T}, g::Sn{N})
    σ = g.σ
    n = size(M, 1)
    for i in 1:n, j in 1:n
        if i == σ[j]
            M[i, j] = one(T)
        else
            M[i, j] = zero(T)
        end
    end

    return M
end

defrep{N}(g::Sn{N}) = defrep!(Matrix{Float64}(N, N), g)

function nearest{N, T}(::Type{Sn{N}}, M_defrep::Matrix{T})
    Sn{N}(munkres(-M_defrep'))
end

# Multiplication
function *{N}(g₁::Sn{N}, g₂::Sn{N})
    Sn{N}(g₁.σ[g₂.σ])
end

# Inversion
function inv{N}(g::Sn{N})
    Sn{N}(indexin(1:N, g.σ))
end


text:
using Compat
using Compat.Test

using Yao
using Yao.Zoo
using Yao.Blocks

using QCGANS, GradOptim

# TODO using an ideal discriminator!
witness_vec = zeros(1<<4)
witness_vec[[0, 3, 5, 10, 12, 15].+1] = 1
#loss(reg::DefaultRegister) = - witness_vec .* (reg |> probs)
using StatsBase
loss(reg::DefaultRegister) = kldivergence(witness_vec, (reg |> probs))

function vgradient(witness_vec::Vector, circuit::AbstractBlock, reg0::AbstractRegister, rots::Vector{RotationGate})
    loss_func = () -> -sum(witness_vec.*(copy(reg0) |> copy |> circuit |> probs))
    ptb = perturb(loss_func, rots, π/2)
    (ptb[:,1] - ptb[:,2])/2
end

using BenchmarkTools

n = 4
depth = 4
pairs = [1=>2, 2=>3, 3=>1]
generator = diff_circuit(n, depth, pairs)
rots = collect_rotblocks(generator)

T = randn(1<<n,1<<n)
T = T' + T
Trank = ndims(T)
witness_vec = T*(zero_state(4) |> generator |> probs)
gexact = vgradient(witness_vec, generator, zero_state(n), rots)*Trank
T[1:1<<n+1:1<<(2*n)] .= 0

function loss_func()
    p = zero_state(nqubits(witness_vec)) |> generator |> probs
    -sum(p'*T*p)
end

# gnum = num_gradient(loss_func, rots)

function train(generator::AbstractBlock, g_learning_rate::Real, niter::Int)
    rots = collect_rotblocks(generator)
    for i in 1:niter
        ggrad = vgradient(T*(zero_state(n) |> generator |> probs), generator, zero_state(n), rots)*2
        dispatch!(+, generator, -ggrad.*g_learning_rate)
        println("loss = $(loss_func())")
    end
end

train(generator, 0.1, 200)

using PyPlot
bar(0:1<<n-1, apply!(zero_state(n), generator)|>probs)

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
# Types in Julia need to be loaded in correct order for dependencies to be resolved.
# This utility script determines the correct order for including generated models.
#
# It parses each model and keeps a table of the other models it refers to.
# At the end it walks through the dependency tree and outputs include statements in the correct order.

using Dates

function typedeps(file::String)
    modulename = split(basename(file), ".")[1]
    contents = read(file, String)
    wrapped = """module $modulename
        $contents
    end"""
    typedeps(Meta.parse(wrapped))
end

function findtype(X::Expr, what::Symbol)
    (X.head === what) && (return X)
    for x in X.args
        isa(x, Expr) && (x.head === what) && return x
    end
    error("no $what in expression")
end

function isbasetype(DT)
    try
        T = eval(DT)
        (T <: Number) || (T <: String) || (T === Any) || (T === DateTime)
    catch
        false
    end
end

function pushdeps(deps, DT)
    if isa(DT, Expr) && DT.args[1] === :Vector
        pushdeps(deps, DT.args[2])
    elseif isa(DT, Expr) && DT.args[1] === :Dict
        pushdeps(deps, DT.args[2])
        pushdeps(deps, DT.args[3])
    elseif !isbasetype(DT)
        @info("    depends on $DT")
        push!(deps, DT)
    end
end

function typedeps(M::Expr)
    mod = findtype(M, :module)
    modblock = findtype(mod, :block)
    typ = findtype(modblock, :struct)
    typedecl = findtype(typ, :<:)
    typeblock = findtype(typ, :block)

    typename = typedecl.args[1]
    @assert typedecl.args[2] === :SwaggerModel

    deps = Vector{Symbol}()
    for d in typeblock.args
        if isa(d, Expr) && (d.head === :(::))
            nullable_type = d.args[2]
            if (nullable_type.head === :curly) && (nullable_type.args[1] === :Union) && (length(nullable_type.args) === 3) && (:Nothing in nullable_type.args)
                pushdeps(deps, first(filter(x->(x !== :Nothing) && (x !== :Union), nullable_type.args)))
            end
        end
    end
    (typename, deps)
end

srcdir(folder) = joinpath(folder, "src")
fullsrcpath(folder, file) = joinpath(srcdir(folder), file)

function gentypetree(folder::String)
    @info("reading $folder/src/model_*.jl")
    TT = Dict{Symbol, Vector{Symbol}}()
    TF = Dict{Symbol,String}()
    for file in readdir(srcdir(folder))
        if startswith(file, "model_")
            @info("parsing $file")
            typename, deps = typedeps(fullsrcpath(folder, file))
            TT[typename] = deps
            TF[typename] = file
        end
    end
    TT, TF
end

satisfied(typename, TT, generated) = (isempty(TT[typename]) || (typename in generated))

function gen(typename, TT, TF, generated, genstack, io=STDOUT)
    (typename in generated) && return
    if !satisfied(typename, TT, generated)
        push!(genstack, typename)
        for T in TT[typename]
            if T == typename
                @info("found recursive type use in $T")
            else
                if T in genstack
                    error("circular type references are not supported, found $T and $typename")
                end
                gen(T, TT, TF, generated, genstack, io)
            end
        end
    end
    println(io, "include(\"", TF[typename], "\")")
    push!(generated, typename)
    (typename in genstack) && delete!(genstack, typename)
    nothing
end

function genincludes(folder::String)
    generated = Set{Symbol}()
    genstack = Set{Symbol}()
    TT, TF = gentypetree(folder)

    open(fullsrcpath(folder, "modelincludes.jl"), "w") do inclfile
        for typename in keys(TT)
            gen(typename, TT, TF, generated, genstack, inclfile)
        end
    end

    nothing
end


text:
using Flux
using Zygote
using CUDA

abstract type NNStructure end

"""
    (nn::NNStructure)(x::AbstractArray{Float32,3})

Make NNStructure able to work with batches.
"""
(nn::NNStructure)(ts::AbstractTrajectoryState) = throw(ErrorException("missing function (::$(typeof(nn)))(::$(typeof(ts)))."))
function (nn::NNStructure)(x::AbstractVector{<:TabularTrajectoryState})
    batch_size = size(x, 3)
    qval = [nn(x[:, :, i]) for i in 1:batch_size]
    hcat(qval...)
end

function (nn::NNStructure)(x::AbstractVector{<:NonTabularTrajectoryState})
    qval = nn.(x)
    return hcat(qval...)
end

include("weighted_graph_gat.jl")
include("geometricflux.jl")
include("cpnn.jl")
include("fullfeaturedcpnn.jl")
include("variableoutputcpnn.jl")

struct ModelNotImplementedError{M} <: Exception
    m::M
    ModelNotImplementedError(m::M) where {M} = new{M}(m)
end

Base.showerror(io::IO, ie::ModelNotImplementedError) = print(io, "Model $(ie.m) not implemented.")

wears_mask(structure::NNStructure) = true
wears_mask(structure) = true # For simpler structures like Flux.Chain


text:
using ArgCheck

include("checks.jl")
include("perf.jl")


text:
using OpenQuantumBase, Test

# === matrix decomposition
v = 1.0 * σx + 2.0 * σy + 3.0 * σz
res = matrix_decompose(v, [σx, σy, σz])
@test isapprox(res, [1.0, 2.0, 3.0])
# === positivity test ===
r = rand(2)
m =
    r[1] * PauliVec[1][2] * PauliVec[1][2]' +
    r[2] * PauliVec[1][1] * PauliVec[1][1]'
@test check_positivity(m)
@test !check_positivity(σx)
@test_logs (:warn, "Input fails the numerical test for Hermitian matrix. Use the upper triangle to construct a new Hermitian matrix.") check_positivity(σ₊)
# == units conversion test ===
@test temperature_2_freq(1e3) ≈ 20.8366176361328 atol = 1e-4 rtol = 1e-4
@test freq_2_temperature(20) ≈ 959.8489324422699 atol = 1e-4 rtol = 1e-4
@test temperature_2_freq(1e3) ≈ 1 /  temperature_2_β(1e3) / 2 / π
@test β_2_temperature(0.47) ≈ 16.251564065921915
# === unitary test ===
u_res = exp(-1.0im * 5 * 0.5 * σx)
@test check_unitary(u_res)
@test !check_unitary([0 1; 0 0])
# === integration test ===
@test OpenQuantumBase.cpvagk((x) -> 1.0, 0, -1, 1)[1] == 0

# == Hamiltonian analysis ===
DH = DenseHamiltonian([(s) -> 1 - s, (s) -> s], [σx, σz], unit=:ħ)
# hfun(s) = (1-s)*real(σx)+ s*real(σz)
# dhfun(s) = -real(σx) + real(σz)
t = [0.0, 1.0]
states = [PauliVec[1][2], PauliVec[1][1]]
res = inst_population(t, states, DH, lvl=1:2)
@test isapprox(res, [[1.0, 0], [0.5, 0.5]])

SH = SparseHamiltonian(
    [(s) -> -(1 - s), (s) -> s],
    [standard_driver(2, sp=true), (0.1 * spσz ⊗ spσi + spσz ⊗ spσz)],
    unit=:ħ,
)
H_check = DenseHamiltonian(
    [(s) -> -(1 - s), (s) -> s],
    [standard_driver(2), (0.1 * σz ⊗ σi + σz ⊗ σz)],
    unit=:ħ,
)
# spdhfun(s) =
#    real(standard_driver(2, sp = true) + (0.1 * spσz ⊗ spσi + spσz ⊗ spσz))
interaction = [spσz ⊗ spσi, spσi ⊗ spσz]
spw, spv = eigen_decomp(SH, [0.5])
w, v = eigen_decomp(H_check, [0.5])
@test w ≈ spw atol = 1e-4
@test isapprox(spv[:, 1, 1], v[:, 1, 1], atol=1e-4) ||
      isapprox(spv[:, 1, 1], -v[:, 1, 1], atol=1e-4)
@test isapprox(spv[:, 2, 1], v[:, 2, 1], atol=1e-4) ||
      isapprox(spv[:, 2, 1], -v[:, 2, 1], atol=1e-4)
# == utility math functions ==
@test log_uniform(1, 10, 3) == [1, 10^0.5, 10]

v = sqrt.([0.4, 0.6])
ρ1 = v*v'
ρ2 = [0.5 0; 0 0.5]
ρ3 = ones(3, 3)/3
@test ρ1 == partial_trace(ρ1 ⊗ ρ2 ⊗ ρ2, [1])
@test ρ2 == partial_trace(ρ1 ⊗ ρ2 ⊗ ρ2, [2])
@test ρ3 ≈ partial_trace(ρ1⊗ρ2⊗ρ3, [2,2,3], [3])
@test_throws ArgumentError partial_trace(ρ1⊗ρ2⊗ρ3, [3,2,3], [3])
@test_throws ArgumentError partial_trace(rand(4,5), [2,2], [1])
@test purity(ρ1) ≈ 1
@test purity(ρ2) == 0.5
@test check_pure_state(ρ1)
@test !check_pure_state(ρ2)
@test !check_pure_state([0.4 0.5; 0.5 0.6])

ρ = PauliVec[1][1]*PauliVec[1][1]'
σ = PauliVec[3][1]*PauliVec[3][1]'
@test fidelity(ρ, σ) ≈ 0.5
@test fidelity(ρ, ρ) ≈ 1
@test check_density_matrix(ρ)
@test !check_density_matrix(σx)

w = [1,1,1,2,3,4,4,5]
@test OpenQuantumBase.find_degenerate(w) == [[1,2,3],[6,7]]
@test isempty(OpenQuantumBase.find_degenerate([1,2,3]))

gibbs = gibbs_state(σz, 12)
@test gibbs[2,2] ≈ 1/(1+exp(-temperature_2_β(12)*2))
@test gibbs[1,1] ≈ 1 - 1/(1+exp(-temperature_2_β(12)*2))

@test low_level_matrix(σz⊗σz+0.1σz⊗σi, 2) == [0.0+0.0im   0.0+0.0im   0.0+0.0im  0.0+0.0im
0.0+0.0im  -0.9+0.0im   0.0+0.0im  0.0+0.0im
0.0+0.0im   0.0+0.0im  -1.1+0.0im  0.0+0.0im
0.0+0.0im   0.0+0.0im   0.0+0.0im  0.0+0.0im]
@test_logs (:warn, "Subspace dimension bigger than total dimension.") low_level_matrix(σz⊗σz+0.1σz⊗σi, 5)


text:
using DataFrames
using CSV
using GLM
using Lathe

df_steam_t = DataFrame(CSV.File("../resources/data3.csv"))
coln_steam = Symbol[]
for i in string.(names(df_steam_t))
    i = replace(strip(i), " " => "_")
    i = replace(i, "-" => "")
    i = replace(i, "/" => "")
    i = replace(i, "(" => "")
    i = replace(i, ")" => "")
    i = replace(i, "%" => "")
    i = replace(i, "," => "")
    i = replace(i, "*" => "")
    i = replace(i, "+" => "")
    println(i)
    push!(coln_steam, Symbol(i))
end

rename!(df_steam_t, coln_steam)
#LP_Turbine_Power_kW
#Condenser_Duty_MMBTUhr
#DAC_Steam_Duty_MMBTUhr
#Total_ST_and_BOP_Auxiliary_Loads_kW
formulasSteamT = FormulaTerm[]
#
fm0st = @formula(LP_Turbine_Power_kW ~ Percent_of_Gas_Turbine_Electrical_Load)
push!(formulasSteamT, fm0st)
fm1st = @formula(Duty_of_CC_Reboiler_MMBTUhr ~ Percent_of_Gas_Turbine_Electrical_Load)
push!(formulasSteamT, fm1st)
fm2st = @formula(DAC_Steam_Duty_MMBTUhr ~ Percent_of_Gas_Turbine_Electrical_Load)
push!(formulasSteamT, fm2st)
fm3st = @formula(Total_ST_and_BOP_Auxiliary_Loads_kW ~ Percent_of_Gas_Turbine_Electrical_Load)
push!(formulasSteamT, fm3st)
coef_steam_df = DataFrame(y = String[], b = Float64[], a = Float64[])
for i in formulasSteamT
    linReg = lm(i, df_steam_t)
    y = string(terms(i)[1])
    b = coef(linReg)[1]
    a = coef(linReg)[2]
    push!(coef_steam_df, (y, b, a))
    println(linReg)
    println(deviance(linReg))
    println(coef(linReg))
end

CSV.write("steam_coeffs_v3.csv", coef_steam_df)

max_stars_repo_name: dthierry/princetonDacLti

max_issues_repo_name: dthierry/princetonDacLti

max_forks_repo_name: dthierry/princetonDacLti


text:
function adjacency_matrix(adj::AbstractMatrix, T::DataType=eltype(adj))
    m, n = size(adj)
    (m == n) || throw(DimensionMismatch("adjacency matrix is not a square matrix: ($m, $n)"))
    T.(adj)
end

function degrees(adj::AbstractMatrix; dir::Symbol=:out)
    if issymmetric(adj)
        d = vec(sum(adj, dims=1))
    else
        if dir == :out
            d = vec(sum(adj, dims=1))
        elseif dir == :in
            d = vec(sum(adj, dims=2))
        elseif dir == :both
            d = vec(sum(adj, dims=1)) + vec(sum(adj, dims=2))
        else
            throw(DomainError(dir, "invalid argument, only accept :in, :out and :both"))
        end
    end
    d
end

degrees(adj::AbstractMatrix, T::DataType; dir::Symbol=:out) = degrees(T.(adj); dir=dir)

"""
    degree_matrix(g[, T]; dir=:out)

Degree matrix of graph `g`. Return a matrix which contains degrees of each vertex in its diagonal.
The values other than diagonal are zeros.

# Arguments

- `g`: should be a adjacency matrix, `FeaturedGraph`, `SimpleGraph`, `SimpleDiGraph` (from Graphs)
    or `SimpleWeightedGraph`, `SimpleWeightedDiGraph` (from SimpleWeightedGraphs).
- `T`: result element type of degree vector; default is the element type of `g` (optional).
- `dir`: direction of degree; should be `:in`, `:out`, or `:both` (optional).

# Examples

```jldoctest
julia> using GraphLaplacians

julia> m = [0 1 1; 1 0 0; 1 0 0];

julia> GraphLaplacians.degree_matrix(m)
3×3 LinearAlgebra.Diagonal{Int64, Vector{Int64}}:
 2  ⋅  ⋅
 ⋅  1  ⋅
 ⋅  ⋅  1
```
"""
function degree_matrix(adj::AbstractMatrix, T::DataType=eltype(adj);
                       dir::Symbol=:out, squared::Bool=false, inverse::Bool=false)
    d = degrees(adj, T, dir=dir)
    squared && (d .= sqrt.(d))
    inverse && (d .= inv.(d); replace!(d, typemax(T)=>zero(T)))
    return Diagonal(T.(d))
end

"""
    laplacian_matrix(g[, T]; dir=:out)

Laplacian matrix of graph `g`.

# Arguments

- `g`: should be a adjacency matrix, `FeaturedGraph`, `SimpleGraph`, `SimpleDiGraph` (from Graphs)
    or `SimpleWeightedGraph`, `SimpleWeightedDiGraph` (from SimpleWeightedGraphs).
- `T`: result element type of degree vector; default is the element type of `g` (optional).
- `dir`: direction of degree; should be `:in`, `:out`, or `:both` (optional).
"""
Graphs.laplacian_matrix(adj::AbstractMatrix, T::DataType=eltype(adj); dir::Symbol=:out) =
    degree_matrix(adj, T, dir=dir) - SparseMatrixCSC(T.(adj))

"""
    normalized_laplacian(g[, T]; dir=:both, selfloop=false)

Normalized Laplacian matrix of graph `g`.

# Arguments

- `g`: should be a adjacency matrix, `FeaturedGraph`, `SimpleGraph`, `SimpleDiGraph` (from Graphs)
    or `SimpleWeightedGraph`, `SimpleWeightedDiGraph` (from SimpleWeightedGraphs).
- `T`: result element type of degree vector; default is the element type of `g` (optional).
- `selfloop`: adding self loop while calculating the matrix (optional).
- `dir`: direction of graph; should be `:in` or `:out` (optional).
"""
function normalized_laplacian(adj::AbstractMatrix, T::DataType=eltype(adj);
                              dir::Symbol=:both, selfloop::Bool=false)
    if dir == :both
        selfloop && (adj += I)
        inv_sqrtD = degree_matrix(adj, T, dir=:both, squared=true, inverse=true)
        return T.(I - inv_sqrtD * adj * inv_sqrtD)
    else
        return T.(I - degree_matrix(adj, T, dir=dir, inverse=true) * adj)
    end
end

@doc raw"""
    scaled_laplacian(g[, T])

Scaled Laplacien matrix of graph `g`,
defined as ``\hat{L} = \frac{2}{\lambda_{max}} L - I`` where ``L`` is the normalized Laplacian matrix.

# Arguments

- `g`: should be a adjacency matrix, `FeaturedGraph`, `SimpleGraph`, `SimpleDiGraph` (from Graphs)
    or `SimpleWeightedGraph`, `SimpleWeightedDiGraph` (from SimpleWeightedGraphs).
- `T`: result element type of degree vector; default is the element type of `g` (optional).
"""
function scaled_laplacian(adj::AbstractMatrix, T::DataType=eltype(adj))
    @assert issymmetric(adj) "scaled_laplacian only works with symmetric matrices"
    E = eigen(Symmetric(Array(adj))).values
    T(2. / maximum(E)) * normalized_laplacian(adj, T) - I
end

"""
    random_walk_laplacian(g[, T]; dir=:out)

Random walk normalized Laplacian matrix of graph `g`.

# Arguments

- `g`: should be a adjacency matrix, `FeaturedGraph`, `SimpleGraph`, `SimpleDiGraph` (from Graphs)
    or `SimpleWeightedGraph`, `SimpleWeightedDiGraph` (from SimpleWeightedGraphs).
- `T`: result element type of degree vector; default is the element type of `g` (optional).
- `dir`: direction of degree; should be `:in`, `:out`, or `:both` (optional).
"""
function random_walk_laplacian(adj::AbstractMatrix, T::DataType=eltype(adj); dir::Symbol=:out)
    d = degrees(adj, dir=dir)
    inv_d = 1 ./ d
    replace!(inv_d, typemax(float(T)) => zero(float(T)))  # avoid degree to be zero
    P = Diagonal(inv_d) * adj
    SparseMatrixCSC(T.(I - P))
end

"""
    signless_laplacian(g[, T]; dir=:out)

Signless Laplacian matrix of graph `g`.

# Arguments

- `g`: should be a adjacency matrix, `FeaturedGraph`, `SimpleGraph`, `SimpleDiGraph` (from Graphs)
    or `SimpleWeightedGraph`, `SimpleWeightedDiGraph` (from SimpleWeightedGraphs).
- `T`: result element type of degree vector; default is the element type of `g` (optional).
- `dir`: direction of degree; should be `:in`, `:out`, or `:both` (optional).
"""
signless_laplacian(adj::AbstractMatrix, T::DataType=eltype(adj); dir::Symbol=:out) =
    degree_matrix(adj, T, dir=dir) + SparseMatrixCSC(T.(adj))


text:
using Test, Plots, JuMP, Ipopt, Clp, LinearAlgebra

u_bound = 1.5
model, obj = TrajectoryOptimization.Dynamics.pendulum!
obj = TrajectoryOptimization.ConstrainedObjective(obj, u_min=-u_bound, u_max=u_bound)

opts = TrajectoryOptimization.SolverOptions()
opts.verbose = true
solver = TrajectoryOptimization.Solver(model,obj,dt=0.1,opts=opts)
U = zeros(model.m,solver.N-1)
results_c, = TrajectoryOptimization.solve(solver, U)
max_c = TrajectoryOptimization.max_violation(results_c)

plot(to_array(results_c.U)')

# costfun = solver.obj.cost
# x = results_c.X[1]
# u = results_c.U[1]
# expansion = taylor_expansion(costfun,x,u)
# Q,R,H,q,r = taylor_expansion(costfun,x,u)
# H

plot(to_array(results_c.λ)[:,1:solver.N-1]')

function λ_second_order_update!(results::SolverIterResults,solver::Solver,verbose::Bool=false)
    n = solver.model.n
    m = solver.model.m
    nm = n+m
    N = solver.N
    p,pI,pE = get_num_constraints(solver)
    p_N,pI_N,pE_N = get_num_terminal_constraints(solver)

    Nz = n*N + m*(N-1) # number of decision variables x, u
    Nu = m*(N-1) # number of control decision variables u

    Np = p*(N-1) + p_N # number of constraints

    Mz = N + N-1 # number of decision vectors
    Mu = N-1 # number of control decision vectors

    B̄ = zeros(Nz,Nu)
    Ā = zeros(Nz,n)

    Q̄ = zeros(Nz,Nz)
    q̄ = zeros(Nz)

    C̄ = zeros(Np,Nz)
    c̄ = zeros(Np)

    λ_tmp = zeros(Np)
    idx_inequality = zeros(Bool,Np)
    active_set = zeros(Bool,Np)

    x0 = solver.obj.x0
    costfun = solver.obj.cost

    for k = 1:N
        x = results.X[k]

        if k != N
            u = results.U[k]
            Q,R,H,q,r = taylor_expansion(costfun,x,u)
        else
            Qf,qf = taylor_expansion(costfun,x)
        end

        # Indices
        idx = ((k-1)*nm + 1):k*nm
        idx2 = ((k-1)*p + 1):k*p
        idx3 = ((k-1)*nm + 1):((k-1)*nm + n)
        idx4 = ((k-1)*nm + n + 1):k*nm
        idx5 = ((k-1)*m + 1):k*m
        idx6 = ((k-1)*p + 1):((k-1)*p + pI)

        # Calculate Ā
        k == 1 ? Ā[idx3,1:n] = 1.0*Matrix(I,n,n) : Ā[idx3,1:n] = prod(results.fdx[1:k-1])

        # Calculate B̄
        if k > 1
            for j = 1:k-1
                idx7 = ((j-1)*m + 1):j*m
                j == k-1 ? B̄[idx3,idx7] = results.fdu[j][1:n,1:m] : B̄[idx3,idx7] = prod(results.fdx[j+1:(k-1)])*results.fdu[j][1:n,1:m]
            end
        end

        # Calculate Q̄, q̄, C̄, c̄
        if k != N
            Q̄[idx,idx] = [Q H';H R]
            q̄[idx] = [q;r]

            G = results.Cx[k]
            H = results.Cu[k]
            C̄[idx2,idx] = [G H]
            c̄[idx2] = results.C[k]

            λ_tmp[idx2] = results.λ[k]
            idx_inequality[idx6] .= true
            active_set[idx2] = results.active_set[k]

            B̄[idx4,idx5] = 1.0*Matrix(I,m,m)
        else
            idx = ((k-1)*nm + 1):Nz
            Q̄[idx,idx] = Qf
            q̄[idx] = qf

            idx2 = ((k-1)*p + 1):Np
            C̄[idx2,idx] = results.Cx[N]
            c̄[idx2] = results.C[N]

            idx6 = ((k-1)*p + 1):((k-1)*p + pI_N)

            λ_tmp[idx2] = results.λ[N]
            idx_inequality[idx6] .= true
            active_set[idx2] .= results.active_set[N]
        end

    end
    N_active_set = sum(active_set)

    # if verbose
    #     # Tests
    #     @test Ā[1:n,:] == Matrix(I,n,n)
    #     k = 3
    #     @test Ā[(k-1)*nm+1:(k-1)*nm+n,:] == prod(results.fdx[1:2])
    #     @test results.fdx[1]*results.fdx[2] == prod(results.fdx[1:2])
    #
    #     k = 7
    #     @test Ā[(k-1)*nm+1:(k-1)*nm+n,:] == prod(results.fdx[1:k-1])
    #     @test Ā[(N-1)*nm+1:Nz,:] == prod(results.fdx[1:N-1])
    #
    #     @test B̄[1:n,1:n] == zeros(n,n)
    #     @test B̄[n+1:nm,1:m] == 1.0*Matrix(I,m,m)
    #     @test B̄[nm+1:nm+n,1:m] == results.fdu[1][1:n,1:m]
    #     @test B̄[nm+n+1:2*nm,m+1:2*m] == 1.0*Matrix(I,m,m)
    #     @test B̄[(N-1)*nm+1:Nz,1:m] == prod(results.fdx[2:N-1])*results.fdu[1]
    #     @test B̄[(N-1)*nm+1:Nz,m+1:2*m] == prod(results.fdx[3:N-1])*results.fdu[2]
    #     @test B̄[(N-1)*nm+1:Nz,(N-2)*m+1:(N-1)*m] == results.fdu[N-1]
    #
    #     k = 1
    #     Q,R,H,q,r = taylor_expansion(costfun,results.X[k],results.U[k])
    #     @test Q̄[1:nm,1:nm] == [Q H'; H R]
    #     @test q̄[1:nm] == [q;r]
    #
    #     k = 13
    #     Q,R,H,q,r = taylor_expansion(costfun,results.X[k],results.U[k])
    #     @test Q̄[(k-1)*nm+1:k*nm,(k-1)*nm+1:k*nm] == [Q H'; H R]
    #     @test q̄[(k-1)*nm+1:k*nm] == [q;r]
    #
    #     k = N
    #     Qf,qf = taylor_expansion(costfun,results.X[k])
    #     @test Q̄[(k-1)*nm+1:Nz,(k-1)*nm+1:Nz] == Qf
    #     @test q̄[(k-1)*nm+1:Nz] == qf
    #
    #     @test C̄[1:p,1:nm] == [results.Cx[1] results.Cu[1]]
    #
    #     k = 9
    #     @test C̄[(k-1)*p+1:k*p,(k-1)*nm+1:k*nm] == [results.Cx[k] results.Cu[k]]
    #     @test C̄[(N-1)*p+1:Np,(N-1)*nm+1:Nz] == results.Cx[N]
    #
    #     @test c̄[1:p] == results.C[1]
    #
    #     k = 17
    #     @test c̄[(k-1)*p+1:k*p] == results.C[k]
    #     @test c̄[(N-1)*p+1:Np] == results.C[N]
    #
    #     @test all(idx_inequality[1:pI] .== true)
    #     @test all(idx_inequality[pI+1:p] .== false)
    #     k = 12
    #
    #     @test all(idx_inequality[(k-1)*p+1:(k-1)*p+pI] .== true)
    #     idx_inequality[(k-1)*p+1:(k-1)*p+pI]
    #     @test all(idx_inequality[(N-1)*p+1:(N-1)*pI_N] .== true)
    # end
    ū = -(B̄'*Q̄*B̄)\(B̄'*(q̄ + Q̄*Ā*x0) + B̄'*C̄'*λ_tmp)

    # println(any(isnan.(q̄)))
    # println(any(isnan.(Ā)))
    # println(any(isnan.(λ_tmp)))
    # println(any(isnan.(ū)))
    # println(any(isnan.(Ā)))
    # println(any(isnan.(B̄)))
    # println(any(isnan.(C̄)))
    # println(any(isnan.(c̄)))
    P = B̄'*Q̄*B̄
    M = q̄ + Q̄*Ā*x0

    # println(rank(P))
    # println(size(P))
    # println(λ_tmp)

    L = 0.5*ū'*P*ū + M'*B̄*ū + λ_tmp'*(C̄*B̄*ū + C̄*Ā*x0 + c̄)

    D = 0.5*M'*B̄*inv(P')*B̄'*M + 0.5*M'*B̄*inv(P')*B̄'*C̄'*λ_tmp + 0.5*λ_tmp'*C̄*B̄*inv(P')*B̄'*M + 0.5*λ_tmp'*C̄*B̄*inv(P')*B̄'*C̄'*λ_tmp - M'*B̄*inv(P)*B̄'*M - M'*B̄*inv(P)*B̄'*C̄'*λ_tmp - λ_tmp'*C̄*B̄*inv(P)*B̄'*M - λ_tmp'*C̄*B̄*inv(P)*B̄'*C̄'*λ_tmp + λ_tmp'*C̄*Ā*x0 + λ_tmp'*c̄
    Q_dual = C̄*B̄*inv(P')*B̄'*C̄' - 2*C̄*B̄*inv(P)*B̄'*C̄'
    q_dual = M'*B̄*inv(P')*B̄'*C̄' - 2*M'*B̄*inv(P)*B̄'*C̄' + x0'*Ā'*C̄' + c̄'
    qq_dual = 0.5*M'*B̄*inv(P')*B̄'*M - M'*B̄*inv(P)*B̄'*M

    DD = 0.5*λ_tmp'*Q_dual*λ_tmp + q_dual*λ_tmp + qq_dual

    @test isapprox(D,L)
    @test isapprox(DD,L)

    # solve QP
    m = JuMP.Model(solver=IpoptSolver(print_level=0))
    # m = JuMP.Model(solver=ClpSolver())

    # @variable(m, u[1:Nu])
    #
    # @objective(m, Min, 0.5*u'*P*u + M'*B̄*u)
    # @constraint(m, con, C̄*B̄*u + C̄*Ā*x0 + c̄ .== 0.)

    # @variable(m, λ[1:Np])
    N_active_set = sum(active_set)

    @variable(m, λ[1:N_active_set])

    # @objective(m, Min, λ'*λ)
    # @constraint(m, con, Q_dual[active_set,active_set]*λ + q_dual[active_set] .== 0.)
    #
    @objective(m, Max, 0.5*λ'*Q_dual[active_set,active_set]*λ + q_dual[active_set]'*λ)
    @constraint(m, con2, λ[idx_inequality[active_set]] .>= 0)
    #
    # # print(m)
    #
    status = JuMP.solve(m)
    #
    if status != :Optimal
        error("QP failed")
    end

    # Solution
    # println("Objective value: ", JuMP.getobjectivevalue(m))
    # println("λ = ", JuMP.getvalue(λ))

    λ_tmp[active_set] = JuMP.getvalue(λ)

    for k = 1:N
        if k != N
            idx = ((k-1)*p + 1):k*p
            results.λ[k] = λ_tmp[idx]#JuMP.getvalue(λ[idx])
            # results.λ[k][1:pI] = max.(0.0,results.λ[k][1:pI])
        else
            idx = ((N-1)*p + 1):Np
            results.λ[k] = λ_tmp[idx]#JuMP.getvalue(λ[idx])
            # results.λ[k][1:pI_N] = max.(0.0,results.λ[k][1:pI_N])
        end

        # tolerance check
        results.λ[k][abs.(results.λ[k]) .< 1e-8] .= 0.
    end
end

λ_second_order_update!(results_c,solver,true)

plot!(to_array(results_c.λ)[:,1:solver.N-1]')

solver = TrajectoryOptimization.Solver(model,obj,dt=0.01,opts=opts)
U = zeros(model.m,solver.N-1)
solver.opts.verbose = true
solver.opts.cost_tolerance = 1e-6
solver.opts.cost_tolerance_intermediate = 1e-6
solver.opts.gradient_norm_tolerance = 1e-6
solver.opts.gradient_norm_tolerance_intermediate = 1e-6
solver.opts.constraint_tolerance = 1e-6
solver.opts.use_second_order_dual_update = true
@time results_c,stats = TrajectoryOptimization.solve(solver, U)

plot(to_array(results_c.λ)[:,1:solver.N-1]')

plot(to_array(results_c.U)[:,1:solver.N-1]')

plot!(log.(stats["cost"] .+ 1000))