text:
using JSON
using MacroTools

abstract type AbstractHGFConfig <: AbstractDict{Symbol, Any} end
abstract type HGFPretrainedConfig <: AbstractHGFConfig end
abstract type HGFConfig <: AbstractHGFConfig end

Base.length(cfg::AbstractHGFConfig) = fieldcount(typeof(cfg))
Base.isempty(::AbstractHGFConfig) = false
Base.getindex(cfg::AbstractHGFConfig, k::Symbol) = getproperty(cfg, k)
Base.iterate(cfg::AbstractHGFConfig) = iterate(cfg, 1)
Base.iterate(cfg::AbstractHGFConfig, i) = i <= length(cfg) ? (fieldname(typeof(cfg), i) => getfield(cfg, i), i+1) : nothing

function load_config(model_name)
  cfg = JSON.parsefile(get_registered_config_path(model_name); dicttype=Dict{Symbol, Any})
  model_type = Symbol(cfg[:model_type])
  return load_config(Val(model_type), cfg)
end

config_type(::Val{model_type}) where {model_type} = error("Unknown model type: $model_type")

function load_config(model_type::Val, cfg)
  cfg_type = config_type(model_type)
  return cfg_type(; cfg...)
end

"""
  `save_config(model_name, config; path=pwd(), config_name=DEFAULT_CONFIG_NAME)`

save the `config` at `<path>/<model_name>/<config_name>`.
"""
function save_config(model_name, config; path=pwd(), config_name=DEFAULT_CONFIG_NAME)
  model_path = joinpath(path, model_name)
  !isdir(model_path) && error("$model_path is not a dir.")
  config_file = joinpath(model_path, config_name)
  open(config_file, "w+") do io
    JSON.print(io, config)
  end
  return config_file
end

struct HGFPretrainedConfigBase <: HGFPretrainedConfig
  # encoder-decoder
  is_encoder_decoder::Bool
  is_decoder::Bool
  add_cross_attention::Bool
  tie_encoder_decoder::Bool
  # sequence generation
  max_length::Int
  min_length::Int
  do_sample::Bool
  early_stopping::Bool
  num_beams::Int
  temeperature::Float64
  top_k::Int
  top_p::Float64
  repetition_penalty::Float64
  length_penalty::Float64
  no_repeat_ngram_size::Int
  bad_words_ids::Nothing
  num_return_sequence::Int
  chunk_size_feed_forward::Int
  # fine-tune task
  architectures::Nothing
  finetuning_task::Nothing
  id2label::Nothing
  label2id::Nothing
  num_labels::Int
  # tokenizer
  prefix::Nothing
  bos_token_id::Nothing
  pad_token_id::Nothing
  eos_token_id::Nothing
  decoder_start_token_id::Nothing

  HGFPretrainedConfigBase() = new(
    false, false, false, false,
    20, 0, false, false, 1, 1.0, 50, 1.0, 1.0, 1.0, 0, nothing, 1, 0,
    nothing, nothing, nothing, nothing, 2,
    nothing, nothing, nothing, nothing, nothing,
  )
end

struct UpdatedConfig{C <: HGFPretrainedConfig, V} <: HGFPretrainedConfig
  config::C
  name::Symbol
  value::V
end

update(cfg::HGFPretrainedConfigBase, name::Symbol, v) = UpdatedConfig(cfg, name, v)
function update(cfg::UpdatedConfig, name::Symbol, v)
  if _name(cfg) == name
    return UpdatedConfig(_config(cfg), name, v)
  else
    return UpdatedConfig(
      update(_config(cfg), name, v),
      _name(cfg), _value(cfg))
  end
end

_name(cfg::UpdatedConfig) = getfield(cfg, :name)
_value(cfg::UpdatedConfig) = getfield(cfg, :value)
_config(cfg::UpdatedConfig) = getfield(cfg, :config)
_base(cfg::UpdatedConfig) = (_base ∘ _config)(cfg)
_base(cfg::HGFPretrainedConfigBase) = cfg

Base.getproperty(cfg::UpdatedConfig, k::Symbol) = k == _name(cfg) ? _value(cfg) : Base.getproperty(_config(cfg), k)

_get_cfg(cfg::UpdatedConfig, k) = k == _name(cfg) ? cfg : _get_cfg(_config(cfg), k)

_current_pair(cfg::UpdatedConfig) = _name(cfg)=>_value(cfg)

_updated(cfg::UpdatedConfig) = tuple(_current_pair(cfg), _updated(_config(cfg))...)
_updated(cfg::HGFPretrainedConfigBase) = ()

Base.iterate(cfg::UpdatedConfig) = iterate(cfg, _name(cfg))
function Base.iterate(cfg::UpdatedConfig, i)
  isnothing(i) && return nothing
  cur = _get_cfg(cfg, i)
  p = _current_pair(cur)
  nc = _config(cur)
  if nc isa UpdatedConfig
    return p, _name(nc)
  else
    return p, nothing
  end
end

function pretrained_config(; kws...)
  global PRETRAIN_CONFIG
  p = PRETRAIN_CONFIG
	for (k, v) in kws
    p = update(p, k, v)
  end
  return p
end

_parent(cfg::HGFConfig) = getfield(cfg, :_parent)

Base.length(cfg::UpdatedConfig) = length(_updated(cfg))
function Base.summary(io::IO, x::UpdatedConfig)
  n = length(x)
  print(io, typeof(x).name)
  print(io, " with ", n, (n==1 ? " entry" : " entries"))
end

const PRETRAIN_CONFIG = HGFPretrainedConfigBase()

Base.length(cfg::HGFConfig) = fieldcount(typeof(cfg)) + length(_updated(_parent(cfg))) - 1
Base.iterate(cfg::HGFConfig) = iterate(cfg, 1)
function Base.iterate(cfg::HGFConfig, i::Int)
  n = fieldcount(typeof(cfg))
  if i < n
    return fieldname(typeof(cfg), i) => getfield(cfg, i), i+1
  else
    pcfg = _parent(cfg)
    if pcfg isa UpdatedConfig
      return iterate(pcfg)
    else
      return nothing
    end
  end
end

Base.iterate(cfg::HGFConfig, i::Union{Symbol, Nothing}) = iterate(_parent(cfg), i)

macro cfgdef(ex)
  @capture(ex, struct T_ fields__ end) || error("only accept struct definition like Base.@kwdef")
  fields_wo_default = map(fields) do fwv
    @capture(fwv, f_ = v_) ? f : fwv
  end
  field_names = map(fields_wo_default) do f
    @capture(f, fn_::ft_) ? fn : f
  end
  @capture(T, Tname_ <: Sname_) || (Tname = T)
  Tname = esc(Tname)
  param_ex = map(fields) do fwv
    @capture(fwv, f_ = v_) ? Expr(:kw, @capture(f, fn_::ft_) ? fn : f, esc(v)) : fwv
  end
  kws = gensym(:kws)
  push!(param_ex, :($(kws)...))
  param_ex = Expr(:parameters, param_ex...)

  return quote
    struct $(esc(T))
      $(esc.(fields_wo_default)...)
      _parent::HGFPretrainedConfig
    end

    function $(Tname)($(param_ex))
      p = pretrained_config(; $(kws)...)
      $(Tname)($(field_names...), p)
    end
  end
end

Base.getproperty(cfg::HGFConfig, k::Symbol) = hasfield(typeof(cfg), k) ? getfield(cfg, k) : getproperty(_parent(cfg), k)

include("./config_bert.jl")
include("./config_gpt2.jl")
include("./config_roberta.jl")


text:
__precompile__()
module easyClipper

    const depsfile = joinpath(@__DIR__, "..", "deps", "deps.jl")

    if isfile(depsfile)
        include(depsfile)
    else
        error("Clipper not build correctly. Please run Pkg.build(\"Clipper\")")
    end

    function __init__()
        check_deps()
    end

    export PolyType, PolyTypeSubject, PolyTypeClip,
           ClipType, ClipTypeIntersection, ClipTypeUnion, ClipTypeDifference, ClipTypeXor,
           PolyFillType, PolyFillTypeEvenOdd, PolyFillTypeNonZero, PolyFillTypePositive, PolyFillTypeNegative,
           JoinType, JoinTypeSquare, JoinTypeRound, JoinTypeMiter,
           EndType, EndTypeClosedPolygon, EndTypeClosedLine, EndTypeOpenSquare, EndTypeOpenRound, EndTypeOpenButt,
           Clip, add_path!, add_paths!, execute, clear!, get_bounds,
           IntPoint, IntRect, orientation, area, pointinpolygon, ClipperOffset,
           PolyNode, execute_pt, contour, ishole, contour, children

    @enum PolyType PolyTypeSubject=0 PolyTypeClip=1

    @enum ClipType ClipTypeIntersection=0 ClipTypeUnion=1 ClipTypeDifference=2 ClipTypeXor=3

    @enum PolyFillType PolyFillTypeEvenOdd=0 PolyFillTypeNonZero=1 PolyFillTypePositive=2 PolyFillTypeNegative=3

    @enum JoinType JoinTypeSquare=0 JoinTypeRound=1 JoinTypeMiter=2

    @enum EndType EndTypeClosedPolygon=0 EndTypeClosedLine=1 EndTypeOpenSquare=2 EndTypeOpenRound=3 EndTypeOpenButt=4

    struct IntPoint
        X::Int64
        Y::Int64
    end

    mutable struct PolyNode{T}
        contour::Vector{T}
        hole::Bool
        open::Bool
        children::Vector{PolyNode{T}}
        parent::PolyNode{T}
        PolyNode{T}(a,b,c) where {T} = new{T}(a,b,c)
        function PolyNode{T}(a,b,c,d) where T
            p = new{T}(a,b,c,d)
            p.parent = p
            return p
        end
        PolyNode{T}(a,b,c,d,e) where {T} = new{T}(a,b,c,d,e)
    end

    Base.convert(::Type{PolyNode{T}}, x::PolyNode{T}) where {T} = x
    function Base.convert(::Type{PolyNode{S}}, x::PolyNode{T}) where {S,T}
        parent(x) !== x && error("must convert a top-level PolyNode (i.e. a PolyTree).")

        pn = PolyNode{S}(convert(Vector{S}, contour(x)), ishole(x), isopen(x))
        pn.children = [PolyNode(y,pn) for y in children(x)]
        pn.parent = pn
    end
    function PolyNode(x::PolyNode, parent::PolyNode{S}) where S
        pn = PolyNode{S}(contour(x), ishole(x), isopen(x))
        pn.children = [PolyNode(y,pn) for y in children(x)]
        pn.parent = parent
        pn
    end

    @inline ishole(x::PolyNode) = x.hole
    @inline Base.isopen(x::PolyNode) = x.open
    @inline contour(x::PolyNode) = x.contour
    @inline children(x::PolyNode) = x.children
    @inline Base.parent(x::PolyNode) = x.parent

    function Base.show(io::IO, node::PolyNode)
        if parent(node) === node
            print(io, "Top-level PolyNode with $(length(children(node))) immediate children.")
        else
            if isopen(node)
                print(io, "Open ")
            else
                print(io, "Closed ")
            end
            print(io, "PolyNode ")
            ishole(node) && print(io, "(hole) ")
            println(io, "with contour:")
            show(io, contour(node))
            println(io, "")
            print(io, "...and $(length(children(node))) immediate children.")
        end
    end

    function Base.show(io::IO, point::IntPoint)
      print(io, "[$(point.X),$(point.Y)]")
    end

    function append!(outputArray::Ptr{Cvoid}, polyIndex::Csize_t, point::IntPoint)
        ourArray = unsafe_pointer_to_objref(outputArray)::Vector{Vector{IntPoint}}

        while (polyIndex + 1) > length(ourArray)
          push!(ourArray, Vector{IntPoint}())
        end

        push!(ourArray[polyIndex + 1], point)
    end

    # private
    function appendpn!(jl_node::Ptr{Cvoid}, point::IntPoint)
        node = unsafe_pointer_to_objref(jl_node)::PolyNode{IntPoint}
        push!(contour(node), point)
    end

    # private
    function newnode(outputTree::Ptr{Cvoid}, ishole::Bool, isopen::Bool)
        tree = unsafe_pointer_to_objref(outputTree)::PolyNode{IntPoint}
        node = PolyNode{IntPoint}(IntPoint[], ishole, isopen, PolyNode{IntPoint}[], tree)
        push!(children(tree), node)
        pointer_from_objref(node)
    end

    #==============================================================#
  	# Static functions
  	#==============================================================#
    function orientation(path::Vector{IntPoint})
        ccall((:orientation, cclipper), Cuchar, (Ptr{IntPoint}, Csize_t),
            path,
            length(path)) == 1 ? true : false
    end

    function area(path::Vector{IntPoint})
        ccall((:area, cclipper), Float64, (Ptr{IntPoint}, Csize_t),
            path,
            length(path))
    end

    function pointinpolygon(pt::IntPoint, path::Vector{IntPoint})
        ccall((:pointinpolygon, cclipper), Cint, (IntPoint, Ptr{IntPoint}, Csize_t),
            pt,
            path,
            length(path))
    end

    #==============================================================#
  	# Clipper object
  	#==============================================================#
    mutable struct Clip
        clipper_ptr::Ptr{Cvoid}

        function Clip()
            clipper = new(ccall((:get_clipper, cclipper), Ptr{Cvoid}, ()))
            finalizer(c -> ccall((:delete_clipper, cclipper), Cvoid, (Ptr{Cvoid},), c.clipper_ptr), clipper)
            clipper
        end
    end

    function add_path!(c::Clip, path::Vector{IntPoint}, polyType::PolyType, closed::Bool)
        ccall((:add_path, cclipper), Cuchar, (Ptr{Cvoid}, Ptr{IntPoint}, Csize_t, Cint, Cuchar),
              c.clipper_ptr,
              path,
              length(path),
              Int(polyType),
              closed) == 1 ? true : false
    end

    function add_paths!(c::Clip, paths::Vector{Vector{IntPoint}}, polyType::PolyType, closed::Bool)
        lengths = Vector{UInt64}()
        for path in paths
            push!(lengths, length(path))
        end

        ccall((:add_paths, cclipper), Cuchar, (Ptr{Cvoid}, Ptr{Ptr{IntPoint}}, Ptr{Csize_t}, Csize_t, Cint, Cuchar),
              c.clipper_ptr,
              paths,
              lengths,
              length(paths),
              Int(polyType),
              closed) == 1 ? true : false
    end

    function execute(c::Clip, clipType::ClipType, subjFillType::PolyFillType, clipFillType::PolyFillType)
        polys = Vector{Vector{IntPoint}}()

        result = ccall((:execute, cclipper), Cuchar, (Ptr{Cvoid}, Cint, Cint, Cint, Any, Ptr{Cvoid}),
                        c.clipper_ptr,
                        Int(clipType),
                        Int(subjFillType),
                        Int(clipFillType),
                        polys,
                        @cfunction(append!, Any, (Ptr{Cvoid}, Csize_t, IntPoint)))

        return result == 1 ? true : false, polys
    end

    function execute_pt(c::Clip, clipType::ClipType, subjFillType::PolyFillType, clipFillType::PolyFillType)
        pt = PolyNode{IntPoint}(IntPoint[], false, false, PolyNode{IntPoint}[])

        result = ccall((:execute_pt, cclipper), Cuchar,
            (Ptr{Cvoid}, Cint, Cint, Cint, Any, Ptr{Cvoid}, Ptr{Cvoid}),
            c.clipper_ptr,
            Int(clipType),
            Int(subjFillType),
            Int(clipFillType),
            pt,
            @cfunction(newnode, Ptr{Cvoid}, (Ptr{Cvoid}, Bool, Bool)),
            @cfunction(appendpn!, Any, (Ptr{Cvoid}, IntPoint)))

        return result == 1 ? true : false, pt
    end

    function clear!(c::Clip)
        ccall((:clear, cclipper), Cvoid, (Ptr{Cvoid},), c.clipper_ptr)
    end

    mutable struct IntRect
        left::Int64
        top::Int64
        right::Int64
        bottom::Int64
    end

    function get_bounds(c::Clip)
        ccall((:get_bounds, cclipper), IntRect, (Ptr{Cvoid}, ), c.clipper_ptr)
    end

    #==============================================================#
  	# ClipperOffset object
  	#==============================================================#
    mutable struct ClipperOffset
        clipper_ptr::Ptr{Cvoid}

        function ClipperOffset(miterLimit::Float64 = 2.0, roundPrecision::Float64 = 0.25)
            clipper = new(ccall((:get_clipper_offset, cclipper), Ptr{Cvoid}, (Cdouble, Cdouble), miterLimit, roundPrecision))
	    finalizer(c -> ccall((:delete_clipper_offset, cclipper), Cvoid, (Ptr{Cvoid},), c.clipper_ptr), clipper)

            clipper
        end
    end

    function add_path!(c::ClipperOffset, path::Vector{IntPoint}, joinType::JoinType, endType::EndType)
        ccall((:add_offset_path, cclipper), Cvoid, (Ptr{Cvoid}, Ptr{IntPoint}, Csize_t, Cint, Cint),
              c.clipper_ptr,
              path,
              length(path),
              Int(joinType),
              Int(endType))
    end

    function add_paths!(c::ClipperOffset, paths::Vector{Vector{IntPoint}}, joinType::JoinType, endType::EndType)
        lengths = Vector{UInt64}()
        for path in paths
            push!(lengths, length(path))
        end

        ccall((:add_offset_paths, cclipper), Cvoid, (Ptr{Cvoid}, Ptr{Ptr{IntPoint}}, Ptr{Csize_t}, Csize_t, Cint, Cint),
              c.clipper_ptr,
              paths,
              lengths,
              length(paths),
              Int(joinType),
              Int(endType))
    end

    function clear!(c::ClipperOffset)
        ccall((:clear_offset, cclipper), Cvoid, (Ptr{Cvoid},), c.clipper_ptr)
    end

    function execute(c::ClipperOffset, delta::Float64)
        polys = Vector{Vector{IntPoint}}()
        result = ccall((:execute_offset, cclipper), Cvoid, (Ptr{Cvoid}, Cdouble, Any, Ptr{Cvoid}),
                        c.clipper_ptr,
                        delta,
                        polys,
                        @cfunction(append!, Any, (Ptr{Cvoid}, Csize_t, IntPoint)))

        return polys
    end

    include("conversion.jl")
    include("utilities.jl")
end

ext: jl


text:
#=
Counterfactual regret minimizaion
- Chance sampling
- Pruning
- Imperfect recall
=#
using Random, Distributions, Printf
using Profile, BenchmarkTools
using Serialization

export randhands!,
       key,
       Key,
       actions,
       bids,
       normalize,
       win,
       train,
       Botprofile,
       Node

mutable struct Node
    regretsum::Dict{UInt8,Float64}
    strategysum::Dict{UInt8,Float64}
end

const Key = Tuple{Array{UInt8,1},Array{UInt8,1}}

key(hand::Vector{UInt8}, history::Vector{UInt8})::Key = (hand, history)

mutable struct Botprofile
    num_samples::Int
    cards_per_hand::Int
    num_ranks::Int
    max_history::Int
    min_bid_quant::Int
    max_bid_quant::Int
    all_bids::Vector{UInt8}
    policys::Dict{Key,Node}

    function Botprofile(ns, cph, r, mh, minq, maxq)
        if minq < 0
            @printf(
                "WARNING minimum bid quantity (%.0f) must be positive, defaulting to 0.\n",
                minq,
            )
            minq = 0
        end
        if maxq > 2 * cph
            @printf(
                "WARNING maximum bid quantity (%.0f) exceeds total number of cards, defaulting to %.0f\n",
                maxq,
                2 * cph,
            )
            maxq = 2 * cph
        end
        ab = mkbids(minq, maxq, r)
        new(ns, cph, r, mh, minq, maxq, ab, Dict{Key,Node}())
    end
end

function Base.show(io::IO, profile::Botprofile)
    @printf(io, "Samples:                   %.0f\n", profile.num_samples)
    @printf(io, "Cards per hand:            %.0f\n", profile.cards_per_hand)
    @printf(io, "Number of ranks:           %.0f\n", profile.num_ranks)
    @printf(io, "Max history:               %.0f\n", profile.max_history)
    @printf(io, "Minimum bid qunatitiy:     %.0f\n", profile.min_bid_quant)
    @printf(io, "Maximum bid qunatitiy:     %.0f\n", profile.max_bid_quant)
end

function mkbids(min_quant, max_quant, ranks)
    m = UInt8(min_quant)
    n = UInt8(max_quant)
    r = UInt8(ranks)
    vcat([0xa * i + j for j = 0x0:(r-0x1), i = m:n]...)
end

pos(x) = x > 0 ? x : zero(x)

issomething(s) = !isnothing(s)

Base.map(f, dict::AbstractDict) = Dict(k => f(v) for (k, v) in dict)
Base.show(io::IO, x::UInt8) = show(io, Int(x))
Base.show(io::IO, x::Vector{UInt8}) = show(io, map(n -> Int(n), x))

function normalize(xs)
    total = values(xs) |> sum
    n = length(xs)
    total > 0 ? map(x -> x / total, xs) : map(_ -> 1.0 / n, xs)
end

parsebid(n) = divrem(n, 10)

function lastn!(hist, n)
    if length(hist) > n
        popfirst!(hist)
    end
    if hist[1] == 0 # cannot start with a challenge.
        popfirst!(hist)
    end
end

function getstrategy(node, realizationweight)
    strategy = map(pos, node.regretsum) |> normalize
    for (k, v) in strategy
        node.strategysum[k] += realizationweight * v
    end
    strategy
end

function actions(bids, history)
    n = length(history)
    if n > 2 && history[end-2] == 0 && history[end] == 0
        return [1]
    elseif n > 1 && history[end] == 0
        b = findfirst(isequal(history[end-1]), bids) + 1
        return pushfirst!(bids[b:end], 1)
    elseif n > 0
        b = findfirst(isequal(history[end]), bids) + 1
        return pushfirst!(bids[b:end], 0)
    else
        return bids
    end
end

function getnode(profile, hand, history)
    k = key(hand, history)
    if !haskey(profile.policys, k)
        dict()::Dict{UInt8,Float64} =
            Dict(k => 0.0 for k in actions(profile.all_bids, history))
        node = Node(dict(), dict())
        profile.policys[k] = node
    end
    profile.policys[k]
end

function win(history, hands)
    (quant, rank) = history[end-2] |> parsebid
    count(c -> c == rank, hands) >= quant
end

function terminal(history, hands)
    plays = length(history)
    if plays > 2 && history[end] == 1
        return win(history, hands) ? 1 : -1
    end
    return nothing
end

function cfr(profile, hands, history, p1, p2)::Float64
    terminalutility = terminal(history, hands)
    if issomething(terminalutility)
        return -terminalutility
    end
    player = length(history) % 2 + 1
    if p1 <= 0 && p2 <= 0
        return 0.0
    end
    nodeutil = 0.0
    prob = [p1, p2][player]
    node = getnode(profile, hands[player, :], history)
    strategy = getstrategy(node, prob)
    util = Dict{UInt8,Float64}()
    acts = keys(strategy)
    for a in acts
        nexthistory = push!(copy(history), a)
        lastn!(nexthistory, profile.max_history)
        util[a] = player == 1 ?
                  -cfr(profile, hands, nexthistory, p1 * strategy[a], p2) :
                  -cfr(profile, hands, nexthistory, p1, p2 * strategy[a])
        nodeutil += strategy[a] * util[a]
    end
    q = [p2, p1][player]
    for a in acts
        node.regretsum[a] += q * (util[a] - nodeutil)
    end
    nodeutil
end

function randhands!(profile, a::Matrix{UInt8})
    for i in eachindex(a)
        a[i] = rand(0x0:(profile.num_ranks-0x1))
    end
    sort!(a, dims = 2)
    return a
end

function train(profile, filename=nothing)
    n = profile.num_samples
    util = 0.0
    hs = Matrix{UInt8}(undef, (2, profile.cards_per_hand))
    for i = 1:n
        if i % 1000 == 0
            println(i)
        end
        randhands!(profile, hs)
        util += cfr(profile, hs, UInt8[], 1.0, 1.0)
    end
    println(util / n)
    if issomething(filename)
        serialize(filename, profile)
    end
    profile
end

max_stars_repo_name: jeffreyrosenbluth/LP

max_issues_repo_name: jeffreyrosenbluth/LP

max_forks_repo_name: jeffreyrosenbluth/LP


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

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, StatsBase, DataFrames, Plots
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.8"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
begin
	grade = ["A","B","C","D","E"]
	weightVect = Weights([0.1,0.2,0.1,0.2,0.4])

	N = 10^6
	data = sample(grade,weightVect,N)
	[count(i->(i==g),data) for g in grade]/N
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.8"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
import Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using DelimitedFiles
using LinearAlgebra
using Printf
using SparseArrays

using BSON
using HDF5
using Interpolations
using OrdinaryDiffEq, DiffEqCallbacks
using Postprocessing
using UnPack
using RecursiveArrayTools
using Roots
using StaticArrays
using SummationByPartsOperators
import SymPy; sp=SymPy

import FFTW; FFTW.set_num_threads(1)

using PyCall, LaTeXStrings; import PyPlot; plt=PyPlot
inset_locator = pyimport("mpl_toolkits.axes_grid.inset_locator")

cycler = pyimport("cycler").cycler
line_cycler   = (cycler(color=["#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442"]) +
                 cycler(linestyle=["-", "--", "-.", ":", "-", "--", "-."]))
marker_cycler = (cycler(color=["#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442"]) +
                 cycler(linestyle=["none", "none", "none", "none", "none", "none", "none"]) +
                 cycler(marker=["4", "2", "3", "1", "+", "x", "."]))

plt.rc("axes", prop_cycle=line_cycler)
plt.rc("text", usetex=true)
plt.rc("text.latex", preamble="\\usepackage{newpxtext}\\usepackage{newpxmath}\\usepackage{commath}\\usepackage{mathtools}")
plt.rc("font", family="serif", size=18.)
plt.rc("savefig", dpi=100)
plt.rc("legend", loc="best", fontsize="medium", fancybox=true, framealpha=0.5)
plt.rc("lines", linewidth=2.5, markersize=10, markeredgewidth=2.5)

function relaxation!(integrator)
    told = integrator.tprev
    uold = integrator.uprev
    tnew = integrator.t
    unew = integrator.u

    γ = one(tnew)
    terminate_integration = false
    γlo = one(γ)/2
    γhi = 3*one(γ)/2
    energy_old = relaxation_functional(uold, integrator.p)
    sign_condition = (relaxation_functional(γlo, unew, uold, integrator.p)-energy_old) * (relaxation_functional(γhi, unew, uold, integrator.p)-energy_old)
    if iszero(sign_condition)
        # everything is fine
        γ = one(tnew)
    elseif sign_condition > 0
        terminate_integration = true
        @warn "Terminating integration because no solution γ can be found."
    else
        γ = find_zero(g -> relaxation_functional(g, unew, uold, integrator.p)-energy_old, (γlo, γhi), Roots.AlefeldPotraShi())
    end
    if γ < eps(typeof(γ))
        terminate_integration = true
        @warn "Terminating integration because γ=$γ is too small."
    end

#     println(γ)
    @. unew = uold + γ * (unew - uold)
    DiffEqBase.set_u!(integrator, unew)
    if !(tnew ≈ first(integrator.opts.tstops))
        tγ = told + γ * (tnew - told)
        DiffEqBase.set_t!(integrator, tγ)
    end

    if terminate_integration
        terminate!(integrator)
    end

    nothing
end


text:
@testset "get_prenex_normal_form" begin
    @customtest String(get_renamed_quantifiers_form(Formula("(~?{y}*{g}P(f(y, n))) & (?{x}P(x) | *{y}P(y))"))) ==
        "(~?{##1}*{##2}P(f(##1, n)) & (?{##3}P(##3) | *{##4}P(##4)))"
    @customtest String(get_renamed_quantifiers_form(Formula("*{x}?{y}P(x, g(y, f(x))) | ~Q(z) | ~*{x}R(x,y)"))) ==
        "(*{##1}?{##2}P(##1, g(##2, f(##1))) | (~Q(z) | ~*{##3}R(##3, y)))"
end


text:
#############################################################################
# quantum_relative_entropy returns tr(A*(log(A)-log(B))) where A and B
# are positive semidefinite matrices.  Note this function uses logarithm
# base e, not base 2, so return value is in units of nats, not bits.
#
# Quantum relative entropy is convex (jointly) in (A,B). This function
# implements the semidefinite programming approximation given in the reference
# below.  Parameters m and k control the accuracy of this approximation: m is
# the number of quadrature nodes to use and k the number of square-roots to
# take. See reference for more details.
#
# Implementation uses the expression
#    D(A||B) = e'*D_{op} (A \otimes I || I \otimes B) )*e
# where D_{op} is the operator relative entropy and e = vec(Matrix(I, n, n)).
#
# All expressions and atoms are subtypes of AbstractExpr.
# Please read expressions.jl first.
#
#REFERENCE
#   Ported from CVXQUAD which is based on the paper: "Lieb's concavity
#   theorem, matrix geometric means and semidefinite optimization" by Hamza
#   Fawzi and James Saunderson (arXiv:1512.03401)
#############################################################################

const MatrixOrConstant = Union{AbstractMatrix,Constant}

struct QuantumRelativeEntropy1 <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr,AbstractExpr}
    size::Tuple{Int,Int}
    m::Integer
    k::Integer

    function QuantumRelativeEntropy1(
        A::AbstractExpr,
        B::AbstractExpr,
        m::Integer,
        k::Integer,
    )
        children = (A, B)
        if size(A) != size(B)
            throw(DimensionMismatch("A and B must be the same size"))
        end
        n = size(A)[1]
        if size(A) != (n, n)
            throw(DimensionMismatch("A and B must be square"))
        end
        return new(
            :quantum_relative_entropy,
            hash(children),
            children,
            (1, 1),
            m,
            k,
        )
    end
end

struct QuantumRelativeEntropy2 <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int,Int}
    m::Integer
    k::Integer
    B::AbstractMatrix
    J::AbstractMatrix
    K::AbstractMatrix

    function QuantumRelativeEntropy2(
        A::AbstractExpr,
        B::AbstractMatrix,
        m::Integer,
        k::Integer,
        nullspace_tol::Real,
    )
        children = (A,)

        if size(A) != size(B)
            throw(DimensionMismatch("A and B must be the same size"))
        end
        n = size(A)[1]
        if size(A) != (n, n)
            throw(DimensionMismatch("A and B must be square"))
        end
        if norm(B - B') > nullspace_tol
            throw(DomainError(B, "B must be Hermitian"))
        end

        # nullspace of A must contain nullspace of B
        v, U = eigen(Hermitian(B))
        if any(v .< -nullspace_tol)
            throw(DomainError(B, "B must be positive semidefinite"))
        end
        J = U'[v.>nullspace_tol, :]
        K = U'[v.<nullspace_tol, :]

        return new(
            :quantum_relative_entropy,
            hash(children),
            children,
            (1, 1),
            m,
            k,
            B,
            J,
            K,
        )
    end
end

sign(atom::Union{QuantumRelativeEntropy1,QuantumRelativeEntropy2}) = Positive()

function monotonicity(atom::QuantumRelativeEntropy1)
    return (NoMonotonicity(), NoMonotonicity())
end
monotonicity(atom::QuantumRelativeEntropy2) = (NoMonotonicity(),)

function curvature(atom::Union{QuantumRelativeEntropy1,QuantumRelativeEntropy2})
    return ConvexVexity()
end

function evaluate(atom::QuantumRelativeEntropy1)
    A = evaluate(atom.children[1])
    B = evaluate(atom.children[2])
    return quantum_relative_entropy(A, B)
end

function evaluate(atom::QuantumRelativeEntropy2)
    A = evaluate(atom.children[1])
    return quantum_relative_entropy(A, atom.B)
end

function quantum_relative_entropy(
    A::AbstractExpr,
    B::AbstractExpr,
    m::Integer = 3,
    k::Integer = 3,
)
    #println("quantum_relative_entropy general case")
    return QuantumRelativeEntropy1(A, B, m, k)
end

function quantum_relative_entropy(
    A::AbstractExpr,
    B::MatrixOrConstant,
    m::Integer = 3,
    k::Integer = 3,
    nullspace_tol::Real = 1e-6,
)
    #println("quantum_relative_entropy general case")
    return QuantumRelativeEntropy2(A, evaluate(B), m, k, nullspace_tol)
end

function quantum_relative_entropy(
    A::MatrixOrConstant,
    B::AbstractExpr,
    m::Integer = 3,
    k::Integer = 3,
)
    #println("quantum_relative_entropy constant A")
    A = evaluate(A)
    return -quantum_entropy(A, m, k) - trace_logm(B, A, m, k)
end

function quantum_relative_entropy(
    A::MatrixOrConstant,
    B::MatrixOrConstant,
    m::Integer = 0,
    k::Integer = 0,
    nullspace_tol::Real = 1e-6,
)
    #println("quantum_relative_entropy constant A and B")
    A = evaluate(A)
    B = evaluate(B)

    if size(A) != size(B)
        throw(DimensionMismatch("A and B must be the same size"))
    end
    if size(A) != (size(A)[1], size(A)[1])
        throw(DimensionMismatch("A and B must be square"))
    end
    if norm(A - A') > nullspace_tol
        throw(DomainError(A, "A must be Hermitian"))
    end
    if norm(B - B') > nullspace_tol
        throw(DomainError(B, "B must be Hermitian"))
    end

    # need to project down to support of A
    v, U = eigen(Hermitian(A))
    if any(v .< -nullspace_tol)
        throw(DomainError(A, "A must be positive semidefinite"))
    end
    if any(eigvals(Hermitian(B)) .< -nullspace_tol)
        throw(DomainError(B, "B must be positive semidefinite"))
    end
    J = U'[v.>nullspace_tol, :]
    Ap = Hermitian(J * A * J')
    Bp = Hermitian(J * B * J')

    if any(eigvals(Bp) .< nullspace_tol)
        return Inf
    end

    return real(tr(Ap * (log(Ap) - log(Bp))))
end

function conic_form!(atom::QuantumRelativeEntropy1, unique_conic_forms)
    #println("conic_form QuantumRelativeEntropy1")
    if !has_conic_form(unique_conic_forms, atom)
        A = atom.children[1]
        B = atom.children[2]
        m = atom.m
        k = atom.k
        n = size(A)[1]
        eye = Matrix(1.0 * I, n, n)
        e = vec(eye)

        conic_form!(A ⪰ 0, unique_conic_forms)
        conic_form!(B ⪰ 0, unique_conic_forms)

        τ = Variable()
        conic_form!(
            τ in
            RelativeEntropyEpiCone(kron(A, eye), kron(eye, conj(B)), m, k, e),
            unique_conic_forms,
        )

        cache_conic_form!(unique_conic_forms, atom, minimize(τ))
    end
    return get_conic_form(unique_conic_forms, atom)
end

function conic_form!(atom::QuantumRelativeEntropy2, unique_conic_forms)
    #println("conic_form QuantumRelativeEntropy2")
    if !has_conic_form(unique_conic_forms, atom)
        A = atom.children[1]
        B = atom.B
        J = atom.J
        K = atom.K
        m = atom.m
        k = atom.k

        conic_form!(A ⪰ 0, unique_conic_forms)

        if length(K) > 0
            conic_form!(K * A * K' == 0, unique_conic_forms)
            Ap = J * A * J'
            Bp = Hermitian(J * B * J')
            τ = -quantum_entropy(Ap, m, k) - real(tr(Ap * log(Bp)))
        else
            τ = -quantum_entropy(A, m, k) - real(tr(A * log(B)))
        end

        cache_conic_form!(unique_conic_forms, atom, minimize(τ))
    end
    return get_conic_form(unique_conic_forms, atom)
end


text:
# ======udp========= need to be finished
const UDPSize = 65536

macro randPort()
    quote
        UInt16(floor(rand()*10000 + 50000))
    end
end

function doBind(port::Integer)
    udp = UDPSocket()
    try
        bind(udp, getipaddr() isa IPv4 ? IPv4(0) : IPv6(0), @randPort)
    catch err
        return nothing, err
    end

    return udp, nothing
end

mutable struct S5UDPSocket
    sock::Union{UDPSocket, Nothing}
end

mutable struct SSUDPSocket
    sock::Union{UDPSocket, Nothing}
    cipher::Union{Cipher, Nothing}
    nonce::Union{Bytes, Nothing}
end

mutable struct NATmap
    map::Dict{Bytes, Union{UDPSocket, SSUDPSocket, Nothing}}
    timeout::Integer
end

import Base.in
function in(nmap::NATmap, key::Bytes)
    return haskey(nmap.map, key)
end

import Sockets.send
import Sockets.recv

function send(sock::SSUDPSocket, ip::Sockets.IPAddr, port::Integer)
end

function recv(sock::SSUDPSocket)
end

function waitWithTimeout()
end

function gethost(buff::Bytes, shift::Integer)
    host, port = nothing, nothing
    if buff[1+shift] == 0x01
        host = IPv4(ntoh(unsafe_load(Ptr{UInt32}(pointer(buff[2+shift:5+shift])))))
        port = UInt16(buff[6+shift]) << 8 + buff[7+shift]
    elseif buff[1+shift] == 0x03
        len = buff[2+shift]
        host = try
            getaddrinfo(String(buff[3+shift:len+2+shift]))
        catch err
            return nothing, nothing, err
        end
        port = UInt16(buff[len+3+shift]) << 8 + buff[len+4+shift]
    elseif buff[1+shift] == 0x04
        host = IPv6(ntoh(unsafe_load(Ptr{UInt128}(pointer(buff[2+shift:17+shift])))))
        port = UInt16(buff[18+shift]) << 8 + buff[19+shift]
    end

    return host, port, nothing
end

function key(data::Bytes, shift::Integer)
    if data[1+shift] == 0x01
        return data[2+shift:7+shift]
    elseif data[1+shift] == 0x03
        len = data[2+shift]
        return data[3+shift:len+4+shift]
    elseif data[1+shift] == 0x04
        return data[2+shift:19+shift]
    end
end

function udpServer(config::SSConfig, cipher::Cipher)
    udp, err = doBind(config.lisPort)
    server = SSUDPSocket(udp, cipher, zeros(UInt8, cipher.ivlen))
    nmap = NATmap(Dict(UInt8[0,0,0,0,0,0] => nothing), 60)

    while true
        data, err = recv(server)
        begin
            host, port, err = gethost(data, 0)
            thiskey = key(data)
            if thiskey in nmap
            end

            @async waitWithTimeout()
        end
    end
end
# =====udp==========


text:
#!/usr/bin/env julia6
using ReadGlobal, MyStats, Fluids

function dinfo(file::String)
    field = readfield(file)
    nx,ny,nz,_,_,_ = getdimsize()
    fdims = (nx,ny,nz)
    med = tmean(field)
    devi = tstd(field,med)
    imin,imax=  min_max_ind(field)
    fmax = field[imax]
    imaxc = ind2sub(fdims,imax)
    fmin = field[imin]
    iminc = ind2sub(fdims,imin)

    println("Mean: ",med)
    println("Std: ",devi)
    println("Max: ",fmax," at ",imaxc)
    println("Min: ",fmin," at ",iminc)
    return 0
end

function dinfo(file1::String,file2::String,file3::String)
    field = FVector(file1,file2,file3)
    nx,ny,nz,_,_,_ = getdimsize()
    fdims = (nx,ny,nz)
    med = mean(field)
    devi = std(field,mean=med)
    imax =  indmax(mod(field))
    fmax = field[imax]
    imax = ind2sub(fdims,imax)
    imin =  indmin(mod(field))
    fmin = field[imin]
    imin = ind2sub(fdims,imin)

    println("Mean: ",med)
    println("Std: ",devi)
    println("Max: ",fmax," at ",imax)
    println("Min: ",fmin," at ",imin)
    return 0
end
                 #d11           d12           d13           d22   d23    d33
function dinfo(file1::String,file2::String,file3::String,file4,file5,file6)
    field = SymTensor(file1,file4,file6,file2,file3,file5)
    nx,ny,nz,_,_,_ = getdimsize()
    fdims = (nx,ny,nz)
    med = mean(field)
    devi = std(field,mean=med)
    imax =  indmax(mod(field))
    fmax = field[imax]
    imax = ind2sub(fdims,imax)
    imin =  indmin(mod(field))
    fmin = field[imin]
    imin = ind2sub(fdims,imin)

    println("Mean: ",med)
    println("Std: ",devi)
    println("Max: ",fmax," at ",imax)
    println("Min: ",fmin," at ",imin)
    return 0
end

dinfo(ARGS...)


text:
#=
using PETSc
using BaseTestNext
# test timestepping (TS) functions
ST = Float64
=#

@testset "TS" begin
  # solve the equation f =  [u1;u2] = [t + 1; 2t + 1] using forward Euler
  # This should be an exact case

  # the jacobian is always the identity for a case where u = f(t) only
  A = Mat(ST, 2, 2)
  A[1,1] = 1.0
  A[2,2] = 1.0

  # evaluate rhs u_t
  function rhs1(ts::TS, t, u::Vec, f::Vec, ctx)

    f_local = LocalVector(f)

    f_local[1] = 1.0
    f_local[2] = 2.0
    restore(f_local)

    return 0
  end

  ts = TS(ST, PETSc.C.TS_LINEAR, PETSc.C.TSEULER)

  set_rhs_function(ts, NullVec[ST], rhs1)
  set_rhs_jac(ts, A, A, ComputeRHSJacobianConstant)

  set_times(ts, 0.0, 0.1, 10, 50.0)  # limited by number of timesteps

  # hold IC/final solution
  u = Vec(ST, 2)
  u[1] = 1.0
  u[2] = 1.0

  solve!(ts, u)
  #not sure how to test this other than "doesn't explode everywhere"
  ksp = KSP(ts)

  # same as above, but using set_ic
  ts = TS(ST, PETSc.C.TS_LINEAR, PETSc.C.TSEULER)

  set_rhs_function(ts, NullVec[ST], rhs1)
  set_rhs_jac(ts, A, A, ComputeRHSJacobianConstant)

  set_times(ts, 0.0, 0.1, 10, 50.0)  # limited by number of timesteps

  u = Vec(ST, 2)
  u[1] = 1.0
  u[2] = 1.0
  set_ic(ts, u)
  solve!(ts)

#  petscview(u)
  @test u[1] ≈ ST(2.0)
  @test u[2] ≈ ST(3.0)

  # solve the 1D heat equation u_t = u_xx with 2nd order space, linear time
  # homogneous Dirchlet BCs

  dt = 0.05
  nsteps = 10
  dx = 0.2  # actually dx^2
  nx = 5  # number of points in domain (including boundaries
  r = dt/dx

  A_julia = zeros(ST, nx, nx)
  A_julia[1,1] = 1.0
  A_julia[nx,nx] = 1.0

  # interior discretization
  for i=2:(nx-1)
    A_julia[i, i-1] = r
    A_julia[i, i] = 1 - 2*r
    A_julia[i, i+1] = r
  end

  u_julia = ones(ST, nx)
  u_julia[1] = 0.0  # homogenous direchlet
  u_julia[nx] = 0.0

  # solve in julia
  for t = 1:(nsteps+0)
    u_next = A_julia*u_julia
    u_julia[:] = u_next
  end

  # now solve in Petsc

  function rhs_jacobian(ts::TS, t, u, A::Mat, B::Mat, ctx)
    dx = ctx[1]
    nx, ny = size(A)
    A[1,1] = 1.0
    A[nx, nx] = 1.0

    for i=2:(nx-1)
      A[i, i-1] = 1/dx
      A[i, i] = -2/dx
      A[i, i+1] = 1/dx
    end

    return 0
  end

  function rhs_func(ts::TS, t, u::Vec, F::Vec, ctx)
    dx = ctx[1]
    nx = length(u)
    u_arr = LocalVector_readonly(u)
    f_arr = LocalVector(F)

    f_arr[1] = 0.0
    f_arr[nx] = 0.0

    for i=2:(nx-1)
      f_arr[i] = u_arr[i-1]/dx + -2*u_arr[i]/dx + u_arr[i+1]/dx
    end

    restore(u_arr)
    restore(f_arr)

    return 0
  end

  ts = TS(ST, PETSc.C.TS_LINEAR, PETSc.C.TSEULER)

  A = Mat(ST, nx, nx)
  u = Vec(ST, nx)

  # set IC
  for i=2:(nx-1)
    u[i] = 1.0
  end

  ctx1 = (dx,)
  ctx2 = (dx,)
  set_rhs_function(ts, NullVec[ST], rhs_func, ctx1)
  set_rhs_jac(ts, A, A, rhs_jacobian, ctx2)

  set_times(ts, 0.0, dt, nsteps , 50.0)  # limited by number of timesteps

  solve!(ts, u)

  for i=1:nx
    @test u[i] ≈ u_julia[i]
  end

  PETSc.PetscDestroy(ts)
  @test PETSc.isfinalized(ts)

  # now solve the same problem using implicit euler
  dt = 0.05
  nsteps = 10
  dx = 0.2  # actually dx^2
  nx = 5  # number of points in domain (including boundaries
  r = dt/dx

  A_julia = zeros(ST, nx, nx)
  A_julia[1,1] = 1.0
  A_julia[nx,nx] = 1.0

  # interior discretization
  for i=2:(nx-1)
    A_julia[i, i-1] = -r
    A_julia[i, i] = 1 + 2*r
    A_julia[i, i+1] = -r
  end

  u_julia = ones(ST, nx)
  u_julia[1] = 0.0  # homogenous direchlet
  u_julia[nx] = 0.0

  # solve in julia
  for t = 1:(nsteps+0)
    u_next = A_julia\u_julia
    u_julia[:] = u_next
  end

  function lhs_func(ts::TS, t, u::Vec, ut::Vec, F::Vec, ctx)
    dx = ctx[1]
    nx = length(u)
    u_arr = LocalVector_readonly(u)
    ut_arr = LocalVector_readonly(ut)
    f_arr = LocalVector(F)

    f_arr[1] = u_arr[1]
    f_arr[nx] = u_arr[nx]

    for i=2:(nx-1)
      f_arr[i] = ut[i] - (u_arr[i-1]/dx + -2*u_arr[i]/dx + u_arr[i+1]/dx)
    end

    restore(u_arr)
    restore(ut_arr)
    restore(f_arr)

    return 0
  end

  function lhs_jacobian(ts::TS, u::Vec, ut::Vec, a, A::Mat, B::Mat, ctx)
    dx = ctx[1]
    nx, ny = size(A)
    A[1,1] = 1.0
    A[nx, nx] = 1.0

    for i=2:(nx-1)
      A[i, i-1] = -1/dx
      A[i, i] = 2/dx + a
      A[i, i+1] = -1/dx
    end

    AssemblyBegin(A, PETSc.C.MAT_FINAL_ASSEMBLY)
    AssemblyEnd(A, PETSc.C.MAT_FINAL_ASSEMBLY)

    return 0
  end

  ts = TS(ST, PETSc.C.TS_LINEAR, PETSc.C.TSBEULER)

  A = Mat(ST, nx, nx, nz = 3)
  A[1,1] = 1.0
  A[nx, nx] = 1.0
  for i=2:(nx-1)
    A[i, i-1] = 1.0
    A[i, i] = 1.0
    A[i, i+1] = 1.0
  end
  AssemblyBegin(A, PETSc.C.MAT_FINAL_ASSEMBLY)
  AssemblyEnd(A, PETSc.C.MAT_FINAL_ASSEMBLY)
  u = Vec(ST, nx)

  # set IC
  for i=2:(nx-1)
    u[i] = 1.0
  end

  ctx1 = (dx,)
  ctx2 = (dx,)
  set_lhs_function(ts, NullVec[ST], lhs_func, ctx1)
  set_lhs_jac(ts, A, A, lhs_jacobian, ctx2)

  set_times(ts, 0.0, dt, nsteps , 50.0)  # limited by number of timesteps
  AssemblyBegin(u)
  AssemblyEnd(u)
  solve!(ts, u)

  for i=1:nx
    @test u[i] ≈ u_julia[i]
  end

  PETSc.PetscDestroy(ts)

end


text:
using SQLQuery

type NewNode{T} <: SQLQuery.QueryNode
    input::T
    args
end

SQLQuery.QUERYNODE[:newnode] = NewNode

SQLQuery.translatesql(nn::NewNode, offset::Int) = "newnode, offset=$offset"

@sqlquery source |>
newnode(are, you.serious) |>
select(*, some, columns)

@sqlquery source |>
filter(name == 3, bar == "whee")

@sqlquery source |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col)

@sqlquery source |>
filter(name == 3, bar == "whee") |>
select(distinct(name = foo * 3, col))

@sqlquery source |>
filter(name == 3, bar == "whee") |>
select(all(name = foo * 3, col))

@sqlquery source |>
filter(name == 3, bar == "whee") |>
select(all(*))

@sqlquery source |>
filter(name == 3, bar == "whee") |>
select(distinct(*))

@sqlquery source |>
filter(name == 3, bar == "whee") |>
select(distinct(col))

@sqlquery source |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col) |>
orderby(name)

@sqlquery source |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col) |>
orderby(name, col)

@sqlquery source |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col) |>
orderby(desc(name))

@sqlquery source |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col) |>
orderby(name, desc(col))

@sqlquery source |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col) |>
orderby(desc(name), desc(col))

@sqlquery source |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col) |>
orderby(desc(name), asc(col))

@sqlquery source |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col) |>
orderby(asc(name), asc(col))

@sqlquery source |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col) |>
orderby(desc(name))

@sqlquery source |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col) |>
orderby(desc(name)) |>
limit(10)

@sqlquery source |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col) |>
orderby(desc(name)) |>
offset(7) |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col) |>
orderby(desc(name)) |>
offset(7) |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col) |>
orderby(desc(name)) |>
offset(7) |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col) |>
orderby(desc(name)) |>
offset(7) |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col) |>
orderby(desc(name)) |>
offset(7) |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col) |>
orderby(desc(name)) |>
offset(7) |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col) |>
orderby(desc(name)) |>
offset(7) |>
filter(name == 3, bar == "whee") |>
select(name = foo * 3, col) |>
orderby(desc(name)) |>
offset(7)

@sqlquery Artists |>
leftjoin(Songs) |>
leftjoin(Albums) |>
select( song_id = Songs._id,
song_name = Songs.name,
Songs.length,
artist_id = Songs.artist_id,
artist_name = Artists.name,
album_id = Songs.album_id,
album_name = Albums.name)

@sqlquery Artists |>
leftjoin(Songs) |>
leftjoin(Albums) |>
select( song_id = Songs._id,
song_name = Songs.name,
Songs.length,
artist_id = Songs.artist_id,
artist_name = Artists.name,
album_id = Songs.album_id,
album_name = Albums.name) |>
groupby(Songs.length, artist_id)

@sqlquery Artists |>
leftjoin(Songs) |>
leftjoin(Albums) |>
select( song_id = Songs._id,
song_name = Songs.name,
Songs.length,
artist_id = Songs.artist_id,
artist_name = Artists.name,
album_id = Songs.album_id,
album_name = Albums.name) |>
groupby(Songs.length, artist_id) |>
orderby(asc(song_id), desc(Songs.length), desc(artist_id))


text:
module MegTools

using AxisKeys
using Infiltrator
using LinearAlgebra
using MAT
using Measures
using NamedDims
using Statistics

# Paths are based on where this module is
# Once you're done with your script and feel you want to generalise it,
# then the functions are put here or 'include'-ed  here
# E.g. include("../scripts/ReadContEps.jl")

include("IO_Meg.jl")
export load_besa_av, load_cont_epochs, load_BSepochs

include("analysis_recipes.jl")
export average_across_trials, select_channels, baseline_correction, find_peaks, collect_peaks,
    find_mean_amplitude, collect_mean_amps, get_averaged_baseline

include("filters.jl")
export highlow_butterworth_filter

include("trigger_keys.jl")
export load_trigger_values

end

max_stars_repo_name: ElectronicTeaCup/MegTools

max_issues_repo_name: ElectronicTeaCup/MegTools

max_forks_repo_name: ElectronicTeaCup/MegTools


text:
export _run_vd

# mutable struct VDPowerModel <: _PM.AbstractACModel _PM.@pm_fields end

"""
run model for Voltge-Deviation objective with AC Power Flow equations
"""

function _run_vd(file, model_type::_PM.Type, optimizer; kwargs...)
    return _PM.run_model(file, model_type, optimizer, _build_vd; kwargs...)
end

"""
given a JuMP model and a PowerModels network data structure,
builds an Voltge-Deviation formulation of the given data and returns the JuMP model
"""

function _build_vd(pm::_PM.AbstractPowerModel)

    _PM.variable_bus_voltage(pm)
    _PM.variable_gen_power(pm)
    _PM.variable_branch_power(pm)
    _PM.variable_dcline_power(pm, bounded = false) # TODO: why false?

    objective_vd(pm)

    _PM.constraint_model_voltage(pm)

    for i in _PM.ids(pm, :ref_buses)
        _PM.constraint_theta_ref(pm, i)
    end

    for i in _PM.ids(pm, :bus)
        _PM.constraint_power_balance(pm, i)
    end

    for (i, branch) in _PM.ref(pm, :branch)
        _PM.constraint_ohms_yt_from(pm, i)
        _PM.constraint_ohms_yt_to(pm, i)

        _PM.constraint_thermal_limit_from(pm, i)
        _PM.constraint_thermal_limit_to(pm, i)
    end

    for i in _PM.ids(pm, :dcline)
        _PM.constraint_dcline_power_losses(pm, i)
    end
end

function objective_vd(pm::_PM.AbstractPowerModel)
    return JuMP.@objective(pm.model, Min,
        sum((var(pm, :vm, content["element_index"]) - content["value"])^2 for (i, content) in pm.ext[:setpoint_v]))
end


text:
const NUMBER_OF_LETTERS = 26

function rotate(key, input::Char)
    if input ∉ 'a':'z' && input ∉ 'A':'Z'
        return input
    end
    first_letter = islowercase(input) ? 'a' : 'A'
    initial_pos = input - first_letter
    return first_letter + (initial_pos + key) % NUMBER_OF_LETTERS
end

function rotate(key, input::String)
    return map(c -> rotate(key, c), input)
end

macro R13_str(input)
    rotate(13, input)
end

for i in 0:26
    @eval macro $(Symbol("R", i, "_str"))(input)
        rotate($i, input)
    end
end

max_stars_repo_name: vladflore/exercism

max_issues_repo_name: vladflore/exercism

max_forks_repo_name: vladflore/exercism


text:
using BenchmarkTools: @btime
import Random
using Distributions: Distribution, MixtureModel, Normal, LogNormal, pdf, ncomponents, probs, ContinuousUnivariateDistribution
import Plots

gmm = MixtureModel(
    Normal.([-1.0, 0.0, 3.0], # mean vector
            [0.3, 0.5, 1.0]),  # std vector
    [0.25, 0.25, 0.5] # component weights
)

xs  = -2.0:0.01:6.0
Plots.plot(xs, pdf.(gmm, xs), legend=nothing)
Plots.ylabel!("\$f_X(x)\$")
Plots.xlabel!("\$x\$")
Plots.title!("Gaussian mixture PDF")
Plots.savefig(joinpath(base_img, "gaussian_mm.pdf"))


text:
# TODO: This function might not be needed! Without @avx it can be merged with
# the two dimensional case
function downconvert!(
    downconverted_signal::StructArray{Complex{T}, 1},
    signal::StructArray{Complex{T}, 1},
    carrier_replica::StructArray{Complex{T}, 1},
    start_sample::Integer,
    num_samples::Integer
) where T
    ds_re = downconverted_signal.re; ds_im = downconverted_signal.im
    s_re = signal.re; s_im = signal.im
    c_re = carrier_replica.re; c_im = carrier_replica.im
    @avx for i = start_sample:num_samples + start_sample - 1
        ds_re[i] = s_re[i] * c_re[i] + s_im[i] * c_im[i]
        ds_im[i] = s_im[i] * c_re[i] - s_re[i] * c_im[i]
    end
    downconverted_signal
end

function downconvert!(
    downconverted_signal::StructArray{Complex{T}, 2},
    signal::StructArray{Complex{T}, 2},
    carrier_replica::StructArray{Complex{T}, 1},
    start_sample::Integer,
    num_samples::Integer
) where T
    ds_re = downconverted_signal.re; ds_im = downconverted_signal.im
    s_re = signal.re; s_im = signal.im
    c_re = carrier_replica.re; c_im = carrier_replica.im
    @avx for i = start_sample:num_samples + start_sample - 1, j = 1:size(s_re, 2)
        ds_re[i, j] = s_re[i, j] * c_re[i] + s_im[i, j] * c_im[i]
        ds_im[i, j] = s_im[i, j] * c_re[i] - s_re[i, j] * c_im[i]
    end
    downconverted_signal
end

function downconvert!(
    downconverted_signal,
    signal,
    carrier_replica,
    start_sample::Integer,
    num_samples::Integer
)
    sample_range = start_sample:start_sample + num_samples - 1
    downconverted_signal[sample_range,:] .= @view(signal[sample_range,:]) .* conj.(@view(carrier_replica[sample_range]))
    downconverted_signal
end


text:
Pkg.add("Vega")
using Vega

#see also: https://vega.github.io/editor/#/examples/vega/word-cloud

corpus = [
"Julia is a high-level, high-performance dynamic programming language for technical computing, with syntax that is familiar to users of other technical computing environments. It provides a sophisticated compiler, distributed parallel execution, numerical accuracy, and an extensive mathematical function library. Julia’s Base library, largely written in Julia itself, also integrates mature, best-of-breed open source C and Fortran libraries for linear algebra, random number generation, signal processing, and string processing. In addition, the Julia developer community is contributing a number of external packages through Julia’s built-in package manager at a rapid pace. IJulia, a collaboration between the IPython and Julia communities, provides a powerful browser-based graphical notebook interface to Julia.",
"Julia programs are organized around multiple dispatch; by defining functions and overloading them for different combinations of argument types, which can also be user-defined. For a more in-depth discussion of the rationale and advantages of Julia over other systems, see the following highlights or read the introduction in the online manual."
]

corpus2 = ["2d", "mesh", "polygon", "2d", "scatter", "image", "subscene", "2d",
"contour", "Anthony", "2d", "contour", "2d", "contour", "2d", "heatmap", "updating",
"2d", "scatter", "animation", "2d", "text", "align", "rotation", "2d", "text",
"scatter", "3d", "similar", "animated", "meshscatter", "axi s", "mesh", "3d",
"views", "scatter", "axis", "lines", "3d", "meshscatter", "3d", "animated",
"video", "axis", "wireframe", "3d", "surface", "mesh", "cat", "3d", "linesegment",
"mesh", "3d", "glow", "scatter", "3d", "scatter", "axis", "marker", "3d", "animated",
"offset", "lines", "3d", "gif", " axis", "3d", "surface", "text", "volume", "3d",
"documentation", "2d", "heatmap", "documentation", "mesh", "cat", "3d",
"documentation", "texture", "mesh", "cat", "3d", "documentation", "axis",
"mesh", "3d", "documentation", "mesh", "cat", "wireframe", "3d", "documentation",
"wireframe", "3d", " documentation", "wireframe", "3d", "documentation", "surface",
"3d", "documentation", "surface", "image", "3d", "documentation", "surface", "2d",
"lines", "documentation", "meshscatter", "3d", "documentation", "2d", "scatter",
"documentation", "Makie.scatter", "Anthony", "2d", "scatter", "docume ntation",
"Makie.scatter", "Anthony", "2d", "scatter", "documentation", "2d", "scatter",
"VideoStream", "documentation", "linesegment", "scatter", " lines", "3d",
"documentation", "linestyle", "legend", "colorlegend", "2d", "documentation",
"legend", "meshscatter", "3d", "documentation", "VideoStream", "linesegment",
"slices", "volume", "3d layout", "heatmap", "documentation", "layout", "contour",
"slices", "volume", "gui", "animation", "documentation", "layout"]

wc = wordcloud(x = corpus2)
colorscheme!(wc, palette = ("Spectral", 11))


text:
function expandkwargs(kwargs::Vector)
  m = length(kwargs)
  keys = Vector{Symbol}(m)
  vals = Vector{Any}(m)
  for i in 1:m
    keys[i] = kwargs[i][1]
    vals[i] = kwargs[i][2]
  end
  return keys, vals
end

macro attach(obj, fields...)
  block = Expr(:block)
  for f in fields
    e = :($f = $obj.$f)
    push!(block.args, e)
  end
  return esc(:($block))
end

const greekdict = Dict(:σ => :sigma, :ρ => :rho, :η => :eta)

function greekin(g::Symbol, c::AbstractArray)
  if g in c
    return true
  else
    s = get(greekdict, g, :notfound)
    if s == :notfound
      error("an unknown greek letter is given")
    else
      return s in c
    end
  end
end

function greekfind(c::AbstractArray, g::Symbol)
  if g in c
    return findfirst(c, g)
  else
    s = get(greekdict, g, :notfound)
    if s == :notfound
      error("an unknown greek letter is given")
    else
      return findfirst(c, s)
    end
  end
end

function ∫(f::Function, a::Number, b::Number; e::Float64=1e-6)
  val, err = quadgk(f, a, b, reltol=e)
  return val
end

# Ref: http://www.eweb.unex.es/eweb/fisteor/santos/files/rdf.nb
function inverselaplace(F::Function, t::Float64, nterm::Int, meuler::Int)
  a::Float64 = 30.0
  h::Float64 = π/t
  u::Float64 = exp(a/2) / t
  x::Float64 = a / 2t

  c(m::Int)::Int = binomial(meuler, m-1)
  fs(s::Complex{Float64})::Float64 = real(F(s))

  suma::Float64 = fs(x + 0im)/2
  @simd for n in 1:nterm
    @inbounds suma += (-1)^n * fs(x + im*n*h)
  end

  su = Vector{Float64}(meuler+2)
  su[1] = suma
  for k in 1:meuler+1
    n = nterm + k
    @inbounds su[k+1] = su[k] + (-1)^n * fs(x + im*n*h)
  end

  # argsu::Float64 = 0
  argsu1::Float64 = 0
  @simd for j in 1:meuler+1
    # @inbounds argsu = argsu + c(j)*su[j]
    @inbounds argsu1 += c(j)*su[j+1]
  end

  # fun::Float64 = u * argsu / 2^meuler
  fun1::Float64 = u * argsu1 / 2^meuler

  return fun1
end

function spline(f::Function, a::Float64, b::Float64, N::Int; bc="error")::Function
  spl = spline_spl(f, a, b, N, bc)

  fspl(x::Real)::Float64 = evaluate(spl, x)
end

function spline_derivative(f::Function, a::Float64, b::Float64, N::Int;
                           bc="error")::Function
  spl = spline_spl(f, a, b, N, bc)

  f′spl(x::Real)::Float64 = derivative(spl, x)
end

function spline_integral(f::Function, a::Float64, b::Float64, N::Int)::Float64
  spl = spline_spl(f, a, b, N, "zero")
  integrate(spl, a, b)
end

function spline_spl(f::Function, a::Float64, b::Float64, N::Int,
                    bc::AbstractString)::Spline1D
  @assert a < b "invalid arguments for spline"

  Δx::Float64 = (b - a) / N
  𝐱::Vector{Float64} = collect(a : Δx : b)
  𝐲::Vector{Float64} = [f(x) for x in 𝐱]

  Spline1D(𝐱, 𝐲, k=3, bc=bc)
end

function cardano(a₃::Real, a₂::Real, a₁::Real, a₀::Real)
  A₂ = a₂ / a₃
  A₁ = a₁ / a₃
  A₀ = a₀ / a₃

  p = (A₁ - 1/3*A₂^2) / 3
  q = (A₀ - 1/3*A₁*A₂ + 2/27*A₂^3) / 2

  D = - (q^2 + p^3)

  if D ≈ 0 && q ≈ 0
    y = [0]
  elseif D ≈ 0
    y = [-2*∛q, ∛q]
  elseif D > 0
    Θ = atan2(√D, -q)
    R = (q^2 + D)^(1/6) * cos(Θ/3)
    Q = (q^2 + D)^(1/6) * sin(Θ/3)
    y = [2R, -R - √3*Q, -R + √3*Q]
  elseif D < 0
    S = ∛(-q + √(-D))
    T = ∛(-q - √(-D))
    y = [S+T, -1/2*(S+T) + √3/2*im*(S-T), -1/2*(S+T) - √3/2*im*(S-T)]
  end

  x = y - A₂/3
end


text:
@testset "Hiring/Breaking Loops" begin
    a = 0
    @test begin
        @loopinvariant (a >= 0) for i = 1:10
            a += 1
        end
        a == 10
    end

    @test_throws ContractBreachException begin
        @loopinvariant (a >= 5) for i = 1:10
            a -= 1
        end
        a == 10
    end

    a = 0
    b = a
    @test begin
        @loopinvariant (a == 0) (b >= a) for i = 1:10
            b += 1
        end
        b == 10
    end

    a = 0
    @test begin
        @loopinvariant (a == 0) (b >= a) for i = 1:10
            b += 1
        end
        b != 10
    end

    a = 0
    b = a
    @test_throws ContractBreachException begin
        @loopinvariant (a == 0) (b >= a) for i = 1:10
            a += 1
        end
        b == 10
    end
end

@testset "Loop Breach Log" begin

    a = 0
    try
        @loopinvariant (a >= 5) for i = 1:10
            a -= 1
        end
    catch e
        b = IOBuffer()
        showerror(b, e)
        @test String(take!(b)) ==
              "ContractBreachException: Breach on Loop Invariant Expression 'a >= 5'"
    end

    try
        a = 0
        b = 100
        @loopinvariant (b >= a) (a == 0) for i = 1:10
            a += 1
        end
    catch e
        b = IOBuffer()
        showerror(b, e)
        @test String(take!(b)) ==
              "ContractBreachException: Breach on Loop Invariant Expression 'a == 0'"
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "GTK+"
version = v"3.24.10"

sources = [
    "https://download.gnome.org/sources/gtk+/3.24/gtk+-3.24.10.tar.xz" =>
    "35a8f107e2b90fda217f014c0c15cb20a6a66678f6fd7e36556d469372c01b03",
]

# Bash recipe for building across all platforms
# TODO: Theora and Opus once their releases are available
script = raw"""
cd $WORKSPACE/srcdir
cd gtk+-$(version)/
./configure --prefix=$prefix --extra-cflags="-I${prefix}/include" --extra-ldflags="-L${prefix}/lib"
make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, ["libgtk"], :gtk),
    LibraryProduct(prefix, ["libgdk"], :gdk),
    LibraryProduct(prefix, ["libgdk_pixbuf"], :gdk_pixbuf),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # The GLib library provides core non-graphical functionality such as high level data types, Unicode manipulation, and an object and type system to C programs.
    # Glib-related dependencies
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/PCRE-v8.42-2/build_PCRE.v8.42.0.jl",
    "https://github.com/giordano/Yggdrasil/releases/download/Libffi-v3.2.1-0/build_Libffi.v3.2.1.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Libiconv-v1.15-0/build_Libiconv.v1.15.0.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Gettext-v0.19.8-0/build_Gettext.v0.19.8.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Glib-v2.59.0%2B0/build_Glib.v2.59.0.jl",

    # The GdkPixbuf library provides facilities for loading images in a variety of file formats.
    # GdkPixBuf dependencies
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.4/build_Zlib.v1.2.11.jl",
    "https://github.com/giordano/Yggdrasil/releases/download/X11-v1.6.8/build_X11.v1.6.8.jl",
    "https://github.com/JuliaIO/LibpngBuilder/releases/download/v1.0.3/build_libpng.v1.6.37.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/JpegTurbo-v2.0.1-0/build_JpegTurbo.v2.0.1.jl",
    "https://github.com/JuliaIO/LibTIFFBuilder/releases/download/v6/build_libtiff.v4.0.9.jl",
    "https://github.com/ianshmean/SharedMimeInfoBuilder/releases/download/v1.0.0/build_shared-mime-info.v1.10.0.jl",

    # Pango is a library for internationalized text handling.
    # Pango-only dependencies
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/FriBidi-v1.0.5%2B0/build_FriBidi.v1.0.5.jl",
    "https://github.com/giordano/Yggdrasil/releases/download/Pango-v.1.42.4/build_Pango.v1.42.4.jl"

    # ATK is the Accessibility Toolkit. It provides a set of generic interfaces allowing accessibility technologies such as screen readers to interact with a graphical user interface.

    # Gobject Introspection is a framework for making introspection data available to language bindings.
    ### May be satisfied via the GLib build...

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: ianshmean/GTKBuilder

max_issues_repo_name: ianshmean/GTKBuilder

max_forks_repo_name: ianshmean/GTKBuilder


text:
using StructuresKit

#Test torsion calculation for a simply-supported beam with
#twist fixed warping free ends.
#ax=1.0 in., offset uniform load

#Ix Iy Ixy J Cw
sectionProperties = [(100.0, 100.0, 0.0, 10.0, 10.0)]

#E  ν
materialProperties = [(29500, 0.30)]

#ax ay
loadLocation = [(1.0, 0.0)]

#kx kϕ
springStiffness = [(0.0, 0.0)]

#member information
#L dL SectionProperties MaterialProperties LoadLocation SpringStiffness
memberDefinitions = [(30*12,3,1,1,1,1)]

#end boundary conditions
#type=1 u''=v''=ϕ''=0 (simply supported), type=2 u'=v'=ϕ'=0  (fixed), type=3 u''=v''=ϕ''=u'''=v'''=ϕ'''=0 (free end, e.g., a cantilever)
endBoundaryConditions = [1 1]

#supports
#location where u=v=ϕ=0
supports = [0.0 30*12]

#load  qx   qy
uniformLoad = (0.0, 5.0/1000)

z, u, v, ϕ, beamProperties = Beam.solve(memberDefinitions, sectionProperties, materialProperties, loadLocation, springStiffness, endBoundaryConditions, supports, uniformLoad)

T = InternalForces.torsion(z, beamProperties.dm, ϕ, beamProperties.E, beamProperties.G, beamProperties.J, beamProperties.Cw)

maxT = maximum(T)
maxT_theory = (uniformLoad[2] * beamProperties.ax[1] * z[end])/2

Error= abs((maxT - maxT_theory))/ maxT_theory

#accept 1% error from numerical solution
@test Error <= 0.01


text:
using ApproxFun, Base.Test
    import ApproxFun: ChebyshevDirichlet,Ultraspherical,space

@test_approx_eq Fun(exp,ChebyshevDirichlet{1,1})(.1) exp(.1)
@test_approx_eq Fun(Fun(exp,ChebyshevDirichlet{1,1}),Ultraspherical(1))(.1) exp(.1)

d=Interval()
sp=ChebyshevDirichlet{1,1}(d)
B=dirichlet(sp)
D=Derivative(d)
L=D^2+I

@test norm(([B;L]\[1.])-([dirichlet(d);L]\[1.])) <10eps()

f=Fun(t->cos(t)+cos(3t),CosSpace)

@test (f.*f-Fun(t->(cos(t)+cos(3t))^2,CosSpace)).coefficients|>norm <100eps()

f=Fun(exp,Taylor(Circle()))
g=Fun(z->1./(z-.1),Hardy{false}(Circle()))
@test_approx_eq (f(1.)+g(1.)) (exp(1.) + 1./(1-.1))

## Periodic
f=Fun(x->exp(-10sin((x-.1)/2)^2),Laurent)
@test_approx_eq f(.5) (Conversion(space(f),Fourier(domain(f)))*f)(.5)
@test_approx_eq f(.5) Fun(f,Fourier)(.5)

## PiecewiseSPace

x=Fun(identity,[-1.,0.,1.])
sp=space(x)
D=Derivative(sp)

u=[dirichlet(sp);
    D^2]\[1];
u2=[dirichlet();Derivative(Chebyshev())^2]\[1.]
@test_approx_eq u(0.) u2(0.)

x=Fun(identity,[-10.,0.,1.,15.])
sp=space(x)
D=Derivative(sp)

u=[dirichlet(sp);
    D^2-x]\[airyai(-10.)];

@test_approx_eq u(0.) airyai(0.)

s=Fun(sin,[-2.,2.])|>abs
c=Fun(cos,[-2.,2.])|>abs
sc=Fun(x->abs(sin(x))+abs(cos(x)),[-2,-π/2,0,π/2,2])
@test norm(sc-(c+s))<100eps()

# max/min creates breakpoints

x=Fun()
g=4*(x-.2)
f=max(-1,g)
f2=min(f,1)
f3=Fun(x->x<-0.05?-1.0:(x<0.45?4*(x-.2):1),[-1.0;-0.05;0.45;1.0])
@test norm(f2(collect(linspace(-1,1,10)))-f3(collect(linspace(-1,1,10)))) < 2eps()

x=Fun(identity,[im,0.,1.])
@test_approx_eq sqrt(1-x)(0.2im) sqrt(1-0.2im)
@test_approx_eq sqrt(1-x)(0.2) sqrt(1-0.2)

w=2/(sqrt(1-x)*sqrt(1+im*x))
@test_approx_eq w(0.2im) 2/(sqrt(1-0.2im)*sqrt(1+im*(0.2im)))
@test_approx_eq w(0.2) 2/(sqrt(1-0.2)*sqrt(1+im*(0.2)))

## Equivalent spaces

@test norm(Fun(cos,Chebyshev)-Fun(cos,Jacobi(-0.5,-0.5)))<100eps()
@test norm(Fun(cos,Chebyshev)-Fun(cos,JacobiWeight(0,0)))<100eps()
@test norm(Fun(cos,Jacobi(-0.5,-0.5))-Fun(cos,JacobiWeight(0,0))) < 100eps()
@test norm(Fun(cos,Chebyshev)-Fun(cos,JacobiWeight(0,0,Jacobi(-0.5,-0.5))))<100eps()
@test norm(Fun(cos,Jacobi(-0.5,-0.5))-Fun(cos,JacobiWeight(0,0,Jacobi(-0.5,-0.5))))<100eps()

## ContinuousSpace

import ApproxFun: PiecewiseInterval,ContinuousSpace

d=PiecewiseInterval(1.,2.,3.,4.)
S=ContinuousSpace(d)
D=Derivative(S)
u=[ldirichlet(S),D-I]\[exp(1.)]

@test_approx_eq u(1.1) exp(1.1)
@test_approx_eq u(3.4) exp(3.4)
@test_approx_eq last(u) exp(4)

d=PiecewiseInterval(0,1.,1.+im,im,0.)
@test_approx_eq Fun(exp,d)(.1) exp(.1)

## Triple SumSpace

x=Fun()
w=log(1-x)+sqrt(1-x^2)
f=w+x
@test_approx_eq f(0.1) (w(0.1)+0.1)
@test_approx_eq (w+1)(0.1) (w(0.1)+1)
@test_approx_eq (w+x+1)(0.1) (w(0.1)+1.1)
@test_approx_eq ((w+x)+1)(0.1) (w(0.1)+1.1)

## SumSpace bug

dsp=JacobiWeight(1.,0.,Jacobi(0.,1.,[0.,1.]))⊕JacobiWeight(0.5,0.,Jacobi(-0.5,0.5,[0.,1.]))
rsp=Legendre([0.,1.])⊕JacobiWeight(0.5,0.,Jacobi(0.5,0.5,[0.,1.]))
C=Conversion(dsp,rsp)
f=Fun([1.,2.,3.,4.,5.],dsp)
@test_approx_eq f(0.1) (C*f)(0.1)

## Piecewise + Cosntant

Γ=Circle()∪Circle(0.0,0.4)
G=Fun(z->in(z,Γ[2])?[1 0; -1/z 1]:[z 0; 0 1/z],Γ)   # Before the 80 wasn’t specified causing inconsistency
@test_approx_eq (G-I)(1.) (G(1.)-I)

## Previoius seffdault

x=Fun(identity,[-1.,1.])
f=x+sin(2x)*sqrt(1-x^2)
@test_approx_eq f(0.1) 0.1+sin(2*0.1)*sqrt(1-0.1^2)

## Check multiple piecewisesapce

x=Fun(identity,[-3,-2])+Fun(identity,[2,3])
w=sqrt(9-x^2)
f=w+Fun()
@test_approx_eq (f+w)(2.5) 2w(2.5)
@test_approx_eq (f+w)(.5) f(.5)

## Check Jacobi recurrence bug

S=Jacobi(-.5,.5)
f=Fun(exp,S)
@test_approx_eq f(0.1) exp(0.1)

## Check cancel conversion works
x=Fun([0.,1.])
f=exp(x)-1
Fun(f,JacobiWeight(1.,0.,[0.,1.]))

## Hermite
f=Fun(x->x+x^2,Hermite())
@test_approx_eq f(1.) 2.

## Arc exp

z=Fun(identity,Arc(0.,.1,0.,π/2))
@test_approx_eq exp(z)(0.1exp(0.2im)) exp(0.1exp(0.2im))


text:
# Helper Script that reads in the problems in the SDPA sparese file format and recreates and
# stores the problem data as matrices in Julia compatible JLD2 format

# The following is stored in the JLD2 file:
# m: dimension of x
# n: dimension of Fs,X
# nblocks: number of blocks in the diagonal structure of the matrices
# blockvec: vector of numbers that give the sizes on the individual blocks
# c: objective vector
# F: array containing the constraint matrices, F[1]=F0, F[2]=F1, F[m+1]=Fm
# optval: Optimal objective value (extracted from README file)

# by Michael Garstka
# University of Oxford, Control Group
# January 2018
# Updated to Julia 1.0 in April 2019

using FileIO, LinearAlgebra, SparseArrays
# Specify path to .dat-s files
dirPath = "./sdplib/"

# create array of all relevant filenames that contain problems
fileNames = []
for f in filter(x -> endswith(x, "dat-s"), readdir(dirPath))
    f = split(f,".")[1]
    push!(fileNames,f)
end

for file in fileNames
  # open file and read everything into memory
  f = open(dirPath*file*".dat-s");
  lines = readlines(f)
  close(f)

  # initialize variables
  counter = 1
  # size of x
  m = Inf
  # size of symmetric F matrices
  n = Inf
  nblocks = Inf
  blockvec = zeros(1)
  c = zeros(1)
  F = zeros(1)
  currentM = -1
  Fm = zeros(1)

  for ln in lines
    # dont do anything if line contains a comment
    if ln[1] == '"' || ln[1] == "*"
      counter == 1
    else
      # m: number of constraint matrices (in SDPA it starts from 0 -> +1)
      if counter == 1
        m = parse(Int64,strip(ln)) + 1
        # now that m is known, initalize the array that holds all the constraint matrices F
        F = Array{SparseMatrixCSC}(undef, m)

      # nblocks: number of blocks in the diagonal structure of the matrices
      elseif counter == 2
        nblocks = parse(Int64,strip(ln))

      # vector of numbers that give the sizes on the individual blocks
      # negative number indicates diagonal submatrix
      elseif counter == 3
      bvec = split(replace(strip(ln,['{','}','(',')']),"+" =>""))
       blockvec = [parse(Float64, ss) for ss in bvec]
       n = Int(sum(abs.(blockvec)))

       # objective function vector c
      elseif counter == 4
        # try first to split by whitespace as delimiter (also remove certain trailing, ending characters)
        cvec = split(replace(strip(ln,['{','}','(',')']),"+" => ""))
        # otherwise try comma as delimiter
        if length(cvec) == 1
          cvec = split(replace(strip(ln,['{','}','(',')']),"+" => ""),",")
        end
        c = [parse(Float64, ss) for ss in cvec]

      # all other lines contain constraint matrices with one entry per line
      # save them directly as sparse matrix
      else
        # FIXME: Accuracy
        line = [parse(Float64, ss) for ss in split(ln)]
        matno = Int(line[1])
        blkno = Int(line[2])
        i = Int(line[3])
        j = Int(line[4])
        entry = line[5]

        if matno > currentM
          # save previously generated matrix in sparse format
          if currentM >= 0
            F[currentM+1] = sparse(Fm)
          end
          # create new temporary matrix that is filled
          Fm = spzeros(n,n)
          currentM = matno
        end

        if blkno == 1
          iF = i
          jF = j
        else
          shift = Int(sum(abs.(blockvec[1:blkno-1])))
          iF = i+shift
          jF = j+shift
        end
        Fm[iF,jF] = entry
        # ensure symmetry
        if iF != jF
          Fm[jF,iF] = entry
        end
      end
      counter += 1
    end
  end
  # after last line has reached save the last F matrix
  F[currentM+1] = Fm

  # extract solution for objective value from README file
  f = open(dirPath*"README")
  lines = readlines(f)
  close(f)
  lines = lines[20:113]
  ln = filter(s -> occursin(file, s) ,lines)
  if length(ln) == 0
    optVal = "Not provided"
  else
    str = split(ln[1])[4]
    if occursin("primal", str)
      optVal = Inf
    elseif occursin("dual", str)
      optVal = -Inf
    else
    optVal = parse(Float64, str)
    end
  end

  # save to JLD file
  save(dirPath*file*".jld2", "m", m-1,"n",n,"nblocks",nblocks,"blockvec",blockvec,"c",c,"F",F,"optVal",optVal)
  println("Saved problem: $(file) to "*dirPath*file*".jld2, Optimal Value: $(optVal)")
end


text:
using Documenter, PowerModelsSecurityConstrained

makedocs(
    modules = [PowerModelsSecurityConstrained],
    sitename = "PowerModelsSecurityConstrained",
    authors = "Carleton Coffrin",
    pages = [
        "Home" => "index.md",
        "Solver Components" => "components.md",
    ]
)

deploydocs(
    repo = "github.com/lanl-ansi/PowerModelsSecurityConstrained.jl.git",
)


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

using ACE, JuLIP, Test
using ACE: evaluate
using JuLIP.Testing: print_tf

#---

@info("Construct basic and parameterised basis")

r0 = 2.3
rcut = 5.5
rin = 0.0
maxn = 5
species = [:X ]
D = SparsePSHDegree(wL = 1.0)

trans = ACE.PolyTransform(1, r0)
J = ACE.OrthPolys.transformed_jacobi(10, trans, rcut, rin)
P1 = ACE.RPI.PSH1pBasis(J, maxn, D=D, species = species)

basis = RPIBasis(P1, 3, D, maxn)

J5 = ACE.OrthPolys.transformed_jacobi(5, trans, rcut, rin)
P1basic = ACE.RPI.BasicPSH1pBasis(J5)
basic = RPIBasis(P1basic, 3, D, maxn)

#--- first test: make sure the bases are equivalent

@info("Test bases with and without parameters match")
for ntest = 1:30
   local R, Z, z0 = ACE.Random.rand_nhd(12, J, species)
   print_tf(@test ACE.evaluate(basis, R, Z, z0) ≈ ACE.evaluate(basic, R, Z, z0))
end
println()

#--- second test: perturb parameters

@info("Test basis with perturbed parameters doesn't match (duh...)")
params = basis.pibasis.basis1p.C[1]
basis.pibasis.basis1p.C[1] .+= 0.1 * (rand(size(params)...) .- 0.5)
for ntest = 1:30
   local R, Z, z0 = ACE.Random.rand_nhd(12, J, species)
   print_tf(@test !(ACE.evaluate(basis, R, Z, z0) ≈ ACE.evaluate(basic, R, Z, z0)))
end
println()


text:
#####
### seconds_analysis,jl
#####
# requires nothing, so long as PSG_01.jl was run.
# experiment with various factors
using Gadfly

include("DataFrameBinning.jl")
df = readtable("C:/Users/Peter.Mancini/Documents/Datafiles/playersessiongap_post.csv")

σ1 = 1
σ2 = 2
σ15 = 1.5

df[:idx] = find_bins(df[:similarity], range(0.0,0.01,101))
change_mid = sum(counts(sample(df[df[:similarity] .< 0.25,:idx], 10000), 0:25))/2
change_cumsum = cumsum(counts(sample(df[df[:similarity] .< 0.25,:idx], 10000), 0:25))
stay_mid = sum(counts(sample(df[df[:similarity] .> 0.90,:idx], 10000), 75:100))/2
stay_cumsum = cumsum(counts(sample(df[df[:similarity] .> 0.90,:idx], 10000), 75:100))

sim_stay = 90
sim_change= 20

Stay_gap = by(df[df[:similarity] .> sim_stay/100, :], :hour, x->median(x[:gap]))
stay_seconds_bydayhour = by(df[df[:similarity] .> sim_stay/100, :], [:day, :hour], x->quantile(x[:gap],[0.5]))

change_gap = by(df[df[:similarity] .< sim_change/100, :], :hour, x->median(x[:gap]))
change_seconds_bydayhour = by(df[df[:similarity] .< sim_change/100, :], [:day, :hour], x->quantile(x[:gap],[0.5]))

Ns_stay = by(df[df[:similarity] .> sim_stay/100,:], [:day, :hour], x-> size(x)[1])
stay_ymaxes = Array{Float64}(7,24)
stay_ymins = Array{Float64}(7,24)

P1= plot(stay_seconds_bydayhour, x="hour", y="x1", Geom.point, color="day",
  Guide.title(string("Durant c:", sim_change, " s:", sim_stay)),
  Guide.YLabel("Seconds Wait Time"),
  Guide.XLabel("By Hour"))
draw(PNG("Gap_stay_bydayhour.png", 8inch, 6inch), P1)


text:
gr()
default()
default(
    fontfamily="Arial",
    tick_direction=:out,
    guidefontsize=9,
    annotationfontfamily="Arial",
    annotationfontsize=10,
    annotationhalign=:left
)

alphabet = "abcdefghijklmnopqrstuvwxyz"

function add_plot_labels!(plt;offset=0)
    n = length(plt.subplots)
    for i = 1:n
        plot!(plt,subplot=i,title="($(alphabet[i+offset]))")
    end
    plot!(
        titlelocation = :left,
        titlefontsize = 10,
        titlefontfamily = "Helvetica Bold"
    )
end

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
require("CVMMean.jl")
require("CVMAddition.jl")
require("CVMVariance.jl")
require("CVMBigramAddition.jl")
require("CVMBigramMean.jl")
require("Corpora.jl")
require("CommonTools.jl")

using Variance
using BigramMean
using BigramAddition
using Addition
using Mean
using ArrayViews
using CVMAbstractModel
using Corpora
using CommonTools.CVMIntermediates
using Calculus

srand(300)

for (ModelModule, Model, Parameters) in [
        (Variance, Variance.VarianceModel, Variance.VarianceParameters),
        # The bigram based models will fail the gradient check
        # since we are using an approximate version of tanh for speed reasons
        #(BigramMean, BigramMean.BigramMeanModel, BigramMean.BigramMeanParameters),
        #(BigramAddition, BigramAddition.BigramAdditionModel, BigramAddition.BigramAdditionParameters),
        (Mean, Mean.MeanModel, Mean.MeanParameters),
        (Addition, Addition.AdditionModel, Addition.AdditionParameters)]
    println("ModelModule: ", ModelModule,
                " Model: ", Model, " Parameters: ", Parameters)

    model = Model(5, 5.0, 0.0, 0.0, 0.0, 0, 1.0, 10, 1)

    function init_params(nvocab::Int64, params::AbstractParameters)
        params.wordvecs = randn(model.wordvecdim, nvocab) * 0.1
        params.gradients = zeros(model.wordvecdim, nvocab)
    end

    function init_data_token(word, iword, params::AbstractParameters)
        return DataToken(word, iword, [1],
                            view(params.wordvecs, :, iword),
                            view(params.gradients, :, iword))
    end

    intermediates = [CVMIntermediates(model.wordvecdim, 25),
                    CVMIntermediates(model.wordvecdim, 25),
                    CVMIntermediates(model.wordvecdim, 25)]

    params1 = Parameters()
    params2 = Parameters()

    init_params(20, params1)
    init_params(19, params2)

    tokens1 = Token[init_data_token("$i",i,params1) for i=1:20]
    tokens2 = Token[init_data_token("$i",i,params2) for i=1:19]

    s_bi_1 = Token[tokens1[rand(1:20)] for i=1:15]
    s_mono_outer_1 = Token[tokens1[rand(1:20)] for i=1:11]
    s_mono_inner_1 = s_mono_outer_1[3:8]
    s_bi_2 = Token[tokens2[rand(1:19)] for i=1:17]
    s_mono_outer_2 = Token[tokens2[rand(1:19)] for i=1:13]
    s_mono_inner_2 = s_mono_outer_2[3:8]

    get_rand(range) = 5

    function set_copy!(a, b)
        @assert length(a) == length(b)
        @simd for i=1:length(a)
            b[i] = a[i]
        end
    end

    # Accumulate gradients (the analytical way)
    accumulate_gradients_batch!(model, s_mono_outer_1,
                        s_mono_inner_1, Vector{Token}[s_bi_1], intermediates,
                            scaling=1.0)

    function joggle_words(wordvecs)
        reshape(wordvecs, size(params1.wordvecs))
        tmp = deepcopy(params1.wordvecs)
        set_copy!(wordvecs, params1.wordvecs)
        err = objective_gradients_batch(model, s_mono_outer_1,
                        s_mono_inner_1, Vector{Token}[s_bi_1], intermediates,
                            scaling=1.0)
        set_copy!(tmp, params1.wordvecs)
        return err
    end

    function calc_fd_words()
        return Calculus.finite_difference(joggle_words,
                    reshape(deepcopy(params1.wordvecs),length(params1.wordvecs)));
    end

    ana_words = reshape(params1.gradients,length(params1.gradients))
    fd_words = calc_fd_words()
    fd_words_shaped = reshape(fd_words, size(params1.gradients))
    @test_approx_eq_eps ana_words fd_words 1e-6
end

println("All tests were successful")

max_stars_repo_name: ogh/binclusion

max_issues_repo_name: ogh/binclusion

max_forks_repo_name: ogh/binclusion


text:
using Test
using KernelFunctions
using Random: seed!

seed!(42)

dims = (10,5)
X = rand(dims...)
x = rand(dims[1])
s = 3.0
v1 = vcat(3.0,4.0*ones(dims[2]-1))
v2 = vcat(3.0,4.0*ones(dims[1]-1))
P = rand(5,10)
sdims = [1,2,3]
f(x) = sin.(x)

@testset "Transform Test" begin
    ## Test Scale Transform
    @testset "IdentityTransform" begin
        @test KernelFunctions.apply(IdentityTransform(),X)==X
    end
    @testset "ScaleTransform" begin
        t = ScaleTransform(s)
        @test all(KernelFunctions.apply(t,X).==s*X)
        s2 = 2.0
        KernelFunctions.set!(t,s2)
        @test all(t.s.==[s2])
        @test isequal(ScaleTransform(s),ScaleTransform(s))
    end
    ## Test ARD Transform
    @testset "ARDTransform" begin
        vt1 = ARDTransform(v1)
        vt2 = ARDTransform(v2)
        @test all(KernelFunctions.apply(vt1,X,obsdim=1).==v1'.*X)
        @test all(KernelFunctions.apply(vt2,X,obsdim=2).==v2.*X)
        newv1 = rand(5)
        KernelFunctions.set!(vt1,newv1)
        @test all(vt1.v .== newv1)
        @test ARDTransform(s,dims[2]).v == ARDTransform(s*ones(dims[2])).v
        @test_throws DimensionMismatch KernelFunctions.apply(vt1,rand(3,4))
    end
    ## Test LowRankTransform
    @testset "LowRankTransform" begin
        tp = LowRankTransform(P)
        @test all(KernelFunctions.apply(tp,X,obsdim=2).==P*X)
        @test all(KernelFunctions.apply(tp,x).==P*x)
        @test tp.proj == P
        P2 = rand(5,10)
        KernelFunctions.set!(tp,P2)
        @test all(tp.proj.==P2)
        @test_throws AssertionError KernelFunctions.set!(tp,rand(6,10))
        @test_throws DimensionMismatch KernelFunctions.apply(tp,rand(11,3))
    end
    ## Test FunctionTransform
    @testset "FunctionTransform" begin
        tf = FunctionTransform(f)
        KernelFunctions.apply(tf,X,obsdim=1)
        @test all(KernelFunctions.apply(tf,X,obsdim=1).==f(X))
    end
    ## Test SelectTransform
    @testset "SelectTransform" begin
        ts = SelectTransform(sdims)
        @test all(KernelFunctions.apply(ts,X,obsdim=2).==X[sdims,:])
        @test all(KernelFunctions.apply(ts,x).==x[sdims])
        sdims2 = [2,3,5]
        KernelFunctions.set!(ts,sdims2)
        @test all(ts.select.==sdims2)
    end
    ## Test ChainTransform
    @testset "ChainTransform" begin
        t = ScaleTransform(s)
        tp = LowRankTransform(P)
        tf = FunctionTransform(f)
        tchain = ChainTransform([t,tp,tf])
        @test all(KernelFunctions.apply(tchain,X,obsdim=2).==f(P*(s*X)))
        @test all(KernelFunctions.apply(tchain,X,obsdim=2).==
                    KernelFunctions.apply(tf∘tp∘t,X,obsdim=2))
    end
end


text:
using Statistics

for (op, fun) in ((:sum, :sum), (:mean, :mean))
    isdefined(Statistics, op) || continue
    _op = Symbol("_", op)
    __op = Symbol("__", op)

    @eval begin
        Statistics.$op(x::AbstractJaxArray; dims=:, keepdims=_autokeep(dims)) = $_op(nothing, x, keepdims, dims)
        Statistics.$op(T::Type, x::AbstractJaxArray; dims=:, keepdims=_autokeep(dims)) = $_op(T, x, keepdims, dims)

        $_op(T, x, keepdims, dims...) = $__op(T, x, keepdims, _convert_dims(size(x), dims...))
        $__op(T, x, keepdims, dims) = numpy.$fun(x, axis=dims, dtype=T, keepdims=keepdims)
    end
end


text:
#==========================================================
X-gate for qubit #5 on IBM Casablanca
==========================================================#
using LinearAlgebra
using Plots
using FFTW
using DelimitedFiles
using Printf
using Ipopt
using Random

Base.show(io::IO, f::Float64) = @printf(io, "%20.13e", f)
pyplot()

using Juqbox

N = 2 # 3
Nguard = 1 # 0

Ntot = N + Nguard

casename = "x-gate" # base file name (used in optimize-once.jl)

# frequencies (in GHz, will be multiplied by 2*pi to get angular frequencies in the Hamiltonian matrix)
fa = 4.9639697
xa = 0.3215826
rot_freq = [fa] # Rotational frequencies

# Duration tailored for the IBM IQ mixer
dt_IQ = 2.0/9
Npulses = 160 # Must be evenly divisable by 16
T = dt_IQ * Npulses # Nano sec

utarget = Matrix{ComplexF64}(I, Ntot, N)
vtarget = Matrix{ComplexF64}(I, Ntot, N)

# unitary target matrix
if Nguard == 0
    utarget[1,1] = 0.0
    utarget[2,1] = 1.0
    utarget[1,2] = 1.0
    utarget[2,2] = 0.0
elseif Nguard == 1
    utarget[1,1] = 0.0
    utarget[2,1] = 1.0
    utarget[3,1] = 0.0
    utarget[1,2] = 1.0
    utarget[2,2] = 0.0
    utarget[3,2] = 0.0
end

omega1 = Juqbox.setup_rotmatrices([N], [Nguard], rot_freq)

# Compute Ra*utarget
rot1 = Diagonal(exp.(im*omega1*T))

# target in the rotating frame
vtarget = rot1*utarget

# target in lab frame
#vtarget = utarget

# setup ansatz for control functions
Nctrl = 1
Nosc = Nctrl
Nfreq = 1 # number of carrier frequencies

Random.seed!(2456)
# initial parameter guess

# setup carrier frequencies
om = zeros(Nosc,Nfreq)
# Note: same frequencies for each p(t) (x-drive) and q(t) (y-drive)
println("Carrier frequencies [GHz]: ", om[:,:]./(2*pi))

# setup drift Hamiltonian
number = Diagonal(collect(0:Ntot-1))

H0 = -0.5*(2*pi)*xa* (number*number - number)
println("Drift Hamiltonian/2*pi: ", H0)

# lowering matrix
amat = Bidiagonal(zeros(Ntot),sqrt.(collect(1:Ntot-1)),:U) # standard lowering operator matrix
# raising matrix
adag = transpose(amat) # raising operator matrix

# package the lowering and raising matrices together into an one-dimensional array of two-dimensional arrays
Hsym_ops=[Array(amat + adag)]
Hanti_ops=[Array(amat - adag)]
H0 = Array(H0)

# max parameter amplitude
maxpar = 4.0*(2*pi/T)/Nfreq # Bigger amplitude than for a constant pulse

maxamp = zeros(Nctrl, Nfreq)
maxamp .= maxpar # Same parameter bounds for all ctrl Hamiltonians and frequencies

# Estimate time step
Pmin = 80
nsteps = calculate_timestep(T, H0, Hsym_ops, Hanti_ops, [maxpar], Pmin)

println("Final time T = ", T, " # time steps per min-period, P = ", Pmin, " # time steps: ", nsteps)

# Initial conditions
Ident = Matrix{Float64}(I, Ntot, Ntot)
U0 = Ident[1:Ntot,1:N]

# setup the simulation parameters
params = Juqbox.objparams([N], [Nguard], T, nsteps, Uinit=U0, Utarget=vtarget, Cfreq=om, Rfreq=rot_freq,
                          Hconst=H0, Hsym_ops=Hsym_ops, Hanti_ops=Hanti_ops)

# setup the initial parameter vector, either randomized or from file
startFromScratch = true # true
startFile = "rabi-pcof-opt-alpha-0.5.dat"

if startFromScratch
    D1 = 5 # Number of B-spline coefficients per frequency, sin/cos
    nCoeff = 2*Nosc*Nfreq*D1 # factor '2' is for sin/cos
    pcof0  = maxpar*0.05.*ones(nCoeff)

    println("*** Starting from constant pcof with amplitude ", maxpar*0.05)
else
    # the data on the startfile must be consistent with the setup!
    # use if you want to have initial coefficients read from file
    pcof0 = vec(readdlm(startFile))
    nCoeff = length(pcof0)
    D1 = div(nCoeff, 2*Nosc*Nfreq) # factor '2' is for sin/cos
    nCoeff = 2*Nosc*Nfreq*D1 # just to be safe if the file doesn't contain the right number of elements
    println("*** Starting from B-spline coefficients in file: ", startFile)
end

# min and max coefficient values
useBarrier = true
minCoeff, maxCoeff = assign_thresholds_ctrl_freq(params, D1, maxamp)

zero_start_end!(params, D1, minCoeff, maxCoeff)

#-maxpar*ones(nCoeff);
#maxCoeff = maxpar*ones(nCoeff);

# For ipopt
maxIter = 150 # optional argument
lbfgsMax = 250 # optional argument

println("*** Settings ***")
println("System Hamiltonian coefficients: (fa, xa) =  ", fa, xa)
println("Total number of states, Ntot = ", Ntot, " Total number of guard states, Nguard = ", Nguard)
println("Using B-spline basis functions with carrier wave, # freq = ", Nfreq)
println("Number of coefficients per spline = ", D1, " Total number of parameters = ", nCoeff)
println("Max parameter amplitudes: maxpar = ", maxpar)
println("Tikhonov coefficients: tik0 = ", params.tik0)

# Allocate all working arrays
wa = Juqbox.Working_Arrays(params, nCoeff)
prob = Juqbox.setup_ipopt_problem(params, wa, nCoeff, minCoeff, maxCoeff, maxIter, lbfgsMax, startFromScratch)

# uncomment to run the gradient checker for the initial pcof
# addOption( prob, "derivative_test", "first-order"); # for testing the gradient

println("Initial coefficient vector stored in 'pcof0'")


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module NotebooksUtils

# import Pluto
import PlutoUI
# import Markdown
import Format

import Makie
import Makie.AbstractPlotting
import Makie.AbstractPlotting.MakieLayout
import WGLMakie
import GLMakie

mutable struct Defs
    authors::String     # authors
    desc::String        # description of the notebook
    width::Int64
    html::String        # work area to accumulate html syntaxt in order to save cell space (ouput at the end of a single cell)
end
Defs() = Defs("?Authors?")
Defs(authors) = Defs(authors, "")
Defs(authors, desc) = Defs(authors, desc, 1000, "")

function fix_path(path)
    return replace(abspath(path), "\\"=>"/")
end

"""
    function run(notebook_filename)

Launch Pluto and allow teh user to open a specific notebook.
"""
function run(notebook_filename)
    run(; path=notebook_filename)
end

"""
    function run(; port=nothing, path=nothing, sysimage_file=nothing )

Launch Pluto and allow teh user to open a specific notebook.
Also allow the usage of a sysimage file for faster loading.
"""
function run(; port=nothing, path=nothing, sysimage_file=nothing, auto_detect_sysimage = false )

    @eval begin
        import Pluto

        local file_path = $path
        local sysimage_file_path = $sysimage_file

        if (file_path === nothing)
            @info "Launching Pluto"
        else
            file_path = fix_path(file_path)
            if (!isfile(file_path))
                @warn "Can't find the notebook file [$(file_path)] - trying to locate it in samples"
                base_filename = splitdir(file_path)[2]
                file_path = fix_path(joinpath(splitdir(splitdir(@__DIR__)[1])[1], "samples", "notebooks", base_filename))
                if (!isfile(file_path))
                    @warn "Can't find the notebook file [$(file_path)] in the samples folder - launching Pluto without a file"
                else
                    @info "Launching Notebook [$(file_path)]"
                end
            else
                @info "Launching Notebook [$(file_path)]"
            end
        end

        # handling the sysimage file
        if (sysimage_file_path !== nothing)
            if (!isfile(sysimage_file_path))
                @warn "Can't find the sysimage file [$(sysimage_file_path)] - launching Pluto without a sysimage"
                sysimage_file_path = nothing
            end
        else
            if ($auto_detect_sysimage)
                if Sys.iswindows()
                    if (!isfile("JuliaSysimage.dll"))
                        sysimage_file_path = "JuliaSysimage.dll"
                    end
                elseif Sys.islinux()
                    if (!isfile("JuliaSysimage.so"))
                        sysimage_file_path = "JuliaSysimage.so"
                    end
                end
            end
        end

        local options = Pluto.Configuration.from_flat_kwargs(
            port=$port,
            notebook=file_path,
            launch_browser=true,
            host="127.0.0.1",
            require_secret_for_open_links = false,
            require_secret_for_access = false,
            run_notebook_on_load = true,
            sysimage=sysimage_file_path,
            # banner = "yes",
            )

        session = Pluto.ServerSession(options=options)

        Pluto.run(session)
    end
end

"""
    run_sample(sample_name::String)

Launch Pluto and allow the user to open a specific sample notebook. If a notebook of the same name exists in the current working folder,
it will be opened in Pluto, otherwise, the original sample notebook will be copied to the current folder and be used.
This beheviour will prevent users from updating the original sample notebook.
"""
function run_sample(sample_name::String, edit_original::Bool = false)
    folder, basename = splitdir(sample_name)
    filename = "DummyFile,jl"
    if isempty(folder)
        if (isfile(basename))
            # file already exists in current folder - launch pluto and point to it
            filename = fix_path(abspath(basename))
        else
            # file does not exists in current directory
            file_in_samples = fix_path(joinpath(splitdir(splitdir(@__DIR__)[1])[1], "samples", "notebooks", basename))
            if (isfile(file_in_samples))
                if (edit_original)
                    @warn "You are about to edit the original sample notebook"
                    filename = file_in_samples
                else
                    # file exists - need to be copied into the current folder
                    src_fn = file_in_samples
                    dst_fn = fix_path(abspath(basename))
                    @info "Copying the sample notebook [$basename] from the samples directory to the current one\n$src_fn ->\n$dst_fn"
                    cp(src_fn, dst_fn)
                    filename = dst_fn
                end
            else
                # file does not exist in samples - need to issue an error
                # @error "Can't find a sample notebook named [file_in_samples] - please check again"
                throw(ArgumentError("Can't find a sample notebook named [$file_in_samples] - please check again" ))
            end
        end
    end

    @info "Running sample notebook from [$filename]"
    run(; path=filename)
end

mutable struct VariableInfo
    bond::Any
    html::String
end
function GetVarInfo(bond)
    res = VariableInfo(bond, HTMLFromObj(bond))
    return res
end

struct UISlider
    range::AbstractRange
    default::Number
    dec::Int16
end
UISlider() = UISlider(1:10, 1, -1)
UISlider(r) = UISlider(r, r.start, -1)
UISlider(r, def) = UISlider(r, def, -1)
# UISlider(r, def, format::String) = UISlider(r, def, format)

function Base.show(io::IO, ::MIME"text/html", slider::UISlider)
    if (slider.dec == -1)
        print(io, """
            <input type="range"
            min="$(first(slider.range))"
            step="$(step(slider.range))"
            max="$(last(slider.range))"
            value="$(slider.default)"
            oninput="this.nextElementSibling.value=this.value">
            <output>$(slider.default)</output>""")
    else
        fmt = "%.$(slider.dec)f"
        print(io, """
            <input type="range"
            min="$(first(slider.range))"
            step="$(step(slider.range))"
            max="$(last(slider.range))"
            value="$(slider.default)"
            oninput="this.nextElementSibling.value=parseFloat(this.value).toFixed($(slider.dec))">
            <output>$(Format.cfmt( fmt, slider.default ))</output>""")
    end
end

function SetDefs(defs::Defs)

    @info "I'm in SetDefs"
    ret = PlutoUI.Show(MIME"text/html"(), """

            """)

    #@info "The return type $(typeof(ret))"

    return ret
end

function DefsClearHTML(defs::Defs)
    defs.html = ""
end

function DefsAddHTML(defs::Defs, html::String)
    defs.html = defs.html * html
end

function DefsHTML(defs::Defs)
    return defs.html
end

# we don't need to use this function anymore - we can extract the html of items and use it using the PlutoUI.show command
# function SetHTMLMarkdown(val::Bool)
#     if (val)
#         @info "SetHTMLMarkdown: Allow Raw HTML Tags in Markdown"
#         @eval Markdown.htmlesc(io::IO, s::AbstractString) = print(io,s)
#     else
#         @info "SetHTMLMarkdown: Original process - Raw HTML is not allowed"
#         @eval Markdown.htmlesc(io::IO, s::AbstractString) =
#                 for ch in s
#                 print(io, get(Markdown._htmlescape_chars, ch, ch))
#             end
#     end
# end

"""
    function SetBackend(defs::Defs, be::String)

    this is my first comment try
"""
function SetBackend(defs::Defs, be::String)

    if (be == "Web")
        @info "Makie backend set to WEB (WGLMakie)"
        WGLMakie.activate!()
        AbstractPlotting.__init__()
        AbstractPlotting.inline!(true)
    else
        @info "Makie backend set to STATIC (GLMakie)"
        GLMakie.activate!()
        AbstractPlotting.__init__()
        AbstractPlotting.inline!(true)
    end

end

"""
    InitNotebook(; port=8449)

initialize the JSServe package.
"""
function InitNotebook(; port=8449)
    @eval begin
        try
            import JSServe
            local port = 8449 # the port you want
            JSServe.JSSERVE_CONFIGURATION.listen_port[] = port
            JSServe.JSSERVE_CONFIGURATION.external_url[] = "http://localhost:$(port)"
            JSServe.JSSERVE_CONFIGURATION.content_delivery_url[] = "http://localhost:$(port)"
            return JSServe.Page() # needs to get displayed by Pluto
        catch e
            @warn "Can't initialize the JSServe package\n$e"
        end
    end
end

function HTMLFromObj(obj)
    io = IOBuffer()
	Base.show(io, MIME"text/html"(), obj)
    res = String(take!(io))

    return res
end

# function MDFromString(str)
#     # SetHTMLMarkdown(true)

# 	io = IOBuffer(str)
#     res_md = Markdown.parse(io, )

#     # SetHTMLMarkdown(false)

#     return res_md
# end

function HTMLFix(html::String)
    html = replace(html, "\r" => "")
    html = replace(html, "\n" => "")
    return html
end

function HTMLFloatingBox(items; name="plutoui-genericfloatingbox", header="?? header ??", kwargs...)

    res = ""

    res = res *
          """""" *
          """$header""" *
          """
""" *
          """"""

    for item in items
        item_level = 1
        if     (startswith(item, "@ "))
            item_level = 2
        elseif (startswith(item, "@@ "))
            item_level = 3
        elseif (startswith(item, "@@@ "))
            item_level = 4
        elseif (startswith(item, "@@@@ "))
            item_level = 5
        elseif (startswith(item, "@@@@@ "))
            item_level = 6
        elseif (startswith(item, "@@@@@@ "))
            item_level = 7
        end

        if (startswith(item, "@"))
            item2 = item[item_level+2:end]
        else
            item2 = item
        end

        # println("Item [$item]  Level [$item_level]")

        res = res *
              """
""" *
              """
""" *
              item2 *
              """
""" *
              """
""" *
              """\n"""

        #Input: $(@bind nnnn MySlider(1:100, 10))
    end

    # closing the nav tags
    res = res *
          """""" *
          """
""" *
          """"""

    # add the style for this floating box
    res = res * HTMLFloatingBoxStyle(name; kwargs...)

    return HTMLFix(res)
end

function HTMLFloatingBoxStyle(name::String; right="1rem", top="20rem", width="25%", kwargs...)

    @info "HTMLFloatingBoxStyle: right=$right, top=$top, width=$width"

    res = """"""

    return HTMLFix(res)
end

function HTMLNewDocLayout()

    res = """

    """

    return HTMLFix(res)
end

function HTMLFixTOC()

    res = """