text:
"""
    Module Basis2D
Includes DG core functions for evaluating 2D bases on quads
"""
module Basis2D

export vandermonde2D, grad_vandermonde2D
export chebyshev2D, equidistant_nodes2D, gauss_quad2D
export dPhi2D

import Basis1D
using LinearAlgebra, SparseArrays

"""
    vandermonde2D(p, x, y)
Return the 2D vandermonde matrix - the order 0:p Legendre polynomials
evaluated at 2D nodes (x,y)
"""
function vandermonde2D(p, x, y)
    V = Array{Float64,4}(undef,length(x),length(y),p+1,p+1)
    lx = Basis1D.vandermonde1D(p, x)
    ly = Basis1D.vandermonde1D(p, y)
    for ipy = 1:p+1
        for ipx = 1:p+1
            V[:,:,ipx,ipy] = lx[:,ipx]*(ly[:,ipy]')
        end
    end
    return V
end

"""
    grad_vandermonde2D(p, x, y)
Return the gradient of the order 0:p Legendre polynomials
evaluated at 2D nodes (x,y)
"""
function grad_vandermonde2D(p, x, y)
    dV = Array{Float64,5}(undef, length(x),length(y), p+1,p+1, 2)
    lx = Basis1D.vandermonde1D(p, x)
    ly = Basis1D.vandermonde1D(p, y)
    dlx = Basis1D.grad_vandermonde1D(p, x)
    dly = Basis1D.grad_vandermonde1D(p, y)
    # Derivative wrt x
    for ipy = 1:p+1
        for ipx = 1:p+1
            dV[:,:,ipx,ipy,1] = dlx[:,ipx]*(ly[:,ipy]')
        end
    end
    # Derivative wrt y
    for ipy = 1:p+1
        for ipx = 1:p+1
            dV[:,:,ipx,ipy,2] = lx[:,ipx]*(dly[:,ipy]')
        end
    end
    return dV
end

"""
    gauss_quad2D(p)
Return the 2D Gaussian quadrature points and weights (x,w) of order p on the domain [-1,1]^2
"""
function gauss_quad2D(p)
    (x1D,w1D) = Basis1D.gauss_quad(p)
    n = length(x1D)
    x2D = Array{Float64,3}(undef, 2, n,n)
    w2D = Array{Float64,2}(undef, n,n)
    for iy = 1:n
        for ix = 1:n
            x2D[1,ix,iy] = x1D[ix]
            x2D[2,ix,iy] = x1D[iy]
            w2D[ix,iy] = w1D[ix]*w1D[iy]
        end
    end
    return (x2D, w2D)
end

"""
    chebyshev2D(p)
Return the 2D chebyshev points of order p on the domain [-1,1]^2
"""
function chebyshev2D(p)
    x1D = Basis1D.chebyshev(p)
    n = length(x1D)
    x2D = Array{Float64,3}(undef, 2, n, n)
    for iy = 1:n
        for ix = 1:n
            x2D[1,ix,iy] = x1D[ix]
            x2D[2,ix,iy] = x1D[iy]
        end
    end
    return x2D
end

"""
    equidistant_nodes2D(p)
Return the "order p" 2D points on the domain [-1,1]^2
"""
function equidistant_nodes2D(p)
    x1D = -1:2/p:1
    n = length(x1D)
    x2D = Array{Float64,3}(undef, 2, n, n)
    for iy = 1:n
        for ix = 1:n
            x2D[1,ix,iy] = x1D[ix]
            x2D[2,ix,iy] = x1D[iy]
        end
    end
    return x2D
end

"""
    interpolation_matrix2D(xy_from, xy_to)
Compute an interpolation matrix from a set of 2D points to another set of 2D points.
Assumes that points interpolating from provide enough accuracy (aka - they are
well spaced out and of high enough order), and define a square. Points
interpolating onto can be of any size, but must be defined on this same square.
Interpolation matrix ∈ ℜ^(2D size of xy_to x 2D size of xy_from)
"""
function interpolation_matrix2D(xy_from, xy_to)
    # Create nodal representation of reference bases
    (x_from, y_from) = xy_from
    order = size(x_from,1) - 1 # Assumes order = (size of x_from) - 1
    @assert size(x_from,1) == size(y_from,1)
    n_from = size(x_from,1)*size(y_from,1)
    l_from = vandermonde2D(order, x_from, y_from)

    eye = diagm(0=>ones(n_from))
    V = reshape(l_from, n_from,n_from)
    coeffs_phi = V \ eye

    # Compute reference bases on the output points
    l_to = vandermonde2D(order, xy_to...)
    n_to = prod(size.(xy_to,1))
    V_to = reshape(l_to, n_to,n_from)

    # Construct interpolation matrix
    Interp2D = V_to*coeffs_phi
end

"""
    dPhi2D(xy_from, xy_to)
Compute the gradient of the basis functions defined by 2D points on another
set of 2D points.

Assumes that points interpolating from provide enough accuracy (aka - they are
well spaced out and of high enough order), and define a square. Points
interpolating onto can be of any size, but must be defined on this same square.
dPhi ∈ ℜ^(size of xy_to × size of xy_from × 2)
"""
function dPhi2D(xy_from, xy_to)
    dim = 2
    # Create nodal representation of reference bases
    (x_from, y_from) = xy_from
    order = size(x_from,1) - 1 # Assumes order = (size of x_from) - 1
    @assert size(x_from,1) == size(y_from,1)
    n_from = size(x_from,1)*size(y_from,1)
    l_from = vandermonde2D(order, x_from, y_from)

    eye = diagm(0=>ones(n_from))
    V = reshape(l_from, n_from,n_from)
    coeffs_phi = V \ eye

    # Compute derivative of reference bases on the output points
    dl_to = grad_vandermonde2D(order, xy_to...)
    n_to = prod(size.(xy_to,1))
    dV = reshape(dl_to, n_to,n_from,dim)
    # Construct gradient of phi = dV*coeffs_phi
    dPhi_to = Array{Float64,3}(undef, n_to,n_from,dim)
    for l = 1:dim
        dPhi_to[:,:,l] = dV[:,:,l]*coeffs_phi
    end
    return dPhi_to
end

end

max_stars_repo_name: NoseKnowsAll/DGToolkit

max_issues_repo_name: NoseKnowsAll/DGToolkit

max_forks_repo_name: NoseKnowsAll/DGToolkit


text:
struct GBMDPSolver <: Solver
    mdp_solver::Solver
    updater::Union{Function, Updater}
end

function solve(solver::GBMDPSolver, pomdp::POMDP)
    if isa(solver.updater, Function)
        updater = solver.updater(pomdp)
    else
        updater = solver.updater
    end
    # belief_mdp = GenerativeBeliefMDP(deepcopy(pomdp), updater)
    belief_mdp = MeanRewardBeliefMDP(deepcopy(pomdp),
                                     updater.resample,
                                     updater.max_frac_replaced
                                    )
    return solve(solver.mdp_solver, belief_mdp)
end

function solve(qs::QMDPSolver, bmdp::Union{GenerativeBeliefMDP, MeanRewardBeliefMDP})
    return solve(qs, bmdp.pomdp)
end

ext: jl


text:
module Trapani

using ArgCheck

using WignerFunctions.Risbo: initial_wigner, Cache

function delta_recurrence(l::Integer, m₁::Integer, m₂::Integer, cache::Cache)
    if l < 0 || abs(l) < m₁ || abs(l) < m₂
        return 0
    elseif m₁ < 0 && m₂ < 0
        return delta_recurrence(l, -m₂, -m₁, cache)
    elseif m₁ < 0
        return ((l - m₂) % 2 == 0 ? 1 : -1) * delta_recurrence(l, -m₁, m₂, cache)
    elseif m₂ < 0
        return ((l - m₁) % 2 == 0 ? 1 : -1) * delta_recurrence(l, m₁, -m₂, cache)
    elseif m₂ < m₁
        return ((m₁ - m₂) % 2 == 0 ? 1 : -1) * delta_recurrence(l, m₂, m₁, cache)
    elseif (l=l, m₁=m₁, m₂=m₂) ∈ keys(cache)
        return cache[(l=l, m₁=m₁, m₂=m₂)]
    end
    result = if l == 0 && m₁ == 0 && m₂ == 0
            1
        elseif l == m₁
            (
                sqrt((2l*l - l) / (2 * (l + m₂) * (l + m₂ - 1)))
                * delta_recurrence(l - 1, l - 1, m₂ - 1, cache)
            )
        else
            (
                2m₂ / sqrt((l - m₁) * (l + m₁ + 1))
                * delta_recurrence(l, m₁ + 1, m₂, cache)
                - sqrt((l - m₁ - 1) * (l + m₁ + 2) / ((l - m₁) * (l + m₁ + 1)))
                * delta_recurrence(l, m₁ + 2, m₂, cache)
            )
    end
    cache[(l=l, m₁=m₁, m₂=m₂)] = result
    result
end

function delta_recurrence(l::Integer, m₁::Integer, m₂::Integer)
    delta_recurrence(l, m₁, m₂, Cache{Float64}())
end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Yosys"
version = v"0.11.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/YosysHQ/yosys.git", "360fed8e4d611fa725a4526cf960383b8e6c6e64")
]

dependencies = [
    Dependency("boost_jll"; compat="=1.76.0"), # max gcc7
    Dependency("Readline_jll"; compat="8.1.1"),
    Dependency("Tcl_jll"; compat="8.6.11"),
    Dependency("Zlib_jll"; compat="1.2.11"),
    Dependency("Libffi_jll"; compat="~3.2.2")
]

# Bash recipe for building across all platforms
script = raw"""
cd yosys
if [[ "${target}" == *-apple-* ]] || [[ "${target}" == *-freebsd* ]]; then
    CONFIG=clang
    OS=Darwin
elif [[ "${target}" == *-x86_64-w64-mingw32* ]]; then
    CONFIG=msys2-64
    OS=Windows
elif [[ "${target}" == *-i686-w64-mingw32* ]]; then
    CONFIG=msys2-32
    OS=Windows
else
    CONFIG=gcc
    OS=Linux
fi
# TODO: we need mingw-dlfcn wrappers for windows plugin support
# TODO: ABC does not build on windows
make ENABLE_LIBYOSYS=1 OS=${OS} CONFIG=${CONFIG} PREFIX=${prefix} LIBDIR=${libdir} -j${nproc}
make install ENABLE_LIBYOSYS=1 OS=${OS} PREFIX=${prefix} LIBDIR=${libdir}
if [[ "${target}" == *-apple-* ]]; then
    mv ${prefix}/lib/libyosys.so  ${prefix}/lib/libyosys.dylib
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> arch(p) != "armv6l" && (Sys.isapple(p) || Sys.islinux(p)), supported_platforms(;experimental=true))
platforms = expand_cxxstring_abis(platforms)
# For some reasons, building for CXX03 string ABI doesn't actually work, skip it
filter!(x -> cxxstring_abi(x) != "cxx03", platforms)

# The products that we will ensure are always built
products = Product[
    ExecutableProduct("yosys", :yosys),
    ExecutableProduct("yosys-config", :yosys_config),
    ExecutableProduct("yosys-filterlib", :yosys_filterlib),
    ExecutableProduct("yosys-smtbmc", :yosys_smtbmc),
    ExecutableProduct("yosys-abc", :yosys_abc),
    LibraryProduct("libyosys", :libyosys)
]

# Build the tarballs, and possibly a `build.jl` as well.
# gcc7 constraint from boost
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"7")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    wclogptype, wclogphydrogentype,
    wclogpcontrib, wclogp

const WCLOGP_TABLE = YAML.load(open(
    joinpath(dirname(@__FILE__), "..", "assets", "const", "wclogp.yaml")
))["logP"]

const ALIPH_HETERO = Set([:H, :C, :N, :O, :P, :S, :F, :Cl, :Br, :I])
const AROM_HETERO = Dict(
    :H  => :C18, :C => :C21, :N => :C22, :O => :C23, :S => :C24,
    :F => :C14, :Cl => :C15, :Br => :C16, :I => :C17
)
# TODO:
const P_BLOCK = (:Al, :B, :Si, :Ga, :Ge, :As, :Se, :Sn, :Te, :Pb,
                 :Ne, :Ar, :Kr, :Xe, :Rn)
const D_BLOCK = (:Fe, :Co, :Cu, :Zn, :Tc, :Cd, :Pt, :Au, :Hg, :Gd)

"""
    wclogptype(mol::GraphMol)

Return Wildman-Crippen LogP atom types.
"""
@cache function wclogptype(mol::GraphMol)
    atomtypes = fill(:undef, nodecount(mol))
    atomsymbol_ = atomsymbol(mol)
    charge_ = charge(mol)
    bondorder_ = bondorder(mol)
    heavyatomcount_ = heavyatomcount(mol)
    hcount_ = hcount(mol)
    pielectron_ = pielectron(mol)
    hybridization_ = hybridization(mol)
    isaromatic_ = isaromatic(mol)
    isaromaticbond_ = isaromaticbond(mol)
    # Carbons
    for i in findall(atomsymbol_ .=== :C)
        if hybridization_[i] === :sp3
            # Aliphatic sp3 hybrid (C1-4,8-12,27)
            nbrs = collect(adjacencies(mol, i))
            if !isempty(setdiff(atomsymbol_[nbrs], ALIPH_HETERO))
                atomtypes[i] = :C27  # Adjacent to inorganic atoms
            elseif all(isaromatic_[nbrs] .=== false)
                if all(atomsymbol_[nbrs] .=== :C)
                    # Aliphatic carbon (C1,2)
                    atomtypes[i] = heavyatomcount_[i] < 3 ? :C1 : :C2
                else
                    # Adjacent to heteroatoms (C3,4)
                    atomtypes[i] = heavyatomcount_[i] < 3 ? :C3 : :C4
                end
            else
                # Adjacent to aromatic atoms (C8-12)
                if heavyatomcount_[i] == 1
                    atomtypes[i] = atomsymbol_[nbrs[1]] === :C ? :C8 : :C9
                elseif heavyatomcount_[i] == 2
                    atomtypes[i] = :C10
                elseif heavyatomcount_[i] == 3
                    atomtypes[i] = :C11
                elseif heavyatomcount_[i] == 4
                    atomtypes[i] = :C12
                end
            end
        elseif isaromatic_[i]
            # Aromatic (C13-25)
            if hcount_[i] == 1
                atomtypes[i] = :C18  # Aromatic non-substituted
                continue
            end
            subst = -1
            substbond = -1
            aromcnt = 0
            for (inc, adj) in neighbors(mol, i)
                if isaromaticbond_[inc]
                    aromcnt += 1
                else
                    subst = adj
                    substbond = inc
                end
            end
            if aromcnt == 3
                atomtypes[i] = :C19 # Bridgehead
            elseif !haskey(AROM_HETERO, atomsymbol_[subst])
                atomtypes[i] = :C13 # Inorganic substituent
            elseif isaromatic_[subst]
                atomtypes[i] = :C20 # Aromatic substituent
            elseif bondorder_[substbond] == 2
                atomtypes[i] = :C25 # Double bond substituent
            else
                # Typical substituent (C14-17,21-24)
                atomtypes[i] = AROM_HETERO[atomsymbol_[subst]]
            end
        elseif hybridization_[i] === :sp2
            # Aliphatic sp2 hybrid (C5-6,26)
            arom = 0
            het = 0
            for (inc, adj) in neighbors(mol, i)
                if bondorder_[inc] == 2 && atomsymbol_[adj] !== :C
                    het += 1
                elseif isaromatic_[adj]
                    arom += 1
                end
            end
            if het > 0
                atomtypes[i] = :C5 # Double bond to non-C
            elseif arom > 0
                atomtypes[i] = :C26 # Adjacent to aromatic
            else
                atomtypes[i] = :C6 # Double bond to C
            end
        elseif hybridization_[i] === :sp
            # Aliphatic sp hybrid (C6-7)
            bonds = collect(incidences(mol, i))
            if any(bondorder_[bonds] .== 3)
                atomtypes[i] = :C7 # Alkyne, Nitrile
            else
                atomtypes[i] = :C6 # Allene
            end
        end
        if atomtypes[i] === :undef
            atomtypes[i] = :CS # Not found
        end
    end

    # Nitrogens
    for i in findall(atomsymbol_ .=== :N)
        if hybridization_[i] === :sp3
            if charge_[i] > 0
                # Ammonium (N10, N13)
                if heavyatomcount_[i] == 4
                    atomtypes[i] = :N13  # Quart-ammonium
                elseif heavyatomcount_[i] < 4
                    atomtypes[i] = :N10  # Protonated amine
                end
            elseif charge_[i] == 0
                # Aliphatic amine (N1,2,7)
                if heavyatomcount_[i] == 1
                    atomtypes[i] = :N1
                elseif heavyatomcount_[i] == 2
                    atomtypes[i] = :N2
                elseif heavyatomcount_[i] == 3
                    atomtypes[i] = :N7
                end
            end
        elseif isaromatic_[i]
            # Aromatic amine (N11,12)
            if charge_[i] > 0
                atomtypes[i] = :N12 # Protonated aromatic
            elseif charge_[i] == 0
                atomtypes[i] = :N11 # Unprotonated aromatic
            end
        elseif hybridization_[i] === :sp2
            if pielectron_[i] == 2
                adjs = collect(adjacencies(mol, i))
                if any(isaromatic_[adjs])
                    # Amine adjacent to aromatic (N3,4,8)
                    if heavyatomcount_[i] == 1
                        atomtypes[i] = :N3
                    elseif heavyatomcount_[i] == 2
                        atomtypes[i] = :N4
                    elseif heavyatomcount_[i] == 3
                        atomtypes[i] = :N8
                    end
                else
                    # Aliphatic amine (N1,2,7)
                    if heavyatomcount_[i] == 1
                        atomtypes[i] = :N1
                    elseif heavyatomcount_[i] == 2
                        atomtypes[i] = :N2
                    elseif heavyatomcount_[i] == 3
                        atomtypes[i] = :N7
                    end
                end
            elseif pielectron_[i] == 1
                # Imine (N5,6)
                if charge_[i] > 0
                    atomtypes[i] = :N13  # Quart-ammonium
                elseif charge_[i] == 0
                    if heavyatomcount_[i] == 1
                        atomtypes[i] = :N5
                    elseif heavyatomcount_[i] == 2
                        atomtypes[i] = :N6
                    end
                end
            end
        elseif hybridization_[i] === :sp
            # sp nitrogen (N9,13,14)
            adjs = adjacencies(mol, i)
            adjatoms = [atomsymbol_[adj] for adj in adjs]
            if charge_[i] == 1 && issetequal(adjatoms, [:N, :N])
                atomtypes[i] = :N14  # Other ionized nitrogen (Azide)
            elseif charge_[i] == 1 && issetequal(adjatoms, [:C, :N])
                atomtypes[i] = :N13  # Quart-ammonium (Diazo)
            elseif charge_[i] == 0 && :C in adjatoms
                atomtypes[i] = :N9 # Nitrile
            end
        end
        if atomtypes[i] === :undef
            if charge_[i] < 0
                atomtypes[i] = :N14  # Other ionized nitrogen (Nitride)
            else
                atomtypes[i] = :NS  # Others
            end
        end
    end

    # Oxygens
    for i in findall(atomsymbol_ .=== :O)
        if hybridization_[i] === :sp3
            if hcount_[i] > 0
                atomtypes[i] = :O2  # Alcohol (O2)
            elseif charge_[i] == 0 && heavyatomcount_[i] == 2
                atomtypes[i] = :O3  # Aliphatic ether
            elseif charge_[i] < 0 && heavyatomcount_[i] == 1
                # Oxide (O5,6,7,12)
                adj = pop!(adjacencies(mol, i))
                if atomsymbol_[adj] in (:O, :N)
                    atomtypes[i] = :O5 # O2 or N-oxide
                elseif atomsymbol_[adj] == :S
                    atomtypes[i] = :O6 # S-oxide
                elseif atomsymbol_[adj] == :C
                    cadjs = adjacencies(mol, adj)
                    cadjatoms = [atomsymbol_[cadj] for cadj in cadjs]
                    if sum(cadjatoms .=== :O) == 2
                        atomtypes[i] = :O12 # Acid (O12)
                    end
                end
                if atomtypes[i] === :undef
                    atomtypes[i] = :O7 # Other oxide
                end
            end
        elseif isaromatic_[i]
            atomtypes[i] = :O1  # Aromatic oxygen (O1)
        elseif hybridization_[i] === :sp2 && charge_[i] == 0
            if heavyatomcount_[i] == 2
                adjs = collect(adjacencies(mol, i))
                if any(isaromatic_[adjs])
                    atomtypes[i] = :O4  # Aromatic ether
                else
                    atomtypes[i] = :O3  # Aliphatic ether
                end
            elseif pielectron_[i] == 2
                atomtypes[i] = :O2  # Alcohol
            else
                # Carbonyl
                adj = pop!(adjacencies(mol, i))
                cadjs = adjacencies(mol, adj)
                pop!(cadjs, i)
                cadjs = collect(cadjs)
                if isaromatic_[adj]
                    atomtypes[i] = :O8  # Aromatic carbonyl
                elseif all(atomsymbol_[cadjs] .!== :C)
                    atomtypes[i] = :O11  # Carbonyl heteroatom
                elseif any(isaromatic_[cadjs])
                    atomtypes[i] = :O10  # Carbonyl aromatic
                else
                    atomtypes[i] = :O9  # Carbonyl aliphatic
                end
            end
        end
        if atomtypes[i] === :undef
            atomtypes[i] = :OS  # Others
        end
    end

    # Others
    for i in findall(atomtypes .=== :undef)
        if atomsymbol_[i] in (:F, :Cl, :Br, :I)
            if charge_[i] == 0
                atomtypes[i] = atomsymbol_[i]
            else
                atomtypes[i] = :Hal  # Ionic halogens
            end
        elseif atomsymbol_[i] === :P
            atomtypes[i] = :P
        elseif atomsymbol_[i] === :S
            if isaromatic_[i]
                atomtypes[i] = :S3  # Aromatic sulfur
            elseif charge_[i] === 0
                atomtypes[i] = :S1  # Aliphatic sulfur
            else
                atomtypes[i] = :S2  # Ionic sulfur
            end
        elseif atomsymbol_[i] === :H
            adjs = adjacencies(mol, i)
            if isempty(adjs)
                atomtypes[i] = :HS  # Proton or hydride (But not meaningful)
            else
                adj = pop!(adjs)
                atomtypes[i] = wclogphydrogentype(mol, adj)
            end
        elseif atomsymbol_[i] in P_BLOCK
            atomtypes[i] = :Me1
        elseif atomsymbol_[i] in D_BLOCK
            atomtypes[i] = :Me2
        end
    end
    return atomtypes
end

function wclogphydrogentype(mol::GraphMol, i)
    atomsymbol_ = atomsymbol(mol)
    hybridization_ = hybridization(mol)
    isaromatic_ = isaromatic(mol)
    if atomsymbol_[i] in (:H, :C)
        return :H1  # Hydrocarbon
    elseif atomsymbol_[i] === :N
        return :H3  # Amine
    elseif atomsymbol_[i] === :O
        adjs = adjacencies(mol, i)
        heavy = [adj for adj in adjs if atomsymbol_[adj] !== :H]
        isempty(heavy) && return :H2  # Alcohol (H2O)
        a = pop!(heavy)
        if atomsymbol_[a] === :N
            return :H3  # Hydroxyamine
        elseif atomsymbol_[a] in (:O, :S)
            return :H4  # Peroxide, sulfoxide
        elseif atomsymbol_[a] === :C
            if hybridization_[a] === :sp2 && !isaromatic_[a]
                return :H4  # Acid
            else
                return :H2  # Alcohol
            end
        else
            return :H2  # Alcohol (Other hydroxyl)
        end
    else
        return :H2  # Other hydrate
    end
end

function wclogpcontrib(mol::GraphMol)
    contrib = zeros(nodecount(mol))
    atomsymbol_ = atomsymbol(mol)
    hcount_ = hcount(mol)
    atomtypes = wclogptype(mol)
    for i in 1:nodecount(mol)
        cont = get(WCLOGP_TABLE, string(atomtypes[i]), 0)
        if atomsymbol_[i] == :H
            contrib[i] = 0  # avoid double count
        elseif hcount_[i] == 0
            contrib[i] = cont
        else
            htype = wclogphydrogentype(mol, i)
            hcont = WCLOGP_TABLE[string(htype)] * hcount_[i]
            contrib[i] = cont + hcont
        end
    end
    return contrib
end

"""
    wclogp(mol::GraphMol) -> Float64

Return predicted logP value calculated by using Wildman and Crippen method.

# Reference

1. Wildman, S. A. and Crippen, G. M. (1999). Prediction of Physicochemical
   Parameters by Atomic Contributions. Journal of Chemical Information and
   Modeling, 39(5), 868–873. https://doi.org/10.1021/ci990307l
"""
wclogp(mol::GraphMol; digits=2
    ) = round(reduce(+, wclogpcontrib(mol); init=0), digits=digits)


text:
using WikitextParser
using Test

@testset "WikitextParser.jl" begin
    # Write your own tests here.
    include("test-wiktionary.jl")
end

@testset "parentheses" begin
    wt = wikitext(parentheses=true)
    @test parse(wt, "(balanced parenthesis)") ==
        aline(TokenPair(:paren,AbstractToken[l"balanced",delim" ", l"parenthesis"]))

    @test parse(wt, "(unclosed parenthesis") ==
        aline(Token(:paren,"("),l"unclosed",delim" ", l"parenthesis")

    @test parse(wt, "unopened parenthesis)") ==
        aline(l"unopened",delim" ", l"parenthesis",Token(:paren,")"))

    @test parse(wt, ")") == aline(Token(:paren,")"))

    @test parse(wt, "{{t|(balanced parenthesis)}}") ==
        aline(Template("t",[""=>aline(TokenPair(
            :paren,AbstractToken[l"balanced", delim" ", l"parenthesis"]))]))
end

@testset "wikitext links" begin
    @test parse(wiki_link(AnyChar()), "[[Test|1]]") ==
        WikiLink("wikt:de", "Test", "", "1")

    @test parse(wiki_link(), "[[Test]]") ==
        WikiLink("wikt:de", "Test", "", "")

    @test parse(wiki_link(), "[[Texas|Lone Star State]]") ==
        WikiLink("wikt:de", "Texas", "", "Lone Star State")
    # @test parse(wiki_link, "[[Texas#Racism|Lone Star State]]") ==
    # (target="Texas",section="Racism",title="Lone Star State")

    ## label should be parsed for styling info
    @test parse(wiki_link(), "[[File:wiki.png|thumb|Wikipedia logo]]") ==
        WikiLink("File", "wiki.png", "", "thumb|Wikipedia logo")

    @test parse(wiki_external_link,
                "[http://www.bictxt.de/Schule]") ==
                    TokenPair(:link,[T"http"protocol, T"://"delimiter, T"www.bictxt.de"domain,
                                     T"/"delimiter,T"Schule"capitalized])
    ## todo: bare urls

    @test parse(wikitext(parentheses=true), "([[Reich]]e)") ==
        aline(TokenPair(:paren,AbstractToken[WikiLink("wikt:de","Reich","","Reiche")]))
end

import WikitextParser: expand_numbers
@testset "number references" begin
    @test parse(expand_numbers,"[1]") == ["1"]
    @test parse(expand_numbers,"[1-3]") == ["1","2","3"]
    @test parse(expand_numbers,"[1-3,10]") == ["1","2","3","10"]
end

@testset "html" begin
    @test parse(WikitextParser.attributes, "class=\"a\" class=\"b\"") == [Token(:class,"a"), Token(:class,"b")]
    @test parse(WikitextParser.attributes, "size=10% class=\"b\"") == [Token(:size,"10%"), Token(:class,"b")]
    @test parse(WikitextParser.attributes, "class=simple") == [Token(:class,"simple")]
    @test parse(wikitext(), "a") == aline(Node(:font,[],aline(l"a")))
    wtp = wikitext();
    ## TODO: profile and optimize
    x = wtp("<ref name=\"Georges\">{{Ref-Georges}}</ref>")
    @test x == [ Line(Token[],AbstractToken[
        Node(:ref, Token[ Token(:name,"Georges") ],
             [ Line(Token[], AbstractToken[Template{Token,AbstractToken}(
                 "Ref-Georges",[])] ) ] ) ] ) ]
    @test string(x)=="<ref name=\"Georges\">{{Ref-Georges}}</ref>"
end

@testset "wikitext templates" begin
    @test parse(wikitext(),"{{}}")==[Line(NamedString[], AbstractToken[Template("")])]
    @test parse(wikitext(),"|")==[Line(NamedString[], AbstractToken[Token(:delimiter,"|")])]
    @test parse(wikitext(), "{{test}}")[1] == Line(LineContent[Template("test")])
    parse(wikitext(), "{{Ü-Tabelle|Ü-links|Ü-rechts}}") ==
        aline(Template("Ü-Tabelle",[""=>aline(l"Ü",delim"-",l"links"),
                                    ""=>aline(l"Ü",delim"-",l"rechts")]))
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "OpenFAST"
version = v"3.0.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/OpenFAST/openfast/archive/refs/tags/v$(version).tar.gz", "9af57af054e4128b6e257a76da368dc4ad0c7fbb2b22d51fc7ea63cdf999c530"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/openfast-*/

if [[ "${target}" == *-mingw* ]]; then
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-lowercase-windows-include.patch
fi

mkdir build && cd build

cmake .. \
-DCMAKE_INSTALL_PREFIX=${prefix} \
-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
-DCMAKE_BUILD_TYPE=Debug \
-DBUILD_SHARED_LIBS=ON \
-DBUILD_FASTFARM=OFF \
-DORCA_DLL_LOAD=OFF \
-DOPENMP=ON \
-DBLAS_LIBRARIES="${libdir}/libopenblas.${dlext}" \
-DLAPACK_LIBRARIES="${libdir}/libopenblas.${dlext}"

#WARNING: compiling this locally can go crazy and lock up your machine, using only 2 jobs and Debug version to make it behave
make -j2
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line

platforms = expand_gfortran_versions(supported_platforms())

#remove aarch64-linux-musl from platforms, this platform does not currently have IEEE_ARITHMETIC enabled for gfortran under the current configure set up
#see also: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=100662
filter!(p -> !(arch(p) == "aarch64" && Sys.islinux(p) && libc(p) == "musl"), platforms)

#remove aarch-apple-darwin platforms, same issue as aarch64-linux-musl
filter!(p -> !(Sys.isapple(p) && arch(p) == "aarch64"), platforms)

#filter arm platforms
filter!(p -> arch(p) ∉ ("armv6l", "armv7l"), platforms)

#filter windows platforms - MinGW is not well supported relative to MSVC
filter!(p -> !Sys.iswindows(p), platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("hydrodyn_driver", :hydrodyn_driver),
    LibraryProduct("libaerodynlib", :libaerodynlib),
    LibraryProduct("libelastodynlib", :libelastodynlib),
    LibraryProduct("libopenfastlib", :libopenfastlib),
    LibraryProduct("libmaplib", :libmaplib),
    LibraryProduct("libbeamdynlib", :libbeamdynlib),
    LibraryProduct("libnwtclibs", :libnwtclibs),
    ExecutableProduct("inflowwind_driver", :inflowwind_driver),
    LibraryProduct("libversioninfolib", :libversioninfolib),
    LibraryProduct("libuaaerolib", :libuaaerolib),
    ExecutableProduct("orca_driver", :orca_driver),
    LibraryProduct("libmoordynlib", :libmoordynlib),
    LibraryProduct("libmapcpplib", :libmapcpplib),
    LibraryProduct("libwdlib", :libwdlib),
    ExecutableProduct("feam_driver", :feam_driver),
    ExecutableProduct("unsteadyaero_driver", :unsteadyaero_driver),
    LibraryProduct("libextptfm_mckflib", :libextptfm_mckflib),
    LibraryProduct("libscdataexlib", :libscdataexlib),
    ExecutableProduct("turbsim", :turbsim),
    LibraryProduct("libfvwlib", :libfvwlib),
    LibraryProduct("libifwlib", :libifwlib),
    ExecutableProduct("openfast", :openfast),
    LibraryProduct("libafinfolib", :libafinfolib),
    LibraryProduct("liborcaflexlib", :liborcaflexlib),
    ExecutableProduct("aerodyn_driver", :aerodyn_driver),
    LibraryProduct("libfoamfastlib", :libfoamfastlib),
    LibraryProduct("libfeamlib", :libfeamlib),
    LibraryProduct("libsctypeslib", :libsctypeslib),
    LibraryProduct("libsubdynlib", :libsubdynlib),
    ExecutableProduct("subdyn_driver", :subdyn_driver),
    LibraryProduct("libopenfoamtypeslib", :libopenfoamtypeslib),
    LibraryProduct("libaeroacoustics", :libaeroacoustics),
    ExecutableProduct("dwm_driver_wind_farm", :dwm_driver_wind_farm),
    LibraryProduct("libscdataextypeslib", :libscdataexttypeslib),
    ExecutableProduct("moordyn_driver", :moordyn_driver),
    LibraryProduct("libhydrodynlib", :libhydrodynlib),
    LibraryProduct("libopenfast_postlib", :libopenfast_postlib),
    LibraryProduct("libscfastlib", :libscfastlib),
    LibraryProduct("libopenfast_prelib", :libopenfast_prelib),
    LibraryProduct("libaerodyn14lib", :libaerodyn14lib),
    ExecutableProduct("beamdyn_driver", :beamdyn_driver),
    ExecutableProduct("servodyn_driver", :servodyn_driver),
    LibraryProduct("libawaelib", :libawaelib),
    LibraryProduct("libicedynlib", :libicedynlib),
    LibraryProduct("libicefloelib", :libicefloelib),
    LibraryProduct("libservodynlib", :libservodynlib)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("OpenBLAS32_jll")
    #Dependency(PackageSpec(name="OpenBLAS_jll", uuid="4536629a-c528-5b80-bd46-f80d51c5b363"))
    Dependency(PackageSpec(name="LAPACK_jll", uuid="51474c39-65e3-53ba-86ba-03b1b862ec14"))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"6")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
# using Julog

@testset "Extended Built-ins" begin

@test resolve(@julog(atom(X)), Clause[])[1] == true
@test resolve(@julog(nonvar(X)), Clause[])[1] == false
@test resolve(@julog(ground(X)), Clause[])[1] == false
@test resolve(@julog(ground(atomic(X))), Clause[])[1] == false
@test resolve(@julog(ground(atomic(x))), Clause[])[1] == true
@test resolve(@julog(arg(1, f(a, b), a)), Clause[])[1] == true
@test resolve(@julog(functor(f(a, b), f, 2)), Clause[])[1] == true
@test resolve(@julog(functor(f(a, b), F, A)), Clause[])[1] == true
@test resolve(@julog(univ(Term, [baz, foo(1)])), Clause[])[1] == true
@test resolve(@julog(univ(foo(hello, X), List)), Clause[])[1] == true

end

# goal = @julog arg(1, f(a,b), a)
# res = resolve(goal, Clause[])


text:
module PositionVelocityTime
    using CoordinateTransformations, DocStringExtensions, Geodesy, GNSSDecoder, GNSSSignals, LinearAlgebra, Parameters
    using Unitful: s, Hz

    export  calc_PVT,
            is_sat_healthy_and_decodable,
            PVTSolution,
            DOP,
            sat_position_rotation,
            sat_position_ECI_2_ECEF,
            sat_position_ECEF,
            user_position,
            SatelliteState,
            get_num_used_sats,
            get_gdop,
            get_pdop,
            get_hdop,
            get_vdop,
            get_tdop

    """
    Struct of decoder, code- and carrierphase of satellite vehicle
    """
    @with_kw struct SatelliteState{CP <: Real}
        decoder_state::GNSSDecoderState
        code_phase::CP
        carrier_phase::CP = 0
    end

    """
    DOP, stores the Dilution of Precision values
    """
    @with_kw struct DOP
        GDOP::Union{Nothing, Float64} = nothing
        PDOP::Union{Nothing, Float64} = nothing
        VDOP::Union{Nothing, Float64} = nothing
        HDOP::Union{Nothing, Float64} = nothing
        TDOP::Union{Nothing, Float64} = nothing
    end

    """
    PVTSolution, with used sats and el and az.
    """
    @with_kw struct PVTSolution
        pos::ECEF = ECEF([0,0,0])
        receiver_time_correction::Float64 = 0
        DOP_val::DOP = DOP()
        #following is added for satllite observation
        used_sats::AbstractVector{Int64} = []
        elevation_sats::AbstractVector{Float64} = []
        azimuth_sats::AbstractVector{Float64} = []
    end

    function PVTSolution(pvt_old::PVTSolution, used_sats, elevation_sats,azimuth_sats)
        PVTSolution(pvt_old.pos, pvt_old.receiver_time_correction, pvt_old.DOP_val, used_sats, elevation_sats, azimuth_sats)
    end

    function get_num_used_sats(pvt_solution::PVTSolution)
        length(pvt_solution.used_sats)
    end

    #Get methods for single DOP values
    function get_gdop(pvt_sol::PVTSolution)
        return pvt_sol.DOP_val.GDOP
    end

    function get_pdop(pvt_sol::PVTSolution)
        return pvt_sol.DOP_val.PDOP
    end

    function get_vdop(pvt_sol::PVTSolution)
        return pvt_sol.DOP_val.VDOP
    end

    function get_hdop(pvt_sol::PVTSolution)
        return pvt_sol.DOP_val.HDOP
    end

    function get_tdop(pvt_sol::PVTSolution)
        return pvt_sol.DOP_val.TDOP
    end

    """
    Calculates elevation and azimuth of satellite

    $SIGNATURES
    user_pos_ecef: user position in ecef coordinates
    sat_pos_ecef: satellite position in ecef coordinates

    This function calculates the elevation and azimuth of the satellite in degree
    """

    function get_elevation_azimuth(user_pos_ecef::ECEF, sat_pos_ecef::ECEF)
        sat_enu = ENUfromECEF(user_pos_ecef, wgs84)(sat_pos_ecef)
        sat_enu_sph = SphericalFromCartesian()(sat_enu)
        az = rad2deg(sat_enu_sph.θ)
        el = rad2deg(sat_enu_sph.ϕ)
        return el, az
    end

    """
    Calculates ECEF position of user

    $SIGNATURES
    ´sat_state´: satellite state, combining decoded data, code- and carrierphase

    ´min_elevation´: The minimum elevation for a satellite, that can be used for a positioning
                     Hint: The elevation starts here with 0° at the horizon and has its maximum with 90° in the zenith.
                     min_elevation = -100, deactivates the  elevation filter

    This function calculates the position of the user in ECEF coordinates
        incl. earth rotation correction

    The implementation follows IS-GPS-200K Table 20-IV.
    """
    function calc_PVT(
        satellite_states::AbstractVector{SatelliteState{Float64}},
        prev_pvt::PVTSolution = PVTSolution(),
        min_elevation::Int64 = -100 #elevation filter is normally turned off
        )

        c = satellite_states[1].decoder_state.constants.c
        usable_satellite_states = filter(x -> is_sat_healthy_and_decodable(x.decoder_state), satellite_states)
        length(usable_satellite_states) >= 4 || throw(SmallData("Not enough usable SV Data"))
        sv_positions = map(x -> sat_position_ECEF(x), usable_satellite_states)

        if get_num_used_sats(prev_pvt) == 0
            pseudoranges = pseudo_ranges(usable_satellite_states)
            pvt_sol = user_position([0,0,0], sv_positions, pseudoranges)
            prev_pos = pvt_sol.pos
        else
            prev_pos = prev_pvt.pos
        end

        # Remove Satellite with low elevation angle
        if min_elevation >= 0 && min_elevation < 90 # min elevaiton filter
            out = map( x-> get_elevation_azimuth(prev_pos, ECEF(x[1],x[2],x[3])), sv_positions)
            el, az = map(x->getindex.(out,x), 1:2)
            el_filter = el .>= min_elevation
            usable_satellite_states = usable_satellite_states[el_filter]
            el_save = el[el_filter]
            az_save = az[el_filter]
            PRN_save = map(x -> x.decoder_state.PRN, usable_satellite_states)
        end

        length(usable_satellite_states) >= 4 || throw(SmallData("Not enough usable SV Data, with elevation filter."))

        #compensation of earth rotation within TOF
        dt =  map(x -> norm(x-prev_pos)/c, sv_positions)
        sv_positions = map(i -> sat_position_rotation(usable_satellite_states[i], dt[i]), 1:length(usable_satellite_states))

        pseudoranges = pseudo_ranges(usable_satellite_states)
        pvt = user_position(prev_pos, sv_positions, pseudoranges)
        pvt = PVTSolution(pvt, PRN_save, el_save, az_save)
        return pvt
    end

    include("pseudo_range.jl")
    include("user_position.jl")
    include("sat_position.jl")
    include("sv_time.jl")
    include("errors.jl")
end #end module


text:
immutable DummyPotential{V} <: AtomicDFT.AbstractPotential end

@testset "Check expressions by creating simple calculator" begin
    calc = a::DummyPotential -> typeof(a).parameters[1]
    (::typeof(calc))(a::AtomicDFT.CoeffPotential) = a.coefficient * calc(a.potential)
    (::typeof(calc))(a::AtomicDFT.PotentialExponent) = calc(a.potential) ^ a.exponent
    (::typeof(calc))(a::AtomicDFT.PotentialFold) = mapreduce(calc, a.op, a._)
    (::typeof(calc))(a::AtomicDFT.PotentialFunction) = (a.func).(calc(a.potential))
    (::typeof(calc))(a::AtomicDFT.ConstantPotential) = a._

    const zero = DummyPotential{0}()
    const one = DummyPotential{1}()
    const two = DummyPotential{2}()

    @test calc(zero) == 0
    @test calc(one) == 1
    @test calc(two) == 2

    @test calc(1zero) == 0
    @test calc(1one) == 1
    @test calc(1two) == 2

    @test calc(2zero) == 0
    @test calc(2one) == 2
    @test calc(two * 2) == 4
    @test calc(two // 2) == 1
    @test calc(two / 1.5) ≈ 2 / 1.5

    @test calc(one + two) == 3
    @test calc(zero + one) == 1
    @test calc(zero + two + one) == 3
    @test calc(one + 1two) == 3
    @test calc(one + 3two) == 7

    @test calc(one + cos(3two)) ≈ 1 + cos(6)
    @test calc(one - cos(3two)) ≈ 1 - cos(6)

    @test calc(cos(two)^2) ≈ cos(2)^2

    @test calc(two * two) == 4
    @test calc(one // two // two) == 1//2//2
    @test calc(cos(two) * sin(one)) ≈ cos(2) * sin(1)
end


text:
module CubedSphereDiscreteModelsTestsMPI
  using PartitionedArrays
  using Test
  using FillArrays
  using Gridap
  using GridapGeosciences

  include("../CubedSphereDiscreteModelsTests.jl")
  include("../ConvergenceAnalysisTools.jl")
  function main(parts)
    num_refs=[2,3,4,5]
    hs=[2.0/2^n for n in num_refs]
    model_args_series=zip(Fill(parts,length(num_refs)),num_refs)
    a,b,s=convergence_discrete_cubed_sphere_surface(hs,model_args_series,4)
    @test round(s,digits=1) ≈ 5.6
  end
  prun(main,mpi,4)
end #module

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
#=
Data aggregation of any kind. Includes handling of missing values as well.
=#
#########################################################################
# missing values handling
#########################################################################
export missing_weights, missing_val

function nomissing(da::AbstractArray{Union{T,Missing},N}, check = any(ismissing, da)) where {T,N}
    check && error("array contains missing values")
    return Array{T,N}(da)
end
nomissing(da::AbstractArray{<:Union{Real, Dates.TimeType}}, args...) = da
nomissing(da::ClimArray) = ClimArray(nomissing(da.data), da.dims, da.refdims, da.name, da.attrib)

"""
    missing_weights(A::ClimArray, val = missing_val(A)) → B, W
Generate a new array `B` with values like `A`, but with `A`'s `missing` values replaced
with `val`. Also generate an array of weights, which has the value 0 when `A` had `missing`,
and the value `1` otherwise.

The output of this function should be used in conjunction with any of ClimateBase.jl
aggregating functions like `spacemean, timemean, ...`, when your data have `missing`
values which you want to _completely skip_ during the aggregation process.

This function returns `A, nothing` if `A` has no `missing` values.
"""
function missing_weights(A::ClimArray{Union{T, Missing}}, val = missing_val(A)) where {T}
    B = zeros(T, size(A))
    W = ones(T, size(A))
    missing_idxs = findall(ismissing, A)
    notmissing_idxs = findall(!ismissing, A)
    B[missing_idxs] .= val
    B[notmissing_idxs] .= A.data[notmissing_idxs]
    W[missing_idxs] .= 0
    return ClimArray(B, dims(A); name = A.name, attrib = A.attrib),
           ClimArray(W, dims(A); name = "weights_for_missing")
end
missing_weights(A::ClimArray{<:Number}, val = nothing) = A, nothing

"""
    missing_val(A)
Return the value that represents "missing" data in `A`, according to `A`'s metadata.
If `A` does not have the `_FillValue` metadata, return 0 instead.
"""
function missing_val(A)
    if A.attrib isa Dict
        return get(A.attrib, "_FillValue", 0)
    else
        return 0
    end
end

#########################################################################
# Aggregation of data, dropagg missings, dimensions, etc.
#########################################################################
# TODO: Define dropagg with weights as well (simple to do, same as timeagg but
# without the complexity of temporal averaging)

export dropagg, nomissing, collapse, drop

"""
    dropagg(f, A, d [, W])
Apply statistics/aggregating function `f` (e.g. `sum` or `mean`) on array `A` across
dimension(s) `d` and drop the corresponding dimension(s) from the result
(Julia inherently keeps singleton dimensions).

If `A` is one dimensional, `dropagg` will return the single number of applying `f(A)`.

Optionally you can provide statistical weights in the form of an array `W`.
`W` must have same size as `A`. An exception is when `d` is only a single
dimension, e.g. `d = Lat`; then `W` is also allowed to be a single vector with
length the same as `dims(A, d)`.
"""
function dropagg(f, A, dims)
    length(size(A)) == 1 && return f(A)
    return dropdims(f(A; dims = dims); dims = dims)
end
# this method is necessary because of "reshaping" happening
# in DimensionalData.jl...
function dropagg(f, A::AbDimArray, dims)
    length(size(A)) == 1 && return f(A)
    r = dropdims(f(A; dims = dims); dims = dims)
    DimensionalData.rebuild(r, Array(r.data))
end

dropagg(f, A::ClimArray, d, w::Nothing) = dropagg(f, A, d)

function dropagg(f, A::ClimArray, d, W)
    odims = otherdims(A, d)
    oidxs = otheridxs(A, d)
    D = length(odims)
    if D == 0 # operation output is a single number
        return f(A, weights(W))
    elseif size(A) == size(W)
        r = map(i -> f(view(A, i), weights(view(W, i))), oidxs)
    elseif length(size(W)) == 1 && length(W) == length(dims(A, d))
        fw = weights(W)
        r = map(i -> f(view(A, i), fw), oidxs)
    else
        error("Given weights `W` have invalid form.")
    end
    return ClimArray(r, odims; name = A.name)
end

"""
    collapse(f, A, dim)
Reduce `A` towards dimension `dim` using the collapsing function `f` (e.g. `mean`).
This means that `f` is applied across all other dimensions of `A`, each of which are
subsequently dropped, leaving only the collapsed result of `A` vs. the remaining dimension.
"""
collapse(f, A, dim) = dropagg(f, A, otherdims(A, dim))

#########################################################################
# Other dimensions
#########################################################################
export otherdims
export otheridxs

"""
    otheridxs(A::ClimArray, Dim)
Return an iterator of indices, that when used can access all indices of `A` *except* those
belonging to dimension(s) `Dim`.

For example, if `A` has dims `(Lon, Lat, Time)` you can get all timeseries of `A`:
```julia
for i in otheridxs(A, Time())
    x = A[i...] # this is a timeseries (Vector) for each lon × lat combination
end
```
or all time+latitude slices with
```julia
for i in otheridxs(A, (Time(), Lat()))
    x = A[i...] # matrix of time × latitude slice for each longitude
end
```
(notice that splatting `i...` is necessary because `otheridxs` generates tuples)
"""
function otheridxs(A, D)
    z = otherdims(A, D)
    az = DimensionalData.basetypeof.(z)
    iters = [(Dim(i) for i in 1:size(A, Dim)) for Dim in az]
    return Iterators.product(iters...)
end

# This is the version @rafaqz suggested, but has worse performance...
# otheridxs(A, D) = map(identity, DimensionalData.dimwise_generators(otherdims(A, D)))

#########################################################################
# Dimensionwise
#########################################################################
#=
"""
    dimwise(f, A, B)
Apply bivariate function `f` (e.g. multiplication `*`) across the (first) matching dimension
of `A` and `B`, where `B` has only one dimension.
For normal Julia arrays match is done on the first matching dimension size.

For example, if `A` has dimensions (Lon, Lat, Time) and `B` has dimension (Lat,)
(and of course both Lat are actually the same thing),
then the result would be an array `C` with same structure as `A` but with values
```julia
@. C[:, i, :] = f(A[:, i, :], B[i])
```
for every `i` across the matching dimension.

Useful when wanting to weight a field `A` based on its latitude or time.
"""
function dimwise(f, A::AbDimArray, B::AbDimArray)
    @assert length(size(B)) == 1 "3rd argument of dimwise must be a vector"
    m = find_matching_dim(A, B)
    # Array(data) necessary because reshapes somewhow fail to play well with the code...
    C = dimwise(f, Array(data(A)), Array(data(B)), m)
    return DimensionalData.basetypeof(A)(C, dims(A))
end

dimwise(f, A::AbstractArray, B, m::Int = find_matching_dim(A, B)) =
dimwise!(copy(A), f, A, B, m)

function dimwise!(C, f, A::AbstractArray{T, N}, B::AbstractVector, m::Int) where {T, N}
    shape = ntuple(i -> (i==m) ? length(B) : 1, N)
    b = reshape(B, shape)
    C .= f.(A, b)
    return C
end

function find_matching_dim(A::AbDimArray, B::AbDimArray)::Int
    dimB = dims(B)[1]
    m = findfirst(isequal(dimB), dims(A))
    isnothing(m) && error("A and B have no matching dimensions.")
    return m
end

function find_matching_dim(A::AbstractArray, B::AbstractVector)::Int
    m = findfirst(isequal(length(B)), size(A))
    isnothing(m) && error("A and B have no matching dimensions.")
    return m
end

=#
export dimwise


text:
__precompile__(true)
module KeysightInstruments
# Based on user's guide for M31XXA/M33XXA Digitizers:
# http://literature.cdn.keysight.com/litweb/pdf/M3100-90002.pdf?id=2796080
# and for M32XXA/M33XXA AWGs:
# http://literature.cdn.keysight.com/litweb/pdf/M3201-90001.pdf?id=2787170

"""
Load libraries: First need to install Keysight SD1 software from
http://www.keysight.com/main/software.jspx?cc=US&lc=eng&nid=-33786.536905585&id=2784055&pageMode=CV
"""
const klib = is_windows() ? "SD1core.dll" : "libSD1core.so"

# Keysight Constants
include("KeysightConstants.jl")
# SD_Module Functions
include("SD_Module.jl")
# SD_AIN Functions
include("SD_AIN.jl")
# SD_AOU Functions
include("SD_AOU.jl")
# SD_Wave Functions
include("SD_Wave.jl")
# SD_HVI Functions
include("SD_HVI.jl")

# Keysight Errors
include("KeysightErrors.jl")
end # module


text:
@testset "N" begin
    @testset "namfrm" begin
        @test namfrm("J2000") == 1
        @test_throws SpiceError namfrm("")
    end
    @testset "ncpos" begin
        string = "BOB, JOHN, TED, AND MARTIN    "
        chars = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
        @test ncpos(string, chars, 1) == 4
        @test ncpos(string, chars, 5) == 5
        @test ncpos(string, chars, 6) == 10
        @test ncpos(string, chars, 11) == 11
        @test ncpos(string, chars, 12) == 15
        @test ncpos(string, chars, 16) == 16
        @test ncpos(string, chars, 17) == 20
        @test ncpos(string, chars, 21) == 27
        @test ncpos(string, chars, 28) == 28
        @test ncpos(string, chars, 29) == 29
        @test ncpos(string, chars, 30) == 30
        @test ncpos(string, chars, -11) == 4
        @test ncpos(string, chars, 0) == 4
        @test ncpos(string, chars, 31) == 0
        @test ncpos(string, chars, 123) == 0
    end
    @testset "ncposr" begin
        string = "BOB, JOHN, TED, AND MARTIN...."
        chars = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
        @test ncposr(string, chars, 30) == 30
        @test ncposr(string, chars, 29) == 29
        @test ncposr(string, chars, 28) == 28
        @test ncposr(string, chars, 27) == 27
        @test ncposr(string, chars, 26) == 20
        @test ncposr(string, chars, 19) == 16
        @test ncposr(string, chars, 15) == 15
        @test ncposr(string, chars, 14) == 11
        @test ncposr(string, chars, 10) == 10
        @test ncposr(string, chars, 9) == 5
        @test ncposr(string, chars, 4) == 4
        @test ncposr(string, chars, 3) == 0
        @test ncposr(string, chars, 0) == 0
        @test ncposr(string, chars, -4) == 0
        @test ncposr(string, chars, 31) == 30
        @test ncposr(string, chars, 123) == 30
    end
    @testset "nearpt" begin
        a, b, c = 1.0, 2.0, 3.0
        point = [3.5, 0.0, 0.0]
        pnear, alt = nearpt(point, a, b, c)
        expected_pnear = [1.0, 0.0, 0.0]
        expected_alt = 2.5
        @test alt ≈ expected_alt
        @test pnear ≈ expected_pnear
        @test_throws SpiceError nearpt(point, -1.0, b, c)
        @test_throws ArgumentError nearpt([1.0, 2.0], a, b, c)
    end
    @testset "npedln" begin
        linept = [1.0e6, 2.0e6, 3.0e6]
        a, b, c = 7.0e5, 7.0e5, 6.0e5
        linedr = [-4.472091234e-1, -8.944182469e-1, -4.472091234e-3]
        pnear, dist = npedln(a, b, c, linept, linedr)
        expected_pnear = [-1633.3111, -3266.6222, 599991.83]
        expected_dist = 2389967.9
        @test dist ≈ expected_dist
        @test pnear ≈ expected_pnear
        @test_throws SpiceError npedln(-1.0, b, c, linept, linedr)
        @test_throws ArgumentError npedln(a, b, c, linept, linedr[1:2])
        @test_throws ArgumentError npedln(a, b, c, linept[1:2], linedr)
    end
    @testset "npelpt" begin
        center = [1.0, 2.0, 3.0]
        smajor = [3.0, 0.0, 0.0]
        sminor = [0.0, 2.0, 0.0]
        point = [-4.0, 2.0, 1.0]
        expected_pnear = [-2.0, 2.0, 3.0]
        expected_dist = 2.8284271
        ellipse = cgv2el(center, smajor, sminor)
        pnear, dist = npelpt(point, ellipse)
        @test dist ≈ expected_dist
        @test pnear ≈ expected_pnear
        @test_throws ArgumentError npelpt(point[1:2], ellipse)
    end
    @testset "nplnpt" begin
        linept = [1.0, 2.0, 3.0]
        linedr = [0.0, 1.0, 1.0]
        point = [-6.0, 9.0, 10.0]
        pnear, dist = nplnpt(linept, linedr, point)
        expected_pnear = [1.0, 9.0, 10.0]
        expected_dist = 7.0
        @test dist ≈ expected_dist
        @test pnear ≈ expected_pnear
        @test_throws SpiceError nplnpt(linept, zeros(3), point)
        @test_throws ArgumentError nplnpt(linept[1:2], linedr, point)
        @test_throws ArgumentError nplnpt(linept, linedr[1:2], point)
        @test_throws ArgumentError nplnpt(linept, linedr, point[1:2])
    end
    @testset "nvc2pl" begin
        normal = [1.0, 1.0, 1.0]
        constant = 23.0
        expected_constant = 13.279056
        expected_normal = [0.57735027, 0.57735027, 0.57735027]
        plane = nvc2pl(normal, constant)
        @test collect(plane.normal) ≈ expected_normal
        @test plane.constant ≈ expected_constant
        @test_throws ArgumentError nvc2pl(normal[1:2], constant)
        @test_throws SpiceError nvc2pl(zeros(3), constant)
    end
    @testset "nvp2pl" begin
        normal = [1.0, 1.0, 1.0]
        point = [1.0, 4.0, 9.0]
        expected_constant = 8.0829038
        expected_normal = [0.57735027, 0.57735027, 0.57735027]
        plane = nvp2pl(normal, point)
        @test collect(plane.normal) ≈ expected_normal
        @test plane.constant ≈ expected_constant
        @test_throws ArgumentError nvp2pl(normal[1:2], point)
        @test_throws ArgumentError nvp2pl(normal, point[1:2])
        @test_throws SpiceError nvp2pl(zeros(3), point)
    end
end


text:
import ConstrainedSystems: init

"""
    ODEFunctionList

## Constructor

`ODEFunctionList(state=nothing,constraint=nothing,ode_rhs=nothing,bc_rhs=nothing,constraint_force=nothing,bc_op=nothing,lin_op=nothing)`

Supply functions and data types for a constrained ODE system. The specified
functions must provide the various parts that comprise the constrained ODE system

``
\\dfrac{dy}{dt} = Ly - B_1 z + r_1(y,t)
``

``
B_2 y = r_2(t)
``

The functions can be in in-place or out-of-place form, but they must all be consistently
of the same form.

- `state =` specifies a prototype of the state vector ``y``.

- `constraint = ` specifies a prototype of the constraint force vector ``z``. If there are no constraints, this can be omitted.

- `ode_rhs =` specifies the right-hand side of the ODEs, ``r_1``. The in-place form of the function is `r1(dy,y,sys::ILMSystem,t)`, the state vector `y`, IL system `sys`, time `t`, and returning ``dy/dt``. The out-of-place form is `r1(y,sys,t)`.

- `bc_rhs = ` specifies the right-hand side of the boundary conditions, ``r_2``. If there are no constraints, this can be omitted. The in-place form is `r2(dz,sys,t)`, returning `dz`, the part of the boundary constraint not dependent on the state vector. The out-of-place form is `r2(sys,t)`.

- `constraint_force = ` supplies the constraint force term in the ODEs, ``B_1 z``. If there are no constraints, this can be omitted. The in-place form is `B1(dy,z,sys)`, returning the contribution to ``dy/dt`` (with the sign convention shown in the equations above) and the out-of-place form is `B1(z,sys)`.

- `bc_op = ` supplies the left-hand side of the boundary constraint, ``B_2 y``. If there are no constraints, this can be omitted. The in-place form is `B2(dz,y,sys)`, the out-of-place form is `B2(y,sys)`.

- `lin_op = ` is optional and specifies a linear operator on the state vector ``L``, to be treated with an exponential integral (i.e., integrating factor) in the time marching. (Alternatively, this part can simply be included in `r_1`). It should have an associated `mul!` operation that acts upon the state vector.

"""
struct ODEFunctionList{RT,LT,BRT,CFT,BOT,ST,CT}
    ode_rhs :: RT
    lin_op :: LT
    bc_rhs :: BRT
    constraint_force :: CFT
    bc_op :: BOT
    state :: ST
    constraint :: CT
end

_has_lin_op(f::ODEFunctionList{RT,LT}) where {RT,LT} = true
_has_lin_op(f::ODEFunctionList{RT,Nothing}) where {RT} = false

function ODEFunctionList(;ode_rhs=nothing,lin_op=nothing,bc_rhs=nothing,constraint_force=nothing,bc_op=nothing,state=nothing,constraint=nothing)

    # Audit the supplied information to make sure it is consistent
    !(state isa Nothing)  || error("need to supply a state vector")
    !(ode_rhs isa Nothing)  || error("need to supply a rhs function for the ODE")
    constraint_stuff = [bc_rhs,constraint_force,bc_op,constraint]
    #!any(i -> i isa Nothing,constraint_stuff) || error("incomplete set of functions for constrained system")

    ODEFunctionList{typeof(ode_rhs),typeof(lin_op),typeof(bc_rhs),typeof(constraint_force),
                    typeof(bc_op),typeof(state),typeof(constraint)}(ode_rhs,lin_op,bc_rhs,constraint_force,bc_op,state,constraint)
end

function ImmersedLayers.ConstrainedODEFunction(sys::ILMSystem{true,SCA,0}) where {SCA}
    @unpack extra_cache = sys
    @unpack f = extra_cache

    _constrained_ode_function(f.lin_op,f.ode_rhs;_func_cache=zeros_sol(sys))
end

function ImmersedLayers.ConstrainedODEFunction(sys::ILMSystem{true})
    @unpack extra_cache = sys
    @unpack f = extra_cache

    _constrained_ode_function(f.lin_op,f.ode_rhs,f.bc_rhs,f.constraint_force,
                           f.bc_op;_func_cache=zeros_sol(sys))
end

function ImmersedLayers.ConstrainedODEFunction(sys::ILMSystem{false})
    @unpack extra_cache = sys
    @unpack f = extra_cache

    rhs! = ConstrainedSystems.r1vector(state_r1 = f.ode_rhs,
                                       aux_r1 = body_rhs!)

    _constrained_ode_function(f.lin_op,rhs!,f.bc_rhs,f.constraint_force,
                              f.bc_op;_func_cache=zeros_sol(sys),
                                      param_update_func=update_system!)
end

@inline _constrained_ode_function(lin_op,args...;kwargs...) =
    ConstrainedODEFunction(args...,lin_op;kwargs...)

@inline _constrained_ode_function(::Nothing,args...;kwargs...) =
    ConstrainedODEFunction(args...;kwargs...)

"""
    zeros_sol(sys::ILMSystem)

Return a zeroed version of the solution vector.
"""
function zeros_sol(sys::ILMSystem{true,SCA,0}) where {SCA}
    @unpack extra_cache = sys
    @unpack f = extra_cache
    return solvector(state=zero(f.state))
end

function zeros_sol(sys::ILMSystem{true})
    @unpack extra_cache = sys
    @unpack f = extra_cache
    return solvector(state=zero(f.state),constraint=zero(f.constraint))
end

function zeros_sol(sys::ILMSystem{false})
    @unpack motions, extra_cache, base_cache = sys
    @unpack bl  = base_cache
    @unpack f = extra_cache
    return solvector(state=zero(f.state),
                     constraint=zero(f.constraint),
                     aux_state=zero(motion_state(bl,motions)))
end

"""
    init_sol(sys::ILMSystem)

Return the initial solution vector.
"""
function init_sol(sys::ILMSystem{true,SCA,0}) where {SCA}
    @unpack extra_cache = sys
    @unpack f = extra_cache
    return solvector(state=zero(f.state))
end

function init_sol(sys::ILMSystem{true})
    @unpack extra_cache = sys
    @unpack f = extra_cache
    return solvector(state=zero(f.state),constraint=zero(f.constraint))
end

function init_sol(sys::ILMSystem{false})
    @unpack motions, extra_cache, base_cache = sys
    @unpack bl  = base_cache
    @unpack f = extra_cache
    return solvector(state=zero(f.state),
                     constraint=zero(f.constraint),
                     aux_state=motion_state(bl,motions))
end

function body_rhs!(dx::Vector{T},x::Vector{T},sys::ILMSystem,t::Real) where {T<:Real}
  @unpack motions, base_cache = sys
  @unpack bl = base_cache
  length(dx) == length(x) || error("wrong length for vector")
  dx .= motion_velocity(bl,motions,t)
  return dx
end

_norm_sq(u) = dot(u,u)
_norm_sq(u::ConstrainedSystems.ArrayPartition) = sum(_norm_sq,u.x)
state_norm(u,t) = sqrt(_norm_sq(u))

"""
    ConstrainedSystems.init(u0,tspan,sys::ILMSystem,[alg=ConstrainedSystems.LiskaIFHERK()])

Initialize the integrator for a time-varying immersed-layer system of PDEs,
described in `sys`.
"""
function init(u0,tspan,sys::ILMSystem;alg=ConstrainedSystems.LiskaIFHERK(),kwargs...)
    @unpack timestep_func, phys_params, extra_cache,base_cache = sys
    @unpack g = base_cache
    fode = ConstrainedODEFunction(sys)

    prob = ODEProblem(fode,u0,tspan,sys)
    dt_calc = timestep_func(g,phys_params)
    return init(prob, alg,dt=dt_calc,internal_norm=state_norm,kwargs...)
end


text:
#=
    Tiled Cholesky factorization in pure Julia. The serial performance
    essentially comes from the `TriangularSolve` and `LoopVectorization`
    packages. The parallelizatino is handled by `DataFlowTask`s.
=#

cholesky(A::Matrix,args...) = cholesky!(deepcopy(A),args...)

function cholesky!(A::Matrix,s=TILESIZE[],tturbo::Val{T}=Val(false)) where {T}
    _cholesky!(PseudoTiledMatrix(A,s),tturbo)
end

# tiled cholesky factorization
function _cholesky!(A::PseudoTiledMatrix,tturbo::Val{T}=Val(false)) where {T}
    m,n = size(A) # number of blocks
    for i in 1:m
        # _chol!(A[i,i],UpperTriangular)
        @dspawn _chol!(A[i,i],UpperTriangular,tturbo) (A[i,i],) (RW,)
        Aii = A[i,i]
        U = UpperTriangular(Aii)
        L = adjoint(U)
        for j in i+1:n
            Aij = A[i,j]
            # TriangularSolve.ldiv!(L,Aij,tturbo)
            @dspawn TriangularSolve.ldiv!(L,Aij,tturbo) (Aii,Aij) (R,RW)
        end
        for j in i+1:m
            Aij = A[i,j]
            for k in j:n
                # TODO: for k = j, only the upper part needs to be updated,
                # dividing the cost of that operation by two
                Ajk = A[j,k]
                Aji = adjoint(Aij)
                Aik = A[i,k]
                # schur_complement!(Ajk,Aji,Aik,tturbo)
                @dspawn schur_complement!(Ajk,Aji,Aik,tturbo) (Ajk,Aij,Aik) (RW,R,R)
            end
        end
    end
    # wait for all computations before returning
    DataFlowTasks.sync()
    return Cholesky(A.data,'U',zero(BlasInt))
end

# a fork-join approach for comparison with the data-flow parallelism
function _cholesky_forkjoin!(A::PseudoTiledMatrix,tturbo::Val{T}=Val(false)) where {T}
    m,n = size(A) # number of blocks
    for i in 1:m
        _chol!(A[i,i],UpperTriangular,tturbo)
        Aii = A[i,i]
        U = UpperTriangular(Aii)
        L = adjoint(U)
        Threads.@threads for j in i+1:n
            Aij = A[i,j]
            TriangularSolve.ldiv!(L,Aij,tturbo)
        end
        # spawn m*(m+1)/2 tasks and sync them at the end
        @sync for j in i+1:m
            Aij = A[i,j]
            for k in j:n
                Ajk = A[j,k]
                Aji = adjoint(Aij)
                Aik = A[i,k]
                Threads.@spawn schur_complement!(Ajk,Aji,Aik,tturbo)
            end
        end
    end
    return Cholesky(A.data,'U',zero(Int32))
end

# Modified from the generic version from LinearAlgebra (MIT license).
function _chol!(A::AbstractMatrix{<:Real}, ::Type{UpperTriangular},tturbo::Val{T}=Val(false)) where {T}
    Base.require_one_based_indexing(A)
    n = LinearAlgebra.checksquare(A)
    @inbounds begin
        for k = 1:n
            Akk = A[k,k]
            for i = 1:k - 1
                Akk -= A[i,k]*A[i,k]
            end
            A[k,k] = Akk
            Akk, info = _chol!(Akk, UpperTriangular)
            if info != 0
                return UpperTriangular(A), info
            end
            A[k,k] = Akk
            AkkInv = inv(Akk')
            if T
                @tturbo warn_check_args=false for j = k + 1:n
                    for i = 1:k - 1
                        A[k,j] -= A[i,k]*A[i,j]
                    end
                end
                @tturbo warn_check_args=false for j in k+1:n
                    A[k,j] = AkkInv*A[k,j]
                end
            else
                @turbo warn_check_args=false for j = k + 1:n
                    for i = 1:k - 1
                        A[k,j] -= A[i,k]*A[i,j]
                    end
                end
                @turbo warn_check_args=false for j in k+1:n
                    A[k,j] = AkkInv*A[k,j]
                end
            end
        end
    end
    return UpperTriangular(A), convert(Int32, 0)
end
## Numbers
function _chol!(x::Number, uplo)
    rx = real(x)
    rxr = sqrt(abs(rx))
    rval =  convert(promote_type(typeof(x), typeof(rxr)), rxr)
    rx == abs(x) ? (rval, convert(Int32, 0)) : (rval, convert(Int32, 1))
end


text:
#########################################################################
# Project #8: DIIS extrapolation for the SCF procedure.
# - P. Pulay, Chem. Phys. Lett. 73, 393-398 (1980).
# - P. Pulay, J. Comp. Chem. 3, 556-560 (1982).
# - T. P. Hamilton and P. Pulay, J. Chem. Phys. 84, 5728-5734 (1986).
# - C. David Sherrill. "Some comments on accellerating convergence of iterative
#      sequences using direct inversion of the iterative subspace (DIIS)".
#      Available at: vergil.chemistry.gatech.edu/notes/diis/diis.pdf. (1998)
#########################################################################

function diis(F, Fstack, DIISstack, s_minushalf, s, P, iteri)
   # orbital gradient in AO basis: F Di Si - S Di Fi
   # better choice:     (S^-0.5)' (F Di Si - S Di Fi) S^-0.5
   DIIS_residual = s * P * F
   DIIS_residual = DIIS_residual' - DIIS_residual
   DIIS_residual = s_minushalf * DIIS_residual * s_minushalf
   if iteri > 1
      push!(Fstack,F)
      push!(DIISstack, DIIS_residual)
   end
   ΔDIIS = norm(DIIS_residual)

   if iteri > 2
      if iteri > 15 # Limit the storage of DIIS arrays
         popfirst!(DIISstack)
         popfirst!(Fstack)
      end
      dim_B = length(Fstack) + 1 # 1 dim for the Lagrange multiplier
      B = zeros(dim_B,dim_B)
      B[end,:]   .= -1
      B[:,  end] .= -1
      B[end,end]  =  0
      for i in eachindex(Fstack), j in eachindex(Fstack)
         B[i,j] = dot(DIISstack[i], DIISstack[j]) # Gramian matrix
      end
      # Solve Lagrange equation of Pulay
      Pulay_rhs = zeros(dim_B)
      Pulay_rhs[end] = -1
      coef_Pulay = B \ Pulay_rhs
      F = zeros(size(F))
      for (Fi,c) in zip(Fstack,coef_Pulay[1:end-1]) # skip Lagrange mult.
         F += c * Fi
      end
   end
   F, Fstack, DIISstack
end


text:
# test leap-frogging vortex rings in 2d
include("vm2d.jl")

# create vortex ring geometry
radius1 = 0.1
radius2 = 0.1
spacing = 0.05
points = [[radius1,0.0,0.0],[-radius1,0.0,0.0],[radius2,0.0,spacing],[-radius2,0.0,spacing]]

# set vectorial vortex strengths
gamma = 0.1
gamma_upper = [0,gamma,0]
gamma_lower = [0,-gamma,0]
gammas = [gamma_upper,gamma_lower,gamma_upper,gamma_lower]

# set initial velocities
v0 = 0.0
vs = [[v0,0.0,0.0],[v0,0.0,0.0],[0.0,0.0,0.0],[0.0,0.0,0.0]]

particles = place(points,gammas,vs)

# set simulation variables
name = "vortexrings"
global currenttime
currenttime = 0.0
timestep = 0.01
numtimesteps = 5 # try 150

# animate!

# set up plot
plt.plot()
plt.hot()
plt.axis("equal")
for particle in particles
    plt.scatter(particle.x[1],particle.x[2],c=particle.gamma[3])
end
plt.title("time: "*string(round(currenttime; digits=3)))
if ~isdir("./sim")
    mkdir("./sim")
end
# animate
for i in 1:numtimesteps
    ## advance timestep
    advance(particles,timestep,Uinf)
    currenttime += timestep
    ## evaluate velocity field
    zs = range(-radius1, stop=radius1, length=6)
    ys = [0.0] #range(-2*radius1,stop=2*radius1,length=20)
    xs = range(particles[1].x[1] - 4*spacing,
            stop = particles[1].x[1] + 4*spacing,
            length = 11)
    (x,y,z,u,v,w) = vfield(particles, xs, ys, zs)
    ## update plot
    name = "getworking"*string(i)
    plot_2d(particles, "zx"; velocity=true, xs=xs, ys=ys, zs=zs, title="time: "*string(round(currenttime; digits=3)), save=true, name=name)
    # plt.quiver(x,y,u,v)
    plt.pause(0.005)
end


text:
using Revise, MathematicalSystems,Reachability, LinearAlgebra, HybridSystems
using LazySets;
using Reachability: solve
using Plots
using TaylorIntegration

@taylorize function spiking_neuron!(t, x, dx)
    local a = 0.02
    local b = 0.2
    local I = 40
    dx[1] = 0.04*x[1]^2 + 5*x[1] + 140 - x[2] + I
    dx[2] = a*(b*x[1] - x[2])
end

function spiking_neurons()

    automaton = LightAutomaton(1) # one node

    inv_one = HPolyhedron([HalfSpace([-1.0, 0.0], 0.2),    #x[1] >= -0.2
                           HalfSpace([0.0, -1.0], 70.0)])  #x[2] >= -70

    m1 = ConstrainedBlackBoxContinuousSystem(spiking_neuron!, 2, inv_one)

    modes = [m1]

    add_transition!(automaton, 1, 1, 1)

    guard_alpha = HPolyhedron([HalfSpace([1.0, 0.0], 30.0),   # x[1] >= 30
                               HalfSpace([-1.0, 0.0], -30.0)]) # x[1] <= 30

    A = [1.0 0.0;0.0 1.0]
    b = [-65.0, 8.0]
    t1 = ConstrainedAffineMap(A,b,guard_alpha)
    #resetmaps
    resetmaps = [t1]
    # switching
    switching = [AutonomousSwitching()]

    ℋ = HybridSystem(automaton, modes, resetmaps, switching)
    #initial condition in mode one
    X0 = Hyperrectangle(low = [-0.2, -65.0], high = [0.2, -60.0])

    initial_condition = [(1, X0)]

    system = InitialValueProblem(ℋ, initial_condition)

    options = Options(:mode=>"reach", :T=>100.0, :plot_vars=>[1, 2], :project_reachset=>false)

    return (system, options)
end

problem, options = spiking_neurons();

options = Options(:mode=>"reach", :T=>100.0, :plot_vars=>[1, 2], :project_reachset=>false, :verbosity => "info")

@time sol_TMJets = solve(problem, options, TMJets(:orderT=>5, :orderQ=>2, :abs_tol=>1e-10),LazyDiscretePost(:check_invariant_intersection=>true))

plot(sol_TMJets, use_subindices=false, alpha=.5, color=:lightblue)


text:
using LiveServer

sw = LiveServer.SimpleWatcher()

# source files to check for updates to avoid infinite loop
EXAMPLEDIR = joinpath(@__DIR__, "src", "examples")
EXAMPLES_jl = [f for f in readdir(EXAMPLEDIR) if endswith(f, ".jl")]
SOURCES = [
    abspath(joinpath(@__DIR__, "make.jl"))
    abspath(joinpath(@__DIR__, "src", "index.md"))
    abspath(joinpath(@__DIR__, "src", "prerequisites.md"))
    abspath(joinpath(@__DIR__, "src", "functions.md"))
    [abspath(joinpath(EXAMPLEDIR, example)) for example in EXAMPLES_jl]
]

append!(sw.watchedfiles, LiveServer.WatchedFile.(SOURCES))

function callback(x)
    # only trigger for source files to avoid infinite loop
    if x in SOURCES
        include(abspath(joinpath(@__DIR__, "make.jl")))
        LiveServer.file_changed_callback(x)
    end
end

callback(abspath(joinpath(@__DIR__, "make.jl"))) # make sure files exist
LiveServer.set_callback!(sw, callback)

LiveServer.serve(sw; dir = abspath(joinpath(@__DIR__, "build")), verbose=true)