text:
type IIDMInput <: AbstractInput end

## Read inputs
function read_input(intput_type::T, instance_path::String) where T<:Type{IIDMInput}
  data = parse_file(instance_path)

  ## Cllecting bus information
  bus = Dict{String, Dict{String, Any}}()
  ind_bus = 1
  for substation in root(data)["substation"]
    substationname = attribute(substation, "id")
    busind = ind_bus
    ind_bus += 1

    voltagelvl = substation["voltageLevel"][1]
    length(substation["voltageLevel"])==1 || warn("Substation $substationname has several voltageLevel, accounting for $(attribute(voltagelvl, "id")) only")

    busname = attribute(voltagelvl, "id")
    length(voltagelvl["busBreakerTopology"][1]["bus"]) == 1 || warn("Substation has several buses, accounting for $busname only..")
    basekV = parse(attribute(voltagelvl, "nominalV"))
    minV = parse(attribute(voltagelvl, "lowVoltageLimit"))
    maxV = parse(attribute(voltagelvl, "highVoltageLimit"))

    busBreakerTopology = voltagelvl["busBreakerTopology"][1]["bus"][1]
    V, θ = parse(attribute(busBreakerTopology, "v")), parse(attribute(busBreakerTopology, "angle"))*π/180

    bus[busname] = Dict{String, Any}()
    bus[busname]["Volt"] = IIDMVolt(busname, busind, substationname, basekV, minV, maxV)
    bus[busname]["MeasureVolt"] = IIDMVoltMeasure(V*exp(im*θ), NaN)

    for cur_load in voltagelvl["load"]
      load_id = attribute(cur_load, "id")
      load0 = parse(attribute(cur_load, "p0")) + im*parse(attribute(cur_load, "q0"))
      load = parse(attribute(cur_load, "p")) + im*parse(attribute(cur_load, "q"))

      abs(load-load0) < 1e-5 || warn("$busname $load_id : load - load0 = $(load-load0)")

      bus[busname]["LOAD_$load_id"] = IIDMLoad(load_id, load)
    end

    for cur_gen in voltagelvl["generator"]
      gen_id = attribute(cur_gen, "id")
      minP = parse(attribute(cur_gen, "minP"))
      maxP = parse(attribute(cur_gen, "maxP"))
      targetV = parse(attribute(cur_gen, "targetV"))
      volt_regulator_on = parse(attribute(cur_gen, "voltageRegulatorOn"))

      attribute(cur_gen["property"][1], "value") == "FUEL" || warn("generator $gen_id of $busname not of expected type FUEL.")
      capabilityCurve = cur_gen["reactiveCapabilityCurve"][1]
      minP = parse(attribute(pt, "p"))
      minQminP, maxQminP = parse(attribute(pt, "minQ")), parse(attribute(pt, "maxQ"))
      maxP = parse(attribute(pt, "p"))
      minQmaxP, maxQmaxP = parse(attribute(pt, "minQ")), parse(attribute(pt, "maxQ"))
      S_bounds = [(minP+im*minQminP, minP+im*maxQminP), (maxP+im*minQmaxP, maxP+im*maxQmaxP)]

      bus[busname]["GEN_$gen_id"] = IIDMGenerator(gen_id, minP, maxP, targetV, volt_regulator_on, S_bounds)
    end
  end

  ## Collecting link information
  link = Dict{Link, Dict{String, Any}}()
  for line in root(data)["line"]
    orig = attribute(line, "voltageLevelId1")
    dest = attribute(line, "voltageLevelId2")
    linkname = Link(orig, dest)
    lineid = attribute(line, "id")
    r = parse(Float64, attribute(line, "r"))
    x = parse(Float64, attribute(line, "x"))
    G1 = parse(Float64, attribute(line, "g1")) + im*parse(Float64, attribute(line, "b1"))
    G2 = parse(Float64, attribute(line, "g2")) + im*parse(Float64, attribute(line, "b2"))
    S1 = parse(Float64, attribute(line, "p1")) + im*parse(Float64, attribute(line, "q1"))
    S2 = parse(Float64, attribute(line, "p2")) + im*parse(Float64, attribute(line, "q2"))
    basekV_orig = bus["$orig"]["Volt"].nomV
    basekV_dest = bus["$dest"]["Volt"].nomV
    lim = parse(Float64, attribute(line["currentLimits1"][1], "permanentLimit"))

    if !haskey(link, linkname)
      link[linkname] = Dict{String, Any}()
    end
    link[linkname][lineid] = IIDMLine_π(r, x, G1, G2, basekV_orig, basekV_dest, lim)
    link[linkname]["Measure_$lineid"] = IIDMLine_πMeasure(S1, S2)
  end

  ds = DataSource(bus, link)

  ## Building GridStructure
  node_linksin = node_linksout = Dict{String, Set{Link}}()
  node_vars = Dict{String, Dict{String, Variable}}()
  link_vars = Dict{Link, Dict{String, Variable}}()
  gs = GridStructure("BaseCase", node_linksin, node_linksout)

  node_formulations = Dict{String, Dict{Tuple{Type, String}, Symbol}}()
  link_formulations = Dict{Link, Dict{Tuple{Type, String}, Symbol}}()
  mp = MathematicalProgramming(node_formulations, link_formulations, node_vars, link_vars)
  return OPFProblems("BaseCase"=>Scenario(ds, gs, mp))
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
using Interp1d
using Test

t = time()

@testset "Interp1d.jl" begin
    # Write your tests here.

    for file in filter(f -> endswith(f, ".jl"), readdir(@__DIR__))
        file in ["runtests.jl"] && continue

        @testset "$(file)" begin
            t = time()
            include(file)
            println("$(file) took $(round(time() - t; digits = 1)) seconds.")
        end
    end
end

ext: jl


text:
show(io::IO, g::Graph) = print(io, "{$(g.nn),$(g.na)} Graph")
show(io::IO, a::Arc) = print(io, "Arc $(a.u) -$(a.key)-> $(a.v)")
show(io::IO, ait::ArcIter) = print(io, "ArcIter $(ait.len)")
function show(io::IO, path::Path)
    if !isempty(path)
        res = "Path {$(path.nodelist[1])}" # first
        for k in eachindex(path.keylist)
            res *= " -$(path.keylist[k])-> {$(path.nodelist[k + 1])}"
        end
        print(io, res)
    else
        res = "Empty Path"
        print(io, res)
    end
end


text:
Δ=ToeplitzOperator([1.],[0.,1.])
@manipulate for n=2:10
    K=SymTriOperator(rand(n),rand(n-1))
    eigvals(Δ+K)|>ApproxFun.plot
end


text:
module JuliaLS

if VERSION < v"1.0.0"
    error("Julia language server only works with julia 1.0.0+")
end

include("main.jl")

function julia_main()::Cint
    try
        main(ARGS)
    catch
        Base.display_error(err, catch_backtrace())
        return 1
    end
    return 0
end

end # module

max_stars_repo_name: wangl-cc/coc-julia

max_issues_repo_name: wangl-cc/coc-julia

max_forks_repo_name: wangl-cc/coc-julia


text:
using Morsel;
using SQLite;

db = SQLiteDB("quotes.db");
try
  res = query(db,"select count(*) from quotes");
catch
  error("Can't find quotes table");
end
nq = res[1][1];

@assert nq > 0;

htmlH = """<!DOCTYPE html>

<title>Julia Quotes Server<</title>

""";

htmlF= "";

app = Morsel.app();
route(app, POST | PUT, "/") do req, res
 htmlB = "
Posting data to the server is not allowed
";
 string(htmlH,htmlB,htmlF);
end

get(app, "/") do req, res
  qid = rand(1:nq);
  sql = "select q.author,c.catname,q.quoname from quotes q ";
  sql *= "join categories c on q.cid = c.id and q.id = $qid";
  res = query(db,sql);
  author  = res[1][1];
  catname = res[2][1];
  quotext = res[3][1];
  htmlB = "
$quotext
$author ($catname)
";
  string(htmlH,htmlB,htmlF);
end

max_stars_repo_path: Module 3/Chapter09/quotes-srv.jl

max_issues_repo_path: Module 3/Chapter09/quotes-srv.jl

max_forks_repo_path: Module 3/Chapter09/quotes-srv.jl


text:
using SparseArrays: SparseMatrixCSC, sparse
using LinearAlgebra: mul!
using Combinatorics: permutations
using LinearAlgebra: rank, svd, Diagonal

"""
`struct SymmetricBasis`

### Constructors

Option 1: pass a `OneParticleBasis`
```julia
SymmetricBasis(φ, symgrp, basis1p, Bsel)
SymmetricBasis(φ, basis1p, Bsel)   # uses default symgrp = O3()
```
will first construct a `PIBasis` from these inputs and then call the second
constructor.

Option 1: pass a `PIBasis`
```julia
SymmetricBasis(φ, symgrp, pibasis)
SymmetricBasis(φ, pibasis)
```
If the PIbasis is already available, this directly constructs a
resulting SymmetricBasis; all possible permutation-invariant basis functions
will be symmetrised and then reduced to a basis (rather than spanning set)
"""
mutable struct SymmetricBasis{PIB, PROP, SYM, REAL} <: ACEBasis
   pibasis::PIB
   A2Bmap::SparseMatrixCSC{PROP, Int}
   symgrp::SYM
   real::REAL
end

Base.length(basis::SymmetricBasis{PIB, PROP}) where {PIB, PROP} =
      size(basis.A2Bmap, 1)

# -------- FIO

==(B1::SymmetricBasis, B2::SymmetricBasis) =
      ( (B1.pibasis == B2.pibasis) &&
        (B1.A2Bmap == B2.A2Bmap) &&
        (B1.real == B2.real) )

write_dict(B::SymmetricBasis{PIB, PROP}) where {PIB, PROP} =
      Dict( "__id__" => "ACE_SymmetricBasis",
            "pibasis" => write_dict(B.pibasis),
            "A2Bmap" => write_dict(B.A2Bmap),
            "symgrp" => write_dict(B.symgrp),
            "isreal" => (B.real == Base.real) )

read_dict(::Val{:ACE_SymmetricBasis}, D::Dict) =
      SymmetricBasis(read_dict(D["pibasis"]),
                     read_dict(D["A2Bmap"]),
                     read_dict(D["symgrp"]),
                     (D["isreal"] ? Base.real : Base.identity) )
# --------

SymmetricBasis(φ::AbstractProperty,
               basis1p::OneParticleBasis,
               Bsel::AbstractBasisSelector;
               kwargs...) =
      SymmetricBasis(φ, basis1p, O3(), Bsel; kwargs...)

SymmetricBasis(φ::AbstractProperty, pibasis; kwargs...) =
      SymmetricBasis(φ, O3(), pibasis; kwargs...)

SymmetricBasis(φ::AbstractProperty,
               basis1p::OneParticleBasis,
               symgrp::SymmetryGroup,
               Bsel::AbstractBasisSelector;
               isreal=isrealB(φ), kwargs...) =
      SymmetricBasis(φ, symgrp,
                     PIBasis(basis1p, symgrp, Bsel;
                             isreal = isrealAA(φ), kwargs..., property = φ);
                     isreal=isreal)

SymmetricBasis(φ::AbstractProperty, symgrp::SymmetryGroup, pibasis::PIBasis;
               isreal=false) =
      SymmetricBasis(φ, symgrp, pibasis, isreal ? Base.real : Base.identity)

function SymmetricBasis(φ::TP, symgrp::SymmetryGroup, pibasis::PIBasis,
                        _real) where {TP <: AbstractProperty}

   # AA index -> AA spec
   AAspec = get_spec(pibasis)

   # construct the reverse mapping AA spec -> iAA and A spec -> iA
   invAAspec = Dict{Any, Int}()
   for (iAA, AA) in enumerate(AAspec)
      invAAspec[AA] = iAA
   end
   invAspec = Dict{Any, Int}()
   for (iA, A) in enumerate(get_spec(pibasis.basis1p))
      invAspec[A] = iA
   end

   # allocate the datastructure that computes and caches the
   # coupling coefficients
   # TODO: should this be stored with the basis?
   #       or maybe written to a file on disk? and then flushed every time
   #       we finish with a basis construction???
   # TODO: for sure this needs to become a function of the symmetry group?
   rotc = Rot3DCoeffs(φ, Float64)
   # allocate triplet format
   TCC = coco_type(TP)
   Irow, Jcol, vals = Int[], Int[], TCC[]
   # count the number of PI basis functions = number of rows
   idxB = 0

   # loop through AA basis, but skip most of them ...
   for (iAA, AA) in enumerate(AAspec)
      # determine whether we need to compute coupling coefficients for this
      # basis function or whether it will be included in a different
      # coco computation?
      if !is_refbasisfcn(symgrp, AA)
         continue
      end
      # compute the cocos
      U, AAcols = coupling_coeffs(symgrp, AA, rotc)

      # loop over the rows of U -> each specifies a basis function which we now
      # need to incorporate into the basis
      for irow = 1:size(U, 1)
         idxB += 1
         # loop over the columns of U / over brows
         for (icol, bcol) in enumerate(AAcols)
            # put bcol into the correct order
            bcol_ordered = _get_ordered(bcol, invAspec)
            # this is a subtle step: bcol and bcol_ordered are equivalent
            # permutation-invariant basis functions. This means we will
            # add the same PI basis function several times, but in the call to
            # `sparse` the values will just be added.
            if !haskey(invAAspec, bcol_ordered)
               @warn("bcol_ordered not in AA-spec")
               @show bcol_ordered
               # @show degree(bcol_ordered, NaiveTotalDegree(), pibasis.basis1p)
               error("bcol_ordered not in AA-spec")
            end
            idxAA = invAAspec[bcol_ordered]
            push!(Irow, idxB)
            push!(Jcol, idxAA)
            push!(vals, U[irow, icol])
         end
      end
   end
   # TODO: filter and throw out everything that hasn't been used!!
   # create CSC: [   triplet    ]  nrows   ncols
   A2Bmap = sparse(Irow, Jcol, vals, idxB, length(AAspec))
   basis = SymmetricBasis(pibasis, A2Bmap, symgrp, _real)
   # clean up a bit, i.e. remove AA basis functions that we don't need
   # to evaluate the symmetric basis
   clean_pibasis!(basis)
   return basis
end

function SymmetricBasis(pibasis, A2Bmap, symgrp, _real)
   PROP = _real(eltype(A2Bmap))
   B_pool = VectorPool{PROP}()
   return SymmetricBasis(pibasis, A2Bmap, symgrp, _real, B_pool)
end

"""
produce the ordered tuple defining the AA basis function uniquely.
"""
function _get_ordered(bb, invAspec)
   iAs = [invAspec[b] for b in bb]
   return bb[ sortperm(iAs, rev = true) ]
end

# ------------------- exporting the basis spec

# this doesn't provide the "full" specification, just collects
# the n and l but not the m or coupling coefficients.

function get_spec(basis::SymmetricBasis)
   spec_AA = get_spec(basis.pibasis)
   spec_B = []
   for iB = 1:length(basis)
      iAA = findfirst(norm.(basis.A2Bmap[iB, :]) .!= 0)
      push!(spec_B, get_sym_spec(basis.symgrp, spec_AA[iAA]))
   end
   return identity.(spec_B)
end

# ---------------- sparsification and cleanup codes

"""
`sparsify!(basis::SymmetricBasis; del = ..., keep = ...)`

sparsify the symmetric basis by either specifyin which basis functions to
keep or which ones to delete.
"""
function sparsify!(basis::SymmetricBasis;
                  del::Union{Nothing, AbstractVector{<: Integer}} = nothing,
                  keep::Union{Nothing, AbstractVector{<: Integer}} = nothing)
   if ( (del == nothing && keep == nothing ) ||
        (del != nothing && keep != nothing ) )
      error("sparsify!: must provide either del or keep kwarg but not both")
   end
   if del != nothing
      keep = setdiff(1:length(basis), del)
   end
   basis.A2Bmap = basis.A2Bmap[keep, :]
   clean_pibasis!(basis; atol = 0.0)
end

# this is also used for Filtering the AA basis if there are zero-rows in the A2B map

"""
Remove elements of the AA basis, when there are zero-rows in the A2B map, i.e.
when some AA basis elements are simply not required to evaluate the
symmetric basis.
"""
function clean_pibasis!(basis::SymmetricBasis; atol = 0.0)
   # get the zero-columns
   Inz = sort( findall(x -> x > atol, sum(norm, basis.A2Bmap, dims=1)[:]) )
   if length(Inz) < size(basis.A2Bmap, 2)
      # remove those columns from the A2Bmap
      sparsify!(basis.pibasis, Inz)
      # remove those columns from the A2Bmap
      basis.A2Bmap = basis.A2Bmap[:, Inz]
   end
   clean_1pbasis!(basis.pibasis)
   return basis
end

# ---------------- A modified sparse matmul

# TODO: move this stuff all to aux?

using SparseArrays: AbstractSparseMatrixCSC,
				        nonzeros, rowvals, nzrange

using LinearAlgebra: Transpose

function genmul!(C, A::AbstractSparseMatrixCSC, B, mulop)
    size(A, 2) == size(B, 1) || throw(DimensionMismatch())
    size(A, 1) == size(C, 1) || throw(DimensionMismatch())
    size(B, 2) == size(C, 2) || throw(DimensionMismatch())
    nzv = nonzeros(A)
    rv = rowvals(A)
    fill!(C, zero(eltype(C)))
    for k in 1:size(C, 2)
        @inbounds for col in 1:size(A, 2)
            αxj = B[col,k]
            for j in nzrange(A, col)
                C[rv[j], k] += mulop(nzv[j], αxj)
            end
        end
    end
    return C
end

function genmul!(C, xA::Transpose{<:Any,<:AbstractSparseMatrixCSC}, B, mulop)
   A = xA.parent
   size(A, 2) == size(C, 1) || throw(DimensionMismatch())
   size(A, 1) == size(B, 1) || throw(DimensionMismatch())
   size(B, 2) == size(C, 2) || throw(DimensionMismatch())
   nzv = nonzeros(A)
   rv = rowvals(A)
   fill!(C, zero(eltype(C)))
   for k in 1:size(C, 2)
       @inbounds for col in 1:size(A, 2)
           tmp = zero(eltype(C))
           for j in nzrange(A, col)
               tmp += mulop(nzv[j], B[rv[j],k])
           end
           C[col,k] += tmp
       end
   end
   return C
end

function genmul(A, B, mulop)
   T = typeof(mulop(A[1], B[1]))
   C = Array{T}(undef, (size(A, 1), size(B)[2:end]...) )
   return genmul!(C, A, B, mulop)
end

# ---------------- Evaluation code

function evaluate(basis::SymmetricBasis, cfg::UConfig)
   AA = evaluate(basis.pibasis, cfg)
   B = evaluate(basis, AA)
   release!(AA)
   return B
end

# NOTE: Nasty and completely not understood type instability here
function evaluate!(B, basis::SymmetricBasis, cfg::UConfig)
   AA = evaluate(basis.pibasis, cfg)
   evaluate!(B, basis, AA)
   release!(AA)
   return B
end

evaluate(basis::SymmetricBasis, AA::AbstractVector{<: Number}) =
      genmul(basis.A2Bmap, AA, (a, b) -> basis.real(a * b))

evaluate!(B, basis::SymmetricBasis, AA::AbstractVector{<: Number}) =
      genmul!(B, basis.A2Bmap, AA, (a, b) -> basis.real(a * b))

# ---------------- gradients

# args... could be nothing or sym
function evaluate_ed(basis::SymmetricBasis, cfg::UConfig, args...)
   AA, dAA = evaluate_ed(basis.pibasis, cfg, args...)
   B, dB = evaluate_ed(basis, AA, dAA)
   release!(AA)
   release!(dAA)
   return B, dB
end

function evaluate_ed(basis::SymmetricBasis, AA::AbstractVector{<: Number}, dAA)
   B = genmul(basis.A2Bmap, AA, (a, b) -> basis.real(a * b))
   dB = genmul(basis.A2Bmap, dAA, (a, b) -> basis.real( ACE.coco_o_daa(a, b) ))
   return B, dB
end

evaluate_d(basis::SymmetricBasis, args...) = evaluate_ed(basis, args...)[2]

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

function scaling(basis::SymmetricBasis, p)
   wwpi = scaling(basis.pibasis, p)
   wwrpi = abs2.(norm.(basis.A2Bmap)) * abs2.(wwpi)
   return sqrt.(wwrpi)
end


text:
using Test
using MechGluecode
using MechanicalUnits
using MechanicalUnits: Unitfu.promote_to_derived
@import_expand(m, s, N) # Takes long if defined late)
using DifferentialEquations

@testset "Unitful time with unitful function" begin
    u0 = 30.0N
    tspan = (0.0, 10.0)s
    prob1 = ODEProblem((du,u,t,p) -> (du[1] = -0.2∙s⁻¹ * u[1]), [u0], tspan)
    prob2 = ODEProblem((u,t,p)    ->         (-0.2∙s⁻¹ * u[1]), u0, tspan)
    prob3 = ODEProblem((u,t,p)    ->         [-0.2∙s⁻¹* u[1]], [u0],tspan)
    so1 =  solve(prob1, Euler(), dt = 1s)
    so2 =  solve(prob2, Midpoint())
    so3 =  solve(prob3, Tsit5())
    plot(so1, denseplot = false, marker = true)
    plot(so2, denseplot = false, marker = true)
    plot(so3, denseplot = false, marker = true)
    sf1(t) = so1(t)[1]
    sf2(t) = so2(t)[1]
    sf3(t) = so3(t)[1]
    plot([sf1, sf2, sf3], xlims =(0,10)s, ylims =(0,40)N)
    plot([t-> sf1(t) / sf3(t), t-> sf2(t) / sf3(t)], xlims =(0,10)s, ylims = (0.5,1.1), minorgrid=true, grid = (:y, :olivedrab, :dot, 1, 0.9))

end
@testset "Unitful time with unitless function" begin
    u0 = 30.0
    tspan = (0.0, 10.0)s
    prob1 = ODEProblem((du,u,t,p) -> (du[1] = -0.2∙s⁻¹ * u[1]), [u0], tspan)
    prob2 = ODEProblem((u,t,p)    ->         (-0.2∙s⁻¹ * u[1]), u0, tspan)
    prob3 = ODEProblem((u,t,p)    ->         [-0.2∙s⁻¹* u[1]], [u0],tspan)
    so1 =  solve(prob1, Euler(), dt = 1s)
    so2 =  solve(prob2, Midpoint())
    so3 =  solve(prob3, Tsit5())
    # don't work plot(so1, denseplot = false, marker = true)
    # don't work plot(so2, denseplot = false, marker = true)
    # don't work plot(so3, denseplot = false, marker = true)
    sf1(t) = so1(t)[1]
    sf2(t) = so2(t)[1]
    sf3(t) = so3(t)[1]
    plot([sf1, sf2, sf3], xlims =(0,10)s, ylims =(0,40))
    plot([t-> sf1(t) / sf3(t), t-> sf2(t) / sf3(t)], xlims =(0,10)s, ylims = (0.5,1.1), minorgrid=true, grid = (:y, :olivedrab, :dot, 1, 0.9))
end

@testset "Plot ODE quantity" begin
    # Don't give us kN from from solution in this case, just defaults.
    promote_to_derived()
    f1 = (y, p, t) -> -1.5y / 0.3s
    u0 = 1.5N
    tspan = (0.0s, 1.0s)
    prob = ODEProblem(f1, u0, tspan)
    sol = solve(prob)
    plot(sol)
end

@testset "Plot ODE system same quantities" begin
    # This is not a good example, as the units are incorrect and have no
    # known physical interpretation
    σ = 10/s
    ρ = 28m
    β = (8 / 3)m
    function lorenz!(du, u, p, t)
        du[1] = (u[2] - u[1])∙σ
        du[2] = (u[1] * (ρ - u[3])) / (m∙s) - u[2]/s
        du[3] =  u[1] * u[2] / (m∙s)    - β ∙ u[3]/(m∙s)
    end
    u0 = [1.0, 0.0, 0.0]m
    tspan = (0.0,100.0)s
    prob = ODEProblem(lorenz!,u0,tspan)
    sol = solve(prob)
    # Plot all variables vs time
    plot(sol)
    # Plot in phase space
    plot(sol, vars=(1,2,3))
end
#=
@testset "Plot ODE system single pendulum" begin
    # Second order non-linear, undamped pendulum:
    # mLθ'' + mg∙sin(θ)
    # where m is mass, unit kg
    #       L is length, unit m
    #       θ is angle to vertical, strictly no unit
    #       g is gravity's acceleration
    # Acceleration
    # θ'' + k * sin(θ) = 0
    #    where k = g / L , unit(k) = s⁻²
    # Reduced to first order non-linear system by defining:
    # θ = y[1]
    # θ' = y[2]
    # It follows that
    # θ'' = y'[2]
    # and
    # y'[1] = y[2]
    # y'[2] = -k * sin(y[1])
    # To solve numerically, we define the derivative.
    # For efficiency, it's evaluated in place, i.e. dy changes at every call.
    function f!(dy, y, k, t)
        dy[1] = y[2]
        dy[2] = -k ∙ sin(y[1])
    end
    # Fixed parameter
    k = upreferred(g / 1.0m)
    # Initial conditions
    θ0 = 0.1
    θ0_vel = 0.0s⁻¹
    θ0_acc = -k *sin(θ0)

    y0 = [θ0, θ0_vel]
    y0 = ArrayPartition([θ0_vel, θ0_acc])
    tspan = (0.0, 100.0)s
    prob = ODEProblem(f!, y0, tspan, k)
    sol = solve(prob)
    # Plot all variables vs time
    plot(sol)
    # Plot in phase space
    plot(solve(prob), vars=(1,2))
end
=#


text:
# from Polynomials.jl (Poly)

#=
poly = (x-1)^7
polycoeffs = [1.,-7.,21.,-35.,35.,-21.,7.,-1.]
xx = sqrt(2.5)
bigxx = sqrt(big(5)/2)
poly=Poly(polycoeffs)
poly1=polyval(xx,poly)
poly2=@Horner(xx,1.,-7.,21.,-35.,35.,-21.,7.,-1.)
poly3=@HornerFMA(xx,1.,-7.,21.,-35.,35.,-21.,7.,-1.)
poly4=parse(Float64, string(polyval(bigxx,poly))) # Float64(polyval(bigxx,poly))

poly4== -0.02238511157128063
poly3== -0.02238511157129503
poly2== -0.022385111571300076
poly1== -0.022385111571300076
=#

"""
* `FastPoly{T<:Number}(a::Vector)`: Construct a polynomial from its coefficients, lowest order first. That is if `p=a_n x^n + ... + a_2 x^2 + a_1 x^1 + a_0`, we construct this through `Poly([a_0, a_1, ..., a_n])`.
Example:
```
FastPoly([1,0,3,4])    # FastPoly(1 + 3x^2 + 4x^3)
```
An optional variable parameter can be added:
```
FastPoly([1,2,3], :s)       # FastPoly(1 + 2s + 3s^2)
```
"""
immutable FastPoly{T<:Real}
    a::Vector{T}
    var::Symbol

    function FastPoly{S}(a::Vector{S}, var)
        # if a == [] we replace it with a = [0]
        if length(a) == 0
            return new{S}(zeros(S,1), Symbol(var))
        else
            # determine the last nonzero element and truncate a accordingly
            a_last = max(1,findlast(a))
            new{S}(a[1:a_last], Symbol(var))
        end
    end
end

FastPoly{S}(a::Vector{S}, var) = FastPoly{S}(a,var)
FastPoly{S}(a::Vector{S}) = FastPoly{S}(a,:x)

#=
   via Julia High Performance by Avik Sengupta (pp 59-62)
=#
# from base
macro Horner(x, p...)
    ex = esc(p[end])
    for i = length(p)-1:-1:1
        ex = :(muladd(t, $ex, $(esc(p[i]))))
    end
    Expr(:block, :(t = $(esc(x))), ex)
end

polyvalHorner(x::Real) = @Horner(x, coeffs)

# from base, with fma instead of muladd
macro HornerFMA(x, p...)
    ex = esc(p[end])
    for i = length(p)-1:-1:1
        ex = :(fma(t, $ex, $(esc(p[i]))))
    end
    Expr(:block, :(t = $(esc(x))), ex)
end

polyvalHornerFMA(x::Union{Float64,Float32}) = @HornerFMA(x, coeffs)

#=
https://github.com/JuliaLang/julia/pull/2987/commits/9c24795ac2918df7b8cba9bb129db8c535d321e8
+# evaluate p[1] + x * (p[2] + x * (....)), i.e. a polynomial via Horner's rule
 +macro horner(x, p...)
 +    ex = p[end]
 +    for i = length(p)-1:-1:1
 +        ex = :($(p[i]) + $x * $ex)
 +    end
 +    ex
 +end
 +
 +# Compute the inverse of the error function: erf(erfinv(x)) == x,
 +# using the rational approximants tabulated in:
 +#     J. M. Blair, C. A. Edwards, and J. H. Johnson, "Rational Chebyshev
 +#     approximations for the inverse of the error function," Math. Comp. 30,
 +#     pp. 827--830 (1976).
 +#         http://dx.doi.org/10.1090/S0025-5718-1976-0421040-7
 +#         http://www.jstor.org/stable/2005402
 +function erfinv(x::Float64)
 +    a = abs(x)
 +    if a >= 1.0
 +        if x == 1.0
 +            return inf(Float64)
 +        elseif x == -1.0
 +            return -inf(Float64)
 +        end
 +        throw(DomainError())
 +    elseif a <= 0.75 # Table 17 in Blair et al.
 +        t = x*x - 0.5625
 +        return x * @horner(t, 0.16030_49558_44066_229311e2,
 +                             -0.90784_95926_29603_26650e2,
 +                              0.18644_91486_16209_87391e3,
 +                             -0.16900_14273_46423_82420e3,
 +                              0.65454_66284_79448_7048e2,
 +                             -0.86421_30115_87247_794e1,
 +                              0.17605_87821_39059_0) /
 +                   @horner(t, 0.14780_64707_15138_316110e2,
 +                             -0.91374_16702_42603_13936e2,
 +                              0.21015_79048_62053_17714e3,
 +                             -0.22210_25412_18551_32366e3,
 +                              0.10760_45391_60551_23830e3,
 +                             -0.20601_07303_28265_443e2,
 +                              0.1e1)
 =#
#=
Base:math.jl

# evaluate p[1] + x * (p[2] + x * (....)), i.e. a polynomial via Horner's rule
macro horner(x, p...)
    ex = esc(p[end])
    for i = length(p)-1:-1:1
        ex = :(muladd(t, $ex, $(esc(p[i]))))
    end
    Expr(:block, :(t = $(esc(x))), ex)
end

# Evaluate p[1] + z*p[2] + z^2*p[3] + ... + z^(n-1)*p[n].  This uses
# Horner's method if z is real, but for complex z it uses a more
# efficient algorithm described in Knuth, TAOCP vol. 2, section 4.6.4,
# equation (3).
macro evalpoly(z, p...)
    a = :($(esc(p[end])))
    b = :($(esc(p[end-1])))
    as = []
    for i = length(p)-2:-1:1
        ai = Symbol("a", i)
        push!(as, :($ai = $a))
        a = :(muladd(r, $ai, $b))
        b = :($(esc(p[i])) - s * $ai) # see issue #15985 on fused mul-subtract
    end
    ai = :a0
    push!(as, :($ai = $a))
    C = Expr(:block,
             :(x = real(tt)),
             :(y = imag(tt)),
             :(r = x + x),
             :(s = muladd(x, x, y*y)),
             as...,
             :(muladd($ai, tt, $b)))
    R = Expr(:macrocall, Symbol("@horner"), :tt, map(esc, p)...)
    :(let tt = $(esc(z))
          isa(tt, Complex) ? $C : $R
      end)
end

  http://stackoverflow.com/questions/28077057/julia-evalpoly-macro-with-varargs

  from Erik Schnetter
  c=[1,2,3,4]
  @eval begin
    function f()
        @evalpoly 0.5 $(c...)
    end
end
=#
#=
Investigation of different methods of fast polynomial evaluation
Gavin S. Reynolds
September 2010
MSc in High Performance Computing
The University of Edinburgh
Year of Presentation: 2010

The Optimised Brute Force function pre-calculates the value of x^2 and use it where possible in
building up the powers for each monomial. For example, if a monomial contains x^7, instead of
7 multiplications of x by x, the function would multiply x^2 by x^2three times then multiply it by x
once. In this example, this would save three multiplications.

function hornerform(x,coeffs)
    n = length(coeffs)
    z = coeffs[1]*x
    for i in 2:(n-1)
         z = (z + coeffs[i])*x
    end
    z + coeffs[n]
end

function bruteforce(x,coeffs)
    n = length(coeffs)
    monomials = zeros(typeof(x),n)
    for i in 1:(n-1)
        monomials[i] = coeffs[i] * (x^(n-i))
    end
    monomials[n] = coeffs[n]
    sum(monomials)
end

function optforce(x,coeffs)
    n = length(coeffs)
    monomials = zeros(typeof(x),n)
    x2 = x*x
    for i in 1:(n-1)
        monomials[i] = coeffs[i]
        numsteps = (n-i) >> 1 # fld(n-i,2)
        for k in 1:numsteps
           monomials[k] *= x2
        end
        if isodd(n-i)
            monomials[i] *= x
        end
    end
    monomials[n] = coeffs[n]
    sum(monomials)
end

end
=#


text:
include("header.jl")
using SpecialFunctions

@testset "unary" begin

    function frand(f,t,d...)
        r = rand(t,d...) .* t(0.5) .+ t(0.25)
        if in(f,(acosh,asec))
            return 1 ./ r
        else
            return r
        end
    end

    bcast(f)=(x->broadcast(f,x))

    unary_fns = Any[]
    for f in Knet.unary_ops
        if isa(f,Tuple); f=f[2]; end
        push!(unary_fns, eval(Meta.parse(f)))
    end

    for f in unary_fns
        #@show f
        bf = bcast(f)
        for t in (Float32, Float64)
            #@show f,t
            sx = frand(f,t)
            @test isa(f(sx),t)
            @test gradcheck(f, sx)
            for n in (1,(1,1),2,(2,1),(1,2),(2,2))
                #@show f,t,n
                ax = frand(f,t,n)
                @test gradcheck(bf, ax)
                if gpu() >= 0
                    gx = KnetArray(ax)
                    cy = bf(ax)
                    gy = bf(gx)
                    @test isapprox(cy,Array(gy))
                    @test gradcheck(bf, gx)
                end
            end
        end
    end
end

nothing


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

@testset "parameter_mapping" begin
end


text:
"""
    TwoLevelChain(N_atoms::Int; Γ = 0.0)

Creates a chain of `N_atoms` two-level atoms with spontaneous decay rate `Γ`.
The coupling rate is set to 1, but can be modified by scaling the Hamiltonian and coupling operator.
Uses `NLevelBasis(2)` as the single atom basis.
"""
function TwoLevelChain(N_atoms::Int; Γ = 0.0)
	singleAtomBasis = NLevelBasis(2)
	σ_WG_i, σ_GW_i  = transition(singleAtomBasis, 2, 1), transition(singleAtomBasis, 1, 2)

	if N_atoms == 1
		groundstateAtoms = nlevelstate(singleAtomBasis, 1)
		H_Atom      = 0.0 * one(singleAtomBasis)
		Dissipators = iszero(Γ) ? [] : [√(Γ) * σ_GW_i]

		return H_Atom, Dissipators, σ_GW_i, groundstateAtoms
	end

	multiAtomBasis = reduce(⊗, repeat([singleAtomBasis], N_atoms))

	groundstateAtoms = reduce(⊗, repeat([nlevelstate(singleAtomBasis, 1)], N_atoms))

	σ_col = sum(i -> embed(multiAtomBasis, i, σ_GW_i), [1:N_atoms;])

	Dissipators = iszero(Γ) ? [] : [embed(multiAtomBasis, i, √(Γ) * σ_GW_i) for i ∈ 1:N_atoms]

	H_Atoms = -0.5im * sum(i -> sum(
			j -> embed(multiAtomBasis, i, σ_WG_i) * embed(multiAtomBasis, j, σ_GW_i),
			[1:i-1;]
		), [2:N_atoms;]
	)
	H_Atoms = (H_Atoms + dagger(H_Atoms))

	return H_Atoms, Dissipators, σ_col, groundstateAtoms
end


text:
#=
The Computer Language Benchmarks Game
 https://salsa.debian.org/benchmarksgame-team/benchmarksgame/

 direct transliteration of the swift#3 program by Ralph Ganszky and Daniel Muellenborn

 modified for Julia 1.0 by Simon Danisch
 tweaked for performance by maltezfaria and Adam Beckmeyer
=#

push!(LOAD_PATH, "MD5");

using Base.Cartesian
using Base64.Base64
using MD5
# include("MD5/MD5.jl")

# Calculate the byte to print for a given vector of 8 real numbers cr
# and a given imaginary component ci. This function should give the
# same result whether prune is true or false but may be faster or
# slower depending on the input.
function mand8(cr, ci, prune)
    zr = zi = tr = ti = @ntuple 8 _-> 0.0

    # In cases where the last call to mand8 resulted in 0x00, the next
    # call is much more likely to result in 0x00, so it's worth it to
    # check several times if the calculation can be aborted
    # early. Otherwise, the relatively costly check can be eliminated.
    if prune
        for _=1:10
            for _=1:5
                zr, zi, tr, ti = calc_sum(cr, ci, zr, zi , tr, ti)
            end
            # If all values are already > 4.0, we already know no
            # values will be <= 4.0 at end of iterations. Reduction
            # with binary & is used here instead of all function to
            # avoid branches caused by short-circuiting behavior.
            foldl(&, @ntuple(8, k-> tr[k] + ti[k] > 4.0)) && return 0x00
        end
    else
        for _=1:50
            zr, zi, tr, ti = calc_sum(cr, ci, zr, zi , tr, ti)
        end
    end

    # For each value <= 4.0, the k-th bit (counting from left) is set to 1
    foldl(|, @ntuple(8, k-> ifelse(tr[k] + ti[k] <= 4.0, 0x80 >> (k - 1), 0x00)))
end

# Single iteration of mandelbrot calculation for vector r of real
# components and vector i or imaginary components.
@inline function calc_sum(cr, ci, zr, zi, tr, ti)
    # Using broadcasting (r2 = r .* r) generates operations on llvm
    # <8 x double> vectors even with --cpu-target=core2 (widest simd
    # register on core2 is <2 x double>). @ntuple results in better
    # codegen (uses <2 x double>).
    zi = @ntuple 8 k-> muladd(zr[k] + zr[k], zi[k], ci)
    zr = @ntuple 8 k-> tr[k] - ti[k] + cr[k]
    tr = @ntuple 8 k-> zr[k] * zr[k]
    ti = @ntuple 8 k-> zi[k] * zi[k]
    zr, zi, tr, ti
end

# Write n by n portable bitmap image of mandelbrot set to io
function mandelbrot(io, n_raw)
    n = (n_raw + 7) ÷ 8 * 8

    # Precalculate real coordinates to check
    xvals = Float64[2i/n - 1.5 for i=0:n-1]
    # Precalculate imaginary coordinates to check
    yvals = Float64[2i/n - 1.0 for i=0:n-1]

    # Create a vector of bytes to output
    out = Vector{UInt8}(undef, n * n ÷ 8)
    # For each row (each imaginary coordinate), spawn a thread to fill
    # out values. At small values of n, this is too fine-grained of
    # parallelism to really be efficient, but it works well for large n.
    @sync for y=1:n
        # Threads.@spawn allows dynamic scheduling instead of static scheduling
        # of Threads.@threads macro.
        # See <https://github.com/JuliaLang/julia/issues/21017>.
        Threads.@spawn @inbounds begin
            ci = yvals[y]
            startofrow = (y - 1) * n ÷ 8
            # The first iteration within a row will generally return 0x00
            prune = true
            for x=1:8:n
                # Calculate whether the (x:x+7)-th real coordinates with
                # the y-th imaginary coordinate belong to the
                # mandelbrot set.
                byte = mand8(@ntuple(8, k-> xvals[x+k-1]), ci, prune)
                out[startofrow + x÷8 + 1] = byte
                prune = byte == 0x00
            end
        end
    end

    println("P4\n$n $n")
    println(bytes2hex(MD5.md5("test")))
end

isinteractive() || mandelbrot(stdout, parse(Int, ARGS[1]))


text:
#=
kerneldecomposition:
- Julia version:
- Author: qd4314
- Date: 2019-03-19
=#

function metacheck(N)
    for n=1:N
        checkexpansionHG(n)
    end
end

function checkexpansionHG(n,g = 0.3)
    ntheta = 1000
    theta = range(0,stop=2*pi,length=ntheta)
	hg(g, mu) = 1/(4*pi)*(1-g^2)./(1+g^2-2*g*mu).^(3/2)
    y = hg.(g,cos.(theta))
	hgcoeffs = [(2*l+1)/(4*pi)*g^l for l = 0:n]
    v = sf_legendre_Pl_array.(n,cos.(theta))
    yapprox = [sum(dot(v[i],hgcoeffs)) for i=1:ntheta]
    println(norm(yapprox .- y) / norm(y))
end

function metadecomposition(N)
    err = zeros(N)
    for n=1:N
		_,_,_,err[n] = kerneldecomposition(n);
    end
    return(err)
end

function kerneldecomposition(KS::KineticSolver,order=20)
    g = 0.01
    HG = computeScatteringKernelHenyeyGreenstein(KS,g)

    index(l,m) = l*(order+1)+m+1

    Ωs = KS.Q.pointsxyz
    nq = KS.nquadpoints
    weights = KS.Q.weights

    Norder = (order+1)*(order+1)
    S = zeros(Norder,Norder)
    counter = 1
    for l=0:order
        for m=-l:+l
            S[counter,counter] =  g^l
            counter += 1
        end
    end

    Y = zeros(order +1,2*(order +1)+1,nq)
    YY = zeros(Norder,nq)
    @polyvar x y z
    counter = 1
    for l=0:order
        for m=-l:l
            sphericalh = ylm(l,m,x,y,z)
            for q = 1 : nq
                Y[l+1,m+l+1,q] = sphericalh(x=>Ωs[q,1],y=>Ωs[q,2],z=>Ωs[q,3])
                YY[counter,q] =  sphericalh(x=>Ωs[q,1],y=>Ωs[q,2],z=>Ωs[q,3])
            end
            counter += 1
        end
    end
    normY = zeros(Norder);
    for i = 1:Norder
        for k = 1:nq
			normY[i] = normY[i] + YY[i,k]^2
        end
    end
    normY = sqrt.(normY)

    Mat = zeros(nq,nq)
    O = zeros(nq,Norder)
    M = zeros(Norder,nq)
    for i=1:nq
        for j=1:nq
            val = 0
            counter = 1
            for l=0:order
                for m=-l:l
                    val +=g^l*Y[l+1,m+l+1,i] * Y[l+1,m+l+1,j]*weights[j]
					O[i,counter] = YY[counter,i] # / normY[counter]
					M[counter,j] = YY[counter,j]*weights[j]# / normY[counter]
                    counter += 1
                end
            end
            Mat[i,j] = val
        end
    end

	v = O*S*M*ones(nq)

	counter = 1
	for q=1:nq
		O[q,:] /= v[q]
	end

	#println(maximum(abs.(v.-HG*ones(nq))))

	mystring = string(norm((O*S*(M) .-HG))/norm(HG))
	println("Error for HG - Prod(Matrixdecomposition)= " * mystring)
	err = norm(O*S*M .- HG) / norm(HG)
    return O,S,M,err

end

function kerneldecomposition(order=20)
    KS =KineticSolver("config.txt")
    kerneldecomposition(KS,order)
end

max_stars_repo_name: camminady/SN

max_issues_repo_name: camminady/SN

max_forks_repo_name: camminady/SN


text:
import MathOptInterface
const MOI = MathOptInterface
const MOIB = MOI.Bridges
const MOIU = MOI.Utilities
const MOIF = MOI.FileFormats

function load_problem!(m::MOI.AbstractOptimizer, fname::String)
    MOI.empty!(m)

    # Read problem from MPS file
    mps = MOIF.MPS.Model()
    MOI.read_from_file(mps, fname)

    # Remove all integer constraints
    l = MOI.get(mps, MOI.ListOfConstraints())
    for cidx in MOI.get(mps, MOI.ListOfConstraintIndices{MOI.SingleVariable, MOI.Integer}())
        MOI.delete(mps, cidx)
    end
    for cidx in MOI.get(mps, MOI.ListOfConstraintIndices{MOI.SingleVariable, MOI.ZeroOne}())
        MOI.delete(mps, cidx)
        # Add lower and upper bounds
        v = MOI.VariableIndex(cidx.value)
        MOI.add_constraint(mps, MOI.SingleVariable(v), MOI.GreaterThan(0.0))
        MOI.add_constraint(mps, MOI.SingleVariable(v), MOI.LessThan(1.0))
    end

    # Set the objective to be
    Tobj = MOI.get(mps, MOI.ObjectiveFunctionType())
    if Tobj != MOI.ScalarAffineFunction{Float64}
        fobj = MOI.get(mps, MOI.ObjectiveFunction{Tobj}())
        # convert to ScalarAffine
        fobj_ = MOI.ScalarAffineFunction{Float64}(fobj)
        MOI.set(mps, MOI.ObjectiveFunction{typeof(fobj_)}(), fobj_)
    end

    # Copy model
    MOI.copy_to(m, mps, copy_names=false)

    return m
end

function run(model::MOI.AbstractOptimizer)
    println("Solver name: ", MOI.get(model, MOI.SolverName()))
    try
        tsolve = @elapsed MOI.optimize!(model)

        # TODO: export solution
        # We can't do that since MOI's MPS reader may not be consistent between runs
        println("Solve time: ", tsolve)
        println("Termination status: ", MOI.get(model, MOI.TerminationStatus()))
        rcount = MOI.get(model, MOI.ResultCount())
        println("Result count: ", rcount)
        if rcount > 0
            println("Primal status: ", MOI.get(model, MOI.PrimalStatus()))
            println("Dual status: ", MOI.get(model, MOI.DualStatus()))
        end
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
    end
end

max_stars_repo_name: mtanneau/LPBenchmarks

max_issues_repo_name: mtanneau/LPBenchmarks

max_forks_repo_name: mtanneau/LPBenchmarks


text:
using PyCall, FileIO, ImageCore, ImageTransformations, CUDA, Tar, SHA
import Knet
using Knet.Layers21, Knet.Ops21
using Knet.Train20: param
using Knet.Ops20: softmax
typename(p::PyObject) = pytypeof(p).__name__
tf = pyimport("tensorflow")
include("../Models21.jl")

# Model translator
function keras2knet(p::PyObject; xtest = nothing, nlayers = -1)
    layers = nlayers >= 1 ? p.layers[1:nlayers] : copy(p.layers)
    model = Block(; name=p.name)
    while !isempty(layers)
        layer1 = popfirst!(layers)
        ltype = typename(layer1)
        # This could also pop one or more layers:
        layer2 = eval(Meta.parse("keras_$ltype"))(layer1, layers)
        layer2 !== nothing && push!(model, layer2)
    end
    if xtest !== nothing # x = rand(Float32,1,224,224,3)
        y1 = p.predict(xtest)
        y2 = model(xtest)
        @assert vec(Array(y1)) ≈ vec(Array(y2))
    end
    return model
end

keras_InputLayer(l, layers) = nothing

keras_ZeroPadding2D(l, layers) = ZeroPad(l.padding)

function keras_Conv2D(l, layers)
    # channels_last corresponds to inputs with shape nhwc (default)
    # channels_first corresponds to inputs with shape nchw
    # filters in channels_last are in hwxy order
    # knet default is whcn for images and whxy for filters
    w = param(permutedims(l.weights[1].numpy(),(2,1,3,4)))
    bias = (l.bias === nothing ? nothing : param(reshape(l.bias.numpy(), (1,1,:,1))))
    normalization, activation = keras_normact(layers)
    channelmajor = false
    crosscorrelation = true
    dilation = l.dilation_rate
    groups = l.groups
    padding = (l.padding == "valid" ? 0 : l.padding == "same" ? (size(w)[1:2] .- 1) .÷ 2 : l.padding)
    stride = l.strides
    Conv(w; activation, normalization, bias, channelmajor, crosscorrelation, dilation, groups, padding, stride)
end

function keras_DepthwiseConv2D(l, layers)
    w = l.weights[1].numpy()
    @assert size(w,4) == 1
    w = param(permutedims(w, (2,1,4,3))) ### (3,3,32,1) => (3,3,1,32)
    bias = (l.bias === nothing ? nothing : param(reshape(l.bias.numpy(), (1,1,:,1))))
    normalization, activation = keras_normact(layers)
    channelmajor = false
    crosscorrelation = true
    dilation = l.dilation_rate
    groups = size(w,4) ###
    padding = (l.padding == "valid" ? 0 : l.padding == "same" ? (size(w)[1:2] .- 1) .÷ 2 : l.padding)
    stride = l.strides
    Conv(w; activation, normalization, bias, channelmajor, crosscorrelation, dilation, groups, padding, stride)
end

function keras_normact(layers)
    if !isempty(layers) && typename(layers[1]) == "BatchNormalization"
        normalization = keras_BatchNormalization(popfirst!(layers), layers)
    else
        normalization = nothing
    end
    if !isempty(layers) && typename(layers[1]) == "ReLU"
        activation = keras_ReLU(popfirst!(layers), layers)
    else
        activation = nothing
    end
    normalization, activation
end

function keras_ReLU(l, layers)
    max_value = (l.max_value === nothing || l.max_value[] === nothing ? Inf : l.max_value[])
    threshold = l.threshold[]
    negative_slope = l.negative_slope[]
    if threshold == negative_slope == 0
        Op(relu; max_value)
    else
        Op(relu; max_value, threshold, negative_slope)
    end
end

function keras_BatchNormalization(l, layers)
    bparam(x) = param(reshape(x.numpy(), (1,1,:,1)))
    mean = bparam(l.moving_mean).value
    var = bparam(l.moving_variance).value
    bias = bparam(l.beta)
    scale = bparam(l.gamma)
    epsilon = l.epsilon
    update = 1-l.momentum
    BatchNorm(; mean, var, bias, scale, epsilon, update)
end

function keras_GlobalAveragePooling2D(l, layers)
    x->begin #GlobalAveragePooling2D
        y = pool(x; mode=mean, window=typemax(Int))
        reshape2d(y)
    end
end

function keras_Reshape(l, layers)
    Op(reshape, (l.target_shape..., :))
end

function keras_Dropout(l, layers)
    Op(dropout, l.rate)
end

function keras_Activation(l, layers)
    if l.activation.__name__ == "softmax"
        Op(softmax; dims=1)
    elseif l.activation.__name__ == "swish"
        swish
    else
        error("Activation $(l.activation.__name__) not yet implemented.")
    end
end

function keras_Dense(l, layers)
    w = param(permutedims(l.weights[1].numpy()))
    bias = param(l.bias.numpy())
    activation = keras_Activation(l, layers)
    Linear(w; bias, activation)
end

function keras_Add(l, layers)
    @warn "Add not defined" maxlog=1
    Op("Add")
end

function keras_Rescaling(l, layers)
    @warn "Rescaling not defined" maxlog=1
    Op("Rescaling")
end

function keras_TFOpLambda(l, layers)
    @warn "TFOpLambda not defined" maxlog=1
    Op("TFOpLambda")
end

function keras_Multiply(l, layers)
    @warn "Multiply not defined" maxlog=1
    Op("Multiply")
end

function keras_Flatten(l, layers)
    @warn "Flatten not defined" maxlog=1
    Op("Flatten")
end

function keras_Normalization(l, layers)
    @warn "Normalization not defined" maxlog=1
    Op("Normalization")
end

### Test and import utils

function kerastest(img="ILSVRC2012_val_00000001.JPEG")
    global pf, jf, px, jx
    pf = tf.keras.applications.MobileNet()
    jf = keras2knet(pf)
    px = imagenet_preprocess(img; normalize=mobilenet_normalize_tf, format="nhwc", atype=Array{Float32}) # 1,224,224,3
    jx = Knet.atype(permutedims(px,(3,2,4,1))) # 224,224,3,1
    @time py = pf(px).numpy() |> vec
    @time jy = jf(jx) |> Array |> vec
    @show jy ≈ py
    idx = sortperm(jy, rev=true)
    [ idx jy[idx] imagenet_labels()[idx] ]
end

function tf_layer_outputs(model, img="ILSVRC2012_val_00000001.JPEG"; normalize=identity, resolution=224)
    model.trainable = false
    x = imagenet_preprocess(img; normalize, format="nhwc", atype=Array{Float32}, resolution) # 1,224,224,3
    outputnames = string.(l.output for l in model.layers)
    outputs = []
    for l in Knet.progress(model.layers)
        if l.input isa PyObject
            x = l(x)
        else
            xs = map(string.(l.input)) do s
                outputs[findfirst(isequal(s), outputnames)]
            end
            x = l(xs)
        end
        push!(outputs, x)
    end
    (x->x.numpy()).(outputs)
end

function tf_layer_outputs_slow(model, img="ILSVRC2012_val_00000001.JPEG"; normalize=identity, resolution=224)
    model.trainable = false
    input = imagenet_preprocess(img; normalize, format="nhwc", atype=Array{Float32}, resolution) # 1,224,224,3
    # layer_output = base_model.get_layer(layer_name).output
    models = []
    for l in Knet.progress(model.layers)
        push!(models, tf.keras.models.Model(model.input, outputs=[l.output]))
    end
    outputs = []
    for m in Knet.progress(models)
        try
            push!(outputs, m.predict(input))
        catch e
            push!(outputs, nothing)
        end
    end
    return outputs
end

function tf_layer_output(model, layer=length(model.layers), img="ILSVRC2012_val_00000001.JPEG"; normalize=identity, resolution=224)
    model.trainable = false
    input = imagenet_preprocess(img; normalize, format="nhwc", atype=Array{Float32}, resolution) # 1,224,224,3
    m = tf.keras.models.Model(model.input, outputs=[model.layers[layer].output])
    m.predict(input)
end

function test_mobilenet_v2(file="mobilenet_v2_100_224_tf.jld2", img="ILSVRC2012_val_00000001.JPEG")
    pf = tf.keras.applications.MobileNetV2() # tf.keras.applications.mobilenet_v2.MobileNetV2()
    px = imagenet_preprocess(img; normalize=mobilenet_normalize_tf, format="nhwc", atype=Array{Float32}, resolution=224) # 1,224,224,3
    jf = MobileNet("mobilenet_v2_100_224_tf")
    jx = Knet.atype(permutedims(px,(3,2,4,1))) # 224,224,3,1
    jf(jx)
    setweights!(jf, file)
    isapprox(vec(pf(px).numpy()), softmax(vec(jf(jx))))
end

function import_mobilenet_tf(
    ; width=1, resolution=224, version=1,
    img="ILSVRC2012_val_00000001.JPEG",
)
    # global pf, jf0, jf, jx, jy, px, py, w
    @assert width in (1, 0.75, 0.5, 0.25)
    @assert resolution in (224, 192, 160, 128)
    alpha = width
    input_shape = (resolution, resolution, 3)

    if version == 1
        pf == tf.keras.applications.MobileNet(; alpha, input_shape)
        block=DWConv
        layout=mobilenet_v1_layout
        output=1024
    elseif version == 2
        pf == tf.keras.applications.MobileNetV2() # TODO: (; alpha, input_shape)
        block=MBConv
        layout=mobilenet_v2_layout
        output=1280
    end

    save_atype = Knet.array_type[]; Knet.array_type[] = Array{Float32}
    px = imagenet_preprocess(img; normalize=mobilenet_normalize_tf, format="nhwc", atype=Array{Float32}, resolution) # 1,224,224,3
    jx = Knet.atype(permutedims(px,(3,2,4,1))) # 224,224,3,1

    jf = MobileNet(; width, resolution, block, layout, output)
    jf(jx) # init weights
    jf0 = keras2knet(pf)
    w = getweights(jf0)
    if version == 1
        # MobileNetV1 ends with conv instead of linear
        w[end] = reshape(w[end],1000)
        w[end-1] = permutedims(reshape(w[end-1],(:,1000)))
    end
    setweights!(jf, w)

    @time py = pf(px).numpy() |> vec
    @time jy = jf(jx) |> vec |> softmax
    @assert @show jy ≈ py       # TODO: this does not match python script, preprocessing?
    Knet.array_type[] = save_atype

    alpha100 = round(Int, 100*width)
    model = "mobilenet_v$(version)_$(alpha100)_$(resolution)"
    saveweights("$(model).jld2", jf)
    run(`tar cf $(model).tar $(model).jld2`)
    sha1 = Tar.tree_hash("$(model).tar")
    run(`gzip $(model).tar`)
    sha2 = open("$(model).tar.gz") do f; bytes2hex(sha256(f)); end
    idx = sortperm(py, rev=true)
    pred = [ idx py[idx] imagenet_labels()[idx] ]
    display(pred[1:5,:]); println()
    @info "git-tree-sha1 = \"$sha1\""
    @info "sha256 = \"$sha2\""
end

function import_mobilenet_v2(
    ; width=1, resolution=224,
    tfmodel = tf.keras.applications.MobileNetV2,
    img="ILSVRC2012_val_00000001.JPEG",
)
    # global pf, jf0, jf, jx, jy, px, py, w
    save_atype = Knet.array_type[]; Knet.array_type[] = Array{Float32}
    @assert width in (1, 0.75, 0.5, 0.25)
    @assert resolution in (224, 192, 160, 128)
    px = imagenet_preprocess(img; normalize=mobilenet_normalize_tf, format="nhwc", atype=Array{Float32}, resolution) # 1,224,224,3
    jx = Knet.atype(permutedims(px,(3,2,4,1))) # 224,224,3,1
    pf = tfmodel()  # ; alpha=width, input_shape=(resolution,resolution,3))
    jf0 = keras2knet(pf)
    w = getweights(jf0)
    jf = MobileNet(; width, resolution, block=MBConv, layout=mobilenet_v2_layout, output=1280)
    jf(jx) # init weights
    setweights!(jf, w)
    @time py = pf(px).numpy() |> vec
    @time jy = jf(jx) |> vec |> softmax
    @assert @show jy ≈ py
    alpha = round(Int, 100*width)
    model = "mobilenet_v2_$(alpha)_$(resolution)"
    saveweights("$(model).jld2", jf)
    run(`tar cf $(model).tar $(model).jld2`)
    sha1 = Tar.tree_hash("$(model).tar")
    run(`gzip $(model).tar`)
    sha2 = open("$(model).tar.gz") do f; bytes2hex(sha256(f)); end
    idx = sortperm(py, rev=true)
    @info "git-tree-sha1 = \"$sha1\""
    @info "sha256 = \"$sha2\""
    Knet.array_type[] = save_atype
    [ idx py[idx] imagenet_labels()[idx] ]
end

function import_mobilenet_v3(   # WIP
    ;
    tfmodel = tf.keras.applications.MobileNetV3Small,
    layout = mobilenet_v3_small_layout,
    img="ILSVRC2012_val_00000001.JPEG",
    width=1,
    resolution=224)
    global pf, jf0, jf, jx, jy, px, py, w0
    save_atype = Knet.array_type[]; Knet.array_type[] = Array{Float32}
    @assert width in (1, 0.75, 0.5, 0.25)
    @assert resolution in (224, 192, 160, 128)
    # px = imagenet_preprocess(img; normalization="tf", format="nhwc", atype=Array{Float32}, resolution) # 1,224,224,3
    px = imagenet_preprocess(img; format="nhwc", atype=Array{Float32}, resolution) # 1,224,224,3
    px = 255 .* px #   x = layers.Rescaling(scale=1. / 127.5, offset=-1.)(x)
    # jx = Knet.atype(permutedims(px,(3,2,4,1))) # 224,224,3,1
    jx = imagenet_preprocess(img; format="whcn", atype=Array{Float32}, resolution) # 224,224,3,1
    jx = 2 .* jx .- 1
    pf = tfmodel()  # ; alpha=width, input_shape=(resolution,resolution,3))
    jf0 = keras2knet(pf)
    w0 = getweights(jf0)
    # MobileNet_tf ends with conv instead of linear
    w0[end-3] = permutedims(reshape2d(w0[end-3]))
    w0[end-2] = vec(w0[end-2])
    w0[end-1] = permutedims(reshape2d(w0[end-1]))
    w0[end] = vec(w0[end])
    jf = MobileNet(; width=1, resolution=224, input=16, output=(576,1024), classes=1000, activation=hardswish, tfpadding=true, bnupdate=0.001, bnepsilon=0.001, dropout=0.2, block=mobilenet_v3_block, layout)
    jf(jx) # init weights
    setweights!(jf, w0)
    @time py = pf(px).numpy() |> vec
    @time jy = jf(jx) |> vec |> softmax
    @show jy ≈ py
    # alpha = round(Int, 100*width)
    # model = "mobilenet_v2_$(alpha)_$(resolution)"
    # saveweights("$(model).jld2", jf)
    # run(`tar cf $(model).tar $(model).jld2`)
    # sha1 = Tar.tree_hash("$(model).tar")
    # run(`gzip $(model).tar`)
    # sha2 = open("$(model).tar.gz") do f; bytes2hex(sha256(f)); end
    # @info "git-tree-sha1 = \"$sha1\""
    # @info "sha256 = \"$sha2\""
    pidx = sortperm(py, rev=true)
    jidx = sortperm(jy, rev=true)
    Knet.array_type[] = save_atype
    [ pidx py[pidx] imagenet_labels()[pidx] jidx jy[jidx] imagenet_labels()[jidx] ]
end

function import_efficientnet(n=0)
    # to test with outs = tf_layer_outputs(pf; normalize=(x->(255 .* x)), resolution=240)
    tfmodel = getproperty(tf.keras.applications, "EfficientNetB$n")
    img="ILSVRC2012_val_00000001.JPEG"
    global pf, jf0, jf, jx, jy, px, py, w0
    model = lowercase(tfmodel.__name__)
    opts = efficientnet_models[model]
    save_atype = Knet.array_type[]; Knet.array_type[] = Array{Float32}
    px = imagenet_preprocess(img; normalize=(x->(255*x)), format="nhwc", atype=Array{Float32}, opts.resolution) # 1,224,224,3
    jx = imagenet_preprocess(img; normalize=efficientnet_normalize, format="whcn", atype=Array{Float32}, opts.resolution) # 224,224,3,1
    pf = tfmodel()  # ; alpha=width, input_shape=(resolution,resolution,3))
    jf0 = keras2knet(pf)
    w0 = getweights(jf0)
    jf = EfficientNet(; opts...)
    jf(jx) # init weights
    setweights!(jf, w0)
    @info "Calling tf model"
    @time py = pf(px).numpy() |> vec
    @info "Calling julia model"
    @time jy = jf(jx) |> vec |> softmax
    @show jy ≈ py
    saveweights("$(model).jld2", jf)
    run(`tar cf $(model).tar $(model).jld2`)
    sha1 = Tar.tree_hash("$(model).tar")
    run(`gzip $(model).tar`)
    sha2 = open("$(model).tar.gz") do f; bytes2hex(sha256(f)); end
    @info "git-tree-sha1 = \"$sha1\""
    @info "sha256 = \"$sha2\""
    pidx = sortperm(py, rev=true)
    jidx = sortperm(jy, rev=true)
    Knet.array_type[] = save_atype
    [ pidx py[pidx] imagenet_labels()[pidx] jidx jy[jidx] imagenet_labels()[jidx] ]
end


text:
################################################################
#
# representation of words as vector of pointers to generators
#
################################################################

@inline Base.push!(pw::Word, i) = push!(pw.gens_ptrs, i)
@inline Base.getindex(pw::Word, n) = pw.gens_ptrs[n]
@inline Base.setindex!(pw::Word, val, n) = pw.gens_ptrs[n] = val

@inline Base.iterate(pw::Word) = iterate(pw.gens_ptrs)
@inline Base.iterate(pw::Word, s) = iterate(pw.gens_ptrs, s)
@inline Base.length(pw::Word) = length(pw.gens_ptrs)
@inline Base.eltype(::Word{I}) where I = I
@inline Base.empty!(w::Word) = empty!(w.gens_ptrs)

@inline function Base.:(==)(pw1::Word, pw2::Word)
    return pw1.gens_ptrs == pw2.gens_ptrs
end

Base.hash(pw::Word, h::UInt) = hash(Word, hash(pw.gens_ptrs, h))

@doc doc"""
    Word(orb, gens_inv, pt[, op=^])
Return a `Word` `w` consising of pointers to inverses of generators `gens_inv`.

The group element `r` such that `first(orb)^r = pt` can be reconstructed via
`prod(gens_inv[i] for i in reverse(w))` e.g. if
 * `g = gens[1]*gens[3]*gens[2]` and
 * `δ = first(orb)^g`
then the result _may be_ `Word([2,3,1])`.
"""
function Word(gens_inv::AbstractVector{<:GroupElement}, orb::AbstractOrbit{T, I}, pt, op=^) where {T, I<:Integer}
    return append!(Word(I[]), gens_inv, orb, pt, op)
end

function Base.append!(w::Word, gens_inv, orb, pt, op=^)
    γ = pt
    while γ ≠ first(orb)
        idx = orb[γ]
        push!(w, idx)
        γ = op(γ, gens_inv[idx]) #assume that this is already the inverse!
    end
    return w
end

@doc doc"""
    (w::Word)(gens[, init])
Computes the evaluation of group word `w` as a group element in generators `gens`.

Optional `init` element is by default the group identity.
"""
@inline function (pw::Word)(gens::AbstractVector{<:GroupElement}, init=one(first(gens)))
    res = init
    @inbounds for i in pw
        res = GroupsCore.mul!(res, res, gens[i])
    end
    return res
end

@doc doc"""
    representative(orb, gens, pt[, ^])
Return a representative of the coset of `Stab(orb)` which takes `first(orb)`
to `pt` i.e. an element `g` of `⟨gens⟩` such that `first(orb)^g = pt`.
"""
function representative(gens::AbstractVector{<:GroupElement}, orb::AbstractOrbit{I, <:Integer}, pt::I, op=^) where I<:Integer
    gens_inv = inv.(gens)
    word = Word(gens_inv, orb, pt, op)
    g = word(gens_inv)
    return inv(g)
end


text:
using PCAs
using Test

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


text:
mockConfig = GraffConfig("http://mock", "9000", "", "", "")
mockConfig.sessionId = "TestSession"
mockConfig.robotId = "TestRobot"

@testset "Robot API" begin
    # Arrange
    mockResponse = HTTP.Response(200, "{\"id\": \"TestSession\",\"description\": \"\",\"robotId\": \"\",\"userId\": \"\",\"initialPoseType\": \"\",\"nodeCount\": 0,\"shouldInitialize\": \"true\",\"createdTimestamp\": \"\",\"lastUpdatedTimestamp\": \"\", \"lastSolvedTimestamp\": \"\", \"isSolverEnabled\": 1, \"links\": {\"self\": \"https://api.synchronysandbox.com/mocking\"}}")
    mockListResponse = HTTP.Response(200, "{\"sessions\": [{\"id\": \"\", \"links\": {\"self\": \"https://api.synchronysandbox.com/mocking\"}}], \"links\": {\"self\": \"https://api.synchronysandbox.com/mocking\"}}")
    mockErrorResponse = HTTP.Response(403)
    mockRequest = SessionDetailsRequest("mockId", "desc", "Pose2", true)
    robotConfig = Dict{String, String}()

    # Generate our mocks for the actual REST calls
    sendRequestMock = @patch _sendRestRequest(synchronyConfig::GraffConfig, verbFunction, url::String; data::String="", headers::Dict{String, String}=Dict{String, String}(), debug::Bool=false) = mockResponse
    sendRequestListMock = @patch _sendRestRequest(synchronyConfig::GraffConfig, verbFunction, url::String; data::String="", headers::Dict{String, String}=Dict{String, String}(), debug::Bool=false) = mockListResponse
    sendRequestErrorMock = @patch _sendRestRequest(synchronyConfig::GraffConfig, verbFunction, url::String; data::String="", headers::Dict{String, String}=Dict{String, String}(), debug::Bool=false) = mockErrorResponse

    sendRequestNodesMock = @patch _sendRestRequest(synchronyConfig::GraffConfig, verbFunction, url::String; data::String="", headers::Dict{String, String}=Dict{String, String}(), debug::Bool=false) = mockResponse

    # Success criteria
    apply(sendRequestListMock) do
        @testset "getSessions" begin
            # Act
            callResponse = getSessions("")
            # Assert
            @test length(callResponse.sessions) == 1
        end
        @testset "isSessionExisting" begin
            # Act
            callResponse = isSessionExisting("", "")
            # Assert
            @test callResponse == true
        end
    end
    apply(sendRequestMock) do
        @testset "getSession" begin
            # Act
            callResponse = getSession("", "")
            # Assert
            @test callResponse.id == "TestSession"
        end
        # @testset "updateSession") do
        #     # Act
        #     callResponse = addSession("", mockRequest)
        #     # Assert
        #     @fact callResponse.name --> "TestRobot"
        # end
        @testset "deleteSession" begin
            # Act & Assert
            callResponse = deleteSession("", "")
        end
        @testset "addSession" begin
            # Act
            callResponse = addSession("", mockRequest)
            # Assert
            @test callResponse.id == "TestSession"
        end
    end

    # apply(sendRequestMock) do
    #     @testset "getNodes") do
    #         # Act
    #         callResponse = getRobotConfig("GearsAD")
    #         # Assert
    #         @fact callResponse["name"] --> "TestRobot"
    #     end
    #     @testset "updateRobotConfig") do
    #         # Act
    #         callResponse = updateRobotConfig("GearsAD", robotConfig)
    #         # Assert
    #         @fact callResponse["name"] --> "TestRobot"
    #     end
    #
    # end

    # Failure mode
    # apply(sendRequestErrorMock) do
    #     @testset "getRobots - Failure Mode") do
    #         # Act & Assert
    #         @fact_throws ErrorException getRobots(    #     end
    #     @testset "getRobot - Failure Mode") do
    #         # Act & Assert
    #         @fact_throws ErrorException getRobot("TestRobot")
    #     end
    #     @testset "addRobot - Failure Mode") do
    #         # Act & Assert
    #         @fact_throws ErrorException addRobot(mockRequest)
    #     end
    #     @testset "deleteRobot - Failure Mode") do
    #         # Act & Assert
    #         @fact_throws ErrorException deleteRobot("GearsAD")
    #     end
    #     @testset "getRobotConfig - Failure Mode") do
    #         # Act & Assert
    #         @fact_throws ErrorException getRobotConfig("GearsAD")
    #     end
    #     @testset "updateRobotConfig - Failure Mode") do
    #         # Act & Assert
    #         @fact_throws ErrorException updateRobotConfig("GearsAD", robotConfig)
    #     end
    # end
end


text:
using Pkg; Pkg.activate(".")
using HeatEqCalc
sol = HeatEqCalc.calc_sol()
HeatEqCalc.save_sol(sol)

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public