text:
#2021/6/23
#Edited by shwh
#This code can be solved Euler Eq.
#conversion from png to gif is conducted by python
using Plots

#parameter
tmax = 300.0
dt = 0.001
xmin, xmax = 0.0, 1.0
xmid = (xmax-xmin)/2
xsp = 100
dx = (xmax-xmin)/xsp
ul, rhol, pl = 0.0, 1.0, 1.0
ur, rhor, pr = 0.0, 0.125, 0.1
gamma = 1.4

xp = collect(0.0:dx:xmax) #for plot
#setup Initial
function setup()
    #global x, Qc, qf, bol, bor
    x = collect(0.0:dx:xmax)
    append!(x, xmax+dx)
    insert!(x, 1, -dx)
    Qc = zeros(length(x),3) #Qc = [rho, rhou, e]
    qf = zeros(length(x),3) #qf = [u,rho, p]
    bol = zeros(3)
    bor = zeros(3)

    for i = 1:length(x)
        if x[i] <= xmid
            qf[i,1] = ul
            qf[i,2] = rhol
            qf[i,3] = pl
            Qc[i,1] = rhol
            Qc[i,2] = rhol*ul
            Qc[i,3] = pl/(gamma-1)+0.5*rhol*ul^2
        else
            qf[i,1] = ur
            qf[i,2] = rhor
            qf[i,3] = pr
            Qc[i,1] = rhor
            Qc[i,2] = rhor*ur
            Qc[i,3] = pr/(gamma-1)+0.5*rhor*ur^2
        end
    end
    for i = 1:3
        bol[i] = Qc[1,i]
        bor[i] = Qc[end,i]
    end
    return x, Qc, qf, bol, bor
end

#calc Ap and Am
function A_pm(ite, Qc, qf)
    #H = (e+p)/rho
    H = (Qc[ite,3]+qf[ite,3])/Qc[ite,1]
    u = qf[ite,1]
    c = sqrt((gamma-1)*(H-0.5*u^2))
    b1 = 0.5*u^2*(gamma-1)/c^2
    b2 = (gamma-1)/c^2

    R = [1.0 1.0 1.0; u-c u u+c; H-u*c 0.5*u^2 H+u*c]
    Rinv = [0.5*(b1+u/c) -0.5*(1/c+b2*u) 0.5*b2; 1-b1 b2*u -b2; 0.5*(b1-u/c) 0.5*(1/c-b2*u) 0.5*b2]

    lam = [(u-c) 0.0 0.0; 0.0 u 0.0; 0.0 0.0 (u+c)]
    lam_abs = [abs(u-c) 0.0 0.0; 0.0 abs(u) 0.0; 0.0 0.0 abs(u+c)]

    return R, Rinv, lam, lam_abs
end

#calc flux by fvs
function fvs(x, Qc, qf)
    #F(i+1/2) = F[i]
    Fp = zeros(length(x),3)
    for i = 1:length(x)-1
        R, Rinv, lam, lam_abs = A_pm(i, Qc, qf)
        Ap =  R*(lam+lam_abs)*Rinv

        R, Rinv, lam, lam_abs = A_pm(i+1, Qc, qf)
        Am = R*(lam-lam_abs)*Rinv

        Fp[i,:] = 0.5*(Ap*Qc[i,:] + Am*Qc[i+1,:])
    end
    return Fp
end

#calc Res = Fp[i]-Fp[i-1]
function calc_Res(Fp,x)
    Res = zeros(length(x), 3)
    for i = 2:length(x)-1
        Res[i,:] = Fp[i,:] - Fp[i-1,:]
    end
    return Res
end

#calc Qc
function calc_Qc(Res,x,Qc)
    lo_R = copy(Res)
    for i = 2:length(x)-1
        Qc[i,:] = Qc[i,:] - dt/dx*lo_R[i,:]
    end
    return Qc
end

function bound(Qc,bol,bor)
    for i = 1:3
        Qc[1,i] = 2*bol[i]-Qc[2,i]
        Qc[end,i] = Qc[end-1,i]
    end
    return Qc
end

function qf2Qc(qf,x)
    lo_Qc = zeros(length(x),3)
    for i = 1:length(x)
        lo_Qc[i,1] = qf[i,2]
        lo_Qc[i,2] = qf[i,2]*qf[i,1]
        lo_Qc[i,3] = (qf[i,3]/(gamma-1)+0.5*qf[i,2]*qf[i,1]^2)
    end
    return lo_Qc
end

function Qc2qf(Qc,x)
    lo_qf = zeros(length(x),3)
    for i = 1:length(x)
        lo_qf[i,1] = Qc[i,2]/Qc[i,1]
        lo_qf[i,2] = Qc[i,1]
        lo_qf[i,3] = (gamma-1)*(Qc[i,3]-0.5*Qc[i,1]*((Qc[i,2]/Qc[i,1])^2))
    end
    return lo_qf
end

function main()
    x, Qc, qf, bol, bor = setup()
    plist = zeros(length(x)-2)
    ulist = zeros(length(x)-2)
    rholist = zeros(length(x)-2)
    for i = 1:length(x)-2
        ulist[i] = qf[i+1,1]
        rholist[i] = qf[i+1,2]
        plist[i] = qf[i+1,3]
    end
    p1 = plot(xp,ulist,xlims = (0.0,1.0), ylims = (0.0,1.0),lw = 3, color = "red", label = "u", title = "0.0s")
    p2 = plot(xp, rholist, xlims = (0.0,1.0), ylims = (0.0,1.0), lw = 3, color = "blue", label = "rho")
    p3 = plot(xp, plist, xlims = (0.0,1.0), ylims = (0.0,1.0), lw = 3, color = "green", label = "p")
    plot(p1, p2, p3, layout=(3,1), size = (1600,900))
    savefig("urhop_"*string(Int(0))*".png")

    for k = 1:tmax
        Fp = fvs(x,Qc,qf)
        Res = calc_Res(Fp,x)
        Qc = calc_Qc(Res,x,Qc)
        Qc = bound(Qc,bol,bor)
        qf = Qc2qf(Qc,x)
        for i = 1:length(x)-2
            ulist[i] = qf[i+1,1]
            rholist[i] = qf[i+1,2]
            plist[i] = qf[i+1,3]
        end
        if k%10 == 0
            p1 = plot(xp,ulist,xlims = (0.0,1.0), ylims = (0.0,1.0),lw = 3, color = "red", label = "u",title = (string(k*dt)*"s"))
            p2 = plot(xp, rholist, xlims = (0.0,1.0), ylims = (0.0,1.0), lw = 3, color = "blue", label = "rho")
            p3 = plot(xp, plist, xlims = (0.0,1.0), ylims = (0.0,1.0), lw = 3, color = "green", label = "p")
            plot(p1, p2, p3, layout=(3,1), size = (1600,900))
            savefig("urhop_"*string(Int(k))*".png")
        end
    end
    return plist
end

@time main()


text:
using OpenTelemetryUber

using Documenter

open(joinpath(@__DIR__, "..", "README.md"), "r") do f_src
    open(joinpath(@__DIR__, "src", "index.md"), "w") do f_dest
        s_dest = read(f_src, String)
        s_dest = replace(s_dest, "<!-- ```@raw html -->" => "```@raw html")
        s_dest = replace(s_dest, "<!-- ``` -->" => "```")
        s_dest = replace(s_dest, "```julia" => "```@example")
        write(f_dest, s_dest)
    end
end

makedocs(
    modules = [OpenTelemetryUber],
    format = Documenter.HTML(
        prettyurls = true,
        analytics = "G-YWC6SQHSVM",
        assets = ["assets/favicon.ico"],
    ),
    sitename = "OpenTelemetry.jl",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "Tutorial" => Any[
            "View Metrics in Prometheus"=>"tutorials/View_Metrics_in_Prometheus/index.md",
            "View Metrics in Prometheus through Open Telemetry Collector"=>"tutorials/View_Metrics_in_Prometheus_through_Open_Telemetry_Collector/index.md",
            "View Metrics/Traces/Logs in Grafana/ElasticAPM"=>"tutorials/View_Logs_Traces,_and_Metrics_Together_in_Grafana_ElasticAPM/index.md",
        ],
        "Tips for Developers" => "tips.md",
        "FAQ" => "FAQ.md",
        "Design" => Any[
            "OpenTelemetryAPI.md",
            "OpenTelemetrySDK.md",
            "OpenTelemetryExporterOtlpProtoGrpc.md",
            "OpenTelemetryExporterPrometheus.md",
            "OpenTelemetryProto.md",
            "OpenTelemetryInstrumentationBase.md",
            "OpenTelemetryInstrumentationDownloads.md",
            "OpenTelemetryInstrumentationDistributed.md",
            "OpenTelemetryInstrumentationHTTP.md",
            "OpenTelemetryInstrumentationGenie.md",
        ],
    ],
)

deploydocs(repo = "github.com/oolong-dev/OpenTelemetry.jl.git", target = "build")

ext: jl


text:
using Test
using EconLH.ProductionFunctionsLH

pf = ProductionFunctionsLH;

function prod_fct_test(fS)
    println(fS)
    n = n_inputs(fS);
    T = length(productivities(fS));
    if T == 1
        xM = LinRange(1.2, 0.4, n)';
    else
        xM = LinRange(0.5, 10.5, T) * LinRange(1.2, 0.4, n)';
    end

    yV = output(fS, xM);
    @test size(yV) == (T,)
    @test all(yV .> 0.0)

    if constant_returns(fS)
        y2V = output(fS, xM .* 2.0);
        @test all(isapprox(yV .* 2.0, y2V, atol = 1e-4));
        for t = 1 : T
            xMin, xMax = extrema(xM[t,:]);
            @test yV[t] / productivities(fS)[t] > xMin
            @test yV[t] / productivities(fS)[t] < xMax
        end
    end

    # Marginal products
    mpM = mproducts(fS, xM);
    dx = 0.00001;
    for j = 1 : n
      x2M = copy(xM);
      x2M[:,j] = x2M[:,j] .+ dx;
      y2V = output(fS, x2M);
      derivV = (y2V .- yV) ./ dx;
      @test  all(abs.(derivV .- mpM[:,j]) ./ max.(0.1, mpM[:,j]) .< 1e-3)
      if !all(abs.(derivV .- mpM[:,j]) ./ max.(0.1, mpM[:,j]) .< 1e-3)
        println(derivV)
        println(mpM[:, j])
      end
    end
end

function output_direct_test(substElast)
    @testset "Direct output computation" begin
        n = 3;
        T = 4;
        alphaV = (1 : n) ./ n;
        xM = (1:T) .+ (1:n)' .+ 0.5;
        yV = ces_output(xM, alphaV, substElast);
        @test all(yV .> 0.0);
        @test size(yV) == (T,);

        # One set of inputs
        y = ces_output(xM[2,:], alphaV, substElast);
        @test isapprox(y, yV[2]);
    end
end

@testset "Production Functions" begin
    n = 3;
    substElastV = (0.5, 0.99, 1.0, 1.01, 3.0);
    TV = (1, 5);
    for T in TV
        for substElast in substElastV
            fS = pf.make_test_ces(T, n, substElast);
            prod_fct_test(fS);
            output_direct_test(substElast);
        end
    end

end

# -------------


text:
module Smolyak

using Iterators

export SmolyakGrid, build_B, build_B!, build_grid,
       SmolyakInterp, interpolate, find_theta, update_theta

include("util.jl")
include("grid.jl")
include("interp.jl")

end # module

max_issues_repo_name: EconForge/Smolyak

max_forks_repo_name: EconForge/Smolyak


text:
@testset "Write to IOBuffer" begin
    csv_string = chomp(read(joinpath(dir, "test_basic.csv"), String))
    df = CSV.read(IOBuffer(csv_string))
    io = IOBuffer()
    CSV.write(io, df)
    written = chomp(String(take!(io)))
    @test written == csv_string
end


text:
# ------------
# Construction
# ------------

export seq, constantly, repeatedly, iterate, concat

seq(xs::List) = xs
seq(xs::Array) = isempty(xs) ? list() : xs[1]:seq(xs[2:end])

seq(xs::Tuple) = seq(collect(xs))

seq(itr) = seq(itr, start(itr))
seq(itr, state) =
  @lazy done(itr, state) ? [] :
    begin x, state = next(itr, state)
      prepend(x, seq(itr, state))
    end

constantly(x) = @lazy x:constantly(x)
constantly(n, x) = @>> constantly(x) take(n)

cycle(xs) = @lazy xs * cycle(xs)

repeatedly(f) = @lazy f():repeatedly(f)
repeatedly(n, f) = @>> repeatedly(f) take(n)

iterate(f, v) = @lazy v:iterate(f, f(v))

range(x, y, step=1) =
  @lazy x <= y ? (x:range(x+step, y, step)) : []
range(x=1) = # Optimisation for y=Inf
  @lazy x:range(x+1)

concat(xs::List, ys::List) =
  @lazy isempty(xs) ? ys : (first(xs):concat(tail(xs), ys))

*(xs::List, ys::List) = concat(xs, ys)
*(xs::BitArray, ys::List) = concat(seq(xs), ys)
*(xs::AbstractArray, ys::List) = concat(seq(xs), ys)
*(xs::List, ys::AbstractArray) = concat(xs, seq(ys))

# ------------
# Manipulation
# ------------

import Base: length, map, reduce, filter, reverse

if VERSION >= v"0.6.0-dev.1015"
  import Base.Iterators: drop, take
elseif VERSION >= v"0.4.0"
  import Base: drop, take
end

export riffle, interpose, take, drop, takelast, droplast, takenth, takewhile, dropwhile,
       lazymap, reductions, remove, dorun, foreach, distinct,
       groupby, partition, partitionby, splitat, splitby, flatten

riffle(ls...) = riffle(map(seq, ls)...)
riffle(ls::List...) =
  @lazy any(isempty, ls) ? [] :
    seq(map(first, ls)) * riffle(map(tail, ls)...)

interpose(xs, args...) = interpose(seq(xs), args...)
interpose(xs::List, y, n = 1) =
  @lazy isempty(xs) ? [] :
       take(n, xs) * (isempty(drop(n, xs)) ? [] :
                        prepend(y, interpose(drop(n, xs), y, n)))

length(l::List) = isempty(l) ? 0 : 1 + length(tail(l))

Base.endof(l::List) = error("Cant use `end` with List.")

take(n::Integer, l::List) =
  @lazy n <= 0 || isempty(l) ? [] : prepend(first(l), take(n-1, tail(l)))

drop(n::Integer, l::List) =
  @lazy n <= 0 ? l : drop(n-1, tail(l))

takelast(n::Integer, l::List) =
  @lazy isempty(drop(n, l)) ? l : takelast(n, tail(l))

Base.last(l::List) = @>> l takelast(1) first

droplast(n::Integer, l::List) =
  map((x,_)->x, l, drop(n,l))

takenth(n::Integer, l::List) =
  @lazy isempty(l) ? [] : first(drop(n-1,l)):takenth(n, drop(n, l))

for f in [:take :drop :takelast :droplast :takenth]
    # This avoid the ambiguity with the base versions
    @eval $f(l::List, n::Int) = $f(n, l)
    @eval $f(l::List, n::Integer) = $f(n, l)
end

takewhile(pred::Function, l::List) =
  @lazy isempty(l) || !pred(first(l)) ? [] : first(l):takewhile(pred, tail(l))

dropwhile(pred::Function, l::List) =
  @lazy isempty(l) || !pred(first(l)) ? l : dropwhile(pred, tail(l))

mapply(f::@compat(Union{Function, DataType}), ls...) =
  @lazy any(isempty, ls) ? [] : prepend(f(map(first, ls)...), mapply(f, map(tail, ls)...))

# Resolves amibguity error
map(f::Function, ls::List...) = mapply(f, ls...)
map(f::DataType, ls::List...) = mapply(f, ls...)

lazymap(f::@compat(Union{Function, DataType}), ls...) = map(f, map(seq, ls)...)

@rec reduce(f::Function, v, xs::List) =
  isempty(xs) ? v : reduce(f, f(v, first(xs)), tail(xs))

reduce(f::Function, xs::List) =
  isempty(xs) ? f() : reduce(f, first(xs), tail(xs))

reductions(f::Function, v, xs::List) =
  @lazy isempty(xs) ? [] : v:reductions(f, f(v, first(xs)), tail(xs))

reductions(f::Function, xs::List) =
  @lazy isempty(xs) ? [] : reductions(f, first(xs), tail(xs))

filter(f::Function, xs::List) =
  @lazy isempty(xs) ? [] :
        f(first(xs)) ? (first(xs):filter(f, tail(xs))) :
        filter(f, tail(xs))

remove(f::Function, xs::List) = filter(x->!f(x), xs)

reverse(xs::List) = reduce((xs, x)->x:xs, list(), xs)

distinct(xs::List) = distinct(xs, Set())

distinct(xs::List, seen::Set) =
  @lazy isempty(xs) ? [] :
    first(xs) in seen ?
      distinct(tail(xs), seen) :
      first(xs):distinct(tail(xs), push!(seen, first(xs)))

function groupby(f, xs::List)
  groups = Dict()
  for x in xs
    k = f(x)
    groups[k] = x:get(groups, k, list())
  end
  return groups
end

partition(n, xs::List; step = n, pad = nothing) =
  @lazy isempty(xs) ? [] :
    @with (l = take(n, xs), len = length(l)),
      len < n ?
        (pad == nothing ? [] : list(l * take(n-len, pad))) :
        (l:partition(n, drop(step, xs); step = step, pad = pad))

partitionby(f, xs::List) =
  @lazy isempty(xs) ? [] :
    @with (x = first(xs), v = f(x), run = takewhile(x->f(x)==v, tail(xs))),
      prepend(x,run):partitionby(f, @lazy drop(length(run)+1, xs))

splitat(n, xs::List) = (take(n, xs), drop(n, xs))

splitby(p::Function, xs::List) = takewhile(p, xs), dropwhile(p, xs)

walk(inner, outer, xs::List) = @>> xs map(inner) outer
walk(inner, outer, x) = outer(x)

prewalk(f, xs)  = walk(x -> prewalk(f, x), identity, f(xs))
postwalk(f, xs) = walk(x -> postwalk(f, x), f, xs)

flatten(x) = list(x)
flatten(xs::List) = reduce((xs, x) -> xs*flatten(x), list(), xs)

# ----------
# Predicates
# ----------

import Base: any, all

==(xs::List, ys::List) =
  isempty(xs) == isempty(ys) &&
    (isempty(xs) || first(xs) == first(ys) && tail(xs) == tail(ys))

if isdefined(Base, :Predicate)
  any(f::Base.Predicate, xs::List) = @>> xs map(f) any
  all(f::Base.Predicate, xs::List) = @>> xs map(f) all
else
  any(f, xs::List) = @>> xs map(f) any
  all(f, xs::List) = @>> xs map(f) all
end

@rec any(xs::List) = isempty(xs) ? false : first(xs) || any(tail(xs))
any(::typeof(@functorize(identity)), xs::List) = any(xs)

@rec all(xs::List) = isempty(xs) ? true : first(xs) && all(tail(xs))
all(::typeof(@functorize(identity)), xs::List) = all(xs)


text:
#=
# Examples
=#

import Pasteee

#=
To access the Paste.ee API, you need an Application Key.
You can create yours at the following link:
<https://paste.ee/account/api>
(after logging in to [Paste.ee](https://paste.ee/))
In that page, first create an Application,
then click ACTIONS -> AUTHORIZATION PAGE, and then SUBMIT.
You will see the Application Key.
=#

const appkey = ENV["PASTEEE_APPKEY"]; nothing #hide

#=
In the following examples I assume that you have assigned your
Application Key to the `appkey` variable.
=#

# Create a paste that expires in one hour.

id = Pasteee.paste(appkey, "Hola mundo"; expiration="3600")

# The paste can be retrieved using the returned `id`.

paste = Pasteee.get(appkey, id)

# Paste.ee pastes are organized in sections.
# Here we retrieve the contents of section number 1 of the paste we just created.

paste["sections"][1]["contents"]

# Delete the paste.

Pasteee.delete(appkey, id)


text:
function init_pseudo_observable_mappings!(m::Model1002)

    pseudo_names = [:y_t, :y_f_t, :NaturalRate, :π_t, :OutputGap, :ExAnteRealRate, :LongRunInflation,
                    :MarginalCost, :Wages, :FlexibleWages, :Hours, :FlexibleHours, :z_t,
                    :Expected10YearRateGap, :NominalFFR, :Expected10YearRate, :Expected10YearNaturalRate,
                    :ExpectedNominalNaturalRate, :NominalRateGap, :LaborProductivityGrowth]

    # Create PseudoObservable objects
    pseudo = OrderedDict{Symbol,PseudoObservable}()
    for k in pseudo_names
        pseudo[k] = PseudoObservable(k)
    end

    # Fill in names and reverse transforms
    pseudo[:y_t].name = "Output Growth"
    pseudo[:y_t].longname = "Output Growth Per Capita"

    pseudo[:y_f_t].name = "Flexible Output Growth"
    pseudo[:y_f_t].longname = "Output that would prevail in a flexible-price economy."

    pseudo[:NaturalRate].name = "Real Natural Rate"
    pseudo[:NaturalRate].longname = "The real interest rate that would prevail in a flexible-price economy."
    pseudo[:NaturalRate].rev_transform = quartertoannual

    pseudo[:π_t].name = "Inflation"
    pseudo[:π_t].longname = "Inflation"
    pseudo[:π_t].rev_transform = quartertoannual

    pseudo[:OutputGap].name = "Output Gap"
    pseudo[:OutputGap].longname = "Output Gap"

    pseudo[:ExAnteRealRate].name = "Ex Ante Real Rate"
    pseudo[:ExAnteRealRate].longname = "Ex Ante Real Rate"
    pseudo[:ExAnteRealRate].rev_transform = quartertoannual

    pseudo[:LongRunInflation].name = "Long Run Inflation"
    pseudo[:LongRunInflation].longname = "Long Run Inflation"
    pseudo[:LongRunInflation].rev_transform = quartertoannual

    pseudo[:MarginalCost].name = "Marginal Cost"
    pseudo[:MarginalCost].longname = "Marginal Cost"

    pseudo[:Wages].name = "Wages"
    pseudo[:Wages].longname = "Wages"

    pseudo[:FlexibleWages].name = "Flexible Wages"
    pseudo[:FlexibleWages].longname = "Wages that would prevail in a flexible-wage economy"

    pseudo[:Hours].name = "Hours"
    pseudo[:Hours].longname = "Hours"

    pseudo[:FlexibleHours].name     = "Flexible Hours"
    pseudo[:FlexibleHours].longname = "Flexible Hours"

    pseudo[:z_t].name     = "z_t"
    pseudo[:z_t].longname = "z_t"

    pseudo[:Expected10YearRateGap].name     = "Expected 10-Year Rate Gap"
    pseudo[:Expected10YearRateGap].longname = "Expected 10-Year Rate Gap"
    pseudo[:Expected10YearRateGap].rev_transform = quartertoannual

    pseudo[:NominalFFR].name     = "Nominal FFR"
    pseudo[:NominalFFR].longname = "Nominal FFR at an annual rate"
    pseudo[:NominalFFR].rev_transform = quartertoannual

    pseudo[:Expected10YearRate].name     = "Expected 10-Year Rate"
    pseudo[:Expected10YearRate].longname = "Expected 10-Year Interest Rate"
    pseudo[:Expected10YearRate].rev_transform = quartertoannual

    pseudo[:Expected10YearNaturalRate].name     = "Expected 10-Year Natural Rate"
    pseudo[:Expected10YearNaturalRate].longname = "Expected 10-Year Natural Rate of Interest"
    pseudo[:Expected10YearNaturalRate].rev_transform = quartertoannual

    pseudo[:ExpectedNominalNaturalRate].name     = "Expected Nominal Natural Rate"
    pseudo[:ExpectedNominalNaturalRate].longname = "Natural Rate + Expected Inflation"
    pseudo[:ExpectedNominalNaturalRate].rev_transform = quartertoannual

    pseudo[:NominalRateGap].name     = "Nominal Rate Gap"
    pseudo[:NominalRateGap].longname = "Nominal FFR - Nominal Natural Rate"
    pseudo[:NominalRateGap].rev_transform = quartertoannual

    pseudo[:LaborProductivityGrowth].name     = "Labor Productivity Growth"
    pseudo[:LaborProductivityGrowth].longname = "Labor Productivity Growth Rate"
    pseudo[:LaborProductivityGrowth].rev_transform = quartertoannual

    # Add to model object
    m.pseudo_observable_mappings = pseudo
end


text:
module Printf
using Base.Grisu
export @printf, @sprintf

### printf formatter generation ###

function gen(s::String)
    args = {}
    blk = Expr(:block, :(local neg, pt, len, exp))
    for x in parse(s)
        if isa(x,String)
            push!(blk.args, :(write(out, $(length(x)==1 ? x[1] : x))))
        else
            c = lowercase(x[end])
            f = c=='f' ? gen_f :
                c=='e' ? gen_e :
                c=='g' ? gen_g :
                c=='c' ? gen_c :
                c=='s' ? gen_s :
                c=='p' ? gen_p :
                         gen_d
            arg, ex = f(x...)
            push!(args, arg)
            push!(blk.args, ex)
        end
    end
    push!(blk.args, :nothing)
    return args, blk
end

### printf format string parsing ###

function parse(s::String)
    # parse format string in to stings and format tuples
    list = {}
    i = j = start(s)
    while !done(s,j)
        c, k = next(s,j)
        if c == '%'
            isempty(s[i:j-1]) || push!(list, s[i:j-1])
            flags, width, precision, conversion, k = parse1(s,k)
            '\'' in flags && error("printf format flag ' not yet supported")
            conversion == 'a'    && error("printf feature %a not yet supported")
            conversion == 'n'    && error("printf feature %n not supported")
            push!(list, conversion == '%' ? "%" : (flags,width,precision,conversion))
            i = j = k
        else
            j = k
        end
    end
    isempty(s[i:end]) || push!(list, s[i:end])
    # coalesce adjacent strings
    i = 1
    while i < length(list)
        if isa(list[i],String)
            for j = i+1:length(list)
                if !isa(list[j],String)
                    j -= 1
                    break
                end
                list[i] *= list[j]
            end
            splice!(list,i+1:j)
        end
        i += 1
    end
    return list
end

## parse a single printf specifier ##

# printf specifiers:
#   %                       # start
#   (\d+\$)?                # arg (not supported)
#   [\-\+#0' ]*             # flags
#   (\d+)?                  # width
#   (\.\d*)?                # precision
#   (h|hh|l|ll|L|j|t|z|q)?  # modifier (ignored)
#   [diouxXeEfFgGaAcCsSp%]  # conversion

next_or_die(s::String, k) = !done(s,k) ? next(s,k) :
    error("invalid printf format string: ", repr(s))

function parse1(s::String, k::Integer)
    j = k
    width = 0
    precision = -1
    c, k = next_or_die(s,k)
    # handle %%
    if c == '%'
        return "", width, precision, c, k
    end
    # parse flags
    while c in "#0- + '"
        c, k = next_or_die(s,k)
    end
    flags = ascii(s[j:k-2])
    # parse width
    while '0' <= c <= '9'
        width = 10*width + c-'0'
        c, k = next_or_die(s,k)
    end
    # parse precision
    if c == '.'
        c, k = next_or_die(s,k)
        if '0' <= c <= '9'
            precision = 0
            while '0' <= c <= '9'
                precision = 10*precision + c-'0'
                c, k = next_or_die(s,k)
            end
        end
    end
    # parse length modifer (ignored)
    if c == 'h' || c == 'l'
        prev = c
        c, k = next_or_die(s,k)
        if c == prev
            c, k = next_or_die(s,k)
        end
    elseif c in "Ljqtz"
        c, k = next_or_die(s,k)
    end
    # validate conversion
    if !(c in "diouxXDOUeEfFgGaAcCsSpn")
        error("invalid printf format string: ", repr(s))
    end
    # TODO: warn about silly flag/conversion combinations
    flags, width, precision, c, k
end

### printf formatter generation ###

function special_handler(flags::ASCIIString, width::Int)
    @gensym x
    blk = Expr(:block)
    pad = '-' in flags ? rpad : lpad
    pos = '+' in flags ? "+" :
          ' ' in flags ? " " : ""
    abn = quote
        isnan($x) ? $(pad("NaN", width)) :
         $x < 0   ? $(pad("-Inf", width)) :
                    $(pad("$(pos)Inf", width))
    end
    ex = :(isfinite($x) ? $blk : write(out, $abn))
    x, ex, blk
end

function pad(m::Int, n, c::Char)
    if m <= 1
        :($n > 0 && write(out,$c))
    else
        @gensym i
        quote
            $i = $n
            while $i > 0
                write(out,$c)
                $i -= 1
            end
        end
    end
end

function print_fixed(out, precision)
    pdigits = pointer(DIGITS)
    ndigits = LEN[1]
    pt = POINT[1]
    if pt <= 0
        # 0.0dddd0
        write(out, '0')
        write(out, '.')
        precision += pt
        while pt < 0
            write(out, '0')
            pt += 1
        end
        write(out, pdigits, ndigits)
        precision -= ndigits
    elseif ndigits <= pt
        # dddd000.000000
        write(out, pdigits, ndigits)
        while ndigits < pt
            write(out, '0')
            ndigits += 1
        end
        write(out, '.')
    else # 0 < pt < ndigits
        # dd.dd0000
        ndigits -= pt
        write(out, pdigits, pt)
        write(out, '.')
        write(out, pdigits+pt, ndigits)
        precision -= ndigits
    end
    while precision > 0
        write(out, '0')
        precision -= 1
    end
end

function print_exp(out, exp)
    write(out, exp < 0 ? '-' : '+')
    exp = abs(exp)
    d = div(exp,100)
    d > 0 && write(out, char('0'+d))
    exp = rem(exp,100)
    write(out, char('0'+div(exp,10)))
    write(out, char('0'+rem(exp,10)))
end

function gen_d(flags::ASCIIString, width::Int, precision::Int, c::Char)
    # print integer:
    #  [dDiu]: print decimal digits
    #  [o]:    print octal digits
    #  [x]:    print hex digits, lowercase
    #  [X]:    print hex digits, uppercase
    #
    # flags:
    #  (#): prefix hex with 0x/0X; octal leads with 0
    #  (0): pad left with zeros
    #  (-): left justify
    #  ( ): precede non-negative values with " "
    #  (+): precede non-negative values with "+"
    #
    x, ex, blk = special_handler(flags,width)
    # interpret the number
    prefix = ""
    if lowercase(c)=='o'
        f = '#' in flags ? :int_0ct : :int_oct
        push!(blk.args, :(($f)($x)))
    elseif c=='x'
        '#' in flags && (prefix = "0x")
        push!(blk.args, :(int_hex($x)))
    elseif c=='X'
        '#' in flags && (prefix = "0X")
        push!(blk.args, :(int_HEX($x)))
    else
        push!(blk.args, :(int_dec($x)))
    end
    push!(blk.args, :(neg = NEG[1]))
    push!(blk.args, :(pt  = POINT[1]))
    # calculate padding
    width -= length(prefix)
    space_pad = width > max(1,precision) && '-' in flags ||
                precision < 0 && width > 1 && !('0' in flags) ||
                precision >= 0 && width > precision
    padding = nothing
    if precision < 1; precision = 1; end
    if space_pad
        if '+' in flags || ' ' in flags
            width -= 1
            if width > precision
                padding = :($width-(pt > $precision ? pt : $precision))
            end
        else
            if width > precision
                padding = :($width-neg-(pt > $precision ? pt : $precision))
            end
        end
    end
    # print space padding
    if padding != nothing && !('-' in flags)
        push!(blk.args, pad(width-precision, padding, ' '))
    end
    # print sign
    '+' in flags ? push!(blk.args, :(write(out, neg?'-':'+'))) :
    ' ' in flags ? push!(blk.args, :(write(out, neg?'-':' '))) :
                   push!(blk.args, :(neg && write(out, '-')))
    # print prefix
    for ch in prefix
        push!(blk.args, :(write(out, $ch)))
    end
    # print zero padding & leading zeros
    if space_pad && precision > 1
        push!(blk.args, pad(precision-1, :($precision-pt), '0'))
    elseif !space_pad && width > 1
        zeros = '+' in flags || ' ' in flags ? :($(width-1)-pt) : :($width-neg-pt)
        push!(blk.args, pad(width-1, zeros, '0'))
    end
    # print integer
    push!(blk.args, :(write(out, pointer(DIGITS), pt)))
    # print padding
    if padding != nothing && '-' in flags
        push!(blk.args, pad(width-precision, padding, ' '))
    end
    # return arg, expr
    :(($x)::Real), ex
end

function gen_f(flags::ASCIIString, width::Int, precision::Int, c::Char)
    # print to fixed trailing precision
    #  [fF]: the only choice
    #
    # flags
    #  (#): always print a decimal point
    #  (0): pad left with zeros
    #  (-): left justify
    #  ( ): precede non-negative values with " "
    #  (+): precede non-negative values with "+"
    #
    x, ex, blk = special_handler(flags,width)
    # interpret the number
    if precision < 0; precision = 6; end
    push!(blk.args, :(fix_dec($x,$precision)))
    push!(blk.args, :(neg = NEG[1]))
    push!(blk.args, :(pt  = POINT[1]))
    push!(blk.args, :(len = LEN[1]))
    # calculate padding
    padding = nothing
    if precision > 0 || '#' in flags
        width -= precision+1
    end
    if '+' in flags || ' ' in flags
        width -= 1
        if width > 1
            padding = :($width-(pt > 0 ? pt : 1))
        end
    else
        if width > 1
            padding = :($width-(pt > 0 ? pt : 1)-neg)
        end
    end
    # print space padding
    if padding != nothing && !('-' in flags) && !('0' in flags)
        push!(blk.args, pad(width-1, padding, ' '))
    end
    # print sign
    '+' in flags ? push!(blk.args, :(write(out, neg?'-':'+'))) :
    ' ' in flags ? push!(blk.args, :(write(out, neg?'-':' '))) :
                   push!(blk.args, :(neg && write(out, '-')))
    # print zero padding
    if padding != nothing && !('-' in flags) && '0' in flags
        push!(blk.args, pad(width-1, padding, '0'))
    end
    # print digits
    if precision > 0
        push!(blk.args, :(print_fixed(out,$precision)))
    else
        push!(blk.args, :(write(out, pointer(DIGITS), len)))
        push!(blk.args, :(while pt >= (len+=1) write(out,'0') end))
        '#' in flags && push!(blk.args, :(write(out, '.')))
    end
    # print space padding
    if padding != nothing && '-' in flags
        push!(blk.args, pad(width-1, padding, ' '))
    end
    # return arg, expr
    :(($x)::Real), ex
end

function gen_e(flags::ASCIIString, width::Int, precision::Int, c::Char)
    # print float in scientific form:
    #  [e]: use 'e' to introduce exponent
    #  [E]: use 'E' to introduce exponent
    #
    # flags:
    #  (#): always print a decimal point
    #  (0): pad left with zeros
    #  (-): left justify
    #  ( ): precede non-negative values with " "
    #  (+): precede non-negative values with "+"
    #
    x, ex, blk = special_handler(flags,width)
    # interpret the number
    if precision < 0; precision = 6; end
    ndigits = min(precision+1,BUFLEN-1)
    push!(blk.args, :(ini_dec($x,$ndigits)))
    push!(blk.args, :(neg = NEG[1]))
    push!(blk.args, :(exp = POINT[1]-1))
    expmark = c=='E' ? "E" : "e"
    if precision==0 && '#' in flags
        expmark = string(".",expmark)
    end
    # calculate padding
    padding = nothing
    width -= precision+length(expmark)+(precision>0)+4
    # 4 = leading + expsign + 2 exp digits
    if '+' in flags || ' ' in flags
        width -= 1 # for the sign indicator
        if width > 0
            padding = :($width-((exp<=-100)|(100<=exp)))
        end
    else
        if width > 0
            padding = :($width-((exp<=-100)|(100<=exp))-neg)
        end
    end
    # print space padding
    if padding != nothing && !('-' in flags) && !('0' in flags)
        push!(blk.args, pad(width, padding, ' '))
    end
    # print sign
    '+' in flags ? push!(blk.args, :(write(out, neg?'-':'+'))) :
    ' ' in flags ? push!(blk.args, :(write(out, neg?'-':' '))) :
                    push!(blk.args, :(neg && write(out, '-')))
    # print zero padding
    if padding != nothing && !('-' in flags) && '0' in flags
        push!(blk.args, pad(width, padding, '0'))
    end
    # print digits
    push!(blk.args, :(write(out, DIGITS[1])))
    if precision > 0
        push!(blk.args, :(write(out, '.')))
        push!(blk.args, :(write(out, pointer(DIGITS)+1, $(ndigits-1))))
        if ndigits < precision+1
            n = precision+1-ndigits
            push!(blk.args, pad(n, n, '0'))
        end
    end
    for ch in expmark
        push!(blk.args, :(write(out, $ch)))
    end
    push!(blk.args, :(print_exp(out, exp)))
    # print space padding
    if padding != nothing && '-' in flags
        push!(blk.args, pad(width, padding, ' '))
    end
    # return arg, expr
    :(($x)::Real), ex
end

function gen_c(flags::ASCIIString, width::Int, precision::Int, c::Char)
    # print a character:
    #  [cC]: both the same for us (Unicode)
    #
    # flags:
    #  (0): pad left with zeros
    #  (-): left justify
    #
    @gensym x
    blk = Expr(:block, :($x = char($x)))
    if width > 1 && !('-' in flags)
        p = '0' in flags ? '0' : ' '
        push!(blk.args, pad(width-1, :($width-charwidth($x)), p))
    end
    push!(blk.args, :(write(out, $x)))
    if width > 1 && '-' in flags
        push!(blk.args, pad(width-1, :($width-charwidth($x)), ' '))
    end
    :(($x)::Integer), blk
end

function gen_s(flags::ASCIIString, width::Int, precision::Int, c::Char)
    # print a string:
    #  [sS]: both the same for us (Unicode)
    #
    # flags:
    #  (0): pad left with zeros
    #  (-): left justify
    #
    @gensym x
    blk = Expr(:block)
    if width > 0
        if !('#' in flags)
            push!(blk.args, :($x = string($x)))
        else
            push!(blk.args, :($x = repr($x)))
        end
        if !('-' in flags)
            push!(blk.args, pad(width, :($width-strwidth($x)), ' '))
        end
        push!(blk.args, :(write(out, $x)))
        if '-' in flags
            push!(blk.args, pad(width, :($width-strwidth($x)), ' '))
        end
    else
        if !('#' in flags)
            push!(blk.args, :(print(out, $x)))
        else
            push!(blk.args, :(show(out, $x)))
        end
    end
    :(($x)::Any), blk
end

# TODO: faster pointer printing.

function gen_p(flags::ASCIIString, width::Int, precision::Int, c::Char)
    # print pointer:
    #  [p]: the only option
    #
    @gensym x
    blk = Expr(:block)
    ptrwidth = WORD_SIZE>>2
    width -= ptrwidth+2
    if width > 0 && !('-' in flags)
        push!(blk.args, pad(width, width, ' '))
    end
    push!(blk.args, :(write(out, '0')))
    push!(blk.args, :(write(out, 'x')))
    push!(blk.args, :(write(out, bytestring(hex(unsigned($x), $ptrwidth)))))
    if width > 0 && '-' in flags
        push!(blk.args, pad(width, width, ' '))
    end
    :(($x)::Ptr), blk
end

function gen_g(flags::ASCIIString, width::Int, precision::Int, c::Char)
    error("printf \"%g\" format specifier not implemented")
end

### core unsigned integer decoding functions ###

macro handle_zero(ex)
    quote
        if $(esc(ex)) == 0
            POINT[1] = 1
            DIGITS[1] = '0'
            return
        end
    end
end

function decode_oct(x::Unsigned)
    @handle_zero x
    POINT[1] = i = div((sizeof(x)<<3)-leading_zeros(x)+2,3)
    while i > 0
        DIGITS[i] = '0'+(x&0x7)
        x >>= 3
        i -= 1
    end
end

function decode_0ct(x::Unsigned)
    # doesn't need special handling for zero
    POINT[1] = i = div((sizeof(x)<<3)-leading_zeros(x)+5,3)
    while i > 0
        DIGITS[i] = '0'+(x&0x7)
        x >>= 3
        i -= 1
    end
end

function decode_dec(x::Unsigned)
    @handle_zero x
    POINT[1] = i = Base.ndigits0z(x)
    while i > 0
        DIGITS[i] = '0'+rem(x,10)
        x = div(x,10)
        i -= 1
    end
end

function decode_hex(x::Unsigned, symbols::Array{Uint8,1})
    @handle_zero x
    POINT[1] = i = (sizeof(x)<<1)-(leading_zeros(x)>>2)
    while i > 0
        DIGITS[i] = symbols[(x&0xf)+1]
        x >>= 4
        i -= 1
    end
end

const hex_symbols = "0123456789abcdef".data
const HEX_symbols = "0123456789ABCDEF".data

decode_hex(x::Unsigned) = decode_hex(x,hex_symbols)
decode_HEX(x::Unsigned) = decode_hex(x,HEX_symbols)

function decode(b::Int, x::BigInt)
    neg = NEG[1] = x.size < 0
    pt = POINT[1] = Base.ndigits(x, abs(b))
    length(DIGITS) < pt+1 && resize!(DIGITS, pt+1)
    neg && (x.size = -x.size)
    ccall((:__gmpz_get_str, :libgmp), Ptr{Uint8},
          (Ptr{Uint8}, Cint, Ptr{BigInt}), DIGITS, b, &x)
    neg && (x.size = -x.size)
end

function decode_0ct(x::BigInt)
    neg = NEG[1] = x.size < 0
    DIGITS[1] = '0'
    if x.size == 0
        POINT[1] = 1
        return
    end
    pt = POINT[1] = Base.ndigits0z(x, 8) + 1
    length(DIGITS) < pt+1 && resize!(DIGITS, pt+1)
    neg && (x.size = -x.size)
    p = convert(Ptr{Uint8}, DIGITS) + 1
    ccall((:__gmpz_get_str, :libgmp), Ptr{Uint8},
          (Ptr{Uint8}, Cint, Ptr{BigInt}), p, 8, &x)
    neg && (x.size = -x.size)
end

### decoding functions directly used by printf generated code ###

# int_*(x)   => fixed precision, to 0th place, filled out
# fix_*(x,n) => fixed precision, to nth place, not filled out
# ini_*(x,n) => n initial digits, filled out

# alternate versions:
#   *_0ct(x,n) => ensure that the first octal digits is zero
#   *_HEX(x,n) => use uppercase digits for hexadecimal

## "int" decoding functions ##
#
# - sets neg[1]
# - sets point[1]
# - implies len[1] = point[1]
#

int_oct(x::Unsigned) = (NEG[1]=false; decode_oct(x))
int_0ct(x::Unsigned) = (NEG[1]=false; decode_0ct(x))
int_dec(x::Unsigned) = (NEG[1]=false; decode_dec(x))
int_hex(x::Unsigned) = (NEG[1]=false; decode_hex(x))
int_HEX(x::Unsigned) = (NEG[1]=false; decode_HEX(x))

macro handle_negative()
    quote
        if $(esc(:x)) < 0
            NEG[1] = true
            $(esc(:x)) = oftype($(esc(:x)),-$(esc(:x)))
        else
            NEG[1] = false
        end
    end
end

int_oct(x::Integer) = (@handle_negative; decode_oct(unsigned(x)))
int_0ct(x::Integer) = (@handle_negative; decode_0ct(unsigned(x)))
int_dec(x::Integer) = (@handle_negative; decode_dec(unsigned(x)))
int_hex(x::Integer) = (@handle_negative; decode_hex(unsigned(x)))
int_HEX(x::Integer) = (@handle_negative; decode_HEX(unsigned(x)))

int_oct(x::BigInt) = decode(8, x)
int_0ct(x::BigInt) = decode_0ct(x)
int_dec(x::BigInt) = decode(10, x)
int_hex(x::BigInt) = decode(16, x)
int_HEX(x::BigInt) = decode(-16, x)

const SmallFloatingPoint = Union(Float64,Float32,Float16)

function int_dec(x::SmallFloatingPoint)
    if x == 0.0
        NEG[1] = false
        POINT[1] = 1
        DIGITS[1] = '0'
        return
    end
    @grisu_ccall x Grisu.FIXED 0
    if LEN[1] == 0
        NEG[1] = false
        POINT[1] = 1
        DIGITS[1] = '0'
    else
        for i = LEN[1]+1:POINT[1]
            DIGITS[i] = '0'
        end
    end
end

int_oct(x::Real) = int_oct(integer(x)) # TODO: real float decoding.
int_0ct(x::Real) = int_0ct(integer(x)) # TODO: real float decoding.
int_dec(x::Real) = int_dec(float(x))
int_hex(x::Real) = int_hex(integer(x)) # TODO: real float decoding.
int_HEX(x::Real) = int_HEX(integer(x)) # TODO: real float decoding.

## fix decoding functions ##
#
# - sets neg[1]
# - sets point[1]
# - sets len[1]; if less than point[1], trailing zeros implied
#

fix_dec(x::Integer, n::Int) = (int_dec(x); LEN[1]=POINT[1])
fix_dec(x::Real, n::Int) = fix_dec(float(x),n)

function fix_dec(x::SmallFloatingPoint, n::Int)
    if n > BUFLEN-1; n = BUFLEN-1; end
    @grisu_ccall x Grisu.FIXED n
    if LEN[1] == 0
        NEG[1] = false
        POINT[1] = 1
        DIGITS[1] = '0'
    end
end

## ini decoding functions ##
#
# - sets neg[1]
# - sets point[1]
# - implies len[1] = n (requested digits)
#

function ini_dec(x::Unsigned, n::Int)
    k = ndigits(x)
    if k <= n
        POINT[1] = k
        for i = k:-1:1
            DIGITS[i] = '0'+rem(x,10)
            x = div(x,10)
        end
        for i = k+1:n
            DIGITS[i] = '0'
        end
    else
        p = Base.powers_of_ten[k-n+1]
        r = rem(x,p)
        if r >= (p>>1)
            x += p
            if x >= Base.powers_of_ten[k+1]
                p *= 10
                k += 1
            end
        end
        POINT[1] = k
        x = div(x,p)
        for i = n:-1:1
            DIGITS[i] = '0'+rem(x,10)
            x = div(x,10)
        end
    end
end

ini_dec(x::Integer, n::Int) = (@handle_negative; ini_dec(unsigned(x),n))
ini_dec(x::Real, n::Int) = ini_dec(float(x),n)

function ini_dec(x::SmallFloatingPoint, n::Int)
    if x == 0.0
        POINT[1] = 1
        NEG[1] = signbit(x)
        ccall(:memset, Ptr{Void}, (Ptr{Void}, Cint, Csize_t), DIGITS, '0', n)
    else
        @grisu_ccall x Grisu.PRECISION n
    end
end

function ini_dec(x::BigInt, n::Int)
    if x.size == 0
        POINT[1] = 1
        NEG[1] = false
        ccall(:memset, Ptr{Void}, (Ptr{Void}, Cint, Csize_t), DIGITS, '0', n)
    else
        d = Base.ndigits0z(x)
        if d <= n
            int_dec(x)
            d == n && return
            p = convert(Ptr{Void}, DIGITS) + POINT[1]
            ccall(:memset, Ptr{Void}, (Ptr{Void}, Cint, Csize_t), p, '0', n - POINT[1])
        else
            int_dec(iround(x/big(10)^(d-n)))
            POINT[1] = d
        end
    end
end

### external printf interface ###

is_str_expr(ex) =
    isa(ex,Expr) && (ex.head == :string || (ex.head == :macrocall && isa(ex.args[1],Symbol) &&
    endswith(string(ex.args[1]),"str")))

function _printf(macroname, io, fmt, args)
    isa(fmt, String) || error("$macroname: format must be a plain static string (no interpolation or prefix)")
    sym_args, blk = gen(fmt)
    if length(sym_args) != length(args)
        error("$macroname: wrong number of arguments")
    end
    for i = length(args):-1:1
        var = sym_args[i].args[1]
        unshift!(blk.args, :($var = $(esc(args[i]))))
    end
    unshift!(blk.args, :(out = $io))
    blk
end

macro printf(args...)
    !isempty(args) || error("@printf: called with zero arguments")
    if isa(args[1], String) || is_str_expr(args[1])
        _printf("@printf", :STDOUT, args[1], args[2:end])
    else
        (length(args) >= 2 && (isa(args[2], String) || is_str_expr(args[2]))) ||
            error("@printf: first or second argument must be a format string")
        _printf("@printf", esc(args[1]), args[2], args[3:end])
    end
end

macro sprintf(args...)
    !isempty(args) || error("@sprintf: called with zero arguments")
    isa(args[1], String) || is_str_expr(args[1]) ||
        error("@sprintf: first argument must be a format string")
    :(sprint(io->$(_printf("@sprintf", :io, args[1], args[2:end]))))
end

end # module

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia


text:
module TestScikitLearn

# using Revise
using MLJBase
using Test
using LinearAlgebra

import MLJModels
import ScikitLearn
using MLJModels.ScikitLearn_
using CategoricalArrays

## CLASSIFIERS
@test_logs (:warn,"kernel parameter is not valid, setting to default=\"rbf\" \n") SVMClassifier(kernel="wrong")
@test_logs (:warn,"penalty parameter is not valid, setting to default=\"l2\" \n") SVMLClassifier(penalty="wrong")
@test_logs (:warn,"loss parameter is not valid, setting to default=\"epsilon_insensitive\" \n") SVMLRegressor(loss="wrong")

plain_classifier = SVMClassifier()
nu_classifier = SVMNuClassifier()
linear_classifier = SVMLClassifier(max_iter=10000)

# test preservation of categorical levels:
task = load_iris();
X, y = X_and_y(task)
train, test = partition(eachindex(y), 0.6) # levels of y are split across split

fitresultC, cacheC, reportC = MLJBase.fit(plain_classifier, 1,
                                          selectrows(X, train), y[train]);
fitresultCnu, cacheCnu, reportCnu = MLJBase.fit(nu_classifier, 1,
                                          selectrows(X, train), y[train]);
fitresultCL, cacheCL, reportCL = MLJBase.fit(linear_classifier, 1,
                                          selectrows(X, train), y[train]);
pcpred = predict(plain_classifier, fitresultC, selectrows(X, test));
nucpred = predict(nu_classifier, fitresultCnu, selectrows(X, test));
lcpred = predict(linear_classifier, fitresultCL, selectrows(X, test));

@test levels(pcpred) == levels(y[train])
@test levels(nucpred) == levels(y[train])
@test levels(lcpred) == levels(y[train])

# test with linear data:
x1 = randn(3000);
x2 = randn(3000);
x3 = randn(3000);
X = (x1=x1, x2=x2, x3=x3);
y = x1 - x2 -2x3;
ycat = map(y) do η
    η > 0 ? "go" : "stop"
end |> categorical;
train, test = partition(eachindex(ycat), 0.8);
fitresultC, cacheC, reportC = MLJBase.fit(plain_classifier, 1,
                                          selectrows(X, train), ycat[train]);
fitresultCnu, cacheCnu, reportCnu = MLJBase.fit(nu_classifier, 1,
                                          selectrows(X, train), ycat[train]);
fitresultCL, cacheCL, reportCL = MLJBase.fit(linear_classifier, 1,
                                          selectrows(X, train), ycat[train]);
pcpred = predict(plain_classifier, fitresultC, selectrows(X, test));
nucpred = predict(nu_classifier, fitresultCnu, selectrows(X, test));
lcpred = predict(linear_classifier, fitresultCL, selectrows(X, test));
@test sum(pcpred .!= ycat[test])/length(ycat) < 0.05
@test sum(nucpred .!= ycat[test])/length(ycat) < 0.05
@test sum(lcpred .!= ycat[test])/length(ycat) < 0.05

## REGRESSORS

plain_regressor = SVMRegressor()
nu_regressor = SVMNuRegressor()
linear_regressor = SVMLRegressor(max_iter=10000)

# test with linear data:
fitresultR, cacheR, reportR = MLJBase.fit(plain_regressor, 1,
                                          selectrows(X, train), y[train]);
fitresultRnu, cacheRnu, reportRnu = MLJBase.fit(nu_regressor, 1,
                                                selectrows(X, train), y[train]);
fitresultRL, cacheRL, reportRL = MLJBase.fit(linear_regressor, 1,
                                             selectrows(X, train), y[train]);

rpred = predict(plain_regressor, fitresultR, selectrows(X, test));
nurpred = predict(nu_regressor, fitresultRnu, selectrows(X, test));
Lrpred = predict(linear_regressor, fitresultRL, selectrows(X, test));

@test norm(rpred - y[test])/sqrt(length(y)) < 0.2
@test norm(nurpred - y[test])/sqrt(length(y)) < 0.2
@test norm(Lrpred - y[test])/sqrt(length(y)) < 0.2

info(SVMClassifier)
info(SVMNuClassifier)
info(SVMLClassifier)
info(SVMRegressor)
info(SVMNuRegressor)
info(SVMLRegressor)

end
true


text:
struct LOCGBuffersSimple{Tx,TQ}
    x::Tx
    Ax::Tx
    x_prev::Tx
    Ax_prev::Tx
    r::Tx
    Ar::Tx
    c::Tx
    X::TQ
    Q::TQ
end
function LOCGBuffersSimple(x::Vector{T}) where T
    Q = Matrix{T}(length(x), 3)
    if all((x)->(x == zero(T)), x)
        x = rand(T, length(x))
    else
        x = copy(x)
    end
    LOCGBuffersSimple{Vector{T}, Matrix{T}}(x, similar(x), rand(T, length(x)), similar(x), similar(x), similar(x), similar(x), Q, similar(Q))
end
function LOCGBuffersSimple{Tx,TQ}(x::Tx, Q::TQ) where {Tx<:AbstractVector, TQ<:AbstractMatrix}
    length(x) == size(Q,1) || throw("Length of x is not equal to size(Q,2).")
    size(Q,2) == 3 || throw("size(Q,2) is not equal to 3.")
    if all((x)->(x == zero(eltype(x))), x)
        x = rand(eltype(x), length(x))
    else
        x = copy(x)
    end
    LOCGBuffersSimple{Tx,TQ}(x, similar(x), rand(eltype(x), length(x)), similar(x), similar(x), similar(x), similar(x), Q, similar(Q))
end

struct LOCGBuffersGeneral{Tx,TQ}
    x::Tx
    Ax::Tx
    Bx::Tx
    x_prev::Tx
    Ax_prev::Tx
    Bx_prev::Tx
    r::Tx
    Ar::Tx
    Br::Tx
    c::Tx
    X::TQ
    Q::TQ
end
function LOCGBuffersGeneral(x::Vector{T}) where T
    Q = Matrix{T}(length(x), 3)
    if all((x)->(x == zero(T)), x)
        x = rand(T, length(x))
    else
        x = copy(x)
    end
    LOCGBuffersGeneral{Vector{T}, Matrix{T}}(x, similar(x), similar(x), rand(T, length(x)), similar(x), similar(x), similar(x), similar(x), similar(x), similar(x), Q, similar(Q))
end
function LOCGBuffersGeneral{Tx,TQ}(x::Tx, Q::TQ) where {Tx<:AbstractVector, TQ<:AbstractMatrix}
    length(x) == size(Q,1) || throw("Length of x is not equal to size(Q,2).")
    size(Q,2) == 3 || throw("size(Q,2) is not equal to 3.")
    if all((x)->(x == zero(eltype(x))), x)
        x = rand(eltype(x), length(x))
    else
        x = copy(x)
    end
    LOCGBuffersGeneral{Tx,TQ}(x, similar(x), similar(x), rand(eltype(x), length(x)), similar(x), similar(x), similar(x), similar(x), similar(x), similar(x), Q, similar(Q))
end


text:
struct GeneDataView{G <: AbstractGene} <: AbstractArray{G, 1}
    parent::Chromosome{G}
    indices::Vector{UInt}
    property::Symbol
end

function Base.getproperty(gene::G, name::Symbol) where {G <: AbstractGene}
    if name in (:parent, :index)
        return getfield(gene, name)
    elseif name in propertynames(gene)
        return gene.parent.genedata[gene.index, name]
    else
        return missing
    end
end

function Base.getproperty(gv::GeneDataView{G}, name::Symbol) where {G <: AbstractGene}
    return getfield(gv, name)
end

function Base.getproperty(genes::AbstractArray{G, 1}, name::Symbol) where {G <: AbstractGene}
    # @assert all(getfield.(genes, :parent) .== Ref(getfield(genes[1], :parent))) "Not all members of `genes` come from the same parent"
    if name == :index
        return [getfield(g, :index) for g in genes]
    elseif name == :parent
        return getfield(genes[1], name)
    else
        return GeneDataView(genes[1].parent, genes.index, name)
    end
end

function Base.setproperty!(gene::G, name::Symbol, x::T) where {G <: AbstractGene, T}
    if hasproperty(gene.parent.genedata, name)
        gene.parent.genedata[gene.index, name] = x
    else
        s = size(gene.parent.genedata, 1)
        gene.parent.genedata[!, name] = Vector{Union{Missing, T}}(missing, s)
        gene.parent.genedata[gene.index, name] = x
    end
    return x
end

Base.size(gv::GeneDataView) = size(gv.indices)
Base.getindex(gv::GeneDataView, i::Int) = getproperty(gv.parent.genes[gv.indices[i]], gv.property)
Base.getindex(gv::GeneDataView, I::AbstractArray) = GeneDataView(gv.parent, gv.indices[I], gv.property)
Base.setindex!(gv::GeneDataView, v, i::Int) = (Base.setproperty!(gv.parent.genes[i], gv.property, v))
Base.similar(gv::GeneDataView{Gene}) = GeneDataView(gv.parent, gv.indices, gv.property)
Base.copy(gv::GeneDataView{Gene}) = GeneDataView(gv.parent, gv.indices, gv.property)
Base.view(gv::GeneDataView{Gene}, I) = GeneDataView(gv.parent, gv.indices[I], gv.property)

function Base.fill!(gv::GeneDataView{Gene}, x)
    if hasproperty(gv.parent.genedata, gv.property)
        xT = convert(eltype(gv.parent.genedata[!, gv.property]), x)
    else
        xT = x
    end
    for I in gv.indices
        @inbounds gv[I] = x
    end
    gv
end


text:
using MetaImageFormat, FileIO, Images, AxisArrays
using Test

# TODO: change the next to lines to `load("example.mhd")` after registry
f = File(format"MetaImage", "example.mhd")
img = MetaImageFormat.load(f)

@test pixelspacing(img) == (25, 20, 40)
@test axisnames(img) == (:R, :A, :I)
@test eltype(img) == N0f16
@test size(img) == (8, 9, 5)

@info "Warning \"AnatomicalOrientation string unrecognized, got ???. Using xyz instead.\" is expected"
f2 = File(format"MetaImage", "example_dup_axes.mhd")
img2 = MetaImageFormat.load(f2)
@test parent(img) == parent(img2)
@test axisnames(img2) == (:x, :y, :z)


text:
function Base.show(io::IO, x::EXPR, d = 0, er = false)
    T = typof(x)
    c =  T === ErrorToken || er ? :red : :normal
    if isidentifier(x)
        printstyled(io, " "^d, valof(x), "  ", x.fullspan, "(", x.span, ")", color = :yellow)
        bindingof(x) !== nothing && printstyled(io, " $(bindingof(x).name)", bindingof(x).t === nothing ? "" : "::", color = :blue)
        refof(x) !== nothing && printstyled(io, " * ", color = :red)
        println(io)
    elseif isoperator(x)
        printstyled(io, " "^d, "OP: ", kindof(x), "  ", x.fullspan, "(", x.span, ")\n", color = c)
    elseif iskw(x)
        printstyled(io, " "^d, kindof(x), "  ", x.fullspan, "(", x.span, ")\n", color = :magenta)
    elseif ispunctuation(x)
        if kindof(x) == Tokens.LPAREN
            printstyled(io, " "^d, "(\n", color = c)
        elseif kindof(x) == Tokens.RPAREN
            printstyled(io, " "^d, ")\n", color = c)
        elseif kindof(x) == Tokens.LSQUARE
            printstyled(io, " "^d, "[\n", color = c)
        elseif kindof(x) == Tokens.RSQUARE
            printstyled(io, " "^d, "]\n", color = c)
        elseif kindof(x) == Tokens.COMMA
            printstyled(io, " "^d, ",\n", color = c)
        else
            printstyled(io, " "^d, "PUNC: ", kindof(x), "  ", x.fullspan, "(", x.span, ")\n", color = c)
        end
    elseif isliteral(x)
        printstyled(io, " "^d, "$(kindof(x)): ", valof(x), "  ", x.fullspan, "(", x.span, ")\n", color = c)
    else
        printstyled(io, " "^d, T, "  ", x.fullspan, "(", x.span, ")", color = c)
        scopeof(x) !== nothing && printstyled(io, " new scope", color = :green)
        bindingof(x) !== nothing && printstyled(io, " $(bindingof(x).name)", bindingof(x).t === nothing ? "" : "::", color = :blue)
        refof(x) !== nothing && printstyled(io, " * ", color = :red)
        println(io)
        x.args === nothing && return
        for a in x.args
            show(io, a, d + 1, er)
        end
    end
end

function Base.show(io::IO, scope::Scope)
    println(io, parentof(scope) === nothing ? "Root scope:" : "Scope:")
    println(io, scope.names isa Dict ? string("[", join(collect(keys(scope.names)), ","), "]") : "[]")
    println(io, scope.modules isa Dict ? string("[", join(collect(keys(scope.modules)), ","), "]") : "[]")
end


text:
function FindEntry(r)
    counter = 0
    while isborder(r, Nord)
        movementsbyfixmoves(r, West, counter)
        counter += 1
        if !isborder(r, Nord)
            break
        end
        movementsbyfixmoves(r, Ost, counter)
        counter += 1
    end
end

function movementsbyfixmoves(r, side, moves)
    while moves != 0
        move!(r, side)
        moves -= 1
    end
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Composition (Modia3D/Composition/_module.jl)
#

function build_tree!(scene::Scene, world::Object3D)::NOTHING
   tree                = scene.tree
   stack               = scene.stack
   allCollisionElements = scene.allCollisionElements
   empty!(tree)
   empty!(stack)
   empty!(scene.allVisuElements)
   empty!(allCollisionElements)

   assign_Visu_CutJoint_Dynamics!(scene, world, world)
   append!(stack, world.children)
   while length(stack) > 0
      obj = pop!(stack)
      assign_Visu_CutJoint_Dynamics!(scene, obj, world)
      if canCollide(obj)
        push!(allCollisionElements, obj)
      end
      push!(tree, obj)
      # Visit children of frame
      append!(stack, obj.children)
   end
   (length(allCollisionElements) > 1) ? (scene.collide = true) : (scene.collide = false)
   return nothing
end

insert_and_dedup!(v::Vector, x) = (splice!(v, searchsorted(v,x), x); v)

function addIndicesOfCutJointsToSuperObj(scene::Scene)
  tmp = collect(values(scene.noCPairsHelp))
  for i=1:length(tmp)
    if length(tmp[i]) == 2
      insert_and_dedup!(scene.noCPairs[minimum(tmp[i])], maximum(tmp[i]))
    else
      error("...from addIndicesOfCutJointsToSuperObj: problems with amount of cut joints")
    end
  end
  return nothing
end

function createAABB_noCPairs(scene::Scene, superObjsRow::SuperObjsRow)
  if length(superObjsRow.superObjCollision.superObj) > 0 && !isempty(superObjsRow.noCPair)
      push!(scene.noCPairs, superObjsRow.noCPair)
  else
    push!(scene.noCPairs, [0])
  end

  if length(superObjsRow.superObjCollision.superObj) > 0
    AABBrow = [Basics.BoundingBox() for i = 1:length(superObjsRow.superObjCollision.superObj)]
    push!(scene.AABB, AABBrow)
  else
    push!(scene.AABB, [])
  end
  return nothing
end

function changeParentToRootObj(newParent::Object3D, obj::Object3D)
  # Save elements of obj
  child_r_rel  = obj.r_rel
  child_R_rel  = obj.R_rel
  parent_r_rel = obj.parent.r_rel
  parent_R_rel = obj.parent.R_rel
  parent = obj.parent

  child_r_abs  = obj.r_abs
  child_R_abs  = obj.R_abs
  parent_r_abs = obj.parent.r_abs
  parent_R_abs = obj.parent.R_abs

  r_new = parent_r_rel + parent_R_rel' * child_r_rel
  R_new = child_R_rel * parent_R_rel
  obj.r_rel = r_new
  obj.R_rel = R_new

  # Reverse obj, so that newParent is the new parent
  obj.parent = newParent
  push!(newParent.children, obj)
  return nothing
end

# the indices of super objects, which can't collide, are stored in a list
function fillStackOrBuffer!(scene::Scene, superObj::SuperObjsRow, obj::Object3D, rootSuperObj::Object3D)::NOTHING
  n_children =  length(obj.children)
  help = fill(false, n_children)

  for i = 1:n_children
    child = obj.children[i]
    if isNotWorld(child)
      if isNotFixed(child)
        push!(scene.buffer, child)
        child.computeAcceleration = true
        if !child.joint.canCollide    # !isFree(child) &&  !( typeof(child.joint) <: Modia3D.AbstractPrismatic )
          push!(superObj.noCPair, length(scene.buffer))
        end
      else
        push!(scene.stack, child)
        assignAccVelo(scene.treeAccVelo, child)
        if !(obj == rootSuperObj)
          changeParentToRootObj(rootSuperObj, child)
          help[i] = true
  end; end; end; end

  if !isempty(obj.children)
    deleteat!(obj.children,help)
  end
  return nothing
end

# it builds the elements which may collide
# to reduce the amount of collision pairs, some assumptions are made:
#   elements which are rigidly attached, can't collide
#     these elements form together a super object
#   elements which are directly connected with a joint can't collide
#     these elements are excluded from the collision list
function build_superObjs!(scene::Scene, world::Object3D)::NOTHING
  if !scene.initSuperObj
  stack = scene.stack
  buffer = scene.buffer
  treeAccVelo  = scene.treeAccVelo
  empty!(stack)
  empty!(buffer)
  empty!(scene.allVisuElements)
  empty!(treeAccVelo)

  world.computeAcceleration = true
  push!(buffer, world)
  actPos = 1
  nPos   = 1

  hasOneCollisionSuperObj  = false
  hasMoreCollisionSuperObj = false

  while actPos <= nPos
    superObjsRow = SuperObjsRow()
    AABBrow      = Array{Basics.BoundingBox,1}()
    rootSuperObj = buffer[actPos]
    assign_Visu_CutJoint_Dynamics!(scene, rootSuperObj, world)
    if rootSuperObj != world
      assignAll(scene, superObjsRow, rootSuperObj, world, actPos)
      push!(treeAccVelo, rootSuperObj)
    end
    fillStackOrBuffer!(scene, superObjsRow, rootSuperObj, rootSuperObj)

    while length(stack) > 0
      frameChild = pop!(stack)
      assign_Visu_CutJoint_Dynamics!(scene, frameChild, world)
      assignAll(scene, superObjsRow, frameChild, world, actPos)
      fillStackOrBuffer!(scene, superObjsRow, frameChild, rootSuperObj)
    end

    if length(superObjsRow.superObjCollision.superObj) > 0 && hasOneCollisionSuperObj == true
      hasMoreCollisionSuperObj = true
    elseif length(superObjsRow.superObjCollision.superObj) > 0
      hasOneCollisionSuperObj = true
    end

    createAABB_noCPairs(scene, superObjsRow)
    push!(scene.superObjs, superObjsRow)
    nPos = length(buffer)
    actPos += 1
  end
  addIndicesOfCutJointsToSuperObj(scene)

#=
  println("superObjRow.superObjCollision.superObj")
  for superObjRow in scene.superObjs
  println("[")
  for a in superObjRow.superObjCollision.superObj
    println(ModiaMath.fullName(a))
  end
  println("]")
  println(" ")
  end
=#

  hasMoreCollisionSuperObj ? (scene.collide = true) : (scene.collide = false)
  scene.initSuperObj = true
  end
  return nothing
end

getAssembly(assembly::Modia3D.AbstractAssembly) = assembly
getAssembly(dummy)                              = nothing

function build_SignalObject3DConnections!(assembly::Modia3D.AbstractAssembly)
  scene     = assembly._internal.scene
  uniqueSignals       = scene.uniqueSignals
  uniqueForceTorques  = scene.uniqueForceTorques
  potentialVarInput   = scene.potentialVarInput
  potentialVarOutput  = scene.potentialVarOutput
  flowVarInput        = scene.flowVarInput
  flowVarOutput       = scene.flowVarOutput

  empty!(uniqueSignals)
  empty!(uniqueForceTorques)
  empty!(potentialVarInput)
  empty!(potentialVarOutput)
  empty!(flowVarInput)
  empty!(flowVarOutput)

  assemblyStack = Modia3D.AbstractAssembly[]
  push!(assemblyStack, assembly)
  while length(assemblyStack) > 0
    assembly = pop!(assemblyStack)

    # for checking if there are any sub assemblies
    for v in fieldnames(typeof(assembly))
      subAssembly = getAssembly(getfield(assembly,v))
      if typeof(subAssembly) != NOTHING && typeof(subAssembly) != Modia3D.Composition.Part
        push!(assemblyStack, subAssembly)
    end; end

    # for signals
    for signal in assembly._internal.uniqueSignals
      if !in(signal, uniqueSignals)
        push!(uniqueSignals, signal)
    end; end

    # for forces and torques
    for forceTorque in assembly._internal.uniqueForceTorques
      if !in(forceTorque, uniqueForceTorques)
        push!(uniqueForceTorques, forceTorque)
    end; end

    append!(potentialVarInput, assembly._internal.potentialVarInput)
    append!(potentialVarOutput, assembly._internal.potentialVarOutput)
    append!(flowVarInput, assembly._internal.flowVarInput)
    append!(flowVarOutput, assembly._internal.flowVarOutput)
  end
  return nothing
end
build_SignalObject3DConnections!(dummy) = nothing

"""
    visualizeAssembly!(assembly::Modia3D.AbstractAssembly)

Visualize the `assembly` defined with macro [`Modia3D.@assembly`](@ref)
in its initial configuration (but without simulating it).
"""
function visualizeAssembly!(assembly::Modia3D.AbstractAssembly)::NOTHING
   Modia3D.initAnalysis!(assembly)
   Modia3D.updatePosition!(assembly)
   Modia3D.visualize!(assembly,0.0)
   Modia3D.visualize!(assembly,1.0)
   Modia3D.closeAnalysis!(assembly)
   return nothing
end

function visualizeWorld!(world::Object3D; sceneOptions = SceneOptions())::NOTHING
   scene = Scene(sceneOptions)
   scene.analysis = ModiaMath.KinematicAnalysis
   initAnalysis!(world, scene)
   updatePosition!(world)
   visualize!(Modia3D.renderer[1], 0.0)
   visualize!(Modia3D.renderer[1], 1.0)
   closeAnalysis!(scene)
   return nothing
end

function chooseAndBuildUpTree(world::Object3D, scene::Scene)
  # Build tree for optimized structure or standard structure
  # collision handling is only available for optimized structure
  if scene.options.useOptimizedStructure
     build_superObjs!(scene, world)
     if scene.options.enableContactDetection && scene.collide
        initializeContactDetection!(world, scene)
        append!(scene.allVisuElements, world.contactVisuObj1)
        append!(scene.allVisuElements, world.contactVisuObj2)
        append!(scene.allVisuElements, world.supportVisuObj1A)
        append!(scene.allVisuElements, world.supportVisuObj1B)
        append!(scene.allVisuElements, world.supportVisuObj1C)
        append!(scene.allVisuElements, world.supportVisuObj2A)
        append!(scene.allVisuElements, world.supportVisuObj2B)
        append!(scene.allVisuElements, world.supportVisuObj2C)
     end
     initializeMassComputation!(scene)
  else
     build_tree!(scene, world)
     if scene.options.enableContactDetection
        error("Collision handling is only possible with the optimized structure. Please set useOptimizedStructure = true in Modia3D.SceneOptions.")
     end
  end
  if scene.visualize
     initializeVisualization(Modia3D.renderer[1], scene.allVisuElements)
  end
  scene.initAnalysis = true
  return nothing
end

function initAnalysis!(world::Object3D, scene::Scene)
   # Initialize spanning tree and visualization (if visualization desired and visual elements present)
  println("initAnalysis!(world::Object3D, scene::Scene)")
  chooseAndBuildUpTree(world, scene)
  return nothing
end

function initAnalysis!(assembly::Modia3D.AbstractAssembly;
                       analysis::ModiaMath.AnalysisType=ModiaMath.KinematicAnalysis)
   # Add visual elements to scene.allVisuElements
   if typeof(assembly._internal.referenceObject3D) == NOTHING
      error("\nError message from Modia3D.initAnalysis!(..):\n",
            typeof(assembly), " has no reference frame.")
   end
   #println("initAnalysis!(assembly::Modia3D.AbstractAssembly;                        analysis::ModiaMath.AnalysisType=ModiaMath.KinematicAnalysis)")
   # Construct Scene(..) object
   world = assembly._internal.referenceObject3D
   so    = assembly._internal.sceneOptions
   sceneOptions::SceneOptions = typeof(so) == NOTHING ? SceneOptions() : so
   scene = Scene(sceneOptions)
   scene.analysis = analysis
   assembly._internal.scene = scene

   chooseAndBuildUpTree(world, scene)
   return nothing
end

# Error Messages
function assertInitAnalysis(assembly::Modia3D.AbstractAssembly, functionName::String)
  if !assembly._internal.scene.initAnalysis
    error("\nError message from ", functionName, ":\n",
          "Call Modia3D.initAnalysis!(..) first.")
  end
end

function errorMessageCollision(functionName::String)
  error("\nError message from ", functionName, ":\n",
        "There are too less solids which can collide.")
end

function visualize!(assembly::Modia3D.AbstractAssembly, time::Float64)
  assertInitAnalysis(assembly, "Modia3D.visualize!(..)")
  scene = assembly._internal.scene
  if scene.visualize
    visualize!(Modia3D.renderer[1], time)
  end
end

function visualize!(scene::Scene, time::Float64)
  if scene.visualize
    visualize!(Modia3D.renderer[1], time)
  end
end

function getDistances!(assembly::Modia3D.AbstractAssembly)
  assertInitAnalysis(assembly, "Modia3D.getDistances!(..)")
  if assembly._internal.scene.collide
    ch = assembly._internal.scene.options.contactDetection
    world = assembly._internal.referenceObject3D
    getDistances!(ch, world)
  else
    errorMessageCollision("Modia3D.getDistances!(..)")
  end
end

function setComputationFlag(assembly::Modia3D.AbstractAssembly)
  assertInitAnalysis(assembly, "Modia3D.setComputationFlag!(..)")
  if assembly._internal.scene.collide
    ch = assembly._internal.scene.options.contactDetection
    setComputationFlag(ch)
  else
    errorMessageCollision("Modia3D.setComputationFlag!(..)")
  end
end

function selectContactPairsWithEvent!(assembly::Modia3D.AbstractAssembly)
  assertInitAnalysis(assembly, "Modia3D.selectContactPairsWithEvent!(..)")
  if assembly._internal.scene.collide
    ch = assembly._internal.scene.options.contactDetection
    world = assembly._internal.referenceObject3D
    selectContactPairsWithEvent!(nothing, ch, world)
  else
    errorMessageCollision("Modia3D.selectContactPairsWithEvent!(..)")
  end
end

function selectContactPairsNoEvent!(assembly::Modia3D.AbstractAssembly)
  assertInitAnalysis(assembly, "Modia3D.selectContactPairsNoEvent!(..)")
  if assembly._internal.scene.collide
    ch = assembly._internal.scene.options.contactDetection
    world = assembly._internal.referenceObject3D
    selectContactPairsNoEvent!(nothing, ch, world)
  else
    errorMessageCollision("Modia3D.selectContactPairsNoEvent!(..)")
  end
end

function performAnalysis!(assembly::Modia3D.AbstractAssembly, time::Float64)
   assertInitAnalysis(assembly, "performAnalysis!(..)")
   scene = assembly._internal.scene

   for frame in scene.tree
      parent = frame.parent
      if !(frame.r_rel ≡ ModiaMath.ZeroVector3D)
         frame.r_abs = parent.r_abs + parent.R_abs'*frame.r_rel
      end
      if !(frame.R_rel ≡ ModiaMath.NullRotation)
         frame.R_abs = frame.R_rel*parent.R_abs
      end
   end

   if scene.visualize
      visualize!(Modia3D.renderer[1], time)
   end
end

function closeContactDetection!(assembly::Modia3D.AbstractAssembly)
  assertInitAnalysis(assembly, "Modia3D.closeContactDetection!(..)")
  if assembly._internal.scene.collide
    ch = assembly._internal.scene.options.contactDetection
    closeContactDetection!(ch)
  else
    errorMessageCollision("Modia3D.closeContactDetection!(..)")
  end
end

function closeAnalysis!(scene::Scene)
     # Close Visualisation
     closeVisualization(Modia3D.renderer[1])
     empty!(scene.allVisuElements)
     scene.visualize = false
     # Close Collision detection
     if scene.collide
        closeContactDetection!(scene.options.contactDetection)
     end
     empty!(scene.stack)
     empty!(scene.buffer)
     empty!(scene.superObjs)
     empty!(scene.noCPairs)
     scene.collide = false

     scene.initSuperObj = false

     # Close signals and Forces and Torques
     empty!(scene.uniqueSignals)
     empty!(scene.uniqueForceTorques)
     empty!(scene.potentialVarInput)
     empty!(scene.potentialVarOutput)
     empty!(scene.flowVarInput)
     empty!(scene.flowVarOutput)

     # Close Spanning tree
     empty!(scene.tree)

     # Close Analysis
     scene.initAnalysis = false
end

closeAnalysis!(assembly::Modia3D.AbstractAssembly) = typeof(assembly._internal.scene) != NOTHING ? closeAnalysis!(assembly._internal.scene) : nothing


text:
function solve()
    s = 0
    n = readline()
    for c in n
        s += parse(Int, c)
    end
    s % 9 == 0 ? "Yes" : "No"
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
function codegen(::Type{Val{:iterate}}, ::Type{MH}, job::BasicMCJob)
  result::Expr
  update = []
  noupdate = []
  body = []

  vform = variate_form(job.pstate)
  if vform != Univariate && vform != Multivariate
    error("Only univariate or multivariate parameter states allowed in MH code generation")
  end

  if job.tuner.verbose
    push!(body, :(_job.sstate.tune.proposed += 1))
  end

  push!(body, :(_job.sstate.proposal = _job.sampler.setproposal(_job.pstate)))

  push!(body, :(_job.sstate.pstate.value =  rand(_job.sstate.proposal)))

  push!(body, :(_job.parameter.logtarget!(_job.sstate.pstate)))

  push!(body, :(_job.sstate.ratio = _job.sstate.pstate.logtarget-_job.pstate.logtarget))
  if !(job.sampler.symmetric)
    push!(body, :(_job.sstate.ratio -= logpdf(_job.sstate.proposal, _job.sstate.pstate.value)))
    if !(job.sampler.normalised)
      push!(body, :(_job.sstate.ratio -= lognormalise(_job.sstate.proposal)))
    end
    push!(body, :(_job.sstate.proposal = _job.sampler.setproposal(_job.sstate.pstate)))
    push!(body, :(_job.sstate.ratio += logpdf(_job.sstate.proposal, _job.pstate.value)))
    if !(job.sampler.normalised)
      push!(body, :(_job.sstate.ratio += lognormalise(_job.sstate.proposal)))
    end
  end

  if vform == Univariate
    push!(update, :(_job.pstate.value = _job.sstate.pstate.value))
  elseif vform == Multivariate
    push!(update, :(_job.pstate.value = copy(_job.sstate.pstate.value)))
  end
  push!(update, :(_job.pstate.logtarget = _job.sstate.pstate.logtarget))
  if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
    push!(update, :(_job.pstate.loglikelihood = _job.sstate.pstate.loglikelihood))
  end
  if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
    push!(update, :(_job.pstate.logprior = _job.sstate.pstate.logprior))
  end
  if in(:accept, job.outopts[:diagnostics])
    push!(update, :(_job.pstate.diagnosticvalues[1] = true))
    push!(noupdate, :(_job.pstate.diagnosticvalues[1] = false))
  end
  if job.tuner.verbose
    push!(update, :(_job.sstate.tune.accepted += 1))
  end

  push!(body, Expr(:if, :(_job.sstate.ratio > 0 || (_job.sstate.ratio > log(rand()))), Expr(:block, update...), noupdate...))

  if job.tuner.verbose
    fmt_iter = format_iteration(ndigits(job.range.burnin))
    fmt_perc = format_percentage()

    push!(body, :(
      if _job.sstate.tune.totproposed <= _job.range.burnin && mod(_job.sstate.tune.proposed, _job.tuner.period) == 0
        rate!(_job.sstate.tune)
        println(
          "Burnin iteration ",
          $(fmt_iter)(_job.sstate.tune.totproposed),
          " of ",
          _job.range.burnin,
          ": ",
          $(fmt_perc)(100*_job.sstate.tune.rate),
          " % acceptance rate"
        )
        reset_burnin!(_job.sstate.tune)
      end
    ))
  end

  if !job.plain
    push!(body, :(produce()))
  end

  @gensym _iterate

  result = quote
    function $_iterate(_job::BasicMCJob)
      $(body...)
    end
  end

  result
end


text:
# Merging of networks and stations

"""
    append!(sxml1, sxml2) -> sxml1

Add all the networks in `sxml2` to `sxml1`, duplicating any networks, stations or
channels that appear in both, even if they are identical.

Note that the `Network`s in the `network` field of `sxml2` will only be
**referred to** in `sxml1`s `network` field.  This means that changes to
`sxml2.network` will also be reflected in `sxml1.network`, even after the
`append!` call.  To avoid this, use `append!(sxml1, deepcopy(sxml2))`.

To return a copy of `sxml1` with the networks of `sxml2` added, whilst leaving
the original `sxml1` unmodified, use `append!(deepcopy(sxml1), sxml2)`.

See also:
[`merge!`](@ref Base.merge!(::StationXML.FDSNStationXML, ::StationXML.FDSNStationXML)),
[`merge`](@ref Base.merge(::StationXML.FDSNStationXML, ::StationXML.FDSNStationXML))
"""
function Base.append!(sxml1::FDSNStationXML, sxml2::FDSNStationXML)
    append!(sxml1.network, sxml2.network)
    sxml1
end

"""
    merge!(sxml1, sxml2[, ...]; warn=true) -> sxml1

Merge two or more [`FDSNStationXML`](@ref) objects into `sxml1`, adding unique
networks, stations and channels so that all are contained in a single object.
Networks are distinguished by their codes; stations and channels are distinguished
by their codes and start and end dates; channels are also distinguished by
their location codes.  Networks, stations and channels which are
duplicated between `sxml1`, `sxml2` (and `sxml3...`) are not duplicated in the
final merged `FDSNStationXML` object.

The implementation always takes the `FDSNStationXML` base information from the first
object.  That is, everything is determined by `sxml1`, apart from the `network` field,
meaning e.g. `source`, `sender`, `created`, etc.

If `warn` is `true` (the default), then warnings are printed when merging
networks or stations with the same code but otherwise different data (e.g.,
channels active at the same time with the same codes, but different azimuths).
If you want to duplicate all networks, stations or channels, use
[`append!`](@ref Base.append!(::StationXML.FDSNStationXML, ::StationXML.FDSNStationXML))
instead.

# Example
Find stations using a web request; `FDSNStationXML` objects are stored in the
stations' `.meta` structure.  We can then merge the station information together.
```
julia> using SeisRequests

julia> stas = get_stations(code="XM.C*..HHZ", starttime="2012-01-01", endtime="2014-01-01")
8-element Vector{Seis.GeogStation{Float64}}:
 Station: XM.C01E..HHZ, lon: 38.385052, lat: 7.22327, dep: 0.0, elev: 1831.0, azi: 0.0, inc: 0.0, meta: 5 keys
 Station: XM.C02E..HHZ, lon: 38.43224, lat: 7.14882, dep: 0.0, elev: 1722.0, azi: 0.0, inc: 0.0, meta: 5 keys
 Station: XM.C03E..HHZ, lon: 38.484409, lat: 7.19804, dep: 0.0, elev: 1759.0, azi: 0.0, inc: 0.0, meta: 5 keys
 Station: XM.C04E..HHZ, lon: 38.476398, lat: 7.21787, dep: 0.0, elev: 1782.0, azi: 0.0, inc: 0.0, meta: 5 keys
 Station: XM.C04E..HHZ, lon: 38.476398, lat: 7.21787, dep: 0.0, elev: 1782.0, azi: 0.0, inc: 0.0, meta: 5 keys
 Station: XM.C05E..HHZ, lon: 38.30228, lat: 7.25978, dep: 0.0, elev: 1946.0, azi: 0.0, inc: 0.0, meta: 5 keys
 Station: XM.C06E..HHZ, lon: 38.372139, lat: 7.13389, dep: 0.0, elev: 1775.0, azi: 0.0, inc: 0.0, meta: 5 keys
 Station: XM.C07E..HHZ, lon: 38.444, lat: 7.26028, dep: 0.0, elev: 1710.0, azi: 0.0, inc: 0.0, meta: 5 keys

julia> sxml = merge(stas.meta.stationxml...) # Merge all separate FDSNStationXML objects
StationXML.FDSNStationXML
  source: String "IRIS-DMC"
  sender: String "IRIS-DMC"
  module_name: String "IRIS WEB SERVICE: fdsnws-station | version: 1.1.47"
  module_uri: String "http://service.iris.edu/fdsnws/station/1/query?starttime=2012-01-01T00:00:00&endtime=2014-01-01T00:00:00&network=XM&station=C*&location=&channel=HHZ&level=channel&nodata=204"
  created: Dates.DateTime
  network: Array{StationXML.Network}((1,))
  schema_version: String "1.1"

julia> channel_codes(sxml) # List of all channel codes in the StationXML
8-element Vector{String}:
 "XM.C01E..HHZ"
 "XM.C02E..HHZ"
 "XM.C03E..HHZ"
 "XM.C04E..HHZ"
 "XM.C04E..HHZ"
 "XM.C05E..HHZ"
 "XM.C06E..HHZ"
 "XM.C07E..HHZ"
```

See also:
[`append!`](@ref Base.append!(::StationXML.FDSNStationXML, ::StationXML.FDSNStationXML))
[`merge`](@ref Base.merge(::StationXML.FDSNStationXML, ::StationXML.FDSNStationXML))
"""
function Base.merge!(sxml1::FDSNStationXML, sxml2::FDSNStationXML; warn=true)
    _merge!(nothing, sxml1, sxml2, warn)
    sxml1
end

function Base.merge!(sxml1::FDSNStationXML, sxml2::FDSNStationXML,
        sxmls::Vararg{FDSNStationXML}; warn=true)
    foldl((x, y) -> merge!(x, y; warn=warn), (sxml1, sxml2, sxmls...))
end

"""
    merge(sxml1, sxml2[, ...]; warn=true) -> sxml

Merge two or more [`FDSNStationXML`](@ref) objects together, adding unique
networks, stations and channels so that all are contained in a single object.

See [`merge!`](@ref Base.merge!(::FDSNStationXML, ::FDSNStationXML)) for
more details of the criteria by which networks, stations and channels
are merged.

See also:
[`append!`](@ref Base.append!(::StationXML.FDSNStationXML, ::StationXML.FDSNStationXML))
"""
Base.merge(sxml1::FDSNStationXML, sxml2::FDSNStationXML; warn=true) =
    merge!(deepcopy(sxml1), sxml2; warn=warn)

function Base.merge(sxml1::FDSNStationXML, sxml2::FDSNStationXML,
        sxmls::Vararg{FDSNStationXML}; warn=true)
    merge!(deepcopy(sxml1), sxml2, sxmls...; warn=warn)
end

_items_for_merging(x::FDSNStationXML, y::FDSNStationXML) = x.network, y.network
_items_for_merging(x::Network, y::Network) = x.station, y.station

"""
Merge two collections `x1` and `x2` together into `x1`.  `x1` is a field of `parent`.
"""
function _merge!(parent, x1, x2, warn)
    y1s, y2s = _items_for_merging(x1, x2)
    # Find network or station codes
    codes1 = (y.code for y in y1s)
    codes2 = (y.code for y in y2s)
    # Items with the same code but which do not overlap in time are
    # considered to be separate (e.g., reusing of temporary network codes;
    # stations which move around in time) and so we don't merge stations or
    # channels.
    # Use `Dict{Int,eltype(y1s)}` rather than `Set{eltype(y1s)}` since this
    # only requires hashing an `Int` rather than a potentially complex type
    # and is therefore be faster.
    different_y2s = Dict{Int,eltype(y1s)}()
    for code in intersect(codes1, codes2)
        # Only iterate over all matching codes to avoid searching through every
        # network/station for both `x1` and `x2`
        i1s = findall(x -> x.code == code, y1s)
        i2s = findall(x -> x.code == code, y2s)
        for i1 in i1s
            y1 = y1s[i1]
            for i2 in i2s
                y2 = y2s[i2]
                # No need to copy identical networks/stations
                if y1 == y2
                    continue
                end
                # Assume that networks/stations with overlapping times are really
                # the same, but with different contents
                if _time_ranges_overlap(y1, y2)
                    _merge!(x1, y1s[i1], y2s[i2], warn)
                else
                    different_y2s[i2] = y2
                end
            end
        end
    end
    append!(y1s, values(different_y2s))

    # Append the items only in x2 to x1
    for code in setdiff(codes2, codes1)
        i2 = findfirst(x -> x.code == code, y2s)
        # deepcopy because the collections contain arrays and it may be unexpected
        # that changing them later in x1 would affect elements of x2
        push!(y1s, deepcopy(y2s[i2]))
    end
    x1
end

"""
Merge two stations together, collecting unique channels into `sta1`.
"""
function _merge!(network::Network, sta1::Station, sta2::Station, warn)
    cha1, cha2 = sta1.channel, sta2.channel
    # Get set of location and channel codes
    codes1 = ((c.code, c.location_code) for c in cha1)
    codes2 = ((c.code, c.location_code) for c in cha2)
    # Channels with the same code but which do not overlap in time are
    # considered to be separate.
    # Use a Dict so we don't need to keep track of whether we have already
    different_c2_channels = Dict{Int,Channel}()
    for (code, loc) in intersect(codes1, codes2)
        i1s = findall(x -> x.code == code && x.location_code == loc, cha1)
        i2s = findall(x -> x.code == code && x.location_code == loc, cha2)
        @assert !(isempty(i1s) || isempty(i2s)) "unexpectedly empty set of channels"
        for i1 in i1s
            c1 = cha1[i1]
            for i2 in i2s
                c2 = cha2[i2]
                # Skip identical channels
                if c1 == c2
                    continue
                end
                # There is a time overlap, so these must be the same channel
                if _time_ranges_overlap(c1, c2)
                    if warn
                        chan_code = join((network.code, sta1.code, loc, code), '.')
                        timerange1 = string(c1.start_date, " - ", c1.end_date)
                        timerange2 = string(c2.start_date, " - ", c2.end_date)
                        @warn("""
                            two channels with code '$(chan_code)' overlap in time but are not the same; only keeping the first
                            Time ranges:
                                $(timerange1)
                                $(timerange2)
                            """)
                    end
                # These are different in time, so we need to add c2
                else
                    different_c2_channels[i2] = c2
                end
            end
        end
    end
    append!(cha1, values(different_c2_channels))

    # Append the channels only in cha2 to cha1
    for (code, loc) in setdiff(codes2, codes1)
        i2 = findfirst(x -> x.code == code && x.location_code == loc, cha2)
        push!(cha1, cha2[i2])
    end
    sta1
end

"""
    _time_ranges_overlap(x, y) -> ::Bool

Return `true` if the two networks, stations or channels `x` and `y` overlap
by more than one point in time.

Note that typically the end time of a period of channel operation is the same
as the start time of the next period of operation, hence even when there is no
intended overlap this creates a single point in time where the two ranges
overlap.  Therefore we require two points to overlap in order for a true
overlap to be decided.

Start and end dates are permitted to be `missing` (absent), in which case it
is assumed that a missing start date corresponds to a station always operating
into the past, and that a missing end date corresponds to a station always
operating into the future.  Hence two stations which both have missing end
dates but set start dates are assumed to overlap, and similarly for a pair
with different end dates but both with missing start dates
"""
function _time_ranges_overlap(x, y)
    b1, e1 = x.start_date, x.end_date
    b2, e2 = y.start_date, y.end_date
    nmissings = count(x -> x === missing, (b1, e1, b2, e2))

    # If all missing or only one non-missing, assume overlap
    if nmissings in (3, 4)
        return true
    elseif nmissings == 2
        # One station has no date information; assume overlap
        if (b1 === missing && e1 === missing) || (b2 === missing && e2 === missing)
            return true
        # Pairs of start or end times are missing; assume overlap
        elseif (b1 === missing && b2 === missing) || (e1 === missing && e2 === missing)
            return true
        # 'Cross terms' are missing, so only need to check one pair
        elseif (b1 === missing && e2 === missing && b2 >= e1) ||
                (b2 === missing && e1 === missing && b1 >= e2)
            return false
        # We have 'cross terms' but there is an overlap
        else
            return true
        end
    # Just convert the single missing value to something we can work with then
    elseif nmissings === 1
        b1 === missing && (b1 = typemin(b2))
        b2 === missing && (b2 = typemin(b1))
        e1 === missing && (e1 = typemax(e2))
        e2 === missing && (e2 = typemax(e1))
    end

    # No missing values now
    # x before y and no overlap
    if e1 <= b2
        return false
    # y before x and no overlap
    elseif e2 <= b1
        return false
    else
        return true
    end
end


text:
using ITensors

N = 4
nmps = 3
cutoff = 1e-8
s = siteinds("S=1/2", N)
ψs = [randomMPS(s; linkdims=2) for _ in 1:nmps]
ρs = [outer(ψ, ψ; cutoff) for ψ in ψs]
ρ = sum(ρs; cutoff)

ψs_full = prod.(ψs)
ρs_full = [ψ'dag(ψ) for ψ in ψs_full]
ρ_full = sum(ρs_full)
@show norm(ρ_full - prod(ρ))


text:
lower(str::AbstractString) = lowercase(strip(str))

function has_ssh_private_key(; env::AbstractDict=ENV)
    value = strip(get(env, PRIVATE_SSH_ENVVAR, ""))

    x1 = haskey(env, PRIVATE_SSH_ENVVAR)
    x2 = value isa AbstractString
    x3 = !isempty(value)
    x4 = value != "false"

    return x1 && x2 && x3 && x4
end

function api_retry(f::Function)
    delays = ExponentialBackOff(; n=10, max_delay=30.0)
    return retry(f; delays)()
end

function _max(x::Union{Nothing,Any}, y)
    if x === nothing
        return y
    else
        return max(x, y)
    end
end

function with_temp_dir(f::Function; kwargs...)
    original_dir = pwd()
    tmp_dir = mktempdir(; kwargs...)
    atexit(() -> rm(tmp_dir; force=true, recursive=true))

    cd(tmp_dir)
    result = f(tmp_dir)
    cd(original_dir)

    rm(tmp_dir; force=true, recursive=true)
    return result
end

function add_compat_section!(project::Dict)
    if !haskey(project, "compat")
        project["compat"] = Dict{Any,Any}()
    end

    return project
end

function get_random_string()
    randint = lpad(string(rand(UInt32)), 11, "0")
    return string(utc_to_string(now_localzone()), "-", randint)
end


text:
function space(::SiteType"Electron";
               conserve_qns = false,
               conserve_sz = conserve_qns,
               conserve_nf = conserve_qns,
               conserve_nfparity = conserve_qns,
               qnname_sz = "Sz",
               qnname_nf = "Nf",
               qnname_nfparity = "NfParity",
               # Deprecated
               conserve_parity = nothing)
  if !isnothing(conserve_parity)
    conserve_nfparity = conserve_parity
  end
  if conserve_sz && conserve_nf
    return [
       QN((qnname_nf,0,-1),(qnname_sz, 0)) => 1
       QN((qnname_nf,1,-1),(qnname_sz,+1)) => 1
       QN((qnname_nf,1,-1),(qnname_sz,-1)) => 1
       QN((qnname_nf,2,-1),(qnname_sz, 0)) => 1
      ]
  elseif conserve_nf
    return [
       QN(qnname_nf,0,-1) => 1
       QN(qnname_nf,1,-1) => 2
       QN(qnname_nf,2,-1) => 1
      ]
  elseif conserve_sz
    return [
       QN((qnname_sz, 0),(qnname_nfparity,0,-2)) => 1
       QN((qnname_sz,+1),(qnname_nfparity,1,-2)) => 1
       QN((qnname_sz,-1),(qnname_nfparity,1,-2)) => 1
       QN((qnname_sz, 0),(qnname_nfparity,0,-2)) => 1
      ]
  elseif conserve_nfparity
    return [
       QN(qnname_nfparity,0,-2) => 1
       QN(qnname_nfparity,1,-2) => 2
       QN(qnname_nfparity,0,-2) => 1
      ]
  end
  return 4
end

state(::SiteType"Electron",::StateName"Emp")  = 1
state(::SiteType"Electron",::StateName"Up")   = 2
state(::SiteType"Electron",::StateName"Dn")   = 3
state(::SiteType"Electron",::StateName"UpDn") = 4
state(st::SiteType"Electron",::StateName"0")    = state(st,StateName("Emp"))
state(st::SiteType"Electron",::StateName"↑")    = state(st,StateName("Up"))
state(st::SiteType"Electron",::StateName"↓")    = state(st,StateName("Dn"))
state(st::SiteType"Electron",::StateName"↑↓")   = state(st,StateName("UpDn"))

function op!(Op::ITensor,
             ::OpName"Nup",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>2,s=>2] = 1.0
  Op[s'=>4,s=>4] = 1.0
end

function op!(Op::ITensor,
             ::OpName"Ndn",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>3,s=>3] = 1.0
  Op[s'=>4,s=>4] = 1.0
end

function op!(Op::ITensor,
             ::OpName"Nupdn",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>4,s=>4] = 1.0
end

function op!(Op::ITensor,
             ::OpName"Ntot",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>2,s=>2] = 1.0
  Op[s'=>3,s=>3] = 1.0
  Op[s'=>4,s=>4] = 2.0
end

function op!(Op::ITensor,
             ::OpName"Cup",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>1,s=>2] = 1.0
  Op[s'=>3,s=>4] = 1.0
end

function op!(Op::ITensor,
             ::OpName"Cdagup",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>2,s=>1] = 1.0
  Op[s'=>4,s=>3] = 1.0
end

function op!(Op::ITensor,
             ::OpName"Cdn",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>1,s=>3] = 1.0
  Op[s'=>2,s=>4] = -1.0
end

function op!(Op::ITensor,
             ::OpName"Cdagdn",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>3,s=>1] = 1.0
  Op[s'=>4,s=>2] = -1.0
end

function op!(Op::ITensor,
             ::OpName"Aup",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>1,s=>2] = 1.0
  Op[s'=>3,s=>4] = 1.0
end

function op!(Op::ITensor,
             ::OpName"Adagup",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>2,s=>1] = 1.0
  Op[s'=>4,s=>3] = 1.0
end

function op!(Op::ITensor,
             ::OpName"Adn",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>1,s=>3] = 1.0
  Op[s'=>2,s=>4] = 1.0
end

function op!(Op::ITensor,
             ::OpName"Adagdn",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>3,s=>1] = 1.0
  Op[s'=>2,s=>4] = 1.0
end

function op!(Op::ITensor,
             ::OpName"F",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>1,s=>1] = +1.0
  Op[s'=>2,s=>2] = -1.0
  Op[s'=>3,s=>3] = -1.0
  Op[s'=>4,s=>4] = +1.0
end

function op!(Op::ITensor,
             ::OpName"Fup",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>1,s=>1] = +1.0
  Op[s'=>2,s=>2] = -1.0
  Op[s'=>3,s=>3] = +1.0
  Op[s'=>4,s=>4] = -1.0
end

function op!(Op::ITensor,
             ::OpName"Fdn",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>1,s=>1] = +1.0
  Op[s'=>2,s=>2] = +1.0
  Op[s'=>3,s=>3] = -1.0
  Op[s'=>4,s=>4] = -1.0
end

function op!(Op::ITensor,
             ::OpName"Sz",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>2,s=>2] = +0.5
  Op[s'=>3,s=>3] = -0.5
end

op!(Op::ITensor,
    ::OpName"Sᶻ",
    st::SiteType"Electron",
    s::Index) = op!(Op,OpName("Sz"),st,s)

function op!(Op::ITensor,
             ::OpName"Sx",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>2,s=>3] = 0.5
  Op[s'=>3,s=>2] = 0.5
end

op!(Op::ITensor,
    ::OpName"Sˣ",
    st::SiteType"Electron",
    s::Index) = op!(Op,OpName("Sx"),st,s)

function op!(Op::ITensor,
             ::OpName"S+",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>2,s=>3] = 1.0
end

op!(Op::ITensor,
    ::OpName"S⁺",
    st::SiteType"Electron",
    s::Index) = op!(Op,OpName("S+"),st,s)
op!(Op::ITensor,
    ::OpName"Sp",
    st::SiteType"Electron",
    s::Index) = op!(Op,OpName("S+"),st,s)
op!(Op::ITensor,
    ::OpName"Splus",
    st::SiteType"Electron",
    s::Index) = op!(Op,OpName("S+"),st,s)

function op!(Op::ITensor,
             ::OpName"S-",
             ::SiteType"Electron",
             s::Index)
  Op[s'=>3,s=>2] = 1.0
end

op!(Op::ITensor,
    ::OpName"S⁻",
    st::SiteType"Electron",
    s::Index) = op!(Op,OpName("S-"),st,s)
op!(Op::ITensor,
    ::OpName"Sm",
    st::SiteType"Electron",
    s::Index) = op!(Op,OpName("S-"),st,s)
op!(Op::ITensor,
    ::OpName"Sminus",
    st::SiteType"Electron",
    s::Index) = op!(Op,OpName("S-"),st,s)

has_fermion_string(::OpName"Cup", ::SiteType"Electron") = true
has_fermion_string(::OpName"Cdagup", ::SiteType"Electron") = true
has_fermion_string(::OpName"Cdn", ::SiteType"Electron") = true
has_fermion_string(::OpName"Cdagdn", ::SiteType"Electron") = true


text:
module CameraGeometry

# package code goes here
using Images

end # module


text:
using ChainRulesCore

@doc raw"
    einsum_grad(ixs, xs, iy, size_dict, cdy, i)

return the gradient of the result of evaluating the `EinCode` w.r.t
the `i`th tensor in `xs`. `cdy` is the result of applying the `EinCode`
to the `xs`.

# example
```jldoctest; setup = :(using OMEinsum)
julia> using OMEinsum: einsum_grad, get_size_dict

julia> a, b = rand(2,2), rand(2,2);

julia> c = einsum(EinCode((('i','j'),('j','k')), ('i','k')), (a,b));

julia> sd = get_size_dict((('i','j'),('j','k')), (a,b));

julia> einsum_grad((('i','j'),('j','k')), (a,b), ('i','k'), sd, c, 1) ≈ c * transpose(b)
true
```
"
function einsum_grad(ixs, xs, iy, size_dict, cdy, i)
    nixs = TupleTools.insertat(ixs, i, (iy,))
    nxs  = TupleTools.insertat( xs, i, (cdy,))
    niy = ixs[i]
    y = einsum(EinCode(nixs, niy), nxs, size_dict)
    try
        conj!(y)
    catch e
        y = conj(y)
    end
    typeof(y) == typeof(xs[i]) && return y
    xs[i] isa Array{<:Real} && return convert(typeof(xs[i]), real(y))
    convert(typeof(xs[i]), y)
end

function ChainRulesCore.rrule(::typeof(einsum), code::EinCode{ixs, iy}, xs::NTuple{N,T where T}, size_dict) where {N, ixs, iy}
    y = einsum(code, xs, size_dict)
    function einsum_pullback(dy)
        dxs = ChainRulesCore.@thunk ntuple(i -> einsum_grad(ixs, xs, iy, size_dict, map(conj, dy), i), N)
        return (NoTangent(), NoTangent(), dxs, NoTangent())
    end
    einsum_pullback(::NoTangent) = (NoTangent(), NoTangent(), NoTangent(), NoTangent())
    return y, einsum_pullback
end

function dynamic_einsum_grad(ixs, xs, iy, size_dict, cdy, i)
    nixs = TupleTools.insertat(ixs, i, (iy,))
    nxs  = TupleTools.insertat( xs, i, (cdy,))
    niy = ixs[i]
    y = dynamic_einsum(nixs, nxs, niy, size_dict)
    try
        conj!(y)
    catch e
        y = conj(y)
    end
    typeof(y) == typeof(xs[i]) && return y
    xs[i] isa Array{<:Real} && return convert(typeof(xs[i]), real(y))
    convert(typeof(xs[i]), y)
end

function ChainRulesCore.rrule(::typeof(dynamic_einsum), ixs, xs::NTuple{N,T where T}, iy, size_dict) where {N}
    y = dynamic_einsum(ixs, xs, iy, size_dict)
    function einsum_pullback(dy)
        dxs = ChainRulesCore.@thunk ntuple(i -> einsum_grad(ixs, xs, iy, size_dict, map(conj, dy), i), N)
        return (NoTangent(), NoTangent(), dxs, NoTangent(), NoTangent())
    end
    einsum_pullback(::NoTangent) = (NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent())
    return y, einsum_pullback
end