text:
module VsI18n

using SimpleI18n

using VsRecorderBase: download_ocr_language, enum_to_string
using ..VsRecorder: datapath

export GlobalI18nContext, GameLanguage
export get_i18n, download_all_ocr_languages, all_ocr_languages,
    get_game_language

const GlobalI18nContext = Ref{Union{Nothing, SimpleI18n.I18nContext}}(nothing)

function __init__()
    GlobalI18nContext[] = SimpleI18n.setup(
        datapath("locales"),
        "en"
    )
end

@enum GameLanguage begin
    EN
    ZHS
    ZHT
    JA
    KO
    ES
    FR
    DE
    IT
end

const OCR_LANGUAGES = Dict(
    EN => "eng",
    ZHS => "chi_sim",
    ZHT => "chi_tra",
    JA => "jpn",
    ES => "spa",
    FR => "fra",
    DE => "deu",
    IT => "ita",
    KO => "kor",
)

function get_game_language(locale::AbstractString = SimpleI18n.get_system_language())
    lang = lowercase(locale)
    t = split(lang, '_', limit = 2)
    lang, region = if length(t) < 2
        t[1], ""
    else
        t
    end
    if lang ∈ ("zh", "chi")
        if region ∈ ("hk", "tw", "hant", "tra")
            ZHT
        else
            ZHS
        end
    elseif lang == "zht"
        ZHT
    elseif lang == "zhs"
        ZHS
    elseif lang ∈ ("ja", "jpn")
        JA
    elseif lang ∈ ("es", "spa")
        ES
    elseif lang ∈ ("fr", "fra")
        FR
    elseif lang ∈ ("de", "deu")
        DE
    elseif lang ∈ ("it", "ita")
        IT
    elseif lang ∈ ("ko", "kor")
        KO
    else  # defaultly fallback to EN
        EN
    end
end
get_game_language(x::GameLanguage) = x

get_code(lang::AbstractString) = SimpleI18n.parse_locale_name(lang)
get_code(lang::GameLanguage) = get_code(enum_to_string(lang))

function download_all_ocr_languages()
    for lang in values(OCR_LANGUAGES)
        download_ocr_language(lang) || return false
        if lang ∈ ("chi_sim", "chi_tra", "jpn")
            download_ocr_language(lang * "_vert") || return false
        end
    end
    true
end

function get_i18n(path...; lang = nothing, default = nothing)
    key = join(path, '.')
    i18n(GlobalI18nContext[], key, isnothing(default) ? key : default; language = lang)
end

end


text:
import PredictMD
import Test

if PredictMD.is_travis_ci_on_linux()
    Test.@test( PredictMD.is_filesystem_root("/") )
    Test.@test_throws(ErrorException, PredictMD.find_package_directory("/"))
end

ext: jl


text:
# map-reduce

const AbstractArrayOrBroadcasted = Union{AbstractArray,Broadcast.Broadcasted}

# GPUArrays' mapreduce methods build on `Base.mapreducedim!`, but with an additional
# argument `init` value to avoid eager initialization of `R` (if set to something).
mapreducedim!(f, op, R::AbstractGPUArray, A::AbstractArrayOrBroadcasted;
              init=nothing) = error("Not implemented") # COV_EXCL_LINE
# resolve ambiguities
Base.mapreducedim!(f, op, R::AbstractGPUArray, A::AbstractArray) = mapreducedim!(f, op, R, A)
Base.mapreducedim!(f, op, R::AbstractGPUArray, A::Broadcast.Broadcasted) = mapreducedim!(f, op, R, A)

neutral_element(op, T) =
    error("""GPUArrays.jl needs to know the neutral element for your operator `$op`.
             Please pass it as an explicit argument to (if possible), or register it
             globally your operator by defining `GPUArrays.neutral_element(::typeof($op), T)`.""")
neutral_element(::typeof(Base.:(|)), T) = zero(T)
neutral_element(::typeof(Base.:(+)), T) = zero(T)
neutral_element(::typeof(Base.add_sum), T) = zero(T)
neutral_element(::typeof(Base.:(&)), T) = one(T)
neutral_element(::typeof(Base.:(*)), T) = one(T)
neutral_element(::typeof(Base.mul_prod), T) = one(T)
neutral_element(::typeof(Base.min), T) = typemax(T)
neutral_element(::typeof(Base.max), T) = typemin(T)

# resolve ambiguities
Base.mapreduce(f, op, A::AbstractGPUArray, As::AbstractArrayOrBroadcasted...;
               dims=:, init=nothing) = _mapreduce(f, op, A, As...; dims=dims, init=init)
Base.mapreduce(f, op, A::Broadcast.Broadcasted{<:AbstractGPUArrayStyle}, As::AbstractArrayOrBroadcasted...;
               dims=:, init=nothing) = _mapreduce(f, op, A, As...; dims=dims, init=init)

function _mapreduce(f, op, As...; dims, init)
    # mapreduce should apply `f` like `map` does, consuming elements like iterators.
    bc = if allequal(size.(As)...)
        Broadcast.instantiate(Broadcast.broadcasted(f, As...))
    else
        # TODO: can we avoid the reshape + view?
        indices = LinearIndices.(As)
        common_length = minimum(length.(indices))
        Bs = map(As) do A
            view(reshape(A, length(A)), 1:common_length)
        end
        Broadcast.instantiate(Broadcast.broadcasted(f, Bs...))
    end

    # figure out the destination container type by looking at the initializer element,
    # or by relying on inference to reason through the map and reduce functions.
    if init === nothing
        ET = Broadcast.combine_eltypes(bc.f, bc.args)
        ET = Base.promote_op(op, ET, ET)
        (ET === Union{} || ET === Any) &&
            error("mapreduce cannot figure the output element type, please pass an explicit init value")

        init = neutral_element(op, ET)
    else
        ET = typeof(init)
    end

    sz = size(bc)
    red = ntuple(i->(dims==Colon() || i in dims) ? 1 : sz[i], length(sz))
    R = similar(bc, ET, red)
    mapreducedim!(identity, op, R, bc; init=init)

    if dims==Colon()
        @allowscalar R[]
    else
        R
    end
end

Base.any(A::AbstractGPUArray{Bool}) = mapreduce(identity, |, A)
Base.all(A::AbstractGPUArray{Bool}) = mapreduce(identity, &, A)

Base.any(f::Function, A::AbstractGPUArray) = mapreduce(f, |, A)
Base.all(f::Function, A::AbstractGPUArray) = mapreduce(f, &, A)
Base.count(pred::Function, A::AbstractGPUArray) = mapreduce(pred, +, A; init = 0)

Base.:(==)(A::AbstractGPUArray, B::AbstractGPUArray) = Bool(mapreduce(==, &, A, B))

# avoid calling into `initarray!`
Base.sum!(R::AbstractGPUArray, A::AbstractGPUArray) = Base.reducedim!(Base.add_sum, R, A)
Base.prod!(R::AbstractGPUArray, A::AbstractGPUArray) = Base.reducedim!(Base.mul_prod, R, A)
Base.maximum!(R::AbstractGPUArray, A::AbstractGPUArray) = Base.reducedim!(max, R, A)
Base.minimum!(R::AbstractGPUArray, A::AbstractGPUArray) = Base.reducedim!(min, R, A)

LinearAlgebra.ishermitian(A::AbstractGPUMatrix) = mapreduce(==, &, A, adjoint(A))


text:
module QJuliaPCG

using QJuliaInterface
using QJuliaEnums
using QJuliaBlas
using QJuliaReduce
using QJuliaSolvers

using LinearAlgebra
using Printf

using Plots

##########
# Reference: H. Van der Vorst, Q. Ye, "Residual replacement strategies for Krylov subspace iterative methods for the convergence of true residuals", 1999
##########

norm2    = QJuliaReduce.gnorm2
rdot     = QJuliaReduce.reDotProduct
cpy      = QJuliaBlas.cpy

@inline function MatPrecon(out::AbstractArray, inp::AbstractArray, outSloppy::AbstractArray, inpSloppy::AbstractArray, K::Function)

	if pointer_from_objref(out) == pointer_from_objref(inp); return; end #nothing to do

    outSloppy .=@. 0.0
    cpy(inpSloppy, inp)       #noop for the alias refs
	K(outSloppy, inpSloppy)
    cpy(out, outSloppy)       #noop for the alias refs

end

function solver(x::AbstractArray, b::AbstractArray, Mat::Any, MatSloppy::Any, param::QJuliaSolvers.QJuliaSolverParam_qj, K::Function, extra_args...)

    is_preconditioned = param.inv_type_precondition != QJuliaEnums.QJULIA_INVALID_INVERTER

    solver_name = is_preconditioned == false ? "CG" : "PCG"

    println("Running ", solver_name ," solver.")

    if is_preconditioned == true; println("Preconditioner  :: ", param.inv_type_precondition); end

    if (param.maxiter == 0)
      if param.use_init_guess == false; x .=@. 0.0; end
      return
    end #if param.maxiter == 0

    mixed = (param.dtype_sloppy != param.dtype)

    if mixed == true; println("Running mixed precision solver.");end

    local r    = Vector{param.dtype}(undef, length(x))
    local y    = Vector{param.dtype}(undef, length(x))
    # aux high precision vector
    local yaux = ones(param.dtype, length(y))
    local yaux2 = zeros(param.dtype, length(y))
    # sloppy residual vector
    local rSloppy    = mixed == true ? zeros(param.dtype_sloppy, length(b)) : r
    local rSloppyOld = zeros(param.dtype_sloppy, length(rSloppy))
    # search vector and Ap result
    local p       = zeros(param.dtype_sloppy, length(rSloppy))
    local s       = zeros(param.dtype_sloppy, length(rSloppy))
    # iterated sloppy solution vector
    local xSloppy = param.use_sloppy_partial_accumulator == true ? x : zeros(param.dtype_sloppy, length(rSloppy))
    # for the preconditioner
    local u       = is_preconditioned == true ? zeros(param.dtype_sloppy, length(rSloppy)) : rSloppy
    local rPre    = param.dtype_precondition != param.dtype_sloppy ? zeros(param.dtype_precondition, length(rSloppy)) : rSloppy
    local pPre    = param.dtype_precondition != param.dtype_sloppy ? zeros(param.dtype_precondition, length(rSloppy)) : u

    local plot_data = Array{Float64, 2}(undef, param.maxiter, 4)

    Precond(out, inp) = MatPrecon(out, inp, pPre, rPre,K)

    b2 = norm2(b)  #Save norm of b
    r2 = 0.0; r2_old = 0.0     #if zero source then we will exit immediately doing no work
    #
    ϵ     = eps(param.dtype_sloppy) / 2.0; ϵh = eps(param.dtype) / 2.0
    deps  = sqrt(ϵ); dfac = 1.1; nfact = 10.0*sqrt(Float64(length(r)))
    xnorm = 0.0; ppnorm = 0.0; Anorm = 0.0
	#
    y .=@. 1.0
    Mat(y, yaux)
    yaux .=@. abs.(y)
    Anorm = findmax(yaux)[1]
    # Relupdates parameters:
    println("Estimated matrix norm is ", Anorm)

    if param.use_init_guess == true
      #r = b - Ax0 <- real
      Mat(r, x)
      r .=@. b - r
      r2 = norm2(r)
      y .=@. x
    else
      r2 = b2;
      r .=@. b
      y .=@. 0.0
    end
	#
    xSloppy .=@. 0.0
    # Relupdates parameters:
    rUpdate = 0
    rNorm   = sqrt(r2)
    #dinit   = ϵh*(rNorm+nfact*Anorm*xnorm)
    dinit   = ϵh*(rNorm+(nfact+1)*Anorm*xnorm)
    dk      = dinit
    #
    cpy(rSloppy, r)
    #
    Precond(u, rSloppy)
    # initialize CG parameters
    β = 0.0; γold = 0
    γ = is_preconditioned == true ? rdot(rSloppy, u) : r2
    #iteration counters
    k = 0; converged = false

    # if invalid residual then convergence is set by iteration count only
    stop = b2*param.tol*param.tol
    println(solver_name," : Initial residual = ", sqrt(r2))
    println(solver_name," : Initial preconditioned residual ", norm(u))

    resIncrease = 0; resIncreaseTotal = 0;  relUpdates = 0

    updateR::Bool = false

    # Main loop:
    max_res_incr = 32
    aux_rel_trigger = max_res_incr
    while (k < param.maxiter && converged == false && aux_rel_trigger > 0 )
      # Update search vector
      p  .=@. u + β*p
      #
      MatSloppy(s, p)
      #
      γold = γ
      #
      pAp = rdot(p, s); ppnorm = norm2(p)
      #
      α = γ / pAp
      #
      xSloppy    .=@. xSloppy + α*p
	  #
      rSloppyOld .=@. rSloppy
      rSloppy    .=@. rSloppy - α*s
      rSloppyOld .=@. rSloppy - rSloppyOld

      # preconditioned residual
      Precond(u, rSloppy)

      r2_old = r2
      # compute remaining dot products
      r2   = norm2(rSloppy)
      γ    = rdot(u, rSloppy)
      γaux = rdot(u, rSloppyOld)
      γnew = γaux >= 0.0 ? γaux : γ

      β    = γnew / γold

      rNormOld = rNorm
      rNorm   = sqrt(r2)

      if(rNormOld < rNorm); aux_rel_trigger -= 1; end

      # xSloppy .= xSloppy + α*p <=> norm2(xSloppy) = norm2(xSloppy) + α*α*norm2(p)
      xnorm = xnorm + α*α*ppnorm
      dkm1  = dk
      #dk    = dkm1 + ϵ*rNorm+ϵh*abs(α)*Anorm*sqrt(xnorm)
      dk  = dkm1 + ϵ*sqrt(rNorm + (Anorm*sqrt(xnorm) + (nfact+4.0)*abs(α)*Anorm*sqrt(ppnorm)))

      updateR =  ( ((dkm1 <= deps*sqrt(r2_old)) && ((dk > (deps * rNorm)))) && (dk > dfac * dinit) )
	  #updateR = ( (dkm1 <= (deps*sqrt(r2_old))) && ((dk > (deps * rNorm))) )

      if updateR == true
        println("Do reliable update.")
        x .=@. xSloppy
        y .=@. x + y
        Mat(r, y)
        r .=@. b - r
        # Reset sloppy residual and solution vectors
        cpy(rSloppy, r)
        xSloppy .=@. 0.0
        # preconditioned residual
        Precond(u, rSloppy)
        #
        r2 = norm2(r)
		#
        γ  = rdot(rSloppy, u)
        # Reorthogonalize previous search direction against the residual vector
        rp = rdot(rSloppy, p) / r2
        p  .=@. p - rp*rSloppy
        # Recompute β after reliable update
        β  = γ / γold
		# Reset reliable parameters
        dinit = ϵh*(sqrt(r2) + Anorm*sqrt(norm2(y)))
        dk = dinit; xnorm = 0.0

        if sqrt(γ) > rNorm
          resIncrease      += 1
          resIncreaseTotal += 1
          println("Update residual is higher than iterative residual..", sqrt(γ), " the previous iter residual is ", rNorm)
        else
          resIncrease       = 0
        end

        relUpdates += 1
      end
      # Check convergence:
      converged = (γ > stop) ? false : true
      # Update iter index
      k += 1

      @printf("%s: %d iteration, iter residual: %le  \n", solver_name, k, sqrt(γ) / sqrt(b2))

    end #while

    x .= @. xSloppy
    y .= @. x + y

    if (param.compute_true_res == true)
      Mat(r, y)

      r .=@. b - r
      r2 = norm2(r)

      param.true_res = sqrt(r2 / b2)
      println(solver_name, ": converged after ", k , "  iterations, relative residual: true = ", sqrt(r2), " after reliable updates number ", relUpdates)

    end #if (param.compute_true_res == true)

    x .= @. y

end #solver

end #QJuliaPCG

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
"""
    precompile_activator(package_path)

Activates precompile of a package by setting `should_precompile = true`

`package_path = pathof_noload(package_name)`
"""
function precompile_activator(package_path::AbstractString)
    includer_file = "$(dirname(package_path))/precompile_includer.jl"
    if !isfile(includer_file)
        error("$includer_file doesn't exists")
    else
        file_text = Base.read(includer_file, String)
        regp = r"(should_precompile\s*=\s*)(true|false)"
        m = match(regp, file_text)
        if m !== nothing
            file_text = replace(file_text, regp => s"\1 true")
            Base.write(includer_file, file_text)
        else
            error("\"should_precompile = ...\" doesn't exist")
        end
    end
end

"""
    precompile_deactivator(package_path)

Deactivates precompile of a package by setting `should_precompile = false`

`package_path = pathof_noload(package_name)`
"""
function precompile_deactivator(package_path::AbstractString)
    includer_file = "$(dirname(package_path))/precompile_includer.jl"
    if !isfile(includer_file)
        error("$includer_file doesn't exists")
    else
        file_text = Base.read(includer_file, String)
        regp = r"(should_precompile\s*=\s*)(true|false)"
        m = match(regp, file_text)
        if m !== nothing
            file_text = replace(file_text, regp => s"\1 false")
            Base.write(includer_file, file_text)
        else
            error("\"should_precompile = ...\" doesn't exist")
        end
    end
end


text: not_implemented_for(x) = error("Not implemented for $x")


text:
export AbstractMStorage, MAX_PREGENERATED_MSTORAGE_SIZE

abstract type AbstractMStorage end

const MAX_PREGENERATED_MSTORAGE_SIZE = 16

import Base: length

"""
    @MStorage(n::Int)

Helper function to generate tuple-like structure MStorageN, but with mutable fields and empty constructor.
It is possible then to take a `snapshot(::MStorage)` which returns a tuple with the same types and values from storage.
Some operators and observables use pregenerated `MStorage` to instatiate uninitialized
mutable storage in case when stream is allowed to not emit any values before completion.

# Generated structure layout

```
struct MStorageN{V1, V2, ..., VN}
    v1 :: V1
    v2 :: V2
    ...
    vn :: VN
end
```

See also: [`setstorage!`](@ref)
"""
macro MStorage(n::Int)
    name   = Symbol(:MStorage, n)
    types  = map(i -> Symbol(:V, i), 1:n) # types  = [ V1, V2, ..., V2 ]
    fields = map(i -> Symbol(:v, i), 1:n) # fields = [ v1, v2, ..., vn ]

    # Generates structure with the following layout:
    #
    # struct MStorageN{V1, V2, ..., VN} <: AbstractMStorage
    #     v1 :: V1
    #     v2 :: V2
    #     ...
    #     vn :: VN
    # end
    structure = Expr(:struct, true,                                   # <-- true stands for mutable
        Expr(:(<:), Expr(:curly, name, types...), :AbstractMStorage), # $name{types...} <: AbstractMStorage
        Expr(:block,                                                  # begin
            map(z -> Expr(:(::), z[1], z[2]), zip(fields, types))..., # v1::V1, v2::V2, ..., vn::VN
            Expr(:(=),                                                # constructor() = new()
                Expr(:where, Expr(:call, Expr(:curly, name, types...)), types...),
                Expr(:block, Expr(:call, :new))
            )
        )
    )

    # Generates function `snapshot`
    #
    # snapshot(storage::$name) = (storage.v1, storage.v2, ..., storage.vn)
    snapshot = Expr(:(=),
        Expr(:call, :snapshot, Expr(:(::), :s, name)),
        Expr(:block, Expr(:tuple, map(f -> Expr(:(.), :s, QuoteNode(f)), fields)...))
    )

    # Generates function `length`
    #
    # length(storage::$name) = $n
    length = Expr(:(=), Expr(:call, :length, Expr(:(::), :s, name)), Expr(:block, n))

    # Generates function `mstorage`
    #
    # `mstorage(::Val{$n}, ::Type{ <: Tuple{V1, V2, ..., VN} }) = $name{V1, V2, ..., VN}()`
    mstorage = Expr(:(=),
        Expr(:where,
            Expr(:call, :mstorage,
                Expr(:(::), Expr(:curly, :Type, Expr(:curly, :Tuple, types...))),
                Expr(:(::), Expr(:curly, :Val, n))
            ),
            types...
        ),
        Expr(:block, Expr(:call, Expr(:curly, name, types...)))
    )

    # Generates function `setstorage!`
    #
    # function setstorage!(s::$name, v, index::Int)
    #     if index === 1
    #         s.v1 = v
    #         return nothing
    #     end
    #     ...
    # end
    setstorage = Expr(:function,
        Expr(:call, :setstorage!, Expr(:(::), :s, name), :v, Expr(:(::), :index, :Int)),
        Expr(:block,
            map(i -> Expr(:if,
                Expr(:call, :(===), :index, i),
                Expr(:block,
                    Expr(:(=), Expr(:(.), :s, QuoteNode(fields[i])), :v),
                    Expr(:return, :nothing)
                )
            ), 1:n)...
        )
    )

    output = quote
        $structure
        $snapshot
        $length
        $mstorage
        $setstorage
    end

    return esc(output)
end

"""
    setstorage!(s, v, ::Val{I}) where I

This function can be used to set a new value `v` for storage `s` with a given value `v` and index `I`.
Using parametrized `Val{I}` for indexing ensures for index to be resolved at compile-time and if-else branch optimization.

See also: [`@MStorage`](@ref)
"""
function setstorage!(s::S, v, ::Val{I}) where { S <: AbstractMStorage, I }
    if I === 1
        s.v1 = v
    elseif I === 2
        s.v2 = v
    elseif I === 3
        s.v3 = v
    elseif I === 4
        s.v4 = v
    elseif I === 5
        s.v5 = v
    elseif I === 6
        s.v6 = v
    elseif I === 7
        s.v7 = v
    elseif I === 8
        s.v8 = v
    elseif I === 9
        s.v9 = v
    elseif I === 10
        s.v10 = v
    elseif I === 11
        s.v11 = v
    elseif I === 12
        s.v12 = v
    elseif I === 13
        s.v13 = v
    elseif I === 14
        s.v14 = v
    elseif I === 15
        s.v15 = v
    elseif I === 16
        s.v16 = v
    end
end

mstorage(::Type{T}, ::Val{N}) where { N, T <: NTuple{N, Any} } = Vector{Any}(undef, N)
snapshot(s::S) where { S <: Vector } = tuple(s...)

setstorage!(s::S, v::T, ::Val{I}) where { T, S <: Vector{T}, I } = s[I] = v
setstorage!(s::S, v::T, I::Int)   where { T, S <: Vector{T} }    = s[I] = v

macro GenerateMStorages()
    output = quote end
    for i in 1:MAX_PREGENERATED_MSTORAGE_SIZE
        output = quote
            $output
            @MStorage($i)
        end
    end
    return esc(output)
end

@GenerateMStorages()

_staticlength(::Type{T}) where { N, T <: NTuple{N, Any} } = Val{N}()

getmstorage(::Type{T}) where T = mstorage(T, _staticlength(T))


text:
using Mimi

@defcomp emissions begin
    regions = Index() # regions index must be specified

    E = Variable(index = [time, regions]) # Total greenhouse gas emissions
    E_Global = Variable(index = [time]) # Global emissions (sum of regional emissions)
    sigma = Parameter(index = [time, regions]) # Emissions output ratio
    YGROSS = Parameter(index = [time, regions]) # Gross output

    function run_timestep(p,v,d,t)
        # Define an equation for E
        for r in d.regions
            v.E[t,r] = p.YGROSS[t,r] * p.sigma[t,r]
        end

        # Define an equation for E_Global
        for r in d.regions
            v.E_Global[t] = sum(v.E[t,:])
        end
    end
end

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
export to_basictypes

"""
    to_basictypes(block::AbstractBlock{N}) where N

convert gates to basic types

    * ChainBlock
    * PutBlock
    * PrimitiveBlock
"""
function to_basictypes end

to_basictypes(block::PrimitiveBlock) = block
function to_basictypes(block::AbstractBlock{N}) where {N}
    throw(NotImplementedError(:to_basictypes, typeof(block)))
end

function to_basictypes(block::RepeatedBlock{N}) where {N}
    chain(N, map(i -> put(N, i => content(block)), block.locs))
end

to_basictypes(block::CachedBlock) = content(block)
function to_basictypes(block::Concentrator{N,<:PrimitiveBlock}) where {N}
    put(N, block.locs => content(block))
end
function to_basictypes(block::Concentrator{N}) where {N}
    to_basictypes(map_address(content(block), AddressInfo(N, [block.locs...])))
end
function to_basictypes(block::Concentrator{N,<:Measure}) where {N}
    map_address(content(block), AddressInfo(N, [block.locs...]))
end
to_basictypes(block::Daggered) where {N} = Daggered(block.content)
to_basictypes(block::Scale) where {N} = Scale(block.alpha, block.content)
to_basictypes(block::KronBlock{N}) where {N} = chain(N, [put(N, i => block[i]) for i in block.locs])
to_basictypes(block::Union{Add,PutBlock,ChainBlock,ControlBlock}) where {N} = block


text:
abstract type Integrator{DT} end

struct ExplicitEuler{DT, ET <: Equation{DT}} <: Integrator{DT}
    equ::ET
    Δt::DT
    ẋ::Vector{DT}

    function ExplicitEuler(equ::ET, Δt::DT) where {DT, ET <: Equation{DT}}
        new{DT,ET}(equ, Δt, zeros(DT, ndims(equ)))
    end
end

function (int::ExplicitEuler)(x₀::AbstractVector, x₁::AbstractVector)
    int.equ.f(int.ẋ, x₀)
    x₁ .= x₀ .+ int.Δt .* int.ẋ
end

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# StochasticPowerModels.jl                                                     #
# An extention package of PowerModels(Distribution).jl for Stochastic (Optimal)#
# Power Flow                                                                   #
# See http://github.com/timmyfaraday/StochasticPowerModels.jl                  #
################################################################################

module StochasticPowerModels

    # import pkgs
    import InfrastructureModels
    import Ipopt
    import JuMP
    import KernelDensity
    import Memento
    import PolyChaos
    import PowerModels

    # import types
    import PowerModels: AbstractPowerModel, AbstractACRModel, AbstractIVRModel

    # pkgs const
    const _IM = InfrastructureModels
    const _KDE = KernelDensity
    const _PCE = PolyChaos
    const _PM = PowerModels
    const _SPM = StochasticPowerModels

    # memento logger
    function __init__()
        global _LOGGER = Memento.getlogger(PowerModels)
    end

    # const
    const nw_id_default = 1

    # funct
    sorted_nw_ids(pm) = sort(collect(_PM.nw_ids(pm)))

    # paths
    const BASE_DIR = dirname(@__DIR__)

    # include
    include("core/constraint_template.jl")
    include("core/objective.jl")
    include("core/variable.jl")

    include("form/acr.jl")
    include("form/iv.jl")

    include("prob/sopf_acr.jl")
    include("prob/sopf_iv.jl")

    include("util/util.jl")

    # export
    export BASE_DIR

    export run_sopf_iv, run_sopf_acr

    export build_stochastic_data
    export pce_coeff, sample, density, print_summary
end


text:
"""
"""
struct LocalToGlobalPosNegArray{T,M,N,L,V,W} <: AbstractArray{Array{T,M},N}
  lid_to_gid::L
  gid_to_val_pos::V
  gid_to_val_neg::W

  @doc """
  """
  function LocalToGlobalPosNegArray(
    lid_to_gid::AbstractArray{<:AbstractArray},
    gid_to_val_pos::AbstractArray,
    gid_to_val_neg::AbstractArray)

    L = typeof(lid_to_gid)
    V = typeof(gid_to_val_pos)
    W = typeof(gid_to_val_neg)
    T = eltype(V)
    M = ndims(eltype(L))
    N = ndims(L)
    new{T,M,N,L,V,W}(lid_to_gid,gid_to_val_pos,gid_to_val_neg)
  end
end

size(a::LocalToGlobalPosNegArray) = size(a.lid_to_gid)

function IndexStyle(::Type{<:LocalToGlobalPosNegArray{T,M,N,L}}) where {T,M,N,L}
  IndexStyle(L)
end

function array_cache(a::LocalToGlobalPosNegArray)
  gids = testitem(a.lid_to_gid)
  T = eltype(a.gid_to_val_pos)
  r = zeros(T,size(gids))
  c = CachedArray(r)
  cl = array_cache(a.lid_to_gid)
  (cl,c)
end

function getindex!(cache,a::LocalToGlobalPosNegArray,i::Integer...)
  (cl,c) = cache
  gids = getindex!(cl,a.lid_to_gid,i...)
  setsize!(c,size(gids))
  r = c.array
  for (i,gid) in enumerate(gids)
    if gid > 0
      r[i] = a.gid_to_val_pos[gid]
    elseif gid < 0
      r[i] = a.gid_to_val_neg[-gid]
    else
      @unreachable "Only positive or negative indices allowed, not zero."
    end
  end
  r
end

function getindex(a::LocalToGlobalPosNegArray,i::Integer...)
  cache = array_cache(a)
  getindex!(cache,a,i...)
end


text:
@testset "Stats (TF,TF-IDF,BM25)" begin

    doc1 = "a a a sample text text"
    doc2 = "another example example text text"
    doc3 = ""
    doc4 = "another another text text text text"

    crps = Corpus(map(StringDocument, [doc1, doc2, doc3, doc4]))

    update_lexicon!(crps)
    m = DocumentTermMatrix(crps)

    T = StringAnalysis.DEFAULT_FLOAT_TYPE

    approx_eq(m1::AbstractMatrix{T}, m2::AbstractMatrix{T}; tol=1e-6
             ) where T = begin
        t = all(abs.(m1-m2) .<= tol)
        return t
    end

    max_tol = 1e-5
    # TF
    # Terms are in alphabetical ordering
    # (correctweights uses older rows are documents format, hence the transpose)
    correctweights = T.([0.408248  0.0       0.0  0.0;
                         0.57735   0.632456  0.0  0.816497;
                         0.707107  0.0       0.0  0.0;
                         0.0       0.447214  0.0  0.57735;
                         0.0       0.632456  0.0  0.0])

    myweights = tf(m)
    @test approx_eq(myweights, correctweights, tol=max_tol)

    myweights = tf(dtm(m))
    @test approx_eq(myweights, correctweights, tol=max_tol)
    @test typeof(myweights) <: SparseMatrixCSC

    myweights = tf(Matrix(dtm(m)))
    @test isnan(sum(myweights)) == 0
    @test approx_eq(myweights, correctweights, tol=max_tol)
    @test typeof(myweights) <: Matrix

    myweights = float(dtm(crps));
    tf!(myweights)
    @test approx_eq(T.(myweights), correctweights, tol=max_tol)
    @test typeof(myweights) <: SparseMatrixCSC
    @test eltype(myweights) == typeof(1.0)

    myweights = float(Matrix(dtm(crps)));
    tf!(myweights)
    @test approx_eq(T.(myweights), correctweights, tol=max_tol)
    @test typeof(myweights) <: Matrix
    @test eltype(myweights) == typeof(1.0)

    @test_throws MethodError tf!(DocumentTermMatrix(crps))
    @test tf!(DocumentTermMatrix{T}(crps)) isa SparseMatrixCSC{T, Int}

    # TF-IDF
    # Terms are in alphabetical ordering
    # (correctweights uses older rows are documents format, hence the transpose)
    correctweights = T.([0.691224  0.0       0.0  0.0;
                         0.57735   0.632456  0.0  0.816497;
                         1.19724   0.0       0.0  0.0;
                         0.0       0.575869  0.0  0.743444;
                         0.0       1.07084   0.0  0.0])
    myweights = tf_idf(m)
    @test approx_eq(myweights, correctweights, tol=max_tol)

    myweights = tf_idf(dtm(m))
    @test approx_eq(myweights, correctweights, tol=max_tol)
    @test typeof(myweights) <: SparseMatrixCSC

    myweights = tf_idf(Matrix(dtm(m)))
    @test isnan(sum(myweights)) == 0
    @test approx_eq(myweights, correctweights, tol=max_tol)
    @test typeof(myweights) <: Matrix

    myweights = float(dtm(crps));
    tf_idf!(myweights)
    @test approx_eq(T.(myweights), correctweights, tol=max_tol)
    @test typeof(myweights) <: SparseMatrixCSC
    @test eltype(myweights) == typeof(1.0)

    myweights = float(Matrix(dtm(crps)));
    tf_idf!(myweights)
    @test approx_eq(T.(myweights), correctweights, tol=max_tol)
    @test typeof(myweights) <: Matrix
    @test eltype(myweights) == typeof(1.0)

    @test_throws MethodError tf_idf!(DocumentTermMatrix(crps))
    @test tf_idf!(DocumentTermMatrix{T}(crps)) isa SparseMatrixCSC{T, Int}

    # Terms are in alphabetical ordering
    # (correctweights uses older rows are documents format, hence the transpose)
    correctweights = T.([0.685309  0.0       0.0  0.0;
                         0.542113  0.654905  0.0  0.713275;
                         1.08029   0.0       0.0  0.0;
                         0.0       0.637042  0.0  0.69807;
                         0.0       1.10885   0.0  0.0])
    myweights = bm_25(m)
    @test approx_eq(myweights, correctweights, tol=max_tol)

    myweights = bm_25(dtm(m))
    @test approx_eq(myweights, correctweights, tol=max_tol)
    @test typeof(myweights) <: SparseMatrixCSC

    myweights = bm_25(Matrix(dtm(m)))
    @test isnan(sum(myweights)) == 0
    @test approx_eq(myweights, correctweights, tol=max_tol)
    @test typeof(myweights) <: Matrix

    myweights = float(dtm(crps));
    bm_25!(myweights)
    @test approx_eq(T.(myweights), correctweights, tol=max_tol)
    @test typeof(myweights) <: SparseMatrixCSC
    @test eltype(myweights) == typeof(1.0)

    myweights = float(Matrix(dtm(crps)));
    bm_25!(myweights)
    @test approx_eq(T.(myweights), correctweights, tol=max_tol)
    @test typeof(myweights) <: Matrix
    @test eltype(myweights) == typeof(1.0)

    @test_throws MethodError bm_25!(DocumentTermMatrix(crps))
    @test bm_25!(DocumentTermMatrix{T}(crps)) isa SparseMatrixCSC{T, Int}
end


text:
## Extend convert to handle JuMP containers
# Array -> SparseAxisArray
function Base.convert(::Type{JuMPC.SparseAxisArray}, arr::Array)
    data = Dict(Tuple(k) => arr[k] for k in CartesianIndices(arr))
    return JuMPC.SparseAxisArray(data)
end

# DenseAxisArray -> SparseAxisArray
function Base.convert(::Type{JuMPC.SparseAxisArray},
                      arr::JuMPC.DenseAxisArray)
    data = Dict(k.I => arr[k] for k in keys(arr))
    return JuMPC.SparseAxisArray(data)
end

# Convert numbers to jump objects
function Base.convert(::Type{JuMP.AbstractJuMPScalar}, c::Number)
    return zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}) + c
end

# Make workaround for keys of containers
_keys(a::JuMPC.SparseAxisArray) = keys(a.data)
_keys(a::AbstractArray) = keys(a)

# Tempororay hack for map of DenseAxisArray
function Base.map(f, a::JuMPC.DenseAxisArray)
    data = map(f, a.data)
    return JuMPC.DenseAxisArray(data, axes(a)...)
end

## Define functions to convert a JuMP array into a vector (need for @BDconstraint)
# AbstractArray
function _make_vector(arr::AbstractArray{T})::Vector{T} where {T}
    if isempty(arr)
        return T[]
    else
        return reduce(vcat, arr)
    end
end

# Array (do nothing)
function _make_vector(arr::Vector{T})::Vector{T} where {T}
    return arr
end

## Define function to vectorize abstractarrays in a consistent way
# Vector
function _make_ordered_vector(
    arr::Vector{T}
    )::Vector{T} where {T}
    return arr
end

# Array
function _make_ordered_vector(
    arr::Union{Array{T}, JuMPC.DenseAxisArray{T}}
    )::Vector{T} where{T}
    return reduce(vcat, arr)
end

# SparseAxisArray
function _make_ordered_vector(
    arr::JuMPC.SparseAxisArray{T}
    )::Vector{T} where {T}
    indices = Collections._get_indices(arr)
    return Collections._make_ordered(arr, indices)
end

# Fallback
function _make_ordered_vector(x::T)::T where {T}
    return x
end

## Define efficient function to check if all elements in array are equal
# method found at https://stackoverflow.com/questions/47564825/check-if-all-the-elements-of-a-julia-array-are-equal/47578613
@inline function _allequal(x::AbstractArray)::Bool
    length(x) < 2 && return true
    e1 = first(x)
    for value in x
        value == e1 || return false
    end
    return true
end

## Convert JuMP variable info to only use Float64
# Change needed
function _make_float_info(info::JuMP.VariableInfo
    )::JuMP.VariableInfo{Float64, Float64, Float64, Float64}
    return JuMP.VariableInfo{Float64, Float64, Float64, Float64}(
        info.has_lb, info.lower_bound, info.has_ub, info.upper_bound,
        info.has_fix, info.fixed_value, info.has_start, info.start,
        info.binary, info.integer)
end

# No change needed
function _make_float_info(
    info::JuMP.VariableInfo{Float64, Float64, Float64, Float64}
    )::JuMP.VariableInfo{Float64, Float64, Float64, Float64}
    return info
end


text:
#=  Defines the actuatorSystem Type that converts control signals from the
    ControlSystem type to commands for the individual actuators
=#
abstract type AbstractActuatorSystem <: DiscreteTimeSystem

########### Template function ###########################
#mutable struct ActuatorSystem{T} <: DiscreteTimeSystem
#    Ts::T #Sample Time
#    Tprev::T #previous evaluation time
#end


text:
export clip

function clip(x, clipval)
    x = (abs(x) > clipval) ? clipval * sign(x) : x
end

function update!(var::Variable, lr)
    # update single Variable
    @. var.value -= lr * setNanInfZero(var.delta)
end

function update!(vars::Vector{Variable}, lr)
    # update multi Variables
    for var in vars
        update!(var, lr)
    end
end

function zerograds!(parameters)
    for var in parameters
        if var.delta !== nothing
            var.delta .= 0.0
        end
    end
end

function zerograds!(O::Optimizer)
    for xparam in O.xparams
        c , θ = xparam
        if θ.delta !== nothing
            θ.delta .= 0.0
        end
    end
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
mutable struct data_type
    opt_val
    have_sol
end

function block_compact_POP(x,f,g,h,k,r)
    n=length(x)
    m=length(g)
    l=length(h)

    dg = Array{UInt16}(undef, m)
    dh = Array{UInt16}(undef,l)

    rf=binomial(k+n,n)
    rg=Array{UInt16}(undef, m)
    rh=Array{UInt16}(undef,l)

    v=monomials(x, 0:k)
    v=reverse(v)

    Uf=union(2*v,monomials(f))
    Ug=Array{Any}(undef, m)
    Uh=Array{Any}(undef, l)

    Tf=v*v'
    Tg=Array{Any}(undef, m)
    Th=Array{Any}(undef, l)

    length_ind_block_f=Array{UInt16}(undef, 1)
    length_ind_block_g=Array{UInt16}(undef, m)
    length_ind_block_h=Array{UInt16}(undef, l)

    ind_block_f=Array{Any}(undef, 1)
    ind_block_g=Array{Any}(undef, m)
    ind_block_h=Array{Any}(undef, l)

    for i = 1:m
        dg[i] = ceil(UInt16,maxdegree(g[i])/2)
        rg[i]=binomial(k-dg[i]+n,n)
        Uf=union(Uf,monomials(g[i]))
        Tg[i]=g[i].*v[1:rg[i]]*v[1:rg[i]]'
    end

    for i = 1:l
        dh[i] = ceil(UInt16,maxdegree(h[i])/2)
        rh[i] = binomial(k-dh[i]+n,n)
        Uf=union(Uf,monomials(h[i]))
        Th[i]=h[i].*v[1:rh[i]]*v[1:rh[i]]'
    end

    for t=1:r
        U=Uf
        for i=1:m
            if t>1 && isempty(Ug[i])==0
                U=union(U,Ug[i])
            end
        end
        for i=1:l
            if t>1 && isempty(Uh[i])==0
                U=union(U,Uh[i])
            end
        end

        Cf=zeros(rf,rf)
        Cf[findall(p -> p in U, Tf)].=1
        Cf[diagind(Cf)] .= 0

        Cg=Array{Any}(undef, m)
        Ch=Array{Any}(undef, l)

        for i=1:m
            Cg[i]=zeros(rg[i],rg[i])

            for p=1:rg[i]
                for q=1:rg[i]
                    if isintersecting(monomials(Tg[i][p,q]),U)==1
                        Cg[i][p,q]=1
                    end
                end
            end
            Cg[i][diagind(Cg[i])] .= 0
        end

        for i=1:l
            Ch[i]=zeros(rh[i],rh[i])

            for p=1:rh[i]
                for q=1:rh[i]
                    if isintersecting(monomials(Th[i][p,q]),U)==1
                        Ch[i][p,q]=1
                    end
                end
            end
            Ch[i][diagind(Ch[i])].= 0
        end

        ind_block_f=connected_components(SimpleGraph(Cf))
        length_ind_block_f=length(ind_block_f)

        for i=1:m
            ind_block_g[i]=connected_components(SimpleGraph(Cg[i]))
            length_ind_block_g[i]=length(ind_block_g[i])
        end

        for i=1:l
            ind_block_h[i]=connected_components(SimpleGraph(Ch[i]))
            length_ind_block_h[i]=length(ind_block_h[i])
        end

        if r >=2
            for j=1:length_ind_block_f
                Uf=union(Uf,vec(v[ind_block_f[j]]*v[ind_block_f[j]]'))
            end

            for i=1:m
                for j=1:length_ind_block_g[i]
                    mono=monomials(g[i])
                    Ug[i]=vec(mono[1]*v[ind_block_g[i][j]]*v[ind_block_g[i][j]]')
                    for p=2:length(mono)
                        Ug[i]=union(Ug[i],vec(mono[p]*v[ind_block_g[i][j]]*v[ind_block_g[i][j]]'))
                    end
                end
            end

            for i=1:l
                for j=1:length_ind_block_h[i]
                    mono=monomials(h[i])
                    Uh[i]=vec(mono[1]*v[ind_block_h[i][j]]*v[ind_block_h[i][j]]')
                    for p=2:length(mono)
                        Uh[i]=union(Uh[i],vec(mono[p]*v[ind_block_h[i][j]]*v[ind_block_h[i][j]]'))
                    end
                end
            end
        end

    end

    model = Model(with_optimizer(Mosek.Optimizer, QUIET=true))

    wSOS=0

    G=Array{Any}(undef, length_ind_block_f)
    for i=1:length_ind_block_f
        G[i]= @variable(model, [1:length(ind_block_f[i]), 1:length(ind_block_f[i])],PSD)
        wSOS +=v[ind_block_f[i]]'*G[i]*v[ind_block_f[i]]
    end

    Gg=Array{Any}(undef, m)
    for j=1:m
        Gg[j]=Array{Any}(undef, length_ind_block_g[j])
        for i=1:length_ind_block_g[j]
            Gg[j][i]= @variable(model, [1:length(ind_block_g[j][i]), 1:length(ind_block_g[j][i])],PSD)
            wSOS += g[j]*v[ind_block_g[j][i]]'*Gg[j][i]*v[ind_block_g[j][i]]
        end
    end

    Gh=Array{Any}(undef, l)
    length_Gh=Array{Any}(undef, l)
    vh=Array{Any}(undef, l)
    for j=1:l
        Gh[j]=Array{Any}(undef, length_ind_block_h[j])
        length_Gh[j]=Array{Any}(undef, length_ind_block_h[j])
        vh[j]=Array{Any}(undef, length_ind_block_h[j])
        for i=1:length_ind_block_h[j]
            vh[j][i]= unique(vec(v[ind_block_h[j][i]]*v[ind_block_h[j][i]]'))
            length_Gh[j][i]= length(vh[j][i])
            Gh[j][i]=  @variable(model, [1:length_Gh[j][i]])
            wSOS += h[j]*vh[j][i]'*Gh[j][i]
        end
    end

    @variable(model, lambda)

    @constraint(model, coefficients(f-lambda - wSOS).== 0)

    @objective(model, Max, lambda)

    optimize!(model)

    opt_val = value(lambda)

    println("Index of block for f : ", ind_block_f)
    println("------------------------------------")
    println("Indices of block for g : ", ind_block_g)
    println("------------------------------------")
    println("Indices of block for h ; ", ind_block_h)
    println("------------------------------------")
    println("Termination status = ", termination_status(model))

    println("Optimal value = ",opt_val)

    have_sol=0

    println("------------------------------------")
    G0=zeros(rf,rf)
    for i=1:length_ind_block_f
        G0[ind_block_f[i],ind_block_f[i]]=value.(G[i])
    end

    rk=rf-rank(G0, 1e-4)

    println("Rank of moment matrix = ", rk)

    if rk >0
        # extraction of Henrion and Lasserre
        F = eigen(G0)
        V = F.vectors

        Ix=sortperm(real.(F.values))

        V=V[:,Ix[1:rk]]
        V=Matrix(V')
        V= rref_with_pivots!(V,1e-3);
        U=V[1]

        U=Matrix(U')
        # Figure out multiplying matrices using YALMIP code
        w=v[V[2]];
        N=zeros(length(V[2]),rk,n)
        flag=1
        for i in 1:n
            xw=x[i]*w
            kk=indexin(xw,v)
            if nothing in kk
                flag=0
                break
            else
                N[:,:,i]=U[kk,:]
            end
        end

        if flag==1

            # Create random convex combination
            rands = rand(n,1);rands = rands/sum(rands);
            M = zeros(length(V[2]),rk);
            for i in 1:n
                M+=rands[i]*N[:,:,i];
            end

            F= schur(M);
            L=F.Z
            # Extract solution
            for i in 1:rk
                atom=[]
                for j = 1:n
                    atom=[atom;L[:,i]'*N[:,:,j]*L[:,i]];
                end
                println("------------------------------------")
                println("atom ",i," = ",atom)
                flag=1
                check=opt_val-polynomial(f)(x => atom)
                println("check lower bound  = ",check)
                if abs(check)>1e-2
                    flag=0
                end
                for j in 1:m
                    check=polynomial(g[j])(x => atom)
                    println("check inequality ",j," = ",check)
                    if check<-1e-3
                        flag=0
                    end
                end

                for j in 1:l
                    check=polynomial(h[j])(x => atom)
                    println("check equality ",j," = ",check)
                    if abs(check)>1e-3
                        flag=0
                    end
                end
                if flag ==1
                    have_sol=1
                    sol=atom
                    println("####################################")
                    println("Solution = ",sol)
                    println("####################################")
                end

            end
        end
    end
    data=data_type(opt_val,have_sol)

return  data
end

function block_noncompact_POP(x,f,g,h,eps,k,r)
    n=length(x)
    m=length(g)
    l=length(h)

    if m==0 && l==0
        d = ceil(UInt16,maxdegree(f)/2)
    else
        d = floor(UInt16,maxdegree(f)/2) + 1
    end
    dg = Array{UInt16}(undef, m)
    dh = Array{UInt16}(undef,l)

    rf=binomial(k+d+n,n)
    rg=Array{UInt16}(undef, m)
    rh=Array{UInt16}(undef,l)

    v=monomials(x, 0:d+k)
    v=reverse(v)

    theta=1+x'*x

    Uf=union(2*v,monomials(theta^k*f))
    Ug=Array{Any}(undef, m)
    Uh=Array{Any}(undef, l)

    Tf=v*v'
    Tg=Array{Any}(undef, m)
    Th=Array{Any}(undef, l)

    length_ind_block_f=Array{UInt16}(undef, 1)
    length_ind_block_g=Array{UInt16}(undef, m)
    length_ind_block_h=Array{UInt16}(undef, l)

    ind_block_f=Array{Any}(undef, 1)
    ind_block_g=Array{Any}(undef, m)
    ind_block_h=Array{Any}(undef, l)

    for i = 1:m
        dg[i] = ceil(UInt16,maxdegree(g[i])/2)
        rg[i]=binomial(k+d-dg[i]+n,n)
        Uf=union(Uf,monomials(g[i]))
        Tg[i]=g[i].*v[1:rg[i]]*v[1:rg[i]]'
    end

    for i = 1:l
        dh[i] = ceil(UInt16,maxdegree(h[i])/2)
        rh[i] = binomial(k+d-dh[i]+n,n)
        Uf=union(Uf,monomials(h[i]))
        Th[i]=h[i].*v[1:rh[i]]*v[1:rh[i]]'
    end

    for t=1:r
        U=Uf
        for i=1:m
            if t>1 && isempty(Ug[i])==0
                U=union(U,Ug[i])
            end
        end
        for i=1:l
            if t>1 && isempty(Uh[i])==0
                U=union(U,Uh[i])
            end
        end

        Cf=zeros(rf,rf)
        Cf[findall(p -> p in U, Tf)].=1
        Cf[diagind(Cf)] .= 0

        Cg=Array{Any}(undef, m)
        Ch=Array{Any}(undef, l)

        for i=1:m
            Cg[i]=zeros(rg[i],rg[i])

            for p=1:rg[i]
                for q=1:rg[i]
                    if isintersecting(monomials(Tg[i][p,q]),U)==1
                        Cg[i][p,q]=1
                    end
                end
            end
            Cg[i][diagind(Cg[i])] .= 0
        end

        for i=1:l
            Ch[i]=zeros(rh[i],rh[i])

            for p=1:rh[i]
                for q=1:rh[i]
                    if isintersecting(monomials(Th[i][p,q]),U)==1
                        Ch[i][p,q]=1
                    end
                end
            end
            Ch[i][diagind(Ch[i])].= 0
        end

        ind_block_f=connected_components(SimpleGraph(Cf))
        length_ind_block_f=length(ind_block_f)

        for i=1:m
            ind_block_g[i]=connected_components(SimpleGraph(Cg[i]))
            length_ind_block_g[i]=length(ind_block_g[i])
        end

        for i=1:l
            ind_block_h[i]=connected_components(SimpleGraph(Ch[i]))
            length_ind_block_h[i]=length(ind_block_h[i])
        end

        if r >=2
            for j=1:length_ind_block_f
                Uf=union(Uf,vec(v[ind_block_f[j]]*v[ind_block_f[j]]'))
            end

            for i=1:m
                for j=1:length_ind_block_g[i]
                    mono=monomials(g[i])
                    Ug[i]=vec(mono[1]*v[ind_block_g[i][j]]*v[ind_block_g[i][j]]')
                    for p=2:length(mono)
                        Ug[i]=union(Ug[i],vec(mono[p]*v[ind_block_g[i][j]]*v[ind_block_g[i][j]]'))
                    end
                end
            end

            for i=1:l
                for j=1:length_ind_block_h[i]
                    mono=monomials(h[i])
                    Uh[i]=vec(mono[1]*v[ind_block_h[i][j]]*v[ind_block_h[i][j]]')
                    for p=2:length(mono)
                        Uh[i]=union(Uh[i],vec(mono[p]*v[ind_block_h[i][j]]*v[ind_block_h[i][j]]'))
                    end
                end
            end
        end

    end

    model = Model(with_optimizer(Mosek.Optimizer, QUIET=true))

    wSOS=0

    G=Array{Any}(undef, length_ind_block_f)
    for i=1:length_ind_block_f
        G[i]= @variable(model, [1:length(ind_block_f[i]), 1:length(ind_block_f[i])],PSD)
        wSOS += v[ind_block_f[i]]'*G[i]*v[ind_block_f[i]]
    end

    Gg=Array{Any}(undef, m)
    for j=1:m
        Gg[j]=Array{Any}(undef, length_ind_block_g[j])
        for i=1:length_ind_block_g[j]
            Gg[j][i]= @variable(model, [1:length(ind_block_g[j][i]), 1:length(ind_block_g[j][i])],PSD)
            wSOS += g[j]*v[ind_block_g[j][i]]'*Gg[j][i]*v[ind_block_g[j][i]]
        end
    end

    Gh=Array{Any}(undef, l)
    length_Gh=Array{Any}(undef, l)
    vh=Array{Any}(undef, l)
    for j=1:l
        Gh[j]=Array{Any}(undef, length_ind_block_h[j])
        length_Gh[j]=Array{Any}(undef, length_ind_block_h[j])
        vh[j]=Array{Any}(undef, length_ind_block_h[j])
        for i=1:length_ind_block_h[j]
            vh[j][i]= unique(vec(v[ind_block_h[j][i]]*v[ind_block_h[j][i]]'))
            length_Gh[j][i]= length(vh[j][i])
            Gh[j][i]=  @variable(model, [1:length_Gh[j][i]])
            wSOS +=h[j]*vh[j][i]'*Gh[j][i]
        end
    end

    @variable(model, lambda)

    @constraint(model, coefficients(theta^k*(f-lambda+eps*theta^d) - wSOS).== 0)

    @objective(model, Max, lambda)

    optimize!(model)

    opt_val = value(lambda)

    println("Index of block for f : ", ind_block_f)
    println("------------------------------------")
    println("Indices of block for g : ", ind_block_g)
    println("------------------------------------")
    println("Indices of block for h ; ", ind_block_h)
    println("------------------------------------")
    println("Termination status = ", termination_status(model))

    println("Optimal value = ",opt_val)

return opt_val
end

function block_noncompact_POP_with_lower_bound(x,f,g,h,eps,k,r)

    n=length(x)
    m=length(g)
    l=length(h)

    d = ceil(UInt16,maxdegree(f)/2)

    dg = Array{UInt16}(undef, m)
    dh = Array{UInt16}(undef,l)

    rf=binomial(k+d+n,n)
    rg=Array{UInt16}(undef, m)
    rh=Array{UInt16}(undef,l)

    v=monomials(x, 0:d+k)
    v=reverse(v)

    theta=1+x'*x

    Uf=union(2*v,monomials(theta^k*f))
    Ug=Array{Any}(undef, m)
    Uh=Array{Any}(undef, l)

    Tf=v*v'
    Tg=Array{Any}(undef, m)
    Th=Array{Any}(undef, l)

    length_ind_block_f=Array{UInt16}(undef, 1)
    length_ind_block_g=Array{UInt16}(undef, m)
    length_ind_block_h=Array{UInt16}(undef, l)

    ind_block_f=Array{Any}(undef, 1)
    ind_block_g=Array{Any}(undef, m)
    ind_block_h=Array{Any}(undef, l)

    for i = 1:m
        dg[i] = ceil(UInt16,maxdegree(g[i])/2)
        rg[i]=binomial(k+d-dg[i]+n,n)
        Uf=union(Uf,monomials(g[i]))
        Tg[i]=g[i].*v[1:rg[i]]*v[1:rg[i]]'
    end

    for i = 1:l
        dh[i] = ceil(UInt16,maxdegree(h[i])/2)
        rh[i] = binomial(k+d-dh[i]+n,n)
        Uf=union(Uf,monomials(h[i]))
        Th[i]=h[i].*v[1:rh[i]]*v[1:rh[i]]'
    end

    for t=1:r
        U=Uf
        for i=1:m
            if t>1 && isempty(Ug[i])==0
                U=union(U,Ug[i])
            end
        end
        for i=1:l
            if t>1 && isempty(Uh[i])==0
                U=union(U,Uh[i])
            end
        end

        Cf=zeros(rf,rf)
        Cf[findall(p -> p in U, Tf)].=1
        Cf[diagind(Cf)] .= 0

        Cg=Array{Any}(undef, m)
        Ch=Array{Any}(undef, l)

        for i=1:m
            Cg[i]=zeros(rg[i],rg[i])

            for p=1:rg[i]
                for q=1:rg[i]
                    if isintersecting(monomials(Tg[i][p,q]),U)==1
                        Cg[i][p,q]=1
                    end
                end
            end
            Cg[i][diagind(Cg[i])] .= 0
        end

        for i=1:l
            Ch[i]=zeros(rh[i],rh[i])

            for p=1:rh[i]
                for q=1:rh[i]
                    if isintersecting(monomials(Th[i][p,q]),U)==1
                        Ch[i][p,q]=1
                    end
                end
            end
            Ch[i][diagind(Ch[i])].= 0
        end

        ind_block_f=connected_components(SimpleGraph(Cf))
        length_ind_block_f=length(ind_block_f)

        for i=1:m
            ind_block_g[i]=connected_components(SimpleGraph(Cg[i]))
            length_ind_block_g[i]=length(ind_block_g[i])
        end

        for i=1:l
            ind_block_h[i]=connected_components(SimpleGraph(Ch[i]))
            length_ind_block_h[i]=length(ind_block_h[i])
        end

        if r >=2
            for j=1:length_ind_block_f
                Uf=union(Uf,vec(v[ind_block_f[j]]*v[ind_block_f[j]]'))
            end

            for i=1:m
                for j=1:length_ind_block_g[i]
                    mono=monomials(g[i])
                    Ug[i]=vec(mono[1]*v[ind_block_g[i][j]]*v[ind_block_g[i][j]]')
                    for p=2:length(mono)
                        Ug[i]=union(Ug[i],vec(mono[p]*v[ind_block_g[i][j]]*v[ind_block_g[i][j]]'))
                    end
                end
            end

            for i=1:l
                for j=1:length_ind_block_h[i]
                    mono=monomials(h[i])
                    Uh[i]=vec(mono[1]*v[ind_block_h[i][j]]*v[ind_block_h[i][j]]')
                    for p=2:length(mono)
                        Uh[i]=union(Uh[i],vec(mono[p]*v[ind_block_h[i][j]]*v[ind_block_h[i][j]]'))
                    end
                end
            end
        end

    end

    model = Model(with_optimizer(Mosek.Optimizer, QUIET=true))

    wSOS=0

    G=Array{Any}(undef, length_ind_block_f)
    for i=1:length_ind_block_f
        G[i]= @variable(model, [1:length(ind_block_f[i]), 1:length(ind_block_f[i])],PSD)
        wSOS += v[ind_block_f[i]]'*G[i]*v[ind_block_f[i]]
    end

    Gg=Array{Any}(undef, m)
    for j=1:m
        Gg[j]=Array{Any}(undef, length_ind_block_g[j])
        for i=1:length_ind_block_g[j]
            Gg[j][i]= @variable(model, [1:length(ind_block_g[j][i]), 1:length(ind_block_g[j][i])],PSD)
            wSOS += g[j]*v[ind_block_g[j][i]]'*Gg[j][i]*v[ind_block_g[j][i]]
        end
    end

    Gh=Array{Any}(undef, l)
    length_Gh=Array{Any}(undef, l)
    vh=Array{Any}(undef, l)
    for j=1:l
        Gh[j]=Array{Any}(undef, length_ind_block_h[j])
        length_Gh[j]=Array{Any}(undef, length_ind_block_h[j])
        vh[j]=Array{Any}(undef, length_ind_block_h[j])
        for i=1:length_ind_block_h[j]
            vh[j][i]= unique(vec(v[ind_block_h[j][i]]*v[ind_block_h[j][i]]'))
            length_Gh[j][i]= length(vh[j][i])
            Gh[j][i]=  @variable(model, [1:length_Gh[j][i]])
            wSOS +=h[j]*vh[j][i]'*Gh[j][i]
        end
    end

    @variable(model, lambda)

    @constraint(model, coefficients(theta^k*(f-lambda+eps*theta^d) - wSOS).== 0)

    @objective(model, Max, lambda)

    optimize!(model)

    opt_val = value(lambda)

    println("Index of block for f : ", ind_block_f)
    println("------------------------------------")
    println("Indices of block for g : ", ind_block_g)
    println("------------------------------------")
    println("Indices of block for h ; ", ind_block_h)
    println("------------------------------------")
    println("Termination status = ", termination_status(model))

    println("Optimal value = ",opt_val)

return opt_val
end

function adding_spherical_constraints(x,g,h,k,r)

    # small parameter
    eps = 1e-1

    n=length(x)
    m=length(g)
    l=length(h)

    dg = Array{UInt16}(undef, m)
    dh = Array{UInt16}(undef,l)

    for i = 1:m
        dg[i] = ceil(UInt16,maxdegree(g[i])/2)
    end

    for i = 1:l
        dh[i] = ceil(UInt16,maxdegree(h[i])/2)
    end

    d_max=maximum(union(union(dg,dh),1))

    # Define centers and square of radius
    a0=zeros(Float64,(n, 1)); a = Matrix{Float64}(I, n, n)

    println("Determine L0:")

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

    # Polynomial to optimize
    f = [x-a0]'*[x-a0]

    if m==0
        g=[f]
    else
        g=[g;f]
    end

    L0=block_noncompact_POP_with_lower_bound(x,f,g,h,eps,k,r)

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

    println("L0 = ", L0)

    println("====================================")

    # Define omegat, t=0,...,n
    omega0 = 0; omega = zeros(n)

    #inequalities polynomial
    g[end]=L0-f

    println("Determine omega",0,":")

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

    data=block_compact_POP(x,f,g,h,k+d_max,r)

    omega0=data.opt_val

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

    println("omega",0," = ", omega0)

    println("====================================")

    #equalities polynomial
    if l==0
        h=[omega0-f]
    else
        h=[h;omega0-f]
    end

    #inequalities polynomial
    g=g[1:end-1]

    if data.have_sol==0
        for t=1:n

            println("Determine omega",t,":")

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

            if t>1
                #equalities polynomial
                h = [h;omega[t-1]-f] ; l=length(h)
            end

            f=[x-a[:,t]]'*[x-a[:,t]]

            data=block_compact_POP(x,f,g,h,k+d_max,r)
            omega[t]=data.opt_val
            println("------------------------------------")
            println("omega",t," = ", omega[t])

            println("====================================")
            if data.have_sol==1
                break
            end
        end
    end
end

function isintersecting(v1, v2)
    flag=0
    for v in v1
        if v in v2
            flag=1
            break
        end
    end
return flag
end


text:
function get_ap(obj::NRLMSISE00InputParameters, arg0::AbsoluteDate)
    return jcall(obj, "getAp", Vector{jdouble}, (AbsoluteDate,), arg0)
end

function get_average_flux(obj::NRLMSISE00InputParameters, arg0::AbsoluteDate)
    return jcall(obj, "getAverageFlux", jdouble, (AbsoluteDate,), arg0)
end

function get_daily_flux(obj::NRLMSISE00InputParameters, arg0::AbsoluteDate)
    return jcall(obj, "getDailyFlux", jdouble, (AbsoluteDate,), arg0)
end

function get_max_date(obj::NRLMSISE00InputParameters)
    return jcall(obj, "getMaxDate", AbsoluteDate, ())
end

function get_min_date(obj::NRLMSISE00InputParameters)
    return jcall(obj, "getMinDate", AbsoluteDate, ())
end


text:
using Compat, ClassicalLaminateTheory

forcesAndMoments =
 @compat Dict(
    :Nx => 2000,
    :Ny => 0,
    :Nxy => 0,
    :Mx => 0,
    :My => 0,
    :Mxy => 0)

materialProperties =@compat Dict(
  :E1 => 20.0*10^6,
  :E2 => 1.5*10^6,
  :G12 => 1.0*10^6,
  :ν12 => .34,
  :alpha1 => 0.2e-6,
  :alpha2 => 0.2e-4,
  :beta1 => 0.2e-6,
  :beta2 => 0.2e-4,
  :thickness => 0.0075)

laminateProperties =@compat Dict(
  :nplies => 2,
  :repeats => 1,
  :symmetric => false,
  :orientation => [0, 90],
  :deltaTemp => 0.0,
  :deltaMoisture => 0
)

println()
createLaminate!(laminateProperties, materialProperties, forcesAndMoments)
show(laminateProperties)

println("ABD matrix:")
round(laminateProperties[:ABD], 3) |> display

println()
println("Total forces and moments:")
round(laminateProperties[:fTotal], 7) |> display

[laminateProperties[:f][1], laminateProperties[:f][2]]


text:
using Mmap
using ImageBinarization
using Images
using Colors
using Dates

log(msg::String) = println("$(Dates.now()) $msg")

dirs = ["/run/media/jakub/3315-3BCD/trainset3/circles", "/run/media/jakub/3315-3BCD/trainset3/triangles/", "/run/media/jakub/3315-3BCD/trainset3/rectangles/"] # input dirs

function toBitMatrix(img::Matrix)
    noalpha = RGB{N0f8}.(img)
    binarized = binarize(noalpha, Otsu())

    Bool.(binarized)
end

toVector(m::BitMatrix) = reshape(m, length(m), 1)[1:end]

function numberOfFiles(dir::String)
    output = readchomp(`src/numberOfFilesInDir.sh $dir`) # faster than readdir: does not collect file names

    return parse(Int64, output)
end

vectorLen = 30*30 # the images have these dimensions
totalNumberOfFiles = sum(numberOfFiles.(dirs))
println(totalNumberOfFiles)

dataFile = open("/home/jakub/Dokumenty/skola/siete/neural/trainset/trainset4", "w+") # output data file
dataMap = mmap(dataFile, BitMatrix, (vectorLen, totalNumberOfFiles))

targetsFile = open("/home/jakub/Dokumenty/skola/siete/neural/trainset/targets4", "w+") # output targets file
targetsMap = mmap(targetsFile, BitMatrix, (3, totalNumberOfFiles))

columnIndex = 1

for (i, dir) in enumerate(dirs)
    target = Bool.(==(i, j) for j in 1:3)

    log("Entering directory $dir")
    for file in readdir(dir, join = true)
        global columnIndex

        log("Loading file $file")
        img = load(file)
        log("Converting")
        img = toBitMatrix(img)
        img = toVector(img)

        log("Adding to matrix")
        dataMap[:, columnIndex] = img
        targetsMap[:, columnIndex] = target

        columnIndex += 1
    end
end

max_stars_repo_name: achjaj/shape-recognition

max_issues_repo_name: achjaj/shape-recognition

max_forks_repo_name: achjaj/shape-recognition


text:
############related to gamma function
#see https://math.stackexchange.com/questions/1441753/approximating-the-digamma-function
#and https://en.wikipedia.org/wiki/Digamma_function
function cudaapproxdigamma(x::Real)
  adj = 0.0f0
  ψ = x
  #the polynomial is only accurate for large x
  for i ∈ 1:3
    adj -= 1f0/ψ
    ψ += 1f0
  end

  ψ = (CUDA.log(ψ) - 1f0 / 2f0 / ψ - 1f0 / 12f0 / (ψ * ψ) +
    1f0/120f0 * CUDA.pow(ψ, -4) - 1f0/252f0 * CUDA.pow(ψ, -6) +
    1f0/240f0 * CUDA.pow(ψ, -8) - 5f0/660f0 * CUDA.pow(ψ, -10) +
    691f0/32760f0 * CUDA.pow(ψ, -12) - 1f0/12f0 * CUDA.pow(ψ, -14)) + adj

  return ψ
end

#runs a customized gamma log pdf that uses a gpu array when available
#adapted from StatsFuns
#https://github.com/JuliaStats/StatsFuns.jl/blob/master/src/distrs/gamma.jl
#and CUDA code
#GPU vectorized log gamma pdf
function vecgammalogpdf(k::TV, θ::TV, x::TV)::TV where TV<:CUDA.CuVector
  Γs_part = cudagammalogpdf_part.(k,θ,x)

  return Γs_part .- cudaveclgamma(k)
end

#does everything except compute the log gamma function
#adapted from StatsFuns
#https://github.com/JuliaStats/StatsFuns.jl/blob/master/src/distrs/gamma.jl
cudagammalogpdf_part(k,θ,x) = - k * CUDA.log(θ) + (k - 1f0) * CUDA.log(x) - x / θ

#vectorized log gamma function
cudaveclgamma(k) = (CUDA.lgamma).(k)

#define the adjoint
Zygote.@adjoint cudaveclgamma(k) = cudaveclgamma(k), y->(y .* cudaapproxdigamma.(k), )
Zygote.refresh()

#lightly optimized cpu version for testing
function vecgammalogpdf(k::Any, θ::Any, x::TV)::TV where TV<:AbstractVector

  Γs::TV = -(SpecialFunctions.loggamma).(k) .- k .* log.(θ) .+ (k .- 1f0) .* log.(x) .- x ./ θ

  return Γs
end

slowvecgammalogpdf(k, θ, x) = ((kᵢ, θᵢ, xᵢ)->logpdf(Gamma(kᵢ,θᵢ), xᵢ)).(k, θ, x)

#test function for vecgammalogpdf
#test it
function testgammagrad(N::Int)
  τ = abs.(randn(Float32, N)) .+ .001f0
  absμ = abs.(randn(Float32, N)) .+ .001f0

  θ = 1.0f0 ./ τ
  k = absμ .* τ
  v = rand(Float32, N)

  cuθ = θ |> gpu
  cuk = k |> gpu
  cuv = v |> gpu

  @info "testing function values (N=$N) cpu: $(sum(vecgammalogpdf(k, θ, v))) " *
    "gpu: $(sum(vecgammalogpdf(cuk, cuθ, cuv)))"
  @info "Additional check of value: $(sum(slowvecgammalogpdf(k, θ, v)))"

  @info "cpu version gradient N=$N"
  if N ≤ 10
    gs = gradient((k,θ)->sum(vecgammalogpdf(k, θ, v)),k,θ)
    @info "gs: $gs"
  else
    gs = @btime gradient((k,θ)->sum(vecgammalogpdf(k, θ, $v)),$k,$θ)
  end

  @info "gpu version gradient N=$N"
  if N ≤ 10
    cugs = gradient((cuk,cuθ)->sum(vecgammalogpdf(cuk, cuθ, cuv)),cuk,cuθ)
    @info "cugs: $cugs"
  else
    cugs = @btime gradient((cuk,cuθ)->sum(vecgammalogpdf(cuk, cuθ, $cuv)),$cuk,$cuθ)
  end
end

###log normal dpf
#sortof a cuda kernal
using Revise, BenchmarkTools, Distributions, Zygote, Flux, CUDA, CUDA
function cudanormlogprop(μᵢ::Real, logτᵢ::Real, xᵢ::Real)
  μMxᵢ = μᵢ-xᵢ
  return -1f0 * μMxᵢ*μMxᵢ*CUDA.exp(logτᵢ) + logτᵢ
end

function vecnormlogprop(μ::TV, logτ::TV, x::TV)::TV where TV<:AbstractVector
  uMx = μ .- x
  return -1f0 .* uMx .* uMx .* exp.(logτ) .+ logτ
end

#Zygote.@adjoint cudaveclgamma(k) = cudaveclgamma(k), y->(y .* cudaapproxdigamma.(k), )
#Zygote.refresh()
function vecnormlogprop(μ::TV, logτ::TV, x::TV)::TV where TV<:CUDA.CuVector
  return cudanormlogprop.(μ,logτ,x)
end
cudanormlogpdffromprop(x::Real) = 0.5f0*x - 0.5f0*log(2f0*3.1415927f0)
slowvecnormlogpdf(μ, σ, x) = ((μᵢ, σᵢ, xᵢ)->logpdf(Normal(μᵢ,σᵢ), xᵢ)).(μ, σ, x)

function vecnormlogpdf(μ::TV, logτ::TV, x::TV)::TV where TV<:CUDA.CuVector
  return cudanormlogpdffromprop.(vecnormlogprop(μ, logτ, x))
end

function vecnormlogpdf(μ::TV, logτ::TV, x::TV)::TV where TV<:AbstractVector
  return 0.5f0 .* vecnormlogprop(μ, logτ, x) .- 0.5f0 .* log.(2f0 .* 3.1415927f0)
end

function testnormgrad(N::Int)
  τ = log.(abs.(randn(Float32, N)) .+ .001f0)
  μ = abs.(randn(Float32, N)) .+ .001f0

  σ = (1.0f0 ./ exp.(τ)).^0.5f0
  v = rand(Float32, N)

  cuμ = μ |> gpu
  cuτ = τ |> gpu
  cuσ = σ |> gpu
  cuv = v |> gpu

  @info "testing prop values (N=$N) cpu: $(sum(vecnormlogprop(μ, τ, v))) " *
    "gpu: $(sum(vecnormlogprop(cuμ, cuτ, cuv)))"
  @info "corrected cpu values: $(sum(vecnormlogpdf(μ, τ, v)))"
  @info "corrected gpu values: $(sum(vecnormlogpdf(cuμ, cuτ, cuv)))"
  @info "Additional check of value: $(sum(slowvecnormlogpdf(μ, σ, v)))"

  @info "cpu version gradient N=$N"
  if N ≤ 10
    gs = gradient((μ, τ)->sum(vecnormlogprop(μ, τ, v)),μ, τ)
    @info "gs: $gs"
  else
    gs = @btime gradient((μ,τ)->sum(vecnormlogprop(μ, τ, $v)),$μ,$τ)
  end

  @info "gpu version gradient N=$N"
  if N ≤ 10
    cugs = gradient((cuμ,cuτ)->sum(vecnormlogprop(cuμ, cuτ, cuv)),cuμ,cuτ)
    @info "cugs: $cugs"
  else
    cugs = @btime gradient((cuμ,cuτ)->sum(vecnormlogprop(cuμ, cuτ, $cuv)),$cuμ,$cuτ)
  end
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# =============
# weight matrix
# =====================

# TODO would probably be better if the matrix stores the adjacency list instead of the key
# TODO These methods are not typestable

struct SimpleValueMatrix{T, K <: Union{Nothing, Val}, G <: AbstractSimpleValueGraph} <: AbstractMatrix{T}
    g::G
end

function weights(g::AbstractSimpleValueGraph{V, E_VAL}) where {V, E_VAL}
    SimpleValueMatrix{E_VAL, Nothing, typeof(g)}(g)
end

function weights(g::AbstractSimpleValueGraph{V, E_VAL}, key) where {V, E_VAL <: TupleOrNamedTuple}
    T = eltype(eltype(g.edgevals[key]))
    SimpleValueMatrix{T, Val{key}, typeof(g)}(g)
end

function getindex(A::SimpleValueMatrix{T, Nothing}, s::Integer, d::Integer) where {T}
    result = get_edgeval(A.g, s, d)
    return ifelse(result == nothing, default_zero_edgeval(T), result)
end

function getindex(A::SimpleValueMatrix{T, Val{key}}, s::Integer, d::Integer) where {T, key}
    result = get_edgeval(A.g, s, d, key)
    return ifelse(result == nothing, default_zero_edgeval(T), result)
end

function size(A::SimpleValueMatrix)
    n = Int(nv(A.g))
    return (n, n)
end

function replace_in_print_matrix(A::SimpleValueMatrix, s::Integer, d::Integer, str::AbstractString)
    has_edge(A.g, s, d) ? str : Base.replace_with_centered_mark(str)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

typealias MortarElements2D Union{Seg2, Seg3}

function newton(f, df, x; tol=1.0e-6, max_iterations=10)
    for i=1:max_iterations
        dx = -f(x)/df(x)
        x += dx
        if norm(dx) < tol
            return x
        end
    end
    error("Newton iteration did not converge in $max_iterations iterations")
end

function cross2(a, b)
    cross([a; 0], [b; 0])[3]
end

function get_slave_elements(problem::Problem)
    filter(el -> haskey(el, "master elements"), get_elements(problem))
end

function project_from_master_to_slave{E<:MortarElements2D}(slave_element::Element{E}, x2, time)
    x1_ = slave_element["geometry"](time)
    n1_ = slave_element["normal"](time)
    x1(xi1) = vec(get_basis(slave_element, [xi1], time))*x1_
    dx1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*x1_
    n1(xi1) = vec(get_basis(slave_element, [xi1], time))*n1_
    dn1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*n1_
    R(xi1) = cross2(x1(xi1)-x2, n1(xi1))
    dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1))
    xi1 = nothing
    try
        xi1 = newton(R, dR, 0.0)
    catch
        warn("projection from master to slave failed with following arguments:")
        warn("slave element x1: $x1_")
        warn("slave element n1: $n1_")
        warn("master element x2: $x2")
        warn("time: $time")
        len = norm(x1_[2] - x1_[1])
        midpnt = mean(x1_)
        dist = norm(midpnt - x2)
        distval = dist/len
        warn("midpoint of slave element: $midpnt")
        warn("length of slave element: $len")
        warn("distance between midpoint of slave element and x2: $dist")
        warn("charasteristic measure: $distval")
        rethrow()
    end
    return xi1
end

function project_from_slave_to_master{E<:MortarElements2D}(master_element::Element{E}, x1, n1, time)
    x2_ = master_element["geometry"](time)
    x2(xi2) = vec(get_basis(master_element, [xi2], time))*x2_
    dx2(xi2) = vec(get_dbasis(master_element, [xi2], time))*x2_
    cross2(a, b) = cross([a; 0], [b; 0])[3]
    R(xi2) = cross2(x2(xi2)-x1, n1)
    dR(xi2) = cross2(dx2(xi2), n1)
    xi2 = newton(R, dR, 0.0)
    return xi2
end

function calculate_normals(elements, time, ::Type{Val{1}}; rotate_normals=false)
    tangents = Dict{Int64, Vector{Float64}}()
    for element in elements
        conn = get_connectivity(element)
        #X1 = element("geometry", time)
        #dN = get_dbasis(element, [0.0], time)
        #tangent = vec(sum([kron(dN[:,i], X1[i]') for i=1:length(X1)]))
        tangent = vec(element([0.0], time, Val{:Jacobian}))
        for nid in conn
            if haskey(tangents, nid)
                tangents[nid] += tangent
            else
                tangents[nid] = tangent
            end
        end
    end

    Q = [0.0 -1.0; 1.0 0.0]
    normals = Dict{Int64, Vector{Float64}}()
    S = collect(keys(tangents))
    for j in S
        tangents[j] /= norm(tangents[j])
        normals[j] = Q*tangents[j]
    end

    if rotate_normals
        for j in S
            normals[j] = -normals[j]
        end
    end

    return normals, tangents
end

function calculate_normals!(elements, time, ::Type{Val{1}}; rotate_normals=false)
    normals, tangents = calculate_normals(elements, time, Val{1}; rotate_normals=rotate_normals)
    for element in elements
        conn = get_connectivity(element)
        update!(element, "normal", time => [normals[j] for j in conn])
        update!(element, "tangent", time => [tangents[j] for j in conn])
    end
end

function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Type{Val{false}})

    props = problem.properties
    field_dim = get_unknown_field_dimension(problem)
    field_name = get_parent_field_name(problem)
    slave_elements = get_slave_elements(problem)

    # 1. calculate nodal normals and tangents for slave element nodes j ∈ S
    normals, tangents = calculate_normals(slave_elements, time, Val{1};
                                          rotate_normals=props.rotate_normals)
    update!(slave_elements, "normal", time => normals)
    update!(slave_elements, "tangent", time => tangents)

    # 2. loop all slave elements
    for slave_element in slave_elements

        nsl = length(slave_element)
        X1 = slave_element("geometry", time)
        n1 = slave_element("normal", time)

        # 3. loop all master elements
        for master_element in slave_element("master elements", time)

            nm = length(master_element)
            X2 = master_element("geometry", time)

            # 3.1 calculate segmentation
            xi1a = project_from_master_to_slave(slave_element, X2[1], time)
            xi1b = project_from_master_to_slave(slave_element, X2[2], time)
            xi1 = clamp([xi1a; xi1b], -1.0, 1.0)
            l = 1/2*abs(xi1[2]-xi1[1])
            isapprox(l, 0.0) && continue # no contribution in this master element

            # 3.2. bi-orthogonal basis
            De = zeros(nsl, nsl)
            Me = zeros(nsl, nsl)
            Ae = zeros(nsl, nsl)
            if props.dual_basis
                for ip in get_integration_points(slave_element, 3)
                    detJ = slave_element(ip, time, Val{:detJ})
                    w = ip.weight*detJ*l
                    xi = ip.coords[1]
                    xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
                    N1 = vec(get_basis(slave_element, xi_s, time))
                    De += w*diagm(N1)
                    Me += w*N1*N1'
                end
                Ae = De*inv(Me)
            else
                Ae = eye(nsl)
            end

            # 3.3. loop integration points of one integration segment and calculate
            # local mortar matrices
            fill!(De, 0.0)
            fill!(Me, 0.0)
            ge = zeros(field_dim*nsl)
            for ip in get_integration_points(slave_element, 2)
                detJ = slave_element(ip, time, Val{:detJ})
                w = ip.weight*detJ*l
                xi = ip.coords[1]
                xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
                N1 = vec(get_basis(slave_element, xi_s, time))
                Phi = Ae*N1
                # project gauss point from slave element to master element in direction n_s
                X_s = N1*X1 # coordinate in gauss point
                n_s = N1*n1 # normal direction in gauss point
                xi_m = project_from_slave_to_master(master_element, X_s, n_s, time)
                N2 = vec(get_basis(master_element, xi_m, time))
                X_m = N2*X2
                De += w*Phi*N1'
                Me += w*Phi*N2'
                if props.adjust
                    haskey(slave_element, "displacement") || continue
                    haskey(master_element, "displacement") || continue
                    norm(mean(X1) - X2[1]) / norm(X1[2] - X1[1]) < props.distval || continue
                    norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) < props.distval || continue
                    u1 = slave_element("displacement", time)
                    u2 = master_element("displacement", time)
                    x_s = X_s + N1*u1
                    x_m = X_m + N2*u2
                    ge += w*vec((x_m-x_s)*Phi')
                end
            end

            # add contribution to contact virtual work
            sdofs = get_gdofs(problem, slave_element)
            mdofs = get_gdofs(problem, master_element)

            for i=1:field_dim
                lsdofs = sdofs[i:field_dim:end]
                lmdofs = mdofs[i:field_dim:end]
                add!(problem.assembly.C1, lsdofs, lsdofs, De)
                add!(problem.assembly.C1, lsdofs, lmdofs, -Me)
                add!(problem.assembly.C2, lsdofs, lsdofs, De)
                add!(problem.assembly.C2, lsdofs, lmdofs, -Me)
            end
            add!(problem.assembly.g, sdofs, ge)

        end # master elements done

    end # slave elements done, contact virtual work ready

end


text:
p1 = SVector{2,Float64}(0.0, 0.0)
p2 = SVector{2,Float64}(1.0, 0.0)
p3 = SVector{2,Float64}(1.0, 1.0)
p4 = SVector{2,Float64}(0.0, 1.0)
v = [p1, p2, p3, p4]

@testset "rotationally symmetric surface mesh creation" begin
    mesh_out = obj_from_point_sequence(v, 4)
    @test length(mesh_out.point) == 10
end


text:
module MyExample

greet() = print("Hello World!")

include("extra_file.jl")

export my_f

end


text:
#
println("Ca) Test of the Einstein module with ASF from an internally generated multiplet.")
#
setDefaults("print summary: open", "zzz-Einstein.sum")
setDefaults("unit: energy", "Kayser")
grid=JAC.Radial.Grid(true)
## grid = Radial.Grid(Radial.Grid(true), rnt = 2.0e-5,h = 1.0e-2, hp = 1.0e-2, NoPoints = 2000)
if  true
    asfSettings   = AsfSettings(AsfSettings(), scField=Basics.DFSField())
    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(26.), asfSettings=asfSettings,
                            configs=[Configuration("[Ne] 3s 3p^6"), Configuration("[Ne] 3s^2 3p^4 3d"), Configuration("[Ne] 3s^2 3p^5")],
                            propertySettings=[ Einstein.Settings([E1,M1,E2], true, LineSelection(false), 0., 0., 10000. ) ])

    wb = perform(wa; output=true)
    #
elseif  true
    ## grid = Radial.Grid(Radial.Grid(true), rnt = 4.0e-6, h = 1.0e-2, hp = 0., rbox = 2.0)
    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(36.),
                            configs=[Configuration("1s^2 2s"), Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s 2p^2")],
                            propertySettings=[ Einstein.Settings([E1], true, LineSelection(true, indexPairs=[(13,0), (15,0)]), 0., 0., 10000. )] )
                            ## Einstein.Settings([M2], true, LineSelection(true, indexPairs=[(13,2), (15,5), (15,4)]), 0., 0., 10000. ) )

    wb = perform(wa; output=true)
    #
elseif  false
    ## grid = Radial.Grid(Radial.Grid(true), rnt = 4.0e-6, h = 1.0e-2, hp = 0., rbox = 2.0)
    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(36.),
                            configs=[Configuration("1s^2 2s^2")],
                            ## configs=[Configuration("1s^2 2s"), Configuration("1s^2 2p"), Configuration("1s^2 3d"), Configuration("1s^2 4f")],
                            ##x configs=[Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s^2 2s")],
                            propertySettings=[ Einstein.Settings([E1, M1, E2], true, LineSelection(), 0., 0., 10000. ) ])
                            ##                 Einstein.Settings([M2], true, LineSelection(true, indexPairs=[(13,2), (15,5), (15,4)]), 0., 0., 10000. ) )

    wb         = perform(wa; output=true)
    frozenOrbs = wb["multiplet:"].levels[1].basis.orbitals
    #
elseif  true
    ## grid = Radial.Grid(Radial.Grid(true), rnt = 4.0e-6, h = 1.0e-2, hp = 0., rbox = 2.0)
    ## asfSettings   = AsfSettings(AsfSettings(), frozenSubshells=[Subshell("1s_1/2"), Subshell("2s_1/2"), Subshell("2p_1/2"), Subshell("2p_3/2"), Subshell("3d_3/2"), Subshell("3d_5/2")],
    ## asfSettings   = AsfSettings(AsfSettings(), frozenSubshells=[Subshell("1s_1/2"), Subshell("2s_1/2")],
    ##                                            startScfFrom=StartFromPrevious(frozenOrbs), generateScf=true)
    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(6.), ## asfSettings=asfSettings,
                            properties=[JAC.EinsteinX()],
                            ##x configs=[Configuration("1s^2 2s"), Configuration("1s^2 2p"), Configuration("1s^2 3d"), Configuration("1s^2 4f")],
                            configs=[Configuration("1s^2 2s^2"), Configuration("1s^2 2s 3p")],
                            einsteinSettings=Einstein.Settings([E1], true, LineSelection(), 0., 0., 10000. ) )
                            ## einsteinSettings=Einstein.Settings([E1,M1,E2], true, LineSelection(true, indexPairs=[(3,1), (5,1)]), 0., 0., 10000. ) )

    wb = perform(wa; output=true)
    #
end
setDefaults("print summary: close", "")


text:
module TensorNetAlgs

using Dates
using JLD2
using FileIO
#using Distributed

using MatrixProductStateTools
using SymTensors
using LinearAlgebra
using TensorOperations
using KrylovKit

#using SparseArrays

include("apply.jl")
include("dmrg.jl")
include("dmrg_envio.jl")
#include("idmrg.jl")
include("tdvp.jl")

export initialenv
export _initialenv_tempfiles

export dmrg1sitesweep!
export dmrg2sitesweep!
export dmrg2sitesweep_envio!
export dmrg2sitesweep_envasyncio!
export idmrg2site
export tdvp1sitesweep!
export tdvp2sitesweep!

end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## dummy stub for https://github.com/JuliaBinaryWrappers/Zlib_jll.jl
baremodule Zlib_jll
using Base, Libdl
Base.Experimental.@compiler_options compile=min optimize=0 infer=false

const PATH_list = String[]
const LIBPATH_list = String[]

export libz

# These get calculated in __init__()
const PATH = Ref("")
const LIBPATH = Ref("")
artifact_dir = ""
libz_handle = C_NULL
libz_path = ""

if Sys.iswindows()
    const libz = "libz.dll"
elseif Sys.isapple()
    const libz = "@rpath/libz.1.dylib"
else
    const libz = "libz.so.1"
end

function __init__()
    global libz_handle = dlopen(libz)
    global libz_path = dlpath(libz_handle)
    global artifact_dir = dirname(Sys.BINDIR)
    LIBPATH[] = dirname(libz_path)
    push!(LIBPATH_list, LIBPATH[])
end

# JLLWrappers API compatibility shims.  Note that not all of these will really make sense.
# For instance, `find_artifact_dir()` won't actually be the artifact directory, because
# there isn't one.  It instead returns the overall Julia prefix.
is_available() = true
find_artifact_dir() = artifact_dir
dev_jll() = error("stdlib JLLs cannot be dev'ed")
best_wrapper = nothing
get_libz_path() = libz_path

end  # module Zlib_jll

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
abstract type AbstractGPModel{T<:Real,L<:AbstractLikelihood,I<:AbstractInference,N} end

@traitdef IsFull{X}
@traitdef IsMultiOutput{X}
@traitdef IsSparse{X}

Base.eltype(::AbstractGPModel{T}) where {T} = T
data(m::AbstractGPModel) = m.data
likelihood(m::AbstractGPModel) = m.likelihood
inference(m::AbstractGPModel) = m.inference
n_latent(::AbstractGPModel{<:Real,<:AbstractLikelihood,<:AbstractInference,N}) where {N} = N
n_output(m::AbstractGPModel) = 1
@traitfn nFeatures(m::TGP) where {TGP <: AbstractGPModel; !IsSparse{TGP}} = nSamples(m)
@traitfn nFeatures(m::TGP) where {TGP <: AbstractGPModel; IsSparse{TGP}} = m.nFeatures

getf(m::AbstractGPModel) = m.f
getf(m::AbstractGPModel, i::Int) = m.f[i]
Base.getindex(m::AbstractGPModel, i::Int) = getf(m, i)

n_iter(m::AbstractGPModel) = n_iter(inference(m))

is_stochastic(m::AbstractGPModel) = is_stochastic(inference(m))

batchsize(m::AbstractGPModel) = batchsize(inference(m))

is_trained(m::AbstractGPModel) = m.trained
set_trained!(m::AbstractGPModel, status::Bool) = m.trained = status

verbose(m::AbstractGPModel) = m.verbose

post_step!(::AbstractGPModel, ::Any) = nothing

function Random.rand!(
    model::AbstractGPModel, A::DenseArray{T}, X::AbstractVector
) where {T<:Real}
    return rand!(MvNormal(predict_f(model, X; cov=true, diag=false)...), A)
end

function Random.rand(model::AbstractGPModel, X::AbstractMatrix, n::Int=1)
    return rand(model, KernelFunctions.RowVecs(X), n)
end

function Random.rand(model::AbstractGPModel{T}, X::AbstractVector, n::Int=1) where {T<:Real}
    if n_latent(model) == 1
        rand!(model, Array{T}(undef, length(X), n), X)
    else
        @error "Sampling not implemented for multiple output GPs"
    end
end

pr_covs(model::AbstractGPModel, x) = pr_cov.(model.f, Ref(x))

pr_means(model::AbstractGPModel) = pr_mean.(model.f)
pr_means(model::AbstractGPModel, X::AbstractVector) = pr_mean.(model.f, Ref(X))

setpr_means!(model::AbstractGPModel, pr_means) = setpr_mean!.(model.f, pr_means)

means(model::AbstractGPModel) = mean.(model.f)

covs(model::AbstractGPModel) = cov.(model.f)

kernels(model::AbstractGPModel) = kernel.(model.f)

setkernels!(model::AbstractGPModel, kernels) = setkernel!.(model.f, kernels)

setZs!(model::AbstractGPModel, Zs) = setZ!.(model.f, Zs)


text:
using Test, Random, LSHFunctions

include(joinpath("..", "utils.jl"))

#========================
Tests
========================#

@testset "Monte Carlo function hashing tests" begin
    Random.seed!(RANDOM_SEED)

    @testset "Construct MonteCarloHash" begin
        # Hash L^2([0,1]) over cosine similarity
        # Sampler μ() chooses a point in [0,1] uniformly at random
        μ() = rand()
        hashfn = MonteCarloHash(cossim, μ)

        @test n_hashes(hashfn) == 1
        @test hashfn.μ == μ
        @test similarity(hashfn) == cossim
        @test hashtype(hashfn) == hashtype(LSHFunction(cossim))

        # Hash L^1([0,1]) over L^1 distance
        hashfn = MonteCarloHash(L1, μ, 10)

        @test n_hashes(hashfn) == 10
        @test similarity(hashfn) == L1
        @test hashtype(hashfn) == hashtype(LSHFunction(ℓ1))
    end

    @testset "Pass invalid similarity function" begin
        # Attempt to construct a MonteCarloHash with an invalid
        # similarity function.
        μ() = rand()

        @test_throws(ErrorException, MonteCarloHash(ℓ1, μ))
        @test_throws(ErrorException, MonteCarloHash(ℓ2, μ))
        @test_throws(ErrorException, MonteCarloHash(jaccard, μ))

        # Construct a hash function (with valid similarity) in the same
        # manner in case the MonteCarloHash API ever changes.
        # This ensures that we won't forget to update these tests.
        _ = MonteCarloHash(L2, μ)
    end

    #==========
    Cosine similarity hashing
    ==========#
    @testset "Hash cosine similarity with trivial inputs" begin
        # Hash over cosine similarity between two functions with cosine similarity
        # zero. The collision rate should be close to 50%.
        f(x) = (0.0 ≤ x ≤ 0.5) ? 1.0 : 0.0;
        g(x) = (0.0 ≤ x ≤ 0.5) ? 0.0 : 1.0;
        hashfn = MonteCarloHash(cossim, rand, 1024)

        @test embedded_similarity(hashfn, f, g) == 0.0

        hf, hg = hashfn(f), hashfn(g)
        @test 0.45 ≤ mean(hf .== hg) ≤ 0.55
    end

    @testset "Hash cosine similarity with nontrivial inputs" begin
        # Test hashing on cosine similarity using step functions defined
        # over [0,N]. The functions are piecewise constant on the intervals
        # [i,i+1], i = 1, ..., N.
        N = 10
        μ() = N * rand()
        hashfn = MonteCarloHash(cossim, μ, 2048)

        @test let success = true, ii = 1
            while success && ii ≤ 128
                f, f_steps = create_step_function(N)
                g, g_steps = create_step_function(N)

                prob = collision_probability(hashfn, cossim(f_steps, g_steps); n_hashes=1)
                hf, hg = hashfn(f), hashfn(g)

                success &= (prob-0.05 ≤ mean(hf .== hg) ≤ prob+0.05)
                ii += 1
            end
            success
        end
    end

    #==========
    L^p distance hashing
    ==========#
    @testset "Hash L^p distance" begin
        # Use step functions to test distances between functions again, since
        # you can map them into R^N isomorphically
        N = 4
        μ() = N * rand()
        hashfn = MonteCarloHash(L1, μ, 1024; volume = N)

        @test let success = true, ii = 1
            while success && ii ≤ 128
                f, f_steps = create_step_function(N)
                g, g_steps = create_step_function(N)

                # Hash collision rate should be close to the probability of
                # collision
                prob = collision_probability(hashfn, ℓ1(f_steps, g_steps); n_hashes=1)
                hf, hg = hashfn(f), hashfn(g)

                success &= (prob-0.05 ≤ mean(hf .== hg) ≤ prob+0.05)
                ii += 1
            end
            success
        end
    end
end


text:
using Luxor, Test, Random
Random.seed!(42)
function circlecircletangent_test(fname)
    Drawing(800, 800, fname)
    origin()
    background("white")
    setopacity(0.7)
    sethue("black")
    # overlapping circles will return four Os
    c1center = c2center = O
    c1radius = 50
    c2radius = 100
    p1, p2, p3, p4 = circlecircleoutertangents(c1center, c1radius, c2center, c2radius)
    @test all(pt -> isequal(pt, O), (p1, p2, p3, p4)) == true
    for i in 1:30
        randomhue()
        c1center = rand(BoundingBox() * 0.5)
        c2center = rand(BoundingBox() * 0.75)
        c1radius = rand(5:30)
        c2radius = rand(30:50)
        circle(c1center, c1radius, :fill)
        circle(c2center, c2radius, :fill)
        pt1, pt2, pt3, pt4 = circlecircleoutertangents(c1center, c1radius, c2center, c2radius)
        poly([pt1, pt2, pt4, pt3], :fill)
        sethue("black")
        circle(c2center, c2radius, :stroke)
        circle(c1center, c1radius, :stroke)
    end
    @test finish() == true
    println("...finished test: output in $(fname)")
    preview()
end

function circlecircletangent_test2(fname)
    Drawing(800, 800, fname)
    origin()
    background("white")
    setline(0.5)
    sethue("black")
    # overlapping circles will return four Os
    c1center = c2center = O
    c1radius = 50
    c2radius = 100
    p1, p2, p3, p4 = circlecircleinnertangents(c1center, c1radius, c2center, c2radius)
    @test all(pt -> isequal(pt, O), (p1, p2, p3, p4)) == true

    for θ in range(0, step=2π/12, length=12)
        @layer begin
            randomhue()
            c1center = polar(100, θ)
            c2center = polar(300, θ)
            c1radius = 20
            c2radius = rand(50:70)
            circle(c1center, c1radius, :fill)
            circle(c2center, c2radius, :fill)

            pt1, pt2, pt3, pt4 = circlecircleinnertangents(c1center, c1radius, c2center, c2radius)

            if !all(pt -> isequal(pt, O), (pt1, pt2, pt3, pt4))
                sethue("green")
                circle(c2center, c2radius, :stroke)
                circle(c1center, c1radius, :stroke)
                circle.((pt1, pt2, pt3, pt4), 4, :fill)
                sethue("black")
                arrow(pt1, pt2, linewidth=2)
                arrow(pt3, pt4, linewidth=2)
            end

            pt1, pt2, pt3, pt4 = circlecircleoutertangents(c1center, c1radius, c2center, c2radius)

            if !all(pt -> isequal(pt, O), (pt1, pt2, pt3, pt4))
                sethue("red")
                arrow(pt1, pt2, linewidth=2)
                arrow(pt3, pt4, linewidth=2)
                circle.((pt1, pt2, pt3, pt4), 3, :fill)
            end
        end
    end
    @test finish() == true
    println("...finished test: output in $(fname)")
end

circlecircletangent_test("circlecircletangent-outer.png")
circlecircletangent_test2("circlecircletangent-inner.png")


text:
import Test: @test, @test_throws, @testset, @test_broken
import Luna: PhysData

@testset "All" begin
@testset "Exceptions" begin
    @test_throws DomainError PhysData.ref_index(:Hello, 800e-9)
end

@testset "refractive indices" begin
    @test PhysData.ref_index(:He, 800e-9) ≈ 1.000031950041203
    @test PhysData.ref_index(:He, 800e-9, 10) ≈ 1.0003180633169397
    @test PhysData.ref_index(:SiO2, 800e-9, lookup=false) ≈ 1.4533172548587419
    @test PhysData.ref_index(:SiO2, 400e-9, lookup=false) ≈ 1.4701161185594052
    @test PhysData.ref_index(:SiO2, PhysData.eV_to_m(6)) ≈ 1.543
    @test PhysData.ref_index(:SiO2, PhysData.eV_to_m(1.455)) ≈ 1.45248
    @test PhysData.ref_index(:SiO2, PhysData.eV_to_m(0.91018)) ≈ 1.44621
    @test real(PhysData.ref_index(:SiO2, PhysData.eV_to_m(121.6))) ≈ 0.9865
    @test imag(PhysData.ref_index(:SiO2, PhysData.eV_to_m(121.6))) ≈ 0.0085
end

@testset "Function equivalence" begin
    @test PhysData.ref_index_fun(:SiO2)(800e-9) == PhysData.ref_index(:SiO2, 800e-9)
    @test PhysData.ref_index_fun(:He)(800e-9) == PhysData.ref_index(:He, 800e-9)
end

@testset "Dispersion" begin
    @test PhysData.dispersion(2, :SiO2, 800e-9, lookup=false) ≈ 3.61619983e-26
    @test isapprox(PhysData.dispersion(2, :He, 800e-9), 9.373942337550116e-31, rtol=1e-5)
    @test isapprox(PhysData.dispersion(2, :He, 800e-9, 10), 9.33043805928079e-30, rtol=1e-5)
end

@testset "glasses" begin
    for g in PhysData.glass
        @test (g == :SiO2 ?
                !isreal(PhysData.ref_index(g, 800e-9)) : # SiO2 ref index is complex
                isreal(PhysData.ref_index(g, 800e-9)))
    end
    @test isreal(PhysData.ref_index(:SiO2, 800e-9, lookup=false))
end

@testset "gases" begin
    for g in PhysData.gas
        @test isreal(PhysData.ref_index(g, 800e-9))
        @test isreal(PhysData.ref_index(g, 200e-9))
        @test isreal(PhysData.ref_index(g, 800e-9, 10))
        @test isreal(PhysData.ref_index(g, 200e-9, 10))
    end
end

@testset "Nonlinear coefficients" begin
    @test_broken PhysData.χ3_gas(:He, 1) ≈ 1.2617371645226101e-27
    @test PhysData.χ3_gas(:Ar, 1) ≈ 2.964158749949189e-26
    @test_broken PhysData.n2_gas(:He, 1) ≈ 3.5647819877255427e-25
    @test PhysData.n2_gas(:He, 2) ≈ 7.125642138007481e-25
    @test_broken PhysData.n2_gas.(:He, [1, 2]) ≈ [3.5647819877255427e-25, 7.125642138007481e-25]
    @test_broken PhysData.n2_gas.([:He, :Ne], 1) ≈ [3.5647819877255427e-25, 6.416061508801999e-25]
    for gas in PhysData.gas[2:end] # Don't have γ3 for Air
        @test isreal(PhysData.n2_gas(gas, 1))
    end
end

@testset "Density" begin
    # compare to refprop
    # room temperature
    @test isapprox(PhysData.density(:Ar, 0.002, 294.0), 4.927169444649429e22, rtol=2e-15)
    @test isapprox(PhysData.density(:Ar, 0.02, 294.0), 4.9272289924653735e23, rtol=9e-15)
    @test isapprox(PhysData.density(:Ar, 2.0, 294.0), 4.933761614600933e25, rtol=4e-13)
    @test isapprox(PhysData.density(:Ar, 40.0, 294.0), 1.0101579129300146e27, rtol=2e-11)
    @test isapprox(PhysData.density(:Ar, 400.0, 294.0), 9.270757850984163e27, rtol=3e-10)
    @test_broken isapprox(PhysData.density(:He, 0.002, 294.0), 4.927180563885407e22, rtol=7e-16)
    @test_broken isapprox(PhysData.density(:He, 0.02, 294.0), 4.927137517453137e23, rtol=2e-16)
    @test_broken isapprox(PhysData.density(:He, 2.0, 294.0), 4.9224080868066745e25, rtol=2e-13)
    @test_broken isapprox(PhysData.density(:He, 40.0, 294.0), 9.66755645602771e26, rtol=6e-16)
    @test_broken isapprox(PhysData.density(:He, 400.0, 294.0), 8.309132317978155e27, rtol=1e-13)
    @test isapprox(PhysData.density(:Xe, 0.002, 294.0), 4.927238737907382e22, rtol=7e-16)
    @test isapprox(PhysData.density(:Xe, 0.02, 294.0), 4.92771934762324e23, rtol=2e-16)
    @test isapprox(PhysData.density(:Xe, 2.0, 294.0), 4.981603290885485e25, rtol=9e-16)
    @test isapprox(PhysData.density(:Xe, 20.0, 294.0), 5.5912074541820115e26, rtol=8e-15)
    @test isapprox(PhysData.density(:Xe, 40.0, 294.0), 1.3512153640360566e27, rtol=7e-16)
    @test isapprox(PhysData.density(:Xe, 60.0, 294.0), 3.218616771878294e27, rtol=6e-15)
    @test isapprox(PhysData.density(:N2, 0.002, 294.0), 4.927165137695148e22, rtol=7e-14)
    @test isapprox(PhysData.density(:N2, 0.02, 294.0), 4.927185920476881e23, rtol=7e-13)
    @test isapprox(PhysData.density(:N2, 2.0, 294.0), 4.929427604602592e25, rtol=6e-11)
    @test isapprox(PhysData.density(:N2, 40.0, 294.0), 9.909239438516337e26, rtol=3e-10)
    @test isapprox(PhysData.density(:N2, 400.0, 294.0), 7.890787516580131e27, rtol=1e-7)
    @test isapprox(PhysData.density(:H2, 0.002, 294.0), 4.9271795187225075e22, rtol=9e-16)
    @test isapprox(PhysData.density(:H2, 0.02, 294.0), 4.927127065719356e23, rtol=6e-16)
    @test isapprox(PhysData.density(:H2, 2.0, 294.0), 4.921361734871252e25, rtol=4e-13)
    @test isapprox(PhysData.density(:H2, 40.0, 294.0), 9.624525613497868e26, rtol=6e-16)
    @test isapprox(PhysData.density(:H2, 400.0, 294.0), 7.840662987579036e27, rtol=6e-16)
    # temperature range
    @test_broken isapprox(PhysData.density(:He, 10.0, 10.0), 9.218927214187151e27, rtol=6e-16)
    @test_broken isapprox(PhysData.density(:He, 10.0, 100.0), 7.142310177179197e26, rtol=6e-16)
    @test_broken isapprox(PhysData.density(:He, 10.0, 1000.0), 7.234660811823096e25, rtol=6e-13)
    @test isapprox(PhysData.density(:Xe, 10.0, 170.0), 1.336261228046876e28, rtol=6e-16)
    @test isapprox(PhysData.density(:Xe, 10.0, 700.0), 1.0361035050644844e26, rtol=6e-13)
    @test isapprox(PhysData.density(:N2, 10.0, 100.0), 1.4849593021049305e28, rtol=6e-7)
    @test isapprox(PhysData.density(:N2, 10.0, 1000.0), 7.217962021789807e25, rtol=6e-10)
end

@testset "Density spline" begin
    P = range(0, 10, length=128)
    Plow = range(0, 0.1, length=128)
    Pfine = range(0, 0.01, length=128)
    for g in PhysData.gas
        dens = PhysData.densityspline(g, Pmax=maximum(P))
        @test all(dens.(P) .≈ PhysData.density.(g, P))
        dens = PhysData.densityspline(g, Pmax=maximum(Plow))
        @test all(dens.(Plow) .≈ PhysData.density.(g, Plow))
        @test all(isapprox.(dens.(Pfine), PhysData.density.(g, Pfine), rtol=1e-6))
    end
end
end


text:
export wrand, wrandn, wrandnp, wsobol, sobol, soboln, sobolnp, wgrid, xgrid, near

##### random

"""
    wrand(d, k) = Weighted(rand(d, k))
    wrand(T, d, k)

Uniformly distributed vectors in `[0,1]^d`, as columns of a `Weighted{Matrix}`
which knows to `clamp` them to this box. `T = Float64` by default.

Keyword `weights=true` gives weights ∝ `1 .+ rand(k)` rather than constant.

    wrandn(d, k) = scale .* Weighted(randn(d, k))

Normally distributed `d`-vectors, of mean zero and std. dev. `scale=1` by default.
Keyword `max=10` clamps absolute values to be less than this.

    wrandnp(d, k)

Absolute value of normally distributed...
"""
wrand(d::Int, k::Int=1; kw...) = wrand(Float64,d,k; kw...)
wrand(T::Type, d::Int, k::Int=1; weights=true) = Weighted(rand(T,d,k), weights ? 1 .+ rand(T,k) : ones(T,k) ,0,1)

@doc @doc(wrand)
wrandn(d::Int, k::Int=1; kw...) = wrandn(Float64,d,k; kw...)
function wrandn(T::Type, d::Int, k::Int=1; scale=one(T), weights=true, max=typemax(T))
    arr = clamp.(scale .* randn(T,d,k), -max, max)
    opt = max==Inf ? WeightOpt() : WeightOpt(clamp=true, lo=-max, hi=max)
    Weighted(arr, normalise(weights ? 1 .+ rand(T,k) : ones(T,k)), opt)
end

@doc @doc(wrand)
wrandnp(d::Int, k::Int=1; kw...) = wrandnp(Float64,d,k; kw...)
function wrandnp(T::Type, d::Int, k::Int=1; scale=one(T), weights=true, max=typemax(T))
    arr = clamp.(abs.(scale .* randn(T,d,k)), 0, max)
    opt = WeightOpt(clamp=true, lo=0, hi=max)
    Weighted(arr, normalise(weights ? 1 .+ rand(T,k) : ones(T,k)), opt)
end

for f in [:wrand, :wrandn, :wrandnp]
    @eval function $f(::Union{typeof(big), Type{BigFloat}}, d::Int, k::Int=1; kw...)
        r64 = $f(d,k;kw...)
        Weighted(big.(r64.array), big.(r64.weights), r64.opt) |> normalise!
    end
end

##### sub-random

import Sobol

"""
    sobol(d, k)

First `k` entries from the `d`-dimensional Sobol sequence, as columns of a `Weighted{Matrix}`.
"""
function sobol(d::Int, k::Int=1)
    s = Sobol.SobolSeq(d)
    cols = reduce(hcat, [Sobol.next!(s) for i=1:k])
    Weighted(cols, ones(k), 0,1)
end

using SpecialFunctions: erfinv

"""
    soboln(d, k)

Normally distributed `d`-vectors, of mean zero and std. dev. 1, generated from Sobol sequence,
as columns of a `Weighted{Matrix}`.
"""
function soboln(d::Int, k::Int=1; scale=1)
    boxed = sobol(d,k)
    Weighted(scale .* erfinv.( -1 .+ 2 .* array(boxed)), boxed.weights)
end

function sobolnp(d::Int, k::Int=1; scale=π)
    boxed = sobol(d,k)
    Weighted(scale .* erfinv.(boxed.array), boxed.weights, 0, Inf)
end

sobol(T::Type, d::Int, k::Int=1) = Weighted(map(T, sobol(d,k).array), ones(T,k), 0,1)
soboln(T::Type, d::Int, k::Int=1) = Weighted(map(T, soboln(d,k).array), ones(T,k))
sobolnp(T::Type, d::Int, k::Int=1) = Weighted(map(T, sobolnp(d,k).array), ones(T,k))

##### ordered!

wgrid(d::Int, range::AbstractRange) = Weighted(xgrid(d,range), extrema(range)...)

##### not Weighted

"""
    xgrid(d, 0:0.1:5)

Gives the matrix whose colums are `d`-vectors, forming a grid of the given range in all dimensions.
"""
xgrid(d::Int, n::Int) = xgrid(d, 0:n)

function (xgrid(d::Int, range::AbstractRange{T}; big=false)::Matrix{T}) where T
    d==1 && return Matrix(collect(range)')

    n = length(range)
    !big && n^d > 4_000_000 && @error "xgrid is happiest with n^d = < 4×10^6, you have $n^$d... keyword big=true will disable this check"

    intcols = zeros(Int,d,n^d)
    @inbounds for i=1:n^d
        @views digits!(intcols[:,i], i-1; base=n)
    end

    range[1] .+ step(range) .* intcols
end

function near(x::Matrix, y, dist::Real; verbose=false)
    yesno = minimum( pairwise2(x,array(y)) ,dims=2) .<= (dist^2)
    sum(yesno)==0 && begin @warn("no nearby points, returning complete matrix"); return x end
    verbose && @info string("kept ",sum(yesno)," / ",size(x,2)," columns: ",round(100*sum(yesno)/size(x,2),2)," percent")
    x[:,vec(yesno)]
end

using Distances

"""
    pairwise2(x, y=x) = Distances.pairwise(SqEuclidean(), x, y)

Resulting `mat[i,j]` is distance sqared from `x[:,i]` to `y[:,j]`. Implementation varies.
"""
pairwise2(x::Matrix, y::Matrix) = pairwise(SqEuclidean(), x, y, dims=2)
pairwise2(x::Matrix) = pairwise(SqEuclidean(), x)

pairwise2(x::AbsMat, y::AbsMat) = diag(x'*x) .+ diag(y'*y)' .- 2x'*y
function pairwise2(x::AbsMat) # = diag(x'*x) .+ diag(x'*x)' .- 2x'*x
    mat = x'*x
    vec = diag(mat)
    vec .+ vec' .- 2 .* mat
end