text:
export @const_gate,
    ConstantGate,
    ConstGate,
    X,
    Y,
    Z,
    H,
    I2,
    SWAP,
    T,
    XGate,
    YGate,
    ZGate,
    HGate,
    I2Gate,
    SWAPGate,
    TGate

"""
    ConstGate

A module contains all constant gate definitions.
"""
module ConstGate
import ..YaoBlocks
export ConstantGate, PauliGate

"""
    ConstantGate{N, T} <: PrimitiveBlock{N, T}

Abstract type for constant gates. Constant gates are
quantum gates with constant value, e.g Pauli gates,
T gate, Hadmard gates, etc.
"""
abstract type ConstantGate{N} <: YaoBlocks.PrimitiveBlock{N} end

include("const_gate_tools.jl")
include("const_gate_gen.jl")

const PauliGate = Union{I2Gate,XGate,YGate,ZGate}

YaoBlocks.cache_key(x::ConstantGate) = 0x1

struct IGate{N} <: ConstantGate{N} end
YaoBlocks.mat(::Type{T}, ::IGate{N}) where {T,N} = IMatrix{1 << N,T}()

YaoBase.ishermitian(::IGate) = true
YaoBase.isunitary(::IGate) = true
end

# import some frequently-used objects
import .ConstGate:
    XGate,
    YGate,
    ZGate,
    HGate,
    I2Gate,
    SWAPGate,
    TGate,
    X,
    Y,
    Z,
    H,
    I2,
    SWAP,
    T,
    @const_gate,
    ConstantGate,
    PauliGate
occupied_locs(::I2Gate) = ()


text:
struct OffsetPrecalc{T,N,C,B,R,X,M,P<:AbstractStridedPointer{T,N,C,B,R,X,M},I} <: AbstractStridedPointer{T,N,C,B,R,X,M}
    ptr::P
    precalc::I
end
@inline Base.pointer(ptr::OffsetPrecalc) = pointer(getfield(ptr, :ptr))
@inline Base.similar(ptr::OffsetPrecalc, p::Ptr) = OffsetPrecalc(similar(getfield(ptr, :ptr), p), getfield(ptr, :precalc))
# @inline pointerforcomparison(p::OffsetPrecalc) = pointerforcomparison(getfield(p, :ptr))
# @inline pointerforcomparison(p::OffsetPrecalc, i::Tuple) = pointerforcomparison(p.ptr, i)
@inline offsetprecalc(x, ::Any) = x
@inline offsetprecalc(x::StridedBitPointer, ::Val) = x
# @inline pointerforcomparison(p::AbstractStridedPointer) = pointer(p)
# @inline pointerforcomparison(p::AbstractStridedPointer, i) = gep(p, i)
@inline ArrayInterface.offsets(p::OffsetPrecalc) = offsets(getfield(p, :ptr))

@inline Base.strides(p::OffsetPrecalc) = strides(getfield(p, :ptr))
@inline ArrayInterface.strides(p::OffsetPrecalc) = strides(getfield(p, :ptr))

@inline function similar_no_offset(sptr::OffsetPrecalc{T}, ptr::Ptr{T}) where {T}
    OffsetPrecalc(similar_no_offset(getfield(sptr, :ptr), ptr), getfield(sptr, :precalc))
end

"""

Basically:

if I ∈ [3,5,7,9]
    c[(I - 1) >> 1]
else
    b * I
end

because

c = b .* [3, 5, 7, 9]

"""
@generated function lazymul(::StaticInt{I}, b, c::Tuple{Vararg{Any,N}}) where {I,N}
    Is = (I - 1) >> 1
    ex = if (isodd(I) && 1 ≤ Is ≤ N) && (c.parameters[Is] !== nothing)
        Expr(:call, GlobalRef(Core, :getfield), :c, Is, false)
    elseif (I ∈ (6, 10) && (I >> 2 ≤ N)) && (c.parameters[I >> 2] !== nothing)
        Expr(:call, :lazymul, Expr(:call, Expr(:curly, :StaticInt, 2)), Expr(:call, GlobalRef(Core, :getfield), :c, I >> 2, false))
    else
        Expr(:call, :lazymul, Expr(:call, Expr(:curly, :StaticInt, I)), :b)
    end
    Expr(:block, Expr(:meta,:inline), ex)
end
@inline lazymul(a, b, c) = lazymul(a, b)
@inline lazymul(a::StaticInt, b, ::Nothing) = lazymul(a, b)

_unwrap(::Type{StaticInt{N}}) where {N} = N
_unwrap(_) = nothing
# descript is a tuple of (unrollfactor) for each ind; if it shouldn't preallocate, unrollfactor may be set to 1
function precalc_quote_from_descript(descript, contig, X)
    precalc = Expr(:tuple)
    anyprecalcs = anydynamicprecals = false
    pstrideextracts = Expr(:block)
    for (i,uf) ∈ enumerate(descript)
        if i == contig || uf < 3
            push!(precalc.args, nothing)
        else
            t = Expr(:tuple)
            Xᵢ = _unwrap(X[i])
            anyprecalcs = true
            if Xᵢ === nothing
                anydynamicprecals = true
                pstride_i = Symbol(:pstride_, i)
                push!(pstrideextracts.args, Expr(:(=), pstride_i, Expr(:call, GlobalRef(Core, :getfield), :pstride, i, false)))
                foreach(u -> push!(t.args, Expr(:call, :vmul_fast, u, pstride_i)), 3:2:uf)
            else
                foreach(u -> push!(t.args, u * Xᵢ), 3:2:uf)
            end
            push!(precalc.args, t)
        end
    end
    q = Expr(:block, Expr(:meta,:inline))
    if anydynamicprecals
        push!(q.args, :(pstride = strides(p)))
        push!(q.args, pstrideextracts)
    end
    if anyprecalcs
        push!(q.args, Expr(:call, :OffsetPrecalc, :p, precalc))
    else
        push!(q.args, :p)
    end
    q
end
@generated function offsetprecalc(p::AbstractStridedPointer{T,N,C,B,R,X,O}, ::Val{descript}) where {T,N,C,B,R,X,O,descript}
    precalc_quote_from_descript(descript, C, X.parameters)
end

@inline tdot(ptr::OffsetPrecalc{T}, a, b) where {T} = tdot(pointer(ptr), a, b, getfield(ptr, :precalc))
@inline tdot(ptr::OffsetPrecalc, ::Tuple{}, ::Tuple{}) = (pointer(ptr), Zero())

ext: jl


text:
pts = invariantize([rand(3) for i = 1:15])
triang = triangulate(pts)

to_approx =  transferoperator(pts, TriangulationBinning(), ExactIntersection())
to_exact = transferoperator(pts, TriangulationBinning(), ApproximateIntersection())

@test to_approx isa TransferOperatorTriangulationApprox
@test to_exact isa TransferOperatorTriangulationExact

@test all(sum(to_approx.transfermatrix, dims = 2) .≈ 1.0)
@test all(sum(to_exact.transfermatrix, dims = 2) .≈ 1.0)


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

mutable struct IoK8sApiCoreV1PodSecurityContext <: SwaggerModel
    fsGroup::Any # spec type: Union{ Nothing, Int64 } # spec name: fsGroup
    runAsGroup::Any # spec type: Union{ Nothing, Int64 } # spec name: runAsGroup
    runAsNonRoot::Any # spec type: Union{ Nothing, Bool } # spec name: runAsNonRoot
    runAsUser::Any # spec type: Union{ Nothing, Int64 } # spec name: runAsUser
    seLinuxOptions::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SELinuxOptions } # spec name: seLinuxOptions
    supplementalGroups::Any # spec type: Union{ Nothing, Vector{Int64} } # spec name: supplementalGroups
    sysctls::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1Sysctl} } # spec name: sysctls

    function IoK8sApiCoreV1PodSecurityContext(;fsGroup=nothing, runAsGroup=nothing, runAsNonRoot=nothing, runAsUser=nothing, seLinuxOptions=nothing, supplementalGroups=nothing, sysctls=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PodSecurityContext, Symbol("fsGroup"), fsGroup)
        setfield!(o, Symbol("fsGroup"), fsGroup)
        validate_property(IoK8sApiCoreV1PodSecurityContext, Symbol("runAsGroup"), runAsGroup)
        setfield!(o, Symbol("runAsGroup"), runAsGroup)
        validate_property(IoK8sApiCoreV1PodSecurityContext, Symbol("runAsNonRoot"), runAsNonRoot)
        setfield!(o, Symbol("runAsNonRoot"), runAsNonRoot)
        validate_property(IoK8sApiCoreV1PodSecurityContext, Symbol("runAsUser"), runAsUser)
        setfield!(o, Symbol("runAsUser"), runAsUser)
        validate_property(IoK8sApiCoreV1PodSecurityContext, Symbol("seLinuxOptions"), seLinuxOptions)
        setfield!(o, Symbol("seLinuxOptions"), seLinuxOptions)
        validate_property(IoK8sApiCoreV1PodSecurityContext, Symbol("supplementalGroups"), supplementalGroups)
        setfield!(o, Symbol("supplementalGroups"), supplementalGroups)
        validate_property(IoK8sApiCoreV1PodSecurityContext, Symbol("sysctls"), sysctls)
        setfield!(o, Symbol("sysctls"), sysctls)
        o
    end
end # type IoK8sApiCoreV1PodSecurityContext

const _property_map_IoK8sApiCoreV1PodSecurityContext = Dict{Symbol,Symbol}(Symbol("fsGroup")=>Symbol("fsGroup"), Symbol("runAsGroup")=>Symbol("runAsGroup"), Symbol("runAsNonRoot")=>Symbol("runAsNonRoot"), Symbol("runAsUser")=>Symbol("runAsUser"), Symbol("seLinuxOptions")=>Symbol("seLinuxOptions"), Symbol("supplementalGroups")=>Symbol("supplementalGroups"), Symbol("sysctls")=>Symbol("sysctls"))
const _property_types_IoK8sApiCoreV1PodSecurityContext = Dict{Symbol,String}(Symbol("fsGroup")=>"Int64", Symbol("runAsGroup")=>"Int64", Symbol("runAsNonRoot")=>"Bool", Symbol("runAsUser")=>"Int64", Symbol("seLinuxOptions")=>"IoK8sApiCoreV1SELinuxOptions", Symbol("supplementalGroups")=>"Vector{Int64}", Symbol("sysctls")=>"Vector{IoK8sApiCoreV1Sysctl}")
Base.propertynames(::Type{ IoK8sApiCoreV1PodSecurityContext }) = collect(keys(_property_map_IoK8sApiCoreV1PodSecurityContext))
Swagger.property_type(::Type{ IoK8sApiCoreV1PodSecurityContext }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1PodSecurityContext[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PodSecurityContext }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PodSecurityContext[property_name]

function check_required(o::IoK8sApiCoreV1PodSecurityContext)
    true
end

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


text:
module CartesianGP

include("BitString.jl")
include("Func.jl")
include("Parameters.jl")
include("Node.jl")
include("Chromosome.jl")
include("Goal.jl")
include("Fitness.jl")
include("Mutate.jl")
include("Execute.jl")
include("Evolution.jl")
include("Utilities.jl")

end


text:
""" Graphviz drawing of categories, functors, diagrams, and related structures.
"""
module GraphvizCategories
export to_graphviz, to_graphviz_property_graph

using ...Syntax, ...Present, ...Theories
using ...CategoricalAlgebra, ...Graphs, ..GraphvizGraphs
import ..Graphviz
import ..GraphvizGraphs: to_graphviz, to_graphviz_property_graph

# Presentations
###############

to_graphviz(pres::Presentation, args...; kw...) =
  to_graphviz(to_graphviz_property_graph(pres, args...; kw...))

function to_graphviz_property_graph(pres::Presentation, Ob::Symbol, Hom::Symbol;
    prog::AbstractString="neato", graph_attrs::AbstractDict=Dict(),
    node_attrs::AbstractDict=Dict(), edge_attrs::AbstractDict=Dict())
  pg = PropertyGraph{Any}(; prog = prog,
    graph = graph_attrs,
    node = merge!(Dict(:shape => "ellipse", :margin => "0"), node_attrs),
    edge = edge_attrs,
  )

  obs = generators(pres, Ob)
  add_vertices!(pg, length(obs))
  for (v, ob) in enumerate(obs)
    set_vprop!(pg, v, :label, string(first(ob)))
  end

  homs = generators(pres, Hom)
  add_edges!(pg, map(f -> generator_index(pres, first(gat_type_args(f))), homs),
                 map(f -> generator_index(pres, last(gat_type_args(f))), homs))
  for (e, hom) in enumerate(homs)
    set_eprop!(pg, e, :label, string(first(hom)))
  end
  pg
end

# Categories
############

function to_graphviz_property_graph(pres::Presentation{Category}; kw...)
  to_graphviz_property_graph(pres, :Ob, :Hom; kw...)
end

function to_graphviz_property_graph(pres::Presentation{MCategory};
    tight_attrs::AbstractDict=Dict(:arrowhead => "empty"), kw...)
  pg = to_graphviz_property_graph(pres, :Ob, :Hom; kw...)
  for tight_hom in generators(pres, :Tight)
    e = generator_index(pres, hom(tight_hom))
    set_eprops!(pg, e, tight_attrs)
  end
  pg
end

to_graphviz(X::AbstractElements; kw...) =
  to_graphviz(to_graphviz_property_graph(X; kw...))

function to_graphviz_property_graph(X::AbstractElements;
    prog::AbstractString="neato", graph_attrs::AbstractDict=Dict(),
    node_attrs::AbstractDict=Dict(), edge_attrs::AbstractDict=Dict(),
    node_labels::Bool=false, edge_labels::Bool=false)
  pg = PropertyGraph{Any}(; prog = prog,
    graph = graph_attrs,
    node = merge!(Dict(:shape => "box", :width => "0", :height => "0",
                       :margin => "0.025"), node_attrs),
    edge = edge_attrs)

  add_vertices!(pg, nparts(X, :El))
  add_edges!(pg, X[:src], X[:tgt])

  for v in parts(X, :El)
    vlabel = X[X[v, :πₑ], :nameo]
    set_vprop!(pg, v, :label, node_labels ? "$v:$vlabel" : "$vlabel")
  end
  for e in parts(X, :Arr)
    elabel = X[X[e, :πₐ], :nameh]
    set_eprop!(pg, e, :label, edge_labels ? "$e:$elabel" : "$elabel")
  end
  pg
end

# Diagrams
##########

to_graphviz(F::FinFunctor; kw...) =
  to_graphviz(to_graphviz_property_graph(F; kw...))

function to_graphviz_property_graph(F::FinFunctor; kw...)
  g = graph(dom(F))
  pg = to_graphviz_property_graph(g; kw...)
  for v in vertices(g)
    lᵥ = g[v, :vname]
    tᵥ = ob_map(F, v)
    set_vprop!(pg, v, :label, string(lᵥ, ":", tᵥ))
  end
  for e in edges(g)
    tₑ = hom_map(F, e)
    set_eprop!(pg, e, :label, string(tₑ))
  end
  pg
end

# Schemas
#########

function to_graphviz_property_graph(pres::Presentation{Schema}; kw...)
  pg = to_graphviz_property_graph(pres, :Ob, :Hom; kw...)
  ob_vertices = vertices(pg)
  hom_edges = edges(pg)

  attr_types = generators(pres, :AttrType)
  attr_vertices = add_vertices!(pg, length(attr_types))
  for (v, attr_type) in zip(attr_vertices, attr_types)
    set_vprops!(pg, v, xlabel=string(first(attr_type)), shape="point")
  end

  attrs = generators(pres, :Attr)
  attr_edges = add_edges!(pg,
    map(attr -> ob_vertices[generator_index(pres, dom(attr))], attrs),
    map(attr -> attr_vertices[generator_index(pres, codom(attr))], attrs))
  for (e, attr) in zip(attr_edges, attrs)
    set_eprop!(pg, e, :label, string(first(attr)))
  end
  pg
end

end

max_forks_repo_name: epatters/CompCat


text:
##############################################################################
# From client
@dict_readable struct Capabilities
    dynamicRegistration::Union{Bool,Missing}
end

const ResourceOperationKind = String
const ResourceOperationKinds = ("create", "rename", "delete")
const FailureHandlingKind = String
const FailureHandlingKinds = ("abort", "transactional", "undo", "textOnlyTransactional")

@dict_readable struct WorkspaceEditCapabilities
    documentChanges::Union{Bool,Missing}
    resourceOperations::Union{Vector{ResourceOperationKind},Missing}
    failureHandling::Union{FailureHandlingKind,Missing}
end

@dict_readable struct SymbolKindCapabilities
    valueSet::Union{Vector{Int},Missing}
end

@dict_readable mutable struct SymbolCapabilities
    dynamicRegistration::Union{Bool,Missing}
    symbolKind::Union{SymbolKindCapabilities,Missing}
end

@dict_readable struct WorkspaceClientCapabilities
    applyEdit::Union{Bool,Missing}
    workspaceEdit::Union{WorkspaceEditCapabilities,Missing}
    didChangeConfiguration::Union{Capabilities,Missing}
    didChangeWatchedFiles::Union{Capabilities,Missing}
    symbol::Union{SymbolCapabilities,Missing}
    executeCommand::Union{Capabilities,Missing}
    workspaceFolders::Union{Bool,Missing}
    configuration::Union{Bool,Missing}
end

@dict_readable struct SynchronizationCapabilities
    dynamicRegistration::Union{Bool,Missing}
    willSave::Union{Bool,Missing}
    willSaveWaitUntil::Union{Bool,Missing}
    didSave::Union{Bool,Missing}
end

@dict_readable struct CompletionItemCapabilities
    snippetSupport::Union{Bool,Missing}
    commitCharactersSupport::Union{Bool,Missing}
    documentationFormat::Union{Vector{String},Missing}
    deprecatedSupport::Union{Bool,Missing}
    preselectSupport::Union{Bool,Missing}
end

@dict_readable struct CompletionItemKindCapabilities
    valueSet::Union{Vector{Int},Missing}
end

@dict_readable struct CompletionCapabilities
    dynamicRegistration::Union{Bool,Missing}
    completionItem::Union{CompletionItemCapabilities,Missing}
    completionItemKind::Union{CompletionItemKindCapabilities,Missing}
    contextSupport::Union{Bool,Missing}
end

@dict_readable struct HoverCapabilities
    dynamicRegistration::Union{Bool,Missing}
    contentFormat::Union{Vector{String},Missing}
end

@dict_readable struct ParameterInformationCapabilities
    labelOffsetSupport::Union{Bool,Missing}
end

@dict_readable struct SignatureInformationCapabilities
    documentationFormat::Union{Vector{String},Missing}
    parameterInformation::Union{ParameterInformationCapabilities,Missing}
end

@dict_readable struct SignatureCapabilities
    dynamicRegistration::Union{Bool,Missing}
    signatureInformation::Union{SignatureInformationCapabilities,Missing}
end

@dict_readable struct DocumentSymbolCapabilities
    dynamicRegistration::Union{Bool,Missing}
    symbolKind::Union{SymbolKindCapabilities,Missing}
    hierarchicalDocumentSymbolSupport::Union{Bool,Missing}
end

@dict_readable struct DeclarationCapabilities
    dynamicRegistration::Union{Bool,Missing}
    linkSupport::Union{Bool,Missing}
end

@dict_readable struct DefinitionCapabilities
    dynamicRegistration::Union{Bool,Missing}
    linkSupport::Union{Bool,Missing}
end

@dict_readable struct TypeDefinitionCapabilities
    dynamicRegistration::Union{Bool,Missing}
    linkSupport::Union{Bool,Missing}
end

@dict_readable struct ImplementationCapabilities
    dynamicRegistration::Union{Bool,Missing}
    linkSupport::Union{Bool,Missing}
end

@dict_readable struct CodeActionKindCapabilities
    valueSet::Vector{String}
end

@dict_readable struct CodeActionLiteralCapabilities
    codeActionKind::CodeActionKindCapabilities
end

@dict_readable struct CodeActionCapabilities
    dynamicRegistration::Union{Bool,Missing}
    codeActionLiteralSupport::Union{CodeActionLiteralCapabilities,Missing}
end

@dict_readable struct RenameCapabilities
    dynamicRegistration::Union{Bool,Missing}
    prepareSupport::Union{Bool,Missing}
end

@dict_readable struct PublishDiagnosticsCapabilities
    relatedInformation::Union{Bool,Missing}
end

@dict_readable struct FoldingRangeCapabilities
    dynamicRegistration::Union{Bool,Missing}
    rangeLimit::Union{Int,Missing}
    lineFoldingOnly::Union{Bool,Missing}
end

@dict_readable struct TextDocumentClientCapabilities
    synchronization::Union{SynchronizationCapabilities,Missing}
    completion::Union{CompletionCapabilities,Missing}
    hover::Union{HoverCapabilities,Missing}
    signatureHelp::Union{SignatureCapabilities,Missing}
    references::Union{Capabilities,Missing}
    documentHighlight::Union{Capabilities,Missing}
    documentSymbol::Union{DocumentSymbolCapabilities,Missing}
    formatting::Union{Capabilities,Missing}
    rangeFormatting::Union{Capabilities,Missing}
    onTypeFormatting::Union{Capabilities,Missing}
    declaration::Union{DeclarationCapabilities,Missing}
    definition::Union{DefinitionCapabilities,Missing}
    typeDefinition::Union{TypeDefinitionCapabilities,Missing}
    implementation::Union{ImplementationCapabilities,Missing}
    codeAction::Union{CodeActionCapabilities,Missing}
    codeLens::Union{Capabilities,Missing}
    documentLink::Union{Capabilities,Missing}
    colorProvider::Union{Capabilities,Missing}
    rename::Union{RenameCapabilities,Missing}
    publishDiagnostics::Union{PublishDiagnosticsCapabilities,Missing}
    foldingRange::Union{FoldingRangeCapabilities,Missing}
end

@dict_readable struct ClientCapabilities
    workspace::Union{WorkspaceClientCapabilities,Missing}
    textDocument::Union{TextDocumentClientCapabilities,Missing}
    window::Union{Any,Missing}
    experimental::Union{Any,Missing}
end

struct InitializeParams
    processId::Union{Int,Nothing}
    rootPath::Union{DocumentUri,Nothing,Missing}
    rootUri::Union{DocumentUri,Nothing}
    initializationOptions::Union{Any,Missing}
    capabilities::ClientCapabilities
    trace::Union{String,Missing}
    workspaceFolders::Union{Vector{WorkspaceFolder},Nothing,Missing}
end

# Requires handwritten implementaiton to account for 3-part Unions
function InitializeParams(dict::Dict)
    InitializeParams(Int(dict["processId"]),
    !haskey(dict, "rootPath") ? missing : dict["rootPath"] === nothing ? nothing : DocumentUri(dict["rootPath"]),
    dict["rootUri"] === nothing ? nothing : DocumentUri(dict["rootUri"]),
    if haskey(dict, "initializationOptions") dict["initializationOptions"] else missing end,
    ClientCapabilities(dict["capabilities"]),
    if haskey(dict, "trace") String(dict["trace"]) else missing end,
    !haskey(dict, "workspaceFolders") ? missing : dict["workspaceFolders"] === nothing ? nothing : WorkspaceFolder.(dict["workspaceFolders"]))
end
##############################################################################

##############################################################################
# Server Response
struct CompletionOptions <: Outbound
    resolveProvider::Union{Bool,Missing}
    triggerCharacters::Union{Vector{String},Missing}
end

struct SignatureHelpOptions <: Outbound
    triggerCharacters::Union{Vector{String},Missing}
end

const CodeActionKind = String
const CodeActionKinds = ("", "quickfix", "refactor", "refactor.extract", "refactor.inline", "source", "source.organiseImports")

struct CodeActionOptions <: Outbound
    codeActionKinds::Union{Vector{CodeActionKind},Missing}
end

struct CodeLensOptions <: Outbound
    resolveProvider::Union{Bool,Missing}
end

struct DocumentOnTypeFormattingOptions <: Outbound
    firstTriggerCharacter::String
    moreTriggerCharacters::Union{Vector{String},Missing}
end

struct RenameOptions <: Outbound
    prepareProvider::Union{Bool,Missing}
end

struct DocumentLinkOptions <: Outbound
    resolveProvider::Union{Bool,Missing}
end

struct ExecuteCommandOptions <: Outbound
    commands::Vector{String}
end

struct SaveOptions <: Outbound
    includeText::Union{Bool,Missing}
end

struct ColorProviderOptions <: Outbound end

struct FoldingRangeProviderOptions <: Outbound end

const TextDocumentSyncKind = Int
const TextDocumentSyncKinds = Dict("None" => 0, "Full" => 1, "Incremental" => 2)

struct TextDocumentSyncOptions <: Outbound
    openClose::Union{Bool,Missing}
    change::Union{TextDocumentSyncKind,Missing}
    willSave::Union{Bool,Missing}
    willSaveWaitUntil::Union{Bool,Missing}
    save::Union{SaveOptions,Missing}
end

struct WorkspaceFoldersOptions <: Outbound
    supported::Union{Bool,Missing}
    changeNotifications::Union{Bool,String,Missing}
end

struct WorkspaceOptions <: Outbound
    workspaceFolders::Union{WorkspaceFoldersOptions,Missing}
end

struct ServerCapabilities <: Outbound
    textDocumentSync::Union{TextDocumentSyncOptions,TextDocumentSyncKind,Missing}
    hoverProvider::Union{Bool,Missing}
    completionProvider::Union{CompletionOptions,Missing}
    signatureHelpProvider::Union{SignatureHelpOptions,Missing}
    definitionProvider::Union{Bool,Missing}
    typeDefinitionProvider::Union{Bool,Missing}
    implementationProvider::Union{Bool,Missing}
    referencesProvider::Union{Bool,Missing}
    documentHighlightProvider::Union{Bool,Missing}
    documentSymbolProvider::Union{Bool,Missing}
    workspaceSymbolProvider::Union{Bool,Missing}
    codeActionProvider::Union{Bool,CodeActionOptions,Missing}
    codeLensProvider::Union{CodeLensOptions,Missing}
    documentFormattingProvider::Union{Bool,Missing}
    documentRangeFormattingProvider::Union{Bool,Missing}
    documentOnTypeFormattingProvider::Union{DocumentOnTypeFormattingOptions,Missing}
    renameProvider::Union{Bool,RenameOptions,Missing}
    documentLinkProvider::Union{DocumentLinkOptions,Missing}
    colorProvider::Union{Bool,ColorProviderOptions,Missing}
    foldingRangeProvider::Union{Bool,FoldingRangeProviderOptions,Missing}
    declarationProvider::Union{Bool,Missing}
    executeCommandProvider::Union{ExecuteCommandOptions,Missing}
    workspace::Union{WorkspaceOptions,Missing}
    experimental::Union{Any,Missing}
end

struct InitializeResult <: Outbound
    capabilities::ServerCapabilities
end

##############################################################################
struct InitializedParams end


text:
#########
#File: c:\Users\digan\Dropbox\Dynamic_Networks\repos\ScoreDrivenExponentialRandomGraphs\_research\analysis_for_paper_revision\new_application_\wiki_load_change_stats_and_estimate.jl
#Project: c:\Users\digan\Dropbox\Dynamic_Networks\repos\ScoreDrivenExponentialRandomGraphs\_research\analysis_for_paper_revision\new_application_
#Created Date: Friday April 23rd 2021
#Author: Domenico Di Gangi,  <digangidomenico@gmail.com>
#-----
#Last Modified: Thursday June 3rd 2021 5:18:21 pm
#Modified By:  Domenico Di Gangi
#-----
#Description:
#-----
########

#%%
using DrWatson
@quickactivate "ScoreDrivenExponentialRandomGraphs"
DrWatson.greet()

using Distributed

begin
using TableView
using Blink
viewtab(df) = body!(Window(), showtable(df))
using Logging
using ScoreDrivenERGM
import ScoreDrivenERGM:StaticNets,DynNets, Utilities, ErgmRcall
using LinearAlgebra
using CSV
using DataFrames
using Dates
using PyPlot
pygui(true)

end
#%%
#load data
@elapsed dfEst = collect_results( datadir("reddit_hyperlinks", "ch_stats_present"))
viewtab(dfEst[Not(:ch_stats)])
begin
row = dfEst[dfEst.f .== "w",:][1,:]
T = length(row.ch_stats)
model = row.model
t0_train = 1
tend_train = T
N = row.n_nodes_T[t0_train:tend_train]
obsT = row.ch_stats[t0_train:tend_train]

ENV["JULIA_DEBUG"] = ScoreDrivenERGM

res_est = DynNets.estimate_and_filter(model, N, obsT; show_trace = true)
fig, ax = DynNets.plot_filtered(model, N, res_est[3])

# res_conf = DynNets.conf_bands_given_SD_estimates(model, N, obsT, DynNets.unrestrict_all_par(model,  res_est.vEstSdResParAll), res_est.ftot_0, [[0.975, 0.025]];  offset=0, plotFlag=true, parUncMethod = "WHITE-MLE", winsorProp = 0.25)
# res_conf.ax[1].set_title("college msgs $(model.staticModel.ergmTermsString)")

end
estSS_T = DynNets.estimate_single_snap_sequence(model, obsT)
fig, ax = DynNets.plot_filtered(model, N, estSS_T, ax=ax, lineType = ".", lineColor="r")

using Statistics
var(res_est.fVecT_filt[:, 100:end], dims = 2)

DynNets.plot_filtered(model, N, estParSS_T; lineType = ".", lineColor = "k", fig = res_conf.fig, ax=res_conf.ax, gridFlag=false)

DynNets.plot_filtered(model, row.n_nodes_T, estSS_T; lineType = ".", lineColor = "k")

DynNets.static_estimate(model, obsT)

vEstSdResPar = res_est.vEstSdResParAll
fVecT_filt , target_fun_val_T, sVecT_filt = DynNets.score_driven_filter(model, N, obsT,  vEstSdResPar, model.indTvPar;ftot_0 = res_est.ftot_0)

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
@testset "Predicates" begin
    # ispalindromic
    @test  ispalindromic(dna"")
    @test !ispalindromic(dna"A")
    @test !ispalindromic(dna"C")
    @test  ispalindromic(dna"AT")
    @test  ispalindromic(dna"CG")
    @test !ispalindromic(dna"AC")
    @test !ispalindromic(dna"TT")
    @test  ispalindromic(dna"ANT")
    @test  ispalindromic(dna"ACGT")
    @test !ispalindromic(dna"ACNT")

    @test  ispalindromic(mer"ACGT"dna)
    @test !ispalindromic(mer"CACG"dna)

    @test  ispalindromic(rna"")
    @test !ispalindromic(rna"A")
    @test !ispalindromic(rna"C")
    @test  ispalindromic(rna"AU")
    @test  ispalindromic(rna"CG")
    @test !ispalindromic(rna"AC")
    @test !ispalindromic(rna"UU")
    @test  ispalindromic(rna"ANU")
    @test  ispalindromic(rna"ACGU")
    @test !ispalindromic(rna"ACNU")

    @test  ispalindromic(mer"ACGU"rna)
    @test !ispalindromic(mer"CACG"rna)

    @test_throws Exception ispalindromic(aa"PQ")

    # hasambiguity
    @test !hasambiguity(dna"")
    @test !hasambiguity(dna"A")
    @test  hasambiguity(dna"N")
    @test !hasambiguity(dna"ACGT")
    @test  hasambiguity(dna"ANGT")

    @test !hasambiguity(mer"ACGT"dna)

    @test !hasambiguity(rna"")
    @test !hasambiguity(rna"A")
    @test  hasambiguity(rna"N")
    @test !hasambiguity(rna"ACGU")
    @test  hasambiguity(rna"ANGU")

    @test !hasambiguity(mer"ACGU"rna)

    @test !hasambiguity(aa"")
    @test !hasambiguity(aa"A")
    @test !hasambiguity(aa"P")
    @test  hasambiguity(aa"B")
    @test  hasambiguity(aa"X")
    @test !hasambiguity(aa"ARNDCQEGHILKMFPSTWYVOU")
    @test  hasambiguity(aa"ARXDCQEGHILKMFPSTWYVOU")

    # isrepetitive
    @test  isrepetitive(dna"")
    @test !isrepetitive(dna"", 1)
    @test  isrepetitive(dna"A")
    @test  isrepetitive(dna"A", 1)
    @test  isrepetitive(dna"AAA")
    @test !isrepetitive(dna"ACGT", 2)
    @test  isrepetitive(dna"AAGT", 2)
    @test  isrepetitive(dna"ACCG", 2)
    @test  isrepetitive(dna"ACGG", 2)
    @test !isrepetitive(dna"ACGTCCGT", 3)
    @test  isrepetitive(dna"ACGCCCGT", 3)

    @test !isrepetitive(mer"ACGT"dna, 2)
    @test  isrepetitive(mer"ACCT"dna, 2)
    @test !isrepetitive(mer"ACCT"dna, 3)
    @test  isrepetitive(mer"ACCCT"dna, 2)
    @test  isrepetitive(mer"ACCCT"dna, 3)

    @test  isrepetitive(rna"")
    @test !isrepetitive(rna"", 1)
    @test  isrepetitive(rna"A")
    @test  isrepetitive(rna"A", 1)
    @test  isrepetitive(rna"AAA")
    @test !isrepetitive(rna"ACGU", 2)
    @test  isrepetitive(rna"AAGU", 2)
    @test  isrepetitive(rna"ACCG", 2)
    @test  isrepetitive(rna"ACGG", 2)
    @test !isrepetitive(rna"ACGUCCGU", 3)
    @test  isrepetitive(rna"ACGCCCGU", 3)

    @test !isrepetitive(mer"ACGU"rna, 2)
    @test  isrepetitive(mer"ACCU"rna, 2)
    @test !isrepetitive(mer"ACCU"rna, 3)
    @test  isrepetitive(mer"ACCCU"rna, 2)
    @test  isrepetitive(mer"ACCCU"rna, 3)

    @test  isrepetitive(aa"")
    @test !isrepetitive(aa"PGQQ")
    @test  isrepetitive(aa"PGQQ", 2)
    @test !isrepetitive(aa"PPQQ", 3)
    @test  isrepetitive(aa"PPPQQ", 3)
end


text:
using PovertyAndInequalityMeasures
using DataFrames
using Test
using CSV
#
# These tests mostly try to replicate examples from
# World Bank 'Handbook on poverty and inequality'
# by Haughton and Khandker
# http://documents.worldbank.org/curated/en/488081468157174849/Handbook-on-poverty-and-inequality
#

# ...
primitive type Amount <: AbstractFloat 64 end

const  TOL = 0.00001; # for nearly equal

"
This just creates an array which is `times` vcat copies of `a`
"
function vcn( a :: Array{Float64}, times :: Int64 )
    nrows = size( a )[1]
    ncols = size( a )[2]
    newrows = nrows*times
    out = zeros( Float64, newrows, ncols )
    p = 0
    for row in 1:nrows
        for i in 1:times
            p += 1
            out[p,:] .= a[row,:]
        end
    end
    out
end

NONE=Array{Symbol,1}();

function removeIgnored( d :: Dict{ Symbol, <:Any}, ignore::Array{Symbol,1} ):: Dict{ Symbol, <:Any}
    for i in ignore
        if haskey( d, i )
            delete!( d, i )
        end
    end
    return d;
end

"
element - by element compare of the type of dicts we use for poverty and inequality output
"
function comparedics( left :: Dict{ Symbol, <:Any}, right :: Dict{ Symbol, <:Any}, ignore::Array{Symbol,1} = NONE ) :: Bool
    left  = removeIgnored( left, ignore )
    right = removeIgnored( right, ignore )
    lk = keys( left )
    if lk != keys( right )
        return false
    end
    for k in lk
        # try catch here in case types are way off
        try
            if !( left[k] ≈ right[k] )
                l = left[k]
                r = right[k]
                print( "comparison failed '$k' : left = $l right = $r")
                return false
            end
       catch e
            return false
       end
    end
    return true
end

@testset "WB Chapter 4 - Poverty " begin
    country_a = [ 1.0 100; 1.0 100; 1 150; 1 150 ]
    country_b = copy( country_a )
    country_b[1:2,2] .= 124
    country_c = [ 1.0 100; 1.0 110; 1 150; 1 160 ]
    #
    # Ch 4 doesn't discuss weighting issues, so
    # we'll add some simple checks for that.
    # a_2 and b_2 should be the same as country_a and _b,
    # but with 2 obs of weight 2 rather than 4 of weight 1
    #
    country_a_2 = [2.0 100; 2.0 150 ]
    country_b_2 = [2.0 124; 2.0 150 ]
    # d should be a big version of a and also produce same result
    country_d = vcn( country_a, 50 )
    # attempt to blow things up with huge a clone
    country_d = vcn( country_c, 100_000 )

    # very unbalanced copy of dataset 1 with 10,000 weight1 1:2 and 2 weight 10,000 7:10
    country_e = vcn( country_c[1:2,:], 10_000 )
    cx = copy(country_c[3:4,:])
    cx[:,1] .= 10_000
    country_e = vcat( country_e, cx )

    line = 125.0
    growth = 0.05

    country_a_pov = make_poverty( country_a, line, growth )
    print("country A " );println( country_a_pov )
    country_a_2_pov = make_poverty( country_a_2, line,growth )
    country_b_pov = make_poverty( country_b, line, growth )
    country_c_pov = make_poverty( country_c, line, growth )
    print("country C " );println( country_c_pov )
    country_d_pov = make_poverty( country_d, line, growth )
    print("country D " );println( country_d_pov )
    country_e_pov = make_poverty( country_e, line, growth )
    print("country E " );println( country_e_pov )

    @test povs_equal( country_a_pov, country_a_2_pov )
    @test povs_equal( country_c_pov, country_e_pov )

    # @test povs_equal( country_c_pov, country_d_pov )
    # test dataframes same as A
    country_a_df = DataFrame( weight=[2.0, 2.0], income=[100.0, 150])
    country_a_pov_df = make_poverty( country_a_df, line, growth, :weight, :income )
    @test povs_equal( country_a_pov_df, country_a_pov )

    # numbers from WP ch. 4
    @test country_a_pov.headcount ≈ 0.5
    @test country_b_pov.headcount ≈ 0.5
    @test country_b_pov.gap ≈ 1.0/250.0
    @test country_c_pov.watts ≈ 0.0877442307
    # some of these are hand-calculations, for from Ada version
    @test isapprox( country_c_pov.gap, 0.080000, atol = TOL )
    @test isapprox( country_c_pov.foster_greer_thorndyke[ 1 ], 0.5000000, atol = TOL ) # pov level
    @test isapprox( country_c_pov.foster_greer_thorndyke[ 2 ], 0.1984059, atol = TOL )
    @test isapprox( country_c_pov.foster_greer_thorndyke[ 3 ], 0.0800000, atol = TOL )
    @test isapprox( country_c_pov.foster_greer_thorndyke[ 4 ], 0.0327530, atol = TOL )
    @test isapprox( country_c_pov.foster_greer_thorndyke[ 5 ], 0.0136000, atol = TOL )
    @test isapprox( country_c_pov.foster_greer_thorndyke[ 6 ], 0.0057192, atol = TOL )
    @test isapprox( country_c_pov.sen , 0.0900000, atol = TOL )
    @test isapprox( country_c_pov.shorrocks , 0.0625000, atol = TOL )
    @test isapprox( country_c_pov.watts , 0.0877442, atol = TOL )
    @test isapprox( country_c_pov.time_to_exit , 1.7548846, atol = TOL )
    @test isapprox( country_c_pov.gini_amongst_poor , 0.0238095, atol = TOL )
    @test isapprox( country_c_pov.poverty_gap_gini , 0.5625000, atol = TOL )
end # poverty testset

#
# reproduce WB Table 6.3 with various cominations ofweights & data
# note table has errors:
#
#    1. uses log10 not ln for  theil##
#    2. has N in wrong place for ge(2) - outside bracket
#
@testset "WB Chapter 6 - Inequality " begin
    c1 = [1.0 10; 1 15; 1 20; 1 25; 1 40; 1 20; 1 30; 1 35; 1 45; 1 90 ]
    # these next are copies of c1 intended
    # to check we haven't screwed up the weighting
    c2 = vcn( c1, 2 )
    c3 = copy( c1 )
    c3[:,1] .= 10_000.0
    c4 = copy( c1 )
    c4[:,1] .= 2.0
    # very unbalanced copy of dataset 1 with 100,000 weight1 1:6 and 4 weight 100,000 7:10
    c64k = vcn( c1[1:6,:], 100_000 )
    cx = copy(c1[7:10,:])
    cx[:,1] .= 100_000
    c64k = vcat( c64k, cx )
    iq1 = make_inequality( c1 )
    iq2 = make_inequality( c2 )
    iq3 = make_inequality( c3 )
    iq4 = make_inequality( c4 )
    iq64k = make_inequality( c64k )
    # weighting and multiplying should make no difference
    println( "iq1");println( iq1 )
    println( "iq2");println( iq2 )
    println( "iq3");println( iq3 )
    println( "iq64k");println( iq64k )

    # test from dataframe
    cdf = DataFrame(income=c1[:,2], weight=c1[:,1])
    iqdf = make_inequality( cdf, :weight, :income)
    @test ineqs_equal( iqdf, iq1 )
    @test ineqs_equal( iq1 , iq64k, include_populations = false )
    @test ineqs_equal( iq1 , iq2, include_populations = false )
    @test ineqs_equal( iq1 , iq3, include_populations = false )
    @test ineqs_equal( iq1 , iq4, include_populations = false )

    @test isapprox( iq1.gini , 0.3272727, atol = TOL )
    @test isapprox( iq1.theil_l, 0.1792203, atol = TOL )
    @test isapprox( iq1.theil_t,  0.1830644, atol = TOL )
    @test isapprox( iq1.generalised_entropy[ 1 ], 0.1883288, atol = TOL )
    @test isapprox( iq1.generalised_entropy[ 2 ], 0.1954897, atol = TOL )
    @test isapprox( iq1.generalised_entropy[ 3 ], 0.2047211, atol = TOL )
    @test isapprox( iq1.generalised_entropy[ 4 ], 0.2162534, atol = TOL )
    @test isapprox( iq1.generalised_entropy[ 5 ], 0.2303812, atol = TOL )
    @test isapprox( iq1.generalised_entropy[ 6 ], 0.2474728, atol = TOL )
    @test isapprox( iq1.atkinson[ 1 ], 0.0446396, atol = TOL )
    @test isapprox( iq1.atkinson[ 2 ], 0.0869155, atol = TOL )
    @test isapprox( iq1.atkinson[ 3 ], 0.1267328, atol = TOL )
    @test isapprox( iq1.atkinson[ 4 ], 0.1640783, atol = TOL )
    @test isapprox( iq1.atkinson[ 5 ], 0.1989991, atol = TOL )
    @test isapprox( iq1.atkinson[ 6 ], 0.2315817, atol = TOL )
    @test isapprox( iq1.atkinson[ 7 ], 0.2619332, atol = TOL )
    @test isapprox( iq1.atkinson[ 8 ], 0.2901688, atol = TOL )
    @test isapprox( iq1.atkinson[ 9 ], 0.3164032, atol = TOL )
    @test isapprox( iq1.hoover, 0.2363636, atol = TOL )
    print( iq1 )
end # inequality testset

hbai_dir = "/mnt/data/hbai/tab/"
if isdir(hbai_dir)
	# HBAI example if available
	    # load each year & jam varnames to lower case
    hbai =  CSV.File("$(hbai_dir)h1819.tab")|>DataFrame
    lcnames = Symbol.(lowercase.(string.(names(hbai))))
    rename!(hbai, lcnames)
    # make scottish subset
    positives = hbai[(hbai.s_oe_ahc .> 0.0 ),:]

    regions = Vector{InequalityMeasures}(undef,0)
    nations = Vector{InequalityMeasures}(undef,0)
    scot = positives[(positives.gvtregn .== 12),:]
    wal  = positives[(positives.gvtregn .== 11),:]
    nire = positives[(positives.gvtregn .== 13),:]
    eng = positives[(positives.gvtregn .< 11),:]

    for reg in 1:13
    	println("region $reg ")
    	if reg != 3
	    	rd = positives[(positives.gvtregn .== reg),:]
	    	ineq_ahc = make_inequality( rd,:gs_newpp,:s_oe_ahc )
	    	push!(regions, ineq_ahc )
	     end
    end
    #  theil decomp isn't exact if <0 incomes included
	uk_ineq_ahc = make_inequality( positives,:gs_newpp,:s_oe_ahc )

	sco_ineq_ahc = make_inequality( scot,:gs_newpp,:s_oe_ahc )
	push!( nations, sco_ineq_ahc )

	wal_ineq_ahc = make_inequality( wal,:gs_newpp,:s_oe_ahc )
	push!( nations, wal_ineq_ahc )

	nire_ineq_ahc = make_inequality( nire,:gs_newpp,:s_oe_ahc )
	push!( nations, nire_ineq_ahc )

	eng_ineq_ahc = make_inequality( eng,:gs_newpp,:s_oe_ahc )
	push!( nations, eng_ineq_ahc )

	dt_regions = add_decomposed_theil( uk_ineq_ahc, regions )
	dt_nations = add_decomposed_theil( uk_ineq_ahc, nations )

    pvline_ahc = positives.mdoeahc[1]*0.6
    pvline_bhc = positives.mdoebhc[1]*0.6

	eng_pov_ahc = make_poverty( eng, pvline_ahc, 0.02, :gs_newpp,:s_oe_ahc )
	sco_pov_ahc = make_poverty( scot, pvline_ahc, 0.02, :gs_newpp,:s_oe_ahc )

end


text:
using Revise
using gtpmip

graph = solve_hfg("test/pddl/hdomain.pddl", "test/pddl/hp1.pddl"; max_levels=100)

# 1


text:
@testset "'MBESS' .............................. " begin

  ss = SampleSpace(5:20, 60)
  prior = p -> Distributions.pdf(Distributions.Beta(2 , 5), p)
  params = MBESS(ss, .2, .3, prior, a = 1, b = 1)
  design = optimaldesign(10, params, solver)

  @test expectedpower(design, prior) > .2

  params = MBESS(ss, .2, .3, prior, a = Inf, b = Inf, targetpower = 0.8)
  design = optimaldesign(10, params, solver)

  z =  z = quadgk(prior, mcrv(parameters(design)), 1)[1]
  tmp = fzero(p -> power(design, p) - .8, 0, 1)
  @test quadgk(prior, tmp, 1)[1]/z > 0.8

end

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
using DiffEqProblemLibrary.DDEProblemLibrary
DDEProblemLibrary.importddeproblems()

using PGFPlotsX

function plot_2delays_analytic(ts)
    prob = DDEProblemLibrary.prob_dde_constant_2delays_scalar
    xmin, xmax = extrema(ts)

    # plot analytic solution
    @pgf ax = Axis(
        {
            xlabel = raw"$t$", ylabel = raw"$x(t)$",
            no_marks, grid = "major",
            xmin = xmin, xmax = xmax
        },
    PlotInc({ thick }, Table(; x = ts, y = prob.f.analytic.(prob.u0, prob.h, prob.p, ts)))
    )

    ax
end

pgfsave(joinpath(@__DIR__, "..", "figures", "generated", "discontinuities_2delays.tex"),
        plot_2delays_analytic(range(0; stop = 1, length = 10_000));
        include_preamble = false)


text:
type FiveDriftingArms <: BanditEnvironment
end

function arms(env::FiveDriftingArms, trial::Int64)
  [
    BernoulliArm(((0 + trial) % 100) / 100.0),
    BernoulliArm(((10 + trial) % 100) / 100.0),
    BernoulliArm(((20 + trial) % 100) / 100.0),
    BernoulliArm((100 - ((1 + trial) % 100)) / 100.0),
    BernoulliArm((100 - trial % 100) / 100.0)
   ]
end

function n_arms(env::FiveDriftingArms)
  5
end

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook


text:
include("../src/QViM.jl")

using .QViM
using LinearAlgebra

quip1 = @quip [
    X >> 1;
    X >> 2 | 1;
]

quip2 = @quip [
    X >> 1;
    X >> 2 | !1;
]

quip3 = @quip [
    X >> 2 | 1;
]

quip4 = @quip [
    1 => C[1]     ;
    X >> 2 | !C[1];
]

quips = [quip1, quip2, quip3, quip4]
for quip in quips
    vm = QVM(quip, 2)
    execute!(vm)
    for (i, j) in enumerate([norm(x)^2 for x in vm.wfn])
        println("$i: $j")
    end
    println()
end


text:
####
#### LinearRegressor
####

"""
    LinearRegressor(; bias::Bool=true)

$LINEAR_DESCR

# Keyword Parameters

- `bias::Bool=true`: if true includes a bias term else fits without bias term.
"""
@mlj_model mutable struct LinearRegressor <: MMI.Deterministic
    bias::Bool = true
end

struct LinearFitresult{T, F<:Real, M<:AbstractArray{F}} <: MMI.MLJType
    sol_matrix::M
    bias::Bool
    names::T
end

_convert(common_type, x::AbstractVector) = convert(AbstractVector{common_type}, x)
_convert(common_type, x::AbstractMatrix) = convert(AbstractMatrix{common_type}, MMI.matrix(x))
matrix_(X::AbstractVector) = X
matrix_(X) = MMI.matrix(X)
_names(y::AbstractVector) = nothing
_names(Y) = collect(MMI.schema(Y).names)

function _matrix(X, target)
    Xmatrix_ = MMI.matrix(X)
    Y_ = matrix_(target)
    common_type = promote_type(eltype(Xmatrix_), eltype(Y_))
    Xmatrix = _convert(common_type, Xmatrix_)
    Y = _convert(common_type, Y_)
    return Xmatrix, Y, _names(target)
end

function MMI.fit(model::LinearRegressor, verbosity::Int, X, y)
    Xmatrix, y_, target_header= _matrix(X, y)
    θ = MS.llsq(Xmatrix, y_; bias=model.bias)
    fitresult = LinearFitresult(θ, model.bias, target_header)
    report = NamedTuple()
    cache = nothing
    return fitresult, cache, report
end

function _regressor_fitted_params(fr::LinearFitresult{Nothing, <:Real, <:AbstractVector})
    return (
        coefficients=fr.sol_matrix[1:end-Int(fr.bias)],
        intercept=ifelse(fr.bias, fr.sol_matrix[end], nothing)
    )
end

function _regressor_fitted_params(fr::LinearFitresult{<:Vector, <:Real, <:AbstractMatrix})
    return (
        coefficients=fr.sol_matrix[1:end-Int(fr.bias), :],
        intercept=fr.bias ? fr.sol_matrix[end, :] : nothing
    )
end

function MMI.fitted_params(::LinearRegressor, fr)
    return _regressor_fitted_params(fr)
end

function _predict_regressor(
    fr::LinearFitresult{Nothing, <:Real, <:AbstractVector},
    Xmat_new::AbstractMatrix,
    prototype
)
    if fr.bias
        return @views Xmat_new * fr.sol_matrix[1:end-1] .+ transpose(fr.sol_matrix[end])
    else
        return @views Xmat_new * fr.sol_matrix[1:end-1]
    end
end

function _predict_regressor(
    fr::LinearFitresult{<:Vector, <:Real, <:AbstractMatrix},
    Xmat_new::AbstractMatrix,
    prototype
)
    if fr.bias
        return MMI.table(
            Xmat_new * @view(fr.sol_matrix[1:end-1, :]) .+ transpose(
                @view(fr.sol_matrix[end, :])
            );
            names=fr.names,
            prototype=prototype
        )
    else
        return MMI.table(
            Xmat_new * @view(fr.sol_matrix[1:end-1, :]);
            names=fr.names,
            prototype=prototype
        )
    end
end

function MMI.predict(::LinearRegressor, fr, Xnew)
    Xmat_new = MMI.matrix(Xnew)
    return _predict_regressor(fr, Xmat_new, Xnew)
end

metadata_model(
    LinearRegressor,
    input=Table(Continuous),
    target=Union{Table(Continuous), AbstractVector{Continuous}},
    weights=false,
    descr=LINEAR_DESCR,
    path="$(PKG).LinearRegressor"
)

####
#### RidgeRegressor
####

_check_typeof_lambda(x)= x isa AbstractVecOrMat || (x isa Real && x ≥ 0)

"""
    RidgeRegressor(; lambda::Union{Real, AbstractVecOrMat}=1.0, bias::Bool=true)

$RIDGE_DESCR

# Keyword Parameters

- `lambda::Union{Real, AbstractVecOrMat}=1.0`: non-negative parameter for the
    regularization strength.
- `bias::Bool=true`: if true includes a bias term else fits without bias term.
"""
@mlj_model mutable struct RidgeRegressor <: MMI.Deterministic
    lambda::Union{Real, AbstractVecOrMat} = 1.0::(_check_typeof_lambda(_))
    bias::Bool = true
end

function MMI.fit(model::RidgeRegressor, verbosity::Int, X, y)
    Xmatrix, y_, target_header = _matrix(X, y)
    θ = MS.ridge(Xmatrix, y_, model.lambda; bias=model.bias)
    fitresult = LinearFitresult(θ, model.bias, target_header)
    report = NamedTuple()
    cache = nothing
    return fitresult, cache, report
end

function MMI.fitted_params(::RidgeRegressor, fr)
    return _regressor_fitted_params(fr)
end

function MMI.predict(::RidgeRegressor, fr, Xnew)
    Xmat_new = MMI.matrix(Xnew)
    return _predict_regressor(fr, Xmat_new, Xnew)
end

metadata_model(
    RidgeRegressor,
    input=Table(Continuous),
    target=Union{Table(Continuous), AbstractVector{Continuous}},
    weights=false,
    descr=RIDGE_DESCR,
    path="$(PKG).RidgeRegressor"
)


text:
# A manufactured solution of a circular wind with constant angular velocity
# on a planetary-sized non-conforming mesh

using OrdinaryDiffEq
using Trixi
using LinearAlgebra

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 1.4
equations = CompressibleEulerEquations3D(gamma)

function initial_condition_circular_wind(x, t, equations::CompressibleEulerEquations3D)
  radius_earth = 6.371229e6
  p = 1e5
  rho = 1.0
  v1 = -10 * x[2] / radius_earth
  v2 = 10 * x[1] / radius_earth
  v3 = 0.0

  return prim2cons(SVector(rho, v1, v2, v3, p), equations)
end

@inline function source_terms_circular_wind(u, x, t, equations::CompressibleEulerEquations3D)
  radius_earth = 6.371229e6
  rho = 1.0

  du1 = 0.0
  du2 = -rho * (10 / radius_earth) * (10 * x[1] / radius_earth)
  du3 = -rho * (10 / radius_earth) * (10 * x[2] / radius_earth)
  du4 = 0.0
  du5 = 0.0

  return SVector(du1, du2, du3, du4, du5)
end

function indicator_test(u::AbstractArray{<:Any,5},
                        mesh, equations, dg::DGSEM, cache;
                        kwargs...)
  alpha = zeros(Int, nelements(dg, cache))

  for element in eachelement(dg, cache)
    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, i, j, k, element)
      lambda, phi, r = cart_to_sphere(x)
      if 0.22 < lambda < 3.3 && 0.45 < phi < 1.3
        alpha[element] = 1
      end
    end
  end

  return alpha
end

function cart_to_sphere(x)
  r = norm(x)
  lambda = atan(x[2], x[1])
  if lambda < 0
    lambda += 2 * pi
  end
  phi = asin(x[3] / r)

  return lambda, phi, r
end

function Trixi.get_element_variables!(element_variables, indicator::typeof(indicator_test), ::AMRCallback)
  return nothing
end

initial_condition = initial_condition_circular_wind

boundary_conditions = Dict(
  :inside  => boundary_condition_slip_wall,
  :outside => boundary_condition_slip_wall
)

# The speed of sound in this example is 374 m/s.
surface_flux = FluxLMARS(374)
# Note that a free stream is not preserved if N < 2 * N_geo, where N is the
# polydeg of the solver and N_geo is the polydeg of the mesh.
# However, the FSP error is negligible in this example.
solver = DGSEM(polydeg=6, surface_flux=surface_flux)

# Use this variable name because it will be doubled in each iteration by the `convergence_test` function.
# It will not be increased if the variable `initial_refinement_level` is used below.
cells_per_dimension = (4, 2)

# Use this variable name for the polydeg of the mesh. It will also be doubled in each iteration.
# Using this variable name will not increase `cells_per_dimension` in each iteration.
initial_refinement_level = 2

trees_per_cube_face = (8, 4) # Use cells_per_dimension here for a convergence test

# Mesh polydeg of 3 (half the solver polydeg) to ensure FSP (see above).
mesh = Trixi.P4estMeshCubedSphere(trees_per_cube_face..., 6.371229e6, 30000.0,
                                  polydeg=initial_refinement_level)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    source_terms=source_terms_circular_wind,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 10.0) # Use a small tspan for performance reasons
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

amr_controller = ControllerThreeLevel(semi, indicator_test,
                                      base_level=0,
                                      max_level=1, max_threshold=0.6)
amr_callback = AMRCallback(semi, amr_controller,
                           interval=0, # Only initial refinement
                           adapt_initial_condition=true,
                           adapt_initial_condition_only_refine=true)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback,
                        save_solution,
                        amr_callback)

###############################################################################
# run the simulation

# Use a Runge-Kutta method with automatic (error based) time step size control
# Enable threading of the RK method for better performance on multiple threads
sol = solve(ode, RDPK3SpFSAL49(thread=OrdinaryDiffEq.True()), abstol=1.0e-6, reltol=1.0e-6,
            save_everystep=false, callback=callbacks);

summary_callback() # print the timer summary


text:
println("Running runtests.jl ...")
using ShapeModels, MultivariateStats, FactCheck, HDF5
FactCheck.setstyle(:compact)

shouldtest(f, a) = length(ARGS) == 0 || in(a, ARGS) ? facts(f, a) : nothing
shouldtestcontext(f, a) = length(ARGS) < 2 || a == ARGS[2] ? context(f, a) : nothing

macro throws_pred(ex) FactCheck.throws_pred(ex) end

landmarks = ShapeModels.examplelandmarks()

shouldtest("basic") do
	a = PCAShapeModel(landmarks)
	@fact size(projection(a.pca)) --> (256,8)
	@fact size(principalvars(a.pca)) --> (8,)
	@fact indim(a.pca) --> 256
	@fact outdim(a.pca) --> 8
end

shouldtest("utils") do
	a = PCAShapeModel(landmarks)
	m = meanshape(a)
	@fact size(m) --> (2,128)
	@fact maximum(abs(mean(m,2))) --> less_than(1e-6)
end

shouldtest("shape") do
	a = PCAShapeModel(landmarks)

    @fact size(shape(a, zeros(nmodes(a)))) --> (2,128)
    if VERSION.minor == 4
        @fact size(a(zeros(nmodes(a)))) --> (2,128)
    end
end

shouldtest("modes") do
	a = PCAShapeModel(landmarks)
	@fact size(modeshapes(a, 1)) --> (2,128,10)
end

println("  ... finished runtests.jl")


text:
include("Tree.jl");

x3 = Variable(UInt(3), semantic=[1.,4.,-2.]);

tree = DivNode(
    Constant(2.) ,
    AddNode(
        DivNode(
            Constant(3.) ,
            SubNode(
                Constant(1.) ,
                DivNode(
                    Constant(2.),
                    x3
                )
            )
        ) ,
        Constant(1.)
    )
);

println("Figure 1: ", string(tree));

targets = [5.,4.,1.];
evaluateTree(tree);

initialEquation = NodeEquation(1., targets, 0., -1.);

calculateEquations!(tree,initialEquation);

println("Table 1:");

nodes, = iterateTree(tree);
for numNode in 1:length(nodes)
    println("-----------------------------------------------------------------------------------------")
    if (isa(nodes[numNode], NonTerminal))
        println("Node ", numNode, ": ", nodes[numNode].name);
    else
        println("Node ", numNode, ": ", string(nodes[numNode]));
    end;
    showEquation(nodes[numNode].equation, nodes[numNode].semantic);
end;

println("-----------------------------------------------------------------------------------------")
println("-----------------------------------------------------------------------------------------")
println("-----------------------------------------------------------------------------------------")

x1 = Variable(UInt(1), semantic=[3.,2.,0.]);
x2 = Variable(UInt(2), semantic=[0.,1.,0.]);
x3 = Variable(UInt(3), semantic=[1.,0.,3.]);
# x4 = Variable(UInt(4), semantic=[2.,0.,4.]);

tree = DivNode(
    Constant(2.) ,
    AddNode(
        DivNode(
            x1,
            Constant(2.)
        ) ,
        MulNode(
            x2,
            SubNode(
                MulNode(
                    x3,
                    # DivNode(
                        Constant(2.) ,
                        # x4
                    # ) ,
                ) ,
                Constant(1.)
            )
        )
    )
);
println("Figure 3: ", string(tree));

evaluateTree(tree);

targets = [5.,4.,1.];
initialEquation = NodeEquation(1., targets, 0., -1.);

calculateEquations!(tree,initialEquation);

println("Table 3: ");

nodes, = iterateTree(tree);
for numNode in 1:length(nodes)
    println("-----------------------------------------------------------------------------------------")
    if (isa(nodes[numNode], NonTerminal))
        println("Node ", numNode, ": ", nodes[numNode].name);
    else
        println("Node ", numNode, ": ", string(nodes[numNode]));
    end;
    showEquation(nodes[numNode].equation, nodes[numNode].semantic);
end;


text:
export ZeroOne
module ZeroOne

export zeroone

function zeroone(s::AbstractString, p = 10^9 + 7)
    nvariants = 1
    nones = 0
    ninvol = 0
    for c in s
        if c === '1'
            nones = mod(nones + nvariants, p)
        elseif c === '0'
            ninvol = mod(ninvol + nones, p)
        elseif c === '?'
            ninvol = mod(ninvol * 2 + nones, p)
            nones = mod(nones * 2 + nvariants, p)
            nvariants = mod(nvariants * 2, p)
        end
    end
    ninvol
end

end # module


text:
### Simple Interfaces (allocate output)
### Double precision
## Type-1

"""
    nufft1d1(xj      :: Array{Float64} or Array{Float32},
             cj      :: Array{ComplexF64} or Array{ComplexF32},
             iflag   :: Integer,
             eps     :: Real,
             ms      :: Integer;
             kwargs...
            ) -> Array{ComplexF64} or Array{ComplexF32}

Compute type-1 1D complex nonuniform FFT.
This computes, to relative precision eps, via a fast algorithm:

              nj
    f(k1) =  SUM c[j] exp(+/-i k1 x(j))  for -ms/2 <= k1 <= (ms-1)/2
             j=1
 # Inputs
  - `xj`      locations of nonuniform sources on interval [-3pi,3pi), length nj
  - `cj`      length-nj complex vector of source strengths. If length(cj)>nj,
            expects a stack of vectors (eg, a nj*ntrans matrix) each of which is
            transformed with the same source locations.
  - `iflag`   if >=0, uses + sign in exponential, otherwise - sign.
  - `eps`     relative precision requested (generally between 1e-15 and 1e-1)
  - `ms`      number of Fourier modes computed, may be even or odd;
            in either case, mode range is integers lying in [-ms/2, (ms-1)/2]
  - kwargs  (optional). See `nufft_opts` and https://finufft.readthedocs.io/en/latest/opts.html
 # Output
  - size `(ms,)` complex vector of Fourier coefficients f, or, if
            `ntrans>1`, matrix of size `(ms,ntrans)`.

"""
function nufft1d1(xj::Array{T},
                  cj::Array{Complex{T}},
                  iflag::Integer,
                  eps::Real,
                  ms::Integer;
                  kwargs...) where T <: finufftReal
    valid_setpts(1,1,xj)
    ntrans = valid_ntr(xj,cj)
    fk = Array{Complex{T}}(undef, ms, ntrans)
    checkkwdtype(T; kwargs...)
    nufft1d1!(xj, cj, iflag, eps, fk; kwargs...,dtype=T)
    return fk
end

"""
    nufft2d1(xj      :: Array{Float64} or Array{Float32}
             yj      :: Array{Float64} or Array{Float32},
             cj      :: Array{ComplexF64} or Array{ComplexF32},
             iflag   :: Integer,
             eps     :: Real,
             ms      :: Integer,
             mt      :: Integer;
             kwargs...
            ) -> Array{ComplexF64}

Compute type-1 2D complex nonuniform FFT.
This computes, to relative precision eps, via a fast algorithm:

                  nj
    f[k1,k2] =   SUM  c[j] exp(+-i (k1 x[j] + k2 y[j]))
                 j=1

    for -ms/2 <= k1 <= (ms-1)/2,  -mt/2 <= k2 <= (mt-1)/2.

 # Inputs
  -  `xj`,`yj`   coordinates of nonuniform sources on the square [-3pi,3pi)^2,
            each a length-nj vector
  -  `cj`      length-nj complex vector of source strengths. If length(cj)>nj,
            expects a stack of vectors (eg, a nj*ntrans matrix) each of which is
            transformed with the same source locations.
  -  `iflag`   if >=0, uses + sign in exponential, otherwise - sign.
   - `eps`     relative precision requested (generally between 1e-15 and 1e-1)
   - `ms`,`mt`   number of Fourier modes requested in x & y; each may be even or odd.
            In either case the mode range is integers lying in [-m/2, (m-1)/2]
  -  kwargs  (optional), see `nufft_opts` and https://finufft.readthedocs.io/en/latest/opts.html
 # Output
  - size `(ms,mt)` complex matrix of Fourier coefficients f
            (ordering given by opts.modeord in each dimension; `ms` fast, `mt` slow),
            or, if `ntrans>1`, a array of size `(ms,mt,ntrans)`.
"""
function nufft2d1(xj      :: Array{T},
                  yj      :: Array{T},
                  cj      :: Array{Complex{T}},
                  iflag   :: Integer,
                  eps     :: Real,
                  ms      :: Integer,
                  mt      :: Integer;
                  kwargs...) where T <: finufftReal
    valid_setpts(1,2,xj,yj)
    ntrans = valid_ntr(xj,cj)
    fk = Array{Complex{T}}(undef, ms, mt, ntrans)
    checkkwdtype(T; kwargs...)
    nufft2d1!(xj, yj, cj, iflag, eps, fk;kwargs...,dtype=T)
    return fk
end

"""
    nufft3d1(xj      :: Array{Float64} or Array{Float32},
             yj      :: Array{Float64} or Array{Float32},
             zj      :: Array{Float64} or Array{Float32},
             cj      :: Array{ComplexF64} or Array{ComplexF32},
             iflag   :: Integer,
             eps     :: Real,
             ms      :: Integer,
             mt      :: Integer,
             mu      :: Integer;
             kwargs...
            ) -> Array{ComplexF64}

Compute type-1 3D complex nonuniform FFT.
This computes, to relative precision eps, via a fast algorithm:

                      nj
    f[k1,k2,k3] =    SUM  c[j] exp(+-i (k1 x[j] + k2 y[j] + k3 z[j]))
                     j=1

    for -ms/2 <= k1 <= (ms-1)/2,  -mt/2 <= k2 <= (mt-1)/2,
        -mu/2 <= k3 <= (mu-1)/2.

 # Inputs
 -   `xj`,`yj`,`zj` coordinates of nonuniform sources on the cube [-3pi,3pi)^3,
             each a length-nj vector
 -   `cj`       length-nj complex vector of source strengths. If length(cj)>nj,
             expects a stack of vectors (eg, a nj*ntrans matrix) each of which is
             transformed with the same source locations.
 -   `iflag`    if >=0, uses + sign in exponential, otherwise - sign.
  -  `eps`      relative precision requested (generally between 1e-15 and 1e-1)
 -   `ms`,`mt`,`mu` number of Fourier modes requested in x,y and z; each may be
             even or odd.
             In either case the mode range is integers lying in [-m/2, (m-1)/2]
  -  kwargs  (optional). See `nufft_opts` and https://finufft.readthedocs.io/en/latest/opts.html
 # Output
  - size `(ms,mt,mu)` complex array of Fourier coefficients f
            (ordering given by opts.modeord in each dimension; `ms` fastest, `mu`
            slowest), or, if `ntrans>1`, a 4D array of size `(ms,mt,mu,ntrans)`.
"""
function nufft3d1(xj      :: Array{T},
                  yj      :: Array{T},
                  zj      :: Array{T},
                  cj      :: Array{Complex{T}},
                  iflag   :: Integer,
                  eps     :: Real,
                  ms      :: Integer,
                  mt      :: Integer,
                  mu      :: Integer;
                  kwargs...)  where T <: finufftReal
    valid_setpts(1,3,xj,yj,zj)
    ntrans = valid_ntr(xj,cj)
    fk = Array{Complex{T}}(undef, ms, mt, mu, ntrans)
    checkkwdtype(T; kwargs...)
    nufft3d1!(xj, yj, zj, cj, iflag, eps, fk;kwargs...,dtype=T)
    return fk
end

## Type-2

"""
    nufft1d2(xj      :: Array{Float64} or Array{Float32},
             iflag   :: Integer,
             eps     :: Real,
             fk      :: Array{ComplexF64} or Array{ComplexF32};
             kwargs...
            ) -> Array{ComplexF64}

Compute type-2 1D complex nonuniform FFT.
This computes, to relative precision eps, via a fast algorithm:

    c[j] = SUM   f[k1] exp(+/-i k1 x[j])      for j = 1,...,nj
            k1
     where sum is over -ms/2 <= k1 <= (ms-1)/2.

 # Input
  -  `xj`      location of nonuniform targets on interval [-3pi,3pi), length nj
   - `fk`      complex Fourier coefficients. If a vector, length sets `ms`
            (with mode ordering given by opts.modeord). If a matrix, each
            column is transformed with the same nonuniform targets.
  -  `iflag`   if >=0, uses + sign in exponential, otherwise - sign.
  -  `eps`     relative precision requested (generally between 1e-15 and 1e-1)
  -  kwargs  (optional). See `nufft_opts` and https://finufft.readthedocs.io/en/latest/opts.html
 # Output
  - complex `(nj,)` vector c of answers at targets, or,
           if `ntrans>1`, matrix of size `(nj,ntrans)`.
"""
function nufft1d2(xj      :: Array{T},
                  iflag   :: Integer,
                  eps     :: Real,
                  fk      :: Array{Complex{T}};
                  kwargs...) where T <: finufftReal
    (nj, nk) = valid_setpts(2,1,xj)
    (ms, ntrans) = get_nmodes_from_fk(1,fk)
    cj = Array{Complex{T}}(undef, nj, ntrans)
    checkkwdtype(T; kwargs...)
    nufft1d2!(xj, cj, iflag, eps, fk;kwargs...,dtype=T)
    return cj
end

"""
    nufft2d2(xj      :: Array{Float64} or Array{Float32},
             yj      :: Array{Float64} or Array{Float32},
             iflag   :: Integer,
             eps     :: Real,
             fk      :: Array{ComplexF64} or Array{ComplexF32};
             kwargs...
            ) -> Array{ComplexF64}

Compute type-2 2D complex nonuniform FFT.
This computes, to relative precision eps, via a fast algorithm:

    c[j] =  SUM   f[k1,k2] exp(+/-i (k1 x[j] + k2 y[j]))  for j = 1,..,nj
           k1,k2
     where sum is over -ms/2 <= k1 <= (ms-1)/2, -mt/2 <= k2 <= (mt-1)/2,

 # Inputs
   -  `xj`,`yj`   coordinates of nonuniform targets on the square [-3pi,3pi)^2,
            each a vector of length nj
  -  `fk`      complex Fourier coefficient matrix, whose size determines (ms,mt).
            (Mode ordering given by opts.modeord, in each dimension.)
            If a 3D array, 3rd dimension sets `ntrans`, and each of `ntrans`
            matrices is transformed with the same nonuniform targets.
  -  `iflag`   if >=0, uses + sign in exponential, otherwise - sign.
  -  `eps`     relative precision requested (generally between 1e-15 and 1e-1)
 -   kwargs  (optional). See `nufft_opts` and https://finufft.readthedocs.io/en/latest/opts.html
 # Output
   - complex size `(nj,)` vector c of answers at targets, or,
            if `ntrans>1`, matrix of size `(nj,ntrans)`.

"""
function nufft2d2(xj      :: Array{T},
                  yj      :: Array{T},
                  iflag   :: Integer,
                  eps     :: Real,
                  fk      :: Array{Complex{T}};
                  kwargs...) where T <: finufftReal
    (nj, nk) = valid_setpts(2,2,xj,yj)
    (ms, mt, ntrans) = get_nmodes_from_fk(2,fk)
    cj = Array{Complex{T}}(undef, nj, ntrans)
    checkkwdtype(T; kwargs...)
    nufft2d2!(xj, yj, cj, iflag, eps, fk;kwargs...,dtype=T)
    return cj
end

"""
    nufft3d2(xj      :: Array{Float64} or Array{Float32},
             yj      :: Array{Float64} or Array{Float32},
             zj      :: Array{Float64} or Array{Float32},
             iflag   :: Integer,
             eps     :: Real,
             fk      :: Array{ComplexF64} or Array{ComplexF32};
             kwargs...
            ) -> Array{ComplexF64}

Compute type-2 3D complex nonuniform FFT.
This computes, to relative precision eps, via a fast algorithm:

    c[j] =   SUM   f[k1,k2,k3] exp(+/-i (k1 x[j] + k2 y[j] + k3 z[j]))
           k1,k2,k3
                           for j = 1,..,nj
     where sum is over -ms/2 <= k1 <= (ms-1)/2, -mt/2 <= k2 <= (mt-1)/2,
                      -mu/2 <= k3 <= (mu-1)/2.

 # Inputs
  -  `xj`,`yj`,`zj` coordinates of nonuniform targets on the cube [-3pi,3pi)^3,
             each a vector of length nj
  -  `fk`       complex Fourier coefficient array, whose size sets `(ms,mt,mu)`.
             (Mode ordering given by opts.modeord, in each dimension.)
             If a 4D array, 4th dimension sets `ntrans`, and each of `ntrans`
             3D arrays is transformed with the same nonuniform targets.
  -  `iflag`    if >=0, uses + sign in exponential, otherwise - sign.
 -   `eps`      relative precision requested (generally between 1e-15 and 1e-1)
   - kwargs   (optional). See `nufft_opts` and https://finufft.readthedocs.io/en/latest/opts.html
 #  Output
   - complex vector c of size `(nj,)` giving answers at targets, or,
             if `ntrans>1`, matrix of size `(nj,ntrans)`.
"""
function nufft3d2(xj      :: Array{T},
                  yj      :: Array{T},
                  zj      :: Array{T},
                  iflag   :: Integer,
                  eps     :: Real,
                  fk      :: Array{Complex{T}};
                  kwargs...) where T <: finufftReal
    (nj, nk) = valid_setpts(2,3,xj,yj,zj)
    (ms, mt, mu, ntrans) = get_nmodes_from_fk(3,fk)
    cj = Array{Complex{T}}(undef, nj, ntrans)
    checkkwdtype(T; kwargs...)
    nufft3d2!(xj, yj, zj, cj, iflag, eps, fk;kwargs...,dtype=T)
    return cj
end

## Type-3

"""
    nufft1d3(xj      :: Array{Float64} or Array{Float32},
             cj      :: Array{ComplexF64} or Array{ComplexF32},
             iflag   :: Integer,
             eps     :: Real,
             sk      :: Array{Float64} or Array{Float32};
             kwargs...
            ) -> Array{ComplexF64}

Compute type-3 1D complex nonuniform FFT.
This computes, to relative precision eps, via a fast algorithm:

             nj
    f[k]  =  SUM   c[j] exp(+-i s[k] x[j]),      for k = 1, ..., nk
             j=1

 # Inputs
 -   `xj`       locations of nonuniform sources on R (real line), length-nj vector.
 -   `cj`       length-nj complex vector of source strengths. If length(cj)>nj,
             expects a size `(nj,ntrans)` matrix each column of which is
             transformed with the same source and target locations.
 -    `iflag`    if >=0, uses + sign in exponential, otherwise - sign.
 - `eps`      relative precision requested (generally between 1e-15 and 1e-1)
 -  `sk`       frequency locations of nonuniform targets on R, length-nk vector.
  -  kwargs   (optional). See `nufft_opts` and https://finufft.readthedocs.io/en/latest/opts.html
 #  Output
  - complex vector f size '(nk,)` of values at targets, or, if `ntrans>1`,
             a matrix of size `(nk,ntrans)`
"""
function nufft1d3(xj      :: Array{T},
                  cj      :: Array{Complex{T}},
                  iflag   :: Integer,
                  eps     :: Real,
                  sk      :: Array{T};
                  kwargs...) where T <: finufftReal
    (nj, nk) = valid_setpts(3,1,xj,T[],T[],sk)
    ntrans = valid_ntr(xj,cj)
    fk = Array{Complex{T}}(undef, nk, ntrans)
    checkkwdtype(T; kwargs...)
    nufft1d3!(xj, cj, iflag, eps, sk, fk;kwargs...,dtype=T)
    return fk
end

"""
    nufft2d3(xj      :: Array{Float64} or Array{Float32},
             yj      :: Array{Float64} or Array{Float32},
             cj      :: Array{ComplexF64} or Array{ComplexF32},
             iflag   :: Integer,
             eps     :: Real,
             sk      :: Array{Float64} or Array{Float32},
             tk      :: Array{Float64} or Array{Float32};
             kwargs...
            ) -> Array{ComplexF64}

Compute type-3 2D complex nonuniform FFT.
This computes, to relative precision eps, via a fast algorithm:

             nj
    f[k]  =  SUM   c[j] exp(+-i (s[k] x[j] + t[k] y[j])),  for k = 1, ..., nk
             j=1

 # Inputs
  -  `xj`,`yj`    coordinates of nonuniform sources in R^2, each a length-nj vector.
   - `cj`    complex vector `(nj,)` of source strengths. If length(cj)>nj,
             expects a `(nj,ntrans)` matrix, each column of which is
             transformed with the same source and target locations.
 -   `iflag`    if >=0, uses + sign in exponential, otherwise - sign.
 -   `eps`      relative precision requested (generally between 1e-15 and 1e-1)
  -  `sk`,`tk`    frequency coordinates of nonuniform targets in R^2,
             each a length-nk vector.
 -   kwargs   (optional). See `nufft_opts` and https://finufft.readthedocs.io/en/latest/opts.html
 # Output
  - complex vector size `(nk,)` of values at targets, or, if `ntrans>1`,
             a matrix of size `(nk,ntrans)`
"""
function nufft2d3(xj      :: Array{T},
                  yj      :: Array{T},
                  cj      :: Array{Complex{T}},
                  iflag   :: Integer,
                  eps     :: Real,
                  sk      :: Array{T},
                  tk      :: Array{T};
                  kwargs...) where T <: finufftReal
    (nj, nk) = valid_setpts(3,2,xj,yj,T[],sk,tk)
    ntrans = valid_ntr(xj,cj)
    fk = Array{Complex{T}}(undef, nk, ntrans)
    checkkwdtype(T; kwargs...)
    nufft2d3!(xj, yj, cj, iflag, eps, sk, tk, fk;kwargs...,dtype=T)
    return fk
end

"""
    nufft3d3(xj      :: Array{Float64} or Array{Float32},
             yj      :: Array{Float64} or Array{Float32},
             zj      :: Array{Float64} or Array{Float32},
             cj      :: Array{ComplexF64} or Array{ComplexF32},
             iflag   :: Integer,
             eps     :: Real,
             sk      :: Array{Float64} or Array{Float32},
             tk      :: Array{Float64} or Array{Float32},
             uk      :: Array{Float64} or Array{Float32};
             kwargs...
            ) -> Array{ComplexF64}

Compute type-3 3D complex nonuniform FFT.
This computes, to relative precision eps, via a fast algorithm:

             nj
    f[k]  =  SUM   c[j] exp(+-i (s[k] x[j] + t[k] y[j] + u[k] z[j])),
             j=1
                             for k = 1, ..., nk

 # Inputs
  -  `xj`,`yj`,`zj` coordinates of nonuniform sources in R^3, each a length-nj vector.
 -   `cj`     complex `(nj,)` vector of source strengths. If length(cj)>nj,
             expects a '(nj,ntrans)' matrix, each column of which is
             transformed with the same source and target locations.
 -   `iflag`    if >=0, uses + sign in exponential, otherwise - sign.
 -   `eps`      relative precision requested (generally between 1e-15 and 1e-1)
  -  `sk`,`tk,`uk` frequency coordinates of nonuniform targets in R^3,
             each a length-nk vector.
  -  kwargs   (optional). See `nufft_opts` and https://finufft.readthedocs.io/en/latest/opts.html
 #  Output
  - size `(nk,)` complex vector f values at targets, or, if `ntrans>1`,
             a matrix of size `(nk,ntrans)`
"""
function nufft3d3(xj      :: Array{T},
                  yj      :: Array{T},
                  zj      :: Array{T},
                  cj      :: Array{Complex{T}},
                  iflag   :: Integer,
                  eps     :: Real,
                  sk      :: Array{T},
                  tk      :: Array{T},
                  uk      :: Array{T};
                  kwargs...) where T <: finufftReal
    (nj, nk) = valid_setpts(3,3,xj,yj,zj,sk,tk,uk)
    ntrans = valid_ntr(xj,cj)
    fk = Array{Complex{T}}(undef, nk, ntrans)
    checkkwdtype(T; kwargs...)
    nufft3d3!(xj, yj, zj, cj, iflag, eps, sk, tk, uk, fk;kwargs...,dtype=T)
    return fk
end

### Direct interfaces (No allocation)
### Double precision
## 1D

"""
    nufft1d1!(xj      :: Array{Float64} or Array{Float32},
              cj      :: Array{ComplexF64} or Array{ComplexF32},
              iflag   :: Integer,
              eps     :: Real,
              fk      :: Array{ComplexF64} or Array{ComplexF32};
              kwargs...
            )

Compute type-1 1D complex nonuniform FFT. Output written to `fk`. See `nufft1d1`.
"""
function nufft1d1!(xj      :: Array{T},
                   cj      :: Array{Complex{T}},
                   iflag   :: Integer,
                   eps     :: Real,
                   fk      :: Array{Complex{T}};
                   kwargs...) where T <: finufftReal
    valid_setpts(1,1,xj)
    ntrans = valid_ntr(xj,cj)
    (ms, ntrans_fk) = get_nmodes_from_fk(1,fk)

    checkkwdtype(T; kwargs...)
    plan = finufft_makeplan(1,[ms;],iflag,ntrans,eps;kwargs...)
    finufft_setpts(plan,xj)
    finufft_exec!(plan,cj,fk)
    ret = finufft_destroy(plan)
    check_ret(ret)
end

"""
    nufft1d2!(xj      :: Array{Float64} or Array{Float32},
              cj      :: Array{ComplexF64} or Array{ComplexF32},
              iflag   :: Integer,
              eps     :: Real,
              fk      :: Array{ComplexF64} or Array{ComplexF32};
              kwargs...
            )

Compute type-2 1D complex nonuniform FFT. Output written to `cj`. See `nufft1d2`.
"""
function nufft1d2!(xj      :: Array{T},
                   cj      :: Array{Complex{T}},
                   iflag   :: Integer,
                   eps     :: Real,
                   fk      :: Array{Complex{T}};
                   kwargs...) where T <: finufftReal
    (nj, nk) = valid_setpts(2,1,xj)
    (ms, ntrans) = get_nmodes_from_fk(1,fk)

    checkkwdtype(T; kwargs...)
    plan = finufft_makeplan(2,[ms;],iflag,ntrans,eps;kwargs...)
    finufft_setpts(plan,xj)
    finufft_exec!(plan,fk,cj)
    ret = finufft_destroy(plan)
    check_ret(ret)
end

"""
    nufft1d3!(xj      :: Array{Float64} or Array{Float32},
              cj      :: Array{ComplexF64} or Array{ComplexF32},
              iflag   :: Integer,
              eps     :: Real,
              sk      :: Array{Float64} or Array{Float32},
              fk      :: Array{ComplexF64} or Array{ComplexF32};
              kwargs...
             )

Compute type-3 1D complex nonuniform FFT. Output written to `fk`. See `nufft1d3`.
"""
function nufft1d3!(xj      :: Array{T},
                   cj      :: Array{Complex{T}},
                   iflag   :: Integer,
                   eps     :: Real,
                   sk      :: Array{T},
                   fk      :: Array{Complex{T}};
                   kwargs...) where T <: finufftReal
    (nj, nk) = valid_setpts(3,1,xj,T[],T[],sk)
    ntrans = valid_ntr(xj,cj)

    checkkwdtype(T; kwargs...)
    plan = finufft_makeplan(3,1,iflag,ntrans,eps;kwargs...)
    finufft_setpts(plan,xj,T[],T[],sk)
    finufft_exec!(plan,cj,fk)
    ret = finufft_destroy(plan)
    check_ret(ret)
end

## 2D

"""
    nufft2d1!(xj      :: Array{Float64} or Array{Float32},
              yj      :: Array{Float64} or Array{Float32},
              cj      :: Array{ComplexF64} or Array{ComplexF32},
              iflag   :: Integer,
              eps     :: Real,
              fk      :: Array{ComplexF64} or Array{ComplexF32};
              kwargs...
            )

Compute type-1 2D complex nonuniform FFT. Output written to `fk`. See `nufft2d1`.
"""
function nufft2d1!(xj      :: Array{T},
                   yj      :: Array{T},
                   cj      :: Array{Complex{T}},
                   iflag   :: Integer,
                   eps     :: Real,
                   fk      :: Array{Complex{T}};
                   kwargs...) where T <: finufftReal
    valid_setpts(1,2,xj,yj)
    ntrans = valid_ntr(xj,cj)
    (ms, mt, ntrans_fk) = get_nmodes_from_fk(2,fk)
    @assert ntrans==ntrans_fk

    checkkwdtype(T; kwargs...)
    plan = finufft_makeplan(1,[ms;mt],iflag,ntrans,eps;kwargs...)
    finufft_setpts(plan,xj,yj)
    finufft_exec!(plan,cj,fk)
    ret = finufft_destroy(plan)
    check_ret(ret)
end

"""
    nufft2d2!(xj      :: Array{Float64} or Array{Float32},
              yj      :: Array{Float64} or Array{Float32},
              cj      :: Array{ComplexF64} or Array{ComplexF32},
              iflag   :: Integer,
              eps     :: Real,
              fk      :: Array{ComplexF64} or Array{ComplexF32};
              kwargs...
            )

Compute type-2 2D complex nonuniform FFT. Output written to `cj`. See `nufft2d2`.
"""
function nufft2d2!(xj      :: Array{T},
                   yj      :: Array{T},
                   cj      :: Array{Complex{T}},
                   iflag   :: Integer,
                   eps     :: Real,
                   fk      :: Array{Complex{T}};
                   kwargs...) where T <: finufftReal
    (nj, nk) = valid_setpts(1,2,xj,yj)
    (ms, mt, ntrans) = get_nmodes_from_fk(2,fk)

    checkkwdtype(T; kwargs...)
    plan = finufft_makeplan(2,[ms;mt],iflag,ntrans,eps;kwargs...)
    finufft_setpts(plan,xj,yj)
    finufft_exec!(plan,fk,cj)
    ret = finufft_destroy(plan)
    check_ret(ret)
end

"""
    nufft2d3!(xj      :: Array{Float64} or Array{Float32},
              yj      :: Array{Float64} or Array{Float32},
              cj      :: Array{ComplexF64} or Array{ComplexF32},
              iflag   :: Integer,
              eps     :: Real,
              sk      :: Array{Float64} or Array{Float32},
              tk      :: Array{Float64} or Array{Float32},
              fk      :: Array{ComplexF64} or Array{ComplexF32};
              kwargs...
             )

Compute type-3 2D complex nonuniform FFT. Output written to 'fk'. See `nufft2d3`.
"""
function nufft2d3!(xj      :: Array{T},
                   yj      :: Array{T},
                   cj      :: Array{Complex{T}},
                   iflag   :: Integer,
                   eps     :: Real,
                   sk      :: Array{T},
                   tk      :: Array{T},
                   fk      :: Array{Complex{T}};
                   kwargs...) where T <: finufftReal
    (nj, nk) = valid_setpts(3,2,xj,yj,T[],sk,tk)
    ntrans = valid_ntr(xj,cj)

    checkkwdtype(T; kwargs...)
    plan = finufft_makeplan(3,2,iflag,ntrans,eps;kwargs...)
    finufft_setpts(plan,xj,yj,T[],sk,tk)
    finufft_exec!(plan,cj,fk)
    ret = finufft_destroy(plan)
    check_ret(ret)
end

## 3D

"""
    nufft3d1!(xj      :: Array{Float64} or Array{Float32},
              yj      :: Array{Float64} or Array{Float32},
              zj      :: Array{Float64} or Array{Float32},
              cj      :: Array{ComplexF64} or Array{ComplexF32},
              iflag   :: Integer,
              eps     :: Real,
              fk      :: Array{ComplexF64} or Array{ComplexF32};
              kwargs...
            )

Compute type-1 3D complex nonuniform FFT. Output written to `fk`. See `nufft3d1`.
"""
function nufft3d1!(xj      :: Array{T},
                   yj      :: Array{T},
                   zj      :: Array{T},
                   cj      :: Array{Complex{T}},
                   iflag   :: Integer,
                   eps     :: Real,
                   fk      :: Array{Complex{T}};
                   kwargs...) where T <: finufftReal
    valid_setpts(1,3,xj,yj,zj)
    ntrans = valid_ntr(xj,cj)
    (ms, mt, mu, ntrans_fk) = get_nmodes_from_fk(3,fk)
    @assert ntrans == ntrans_fk

    checkkwdtype(T; kwargs...)
    plan = finufft_makeplan(1,[ms;mt;mu],iflag,ntrans,eps;kwargs...)
    finufft_setpts(plan,xj,yj,zj)
    finufft_exec!(plan,cj,fk)
    ret = finufft_destroy(plan)
    check_ret(ret)
end

"""
    nufft3d2!(xj      :: Array{Float64} or Array{Float32},
              yj      :: Array{Float64} or Array{Float32},
              zj      :: Array{Float64} or Array{Float32},
              cj      :: Array{ComplexF64} or Array{ComplexF32},
              iflag   :: Integer,
              eps     :: Real,
              fk      :: Array{ComplexF64} or Array{ComplexF32};
              kwargs...
            )

Compute type-2 3D complex nonuniform FFT. Output written to `cj`. See `nufft3d2`.
"""
function nufft3d2!(xj      :: Array{T},
                   yj      :: Array{T},
                   zj      :: Array{T},
                   cj      :: Array{Complex{T}},
                   iflag   :: Integer,
                   eps     :: Real,
                   fk      :: Array{Complex{T}};
                   kwargs...) where T <: finufftReal
    (nj, nk) = valid_setpts(2,3,xj,yj,zj)
    (ms, mt, mu, ntrans) = get_nmodes_from_fk(3,fk)

    checkkwdtype(T; kwargs...)
    plan = finufft_makeplan(2,[ms;mt;mu],iflag,ntrans,eps;kwargs...)
    finufft_setpts(plan,xj,yj,zj)
    finufft_exec!(plan,fk,cj)
    ret = finufft_destroy(plan)
    check_ret(ret)
end

"""
    nufft3d3!(xj      :: Array{Float64} or Array{Float32},
              yj      :: Array{Float64} or Array{Float32},
              zj      :: Array{Float64} or Array{Float32},
              cj      :: Array{ComplexF64} or Array{ComplexF32},
              iflag   :: Integer,
              eps     :: Real,
              sk      :: Array{Float64} or Array{Float32},
              tk      :: Array{Float64} or Array{Float32},
              uk      :: Array{Float64} or Array{Float32},
              fk      :: Array{ComplexF64} or Array{ComplexF32};
              kwargs...
             )

Compute type-3 3D complex nonuniform FFT. Output written to `fk`. See `nufft3d3`.
"""
function nufft3d3!(xj      :: Array{T},
                   yj      :: Array{T},
                   zj      :: Array{T},
                   cj      :: Array{Complex{T}},
                   iflag   :: Integer,
                   eps     :: Real,
                   sk      :: Array{T},
                   tk      :: Array{T},
                   uk      :: Array{T},
                   fk      :: Array{Complex{T}};
                   kwargs...) where T <: finufftReal
    (nj, nk) = valid_setpts(3,3,xj,yj,zj,sk,tk,uk)
    ntrans = valid_ntr(xj,cj)

    checkkwdtype(T; kwargs...)
    plan = finufft_makeplan(3,3,iflag,ntrans,eps;kwargs...)
    finufft_setpts(plan,xj,yj,zj,sk,tk,uk)
    finufft_exec!(plan,cj,fk)
    ret = finufft_destroy(plan)
    check_ret(ret)
end


text:
function evaluate_to_array(backend::ArrayBackend, x::Number, target_dimensions; seed=nothing)
    evaluate(backend, convert(FieldExpr, x), target_dimensions)
end

function evaluate_to_array(backend::ArrayBackend, fexpr::FieldExpr, target_dimensions; seed=nothing)
    af = wrapped_array_function([fexpr => target_dimensions], [])
    cf = compile_wrapped_array_function(backend, af)
    cf(seed=seed)
end

function evaluate(backend::ArrayBackend, fexpr::FieldExpr, target_dimensions=nothing; seed=nothing)

    if target_dimensions === nothing
        target_dimensions = collect(known_dimensions_expr(fexpr))
    end

    if !is_stochastic(fexpr)
        dims = target_dimensions
    else
        dims = [target_dimensions..., STOCHASTIC_DIMENSION]
    end

    ffunc = field_function([fexpr => dims], [])
    cffunc = compile_field_function(backend, ffunc)
    cffunc()
end

evaluate(args...; kwds...) = evaluate(default_array_backend(), args...; kwds...)

function known_field{T <: KnownDimension}(x::BackendArray, dimensions::Array{T, 1})
    @assert ndims(x) == length(dimensions)
    sdims = find(dimensions .== STOCHASTIC_DIMENSION)
    if length(sdims) == 0
        regular_dims = dimensions
        stochastic_info = StochasticInfo()
    else
        regular_dims = [dimensions[1:sdims[1]-1]..., dimensions[sdims[1]+1:end]...]
        stochastic_info = StochasticInfo(size(x, sdims[1]))
    end
    KnownField(regular_dims, x, stochastic_info)
end


text:
function cstrt_datta_ngon(n)
if mod(log2(n),1)==0 && n>=16
    # initialization
    d = atan(tan(2*pi/n)*tan(pi/n))-asin(sin(2*pi/n)^2/sqrt(2+2*cos(2*pi/n)*cos(4*pi/n)));
    x = zeros(Int(3*n/8));
    for i = 1:Int(3*n/8)
        if mod(i,3) == 2
            x[i] = 2*(pi/n + (-1)^(i-1)*d);
        else
            x[i] = pi/n + (-1)^(i-1)*d;
        end
    end
    p = zeros(Int(3*n/8-1)); q = zeros(Int(3*n/8-1));
    p[1] = sin(x[1]); q[1] = cos(x[1]);
    for i = 2:Int(3*n/8-1)
        p[i] = p[i-1] - (-1)^i*sin(sum(x[1:i]));
        q[i] = q[i-1] - (-1)^i*cos(sum(x[1:i]));
    end
    # construction
    a = zeros(n-1); b = zeros(n-1);
    a[Int(n/2-1)] = p[1]; a[Int(n/2+1)] = -a[Int(n/2-1)];
    b[Int(n/2-1)] = q[1]; b[Int(n/2+1)] = b[Int(n/2-1)];
    a[n-1] = a[Int(n/2-1)]-sin(x[1]+x[2]/2); a[1] = -a[n-1];
    b[n-1] = b[Int(n/2-1)]-cos(x[1]+x[2]/2); b[1] = b[n-1];
    a[n-2] = p[2]; a[2] = -a[n-2];
    b[n-2] = q[2]; b[2] = b[n-2];
    for i = 1:Int(n/8-1)
        a[Int(mod(n/2-1+i*(3*n/2-2),n))] = p[3*i+1];
        a[Int(n-mod(n/2-1+i*(3*n/2-2),n))] = -a[Int(mod(n/2-1+i*(3*n/2-2),n))];
        b[Int(mod(n/2-1+i*(3*n/2-2),n))] = q[3*i+1];
        b[Int(n-mod(n/2-1+i*(3*n/2-2),n))] = b[Int(mod(n/2-1+i*(3*n/2-2),n))];
        a[Int(mod(n-1+i*(3*n/2-2),n))] = a[Int(mod(n/2-1+i*(3*n/2-2),n))]-(-1)^i*sin(sum(x[1:3*i+1])+x[3*i+2]/2);
        a[Int(n-mod(n-1+i*(3*n/2-2),n))] = -a[Int(mod(n-1+i*(3*n/2-2),n))];
        b[Int(mod(n-1+i*(3*n/2-2),n))] = b[Int(mod(n/2-1+i*(3*n/2-2),n))]-(-1)^i*cos(sum(x[1:3*i+1])+x[3*i+2]/2);
        b[Int(n-mod(n-1+i*(3*n/2-2),n))] = b[Int(mod(n-1+i*(3*n/2-2),n))];
        a[Int(mod(n-2+i*(3*n/2-2),n))] = p[3*i+2];
        a[Int(n-mod(n-2+i*(3*n/2-2),n))] = -a[Int(mod(n-2+i*(3*n/2-2),n))];
        b[Int(mod(n-2+i*(3*n/2-2),n))] = q[3*i+2];
        b[Int(n-mod(n-2+i*(3*n/2-2),n))] = b[Int(mod(n-2+i*(3*n/2-2),n))];
        a[Int(mod(n+i*(3*n/2-2),n))] = p[3*i];
        a[Int(n-mod(n+i*(3*n/2-2),n))] = -a[Int(mod(n+i*(3*n/2-2),n))];
        b[Int(mod(n+i*(3*n/2-2),n))] = q[3*i];
        b[Int(n-mod(n+i*(3*n/2-2),n))] = b[Int(mod(n+i*(3*n/2-2),n))];
    end
    a[Int(n/2)] = 0; b[Int(n/2)] = 1;
    return(a,b)
end
end

max_stars_repo_name: cbingane/optigon

max_issues_repo_name: cbingane/optigon

max_forks_repo_name: cbingane/optigon


text:
"""
$(SIGNATURES)

Convenience function to call the base markdown to html converter on "simple" strings (i.e. strings
that don't need to be further considered and don't contain anything else than markdown tokens).
The boolean `stripp` indicates whether to remove the inserted `<p>` and `</p>` by the base markdown
processor, this is relevant for things that are parsed within latex commands etc.
"""
function md2html(ss::AS; stripp::Bool=false)::AS
    # if there's nothing, return that...
    isempty(ss) && return ss
    # Use Julia's Markdown parser followed by Julia's MD->HTML conversion
    partial = ss |> fix_inserts |> Markdown.parse
    # take over from the parsing of indented blocks
    for (i, c) in enumerate(partial.content)
        c isa Markdown.Code || continue
        partial.content[i] = Markdown.Paragraph(Any[c.code])
    end
    # Use Julia's MD->HTML conversion
    partial = partial |> Markdown.html
    # Markdown.html transforms {{ with HTML entities but we don't want that
    partial = replace(partial, r"{{" => "{{")
    partial = replace(partial, r"}}" => "}}")
    # In some cases, base converter adds
...
\n which we might not want
    stripp || return partial
    if startswith(partial, "
") && endswith(partial, "
\n")
        partial = chop(partial, head=3, tail=5)
    end
    return partial
end

"""
$(SIGNATURES)

Convenience function to check if `idx` is smaller than the length of `v`, if it is, then return the starting point of `v[idx]` (via `from`), otherwise return `BIG_INT`.
"""
from_ifsmaller(v::Vector, idx::Int, len::Int)::Int = (idx > len) ? BIG_INT : from(v[idx])

"""
$(SIGNATURES)

Since divs are recursively processed, once they've been found, everything inside them needs to be
deactivated and left for further re-processing to avoid double inclusion.
"""
function deactivate_divs(blocks::Vector{OCBlock})::Vector{OCBlock}
    active_blocks = ones(Bool, length(blocks))
    for (i, β) ∈ enumerate(blocks)
        fromβ, toβ = from(β), to(β)
        active_blocks[i] || continue
        if β.name == :DIV
            innerblocks = findall(b -> (fromβ < from(b) < toβ), blocks)
            active_blocks[innerblocks] .= false
        end
    end
    return blocks[active_blocks]
end

"""
$(SIGNATURES)

The insertion token have whitespaces around them: ` ##FDINSERT## `, this mostly helps but causes
a problem when combined with italic or bold markdown mode since `_blah_` works but not `_ blah _`.
This function looks for any occurrence of `[\\*_] ##FDINSERT##` or the opposite and removes the
extraneous whitespace.
"""
fix_inserts(s::AS)::String =
    replace(replace(s, r"([\*_]) ##FDINSERT##" => s"\1##FDINSERT##"),
                       r"##FDINSERT## ([\*_])" => s"##FDINSERT##\1")

"""
$(SIGNATURES)

Takes a list of tokens and deactivate tokens that happen to be in a multi-line
md-def. Used in [`convert_md`](@ref).
"""
function preprocess_candidate_mddefs!(tokens::Vector{Token})
    isempty(tokens) && return nothing
    # process:
    # 1. find a MD_DEF_OPEN token
    # 2. Look for the first LINE_RETURN (proper)
    # 3. try to parse the content with Meta.parse. If it completely fails,
    # error, otherwise consider the  span of the first ok expression and
    # discard all tokens within its span leaving effectively just the
    # opening MD_DEF_OPEN and closing LINE_RETURN
    from_to_list = Pair{Int,Int}[]
    i = 0
    while i < length(tokens)
        i += 1
        τ  = tokens[i]
        τ.name == :MD_DEF_OPEN || continue
        # look ahead stopping with the first LINE_RETURN
        j = findfirst(τc -> τc.name ∈ (:LINE_RETURN, :EOS), tokens[i+1:end])
        if j !== nothing
            e = i+j
            # discard if it's a single line
            any(τx -> τx.name == :LR_INDENT, tokens[i:e]) || continue
            push!(from_to_list, (i => e))
        end
    end
    remove = Int[]
    s = str(first(tokens))
    for (i, j) in from_to_list
        si = nextind(s, to(tokens[i]))
        sj = prevind(s, from(tokens[j]))
        sc = subs(s, si, sj)
        ex, pos = Meta.parse(sc, 1)

        # find where's the first character after the effective end of the
        # definition (right-strip)
        str_expr = subs(s, si, si + prevind(sc, pos))
        stripped = strip(str_expr)
        start_id = findfirst(c -> c == stripped[1], str_expr)
        last_id  = prevind(str_expr, start_id + length(stripped))

        # find the first token after si + next_id - 1 and discard what's before
        c = findfirst(k -> from(tokens[k]) > si + last_id - 1, i+1:j)
        # we discard everything between i+1 and i+c-1 (we want to keep i+c)
        if !isnothing(c)
            append!(remove, i+1:i+c-1)
        end
    end
    deleteat!(tokens, remove)
    return nothing
end


text:
for matrix view
using Gadfly
M = rand(10,10)
spy(M)

max_stars_repo_name: xwghua/DRiP-public

max_issues_repo_name: xwghua/DRiP-public

max_forks_repo_name: xwghua/DRiP-public