text:
module V3

using DataFrames
using DocStringExtensions
using Printf
using UUIDs
using ...STMMRunner
export run_mstm

abstract type MSTMOutput end

struct NearField
    spheres::DataFrame
    field::DataFrame
end

struct MSTMFixedOutput <: MSTMOutput
    q_ext_unpolarized::Float64
    q_abs_unpolarized::Float64
    q_scat_unpolarized::Float64
    asymmetric_parameter_unpolarized::Float64
    q_ext_parallel::Float64
    q_abs_parallel::Float64
    q_scat_parallel::Float64
    q_ext_perpendicular::Float64
    q_abs_perpendicular::Float64
    q_scat_perpendicular::Float64
    scattering_matrix::DataFrame
    near_field::Union{NearField,Nothing}
end

struct MSTMRandomOutput <: MSTMOutput
    q_ext::Float64
    q_abs::Float64
    q_scat::Float64
    asymmetric_parameter::Float64
    scattering_matrix::DataFrame
    scattering_matrix_expansion_coefficients::DataFrame
    t_matrix::Union{DataFrame,Nothing}
end

"""
Use the given configuration to run MSTM v3.

$(SIGNATURES)

- If `keep = true`, the working directory will not be removed after the run.
- `mstm_exe_name` specifies the name or path of your compiled MSTM v3 executable.
"""
function run_mstm(cfg::STMMConfig; keep::Bool = false, mstm_exe_name::String = "mstm3")
    current_dir = pwd()

    id = string(uuid1())
    if cfg.working_directory == ""
        dir = joinpath(pwd(), "mstm_tmp_$id")
    else
        dir = cfg.working_directory
    end
    results = nothing
    proc = nothing

    try
        if !ispath(dir)
            @debug "[Run MSTM] Creating workding directory"
            mkdir(dir)
        end

        @debug "[Run MSTM] Entering working directory"
        cd(dir)

        @debug "[Run MSTM] Writing MSTM input file"
        write_input(cfg)

        @debug "[Run MSTM] Running MSTM..."
        proc = open(
            `mpiexec -n $(cfg.number_processors) $(mstm_exe_name)`,
            stdout;
            write = true,
        )
        wait(proc)

        @debug "[Run MSTM] Collecting MSTM output"
        results = collect_output(cfg)
    catch e
        @error "[Run MSTM] $e"
        @error "[Run MSTM] Runtime error occurred"
    finally
        @debug "[Run MSTM] Leaving working directory"
        cd(current_dir)

        if !isnothing(proc)
            @debug "[Run MSTM] Killing MSTM process"
            kill(proc)
            sleep(0.2)
            kill(proc) # We need two SIGTERMs to actually kill the process.
            sleep(0.2)
            if process_running(proc)
                @warn "Failed to kill, please retry manually."
            end
        end

        if !keep && occursin(String(id), dir)
            @assert ispath(dir)
            @debug "[Run MSTM] Removing temporary directory"
            rm(dir; force = true, recursive = true)
        end

        return results
    end
end

function format_sphere(sphere::Sphere)::String
    radius_and_center = @sprintf "%.4e %.4e %.4e %.4e " sphere.r sphere.x sphere.y sphere.z

    if !isempty(sphere.t_matrix)
        return radius_and_center * sphere.t_matrix
    else
        return radius_and_center *
               (@sprintf "%.4e %.4e %.4e %.4e" sphere.m.re sphere.m.im sphere.β.re sphere.β.im)
    end
end

# We set `write_sphere_data` to `0` since we do not need them.
# We set `near_field_output_data` to `2` (the full set).
# The manual uses `near_field_translation_distance`, which should be `near_field_distance`.
function write_input(cfg::STMMConfig)
    inp = """
number_spheres
$(length(cfg.spheres))
output_file
$(cfg.output_file)
run_print_file
$(cfg.run_print_file)
write_sphere_data
0
medium_real_ref_index
$(cfg.layers[1].m.re)
medium_imag_ref_index
$(cfg.layers[1].m.im)
medium_real_chiral_factor
$(cfg.layers[1].β.re)
medium_imag_chiral_factor
$(cfg.layers[1].β.im)
target_euler_angles_deg
$(join(cfg.target_euler_angles_deg, ","))
mie_epsilon
$(cfg.mie_epsilon)
translation_epsilon
$(cfg.translation_epsilon)
solution_epsilon
$(cfg.solution_epsilon)
max_number_iterations
$(cfg.max_iterations)
store_translation_matrix
$(Int(cfg.store_translation_matrix))
sm_number_processors
$(cfg.sm_number_processors)
near_field_distance
$(cfg.near_field_translation_distance)
iterations_per_correction
$(cfg.iterations_per_correction)
min_scattering_angle_deg
$(cfg.θ_min)
max_scattering_angle_deg
$(cfg.θ_max)
min_scattering_plane_angle_deg
$(cfg.ϕ_min)
max_scattering_plane_angle_deg
$(cfg.ϕ_max)
delta_scattering_angle_deg
$(cfg.Δθ)
incident_or_target_frame
$(Int(cfg.frame))
normalize_scattering_matrix
$(Int(cfg.normalize_scattering_matrix))
azimuth_average_scattering_matrix
$(Int(cfg.azimuthal_average))
gaussian_beam_constant
$(cfg.gaussian_beam_constant)
gaussian_beam_focal_point
$(join(cfg.gaussian_beam_focal_point, ","))
fixed_or_random_orientation
$(Int(cfg.orientation))
incident_azimuth_angle_deg
$(cfg.α)
incident_polar_angle_deg
$(cfg.β)
calculate_scattering_coefficients
$(Int(cfg.calculate_scattering_coefficients))
scattering_coefficient_file
$(cfg.scattering_coefficient_file)
track_iterations
$(Int(cfg.track_iterations))
calculate_near_field
$(Int(cfg.calculate_near_field))
near_field_plane_coord
$(cfg.near_field_plane_coord)
near_field_plane_position
$(cfg.near_field_plane_position)
near_field_plane_vertices
$(join(cfg.near_field_plane_vertices, ","))
spacial_step_size
$(cfg.near_field_step_size)
polarization_angle_deg
$(cfg.polarization_angle_deg)
near_field_output_file
$(cfg.near_field_output_file)
near_field_output_data
2
plane_wave_epsilon
$(cfg.plane_wave_epsilon)
calculate_t_matrix
$(Int(cfg.calculate_t_matrix))
t_matrix_file
$(cfg.t_matrix_file)
t_matrix_convergence_epsilon
$(cfg.t_matrix_epsilon)
sphere_sizes_and_positions
$(join(map(format_sphere, cfg.spheres),"\n"))
end_of_options
"""

    open("mstm.inp", "w") do io
        write(io, inp)
    end
end

function read_ints(s::AbstractString)::Vector{Int}
    return map(x -> parse(Int, x), filter(!isempty, split(s, r"\s+")))
end

function read_floats(s::AbstractString)::Vector{Float64}
    return map(x -> parse(Float64, x), filter(!isempty, split(s, r"\s+")))
end

function collect_output(cfg::STMMConfig)::MSTMOutput
    @assert isfile(cfg.output_file)

    out = split(read(open(cfg.output_file), String), "\n")
    i = 1

    if cfg.orientation == FixedOrientation
        q = zeros(10)
    else
        q = zeros(4)
    end

    θ = Float64[]
    ϕ = Float64[]
    s11 = Float64[]
    s12 = Float64[]
    s13 = Float64[]
    s14 = Float64[]
    s21 = Float64[]
    s22 = Float64[]
    s23 = Float64[]
    s24 = Float64[]
    s31 = Float64[]
    s32 = Float64[]
    s33 = Float64[]
    s34 = Float64[]
    s41 = Float64[]
    s42 = Float64[]
    s43 = Float64[]
    s44 = Float64[]

    while !occursin("calculation results for run", out[i])
        i += 1
    end

    # Skip header
    i += 2

    if cfg.orientation == FixedOrientation
        # Fixed orientation

        # Read unpolarized efficiencies
        q[1:4] = read_floats(out[i+1])
        i += 2

        # Read parallel efficiencies
        q[5:7] = read_floats(out[i+1])
        i += 2

        # Read perpendicular efficiencies
        q[8:10] = read_floats(out[i+1])
        i += 4

        # Read scattering matrix
        has_phi = !cfg.azimuthal_average
        while !isempty(out[i])
            v = read_floats(out[i])
            i += 1

            push!(θ, v[1])
            push!(ϕ, has_phi ? v[2] : 0.0)

            push!(s44, v[end])
            push!(s43, v[end-1])
            push!(s42, v[end-2])
            push!(s41, v[end-3])
            push!(s34, v[end-4])
            push!(s33, v[end-5])
            push!(s32, v[end-6])
            push!(s31, v[end-7])
            push!(s24, v[end-8])
            push!(s23, v[end-9])
            push!(s22, v[end-10])
            push!(s21, v[end-11])
            push!(s14, v[end-12])
            push!(s13, v[end-13])
            push!(s12, v[end-14])
            push!(s11, v[end-15])
        end
        scattering_matrix = DataFrame(;
            θ,
            ϕ,
            s11,
            s12,
            s13,
            s14,
            s21,
            s22,
            s23,
            s24,
            s31,
            s32,
            s33,
            s34,
            s41,
            s42,
            s43,
            s44,
        )

        if cfg.calculate_near_field && isfile(cfg.near_field_output_file)
            nf = split(read(open(cfg.near_field_output_file), String), "\n")

            # Read near field grid size
            Nx, Ny = read_ints(nf[1])
            Ns = read_ints(nf[2])[1]

            sx = Float64[]
            sy = Float64[]
            sz = fill(cfg.near_field_plane_position, Ns)
            r = Float64[]
            for j = 1:Ns
                v = read_floats(nf[j+2])
                push!(sx, v[1])
                push!(sy, v[2])
                push!(r, v[3])
            end

            if cfg.near_field_plane_coord == 1
                spheres = DataFrame(; x = sz, y = sx, z = sy, r)
            elseif cfg.near_field_plane_coord == 2
                spheres = DataFrame(; x = sy, y = sz, z = sx, r)
            else
                @assert cfg.near_field_plane_coord == 3
                spheres = DataFrame(; x = sx, y = sy, z = sz, r)
            end

            # Read near field data
            E² = Float64[]
            Ex = ComplexF64[]
            Ey = ComplexF64[]
            Ez = ComplexF64[]
            Hx = ComplexF64[]
            Hy = ComplexF64[]
            Hz = ComplexF64[]
            X = Float64[]
            Y = Float64[]
            Z = fill(cfg.near_field_plane_position, Nx * Ny)
            for j = 1:Nx*Ny
                v = read_floats(nf[j+2+Ns])

                push!(X, v[1])
                push!(Y, v[2])

                push!(Ex, complex(v[3], v[4]))
                push!(Ey, complex(v[5], v[6]))
                push!(Ez, complex(v[7], v[8]))

                push!(Hx, complex(v[9], v[10]))
                push!(Hy, complex(v[11], v[12]))
                push!(Hz, complex(v[13], v[14]))
            end

            if cfg.near_field_plane_coord == 1
                x, y, z = Z, X, Y
            elseif cfg.near_field_plane_coord == 2
                x, y, z = Y, Z, X
            else
                @assert cfg.near_field_plane_coord == 3
                x, y, z = X, Y, Z
            end

            field = DataFrame(; x, y, z, Ex, Ey, Ez, Hx, Hy, Hz)
            near_field = NearField(spheres, field)
        else
            near_field = nothing
        end

        return MSTMFixedOutput(q..., scattering_matrix, near_field)
    else
        # Random orientation

        # Read efficiencies
        q[1:4] = read_floats(out[i+1])
        i += 3

        # Read scattering matrix
        Nθ = Int(round((cfg.θ_max - cfg.θ_min) / cfg.Δθ)) + 1
        for j = 1:Nθ
            v = read_floats(out[i+j])

            push!(θ, v[1])
            push!(ϕ, 0.0)

            # Refer to Mishchenko and Yurkin (2017) for the reciprocity relation
            push!(s11, v[2])
            push!(s12, v[3])
            push!(s13, v[4])
            push!(s14, v[5])
            push!(s21, v[3])
            push!(s22, v[6])
            push!(s23, v[7])
            push!(s24, v[8])
            push!(s31, -v[4])
            push!(s32, -v[7])
            push!(s33, v[9])
            push!(s34, v[10])
            push!(s41, v[5])
            push!(s42, v[8])
            push!(s43, -v[10])
            push!(s44, v[11])
        end
        scattering_matrix = DataFrame(;
            θ,
            ϕ,
            s11,
            s12,
            s13,
            s14,
            s21,
            s22,
            s23,
            s24,
            s31,
            s32,
            s33,
            s34,
            s41,
            s42,
            s43,
            s44,
        )
        i += Nθ + 3

        # Read scattering matrix expansion coefficients
        w = Int[]
        a11 = Float64[]
        a12 = Float64[]
        a13 = Float64[]
        a14 = Float64[]
        a22 = Float64[]
        a23 = Float64[]
        a24 = Float64[]
        a32 = Float64[]
        a33 = Float64[]
        a34 = Float64[]
        a44 = Float64[]

        while !isempty(out[i])
            v = read_floats(out[i])
            i += 1

            push!(w, Int(v[1]))
            push!(a11, v[2])
            push!(a22, v[3])
            push!(a33, v[4])
            push!(a23, v[5])
            push!(a32, v[6])
            push!(a44, v[7])
            push!(a12, v[8])
            push!(a34, v[9])
            push!(a13, v[10])
            push!(a24, v[11])
            push!(a14, v[12])
        end
        scattering_matrix_expansion_coefficients =
            DataFrame(; w, a11, a12, a13, a14, a22, a23, a24, a32, a33, a34, a44)

        return MSTMRandomOutput(
            q...,
            scattering_matrix,
            scattering_matrix_expansion_coefficients,
            nothing,
        )
    end
end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "CyrusSASL"
version = v"2.1.27"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/cyrusimap/cyrus-sasl/releases/download/cyrus-sasl-2.1.27/cyrus-sasl-2.1.27.tar.gz", "26866b1549b00ffd020f188a43c258017fa1c382b3ddadd8201536f72efb05d5"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/cyrus-sasl-*/
if [[ "${target}" == *-mingw* ]]; then
    # Patches from
    # https://github.com/msys2/MINGW-packages/tree/72dac6ecd4dbbe2cbb29ebbb355b6742a9103150/mingw-w64-cyrus-sasl
    atomic_patch -p1 ../patches/02-exeext.patch
    atomic_patch -p1 ../patches/03-fix-plugins.patch
    atomic_patch -p1 ../patches/04-manpage-paths.patch
    atomic_patch -p1 ../patches/14-MinGW-w64-add-LIBSASL_API-to-function-definitions.patch
    atomic_patch -p1 ../patches/15-MinGW-w64-define-LIBSASL_EXPORTS_eq_1-for-sasldb.patch
    atomic_patch -p1 ../patches/16-MinGW-w64-define-WIN32_LEAN_AND_MEAN-avoiding-handle_t-redef.patch
    atomic_patch -p1 ../patches/17-MinGW-w64-define-S_IRUSR-and-S_IWUSR.patch
    atomic_patch -p1 ../patches/19-paths-relocation.patch
    atomic_patch -p1 ../patches/20-mingw-tchar.patch
    atomic_patch -p1 ../patches/21-fix-getopt-guard.patch

    # Remove incompatible typedef
    atomic_patch -p1 ../patches/30-remove-extra-incompatible-typedef.patch
fi
atomic_patch -p1 ../patches/macos-shared-lib-extension.patch
autoreconf -vi
if [[ "${target}" == *-mingw* ]]; then
    # Copy the right header file for Windows into `include/`...
    cp win32/include/md5global.h include/md5global.h

    # ...and don't regenerate a wrong one with `make.  This patch needs to be
    # applied _after_ autoreconf, which seems to somehow revert the changes :-(
    atomic_patch -p1 ../patches/31-do-not-make-mdf5global_h.patch
fi
./configure --prefix=${prefix} \
    --build=${MACHTYPE} --host=${target} \
    --with-openssl=${prefix} \
    --with-sqlite3=${prefix} \
    --oldincludedir=${includedir} \
    --enable-ntlm \
    --disable-gssapi \
    --with-dblib=gdbm \
    --disable-static \
    --disable-ldapdb \
    --without-saslauthd \
    --without-pwcheck \
    --without-des \
    --without-authdaemond \
    --disable-sample \
    --with-plugindir=${prefix}/lib/sasl2 \
    --with-configdir=${prefix}/etc/sasl2:${prefix}/etc/sasl:${prefix}/lib/sasl2
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libanonymous", :libanonymous, "lib/sasl2"),
    LibraryProduct("libplain", :libplan, "lib/sasl2"),
    LibraryProduct("libscram", :libscram, "lib/sasl2"),
    LibraryProduct("libotp", :libotp, "lib/sasl2"),
    LibraryProduct("libdigestmd5", :libdigestmd5, "lib/sasl2"),
    LibraryProduct("libcrammd5", :libcrammd5, "lib/sasl2"),
    LibraryProduct("libntlm", :libntlm, "lib/sasl2"),
    LibraryProduct("libsasl2", :libsasl2)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="OpenSSL_jll", uuid="458c3c95-2e84-50aa-8efc-19380b2a3a95"))
    Dependency(PackageSpec(name="SQLite_jll", uuid="76ed43ae-9a5d-5a62-8c75-30186b810ce8"))
    Dependency(PackageSpec(name="Gdbm_jll", uuid="54ca2031-c8dd-5cab-9ed4-295edde1660f"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

ext: jl

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function compute_xi_t!(gauge_t::Gauge, bulkconstrain::BulkConstrained,
                       bulkevol::BulkEvolved, deriv::BulkDeriv, gauge::Gauge,
                       cache::HorizonCache, sys::System{Outer}, ::ConstantGauge)
    xi_t = getxi(gauge_t)
    fill!(xi_t, 0)
    nothing
end

function compute_xi_t!(gauge_t::Gauge, bulkconstrain::BulkConstrained,
                       bulkevol::BulkEvolved, deriv::BulkDeriv, gauge::Gauge,
                       cache::HorizonCache, sys::System{Outer}, gaugecondition::Advect_xi)
    _, Nx, Ny = size(sys)

    Dx  = sys.Dx
    Dy  = sys.Dy

    xiGF = getxi(gauge)
    xi_t = getxi(gauge_t)

    vx   = gaugecondition.xi_vx
    vy   = gaugecondition.xi_vy

    @fastmath @inbounds Threads.@threads for j in 1:Ny
        @inbounds for i in 1:Nx
            xi_x  = Dx(xiGF, 1,i,j)
            xi_y  = Dy(xiGF, 1,i,j)

            xi_t[1,i,j] = -vx * xi_x - vy * xi_y
        end
    end

    nothing
end

#= Solving for xi_t

this is a 2D PDE of the type

  (axx Dxx + ayy Dyy + Axy Dxy + bx Dx + by Dy + c Id) xi_t = -S

we first build each operator (Dxx, Dyy, etc) through a Kronecker product and
then overwrite each of them with the corresponding coefficient. we then sum them
all up and solve the linear system.

we use SparseMatrices since these are finite differencing operators. also, there
is no reason to use the same operators as the ones we use elsewhere. here, as
we're solving a very large linear system, it may be better to use second order
accurate operators since the resulting matrix is smaller and therefore much
faster to invert.

since this is a gauge condition, i think it shouldn't affect the overall
convergence order.

=#
function compute_xi_t!(gauge_t::Gauge, bulkconstrain::BulkConstrained,
                       bulkevol::BulkEvolved, deriv::BulkDeriv, gauge::Gauge,
                       cache::HorizonCache, sys::System{Outer},
                       gaugecondition::ConstantAH)
    _, Nx, Ny = size(sys)
    bulk = Bulk(bulkevol, bulkconstrain)

    Du  = sys.Du
    Dx  = sys.Dx
    Dxx = sys.Dxx
    Dy  = sys.Dy
    Dyy = sys.Dyy

    interp = sys.uinterp

    B1_uAH      = cache.bulkhorizon.B1_uAH
    B2_uAH      = cache.bulkhorizon.B2_uAH
    G_uAH       = cache.bulkhorizon.G_uAH
    phi_uAH     = cache.bulkhorizon.phi_uAH
    S_uAH       = cache.bulkhorizon.S_uAH
    Fx_uAH      = cache.bulkhorizon.Fx_uAH
    Fy_uAH      = cache.bulkhorizon.Fy_uAH
    Sd_uAH      = cache.bulkhorizon.Sd_uAH
    B1d_uAH     = cache.bulkhorizon.B1d_uAH
    B2d_uAH     = cache.bulkhorizon.B2d_uAH
    Gd_uAH      = cache.bulkhorizon.Gd_uAH
    phid_uAH    = cache.bulkhorizon.phid_uAH
    A_uAH       = cache.bulkhorizon.A_uAH
    Du_B1_uAH   = cache.bulkhorizon.Du_B1_uAH
    Du_B2_uAH   = cache.bulkhorizon.Du_B2_uAH
    Du_G_uAH    = cache.bulkhorizon.Du_G_uAH
    Du_phi_uAH  = cache.bulkhorizon.Du_phi_uAH
    Du_S_uAH    = cache.bulkhorizon.Du_S_uAH
    Du_Fx_uAH   = cache.bulkhorizon.Du_Fx_uAH
    Du_Fy_uAH   = cache.bulkhorizon.Du_Fy_uAH
    Du_Sd_uAH   = cache.bulkhorizon.Du_Sd_uAH
    Du_B1d_uAH  = cache.bulkhorizon.Du_B1d_uAH
    Du_B2d_uAH  = cache.bulkhorizon.Du_B2d_uAH
    Du_Gd_uAH   = cache.bulkhorizon.Du_Gd_uAH
    Du_A_uAH    = cache.bulkhorizon.Du_A_uAH
    Duu_B1_uAH  = cache.bulkhorizon.Duu_B1_uAH
    Duu_B2_uAH  = cache.bulkhorizon.Duu_B2_uAH
    Duu_G_uAH   = cache.bulkhorizon.Duu_G_uAH
    Duu_S_uAH   = cache.bulkhorizon.Duu_S_uAH
    Duu_Fx_uAH  = cache.bulkhorizon.Duu_Fx_uAH
    Duu_Fy_uAH  = cache.bulkhorizon.Duu_Fy_uAH
    Duu_A_uAH   = cache.bulkhorizon.Duu_A_uAH

    axx         = cache.axx
    ayy         = cache.ayy
    axy         = cache.axy
    bx          = cache.bx
    by          = cache.by
    cc          = cache.cc
    b_vec       = cache.b_vec

    Dx_2D       = cache.Dx_2D
    Dxx_2D      = cache.Dxx_2D
    Dy_2D       = cache.Dy_2D
    Dyy_2D      = cache.Dyy_2D
    Dxy_2D      = cache.Dxy_2D
    _Dx_2D      = cache._Dx_2D
    _Dxx_2D     = cache._Dxx_2D
    _Dy_2D      = cache._Dy_2D
    _Dyy_2D     = cache._Dyy_2D
    _Dxy_2D     = cache._Dxy_2D

    xi_t = getxi(gauge_t)

    uAH   = gaugecondition.u_AH
    kappa = gaugecondition.kappa

    u2 = uAH * uAH
    u3 = uAH * uAH * uAH
    u4 = uAH * uAH * uAH * uAH

    # take u-derivatives of Sd, B1d, B2d and Gd. these are not needed in the
    # nested system, so they haven't been computed before. since we need them
    # here, compute them now
    @sync begin
        @spawn mul!(deriv.Du_Sd,  Du,  bulkconstrain.Sd)
        @spawn mul!(deriv.Du_B1d, Du,  bulkconstrain.B1d)
        @spawn mul!(deriv.Du_B2d, Du,  bulkconstrain.B2d)
        @spawn mul!(deriv.Du_Gd,  Du,  bulkconstrain.Gd)
    end

    # interpolate bulk functions (and u-derivatives) to the u=uAH surface
    @inbounds Threads.@threads for j in 1:Ny
        @inbounds for i in 1:Nx
            B1_uAH[1,i,j]       = interp(view(bulk.B1,  :,i,j))(uAH)
            B2_uAH[1,i,j]       = interp(view(bulk.B2,  :,i,j))(uAH)
            G_uAH[1,i,j]        = interp(view(bulk.G,   :,i,j))(uAH)
            phi_uAH[1,i,j]      = interp(view(bulk.phi, :,i,j))(uAH)
            S_uAH[1,i,j]        = interp(view(bulk.S,   :,i,j))(uAH)
            Fx_uAH[1,i,j]       = interp(view(bulk.Fx,  :,i,j))(uAH)
            Fy_uAH[1,i,j]       = interp(view(bulk.Fy,  :,i,j))(uAH)
            Sd_uAH[1,i,j]       = interp(view(bulk.Sd,  :,i,j))(uAH)
            B1d_uAH[1,i,j]      = interp(view(bulk.B1d, :,i,j))(uAH)
            B2d_uAH[1,i,j]      = interp(view(bulk.B2d, :,i,j))(uAH)
            Gd_uAH[1,i,j]       = interp(view(bulk.Gd,  :,i,j))(uAH)
            phid_uAH[1,i,j]     = interp(view(bulk.phid,:,i,j))(uAH)
            A_uAH[1,i,j]        = interp(view(bulk.A,   :,i,j))(uAH)

            Du_B1_uAH[1,i,j]    = interp(view(deriv.Du_B1,  :,i,j))(uAH)
            Du_B2_uAH[1,i,j]    = interp(view(deriv.Du_B2,  :,i,j))(uAH)
            Du_G_uAH[1,i,j]     = interp(view(deriv.Du_G,   :,i,j))(uAH)
            Du_phi_uAH[1,i,j]   = interp(view(deriv.Du_phi, :,i,j))(uAH)
            Du_S_uAH[1,i,j]     = interp(view(deriv.Du_S,   :,i,j))(uAH)
            Du_Fx_uAH[1,i,j]    = interp(view(deriv.Du_Fx,  :,i,j))(uAH)
            Du_Fy_uAH[1,i,j]    = interp(view(deriv.Du_Fy,  :,i,j))(uAH)
            Du_Sd_uAH[1,i,j]    = interp(view(deriv.Du_Sd,  :,i,j))(uAH)
            Du_B1d_uAH[1,i,j]   = interp(view(deriv.Du_B1d, :,i,j))(uAH)
            Du_B2d_uAH[1,i,j]   = interp(view(deriv.Du_B2d, :,i,j))(uAH)
            Du_Gd_uAH[1,i,j]    = interp(view(deriv.Du_Gd,  :,i,j))(uAH)
            Du_A_uAH[1,i,j]     = interp(view(deriv.Du_A,   :,i,j))(uAH)

            Duu_B1_uAH[1,i,j]   = interp(view(deriv.Duu_B1,  :,i,j))(uAH)
            Duu_B2_uAH[1,i,j]   = interp(view(deriv.Duu_B2,  :,i,j))(uAH)
            Duu_G_uAH[1,i,j]    = interp(view(deriv.Duu_G,   :,i,j))(uAH)
            Duu_S_uAH[1,i,j]    = interp(view(deriv.Duu_S,   :,i,j))(uAH)
            Duu_Fx_uAH[1,i,j]   = interp(view(deriv.Duu_Fx,  :,i,j))(uAH)
            Duu_Fy_uAH[1,i,j]   = interp(view(deriv.Duu_Fy,  :,i,j))(uAH)
            Duu_A_uAH[1,i,j]    = interp(view(deriv.Duu_A,   :,i,j))(uAH)
        end
    end

    ind2D  = LinearIndices(B1_uAH[1,:,:])

    # coefficients of the derivative operators
    @fastmath @inbounds Threads.@threads for j in 1:Ny
        @inbounds for i in 1:Nx
            idx   = ind2D[i,j]

            xi    = gauge.xi[1,i,j]
            xi_x  = Dx(gauge.xi, 1,i,j)
            xi_y  = Dy(gauge.xi, 1,i,j)
            xi_xx = Dxx(gauge.xi, 1,i,j)
            xi_yy = Dyy(gauge.xi, 1,i,j)
            xi_xy = Dx(Dy, gauge.xi, 1,i,j)

            B1    = B1_uAH[1,i,j]
            B2    = B2_uAH[1,i,j]
            G     = G_uAH[1,i,j]
            phi   = phi_uAH[1,i,j]
            S     = S_uAH[1,i,j]
            Fx    = Fx_uAH[1,i,j]
            Fy    = Fy_uAH[1,i,j]
            Sd    = Sd_uAH[1,i,j]
            B1d   = B1d_uAH[1,i,j]
            B2d   = B2d_uAH[1,i,j]
            Gd    = Gd_uAH[1,i,j]
            phid  = phid_uAH[1,i,j]
            A     = A_uAH[1,i,j]

            # r derivatives

            B1p        = -u2 * Du_B1_uAH[1,i,j]
            B2p        = -u2 * Du_B2_uAH[1,i,j]
            Gp         = -u2 * Du_G_uAH[1,i,j]
            phip       = -u2 * Du_phi_uAH[1,i,j]
            Sp         = -u2 * Du_S_uAH[1,i,j]
            Fxp        = -u2 * Du_Fx_uAH[1,i,j]
            Fyp        = -u2 * Du_Fy_uAH[1,i,j]
            Sdp        = -u2 * Du_Sd_uAH[1,i,j]
            B1dp       = -u2 * Du_B1d_uAH[1,i,j]
            B2dp       = -u2 * Du_B2d_uAH[1,i,j]
            Gdp        = -u2 * Du_Gd_uAH[1,i,j]
            Ap         = -u2 * Du_A_uAH[1,i,j]

            B1pp       = 2*u3 * Du_B1_uAH[1,i,j]  + u4 * Duu_B1_uAH[1,i,j]
            B2pp       = 2*u3 * Du_B2_uAH[1,i,j]  + u4 * Duu_B2_uAH[1,i,j]
            Gpp        = 2*u3 * Du_G_uAH[1,i,j]   + u4 * Duu_G_uAH[1,i,j]
            Spp        = 2*u3 * Du_S_uAH[1,i,j]   + u4 * Duu_S_uAH[1,i,j]
            Fxpp       = 2*u3 * Du_Fx_uAH[1,i,j]  + u4 * Duu_Fx_uAH[1,i,j]
            Fypp       = 2*u3 * Du_Fy_uAH[1,i,j]  + u4 * Duu_Fy_uAH[1,i,j]
            App        = 2*u3 * Du_A_uAH[1,i,j]   + u4 * Duu_A_uAH[1,i,j]

            # x and y derivatives

            B1_x    = Dx(B1_uAH,   1,i,j)
            B2_x    = Dx(B2_uAH,   1,i,j)
            G_x     = Dx(G_uAH,    1,i,j)
            phi_x   = Dx(phi_uAH,  1,i,j)
            S_x     = Dx(S_uAH,    1,i,j)
            Fx_x    = Dx(Fx_uAH,   1,i,j)
            Fy_x    = Dx(Fy_uAH,   1,i,j)
            Sd_x    = Dx(Sd_uAH,   1,i,j)
            B1d_x   = Dx(B1d_uAH,  1,i,j)
            B2d_x   = Dx(B2d_uAH,  1,i,j)
            Gd_x    = Dx(Gd_uAH,   1,i,j)
            A_x     = Dx(A_uAH,    1,i,j)

            B1_y    = Dy(B1_uAH,   1,i,j)
            B2_y    = Dy(B2_uAH,   1,i,j)
            G_y     = Dy(G_uAH,    1,i,j)
            phi_y   = Dy(phi_uAH,  1,i,j)
            S_y     = Dy(S_uAH,    1,i,j)
            Fx_y    = Dy(Fx_uAH,   1,i,j)
            Fy_y    = Dy(Fy_uAH,   1,i,j)
            Sd_y    = Dy(Sd_uAH,   1,i,j)
            B1d_y   = Dy(B1d_uAH,  1,i,j)
            B2d_y   = Dy(B2d_uAH,  1,i,j)
            Gd_y    = Dy(Gd_uAH,   1,i,j)
            A_y     = Dy(A_uAH,    1,i,j)

            B1p_x   = -u2 * Dx(Du_B1_uAH, 1,i,j)
            B2p_x   = -u2 * Dx(Du_B2_uAH, 1,i,j)
            Gp_x    = -u2 * Dx(Du_G_uAH,  1,i,j)
            Sp_x    = -u2 * Dx(Du_S_uAH,  1,i,j)
            Fxp_x   = -u2 * Dx(Du_Fx_uAH, 1,i,j)
            Fyp_x   = -u2 * Dx(Du_Fy_uAH, 1,i,j)
            Ap_x    = -u2 * Dx(Du_A_uAH,  1,i,j)

            B1p_y   = -u2 * Dy(Du_B1_uAH, 1,i,j)
            B2p_y   = -u2 * Dy(Du_B2_uAH, 1,i,j)
            Gp_y    = -u2 * Dy(Du_G_uAH,  1,i,j)
            Sp_y    = -u2 * Dy(Du_S_uAH,  1,i,j)
            Fxp_y   = -u2 * Dy(Du_Fx_uAH, 1,i,j)
            Fyp_y   = -u2 * Dy(Du_Fy_uAH, 1,i,j)
            Ap_y    = -u2 * Dy(Du_A_uAH,  1,i,j)

            Fy_xx   = Dxx(Fy_uAH,  1,i,j)
            A_xx    = Dxx(A_uAH,   1,i,j)

            Fx_yy   = Dyy(Fx_uAH,  1,i,j)
            A_yy    = Dyy(A_uAH,   1,i,j)

            Fx_xy   = Dx(Dy, Fx_uAH, 1,i,j)
            Fy_xy   = Dx(Dy, Fy_uAH, 1,i,j)
            A_xy    = Dx(Dy, A_uAH,  1,i,j)

            vars =  (
                kappa, xi, xi_x, xi_y, xi_xx, xi_yy, xi_xy,
                B1   , B2   , G   , phi  , S    , Fx    , Fy    , Sd ,  B1d  , B2d  , Gd,  phid, A   ,
                B1p  , B2p  , Gp  , phip , Sp   , Fxp   , Fyp   , Sdp,  B1dp , B2dp , Gdp,       Ap  ,
                B1pp , B2pp , Gpp ,        Spp  , Fxpp  , Fypp  ,                                App ,
                B1_x , B2_x , G_x , phi_x, S_x  , Fx_x  , Fy_x  , Sd_x, B1d_x, B2d_x, Gd_x,      A_x ,
	        B1_y , B2_y , G_y , phi_y, S_y  , Fx_y  , Fy_y  , Sd_y, B1d_y, B2d_y, Gd_y,      A_y ,
                B1p_x, B2p_x, Gp_x,        Sp_x , Fxp_x , Fyp_x ,                                Ap_x,
                B1p_y, B2p_y, Gp_y,        Sp_y , Fxp_y , Fyp_y ,                                Ap_y,
                                                          Fy_xx ,                                A_xx,
                                                  Fx_yy ,                                        A_yy,
                                                  Fx_xy , Fy_xy ,                                A_xy
            )

            a11, a22, a12, b1, b2, c, SS = xi_t_eq_coeff(vars, sys.gridtype)

            axx[idx]   = a11
            ayy[idx]   = a22
            axy[idx]   = a12
            bx[idx]    = b1
            by[idx]    = b2
            cc[idx]    = c

            b_vec[idx] = -SS
        end
    end

    axxId = Diagonal(axx)
    ayyId = Diagonal(ayy)
    axyId = Diagonal(axy)
    bxId  = Diagonal(bx)
    byId  = Diagonal(by)
    ccId  = Diagonal(cc)

    # build the differential operators by multiplying with the coefficients
    # computed in the loop above (note that Dxx_2D, Dyy_2D, etc, are never
    # overwritten)
    mul!(_Dxx_2D, axxId, Dxx_2D)
    mul!(_Dyy_2D, ayyId, Dyy_2D)
    mul!(_Dxy_2D, axyId, Dxy_2D)
    mul!(_Dx_2D,  bxId,  Dx_2D)
    mul!(_Dy_2D,  byId,  Dy_2D)

    # build actual operator to be inverted
    A_mat = _Dxx_2D + _Dyy_2D + _Dxy_2D + _Dx_2D + _Dy_2D + ccId

    # solve system
    A_fact = lu(A_mat)
    ldiv!(A_fact, b_vec)

    copyto!(xi_t, b_vec)

    nothing
end


text:
module ServoWithRampAndPrismatic

using ModiaLang
using Unitful
import Modia3D

# ModiaLang models
include("$(ModiaLang.path)/models/Blocks.jl")
include("$(ModiaLang.path)/models/Electric.jl")
include("$(ModiaLang.path)/models/Translational.jl")

import Modia3D
using  Modia3D.ModiaInterface

ControllerPrism = Model(
    # Interface
    refLoadForce     = input,
    actualMotorPos   = input,
    actualMotorSpeed = input,
    refForce         = output,
    # Components
    gainOuter     = Gain,
    posFeedback   = Feedback,
    gainInner     = Gain,
    speedFeedback = Feedback,
    PI = PI,

    connect = :[
        (refLoadForce, gainOuter.u)
        (gainOuter.y,  posFeedback.u1)
        (actualMotorPos, posFeedback.u2)
        (posFeedback.y,  gainInner.u)
        (gainInner.y,    speedFeedback.u1)
        (actualMotorSpeed, speedFeedback.u2)
        (speedFeedback.y,  PI.u)
        (PI.y, refForce)
    ]
)

DrivePrism = Model(
    # Interface
    refForce    = input,
    actualSpeed = output,
    actualPos   = output,
    flange      = TranslationalFlange,
    # Components
    force        = UnitlessForce,
    speedSensor  = UnitlessVelocitySensor,
    posSensor    = UnitlessPositionSensor,
#    motorInertia = Inertia,
#    idealGear    = IdealGear,

    connect = :[
        (refForce, force.f)
        (force.flange, speedSensor.flange, posSensor.flange, flange)
        #(motorInertia.flange_b, idealGear.flange_a)
        #(idealGear.flange_b, flange)
        (speedSensor.v, actualSpeed)
        (posSensor.s, actualPos)
    ]
)

ServoPrism = Model(
    refLoadForce = input,
    flange       = TranslationalFlange,
    ppi          = ControllerPrism,
    drive        = DrivePrism,

    connect = :[
        # inputs of ppi
        (refLoadForce, ppi.refLoadForce)
        (ppi.actualMotorSpeed, drive.actualSpeed)
        (ppi.actualMotorPos,   drive.actualPos)
        # output of ppi --> input of motor
        (ppi.refForce, drive.refForce)
        # Flange of drive
        (drive.flange, flange)
    ]
)

arm_joint_1_obj = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/arm_joint_1.obj")

m = 2
axis = 3
vmat1 = VisualMaterial(color="LightBlue", transparency=0.5)
m1=1.390
translation1 = [0.033, 0, 0]
rotation1 = Modia3D.rot1(180u"°")

k1 = 50.0
k2 = 0.1
T2 = 1.0
gearRatio = 156.0
motorInertia = (0.0000135 + 0.000000409)*u"kg*m^2"
J1 = 0.32

servoParameters = Map(
                      ppi = Map(
                          gainOuter = Map(k=gearRatio),
                          gainInner = Map(k=k1),
                          PI        = Map(k=k2, T=T2)
                      ),
                      drive = Map(
                          #motorInertia = Map(J=motorInertia),
                          #idealGear = Map(ratio=gearRatio)
                      )
                  )

TestServo = Model(
    world = Object3D(feature=Scene()),
    body  = Object3D(feature=Solid(shape=FileMesh(filename=arm_joint_1_obj),
                                   massProperties=MassPropertiesFromShapeAndMass(mass=m1),
                                   visualMaterial=vmat1)),
    obj2  = Object3D(parent=:body,
                     translation=translation1,
                     rotation=rotation1),

    rev = PrismaticWithFlange(obj1=:world, obj2=:obj2, axis=axis, s=Var(init=0.0), v=Var(init=0.0)),

    ramp  = Ramp | Map(duration=1.18u"s", height=2.95),
    servo = ServoPrism | servoParameters,
    connect = :[
        (ramp.y, servo.refLoadForce)
        (servo.flange, rev.flange)
    ]
)

servo = @instantiateModel(buildModia3D(TestServo), unitless=true, logCode=true, log=true)

stopTime = 4.0
tolerance = 1e-6
requiredFinalStates = [-269.9911219480055, 1594.1933857866297, 25436.042732813803]
simulate!(servo, stopTime=stopTime, tolerance=tolerance, log=true, logStates=true, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plotVariables = [("ramp.y", "rev.s"); "rev.v"; "servo.ppi.PI.x"; "servo.ppi.refForce"]
plot(servo, plotVariables, figure=1)

end


text:
using IntervalArithmetic, IntervalOptimisation
using Test
using IntervalOptimisation: numeric_type

@testset "IntervalOptimisation tests" begin
    @testset "numeric_type" begin
        x = -10..10
        big_x = big(x)
        @test numeric_type(x) == Float64
        @test numeric_type(big_x) == BigFloat
        @test numeric_type(IntervalBox(x, x)) == Float64
        @test numeric_type(IntervalBox(big_x, big_x)) == BigFloat
    end

    @testset "Minimise in 1D using default data structure i.e HeapedVector" begin
        global_min, minimisers = minimise(x->x, -10..10)
        @test global_min ⊆ -10 .. -9.999
        @test length(minimisers) == 1
        @test minimisers[1] ⊆ -10 .. -9.999

        # same but with maximise
        global_max, maximisers = maximise(x->x, -10..10)
        @test global_max ⊆ 9.999 .. 10
        @test length(maximisers) == 1
        @test maximisers[1] ⊆ 9.999 .. 10

        # same but with BigFloats
        global_min, minimisers = minimise(x->x, -big(10.0)..big(10.0))
        @test global_min ⊆ -10 .. -9.999
        @test length(minimisers) == 1
        @test minimisers[1] ⊆ -10 .. -9.999

        global_min, minimisers = minimise(x->x^2, -10..11, tol = 1e-10)
        @test global_min ⊆ 0..1e-20
        @test length(minimisers) == 1
        @test minimisers[1] ⊆ -0.1..0.1

        global_min, minimisers = minimise(x->(x^2-2)^2, -10..11)
        @test global_min ⊆ 0..1e-7
        @test length(minimisers) == 2
        @test sqrt(2) ∈ minimisers[2]
    end

    for Structure in (SortedVector, HeapedVector)

        @testset "Minimise in 1D using SortedVector" begin
            global_min, minimisers = minimise(x->x, -10..10, structure = Structure)
            @test global_min ⊆ -10 .. -9.999
            @test length(minimisers) == 1
            @test minimisers[1] ⊆ -10 .. -9.999

            # same but with maximise
            global_max, maximisers = maximise(x->x, -10..10, structure = Structure)
            @test global_max ⊆ 9.999 .. 10
            @test length(maximisers) == 1
            @test maximisers[1] ⊆ 9.999 .. 10

            # same but with BigFloats
            global_min, minimisers = minimise(x->x, -big(10.0)..big(10.0), structure = Structure)
            @test global_min ⊆ -10 .. -9.999
            @test length(minimisers) == 1
            @test minimisers[1] ⊆ -10 .. -9.999

            global_min, minimisers = minimise(x->x^2, -10..11, tol=1e-10, structure = Structure)
            @test global_min ⊆ 0..1e-20
            @test length(minimisers) == 1
            @test minimisers[1] ⊆ -0.1..0.1

            global_min, minimisers = minimise(x->(x^2-2)^2, -10..11, structure = Structure)
            @test global_min ⊆ 0..1e-7
            @test length(minimisers) == 2
            @test sqrt(2) ∈ max(minimisers[1], minimisers[2])
        end

        @testset "Discontinuous function in 1D" begin

            H(x) = (sign(x) + 1) / 2   # Heaviside function except at 0, where H(0) = 0.5
            global_min, minimisers = minimise(x -> abs(x) + H(x) - 1, -10..11, tol=1e-5, structure = Structure)
            @test global_min ⊆ -1 .. -0.9999
            @test length(minimisers) == 1
            @test 0 ∈ minimisers[1]
            @test diam(minimisers[1]) <= 1e-5
        end

        @testset "Smooth function in 2D" begin
            global_min, minimisers = minimise( X -> ( (x,y) = X; x^2 + y^2 ), (-10..10) × (-10..10), structure = Structure )
            @test global_min ⊆ 0..1e-7
            @test all(X ⊆ (-1e-3..1e3) × (-1e-3..1e-3) for X in minimisers)

            # same but with maximise
            global_max, maximisers = maximise( X -> ( (x,y) = X; x^2 + y^2 ), (-10..10) × (-10..10), structure = Structure )
            @test global_max ⊆ 199.9..200
            m = (9.99..10)
            @test all(X ⊆ m × m || X ⊆ -m × m || X ⊆ m × -m || X ⊆ -m × -m for X in maximisers)

            # same but with BigFloats
            global_min, minimisers = minimise( X -> ( (x,y) = X; x^2 + y^2 ), (-big(10.0)..big(10.0)) × (-big(10.0)..big(10.0)), structure = Structure )
            @test global_min ⊆ 0..1e-7
            @test all(X ⊆ big(-1e-3..1e3) × big(-1e-3..1e-3) for X in minimisers)

        end

        @testset "Non smooth function in 2D" begin
            global_min, minimisers = minimise( X -> ( (x,y) = X; abs(x-y) + abs(x) ), (-1..1) × (-1..1), structure = Structure )
            @test global_min == 0..0 # it must return the exact minimum (if mid is used to find the first global_min)
            @test all(X ⊆ (-2e-3.. 2e-3) × (-2e-3..2e-3) for X in minimisers)
        end

    end

end


text:
module TestExamplesTreeMesh3DPart1

using Test
using Trixi

include("test_trixi.jl")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "TreeMesh3D Part 1" begin

# Run basic tests
@testset "Examples 3D" begin
  # Compressible Euler
  include("test_tree_3d_euler.jl")
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end # TreeMesh3D Part 1

end #module


text:
#function logsumexp(x)
#    max_x = maximum(x)
#    exp_x = exp.(x .- max_x)
#    log(sum(exp_x)) .+ max_x
#end

function mbar_equation!(F, f_k, K, N_max, N_k, u_kln, idx, log_wi_jn, log_term, FlogN)
    for k = 1:K
        FlogN[k] = log(Float64(N_k[k])) + f_k[k]
    end
    for i = 1:K
        for j = 1:K
            for n = 1:N_k[j]
                ui_xjn = u_kln[j, i, n];
                for k = 1:K
                    uk_xjn = u_kln[j, k, n];
                    log_term[k] = FlogN[k] - (uk_xjn - ui_xjn);
                end
                log_wi_jn[j, n] = -logsumexp(log_term);
            end
        end
        F[i] = f_k[i] + logsumexp(log_wi_jn[idx])
    end
    F[:] = F[:] .- F[1]
end

"""
    mbar(u_kl; ftol=1e-8, iterations=10^2) -> F

Estimates the free energy differences of umbrella-windowed systems by using the Multistate Bennet Acceptance Ratio Method (MBAR).
Let K be # of umbrellas or different ensembles. `u_kl` is a K x K Array whose elements are reduced bias-factor or
potential energy of umbrella (or ensemble) simulation data k evaluated at umbrella (or ensemble) l.

Returns (dimensionless) free energies of umbrella-windows or ensembles `f_k`.

# References
```
M. R. Shirts and J. D. Chodera, J. Chem. Phys. 129, 124105 (2008).
```
"""
function mbar(u_kl; ftol=1e-8, iterations=10^2)
    # K: number of umbrella windows
    K, L = size(u_kl)

    # N_k: number of data in k-th umbrella window
    N_k = zeros(Int64, K)
    for k = 1:K
        N_k[k] = length(u_kl[k, 1])
    end
    N_max = maximum(N_k)

    # conversion from array of array (u_kl) to array (u_kln)
    u_kln = zeros(Float64, K, K, N_max)
    for k = 1:K
        for l = 1:K
            u_kln[k, l, 1:N_k[k]] = u_kl[k, l]
        end
    end

    dummy_kn = zeros(Int64, K, N_max)
    for k = 1:K
      dummy_kn[k, 1:N_k[k]] .= 1
    end
    idx = dummy_kn .> 0;

    # solve the MBAR equation
    log_wi_jn = zeros(Float64, K, N_max)
    log_term = zeros(Float64, K)
    FlogN = zeros(Float64, K)
    # f!(F, x) = mbar_equation!(F, x, K, N_max, N_k, u_kln, idx, log_wi_jn, log_term, FlogN)
    # x_init = zeros(Float64, K)
    # sol = nlsolve(f!, x_init, autodiff=:forward, ftol=ftol, iterations=iterations)
    # #sol = nlsolve(f!, x_init, method = :anderson, ftol=ftol, iterations=iterations)
    # f_k = sol.zero
    # f_k .- f_k[1]

    F = zeros(Float64, K)
    x_init = zeros(Float64, K)
    for n = 1:1000
        mbar_equation!(F, x_init, K, N_max, N_k, u_kln, idx, log_wi_jn, log_term, FlogN)
        x_init = x_init - F
        x_init = x_init .- x_init[1]
    end
    f_k = x_init
end


text:
using SpecialFunctions
"""
    gammaf()
Variables: 0.5 <= x <= 5;
Objective: Minimize gamma(x)
"""

function gammaf()
    best_known_objective = gamma(1.4616321)
    cplex_optimizer =
        JuMP.optimizer_with_attributes(CPLEX.Optimizer, "CPX_PARAM_SCRIND" => 0)
    error_tolerances = [NaN64, 1e-3]
    # base partition holds the inflection points and the end-points
    milp_statistics = zeros(length(error_tolerances), 5)
    lp_statistics = zeros(length(error_tolerances), 5)

    base_partition = [0.5, 5]
    # functions
    f = x -> gamma(x)

        # solve MILP for various error tolerances
        for i = 1:length(error_tolerances)
            error_tolerance = error_tolerances[i]
            # construct MILP relaxations
            milp = Model(cplex_optimizer)
            var, x_index, y_index = add_milp_relaxation(
                    f,
                    base_partition,
                    error_tolerance,
                    milp,
                    name = "1",
                )
            @objective(milp, Min, var[y_index])
            _, solve_time, solve_bytes_alloc, sec_in_gc = @timed optimize!(milp)
            relaxation_objective = objective_value(milp)
            @assert relaxation_objective ≤ best_known_objective
            relative_gap =
                abs(best_known_objective - relaxation_objective) / (abs(relaxation_objective) + 1e-4)
            milp_statistics[i, :] = [
                error_tolerance,
                best_known_objective,
                relaxation_objective,
                relative_gap,
                round(solve_time, digits = 2),
            ]
        end

        for i = 1:length(error_tolerances)
            error_tolerance = error_tolerances[i]
            # construct LP relaxations
            lp = Model(cplex_optimizer)
            var, x_index, y_index = add_lp_relaxation(
                    f,
                    base_partition,
                    error_tolerance,
                    lp,
                    name = "1",
                )
            @objective(lp, Min, var[y_index])
            _, solve_time, solve_bytes_alloc, sec_in_gc = @timed optimize!(lp)
            relaxation_objective = objective_value(lp)
            @assert relaxation_objective ≤ best_known_objective
            relative_gap =
                abs(best_known_objective - relaxation_objective) / abs(relaxation_objective)
            lp_statistics[i, :] = [
                error_tolerance,
                best_known_objective,
                relaxation_objective,
                relative_gap,
                round(solve_time, digits = 2),
            ]
        end

        open("csv/gammaf_milp.csv", "w") do io
            writedlm(io, ["error_tol" "g_opt" "relaxation" "relative_gap" "time"], ',')
            writedlm(io, milp_statistics, ",")
        end

        open("csv/gammaf_lp.csv", "w") do io
            writedlm(io, ["error_tol" "g_opt" "relaxation" "relative_gap" "time"], ',')
            writedlm(io, lp_statistics, ",")
        end

end

max_stars_repo_name: sujeevraja/univariate-experiments

max_issues_repo_name: sujeevraja/univariate-experiments

max_forks_repo_name: sujeevraja/univariate-experiments


text:
using ITensors
using ITensorGaussianMPS
using LinearAlgebra

# Electrons

# Half filling
Nx, Ny = 6, 3
N = Nx * Ny
Nf = N

@show Nx, Ny
@show N, Nf

# Maximum MPS link dimension
_maxlinkdim = 1_000

@show _maxlinkdim

# DMRG cutoff
_cutoff = 1e-5

# Hopping
t = 1.0

# Electron-electon on-site interaction
U = 4.0

@show t, U

lattice = square_lattice(Nx, Ny; yperiodic=true)

# Make the free fermion Hamiltonian for the up spins
os_up = OpSum()
for b in lattice
  os_up .+= -t, "Cdagup", b.s1, "Cup", b.s2
  os_up .+= -t, "Cdagup", b.s2, "Cup", b.s1
end

# Make the free fermion Hamiltonian for the down spins
os_dn = OpSum()
for b in lattice
  os_dn .+= -t, "Cdagdn", b.s1, "Cdn", b.s2
  os_dn .+= -t, "Cdagdn", b.s2, "Cdn", b.s1
end

# Hopping Hamiltonian with 2*N spinless fermions,
# alternating up and down spins
h = hopping_hamiltonian(os_up, os_dn)

# Get the Slater determinant
Φ = slater_determinant_matrix(h, Nf)

println()
println("Exact free fermion energy: ", tr(Φ'h * Φ))
println()

# Create an MPS from the slater determinant.
# In this example we are turning of spin conservation.
s = siteinds("Electron", N; conserve_qns=true, conserve_sz=false)
println("Making free fermion starting MPS")
@time ψ0 = slater_determinant_to_mps(
  s, Φ; eigval_cutoff=1e-4, cutoff=_cutoff, maxdim=_maxlinkdim
)
@show maxlinkdim(ψ0)

os = os_up + os_dn
for n in 1:N
  os .+= U, "Nupdn", n
end
H = MPO(os, s)

# Random starting state
ψr = randomMPS(s, n -> n ≤ Nf ? (isodd(n) ? "↑" : "↓") : "0")

println("\nRandom starting state energy")
@show flux(ψr)
@show inner(ψr, H, ψr)

println("\nFree fermion MPS starting state energy")
@show flux(ψ0)
@show inner(ψ0, H, ψ0)

println("\nStart from random product state")
sweeps = Sweeps(10)
setmaxdim!(sweeps, 10, 20, 100, 200, _maxlinkdim)
setcutoff!(sweeps, _cutoff)
setnoise!(sweeps, 1e-7, 1e-8, 1e-10, 0.0)
@time dmrg(H, ψr, sweeps)

println("\nStart from free fermion state")
sweeps = Sweeps(10)
setmaxdim!(sweeps, _maxlinkdim)
setcutoff!(sweeps, _cutoff)
@time dmrg(H, ψ0, sweeps)

nothing


text:
using OpenSCADPolyhedron
using OpenSCADIcosahedron

xscale = 1.15470053838
hscale = 1.51152262815
l1 = .3 # length toward center of triangle in plane from edge
l2 = .1 # overlap length of plexiglass over outer support (must be less than l1)
h = .3  # total height of inner rib
d = .15  # thickness of outer support (so depth of plexiglass traingles behind outer faces)
alpha = 1.2059325 # = half the dihedral angle of an icosahedron ~ 69 degrees

# FIGURE 1:
# Cross section parameters and points layout
#
#              l1
#    _____________________
# d |                    / |
#   |________           /  |
#      l2    |         /   | h
#            |________/    |

#=
icosahedronPoints = [
  Point3{Float64}([-1., -0.57735, 1.51152]),
  Point3{Float64}([ 0.,  1.1547,  1.51152]),
  Point3{Float64}([ 1., -0.57735, 1.51152]),
]
=#

# generates the points that define the triangle crosssection in the xz plane
# points are defined counter clockwise starting with top right point in Figure 1
xpoints = [
  Point3{Float64}([xscale,              0, hscale]),
  Point3{Float64}([xscale - 2l1,        0, hscale]),
  Point3{Float64}([xscale - 2l1,        0, hscale - d]),
  Point3{Float64}([xscale - 2(l1 - l2), 0, hscale - d]),
  Point3{Float64}([xscale - 2(l1 - l2), 0, hscale - h]),
  Point3{Float64}([xscale - 2*h/tan(alpha), 0, hscale - h]),
]
temppoly1 = Polyhedron{Float64}(xpoints)

strutTransforms = Function[
  p::Point3{Float64} -> p,
  p::Point3{Float64} -> Point3{Float64}([ -0.5 -0.866025403784 0.; 0.866025403784 -0.5 0.; 0. 0. 1.] * p.x),
]

triangleTransforms = Function[
  p::Point3{Float64} -> Point3{Float64}([ 0.866025403784  0.5  0.; -0.5  0.866025403784  0.; 0. 0. 1.] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ 0.             -1.   0.;  1.   0.              0.; 0. 0. 1.] * p.x),
  p::Point3{Float64} -> Point3{Float64}([-0.866025403784  0.5  0.; -0.5 -0.866025403784  0.; 0. 0. 1.] * p.x),
]

temppoly2 = generatePolyhedron(temppoly1, strutTransforms, Dict{(Int, Int), (Int, Int)}())
append!(temppoly2, Int[1, 2, 8, 7])
append!(temppoly2, Int[2, 3, 9, 8])
append!(temppoly2, Int[3, 4, 10, 9])
append!(temppoly2, Int[4, 5, 11, 10])
append!(temppoly2, Int[5, 6, 12, 11])
#append!(temppoly2, Int[])
temppoly2symmetries = [
  (2,  1) => (1,  7),
  (2,  2) => (1,  8),
  (2,  3) => (1,  9),
  (2,  4) => (1, 10),
  (2,  5) => (1, 11),
  (2,  6) => (1, 12),
  (3,  7) => (1,  1),
  (3,  8) => (1,  2),
  (3,  9) => (1,  3),
  (3, 10) => (1,  4),
  (3, 11) => (1,  5),
  (3, 12) => (1,  6),
  (3,  1) => (2,  7),
  (3,  2) => (2,  8),
  (3,  3) => (2,  9),
  (3,  4) => (2, 10),
  (3,  5) => (2, 11),
  (3,  6) => (2, 12),
]

temppoly3 = generatePolyhedron(temppoly2, triangleTransforms, temppoly2symmetries)
#println(stringify(temppoly3))

layers = Set{Array{Int,1}}()
push!(layers, Int[13, 7, 1])
push!(layers, Int[18, 12, 6])

symmetries = expandSymmetries(icosahedronSymmetries, icosahedronFace, layers)
lampahedron = generatePolyhedron(temppoly3, icosahedronTransforms, symmetries)

println(stringify(lampahedron))

max_stars_repo_name: erikstrand/OpenSCADPolyhedron

max_issues_repo_name: erikstrand/OpenSCADPolyhedron

max_forks_repo_name: erikstrand/OpenSCADPolyhedron


text:
using Test
using Animations
using Colors

@testset "keyframes" begin

    kf1 = Keyframe(1, 5.0)
    kf2 = Keyframe(3, 10.0)
    kf3 = Keyframe(6, 0.0)

    animation = Animation(
        [kf1, kf2, kf3],
        [sineio(), linear()]
    )

    @test timestamps(animation) == [1, 3, 6]
    @test keyvalues(animation) == [5.0, 10.0, 0.0]
end

@testset "animation creation" begin
    anim = Animation(
        1, 5,
        2, 10,
        saccadic(2),
        3, 20,
        sineio(),
        4, 100
    )

    @test keyvalues(anim) == [5, 10, 20, 100]
    @test timestamps(anim) == [1, 2, 3, 4]
    @test easings(anim) == [linear(), saccadic(2), sineio()]
end

@testset "animation creation rel values" begin
    anim = Animation(
        0, 3,
        rel(2), 4,
        rel(1), 3
    )
    @test timestamps(anim) == [0, 2, 3]
end

@testset "animation creation defaulteasing" begin
    e = polyin(3; n=3, yoyo=true)
    anim = Animation(
        1, 5,
        2, 10,
        saccadic(2),
        3, 20,
        sineio(),
        4, 100;
        defaulteasing=e
    )

    @test keyvalues(anim) == [5, 10, 20, 100]
    @test timestamps(anim) == [1, 2, 3, 4]
    @test easings(anim) == [e, saccadic(2), sineio()]
end

@testset "temporally changed animations" begin
    anim = Animation(
        1, 10,
        2, 20,
        3, 30,
    )

    anim_later = anim + 1
    @test timestamps(anim_later) == [2, 3, 4]
    anim_earlier = anim - 1
    @test timestamps(anim_earlier) == [0, 1, 2]
    anim_longer = anim * 2
    @test timestamps(anim_longer) == [2, 4, 6]
    anim_shorter = anim / 2
    @test timestamps(anim_shorter) == [0.5, 1, 1.5]
end

@testset "vector interpolate" begin

    kf1 = Keyframe(0, [0.0, 0.0, 0.0])
    kf2 = Keyframe(1, [1.0, 2.0, 3.0])

    animation = Animation(
        [kf1, kf2],
        sineio()
    )

    @test animation(0) == [0.0, 0.0, 0.0]
    @test animation(0.5) == [0.5, 1.0, 1.5]
    @test animation(1) == [1.0, 2.0, 3.0]
end

@testset "color interpolate" begin

    c1 = RGB(0.0, 0, 0)
    c2 = RGB(1, 0.5, 0.3)
    kf1 = Keyframe(0, c1)
    kf2 = Keyframe(1, c2)

    animation = Animation(
        [kf1, kf2],
        sineio()
    )

    @test animation(0) == c1
    @test animation(0.5) == Animations.Colors.weighted_color_mean(0.5, c1, c2)
    @test animation(1) == c2
end

@testset "string steps" begin

    kf1 = Keyframe(0, "first")
    kf2 = Keyframe(1, "second")
    kf3 = Keyframe(2, "third")

    animation = Animation(
        [kf1, kf2, kf3],
        noease()
    )

    @test animation(0) == "first"
    @test animation(0.5) == "first"
    @test animation(1) == "second"
    @test animation(1.5) == "second"
    @test animation(2) == "third"
end

@testset "yoyo repeat" begin

    animation = Animation(
        [0, 3],
        [0.0, 1.0],
        linear(n=3, yoyo=true)
    )

    results = [0.0, 0.25, 0.5, 0.75, 1.0, 0.75, 0.5, 0.25, 0.0, 0.25, 0.5, 0.75, 1.0]
    @test animation.(collect(0:0.25:3)) == results

end

@testset "even yoyo" begin
    animation = Animation(
        [1, 3],
        [0.0, 1.0],
        linear(n=2, yoyo=true)
    )

    @test animation.(0:4) == [0, 0, 1, 0, 0]
end

@testset "mixed easing" begin
    animation = Animation(
        [0, 1],
        [0, 1],
        mixed(sineio(), polyin(2))
    )
    @test animation(0) == 0
    @test animation(1) == 1
end

@testset "multiplied easing" begin
    animation = Animation(
        [0, 1],
        [0, 1],
        multiplied(sineio(), polyin(2))
    )
    @test animation(0) == 0
    @test animation(1) == 1
end

@testset "eased easing" begin
    animation = Animation(
        [0, 1],
        [0, 1],
        eased(sineio(), polyin(2), saccadic(2))
    )
    @test animation(0) == 0
    @test animation(1) == 1
end

@testset "loop" begin
    animation = Animation(
        0, 0.0,
        sineio(n=2, yoyo=true),
        2, 1.0
    )
    loop = Loop(animation, 0, 1, 2)
    @test loop(0) == 0
    @test loop(1) == 1
    @test loop(2) == 0
    @test loop(2.5) == 0
    @test loop(3) == 0
    @test loop(4) == 1
    @test loop(5) == 0
    @test loop(7) == 0 # two loops done
end

@testset "sequence" begin
    anim1 = Animation(
        0, 0.0,
        1, 1.0,
    )
    anim2 = Animation(
        0, 1.0,
        1, 0.0
    )
    sequence = Sequence([anim1, anim2], 1, 1)
    @test sequence(0) == 0.0
    @test sequence(1) == 0.0
    @test sequence(1.5) == 0.5
    @test sequence(2) == 1.0
    @test sequence(2.5) == 1.0
    @test sequence(3) == 1.0
    @test sequence(3.5) == 0.5
    @test sequence(4) == 0.0
    @test sequence(5) == 0.0
end

@testset "sequence with loop" begin
    anim1 = Animation(
        0, 0.0,
        1, 1.0,
    )
    anim2 = Animation(
        0, 1.0,
        1, 0.0
    )
    loop = Loop(anim2, 0, 1, 2)
    sequence = Sequence([anim1, loop], 1, 1)
    @test sequence(0) == 0.0
    @test sequence(1) == 0.0
    @test sequence(1.5) == 0.5
    @test sequence(2) == 1.0
    @test sequence(2.5) == 1.0
    # start of loop
    @test sequence(3) == 1.0
    @test sequence(3.5) == 0.5
    @test sequence(4) == 0.0
    @test sequence(4.5) == 0.0
    @test sequence(5) == 1.0
    @test sequence(5.5) == 0.5
    @test sequence(6) == 0.0
    @test sequence(7) == 0.0
end

@testset "animate_async" begin

    anim = Animation([0, 2], [5.0, 10.0])

    loopanim = Animation(
        0, 4.0,
        sineio(n=2, yoyo=true),
        2, 9.0
    )
    loop = Loop(loopanim, 0, 1, 2)

    anim1 = Animation(
        0, 3.0,
        1, 1.0,
    )
    anim2 = Animation(
        0, 1.0,
        1, 8.0
    )
    sequence = Sequence([anim1, anim2], 1, 1)

    vals = [0.0, 0.0, 0.0]
    last_t = Ref(0.0)
    animate_async(anim, loop, sequence; duration = 0.1, fps = 30) do t, an, lo, se
        last_t[] = t
        vals[1] = an
        vals[2] = lo
        vals[3] = se
    end |> wait

    t = last_t[]
    @test vals == [anim(t), loop(t), sequence(t)]

    starttime = time()
    timeref = Ref(time())
    # start long loop to then stop it prematurely
    animtask = animate_async(Loop(anim, 0, 0, 1000); duration = 5) do t, lo
        timeref[] = time()
    end
    stop(animtask)
    # give the task the opportunity to keep running
    sleep(0.5)

    @test timeref[] - starttime < 0.1 # stupid way to test that the loop didn't run much after stopping
end


text:
##
# Begin general norm routines
##
function normsq_l2(x::Matrix{T}, grid1::Vector, grid2::Vector) where {T<:Real}
    v0 = zero(T)
    sz1, sz2 = size(x)
    for j in 2:sz2
        dx = grid2[j] - grid2[j-1]
        for i in 2:sz1
            v0 += abs2(x[i, j]) * dx * (grid1[i] - grid1[i-1])
        end
    end

    v0
end
function normsq_l2(x::Matrix{T}, dx::Real, dt::Real) where {T<:Real}
    v0 = zero(T)
    sz1, sz2 = size(x)
    for j in 2:sz2
        for i in 2:sz1
            v0 += abs2(x[i, j])
        end
    end

    v0 * dt * dx
end
normsq_l2(x::Matrix, grid1::AbstractRange, grid2::AbstractRange) = normsq_l2(x, step(grid1), step(grid2))
norm_l2(x::Matrix, grid1, grid2) = sqrt(normsq_l2(x, grid1, grid2))

function norm_l1(x::Matrix{T}, grid1::Vector, grid2::Vector) where {T<:Real}
    v0 = zero(T)
    sz1, sz2 = size(x)
    for j in 2:sz2
        dx = grid2[j] - grid2[j-1]
        for i in 2:sz1
            v0 += abs(x[i, j]) * dx * (grid1[i] - grid1[i-1])
        end
    end

    v0
end
function norm_l1(x::Matrix{T}, dx::Real, dt::Real) where {T<:Real}
    v0 = zero(T)
    sz1, sz2 = size(x)
    for j in 2:sz2
        for i in 2:sz1
            v0 += abs(x[i, j])
        end
    end

    v0 * dx * dt
end
norm_l1(x::AbstractMatrix, grid1::AbstractRange, grid2::AbstractRange) = norm_l1(x, step(grid1), step(grid2))
##
# End general norm routines
##

export norm_l2, norm_w21_only, norm_w22_only, norm_w21, norm_w22, norm_null

function normsq_l2(x::Vector{T}, grid::Vector) where {T<:Real}
    gnext = first(grid)
    v0 = zero(T)
    for i in Iterators.drop(eachindex(x), 1)
        gcurr, gnext = gnext, grid[i]
        v0 += abs2(x[i]) * (gnext - gcurr)
    end
    v0
end
function normsq_l2(x::Vector{T}, dx::Real) where {T<:Real}
    sum(abs2, Iterators.drop(x, 1))*dx
end
norm_l2(x::Vector, grid::Vector) = sqrt(normsq_l2(x, grid))
norm_l2(x::Vector, dx::Real) = sqrt(normsq_l2(x, dx))

##
# Begin discrete Sobolev norm routines
##

function normsq_w21_only(x::Vector{T}, grid::Vector) where {T<:Real}
    xnext, gnext, v0 = first(x), first(grid), zero(T)
    for i in Iterators.drop(eachindex(x), 1)
        xcurr, xnext, gcurr, gnext = xnext, x[i], gnext, grid[i]
        v0 += abs2(xnext - xcurr) / (gnext - gcurr)
    end
    v0
end
# Below implementation is much faster for large (i.e. >512 points.)
function normsq_w21_only(x::Vector{T}, dx::Real) where {T<:Real}
    xnext, v0 = first(x), zero(T)
    for i in Iterators.drop(eachindex(x), 1)
        xcurr, xnext = xnext, x[i]
        v0 += abs2(xnext - xcurr)
    end
    v0 / dx
end

norm_w21_only(x::Vector, grid::Vector) = sqrt(normsq_w21_only(x, grid))
norm_w21_only(x::Vector, dx::Real) = sqrt(normsq_w21_only(x, dx))

normsq_w21(x::Vector, grid::Vector) = normsq_l2(x, grid) + normsq_w21_only(x, grid)
normsq_w21(x::Vector, dx::Real) = normsq_l2(x, dx) + normsq_w21_only(x, dx)

norm_w21(x::Vector, grid::Vector) = sqrt(normsq_w21(x, grid))
norm_w21(x::Vector, dx::Real) = sqrt(normsq_w21(x, dx))

function normsq_w22_only(x::Vector{T}, grid::Vector) where {T<:Real}
    #   Note: to handle the general case with full accuracy, a more involved algorithm
    #   is required, leading this to be quite a bit more computationally intensive...
    nx = length(x)
    #   v0 = zero(T)
    if nx < 3
        return NaN
        #     return v0
    end
    xprev, xcurr, xnext, dx1, dx2 = x[1], x[2], x[3], grid[2] - grid[1], grid[3] - grid[2]
    rdx1, rdx2, rdxs = 1 / dx1, 1 / dx2, 1/(dx1 + dx2)
    v0 = 4*abs2(rdx1 * rdxs * xprev - rdx1 * rdx2 * xcurr + rdx2 * rdxs * xnext)
    if nx == 3
        return v0 * dx1
    end
    for i in 3:(nx - 1)
        xprev, xcurr, xnext, dx1, dx2 = xcurr, xnext, x[i + 1], dx2, grid[i + 1] - grid[i]
        rdx1, rdx2, rdxs = rdx2, 1 / dx2, 1 / (dx1 + dx2)
        v0 += 4*abs2(rdx1 * rdxs * xprev - rdx1 * rdx2 * xcurr + rdx2 * rdxs * xnext)
    end
    v0 * dx1
end
function normsq_w22_only(x::Vector, dx::Real)
    nx = length(x)
    if nx < 3
        return NaN
    end
    xprev, xcurr, xnext = x[1], x[2], x[3]
    v0 = abs2(xprev - 2*xcurr + xnext)
    if nx == 3
        return v0 / dx^3
    end
    for i in 3:(nx - 1)
        xprev, xcurr, xnext = xcurr, xnext, x[i + 1]
        v0 += abs2(xprev - 2*xcurr + xnext)
    end
    v0 / dx^3
end

norm_w22_only(x::Vector, grid::Vector) = sqrt(normsq_w22_only(x, grid))
norm_w22_only(x::Vector, dx::Real) = sqrt(normsq_w22_only(x, dx))

normsq_w22(x::Vector, grid::Vector) = normsq_w21(x, grid) + normsq_w22_only(x, grid)
normsq_w22(x::Vector, dx::Real) = normsq_w21(x, dx) + normsq_w22_only(x, dx)

norm_w22(x::Vector, grid::Vector) = sqrt(normsq_w22(x, grid))
norm_w22(x::Vector, dx::Real) = sqrt(normsq_w22(x, dx))

#TODO: Either generalize norms to work against any AbstractArray, or call this
# exportedVectorNorms.
const exportedNorms = [norm_l2,
                       normsq_w21, norm_w21, normsq_w21_only, norm_w21_only,
                       normsq_w22, norm_w22, normsq_w22_only, norm_w22_only
                       ]
const exportedNormSymbols = [:norm_l2,
                             :normsq_w21, :norm_w21, :normsq_w21_only, :norm_w21_only,
                             :normsq_w22, :norm_w22, :normsq_w22_only, :norm_w22_only
                             ]
##
# End discrete Sobolev norm routines
##

## Define functions on AbstractRanges using definitions with constant grid spacing
for fn in exportedNormSymbols
    @eval ($fn)(x::Vector, grid::AbstractRange) = ($fn)(x, step(grid))
    for T in [Float32, Float64, BigFloat]
        @eval precompile($fn, (Vector{$T}, $T))
        @eval precompile($fn, (Vector{$T}, Vector{$T}))
    end
end


text:
nums = parse.(Int, readlines("day10/input.txt"))

push!(nums, 0)
push!(nums, 3+maximum(nums))
sort!(nums)

println(nums)

dp = Dict()
function dp_proxy(i)
    if i in keys(dp)
        return dp[i]
    else
        result = connect_suffix(i)
        dp[i] = result
        return result
    end
end

function connect_suffix(i)
    global nums
    num_ways = 0
    if i == length(nums)
        return 1
    end
    if i + 1 <= length(nums) && (nums[i]+1 <= nums[i + 1] <= nums[i]+3)
        num_ways += dp_proxy(i+1)
    end
    if i + 2 <= length(nums) && (nums[i]+1 <= nums[i + 2] <= nums[i]+3)
        num_ways += dp_proxy(i+2)
    end
    if i + 3 <= length(nums) && (nums[i]+1 <= nums[i + 3] <= nums[i]+3)
        num_ways += dp_proxy(i+3)
    end
    return num_ways
end

println(connect_suffix(1))

# holy hell DP! that was exciting!


text:
module TestUTF8

include("preamble.jl")

using SymSpellChecker: add_edits!

@testset "add_edits! should generate correct delete sequence for ascii" begin
    considered_deletes = Set{String}()
    candidates = Set{String}()
    candidate = "world"
    candidate_len = length(candidate)
    add_edits!(considered_deletes, candidates, candidate, candidate_len)
    @test length(candidates) == 5
    @test candidates == Set(["orld", "wrld", "wold", "word", "worl"])
end

@testset "add_edits! should generate correct delete sequence for utf-8" begin
    considered_deletes = Set{String}()
    candidates = Set{String}()
    candidate = "земля"
    candidate_len = length(candidate)
    add_edits!(considered_deletes, candidates, candidate, candidate_len)
    @test length(candidates) == 5
    @test candidates == Set(["емля", "змля", "зеля", "земя", "земл"])
end

end # module


text:
export humanize!, humanize

using ARFIMA, Statistics
const φ0 = ARFIMA.SVector(-0.5, -1.5)

"""
    humanize!(notes, property, σ, noise = :ARFIMA; kwargs...)
Humanize given `notes` by adding noise of standard deviation `σ` to their `property`,
typically either `:position` or `:velocity`. Research[^Datseris2019] suggests that `σ`
should be around 40 (for `:position`) but that depends on the BPM, and around 10 for `:velocity`.

The `noise` argument decides the type of noise:
* `:ARFIMA` uses ARFIMA.jl and attempts to generate a power-law correlated (pink) noise.
  Keywords `d = 0.25, φ = SVector(-0.5, -1.5)` are propagated to function `arfima`.
* `:white` plain ol' white noise.

Use `humanize` for a non-modifying version.

[^Datseris2019]: Datseris, G., et al. [Microtiming Deviations and Swing Feel in Jazz. Sci Rep 9, 19824 (2019)](https://doi.org/10.1038/s41598-019-55981-3)
"""
function humanize!(n::Notes, property, σ, noise = :ARFIMA; d=0.25, φ=φ0)
    N = length(n)
    if noise == :ARFIMA
        ξ = arfima(N, 1.0, d, φ)
    elseif noise == :white
        ξ = randn(N)
    else
        error("Unrecognized noise type")
    end
    ξ .*= σ/std(ξ)
    for j in 1:N
        setproperty!(n[j], property, getproperty(n[j], property) + round(Int, ξ[j]))
    end
    return n
end

humanize(notes, args...; kwargs...) = humanize!(copy(notes), args...; kwargs...)


text:
using JuMP
using Ipopt
using PowerModels
using StochasticPowerModels

# constants
const PM = PowerModels
const SPM = StochasticPowerModels

# solvers
ipopt_solver = Ipopt.Optimizer

# input
deg  = 2
aux  = true
red  = false
case = "case5_spm.m"

# data
file  = joinpath(BASE_DIR, "test/data/matpower", case)

#-----------------------------------
# run the convenience functions for stochastic OPF for IVR and ACR
result_ivr = run_sopf_iv(file, PM.IVRPowerModel, ipopt_solver, aux=aux, deg=deg, red=red)
result_acr = run_sopf_acr(file, PM.ACRPowerModel, ipopt_solver, aux=aux, deg=deg)

@assert result_ivr["termination_status"] == PM.LOCALLY_SOLVED
@assert result_acr["termination_status"] == PM.LOCALLY_SOLVED
#the optimal objective values (expectation) are
obj_ivr = result_ivr["objective"]
obj_acr = result_acr["objective"]
@assert obj_ivr ≈ obj_acr

# print variables for all polynomial indices k
SPM.print_summary(result_ivr["solution"])

# print variables for a specific index k
k=1
SPM.print_summary(result_ivr["solution"]["nw"]["$k"])

# get polynomial chaos coefficients for specific component
pg_coeff = pce_coeff(result_ivr, "gen", 1, "pg")

# obtain 10 samples of the generator active power output variable
pg_sample = sample(result_ivr, "gen", 1, "pg"; sample_size=10)

# obtain an kernel density estimate of the generator active power output variable
pg_density = density(result_ivr, "gen", 1, "pg"; sample_size=10)

#-----------------------------------
# alternatively, you can first read in PowerModels dict,
# from a file with stochastic data extensions:
data  = PM.parse_file(file)

result_ivr2 = run_sopf_iv(data, PM.IVRPowerModel, ipopt_solver; aux=aux, deg=deg)
@assert result_ivr2["termination_status"] == PM.LOCALLY_SOLVED
obj_ivr2 = result_ivr2["objective"]

#-----------------------------------
# finally, we can also build the multinetwork dictionary here.
# we run the stochastic replicate function, introduces polynomial index k
sdata = SPM.build_stochastic_data(data, deg)

# run the reduced IVR with auxiliary variables
result_ivr3 = PM.run_model(sdata, PM.IVRPowerModel, ipopt_solver, SPM.build_sopf_iv_with_aux; multinetwork=true, solution_processors=[PM.sol_data_model!])
@assert result_ivr3["termination_status"] == PM.LOCALLY_SOLVED
obj_ivr3 = result_ivr3["objective"]

@assert obj_ivr ≈ obj_ivr2 ≈ obj_ivr3


text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

# sensitive dependence
Random.seed!(5)
fig = figure(figsize=(figx/2, figy))
ax = gca()
u0 = [10,10.0,10]
lo = Systems.lorenz([10,10.0,10])

for i in 1:3
    u = u0 .+ i*1e-3
    tr = trajectory(lo, 15, u)
    plot(0:0.01:15, tr[:, 1]; c = COLORS[[1,4,3][i]])
end

xlabel("\$t\$",labelpad = -15)
yticks(-15:10:15)
xticks(0:5:15)
ylabel("\$x\$",labelpad = -20)
tight_layout(pad = 0.25)
wsave(plotsdir("3","sensitive"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
module Buttons

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export btn, btngroup, Btn

register_normal_element("q__btn", context = @__MODULE__)
register_normal_element("q__btn__group", context = @__MODULE__)

function btn( label::String = "",
              args...;
              wrap::Function = StippleUI.DEFAULT_WRAPPER,
              content::Union{String,Vector,Function} = "",
              kwargs...)

  wrap() do
    q__btn([isa(content, Function) ? content() : join(content)],
            args...;
            attributes([:label => label, kwargs...], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
  end
end

function btn( content::Union{Function,Vector},
              label::String = "",
              args...;
              wrap::Function = StippleUI.DEFAULT_WRAPPER,
              kwargs...)
  btn(label, args...; wrap = wrap, content = content, kwargs...)
end

function btn( args...;
              label::String = "",
              content::Union{String,Vector,Function} = "",
              wrap::Function = StippleUI.DEFAULT_WRAPPER,
              kwargs...)
  btn(label, args...; wrap = wrap, content = content, kwargs...)
end

mutable struct Btn
  label
  args
  wrap
  content
  kwargs

  Btn(label::String = "",
      args...;
      wrap::Function = StippleUI.DEFAULT_WRAPPER,
      content::Union{String,Vector,Function} = "",
      kwargs...) = new(fieldname, args, wrap, mask, kwargs)
end

Base.string(bt::Btn) = string(bt.label, bt.args...; bt.wrap, bt.content, bt.kwargs...)

function btngroup(args...;
                    wrap::Function = StippleUI.DEFAULT_WRAPPER,
                    kwargs...)
  wrap() do
    Genie.Renderer.Html.q__btn__group(args...; kwargs...)
  end
end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ a9ed5ba5-9f59-4bbf-b35c-4a5a36d0b6e0
using DelimitedFiles

# ╔═╡ 291ffe60-5c63-415a-a5f7-1e40e59e2cf9
# easier way with IterTools
using IterTools

# ╔═╡ 748af037-7d04-4d73-9aaa-575c48fac5f3
depths = readdlm("inputs/day_01.txt")[:, 1]

# ╔═╡ 844a02ca-0744-4523-99a8-37918bff34b3
md"""
# Part 1

To do this, count the number of times a depth measurement increases from the previous measurement. (There is no measurement before the first measurement.) In the example above, the changes are as follows.
"""

# ╔═╡ 6467bbda-bf07-4ee0-81a3-ae4511c6e8eb
# the ugly way
begin
deeper_counts = 0

for (i, depth) in enumerate(depths[1:end-1])
	if depth < depths[i + 1]
		global deeper_counts += 1
	end
end

deeper_counts
end

# ╔═╡ 3407599f-b7f6-4b02-a36b-5792e4111be6
# the neat way with broadcasting
sum((depths[2:end] .- depths[1:end-1]) .> 0)

# ╔═╡ 6da7f6d6-755b-43cf-aac2-b99299222955
md"""Part 2
Considering every single measurement isn't as useful as you expected: there's just too much noise in the data.

Instead, consider sums of a three-measurement sliding window. Again considering the above example:

```
199  A
200  A B
208  A B C
210    B C D
200  E   C D
207  E F   D
240  E F G
269    F G H
260      G H
263        H
```

Start by comparing the first and second three-measurement windows. The measurements in the first window are marked A (199, 200, 208); their sum is 199 + 200 + 208 = 607. The second window is marked B (200, 208, 210); its sum is 618. The sum of measurements in the second window is larger than the sum of the first, so this first comparison increased.

Your goal now is to count the number of times the sum of measurements in this sliding window increases from the previous sum. So, compare A with B, then compare B with C, then C with D, and so on. Stop when there aren't enough measurements left to create a new three-measurement sum.
"""

# ╔═╡ 62fa6097-336b-4ee3-aa2e-7536ac7e66a8
# the purist way
function rolling_counter(depths)
	rolling_depths = []
	current_depths = []

	for depth in depths
		push!(current_depths, depth)

		# in case we still haven't filled our current_depths array
		length(current_depths) < 3 && continue
		push!(rolling_depths, sum(current_depths))
		popfirst!(current_depths)
	end

	sum((rolling_depths[2:end] .- rolling_depths[1:end-1]) .> 0)
end

# ╔═╡ 60bbbc71-0b56-4123-8a28-2fb83368e826
rolling_counter(depths)

# ╔═╡ 82dd70b8-7f72-445f-bd4d-78225ddaa9f9
rolling_depths = sum.(IterTools.partition(depths, 3, 1))

# ╔═╡ 6ca1b841-c5f9-4dc7-8ef5-521cadfdf83b
sum((rolling_depths[2:end] .- rolling_depths[1:end-1]) .> 0)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"

[compat]
IterTools = "~1.4.0"
"""

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

julia_version = "1.7.0"
manifest_format = "2.0"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.IterTools]]
git-tree-sha1 = "fa6287a4469f5e048d763df38279ee729fbd44e5"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.4.0"

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

# ╔═╡ Cell order:
# ╠═a9ed5ba5-9f59-4bbf-b35c-4a5a36d0b6e0
# ╠═748af037-7d04-4d73-9aaa-575c48fac5f3
# ╠═844a02ca-0744-4523-99a8-37918bff34b3
# ╠═6467bbda-bf07-4ee0-81a3-ae4511c6e8eb
# ╠═3407599f-b7f6-4b02-a36b-5792e4111be6
# ╠═6da7f6d6-755b-43cf-aac2-b99299222955
# ╠═62fa6097-336b-4ee3-aa2e-7536ac7e66a8
# ╠═60bbbc71-0b56-4123-8a28-2fb83368e826
# ╠═291ffe60-5c63-415a-a5f7-1e40e59e2cf9
# ╠═82dd70b8-7f72-445f-bd4d-78225ddaa9f9
# ╠═6ca1b841-c5f9-4dc7-8ef5-521cadfdf83b
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
abstract type ActivationFunction end

struct Id <: ActivationFunction end

(::Id)(x) = x

struct ReLU <: ActivationFunction end

(::ReLU)(x) = max.(x, zero(eltype(x)))

struct Sigmoid <: ActivationFunction end

(::Sigmoid)(x) = @. 1 / (1 + exp(-x))

struct Tanh <: ActivationFunction end

(::Tanh)(x) = tanh.(x)


text:
using RecipesBase

export partitionplot

#-----------------------------------------------------------------# partitionplot
@userplot struct PartitionPlot
    args
end

getvalue(x) = x
getvalue(x::DataValues.DataValue) = get(x)

indextype(x::Type) = x
indextype(x::Type{DataValues.DataValue{T}}) where {T} = T

@recipe function f(o::PartitionPlot; nparts = 100, stat = Extrema(), by = nothing, dropmissing = false)
    t = o.args[1]
    sel_x = o.args[2]
    if length(o.args) == 3
        sel_y = o.args[3]
        T = indextype.(fieldtype(eltype(t), IndexedTables.colindex(t, sel_x)))
        s = dropmissing ?
            FTSeries(IndexedPartition(T, stat, nparts); filter = x->all(!isnull, x), transform = x->getvalue.(x)) :
            FTSeries(IndexedPartition(T, stat, nparts))
        if by == nothing
            reduce(s, t; select = (sel_x, sel_y))
        else
            out = collect(groupreduce(s, t, by; select = (sel_x, sel_y)))
            for i in 1:length(out)
                @series begin
                    label --> OnlineStats.name(stat,false,false) * " of $sel_x ($(out[i][1]))"
                    out[i][2]
                end
            end
        end
    elseif length(o.args) == 2
        s = dropmissing ?
            FTSeries(Partition(stat, nparts); filter = !isnull, transform = getvalue) :
            FTSeries(Partition(stat, nparts))
        if by == nothing
            reduce(s, t; select = sel_x)
        else
            out = groupreduce(s, t, by; select = sel_x)
            for i in 1:length(out)
                @series begin
                    label --> OnlineStats.name(stat,false,false) * " of $sel_x ($(out[i][1]))"
                    out[i][2]
                end
            end
        end
    end
end

"""
    partitionplot(table, y;    stat=Extrema(), nparts=100, by=nothing, dropmissing=false)
    partitionplot(table, x, y; stat=Extrema(), nparts=100, by=nothing, dropmissing=false)

Plot a summary of variable `y` against `x` (`1:length(y)` if not specified).  Using `nparts`
approximately-equal sections along the x-axis, the data in `y` over each section is
summarized by `stat`.
"""
partitionplot


text:
@testset "Test MPS register initialisation" begin
    @test begin
        qreg = MPSQuantumRegister{Int}(3, "000")
        to_str(qreg) == "(1|000>)"
    end
    @test begin
        qreg = MPSQuantumRegister{Int}(3, "100")
        to_str(qreg) == "(1|100>)"
    end
    @test begin
        qreg = MPSQuantumRegister{Int}(3, "001")
        to_str(qreg) == "(1|001>)"
    end
end

@testset "Test convert MPS register to FullStateRegister" begin
    @test begin
        qreg = MPSQuantumRegister{Int}(3, "100")
        qreg_ref = FullStateQuantumRegister{Int}(3, "100")
        convert(FullStateQuantumRegister{Int}, qreg) == qreg_ref
    end
    @test begin
        qreg = MPSQuantumRegister{Int}(3, "100")
        qreg_ref = FullStateQuantumRegister{ComplexF64}(3, "100")
        convert(FullStateQuantumRegister{ComplexF64}, qreg) == qreg_ref
    end
end

@testset "Test applying 1 qubit gates on mps quantum register" begin
    @test begin
        qreg = MPSQuantumRegister{ComplexF64}(3, "000")
        apply_1qubit!(qreg, Gates.x, 3)
        to_str(qreg) == to_str(FullStateQuantumRegister{ComplexF64}(3, "001"))
    end
end

@testset "Test conversion from two qubit gate to mpo and back" begin
    @test begin
        mpo = two_qubit_gate_to_mpo(Gates.cx)
        gate = mpo_to_two_qubit_gate(mpo)
        gate ≈ Gates.cx
    end
    @test begin
        mpo = two_qubit_gate_to_mpo(Gates.rk(4))
        gate = mpo_to_two_qubit_gate(mpo)
        gate ≈ Gates.rk(4)
    end
    @test begin
        gate = rand(ComplexF64, 4, 4)
        gate = gate + gate'
        mpo = two_qubit_gate_to_mpo(gate)
        gate_from_mpo = mpo_to_two_qubit_gate(mpo)
        gate ≈ gate_from_mpo
    end
end

@testset "Test applying 2 qubit gates on mps quantum register" begin
    @test begin
        qreg = MPSQuantumRegister{ComplexF64}(3, "100")
        apply_2qubit!(qreg, Gates.cx, 1, 2)
        convert(FullStateQuantumRegister{ComplexF64}, qreg) ≈ FullStateQuantumRegister{ComplexF64}(3, "110")
    end
    @test begin
        qreg = MPSQuantumRegister{ComplexF64}(3, "100")
        apply_2qubit!(qreg, Gates.cx, 2, 1)
        convert(FullStateQuantumRegister{ComplexF64}, qreg) ≈ FullStateQuantumRegister{ComplexF64}(3, "100")
    end
end

@testset "Test MPS compression on full MPS register" begin
    @test begin
        qreg = MPSQuantumRegister{ComplexF64}(3, "000")
        apply_1qubit!(qreg, Gates.h, 1)
        apply_2qubit!(qreg, Gates.cx, 1, 2)
        apply_2qubit!(qreg, Gates.cx, 2, 3)
        qreg_state = convert(FullStateQuantumRegister{ComplexF64}, qreg)
        compress!(qreg)
        convert(FullStateQuantumRegister{ComplexF64}, qreg) ≈ qreg_state
    end

    @testset "On larger system with random prep circuit" begin
        qreg_mps = MPSQuantumRegister{ComplexF64}(10, "0"^10)
        qreg_full = FullStateQuantumRegister{ComplexF64}(10, "0"^10)
        qc = simple_state_prep_circuit(10, 8)
        execute!(qc, qreg_mps)
        execute!(qc, qreg_full)

        @test begin
            # check that before compression the full state equals mps state
            qreg_full ≈ convert(FullStateQuantumRegister{ComplexF64}, qreg_mps)
        end

        @test begin
            # check states still equal after compression
            compress!(qreg_mps)
            qreg_full ≈ convert(FullStateQuantumRegister{ComplexF64}, qreg_mps)
        end
    end
end

@testset "Test execute function with compression" begin
    qreg1 = MPSQuantumRegister{ComplexF64}(3, "000")
    qreg2 = MPSQuantumRegister{ComplexF64}(3, "000")
    qc = ghz_circuit(3)
    # add some additional unecessary gates that should cancel out
    push!(qc.ops, Gates.GateOperator([1, 2], Gates.cx))
    push!(qc.ops, Gates.GateOperator([1, 2], Gates.cx))
    push!(qc.ops, Gates.GateOperator([2, 3], Gates.cx))
    push!(qc.ops, Gates.GateOperator([2, 3], Gates.cx))
    push!(qc.ops, Gates.GateOperator([1, 3], Gates.cx))
    push!(qc.ops, Gates.GateOperator([1, 3], Gates.cx))
    execute!(qc, qreg1) # with no compression
    execute!(qc, qreg2, 1) # with compression after each gate application

    @test begin
        # test two registers match
        convert(FullStateQuantumRegister{ComplexF64}, qreg1) ≈
        convert(FullStateQuantumRegister{ComplexF64}, qreg2)
    end

    @test begin
        # test that bond dimenions are smaller
        all(size(qreg1.state[2])[2:end] .> size(qreg2.state[2])[2:end])
    end
end


text:
using Glob
using FileIO, JLD2

export ItemCollection, search, details

"""
    Item
An Item in the experiment. Contains the parsed arguments.
"""
struct Item
    folder_str::String
    parsed_args::Dict
    hash_keys
end

function Item(settings_file::AbstractString)
    dict = FileIO.load(settings_file)
    if "settings_dict" ∈ keys(dict)
        dict = dict["settings_dict"]
    end
    parsed_args = dict["parsed_args"]
    used_keys = dict["used_keys"]
    return Item(settings_file, parsed_args, used_keys)
end

"""
    ItemCollection
A collection of items. Mostly helpful, but not really used yet.
"""
mutable struct ItemCollection
    items::Array{Item,1}
end

function ItemCollection(dir::AbstractString; settings_file="settings.jld2")
    dir_list = glob(joinpath(dir, "*", settings_file))
    items = Array{Item,1}()
    for p in dir_list
        append!(items, [Item(p)])
    end

    return ItemCollection(items)
end

"""
    search

Search for specific entries, or a number of entries.
"""
function search(itemCollection::ItemCollection, search_dict)

    dict_keys = keys(search_dict)
    found_items = Array{Item, 1}()
    hash_codes = Array{UInt64, 1}()
    save_dirs = Array{String, 1}()
    for (item_idx, item) in enumerate(itemCollection.items)
        if search_dict == filter(k->((k[1] in dict_keys)), item.parsed_args)
            push!(found_items, item)
            push!(hash_codes, item.parsed_args["_HASH"])
            push!(save_dirs, item.parsed_args["_SAVE"])
        end
    end
    return hash_codes, save_dirs, found_items

end

search(dir::AbstractString, search_dict; settings_file="settings.jld2") =
    search(ItemCollection(dir; settings_file=settings_file), search_dict)

"""
    details

get details of the pointed directory
"""
function details(itemCollection::ItemCollection)
    # itemCollection = ItemCollection(dir; settings_file=settings_file)

    for item in itemCollection.items
        args = ["$(k):$(item.parsed_args[k])" for k in item.hash_keys]
        println(item.parsed_args["_SAVE"], " ", join(args, " "))
    end
end

details(dir::AbstractString; settings_file="settings.jld2") =
    details(ItemCollection(dir; settings_file=settings_file))

import Base.-
"""
    -(l::Dict{K,T}, r::Dict{K,T}) where {K, T}

Get the difference between two dictionaries. Helper function for diff
"""
function -(l::Dict{K, V}, r::Dict{K, V}) where {K, V}
    ret = Dict{K, V}()
    for key in keys(l)
        if l[key] != r[key]
            ret[key] = (l[key], r[key])
        end
    end
    return ret
end

-(l::Item, r::Item) = l.parsed_args - r.parsed_args

"""
    diff

get difference of the list of items.
"""
function Base.diff(items::Array{Reproduce.Item, 1};
              exclude_keys::Union{Array{String,1}, Array{Symbol,1}} = Array{String, 1}(),
              exclude_parse_values::Bool=true)
    kt = keytype(items[1].parsed_args)
    exclude_keys = kt.(exclude_keys)
    if exclude_parse_values == true
        exclude_keys = [exclude_keys; kt.([HASH_KEY, SAVE_NAME_KEY, GIT_INFO_KEY])]
    end
    diff_parsed = Dict{kt, Array}()
    for item in items
        tmp_dict = items[1] - item
        for key in filter((k)->k ∉ exclude_keys, keys(tmp_dict))
            if key ∉ keys(diff_parsed)
                diff_parsed[key] = Array{typeof(items[1].parsed_args[key]), 1}()
            end
            if tmp_dict[key][1] ∉ diff_parsed[key]
                push!(diff_parsed[key], tmp_dict[key][1])
            end
            if tmp_dict[key][2] ∉ diff_parsed[key]
                push!(diff_parsed[key], tmp_dict[key][2])
            end
        end
    end
    for key in keys(diff_parsed)
        sort!(diff_parsed[key])
    end
    return diff_parsed
end

Base.diff(itemCollection::ItemCollection; kwargs...) =
    Base.diff(itemCollection.items; kwargs...)

Base.diff(base_dir::AbstractString; settings_file="settings.jld2", kwargs...) =
    Base.diff(ItemCollection(base_dir; settings_file=settings_file); kwargs...)


text:
using MarketData

facts("Moving Averages") do

    context("sma") do
        @fact sma(cl, 10).values[1]      --> roughly(98.782, atol=.01)  # TTR value 98.782
        @fact sma(cl, 10).values[2]      --> roughly(97.982, atol=.01)  # TTR value 97.982
        @fact sma(cl, 10).values[3]      --> roughly(98.388, atol=.01)  # TTR value 98.388
        @fact sma(cl, 10).values[490]    --> roughly(21.266, atol=.01)  # TTR value 21.266
        @fact sma(cl, 10).values[491]    --> roughly(21.417, atol=.01)  # TTR value 21.417
        @fact sma(cl, 10).timestamp[491] --> Date(2001,12,31)
    end

    context("ema") do
        @fact ema(cl, 10).values[1]      --> roughly(98.78200, atol=.01)   # TTR value 98.78200
        @fact ema(cl, 10).values[2]      --> roughly(99.71982, atol=.01)   # TTR value 99.71982
        @fact ema(cl, 10).values[3]      --> roughly(100.96349, atol=.01)  # TTR value 100.96349
        @fact ema(cl, 10).values[490]    --> roughly(21.58580, atol=.01)   # TTR value 21.58580
        @fact ema(cl, 10).values[491]    --> roughly(21.64293, atol=.01)   # TTR value 21.64293
        @fact ema(cl, 10).timestamp[491] --> Date(2001,12,31)
    end

    context("ema wilder true") do
        @fact ema(cl, 10, wilder=true).values[1]      --> roughly(98.7820, atol=.01)  # TTR value 98.7820
        @fact ema(cl, 10, wilder=true).values[2]      --> roughly(99.2978, atol=.01)  # TTR value 99.2978
        @fact ema(cl, 10, wilder=true).values[3]      --> roughly(100.0240, atol=.01)  # TTR value 100.0240
        @fact ema(cl, 10, wilder=true).values[490]    --> roughly(21.34556, atol=.01)  # TTR value 21.34556
        @fact ema(cl, 10, wilder=true).values[491]    --> roughly(21.40100, atol=.01)  # TTR value 21.40100
        @fact ema(cl, 10, wilder=true).timestamp[491] --> Date(2001,12,31)
    end

    context("Array dispatch is correct too") do
        @fact sma(cl.values, 10)[1]                --> roughly(98.782, atol=.01)   # same values as above TimeArray dispatch
        @fact sma(cl.values, 10)[2]                --> roughly(97.982, atol=.01)
        @fact sma(cl.values, 10)[3]                --> roughly(98.388, atol=.01)
        @fact sma(cl.values, 10)[490]              --> roughly(21.266, atol=.01)
        @fact sma(cl.values, 10)[491]              --> roughly(21.417, atol=.01)
        @fact ema(cl.values, 10)[1]                --> roughly(98.78200, atol=.01)
        @fact ema(cl.values, 10)[2]                --> roughly(99.71982, atol=.01)
        @fact ema(cl.values, 10)[3]                --> roughly(100.96349, atol=.01)
        @fact ema(cl.values, 10)[490]              --> roughly(21.58580, atol=.01)
        @fact ema(cl.values, 10)[491]              --> roughly(21.64293, atol=.01)
        @fact ema(cl.values, 10, wilder=true)[1]   --> roughly(98.7820, atol=.01)
        @fact ema(cl.values, 10, wilder=true)[2]   --> roughly(99.2978, atol=.01)
        @fact ema(cl.values, 10, wilder=true)[3]   --> roughly(100.0240, atol=.01)
        @fact ema(cl.values, 10, wilder=true)[490] --> roughly(21.34556, atol=.01)
        @fact ema(cl.values, 10, wilder=true)[491] --> roughly(21.40100, atol=.01)
    end
end


text:
# L2,1 norm/Sum of norms of columns or rows (times a constant)

export NormL21

"""
**Sum-of-``L_2`` norms**

    NormL21(λ=1, dim=1)

Returns the function
```math
f(X) = λ⋅∑_i\\|x_i\\|
```
for a nonnegative `λ`, where ``x_i`` is the ``i``-th column of ``X`` if `dim == 1`, and the ``i``-th row of ``X`` if `dim == 2`.
In words, it is the sum of the Euclidean norms of the columns or rows.
"""
struct NormL21{R <: Real, I <: Integer} <: ProximableFunction
    lambda::R
    dim::I
    function NormL21{R,I}(lambda::R, dim::I) where {R <: Real, I <: Integer}
        if lambda < 0
            error("parameter λ must be nonnegative")
        else
            new(lambda, dim)
        end
    end
end

is_convex(f::NormL21) = true

NormL21(lambda::R=1, dim::I=1) where {R <: Real, I <: Integer} = NormL21{R, I}(lambda, dim)

function (f::NormL21)(X::AbstractArray{T, 2}) where {R, T <: RealOrComplex{R}}
    nslice = R(0)
    n21X = R(0)
    if f.dim == 1
        for j = 1:size(X, 2)
            nslice = R(0)
            for i = 1:size(X, 1)
                nslice += abs(X[i, j])^2
            end
            n21X += sqrt(nslice)
        end
    elseif f.dim == 2
        for i = 1:size(X, 1)
            nslice = R(0)
            for j = 1:size(X, 2)
                nslice += abs(X[i, j])^2
            end
            n21X += sqrt(nslice)
        end
    end
    return f.lambda * n21X
end

function prox!(Y::AbstractArray{T, 2}, f::NormL21, X::AbstractArray{T, 2}, gamma::Real=1) where {R, T <: RealOrComplex{R}}
    gl = gamma * f.lambda
    nslice = R(0)
    n21X = R(0)
    if f.dim == 1
        for j = 1:size(X, 2)
            nslice = R(0)
            for i = 1:size(X, 1)
                nslice += abs(X[i, j])^2
            end
            nslice = sqrt(nslice)
            scal = 1 - gl / nslice
            scal = scal <= 0 ? R(0) : scal
            for i = 1:size(X, 1)
                Y[i, j] = scal * X[i, j]
            end
            n21X += scal * nslice
        end
    elseif f.dim == 2
        for i = 1:size(X, 1)
            nslice = R(0)
            for j = 1:size(X, 2)
                nslice += abs(X[i, j])^2
            end
            nslice = sqrt(nslice)
            scal = 1-gl/nslice
            scal = scal <= 0 ? R(0) : scal
            for j = 1:size(X, 2)
                Y[i, j] = scal * X[i, j]
            end
            n21X += scal * nslice
        end
    end
    return f.lambda * n21X
end

fun_name(f::NormL21) = "sum of Euclidean norms"
fun_dom(f::NormL21) = "AbstractArray{Real,2}, AbstractArray{Complex, 2}"
fun_expr(f::NormL21) = "x ↦ λsum(||x_i||)"
fun_params(f::NormL21) = "λ = $(f.lambda), dim = $(f.dim)"

function prox_naive(f::NormL21, X::AbstractArray{T,2}, gamma::Real=1.0) where T <: RealOrComplex
    Y = max.(0, 1 .- f.lambda * gamma ./ sqrt.(sum(abs.(X).^2, dims=f.dim))) .* X
    return Y, f.lambda * sum(sqrt.(sum(abs.(Y).^2, dims=f.dim)))
end


text:
# promotion methods for CircuitParams
for p in (:Sparams, :Yparams, :Zparams, :ABCDparams)
    @eval promote_rule(::Type{($p){T}}, ::Type{($p){S}}) where {T<:Real, S<:Real} =
        ($p){promote_type(T, S)}
end

# Conversion between general nPort parameters:
for p in (:Sparams, :Yparams, :Zparams, :ABCDparams)
    @eval convert(::Type{$p}, s::NetworkData{S, $p{T}}) where
        {S<:Real, T<:Real} = s
end

convert(::Type{Zparams}, s::NetworkData{S, Sparams{T}}) where
    {S<:Real, T<:Real} = _S_to_Z(s)
convert(::Type{Yparams}, s::NetworkData{S, Sparams{T}}) where
    {S<:Real, T<:Real} = _S_to_Y(s)
convert(::Type{Sparams}, z::NetworkData{S, Zparams{T}}) where
    {S<:Real, T<:Real} = _Z_to_S(z)
convert(::Type{Sparams}, y::NetworkData{S, Yparams{T}}) where
    {S<:Real, T<:Real} = _Y_to_S(y)
convert(::Type{Zparams}, y::NetworkData{S, Yparams{T}}) where
    {S<:Real, T<:Real} =
    NetworkData(y.ports, y.frequency,
    [Zparams(inv(y.params[n].data)) for n in 1:y.nPoint])
convert(::Type{Yparams}, z::NetworkData{S, Zparams{T}}) where
    {S<:Real, T<:Real} =
    NetworkData(z.ports, z.frequency,
    [Yparams(inv(z.params[n].data)) for n in 1:z.nPoint])

"""
Helper matrices related to reference impedances
"""
Z_ref(ports::Vector{Port}) = diagm(impedances(ports), 0)
Z_ref(D::NetworkData) = Z_ref(D.ports)
G_ref(ports::Vector{Port}) = diagm(1./sqrt.(abs.(impedances(ports))), 0)
G_ref(D::NetworkData) = G_ref(D.ports)

"""
Conversion from S-parameters to Z-parameters
"""
function _S_to_Z(s::NetworkData{S, Sparams{T}}) where {S<:Real, T<:Real}
    E = eye(Complex{T}, s.nPort)
    _Z_ref, _G_ref = Z_ref(s), G_ref(s)
    params = [Zparams(inv(_G_ref) * inv(E - s.params[n].data) *
        (E + s.params[n].data) * _Z_ref * _G_ref) for n in 1:s.nPoint]
    return NetworkData(s.ports, s.frequency, params)
end

"""
Conversion from S-parameters to Y-parameters
"""
function _S_to_Y(s::NetworkData{S, Sparams{T}}) where {S<:Real, T<:Real}
    E = eye(Complex{T}, s.nPort)
    _Z_ref, _G_ref = Z_ref(s), G_ref(s)
    params = [Yparams(inv(_G_ref) * inv(_Z_ref) *
        inv(E + s.params[n].data) * (E - s.params[n].data) *
        _G_ref) for n in 1:s.nPoint]
    return NetworkData(s.ports, s.frequency, params)
end

"""
Conversion from Z-parameters to S-parameters
"""
function _Z_to_S(z::NetworkData{S, Zparams{T}}) where {S<:Real, T<:Real}
    _Z_ref, _G_ref = Z_ref(z), G_ref(z)
    params = [Sparams(_G_ref * (z.params[n].data - _Z_ref) *
        inv(z.params[n].data + _Z_ref) * inv(_G_ref)) for n in 1:z.nPoint]
    return NetworkData(z.ports, z.frequency, params)
end

"""
Conversion from Y-parameters to S-parameters
"""
function _Y_to_S(y::NetworkData{S, Yparams{T}}) where {S<:Real, T<:Real}
    E = eye(Complex{T}, y.nPort)
    _Z_ref, _G_ref = Z_ref(y), G_ref(y)
    params = [Sparams(_G_ref * (E - _Z_ref * y.params[n].data) *
        inv(E + _Z_ref * y.params[n].data) * inv(_G_ref)) for n in 1:y.nPoint]
    return NetworkData(y.ports, y.frequency, params)
end

# for two-port parameters
convert(::Type{Sparams}, abcd::NetworkData{S, ABCDparams{T}}) where
    {S<:Real, T<:Real} = _ABCD_to_S(abcd)
convert(::Type{ABCDparams}, s::NetworkData{S, Sparams{T}}) where
    {S<:Real, T<:Real} = _S_to_ABCD(s)

"""
Conversion from ABCD-parameters to S-parameters
"""
function _ABCD_to_S(abcd::NetworkData{S, ABCDparams{T}}) where {S<:Real, T<:Real}
    Z0 = impedances(abcd)[1]
    s = [begin
        (A, B, C, D) = (abcd[(1, 1), n], abcd[(1, 2), n], abcd[(2, 1), n], abcd[(2, 2), n])
        S11 = (A + B / Z0 - C * Z0 - D) / (A + B / Z0 + C * Z0 + D)
        S12 = 2 * (A * D - B * C) / (A + B / Z0 + C * Z0 + D)
        S21 = 2 / (A + B / Z0 + C * Z0 + D)
        S22 = (-A + B / Z0 - C * Z0 + D) / (A + B / Z0 + C * Z0 + D)
        Sparams([S11 S12; S21 S22])
        end for n in 1:abcd.nPoint]
    return NetworkData(abcd.ports, abcd.frequency, s)
end

"""
Conversion from S-parameters to ABCD-parameters
"""
function _S_to_ABCD(s::NetworkData{S, Sparams{T}}) where {S<:Real, T<:Real}
    if s.nPort != 2
        error("Error: ABCD-parameters are defined only for 2-port networks")
    end
    if is_uniform(s) == false
        error("Error: The port impedances of a ABCDparams must be uniform")
    end
    Z0 = impedances(s)[1]
    abcd = [begin
    (S11, S12, S21, S22) = (s[(1, 1), n], s[(1, 2), n], s[(2, 1), n], s[(2, 2), n])
    A = ((1 + S11) * (1 - S22) + S12 * S21) / (2 * S21)
    B = Z0 * ((1 + S11) * (1 + S22) - S12 * S21) / (2 * S21)
    C = 1 / Z0 * ((1 - S11) * (1 - S22) - S12 * S21) / (2 * S21)
    D = ((1 - S11) * (1 + S22) + S12 * S21) / (2 * S21)
    ABCDparams([A B; C D])
    end for n in 1:s.nPoint]  # converting S-parameters into ABCD parameters
    return NetworkData(s.ports, s.frequency, abcd)
end