text:
function year_to_short_string(y)
    y = y % 10
    return "$y"
end

function number_to_char(n::Number)::Char
    if n >= 0 && n <= 9
        return Char(codepoint('0') + n)
    else
        return Char(codepoint('A') + n - 10)
    end
end

function char_to_number(c::Char)::UInt32
    if c >= '0' && c <= '9'
        return codepoint(c) - codepoint('0')
    else
        return codepoint(c) - codepoint('A') + 10
    end
end

struct IGCShortFilename
    date::Date
    manufacturer::Char
    serial_id::String
    flight_of_day::UInt64
end

struct IGCLongFilename
    date::Date
    manufacturer::String
    serial_id::String
    flight_of_day::UInt64
end

"""
    filename(f::IGCShortFilename)::String

Return a short filename as string
YMDCXXXF.IGC
    Y: year; value 0 to 9, cycling every 10 years
    M: month; value 1 to 9 then A for 10, B=11, C=12
    D: day; value 1 to 9 then A=10, B=11, C=12, D=13, E=14, F=15, G=16, H=17, I=18, J=19,
        K=20, L=21, M=22, N=23, O=24, P=25, Q=26, R=27, S=28, T=29, U=30, V=31.
    C: manufacturer's single-letter IGC identifier
    XXX: unique flight recorder Serial ID (S/ID); 3 alphanumeric characters

Example
56HCXXX2.IGC
"""
function filename(f::IGCShortFilename)::String
    d = f.date
    Y = year_to_short_string(Dates.year(d))
    M = number_to_char(Dates.month(d))
    D = number_to_char(Dates.day(d))
    C = f.manufacturer[1]
    XXX = f.serial_id[1:3]
    F = number_to_char(f.flight_of_day)
    ext = ".IGC"
    return Y * M * D * C * XXX * F * ext
end

function decade(y)
    return y - y % 10
end

function parse(::Type{IGCShortFilename}, s::String; year=Dates.year(Dates.today()))
    year = decade(year)
    pattern = r"(.)(.)(.)(.)(...)(.).IGC"
    m = match(pattern, s)
    Y = year + parse(Int, m[1])
    M = char_to_number(m[2][1])
    D = char_to_number(m[3][1])
    C = m[4][1]
    XXX = m[5]
    F = char_to_number(m[6][1])
    return IGCShortFilename(
        Date(Y, M, D),
        C,
        XXX,
        F
    )
end

"""
    filename(f::IGCLongFilename)::String

Return a long filename as string
YYYY-MM-DD-MMM-XXX-FF.IGC
    YYYY: year
    MM: month
    DD: day
    MMM: manufacturer
    XXX: unique flight recorder Serial ID (S/ID); 3 alphanumeric characters
    FF: flight of the day

Example
2015-06-17-MMM-XXX-02.IGC
"""
function filename(f::IGCLongFilename)::String
    d = f.date
    YYYY = lpad(Dates.year(d) % 10000, 4, "0")
    MM = lpad(Dates.month(d), 2, "0")
    DD = lpad(Dates.day(d), 2, "0")
    MMM = f.manufacturer[1:3]
    XXX = f.serial_id[1:3]
    FF = lpad(f.flight_of_day % 100, 2, "0")
    return "$YYYY-$MM-$DD-$MMM-$XXX-$FF.IGC"
end

function parse(::Type{IGCLongFilename}, s::String)
    pattern = r"(\d\d\d\d)-(\d\d)-(\d\d)-(.{1,3})-(.{1,3})-(\d\d).IGC"
    m = match(pattern, s)
    YYYY = parse(Int, m[1])
    MM = parse(Int, m[2])
    DD = parse(Int, m[3])
    MMM = m[4]
    XXX = m[5]
    FF = parse(Int, m[6])
    return IGCLongFilename(
        Date(YYYY, MM, DD),
        MMM,
        XXX,
        FF
    )
end


text:
using Smolyak

const verbose = true
const print_str = "mean abs diff: %.3e\nmax abs diff: %.3e\nmin abs diff: %.3e"
@eval print_func(x, y, z) = @printf($print_str, x, y, z)

# L2-norm squared
f1(x) = squeeze(sum(x .^ 2, 2), 2)
f1_p(x) = 2 * x

f2(x) = sin(x[:,1]) + cos(x[:, 2])
f2_p(x) = [cos(x[:, 1]) -sin(x[:, 2])]

function test_interp2d_derivs(d::Int, mu::Int, f::Function, f_prime::Function, bds::Real=1)
    sg = SmolyakGrid(d, mu, -bds, bds)

    f_on_grid = f(sg.grid)
    si = SmolyakInterp(sg, f_on_grid)

    srand(42)
    test_points = randn(100, d)

    # Make sure it is bounded by -bds, bds
    test_points = (bds - 0.05) * test_points/maximum(abs(test_points))

    true_vals = f(test_points)
    true_vals_prime = f_prime(test_points)
    i_vals = si(test_points)

    mean_ad = mean(abs(i_vals - true_vals))
    max_ad = maximum(abs(i_vals - true_vals))
    min_ad = minimum(abs(i_vals - true_vals))

    if verbose
        println("Interpolation results\n" * "#" ^ 21)
        print_func(mean_ad, max_ad, min_ad)
    end

    # return i_vals_prime
end

test_interp2d_derivs(15, 3, f2, f2_p, -2)

@time test_interp2d_derivs(15, 3, f2, f2_p, -2)

ext: jl

max_issues_repo_name: EconForge/Smolyak

max_forks_repo_name: EconForge/Smolyak


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _aitoff{T<:AbstractFloat}(l::T, b::T)
    l > 180.0 && (l -= 360.0)
    alpha2 = deg2rad(l/2.0)
    delta = deg2rad(b)
    r2 = sqrt(2.0)
    f = 2.0*r2/pi
    cdec = cos(delta)
    denom = sqrt(1.0 + cdec*cos(alpha2))
    return rad2deg(cdec*sin(alpha2)*2.0*r2/denom/f), rad2deg(sin(delta)*r2/denom/f)
end

"""
    aitoff(l, b) -> x, y

### Purpose ###

Convert longitude `l` and latitude `b` to `(x, y)` using an Aitoff projection.

### Explanation ###

This function can be used to create an all-sky map in Galactic coordinates with
an equal-area Aitoff projection.  Output map coordinates are zero longitude
centered.

### Arguments ###

* `l`: longitude, scalar or vector, in degrees.
* `b`: latitude, number of elements as `l`, in degrees.

Coordinates can be given also as a 2-tuple `(l, b)`.

### Output ###

2-tuple `(x, y)`.

* `x`: x coordinate, same number of elements as `l`.  `x` is normalized to be in
  \$[-180, 180]\$.
* `y`: y coordinate, same number of elements as `l`.  `y` is normalized to be in
  \$[-90, 90]\$.

### Example ###

Get \$(x ,y)\$ Aitoff coordinates of Sirius, whose Galactic coordinates are
\$(227.23, -8.890)\$.

``` julia
x, y = aitoff(227.23, -8.890)
# => (-137.92196683723276,-11.772527357473054)
```

### Notes ###

See AIPS memo No. 46
(ftp://ftp.aoc.nrao.edu/pub/software/aips/TEXT/PUBL/AIPSMEMO46.PS), page 4, for
details of the algorithm.  This version of `aitoff` assumes the projection is
centered at b=0 degrees.

Code of this function is based on IDL Astronomy User's Library.
"""
aitoff(l::Real, b::Real) = _aitoff(promote(float(l), float(b))...)

aitoff(lb::Tuple{Real, Real}) = aitoff(lb...)

function aitoff{L<:Real,B<:Real}(l::AbstractArray{L}, b::AbstractArray{B})
    @assert length(l) == length(b)
    typel = typeof(float(one(L)))
    x = similar(l, typel)
    y = similar(b, typel)
    for i in eachindex(l)
        x[i], y[i] = aitoff(l[i], b[i])
    end
    return x, y
end


text:
using ...semiclassical
import ...semiclassical: State

"""
    stochastic.schroedinger_semiclassical(tspan, state0, fquantum, fclassical[; fout, ...])
Integrate time-dependent Schrödinger equation coupled to a classical system.
# Arguments
* `tspan`: Vector specifying the points of time for which the output should
        be displayed.
* `state0`: Initial semi-classical state [`semiclassical.State`](@ref).
* `fquantum`: Function `f(t, psi, u) -> H` returning the time and or state
        dependent Hamiltonian.
* `fclassical`: Function `f(t, psi, u, du)` calculating the possibly time and
        state dependent derivative of the classical equations and storing it
        in the vector `du`.
* `fstoch_quantum=nothing`: Function `f(t, psi, u) -> Hs` that returns a vector
        of operators corresponding to the stochastic terms of the Hamiltonian.
        NOTE: Either this function or `fstoch_classical` has to be defined.
* `fstoch_classical=nothing`: Function `f(t, psi, u, du)` that calculates the
        stochastic terms of the derivative `du`.
        NOTE: Either this function or `fstoch_quantum` has to be defined.
* `fout=nothing`: If given, this function `fout(t, state)` is called every time
        an output should be displayed. ATTENTION: The given state is neither
        normalized nor permanent!
* `noise_processes=0`: Number of distinct quantum noise processes in the equation.
        This number has to be equal to the total number of noise operators
        returned by `fstoch`. If unset, the number is calculated automatically
        from the function output.
        NOTE: Set this number if you want to avoid an initial calculation of
        the function output!
* `noise_prototype_classical=nothing`: The equivalent of the optional argument
        `noise_rate_prototype` in `StochasticDiffEq` for the classical
        stochastic function `fstoch_classical` only. Must be set for
        non-diagonal classical noise or combinations of quantum and classical
        noise. See the documentation for details.
* `normalize_state=false`: Specify whether or not to normalize the state after
        each time step taken by the solver.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function schroedinger_semiclassical(tspan, state0::S, fquantum::Function,
                fclassical::Function; fstoch_quantum::Union{Nothing, Function}=nothing,
                fstoch_classical::Union{Nothing, Function}=nothing,
                fout::Union{Function,Nothing}=nothing,
                noise_processes::Int=0,
                noise_prototype_classical=nothing,
                normalize_state::Bool=false,
                kwargs...) where {B<:Basis,T<:Ket{B},S<:State{B,T}}
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    dschroedinger_det(t, state::S, dstate::S) =
            semiclassical.dschroedinger_dynamic(t, state, fquantum, fclassical, dstate)

    if isa(fstoch_quantum, Nothing) && isa(fstoch_classical, Nothing)
        throw(ArgumentError("No stochastic functions provided!"))
    end

    x0 = Vector{eltype(state0)}(undef, length(state0))
    recast!(state0, x0)
    state = copy(state0)
    dstate = copy(state0)

    if noise_processes == 0
        n = 0
        if isa(fstoch_quantum, Function)
            fs_out = fstoch_quantum(0.0, state0.quantum, state0.classical)
            n += length(fs_out)
        end
    else
        n = noise_processes
    end

    if n > 0 && isa(fstoch_classical, Function)
        if isa(noise_prototype_classical, Nothing)
            throw(ArgumentError("noise_prototype_classical must be set for combinations of quantum and classical noise!"))
        end
    end

    if normalize_state
        len_q = length(state0.quantum)
        function norm_func(u::Vector, t, integrator)
            u .= [normalize!(u[1:len_q]); u[len_q+1:end]]
        end
        ncb = DiffEqCallbacks.FunctionCallingCallback(norm_func;
                 func_everystep=true,
                 func_start=false)
    else
        ncb = nothing
    end

    dschroedinger_stoch(dx::DiffArray, t, state::S, dstate::S, n::Int) =
            dschroedinger_stochastic(dx, t, state, fstoch_quantum, fstoch_classical, dstate, n)

    integrate_stoch(tspan_, dschroedinger_det, dschroedinger_stoch, x0, state, dstate, fout, n;
                    noise_prototype_classical = noise_prototype_classical,
                    ncb=ncb,
                    kwargs...)
end

"""
    stochastic.master_semiclassical(tspan, rho0, H, Hs, J; <keyword arguments>)
Time-evolution according to a stochastic master equation.
For dense arguments the `master` function calculates the
non-hermitian Hamiltonian and then calls master_nh which is slightly faster.
# Arguments
* `tspan`: Vector specifying the points of time for which output should
        be displayed.
* `rho0`: Initial density operator. Can also be a state vector which is
        automatically converted into a density operator.
* `fquantum`: Function `f(t, rho, u) -> (H, J, Jdagger)` or
        `f(t, rho, u) -> (H, J, Jdagger, rates)` giving the deterministic
        part of the master equation.
* `fclassical`: Function `f(t, rho, u, du)` that calculates the classical
        derivatives `du`.
* `fstoch_quantum=nothing`: Function `f(t, rho, u) -> C, Cdagger`
        that returns the stochastic operator for the superoperator of the form
        `C[i]*rho + rho*Cdagger[i]`.
* `fstoch_classical=nothing`: Function `f(t, rho, u, du)` that calculates the
        stochastic terms of the derivative `du`.
* `rates=nothing`: Vector or matrix specifying the coefficients (decay rates)
        for the jump operators. If nothing is specified all rates are assumed
        to be 1.
* `fout=nothing`: If given, this function `fout(t, rho)` is called every time
        an output should be displayed. ATTENTION: The given state rho is not
        permanent! It is still in use by the ode solver and therefore must not
        be changed.
* `noise_processes=0`: Number of distinct quantum noise processes in the equation.
        This number has to be equal to the total number of noise operators
        returned by `fstoch`. If unset, the number is calculated automatically
        from the function output.
        NOTE: Set this number if you want to avoid an initial calculation of
        the function output!
* `noise_prototype_classical=nothing`: The equivalent of the optional argument
        `noise_rate_prototype` in `StochasticDiffEq` for the classical
        stochastic function `fstoch_classical` only. Must be set for
        non-diagonal classical noise or combinations of quantum and classical
        noise. See the documentation for details.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master_semiclassical(tspan, rho0::S,
                fquantum::Function, fclassical::Function;
                fstoch_quantum::Union{Function, Nothing}=nothing,
                fstoch_classical::Union{Function, Nothing}=nothing,
                rates::DecayRates=nothing,
                fout::Union{Function,Nothing}=nothing,
                noise_processes::Int=0,
                noise_prototype_classical=nothing,
                nonlinear::Bool=true,
                kwargs...) where {B<:Basis,T<:Operator{B,B},S<:State{B,T}}

    tmp = copy(rho0.quantum)
    if isa(fstoch_quantum, Nothing) && isa(fstoch_classical, Nothing)
        throw(ArgumentError("No stochastic functions provided!"))
    end

    if noise_processes == 0
        n = 0
        if isa(fstoch_quantum, Function)
            fq_out = fstoch_quantum(0, rho0.quantum, rho0.classical)
            n += length(fq_out[1])
        end
    else
        n = noise_processes
    end

    if n > 0 && isa(fstoch_classical, Function)
        if isa(noise_prototype_classical, Nothing)
            throw(ArgumentError("noise_prototype_classical must be set for combinations of quantum and classical noise!"))
        end
    end

    dmaster_determ(t, rho::S, drho::S) =
            semiclassical.dmaster_h_dynamic(t, rho, fquantum, fclassical, rates, drho, tmp)

    dmaster_stoch(dx::DiffArray, t, rho::S,
                    drho::S, n::Int) =
        dmaster_stoch_dynamic(dx, t, rho, fstoch_quantum, fstoch_classical, drho, n)

    integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch, rho0, fout, n;
                noise_prototype_classical=noise_prototype_classical,
                kwargs...)
end
master_semiclassical(tspan, psi0::State{B,T}, args...; kwargs...) where {B<:Basis,T<:Ket{B}} =
        master_semiclassical(tspan, dm(psi0), args...; kwargs...)

# Derivative functions
function dschroedinger_stochastic(dx::Vector, t,
        state::S, fstoch_quantum::Function, fstoch_classical::Nothing,
        dstate::S, ::Int) where {B<:Basis,T<:Ket{B},S<:State{B,T}}
    H = fstoch_quantum(t, state.quantum, state.classical)
    recast!(dx, dstate)
    QO_CHECKS[] && check_schroedinger(state.quantum, H[1])
    dschroedinger(state.quantum, H[1], dstate.quantum)
    recast!(dstate, dx)
end
function dschroedinger_stochastic(dx::Matrix,
        t, state::S, fstoch_quantum::Function,
        fstoch_classical::Nothing, dstate::S, n::Int) where {B<:Basis,T<:Ket{B},S<:State{B,T}}
    H = fstoch_quantum(t, state.quantum, state.classical)
    for i=1:n
        dx_i = @view dx[:, i]
        recast!(dx_i, dstate)
        QO_CHECKS[] && check_schroedinger(state.quantum, H[i])
        dschroedinger(state.quantum, H[i], dstate.quantum)
        recast!(dstate, dx_i)
    end
end
function dschroedinger_stochastic(dx::DiffArray, t,
            state::S, fstoch_quantum::Nothing, fstoch_classical::Function,
            dstate::S, ::Int) where {B<:Basis,T<:Ket{B},S<:State{B,T}}
    dclassical = @view dx[length(state.quantum)+1:end, :]
    fstoch_classical(t, state.quantum, state.classical, dclassical)
end
function dschroedinger_stochastic(dx::Matrix, t, state::S, fstoch_quantum::Function,
            fstoch_classical::Function, dstate::S, n::Int) where {B<:Basis,T<:Ket{B},S<:State{B,T}}
    dschroedinger_stochastic(dx, t, state, fstoch_quantum, nothing, dstate, n)

    dx_i = @view dx[length(state.quantum)+1:end, n+1:end]
    fstoch_classical(t, state.quantum, state.classical, dx_i)
end

function dmaster_stoch_dynamic(dx::Vector, t,
            state::S, fstoch_quantum::Function,
            fstoch_classical::Nothing, dstate::S, ::Int) where {B<:Basis,T<:Operator{B,B},S<:State{B,T}}
    result = fstoch_quantum(t, state.quantum, state.classical)
    QO_CHECKS[] && @assert length(result) == 2
    C, Cdagger = result
    QO_CHECKS[] && check_master_stoch(state.quantum, C, Cdagger)
    recast!(dx, dstate)
    QuantumOpticsBase.mul!(dstate.quantum,C[1],state.quantum)
    QuantumOpticsBase.mul!(dstate.quantum,state.quantum,Cdagger[1],true,true)
    dstate.quantum.data .-= tr(dstate.quantum)*state.quantum.data
    recast!(dstate, dx)
end
function dmaster_stoch_dynamic(dx::Matrix, t,
            state::S, fstoch_quantum::Function,
            fstoch_classical::Nothing, dstate::S, n::Int) where {B<:Basis,T<:Operator{B,B},S<:State{B,T}}
    result = fstoch_quantum(t, state.quantum, state.classical)
    QO_CHECKS[] && @assert length(result) == 2
    C, Cdagger = result
    QO_CHECKS[] && check_master_stoch(state.quantum, C, Cdagger)
    for i=1:n
        dx_i = @view dx[:, i]
        recast!(dx_i, dstate)
        QuantumOpticsBase.mul!(dstate.quantum,C[i],state.quantum)
        QuantumOpticsBase.mul!(dstate.quantum,state.quantum,Cdagger[i],true,true)
        dstate.quantum.data .-= tr(dstate.quantum)*state.quantum.data
        recast!(dstate, dx_i)
    end
end
function dmaster_stoch_dynamic(dx::DiffArray, t,
            state::S, fstoch_quantum::Nothing,
            fstoch_classical::Function, dstate::S, ::Int) where {B<:Basis,T<:Operator{B,B},S<:State{B,T}}
    dclassical = @view dx[length(state.quantum)+1:end, :]
    fstoch_classical(t, state.quantum, state.classical, dclassical)
end
function dmaster_stoch_dynamic(dx::Matrix, t,
            state::S, fstoch_quantum::Function,
            fstoch_classical::Function, dstate::S, n::Int) where {B<:Basis,T<:Operator{B,B},S<:State{B,T}}
    dmaster_stoch_dynamic(dx, t, state, fstoch_quantum, nothing, dstate, n)

    dx_i = @view dx[length(state.quantum)+1:end, n+1:end]
    fstoch_classical(t, state.quantum, state.classical, dx_i)
end

function integrate_master_stoch(tspan, df::Function, dg::Function,
                        rho0::State{B,T}, fout::Union{Nothing, Function},
                        n::Int;
                        kwargs...) where {B<:Basis,T<:Operator{B,B}}
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    x0 = Vector{eltype(rho0)}(undef, length(rho0))
    recast!(rho0, x0)
    state = copy(rho0)
    dstate = copy(rho0)
    integrate_stoch(tspan_, df, dg, x0, state, dstate, fout, n; kwargs...)
end

function recast!(state::State, x::SubArray)
    N = length(state.quantum)
    copyto!(x, 1, state.quantum.data, 1, N)
    copyto!(x, N+1, state.classical, 1, length(state.classical))
    x
end
function recast!(x::SubArray, state::State)
    N = length(state.quantum)
    copyto!(state.quantum.data, 1, x, 1, N)
    copyto!(state.classical, 1, x, N+1, length(state.classical))
end


text:
module ScoreAD

using ForwardDiff: Dual, value, partials

"""
    score_AD(p)

When `p` is real, return `1` (with the same type).

When `p` is a primitive for automatic differentiation, return ``log ∂p``, ie the
score function.

This function can be used for gradient estimation using simulation. See the
references in the README.
"""
score_AD(p::Real) = one(p)

score_AD(p::Dual{T,V}) where {T,V} = Dual{T}(one(V), partials(log(p)))

"""
    score_AD_log(ℓ)

Similar to [`score_AD`](@ref), but uses the log likelihood, ie
```julia
score_AD(p) == score_AD_log(log(p))
```
"""
score_AD_log(ℓ::Real) = one(ℓ)

score_AD_log(ℓ::Dual{T,V}) where {T,V} = Dual{T}(one(V), partials(ℓ))

"""
    reject_nonfinite(ℓ)

When the likelihood `ℓ` (or its derivative, when applicable) is not finite,
replace the value with `-Inf`.

Always return the same type as `ℓ`.

# Usage

This useful when calculations using automatic differentiation of a score
function in a (finite) sample, which may deliver non-finite derivatives for
extreme values. In this case, reject the parameter.
"""
@inline reject_nonfinite(ℓ::V) where {V <: Real} =
    ifelse(isfinite(ℓ), ℓ, V(-Inf))

@inline function reject_nonfinite(ℓ::Dual{T,V}) where {T,V}
    p = partials(ℓ)
    if isfinite(value(ℓ)) && all(isfinite, p)
        ℓ
    else
        Dual{T}(V(-Inf), p)
    end
end

end # module


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7, x8, x9, x10)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9", "x10"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7, x8, x9, x10]);F = [8980*x1^2*x2^2+3960*x1^2*x2*x4+264*x1^2*x2*x5+900*x1^2*x4^2+120*x1^2*x4*x5+4*x1^2*x5^2+3600*x1^2*x7^2-5280*x1^2*x7*x9+1936*x1^2*x9^2+2772*x1*x10*x2*x5-10880*x1*x10*x2*x9+1260*x1*x10*x4*x5+84*x1*x10*x5^2-10336*x1*x2^2*x9-7616*x1*x2*x4*x9+2992*x1*x2*x5*x6-8052*x1*x2*x7*x8-3944*x1*x2*x7*x9-11520*x1*x3*x7^2+8448*x1*x3*x7*x9-3660*x1*x4*x7*x8-244*x1*x5*x7*x8-11520*x1*x5*x7*x9+8448*x1*x5*x9^2+1560*x1*x7^2*x9-1144*x1*x7*x9^2+1681*x10^4+410*x10^2*x2*x7-4264*x10^2*x3*x8+441*x10^2*x5^2+5909*x10^2*x8^2+6400*x10^2*x9^2-13110*x10*x2*x6*x8-2760*x10*x2*x7*x8+12160*x10*x2*x9^2-7038*x10*x3*x6*x8+8960*x10*x4*x9^2-3520*x10*x5*x6*x9-2562*x10*x5*x7*x8-10212*x10*x6^2*x8+4640*x10*x7*x9^2+9025*x2^2*x6^2+3800*x2^2*x6*x7+425*x2^2*x7^2+5776*x2^2*x9^2+9690*x2*x3*x6^2+2040*x2*x3*x6*x7-520*x2*x3*x7*x8+8512*x2*x4*x9^2-3344*x2*x5*x6*x9+14060*x2*x6^3+2960*x2*x6^2*x7+140*x2*x7*x8^2+4408*x2*x7*x9^2+2601*x3^2*x6^2+9216*x3^2*x7^2+2704*x3^2*x8^2+18432*x3*x5*x7*x9+7548*x3*x6^3-2496*x3*x7^2*x9-1456*x3*x8^3+3136*x4^2*x9^2-2464*x4*x5*x6*x9+3248*x4*x7*x9^2+484*x5^2*x6^2+9216*x5^2*x9^2-1276*x5*x6*x7*x9-2496*x5*x7*x9^2+5476*x6^4+3721*x7^2*x8^2+1010*x7^2*x9^2+196*x8^4-1886*x1*x10^2-230*x1*x2*x7+9636*x1*x2*x9+2392*x1*x3*x8+4380*x1*x4*x9+4800*x1*x5*x7-3228*x1*x5*x9-644*x1*x8^2-738*x10^2*x7+12420*x10*x5*x8+3066*x10*x5*x9-17100*x2*x5*x6-3600*x2*x5*x7-90*x2*x7^2-9180*x3*x5*x6-7680*x3*x5*x7+936*x3*x7*x8-7680*x5^2*x9-13320*x5*x6^2+1040*x5*x7*x9-252*x7*x8^2-8906*x7*x8*x9+529*x1^2+414*x1*x7+9700*x5^2+81*x7^2+5329*x9^2, 17960*x1^2*x2+3960*x1^2*x4+264*x1^2*x5+2772*x1*x10*x5-10880*x1*x10*x9-20672*x1*x2*x9-7616*x1*x4*x9+2992*x1*x5*x6-8052*x1*x7*x8-3944*x1*x7*x9+410*x10^2*x7-13110*x10*x6*x8-2760*x10*x7*x8+12160*x10*x9^2+18050*x2*x6^2+7600*x2*x6*x7+850*x2*x7^2+11552*x2*x9^2+9690*x3*x6^2+2040*x3*x6*x7-520*x3*x7*x8+8512*x4*x9^2-3344*x5*x6*x9+14060*x6^3+2960*x6^2*x7+140*x7*x8^2+4408*x7*x9^2-230*x1*x7+9636*x1*x9-17100*x5*x6-3600*x5*x7-90*x7^2, -11520*x1*x7^2+8448*x1*x7*x9-4264*x10^2*x8-7038*x10*x6*x8+9690*x2*x6^2+2040*x2*x6*x7-520*x2*x7*x8+5202*x3*x6^2+18432*x3*x7^2+5408*x3*x8^2+18432*x5*x7*x9+7548*x6^3-2496*x7^2*x9-1456*x8^3+2392*x1*x8-9180*x5*x6-7680*x5*x7+936*x7*x8, 3960*x1^2*x2+1800*x1^2*x4+120*x1^2*x5+1260*x1*x10*x5-7616*x1*x2*x9-3660*x1*x7*x8+8960*x10*x9^2+8512*x2*x9^2+6272*x4*x9^2-2464*x5*x6*x9+3248*x7*x9^2+4380*x1*x9, 264*x1^2*x2+120*x1^2*x4+8*x1^2*x5+2772*x1*x10*x2+1260*x1*x10*x4+168*x1*x10*x5+2992*x1*x2*x6-244*x1*x7*x8-11520*x1*x7*x9+8448*x1*x9^2+882*x10^2*x5-3520*x10*x6*x9-2562*x10*x7*x8-3344*x2*x6*x9+18432*x3*x7*x9-2464*x4*x6*x9+968*x5*x6^2+18432*x5*x9^2-1276*x6*x7*x9-2496*x7*x9^2+4800*x1*x7-3228*x1*x9+12420*x10*x8+3066*x10*x9-17100*x2*x6-3600*x2*x7-9180*x3*x6-7680*x3*x7-15360*x5*x9-13320*x6^2+1040*x7*x9+19400*x5, 2992*x1*x2*x5-13110*x10*x2*x8-7038*x10*x3*x8-3520*x10*x5*x9-20424*x10*x6*x8+18050*x2^2*x6+3800*x2^2*x7+19380*x2*x3*x6+2040*x2*x3*x7-3344*x2*x5*x9+42180*x2*x6^2+5920*x2*x6*x7+5202*x3^2*x6+22644*x3*x6^2-2464*x4*x5*x9+968*x5^2*x6-1276*x5*x7*x9+21904*x6^3-17100*x2*x5-9180*x3*x5-26640*x5*x6, 7200*x1^2*x7-5280*x1^2*x9-8052*x1*x2*x8-3944*x1*x2*x9-23040*x1*x3*x7+8448*x1*x3*x9-3660*x1*x4*x8-244*x1*x5*x8-11520*x1*x5*x9+3120*x1*x7*x9-1144*x1*x9^2+410*x10^2*x2-2760*x10*x2*x8-2562*x10*x5*x8+4640*x10*x9^2+3800*x2^2*x6+850*x2^2*x7+2040*x2*x3*x6-520*x2*x3*x8+2960*x2*x6^2+140*x2*x8^2+4408*x2*x9^2+18432*x3^2*x7+18432*x3*x5*x9-4992*x3*x7*x9+3248*x4*x9^2-1276*x5*x6*x9-2496*x5*x9^2+7442*x7*x8^2+2020*x7*x9^2-230*x1*x2+4800*x1*x5-738*x10^2-3600*x2*x5-180*x2*x7-7680*x3*x5+936*x3*x8+1040*x5*x9-252*x8^2-8906*x8*x9+414*x1+162*x7, -8052*x1*x2*x7-3660*x1*x4*x7-244*x1*x5*x7-4264*x10^2*x3+11818*x10^2*x8-13110*x10*x2*x6-2760*x10*x2*x7-7038*x10*x3*x6-2562*x10*x5*x7-10212*x10*x6^2-520*x2*x3*x7+280*x2*x7*x8+5408*x3^2*x8-4368*x3*x8^2+7442*x7^2*x8+784*x8^3+2392*x1*x3-1288*x1*x8+12420*x10*x5+936*x3*x7-504*x7*x8-8906*x7*x9, -5280*x1^2*x7+3872*x1^2*x9-10880*x1*x10*x2-10336*x1*x2^2-7616*x1*x2*x4-3944*x1*x2*x7+8448*x1*x3*x7-11520*x1*x5*x7+16896*x1*x5*x9+1560*x1*x7^2-2288*x1*x7*x9+12800*x10^2*x9+24320*x10*x2*x9+17920*x10*x4*x9-3520*x10*x5*x6+9280*x10*x7*x9+11552*x2^2*x9+17024*x2*x4*x9-3344*x2*x5*x6+8816*x2*x7*x9+18432*x3*x5*x7-2496*x3*x7^2+6272*x4^2*x9-2464*x4*x5*x6+6496*x4*x7*x9+18432*x5^2*x9-1276*x5*x6*x7-4992*x5*x7*x9+2020*x7^2*x9+9636*x1*x2+4380*x1*x4-3228*x1*x5+3066*x10*x5-7680*x5^2+1040*x5*x7-8906*x7*x8+10658*x9, 2772*x1*x2*x5-10880*x1*x2*x9+1260*x1*x4*x5+84*x1*x5^2+6724*x10^3+820*x10*x2*x7-8528*x10*x3*x8+882*x10*x5^2+11818*x10*x8^2+12800*x10*x9^2-13110*x2*x6*x8-2760*x2*x7*x8+12160*x2*x9^2-7038*x3*x6*x8+8960*x4*x9^2-3520*x5*x6*x9-2562*x5*x7*x8-10212*x6^2*x8+4640*x7*x9^2-3772*x1*x10-1476*x10*x7+12420*x5*x8+3066*x5*x9];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
module infopunk_tools

function cat(x)
	print(read(x, String))
end

end # module


text:
export loaddata_analyze, savedata_analyze, read_analyze_data, read_analyze_hdr

#####################
#  Analyze export
#####################

baremodule ElementType
  FLOAT32 = 16
  INT16 = 4
  RGBA = 2304
  RGB = 128
end

function _bitpix(elemType)
  if elemType == ElementType.FLOAT32
    return 32
  elseif elemType == ElementType.INT16
    return 16
  elseif elemType == ElementType.RGBA
    return 32
  elseif elemType == ElementType.RGB
    return 24
  else
    error("elemType ", elemType, " not yet supported")
  end
end

function savedata_analyze(filename::AbstractString, image::ImageMeta{T};
                          kargs...) where {T<:ColorTypes.TransparentColor}
  rgbimage = convert(RGB{N0f8},image)
  savedata_analyze{T<:TransparentColor}(filename, rgbimage; kargs...)
end

function savedata_analyze(filename::AbstractString, im::ImageMeta; verbose = false, saveas16bit=false, perm=nothing, permRGBData=false)
  filenamebase, ext = splitext(filename)

  im[:permRGBData] = permRGBData

  if get(im,"dim",3) < 3
    warn("voxel size will be wrong for 1D/2D data")
  end
  T = eltype(im)
  if T <: AbstractFloat
    c = arraydata(im)
    im[:cmin] = minimum(c)
    im[:cmax] = maximum(c)
    if saveas16bit
        im[:scalingFactor] = 2^14 / maximum(c)
        c *= im[:scalingFactor]
        c = round(Int16,c)
        im[:datatype] = ElementType.INT16
    else
        im[:scalingFactor] = 1
        if eltype(c) != Float32
          #c = map(Float32,c) #this does not work for franzi
          c = AxisArray(map(f, c.data), c.axes...)
        end
        im[:datatype] = ElementType.FLOAT32
    end
    im[:size] = [size(c,i) for i=1:4]
    im[:ndims] = ndims(c)
  else
      c = arraydata(im).data
      im[:cmin] = 0
      im[:cmax] = 1
      im[:scalingFactor] = 1
    if T <: RGB
      # This is what is returned by colorize or blend
      im[:datatype] = ElementType.RGB
      im[:size] = [size(c,i) for i=1:4]
      im[:ndims] = ndims(c)

    elseif T == UInt8
      # This is what is returned by fusedVolumes
      im[:datatype] = ElementType.RGB
      c = reorderColor(c)
      im[:size] = [size(c,i) for i=1:4]
      im[:ndims] = ndims(c)

    else
      error("Colored data has to be passed as ...")
    end
    if permRGBData

      if ndims(c) == 3
        c = reinterpret(UInt8, c, (3, size(c,1), size(c,2), size(c,3)))
        c = permutedims(c,[2,3,1,4])
      else
        c = reinterpret(UInt8, c, (3, size(c,1), size(c,2), size(c,3), size(c,4)))
        c = permutedims(c,[2,3,1,4,5])
      end
    end
  end

  header = properties(im)
  header[:pixelspacing] = collect(ustrip.(pixelspacing(im)))
  header[:offset] = collect(imcenter(im))

  if ext == ".nii"
    im[:nii_file] = true
    im[:vox_offset] = 352

    open(filenamebase*".nii","w") do fd
      write_analyze_hdr(fd, header)
      write(fd, Int32(0))
      write(fd, c)
    end
  else
    im[:nii_file] = false
    im[:vox_offset] = 0
    # write data
    open(filenamebase*".img","w") do fd
      write(fd, c)
    end

    # write Header
    open(filenamebase*".hdr","w") do fd
      write_analyze_hdr(fd, header)
    end
  end

end

function loaddata_analyze(filename::AbstractString)
  filenamebase, ext = splitext(filename)

  header = nothing
  c = nothing

  if ext == ".nii"
    c = open(filenamebase*".nii","r") do fd
      header = read_analyze_hdr(fd)
      read(fd,Int32)
      read_analyze_data(fd,header)
    end
  else
    open(filenamebase*".hdr","r") do fd
      header = read_analyze_hdr(fd)
    end

    c = open(filenamebase*".img","r") do fd
      read_analyze_data(fd,header)
    end
  end

  pixspacing = header[:pixelspacing]*1000u"mm"
  offset = header[:offset]*1000u"mm" .- 0.5.*pixspacing.*header[:size][1:3] .+ 0.5.*pixspacing
  if ndims(c) == 4
    im = AxisArray(c, (:x,:y,:z,:time), tuple(pixspacing...,one(typeof(header[:pixelspacing][1]))*u"s"),
                                   tuple(offset..., zero(typeof(header[:offset][1]))*u"s"))
  else
    im = AxisArray(c, (:x,:y,:z), tuple(pixspacing...), tuple(offset... ))
  end

  delete!(header, "pixelspacing")
  delete!(header, "timedim")

  imMeta = ImageMeta(im, header)

  return imMeta
end

function read_analyze_data(fd,header)
  if header[:datatype] == ElementType.RGB
    # colored data
    s = header[:size]
    if header[:permRGBData] == false
      c = read!(fd, Array{RGB{N0f8}}(undef, header[:size]...))
    else
      if header[:dim] == 3
        @info "3D data"
        c = read!(fd, Array{Int8}(undef, s[1], s[2], 3, s[3]))
	#c = read(fd, UInt8, s[1], s[2], 3, s[3])
        c = permutedims(c,[3,1,2,4])
        c = reinterpret(RGB{N0f8}, c, (size(c,2),size(c,3),size(c,4)))
      else
        @info "4D data"
        c = read!(fd, Array{Int8}(undef, s[1], s[2], 3, s[3], s[4]))
        #c = read(fd, UInt8, s[1], s[2], 3, s[3], s[4])
        c = permutedims(c,[3,1,2,4,5])
        c = reinterpret(RGB{N0f8}, c, (size(c,2),size(c,3),size(c,4),size(c,5)))
      end
    end
  else
    if header[:datatype] == ElementType.FLOAT32
      c = read!(fd, Array{Float32}(undef, header[:size]...))
    else
      c = map(Float32, read!(fd, Array{Int16}(undef, header[:size]...)) )
    end

    if header[:scalingFactor] > 0
      c *= (1/header[:scalingFactor])
    end
  end
  return c
end

function write_analyze_hdr(fd, header)

    s = header[:pixelspacing].*1000

    write(fd, Int32(348)) #sizeof_hdr::Int32

    data_type = "MPI"
    write(fd, data_type)
    write(fd, zeros(Int8,10-length(data_type))) #char data_type[10];

    write(fd, zeros(Int8,18)) #char db_name[18];
    write(fd, Int32(16384)) #  int extents;
    write(fd, Int16(0)) #short int session_error;
    write(fd, 'r') #char regular;
    write(fd, '0') #char hkey_un0;

    write(fd, Int16(header[:ndims])) #short int dim[8];
    write(fd, Int16(header[:size][1]))
    write(fd, Int16(header[:size][2]))
    write(fd, Int16(header[:size][3]))
    write(fd, Int16(header[:size][4]))
    write(fd, ones(Int16,3))

    write(fd, "mm  ") #vox_units	uchar[4]		spatial units of measure for a voxel
    write(fd, "mmol(Fe)") # cal_units	uchar[8]	name of the calibration unit

    write(fd, zeros(Int16,1)) # unused1	Int16

    write(fd, Int16(header[:datatype])) #short int datatype; -> Float32
    write(fd, Int16(_bitpix(header[:datatype]))) #short int bitpix;
    write(fd, Int16(0)) #short int dim_un0;

    write(fd, Float32(0)) # number of spatial dims
    write(fd, Float32(s[1]))
    write(fd, Float32(s[2]))
    write(fd, Float32(s[3]))
    write(fd, zeros(Float32,4)) # float pixdim[8];

    write(fd, Float32(header[:vox_offset])) #float vox_offset;
    #float funused1;
    #float funused2;

    off = ustrip.(uconvert.(Unitful.mm,header[:offset]))
    write(fd, Float32(off[1]))
    write(fd, Float32(off[2]))
    write(fd, Float32(off[3])) #float funused3;
    write(fd, Float32(header[:scalingFactor])) #float cal_max;  HERE WE STORE THE SCALING FACTOR!
    write(fd, Float32(0)) #float cal_min;
    write(fd, Float32(header[:permRGBData])) #float compressed;
    write(fd, Float32(0)) #float verified;
    #write(fd, round(Int32,header[:cmin]))
    #write(fd, round(Int32,header[:cmax]))#int glmax,glmin;
    write(fd, zero(Int32))
    write(fd, Int32(typemax(Int32)))#int glmax,glmin;

    description = Vector{UInt8}(header[:experimentDescription])
    description = length(description) > 9 ? description[1:9] : description
    @debug "" description
    write(fd, description)
    write(fd, zeros(Int8,80-length(description))) #char descrip[80];

    write(fd, zeros(Int8,24)) #char aux_file[24];
    write(fd, Int8(0)) #char orient;
    write(fd, zeros(Int8,10)) #char originator[10];
    write(fd, zeros(Int8,10)) #char generated[10];

    scannum = Vector{UInt8}(string(header[:experimentNumber]))
    scannum = length(scannum) > 9 ? scannum[1:9] : scannum
    @debug "" scannum
    write(fd, scannum)
    write(fd, zeros(Int8,10-length(scannum))) #char scannum[10];

    subjName = Vector{UInt8}(header[:experimentSubject])
    subjName = length(subjName) > 9 ? subjName[1:9] : subjName
    @debug "" subjName
    write(fd, subjName)
    write(fd, zeros(Int8,10-length(subjName))) #char patient_id[10];

    dateStr, timeStr = split("$(header[:acqStartTime])","T")
    dateStr = prod(split(dateStr,"-"))
    timeStr = split(timeStr,".")[1]
    timeStr = prod(split(timeStr,":"))

    date = Vector{UInt8}(dateStr)
    date = length(date) > 9 ? date[1:9] : date
    write(fd, date)
    write(fd, zeros(Int8,10-length(date))) #char exp_date[10];

    time = Vector{UInt8}(timeStr)
    time = length(time) > 9 ? time[1:9] : time
    write(fd, time)
    write(fd, zeros(Int8,10-sizeof(time))) #char exp_time[10];

    write(fd, zeros(Int8,3)) #char hist_un0[3];
    write(fd, Int32(0)) #int views
    write(fd, Int32(0)) #  int vols_added;
    write(fd, Int32(0)) #  int start_field;
    write(fd, Int32(0)) #  int field_skip;
    write(fd, Int32(0))
    write(fd, Int32(0)) #  int omax, omin;
    write(fd, Int32(0)) #  int smax, smin;
    if !header[:nii_file]
      write(fd, UInt32(0x6E693100))
    else
      write(fd, UInt32(0x6E2B3100))
    end
end

function read_analyze_hdr(fd)
    header = Dict{Symbol,Any}()

    read(fd, Int32) #sizeof_hdr::Int32

    data_type = String( read!(fd, Array{UInt8}(undef, 10)) ) #char data_type[10];
    header[:data_type] = data_type

    read!(fd, Array{Int8}(undef, 18)) #char db_name[18];
    read(fd, Int32) #  int extents;
    read(fd, Int16) #short int session_error;
    read(fd, Int8) #char regular;
    read(fd, Int8) #char hkey_un0;

    dim = read(fd, Int16)
    header[:dim] = dim
    sizes = read!(fd, Array{Int16}(undef, 7)) #short int dim[8];
    header[:size] = sizes[1:dim]

    read!(fd, Array{Int8}(undef, 4)) #vox_units	uchar[4]		spatial units of measure for a voxel
    read!(fd, Array{Int8}(undef, 8)) # cal_units	uchar[8]									name of the calibration unit

    read(fd, Int16) # unused1	Int16

    header[:datatype] = read(fd, Int16) #short int datatype; -> Float32
    read(fd, Int16) #short int bitpix;
    read(fd, Int16) #short int dim_un0;

    read(fd, Float32)
    voxelSize = read!(fd, Array{Float32}(undef, 7)) # float pixdim[8];
    header[:pixelspacing] = voxelSize[1:3] ./ 1000 # this is restricted to 3D...

    if header[:dim] == 4
      header[:timedim] = 4
    else
      header[:timedim] = 0 #This is likely not always correct
    end

    read(fd, Float32) #float vox_offset;
    #float funused1;
    #float funused2;
    off = read!(fd, Array{Float32}(undef,3)) #float funused3;
    header[:offset] = off
    if header[:data_type][1:3]!="MPI"
      header[:offset][:] .= 0.0
    end
    header[:scalingFactor] = read(fd, Float32) #float cal_max; #HERE WE READ THE SCALING FACTOR
    read(fd, Float32) #float cal_min;
    header[:permRGBData] = (read(fd, Float32) != 0.0) #float compressed;
    read(fd, Float32) #float verified;
    cmin = read(fd, Int32)
    header[:cmin] = cmin
    cmax = read(fd, Int32)
    header[:cmax] = cmax #int glmax,glmin;

    descr = String( read!(fd, Array{UInt8}(undef, 80)) ) #char descrip[80];
    header[:experimentDescription] = descr

    read!(fd, Array{UInt8}(undef, 24)) #char aux_file[24];
    read(fd, UInt8) #char orient;
    read!(fd, Array{UInt8}(undef, 10)) #char originator[10];
    read!(fd, Array{UInt8}(undef, 10)) #char generated[10];

    scannum =  String( read!(fd, Array{UInt8}(undef, 10)) ) #char scannum[10];
    header[:experimentNumber] = scannum

    patient_id =  String( read!(fd, Array{UInt8}(undef, 10)) ) #char patient_id[10];
    header[:experimentSubject] = patient_id

    exp_date =  String( read!(fd, Array{UInt8}(undef, 10)) ) #char exp_date[10];
    header[:date] = exp_date

    exp_time =  String( read!(fd, Array{UInt8}(undef, 10)) ) #char exp_time[10];
    header[:time] = exp_time

    read!(fd, Array{UInt8}(undef, 3)) #char hist_un0[3];
    read(fd, Int32) #int views
    read(fd, Int32) #  int vols_added;
    read(fd, Int32) #  int start_field;
    read(fd, Int32) #  int field_skip;
    read(fd, Int32)
    read(fd, Int32) #  int omax, omin;
    read(fd, Int32)
    read(fd, Int32) #  int smax, smin;

  return header
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# extensions to Core types to add features in Base

# hook up VecElement constructor to Base.convert
VecElement{T}(arg) where {T} = VecElement{T}(convert(T, arg))
convert(::Type{T}, arg)  where {T<:VecElement} = T(arg)
convert(::Type{T}, arg::T) where {T<:VecElement} = arg

# ## dims-type-converting Array constructors for convenience
# type and dimensionality specified, accepting dims as series of Integers
Vector{T}(::UndefInitializer, m::Integer) where {T} = Vector{T}(undef, Int(m))
Matrix{T}(::UndefInitializer, m::Integer, n::Integer) where {T} = Matrix{T}(undef, Int(m), Int(n))
Array{T,N}(::UndefInitializer, d::Vararg{Integer,N}) where {T,N} = Array{T,N}(undef, convert(Tuple{Vararg{Int}}, d))
# type but not dimensionality specified, accepting dims as series of Integers
Array{T}(::UndefInitializer, m::Integer) where {T} = Array{T,1}(undef, Int(m))
Array{T}(::UndefInitializer, m::Integer, n::Integer) where {T} = Array{T,2}(undef, Int(m), Int(n))
Array{T}(::UndefInitializer, m::Integer, n::Integer, o::Integer) where {T} = Array{T,3}(undef, Int(m), Int(n), Int(o))
Array{T}(::UndefInitializer, d::Integer...) where {T} = Array{T}(undef, convert(Tuple{Vararg{Int}}, d))
# dimensionality but not type specified, accepting dims as series of Integers
Vector(::UndefInitializer, m::Integer) = Vector{Any}(undef, Int(m))
Matrix(::UndefInitializer, m::Integer, n::Integer) = Matrix{Any}(undef, Int(m), Int(n))
# Dimensions as a single tuple
Array{T}(::UndefInitializer, d::NTuple{N,Integer}) where {T,N} = Array{T,N}(undef, convert(Tuple{Vararg{Int}}, d))
Array{T,N}(::UndefInitializer, d::NTuple{N,Integer}) where {T,N} = Array{T,N}(undef, convert(Tuple{Vararg{Int}}, d))
# empty vector constructor
Vector() = Vector{Any}(undef, 0)

# Array constructors for nothing and missing
# type and dimensionality specified
Array{T,N}(::Nothing, d...) where {T,N} = fill!(Array{T,N}(undef, d...), nothing)
Array{T,N}(::Missing, d...) where {T,N} = fill!(Array{T,N}(undef, d...), missing)
# type but not dimensionality specified
Array{T}(::Nothing, d...) where {T} = fill!(Array{T}(undef, d...), nothing)
Array{T}(::Missing, d...) where {T} = fill!(Array{T}(undef, d...), missing)

max_stars_repo_name: EricForgy/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
using SynapticFilter
using BenchmarkTools, Test

const RUN_BENCHMARKS = false # optional intermediate benchmarks

@testset "SynapticFilter.jl" begin
    include("models.jl")
    include("filters.jl")
    include("performance.jl")
end


text:
module BismarkSummary

# package code goes here
export parse_bismark_report
export get_bismark_report_filenames_from_plite_pipeline
export parse_bismark_reports
export add_report_info_to_datasource!
export make_bismark_report
export bismark_report
include(Pkg.dir("BismarkSummary","src","bismark-report.jl"))

export make_coverage_stats_table
include(Pkg.dir("BismarkSummary","src","bismark-coverage.jl"))
end # module


text:
const seed = 420  # Random seed to use for all pseudorandom number generators.

datatuple(A) = NamedTuple{propertynames(A)}(Array(data(a)) for a in A)

const T₀ = 9.85
const S₀ = 35.0

function get_output_tuple(output, iter, tuplename)
    file = jldopen(output.filepath, "r")
    output_tuple = file["timeseries/$tuplename/$iter"]
    close(file)
    return output_tuple
end

function run_thermal_bubble_regression_tests(arch)
    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 100, 100, 100
    Δt = 6

    model = BasicModel(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), architecture=arch, ν=4e-2, κ=4e-2,
                       coriolis=FPlane(f=1e-4))

    model.tracers.T.data.parent .= T₀
    model.tracers.S.data.parent .= S₀

    # Add a cube-shaped warm temperature anomaly that takes up the middle 50%
    # of the domain volume.
    i1, i2 = round(Int, Nx/4), round(Int, 3Nx/4)
    j1, j2 = round(Int, Ny/4), round(Int, 3Ny/4)
    k1, k2 = round(Int, Nz/4), round(Int, 3Nz/4)
    model.tracers.T.data[i1:i2, j1:j2, k1:k2] .+= 0.01

    nc_writer = NetCDFOutputWriter(dir=".",
                                   prefix="thermal_bubble_regression_",
                                   frequency=10, padding=2)

    # Uncomment to include a NetCDF output writer that produces the regression.
    # push!(model.output_writers, nc_writer)

    time_step!(model, 10, Δt)

    u = read_output(nc_writer, "u", 10)
    v = read_output(nc_writer, "v", 10)
    w = read_output(nc_writer, "w", 10)
    T = read_output(nc_writer, "T", 10)
    S = read_output(nc_writer, "S", 10)

    field_names = ["u", "v", "w", "T", "S"]
    fields = [model.velocities.u, model.velocities.v, model.velocities.w, model.tracers.T, model.tracers.S]
    fields_gm = [u, v, w, T, S]
    for (field_name, φ, φ_gm) in zip(field_names, fields, fields_gm)
        φ_min = minimum(Array(data(φ)) - φ_gm)
        φ_max = maximum(Array(data(φ)) - φ_gm)
        φ_mean = mean(Array(data(φ)) - φ_gm)
        φ_abs_mean = mean(abs.(Array(data(φ)) - φ_gm))
        φ_std = std(Array(data(φ)) - φ_gm)
        @info(@sprintf("Δ%s: min=%.6g, max=%.6g, mean=%.6g, absmean=%.6g, std=%.6g\n",
                       field_name, φ_min, φ_max, φ_mean, φ_abs_mean, φ_std))
    end

    # Now test that the model state matches the regression output.
    @test all(Array(data(model.velocities.u)) .≈ u)
    @test all(Array(data(model.velocities.v)) .≈ v)
    @test all(Array(data(model.velocities.w)) .≈ w)
    @test all(Array(data(model.tracers.T))    .≈ T)
    @test all(Array(data(model.tracers.S))    .≈ S)
end

function run_rayleigh_benard_regression_test(arch)

    #####
    ##### Parameters
    #####
          α = 2                 # aspect ratio
          n = 1                 # resolution multiple
         Ra = 1e6               # Rayleigh number
    Nx = Ny = 8n * α            # horizontal resolution
    Lx = Ly = 1.0 * α           # horizontal extent
         Nz = 16n               # vertical resolution
         Lz = 1.0               # vertical extent
         Pr = 0.7               # Prandtl number
          a = 1e-1              # noise amplitude for initial condition
         Δb = 1.0               # buoyancy differential

    # Rayleigh and Prandtl determine transport coefficients
    ν = sqrt(Δb * Pr * Lz^3 / Ra)
    κ = ν / Pr

    #####
    ##### Model setup
    #####

    # Force salinity as a passive tracer (βS=0)
    S★(x, z) = exp(4z) * sin(2π/Lx * x)
    FS(i, j, k, grid, time, U, Φ, params) = 1/10 * (S★(grid.xC[i], grid.zC[k]) - Φ.S[i, j, k])

    model = Model(
               architecture = arch,
                       grid = RegularCartesianGrid(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz)),
                    closure = ConstantIsotropicDiffusivity(ν=ν, κ=κ),
                   buoyancy = BuoyancyTracer(),
        boundary_conditions = BoundaryConditions(T=HorizontallyPeriodicBCs(
                                top=BoundaryCondition(Value, 0.0), bottom=BoundaryCondition(Value, Δb))),
                    forcing = Forcing(FS=FS)
    )

    ArrayType = typeof(model.velocities.u.data.parent)  # The type of the underlying data, not the offset array.
    Δt = 0.01 * min(model.grid.Δx, model.grid.Δy, model.grid.Δz)^2 / ν

    spinup_steps = 1000
      test_steps = 100

    output_U(model) = datatuple(model.velocities)
    output_Φ(model) = datatuple(model.tracers)
    output_G(model) = datatuple(model.timestepper.Gⁿ)
    outputfields = Dict(:U=>output_U, :Φ=>output_Φ, :G=>output_G)

    prefix = "data_rayleigh_benard_regression"
    outputwriter = JLD2OutputWriter(model, outputfields; dir=".", prefix=prefix,
                                    frequency=test_steps, including=[])

    #####
    ##### Initial condition and spinup steps for creating regression test data
    #####

    #=
    @warn ("Generating new data for the Rayleigh-Benard regression test.
           New regression test data generation will fail unless the JLD2
           file $prefix.jld2 is manually deleted.")

    ξ(z) = a * rand() * z * (Lz + z) # noise, damped at the walls
    b₀(x, y, z) = (ξ(z) - z) / Lz
    set_ic!(model, T=b₀)

    time_step!(model, spinup_steps-test_steps, Δt)
    push!(model.output_writers, outputwriter)
    time_step!(model, 2test_steps, Δt)
    =#

    #####
    ##### Regression test
    #####

    # Load initial state
    u₀, v₀, w₀ = get_output_tuple(outputwriter, spinup_steps, :U)
    T₀, S₀ = get_output_tuple(outputwriter, spinup_steps, :Φ)
    Gu, Gv, Gw, GT, GS = get_output_tuple(outputwriter, spinup_steps, :G)

    data(model.velocities.u) .= ArrayType(u₀)
    data(model.velocities.v) .= ArrayType(v₀)
    data(model.velocities.w) .= ArrayType(w₀)
    data(model.tracers.T)    .= ArrayType(T₀)
    data(model.tracers.S)    .= ArrayType(S₀)

    data(model.timestepper.Gⁿ.Gu) .= ArrayType(Gu)
    data(model.timestepper.Gⁿ.Gv) .= ArrayType(Gv)
    data(model.timestepper.Gⁿ.Gw) .= ArrayType(Gw)
    data(model.timestepper.Gⁿ.GT) .= ArrayType(GT)
    data(model.timestepper.Gⁿ.GS) .= ArrayType(GS)

    model.clock.iteration = spinup_steps
    model.clock.time = spinup_steps * Δt
    length(model.output_writers) > 0 && pop!(model.output_writers)

    # Step the model forward and perform the regression test
    time_step!(model, test_steps, Δt; init_with_euler=false)

    u₁, v₁, w₁ = get_output_tuple(outputwriter, spinup_steps+test_steps, :U)
    T₁, S₁ = get_output_tuple(outputwriter, spinup_steps+test_steps, :Φ)

    field_names = ["u", "v", "w", "T", "S"]
    fields = [model.velocities.u, model.velocities.v, model.velocities.w, model.tracers.T, model.tracers.S]
    fields_gm = [u₁, v₁, w₁, T₁, S₁]
    for (field_name, φ, φ_gm) in zip(field_names, fields, fields_gm)
        φ_min = minimum(Array(data(φ)) - φ_gm)
        φ_max = maximum(Array(data(φ)) - φ_gm)
        φ_mean = mean(Array(data(φ)) - φ_gm)
        φ_abs_mean = mean(abs.(Array(data(φ)) - φ_gm))
        φ_std = std(Array(data(φ)) - φ_gm)
        @info(@sprintf("Δ%s: min=%.6g, max=%.6g, mean=%.6g, absmean=%.6g, std=%.6g\n",
                       field_name, φ_min, φ_max, φ_mean, φ_abs_mean, φ_std))
    end

    # Now test that the model state matches the regression output.
    @test all(Array(data(model.velocities.u)) .≈ u₁)
    @test all(Array(data(model.velocities.v)) .≈ v₁)
    @test all(Array(data(model.velocities.w)) .≈ w₁)
    @test all(Array(data(model.tracers.T))    .≈ T₁)
    @test all(Array(data(model.tracers.S))    .≈ S₁)
    return nothing
end

@testset "Regression" begin
    println("Running regression tests...")

    for arch in archs
        @testset "Thermal bubble [$(typeof(arch))]" begin
            println("  Testing thermal bubble regression [$(typeof(arch))]")
            run_thermal_bubble_regression_tests(arch)
        end

        @testset "Rayleigh–Bénard tracer [$(typeof(arch))]" begin
            println("  Testing Rayleigh–Bénard tracer regression [$(typeof(arch))]")
            run_rayleigh_benard_regression_test(arch)
        end
    end
end


text:
include("../utils/format/Heading.jl")

# float32/64
float_64= 1e-2
float_32= 1f-2
println(typeof(float_64), " ", typeof(float_32))
#: Float64 Float32

# underscore as digit separator
println(100_000, " ", 0.00_00_5)
#: 100000 5.0e-5

# flotiing point has two zeros - positive and negative
# both are `=` same but `bitstring` different
println(0.0 == -0.0)
#: true
println(bitstring(0.0) == bitstring(-0.0))
#: false

Heading.main("Special floating point values")

println(-5/0)
println(-5/Inf)
#: -Inf
#: -0.0

println(Inf+500)
#: Inf

# as Inf+x (or) Inf-x (or) Inf + Inf all equal to Inf
# Information is lost. So, Inf - Inf cannot be exactly 0
println(Inf - Inf)
#: NaN

# epsilon: distance between 1.0 and the next larger
# representable floating-point value i.e smallest distance
# possible from 1.0
println(eps(Float16))
println(eps(Float32))
println(eps(Float64))
println(eps()) # same as 64
#: 0.000977
#: 1.1920929e-7
#: 2.220446049250313e-16
#: 2.220446049250313e-16

println(typemax(BigFloat))
#: Inf


text:
module FoldsCatwalk

export CatwalkEx, OptimizeInner

using Catwalk
using Transducers
using Transducers:
    Executor,
    R_,
    Transducer,
    complete,
    extract_transducer,
    inner,
    next,
    start,
    xform

"""
    CatwalkEx(batchsize)

JIT compile a fold for every `batchsize` iteration. See [`OptimizeInner`](@ref)
for how to specify what to JIT.
"""
struct CatwalkEx <: Executor
    batchsize::Int
end
CatwalkEx() = CatwalkEx(1_000_000)

Transducers.maybe_set_simd(exc::CatwalkEx, _) = exc

"""
    OptimizeInner([factory = Catwalk.JIT])

A marker for instructing `CatwalkEx` to optimize transducers and reducing
function inner to this position. Optionally, it takes a callable to create a
Catwalk JIT object.

```
xs |> Map(f) |> OptimizeInner() |> Map(g)
                                   ------
                                  JIT target
```

It supports the type instability caused by the outer transducers (such as
`Map(f)`), collection (`xf`), and/or the reducing function output.
"""
struct OptimizeInner{Factory}
    factory::Factory
end

OptimizeInner(::Type{T}) where {T} = OptimizeInner{Type{T}}(T)
OptimizeInner() = OptimizeInner(Catwalk.JIT)

(f::OptimizeInner{Factory})(outer::Outer) where {Factory,Outer} =
    OptimizeInnerFoldable{Factory,Outer}(f.factory, outer)

struct OptimizeInnerFoldable{Factory,Outer}
    factory::Factory
    outer::Outer
end

function Transducers.transduce(xf, rf, init, coll, exc::CatwalkEx)
    xfi, opt = extract_transducer(xf(coll))
    opt::OptimizeInnerFoldable
    xfo, itr = extract_transducer(opt.outer)
    rf0 = reducingfunction(xfi ∘ xfo, rf)
    acc = start(rf0, init)
    y = iterate(itr)
    y === nothing && return complete(rf0, acc)
    acc = next(rf0, acc, first(y))
    state = last(y)
    batchsize = exc.batchsize
    @assert batchsize > 0
    jit = opt.factory()
    while true
        Catwalk.step!(jit)
        rf1 = reducingfunction(xfi ∘ OptimizeXF(Catwalk.ctx(jit)) ∘ xfo, rf)
        acc, state, done = onebatch(rf1, acc, itr, state, batchsize)
        done && return complete(rf1, acc)
    end
end

function onebatch(rf::RF, acc, itr, state, counter) where {RF}
    jitctx = xform(inner(rf)).jitctx
    while counter != 0
        y = iterate(itr, state)
        y === nothing && return acc, state, true
        state = last(y)
        val = next(rf, acc, first(y), jitctx)
        val isa Reduced && return val, state, true
        acc = val
        counter -= 1
    end
    return acc, state, false
end

struct OptimizeXF{C} <: Transducer
    jitctx::C
end

@inline @jit next mapped function Transducers.next(rf::R_{Map}, acc, input, jitctx)
    mapped = xform(rf).f(input)
    next(inner(rf), acc, mapped)
end

function Transducers.next(rf::R_{OptimizeXF}, acc, input)
    next(inner(rf), acc, input)
end

end # module FoldsCatwalk


text:
"""
	stationarydist(m::MarkovProcess)

Phân phối dừng của quá trình Markov
"""
function stationarydist(m::MarkovProcess)
	G = m.G
	A = vcat(transpose(G)[1:end-1, :], ones(1, size(G, 1)))
	B = vcat(zeros(size(G, 1) - 1), 1)
	inv(A) * B
end

"""
	imbeddedchain(m::MarkovProcess)

Xích Markov nhúng của một quá trình Markov
"""
function imbeddedchain(m::MarkovProcess)
	λ = -diag(m.G)
	p = m.G ./ λ - I
	p[diagind(p)] .= 0
	p[diagind(p)] = 1 .- sum(p; dims=2)
	MarkovChain(p)
end

"""
	transition_matrix(m::MarkovProcess, t)

Tạo ma trận chuyển tại thời điểm t
"""
function transition_matrix(m::MarkovProcess, t)
	D, Q = eigen(m.G)
	Q * diagm(exp.(D * t)) * inv(Q)
end

function Base.show(io::IO, m::MarkovProcess)
	display("Quá trình Markov với ma trận sinh")
	display(m.G)
end


text:
# using Revise
using TopOpt, LinearAlgebra, StatsFuns
using StatsFuns: logsumexp
# using Makie
# using TopOpt.TopOptProblems.Visualization: visualize

E = 1.0 # Young’s modulus
v = 0.3 # Poisson’s ratio
f = 1.0 # downward force
rmin = 3.0

problems = Any[
    PointLoadCantilever(Val{:Linear}, (160, 40), (1.0, 1.0), E, v, f),
    HalfMBB(Val{:Linear}, (60, 20), (1.0, 1.0), E, v, f),
    LBeam(Val{:Linear}, Float64),
]
problem_names = ["Cantilever beam", "Half MBB beam", "L-beam"]

for i in 1:length(problems)
    println(problem_names[i])
    problem = problems[i]
    # Parameter settings
    V = 0.5 # volume fraction
    xmin = 0.0001 # minimum density
    steps = 40 # maximum number of penalty steps, delta_p0 = 0.1
    convcriteria = KKTCriteria()
    solver = FEASolver(Direct, problem; xmin=xmin)
    x0 = fill(1.0, length(solver.vars))
    for p in [1.0, 2.0, 3.0]
        #penalty = TopOpt.PowerPenalty(1.0)
        global penalty = TopOpt.PowerPenalty(p)
        # Define a finite element solver
        solver = FEASolver(Direct, problem; xmin=xmin, penalty=penalty)
        # Define compliance objective
        global stress = TopOpt.MicroVonMisesStress(solver)
        global filter = DensityFilter(solver; rmin=rmin)
        global volfrac = Volume(problem, solver)

        obj = x -> volfrac(filter(x)) - V
        constr = x -> begin
            s = stress(filter(x))
            thr = 10
            vcat((s .- thr) / 100, logsumexp(s) - log(length(s)) - thr)
        end
        alg = PercivalAlg()
        options = PercivalOptions()
        optimizer = Optimizer(obj, constr, x0, alg; options=options)
        # Define continuation SIMP optimizer
        simp = SIMP(optimizer, solver, p)
        # Solve
        global result = simp(x0)
        x0 = result.topology
    end
    #visualize(
    #    problem; topology = penalty.(filter(result.topology)), default_exagg_scale = 0.07,
    #    scale_range = 10.0, vector_linewidth = 3, vector_arrowsize = 0.5,
    #)
end


text:
@info "starting 1D1V model"
TD = 30000 # eV
Te = 10000 # eV

D = species.D
e = species.e

function HotCarrier(T)
    Kb = 8.617333262145e-5
    function P(x,v)
        v_ = sqrt(sum(v .^2))
        exp(-v_/(Kb*T))
    end
end

D_D = Distribution(HotCarrier(TD), D)
G = Geometry()

plasma = ElectrostaticPlasma([D_D], G)

sol1 = Plasma.solve(plasma, dim=1, GPU=false, inner_layers=32, time_ub=100.0, ub=10.0; strategy=QuadratureTraining())
sol2 = Plasma.solve(plasma, dim=1, GPU=false, inner_layers=32, time_ub=100.0, ub=10.0; strategy=StochasticTraining(200))

@save "1d1v_quadrature.bson" sol1
@save "1d1v_quadrature.bson" sol2

Plasma.plot(sol1)
Plasma.plot(sol2)


text:
module SurrogateModelOptim

export  smoptimize,
        surrogate_model,
        model_infill,
        scaled_LHC_sampling_plan

using LatinHypercubeSampling
import LatinHypercubeSampling.Continuous, LatinHypercubeSampling.Categorical
using ScatteredInterpolation
using BlackBoxOptim
using Parameters
using Distances
using Statistics
using StatsBase
using LinearAlgebra
using Distributed
using StaticArrays
import NearestNeighbors: KDTree, knn
using Printf

include("types.jl")
include("LHC_sampling_plan.jl")
include("model_infill_utilities.jl")
include("model_infill.jl")
include("surrogate_model_utilities.jl")
include("surrogate_model.jl")
include("smoptimize.jl")

end # module


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

function SCIPincludeConshdlrKnapsack(scip)
    ccall((:SCIPincludeConshdlrKnapsack, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPcreateConsKnapsack(scip, cons, name, nvars, vars, weights, capacity, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsKnapsack, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Clonglong}, Clonglong, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, nvars, vars, weights, capacity, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicKnapsack(scip, cons, name, nvars, vars, weights, capacity)
    ccall((:SCIPcreateConsBasicKnapsack, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Clonglong}, Clonglong), scip, cons, name, nvars, vars, weights, capacity)
end

function SCIPaddCoefKnapsack(scip, cons, var, weight)
    ccall((:SCIPaddCoefKnapsack, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Clonglong), scip, cons, var, weight)
end

function SCIPgetCapacityKnapsack(scip, cons)
    ccall((:SCIPgetCapacityKnapsack, libscip), Clonglong, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPchgCapacityKnapsack(scip, cons, capacity)
    ccall((:SCIPchgCapacityKnapsack, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Clonglong), scip, cons, capacity)
end

function SCIPgetNVarsKnapsack(scip, cons)
    ccall((:SCIPgetNVarsKnapsack, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetVarsKnapsack(scip, cons)
    ccall((:SCIPgetVarsKnapsack, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetWeightsKnapsack(scip, cons)
    ccall((:SCIPgetWeightsKnapsack, libscip), Ptr{Clonglong}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetDualsolKnapsack(scip, cons)
    ccall((:SCIPgetDualsolKnapsack, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetDualfarkasKnapsack(scip, cons)
    ccall((:SCIPgetDualfarkasKnapsack, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetRowKnapsack(scip, cons)
    ccall((:SCIPgetRowKnapsack, libscip), Ptr{SCIP_ROW}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPsolveKnapsackExactly(scip, nitems, weights, profits, capacity, items, solitems, nonsolitems, nsolitems, nnonsolitems, solval, success)
    ccall((:SCIPsolveKnapsackExactly, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cint, Ptr{Clonglong}, Ptr{Cdouble}, Clonglong, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{UInt32}), scip, nitems, weights, profits, capacity, items, solitems, nonsolitems, nsolitems, nnonsolitems, solval, success)
end

function SCIPsolveKnapsackApproximately(scip, nitems, weights, profits, capacity, items, solitems, nonsolitems, nsolitems, nnonsolitems, solval)
    ccall((:SCIPsolveKnapsackApproximately, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cint, Ptr{Clonglong}, Ptr{Cdouble}, Clonglong, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}), scip, nitems, weights, profits, capacity, items, solitems, nonsolitems, nsolitems, nnonsolitems, solval)
end

function SCIPseparateKnapsackCuts(scip, cons, sepa, vars, nvars, weights, capacity, sol, usegubs, cutoff, ncuts)
    ccall((:SCIPseparateKnapsackCuts, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_SEPA}, Ptr{Ptr{SCIP_VAR}}, Cint, Ptr{Clonglong}, Clonglong, Ptr{SCIP_SOL}, UInt32, Ptr{UInt32}, Ptr{Cint}), scip, cons, sepa, vars, nvars, weights, capacity, sol, usegubs, cutoff, ncuts)
end

function SCIPseparateRelaxedKnapsack(scip, cons, sepa, nknapvars, knapvars, knapvals, valscale, rhs, sol, cutoff, ncuts)
    ccall((:SCIPseparateRelaxedKnapsack, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_SEPA}, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Cdouble, Cdouble, Ptr{SCIP_SOL}, Ptr{UInt32}, Ptr{Cint}), scip, cons, sepa, nknapvars, knapvars, knapvals, valscale, rhs, sol, cutoff, ncuts)
end


text:
include("ResultTests.jl")
include("TestFuncTests.jl")
include("RulesTests.jl")


text:
using DataDeps

register(DataDep(
    "iris",
    """
    Dataset: iris
    Author: R. A. Fisher (donated by Michael Marshall)
    Website: https://archive.ics.uci.edu/ml/datasets/iris

    Observations: 150
    Features:     4
    Classes:      3
    """,
    "https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data",
    "6f608b71a7317216319b4d27b4d9bc84e6abd734eda7872b71a458569e2656c0",
    post_fetch_method = (path -> begin
        MVDA.process_dataset(path,
            header=false,
            target_index=5,
            feature_indices=1:4,
            ext=".csv",
        )
    end),
))


text:
include("det_loocv.jl")
using Test

#####
##### test det_loocv
#####
#the inremental_cholesky is just a wrapper type that allows one to easily extend the cholesky
#decomposition without allocation more space at each iteration, because a fixed amount of space
#is allocated at the beginning, and what you get to use is a view of the top left n x n corner
#println("test det_loocv.jl...")
n=5
A = rand(n, n); A = A'*A + UniformScaling(4.0)
#IC = incremental_cholesky!(A, n-1) #make it n-1 just for kicks
IC = cholesky(A)

vals = zeros(1, n)
vals1 = zeros(1, n)
for i=1:n
    subdet = det(A[[1:i-1;i+1:end], [1:i-1;i+1:end]])
    vals[i] = subdet
    #println("actual subdets $i: ", subdet)
end

det_original = det(IC)
for i = 1:n
    subdet = det_loocv(IC, det_original, i)
    vals1[i] = subdet
    #println("fast subdets $i: ", subdet)
end

@testset "det_loocv.jl" begin
    @test norm(vals-vals1)<1e-10
end

#####
##### test det_XΣX_loocv
#####
#println("test det_XΣX_loocv...")

n = 20
m = 3
X = rand(n, m)
A = (w = rand(n, n);w'*w + UniformScaling(5))
cholA = cholesky(A)
detA = det(cholA)
bigmat = [A X; X' zeros(m, m)]
#for i = 1:n
#    println("actual subdet $i: ", det(bigmat[[1:i-1;i+1:end], [1:i-1;i+1:end]]))
#end
vals = zeros(1, n)
vals1 = zeros(1, n)
for i = 1:n
    X_minus_i = X[[1:i-1;i+1:end], :]
    sigma_minus_i = A[[1:i-1;i+1:end], [1:i-1;i+1:end]]
    vals[i] = det(X_minus_i'*(sigma_minus_i\X_minus_i))
    #println("actual subdet $i: ", det(X_minus_i'*(sigma_minus_i\X_minus_i)))
end
for i =1:n
    subdet =  det_XΣX_loocv(X, cholA, detA, i)
    vals1[i] = subdet
    #println("fast subdet $i: ", subdet)
end

@testset "det_XΣX_loocv" begin
    @test norm(vals-vals1)<1e-10
end

#####
##### test logdet_loocv_IC
#####
#println("test logdet_loocv_IC...")
n=5
A = rand(n, n); A = A'*A + UniformScaling(4.0)
B = copy(A)
IC = incremental_cholesky!(B, n)
#IC = cholesky(A)
vals = zeros(1, n)
for i=1:n
    sublogdet = log(det(A[[1:i-1;i+1:end], [1:i-1;i+1:end]]))
    vals[i] = sublogdet
    #println("actual sublogdets $i: ", sublogdet)
end
vals1 = zeros(1, n)
logdet_original = log(det(IC))
for i = 1:n
    sublogdet = logdet_loocv_IC(IC, logdet_original, i)
    vals1[i] = sublogdet
    #println("fast sublogdet $i: ", sublogdet)
end
@testset "logdet_loocv_IC" begin
    @test norm(vals-vals1)<1e-10
end

#####
##### test logdet_loocv
#####
#println("test logdet_loocv...")
n=5
#A = rand(n, n); A = A'*A + UniformScaling(1.0)
#println("A: ", A)
val = zeros(1, n)
chol = cholesky(A)
for i=1:n
    sublogdet = log(det(A[[1:i-1;i+1:end], [1:i-1;i+1:end]]))
    vals[i] = sublogdet
    #println("actual sublogdets $i: ", sublogdet)
end

val1 =zeros(1, n)
logdet_original = log(det(chol))
for i = 1:n
    sublogdet = logdet_loocv(chol, logdet_original, i)
    vals1[i] = sublogdet
    #println("fast sublogdet $i: ", sublogdet)
end

@testset "logdet_loocv" begin
    @test norm(vals-vals1)<1e-10
end

#####
##### test logdet_XΣX_loocv
#####
#println("Running test logdet_XΣX_loocv...")
n = 30
m = 10
X = rand(n, m)
A = (w = rand(n, n);w'*w + UniformScaling(5))
cholA = cholesky(A)
logdetA = logdet(cholA)
bigmat = [A X; X' zeros(m, m)]
#for i = 1:n
#    println("actual subdet $i: ", det(bigmat[[1:i-1;i+1:end], [1:i-1;i+1:end]]))
#end
vals = zeros(1, n)
@time begin
for i = 1:n
    X_minus_i = X[[1:i-1;i+1:end], :]
    sigma_minus_i = A[[1:i-1;i+1:end], [1:i-1;i+1:end]]
    vals[i] = log(abs(det(X_minus_i'*(sigma_minus_i\X_minus_i))))
    #println("actual sublogdet $i: ", log(abs(det(X_minus_i'*(sigma_minus_i\X_minus_i)))))
end
end
vals1 = zeros(1, n)
@time begin
for i =1:n
    sublogdet =  logdet_XΣX_loocv(X, cholA, logdetA, i)
    vals1[i] = sublogdet
    #println("fast sublogdet $i: ", sublogdet)
end
end
println("norm of errors: ", norm(vals-vals1))

@testset "logdet_XΣX_loocv" begin
    @test norm(vals-vals1)<1e-10
end

#####
##### test logdet_XΣX_loocv_IC
#####
#println("Running test logdet_XΣX_loocv_IC...")
n = 30
m = 10
X = rand(n, m)
#A = (w = rand(n, n);w'*w + UniformScaling(5))
B = copy(A) #from previous
cholA = incremental_cholesky!(B, size(A)[1])
logdetA = logdet(cholA)
bigmat = [A X; X' zeros(m, m)]
#for i = 1:n
#    println("actual subdet $i: ", det(bigmat[[1:i-1;i+1:end], [1:i-1;i+1:end]]))
#end
vals = zeros(1, n)
@time begin
for i = 1:n
    X_minus_i = X[[1:i-1;i+1:end], :]
    sigma_minus_i = A[[1:i-1;i+1:end], [1:i-1;i+1:end]]
    vals[i] = log(abs(det(X_minus_i'*(sigma_minus_i\X_minus_i))))
    #println("actual sublogdet $i: ", log(abs(det(X_minus_i'*(sigma_minus_i\X_minus_i)))))
end
end
vals1 = zeros(1, n)
@time begin
for i =1:n
    sublogdet =  logdet_XΣX_loocv_IC(X, cholA, logdetA, i)
    vals1[i] = sublogdet
    #println("fast sublogdet $i: ", sublogdet)
end
end
#println("norm of errors: ", norm(vals-vals1))

@testset "logdet_XΣX_loocv_IC" begin
    @test norm(vals-vals1)<1e-10
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module MEngProject
include("./LaminartFunc.jl")
include("./LaminartInitFunc.jl")
include("./LaminartEqImfilter.jl")
include("./LaminartEqImfilterGPU_FFT.jl")
include("./LaminartEqImfilterGPU_FIR.jl")
include("./LaminartEqImfilterGPU_IIR.jl")
include("./LaminartEqConv.jl")
include("./LaminartKernels.jl")
include("./Utils.jl")
include("./Parameters.jl")
export LaminartFunc,
    LaminartInitFunc,
    LaminartEqImfilter,
    LaminartEqConv,
    LaminartKernels,
    Utils,
    Parameters,
    LaminartEqImfilterGPU_FFT,
    LaminartEqImfilterGPU_FIR,
    LaminartEqImfilterGPU_IIR
end


text:
module sacco_comp_pl_examples
# A mixed-enhanced finite-element for the analysis of laminated composite plates
# FERDINANDO AURICCHIO, ELIO SACCO

# Laminated [0/90] cross-ply square simply supported plate. Uniform Distributed
# Loading. Results in Table IV.

# Note: The scaling of the horizontal displacement of the midpoint of the edge
# is wrong in the paper. The power of the aspect ratio should also be 4.
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test(visualize = true)
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule

    a = 250 * phun("mm") # note the changed dimension
    nx = ny = 28
    # high modulus orthotropic graphite=epoxy composite material
    E1 = 100 * phun("GPa")
    E2 = E1 / 25
    G12 = G13 = E2 * 0.5
    nu12 = 0.25
    nu23 = 0.25
    G23 = E2 * 0.2
    q = 10.0 * phun("kilo*Pa")
    thickness = a / 10
    # With these inputs, the Abaqus-verified solution is -1176.8346659735575

    tolerance = a / nx / 100
    CM = CompositeLayupModule

    mater = CM.lamina_material(E1, E2, nu12, G12, G13, G23)
    plies = CM.Ply[
        CM.Ply("ply_0", mater, thickness / 2, 0),
        CM.Ply("ply_90", mater, thickness / 2, 90)
    ]
    mcsys = CM.cartesian_csys((1, 2, 3))
    layup = CM.CompositeLayup("sacco", plies, mcsys)

    fens, fes = T3block(a, a, nx, ny)
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz, 1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz, 1), 6))

    # Apply EBC's
    # Simple support
    # The boundary conditions are a bit peculiar: refer to the paper
    l1 = selectnode(fens; box = Float64[0 0 -Inf Inf 0 0], inflate = tolerance)
    for i in [2, 4]
        setebc!(dchi, l1, true, i)
    end
    l1 = selectnode(fens; box = Float64[a a -Inf Inf 0 0], inflate = tolerance)
    for i in [2, 4]
        setebc!(dchi, l1, true, i)
    end
    l1 = selectnode(fens; box = Float64[-Inf Inf 0 0 0 0], inflate = tolerance)
    for i in [1, 5]
        setebc!(dchi, l1, true, i)
    end
    l1 = selectnode(fens; box = Float64[-Inf Inf a a 0 0], inflate = tolerance)
    for i in [1, 5]
        setebc!(dchi, l1, true, i)
    end
    # Simple support
    l1 = connectednodes(meshboundary(fes))
    for i in [3]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi)

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi)

    lfemm = FEMMBase(IntegDomain(fes, TriRule(1)))
    fi = ForceIntensity(FFlt[0, 0, q, 0, 0, 0])
    F = distribloads(lfemm, geom0, dchi, fi, 2)

    # Solve
    U = K \ F
    scattersysvec!(dchi, U[:])

    # Visualization
    if visualize
        scalars = []
        for nc in 1:3
            fld = fieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = layup.csys)
            push!(scalars, ("m$nc", fld.values))
            fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = layup.csys)
            push!(scalars, ("em$nc", fld.values))
        end
        vtkwrite("sacco_comp_pl_examples-m.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
        scalars = []
        for nc in 1:3
            fld = fieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = layup.csys)
            push!(scalars, ("n$nc", fld.values))
            fld = elemfieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = layup.csys)
            push!(scalars, ("en$nc", fld.values))
        end
        vtkwrite("sacco_comp_pl_examples-n.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
        scalars = []
        for nc in 1:2
            fld = fieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = layup.csys)
            push!(scalars, ("q$nc", fld.values))
            fld = elemfieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = layup.csys)
            push!(scalars, ("eq$nc", fld.values))
        end
        vtkwrite("sacco_comp_pl_examples-q.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])

    end
    true
end

function allrun()
    println("#####################################################")
    println("# test_convergence ")
    test()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
##

using ACE
using Printf, Test, LinearAlgebra, StaticArrays
using ACE: evaluate, evaluate_d, Rn1pBasis, Ylm1pBasis,
      PositionState, Product1pBasis
using Random: shuffle
using ACEbase.Testing: dirfdtest, fdtest, print_tf, test_fio

##

@info "Build a 1p basis from scratch"

maxdeg = 5
r0 = 1.0
rcut = 3.0
maxorder = 3
Bsel = SimpleSparseBasis(maxorder, maxdeg)

trans = PolyTransform(1, r0)   # r -> x = 1/r^2
J = transformed_jacobi(maxdeg, trans, rcut; pcut = 2)   #  J_n(x) * (x - xcut)^pcut
Rn = Rn1pBasis(J)
Ylm = Ylm1pBasis(maxdeg)
B1p = Product1pBasis( (Rn, Ylm) )
ACE.init1pspec!(B1p, Bsel)

nX = 10
Xs = rand(PositionState{Float64}, Rn, nX)
cfg = ACEConfig(Xs)

A = evaluate(B1p, cfg)

@info("test against manual summation")
A1 = sum( evaluate(B1p, X) for X in Xs )
println(@test A1 ≈ A)

@info("test permutation invariance")
for ntest = 1:30
   print_tf(@test A ≈ evaluate(B1p, ACEConfig(shuffle(Xs))))
end
println()

##

# test_fio(Ylm)
# D = ACE.write_dict(Ylm)
# Ylm_ = ACE.read_dict(D)
# Ylm.SH.alp == Ylm_.SH.alp
# Ylm.SH.alp.B ≈ Ylm_.SH.alp.B
# @which ( Ylm.SH.alp == Ylm_.SH.alp )

@info("Test FIO")
for _B in (J, Rn, Ylm, B1p)
   println((@test(all(test_fio(_B)))))
end

##

@info("Ylm1pBasis gradients")
Y = ACE.acquire_B!(Ylm, Xs[1])
dY = ACE.acquire_dB!(Ylm, Xs[1])
println(@test (typeof(dY) == eltype(Ylm.dB_pool.arrays[Base.Threads.threadid()])))
ACE.evaluate!(Y, Ylm, Xs[1])
ACE.evaluate_d!(dY, Ylm, Xs[1])
Y1 = ACE.acquire_B!(Ylm, Xs[1])
dY1 = ACE.acquire_dB!(Ylm, Xs[1])
ACE.evaluate_ed!(Y1, dY1, Ylm, Xs[1])

println(@test (evaluate(Ylm, Xs[1]) ≈ Y))
println(@test (evaluate_d(Ylm, Xs[1]) ≈ dY))
println(@test ((Y ≈ Y1) && (dY ≈ dY1)) )

_vec2X(x) = PositionState{Float64}((rr = SVector{3}(x),))

for ntest = 1:30
   x0 = randn(3)
   c = rand(length(Y))
   F = x -> sum(ACE.evaluate(Ylm, _vec2X(x)) .* c)
   dF = x -> sum(ACE.evaluate_d(Ylm, _vec2X(x)) .* c).rr |> Vector
   print_tf(@test fdtest(F, dF, x0; verbose=false))
end
println()
##

@info("Rn1pBasis gradients")

for ntest = 1:30
   x0 = randn(3)
   c = rand(length(Rn))
   F = x -> sum(ACE.evaluate(Rn, _vec2X(x)) .* c)
   dF = x -> sum(ACE.evaluate_d(Rn, _vec2X(x)) .* c).rr |> Vector
   print_tf(@test fdtest(F, dF, x0; verbose=false))
end
println()

##

@info("Product basis evaluate_ed! tests")

A1 = ACE.acquire_B!(B1p, cfg)
ACE.evaluate!(A1, B1p, cfg)
A2 = ACE.acquire_B!(B1p, cfg)
dA = ACE.acquire_dB!(B1p, cfg)
ACE.evaluate_ed!(A2, dA, B1p, cfg)
println(@test A1 ≈ A2)

println(@test( evaluate_d(B1p, cfg) ≈ dA ))

##
@info("Product basis gradient test")

for ntest = 1:30
   x0 = randn(3)
   c = rand(length(B1p))
   F = x -> sum(ACE.evaluate(B1p, _vec2X(x)) .* c)
   dF = x -> sum(ACE.evaluate_d(B1p, ACEConfig([_vec2X(x)])) .* c).rr |> Vector
   print_tf(@test fdtest(F, dF, x0; verbose=false))
end
println()
##


text:
include("learn_perceptron.jl")
using DataArrays

# create initial weight vector
w_init = randn(3,1);

# Create initial train data
train_data = @data([0.80857 0.83721;
                    0.35714 0.8505;
                    -0.75143 -0.7309;
                    -0.3 0.12625;
                    0.64286 -0.54485]);

target = @data([-1.0; -1.0; -1.0; -1.0;-1.0]);

train_data_biased = add_bias(train_data);

number_of_errors = eval_perceptron(train_data_biased, w_init, target)
println("The errors before starting training: ", number_of_errors)
while (number_of_errors > 0)
  w_init = learn_perceptron(train_data_biased, target, w_init)
  number_of_errors = eval_perceptron(train_data_biased, w_init, target)
end
println("Desired weight vector for a given train set: ")
println(w_init)

max_stars_repo_name: kirnap/Machine-Learning-in-Julia

max_issues_repo_name: kirnap/Machine-Learning-in-Julia

max_forks_repo_name: kirnap/Machine-Learning-in-Julia


text:
# Autogenerated wrapper script for helpme_jll for x86_64-apple-darwin
export libhelpme

using FFTW_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("helpme")
JLLWrappers.@declare_library_product(libhelpme, "@rpath/libhelpme.dylib")
function __init__()
    JLLWrappers.@generate_init_header(FFTW_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libhelpme,
        "lib/libhelpme.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# ╔═╡ 13b71454-28e9-11eb-3dfd-d591460d236a
using Gnuplot

# ╔═╡ c898b298-5d3a-4734-adf2-b69a76f4aede
md"""
# Drude Conductivity
#### Jarv

Finally got around (Nov 2020) to actually plotting the Ashcroft & Mermin freq.dependent Drude response.

"""

# ╔═╡ bc88379c-28e8-11eb-1fa6-87caf0f24fb9
σ(ω; τ=0.5, σ0=1.0) = σ0/(1-im*ω*τ)
# in the Drude theory σ0 is ne^2τ/m - see (1.29) p. 16 Ashcroft & Mermin

# ╔═╡ 2a16bede-28e9-11eb-099f-119bdeceaca1
ω=0:0.1:20

# ╔═╡ 22ef9f21-a008-41ec-827d-b29636ae27ec
md"""
To get the figure as you see in the textbook, you need to plot on a linear-log scale!
"""

# ╔═╡ 90ffdc66-28e9-11eb-1191-89949f00b0b7
begin
	@gp "set xlabel 'Frequency. ω'"
	@gp :- "set ylabel 'Conductivity, σ'"
	@gp :- ω real(σ.(ω)) "w l title 'Drude Real'"
	@gp :- ω imag(σ.(ω)) "w l title 'Drude Imaginary'"
	# pay attention to the logscale!
	@gp :- "set logscale x"
end

# ╔═╡ f039f47c-28ea-11eb-251e-9f923ac538a0
# Without logscle, for comparison
@gp :- "unset logscale x"

# ╔═╡ da70d77e-28ee-11eb-185e-47d5528a1743
@gp :- "set logscale y"

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Gnuplot = "dc211083-a33a-5b79-959f-2ff34033469d"

[compat]
Gnuplot = "~1.4.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "a325370b9dd0e6bf5656a6f1a7ae80755f8ccc46"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.7.2"

[[ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "a851fec56cb73cfdf43762999ec72eff5b86882a"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.15.0"

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "31d0151f5716b655421d9d75b7fa74cc4e744df2"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.39.0"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[Gnuplot]]
deps = ["ColorSchemes", "ColorTypes", "Colors", "DataStructures", "REPL", "ReplMaker", "StatsBase", "StructC14N", "Test"]
git-tree-sha1 = "637ee67f78c4adb222ee11c590bd1d0a1a50ba88"
uuid = "dc211083-a33a-5b79-959f-2ff34033469d"
version = "1.4.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[LogExpFunctions]]
deps = ["ChainRulesCore", "DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "34dc30f868e368f8a17b728a1238f3fcda43931a"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.3"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[ReplMaker]]
deps = ["REPL", "Unicode"]
git-tree-sha1 = "76098218397ec93b925b70ce355144d539b1a8b4"
uuid = "b873ce64-0db9-51f5-a568-4457d8e49576"
version = "0.2.5"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "65fb73045d0e9aaa39ea9a29a5e7506d9ef6511f"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.11"

[[StructC14N]]
deps = ["Test"]
git-tree-sha1 = "932b930747486efdc5f78297480c311d58374d04"
uuid = "d2514e9c-36c4-5b8e-97e2-51e7675c221c"
version = "0.3.0"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─c898b298-5d3a-4734-adf2-b69a76f4aede
# ╠═bc88379c-28e8-11eb-1fa6-87caf0f24fb9
# ╠═13b71454-28e9-11eb-3dfd-d591460d236a
# ╠═2a16bede-28e9-11eb-099f-119bdeceaca1
# ╟─22ef9f21-a008-41ec-827d-b29636ae27ec
# ╠═90ffdc66-28e9-11eb-1191-89949f00b0b7
# ╠═f039f47c-28ea-11eb-251e-9f923ac538a0
# ╠═da70d77e-28ee-11eb-185e-47d5528a1743
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: Frost-group/Pluto-notebooks

max_issues_repo_name: Frost-group/Pluto-notebooks

max_forks_repo_name: Frost-group/Pluto-notebooks


text:
println("Advent of Code 2017, day 5")

f = open("input5")
input = readlines(f)
close(f)

jumps = [parse(Int,x) for x in input]

# Part 1

l = length(jumps)
steps = 0
i = 1
while 1 <= i <= l
  j = jumps[i]
  jumps[i] = j+1
  global i += j
  global steps += 1
end

println("Part 1: ", steps)

# Part 2

jumps = [parse(Int,x) for x in input]

steps = 0
i = 1
while 1 <= i <= l
  j = jumps[i]
  if j >= 3
    jumps[i] = j-1
  else
    jumps[i] = j+1
  end
  global i += j
  global steps += 1
end

println("Part 2: ", steps)

max_stars_repo_name: venanzio/AoC

max_issues_repo_name: venanzio/AoC

max_forks_repo_name: venanzio/AoC


text:
# ------------------------------------------------------------------------------
# Schwarz preconditioner
# ------------------------------------------------------------------------------

"""
```julia
Schwarz(operator)
```

Schwarz diagonal preconditioner for finite element operators

# Arguments:
- `operator`: finite element operator to precondition

# Returns:
- Schwarz preconditioner object

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
mass = GalleryOperator("mass", 4, 4, mesh);

# preconditioner
schwarz = Schwarz(mass);

# verify
println(schwarz)
println(schwarz.operator)

# output
schwarz preconditioner
finite element operator:
2d mesh:
    dx: 1.0
    dy: 1.0

2 inputs:
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    interpolation
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    quadratureweights

1 output:
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    interpolation
```
"""
mutable struct Schwarz <: AbstractPreconditioner
    # data never changed
    operator::Operator

    # data empty until assembled
    operatorinverse::AbstractArray{Float64,2}

    # inner constructor
    Schwarz(operator::Operator) = new(operator)
end

# printing
# COV_EXCL_START
Base.show(io::IO, preconditioner::Schwarz) = print(io, "schwarz preconditioner")
# COV_EXCL_STOP

# ------------------------------------------------------------------------------
# data for computing symbols
# ------------------------------------------------------------------------------

"""
```julia
getoperatorinverse(preconditioner)
```

Compute or retrieve the inverse of the symbol matrix diagonal for a Schwarz
    preconditioner

# Returns:
- Symbol matrix diagonal inverse for the operator

# Example:
```jldoctest
# setup
mesh = Mesh1D(1.0);
diffusion = GalleryOperator("diffusion", 3, 3, mesh);

# preconditioner
schwarz = Schwarz(diffusion)

# note: either syntax works
diagonalinverse = LFAToolkit.getoperatorinverse(schwarz);
diagonalinverse = schwarz.operatorinverse;

# verify
@assert diagonalinverse ≈ [6/7 0; 0 3/4]

# output

```
"""
function getoperatorinverse(preconditioner::Schwarz)
    # assemble if needed
    if !isdefined(preconditioner, :operatorinverse)
        elementmatrix = preconditioner.operator.elementmatrix
        preconditioner.operatorinverse = pinv(Matrix(elementmatrix))
    end

    # return
    return getfield(preconditioner, :operatorinverse)
end

# ------------------------------------------------------------------------------
# get/set property
# ------------------------------------------------------------------------------

function Base.getproperty(preconditioner::Schwarz, f::Symbol)
    if f == :operatorinverse
        return getoperatorinverse(preconditioner)
    else
        return getfield(preconditioner, f)
    end
end

function Base.setproperty!(preconditioner::Schwarz, f::Symbol, value)
    if f == :operator
        throw(ReadOnlyMemoryError()) # COV_EXCL_LINE
    else
        return setfield!(preconditioner, f, value)
    end
end

# ------------------------------------------------------------------------------
# compute symbols
# ------------------------------------------------------------------------------

"""
```julia
computesymbols(preconditioner, ω, θ)
```

Compute or retrieve the symbol matrix for a Schwarz preconditioned operator

# Arguments:
- `preconditioner`: Schwarz preconditioner to compute symbol matrix for
- `ω`:              Nothing
- `θ`:              Fourier mode frequency array (one frequency per dimension)

# Returns:
- Symbol matrix for the Schwarz preconditioned operator

# Example:
```jldoctest
using LinearAlgebra

for dimension in 1:3
    # setup
    mesh = []
    if dimension == 1
        mesh = Mesh1D(1.0);
    elseif dimension == 2
        mesh = Mesh2D(1.0, 1.0);
    elseif dimension == 3
        mesh = Mesh3D(1.0, 1.0, 1.0);
    end
    diffusion = GalleryOperator("diffusion", 3, 3, mesh);

    # preconditioner
    schwarz = Schwarz(diffusion);

    # compute symbols
    A = computesymbols(schwarz, [], π*ones(dimension));

    # verify
    using LinearAlgebra;
    eigenvalues = real(eigvals(A));
    if dimension == 1
        @assert max(eigenvalues...) ≈ 1/7
    elseif dimension == 2
        @assert min(eigenvalues...) ≈ -1/14
    elseif dimension == 3
        @assert min(eigenvalues...) ≈ -0.33928571428571486
    end
end

# output

```
"""
function computesymbols(preconditioner::Schwarz, ω::Array, θ::Array)

    rowmodemap = preconditioner.operator.rowmodemap
    columnmodemap = preconditioner.operator.columnmodemap
    dimension = preconditioner.operator.dimension
    elementmatrix = preconditioner.operator.elementmatrix
    Minv = preconditioner.operatorinverse
    nodecoordinatedifferences = preconditioner.operator.nodecoordinatedifferences
    numberrows, numbercolumns = size(elementmatrix)

    symbolmatrixnodes = zeros(ComplexF64, numberrows, numbercolumns)
    for i = 1:numberrows, j = 1:numbercolumns
        symbolmatrixnodes[i, j] =
            Minv[i, j] *
            ℯ^(im * sum([θ[k] * nodecoordinatedifferences[i, j, k] for k = 1:dimension]))
    end
    symbolmatrixmodes = rowmodemap * symbolmatrixnodes * columnmodemap

    A = computesymbols(preconditioner.operator, θ)
    return I - symbolmatrixmodes * A
end

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


text:
using Rewrite
using Test

@theory! begin
    a => FreeTheory()
    b => FreeTheory()
    c => FreeTheory()
    f => FreeTheory()
    g => FreeTheory()
    h => FreeTheory()
    p => CTheory()
    q => CTheory()
    s => ACTheory()
    t => ACTheory()
end

a = @term(a)
b = @term(b)
c = @term(c)
f(xs...) = @term(f($(xs...)))
g(xs...) = @term(g($(xs...)))
h(xs...) = @term(h($(xs...)))
p(x, y) = @term(p($x, $y))
q(x, y) = @term(q($x, $y))
s(xs...) = @term(s($(xs...)))
t(xs...) = @term(t($(xs...)))
x, y, z = Variable(), Variable(), Variable()

@testset "construction" begin
    include("construction.jl")
end

@testset "match" begin
    include("match.jl")
end

@testset "rewrite" begin
    include("rewrite.jl")
end