text:
_lattice_primitives(::TYPE,angles,constants::NTuple{1}) where TYPE = (Point{TYPE}(constants[1]),)

Base.getindex(lat::Lattice{1},i::CartesianIndex) = getindex(lat,i[1])
Base.getindex(lat::Lattice{1},i) = map(i->lat[i],i)

Base.getindex(lat::Lattice{1,Bravais} ,i::Real)     = lat.origin + i*lat.a*lat.e1
Base.getindex(lat::Lattice{1,Cartesian} ,i::Real)   = lat.origin + i*lat.dx*lat.e1
Base.getindex(lat::Lattice{1,Polar}     ,i::Real)   = lat.origin + lat.r0 + i*lat.dr
Base.getindex(lat::Lattice{1,Spherical} ,i::Real)   = lat.origin + lat.r0 + i*lat.dr

latticeindex(lat::Lattice{1,Bravais}    ,p::Point{1}) = (p.x-lat.origin.x)/lat.a
latticeindex(lat::Lattice{1,Cartesian}  ,p::Point{1}) = (p.x-lat.origin.x)/lat.dx
latticeindex(lat::Lattice{1,Polar}      ,p::Point{1}) = (p.x-lat.origin.x)/lat.dr - 1/2
latticeindex(lat::Lattice{1,Spherical}  ,p::Point{1}) = (p.x-lat.origin.x)/lat.dr - 1/2

################################################################################
# Bravais
function Base.getproperty(lat::Lattice{1,Bravais}, sym::Symbol)
    if sym == :a
        return getfield(lat,:constants)[1]
    elseif sym == :x0
        return getfield(lat,:origin)[1]
    elseif sym == :e1
        return getfield(lat,:e)[1]
    else
        return getfield(lat,sym)
    end
end

function Base.propertynames(::Lattice{1,Bravais}, private=false)
    if private
        return fieldnames(Lattice)
    else
        return (:a, :x0)
    end
end

################################################################################
# Cartesian
function Base.getproperty(lat::Lattice{1,Cartesian}, sym::Symbol)
    if sym == :dx
        return getfield(lat,:constants)[1]
    elseif sym == :x0
        return getfield(lat,:origin)[1]
    elseif sym == :e1
        return getfield(lat,:e)[1]
    else
        return getfield(lat,sym)
    end
end

function Base.propertynames(::Lattice{1,Cartesian},private=false)
    if private
        return fieldnames(Lattice)
    else
        return (:dx, :origin, :x0)
    end
end

################################################################################
# Polar
function Base.getproperty(lat::Lattice{1,Polar}, sym::Symbol)
    if sym == :dr
        return getfield(lat,:constants)[1]
    elseif sym == :r0
        return getfield(lat,:constants)[1]/2
    else
        return getfield(lat,sym)
    end
end

function Base.propertynames(::Lattice{1,Polar},private=false)
    if private
        return fieldnames(Lattice)
    else
        return (:dr, :r0)
    end
end

################################################################################
# Spherical
function Base.getproperty(lat::Lattice{1,Spherical}, sym::Symbol)
    if sym == :dr
        return getfield(lat,:constants)[1]
    elseif sym == :r0
        return getfield(lat,:constants)[1]/2
    else
        return getfield(lat,sym)
    end
end

function Base.propertynames(::Lattice{1,Spherical},private=false)
    if private
        return fieldnames(Lattice)
    else
        return (:dr, :r0)
    end
end

################################################################################
# PRETTY PRINTING
import ..PRINTED_COLOR_NUMBER
import ..PRINTED_COLOR_WARN
import ..PRINTED_COLOR_DARK

function Base.show(io::IO, lat::Lattice{1,TYPE}) where TYPE
    get(io,:tabbed2,false) ? print(io,"\t\t") : nothing
    if TYPE <: Bravais
        printstyled(io,"1D Bravais Lattice ",color=PRINTED_COLOR_DARK)
        print(io,"(a=")
        printstyled(io, fmt("3.4f",lat.a),color=PRINTED_COLOR_NUMBER)
    elseif TYPE <: Cartesian
        printstyled(io,"1D Cartesian Lattice ",color=PRINTED_COLOR_DARK)
        print(io,"(dx=")
        printstyled(io, fmt("3.4f",lat.dx),color=PRINTED_COLOR_NUMBER)
    elseif TYPE <: Polar
        printstyled(io,"1D Polar Lattice ",color=PRINTED_COLOR_DARK)
        print(io,"(dr=")
        printstyled(io, fmt("3.4f",lat.dr),color=PRINTED_COLOR_NUMBER)
    elseif TYPE <: Spherical
        printstyled(io,"1D Spherical Lattice ",color=PRINTED_COLOR_DARK)
        print(io,"(dr=")
        printstyled(io, fmt("3.4f",lat.dr),color=PRINTED_COLOR_NUMBER)
    end
    print(io, ", origin=")
    printstyled(io, fmt("3.4f",lat.origin[1]),color=PRINTED_COLOR_NUMBER)
    print(io,")")
end

################################################################################
# PLOTTING
@recipe function f(lat::Lattice{1},n1=-10,n2=10)
    ylims --> [-1,1]
    lat[n1:n2]
end


text:
using UnitFloat

const M = 300
const N = 100000

function bench_mult_ufloat(M::Int, N::Int)
    x = UFloat.(rand(Float32, (N, 1)))
    y = UFloat.(rand(Float32, (N, 1)))
    @elapsed for _ in 1:M, i in 1:N
        x[i] = x[i] * y[i]
    end
end

function bench_mult_float32(M::Int, N::Int)
    x = rand(Float32, (N, 1))
    y = rand(Float32, (N, 1))
    @elapsed for _ in 1:M, i in 1:N
        x[i] = x[i] * y[i]
    end
end

function bench_mult_float64(M::Int, N::Int)
    x = rand(Float64, (N, 1))
    y = rand(Float64, (N, 1))
    @elapsed for _ in 1:M, i in 1:N
        x[i] = x[i] * y[i]
    end
end

function bench_mult_bfloat(M::Int, N::Int)
    x = BigFloat.(rand(Float64, (N, 1)))
    y = BigFloat.(rand(Float64, (N, 1)))
    @elapsed for _ in 1:M, i in 1:N
        x[i] = x[i] * y[i]
    end
end

function bench_add_ufloat(M::Int, N::Int)
    x = UFloat.(rand(Float32, (N, 1)))
    y = UFloat.(rand(Float32, (N, 1)))
    @elapsed for _ in 1:M, i in 1:N
        x[i] = x[i] + y[i]
    end
end

function bench_add_float32(M::Int, N::Int)
    x = rand(Float32, (N, 1))
    y = rand(Float32, (N, 1))
    @elapsed for _ in 1:M, i in 1:N
        x[i] = x[i] + y[i]
    end
end

function bench_add_float64(M::Int, N::Int)
    x = rand(Float64, (N, 1))
    y = rand(Float64, (N, 1))
    @elapsed for _ in 1:M, i in 1:N
        x[i] = x[i] + y[i]
    end
end

function bench_add_bfloat(M::Int, N::Int)
    x = BigFloat.(rand(Float64, (N, 1)))
    y = BigFloat.(rand(Float64, (N, 1)))
    @elapsed for _ in 1:M, i in 1:N
        x[i] = x[i] + y[i]
    end
end

println("Results for multiplication ($M runs of $N):")
println(" - UFloat:   $(bench_mult_ufloat(M, N) / M)")
println(" - Float32:  $(bench_mult_float32(M, N) / M)")
println(" - Float64:  $(bench_mult_float64(M, N) / M)")
println(" - BigFloat: $(bench_mult_bfloat(M, N) / M)")

println("Results for addition ($M runs of $N):")
println(" - UFloat:   $(bench_add_ufloat(M, N) / M)")
println(" - Float32:  $(bench_add_float32(M, N) / M)")
println(" - Float64:  $(bench_add_float64(M, N) / M)")
println(" - BigFloat: $(bench_add_bfloat(M, N) / M)")

ext: jl


text:
using Graphs, SimpleWeightedGraphs
using Plots, GraphRecipes
using AlphaZero
using GraphNeuralNetworks
using Statistics
import Flux
import AlphaZero.GI

struct GameSpec <: GI.AbstractGameSpec
    gnnGraph
end

function randGraph(numVerticies::Int)
    graph = rand(Float32, numVerticies, numVerticies)
    foreach(enumerate(eachcol(graph))) do (idx, col)
        graph[idx, :] .= col
        graph[idx, idx] = 0
    end
    nodes = collect(1:numVerticies)
    sources = Vector(vec(repeat(nodes, 1, numVerticies - 1)'))
    targets = vcat(map(idx -> filter(val -> val != idx, nodes), nodes)...)
    weights = map(zip(sources, targets)) do (source, target)
        graph[source, target]
    end
    return GNNGraph(sources, targets, weights; ndata = (; x = ones(Float32, 1, numVerticies)))
end

# GameSpec() = GameSpec(randGraph(rand(collect(1:20))))
GameSpec() = GameSpec(randGraph(10))
GameSpec(size::Int) = GameSpec(randGraph(size))

mutable struct GameEnv <: GI.AbstractGameEnv
    gnnGraph
    maskedActions::Vector{Bool} # Masked actions and visitedVerticies can be derived from graph, but are included for clarity
    visitedVerticies::Vector{Int}
    finished::Bool
    specGraph
end
GI.spec(game::GameEnv) = GameSpec(game.specGraph)

function GI.init(spec::GameSpec)
    return GameEnv(spec.gnnGraph, trues(spec.gnnGraph.num_nodes), Vector{Int}(), false, spec.gnnGraph)
end

function GI.set_state!(game::GameEnv, state)
    game.gnnGraph = state.gnnGraph
    game.maskedActions = state.availableActions
    game.visitedVerticies = getPath(state.gnnGraph, game.visitedVerticies)
    any(game.maskedActions) || (game.finished = true)
    return
end

function getPath(gnnGraph, visitedVerticies)
    isempty(visitedVerticies) && return Vector{Int}()
    lastVertex = last(visitedVerticies)
    if isone(count(idx -> idx == lastVertex, gnnGraph.graph[1]))
        push!(visitedVerticies, gnnGraph.graph[2][findfirst(ind -> ind == lastVertex, gnnGraph.graph[1])])
    end
    # sources = gnnGraph.graph[1]
    # targets = gnnGraph.graph[2]
    # madeConnection(vertex) = isone(count(val -> val == vertex, sources))
    # isConnected(vertex) = isone(count(val -> val == vertex, targets))
    # startingVertex = filter(idx -> madeConnection(idx) && !isConnected(idx), unique(sources))
    # isnothing(startingVertex) && (startingVertex = 1)
    # path = [startingVertex]
    # foreach(path) do vertex
    #     if madeConnection(vertex)
    #         push!(path, targets[findfirst(idx -> idx == vertex, sources)])
    #     end
    # end
    return visitedVerticies
end

function Base.hash(gnn::GNNGraph, h::UInt64)
    hash(hash(gnn.graph[1]) + hash(gnn.graph[2]) + hash(gnn.graph[3]) + hash(gnn.ndata) + h)
end

GI.two_players(::GameSpec) = false
GI.actions(a::GameSpec) = collect(range(1, length = a.gnnGraph.num_nodes))
GI.clone(g::GameEnv) = GameEnv(g.gnnGraph, deepcopy(g.maskedActions), deepcopy(g.visitedVerticies), g.finished, g.specGraph)
GI.white_playing(::GameEnv) = true
GI.game_terminated(g::GameEnv) = g.finished
GI.available_actions(g::GameEnv) = collect(range(1, length = g.gnnGraph.num_nodes))[g.maskedActions]
GI.actions_mask(game::GameEnv) = game.maskedActions

function GI.current_state(g::GameEnv)
    return (gnnGraph = g.gnnGraph, availableActions = g.maskedActions)
    # modifiedAdacencyList = deepcopy(g.graph[3])
    # if length(g.visitedVerticies) > 1
    #     foreach(enumerate(g.visitedVerticies[2:end])) do (idx, node)
    #         prevNode = state.path[idx]
    #         # Values are set to zero because SimpleWeightedGraphs discards edges with weights of zero.
    #         modifiedAdacencyList[prevNode, :] .= 0
    #         modifiedAdacencyList[:, node] .= 0
    #         modifiedAdacencyList[node, prevNode] = 0
    #         modifiedAdacencyList[prevNode, node] = spec.fadjlist[prevNode, node]
    #     end
    # end
    # return (ndata = (data = modifiedAdacencyList, x = ones(2)), availableActions = g.maskedActions)
end

function GI.play!(g::GameEnv, vertex::Int)
    adjMatrix = adjacency_matrix(g.gnnGraph)

    maskedActions = deepcopy(g.maskedActions)
    maskedActions[vertex] = false
    if isempty(g.visitedVerticies)
        sources = g.gnnGraph.graph[1]
        targets = g.gnnGraph.graph[2]
        weights = g.gnnGraph.graph[3]
    else
        edgeLength = deepcopy(adjMatrix[last(g.visitedVerticies), vertex])
        adjMatrix[:, vertex] .= 0
        adjMatrix[last(g.visitedVerticies),:] .= 0
        adjMatrix[last(g.visitedVerticies), vertex] = edgeLength

        sources = Vector{Int}()
        targets = Vector{Int}()
        weights = Vector{Float32}()
        foreach(enumerate(eachrow(adjMatrix))) do (i, col)
            foreach(enumerate(col)) do (j, val)
                if !iszero(val)
                    push!(sources, i)
                    push!(targets, j)
                    push!(weights, val)
                end
            end
        end
    end
    graph = GNNGraph(sources, targets, weights; ndata = (; x = ones(Float32, 1, g.gnnGraph.num_nodes)))
    state = (gnnGraph = graph, availableActions = maskedActions)

    push!(g.visitedVerticies, vertex)
    GI.set_state!(g, state)
    return
end

function GI.state_dim(game_spec::GameSpec)
    return size(game_spec.gnnGraph.ndata.x)[1]
end

function GI.white_reward(g::GameEnv)
    isempty(g.visitedVerticies[1:end-1]) && (return 0.0)
    # sources = g.gnnGraph.graph[1]
    # indicies = findall(idx -> idx ∈ g.visitedVerticies[1:end-1], sources)
    return -1 * mean(adjacency_matrix(g.gnnGraph))
end

function GI.heuristic_value(g::GameEnv)
    return GI.white_reward(g)
end

function GI.render(g::GameEnv)
    display(graphplot(adjacency_matrix(g.gnnGraph); curves = false))
end

function GI.graph_state(spec::GameSpec, state)
    return state.gnnGraph
end

function GI.action_string(::GameSpec, a)
    return string(a)
end

function GI.read_state(game::GameSpec)
    nodes = collect(Base.OneTo(game.gnnGraph.num_nodes))
    maskedActions = trues(game.gnnGraph.num_nodes)
    foreach(nodes) do node
        maskedActions[node] = count(vert -> vert == node, game.gnnGraph.graph[2]) == 1 ? false : true
    end
    return (gnnGraph = game.gnnGraph, availableActions = maskedActions)
end

function GI.parse_action(game::GameSpec, input::String)
    try
        p = parse(Int, input)
        return 1 <= p <= game.gnnGraph.num_nodes ? p : nothing
    catch
        return nothing
    end
end

function convert_sample(
    gspec::AlphaZero.Examples.Tsp.GameSpec,
    wp::SamplesWeighingPolicy,
    e::AlphaZero.TrainingSample)

  if wp == CONSTANT_WEIGHT
    w = Float32[1]
  elseif wp == LOG_WEIGHT
    w = Float32[log2(e.n) + 1]
  else
    @assert wp == LINEAR_WEIGHT
    w = Float32[n]
  end
  x = GI.graph_state(gspec, e.s)
  a = GI.actions_mask(GI.init(gspec, e.s))
  p = zeros(size(a))
  p[a] = e.π
  v = [e.z]
  return (; w, x, a, p, v)
end

function AlphaZero.convert_samples(
  gspec::AlphaZero.Examples.Tsp.GameSpec,
  wp::SamplesWeighingPolicy,
  es::Vector{<:AlphaZero.TrainingSample})

ces = [convert_sample(gspec, wp, e) for e in es]
W = Flux.batch((e.w for e in ces))
X = Flux.batch((e.x for e in ces))
X = Flux.batch(Vector{typeof(X[1])}(X))
A = Flux.batch((e.a for e in ces))
P = Flux.batch((e.p for e in ces))
V = Flux.batch((e.v for e in ces))
function f32(arr)
  if typeof(arr) <: Matrix
    return convert(AbstractArray{Float32}, arr)
  else
    return arr
  end
end
return map(f32, (; W, X, A, P, V))
end


text:
import Base: +, -, *, &, |, xor, ==, ~

for op in (:+, :-, :&, :|, :xor, :(==))
    @eval begin
        function ($op)(a::PyenumValues{T}, b::PyenumValues{S}) where {T<:Integer,S<:Integer}
            N = promote_type(T, S)
            ($op)(N(a.value), N(b.value))
        end
        function ($op)(a::PyenumValues{T}, b::PyenumValues{S}) where {T<:String,S<:String}
            print("convert")
            N = promote_type(T, S)
            ($op)(N(a.value), N(b.value))
        end
        function ($op)(a::PyenumValues{T}, b::S) where {T<:Integer,S<:Integer}
            N = promote_type(T, S)
            ($op)(N(a.value), N(b))
        end
        function ($op)(a::T, b::PyenumValues{S}) where {T<:Integer,S<:Integer}
            N = promote_type(T, S)
            ($op)(N(a), N(b.value))
        end
    end
end

function ~(a::PyenumValues{T}) where {T<:Integer}
    ~(T(a.value))
end

function ~(a::PyenumValues{T}) where {T<:String}
    ~(T(a.value))
end

Base.convert(::Type{T1}, x::PyenumValues{T2}) where {T1<:Integer,T2<:Integer} = convert(T1, T2(x.value))
Base.convert(::Type{T1}, x::PyenumValues{T2}) where {T1<:String,T2<:String} = convert(T1, T2(x.value))

max_stars_repo_name: MiguelMarcelino/PyEnum

max_issues_repo_name: MiguelMarcelino/PyEnum

max_forks_repo_name: MiguelMarcelino/PyEnum


text:
include("CenterOfMass.jl")
include("setup_parameters.jl")

plt[:close]("all")

actuator_limit = 3.0
max_robots = 4
csqmi_ratio = 10.0

prior = initialize_prior(circle_ps, resolution, interior_q, masses)

mass = min_mass + rand() * (max_mass - min_mass)
com_p = r_attach * rand_in_circle()
gravity_w = mass*g*point_to_wrench(com_p)
theta = [com_p[1:2]; mass]
println("\nCOM: $(com_p), mass: $(mass)")

robots = Array{Int64,1}()

n_measurement = 10

belief = deepcopy(prior)

for ii = 1:n_measurement
  figure()
  plot_attachment_points(attachment_ps)
  plot_occupied_points(attachment_ps[robots])

  cloud = to_cloud(belief, interior_q)
  scaled_ps = (15*cloud[4,:]'/maximum(cloud[4,:]))
  scatter3D(cloud[1,:]', cloud[2,:]', cloud[3,:]', "z", scaled_ps.^2, "purple", alpha = 0.5)

  scatter3D([com_p[1]], [com_p[2]], [mass], "z", 200, "red", marker="*", alpha=1)
  xlabel("X")
  ylabel("Y")
  zlabel("M")
  title("Iteration: $(ii)")

  println("\nMeasurement: $(ii), robots: $(robots)\n")

  @time @show csqmi_current, robots_current = maximize_csqmi_available_robots(
    robots, belief, circle_ws, attachment_ws, sigma, interior_q, actuator_limit)

  csqmi = csqmi_current
  robots_measurement = robots_current

  if length(robots) < max_robots
    @time @show csqmi_additional, robots_additional, additional_robot =
      maximize_csqmi_additional_robot(robots, belief,
        circle_ws, attachment_ws, sigma, interior_q, actuator_limit)

    if length(robots) ==0 || csqmi_additional / csqmi_current > csqmi_ratio
      println("\nAdding robot $(additional_robot), csqmi: $(csqmi_additional) vs: $(csqmi_current)\n")

      csqmi = csqmi_additional
      robots_measurement = robots_additional
      push!(robots, additional_robot)

      plot_new_point(attachment_ps[additional_robot])
    end
  end

  plot_measurement_points(attachment_ps[robots_measurement])

  applied_w = mean(attachment_ws[robots_measurement])
  total_actuator_limit = length(robots_measurement) * actuator_limit

  applied_f = critical_force_iterative(applied_w, gravity_w,
    circle_ws, total_actuator_limit)

  f_hat = applied_f + sigma * randn()

  println("\nRobots measurement: $(robots_measurement)")
  println("Measurement: $(applied_f)/$(f_hat), total limit: $(total_actuator_limit)\n")

  critical_forces = get_critical_values(circle_ws, applied_w, belief,
    interior_q, total_actuator_limit)

  update_prior!(belief, critical_forces, f_hat, sigma)

  estimate = weighted_average(belief)

  error = normalized_error(estimate - theta)
  println("Estimate: $(estimate), actual: $(theta), normalized error: $(error)")
end

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
using Weave

using Pkg
if haskey(ENV, "GITLAB_CI")
    Pkg.add([PackageSpec(name = x; rev = "master")
             for x in ["CUDAapi", "GPUArrays", "CUDAnative", "NNlib", "CUDAdrv"]])
end

cd(joinpath(@__DIR__, "src")) do
    rm("../build"; force=true, recursive=true)

    # intro tutorial
    weave("intro.jl", out_path="../build", doctype="md2html")
    cp("intro1.png", "../build/intro1.png")
end


text:
import DifferentialEquations;	const DE = DifferentialEquations
"""
	dOmega!(dOmega, Omega, args, t)
Differential equation for the first terms of the Magnus expansion.
#Arguments
- dOmega: 2-dimensional array with the time derivative of the Magnus expansion
terms.
- Omega: 2-dimensional array with the Magnus expansion terms. The nth column is
the nth term of the Magnus expansion: Omega[:,n] = Omega_n.
- args: Array with the many different arguments of the function.
args[1] = order: The number of terms of the Magnus expansion that will be computed.\n
args[2] = A: the matrix of function that defines the linear ODE (dx/dt = A*x).\n
args[3] = S: array of arrays used to store the S operators. ``S[n] = S_n``, where
S[n] is a 2-dimensional array. ``(S[n])[:,j] = S_n^{(j)}``.\n
args[4] = B: array where B[n] = b_n/n!. b_n is the n'th Bernoulli number.\n
args[5] = Snj: array to temporarily store S_n^{(j)}\n
args[6] = temp1 and args[7] = temp2: temporary arrays to store data.\n
- t: time variable.\n
"""
function dOmega!(dOmega, Omega, args, t)
	order = args[1]::Integer
	A = args[2]
	comm! = args[3]
	S = args[4]
	B = args[5]
	Snj = args[6]
	temp = args[7]
	temp2 = args[8]

	dOmega[:, 1] .= A(t)

	S2 = S[2]
	Omega1 = view(Omega, :, 1)
	dOmega1 = view(dOmega, :, 1)
	comm!(Snj, Omega1, dOmega1)
	S2[:, 1] .= Snj
	dOmega[:, 2] .= B[1].*Snj
	for n = 3: order
		Sn = S[n]
		update_S!(Snj, S, n, 1, comm!, Omega, Omega1, dOmega1, temp2)
		temp .= B[1].*Snj
		for j = 2: n - 1
			update_S!(Snj, S, n, j, comm!, Omega, Omega1, dOmega1, temp2)
			if j%2 == 0
				temp .= temp .+ B[j].*Snj
			end
		end
		dOmega[:, n] .= temp
	end
	return nothing
end
"""
	update_S!(Snj, S, n, j, Omega, Omega1, dOmega1, temp2)
Calculate ``S_n^{(j)}`` recursivelly and store the result in S and Snj. The result
is stored in Snj just for efficiency purposes.
#Arguments
See dOmega! documentation.
"""
function update_S!(Snj, S, n, j, comm!, Omega, Omega1, dOmega1, temp2)
	Sn = S[n]
	if j == 1
		Omegam = view(Omega, :, n - 1)
		comm!(Snj, Omegam, dOmega1)
	else
		comm!(Snj, Omega1, view(S[n - 1], :, j - 1))
		for m = 2: n - j
			Omegam = view(Omega, :, m)
			comm!(temp2, Omegam, view(S[n - m], :, j - 1))
			Snj .= Snj .+ temp2
		end
	end
	Sn[:, j] .= Snj
	return nothing
end
"""
	bernoulli(n)
Calculate the n'th Bernoulli number. This function was taken from Rosetta Code
(https://rosettacode.org/wiki/Bernoulli_numbers#Julia).
#Arguments:
- n: the index of the Bernoulli number.
"""
function bernoulli(n)
    A = Vector{Rational{Int64}}(undef, n + 1)
    for m = 0: n
        A[m + 1] = 1//(m + 1)
        for j = m: -1: 1
            A[j] = j*(A[j] - A[j + 1])
        end
    end
    return A[1]
end
"""
	solveOmega(tf, order, A, numdims)
Calculate the Magnus expansion.
#Arguments
- tf: the time at which the Magnus expansion is calculated. It is assumed that
the time evolution goes from 0 to tf.
- order: the number of calculated terms in the Magnus expansion.
- A: the matrix of function that defines the linear ODE (dx/dt = A*x).
- numdims: the number of dimensions of the ODE.
"""
function solveOmega(tf, order, A, comm!, numdims)
	order_aux = order
	order = max(2, order)

	numtype = typeof(tf)
	S = createS(order, numdims, numtype)
	B = zeros(Rational{Int64}, order)
	for j = 1: order
		B[j] = bernoulli(j)//factorial(j)
	end
	B = convert(Array{numtype}, B)
	B = Tuple(B)

	Omega0 = zeros(numtype, numdims, order)
	tspan = (zero(tf), tf)
	temp1 = zeros(numtype, numdims)
	temp2 = zeros(numtype, numdims)
	temp3 = zeros(numtype, numdims)
	args = (order, A, comm!, S, B, temp1, temp2, temp3)
	prob = DE.ODEProblem(dOmega!, Omega0, tspan, args)
	sol = DE.solve(prob, DE.Vern9(), saveat=[tf], reltol=1e-11, abstol=1e-14,
				   maxiters = 10^6)

	if sol.retcode == :Success
		return sol(tf)[:, 1:order_aux]
	else
		error("Calculation of the Magnus series was unsuccessful.")
	end
end
"""
	createS(order, numdims, numtype)
Create the S array, which is an argument of the dOmega! function.
#Arguments:
-numtype: The type of the numbers used in problem.
See solveOmega documentation.
"""
function createS(order, numdims, numtype)
	S = [zeros(numtype, numdims, 1)]
	for n=2:order
		push!(S, zeros(numtype, numdims, n-1))
	end
	return S
end

### TEST ###

max_stars_repo_name: thalesfr/MagnusCorrection

max_issues_repo_name: thalesfr/MagnusCorrection

max_forks_repo_name: thalesfr/MagnusCorrection


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# ReverseMcCormick.jl
# A Reverse McCormick operator library in Julia
# See https://github.com/PSORLab/ReverseMcCormick.jl
#############################################################################
# src/reverse_operators/exponential.jl
# Contains definitions of reverse McCormick operators of exp, exp2, exp10, expm1
# log, log2, log10, log1p
#############################################################################

"""
$(SIGNATURES)

Reverse McCormick operator for `exp`.
"""
function exp_rev(y::MC, x::MC)
    isempty(y) && (return y, y)
    y = y ∩ Interval{Float64}(0.0, Inf)
    x = x ∩ log(y)
    y, x
end

"""
$(SIGNATURES)

Reverse McCormick operator for `exp2`.
"""
function exp2_rev(y::MC, x::MC)
    isempty(y) && (return y, y)
    y = y ∩ Interval{Float64}(0.0, Inf)
    x = x ∩ log2(y)
    y, x
end

"""
$(SIGNATURES)

Reverse McCormick operator for `exp10`.
"""
function exp10_rev(y::MC, x::MC)
    isempty(y) && (return y, y)
    y = y ∩ Interval{Float64}(0.0, Inf)
    x = x ∩ log10(y)
    y, x
end

"""
$(SIGNATURES)

Reverse McCormick operator for `expm1`.
"""
function expm1_rev(y::MC, x::MC)
    isempty(y) && (return y, y)
    y = y ∩ Interval{Float64}(-1.0, Inf)
    x = x ∩ log1p(y)
    y, x
end

"""
$(SIGNATURES)

Reverse McCormick operator for `log`.
"""
function log_rev(y::MC, x::MC)
    isempty(y) && (return y, y)
    x = x ∩ exp(y)
    y, x
end

"""
$(SIGNATURES)

Reverse McCormick operator for `log2`.
"""
function log2_rev(y::MC, x::MC)
    isempty(y) && (return y, y)
    x = x ∩ exp2(y)
    y, x
end

"""
$(SIGNATURES)

Reverse McCormick operator for `log10`.
"""
function log10_rev(y::MC, x::MC)
    isempty(y) && (return y, y)
    x = x ∩ exp10(y)
    y, x
end

"""
$(SIGNATURES)

Reverse McCormick operator for `log1p`.
"""
function log1p_rev(y::MC, x::MC)
    isempty(y) && (return y, y)
    x = x ∩ expm1(y)
    y, x
end


text:
# demonstrates the use of pch=0 for pixels

using SixelGraphics

n = 128
x = 0.1*ones(n)
for i=2:n
  x[i] = 4*x[i-1]*(1-x[i-1])
end

y = x[2:n]
x = x[1:n-1]

sixelplot(x,y, title="Logistic Map", xlab="x[i]", ylab="x[i+1]",
          xlim=[0,1], ylim=[0,1], typ='p', pch=0)


text:
using Documenter, MIPVerify

makedocs(
    modules = [MIPVerify],
    format = :html,
    sitename = "MIPVerify.jl",
    authors = "Vincent Tjeng and contributors.",
    pages = [
        "Home" => "index.md",
        "Tutorials" => "tutorials.md",
        "Finding Adversarial Examples" => [
            "finding_adversarial_examples/single_image.md",
            "finding_adversarial_examples/batch_processing.md"
        ],
        "Importing" => [
            "utils/import_example_nets.md",
            "utils/import_weights.md",
            "utils/import_datasets.md",
        ],
        "Working with Neural Net Parameters" => [
            "net_components/overview.md",
            "net_components/layers.md",
            "net_components/nets.md",
            "net_components/core_ops.md"
        ]
    ],
checkdocs = :exports
)

deploydocs(
    deps = nothing,
    repo = "github.com/vtjeng/MIPVerify.jl.git",
    target = "build",
    make = nothing,
    julia = "0.6",
    osname = "linux"
)


text:
struct AcquisitionBucketStats
    kspace_encode_step_1::Set{UInt16}
    kspace_encode_step_2::Set{UInt16}
    contrast::Set{UInt16}
    slice::Set{UInt16}
    phase::Set{UInt16}
    repetition::Set{UInt16}
    segment::Set{UInt16}
    average::Set{UInt16}
    set::Set{UInt16}
end

struct AcquisitionBucket

    data::Vector{MRD.Acquisition}
    datastats::Vector{AcquisitionBucketStats}
    ref::Vector{MRD.Acquisition}
    refstats::Vector{AcquisitionBucketStats}
    waveforms::Vector{MRD.Waveform}
end

module AcqBuckInternal

using Gadgetron.MRD
struct bundle_meta
    count::UInt64
    header_bytes::UInt64
    data_bytes::UInt64
    trajectory_bytes::UInt64
end

struct waveform_meta
    count::UInt64
    header_bytes::UInt64
    data_bytes::UInt64
end

struct bucket_meta
    data::bundle_meta
    reference::bundle_meta
    data_stats::UInt64
    reference_stats::UInt64
    waveforms::waveform_meta
end

function read_data_as_array(io::IO, data_type::Type{T},shape) where {T}
    res = Matrix{T}(undef, shape...)
    Base.read!(io,res)
    return res
end

function read_waveforms(io::IO, meta::waveform_meta)
    headers = [MRD.read(io,MRD.RawWaveformHeader) for i in 1:meta.count]
    data_arrays = [ read_data_as_array(io,UInt32,(header.number_of_samples, header.channels)) for header in headers ]

    return [MRD.Waveform(MRD.WaveformHeader(header),data) for (header,data) in zip(headers,data_arrays)]
end

function read_acquisitions(io::IO, meta::bundle_meta)
    headers = [MRD.read(io,MRD.RawAcquisitionHeader) for i in 1:meta.count]

    trajectories = [read_data_as_array(io,Float32,(head.number_of_samples,head.trajectory_dimensions)) for head in headers]

    acqs = [read_data_as_array(io,ComplexF32,(head.number_of_samples,head.active_channels)) for head in headers]

    return [MRD.Acquisition(MRD.AcquisitionHeader(header),data,trajectory) for (header,data,trajectory) in zip(headers,acqs,trajectories)]
end
end

MRD.read(io::IO, ::Type{AcquisitionBucketStats}) = MRD.read_struct(io, AcquisitionBucketStats)

function MRD.read(io::IO, ::Type{AcquisitionBucket})
    meta = MRD.unsafe_read(io, AcqBuckInternal.bucket_meta)

    AcquisitionBucket(
        AcqBuckInternal.read_acquisitions(io,meta.data),
        MRD.read(io,Vector{AcquisitionBucketStats}),
        AcqBuckInternal.read_acquisitions(io,meta.reference),
        MRD.read(io, Vector{AcquisitionBucketStats}),
        AcqBuckInternal.read_waveforms(io,meta.waveforms)
    )
end


text:
# CGT97
# Corless, R. M.; Gianni, P. M. & Trager, B. M.
# A reordered Schur factorization method for zero-dimensional polynomial systems with multiple roots
# Proceedings of the 1997 international symposium on Symbolic and algebraic computation, 1997, 133-140

# MD95
# Manocha, D. & Demmel, J. Algorithms for intersecting parametric and algebraic curves II: multiple intersections
# Graphical Models and Image Processing, Elsevier, 1995, 57, 81-100

using LinearAlgebra # for I
@testset "Section 4.1 MD95" begin
    η = 1e-10
    A = [zeros(10, 1) Matrix(I, 10, 10); zeros(1, 10) 0.5]
    A[10, 11] = 0
    A[10, 1] = η
    @test sort.(SemialgebraicSets.clusterordschur(A, sqrt(eps(Float64)))[2]) == [[1, 2, 3, 4, 5, 6, 7, 8, 9, 10], [11]]
end

@testset "Example 4.1 MD95" begin
    A = [  0  0  0  1  0  0;
           0  0  0  0  1  0;
           0  0  0  0  0  1;
          -1  0  1 -1  0  0;
          -1  0  1 -2 -1  0;
          -1  0  1 -2 -2 -1]
    @test sort.(SemialgebraicSets.clusterordschur(A, sqrt(eps(Float64)))[2]) == [[2], [1, 5, 6]]
end

function testelements(X, Y; atol=Base.rtoldefault(Float64), kwargs...)
    @test length(X) == length(Y)
    for y in Y
        @test any(x -> isapprox(x, y; atol=atol, kwargs...), X)
    end
end
function testelementstypes(X, Y; kwargs...)
    testelements(X, Y; kwargs...)
    for T in [Rational{Int}, Float64]
        if X isa FixedVariablesSet
            U = T
        else
            U = float(T)
        end
        V = MultivariatePolynomials.changecoefficienttype(X, T)
        testelements(V, Y; kwargs...)
        @test eltype(V) == Vector{U}
        @test collect(V) isa Vector{Vector{U}}
    end
end

# We use a fixed RNG in the tests to decrease nondeterminism. There is still nondeterminism in LAPACK though
using Random
solver = ReorderedSchurMultiplicationMatricesSolver(sqrt(eps(Float64)), MersenneTwister(0))

@testset "Zero-dimensional ideal" begin
    Mod.@polyvar x y z
    V = @set x == y
    @test !iszerodimensional(V)
    @static if VERSION >= v"0.7-"
        @test_throws ErrorException iterate(V)
    else
        @test_throws ErrorException start(V)
    end
    @test_throws ErrorException length(V)
    V = @set 4x^2 == -5x && 3x^3 == 0 solver
    @test V.solver.solver === solver
    @test iszerodimensional(V)
    testelementstypes(V, [[0]])
    V = @set y == x^2 && z == x^3 solver
    @test !iszerodimensional(V)
    V = @set x^3 == 2x*y &&  x^2*y == 2y^2 + x solver
    @test iszerodimensional(V)
    testelementstypes(V, [[0, 0]])
    V = @set x == 1 solver
    @test iszerodimensional(V)
    testelementstypes(V, [[1]])
    V = @set x == 1 && y == 2 solver
    @test iszerodimensional(V)
    testelementstypes(V, [[1, 2]])
    V = @set x == 4 && y^2 == x solver
    @test iszerodimensional(V)
    testelementstypes(V, [[4, 2], [4, -2]])
    V = @set x^2 + x == 6 && y == x+1 solver
    @test iszerodimensional(V)
    testelements(V, [[2, 3], [-3, -2]])
    V = @set x^2 + x == 6 && y^2 == x solver
    @test iszerodimensional(V)
    testelements(V, [[2, √2], [2, -√2]])
end

@testset "Projective zero-dimensional ideal" begin
    Mod.@polyvar x y z
    V = projectivealgebraicset([x - y], solver)
    @test iszerodimensional(V)
    testelementstypes(V, [[1, 1]])
    V = @set x + y == z solver
    V.projective = true
    @test !iszerodimensional(V)
    V = @set y == 2x solver
    V.projective = true
    @test iszerodimensional(V)
    testelementstypes(V, [[1, 2]])
    V = @set x + y == y solver
    V.projective = true
    @test iszerodimensional(V)
    testelementstypes(V, [[0, 1]])
    V = projectivealgebraicset([x + y - x])
    @test iszerodimensional(V)
    testelementstypes(V, [[1, 0]])
end

@testset "Example 5.1 of CGT97" begin
    ɛ = 1e-4
    Iɛ = [1 - ɛ 0
          0     1 + ɛ]
    J = [0 1
         1 0]
    Z = zeros(2, 2)
    A = [Iɛ Z
         Z  J]
    B = [J Z
         Z Iɛ]
    α = 0.219
    testelements(SemialgebraicSets._solvemultiplicationmatrices([A, B], [α, 1-α], ReorderedSchurMultiplicationMatricesSolver{Float64}()), [[1.0, -1.0], [1.0, 1.0], [-1.0, 1.0]]; rtol=1e-7)
end

@testset "Example 4.3 of MD95" begin
    Mod.@polyvar x y
    V = @set x^2 + 4y^4 == 4y^2 && 3x^4 + y^2 == 5x^3 solver
    # This test is tricky because in the schur decomposition, the 4 last eigenvalues are e.g. 3.4e-7, -1.7e-7+3e-7im, -1.7e-7-3e-7im, -6e-16
    # the second and third do not seem that close but when the three first are averaged it is very close to zero.
    @test iszerodimensional(V)
    testelementstypes(V, [[0.66209555, 0.935259169], [0.66209555, -0.935259169], [0.0516329456, -0.025825086], [0.0516329456, 0.025825086], [0, 0]])
end

@testset "Example 5.2 of CGT97" begin
    Mod.@polyvar x y z
    V = @set x^2 + y^2 == 1 && x^3 + (2 + z)*x*y + y^3 == 1 && z^2 == 2 solver
    @test iszerodimensional(V)
    iszd, B = monomialbasis(V.I)
    @test iszd
    @test B == [y^3*z, x*y*z, y^3, y^2*z, x*y, x*z, y^2, y*z, x, y, z, 1]
    testelements(V, [[0, 1, √2], [0, 1, -√2], [1, 0, -√2], [1, 0, √2], [-√2/2, -√2/2, √2], [√2/2, √2/2, -√2]])
end

#@testset "Example 4.4 of MD95 and 5.3 of CGT97" begin
#    Mod.@polyvar x y
#    F = -2-7x+14x^3-7x^5+x^7 + (7-42x^2+35x^4-7x^6)*y + (16+42x-70x^3+21x^5)*y^2 + (-14+70x^2-35x^4)*y^3 + (-20-35x+35x^3)*y^4 + (7-21x^2)*y^5 + (8+7x)*y^6 - y^7 - y^8
#    G = differentiate(F, y)
#    @show G
#    V = @set F == 0 && G == 0 solver
#    @test iszerodimensional(V)
#    sols = [[-3.91298, -1.95065], [-3.23984, -1.56367], [-3.21615, -1.41421], [-3.09474, -1.84776], [3.6497, 1.84776], [3.65578, 1.80399], [2.2928, 1.84776], [2.66119, 1.41421], [2.68379, 1.23369], [2.5673, 0.765367], [0.600779, 1.84776], [0.969172, 1.41421], [-1.85925, -1.41422], [-1.40272, -1.84775], [-2.01312, -0.812102], [-2.01235, -0.765367], [1.26106, 0.265359], [-0.481613, 0.765367], [1.03657, -0.765367], [-0.483778, 0.630692], [-0.0458213, -1.84776], [-1.80194, 0.]]
#    @show length(sols)
#    #testelements(V, sols)
#    @show length(V)
#    for v in collect(V)
#        @show v
#        @show F(variables(F) => v)
#        @show G(variables(G) => v)
#    end
#end

@testset "Complex not yet implemented" begin
    Mod.@polyvar x
    V = @set x^2 == im
    @test_throws MethodError collect(V)
end


text:
# SPDX-License-Identifier: MIT

module FSMs

export AbstractFSM
export states
export arcs
export semiring

export AbstractMutableFSM
export addstate!
export addarc!

export State
export isinit
export isfinal

export Arc

export VectorFSM
export HierarchicalFSM
export MatrixFSM
export UnionMatrixFSM

export gpu
export determinize
export minimize
export renormalize
export transpose
export remove_label

using LinearAlgebra
using SparseArrays
using CUDA
using CUDA.CUSPARSE

using ..Semirings

include("abstractfsm.jl")
include("vectorfsm.jl")
include("hierarchicalfsm.jl")
include("matrixfsm.jl")
include("fsmop.jl")

end

max_stars_repo_name: MartinKocour/HiddenMarkovModel

max_issues_repo_name: MartinKocour/HiddenMarkovModel

max_forks_repo_name: BUTSpeechFIT/ProbabilisticFSMs


text:
using Flux
using WindMixing
using JLD2
using FileIO
using OceanParameterizations
using OrdinaryDiffEq
using Random
using GalacticOptim
using LinearAlgebra

BLAS.set_num_threads(1)

# train_files = [
#     # "wind_-5e-4_cooling_4e-8",
#     # "wind_-1e-3_cooling_4e-8",
#     # "wind_-2e-4_cooling_1e-8",
#     # "wind_-1e-3_cooling_2e-8",
#     # "wind_-5e-4_cooling_1e-8",
#     # "wind_-2e-4_cooling_5e-8",
#     # "wind_-5e-4_cooling_3e-8",
#     "wind_-2e-4_cooling_3e-8",
#     # "wind_-1e-3_cooling_3e-8",
#     # "wind_-1e-3_heating_-4e-8",
#     # "wind_-1e-3_heating_-1e-8",
#     # "wind_-1e-3_heating_-3e-8",
#     # "wind_-5e-4_heating_-5e-8",
#     # "wind_-5e-4_heating_-3e-8",
#     # "wind_-5e-4_heating_-1e-8",
#     # "wind_-2e-4_heating_-5e-8",
#     # "wind_-2e-4_heating_-3e-8",
#     # "wind_-2e-4_heating_-1e-8",
# ]

train_files = ["wind_-2e-4_diurnal_3.5e-8"]

PATH = pwd()

OUTPUT_PATH = joinpath(PATH, "training_output")
OUTPUT_PATH = "D:\\University Matters\\MIT\\CLiMA Project\\OceanParameterizations.jl\\training_output"

VIDEO_PATH = joinpath(PATH, "Output")

EXTRACTED_OUTPUT_PATH = joinpath(PATH, "extracted_training_output")
EXTRACTED_OUTPUT_PATH = OUTPUT_PATH

FILE_NAME = "NDE_18sim_windcooling_windheating_18sim5paramsBFGS_divide1f5_gradient_smallNN_mish_scale_5e-3_rate_1e-4"
FILE_PATH = joinpath(OUTPUT_PATH, "$(FILE_NAME).jld2")

EXTRACTED_FILE_PATH = joinpath(EXTRACTED_OUTPUT_PATH, "$(FILE_NAME)_extracted.jld2")
# @assert !isfile(FILE_PATH)

# PARAMETERS_PATH = joinpath(EXTRACTED_OUTPUT_PATH, "parameter_optimisation_18sim_windcooling_windheating_5params_BFGS_extracted.jld2")

# parameters_file = jldopen(PARAMETERS_PATH)
# mpp_parameters = parameters_file["parameters"]
# close(parameters_file)

# ν₀_initial = 1f-4
# ν₋_initial = 1f-1
# ΔRi_initial = 1f-1
# Riᶜ_initial = 0.25f0
# Pr_initial = 1f0

# mpp_scalings = 1 ./ [ν₀_initial, ν₋_initial, ΔRi_initial, Riᶜ_initial, Pr_initial]

# ν₀, ν₋, ΔRi, Riᶜ, Pr = mpp_parameters ./ mpp_scalings

ν₀ = 1f-4
ν₋ = 1f-1
ΔRi = 1f-1
Riᶜ = 0.25f0
Pr = 1f0

# FILE_PATH_uw = joinpath(PATH, "extracted_training_output", "uw_NN_training_1sim_-1e-3_extracted.jld2")
# FILE_PATH_vw = joinpath(PATH, "extracted_training_output", "vw_NN_training_1sim_-1e-3_extracted.jld2")
# FILE_PATH_wT = joinpath(PATH, "extracted_training_output", "wT_NN_training_1sim_-1e-3_extracted.jld2")

# uw_file = jldopen(FILE_PATH_uw, "r")
# vw_file = jldopen(FILE_PATH_vw, "r")
# wT_file = jldopen(FILE_PATH_wT, "r")

# uw_NN = uw_file["neural_network"]
# vw_NN = vw_file["neural_network"]
# wT_NN = wT_file["neural_network"]

# FILE_PATH_NN = joinpath(PATH, "extracted_training_output",
#                         "NDE_training_mpp_8sim_wind_mixing_cooling_diffusivity_1e-1_Ri_1e-1_weights_divide1f5_gradient_smallNN_scale_5e-3_rate_1e-4_extracted.jld2")

# @assert isfile(FILE_PATH_NN)
# file = jldopen(FILE_PATH_NN, "r")

# uw_NN = file["neural_network/uw"]
# vw_NN = file["neural_network/vw"]
# wT_NN = file["neural_network/wT"]

N_inputs = 96
hidden_units = 400
N_outputs = 31

# weights, re = Flux.destructure(Chain(Dense(N_inputs, 32, swish), Dense(32, 400, swish), Dense(400, N_outputs)))
# weights, re = Flux.destructure(Chain(Dense(N_inputs, hidden_units, swish), Dense(hidden_units, hidden_units, swish), Dense(hidden_units, N_outputs)))
weights, re = Flux.destructure(Chain(Dense(N_inputs, 50, mish), Dense(50, 20, mish), Dense(20, 31)))

uw_NN = re(weights ./ 1f5)
vw_NN = re(weights ./ 1f5)
wT_NN = re(weights ./ 1f5)

gradient_scaling = 5f-3
training_fractions = (T=0.8f0, ∂T∂z=0.8f0, profile=0.5f0)

diurnal = occursin("diurnal", train_files[1])

train_parameters = Dict(
                               "ν₀" => ν₀,
                               "ν₋" => ν₋,
                              "ΔRi" => ΔRi,
                              "Riᶜ" => Riᶜ,
                               "Pr" => Pr,
                                "κ" => 10f0,
    "modified_pacanowski_philander" => true,
            "convective_adjustment" => false,
                   "smooth_profile" => false,
                        "smooth_NN" => false,
                        "smooth_Ri" => false,
                   "train_gradient" => true,
                     "zero_weights" => true,
                 "gradient_scaling" => gradient_scaling,
               "training_fractions" => training_fractions,
                          "diurnal" => diurnal,
    )

# train_epochs = [1]
# train_tranges = [1:9:1153]
# train_iterations = [600]
# train_optimizers = [[ADAM(1e-4)]]

train_epochs = [1]
train_tranges = [1:20:200]
train_iterations = [5]
train_optimizers = [[ADAM(3e-4)]]

timestepper = ROCK4()

function train(FILE_PATH, train_files, train_epochs, train_tranges, train_parameters, train_optimizers, train_iterations, uw_NN, vw_NN, wT_NN, timestepper)
    @info "Writing metadata"
    write_metadata_NDE_training(FILE_PATH, train_files, train_epochs, train_tranges, train_parameters, train_optimizers, uw_NN, vw_NN, wT_NN)
    for i in 1:length(train_epochs)
        @info "iteration $i/$(length(train_epochs)), time range $(train_tranges[i])"
        if train_parameters["modified_pacanowski_philander"]
            uw_NN, vw_NN, wT_NN = train_NDE(uw_NN, vw_NN, wT_NN, train_files, train_tranges[i], timestepper, train_optimizers[i], train_epochs[i], FILE_PATH, i,
                                                         maxiters = train_iterations[i],
                                    modified_pacanowski_philander = train_parameters["modified_pacanowski_philander"],
                                            convective_adjustment = train_parameters["convective_adjustment"],
                                                               ν₀ = train_parameters["ν₀"],
                                                               ν₋ = train_parameters["ν₋"],
                                                              ΔRi = train_parameters["ΔRi"],
                                                              Riᶜ = train_parameters["Riᶜ"],
                                                               Pr = train_parameters["Pr"],
                                                                κ = train_parameters["κ"],
                                                   smooth_profile = train_parameters["smooth_profile"],
                                                        smooth_NN = train_parameters["smooth_NN"],
                                                        smooth_Ri = train_parameters["smooth_Ri"],
                                                   train_gradient = train_parameters["train_gradient"],
                                                     zero_weights = train_parameters["zero_weights"],
                                                #  gradient_scaling = train_parameters["gradient_scaling"],
                                               training_fractions = train_parameters["training_fractions"],
                                                          diurnal = train_parameters["diurnal"],
                                    )
        else
            uw_NN, vw_NN, wT_NN = train_NDE(uw_NN, vw_NN, wT_NN, train_files, train_tranges[i], timestepper, train_optimizers[i], train_epochs[i], FILE_PATH, i,
                                                         maxiters = train_iterations[i],
                                    modified_pacanowski_philander = train_parameters["modified_pacanowski_philander"],
                                            convective_adjustment = train_parameters["convective_adjustment"],
                                                                κ = train_parameters["κ"],
                                                   smooth_profile = train_parameters["smooth_profile"],
                                                        smooth_NN = train_parameters["smooth_NN"],
                                                        smooth_Ri = train_parameters["smooth_Ri"],
                                                   train_gradient = train_parameters["train_gradient"],
                                                     zero_weights = train_parameters["zero_weights"],
                                                #  gradient_scaling = train_parameters["gradient_scaling"],
                                               training_fractions = train_parameters["training_fractions"],
                                                          diurnal = train_parameters["diurnal"],
                                                )
        end
    end
    return uw_NN, vw_NN, wT_NN
end

uw_NN_res, vw_NN_res, wT_NN_res = train(FILE_PATH, train_files, train_epochs, train_tranges, train_parameters, train_optimizers, train_iterations,
                                    uw_NN, vw_NN, wT_NN, timestepper)

# extract_NN(FILE_PATH, EXTRACTED_FILE_PATH, "NDE")

# test_files = [
#     # "wind_-5e-4_cooling_4e-8",
#     # "wind_-1e-3_cooling_4e-8",
#     # "wind_-2e-4_cooling_1e-8",
#     # "wind_-1e-3_cooling_2e-8",
#     # "wind_-5e-4_cooling_1e-8",
#     # "wind_-2e-4_cooling_5e-8",
#     # "wind_-5e-4_cooling_3e-8",
#     "wind_-2e-4_cooling_3e-8",
#     # "wind_-1e-3_cooling_3e-8",
#     # "wind_-1e-3_heating_-4e-8",
#     # "wind_-1e-3_heating_-1e-8",
#     # "wind_-1e-3_heating_-3e-8",
#     # "wind_-5e-4_heating_-5e-8",
#     # "wind_-5e-4_heating_-3e-8",
#     # "wind_-5e-4_heating_-1e-8",
#     # "wind_-2e-4_heating_-5e-8",
#     # "wind_-2e-4_heating_-3e-8",
#     # "wind_-2e-4_heating_-1e-8",
# ]

# animate_training_results(test_files, DATA_NAME,
#                          EXTRACTED_DATA_DIR=EXTRACTED_OUTPUT_PATH, OUTPUT_DIR=VIDEO_PATH)


text:
module WeakDivPerpTests
   using Gridap
   using GridapGeosciences
   using Plots
   using Test

   include("ConvergenceAnalysisTools.jl")

   # Nondivergent velocity field tangent to the sphere
   # Reference:
   #     Lauritzen et. al. GMD 2012
   #     Eq. (18) and (19)
   function W(θϕ)
     θ,ϕ = θϕ
     u = 2*sin(θ)^2*sin(2*ϕ) + 2π/5*cos(ϕ)
     v = 2*sin(2*θ)*cos(ϕ)
     spherical_to_cartesian_matrix(VectorValue(θ,ϕ,1.0))⋅VectorValue(u,v,0)
   end

   """
   Spherical divergence on unit sphere applied to perp(W)
   """
   function divM_perpW(θϕ)
     n_times_W=(θϕ)->(W(θϕ)×normal_unit_sphere(θϕ))
     f=divergence_unit_sphere(n_times_W)
     f(θϕ)
   end

   function compute_error_weak_div_perp(model,order,degree)
     # Setup geometry
     Ω=Triangulation(model)
     dΩ=Measure(Ω,degree)
     n=get_normal_vector(model)

     # Setup H(div) spaces
     reffe_rt = ReferenceFE(raviart_thomas, Float64, order)
     V = FESpace(model, reffe_rt ; conformity=:HDiv)
     U = TrialFESpace(V)

     # H1 spaces
     reffe_lgn = ReferenceFE(lagrangian, Float64, order+1)
     S = FESpace(model, reffe_lgn; conformity=:H1)
     R = TrialFESpace(S)

     # Project the analytical velocity field W onto the H(div) space
     a1(u,v) = ∫(v⋅u)dΩ
     b1(v)   = ∫(v⋅(W∘xyz2θϕ))dΩ
     op      = AffineFEOperator(a1,b1,U,V)
     wh      = solve(op)

     e = (W∘xyz2θϕ)-wh
     println(sqrt(sum(∫(e⋅e)dΩ)))

     # Compute weak div_perp operator
     a2(u,v) = ∫(v*u)dΩ
     b2(v)   = ∫(⟂(∇(v),n)⋅(-wh) )dΩ
     op      = AffineFEOperator(a2,b2,R,S)
     divwh   = solve(op)
     e = divwh-divM_perpW∘xyz2θϕ

     sqrt(sum(∫(e*e)dΩ))
   end
   @time ahs0,ak0errors,as0=convergence_study(compute_error_weak_div_perp,generate_n_values(2),0,4)
   @test round(as0,digits=3) ≈ 2.084

   @time bihs0,bik0errors,bis0=convergence_study(compute_error_weak_div_perp,generate_n_values(2),1,0,4)
   @test round(bis0,digits=3) ≈ 0.947

   @time biqhs0,biqk0errors,biqs0=convergence_study(compute_error_weak_div_perp,generate_n_values(2),2,0,4)
   @test round(biqs0,digits=3) ≈ 2.083

  #  plotd=plot([ahs0,bihs0,biqhs0],[ak0errors,bik0errors,biqk0errors],
  #  xaxis=:log, yaxis=:log,
  #  label=["k=0 analytical map" "k=0 bilinear map" "k=0 biquadratic map"],
  #  shape=:auto,
  #  xlabel="h",ylabel="L2 error norm", legend=:bottomright)

  #  savefig(plotd,"L2_error_weak_div_perp.png")

end


text:
# Script to test spglinecurvy

include("spglinecurvy.jl")

using GenSPGL
using MAT
var = matopen("/home/slim/klensink/.julia/v0.6/GenSPGL/compare/testspglc.mat")

gStep= read(var,"gStep")
tau = read(var,"tau")
lastFv = vec(read(var,"lastFv"))
b = vec(read(var,"b"))
x = vec(read(var,"x"))
g = vec(read(var,"g"))
timeProject = zero(Float64)
A = NLfunForward
opts = spgOptions(  optTol = 1e-5,
                    bpTol = 1e-5,
                    decTol = 1e-5,
                    project = TraceNorm_project,
                    primal_norm = TraceNorm_primal,
                    dual_norm = TraceNorm_dual,
                    proxy = true,
                    ignorePErr = true,
                    iterations = 150,
                    verbosity = 1)

# Avoid anon func
afunT(x) = reshape(x,355,709)
params = Dict{String, Any}([("nr", 40)
                                ("Ind", vec(b) .== 0)
                                ("numr", 355)
                                ("numc", 709)
                                ("funForward", NLfunForward)
                                ("afunT", afunT)
                                ("afun", afun)
                                ("mode", 1)
                                ("ls", 1)
                                ("logical", 0)
                                ("funPenalty", funLS)])
f, x, r, nLine, stepG, lnErr, localProdA = spglinecurvy(A, x, gStep*g, maximum(lastFv),
                                A, funLS, b, project, timeProject, tau, opts, params)


text:
using ArrayLayouts, Test
import ArrayLayouts: MemoryLayout, @_layoutlmul, triangulardata

include("test_layouts.jl")
include("test_muladd.jl")
include("test_ldiv.jl")

struct MyMatrix <: LayoutMatrix{Float64}
    A::Matrix{Float64}
end

Base.getindex(A::MyMatrix, k::Int, j::Int) = A.A[k,j]
Base.setindex!(A::MyMatrix, v, k::Int, j::Int) = setindex!(A.A, v, k, j)
Base.size(A::MyMatrix) = size(A.A)
Base.strides(A::MyMatrix) = strides(A.A)
Base.unsafe_convert(::Type{Ptr{T}}, A::MyMatrix) where T = Base.unsafe_convert(Ptr{T}, A.A)
MemoryLayout(::Type{MyMatrix}) = DenseColumnMajor()

@testset "LayoutMatrix" begin
    A = MyMatrix(randn(5,5))
    for (kr,jr) in ((1:2,2:3), (:,:), (:,1:2), (2:3,:), ([1,2],3:4), (:,[1,2]), ([2,3],:))
        @test A[kr,jr] == A.A[kr,jr]
    end
    b = randn(5)
    for Tri in (UpperTriangular, UnitUpperTriangular, LowerTriangular, UnitLowerTriangular)
        @test ldiv!(Tri(A), copy(b)) ≈ ldiv!(Tri(A.A), copy(b))
        @test lmul!(Tri(A), copy(b)) ≈ lmul!(Tri(A.A), copy(b))
    end

    @test copyto!(MyMatrix(Array{Float64}(undef,5,5)), A) == A
    @test copyto!(Array{Float64}(undef,5,5), A) == A
    @test copyto!(MyMatrix(Array{Float64}(undef,5,5)), A.A) == A
    @test copyto!(view(MyMatrix(Array{Float64}(undef,5,5)),1:3,1:3), view(A,1:3,1:3)) == A[1:3,1:3]
    @test copyto!(view(MyMatrix(Array{Float64}(undef,5,5)),:,:), A) == A
    @test copyto!(MyMatrix(Array{Float64}(undef,3,3)), view(A,1:3,1:3)) == A[1:3,1:3]
    @test copyto!(view(MyMatrix(Array{Float64}(undef,5,5)),:,:), A.A) == A
    @test copyto!(Array{Float64}(undef,3,3), view(A,1:3,1:3)) == A[1:3,1:3]

    @test qr(A).factors ≈ qr(A.A).factors
    @test qr(A,Val(true)).factors ≈ qr(A.A,Val(true)).factors
    @test lu(A).factors ≈ lu(A.A).factors
    @test lu(A,Val(true)).factors ≈ lu(A.A,Val(true)).factors
    @test_throws ErrorException qr!(A)
    @test_throws ErrorException lu!(A)
end

struct MyUpperTriangular{T} <: AbstractMatrix{T}
    A::UpperTriangular{T,Matrix{T}}
end

MyUpperTriangular{T}(::UndefInitializer, n::Int, m::Int) where T = MyUpperTriangular{T}(UpperTriangular(Array{T}(undef, n, m)))
MyUpperTriangular(A::AbstractMatrix{T}) where T = MyUpperTriangular{T}(UpperTriangular(Matrix{T}(A)))
Base.convert(::Type{MyUpperTriangular{T}}, A::MyUpperTriangular{T}) where T = A
Base.convert(::Type{MyUpperTriangular{T}}, A::MyUpperTriangular) where T = MyUpperTriangular(convert(AbstractArray{T}, A.A))
Base.convert(::Type{MyUpperTriangular}, A::MyUpperTriangular)= A
Base.convert(::Type{AbstractArray{T}}, A::MyUpperTriangular) where T = MyUpperTriangular(convert(AbstractArray{T}, A.A))
Base.convert(::Type{AbstractMatrix{T}}, A::MyUpperTriangular) where T = MyUpperTriangular(convert(AbstractArray{T}, A.A))
Base.convert(::Type{MyUpperTriangular{T}}, A::AbstractArray{T}) where T = MyUpperTriangular{T}(A)
Base.convert(::Type{MyUpperTriangular{T}}, A::AbstractArray) where T = MyUpperTriangular{T}(convert(AbstractArray{T}, A))
Base.convert(::Type{MyUpperTriangular}, A::AbstractArray{T}) where T = MyUpperTriangular{T}(A)
Base.getindex(A::MyUpperTriangular, kj...) = A.A[kj...]
Base.getindex(A::MyUpperTriangular, ::Colon, j::AbstractVector) = MyUpperTriangular(A.A[:,j])
Base.setindex!(A::MyUpperTriangular, v, kj...) = setindex!(A.A, v, kj...)
Base.size(A::MyUpperTriangular) = size(A.A)
Base.similar(::Type{MyUpperTriangular{T}}, m::Int, n::Int) where T = MyUpperTriangular{T}(undef, m, n)
Base.similar(::MyUpperTriangular{T}, m::Int, n::Int) where T = MyUpperTriangular{T}(undef, m, n)
Base.similar(::MyUpperTriangular, ::Type{T}, m::Int, n::Int) where T = MyUpperTriangular{T}(undef, m, n)
LinearAlgebra.factorize(A::MyUpperTriangular) = factorize(A.A)

MemoryLayout(::Type{MyUpperTriangular{T}}) where T = MemoryLayout(UpperTriangular{T,Matrix{T}})
triangulardata(A::MyUpperTriangular) = triangulardata(A.A)

@_layoutlmul MyUpperTriangular

@testset "MyUpperTriangular" begin
    A = randn(5,5)
    B = randn(5,5)
    x = randn(5)
    U = MyUpperTriangular(A)

    @test lmul!(U, copy(x)) ≈ U * x
    @test lmul!(U, copy(B)) ≈ U * B

    @test_skip lmul!(U,view(copy(B),collect(1:5),1:5)) ≈ U * B
end


text:
mutable struct ProjMPS
  lpos::Int
  rpos::Int
  nsite::Int
  M::MPS
  LR::Vector{ITensor}
  ProjMPS(M::MPS) = new(0, length(M) + 1, 2, M, Vector{ITensor}(undef, length(M)))
end

nsite(P::ProjMPS) = P.nsite

Base.length(P::ProjMPS) = length(P.M)

function lproj(P::ProjMPS)
  (P.lpos <= 0) && return nothing
  return P.LR[P.lpos]
end

function rproj(P::ProjMPS)
  (P.rpos >= length(P) + 1) && return nothing
  return P.LR[P.rpos]
end

function product(P::ProjMPS, v::ITensor)::ITensor
  if nsite(P) != 2
    error("Only two-site ProjMPS currently supported")
  end

  Lpm = dag(prime(P.M[P.lpos + 1], "Link"))
  !isnothing(lproj(P)) && (Lpm *= lproj(P))

  Rpm = dag(prime(P.M[P.rpos - 1], "Link"))
  !isnothing(rproj(P)) && (Rpm *= rproj(P))

  pm = Lpm * Rpm

  pv = scalar(pm * v)

  Mv = pv * dag(pm)

  return noprime(Mv)
end

#function Base.eltype(P::ProjMPS)
#  elT = eltype(P.M[P.lpos+1])
#  for j = P.lpos+2:P.rpos-1
#    elT = promote_type(elT,eltype(P.M[j]))
#  end
#  if !isnull(lproj(P))
#    elT = promote_type(elT,eltype(lproj(P)))
#  end
#  if !isnull(rproj(P))
#    elT = promote_type(elT,eltype(rproj(P)))
#  end
#  return elT
#end

(P::ProjMPS)(v::ITensor) = product(P, v)

#function Base.size(P::ProjMPS)::Tuple{Int,Int}
#  d = 1
#  if P.lpos > 0
#    d *= dim(linkind(M,P.lpos))
#  end
#  for j=P.lpos+1:P.rpos-1
#    d *= dim(siteind(P.M,j))
#  end
#  if P.rpos-1 < N
#    d *= dim(linkind(M,P.rpos-1))
#  end
#  return (d,d)
#end

function makeL!(P::ProjMPS, psi::MPS, k::Int)
  while P.lpos < k
    ll = P.lpos
    if ll <= 0
      P.LR[1] = psi[1] * dag(prime(P.M[1], "Link"))
      P.lpos = 1
    else
      P.LR[ll + 1] = P.LR[ll] * psi[ll + 1] * dag(prime(P.M[ll + 1], "Link"))
      P.lpos += 1
    end
  end
end

function makeR!(P::ProjMPS, psi::MPS, k::Int)
  N = length(P.M)
  while P.rpos > k
    rl = P.rpos
    if rl >= N + 1
      P.LR[N] = psi[N] * dag(prime(P.M[N], "Link"))
      P.rpos = N
    else
      P.LR[rl - 1] = P.LR[rl] * psi[rl - 1] * dag(prime(P.M[rl - 1], "Link"))
      P.rpos -= 1
    end
  end
end

function position!(P::ProjMPS, psi::MPS, pos::Int)
  makeL!(P, psi, pos - 1)
  makeR!(P, psi, pos + nsite(P))

  #These next two lines are needed
  #when moving lproj and rproj backward
  P.lpos = pos - 1
  return P.rpos = pos + nsite(P)
end


text:
using Pages, Mustache, PlotlyJS

Pages.start()

Endpoint("/examples/plot.ly") do request::Request
    d = Dict("host" => request.headers["Host"])
    template = Mustache.template_from_file(joinpath(dirname(@__FILE__),"PagesJL.html"))
    render(template,d)
end

function example_plotly()

n = 10
data = [scatter(;x = 1:n,y = rand(n))]
layout = Layout(;title = "Create a new <div> and plot a single trace",
    width = 600,
    height = 300
)
id = "Plot1"
Pages.add(Pages.Element(id=id))
Pages.broadcast("script","""Plotly.newPlot("$(id)",$(data),$(layout));""")

trace1 = scatter(;x = 1:n, y = rand(n))
trace2 = scatter(;x = 1:n, y = rand(n))
data = [trace1, trace2]
layout = Layout(;title = "Create a new <div> and plot an array of traces",
    width = 600,
    height = 300
)
id = "Plot2"
Pages.add(Pages.Element(id=id))
Pages.broadcast("script","""Plotly.newPlot("$(id)",$(data),$(layout));""")

trace1 = scatter(;
    x = [1, 2, 3, 4],
    y = [10, 15, 13, 17],
    mode = "markers")

trace2 = scatter(;
    x = [2, 3, 4, 5],
    y = [16, 5, 11, 10],
    mode = "lines")

trace3 = scatter(;
    x = [1, 2, 3, 4],
    y = [12, 9, 15, 12],
    mode = "lines+markers")

data = [trace1, trace2, trace3]

layout = Layout(;
    title = "Line and Scatter Plot",
    height = 400,
    width = 480)

id = "Plot3"
Pages.add(Pages.Element(id=id))
Pages.broadcast("script","""Plotly.newPlot("$(id)",$(data));""")

country = ["Switzerland (2011)", "Chile (2013)", "Japan (2014)",
           "United States (2012)", "Slovenia (2014)", "Canada (2011)",
           "Poland (2010)", "Estonia (2015)", "Luxembourg (2013)",
           "Portugal (2011)"]

votingPop = [40, 45.7, 52, 53.6, 54.1, 54.2, 54.5, 54.7, 55.1, 56.6]
regVoters = [49.1, 42, 52.7, 84.3, 51.7, 61.1, 55.3, 64.2, 91.1, 58.9]

# notice use of `attr` function to make nested attributes
trace1 = scatter(;x=votingPop, y=country, mode="markers",
                  name="Percent of estimated voting age population",
                  marker=attr(color="rgba(156, 165, 196, 0.95)",
                              line_color="rgba(156, 165, 196, 1.0)",
                              line_width=1, size=16, symbol="circle"))

trace2 = scatter(;x=regVoters, y=country, mode="markers",
                  name="Percent of estimated registered voters")
# also could have set the marker props above by using a dict
trace2["marker"] = Dict(:color => "rgba(204, 204, 204, 0.95)",
                       :line => Dict(:color=> "rgba(217, 217, 217, 1.0)",
                                     :width=> 1),
                       :symbol => "circle",
                       :size => 16)

data = [trace1, trace2]

layout = Layout(Dict{Symbol,Any}(:paper_bgcolor => "rgb(254, 247, 234)",
                                 :plot_bgcolor => "rgb(254, 247, 234)");
                title="Votes cast for ten lowest voting age population in OECD countries",
                width=600, height=600, hovermode="closest",
                margin=Dict(:l => 140, :r => 40, :b => 50, :t => 80),
                xaxis=attr(showgrid=false, showline=true,
                           linecolor="rgb(102, 102, 102)",
                           titlefont_font_color="rgb(204, 204, 204)",
                           tickfont_font_color="rgb(102, 102, 102)",
                           autotick=false, dtick=10, ticks="outside",
                           tickcolor="rgb(102, 102, 102)"),
                legend=attr(font_size=10, yanchor="middle",
                            xanchor="right"),
                )

id = "Plot4"
Pages.add(Pages.Element(id=id))
Pages.broadcast("script","""Plotly.newPlot("$(id)",$(data),$(layout));""")

end

# example_plotly()


text:
using C2Julia
using Clang
using Clang.LibClang

liblsoda_include = joinpath(@__DIR__, "liblsoda", "src")
liblsoda_c = joinpath(@__DIR__, "liblsoda", "src", "lsoda.c")

ctx = DefaultContext()

std_headers = find_std_headers()
compiler_flags = ["-I$liblsoda_include", "-I$LLVM_INCLUDE", ("-I".*std_headers)...]

push!(ctx.trans_units, TranslationUnit(ctx.index,
                                       liblsoda_c,
                                       compiler_flags,
                                       CXTranslationUnit_None))

root_cursor = getcursor(ctx.trans_units[1])
header = spelling(root_cursor)
ctx.children = children(root_cursor)

check_opt = search(root_cursor, "check_opt")[1]

outpath = joinpath(@__DIR__, "lsoda_test.jl")
open(outpath, "w+") do f
    println(f, translate(check_opt).expr)
end


text:
using SuperEnum

@se Colors begin
    RED = 0
    GREEN = 1
    BLUE = 2

end

@se Permissions begin
    R = 1
    W = 2
    X = 16

end

function show()
    color_map = Dict(Colors.RED => "red", Colors.GREEN => "green", Colors.BLUE => "blue")
    a = Colors.GREEN
    if a == Colors.GREEN
        println(join(["green"], " "))
    else

        println(join(["Not green"], " "))
    end
    b = Permissions.R
    if b == Permissions.R
        println(join(["R"], " "))
    else

        println(join(["Not R"], " "))
    end
    @assert(length(color_map) != 0)
end

function main()
    show()
end

main()


text:
import Base.Broadcast: broadcasted

let
    # Should this be Array{Tainted,1} ?
    @test jactestval((Y,X) -> typeof(X), [1], [2]) == Array{Int, 1}
    @test jactestval((Y,X) -> typeof(Y), [1], [2]) == Array{Int, 1}

    @test jactestval((Y,X) -> eltype(X), [1], [2]) == Int
    @test jactestval((Y,X) -> typeof(X[1]), [1], [2]) == Int
    @test jactestval((Y,X) -> typeof(X[1]/1), [1], [2]) == Float64

    f(y,x) = x .+ 1

    @test jactestval((Y,X) -> broadcasted(+, X, 1)[1], [1], [2]) == 3
    @test jactestval(f, [1], [2]) == [3]
    @test sparse(jactestmeta(f, [1], [2])) == sparse([], [], true, 1, 1)

    g(y,x) = y[:] .= x .+ 1
    #g(y,x) = y .= x .+ 1 -- memove

    println("Broadcast timings")
    println("  y .= x")
    # test path of unsafe_copy from Input to Output
    @test @time jacobian_sparsity((y,x) -> y .= x, [1,2,3], [1,2,3]) == sparse([1,2,3], [1,2,3], true)
    println("  y[:] .= x .+ 1")
    @test @time sparse(jactestmeta(g, [1], [2])) == sparse([1], [1], true)
    println("  y[1:2] .= x[2:3]")
    # test path of unsafe_copy from Input to an intermediary
    @test @time jacobian_sparsity((y,x) -> y[1:2] .= x[2:3], [1,2,3], [1,2,3]) == sparse([1,2],[2,3],true, 3,3)

    using LinearAlgebra, SparsityDetection

    function testsparse!(out, x)
        A = Tridiagonal(x[2:end], x, x[1:end-1])
        mul!(out, A, x)
    end
    x = [1:4;]; out = similar(x);
    @test jacobian_sparsity(testsparse!, out, x) == sparse([1,2,1,2,3,2,3,4,3,4],
                                                   [1,1,2,2,2,3,3,3,4,4], true)
end

@testset "BLAS" begin
    function f(out,in)
        A = rand(length(in), length(in))
        out .= A * in
        return nothing
    end

    x = [1.0:10;]
    out = similar(x)
    @test all(jacobian_sparsity(f, out, x) .== 1)
end

@testset "avoid branches in primitive functions with isleaf" begin
    # without the isleaf fix, this would go into an infinite loop
    # fixes issue #30
    x = rand(3)
    y = similar(x)
    function f(y, x)
        for i in 1:length(x)
            y[i] = exp(x[i])
        end
        return nothing
    end
    @test jacobian_sparsity(f, y, x) == sparse([1, 2, 3], [1, 2, 3], true)

    # this example tests that a function that gets tagged also indicates
    # isleaf correctly
    @test jacobian_sparsity(y,x) do y,x
        y .= exp.(x)
    end == sparse([1,2,3],[1,2,3],[1,1,1])
end

using SparsityDetection: Fixed
@testset "Fixed params" begin
    f(y,x,z) = z > 0 ? y .= x : y .= reverse(x)

    @test jacobian_sparsity(f, [1,1,1], [1,2,3], Fixed(1)) == sparse([1,2,3],[1,2,3],true)
    @test jacobian_sparsity(f, [1,1,1], [1,2,3], Fixed(-1)) == sparse([3,2,1],[1,2,3],true)
    @test jacobian_sparsity(f, [1,1,1], [1,2,3], 2) == sparse([1,2,3,3,1],[1,2,3,1,3],true)
end


text:
export subscript, superscript, randomcolor
export lighten_color, darken_color
export record_interaction

"""
    record_interaction(file, figure; framerate = 30, total_time = 10)
Start recording whatever interaction is happening on some `figure` into a video
output in `file` (recommended to end in `".mp4"`).

## Keywords
* `framerate = 30`
* `total_time = 10`: Time to record for, in seconds
# `sleep_time = 1`: Time to call `sleep()` before starting to save.
"""
function record_interaction(file, figure;
        framerate = 30, total_time = 10, sleep_time = 1,
    )
    ispath(dirname(file)) || mkpath(dirname(file))
    sleep(sleep_time)
    framen = framerate*total_time
    record(figure, file; framerate) do io
        for i = 1:framen
            sleep(1/framerate)
            recordframe!(io)
        end
    end
    return
end
record_interaction(figure::Figure, file; kwargs...) =
record_interaction(file, figure; kwargs...)

"""
    subscript(i::Int)
Transform `i` to a string that has `i` as a subscript.
"""
function subscript(i::Int)
    if i < 0
        "₋"*subscript(-i)
    elseif i == 1
        "₁"
    elseif i == 2
        "₂"
    elseif i == 3
        "₃"
    elseif i == 4
        "₄"
    elseif i == 5
        "₅"
    elseif i == 6
        "₆"
    elseif i == 7
        "₇"
    elseif i == 8
        "₈"
    elseif i == 9
        "₉"
    elseif i == 0
        "₀"
    else
        join(subscript.(digits(i)))
    end
end

"""
    superscript(i::Int)
Transform `i` to a string that has `i` as a superscript.
"""
function superscript(i::Int)
    if i < 0
        "⁻"*superscript(-i)
    elseif i == 1
        "¹"
    elseif i == 2
        "²"
    elseif i == 3
        "³"
    elseif i == 4
        "⁴"
    elseif i == 5
        "⁵"
    elseif i == 6
        "⁶"
    elseif i == 7
        "⁷"
    elseif i == 8
        "⁸"
    elseif i == 9
        "⁹"
    elseif i == 0
        "⁰"
    else
        join(superscript.(digits(i)))
    end
end

"""
    randomcolor(args...) = RGBAf0(0.9 .* (rand(), rand(), rand())..., 0.75)
"""
randomcolor(args...) = RGBAf0(0.9 .* (rand(), rand(), rand())..., 0.75)

function colors_from_map(cmap, α, N)
    N == 1 && return [Makie.to_colormap(cmap, 2)[1]]
    cs = [RGBAf0(c.r, c.g, c.b, α) for c in Makie.to_colormap(cmap, N)]
end

function pushupdate!(o::Observable, v)
    push!(o[], v)
    o[] = o[]
    return o
end

"""
    darken_color(c, f = 1.2)
Darken given color `c` by a factor `f`.
If `f` is less than 1, the color is lightened instead.
"""
function darken_color(c, f = 1.2)
    c = to_color(c)
    return RGBAf0(clamp.((c.r/f, c.g/f, c.b/f, c.alpha), 0, 1)...)
end

"""
    lighten_color(c, f = 1.2)
Lighten given color `c` by a factor `f`.
If `f` is less than 1, the color is darkened instead.
"""
function lighten_color(c, f = 1.2)
    c = to_color(c)
    return RGBAf0(clamp.((c.r*f, c.g*f, c.b*f, c.alpha), 0, 1)...)
end

"""
    to_alpha(c, α = 0.75)
Create a color same as `c` but with given alpha channel.
"""
function to_alpha(c, α = 1.2)
    c = to_color(c)
    return RGBAf0(c.r, c.g, c.b, α)
end

##########################################################################################
# Polygon stuff
##########################################################################################
using Makie.GeometryBasics # for using Polygons
export rotate2D, scale, Polygon, Point2f0

translate(p::Polygon, point) = Polygon(decompose(Point2f0, p.exterior) .+ point)

"""
    rotate2D(p::Polygon, θ)
Rotate given polygon counter-clockwise by `θ` (in radians).
"""
function rotate2D(p::Polygon, θ)
    sinφ, cosφ = sincos(θ)
    Polygon(map(
        p -> Point2f0(cosφ*p[1] - sinφ*p[2], sinφ*p[1] + cosφ*p[2]),
        decompose(Point2f0, p.exterior)
    ))
end

"""
    scale(p::Polygon, s)
Scale given polygon by `s`, assuming polygon's "center" is the origin.
"""
scale(p::Polygon, s) = Polygon(decompose(Point2f0, p.exterior) .* Float32(s))


text:
# test_meshwarp()
println("testing meshwarp")
img = reshape(float(collect(1:121) .% 2), 11, 11) # 11x11 checkerboard
triangles = [1 2 3]

tform = [
    1 0 0;
    0 1 0;
    0 0 1
]
src = [
    0.0 0.0;
    0.0 10.0;
    10.0 0.0
]
dst = warp_pts(tform, src)
offset = [0, 0]
img_warped, warped_offset = meshwarp(img, Array(src'), Array(dst'), Array(triangles'), offset)
# @test_approx_eq_eps img_warped img 1e-5
@test warped_offset == [0, 0]

tform = [
    1 0 0;
    0 1 0;
    0 0 1]
src = [
    0.0 0.0;
    0.0 10.0;
    10.0 0.0]
dst = warp_pts(tform, src)
offset = [5, 10]
img_warped, warped_offset = meshwarp(img, Array(src'), Array(dst'), Array(triangles'), offset)
@test img_warped ≈ zeros(11, 11) atol=1e-10
@test warped_offset == [0, 0]

tform = [1 0 0; 0 1 0; 0 0 1]
src = [0.0 0.0; 0.0 20.0; 20.0 0.0]
dst = warp_pts(tform, src)
offset = [5, 10]
img_warped, warped_offset = meshwarp(img, Array(src'), Array(dst'), Array(triangles'), offset)
@test img_warped[6, 11] == 1.0
@test warped_offset == [0, 0]

src = [2.0 2.0; 10.0 2.0; 6.0 10.0]
dst = [4.0 2.0; 8.0 2.0; 6.0 10.0]
offset = [0, 0]
img_warped, warped_offset = meshwarp(img, Array(src'), Array(dst'), Array(triangles'), offset)
@test warped_offset == [4, 2]

# test_poly2mask()
# vertices=[1.5 2.5; 2.5 1.5; 3.5 4.5]
# img, off = poly2mask(vertices[:,1],vertices[:,2])
# """ off should be [1,1] and img should be:
# false  false  false  false  false
# false   true   true  false  false
# false  false   true   true  false
# false  false  false  false  false
# """
# vertices=[1.5 2.4; 2.5 1.5; 3.5 4.5]
# img, off = poly2mask(vertices[:,1],vertices[:,2])
# """ off should be [1,1] and img should be:
# false  false  false  false  false
# false   true  false  false  false
# false  false   true  false  false
# false  false  false  false  false,
# """
# vertices=[1 2; 2.5 1.5; 3.5 4.5]
# img, off = poly2mask(vertices[:,1],vertices[:,2])
# """ corner case: off should be [1,1] and img should be:
# false   true  false  false  false
# false   true   true  false  false
# false  false   true   true  false
# false  false  false  false  false
# """
# p=plot(x=vertices[:,1],y=vertices[:,2])
# display(p)
# img,off


text:
using Statistics
function fit(x; axis = 1) #this function dose only performs calculations.
    if axis == 2
        x = x'
    end
    if ndims(x) == 1
        p = [mean(x), std(x)]
    else p = []
        for i in 1:size(x)[2]
            push!(p, [mean(x[:, i]), std(x[:, i])])
        end
    end
    return p
end
function transform(x, p; axis = 1) #Only transform is performed based on the calculation results of the'fit'function and'fit_transform'function.
    if axis == 2
        x = x'
    end
    if ndims(x) == 1
        t = (x .- p[1]) / p[2]
    else
        t = []
        for i in 1:size(x)[2]
            push!(t, (x[:, i] .- p[i][1]) / p[i][2])
        end
    end
    return t
end
function inverse_transform(x, p; axis = 1)
    if axis == 2
        x = x'
    end
    if ndims(x) == 1
        t = (x * p[2]) .+ p[1]
    else
        t = []
        for i in 1:size(x)[2]
            push!(t, (x[:, i][1] * p[i][2]) .+ p[i][1])
        end
    end
    return t
end
function fit_transform(x; axis = 1) #This function handles both the fit function and the transform function.
    if axis == 2
        x = x'
    end
    if ndims(x) == 1
        p = [mean(x), std(x)]
        t = (x .- p[1]) / p[2]
    else t, p = [], []
        for i in 1:size(x)[2]
            push!(p, [mean(x[:, i]), std(x[:, i])])
            push!(t, (x[:, i] .- p[i][1]) / p[i][2])
        end
    end
    return t, p
end


text:
#####
##### There's no free surface with a rigid lid
#####

explicit_barotropic_pressure_x_gradient(i, j, k, grid, free_surface::Nothing) = zero(eltype(grid))
explicit_barotropic_pressure_y_gradient(i, j, k, grid, free_surface::Nothing) = zero(eltype(grid))


text:
@testset "Wigner3j" begin
    w3j_prealloc = zeros(Cdouble, 5)
    for j2 in 0:2, j3 = 0:2, m2 in -j2:j2, m3 in -j3:j3
        jmax = j2 + j3
        for m1 in -jmax:jmax
            w3j, jmin, jmax = Wigner3j(j2, j3, m1, m2, m3)
            Wigner3j!(w3j_prealloc, j2, j3, m1, m2, m3)
            @test w3j_prealloc[1:(jmax - jmin + 1)] == w3j
        end
    end

    @testset "known values" begin
        for j in 0:20, m in -j:j
            w3j, jmin, jmax = Wigner3j(j, 0, m, -m, 0)
            @test length(w3j) == 1
            @test w3j[1] ≈ (-1)^(j-m)/√(2j+1)
        end
    end
end


text:
using DemoPackageZW
using Test

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


text:
mdpath, resultpath = ARGS
using PkgBenchmark: PkgBenchmark, export_markdown
result = PkgBenchmark.readresults(resultpath)

if mdpath == "-"
    export_markdown(stdout, result)
else
    open(mdpath, write = true) do io
        export_markdown(io, result)
    end
end


text:
using ComplexMixtures
using PDBTools
using Plots
using LaTeXStrings
using EasyFit

function fig() # to simplify globals

# Plot defaults
plot_font = "Computer Modern"
default(
    fontfamily=plot_font,
    linewidth=2.5,
    framestyle=:box,
    label=nothing,
    grid=false,
    palette=:tab10
)
scalefontsizes()
#scalefontsizes(1.3)

# Load system PDB file
system = readPDB("../simulation/equilibrated.pdb")

# Select the atoms corresponding to glycerol and water
glyc = select(system,"resname GLLM")
water = select(system,"water")

# Set the solute structure for Glycerol, used to extract atomic contributions
solute = Selection(glyc,natomspermol=14)

# Load previously computed data
mddf_glyc = load("./mddf_glyc.json")
mddf_water_glyc = load("./mddf_water_glyc.json")

l = @layout [
    grid(2,1)
    grid(2,1)
    grid(2,1)
]
p = plot(layout=l)

# Plot the correlation functions
x = mddf_glyc.d # distances
y = movavg(mddf_glyc.mddf,n=10).x # the mddf (using movavg to smooth noise)
plot!(p,x,y,label="Glycerol-Glycerol",subplot=1)
x = mddf_water_glyc.d
y = movavg(mddf_water_glyc.mddf,n=10).x
plot!(p,x,y,label="Water-Glycerol",subplot=1)
plot!(p,
    ylabel="MDDF",
    xlim=(0,20),subplot=1,
    xticks=:none
)

# Plot the KB integrals
y = movavg(mddf_glyc.kb,n=10).x
plot!(p,x,y,subplot=2)
y = movavg(mddf_water_glyc.kb,n=10).x
plot!(p,x,y,subplot=2)
plot!(p,
    xlabel=L"\textrm{Distance / \AA}",ylabel=L"\textrm{KB~/~cm^3~mol^{-1}}",
    xlim=(0,20),subplot=2
)

# Plot some group contributions to the MDDF. We select the atom names
# corresponding to each type of group of the glycerol molecule.
hydroxyls = ["O1","O2","O3","HO1","HO2","HO3"]
aliphatic = ["C1","C2","C3","H11","H12","H2","H31","H32"]

# Retrieve the contributions of these groups to the MDDFs
hydroxyl_contrib = contrib(solute,mddf_glyc.solvent_atom,hydroxyls)
aliphatic_contrib = contrib(solute,mddf_glyc.solvent_atom,aliphatic)

# Plot group contributions
x = mddf_glyc.d
y = movavg(mddf_glyc.mddf,n=10).x
plot!(p,x,y,label="Total",subplot=3)
y = movavg(hydroxyl_contrib,n=10).x
plot!(p,x,y,label="Hydroxyl contributions",subplot=3)
y = movavg(aliphatic_contrib,n=10).x
plot!(p,x,y,label="Aliphatic contributions",subplot=3)
plot!(p,
    xlim=(1,6),
    ylim=(0,3.8),
    ylabel="MDDF",
    subplot=3,
    xticks=:none
)

x = mddf_water_glyc.d
y = movavg(mddf_water_glyc.mddf,n=10).x
plot!(p,x,y,subplot=4)
hydroxyl_contrib = contrib(solute,mddf_water_glyc.solute_atom,hydroxyls)
aliphatic_contrib = contrib(solute,mddf_water_glyc.solute_atom,aliphatic)

y = movavg(hydroxyl_contrib,n=10).x
plot!(p,x,y,subplot=4)
y = movavg(aliphatic_contrib,n=10).x
plot!(p,x,y,subplot=4)
plot!(p,
    xlim=(1,6),
    ylim=(0,3.8),
    xlabel=L"\textrm{Distance / \AA}",
    ylabel="MDDF",
    subplot=4
)

# 2D maps plot of group contributions

# Glycerol groups
groups = [
    (["O1","HO1"],"OH"), # first hydroxyl
    (["C1","H11","H12"],L"\textrm{CH_2}"), # first CH2
    (["O2","HO2"],"OH"), # second hydroxyl
    (["C2","H2"],"CH"), # CH
    (["C3","H31","H32"],L"\textrm{CH_2}"), # second CH2
    (["O3","HO3"],"OH") # third hydroxyl
]

# Contributions to Glycerol-Glycerol autocorrelation
group_contrib = zeros(length(mddf_glyc.d),length(groups))
for (igroup, group) in pairs(groups)
    group_contrib[:,igroup] .= contrib(solute,mddf_glyc.solute_atom,group[1])
end

idmin = findfirst( d -> d > 1.5, mddf_glyc.d)
idmax = findfirst( d -> d > 3.0, mddf_glyc.d)
labels = [ "OH", L"\textrm{CH_2}", "OH", "CH", L"\textrm{CH_2}", "OH" ]

contourf!(p,
    1:length(groups),
    mddf_glyc.d[idmin:idmax],
    group_contrib[idmin:idmax,:],
    color=cgrad(:tempo),linewidth=1,linecolor=:black,
    colorbar=:none,levels=5,
    xticks=:none,
    ylabel=L"\textrm{Distance/\AA}",
    subplot=5
)

# Water-glycerol interactions (Glycerol contributions)
group_contrib = zeros(length(mddf_glyc.d),length(groups))
for (igroup, group) in pairs(groups)
    group_contrib[:,igroup] .= contrib(solute,mddf_water_glyc.solute_atom,group[1])
end
contourf!(p,
    1:length(groups),
    mddf_glyc.d[idmin:idmax],
    group_contrib[idmin:idmax,:],
    color=cgrad(:tempo),linewidth=1,linecolor=:black,
    colorbar=:none,levels=5,
    xticks=(1:length(groups),labels),xrotation=60,
    ylabel=L"\textrm{Distance/\AA}",
    subplot=6
)
plot!(p,
    xlabel="Glycerol group",
    subplot=6
)

for (y,lab) in [(13.5,"A)"),
                (11.7,"B)"),
                (8.8,"C)"),
                (7.0,"D)"),
                (4.10,"E)"),
                (2.20,"F)")]
    annotate!(p,
        0, y,
        text(lab,plot_font,10),
        subplot=6,
    )
end

annotate!(p,1.1, 3.4, text("Glycerol-Glycerol", plot_font, 10, :left), subplot=3)
annotate!(p,1.1, 3.4, text("Water-Glycerol", plot_font, 10, :left), subplot=4)

annotate!(p, 4, 1.7, text("Glycerol-Glycerol", plot_font, 10, :left), subplot=5)
annotate!(p, 4, 1.7, text("Water-Glycerol", plot_font, 10, :left), subplot=6)

plot!(p,
    size=(550,800),
    leftmargin=10Plots.Measures.mm # adjust margin
)

savefig(p,"./glyc_water.png")
savefig(p,"./figure3.pdf")

return p
end;
fig()


text:
# Simple Advection problem

# build mesh
include("../src/DiscontinuousGalerkin.jl")
using DiscontinuousGalerkin

const CFL = 0.5
const Tend = 1.0

f(::Type{Val{:jac}},u::Vector) = eye(size(u,1))
f(u::Vector) = u

#define max wave speed
max_w_speed(u, f) = 1

f0(x) = sin(2*π*x)

function get_problem(N)
  mesh = Uniform1DFVMesh(N,-2.0,2.0,:PERIODIC,:PERIODIC)
  ConservationLawsProblem(f0,f,CFL,Tend,mesh)
end
#Compile
prob = get_problem(10)

#Run
prob = get_problem(200)

#Generate basis for Discontinuous Galerkin Scheme
basis=legendre_basis(3)

#Optional dgLimiter
limiter! = DGLimiter(prob, basis, Linear_MUSCL_Limiter())

#Solve problem
@time sol = solve(prob, DiscontinuousGalerkinScheme(basis, advection_num_flux; max_w_speed = max_w_speed); TimeIntegrator = SSPRK22(limiter!))

#Plot
using Plots
plot(sol)

max_stars_repo_name: Paulms/DiscontinuousGalerkin

max_issues_repo_name: Paulms/DiscontinuousGalerkin

max_forks_repo_name: Paulms/DiscontinuousGalerkin


text:
module Cases
using DataFrames, CSV

include("../src/DistSystem.jl")
import .DistSystem

export case3_dist, case33, case3_dist_no_dgs

function case3_dist()
    VN = 12500
    Sb = 1e6
    Yb = Sb / (VN^2)
    Vsub = 1.0
    P_limt = 10
    Q_limit = 10
    bus_sub = 1
    nbuses = 3
    VL = 0.90
    VH = 1.05
    sub = DistSystem.Substation(VN, bus_sub, Vsub, P_limt, Q_limit, [0.003, 12, 240])
    data = DataFrame(CSV.File("case3_dist.csv"))

    sys = DistSystem.factory_system(data, VL, VH, sub)

    sys.dgs = [
        DistSystem.DER(3, 0.001, 0.001, 0.05, 1.0, [0.0, 0.03], [-0.03, 0.03], [0.0, 1.0], [0.0026, 10.26, 210]), #DG dispatchable
        DistSystem.DER(2, 0.002, 0.002, 0.05, 1.0, [0.0, 0.03], [-0.03, 0.03], [-1.0, 0.0, 1.0], [0.0026, 10.26, 210]), #DG dispatchable
    ]
    sys.m_load = [0.5, 0.8, 1.0]
    sys.m_new_dg = [-1.0, 0.0, 1]

    return sys, "case3_dist"
end

function case3_dist_no_dgs()
    VN = 12500
    Sb = 1e6
    Yb = Sb / (VN^2)
    Vsub = 1.0
    P_limt = 10
    Q_limit = 10
    bus_sub = 1
    nbuses = 3
    VL = 0.90
    VH = 1.05
    sub = DistSystem.Substation(VN, bus_sub, Vsub, P_limt, Q_limit, [0.003, 12, 240])
    data = DataFrame(CSV.File("case3_dist.csv"))

    sys = DistSystem.factory_system(data, VL, VH, sub)
    sys.m_load = [0.5, 0.8, 1.0]
    sys.m_new_dg = [-1.0, 0.0, 1]

    return sys, "case3_dist_no_dgs"
end

function case33()
    VN = 12660
    Sb = 1e6
    Yb = Sb / (VN^2)
    Vsub = 1.0
    P_limt = 4
    Q_limit = 2.5
    bus_sub = 1
    nbuses = 33
    VL = 0.90
    VH = 1.05
    sub = DistSystem.Substation(VN, bus_sub, Vsub, P_limt, Q_limit, [0.003, 12, 240])
    data = DataFrame(CSV.File("case33.csv"))
    sys = DistSystem.factory_system(data, VL, VH, sub)
    sys.dgs = [
        DistSystem.DER(18, 0.02, 0.02, 0.05, 0.05, [0.0, 0.03], [-0.03, 0.03], [0.0, 1.0], [0.0026, 10.26, 210]), #DG dispatchable
        DistSystem.DER(22, 0.02, 0.02, 0.05, 0.05, [0.0, 0.0], [-0.03, 0.03], [0.0, 1.0], [0.0026, 10.26, 210]), #DG non-active dispatchable
        DistSystem.DER(33, 0.02, 0.02, 0.15, 0.05, [-0.01, 0.01], [-0.01, 0.01], [-1.0, 0.0, 1.0], [0.0026, 10.26, 210]), #ESS
        DistSystem.DER(25, 0.02, 0.02, 0.0, 0.05, [0.0, 0.0], [0.0, 0.0], [0.0, 1.0], [0.0026, 10.26, 210]), #DG non-dispatchable
    ]
    sys.m_load = [0.5, 0.8, 1.0]
    sys.m_new_dg = [-1.0, 0.0, 1]
    return sys, "case33"
end

end


text:
include("calcPSF.jl")

##################
print("Program begins...... ")
tic()
println("Time in...")

###

max_stars_repo_name: xwghua/DRiP-public

max_issues_repo_name: xwghua/DRiP-public

max_forks_repo_name: xwghua/DRiP-public


text:
using AlgebraicDiffEq, DiffEqDevTools, DiffEqProblemLibrary
using Base.Test


text:
module AttachmapTests

using Test
using InplaceArrays.TensorValues
using InplaceArrays.Arrays
using InplaceArrays.Fields
using InplaceArrays.Fields: OtherMockBasis, MockBasis
using FillArrays

p1 = Point(2,2)
p2 = Point(4,2)
p3 = Point(1,3)
p4 = Point(5,2)
x = [p1,p2,p3,p4]

np = length(x)
v = 2.0
d = 2
ndof = 8
wi = 3.0
w = fill(wi,ndof)
f = OtherMockBasis{d}(ndof)
r = MockBasis{d}(v,ndof)
ϕ = lincomb(f,w)

b = attachmap(r,ϕ)
rx = evaluate(r,x)
bx = rx
∇rx = evaluate(∇(r),x)
∇ϕx = evaluate(∇(ϕ),x)
∇bx = zeros(VectorValue{d,Float64},np,ndof)
for i in 1:np
  jacinv = inv(∇ϕx[i])
  for j in 1:ndof
    ∇bx[i,j] = jacinv*∇rx[i,j]
  end
end
test_field(b,x,bx,grad=∇bx)

l = 10
ax = fill(x,l)
af = Fill(f,l)
ar = Fill(r,l)
aw = fill(w,l)
aϕ = lincomb(af,aw)

ab = attachmap(ar,aϕ)
abx = fill(bx,l)
a∇bx = fill(∇bx,l)
test_array_of_fields(ab,ax,abx,grad=a∇bx)

end # module


text:
#--------------------------------------
# This script defines a bunch of
# helper functions to assist with the
# core functionality of PosteriorPlots.jl
#--------------------------------------

#--------------------------------------
# Author: Trent Henderson, 23 July 2021
#--------------------------------------

#----------- Histogram ------------

function histhelper(data::DataFrame, centrality::String, p::Symbol, add_legend::Bool)

    # Compute median

    if centrality == "median"
        m = median(data[!, p])
    else
        m = mean(data[!, p])
    end

    # Set up graphics helpers

    mycolor = cgrad(:blues)[1]
    mycolor2 = cgrad(:blues)[2]

    # Draw plot

    myPlot = plot(data[!, p], seriestype = :histogram, fillalpha = 0.5,
                  xlabel = "Value", ylabel = "", label = "",
                  color = mycolor, title = string(p), size = (800, 800),
                  legend = add_legend)

    plot!([m], seriestype = "vline", color = mycolor2, label = centrality, linewidth = 2.5)

    return myPlot
end

#----------- Density --------------

function denshelper(data::DataFrame, centrality::String, p::Symbol, lower::Float64, upper::Float64, add_legend::Bool)

    Random.seed!(123)

    # Compute median

    thevec = convert(Vector, data[!, p])

    if centrality == "median"
        m = median(thevec)
    else
        m = mean(thevec)
    end

    y = kde(thevec)
    lowerquantile = quantile(thevec, lower)
    upperquantile = quantile(thevec, upper)

    # Set up graphics helpers

    mycolor = theme_palette(:auto).colors.colors[1]

    # Draw plot

    myPlot = plot(thevec, title = string(p), fillalpha = 0.8,
         xlabel = "Value", ylabel = "Density", label = "",
         seriestype = :density, color = mycolor, size = (800, 800))

    plot!(range(lowerquantile, stop = upperquantile, length = 100), thevec -> pdf(y,thevec),
           color =  mycolor, fill = (0, 0.4, mycolor), label = "Credible Interval", legend = add_legend)

    plot!([m], seriestype = "vline", color = mycolor, label = centrality, linewidth = 2.5)

    return myPlot
end

#----------- Soss.jl parser -------

# Helper function to split an arbitrarily sized array returned by `dynamicHMC`
# into separate columns for each parameter. This likely applies to coefficients more often than
# not and potentially intercepts if mixed effects models are used.

function expandcol(df::DataFrame, thecol)

    expandlength =  length(df[1, thecol])
    @assert all(length.(df[!, thecol]) .== expandlength)
    expandnames = [Symbol(thecol, "_", i) for i ∈ 1:expandlength]
    subdf = DataFrame([getindex.(df[:, thecol], i) for i ∈ 1:lastindex(df[1, thecol])], expandnames)
    return hcat(df, subdf)
end

#----------- Quantile maker -------

# Helper function to produce 2 numbers that bound a quantile for a user-specified probability

function generatequantile(x::Float64)

    upperbound = round(((x + 1) / 2), digits = 3) # Solve for upper bound
    lowerbound = round(upperbound - x, digits = 3) # Substitute in for lower bound
    return lowerbound, upperbound
end


text:
using RCall, JLD, DataFrames, DataFramesMeta, GibbsTypePriors, Nemo

accuracy_Cnk_sigma = load("test/graphical_tests/saves_for_graphical_tests/accuracy_Cnk_sigma.jld", "accuracy_Cnk_sigma") |>
    df -> @transform(df, Accuracy = min.(:Accuracy, GibbsTypePriors.RR(22)/GibbsTypePriors.RR(7.) |> accuracy_bits))

accuracy_Cnk_rec_sigma = load("test/graphical_tests/saves_for_graphical_tests/accuracy_Cnk_rec_sigma.jld", "accuracy_Cnk_sigma") |>
        df -> @transform(df, Accuracy = min.(:Accuracy, GibbsTypePriors.RR(22)/GibbsTypePriors.RR(7.) |> accuracy_bits))  |>
        df -> @transform(df, s = Float64.(:s))

R"library(tidyverse)"

R"p = as_tibble($accuracy_Cnk_sigma) %>%
    mutate(type = ifelse(test = n == k, yes = 'Large k', no = 'Small k'), fun = 'Iterative') %>%
    bind_rows(as_tibble($accuracy_Cnk_rec_sigma) %>%
        mutate(type = ifelse(test = n == k, yes = 'Large k', no = 'Small k'), fun = 'Recursive')) %>%
    ggplot(aes(x = n, y = Accuracy, colour = s, group = interaction(s, type))) +
    theme_bw() +
    facet_wrap(~fun) +
    geom_point(aes(shape = type)) +
    geom_line(aes(linetype = type)) +
    scale_linetype(name='') +
    scale_shape(name='') +
    viridis::scale_colour_viridis(name=expression(sigma)) +
    geom_hline(yintercept = 64, linetype = 'dotted', colour = 'red')
    "

R"pdf('test/graphical_tests/figures_graphical_tests/accuracy_Cnk_sigma.pdf', height = 5, width = 5)
    plot(p)
dev.off()"

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
using RDatasets, DataFrames, FixedEffectModels, GLM, Base.Test

x = dataset("plm", "Cigar")
x[:pState] = pool(x[:State])
x[:pYear] = pool(x[:Year])

function glm_helper(formula::Formula, x::DataFrame)
    model_response(ModelFrame(formula, x)) - predict(glm(formula, x, Normal(), IdentityLink()))
end
function glm_helper(formula::Formula, x::DataFrame, wts::Symbol)
    model_response(ModelFrame(formula, x)) - predict(glm(formula, x, Normal(), IdentityLink(), wts = convert(Array{Float64}, x[wts])))
end

test = (
    convert(Array{Float64}, partial_out(Sales + Price ~ NDI, x)),
    convert(Array{Float64}, partial_out(Sales + Price ~ NDI |> pState, x)),
    convert(Array{Float64}, partial_out(Sales + Price ~ 1 |> pState, x)),
    convert(Array{Float64}, partial_out(Sales + Price ~ 1, x)),
    mean(convert(Array{Float64}, partial_out(Sales + Price ~ NDI, x, add_mean = true)), 1),
    mean(convert(Array{Float64}, partial_out(Sales + Price ~ NDI |> pState, x, add_mean = true)), 1),
    mean(convert(Array{Float64}, partial_out(Sales + Price ~ 1 |> pState, x, add_mean = true)), 1),
    mean(convert(Array{Float64}, partial_out(Sales + Price ~ 1, x, add_mean = true)), 1),
    convert(Array{Float64}, partial_out(Sales + Price ~ NDI, x, weight = :Pop)),
    convert(Array{Float64}, partial_out(Sales + Price ~ NDI |> pState, x, weight = :Pop)),
    convert(Array{Float64}, partial_out(Sales + Price ~ 1 |> pState, x, weight = :Pop)),
    convert(Array{Float64}, partial_out(Sales + Price ~ 1, x, weight = :Pop)),
    )

answer = (
    hcat(glm_helper(Sales ~ NDI, x), glm_helper(Price ~ NDI, x)),
    hcat(glm_helper(Sales ~ NDI + pState, x), glm_helper(Price ~ NDI + pState, x)),
    hcat(glm_helper(Sales ~ pState, x), glm_helper(Price ~ pState, x)),
    hcat(glm_helper(Sales ~ 1, x), glm_helper(Price ~ 1, x)),
    hcat(mean(x[:Sales]), mean(x[:Price])),
    hcat(mean(x[:Sales]), mean(x[:Price])),
    hcat(mean(x[:Sales]), mean(x[:Price])),
    hcat(mean(x[:Sales]), mean(x[:Price])),
    hcat(glm_helper(Sales ~ NDI, x, :Pop), glm_helper(Price ~ NDI, x, :Pop)),
    hcat(glm_helper(Sales ~ NDI + pState, x, :Pop), glm_helper(Price ~ NDI + pState, x, :Pop)),
    hcat(glm_helper(Sales ~ pState, x, :Pop), glm_helper(Price ~ pState, x, :Pop)),
    hcat(glm_helper(Sales ~ 1, x, :Pop), glm_helper(Price ~ 1, x, :Pop))
    )

for i in 1:12
    @test_approx_eq_eps test[i] answer[i]	1e-5
end

x[1, :Sales] = NA
x[2, :Price]  = NA
x[5, :Pop]  = NA
x[6, :Pop]  = -1.0
partial_out(Sales + Price ~ 1, x, weight = :Pop)


text: const DB = "db.sqlite"

max_stars_repo_name: Arkoniak/ZulipSnippetBot

max_issues_repo_name: Arkoniak/ZulipSnippetBot

max_forks_repo_name: Arkoniak/ZulipSnippetBot


text:
export GeneralizedMinimalResidualMethod

"""
    GeneralizedMinimalResidualMethod(; M=30, K=10)

# GMRES
This object represents an iterative Krylov method for solving a linear system.
The constructor parameter `M` is the number of steps after which the algorithm
is restarted (if it has not converged), `K` is the maximal number of restarts.
The amount of memory required for the solver state is roughly `(M + 1) * N`, where `N`
is the number of unknowns. This uses the restarted Generalized Minimal
Residual method of Saad and Schultz (1986).

## References
    @article{saad1986gmres,
      title={GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems},
      author={Saad, Youcef and Schultz, Martin H},
      journal={SIAM Journal on scientific and statistical computing},
      volume={7},
      number={3},
      pages={856--869},
      year={1986},
      publisher={SIAM}
    }
"""
struct GeneralizedMinimalResidualMethod <: AbstractKrylovMethod
    "Maximum number of Krylov iterations"
    M::Int
    "Maximum number of restarts"
    K::Int
    function GeneralizedMinimalResidualMethod(; M=30, K=10)
        return new(M, K)
    end
end

mutable struct GMRESCache{M, MP1, MMP1, T, AT} <: AbstractLinearSolverCache
    krylov_basis::NTuple{MP1, AT}
    "Hessenberg matrix"
    H::Matrix{T}
    "rhs of the least squares problem"
    g0::Vector{T}
    "work vector for preconditioning"
    Wvec::AT
end

function cache(
    krylov_alg::GeneralizedMinimalResidualMethod,
    Q::AT,
) where {AT}
    Wvec = similar(Q)
    M = krylov_alg.M
    krylov_basis = ntuple(i -> similar(Q), M + 1)
    FT = eltype(Q)
    H = zeros(M + 1, M)
    g0 = zeros(M + 1)

    return GMRESCache{M, M + 1, M * (M + 1), eltype(Q), AT}(
        krylov_basis,
        H,
        g0,
        Wvec,
    )
end

function LSinitialize!(
    ::GeneralizedMinimalResidualMethod,
    solver::LinearSolver,
    Q,
    Qrhs,
    args...,
)
    linearoperator! = solver.linop!
    pc = solver.pc

    cache = solver.cache
    g0 = cache.g0
    krylov_basis = cache.krylov_basis
    atol = solver.atol

    @assert size(Q) == size(krylov_basis[1])

    # store the initial residual in krylov_basis[1]
    linearoperator!(krylov_basis[1], Q, args...)
    @. krylov_basis[1] = Qrhs - krylov_basis[1]

    # apply left preconditioning
    if isa(pc.pc_side, PCleft)
        PCapply!(pc, cache.Wvec, krylov_basis[1], args...)
        @. krylov_basis[1] = cache.Wvec
    end

    residual_norm = norm(krylov_basis[1])
    threshold = solver.rtol * residual_norm

    converged = false
    if threshold < solver.atol
        converged = true
        return converged, threshold
    end

    fill!(g0, 0)
    g0[1] = residual_norm
    krylov_basis[1] ./= residual_norm

    converged, max(threshold, atol)
end

function LSsolve!(
    krylov_alg::GeneralizedMinimalResidualMethod,
    solver::LinearSolver,
    threshold,
    Q,
    Qrhs,
    args...,
)
    converged = false
    iter = 0
    total_iters = 0
    residual_norm = typemax(eltype(Q))

    while !converged && iter < krylov_alg.K
        converged, cycle_iters, residual_norm = gmres_cycle!(solver, threshold, Q, Qrhs, args...)

        iter += 1
        total_iters += cycle_iters

        # If we blow up, we want to know about it
        if !isfinite(residual_norm)
            error("norm of residual is not finite after $total_iters iterations.")
        end
    end

    (converged, total_iters, residual_norm)
end

function gmres_cycle!(
    solver::LinearSolver,
    threshold,
    Q,
    Qrhs,
    args...,
)
    linearoperator! = solver.linop!
    pc = solver.pc
    cache = solver.cache
    krylov_basis = cache.krylov_basis
    H = cache.H
    g0 = cache.g0
    rtol = solver.rtol
    atol = solver.atol

    converged = false
    residual_norm = typemax(eltype(Q))
    Ω = LinearAlgebra.Rotation{eltype(Q)}([])
    j = 1
    for outer j in 1:solver.krylov_alg.M
        # apply right preconditioning
        if isa(pc.pc_side, PCright)
            PCapply!(pc, cache.Wvec, krylov_basis[j], args...)
            @. krylov_basis[j] = cache.Wvec
        end

        # apply the linear operator
        linearoperator!(krylov_basis[j + 1], krylov_basis[j], args...)

        # apply left preconditioning
        if isa(pc.pc_side, PCleft)
            PCapply!(pc, cache.Wvec, krylov_basis[j + 1], args...)
            @. krylov_basis[j + 1] = cache.Wvec
        end

        # Arnoldi using the Modified Gram Schmidt orthonormalization
        for i in 1:j
            H[i, j] = dot(krylov_basis[j + 1], krylov_basis[i])
            @. krylov_basis[j + 1] -= H[i, j] * krylov_basis[i]
        end
        H[j + 1, j] = norm(krylov_basis[j + 1])
        krylov_basis[j + 1] ./= H[j + 1, j]

        # apply the previous Givens rotations to the new column of H
        @views H[1:j, j:j] .= Ω * H[1:j, j:j]

        # compute a new Givens rotation to zero out H[j + 1, j]
        G, _ = givens(H, j, j + 1, j)

        # apply the new rotation to H and the rhs
        H .= G * H
        g0 .= G * g0

        # compose the new rotation with the others
        Ω = lmul!(G, Ω)

        residual_norm = abs(g0[j + 1])

        if residual_norm < threshold
            converged = true
            break
        end
    end

    # solve the triangular system
    # The commented-out version occasionally introduces a ~150% increase in
    # allocations and a ~5% increase in time. When combined with the commented-
    # out kernel abstraction code below, this becomes a ~800% increase in
    # allocations and a ~200% increase in time.
    # y = NTuple{j}(@views UpperTriangular(H[1:j, 1:j]) \ g0[1:j])
    y = @views UpperTriangular(H[1:j, 1:j]) \ g0[1:j]

    # compose the solution
    # The commented-out version occasionally introduces a ~100% increase in time
    # without any increase in allocations. When combined with the commented-out
    # code above, this becomes a ~800% increase in allocations and a ~200%
    # increase in time.
    # rv_Q = realview(Q)
    # rv_krylov_basis = realview.(krylov_basis)
    # groupsize = 256
    # event = Event(array_device(Q))
    # event = linearcombination!(array_device(Q), groupsize)(
    #     rv_Q,
    #     y,
    #     rv_krylov_basis,
    #     true;
    #     ndrange = length(rv_Q),
    #     dependencies = (event,),
    # )
    # wait(array_device(Q), event)
    for i in 1:j
        @. Q += y[i] * krylov_basis[i]
    end

    # unwind right-preconditioning
    if isa(pc.pc_side, PCright)
        PCapply!(pc, cache.Wvec, Q, args...)
        @. Q = cache.Wvec
    end

    # If not converged, restart by reinitializing with current Q
    converged || LSinitialize!(solver, Q, Qrhs, args...)
    (converged, j, residual_norm)
end

max_stars_repo_name: CliMA/SystemSolvers

max_issues_repo_name: CliMA/SystemSolvers


text:
get_states_url() = "https://cdn-api.co-vin.in/api/v2/admin/location/states"
get_districts_url(state_id::Int) = "https://cdn-api.co-vin.in/api/v2/admin/location/districts/$(state_id)"
find_by_pin_url(pin::Int, date::String) = "https://cdn-api.co-vin.in/api/v2/appointment/sessions/public/findByPin?pincode=$(pin)&date=$(date)"
find_by_district_url(district_id::Int, date::String) ="https://cdn-api.co-vin.in/api/v2/appointment/sessions/public/findByDistrict?district_id=$(district_id)&date=$(date)"
calendar_by_pin_url(pin::Int, date::String) ="https://cdn-api.co-vin.in/api/v2/appointment/sessions/public/calendarByPin?pincode=$(pin)&date=$(date)"
calendar_by_district_url(district_id::Int, date::String) = "https://cdn-api.co-vin.in/api/v2/appointment/sessions/public/calendarByDistrict?district_id=$(district_id)&date=$(date)"
download_certificate_url(beneficiary_reference_id::String) = "https://cdn-api.co-vin.in/api/v2/registration/certificate/public/download?beneficiary_reference_id=$(beneficiary_reference_id)"

function formatted_date(date::Date)
    dt = Dates.format(date, "dd-mm-yyyy")
    return dt
end


text:
struct GaussianExpansion{T}
    n::Vector{Int}
    β::Vector{T}
    B::Vector{T}
end

Base.length(ge::GaussianExpansion) = length(ge.n)

Base.:(==)(a::GaussianExpansion, b::GaussianExpansion) =
    a.n == b.n && a.β == b.β && a.B == b.B

Base.hash(ge::GaussianExpansion, h::UInt) = hash((ge.n, ge.β, ge.B), h)

function (ge::GaussianExpansion{T})(r::AbstractVector) where T
    V = zeros(T, length(r))
    for k in 1:length(ge)
        V += ge.B[k] * r.^(ge.n[k]-2).*exp.(-ge.β[k]*r.^2)
    end
    V
end


text:
function solve(prob::NEqProblem, x, method::Anderson, options::NEqOptions)
    if !(mstyle(prob) === InPlace())
        throw(ErrorException("solve() not defined for OutOfPlace() with Anderson"))
    end
    function fixedfromnleq(F, x)
        F .= value(prob, F, x) .+ x
    end
    fixedpoint!(
        fixedfromnleq,
        x,
        method;
        # kwargs
        Gx = similar(x),
        Fx = similar(x),
        f_abstol = sqrt(eps(eltype(x))),
        maxiter = options.maxiter,
    )

end

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers