text:
"""
    abstract type AbstractElement{D,T}

Elements given by a fixed interpolation schemes mapping points on the the domain
`D<:AbstractReferenceShape` (of singleton type) to a value of type `T`.

Instances `el` of `AbstractElement` are expected to implement:
- `el(x̂)`: evaluate the interpolation scheme at the (reference) coordinate `x̂
  ∈ D`.
- `jacobian(el,x̂)` : evaluate the jacobian matrix of the interpolation at the
  (reference) coordinate `x ∈ D`.

!!! note
    For performance reasons, both `el(x̂)` and `jacobian(el,x̂)` should
    take as input a `StaticVector` and output a static vector or static array.
"""
abstract type AbstractElement{D,T} end

function (el::AbstractElement)(x)
    abstractmethod(el)
end

function jacobian(el,x)
    abstractmethod(el)
end

function WavePropBase.normal(el::AbstractElement, u)
    @assert u ∈ domain(el)
    jac = jacobian(el, u)
    normal(jac)
end

domain(::SType{<:AbstractElement{D}}) where {D <: AbstractReferenceShape} = D()

return_type(el::AbstractElement{D,T}) where {D,T} = T

domain_dimension(t::SType{<:AbstractElement}) = dimension(domain(t))

range_dimension(el::AbstractElement{R,T})        where {R,T} = length(T)

range_dimension(t::Type{<:AbstractElement{R,T}}) where {R,T} = length(T)

"""
    struct LagrangeElement{D,Np,T} <: AbstractElement{D,T}

Standard element over `D <: AbstractReferenceShape` commonly used in finite
element methods. The underlying polynomial space is [`Pk{D,K}`](@ref), and its
interpolant maps the `Np` `reference_nodes` in `D` to `Np` values of type `T`
stored in the field `vals`.
"""
struct LagrangeElement{D,Np,T} <: AbstractElement{D,T}
    vals::SVector{Np,T}
end

vals(el::LagrangeElement) = el.vals

"""
    reference_nodes(el::LagrangeElement)

Return the reference nodes on `domain(el)` used for the polynomial
interpolation. The function values on these nodes completely determines the
interpolating polynomial.

We use the same convention as `gmsh` for defining the reference nodes and their
order (see [node
ordering](https://gmsh.info/doc/texinfo/gmsh.html#Node-ordering) on `gmsh`
documentation).
"""
function reference_nodes end

# a contructor which infers the extra information from nodes
function LagrangeElement{R,K}(vals::SVector{Np,T}) where {R,K,T,Np}
    LagrangeElement{R,K,T}(vals)
end

# a contructor which infers the extra information from nodes
function LagrangeElement{D}(nodes::SVector{Np,T}) where {D,Np,T}
    LagrangeElement{D,Np,T}(nodes)
end

# constructor which converts each entry to a Point, and then creates an SVector
# of that.
function LagrangeElement{D}(vals) where {D}
    vals = SVector.(vals)
    LagrangeElement{D}(SVector(vals))
end

# a convenience constructor to allow things like LagrangeLine(a,b) instead of LagrangeLine((a,b))
function LagrangeElement{D}(nodes...) where {D}
    LagrangeElement{D}(nodes)
end

"""
    degree(el::LagrangeElement)

The polynomial degree of the element. A `LagrangeElement` of degree `K` and
domain `D` belongs to the space [`Pk{D,K}`](@ref).
"""
function degree(::SType{LagrangeElement{D,Np}})::Int where {D,Np}
    if D == ReferenceLine
        return Np-1
    elseif D == ReferenceTriangle
        K   = (-3 + sqrt(1+8*Np))/2
        return K
    elseif D == ReferenceTetrahedron
        notimplemented()
    elseif D == ReferenceSquare
        return sqrt(Np) - 1
    elseif D == ReferenceCube
        return Np^(1/3) - 1
    else
        notimplemented()
    end
end

function polynomial_space(::SType{LagrangeElement{D,Np}}) where {D,Np}
    K = degree(LagrangeElement{D,Np})
    Pk{D,K}()
end

"""
    const LagrangePoint{N,T} = LagrangeElement{ReferencePoint,1,SVector{N,T}}
"""
const LagrangePoint{N,T}        = LagrangeElement{ReferencePoint,1,SVector{N,T}}

"""
    const LagrangeLine = LagrangeElement{ReferenceLine}
"""
const LagrangeLine        = LagrangeElement{ReferenceLine}

"""
    const LagrangeTriangle = LagrangeElement{ReferenceTriangle}
"""
const LagrangeTriangle    = LagrangeElement{ReferenceTriangle}

"""
    const LagrangeTetrahedron = LagrangeElement{ReferenceTetrahedron}
"""
const LagrangeTetrahedron = LagrangeElement{ReferenceTetrahedron}

"""
    const LagrangeSquare = LagrangeElement{ReferenceSquare}
"""
const LagrangeSquare   = LagrangeElement{ReferenceSquare}

#=
Hardcode some basic elements.
TODO: Eventually this could/should be automated, at least for the LagrangeElements.
=#

# P1 for ReferenceLine
function (el::LagrangeLine{2})(u)
    @assert u ∈ ReferenceLine()
    v = vals(el)
    v[1] + (v[2] - v[1]) * u[1]
end

function jacobian(el::LagrangeLine{2}, u)
    @assert u ∈ ReferenceLine()
    v = vals(el)
    return hcat(v[2] - v[1])
end

# P2 for ReferenceLine
function (el::LagrangeLine{3})(u)
    @assert u ∈ domain(el)
    v = vals(el)
    v[1] + (4 * v[3] - 3 * v[1] - v[2]) * u[1]  + 2 * (v[2] + v[1] - 2 * v[3]) * u[1]^2
end

function jacobian(el::LagrangeLine{3}, u)
    @assert u ∈ domain(el)
    v = vals(el)
    hcat(4 * v[3] - 3 * v[1] - v[2] + 4 * (v[2] + v[1] - 2 * v[3]) * u[1])
end

# P1 for ReferenceTriangle
function (el::LagrangeTriangle{3})(u)
    @assert u ∈ domain(el)
    v = vals(el)
    v[1] + (v[2] - v[1]) * u[1] + (v[3] - v[1]) * u[2]
    # v[1]*(1-u[1]-u[2]) + v[2]*(u[1]) + v[3]*u[2]
end

function jacobian(el::LagrangeTriangle{3}, u)
    @assert u ∈ domain(el)
    v   = vals(el)
    jac = hcat( v[2] - v[1],
                v[3] - v[1])
    return jac
end

# P2 for ReferenceTriangle
function (el::LagrangeElement{ReferenceTriangle,6})(u)
    @assert u ∈ domain(el)
    v = vals(el)
    return (1+u[2]*(-3+2u[2])+u[1]*(-3+2u[1]+4u[2]))*v[1] +
           u[1]*(-v[2]+u[1]*(2v[2]-4v[4])+4v[4]+u[2]*(-4v[4]+4v[5]-4v[6])) +
           u[2]*(-v[3]+u[2]*(2v[3]-4v[6])+4v[6])
end

function jacobian(el::LagrangeElement{ReferenceTriangle,6}, u)
    @assert u ∈ domain(el)
    v = vals(el)
    hcat(
        (-3+4u[1]+4u[2])*v[1] - v[2] + u[1]*(4v[2]-8v[4]) + 4v[4] + u[2]*(-4v[4]+4v[5]-4v[6]),
        (-3+4u[1]+4u[2])*v[1] - v[3] + u[2]*(4v[3]-8v[6]) + u[1]*(-4v[4]+4v[5]-4v[6]) + 4v[6]
    )
end

# P1 for ReferenceSquare
function (el::LagrangeElement{ReferenceSquare,4})(u)
    @assert u ∈ domain(el)
    v = vals(el)
    v[1] + (v[2] - v[1]) * u[1] + (v[4] - v[1]) * u[2] + (v[3] + v[1] - v[2] - v[4]) * u[1] * u[2]
end

function jacobian(el::LagrangeElement{ReferenceSquare,4}, u)
    @assert u ∈ domain(el)
    v = vals(el)
    hcat(
        ((v[2] - v[1]) + (v[3] + v[1] - v[2] - v[4]) * u[2]),
        ((v[4] - v[1]) + (v[3] + v[1] - v[2] - v[4]) * u[1])
    )
end

# P1 for ReferenceTetrahedron
function (el::LagrangeElement{ReferenceTetrahedron,4})(u)
    @assert u ∈ domain(el)
    v = vals(el)
    v[1] + (v[2] - v[1]) * u[1] + (v[3] - v[1]) * u[2] + (v[4] - v[1]) * u[3]
end

function jacobian(el::LagrangeElement{ReferenceTetrahedron,4}, u)
    @assert u ∈ domain(el)
    v = vals(el)
    hcat(
        (v[2] - v[1]),
        (v[3] - v[1]),
        (v[4] - v[1])
    )
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
using CuArrays, CUDAnative, GPUArrays
using NiLang, NiLang.AD

"""
A reversible swap kernel for GPU for SWAP gate in quantum computing.
See the irreversible version for comparison

http://tutorials.yaoquantum.org/dev/generated/developer-guide/2.cuda-acceleration/
"""
@i @inline function swap_kernel(state::AbstractVector{T}, mask1, mask2) where T
    @invcheckoff b ← (blockIdx().x-1) * blockDim().x + threadIdx().x
    @invcheckoff if (b < length(state), ~)
        if (b&mask1==0 && b&mask2==mask2, ~)
            NiLang.SWAP(state[b+1], state[b ⊻ (mask1|mask2) + 1])
        end
    end
end

# TODO: support ::Type like argument.
"""
SWAP gate in quantum computing.
"""
@i function instruct!(state::CuVector, gate::Val{:SWAP}, locs::Tuple{Int,Int})
    mask1 ← 1 << (tget(locs, 1)-1)
    mask2 ← 1 << (tget(locs, 2)-1)
    XY ← GPUArrays.thread_blocks_heuristic(length(state))
    @cuda threads=tget(XY,1) blocks=tget(XY,2) swap_kernel(state, mask1, mask2)
end

using Test
@testset "swap gate" begin
    v = cu(randn(128))
    v1 = instruct!(copy(v), Val(:SWAP), (3,4))[1]
    v2 = instruct!(copy(v1), Val(:SWAP), (3,4))[1]
    v3 = (~instruct!)(copy(v1), Val(:SWAP), (3,4))[1]
    @test !(v ≈ v1)
    @test v ≈ v2
    @test v ≈ v3
end

@i function loss(out!, state::CuVector)
    instruct!(state, Val(:SWAP), (3,4))
    out! += state[4]
end

loss(0.0, CuArray(randn(128)))
Grad(loss)(Loss(0.0), CuArray(randn(128)))

####################### A different loss ###############
@i function loss(out!, state::CuVector, target::CuVector)
    instruct!(state, Val(:SWAP), (3,4))
    out! += state' * target
end

# requires defining a new primitive, we don't how to parallelize a CUDA program automatically yet.
using LinearAlgebra: Adjoint
function (_::MinusEq{typeof(*)})(out!::GVar, x::Adjoint{<:Any, <:CuVector{<:GVar}}, y::CuVector{<:GVar})
    chfield(out!, value, value(out!)-(value.(x) * value.(y))[]),
    chfield.(parent(x), grad, grad.(parent(x)) .+ grad(out!)' .* conj.(value.(y)))',
    chfield.(y, grad, grad.(y) .+ grad(out!) .* conj.(value.(x')))
end

function (_::PlusEq{typeof(*)})(out!::GVar, x::Adjoint{<:Any, <:CuVector{<:GVar}}, y::CuVector{<:GVar})
    chfield(out!, value, value(out!)+(value.(x) * value.(y))[]),
    chfield.(parent(x), grad, grad.(parent(x)) .- grad(out!)' .* conj.(value.(y)))',
    chfield.(y, grad, grad.(y) .- grad(out!) .* conj.(value.(x')))
end

function (_::PlusEq{typeof(*)})(out!, x, y)
    out! += x * y
    out!, x, y
end

function (_::MinusEq{typeof(*)})(out!, x, y)
    out! -= x * y
    out!, x, y
end

loss(0.0, CuArray(randn(128)), CuArray(randn(128)))
loss'(Loss(0.0), CuArray(randn(128)), CuArray(randn(128)))

\begin{minipage}{.88\columnwidth}
\begin{lstlisting}
using CuArrays, CUDAnative, GPUArrays
using NiLang, NiLang.AD

@i @inline function swap_kernel(state, mask1, mask2)
    @invcheckoff begin
        b ← (blockIdx().x-1) *
            blockDim().x + threadIdx().x - 1
        if (b < length(state), ~)
            if (b&mask1==0 && b&mask2==mask2, ~)
                SWAP(state[b+1], state[b ⊻
                    (mask1|mask2) + 1])
            end
        end
        b → (blockIdx().x-1) *
            blockDim().x + threadIdx().x - 1
    end
end
\end{lstlisting}
\end{minipage}

This kernel function simulates the SWAP gate in quantum computing.
Here, one must use the macro \texttt{@invcheckoff} to turn off the reversibility checks. It is necessary because the possible error thrown in a kernel function can not be handled on a CUDA kernel.
One can launch this kernel function to GPUs with a single macro \texttt{@cuda}, as shown in the following using case.

\begin{minipage}{\columnwidth}
\begin{lstlisting}[multicols=2]
julia> @i function instruct!(state::CuVector,
            gate::Val{:SWAP}, locs::Tuple{Int,Int})
           mask1 ← 1 << (tget(locs, 1)-1)
           mask2 ← 1 << (tget(locs, 2)-1)
           XY ← GPUArrays.thread_blocks_heuristic(
                length(state))
           @cuda threads=tget(XY,1) blocks=tget(XY,
                2) swap_kernel(state, mask1, mask2)
       end

julia> instruct!(CuArray(randn(8)),
            Val(:SWAP), (1,3))[1]
8-element CuArray{Float64,1,Nothing}:
 -0.06956048379200473
 -0.6464176838567472
 -0.06523362834285944
 -0.7314356941903547
  1.512329204247244
  0.9773772766637732
  1.6473223915215722
 -1.0631789613639087
\end{lstlisting}
\end{minipage}

One can also write kernels with KernelAbstaction. It solves many compatibility issues related to different function calls on GPU and CPU.

\begin{minipage}{.88\columnwidth}
\begin{lstlisting}
@i @kernel function swap_kernel2(state, mask1, mask2)
    @invcheckoff begin
        b ← @index(Global)
        if (b < length(state), ~)
            if (b&mask1==0 && b&mask2==mask2, ~)
                SWAP(state[b+1], state[b ⊻
                    (mask1|mask2) + 1])
            end
        end
        b → @index(Global)
    end
end
\end{lstlisting}
\end{minipage}

We can use the macro \texttt{@launchkernel} to launch a kernel.
The first parameter is a device.
The second parameter is the block size.
The third parameter is the number of threads.
The last parameter is a kernel function call to be launched.

\begin{minipage}{\columnwidth}
\begin{lstlisting}[multicols=2]
julia> @i function instruct!(state::CuVector,
            gate::Val{:SWAP}, locs::Tuple{Int,Int})
           mask1 ← 1 << (tget(locs, 1)-1)
           mask2 ← 1 << (tget(locs, 2)-1)
           XY ← GPUArrays.thread_blocks_heuristic(
                length(state))
           @launchkernel CUDA() 256 length(out!
                ) swap_kernel2(state, mask1, mask2)
       end

julia> instruct!(CuArray(randn(8)),
            Val(:SWAP), (1,3))[1]
8-element CuArray{Float64,1,Nothing}:
  2.1492759883720525
  2.326837084303501
  1.4587667131427016
 -1.3273806428138293
 -0.03975355575683114
 -0.10763082744447787
 -1.7111718557581195
 -0.47922613687722704
\end{lstlisting}
\end{minipage}

ext: jl


text:
using Revise, WaltonForaging
ispath("/home/beatriz/Documents/") ? (main_path ="/home/beatriz/Documents/") : (main_path = "/Users/dariosarra/Documents/Lab/Walton/WaltonForaging")
fig_dir = joinpath(main_path,"Figures")
# pokes = DataFrame(CSV.read(joinpath(main_path,"Miscellaneous/Poke_data.csv")))#, DataFrame)
# pokes2 = DataFrame(CSV.read(joinpath(main_path,"Miscellaneous/New_Poke_data.csv")))#, DataFrame)
pokes = CSV.read(joinpath(main_path,"Miscellaneous/States_Poke_data.csv"), DataFrame)
##
Exp = "DAphotometry"#"5HTPharma"
FileList = readdir(joinpath(main_path,Exp,"RawData"))
filter!(r -> ismatch(r".txt",r), FileList)
session = joinpath(main_path,Exp,"RawData",FileList[1])
#=
Depending on when the p line is written, during forge port poke or during travel port poke,
T  = poke time required to get reward (retrospective?), poke time required to travel
    (initial value to threshold 2000)
IFT = always tells you the initial average of the dist for which the richness is taken
AFT = during reward collection average dist for that specific reward (the updated average),
    during travel for how long in the previous forage before leaving
There is a timeout for going back to travel if animal don't engage in foraging
=#
##
session
pokes = process_raw_session(session; observe = true)
combine(groupby(pokes,:Patch), :RewardAvailable => x -> sum(.!ismissing.(x)))
bouts = process_bouts(pokes; observe = true)
##
findfirst(FileList.== "fp14-2019-05-30-172855.txt")
AllBouts = DataFrame()
for file in FileList[1:end]
    session = joinpath(main_path,Exp,"RawData",file)
    pokes = process_raw_session(session; observe = false)
    bouts = process_bouts(pokes; observe = false)
    # path = joinpath(main_path,Exp,"Processed",replace(file,"txt"=>"csv"))
    # CSV.write(path, bouts)
    if isempty(AllBouts)
        AllBouts = bouts
        allowmissing!(AllBouts)
    else
        append!(AllBouts,bouts)
    end
    println(file)
end
path = joinpath(main_path,Exp,"Processed","AllBouts.csv")
CSV.write(path, AllBouts)
##
nrow(AllBouts)
check_idx = findfirst(AllBouts.SubjectID .== "fp17" .&&
    AllBouts.Startdate .== "2019/06/2015:06:16" .&&
    AllBouts.Patch .== 53)

open_html_table(AllBouts[check_idx-10:check_idx+10,:])
findfirst(ismissing(AllBouts.In))
##
test = filter(r -> r.State == "Forage" && r.GiveUp && !ismissing(r.ForageTime_Sum), AllBouts)
test.Richness = categorical(test.Richness)
levels!(test.Richness,["poor","medium","rich"])
test.Travel = categorical(test.Travel)
levels!(test.Travel,["Short","Long"])
df1 = combine(groupby(test,[:Travel,:Richness,:SubjectID]),:ForageTime_Sum .=> mean => :ForageTime)
df2 = combine(groupby(df1,[:Travel,:Richness]),:ForageTime .=> [mean, sem])
@df df2 scatter(:Richness,:ForageTime_mean,group = :Travel,yerror=:ForageTime_sem,
    ylims=(1000,3000), legend = :outerright,
    ylabel = "Average forage time \n from last reward at leaving", xlabel = "Patch richness",
    legendtitle = "Travel type")
savefig("/Users/dariosarra/Documents/Lab/Walton/WaltonForaging/DAphotometry/CoarseAnalysis_LastBout.pdf")
##
test = filter(r -> r.State == "Forage", AllBouts)
test.Richness = categorical(test.Richness)
levels!(test.Richness,["poor","medium","rich"])
test.Travel = categorical(test.Travel)
levels!(test.Travel,["Short","Long"])
df1 = combine(groupby(test,[:Travel,:Richness,:Patch,:SubjectID,:Startdate]),:ForageTime_Sum .=>
    (x->sum(skipmissing(x))) => :ForageTime)
df2 = combine(groupby(df1,[:Travel,:Richness,:SubjectID]),:ForageTime .=> mean => :ForageTime)
df3 = combine(groupby(df2,[:Travel,:Richness]),:ForageTime .=> [mean, sem])
@df df3 scatter(:Richness,:ForageTime_mean, group = :Travel, yerror =:ForageTime_sem,
    ylims = (2000,9000), legend = :outerright,
    ylabel = "Average forage time \n from last reward at leaving", xlabel = "Patch richness",
    legendtitle = "Travel type")
savefig("/Users/dariosarra/Documents/Lab/Walton/WaltonForaging/DAphotometry/CoarseAnalysis2_FullPatch.pdf")


text:
function reset!(sim::Simulation)
    sim.t = 0
    sim.queue = Vector{Function}()

    empty!(sim.stt.current_stock)
    empty!(sim.stt.dispatch_queue)

    empty!(sim.acc.executed_dispatch)
    empty!(sim.acc.inventory_levels)
    empty!(sim.acc.demand_fulfillment)

    return sim
end

function initiate!(sim::Simulation; reset::Bool=true)
    reset && reset!(sim) # reset simulation
    if (isempty(sim.stt) || isempty(sim.acc))
        append!(sim.stt.current_stock, sim.libs.init_stock)
        append!(sim.acc.inventory_levels,
            insertcols!(
                sim.libs.init_stock |> copy, :id, :periode => fill(sim.t), after=true
            )
        )
    else
        println("Simulation already initiated!")
    end
    return sim
end

function schedule!(sim::Simulation)
    sim.t < sim.params.T && begin
        push!(sim.queue, plan!)
        push!(sim.queue, transport!)
        push!(sim.queue, fulfill!)
    end

    return nothing
end

function run!(sim::Simulation)
    (isempty(sim.stt) || isempty(sim.acc)) && error("Simulation uninitiated!")

    println("Starting simulation")
    start = time()
    while sim.t < sim.params.T
        if sim.fixed
            sim.params.H = min(sim.params.H, sim.params.T - sim.t)
        end
        println(sim) # verbose
        schedule!(sim)
        while !isempty(sim.queue)
            fn = popfirst!(sim.queue)
            sim |> fn
        end
        sim.t += 1
    end
    stop = time()

    println("Simulation stopped. Total duration: $(stop-start)")

    return nothing
end


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

using Documenter, Sunny
# Need to import libraries which expose functions hidden behind Requires.jl
using GLMakie

makedocs(
    sitename="SU(N) Spin Simulations",
    pages = [
        "index.md",
        "Getting Started" => "getting-started.md",
        "Examples" => "examples.md",
        "Library" => "library.md",
        "Structure Factor Calculations" => "structure-factor.md",
        "Internals" => "internals.md",
    ],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    )
)

deploydocs(
    repo = "github.com/SunnySuite/Sunny.jl.git",
    devbranch = "main",
)


text:
#!/usr/bin/env julia --project
using ClimateMachine
ClimateMachine.init(parse_clargs = true)

using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.ConfigTypes
using ClimateMachine.Diagnostics
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.TurbulenceClosures
using ClimateMachine.SystemSolvers: ManyColumnLU
using ClimateMachine.Mesh.Filters
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Interpolation
using Thermodynamics.TemperatureProfiles
using Thermodynamics: total_energy, air_density
using ClimateMachine.VariableTemplates

using Distributions: Uniform
using LinearAlgebra
using StaticArrays

using CLIMAParameters
using CLIMAParameters.Planet:
    R_d, day, grav, cp_d, planet_radius, Omega, kappa_d, MSLP
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

import CLIMAParameters
CLIMAParameters.Planet.Omega(::EarthParameterSet) = 0.0
CLIMAParameters.Planet.planet_radius(::EarthParameterSet) = 6.371e6 / 125.0
CLIMAParameters.Planet.MSLP(::EarthParameterSet) = 1e5

function init_nonhydrostatic_gravity_wave!(problem, bl, state, aux, localgeo, t)
    FT = eltype(state)

    # grid
    φ = latitude(bl, aux)
    λ = longitude(bl, aux)
    z = altitude(bl, aux)

    # parameters
    param_set = parameter_set(bl)
    _grav::FT = grav(param_set)
    _cp::FT = cp_d(param_set)
    _Ω::FT = Omega(param_set)
    _a::FT = planet_radius(param_set)
    _R_d::FT = R_d(param_set)
    _kappa::FT = kappa_d(param_set)
    _p_eq::FT = MSLP(param_set)

    N::FT = 0.01
    u_0::FT = 0.0
    G::FT = _grav^2 / N^2 / _cp
    T_eq::FT = 300
    Δθ::FT = 0.0
    d::FT = 5e3
    λ_c::FT = 2 * π / 3
    φ_c::FT = 0
    L_z::FT = 20e3

    # initial velocity profile (we need to transform the vector into the Cartesian
    # coordinate system)
    u_sphere = SVector{3, FT}(u_0 * cos(φ), 0, 0)
    u_init = sphr_to_cart_vec(bl.orientation, u_sphere, aux)

    # background temperature
    T_s::FT =
        G +
        (T_eq - G) *
        exp(-u_0 * N^2 / 4 / _grav^2 * (u_0 + 2 * _Ω * _a) * (cos(2 * φ) - 1))
    T_b::FT = G * (1 - exp(N^2 / _grav * z)) + T_s * exp(N^2 / _grav * z)

    # pressure
    p_s::FT =
        _p_eq *
        exp(u_0 / 4 / G / _R_d * (u_0 + 2 * _Ω * _a) * (cos(2 * φ) - 1)) *
        (T_s / T_eq)^(1 / _kappa)
    p::FT = p_s * (G / T_s * exp(-N^2 / _grav * z) + 1 - G / T_s)^(1 / _kappa)

    # background potential temperature
    θ_b::FT = T_b * (_p_eq / p)^_kappa

    # potential temperature perturbation
    r::FT = _a * acos(sin(φ_c) * sin(φ) + cos(φ_c) * cos(φ) * cos(λ - λ_c))
    s::FT = d^2 / (d^2 + r^2)
    θ′::FT = Δθ * s * sin(2 * π * z / L_z)

    # temperature perturbation
    T′::FT = θ′ * (p / _p_eq)^_kappa

    # temperature
    T::FT = T_b + T′

    # density
    ρ = air_density(param_set, T_b, p)

    # potential & kinetic energy
    e_pot = gravitational_potential(bl.orientation, aux)
    e_kin::FT = 0.5 * sum(abs2.(u_init))

    state.ρ = ρ
    state.ρu = ρ * u_init
    state.energy.ρe = ρ * total_energy(param_set, e_kin, e_pot, T)

    nothing
end

function config_nonhydrostatic_gravity_wave(FT, poly_order, resolution)
    # Set up a reference state for linearization of equations
    temp_profile_ref =
        DecayingTemperatureProfile{FT}(param_set, FT(300), FT(100), FT(27.5e3))
    ref_state = HydrostaticState(temp_profile_ref)

    domain_height::FT = 10e3               # distance between surface and top of atmosphere (m)

    # Set up the atmosphere model
    exp_name = "NonhydrostaticGravityWave"

    physics = AtmosPhysics{FT}(
        param_set;
        ref_state = ref_state,
        turbulence = ConstantKinematicViscosity(FT(0)),
        moisture = DryModel(),
    )

    model = AtmosModel{FT}(
        AtmosGCMConfigType,
        physics;
        init_state_prognostic = init_nonhydrostatic_gravity_wave!,
        source = (Gravity(),),
    )

    config = ClimateMachine.AtmosGCMConfiguration(
        exp_name,
        poly_order,
        resolution,
        domain_height,
        param_set,
        init_nonhydrostatic_gravity_wave!;
        model = model,
    )

    return config
end

function config_diagnostics(FT, driver_config)
    interval = "40000steps" # chosen to allow a single diagnostics collection

    _planet_radius = FT(planet_radius(param_set))

    info = driver_config.config_info
    boundaries = [
        FT(-90.0) FT(-180.0) _planet_radius
        FT(90.0) FT(180.0) FT(_planet_radius + info.domain_height)
    ]
    resolution = (FT(10), FT(10), FT(100)) # in (deg, deg, m)
    interpol = ClimateMachine.InterpolationConfiguration(
        driver_config,
        boundaries,
        resolution,
    )

    dgngrp = setup_atmos_default_diagnostics(
        AtmosGCMConfigType(),
        interval,
        driver_config.name,
        interpol = interpol,
    )

    return ClimateMachine.DiagnosticsConfiguration([dgngrp])
end

function main()
    # Driver configuration parameters
    FT = Float64                             # floating type precision
    poly_order = 5                           # discontinuous Galerkin polynomial order
    n_horz = 8                               # horizontal element number
    n_vert = 4                               # vertical element number
    timestart = FT(0)                        # start time (s)
    timeend = FT(3600)                       # end time (s)

    # Set up driver configuration
    driver_config =
        config_nonhydrostatic_gravity_wave(FT, poly_order, (n_horz, n_vert))

    # Set up experiment
    CFL = FT(0.4)
    solver_config = ClimateMachine.SolverConfiguration(
        timestart,
        timeend,
        driver_config,
        Courant_number = CFL,
        CFL_direction = HorizontalDirection(),
    )

    # Set up diagnostics
    dgn_config = config_diagnostics(FT, driver_config)

    # Set up user-defined callbacks
    filterorder = 64
    filter = ExponentialFilter(solver_config.dg.grid, 0, filterorder)
    cbfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            AtmosFilterPerturbations(driver_config.bl),
            solver_config.dg.grid,
            filter,
            state_auxiliary = solver_config.dg.state_auxiliary,
        )
        nothing
    end

    # Run the model
    result = ClimateMachine.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        user_callbacks = (cbfilter,),
        check_euclidean_distance = true,
    )
end

main()


text: @test Id(2) == Symp([1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1])


text:
# using Optim
# using PyPlot
# using FFTW
# using ACME
# using ElectricalEngineering

# include("Load_LTSpice_Net.jl")
# include("RunAnalysis.jl")
# include("SPICE2Matrix.jl")
# include("ToroidOptimizer.jl")

"""
This function takes in similar inputs to DesignDriveFilter. It will iterate that function with varying matching impedances to minimize the expected drift in the current.
"""
function DesignDriveFilter_OptimDrift(
    LDrive,
    RDrive,
    TargetZ,
    DriveFreq;
    CDrive = 1e6,
    NumDriveElements = 1,
    WireDiam = 2e-3,
    WireFillFac = 0.75,
    PlotSrcR = TargetZ,
    PlotFTs = true,
    VSrc = 2,
    DetermineFreq = false,
    AddNotchFreq = nothing,
    FilterZ = TargetZ,
    RDampVal = FilterZ
)

    function WrapperFunk(Zin)
    DriveFreq, CurrentVec, Results, SPICE_DF,inputs,InputList,FreqList = DesignDriveFilter(
        LDrive,RDrive, Zin,DriveFreq;
        CDrive = CDrive,
        NumDriveElements = NumDriveElements,
        WireDiam = WireDiam,
        DetermineFreq=DetermineFreq,
        AddNotchFreq = AddNotchFreq,
        FilterZ = Zin)

        CircModDesign.DetermineComponentsTempCoeffs(SPICE_DF,InputList,1,DriveFreq,"LDrive")

        return SPICE_DF.DriftCoeff[findfirst(isequal("CDrive"),SPICE_DF.Name)]
    end

    MinZ = optimize(WrapperFunk,10,15)

    DriveFreq, CurrentVec, Results, SPICE_DF,inputs,InputList,FreqList = DesignDriveFilter(
        LDrive,RDrive, MinZ.minimizer,DriveFreq;
        CDrive = CDrive,
        WireDiam = WireDiam,
        DetermineFreq=DetermineFreq,
        AddNotchFreq = AddNotchFreq,
        FilterZ = MinZ.minimizer)

        DetermineComponentsTempCoeffs(SPICE_DF,InputList,1,DriveFreq,"LDrive")
    return MinZ.minimizer, DriveFreq, CurrentVec, Results, SPICE_DF,inputs,InputList,FreqList
end

"""
This function takes in  the following parameters:

    LDrive, drive coil inductance (H)

    RDrive, Drive coil resistance (Ω)

    TargetZ, Amplifier/target impedance to match to  (Ω)

    DriveFreq; drive frequency (Hz)

    CDrive = 1e6, Any series capacitance with drive coil

    NumDriveElements = 1, The number of series copies of the L,C,R for example if there are 6 pancakes of LDrive, RDrive, Cdrive, then this would be 6

    WireDiam = 2e-3, Wire diameter for toroids

    WireFillFac = 0.75, If litz wire, this is the ratio of copper to air

    Example (Copy/paste):
    RDrive = 400e-3
    LDrive = 400e-6
    Amp_Z = 4
    DriveFreq = 25e3

    DesignDriveFilter(LDrive,RDrive,Amp_Z,DriveFreq)

"""
function DesignDriveFilter(
    LDrive,
    RDrive,
    TargetZ,
    DriveFreq;
    CDrive = 1e6,
    NumDriveElements = 1,
    WireDiam = 2e-3,
    WireFillFac = 0.75,
    PlotSrcR = TargetZ,
    PlotFTs = true,
    VSrc = 2,
    DetermineFreq = false,
    AddNotchFreq = nothing,
    FilterZ = TargetZ,
    RDampVal = FilterZ
)

    ZeroVal = 1e-9 #There are issues with zero-valued components
    ωDr = 2 * π * DriveFreq
    ZDrive =
        RDrive * NumDriveElements +
        im * ωDr * LDrive * NumDriveElements +
        NumDriveElements ./ (im * 2 * pi * DriveFreq * CDrive) #Equivalent complex impednace of the load

    Reactance_Load =
        ωDr * LDrive * NumDriveElements - NumDriveElements ./ (ωDr * CDrive)
    Reactance_Load = round(Reactance_Load; sigdigits = 3)
    println("The reactance of the load is: $(round(Reactance_Load;sigdigits=3)) Ω")
    if ~DetermineFreq
        if (Reactance_Load > 0)

            SerCap,CParAct = ImpMatch_LLoad(TargetZ, ZDrive, Reactance_Load,ωDr)

            LTee_2 = ZeroVal
            LTee_2_ESR = ZeroVal
            LTee_1 = ZeroVal
            LTee_1_ESR = ZeroVal

        elseif (Reactance_Load < 0)

            LTee_2,LTee_2_ESR,CParAct = ImpMatch_CLoad(TargetZ, ZDrive, Reactance_Load,ωDr)

            LTee_1 = ZeroVal
            LTee_1_ESR = ZeroVal
            SerCap = ZeroVal

        elseif (Reactance_Load == 0)
            matchRatio = real(TargetZ) / real(ZDrive)
            Q = sqrt(matchRatio - 1)
            Xs = Q * real(ZDrive)
            LSer2 = Xs / (ωDr)
            LTee_2 = LSer2
            LTee_1 = ZeroVal
            LTee_1_ESR = ZeroVal
            CParAct = findResPair((1 + Q^(-2)) * LSer2, DriveFreq)
            LTee_2_Geom =
                ToroidOptimizer(WireDiam, LTee_2; CuFillFactor = WireFillFac)
            LTee_2_ESR = LTee_2_Geom.DCore.Resistance
            SerCap = 1e6
        end
    else

        matchRatio = real(TargetZ) / real(ZDrive)
        Q = √(matchRatio - 1)
        Xs = Q * real(ZDrive)

        DriveFreq = findFilterFreq(Q,real(ZDrive), LDrive * NumDriveElements, CDrive / NumDriveElements)
        ωDr = 2*π*DriveFreq
        println("Determined Q to be $(round(Q))")
        println("Determined drive freq. to be $(round(DriveFreq)) Hz")
        ZSer = 1*im * ωDr*LDrive * NumDriveElements -  1*im /(ωDr* CDrive / NumDriveElements) +RDrive* NumDriveElements
        YSer = 1/ZSer

        CParAct = abs.(imag(YSer))/ωDr
        LTee_2 = ZeroVal
        LTee_2_ESR = ZeroVal
        LTee_1 = ZeroVal
        LTee_1_ESR = ZeroVal
        SerCap = 1e6
    end
    println("L Tee 1 =  $(round(LTee_1*1e6;sigdigits=3))μH ")
    println("L Tee 2 =  $(round(LTee_2*1e6;sigdigits=3))μH ")
    println("CParAct =  $(round(CParAct*1e6;sigdigits=3))μF ")
    println("SerCap =  $(round(SerCap*1e6;sigdigits=3))μF ")

    (LFilt, CFilt) = Butterworth_2(FilterZ, DriveFreq)
    LFiltMatch_C = findResPair(LFilt, DriveFreq)
    CFiltMatch_L = findResPair(CFilt, DriveFreq)

    LFilt1_Geom =
            ToroidOptimizer(WireDiam, LFilt; CuFillFactor = WireFillFac)
    LFilt1_ESR = LFilt1_Geom.DCore.Resistance

    LFilt2_Geom =
            ToroidOptimizer(WireDiam, CFiltMatch_L; CuFillFactor = WireFillFac)
    LFilt2_ESR = LFilt2_Geom.DCore.Resistance
    println(
        "LFilt =  $(round(LFilt*1e6;sigdigits=3))μH matched with: LFiltMatch_C =  $(round(LFiltMatch_C*1e6;sigdigits=3))μF ",
    )
    println(
        "CFilt =  $(round(CFilt*1e6;sigdigits=3))μF matched with: CFiltMatch_L = $(round(CFiltMatch_L*1e6;sigdigits=3))μH ",
    )

    if ~(AddNotchFreq===nothing)
        LNotch, CNotch, LNotch_Tune = makeNotchSection(AddNotchFreq, DriveFreq, TargetZ;LVal = 100e-6)
        LNotch_Geom =
                ToroidOptimizer(WireDiam, LNotch; CuFillFactor = WireFillFac)
        LNotch_ESR = LNotch_Geom.DCore.Resistance
        LNotch_Tune_Geom =
                ToroidOptimizer(WireDiam, LNotch_Tune; CuFillFactor = WireFillFac)
        LNotch_Tune_ESR = LNotch_Tune_Geom.DCore.Resistance
    else
        LNotch = 10
        CNotch = 1e-12
        LNotch_Tune = 1e-12
        LNotch_ESR = 1e-3
        LNotch_Tune_ESR = 1e-3
    end

    CurrentVec, Results, SPICE_DF,inputs,InputList,FreqList = CircModel_SPICE(DriveFreq, VSrc,
    1e-3,
    RDrive,
    NumDriveElements,
    LDrive,
    CDrive,
    SerCap,
    LTee_2,
    LTee_2_ESR,
    LTee_1,
    LTee_1_ESR,
    CParAct,
    LFilt,
    LFilt2_ESR,
    LFilt1_ESR,
    LFiltMatch_C,
    CFilt,
    CFiltMatch_L,
    LNotch,
    LNotch_ESR,
    CNotch,
    LNotch_Tune,
    LNotch_Tune_ESR;
    PlotOn = PlotFTs,
    RDampVal = RDampVal)

    return DriveFreq, CurrentVec, Results, SPICE_DF,inputs,InputList,FreqList

end

function CircModel_SPICE(DriveFreq, VSrc,
    PlotSrcR,
    RDrive,
    NumDriveElements,
    LDrive,
    CDrive,
    SerCap,
    LTee_2,
    LTee_2_ESR,
    LTee_1,
    LTee_1_ESR,
    CParAct,
    LFilt,
    LFilt2_ESR,
    LFilt1_ESR,
    LFiltMatch_C,
    CFilt,
    CFiltMatch_L,
    LNotch,
    LNotch_ESR,
    CNotch,
    LNotch_Tune,
    LNotch_Tune_ESR;
    PlotOn = true,
    FreqList = 1000:10:100e3,
    ArchetypeNetFileName = nothing,
    RDampVal = nothing)

    if ArchetypeNetFileName === nothing

        ArchetypeNetFileName = joinpath(dirname(pathof(CircModDesign)),"Filter_Archetype_Damped_2.net")

    end
    SPICE_DF,NodeList,InputList,NumVSources = SPICE2Matrix(ArchetypeNetFileName)
    UpdateElementVal!(SPICE_DF,"RSrc",PlotSrcR)

    UpdateElementVal!(SPICE_DF,"LNotch_Tune",LNotch_Tune)
    UpdateElementVal!(SPICE_DF,"LNotch",LNotch)
    UpdateElementESR!(SPICE_DF,"LNotch_Tune",LNotch_Tune_ESR)
    UpdateElementESR!(SPICE_DF,"LNotch",LNotch_ESR)
    UpdateElementVal!(SPICE_DF,"CNotch",CNotch)

    UpdateElementVal!(SPICE_DF,"LFilt1",LFilt)
    UpdateElementVal!(SPICE_DF,"LFilt2",LFilt)
    UpdateElementVal!(SPICE_DF,"CLFilt_C1",LFiltMatch_C)
    UpdateElementVal!(SPICE_DF,"CLFilt_C2",LFiltMatch_C)
    UpdateElementESR!(SPICE_DF,"LFilt1",LFilt1_ESR)

    UpdateElementVal!(SPICE_DF,"CFilt",CFilt)
    UpdateElementVal!(SPICE_DF,"LCFilt_L",CFiltMatch_L)
    UpdateElementESR!(SPICE_DF,"LCFilt_L",LFilt2_ESR)

    UpdateElementVal!(SPICE_DF,"RDrive",RDrive)
    UpdateElementVal!(SPICE_DF,"LDrive",NumDriveElements*LDrive)
    UpdateElementVal!(SPICE_DF,"CDrive",CDrive/NumDriveElements)

    UpdateElementVal!(SPICE_DF,"CSer",SerCap)
    UpdateElementVal!(SPICE_DF,"CPar",CParAct)

    UpdateElementVal!(SPICE_DF,"LTee1",LTee_1)
    UpdateElementESR!(SPICE_DF,"LTee1",LTee_1_ESR)
    UpdateElementVal!(SPICE_DF,"LTee2",LTee_2)
    UpdateElementESR!(SPICE_DF,"LTee2",LTee_2_ESR)

    if ~(RDampVal===nothing)

        UpdateElementVal!(SPICE_DF,"LDamp",LFilt)
        UpdateElementESR!(SPICE_DF,"LDamp",LFilt1_ESR)
        UpdateElementVal!(SPICE_DF,"CDamp",LFiltMatch_C)
        UpdateElementVal!(SPICE_DF,"RDamp2",RDampVal)
        UpdateElementVal!(SPICE_DF,"RDamp1",RDampVal)

    end

    inputs = zeros(length(InputList))
    inputs[end-(NumVSources-1):end] .= 1

    ResultNodeNames = vcat("V(".*NodeList.*")", "I(".*(InputList[end-(NumVSources-1):end]).*")")

    Results = [ abs.(inv(SPICE_DF2Matrix_ω(SPICE_DF,2*π*FreqList[i],InputList))*inputs) for i in 1:length(FreqList)]
    Results = hcat(Results...)
    ResDict = Dict(ResultNodeNames[1] => Results[1,:])
    for i in 2:length(ResultNodeNames)

         merge!(ResDict,Dict(ResultNodeNames[i] => Results[i,:]))

    end

    CurrentVec = VSrc.* plotACElCurrent(SPICE_DF,FreqList,Results,"LDrive")

    if PlotOn
        pygui(true)
        semilogy(FreqList,abs.(CurrentVec[:]))

        xlabel("Frequency, Hz")
        ylabel("Current in LDrive")
    end

    getElementCurrents(SPICE_DF,Results,DriveFreq)
    return CurrentVec, Results, SPICE_DF,inputs,InputList,FreqList
end

function ImpMatch_LLoad(TargetZ, ZDrive, Reactance_Load,ωDr)
    println("Load is inductive")
    EquivSerL = Reactance_Load / ωDr
    println(
        "Load appears to be a inductor with a value of $(round(EquivSerL*1e6;sigdigits=3))μH and a ESR of $(real(ZDrive))Ω ",
    )

    matchRatio = real(TargetZ) / real(ZDrive)
    Q = sqrt(matchRatio - 1)
    Xs = Q * real(ZDrive) #Target reactance

    X_SerCap = Reactance_Load - Xs #Reactance of series capacitor
    SerCap = 1 / (ωDr * X_SerCap) #Series Capacitor value
    CParAct = Q / (ωDr * TargetZ)

    return SerCap,CParAct

end

function ImpMatch_CLoad(TargetZ, ZDrive, Reactance_Load,ωDr)
    println("Load is capacitive")
        EquivSerC = abs.(1 / (Reactance_Load * ωDr))
        println(
            "Load appears to be a capacitor with a value of $(round(EquivSerC*1e6;sigdigits=3))μF and a ESR of $(real(ZDrive))Ω ",
        )
        LSer = abs.(Reactance_Load) / (ωDr)
        matchRatio = real(TargetZ) / real(ZDrive)
        Q = sqrt(matchRatio - 1)
        Xs = Q * real(ZDrive)
        LSer2 = Xs / (ωDr)
        LTee_2 = LSer2 + LSer

        CParAct = findResPair((1 + Q^(-2)) * LSer2, DriveFreq)
        LTee_2_Geom =
            ToroidOptimizer(WireDiam, LTee_2; CuFillFactor = WireFillFac)
        LTee_2_ESR = LTee_2_Geom.DCore.Resistance
        println("Q = $Q")

        return LTee_2,LTee_2_ESR,CParAct

end

"""
This function determines the ideal frequency of a resonant RLC such that it can be matched to a different impedance. It assumes there is already an L, C, and R, so at some frequency it the real part of the admittance will be the matched admittance. This admittance is made fully real with a shunt reactive element

"""
function findFilterFreq(Q,RLoad,L,C)
f = ((Q*RLoad + √(Q^2*RLoad^2+4*L/C)) / (4*π*L) , (Q*RLoad - √(Q^2*RLoad^2+4*L/C)) / (4*π*L) )

return maximum(f)
end

"""
This function makes a filter section that notches out a frequency (NotchFreq), but is apparently real at another (DriveFreq)
    the inputs are:
        * NotchFreq: the freq that will be shunted out
        * DriveFreq: The freq that should remain real
        * RMatch: The resistance the drive coil is tuned to
        * LVal is a kwarg because either the inductor or capacitor value should be pre-defined
The outputs are:
    LVal(Henries)
    CVal(Farads)
    LTune(Henries)
"""
function makeNotchSection(NotchFreq, DriveFreq, RMatch;LVal = 100e-6)
CVal = findResPair(LVal, NotchFreq)
Resid_Reactance = -1*(imag(Par(Z_Cap(CVal,DriveFreq)+Z_Ind(LVal,DriveFreq),RMatch)))
LTune = Resid_Reactance/(2*π*DriveFreq)
    return LVal,CVal,LTune
end

"""

Calculates parallel impedances.
If given two inputs, it assumes both are inputs
If given an array, it takes the parallel combo of it all

Eg.
ZList = [1 + 1*im, 1]
Par(ZList[1], ZList[2])
Par(ZList)
"""
function Par(Z1, Z2)
    return 1 / (1 / Z1 + 1 / Z2)
end
function Par(Z::Array)
    Y = 1 ./Z
    YTotal = sum(Y)
    Zeff = 1/YTotal
    return Zeff
end

"""
Parallel component lists
Input is Array{Tuple{Number,String},1}
e.g.
Freq = 25e3
CompList = [(3, "R"),(30e-6, "L"),(200e-9,"C")]
Par(CompList,Freq)

"""
function Par(C::Array{Tuple{Real,String},1},freq)

    Z = Complex.(zeros(length(C)))
    for i in 1:length(C)
        if (lowercase(C[i][2])=="r") | (lowercase(C[i][2])=="resistor")
            Z[i] = C[i][1]
        elseif (lowercase(C[i][2])=="l") | (lowercase(C[i][2])=="inductor" )
            Z[i] = Z_Ind(C[i][1],freq)
        elseif (lowercase(C[i][2])=="c") | (lowercase(C[i][2])=="capacitor" )
            Z[i] = Z_Cap(C[i][1],freq)
        elseif lowercase(C[i][2])=="z"
            Z[i] = C[i][1]
        else
            error("Unknown component")
        end

    end

    Y = 1 ./Z
    YTotal = sum(Y)
    Zeff = 1/YTotal
    return Zeff
end

"""
Converts from a polar coordinate (Mag, phase) to cartesian
if you enter a tuple or array it assumes first is mag.

"""
function Polar2Cart(Mag, Phase)
    Real = real.(Mag * exp(-im * Phase / 360 * 2 * pi))
    Imag = imag.(Mag * exp(-im * Phase / 360 * 2 * pi))
    return (Real, Imag)
end

function Polar2Cart(MagPhaseTup::Union{Tuple, Array})
    Polar2Cart(MagPhaseTup[1], MagPhaseTup[2])
end

function Z_Cap(C, f)
    return (1 ./ (im * 2 * pi .* f .* C))
end
function Z_Ind(L, f)
    return im * 2 * pi .* f .* L
end

"""
This function takes in some reactance element (E) and a frequency, and finds the pair

"""
function findResPair(E, f)
    ##Finds the resonant pair for a reactive element, E
    ## ωL = 1/(ωC) → L = 1/(ω^2C),and vice versa
    return 1 ./ ((2 * pi * f) .^ 2 * E)
end

function Butterworth_2(Z, f)
    #Sqrt(L/C) = Z ; f = 1/(2π*sqrt(LC))
    #L = Z^2*C
    #f = 1/(2π*Z C) → C = 1/(2*pi*f*Z)
    #Butterworth 2nd order = [LX = sqrt(2)Z, CX = sqrt(2)Z]
    C = 1 / (2 * pi * f) / Z * sqrt(2)
    L = 1 / (2 * pi * f) * Z * sqrt(2)
    return (L, C)
end

function PlotPhasor(
    Impedance::Tuple,
    PrevImpedance;
    Color = nothing,
    NumSteps = 1,
    Admittance = false,
    HeadWidth = 0.01,
)

    ## This function takes in an impedance, and plots it using 'phasor()'
    #optionally, it takes in a previous impedance which can be added in series or parallel
    #It then plots the line using PrevImpedance as a base
    #If Admittance is true, it means the inputs are Admittance not impedance, so parallel and series are switched
    NewImp = 0 + 0im
    dx = 0
    dy = 0
    for I = 1:NumSteps
        if I > 1
            #    Color = "gray"
        end
        if lowercase(Impedance[2]) == "series"
            if !Admittance
                NewImp = PrevImpedance + Impedance[1] / NumSteps
            else
                NewImp = Par(PrevImpedance, Impedance[1] * NumSteps)
            end
        elseif lowercase(Impedance[2]) == "parallel"
            if !Admittance
                NewImp = Par(PrevImpedance, Impedance[1] * NumSteps)
            else
                NewImp = PrevImpedance + Impedance[1] / NumSteps
            end
        else
            error("must be series or parallel")
        end

        RealParts = [real(NewImp), real(PrevImpedance)]
        ImagParts = [imag(NewImp), imag(PrevImpedance)]

        dx = -(RealParts[2] - RealParts[1])
        dy = -(ImagParts[2] - ImagParts[1])
        if Color !== nothing

            PyPlot.plot(RealParts, ImagParts, color = "gray")

            if I == 1
                phasor(PrevImpedance, color = Color)
            end
        else
            PyPlot.plot(RealParts, ImagParts)
            if I == 1
                phasor(PrevImpedance)
            end
        end
        PrevImpedance = NewImp
    end
    # arrow(real(NewImp)-dx,imag(NewImp)-dy,dx,dy,head_width=HeadWidth,fc="gray",ec="k",alpha=1.,width = 0.00001,head_length= 2*HeadWidth)

    xlabel("Real")
    ylabel("Imag.")
    return NewImp
end

"""
This function takes in an array of tuples structures such as:

    [(<Complex impedance #1> , <"series" or "Parallel">),
    (<Complex impedance #2> , <"series" or "Parallel">)
    .
    .
    .
    (<Complex impedance #n> , <"series" or "Parallel">)]

    e.g.

    ZList = [(1+1*im, "series"),
             (1+5*im, "series"),
             (5+1*im, "parallel"),]
    PlotImpedanceTransformList(ZList)
"""
function PlotImpedanceTransformList(ImpList; InitialImp = nothing, ArrHeadWidth =1)
    ColorList = ["r", "g", "b", "magenta", "teal"]
    L = zeros(length(ImpList))
    for I = 1:length(ImpList)
        L[I] = abs(ImpList[I][1])
    end
    LongestVector = maximum(L)
    HeadWidth = LongestVector / ArrHeadWidth
    ##Impedance Plotting
    if InitialImp !== nothing
        PrevImp = InitialImp
        StartVal = 1
    else
        PrevImp = ImpList[1][1]
        StartVal = 2
    end
    for I = StartVal:length(ImpList)
        PrevImp = PlotPhasor(
            ImpList[I],
            PrevImp;
            NumSteps = 500,
            Color = ColorList[I],
            HeadWidth = HeadWidth,
        )
    end
    phasor(PrevImp, color = ColorList[length(ImpList)+1])
    title("Impedance")
    ##Admittance Plotting
    figure()

    L = zeros(length(ImpList))
    for I = 1:length(ImpList)
        L[I] = abs(1 / ImpList[I][1])
    end
    LongestVector = maximum(L)
    HeadWidth = LongestVector / ArrHeadWidth
    if InitialImp !=0 nothing
        PrevImp = 1 / InitialImp
        StartVal = 1
    else
        PrevImp = 1 / ImpList[1][1]
        StartVal = 2
    end
    for I = StartVal:length(ImpList)
        PrevImp = PlotPhasor(
            (1 / ImpList[I][1], ImpList[I][2]),
            PrevImp;
            NumSteps = 500,
            Color = ColorList[I],
            Admittance = true,
            HeadWidth = HeadWidth,
        )
    end
    phasor(PrevImp, color = ColorList[length(ImpList)+1])
    title("Admittance")
end

function CircModel_MatchingTFilt(DriveFreq, VSrc,
    PlotSrcR,
    RDrive,
    NumDriveElements,
    LDrive,
    CDrive,
    SerCap,
    LTee_2,
    LTee_2_ESR,
    LTee_1,
    LTee_1_ESR,
    CParAct,
    LFilt,
    LFilt2_ESR,
    LFilt1_ESR,
    LFiltMatch_C,
    CFilt,
    CFiltMatch_L;
    PlotOn = false

)

    circ = @circuit begin
        j_in = voltagesource()
        rs = resistor(PlotSrcR)
        ESR = resistor(RDrive)
        ESL = inductor(NumDriveElements * LDrive)
        ESC = capacitor(CDrive / NumDriveElements)
        MatchCSer = capacitor(SerCap)
        LTee_2_mod = inductor(LTee_2)
        LTee_2_ESRmod = resistor(LTee_2_ESR)
        LTee_1_mod = inductor(LTee_1)
        LTee_1_ESRmod = resistor(LTee_1_ESR)
        CParAct_mod = capacitor(CParAct)
        LFilt_mod = inductor(LFilt)
        LFiltESR_mod = resistor(LFilt1_ESR)
        LFilt_mod2 = inductor(LFilt)
        LFiltESR_mod2 = resistor(LFilt1_ESR)
        LFiltMatch_C_mod2 = capacitor(LFiltMatch_C)
        CFilt_mod = capacitor(CFilt)
        CFiltMatch_L_mod = inductor(CFiltMatch_L)
        LFilt2_ESR_mod = resistor(LFilt2_ESR)
        LFiltMatch_C_mod = capacitor(LFiltMatch_C)
        i_out = currentprobe()

        j_in[+] ⟷ rs[1] #\longleftrightarrow
        j_in[-] ⟷ gnd
        rs[2] ⟷ LFiltESR_mod[1]
        LFiltESR_mod[2] ⟷ LFilt_mod[1]
        LFilt_mod[2] ⟷ LFiltMatch_C_mod[1]
        LFiltMatch_C_mod[2] ⟷ CFilt_mod[1]
        LFiltMatch_C_mod[2] ⟷ CFiltMatch_L_mod[1]
        CFiltMatch_L_mod[2] ⟷ LFilt2_ESR_mod[1]
        LFilt2_ESR_mod[2] ⟷ gnd
        CFilt_mod[2] ⟷ gnd
        LFiltMatch_C_mod[2] ⟷ LFiltESR_mod2[1]
        LFiltESR_mod2[2] ⟷ LFilt_mod2[1]
        LFilt_mod2[2] ⟷ LFiltMatch_C_mod2[1]
        LFiltMatch_C_mod2[2] ⟷ LTee_1_mod[1]
        LTee_1_mod[2] ⟷ LTee_1_ESRmod[1]
        LTee_1_ESRmod[2] ⟷ CParAct_mod[1]
        CParAct_mod[2] ⟷ gnd
        LTee_1_ESRmod[2] ⟷ LTee_2_mod[1]
        LTee_2_mod[2] ⟷ LTee_2_ESRmod[1]
        LTee_2_ESRmod[2] ⟷ ESC[1]
        ESC[2] ⟷ MatchCSer[1]
        MatchCSer[2] ⟷ ESL[1]
        ESL[2] ⟷ ESR[1]
        ESR[2] ⟷ i_out[+]
        i_out[-] ⟷ gnd
    end
    println("Starting to model")
    fs_model = 1e6
    model = DiscreteModel(circ, 1 / fs_model)
    FreqList = 1000:100:100e3
    NumPeriods = 1000
    df = FreqList[2] - FreqList[1]
    Mags = zeros(length(FreqList), 1)
    WindowHanning(N) =
        0.5 .- 0.5 .* cos.(2 .* pi .* collect(0:(N-1)) ./ (N - 1))
    println("Starting to run sim")
    for F in FreqList
        y = run!(
            model,
            [
                VSrc .* sin(2π * F / fs_model * n) for c = 1:1,
                n = 1:NumPeriods*fs_model/F
            ],
        )
        y = y[:] .* WindowHanning(length(y))[:]
        FFT_y = abs.(fft(y)) / length(y) .* 4

        Mags[Int((F - FreqList[1]) / df)+1] = maximum(FFT_y)
    end

    y2 = run!(
        model,
        [
            VSrc .* sin(2π * DriveFreq / fs_model * n) for c = 1:1,
            n = 1:NumPeriods*fs_model/DriveFreq
        ],
    )
    y2 = y2[:] .* WindowHanning(length(y2))[:]
    FFT_y2 = abs.(fft(y2)) / length(y2) .* 4

    AmpsPerVolt = maximum(FFT_y2)

    if PlotOn
        semilogy(FreqList, Mags[:])
        xlabel("Frequency")
        ylabel("Drive Current")
        semilogy(DriveFreq, AmpsPerVolt,"r*")
    end

    println("Tx current per $VSrc Volts is $(round(AmpsPerVolt,sigdigits=3)) Amps")
    return AmpsPerVolt
end

function findMinima_Bounds(x,y;StartVal=nothing,StopVal=nothing)

    if !(StartVal===nothing)
        StartIndex = minimum(findall(x .>= StartVal))
    else
        StartIndex = 1
    end

    if !(StopVal===nothing)
        StopIndex = maximum(findall(x .<= StopVal))
    else
        StopIndex = length(x)
    end

    xSubset = x[StartIndex:StopIndex]
    ySubset = y[StartIndex:StopIndex]
    FirstDiff = diff(ySubset)
    SecondDiff = diff(FirstDiff)
    Tmp = SecondDiff .< 0
    push!(Tmp,1)
    FirstDiff[Tmp] .= 1000 #Only look at values where the second derivative is positive

    minIndex = findfirst(isequal(minimum(abs.(FirstDiff))),abs.(FirstDiff))

    Min_x = xSubset[minIndex]
    Min_y = ySubset[minIndex]
    return Min_x, Min_y

end


text:
### A Pluto.jl notebook ###
# v0.17.4

using Markdown
using InteractiveUtils

# ╔═╡ 59ed6c40-8025-11ec-0ace-d5b826d81255
begin
	bigbreak = HTML("<br>"^4);
	using PlutoUI
	using Test
end

# ╔═╡ c3b89eff-c3e6-47d0-acc8-661dde65305c
md"""
# Data structures
Once we start working with many pieces of data at once, it will be convenient for us to store data in structures like arrays or dictionaries (rather than just relying on variables).

Types of data structures covered:
1. Tuples and NamedTuples
2. Dictionaries
3. Arrays

As an overview, tuples and arrays are both ordered sequences of elements (so we can index into them). Dictionaries and arrays are both mutable. We'll explain this more below!

"""

# ╔═╡ a58c0b6f-c5f4-43e9-87ad-4eca62157ce4
md"""
## Tuples
We can create a tuple by enclosing an ordered collection of elements in ( ).

Syntax:
```julia
(item1, item2, ...)
```
"""

# ╔═╡ 30233243-ed9d-4bfd-a9af-fcde96677998
myfavoriteanimals = ("penguins", "cats", "sugargliders")

# ╔═╡ 15717659-eb1a-40d4-abfc-7d680f4722c8
md"We can index into this tuple,"

# ╔═╡ 9cd7ec78-75dc-4cfc-a4ef-2fd8e9a8d40c
myfavoriteanimals[1]

# ╔═╡ dae03614-aa84-483a-8e54-5f493a5a1b82
md"but since tuples are immutable, we can't update it"

# ╔═╡ 6d1bb7e6-6cff-4cee-873d-ccc82f4f0be8
myfavoriteanimals[1] = "otters"

# ╔═╡ d85c4c9a-3470-4b5a-95af-f2e9c2b04c00
md"""## Now in 1.6: NamedTuples
As you might guess, NamedTuples are just like Tuples except that each element additionally has a name! They have a special syntax using = inside a tuple:
```julia
(name1 = item1, name2 = item2, ...)
```
_Special case, for just one entry, use `(; name1 = item1)`_
"""

# ╔═╡ a8e3433a-a9d3-4081-a8fe-5fa8732c8d71
yourfavoriteanimals = (bird = "penguins", mammal = "cats", marsupial = "sugargliders")

# ╔═╡ 9bde2570-6ce5-4ad7-8b0b-dd7889f8bea1
md"
Like regular Tuples, NamedTuples are ordered, so we can retrieve their elements via indexing:"

# ╔═╡ f425304c-29c1-40eb-b606-9dc39a62e252
yourfavoriteanimals[1]

# ╔═╡ 9db7047d-d900-403c-9a10-85328d5faf44
md"They also add the special ability to access values by their name:"

# ╔═╡ 560ddf25-33e0-4075-b1f4-475f308c99a7
yourfavoriteanimals.bird

# ╔═╡ ff1774ca-7a03-474a-80c7-ff2bd9f38e61
bigbreak

# ╔═╡ e13cae8d-e8ad-4ed8-a844-e1cedb452493
md""" ## Dictionaries
If we have sets of data related to one another, we may choose to store that data in a dictionary. We can create a dictionary using the Dict() function, which we can initialize as an empty dictionary or one storing key, value pairs.

Syntax:
```julia
Dict(key1 => value1, key2 => value2, ...)
```
A good example is a contacts list, where we associate names with phone numbers.
"""

# ╔═╡ 079ece15-7872-4e24-b762-418152f232e6
myphonebook = Dict("Jenny" => "867-5309", "Ghostbusters" => "555-2368")

# ╔═╡ 37887043-b092-4d0f-bbc7-4288c4c761cd
md"""
In this example, each name and number is a "key" and "value" pair. We can grab Jenny's number (a value) using the associated key
"""

# ╔═╡ 07fb97a3-60f0-40e5-befe-5d27a0b5d597
myphonebook["Jenny"]

# ╔═╡ 5efd2a9e-7bc9-424e-812d-2fd554677b55
md"We can add another entry to this dictionary as follows"

# ╔═╡ 6647ab87-0ff0-4d61-ad10-897942b461dc
myphonebook["Kramer"] = "555-FILK"

# ╔═╡ 0b438e45-fa85-435b-9319-741379f8bf2d
myphonebook

# ╔═╡ a9ef31b2-9c59-4ed4-ba78-173a383211e2
md"We can delete Kramer from our contact list - and simultaneously grab his number - by using pop!"

# ╔═╡ d869d9fd-703b-475b-b9c2-337cbc1b967a
pop!(myphonebook, "Kramer")

# ╔═╡ 09eb599f-f454-4167-8b63-002e6b3b1ef9
myphonebook

# ╔═╡ 5c30cf8d-1295-4e42-ba93-2b49b005455c
md"Unlike tuples and arrays, dictionaries are not ordered. So, we can't index into them."

# ╔═╡ d96442d4-eaa9-4339-b55f-e48e39169e48
myphonebook[1]

# ╔═╡ 8e1fac94-4a44-499b-a49b-e4761549b166
md"
In the example above, julia thinks you're trying to access a value associated with the key 1."

# ╔═╡ 0536d875-6f77-4834-9fa4-a0a31f7fa09a
bigbreak

# ╔═╡ d54f80fa-750b-494b-be58-933d375fb8c2
md"""
## Arrays
Unlike tuples, arrays are mutable. Unlike dictionaries, arrays contain ordered collections.
We can create an array by enclosing this collection in `[ ]`.

Syntax:
```julia
[item1, item2, ...]
```
For example, we might create an array to keep track of my friends

"""

# ╔═╡ 8bd0408e-7b00-419c-930c-8e9361ad1574
myfriends = ["Ted", "Robyn", "Barney", "Lily", "Marshall"]

# ╔═╡ 3b2e4202-19be-4de5-92ae-b3911af42f7f
typeof(myfriends)

# ╔═╡ 3b3cab5a-c5e2-4587-9aae-a97cbc15f238
md"The 1 in `Array{String,1}` means this is a one dimensional vector. An `Array{String,2}` would be a 2d matrix, etc. The String is the type of each element."

# ╔═╡ c94ca702-fe13-468f-938a-be9b16a879e2
md"or to store a sequence of numbers"

# ╔═╡ b577d259-b564-47c8-b732-a0af01f30014
fibonacci = [1, 1, 2, 3, 5, 8, 13]

# ╔═╡ 24dea8be-3930-4049-9a44-b8151e164d91
typeof(fibonacci)

# ╔═╡ 7701def7-26b9-4f53-90a9-3eec7be70a16
eltype(fibonacci)

# ╔═╡ d1571780-39a7-4d8d-9ee0-b9770e235bdd
md"we can also mix types"

# ╔═╡ 98457248-d957-4376-8086-6b5866a743c9
mixture = [1, 1, 2, 3, "Ted", "Robyn"]

# ╔═╡ 0553623d-5416-4cf4-8cd8-a836dbe5c41f
typeof(mixture)

# ╔═╡ cf556b7b-1d94-43d7-99f8-0a134bd9865b
md"Once we have an array, we can grab individual pieces of data from inside that array by indexing into the array. For example, if we want the third friend listed in myfriends, we write"

# ╔═╡ 46bdba66-e0fa-4631-b0a1-9d532aa25e89
myfriends[1]

# ╔═╡ 6debf09f-ff30-4e30-93b3-8872c258908c
md"
We can use indexing to edit an existing element of an array"

# ╔═╡ 28821165-3106-4983-be9f-3155cbc4ba53
myfriends[3] = "Baby Bop"

# ╔═╡ b221fbd3-c519-4c9c-8c98-bcee6fae796d
bigbreak

# ╔═╡ 3b5b7dcf-6a56-49cb-9dd8-1b2aba81ef87
md"""
### Yes, Julia is 1-based indexing, not 0-based like Python.
_(I share your pain)._

"""

# ╔═╡ 8f398c3a-9def-4f41-a171-45d38527c83c
bigbreak

# ╔═╡ a39e4255-43ee-4ed3-8282-fd25ca0b8d53
md"""
Let's continue anyway.

We can also edit the array by using the push! and pop! functions. push! adds an element to the end of an array and pop! removes the last element of an array.

We can add another number to our fibonnaci sequence
"""

# ╔═╡ 214cde58-aa66-4499-9a5e-5a6a36c22115
push!(fibonacci, 21)

# ╔═╡ 7aa62fe2-9c5e-4640-8960-21a9084e436f
md"and then remove it"

# ╔═╡ 687f65ad-2320-4cab-b576-0b94de505c10
pop!(fibonacci)

# ╔═╡ 55e690c8-14c2-46f8-bd3c-f64e34a8d74e
fibonacci

# ╔═╡ 58fbd165-41f3-4cac-aaa3-e498e4621d1d
md"So far I've given examples of only 1D arrays of scalars, but arrays can have an arbitrary number of dimensions and can also store other arrays.

For example, the following are arrays of arrays:"

# ╔═╡ a6fc5004-5ef2-409c-9b6a-5e7e7efb25ad
favorites = [["koobideh", "chocolate", "eggs"],["penguins", "cats", "sugargliders"]]

# ╔═╡ c4d8af83-f7f9-47a0-9023-fe46e6bc1109
favorites[2][1]

# ╔═╡ 6229a275-b52d-4c81-a05b-f22d15b0cf32
numbers = [[1, 2, 3], [4, 5], [6, 7, 8, 9]]

# ╔═╡ 53cd90d5-8a70-43a6-b72e-d61185c0b7aa
typeof(numbers)

# ╔═╡ 4fdc2f91-dc97-4c58-9cc1-486d30ad7890
bigbreak

# ╔═╡ 08924e9f-27b8-444a-b9af-9fa34cf1d61b
md"Below are examples of 2D and 3D arrays populated with random values."

# ╔═╡ daa23583-9b16-4a76-bde8-9324f3dc1908
rand(4,3)

# ╔═╡ 12e389d6-9ac2-49b5-a678-cf71787541eb
rand(4,3,2)

# ╔═╡ e791a417-a00a-479e-b433-3773dd74a48e
bigbreak

# ╔═╡ b18f9af8-130d-4bac-ba8f-230e89d229b7
X = [1 2; 2 3 ]

# ╔═╡ 26de3d3b-9229-4c34-a7d1-e0c89ef6f843
Z = [i+j for i in 1:10, j in 1:10]

# ╔═╡ 30933711-723e-4b8f-8be9-3928635e068e
Z[3:5, 2:2:6]

# ╔═╡ eb9a5d9f-9f39-401c-aa35-102414ed99de
bigbreak

# ╔═╡ 3280db88-0d68-4426-9f89-0fd3914e6ffa
md"Be careful when you want to copy arrays!"

# ╔═╡ 4fa1d5b9-36a0-491e-a034-c887ee7338d8
fibonacci

# ╔═╡ 68ce440d-53f8-4312-8911-a7348e602a1f
somenumbers = copy(fibonacci)

# ╔═╡ 1be992a6-fc03-437f-9dc1-6072f0b08327
somenumbers[1] = 404

# ╔═╡ 36086001-a2bd-4648-b2b3-3d526f0dcf71
fibonacci

# ╔═╡ 5ece33d6-4622-4385-b84c-639bd7571915
md"Editing somenumbers caused fibonacci to get updated as well!

In the above example, we didn't actually make a copy of fibonacci. **We just created a new way to access the entries in the array bound to fibonacci.**

If we'd like to make a copy of the array bound to fibonacci, we can use the copy function."

# ╔═╡ ce690520-a919-4e3d-a201-ec2baedf913d
let
	a = [[1,2],1]

	b = deepcopy(a)
	b[1][1] = 10
	b[2] = 2
	a
end

# ╔═╡ 28263f4f-37b0-43fd-9eab-6169d82b0581
bigbreak

# ╔═╡ 421ab876-efa5-4a7d-9d33-5b9eca5f528e
begin
	# let's restore the original array
	fibonacci[1] = 1
	fibonacci
end

# ╔═╡ 142b9281-196f-4d3e-ae7d-fbc1011255f6
somemorenumbers = copy(fibonacci)

# ╔═╡ 91e7e129-c7b5-4e1f-9025-bdac4d569b46
somemorenumbers[1] = 404

# ╔═╡ 0e20de93-a6f1-415c-a08e-53de67065bc4
fibonacci

# ╔═╡ be7bad99-185a-43ed-9c8d-33e5abcdbbd0

# ╔═╡ 40b0919d-64a5-4a89-95fb-5266c2c07260
md"In this last example, fibonacci was not updated. Therefore we see that the arrays bound to somemorenumbers and fibonacci are distinct."

# ╔═╡ ec0fa200-0c79-467a-af5c-97e79f37253d
let
	function gradientdecent(f, x0)

	end

	x0 = [1., 2.]
	gradientdecent( nothing, x0)
end

# ╔═╡ 55633c9d-2c1b-40c5-9475-48f18b1f8d77
bigbreak

# ╔═╡ 69e96bb3-fc25-4f44-89cd-8f2aeb7065fe
md"Variables are just associated to values!"

# ╔═╡ 9d4aeb68-6649-42f6-8a93-997a136499d4
begin
	A = [1,2,3]

	A_orig = Ref(A)  # points to the array [1,2,3]

	A .= (2 * A)

end

# ╔═╡ 2a3d5a98-751c-48c9-b102-89bc40e74251
typeof(2*A)

# ╔═╡ 457bc2c7-9aae-4b28-b216-ef9454227651
A

# ╔═╡ f35046c4-4a6e-42e5-b691-78b6ef3b0d27
md"The original array (might) still exists and we wasted memory!"

# ╔═╡ 8af669a5-4f1f-4f0e-9de7-030e3196f96b
A_orig[]

# ╔═╡ 32035322-60ec-4ed7-a340-84c6457a2464
md"""
#### Why did we not overwrite the original array?
Imagine, in the first step we did
```julia
A = [1, 2, 3]
```

$(LocalResource("images/var_1.svg", :width => 800))
After the command
```julia
A = 2*A
```
the variable is bound to a new value. **The old value remains unchanged!**

$(LocalResource("images/var_2.svg", :width => 800))

"""

# ╔═╡ 2c04e557-4afa-4acc-ace1-45f994aa3ed3
bigbreak

# ╔═╡ 83c37801-d248-442a-a360-829f740f5cd9
md"""
## Exercises
3.1
Create an array, `x` with the number 1 to 20 and add `+ 10` to all odd entries.

"""

# ╔═╡ 4170966b-cbb4-4b7b-b6e0-528f63a1c019
begin
	x = []
	# ...
end

# ╔═╡ f720b6ad-eea0-4c07-b528-822b80305000
md"""_The current solution is $( (x == [i + (i%2==1)*10 for i in 1:20]) ? "correct" : "false")._"""

# ╔═╡ 7831cafd-8ac8-41be-ad0d-a5abeeb66276
bigbreak

# ╔═╡ fc60372a-f853-4005-bc40-64e9d0615a32
md"""
3.2

Try to add `"Emergency"` as key to `ourphonebook` with the value `string(911)` with the following code
```julia
ourphonebook["Emergency"] = 112
```
Why doesn't this work?
"""

# ╔═╡ 89b74238-dd21-4063-bdb8-d64adcc2184c
ourphonebook = Dict("Jenny" => "867-5309", "Ghostbusters" => "555-2368")

# ╔═╡ b8c1cdf5-550e-49be-83b0-2a01fc583886
# enter you code here

# ╔═╡ 1b7789aa-1969-486c-a733-6e3d10d135bc

# ╔═╡ f42fa416-47d1-49b6-9d4b-365740908e03
md"How can you change the definition of `ourphonebook` such that you can add and integer value?"

# ╔═╡ 4f191e04-1972-47af-b1d9-2de4857833bb
bigbreak

# ╔═╡ 342337cc-6202-45f7-bf5d-6688b770ac75
md"""
!!! hint "Hint"
      Check out `typeof()` for these dictionaries. You can use `Dict{K,V}` instead of the generic `Dict` to create a dictionary with key type `T` and value type `V`. What should `V` be?
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[compat]
PlutoUI = "~0.7.32"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "0b5cfbb704034b5b4c1869e36634438a047df065"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.1"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "ae6145ca68947569058866e443df69587acc1806"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.32"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─59ed6c40-8025-11ec-0ace-d5b826d81255
# ╟─c3b89eff-c3e6-47d0-acc8-661dde65305c
# ╟─a58c0b6f-c5f4-43e9-87ad-4eca62157ce4
# ╠═30233243-ed9d-4bfd-a9af-fcde96677998
# ╟─15717659-eb1a-40d4-abfc-7d680f4722c8
# ╠═9cd7ec78-75dc-4cfc-a4ef-2fd8e9a8d40c
# ╟─dae03614-aa84-483a-8e54-5f493a5a1b82
# ╠═6d1bb7e6-6cff-4cee-873d-ccc82f4f0be8
# ╟─d85c4c9a-3470-4b5a-95af-f2e9c2b04c00
# ╠═a8e3433a-a9d3-4081-a8fe-5fa8732c8d71
# ╟─9bde2570-6ce5-4ad7-8b0b-dd7889f8bea1
# ╠═f425304c-29c1-40eb-b606-9dc39a62e252
# ╟─9db7047d-d900-403c-9a10-85328d5faf44
# ╠═560ddf25-33e0-4075-b1f4-475f308c99a7
# ╟─ff1774ca-7a03-474a-80c7-ff2bd9f38e61
# ╟─e13cae8d-e8ad-4ed8-a844-e1cedb452493
# ╠═079ece15-7872-4e24-b762-418152f232e6
# ╟─37887043-b092-4d0f-bbc7-4288c4c761cd
# ╠═07fb97a3-60f0-40e5-befe-5d27a0b5d597
# ╟─5efd2a9e-7bc9-424e-812d-2fd554677b55
# ╠═6647ab87-0ff0-4d61-ad10-897942b461dc
# ╠═0b438e45-fa85-435b-9319-741379f8bf2d
# ╟─a9ef31b2-9c59-4ed4-ba78-173a383211e2
# ╠═d869d9fd-703b-475b-b9c2-337cbc1b967a
# ╠═09eb599f-f454-4167-8b63-002e6b3b1ef9
# ╟─5c30cf8d-1295-4e42-ba93-2b49b005455c
# ╠═d96442d4-eaa9-4339-b55f-e48e39169e48
# ╟─8e1fac94-4a44-499b-a49b-e4761549b166
# ╟─0536d875-6f77-4834-9fa4-a0a31f7fa09a
# ╟─d54f80fa-750b-494b-be58-933d375fb8c2
# ╠═8bd0408e-7b00-419c-930c-8e9361ad1574
# ╟─3b2e4202-19be-4de5-92ae-b3911af42f7f
# ╟─3b3cab5a-c5e2-4587-9aae-a97cbc15f238
# ╟─c94ca702-fe13-468f-938a-be9b16a879e2
# ╠═b577d259-b564-47c8-b732-a0af01f30014
# ╠═24dea8be-3930-4049-9a44-b8151e164d91
# ╠═7701def7-26b9-4f53-90a9-3eec7be70a16
# ╟─d1571780-39a7-4d8d-9ee0-b9770e235bdd
# ╠═98457248-d957-4376-8086-6b5866a743c9
# ╠═0553623d-5416-4cf4-8cd8-a836dbe5c41f
# ╟─cf556b7b-1d94-43d7-99f8-0a134bd9865b
# ╠═46bdba66-e0fa-4631-b0a1-9d532aa25e89
# ╟─6debf09f-ff30-4e30-93b3-8872c258908c
# ╠═28821165-3106-4983-be9f-3155cbc4ba53
# ╟─b221fbd3-c519-4c9c-8c98-bcee6fae796d
# ╟─3b5b7dcf-6a56-49cb-9dd8-1b2aba81ef87
# ╟─8f398c3a-9def-4f41-a171-45d38527c83c
# ╟─a39e4255-43ee-4ed3-8282-fd25ca0b8d53
# ╠═214cde58-aa66-4499-9a5e-5a6a36c22115
# ╟─7aa62fe2-9c5e-4640-8960-21a9084e436f
# ╠═687f65ad-2320-4cab-b576-0b94de505c10
# ╠═55e690c8-14c2-46f8-bd3c-f64e34a8d74e
# ╟─58fbd165-41f3-4cac-aaa3-e498e4621d1d
# ╠═a6fc5004-5ef2-409c-9b6a-5e7e7efb25ad
# ╠═c4d8af83-f7f9-47a0-9023-fe46e6bc1109
# ╠═6229a275-b52d-4c81-a05b-f22d15b0cf32
# ╠═53cd90d5-8a70-43a6-b72e-d61185c0b7aa
# ╟─4fdc2f91-dc97-4c58-9cc1-486d30ad7890
# ╟─08924e9f-27b8-444a-b9af-9fa34cf1d61b
# ╠═daa23583-9b16-4a76-bde8-9324f3dc1908
# ╠═12e389d6-9ac2-49b5-a678-cf71787541eb
# ╟─e791a417-a00a-479e-b433-3773dd74a48e
# ╠═b18f9af8-130d-4bac-ba8f-230e89d229b7
# ╠═26de3d3b-9229-4c34-a7d1-e0c89ef6f843
# ╠═30933711-723e-4b8f-8be9-3928635e068e
# ╟─eb9a5d9f-9f39-401c-aa35-102414ed99de
# ╟─3280db88-0d68-4426-9f89-0fd3914e6ffa
# ╠═4fa1d5b9-36a0-491e-a034-c887ee7338d8
# ╠═68ce440d-53f8-4312-8911-a7348e602a1f
# ╠═1be992a6-fc03-437f-9dc1-6072f0b08327
# ╠═36086001-a2bd-4648-b2b3-3d526f0dcf71
# ╟─5ece33d6-4622-4385-b84c-639bd7571915
# ╠═ce690520-a919-4e3d-a201-ec2baedf913d
# ╟─28263f4f-37b0-43fd-9eab-6169d82b0581
# ╠═421ab876-efa5-4a7d-9d33-5b9eca5f528e
# ╠═142b9281-196f-4d3e-ae7d-fbc1011255f6
# ╠═91e7e129-c7b5-4e1f-9025-bdac4d569b46
# ╠═0e20de93-a6f1-415c-a08e-53de67065bc4
# ╠═be7bad99-185a-43ed-9c8d-33e5abcdbbd0
# ╟─40b0919d-64a5-4a89-95fb-5266c2c07260
# ╠═ec0fa200-0c79-467a-af5c-97e79f37253d
# ╟─55633c9d-2c1b-40c5-9475-48f18b1f8d77
# ╟─69e96bb3-fc25-4f44-89cd-8f2aeb7065fe
# ╠═9d4aeb68-6649-42f6-8a93-997a136499d4
# ╠═2a3d5a98-751c-48c9-b102-89bc40e74251
# ╠═457bc2c7-9aae-4b28-b216-ef9454227651
# ╟─f35046c4-4a6e-42e5-b691-78b6ef3b0d27
# ╠═8af669a5-4f1f-4f0e-9de7-030e3196f96b
# ╟─32035322-60ec-4ed7-a340-84c6457a2464
# ╟─2c04e557-4afa-4acc-ace1-45f994aa3ed3
# ╟─83c37801-d248-442a-a360-829f740f5cd9
# ╠═4170966b-cbb4-4b7b-b6e0-528f63a1c019
# ╟─f720b6ad-eea0-4c07-b528-822b80305000
# ╟─7831cafd-8ac8-41be-ad0d-a5abeeb66276
# ╟─fc60372a-f853-4005-bc40-64e9d0615a32
# ╠═89b74238-dd21-4063-bdb8-d64adcc2184c
# ╠═b8c1cdf5-550e-49be-83b0-2a01fc583886
# ╟─1b7789aa-1969-486c-a733-6e3d10d135bc
# ╟─f42fa416-47d1-49b6-9d4b-365740908e03
# ╟─4f191e04-1972-47af-b1d9-2de4857833bb
# ╟─342337cc-6202-45f7-bf5d-6688b770ac75
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: slides/2 - Data structures.jl

max_stars_repo_name: SteffenPL/Julia-for-mathematicians

max_issues_repo_path: slides/2 - Data structures.jl

max_issues_repo_name: SteffenPL/Julia-for-mathematicians

max_forks_repo_path: slides/2 - Data structures.jl

max_forks_repo_name: SteffenPL/Julia-for-mathematicians


text:
module Fields

export classfields, FieldDescriptor

using JavaCall.JNI
using JavaCall.Conversions
using JavaCall.Core
using JavaCall.Reflection: Classes, Modifiers

#=
Struct to hold information about java fields
that can be accessed by the user.
This is not a replecament to java.lang.reflect.Field
as it only holds essential information for
code generation

name:       Name of the field

type:       Descriptor for the type of the field

modifiers:  Information about the field modifiers
=#
struct FieldDescriptor
    name::String
    type::Classes.ClassDescriptor
    modifiers::Modifiers.ModifiersDescriptor
end

function Base.show(io::IO, f::FieldDescriptor)
    print(
        io,
        "FieldDescriptor{name: ", f.name,
        ", type: ", f.type,
        ", modifiers: ", f.modifiers,
        "}"
    )
end

function Base.:(==)(x::FieldDescriptor, y::FieldDescriptor)
    x.name == y.name &&
    x.type == y.type &&
    x.modifiers == y.modifiers
end

function descriptorfromfield(field::jobject)
    name = callinstancemethod(field, :getName, Symbol("java.lang.String"), [])
    type = callinstancemethod(field, :getType, Symbol("java.lang.Class"), [])
    FieldDescriptor(
        convert_to_string(String, name),
        Classes.descriptorfromclass(type),
        Modifiers.methodmodifiers(field))
end

classfields(classname::Symbol) = classfields(Classes.findclass(classname))

function classfields(classdescriptor::Classes.ClassDescriptor)
    array = convert_to_vector(Vector{jobject}, callinstancemethod(
        classdescriptor.jniclass,
        :getFields,
        Vector{Symbol("java.lang.reflect.Field")},
        []))
    map(descriptorfromfield, array)
end

end


text:
module SharedDicts

using Base: KeyIterator, ValueIterator

import Base: getindex, haskey, in, start, done, next, isempty, length
export SharedDict, share, getindex, haskey, in, start, done, next, isempty, length

include("shareddict.jl")

end # module


text:
# ------------------------------------------------------------
# Mimic inheritance of common fields
# ------------------------------------------------------------

function Base.getproperty(x::AbstractStructFRSpace, name::Symbol)
    options =
        union(fieldnames(PSpace1D), fieldnames(PSpace2D), fieldnames(KitBase.CSpace2D))

    if name in options
        return getfield(x.base, name)
    else
        return getfield(x, name)
    end
end

function Base.getproperty(x::AbstractUnstructFRSpace, name::Symbol)
    if name in fieldnames(UnstructPSpace)
        return getfield(x.base, name)
    else
        return getfield(x, name)
    end
end

function Base.propertynames(x::AbstractStructFRSpace, private::Bool = false)
    public = fieldnames(typeof(x))
    true ?
    (
        (
            public ∪
            union(fieldnames(PSpace1D), fieldnames(PSpace2D), fieldnames(KitBase.CSpace2D))
        )...,
    ) : public
end

function Base.propertynames(x::AbstractUnstructFRSpace, private::Bool = false)
    public = fieldnames(typeof(x))
    true ? ((public ∪ fieldnames(UnstructPSpace))...,) : public
end

# ------------------------------------------------------------
# Accuracy analysis
# ------------------------------------------------------------

function L1_error(u::T, ue::T, Δx) where {T<:AbstractArray}
    return sum(abs.(u .- ue) .* Δx)
end

function L2_error(u::T, ue::T, Δx) where {T<:AbstractArray}
    return sqrt(sum((abs.(u .- ue) .* Δx) .^ 2))
end

function L∞_error(u::T, ue::T, Δx) where {T<:AbstractArray}
    return maximum(abs.(u .- ue) .* Δx)
end


text:
using SumOfSquares
using DynamicPolynomials
using MosekTools

# Create symbolic variables (not JuMP decision variables)
@polyvar x1 x2

# Create a Sum of Squares JuMP model with the Mosek solver
model = SOSModel(Mosek.Optimizer)

# Create a JuMP decision variable for the lower bound
@variable(model, γ)

# f(x) is the Goldstein-Price function
f1 = x1+x2+1
f2 = 19-14*x1+3*x1^2-14*x2+6*x1*x2+3*x2^2
f3 = 2*x1-3*x2
f4 = 18-32*x1+12*x1^2+48*x2-36*x1*x2+27*x2^2

f = (1+f1^2*f2)*(30+f3^2*f4)

# Constraints f(x) - γ to be sum of squares
@constraint(model, f >= γ)

@objective(model, Max, γ)

optimize!(model)

# The lower bound found is 3
println(objective_value(model))


text:
################################################################################
#
#  Functions for completion map
#
################################################################################

function image(f::CompletionMap, a::nf_elem)
  Qx = parent(parent(a).pol)
  return evaluate(Qx(a), f.prim_img)
end

function preimage(f::CompletionMap{LocalField{qadic, EisensteinLocalField}, LocalFieldElem{qadic, EisensteinLocalField}}, a::LocalFieldElem{qadic, EisensteinLocalField})
  Kp = codomain(f)
  @assert Kp === parent(a)
  Qq = base_field(Kp)
  Qpx = parent(defining_polynomial(Qq))
  coeffs = Vector{nf_elem}()
  for i = 0:degree(a.data)
    as_pol = Qpx(coeff(a.data, i))
    as_fmpq_poly = map_coefficients(lift, as_pol)
    push!(coeffs, evaluate(as_fmpq_poly, f.inv_img[1]))
  end
  K = domain(f)
  Kx = PolynomialRing(K, "x")[1]
  r = Kx(coeffs)
  return evaluate(r, f.inv_img[2])
end

function preimage(f::CompletionMap{LocalField{padic, EisensteinLocalField}, LocalFieldElem{padic, EisensteinLocalField}}, a::LocalFieldElem{padic, EisensteinLocalField})
  @assert codomain(f) === parent(a)
  return evaluate(map_coefficients(lift, a.data), f.inv_img[2])
end

function preimage(f::CompletionMap{FlintQadicField, qadic}, a::qadic)
  Kp = codomain(f)
  @assert Kp == parent(a)
  Qpx = parent(defining_polynomial(Kp))
  as_pol = Qpx(a)
  as_fmpq_poly = map_coefficients(lift, as_pol)
  return evaluate(as_fmpq_poly, f.inv_img[1])
end

function prime(f::CompletionMap)
  return f.P
end

################################################################################
#
#  Lifting
#
################################################################################

function _lift(a::nf_elem, f::fmpz_poly, prec::Int, P::NfOrdIdl)
  i = prec*ramification_index(P)
  chain = [i]
  while i > 2
    i = div(i+1, 2)
    push!(chain, i)
  end
  push!(chain, 2)
  der = derivative(f)
  for i in length(chain):-1:1
    ex, r = divrem(chain[i], ramification_index(P))
    if r > 0
      ex += 1
    end
    minP2i = minimum(P)^ex
    a = a - f(a)//der(a)
    a = mod(a, minP2i)
  end
  return a
end

function _increase_precision(a::nf_elem, f::fmpz_poly, prec::Int, new_prec::Int, P::NfOrdIdl)
  i = new_prec
  chain = [new_prec]
  while i > prec
    i = div(i+1, 2)
    pushfirst!(chain, max(i, prec))
  end
  der = derivative(f)
  for i = 2:length(chain)
    ex, r = divrem(2^chain[i], ramification_index(P))
    if r > 0
      ex += 1
    end
    minP2i = minimum(P)^ex
    a = a - f(a)//der(a)
    a = mod(a, minP2i)
  end
  return a
end

################################################################################
#
#  Generic completion
#
################################################################################

function generic_completion(K::AnticNumberField, P::NfOrdIdl, precision::Int = 64)
  OK = order(P)
  @assert isprime(P)
  @assert nf(OK) == K
  f = degree(P)
  e = ramification_index(P)
  prec_padics = div(precision, e)
  Qp = PadicField(minimum(P), prec_padics)
  Zp = maximal_order(Qp)
  Qq, gQq = QadicField(minimum(P), f, prec_padics)
  Qqx, gQqx = PolynomialRing(Qq, "x")
  q, mq = ResidueField(Qq)
  F, mF = ResidueFieldSmall(OK, P)
  mp = find_morphism(q, F)
  g = gen(q)
  gq_in_K = (mF\(mp(g))).elem_in_nf
  Zx = PolynomialRing(FlintZZ, "x")[1]
  pol_gq = map_coefficients(lift,  defining_polynomial(Qq))
  gq_in_K = _lift(gq_in_K, pol_gq, precision, P)
  #@assert mF(OK(gq_in_K)) == mp(g)
  u = uniformizer(P).elem_in_nf
  pows_gq = powers(gq_in_K, f-1)
  els = Vector{nf_elem}()
  el = one(K)
  for i = 1:e
    for j = 1:f
      push!(els, el*pows_gq[j])
    end
    mul!(el, el, u)
  end
  append!(els,  map(elem_in_nf, basis(P^precision, copy = false)))
  MK = basis_matrix(els, FakeFmpqMat)
  bK = basis_matrix(nf_elem[u^e, gen(K)], FakeFmpqMat)
  d = lcm(denominator(MK, copy = false), denominator(bK, copy = false))
  if d != denominator(MK, copy = false)
    mul!(MK.num, mK.num, divexact(d, denominator(MK, copy = false)))
  end
  if d != denominator(bK, copy = false)
    mul!(bK.num, bK.num, divexact(d, denominator(bK, copy = false)))
  end
  MZp = map_entries(Zp, MK.num)
  bZp = map_entries(Zp, bK.num)
  fl, xZp = can_solve_with_solution(MZp, bZp, side = :left)
  @assert fl
  coeffs_eisenstein = Vector{qadic}(undef, e+1)
  for i = 1:e
    coeff = Qq()
    for j = 0:f-1
      coeff += (gQq^j)*xZp[1, j+1+(i-1)*f].x
    end
    coeffs_eisenstein[i] = -coeff
  end
  coeffs_eisenstein[e+1] = one(Qq)
  pol_gen = Qqx(coeffs_eisenstein)
  Kp, gKp = eisenstein_extension(pol_gen, "a")
  img_prim_elem = Vector{qadic}(undef, e)
  for i = 1:e
    coeff = Qq()
    for j = 0:f-1
      coeff += (gQq^j)*xZp[2, j+1+(i-1)*f].x
    end
    img_prim_elem[i] = coeff
  end
  img = Kp(Qqx(img_prim_elem))
  completion_map = CompletionMap(K, Kp, img, (gq_in_K, u), precision)
  completion_map.P = P
  return Kp, completion_map
end

function setprecision!(f::CompletionMap{LocalField{qadic, EisensteinLocalField}, LocalFieldElem{qadic, EisensteinLocalField}}, new_prec::Int)
  if new_prec < f.precision
    K = domain(f)
    setprecision!(K, new_prec)
    setprecision!(base_field(K), new_prec)
    setprecision!(f.prim_img, new_prec)
  else
    #I need to increase the precision of the data
    P = prime(f)
    f = inertia_degree(P)
    e = ramification_index(P)
    Kp = codomain(f)
    gq, u = f.inv_img
    Zx = PolynomialRing(FlintZZ, "x")[1]
    pol_gq = lift(Zx, defining_polynomial(q))
    gq = _increase_precision(gq, pol_gq, f.precision, new_prec, P)
    f.inv_img[1] = gq
    Kp = codomain(f)
    Qq = base_field(Kp)
    setprecision!(Qq, new_prec)
    Qqx = PolynomialRing(Qq, "x")[1]
    Qp = PadicField(prime(Kp), new_prec)
    Qpx = PolynomialRing(Qp, "x")
    ex, r = divrem(precision, e)
    if r > 0
      ex += 1
    end
    modulo = minimum(P)^ex
    pows_gq = powers_mod(gq_in_K, f, modulo)
    els = Vector{nf_elem}()
    el = one(K)
    for i = 1:e
      for j = 1:f
        push!(els, mod(el*pows_gq[j], modulo))
      end
      mul!(el, el, u)
    end
    MK = basis_matrix(els)
    bK = basis_matrix(nf_elem[_powermod(u, e, modulo), gen(K)])
    MQp = map_entries(Qp, MK)
    bQp = map_entries(Qp, bK)
    fl, xQp = can_solve_with_solution(MQp, bQp, side = :left)
    @assert fl
    coeffs_eisenstein = Vector{qadic}(undef, e+1)
    for i = 1:e
      coeff = Qq()
      for j = 0:f-1
        coeff -= (gQq^j)*xQp[1, j+1+(i-1)*f]
      end
      coeffs_eisenstein[i] = coeff
    end
    coeffs_eisenstein[e+1] = one(Qq)
    pol_gen = Qqx(coeffs_eisenstein)
    Kp.defining_polynomial = pol_gen
    img_prim_elem = Vector{qadic}(undef, e)
    for i = 1:e
      coeff = Qq()
      for j = 0:f-1
        coeff += (gQq^j)*xQp[2, j+1+(i-1)*f]
      end
      img_prim_elem[i] = coeff
    end
    f.prim_img = Kp(Qqx(img_prim_elem))
  end
  return nothing
end

################################################################################
#
#   Totally ramified case
#
################################################################################

function totally_ramified_completion(K::AnticNumberField, P::NfOrdIdl, precision::Int = 64)
  @assert precision > 0
  OK = order(P)
  @assert isprime(P)
  @assert nf(OK) == K
  @assert isone(degree(P))
  e = ramification_index(P)
  Qp = PadicField(minimum(P), precision)
  Zp = maximal_order(Qp)
  Zx = FlintZZ["x"][1]
  Qpx = PolynomialRing(Qp, "x")[1]
  u = uniformizer(P).elem_in_nf
  pows_u = powers(u, e-1)
  bK = basis_matrix(nf_elem[u*pows_u[end], gen(K)], FakeFmpqMat)
  append!(pows_u, map(elem_in_nf, basis(P^precision, copy = false)))
  MK = basis_matrix(pows_u, FakeFmpqMat)
  d = lcm(denominator(MK, copy = false), denominator(bK, copy = false))
  if d != denominator(MK, copy = false)
    mul!(MK.num, mK.num, divexact(d, denominator(MK, copy = false)))
  end
  if d != denominator(bK, copy = false)
    mul!(bK.num, bK.num, divexact(d, denominator(bK, copy = false)))
  end
  MZp = map_entries(Zp, MK.num)
  bZp = map_entries(Zp, bK.num)
  fl, xZp = can_solve_with_solution(MZp, bZp, side = :left)
  @assert fl
  coeffs_eisenstein = Vector{padic}(undef, e+1)
  for i = 1:e
    coeffs_eisenstein[i] = -xZp[1, i].x
  end
  coeffs_eisenstein[e+1] = one(Qp)
  pol_gen = Qpx(coeffs_eisenstein)
  Kp, gKp = eisenstein_extension(pol_gen, "a")
  img_prim_elem = Vector{padic}(undef, e)
  for i = 1:e
    img_prim_elem[i] = xZp[2, i].x
  end
  img = Kp(Qpx(img_prim_elem))
  if Nemo.precision(img) < Nemo.precision(Kp)
    img = newton_lift(Zx(defining_polynomial(K)), img, Nemo.precision(img), Nemo.precision(Kp))
  end
  completion_map = CompletionMap(K, Kp, img, u, precision)
  completion_map.P = P
  return Kp, completion_map
end

function setprecision!(f::CompletionMap{LocalField{padic, EisensteinLocalField}, LocalFieldElem{padic, EisensteinLocalField}}, new_prec::Int)
  if new_prec < f.precision
    K = domain(f)
    setprecision!(K, new_prec)
    setprecision!(base_field(K), new_prec)
    setprecision!(f.prim_img, new_prec)
  else
    #I need to increase the precision of the data
    P = prime(f)
    e = ramification_index(P)
    u = f.inv_img[2]
    Kp = codomain(f)
    ex, r = divrem(new_prec, ramification_index(P))
    if r > 0
      ex += 1
    end
    Qp = PadicField(prime(Kp), div(new_prec, e)+1)
    Zp = maximal_order(Qp)
    Qpx = PolynomialRing(Qp, "x")
    pows_u = powers(u, e-1)
    bK = basis_matrix(nf_elem[u*pows_u[end], gen(K)])
    append!(pows_u, map(elem_in_nf, basis(P^new_prec, copy = false)))
    MK = basis_matrix(pows_u)
    MQp = map_entries(Zp, MK)
    bQp = map_entries(Zp, bK)
    fl, xZp = can_solve_with_solution(MZp, bZp, side = :left)
    @assert fl
    coeffs_eisenstein = Vector{padic}(undef, e+1)
    for i = 1:e
      coeffs_eisenstein[i] = -xZp[1, i].x
    end
    coeffs_eisenstein[e+1] = one(Qp)
    pol_gen = Qpx(coeffs_eisenstein)
    Kp.defining_polynomial = pol_gen
    Kp.precision = new_prec
    #Should I update the traces too?
    img_prim_elem = Vector{padic}(undef, e)
    for i = 1:e
      img_prim_elem[i] = xZp[2, i].x
    end
    img = Kp(Qpx(img_prim_elem))
    f.prim_img = img
  end
  return nothing
end

################################################################################
#
#  Unramified case
#
################################################################################

function unramified_completion(K::AnticNumberField, P::NfOrdIdl, precision::Int = 64)
  OK = order(P)
  @assert isprime(P)
  @assert nf(OK) == K
  @assert isone(ramification_index(P))
  f = degree(P)
  p = minimum(P)
  Qq, gQq = QadicField(p, f, precision)
  Qp = PadicField(p, precision)
  Zp = maximal_order(Qp)
  q, mq = ResidueField(Qq)
  F, mF = ResidueFieldSmall(OK, P)
  mp = find_morphism(q, F)
  g = gen(q)
  gq_in_K = (mF\(mp(g))).elem_in_nf
  Zx = PolynomialRing(FlintZZ, "x")[1]
  pol_gq = lift(Zx, defining_polynomial(q))
  gq_in_K = _lift(gq_in_K, pol_gq, precision, P)
  #To compute the image of the primitive element, we use linear algebra if p is an index divisor
  #Hensel lifting otherwise
  if !isdefining_polynomial_nice(K) || isindex_divisor(OK, minimum(P))
    els = powers(gq_in_K, f-1)
    append!(els, map(elem_in_nf, basis(P^precision)))
    MK = basis_matrix(els, FakeFmpqMat)
    bK = basis_matrix(nf_elem[gen(K)], FakeFmpqMat)
    d = lcm(denominator(MK, copy = false), denominator(bK, copy = false))
    if d != denominator(MK, copy = false)
      mul!(MK.num, mK.num, divexact(d, denominator(MK, copy = false)))
    end
    if d != denominator(bK, copy = false)
      mul!(bK.num, bK.num, divexact(d, denominator(bK, copy = false)))
    end
    MZp = map_entries(Zp, MK.num)
    bZp = map_entries(Zp, bK.num)
    fl, xZp = can_solve_with_solution(MZp, bZp, side = :left)
    @assert fl
    img = Qq()
    for j = 0:f-1
      img += (gQq^j)*xZp[1, j+1].x
    end
  else
    el = mq\(mp\(mF(OK(gen(K)))))
    img = newton_lift(Zx(K.pol), el)
  end
  completion_map = CompletionMap(K, Qq, img, gq_in_K, precision)
  completion_map.P = P
  return Qq, completion_map
end

function setprecision!(f::CompletionMap{FlintQadicField, qadic}, new_prec::Int)
  Kp = codomain(f)
  setprecision!(Kp, new_prec)
  if new_prec < f.precision
    setprecision!(f.prim_img, new_prec)
  else
    P = prime(f)
    f = inertia_degree(P)
    gq, u = f.inv_img
    Zx = PolynomialRing(FlintZZ, "x")[1]
    q, mq = ResidueField(Kp)
    pol_gq = lift(Zx, defining_polynomial(q))
    gq = _increase_precision(gq, pol_gq, f.precision, new_prec, P)
    f.inv_img[1] = gq
    setprecision!(Qq, new_prec)
    #To increase the precision of the image of the primitive element, I use Hensel lifting
    f.prim_img = newton_lift(Zx(defining_polynomial(domain(f))), f.prim_img, new_prec, f.precision)
  end
  f.precision = new_prec
  return nothing
end


text:
#!/usr/bin/julia

using Pkg

installed() = Pkg.installed()

if !haskey(installed(),"ProjectAssistant")
    Pkg.add("ProjectAssistant")
    Pkg.update()
end

using ProjectAssistant
project_activate(abspath(@__DIR__,"../"))


text:
using AbstractAlgebra

# Function taken from the elimination project, to insert a reference
#------------------------------------------------------------------------------

function dict_to_poly(dict_monom::Dict{Array{Int, 1}, <: RingElem}, poly_ring::MPolyRing)
    builder = MPolyBuildCtx(poly_ring)
    for (monom, coef) in pairs(dict_monom)
        push_term!(builder, poly_ring.base_ring(coef), monom)
    end
    return finish(builder)
end

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

function parent_ring_change(poly::MPolyElem, new_ring::MPolyRing)
    """
    Converts a polynomial to a different polynomial ring
    Input
      - poly - a polynomial to be converted
      - new_ring - a polynomial ring such that every variable name
          appearing in poly appears among the generators
    Output: a polynomial in new_ring "equal" to poly
    """
    old_ring = parent(poly)
    # construct a mapping for the variable indices
    var_mapping = Array{Any, 1}()
    for u in symbols(old_ring)
        push!(
            var_mapping,
            findfirst(v -> (string(u) == string(v)), symbols(new_ring))
        )
    end
    builder = MPolyBuildCtx(new_ring)
    one_new = one(new_ring.base_ring)
    for term in zip(exponent_vectors(poly), coeffs(poly))
        exp, coef = term
        # detecting Singular.n_unknown
        if (typeof(coef) <: Singular.n_unknown) && !(typeof(one_new) <: Singular.n_unknown)
            coef = Singular.libSingular.julia(Singular.libSingular.cast_number_to_void(coef.ptr))
        end
        new_exp = [0 for _ in gens(new_ring)]
        for i in 1:length(exp)
            if exp[i] != 0
                if var_mapping[i] == nothing
                    throw(Base.ArgumentError("The polynomial contains a variable not present in the new ring $poly"))
                else
                    new_exp[var_mapping[i]] = exp[i]
                end
            end
        end
        push_term!(builder, new_ring.base_ring(coef), new_exp)
    end
    return finish(builder)
end

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

function str_to_var(s::String, ring::MPolyRing)
    ind = findfirst(v -> (string(v) == s), symbols(ring))
    if ind == nothing
        throw(Base.KeyError("Variable $s is not found in ring $ring"))
    end
    return gens(ring)[ind]
end

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

function var_to_str(v::MPolyElem)
    ind = findfirst(vv -> vv == v, gens(parent(v)))
    return string(symbols(parent(v))[ind])
end

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

function switch_ring(v::MPolyElem, ring::MPolyRing)
    """
    For a variable v, returns a variable in ring with the same name
    """
    ind = findfirst(vv -> vv == v, gens(parent(v)))
    return str_to_var(string(symbols(parent(v))[ind]), ring)
end

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

function extract_coefficients(poly::P, variables::Array{P, 1}) where P <: MPolyElem
    """
    Intput:
        poly - multivariate polynomial
        variables - a list of variables from the generators of the ring of p
    Output:
        dictionary with keys being tuples of length len(variables) and values being
        polynomials in the variables other than variables which are the coefficients
        at the corresponding monomials (in a smaller polynomial ring)
    """
    var_to_ind = Dict([(v, findfirst(e -> (e == v), gens(parent(poly)))) for v in variables])
    indices = [var_to_ind[v] for v in variables]

    coeff_vars = filter(v -> !(var_to_str(v) in map(var_to_str, variables)), gens(parent(poly)))
    new_ring, new_vars = PolynomialRing(base_ring(parent(poly)), map(var_to_str, coeff_vars))
    coeff_var_to_ind = Dict([(v, findfirst(e -> (e == v), gens(parent(poly)))) for v in coeff_vars])
    FieldType = typeof(one(base_ring(new_ring)))

    result = Dict{Array{Int, 1}, Dict{Array{Int, 1}, FieldType}}()

    for (monom, coef) in zip(exponent_vectors(poly), coeffs(poly))
        var_slice = [monom[i] for i in indices]
        if !haskey(result, var_slice)
            result[var_slice] = Dict{Array{Int, 1}, FieldType}()
        end
        new_monom = [0 for _ in 1:length(coeff_vars)]
        for i in 1:length(new_monom)
            new_monom[i] = monom[coeff_var_to_ind[coeff_vars[i]]]
        end
        result[var_slice][new_monom] = coef
    end

    return Dict(k => dict_to_poly(v, new_ring) for (k, v) in result)
end

max_stars_repo_name: pogudingleb/DifferentialAlgebra

max_issues_repo_name: pogudingleb/DifferentialAlgebra

max_forks_repo_name: pogudingleb/DifferentialAlgebra


text:
function qap(n, W, D, ::Val{:JuMP})
    model = JuMP.Model(CBLS.Optimizer)

    @variable(model, 1 ≤ X[1:n] ≤ n, Int)
    @constraint(model, X in AllDifferent())

    Σwd = p -> sum(sum(W[p[i], p[j]] * D[i, j] for j in 1:n) for i in 1:n)

    @objective(model, Min, ScalarFunction(Σwd))

    return model, X
end

"""
    qap(n, weigths, distances; modeler = :JuMP)

Modelize an instance of the Quadractic Assignment Problem with
- `n`: number of both facilities and locations
- `weights`: `Matrix` of the weights of each pair of facilities
- `distances`: `Matrix` of distances between locations
- `modeler`: Default to `:JuMP`. No other modeler available for now.

# From Wikipedia
There are a set of `n` facilities and a set of `n` locations. For each pair of locations, a distance is specified and for each pair of facilities a `weight` or flow is specified (e.g., the amount of supplies transported between the two facilities). The problem is to assign all facilities to different locations with the goal of minimizing the sum of the `distances` multiplied by the corresponding flows.
"""
qap(n, weigths, distances; modeler = :JuMP) = qap(n, weigths, distances, Val(modeler))