text:
Base.@kwdef struct SigMFGlobal
    datatype::String = "i"
    sample_rate::Union{Float64, Nothing} = nothing
    version::String = _sigmf_default_version
    num_channels::Union{UInt, Nothing} = nothing
    sha512::Union{String, Nothing} = nothing
    offset::Union{UInt, Nothing} = nothing
    description::Union{String, Nothing} = nothing
    author::Union{String, Nothing} = nothing
    meta_doi::Union{String, Nothing} = nothing
    data_doi::Union{String, Nothing} = nothing
    recorder::Union{String, Nothing} = nothing
    license::Union{String, Nothing} = nothing
    hw::Union{String, Nothing} = nothing
    dataset::Union{String, Nothing} = nothing
    trailing_bytes::Union{String, Nothing} = nothing
    metadata_only::Union{Bool, Nothing} = nothing
    geolocation::Union{String, Nothing} = nothing # TODO: GEOJSON
    extensions::Union{Dict{String,String}, Nothing} = nothing
    collection::Union{String, Nothing} = nothing
end

Base.@kwdef struct SigMFCaptures
    sample_start::UInt = 0
    global_index::Union{UInt, Nothing} = nothing
    header_bytes::Union{UInt, Nothing} = nothing
    frequency::Union{Float64, Nothing} = nothing
    datetime::Union{String, Nothing} = nothing
end

Base.@kwdef struct SigMFAnnotations
    sample_start::UInt = 0
    sample_count::Union{UInt, Nothing} = nothing
    generator::Union{String, Nothing} = nothing
    label::Union{String, Nothing} = nothing
    comment::Union{String, Nothing} = nothing
    freq_lower_edge::Union{Float64, Nothing} = nothing
    freq_upper_edge::Union{Float64, Nothing} = nothing
    signal_detail_extension::Union{SignalDetailExtension, Nothing} = nothing
    signal_emitter_extension::Union{SignalEmitterExtension, Nothing} = nothing
end

Base.@kwdef struct SigMFMeta
    meta_global::SigMFGlobal = SigMFGlobal()
    meta_captures::Vector{SigMFCaptures} = SigMFCaptures[]
    meta_annotations::Vector{SigMFAnnotations} = SigMFAnnotations[]
end

Base.@kwdef struct SigMFCollection
    version::String = _sigmf_default_version
    description::Union{String, Nothing} = nothing
    author::Union{String, Nothing} = nothing
    collection_doi::Union{String, Nothing} = nothing
    license::Union{String, Nothing} = nothing
    extensions::Union{Dict{String,String}, Nothing} = nothing
    streams::Vector{Dict{String,String}} = Dict{String,String}[]
end

struct SigMFFile{T}
    meta::SigMFMeta
    data::T
end

## Struct type mapping for JSON 3

StructTypes.StructType(::Type{SigMFGlobal}) = StructTypes.Struct()

function StructTypes.names(::Type{SigMFGlobal})
    ((:datatype,        Symbol("core:datatype")),
     (:sample_rate,     Symbol("core:sample_rate")),
     (:version,         Symbol("core:version")),
     (:num_channels,    Symbol("core:num_channels")),
     (:sha512,          Symbol("core:sha512")),
     (:offset,          Symbol("core:offset")),
     (:description,     Symbol("core:description")),
     (:author,          Symbol("core:author")),
     (:meta_doi,        Symbol("core:meta_doi")),
     (:data_doi,        Symbol("core:data_doi")),
     (:recorder,        Symbol("core:recorder")),
     (:license,         Symbol("core:license")),
     (:hw,              Symbol("core:hw")),
     (:dataset,         Symbol("core:dataset")),
     (:trailing_bytes,  Symbol("core:trailing_bytes")),
     (:metadata_only,   Symbol("core:metadata_only")),
     (:geolocation,     Symbol("core:geolocation")),
     (:extensions,      Symbol("core:extensions")),
     (:collection,      Symbol("core:collection")))
end

StructTypes.StructType(::Type{SigMFCaptures}) = StructTypes.Struct()

function StructTypes.names(::Type{SigMFCaptures})
    ((:sample_start, Symbol("core:sample_start")),
    (:global_index, Symbol("core:global_index")),
    (:frequency, Symbol("core:frequency")),
    (:datetime, Symbol("core:datetime")))
end

StructTypes.StructType(::Type{SigMFAnnotations}) = StructTypes.Struct()

function StructTypes.names(::Type{SigMFAnnotations})
    ((:sample_start, Symbol("core:sample_start")),
    (:sample_count, Symbol("core:sample_count")),
    (:generator, Symbol("core:generator")),
    (:label, Symbol("core:label")),
    (:comment, Symbol("core:comment")),
    (:freq_lower_edge, Symbol("core:freq_lower_edge")),
    (:freq_upper_edge, Symbol("core:freq_upper_edge")),
    (:signal_detail_extension, Symbol("signal:detail")),
    (:signal_emitter_extension, Symbol("signal:emitter")))
end

StructTypes.StructType(::Type{SigMFMeta}) = StructTypes.Struct()

function StructTypes.names(::Type{SigMFMeta})
    ((:meta_global,  Symbol("global")),
    (:meta_captures,Symbol("captures")),
    (:meta_annotations,Symbol("annotations")))
end


text:
struct Configuration
    p::Ptr{configuration_t}
    key::id_t
    type::configuration_type_bitset_t
end

function Configuration(p::Ptr{configuration_t}, key::id_t)
    return Configuration(p, key, configuration_type(p, key))
end

function Configuration(ctl::Control)
    p = control_configuration(ctl)
    return Configuration(p, configuration_root(p))
end

Base.unsafe_convert(::Type{Ptr{configuration_t}}, c::Configuration) = c.p

function subkey(c::Configuration, name::String)
    c.type == configuration_type_map || error("Configuration $(c.key) is not a map type")
    if configuration_map_has_subkey(c, c.key, name)
        return configuration_map_at(c, c.key, name)
    else
        return nothing
    end
end

type(c::Configuration, key::id_t) = configuration_type(c,key)
type(c) = configuration_type(c,c.key)

function Base.getproperty(c::Configuration, v::Symbol)
    if v ∈ fieldnames(Configuration)
        return getfield(c, v)
    end
    s = String(v)
    k = subkey(c,s)
    k === nothing && error("type Configuration has no field $(String(v))")
    if type(c,k) == configuration_type_value
        if configuration_value_is_assigned(c, k)
            return configuration_value_get(c, k)
        else
            return nothing
        end
    end
    return Configuration(c.p, k)
end

function Base.propertynames(c::Configuration)
    return Symbol.(keys(c))
end
function description(c::Configuration, key::id_t)
    return configuration_description(c, key)
end
description(c::Configuration) = description(c, c.key)

function keys(c::Configuration)
    if c.type == configuration_type_map
        v = Vector{String}()
        for k ∈ 0:(configuration_map_size(c, c.key) - 1)
            push!(v, configuration_map_subkey_name(c,c.key, k))
        end
        return v
    else
        return nothing
    end
end

ext: jl


text:
using Commas
using HMMs
using PyPlot

spx = Commas.readcomma( joinpath( homedir(), "myenvironment", "dev", "SPXHMM", "data", "SPX" ) )
y = log.(1 .+ spx[:vwretd])

parameters = zeros(3)
HMMs.fit_mle!( HMMs.GenTDist, parameters, y, ones(length(y)), Dict{Symbol,Any}() )

delta = 0.0002
bins = collect(-0.01:delta:0.01)

clf()
hist( y, bins=bins )

midpoints = (bins[1:end-1] + bins[2:end])/2
densities = delta*length(y).*pdf.( HMMs.GenTDist( parameters... ), midpoints )
plot( midpoints, densities )

delta = 0.001
bins = collect(-0.25:delta:0.25)

clf()
hist( y, bins=bins, log=true )

midpoints = (bins[1:end-1] + bins[2:end])/2
densities = delta*length(y).*pdf.( HMMs.GenTDist( parameters... ), midpoints )
plot( midpoints, densities )

mu = mean(y)
sigma = std(y)
abs(mean(((y .- mu)./sigma).^3))^(1/3)

symy = [y; 2*mu .- y]
symparameters = zeros(3)
HMMs.fit_mle!( HMMs.GenTDist, symparameters, symy, ones(length(symy)), Dict{Symbol,Any}() )
parameters

t = HMMs.GenTDist( hmm.stateparameters[:,2]... )
using PyPlot
delta = 0.0005
bins = collect(-0.1:delta:0.1)
midpoints = (bins[1:end-1] + bins[2:end])/2
tdensities = pdf.( t, midpoints )
clf()
semilogy( midpoints, tdensities )
mu = hmm.stateparameters[1,2]
sigma = hmm.stateparameters[2,2]
nu = hmm.stateparameters[3,2]
densities = pdf.( Normal( mu, sigma*sqrt(nu/(nu-2)) ), midpoints )
plot( midpoints, densities )


text:
using ArcadeLearningEnvironment

episodes = 50

game = Game("space_invaders", 101)
println("repeat_action_probability: ", get_float(game, "repeat_action_probability"))
println("frame_skip: ", get_int(game, "frame_skip"))
println("color_averaging: ", get_bool(game, "color_averaging"))
println("record_screen_dir: ", get_string(game, "record_screen_dir"))
println("Actions: ", game.actions)
println("Inputs: ", typeof(get_inputs(game)))
println("RGB: ", typeof(get_rgb(game)))
draw(game)
close!(game)


text:
using FactCheck

include("mock/mock.jl")
include("mock/patch.jl")
include("registered_fixtures.jl")
include("fact_check_support.jl")
include("file_fixtures.jl")
include("fixture.jl")
include("matchers.jl")

# Throws errors when a @fact failed a test.
exitstatus()


text:
function  compute()
	for i=length(triangle)-1:-1:1
		for j=1: length(triangle[i])
			triangle[i][j]+=max(triangle[i+1][j],triangle[i+1][j+1])
		end
	end

	return triangle[1][1]

end

triangle = [  # Mutable
	[75],
	[95,64],
	[17,47,82],
	[18,35,87,10],
	[20, 4,82,47,65],
	[19, 1,23,75, 3,34],
	[88, 2,77,73, 7,63,67],
	[99,65, 4,28, 6,16,70,92],
	[41,41,26,56,83,40,80,70,33],
	[41,48,72,33,47,32,37,16,94,29],
	[53,71,44,65,25,43,91,52,97,51,14],
	[70,11,33,28,77,73,17,78,39,68,17,57],
	[91,71,52,38,17,14,91,43,58,50,27,29,48],
	[63,66, 4,68,89,53,67,30,73,16,69,87,40,31],
	[ 4,62,98,27,23, 9,70,98,73,93,38,53,60, 4,23],
]

println(compute())


text:
using GaussianBasis
using Molecules
using TensorOperations

function gradient(wfn::RHF, eri_type::Fermi.Integrals.Chonky)

    # Following eq. on C.3. Szabo & Ostlund
    atoms = wfn.molecule.atoms
    Natoms = length(atoms)
    bset = BasisSet(wfn.orbitals.basis, atoms)
    nbas = bset.nbas

    ∂E = zeros(Natoms, 3)

    @views Co = wfn.orbitals.C[:,1:wfn.ndocc]

    P = 2.0 * Co * Co'
    Q = 2.0 * Co * diagm(wfn.orbitals.eps[1:wfn.ndocc]) * Co'

    # Preallocate arrays
    ∂H = zeros(nbas, nbas, 3)
    ∂S = zeros(nbas, nbas, 3)
    ∂ERI = zeros(nbas, nbas, nbas, nbas, 3)

    # Intermediate auxiliary arrays
    AUX_ERI = zeros(nbas, nbas, nbas, nbas)
    AUX = zeros(nbas, nbas)

    for a in eachindex(atoms)

        ∂H .= 0
        ∂S .= 0
        ∂ERI .= 0

        # Nuclear repulsion
        ∂E[a, :] .= Molecules.∇nuclear_repulsion(atoms, a)

        # Store kinetic into H and nuclear into S
        GaussianBasis.∇kinetic!(∂H, bset, a)
        GaussianBasis.∇nuclear!(∂S, bset, a)

        ∂H .+= ∂S
        ∂S .= 0

        # Now use S for overlap
        GaussianBasis.∇overlap!(∂S, bset, a)

        GaussianBasis.∇ERI_2e4c!(∂ERI, bset, a)

        for q in 1:3
            @views vH = ∂H[:,:,q]
            AUX .= P .* vH
            ∂E[a, q]  = sum(AUX)

            @views vS = ∂S[:,:,q]
            AUX .= Q .* vS
            ∂E[a, q] -= sum(AUX)

            @views ERIq = ∂ERI[:,:,:,:,q]
            permutedims!(AUX_ERI, ERIq, (1,4,3,2))
            AUX_ERI .*= -0.5
            AUX_ERI .+= ERIq

            @tensoropt X = 0.5*P[μ,ν]*P[λ,σ]*AUX_ERI[μ,ν,σ,λ]
            ∂E[a,q] += X
        end
    end

    return ∂E
end

function gradient(wfn::RHF)

    # Following eq. on C.3. Szabo & Ostlund
    atoms = wfn.molecule.atoms
    Natoms = length(atoms)
    bset = BasisSet(wfn.orbitals.basis, atoms)
    nbas = bset.nbas

    ∂E = zeros(Natoms, 3)

    @views Co = wfn.orbitals.C[:,1:wfn.ndocc]

    P = 2.0 * Co * Co'
    Q = 2.0 * Co * diagm(wfn.orbitals.eps[1:wfn.ndocc]) * Co'

    # Preallocate arrays
    ∂H = zeros(nbas, nbas, 3)
    ∂S = zeros(nbas, nbas, 3)

    # Intermediate auxiliary arrays
    AUX = zeros(nbas, nbas)

    for a in eachindex(atoms)

        ∂H .= 0
        ∂S .= 0

        # Nuclear repulsion
        ∂E[a, :] .= Molecules.∇nuclear_repulsion(atoms, a)

        # Store kinetic into H and nuclear into S
        GaussianBasis.∇kinetic!(∂H, bset, a)
        GaussianBasis.∇nuclear!(∂S, bset, a)

        ∂H .+= ∂S
        ∂S .= 0

        # Now use S for overlap
        GaussianBasis.∇overlap!(∂S, bset, a)

        idx, xyz... = GaussianBasis.∇sparseERI_2e4c(bset, a)

        for q in 1:3
            @views vH = ∂H[:,:,q]
            AUX .= P .* vH
            ∂E[a, q]  = sum(AUX)

            @views vS = ∂S[:,:,q]
            AUX .= Q .* vS
            ∂E[a, q] -= sum(AUX)

            ∇q = xyz[q]

            for i in eachindex(idx)
                μ,ν,λ,σ = idx[i] .+ 1
                ∇k = ∇q[i]

                if abs(∇k) < 1e-12
                    continue
                end

                # perms is 1 if all indexes are the same (which must not happen)
                # because it implies that the integral is on only one center e.g. (11|11)
                # for which the derivative must be zero
                # In turn, perms will be 2, 4, or 8 depending on the number of permutations
                # allowed. e.g. (11|23) yields 4: (σ != λ) and μν != λσ
                perms = 1 << (μ !== ν) << (σ !== λ) << (Set((μ,ν)) != Set((σ,λ)))

                if perms === 8
                    ∂E[a,q] += (4.0*P[μ,ν]*P[λ,σ] - P[μ,σ]*P[λ,ν] - P[μ,λ]*P[ν,σ])*∇k
                elseif perms === 4
                    ∂E[a,q] += (2.0*P[μ,ν]*P[λ,σ] - 0.5*(P[μ,σ]*P[λ,ν] + P[μ,λ]*P[ν,σ]))*∇k
                else
                    ∂E[a,q] += (P[μ,ν]*P[λ,σ] - 0.25*(P[μ,σ]*P[λ,ν] + P[μ,λ]*P[ν,σ]))*∇k
                end
            end
        end
    end

    return ∂E
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
# Workaround for JuliaLang/julia/pull/28625
if Base.HOME_PROJECT[] !== nothing
  Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

using
  Documenter,
  GeophysicalFlows

makedocs(
 modules = [GeophysicalFlows],
 doctest = false,
   clean = true,
checkdocs = :all,
  format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
 authors = "Navid C. Constantinou and Gregory L. Wagner",
sitename = "GeophysicalFlows.jl",
   pages = Any[
            "Home" => "index.md",
            "Modules" => Any[
              "modules/twodturb.md",
              "modules/barotropicqg.md",
              "modules/multilayerqg.md"
            ],
            "DocStrings" => Any[
            "man/types.md",
            "man/functions.md"]
           ]
)

deploydocs(repo = "github.com/FourierFlows/GeophysicalFlows.jl.git")


text:
baremodule FloatSemantics
    const IEEEHalf   = 0
    const IEEESingle = 1
    const IEEEDouble = 2
    const IEEEQuad   = 3
    const PPCDoubleDouble   = 4
    const x87DoubleExtended = 5
    const bogus = 6
end

baremodule TypeKind
    const void      = 0   # type with no size
    const half      = 1   # 16 bit floating point type
    const float     = 2   # 32 bit floating point type
    const double    = 3   # 64 bit floating point type
    const x86_fp80  = 4   # 80 bit floating point type
    const fp128     = 5   # 128 bit floating point type (112-bit mantissa)
    const ppc_fp128 = 6   # 128 bit floating point type (two 64 bits)
    const label     = 7   # labels
    const integer   = 8   # arbitrary bit width integers
    const func      = 9   # functions
    const struct    = 10  # structures
    const array     = 11  # arrays
    const pointer   = 12  # pointers
    const vector    = 13  # SIMD packed format or other vector type
    const metadata  = 14  # metadata
    const x86_mmx   = 15  # x86 specific multimedia SIMD types
end

baremodule ValueSubclass
    const argument = 0
    const basic_block = 1
    const func = 2
    const global_alias = 3
    const global_variable = 4
    const undef_value = 5
    const block_address = 6
    const const_expr = 7
    const const_aggregate_zero = 8
    const const_data_array = 9
    const const_data_vector = 10
    const const_int = 11
    const const_fp = 12
    const const_array = 13
    const const_struct = 14
    const const_vector = 15
    const const_pointer_null = 16
    const md_node = 17
    const md_string = 18
    const inline_asm = 19
    const pseudo_source_value = 20
    const fixed_stack_pseudo_source_value = 21
    const instruction = 22
end

baremodule Opcode
    # Terminator instructions
    const ret = 1
    const br = 2
    const switch = 3
    const indirectbr = 4
    const invoke = 5
    const unreachable = 7

    # Standard binary operators
    const add = 8
    const fadd = 9
    const sub = 10
    const fsub = 11
    const mul = 12
    const fmul = 13
    const udiv = 14
    const sdiv = 15
    const fdiv = 16
    const urem = 17
    const srem = 18
    const frem = 19

    # Logical operators
    const shl = 20
    const lshr = 21
    const ashr = 22
    const and = 23
    const or = 24
    const xor = 25

    # Memory operators
    const alloca = 26
    const load  = 27
    const store = 28
    const get_element_ptr = 29

    # Cast operators
    const trunc = 30
    const zext = 31
    const sext = 32
    const fptoui = 33
    const fptosi = 34
    const uitofp = 35
    const sitofp = 36
    const fptrunc = 37
    const fpext = 38
    const ptrtoint = 39
    const inttoptr = 40
    const bitcast = 41
    const addrspacecast = 60

    # Other operators
    const icmp = 42
    const fcmp = 43
    const phi = 44
    const call = 45
    const select = 46
    const userop1 = 47
    const userop2 = 48
    const vaarg = 49
    const extract_elem = 50
    const insert_elem = 51
    const shuffle_vector = 52
    const extract_value = 53
    const insert_value = 54

    # Atomic operators
    const fence = 55
    const atomic_cmpxchg = 56
    const atomic_rmw = 57

    # Exception handling operators
    const resume = 58
    const landingpad = 59
end

baremodule IntPred
    const eq  = 32 # equal
    const ne  = 33 # not equal
    const ugt = 34 # unsigned greater than
    const uge = 35 # unsigned greater than or equal
    const ult = 36 # unsinged less than
    const ule = 37 # unsigned less than or equal
    const sgt = 38 # signed greater than
    const sge = 39 # signed greater than or equal
    const slt = 40 # signed less than
    const sle = 41 # signed less than or equal
end

baremodule RealPred
    const predicate_false =  0  # Always false (always folded)
    const oeq =  1    # True if ordered and equal
    const ogt =  2    # True if ordered and greater than
    const oge =  3    # True if ordered and greater than or equal
    const olt =  4    # True if ordered and less than
    const ole =  5    # True if ordered and less than or equal
    const one =  6    # True if ordered and operands are unequal
    const ord =  7    # True if ordered (no nans)
    const uno =  8    # True if unordered: isnan(X) | isnan(Y)
    const ueq =  9    # True if unordered or equal
    const ugt = 10    # True if unordered or greater than
    const uge = 11    # True if unordered, greater than, or equal
    const ult = 12    # True if unordered or less than
    const ule = 13    # True if unordered, less than, or equal
    const une = 14    # True if unordered or not equal
    const predicate_true = 15 # Always true (never folded)
end

baremodule CallingCov
    const c = 0
    const fast = 8
    const cold = 9
    const web_kit_js = 12
    const any_reg = 13
    const x86_stdcall = 64
    const x86_fastcall = 65
end

baremodule Attribute
    const z_ext = 1
    const s_ext = 2
    const no_return = 4
    const in_reg = 8
    const struct_ret = 16
    const no_unwind = 32
    const no_alias = 64
    const by_val = 128
    const nest = 256
    const read_none = 512
    const read_only = 1024
    const no_inline = 2048
    const always_inline = 4096
    const optimize_for_size = 8192
    const stack_protect = 16384
    const stack_protect_req = 32768
    const alignment = 2031616
    const no_capture = 2097152
    const no_red_zone = 4194304
    const no_implicit_float = 8388608
    const naked = 16777216
    const inline_hint = 33554432
    const stack_alignment = 469762048
    const returns_twice = 536870912
    const uw_table = 1073741824
end


text:
@testset "B-coefficients" begin

using SHIPs
using Test, Printf, LinearAlgebra, StaticArrays, BenchmarkTools, Test
using SHIPs: generate_ZKL, generate_ZKL_tuples, SparseSHIP, maxL
using SHIPs.Rotations: ClebschGordan, clebschgordan
using JuLIP.Testing

# function naive_Bcoeff(ll::SVector{4}, mm, cg)
#    coeff = 0.0
#    @assert sum(mm) ≈ 0
#    M = mm[1] + mm[2]
#    for J = max(abs(ll[1]-ll[2]), abs(ll[3]-ll[4])):min(ll[1]+ll[2], ll[3]+ll[4])
#       # @assert clebschgordan(ll[1], mm[1], ll[2], mm[2], J, M) == cg(ll[1], mm[1], ll[2], mm[2], J, M)
#       # @assert clebschgordan(ll[3], mm[3], ll[4], mm[4], J, -M) == cg(ll[3], mm[3], ll[4], mm[4], J, -M)
#       coeff += ( (-1)^(M) *
#                   clebschgordan(ll[1], mm[1], ll[2], mm[2], J, M) *
#                   clebschgordan(ll[3], mm[3], ll[4], mm[4], J, -M) )
#    end
#    return coeff
# end
#
# function naive_Bcoeff(ll::SVector{3}, mm, cg)
#    @assert ( clebschgordan(ll[1], mm[1], ll[2], mm[2], ll[3], -mm[3]) ≈
#               cg(ll[1], mm[1], ll[2], mm[2], ll[3], -mm[3]) )
#    return (-1)^mm[3] * cg(ll[1], mm[1], ll[2], mm[2], ll[3], -mm[3])
# end
#
# function check_Bcoeffs(ll::SVector, cg)
#    allzeros = true
#    for mm in SHIPs._mrange(ll)
#       if abs(mm[end]) > ll[end]; continue; end
#       C = naive_Bcoeff(ll, mm, cg)
#       _C = SHIPs._Bcoeff(ll, mm, cg)
#       if !( (C ≈ _C) || (abs(C-_C) < 1e-12) )
#          return false, allzeros
#       end
#       if !(C ≈ 0);  allzeros = false; end
#    end
#    return true, allzeros
# end

##
@info("Testing _mrange")
ll = SVector(2,3,4)
mrange2 = collect(SHIPs._mrange(ll))
testrg = CartesianIndices( (-2:2, -3:3) )
println(@test all( mpre -> (abs(sum(mpre.I)) > ll[end]) ||
                           (SVector(mpre.I..., -sum(mpre.I)) ∈ mrange2),
                   testrg ))

ll = SVector(4,2,3,4)
mrange3 = collect(SHIPs._mrange(ll))
testrg = CartesianIndices( (-4:4, -2:2, -3:3) )
println(@test all( mpre -> (abs(sum(mpre.I)) > ll[end]) ||
                           (SVector(mpre.I..., -sum(mpre.I)) ∈ mrange3),
                   testrg ))

##
Deg = SparseSHIP(3, 5; wL = 1.0)
Bcoefs = SHIPs.Rotations.CoeffArray()
KL, Nu =  generate_ZKL_tuples(Deg, Bcoefs; filter=false)
KL = KL[1]
_, Nu_filter = generate_ZKL_tuples(Deg, Bcoefs; filter=true)
Izero = Int[]
Iodd = Int[]

Nu3 = Nu[3]
Nu3_filter = Nu_filter[3]

@info("Testing the RI coefficients for Deg = $Deg, 4B")
for (i, zν) in enumerate(Nu3)
   # global Izero, Iodd
   ν = zν.ν
   ll = getfield.(KL[ν], :l)    # SVector([KL[ν[i]].l for i = 1:length(ν)]...)
   # pass, isz = check_Bcoeffs(ll, cg)
   # if isz; push!(Izero, i); end
   Clm = SHIPs.Rotations.single_B(Bcoefs, ll)
   if norm(Clm) == 0; push!(Izero, i); end
   if isodd(sum(ll)); push!(Iodd, i); end
   # print_tf(@test pass)
end
println()
Izodd = union(Izero, Iodd)
@info("""   #Nu = $(length(Nu3)), #Nu_filt = $(length(Nu3_filter)), #(zero or odd) = $(length(Izodd)) """)
println(@test (length(Nu3) == length(Nu3_filter) + length(Izodd)))

##
Deg = SparseSHIP(4, 10; wL = 2.0)
Bcoefs = SHIPs.Rotations.CoeffArray()
KL, Nu =  generate_ZKL_tuples(Deg, Bcoefs; filter=false)
KL = KL[1]
_, Nu_filter = generate_ZKL_tuples(Deg, Bcoefs; filter=true)
Nu4 = Nu[4]
Nu4_filter = Nu_filter[4]
Izero = Int[]
Iodd = Int[]
@info("Testing the RI coefficients for Deg = $Deg, 5B")
for (i, zν) in enumerate(Nu4)
   # global Izero, Iodd
   ν = zν.ν
   ll = SVector([KL[ν[i]].l for i = 1:length(ν)]...)
   # pass, isz = check_Bcoeffs(ll, cg)
   Clm = SHIPs.Rotations.single_B(Bcoefs, ll)
   isz = (norm(Clm) == 0)
   if isz; push!(Izero, i); end
   if isodd(sum(ll)); push!(Iodd, i); end
   # print_tf(@test pass)
end
println()
Izodd = union(Izero, Iodd)
@info("""   #Nu = $(length(Nu4)), #Nu_filt = $(length(Nu4_filter)), #(zero or odd) = $(length(Izodd)) """)
println(@test (length(Nu4) == length(Nu4_filter) + length(Izodd)))

end

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
# # Multiclass classification
#
# As of version `v0.7` and over, the `LSSVClassifier` can handle multiclass classification
# problems.
# It does so by using the so-called _one-vs-one approach_. This is described as follows:
#
# - For a classification problem with ``k`` classes, one trains ``k(k-1)/2`` binary classifiers.
# - Each binary classifier deals with every possible class pair combination. For instance, for a three class classification problem, we need to train three classifiers. The first will deal with the (1, 2) binary classification problem. The second with the (1, 3), while the third classifier will deal with the (2, 3) class pair.
# - We then save all these parameters, and on the prediction step, we ask every classifier to predict a given data instance. Care must be taken to ensure that the class pairs are decoded correctly.
# - A _voting scheme_ is then applied to these predictions. The class that has the largest number of votes when a pediction is carried out is the class that is taken as the predicted class. In the important case of a tie, the predicted class is always that of the lowest index. This is obviously not a great strategy, but it's simple enough to break the tie. See[^1] to read more about this strategy.
#
# Here, we deal with the binary classification problems with the underlying implementation
# of the `LSSVClassifier` and finally pool all of our results together. We then apply
# the voting scheme and we are done.
#
# In this example we will showcase this ability on the very famous and not so difficult
# problem of the [Iris flower dataset.](https://en.wikipedia.org/wiki/Iris_flower_data_set)
#
# This problem is a 3-class classification problem, with 4 features and a total of 150
# samples; 50 for each class.
# It is expected that with a good set of hyperparameters, our implementation should easily
# obtain a 100% accuracy classification rate on this problem.
#
# Let us begin by importing our packages and setting an RNG to ensure reproducibility.
#
using MLJ, MLJBase
using Random
using LeastSquaresSVM

rng = MersenneTwister(957);

# The Iris dataset is so famous, it is included in the `MLJ.jl` package. By using the
# very convenient macro `@load_iris` we can obtain the full dataset, as follows.
X, y = @load_iris;

# Then, we split the dataset. We have said before that the datset is small, so we will
# keep this in mind. A 60-40 train-test split should suffice.
train, test = partition(eachindex(y), 0.6, shuffle=true, rng=rng);

# We will now instantiate our model. Recall that the default kernel is the RBF kernel.
model = LSSVClassifier();

# Now we set some ranges for the hyperparameters and also instantiate a self-tuning model
# Note however that a `repeats` option is now used. What this does is to repeat the
# resampling of the dataset that many times.
# This is just to ensure that the model is learning correctly and not being biased.
r1 = range(model, :σ, lower=1e-2, upper=1e-1);
r2 = range(model, :γ, lower=140, upper=150);
self_tuning_model = TunedModel(
    model=model,
    tuning=Grid(goal=200, rng=rng),
    resampling=StratifiedCV(nfolds=5),
    range=[r1, r2],
    measure=accuracy,
    acceleration=CPUThreads(),
    repeats=10 # Add more resampling to be sure we are not biased
);

# We will use a pipeline to first standardize the data and then feed it to the model.
pipe = @pipeline(Standardizer(), self_tuning_model);
mach = MLJ.machine(pipe, X, y);

# It's time to train. Remeber that if you are going to use multithreading for the
# hyperparameter search, you need to start your `Julia` session with the environment
# variable `JULIA_NUM_THREADS` set to a different number than 1.
MLJ.fit!(mach, rows=train, verbosity=0);

# We want to see the best set of hyperparameters.
fitted_params(mach).deterministic_tuned_model.best_model

# Great! Let us see how good was the model. Let's check with the test set.
results = MLJ.predict(mach, rows=test);
acc = MLJ.accuracy(results, y[test]);
acc

# Fantastic! A 100% accuracy, as was expected.
# As a good measure, let us inspect the confusion matrix of the classification problem.
results = coerce(results, OrderedFactor);
y_ordered = coerce(y[test], OrderedFactor);
confusion_matrix(results, y_ordered)

# ## Conclusions
#
# Again, just as expected, the result is perfect. But do note that this does not mean that
# the implementation or the model is the best there is. There are other models and
# implementations that should be checked along with this one.
#
# As a good rule of thumb, when doing machine learning, always try to train _at least_
# one more model to compare performance.
#
# [^1]: Chih-Wei Hsu and Chih-Jen Lin (2002) ‘A comparison of methods for multiclass support vector machines’, IEEE Transactions on Neural Networks, 13(2), pp. 416. doi: 10.1109/72.991427.
#

max_stars_repo_name: gmagannaDevelop/Elysivm

max_issues_repo_name: gmagannaDevelop/Elysivm

max_forks_repo_name: gmagannaDevelop/Elysivm


text:
"""
Mag2Dpoly

A module to perform magnetic anomaly calculations for 2D polygonal bodies.

# Exports

$(EXPORTS)
"""
module Mag2Dpoly

using DocStringExtensions

export BodySegments2D,MagPolygBodies2D,MagnetizVector
export tmagpolybodies2D,tmagpolybodies2Dgen
export tmagpoly2D,tmagpoly2Dgen

include("magutils.jl")
include("magdatastruct.jl")
include("mag2dpolybodies.jl")

end # module


text:
module KLB

using FileIO

export klbheader, loadklb, loadklb!

add_format(format"KLB", (), ".klb")

function FileIO.load(
    file::File{format"KLB"}
    ;
    bounds::Array{UInt32} = UInt32[0],
    numthreads::Int = Sys.CPU_CORES,
    inplace::AbstractArray = UInt8[0],
    nochecks::Bool = false
    )
    if length(bounds) != 1
        if length(inplace) != 1
            loadklb!(inplace, filename(file), bounds[:,1], bounds[:,2], numthreads=numthreads, nochecks=nochecks)
        else
            loadklb(filename(file), bounds[:,1], bounds[:,2], numthreads=numthreads)
        end
    else
        if length(inplace) != 1
            loadklb!(inplace, filename(file), numthreads=numthreads, nochecks=nochecks)
        else
            loadklb(filename(file), numthreads=numthreads)
        end
    end
end

function FileIO.save(
    file::File{format"KLB"},
    A::AbstractArray
    ;
    numthreads::Int = Sys.CPU_CORES
    )
    writearray(filename(file), A, numthreads)
end

function klbheader( file::AbstractString )
    imagesize = zeros(UInt32, 5)
    blocksize = zeros(UInt32, 5)
    pixelspacing = ones(Float32, 5)
    ktype = Ref{Cint}(-1)
    compressiontype = Ref{Cint}(-1)
    metadata = repeat(" ", 256)

    errid = ccall( (:readKLBheader, "klb"), Cint,
        (Cstring, Ptr{UInt32}, Ref{Cint}, Ptr{Float32}, Ptr{UInt32}, Ref{Cint}, Ptr{Cchar}),
        file, imagesize, ktype, pixelspacing, blocksize, compressiontype, metadata)

    if errid != 0
        error("Could not read KLB header of file '$file'. Error code $errid")
    end

    Dict{AbstractString, Any}(
        "imagesize" => round(Int, imagesize),
        "blocksize" => round(Int, blocksize),
        "pixelspacing" => pixelspacing,
        "metadata" => strip(metadata),
        "datatype" => juliatype(ktype[]),
        "compressiontype" => compressiontype[],
        "spatialorder" => ["x", "y", "z"],
        "colordim" => 4,
        "timedim" => 5)
end

function loadklb(
    file::AbstractString
    ;
    numthreads::Int = Sys.CPU_CORES
    )
    header = klbheader(file)
    A = Array(header["datatype"], header["imagesize"]...)
    loadklb!(A, file, numthreads=numthreads, nochecks=true)
    A
end

function loadklb!(
    A::AbstractArray,
    file::AbstractString
    ;
    numthreads::Int = Sys.CPU_CORES,
    nochecks::Bool = false
    )
    if !nochecks
        header = klbheader(file)
        assert( header["datatype"] == eltype(A) )
        imagesize = header["imagesize"]
        for d in 1:5
            assert( imagesize[d] == size(A, d) )
        end
    end

    errid = ccall( (:readKLBstackInPlace, "klb"), Cint,
        (Cstring, Ptr{Void}, Ref{Cint}, Cint),
        file, A, Ref{Cint}(0), numthreads)

    if errid != 0
        error("Could not read KLB file '$file'. Error code $errid")
    end
end

function loadklb(
    file::AbstractString,
    lower_bounds::Vector{UInt32},
    upper_bounds::Vector{UInt32}
    ;
    numthreads::Int = Sys.CPU_CORES
    )
    header = klbheader(file)
    roisize = 1 + upper_bounds - lower_bounds
    A = Array(header["datatype"], roisize...)
    loadklb!(A, file, lower_bounds, upper_bounds, numthreads=numthreads, nochecks=true)
    A
end

function loadklb!(
    A::AbstractArray,
    file::AbstractString,
    lower_bounds::Vector{UInt32},
    upper_bounds::Vector{UInt32}
    ;
    numthreads::Int = Sys.CPU_CORES,
    nochecks::Bool = false
    )
    lb = lower_bounds - UInt32(1)
    ub = upper_bounds - UInt32(1)

    if !nochecks
        header = klbheader(file)
        assert( header["datatype"] == eltype(A) )
        roisize = 1 + ub - lb
        for d in 1:5
            assert( roisize[d] == size(A, d) )
        end
    end

    errid = ccall( (:readKLBroiInPlace, "klb"), Cint,
        (Cstring, Ptr{Void}, Ptr{UInt32}, Ptr{UInt32}, Cint),
        file, A, lb, ub, numthreads)

    if errid != 0
        error("Could not read KLB file '$file'. Error code $errid")
    end
    A
end

function writearray(
    filepath::AbstractString,
    A::AbstractArray,
    numthreads::Int = Sys.CPU_CORES
    ;
    pixelspacing = C_NULL,
    blocksize = C_NULL,
    compressiontype::Int = 1,
    metadata = C_NULL
    )
    ktype = klbtype(eltype(A))
    imagesize = UInt32[i for i in size(A)]
    while length(imagesize) < 5
        push!(imagesize, 1)
    end

    errid = ccall( (:writeKLBstack, "klb"), Cint,
        (Ptr{Void}, Cstring, Ptr{UInt32}, Cint, Cint, Ptr{Float32}, Ptr{UInt32}, Cint, Ptr{Cchar}),
        A, filepath, imagesize, ktype, numthreads, pixelspacing, blocksize, compressiontype, metadata)

    if errid != 0
        error("Could not write KLB file '$filepath'. Error code $errid")
    end
end

function juliatype( klbtype::Cint )
    if klbtype == 0
        return UInt8
    elseif klbtype == 1
        return UInt16
    elseif klbtype == 2
        return UInt32
    elseif klbtype == 3
        return UInt64
    elseif klbtype == 4
        return Int8
    elseif klbtype == 5
        return Int16
    elseif klbtype == 6
        return Int32
    elseif klbtype == 7
        return Int64
    elseif klbtype == 8
        return Float32
    elseif klbtype == 9
        return Float64
    end
    error( "Unknown or unsupported data type of KLB array: $klbtype" )
end

function klbtype( juliatype::Type )
    if juliatype == UInt8
        return 0
    elseif juliatype == UInt16
        return 1
    elseif juliatype == UInt32
        return 2
    elseif juliatype == UInt64
        return 3
    elseif juliatype == Int8
        return 4
    elseif juliatype == Int16
        return 5
    elseif juliatype == Int32
        return 6
    elseif juliatype == Int64
        return 7
    elseif juliatype == Float32
        return 8
    elseif juliatype == Float64
        return 9
    end
    error( "Unknown or unsupported data type of KLB array: $juliatype" )
end

end # module