text:
"""
    resp_gen(
        pars::DataFrames.DataFrame,
        theta::Vector{Float64},
        design::Matrix{Float64};
        model = "2PL", #1PL, 2PL, 3PL
        parametrization = "at-b", #"at-b, at-ab, at+b, at+ab"
        D = 1,#any number
        method = "classic uniform",
    )

# Description

Generate dichotomous responses from a dataframe of item parameters and a vector of ability points (theta).

# Arguments

- **`pars::DataFrames.DataFrame`** : Required. DataFrame with item parameters (difficulty: b or d, discrimination: a or a1, guessing: c or g).
- **`theta::Vector{Float64}`** : Required. Vector of ability points.
- **`design::Matrix{Float64}`** : Required. The design matrix.
- **`model`**: Optional. Default
- **`parametrization`** : Optional. Default: \"at-b\". Values:  \"at-b\", \"at-ab\", \"at+b\", \"at+ab\". IRT model parametrization. Ex: at-b is ``Da(\\theta-b)``.
- **`D`**: Optional. Default: 1.
- **`method`**: Optional. Default: \"classic uniform\".

# Output

- **`i::Matrix{Float64}`**: Matrix of binary responses.
"""
function resp_gen(
    pars::DataFrames.DataFrame,
    theta::Vector{Float64},
    design::Matrix{Float64};
    model = "2PL", #1PL, 2PL, 3PL
    parametrization = "at-b", #"at-b, at-ab, at+b, at+ab"
    D = 1,#any number (1, 1.6)
    method = "classic uniform",
)

    a, a2, b, c = _desume_pars(pars; model = model, parametrization = parametrization)

    n_items = size(b, 1)
    N = size(theta, 1)
    n_index = Array{Array{Int64,1},1}(undef, n_items)
    i_index = Array{Array{Int64,1},1}(undef, N)
    for n = 1:N
        i_index[n] = findall(design[:, n] .== 1.0)
        if n <= n_items
            n_index[n] = findall(design[n, :] .== 1.0)
        end
    end
    p = Matrix(undef, n_items, N)
    lp = Matrix(undef, n_items, N)
    lq = Matrix(undef, n_items, N)
    for i = 1:n_items
        for n = 1:N
            pr = c[i] + ((1 - c[i]) * (1 / (1 + exp(-(b[i] + a[i] * theta[n])))))
            p[i, n] = pr
            lp[i, n] = log(pr)
            lq[i, n] = log(1 - pr)
        end
    end
    resp = Matrix(undef, n_items, N)
    if method == "cumulated students"
        @fastmath @inbounds for i = 1:n_items
            gapScore = 3
            while gapScore >= 2
                p2 = p[i, :]
                p2 = hcat((1 .- p2), p2)#2
                unif = rand(Uniform(0, 1), N)#5
                n = 1#6
                while n <= N#7
                    csum = p2[n, 1]#8
                    cat = 0#9
                    if design[i, n] == 0
                        resp[i, n] = 0
                        n = n + 1
                    else
                        while csum < unif[n]
                            cat = cat + 1
                            if (cat == 2)
                                break
                            end
                            csum = csum + p2[n, (cat+1)]
                        end
                        resp[i, n] = cat
                        n = n + 1
                    end
                end
                gapScore = abs(sum(resp[i, :]) - sum(p[i, n] for n in n_index[i]))
            end
        end
    end
    if method == "cumulated items"
        @fastmath @inbounds for n = 1:N#4
            #gapScore=3
            #while gapScore>=2
            p2 = p[:, n]
            p2 = hcat((1 .- p2), p2)#2
            unif = rand(Uniform(0, 1), n_items)#5
            StatsBase.samplei =
                StatsBase.sample(collect(1:n_items), n_items, replace = false)
            i = 1#6
            while i <= n_items#7
                csum = p2[StatsBase.samplei[i], 1]#8
                cat = 0#9
                if design[StatsBase.samplei[i], n] == 0
                    resp[StatsBase.samplei[i], n] = 0#missing
                    i = i + 1
                else
                    while csum < unif[StatsBase.samplei[i]]
                        cat = cat + 1
                        if (cat == 2)
                            break
                        end
                        csum = csum + p2[StatsBase.samplei[i], cat+1]
                    end
                    resp[StatsBase.samplei[i], n] = cat
                    i = i + 1
                end
            end
            #gapScore=abs(sum(skipmissing(resp[:,n]))-sum(p[i,n] for i in i_index[n]))
            #if gapScore>=0.5
            #	println("person ",n," gap=",gapScore)
            #end
            #end
        end
    end
    if method == "classic uniform"
        for n = 1:N#4
            #gapScore=2
            #while gapScore>=0.5
            unif = rand(Uniform(0, 1), n_items)#5
            StatsBase.samplei =
                StatsBase.sample(collect(1:n_items), n_items, replace = false)
            i = 1#6
            while i <= n_items#7
                if design[StatsBase.samplei[i], n] == 0
                    resp[StatsBase.samplei[i], n] = 0#missing
                    i = i + 1
                else
                    if unif[StatsBase.samplei[i]] < p[StatsBase.samplei[i], n]
                        resp[StatsBase.samplei[i], n] = 1
                    else
                        resp[StatsBase.samplei[i], n] = 0
                    end
                    i = i + 1
                end
            end
            #gapScore=abs(sum(skipmissing(resp[:,n]))-sum(p[i,n] for i in i_index[n]))
            #if gapScore>=0.5
            #	println("person ",n," gap=",gapScore)
            #end
            #end
        end
    end
    if method == "cumulated pattern"
        for n = 1:N#4+
            println("start person ", n)
            i_index = i_index[n]
            i_index = size(i_index, 1)
            p2 = p[i_index, n]
            p2 = hcat((1 .- p2), p2)#2
            patterns = Vector{Vector{Int64}}(undef, 1)
            nPattern = Vector{Int64}(undef, 1)
            for r = 1:1000
                respn = Vector{Int64}(undef, i_index)
                unif = rand(Uniform(0, 1), i_index)#5
                StatsBase.samplei =
                    StatsBase.sample(collect(1:i_index), i_index, replace = false)
                for i = 1:i_index
                    csum = p2[StatsBase.samplei[i], 1]#8
                    cat = 0#9
                    while csum < unif[StatsBase.samplei[i]]
                        cat = cat + 1
                        if (cat == 2)
                            break
                        end
                        csum = csum + p2[StatsBase.samplei[i], cat+1]
                    end
                    respn[i] = cat
                end
                if r == 1
                    patterns[1] = respn
                    nPattern[1] = 1
                else
                    println(size(patterns, 1))
                    corr = 0
                    for pat = 1:size(patterns, 1)
                        if patterns[pat] == respn
                            corr = pat
                        end
                    end
                    if corr == 0
                        push!(patterns, respn)
                        nPattern = hcat(nPattern, 1)
                    else
                        nPattern[corr] = nPattern[corr] + 1
                    end
                end

            end
            nPatterns = size(patterns, 1)
            probPatterns = nPattern ./ 1000
            println(probPatterns)
            resp[i_index, n] .= StatsBase.sample(patterns, pweights(probPatterns), 1)
            println("end person ", n)
        end
    end
    if method == "classic uniform pattern"
        for n = 1:N#4+
            println("start person ", n)
            i_index = i_index[n]
            i_index = size(i_index, 1)
            patterns = Vector{Vector{Int64}}(undef, 1)
            nPattern = Vector{Int64}(undef, 1)
            for r = 1:1000
                respn = Vector{Int64}(undef, i_index)
                unif = rand(Uniform(0, 1), i_index)#5
                StatsBase.samplei =
                    StatsBase.sample(collect(1:i_index), i_index, replace = false)
                for i = 1:i_index#7
                    if unif[StatsBase.samplei[i]] < p[StatsBase.samplei[i], n]
                        respn[i] = 1
                    else
                        respn[i] = 0
                    end
                end
                if r == 1
                    patterns[1] = respn
                    nPattern[1] = 1
                else
                    println(size(patterns, 1))
                    corr = 0
                    for pat = 1:size(patterns, 1)
                        if patterns[pat] == respn
                            corr = pat
                        end
                    end
                    if corr == 0
                        push!(patterns, respn)
                        nPattern = hcat(nPattern, 1)
                    else
                        nPattern[corr] = nPattern[corr] + 1
                    end
                end
            end
            nPatterns = size(patterns, 1)
            probPatterns = nPattern ./ 1000
            resp[i_index, n] .= StatsBase.sample(patterns, pweights(probPatterns), 1)
            println("end person ", n)
        end
    end
    return resp
end

function sum_pos(x::Vector{Float64})
    x = x[x.>0]
    if size(x, 1) == 0
        return zero(Float64)::Float64
    else
        return sum(x)::Float64
    end
end

function read_settings_file(settings_file::String)
    include(settings_file)
    return Inputs::InputSettings
end


text:
################################################################################
#
#     Hecke.jl : Hecke main file
#
# This file is part of Hecke.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# (C) 2015-2019 Claus Fieker, Tommy Hofmann
# (C) 2020-2021 Claus Fieker, Tommy Hofmann, Carlo Sircana
#
################################################################################

@doc Markdown.doc"""
Hecke is a Julia package for algorithmic algebraic number theory.
For more information please visit

    https://github.com/thofma/Hecke.jl

"""
module Hecke

################################################################################
#
#  Import/export
#
################################################################################

import Base: show, minimum, rand, prod, copy, rand, ceil, round, size, in,
             powermod, ^, getindex, ==, <, >, +, *, /, \, -, !=, getindex,
             setindex!, transpose, getindex, //, floor, max, BigFloat,
             precision, first, StepRange, show, inv, div, divrem, one, zero, iseven, isodd,
             convert, angle, abs2, isless, exponent, isfinite, zeros, rem, min,
             numerator, denominator, exp, maximum, intersect, reduce, sqrt, haskey, merge,
	     powermod

# To make all exported Nemo functions visible to someone using "using Hecke"
# we have to export everything again
# dong it the "import" route, we can pick & choose...

using Requires

using LinearAlgebra, Markdown, InteractiveUtils, Libdl, Distributed, Printf, SparseArrays, Serialization, Random, Pkg, Test

import AbstractAlgebra

import LinearAlgebra: dot, istriu, nullspace, rank, ishermitian

import SparseArrays: nnz

import Serialization: serialize, deserialize

import Random: rand!
using Random: Sampler, SamplerTrivial, GLOBAL_RNG

using RandomExtensions: RandomExtensions, make, Make2, Make3, Make4

import Nemo

# TODO: remove/simplify the following once Nemo has IntegerUnion
# (and the version adding IntegerUnion is required in Project.toml)
if isdefined(Nemo, :IntegerUnion)
  import Nemo.IntegerUnion
else
  const IntegerUnion = Union{Integer, Nemo.fmpz}
end

import Pkg

exclude = [:Nemo, :AbstractAlgebra, :RealField, :zz, :qq, :factor, :call,
           :factors, :parseint, :strongequal, :window, :xgcd, :rows, :cols,
           :can_solve, :set_entry!]

for i in names(Nemo)
  i in exclude && continue
  eval(Meta.parse("import Nemo." * string(i)))
  eval(Expr(:export, i))
end

import Nemo: acb_struct, Ring, Group, Field, NmodRing, nmod, arf_struct,
             elem_to_mat_row!, elem_from_mat_row, gfp_elem, gfp_mat,
             gfp_fmpz_elem, Zmodn_poly, Zmodn_mat, GaloisField,
             GaloisFmpzField, acb_vec, array, acb_vec_clear, force_coerce,
             force_op, fmpz_mod_ctx_struct, divisors

export show, StepRange, domain, codomain, image, preimage, modord, resultant,
       next_prime, ispower, number_field, factor

###############################################################################
#
#   Library initialisation
#
###############################################################################

const pkgdir = joinpath(dirname(pathof(Hecke)), "..")

global const number_field = NumberField

function MaximalOrder
end

global const maximal_order = MaximalOrder

function __init__()
  # Because of serialization/deserialization problems, the base rings would differ otherwise.
  Hecke.Globals.Zx.base_ring = FlintZZ
  Hecke.Globals.Qx.base_ring = FlintQQ

  # Check if were are non-interactive
  bt = Base.process_backtrace(Base.backtrace())
  isinteractive_manual = all(sf -> sf[1].func != :_tryrequire_from_serialized, bt)

  # Respect the -q flag
  isquiet = Bool(Base.JLOptions().quiet)

  show_banner = !isquiet && isinteractive_manual && isinteractive() &&
                !any(x->x.name in ["Oscar"], keys(Base.package_locks)) &&
                get(ENV, "HECKE_PRINT_BANNER", "true") != "false"

  if show_banner
    println("")
    print("Welcome to \n")
    printstyled("
    _    _           _
   | |  | |         | |
   | |__| | ___  ___| | _____
   |  __  |/ _ \\/ __| |/ / _ \\
   | |  | |  __/ (__|   <  __/
   |_|  |_|\\___|\\___|_|\\_\\___|
    ", color = :red)

    println()
    print("Version")
    printstyled(" $VERSION_NUMBER ", color = :green)
    print("... \n ... which comes with absolutely no warranty whatsoever")
    println()
    println("(c) 2015-2021 by Claus Fieker, Tommy Hofmann and Carlo Sircana")
    println()
  end

  #if inNotebook()  # to make toggle work in IJulia
  #  display("text/html", "\$\\require{action}\$")
  #end

  t = create_accessors(AnticNumberField, acb_root_ctx, get_handle())
  global _get_nf_conjugate_data_arb = t[1]
  global _set_nf_conjugate_data_arb = t[2]

  t = create_accessors(AnticNumberField, Dict{Int, acb_roots}, get_handle())
  global _get_nf_conjugate_data_arb_roots = t[1]
  global _set_nf_conjugate_data_arb_roots = t[2]

  t = create_accessors(AnticNumberField,
                       Tuple{Int, nf_elem},
                       get_handle())
  global _get_nf_torsion_units = t[1]
  global _set_nf_torsion_units = t[2]

  t = create_accessors(AnticNumberField, NfOrd, get_handle())

  global _get_maximal_order_of_nf = t[1]
  global _set_maximal_order_of_nf = t[2]

  t = create_accessors(NfOrd, ClassGrpCtx, get_handle())

  global _get_ClassGrpCtx_of_order = t[1]
  global _set_ClassGrpCtx_of_order = t[2]

  t = create_accessors(NfOrd, UnitGrpCtx, get_handle())

  global _get_UnitGrpCtx_of_order = t[1]
  global _set_UnitGrpCtx_of_order = t[2]

  t = create_accessors(AnticNumberField, Dict, get_handle())

  global _get_places_uniformizers = t[1]
  global _set_places_uniformizers = t[2]

  t = create_accessors(AnticNumberField, roots_ctx, get_handle())

  global _get_roots_ctx_of_nf = t[1]
  global _set_roots_ctx_of_nf = t[2]

  t = create_accessors(AnticNumberField, Array, get_handle())

  global _get_cyclotomic_ext_nf = t[1]
  global _set_cyclotomic_ext_nf = t[2]

  t = create_accessors(AnticNumberField, Array, get_handle())

  global _get_automorphisms_nf = t[1]
  global _set_automorphisms_nf = t[2]

  t = create_accessors(AnticNumberField, NfAbsOrd{AnticNumberField, nf_elem}, get_handle())
  global _get_equation_order_of_nf = t[1]
  global _set_equation_order_of_nf = t[2]

  t = create_accessors(SimpleNumField, NfRelOrd, get_handle())

  global _get_maximal_order_of_nf_rel = t[1]
  global _set_maximal_order_of_nf_rel = t[2]

  t = create_accessors(NfOrd, MapClassGrp, get_handle())

  global _get_picard_group = t[1]
  global _set_picard_group = t[2]

  t = create_accessors(NfOrd, MapUnitGrp, get_handle())

  global _get_unit_group_non_maximal = t[1]
  global _set_unit_group_non_maximal = t[2]

  t = create_accessors(AnticNumberField, FacElemMon{AnticNumberField}, get_handle())

  global _get_fac_elem_mon_of_nf = t[1]
  global _set_fac_elem_mon_of_nf = t[2]

  t = create_accessors(NfRel, Array, get_handle())

  global _get_automorphisms_nf_rel = t[1]
  global _set_automorphisms_nf_rel = t[2]

  t = Hecke.create_accessors(AnticNumberField, Dict{Int, Tuple{qAdicRootCtx, Dict{nf_elem, Any}}}, get_handle())
  global _get_nf_conjugate_data_qAdic = t[1]
  global _set_nf_conjugate_data_qAdic = t[2]

  t = Hecke.create_accessors(AnticNumberField, Tuple{Int, Int}, get_handle())
  global _get_nf_signature = t[1]
  global _set_nf_signature = t[2]

  t = Hecke.create_accessors(AnticNumberField, Any, get_handle())
  global _get_nf_prime_data_lifting = t[1]
  global _set_nf_prime_data_lifting = t[2]

  t = Hecke.create_accessors(AnticNumberField, Any, get_handle())
  global _get_nf_norm_relation = t[1]
  global _set_nf_norm_relation = t[2]

  global R = _RealRing()

  global flint_rand_ctx = flint_rand_state()

  @require GAP="c863536a-3901-11e9-33e7-d5cd0df7b904" begin
    include("FieldFactory/fields.jl")
    include("ModAlgAss/GAPMeatAxe.jl")
    #@require Revise="295af30f-e4ad-537b-8983-00126c2a3abe" begin
    #  import .Revise
    #  #Revise.track(Hecke, joinpath(pkgdir, "src/FieldFactory/fields.jl"))
    #  #Revise.track(Hecke, "FieldFactory/abelian_layer.jl")
    #  #Revise.track(Hecke, "FieldFactory/brauer.jl")
    #  #Revise.track(Hecke, "FieldFactory/merge.jl")
    #  #Revise.track(Hecke, "FieldFactory/read_write.jl")
    #end
  end

  @require Polymake="d720cf60-89b5-51f5-aff5-213f193123e7" begin
    include("AlgAssRelOrd/NEQ_polymake.jl")
  end

  resize!(_RealRings, Threads.nthreads())
  for i in 1:Threads.nthreads()
     _RealRings[i] = _RealRing()
  end

end

module Globals
  using Hecke
  const Qx, _ = PolynomialRing(FlintQQ, "x", cached = false)
  const Zx, _ = PolynomialRing(FlintZZ, "x", cached = false)
end

using .Globals

################################################################################
#
#  Verbose printing and custom assertions
#
################################################################################

include("Assertions.jl")

################################################################################
#
#  Deprecations
#
################################################################################

include("Deprecations.jl")

################################################################################
#
#  Setter and getter for objects
#
################################################################################

function _get_maximal_order(K::AnticNumberField)
  return _get_maximal_order_of_nf(K)
end

function _set_maximal_order(K::AnticNumberField, O)
  _set_maximal_order_of_nf(K, O)
end

function _get_nf_equation_order(K::AnticNumberField)
  return _get_equation_order_of_nf(K)::NfAbsOrd{AnticNumberField, nf_elem}
end

function _set_nf_equation_order(K::AnticNumberField, O)
  _set_equation_order_of_nf(K, O)
  return nothing
end

function conjugate_data_arb(K::AnticNumberField)
  try
    c = _get_nf_conjugate_data_arb(K)::acb_root_ctx
    return c
  catch
    c = acb_root_ctx(K.pol)
    _set_nf_conjugate_data_arb(K, c)
    return c::acb_root_ctx
  end
end

function conjugate_data_arb_roots(K::AnticNumberField, p::Int)
  already_set = false
  local c
  try
    c = _get_nf_conjugate_data_arb_roots(K)::Dict{Int, acb_roots}
    already_set = true
  catch
    c = Dict{Int, acb_roots}()
  end

  if already_set && haskey(c, p)
    return c[p]
  end

  #if p > 2^18
  #  Base.show_backtr(STDOUT, backtr())
  #end
  if Nemo.iscyclo_type(K)
    # Use that e^(i phi) = cos(phi) + i sin(phi)
    # Call sincospi to determine these values
    f = get_attribute(K, :cyclo)::Int
    pstart = max(p, 2) # Sometimes this gets called with -1
    local _rall::Vector{Tuple{arb, arb}}
    rreal = arb[]
    rcomplex = Vector{acb}(undef, div(degree(K), 2))
    while true
      R = ArbField(pstart, cached = false)
      # We need to pair them
      _rall = Tuple{arb, arb}[ sincospi(fmpq(2*k, f), R) for k in 1:f if gcd(f, k) == 1]
      if all(x -> radiuslttwopower(x[1], -p) && radiuslttwopower(x[2], -p), _rall)
        CC = AcbField(pstart, cached = false)
        rall = acb[ CC(l[2], l[1]) for l in _rall]
        j = 1
        good = true
        for i in 1:degree(K)
          if ispositive(_rall[i][1])
            rcomplex[j] = rall[i]
            j += 1
          else
            if !isnegative(_rall[i][1])
              # The precision was not large enough to determine the sign of the imaginary part
              good = false
            end
          end
        end
        good && break
      end
      pstart = Int(ceil(1.3 * pstart))
    end

  else
    rootc = conjugate_data_arb(K)
    q = rootc.prec
    while q < p
      refine(rootc)
      q = rootc.prec
    end
    @assert p <= q
    rall = deepcopy(rootc.roots)
    rreal = deepcopy(rootc.real_roots)
    rcomplex = deepcopy(rootc.complex_roots)
  end

  for z in rall
    expand!(z, -p)
  end
  for z in rreal
    expand!(z, -p)
  end
  for z in rcomplex
    expand!(z, -p)
  end
  c[p] = acb_roots(p, rall, rreal, rcomplex)
  if !already_set
    _set_nf_conjugate_data_arb_roots(K, c)
  end
  return c[p]::acb_roots
end

function signature(K::AnticNumberField)
  try
    sig = _get_nf_signature(K)::Tuple{Int, Int}
    return sig
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
  end
  sig = signature(defining_polynomial(K))
  _set_nf_signature(K, sig)
  return sig::Tuple{Int, Int}
end

function _get_prime_data_lifting(K::AnticNumberField)
  try
    c = _get_nf_prime_data_lifting(K)
    return c
  catch
    c = Dict()
    _set_nf_prime_data_lifting(K, c)
    return c
  end
end

################################################################################
#
#  Intermediate backwards compatibility
#
################################################################################

trace(x...) = tr(x...)

################################################################################
#
#  Version number
#
################################################################################

deps = Pkg.dependencies()
if haskey(deps, Base.UUID("3e1990a7-5d81-5526-99ce-9ba3ff248f21"))
  ver = Pkg.dependencies()[Base.UUID("3e1990a7-5d81-5526-99ce-9ba3ff248f21")]
  if occursin("/dev/", ver.source)
    global VERSION_NUMBER = VersionNumber("$(ver.version)-dev")
  else
    global VERSION_NUMBER = VersionNumber("$(ver.version)")
  end
else
  global VERSION_NUMBER = "building"
end

######################################################################
# named printing support
######################################################################

# to use:
# in HeckeMap
#   in the show function, start with @show_name(io, map)
# for other objetcs
#   add @attributes to the struct
#   add @show_name(io, obj) to show
#   optionally, add @show_special(io, obj) as well
# on creation, or whenever, call set_name!(obj, string)
# @show_name will set on printing if bound in the REPL
# moved into AbstractAlgebra

#maps are different - as are number fields

################################################################################
#
#  Abstract map type
#
################################################################################

abstract type HeckeMap <: SetMap end  #needed here for the hasspecial stuff
             #maybe move to Maps?

import AbstractAlgebra: get_attribute, set_attribute!, @show_name, @show_special,
       _get_attributes, _get_attributes!, _is_attribute_storing_type,
       @show_special_elem, @attributes, extra_name, set_name!, find_name,
       @declare_other # remove this once it is removed

import AbstractAlgebra: get_special, set_special # TODO: remove this once Oscar doesn't need it anynmore

# Hecke maps store attributes in the header object
_get_attributes(G::Map{<:Any, <:Any, HeckeMap, <:Any}) = _get_attributes(G.header)
_get_attributes!(G::Map{<:Any, <:Any, HeckeMap, <:Any}) = _get_attributes!(G.header)
_is_attribute_storing_type(::Type{Map{<:Any, <:Any, HeckeMap, <:Any}}) = true

import Nemo: libflint, libantic, libarb  #to be able to reference libraries by full path
                                         #to avoid calling the "wrong" copy

################################################################################
#
#  Custom test function
#
################################################################################

function _adjust_path(x::String)
  if Sys.iswindows()
    return replace(x, "/" => "\\")
  else
    return x
  end
end

function test_module(x, new::Bool = true; long::Bool = false)
   julia_exe = Base.julia_cmd()
   # On Windows, we also allow bla/blub"
   x = _adjust_path(x)
   if x == "all"
     test_file = joinpath(pkgdir, "test", "runtests.jl")
   else
     test_file = joinpath(pkgdir, "test", "$x.jl")
   end

   setup_file = joinpath(pkgdir, "test", "setup.jl")

   if new
     cmd = "using Test; using Hecke; Hecke.assertions(true); long_test = $long; include(\"$(setup_file)\"); include(\"$test_file\");"
     @info("spawning ", `$julia_exe -e \"$cmd\"`)
     proj = Base.active_project()
     run(`$(julia_exe) --project=$(proj) -e $(cmd)`)
   else
     long_test = long
     assertions(true)
     @info("Running tests for $x in same session")
     include(test_file)
     assertions(false)
   end
end

################################################################################
#
#  Verbose time printing
#
################################################################################

macro vtime(args...)
  if length(args) == 2
    msg = string(args[2])
    quote
      if get_verbose_level($(args[1])) >= 1
        local t0 = time_ns()
        local val = $(esc(args[2]))
        println((time_ns()-t0)/1e9, " @ ", $msg)
        val
      else
        local val2 = $(esc(args[2]))
        val2
      end
    end
  elseif length(args) == 3
    msg = string(args[3])
    quote
      if get_verbose_level($(args[1])) >= $(args[2])
        local t0 = time_ns()
        local val = $(esc(args[3]))
        println((time_ns()-t0)/1e9, " @ ", $msg)
        val
      else
        local val2 = $(esc(args[3]))
        val2
      end
    end
  end
end

#usage
# @vtime_add_ellapsed :ClassGroup 2 clg :saturate  s= hnf(a)
# @vtime_add :ClassGroup 2 clg :saturate  0.5
# -> clg.time[:saturate] +=
function _vtime_add(D::Dict, k::Any, v::Any)
  if haskey(D, k)
    D[k] += v
  else
    D[k] = v
  end
end

macro vtime_add(flag, level, var, key, value)
  quote
    if get_verbose_level($flag) >= $level
      _vtime_add($(esc(var)).time, $key, $(esc(value)))
    end
  end
end

macro vtime_add_elapsed(flag, level, var, key, stmt)
  quote
    tm = @elapsed $(esc(stmt))
    if get_verbose_level($flag) >= $level
      _vtime_add($(esc(var)).time, $key, tm)
    end
  end
end

################################################################################
#
#  Exception types
#
################################################################################

mutable struct LowPrecisionCholesky <: Exception end

Base.showerror(io::IO, e::LowPrecisionCholesky) =
    print(io, """
    Negative diagonal in Cholesky decomposition, probably a precision issue""")

mutable struct LowPrecisionLLL <: Exception end

Base.showerror(io::IO, e::LowPrecisionLLL) =
    print(io, """
    Transformation matrix has too large entries relative to precision in LLL""")

mutable struct NotImplemented <: Exception end

Base.showerror(io::IO, ::NotImplemented) =
    print(io, """Not implemented (yet).""")

# what is this function doing here?
function checkbounds(a::Int, b::Int) nothing; end;

################################################################################
add_assert_scope(:PID_Test)
set_assert_level(:PID_Test, 0)

################################################################################
#
#  "Submodules"
#
################################################################################

include("HeckeTypes.jl")
include("NumField/NfRel/Types.jl")
include("AlgAss/Types.jl")
include("AlgAssAbsOrd/Types.jl")
include("AlgAssRelOrd/Types.jl")
include("Grp.jl")
include("Map.jl")
include("GrpAb.jl")
include("Misc.jl")
include("LinearAlgebra.jl")
include("NumField.jl")
include("NumFieldOrd.jl")
include("FunField.jl")
include("Sparse.jl")
include("BigComplex.jl")
include("conjugates.jl")
include("analytic.jl")
include("helper.jl")
include("EllCrv.jl")
include("LargeField.jl")
include("RCF.jl")
include("ModAlgAss.jl")
include("AlgAss.jl")
include("AlgAssAbsOrd.jl")
include("AlgAssRelOrd.jl")
include("LocalField.jl")
include("QuadForm.jl")
include("FieldFactory.jl")
include("../examples/NFDB.jl")

const _RealRings = _RealRing[_RealRing()]

################################################################################
#
#  Download other stuff
#
################################################################################

include("Get.jl")

################################################################################
#
#  Precompilation
#
################################################################################

#precompile(maximal_order, (AnticNumberField, ))

################################################################################
#
#  Object overloading for map types
#
################################################################################

for T in subtypes(Map(HeckeMap))
  (M::T)(a) = image(M, a)
end

################################################################################
#
#  Verbose printing
#
################################################################################

@doc Markdown.doc"""
    vshow(A) -> Nothing

Prints all fields of $A$.
"""
function vshow(A)
  for i in fieldnames(typeof(A))
    if isdefined(A, i)
      print("$i: ")
      println(getfield(A, i), "\n")
    else
      println("$i: Not definied")
    end
  end
end

################################################################################
#
#  Element types for parent types
#
################################################################################

# Nemo only provides element_types for parent objects

elem_type(::Type{FacElemMon{T}}) where {T} = FacElem{elem_type(T), T}

elem_type(::Type{Generic.ResRing{T}}) where {T} = Generic.Res{T}

################################################################################
#
#  Aliases
#
################################################################################

hasroot(a...) = ispower(a...)  # catch all... needs revisiting:
                               #hasroot(poly) != ispower(poly)....

Base.issubset(K::NumField, L::NumField) = issubfield(K, L)[1]
Base.issubset(C::ClassField, B::ClassField) = issubfield(C, B)

################################################################################
#
#  Trace function calls
#
################################################################################

_trace_cache = Dict()

function _trace_call(;cache = _trace_cache, depth = 5, print = false)
  s = Base.stacktrace()[3:3 + depth - 1]
  if !haskey(cache, s)
    if print
      println("\n Trace call ... \n ")
      Base.show_backtrace(stdout, s)
    end
    cache[s] = true
  end
end

function _print_traces(;cache = _trace_cache)
  println(cache)
end

################################################################################
#
# examples
#
################################################################################

function example(s::String)
  Base.include(Main, joinpath(dirname(pathof(Hecke)), "..", "examples", s))
end

function revise(s::String)
  s = joinpath(dirname(pathof(Hecke)), "..", "examples", s)
  Main.Revise.track(Main, s)
end

function system(s::String)
  Base.include(Main, joinpath(dirname(pathof(Hecke)), "..", "system", s))
end

function build()
  system("Build.jl")
end

function doc_init(;path=mktempdir())
  global docsproject = path
  if !isfile(joinpath(docsproject,"Project.toml"))
    cp(joinpath(pkgdir, "docs", "Project.toml"), joinpath(docsproject,"Project.toml"))
  end
  Pkg.activate(docsproject) do
    # we dev all packages with the paths from where they are currently loaded
    Pkg.develop(path=pkgdir)
    Pkg.instantiate()
    Base.include(Main, joinpath(pkgdir, "docs", "Build.jl"))
  end
end

#function doc_update_deps()
#  Pkg.activate(Pkg.update, joinpath(oscardir, "docs"))
#end

function open_doc()
    filename = normpath(pkgdir, "docs", "build", "index.html")
    @static if Sys.isapple()
        run(`open $(filename)`; wait = false)
    elseif Sys.islinux() || Sys.isbsd()
        run(`xdg-open $(filename)`; wait = false)
    elseif Sys.iswindows()
        cmd = get(ENV, "COMSPEC", "cmd.exe")
        run(`$(cmd) /c start $(filename)`; wait = false)
    else
        @warn("Opening files the default application is not supported on this OS.",
              KERNEL = Sys.KERNEL)
    end
end

function build_doc(; doctest=false, strict=false, format=:mkdocs)
  if !isdefined(Main, :Build)
    doc_init()
  end
  Pkg.activate(docsproject) do
    Base.invokelatest(Main.Build.make, Hecke; strict=strict, local_build=true, doctest=doctest, format=format)
  end
  if format == :html
    open_doc()
  elseif format == :mkdocs
    println("""Run `mkdocs serve` inside `../Hecke/docs/` to view the documentation.

            Use `format = :html` for a simplified version of the docs which does
            not require `mkdocs`.
            """)
  end
end

#html_build = Ref(false)
#
#function build_doc(html::Bool = false)
#  _html_build = html_build[]
#  html_build[] = html
#  Base.include(Main, joinpath(dirname(pathof(Hecke)), "..", "docs", "make_local.jl"))
#  html_build[] = _html_build
#end

function percent_P()
  s = Base.active_repl.mistate
  REPL = Base.REPL_MODULE_REF.x

  #from Rafael:
  function print_history(hist)
    println()
    for i = hist.start_idx+1:lastindex(hist.history)
      if hist.modes[i] == :julia
        println('[', i-hist.start_idx, "] ", hist.history[i])
      end
    end
  end
  print_history(REPL.LineEdit.mode(s).hist)
end

export percent_P

#same (copied) as in stdlib/v1.0/InteractiveUtils/src/InteractiveUtils.jl
#difference: names(m, all = true) to also see non-exported variables, aka
# caches...

function varinfo(m::Module=Main, pattern::Regex=r"")
    rows =
        Any[ let value = getfield(m, v)
                 Any[string(v),
                     (value===Base || value===Main || value===Core ? "" : format_bytes(summarysize(value))),
                     summary(value)]
             end
             for v in sort!(names(m, all = true)) if isdefined(m, v) && occursin(pattern, string(v)) ]

    pushfirst!(rows, Any["name", "size", "summary"])

    return Markdown.MD(Any[Markdown.Table(rows, Symbol[:l, :r, :l])])
end
varinfo(pat::Regex) = varinfo(Main, pat)

function print_cache(sym::Vector{Any})
  for f in sym;
    #if f[2] isa Array || f[2] isa Dict || f[2] isa IdDict;
    try
      print(f[1], " ", length(f[2]), "\n");
    catch e
    end
  end
end

function print_cache()
  print_cache(find_cache(Nemo))
  print_cache(find_cache(Nemo.Generic))
  print_cache(find_cache(Hecke))
end

function find_cache(M::Module)
  sym = []
  for a in collect(names(M, all = true))
    d = Meta.parse("$M.$a")
      try
        z = eval(d);
        if isa(z, AbstractArray) || isa(z, AbstractDict)
          push!(sym, (d, z))
        end
    catch e
    end
  end
  return sym
end

protect = [:(Hecke.ASSERT_LOOKUP), :(Hecke.VERBOSE_LOOKUP),
           :(Hecke.ASSERT_SCOPE), :(Hecke.VERBOSE_SCOPE),
           :(Hecke._euler_phi_inverse_maximum),
           :(Hecke.odlyzko_bound_grh),
           :(Hecke.nC), :(Hecke.B1), #part of ECM
           :(Hecke.VERBOSE_PRINT_INDENT),
           :(Hecke._RealRings),
           :(Hecke.protect)] # We need to protect protect itself :)
                             # Otherwise it might emptied and then everything
                             # is emptied.

function clear_cache(sym::Vector{Any})
  for f in sym;
    if f[1] in protect
      continue
    end
    try
      empty!(f[2])
    catch e
    end
  end
end

function clear_cache()
  clear_cache(find_cache(Nemo))
  clear_cache(find_cache(Nemo.Generic))
  clear_cache(find_cache(Hecke))
end

@inline __get_rounding_mode() = Base.MPFR.rounding_raw(BigFloat)

using .NormRel

end # module

ext: jl


text:
using Random, Plots, StatsPlots, Distributions, DataFrames, Printf
# Define calculation of posterior
"""
Posterior Distribution (Gamma dist) for Poisson parameter λ

    Inputs
    ------
    data: data (0 or 1)
    α0: shape parameter of prior (Gamma dist)
    β0: inversed scale parameter of prior (Gamma dist)
    prob: credible interval probability (0 < prob < 1)

    Outputs
    -------
    results: posterior summary stats table (DataFrames)
    α_star: shape parameter of posterior (Gamma dist)
    β: inversed scale parameter (1 / θ) of posterior (Gamma dist)

"""
function poisson_stats(data, α0, β0, prob)
    # posterior
    n = size(data)[1]
    sum_data = sum(data)
    α_star = sum_data + α0
    β_star = n + β0
    θ_star = 1.0 / β_star
    d_posterior = Gamma(α_star, θ_star)
    # mean
    mean_λ = mean(d_posterior)
    # median
    median_λ = median(d_posterior)
    # mode
    mode_λ = mode(d_posterior)
    # Analytic Mode for Beta Dist
    # (a - 1.0) / (a + b - 2.0)
    # std
    sd_λ = std(d_posterior)
    # credible interval
    # ci upper boud
    ci_λ_upper = quantile(d_posterior, prob - (1-prob)/2)
    # ci lower boud
    ci_λ_lower = quantile(d_posterior, (1-prob)/2)
    ci_λ = [ci_λ_lower, ci_λ_upper]

    posterior_stats = reshape([mean_λ, median_λ, mode_λ, sd_λ, ci_λ], 1, 5)
    stats_string = ["Mean", "Median", "Mode", "SD", "CI"]
    results = DataFrame(posterior_stats, stats_string)
    return results, α_star, β_star
end

# Generate moc datas
λ = 3.0
n = 50
Random.seed!(99)
data = rand(Poisson(lam), n)

# Set prior
α0 = 1.0
β0 = 1.0

# Calculate posterior statistics
prob = 0.95
results, α_star, β_star = poisson_stats(data, α0, β0, prob)
print(results)

# Visualize Posterior
θ0, θ_star = 1.0/β0, 1.0/β_star
plot(Gamma(α0, θ0), xlims=(0, 6), ylims=(0, 1.75), xlabel="λ", ylabel="Probability Density", label=@sprintf("prior: Gamma(%.1f, %.1f)",α0, θ0), linestyle=:dash)
plot!(Gamma(α_star, θ_star), label=@sprintf("prior: Gamma(%.1f, %.2f)",α_star, θ_star), linestyle=:solid)


text:
import ..Datasets: dataset, datafile, reader, Mode, Test

rte_init() = register(DataDep(
    "GLUE-RTE",
    """
    Recognizing Textual Entailment (RTE) task (GLUE version)
    """,
    "https://dl.fbaipublicfiles.com/glue/data/RTE.zip",
    "6bf86de103ecd335f3441bd43574d23fef87ecc695977a63b82d5efb206556ee";
    post_fetch_method = fn -> begin
      mv(fn, "RTE.zip")
      DataDeps.unpack("RTE.zip")
      innerdir = "RTE"
      innerfiles = readdir(innerdir)
      mv.(joinpath.(innerdir, innerfiles), innerfiles)
      rm(innerdir)
    end
))

struct RTE <: Dataset end

function dataset(::Type{M}, d::RTE) where M <: Mode
    ds = reader(datafile(M, d))
    header = split(take!(ds), '\t')
    field_num = length(header)
    needed_field = M == Test ? (2,3) : (2,3,4)
    rds = get_channels(String, length(needed_field))

    task = @async begin
        for s ∈ ds
            s = split(s, '\t')
            length(s) != field_num && continue
            for (i, j) ∈ enumerate(needed_field)
                put!(rds[i], s[j])
            end
        end
    end
    for rd ∈ rds
        bind(rd, task)
    end

    rds
end

function trainfile(::RTE)
    datadep"GLUE-RTE/train.tsv"
end

function devfile(::RTE)
    datadep"GLUE-RTE/dev.tsv"
end

function testfile(::RTE)
    datadep"GLUE-RTE/test.tsv"
end

get_labels(::RTE) = ("entailment", "not_entailment")


text:
@testset "$TableType" for TableType in (columntable, rowtable, DataFrame)

    table = TableType((a=[1, 2, 3], b=[4, 5, 6]))

    @test is_transformable(table)

    @testset "apply" begin

        @testset "OneToOne" begin

            T = FakeOneToOneTransform()

            @test isequal(
                FeatureTransforms.apply(table, T),
                TableType((Column1=ones(3), Column2=ones(3)))
            )
            @test isequal(
                FeatureTransforms.apply(table, T; cols=:a),
                TableType((Column1=ones(3), ))
            )
            @test isequal(
                FeatureTransforms.apply(table, T; header=[:x, :y]),
                TableType((x=ones(3), y=ones(3)))
            )
        end

        @testset "OneToMany" begin

            T = FakeOneToManyTransform()

            @test isequal(
                FeatureTransforms.apply(table, T),
                TableType((Column1=ones(3), Column2=ones(3), Column3=ones(3), Column4=ones(3)))
            )
            @test isequal(
                FeatureTransforms.apply(table, T; cols=:a),
                TableType((Column1=ones(3), Column2=ones(3)))
            )
            @test isequal(
                FeatureTransforms.apply(table, T; header=[:x, :y, :p, :q]),
                TableType((x=ones(3), y=ones(3), p=ones(3), q=ones(3)))
            )
        end

        @testset "ManyToOne" begin

            T = FakeManyToOneTransform()

            @test isequal(
                FeatureTransforms.apply(table, T),
                TableType((Column1=ones(3),))
            )
            @test isequal(
                FeatureTransforms.apply(table, T; cols=(:a, :b)),
                TableType((Column1=ones(3),))
            )
            @test isequal(
                FeatureTransforms.apply(table, T; header=[:x]),
                TableType((x=ones(3),))
            )
        end

        @testset "ManyToMany" begin

            T = FakeManyToManyTransform()

            @test isequal(
                FeatureTransforms.apply(table, T),
                TableType((
                    Column1=ones(3), Column2=ones(3), Column3=ones(3),  Column4=ones(3)
                ))
            )
            @test isequal(
                FeatureTransforms.apply(table, T; cols=(:a, :b)),
                TableType((
                        Column1=ones(3), Column2=ones(3), Column3=ones(3),  Column4=ones(3)
                    ))
                )

            @test isequal(
                FeatureTransforms.apply(table, T; header=[:p, :q, :r, :s]),
                TableType((p=ones(3), q=ones(3), r=ones(3), s=ones(3)))
            )
        end

    end

    # TODO: test passing in dims / inds too apply! after addressing
    # https://github.com/invenia/FeatureTransforms.jl/issues/68
    @testset "apply!" begin
        T = FakeOneToOneTransform()

        # Cannot mutate NamedTuple elements
        if isrowtable(table)
            @test_throws MethodError FeatureTransforms.apply!(table, T)
        else
            _table = deepcopy(table)
            FeatureTransforms.apply!(_table, T)
            @test _table == TableType((a=ones(3), b=ones(3)))

            _table = deepcopy(table)
            FeatureTransforms.apply!(_table, T; cols=:a)
            @test _table == TableType((a=[1, 1, 1], b=[4, 5, 6]))

            _table = deepcopy(table)
            @test_broken FeatureTransforms.apply!(_table, T; cols=:b, dims=[1, 2])
            @test_broken _table == TableType((a=[1, 2, 3], b=[1, 1, 6]))
        end
    end

    @testset "apply_append" begin

        @testset "OneToOne" begin

            T = FakeOneToOneTransform()

            @test isequal(
                FeatureTransforms.apply_append(table, T),
                TableType((a=[1, 2, 3], b=[4, 5, 6], Column1=ones(3), Column2=ones(3)))
            )

            @test isequal(
                FeatureTransforms.apply_append(table, T; header=[:x, :y]),
                TableType((a=[1, 2, 3], b=[4, 5, 6], x=ones(3), y=ones(3)))
            )
        end

        @testset "OneToMany" begin

            T = FakeOneToManyTransform()

            @test isequal(
                FeatureTransforms.apply_append(table, T),
                TableType((
                    a=[1, 2, 3], b=[4, 5, 6],
                    Column1=ones(3), Column2=ones(3), Column3=ones(3), Column4=ones(3)
                ))
            )
            @test isequal(
                FeatureTransforms.apply_append(table, T; cols=:a),
                TableType((a=[1, 2, 3], b=[4, 5, 6], Column1=ones(3), Column2=ones(3)))
            )
            @test isequal(
                FeatureTransforms.apply_append(table, T; header=[:x, :y, :p, :q]),
                TableType((
                    a=[1, 2, 3], b=[4, 5, 6], x=ones(3), y=ones(3), p=ones(3), q=ones(3)
                ))
            )
        end

        @testset "ManyToOne" begin

            T = FakeManyToOneTransform()

            @test isequal(
                FeatureTransforms.apply_append(table, T),
                TableType((a=[1, 2, 3], b=[4, 5, 6], Column1=ones(3)))
            )

            @test isequal(
                FeatureTransforms.apply_append(table, T; header=[:x]),
                TableType((a=[1, 2, 3], b=[4, 5, 6], x=ones(3)))
            )
        end

        @testset "ManyToMany" begin

            T = FakeManyToManyTransform()

            @test isequal(
                FeatureTransforms.apply_append(table, T),
                TableType((
                    a=[1, 2, 3], b=[4, 5, 6],
                    Column1=ones(3), Column2=ones(3),
                    Column3=ones(3), Column4=ones(3),
                ))
            )

            @test isequal(
                FeatureTransforms.apply_append(table, T; header=[:p, :q, :r, :s]),
                TableType((
                    a=[1, 2, 3], b=[4, 5, 6], p=ones(3), q=ones(3), r=ones(3), s=ones(3)
                ))
            )
        end

    end

end


text:
# Code used to generate closest correlation matrix problems. Save in Testproblems/ClosestCorr with an accurante
# solution from MOSEK

# Original problem has the following format:
# min_X   1/2 ||X-C||^2
# s.t.    Xii = 1
#         X ⪴ 0

# workspace()
# using JLD, JuMP, Mosek, Distributions

# specify number of problems that start with matrix C as a perturbed correlation matrix and problems where C is a random matrix
nnCorr = 25
nnRand = 25
nn = nnCorr+nnRand
# this function creates a matrix A that slices out the diagonal entries Xii of a vectorized square matrix x=vec(X)
function createDiagonalExtractor(n)
  A = spzeros(n,n^2)
  A[1,1] = 1
  for iii=2:n-1
    col = (iii-1)*(n+1)
    A[iii,col+1] = 1
  end
  A[n,n^2] = 1
  return A
end

# https://stats.stackexchange.com/questions/2746/how-to-efficiently-generate-random-positive-semidefinite-correlation-matrices
# VINE METHOD to generate random correlation matrices
#distributed ~ det(S)^eta [or maybe det(S)^(eta-1), not sure]
function randomCorrMatrix(d, eta)
    β = eta + (d-1)/2;
    P = zeros(d,d);           #storing partial correlations
    S = eye(d);

    for k = 1:d-1
        β = β - 1/2;
        for i = k+1:d
            P[k,i] = rand(Beta(β,β)) # sampling from β
            P[k,i] = (P[k,i]-0.5)*2;     #linearly shifting to [-1, 1]
            p = P[k,i];
            for l = (k-1):-1:1 #converting partial correlation to raw correlation
                p = p * sqrt((1-P[l,i]^2)*(1-P[l,k]^2)) + P[l,i]*P[l,k];
            end
            S[k,i] = p;
            S[i,k] = p;
        end
    end
    return S
end

rng = MersenneTwister(12345)
dirPath = "../DataFiles/Julia/ClosestCorr/"
!ispath(dirPath) && mkdir(dirPath)

 for iii =1:1:nn
  # choose size of problem
  n = rand(rng,6:40)
  eta = abs(randn(rng))
  # create random correlation matrix and perturb it
  if iii<=nnCorr
    C = randomCorrMatrix(n,eta)
    E  = randn(n,n)*1e-1
    C = C + (E + E')/2
  else
    C = randn(rng,n,n)
  end

  c = vec(C)

  isposdef(C) && warn("The perturbed correlation matrix is still pos def.")

  # put problem into solver format
  # min   x'Ix - 1*t
  # s.t.  Ax + s1 = 0
  #       x = s4
  #       -t <= c'*x <=t
  # s1 in 0, s2,s3 in R+, s4 in SDPcone

  n2 = n^2
  m = n+n2

  P = speye(n2)
  q = -vec(C)
  r = 0.5*vec(C)'*vec(C)
  b = [ones(n);zeros(n2)]
  A = createDiagonalExtractor(n)
  Aa = [A; -speye(n2)]
  # specify cone
  Kf = n
  Kl = 0
  Kq = []
  Ks = [n^2]
  # # solve accurately once with mosek
  model = Model(solver=MosekSolver())
  @variable(model, X[1:n,1:n], SDP)
  @variable(model, t)
  x = vec(X)
  @objective(model, Min, t)
  @constraint(model, norm(x-c) <= t)
  @constraint(model, A*x.== b[1:n])
  status = JuMP.solve(model)

  # correct optimal objective value since slightly different problem is solved
  objTrue = getobjectivevalue(model)
  objCorrected = 0.5*objTrue^2
  solTrue = getvalue(x)
  # println("Objective value: ",objTrue)
  # println("Corrected objective value: ", objCorrected)
  # println("x = ", getvalue(x))

  nr = "$(iii)"
  if iii < 10
    nr = "0$(iii)"
  end

  fn = "ClosestCorr"*nr*".jld"
  problemType = "Closest Correlation Matrix"
  problemName = "ClosestCorr"*nr

  JLD.save(dirPath*fn,"n",n,"m",m,"A",Aa,"b",b,"P",P,"q",q,"r",r,"C",C,"objTrue",objCorrected,"solTrue",solTrue,"problemType",problemType,"problemName",problemName,"Kf",Kf,"Kl",Kl,"Kq",Kq,"Ks",Ks)
  println("$(iii)/$(nn) completed!")
end


text:
#include("example_nnet.jl")
using LinearAlgebra
function ComputeValid(predict_func::Function, W, xvalid, yvalid, valid)
    yhat = predict_func(W,[ones(size(xvalid,1)) xvalid],nHidden)
    push!(valid,sum((yhat-yvalid).^2)/length(yvalid))

end

function BS_Select(;validation_array,f_array, gradcalc::Function,xtrain,ytrain,W,nHidden, alist, Blist)
    #Get most recent step and batch sizes
    a=alist[end]
    B=Blist[end]
    n,d=size(xtrain)

    #Get Gradient Angle
    i=rand(1:n,B)
    j=rand(1:n,B)
    _,gi=gradcalc(W, xtrain[i,:],ytrain[i],nHidden)
    _,gj=gradcalc(W,xtrain[j,:],ytrain[j], nHidden)
    costheta=dot(gi,gj)/(norm(gi)*norm(gj))
    #theta=acos(costheta)

    #alpha decreases
    if length(validation_array) >= 4 && validation_array[end-2]>validation_array[end]
        a_adjust = (1 - abs(validation_array[end-2] - validation_array[end])/validation_array[end-2])^(1/4)
        a = a*a_adjust
    end
    push!(alist,a)
    push!(Blist,B)
    print(a," ")

    return a,B
end

function SGDBabysitter(gradcalc::Function, nn_predict::Function, maxIter, nHidden,
                        nParams, xtrain, ytrain, xvalid, yvalid )

    n = size(xtrain,1)
    valid=[]
    flist=[]
    alist=[]
    Blist=[]
    #placehold for initialization
    B=100
    a = 1e-4
    push!(Blist,B)
    push!(alist,a)
    W= randn(nParams,1)

    for t in 1:maxIter

    	i = rand(1:n,B)
    	f,g = gradcalc(W, xtrain[i,:], ytrain[i], nHidden)
    	W = W - a*g
        #Every few iterations, compute validation and append f,g,validation

          if (mod(t-1,round(maxIter/100)) == 0)
              # print("Training iteration = $(t-1) ")
              push!(flist,f)
              ComputeValid(nn_predict, W, xvalid, yvalid, valid)
          end
        if (mod(t-1,round(maxIter/50)) == 0)
            a,B=BS_Select(validation_array=valid,f_array=flist, gradcalc=gradcalc,xtrain=xtrain,ytrain=ytrain,W=W,nHidden=nHidden, alist=alist, Blist=Blist)

          end
    end

    display(plot(1:length(valid), valid))

    return W, valid
end


text:
using BenchmarkTools
using LazilyInitializedFields
import Random
import StableRNGs

const rng = StableRNGs.StableRNG(0)
randstring(args...) = Random.randstring(rng, args...)
rand(args...) = Random.rand(rng, args...)

const SUITE = BenchmarkGroup()
const N = 10^6

mutable struct NormalMutable
    x::Float64
    y::Float64
end
const normal_mutables = [NormalMutable(rand(Float64), rand(Float64)) for i in 1:N]

@lazy mutable struct Mutable
    x::Float64
    @lazy y::Float64
end
const mutables = [Mutable(rand(Float64), rand(Float64)) for i in 1:N]

struct NormalImmutable
    x::Float64
    y::Float64
end
const normal_immutables = [NormalImmutable(rand(Float64), rand(Float64)) for i in 1:N]

@lazy struct Immutable
    x::Float64
    @lazy y::Float64
end
const immutables = [Immutable(rand(Float64), rand(Float64)) for i in 1:N]

@lazy struct ImmutableType{T}
    x::T
    @lazy y::Float64
end
const immutables_type = [ImmutableType(rand(Float64), rand(Float64)) for i in 1:N]

function perf_sum_eager(X)
    s = 0
    @simd for v in X
        s += v.x
    end
    return s
end

function perf_sum_lazy(X)
    s = 0.0
    for v in X
        s += v.y
    end
    return s
end

SUITE["normal mutable - eager"]   = @benchmarkable perf_sum_eager($normal_mutables)
SUITE["mutable - eager"]          = @benchmarkable perf_sum_eager($mutables)
SUITE["normal immutable - eager"] = @benchmarkable perf_sum_eager($normal_immutables)
SUITE["immutable - eager"]        = @benchmarkable perf_sum_eager($immutables)
SUITE["immutable_type - eager"]   = @benchmarkable perf_sum_eager($immutables_type)

SUITE["normal mutable - lazy"]   = @benchmarkable perf_sum_lazy($normal_mutables)
SUITE["mutable - lazy"]          = @benchmarkable perf_sum_lazy($mutables)
SUITE["normal immutable - lazy"] = @benchmarkable perf_sum_lazy($normal_immutables)
SUITE["immutable - lazy"]        = @benchmarkable perf_sum_lazy($immutables)
SUITE["immutable_type - lazy"]   = @benchmarkable perf_sum_lazy($immutables_type)


text:
using ImageInterpLast

using Test

a = zeros(10,10,10,10)
for t = 1:10
    a[:,:,:,t] .= t
end

img0 = interp_last(a, 0.5; cache3d=false, out_type=Float64)
img1 = interp_last(a, 0.5; cache3d=false, out_type=Float32)
img2 = ImgItpLast(a, 0.5, 2, Float64)
img3 = ImgItpLast(a, 0.5, 3, Float64)

for img in (img0,img1,img2,img3)
    @test all(img[:,:,:,1].==1.5)
    @test all(img[1:2,:,:,1].==1.5)
    @test all(img[1:2,:,2,1].==1.5)
    @test all(img[:,:,:,2].==2.5)
    @test all(img[1:2,:,:,2].==2.5)
    @test all(img[1:2,:,2,2].==2.5)
end

coefs = fill(0.25, 10)
coefs[6:10] .= 1.0

img2 = ImgItpLast(a, coefs, 2, Float64)
img3 = ImgItpLast(a, coefs, 3, Float64)

for img in (img2,img3)
    @test all(img[:,:,1:5,1].==1.75)
    @test all(img[1:2,:,1:5,1].==1.75)
    @test all(img[1:2,:,2,1].==1.75)

    @test all(img[:,:,6:10,1].==1.0)
    @test all(img[1:2,:,6:10,1].==1.0)
    @test all(img[1:2,:,7,1].==1.0)

    @test all(img[:,:,1:5,2].==2.75)
    @test all(img[1:2,:,1:5,2].==2.75)
    @test all(img[1:2,:,2,2].==2.75)

    @test all(img[:,:,6:10,2].==2.0)
    @test all(img[1:2,:,6:10,2].==2.0)
    @test all(img[1:2,:,7,2].==2.0)
end

img4 = interp_last(a, coefs; cache3d=false, out_type=Float64)
img5 = interp_last(a, coefs; cache3d=true, out_type=Float64)

@test all(img4[:,:,:,:].==img2[:,:,:,:])
@test all(img5[:,:,:,:].==img3[:,:,:,:])


text:
using SchattenNorms
using Test
using LinearAlgebra

import Random
import LinearAlgebra

Random.seed!(31415926)

@testset "snorm sanity checks" begin
    @test isapprox(snorm([1 0; 0 -1]),sqrt(2.0))
    @test isapprox(snorm([1 0; 0 -1],1.0),2.0)
    @test isapprox(snorm([1 0; 0 -1],Inf),1.0)
end

@testset "snorm" begin
    for i=1:100
        for d=2:4
            R = randn(d,d)
            U = LinearAlgebra.svd(randn(d,d)).U
            V = LinearAlgebra.svd(randn(d,d)).U
            p = sort(abs.(randn(3)) .+ 1.0)
            @test isapprox(snorm(R),snorm(R,2.0))
            @test isapprox(snorm(U*R*V,1),snorm(R,1))
            @test isapprox(snorm(U*R*V),snorm(R,2.0))
            @test isapprox(snorm(U*R*V,Inf),snorm(R,Inf))
            @test isapprox(snorm(U*R*V,p[1]),snorm(R,p[1]))
            @test snorm(R,1) >= snorm(R,2) >= snorm(R,Inf)
            @test snorm(R,p[1]) >= snorm(R,p[2]) >= snorm(R,p[3])
            @test snorm(R,1.0) == snorm(R,1) == trnorm(R) == nucnorm(R)
        end
    end
end

# E
# U,V
# @test isapprox(cbnorm(E),dnorm(E'))
# @test isapprox(cbnorm(E),cbnorm(U*E*V))
# @test isapprox(dnorm(E),dnorm(U*E*V))

include("test-dnorm.jl")


text:
require("../src/quantumoptics.jl")

using quantumoptics

srand(0)

N1 = 2
N2 = 2
N3 = 2
N4 = 2
N5 = 4
N6 = 4
N7 = 32
N = N1*N2*N3*N4*N5*N6*N7

b1 = GenericBasis([N1])
b2 = GenericBasis([N2])
b3 = GenericBasis([N3])
b4 = GenericBasis([N4])
b5 = GenericBasis([N5])
b6 = GenericBasis([N6])
b7 = GenericBasis([N7])

basis_total = CompositeBasis(b1, b2, b3, b4, b5, b6); N = N1*N2*N3*N4*N5*N6
#basis_total = CompositeBasis(b1, b2, b3, b4, b5, b6, b7)

op1 = Operator(b1, rand(N1,N1))
op2 = Operator(b2, rand(N2,N2))

psi = Ket(basis_total, rand(N))

psi2 = Ket(basis_total, rand(N))
psi = tensor(psi, dagger(psi))
@time psi2 = tensor(psi2, dagger(psi2))

l1 = operators_lazy.LazyTensor(basis_total, basis_total, [1], [op1])
l2 = operators_lazy.LazyTensor(basis_total, basis_total, [2], [op2])
l3 = operators_lazy.LazyTensor(basis_total, basis_total, [1,2], [op1,op2])

#@time operators_lazy.gemm2!(Complex(1.), l3, psi.data, Complex(0.), psi2.data)
#@time operators_lazy.gemm2!(Complex(1.), l3, psi.data, Complex(0.), psi2.data)
#println(vecnorm(vec((psi2-full(l3)*psi).data)))
@time operators_lazy.gemm!(Complex(1.), l3, psi.data, Complex(0.), psi2.data)
@time operators_lazy.gemm!(Complex(1.), l3, psi.data, Complex(0.), psi2.data)
println(vecnorm(vec((psi2-full(l3)*psi).data)))

# @time operators_lazy.gemm2!(Complex(1.), psi.data, l3, Complex(0.), psi2.data)
# @time operators_lazy.gemm2!(Complex(1.), psi.data, l3, Complex(0.), psi2.data)
# println(vecnorm(vec((psi2-psi*full(l3)).data)))
@time operators_lazy.gemm!(Complex(1.), psi.data, l3, Complex(0.), psi2.data)
@time operators_lazy.gemm!(Complex(1.), psi.data, l3, Complex(0.), psi2.data)
println(vecnorm(vec((psi2-psi*full(l3)).data)))

# @time operators_lazy.gemm!(Complex(1.), (l2+l3), psi.data, Complex(0.), psi2.data)
# @time operators_lazy.gemm!(Complex(1.), (l2+l3), psi.data, Complex(0.), psi2.data)
# println(vecnorm(vec(((l3+l2)*psi - psi2).data)))

# @time operators_lazy.gemm!(Complex(1.), psi.data, (l2+l3), Complex(0.), psi2.data)
# @time operators_lazy.gemm!(Complex(1.), psi.data, (l2+l3), Complex(0.), psi2.data)
# println(vecnorm(vec((psi*(l3+l2) - psi2).data)))
# println(vecnorm(vec(((l3+l2)*psi - l3*psi - l2*psi).data)))


text:
module HTTP2

const Headers = Vector{Tuple{String,String}}

bytearr(a::Vector{UInt8}) = a
bytearr(cs::Base.CodeUnits{UInt8,String}) = convert(Vector{UInt8}, cs)
bytearr(s::String) = bytearr(codeunits(s))

# package code goes here
include("Frame.jl")
include("Session.jl")

end # module


text:
module TALib

using DataStructures

pkg_name = "TALib"
depsjl = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
println(depsjl)
if isfile(depsjl)
    include(depsjl)
else
    error("$pkg_name not properly installed. Please run Pkg.build(\"$pkg_name\")")
end

include("constants.jl")
include("path.jl")

include("describe.jl")
D_INDICATORS, INDICATORS, D_GROUPS = get_ta_func_constants()
export FunctionDescriptionXML

include("tools.jl")

function _ta_check_success(function_name::String, ret_code::TA_RetCode)
    if isequal(TA_SUCCESS::TA_RetCode, ret_code)
        return true
    else
        error("$function_name function failed with error code $errorCode")
    end
end

for f in [:GetVersionString, :GetVersionMajor, :GetVersionMinor, :GetVersionBuild, :GetVersionDate, :GetVersionTime]
    f_str = string(f)
    f_ta_str = "TA_" * f_str
    @eval begin
        function ($f)()
            unsafe_string(ccall(($f_ta_str, libta_lib), Cstring, ()))
        end
    end

    @eval begin
        export ($f)
    end

end

# ===

for f in [:Initialize, :Shutdown]
    f_str = string(f)
    f_ta_str = "TA_" * f_str
    @eval begin
        function ($f)()
            ta_func = () -> ccall( ($f_ta_str, libta_lib), TA_RetCode, ())
            ret_code = ta_func()
            _ta_check_success($f_str, ret_code)
        end
    end

    @eval begin
        export ($f)
    end
end

# ===

code_generators = ["ta_func_api_code_level0",
    "ta_func_api_code_level1",
    "ta_func_api_code_level2_dataframes",
    "ta_func_api_code_level2_timearrays"]

import_generated = true
for code_generator in code_generators
    filename = joinpath(basepath(), "generated", code_generator * ".jl")
    if !isfile(filename)
        import_generated = false
        @warn("code generation is required for $filename")
        break
    end
end
if !import_generated
    warn("code generation is required to use TALib.jl")
    include("ta_func_api_gen.jl")
end

for code_generator in code_generators
    filename = joinpath("generated", code_generator * ".jl")
    try
        #info("include $filename")
        include(filename)
    catch e
        error(e)
    end
end

for f in INDICATORS
    @eval begin
        export ($f)
    end
end

export D_INDICATORS, INDICATORS, D_GROUPS

end # module


text:
import Base: (^), (*), (~), (/), (|)
ParserOperatorTypes = Union{AbstractToken, AbstractString, Char}

(*)(x, y::AbstractToken) = sSequence(parser(x),y)
(*)(x::AbstractToken, y) = sSequence(x,parser(y))
"""
    (*)(x::Any, y::AbstractToken)
    (*)(x::AbstractToken, y::Any)
    (*)(x::AbstractToken, y::AbstractToken)

Chain parsers in [`sSequence`](@ref).
See also [`@seq`](@ref).
"""
(*)(x::AbstractToken, y::AbstractToken) = sSequence(x,y)

## todo: cuts

"""
    (/)(x::ParserOperatorTypes, y::ParserOperatorTypes)

`Sequence(PositiveLookbehind(x),y)`

```jldoctest
julia> match("is "/"match", "no match is match").offset
13
```

!!! note
    This syntax is reviewed and I your appreciate your comments!
"""
(/)(x::ParserOperatorTypes, y::ParserOperatorTypes) =
    Sequence(PositiveLookbehind(x),y)

"""
    Base.broadcasted(::typeof((&)), x::ValueNotIn, y::ValueNotIn)

Character matchers `m` like `Union{ValueIn,ValueNotIn,T}`, or any
type `T` providing a `ismatch(m::T,c::Char)::Bool` method represent a
"sparse" bitarray for all characters.

!!! note
    Please consider the broadcast API a draft you are invited to comment to.

```jldoctest
julia> CharNotIn("abc") .& CharNotIn("z")
[^abcz] ValueNotIn
::Char

julia> CharIn("abc") .& CharNotIn("c")
[ab] ValueIn
::Char
```
"""
Base.broadcasted(::typeof((&)), x::ValueNotIn, y::ValueNotIn) =
    ValueNotIn(x.pcre*y.pcre, x.sets, y.sets)

Base.broadcasted(::typeof((&)), x::ValueIn, y::ValueNotIn) =
    ValueIn(setdiff(x.sets, y.sets))

Base.broadcasted(::typeof((&)), x::Union{ValueIn,ValueNotIn}, ::AnyValue) =
    x

Base.broadcasted(::typeof((&)), x::NamedParser, y) =
    NamedParser(x.name,x.parser .& y; doc=x.doc)

Base.broadcasted(::typeof((&)), x::MatchedSubSequence, y) =
    MatchedSubSequence(x.parser .& y)

Base.broadcasted(::typeof((&)), x::Transformation, y) =
    Transformation(x.transform, x.parser .& y)

Base.broadcasted(::typeof((|)), x::ValueIn, y::ValueIn) =
    ValueIn(x.pcre*y.pcre, x.sets,y.sets)

(|)(x, y::ParserOperatorTypes) = Either(parser(x),y; simplify=true)
(|)(x::ParserOperatorTypes, y) = Either(x,parser(y); simplify=true)
(|)(x::ParserOperatorTypes, y::ParserOperatorTypes) = Either(x,y; simplify=true)
(|)(x::AbstractToken, y::ParserOperatorTypes) = Either(x,y; simplify=true)

"""
    (|)(x::AbstractToken, y)
    (|)(x, y::AbstractToken)
    (|)(x::AbstractToken, y::AbstractToken)

Operator syntax for `Either(x, y; simplify=true)`.

```jldoctest
julia> 'a' | CharIn("AB") | "bc"
|🗄 Either
├─ a
├─ [AB] ValueIn
└─ bc
::Union{Char, SubString{String}}

```
"""
(|)(x::CombinedParser, y::CombinedParser)

"""
    (|)(x::AbstractToken{T}, default::Union{T,Missing})

Operator syntax for `Optional(x, default=default)`.

```jldoctest
julia> parser("abc") | "nothing"
|🗄 Either
├─ abc
└─ nothing
::SubString{String}

julia> parser("abc") | missing
abc? |missing
::Union{Missing, SubString{String}}

```

"""
(|)(x::AbstractToken{T}, default::Union{T,Missing}) where { T } = Optional(x,default=default)

(|)(x::AbstractToken{Any}, y::AbstractToken) = Either(x, y; simplify=true)

function (|)(x::Char, y::Char)
    ValueIn(tuple(x,y))
end
function (|)(x::ValueIn, y::Char)
    ValueIn(tuple(x.sets...,y))
end

"""
    `(|)(x::Either, T::Type)`

Return new Either with `T` added to result_type(x).
todo: Note that the options array is kept. As a consequence `push!`on result will also push to `x`.
"""
(|)(x::Either, T::Type) =
    Either{Union{result_type(x),T}}(x.options)


text:
@testset "Scaling" begin
    data = (x=collect(1:10),
            y=rand(10) .+ 3,
            z=Symbol.(repeat('a':'e', 2)))

    @testset "Automatic scaling" begin
        xc = concrete_term(term(:x), data, Scale())
        @test xc isa ScaledTerm
        @test xc.scale == std(data.x)
        @test modelcols(xc, data) == data.x ./ std(data.x) == data.x ./ xc.scale

        yc = concrete_term(term(:y), data, Scale())
        @test yc isa ScaledTerm
        @test yc.scale == std(data.y)
        @test modelcols(yc, data) == data.y ./ std(data.y) == data.y ./ yc.scale

        @testset "alternative scale function" begin
            f = mad
            xc = concrete_term(term(:x), data, Scale(f))
            @test xc isa ScaledTerm
            @test xc.scale == f(data.x)
            @test modelcols(xc, data) == data.x ./ f(data.x) == data.x ./ xc.scale
            # why test this? well this makes sure that our tests
            # wouldn't pass if we were using the default scale function
            # in other words, this tests we're actually hitting a different codepath
            @test !isapprox(std(data.x), f(data.x))
        end
    end

    @testset "Manual scaling" begin
        xc = concrete_term(term(:x), data, Scale(5))
        @test xc isa ScaledTerm
        @test xc.scale == 5 != std(data.x)
        @test modelcols(xc, data) == data.x ./ 5 == data.x ./ xc.scale

        yc = concrete_term(term(:y), data, Scale(5))
        @test yc isa ScaledTerm
        @test yc.scale == 5 != std(data.y)
        @test modelcols(yc, data) == data.y ./ 5 == data.y ./ yc.scale
    end

    @testset "Schema hints dict" begin
        sch = schema(data, Dict(:x => Scale(), :y => Scale(2)))
        xc = sch[term(:x)]
        @test xc isa ScaledTerm
        @test !StatsModels.needs_schema(xc)
        @test xc.scale == std(data.x)

        yc = sch[term(:y)]
        @test yc isa ScaledTerm
        @test yc.scale == 2

        @test modelcols(xc, data) == data.x ./ std(data.x)
        @test modelcols(yc, data) == data.y ./ 2
    end

    @testset "scale function" begin
        sch = schema(data)
        x = sch[term(:x)]

        xc = scale(x)
        @test xc isa ScaledTerm
        @test xc.scale == sqrt(x.var)

        xc2 = scale(x, Scale(2))
        xc22 = scale(x, 2)
        @test xc2.scale == xc22.scale == 2

        @testset "scaling vectors" begin
            x = collect(1:5) # int
            y = x * 5.0      # float
            @test scale(x) == x ./ std(x)
            @test scale(y) == y ./ std(y)
            @test scale(maximum, x) == (x ./ 5)
            @test all(scale(y, x) .== 5)

            @test scale!(copy(y)) == scale(y)
            # test mutation
            y2 = copy(y)
            scale!(y2)
            @test y2 == scale(y)

            # this is int on int so you can modify in place
            @test scale!(minimum, copy(x)) == x
            # this converts exactly to int, so you can do it in place
            @test scale!([2 4 6]) == scale([2 4 6])
            @test_throws InexactError scale!(std, [1, 2])
            @test_throws InexactError scale!([1, 2])
            @test_throws MethodError scale!(v -> 1, ["a", "b"])
        end
    end

    @testset "plays nicely with formula" begin
        f = @formula(0 ~ x * y)

        sch = schema(f, data)

        sch_c = schema(f, data, Dict(:x => Scale(), :y => Scale(2)))
        ff_c = apply_schema(f, sch_c)

        mm_c = modelcols(ff_c.rhs, data)
        @test mm_c == hcat(data.x ./ std(data.x),
                           data.y ./ 2,
                           (data.x ./ std(data.x)) .* (data.y ./ 2))

        @test coefnames(ff_c.rhs) ==
              ["x(scaled: 3.03)", "y(scaled: 2)", "x(scaled: 3.03) & y(scaled: 2)"]

        # round-trip schema is empty since needs_schema is false
        sch_2 = schema(ff_c, data)
        @test isempty(sch_2.schema)

        @testset "missing omit" begin
            d = (x=[1, 2, missing, 3], y=[1, missing, 2, 3], z=[missing, 1, 2, 3])
            dd, nonmissing = StatsModels.missing_omit(d, ff_c)
            # z is ignored:
            @test findall(nonmissing) == [1, 4]
            # x and y are handled correclty
            @test dd.x == d.x[nonmissing]
            @test dd.y == d.y[nonmissing]
        end
    end

    @testset "printing" begin
        xc = concrete_term(term(:x), data, Scale())
        @test StatsModels.termsyms(xc) == Set([:x])
        @test "$(xc)" == string_mime(MIME("text/plain"), xc) == "x(scaled: 3.03)"
        @test coefnames(xc) == "x(scaled: 3.03)"
    end

    @testset "categorical term" begin
        z = concrete_term(term(:z), data)
        @test_throws ArgumentError scale(z)
        @test_throws ArgumentError scale(z, Scale())
        zc = scale(z, Scale(0.5))
        @test modelcols(zc, data) == modelcols(z, data) ./ 0.5
        @test StatsModels.width(zc) == StatsModels.width(z)
        @test coefnames(zc) == coefnames(z) .* "(scaled: 0.5)"

        zc2 = scale(z, Scale([1 2 3 4]))
        @test string_mime(MIME("text/plain"), zc2) == "z(scaled: [1 2 3 4])"
        @test modelcols(zc2, data) == modelcols(z, data) ./ [1 2 3 4]
        @test coefnames(zc2) == coefnames(z) .* "(scaled: " .* string.([1, 2, 3, 4]) .* ")"
    end
end


text:
using RoamesGeometry
using StaticArrays
using Displaz
using CoordinateTransformations
using Colors

@static if VERSION >= v"0.7"
    linspace(a,b,c) = range(a, stop = b, length = c)
end

# Make a scene

pole1 = Line(SVector(0.0, 0.0, 0.0), SVector(0.0, 0.0, 8.0))
crossarm1 = Line(SVector(-0.6, 0.0, 7.0), SVector(0.6, 0.0, 7.0))

pole2 = Line(SVector(0.0, 20.0, 0.0), SVector(0.0, 20.0, 8.0))
crossarm2 = Line(SVector(-0.6, 20.0, 7.0), SVector(0.6, 20.0, 7.0))

wire1 = Quadratic(SVector(-0.5, 0.0, 6.8), SVector(-0.5, 10.0, 5.0), SVector(-0.5, 20.0, 6.8))
wire2 = Quadratic(SVector(0.5, 0.0, 6.8), SVector(0.5, 10.0, 5.0), SVector(0.5, 20.0, 6.8))
wire3 = Quadratic(SVector(0.0, 0.0, 8.0), SVector(0.0, 10.0, 6.2), SVector(0.0, 20.0, 8.0))

catwire = Catenary(pi/2, 0.0, 10.0, 6.2, -10.0, 10.0, 28.07275088)

# Combine objects
lines = [pole1, crossarm1, pole2, crossarm2]
quads = [wire1, wire2, wire3]
cats = [catwire]

# plot them
Displaz.plot3d(lines; color = [0.2, 0.2, 0.8])
Displaz.plot3d(quads; color = [0.8, 0.2, 0.1])
Displaz.plot3d(cats; color = [0.8, 0.8, 0.1])

Displaz.plot3d(boundingbox(lines); color = [0.3, 0.3, 0.3])
#Displaz.plot3d(boundingbox(quads); color = [0.3, 0.3, 0.3])
Displaz.plot3d(boundingbox(cats); color = [0.3, 0.3, 0.3])

# Make another bay, test single plotting and transformations

tfm = Translation(SVector(0.0, 20.0, 0.0))

Displaz.plot3d(tfm(pole2), color = [0.2, 0.2, 0.8])
Displaz.plot3d(tfm(crossarm2), color = [0.2, 0.2, 0.8])
Displaz.plot3d(tfm(wire1), color = [0.8, 0.2, 0.1])
Displaz.plot3d(tfm(wire2), color = [0.8, 0.2, 0.1])
Displaz.plot3d(tfm(wire3), color = [0.8, 0.2, 0.1])
Displaz.plot3d(tfm(catwire), color = [0.8, 0.8, 0.1])

# Test the distance functions

dx = 50.0
dy = -30.0
cat = Catenary(pi/2, dx, 10.0 + dy, 6.2, -10.0, 10.0, 10.0)
#cat = Translation(SVector(dx, 0.0, 0.0))(catwire)
points = vec([SVector(x + dx, y + dy, z) for x ∈ linspace(0.0, 4.0, 20), y ∈ linspace(-4.0, 24.0, 140), z ∈ linspace(2.0, 15.0, 130)])
dists = map(p -> distance(cat, p), points)
colors = RGB.(exp.(-1 .* dists), 2 .* dists .* exp.(-1 .* dists), 3 .* (dists .* dists .- 1) .* exp.(-1 .* dists))
Displaz.plot3d(cat, color = [0.8, 0.8, 0.1])
Displaz.plot3d(points, color = colors, label = "distance field")

dx = 50.0
dy = 0.0
cat = Catenary(pi/2, dx, 10.0 + dy, 6.2, -10.0, 10.0, 10.0)
#cat = Translation(SVector(dx, 0.0, 0.0))(catwire)
points = vec([SVector(x + dx, y + dy, z) for x ∈ linspace(0.0, 4.0, 20), y ∈ linspace(-4.0, 24.0, 140), z ∈ linspace(2.0, 15.0, 130)])
dists = map(p -> powerline_distances(cat, p)[2], points)
colors = RGB.(1 ./ (1 .+ (dists .+ 1).*(dists .+ 1)), 1 ./ (1 .+ dists .* dists), 1 ./ (1 .+ (dists .+ 2) .* (dists .+ 2)))
Displaz.plot3d(cat, color = [0.8, 0.8, 0.1])
Displaz.plot3d(points, color = colors, label = "r field")

dx = 50.0
dy = 30.0
cat = Catenary(pi/2, dx, 10.0 + dy, 6.2, -10.0, 10.0, 10.0)
#cat = Translation(SVector(dx, 0.0, 0.0))(catwire)
points = vec([SVector(x + dx, y + dy, z) for x ∈ linspace(0.0, 4.0, 20), y ∈ linspace(-4.0, 24.0, 140), z ∈ linspace(2.0, 15.0, 130)])
dists = map(p -> powerline_distances(cat, p)[3], points)
colors = RGB.(1 ./ (1 .+ 0.25 .* (dists .- 2).*(dists .- 2)), 1 ./ (1 .+ 0.25 .* dists .* dists), 1 ./ (1. .+ 0.25 .* (dists .+ 2).*(dists .+ 2)))
Displaz.plot3d(cat, color = [0.8, 0.8, 0.1])
Displaz.plot3d(points, color = colors, label = "z field")

dx = 50.0
dy = 60.0
cat = Catenary(pi/2, dx, 10.0 + dy, 6.2, -10.0, 10.0, 10.0)
#cat = Translation(SVector(dx, 0.0, 0.0))(catwire)
points = vec([SVector(x + dx, y + dy, z) for x ∈ linspace(0.0, 4.0, 20), y ∈ linspace(-4.0, 24.0, 140), z ∈ linspace(2.0, 15.0, 130)])
dists = map(p -> powerline_distances(cat, p)[4], points)
colors = RGB.(1 ./ (1 .+ 0.02 .* (dists .- 10).*(dists .- 10)), 1 ./ (1 .+ 0.02 .* dists.*dists), 1 ./ (1 .+ 0.02 .* (dists .+ 10).*(dists .+ 10)))
Displaz.plot3d(cat, color = [0.8, 0.8, 0.1])
Displaz.plot3d(points, color = colors, label = "d field")


text:
# circ operator

f(x) = x + 1
g(x) = x - 1

println((f ∘ g)(1))

# Map composition

println(map(first ∘ reverse ∘ uppercase, split("compose Functions Using map")))

# Function chaining by piping

1:10 |> sum |> x -> x^2 |> println

# Dot Syntax

arr = [1, 24, 36]

operatedArr = sqrt.(arr)
println(operatedArr)

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
using SortMerge, Random, Test

# --------------------------------------------------------------------
a1 = [1,2,4,6,7,10]
a2 = [2,3,5,9,10]
j = sortmerge(a1, a2)
@test size(j)[1] == 2
@test j[1][1] == 2
@test j[1][2] == 6
@test j[2][1] == 1
@test j[2][2] == 5
@test a1[countmatch(j, 1) .== 0]  == [1, 4, 6, 7]
@test a2[countmatch(j, 2) .== 0]  == [3, 5, 9]

# --------------------------------------------------------------------
nn = 1_0_000
a1 = rand(MersenneTwister(0), 1.:nn,   nn);
a2 = rand(MersenneTwister(1), 1.:nn, 2*nn);
a1 = unique(a1);

j = sortmerge(a1, a2)
@test sum(abs.(a1[j[1]] .- a2[j[2]])) == 0
cm = countmatch(j, 1); for i in 1:length(cm); @test cm[i] == length(findall(a1[i] .== a2)); end
cm = countmatch(j, 2); for i in 1:length(cm); @test cm[i] == length(findall(a2[i] .== a1)); end

println(a1[countmatch(j, 1) .== 8])
kk = multimatch(j, 1, 8)

println(a1[kk[1]])
println(a2[kk[2]])
@test sum(abs.(a1[kk[1]] .- a2[kk[2]])) == 0

k1 = findall(countmatch(j, 1) .== 0);
k2 = findall(countmatch(j, 2) .== 0);
for i in k1; @test length(findall(a1[i] .== a2)) == 0; end
for i in k2; @test length(findall(a2[i] .== a1)) == 0; end

sort!(a1);
sort!(a2);
j = sortmerge(a1, a2, sorted=true)

@test sum(abs.(a1[j[1]] .- a2[j[2]])) == 0
cm = countmatch(j, 1); for i in 1:length(cm); @test cm[i] == length(findall(a1[i] .== a2)); end
cm = countmatch(j, 2); for i in 1:length(cm); @test cm[i] == length(findall(a2[i] .== a1)); end


text:
#wave = cwts[5]; bc = bcs[3]; β=βs[1]; ave = averagingLengths[1]; n = ns[1]
@testset "Wavelet properties" begin
    bcs =(ContinuousWavelets.PerBoundary(),ContinuousWavelets.SymBoundary(), ContinuousWavelets.ZPBoundary())
    cwts = (cDb2,ContinuousWavelets.dog1,ContinuousWavelets.dog2,ContinuousWavelets.paul1,ContinuousWavelets.morl,ContinuousWavelets.morse)
    βs =(1,4)
    ns = (128, 1382,2039)
    averagingLengths=(0,.5,1,2,4)
    @testset "length $n, with type $wave, bc $bc, β=$β, ave=$(ave)" for n in ns, wave in cwts, bc in bcs, β in βs, ave in averagingLengths
        wav = ContinuousWavelets.wavelet(wave, β=β,boundary=bc,averagingLength = ave)
        Ŵ,ω = (3,1)
        with_logger(ConsoleLogger(stderr, Logging.Error)) do
            Ŵ,ω = ContinuousWavelets.computeWavelets(n, wav)
        end
        if wav.waveType isa Union{<:Morlet, <:Paul, <:Morse}
            nFreq, nSpace = ContinuousWavelets.setn(n, wav)
            fullŴ = [Ŵ; zeros(nFreq-1, size(Ŵ, 2))]
            fullŴ[end-nFreq+1:end,1] = reverse(Ŵ[:,1])
            spaceWaves = ifft(fullŴ,1)
        else
            if bc isa PerBoundary
                nLong = n
            elseif bc isa SymBoundary
                nLong = 2n
            elseif bc isa ZPBoundary
                nLong = 2^(ceil(Int,log2(n)))
            end
            spaceWaves = irfft(Ŵ,nLong,1)
        end
        # since the signals are centered at the boundary, the area that should be zero is in the middle
        nonSupported = spaceWaves[ceil(Int,n/2):end-ceil(Int,n/2),:]
        supported = spaceWaves[[1:floor(Int,n/2); end-floor(Int,n/2):end],:]
        # make sure the wavelets don't have significant support out into the mirrored signal
        if !(bc isa PerBoundary) && !(wave isa Paul) && size(Ŵ,2) > 1
            @test max([norm(nonSupported[:,i],Inf) / norm(supported[:,i],Inf) for i = 2:size(spaceWaves,2)]...) ≤ 1e-2
        elseif wave isa Paul && size(Ŵ,2) > 1 # the Paul wavelet decay is just too slow to avoid bleedover, so best to ignore the averaging size
            @test max([norm(nonSupported[:,i],Inf) / norm(supported[:,i],Inf) for i = 2:size(spaceWaves,2)]...) ≤ 1e-1
        end
        # if the averaging length is 0 and its only averaging at the sizes given, something is wrong
        @test size(Ŵ,2)>1 || ave>0
        if !(wave isa ContOrtho) && size(Ŵ,2) > 1
            # Guaranteed in a very different way for ContOrtho, and no guarantees made for just averaging
            @test max(abs.(Ŵ[end, :])...)  / norm(supported,Inf) ≤ 1e-1
        end
    end
end


text:
function nterms(t)
    if istree(t)
        return reduce(+, map(nterms, arguments(t)), init=0)
    else
        return 1
    end
end

# Soft pivoted
# Note: we call this function with a matrix of Union{SymbolicUtils.Symbolic, Any}
function sym_lu(A; check=true)
    SINGULAR = typemax(Int)
    m, n = size(A)
    F = map(x->x isa Num ? x : Num(x), A)
    minmn = min(m, n)
    p = Vector{LinearAlgebra.BlasInt}(undef, minmn)
    info = 0
    for k = 1:minmn
        kp = k
        amin = SINGULAR
        for i in k:m
            absi = _iszero(F[i, k]) ? SINGULAR : nterms(F[i,k])
            if absi < amin
                kp = i
                amin = absi
            end
        end

        p[k] = kp

        if amin == SINGULAR && !(amin isa Symbolic) && (amin isa Number) && iszero(info)
            info = k
        end

        # swap
        for j in 1:n
            F[k, j], F[kp, j] = F[kp, j], F[k, j]
        end

        for i in k+1:m
            F[i, k] = F[i, k] / F[k, k]
        end
        for j = k+1:n
            for i in k+1:m
                F[i, j] = F[i, j] - F[i, k] * F[k, j]
            end
        end
    end
    check && LinearAlgebra.checknonsingular(info, Val{true}())
    LU(F, p, convert(LinearAlgebra.BlasInt, info))
end

# Given a vector of equations and a
# list of independent variables,
# return the coefficient matrix `A` and a
# vector of constants (possibly symbolic) `b` such that
# A \ b will solve the equations for the vars
function A_b(eqs::AbstractArray, vars::AbstractArray, check)
    exprs = rhss(eqs) .- lhss(eqs)
    if check
        for ex in exprs
            @assert isaffine(ex, vars)
        end
    end
    A = jacobian(exprs, vars)
    b = A * vars - exprs
    A, b
end

"""
$(TYPEDSIGNATURES)

Solve equation(s) `eqs` for a set of variables `vars`.

Assumes `length(eqs) == length(vars)`

Currently only works if all equations are linear. `check` if the expr is linear
w.r.t `vars`.
"""
function solve_for(eq, var; simplify=false, check=true) # scalar case
    # simplify defaults for `false` as canonicalization should handle most of
    # the cases.
    a, b, islinear = linear_expansion(eq, var)
    check && @assert islinear
    islinear || return nothing
    # a * x + b = 0
    x = a \ -b
    simplify ? SymbolicUtils.simplify(x, expand=true) : x
end

function solve_for(eqs::AbstractArray, vars::AbstractArray; simplify=true, check=true)
    length(eqs) == 1 == length(vars) && return [solve_for(eqs[1], vars[1]; simplify=simplify, check=check)]
    A, b = A_b(eqs, vars, check)
    #TODO: we need to make sure that `solve_for(eqs, vars)` contains no `vars`
    sol = _solve(A, b, simplify)
    map(Num, sol)
end

function _solve(A::AbstractMatrix, b::AbstractArray, do_simplify)
    A = SymbolicUtils.simplify.(Num.(A), expand=true)
    b = SymbolicUtils.simplify.(Num.(b), expand=true)
    sol = value.(sym_lu(A) \ b)
    do_simplify ? SymbolicUtils.simplify.(sol, expand=true) : sol
end

LinearAlgebra.ldiv!(A::UpperTriangular{<:Union{Symbolic,Num}}, b::AbstractVector{<:Union{Symbolic,Num}}, x::AbstractVector{<:Union{Symbolic,Num}} = b) = symsub!(A, b, x)
function symsub!(A::UpperTriangular, b::AbstractVector, x::AbstractVector = b)
    LinearAlgebra.require_one_based_indexing(A, b, x)
    n = size(A, 2)
    if !(n == length(b) == length(x))
        throw(DimensionMismatch("second dimension of left hand side A, $n, length of output x, $(length(x)), and length of right hand side b, $(length(b)), must be equal"))
    end
    @inbounds for j in n:-1:1
        _iszero(A.data[j,j]) && throw(SingularException(j))
        xj = x[j] = b[j] / A.data[j,j]
        for i in j-1:-1:1
            sub = _isone(xj) ? A.data[i,j] : A.data[i,j] * xj
            if !_iszero(sub)
                b[i] -= sub
            end
        end
    end
    x
end

LinearAlgebra.ldiv!(A::UnitLowerTriangular{<:Union{Symbolic,Num}}, b::AbstractVector{<:Union{Symbolic,Num}}, x::AbstractVector{<:Union{Symbolic,Num}} = b) = symsub!(A, b, x)
function symsub!(A::UnitLowerTriangular, b::AbstractVector, x::AbstractVector = b)
    LinearAlgebra.require_one_based_indexing(A, b, x)
    n = size(A, 2)
    if !(n == length(b) == length(x))
        throw(DimensionMismatch("second dimension of left hand side A, $n, length of output x, $(length(x)), and length of right hand side b, $(length(b)), must be equal"))
    end
    @inbounds for j in 1:n
        xj = x[j] = b[j]
        for i in j+1:n
            sub = _isone(xj) ? A.data[i,j] : A.data[i,j] * xj
            if !_iszero(sub)
                b[i] -= sub
            end
        end
    end
    x
end

minor(B, j) = B[2:end, 1:size(B,2) .!= j]
function LinearAlgebra.det(A::AbstractMatrix{<:Num}; laplace=true)
    if laplace
        n = LinearAlgebra.checksquare(A)
        if n == 1
            return A[1, 1]
        elseif n == 2
            return A[1, 1] * A[2, 2] - A[1, 2] * A[2, 1]
        else
            return sum((-1)^(1+j) * A[1,j] * det(minor(A,j), laplace=true)
                       for j in axes(A,2))
        end
    else
        if istriu(A) || istril(A)
            return det(UpperTriangular(A))
        end
        return det(lu(A; check = false))
    end
end

function LinearAlgebra.norm(x::AbstractArray{Num}, p::Real=2)
    p = value(p)
    issym = p isa Symbolic
    if !issym && p == 2
        sqrt(sum(x->abs2(x), x))
    elseif !issym && isone(p)
        sum(abs, x)
    elseif !issym && isinf(p)
        mapreduce(abs, max, x)
    else
        sum(x->abs(x)^p, x)^inv(p)
    end
end

"""
    (a, b, islinear) = linear_expansion(t, x)

When `islinear`, return `a` and `b` such that `a * x + b == t`.
"""
function linear_expansion(t, x)
    a, b, islinear = _linear_expansion(t, x)
    x isa Num ? (Num(a), Num(b), islinear) : (a, b, islinear)
end
# _linear_expansion always returns `Symbolic`
function _linear_expansion(t::Equation, x)
    a₂, b₂, islinear = linear_expansion(t.rhs, x)
    islinear || return (a₂, b₂, false)
    a₁, b₁, islinear = linear_expansion(t.lhs, x)
    # t.rhs - t.lhs = 0
    return (a₂ - a₁, b₂ - b₁, islinear)
end
trival_linear_expansion(t, x) = isequal(t, x) ? (1, 0, true) : (0, t, true) # trival case
function _linear_expansion(t, x)
    t = value(t)
    t isa Symbolic || return (0, t, true)
    x = value(x)
    istree(t) || return trival_linear_expansion(t, x)
    isequal(t, x) && return (1, 0, true)

    op, args = operation(t), arguments(t)
    op isa Differential && trival_linear_expansion(t, x)

    if op === (+)
        a₁ = b₁ = 0
        islinear = true
        # (a₁ x + b₁) + (a₂ x + b₂) = (a₁ + a₂) x + (b₁ + b₂)
        for (i, arg) in enumerate(args)
            a₂, b₂, islinear = linear_expansion(arg, x)
            islinear || return (a₁, b₁, false)
            a₁ += a₂
            b₁ += b₂
        end
        return (a₁, b₁, true)
    elseif op === (-)
        @assert length(args) == 1
        a, b, islinear = linear_expansion(args[1], x)
        return (-a, -b, islinear)
    elseif op === (*)
        # (a₁ x + b₁) (a₂ x + b₂) = a₁ a₂ x² + (a₁ b₂ + a₂ b₁) x + b₁ b₂
        a₁ = 0
        b₁ = 1
        islinear = true
        for (i, arg) in enumerate(args)
            a₂, b₂, islinear = linear_expansion(arg, x)
            (islinear && _iszero(a₁ * a₂)) || return (a₁, b₁, false)
            a₁ = a₁ * b₂ + a₂ * b₁
            b₁ *= b₂
        end
        return (a₁, b₁, true)
    elseif op === (^)
        # (a₁ x + b₁)^(a₂ x + b₂) is linear => a₂ = 0 && (b₂ == 1 || a₁ == 0)
        a₂, b₂, islinear = linear_expansion(args[2], x)
        (islinear && _iszero(a₂)) || return (0, 0, false)
        a₁, b₁, islinear = linear_expansion(args[1], x)
        _isone(b₂) && (a₁, b₁, islinear)
        (islinear && _iszero(a₁)) || return (0, 0, false)
        return (0, b₁^b₂, islinear)
    else
        for (i, arg) in enumerate(args)
            a, b, islinear = linear_expansion(arg, x)
            (_iszero(a) && islinear) || return (0, 0, false)
        end
        return (0, t, true)
    end
end