text:
bulk_types = Dict("n" => :ntype,
    "n-type" => :ntype,
    "ntype" => :ntype,
    "p-type" => :ptype,
    "ptype" => :ptype,
    "p" => :ptype  )
unit_conversion = Dict{String, Unitful.Units}(
    "nm" => u"nm", "um" => u"μm", "mm" => u"mm", "cm" => u"cm", "m" => u"m", #length
    "deg" => u"°","rad" => u"rad", #angle
    "V" => u"V", "kV" => u"kV", #Potential
    "K" => u"K", "Kelvic" => u"K", "C" => u"C")

include("Object.jl")
include("Passive.jl")
include("Contacts.jl")
include("Semiconductor.jl")
include("SolidStateDetector.jl")

"""
    SolidStateDetector{T}(filename::AbstractString)::SolidStateDetector{T} where {T <: SSDFloat}

Reads in a config-JSON file and returns an Detector struct which holds all information specified in the config file.
"""
function parse_config_file(filename::AbstractString)::Dict{Any,Any}
    if endswith(filename, ".toml")
        error("currently only .json and .yaml files are supported. We intend to add .toml support in the near future")
    elseif endswith(filename, ".json")
        dicttext = read(filename, String)
        parsed_dict = JSON.parse(dicttext)
    elseif endswith(filename, ".yaml")
        parsed_dict = YAML.load_file(filename)
    else
        error("currently only .json and .yaml files are supported.")
    end
    parsed_dict
end

function yaml2json_convert(filename::String)
    data = YAML.load(open(filename))
    open(replace(filename, ".yaml" => ".json"),"w") do f
            JSON.print(f, data, 4)
    end
end
function yaml2json(directory::String)# or filename
    if endswith(directory, ".yaml")
        yaml2json_convert(directory)
    else
        yamlfiles = filter(x->endswith(x,".yaml"),readdir(directory))
        for filename in yamlfiles
            yaml2json_convert(filename)
        end
    end
end
function SolidStateDetector{T}(filename::AbstractString)::SolidStateDetector{T} where {T <: SSDFloat}
    parsed_dict = parse_config_file(filename)
    return SolidStateDetector{T}(parsed_dict)
end

function SolidStateDetector(T::Type{<:AbstractFloat} = Float32, filename::AbstractString = SSD_examples[:InvertedCoax])::SolidStateDetector{T}
    SolidStateDetector{T}(filename)
end
function SolidStateDetector(filename::AbstractString)::SolidStateDetector{Float32}
    SolidStateDetector{Float32}(filename)
end

@inline in(point::CylindricalPoint, detector::SolidStateDetector) =
    contains(detector, point)

@inline in(point::StaticArray{Tuple{3}, <:Real}, detector::SolidStateDetector) =
    convert(CartesianPoint, point) in detector

@inline in(point::StaticArray{Tuple{3}, <:Quantity}, detector::SolidStateDetector) =
    to_internal_units(u"m", point) in detector

@inline in(point::CartesianPoint, detector::SolidStateDetector) =
    convert(CylindricalPoint, point) in detector

function get_important_points(c::SolidStateDetector{T}, s::Symbol)::Vector{T} where {T <: SSDFloat}
    imp::Vector{T} = []
    for contact in c.contacts
        for g in vcat(contact.geometry_positive,contact.geometry_negative)
            append!(imp, get_important_points(g, Val{s}()))
        end
    end
    return uniq(sort(imp))
end

function is_boundary_point(c::SolidStateDetector, pt::AbstractCoordinatePoint{T}, ax1::Vector{T}, ax2::Vector{T}, ax3::Vector{T})::Tuple{Bool, Real, Int} where {T <: SSDFloat}
    # if false #!(p in c)
    #     return false, T(0), 0
    # else
    for contact in c.contacts
        if in(pt, contact, ax1)
            return true, contact.potential, contact.id
        end
    end
    for external_part in c.external_parts
        if in(pt, external_part, ax1)
            return true, external_part.potential, external_part.id
        end
    end
    return false, 0.0, 0
    # end
end

# function point_type(c::SolidStateDetector{T}, p::CylindricalPoint{T})::Tuple{Symbol,Int} where T
#     for contact in c.contacts
#         if p in contact #|| geom_round(p) in contact #|| in(go_to_nearest_gridpoint(c,p), contact, c.rs) || in(go_to_nearest_gridpoint(c,geom_round(p)), contact, c.rs)
#             return :electrode, contact.id
#         end
#     end
#     sp = is_surface_point(c, p)
#     if sp[1]
#         for contact in c.contacts
#             if in(go_to_nearest_gridpoint(c,p), contact, c.rs) && abs(sum(sp[2])) > 1
#                 println("olo")
#                 return :electrode, contact.id
#             else
#                 return :floating_boundary, 0
#             end
#         end
#     elseif !(p in c)
#         return :outside, -1
#     else
#         return :bulk, -1
#     end
# end
function point_type(c::SolidStateDetector{T}, grid::Grid{T, 3}, p::CylindricalPoint{T})::Tuple{UInt8, Int, CartesianVector{T}} where {T <: SSDFloat}
    surface_normal::CartesianVector{T} = CartesianVector{T}(0, 0, 0) # need undef version for this
    for contact in c.contacts
        if p in contact || geom_round(p) in contact || in(searchsortednearest(grid, p), contact) || in(searchsortednearest(grid, geom_round(p)), contact)
            return CD_ELECTRODE::UInt8, contact.id, surface_normal
        end
    end
    on_surface, surface_normal = is_surface_point_and_normal_vector(c, p)
    if on_surface
        for contact in c.contacts
            if in(searchsortednearest(grid, p), contact) #&& abs(sum(sp[2])) > 1
                return CD_ELECTRODE::UInt8, contact.id, surface_normal
            else
                return CD_FLOATING_BOUNDARY::UInt8, -1, surface_normal
            end
        end
    elseif !(p in c)
        return CD_OUTSIDE::UInt8, -1, surface_normal
    else
        return CD_BULK::UInt8, -1, surface_normal
    end
end

# Point types for charge drift
const CD_ELECTRODE = 0x00
const CD_OUTSIDE = 0x01
const CD_BULK = 0x02
const CD_FLOATING_BOUNDARY = 0x04 # not 0x03, so that one could use bit operations here...

"""
    For charge drift...
"""
function point_type(c::SolidStateDetector{T}, grid::Grid{T, 3}, p::CartesianPoint{T})::Tuple{UInt8, Int, CartesianVector{T}} where {T <: SSDFloat}
    surface_normal::CartesianVector{T} = CartesianVector{T}(0, 0, 0) # need undef version for this
    for contact in c.contacts
        if p in contact #|| geom_round(p) in contact #|| in(go_to_nearest_gridpoint(c,p), contact, c.rs) || in(go_to_nearest_gridpoint(c,geom_round(p)), contact, c.rs)
            return CD_ELECTRODE::UInt8, contact.id, surface_normal
        end
    end
    on_surface, surface_normal = is_surface_point_and_normal_vector(c, p) # surface_normal::CartesianVector{T}
    if on_surface
        for contact in c.contacts
            if in(searchsortednearest(grid, p), contact)
                return CD_ELECTRODE::UInt8, contact.id, surface_normal
            else
                return CD_FLOATING_BOUNDARY::UInt8, -1, surface_normal
            end
        end
    elseif !(p in c)
        return CD_OUTSIDE::UInt8, -1, surface_normal
    else
        return CD_BULK::UInt8, -1, surface_normal
    end
end

# go_to_nearest_gridpoint
function searchsortednearest(grid::Grid{T, 3, :cylindrical}, pt::CylindricalPoint{T})::CylindricalPoint{T} where {T <: SSDFloat}
    idx1::Int = searchsortednearest(grid.axes[1].ticks, pt.r)
    idx2::Int = searchsortednearest(grid.axes[2].ticks, pt.φ)
    idx3::Int = searchsortednearest(grid.axes[3].ticks, pt.z)
    CylindricalPoint{T}(grid.axes[1].ticks[idx1], grid.axes[2].ticks[idx2], grid.axes[3].ticks[idx3])
end
function searchsortednearest(grid::Grid{T, 3, :cartesian}, pt::CartesianPoint{T})::CartesianPoint{T} where {T <: SSDFloat}
    idx1::Int = searchsortednearest(grid.axes[1].ticks, pt.x)
    idx2::Int = searchsortednearest(grid.axes[2].ticks, pt.y)
    idx3::Int = searchsortednearest(grid.axes[3].ticks, pt.z)
    CartesianPoint{T}(grid.axes[1].ticks[idx1], grid.axes[2].ticks[idx2], grid.axes[3].ticks[idx3])
end

function is_surface_point_and_normal_vector(c::SolidStateDetector{T}, p::CylindricalPoint{T})::Tuple{Bool, CartesianVector{T}} where {T <: SSDFloat}
    if !(p in c)
        return false, CartesianPoint{T}(0, 0, 0)
    end
    n::MVector{3,T} = @MVector T[0, 0, 0]
    look_around::Vector{Bool} = [   CylindricalPoint{T}(prevfloat(p.r), p.φ, p.z) in c,
                                    CylindricalPoint{T}(nextfloat(p.r), p.φ, p.z) in c,
                                    CylindricalPoint{T}(p.r, prevfloat(p.φ), p.z) in c,
                                    CylindricalPoint{T}(p.r, nextfloat(p.φ), p.z) in c,
                                    CylindricalPoint{T}(p.r, p.φ, prevfloat(p.z)) in c,
                                    CylindricalPoint{T}(p.r, p.φ, nextfloat(p.z)) in c]
    if !(false in look_around)
        return false , CartesianPoint{T}(n...)
    else
        if (look_around[1]==false) n[1] -= 1 end
        if (look_around[2]==false) n[1] += 1 end
        if (look_around[3]==false) n[2] -= 1 end
        if (look_around[4]==false) n[2] += 1 end
        if (look_around[5]==false) n[3] -= 1 end
        if (look_around[6]==false) n[3] += 1 end
        # println(look_around , " " , n)
        Rα::SMatrix{3,3,T} = @SArray([cos(p.φ) -1*sin(p.φ) 0;sin(p.φ) cos(p.φ) 0;0 0 1])
        # return true, geom_round(CartesianVector{T}((Rα * n)...))
        return true, CartesianVector{T}((Rα * n)...)
    end
end
function is_surface_point_and_normal_vector(c::SolidStateDetector{T}, p::CartesianPoint{T})::Tuple{Bool, CartesianVector{T}} where T
    if !(p in c) # is this necessary?
        return false, CartesianVector{T}(0, 0, 0)
    end
    n::MVector{3,T} = @MVector T[0.0,0.0,0.0]
    look_around::Vector{Bool} = [   CartesianPoint{T}(prevfloat(p.x), p.y, p.z) in c,
                                    CartesianPoint{T}(nextfloat(p.x), p.y, p.z) in c,
                                    CartesianPoint{T}(p.x, prevfloat(p.y), p.z) in c,
                                    CartesianPoint{T}(p.x, nextfloat(p.y), p.z) in c,
                                    CartesianPoint{T}(p.x, p.y, prevfloat(p.z)) in c,
                                    CartesianPoint{T}(p.x, p.y, nextfloat(p.z)) in c]
    if !(false in look_around)
        return false, n
    else
        if (look_around[1] == false) n[1] -= 1 end
        if (look_around[2] == false) n[1] += 1 end
        if (look_around[3] == false) n[2] -= 1 end
        if (look_around[4] == false) n[2] += 1 end
        if (look_around[5] == false) n[3] -= 1 end
        if (look_around[6] == false) n[3] += 1 end
        # println(look_around , " " , n)
        return true, n
    end
end

function get_charge_density(sc::Semiconductor{T}, pt::AbstractCoordinatePoint{T})::T where {T <: SSDFloat}
    get_charge_density(sc.charge_density_model, pt)
end

function get_ρ_and_ϵ(pt::AbstractCoordinatePoint{T}, ssd::SolidStateDetector{T})::Tuple{T, T} where {T <: SSDFloat}
    if in(pt,ssd.semiconductors)
        for sc in ssd.semiconductors
            if in(pt, sc)
                ρ::T = get_charge_density(sc, pt) * elementary_charge
                ϵ::T = sc.material.ϵ_r
                return ρ, ϵ
            end
        end
    elseif in(pt,ssd.passives)
        for ep in ssd.passives
            if pt in ep
                ρ = 0
                ϵ = ep.material.ϵ_r
            end
        end
        return ρ, ϵ
    else
        ρ = 0
        ϵ = ssd.medium.ϵ_r
        return ρ, ϵ
    end
end

function write_grid_to_detector!(ssd::SolidStateDetector{T}, grid::Grid{T, 3, :cylindrical})::Nothing where {T}
    ssd.rs = grid[:r].ticks
    ssd.φs = grid[:φ].ticks
    ssd.zs = grid[:z].ticks
    nothing
end

function set_pointtypes_and_fixed_potentials!(pointtypes::Array{PointType, N}, potential::Array{T, N},
        grid::Grid{T, N, :cylindrical}, ssd::SolidStateDetector{T}; weighting_potential_contact_id::Union{Missing, Int} = missing)::Nothing where {T <: SSDFloat, N}

    channels::Array{Int, 1} = if !ismissing(weighting_potential_contact_id)
        [weighting_potential_contact_id]
    else
        Int[]
    end

    axr::Vector{T} = grid[:r].ticks
    axφ::Vector{T} = grid[:φ].ticks
    axz::Vector{T} = grid[:z].ticks

    for iz in axes(potential, 3)
        z::T = axz[iz]
        for iφ in axes(potential, 2)
            φ::T = axφ[iφ]
            for ir in axes(potential, 1)
                r::T = axr[ir]
                pt::CylindricalPoint{T} = CylindricalPoint{T}( r, φ, z )

                if in(pt, ssd)
                    pointtypes[ ir, iφ, iz ] += pn_junction_bit
                end

            end
        end
    end
    for contact in ssd.contacts
        pot::T = if ismissing(weighting_potential_contact_id)
            contact.potential
        else
            contact.id == weighting_potential_contact_id ? 1 : 0
        end
        contact_gridpoints = paint_object(contact,grid)
        for gridpoint in contact_gridpoints
            potential[ gridpoint... ] = pot
            pointtypes[ gridpoint... ] = zero(PointType)
        end
    end
    nothing
end

function set_pointtypes_and_fixed_potentials!(pointtypes::Array{PointType, N}, potential::Array{T, N},
    grid::Grid{T, N, :cartesian}, ssd::SolidStateDetector{T}; weighting_potential_contact_id::Union{Missing, Int} = missing)::Nothing where {T <: SSDFloat, N}

    channels::Array{Int, 1} = if !ismissing(weighting_potential_contact_id)
        [weighting_potential_contact_id]
    else
        Int[]
    end

    axx::Vector{T} = grid[:x].ticks
    axy::Vector{T} = grid[:y].ticks
    axz::Vector{T} = grid[:z].ticks

    for iz in axes(potential, 3)
        z::T = axz[iz]
        for iy in axes(potential, 2)
            y::T = axy[iy]
            for ix in axes(potential, 1)
                x::T = axx[ix]
                pt::CartesianPoint{T} = CartesianPoint{T}( x, y, z )

                if in(pt, ssd)
                    pointtypes[ ix, iy, iz ] += pn_junction_bit
                end

            end
        end
    end
    for contact in ssd.contacts
        pot::T = if ismissing(weighting_potential_contact_id)
            contact.potential
        else
            contact.id == weighting_potential_contact_id ? 1 : 0
        end
        contact_gridpoints = paint_object(contact,grid)
        for gridpoint in contact_gridpoints
            potential[ gridpoint... ] = pot
            pointtypes[ gridpoint... ] = zero(PointType)
        end
    end
    nothing
end

function json_to_dict(inputfile::String)::Dict
    parsed_json_file = Dict()
    open(inputfile,"r") do f
        global parsed_json_file
        dicttext = readstring(f)
        parsed_json_file = JSON.parse(dicttext)
    end
    return parsed_json_file
end
function bounding_box(d::SolidStateDetector{T})::NamedTuple where T
    (
    r_range = ClosedInterval{T}(0.0,d.crystal_radius),
    φ_range = ClosedInterval{T}(0.0,2π),
    z_range = ClosedInterval{T}(0.0,d.crystal_length)
    )
end

function Grid(  detector::SolidStateDetector{T, :cylindrical};
                init_grid_spacing::Vector{<:Real} = [0.005, deg2rad(5.0), 0.005],
                for_weighting_potential::Bool = false,
                full_2π::Bool = false)::CylindricalGrid{T} where {T}

    init_grid_spacing::Vector{T} = T.(init_grid_spacing)

    important_r_points::Vector{T} = get_important_points(detector, :r)
    important_φ_points::Vector{T} = get_important_points(detector, :φ)
    important_z_points::Vector{T} = get_important_points(detector, :z)

    push!(important_r_points, detector.world.intervals[1].right)
    push!(important_r_points, detector.world.intervals[1].left)
    important_r_points = uniq(sort(important_r_points))
    push!(important_z_points, detector.world.intervals[3].right)
    push!(important_z_points, detector.world.intervals[3].left)
    important_z_points = uniq(sort(important_z_points))
    push!(important_φ_points, detector.world.intervals[2].right)
    push!(important_φ_points, detector.world.intervals[2].left)
    important_φ_points = uniq(sort(important_φ_points))

    # r
    L, R, BL, BR = get_boundary_types(detector.world.intervals[1])
    int_r = Interval{L, R, T}(detector.world.intervals[1].left, detector.world.intervals[1].right)
    ax_r::DiscreteAxis{T, BL, BR} = DiscreteAxis{BL, BR}(int_r, step = init_grid_spacing[1])
    rticks::Vector{T} = merge_axis_ticks_with_important_ticks(ax_r, important_r_points, atol=init_grid_spacing[1]/4)
    ax_r = DiscreteAxis{T, BL, BR}(int_r, rticks)

    # φ
    L, R, BL, BR = get_boundary_types(detector.world.intervals[2])
    int_φ = Interval{L, R, T}(detector.world.intervals[2].left, detector.world.intervals[2].right)
    if full_2π == true || (for_weighting_potential && (detector.world.intervals[2].left != detector.world.intervals[2].right))
        L, R, BL, BR = :closed, :open, :periodic, :periodic
        int_φ = Interval{L, R, T}(0, 2π)
    end
    ax_φ = if int_φ.left == int_φ.right
        DiscreteAxis{T, BL, BR}(int_φ, T[int_φ.left])
    else
        DiscreteAxis{BL, BR}(int_φ, step = init_grid_spacing[2])
    end
    if length(ax_φ) > 1
        φticks::Vector{T} = merge_axis_ticks_with_important_ticks(ax_φ, important_φ_points, atol=deg2rad(init_grid_spacing[2])/4)
        ax_φ = typeof(ax_φ)(int_φ, φticks)
    end
    if isodd(length(ax_φ)) && length(ax_φ) > 1 # must be even
        int_φ = ax_φ.interval
        φticks = ax_φ.ticks
        push!(φticks, geom_round((φticks[end] + φticks[end-1]) * 0.5))
        sort!(φticks)
        ax_φ = typeof(ax_φ)(int_φ, φticks) # must be even
    end
    if length(ax_φ) > 1
        @assert iseven(length(ax_φ)) "CylindricalGrid must have even number of points in φ."
    end

    #z
    L, R, BL, BR = get_boundary_types(detector.world.intervals[3])
    int_z = Interval{L, R, T}(detector.world.intervals[3].left, detector.world.intervals[3].right)
    ax_z::DiscreteAxis{T, BL, BR} = DiscreteAxis{BL, BR}(int_z, step = init_grid_spacing[3])
    zticks::Vector{T} = merge_axis_ticks_with_important_ticks(ax_z, important_z_points, atol=init_grid_spacing[3]/2)
    ax_z = typeof(ax_z)(int_z, zticks)
    if isodd(length(ax_z)) # must be even
        int_z = ax_z.interval
        zticks = ax_z.ticks
        push!(zticks, geom_round((zticks[end] + zticks[end-1]) * 0.5))
        sort!(zticks)
        ax_z = typeof(ax_z)(int_z, zticks) # must be even
    end
    @assert iseven(length(ax_z)) "CylindricalGrid must have even number of points in z."

    return CylindricalGrid{T}( (ax_r, ax_φ, ax_z) )
end

function Grid(  detector::SolidStateDetector{T, :cartesian};
                init_grid_spacing::Vector{<:Real} = [0.001, 0.001, 0.001])::CartesianGrid3D{T} where {T}

    important_x_points::Vector{T} = get_important_points(detector, :x)
    important_y_points::Vector{T} = get_important_points(detector, :y)
    important_z_points::Vector{T} = get_important_points(detector, :z)

    init_grid_spacing::Vector{T} = T.(init_grid_spacing)

    # x
    L, R, BL, BR = get_boundary_types(detector.world.intervals[1])
    int_x = Interval{L, R, T}(detector.world.intervals[1].left, detector.world.intervals[1].right)
    ax_x::DiscreteAxis{T, BL, BR} = DiscreteAxis{BL, BR}(int_x, step = init_grid_spacing[1])
    xticks::Vector{T} = merge_axis_ticks_with_important_ticks(ax_x, important_x_points, atol = init_grid_spacing[1] / 2)
    ax_x = typeof(ax_x)(int_x, xticks)
    if isodd(length(ax_x)) # RedBlack dimension must be of even length
        xticks = ax_x.ticks
        push!(xticks, geom_round((xticks[end] + xticks[end-1]) * 0.5))
        sort!(xticks)
        ax_x = DiscreteAxis{T, BL, BR}(int_x, xticks) # must be even
    end
    @assert iseven(length(ax_x)) "CartesianGrid3D must have even number of points in z."

    # y
    L, R, BL, BR = get_boundary_types(detector.world.intervals[2])
    int_y = Interval{L, R, T}(detector.world.intervals[2].left, detector.world.intervals[2].right)
    ax_y::DiscreteAxis{T, BL, BR} = DiscreteAxis{BL, BR}(int_y, step = init_grid_spacing[2])
    yticks::Vector{T} = merge_axis_ticks_with_important_ticks(ax_y, important_y_points, atol = init_grid_spacing[2] / 2)
    ax_y = typeof(ax_y)(int_y, yticks)

    # z
    L, R, BL, BR = get_boundary_types(detector.world.intervals[3])
    int_z = Interval{L, R, T}(detector.world.intervals[3].left, detector.world.intervals[3].right)
    ax_z::DiscreteAxis{T, BL, BR} = DiscreteAxis{BL, BR}(int_z, step = init_grid_spacing[3])
    zticks::Vector{T} = merge_axis_ticks_with_important_ticks(ax_z, important_z_points, atol = init_grid_spacing[3] / 2)
    ax_z = typeof(ax_z)(int_z, zticks)

    return CartesianGrid3D{T}( (ax_x, ax_y, ax_z) )
end

include("plot_recipes.jl")


text:
function map_socs!(v::Vector{Float64}, cones::ConicSets, a::AuxiliaryData)
    cont = 0
    for (idx, sdp) in enumerate(cones.sdpcone)
        cont += div(sdp.sq_side * (sdp.sq_side + 1), 2)
    end
    sizehint!(a.soc_v, length(cones.socone))
    sizehint!(a.soc_s, length(cones.socone))
    for (idx, soc) in enumerate(cones.socone)
        len = soc.len
        push!(a.soc_s, view(v, cont + 1))
        push!(a.soc_v, view(v, cont + 2:cont + len))
        cont += len
    end
    return nothing
end

ext: jl


text:
#==
    Note on implementation of factorisition types:
    The general strategy is to create a named factorization type which wraps other factorizations.
    This allows us to define `getproperty` methods which return `NamedDimsArray`s around the
    corresponding properties of the parent factorization types.

    NOTE: Maybe we could just have a `Named` type that works with arrays and factorizations
    kind of like `Adjoint`?
==#

struct NamedFactorization{L, T, F<:Factorization{T}} <: Factorization{T}
    fact::F
end

function NamedFactorization{L}(fact::F) where {L, T, F <: Factorization{T}}
    return NamedFactorization{L,T,F}(fact)
end

# Need parent to explicitly call `getfield` to avoid infinitely calling `getproperty`
Base.parent(named::NamedFactorization) = getfield(named, :fact)
Base.size(named::NamedFactorization) = size(parent(named))
Base.propertynames(named::NamedFactorization; kwargs...) = propertynames(parent(named))

# Factorization type specific initial iterate calls
Base.iterate(named::NamedFactorization{L, T, <:LU}) where {L, T} = (named.L, Val(:U))
Base.iterate(named::NamedFactorization{L, T, <:LQ}) where {L, T} = (named.L, Val(:Q))
Base.iterate(named::NamedFactorization{L, T, <:SVD}) where {L, T} = (named.U, Val(:S))
function Base.iterate(
    named::NamedFactorization{L, T, <:Union{QR, LinearAlgebra.QRCompactWY, QRPivoted}}
) where {L, T}
    return (named.Q, Val(:R))
end

# Generic iterate follow up calls
function Base.iterate(named::NamedFactorization, st::Val{D}) where D
    r = iterate(parent(named), st)
    r === nothing && return nothing
    return (getproperty(named, D), last(r))
end

# Convenience constructors
for func in (:lu, :lu!, :lq, :lq!, :svd, :svd!, :qr, :qr!)
    @eval begin
        function LinearAlgebra.$func(nda::NamedDimsArray{L, T}, args...; kwargs...) where {L, T}
            return NamedFactorization{L}($func(parent(nda), args...; kwargs...))
        end
    end
end

# `getproperty` wrappers
## LU
function Base.getproperty(fact::NamedFactorization{L, T, <:LU}, d::Symbol) where {L, T}
    inner = getproperty(parent(fact), d)
    n1, n2 = L
    if d === :L
        return NamedDimsArray{(n1, :_)}(inner)
    elseif d === :U
        return NamedDimsArray{(:_, n2)}(inner)
    elseif d === :P
        perm_matrix_labels = (first(L), first(L))
        return NamedDimsArray{perm_matrix_labels}(inner)
    elseif d === :p
        perm_vector_labels = (first(L),)
        return NamedDimsArray{perm_vector_labels}(inner)
    else
        return inner
    end
end

## LQ
function Base.getproperty(fact::NamedFactorization{L, T, <:LQ}, d::Symbol) where {L, T}
    inner = getproperty(parent(fact), d)
    n1, n2 = L
    if d === :L
        return NamedDimsArray{(n1, :_)}(inner)
    elseif d === :Q
        return NamedDimsArray{(:_, n2)}(inner)
    else
        return inner
    end
end

## svd
function Base.getproperty(fact::NamedFactorization{L, T, <:SVD}, d::Symbol) where {L, T}
    inner = getproperty(parent(fact), d)
    n1, n2 = L
    # Naming based off the SVD visualization on wikipedia
    # # https://en.wikipedia.org/wiki/File:Singular_value_decomposition_visualisation.svg
    if d === :U
        return NamedDimsArray{(n1, :_)}(inner)
    elseif d === :V
        return NamedDimsArray{(n2, :_)}(inner)
    elseif d === :Vt
        return NamedDimsArray{(:_, n2)}(inner)
    else # :S
        return inner
    end
end

## qr
function Base.getproperty(
    fact::NamedFactorization{L, T, F},
    d::Symbol
) where {L, T, F<:Union{QR, LinearAlgebra.QRCompactWY, QRPivoted}}
    inner = getproperty(parent(fact), d)
    n1, n2 = L
    if d === :Q
        return NamedDimsArray{(n1,:_)}(inner)
    elseif d === :R
        return NamedDimsArray{(:_, n2)}(inner)
    elseif F <: QRPivoted && d === :P
        return NamedDimsArray{(n2, n2)}(inner)
    elseif F <: QRPivoted && d === :p
        return NamedDimsArray{(n2,)}(inner)
    else
        return inner
    end
end

function LinearAlgebra.:\(
    fact::NamedFactorization{L,T,F}, nda::NamedDimsArray{W}
) where {L,T,F<:Factorization{T},W}
    n1, n2 = L
    n1 != W[1] && throw(
        DimensionMismatch(
            "Mismatched dimensions with factorization: $L and NamedDimsArray: $W"
        ),
    )
    return NamedDimsArray{(n2,)}(LinearAlgebra.:\(parent(fact), parent(nda)))
end

function LinearAlgebra.:\(
    fact::NamedFactorization{L,T,F}, nda::AbstractVector
) where {L,T,F<:Factorization{T}}
    n1, n2 = L
    return NamedDimsArray{(n2,)}(LinearAlgebra.:\(parent(fact), nda))
end

# Specialised routines for \ often do in-place ops that result in the nameddim populated from B
# Leading to an incorrect named-dim
# We also unname here because that handles wrapper types
for S in (UpperTriangular, LowerTriangular)
    @eval begin
        function LinearAlgebra.:\(
            A::$S{T,<:NamedDimsArray{L}}, B::AbstractVector
        ) where {L,T}
            n1, n2 = L
            return NamedDimsArray{(n2,)}(LinearAlgebra.:\($S(unname(A)), parent(B)))
        end
    end
end

# Diagonal on a nameddim presently loses its nameddimsness. So just pass through for now.
LinearAlgebra.:\(A::Diagonal, B::NamedDimsArray) = LinearAlgebra.:\(A, parent(B))


text:
type SAGAcbsd <: StepData
	A
	y
	m::Int64
	numClasses::Int64
	c::Vector{Int64}
	batchSize::Int64
	getStep::Function
	stepString::AbstractString
	function SAGAcbsd(numVars::Int64,y, stepString::AbstractString,numDp::Int64, batchSize::Int64,c)
		A = sum(y)
		new(A, y, numDp, size(y)[1],c,batchSize, SAGAcbComputation, stepString)
	end
end

function SAGAcb(numVars::Int64, y,batchSize::Int64, c, numDP::Int64)
	SAGAcbsd(numVars,y, "SAGAcb.$batchSize", numDP,batchSize,c)
end

function SAGAcbComputation(x, k, gnum, sd::SAGAcbsd, problem::Problem; outputLevel =0 )

	outputLevel >0 && println("SAGAcbComputation started")

	B=zeros(size(x))

	sumy=zeros(size(x))

	indices = sample(1:sd.m,sd.batchSize;replace=false)
	nummm=0
	for i in indices

		B +=sd.y[sd.c[i]]
		(f,sampleG) = (problem.getSampleFunctionAt(i))(x)
		sd.y[sd.c[i]]=sampleG
		sumy += sampleG
		nummm +=1
		gnum += 1

		#step 3
	end
	outputLevel >0 && println("Acutal step determination")

	step =((1/sd.numClasses)*sd.A) - (1/nummm)*(B - sumy)

	outputLevel >0 && println("Updating current sum")
	sd.A=sd.A-B+sumy
	outputLevel >0 && println("SAGAcbComputation finished")
	(step,gnum)

end


text:
using Gnuplot
x = y = -15:0.33:15
fz(x,y) = sin.(sqrt.(x.^2 + y.^2))./sqrt.(x.^2+y.^2)
fxy = [fz(x,y) for x in x, y in y]
@gsp x y fxy "w l lc 'orange'"
save(term="pngcairo size 800,800", output="./code/plt3d_ex3.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
export AmoebaFiber2D, AmoebaFiber3D, CoamoebaFiber2D, CoamoebaFiber3D, ContourFiber2D,
    ImaginaryFiber2D, ImaginaryFiber3D

abstract type AbstractFiber{T, M, N} end

abstract type AbstractAmoebaFiber{T, N} <: AbstractFiber{T, 2, N} end
"""
An `AmoebaFiber2D` is a representation of the fiber of the amoeba ``\\mathcal{A}_f``
at ``(w_1,w_2)`` where `f` is a bivariate polynomial.
"""
mutable struct AmoebaFiber2D{T, E1, E2} <: AbstractAmoebaFiber{T, 2}
    f_re::SP.Polynomial{T, E1, Nothing}
    f_im::SP.Polynomial{T, E2, Nothing}

    # Tempory storage
    v::NTuple{2, Float64} # exp(w1), exp(w2)
    U::Matrix{T}
    x::SVector{4, Float64}

    # We use this to scale the system to improve the numerics
    exponents::NTuple{2, Matrix{Int}}
    original_coefficients::NTuple{2, Vector{T}}
end

"""
An `AmoebaFiber3D` is a representation of the fiber of the amoeba ``\\mathcal{A}_f``
at ``(w_1, w_2, w_3)`` where `f` is a trivariate polynomial.
"""
mutable struct AmoebaFiber3D{T, E1, E2} <: AbstractAmoebaFiber{T, 3}
    f_re::SP.Polynomial{T, E1, Nothing}
    f_im::SP.Polynomial{T, E2, Nothing}

    # Tempory storage
    v::NTuple{3, Float64} # exp(w1), exp(w2), exp(w3)
    U::Matrix{T}
    x::SVector{6, Float64}

    # We use this to scale the system to improve the numerics
    exponents::NTuple{2, Matrix{Int}}
    original_coefficients::NTuple{2, Vector{T}}
end

abstract type AbstractCoamoebaFiber{T, N} <: AbstractFiber{T, 2, N} end

"""
An `CoamoebaFiber2D` is a representation of the fiber of the coamoeba ``\\mathcal{A}^{\\prime}_f``
at ``(θ_1, θ_2)`` where `f` is a bivariate polynomial.
"""
mutable struct CoamoebaFiber2D{T, E1, E2} <: AbstractCoamoebaFiber{T, 2}
    f_re::SP.Polynomial{T, E1, Nothing}
    f_im::SP.Polynomial{T, E2, Nothing}

    # temporary storage
    sincosθ::NTuple{4, Float64}
    U::Matrix{T}
    x::SVector{4, Float64}
end

mutable struct CoamoebaFiber3D{T, E1, E2} <: AbstractCoamoebaFiber{T, 3}
    f_re::SP.Polynomial{T, E1, Nothing}
    f_im::SP.Polynomial{T, E2, Nothing}

    # temporary storage
    sincosθ::NTuple{6, Float64}
    U::Matrix{T}
    x::SVector{6, Float64}
end

abstract type AbstractContourFiber{T, N} <: AbstractFiber{T, 3, N} end
"""
An `ContourFiber2D` is a representation of the fiber of the contour of the amoeba ``\\mathcal{A}^{\\prime}_f``
at ``(w_1, w_2)`` along a onedimensional affine subspace
where `f` is a bivariate polynomial.
"""
mutable struct ContourFiber2D{T, E1, E2, E3} <: AbstractContourFiber{T, 3}
    f_re::SP.Polynomial{T, E1, Nothing}
    f_im::SP.Polynomial{T, E2, Nothing}
    g::SP.Polynomial{T, E3, Nothing}
    exponents::NTuple{3, Matrix{Int}}
    v_i::Float64
    x::SVector{4, Float64}
    U::Matrix{T}
    w1_fixed::Bool
    original_coefficients::NTuple{3, Vector{T}}
end

abstract type AbstractImaginaryFiber{T, N} <: AbstractFiber{T, 2, N} end

mutable struct ImaginaryFiber2D{T, E1, E2} <: AbstractImaginaryFiber{T, 2}
    f_re::SP.Polynomial{T, E1, Nothing}
    f_im::SP.Polynomial{T, E2, Nothing}

    y::NTuple{2, T}
    U::Matrix{T}
end

mutable struct ImaginaryFiber3D{T, E1, E2} <: AbstractImaginaryFiber{T, 3}
    f_re::SP.Polynomial{T, E1, Nothing}
    f_im::SP.Polynomial{T, E2, Nothing}

    y::NTuple{3, T}
    U::Matrix{T}
end

export AbstractAlgorithm, AbstractGridAlgorithm, ArchTrop, Simple, Greedy, Coarse

abstract type AbstractAlgorithm end
abstract type AbstractGridAlgorithm <: AbstractAlgorithm end

struct ArchTrop <: AbstractGridAlgorithm end
struct Coarse <: AbstractGridAlgorithm end
struct Simple <: AbstractGridAlgorithm end
struct Greedy <: AbstractGridAlgorithm end
struct Polygonal <: AbstractAlgorithm end

export NewtonPolygon
"""
    NewtonPolygon

Representation of a [Newton polygon](https://en.wikipedia.org/wiki/Newton_polygon) of a polynomial.

## Fields
- `polynomial::P`: the polynomial from which the NewtonPolygon was constructed.
- `lattices::Vector{SVector{2,Int}}`: the support of the polynomial `p`.
- `vertices::Vector{Int}`: the vertex indicdes of the Newton polygon.
- `facets::Vector{SVector{2,Int}}`: the facets of the Newton polygon, i.e. it's convex hull.
- `subdivision::Vector{SVector{3,Int}}`: the simplices of the regular subdivision of the
Newton polygon
"""
struct NewtonPolygon
    lattices::Vector{SVector{2, Int}}
    vertices::Vector{Int}
    facets::Vector{SVector{2, Int}}
    subdivision::Vector{Vector{Int}}
end

const Range = StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}

abstract type AbstractGrid{N} end
"""
A structure representing a 2D Grid.
"""
struct Grid2D <: AbstractGrid{2}
    xrange::Range
    yrange::Range
end

"""
A structure representing a 3D Grid.
"""
struct Grid3D <: AbstractGrid{3}
    xrange::Range
    yrange::Range
    zrange::Range
end

abstract type AbstractBitmap{N} end

export Bitmap2D, Bitmap3D
"""
A structure which holds a `BitMatrix` indicating active / inactive pixels and associated
world coordinates `xmin`, `xmax`, `ymin`, `ymax`.
"""
struct Bitmap2D <: AbstractBitmap{2}
    data::BitArray{2}
    grid::Grid2D
end

"""
Structute which holds a `BitArray{3}` indicating active / inactive voxels and associated
world coordinates `xmin`, `xmax`, `ymin`, `ymax`, `zmin`, `zmax`.
"""
struct Bitmap3D <: AbstractBitmap{3}
    data::BitArray{3}
    grid::Grid3D
end

export Tropical, TropicalCurve

"""
A `Tropical` number is an element of the semi-ring (ℝ ∪ {-∞}, ⊕, ⊙).
"""
struct Tropical{T<:Real} <: Number
    val::T
    isinf::Bool

    function Tropical{T}(val::T, isinf::Bool) where {T<:Real}
        new(val, isinf)
    end
    Tropical{T}(t::Tropical) where {T<:Real} = new(convert(T, t.val), t.isinf)
end

const AbstractTropicalPolynomial{T} = MP.AbstractPolynomial{Tropical{T}}

"""
A `TropicalCurve` is a tropical hypersurface of a bivariate polynomial.
"""
struct TropicalCurve{P<:AbstractTropicalPolynomial}
    polynomial::P
    vertices::Vector{SVector{2, Float64}}
    segments::Vector{SVector{2, Int64}}
    halfrays::Vector{Tuple{Int64, SVector{2, Float64}}}
end

export Spine2D, ComponentComplement, PolygonalAmoeba

"""
Represents a component of the complement of the amoeba.
It contains information about the order of the component, whether it is bounded
and one point out of the component.
"""
struct ComponentComplement{N}
    order::SVector{N, Int}
    x::SVector{N, Float64}
    bounded::Bool
end

"""
A `Spine2D` represents the spine of a two-dimensional amoeba. This also contains
informations about the components of the complement and the approximation used to obtain
this spine.
"""
struct Spine2D{P<:AbstractTropicalPolynomial}
    ronkin_polynomial::P
    curve::TropicalCurve{P}
    components_complement::Vector{ComponentComplement{2}}
    amoeba_approximation::Bitmap2D
end

"""
    PolygonalAmoeba

Holding the result of call to [amoeba](@ref) with the `Polygonal()` algorithm.
"""
struct PolygonalAmoeba{P}
    spine::Spine2D{P}
    polygons::Vector{NTuple{4, SVector{2, Float64}}}
    domain::NTuple{4, Float64}
    # qau
    error_bound::Float64
    iterations::Int
end

export Coamoeba
"""
    Coamoeba{N}

Holding a `BitArray` represeting the coamoeba embedded in ``[0,2π]^N``.
"""
struct Coamoeba{N}
    data::BitArray{N}
end

# We introduce some aliases
const Point = SVector{2, Float64}
const Direction = SVector{2, Float64}
const Pillar = Tuple{Point, Direction}


text:
#jl #! format: off
# Line Modeling simulation with [PowerSimulationsDynamics.jl](https://github.com/NREL-SIIP/PowerSimulationsDynamics.jl)

# **Originally Contributed by**: Rodrigo Henriquez and José Daniel Lara

# ## Introduction

# This tutorial will introduce an example of considering dynamic lines in `PowerSimulationsDynamics`.
# Note that this tutorial is for `PowerSimulationsDynamics`.

# This tutorial presents a simulation of a three-bus system, with an infinite bus (represented as a
# voltage source behind an impedance) at bus 1, a one d- one q- machine on bus 2 and an inverter
# of 19 states, as a virtual synchronous machine at bus 3. The perturbation will be the trip of
# two of the three circuits (triplicating its resistance and impedance) of the line that connects
# bus 1 and bus 3. This case also consider a dynamic line model for connection between buses
# 2 and 3. We will compare it against a system without dynamic lines.

# In addition, note that `PowerSimulationsDynamics` will convert ConstantPower loads to RLC
# loads for transient simulations.

# It is recommended to check `Tutorial 1: OMIB` first, since that includes more details and
# explanations on all definitions and functions.

# # Step 1: Package Initialization
using SIIPExamples
using PowerSimulationsDynamics
PSID = PowerSimulationsDynamics
using PowerSystems
using Sundials
using Plots

# # Step 2: Data creation
file_dir = joinpath(
    dirname(dirname(pathof(SIIPExamples))),
    "script",
    "4_PowerSimulationsDynamics_examples",
    "Data",
)
threebus_sys = System(joinpath(file_dir, "threebus_sys.json"));

# In addition, we will create a new copy of the system on which we will simulate the same
# case, but will consider dynamic lines:

threebus_sys_dyn = deepcopy(threebus_sys);

# # Step 3: Create the fault and simulation on the Static Lines system

# First, we construct the perturbation, by properly computing the new Ybus on the system:

#Make a copy of the original system
sys2 = deepcopy(threebus_sys)
#Triplicates the impedance of the line named "BUS 1-BUS 3-i_1"
fault_branches = get_components(ACBranch, sys2)
for br in fault_branches
    if get_name(br) == "BUS 1-BUS 3-i_1"
        br.r = 3 * br.r
        br.x = 3 * br.x
        b_new = (from = br.b.from / 3, to = br.b.to / 3)
        br.b = b_new
    end
end
#Obtain the new Ybus
Ybus_fault = Ybus(sys2).data
#Define Fault: Change of YBus
Ybus_change = NetworkSwitch(
    1.0, #change at t = 1.0
    Ybus_fault, #New YBus
);

# Now, we construct the simulation:

#Time span of our simulation
tspan = (0.0, 30.0)

#Define Simulation
sim = PSID.Simulation(
    PSID.ImplicitModel, #Type of model used
    threebus_sys, #system
    pwd(), #folder to output results
    tspan, #time span
    Ybus_change, #Type of perturbation
)

# We can obtain the initial conditions as:

#Will print the initial states. It also give the symbols used to describe those states.
print_device_states(sim)
#Will export a dictionary with the initial condition values to explore
x0_init = PSID.get_initial_conditions(sim)

# # Step 4: Run the simulation of the Static Lines System

#Run the simulation
PSID.execute!(
    sim, #simulation structure
    IDA(), #Sundials DAE Solver
    dtmax = 0.02, #Maximum step size
)

# # Step 5: Store the solution

series2 = get_voltage_magnitude_series(sim, 102)
zoom = [
    (series2[1][ix], series2[2][ix]) for
    (ix, s) in enumerate(series2[1]) if (s > 0.90 && s < 1.6)
];

# # Step 3.1: Create the fault and simulation on the Dynamic Lines system

# An important aspect to consider is that DynamicLines must not be considered in the computation of the Ybus. First we construct the Dynamic Line, by finding the Line named "BUS 2-BUS 3-i_3", and then adding it to the system.

# get component return the Branch on threebus_sys_dyn named "BUS 2-BUS 3-i_3"
dyn_branch = DynamicBranch(get_component(Branch, threebus_sys_dyn, "BUS 2-BUS 3-i_3"))
# Adding a dynamic line will inmediately remove the static line from the system.
add_component!(threebus_sys_dyn, dyn_branch)

# Similarly, we construct the Ybus fault by creating a copy of the original system, but
# removing the Line "BUS 2-BUS 3-i_3" to avoid considering it in the Ybus:

#Make a copy of the original system
sys3 = deepcopy(threebus_sys);
#Remove Line "BUS 2-BUS 3-i_3"
remove_component!(Line, sys3, "BUS 2-BUS 3-i_3")
#Triplicates the impedance of the line named "BUS 1-BUS 2-i_1"
fault_branches2 = get_components(Line, sys3)
for br in fault_branches2
    if get_name(br) == "BUS 1-BUS 3-i_1"
        br.r = 3 * br.r
        br.x = 3 * br.x
        b_new = (from = br.b.from / 3, to = br.b.to / 3)
        br.b = b_new
    end
end
#Obtain the new Ybus
Ybus_fault_dyn = Ybus(sys3).data
#Define Fault: Change of YBus
Ybus_change_dyn = PowerSimulationsDynamics.NetworkSwitch(
    1.0, #change at t = 1.0
    Ybus_fault_dyn, #New YBus
)

# # Step 4.1: Run the simulation of the Dynamic Lines System

# Now, we construct the simulation:

# Time span of our simulation
tspan = (0.0, 30.0)

# Define Simulation
sim_dyn = PSID.Simulation(
    PSID.ImplicitModel, #Type of model used
    threebus_sys_dyn, #system
    pwd(), #folder to output results
    tspan, #time span
    Ybus_change_dyn, #Type of perturbation
)

# Run the simulation
PSID.execute!(
    sim_dyn, #simulation structure
    IDA(), #Sundials DAE Solver
    dtmax = 0.02, #Maximum step size
)

# We can obtain the initial conditions as:

#Will print the initial states. It also give the symbols used to describe those states.
print_device_states(sim_dyn)
#Will export a dictionary with the initial condition values to explore
x0_init_dyn = PSID.get_initial_conditions(sim_dyn)

# # Step 5.1: Store the solution

series2_dyn = get_voltage_magnitude_series(sim_dyn, 102)
zoom_dyn = [
    (series2_dyn[1][ix], series2_dyn[2][ix]) for
    (ix, s) in enumerate(series2_dyn[1]) if (s > 0.90 && s < 1.6)
];

# # Step 6.1: Compare the solutions:

# We can observe the effect of Dynamic Lines

Plots.plot(series2_dyn, label = "V_gen_dyn")
Plots.plot!(series2, label = "V_gen_st", xlabel = "Time [s]", ylabel = "Voltage [pu]")

# that looks quite similar. The differences can be observed in the zoom plot:

Plots.plot(zoom_dyn, label = "V_gen_dyn")
Plots.plot!(zoom, label = "V_gen_st", xlabel = "Time [s]", ylabel = "Voltage [pu]")


text:
using BBVI
using Test

@testset "BBVI" begin
  # TODO: Remove this test.
  @test BBVI.greet("Arthur") == "Hello, Arthur!"
end


text:
################
#### bitwise and
################

function bitand!(a::ArbUInt, b::ArbUInt)
    # writes to a
    for i in 1:min(length(a.data), length(b.data))
        a.data[i] &= b.data[i]
    end
    resize!(a.data, min(length(a.data), length(b.data)))
    normalize!(a)
end

function bitand!(a::ArbUInt, b::T) where T <: Unsigned
    if sizeof(ArbDigit) >= sizeof(T)
        resize!(a.data, 1)
        a.data[begin] &= ArbDigit(b)
    else
        resize!(a.data, sizeof(T) ÷ sizeof(ArbDigit))
        a.data .&= reinterpret(ArbDigit, [b])
    end
    a
end

function Base.:(&)(a::ArbUInt, b::ArbUInt)
    short, long = length(a.data) <= length(b.data) ? (a,b) : (b,a)
    bitand!(deepcopy(short), long)
end
function Base.:(&)(a::ArbUInt, b::T) where T <: Unsigned
    if sizeof(ArbDigit) >= sizeof(T)
        return ArbUInt(a.data[begin] & b)
    end
    n_data = a.data[begin:(sizeof(T) ÷ sizeof(ArbDigit))]
    bitand!(ArbUInt(n_data), b) # don't have to keep the upper bits!
end
Base.:(&)(a::T, b::ArbUInt) where T <: Unsigned = b & a

###############
#### bitwise or
###############

function bitor!(a::ArbUInt, b::ArbUInt)
    # writes to a
    for i in 1:min(length(a.data), length(b.data))
        a.data[i] |= b.data[i]
    end
    if length(b.data) > length(a.data)
        append!(a.data, @view(b.data[length(a.data)+1:end]))
    end
    a
end

function bitor!(a::ArbUInt, b::T) where T <: Unsigned
    # writes to a
    if sizeof(ArbDigit) >= sizeof(T)
        iszero(a) && push!(a.data, zero(ArbDigit))
        a.data[begin] |= ArbDigit(b)
    else
        mod_size = (sizeof(T)÷sizeof(ArbDigit))
        if length(a.data) < mod_size
            old_length = length(a.data)
            resize!(a.data, mod_size)
            a.data[old_length + 1:end] .= zero(ArbDigit)
        end
        a.data[begin:mod_size] .|= reinterpret(ArbDigit, [b])
    end
    normalize!(a)
end

function Base.:(|)(a::ArbUInt, b::ArbUInt)
    short, long = length(a.data) <= length(b.data) ? (a,b) : (b,a)
    bitor!(deepcopy(long), short) # gotta keep the upper bits
end
Base.:(|)(a::ArbUInt, b::T) where T <: Unsigned = bitor!(deepcopy(a), b)
Base.:(|)(a::T, b::ArbUInt) where T <: Unsigned = bitor!(deepcopy(b), a)

################
#### bitwise xor
################

function bitxor!(a::ArbUInt, b::ArbUInt)
    for i in 1:min(length(a.data), length(b.data))
        a.data[i] ⊻= b.data[i]
    end
    if length(b.data) > length(a.data)
        append!(a.data, @view(b.data[length(a.data)+1:end]))
    end
    normalize!(a)
end

function bitxor!(a::ArbUInt, b::T) where T <: Unsigned
    if sizeof(ArbDigit) >= sizeof(T)
        iszero(a) && push!(a.data, zero(ArbDigit))
        a.data[begin] ⊻= ArbDigit(b)
    else
        mod_size = (sizeof(T)÷sizeof(ArbDigit))
        if length(a.data) < mod_size
            old_length = length(a.data)
            resize!(a.data, mod_size)
            a.data[old_length + 1:end] .= zero(ArbDigit)
        end
        a.data[begin:mod_size] .⊻= reinterpret(ArbDigit, [b])
    end
    normalize!(a)
end

function Base.:(⊻)(a::ArbUInt, b::ArbUInt)
    short, long = length(a.data) <= length(b.data) ? (a,b) : (b,a)
    bitxor!(deepcopy(long), short) # gotta keep the upper bits
end
Base.:(⊻)(a::ArbUInt, b::T) where T <: Unsigned = bitxor!(deepcopy(a), b)
Base.:(⊻)(a::T, b::ArbUInt) where T <: Unsigned = bitxor!(deepcopy(b), a)

#####################
#### bitwise negation
#####################

"""
    ~(a)

Bitwise not. Does not normalize its result.
"""
function bitneg!(a::ArbUInt)
    map!(~, a.data, a.data)
    a # we mustn't normalize here, else a == ~~a doesn't hold anymore.
end
Base.:(~)(a::ArbUInt) = bitneg!(deepcopy(a))

###############
#### bitwise nor
###############

function bitnor!(a::ArbUInt, b::ArbUInt)
    # writes to a
    len_a = length(a.data)
    len_b = length(b.data)
    for i in 1:min(len_a, len_b)
        # a.data[i] ⊽= b.data[i] # FIXME: I think this should be possible?
        a.data[i] = a.data[i] ⊽ b.data[i]
    end
    if len_b != len_a
        if len_b > len_a
            append!(a.data, @view(b.data[len_a+1:end]))
        end
        a_data = @view(a.data[min(len_a,len_b)+1:end])
        map!(~, a_data, a_data)
    end
    a
end

function bitnor!(a::ArbUInt, b::T) where T <: Unsigned
    # writes to a
    if sizeof(ArbDigit) >= sizeof(T)
        iszero(a) && push!(a.data, zero(ArbDigit))
        # a.data[begin] ⊽= ArbDigit(b) # FIXME: I think this should be possible?
        a.data[begin] = a.data[begin] ⊽ ArbDigit(b)
        a_data = @view(a.data[begin+1:end])
        map!(~, a_data, a_data)
    else
        mod_size = (sizeof(T)÷sizeof(ArbDigit))
        if length(a.data) < mod_size
            old_length = length(a.data)
            resize!(a.data, mod_size)
            a.data[old_length + 1:end] .= zero(ArbDigit)
        end
        # a.data[begin:mod_size] .⊽= reinterpret(ArbDigit, [b]) # FIXME: I think this should be possible?
        a.data[begin:mod_size] .= a.data[begin:mod_size] .⊽ reinterpret(ArbDigit, [b])
        a_data = @view(a.data[mod_size+1:end])
        map!(~, a_data, a_data)
    end
    normalize!(a)
end

function Base.nor(a::ArbUInt, b::ArbUInt)
    short, long = length(a.data) <= length(b.data) ? (a,b) : (b,a)
    bitnor!(deepcopy(long), short) # gotta keep the upper bits
end
Base.nor(a::ArbUInt, b::T) where T <: Unsigned = bitnor!(deepcopy(a), b)
Base.nor(a::T, b::ArbUInt) where T <: Unsigned = bitnor!(deepcopy(b), a)

#################
#### bitwise nand
#################

function bitnand!(a::ArbUInt, b::ArbUInt)
    # writes to a
    len_a = length(a.data)
    len_b = length(b.data)
    for i in 1:min(len_a, len_b)
        # a.data[i] ⊼= b.data[i] # FIXME: I think this should be possible?
        a.data[i] = a.data[i] ⊼ b.data[i]
    end
    if len_b > len_a
        append!(a.data, @view(b.data[len_a+1:end]))
    end
    fill!(@view(a.data[min(len_a,len_b)+1:end]), typemax(ArbDigit))
    a
end

function bitnand!(a::ArbUInt, b::T) where T <: Unsigned
    # writes to a
    if sizeof(ArbDigit) >= sizeof(T)
        iszero(a) && push!(a.data, zero(ArbDigit))
        # a.data[begin] ⊼= ArbDigit(b) # FIXME: I think this should be possible?
        a.data[begin] = a.data[begin] ⊼ ArbDigit(b)
        fill!(@view(a.data[begin+1:end]), typemax(ArbDigit))
    else
        mod_size = (sizeof(T)÷sizeof(ArbDigit))
        if length(a.data) < mod_size
            old_length = length(a.data)
            resize!(a.data, mod_size)
            a.data[old_length + 1:end] .= zero(ArbDigit)
        end
        # a.data[begin:mod_size] .⊼= reinterpret(ArbDigit, [b]) # FIXME: I think this should be possible?
        a.data[begin:mod_size] .= a.data[begin:mod_size] .⊼ reinterpret(ArbDigit, [b])
        fill!(@view(a.data[mod_size+1:end]), typemax(ArbDigit))
    end
    normalize!(a)
end

function Base.nand(a::ArbUInt, b::ArbUInt)
    short, long = length(a.data) <= length(b.data) ? (a,b) : (b,a)
    bitnand!(deepcopy(long), short) # gotta keep the upper bits
end
Base.nand(a::ArbUInt, b::T) where T <: Unsigned = bitnand!(deepcopy(a), b)
Base.nand(a::T, b::ArbUInt) where T <: Unsigned = bitnand!(deepcopy(b), a)


text:
Re = 6371e3
ρc = 3996.0 * 1026.0
g  = 9.81
ρ_fw = 1000.0      # kg / m^3
L_fusion = 3.337e5 # J / kg

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


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

using Markdown
using InteractiveUtils

# ╔═╡ 42e620aa-5f4c-11eb-2ebf-85814cf720e7
begin
	using Pkg; Pkg.activate(".") # using the env in /notebooks/day2
	using PlutoUI
	using DSJulia
end

# ╔═╡ 122cffca-5fdc-11eb-3555-b39b818f1116
let
	using JuMP
	using GLPK
end

# ╔═╡ 0963185c-5fdc-11eb-0eed-89d514850353
using BioSequences

# ╔═╡ 24806108-5fdc-11eb-2f19-bb09f836f893
using Sockets

# ╔═╡ 4621c212-5fc7-11eb-2c9d-ad577506420e
using Plots

# ╔═╡ b2c1cef8-5fe5-11eb-20c7-134432196893
using Distributed

# ╔═╡ 31c1e25e-5e53-11eb-2467-9153d30962d5
md"""
# Metaprogramming

The strongest legacy of [Lisp](https://en.wikipedia.org/wiki/Lisp_(programming_language)) in the Julia language is its metaprogramming support. Like Lisp,
Julia represents its own code as a data structure of the language itself. Since code is represented
by objects that can be created and manipulated from within the language, it is possible for a
program to transform and generate its own code. This allows sophisticated code generation without
extra build steps, and also allows true Lisp-style macros operating at the level of abstract syntax trees.
Because all data types and
code in Julia are represented by Julia data structures, powerful reflection
capabilities are available to explore the internals of a program and its types just like any other
data.
"""

# ╔═╡ 4ab33c0e-5e53-11eb-2e63-2dd6f06de3ba
md"""
## Program representation

Every Julia program starts its life as a string:
"""

# ╔═╡ e76d8f04-5e53-11eb-26df-db496622642d
prog = "1 + 1"

# ╔═╡ f361d734-5e53-11eb-3957-61bd0a370b5a
md"""
**What happens next?**

The next step is to parse each string
into an object called an expression, represented by the Julia type `Expr`.
Parsing means taking the input (in this case, a string) and building a data structure – often some kind of parse tree, [abstract syntax tree](https://en.wikipedia.org/wiki/Abstract_syntax_tree) or other hierarchical structure, giving a structural representation of the input while checking for correct syntax.
"""

# ╔═╡ 15aaa5b6-5e54-11eb-067f-9dcd096940b6
ex1 = Meta.parse(prog)

# ╔═╡ 34b6a71e-5e54-11eb-0324-4bffbe85c813
typeof(ex1)

# ╔═╡ 22770a30-5e54-11eb-28db-1359ab1f402f
md"""
`Expr` objects contains two parts:

  * a `Symbol` identifying the kind of expression. A symbol is an [interned string](https://en.wikipedia.org/wiki/String_interning)
    identifier (string interning is a method of storing only one copy of each distinct string value, which must be immutable).

"""

# ╔═╡ 41f9633c-5e54-11eb-323b-e9703c674b0f
ex1.head

# ╔═╡ 48104478-5e54-11eb-2f23-c98c4ccbb763
md"""
  * the expression arguments, which may be symbols, other expressions, or literal values:
"""

# ╔═╡ 538e2054-5e54-11eb-2b9c-451637e1d8ee
ex1.args

# ╔═╡ 6e3beaf8-5e54-11eb-154b-95f5887efcb0
md"""
Expressions may also be constructed directly in [prefix notation](https://en.wikipedia.org/wiki/Polish_notation) (= Polish notation, a mathematical notation in which operators precede their operands, i.e. `+ 1 1` instead of `1 + 1` with infix notation):
"""

# ╔═╡ 7de376a8-5e54-11eb-0356-9b74225b671e
ex2 = Expr(:call, :+, 1, 1)

# ╔═╡ 825b45a6-5e54-11eb-27e6-f3cca8a2f67e
md"""
The two expressions constructed above – by parsing and by direct construction – are equivalent:
"""

# ╔═╡ 8b23cac8-5e54-11eb-2424-5972ea3d2e21
ex1 == ex2

# ╔═╡ 8ecd49ba-5e54-11eb-24b6-a328043b7a07
md"""
**The key point here is that Julia code is internally represented as a data structure that is accessible
from the language itself.**

The `dump` function provides indented and annotated display of `Expr` objects:
"""

# ╔═╡ 985a57f2-5e54-11eb-1f2f-bf23014c6171
@terminal dump(ex2)

# ╔═╡ d66ea502-5e54-11eb-2694-4766a05aed38
md"""
### Symbols

The `:` character has two syntactic purposes in Julia. The first form creates a [`Symbol`](@ref),
an [interned string](https://en.wikipedia.org/wiki/String_interning) used as one building-block
of expressions:
"""

# ╔═╡ db1e9936-5e54-11eb-376f-8ddea65d976c
s = :foo

# ╔═╡ e4f50fda-5e54-11eb-2210-6b8eb531d898
typeof(s)

# ╔═╡ ea891b6e-5e54-11eb-24c9-357a8a65d1c0
md"""
The `Symbol` constructor takes any number of arguments and creates a new symbol by concatenating
their string representations together:
"""

# ╔═╡ efe79c46-5e54-11eb-23c9-3107ea94e5b9
:foo == Symbol("foo")

# ╔═╡ f5fac50e-5e54-11eb-2ff9-cf7e30a4a8ad
Symbol("func", 10)

# ╔═╡ fe35b31e-5e54-11eb-258e-29913f1e8723
Symbol(:var, '_', "sym")

# ╔═╡ 0464f1be-5e55-11eb-34df-ef1217012663
md"""Note that to use `:` syntax, the symbol's name must be a valid identifier.
Otherwise the `Symbol(str)` constructor must be used.

In the context of an expression, symbols are used to indicate access to variables; when an expression
is evaluated, a symbol is replaced with the value bound to that symbol in the appropriate **scope**.
"""

# ╔═╡ abfc3392-5e5a-11eb-149d-ad3dd2755524
md"""
### Quoting

The second syntactic purpose of the `:` character is to create expression objects without using
the explicit `Expr` constructor. This is referred to as *quoting*. The `:` character, followed
by paired parentheses around a single statement of Julia code, produces an `Expr` object based
on the enclosed code. Here is an example of the short form used to quote an arithmetic expression:
"""

# ╔═╡ 917a579c-5e5a-11eb-3fe6-cb0dcad5e4e6
ex = :(a + b * c + 1)

# ╔═╡ 82d3220a-5e5a-11eb-3c49-d5e7afb1af2a
typeof(ex)

# ╔═╡ 67954702-5e5a-11eb-1ab1-e9a285ee208d
ex.args

# ╔═╡ 2ae6ac9c-5e5a-11eb-060c-b3eae1df0493
md"""
### Interpolation

Direct construction of `Expr` objects with value arguments is powerful, but `Expr` constructors
can be tedious compared to "normal" Julia syntax. As an alternative, Julia allows *interpolation* of
literals or expressions into quoted expressions. Interpolation is indicated by a prefix `$`.

In this example, the value of variable `a` is interpolated:
"""

# ╔═╡ 5ef3d532-5e5a-11eb-30b3-a569e0bd736e
a = 1

# ╔═╡ 3bb1f892-5e5a-11eb-0946-e5e977a4e5b2
:($a + b)

# ╔═╡ 5763cca0-5e5a-11eb-1341-3d9bfb47fb58
md"""
Interpolating into an unquoted expression is not supported and will cause a compile-time error:
"""

# ╔═╡ e6a3902a-5e58-11eb-2bd2-4bc779bcebe7
$a + b

# ╔═╡ 11991248-5e55-11eb-2748-992f6fe48620
md"""
The use of `$` for expression interpolation is intentionally reminiscent of string interpolation and command interpolation. Expression interpolation allows convenient, readable programmatic construction of complex Julia expressions.
"""

# ╔═╡ 1c016bec-5e5d-11eb-3687-931da82f8c04
md"""

### Splatting interpolation

Notice that the `$` interpolation syntax allows inserting only a single expression into an
enclosing expression.
Occasionally, you have an array of expressions and need them all to become arguments of
the surrounding expression.
This can be done with the syntax `$(args...)`.

For example, the following code generates a function call where the number of arguments is
determined programmatically:
"""

# ╔═╡ 552ea3bc-5e5d-11eb-0c97-e1e6db526df3
args = [:x, :y, :z]

# ╔═╡ 629c4f9a-5e5d-11eb-2731-bf7457755829
:(f(1, $(args...)))

# ╔═╡ 66c113d0-5e5d-11eb-0aab-b3fa7f98eae5
md"""
### Evaluating expressions

Given an expression object, one can cause Julia to evaluate (execute) it at global scope using `eval`
"""

# ╔═╡ 6fc2fb9c-5e5d-11eb-0f21-6524555a9c42
:(1 + 2)

# ╔═╡ 6dab1d58-5e5d-11eb-3381-a1796f8d14c5
eval(:(1 + 2))

# ╔═╡ b833d8ba-5e5d-11eb-3311-899efb15c1b0
md"""
## Macros

Now that we have an understanding of the basic concepts of code representation in Julia, we can introduce the core concept of this notebook: macros.
Macros provide a method to include generated code in the final body of a program. A macro maps
a tuple of arguments to a returned *expression*, and the resulting expression is compiled directly
rather than requiring a runtime `eval` call. Macro arguments may include expressions,
literal values, and symbols.

In the following examples, we will show that macros allow us to
1. modify code before it runs
2. elegantly add new features or syntax
3. process strings at compile time instead of runtime

### Basics

Here is an extraordinarily simple macro:
"""

# ╔═╡ 169c6aa2-5e5e-11eb-1d84-e7e3b4641e10
macro sayhello()
	return :( println("Hello, world!") )
end

# ╔═╡ ae1525ea-5e5e-11eb-3906-593143776559
@terminal @sayhello

# ╔═╡ c40cebe6-5e5e-11eb-1f7f-37134b7a449f
md"""
Macros have a dedicated character in Julia's syntax: the `@` (at-sign), followed by the unique
name declared in a `macro NAME ... end` block. In this example, the compiler will replace all
instances of `@sayhello` with:
"""

# ╔═╡ f5ab41f4-5e5f-11eb-02f7-efcb3e05e616
:( println("Hello, world!") )

# ╔═╡ 0001a53a-5e60-11eb-18df-d7e9e7dcf43f
md"""
When `@sayhello` is entered in the REPL, the expression executes immediately, thus we only see the evaluation result:
"""

# ╔═╡ bda7fd66-5e5e-11eb-1e06-eb67a0178e14
@terminal @sayhello()

# ╔═╡ 1707320e-5e60-11eb-3f33-798c81663488
md"""Now, consider a slightly more complex macro:"""

# ╔═╡ 1fdf9600-5e60-11eb-30d3-b167811ae47f
macro sayhello(name)
	return :( println("Hello, ", $name) )
end

# ╔═╡ 2dcacdfc-5e60-11eb-3d04-9f2ffb46d0fe
md"""
This macro takes one argument: `name`. When `@sayhello` is encountered, the quoted expression is *expanded* to interpolate the value of the argument into the final expression:
"""

# ╔═╡ 36fc5076-5e60-11eb-04be-bd9f542b1ff3
@terminal @sayhello("Mr. Bond")

# ╔═╡ 0c787198-61a0-11eb-1e30-fbcb5a624076
@terminal @sayhello "Mr. Bond"

# ╔═╡ 4a732eb8-5e60-11eb-197b-0729d0f1be15
md"""
We can view the quoted return expression using the function `macroexpand` or the macro `@macroexpand` (**important note:**
this is an extremely useful tool for debugging macros).
We can see that the `"Mr. Bond"` literal has been interpolated into the expression."""

# ╔═╡ 536addec-5e61-11eb-1a9d-d7f2c190d6c2
@macroexpand @sayhello("Mr. Bond")

# ╔═╡ 855412b2-5e61-11eb-3aba-bb862489414e
md"""

### Hold on: why macros?

We have already seen a function `f(::Expr...) -> Expr` in a previous section. In fact, `macroexpand` is also such a function. So, why do macros exist?

Macros are necessary because they execute when code is parsed, therefore, macros allow the programmer
to generate and include fragments of customized code *before* the full program is run. To illustrate
the difference, consider the following example.

Note that in this case, you need to look at your REPL to see the output. A macro definition is not allowed in the local scope, so we cannot wrap this with our `@terminal` macro.
"""

# ╔═╡ cabb3cdc-5e62-11eb-2479-bb8337c05292
macro twostep(arg)
	println("I execute at parse time. The argument is: ", arg)
	str1 = "I execute at runtime. "
	str2 = "The argument is: "
	message = str1 * str2
	return :(println($message, $arg))
end

# ╔═╡ f33f0cfe-5f52-11eb-30e0-93c4b899aaa7
md"""Note that the computation of message is compiled away if we expand the macro!"""

# ╔═╡ 01feeca0-5e63-11eb-11a0-1b66a92127d1
ex_twostep = @macroexpand @twostep :(1, 2, 3)

# ╔═╡ f1969888-5e65-11eb-18a3-879d2f87b447
@terminal dump(ex_twostep)

# ╔═╡ 0613f620-5e66-11eb-08d9-01dd3a403321
md"""
### Macro invocation

Macros are invoked with the following general syntax:

```julia
@name expr1 expr2 ...
@name(expr1, expr2, ...)
```

"""

# ╔═╡ 29a2c0ee-5e66-11eb-2194-6f1c3cef563a
md"""
Note the distinguishing `@` before the macro name and the lack of commas between the argument
expressions in the first form, and the lack of whitespace after `@name` in the second form. The
two styles should not be mixed. For example, the following syntax is different from the examples
above; it passes the tuple `(expr1, expr2, ...)` as one argument to the macro:

```julia
@name (expr1, expr2, ...)
```
"""

# ╔═╡ a1e86478-5e66-11eb-2ef9-6d460e707013
md"""
### Building an advanced macro

Here is a simplified definition of Julia's `@assert` macro, which checks if an expression is true:

(`... ? ... : ...` is the ternary if-else operator seen in `01-basics.jl` )
"""

# ╔═╡ e61381aa-5e66-11eb-3347-5d26b61e6c17
macro assert(ex)
	return :( $ex ? nothing : throw(AssertionError($(string(ex)))) )
end

# ╔═╡ a43490de-5e67-11eb-253e-f7455ade3ffe
md"""This macro can be used like this:"""

# ╔═╡ c79ad04c-5e67-11eb-1fd9-bf7ff40f0862
md"""
In place of the written syntax, the macro call is expanded at parse time to its returned result. This is equivalent to writing:"""

# ╔═╡ f34e9020-5e67-11eb-1987-e3e6a0969705
1 == 1.0 ? nothing : throw(AssertionError("1 == 1.0"))

# ╔═╡ f81955fe-5e67-11eb-11b9-5f9fd0096a09
1 == 0 ? nothing : throw(AssertionError("1 == 0"))

# ╔═╡ fb772ef4-5e67-11eb-218e-330ba087bcc8
md"""
That is, in the first call, the expression `:(1 == 1.0)` is spliced into the test condition slot,
while the value of `string(:(1 == 1.0))` is spliced into the assertion message slot. The entire
expression, thus constructed, is placed into the syntax tree where the `@assert` macro call occurs.
Then at execution time, if the test expression evaluates to true, then `nothing` is returned,
whereas if the test is false, an error is raised indicating the asserted expression that was false.
Notice that it would not be possible to write this as a function, since only the *value* of the
condition is available and it would be impossible to display the expression that computed it in
the error message.

The actual definition of `@assert` in Julia Base is more complicated. It allows the
user to optionally specify their own error message, instead of just printing the failed expression.
Just like in functions with a variable number of arguments, this is specified with an ellipses
following the last argument:
"""

# ╔═╡ 2df2baa8-5e68-11eb-2f28-9968b7ecbcd5
macro assert(ex, msgs...)
	msg_body = isempty(msgs) ? ex : msgs[1]
	msg = string(msg_body)
	return :($ex ? nothing : throw(AssertionError($msg)))
end

# ╔═╡ b0dfd0dc-5e67-11eb-19cf-5970a2e80bfa
@assert 1 == 1.0

# ╔═╡ c42d2d6a-5e67-11eb-1564-2dcb0e91fb19
@assert 1 == 0

# ╔═╡ 376fe8e4-5e68-11eb-1d7d-c9b1945b133f
md"""
Now `@assert` has two modes of operation depending upon the number of arguments it receives!
If there is only one argument, the tuple of expressions captured by `msgs` will be empty and it
will behave the same as the simpler definition above with only one argument. But now if the user specifies a second argument,
it is printed in the message body instead of the failing expression. You can inspect the result
of a macro expansion with `@macroexpand`.
"""

# ╔═╡ 4750ce04-5e68-11eb-237e-3fed9eb1f4c5
@macroexpand @assert a == b

# ╔═╡ 4dce26dc-5e68-11eb-1ff7-8974f93e3cb8
@macroexpand @assert a==b "a should equal b!"

# ╔═╡ 558920ca-5e68-11eb-3533-d5d596e3884d
md"""

There is yet another case that the actual `@assert` macro handles: what if, in addition to printing
"a should equal b," we wanted to print their values? One might naively try to use string interpolation
in the custom message, e.g., `@assert a==b "a ($a) should equal b ($b)!"`, but this won't work
as expected with the above macro. Can you see why? Recall from string interpolation that an interpolated string is rewritten to a call to `string`. Compare:
"""

# ╔═╡ 6d69106a-5e68-11eb-0d50-e17b4db5db74
typeof(:("a should equal b"))

# ╔═╡ 715464cc-5e68-11eb-3e1a-977db04a7db6
typeof(:("a ($a) should equal b ($b)!"))

# ╔═╡ 7588018e-5e68-11eb-3be2-b97ce124d6bf
@terminal dump(:("a ($a) should equal b ($b)!"))

# ╔═╡ 84771eaa-5e68-11eb-0684-cf1a3118d10b
md"""
So now instead of getting a plain string in `msg_body`, the macro is receiving a full expression that will need to be evaluated in order to display as expected. This can be spliced directly into the returned expression as an argument to the `string` call; see `error.jl` for the complete implementation.

The `@assert` macro makes great use of splicing into quoted expressions to simplify the manipulation of expressions inside the macro body.
"""

# ╔═╡ 43a19ace-5e6b-11eb-005b-b7d2d9a46878
md"""
### Macros and dispatch

Macros, just like Julia functions, are generic. This means they can also have multiple method definitions, thanks to multiple dispatch:
"""

# ╔═╡ 508167c4-5e6b-11eb-2587-b967b79cf74c
macro m end

# ╔═╡ 59f86d98-5e6b-11eb-3259-7312501b70bf
macro m(args...)
	"$(length(args)) arguments"
end

# ╔═╡ 5edf55e2-5e6b-11eb-383e-5965ae694c4c
macro m(x,y)
   	"Two arguments"
end

# ╔═╡ 3dc3aa10-5e6c-11eb-3338-7777b8cb97a5
md"""
However one should keep in mind, that macro dispatch is based on the types of AST
that are handed to the macro, not the types that the AST evaluates to at runtime:"""

# ╔═╡ 34565ba8-5e6c-11eb-2a27-2ff13394bcc8
macro m(::Int)
	"An Integer"
end

# ╔═╡ 65c94e9c-5e6b-11eb-38a8-4d60c692222b
@m "asdl"

# ╔═╡ 42fc5284-5e6c-11eb-0edb-9555a320ec55
@m 1 2

# ╔═╡ 02e4bf64-5e6d-11eb-36b5-fb57c87bff81
@m 3

# ╔═╡ 07f8cc34-5e6d-11eb-36e4-b56b2b92a2e2
x = 2

# ╔═╡ 0bbfce80-5e6d-11eb-07f6-b71086357384
@m x

# ╔═╡ 0e00bff6-5e6d-11eb-342d-bfded7ef8a96
md"""
## Code Generation

When a significant amount of repetitive boilerplate code is required, it is common to generate
it programmatically to avoid redundancy. In most languages, this requires an extra build step,
and a separate program to generate the repetitive code. In Julia, expression interpolation and `eval` allow such code generation to take place in the normal course of program execution.
For example, consider the following custom type
"""

# ╔═╡ 36805374-5e6d-11eb-205c-0524096bdf27
struct MyNumber
    x::Float64
end

# ╔═╡ 4056c9bc-5e6d-11eb-377e-238ef8633c74
md"""for which we want to add a number of methods to. We can do this programmatically in the following loop:"""

# ╔═╡ 49807d8c-5e6d-11eb-0170-f1419837d958
for op = (:sin, :cos, :tan, :log, :exp, :log)
    @eval Base.$op(a::MyNumber) = MyNumber($op(a.x))
end

# ╔═╡ 4ef3f280-5e6d-11eb-2020-b7cd9e6b9c03
md"""and we can now use those functions with our custom type:"""

# ╔═╡ 53f36f36-5e6d-11eb-3f06-6de0eff4d173
y = MyNumber(π)

# ╔═╡ 3428d7b8-5e6d-11eb-32bc-af1df0579e60
sin(y)

# ╔═╡ 18e3753a-5f55-11eb-0a11-47ca6abe9186
md"This will not work since we only have defined `log` and not `log10`."

# ╔═╡ 2f4fe6ce-5e6c-11eb-1744-2f9534c8b6ba
log10(y)

# ╔═╡ cae2546a-5e6b-11eb-3fc2-69f9a15107cf
md"""In this manner, Julia acts as its own preprocessor, and allows code generation from inside the language."""

# ╔═╡ 0ac5096e-61a1-11eb-24ba-13bc0a3377c6

# ╔═╡ efc1357a-61a0-11eb-20f5-a15338080e4c
md"### Example: domain specific languages"

# ╔═╡ 4d9de330-5f55-11eb-0d89-eb81e6c9ffab
md"
Code generation can for instance be used to simplify the creation of a mathematical optimisation problem. In this case, we we'll use the `JuMP` package.
`JuMP` ('Julia for Mathematical Programming') is an open-source modeling language that is embedded in Julia. It allows users to formulate various classes of optimization problems (linear, mixed-integer, quadratic, conic quadratic, semidefinite, and nonlinear) with easy-to-read code. `JuMP` also makes advanced optimization techniques easily accessible from a high-level language.

As a dummy example, let us consider the following linear programming problem:

``\max_{x,y}\,\,x + 2y``

``\text{s.t.}``

``x + y \leq 1``

``0\leq x, y \leq 1``

Which can be transcribed into a `JuMP` model as:
"

# ╔═╡ 92458f54-5f57-11eb-179f-c749ce06f7e2
let
	model = Model(GLPK.Optimizer)
	@variable(model, 0 <= x <= 1)
	@variable(model, 0 <= y <= 1)
	@constraint(model, x + y <= 1)
	@objective(model, Max, x + 2y)
	optimize!(model)
	value(x), value(y), objective_value(model)
end

# ╔═╡ d1222dcc-5f57-11eb-2c87-77a16ea8e65d
md"
Without the macros, the code would be more difficult to read. As an example, check the `macroexpansion` of `@constraint` to see the bunch of code that is generated behind the scenes:
"

# ╔═╡ 0b1f3542-5f58-11eb-1610-75df3e2a4a76
let
	model = Model(GLPK.Optimizer)
	@macroexpand @constraint(model, x + y <= 1)
end

# ╔═╡ 57f243ec-5fe3-11eb-04d9-d776ccba4c0f
md"Finally, this section on code generation ends with two more examples:
1. a macro to repeat a certain expression n-times
2. a macro to repeat a certain expression until a condition is met

Are these macros the best way to tackle the problems at hand? Maybe not, but they do give a nice illustration of how the code generation works.
"

# ╔═╡ 211cf9ca-5fe3-11eb-2596-71b1848753d7
macro dotimes(n, body)
    quote
        for i = 1:$(esc(n))
            $(esc(body))
        end
    end
end

# ╔═╡ 35ccaee2-5fe3-11eb-3bd0-e355b0fd3f71
@macroexpand @dotimes 3 println("hi there")

# ╔═╡ 434669d2-5fe3-11eb-1f4a-27edf1f9e7f6
macro until(condition, block)
    quote
        while true
            $(esc(block))
            if $(esc(condition))
                break
            end
        end
    end
end

# ╔═╡ 9612965e-5fe3-11eb-172a-59cc73a26ab6
PlutoUI.with_terminal() do
	i = 0
	@until i == 10 begin
		i += 1
		println(i)
	end
end

# ╔═╡ 3c2cce20-5fe3-11eb-32e3-f542f9a94a6e
@macroexpand @until j == 10 begin
	j += 1
	println(j)
end

# ╔═╡ 10a8532a-5e60-11eb-3331-974e708cb39d
md"""
## Non-Standard String Literals

String literals prefixed by an identifier are called non-standard string literals, and can have different semantics than un-prefixed string literals. For example:

  * `r"^\s*(?:#|$)"` produces a regular expression object rather than a string
  * `b"DATA\xff\u2200"` is a byte array literal for `[68,65,84,65,255,226,136,128]`.

Perhaps surprisingly, these behaviors are not hard-coded into the Julia parser or compiler. Instead,
they are custom behaviors provided by a general mechanism that anyone can use: prefixed string
literals are parsed as calls to specially-named macros. For example, the regular expression macro
is just the following:
```julia
macro r_str(p)
    Regex(p)
end
```
"""

# ╔═╡ d8a9ad46-5e7a-11eb-1cce-3d4bc49cd332
md"""
That's all. This macro says that the literal contents of the string literal `r"^\s*(?:#|$)"` should
be passed to the `@r_str` macro and the result of that expansion should be placed in the syntax
tree where the string literal occurs. In other words, the expression `r"^\s*(?:#|$)"` is equivalent
to placing the following object directly into the syntax tree:
"""

# ╔═╡ e85b2d3c-5e7a-11eb-3b89-e11bc274f7cf
Regex("^\\s*(?:#|\$)")

# ╔═╡ ec73cd70-5e7a-11eb-0285-6b1fabd1289d
md"""
Not only is the string literal form shorter and far more convenient, but it is also more efficient:
since the regular expression is compiled, which takes time, and the `Regex` object is actually created *when the code is compiled*,
the compilation occurs only once, rather than every time the code is executed. Consider if the
regular expression occurs in a loop:
"""

# ╔═╡ f1d24ba2-5e7a-11eb-010a-f543cfc306b5
PlutoUI.with_terminal() do
	for line ∈ ["first", "sec#nd", "third"]
		m = match(r"#", line)
		if m === nothing
			println("nothing found")
		else
			println("found something")
		end
	end
end

# ╔═╡ f838762c-5e7a-11eb-000a-e963fab5e97d
md"""
Since the regular expression `r"^\s*(?:#|$)"` is compiled and inserted into the syntax tree when
this code is parsed, the expression is only compiled once instead of each time the loop is executed.
In order to accomplish this without macros, one would have to write this loop like this:
"""

# ╔═╡ 149c1c12-5e7b-11eb-297d-6da57f45891b
PlutoUI.with_terminal() do
	re = Regex("#")
	for line ∈ ["first", "sec#nd", "third"]
		m = match(re, line)
		if m === nothing
			println("nothing found")
		else
			println("found something")
		end
	end
end

# ╔═╡ 780da60e-5e7c-11eb-0c1a-55be73da188f
md"""
Moreover, if the compiler could not determine that the regex object was constant over all loops,
certain optimizations might not be possible, making this version still less efficient than the
more convenient literal form above. Of course, there are still situations where the non-literal
form is more convenient: if one needs to interpolate a variable into the regular expression, one
must take this more verbose approach. In the vast majority of use cases, however, regular expressions
are not constructed based on run-time data. In this majority of cases, the ability to write regular
expressions as compile-time values is invaluable.
"""

# ╔═╡ 90a672d0-5f5c-11eb-3a62-1b034cd41e67
md"
Designing your own custom string literals can be done as such:
"

# ╔═╡ b4104d12-5e7c-11eb-3e22-a78b74526129
macro foo_str(str, flag)
    # do stuff
	str, flag
end

# ╔═╡ 20393e86-5e7d-11eb-18e3-613890472903
foo"this is the string"theflag

# ╔═╡ 5069df16-5e7d-11eb-217d-6b740e9b3559
md"""
The first example of a custom macro can be found in every notebook in this course! It is the markdown string literal, which allows the usage of Markdown markup language to prettify these lectures!

```julia
md"I am a Markdown string with glorious **formatting capabilities**."
```
I am a Markdown string with glorious *formatting* **capabilities**.
"""

# ╔═╡ 1a50ff0e-5e7d-11eb-2fc4-cd5c12015751
md"The second example is from the `BioSequences.jl` package. In that bioinformatics package, you can define sequences e.g. DNA and RNA as string literals:
"

# ╔═╡ 0664eb22-5e7d-11eb-07a6-ed35143bf03f
dna"ACGT"

# ╔═╡ a965bd70-5e7c-11eb-13dd-5fe83950af11
md"Repetition and concatenation can be performed pretty straightforward:"

# ╔═╡ fb0c59fa-5f61-11eb-1ab7-bf7b80746aa2
repeat(dna"TTAGGG", 10)

# ╔═╡ a157dafc-5e7c-11eb-105d-4db3ee9dbec6
dna = dna"ACGT" * dna"TGCAA"

# ╔═╡ 0ec9887c-5f5e-11eb-29e5-931881dcb222
md"Other typical string operations such as pushing new values work as you would expect"

# ╔═╡ 1d4aca5a-5f5e-11eb-23c3-d1ddef4e6ee6
push!(dna, DNA_A)

# ╔═╡ c4c2137a-5f5d-11eb-1e06-8f5e2b1f1c0e
md"There exist methods to convert a DNA sequence to its RNA equivalent:"

# ╔═╡ e13fcbbe-5f5d-11eb-2074-fdbc945861e1
rna = convert(LongRNASeq, dna)

# ╔═╡ 2c1308f4-5f5e-11eb-058f-79cbf7dd23c6
md"Note that altough the printout of the DNA and RNA object is different because of different nucleotides. The information content is the same, and as such this statement is true:"

# ╔═╡ fdf46c92-5f5d-11eb-345b-e145c9c2874c
dna.data === rna.data

# ╔═╡ 37da420a-5f61-11eb-3688-09e6188165e3
md"Remember before that we mentioned that you can add a flag to a string literal? In `BioSequences.jl` this has a use case.

If you have a function that generates a sequence, and you want it to create a new sequence each time it is called, then you can add a flag to the end of the sequence literal to dictate behaviour: A flag of 's' means 'static': the sequence will be allocated before code is run, as is the default behaviour. However providing 'd' flag changes the behaviour: 'd' means 'dynamic': the sequence will be allocated whilst the code is running, and not before. So to change foo so as it creates a new sequence each time it is called, simply add the 'd' flag to the sequence literal:"

# ╔═╡ d9f18616-5e8d-11eb-1b1a-3bcc749fb467
function getdna_dynamic()
	s = dna"CTT"d     # 'd' flag appended to the string literal.
	push!(s, DNA_A)
end

# ╔═╡ ea3da4dc-5e8d-11eb-0910-091b926aea58
getdna_dynamic()

# ╔═╡ f13dd22a-5e8d-11eb-1def-137de66123ba
getdna_dynamic()

# ╔═╡ be11ec9a-5f61-11eb-374a-c96489ea582d
md"Output of `getdna_dynamic()` stays the same!"

# ╔═╡ f44d757e-5e8d-11eb-18f1-719d3c7b8688
function getdna_static()
	s = dna"CTT"s     # 's' flag appended to the string literal.
	push!(s, DNA_A)
end

# ╔═╡ fa310fdc-5e8d-11eb-024b-7ded3e213112
getdna_static()

# ╔═╡ 1acb02d4-5e8e-11eb-232c-f756e92c3f97
getdna_static()

# ╔═╡ 2dcb4600-5e8e-11eb-093e-59eb79f0cf20
md"""
Be careful when you are using sequence literals inside of functions, and inside the bodies of things like for loops. And if you use them and are unsure, use the 's' and 'd' flags to ensure the behaviour you get is the behaviour you intend.
"""

# ╔═╡ 2e2601b6-5e94-11eb-3613-eb5fee19b6b7
md"
## Overview of some interesting macros
"

# ╔═╡ 311a4666-5fc5-11eb-2cc0-ef66e5be96e3
md"check if an expression is true"

# ╔═╡ 4622422a-5fc5-11eb-20f2-9bfe466f8f30
@assert true == true

# ╔═╡ 5ba5d6b6-5fc5-11eb-256b-73a954a5db68
md"Integers and floating point numbers with arbitrary precision. This macro exists because promoting a floating point number to a `BigFloat` will keep the approximation error of `Float64`."

# ╔═╡ a47ab60e-5fc5-11eb-368a-d76fc6ee640d
big"0.1"

# ╔═╡ a99546c2-5fc5-11eb-0902-71cb283543c8
@big_str "0.1"

# ╔═╡ af2c7d4e-5fc5-11eb-30d4-a367d2c0db23
big(0.1)

# ╔═╡ 4a335466-5fc6-11eb-2273-312fd3928616
md"There exists a lot more unique string literals than we have shown here, such as html strings, ip address literals, etc."

# ╔═╡ 17f771dc-5fc6-11eb-3b83-ed51f87194e0
html"""<!-- HTML generated using hilite.me --><div style="background: #ffffff; overflow:auto;width:auto;border:solid gray;border-width:.1em .1em .1em .8em;padding:.2em .6em;"><table><tr><td><pre style="margin: 0; line-height: 125%">1</pre></td><td><pre style="margin: 0; line-height: 125%">println(<span style="color: #a31515">"hello world!"</span>)
</pre></td></tr></table></div>
"""

# ╔═╡ 92836b48-5fc6-11eb-19b6-1d7a5c310052
@terminal println("There is nothing like ", ip"127.0.0.1")

# ╔═╡ f26e65a8-5fc6-11eb-3cd7-3f40c12a6abb
md"Find out what function is exactly used in multiple dipatching."

# ╔═╡ 154e2392-5fc7-11eb-250b-e726ab7710e7
@which sin(2.2)

# ╔═╡ 2cfa3952-5fc7-11eb-2e79-95f86aafbf2c
@which sin(2)

# ╔═╡ 4aece42a-5fc7-11eb-291f-99120483339f
md"
There exist some handy macros to analyse execution time and memory allocation.

```julia
using BenchmarkTools
@btime sin(2)
@benchmark sin(2)
@elapsed sin(2)
```
"

# ╔═╡ 295e0028-5fd6-11eb-39e4-b392de1c0de2
md"""
When developing modules, scripts or julia packages, you can use the `@info`, `@warn`, `@error` and `@debug` as logging macros. They are mostly useful in packages, not in notebooks like this.

An example usage would be a warning thrown by an optimisation algorithm to tell you that e.g. the predefined accuracy or tolerance was not reached.

Note, to make `@debug` you need to set an environment variable to mark that you are in debug mode: `ENV["JULIA_DEBUG"] = "all"`
"""

# ╔═╡ 4a9d85ba-5fc7-11eb-1b6e-59ef2259dd24
PlutoUI.with_terminal() do
	@info "Information comes here"
	@error "Error has been found at this exact location"
	@warn "Same. but for a warning"
end

# ╔═╡ 34829ada-5fdc-11eb-263a-cf60a9650556
md"Plot recipes are a nifty thing to make plots for your custom data type:"

# ╔═╡ 5c33798c-5fda-11eb-186f-79450fb59e3b
struct TemperatureMeas
	t1
	t2
	t3
	TemperatureMeas(t1, t2, t3) = new(t1, t2, t3)
end

# ╔═╡ 1ca80d42-5fdb-11eb-153f-9ddb25a4fb6b
@recipe function f(tempmeas::TemperatureMeas)
    xguide --> "time (seconds)"
    yguide --> "temperature (Celsius)"
    [tempmeas.t1, tempmeas.t2, tempmeas.t3] # return the arguments (input data) for the next recipe
end

# ╔═╡ 1a793992-5fdb-11eb-3b38-05544f268a65
plot(TemperatureMeas(Plots.fakedata(50), Plots.fakedata(50), Plots.fakedata(50)))

# ╔═╡ d418dca4-5fe5-11eb-3f1a-a31e5dc5f7e9
md"
Using the `Distributed` module, one can easily transform certain types of code to run distributed (i.e. on multiple cores). Below is an example of how you can transform a for-loop into its distributed version.
"

# ╔═╡ f1b70998-5fdb-11eb-292b-0776f1b03816
PlutoUI.with_terminal() do
	@sync @distributed for i ∈ 1:5
		println(i)
	end
end

# ╔═╡ Cell order:
# ╠═42e620aa-5f4c-11eb-2ebf-85814cf720e7
# ╟─31c1e25e-5e53-11eb-2467-9153d30962d5
# ╟─4ab33c0e-5e53-11eb-2e63-2dd6f06de3ba
# ╠═e76d8f04-5e53-11eb-26df-db496622642d
# ╟─f361d734-5e53-11eb-3957-61bd0a370b5a
# ╠═15aaa5b6-5e54-11eb-067f-9dcd096940b6
# ╠═34b6a71e-5e54-11eb-0324-4bffbe85c813
# ╟─22770a30-5e54-11eb-28db-1359ab1f402f
# ╠═41f9633c-5e54-11eb-323b-e9703c674b0f
# ╟─48104478-5e54-11eb-2f23-c98c4ccbb763
# ╠═538e2054-5e54-11eb-2b9c-451637e1d8ee
# ╟─6e3beaf8-5e54-11eb-154b-95f5887efcb0
# ╠═7de376a8-5e54-11eb-0356-9b74225b671e
# ╟─825b45a6-5e54-11eb-27e6-f3cca8a2f67e
# ╠═8b23cac8-5e54-11eb-2424-5972ea3d2e21
# ╟─8ecd49ba-5e54-11eb-24b6-a328043b7a07
# ╠═985a57f2-5e54-11eb-1f2f-bf23014c6171
# ╟─d66ea502-5e54-11eb-2694-4766a05aed38
# ╠═db1e9936-5e54-11eb-376f-8ddea65d976c
# ╠═e4f50fda-5e54-11eb-2210-6b8eb531d898
# ╟─ea891b6e-5e54-11eb-24c9-357a8a65d1c0
# ╠═efe79c46-5e54-11eb-23c9-3107ea94e5b9
# ╠═f5fac50e-5e54-11eb-2ff9-cf7e30a4a8ad
# ╠═fe35b31e-5e54-11eb-258e-29913f1e8723
# ╟─0464f1be-5e55-11eb-34df-ef1217012663
# ╟─abfc3392-5e5a-11eb-149d-ad3dd2755524
# ╠═917a579c-5e5a-11eb-3fe6-cb0dcad5e4e6
# ╠═82d3220a-5e5a-11eb-3c49-d5e7afb1af2a
# ╠═67954702-5e5a-11eb-1ab1-e9a285ee208d
# ╟─2ae6ac9c-5e5a-11eb-060c-b3eae1df0493
# ╠═5ef3d532-5e5a-11eb-30b3-a569e0bd736e
# ╠═3bb1f892-5e5a-11eb-0946-e5e977a4e5b2
# ╟─5763cca0-5e5a-11eb-1341-3d9bfb47fb58
# ╠═e6a3902a-5e58-11eb-2bd2-4bc779bcebe7
# ╟─11991248-5e55-11eb-2748-992f6fe48620
# ╟─1c016bec-5e5d-11eb-3687-931da82f8c04
# ╠═552ea3bc-5e5d-11eb-0c97-e1e6db526df3
# ╠═629c4f9a-5e5d-11eb-2731-bf7457755829
# ╟─66c113d0-5e5d-11eb-0aab-b3fa7f98eae5
# ╠═6fc2fb9c-5e5d-11eb-0f21-6524555a9c42
# ╠═6dab1d58-5e5d-11eb-3381-a1796f8d14c5
# ╟─b833d8ba-5e5d-11eb-3311-899efb15c1b0
# ╠═169c6aa2-5e5e-11eb-1d84-e7e3b4641e10
# ╠═ae1525ea-5e5e-11eb-3906-593143776559
# ╟─c40cebe6-5e5e-11eb-1f7f-37134b7a449f
# ╠═f5ab41f4-5e5f-11eb-02f7-efcb3e05e616
# ╟─0001a53a-5e60-11eb-18df-d7e9e7dcf43f
# ╠═bda7fd66-5e5e-11eb-1e06-eb67a0178e14
# ╟─1707320e-5e60-11eb-3f33-798c81663488
# ╠═1fdf9600-5e60-11eb-30d3-b167811ae47f
# ╟─2dcacdfc-5e60-11eb-3d04-9f2ffb46d0fe
# ╠═36fc5076-5e60-11eb-04be-bd9f542b1ff3
# ╠═0c787198-61a0-11eb-1e30-fbcb5a624076
# ╟─4a732eb8-5e60-11eb-197b-0729d0f1be15
# ╠═536addec-5e61-11eb-1a9d-d7f2c190d6c2
# ╟─855412b2-5e61-11eb-3aba-bb862489414e
# ╠═cabb3cdc-5e62-11eb-2479-bb8337c05292
# ╟─f33f0cfe-5f52-11eb-30e0-93c4b899aaa7
# ╠═01feeca0-5e63-11eb-11a0-1b66a92127d1
# ╠═f1969888-5e65-11eb-18a3-879d2f87b447
# ╟─0613f620-5e66-11eb-08d9-01dd3a403321
# ╟─29a2c0ee-5e66-11eb-2194-6f1c3cef563a
# ╟─a1e86478-5e66-11eb-2ef9-6d460e707013
# ╠═e61381aa-5e66-11eb-3347-5d26b61e6c17
# ╟─a43490de-5e67-11eb-253e-f7455ade3ffe
# ╠═b0dfd0dc-5e67-11eb-19cf-5970a2e80bfa
# ╠═c42d2d6a-5e67-11eb-1564-2dcb0e91fb19
# ╟─c79ad04c-5e67-11eb-1fd9-bf7ff40f0862
# ╠═f34e9020-5e67-11eb-1987-e3e6a0969705
# ╠═f81955fe-5e67-11eb-11b9-5f9fd0096a09
# ╟─fb772ef4-5e67-11eb-218e-330ba087bcc8
# ╠═2df2baa8-5e68-11eb-2f28-9968b7ecbcd5
# ╟─376fe8e4-5e68-11eb-1d7d-c9b1945b133f
# ╠═4750ce04-5e68-11eb-237e-3fed9eb1f4c5
# ╠═4dce26dc-5e68-11eb-1ff7-8974f93e3cb8
# ╟─558920ca-5e68-11eb-3533-d5d596e3884d
# ╠═6d69106a-5e68-11eb-0d50-e17b4db5db74
# ╠═715464cc-5e68-11eb-3e1a-977db04a7db6
# ╠═7588018e-5e68-11eb-3be2-b97ce124d6bf
# ╟─84771eaa-5e68-11eb-0684-cf1a3118d10b
# ╟─43a19ace-5e6b-11eb-005b-b7d2d9a46878
# ╠═508167c4-5e6b-11eb-2587-b967b79cf74c
# ╠═59f86d98-5e6b-11eb-3259-7312501b70bf
# ╠═5edf55e2-5e6b-11eb-383e-5965ae694c4c
# ╠═65c94e9c-5e6b-11eb-38a8-4d60c692222b
# ╠═42fc5284-5e6c-11eb-0edb-9555a320ec55
# ╟─3dc3aa10-5e6c-11eb-3338-7777b8cb97a5
# ╠═34565ba8-5e6c-11eb-2a27-2ff13394bcc8
# ╠═02e4bf64-5e6d-11eb-36b5-fb57c87bff81
# ╠═07f8cc34-5e6d-11eb-36e4-b56b2b92a2e2
# ╠═0bbfce80-5e6d-11eb-07f6-b71086357384
# ╟─0e00bff6-5e6d-11eb-342d-bfded7ef8a96
# ╠═36805374-5e6d-11eb-205c-0524096bdf27
# ╟─4056c9bc-5e6d-11eb-377e-238ef8633c74
# ╠═49807d8c-5e6d-11eb-0170-f1419837d958
# ╟─4ef3f280-5e6d-11eb-2020-b7cd9e6b9c03
# ╠═53f36f36-5e6d-11eb-3f06-6de0eff4d173
# ╠═3428d7b8-5e6d-11eb-32bc-af1df0579e60
# ╟─18e3753a-5f55-11eb-0a11-47ca6abe9186
# ╠═2f4fe6ce-5e6c-11eb-1744-2f9534c8b6ba
# ╟─cae2546a-5e6b-11eb-3fc2-69f9a15107cf
# ╟─0ac5096e-61a1-11eb-24ba-13bc0a3377c6
# ╟─efc1357a-61a0-11eb-20f5-a15338080e4c
# ╟─4d9de330-5f55-11eb-0d89-eb81e6c9ffab
# ╠═122cffca-5fdc-11eb-3555-b39b818f1116
# ╠═92458f54-5f57-11eb-179f-c749ce06f7e2
# ╟─d1222dcc-5f57-11eb-2c87-77a16ea8e65d
# ╠═0b1f3542-5f58-11eb-1610-75df3e2a4a76
# ╟─57f243ec-5fe3-11eb-04d9-d776ccba4c0f
# ╠═211cf9ca-5fe3-11eb-2596-71b1848753d7
# ╠═35ccaee2-5fe3-11eb-3bd0-e355b0fd3f71
# ╠═434669d2-5fe3-11eb-1f4a-27edf1f9e7f6
# ╠═9612965e-5fe3-11eb-172a-59cc73a26ab6
# ╠═3c2cce20-5fe3-11eb-32e3-f542f9a94a6e
# ╟─10a8532a-5e60-11eb-3331-974e708cb39d
# ╟─d8a9ad46-5e7a-11eb-1cce-3d4bc49cd332
# ╠═e85b2d3c-5e7a-11eb-3b89-e11bc274f7cf
# ╟─ec73cd70-5e7a-11eb-0285-6b1fabd1289d
# ╠═f1d24ba2-5e7a-11eb-010a-f543cfc306b5
# ╟─f838762c-5e7a-11eb-000a-e963fab5e97d
# ╠═149c1c12-5e7b-11eb-297d-6da57f45891b
# ╟─780da60e-5e7c-11eb-0c1a-55be73da188f
# ╟─90a672d0-5f5c-11eb-3a62-1b034cd41e67
# ╠═b4104d12-5e7c-11eb-3e22-a78b74526129
# ╠═20393e86-5e7d-11eb-18e3-613890472903
# ╟─5069df16-5e7d-11eb-217d-6b740e9b3559
# ╟─1a50ff0e-5e7d-11eb-2fc4-cd5c12015751
# ╠═0963185c-5fdc-11eb-0eed-89d514850353
# ╠═0664eb22-5e7d-11eb-07a6-ed35143bf03f
# ╟─a965bd70-5e7c-11eb-13dd-5fe83950af11
# ╠═fb0c59fa-5f61-11eb-1ab7-bf7b80746aa2
# ╠═a157dafc-5e7c-11eb-105d-4db3ee9dbec6
# ╟─0ec9887c-5f5e-11eb-29e5-931881dcb222
# ╠═1d4aca5a-5f5e-11eb-23c3-d1ddef4e6ee6
# ╟─c4c2137a-5f5d-11eb-1e06-8f5e2b1f1c0e
# ╠═e13fcbbe-5f5d-11eb-2074-fdbc945861e1
# ╟─2c1308f4-5f5e-11eb-058f-79cbf7dd23c6
# ╠═fdf46c92-5f5d-11eb-345b-e145c9c2874c
# ╟─37da420a-5f61-11eb-3688-09e6188165e3
# ╠═d9f18616-5e8d-11eb-1b1a-3bcc749fb467
# ╠═ea3da4dc-5e8d-11eb-0910-091b926aea58
# ╠═f13dd22a-5e8d-11eb-1def-137de66123ba
# ╟─be11ec9a-5f61-11eb-374a-c96489ea582d
# ╠═f44d757e-5e8d-11eb-18f1-719d3c7b8688
# ╠═fa310fdc-5e8d-11eb-024b-7ded3e213112
# ╠═1acb02d4-5e8e-11eb-232c-f756e92c3f97
# ╟─2dcb4600-5e8e-11eb-093e-59eb79f0cf20
# ╟─2e2601b6-5e94-11eb-3613-eb5fee19b6b7
# ╟─311a4666-5fc5-11eb-2cc0-ef66e5be96e3
# ╠═4622422a-5fc5-11eb-20f2-9bfe466f8f30
# ╟─5ba5d6b6-5fc5-11eb-256b-73a954a5db68
# ╠═a47ab60e-5fc5-11eb-368a-d76fc6ee640d
# ╠═a99546c2-5fc5-11eb-0902-71cb283543c8
# ╠═af2c7d4e-5fc5-11eb-30d4-a367d2c0db23
# ╟─4a335466-5fc6-11eb-2273-312fd3928616
# ╟─17f771dc-5fc6-11eb-3b83-ed51f87194e0
# ╠═24806108-5fdc-11eb-2f19-bb09f836f893
# ╠═92836b48-5fc6-11eb-19b6-1d7a5c310052
# ╟─f26e65a8-5fc6-11eb-3cd7-3f40c12a6abb
# ╠═154e2392-5fc7-11eb-250b-e726ab7710e7
# ╠═2cfa3952-5fc7-11eb-2e79-95f86aafbf2c
# ╟─4aece42a-5fc7-11eb-291f-99120483339f
# ╟─295e0028-5fd6-11eb-39e4-b392de1c0de2
# ╠═4a9d85ba-5fc7-11eb-1b6e-59ef2259dd24
# ╟─34829ada-5fdc-11eb-263a-cf60a9650556
# ╠═4621c212-5fc7-11eb-2c9d-ad577506420e
# ╠═5c33798c-5fda-11eb-186f-79450fb59e3b
# ╠═1ca80d42-5fdb-11eb-153f-9ddb25a4fb6b
# ╠═1a793992-5fdb-11eb-3b38-05544f268a65
# ╟─d418dca4-5fe5-11eb-3f1a-a31e5dc5f7e9
# ╠═b2c1cef8-5fe5-11eb-20c7-134432196893
# ╠═f1b70998-5fdb-11eb-292b-0776f1b03816


text:
printstyled("Tests complete. Cleaning up...\n", color=:light_green)
flush(out)
close(out)
for fpat in ("*.mseed", "*.SAC", "*.sac", "*.geocsv", "*.xml", "*.h5", "FDSNevq.log")
  try
    files = ls(fpat)
    for f in files
      rm(f)
    end
  catch err
    println("Attempting to delete ", fpat, " threw error: ", err)
  end
end


text:
## get data x, y, property
using DelimitedFiles, PyPlot
d, h = readdlm("mohodepth_GDA_LCC.csv",',', header=true)
x, y = d[:,5], d[:,6]
confidence = d[:,3]
confidence[confidence.<0.1].=0.1
σ = 1 ./confidence
z = d[:,4]
f = figure()
scatter(x, y, c=z, cmap="RdBu")
xlabel("Easting m")
ylabel("Northing m")
cb = colorbar()
cb.ax.set_xlabel("km")
gca().set_aspect(1.)
## make a grid
using NearestNeighbors
xmin, xmax = extrema(x)
ymin, ymax = extrema(y)
delx, dely = 20_000., 20_000
xgrid, ygrid = [x for y in ymin:dely:ymax, x in xmin:delx:xmax], [y for y in ymin:dely:ymax, x in xmin:delx:xmax]
xall = [xgrid[:]'; ygrid[:]']
kdtree = KDTree(xall)
idxs, = nn(kdtree, [x[:]';y[:]'])
imagegrid = NaN .+ zeros(size(xgrid))
noisegrid = copy(imagegrid)
imagegrid[idxs] .= z
noisegrid[idxs] .= σ
figure(figsize=(12,4))
s1 = subplot(121)
imshow(imagegrid, extent=[xmin, xmax, ymin, ymax], origin="lower", aspect=1., cmap="RdBu")
title("moho depth")
xlabel("Easting m")
ylabel("Northing m")
cb = colorbar()
cb.ax.set_xlabel("km")
subplot(122, sharex=s1, sharey=s1)
imshow(noisegrid, extent=[xmin, xmax, ymin, ymax], origin="lower", aspect=1., cmap="viridis")
title("noise on depth")
xlabel("Easting m")
ylabel("Northing m")
cb = colorbar()
cb.ax.set_xlabel("km")
plt.tight_layout()

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
@testset "Test Dummy" begin
  @test 1 == 1
end


text:
ys = cu([3. 3. 4. 4. 5.;
         5. 5. 6. 6. 7.])
us = cu(2*ones(2, 3, 4))
xs = CuArray{Int64}([1 2 3 4;
                     4 2 1 3;
                     3 5 5 3])

∇y_mul = [4. 4. 16. 4. 4.; 4. 4. 16. 4. 4.]
∇y_div = [.25 .25 .0625 .25 .25; .25 .25 .0625 .25 .25]
∇u_mean = cat([.5 .5 .25; .5 .5 .25], [.5 .5 .5; .5 .5 .5],
              [.25 .5 .5; .25 .5 .5], [.5 .25 .25; .5 .25 .25], dims=3)

@testset "cuda/grad" begin
    @testset "pool" begin
        @test Zygote.gradient(x -> sum(sumpool(x, us)), xs) == (nothing,)
        @test Zygote.gradient(x -> sum(sumpool(xs, x)), us) == (ones(2, 3, 4),)

        @test Zygote.gradient(x -> sum(subpool(x, us)), xs) == (nothing,)
        @test Zygote.gradient(x -> sum(subpool(xs, x)), us) == (-ones(2, 3, 4),)

        @test Zygote.gradient(x -> sum(maxpool(x, us)), xs) == (nothing,)
        @test Zygote.gradient(x -> sum(maxpool(xs, x)), us) == (ones(2, 3, 4),)

        @test Zygote.gradient(x -> sum(minpool(x, us)), xs) == (nothing,)
        @test Zygote.gradient(x -> sum(minpool(xs, x)), us) == (ones(2, 3, 4),)

        @test Zygote.gradient(x -> sum(prodpool(x, us)), xs) == (nothing,)
        @test Zygote.gradient(x -> sum(prodpool(xs, x)), us) == (2048*ones(2, 3, 4),)

        @test Zygote.gradient(x -> sum(divpool(x, us)), xs) == (nothing,)
        @test Zygote.gradient(x -> sum(divpool(xs, x)), us) == (-ones(2, 3, 4)/8192,)

        @test Zygote.gradient(x -> sum(meanpool(x, us)), xs) == (nothing,)
        @test Zygote.gradient(x -> sum(meanpool(xs, x)), us) == (∇u_mean,)
    end
end


text:
"""
Base routines used in higher level algorithms implemeted using the
TensorOperations package

product of site tensor A with local operator O and right environment R
"""
function prod_AOR(A::AbstractArray{<:Number, 3}, O::AbstractArray{<:Number, 4},
    R::AbstractArray{<:Number, 3})

    a1, a2, a3 = size(A)
    _, r2, r3 = size(R)
    o1 = size(O, 1)

    AOR = reshape(A, a1*a2, a3)*reshape(R, a3, r2*r3)
    AOR = permutedims(reshape(AOR, a1, a2, r2, r3), [3, 2, 1, 4])
    AOR = reshape(O, o1*a2, r2*a2)*reshape(AOR, r2*a2, a1*r3)
    permutedims(reshape(AOR, o1, a2, a1, r3), [3, 1, 2, 4])

end

"""
product of site tensor A with local operator O and left environment L
"""

function prod_AOL(A::AbstractArray{<:Number, 3}, O::AbstractArray{<:Number, 4},
    L::AbstractArray{<:Number, 3})

    a1, d, a3 = size(A)
    _, l2, l3 = size(L)
    o3 = size(O, 3)

    AOL = reshape(L, a1, l2*l3)'*reshape(A, a1, d*a3)
    AOL = permutedims(reshape(AOL, l2, l3, d, a3), [2, 4, 1, 3])
    AOL = reshape(AOL, l3*a3, l2*d)
    AOL = AOL*reshape(permutedims(O, [1, 4, 2, 3]), l2*d, d*o3)
    permutedims(reshape(AOL, l3, a3, d, o3), [1, 3, 2, 4])

end

"""
product of left environment and AOR from prod_AOR
"""
function prod_LR(L::AbstractArray{<:Number, 3}, R::AbstractArray{<:Number, 4})

    r1, r2, r3, r4 = size(R)
    l3 = size(L, 3)
    reshape(reshape(L, r1*r2, l3)'*reshape(R, r1*r2, r3*r4), l3, r3, r4)

end

function prod_LR!(A0, L::AbstractArray{<:Number, 3}, R::AbstractArray{<:Number, 4})

    r1, r2, r3, r4 = size(R)
    l3 = size(L, 3)
    A0 = reshape(A0, l3, r3*r4)
    Ac_mul_B!(A0, reshape(L, r1*r2, l3), reshape(R, r1*r2, r3*r4))
    #A0 = reshape(A0, l3, r3, r4)

end
"""
product of right environment and AOL from prod_AOL
"""
function prod_LR(L::AbstractArray{<:Number, 4}, R::AbstractArray{<:Number, 3})

    l1, l2, l3, l4 = size(L)
    r3 = size(R, 3)
    reshape(reshape(L, l1*l2, l3*l4)*reshape(R, l3*l4, r3), l1, l2, r3)

end

function prod_LR!(A0, L::AbstractArray{<:Number, 4}, R::AbstractArray{<:Number, 3})

    l1, l2, l3, l4 = size(L)
    r3 = size(R, 3)
    A0 = reshape(A0, l1*l2, r3)
    A_mul_B!(A0, reshape(L, l1*l2, l3*l4), reshape(R, l3*l4, r3))
    #A0 = reshape(A0, l1, l2, r3)

end
"""
product of left environment and larger right environment or site mps
"""
function prod_LR(L::AbstractArray{<:Number, 2}, R::AbstractArray{<:Number, 3})

    r1, r2, r3 = size(R)
    reshape(L*reshape(R, r1, r2*r3), size(L, 1), r2, r3)

end

function prod_LR_co(L::AbstractArray{<:Number, 2}, R::AbstractArray{<:Number, 3})

    r1, r2, r3 = size(R)
    reshape(L.'*reshape(R, r1, r2*r3), size(L, 2), r2, r3)

end

function prod_LR!(A0, L::AbstractArray{<:Number, 2}, R::AbstractArray{<:Number, 3})

    r1, r2, r3 = size(R)
    l1 = size(L, 1)
    A0 = reshape(A0, l1, r2*r3)
    A_mul_B!(A0, L, reshape(R, r1, r2*r3))
    #A0 = reshape(A0, l1, r2, r3)

end

function prod_LR_co!(A0, L::AbstractArray{<:Number, 2}, R::AbstractArray{<:Number, 3})

    r1, r2, r3 = size(R)
    l2 = size(L, 2)
    A0 = reshape(A0, l2, r2*r3)
    At_mul_B!(A0, L, reshape(R, r1, r2*r3))
    #A0 = reshape(A0, l2, r2, r3)

end
"""
product of right environment and larger left environment or site mps
"""
function prod_LR(L::AbstractArray{<:Number, 3}, R::AbstractArray{<:Number, 2})

    l1, l2, l3 = size(L)
    reshape(reshape(L, l1*l2, l3)*R, l1, l2, size(R, 2))

end

function prod_LR!(A0, L::AbstractArray{<:Number, 3}, R::AbstractArray{<:Number, 2})

    l1, l2, l3 = size(L)
    r2 = size(R, 2)
    A0 = reshape(A0, l1*l2, r2)
    A_mul_B!(A0, reshape(L, l1*l2, l3), R)
    #A0 = reshape(A0, l1, l2, r2)

end

"""
update right environment from AOR and site tensor A0
"""
function update_re(A0::AbstractArray{<:Number, 3},
    AOR::AbstractArray{<:Number, 4})

    r1, r2, r3, r4 = size(AOR)
    a1 = size(A0, 1)
    reshape(reshape(AOR, r1*r2, r3*r4)*reshape(A0, a1, r3*r4)', r1, r2, a1)

end

function update_re!(env, A0::AbstractArray{<:Number, 3},
    AOR::AbstractArray{<:Number, 4})

    r1, r2, r3, r4 = size(AOR)
    a1 = size(A0, 1)
    env = reshape(env, r1*r2, a1)
    A_mul_Bc!(env, reshape(AOR, r1*r2, r3*r4), reshape(A0, a1, r3*r4))
    #env = reshape(env, r1, r2, a1)

end
"""
update right environment from A, O and A0
"""
function update_re(A0::AbstractArray{<:Number, 3},
    O::AbstractArray{<:Number, 4}, A::AbstractArray{<:Number, 3},
    FR::AbstractArray{<:Number, 3})

    AOR = prod_AOR(A, O, FR)
    update_re(A0, AOR)

end

function update_re!(env, A0::AbstractArray{<:Number, 3},
    O::AbstractArray{<:Number, 4}, A::AbstractArray{<:Number, 3},
    FR::AbstractArray{<:Number, 3})

    AOR = prod_AOR(A, O, FR)
    update_re!(env, A0, AOR)

end

"""
update left environment from AR and site tensor A0
"""
function update_renv(A0::AbstractArray{<:Number, 3},
    AR::AbstractArray{<:Number, 3})

    a1, d, a2 = size(A0)
    reshape(AR, size(AR, 1), d*a2)*reshape(A0, a1, d*a2)'

end

function update_renv!(env, A0::AbstractArray{<:Number, 3},
    AR::AbstractArray{<:Number, 3})

    a1, d, a2 = size(A0)
    A_mul_Bc!(env, reshape(AR, size(AR, 1), d*a2), reshape(A0, a1, d*a2))

end

"""
update right environment from A and A0
"""
function update_renv(A0::AbstractArray{<:Number, 3},
    A::AbstractArray{<:Number, 3}, FR::AbstractArray{<:Number, 2})

    AR = prod_LR(A, FR)
    update_renv(A0, AR)

end

function update_renv!(env, A0::AbstractArray{<:Number, 3},
    A::AbstractArray{<:Number, 3}, FR::AbstractArray{<:Number, 2})

    AR = prod_LR(A, FR)
    update_renv!(env, A0, AR)

end
"""
update left environment from AOL and site tensor A0, in the case the left
environment is conjugated
"""
function update_le(A0::AbstractArray{<:Number, 3},
    AOL::AbstractArray{<:Number, 4})

    l1, l2, l3, l4 = size(AOL)
    a3 = size(A0, 3)
    reshape(reshape(AOL, l1*l2, l3*l4)'*reshape(A0, l1*l2, a3), l3, l4, a3)

end

function update_le!(env, A0::AbstractArray{<:Number, 3},
    AOL::AbstractArray{<:Number, 4})

    l1, l2, l3, l4 = size(AOL)
    a3 = size(A0, 3)
    env = reshape(env, l3*l4, a3)
    Ac_mul_B!(env, reshape(AOL, l1*l2, l3*l4), reshape(A0, l1*l2, a3))
    #env = reshape(env, l3, l4, a3)

end

"""
update left environment from A, O and A0
"""
function update_le(A0::AbstractArray{<:Number, 3},
    O::AbstractArray{<:Number, 4}, A::AbstractArray{<:Number, 3},
    FL::AbstractArray{<:Number, 3})

    AOL = prod_AOL(A, O, FL)
    update_le(A0, AOL)

end

function update_le!(env, A0::AbstractArray{<:Number, 3},
    O::AbstractArray{<:Number, 4}, A::AbstractArray{<:Number, 3},
    FL::AbstractArray{<:Number, 3})

    AOL = prod_AOL(A, O, FL)
    update_le!(env, A0, AOL)

end

"""
update left environment from AL and site tensor A0
"""
function update_lenv(A0::AbstractArray{<:Number, 3},
    AL::AbstractArray{<:Number, 3})

    a1, d, a3 = size(A0)
    reshape(A0, a1*d, a3)'*reshape(AL, a1*d, size(AL, 3))

end

function update_lenv!(env, A0::AbstractArray{<:Number, 3},
    AL::AbstractArray{<:Number, 3})

    a1, d, a3 = size(A0)
    Ac_mul_B!(env, reshape(A0, a1*d, a3), reshape(AL, a1*d, size(AL, 3)))

end

function update_lenv_co!(env, A0::AbstractArray{<:Number, 3},
    AL::AbstractArray{<:Number, 3})

    a1, d, a3 = size(A0)
    At_mul_B!(env, reshape(AL, a1*d, size(AL, 3)), reshape(conj(A0), a1*d, a3))

end
"""
update left environment from A, and A0
"""
function update_lenv(A0::AbstractArray{<:Number, 3},
    A::AbstractArray{<:Number, 3}, L::AbstractArray{<:Number, 2})

    AL = prod_LR(L, A)
    update_lenv(A0, AL)

end

function update_lenv!(env, A0::AbstractArray{<:Number, 3},
    A::AbstractArray{<:Number, 3}, L::AbstractArray{<:Number, 2})

    AL = prod_LR(L, A)
    update_lenv!(env, A0, AL)

end

function update_lenv_co!(env, A0::AbstractArray{<:Number, 3},
    A::AbstractArray{<:Number, 3}, L::AbstractArray{<:Number, 2})

    AL = prod_LR_co(L, A)
    update_lenv_co!(env, A0, AL)

end
"""
update left environment from AL without conjugation of site tensor A0
"""
function update_lenv_nc(A0::AbstractArray{<:Number, 3},
    AL::AbstractArray{<:Number, 3})

    a1, d, a3 = size(A0)
    reshape(A0, a1*d, a3).'*reshape(AL, a1*d, size(AL, 3))

end

"""
update left environment from A without conjugation of A0
"""
function update_lenv_nc(A0::AbstractArray{<:Number, 3},
    A::AbstractArray{<:Number, 3}, L::AbstractArray{<:Number, 2})

    AL = prod_LR(L, A)
    update_lenv_nc(A0, AL)

end

"""
local operator expectation value assuming mixed conjugation
"""
function local_exp(A0, A, O)

    a1, d, a3 = size(A)
    Ap = reshape(permutedims(A, [2, 1, 3]), d, a1*a3)
    A0p = reshape(permutedims(A0, [2, 1, 3]), d, a1*a3)
    Ap = Ap*A0p'
    trace(O*Ap)

end