text:
module TensorOperations

using TupleTools
using Strided
using Strided: AbstractStridedView, UnsafeStridedView
using LinearAlgebra
using LinearAlgebra: mul!, BLAS.BlasFloat
using LRUCache
using Requires

# Exports
#---------
# export macro API
export @tensor, @tensoropt, @tensoropt_verbose, @optimalcontractiontree, @notensor, @ncon
export @cutensor

export enable_blas, disable_blas, enable_cache, disable_cache, clear_cache, cachesize

# export function based API
export ncon
export tensorcopy, tensoradd, tensortrace, tensorcontract, tensorproduct, scalar
export tensorcopy!, tensoradd!, tensortrace!, tensorcontract!, tensorproduct!

# Convenient type alias
const IndexTuple{N} = NTuple{N,Int}

# An exception type for reporting errors in the index specificatino
struct IndexError{S<:AbstractString} <: Exception
    msg::S
end

# Index notation
#----------------
@nospecialize
include("indexnotation/verifiers.jl")
include("indexnotation/analyzers.jl")
include("indexnotation/preprocessors.jl")
include("indexnotation/ncon.jl")
include("indexnotation/instantiators.jl")
include("indexnotation/postprocessors.jl")
include("indexnotation/parser.jl")
include("indexnotation/poly.jl")
include("indexnotation/optdata.jl")
include("indexnotation/optimaltree.jl")
include("indexnotation/tensormacros.jl")
include("indexnotation/indexordertree.jl")
@specialize

# Implementations
#-----------------
include("implementation/indices.jl")
include("implementation/tensorcache.jl")
include("implementation/stridedarray.jl")
include("implementation/diagonal.jl")

# Functions
#-----------
include("functions/simple.jl")
include("functions/ncon.jl")
include("functions/inplace.jl")

# Global package settings
#-------------------------
# A switch for enabling/disabling the use of BLAS for tensor contractions
const _use_blas = Ref(true)
use_blas() = _use_blas[]
function disable_blas()
    _use_blas[] = false
    return
end
function enable_blas()
    _use_blas[] = true
    return
end

# A cache for temporaries of tensor contractions
const _use_cache = Ref(true)
use_cache() = _use_cache[]

function default_cache_size()
    return min(1<<32, Int(Sys.total_memory())>>2)
end

# methods used for the cache: see implementation/tensorcache.jl for more info
function memsize end
function similar_from_indices end
function similarstructure_from_indices end

taskid() = convert(UInt, pointer_from_objref(current_task()))

const cache = LRU{Any, Any}(; by = memsize, maxsize = default_cache_size())

"""
    disable_cache()

Disable the cache for further use but does not clear its current contents.
Also see [`clear_cache()`](@ref)
"""
function disable_cache()
    _use_cache[] = false
    return
end

"""
    enable_cache(; maxsize::Int = ..., maxrelsize::Real = ...)

(Re)-enable the cache for further use; set the maximal size `maxsize` (as number of bytes)
or relative size `maxrelsize`, as a fraction between 0 and 1, resulting in
`maxsize = floor(Int, maxrelsize * Sys.total_memory())`. Default value is `maxsize = 2^30` bytes, which amounts to 1 gigabyte of memory.
"""
function enable_cache(; maxsize::Int = -1, maxrelsize::Real = 0.0)
    if maxsize == -1 && maxrelsize == 0.0
        maxsize = default_cache_size()
    elseif maxrelsize > 0
        maxsize = max(maxsize, floor(Int, maxrelsize*Sys.total_memory()))
    else
        @assert maxsize >= 0
    end
    _use_cache[] = true
    resize!(cache; maxsize = maxsize)
    return
end

"""
    clear_cache()

Clear the current contents of the cache.
"""
function clear_cache()
    empty!(cache)
    return
end

"""
    cachesize()

Return the current memory size (in bytes) of all the objects in the cache.
"""
cachesize() = cache.currentsize

# Initialization
#-----------------
function __init__()
    resize!(cache; maxsize = default_cache_size())

    @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin
        if CUDA.functional() && CUDA.has_cutensor()
            const CuArray = CUDA.CuArray
            const CublasFloat = CUDA.CUBLAS.CublasFloat
            const CublasReal = CUDA.CUBLAS.CublasReal
            for s in (:handle, :CuDefaultStream, :CuTensorDescriptor, :cudaDataType,
                    :cutensorContractionDescriptor_t, :cutensorContractionFind_t,
                    :cutensorContractionPlan_t,
                    :CUTENSOR_OP_IDENTITY, :CUTENSOR_OP_CONJ, :CUTENSOR_OP_ADD,
                    :CUTENSOR_ALGO_DEFAULT,  :CUTENSOR_WORKSPACE_RECOMMENDED,
                    :cutensorPermutation, :cutensorElementwiseBinary, :cutensorReduction,
                    :cutensorReductionGetWorkspace, :cutensorComputeType,
                    :cutensorGetAlignmentRequirement, :cutensorInitContractionDescriptor,
                    :cutensorInitContractionFind, :cutensorContractionGetWorkspace,
                    :cutensorInitContractionPlan, :cutensorContraction)
                eval(:(const $s = CUDA.CUTENSOR.$s))
            end
            include("implementation/cuarray.jl")
            @nospecialize
            include("indexnotation/cutensormacros.jl")
            @specialize
        end
    end
end

# Some precompile statements
#----------------------------
function _precompile_()
    AVector = Vector{Any}
    for N = 1:8
        @assert precompile(Tuple{typeof(isperm), NTuple{N,Int}})
    end
    @assert precompile(Tuple{typeof(_intersect), Base.BitArray{1}, Base.BitArray{1}})
    @assert precompile(Tuple{typeof(_intersect), Base.BitSet, Base.BitSet})
    @assert precompile(Tuple{typeof(_intersect), UInt128, UInt128})
    @assert precompile(Tuple{typeof(_intersect), UInt32, UInt32})
    @assert precompile(Tuple{typeof(_intersect), UInt64, UInt64})
    @assert precompile(Tuple{typeof(_isemptyset), Base.BitArray{1}})
    @assert precompile(Tuple{typeof(_isemptyset), Base.BitSet})
    @assert precompile(Tuple{typeof(_isemptyset), UInt128})
    @assert precompile(Tuple{typeof(_isemptyset), UInt32})
    @assert precompile(Tuple{typeof(_isemptyset), UInt64})
    @assert precompile(Tuple{typeof(_ncontree!), AVector, Vector{Vector{Int64}}})
    @assert precompile(Tuple{typeof(_setdiff), Base.BitArray{1}, Base.BitArray{1}})
    @assert precompile(Tuple{typeof(_setdiff), Base.BitSet, Base.BitSet})
    @assert precompile(Tuple{typeof(_setdiff), UInt128, UInt128})
    @assert precompile(Tuple{typeof(_setdiff), UInt32, UInt32})
    @assert precompile(Tuple{typeof(_setdiff), UInt64, UInt64})
    @assert precompile(Tuple{typeof(_union), Base.BitArray{1}, Base.BitArray{1}})
    @assert precompile(Tuple{typeof(_union), Base.BitSet, Base.BitSet})
    @assert precompile(Tuple{typeof(_union), UInt128, UInt128})
    @assert precompile(Tuple{typeof(_union), UInt32, UInt32})
    @assert precompile(Tuple{typeof(_union), UInt64, UInt64})
    @assert precompile(Tuple{typeof(_nconmacro), Int, Int, Int})
    @assert precompile(Tuple{typeof(addcost), Power{:χ, Int64}, Power{:χ, Int64}})
    @assert precompile(Tuple{typeof(degree), Power{:x, Int64}})
    @assert precompile(Tuple{typeof(instantiate_contraction), Int, Int, Expr, Int, AVector, AVector, Int})
    @assert precompile(Tuple{typeof(instantiate_generaltensor), Int, Int, Expr, Int, AVector, AVector, Int})
    @assert precompile(Tuple{typeof(instantiate_linearcombination), Int, Int, Expr, Int, AVector, AVector, Int})
    @assert precompile(Tuple{typeof(instantiate), Int, Int, Expr, Int, AVector, AVector, Int})
    @assert precompile(Tuple{typeof(instantiate), Int, Int, Expr, Int, AVector, AVector})
    @assert precompile(Tuple{typeof(instantiate), Expr, Bool, Expr, Int64, AVector, AVector})
    @assert precompile(Tuple{typeof(instantiate), Nothing, Bool, Expr, Bool, AVector, AVector, Bool})
    @assert precompile(Tuple{typeof(instantiate), Symbol, Bool, Expr, Bool, AVector, AVector})
    @assert precompile(Tuple{typeof(instantiate_eltype), Expr})
    @assert precompile(Tuple{typeof(instantiate_scalar), Expr})
    @assert precompile(Tuple{typeof(instantiate_scalar), Float64})
    @assert precompile(Tuple{typeof(disable_blas)})
    @assert precompile(Tuple{typeof(disable_cache)})
    @assert precompile(Tuple{typeof(enable_blas)})
    @assert precompile(Tuple{typeof(enable_cache)})
    @assert precompile(Tuple{typeof(expandconj), Expr})
    @assert precompile(Tuple{typeof(expandconj), Symbol})
    @assert precompile(Tuple{typeof(getallindices), Expr})
    @assert precompile(Tuple{typeof(getallindices), Int})
    @assert precompile(Tuple{typeof(getallindices), Symbol})
    @assert precompile(Tuple{typeof(getindices), Symbol})
    @assert precompile(Tuple{typeof(getindices), Expr})
    @assert precompile(Tuple{typeof(gettensorobject), Int})
    @assert precompile(Tuple{typeof(getlhs), Expr})
    @assert precompile(Tuple{typeof(getrhs), Expr})
    @assert precompile(Tuple{typeof(isindex), Expr})
    @assert precompile(Tuple{typeof(isindex), Symbol})
    @assert precompile(Tuple{typeof(isindex), Int})
    @assert precompile(Tuple{typeof(hastraceindices), Expr})
    @assert precompile(Tuple{typeof(isassignment), Expr})
    @assert precompile(Tuple{typeof(isdefinition), Expr})
    @assert precompile(Tuple{typeof(isgeneraltensor), Expr})
    @assert precompile(Tuple{typeof(istensor), Expr})
    @assert precompile(Tuple{typeof(istensorexpr), Expr})
    @assert precompile(Tuple{typeof(isnconstyle), Array{AVector, 1}})
    @assert precompile(Tuple{typeof(isscalarexpr), Expr})
    @assert precompile(Tuple{typeof(isscalarexpr), Float64})
    @assert precompile(Tuple{typeof(isscalarexpr), LineNumberNode})
    @assert precompile(Tuple{typeof(isscalarexpr), Symbol})
    @assert precompile(Tuple{typeof(istensorexpr), Expr})
    @assert precompile(Tuple{typeof(isgeneraltensor), Expr})
    @assert precompile(Tuple{typeof(decomposetensor), Expr})
    @assert precompile(Tuple{typeof(normalizeindex), Int})
    @assert precompile(Tuple{typeof(normalizeindex), Symbol})
    @assert precompile(Tuple{typeof(normalizeindex), Expr})
    @assert precompile(Tuple{typeof(mulcost), Power{:χ, Int64}, Power{:χ, Int64}})
    @assert precompile(Tuple{typeof(ncontree), Vector{AVector}})
    @assert precompile(Tuple{typeof(optdata), Expr})
    @assert precompile(Tuple{typeof(optdata), Expr, Expr})
    @assert precompile(Tuple{typeof(optimaltree), Vector{AVector}, Base.Dict{Any, Power{:χ, Int64}}})
    @assert precompile(Tuple{typeof(parsecost), Expr})
    @assert precompile(Tuple{typeof(parsecost), Int64})
    @assert precompile(Tuple{typeof(parsecost), Symbol})
    @assert precompile(Tuple{typeof(storeset), Type{Base.BitArray{1}}, AVector, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{Base.BitArray{1}}, Array{Int64, 1}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{Base.BitArray{1}}, Base.Set{Int64}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{Base.BitSet}, AVector, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{Base.BitSet}, Array{Int64, 1}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{Base.BitSet}, Base.Set{Int64}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt128}, AVector, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt128}, Array{Int64, 1}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt128}, Base.Set{Int64}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt32}, AVector, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt32}, Array{Int64, 1}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt32}, Base.Set{Int64}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt64}, AVector, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt64}, Array{Int64, 1}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt64}, Base.Set{Int64}, Int64})
    @assert precompile(Tuple{typeof(tensorify), Expr})
    @assert precompile(Tuple{typeof(extracttensorobjects), Any})
    @assert precompile(Tuple{typeof(_flatten), Expr})
    # @assert precompile(Tuple{typeof(processcontractions), Any, Any, Any})
    @assert precompile(Tuple{typeof(defaultparser), Expr})
    @assert precompile(Tuple{typeof(defaultparser), Any})
    @assert precompile(Tuple{typeof(unique2), AVector})
    @assert precompile(Tuple{typeof(unique2), Array{Int64, 1}})
    @assert precompile(Tuple{typeof(use_blas)})
    @assert precompile(Tuple{typeof(use_cache)})
end
_precompile_()

end # module


text:
using Plots, GLVisualize, GeometryTypes, GLAbstraction
pl_size = if !isempty(GLVisualize.get_screens())
    widths(GLVisualize.current_screen())
else
    (800, 500)
end

description = """
Showing off different line types and marker types,
like 3D meshes, characters and images.
"""

glvisualize(size=pl_size)

ys = Vector[[sin(i*7)/x * cos(x^i) for x in linspace(0.6,5,10)] for i=1:5]
cat = loadasset("cat.obj")
cat = rotationmatrix_y(rad2deg(110)) * cat
marker = Any[:circle, :d, cat, loadasset("foxy.png"), '🐱']
marker = reshape(marker, (1, 5))
colors = reshape(colormap("Blues", 10)[4:8], (5, 1))

plot(
    ys,
    # color = colors,
    line = ([:dot :dash :dashdot :solid :dot], [1 2 3 4 5]),
    marker = marker,
    markersize = [5 10 20 10 15],
    # markercolor = reshape(colormap("RdBu", 5), (5, 1)),
    markerstrokewidth = [1 0.1 0 2 4],
    markerstrokecolor = [:blue :blue :black :blue :white],
)

gui()

ext: jl


text:
include("Solution.jl")
include("Params.jl")

""" Descent Heuristic
1. Select a starting solution.
2. Check the change in the value of the objective function
     for all pairwise exchanges of facilities.
3. If an improving exchange is found,
     the best improving exchange is executed and we go to Step 2.
4. If no improving exchange is found, the algorithm terminates.

Local Search algorithm
S = Initial Solution
While not Terminated
    Explore( N(S) ) ;
    If there is no better neighbor in N(s) Then Stop ;
    S = Select( N(S) ) ;
End While
Return Final solution found (local optima)
"""
function descentHeuristic(sol::Solution; max_iter::Float64=Inf)::Solution
    improved::Bool = true
    num_iter = 0
    while ((improved == true) && (num_iter < max_iter))
        neighborhood = generateSwapNeighbors(sol) # Generate N(S)
        improved = false
        for neighbor in neighborhood    # Explore( N(S) )
            if (neighbor.cost < sol.cost)
                sol = neighbor
                improved = true
            end
        end
        num_iter += 1
    end
    return sol  # return best solution found
end

function cyclicHeuristic(sol::Solution; max_iter::Float64=Inf)::Solution
    improved = false
    neighborhood = generateCyclicNeighbors(sol) # Generate N(S)

    for neighbor in neighborhood    # Explore( N(S) )
        if (neighbor.cost < sol.cost)
            sol = neighbor
            improved = true
        end
    end
    return sol  # return best solution found
end

# Multi-start local search
function multistartLS(params::Params, num_trials::Int; localSearchFunc::Function=descentHeuristic)
    startTime = time()
    results = Dict{String,Any}()
    bestSol = Solution(params) # Random solution
    bestSol = localSearchFunc(bestSol) # Local search
    trial = 1
    while ((trial <= num_trials) && (time() - startTime < params.maxTime) && (bestSol.cost != params.solutionCost))
        sol = Solution(params) # Random solution
        sol = localSearchFunc(sol) # Local Search
        if (sol.cost < bestSol.cost)
            bestSol = sol
        end
        trial += 1
    end
    # Results
    results["time"] = time() - startTime
    results["finalCost"] = bestSol.cost
    results["permutation"] = bestSol.permutation
    results["numTrials"] = trial
    if isnan(params.solutionCost)
        results["gap"] = NaN
    else
        results["gap"] = (bestSol.cost - params.solutionCost) / params.solutionCost
    end
    return bestSol, results
end

max_stars_repo_name: BSAraujo/QAP-Metaheuristics

max_issues_repo_name: BSAraujo/QAP-Metaheuristics

max_forks_repo_name: BSAraujo/QAP-Metaheuristics


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

function niScope_init(resourceName, idQuery, resetDevice, vi)
    ccall((:niScope_init, libniscope), ViStatus, (ViRsrc, ViBoolean, ViBoolean, Ptr{ViSession}), resourceName, idQuery, resetDevice, vi)
end

function niScope_InitWithOptions(resourceName, idQuery, resetDevice, optionString, vi)
    ccall((:niScope_InitWithOptions, libniscope), ViStatus, (ViRsrc, ViBoolean, ViBoolean, ViConstString, Ptr{ViSession}), resourceName, idQuery, resetDevice, optionString, vi)
end

function niScope_close(vi)
    ccall((:niScope_close, libniscope), ViStatus, (ViSession,), vi)
end

function niScope_AutoSetup(vi)
    ccall((:niScope_AutoSetup, libniscope), ViStatus, (ViSession,), vi)
end

function niScope_ConfigureVertical(vi, channelList, vrange, offset, coupling, probeAttenuation, enabled)
    ccall((:niScope_ConfigureVertical, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViReal64, ViInt32, ViReal64, ViBoolean), vi, channelList, vrange, offset, coupling, probeAttenuation, enabled)
end

function niScope_ConfigureChanCharacteristics(vi, channelList, inputImpedance, maxInputFrequency)
    ccall((:niScope_ConfigureChanCharacteristics, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViReal64), vi, channelList, inputImpedance, maxInputFrequency)
end

function niScope_ConfigureHorizontalTiming(vi, minSampleRate, minNumPts, refPosition, numRecords, enforceRealtime)
    ccall((:niScope_ConfigureHorizontalTiming, libniscope), ViStatus, (ViSession, ViReal64, ViInt32, ViReal64, ViInt32, ViBoolean), vi, minSampleRate, minNumPts, refPosition, numRecords, enforceRealtime)
end

function niScope_ConfigureClock(vi, inputClockSource, outputClockSource, clockSyncPulseSource, masterEnabled)
    ccall((:niScope_ConfigureClock, libniscope), ViStatus, (ViSession, ViConstString, ViConstString, ViConstString, ViBoolean), vi, inputClockSource, outputClockSource, clockSyncPulseSource, masterEnabled)
end

function niScope_ExportSignal(vi, signal, signalIdentifier, outputTerminal)
    ccall((:niScope_ExportSignal, libniscope), ViStatus, (ViSession, ViInt32, ViConstString, ViConstString), vi, signal, signalIdentifier, outputTerminal)
end

function niScope_AdjustSampleClockRelativeDelay(vi, delay)
    ccall((:niScope_AdjustSampleClockRelativeDelay, libniscope), ViStatus, (ViSession, ViReal64), vi, delay)
end

function niScope_ConfigureTriggerEdge(vi, triggerSource, level, slope, triggerCoupling, holdoff, delay)
    ccall((:niScope_ConfigureTriggerEdge, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViInt32, ViInt32, ViReal64, ViReal64), vi, triggerSource, level, slope, triggerCoupling, holdoff, delay)
end

function niScope_ConfigureTriggerGlitch(vi, triggerSource, level, width, polarity, glitchCondition, triggerCoupling, holdoff, delay)
    ccall((:niScope_ConfigureTriggerGlitch, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViReal64, ViInt32, ViInt32, ViInt32, ViReal64, ViReal64), vi, triggerSource, level, width, polarity, glitchCondition, triggerCoupling, holdoff, delay)
end

function niScope_ConfigureTriggerHysteresis(vi, triggerSource, level, hysteresis, slope, triggerCoupling, holdoff, delay)
    ccall((:niScope_ConfigureTriggerHysteresis, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViReal64, ViInt32, ViInt32, ViReal64, ViReal64), vi, triggerSource, level, hysteresis, slope, triggerCoupling, holdoff, delay)
end

function niScope_ConfigureTriggerWindow(vi, triggerSource, lowLevel, highLevel, windowMode, triggerCoupling, holdoff, delay)
    ccall((:niScope_ConfigureTriggerWindow, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViReal64, ViInt32, ViInt32, ViReal64, ViReal64), vi, triggerSource, lowLevel, highLevel, windowMode, triggerCoupling, holdoff, delay)
end

function niScope_ConfigureTriggerSoftware(vi, holdoff, delay)
    ccall((:niScope_ConfigureTriggerSoftware, libniscope), ViStatus, (ViSession, ViReal64, ViReal64), vi, holdoff, delay)
end

function niScope_SendSoftwareTriggerEdge(vi, whichTrigger)
    ccall((:niScope_SendSoftwareTriggerEdge, libniscope), ViStatus, (ViSession, ViInt32), vi, whichTrigger)
end

function niScope_ConfigureTriggerImmediate(vi)
    ccall((:niScope_ConfigureTriggerImmediate, libniscope), ViStatus, (ViSession,), vi)
end

function niScope_ConfigureTriggerRunt(vi, triggerSource, lowThreshold, highThreshold, polarity, triggerCoupling, holdoff, delay)
    ccall((:niScope_ConfigureTriggerRunt, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViReal64, ViInt32, ViInt32, ViReal64, ViReal64), vi, triggerSource, lowThreshold, highThreshold, polarity, triggerCoupling, holdoff, delay)
end

function niScope_ConfigureTriggerDigital(vi, triggerSource, slope, holdoff, delay)
    ccall((:niScope_ConfigureTriggerDigital, libniscope), ViStatus, (ViSession, ViConstString, ViInt32, ViReal64, ViReal64), vi, triggerSource, slope, holdoff, delay)
end

function niScope_ConfigureTriggerVideo(vi, triggerSource, enableDcRestore, signalFormat, event, lineNumber, polarity, triggerCoupling, holdoff, delay)
    ccall((:niScope_ConfigureTriggerVideo, libniscope), ViStatus, (ViSession, ViConstString, ViBoolean, ViInt32, ViInt32, ViInt32, ViInt32, ViInt32, ViReal64, ViReal64), vi, triggerSource, enableDcRestore, signalFormat, event, lineNumber, polarity, triggerCoupling, holdoff, delay)
end

function niScope_ConfigureTriggerWidth(vi, triggerSource, level, lowThreshold, highThreshold, polarity, condition, triggerCoupling, holdoff, delay)
    ccall((:niScope_ConfigureTriggerWidth, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViReal64, ViReal64, ViInt32, ViInt32, ViInt32, ViReal64, ViReal64), vi, triggerSource, level, lowThreshold, highThreshold, polarity, condition, triggerCoupling, holdoff, delay)
end

function niScope_ConfigureEqualizationFilterCoefficients(vi, channelList, numberOfCoefficients, coefficients)
    ccall((:niScope_ConfigureEqualizationFilterCoefficients, libniscope), ViStatus, (ViSession, ViConstString, ViInt32, Ptr{ViReal64}), vi, channelList, numberOfCoefficients, coefficients)
end

function niScope_GetEqualizationFilterCoefficients(vi, channel, numberOfCoefficients, coefficients)
    ccall((:niScope_GetEqualizationFilterCoefficients, libniscope), ViStatus, (ViSession, ViConstString, ViInt32, Ptr{ViReal64}), vi, channel, numberOfCoefficients, coefficients)
end

function niScope_GetFrequencyResponse(vi, channel, bufferSize, frequencies, amplitudes, phases, numberOfFrequencies)
    ccall((:niScope_GetFrequencyResponse, libniscope), ViStatus, (ViSession, ViConstString, ViInt32, Ptr{ViReal64}, Ptr{ViReal64}, Ptr{ViReal64}, Ptr{ViInt32}), vi, channel, bufferSize, frequencies, amplitudes, phases, numberOfFrequencies)
end

function niScope_ConfigureAcquisition(vi, acquisitionType)
    ccall((:niScope_ConfigureAcquisition, libniscope), ViStatus, (ViSession, ViInt32), vi, acquisitionType)
end

function niScope_InitiateAcquisition(vi)
    ccall((:niScope_InitiateAcquisition, libniscope), ViStatus, (ViSession,), vi)
end

function niScope_Abort(vi)
    ccall((:niScope_Abort, libniscope), ViStatus, (ViSession,), vi)
end

function niScope_Commit(vi)
    ccall((:niScope_Commit, libniscope), ViStatus, (ViSession,), vi)
end

function niScope_Read(vi, channelList, timeout, numSamples, waveform, wfmInfo)
    ccall((:niScope_Read, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViInt32, Ptr{ViReal64}, Ptr{niScope_wfmInfo}), vi, channelList, timeout, numSamples, waveform, wfmInfo)
end

function niScope_Fetch(vi, channelList, timeout, numSamples, waveform, wfmInfo)
    ccall((:niScope_Fetch, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViInt32, Ptr{ViReal64}, Ptr{niScope_wfmInfo}), vi, channelList, timeout, numSamples, waveform, wfmInfo)
end

function niScope_FetchBinary8(vi, channelList, timeout, numSamples, waveform, wfmInfo)
    ccall((:niScope_FetchBinary8, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViInt32, Ptr{ViInt8}, Ptr{niScope_wfmInfo}), vi, channelList, timeout, numSamples, waveform, wfmInfo)
end

function niScope_FetchBinary16(vi, channelList, timeout, numSamples, waveform, wfmInfo)
    ccall((:niScope_FetchBinary16, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViInt32, Ptr{ViInt16}, Ptr{niScope_wfmInfo}), vi, channelList, timeout, numSamples, waveform, wfmInfo)
end

function niScope_FetchBinary32(vi, channelList, timeout, numSamples, waveform, wfmInfo)
    ccall((:niScope_FetchBinary32, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViInt32, Ptr{ViInt32}, Ptr{niScope_wfmInfo}), vi, channelList, timeout, numSamples, waveform, wfmInfo)
end

function niScope_FetchComplex(vi, channelList, timeout, numSamples, wfm, wfmInfo)
    ccall((:niScope_FetchComplex, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViInt32, Ptr{NIComplexNumber}, Ptr{niScope_wfmInfo}), vi, channelList, timeout, numSamples, wfm, wfmInfo)
end

function niScope_FetchComplexBinary16(vi, channelList, timeout, numSamples, wfm, wfmInfo)
    ccall((:niScope_FetchComplexBinary16, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViInt32, Ptr{NIComplexI16}, Ptr{niScope_wfmInfo}), vi, channelList, timeout, numSamples, wfm, wfmInfo)
end

function niScope_AcquisitionStatus(vi, acquisitionStatus)
    ccall((:niScope_AcquisitionStatus, libniscope), ViStatus, (ViSession, Ptr{ViInt32}), vi, acquisitionStatus)
end

function niScope_ActualNumWfms(vi, channelList, numWfms)
    ccall((:niScope_ActualNumWfms, libniscope), ViStatus, (ViSession, ViConstString, Ptr{ViInt32}), vi, channelList, numWfms)
end

function niScope_ActualMeasWfmSize(vi, arrayMeasFunction, measWaveformSize)
    ccall((:niScope_ActualMeasWfmSize, libniscope), ViStatus, (ViSession, ViInt32, Ptr{ViInt32}), vi, arrayMeasFunction, measWaveformSize)
end

function niScope_ActualRecordLength(vi, recordLength)
    ccall((:niScope_ActualRecordLength, libniscope), ViStatus, (ViSession, Ptr{ViInt32}), vi, recordLength)
end

function niScope_SampleRate(vi, sampleRate)
    ccall((:niScope_SampleRate, libniscope), ViStatus, (ViSession, Ptr{ViReal64}), vi, sampleRate)
end

function niScope_SampleMode(vi, sampleMode)
    ccall((:niScope_SampleMode, libniscope), ViStatus, (ViSession, Ptr{ViInt32}), vi, sampleMode)
end

function niScope_AddWaveformProcessing(vi, channelList, measFunction)
    ccall((:niScope_AddWaveformProcessing, libniscope), ViStatus, (ViSession, ViConstString, ViInt32), vi, channelList, measFunction)
end

function niScope_ClearWaveformProcessing(vi, channelList)
    ccall((:niScope_ClearWaveformProcessing, libniscope), ViStatus, (ViSession, ViConstString), vi, channelList)
end

function niScope_ClearWaveformMeasurementStats(vi, channelList, clearableMeasurementFunction)
    ccall((:niScope_ClearWaveformMeasurementStats, libniscope), ViStatus, (ViSession, ViConstString, ViInt32), vi, channelList, clearableMeasurementFunction)
end

function niScope_ReadMeasurement(vi, channelList, timeout, scalarMeasFunction, result)
    ccall((:niScope_ReadMeasurement, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViInt32, Ptr{ViReal64}), vi, channelList, timeout, scalarMeasFunction, result)
end

function niScope_FetchMeasurement(vi, channelList, timeout, scalarMeasFunction, result)
    ccall((:niScope_FetchMeasurement, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViInt32, Ptr{ViReal64}), vi, channelList, timeout, scalarMeasFunction, result)
end

function niScope_FetchMeasurementStats(vi, channelList, timeout, scalarMeasFunction, result, mean, stdev, minval, maxval, numInStats)
    ccall((:niScope_FetchMeasurementStats, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViInt32, Ptr{ViReal64}, Ptr{ViReal64}, Ptr{ViReal64}, Ptr{ViReal64}, Ptr{ViReal64}, Ptr{ViInt32}), vi, channelList, timeout, scalarMeasFunction, result, mean, stdev, minval, maxval, numInStats)
end

function niScope_FetchArrayMeasurement(vi, channelList, timeout, arrayMeasFunction, measWfmSize, measWfm, wfmInfo)
    ccall((:niScope_FetchArrayMeasurement, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViInt32, ViInt32, Ptr{ViReal64}, Ptr{niScope_wfmInfo}), vi, channelList, timeout, arrayMeasFunction, measWfmSize, measWfm, wfmInfo)
end

function niScope_reset(vi)
    ccall((:niScope_reset, libniscope), ViStatus, (ViSession,), vi)
end

function niScope_self_test(vi, selfTestResult, selfTestMessage)
    ccall((:niScope_self_test, libniscope), ViStatus, (ViSession, Ptr{ViInt16}, Ptr{ViChar}), vi, selfTestResult, selfTestMessage)
end

function niScope_Disable(vi)
    ccall((:niScope_Disable, libniscope), ViStatus, (ViSession,), vi)
end

function niScope_ResetDevice(vi)
    ccall((:niScope_ResetDevice, libniscope), ViStatus, (ViSession,), vi)
end

function niScope_CalSelfCalibrate(vi, channelList, option)
    ccall((:niScope_CalSelfCalibrate, libniscope), ViStatus, (ViSession, ViConstString, ViInt32), vi, channelList, option)
end

function niScope_revision_query(vi, driverRevision, firmwareRevision)
    ccall((:niScope_revision_query, libniscope), ViStatus, (ViSession, Ptr{ViChar}, Ptr{ViChar}), vi, driverRevision, firmwareRevision)
end

function niScope_ProbeCompensationSignalStart(vi)
    ccall((:niScope_ProbeCompensationSignalStart, libniscope), ViStatus, (ViSession,), vi)
end

function niScope_ProbeCompensationSignalStop(vi)
    ccall((:niScope_ProbeCompensationSignalStop, libniscope), ViStatus, (ViSession,), vi)
end

function niScope_CableSenseSignalStart(vi)
    ccall((:niScope_CableSenseSignalStart, libniscope), ViStatus, (ViSession,), vi)
end

function niScope_CableSenseSignalStop(vi)
    ccall((:niScope_CableSenseSignalStop, libniscope), ViStatus, (ViSession,), vi)
end

function niScope_IsDeviceReady(resourceName, channelList, deviceReady)
    ccall((:niScope_IsDeviceReady, libniscope), ViStatus, (ViRsrc, ViConstString, Ptr{ViBoolean}), resourceName, channelList, deviceReady)
end

function niScope_GetChannelName(vi, index, bufferSize, channelString)
    ccall((:niScope_GetChannelName, libniscope), ViStatus, (ViSession, ViInt32, ViInt32, Ptr{ViChar}), vi, index, bufferSize, channelString)
end

function niScope_GetChannelNameFromString(vi, index, bufferSize, name)
    ccall((:niScope_GetChannelNameFromString, libniscope), ViStatus, (ViSession, ViConstString, ViInt32, Ptr{ViChar}), vi, index, bufferSize, name)
end

function niScope_errorHandler(vi, errorCode, errorSource, errorDescription)
    ccall((:niScope_errorHandler, libniscope), ViStatus, (ViSession, ViStatus, Ptr{ViChar}, Ptr{ViChar}), vi, errorCode, errorSource, errorDescription)
end

function niScope_GetError(vi, errorCode, bufferSize, description)
    ccall((:niScope_GetError, libniscope), ViStatus, (ViSession, Ptr{ViStatus}, ViInt32, Ptr{ViChar}), vi, errorCode, bufferSize, description)
end

function niScope_GetErrorMessage(vi, errorCode, bufferSize, errorMessage)
    ccall((:niScope_GetErrorMessage, libniscope), ViStatus, (ViSession, ViStatus, ViInt32, Ptr{ViChar}), vi, errorCode, bufferSize, errorMessage)
end

function niScope_LockSession(vi, callerHasLock)
    ccall((:niScope_LockSession, libniscope), ViStatus, (ViSession, Ptr{ViBoolean}), vi, callerHasLock)
end

function niScope_UnlockSession(vi, callerHasLock)
    ccall((:niScope_UnlockSession, libniscope), ViStatus, (ViSession, Ptr{ViBoolean}), vi, callerHasLock)
end

function niScope_GetAttributeViInt32(vi, channelList, attributeId, value)
    ccall((:niScope_GetAttributeViInt32, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViInt32}), vi, channelList, attributeId, value)
end

function niScope_SetAttributeViInt32(vi, channelList, attributeId, value)
    ccall((:niScope_SetAttributeViInt32, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViInt32), vi, channelList, attributeId, value)
end

function niScope_CheckAttributeViInt32(vi, channelList, attributeId, value)
    ccall((:niScope_CheckAttributeViInt32, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViInt32), vi, channelList, attributeId, value)
end

function niScope_GetAttributeViInt64(vi, channelList, attributeId, value)
    ccall((:niScope_GetAttributeViInt64, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViInt64}), vi, channelList, attributeId, value)
end

function niScope_SetAttributeViInt64(vi, channelList, attributeId, value)
    ccall((:niScope_SetAttributeViInt64, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViInt64), vi, channelList, attributeId, value)
end

function niScope_CheckAttributeViInt64(vi, channelList, attributeId, value)
    ccall((:niScope_CheckAttributeViInt64, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViInt64), vi, channelList, attributeId, value)
end

function niScope_GetAttributeViReal64(vi, channelList, attributeId, value)
    ccall((:niScope_GetAttributeViReal64, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViReal64}), vi, channelList, attributeId, value)
end

function niScope_SetAttributeViReal64(vi, channelList, attributeId, value)
    ccall((:niScope_SetAttributeViReal64, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViReal64), vi, channelList, attributeId, value)
end

function niScope_CheckAttributeViReal64(vi, channelList, attributeId, value)
    ccall((:niScope_CheckAttributeViReal64, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViReal64), vi, channelList, attributeId, value)
end

function niScope_GetAttributeViString(vi, channelList, attributeId, bufSize, value)
    ccall((:niScope_GetAttributeViString, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViInt32, Ptr{ViChar}), vi, channelList, attributeId, bufSize, value)
end

function niScope_SetAttributeViString(vi, channelList, attributeId, value)
    ccall((:niScope_SetAttributeViString, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViConstString), vi, channelList, attributeId, value)
end

function niScope_CheckAttributeViString(vi, channelList, attributeId, value)
    ccall((:niScope_CheckAttributeViString, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViConstString), vi, channelList, attributeId, value)
end

function niScope_GetAttributeViSession(vi, channelList, attributeId, value)
    ccall((:niScope_GetAttributeViSession, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViSession}), vi, channelList, attributeId, value)
end

function niScope_SetAttributeViSession(vi, channelList, attributeId, value)
    ccall((:niScope_SetAttributeViSession, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViSession), vi, channelList, attributeId, value)
end

function niScope_CheckAttributeViSession(vi, channelList, attributeId, value)
    ccall((:niScope_CheckAttributeViSession, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViSession), vi, channelList, attributeId, value)
end

function niScope_GetAttributeViBoolean(vi, channelList, attributeId, value)
    ccall((:niScope_GetAttributeViBoolean, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViBoolean}), vi, channelList, attributeId, value)
end

function niScope_SetAttributeViBoolean(vi, channelList, attributeId, value)
    ccall((:niScope_SetAttributeViBoolean, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViBoolean), vi, channelList, attributeId, value)
end

function niScope_CheckAttributeViBoolean(vi, channelList, attributeId, value)
    ccall((:niScope_CheckAttributeViBoolean, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViBoolean), vi, channelList, attributeId, value)
end

function niScope_ResetAttribute(vi, channelList, attributeId)
    ccall((:niScope_ResetAttribute, libniscope), ViStatus, (ViSession, ViConstString, ViAttr), vi, channelList, attributeId)
end

function niScope_ResetAllAttributes(vi)
    ccall((:niScope_ResetAllAttributes, libniscope), ViStatus, (ViSession,), vi)
end

function niScope_ImportAttributeConfigurationBuffer(vi, sizeInBytes, configuration)
    ccall((:niScope_ImportAttributeConfigurationBuffer, libniscope), ViStatus, (ViSession, ViInt32, Ptr{ViInt8}), vi, sizeInBytes, configuration)
end

function niScope_ExportAttributeConfigurationBuffer(vi, sizeInBytes, configuration)
    ccall((:niScope_ExportAttributeConfigurationBuffer, libniscope), ViStatus, (ViSession, ViInt32, Ptr{ViInt8}), vi, sizeInBytes, configuration)
end

function niScope_ImportAttributeConfigurationFile(vi, filePath)
    ccall((:niScope_ImportAttributeConfigurationFile, libniscope), ViStatus, (ViSession, ViConstString), vi, filePath)
end

function niScope_ExportAttributeConfigurationFile(vi, filePath)
    ccall((:niScope_ExportAttributeConfigurationFile, libniscope), ViStatus, (ViSession, ViConstString), vi, filePath)
end

function niScope_GetScalingCoefficients(vi, channelList, bufferSize, coefficientInfo, numberOfCoefficientSets)
    ccall((:niScope_GetScalingCoefficients, libniscope), ViStatus, (ViSession, ViConstString, ViInt32, Ptr{niScope_coefficientInfo}, Ptr{ViInt32}), vi, channelList, bufferSize, coefficientInfo, numberOfCoefficientSets)
end

function niScope_GetNormalizationCoefficients(vi, channelList, bufferSize, coefficientInfo, numberOfCoefficientSets)
    ccall((:niScope_GetNormalizationCoefficients, libniscope), ViStatus, (ViSession, ViConstString, ViInt32, Ptr{niScope_coefficientInfo}, Ptr{ViInt32}), vi, channelList, bufferSize, coefficientInfo, numberOfCoefficientSets)
end

function niScope_GetStreamEndpointHandle(vi, streamName, writerHandle)
    ccall((:niScope_GetStreamEndpointHandle, libniscope), ViStatus, (ViSession, ViConstString, Ptr{ViUInt32}), vi, streamName, writerHandle)
end


text:
@testset "Symmetric matrices" begin
    m = [
        1 2 1
        2 3 2
        1 2 1
    ]

    @test unit(m) == m
    @test rotate90(m) == m
    @test rotate180(m) == m
    @test rotate270(m) == m
    @test flipdiag(m) == m
    @test flipadiag(m) == m
    @test flipx(m) == m
    @test flipy(m) == m
end

@testset "Non-symmetric matrices" begin
    m = [
        1 2 3
        4 5 6
        7 8 9
    ]

    @test unit(m) == m
    @test rotate90(m) == [
        7 4 1
        8 5 2
        9 6 3
    ]
    @test rotate180(m) == [
        9 8 7
        6 5 4
        3 2 1
    ]
    @test rotate180(m) == [
        9 8 7
        6 5 4
        3 2 1
    ]
    @test rotate270(m) == [
        3 6 9
        2 5 8
        1 4 7
    ]
    @test flipdiag(m) == [
        1 4 7
        2 5 8
        3 6 9
    ]
    @test flipadiag(m) == [
        9 6 3
        8 5 2
        7 4 1
    ]
    @test flipx(m) == [
        7 8 9
        4 5 6
        1 2 3
    ]
    @test flipy(m) == [
        3 2 1
        6 5 4
        9 8 7
    ]
end

@testset "Group operation" begin
    # the '∘' function acts as function composition and group operation at the same time
    # it can be used prevent unneeded computations, e.g., (rotate180 ∘ rotate180)(m) does
    # not actually rotate m twice by 180 degrees but directly returns unit(m)
    m = rand(10, 10)

    for g1 in SquareSymmetries.D4
        for g2 in SquareSymmetries.D4
            @test (g1 ∘ g2)(m) == g1(g2(m))
        end
    end
end

@testset "Inverse elements" begin
    m = rand(10, 10)

    # for all group elements g and all matrices m, g(g^(-1)(m)) = g^(-1)(g(m)) = m
    for g in SquareSymmetries.D4
        @test g(inv(g)(m)) == inv(g)(g(m)) == m
    end

    for g in SquareSymmetries.D4
        @test g ∘ inv(g) == inv(g) ∘ g == unit
    end
end

@testset "Symmetries" begin
    m = rand(10, 10)

    syms = symmetries(m)
    @test length(syms) == 8
    @test syms[1] == unit(m)
    @test syms[2] == rotate90(m)
    @test syms[3] == rotate180(m)
    @test syms[4] == rotate270(m)
    @test syms[5] == flipx(m)
    @test syms[6] == flipy(m)
    @test syms[7] == flipdiag(m)
    @test syms[8] == flipadiag(m)
end


text:
function pp62(data::Dict)

  # Setup basic dimensions of arrays

  # Parse & check FEdict data

  if :struc_el in keys(data)
    struc_el = data[:struc_el]
    @assert typeof(struc_el) <: StructuralElement
  else
    println("No fin_el type specified.")
    return
  end

  ndim = struc_el.ndim
  nst = struc_el.nst

  if struc_el.axisymmetric
    nst = 4
  end

  fin_el = struc_el.fin_el
  @assert typeof(fin_el) <: FiniteElement

  if typeof(fin_el) == Line
    (nels, nn) = mesh_size(fin_el, struc_el.nxe)
  elseif typeof(fin_el) == Triangle || typeof(fin_el) == Quadrilateral
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye)
  elseif typeof(fin_el) == Hexahedron
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye, struc_el.nze)
  else
    println("$(typeof(fin_el)) is not a known finite element.")
    return
  end

  nodof = fin_el.nodof           # Degrees of freedom per node
  ndof = fin_el.nod * nodof      # Degrees of freedom per fin_el

  # Update penalty if specified in FEdict

  penalty = 1e20
  if :penalty in keys(data)
    penalty = data[:penalty]
  end

  tol = 1.0e-10
  if :tol in keys(data)
    tol = copy(data[:tol])
  end

  limit = 100
  if :limit in keys(data)
    limit = copy(data[:limit])
  end

  if :cg_tol in keys(data)
    cg_tol = data[:cg_tol]
  else
    cg_tol = 1.0e-5
  end

  if :cg_limit in keys(data)
    cg_limit = data[:cg_limit]
  else
    cg_limit = Int(1e3)
  end

  # Allocate all arrays

  # Start with arrays to be initialized from FEdict

  if :properties in keys(data)
    prop = zeros(size(data[:properties], 1), size(data[:properties], 2))
    for i in 1:size(data[:properties], 1)
      prop[i, :] = data[:properties][i, :]
    end
  else
    println("No :properties key found in FEdict")
  end

  nf = ones(Int, nodof, nn)
  if :support in keys(data)
    for i in 1:size(data[:support], 1)
      nf[:, data[:support][i][1]] = data[:support][i][2]
    end
  end

  x_coords = zeros(nn)
  if :x_coords in keys(data)
    x_coords = data[:x_coords]
  end

  y_coords = zeros(nn)
  if :y_coords in keys(data)
    y_coords = data[:y_coords]
  end

  z_coords = zeros(nn)
  if :z_coords in keys(data)
    z_coords = data[:z_coords]
  end

  g_coord = zeros(ndim,nn)
  if :g_coord in keys(data)
    g_coord = data[:g_coord]'
  end

  g_num = zeros(Int, fin_el.nod, nels)
  if :g_num in keys(data)
    g_num = reshape(data[:g_num]', fin_el.nod, nels)
  end

  etype = ones(Int, nels)
  if :etype in keys(data)
    etype = data[:etype]
  end

  # All other arrays

  points = zeros(struc_el.nip, ndim)
  g = zeros(Int, ndof)
  fun = zeros(fin_el.nod)
  coord = zeros(fin_el.nod, ndim)
  gamma = zeros(nels)
  jac = zeros(ndim, ndim)
  der = zeros(ndim, fin_el.nod)
  deriv = zeros(ndim, fin_el.nod)
  bee = zeros(nst,ndof)
  km = zeros(ndof, ndof)
  mm = zeros(ndof, ndof)
  gm = zeros(ndof, ndof)
  kg = zeros(ndof, ndof)
  eld = zeros(ndof)
  weights = zeros(struc_el.nip)
  g_g = zeros(Int, ndof, nels)
  num = zeros(Int, fin_el.nod)
  actions = zeros(ndof, nels)
  displacements = zeros(size(nf, 1), ndim)
  gc = ones(ndim)
  dee = zeros(nst,nst)
  sigma = zeros(nst)
  axial = zeros(nels)

  eps = zeros(nst)
  stress = zeros(nst)
  bload =zeros(ndof)
  eload =zeros(ndof)
  erate = zeros(nst)
  evp = zeros(nst)
  devp = zeros(nst)
  m1 = zeros(nst, nst)
  m2 = zeros(nst, nst)
  m3 = zeros(nst, nst)
  flow = zeros(nst, nst)

  formnf!(nodof, nn, nf)
  neq = maximum(nf)

  @assert :qincs in keys(data)
  qincs = data[:qincs]

  loads = zeros(Float64, neq+1)
  bdylds = zeros(Float64, neq+1)
  oldis = zeros(Float64, neq+1)
  totd = zeros(Float64, neq+1)
  p = zeros(nst, struc_el.nip, nels)
  x = zeros(neq+1)
  u = zeros(neq+1)
  xnew = zeros(neq+1)
  diag_precon = zeros(neq+1)
  d = zeros(neq+1)

  # Find global array sizes
  for iel in 1:nels
    geom_rect!(fin_el, iel, x_coords, y_coords, coord, num, struc_el.direction)
    num_to_g!(num, nf, g)
    g_num[:, iel] = num
    g_coord[:, num] = coord'
    g_g[:, iel] = g
  end

  tensor = zeros(nst, struc_el.nip, nels)
  storkm = zeros(ndof, ndof, nels)
  dt = 1.0e15
  ddt = 0.0
  sigm = 0.0
  dsbar = 0.0
  lode_theta = 0.0

  sample!(fin_el, points, weights)

  # element stiffnes integration, storage and conditioner

  for iel in 1:nels
    ddt = 4.0*(1.0+prop[etype[iel], 3])/(3.0*prop[etype[iel], 2])
    if ddt < dt
      dt = ddt
    end
    deemat!(dee, prop[etype[iel], 2], prop[etype[iel], 3])
    num = g_num[:, iel]
    coord = g_coord[:, num]'              # Transpose
    g = g_g[:, iel]
    km = zeros(ndof, ndof)
    for i in 1:struc_el.nip
      shape_der!(der, points, i)
      jac = der*coord
      detm = det(jac)
      jac = inv(jac)
      deriv = jac*der
      beemat!(bee, deriv)
      km += (bee')*dee*bee*detm*weights[i]
    end
    storkm[:,:,iel] = km
    for k in 1:ndof
      diag_precon[g[k]+1] = diag_precon[g[k]+1] + km[k,k]
    end
  end

  diag_precon = 1 ./ diag_precon

  # Read load weightings

  loaded_nodes = 0
  node = Int[]
  val = Array{Float64, 2}
  if :loaded_nodes in keys(data)
    loaded_nodes = size(data[:loaded_nodes], 1)
    node = zeros(Int, loaded_nodes)
    val = zeros(loaded_nodes, size(data[:loaded_nodes][1][2], 2))
  end

  if :loaded_nodes in keys(data)
    for i in 1:loaded_nodes
      node[i] = data[:loaded_nodes][i][1]
      val[i,:] = data[:loaded_nodes][i][2]
    end
  end

  nf1 = deepcopy(nf) + 1

  # Load increment loop

  converged = false
  ptot = 0.0
  iters = 0
  cg_tot = 0
  diag_precon[1] = 0.0
  iy = 0

  for iy in 1:size(qincs, 1)
    ptot = ptot + qincs[iy]
    iters = 0
    bdylds = zeros(Float64, neq+1)
    evpt = zeros(nst, struc_el.nip, nels)
    cg_tot = 0
    cg_converged = false

    # Plastic iteration loop

    while true
      iters += 1
      loads = zeros(Float64, neq+1)
      for i in 1:loaded_nodes
        loads[nf1[:,node[i]]] = val[i,:] * qincs[iy]
      end
      loads += bdylds
      d = diag_precon .* loads
      p = d
      x = zeros(neq+1)
      cg_iters = 0

      # Pcg equation solution
      tmpconvcrit = 0.0

      while true
        cg_iters += 1
        u = zeros(neq+1)

        for iel in 1:nels
          g = g_g[:, iel]
          km = storkm[:, :, iel]
          u[g+1] += km * p[g+1]
        end
        up = dot(loads, d)
        alpha = up ./ dot(p, u)
        xnew = x + p .* alpha
        loads -= u .* alpha
        d = diag_precon .* loads
        beta = dot(loads, d) ./ up
        p = d + p * beta
        #tmpconvcrit =  maximum(abs.(xnew-x))/maximum(abs.(xnew))
        cg_converged = checon(xnew, x, cg_tol)
        #iy < 2 && println([iters cg_iters cg_tot+cg_iters tmpconvcrit])
        if cg_converged || (cg_iters == cg_limit)
          break
        end
      end
      cg_tot += cg_iters
      loads = xnew
      loads[1] = 0.0

      # Check plastic convergence

      converged = checon(loads, oldis, tol)
      if iters == 1
        converged = false
      end
      if converged || (iters == limit)
        bdylds = zeros(Float64, neq+1)
      end

      # Go round the Gauss point

      for iel in 1:nels
        deemat!(dee, prop[etype[iel], 2], prop[etype[iel], 3])
        num = g_num[:, iel]
        coord = g_coord[:, num]'
        g = g_g[:, iel]
        eld = loads[g+1]
        bload = zeros(ndof)
        for i in 1:struc_el.nip
          shape_der!(der, points, i)
          jac = der*coord
          detm = det(jac)
          jac = inv(jac)
          deriv = jac*der
          beemat!(bee, deriv)
          eps = bee*eld
          eps -= evpt[:, i, iel]
          sigma = dee*eps
          stress = sigma + tensor[:, i, iel]
          (sigm, dsbar, lode_theta) = invar(stress, sigm, dsbar, lode_theta)

          # Check whether yield is violated

          f = dsbar - sqrt(3.0)*prop[etype[iel], 1]
          if converged || iters == limit
            devp = deepcopy(stress)
          else
            if f >= 0.0
              dq1 = 0.0
              dq2 = 3.0/2.0/dsbar
              dq3 = 0.0
              (m1, m2, m3) = formm!(stress, m1, m2, m3)
              flow = f*(m1*dq1 + m2*dq2 + m3*dq3)
              erate = flow*stress
              evp = erate*dt
              evpt[:,i,iel] += evp
              devp = dee*evp
            end
          end
          if f >= 0.0 || converged || iters == limit
            eload = devp' * bee
            bload += eload' .* detm .* weights[i]
          end

          # Update the Gauss point stresses

          if converged || iters == limit
            tensor[:,i,iel] = stress
          end
        end

        # Compute the total bodyloads vector

        bdylds[g+1] += bload
        bdylds[1] = 0.0
      end
      if converged || iters == limit
        break
      end
    end
    totd += loads
    totdstr = @sprintf("%+.4e", totd[nf1[2, node[1]]])

    if iy < 10
      println(" $(iy)       $(ptot)    $(totdstr)    $(iters)       $(round.(cg_tot/iters, digits=2))")
    else
      println("$(iy)       $(ptot)    $(totdstr)    $(iters)       $(round.(cg_tot/iters, digits=2))")
    end

  end
  (ptot, totd[nf1[2, node[1]]], iters, round.(cg_tot/iters, 2))
end

function p62a(data)
  lst = [deepcopy(data) for i in 1:length(data[:qincs])];
  for i in 2:length(data[:qincs])
    lst[i][:qincs] = [sum(data[:qincs][1:(i-1)]), data[:qincs][i]]
  end
  np = nprocs()  # determine the number of processes available
  n = length(lst)
  results = Vector{Any}(n)
  local ix = 2
  # function to produce the next work item from the queue.
  # in this case it's just an index.
  nextidx() = (idx=ix; ix+=1; idx)
  @sync begin
      for p=1:np
          if p != myid() || np == 1
              @async begin
                  while true
                      idx = nextidx()
                      if idx > n
                          break
                      end
                      results[idx] = remotecall_fetch(pp62, p, lst[idx])
                  end
              end
          end
      end
  end
  results[2:end]
end

max_stars_repo_path: src/6 Material Nonlinearity/p62a.jl

max_issues_repo_path: src/6 Material Nonlinearity/p62a.jl

max_forks_repo_path: src/6 Material Nonlinearity/p62a.jl


text:
using FlexDates

using DiscreteRanges: DiscreteRange, (..)
using Test
using Dates: Day, Date

@testset "constructors, equality, arithmetic" begin
    E1 = Date(1980, 1, 1)
    E2 = Date(1990, 1, 1)

    D1 = Date(2001, 1, 8)
    D2 = Date(2007, 2, 9)

    F11 = FlexDate{E1, Int16}(D1)
    F12 = FlexDate{E2, Int32}(D1)

    F21 = FlexDate{E1, Int16}(D2)
    F22 = FlexDate{E2, Int32}(D2)

    F3 = FlexDate{E1}(0)
    F4 = FlexDate{E1}(5)
    F5 = FlexDate{E1}(-5)

    @test F11 == F12
    @test F21 == F22
    @test F11 ≠ F22
    @test F5 < F3 < F4
    @test F5 < E1 < F4
    @test F3 == E1

    @test F5 - F3 == Day(-5)
    @test F3 + Day(5) == F4
    @test F3 - Day(5) == F5

    @test length(F3..F4) == 6
    @test length(F4..F3) == 0

    @test_throws InexactError F11 + Day(typemax(Int16) + 1)
end

@testset "FlexDate hashing" begin
    D = Date(2001, 3, 4)
    @test hash(FlexDate{Date(1980,1,1),Int16}(D)) ==
               hash(FlexDate{Date(1990,2,3),Int32}(D))
end

@testset "interval (Int16)" begin
    FD = FlexDate{Date(2005,1,1), Int16}

    d1 = FD(2005, 1, 1)
    d2 = FD(2008, 3, 7)
    d3 = FD(2009, 8, 11)
    d4 = FD(2013, 12, 1)

    @test d2 ∈ d1..d3
    @test (d1..d3) ∩ (d2..d4) == d2..d3
    @test (d1..d3) ∪ (d2..d4) == d1..d4

    @test repr(d2) == "2008-03-07 [2005-01-01 + Int16 days]"

    @test repr(d1..d3) == "2005-01-01..2009-08-11 [2005-01-01 + Int16 days]"

    @test typemax(FD) == FD(typemax(Int16))
    @test typemin(FD) == FD(typemin(Int16))
end

@testset "promotion" begin
    E1 = Date(1980, 1, 1)
    E2 = Date(1980, 5, 30)
    D1 = Date(1970, 2, 7)
    D2 = Date(1991, 4, 9)
    F11 = FlexDate{E1, Int16}(D1)
    F12 = FlexDate{E1, Int32}(D2)
    @test promote(F11, F12) ≡ (FlexDate{E1, Int32}(D1), F12)
    @test promote(F11, FlexDate{E2,Int16}(D2)) ≡ (F11, FlexDate{E1, Int16}(D2))
    @test promote(F11, FlexDate{E2,Int32}(D2)) ≡ (FlexDate{E1, Int32}(D1), F12)
end


text:
module linearstep

include("adam.jl")

end # module


text:
macro manual_union_split(body::Expr, conditions...)
    if !(body.head === :-> && length(body.args) == 2 && body.args[1] == :(()))
        error(
            "`@manual_union_split` is intended to be used with a `do` blockt",
            " with no argumen",
        )
    end
    body = body.args[2]
    ex = foldr(conditions; init = body) do c, ex
        quote
            if $c
                $body
            else
                $ex
            end
        end
    end
    esc(ex)
end

function ithtype(::Type{T}, i::Val) where {T}
    S = foldrunion(T, Val(1)) do S, j
        if j === i
            S
        else
            if j isa Val
                Val(valueof(j) + 1)
            else
                j
            end
        end
    end
    if S isa Type
        return S
    else
        return Union{}
    end
end


text:
using Hop, Test, LinearAlgebra

function getgraphene()
    lat = [1.0 0.5 0.0; 0.0 (√3)/2 0.0; 0.0 0.0 1.0]
    positions = [1/3 2/3; 1/3 2/3; 0.0 0.0]
    graphene = TightBindingModel(lat, positions)
    sethopping!(graphene, 1, 2, [0, 0, 0], -1.0)
    sethopping!(graphene, 2, 1, [1, 0, 0], -1.0)
    sethopping!(graphene, 2, 1, [0, 1, 0], -1.0)
    return graphene
end

let
    graphene = getgraphene()
    @test graphene.norbits == 2
    @test graphene.rlat ≈ [
        2π 0.0 0.0;
        -2π*tan(π/6) 2π/cos(π/6) 0.0;
        0.0 0.0 2π;
    ]
    @test geteigvals(graphene, [0.0, 0.0, 0.0]) ≈ [-3.0, 3.0]
    @test gethamiltonian(graphene, [0.0, 0.0, 0.0])*geteig(graphene, [0.0, 0.0, 0.0])[2][:, 1] ≈
        -3.0*geteig(graphene, [0.0, 0.0, 0.0])[2][:, 1]
    kp = KPath(graphene.rlat, [1 0 1 0; 0 1 0 1; 0 0 0 0], 3)
    bands = getband(graphene, kp)
    @test kp.distances ≈ [0, 2π, 4π, 4π, 6π, 8π]
    @test bands[:, 1:3] ≈ [-3.0 -1.0 -3.0; 3.0 1.0 3.0]
end

let    # Kane-Mele
    lat = [1 0.5 0; 0 (√3)/2 0; 0 0 1]
    positions = [1/3 2/3; 1/3 2/3; 0 0]
    tm = TightBindingModel(lat, positions, spinful=true)
    σ0 = [1 0; 0 1]
    σ1 = [0 1; 1 0]
    σ2 = [0 -im; im 0]
    σ3 = [1 0; 0 -1]
    onsite = 1.0
    t = 1.0
    so = 0.6*t*0.5
    sethopping!(tm, 1, 1, [0, 0, 0], σ0*onsite)
    sethopping!(tm, 2, 2, [0, 0, 0], -σ0*onsite)
    sethopping!(tm, 1, 2, [0, 0, 0], σ0*t)
    sethopping!(tm, 1, 2, [0, -1, 0], σ0*t)
    sethopping!(tm, 1, 2, [-1, 0, 0], σ0*t)
    sethopping!(tm, 1, 1, [0, 1, 0], -im*so*σ3)
    sethopping!(tm, 1, 1, [1, 0, 0], im*so*σ3)
    sethopping!(tm, 1, 1, [1, -1, 0], -im*so*σ3)
    sethopping!(tm, 2, 2, [0, 1, 0], im*so*σ3)
    sethopping!(tm, 2, 2, [1, 0, 0], -im*so*σ3)
    sethopping!(tm, 2, 2, [1, -1, 0], im*so*σ3)
    @test geteigvals(tm, [2/3, 1/3, 0.0]) ≈ [-2.55885, -0.558846, 0.558846, 2.55885] atol=1.0e-5
end

include("floquet.jl")


text:
using CUDA

import CUDA: CuArray

togpudevice(x) = gpu(x)
togpudevice(x, xs...) = (todevice(x), map(todevice, xs)...)
togpudevice(x::Union{Tuple, NamedTuple}) = map(todevice, x)
togpudevice(x::AbstractArray{<:Union{Int, <: NTuple{N, Int} where N}}) = CUDA.CuArray(x)

include("./statistic.jl")
include("./batch_gemm_gpu.jl")
include("./batch_tril_gpu.jl")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ErrorResponse2 <: SwaggerModel
    error::Any # spec type: Union{ Nothing, ErrorDetails } # spec name: error

    function ErrorResponse2(;error=nothing)
        o = new()
        validate_property(ErrorResponse2, Symbol("error"), error)
        setfield!(o, Symbol("error"), error)
        o
    end
end # type ErrorResponse2

const _property_map_ErrorResponse2 = Dict{Symbol,Symbol}(Symbol("error")=>Symbol("error"))
const _property_types_ErrorResponse2 = Dict{Symbol,String}(Symbol("error")=>"ErrorDetails")
Base.propertynames(::Type{ ErrorResponse2 }) = collect(keys(_property_map_ErrorResponse2))
Swagger.property_type(::Type{ ErrorResponse2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ErrorResponse2[name]))}
Swagger.field_name(::Type{ ErrorResponse2 }, property_name::Symbol) =  _property_map_ErrorResponse2[property_name]

function check_required(o::ErrorResponse2)
    true
end

function validate_property(::Type{ ErrorResponse2 }, name::Symbol, val)
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to crop strings.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export left_crop, fit_string_in_field, get_crop_to_fit_string_in_field
export right_crop

"""
    left_crop(str::AbstractString, crop_width::Int)

Left crop from `str` a field with printable width `crop_width`.

This function return two strings:

- The ANSI escape sequence (non-printable string) in the cropped field; and
- The cropped string.
"""
function left_crop(str::AbstractString, crop_width::Int)
    buf_ansi = IOBuffer()
    buf_str = IOBuffer(sizehint = floor(Int, sizeof(str) - crop_width))
    state = :text

    for c in str
        if crop_width ≤ 0
            write(buf_str, c)
        else
            state = _process_string_state(c, state)

            if state == :text
                crop_width -= textwidth(c)

                # If `crop_width` is negative, then it means that we have a
                # character that occupies more than 1 character. In this case,
                # we fill the string with space.
                if crop_width < 0
                    write(buf_str, " "^(-crop_width))
                    crop_width = 0
                end
            else
                write(buf_ansi, c)
            end
        end
    end

    return String(take!(buf_ansi)), String(take!(buf_str))
end

"""
    fit_string_in_field(str::AbstractString, field_width::Int; kwargs...)

Crop the string `str` to fit it in a field with width `field_width`.

# Keywords

- `add_continuation_char::Bool`: If `true`, a continuation character is added to
    the cropped string end. Notice that the returned string always has a
    printable width of `field_size`. (**Default** = `true`)
- `add_space_in_continuation_char::Bool`: If `true`, a space is added before the
    continuation character if `crop_size` is `:right`, or after the continuation
    character if `crop_side` if `:left`. (**Default** = `false`)
- `continuation_char::Char`: The continuation character to add if the string is
    cropped. (**Default** = `…`)
- `crop_side::Symbol`: Select from which side the characters must be removed to
    fit the string into the field. It can be `:right` or `:left`.
    (**Default** = `:right`)
- `field_margin::Int`: Consider an additional margin in the field if it must be
    cropped. (**Default** = 0)
- `keep_ansi::Bool`: If `true`, the ANSI escape sequences found in the cropped
    part will be kept. (**Default** = `true`)
- `printable_string_width::Int`: Provide the printable string width to reduce
    the computational burden. If this parameters is lower than 0, the printable
    width is compute internally. (**Default** = -1)
"""
function fit_string_in_field(
    str::AbstractString,
    field_width::Int;
    add_continuation_char::Bool = true,
    add_space_in_continuation_char::Bool = false,
    continuation_char::Char = '…',
    crop_side::Symbol = :right,
    field_margin::Int = 0,
    keep_ansi::Bool = true,
    printable_string_width::Int = -1
)
    str_width = printable_string_width < 0 ?
        printable_textwidth(str) :
        printable_string_width

    crop = get_crop_to_fit_string_in_field(
        str,
        field_width - field_margin;
        add_continuation_char,
        add_space_in_continuation_char,
        continuation_char,
        printable_string_width = str_width
    )

    if crop ≤ field_margin
        return str
    end

    cont_str = add_continuation_char ? string(continuation_char) : ""

    # Crop from the right
    # ==========================================================================

    if crop_side == :right
        cropped_str, ansi = right_crop(
            str,
            crop;
            keep_escape_seq = keep_ansi,
            printable_string_width = str_width
        )

        if add_space_in_continuation_char
            return cropped_str * " " * cont_str * ansi
        else
            return cropped_str * cont_str * ansi
        end

    # Crop from the left
    # ==========================================================================

    else
        ansi, cropped_str = left_crop(str, crop)

        result = keep_ansi ? ansi : ""

        if add_space_in_continuation_char
            result *= cont_str * " " * cropped_str
        else
            result *= cont_str * cropped_str
        end

        return result
    end
end

"""
    get_crop_to_fit_string_in_field(str::AbstractString, field_width::Int; kwargs...)

Get the number of printable characters to be cropped so that the string `str`
can fit in the field with size `field_width`.

# Keywords

- `add_continuation_char::Bool`: If `true`, a continuation character is added to
    the cropped string end. Notice that the returned string always has a
    printable width of `field_size`. (**Default** = `true`)
- `add_space_in_continuation_char::Bool`: If `true`, a space is added before the
    continuation character if `crop_size` is `:right`, or after the continuation
    character if `crop_side` if `:left`. (**Default** = `false`)
- `continuation_char::Char`: The continuation character to add if the string is
    cropped. (**Default** = `…`)
- `printable_string_width::Int`: Provide the printable string width to reduce
    the computational burden. If this parameters is lower than 0, the printable
    width is compute internally. (**Default** = -1)
"""
function get_crop_to_fit_string_in_field(
    str::AbstractString,
    field_width::Int;
    add_continuation_char::Bool = true,
    add_space_in_continuation_char::Bool = false,
    continuation_char::Char = '…',
    printable_string_width::Int = -1
)

    str_width = printable_string_width < 0 ?
        printable_textwidth(str) :
        printable_string_width

    Δ = str_width - field_width

    # If the field is larger than the string, then we do not need to crop.
    if Δ ≤ 0
        return 0
    end

    # If the user is asking for the continuation char, then we must crop the
    # string to account for the continuation char.
    cont_str = ""

    if add_continuation_char
        cont_str_width = textwidth(continuation_char)

        if add_space_in_continuation_char
            cont_str_width += 1
        end

        Δ += cont_str_width

        # If we are left with no space, then we must crop the entire string.
        if Δ > str_width
            return str_width
        end
    end

    return Δ
end

"""
    right_crop(str::AbstractString, crop_width::Int; kwargs...)

Right crop from `str` a field with printable width `crop_width`.

This function return two strings:

- The cropped string; and
- The ANSI escape sequence (non-printable string) in the cropped field.

# Keyword

- `keep_escape_seq::Bool`: If `false`, then the ANSI escape sequence in the
    cropped field will not be computed. In this case, the second argument
    returned is always empty. (**Default** = `true`)
- `printable_string_width::Int`: Provide the printable string width to reduce
    the computational burden. If this parameters is lower than 0, the printable
    width is compute internally. (**Default** = -1)

!!! note
    If the keyword `keep_escape_seq` is set to `true`, then all the string must
    be processed, which can lead to a substantial increase in computational
    time.
"""
function right_crop(
    str::AbstractString,
    crop_width::Int;
    keep_escape_seq::Bool = true,
    printable_string_width::Int = -1
)
    buf_ansi = IOBuffer()
    buf_str = IOBuffer(sizehint = floor(Int, max(0, sizeof(str) - crop_width)))

    str_width = printable_string_width < 0 ?
        printable_textwidth(str) :
        printable_string_width

    remaining_chars = str_width - crop_width
    state = :text

    for c in str
        state = _process_string_state(c, state)

        if remaining_chars ≤ 0
            !keep_escape_seq && break
            state != :text && write(buf_ansi, c)
        else
            if state == :text
                Δ = textwidth(c)
                remaining_chars -= Δ

                # If `remaining_chars` is negative, then it means that we have a
                # character that occupies more than 1 character. In this case,
                # we fill the string with space.
                if remaining_chars < 0
                    write(buf_str, " "^(-remaining_chars))
                    remaining_chars = 0
                else
                    write(buf_str, c)
                end
            else
                write(buf_str, c)
            end
        end
    end

    return String(take!(buf_str)), String(take!(buf_ansi))
end


text:
indices = 0:49
rates = [784]
d = 50
max_iter = 100

fw_arguments = (
    line_search=FrankWolfe.MonotonousNonConvexStepSize(),
    max_iteration=max_iter,
    print_iter=max_iter / 10,
    verbose=true,
    #lazy=true,
)


text:
"""
$(SIGNATURES)

SecondaryCodeOrBitDetector to detect the secondary code or bits shift.
"""
struct SecondaryCodeOrBitDetector
    buffer::UInt64
    length::Int
    found::Bool
end

function SecondaryCodeOrBitDetector()
    SecondaryCodeOrBitDetector(0, 0, false)
end

@inline get_buffer(sc_bit_detector::SecondaryCodeOrBitDetector) = sc_bit_detector.buffer
@inline found(sc_bit_detector::SecondaryCodeOrBitDetector) = sc_bit_detector.found
@inline length(sc_bit_detector::SecondaryCodeOrBitDetector) = sc_bit_detector.length

"""
$(SIGNATURES)

Finds the secondary code or the bit shift.
"""
function find(::Type{S}, detector, prompt_correlator) where S <: AbstractGNSSSystem
    found(detector) && return detector
    prns = get_buffer(detector) << 1 + UInt64(real(prompt_correlator) > 0)
    num_prns = length(detector) + 1
    bit_found = is_upcoming_integration_new_bit(S, prns, num_prns)
    SecondaryCodeOrBitDetector(prns, num_prns, bit_found)
end


text:
###################################
## Example code
## for Bayesian Gaussin Mixture Model

using PyPlot, PyCall
push!(LOAD_PATH,".")
import GaussianMixtureModel

"""
Visualize data & estimation in 2D space.
"""
function visualize_2D(X::Matrix{Float64}, S::Matrix{Float64}, S_est::Matrix{Float64}, text)
    cmp = get_cmap("jet")

    K1 = size(S, 1)
    K2 = size(S_est, 1)
    col1 = [pycall(cmp.o, PyAny, Int(round(val)))[1:3] for val in linspace(0,255,K1)]
    col2 = [pycall(cmp.o, PyAny, Int(round(val)))[1:3] for val in linspace(0,255,K2)]

    f, (ax1, ax2) = subplots(1,2,num=text)
    f[:clf]()
    f, (ax1, ax2) = subplots(1,2,num=text)

    for k in 1 : K1
        ax1[:scatter](X[1, S[k,:].==1], X[2, S[k,:].==1], color=col1[k])
    end
    ax1[:set_title]("truth")

    for k in 1 : K2
        ax2[:scatter](X[1, S_est[k,:].==1], X[2, S_est[k,:].==1], color=col2[k])
    end

    ax2[:set_title]("estimation")
end

"""
Run a test script for 2D data clustering.
"""
function test_2D()
    ## set model
    D = 2 # data dimension
    K = 4 #  number of mixture components
    alpha = 100.0 * ones(K)
    beta = 0.1
    m = zeros(D)
    nu = D + 1.0
    W = eye(D)
    cmp = [GaussianMixtureModel.GW(beta, m, nu, W) for _ in 1 : K]
    bgmm = GaussianMixtureModel.BGMM(D, K, alpha, cmp)

    ## generate data
    N = 300
    gmm = GaussianMixtureModel.sample_GMM(bgmm)
    X, S = GaussianMixtureModel.sample_data(gmm, N)

    ## inference
    max_iter = 100
    tic()
    S_est, post_bgmm, VB = GaussianMixtureModel.learn_VI(X, bgmm, max_iter)
    #S_est, post_bgmm, VB = GaussianMixtureModel.learn_GS(X, bgmm, max_iter)
    #S_est, post_bgmm, VB = GaussianMixtureModel.learn_CGS(X, bgmm, max_iter)
    toc()

    ## plot
    visualize_2D(X, S, GaussianMixtureModel.winner_takes_all(S_est), "2D plot")

    # VB check
    figure("ELBO")
    clf()
    plot(VB)
    ylabel("ELBO")
    xlabel("iterations")
    show()
end

test_2D()


text:
function newton!(iters::Vector{Int64},Q::Vector{Float64},A::Vector{Float64},
    children::Vector{Int64},Qold::Vector{Float64},Aold::Vector{Float64},W::Vector{Float64},
    beta::Vector{Float64},A0::Vector{Float64},c0::Vector{Float64},rho::Float64,
    f::Function,J::Function,maxiter::Int16,maxval::Float64,epsJ::Float64,epsN::Float64)

    numchildren = length(children);

    # initial guess state vector (flow rate/area pairs)
    x0 = zeros(2*numchildren+2);
    xx = zeros(2*numchildren+2);
    x0[1] = Qold[1];
    x0[2] = Aold[1];
    for i = 1:numchildren
        x0[2*i+1] = Qold[i+1];
        x0[2*i+2] = Aold[i+1];
    end

    # setup for iterations
    xx .= x0;
    x = zeros(length(x0));
    N = 1; # iteration counter
    xn = zeros(length(x0));

    # allocation for objective function/Jacobian
    if numchildren == 1
        ff = zeros(4);
        JJ = zeros(4,4);
    elseif numchildren == 2
        ff = zeros(6);
        JJ = zeros(6,6);
    elseif numchildren == 3
        ff = zeros(8);
        JJ = zeros(8,8);
    end

    # Newton iterations
    while N <= maxiter
        # determine Jacobian, check invertibility
        # println("Interior Jacobian:")
        J(JJ,xx,beta,rho);
        if abs(det(JJ)) < epsJ
            error("Newton Jacobian is singular.");
        end
        # update state vector
        # println("State update:")
        f(ff,xx,beta,A0,rho,c0,W);
        xn .= xx .- inv(JJ)*ff
        # check if sufficiently close to root
        # println("Root checking:")
        f(ff,xn,beta,A0,rho,c0,W);
        if norm(ff) <= epsN
            x = xn;
            # println(x)
            iters[1] += N;
            break
        end
        # check for divergence
        f(ff,xn,beta,A0,rho,c0,W)
        if norm(ff) >= maxval
            iters[1] += N;
            error("Newton iteration diverged.");
        end
        N+=1;
        xx = xn;
        if N == maxiter
            error("Newton iteration failed to converge.");
        end
    end

    # update junction using converged state vector
    # println("Junction update:")
    Q[1] = x[1];
    A[1] = x[2];
    for i = 1:numchildren
        Q[i+1] = x[2*i+1];
        A[i+1] = x[2*i+2];
    end
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
@testset "BatchMvNormal" begin

    # for the BatchMvNormal tests to work BatchMvNormals have to be functors!
    Flux.@functor ConditionalDists.BatchDiagMvNormal
    Flux.@functor ConditionalDists.BatchScalMvNormal

    # BatchDiagMvNormal
    xlen = 10
    batch = 3
    rtol = atol = 1e-6
    μs = rand(Float32, xlen, batch)
    σs = rand(Float32, xlen, batch)
    xs = rand(Float32, xlen, batch)

    d = BatchMvNormal(μs,σs)
    @test eltype(d) == eltype(μs)
    @test Distributions.params(d) == (μs,σs)
    @test all(mean(d) .== μs)
    @test all(var(d) .== σs .^2)

    llh = logpdf(d|>gpu, xs|>gpu) |> cpu
    testds  = [MvNormal(m,s) for (m,s) in zip(eachcol(μs), eachcol(σs))]
    testllh = Float32.([logpdf(p,x) for (p,x) in zip(testds,eachcol(xs))])
    @test all(isapprox.(llh, testllh, rtol=rtol, atol=atol))

    f(x, μ, σ) = sum(logpdf(BatchMvNormal(μ, σ) |> gpu, x |> gpu))
    @test_nowarn Flux.gradient(f, xs, μs, σs)

    f(μ, σ) = sum(rand(BatchMvNormal(μ,σ) |> gpu))
    @test_nowarn Flux.gradient(f, μs, σs)

    # BatchScalMvNormal
    μs = rand(Float32, xlen, batch)
    ss = rand(Float32, batch)
    xs = rand(Float32, xlen, batch)
    d = BatchMvNormal(μs,ss)

    llh = logpdf(d|>gpu,xs|>gpu) |> cpu
    testds  = [MvNormal(m,s) for (m,s) in zip(eachcol(μs),ss)]
    testllh = Float32.([logpdf(p,x) for (p,x) in zip(testds,eachcol(xs))])
    @test all(isapprox.(llh, testllh, rtol=rtol, atol=atol))

    f(x, μ, σ) = sum(logpdf(BatchMvNormal(μ|>gpu, σ|>gpu), x|>gpu))
    @test_nowarn Flux.gradient(f, xs, μs, ss)

    f(μ, σ) = sum(rand(BatchMvNormal(μ,σ)|>gpu))
    @test_nowarn Flux.gradient(f, μs, ss)

end


text:
using Gadgetron
using Test

@testset "AcquisitionSerialization" begin

    acq = MRD.Acquisition(MRD.AcquisitionHeader(),zeros(ComplexF32,192,12))

    io = IOBuffer()
    MRD.write(io,acq )
    seekstart(io)
    acq2 = MRD.read(io,MRD.Acquisition)

    @test acq.header == acq2.header
    @test acq.data == acq2.data
    @test acq == acq2
end


text:
#=
using Test
using LatticeTools

using Combinatorics
using LinearAlgebra
using YAML

@testset "Group" begin
    @testset "Exceptions" begin
        @test_throws ArgumentError FiniteGroup([1 1; 1 2])
        @test_throws ArgumentError FiniteGroup([1 2 3; 2 1 1; 3 1 1])
        @test_throws ArgumentError FiniteGroup([2 3; 3 2])

        # Right Bol Loop 8.1.4.0 of order 8
        # http://ericmoorhouse.org/pub/bol/htmlfiles8/8_1_4_0.html
        @test_throws ArgumentError FiniteGroup([
            1 2 3 4 5 6 7 8;
            2 8 6 1 7 3 5 4;
            3 7 8 6 1 4 2 5;
            4 1 7 8 6 5 3 2;
            5 6 1 7 8 2 4 3;
            6 3 4 5 2 8 1 7;
            7 5 2 3 4 1 8 6;
            8 4 5 2 3 7 6 1;
        ])
    end
    @testset "FiniteGroup-Abelian" begin
        @test_throws ArgumentError FiniteGroup([1 1 1; 1 1 1])
        @test_throws ArgumentError FiniteGroup([1 1 1; 1 1 1; 1 1 1])

        # Example: Z₃
        mtab = [
            1 2 3;
            2 3 1;
            3 1 2;
        ]
        group = FiniteGroup(mtab)

        @test group_order(group) == 3
        @test group_order(group, 1) == 1
        @test group_order(group, 2) == 3
        @test group_order(group, 3) == 3
        @test group_order(group, [1,2,3]) == [1,3,3]
        @test period_length(group, 1) == 1
        @test period_length(group, 2) == 3
        @test period_length(group, 3) == 3
        @test period_length(group, [1,2,3]) == [1,3,3]
        @test isabelian(group)
        @test group_multiplication_table(group) == mtab
        @test element(group, 2) == 2
        @test element(group, 1:2) == 1:2
        @test elements(group) == 1:3
        @test element_name(group, 2) == "2"
        @test element_name(group, 1:2) == ["1", "2"]
        @test element_names(group) == ["1", "2", "3"]
        @test_throws BoundsError element(group, 5)
        @test_throws BoundsError element_name(group, 5)

        gp = group_product(group)
        for i in 1:3, j in 1:3
            @test group_product(group, i, j) == mtab[i, j]
            @test gp(i, j) == mtab[i, j]
        end
        @test group_product(group, 2, BitSet([1,2])) == BitSet([2,3])
        @test group_product(group, BitSet([1,2]), 2) == BitSet([2,3])
        @test group_product(group, BitSet([1,2]), BitSet([1,2])) == BitSet([1,2,3])

        @test gp(2, BitSet([1,2])) == BitSet([2,3])
        @test gp(BitSet([1,2]), 2) == BitSet([2,3])
        @test gp(BitSet([1,2]), BitSet([1,2])) == BitSet([1,2,3])

        @test group_inverse(group, 1) == 1
        @test group_inverse(group, 2) == 3
        @test group_inverse(group, 3) == 2
        @test group_inverse(group, [1,2]) == [1,3]

        ginv = group_inverse(group)
        @test ginv(1) == 1
        @test ginv(2) == 3
        @test ginv(3) == 2
        @test ginv([1,2]) == [1,3]

        @test generate_subgroup(group, 1) == BitSet([1])
        @test generate_subgroup(group, 2) == BitSet([1,2,3])
        @test generate_subgroup(group, [1,2]) == BitSet([1,2,3])

        @test issubgroup(group, Set([1]))
        @test !issubgroup(group, Set([1,2]))
        @test minimal_generating_set(group) == [2]

        @test group_multiplication_table([[1 0; 0 1], [1 0; 0 -1]]) == [1 2; 2 1]

        @test ishomomorphic(1:3, group; product=gp)
        @test !ishomomorphic(1:2, group; product=gp)
    end

    @testset "FiniteGroup-Nonabelian" begin
        # C3v, element#2 = C3, element#4 = σᵥₐ
        # (non-abelian)
        group = FiniteGroup([
            1 2 3 4 5 6;
            2 3 1 6 4 5;
            3 1 2 5 6 4;
            4 5 6 1 2 3;
            5 6 4 3 1 2;
            6 4 5 2 3 1;
        ])
        generators = minimal_generating_set(group)
        @test generators == [2, 4]
        @test generate_subgroup(group, generators) == BitSet(1:6) # completely generates

        ϕ = [1, 3, 4, 5, 2, 6] # group isomorphism
        mtab1 = group_multiplication_table(group)
        mtab2 = zeros(Int, (6,6))
        for x in 1:6, y in 1:6
            # ϕ(x)⋅ϕ(y) = ϕ(x⋅y)
            mtab2[ϕ[x], ϕ[y]] = ϕ[mtab1[x,y]]
        end
        group2 = FiniteGroup(mtab2)

        # ϕ: group  →  group2
        #       x   ↦  ϕ(x)
        ϕ2 = group_isomorphism(group, group2)
        mtab3 = zeros(Int, (6,6))
        for x in 1:6, y in 1:6
            # ϕ(x)⋅ϕ(y) = ϕ(x⋅y)
            mtab3[ϕ2[x], ϕ2[y]] = ϕ2[mtab1[x,y]]
        end
        @test !isnothing(group_isomorphism(group2, FiniteGroup(mtab3)))  # ϕ and ϕ2 are equivalent

        for ϕ in permutations(2:6)
            ϕ = vcat([1], ϕ)
            mtab1 = group_multiplication_table(group)
            mtab2 = zeros(Int, (6,6))
            for x in 1:6, y in 1:6
                # ϕ(x)⋅ϕ(y) = ϕ(x⋅y)
                mtab2[ϕ[x], ϕ[y]] = ϕ[mtab1[x,y]]
            end
            group2 = FiniteGroup(mtab2)
            @test !isnothing(group_isomorphism(group, group2))
        end
    end

    @testset "minimal_generating_set" begin
        for i in 1:32
            psym = PointSymmetryDatabase.get(i)
            group = psym.group
            mgs = minimal_generating_set(group)
            @test length(mgs) <= length(psym.generators)
            @test generate_subgroup(group, mgs) == BitSet(1:group_order(group))
        end
    end

    @testset "group isomorphism" begin
        # 4
        group1 = FiniteGroup([
            1 2 3 4;
            2 1 4 3;
            3 4 2 1;
            4 3 1 2;
        ])
        group1p= FiniteGroup([
            1 2 3 4;
            2 3 4 1;
            3 4 1 2;
            4 1 2 3;
        ])
        @test !isnothing(group_isomorphism(group1, group1p))
        # 2/m
        group2 = FiniteGroup([
            1 2 3 4;
            2 1 4 3;
            3 4 1 2;
            4 3 2 1;
        ])
        @test isnothing(group_isomorphism(group1, group2))
    end # @testset "group isomorphism" begin

    @testset "group isomorphism 2" begin
        # D2h
        group1 = FiniteGroup([
            1 2 3 4 5 6 7 8;
            2 1 4 3 6 5 8 7;
            3 4 1 2 7 8 5 6;
            4 3 2 1 8 7 6 5;
            5 6 7 8 1 2 3 4;
            6 5 8 7 2 1 4 3;
            7 8 5 6 3 4 1 2;
            8 7 6 5 4 3 2 1;
        ])
        group2 = FiniteGroup([
            1 2 3 4 5 6 7 8;
            2 1 4 3 6 5 8 7;
            3 4 2 1 7 8 6 5;
            4 3 1 2 8 7 5 6;
            5 6 7 8 1 2 3 4;
            6 5 8 7 2 1 4 3;
            7 8 6 5 3 4 2 1;
            8 7 5 6 4 3 1 2;
        ])
        @test isnothing(group_isomorphism(group1, group2))
    end

    @testset "group isomorphism 3" begin
        tsym1 = FiniteTranslationSymmetry([3 0; 0 3])
        tsym2 = FiniteTranslationSymmetry([3 1; 0 3])
        @test isnothing(group_isomorphism(tsym1.group, tsym2.group))
    end

    # TODO
    # Are there two non-isomorphic groups with
    # - Same conjugacy classes
    # - Same period lengths
end
=#


text:
"""
	gmtconvert(cmd0::String="", arg1=nothing, kwargs...)

Convert, Paste, and/or Extract columns from data tables

Full option list at [`gmtconvert`]($(GMTdoc)gmtconvert.html)

Parameters
----------

- **A** | **h_cat** :: [Type => Str | []]

    The records from the input files should be pasted horizontally, not appended vertically [Default].
    ($(GMTdoc)gmtconvert.html#a)
- **C** | **n_records** :: [Type => Str]  ``Arg = [+lmin][+umax][+i]``

    Only output segments whose number of records matches your given criteria:
    ($(GMTdoc)gmtconvert.html#c)
- **D** | **dump** :: [Type => Str | []]   ``Arg = [template[+oorig]]``

    For multiple segment data, dump each segment to a separate output file.
    ($(GMTdoc)gmtconvert.html#d)
- **E** | **first_last** :: [Type => Str | []]   ``Arg = [f|l|m|Mstride]``

    Only extract the first and last record for each segment of interest.
    ($(GMTdoc)gmtconvert.html#e)
- **F** | **conn_method** :: [Type => Str | []]   ``Arg = [c|n|r|v][refpoint]``

    Alter the way points are connected (by specifying a scheme) and data are grouped (by specifying a method).
    ($(GMTdoc)gmtconvert.html#f)
- **I** | **invert** | **reverse** :: [Type => Str | Bool]      ``Arg = [tsr]``

    Invert the order of items, i.e., output the items in reverse order, starting with the last
    and ending up with the first item.
    ($(GMTdoc)gmtconvert.html#i)
- **L** | **list_only** :: -[Type => Bool]

    Only output a listing of all segment header records and no data records.
    ($(GMTdoc)gmtconvert.html#l)
- **N** | **sort** :: [Type => Str | Number]      ``Arg = [-|+]col``

    Numerically sort each segment based on values in column col.
    ($(GMTdoc)gmtconvert.html#n)
- **Q** | **select_num** :: [Type => Str]      ``Arg =  [~]selection``

    Only write segments whose number is included in ``selection`` and skip all others.
    ($(GMTdoc)gmtconvert.html#q)
- **S** | **select_hdr** :: [Type => Str]      ``Arg =  [~]”search string” or [~]/regexp/[i]``

    Only output those segments whose header record contains the specified text string.
    ($(GMTdoc)gmtconvert.html#s)
- **T** | **suppress ** :: [Type => Str | []]    ``Arg = [h|d]``

    Suppress the writing of certain records on output. Append h to suppress segment headers
    [Default] or d to suppress duplicate data records. Use T=:hd to suppress both types of records.
    ($(GMTdoc)gmtconvert.html#t)
- **W** | **word2num** :: [Type => Str | []]      ``Arg = [+n]``

    Attempt to gmtconvert each word in the trialing text to a number and append such values
    to the numerical output columns.
    ($(GMTdoc)gmtconvert.html#w)
- **Z** | **range** :: [Type => Str | []]      ``Arg =  [first][:last]``

    Limit output to the specified record range. If first is not set it defaults to record 0
    (very first record) and if last is not set then it defaults to the very last record.
    ($(GMTdoc)gmtconvert.html#z)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_o)
- $(GMT.opt_s)
- $(GMT.opt_swap_xy)
"""
function gmtconvert(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("gmtconvert", cmd0, arg1)

	d = KW(kwargs)

	cmd = parse_common_opts(d, "", [:V_params :b :d :e :f :g :h :i :o :s :yx])
	cmd = parse_these_opts(cmd, d, [[:A :h_cat], [:C :n_records], [:D :dump], [:E :first_last], [:F :conn_method],
		[:I :invert :reverse], [:L :list_only], [:N :sort], [:Q :select_num], [:S :select_hdr],
		[:T :suppress], [:W :word2num], [:Z :range]])

	common_grd(d, cmd0, cmd, "gmtconvert ", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
gmtconvert(arg1; kw...) = gmtconvert("", arg1; kw...)


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

using Markdown
using InteractiveUtils

# ╔═╡ f29303e2-1d66-11eb-3313-f335871e9674
using Pkg, DrWatson

# ╔═╡ f2935b8a-1d66-11eb-02e8-3dacec840928
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using LinearAlgebra, StatsBase, DataFrames
end

# ╔═╡ 7877839e-1d66-11eb-0a3c-5953d222ccbc
md"## Listing 1.8"

# ╔═╡ f293d704-1d66-11eb-36d8-4bb2f8dff940
md"##### Transition probability matrix."

# ╔═╡ f2a01cf8-1d66-11eb-3573-d91ed64fb7bc
P = [0.5 0.4 0.1;
     0.3 0.2 0.5;
     0.5 0.3 0.2];

# ╔═╡ f2a0c444-1d66-11eb-0f37-99e69755fe8f
md"##### First way."

# ╔═╡ f2a67846-1d66-11eb-022c-6f0a51c33a42
piProb1 = (P^100)[1,:]

# ╔═╡ f2b43238-1d66-11eb-3713-c5c8f2d4af76
md"##### Second way."

# ╔═╡ f2b4d63e-1d66-11eb-3d9f-65ba9b1ac506
begin
	A = vcat((P' - I)[1:2,:],ones(3)')
	b = [0 0 1]'
	piProb2 = (A\b)[:, 1]
end

# ╔═╡ f2c17bc8-1d66-11eb-0569-9b84c6b33cca
md"##### Third way."

# ╔═╡ f2c2e8fa-1d66-11eb-3421-3d76cd8c9a73
begin
	eigVecs = eigvecs(copy(P'))
	highestVec = eigVecs[:,findmax(abs.(eigvals(P)))[2]]
	piProb3 = Array{Float64}(highestVec)/norm(highestVec,1)
end

# ╔═╡ f2d180ba-1d66-11eb-216b-e3e640811323
md"##### Fourth way."

# ╔═╡ f2da9662-1d66-11eb-1f95-3da69865d370
begin
numInState = zeros(Int,3)
	state = 1
	N = 10^6
	for t in 1:N
		numInState[state] += 1
		global state = sample(1:3,weights(P[state,:]))
	end
	piProb4 = numInState/N
end

# ╔═╡ f2dbd9d2-1d66-11eb-12b9-af2d7f1f73da
DataFrame(:piProb1 => piProb1, :piProb2 => piProb2, :piProb3 => piProb3, :piProb4 => piProb4)

# ╔═╡ f2e37750-1d66-11eb-3f4e-e7b249fb2052
md"## End of listing 1.8"

# ╔═╡ Cell order:
# ╟─7877839e-1d66-11eb-0a3c-5953d222ccbc
# ╠═f29303e2-1d66-11eb-3313-f335871e9674
# ╠═f2935b8a-1d66-11eb-02e8-3dacec840928
# ╟─f293d704-1d66-11eb-36d8-4bb2f8dff940
# ╠═f2a01cf8-1d66-11eb-3573-d91ed64fb7bc
# ╟─f2a0c444-1d66-11eb-0f37-99e69755fe8f
# ╠═f2a67846-1d66-11eb-022c-6f0a51c33a42
# ╟─f2b43238-1d66-11eb-3713-c5c8f2d4af76
# ╠═f2b4d63e-1d66-11eb-3d9f-65ba9b1ac506
# ╟─f2c17bc8-1d66-11eb-0569-9b84c6b33cca
# ╠═f2c2e8fa-1d66-11eb-3421-3d76cd8c9a73
# ╟─f2d180ba-1d66-11eb-216b-e3e640811323
# ╠═f2da9662-1d66-11eb-1f95-3da69865d370
# ╠═f2dbd9d2-1d66-11eb-12b9-af2d7f1f73da
# ╟─f2e37750-1d66-11eb-3f4e-e7b249fb2052


text:
using SCS
const SCS_PRIMAL_FACTORY =
    MOI.OptimizerWithAttributes(SCS.Optimizer, MOI.Silent() => true)
const SCS_DUAL_FACTORY = dual_optimizer(SCS_PRIMAL_FACTORY)
const SCS_PRIMAL_OPT = MOI.instantiate(SCS_PRIMAL_FACTORY)
const SCS_DUAL_OPT = MOI.instantiate(SCS_DUAL_FACTORY)

push!(primal_power_cone_factory, SCS_PRIMAL_FACTORY)
push!(dual_power_cone_factory, SCS_DUAL_FACTORY)
push!(dual_power_cone_optimizer, SCS_DUAL_OPT)
push!(primal_power_cone_optimizer, SCS_PRIMAL_OPT)

@testset "SCS Exponential Cone Problems" begin
    list_of_exp_problems = [exp1_test, exp2_test]
    test_strong_duality(list_of_exp_problems, SCS_PRIMAL_FACTORY)
end

@testset "SCS Power Cone Problems" begin
    list_of_pow_problems = [pow1_test, pow2_test]
    test_strong_duality(
        list_of_pow_problems,
        SCS_PRIMAL_FACTORY;
        atol = 1e-3,
        rtol = 1e-3,
    )
end


text:
using Test
using PotentialFitting
using PotentialCalculation
using PotentialDB

using ScikitLearn

@sk_import linear_model: LinearRegression

data=load_data_file(normpath(joinpath(dirname(pathof(PotentialFitting)),"../test", "data", "test.jld")))

m1=MoleculeIdenticalInformation{AtomOnlySymbol}(data["cluster1"].atoms)
m2=MoleculeIdenticalInformation{AtomOnlySymbol}(data["cluster2"].atoms)

topo=[]
push!(topo, PairPotentialTopology(LennardJones(),1,6))
push!(topo, PairPotentialTopology(LennardJones(), 5,6))
push!(topo, PairPotentialTopology(GeneralPowers(),4,6))
push!(topo, PairPotentialTopology(LennardJones(),2,6))
push!(topo, ClusterPotentialTopology(GeneralAngle([-6, -12], [0, 1]), 1,3,6))

mpp1=MoleculePairPotential(m1,m2, LennardJones())
mpp=MoleculePairPotential(m1,m2)

mpp.topology = topo

fdata = FitData(mpp, data["Points"], data["Energy"])

model = LinearRegression()

setweight_e_more!(fdata, 0, 1500)
setweight_e_less!(fdata,4,80)
setweight_e_less!(fdata,16,120)

fit_potential!(model, mpp, fdata)
rmsd(data["Points"],data["Energy"],mpp)

pl = plot_compare(data["Points"][:,1], data["Energy"][:,1], mpp);

@testset "Fit" begin
    r = defaultregistry()
    data = loadpotential(r, "7")

    m1=MoleculeIdenticalInformation{AtomOnlySymbol}(data["cluster1"].atoms)
    m2=MoleculeIdenticalInformation{AtomOnlySymbol}(data["cluster2"].atoms)
    mpp = MoleculePairPotential(m1,m2, GeneralPowers(-6,-8,-10,-12,-7,-9,-11))

    fdata = FitData(mpp, data, data)
    setweight_e_more!(fdata, 0.1, 1000)
    setweight_e_more!(fdata, 0.01, 5000)
    setweight_e_less!(fdata, 2, -50)
    setweight_e_less!(fdata, 4, -75)
    setweight_e_less!(fdata, 16, -100)

    model = LinearRegression(fit_intercept=false)
    fit_potential!(model, mpp, fdata)

    @test rmsd(data["Points"], data["Energy"], mpp, emax=-10) < 10
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""TopologySpreadConstraint specifies how to spread matching pods among the given topology.

    IoK8sApiCoreV1TopologySpreadConstraint(;
        labelSelector=nothing,
        maxSkew=nothing,
        topologyKey=nothing,
        whenUnsatisfiable=nothing,
    )

    - labelSelector::IoK8sApimachineryPkgApisMetaV1LabelSelector : LabelSelector is used to find matching pods. Pods that match this label selector are counted to determine the number of pods in their corresponding topology domain.
    - maxSkew::Int32 : MaxSkew describes the degree to which pods may be unevenly distributed. It's the maximum permitted difference between the number of matching pods in any two topology domains of a given topology type. For example, in a 3-zone cluster, MaxSkew is set to 1, and pods with the same labelSelector spread as 1/1/0: | zone1 | zone2 | zone3 | |   P   |   P   |       | - if MaxSkew is 1, incoming pod can only be scheduled to zone3 to become 1/1/1; scheduling it onto zone1(zone2) would make the ActualSkew(2-0) on zone1(zone2) violate MaxSkew(1). - if MaxSkew is 2, incoming pod can be scheduled onto any zone. It's a required field. Default value is 1 and 0 is not allowed.
    - topologyKey::String : TopologyKey is the key of node labels. Nodes that have a label with this key and identical values are considered to be in the same topology. We consider each <key, value> as a \"bucket\", and try to put balanced number of pods into each bucket. It's a required field.
    - whenUnsatisfiable::String : WhenUnsatisfiable indicates how to deal with a pod if it doesn't satisfy the spread constraint. - DoNotSchedule (default) tells the scheduler not to schedule it - ScheduleAnyway tells the scheduler to still schedule it It's considered as \"Unsatisfiable\" if and only if placing incoming pod on any topology violates \"MaxSkew\". For example, in a 3-zone cluster, MaxSkew is set to 1, and pods with the same labelSelector spread as 3/1/1: | zone1 | zone2 | zone3 | | P P P |   P   |   P   | If WhenUnsatisfiable is set to DoNotSchedule, incoming pod can only be scheduled to zone2(zone3) to become 3/2/1(3/1/2) as ActualSkew(2-1) on zone2(zone3) satisfies MaxSkew(1). In other words, the cluster can still be imbalanced, but scheduler won't make it *more* imbalanced. It's a required field.
"""
mutable struct IoK8sApiCoreV1TopologySpreadConstraint <: SwaggerModel
    labelSelector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: labelSelector
    maxSkew::Any # spec type: Union{ Nothing, Int32 } # spec name: maxSkew
    topologyKey::Any # spec type: Union{ Nothing, String } # spec name: topologyKey
    whenUnsatisfiable::Any # spec type: Union{ Nothing, String } # spec name: whenUnsatisfiable

    function IoK8sApiCoreV1TopologySpreadConstraint(;labelSelector=nothing, maxSkew=nothing, topologyKey=nothing, whenUnsatisfiable=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1TopologySpreadConstraint, Symbol("labelSelector"), labelSelector)
        setfield!(o, Symbol("labelSelector"), labelSelector)
        validate_property(IoK8sApiCoreV1TopologySpreadConstraint, Symbol("maxSkew"), maxSkew)
        setfield!(o, Symbol("maxSkew"), maxSkew)
        validate_property(IoK8sApiCoreV1TopologySpreadConstraint, Symbol("topologyKey"), topologyKey)
        setfield!(o, Symbol("topologyKey"), topologyKey)
        validate_property(IoK8sApiCoreV1TopologySpreadConstraint, Symbol("whenUnsatisfiable"), whenUnsatisfiable)
        setfield!(o, Symbol("whenUnsatisfiable"), whenUnsatisfiable)
        o
    end
end # type IoK8sApiCoreV1TopologySpreadConstraint

const _property_map_IoK8sApiCoreV1TopologySpreadConstraint = Dict{Symbol,Symbol}(Symbol("labelSelector")=>Symbol("labelSelector"), Symbol("maxSkew")=>Symbol("maxSkew"), Symbol("topologyKey")=>Symbol("topologyKey"), Symbol("whenUnsatisfiable")=>Symbol("whenUnsatisfiable"))
const _property_types_IoK8sApiCoreV1TopologySpreadConstraint = Dict{Symbol,String}(Symbol("labelSelector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("maxSkew")=>"Int32", Symbol("topologyKey")=>"String", Symbol("whenUnsatisfiable")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1TopologySpreadConstraint }) = collect(keys(_property_map_IoK8sApiCoreV1TopologySpreadConstraint))
Swagger.property_type(::Type{ IoK8sApiCoreV1TopologySpreadConstraint }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1TopologySpreadConstraint[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1TopologySpreadConstraint }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1TopologySpreadConstraint[property_name]

function check_required(o::IoK8sApiCoreV1TopologySpreadConstraint)
    (getproperty(o, Symbol("maxSkew")) === nothing) && (return false)
    (getproperty(o, Symbol("topologyKey")) === nothing) && (return false)
    (getproperty(o, Symbol("whenUnsatisfiable")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1TopologySpreadConstraint }, name::Symbol, val)
end