text:
module CBSOData4

import HTTP
import JSON
import Tables

export get_catalogs, get_tables, get_meta, ODataTable_long

BASE_URL = "http://beta-odata4.cbs.nl"
CATALOG = "CBS"

function get_odata(url)
    result = HTTP.get(url)
    data = JSON.parse(String(result.body))
end

"""
    get_catalogs(base)

Get all the catalogs as a list of dicts. The `Identifier` can be used as the
catalog for queries about datasets in a certain catalog.  If no parameters are
given, the data is retrieved from the (CBS)[https://www.cbs.nl] OData4 portal.

The optional parameter is:

  - base::String, basename of the OData4 server
"""
function get_catalogs(base=BASE_URL)
    url = base * "/" * "/Catalogs"
    return get_odata(url)["value"]
end

"""
    get_tables(base, catalog)

Get all the tables in the catalog as a list of dicts. The `Identifier` can be
used to get the actual data. If no parameters are given, the data is retrieved
from the (CBS)[https://www.cbs.nl] OData4 portal.

The optional parameters are:

  - base::String, basename of the OData4 server
  - catalog::String, path part of the OData4 server for catalog information

"""
function get_tables(base=BASE_URL, catalog=CATALOG)
    url = base * "/" * catalog * "/Datasets"
    return get_odata(url)["value"]
end

"""
    get_meta(table; base, catalog)

    Get the metadata for a given `table`. The result is a dict with the keys `"Dimensions"`, `"MeasureCodes"`, `"Properties"` and all the classifications (seperate codes and groups).

Optional parameters are:

  - base::String, basename of the OData4 server
  - catalog::String, path part of the OData4 server for regular data
"""
function get_meta(table; base=BASE_URL, catalog=CATALOG)
    tableurl = base * "/" * catalog * "/" * table
    metalinks = get_odata(tableurl)["value"]
    meta = Dict{String,Any}()
    for link in metalinks
        name = link["name"]
        if name in ("Observations", "Properties")
            continue
        end
        linkurl = tableurl * "/" * name
        meta[name] = get_odata(linkurl)["value"]
    end
    propurl = tableurl * "/Properties"
    meta["Properties"] = get_odata(propurl)
    return meta
end

function get_block(url)
    data = get_odata(url)
    nextlink = "@odata.nextLink" in keys(data) ? data["@odata.nextLink"] : nothing
    return (data["value"], nextlink)
end

function get_table(table; base=BASE_URL, catalog=CATALOG, columns=[], filter="")
    if length(columns) == 0
        select = ""
    else
        select = "&\$select=" * join(columns, ",")
    end
    if length(filter) > 0
        filter = "&\$filter=" * HTTP.escape(filter)
    end
    url = base * "/" * catalog * "/" * table * "/Observations?" * select * filter
    data = get_odata(url)
    nextlink = "@odata.nextLink" in keys(data) ? data["@odata.nextLink"] : nothing
    return (data["value"], nextlink)
end

function get_count(table; base=BASE_URL, catalog=CATALOG, columns=[], filter="")
    if length(columns) == 0
        select = ""
    else
        select = "&\$select=" * join(columns, ",")
    end
    if length(filter) > 0
        filter = "&\$filter=" * HTTP.escape(filter)
    end
    url = base * "/" * catalog * "/" * table * "/Observations?\$count=true&\$top=1" * select * filter
    data = get_odata(url)
    return data["@odata.count"]
end

struct ODataTable_long
    table::String
    block::Vector{Any}
    nextlink::Union{Nothing, String}
    names::Vector{Symbol}
    types::Vector{Type}
    length::Int
end

function firstline(t::ODataTable_long)
    datablock = getfield(t, :block)
    nextlink = getfield(t, :nextlink)
    row = 1
    return (OdataRow_long(datablock[row], t), (datablock, nextlink, row))
end

function nextline(t::ODataTable_long, status)
    datablock, nextlink, row = status
    if row == length(datablock) && !isnothing(nextlink)
        datablock, nextlink = get_block(nextlink)
        row = 0
    end
    row += 1
    if row > length(datablock)
        return nothing
    else
        return (OdataRow_long(datablock[row], t), (datablock, nextlink, row))
    end
end

struct OdataRow_long <: Tables.AbstractRow
    row::Dict{String, Any}
    table::ODataTable_long
end

Tables.istable(::ODataTable_long) = true
names(t::ODataTable_long) = getfield(t, :names)
Tables.columnnames(t::ODataTable_long) = getfield(t, :names)
types(t::ODataTable_long) = getfield(t, :types)
Tables.schema(t::ODataTable_long) = Tables.Schema(names(t), types(t))

Tables.rowaccess(::ODataTable_long) = true
Tables.rows(t::ODataTable_long) = t
Base.eltype(t::ODataTable_long) = OdataRow_long
Base.length(t::ODataTable_long) = getfield(t, :length)
Base.iterate(t::ODataTable_long) = firstline(t)
Base.iterate(t::ODataTable_long, st) = nextline(t, st)

Tables.getcolumn(r::OdataRow_long, s::String) = something(getfield(r, :row)[s], missing)
Tables.getcolumn(r::OdataRow_long, i::Int) = something(collect(values(getfield(r, :row)))[getfield(getfield(r, :table), :names)[i]], missing)
Tables.getcolumn(r::OdataRow_long, s::Symbol) = something(getfield(r, :row)[String(s)], missing)

Tables.columnnames(r::OdataRow_long) = collect(keys(getfield(r, :row)))
Base.NamedTuple(r::OdataRow_long) = NamedTuple{Tuple(Symbol.(keys(getfield(r, :row))))}(values(getfield(r, :row)))

"""
    ODataTable_long(table; columns, filters, base, catalog)

Get a `Tables` object for the dataset defined by `table`.
Optional parameters are:

  - columns::Vector{String}, exact columnnames to select, empty for all columns
  - filters::String, OData4 specification for row filter
  - base::String, basename of the OData4 server
  - catalog::String, path part of the OData4 server for catalog information

# Examples
```julia
tbl = CBSOData4.ODataTable_long("82931NED", filter="Measure in ['T001036']")
```
"""
function ODataTable_long(table; columns=[], filter="", base=BASE_URL, catalog=CATALOG)
    recordcount = get_count(table, base=base, catalog=catalog, columns=columns, filter=filter)
    nextblock, nextlink = get_table(table, base=base, catalog=catalog, columns=columns, filter=filter)
    meta = get_meta(table)
    names = Symbol.(collect(keys(nextblock[1])))
    types = [typeof(nextblock[1][String(name)]) for name in names]
    return ODataTable_long(table, nextblock, nextlink, names, types, recordcount)
end

end # module


text:
using NMF
using Parameters

@with_kw struct CoherenceNMF
  normalize::Bool = true
  maxiter::Int = 2000
  tol::Float64 = 1e-4
  warn_converge::Bool = true
  reuse_solution::Bool = false
end
CoherenceMethod(::Type{Val{:nmf}},params) = CoherenceNMF(;params...)

Base.eltype(::CoherenceNMF,::Type) = Float64
Base.eltype(::CoherenceNMF,::AbstractArray) = Float64

function with_method(foreach_window,method::CoherenceNMF,K)
  convergence_count,total_count = 0,0
  Winit,Hinit = fill(0.0,(0,0)),fill(0.0,(0,0))
  method_ = NMF.ALSPGrad{Float64}(tol=method.tol,
                                 maxiter=method.maxiter)

  foreach_window() do window
    x_t = reshape(window,size(window,1)*size(window,2),:)
    k = min(size(x_t,1),K)
    nmf_input = abs.(x_t) .+ method.tol/2

    if isempty(Winit) || !method.reuse_solution
      Winit,Hinit = NMF.nndsvd(nmf_input,k)
    end

    solution = NMF.solve!(method_,nmf_input,Winit,Hinit)
    total_count += 1

    # solution = nnmf(abs.(x_t) .+ method.tol/2,k,init=:nndsvd,
    #                 tol=method.tol,
    #                 maxiter=method.maxiter)

    convergence_count += !solution.converged

    y = reshape(solution.H,:,size(window)[3:end]...)
    permutedims(y,[2:ndims(window)-1;1])
  end

  if convergence_count > 0 && method.warn_converge
    percent = round(100convergence_count / total_count,digits=2)
    if percent > 0
      @info("$percent% of frames (n = $convergence_count) failed to fully "*
            "converge to a solution.")
    else
      @info("<0.01% of frames (n = $convergence_count) failed to fully "*
            "converge to a solution.")
    end
  end
end

ext: jl

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
using Kalman
using Base.Test

include("bkftest.jl")
include("ekftest.jl")
include("ukf/augmented.jl")

if VERSION >= v"0.4.0-"
    include("nullabletest.jl")
end


text:
"""
Yaroslav Salii, 2020

This script reads every .sop (SOPLIB06) from `./sopDir`,
and transforms them into .sop (TSPLIB) in `./outDir`
CAVEAT: this script does not check if precedence constraints are *transitively closed*,
although is mandated by TSPLIB standard. SOPLIB06 is *not* transitively closed.
"""

using DelimitedFiles #just to read the ingress SOPLIB06 .sop files

sopDir = "SOPLIB"
outDir = "t-SOPLIB"
outPath = joinpath(".",outDir)
rm(outPath, recursive=true, force=true) #trash the previous output
mkpath(outPath) #create the output directory

#get all filenames from ./$sopDir that end in "sop"
sops = filter(s -> s[end-2:end]=="sop", readdir(sopDir))
#now just stick all the working code into the *for* loop. Feels so dirty!
for ifName in sops
    ifPath = joinpath(".",sopDir,ifName)
    rawA = readdlm(ifPath,'\t') #SOPLIB06 are CSVs with '\t';
    #Transform: first drop the empty column (the last one, from dangling delimiters)
    A = map(Int,rawA[1:end,1:end-1]);
    #Transform: add the missing TSPLIB-SOP quirk, depot to terminal = 1000000
    A[1,end] = 1000000
    #Transform: reduce each row by concat-space
    B=[ reduce((s, z) -> "$s $z",map(string,A[row,1:end])) for row=1:size(A,1)]
    #Transform: intersperse the rows with newlines '\n'
    C = reduce( (s,z) -> s*'\n'*z,B)

    #Prep. output: collect TSPLIB-like stats in a string
    comment = "Part of SOPLIB06 by Roberto Montemanni"
    dim = size(A,1) #how many vertices, including depot and terminal
    #careful: make sure $ifName is ONLY file name, no directories, no paths
    #unimplemented quirk: three spaces' indent before each row
    legend ="""
    NAME: $ifName
    TYPE: SOP
    DIMENSION: $dim
    EDGE_WEIGHT_TYPE: EXPLICIT
    EDGE_WEIGHT_FORMAT: FULL_MATRIX
    EDGE_WEIGHT_SECTION
    $dim
    $C
    EOF
    """

    # Writing prep
    ofName = ifName
    outFilePath = joinpath(outPath,ofName)

    # Write $legend into $outFilePath
    open(outFilePath,"w") do outFile
        write(outFile,legend)
    end
end


text:
using StatsBase
using DataFrames
using CSV
include("SupportFunction.jl")
include("DataLoad.jl")
println("Leeeeroooy Jenkins")

path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/Speed Benchmark Proffs/New/"
cd(path)

minValArr = zeros(6,3)
minIndexArr = zeros(6,3)

for i = 1:6
    lookUpVal = [1,101,201,301,401,501]
    normData = CSV.read(string(lookUpVal[i])*"bestNorm.csv",delim=",", nullable = false)
    warmstartData = CSV.read(string(lookUpVal[i])*"bestWarm.csv",delim=",", nullable = false)
    mathData = CSV.read(string(lookUpVal[i])*"bestMath.csv",delim=",", nullable = false)

    normMin = findmin(Array(normData))
    warmMin = findmin(Array(warmstartData))
    mathMin = findmin(Array(mathData))

    if normMin[2] > 13
        normIndex = normMin[2]%13
        normMinVal = normMin[1]
    else
        normIndex = normMin[2]
        normMinVal = normMin[1]
    end

    if warmMin[2] > 13
        warmIndex = warmMin[2]%13
        warmMinVal = warmMin[1]
    else
        warmIndex = warmMin[2]
        warmMinVal = warmMin[1]
    end

    if mathMin[2] > 13
        mathIndex = mathMin[2]%13
        mathMinVal = mathMin[1]
    else
        mathIndex = mathMin[2]
        mathMinVal = mathMin[1]
    end

    minValArr[i,1] = normMinVal
    minValArr[i,2] = warmMinVal
    minValArr[i,3] = mathMinVal

    minIndexArr[i,1] = normIndex
    minIndexArr[i,2] = warmIndex
    minIndexArr[i,3] = warmIndex
end

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
using JLD2, Test, FileIO

using CodecZlib, CodecBzip2, CodecLz4

# This is for testing the different syntax versions as well as the library
@testset "Compression with CodecZlib" begin
    fn = joinpath(mktempdir(), "test.jld2")

    randomdata = repeat(rand(2000), 10)
    @save fn {compress=ZlibCompressor()} randomdata
    r = jldopen(f -> f["randomdata"], fn, "r")
    @test r == randomdata

    jldopen(fn, "w"; compress=ZlibCompressor()) do f
        f["randomdata"] = randomdata
    end
    @test load(fn, "randomdata") == randomdata

    jldopen(fn, "w") do f
        write(f, "randomdata", randomdata; compress=ZlibCompressor())
    end
    @test load(fn, "randomdata") == randomdata
end

@testset "Compression with CodecBzip2" begin
    fn = joinpath(mktempdir(), "test.jld2")

    randomdata = repeat(rand(2000), 10)
    @save fn {compress=Bzip2Compressor()} randomdata

    r = jldopen(f -> f["randomdata"], fn, "r")
    @test r == randomdata
end

@testset "Compression with CodecLz4" begin
    fn = joinpath(mktempdir(), "test.jld2")

    randomdata = repeat(rand(2000), 10)
    @save fn {compress=LZ4FrameCompressor()} randomdata

    r = jldopen(f -> f["randomdata"], fn, "r")
    @test r == randomdata
end

@testset "Verify Correctness of Compressor" begin
    fn = joinpath(mktempdir(), "test.jld2")

    @test_throws ArgumentError jldopen(fn, "w"; compress = 42)

    jldopen(fn, "w") do f
        @test_throws ArgumentError write(f, "x", zeros(10); compress = π)
    end
end


text:
using ArgParse

using PkgTemplates
import PkgTemplates: Template
using PkgTemplates: priority, prehook, hook, posthook
include("myplugins.jl")

t_with_jupyter = Template(;
    dir = pwd(),
    julia = v"1.6",
    plugins = [
        License(; name = "MIT"),
        Git(;
            ignore = [
                "docs/build/",
                "docs/site/",
                ".ipynb_checkpoints",
                "*.ipynb",
                "*.gif",
                ".DS_Store",
            ],
            manifest = false,
            ssh = true
        ),
        GitHubActions(;
            extra_versions = ["1.6", "1.7", "nightly"]
        ),
        Documenter{GitHubActions}(),
        Readme(;
            inline_badges = true,
            badge_order = DataType[
                GitHubActions,
                Documenter{GitHubActions},
            ]
        ),
        Dockerfile(with_jupyter = true),
        DockerCompose(with_jupyter = true),
        Makefile(),
        DevContainer(),
        Jupytext(),
        PlaygroundNotebook(),
        PlaygroundPluto(),
        VSCodeExtensions(),
    ]
)

t = Template(;
    dir = pwd(),
    julia = v"1.6",
    plugins = [
        License(; name = "MIT"),
        Git(;
            ignore = [
                "docs/build/",
                "docs/site/",
                #".ipynb_checkpoints",
                #"*.ipynb",
                "*.gif",
                ".DS_Store",
            ],
            manifest = false, ssh = true),
        GitHubActions(;
            extra_versions = ["1.6", "1.7", "nightly"]
        ),
        Documenter{GitHubActions}(),
        Readme(;
            inline_badges = true,
            badge_order = DataType[
                GitHubActions,
                Documenter{GitHubActions},
            ]
        ),
        Dockerfile(with_jupyter = false),
        DockerCompose(with_jupyter = false),
        Makefile(),
        DevContainer(),
        VSCodeExtensions(),
    ]
)

function main()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "pkgname"
        help = "Pkg name"
        required = true
        "--with-jupyter"
        help = "create template for Jupyter workflow"
        action = :store_true
    end
    parsed_args = parse_args(ARGS, s)
    pkgname = parsed_args["pkgname"]
    if !endswith(pkgname, ".jl")
        pkgname *= ".jl"
    end
    if parsed_args["with-jupyter"]
        t_with_jupyter(pkgname)
    else
        t(pkgname)
    end
    nojl = split(pkgname, ".jl")[begin]
    @info "Rename $nojl => $pkgname"
    mv(nojl, pkgname)
end

main()


text:
@time @testset "utils" begin
    @testset "check_sizes" begin
        a = Octavian.StaticInt{1}()
        b = Octavian.StaticInt{2}()
        @test Octavian._select(a, a) === a
        @test Octavian._select(a, 1) === a
        @test Octavian._select(1, a) === a
        @test Octavian._select(1, 1) === 1
        @test_throws AssertionError Octavian.matmul_sizes(rand(3,2), rand(3,4), rand(5,2))
        @test_throws AssertionError Octavian.matmul_sizes(rand(3,2), rand(3,4), rand(4,5))
        @test_throws AssertionError Octavian.matmul_sizes(rand(3,2), rand(5,4), rand(4,2))
    end
end


text:
module FLOWMath

include("quadrature.jl")
export trapz

include("roots.jl")
export brent

include("smooth.jl")
export abs_smooth
export ksmax, ksmin
export ksmax_adaptive, ksmin_adaptive
export sigmoid
export sigmoid_blend
export cubic_blend
export quintic_blend

include("interpolate.jl")
export Akima
export derivative
export gradient
export akima
export linear
export interp2d
export interp3d
export interp4d

end # module


text:
# Convert all problems from JLD format into MATLAB's MAT format

# using JLD, MAT
fromPath = "../DataFiles/Julia/"
toPath = "../DataFiles/MATLAB/"
existingFolders = readdir(toPath)
problemTypes = []
for f in filter(x -> !startswith(x, "."), readdir(fromPath))
    f = split(f,".")[1]
    push!(problemTypes,String(f))
 end
 filter!(x->in(x,["RobustControl"]),problemTypes)

for pType in problemTypes
    println(">>> Current ProblemType: $(pType).")
    subDirPath = fromPath * pType *"/"

    problems = []
    for f in filter(x -> endswith(x, ".jld"), readdir(subDirPath))
        f = split(f,".")[1]
        push!(problems,String(f))
    end

    for p in problems
      # load data from JLD
      data = JLD.load("$(subDirPath)"*"$(p).jld")

      # check if folder already exists in MATLAB dir, otherwise create new
      if !in(pType,existingFolders)
        mkdir(toPath*pType*"/")
        push!(existingFolders,pType)
        end
        mfile = matopen(toPath*pType*"/"*p*".mat", "w")
      # iterate over all keys available
      for k in collect(keys(data))
        write(mfile,k,data[k])
      end
      close(mfile)
      println(" "^8*"$(p) completed!")

    end

end


text:
module HeresValidator

using SearchLight, SearchLight.Validation

function not_empty(field::Symbol, m::T, args::Vararg{Any})::ValidationResult where {T<:AbstractModel}
  isempty(getfield(m, field)) && return ValidationResult(invalid, :not_empty, "should not be empty")

  ValidationResult(valid)
end

end

max_stars_repo_name: huaxk/rlservice

max_issues_repo_name: huaxk/rlservice

max_forks_repo_name: huaxk/rlservice


text:
# Various color scales.

"""
    function weighted_color_mean(cs::AbstractArray{Lab{T},1},
                                 ws::AbstractArray{S,1}) where {S <: Number,T}
Return the mean of Lab colors `cs` as weighted by `ws`.
"""
function weighted_color_mean(
        cs::AbstractArray{Lab{T},1}, ws::AbstractArray{S,1}) where {S <: Number,T}
    l = 0.0
    a = 0.0
    b = 0.0
    sumws = sum(ws)
    for (c, w) in zip(cs, ws)
        w /= sumws
        l += w * c.l
        a += w * c.a
        b += w * c.b
    end
    Lab(l, a, b)
end

# Discrete scales

"""
    lab_rainbow(l, c, h0, n)

Generate `n` colors in the LCHab colorspace by using a fixed
luminance `l` and chroma `c`, and varying the hue, starting at `h0`.
"""
lab_rainbow(l, c, h0, n) = [LCHab(l, c, h0 + 360.0 * (i - 1) / n) for i in 1:n]

"""
    luv_rainbow(l, c, h0, n)

Generate `n` colors in the LCHuv colorspace by using a fixed
luminance `l` and chroma `c`, and varying the hue, starting at `h0`.
"""
luv_rainbow(l, c, h0, n) = [LCHuv(l, c, h0 + 360.0 * (i - 1) / n) for i in 1:n]

# Helpful for Experimenting
function plot_color_scale(colors::Vector{T}) where T <: Color
    println(colors)
    canvas(UnitBox(length(colors), 1)) <<
            (compose([rectangle(i-1, 0, 1, 1) << fill(c)
                      for (i, c) in enumerate(colors)]...) << stroke(nothing))
end

# Continuous scales

"""
    function lab_gradient(cs::Color...)

Generate a function `f(p)` that creates a gradient between n≥2 colors, where `0≤p≤1`.
If you have a collection of colors, then use the splatting operator `...`:
```julia
f = Scale.lab_gradient(range(HSV(0,1,1), stop=HSV(250,1,1), length=100)...)
```
Function `f` can be used like so: `Scale.color_continuous(colormap=f)`.
"""
function lab_gradient(cs::Color...)
    length(cs) < 2 && error("Two or more colors are needed for gradients")

    cs_lab = [convert(Lab, c) for c in cs]
    n = length(cs_lab)
    function f(p::Float64)
        @assert 0.0 <= p <= 1.0
        i = 1 + min(n - 2, max(0, floor(Int, p*(n-1))))
        w = p*(n-1) + 1 - i
        weighted_color_mean([cs_lab[i], cs_lab[i+1]], [1.0 - w, w])
    end
    f
end

"""
    function lab_gradient(cs...)

Can be applied to other types, e.g. `Scale.lab_gradient("blue","ghostwhite","red")`
"""
lab_gradient(cs...) = lab_gradient(Gadfly.parse_colorant(cs)...)

"""
    function lchabmix(c0_, c1_, r, power)
"""
function lchabmix(c0_, c1_, r, power)
    c0 = convert(LCHab, c0_)
    c1 = convert(LCHab, c1_)
    w = r^power
    lspan = c1.l - c0.l
    cspan = c1.c - c0.c
    hspan = c1.h - c0.h
    return LCHab(c0.l + w*lspan, c0.c + w*cspan, c0.h + w*hspan)
end


text:
# coding: utf-8

# In[1]:

#check if a string is a palindrome
#reversing a string is extremely easy
#at first i was thinking about using pop function
#however, for an empty string, pop function is still functional
#it would not stop itself
#it would show errors of popping from empty list
#so i have to use the classic way, slicing
#each time we slice the final letter
#to remove the unnecessary marks, we have to use regular expression
function palindrome_slicing(rawtext)

    #regex is very different in julia
    #we need to use match attribute to get string
    text=[i.match for i in eachmatch(r"[a-zA-Z0-9]",lowercase(rawtext))]

    #use non recursive slicing
    if text[length(text):-1:1]==text

        return true

    else

        return false

    end

end

# In[2]:

function palindrome_recursion(text)

    #this is the base case, when string is empty
    #we just return empty
    if isempty(text)

        return text

    else

        #such a pain without a python style list+list
        return cat([text[end]],
            palindrome_recursion(text[1:end-1]),
            dims=1)

    end

end

# In[3]:

function recursion_check(rawtext)

    #this part is the regular expression to get only characters
    #and format all of em into lower cases
    text=[i.match for i in eachmatch(r"[a-zA-Z0-9]",lowercase(rawtext))]

    if palindrome_recursion(text)==text

        return true

    else

        return false

    end

end

# In[4]:

#or creating a new list
#using loop to append popped items from the old list
function palindrome_list(rawtext)

    text=[i.match for i in eachmatch(r"[a-zA-Z0-9]",lowercase(rawtext))]

    new=String[]

    for i in 1:length(text)

        #be careful with append
        #append will create char type instead of string type
        push!(new,text[end-i+1])

    end

    if new==text

        return true

    else

        return false

    end

end

# In[5]:

#alternatively, we can use deque
#we pop from both sides to see if they are equal
#if not return false
#note that the length of string could be an odd number
#we wanna make sure the pop should take action while length of deque is larger than 1
function palindrome_deque(rawtext)

    text=[i.match for i in eachmatch(r"[a-zA-Z0-9]",lowercase(rawtext))]

    while length(text)>=1

        if pop!(text)!=popfirst!(text)

            return false

        end

    end

    return true

end

# In[6]:

#0.257236 seconds (636.24 k allocations: 31.108 MiB, 2.37% gc time)
@time begin

    recursion_check("Evil is a deed as I live!")

end

# In[7]:

#0.058568 seconds (87.67 k allocations: 4.231 MiB)
@time begin

    palindrome_slicing("Evil is a deed as I live!")

end

# In[8]:

#0.084145 seconds (140.27 k allocations: 6.717 MiB, 8.51% gc time)
@time begin

    palindrome_list("Evil is a deed as I live!")

end

# In[9]:

#0.062403 seconds (74.95 k allocations: 3.500 MiB)
@time begin

    palindrome_deque("Evil is a deed as I live!")

end

# In[10]:

#in julia, the fastest is still the slicing approach

max_stars_repo_path: palindrome checker 4 methods.jl

max_issues_repo_path: palindrome checker 4 methods.jl

max_forks_repo_path: palindrome checker 4 methods.jl


text:
function tplot(P, df, args...;
			   figure = (;),
			   axis_attr = (;),
			   legend_attr = (;),
			   title = nothing,
			   attr_var_pairs...)

	fig = Figure(; figure...)

    # 1. Preparations
    spec = specification(P, df, args, attr_var_pairs)

    # 2. Figure out positions of axis and legend
	if !haskey(legend_attr, :position)
		legend_attr = (; position = :top, legend_attr...)
	end

	legend_position = legend_attr.position

	@unpack ax_pos, legs_pos =
			outer_legend_position(fig, has_legend_or_colorbar(spec), legend_position)

	draw!(ax_pos, legs_pos, spec, attr_var_pairs; legend_attr, axis_attr)

	# 5. Add Title
	if !isnothing(title)
		Label(fig[0,:], title, tellwidth = false, tellheight = true)
	end

	(; fig)
end

function draw!(ax_pos, legs_pos, spec, attr_var_pairs; legend_attr = (;), axis_attr = (;))

	# 3. Draw Axis/Axes
	draw_axis!(ax_pos, spec, axis_attr)

	# 4. Draw Legend/Colorbar
	add_legend(legs_pos, spec; legend_attr...)

	nothing
end

function lplot(args...; kwargs...)
	@unpack fig = tplot(args...; kwargs...)
	fig
end

max_stars_repo_name: greimel/TabularMakie

max_issues_repo_name: greimel/TabularMakie

max_forks_repo_name: greimel/TabularMakie


text:
using PackageEvaluator

function testAll(limit = Inf, writeJSON=true, writeHTML=false)

  if writeHTML
    summary = open("index.html","w")

    # Setup the index file by copying in a header from a file in
    # the extra/ folder
    preindex_path = joinpath(Pkg.dir("PackageEvaluator"),"extra","preindex.inc")
    preindex = open(preindex_path,"r")
    write(summary, readall(preindex))
    close(preindex)
  end

  # Walk through each package in METADATA (assume updated)
  available_pkg = Pkg.available()
  done = 0
  for pkg_name in available_pkg
    println("##### Current package: $pkg_name")

    # Check to see if already have JSON
    if writeJSON && !writeHTML && isfile(joinpath(Pkg.dir("PackageEvaluator"),"extra","$(pkg_name).json"))
      println("      ?????? Already have JSON, skipping")
      continue
    end

    # Run any preprocessing
    exceptions_before(pkg_name)

    # Run PackageEvaluator
    features = nothing
    try
      features = evalPkg(pkg_name, true)  # add and remove it
    catch
      # Couldn't process package, about!
      println("      !!!!!! evalPkg failed")
      continue
    end

    # Run any postprocessing
    exceptions_after(pkg_name)

    # Write a row
    if writeHTML
      write(summary, "
\n")
      # Package
      write(summary, "  $pkg_name	\n")
      # Repo
      write(summary, "  $(pkg_name).jl	\n")
      # License
      if features[:LICENSE] == "Unknown"
        write(summary, "  ?	")
      else
        write(summary, "  $(features[:LICENSE])	\n")
      end

      # Testing
      write(summary, "  	")
      details = getDetailsString(pkg_name, features)
      write(summary, "  $details	\n")

      # Pkg req
      if features[:REQUIRE_VERSION]
        write(summary, "  	")
      else
        write(summary, "  	")
      end
      # META req
      if features[:REQUIRES_OK]
        write(summary, "  	")
      else
        write(summary, "  	")
      end
      # TravisCI
      if features[:TEST_TRAVIS]
        write(summary, "  	")
      else
        write(summary, "  	")
      end

      write(summary, "
\n")
    end

    # WriteJSON
    if writeJSON
      json_str = featuresToJSON(pkg_name, features)
      json_fp = open(joinpath(Pkg.dir("PackageEvaluator"),"extra","$(pkg_name).json"),"w")
      write(json_fp, json_str)
      close(json_fp)
    end

    # Limit number of packages to test
    done = done + 1
    if done >= limit
      break
    end
  end

  if writeHTML
    write(summary, "

")
    close(summary)
  end
end

###############################################################################
# exceptions_before
# Any "special" testing commands to be done
function exceptions_before(pkg_name)
  if pkg_name == "JuMP"
    Pkg.add("Cbc")
    Pkg.add("Clp")
  end
end

# exceptions_after
# Any cleanup for packages that had exceptions
function exceptions_after(pkg_name)
  if pkg_name == "JuMP"
    Pkg.rm("Cbc")
    Pkg.rm("Clp")
  end
end

testAll(1,true,false)


text:
# Computes TTV coefficients for first-order eccentricity
# solution from Agol & Deck (2015).  Please cite this paper
# if you make use of this in your research.

include("laplace_wisdom.jl")

u(gamma::T,c1::T,c2::T) where {T<:Number}= ((3+gamma*gamma)*c1+2*gamma*c2)/(gamma*gamma*(1-gamma*gamma))
# m=+/-1
v(z::T,d1::T,d2::T,m::Integer) where {T<:Number}= ((m*(1-z*z)+6*z)*d1+(2+z*z)*d2)/(z*(1-z*z)*(z+m)*(z+2*m))

function ttv_succinct!(jmax::Integer,alpha::T,f1::Array{T,2},f2::Array{T,2}) where T<:Real

  # See simple_solution.pdf 7/16/2015

  b=zeros(T,jmax+2,3)  #check memory allocation >>>>>>>>>>>>
  #println("Computing Laplace Coeffcients: ",alpha)
  for i=0:2
    for j=0:jmax
      b[j+1,i+1]=laplace_wisdom(1//2,j,i,alpha)/alpha^i
    end
  end
  #println("Done with Laplace Coeffcients: ",alpha)

  sqrtalpha = sqrt(alpha)

  # Loop over j:
  @inbounds for j=0:jmax
    # \delta_{j1} (this is indirect coefficient which is only needed for j=1)
    dj1 = j==1 ? 1.0 : 0.0

    # Compute dimensionless frequencies (equation 30):
    beta = j*(1-alpha*sqrtalpha)
    kappa =  beta /  (alpha*sqrtalpha)

    # Compute disturbing function coefficients (equation 31):
    A_j00 = b[j+1,1]
    A_j10 =  alpha* b[j+1,2]
    A_j01 = -(A_j10 + A_j00)
    A_j20 =  alpha*alpha * b[j+1,3]
    A_j11 = -(2*A_j10 + A_j20)
    A_j02 = 2*A_j00 + 4*A_j10 + A_j20
    jd=convert(eltype(alpha),j)
    # Inner planet coefficients, in order k=0,-1,1,-2,2 (see Table 1):
    if j >=2
      f1[j+1,1]=alpha*u(beta          ,jd*(    A_j00-alpha*dj1),A_j10-alpha*dj1)
      f1[j+1,2]=alpha*u(beta-1.0      ,jd*(-jd*A_j00-0.5*A_j10+1.5*alpha*dj1),-jd*A_j10-0.5*A_j20+alpha*dj1)
      f1[j+1,3]=alpha*u(beta+1.0      ,jd*( jd*A_j00-0.5*A_j10-0.5*alpha*dj1), jd*A_j10-0.5*A_j20-alpha*dj1)
      f1[j+1,4]=alpha*u(beta-alpha*sqrtalpha,jd*( jd*A_j00-0.5*A_j01-2.0*alpha*dj1), jd*A_j10-0.5*A_j11-2.0*alpha*dj1)
      f1[j+1,5]=alpha*u(beta+alpha*sqrtalpha,jd*(-jd*A_j00-0.5*A_j01),-jd*A_j10-0.5*A_j11)
    else
      if j==0
        f1[j+1,4]=alpha*u(beta-alpha*sqrtalpha,jd*( jd*A_j00-0.5*A_j01-2.0*alpha*dj1), jd*A_j10-0.5*A_j11-2.0*alpha*dj1)
      else
        f1[j+1,1]=alpha*u(beta          ,jd*(    A_j00-alpha*dj1),A_j10-alpha*dj1)
        f1[j+1,2]=alpha*u(beta-1.0      ,jd*(-jd*A_j00-0.5*A_j10+1.5*alpha*dj1),-jd*A_j10-0.5*A_j20+alpha*dj1)
        f1[j+1,3]=alpha*u(beta+1.0      ,jd*( jd*A_j00-0.5*A_j10-0.5*alpha*dj1), jd*A_j10-0.5*A_j20-alpha*dj1)
        f1[j+1,4]=alpha*u(beta-alpha*sqrtalpha,jd*( jd*A_j00-0.5*A_j01-2.0*alpha*dj1), jd*A_j10-0.5*A_j11-2.0*alpha*dj1)
      end
    end
    # Add in the k=\pm 1 coefficients (note that d1 & d2 are the same as c1 & c2 for k=0):
    if j >= 1
      f1[j+1,2]=f1[j+1,2]+alpha*v(beta,jd*(A_j00-alpha*dj1),A_j10-alpha*dj1,-1)
      f1[j+1,3]=f1[j+1,3]+alpha*v(beta,jd*(A_j00-alpha*dj1),A_j10-alpha*dj1, 1)
    end

  # Now for the outer planet:
  # Outer planet coefficients, in order k=0,-2,2,-1,1 (see Table 1):
    one_over_alpha_squared = 1/(alpha*alpha)
    if j >= 2
      f2[j+1,1]=u(kappa,-jd*(A_j00-dj1*one_over_alpha_squared),A_j01-dj1*one_over_alpha_squared)
      f2[j+1,2]=u(kappa-1,-jd*(jd*A_j00-0.5*A_j01-0.5*dj1*one_over_alpha_squared),jd*A_j01-0.5*A_j02-dj1*one_over_alpha_squared)
      f2[j+1,3]=u(kappa+1,-jd*(-jd*A_j00-0.5*A_j01+1.5*dj1*one_over_alpha_squared),-jd*A_j01-0.5*A_j02+dj1*one_over_alpha_squared)
      f2[j+1,4]=u(kappa-1/(alpha*sqrtalpha),-jd*(-jd*A_j00-0.5*A_j10),-jd*A_j01-0.5*A_j11)
    else
      if j == 1
        f2[j+1,1]=u(kappa,-jd*(A_j00-dj1*one_over_alpha_squared),A_j01-dj1*one_over_alpha_squared)
        f2[j+1,2]=u(kappa-1,-jd*(jd*A_j00-0.5*A_j01-0.5*dj1*one_over_alpha_squared),jd*A_j01-0.5*A_j02-dj1*one_over_alpha_squared)
        f2[j+1,3]=u(kappa+1,-jd*(-jd*A_j00-0.5*A_j01+1.5*dj1*one_over_alpha_squared),-jd*A_j01-0.5*A_j02+dj1*one_over_alpha_squared)
      end
    end
    f2[j+1,5]=u(kappa+1/(alpha*sqrtalpha),-jd*(jd*A_j00-0.5*A_j10-2.0*dj1*one_over_alpha_squared),jd*A_j01-0.5*A_j11-2.0*dj1*one_over_alpha_squared)
  # Add in the k=\pm 2 coefficients (note that d1 & d2 are the same as c1 & c2 for k=0):
    if j >= 1
      f2[j+1,2]=f2[j+1,2]+v(kappa,-jd*(A_j00-dj1*one_over_alpha_squared),A_j01-dj1*one_over_alpha_squared,-1)
      f2[j+1,3]=f2[j+1,3]+v(kappa,-jd*(A_j00-dj1*one_over_alpha_squared),A_j01-dj1*one_over_alpha_squared, 1)
    end
  # That's it!
  end #end of loop
  return
end

# Computes Taylor series expansion for when period is fixed. code exists somewhere
function ttv_succinct!(jmax::Integer,alpha::Number,f1::Array{Float64,2},f2::Array{Float64,2},b::Array{Float64,2},alpha0::Number,b0::Array{Float64,2})

  # See simple_solution.pdf 7/16/2015

  # Fourth-order Taylor expansion approximation of Laplace coefficients:
  dalpha = alpha-alpha0
  for i=0:2
    for j=0:jmax
      b[j+1,i+1]=b0[i+1,j+1]+dalpha*(b0[i+2,j+1]+0.5*dalpha*(b0[i+3,j+1]+dalpha/3.0*(b0[i+4,j+1]+dalpha*0.25*b0[i+5,j+1])))
    end
  end

  sqrtalpha = sqrt(alpha)

  # Loop over j:
  @inbounds for j=0:jmax
    # \delta_{j1} (this is indirect coefficient which is only needed for j=1)
    dj1 = j==1 ? 1.0 : 0.0

    # Compute dimensionless frequencies (equation 30):
    beta = j*(1-alpha*sqrtalpha)
    kappa =  beta /  (alpha*sqrtalpha)

    # Compute disturbing function coefficients (equation 31):
    A_j00 = b[j+1,1]
    A_j10 =  alpha* b[j+1,2]
    A_j01 = -(A_j10 + A_j00)
    A_j20 =  alpha*alpha * b[j+1,3]
    A_j11 = -(2*A_j10 + A_j20)
    A_j02 = 2*A_j00 + 4*A_j10 + A_j20
    jd=convert(eltype(alpha),j)
    # Inner planet coefficients, in order k=0,-1,1,-2,2 (see Table 1):
    if j >=2
      f1[j+1,1]=alpha*u(beta          ,jd*(    A_j00-alpha*dj1),A_j10-alpha*dj1)
      f1[j+1,2]=alpha*u(beta-1.0      ,jd*(-jd*A_j00-0.5*A_j10+1.5*alpha*dj1),-jd*A_j10-0.5*A_j20+alpha*dj1)
      f1[j+1,3]=alpha*u(beta+1.0      ,jd*( jd*A_j00-0.5*A_j10-0.5*alpha*dj1), jd*A_j10-0.5*A_j20-alpha*dj1)
      f1[j+1,4]=alpha*u(beta-alpha*sqrtalpha,jd*( jd*A_j00-0.5*A_j01-2.0*alpha*dj1), jd*A_j10-0.5*A_j11-2.0*alpha*dj1)
      f1[j+1,5]=alpha*u(beta+alpha*sqrtalpha,jd*(-jd*A_j00-0.5*A_j01),-jd*A_j10-0.5*A_j11)
    else
      if j==0
        f1[j+1,4]=alpha*u(beta-alpha*sqrtalpha,jd*( jd*A_j00-0.5*A_j01-2.0*alpha*dj1), jd*A_j10-0.5*A_j11-2.0*alpha*dj1)
      else
        f1[j+1,1]=alpha*u(beta          ,jd*(    A_j00-alpha*dj1),A_j10-alpha*dj1)
        f1[j+1,2]=alpha*u(beta-1.0      ,jd*(-jd*A_j00-0.5*A_j10+1.5*alpha*dj1),-jd*A_j10-0.5*A_j20+alpha*dj1)
        f1[j+1,3]=alpha*u(beta+1.0      ,jd*( jd*A_j00-0.5*A_j10-0.5*alpha*dj1), jd*A_j10-0.5*A_j20-alpha*dj1)
        f1[j+1,4]=alpha*u(beta-alpha*sqrtalpha,jd*( jd*A_j00-0.5*A_j01-2.0*alpha*dj1), jd*A_j10-0.5*A_j11-2.0*alpha*dj1)
      end
    end
    # Add in the k=\pm 1 coefficients (note that d1 & d2 are the same as c1 & c2 for k=0):
    if j >= 1
      f1[j+1,2]=f1[j+1,2]+alpha*v(beta,jd*(A_j00-alpha*dj1),A_j10-alpha*dj1,-1)
      f1[j+1,3]=f1[j+1,3]+alpha*v(beta,jd*(A_j00-alpha*dj1),A_j10-alpha*dj1, 1)
    end

  # Now for the outer planet:
  # Outer planet coefficients, in order k=0,-2,2,-1,1 (see Table 1):
    one_over_alpha_squared = 1/(alpha*alpha)
    if j >= 2
      f2[j+1,1]=u(kappa,-jd*(A_j00-dj1*one_over_alpha_squared),A_j01-dj1*one_over_alpha_squared)
      f2[j+1,2]=u(kappa-1,-jd*(jd*A_j00-0.5*A_j01-0.5*dj1*one_over_alpha_squared),jd*A_j01-0.5*A_j02-dj1*one_over_alpha_squared)
      f2[j+1,3]=u(kappa+1,-jd*(-jd*A_j00-0.5*A_j01+1.5*dj1*one_over_alpha_squared),-jd*A_j01-0.5*A_j02+dj1*one_over_alpha_squared)
      f2[j+1,4]=u(kappa-1/(alpha*sqrtalpha),-jd*(-jd*A_j00-0.5*A_j10),-jd*A_j01-0.5*A_j11)
    else
      if j == 1
        f2[j+1,1]=u(kappa,-jd*(A_j00-dj1*one_over_alpha_squared),A_j01-dj1*one_over_alpha_squared)
        f2[j+1,2]=u(kappa-1,-jd*(jd*A_j00-0.5*A_j01-0.5*dj1*one_over_alpha_squared),jd*A_j01-0.5*A_j02-dj1*one_over_alpha_squared)
        f2[j+1,3]=u(kappa+1,-jd*(-jd*A_j00-0.5*A_j01+1.5*dj1*one_over_alpha_squared),-jd*A_j01-0.5*A_j02+dj1*one_over_alpha_squared)
      end
    end
    f2[j+1,5]=u(kappa+1/(alpha*sqrtalpha),-jd*(jd*A_j00-0.5*A_j10-2.0*dj1*one_over_alpha_squared),jd*A_j01-0.5*A_j11-2.0*dj1*one_over_alpha_squared)
  # Add in the k=\pm 2 coefficients (note that d1 & d2 are the same as c1 & c2 for k=0):
    if j >= 1
      f2[j+1,2]=f2[j+1,2]+v(kappa,-jd*(A_j00-dj1*one_over_alpha_squared),A_j01-dj1*one_over_alpha_squared,-1)
      f2[j+1,3]=f2[j+1,3]+v(kappa,-jd*(A_j00-dj1*one_over_alpha_squared),A_j01-dj1*one_over_alpha_squared, 1)
    end
  # That's it!
  end
  return
end

max_stars_repo_name: bmlindor/TTVFaster

max_issues_repo_name: bmlindor/TTVFaster

max_forks_repo_name: bmlindor/TTVFaster


text:
###
# This file contains BLAS/BandedMatrix overrides for operators
# that depend on the structure of BandedMatrix
####

# default copy is to loop through
# override this for most operators.
function default_BandedMatrix(S::Operator)
    Y=BandedMatrix{eltype(S)}(undef, size(S), bandwidths(S))

    for j=1:size(S,2),k=colrange(Y,j)
        @inbounds inbands_setindex!(Y,S[k,j],k,j)
    end

    Y
end

# default copy is to loop through
# override this for most operators.
function default_RaggedMatrix(S::Operator)
    data=Array{eltype(S)}(undef, 0)
    cols=Array{Int}(undef, size(S,2)+1)
    cols[1]=1
    for j=1:size(S,2)
        cs=colstop(S,j)
        K=cols[j]-1
        cols[j+1]=cs+cols[j]
        resize!(data,cols[j+1]-1)

        for k=1:cs
            data[K+k]=S[k,j]
        end
    end

    RaggedMatrix(data,cols,size(S,1))
end

function default_Matrix(S::Operator)
    n, m = size(S)
    if isinf(n) || isinf(m)
        error("Cannot convert $S to a Matrix")
    end

    eltype(S)[S[k,j] for k=1:n, j=1:m]
end

# The diagonal of the operator may not be the diagonal of the sub
# banded matrix, so the following calculates the row of the
# Banded matrix corresponding to the diagonal of the original operator

diagindshift(S,kr,jr) = first(kr)-first(jr)
diagindshift(S::SubOperator) = diagindshift(S,parentindices(S)[1],parentindices(S)[2])

#TODO: Remove
diagindrow(S,kr,jr) = bandwidth(S,2)+first(jr)-first(kr)+1
diagindrow(S::SubOperator) = diagindrow(S,parentindices(S)[1],parentindices(S)[2])

#####
# Conversions
#####

function BandedMatrix(S::SubOperator{T,ConcreteConversion{Chebyshev{DD,RR},Ultraspherical{Int,DD,RR},T},
                              Tuple{UnitRange{Int},UnitRange{Int}}}) where {T,DD,RR}
    # we can assume order is 1
    ret = BandedMatrix{eltype(S)}(undef, size(S), bandwidths(S))
    kr,jr = parentindices(S)
    dg = diagindshift(S)

    @assert -bandwidth(ret,1) ≤ dg ≤ bandwidth(ret,2)-2

    ret[band(dg)] .= 0.5
    ret[band(dg+1)] .= 0.0
    ret[band(dg+2)] .= -0.5

    # correct first entry
    if 1 in kr && 1 in jr
        ret[1,1] = 1.0
    end

    ret
end

function BandedMatrix(V::SubOperator{T,ConcreteConversion{Ultraspherical{LT,DD,RR},Ultraspherical{LT,DD,RR},T},
                                                                  Tuple{UnitRange{Int},UnitRange{Int}}}) where {T,LT,DD,RR}

    n,m = size(V)
    V_l, V_u = bandwidths(V)
    ret = BandedMatrix{eltype(V)}(undef, (n,m), (V_l,V_u))
    kr,jr = parentindices(V)
    dg = diagindshift(V)

    λ = order(rangespace(parent(V)))
    c = λ-one(T)

    # need to drop columns

    1-n ≤ dg ≤ m-1 && (ret[band(dg)] .= c./(jr[max(0,dg)+1:min(n+dg,m)] .- 2 .+ λ))
    1-n ≤ dg+1 ≤ m-1 && (ret[band(dg+1)] .= 0)
    1-n ≤ dg+2 ≤ m-1 && (ret[band(dg+2)] .= c./(2 .- λ .- jr[max(0,dg+2)+1:min(n+dg+2,m)]))

    ret
end

#####
# Derivatives
#####

function BandedMatrix(S::SubOperator{T,ConcreteDerivative{Chebyshev{DD,RR},K,T},
                                                     Tuple{UnitRange{Int},UnitRange{Int}}}) where {T,K,DD,RR}

    n,m = size(S)
    ret = BandedMatrix{eltype(S)}(undef, (n,m), bandwidths(S))
    kr,jr = parentindices(S)
    dg = diagindshift(S)

    D = parent(S)
    k = D.order
    d = domain(D)

    C=convert(T,pochhammer(one(T),k-1)/2*(4/(complexlength(d)))^k)

    # need to drop columns

    if 1-n ≤ dg+k ≤ m-1
        ret[band(dg+k)] .= C.*(jr[max(0,dg+k)+1:min(n+dg+k,m)] .- one(T))
    end

    ret
end

function BandedMatrix(S::SubOperator{T,ConcreteDerivative{Ultraspherical{LT,DD,RR},K,T},
                                                  Tuple{UnitRange{Int},UnitRange{Int}}}) where {T,K,DD,RR,LT}
    n,m = size(S)
    ret = BandedMatrix{eltype(S)}(undef, (n,m), bandwidths(S))
    kr,jr = parentindices(S)
    dg = diagindshift(S)

    D = parent(S)
    k = D.order
    λ = order(domainspace(D))
    d = domain(D)

    C = convert(T,pochhammer(one(T)*λ,k)*(4/(complexlength(d)))^k)
    ret[band(dg+k)] .= C

    ret
end

####
# Toeplitz/Hankel
####

# αn,α0,αp give the constant for the negative, diagonal and positive
# entries.  The usual case is 1,2,1
function toeplitz_axpy!(αn,α0,αp,neg,pos,kr,jr,ret)
    dat=ret.data
    m=size(dat,2)

    dg=diagindrow(ret,kr,jr)

    # diagonal
    if dg ≥ 1
        α0p=α0*pos[1]
        @simd for j=1:m
            @inbounds dat[dg,j]+=α0p
        end
    end

    # positive entries
    for k=2:min(length(pos),dg)
        αpp=αp*pos[k]
        @simd for j=1:m
            @inbounds dat[dg-k+1,j]+=αpp
        end
    end

    # negative entries
    for k=1:min(length(neg),size(dat,1)-dg)
        αnn=αn*neg[k]
        @simd for j=1:m
            @inbounds dat[dg+k,j]+=αnn
        end
    end

    ret
end

# this routine is for when we want a symmetric toeplitz
function sym_toeplitz_axpy!(αn,α0,αp,cfs,kr,jr,ret)
    dg=diagindrow(ret,kr,jr)

    dat=ret.data
    m=size(dat,2)
    # diagonal
    if dg ≥ 1
        α0p=α0*cfs[1]
        @simd for j=1:m
            @inbounds dat[dg,j]+=α0p
        end
    end

    # positive entries
    for k=2:min(length(cfs),dg)
        αpp=αp*cfs[k]
        @simd for j=1:m
            @inbounds dat[dg-k+1,j]+=αpp
        end
    end

    # negative entries
    for k=2:min(length(cfs),size(dat,1)-dg+1)
        αnn=αn*cfs[k]
        @simd for j=1:m
            @inbounds dat[dg+k-1,j]+=αnn
        end
    end

    ret
end

toeplitz_axpy!(α,neg,pos,kr,jr,ret) =
    toeplitz_axpy!(α,α,α,neg,pos,kr,jr,ret)

sym_toeplitz_axpy!(α0,αp,cfs,kr,jr,ret) =
    sym_toeplitz_axpy!(αp,α0,αp,cfs,kr,jr,ret)

function hankel_axpy!(α,cfs,kr,jr,ret)
    dat=ret.data

    st=stride(dat,2)-2
    mink=first(kr)+first(jr)-1

    N=length(dat)

    # dg gives the row corresponding to the diagonal of the original operator
    dg=diagindrow(ret,kr,jr)

    # we need the entry where the first entry is written to
    # this is going to be a shift of the diagonal of the true operator
    dg1=dg+first(kr)-first(jr)

    for k=mink:min(length(cfs),size(dat,1)+mink-dg1)
        dk=k-mink+dg1
        nk=k-mink
        αc=α*cfs[k]
        @simd for j=dk:st:min(dk+nk*st,N)
            @inbounds dat[j]+=αc
        end
    end

    ret
end


text:
@testset "Bronstein Problems" begin
    (x, ) = @variables x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Manuel*Bronstein - Symbolic*Integration*Tutorial*(1998)=#

    #= ::Section::Closed:: =#
    #=2  Algebraic*Functions=#

    @test_int [(2*x^8 + 1)*(sqrt(x^8 + 1)/(x^17 + 2*x^9 + x)), x, 6, -(1/(4*sqrt(1 + x^8))) - (1/4)*atanh(sqrt(1 + x^8))]
    @test_int [1/(1 + x^2), x, 1, atan(x)]
    @test_int [sqrt(x^8 + 1)/(x*(x^8 + 1)), x, 3, (-(1/4))*atanh(sqrt(1 + x^8))]
    @test_int [x/sqrt(1 - x^3), x, 3, (2*sqrt(1 - x^3))/(1 + sqrt(3) - x) - (3^(1/4)*sqrt(2 - sqrt(3))*(1 - x)*sqrt((1 + x + x^2)/(1 + sqrt(3) - x)^2)*Elliptic.E(asin((1 - sqrt(3) - x)/(1 + sqrt(3) - x)), -7 - 4*sqrt(3)))/(sqrt((1 - x)/(1 + sqrt(3) - x)^2)*sqrt(1 - x^3)) + (2*sqrt(2)*(1 - x)*sqrt((1 + x + x^2)/(1 + sqrt(3) - x)^2)*Elliptic.F(asin((1 - sqrt(3) - x)/(1 + sqrt(3) - x)), -7 - 4*sqrt(3)))/(3^(1/4)*sqrt((1 - x)/(1 + sqrt(3) - x)^2)*sqrt(1 - x^3))]
    @test_int [1/(x*sqrt(1 - x^3)), x, 3, (-(2/3))*atanh(sqrt(1 - x^3))]
    @test_int [x/sqrt(x^4 + 10*x^2 - 96*x - 71), x, 1, (1/8)*log(10001 + 3124*x^2 - 1408*x^3 + 54*x^4 - 128*x^5 + 20*x^6 + x^8 + sqrt(-71 - 96*x + 10*x^2 + x^4)*(781 - 528*x + 27*x^2 - 80*x^3 + 15*x^4 + x^6))]

    #= ::Section::Closed:: =#
    #=3  Elementary*Functions=#

    @test_int [(x - tan(x))/tan(x)^2, x, 6, -(x^2/2) - x*cot(x)]
    @test_int [1 + x*tan(x) + tan(x)^2, x, 7, (im*x^2)/2 - x*log(1 + ℯ^(2*im*x)) + (1/2)*im*PolyLog(2, -ℯ^(2*im*x)) + tan(x)]
    @test_int [sin(x)/x, x, 1, SinIntegral(x)]
    @test_int [(3*(x + ℯ^x)^(1/3) + (2*x^2 + 3*x)*ℯ^x + 5*x^2)/(x*(x + ℯ^x)^(1/3)), x, 8, 3*x*(ℯ^x + x)^(2/3) + 3*log(x)]

    @test_int [1/x + (1 + 1/x)/(x + log(x))^(3/2), x, 2, log(x) - 2/sqrt(x + log(x))]
    @test_int [(log(x)^2 + 2*x*log(x) + x^2 + (x + 1)*sqrt(x + log(x)))/(x*log(x)^2 + 2*x^2*log(x) + x^3), x, -3, log(x) - 2/sqrt(x + log(x))]

    @test_int [(2*log(x)^2 - log(x) - x^2)/(log(x)^3 - x^2*log(x)), x, 6, (-(1/2))*log(x - log(x)) + (1/2)*log(x + log(x)) + LogIntegral(x)]
    #= [log(1 + ℯ^x)^(1/3)/(1 + log(1 + ℯ^x)), x, 0, CannotIntegrate(log(1 + ℯ^x)^(1/3)/(1 + log(1 + ℯ^x)), x)] =#
    #= [((x^2 + 2*x + 1)*sqrt(x + log(x)) + (3*x + 1)*log(x) + 3*x^2 + x)/((x*log(x) + x^2)*sqrt(x + log(x)) + x^2*log(x) + x^3), x, 0, 2*sqrt(x + log(x)) + 2*log(x + sqrt(x + log(x)))] =#

    #= ::Title:: =#
    #=Manuel*Bronstein - Symbolic*Integration*im; Transcendental*FunctionsTutorial*(2005)=#

    #= ::Section::Closed:: =#
    #=2.8  Rioboo's*Algorithm*for*Real*Rational*Functions=#

    @test_int [(x^4 - 3*x^2 + 6)/(x^6 - 5*x^4 + 5*x^2 + 4), x, 1, -atan(sqrt(3) - 2*x) + atan(sqrt(3) + 2*x) + atan((1/2)*x*(1 - 3*x^2 + x^4))]
    end

max_stars_repo_path: test/0 Independent test suites/Bronstein Problems.jl


text:
module AbstractFFTs

import ChainRulesCore

export fft, ifft, bfft, fft!, ifft!, bfft!,
       plan_fft, plan_ifft, plan_bfft, plan_fft!, plan_ifft!, plan_bfft!,
       rfft, irfft, brfft, plan_rfft, plan_irfft, plan_brfft,
       fftshift, ifftshift, Frequencies, fftfreq, rfftfreq

include("definitions.jl")
include("chainrules.jl")

end # module


text:
using DecomposedMDPSolver
using POMDPs, POMDPModels, POMDPSimulators, POMDPModelTools
using LocalFunctionApproximation, GridInterpolations, LocalApproximationPolicyEvaluation
using Flux
using Random
using BSON

# Step 1 - Setup the gridworld problem
action_probability(mdp::SimpleGridWorld, s, a) = 0.25
POMDPs.convert_s(::Type{AbstractArray}, s::GWPos, mdp::SimpleGridWorld) = SVector{4, Float64}(1., s[1], s[2], s[1]*s[2])
POMDPs.gen(mdp::SimpleGridWorld, s, a, rng = Random.GLOBAL_RNG) = (sp =rand(transition(mdp, s, a )), r=reward(mdp, s, a))
POMDPs.initialstate(mdp::SimpleGridWorld) = rand(initialstate_distribution(mdp))
g_size = (9,9)
g = SimpleGridWorld(size = g_size, rewards = Dict(GWPos(g_size...) => 1, GWPos(1,1) => 0), tprob = 1., discount=1)

# Step 2 - Solve the problem exactly using local approximation
grid = RectangleGrid([1:g_size[1] ...], [1:g_size[2]...])
interp = LocalGIFunctionApproximator(grid)
solver = LocalPolicyEvalSolver(interp, action_probability, is_mdp_generative = true, n_generative_samples = 1,  verbose = true, max_iterations = 2000, belres = 1e-6)
dp_policy = solve(solver, g)
render(g, (s=GWPos(1,1), r=1), color = (s) -> 20. *(value(dp_policy, s) - 0.5))

# Step 3 - Solve the problem using a MC policy eval with function approximation
V_network = Chain((x) -> x .- 5.f0 ./ 5.f0, Dense(2,32, relu), Dense(32, 1, sigmoid))
render(g, (s=GWPos(1,1), r=1), color = (s) -> 20. *(V_network(convert_s(Array{Float64, 1}, s, g))[1] - 0.5))

mc_policy_eval!(g, V_network, action_probability, iterations=100)

render(g, (s=GWPos(1,1), r=1), color = (s) -> 20. *(V_network(convert_s(Array{Float64, 1}, s, g))[1] - 0.5))

new_V = BSON.load("log/_4/last_model.bson")[:V]
render(g, (s=GWPos(1,1), r=1), color = (s) -> 20. *(new_V(convert_s(Array{Float64, 1}, s, g))[1] - 0.5))


text:
using Pkg
Pkg.develop(PackageSpec(path=pwd()))
Pkg.instantiate()

using Documenter, MolecularGraph

makedocs(
    sitename="MolecularGraph.jl",
    pages = [
        "Home" => "index.md",
        "MolecularGraph" => [
            "Molecule" => "moleculargraph/molecule.md",
            "I/O" => "moleculargraph/io.md",
            "Structure drawing" => "moleculargraph/draw.md",
            "Chemical properties" => "moleculargraph/properties.md",
            "Preprocessing" => "moleculargraph/preprocess.md",
            "Molecular mass/weight" => "moleculargraph/mass.md",
            "Structure match" => "moleculargraph/structurematch.md",
            "Molecular queries" => "moleculargraph/query.md"
        ],
        "MolecularGraph.Graph" => "graph.md",
        "MolecularGraph.Geometry" => "geometry.md"
    ],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    )
)

deploydocs(
    repo="github.com/mojaie/MolecularGraph.jl.git",
    push_preview=true
)


text:
# QuObj Tests
using Schrodinger
using Compat, Compat.Test, Compat.LinearAlgebra, Compat.SparseArrays
println("Testing Math...")

@testset "QuObj Basic Math" begin
# Build a few different variable for testing
g = basis(2,0)
e1 = dense(basis(2,1))
ρ = maxmixed(4)
α = rand(ComplexF64)
ψα = coherent(15,α)
β = rand(ComplexF64)
ψβ = coherent(15,β)
σ = dense(create(4))
a4 = destroy(4)
a2 = destroy(2)
adag2 = create(2)
plus = normalize!(g+1.0im*basis(2,1))

@testset "QuObj/Number Algebra" begin
    # Test algebra between QuObjects and numbers
    @test +g == g
    @test -g == 0-g
    @test +e1 == e1
    @test -e1 == 0-e1
    @test +(g') == g'
    @test -(g') == 0-g'
    @test +(e1') == e1'
    @test -(e1') == 0-e1'
    @test +ρ == ρ
    @test -ρ == 0-ρ
    @test +σ == σ
    @test -σ == 0-σ
    @test g+1 == 1+g == Ket([2,1])
    @test e1-1 == Ket([-1,0])
    @test 1-e1 == Ket([1,0])
    @test g'+1 == 1+g' == Bra([2,1])
    @test e1'-1 == Bra([-1,0])
    @test 1-e1' == Bra([1,0])
    @test ρ+1 == 1+ρ == ρ+qeye(4)
    @test ρ-0.25 == 0.25-ρ == qzero(4)
    @test σ+2 == 2+σ == σ+2*qeye(4)
    @test σ-2 == -(2-σ) == σ-2*qeye(4)
    @test 2*g == g*2 == Ket(sparse([2,0]))
    @test 2*e1 == e1*2 == Ket([0,2])
    @test 2*(g') == (g')*2 == Bra(sparse([2,0]))
    @test 2*(e1') == (e1')*2 == Bra([0,2])
    @test 4*ρ == ρ*4 == qeye(4)
    @test 2σ == σ*2 == σ+σ
    @test g/2 == Ket(sparse([0.5,0]))
    @test e1/1 == e1
    @test (g')/2 == Bra(sparse([0.5,0]))
    @test (e1')/1 == e1'
    @test_throws ArgumentError 2/g
    @test_throws ArgumentError 2/e1
    @test_throws ArgumentError 2/(g')
    @test_throws ArgumentError 2/(e1')
    if VERSION > v"0.7.0-"
        @test (x=(g/0);isnan(x[2])&&isinf(x[1])) # see julia PR #22715
    end
    @test (x=(e1/0);isnan(x[1])&&isinf(x[2]))
    @test ρ/2 == qeye(4)/8
    @test (x=ρ/0;isnan(x[2,1])&&isinf(x[1,1]))
    @test (x=σ/0;isnan(x[1,1])&&isinf(x[2,1]))
    @test_throws ArgumentError 3/ρ
    @test_throws ArgumentError 3/σ
    @test_throws ArgumentError plus^3
    @test_throws ArgumentError plus^3.3
    @test ρ^3 == 0.015625*qeye(4)
    @test σx^2 == σy^2 == σz^2 == σ0
    if VERSION < v"0.7.0-"
        @test_throws MethodError data(ρ^2.5)
    else
        @test data(ρ^2.5) == data(ρ)^2.5
    end
    @test data(dense(ρ)^2.5) == full(ρ)^2.5
    @test (σ^2)[4,2] == √(2)*√(3)
    @test_broken data(σ^2.5) == data(σ)^2.5
end

@testset "QuObj/QuObj Algebra" begin
    @test g+g == 2g
    @test g-g == Ket(sparse([0,0]))
    @test e1+e1+e1 == 3*e1
    @test g+e1 == Ket([1,1])
    @test g-e1 == Ket([1,-1])
    @test g'+g' == 2g'
    @test e1'+e1'+e1' == 3*e1'
    @test g'+e1' == Bra([1,1])
    @test g'-e1' == Bra([1,-1])
    @test ρ+ρ == 2ρ
    @test σ+σ+σ == 3σ
    @test σ-σ == qzero(4)
    @test ρ-ρ == qzero(4)
    @test data(ρ+σ) == data(ρ)+data(σ)
    @test data(ρ-σ) == data(ρ)-data(σ)
    @test_throws ArgumentError g+σ
    @test_throws ArgumentError g-σ
    @test ρ*ρ == ρ^2
    @test σ*σ == σ^2
    @test σ*a4 ≈ numberop(4)
    @test a4'*a4 ≈ numberop(4)
    @test σ*σ' ≈ numberop(4)
    @test a4'*σ' ≈ numberop(4)
    @test adag2*g == adag2*Bra(g)'== e1
    @test adag2'*e1 == adag2'*Bra(e1)' == g
    @test g'*a2 == Bra(g)*a2 == e1'
    @test g'*adag2' == Bra(g)*adag2' == e1'
    @test g'*g == Bra(e1)*e1 == 1
    @test dot(plus,plus) ≈ 1
    @test dot(ψα,ψβ) ≈ exp(-(abs2(α)+abs2(β))/2+α'*β)
    @test plus*plus' == plus*Bra(plus) ≈ normalize!(σ0 + σy)
    @test_throws ArgumentError g*e1
    @test_throws ArgumentError g/e1
end

@testset "Operator Inner Product" begin
    @test inner(σx,σx) == 2
    @test inner(σx,σy) == 0
    @test inner(σx,σz) == 0
    @test inner(ρ*projectorop(size(ρ,1),1:2),ρ) == 2/size(ρ,1)^2
    @test inner(Operator(ψα),Operator(ψβ)) ≈ exp(-abs2(α-β))
    @test_throws DimensionMismatch inner(Operator(ψα),ρ)
end

@testset "Tensor Product" begin
    @test g ⊗ g   == Ket(Bra(g)  ⊗ Bra(g)) == Ket(sparse([1,0,0,0]),(2,2))
    @test g ⊗ e1  == Ket(Bra(g)  ⊗ Bra(e1)) == Ket(sparse([0,1,0,0]),(2,2))
    @test e1 ⊗ g  == Ket(Bra(e1) ⊗ Bra(g)) == Ket(sparse([0,0,1,0]),(2,2))
    @test e1 ⊗ e1 == Ket(Bra(e1) ⊗ Bra(e1)) == Ket([0,0,0,1],(2,2))
    @test data(Operator(g) ⊗ σ) == [data(σ) zeros(4,4); zeros(4,8)]
    @test data(ρ ⊗ Operator(g)) == (A=Diagonal(fill(0.25,8));A.diag[2:2:8].=0;A)
    @test_throws ArgumentError e1 ⊗ ρ
    @test_throws ArgumentError Bra(e1) ⊗ σ
end

@testset "Transposition" begin
    @test g' == Bra([1,0])
    @test (g+1im*e1)' == Bra([1,-1im])
    @test conj((g+1im*e1)) == Ket([1,-1im])
    @test Bra(g)' == Ket([1,0])
    @test Bra([1,1im])' == Ket([1,-1im])
    @test conj(Bra([1,1im])) == Bra([1,-1im])
    @test ρ' == transpose(ρ) == ρ
    @test σ' == transpose(σ)
end

@testset "Math Functions" begin
    for f in [sqrt,log,exp], A in [Operator(coherent(4,1.2+3im)),ρ,create(4)+destroy(4)]
        @test data(f(A)) == f(full(A))
    end
end

@testset "Misc Functions" begin
    obj = [g,e1,plus,Bra(plus),σ,ρ,Operator(coherent(4,1.2+3im))]
    for A in obj
        for f in [real,imag,abs,abs2]
            @test f(A) == f.(data(A))
        end
        @static if VERSION > v"0.7.0-"
            @test data(round(A,digits=4)) == round.(data(A),digits=4)
        else
            @test data(round(A,4)) == round.(data(A),4)
        end
    end
end

@testset "Eigen" begin
    E = 3.3; Δ = 0.25
    O = dense(E*σz+Δ*σx)
    @test eigvals(O) == [-√(E^2+Δ^2), √(E^2+Δ^2)]
    @test eigen(O)[2][1] ≈ normalize!(Ket([-(E-√(E^2+Δ^2))/Δ,-1]))
    a = destroy(20) ⊗ qeye(2)
    s = qeye(20) ⊗ σ₋
    fc = 5; fa = 6; g = 0.025;
    H = fc*(a'a+0.5) + fa*(s's-0.5) + g*(a's + a*s')
    function En(l)
        n,q,Δ = (l-1)÷2, 1-2*(l%2), fa-fc
        return l==0 ? -Δ/2 : (n+1)*fc + 0.5q*√(4g^2*(n+1)+Δ^2)
    end
    θ(n) = 0.5*atan(2g*√(n+1)/(fa-fc))
    D,V = eigs(H,nev=10)
    @test D ≈ En.(0:9)
    @test abs(inner(V[6], cos(θ(2))*ket((3,0),(20,2)) - sin(θ(2))*ket((2,1),(20,2)))) ≈ 1
end

end


text:
using Pkg
Pkg.activate("Queller", io=devnull)
Pkg.instantiate(io=devnull)

using Queller
main()

max_stars_repo_name: mvmorin/queller-bot

max_issues_repo_name: mvmorin/queller-bot

max_forks_repo_name: mvmorin/queller-bot


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.api.concept.type;
#
# import typedb.client.api.TypeDBTransaction;
# import typedb.client.api.concept.thing.Relation;
#
# import javax.annotation.CheckReturnValue;
# import javax.annotation.Nullable;
# import java.util.stream.Stream;
#
# public interface RelationType extends ThingType {
#
#     @Override
#     default boolean isRelationType() {
#         return true;
#     }
#
#     @Override
#     RelationType.Remote asRemote(TypeDBTransaction transaction);
#
#     interface Remote extends ThingType.Remote, RelationType {
#
#         @CheckReturnValue
#         Relation create();
#
#         @Override
#         @CheckReturnValue
#         Stream<? extends Relation> getInstances();
#
#         @CheckReturnValue
#         Stream<? extends RoleType> getRelates();
#
#         @Nullable
#         @CheckReturnValue
#         RoleType getRelates(String roleLabel);
#
#         void setRelates(String roleLabel);
#
#         void setRelates(String roleLabel, String overriddenLabel);
#
#         void unsetRelates(String roleLabel);
#
#         @Override
#         @CheckReturnValue
#         Stream<? extends RelationType> getSubtypes();
#
#         void setSupertype(RelationType superRelationType);
#     }
# }


text:
using TimeZones: timezone_names
using Dates: Millisecond

names = timezone_names()
@test length(names) >= 436
@test isa(names, Vector{String})
@test issorted(names)

# Make sure that extra time zones exist from "deps/tzdata/utc".
# If tests fail try rebuilding the package: Pkg.build("TimeZones")
@test "UTC" in names
@test "GMT" in names

timezones = all_timezones()
@test length(timezones) == length(names)
@test isa(timezones, Vector{TimeZone})

timezones = timezones_from_abbr("EET")
@test length(timezones) >= 29
@test isa(timezones, Vector{TimeZone})

# Note: Unlike time zone names it is possible, although unlikely, for the number of
# abbreviations to decrease over time. A number of abbreviations were eliminated in 2016
# instead of using invented or obsolete alphanumeric time zone abbreviations. Since the
# number of abbreviates is hard to predict we'll avoid testing the number of abbreviations.
abbrs = timezone_abbrs()
@test isa(abbrs, Vector{String})
@test issorted(abbrs)

wpg = first(compile("America/Winnipeg", tzdata["northamerica"]))
apia = first(compile("Pacific/Apia", tzdata["australasia"]))
paris = first(compile("Europe/Paris", tzdata["europe"]))

@testset "next_transition_instant" begin
    @testset "non-existent" begin
        local zone = FixedTimeZone("CST", -6 * 3600)

        instant = next_transition_instant(ZonedDateTime(2018, 1, 1, wpg))
        expected_instant = ZonedDateTime(DateTime(2018, 3, 11, 8), wpg, zone)
        expected_valid = ZonedDateTime(2018, 3, 11, 3, wpg)

        @test instant === expected_instant
        @test instant == expected_valid
        @test instant !== expected_valid
        @test instant + Millisecond(0) === expected_valid
    end

    @testset "ambiguous" begin
        local zone = FixedTimeZone("CDT", -6 * 3600, 3600)

        instant = next_transition_instant(ZonedDateTime(2018, 6, 1, wpg))
        expected_instant = ZonedDateTime(DateTime(2018, 11, 4, 7), wpg, zone)
        expected_valid = ZonedDateTime(2018, 11, 4, 1, wpg, 2)

        @test instant === expected_instant
        @test instant == expected_valid
        @test instant !== expected_valid
        @test instant + Millisecond(0) === expected_valid
    end

    @testset "upcoming" begin
        local patch = @patch now(tz::TimeZone) = ZonedDateTime(2000, 1, 1, tz)
        apply(patch) do
            @test next_transition_instant(wpg) == ZonedDateTime(2000, 4, 2, 3, wpg)
        end
    end

    @testset "fixed time zone" begin
        @test next_transition_instant(ZonedDateTime(2018, 1, 1, tz"UTC")) === nothing
    end

    @testset "no future transition" begin
        # Determine the last transition instant for the time zone
        t = wpg.transitions[end]
        last_trans_instant = ZonedDateTime(t.utc_datetime, wpg, t.zone)

        @test next_transition_instant(last_trans_instant) !== nothing
        @test next_transition_instant(last_trans_instant + Millisecond(1)) === nothing
    end
end

@testset "show_next_transition" begin
    @testset "non-existent" begin
        @test sprint(show_next_transition, ZonedDateTime(2018, 1, 1, wpg)) ==
            """
            Transition Date:   2018-03-11
            Local Time Change: 02:00 → 03:00 (Forward)
            Offset Change:     UTC-6/+0 → UTC-6/+1
            Transition From:   2018-03-11T01:59:59.999-06:00 (CST)
            Transition To:     2018-03-11T03:00:00.000-05:00 (CDT)
            """
    end

    @testset "ambiguous" begin
         @test sprint(show_next_transition, ZonedDateTime(2018, 6, 1, wpg)) ==
            """
            Transition Date:   2018-11-04
            Local Time Change: 02:00 → 01:00 (Backward)
            Offset Change:     UTC-6/+1 → UTC-6/+0
            Transition From:   2018-11-04T01:59:59.999-05:00 (CDT)
            Transition To:     2018-11-04T01:00:00.000-06:00 (CST)
            """
    end

    @testset "standard offset change" begin
        @test sprint(show_next_transition, ZonedDateTime(2011, 12, 1, apia)) ==
            """
            Transition Date:   2011-12-30
            Local Time Change: 00:00 → 00:00 (Forward)
            Offset Change:     UTC-11/+1 → UTC+13/+1
            Transition From:   2011-12-29T23:59:59.999-10:00 (SDT)
            Transition To:     2011-12-31T00:00:00.000+14:00 (WSDT)
            """
    end

    @testset "dst offset change" begin
        @test sprint(show_next_transition, ZonedDateTime(1945, 4, 1, paris)) ==
            """
            Transition Date:   1945-04-02
            Local Time Change: 02:00 → 03:00 (Forward)
            Offset Change:     UTC+0/+1 → UTC+0/+2
            Transition From:   1945-04-02T01:59:59.999+01:00 (WEST)
            Transition To:     1945-04-02T03:00:00.000+02:00 (WEMT)
            """
    end

    @testset "upcoming" begin
        local patch = @patch now(tz::TimeZone) = ZonedDateTime(2000, 1, 1, tz)
        apply(patch) do
            @test occursin("2000-04-02", sprint(show_next_transition, wpg))
        end
    end

    @testset "fixed time zone" begin
        msg = "No transitions exist in time zone UTC"
        instant = ZonedDateTime(2019 ,1, 1, tz"UTC")

        io = IOBuffer()
        @test_logs (:warn, msg) show_next_transition(io, instant)
        @test isempty(read(seekstart(io), String))
    end

    @testset "no future transition" begin
        msg = "No transition exists in America/Winnipeg after: 2038-03-14T01:59:59.999-06:00"
        instant = ZonedDateTime(wpg.cutoff - Millisecond(1), wpg; from_utc=true)

        io = IOBuffer()
        @test_logs (:warn, msg) show_next_transition(io, instant)
        @test isempty(read(seekstart(io), String))
    end
end


text:
using Test
using SimpleTensorNetworks
using Compose, Viznet

@testset "viz" begin
    tn = TensorNetwork([LabeledTensor(randn(2,2), ['a', 'b']), LabeledTensor(randn(2,2), ['b', 'c']), LabeledTensor(randn(2,2), ['a', 'c'])])
    @test viz_tnet(tn) isa Context
end


text:
# Stores data in a simple dict in memory
struct DictStore <: AbstractStore
  name::String
  a::Dict{String,Vector{UInt8}}
  subdirs::Dict{String,DictStore}
end
function DictStore(name)
  isempty(name) && (name="data")
  a = Dict{String,Vector{UInt8}}()
  subdirs = Dict{String,DictStore}()
  DictStore(name, a, subdirs)
end
DictStore() = DictStore("")
Base.show(io::IO,d::DictStore) = print(io,"Dictionary Storage")

storagesize(d::DictStore) = sum(i->i[1] ∉ (".zattrs",".zarray") ? sizeof(i[2]) : zero(sizeof(i[2])),d.a)
zname(s::DictStore) = s.name

Base.getindex(d::DictStore,i::AbstractString) = get(d.a,i,nothing)
Base.setindex!(d::DictStore,v,i::AbstractString) = d.a[i] = v

Base.delete!(d::DictStore,i::AbstractString) = delete!(d.a,i)

subdirs(d::DictStore) = keys(d.subdirs)
Base.keys(d::DictStore) = keys(d.a)
newsub(d::DictStore, n) = d.subdirs[n] = DictStore(n)
getsub(d::DictStore, n) = d.subdirs[n]

path(d::DictStore) = ""


text:
"""
    This provides units and creates a non-dimensionalization object
"""
module Units
using Unitful
import Unitful: superscript
using Parameters

import Base.show
using GeoParams: AbstractMaterialParam, AbstractMaterialParamsStruct, AbstractPhaseDiagramsStruct, PerpleX_LaMEM_Diagram

# Define additional units that are useful in geodynamics
@unit    Myrs  "Myrs"   MillionYears    1000000u"yr"    false

function __init__()
    Unitful.register(Units)
end

Unitful.register(Units)

# Define a number of useful units
const km    = u"km"
const m     = u"m"
const cm    = u"cm"
const mm    = u"mm"
const Myrs  = u"Myrs"
const yr    = u"yr"
const s     = u"s"
const kg    = u"kg"
const g     = u"g"
const Pa    = u"Pa"
const MPa   = u"MPa"
const kbar  = u"kbar"
const Pas   = u"Pa*s"
const K     = u"K"
const C     = u"°C"
const mol   = u"mol"
const kJ    = u"kJ"
const J     = u"J"
const Watt  = u"W"
const μW    = u"μW"

export
    km, m, cm, mm, Myrs, yr, s, MPa, Pa, kbar, Pas, K, C, g, kg, mol, J, kJ, Watt, μW,
    GeoUnit, GeoUnits, GEO_units, SI_units, NO_units, AbstractGeoUnits,
    Nondimensionalize, Nondimensionalize!, Dimensionalize, Dimensionalize!,
    superscript, upreferred, GEO, SI, NONE, isDimensional, Value, NumValue, Unit

"""
AbstractGeoUnits

Abstract supertype for geo units.
"""
abstract type AbstractGeoUnits{TYPE} end

abstract type AbstractUnitType end

struct GEO <: AbstractUnitType end
struct SI  <: AbstractUnitType end
struct NONE<: AbstractUnitType end

"""
    Structure that holds a GeoUnit parameter and their dimensions

    Having that is useful, as non-dimensionalization removes the units from a number
    and we thus no longer know how to transfer it back to the correct units.

"""
mutable struct GeoUnit
    val                 # the actual value ()
    unit :: Unitful.FreeUnits
end

GeoUnit(v::Unitful.Quantity)            =   GeoUnit(v, unit(v))     # store the units
GeoUnit(v::Number)                      =   GeoUnit(v, NoUnits)     # in case we just have a number with no units
GeoUnit(v::Array)                       =   GeoUnit(v, NoUnits)     # array, no units
GeoUnit(v::Array{Unitful.Quantity})     =   GeoUnit(v, unit.(v))    # with units
GeoUnit(v::StepRange)                   =   GeoUnit(v, unit.(v))    # with units
Value(v::GeoUnit)                       =   v.val                   # get value of GeoUnit
NumValue(v::GeoUnit)                    =   ustrip(v.val)           # numeric value, with no units
Unit(v::GeoUnit)                        =   v.unit                  # extract unit

Base.convert(::Type{Float64}, v::GeoUnit)       =   v.val
Base.convert(::Type{GeoUnit}, v::Quantity)      =   GeoUnit(v)
Base.convert(::Type{GeoUnit}, v::Number)        =   GeoUnit(v, NoUnits)
Base.convert(::Type{GeoUnit}, v::Vector)        =   GeoUnit(v, unit(v[1]))
Base.convert(::Type{GeoUnit}, v::Array)         =   GeoUnit(v, unit(v[1]))
Base.convert(::Type{GeoUnit}, v::StepRangeLen)  =   GeoUnit(v, unit(v[1]))
Base.convert(::Type{GeoUnit}, v::StepRange)     =   GeoUnit(v, unit(v[1]))

# define a few basic routines so we can easily operate with GeoUnits
Base.show(io::IO, x::GeoUnit)   =   println(x.val)
Base.length(v::GeoUnit)         =   length(v.val)
Base.size(v::GeoUnit)           =   size(v.val)

Base.:*(x::GeoUnit, y::Number)  = x.val*y
Base.:+(x::GeoUnit, y::Number)  = x.val+y
Base.:/(x::GeoUnit, y::Number)  = x.val/y
Base.:-(x::GeoUnit, y::Number)  = x.val-y

Base.:*(x::GeoUnit, y::GeoUnit) = x.val*y.val
Base.:+(x::GeoUnit, y::GeoUnit) = x.val+y.val
Base.:/(x::GeoUnit, y::GeoUnit) = x.val/y.val
Base.:-(x::GeoUnit, y::GeoUnit) = x.val-y.val

Base.:*(x::Number, y::GeoUnit)  = y.val*x
Base.:+(x::Number, y::GeoUnit)  = y.val+x
Base.:/(x::Number, y::GeoUnit)  = x/y.val
Base.:-(x::Number, y::GeoUnit)  = x-y.val

Base.:*(x::GeoUnit, y::Quantity)  = GeoUnit(x.val*y, x.unit)
Base.:+(x::GeoUnit, y::Quantity)  = GeoUnit(x.val+y, x.unit)
Base.:/(x::GeoUnit, y::Quantity)  = GeoUnit(x.val/y, x.unit)
Base.:-(x::GeoUnit, y::Quantity)  = GeoUnit(x.val-y, x.unit)

Base.:*(x::GeoUnit, y::Array)   = GeoUnit(x.val*y, x.unit)
Base.:/(x::GeoUnit, y::Array)   = GeoUnit(x.val/y, x.unit)
Base.:+(x::GeoUnit, y::Array)   = GeoUnit(x.val+y, x.unit)
Base.:-(x::GeoUnit, y::Array)   = GeoUnit(x.val-y, x.unit)

Base.:*(x::Array, y::GeoUnit)   = GeoUnit(x*y.val, y.unit)
Base.:/(x::Array, y::GeoUnit)   = GeoUnit(x/y.val, y.unit)
Base.:+(x::Array, y::GeoUnit)   = GeoUnit(x +y.val, y.unit)
Base.:-(x::Array, y::GeoUnit)   = GeoUnit(x -y.val, y.unit)

Base.:*(x::GeoUnit, y::Vector)   = GeoUnit(x.val*y, x.unit)
Base.:/(x::GeoUnit, y::Vector)   = GeoUnit(x.val/y, x.unit)
Base.:+(x::GeoUnit, y::Vector)   = GeoUnit(x.val+y, x.unit)
Base.:-(x::GeoUnit, y::Vector)   = GeoUnit(x.val-y, x.unit)

Base.:*(x::GeoUnit, y::StepRange)   = GeoUnit(x.val*y, x.unit)
Base.:/(x::GeoUnit, y::StepRange)   = GeoUnit(x.val/y, x.unit)
Base.:+(x::GeoUnit, y::StepRange)   = GeoUnit(x.val+y, x.unit)
Base.:-(x::GeoUnit, y::StepRange)   = GeoUnit(x.val-y, x.unit)

Base.getindex(x::GeoUnit, i::Int64, j::Int64, k::Int64) = x.val[i,j,k]
Base.getindex(x::GeoUnit, i::Int64, j::Int64) = x.val[i,j]
Base.getindex(x::GeoUnit, i::Int64) = x.val[i]

Base.setindex!(x::GeoUnit, v::Any, i::Int64, j::Int64, k::Int64) = x.val[i,j,k] = v
Base.setindex!(x::GeoUnit, v::Any, i::Int64, j::Int64) = x.val[i,j] = v
Base.setindex!(x::GeoUnit, v::Any, i::Int64) = x.val[i] = v

"""
    GeoUnits

    Structure that holds parameters used for non-dimensionalization
"""
@with_kw_noshow struct GeoUnits{TYPE}
    # Selectable input parameters
    temperature     =   1               #   Characteristic temperature  [C or K]
    length          =   1               #   Characteristic length unit  [km, m or -]
    stress          =   1               #   Characteristic stress unit  [MPa, Pa or -]
    time            =   1               #   Characteristic time unit    [Myrs, s pr -]
    viscosity       =   1

    # primary characteristic units
    K               =   1;                      # temperature in SI units for material parameter scaling
    s               =   1;                      # time in SI units for material parameter scaling
    m               =   1;                      # length in SI units for material parameter scaling
    Pa              =   1;                      # stress in SI units
    kg              =   upreferred(Pa*m*s^2)    # compute mass from pascal. Note: this may result in very large values

    # Main SI units (used later for nondimensionalization)
    Length          =   m
    Mass            =   kg
    Time            =   s
    Temperature     =   K
    Amount          =   1mol
    Second          =   s

    # Not defined, as they are not common in geodynamics:
    #Current
    #Luminosity

    # Derived units
    N               =   kg*m/s^2        # Newton
    J               =   N*m             # Joule
    W               =   J/s             # Watt
    area            =   m^2             # area   in SI units for material parameter scaling
    volume          =   m^3             # volume
    velocity        =   m/s
    density         =   kg/m^3
    acceleration    =   m/s^2
    force           =   kg*m/s^2
    strainrate      =   1/s
    heatcapacity    =   J/kg/K
    conductivity    =   W/m/K

    # Helpful
    SecYear         =   3600*24*365.25
    Myrs            =   1e6
    cmYear          =   SecYear*100     # to transfer m/s -> cm/yr

end

"""
    GEO_units(;length=1000km, temperature=1000C, stress=10MPa, viscosity=1e20Pas)

Creates a non-dimensionalization object using GEO units.

GEO units implies that upon dimensionalization, `time` will be in `Myrs`, `length` in `km`, stress in `MPa`, etc.
which is more convenient for typical geodynamic simulations than SI units
The characteristic values given as input can be in arbitrary units (`km` or `m`), provided the unit is specified.

# Examples:
```julia-repl
julia> CharUnits = GEO_units()
Employing GEO units
Characteristic values:
         length:      1000 km
         time:        0.3169 Myrs
         stress:      10 MPa
         temperature: 1000.0 °C
julia> CharUnits.velocity
1.0e-7 m s⁻¹
```
If we instead have a crustal-scale simulation, it is likely more appropriate to use a different characteristic `length`:
```julia-repl
julia> CharUnits = GEO_units(length=10km)
Employing GEO units
Characteristic values:
         length:      10 km
         time:        0.3169 Myrs
         stress:      10 MPa
         temperature: 1000.0 °C
```
"""
function GEO_units(;length=1000km, temperature=1000C, stress=10MPa, viscosity=1e20Pas)

    if unit(temperature)==NoUnits;  temperature = temperature*C;        end
    if unit(length)==NoUnits;       length      = length*u"km";         end
    if unit(stress)==NoUnits;       stress      = stress*u"MPa";        end
    if unit(viscosity)==NoUnits;    viscosity   = viscosity*u"Pa*s";    end

    T       =   uconvert(C,      temperature)
    Le      =   uconvert(km,     length);
    Sigma   =   uconvert(MPa,    stress)
    Eta     =   uconvert(Pas,    viscosity)

    T_SI    =   uconvert(K,      T);
    Le_SI   =   uconvert(m,      Le);
    Sigma_SI=   uconvert(Pa,     Sigma)
    Time_SI =   Eta/Sigma_SI;
    t       =   uconvert(Myrs,   Time_SI)

    GeoUnits{GEO}(length=Le, temperature=T,      stress=Sigma,  viscosity=Eta, time=t,
                  m=Le_SI,   K=T_SI,   Pa=Sigma_SI,  s=Time_SI)
end

"""
    CharUnits = SI_units(length=1000m, temperature=1000K, stress=10Pa, viscosity=1e20)

Specify the characteristic values using SI units

# Examples:
```julia-repl
julia> CharUnits = SI_units(length=1000m)
Employing SI units
Characteristic values:
         length:      1000 m
         time:        1.0e19 s
         stress:      10 Pa
         temperature: 1000.0 K
```
Note that the same can be achieved if the input is given in `km`:
```julia-repl
julia> CharUnits = SI_units(length=1km)
```
"""
function SI_units(;length=1000m, temperature=1000K, stress=10Pa, viscosity=1e20Pas)

    if unit(temperature)==NoUnits;  temperature = temperature*K;     end
    if unit(length)==NoUnits;       length      = length*u"m";         end
    if unit(stress)==NoUnits;       stress      = stress*u"Pa";        end
    if unit(viscosity)==NoUnits;    viscosity   = viscosity*u"Pa*s";    end

    T       =   uconvert(K,     temperature)
    Le      =   uconvert(m,     length);
    Sigma   =   uconvert(Pa,    stress)
    Eta     =   uconvert(Pas,   viscosity)

    T_SI    =   uconvert(K,      T);
    Le_SI   =   uconvert(m,      Le);
    Sigma_SI=   uconvert(Pa,     Sigma)
    Time_SI =   Eta/Sigma_SI;
    t       =   uconvert(s,      Time_SI)

    GeoUnits{SI}(length=Le, temperature=T,      stress=Sigma,  viscosity=Eta, time=t,
                 m=Le_SI,   K=T_SI,   Pa=Sigma_SI,  s=Time_SI)
end

"""
    CharUnits = NO_units(length=1, temperature=1, stress=1, viscosity=1)

Specify the characteristic values in non-dimensional units

# Examples:
```julia-repl
julia> using GeoParams;
julia> CharUnits = NO_units()
Employing NONE units
Characteristic values:
         length:      1
         time:        1.0
         stress:      1
         temperature: 1.0
```
"""
function NO_units(;length=1, temperature=1, stress=1, viscosity=1)

    if unit(temperature)!=NoUnits;  error("temperature should not have units")    end
    if unit(length)!=NoUnits;       error("length should not have units")    end
    if unit(stress)!=NoUnits;       error("stress should not have units")    end
    if unit(viscosity)!=NoUnits;    error("viscosity should not have units")    end

    T       =   temperature
    Le      =   length;
    Sigma   =   stress
    Eta     =   viscosity
    Time    =   Eta/Sigma;

    GeoUnits{NONE}(length=Le, temperature=T, stress=Sigma, viscosity=Eta, time=Time,
                 m=Le, K=T, Pa=Sigma, s=Time)
end

"""
    Nondimensionalize(param, CharUnits::GeoUnits{TYPE})

Nondimensionalizes `param` using the characteristic values specified in `CharUnits`

# Example 1
```julia-repl
julia> using GeoParams;
julia> CharUnits =   GEO_units();
julia> v         =   3cm/yr
3 cm yr⁻¹
julia> v_ND      =   Nondimensionalize(v, CharUnits)
0.009506426344208684
```
# Example 2
In geodynamics one sometimes encounters more funky units
```julia-repl
julia> CharUnits =   GEO_units();
julia> A         =   6.3e-2MPa^-3.05*s^-1
0.063 MPa⁻³·⁰⁵ s⁻¹
julia> A_ND      =   Nondimensionalize(A, CharUnits)
7.068716262102384e14
```

In case you are interested to see how the units of `A` look like in different units, use this function from the [Unitful](https://github.com/PainterQubits/Unitful.jl) package:
```julia-repl
julia> uconvert(u"Pa^-3.05*s^-1",A)
3.157479571851836e-20 Pa⁻³·⁰⁵
```
and to see it decomposed in the basic `SI` units of length, mass and time:
```julia-repl
julia> upreferred(A)
3.1574795718518295e-20 m³·⁰⁵ s⁵·¹ kg⁻³·⁰⁵
```
"""
function Nondimensionalize(param, g::GeoUnits{TYPE}) where {TYPE}
    if typeof(param) == String
        param_ND = param # The parameter is a string, cannot be nondimensionalized
    elseif unit(param)!=NoUnits
        dim         =   Unitful.dimension(param);                   # Basic SI units
        char_val    =   1.0;
        foreach((typeof(dim).parameters[1])) do y
            val = upreferred(getproperty(g, Unitful.name(y)))       # Retrieve the characteristic value from structure g
            pow = Float64(y.power)                                  # power by which it should be multiplied
            char_val *= val^pow                                     # multiply characteristic value
        end
        param_ND = upreferred.(param)/char_val
    else
        param_ND = param # The parameter has no units, so there is no way to determine how to nondimensionize it
    end
    return param_ND
end

function Nondimensionalize(param::Array, g::GeoUnits{TYPE}) where {TYPE}
    if unit(param[1])!=NoUnits
        dim         =   Unitful.dimension.(param);                   # Basic SI units
        char_val    =   1.0;
        foreach((typeof(dim[1]).parameters[1])) do y
            val = upreferred(getproperty(g, Unitful.name(y)))       # Retrieve the characteristic value from structure g
            pow = Float64(y.power)                                  # power by which it should be multiplied
            char_val *= val^pow                                     # multiply characteristic value
        end
        param_ND = upreferred.(param)/char_val
    else
        param_ND = param # The parameter has no units, so there is no way to determine how to nondimensionize it
    end
    return param_ND
end

"""
    Nondimensionalize!(param::GeoUnit, CharUnits::GeoUnits{TYPE})

Nondimensionalizes `param` (given as GeoUnit) using the characteristic values specified in `CharUnits` in-place

# Example 1
```julia-repl
julia> using GeoParams;
julia> CharUnits =   GEO_units();
julia> v         =   GeoUnit(3cm/yr)
3 cm yr⁻¹
julia> Nondimensionalize!(v, CharUnits)
0.009506426344208684
```
# Example 2
```julia-repl
julia> CharUnits =   GEO_units();
julia> A         =   GeoUnit(6.3e-2MPa^-3.05*s^-1)
0.063 MPa⁻³·⁰⁵ s⁻¹
julia> A_ND      =   Nondimensionalize(A, CharUnits)
7.068716262102384e14
```
"""
function Nondimensionalize!(param::GeoUnit, g::GeoUnits{TYPE}) where {TYPE}
    if unit.(param.val)!=NoUnits
        dim         =   Unitful.dimension(param.val[1]);                   # Basic SI units
        char_val    =   1.0;
        foreach((typeof(dim).parameters[1])) do y
            val = upreferred(getproperty(g, Unitful.name(y)))       # Retrieve the characteristic value from structure g
            pow = Float64(y.power)                                  # power by which it should be multiplied
            char_val *= val^pow                                     # multiply characteristic value
        end
        param.val = upreferred.(param.val)/char_val;
    else
        param = param # The parameter has no units, so there is no way to determine how to nondimensionize it
    end
end

function Nondimensionalize!(param::String, g::GeoUnits{TYPE}) where {TYPE}
    param_ND = param
    return nothing
end

function Nondimensionalize!(param::Array, g::GeoUnits{TYPE}) where {TYPE}

    if unit.(param[1])!=NoUnits
        dim         =   Unitful.dimension(param[1]);                   # Basic SI units
        char_val    =   1.0;
        foreach((typeof(dim).parameters[1])) do y
            val = upreferred(getproperty(g, Unitful.name(y)))       # Retrieve the characteristic value from structure g
            pow = Float64(y.power)                                  # power by which it should be multiplied
            char_val *= val^pow                                     # multiply characteristic value
        end
        param = upreferred.(param)/char_val;
    else
        param = param # The parameter has no units, so there is no way to determine how to nondimensionize it
    end

end

"""
    Nondimensionalize!(MatParam::AbstractMaterialParam, CharUnits::GeoUnits{TYPE})

Non-dimensionalizes a material parameter structure (e.g., Density, CreepLaw)

"""
function Nondimensionalize!(MatParam::AbstractMaterialParam, g::GeoUnits{TYPE}) where {TYPE}
    for param in fieldnames(typeof(MatParam))
        if typeof(getfield(MatParam, param))==GeoUnit
            z=getfield(MatParam, param)
            Nondimensionalize!(z, g)
            setfield!(MatParam, param, z)
        end
    end
end

"""
    Nondimensionalize!(phase_mat::MaterialParams, g::GeoUnits{TYPE})

Nondimensionalizes all fields within the Material Parameters structure that contain material parameters
"""
function Nondimensionalize!(phase_mat::AbstractMaterialParamsStruct, g::GeoUnits{TYPE}) where {TYPE}

    for param in fieldnames(typeof(phase_mat))
        fld = getfield(phase_mat, param)
        if ~isnothing(fld)
            if typeof(fld[1]) <: AbstractPhaseDiagramsStruct

                # in case we employ a phase diagram
                temp = PerpleX_LaMEM_Diagram(fld[1].Name, CharDim = g)

                setfield!(phase_mat, param, (temp,))
            else
                # otherwise non-dimensionalize
                for i=1:length(fld)
                    if typeof(fld[i]) <: AbstractMaterialParam
                        Units.Nondimensionalize!(fld[i],g)
                    end
                end
            end
        end
    end
    phase_mat.Nondimensional = true
end

"""
    Dimensionalize(param, param_dim::Unitful.FreeUnits, CharUnits::GeoUnits{TYPE})

Dimensionalizes `param` into the dimensions `param_dim` using the characteristic values specified in `CharUnits`.

# Example
```julia-repl
julia> CharUnits =   GEO_units();
julia> v_ND      =   Nondimensionalize(3cm/yr, CharUnits)
0.031688087814028945
julia> v_dim     =   Dimensionalize(v_ND, cm/yr, CharUnits)
3.0 cm yr⁻¹
```

"""
function Dimensionalize(param_ND, param_dim::Unitful.FreeUnits, g::GeoUnits{TYPE}) where {TYPE}

    dim         =   Unitful.dimension(param_dim);                   # Basic SI units
    char_val    =   1.0;
    foreach((typeof(dim).parameters[1])) do y
        val = upreferred(getproperty(g, Unitful.name(y)))       # Retrieve the characteristic value from structure g
        pow = Float64(y.power)                                  # power by which it should be multiplied
        char_val *= val^pow                                     # multiply characteristic value
    end
    param = uconvert.(param_dim, param_ND*char_val)

    return param
end

"""
    Dimensionalize!(param::GeoUnit, CharUnits::GeoUnits{TYPE})

Dimensionalizes `param` again to the values that it used to have using the characteristic values specified in `CharUnits`.

# Example
```julia-repl
julia> CharUnits =   GEO_units();
julia> x = GeoUnit(3cm/yr)
julia> Nondimensionalize!(x, CharUnits)
julia> Dimensionalize!(x, CharUnits)
3.0 cm yr⁻¹
```

"""
function Dimensionalize!(param::GeoUnit, g::GeoUnits{TYPE}) where {TYPE}

    dim         =   Unitful.dimension(param.unit);                   # Basic SI units
    char_val    =   1.0;
    foreach((typeof(dim).parameters[1])) do y
        val = upreferred(getproperty(g, Unitful.name(y)))       # Retrieve the characteristic value from structure g
        pow = Float64(y.power)                                  # power by which it should be multiplied
        char_val *= val^pow                                     # multiply characteristic value
    end
    param.val = uconvert.(param.unit, param.val*char_val)

end

"""
    Dimensionalize!(MatParam::AbstractMaterialParam, CharUnits::GeoUnits{TYPE})

Dimensionalizes a material parameter structure (e.g., Density, CreepLaw)

"""
function Dimensionalize!(MatParam::AbstractMaterialParam, g::GeoUnits{TYPE}) where {TYPE}

    for param in fieldnames(typeof(MatParam))
        if typeof(getfield(MatParam, param))==GeoUnit
            z=getfield(MatParam, param)
            Dimensionalize!(z, g)
            setfield!(MatParam, param, z)
        end
    end

end

"""
    isDimensional(MatParam::AbstractMaterialParam)

`true` if MatParam is in dimensional units.
"""
function isDimensional(MatParam::AbstractMaterialParam)
    isDim = false;
    for param in fieldnames(typeof(MatParam))
        if typeof(getfield(MatParam, param))==GeoUnit
            z=getfield(MatParam, param)
            if unit(z.val)!=NoUnits
                isDim=true;
            end
        end
    end
    return isDim
end

"""
    Dimensionalize!(phase_mat::MaterialParams, g::GeoUnits{TYPE})

Dimensionalizes all fields within the Material Parameters structure that contain material parameters
"""
function Dimensionalize!(phase_mat::AbstractMaterialParamsStruct, g::GeoUnits{TYPE}) where {TYPE}

    for param in fieldnames(typeof(phase_mat))
        fld = getfield(phase_mat, param)
        if ~isnothing(fld)
            for i=1:length(fld)
                if typeof(fld[i]) <: AbstractMaterialParam
                    Units.Dimensionalize!(fld[i],g)
                end
            end
        end
    end
    phase_mat.Nondimensional = false

end

# Define a view for the GEO_Units structure
function show(io::IO, g::GeoUnits{TYPE})  where {TYPE}
    print(io, "Employing $TYPE units \n",
              "Characteristic values: \n",
              "         length:      $(g.length)\n",
              "         time:        $(round(ustrip(g.time),digits=4)) $(unit(g.time))\n",
              "         stress:      $(g.stress)\n",
              "         temperature: $(Float64(g.temperature))\n")
end

# This replaces the viewer of the Unitful package, such that the printing of units is done as floats (better)
function Unitful.superscript(i::Rational{Int64}; io=nothing)
    string(superscript(float(i)))
    if io === nothing
        iocontext_value = nothing
    else
        iocontext_value = get(io, :fancy_exponent, nothing)
    end
    if iocontext_value isa Bool
        fancy_exponent = iocontext_value
    else
        v = get(ENV, "UNITFUL_FANCY_EXPONENTS", Sys.isapple() ? "true" : "false")
        t = tryparse(Bool, lowercase(v))
        fancy_exponent = (t === nothing) ? false : t
    end
    if fancy_exponent
        return superscript(float(i))
    else
        return  "^" * string(float(i))
    end

end

Unitful.superscript(i::Float64) = map(repr(i)) do c
    c == '-' ? '\u207b' :
    c == '1' ? '\u00b9' :
    c == '2' ? '\u00b2' :
    c == '3' ? '\u00b3' :
    c == '4' ? '\u2074' :
    c == '5' ? '\u2075' :
    c == '6' ? '\u2076' :
    c == '7' ? '\u2077' :
    c == '8' ? '\u2078' :
    c == '9' ? '\u2079' :
    c == '0' ? '\u2070' :
    c == '.' ? '\u0387' :
    error("unexpected character")
end

end


text:
import ..UncertainValues:
	TheoreticalDistributionScalarValue,
    AbstractUncertainTwoParameterScalarValue,
    AbstractUncertainThreeParameterScalarValue,
    UncertainScalarTheoreticalTwoParameter,
    UncertainScalarTheoreticalThreeParameter,
    UncertainScalarNormallyDistributed,
    UncertainScalarUniformlyDistributed,
    UncertainScalarBetaDistributed,
    UncertainScalarBetaPrimeDistributed,
    UncertainScalarBetaBinomialDistributed,
    UncertainScalarBinomialDistributed,
    UncertainScalarGammaDistributed,
    UncertainScalarFrechetDistributed
import Distributions

# Resample for generic
resample(uv::TheoreticalDistributionScalarValue) = rand(uv.distribution)
resample(uv::TheoreticalDistributionScalarValue, n::Int) = rand(uv.distribution, n)

# Custom resample methods for each type of uncertain scalars based on
# distributions (in case we want to implement custom sampling for some of them)
resample(uv::UncertainScalarTheoreticalThreeParameter) = rand(uv.distribution)
resample(uv::UncertainScalarTheoreticalTwoParameter) = rand(uv.distribution)
resample(uv::UncertainScalarNormallyDistributed) = rand(uv.distribution)
resample(uv::UncertainScalarUniformlyDistributed) = rand(uv.distribution)
resample(uv::UncertainScalarBetaDistributed) = rand(uv.distribution)
resample(uv::UncertainScalarBetaPrimeDistributed) = rand(uv.distribution)
resample(uv::UncertainScalarBetaBinomialDistributed) = rand(uv.distribution)
resample(uv::UncertainScalarGammaDistributed) = rand(uv.distribution)
resample(uv::UncertainScalarFrechetDistributed) = rand(uv.distribution)
resample(uv::UncertainScalarBinomialDistributed) = rand(uv.distribution)

resample(uv::UncertainScalarTheoreticalThreeParameter, n::Int) = rand(uv.distribution, n)
resample(uv::UncertainScalarTheoreticalTwoParameter, n::Int) = rand(uv.distribution, n)
resample(uv::UncertainScalarNormallyDistributed, n::Int) = rand(uv.distribution, n)
resample(uv::UncertainScalarUniformlyDistributed, n::Int) = rand(uv.distribution, n)
resample(uv::UncertainScalarBetaDistributed, n::Int) = rand(uv.distribution, n)
resample(uv::UncertainScalarBetaPrimeDistributed, n::Int) = rand(uv.distribution, n)
resample(uv::UncertainScalarBetaBinomialDistributed, n::Int) = rand(uv.distribution, n)
resample(uv::UncertainScalarGammaDistributed, n::Int) = rand(uv.distribution, n)
resample(uv::UncertainScalarFrechetDistributed, n::Int) = rand(uv.distribution, n)
resample(uv::UncertainScalarBinomialDistributed, n::Int) = rand(uv.distribution, n)

resample(x::Distributions.Truncated) = rand(x)
resample(x::Distributions.Truncated, n::Int) = rand(x, n)

resample(x::Distributions.Distribution) = rand(x)
resample(x::Distributions.Distribution, n::Int) = rand(x, n)


text:
function build_couplings(; polaron=false)
    σ₊ = [0.0 1;0 0]
    σ₋ = [0.0 0;1 0]
    if polaron == false
        ConstantCouplings(["10ZII", "IZI", "IIZ"], unit=:ħ)
    else
        ConstantCouplings([σ₊ ⊗ σi ⊗ σi, σ₋ ⊗ σi ⊗ σi, σi ⊗ σz ⊗ σi, σi ⊗ σi ⊗ σz], unit=:ħ)
    end
end

function interp_couplings(proj, inds, sm, polaron, stage)
    if polaron == false
        ops = [[op[i][1:4, 1:4] for op in proj.op] for i in 1:3]
    else
        spfun, = build_s(stage)
        As = A.(spfun.(proj.s))
        ops1 = [As .* [op[i][1:4, 1:4] for op in proj.op] for i in 1:2]
        ops2 = [[op[i][1:4, 1:4] for op in proj.op] for i in 3:4]
        ops = vcat(ops1, ops2)
    end
    ops = [interp_single_coupling(proj.s, x, inds, sm) for x in ops]
    if polaron == false
        [CustomCouplings(ops, unit=:ħ), CustomCouplings([ops[1]], unit=:ħ)]
    else
        [CustomCouplings(ops[1:2], unit=:ħ), CustomCouplings(ops[3:4], unit=:ħ)]
    end
end

function interp_single_coupling(s, op, ::Nothing, ::Any)
    coupling = function (s, itp)
        res = Array{Float64,2}(undef, 4, 4)
        for i = 1:4
            for j = 1:4
                res[i, j] = itp(i, j, s)
            end
        end
        res
    end
    itp = construct_interpolations(s, op)
    (s) -> coupling(s, itp)
end

function interp_single_coupling(s_axis, op, inds, sm)
    coupling = function (s, itp)
        res = Array{Float64,2}(undef, 4, 4)
        for i = 1:4
            for j = 1:4
                res[i, j] = itp(i, j, s)
            end
        end
        res
    end
    itp1 = construct_interpolations(s_axis[1:inds], op[1:inds])
    itp2 = construct_interpolations(s_axis[inds + 1:end], op[inds + 1:end])
    function (s)
        if s <= sm
            coupling(s, itp1)
        else
            coupling(s, itp2)
        end
    end
end

function rotate_couplings(coupling, v, polaron)
    cmat = coupling(0)
    cmat = [real.(v' * m * v) for m in cmat]
    if polaron == false
        [ConstantCouplings(cmat, unit=:ħ), ConstantCouplings([cmat[1]], unit=:ħ)]
    else
        [ConstantCouplings(A(sp) * cmat[1:2], unit=:ħ), ConstantCouplings(cmat[3:4], unit=:ħ)]
    end
end


text:
#=
  Random walk model

  from Foundations of Probabilistic Programming,
  page 3

  The distribution of u and v are about the same as this
ulia> show_var_dist_pct(rand(DiscreteUniform(-100,100),10000),10)

Distributions of variable (num:10)
-76.00000 =>      67  0.006700
-38.00000 =>      65  0.006500
17.00000 =>      64  0.006400
59.00000 =>      64  0.006400
81.00000 =>      63  0.006300
24.00000 =>      62  0.006200
74.00000 =>      61  0.006100
-6.00000 =>      61  0.006100
-100.00000 =>      60  0.006000
-95.00000 =>      60  0.006000

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

data = []
@model function random_walk()
    u ~ DiscreteUniform(-10,10)
    v ~ DiscreteUniform(-10,10)
    function step(tu=0,tv=0)
        x ~ flip()
        y ~ flip()
        tu = tu + (x-y)
        tv = tv + (x+y-1)
        return tu,tv
    end
    # println("\nSTART")
    datad = [[0,0]]
    u,v = step(0,0)
    push!(datad,[u,v])
    c = 0
    while u != 0 || v != 0
        u,v = step(u,v)
        push!(datad,[u,v])
        c += 1
        # Note: We must restrict the walk, otherwise
        #       it'll blow up...
        if u == 0 && v == 0
            println("BACK TO 0,0!")
        end
        if c < 10
            push!(data,datad)
        else
            break
        end

    end
end

model = random_walk()

num_chns = 4
chns = sample(model, MH(), MCMCThreads(), 1000, num_chns)
# chns = sample(model, MH(), 10_000)

display(chns)
show_var_dist_pct(chns, :x)
show_var_dist_pct(chns, :y)
show_var_dist_pct(chns, :u,20)
show_var_dist_pct(chns, :v,20)

# for d in data
#    println(d)
# end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
__precompile__(false)  # Task cannot be serialized

module JobSchedulers

using Base.Threads
using Dates, DataFrames, JSON
using JLD2
using Pipelines

include("jobs.jl")
export Job, result

include("thread_utils.jl")

include("scheduler.jl")
export B, KB, MB, GB, TB
export submit!, cancel!
export QUEUING, RUNNING, DONE, FAILED, CANCELLED

include("pretty_print.jl")
export queue, all_queue, json_queue

include("query.jl")
export job_query_by_id, job_query

include("control.jl")
export scheduler_start, scheduler_stop, scheduler_status
export set_scheduler_update_second,
set_scheduler_max_cpu,
set_scheduler_max_mem,
set_scheduler_max_job

include("backup.jl")
export set_scheduler_backup, backup

include("compat_pipelines.jl")
export close_in_future

function __init__()
    # initiating THREAD_POOL
    c = Channel{Int}(nthreads() - 1)
    THREAD_POOL[] = c
    foreach(i -> put!(c, i), 2:nthreads())  # the thread 1 is reserved for JobScheduler, when nthreads > 2
end

end


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("utilities.jl")

include("../model.jl")

@testset "Indicator bridge" begin
    # linear problem with indicator constraint
    # similar to indicator1_test with reversed z1
    # max  2x1 + 3x2
    # s.t. x1 + x2 <= 10
    #      z1 == 0 ==> x2 <= 8
    #      z2 == 1 ==> x2 + x1/5 <= 9
    #      (1-z1) + z2 >= 1 <=> z2 - z1 >= 0
    model = MOIU.MockOptimizer(Model{Float64}());
    config = MOIT.TestConfig()

    x1 = MOI.add_variable(model)
    x2 = MOI.add_variable(model)
    z1 = MOI.add_variable(model)
    z2 = MOI.add_variable(model)

    vc1 = MOI.add_constraint(model, z1, MOI.ZeroOne())
    @test vc1.value == z1.value
    vc2 = MOI.add_constraint(model, z2, MOI.ZeroOne())
    @test vc2.value == z2.value
    f1 = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z1)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x2)),
        ],
        [0.0, 0.0]
    )
    iset1 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}(MOI.LessThan(8.0))

    bridge = MOIB.bridge_constraint(MOIB.IndicatorActiveOnFalseBridge, model, f1, iset1)

    z1comp = bridge.variable_index
    @test MOI.get(model, MOI.ConstraintFunction(), bridge.zero_one_cons) == MOI.SingleVariable(z1comp)
    @test MOI.get(model, MOI.ConstraintSet(), bridge.disjunction_cons) == MOI.EqualTo(1.0)
    disjunction_cons = MOI.get(model, MOI.ConstraintFunction(), bridge.disjunction_cons)
    for t in disjunction_cons.terms
        @test t.variable_index == z1 || t.variable_index == z1comp
        @test t.coefficient ≈ 1.0
    end
end