text:
#
# @description Implements the ECB (Electronic Codebook) mode of operation algorithm.
#
# @author Alejandro Claro (alejandro.claro@gmail.com)
#
# Copyright 2017 All rights reserved.
# Use of this source code is governed by a MIT-style license that can be found in the LICENSE file.
#
export EcbModeCipher
export block_size, key_size, reset!, encrypt_next_blocks!, encrypt_last_block!, decrypt_next_blocks!, decrypt_last_block!

# @description Defines the CBC cipher mode data struture.
struct EcbModeCipher <: MultiBlockCipher
  block_cipher::BlockCipher
  padder::PaddingMethod
  partial_block::Vector{UInt8}

  # @description Initialize the ECB cipher.
  #
  # @param {BlockCipher}   block_cipher The block cipher.
  # @param {PaddingMethod} padder       The padding method.
  function CbcModeCipher(block_cipher::BlockCipher, padder::PaddingMethod)
    return new(block_cipher, padder, UInt8[]))
  end
end

# @description Gets the cipher block size.
#
# @param {EcbModeCipher} self The cipher data struture.
#
# @return {Integer} The block size in bytes.
function block_size(self::EcbModeCipher)
  return block_size(self.block_cipher)
end

# @description Gets the cipher key size.
#
# @param {EcbModeCipher} self The cipher data struture.
#
# @return {Integer} The key size in bytes.
function key_size(self::EcbModeCipher)
  return key_size(self.block_cipher)
end

# @description Resets the state of the CBC data structure.
#
# @param {EcbModeCipher} self The hash algorithm.
#
# @return The reference to the resetted cipher data structure.
function reset!(self::EcbModeCipher)
  empty!(self.partial_block)
  return self
end

# @description Continues a multiple-part encryption operation processing another data part.
#
# @param {EcbModeCipher} self The cipher data struture.
# @param {Vector{UInt8}} data The message part to encrypt.
#
# @return {Vector{UInt8}} The resulting ciphertext blocks. It can be a empty vector if the input data is too short to
# result in a new block.
function encrypt_next_blocks!(self::EcbModeCipher, data::Vector{UInt8})
  data   = vcat(self.partial_block, data)
  step   = block_size(self)
  len    = length(data)
  result = UInt8[]

  for i in 1 : step : (len - 1)
    ciphertext = encrypt(self.block_cipher, data[i : (i + step - 1)])
    append!(result, ciphertext)
  end

  empty!(self.partial_block)
  append!(self.partial_block, data[(end - mod(len, step) + 1) : end])

  return result
end

# @description Encrypts the given last plaintext block.
#
# @param {EcbModeCipher} self  The cipher data struture.
# @param {Vector{UInt8}} block The last block from the message to encrypt.
#
# @return {Vector{UInt8}} The resulting last ciphertext blocks. This could be more than the size of a block.
function encrypt_last_block!(self::EcbModeCipher, block::Vector{UInt8})
  result       = encrypt_next_blocks!(self, block)
  padded_block = pad(self.padder, self.partial_block, block_size(self))
  ciphertext   = encrypt(self.block_cipher, padded_block)

  append!(result, ciphertext)
  reset!(self)

  return result
end

# @description Continues a multiple-part decryption operation processing another data part.
#
# @param {EcbModeCipher} self  The cipher data struture.
# @param {Vector{UInt8}} data  The message part to decrypt.
#
# @return {Vector{UInt8}} The resulting plaintext blocks. It can be a empty vector if the input data is too short to
# result in a new block.
function decrypt_next_blocks!(self::EcbModeCipher, data::Vector{UInt8})
  data   = vcat(self.partial_block, data)
  step   = block_size(self)
  len    = length(data)
  result = UInt8[]

  for i in 1 : step : (len - 1)
    ciphertext = data[i : (i + step - 1)]
    plaintext  = decrypt(self.block_cipher, ciphertext)
    append!(result, plaintext)
  end

  empty!(self.partial_block)
  append!(self.partial_block, data[(end - mod(len, step) + 1) : end])

  return result
end

# @description Decrypts the given last ciphertext block.
#
# @param {EcbModeCipher} self  The cipher data struture.
# @param {Vector{UInt8}} block The last block from the message to decrypt.
#
# @return {Vector{UInt8}} The resulting last plaintext block. This could be less than the size of a block.
function decrypt_last_block!(self::EcbModeCipher, block::Vector{UInt8})
  result = decrypt_next_blocks!(self, block)

  if !isempty(self.partial_block)
    throw(BlockSizeError())
  end

  result = unpad(self.padder, result)
  reset!(self)

  return result
end


text:
using Images: mapwindow, Fill, label_components
using StatsBase: counts

function part1(input)
    heights = reduce(hcat, (parse.(Int, split(line, "")) for line in eachline(input)))
    return mapwindow(heights, (3, 3), border = Fill(9)) do x
        x[2,2] < min(x[1,2], x[3,2], x[2,1], x[2,3]) ? 1 + x[2,2] : 0
    end |> sum
end

function part2(input)
    heights = reduce(hcat, (parse.(Int, split(line, "")) for line in eachline(input)))
    return prod(sort(counts(label_components(heights .< 9))[2:end], rev=true)[1:3])
end

ext: jl


text:
module CAM
using MLStyle
using MLStyle.Infras: @format, mangle
using CamJulia.Runtime
using CamJulia

@use Enum # replace `S() => ...` with `S => ...`

# held here until merged into MLStyle master.
def_pattern(CAM,
    predicate = (@λ begin
        :[$_ for $_ in $_] -> true
        _ -> false
    end),
    rewrite = (tag, case, mod) -> begin
        @match case begin
            :[$expr for $iter in $seq] => begin
                function (body)
                    iter_var = mangle(mod)
                    produced_elt_var = mangle(mod)
                    inner_test_var = mangle(mod)
                    seq_var = mangle(mod)
                    decons_elt = mk_pattern(iter_var, expr, mod)(iter)
                    decons_seq = mk_pattern(seq_var, seq, mod)(body)
                    @format [
                        seq_var, iter_var, inner_test_var, produced_elt_var,
                        body, tag, decons_elt, decons_seq, push!] quote
                        if tag isa Vector
                            let seq_var = [], inner_test_var = true
                                for iter_var in tag
                                    produced_elt_var = decons_elt
                                    if produced_elt_var === failed
                                        inner_test_var = false
                                        break
                                    end
                                    push!(seq_var, produced_elt_var)
                                end
                                if inner_test_var
                                    decons_seq
                                else
                                    failed
                                end
                            end
                        else
                            failed
                        end
                    end
                end
            end
        end

    end
)

export IR, Let, LetRec, If, While, Mutate, Fun
export App, Var, Block, Join, Proj
export BigIntConst, IntConst, DoubleConst
export StrConst, ChConst, BoolConst, NilConst
export Internal, Located, SymConst

@data IR begin
    Let(String, IR, IR)
    LetRec(Vector{Tuple{String, IR}}, IR)
    If(IR, IR, IR)
    While(IR, IR)
    Mutate(String, IR)
    Fun(Vector{String}, IR)
    App(IR, Vector{IR})
    Var(String)
    Block(Vector{IR})
    Join(Vector{IR})
    Proj(IR, IR)

    BigIntConst(IdrisIntegerT)
    IntConst(Int)
    DoubleConst(Float64)
    StrConst(String)
    ChConst(Char)
    BoolConst(Bool)
    SymConst(String)
    NilConst()

    Internal(String)
    Located(fname::Union{Nothing, String}, line::Int, column::Int, value::IR)
end

@active Sym(x :: String) begin
    Symbol(x)
end

@active Julia(x) begin
    ir_to_julia(x)
end

export ir_to_julia

@active LetToDef(tp) begin
    @match tp begin
        (Sym(fname), Fun([Sym(arg) for arg in args], Julia(body))) =>
        :($fname($(args...), ) = $body)
    end
end

@generated function projection(x, i :: Int)
    x <: Tuple ? :(x[i]) : :(if i === 1; x else throw("internal runtime error") end)
end

ir_to_julia(ir::IR) =
    @match ir begin
        SymConst(Sym(sym)) => QuoteNode(sym)
        Var(Sym(sym)) => sym
        # Let(Sym(fn_name), Fun([Sym(arg) for arg in args], Julia(fn_body)), Julia(body)) =>
        #     :(let $fn_name($(args...), ) = $fn_body; $body end)
        Let(Sym(s), Julia(value), Julia(body)) =>
            :(let $s = $value; $body end)

        # letrec: bindings must be all functions
        LetRec([LetToDef(tp) for tp in seq], Julia(value)) =>
            Expr(:block, seq..., value)
        If(Julia(cond), Julia(iftrue), Julia(iffalse)) =>
            :(if $cond; $iftrue else $iffalse end)
        While(Julia(cond), Julia(body)) =>
            :(while $cond; $body end)
        Mutate(Sym(s), Julia(value)) =>
            :($s = $value)
        Fun([Sym(s) for s in args], Julia(body)) =>
            :(function ($(args...), ) $body end)

        App(Internal("builtin-get_module"), [StrConst(c)]) =>
            let paths = map(Symbol, split(c, ".")),
                imp = @match paths begin
                    [] => throw("no module accessing path given! requires a path like `a.b.c`!")
                    [x] => x
                    _ => Expr(:., paths)
                end
                quote
                    @eval import $imp
                    $imp
                end
            end

        App(Julia(fn), [Julia(arg) for arg in args]) =>
            :($fn($(args...), ))

        Block([Julia(expr) for expr in seq]) =>
            :(begin $(seq...) end)
        Join([Julia(elt) for elt in elts]) =>
            :($(elts...), )
        Proj(Julia(major), Julia(ith)) =>
            :($projection($major, $ith + 1))
            # :($major[$ith + 1])

        StrConst(c) => foldr(c, init=CamJulia.IdrisList.IdrisNil{Char}()) do each, prev
            each ^ prev
        end

        BigIntConst(c) || IntConst(c) || DoubleConst(c) ||
        # TODO: str const needs to be unescaped, e.g., '\\a' -> '\a'
        ChConst(c) || BoolConst(c) => c

        NilConst => :nothing # kind of tricky for Julia use :nothing to represent nothing in ASTs

        Internal(s) => rt_support[s]
        Located(fname, lineno, _, Julia(value)) => # colno doesn't have an effect in current julia.
            let lineno = LineNumberNode(lineno, fname)
                :(begin $lineno; $value end)
            end
        a => throw(a)
    end
end

max_stars_repo_name: thautwarm/idris-cam

max_issues_repo_name: thautwarm/idris-cam

max_forks_repo_name: thautwarm/idris-cam


text:
for i = 0 : (x = 97; 0) end
println(x - 55)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
include((@__DIR__)*"/../examples/optimal_control.jl")

@test termination_status(graph) == MOI.LOCALLY_SOLVED || termination_status(graph) == MOI.ITERATION_LIMIT


text:
# Julia wrapper for header: xc.h
# Automatically generated using Clang.jl

function xc_reference()
    ccall((:xc_reference, libxc), Cstring, ())
end

function xc_reference_doi()
    ccall((:xc_reference_doi, libxc), Cstring, ())
end

function xc_version(major, minor, micro)
    ccall((:xc_version, libxc), Cvoid, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), major, minor, micro)
end

function xc_version_string()
    ccall((:xc_version_string, libxc), Cstring, ())
end

function xc_func_reference_get_ref(reference)
    ccall((:xc_func_reference_get_ref, libxc), Cstring, (Ptr{func_reference_type},), reference)
end

function xc_func_reference_get_doi(reference)
    ccall((:xc_func_reference_get_doi, libxc), Cstring, (Ptr{func_reference_type},), reference)
end

function xc_func_reference_get_bibtex(reference)
    ccall((:xc_func_reference_get_bibtex, libxc), Cstring, (Ptr{func_reference_type},), reference)
end

function xc_func_info_get_number(info)
    ccall((:xc_func_info_get_number, libxc), Cint, (Ptr{xc_func_info_type},), info)
end

function xc_func_info_get_kind(info)
    ccall((:xc_func_info_get_kind, libxc), Cint, (Ptr{xc_func_info_type},), info)
end

function xc_func_info_get_name(info)
    ccall((:xc_func_info_get_name, libxc), Cstring, (Ptr{xc_func_info_type},), info)
end

function xc_func_info_get_family(info)
    ccall((:xc_func_info_get_family, libxc), Cint, (Ptr{xc_func_info_type},), info)
end

function xc_func_info_get_flags(info)
    ccall((:xc_func_info_get_flags, libxc), Cint, (Ptr{xc_func_info_type},), info)
end

function xc_func_info_get_references(info, number)
    ccall((:xc_func_info_get_references, libxc), Ptr{func_reference_type}, (Ptr{xc_func_info_type}, Cint), info, number)
end

function xc_func_info_get_n_ext_params(info)
    ccall((:xc_func_info_get_n_ext_params, libxc), Cint, (Ptr{xc_func_info_type},), info)
end

function xc_func_info_get_ext_params_name(p, number)
    ccall((:xc_func_info_get_ext_params_name, libxc), Cstring, (Ptr{xc_func_info_type}, Cint), p, number)
end

function xc_func_info_get_ext_params_description(info, number)
    ccall((:xc_func_info_get_ext_params_description, libxc), Cstring, (Ptr{xc_func_info_type}, Cint), info, number)
end

function xc_func_info_get_ext_params_default_value(info, number)
    ccall((:xc_func_info_get_ext_params_default_value, libxc), Cdouble, (Ptr{xc_func_info_type}, Cint), info, number)
end

function xc_functional_get_number(name)
    ccall((:xc_functional_get_number, libxc), Cint, (Cstring,), name)
end

function xc_functional_get_name(number)
    ccall((:xc_functional_get_name, libxc), Cstring, (Cint,), number)
end

function xc_family_from_id(id, family, number)
    ccall((:xc_family_from_id, libxc), Cint, (Cint, Ptr{Cint}, Ptr{Cint}), id, family, number)
end

function xc_number_of_functionals()
    ccall((:xc_number_of_functionals, libxc), Cint, ())
end

function xc_maximum_name_length()
    ccall((:xc_maximum_name_length, libxc), Cint, ())
end

function xc_available_functional_numbers(list)
    ccall((:xc_available_functional_numbers, libxc), Cvoid, (Ptr{Cint},), list)
end

function xc_available_functional_numbers_by_name(list)
    ccall((:xc_available_functional_numbers_by_name, libxc), Cvoid, (Ptr{Cint},), list)
end

function xc_available_functional_names(list)
    ccall((:xc_available_functional_names, libxc), Cvoid, (Ptr{Cstring},), list)
end

function xc_func_alloc()
    ccall((:xc_func_alloc, libxc), Ptr{xc_func_type}, ())
end

function xc_func_init(p, functional, nspin)
    ccall((:xc_func_init, libxc), Cint, (Ptr{xc_func_type}, Cint, Cint), p, functional, nspin)
end

function xc_func_end(p)
    ccall((:xc_func_end, libxc), Cvoid, (Ptr{xc_func_type},), p)
end

function xc_func_free(p)
    ccall((:xc_func_free, libxc), Cvoid, (Ptr{xc_func_type},), p)
end

function xc_func_get_info(p)
    ccall((:xc_func_get_info, libxc), Ptr{xc_func_info_type}, (Ptr{xc_func_type},), p)
end

function xc_func_set_dens_threshold(p, t_dens)
    ccall((:xc_func_set_dens_threshold, libxc), Cvoid, (Ptr{xc_func_type}, Cdouble), p, t_dens)
end

function xc_func_set_zeta_threshold(p, t_zeta)
    ccall((:xc_func_set_zeta_threshold, libxc), Cvoid, (Ptr{xc_func_type}, Cdouble), p, t_zeta)
end

function xc_func_set_sigma_threshold(p, t_sigma)
    ccall((:xc_func_set_sigma_threshold, libxc), Cvoid, (Ptr{xc_func_type}, Cdouble), p, t_sigma)
end

function xc_func_set_tau_threshold(p, t_tau)
    ccall((:xc_func_set_tau_threshold, libxc), Cvoid, (Ptr{xc_func_type}, Cdouble), p, t_tau)
end

function xc_func_set_ext_params(p, ext_params)
    ccall((:xc_func_set_ext_params, libxc), Cvoid, (Ptr{xc_func_type}, Ptr{Cdouble}), p, ext_params)
end

function xc_func_set_ext_params_name(p, name, par)
    ccall((:xc_func_set_ext_params_name, libxc), Cvoid, (Ptr{xc_func_type}, Cstring, Cdouble), p, name, par)
end

function xc_lda(p, np, rho, zk, vrho, v2rho2, v3rho3, v4rho4)
    ccall((:xc_lda, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, zk, vrho, v2rho2, v3rho3, v4rho4)
end

function xc_gga(p, np, rho, sigma, zk, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2, v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3, v4rho4, v4rho3sigma, v4rho2sigma2, v4rhosigma3, v4sigma4)
    ccall((:xc_gga, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, zk, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2, v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3, v4rho4, v4rho3sigma, v4rho2sigma2, v4rhosigma3, v4sigma4)
end

function xc_mgga(p, np, rho, sigma, lapl_rho, tau, zk, vrho, vsigma, vlapl, vtau, v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, v2lapl2, v2lapltau, v2tau2, v3rho3, v3rho2sigma, v3rho2lapl, v3rho2tau, v3rhosigma2, v3rhosigmalapl, v3rhosigmatau, v3rholapl2, v3rholapltau, v3rhotau2, v3sigma3, v3sigma2lapl, v3sigma2tau, v3sigmalapl2, v3sigmalapltau, v3sigmatau2, v3lapl3, v3lapl2tau, v3lapltau2, v3tau3, v4rho4, v4rho3sigma, v4rho3lapl, v4rho3tau, v4rho2sigma2, v4rho2sigmalapl, v4rho2sigmatau, v4rho2lapl2, v4rho2lapltau, v4rho2tau2, v4rhosigma3, v4rhosigma2lapl, v4rhosigma2tau, v4rhosigmalapl2, v4rhosigmalapltau, v4rhosigmatau2, v4rholapl3, v4rholapl2tau, v4rholapltau2, v4rhotau3, v4sigma4, v4sigma3lapl, v4sigma3tau, v4sigma2lapl2, v4sigma2lapltau, v4sigma2tau2, v4sigmalapl3, v4sigmalapl2tau, v4sigmalapltau2, v4sigmatau3, v4lapl4, v4lapl3tau, v4lapl2tau2, v4lapltau3, v4tau4)
    ccall((:xc_mgga, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, lapl_rho, tau, zk, vrho, vsigma, vlapl, vtau, v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, v2lapl2, v2lapltau, v2tau2, v3rho3, v3rho2sigma, v3rho2lapl, v3rho2tau, v3rhosigma2, v3rhosigmalapl, v3rhosigmatau, v3rholapl2, v3rholapltau, v3rhotau2, v3sigma3, v3sigma2lapl, v3sigma2tau, v3sigmalapl2, v3sigmalapltau, v3sigmatau2, v3lapl3, v3lapl2tau, v3lapltau2, v3tau3, v4rho4, v4rho3sigma, v4rho3lapl, v4rho3tau, v4rho2sigma2, v4rho2sigmalapl, v4rho2sigmatau, v4rho2lapl2, v4rho2lapltau, v4rho2tau2, v4rhosigma3, v4rhosigma2lapl, v4rhosigma2tau, v4rhosigmalapl2, v4rhosigmalapltau, v4rhosigmatau2, v4rholapl3, v4rholapl2tau, v4rholapltau2, v4rhotau3, v4sigma4, v4sigma3lapl, v4sigma3tau, v4sigma2lapl2, v4sigma2lapltau, v4sigma2tau2, v4sigmalapl3, v4sigmalapl2tau, v4sigmalapltau2, v4sigmatau3, v4lapl4, v4lapl3tau, v4lapl2tau2, v4lapltau3, v4tau4)
end

function xc_lda_exc(p, np, rho, zk)
    ccall((:xc_lda_exc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, zk)
end

function xc_gga_exc(p, np, rho, sigma, zk)
    ccall((:xc_gga_exc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, zk)
end

function xc_mgga_exc(p, np, rho, sigma, lapl, tau, zk)
    ccall((:xc_mgga_exc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, lapl, tau, zk)
end

function xc_lda_exc_vxc(p, np, rho, zk, vrho)
    ccall((:xc_lda_exc_vxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, zk, vrho)
end

function xc_gga_exc_vxc(p, np, rho, sigma, zk, vrho, vsigma)
    ccall((:xc_gga_exc_vxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, zk, vrho, vsigma)
end

function xc_mgga_exc_vxc(p, np, rho, sigma, lapl, tau, zk, vrho, vsigma, vlapl, vtau)
    ccall((:xc_mgga_exc_vxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, lapl, tau, zk, vrho, vsigma, vlapl, vtau)
end

function xc_lda_vxc(p, np, rho, vrho)
    ccall((:xc_lda_vxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, vrho)
end

function xc_gga_vxc(p, np, rho, sigma, vrho, vsigma)
    ccall((:xc_gga_vxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, vrho, vsigma)
end

function xc_mgga_vxc(p, np, rho, sigma, lapl, tau, vrho, vsigma, vlapl, vtau)
    ccall((:xc_mgga_vxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, lapl, tau, vrho, vsigma, vlapl, vtau)
end

function xc_lda_exc_vxc_fxc(p, np, rho, zk, vrho, v2rho2)
    ccall((:xc_lda_exc_vxc_fxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, zk, vrho, v2rho2)
end

function xc_gga_exc_vxc_fxc(p, np, rho, sigma, zk, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2)
    ccall((:xc_gga_exc_vxc_fxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, zk, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2)
end

function xc_mgga_exc_vxc_fxc(p, np, rho, sigma, lapl, tau, zk, vrho, vsigma, vlapl, vtau, v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, v2lapl2, v2lapltau, v2tau2)
    ccall((:xc_mgga_exc_vxc_fxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, lapl, tau, zk, vrho, vsigma, vlapl, vtau, v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, v2lapl2, v2lapltau, v2tau2)
end

function xc_lda_vxc_fxc(p, np, rho, vrho, v2rho2)
    ccall((:xc_lda_vxc_fxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, vrho, v2rho2)
end

function xc_gga_vxc_fxc(p, np, rho, sigma, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2)
    ccall((:xc_gga_vxc_fxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2)
end

function xc_mgga_vxc_fxc(p, np, rho, sigma, lapl, tau, vrho, vsigma, vlapl, vtau, v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, v2lapl2, v2lapltau, v2tau2)
    ccall((:xc_mgga_vxc_fxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, lapl, tau, vrho, vsigma, vlapl, vtau, v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, v2lapl2, v2lapltau, v2tau2)
end

function xc_lda_fxc(p, np, rho, v2rho2)
    ccall((:xc_lda_fxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, v2rho2)
end

function xc_gga_fxc(p, np, rho, sigma, v2rho2, v2rhosigma, v2sigma2)
    ccall((:xc_gga_fxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, v2rho2, v2rhosigma, v2sigma2)
end

function xc_mgga_fxc(p, np, rho, sigma, lapl, tau, v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, v2lapl2, v2lapltau, v2tau2)
    ccall((:xc_mgga_fxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, lapl, tau, v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, v2lapl2, v2lapltau, v2tau2)
end

function xc_lda_exc_vxc_fxc_kxc(p, np, rho, zk, vrho, v2rho2, v3rho3)
    ccall((:xc_lda_exc_vxc_fxc_kxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, zk, vrho, v2rho2, v3rho3)
end

function xc_gga_exc_vxc_fxc_kxc(p, np, rho, sigma, zk, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2, v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3)
    ccall((:xc_gga_exc_vxc_fxc_kxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, zk, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2, v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3)
end

function xc_mgga_exc_vxc_fxc_kxc(p, np, rho, sigma, lapl, tau, zk, vrho, vsigma, vlapl, vtau, v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, v2lapl2, v2lapltau, v2tau2, v3rho3, v3rho2sigma, v3rho2lapl, v3rho2tau, v3rhosigma2, v3rhosigmalapl, v3rhosigmatau, v3rholapl2, v3rholapltau, v3rhotau2, v3sigma3, v3sigma2lapl, v3sigma2tau, v3sigmalapl2, v3sigmalapltau, v3sigmatau2, v3lapl3, v3lapl2tau, v3lapltau2, v3tau3)
    ccall((:xc_mgga_exc_vxc_fxc_kxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, lapl, tau, zk, vrho, vsigma, vlapl, vtau, v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, v2lapl2, v2lapltau, v2tau2, v3rho3, v3rho2sigma, v3rho2lapl, v3rho2tau, v3rhosigma2, v3rhosigmalapl, v3rhosigmatau, v3rholapl2, v3rholapltau, v3rhotau2, v3sigma3, v3sigma2lapl, v3sigma2tau, v3sigmalapl2, v3sigmalapltau, v3sigmatau2, v3lapl3, v3lapl2tau, v3lapltau2, v3tau3)
end

function xc_lda_vxc_fxc_kxc(p, np, rho, vrho, v2rho2, v3rho3)
    ccall((:xc_lda_vxc_fxc_kxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, vrho, v2rho2, v3rho3)
end

function xc_gga_vxc_fxc_kxc(p, np, rho, sigma, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2, v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3)
    ccall((:xc_gga_vxc_fxc_kxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2, v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3)
end

function xc_mgga_vxc_fxc_kxc(p, np, rho, sigma, lapl, tau, vrho, vsigma, vlapl, vtau, v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, v2lapl2, v2lapltau, v2tau2, v3rho3, v3rho2sigma, v3rho2lapl, v3rho2tau, v3rhosigma2, v3rhosigmalapl, v3rhosigmatau, v3rholapl2, v3rholapltau, v3rhotau2, v3sigma3, v3sigma2lapl, v3sigma2tau, v3sigmalapl2, v3sigmalapltau, v3sigmatau2, v3lapl3, v3lapl2tau, v3lapltau2, v3tau3)
    ccall((:xc_mgga_vxc_fxc_kxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, lapl, tau, vrho, vsigma, vlapl, vtau, v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, v2lapl2, v2lapltau, v2tau2, v3rho3, v3rho2sigma, v3rho2lapl, v3rho2tau, v3rhosigma2, v3rhosigmalapl, v3rhosigmatau, v3rholapl2, v3rholapltau, v3rhotau2, v3sigma3, v3sigma2lapl, v3sigma2tau, v3sigmalapl2, v3sigmalapltau, v3sigmatau2, v3lapl3, v3lapl2tau, v3lapltau2, v3tau3)
end

function xc_lda_kxc(p, np, rho, v3rho3)
    ccall((:xc_lda_kxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, v3rho3)
end

function xc_gga_kxc(p, np, rho, sigma, v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3)
    ccall((:xc_gga_kxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3)
end

function xc_mgga_kxc(p, np, rho, sigma, lapl, tau, v3rho3, v3rho2sigma, v3rho2lapl, v3rho2tau, v3rhosigma2, v3rhosigmalapl, v3rhosigmatau, v3rholapl2, v3rholapltau, v3rhotau2, v3sigma3, v3sigma2lapl, v3sigma2tau, v3sigmalapl2, v3sigmalapltau, v3sigmatau2, v3lapl3, v3lapl2tau, v3lapltau2, v3tau3)
    ccall((:xc_mgga_kxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, lapl, tau, v3rho3, v3rho2sigma, v3rho2lapl, v3rho2tau, v3rhosigma2, v3rhosigmalapl, v3rhosigmatau, v3rholapl2, v3rholapltau, v3rhotau2, v3sigma3, v3sigma2lapl, v3sigma2tau, v3sigmalapl2, v3sigmalapltau, v3sigmatau2, v3lapl3, v3lapl2tau, v3lapltau2, v3tau3)
end

function xc_lda_lxc(p, np, rho, v4rho4)
    ccall((:xc_lda_lxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, v4rho4)
end

function xc_gga_lxc(p, np, rho, sigma, v4rho4, v4rho3sigma, v4rho2sigma2, v4rhosigma3, v4sigma4)
    ccall((:xc_gga_lxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, v4rho4, v4rho3sigma, v4rho2sigma2, v4rhosigma3, v4sigma4)
end

function xc_mgga_lxc(p, np, rho, sigma, lapl, tau, v4rho4, v4rho3sigma, v4rho3lapl, v4rho3tau, v4rho2sigma2, v4rho2sigmalapl, v4rho2sigmatau, v4rho2lapl2, v4rho2lapltau, v4rho2tau2, v4rhosigma3, v4rhosigma2lapl, v4rhosigma2tau, v4rhosigmalapl2, v4rhosigmalapltau, v4rhosigmatau2, v4rholapl3, v4rholapl2tau, v4rholapltau2, v4rhotau3, v4sigma4, v4sigma3lapl, v4sigma3tau, v4sigma2lapl2, v4sigma2lapltau, v4sigma2tau2, v4sigmalapl3, v4sigmalapl2tau, v4sigmalapltau2, v4sigmatau3, v4lapl4, v4lapl3tau, v4lapl2tau2, v4lapltau3, v4tau4)
    ccall((:xc_mgga_lxc, libxc), Cvoid, (Ptr{xc_func_type}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, np, rho, sigma, lapl, tau, v4rho4, v4rho3sigma, v4rho3lapl, v4rho3tau, v4rho2sigma2, v4rho2sigmalapl, v4rho2sigmatau, v4rho2lapl2, v4rho2lapltau, v4rho2tau2, v4rhosigma3, v4rhosigma2lapl, v4rhosigma2tau, v4rhosigmalapl2, v4rhosigmalapltau, v4rhosigmatau2, v4rholapl3, v4rholapl2tau, v4rholapltau2, v4rhotau3, v4sigma4, v4sigma3lapl, v4sigma3tau, v4sigma2lapl2, v4sigma2lapltau, v4sigma2tau2, v4sigmalapl3, v4sigmalapl2tau, v4sigmalapltau2, v4sigmatau3, v4lapl4, v4lapl3tau, v4lapl2tau2, v4lapltau3, v4tau4)
end

function xc_gga_ak13_get_asymptotic(homo)
    ccall((:xc_gga_ak13_get_asymptotic, libxc), Cdouble, (Cdouble,), homo)
end

function xc_gga_ak13_pars_get_asymptotic(homo, ext_params)
    ccall((:xc_gga_ak13_pars_get_asymptotic, libxc), Cdouble, (Cdouble, Ptr{Cdouble}), homo, ext_params)
end

function xc_hyb_exx_coef(p)
    ccall((:xc_hyb_exx_coef, libxc), Cdouble, (Ptr{xc_func_type},), p)
end

function xc_hyb_cam_coef(p, omega, alpha, beta)
    ccall((:xc_hyb_cam_coef, libxc), Cvoid, (Ptr{xc_func_type}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), p, omega, alpha, beta)
end

function xc_nlc_coef(p, nlc_b, nlc_C)
    ccall((:xc_nlc_coef, libxc), Cvoid, (Ptr{xc_func_type}, Ptr{Cdouble}, Ptr{Cdouble}), p, nlc_b, nlc_C)
end


text:
type Texture
    ptr::Ptr{Void}

    function Texture(ptr::Ptr{Void})
        t = new(ptr)
        # finalizer(t, destroy)
        # t
    end
end

Texture() = Texture(C_NULL)
function Texture(width::Integer, height::Integer)
    Texture(ccall((:sfTexture_create, libcsfml_graphics), Ptr{Void}, (UInt32, UInt32,), width, height))
end

function Texture(filename::AbstractString)
    Texture(ccall((:sfTexture_createFromFile, libcsfml_graphics), Ptr{Void}, (Ptr{Cchar}, Ptr{Void},), filename, C_NULL))
end

function Texture(image::Image)
    Texture(ccall((:sfTexture_createFromImage, libcsfml_graphics), Ptr{Void}, (Ptr{Void}, Ptr{Void},), image.ptr, C_NULL))
end

function copy(texture::Texture)
    return Texture(ccall((:sfTexture_copy, libcsfml_graphics), Ptr{Void}, (Ptr{Void},), texture.ptr))
end

function destroy(texture::Texture)
    # println("Destroy $(texture.ptr)")
    ccall((:sfTexture_destroy, libcsfml_graphics), Void, (Ptr{Void},), texture.ptr)
end

function get_size(texture::Texture)
    return ccall((:sfTexture_getSize, libcsfml_graphics), Vector2u, (Ptr{Void},), texture.ptr)
end

function copy_to_image(texture::Texture)
    return Image(ccall((:sfTexture_copyToImage, libcsfml_graphics), Ptr{Void}, (Ptr{Void},), texture.ptr))
end

function update(texture::Texture, window::RenderWindow, width::Integer=0, height::Integer=0)
    ccall((:sfTexture_updateFromWindow, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void}, UInt32, UInt32,), texture.ptr, window.ptr, width, height)
end

function update(texture::Texture, image::Image, x_off::Integer=0, y_off::Integer=0)
    ccall((:sfTexture_updateFromImage, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void}, UInt32, UInt32,), texture.ptr, image.ptr, x_off, y_off)
end

function set_smooth(texture::Texture, smooth::Bool)
    ccall((:sfTexture_setSmooth, libcsfml_graphics), Void, (Ptr{Void}, Int32,), texture.ptr, smooth)
end

function is_smooth(texture::Texture)
    ccall((:sfTexture_isSmooth, libcsfml_graphics), Bool, (Ptr{Void},), texture.ptr)
end

function set_repeated(texture::Texture, repeated::Bool)
    ccall((:sfTexture_setRepeated, libcsfml_graphics), Void, (Ptr{Void}, Int32), texture.ptr, repeated)
end

function is_repeated(texture::Texture)
    ccall((:sfTexture_isRepeated, libcsfml_graphics), Bool, (Ptr{Void},), texture.ptr)
end


text:
##############################  Gradients  ###################################
function grad_potential_energy(r::Atom, p::SNAP)
    println("Not yet implmented.")
    return 0.0
end

function grad_force(r::Atom, p::SNAP)
    println("Not yet implemented.")
    return 0.0
end

function grad_virial(r::Atom, p::SNAP)
    println("Not yet implemented")
    return 0.0
end


text:
"""
    Revise.WatchList

A struct for holding files that live inside a directory.
Some platforms (OSX) have trouble watching too many files. So we
watch parent directories, and keep track of which files in them
should be tracked.

Fields:
- `timestamp`: mtime of last update
- `trackedfiles`: Set of filenames, generally expressed as a relative path
"""
mutable struct WatchList
    timestamp::Float64         # unix time of last revision
    trackedfiles::Dict{String,PkgId}
end

const DocExprs = Dict{Module,Vector{Expr}}
const ExprsSigs = OrderedDict{RelocatableExpr,Union{Nothing,Vector{Any}}}
const DepDictVals = Tuple{Module,RelocatableExpr}
const DepDict = Dict{Symbol,Set{DepDictVals}}

Base.setindex!(ex_sigs::ExprsSigs, val, ex::Expr) = setindex!(ex_sigs, val, RelocatableExpr(ex))

function Base.show(io::IO, exsigs::ExprsSigs)
    compact = get(io, :compact, false)
    if compact
        n = 0
        for (rex, sigs) in exsigs
            sigs === nothing && continue
            n += length(sigs)
        end
        print(io, "ExprsSigs(<$(length(exsigs)) expressions>, <$n signatures>)")
    else
        print(io, "ExprsSigs with the following expressions: ")
        for def in keys(exsigs)
            print(io, "\n  ")
            Base.show_unquoted(io, RelocatableExpr(unwrap(def)), 2)
        end
    end
end

"""
    ModuleExprsSigs

For a particular source file, the corresponding `ModuleExprsSigs` is a mapping
`mod=>exprs=>sigs` of the expressions `exprs` found in `mod` and the signatures `sigs`
that arise from them. Specifically, if `mes` is a `ModuleExprsSigs`, then `mes[mod][ex]`
is a list of signatures that result from evaluating `ex` in `mod`. It is possible that
this returns `nothing`, which can mean either that `ex` does not define any methods
or that the signatures have not yet been cached.

The first `mod` key is guaranteed to be the module into which this file was `include`d.

To create a `ModuleExprsSigs` from a source file, see [`Revise.parse_source`](@ref).
"""
const ModuleExprsSigs = OrderedDict{Module,ExprsSigs}

if VERSION >= v"1.5.0-DEV.472"
    Base.typeinfo_prefix(io::IO, mexs::ModuleExprsSigs) = string(typeof(mexs).name), true
else
    Base.typeinfo_prefix(io::IO, mexs::ModuleExprsSigs) = string(typeof(mexs).name)
end

"""
    fm = ModuleExprsSigs(mod::Module)

Initialize an empty `ModuleExprsSigs` for a file that is `include`d into `mod`.
"""
ModuleExprsSigs(mod::Module) = ModuleExprsSigs(mod=>ExprsSigs())

Base.isempty(fm::ModuleExprsSigs) = length(fm) == 1 && isempty(first(values(fm)))

"""
    FileInfo(mexs::ModuleExprsSigs, cachefile="")

Structure to hold the per-module expressions found when parsing a
single file.
`mexs` holds the [`Revise.ModuleExprsSigs`](@ref) for the file.

Optionally, a `FileInfo` can also record the path to a cache file holding the original source code.
This is applicable only for precompiled modules and `Base`.
(This cache file is distinct from the original source file that might be edited by the
developer, and it will always hold the state
of the code when the package was precompiled or Julia's `Base` was built.)
When a cache is available, `mexs` will be empty until the file gets edited:
the original source code gets parsed only when a revision needs to be made.

Source cache files greatly reduce the overhead of using Revise.
"""
struct FileInfo
    modexsigs::ModuleExprsSigs
    cachefile::String
end
FileInfo(fm::ModuleExprsSigs) = FileInfo(fm, "")

"""
    FileInfo(mod::Module, cachefile="")

Initialze an empty FileInfo for a file that is `include`d into `mod`.
"""
FileInfo(mod::Module, cachefile::AbstractString="") = FileInfo(ModuleExprsSigs(mod), cachefile)

FileInfo(fm::ModuleExprsSigs, fi::FileInfo) = FileInfo(fm, fi.cachefile)

function Base.show(io::IO, fi::FileInfo)
    print(io, "FileInfo(")
    for (mod, exsigs) in fi.modexsigs
        show(io, mod)
        print(io, "=>")
        show(io, exsigs)
        print(io, ", ")
    end
    if !isempty(fi.cachefile)
        print(io, "with cachefile ", fi.cachefile)
    end
    print(io, ')')
end

"""
    PkgData(id, path, fileinfos::Dict{String,FileInfo})

A structure holding the data required to handle a particular package.
`path` is the top-level directory defining the package,
and `fileinfos` holds the [`Revise.FileInfo`](@ref) for each file defining the package.

For the `PkgData` associated with `Main` (e.g., for files loaded with [`includet`](@ref)),
the corresponding `path` entry will be empty.
"""
mutable struct PkgData
    info::PkgFiles
    fileinfos::Vector{FileInfo}
end

PkgData(id::PkgId, path) = PkgData(PkgFiles(id, path), FileInfo[])
PkgData(id::PkgId, ::Nothing) = PkgData(id, "")
function PkgData(id::PkgId)
    bp = basepath(id)
    if !isempty(bp)
        bp = normpath(bp)
    end
    PkgData(id, bp)
end

# Abstraction interface for PkgData
Base.PkgId(pkgdata::PkgData) = PkgId(pkgdata.info)
CodeTracking.basedir(pkgdata::PkgData) = basedir(pkgdata.info)
CodeTracking.srcfiles(pkgdata::PkgData) = srcfiles(pkgdata.info)

function fileindex(info, file::AbstractString)
    for (i, f) in enumerate(srcfiles(info))
        f == file && return i
    end
    return nothing
end

function hasfile(info, file)
    if isabspath(file)
        file = relpath(file, info)
    end
    fileindex(info, file) !== nothing
end

function fileinfo(pkgdata::PkgData, file::AbstractString)
    i = fileindex(pkgdata, file)
    i === nothing && error("file ", file, " not found")
    return pkgdata.fileinfos[i]
end
fileinfo(pkgdata::PkgData, i::Integer) = pkgdata.fileinfos[i]

function Base.push!(pkgdata::PkgData, pr::Pair{<:AbstractString,FileInfo})
    push!(srcfiles(pkgdata), pr.first)
    push!(pkgdata.fileinfos, pr.second)
    return pkgdata
end

function Base.show(io::IO, pkgdata::PkgData)
    compact = get(io, :compact, false)
    print(io, "PkgData(")
    if compact
        print(io, '"', pkgdata.info.basedir, "\", ")
        nexs, nsigs, nparsed = 0, 0, 0
        for fi in pkgdata.fileinfos
            thisnexs, thisnsigs = 0, 0
            for (mod, exsigs) in fi.modexsigs
                for (rex, sigs) in exsigs
                    thisnexs += 1
                    sigs === nothing && continue
                    thisnsigs += length(sigs)
                end
            end
            nexs += thisnexs
            nsigs += thisnsigs
            if thisnexs > 0
                nparsed += 1
            end
        end
        print(io, nparsed, '/', length(pkgdata.fileinfos), " parsed files, ", nexs, " expressions, ", nsigs, " signatures)")
    else
        show(io, pkgdata.info.id)
        println(io, ", basedir \"", pkgdata.info.basedir, ':')
        for (f, fi) in zip(pkgdata.info.files, pkgdata.fileinfos)
            print(io, "  \"", f, "\": ")
            show(IOContext(io, :compact=>true), fi)
            print(io, '\n')
        end
    end
end

struct GitRepoException <: Exception
    filename::String
end

function Base.showerror(io::IO, ex::GitRepoException)
    print(io, "no repository at ", ex.filename, " to track stdlibs you must build Julia from source")
end

"""
    thunk = TaskThunk(f, args)

To facilitate precompilation and reduce latency, we avoid creation of anonymous thunks.
`thunk` can be used as an argument in `schedule(Task(thunk))`.
"""
struct TaskThunk
    f          # deliberately untyped
    args       # deliberately untyped
end

@noinline (thunk::TaskThunk)() = thunk.f(thunk.args...)


text:
#
# This file is part of the Simulate.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#

"""
    𝐶

`𝐶` (𝐶 = `\\itC`+`tab`) is the default simulation clock. If you do one
simulation at a time, you can use it for time keeping.

# Examples

```jldoctest
julia> using Simulate

julia> reset!(𝐶)
"clock reset to t₀=0.0, sampling rate Δt=0.0."

julia> 𝐶  # default clock
Clock thread 1 (+ 0 ac): state=Simulate.Idle(), t=0.0 , Δt=0.0 , prc:0
  scheduled ev:0, cev:0, sampl:0
```
"""
const 𝐶 = Clock()

# ----------------------------------------------------
# setting and getting clock parameters
# ----------------------------------------------------

"""
    setUnit!(clk::Clock, new::FreeUnits)

set a clock to a new time unit in `Unitful`. If necessary convert
current clock times to the new unit.

# Arguments
- `clk::Clock`
- `new::FreeUnits`: new is one of `ms`, `s`, `minute` or `hr` or another Unitful
    `Time` unit.

# Examples
```jldoctest
julia> using Simulate, Unitful

julia> import Unitful: Time, s, minute, hr

julia> c = Clock(t0=60)     # setup a new clock with t0=60
Clock thread 1 (+ 0 ac): state=Simulate.Undefined(), t=60.0 , Δt=0.0 , prc:0
  scheduled ev:0, cev:0, sampl:0

julia> tau(c) # current time is 60.0 NoUnits
60.0

julia> setUnit!(c, s)       # set clock unit to Unitful.s
60.0 s

julia> tau(c) # current time is now 60.0 s
60.0 s

julia> setUnit!(c, minute)  # set clock unit to Unitful.minute
1.0 minute

julia> tau(c)               # current time is now 1.0 minute
1.0 minute

julia> isa(tau(c), Time)
true

julia> uconvert(s, tau(c))  # ... which can be converted to other time units
60.0 s

julia> tau(c).val           # it has a value of 1.0
1.0

julia> c.time               # internal clock time is set to 1.0 (a Float64)
1.0

julia> c.unit               # internal clock unit is set to Unitful.minute
minute
```
"""
function setUnit!(clk::Clock, new::FreeUnits)
    if isa(1new, Time)
        if clk.unit == new
            println("clock is already set to $new")
        elseif clk.unit == NoUnits
            clk.unit = new
        else
            old = clk.unit
            clk.unit = new
            fac = uconvert(new, 1*old).val
            clk.time *= fac
            clk.end_time *= fac
            clk.tev *= fac
            clk.Δt *= fac
            clk.tn *= fac
        end
    else
        clk.unit = NoUnits
    end
    tau(clk)
end

"""
```
tau(clk::Clock=𝐶)
```
Return the current simulation time.

# Examples

```jldoctest
julia> using Simulate

julia> reset!(𝐶)
"clock reset to t₀=0.0, sampling rate Δt=0.0."
julia> tau() # gives the central time
0.0
```
"""
tau(clk::Clock=𝐶) = clk.unit == NoUnits ? clk.time : clk.time*clk.unit

"""
```
sync!(clk::Clock, to::Clock=𝐶)
```
Force a synchronization of two clocks. Change all registered times of
`clk` accordingly. Convert or force clk.unit to to.unit.
"""
function sync!(clk::Clock, to::Clock=𝐶)
    if (clk.unit == NoUnits) | (clk.unit == to.unit)
        fac = 1
    elseif to.unit == NoUnits
        println(stderr, "Warning: deleted time unit without conversion")
        fac = 1
    else
        fac = uconvert(to.unit, 1clk.unit).val
    end
    Δt = to.time - clk.time*fac
    clk.time = clk.time*fac + Δt
    clk.unit = to.unit
    clk.tn  = clk.tn*fac + Δt
    clk.tev  = clk.tev*fac + Δt
    clk.end_time = clk.end_time*fac + Δt
    clk.Δt = to.Δt
    evq = PriorityQueue{DiscreteEvent,Float64}()
    for (ev, t) ∈ pairs(clk.sc.events)
        evq[ev] = t*fac + Δt
    end
    clk.sc.events = evq
    clk
end

"""
```
reset!(clk::Clock, Δt::Number=0; t0::Number=0, hard::Bool=true, unit=NoUnits)
```
reset a clock

# Arguments
- `clk::Clock`
- `Δt::Number=0`: time increment
- `t0::Float64=0` or `t0::Time`: start time
- `hard::Bool=true`: time is reset, all scheduled events and sampling are
    deleted. If hard=false, then only time is reset, event and
    sampling times are adjusted accordingly.
- `unit=NoUnits`: the Time unit for the clock after reset. If a `Δt::Time` is
    given, its Time unit goes into the clock Time unit. If only t0::Time is given,
    its Time unit goes into the clock time unit.

# Examples

```jldoctest
julia> using Simulate, Unitful

julia> import Unitful: s

julia> c = Clock(1s, t0=60s)
Clock thread 1 (+ 0 ac): state=Simulate.Undefined(), t=60.0 s, Δt=1.0 s, prc:0
  scheduled ev:0, cev:0, sampl:0

julia> reset!(c)
"clock reset to t₀=0.0, sampling rate Δt=0.0."

julia> c
Clock thread 1 (+ 0 ac): state=Simulate.Idle(), t=0.0 , Δt=0.0 , prc:0
  scheduled ev:0, cev:0, sampl:0
```
"""
function reset!(clk::Clock, Δt::Number=0;
                t0::Number=0, hard::Bool=true, unit=NoUnits)
    if  isa(1unit, Time)
        Δt = isa(Δt, Time) ? uconvert(unit, Δt).val : Δt
        t0 = isa(t0, Time) ? uconvert(unit, t0).val : t0
    elseif isa(Δt, Time)
        unit = Unitful.unit(1Δt)
        Δt = Δt.val
        t0 = isa(t0, Time) ? uconvert(unit, t0).val : t0
    elseif isa(t0, Time)
        unit = Unitful.unit(t0)
        t0 = t0.val
    else
        nothing
    end
    if hard
        clk.state = Idle()
        clk.time = t0
        clk.unit = unit
        clk.tn = t0
        clk.tev = t0
        clk.end_time = t0
        clk.evcount = 0
        clk.scount = 0
        clk.Δt = Δt
        clk.sc.events = PriorityQueue{DiscreteEvent,Float64}()
        clk.sc.cevents = DiscreteCond[]
        clk.processes = Dict{Any, Prc}()
        clk.sc.samples = Sample[]
    else
        sync!(clk, Clock(Δt, t0=t0, unit=unit))
    end
    "clock reset to t₀=$(float(t0*unit)), sampling rate Δt=$(float(Δt*unit))."
end

"""
    tadjust(clk::Clock, t::Unitful.Time) :: Float64

Adjust/convert `t` given according to clock settings and return a Float64 value.
"""
function tadjust(clk::Clock, t::Unitful.Time) :: Float64
    if clk.unit == NoUnits
        println(stderr, "Warning: clock has no time unit, ignoring units")
        return t.val
    else
        return uconvert(clk.unit, t).val
    end
end

"""
```
sample_time!([clk::Clock], Δt::Number)
```
set the clock's sample rate starting from now (`tau(clk)`).

# Arguments
- `clk::Clock`: if not supplied, set the sample rate on 𝐶,
- `Δt::Number`: sample rate, time interval for sampling
"""
function sample_time!(clk::Clock, Δt::Number)
    clk.Δt = Δt isa Unitful.Time ? tadjust(clk, Δt) : Δt
    clk.tn = clk.time + clk.Δt
end
sample_time!(Δt::Number) = sample_time!(𝐶, Δt)

"""
```
periodic!([clk::Clock], ex::Union{Expr, Fun}, Δt::Number=clk.Δt;
        scope::Module=Main, spawn=false)
```
Register a function or expression for periodic execution at the clock`s sample rate.

# Arguments
- `clk::Clock`: if not supplied, it samples on 𝐶,
- `ex::Union{Expr, Fun}`: an expression or function,
- `Δt::Number=clk.Δt`: set the clock's sampling rate, if no Δt is given, it takes
    the current sampling rate, if that is 0, it calculates one,
- `scope::Module=Main`: optional, an evaluation scope for a given expression.
"""
function periodic!(clk::Clock, ex::Union{Expr, Fun}, Δt::Number=clk.Δt;
                 scope::Module=Main, spawn=false)
    clk.Δt = Δt == 0 ? scale(clk.end_time - clk.time)/100 : Δt
    assign(clk, Sample(ex, scope), spawn ? spawnid(clk) : 0)
end
periodic!(ex::Union{Expr, Fun}, Δt::Number=𝐶.Δt; kw...) = periodic!(𝐶, ex, Δt; kw...)

"Is a Clock busy?"
busy(clk::Clock) = clk.state == Busy()

# ----------------------------------------------------
# step! transition functions for Clocks
# ----------------------------------------------------

"""
    step!(clk::Clock, ::Undefined, ::Init)

initialize a clock.
"""
function step!(clk::Clock, ::Undefined, ::Init)
    clk.state = Idle()
end

"""
    step!(clk::Clock, ::Undefined, σ::Union{Step,Run})

if uninitialized, initialize and then Step or Run.
"""
function step!(clk::Clock, ::Undefined, σ::Union{Step,Run})
    step!(clk, clk.state, Init(0))
    step!(clk, clk.state, σ)
end

"""
    setTimes(clk::Clock)

set clock times for next event or sampling action. The internal clock times
`clk.tev` and `clk.tn` must always be set to be at least `clk.time`.
"""
function setTimes(clk::Clock)
    if !isempty(clk.sc.events)
        clk.tev = nextevtime(clk)
        clk.tn = clk.Δt > 0 ? clk.time + clk.Δt : clk.time
    else
        clk.tn = clk.Δt > 0 ? clk.time + clk.Δt : clk.time
        clk.tev = clk.tn
    end
end

"""
    step!(c::Clock, ::Union{Idle,Busy,Halted}, ::Step)

step forward to next tick or scheduled event.

At a tick evaluate 1) all sampling functions or expressions, 2) all conditional
events, then 3) if an event is encountered, trigger the event.

The internal clock times `c.tev` and `c.tn` are always at least `c.time`.
"""
step!(c::Clock, ::Union{Idle,Halted}, ::Step) = do_step!(c)

"""
    do_run(c::Clock, Δt::Float64)

Run a clock for a time Δt.
"""
function do_run!(c::Clock, Δt::Float64)
    c.end_time = c.time + Δt
    c.evcount = 0
    c.scount = 0
    setTimes(c)
    while any(i->(c.time < i ≤ c.end_time), (c.tn, c.tev))
        do_step!(c)
        if c.state == Halted()
            return c.end_time
        end
    end
    c.time = c.end_time
end

"""
    step!(clk::Clock, ::Idle, σ::Run)

Run a simulation for a given duration.

The duration is given with `Run(duration)`. Call scheduled events and evaluate
sampling expressions at each tick in that timeframe.
"""
function step!(clk::Clock, ::Idle, σ::Run)
    tend = do_run!(clk, σ.duration)
    if clk.state == Halted()
        return
    end
    # catch remaining events scheduled for the end_time
    while (length(clk.sc.events) ≥ 1) && (nextevtime(clk) ≤ tend + Base.eps(tend)*10)
        step!(clk, clk.state, Step())
        tend = nextfloat(tend)
    end
    clk.time = clk.end_time

    "run! finished with $(clk.evcount) clock events, $(clk.scount) sample steps, simulation time: $(clk.time)"
end

"""
    step!(clk::Clock, ::Busy, ::Stop)

Stop the clock.
"""
function step!(clk::Clock, ::Busy, ::Stop)
    clk.state = Halted()
    "Halted after $(clk.evcount) events, simulation time: $(clk.time)"
end

"""
    step!(clk::Clock, ::Halted, ::Resume)

Resume a halted clock.
"""
function step!(clk::Clock, ::Halted, ::Resume)
    clk.state = Idle()
    step!(clk, clk.state, Run(clk.end_time - clk.time))
end

"""
    step!(clk::Clock, q::ClockState, σ::ClockEvent)

catch all step!-function.
"""
function step!(clk::Clock, q::ClockState, σ::ClockEvent)
    println(stderr, "Warning: undefined transition ",
            "$(typeof(clk)), ::$(typeof(q)), ::$(typeof(σ)))\n",
            "maybe, you should reset! the clock!")
end

# --------------------------------------------------------
# Clock user interface functions
# --------------------------------------------------------

"""
    run!(clk::Clock, duration::Number)

Run a simulation for a given duration. Call scheduled events and evaluate
sampling expressions at each tick in that timeframe.
"""
function run!(clk::Clock, duration::Number)
    duration = duration isa Unitful.Time ? tadjust(clk, duration) : duration
    step!(clk, clk.state, Run(duration))
end

"""
    incr!(clk::Clock)

Take one simulation step, execute the next tick or event.
"""
incr!(clk::Clock) = step!(clk, clk.state, Step())

"""
    stop!(clk::Clock)

Stop a running simulation.
"""
stop!(clk::Clock) = step!(clk, clk.state, Stop())

"""
    resume!(clk::Clock)

Resume a halted simulation.
"""
resume!(clk::Clock) = step!(clk, clk.state, Resume())

"""
    init!(clk::Clock)

initialize a clock.
"""
init!(clk::Clock) = step!(clk, clk.state, Init(""))


text:
using DifferentialEquations
using NLsolve

function f(du, u, p,t )
    du[1] = 2- 2u[1]
    du[2] = u[1] - 4u[2]
end

u0 = zeros(2)

# Build problem
prob = SteadyStateProblem(f,u0)
sol = solve(prob,SSRootfind())

sol = solve(prob,SSRootfind(nlsolve = (f,u0,abstol) -> (res=NLsolve.nlsolve(f,u0,autodiff=:forward,method=:newton,iterations=Int(1e6),ftol=abstol);res.zero) ))


text:
function DeterministicMetadata(ts::AbstractDeterministic)
    return DeterministicMetadata(
        get_name(ts),
        get_resolution(ts),
        get_initial_timestamp(ts),
        get_interval(ts),
        get_count(ts),
        get_uuid(ts),
        get_horizon(ts),
        typeof(ts),
        get_scaling_factor_multiplier(ts),
    )
end

function serialize(::Type{<:T}) where {T <: AbstractDeterministic}
    # This currently cannot be done for all InfrastructureSystemsTypes.
    # Some are encoded directly as strings.
    @debug "serialize" _group = LOG_GROUP_SERIALIZATION T
    data = Dict{String, Any}()
    add_serialization_metadata!(data, T)
    return data
end

function deserialize_to_dict(::Type{T}, data::Dict) where {T <: DeterministicMetadata}
    # This is custom because of time_series_type.
    # Duplicated from src/serialization.jl
    vals = Dict{Symbol, Any}()
    for (field_name, field_type) in zip(fieldnames(T), fieldtypes(T))
        val = data[string(field_name)]
        if val isa Dict && haskey(val, METADATA_KEY)
            metadata = get_serialization_metadata(val)
            if haskey(metadata, FUNCTION_KEY)
                vals[field_name] = deserialize(Function, val)
            else
                type = get_type_from_serialization_metadata(metadata)
                if field_name == :time_series_type
                    vals[field_name] = type
                else
                    vals[field_name] = deserialize(type, val)
                end
            end
        else
            vals[field_name] = deserialize(field_type, val)
        end
    end
    return vals
end


text:
"""
    HubbardReal1D(address; u=1.0, t=1.0)

Implements a one-dimensional Bose Hubbard chain in real space.

```math
\\hat{H} = -t \\sum_{\\langle i,j\\rangle} a_i^† a_j + \\frac{u}{2}\\sum_i n_i (n_i-1)
```

# Arguments

* `address`: the starting address, defines number of particles and sites.
* `u`: the interaction parameter.
* `t`: the hopping strength.

# See also

* [`HubbardMom1D`](@ref)
* [`ExtendedHubbardReal1D`](@ref)

"""
struct HubbardReal1D{TT,A<:AbstractFockAddress,U,T} <: AbstractHamiltonian{TT}
    add::A
end

function HubbardReal1D(addr; u=1.0, t=1.0)
    U, T = promote(float(u), float(t))
    return HubbardReal1D{typeof(U),typeof(addr),U,T}(addr)
end

function Base.show(io::IO, h::HubbardReal1D)
    print(io, "HubbardReal1D($(h.add); u=$(h.u), t=$(h.t))")
end

function starting_address(h::HubbardReal1D)
    return getfield(h, :add)
end

LOStructure(::Type{<:HubbardReal1D{<:Real}}) = IsHermitian()

Base.getproperty(h::HubbardReal1D, s::Symbol) = getproperty(h, Val(s))
Base.getproperty(h::HubbardReal1D{<:Any,<:Any,U}, ::Val{:u}) where U = U
Base.getproperty(h::HubbardReal1D{<:Any,<:Any,<:Any,T}, ::Val{:t}) where T = T
Base.getproperty(h::HubbardReal1D, ::Val{:add}) = getfield(h, :add)

function num_offdiagonals(::HubbardReal1D, address::BoseFS)
    return 2 * num_occupied_modes(address)
end

"""
    bose_hubbard_interaction(address)

Return Σ_i *n_i* (*n_i*-1) for computing the Bose-Hubbard on-site interaction (without the
*U* prefactor.)

# Example

```jldoctest
julia> Hamiltonians.bose_hubbard_interaction(BoseFS{4,4}((2,1,1,0)))
2
julia> Hamiltonians.bose_hubbard_interaction(BoseFS{4,4}((3,0,1,0)))
6
```
"""
function bose_hubbard_interaction(b::BoseFS{<:Any,<:Any,A}) where A
    return bose_hubbard_interaction(Val(num_chunks(A)), b)
end

@inline function bose_hubbard_interaction(_, b::BoseFS)
    result = 0
    for (n, _, _) in occupied_modes(b)
        result += n * (n - 1)
    end
    return result
end

@inline function bose_hubbard_interaction(::Val{1}, b::BoseFS)
    # currently this ammounts to counting occupation numbers of modes
    chunk = chunks(b.bs)[1]
    matrixelementint = 0
    while !iszero(chunk)
        chunk >>>= (trailing_zeros(chunk) % UInt) # proceed to next occupied mode
        bosonnumber = trailing_ones(chunk) # count how many bosons inside
        # surpsingly it is faster to not check whether this is nonzero and do the
        # following operations anyway
        chunk >>>= (bosonnumber % UInt) # remove the counted mode
        matrixelementint += bosonnumber * (bosonnumber - 1)
    end
    return matrixelementint
end

function diagonal_element(h::HubbardReal1D, address::BoseFS)
    h.u * bose_hubbard_interaction(address) / 2
end

function get_offdiagonal(h::HubbardReal1D, add::BoseFS, chosen)
    naddress, onproduct = hopnextneighbour(add, chosen)
    return naddress, - h.t * sqrt(onproduct)
end


text:
import .Utils: pow_of_minus_one

####
#### OpTerm{AbstractFermiOp}
####

const FermiTerm = OpTerm{FermiOp}
const FermiSum = OpSum{FermiOp}
const AFermiTerm = OpTerm{<:AbstractFermiOp}
const DenseFermiTerm = DenseOpTerm{<:AbstractFermiOp}
const AFermiSum = OpSum{<:AbstractFermiOp}

macro fermi_str(str)
    return _op_term_macro_helper(FermiTerm, str)
end

####
#### Predicates
####

function IsApprox.ishermitian(ft::AFermiTerm, approx_test::AbstractApprox=Equal())
    for x in op_string(ft)
        if FermiOps.is_raise_lower(x)
            return false
        end
    end
    return IsApprox.ishermitian(ft.coeff, approx_test)
end

####
#### Algebra / mathematical operations
####

"""
    count_ladder_ops(v::AbstractVector)

Return the number of explicit raising and lowering operators in `v`.
"""
count_ladder_ops(v::AbstractVector) = count(FermiOps.is_raise_lower, v)

function fermi_phase_count(v1::AbstractVector, v2::AbstractVector)
    length(v1) == length(v2) || throw(DimensionMismatch("v1 an v2 differ in length"))
    _count = count_ladder_ops(v1)
    phase_count = 0
    @inbounds for i in eachindex(v1)
        if is_raise_lower(v1[i])
            _count -= 1
        end
        if !iseven(_count) && is_raise_lower(v2[i])
            phase_count += 1
        end
    end
    return phase_count
end

function AbstractOps.compute_phase(_::AbstractOps.DummyPhaseData, s1::AbstractArray{<:FermiOp}, s2::AbstractArray{<:FermiOp})
    phase_count = fermi_phase_count(s1, s2)
    return pow_of_minus_one(phase_count)
end

function Base.:^(ft::AFermiTerm, n::Integer)
    n < 0 && throw(DomainError(n))
    n == 0 && return one(ft)
    n == 1 && return ft
    ## Any +,-,0 sends the entire string to zero
    if any(x -> x === Raise || x === Lower || x === Zero, op_string(ft))
        return zero(ft)
    end
    ## E, N, I, are idempotent
    return strip_typeof(ft)(copy(op_string(ft)), ft.coeff^n)
end

#count_bodies() = count_bodies(op_string(ft))

function count_bodies(ft::OpTerm)
    n = sum(FermiOps.count_raise_lower, op_string(ft))
    return n ÷ 2
end

# This may account correctly for phase
# Looks like each excitation gives a sign flip. That is, each
# -,+ pair gives a sign flip. So odd number of excitations gives sign flip
# No flip: number op, double excitation, coulomb interact.
# a flip: excitation op, number-excitation
function Base.adjoint(ft::AFermiTerm)
    n_excitations = count_ladder_ops(op_string(ft)) ÷ 2
    new_coeff = iseven(n_excitations) ? conj(ft.coeff) : -conj(ft.coeff)
    return strip_typeof(ft)(adjoint.(op_string(ft)), new_coeff)
end

## TODO: improve efficiency
function Base.adjoint(ft::AFermiSum)
    terms = [adjoint(x) for x in ft]
    return strip_typeof(ft)(terms)
end


text:
function train(Xs::Signal, Ys::Signal, net::Network; iter=50, batch_size=10)
    errors = Float32[]
    for t in 1:iter
        batch_num = floor(Int64, (size(Xs, 2) / batch_size)) - 1
        Y = Ys[:, 1:2]
        for i in 1:batch_num
            ids = [1:batch_size;] + (i * batch_size)
            X   = Xs[:, ids]
            Y   = Ys[:, ids]

            backprop!(net, X, Y)
            update!(net)

            if (i % 500) == 0
                report(net, errors, Y, ids)
           end
        end
    end
    errors
end

function test(Xs::Signal, Ys::Signal, net::Network; iter=10, batch_size=10)
    miss = 0
    for t in 1:iter
        ids = [1:batch_size;] + (t * batch_size)
        X   = Xs[:, ids]
        Y   = Ys[:, ids]

        Z = forward!(net, X)

        for n in 1:batch_size
            ys, yn = findmax(Y[:, n])
            zs, zn = findmax(Z[:, n])
            if yn != zn
                miss += 1
            end
        end
    end
    1.0 - (miss ./ (iter * batch_size))
end

function report(net, errors, Y, ids)
    l = net.layers[end]
    cost = l.loss(l, Y, l.Z)
    Base.push!(errors, cost)
    println("$(ids[1]) ~ $(ids[end])", " ",
            printfp("cost", cost),     "\t\t|\t",
            report(net.layers[1]),     "\t|\t",
            report(net.layers[2])
    )
end

function report(l::Layer)
    join([printfp("E", l.ΔE),
          printfp("W", l.W),
          printfp("Z", l.Z),
          printfp("X", l.X)], "  ")
end

function report(l::SimpleRecurrentLayer)
    join([printfp("E",  l.ΔE),
          printfp("W",  l.W),
          printfp("Wh", l.Wh),
          printfp("Z",  l.Z),
          printfp("X",  l.X)], "  ")
end

function printfp(label::AbstractString, x::FloatingPoint)
    join([label, ":", @sprintf("%.2f", x)])
end

function printfp(label::AbstractString, x::Matrix{Float32})
    printfp(label, sum(x))
end


text:
⊆(x::Complex{Interval{T}}, y::Complex{Interval{T}}) where T = real(x) ⊆ real(y) && imag(x) ⊆ imag(y)

function ^(x::Complex{Interval{T}}, n::Integer) where {T}
    if n < 0
        return inv(x)^n
    end

    return Base.power_by_squaring(x, n)
end

function ^(x::Complex{Interval{T}}, y::Real) where {T}
    return exp(y*log(x))
end

function ^(x::Complex{Interval{T}}, y::Complex) where {T}
    return exp(y*log(x))
end

function ssqs(x::T, y::T,RND::RoundingMode) where T<:AbstractFloat
    k::Int = 0
    ρ = +(*(x,x,RND),*(y,y,RND),RND)
    if !isfinite(ρ) && (isinf(x) || isinf(y))
        ρ = convert(T, Inf)
    elseif isinf(ρ) || (ρ==0 && (x!=0 || y!=0)) || ρ<nextfloat(zero(T))/(2*eps(T)^2)
        m::T = max(abs(x), abs(y))
        k = m==0 ? m : exponent(m)
        xk, yk = ldexp(x,-k), ldexp(y,-k)
        ρ = +(*(xk,xk,RND),*(yk,yk,RND),RND)
    end
    ρ, k
end

# function ssqs(x::Interval{T}, y::Interval{T}) where T<:AbstractFloat
#     x = abs(x); y = abs(y)
#     ρl,kl = ssqs(inf(x),inf(y),RoundDown)
#     ρu,ku = ssqs(sup(x),sup(y),RoundUp)
#     ρl, kl, ρu, ku
# end

function sqrt_realpart(x::T, y::T,RND::RoundingMode) where T<:AbstractFloat
    # @assert x ≥ 0
    ρ, k = ssqs(x,y,RND)
    if isfinite(x) ρ = +(ldexp(x,-k),sqrt(ρ,RND),RND) end

    if isodd(k)
        k = div(k-1,2)
    else
        k = div(k,2)-1
        ρ *= 2
    end
    ρ = ldexp(sqrt(ρ,RND),k) #sqrt((abs(z)+abs(x))/2) without over/underflow
end

function sqrt(z::Complex{Interval{T}}) where T<:AbstractFloat
    x, y = reim(z)

    (inf(x) < 0 && inf(y) < 0 && sup(y) >= 0) && error("Interval lies across branch cut")
    x == y == 0 && return Complex(zero(x),y)

    (inf(x) < 0 && sup(x) > 0) && return sqrt((inf(x)..0) + im*y) ∪ sqrt((0..sup(x)) + im*y)

    absx = abs(x)

    absy= abs(y)
    ρl = sqrt_realpart(inf(absx),inf(absy),RoundDown)
    ρu = sqrt_realpart(sup(absx),sup(absy),RoundUp)

    ηu = sup(y)
    if isfinite(ηu)
        if ηu >= 0
            ηu_re = sqrt_realpart(inf(absx),ηu,RoundDown)
        else
            ηu_re = sqrt_realpart(sup(absx),ηu,RoundUp)
        end
        ηu = ηu_re ≠ 0 ?  /(ηu,2ηu_re,RoundUp) : zero(T)
    end

    ηl = inf(y)
    if ηl >= 0
        ηl_re = sqrt_realpart(sup(absx),ηl,RoundUp)
    else
        ηl_re = sqrt_realpart(inf(absx),ηl,RoundDown)
    end
    ηl = ηl_re ≠ 0 ? /(ηl,2ηl_re,RoundDown) : zero(T)

    ξ = ρl..ρu
    η = ηl..ηu

    if mid(x)<0
        ξ = flipsign(η,mid(η))
        η = flipsign((ρl..ρu),mid(η))
    end

    Complex(ξ,η)
end

function log(z::Complex{T}) where T<:Interval
    ρ = abs(z)
    θ = angle(z)

    return log(ρ) + im * θ
end

function abs2(z::Complex{T}) where T<:Interval
    return real(z)^2 + imag(z)^2
end

function abs(z::Complex{T}) where T<:Interval
    return sqrt(abs2(z))
end
#
# # \left( |x|^p \right)^{1/p}.
# function norm(z::Complex{T}, p=2) where T<:Interval
#     return (abs(z)^(p))^(1 / p)
# end


text:
struct Camera
    origin::AbstractVec
    lower_left_corner::AbstractVec
    horizontal::AbstractVec
    vertical::AbstractVec
    u::AbstractVec
    v::AbstractVec
    lens_radius::AbstractFloat

    function Camera(
        look_from::AbstractVec,
        look_to::AbstractVec,
        up::AbstractVec,
        vfov::AbstractFloat,
        aspect::AbstractFloat,
        aperture::AbstractFloat,
        focus_dist::AbstractFloat)

        Θ::AbstractFloat = vfov * π / 180
        half_height::AbstractFloat = tan(Θ / 2)
        half_width::AbstractFloat = aspect * half_height
        w::AbstractVec = normalize(look_from - look_to)
        u::AbstractVec = normalize(cross(up, w))
        v::AbstractVec = cross(w, u)
        lower_left_corner::AbstractVec = look_from - focus_dist * (half_width * u + half_height * v + w)
        horizontal::AbstractVec = 2.0 * focus_dist * half_width * u
        vertical::AbstractVec = 2.0 * focus_dist * half_height * v
        lens_radius::AbstractFloat = aperture / 2

        new(look_from, lower_left_corner, horizontal, vertical, u, v, lens_radius)
    end
end

function get_ray(camera::Camera, s::AbstractFloat, t::AbstractFloat)::Ray
    rd::AbstractVec = camera.lens_radius * random_in_unit_disk()
    offset::AbstractVec = camera.u * x(rd) + camera.v * y(rd)
    Ray(camera.origin + offset, camera.lower_left_corner + s * camera.horizontal + t * camera.vertical - (camera.origin + offset))
end

max_stars_repo_name: goooxu/JuliaRayTracing

max_issues_repo_name: goooxu/JuliaRayTracing

max_forks_repo_name: goooxu/JuliaRayTracing


text:
# # Term sparsity

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Sparsity/term_sparsity.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Sparsity/term_sparsity.ipynb)
# **Adapted from**: Example 3.5 of [WML20b]
#
# [WML20a] Wang, Jie, Victor Magron, and Jean-Bernard Lasserre.
# *TSSOS: A Moment-SOS hierarchy that exploits term sparsity*.
# arXiv preprint arXiv:1912.08899 (2020).
#
# [WML20b] Wang, Jie, Victor Magron, and Jean-Bernard Lasserre.
# *Chordal-TSSOS: a moment-SOS hierarchy that exploits term sparsity with chordal extension*.
# arXiv preprint arXiv:2003.03210 (2020).

using Test #src
using DynamicPolynomials
@polyvar x[1:3]

# We would like to find the minimum value of the polynomial

poly = x[1]^2 - 2x[1]*x[2] + 3x[2]^2 - 2x[1]^2*x[2] + 2x[1]^2*x[2]^2 - 2x[2]*x[3] + 6x[3]^2 + 18x[2]^2*x[3] - 54x[2]*x[3]^2 + 142x[2]^2*x[3]^2

# The minimum value of the polynomial can be found to be zero.

import CSDP
solver = CSDP.Optimizer
using SumOfSquares
function sos_min(sparsity)
    model = Model(solver)
    @variable(model, t)
    @objective(model, Max, t)
    con_ref = @constraint(model, poly - t in SOSCone(), sparsity = sparsity)
    optimize!(model)
    return value(t), moment_matrix(con_ref)
end

bound, ν = sos_min(NoSparsity())
@test bound ≈ 0 atol=1e-6 #src
bound

# We find the corresponding minimizer `(0, 0, 0)` by matching the moments
# of the moment matrix with a dirac measure centered at this minimizer.

extractatoms(ν, 1e-6)

# We can see below that the basis contained 6 monomials hence we needed to use 6x6 PSD matrix variables.

@test ν.basis.monomials == [x[1]*x[2], x[2]*x[3], x[1], x[2], x[3], 1] #src
ν.basis

# Using the monomial/term sparsity method of [WML20a] based on cluster completion, we find the same bound.

bound, ν = sos_min(MonomialSparsity())
@test bound ≈ 0 atol=1e-6 #src
bound

# Which is not suprising as no sparsity reduction could be performed.

@test length(ν.sub_moment_matrices) == 1 #src
@test ν.sub_moment_matrices[1].basis.monomials == [x[1]*x[2], x[2]*x[3], x[1], x[2], x[3], 1] #src
[sub.basis for sub in ν.sub_moment_matrices]

# Using the monomial/term sparsity method of [WML20b] based on chordal completion, the lower bound is smaller than 0.

bound, ν = sos_min(MonomialSparsity(ChordalCompletion()))
@test bound ≈ -0.00355 rtol=1e-3 #src
bound

# However, this bound was obtained with an SDP with 4 matrices of size 3x3.

@test length(ν.sub_moment_matrices) == 4                                  #src
@test ν.sub_moment_matrices[1].basis.monomials == [x[2]*x[3], x[2], x[3]] #src
@test ν.sub_moment_matrices[2].basis.monomials == [x[2]*x[3], x[2], 1]    #src
@test ν.sub_moment_matrices[3].basis.monomials == [x[1], x[2], 1]         #src
@test ν.sub_moment_matrices[4].basis.monomials == [x[1]*x[2], x[1], 1]    #src
[sub.basis for sub in ν.sub_moment_matrices]


text:
using LinearAlgebra
dt = 0.002
function flow(u::Array{Float64,2},s::Array{Float64,1})
    sigma, rho, beta = s
    x, y, z = view(u,1,:), view(u,2,:), view(u,3,:)
    return [sigma.*(y - x)  x.*(rho .- z) - y  x.*y - beta.*z]'
end
function flow(u::Array{Float64,1},s::Array{Float64,1})
    sigma, rho, beta = s
    x, y, z = u[1], u[2], u[3]
    return [sigma*(y - x), x*(rho - z) - y, x*y - beta*z]
end
function mag_flow(u::Array{Float64,2},s::Array{Float64,1})
    d, n = size(u)
    mags = zeros(n)
    for i = 1:n
    	mags[i] = norm(flow(u[:,i],s))
    end
    return mags
end
function mag_flow(u::Array{Float64,1},s::Array{Float64,1})
    return norm(flow(u,s))
end
function dmag_flow(u::Array{Float64,2},s::Array{Float64,1})
    sigma, rho, beta = s
    d, n = size(u)
    db = zeros(d, n)
    for i = 1:n
    	x, y, z = u[1,i],u[2,i],u[3,i]
    	term = [sigma*(y - x), x*(rho - z) - y,
    		x*y - beta*z]
    	b = norm(term)
    	db[1,i] = dot([-sigma, (rho - z), y], term./b)
    	db[2,i] = dot([sigma, -1, x], term./b)
    	db[3,i] = dot([0., -x, -beta], term./b)
    end
    return db
end
function dmag_flow(u::Array{Float64,1},s::Array{Float64,1})
    sigma, rho, beta = s
    x, y, z = u[1], u[2], u[3]

    term = [sigma*(y - x), x*(rho - z) - y,
    		x*y - beta*z]
    b = norm(term)
    dxb2 = dot([-sigma, (rho - z), y], term./b)
    dyb2 = dot([sigma, -1, x], term./b)
    dzb2 = dot([0., -x, -beta], term./b)
    return [dxb2, dyb2, dzb2]
end
function dunit_flow(u::Array{Float64,1},s::Array{Float64,1})
    sigma, rho, beta = s
    x, y, z = u[1], u[2], u[3]
    F = [sigma*(y - x), x*(rho - z) - y, x*y - beta*z]
    b2 = dot(F,F)
    db = dmag_flow(u, s)
    dF = -(db'/b2).*[F F F]
    dF = dF .+ (1.0/sqrt(b2))*dflow(u,s)
    return dF
end

function step(u0, s, n)
    u_trj = zeros(3, n+1)
    u_trj[:,1] = u0
    n = n+1
    for i = 2:n
        u = u_trj[:,i-1]
        k1 = dt*flow(u, s)
        k2 = dt*flow(u .+ k1, s)
        u_trj[:,i] = @. u + (k1 + k2)/2
    end
    return u_trj
end
function next(u, s)
    k1 = dt*flow(u, s)
    k2 = dt*flow(u .+ k1, s)
    u_next = @. u + (k1 + k2)/2
    return u_next
end
function dflow(u::Array{Float64,2}, s::Array{Float64,1})
    sigma, rho, beta = s
    d, n = size(u)
    x = view(u,1,:)
    y = view(u,2,:)
    z = view(u,3,:)
    du = zeros(n, d, d)
    @. du[:,1,1] = -sigma
    @. du[:,1,2] = sigma
    @. du[:,2,1] = (rho - z)
    @. du[:,2,2] = -1.0
    @. du[:,2,3] = -x
    @. du[:,3,1] = y
    @. du[:,3,2] = x
    @. du[:,3,3] = -beta
    return reshape([du[:,:,1]'; du[:,:,2]';
                	du[:,:,3]'], d, d, n)
end
function dflow(u::Array{Float64,1},s::Array{Float64,1})
    du = zeros(3,3)
    x, y, z = u[1], u[2], u[3]
    sigma, rho, beta = s
    du[1,1] = -sigma
    du[1,2] = sigma
    du[2,1] = (rho - z)
    du[2,2] = -1.0
    du[2,3] = -x
    du[3,1] = y
    du[3,2] = x
    du[3,3] = -beta
    return du
end
function dstep(u::Array{Float64,2}, s::Array{Float64,1})
    d, n = size(u)
    du = zeros(d, d, n)
    dui = zeros(d,d)
    for i = 1:n
        ui = u[:,i]
        k1 = dt*flow(ui, s)
        dk1 = dt*dflow(ui, s)
        dk2 = dt*dflow(ui .+ k1, s)
        dui .= 1.0*I(d) .+ 0.5*(dk1 .+ dk2*(1.0*I(d) .+ dk1))
        du[:,:,i] .= dui
    end
    return du
end
function dstep(u::Array{Float64,1}, s::Array{Float64,1})
    d = 3
    du = zeros(d, d)
    k1 = dt*flow(u, s)
    dk1 = dt*dflow(u, s)
    dk2 = dt*dflow(u .+ k1, s)
    du = 1.0*I(d) .+ 0.5*(dk1 .+ dk2*(1.0*I(d) .+ dk1))
    return du
end

function d2flow(u::Array{Float64,2}, s::Array{Float64,1})
    n = size(u)[2]
    ddu = zeros(3,3,3,n)
    for i = 1:n
        ddu[2,3,1,i] = -1.0
        ddu[3,2,1,i] = 1.0
        ddu[3,1,2,i] = 1.0
        ddu[2,1,3,i] = -1.0
    end
    return ddu
end
function d2flow(u::Array{Float64,1}, s::Array{Float64,1})
    ddu = zeros(3,3,3)
    ddu[2,3,1] = -1.0
    ddu[3,2,1] = 1.0
    ddu[3,1,2] = 1.0
    ddu[2,1,3] = -1.0
    return ddu
end
function d2step(u::Array{Float64,2}, s::Array{Float64,1})
    d, n = size(u)
    ddu = zeros(d, d, d, n)
    k1 = dt*flow(u, s)
    dk1 = dt*dflow(u,s)
    dk2 = dt*dflow(u .+ k1, s)
    d2 = dt*d2flow(u, s)
    d2p = permutedims(d2,[1, 3, 2, 4])
    A = zeros(d,d,d)
    for k = 1:n
        for j = 1:d
        	A[:,:,j] .= d2p[:,:,j,k]*(1.0*I(d) .+ dk1[:,:,k])
        end
        A = permutedims(A, [1,3,2])
        for i = 1:d
        		ddu[:,:,i,k] = 0.5*(d2[:,:,i,k] .+
        				dk2[:,:,k]*d2[:,:,i,k] .+ A[:,:,i]*
        				(1.0*I(d) .+ dk1[:,:,k]))
        end
    end
    return ddu
end
function d2step(u::Array{Float64,1}, s::Array{Float64,1})
    ddu = zeros(3, 3, 3)
    k1 = dt*flow(u, s)
    dk1 = dt*dflow(u,s)
    dk2 = dt*dflow(u .+ k1, s)
    d2 = dt*d2flow(u, s)
    d2p = permutedims(d2,[1, 3, 2])
    A = zeros(3,3,3)
    for j = 1:3
        A[:,:,j] .= d2p[:,:,j]*(1.0*I(3) .+ dk1)
    end
    A = permutedims(A, [1,3,2])
    for i = 1:3
        ddu[:,:,i] = 0.5*(d2[:,:,i] .+ dk2*d2[:,:,i] .+ A[:,:,i]*
        					  (1.0*I(3) .+ dk1))
    end

    return ddu
end
function pert(u::Array{Float64,2},s::Array{Float64,1})
    n = size(u)[2]
    dsflow(x) = [0., x, 0.]
    pp = zeros(3, n)
    for i = 1:n
    	ui = u[:,i]
    	k1 = dt*flow(ui, s)
    	dk2 = dt*dflow(ui .+ k1,s)
    	dsk1 = dt*dsflow(ui[1])
    	dsk2 = dt*dsflow(ui[1] .+ k1[1]) .+
    		   dk2*dsk1
    	pp[:,i] = 0.5*(dsk1 + dsk2)
    end
    return pp
end
function pert(u::Array{Float64,1},s::Array{Float64,1})
    dsflow(x) = [0., x, 0.]
    k1 = dt*flow(u, s)

    dk2 = dt*dflow(u .+ k1,s)
    dsk1 = dt*dsflow(u[1])
    dsk2 = dt*dsflow(u[1] .+ k1[1]) .+
    		   dk2*dsk1
    pp = 0.5*(dsk1 + dsk2)
    return pp
end
function dpert_next(u_arr::Array{Float64,2},s::Array{Float64,1})
    d, n = size(u_arr)
    dpp = zeros(3,3,n)
    ddsflow = zeros(3,3)
    ddsflow[2,1] = 1.0
    dsflow(x) = [0., x, 0.]
    for j = 1:n
    	u = u_arr[:,j]
    	k1 = dt*flow(u, s)
    	dk1 = dt*dflow(u,s)
    	dk2 = dt*dflow(u .+ k1,s)
    	ddk2 = dt*d2flow(u .+ k1, s)
    	dsk1 = dt*dsflow(u[1])
    	ddsk1 = dt*ddsflow
    	d2 = dt*d2flow(u, s)
        d2p = permutedims(d2,[1, 3, 2])
        A = zeros(3,3,3)
        for j = 1:3
        	A[:,:,j] .= d2p[:,:,j]*(1.0*I(3) .+ dk1)
        end
        A = permutedims(A, [1,3,2])
    	ddsk2 = dt*ddsflow*(1.0*I(3) .+ dk1) .+
    		dk2*ddsk1
    	for i = 1:3
    		ddsk2[:,i] .+= A[:,:,i]*dsk1
    	end
    	dpp[:,:,j] .= 0.5*(ddsk1 + ddsk2)

      end
    return dpp
end
function dpert_next(u::Array{Float64,1},s::Array{Float64,1})
    dpp = zeros(3,3)
    ddsflow = zeros(3,3)
    ddsflow[2,1] = 1.0

    k1 = dt*flow(u, s)
    dk1 = dt*dflow(u,s)

    dk2 = dt*dflow(u .+ k1,s)
    ddk2 = dt*d2flow(u .+ k1, s)

    dsflow(x) = [0., x, 0.]
    dsk1 = dt*dsflow(u[1])
    ddsk1 = dt*ddsflow

    d2 = dt*d2flow(u, s)
    d2p = permutedims(d2,[1, 3, 2])
    A = zeros(3,3,3)
    for j = 1:3
        A[:,:,j] .= d2p[:,:,j]*(1.0*I(3) .+ dk1)
    end
    A = permutedims(A, [1,3,2])

    ddsk2 = dt*ddsflow*(1.0*I(3) .+ dk1) .+
    		dk2*ddsk1
    for j = 1:3
    	ddsk2[:,j] .+= A[:,:,j]*dsk1
    end
    dpp .= 0.5*(ddsk1 + ddsk2)
    return dpp
end
function dpert(u::Array{Float64, 1}, s::Array{Float64,1})
    du = dstep(u,s)
    dpp_next = dpert_next(u,s)
    dpp = du'\dpp_next'
    return dpp'
end
function dpert(u::Array{Float64, 2}, s::Array{Float64,1})
    du = dstep(u,s)
    dpp_next = dpert_next(u,s)
    d, n = size(u)
    dpp = zeros(3,3,n)
    for i = 1:n
    	dpp[:,:,i] .= du[:,:,i]'\dpp_next[:,:,i]'
    	dpp[:,:,i] .= dpp[:,:,i]'
    end
    return dpp
end

### AD functions
function lorenz63_rhs_ad(du, u, s, t)
    du[1] = s[1]*(u[2] - u[1])
    du[2] = u[1]*(s[2] - u[3]) - u[2]
    du[3] = u[1]*u[2] - s[3]*u[3]
end
function lorenz63_ad(u0, s, n)
    t = n*dt
    prob = ODEProblem(lorenz63_rhs_ad, u0, (0.,t), s)
    sol = Array(solve(prob, Tsit5(), saveat=dt))
    return sol[:,end]
end
function obj_fun(u0, s)
    prob = ODEProblem(lorenz63_rhs_ad, u0, (0.,1.1), s)
    #_prob = remake(prob,u0=u0,p=s)
    sol = solve(prob, Tsit5(), saveat=0.005)
    sum(sol[3,:])/size(sol,2)
end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
using NumericExtensions
using Base.Test

## data

a1 = 2 * rand(8) + 1.0
a2 = 2 * rand(8, 7) + 1.0
a3 = 2 * rand(8, 7, 6) + 1.0
a4 = 2 * rand(8, 7, 6, 5) + 1.0

b1 = 2 * rand(8) - 1.0
b2 = 2 * rand(8, 7) - 1.0
b3 = 2 * rand(8, 7, 6) - 1.0
b4 = 2 * rand(8, 7, 6, 5) - 1.0

c1 = 2 * rand(8) - 1.0
c2 = 2 * rand(8, 7) - 1.0
c3 = 2 * rand(8, 7, 6) - 1.0
c4 = 2 * rand(8, 7, 6, 5) - 1.0

ua1 = view(a1, 1:6)
ub1 = view(b1, 2:7)
va1 = view(a1, 1:2:7)
vb1 = view(b1, 3:1:6)

ua2 = view(a2, 1:6, 1:5)
ub2 = view(b2, 2:7, 2:6)
va2 = view(a2, 1:2:7, 1:5)
vb2 = view(b2, 3:1:6, 2:6)

ua3 = view(a3, 1:6, 1:5, 1:3)
ub3 = view(b3, 2:7, 2:6, 3:5)
va3 = view(a3, 1:2:7, 1:5, 1:3)
vb3 = view(b3, 3:1:6, 2:6, 3:5)

ua4 = view(a4, 1:6, 1:5, 1:3, 1:4)
ub4 = view(b4, 2:7, 2:6, 3:5, 2:5)
va4 = view(a4, 1:2:7, 1:5, 1:3, 1:4)
vb4 = view(b4, 3:1:6, 2:6, 3:5, 2:5)

## test cases

arrs_a = {a1, a2, a3, a4, ua1, va1, ua2, va2, ua3, va3, ua4, va4}
arrs_b = {b1, b2, b3, b4, ub1, vb1, ub2, vb2, ub3, vb3, ub4, vb4}

## unary cases

println("  -- unary functions")

unaryfs = [(Negate, -, negate!),
           (AbsFun, abs, abs!),
           (Abs2Fun, abs2, abs2!),
           (SqrtFun, sqrt, sqrt!),
           (FloorFun, floor, floor!),
           (CeilFun, ceil, ceil!),
           (ExpFun, exp, exp!),
           (LogFun, log, log!)]

for a in arrs_a, (Op, vf, vf!) in unaryfs
    r = vf(copy(a))
    @test_approx_eq map(Op(), a) r

    y = zeros(size(r))
    vf!(y, a)
    @test_approx_eq y r

    ac = copy(a)
    vf!(ac)
    @test_approx_eq ac r
end

## binary cases

println("  -- binary functions")

binaryfs = [(Add, +, add!),
            (Subtract, -, subtract!),
            (Multiply, .*, multiply!),
            (Divide, ./, divide!)]

for (a, b) in zip(arrs_a, arrs_b), (Op, vf, vf!) in binaryfs
    r = vf(copy(b), copy(a))
    @test_approx_eq map(Op(), b, a) r
    @test_approx_eq map(Op(), a, 2.0) vf(copy(a), 2.0)
    @test_approx_eq map(Op(), 2.0, b) vf(2.0, copy(b))

    y = zeros(size(r))
    vf!(y, b, a)
    @test_approx_eq y r

    bc = copy(b)
    vf!(bc, a)
    @test_approx_eq bc r

    bc = copy(b)
    vf!(bc, 2.0)
    @test_approx_eq bc vf(copy(b), 2.0)
end

## mapdiff cases

println("  -- mapdiff functions")

for (a, b) in zip(arrs_a, arrs_b)
    ac = copy(a)
    bc = copy(b)
    v = 2.6

    @test_approx_eq absdiff(a, b) abs(ac - bc)
    @test_approx_eq absdiff(a, v) abs(ac - v)
    @test_approx_eq absdiff(v, b) abs(v - bc)

    @test_approx_eq sqrdiff(a, b) abs2(ac - bc)
    @test_approx_eq sqrdiff(a, v) abs2(ac - v)
    @test_approx_eq sqrdiff(v, b) abs2(v - bc)

    @test_approx_eq absdiff!(zeros(size(a)), a, b) abs(ac - bc)
    @test_approx_eq absdiff!(zeros(size(a)), a, v) abs(ac - v)
    @test_approx_eq absdiff!(zeros(size(a)), v, b) abs(v - bc)

    @test_approx_eq sqrdiff!(zeros(size(a)), a, b) abs2(ac - bc)
    @test_approx_eq sqrdiff!(zeros(size(a)), a, v) abs2(ac - v)
    @test_approx_eq sqrdiff!(zeros(size(a)), v, b) abs2(v - bc)
end

## ternary cases

println("  -- ternary functions")

@test_approx_eq map(FMA(), a1, b1, c1) a1 + b1 .* c1

v = 2.5

for (a, b) in zip(arrs_a, arrs_b)
    ac = copy(a)
    bc = copy(b)

    @test_approx_eq map(FMA(), v, a, b) v + ac .* bc
    @test_approx_eq map(FMA(), a, v, b) ac + v .* bc
    @test_approx_eq map(FMA(), a, b, v) ac + bc .* v

    fma!(ac, a, b)
    @test_approx_eq ac a + a .* b

    ac = copy(a)
    fma!(ac, b, v)
    @test_approx_eq ac a + b .* v
end


text:
using InteractiveUtils

function nonperiodic(bd::Billiard)
    toplot = Obstacle{eltype(bd)}[]
    for obst in bd
        typeof(obst) <: PeriodicWall && continue
        push!(toplot, obst)
    end
    return toplot
end

periodicerror() = throw(ArgumentError(
"The billiard must be periodic, i.e. has at least two `PeriodicWall` obstacles."
))

"""
    plot(bd::Billiard; ax = (figure(); gca()))
Plot all obstacles in `bd` using the default arguments, set
`xlim` and `ylim` to be 20% larger than `cellsize` and
set the axis aspect ratio to equal.

    plot(bd::Billiard, xmin, ymin, xmax, ymax;
         hexagonal = false, ax = (figure(); gca()))
Plot the given **periodic** billiard `bd`, repeatedly
plotting from `(xmin, ymin)` to `(xmax, ymax)`.
Works for either rectangular periodic billiards, or hexagonal ones. Use
keyword `hexagonal` to denote which one you want.

    plot(bd::Billiard, xt::Vector, yt::Vector; hexagonal = false,
         ax = (figure(); gca()), plot_orbit = true, orbit_color = "C0")
Plot the given billiard `bd` and an acompanying orbit resulting from
[`timeseries!`](@ref) using the limits defined by `xt` and `yt`.
Works for both periodic and normal billiards.
"""
function plot(bd::Billiard) end

function plot(bd::Billiard{T};
    ax = (PyPlot.figure(); PyPlot.gca())) where {T}
    PyPlot.sca(ax)
    for obst in bd; plot(obst); end
    xmin, ymin, xmax, ymax = cellsize(bd)
    dx = xmax - xmin; dy = ymax - ymin
    ax.set_aspect("equal")
    if !isinf(xmin) && !isinf(xmax)
        PyPlot.xlim(xmin - 0.1dx, xmax + 0.1dx)
    end
    if !isinf(ymin) && !isinf(ymax)
        PyPlot.ylim(ymin - 0.1dy, ymax + 0.1dy)
    end
    return nothing
end

function plot(bd::Billiard, xmin, ymin, xmax, ymax;
    hexagonal = false, ax = (PyPlot.figure(); PyPlot.gca()))

    isperiodic(bd) || periodicerror()
    PyPlot.sca(ax)

    n = count(x -> typeof(x) <: PeriodicWall, bd)
    if hexagonal
        n != 6 && throw(ArgumentError("Hexagonally periodic billiards have "*
        "exactly 6 periodic walls arranged as a perfect hexagon. Use the "*
        "function `billiard_polygon(6, r, R, setting = \"periodic\")`."))
        plot_periodic_hexagon(bd, xmin, ymin, xmax, ymax)
    else
        n ∉ (2, 4) && throw(ArgumentError(
        "Rectangular periodic billiards must have "*
        "exactly 2 or 4 periodic walls."))
        plot_periodic_rectangle(bd, xmin, ymin, xmax, ymax)
    end

    PyPlot.xlim(xmin, xmax)
    PyPlot.ylim(ymin, ymax)
    return nothing
end

function plot(bd, xt::AbstractVector, yt::AbstractVector;
    hexagonal = false, ax = (PyPlot.figure(); PyPlot.gca()),
    plot_orbit = true, orbit_color = "C0")

    if isperiodic(bd)
        xmin = minimum(xt); xmax = maximum(xt)
        ymin = minimum(yt); ymax = maximum(yt)

        plot(bd, xmin, ymin, xmax, ymax; hexagonal = hexagonal, ax = ax)
    else
        plot(bd; ax = ax)
    end

    if plot_orbit
        PyPlot.sca(ax)
        ax.scatter(xt[1], yt[1], color = "gray", s = 20.0, zorder = 99)
        ax.plot(xt, yt, color = orbit_color, zorder = 1)
    end

    cellxmin, cellymin, cellxmax, cellymax = cellsize(bd)
    dcx = cellxmax - cellxmin
    dcy = cellymax - cellymin

    if isperiodic(bd)
        if hexagonal
            PyPlot.xlim(xmin - dcx/2, xmax + dcx/2)
            PyPlot.ylim(ymin - dcy/2, ymax + dcy/2)
            PyPlot.gca().set_aspect("equal")
        else
            rr = xmin < 0 ? (div(xmin, dcx)-1) : (div(xmin, dcx))
            ll = ymin < 0 ? (div(ymin, dcy)-1) : (div(ymin, dcy))

            PyPlot.xlim(cellxmin + rr*dcx,
            cellxmin + (div(xmax, dcx)+1)*dcx)
            PyPlot.ylim(cellymin + ll*dcy,
            cellymin + (div(ymax, dcy)+1)*dcy)
            PyPlot.gca().set_aspect("equal")
        end
    end
    return nothing
end

function plot_periodic_rectangle(bd, xmin, ymin, xmax, ymax)

    # Cell limits:
    cellxmin, cellymin, cellxmax, cellymax = cellsize(bd)
    dcx = cellxmax - cellxmin
    dcy = cellymax - cellymin

    # Find displacement vectors
    dx = (floor((xmin - cellxmin)/dcx):1:ceil((xmax - cellxmax)/dcx)) * dcx
    dy = (floor((ymin - cellymin)/dcy):1:ceil((ymax - cellymax)/dcy)) * dcy

    # Plot displaced Obstacles
    toplot = nonperiodic(bd)
    for x in dx
        for y in dy
            disp = SVector(x,y)
            for obst in toplot
                plot(translate(obst, disp))
            end
        end
    end
end

function plot_periodic_hexagon(bd, xmin, ymin, xmax, ymax)

    v = 1
    while !(typeof(bd[v]) <: PeriodicWall); v += 1; end
    space = norm(bd[v].sp - bd[v].ep)*√3

    basis_a = space*SVector(0.0, 1.0)
    basis_b = space*SVector(√3/2, 1/2)
    basis_c = space*SVector(√3, 0.0)

    PyPlot.xlim(xmin - space/2, xmax + space/2)
    PyPlot.ylim(ymin, ymax)
    PyPlot.gca().set_aspect("equal")

    # Cell limits:
    cellxmin, cellymin, cellxmax, cellymax = cellsize(bd)
    dcx = cellxmax - cellxmin
    dcy = cellymax - cellymin
    jmin = Int((ymin - cellymin - dcy/2)÷space) - 1
    jmax = Int((ymax - cellymax + dcy/2)÷space) + 1
    imin = Int((xmin - cellxmin - dcx/2)÷(√3*space)) - 1
    imax = Int((xmax - cellxmax + dcx/2)÷(√3*space)) + 1

    obstacles = nonperiodic(bd)
    for d in obstacles
        for j ∈ jmin:jmax
            for i ∈ imin:imax
                plot(translate(d, j*basis_a + i*basis_c))
                plot(translate(d, j*basis_a + i*basis_c + basis_b))
            end
        end
    end
end


text:
txt = """Given\$a\$txt\$file\$of\$many\$lines,\$where\$fields\$within\$a\$line\$
are\$delineated\$by\$a\$single\$'dollar'\$character,\$write\$a\$program
that\$aligns\$each\$column\$of\$fields\$by\$ensuring\$that\$words\$in\$each\$
column\$are\$separated\$by\$at\$least\$one\$space.
Further,\$allow\$for\$each\$word\$in\$a\$column\$to\$be\$either\$left\$
justified,\$right\$justified,\$or\$center\$justified\$within\$its\$column."""

# left/right/center justification of strings:
ljust(s, width) = s * " "^max(0, width - length(s))
rjust(s, width) = " "^max(0, width - length(s)) * s
function center(s, width)
  pad = width - length(s)
  if pad <= 0
    return s
  else
    pad2 = div(pad, 2)
    return " "^pad2 * s * " "^(pad - pad2)
  end
end

parts = [split(rstrip(line, '$'), '$') for line in split(txt, '\n')]

max_widths = zeros(Int, maximum(length, parts))
for line in parts
  for (i, word) in enumerate(line)
    max_widths[i] = max(max_widths[i], length(word))
  end
end
max_widths += 1 # separate cols by at least one space

for (label, justify) in (("Left", ljust), ("Right",rjust), ("Center",center))
  println(label, " column-aligned output:")
  for line in parts
    for (j, word) in enumerate(line)
      print(justify(word, max_widths[j]))
    end
    println()
  end
  println("-"^sum(max_widths))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module Theon

using GeometryTypes
using JSON

end # module


text:
using Test
# using Revise
using PolynomialBasis

PB = PolynomialBasis

function allequal(v1, v2)
    return all(v1 .≈ v2)
end

function allequal(v1, v2, tol)
    np = length(v1)
    f = length(v2) == np
    return f && all([isapprox(v1[i], v2[i], atol = tol) for i = 1:np])
end

coeffs = reshape([1.0], 1, 1)
b0 = PB.LagrangeBasis(0)
b1 = PB.LagrangeBasis(1)
tp1 = PB.LagrangeTensorProductBasis(2, 1)

@test_throws AssertionError PB.InterpolatingPolynomial(coeffs, b1)
ip = PB.InterpolatingPolynomial(coeffs, b0)

@test allequal(ip.coeffs, coeffs)

float_type = typeof(0.0)
ip = PB.InterpolatingPolynomial(float_type, 2, b1)
@test size(ip.coeffs) == (2, 2)

ip = PB.InterpolatingPolynomial(3, tp1)
@test size(ip.coeffs) == (3, 4)

ip = PB.InterpolatingPolynomial(1, 1, 2; start = 0.0)
@test allequal(ip.basis(0.0), [1.0, 0.0, 0.0])
@test allequal(ip.basis(0.5), [0.0, 1.0, 0.0])
@test allequal(ip.basis(1.0), [0.0, 0.0, 1.0])

v = [1.0 2.0 3.0]
PB.update!(ip, v)
@test allequal(ip.coeffs, v)
@test_throws DimensionMismatch PB.update!(ip, 1:4)

ip = PB.InterpolatingPolynomial(1, 2, 2)
PB.update!(ip, repeat([1.0, 2.0, 3.0], 3))
@test ip(0.0, -1.0) ≈ 1.0
@test ip(1.0, 0.0) ≈ 2.0
@test ip(-1.0, 1.0) ≈ 3.0
@test ip([-1.0, 1.0]...) ≈ 3.0

ip = PB.InterpolatingPolynomial(2, 2, 2)
coeffs = [
    1.0 1.0 1.0 2.0 2.0 2.0 3.0 3.0 3.0
    5.0 5.0 5.0 7.0 7.0 7.0 9.0 9.0 9.0
]
PB.update!(ip, coeffs)
@test allequal(ip(-1.0, 0.0), [1.0, 5.0])
@test allequal(ip(0.0, -1.0), [2.0, 7.0])
@test allequal(ip(1.0, 1.0), [3.0, 9.0])

P = PB.InterpolatingPolynomial(1, 1, 2)
f(x) = 3x^2 + 2x - 5
df(x) = 6x + 2
coeffs = f.(P.basis.points)
PB.update!(P, coeffs)
@test all(PB.gradient(P, [-1.0]) .≈ [df(-1.0)])
@test all(PB.gradient(P, [0.19]) .≈ [df(0.19)])

P = PB.InterpolatingPolynomial(1, 2, 2)
f(x, y) = 3x^2 * y + 2y^2 + 2x * y
dfx(x, y) = 6 * x * y + 2 * y
dfy(x, y) = 3x^2 + 4y + 2x
coeffs = [f((P.basis.points[:, i])...) for i = 1:9]
PB.update!(P, coeffs)
@test allequal(PB.gradient(P, [0.1, 0.2]), [dfx(0.1, 0.2) dfy(0.1, 0.2)])
@test allequal(PB.gradient(P, 1, [0.7, -0.3]), dfx(0.7, -0.3))
@test allequal(PB.gradient(P, 2, [-0.5, 0.75]), dfy(-0.5, 0.75))
@test allequal(PB.gradient(P, [0.1, 0.45]), [dfx(0.1, 0.45) dfy(0.1, 0.45)])

poly = PB.InterpolatingPolynomial(1, 1, 3)
f(x) = 7x^3 + 2x^2
coeffs = f.(poly.basis.points)
update!(poly, coeffs)
testh(x) = 42x + 4
@test all([
    allequal(PB.hessian(poly, p)[1], testh(p)) for p in poly.basis.points
])

poly = PB.InterpolatingPolynomial(1, 2, 2)
f(x) = 12x[1]^2 + 7x[1] * x[2] + 11x[2]^2
h = [24, 7, 22]
coeffs = mapslices(f, Array(poly.basis.points), dims = 1)
update!(poly, coeffs)
p = poly.basis.points
@test all([allequal(PB.hessian(poly, p[:, i]), h, 1e2eps()) for i = 1:9])

poly = PB.InterpolatingPolynomial(1, 2, 3)
f(x) = x[1]^3 + 2x[1]^2 * x[2] + 18.0
tp3h(x) = [6x[1] + 4x[2], 4x[1], 0.0]
p = poly.basis.points
coeffs = mapslices(f, Array(p), dims = 1)
update!(poly, coeffs)
@test all([
    allequal(PB.hessian(poly, p[:, i]), tp3h(p[:, i]), 1e3eps()) for i = 1:16
])


text:
doc"""
Stability Counter
"""
type StabilityCounter
    counter::Int64
    ishalt::Function
    update::Function

    function StabilityCounter(halt=150::Int64, abs_tol=1e-12, rel_tol=1e-09)
        this = new()
        this.counter = 0
        this.ishalt = () -> this.counter >= halt

        function update(ΔE::Float64)
            if ΔE <= abs_tol #abs(ΔE) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)
                this.counter += 1
            else
                this.counter = 0
            end
        end
        this.update = update
        return this
    end
end


text:
@testset "Q+μD" begin
    # This was removed (not used)
end


text:
""" Exports Operational Model Results to a path"""
function write_model_results(results::Dict, save_path::String)

    if !isdir(save_path)
        @error("Specified path is not valid. Run write_results to save results.")
    end

    for (k,v) in results[:vars]
        file_path = joinpath(save_path, "$(k).feather")
        Feather.write(file_path, v)
    end

    write_optimizer_log(results[:optimizer_log], save_path)
    println("Files written to $save_path folder.")

    return

end

function write_optimizer_log(optimizer_log::Dict{Symbol, Any}, save_path::AbstractString)

    optimizer_log[:termination_status] = Int(optimizer_log[:termination_status])
    optimizer_log[:primal_status] = Int(optimizer_log[:primal_status])
    optimizer_log[:dual_status] = Int(optimizer_log[:dual_status])
    optimizer_log[:solve_time] = optimizer_log[:solve_time]

    df = DataFrames.DataFrame(optimizer_log)
    file_path = joinpath(save_path,"optimizer_log.feather")
    Feather.write(file_path, df)

    return

end


text:
let
    using Random
    Random.seed!(0)
    p = 3
    n = 12
    Ut = randn(p,n)
    Vt = randn(p,n)
    x = randn(n)
    A = egrss.full(Ut,Vt)
    yref = A*x
    y = egrss.symv(Ut,Vt,x);
    @test ( isapprox(y,yref) )
end

max_stars_repo_name: martinandersen/EGRSS

max_issues_repo_name: martinandersen/EGRSS

max_forks_repo_name: martinandersen/EGRSS


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

function SCIPpresolComp(elem1, elem2)
    ccall((:SCIPpresolComp, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPpresolCompName(elem1, elem2)
    ccall((:SCIPpresolCompName, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPpresolGetData(presol)
    ccall((:SCIPpresolGetData, libscip), Ptr{SCIP_PRESOLDATA}, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolSetData(presol, presoldata)
    ccall((:SCIPpresolSetData, libscip), Cvoid, (Ptr{SCIP_PRESOL}, Ptr{SCIP_PRESOLDATA}), presol, presoldata)
end

function SCIPpresolGetName(presol)
    ccall((:SCIPpresolGetName, libscip), Cstring, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolGetDesc(presol)
    ccall((:SCIPpresolGetDesc, libscip), Cstring, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolGetPriority(presol)
    ccall((:SCIPpresolGetPriority, libscip), Cint, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolGetMaxrounds(presol)
    ccall((:SCIPpresolGetMaxrounds, libscip), Cint, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolGetTiming(presol)
    ccall((:SCIPpresolGetTiming, libscip), SCIP_PRESOLTIMING, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolSetTiming(presol, timing)
    ccall((:SCIPpresolSetTiming, libscip), Cvoid, (Ptr{SCIP_PRESOL}, SCIP_PRESOLTIMING), presol, timing)
end

function SCIPpresolIsInitialized(presol)
    ccall((:SCIPpresolIsInitialized, libscip), UInt32, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolGetSetupTime(presol)
    ccall((:SCIPpresolGetSetupTime, libscip), Cdouble, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolGetTime(presol)
    ccall((:SCIPpresolGetTime, libscip), Cdouble, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolGetNFixedVars(presol)
    ccall((:SCIPpresolGetNFixedVars, libscip), Cint, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolGetNAggrVars(presol)
    ccall((:SCIPpresolGetNAggrVars, libscip), Cint, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolGetNChgVarTypes(presol)
    ccall((:SCIPpresolGetNChgVarTypes, libscip), Cint, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolGetNChgBds(presol)
    ccall((:SCIPpresolGetNChgBds, libscip), Cint, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolGetNAddHoles(presol)
    ccall((:SCIPpresolGetNAddHoles, libscip), Cint, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolGetNDelConss(presol)
    ccall((:SCIPpresolGetNDelConss, libscip), Cint, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolGetNAddConss(presol)
    ccall((:SCIPpresolGetNAddConss, libscip), Cint, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolGetNUpgdConss(presol)
    ccall((:SCIPpresolGetNUpgdConss, libscip), Cint, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolGetNChgCoefs(presol)
    ccall((:SCIPpresolGetNChgCoefs, libscip), Cint, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolGetNChgSides(presol)
    ccall((:SCIPpresolGetNChgSides, libscip), Cint, (Ptr{SCIP_PRESOL},), presol)
end

function SCIPpresolGetNCalls(presol)
    ccall((:SCIPpresolGetNCalls, libscip), Cint, (Ptr{SCIP_PRESOL},), presol)
end