text:
struct TypeID{typeid} end
TypeID(typeid) = TypeID{Symbol(typeid)}()
show(io::IO, ::TypeID{typeid}) where {typeid} = print(io, "Genser.TypeID(", string(typeid), ")")
print(io::IO, ::TypeID{typeid}) where {typeid} = print(io, typeid)

# Main to converter
function togenser(::Type{GenserDataType}, v::V) where V
    T = gensertypefor(V)
    togenser(T, v)
end

# Main from converter
function fromgenser(::Type{T}, v::GenserDataType) where T
    throw(ArgumentError("unable to convert from $(typeof(v)) to $(T)"))
end

# Default convert_to_type
function convert_to_type(type, v::GenserDataType) where T
    if type isa TypeID
        fromgenser(Any, v)
    else
        fromgenser(type, v)
    end
end

togenser(::Type{T}, v::T) where T <: GenserDataType = v

Base.convert(::Type{T}, v::T) where T <:GenserDataType = v
Base.convert(::Type{T}, v) where T <:GenserDataType = togenser(T, v)
Base.convert(::Type{T}, v::V) where T where V <:GenserDataType = fromgenser(T, v)
Base.convert(::Type{Any}, v::V) where V <:GenserDataType = fromgenser(Any, v)

togenser(v) = togenser(GenserDataType, v)
fromgenser(v::V) where V <: GenserDataType = fromgenser(Any, v)
fromgenser(::Type{T}, v::GenserValue{T}) where T = v.value

# Any
fromgenser(::Type{Any}, v::GenserAny) = fromgenser(v.value)

# Nothing
togenser(T::Type{<:GenserNothingValue}, v) = T()
fromgenser(::Type{>:Nothing}, v::GenserNothingValue) = nothing
fromgenser(::Type{Any}, v::GenserNothingValue) = nothing

# Integer types
for (T, GT) = (
    Int8 => GenserInt8,
    UInt8 => GenserUInt8,
    Int16 => GenserInt16,
    UInt16 => GenserUInt16,
    Int32 => GenserInt32,
    UInt32 => GenserUInt32,
    Int64 => GenserInt64,
    UInt64 => GenserUInt64,
    Int128 => GenserInt128,
    UInt128 => GenserUInt128,
    Bool => GenserBool,
    BigInt => GenserBigInt,
    Float16 => GenserFloat16,
    Float32 => GenserFloat32,
    Float64 => GenserFloat64,
    BigFloat => GenserBigFloat,
    Rational => GenserRational,
)

    @eval togenser(::Type{$GT}, v::$GT) = v
    @eval togenser(::Type{$GT}, v::V) where V <: Number = $GT($T(v))
    @eval togenser(::Type{$GT}, v::V) where V <: AbstractChar = $GT($T(v))
    if T == Rational
        @eval togenser(::Type{$GT}, v::V) where V <: AbstractString = $GT($T(parse(Float64, v)))
    else
        @eval togenser(::Type{$GT}, v::V) where V <: AbstractString = $GT(parse($T, v))
    end

    @eval fromgenser(T::Type{<:Number}, v::$GT) = T(v.value)
    @eval fromgenser(T::Type{Char}, v::$GT) = T(v.value)
    @eval fromgenser(T::Type{>:$T}, v::$GT) = v.value
end

# Char
togenser(::Type{GenserChar}, v::GenserChar) = v
togenser(::Type{GenserChar}, v::Char) = GenserChar(v)
togenser(::Type{GenserChar}, v::V) where V <: Number = GenserChar(Char(v))
togenser(::Type{GenserChar}, v::V) where V <: AbstractString = begin
    @assert length(v) == 1 "cannot convert a multi-character or empty string to a char"
    GenserChar(v[1])
end

fromgenser(::Type{Char}, v::GenserChar) = v.value
fromgenser(::Type{AbstractString}, v::GenserChar) = Base.string(v.value)
fromgenser(::Type{String}, v::GenserChar) = Base.string(v.value)
fromgenser(T::Type{<:Number}, v::GenserChar) = T(v.value)
fromgenser(::Type{>:Char}, v::GenserChar) = v.value

# Strings
togenser(T::Type{<:GenserStringValue}, v) = T(string(v))

fromgenser(::Type{AbstractString}, v::V) where V <: GenserStringValue = v.value
fromgenser(T::Type{<:Number}, v::V) where V <:GenserStringValue = parse(T, v.value)
fromgenser(T::Type{<:AbstractChar}, v::V) where V <: GenserStringValue = begin
    @assert length(v.value) == 1 "cannot convert a multi-character or empty string to a $T"
    T(v.value[1])
end
fromgenser(::Type{Symbol}, v::V) where V <:GenserStringValue = Symbol(v.value)
fromgenser(::Type{>:AbstractString}, v::V) where V <: GenserStringValue = v.value
fromgenser(::Type{>:String}, v::V) where V <: GenserStringValue = string(v.value)
fromgenser(T::Type{<:AbstractString}, v::V) where V <: GenserStringValue = begin
    T = T.isconcretetype ? T : typeof(v.value)
    T(v.value)
end
fromgenser(::Type{Base.UUID}, v::V) where V <: GenserStringValue = UUID(v.value)
fromgenser(::Type{Vector{UInt8}}, v::V) where V <: GenserStringValue = Vector{UInt8}(v.value)

# Binary
togenser(T::Type{<:GenserBinaryValue}, v::Vector{UInt8}) = T(v)
togenser(T::Type{<:GenserBinaryValue}, v::Array) = begin
    v = map(UInt8, v)
    T(v)
end
togenser(::Type{GenserBinary}, v::AbstractString) = begin
    v = Vector{UInt8}(v)
    GenserBinary(v)
end
togenser(::Type{GenserBinaryValue{E}}, v::AbstractString) where E = begin
    C = genser_converter_for_encoding(E)
    v = C.decode(v)
    GenserBinaryValue{E}(v)
end

fromgenser(::Type{Vector{UInt8}}, v::V) where {V <: GenserBinaryValue} = v.value
fromgenser(T::Type{<:Vector}, v::V) where {V <: GenserBinaryValue} = begin
    v = map(eltype(T), (v.value))
    v
end
fromgenser(T::Type{<:Matrix}, v::V) where {V <: GenserBinaryValue} = begin
    v = map(eltype(T), (v.value))
    # Unknown dimensions...
    hcat(v)
end
fromgenser(::Type{AbstractString}, v::GenserBinary) = String(v.value)
fromgenser(T::Type{<:AbstractString}, v::GenserBinary) = T(v.value)
fromgenser(::Type{AbstractString}, v::GenserBinaryValue{E}) where E = begin
    C = genser_converter_for_encoding(E)
    v = C.encode(v.value)
    v
end
fromgenser(T::Type{<:AbstractString}, v::GenserBinaryValue{E}) where E = begin
    C = genser_converter_for_encoding(E)
    v = C.encode(v.value)
    T(v)
end
fromgenser(::Type{UUID}, v::V) where {V <: GenserBinaryValue} = begin
    @assert length(v.value) == 16 "cannot convert a binary to a uuid when its length is not  16 bytes"
    v = ntoh(reinterpret(UInt128, v.value)[1])
    UUID(v)
end
fromgenser(::Type{>:Vector{UInt8}}, v::V) where {V <: GenserBinaryValue} = v.value

# UUID
togenser(::Type{GenserUUID}, v::UUID) = GenserUUID(v)
togenser(::Type{GenserUUID}, v::AbstractString) = GenserUUID(UUID(v))
togenser(::Type{GenserUUID}, v::Vector{UInt8}) = begin
    @assert length(v) == 16 "cannot convert a binary to a uuid when its length is not 16 bytes"
    v = hton(reinterpret(UInt128, v)[1])
    GenserUUID(UUID(v))
end
togenser(::Type{GenserUUID}, v::UInt128) = GenserUUID(UUID(v))

# fromgenser(::Type{String}, v::GenserUUID) = string(v.value)
fromgenser(T::Type{AbstractString}, v::GenserUUID) = string(v.value)
fromgenser(T::Type{<:AbstractString}, v::GenserUUID) = T(string(v.value))
fromgenser(::Type{T}, v::GenserUUID) where T <: AbstractString = T(string(v.value))
fromgenser(::Type{Vector{UInt8}}, v::GenserUUID) = [reinterpret(UInt8, [hton(v.value.value)])...]
fromgenser(::Type{UInt128}, v::GenserUUID) = v.value.value
fromgenser(::Type{>:UUID}, v::GenserUUID) = v.value
fromgenser(::Type{<:UUID}, v::GenserUUID) = v.value

# Symbol
togenser(::Type{GenserSymbol}, v::Symbol) = GenserSymbol(v)
togenser(::Type{GenserSymbol}, v) = GenserSymbol(Symbol(v))

fromgenser(::Type{String}, v::GenserSymbol) = String(v.value)
fromgenser(::Type{AbstractString}, v::GenserSymbol) = String(v.value)
fromgenser(T::Type{<:AbstractString}, v::GenserSymbol) = T(String(v.value))
fromgenser(::Type{>:Symbol}, v::GenserSymbol) = v.value

# Sequences
togenser(::Type{GenserSequence{T}}, vs::AbstractVector{T}) where T <: GenserDataType = GenserSequence{T}(vs)
togenser(::Type{GenserSequence{T}}, vs::AbstractVector) where T <: GenserDataType = begin
    newvs = map(vs) do v
        togenser(T, v)
    end

    GenserSequence{T}(newvs)
end
togenser(T::Type{<: GenserSequence}, vs::AbstractArray) = togenser(T, vs[:])
togenser(T::Type{<: GenserSequence}, v::AbstractString) = togenser(T, [v...])
togenser(T::Type{GenserSequence}, vs) = begin
    vs = Vector(vs)
    togenser(T, vs)
end

fromgenser(T::Type{<: AbstractVector}, v::GenserSequence) = begin
    E = eltype(T)
    vs = map(v.value) do subv
        fromgenser(E, subv)
    end

    T(vs)
end

fromgenser(T::Type{>: Vector{E}}, v::GenserSequence) where E = begin
    vs = map(v.value) do subv
        fromgenser(E, subv)
    end

    vs::T
end

fromgenser(::Type{Any}, v::GenserSequence) = begin
    vs = map(v.value) do subv
        fromgenser(Any, subv)
    end

    vs
end

fromgenser(T::Type{<: AbstractSet}, v::GenserSequence) = begin
    E = eltype(T)
    vs = map(v.value) do subv
        fromgenser(E, subv)
    end

    T = T.isconcretetype ? T : Set
    T(vs)
end

convert_to_type(m, v::GenserSequence) = begin
    map(v.value) do subv
        convert_to_type(m, subv)
    end
end

# Sets
togenser(::Type{GenserSet{T}}, vs::AbstractSet{T}) where T <: GenserDataType = GenserSet{T}(vs)
togenser(::Type{GenserSet{T}}, vs::AbstractSet) where T <: GenserDataType = begin
    newvs = map([vs...]) do v
        togenser(T, v)
    end

    GenserSet{T}(Set(newvs))
end
togenser(T::Type{<:GenserSet}, v) = togenser(T, Set(v))

fromgenser(T::Type{Any}, v::GenserSet) = begin
    vs = map([v.value...]) do subv
        fromgenser(Any, subv)
    end

    Set(vs)
end

fromgenser(T::Type{<: AbstractSet}, v::GenserSet) = begin
    E = eltype(T)
    vs = map([v.value...]) do subv
        fromgenser(E, subv)
    end

    T = T.isconcretetype ? T : Set
    T(vs)
end

fromgenser(T::Type{<: AbstractVector}, v::GenserSet) = begin
    E = eltype(T)
    vs = map([v.value...]) do subv
        fromgenser(E, subv)
    end

    T(vs)
end

fromgenser(T::Type{>: Set{E}}, v::GenserSet) where E = begin
    vs = map([v.value...]) do subv
        fromgenser(E, subv)
    end
    Set(vs)::T
end

convert_to_type(m, v::GenserSet) = begin
    newv = map([v.value...]) do subv
        convert_to_type(m, subv)
    end

    Set(newv)
end

# Tuples
togenser(::Type{GenserTuple{T}}, v::T) where T = GenserTuple{T}(v)
togenser(::Type{GenserTuple{T}}, v) where T = GenserTuple{T}(T(v))

fromgenser(T::Type{<:Tuple}, v::GenserTuple) = begin
    @assert length(v.value) >= length(T.types)

    vals = []
    for (pos, E) in pairs(T.types)
        push!(vals, fromgenser(E, v.value[pos]))
    end

    T(vals)
end

fromgenser(T::Type{<:AbstractVector}, v::GenserTuple) = begin
    E = eltype(T)
    v = Base.map(v.value) do subv
        fromgenser(E, subv)
    end

    T([v...])
end

convert_to_type(m, v::GenserTuple) = begin
    v = Base.map(v.value) do subv
        convert_to_type(m, subv)
    end

    v
end

fromgenser(T::Type{Any}, v::GenserTuple) = begin
    v = Base.map(v.value) do subv
        fromgenser(Any, subv)
    end

    v
end

# Dicts
togenser(::Type{GenserDict{K,V}}, vs::AbstractDict{K,V}) where K <: GenserDataType where V <: GenserDataType = GenserDict{K,V}(vs)
togenser(::Type{GenserDict{K,V}}, vs::AbstractDict) where K <: GenserDataType where V <: GenserDataType = begin
    pairs = []
    for (k,v) in vs
        k = togenser(K, k)
        v = togenser(V, v)
        push!(pairs, k => v)
    end
    newvs = Dict(pairs)
    GenserDict{K,V}(newvs)
end

togenser(::Type{GenserDict{K,V}}, vs::T) where {K,V,T} = begin
    pairs = []
    for k = fieldnames(T)
        v = getfield(vs, k)
        k = togenser(K, k)
        v = togenser(V, v)
        push!(pairs, k => v)
    end
    newvs = Dict(pairs)
    GenserDict{K,V}(newvs)
end

fromgenser(T::Type{<:AbstractDict}, vs::GenserDict{K,V}) where K where V = begin
    TK = keytype(T)
    TV = valtype(T)
    pairs = []
    for (k,v) in vs.value
        k = fromgenser(TK, k)
        v = fromgenser(TV, v)
        push!(pairs, k => v)
    end
    T(pairs)
end

fromgenser(T::Type{>:Dict{DK, DV}}, vs::GenserDict{K,V}) where {DK, DV, K, V} = begin
    pairs = []
    for (k,v) in vs.value
        k = fromgenser(DK, k)
        v = fromgenser(DV, v)
        push!(pairs, k => v)
    end
    Dict(pairs)
end

fromgenser(T::Type{<:NamedTuple}, vs::GenserDict{K,V}) where {K, V} = begin
    args = []
    for (key, type) = zip(fieldnames(T), T.types)
        gkey = togenser(K, key)
        if haskey(vs.value, gkey)
            push!(args, fromgenser(type, vs.value[gkey]))
        else
            # Attempt to convert from nothing
            push!(args, fromgenser(type, GenserUndefined()))
        end
    end
    T(args)
end

fromgenser(T::Type{Any}, vs::GenserDict{K,V}) where {K, V} = begin
    pairs = []
    for (k,v) in vs.value
        k = fromgenser(Any, k)
        v = fromgenser(Any, v)
        push!(pairs, k => v)
    end
    Dict(pairs)
end

fromgenser(::Type{T}, v::GenserDict{K,V}) where {T,K,V} = begin
    args = []
    for (key, type) = zip(fieldnames(T), T.types)
        gkey = togenser(K, key)
        if haskey(v.value, gkey)
            push!(args, fromgenser(type, v.value[gkey]))
        else
            # Attempt to convert from nothing
            push!(args, fromgenser(type, GenserUndefined()))
        end
    end

    # TODO: Need a way to set an instantiation strategy for the type
    T(args...)
end

convert_to_type(m, vs::GenserDict{K,V}) where {K, V} = begin
    pairs = []
    for (k,v) in vs.value
        k = convert_to_type(m, k)
        v = convert_to_type(m, v)
        push!(pairs, k => v)
    end
    Dict(pairs)
end

# Records
togenser(::Type{GenserRecord{T}}, v::AbstractDict) where {T <: NamedTuple} = begin
    args = []
    for (key, type) = zip(fieldnames(T), T.types)
        if haskey(v, key)
            push!(args, togenser(type, v[key]))
        elseif haskey(v, string(key))
            push!(args, togenser(type, v[string(key)]))
        else
            # Attempt to convert from nothing
            push!(args, togenser(type, nothing))
        end
    end

    GenserRecord(T(args))
end
togenser(::Type{GenserRecord{T}}, v::T) where T = GenserRecord{T}(v)
togenser(::Type{GenserRecord{T}}, v::V) where T where V = begin
    args = []
    for (key, type) = zip(fieldnames(T), T.types)
        if hasfield(V, key)
            push!(args, togenser(type, getfield(v, key)))
        else
            # Attempt to convert from nothing
            push!(args, togenser(type, nothing))
        end
    end

    GenserRecord(T(args))
end

fromgenser(T::Type{<:NamedTuple}, v::GenserRecord) = begin
    args = []
    for (key, type) = zip(fieldnames(T), T.types)
        if hasfield(typeof(v.value), key)
            push!(args, fromgenser(type, getfield(v.value, key)))
        else
            # Attempt to convert from nothing
            push!(args, fromgenser(type, GenserUndefined()))
        end
    end

    T(args)
end

fromgenser(T::Type{<:AbstractDict}, v::GenserRecord) = begin
    TK = typeof(T) == UnionAll ? Symbol : keytype(T)
    TV = typeof(T) == UnionAll ? Any : valtype(T)
    entries = []
    for (k,v) = Base.pairs(v.value)
        k = fromgenser(TK, GenserSymbol(k))
        v = fromgenser(TV, v)
        push!(entries, k => v)
    end
    T = T === AbstractDict ? Dict : T
    T(entries)
end

fromgenser(T::Type{Any}, v::GenserRecord) = begin
    Base.map(v.value) do subv
        fromgenser(Any, subv)
    end
end

fromgenser(::Type{T}, v::GenserRecord) where T = begin
    args = []
    for (key, type) = zip(fieldnames(T), T.types)
        if hasfield(typeof(v.value), key)
            push!(args, fromgenser(type, getfield(v.value, key)))
        else
            # Attempt to convert from nothing
            push!(args, fromgenser(type, GenserUndefined()))
        end
    end

    # TODO: Need a way to set an instantiation strategy for the type
    T(args...)
end

convert_to_type(m, v::GenserRecord) = begin
    Base.map(v.value) do subv
        convert_to_type(m, subv)
    end
end

# Optionals
togenser(::Type{GenserOptional{T}}, v::T) where T = GenserOptional{T}(v)
togenser(::Type{GenserOptional{T}}, v::Nothing) where T = GenserOptional{T}(GenserUndefined())
togenser(::Type{GenserOptional{T}}, v) where T = GenserOptional{T}(togenser(T, v))

fromgenser(T::Type{>: Nothing}, v::V) where V <: GenserOptional = begin
    if v.value isa GenserNothingValue
        nothing
    elseif typeof(T) == Union
        fromgenser(unoptionalize(T), v.value)
    else
        fromgenser(T, v.value)
    end
end

convert_to_type(m, v::GenserOptional) = convert_to_type(m, v.value)

# Variants
togenser(::Type{GenserVariant{T}}, v::T) where T = GenserVariant{T}(v)
togenser(::Type{GenserVariant{T}}, v::V) where {T, V} = begin
    if V <: T
        return GenserVariant{T}(v)
    end

    GenserVariant{T}(togenser(v))
end

fromgenser(::Type{T}, v::V) where T where V <: GenserVariant = begin
    # TODO: Fix this
    fromgenser(T, v.value)
end
convert_to_type(m, v::GenserVariant) = convert_to_type(m, v.value)

# Any
togenser(::Type{GenserAny}, v::T) where T <: GenserDataType = GenserAny(v)
togenser(::Type{GenserAny}, v::T) where T = GenserAny(togenser(v))
fromgenser(::Type{T}, v::GenserAny) where T = fromgenser(T, v.value)


text:
using Compat
using ITensors
using Test
using LinearAlgebra

@testset "Threading" begin
  blas_num_threads = Compat.get_num_threads()
  strided_num_threads = ITensors.NDTensors.Strided.get_num_threads()

  BLAS.set_num_threads(1)
  ITensors.NDTensors.Strided.set_num_threads(1)

  @testset "Threaded contraction" begin
    i = Index([QN(0) => 500, QN(1) => 500])
    A = randomITensor(i', dag(i))

    using_threaded_blocksparse = ITensors.disable_threaded_blocksparse()

    R = A' * A

    ITensors.enable_threaded_blocksparse()

    Rthreaded = A' * A

    @test R ≈ Rthreaded

    #ITensors.disable_threaded_blocksparse()
    #time = @elapsed B = A' * A
    #ITensors.enable_threaded_blocksparse()
    #time_threaded = @elapsed B = A' * A
    #@test time > time_threaded

    if !using_threaded_blocksparse
      ITensors.disable_threaded_blocksparse()
    end
  end

  @testset "Contraction resulting in no blocks with threading bug" begin
    i = Index([QN(0) => 1, QN(1) => 1])
    A = emptyITensor(i', dag(i))
    B = emptyITensor(i', dag(i))
    A[i' => 1, i => 1] = 11.0
    B[i' => 2, i => 2] = 22.0

    using_threaded_blocksparse = ITensors.disable_threaded_blocksparse()
    C1 = A' * B
    ITensors.enable_threaded_blocksparse()
    C2 = A' * B
    if !using_threaded_blocksparse
      ITensors.disable_threaded_blocksparse()
    end

    @test nnzblocks(C1) == 0
    @test nnzblocks(C2) == 0
    @test nnz(C1) == 0
    @test nnz(C2) == 0
    @test C1 ≈ C2
  end

  @testset "Bug fixed in threaded block sparse" begin
    maxdim = 10
    nsweeps = 2
    outputlevel = 0
    cutoff = 0.0
    Nx = 4
    Ny = 2
    U = 4.0
    t = 1.0
    N = Nx * Ny
    sweeps = Sweeps(nsweeps)
    maxdims = min.(maxdim, [100, 200, 400, 800, 2000, 3000, maxdim])
    maxdim!(sweeps, maxdims...)
    cutoff!(sweeps, cutoff)
    noise!(sweeps, 1e-6, 1e-7, 1e-8, 0.0)
    sites = siteinds("Electron", N; conserve_qns=true)
    lattice = square_lattice(Nx, Ny; yperiodic=true)
    ampo = OpSum()
    for b in lattice
      ampo .+= -t, "Cdagup", b.s1, "Cup", b.s2
      ampo .+= -t, "Cdagup", b.s2, "Cup", b.s1
      ampo .+= -t, "Cdagdn", b.s1, "Cdn", b.s2
      ampo .+= -t, "Cdagdn", b.s2, "Cdn", b.s1
    end
    for n in 1:N
      ampo .+= U, "Nupdn", n
    end
    H = MPO(ampo, sites)
    Hsplit = splitblocks(linkinds, H)
    state = [isodd(n) ? "↑" : "↓" for n in 1:N]
    ψ0 = productMPS(sites, state)
    using_threaded_blocksparse = ITensors.enable_threaded_blocksparse()
    energy, _ = dmrg(H, ψ0, sweeps; outputlevel=outputlevel)
    energy_split, _ = dmrg(Hsplit, ψ0, sweeps; outputlevel=outputlevel)
    @test energy_split ≈ energy
    if !using_threaded_blocksparse
      ITensors.disable_threaded_blocksparse()
    end
  end

  BLAS.set_num_threads(blas_num_threads)
  ITensors.NDTensors.Strided.set_num_threads(strided_num_threads)
end

ext: jl


text:
"""
## rod_mm!

This subroutine forms the consistent mass matrix of a 1-d "rod" fin_el.

### Function
```julia
rod_mm!(km, ell)
```

### Arguments
```julia
* mm::Matrix{Float64}       : Element mass matrix (Updated)
* ell::Float64              : Element length
```
"""
function rod_mm!(mm::Matrix{Float64}, ell::Float64)
  #
  # This subroutine forms the consistent mass matrix of a 1-d "rod" fin_el.
  #
 mm = [1.0/3.0 1.0/6.0; 1.0/6.0 1.0/3.0]
 mm *= ell
end


text:
#
# Correctness Tests
#

using RDatasets
using Test

my_tests = ["packages.jl",
            "datasets.jl",
            "dataset.jl"]

@testset "$my_test" for my_test in my_tests
    @info "Running $(my_test)..."
    include(my_test)
end

max_stars_repo_name: blockchainhelppro/Wadren-Analysis

max_issues_repo_name: blockchainhelppro/Wadren-Analysis

max_forks_repo_name: blockchainhelppro/Wadren-Analysis


text:
"""
```
transform_data(m::AbstractDSGEModel, levels::DataFrame; cond_type::Symbol = :none,
    verbose::Symbol = :low)
```

Transform data loaded in levels and order columns appropriately for the DSGE
model. Returns DataFrame of transformed data.

The DataFrame `levels` is output from `load_data_levels`. The series in levels are
transformed as specified in `m.observable_mappings`.

- To prepare for per-capita transformations, population data are filtered using
  `hpfilter`. The series in `levels` to use as the population series is given by
  the `population_mnemonic` setting. If `use_population_forecast(m)`, a
  population forecast is appended to the recorded population levels before the
  filtering. Both filtered and unfiltered population levels and growth rates are
  added to the `levels` data frame.
- The transformations are applied for each series using the `levels` DataFrame
  as input.

Conditional data (identified by `cond_type in [:semi, :full]`) are handled
slightly differently: If `use_population_forecast(m)`, we drop the first period
of the population forecast because we treat the first forecast period
`date_forecast_start(m)` as if it were data. We also only apply transformations
for the observables given in `cond_full_names(m)` or `cond_semi_names(m)`.
"""
function transform_data(m::AbstractDSGEModel, levels::DataFrame;
                        cond_type::Symbol = :none, verbose::Symbol = :low)

    n_obs, _ = size(levels)

    # Step 1: HP filter (including population forecasts, if they're being used)
    population_mnemonic = parse_population_mnemonic(m)[1]
    if !isnull(population_mnemonic)
        population_forecast_levels = if use_population_forecast(m)
            read_population_forecast(m; verbose = verbose)
        else
            DataFrame()
        end

        population_data, _ = transform_population_data(levels, population_forecast_levels,
                                                       get(population_mnemonic);
                                                       verbose = verbose,
                                                       use_hpfilter = hpfilter_population(m))

        levels = join(levels, population_data, on = :date, kind = :left)
        name_maps = [s => t for (s,t) = zip([:filtered_population_recorded, :dlfiltered_population_recorded, :dlpopulation_recorded],
                [:filtered_population, :filtered_population_growth, :unfiltered_population_growth])]
        rename!(levels, name_maps)
    end

    # Step 2: apply transformations to each series
    transformed = DataFrame()
    transformed[!, :date] = levels[:, :date]

    data_transforms = collect_data_transforms(m)

    for series in keys(data_transforms)
        println(verbose, :high, "Transforming series $series...")
        f = data_transforms[series]
        transformed[!, series] = f(levels)
    end

    sort!(transformed, [:date])

    return transformed
end

function collect_data_transforms(m; direction=:fwd)

    data_transforms = OrderedDict{Symbol,Function}()

    # Parse vector of observable mappings into data_transforms dictionary
    for obs in keys(m.observable_mappings)
        data_transforms[obs] = getfield(m.observable_mappings[obs], Symbol(string(direction) * "_transform"))
    end

    data_transforms
end

"""
```
transform_population_data(population_data, population_forecast,
    population_mnemonic; verbose = :low)
```

Load, HP-filter, and compute growth rates from population data in
levels. Optionally do the same for forecasts.

### Inputs

- `population_data`: pre-loaded DataFrame of historical population data
  containing the columns `:date` and `population_mnemonic`. Assumes this is
  sorted by date.
- `population_forecast`: pre-loaded `DataFrame` of population forecast
  containing the columns `:date` and `population_mnemonic`
- `population_mnemonic`: column name for population series in `population_data`
  and `population_forecast`

### Keyword Arguments

- `verbose`: one of `:none`, `:low`, or `:high`
- `use_hpfilter`: whether to HP filter population data and forecast. See `Output` below.
- `pad_forecast_start::Bool`: Whether you want to re-size
the population_forecast such that the first index is one quarter ahead of the last index
of population_data. Only set to false if you have manually constructed population_forecast
to artificially start a quarter earlier, so as to avoid having an unnecessary missing first entry.

### Output

Two dictionaries containing the following keys:

- `population_data_out`:
  + `:filtered_population_recorded`: HP-filtered historical population series (levels)
  + `:dlfiltered_population_recorded`: HP-filtered historical population series (growth rates)
  + `:dlpopulation_recorded`: Non-filtered historical population series (growth rates)

- `population_forecast_out`:
  + `:filtered_population_forecast`: HP-filtered population forecast series (levels)
  + `:dlfiltered_population_forecast`: HP-filtered population forecast series (growth rates)
  + `:dlpopulation_forecast`: Non-filtered historical population series (growth rates)

If `population_forecast_file` is not provided, the r\"*forecast\" fields will be
empty. If `use_hpfilter = false`, then the r\"*filtered*\" fields will be
empty.
"""
function transform_population_data(population_data::DataFrame, population_forecast::DataFrame,
                                   population_mnemonic::Symbol; verbose = :low,
                                   use_hpfilter::Bool = true,
                                   pad_forecast_start::Bool = false)

    # Unfiltered population data
    population_recorded = population_data[:, [:date, population_mnemonic]]

    # Make sure first period of unfiltered population forecast is the first forecast quarter
    if !isempty(population_forecast) && !pad_forecast_start
        last_recorded_date = population_recorded[end, :date]
        if population_forecast[1, :date] <= last_recorded_date
            last_recorded_ind   = findall(population_forecast[!,:date] .== last_recorded_date)[1]
            population_forecast = population_forecast[(last_recorded_ind+1):end, :]
        end
        @assert subtract_quarters(population_forecast[1, :date], last_recorded_date) == 1
    end

    # population_all: full unfiltered series (including forecast)
    population_all = if isempty(population_forecast)
        population_recorded[!,population_mnemonic]
    else
        pop_all = vcat(population_recorded, population_forecast)
        pop_all[!,population_mnemonic]
    end

    # HP filter
    if use_hpfilter
        population_all = convert(Array{Union{Float64, Missing}}, population_all)
        filtered_population, _ = hpfilter(population_all, 1600)
    end

    # Output dictionary for population data
    population_data_out = DataFrame()
    population_data_out[!,:date] = convert(Array{Date}, population_recorded[!,:date])
    population_data_out[!,:dlpopulation_recorded] = difflog(population_recorded[!,population_mnemonic])

    n_population_forecast_obs = size(population_forecast,1)

    if use_hpfilter
        filt_pop_recorded = filtered_population[1:end-n_population_forecast_obs]
        population_data_out[!,:filtered_population_recorded] = filt_pop_recorded
        population_data_out[!,:dlfiltered_population_recorded] = difflog(filt_pop_recorded)
    end

    # Output dictionary for population forecast
    population_forecast_out = DataFrame()
    if n_population_forecast_obs > 0
        population_forecast_out[!,:date] = convert(Array{Date}, population_forecast[!,:date])
        population_forecast_out[!,:dlpopulation_forecast] = difflog(population_forecast[!,population_mnemonic])

        if use_hpfilter
            filt_pop_fcast = filtered_population[end-n_population_forecast_obs:end]
            population_forecast_out[!,:filtered_population_forecast]   = filt_pop_fcast[2:end]
            population_forecast_out[!,:dlfiltered_population_forecast] = difflog(filt_pop_fcast)[2:end]
        end
    end

    return population_data_out, population_forecast_out
end


text:
# # Tutorial 1
#
# ## Geometry basics, elementary entities, physical groups
#
#
# The Julia API is entirely defined in the `gmsh.jl` module (which contains the
# full documentation of all the functions in the API):
import gmsh
#
# Before using any functions in the Julia API, Gmsh must be initialized:
gmsh.initialize()
#
# Next we add a new model named "t1" (if `gmsh.model.add()` is not called a new
# unnamed model will be created on the fly, if necessary):
gmsh.model.add("t1")
#
# The Julia API provides direct access to each supported geometry (CAD)
# kernel. The built-in kernel is used in this first tutorial: the corresponding
# API functions have the `gmsh.model.geo` prefix.
#
# The first type of `elementary entity` in Gmsh is a `Point`. To create a point
# with the built-in CAD kernel, the Julia API function is
# gmsh.model.geo.addPoint():
# - the first 3 arguments are the point coordinates (x, y, z)
# - the next (optional) argument is the target mesh size close to the point
# - the last (optional) argument is the point tag (a stricly positive integer
#   that uniquely identifies the point)
lc = 1e-2
gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
#
# The distribution of the mesh element sizes will be obtained by interpolation
# of these mesh sizes throughout the geometry. Another method to specify mesh
# sizes is to use general mesh size Fields (see `t10.jl`).
# If no target mesh size of provided, a default uniform coarse size will be used
# for the model, based on the overall model size.
#
# We can then define some additional points. All points should have different
# tags:
gmsh.model.geo.addPoint(.1, 0,  0, lc, 2)
gmsh.model.geo.addPoint(.1, .3, 0, lc, 3)
#
# If the tag is not provided explicitly, a new tag is automatically created, and
# returned by the function:
p4 = gmsh.model.geo.addPoint(0, .3, 0, lc)
#
# Curves are Gmsh's second type of elementery entities, and, amongst curves,
# straight lines are the simplest. The API to create straight line segments with
# the built-in kernel follows the same conventions: the first 2 arguments are
# point tags (the start and end points of the line), and the last (optional one)
# is the line tag.
#
# In the commands below, for example, the line 1 starts at point 1 and ends at
# point 2.
#
# Note that curve tags are separate from point tags - hence we can reuse tag `1`
# for our first curve. And as a general rule, elementary entity tags in Gmsh
# have to be unique per geometrical dimension.
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(3, 2, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)
#
# The third elementary entity is the surface. In order to define a simple
# rectangular surface from the four curves defined above, a curve loop has first
# to be defined. A curve loop is defined by an ordered list of connected curves,
# a sign being associated with each curve (depending on the orientation of the
# curve to form a loop). The API function to create curve loops takes a list
# of integers as first argument, and the curve loop tag (which must be unique
# amongst curve loops) as the second (optional) argument:
gmsh.model.geo.addCurveLoop([4, 1, -2, 3], 1)
#
# We can then define the surface as a list of curve loops (only one here,
# representing the external contour, since there are no holes--see `t4.jl` for
# an example of a surface with a hole):
gmsh.model.geo.addPlaneSurface([1], 1)
#
# Before they can be meshed (and, more generally, before they can be used by API
# functions outside of the built-in CAD kernel functions), the CAD entities must
# be synchronized with the Gmsh model, which will create the relevant Gmsh data
# structures. This is achieved by the gmsh.model.geo.synchronize() API call for
# the built-in CAD kernel. Synchronizations can be called at any time, but they
# involve a non trivial amount of processing; so while you could synchronize the
# internal CAD data after every CAD command, it is usually better to minimize
# the number of synchronization points.
gmsh.model.geo.synchronize()
#
# At this level, Gmsh knows everything to display the rectangular surface 1 and
# to mesh it. An optional step is needed if we want to group elementary
# geometrical entities into more meaningful groups, e.g. to define some
# mathematical ("domain", "boundary"), functional ("left wing", "fuselage") or
# material ("steel", "carbon") properties.
#
# Such groups are called "Physical Groups" in Gmsh. By default, if physical
# groups are defined, Gmsh will export in output files only mesh elements that
# belong to at least one physical group. (To force Gmsh to save all elements,
# whether they belong to physical groups or not, set the `Mesh.SaveAll' option
# to 1.) Physical groups are also identified by tags, i.e. stricly positive
# integers, that should be unique per dimension (0D, 1D, 2D or 3D). Physical
# groups can also be given names.
#
# Here we define a physical curve that groups the left, bottom and right curves
# in a single group (with prescribed tag 5); and a physical surface with name
# "My surface" (with an automatic tag) containing the geometrical surface 1:
gmsh.model.addPhysicalGroup(1, [1, 2, 4], 5)
ps = gmsh.model.addPhysicalGroup(2, [1])
gmsh.model.setPhysicalName(2, ps, "My surface")
#
# We can then generate a 2D mesh
gmsh.model.mesh.generate(2)
#
# ... and save it to disk
gmsh.write("t1.msh")
#
# Remember that by default, if physical groups are defined, Gmsh will export in
# the output mesh file only those elements that belong to at least one physical
# group. To force Gmsh to save all elements, you can use
#
# `gmsh.option.setNumber("Mesh.SaveAll", 1)`
#
# By default, Gmsh saves meshes in the latest version of the Gmsh mesh file
# format (the `MSH` format). You can save meshes in other mesh formats by
# specifying a filename with a different extension. For example
#
# `gmsh.write("t1.unv")`
#
# will save the mesh in the UNV format. You can also save the mesh in older
# versions of the MSH format: simply set
#
# `gmsh.option.setNumber("Mesh.MshFileVersion", x)`
#
# for any version number `x`. As an alternative, you can also not specify the
# format explicitly, and just choose a filename with the `.msh2` or `.msh4`
# extension.
#
# Note that starting with Gmsh 3.0, models can be built using other geometry
# kernels than the default "built-in" kernel. To use the OpenCASCADE CAD kernel
# instead of the built-in kernel, you should use the functions with the
# `gmsh.model.occ` prefix.
#
# Different CAD kernels have different features. With OpenCASCADE, instead of
# defining the surface by successively defining 4 points, 4 curves and 1 curve
# loop, one can define the rectangular surface directly with
#
# `gmsh.model.occ.addRectangle(.2, 0, 0, .1, .3)`
#
# After synchronization with the Gmsh model with
#
# `gmsh.model.occ.synchronize()`
#
# the underlying curves and points could be accessed with
# `gmsh.model.getBoundary()`.
#
# This should be called when you are done using the Gmsh Julia API:
gmsh.finalize()


text:
function mark_kross_diag!(r)

    function inverse_side_2(side, k)
        return HorizonSide((Int(side) + k) % 4)
    end

    function mark_diag(s)
        k = 0
        while (!isborder(r, s) && !isborder(r, inverse_side_2(s, 1)))
            move!(r, s)
            move!(r, inverse_side_2(s, 1))
            putmarker!(r)
            k += 1
        end
        for i in 1:k
            move!(r, inverse_side_2(s, 2))
            move!(r, inverse_side_2(s, 3))
        end
    end

    for i in [Nord, Ost, Sud, West]
        mark_diag(i)
    end
    putmarker!(r)
end

max_stars_repo_name: dmitrii-a-piskunov/mirea-progs

max_issues_repo_name: dmitrii-a-piskunov/mirea-progs

max_forks_repo_name: dmitrii-a-piskunov/mirea-progs


text:
using Test, LinearAlgebra, BenchmarkTools, Parameters, Revise

N = 64
U0 = randn(N,N,N)

U1 = similar(U0)
U2 = similar(U0)

x = LinRange(-1,1,N)
y = x'

@time z = reshape(x,1,1,length(x))

V(x,y,z,t) = x^2 + y^2 + z^2

@time @. U0 += V(x,y,z,0.)

@time @. U0 = V(x,y,z,0.)

@time @. U0 += V(x,y,z,0.)

@time @. U0 = U0 + U1

function Vadd!(V1,V0,V,t)
    @. V1 = V0 + V(x,y,z,t)
    return nothing
end

@time Vadd!(U1,U0,V,0.)

# function Valloc!(V1,V,t)
#     @. V1 = V(x,y,z,t)
#     return nothing
# end

# @btime Valloc!(V1,V,0.)

# macro

using FourierGPE
import FourierGPE.V

# static potential array
struct VStatic{N} <: UserParams
    V0::Array{Complex{Float64},N}
end

V0(x,y,z) = 0.25*(x^2 + y^2 + 4*z^2)
p0 = VStatic(V0.(x,y,z) .|> complex)

# dynamical potential function
@with_kw mutable struct Potential <: UserParams
    V::Expr = :( 0.5*(x^2 + y^2 + 4*z^2) )
end

macro potential(V,p)
    @eval $p = Potential($V)
end

macro staticpotential(V,p)
    @eval $p = VStatic($V.(x,y,z) .|> complex)
end

#TODO change sim to have V0, V, Vt fields
macro potentials(sim,p,q)
    @eval Sim($sim,V0 = $p,V = $q)
end

pot = Potential()
ex = pot.V

@potential :(0.25*(x^2 + y^2 + 4*z^2)) p1

@staticpotential V0 p2

#TODO V0 should just be a regular array
@test p2 == p0

p0 .≈ p2


text:
# Copyright 2020 Paul Melis (paul.melis@surf.nl)
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Printf
using StaticArrays

include("halfedge.jl")

# Some utility functions to handle Julia's 1-based indexing
function get_vertex(V, idx)
    s = 1+3*(idx-1)
    return SVector(V[s], V[s+1], V[s+2])
end

function set_vertex(V, idx, v)
    s = 1+3*(idx-1)
    V[s:s+2] = v
end

function time_subdivide(vertices::Array, loop_start::Array, loop_total::Array, loops::Array)
    @timev subdivide(vertices, loop_start, loop_total, loops)
end

function subdivide(vertices::Array, loop_start::Array, loop_total::Array, loops::Array)
    #println(vertices)
    #println(loop_start)
    #println(loop_total)
    #println(loops)

    t0 = time()

    num_vertices, num_faces, num_edges, vertex_start_edges, face_start_edges, polygon_edges = build(vertices, loop_start, loop_total, loops)

    #println(num_vertices)
    #println(num_faces)
    #println(num_edges)
    #println(vertex_start_edges)
    #println(face_start_edges)
    #println(polygon_edges)

    t1 = time()
    @printf("(Julia) Building half edges done in %.3fms\n", 1000*(t1-t0))

    println("(Julia) Input: $(num_vertices) vertices, $(num_faces) polygons, $(num_edges) polygon edges")

    # One new vertex for each input face, one new vertex for each edge
    output_num_vertices = num_vertices + num_faces + num_edges
    # 1 .. NV           Original input vertices (initially, are overwritten later on)
    # NV+1 .. NV+NF     New face points
    # NV+NF+1 .. end    New edge points
    output_vertices = Array{Float32}(undef, 3*output_num_vertices)
    # Copy original input vertex positions, to be modified later on
    output_vertices[1:3*num_vertices] = vertices

    # Output is always all quads, as each input face is split into n quads,
    # where n is the number of vertices in the face
    output_num_quads = sum(loop_total)
    output_loop_start = collect(range(UInt32(1), step=UInt32(4), stop=UInt32(output_num_quads*4)))
    output_loop_total = fill(UInt32(4), output_num_quads)
    output_loops = Array{UInt32}(undef, 4*output_num_quads)

    # Wether a vertex is on the boundary
    is_boundary_vertex = zeros(Bool, num_vertices)

    println("(Julia) Output: $(output_num_vertices) vertices, $(output_num_quads) quads")

    function face_point_index(fi)
        return num_vertices + fi
    end

    function edge_point_index(ei)
        return num_vertices + num_faces + ei
    end

    #
    # Subdivide
    #

    # Add new face points: average of existing original face vertices
    for fi = 1:num_faces
        sum = SVector{3, Float32}(0.0, 0.0, 0.0)
        n = 0

        he = start = face_start_edges[fi]
        while true
            sum += get_vertex(vertices, he.source)
            n += 1

            he = he.next
            if he == start break end
        end

        #println("face point ", fi, " ", sum/n)

        set_vertex(output_vertices, face_point_index(fi), sum / n)

    end

    #println("FV ", face_vertices)
    #println("OV ", output_vertices)

    # Add new edge points: average of edge endpoints and neighbouring face points

    for ei = 1:num_edges

        he = polygon_edges[ei]

        if he.sibling != nothing

            edge_point = (
                get_vertex(output_vertices, face_point_index(he.face))
                +
                get_vertex(output_vertices, face_point_index(he.sibling.face))
                +
                get_vertex(output_vertices, he.source)
                +
                get_vertex(output_vertices, he.target)
            ) / 4

        else

            # Boundary rule
            edge_point = 0.5*(get_vertex(output_vertices, he.source) + get_vertex(output_vertices, he.target))

            is_boundary_vertex[he.source] = is_boundary_vertex[he.target] = true

        end

        #println("edge point ", ei, " ", edge_point)

        set_vertex(output_vertices, edge_point_index(ei), edge_point)
    end

    # Move original input vertices to new positions

    for vi = 1:num_vertices

        if !haskey(vertex_start_edges, vi)
            # Skip unconnected vertices
            continue
        end

        P = get_vertex(output_vertices, vi)

        if is_boundary_vertex[vi]

            he = vertex_start_edges[vi]

            Pp = (
                get_vertex(output_vertices, he.target)
                +
                6*P
                +
                get_vertex(output_vertices, he.prev.source)
            ) / 8

            set_vertex(output_vertices, vi, Pp)

        else

            F_sum = SVector{3, Float32}(0, 0, 0)
            R_sum = SVector{3, Float32}(0, 0, 0)
            n = 0

            # R = average of edge midpoints
            #   = 1/n * sum_i(0.5*(P+Q_i))
            #   = 1/n * (0.5*n*P + 0.5*sum_i(Q_i))
            #   = 0.5*P + 1/2n * sum_i(Q_i)

            he = start = vertex_start_edges[vi]
            while true
                F_sum += get_vertex(output_vertices, face_point_index(he.face))
                R_sum += get_vertex(output_vertices, he.target)
                n += 1

                @assert he.sibling != nothing

                he = he.sibling.next
                if he == start break end
            end

            F = F_sum / n
            R = R_sum / (2*n) + 0.5f0 * P

            set_vertex(output_vertices, vi, (F + 2*R + (n-3)*P) / n)

        end
    end

    # Create new face loops for the quads for each subdivided original face

    offs = 1
    ofi = 1

    for fi = 1:num_faces

        he = start = face_start_edges[fi]
        while true
            output_loop_start[ofi] = offs
            output_loops[offs] = edge_point_index(he.index)
            output_loops[offs+1] = he.target
            output_loops[offs+2] = edge_point_index(he.next.index)
            output_loops[offs+3] = face_point_index(he.face)
            offs += 4
            ofi += 1

            he = he.next
            if he == start break end
        end
    end

    t2 = time()
    @printf("(Julia) Subdivision done in %.3fms\n", 1000*(t2-t1))

    return output_vertices, output_loop_start, output_loop_total, output_loops
end

max_stars_repo_name: paulmelis/blender-julia-test

max_issues_repo_name: paulmelis/blender-julia-test

max_forks_repo_name: paulmelis/blender-julia-test


text:
using TensorTrains
export stsolve, inexstsolve, htsolve

function stsolve(L::Int64, D::Int64, residual::Function,
            normest::Float64, invnormest::Float64,
            ɛ::Float64, Θ::Float64 = 0.25,
            ux::Tensor = Vector{Tensor3}(), αx::Float64 = 0.,
            maxiter::Int64 = typemax(Int64))

  μ = 2 / (1/invnormest + normest);
  κ = normest * invnormest;
  ρ = (κ - 1.) / (κ + 1.);

  if length(ux) > 0
    u = ux
  else
    u = Tensor(2^D, L, 0)
  end
  res = norm!(residual(u))

  udiff = (1. - ρ)/(ρ*normest) * res
      # ensures that α is not decreased in first iteration

  if αx > 0.
    α = αx
  else
    α = .5 * μ * res
  end

  println("κ = ", κ, ", ρ = ", ρ, ", μ = ", μ)

  reslist = Float64[]

  i = 0
  while res > ɛ && i < maxiter
    i += 1
    u0 = deepcopy(u)

    r = residual(u)

    println("   res ranks ", ranks(r))
    println("   max ", maximum(ranks(r)))

    res = norm!(r)
    push!(reslist, res)

    if res <= ɛ
      break
    elseif udiff <= .95 * (1-ρ)/(ρ*normest) * res
      α *= Θ
    end

    add!(u, -μ, r)
    softthresh!(u, α)

    add!(u0, -1., u)
    udiff = norm!(u0)

    println(i, ": α = ", α, ", res ", res, ", udiff ", udiff,
          "  max rank ", maximum(ranks(u)))
  end

  return u, α, reslist
end

function inexstsolve(L::Int64, D::Int64, residual::Function,
            normest::Float64, invnormest::Float64,
            ɛ::Float64, Θ::Float64 = 0.25, ν::Float64 = .95)
  μ = 2 / (1/invnormest + normest);
	κ = normest * invnormest;
	ρ = (κ - 1.) / (κ + 1.);

  τ1 = .5*(1. - ρ)/(3. - ρ)
  τ2 = .25*(1. - ρ)

  ρ̂ = (ρ + τ2*(1. + ρ)/(1. - τ2))

  u = Tensor(2^D, L, 0)
  r0, _ = residual(u, 0.)
  res = norm!(r0)

  α = .5 * μ * res

  B = (1. - ρ)*(1. - τ1)*ν / (normest * (1. + τ2) * ρ̂)
  M = min( (1. - τ1)*τ2*B / (μ * (1 + τ1 + normest*B)),
      (1. - τ1)^2 * ν*ρ*τ2 / (μ * (ρ*(1. + τ1)*(1. + τ2) + ν*(1. - τ1)*(1. - ρ))) )

  udiff = (B / ν) * res

  println("κ = ", κ, ", ρ = ", ρ, ", μ = ", μ)

  reslist = Float64[]
  uranklist = Int64[]
  rranklist = Int64[]
  iranklist = Int64[]

  η = τ1 * res
  i = 0

  while true
    i += 1
    u0 = deepcopy(u)

    r, imax = residual(u, η)
    res = norm!(r)

    println("   η = ", η, ": max res rank ",
            maximum(ranks(r)), "  internal max ", imax)

    while η > τ1 * res
      η *= .5
      r, imax = residual(u, η)
      res = norm!(r)

      println("   (1) η = ", η, ": max res rank ",
              maximum(ranks(r)), "  internal max ", imax)
    end

    push!(reslist, res + η)
    push!(uranklist, maximum(ranks(u)))
    push!(rranklist, maximum(ranks(r)))
    push!(iranklist, maximum(imax))

    if res + η <= ɛ
      break
    elseif udiff <= B * res
      α *= Θ
      η = τ1 * res
    end

    add!(u, -μ, r)
    softthresh!(u, α)

    diff = add(1., u0, -1., u)
    udiff = norm!(diff)

    println(i, ": α = ", α, ", res ", res, ", udiff ", udiff,
      "  max rank ", maximum(ranks(u)))

    while η > (τ2 / μ) * udiff && η > M * res
      η *= .5
      r, imax = residual(u, η)
      res = norm!(r)

      println("   (2) η = ", η, ": max res rank ",
              maximum(ranks(r)), "  internal max ", imax)

      add!(u, -μ, r)
      softthresh!(u, α)

      diff = add(1., u0, -1., u)
      udiff = norm!(diff)

      println(i, ": α = ", α, ", res ", res, ", udiff ", udiff,
        "  max rank ", maximum(ranks(u)))
    end

  end
  return u, (reslist, uranklist, rranklist, iranklist)
end

function htsolve(L::Int64, D::Int64, residual::Function,
      normest::Real, invnormest::Real, ɛ0::Real,
      u = Tensor(2^D, L, 0), κP::Real = sqrt(L-1.), δ0::Real = 0.)
  α = .05
  θ = .5
  β = .5

  κ1 = 1/(1 + κP*(1+α))
  κ2 = κP*(1+α)*κ1

  ω = 2 / (1/invnormest + normest)
  ρ = (normest*invnormest - 1.) / (normest*invnormest + 1.)

  J = 0
  while ρ^J * (1. + (ω + β)*J) > θ*κ1
    J += 1
  end

  if δ0 > 0.
    ε = δ0
  else
    ε = invnormest * norm(residual(u, 0.)[1])
  end

  reslist = Vector{Float64}[]
  uranklist = Int64[]
  rranklist = Int64[]
  iranklist = Int64[]

  k = 0
  while ε > ɛ0
    k += 1

    for j = 0:J
      η = ρ^(j+1) * ε

      r, mr = residual(u, η)
      rnorm = norm(r)

      push!(reslist, rnorm + η)
      push!(uranklist, maximum(ranks(u)))
      push!(rranklist, maximum(ranks(r)))
      push!(iranklist, maximum(imax))

      println(k, ", ", j, ": ", rnorm + η, "  (", mr, ")")

      if invnormest*(rnorm+β*η) <= θ*κ1*ε
        break
      end

      add!(u, -ω, r)
      svdtrunc!(u, svd!(u), β*η)

      println("   ", ranks(u)')
    end

    ε *= θ
    svdtrunc!(u, svd!(u), κ2*ε)
  end

  return u, (reslist, uranklist, rranklist, iranklist)
end


text:
using Pkg

handle("ping") do
  "pong"
end

handle("args") do args
  if length(filter!(!isempty, args)) !== 0
    append!(Base.ARGS, args)
  end
end

handle("echo") do data
  data
end

handle("cd") do path
  hideprompt() do
    cd(path)
    @info "Changed working directory to `$path`"
  end
end

handle("activateProject") do dir
  hideprompt() do
    Pkg.activate(dir)
  end
end

handle("activateParentProject") do dir
  hideprompt() do
    if (path = find_project_file(dir)) === nothing
      @warn "No project file found for `$dir`"
      return
    end
    Pkg.activate(path)
  end
end

function find_project_file(dir)
  while true
    next_dir = dirname(dir)
    (
      next_dir == dir || # ensure to escape infinite recursion
      isempty(dir)       # reached to the system root
    ) && return nothing
    path = joinpath(dir, "Project.toml")
    isfile(path) && return path
    dir = next_dir
  end
end

handle("activateDefaultProject") do
  hideprompt() do
    Pkg.activate()
  end
end

handle("evalsimple") do code
  Core.eval(Main, Meta.parse(code))
end

handle("exit") do
  exit()
end

handle("packages") do
  finddevpackages()
end

const PlotPaneEnabled = Ref(true)

handle("enableplotpane") do enable
  if enable
    Media.setdisplay(Media.Graphical, PlotPane())
    PlotPaneEnabled[] = true
  else
    Media.unsetdisplay(Media.Graphical)
    PlotPaneEnabled[] = false
  end
  nothing
end

handle("cancelCallback") do args...
  # TODO: Do something sensible here.
  # Until then it's better to silently fail than spam the REPL with "unrecognised
  # message" warnings.
end

using Pkg: status
handle("reportinfo") do
  io = IOBuffer()
  versioninfo(io)
  println(io)

  old = stdout
  rd, wr = redirect_stdout()
  Pkg.status()
  redirect_stdout(old)
  close(wr)
  println(io, String(read(rd)))

  String(take!(io))
end


text:
using IntervalCensored
using Test

@testset "IntervalCensored.jl" begin
    # Write your tests here.
end

max_stars_repo_name: abikoushi/IntervalCensored

max_issues_repo_name: abikoushi/IntervalCensored

max_forks_repo_name: abikoushi/IntervalCensored


text:
export selector, clearconsole, input, structure, @profiler

"""
    selector([xs...]) -> x

Allow the user to select one of the `xs`.

`xs` should be an iterator of strings. Currently there is no fallback in other
environments.
"""
selector(xs) = Main.Atom.selector(xs)

"""
    clearconsole()

Clear the console if Juno is used; does nothing otherwise.
"""
clearconsole() = isactive() && Main.Atom.clearconsole()

"""
    input(prompt = "") -> "..."

Prompt the user to input some text, and return it. Optionally display a prompt.
"""
input(prompt = "") = (print(prompt); isactive() ? Main.Atom.input() : readline())

"""
    notify(msg::AbstractString)

Display `msg` as an OS specific notification.

Useful for signaling the end of a long running computation or similar. This
disregards the `Notifications` setting in `julia-client`. Falls back to
`@info(msg)` in other environments.
"""
notify(msg::AbstractString) = isactive() ? Main.Atom.sendnotify(msg) : @info(msg)

"""
    notification(message::AbstractString; kind = :Info, options = Dict())

Adds an notification via Atom's builtin notification system.
`message` and `options` serves same as described in [`NotificationManager` API](https://flight-manual.atom.io/api/v1.46.0/NotificationManager/#instance-addSuccess),
while `kind` specifies which kind of notification will be used as follows:
```js
atom.notifications[`add\${kind}`](message, options)
```

Falls back to `@info(message)` in other environments.
"""
function notification(message::AbstractString; kind = :Info, options = Dict())
  if isactive()
    Main.Atom.sendnotification(message; kind = kind, options = options)
  else
    kind = string(kind)
    occursin(r"warn"i, kind) ? @warn(msg) :
    occursin(r"error"i, kind) ? @error(msg) :
    @info(msg)
  end
end

"""
    plotsize()

Get the size of Juno's plot pane in `px`. Returns `[100, 100]` as a fallback value.
"""
plotsize() = isactive() ? Main.Atom.plotsize() : [100, 100]

"""
    syntaxcolors(selectors = Atom.SELECTORS)::Dict{String, UInt32}

Get the colors used by the current Atom theme.
`selectors` should be a `Dict{String, Vector{String}}` which assigns a css
selector (e.g. `syntax--julia`) to a name (e.g. `variable`).
"""
syntaxcolors(selectors) = isactive() ? Main.Atom.syntaxcolors(selectors) : Dict{String, UInt32}()
syntaxcolors() = isactive() ? Main.Atom.syntaxcolors() : Dict{String, UInt32}()

"""
    profiler(data=Profile.fetch(); lidict=nothing, C=false, combine=true, recur=:off, pruned=[])

Show profile information as an in-editor flamechart.
Any keyword argument that [`FlameGraphs.flamegraph`](@ref) can take could be given.
"""
profiler(args...; kwargs...) = isactive() && Main.Atom.Profiler.profiler(args...; kwargs...)

"""
    @profiler exp [kwargs...]

Clear currently collected profile traces, profile the provided expression and show
  it via `Juno.profiler()`.

Any keyword argument that [`FlameGraphs.flamegraph`](@ref) can take could be given
  as optional arugments `kwargs...`

```julia
# profile a function call
@profiler fname(fargs)

# include ccalls and compress recursive calls
@profiler fname(fargs) C = true recur = :flat
```
"""
macro profiler(exp, kwargs...)
  quote
    let
      $(Profile).clear()
      res = $(Profile).@profile $(esc(exp))
      profiler(; $(map(esc, kwargs)...))
      res
    end
  end
end

"""
    structure(x)

Display `x`'s underlying representation, rather than using its normal display
method.

For example, `structure(:(2x+1))` displays the `Expr` object with its
`head` and `args` fields instead of printing the expression.
"""
structure(args...) = Main.Atom.structure(args...)


text:
module Isotropic

export bulk_modulus,
    young_modulus, lame1st, shear_modulus, poisson_ratio, longitudinal_modulus

const ALLOWED_KEYS = (:K, :E, :λ, :G, :ν, :M)

bulk_modulus(; kwargs...) = bulk_modulus(NamedTuple(kwargs))
bulk_modulus((E, λ)::NamedTuple{(:E, :λ)}) = (E + 3λ + _R(E, λ)) / 6
bulk_modulus((E, G)::NamedTuple{(:E, :G)}) = E * G / 3(3G - E)
bulk_modulus((E, ν)::NamedTuple{(:E, :ν)}) = E / 3(1 - 2ν)
bulk_modulus((E, M)::NamedTuple{(:E, :M)}) = (3M - E + _S(E, M)) / 6
bulk_modulus((λ, G)::NamedTuple{(:λ, :G)}) = λ + 2G / 3
bulk_modulus((λ, ν)::NamedTuple{(:λ, :ν)}) = λ * (1 + ν) / 3ν
bulk_modulus((λ, M)::NamedTuple{(:λ, :M)}) = (M + 2λ) / 3
bulk_modulus((G, ν)::NamedTuple{(:G, :ν)}) = 2G * (1 + ν) / 3(1 - 2ν)
bulk_modulus((G, M)::NamedTuple{(:G, :M)}) = M - 4G / 3
bulk_modulus((ν, M)::NamedTuple{(:ν, :M)}) = M * (1 + ν) / 3(1 - ν)
function bulk_modulus(x::NamedTuple)
    _checkkeys(x)
    return haskey(x, :K) ? x[:K] : bulk_modulus(_reverse(x))
end

young_modulus(; kwargs...) = young_modulus(NamedTuple(kwargs))
young_modulus((K, λ)::NamedTuple{(:K, :λ)}) = 9K * (K - λ) / (3K - λ)
young_modulus((K, G)::NamedTuple{(:K, :G)}) = 9K * G / (3K + G)
young_modulus((K, ν)::NamedTuple{(:K, :ν)}) = 3K * (1 - 2ν)
young_modulus((K, M)::NamedTuple{(:K, :M)}) = 9K * (M - K) / (3K + M)
young_modulus((λ, G)::NamedTuple{(:λ, :G)}) = G * (3λ + 2G) / (λ + G)
young_modulus((λ, ν)::NamedTuple{(:λ, :ν)}) = λ * (1 + ν) * (1 - 2ν) / ν
young_modulus((λ, M)::NamedTuple{(:λ, :M)}) = (M - λ) * (M + 2λ) / (M + λ)
young_modulus((G, ν)::NamedTuple{(:G, :ν)}) = 2G * (1 + ν)
young_modulus((G, M)::NamedTuple{(:G, :M)}) = G * (3M - 4G) / (M - G)
young_modulus((ν, M)::NamedTuple{(:ν, :M)}) = M * (1 + ν) * (1 - 2ν) / (1 - ν)
function young_modulus(x::NamedTuple)
    _checkkeys(x)
    return haskey(x, :E) ? x[:E] : young_modulus(_reverse(x))
end

lamé1st(; kwargs...) = lamé1st(NamedTuple(kwargs))
lamé1st((K, E)::NamedTuple{(:K, :E)}) = (9K^2 - 3K * E) / (9K - E)
lamé1st((K, G)::NamedTuple{(:K, :G)}) = K - 2G / 3
lamé1st((K, ν)::NamedTuple{(:K, :ν)}) = 3K * ν / (1 + ν)
lamé1st((K, M)::NamedTuple{(:K, :M)}) = (3K - M) / 2
lamé1st((E, G)::NamedTuple{(:E, :G)}) = G * (E - 2G) / (3G - E)
lamé1st((E, ν)::NamedTuple{(:E, :ν)}) = E * ν / (1 + ν) / (1 - 2ν)
lamé1st((E, M)::NamedTuple{(:E, :M)}) = (M - E + _S(E, M)) / 4
lamé1st((G, ν)::NamedTuple{(:G, :ν)}) = 2G * ν / (1 - 2ν)
lamé1st((G, M)::NamedTuple{(:G, :M)}) = M - 2G
lamé1st((ν, M)::NamedTuple{(:ν, :M)}) = M * ν / (1 - ν)
function lamé1st(x::NamedTuple)
    _checkkeys(x)
    return haskey(x, :λ) ? x[:λ] : lamé1st(_reverse(x))
end
const lame1st = lamé1st

shear_modulus(; kwargs...) = shear_modulus(NamedTuple(kwargs))
shear_modulus((K, E)::NamedTuple{(:K, :E)}) = 3K * E / (9K - E)
shear_modulus((K, λ)::NamedTuple{(:K, :λ)}) = 3(K - λ) / 2
shear_modulus((K, ν)::NamedTuple{(:K, :ν)}) = 3K * (1 - 2ν) / 2(1 + ν)
shear_modulus((K, M)::NamedTuple{(:K, :M)}) = 3(M - K) / 4
shear_modulus((E, λ)::NamedTuple{(:E, :λ)}) = (E - 3λ + _R(E, λ)) / 4
shear_modulus((E, ν)::NamedTuple{(:E, :ν)}) = E / 2(1 + ν)
shear_modulus((E, M)::NamedTuple{(:E, :M)}) = (3M + E - _S(E, M)) / 8
shear_modulus((λ, ν)::NamedTuple{(:λ, :ν)}) = λ * (1 - 2ν) / 2ν
shear_modulus((λ, M)::NamedTuple{(:λ, :M)}) = (M - λ) / 2
shear_modulus((ν, M)::NamedTuple{(:ν, :M)}) = M * (1 - 2ν) / 2(1 - ν)
function shear_modulus(x::NamedTuple)
    _checkkeys(x)
    return haskey(x, :G) ? x[:G] : shear_modulus(_reverse(x))
end

poisson_ratio(; kwargs...) = poisson_ratio(NamedTuple(kwargs))
poisson_ratio((K, E)::NamedTuple{(:K, :E)}) = (3K - E) / 6K
poisson_ratio((K, λ)::NamedTuple{(:K, :λ)}) = λ / (3K - λ)
poisson_ratio((K, G)::NamedTuple{(:K, :G)}) = (3K - 2G) / 2(3K + G)
poisson_ratio((K, M)::NamedTuple{(:K, :M)}) = (3K - M) / (3K + M)
poisson_ratio((E, λ)::NamedTuple{(:E, :λ)}) = 2λ / (E + λ + _R(E, λ))
poisson_ratio((E, G)::NamedTuple{(:E, :G)}) = E / 2G - 1
poisson_ratio((E, M)::NamedTuple{(:E, :M)}) = (E - M + _S(E, M)) / 4M
poisson_ratio((λ, G)::NamedTuple{(:λ, :G)}) = λ / 2(λ + G)
poisson_ratio((λ, M)::NamedTuple{(:λ, :M)}) = λ / (M + λ)
poisson_ratio((G, M)::NamedTuple{(:G, :M)}) = (M - 2G) / 2(M - G)
function poisson_ratio(x::NamedTuple)
    _checkkeys(x)
    return haskey(x, :ν) ? x[:ν] : poisson_ratio(_reverse(x))
end

longitudinal_modulus(; kwargs...) = longitudinal_modulus(NamedTuple(kwargs))
longitudinal_modulus((K, E)::NamedTuple{(:K, :E)}) = 3K * (3K + E) / (9K - E)
longitudinal_modulus((K, λ)::NamedTuple{(:K, :λ)}) = 3K - 2λ
longitudinal_modulus((K, G)::NamedTuple{(:K, :G)}) = K + 4G / 3
longitudinal_modulus((K, ν)::NamedTuple{(:K, :ν)}) = 3K * (1 - ν) / (1 + ν)
longitudinal_modulus((E, λ)::NamedTuple{(:E, :λ)}) = (E - λ + _R(E, λ)) / 2
longitudinal_modulus((E, G)::NamedTuple{(:E, :G)}) = G * (4G - E) / (3G - E)
longitudinal_modulus((E, ν)::NamedTuple{(:E, :ν)}) = E * (1 - ν) / (1 + ν) / (1 - 2ν)
longitudinal_modulus((λ, G)::NamedTuple{(:λ, :G)}) = λ + 2G
longitudinal_modulus((λ, ν)::NamedTuple{(:λ, :ν)}) = λ * (1 - ν) / ν
longitudinal_modulus((G, ν)::NamedTuple{(:G, :ν)}) = 2G * (1 - ν) / (1 - 2ν)
function longitudinal_modulus(x::NamedTuple)
    _checkkeys(x)
    return haskey(x, :M) ? x[:M] : longitudinal_modulus(_reverse(x))
end

# These are helper functions and should not be exported!
_R(E, λ) = sqrt(E^2 + 9λ^2 + 2E * λ)
_S(E, M) = sqrt(E^2 + 9M^2 - 10E * M)  # FIXME: ±S

_reverse(x::NamedTuple) = (; zip(reverse(propertynames(x)), reverse(values(x)))...)

_checkkeys(x) = @assert all(key ∈ ALLOWED_KEYS for key in propertynames(x))

end


text:
@syms x

function generate_rand_poly(x; min_deg=0, max_deg=20, sparcity=0.5, rational=0.25)
    n = rand(min_deg:max_deg)
    p = rand() < rational ? Poly(Rational, x) : Poly(x)
    r = rand(1:10) // rand(1:10)
    for i = 0:n
        if rand() > sparcity || i == n
            s = rand() < 0.5 ? +1 : -1
            p[i] = (p.T == Rational ? s * r * rand(0:10)//rand(1:5) : s * 10 * randn())
        end
    end
    p
end

function test_eq(x, f; n=10, abstol=1e-8, min_deg=0, max_deg=20, sparcity=0.5, rational=0.25)
    k = 0
    for i = 1:n
        p = generate_rand_poly(x; max_deg, sparcity, rational)
        q = generate_rand_poly(x; min_deg, max_deg, sparcity, rational)
        try
            r = f(p, q)
            if abs(maximum(extract_coef(r))) > abstol
                @warn "+/- mismatch: $r"
                println(p, ", ", q)
            else
                k += 1
            end
        catch e
            println(e)
            println(p, ", ", q)
        end
    end
    @info "$k ok!"
end

test_addsub(x) = test_eq(x, (p,q)->(p+q)-p-q)
test_muldiv(x) = test_eq(x, (p,q)->p-(p/q)*q-(p%q))
test_gcd(x) = test_eq(x, (p,q)->p%gcd_naive(p,q)+q%gcd_naive(p,q); max_deg=5, rational=1.0)

function test_deriv(x; n=10, min_deg=0, max_deg=20, sparcity=0.5, rational=1)
    k = 0
    for i = 1:n
        p = generate_rand_poly(x; min_deg, max_deg, sparcity, rational)
        try
            q = derivative(p)
            Δp = p - integrate(q)
            if degree(Δp) > 0
                @warn "deriv/integrate mismatch: $Δp"
            else
                k += 1
            end
        catch e
            println(e)
        end
    end
    @info "$k ok!"
end

function test_shubert(x; n=10, min_deg=1, max_deg=3, sparcity=0.5, rational=1)
    k = 0
    for i = 1:n
        p1 = generate_rand_poly(x; min_deg, max_deg, sparcity, rational)
        p2 = generate_rand_poly(x; min_deg, max_deg, sparcity, rational)
        p3 = generate_rand_poly(x; min_deg, max_deg, sparcity, rational)
        _, p1 = integer_poly(p1)
        _, p2 = integer_poly(p2)
        _, p3 = integer_poly(p3)
        p = p1 * p2 * p3
        try
            println('(',p1,") * (",p2,") * (",p3,')')
            println("\t=> ", factor_schubert_kronecker(p))
            k += 1
        catch e
            println(e)
        end
    end
    @info "$k ok!"
end

function test_factor(x)
    ps = [
        gen_poly(x,  64, 56, 14, 1),
        gen_poly(x,  -12, 53, -57, 18),
        gen_poly(x,  273, -86, -73, 6),
        gen_poly(x,  -2, -1, 4, -1, 6),
        gen_poly(x,  -30, -33, 8, -11, 6),
        gen_poly(x,  -2, 7, 20, -24, -6, 5),
        gen_poly(x,  8, -38, 27, 47, -11, 15),
        poly(x^4 - 4, x),
        poly(x^4 -8x^2 - 9, x),
        poly(6x^4 - 7x^3 + 5x^2 - 20x + 17, x),
        poly(x^6 - 1, x),
        poly(x^6 + 1, x),
        poly(x^5 + x + 1, x),
        gen_poly(x,  -35, 11, 6),
        gen_poly(x,  -16, 25),
        gen_poly(x,  4, -13, 24, -19, 6),
        gen_poly(x,  -2, 3, 5, 4, 2, 2, 2),
        gen_poly(x,  -27, 117, -90, 90, -48, 8),
        gen_poly(x,  2, 9, 25, 35, 39, 30),
        gen_poly(x,  -1, 1, 2, -2, -1, 1),
        gen_poly(x,  1, 2, 1, 1, 2, 1, 1, 2, 2, 2, 1),
        poly(x^8 - 4x^6 + 16x^2 - 16, x),
        gen_poly(x,  -1, 0, 3, 2, -3, -6, 6, 3, -2, -3, 1),
    ]

    k = 0
    for p = ps
        try
            printstyled(p, '\n'; color=:green)
            f = factor(p; verbose=true)
            println("\t", f)
            k += 1
            # q = expand(*(f...))
            # if iszero(p - q)
            #     @info "OK!"
            #     k += 1
            # else
            #     @warn "mismatch!"
            # end
        catch e
            println(e)
        end
    end
    @info "$k ok!"
end


text:
using GeometricIntegrators
using Weave

ENV["GKSwstype"] = "100"

tutorial_path = joinpath(dirname(pathof(GeometricIntegrators)), "../docs/tutorial", "tutorial.jmd")
build_path = joinpath(dirname(pathof(GeometricIntegrators)), "../docs/build/tutorial")

weave(tutorial_path,
  out_path=build_path,
  doctype = "md2html")


text:
using Juno
using Test

@test Juno.isactive() == false

@test Juno.plotsize() == [100, 100]

@test_throws Exception selector(["foo", "bar", "baz"])

let i = 0, x
  x = @progress for _ = 1:100
    i += 1
  end
  @test i == 100
  @test x == nothing
end

let i = 0, r = -50:10:50, x
  x = @progress for _ in r
    i += 1
  end
  @test i == 11
  @test x == nothing
end

let i = 0, x
  x = @progress "named" for _ = 1:100
    i += 1
  end
  @test i == 100
  @test x == nothing
end

let i = 0, j = 0, x
  x = @progress for _ = 1:10, __ = 1:20
    i += 1
  end
  @test i == 200
  @test x == nothing
end

let i = 0, j = 0, x
  bar = "bar"
  x = @progress "foo $bar" for _ = 1:10
    i += 1
  end
  @test i == 10
  @test x == nothing
end

let x,y
  x = @progress y = [i+3j for i=1:3, j=1:4]
  @test y == reshape(4:15,3,4)
  @test x == y
end

let a = [], x
  x = @progress for i=1:3, j=[-5,-2,-1,8]
      j > 0 && continue
      push!(a,(i,j))
      i > 1 && break
  end
  @test a == [(1,-5),(1,-2),(1,-1),(2,-5)]
  @test x == nothing
end

@test Juno.notify("hi") == nothing


text:
const turing_key_map = Dict(
    "acceptance_rate" => "mean_tree_accept",
    "hamiltonian_energy" => "energy",
    "hamiltonian_energy_error" => "energy_error",
    "is_adapt" => "tune",
    "max_hamiltonian_energy_error" => "max_energy_error",
    "n_steps" => "tree_size",
    "numerical_error" => "diverging",
    "tree_depth" => "depth",
)
const stan_key_map = Dict(
    "accept_stat__" => "accept_stat",
    "divergent__" => "diverging",
    "energy__" => "energy",
    "lp__" => "lp",
    "n_leapfrog__" => "n_leapfrog",
    "stepsize__" => "stepsize",
    "treedepth__" => "treedepth",
)
const stats_key_map = merge(turing_key_map, stan_key_map)

"""
    reshape_values(x::AbstractArray) -> AbstractArray

Convert from `MCMCChains` variable values with dimensions `(ndraw, size..., nchain)` to
ArviZ's expected `(nchain, ndraw, size...)`.
"""
reshape_values(x::AbstractArray{T,N}) where {T,N} = permutedims(x, (N, 1, 2:(N - 1)...))

headtail(x) = x[1], x[2:end]

function split_locname(name)
    name = replace(name, r"[\[,]" => '.')
    name = replace(name, ']' => "")
    name, loc = headtail(split(name, '.'))
    length(loc) == 0 && return name, ()
    loc = tryparse.(Int, loc)
    Nothing <: eltype(loc) && return name, ()
    return name, tuple(loc...)
end

function varnames_locs_dict(loc_names, loc_str_to_old)
    vars_to_locs = Dict()
    for loc_name in loc_names
        var_name, loc = split_locname(loc_name)
        if var_name ∉ keys(vars_to_locs)
            vars_to_locs[var_name] = ([loc_str_to_old[loc_name]], [loc])
        else
            push!(vars_to_locs[var_name][1], loc_str_to_old[loc_name])
            push!(vars_to_locs[var_name][2], loc)
        end
    end
    return vars_to_locs
end

function attributes_dict(chns::Chains)
    info = delete(chns.info, :hashedsummary)
    return Dict{String,Any}((string(k), v) for (k, v) in pairs(info))
end
attributes_dict(::Nothing) = Dict()

function section_dict(chns::Chains, section)
    ndraws, _, nchains = size(chns)
    loc_names_old = getfield(chns.name_map, section) # old may be Symbol or String
    loc_names = string.(loc_names_old)
    loc_str_to_old =
        Dict(name_str => name_old for (name_str, name_old) in zip(loc_names, loc_names_old))
    vars_to_locs = varnames_locs_dict(loc_names, loc_str_to_old)
    vars_to_arrays = Dict{String,Array}()
    for (var_name, names_locs) in vars_to_locs
        loc_names, locs = names_locs
        max_loc = maximum(hcat([[loc...] for loc in locs]...); dims = 2)
        ndim = length(max_loc)
        sizes = tuple(max_loc...)

        oldarr = reshape_values(replacemissing(Array(chns.value[:, loc_names, :])))
        if ndim == 0
            arr = dropdims(oldarr; dims = 3)
        else
            arr = Array{Union{typeof(NaN),eltype(oldarr)}}(undef, nchains, ndraws, sizes...)
            fill!(arr, NaN)
            for i in eachindex(locs)
                arr[:, :, locs[i]...] = oldarr[:, :, i]
            end
        end
        vars_to_arrays[var_name] = arr
    end
    return vars_to_arrays
end

function chains_to_dict(
    chns::Chains;
    ignore = String[],
    section = :parameters,
    rekey_fun = identity,
)
    section in sections(chns) || return Dict()
    chns_dict = section_dict(chns, section)
    removekeys!(chns_dict, ignore)
    return rekey_fun(chns_dict)
end
chains_to_dict(::Nothing; kwargs...) = nothing

"""
    convert_to_inference_data(obj::Chains; group = :posterior, kwargs...) -> InferenceData

Convert the chains `obj` to an [`InferenceData`](@ref) with the specified `group`.

Remaining `kwargs` are forwarded to [`from_mcmcchains`](@ref).
"""
function convert_to_inference_data(chns::Chains; group = :posterior, kwargs...)
    group = Symbol(group)
    group == :posterior && return from_mcmcchains(chns; kwargs...)
    return from_mcmcchains(; group => chns)
end

@doc doc"""
    from_mcmcchains(posterior::Chains; kwargs...) -> InferenceData
    from_mcmcchains(; kwargs...) -> InferenceData
    from_mcmcchains(
        posterior::Chains,
        posterior_predictive::Any,
        predictions::Any,
        log_likelihood::Any;
        kwargs...
    ) -> InferenceData

Convert data in an `MCMCChains.Chains` format into an [`InferenceData`](@ref).

Any keyword argument below without an an explicitly annotated type above is allowed, so long
as it can be passed to [`convert_to_inference_data`](@ref).

# Arguments

- `posterior::Chains`: Draws from the posterior

# Keywords

- `posterior_predictive::Any=nothing`: Draws from the posterior predictive distribution or
    name(s) of predictive variables in `posterior`
- `predictions::Any=nothing`: Out-of-sample predictions for the posterior.
- `prior::Any=nothing`: Draws from the prior
- `prior_predictive::Any=nothing`: Draws from the prior predictive distribution or name(s)
    of predictive variables in `prior`
- `observed_data::Dict{String,Array}=nothing`: Observed data on which the `posterior` is
    conditional. It should only contain data which is modeled as a random variable. Keys are
    parameter names and values.
- `constant_data::Dict{String,Array}=nothing`: Model constants, data included in the model
    which is not modeled as a random variable. Keys are parameter names and values.
- `predictions_constant_data::Dict{String,Array}=nothing`: Constants relevant to the model
     predictions (i.e. new `x` values in a linear regression).
- `log_likelihood::Any=nothing`: Pointwise log-likelihood for the data. It is recommended
     to use this argument as a dictionary whose keys are observed variable names and whose
     values are log likelihood arrays.
- `log_likelihood::String=nothing`: Name of variable in `posterior` with log likelihoods
- `library=MCMCChains`: Name of library that generated the chains
- `coords::Dict{String,Vector}=nothing`: Map from named dimension to named indices
- `dims::Dict{String,Vector{String}}=nothing`: Map from variable name to names of its
    dimensions

# Returns

- `InferenceData`: The data with groups corresponding to the provided data
"""
from_mcmcchains

function from_mcmcchains(
    posterior,
    posterior_predictive,
    predictions,
    log_likelihood;
    library = MCMCChains,
    kwargs...,
)
    kwargs = convert(Dict, merge((; dims = nothing), kwargs))
    library = string(library)
    rekey_fun = d -> rekey(d, stats_key_map)

    # Convert chains to dicts
    post_dict = chains_to_dict(posterior)
    stats_dict = chains_to_dict(posterior; section = :internals, rekey_fun = rekey_fun)
    stats_dict = enforce_stat_types(stats_dict)

    all_idata = InferenceData()
    for (group, group_data) in [
        :posterior_predictive => posterior_predictive,
        :predictions => predictions,
        :log_likelihood => log_likelihood,
    ]
        group_data === nothing && continue
        if group_data isa Union{Symbol,String}
            group_data = [string(group_data)]
        end
        if group_data isa Union{AbstractVector{Symbol},NTuple{N,Symbol} where {N}}
            group_data = map(string, group_data)
        end
        if group_data isa Union{AbstractVector{String},NTuple{N,String} where {N}}
            group_data = popsubdict!(post_dict, group_data)
        end
        group_dataset = convert_to_dataset(group_data; library = library, kwargs...)
        setattribute!(group_dataset, "inference_library", library)
        concat!(all_idata, InferenceData(; group => group_dataset))
    end

    attrs = attributes_dict(posterior)
    attrs = merge(attrs, Dict("inference_library" => library))
    kwargs = convert(Dict, merge((; attrs = attrs, dims = nothing), kwargs))
    post_idata = _from_dict(post_dict; sample_stats = stats_dict, kwargs...)
    concat!(all_idata, post_idata)
    return all_idata
end
function from_mcmcchains(
    posterior = nothing;
    posterior_predictive = nothing,
    predictions = nothing,
    prior = nothing,
    prior_predictive = nothing,
    observed_data = nothing,
    constant_data = nothing,
    predictions_constant_data = nothing,
    log_likelihood = nothing,
    library = MCMCChains,
    kwargs...,
)
    kwargs = convert(Dict, merge((; dims = nothing, coords = nothing), kwargs))

    all_idata = from_mcmcchains(
        posterior,
        posterior_predictive,
        predictions,
        log_likelihood;
        library = library,
        kwargs...,
    )

    if prior !== nothing
        pre_prior_idata = convert_to_inference_data(
            prior;
            posterior_predictive = prior_predictive,
            library = library,
            kwargs...,
        )
        prior_idata = rekey(
            pre_prior_idata,
            Dict(
                :posterior => :prior,
                :posterior_predictive => :prior_predictive,
                :sample_stats => :sample_stats_prior,
            ),
        )
        concat!(all_idata, prior_idata)
    end

    for (group, group_data) in [
        :observed_data => observed_data,
        :constant_data => constant_data,
        :predictions_constant_data => predictions_constant_data,
    ]
        group_data === nothing && continue
        group_dict = convert(Dict, group_data)
        group_dataset =
            convert_to_constant_dataset(group_dict; library = library, kwargs...)
        concat!(all_idata, InferenceData(; group => group_dataset))
    end

    return all_idata
end

"""
    from_cmdstan(posterior::Chains; kwargs...) -> InferenceData

Call [`from_mcmcchains`](@ref) on output of `CmdStan`.
"""
function from_cmdstan(posterior::Chains; kwargs...)
    return from_mcmcchains(posterior; library = "CmdStan", kwargs...)
end


text:
abstract type AbstractExpression end

include("variable.jl")
include("expressions/expression.jl")
include("terms/term.jl")
include("problem.jl")


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export proj
"""
`proj(get_iolabel::Function, fst::WFST)`

Projects input labels to output labels (or viceversa).
The function `get_iolabel` should either be the function `get_ilabel` or `get_olabel`.

```julia
julia> A = linearfst(["a","b","c"],[1,2,3],ones(3))
WFST #states: 4, #arcs: 3, #isym: 3, #osym: 3
|1/1.0|
a:1/1.0 → (2)
(2)
b:2/1.0 → (3)
(3)
c:3/1.0 → (4)
((4/1.0))

julia> proj(get_ilabel, A)
WFST #states: 4, #arcs: 3, #isym: 3, #osym: 3
|1/1.0|
a:a/1.0 → (2)
(2)
b:b/1.0 → (3)
(3)
c:c/1.0 → (4)
((4/1.0))

julia> proj(get_olabel, A)
WFST #states: 4, #arcs: 3, #isym: 3, #osym: 3
|1/1.0|
1:1/1.0 → (2)
(2)
2:2/1.0 → (3)
(3)
3:3/1.0 → (4)
((4/1.0))

```
"""
function proj(get_iolabel::Function, fst::WFST)
  if !( (get_iolabel == get_ilabel) || (get_iolabel == get_olabel) )
    throw(ErrorException("Input function should be either `ilabel` or `olabel`."))
  end
  pfst = arcmap(proj, fst, get_iolabel)
  if get_iolabel == get_ilabel
    pfst = WFST(pfst.states,pfst.initials,pfst.finals,pfst.isym,pfst.isym)
  else
    pfst = WFST(pfst.states,pfst.initials,pfst.finals,pfst.osym,pfst.osym)
  end
  return pfst
end


text:
# ---
# title: 1632. Rank Transform of a Matrix
# id: problem1632
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Greedy, Union Find
# link: <https://leetcode.com/problems/rank-transform-of-a-matrix/description/>
# hidden: true
# ---
#
# Given an `m x n` `matrix`, return _a new matrix_`answer`
# _where_`answer[row][col]` _is the_ _ **rank** of _`matrix[row][col]`.
#
# The **rank** is an **integer** that represents how large an element is
# compared to other elements. It is calculated using the following rules:
#
#   * The rank is an integer starting from `1`.
#   * If two elements `p` and `q` are in the **same row or column** , then:
#     * If `p < q` then `rank(p) < rank(q)`
#     * If `p == q` then `rank(p) == rank(q)`
#     * If `p > q` then `rank(p) > rank(q)`
#   * The **rank** should be as **small** as possible.
#
# It is guaranteed that `answer` is unique under the given rules.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/18/rank1.jpg)
#
#
#
#     Input: matrix = [[1,2],[3,4]]
#     Output: [[1,2],[2,3]]
#     Explanation:
#     The rank of matrix[0][0] is 1 because it is the smallest integer in its row and column.
#     The rank of matrix[0][1] is 2 because matrix[0][1] > matrix[0][0] and matrix[0][0] is rank 1.
#     The rank of matrix[1][0] is 2 because matrix[1][0] > matrix[0][0] and matrix[0][0] is rank 1.
#     The rank of matrix[1][1] is 3 because matrix[1][1] > matrix[0][1], matrix[1][1] > matrix[1][0], and both matrix[0][1] and matrix[1][0] are rank 2.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/18/rank2.jpg)
#
#
#
#     Input: matrix = [[7,7],[7,7]]
#     Output: [[1,1],[1,1]]
#
#
# **Example 3:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/18/rank3.jpg)
#
#
#
#     Input: matrix = [[20,-21,14],[-19,4,19],[22,-47,24],[-19,4,19]]
#     Output: [[4,2,3],[1,3,4],[5,1,6],[1,3,4]]
#
#
# **Example 4:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/18/rank4.jpg)
#
#
#
#     Input: matrix = [[7,3,6],[1,4,5],[9,8,2]]
#     Output: [[5,1,4],[1,2,3],[6,3,1]]
#
#
#
#
# **Constraints:**
#
#   * `m == matrix.length`
#   * `n == matrix[i].length`
#   * `1 <= m, n <= 500`
#   * `-109 <= matrix[row][col] <= 109`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Optim

#cumulpdf(x,alpha,xmax) =  x^(1-alpha) - xmax^(1-alpha)

cumulpdf(x,a,xmax) = (x^(1-a)-xmax^(1-a)) / (1-xmax^(1-a))

function loghistogram(ts, nb=20)
# histogram of log-transformed values
   his = fit(Histogram, ts, nbins=nb)
   x1 = his.edges[1][2:end]
   ind = findall(x->x != 0, his.weights)  # remove bins with entry 0
   w = his.weights[ind]
   w = w ./ float(sum(w))
   x1 = x1[ind]
   y1 = log.(w) .- x1
   ex1 = exp.(x1)
   ey1 = exp.(y1)
   x1,y1,ex1,ey1
end

function log2histogram1(ts)
# histogram with logarithmic binning
# bin positions at powers of two, 2,4,8,16,...
    maxt = maximum(ts)
    maxind = ceil(Int,log2(maxt))
    bins = 2.0 .^ (1:maxind)   # bins
    nbins = length(bins)
    w = zeros(nbins)        # weights
    for i=1:length(ts)
       t = ts[i]
       if t == 1
           w[1] += 1
       else
           ind = ceil(Int,log2(t))
           w[ind] +=1
       end
    end
    ind = findall(x->x != 0, w)  # remove bins with entry 0
    w = w[ind]
    bins = bins[ind]
    distr = w ./ bins  # get distribution by normalization
    distr = distr / sum(distr)

    bins, distr
end

function log2histogram(ts)
# histogram with logarithmic binning
# bin positions at powers of two, 1,2,4,8,16,...
    maxt = maximum(ts)
    maxind = ceil(Int,log2(maxt))
    bins = 2.0 .^ (0:maxind)   # bins
    nbins = length(bins)
    weight = zeros(nbins)        # weights
    widths = zeros(nbins)        # widths of bins
    for t in ts
       ind = ceil(Int,log2(t)) + 1
       weight[ind] +=1
    end
    widths[1] = 0.5
    widths[2:end] = bins[1:end-1]
    ind = findall(x->x != 0, weight)  # remove bins with entry 0
    weight = weight[ind]
    widths = widths[ind]
    bins   = bins[ind]
    distr = weight ./ widths  # get distribution by normalization
    distr = distr / sum(distr)

    bins, distr
end

# log likelihood
like(mu,g,nmax) = log((1-mu)/(nmax^(1-mu)-1)) - mu*g

like_inf(mu,g) = log(mu-1) - mu*g

function standard_err(M,mu,nmax)
# standard error of exponent
   mu1 = mu-1.0
   denom = 1-nmax^mu1
   1.0 / sqrt(1.0/(mu1*mu1) - nmax^mu1 * log(log(nmax))/(denom*denom)) /
   sqrt(M)
end

function KS_test(C,n,mu,nmax)
    # Kolmogorov-Smirnov test
   diff = abs.(cumulpdf.(C,mu,nmax) .- (1:n)/n)
   maximum(diff)
end

function estimate(C)
   nmax = maximum(C)
   n = length(C)
   estimate(C,n,nmax)
end

function estimate(C,n,nmax)
   # estimate fitting values for truncated power law
   # critical exponent, standard deviation Kolmogorov-Smirnov distance
   g = sum(log.(C)) / n   # log(geometric mean)
   #g = sum(log.(2*C)) / n   # try correction from Clauset
   # find maximum using Brent's method
   mu = optimize(x->-like(x,g,nmax), 1.01, 3.0).minimizer
  # mu = optimize(x->-like(x,g,nmax), 1.01, 3.0, GoldenSection()).minimizer
   st = standard_err(n,mu,nmax)
   #d = KS_test(C,n,mu,nmax)
   #mu, st, d
   mu, st
end

function randPL(n,mu,nmax)
# generate sample of random numbers taken from a truncated power law
   mu1 = 1.0 - mu
   r = (1.0 .- rand(n) .* (1.0 - nmax^mu1)) .^ (1.0/mu1)
   sort(r, rev=true)
end


text:
## Landcover variables loading functions

# Landcover variables loading
function landcover(layers::Vector{Int64}; resolution::AbstractString="10", path::AbstractString="assets")
    ## Get file paths
    # List files in path
    lc_files = readdir("$(path)")
    # Filter for selected resolution
    filter!(x -> occursin.("lc_", x), lc_files)
    # Create paths for selected layers only
    paths = joinpath.(path, lc_files)[layers]

    ## Load data
    # Read raster dataset
    datasets = [ArchGDAL.read(p) for p in paths]
    # Read raster values from band 1 (only band in this case)
    values = [ArchGDAL.read(d, 1) for d in datasets]
    # Convert from UInt8 to Float64
    values_int = [convert.(Float64, v) for v in values]

    ## Preliminary manipulations
    # Reverse rows and permute dimensions
    landcover_mat = [permutedims(v[:,end:-1:1]) for v in values_int]
    # Replace 255 (default no data values) by NaN
    [replace!(l, 255 => NaN) for l in landcover_mat]

    # Fill missing latitudes with NaNs (latitude extent is only (-60,80) for landcover data)
    nlat, nlon = size(landcover_mat[1])
    slim, nlim = abs(-60), 80
    res = Int64(nlat/(nlim+slim))
    south_nans = fill(NaN, ((90-slim)*res, nlon))
    north_nans = fill(NaN, ((90-nlim)*res, nlon))
    landcover_grids = [vcat(south_nans, l, north_nans) for l in landcover_mat]

    # Convert to SimpleSDMLayers
    landcover_layers = SimpleSDMPredictor.(landcover_grids, -180.0, 180.0, -90.0, 90.0)

    return landcover_layers
end

"""
    landcover(layer::Int64; x...)

Return a single landcover layer.
"""
landcover(layer::Int64; x...) = landcover([layer]; x...)[1]

"""
    landcover(layers::UnitRange{Int64}; x...)

Return a range of landcover layers.
"""
landcover(layers::UnitRange{Int64}; x...) = landcover(collect(layers); x...)


text:
push!(LOAD_PATH,"../src/")

using Documenter, BayesMultiFit

makedocs(sitename="BayesMultiFit.jl Documentation",
modules=[BayesMultiFit]
)

deploydocs(;
    repo="github.com/LidkeLab/BayesMultiFit.jl.git"
)


text:
struct FreeTimeConstraints <: Constraints
	n::Int
	ineq
end

function free_time_constraints(T)
	n = T - 2
	ineq = (1:0)

	return FreeTimeConstraints(n, ineq)
end

function constraints!(c, Z, con::FreeTimeConstraints, model, idx, h, T)
	for t = 1:T-2
		h = view(Z,idx.u[t])[end]
		h⁺ = view(Z,idx.u[t + 1])[end]
		c[t] = h⁺ - h
	end

    return nothing
end

function constraints_jacobian!(∇c, Z, con::FreeTimeConstraints, model, idx, h, T)
	shift = 0
	for t = 1:T-2
		r_idx = (t:t)

		c_idx = (idx.u[t][end]:idx.u[t][end])
		len = 1
		∇c[shift + len] = -1.0
		shift += 1

		c_idx = (idx.u[t + 1][end]:idx.u[t + 1][end])
		len = 1
		∇c[shift + len] = 1.0
		shift += 1
	end

	return nothing
end

function constraints_sparsity(con::FreeTimeConstraints, model, idx, T;
	shift_row = 0, shift_col = 0)
    row = []
    col = []

	for t = 1:T-2
		r_idx = shift_row .+ (t:t)

		c_idx = shift_col .+ (idx.u[t][end]:idx.u[t][end])
		row_col!(row, col, r_idx, c_idx)

		c_idx = shift_col .+ (idx.u[t + 1][end]:idx.u[t + 1][end])
		row_col!(row, col, r_idx, c_idx)
	end
    return collect(zip(row, col))
end

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
using RCall

function snap_pour_point(flowaccu::GeoArray, shp, outfile = "outlet_snaped.shp"; ngrid = 15)
    @time res = clip_point(flowaccu, shp; ngrid = ngrid);
    lst = res["data"];
    shp2 = shp[:, setdiff(1:15, 13)]

    R"""
    library(data.table)
    library(purrr)
    # library()
    shp = $shp2
    # res = map($lst, ~data.table(.x)[value > 100]) %>% set_names(seq_along(.))
    res = $lst %>% set_names(seq_along(.))
    save(res, shp, file = "flowaccu.rda")

    df = snap_pour_point(res, shp)
    # "I:/ChinaBasins/shp/chinaRunoff_stationInfo_Yangtze_sp87_snaped.shp"
    sp2::df2sp(df) %>% sp2::write_shp($outfile)
    # $lst[[1]] %>% data.table()
    """
end

export snap_pour_point


text:
module LightGraphsFlows

using LightGraphs
const lg = LightGraphs

using SimpleTraits: @traitfn, @traitimpl
import SimpleTraits

using JuMP

using SparseArrays: spzeros, sparse, sparsevec
using Markdown: @doc_str

import Base: getindex, size, transpose, adjoint

include("maximum_flow.jl")
include("edmonds_karp.jl")
include("dinic.jl")
include("boykov_kolmogorov.jl")
include("push_relabel.jl")
include("multiroute_flow.jl")
include("kishimoto.jl")
include("ext_multiroute_flow.jl")
include("mincost.jl")
include("mincut.jl")

export
maximum_flow, EdmondsKarpAlgorithm, DinicAlgorithm, BoykovKolmogorovAlgorithm, PushRelabelAlgorithm,
multiroute_flow, KishimotoAlgorithm, ExtendedMultirouteFlowAlgorithm, mincost_flow

end


text:
# We mark `isless()` as `@pure` because it's used in computing the correct
# (sorted) monomial order during type promotion, so its results must be
# assumed to be hyper-pure. Overloading this method in user code would be
# a bad idea.
@pure function isless(::Type{V1}, ::Type{V2}) where {V1 <: Variable, V2 <: Variable}
    name(V1) < name(V2)
end

function promote_rule(::Type{V1}, ::Type{V2}) where {V1 <: Variable, V2 <: Variable}
    if V1 < V2
        Monomial{(V1(), V2()), 2}
    else
        Monomial{(V2(), V1()), 2}
    end
end

function promote_rule(::Type{M}, ::Type{V}) where {V <: Variable, M <: Monomial}
    promote_rule(Monomial{(V(),), 1}, M)
end
promote_rule(V::Type{<:Variable}, M::Type{<:Monomial}) = promote_rule(M, V)

function promote_rule(::Type{Term{T, M}}, ::Type{V}) where {V <: Variable, T, M <: Monomial}
    Term{T, promote_type(V, M)}
end
promote_rule(V::Type{<:Variable}, T::Type{<:Term}) = promote_rule(T, V)

function _promote_monomial_noncommutative(::Type{Monomial{V1, N1}}, ::Type{Monomial{V2, N2}}) where {N1, V1, N2, V2}
    if V1 > V2
        _promote_monomial(Monomial{V2, N2}, Monomial{V1, N1})
    else
        varnames = shortest_common_supersequence(name.(V1), name.(V2))
        vars = Expr(:tuple, [Expr(:call, Expr(:curly, :Variable, Expr(:quote, n))) for n in varnames]...)
        quote
            Monomial{$vars, $(length(varnames))}
        end
    end
end

@pure merge(v1::Tuple{Vararg{Variable}}) = v1

@pure function merge(v1::Tuple{Vararg{Variable}}, v2::Variable)
    for i in 1:length(v1)
        if v2 === v1[i]
            return v1
        elseif v2 > v1[i]
            return tuple(v1[1:i-1]..., v2, v1[i:end]...)
        end
    end
    return tuple(v1..., v2)
end

@pure function merge(v1::Tuple{Vararg{Variable}}, v2::Variable, vars::Variable...)
    merge(merge(v1, v2), vars...)
end

@pure function promote_rule(::Type{<:Monomial{V1}}, ::Type{<:Monomial{V2}}) where {V1, V2}
    vars = merge(V1, V2...)
    Monomial{vars, length(vars)}
end

function promote_rule(::Type{Term{T, M2}}, ::Type{M1}) where {M1 <: Monomial, T, M2 <: Monomial}
    Term{T, promote_type(M1, M2)}
end
promote_rule(M::Type{<:Monomial}, T::Type{<:Term}) = promote_rule(T, M)

promote_rule(::Type{Term{T1, M1}}, ::Type{Term{T2, M2}}) where {T1, M1 <: Monomial, T2, M2 <: Monomial} = Term{promote_type(T1, T2), promote_type(M1, M2)}

function promote_rule(::Type{<:Polynomial{<:Any, T2}}, ::Type{T1}) where {T1 <: TermLike, T2 <: Term}
    T = promote_type(T1, T2)
    Polynomial{coefficienttype(T), T, Vector{T}}
end
promote_rule(T::Type{Term{T1, M1}}, P::Type{<:Polynomial}) where {T1, M1 <: Monomial} = promote_rule(P, T)
promote_rule(T::Type{<:Monomial}, P::Type{<:Polynomial}) = promote_rule(P, T)
promote_rule(T::Type{<:Variable}, P::Type{<:Polynomial}) = promote_rule(P, T)

function promote_rule(::Type{<:Polynomial{<:Any, T1}}, ::Type{<:Polynomial{<:Any, T2}}) where {T1 <: TermLike, T2 <: TermLike}
    T = promote_type(T1, T2)
    Polynomial{coefficienttype(T), T, Vector{T}}
end

function MP.promote_rule_constant(::Type{S}, ::Type{V}) where {S, V <: Variable}
    Term{S, Monomial{(V(),), 1}}
end

function MP.promote_rule_constant(::Type{S}, ::Type{M}) where {S, M <: Monomial}
    Term{S, M}
end

MP.promote_rule_constant(::Type{S}, t::Type{Term{T, M}}) where {T, M <: Monomial, S} = Term{promote_type(S, T), M}

function MP.promote_rule_constant(::Type{S}, ::Type{<:Polynomial{<:Any, T}}) where {S, T <: TermLike}
    R = promote_type(S, T)
    Polynomial{coefficienttype(R), R, Vector{R}}
end


text:
function eigen(mpo::MPO)
    F = eigen(mpo2mat(mpo));
    return F.values
end

max_stars_repo_name: clarazen/BigMat

max_issues_repo_name: clarazen/BigMat

max_forks_repo_name: clarazen/BigMat


text:
using AlternatingProjections, FFTW
# using LinearAlgebra
using Test

@testset "AlternatingProjections.jl" begin
    # Write your own tests here.

    S = ConstrainedBySupport([true, false,true])
    A = ConstrainedByAmplitude([1, sqrt(2), 5])
    x = [1, 2, 3]
    y = [2im, -2. + 2im, 6 - 8im]
    @test project(x, S) == [1, 0, 3]
    @test project(y, A) ≈  [im, -1 + im, 3 - 4im]

    xr=copy(x)
    reflect!(xr,x,S)
    @test xr == [1, -2, 3]

    y = zeros(ComplexF64,10,10)  # for ComplexF64 increas the number of iterations
    y[1:5,1:5] = randn(ComplexF32, 5,5)
    Y = fft(y)
    # z = apsolve(abs.(y),abs.(Y), GS, maxit =3000,maxϵ = 1e-18)
    # @test abs.(z) ≈ abs.(y)
    # @test abs.(fft(z)) ≈  abs.(Y)

    # p= PR(abs.(y),abs.(Y))
    # gs=AP(3000,1e-6)
    # z= solve(p, one(y),gs)
    # @test abs.(z) ≈ abs.(y)
    # z, h, snap = solve(p, one(y),gs, true, [1, 5, 10])
    # @test length(h) == 3000
    # @test size(snap) == (10,10,3)

end

@testset "FourierTransformedSet" begin
    a = ConstrainedByAmplitude(ones(5,5))
    A = FourierTransformedSet(a);
    x = rand(Complex{Float64}, 5,5);
    x = project!(x,a)
    @test abs.(x) ≈ ones(5,5)

    y = similar(x);
    project!(y,x,A)
    @test abs.(fft(y)) ≈ ones(5,5)

    @test AlternatingProjections.getelement(a) == ones(5,5)
    B = zeros(Complex,(5,5)); B[1] +=1;
    @test AlternatingProjections.getelement(A) ≈ B

end

@testset "Amplitude and Shape Constraints" begin
    a = [ 8 - i^2 - j^2 for i in -2:2, j in -2:2]
    m = [ i >= 0 && j >= 0 for i in -2:2, j in -2:2]

    A = ConstrainedByAmplitude(a)
    B = ConstrainedByAmplitudeMasked(a,m)
    C = ConstrainedByShape(a)
    D = ConstrainedByShapeMasked(a, m)

    x = randn(ComplexF64, 5,5)
    y = project(x,A)
    z = project(x, B)
    v = project(x,C)
    w = project(x, D)

    @test abs.(y) ≈ a
    @test abs.(z)[m] ≈ a[m]
    @test abs.(v)/sum(abs.(v)) ≈ a /sum(a)

    s = abs.(x)[D.mask]' * D.amp[D.mask] / D.n

    @test abs.(w[m])/s ≈ a[m]

end

@testset "ConstrainedByShapeSaturated" begin
    a = map(x -> 1-x^2, -1:0.1:1)
    meas = 2 * a
    measset = ConstrainedByShapeSaturated(meas, meas .<=1)
    meassat = copy(meas)
    meassat[meassat .> 1] .= 1

    p = project(a, measset)
    @test p == a

    p2 = project(meas, measset)
    @test p2 == meas

end

@testset "solve" begin

    # Quick Gerchberg-Saxton
    n = 50
    m = 10

    amp = 0.35
    x = zeros(ComplexF64, n, n); x[1:m,1:m] .= exp.(amp * 2π * im * randn(Float64, m, m))
    a = abs.(x)
    X = fft(x)
    A = abs.(X)
    aset = ConstrainedByAmplitude(a)
    Aset = FourierTransformedSet(ConstrainedByAmplitude(A))
    problem = TwoSetsFP(aset,Aset)
    # start with a good initial guess
    sol = solve(problem, APparam(), x⁰ = x + 1.5 * randn(ComplexF64, n, n), maxit = 1000)

    function testsolution()
        phasediff = -pi .+ mod2pi.(pi .+ angle.(sol[1][1:m,1:m]) .- angle.(x[1:m,1:m]))
        phasediff .-= sum(phasediff)/m^2
        phaserms = sqrt(sum(abs2,phasediff))/m
        #phase difference can be visualised
        # heatmap(phasediff)
        println("phase rms = $phaserms")
        @test phaserms < 1e-6
    end

    testsolution()

    # and the same problem solved with DR method (it requires more iterations and other starting point)
    sol = solve(problem, DRparam(), x⁰ = x - 1.5 * randn(ComplexF64, n, n), maxit = 3000)

    testsolution()

    # and with DRAP algoritm
    drap = DRAPparam(missing,500,missing,true,[1],0.1)
    # sol = solve(problem, drap, maxit=1000);
    sol = solve(problem, drap, x⁰ = x + 1.5 * randn(ComplexF64, n, n), maxit=1000);
    testsolution()
end


text:
module EikonalInv
export EikonalInvParam
export getEikonalInvParam

using jInv.Mesh
using jInv.Utils
using FactoredEikonalFastMarching
using MAT

import jInv.ForwardShare.getData
import jInv.ForwardShare.getSensTMatVec
import jInv.ForwardShare.getSensMatVec

import jInv.ForwardShare.ForwardProbType

function getFieldsFileName()
	tfilename = string("tempEikFields_worker",myid(),".mat");
end

"""
type EikonalInvParam

A type to be used with jInv for solving the inverse eikonal equation (e.g., travel time tomography).

Fields:

	Mesh 			:: RegularMesh
	Sources 		:: SparseMatrixCSC  		  - The matrix of point source locations of size N x num_sources
	Receivers 		:: SparseMatrixCSC  		- The matrix of point reciever locations of size N x num_recievers
	HO          	::Bool              		- flag for higher-order discretization
	eikonalParams	::Array{EikonalParam}   - An array of forward problem parameters. See FactoredEikonalFastMarching for EikonalParam.
  useFilesForFields :: Bool             - flag for writing fields to disk

Constructor:

	getEikonalInvParam
"""
type EikonalInvParam <: ForwardProbType
	Mesh      		  :: RegularMesh
  Sources			  :: SparseMatrixCSC
  Receivers		  :: SparseMatrixCSC
	HO        		  :: Bool
	eikonalParams	  :: Array{EikonalParam}
	useFilesForFields :: Bool
end

import jInv.ForwardShare.getSensMatSize
getSensMatSize(pFor::EikonalInvParam) = (size(pFor.Receivers,2)*size(pFor.Sources,2),prod(pFor.Mesh.n+1))

"""
function EikonalInv.getEikonalInvParam

Constructor for EikonalInvParam

Input:

	Mesh::RegularMesh
	Sources::SparseMatrixCSC      - The matrix of point source locations of size N x num_sources
	Receivers::SparseMatrixCSC    - The matrix of point reciever locations of size N x num_recievers
	HO::Bool=false							  - flag for higher-order discretization
	useFilesForFields::Bool=false - flag for writing fields to disk
"""
function getEikonalInvParam(Mesh::RegularMesh,Sources::SparseMatrixCSC,Receivers::SparseMatrixCSC,HO::Bool=false,useFilesForFields::Bool = false)
	## This function does not use the parallel mechanism of jInv.
	return EikonalInvParam(Mesh,Sources,Receivers,HO,Array{EikonalParam}(0),useFilesForFields);
end

function getEikonalInvParam(Mesh::RegularMesh,Sources::SparseMatrixCSC,Receivers::SparseMatrixCSC,HO::Bool,numWorkers::Int64,useFilesForFields::Bool = false)
	## This function does use the parallel mechanism of jInv (i.e., returns a RemoteChannel), even if numWorkers=1.
	if numWorkers > nworkers()
		numWorkers = nworkers();
	end
	SourcesSubInd = Array{Array{Int64,1}}(numWorkers);
	ActualWorkers = workers();
	if numWorkers < nworkers()
		ActualWorkers = ActualWorkers[1:numWorkers];
	end
	continuationDivision = [1;numWorkers+1];

	pFor   = Array{RemoteChannel}(numWorkers)
	i = 1; nextidx() = (idx=i; i+=1; idx)

	nsrc  = size(Sources,2);
	# send out jobs
	@sync begin
		for p = ActualWorkers
			@async begin
				while true
					idx = nextidx()
					if idx > numWorkers
						break
					end
					I_p = getSourcesIndicesOfKthWorker(numWorkers,idx,nsrc);
					SourcesSubInd[idx] = I_p;
					# find src and rec on mesh
					pFor[idx]  = initRemoteChannel(getEikonalInvParam,p, Mesh, Sources[:,I_p], Receivers, HO,useFilesForFields);
					wait(pFor[idx])
				end
			end
		end
	end
	return pFor,continuationDivision,SourcesSubInd# Array of Remote Refs
end

import jInv.Utils.clear!
function clear!(pFor::EikonalInvParam)
for k = 1:length(pFor.eikonalParams)
	clear!(pFor.eikonalParams[k]);
end
return pFor;
end

include("./SeismicUtils/SeismicUtils.jl")
include("getData.jl")
include("getSensMatVec.jl")
include("getSensTMatVec.jl")

end


text:
using .StanSample

include("stan_plot_model_coef.jl")
include("stan_plot_models.jl")
include("stan_trankplot.jl")


text:
function foo(x)
  str = x("world")
  println("hello $(str)!")
end
foo(y -> "blue $y") # prints "hello blue world"

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode