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

using BAT
using Test

using Distributions, StatsBase, ValueShapes
using LinearAlgebra
using ForwardDiff, Zygote, DistributionsAD

@testset "zygote" begin
    posterior = ZygoteAD()|BAT.example_posterior_with_dirichlet()

    target, trafo = bat_transform(PriorToGaussian(), posterior)
    v = bat_initval(target).result
    f = logdensityof(target)
    gradlogp_zg = valgradof(f)
    gradlogp_fd = valgradof(logdensityof(bat_transform(PriorToGaussian(), ForwardDiffAD() | BAT.example_posterior_with_dirichlet()).result))
    @test @inferred(gradlogp_zg(v))[1] ≈ gradlogp_fd(v)[1]
    @test @inferred(gradlogp_zg(v))[2] ≈ gradlogp_fd(v)[2]

    utarget, utrafo = bat_transform(PriorToGaussian(), unshaped(posterior))
    x = bat_initval(utarget).result
    @test @inferred(Zygote.jacobian(inverse(utrafo), x)[1]) ≈ ForwardDiff.jacobian(inverse(utrafo), x)

    hmc_sampling_alg = MCMCSampling(
        mcalg = HamiltonianMC(),
        nchains = 2,
        nsteps = 100,
        init = MCMCChainPoolInit(init_tries_per_chain = 2, nsteps_init = 5),
        burnin = MCMCMultiCycleBurnin(nsteps_per_cycle = 100, max_ncycles = 1),
        strict = false
    )

    @test bat_sample(posterior, hmc_sampling_alg).result isa DensitySampleVector
end


text:
__precompile__(true)

module GGPlots

using Reexport
@reexport using StatPlots
@reexport using DataFrames
import Base: +

export
    ggplot,
    geom_point,
    geom_smooth

# -------------------------------------------------------------------------
# geoms

abstract Geom

const _point_keys = KW(:color => :marker_z)
immutable Point <: Geom
    kw::KW
end
function Point(; kw...)
    g = Point(KW(:seriestype => :scatter))
    for (k,v) in kw
        g.kw[get(_point_keys, k, k)] = v
    end
    g
end
const geom_point = Point

const _smooth_keys = KW(:method => :smooth)
const _smooth_vals = KW(:lm => true)
immutable Smooth <: Geom
    kw::KW
end
function Smooth(; kw...)
    g = Smooth(KW())
    for (k,v) in kw
        g.kw[get(_smooth_keys, k, k)] = get(_smooth_vals, v, v)
    end
    g
end
const geom_smooth = Smooth

# -------------------------------------------------------------------------
# ggplot

immutable GGPlot
    df::AbstractDataFrame
    args::Tuple
    kw::KW
end
function ggplot(df::AbstractDataFrame, args...; kw...)
    GGPlot(df, args, KW(kw))
end

function +(gg::GGPlot, geom::Geom)
    merge!(gg.kw, geom.kw)
    gg
end

StatPlots.plot(gg::GGPlot) = plot(gg.df, gg.args...; gg.kw...)
Base.display(gg::GGPlot) = display(plot(gg))

end # module

ext: jl


text:
using RigidBodyDynamicsDiff
using RigidBodyDynamics
using ForwardDiff
using DiffResults
using Random
using Test

using RigidBodyDynamics: velocity_to_configuration_derivative_jacobian
using RigidBodyDynamics: configuration_derivative_to_velocity_jacobian

include(joinpath(@__DIR__, "..", "test", "forwarddiff_compatible.jl"))

struct MyTag end

@testset "differential (nthreads = $(Threads.nthreads()))" begin
    Random.seed!(1)

    T = Float64
    urdf = joinpath(dirname(pathof(RigidBodyDynamics)), "..", "test", "urdf", "atlas.urdf")
    mechanism = parse_urdf(urdf, scalar_type=T, root_joint_type=QuaternionFloating{T}())
    state = MechanismState{T}(mechanism)
    result = DynamicsResult{T}(mechanism)
    rand!(state)
    q = configuration(state)
    v = velocity(state)
    v̇ = rand!(similar(velocity(state)))
    τ = rand!(similar(velocity(state)))
    nv = num_velocities(mechanism)
    statecache = StateCache(mechanism)
    resultcache = DynamicsResultCache(mechanism)
    normalize = Ref(false)

    configs = []
    push!(configs, (
        "mass_matrix",
        mass_matrix_differential!,
        forwarddiff_compatible(mass_matrix!, statecache, resultcache, normalize),
        zeros(nv * nv)
    ))
    push!(configs, (
        "dynamics_bias",
        dynamics_bias_differential!,
        forwarddiff_compatible(dynamics_bias!, statecache, resultcache, v, normalize),
        zeros(nv)
    ))
    push!(configs, (
        "inverse_dynamics",
        (dest, state, cache) -> inverse_dynamics_differential!(dest, state, v̇, cache),
        forwarddiff_compatible(inverse_dynamics!, statecache, resultcache, v, v̇, normalize),
        zeros(nv)
    ))
    push!(configs, (
        "dynamics",
        (dest, state, cache) -> dynamics_differential!(dest, result, state, τ, cache),
        forwarddiff_compatible(dynamics!, statecache, resultcache, v, τ, normalize),
        zeros(nv)
    ))

    for (name, differentialfun, fdfun, out) in configs
        @testset "$name" begin
            diffcache = DifferentialCache{MyTag}(mechanism)

            fdjacresult = DiffResults.JacobianResult(out, q)
            fdconfig = ForwardDiff.JacobianConfig(fdfun, out, q, ForwardDiff.Chunk(1))
            differential = fill(T(NaN), length(out), nv)
            differentialfun(differential, state, diffcache)

            # Without normalization:
            normalize[] = false
            ForwardDiff.jacobian!(fdjacresult, fdfun, out, q, fdconfig)
            fdjac = DiffResults.jacobian(fdjacresult)
            velocity_to_configuration_derivative_jacobian(state)
            @test differential ≈ fdjac * velocity_to_configuration_derivative_jacobian(state) atol = 1e-8

            # Without normalization, the normal component in dM is zero (since we're only working
            # in the tangent plane), but the ForwardDiff version *does* have a normal component,
            # so the following is not true:
            #   @test dM * configuration_derivative_to_velocity_jacobian(state) ≈ fdjac atol = 1e-10
            # But with normalization:

            normalize[] = true
            ForwardDiff.jacobian!(fdjacresult, fdfun, out, q, fdconfig)
            fdjac = DiffResults.jacobian(fdjacresult)
            @test differential * configuration_derivative_to_velocity_jacobian(state) ≈ fdjac atol = 1e-8
            @test differential ≈ fdjac * velocity_to_configuration_derivative_jacobian(state) atol = 1e-8
        end
    end
end

@testset "benchmarks" begin
    include(joinpath(@__DIR__, "..", "perf", "runbenchmarks.jl"))
end


text:
""" Seconds since the first sample """
const TimeStamp = Second

# TODO: Support custom fields (wrap NamedTuple)
struct Process
    timestamp::TimeStamp
    pid::Int
    user::Symbol
    priority::Int                                               # PR
    nice::Int                                                   # NI
    virtual_memory_bytes::Int                                   # VIRT
    resident_memory_bytes::Int                                  # RES
    shared_memory_bytes::Int                                    # SHR
    state::Symbol                                               # S
    cpu_usage::Float64                                          # %CPU
    memory_usage::Float64                                       # %MEM
    cpu_time::Float64                                           # TIME+
    command::String
end

struct TimeOfDay
    hour::Int8
    minute::Int8
    second::Int8
end

Base.:(==)(a::TimeOfDay, b::TimeOfDay) =
    a.hour == b.hour && a.minute == b.minute && a.second == b.second

function Base.:<(a::TimeOfDay, b::TimeOfDay)
    a.hour < b.hour && return true
    a.hour > b.hour && return false
    a.minute < b.minute && return true
    a.minute > b.minute && return false
    return a.second < b.second
end

as_seconds(t::TimeOfDay) = 60 * (60 * t.hour + t.minute) + t.second

struct SamplesIterator{IO}
    io::IO
end

mutable struct ProcessIterator{IO}
    timestamp::TimeStamp
    io::IO
    done::Bool
end

struct Sample
    timestamp::TimeStamp
    timeofday::TimeOfDay
    processes::ProcessIterator
end

Base.eltype(::Type{<:SamplesIterator}) = Sample
Base.eltype(::Type{<:ProcessIterator}) = Process

Base.IteratorSize(::Type{<:SamplesIterator}) = Base.SizeUnknown()
Base.IteratorSize(::Type{<:ProcessIterator}) = Base.SizeUnknown()

function Base.iterate(it::SamplesIterator, state = nothing)
    io = it.io

    if state !== nothing
        if !state.pit.done
            while !all(isspace, readline(io))
            end
            state.pit.done = true
        end
    end

    ln = readline(io)
    if all(isspace, ln) && eof(io)
        return nothing
    end

    m = match(r"^top *- *([0-9]{2}):([0-9]{2}):([0-9]{2}) ", ln)
    if m === nothing
        error("unexpected first line of summary: ", ln)
    end
    timeofday = TimeOfDay(parse(Int8, m[1]), parse(Int8, m[2]), parse(Int8, m[3]))

    if state === nothing
        seconds = 0
    else
        seconds = state.seconds
        if timeofday < state.timeofday
            # Assuming the delay (as in `top -d $delay`) is < 1 day
            # TODO: verify this using the uptime information
            rest = 24 * 60 * 60 - as_seconds(state.timeofday)
            @assert rest >= 0
            seconds += rest
            seconds += as_seconds(timeofday)
        else
            seconds += as_seconds(timeofday) - as_seconds(state.timeofday)
        end
    end

    while !all(isspace, readline(io))
    end

    ln = readline(io)
    m = match(
        r"^ +PID +USER +PR +NI +VIRT +RES +SHR +S +%CPU +%MEM +TIME\+ +COMMAND *$",
        ln,
    )
    if m === nothing
        eof(io) && return nothing
        error("unexpected first line of table: ", ln)
    end

    timestamp = Second(seconds)
    pit = ProcessIterator(timestamp, io, false)

    return (Sample(timestamp, timeofday, pit), (; seconds, timeofday, pit))
end

function Base.iterate(it::ProcessIterator, _ignored = nothing)
    timestamp = it.timestamp
    io = it.io
    it.done && return nothing

    ln = readline(io)
    cols = split(ln, limit = 12)
    if length(cols) < 11
        it.done = true
        return nothing
    end

    let (pid, user, pr, ni, virt, res, shr, s, pcpu, pmem, time) = cols
        command = get(cols, 12, "")
        args = (
            tryparse(Int, pid),
            Symbol(user),
            pr == "rt" ? typemax(Int) : tryparse(Int, pr),
            tryparse(Int, ni),
            tryparse_memory(virt),
            tryparse_memory(res),
            tryparse_memory(shr),
            Symbol(s),
            tryparse(Float64, pcpu),
            tryparse(Float64, pmem),
            tryparse_time(time),
            command,
        )
        if all(!isnothing, args)
            return (Process(timestamp, map(something, args)...), nothing)
        end
        cols = join(findall(isnothing, args), ", ", " and ")
        error("unexpected value(s) in column(s) $cols of line: ", ln)
    end
end

function tryparse_memory(str)
    unit = lowercase(str[end])
    endat = lastindex(str) - 1
    if unit == 'k'
        scaling = 2^10
    elseif unit == 'm'
        scaling = 2^20
    elseif unit == 'g'
        scaling = 2^30
    elseif unit == 't'
        scaling = 2^40
    elseif unit == 'p'
        scaling = 2^50
    else
        scaling = 2^10  # KiB
        endat = lastindex(str)
    end
    v = tryparse(Float64, view(str, 1:endat))
    v === nothing && return v
    return trunc(Int, v * scaling)
end

function tryparse_time(str)
    m = match(
        r"""
        ^(?:
            (?<minute>
                [0-9]*
            ):
        )?
        (?<second>
            [0-9]+
            (?:
                \.[0-9]*
            )?
        )$
        """x,
        str,
    )
    m === nothing && return nothing
    return parse(Int, m[:minute]) * 60 + parse(Float64, m[:second])
end


text:
# ---
# title: 240. Search a 2D Matrix II
# id: problem240
# author: woclass
# date: 2020-10-29
# difficulty: Medium
# categories: Binary Search, Divide and Conquer
# link: <https://leetcode.com/problems/search-a-2d-matrix-ii/description/>
# hidden: true
# ---
#
# Write an efficient algorithm that searches for a `target` value in an `m x n`
# integer `matrix`. The `matrix` has the following properties:
#
#   * Integers in each row are sorted in ascending from left to right.
#   * Integers in each column are sorted in ascending from top to bottom.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/24/searchgrid2.jpg)
#
#
#
#     Input: matrix = [[1,4,7,11,15],[2,5,8,12,19],[3,6,9,16,22],[10,13,14,17,24],[18,21,23,26,30]], target = 5
#     Output: true
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/24/searchgrid.jpg)
#
#
#
#     Input: matrix = [[1,4,7,11,15],[2,5,8,12,19],[3,6,9,16,22],[10,13,14,17,24],[18,21,23,26,30]], target = 20
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `m == matrix.length`
#   * `n == matrix[i].length`
#   * `1 <= n, m <= 300`
#   * `-109 <= matix[i][j] <= 109`
#   * All the integers in each row are **sorted** in ascending order.
#   * All the integers in each column are **sorted** in ascending order.
#   * `-109 <= target <= 109`
#
#
## @lc code=start
using LeetCode

function search_matrix(matrix::Matrix{Int}, target::Int)::Bool
    row, col = size(matrix)
    i, j = row, 1

    while (i > 0 && j <= col)
        last = matrix[i, j]

        if (last == target)
            return true
        elseif last > target
            i -= 1
        else
            j += 1
        end
    end

    return false
end
## @lc code=end


text:
@defcomp damages begin
    regions = Index()

    DAMFRAC = Variable(index=[time, regions]) # Damages as % of GDP

    TATM = Parameter(index=[time]) # Increase temperature of atmosphere (degrees C from 1900)
    SLRDAMAGES = Parameter(index=[time, regions])
    a1 = Parameter(index=[regions]) # Damage intercept
    a2 = Parameter(index=[regions]) # Damage quadratic term
    a3 = Parameter(index=[regions]) # Damage exponent

    function run_timestep(p, v, d, t)

        #Define function for DAMFRAC
        for r in d.regions
            v.DAMFRAC[t,r] = (((p.a1[r] * p.TATM[t]) + (p.a2[r] * p.TATM[t]^p.a3[r])) / 100) + (p.SLRDAMAGES[t,r] / 100)
        end
    end
end

max_stars_repo_name: fdennig/RICEupdate

max_issues_repo_name: fdennig/RICEupdate

max_forks_repo_name: fdennig/RICEupdate


text:
# Convenience functions for working with MPI
using MPI

"""
Initialize MPI. Must be called before doing any non-trivial MPI work
(even in the single-process case). Unlike the MPI.Init() function,
this can be called multiple times.
"""
function mpi_ensure_initialized()
    # MPI Thread level 3 means that the environment is multithreaded, but that only
    # one thread will call MPI at once
    # see https://www.open-mpi.org/doc/current/man3/MPI_Init_thread.3.php#toc7
    # TODO look more closely at interaction between MPI and threads
    MPI.Initialized() || MPI.Init_thread(MPI.ThreadLevel(3))
end

"""
Number of processors used in MPI. Can be called without ensuring initialization.
"""
mpi_nprocs(comm=MPI.COMM_WORLD)  = (mpi_ensure_initialized(); MPI.Comm_size(comm))
mpi_master() = (mpi_ensure_initialized(); MPI.Comm_rank(MPI.COMM_WORLD) == 0)

mpi_sum( arr, comm::MPI.Comm)  = MPI.Allreduce( arr, +, comm)
mpi_sum!(arr, comm::MPI.Comm)  = MPI.Allreduce!(arr, +, comm)
mpi_min( arr, comm::MPI.Comm)  = MPI.Allreduce( arr, min, comm)
mpi_min!(arr, comm::MPI.Comm)  = MPI.Allreduce!(arr, min, comm)
mpi_max( arr, comm::MPI.Comm)  = MPI.Allreduce( arr, max, comm)
mpi_max!(arr, comm::MPI.Comm)  = MPI.Allreduce!(arr, max, comm)
mpi_mean(arr, comm::MPI.Comm)  = mpi_sum(arr, comm) ./ mpi_nprocs(comm)
mpi_mean!(arr, comm::MPI.Comm) = (mpi_sum!(arr, comm); arr ./= mpi_nprocs(comm))

"""
Splits an iterator evenly between the processes of `comm` and returns the part handled
by the current process.
"""
function mpi_split_iterator(itr, comm)
    nprocs = mpi_nprocs(comm)
    @assert nprocs <= length(itr)
    split_evenly(itr, nprocs)[1 + MPI.Comm_rank(comm)]  # MPI ranks are 0-based
end


text:
using Test

@testset "running tests for module BellComm.jl" begin

    println("running BellComm.jl unit tests:")
    for test in readdir("./test/BellComm/unit/")
        # run only julia files in test directory
        if occursin(r"^.*\.jl$", test)
            println("./BellComm/unit/$test")
            include("./unit/$test")
        end
    end

    println("running BellComm.jl integration tests:")
    for test in readdir("./test/BellComm/integration/")
        if occursin(r"^.*\.jl$", test)
            println("./BellComm/integration/$test")
            include("./integration/$test")
        end
    end
end


text:
include("iris_tools.jl") # IrisTools

using IrisTools #

# fname # ERROR: UndefVarError: fname not defined
println(species) # ASCIIString["I. setosa","I. versicolor","I. virginica"]

import IrisTools.fname
println(fname) # iris.txt

iris1 = Iris(5.1, 3.5, 1.4, 0.2, "I. setosa")
println(to_string(iris1))
# This I. setosa has measurements:
#     Sepal length: 5.1
#     Sepal width: 3.5
#     Petal length: 1.4
#     Petal width: 0.2

# reloading:
workspace()

include("iris_tools.jl")
using IrisTools

max_stars_repo_path: Section 2/use_iris_tools.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 2/use_iris_tools.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 2/use_iris_tools.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
@testset ExtendedTestSet "Annealing" begin
    @testset ExtendedTestSet "Edge identity" begin
        for anneal_fn in (anneal_const, anneal_cosine, anneal_exp, anneal_linear)
            @test anneal_fn(0, 0.1, 1) ≈ 0.1
            @test anneal_fn(1, 0.1, 1) ≈ 1
        end
    end
    @testset ExtendedTestSet "Exponential annealing undefined from 0" begin
        @test_throws DomainError anneal_exp(0, 0, 1)
    end
end

@testset ExtendedTestSet "Functional metrics" begin
    @test accuracy([1 0; 0 1], [1 0; 0 1]) == 1
    @test accuracy([0.9 0.1; 0.1 0.9], [1 0; 0 1]) == 1

end


text:
module density_demo

include("util.jl")
export
    xlogx,
    neumann_entropy,
    binary_entropy,
    thermal_logZ,
    latent_density,
    sigmoid,
    trdistance,
    fidelity

include("adjacency_map.jl")
export
    lattice_adjacency_map,
    lattice_adjacency_map_horz,
    lattice_adjacency_map_vert,
    lattice_map_1d

include("construct_hamiltonian.jl")
export
    construct_hamiltonian,
    construct_1d_hamiltonian,
    construct_2d_hamiltonian

include("thermal_states.jl")
export
    thermal_state,
    thermal_logZ

include("construct_circuit.jl")
export
    create_cgc

end # module

max_stars_repo_name: cmendl/density-matrix-bloch-qcircuit

max_issues_repo_name: cmendl/density-matrix-bloch-qcircuit

max_forks_repo_name: cmendl/density-matrix-bloch-qcircuit


text:
module WGLMakie

using Hyperscript
using JSServe
using Observables
using AbstractPlotting
using Colors
using ShaderAbstractions
using LinearAlgebra
using GeometryBasics
using ImageMagick
using FreeTypeAbstraction
using StaticArrays

using JSServe: Session
using JSServe: @js_str, onjs, Dependency, App
using JSServe.DOM

using ShaderAbstractions: VertexArray, Buffer, Sampler, AbstractSampler
using ShaderAbstractions: InstancedProgram

import AbstractPlotting.FileIO
using AbstractPlotting: get_texture_atlas, glyph_uv_width!, SceneSpace, Pixel
using AbstractPlotting: attribute_per_char, glyph_uv_width!, layout_text
using AbstractPlotting: apply_transform, transform_func_obs

struct WebGL <: ShaderAbstractions.AbstractContext end
struct WGLBackend <: AbstractPlotting.AbstractBackend end
#["https://unpkg.com/three@0.123.0/build/three.min.js"
const THREE = Dependency(:THREE,
                                 ["https://cdn.jsdelivr.net/gh/mrdoob/three.js/build/three.js"])

const WGL = Dependency(:WGLMakie, [joinpath(@__DIR__, "wglmakie.js")])
const WEBGL = Dependency(:WEBGL, [joinpath(@__DIR__, "WEBGL.js")])

include("three_plot.jl")
include("serialization.jl")
include("events.jl")
include("particles.jl")
include("lines.jl")
include("meshes.jl")
include("imagelike.jl")
include("display.jl")

function activate!()
    b = WGLBackend()
    AbstractPlotting.register_backend!(b)
    AbstractPlotting.current_backend[] = b
    AbstractPlotting.set_glyph_resolution!(AbstractPlotting.Low)
    return
end

const TEXTURE_ATLAS_CHANGED = Ref(false)

function __init__()
    # Activate WGLMakie as backend!
    activate!()
    browser_display = JSServe.BrowserDisplay() in Base.Multimedia.displays
    AbstractPlotting.inline!(!browser_display)
    # We need to update the texture atlas whenever it changes!
    # We do this in three_plot!
    AbstractPlotting.font_render_callback!() do sd, uv
        TEXTURE_ATLAS_CHANGED[] = true
    end
end

for name in names(AbstractPlotting)
    if name !== :Button && name !== :Slider
        @eval import AbstractPlotting: $(name)
        @eval export $(name)
    end
end

end # module


text:
export KooshballTrajectory, kooshballNodes, kooshballDensity, kooshBallTimes

"""
    KooshballTrajectory(numProfiles, numSamplingPerProfile
                  ; TE::Float64=0.0
                  , AQ::Float64=1.e-3
                  , kargs...)

returns a 3d kooshball trajectory.

# Arguments
* `numProfiles::Int64`            - number of profiles
* `numSamplingPerProfile::Int64`  - number of sampling points per profile
* (`TE::Float64=0.0`)             - echo time in s
* (`AQ::Float64=1.e-3`)           - readout duration in s (per profile)
"""
function KooshballTrajectory(::Type{T}, numProfiles, numSamplingPerProfile
                  ; TE=0.0
                  , AQ=1.e-3
                  , kargs...) where T
  nodes = kooshballNodes(T, numProfiles, numSamplingPerProfile; numSlices=numSlices)
  times = kooshballTimes(T, numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ)
  return  Trajectory("Kooshball", nodes, times, TE, AQ, numProfiles, numSamplingPerProfile, 1, false, true)
end

function kooshballNodes(::Type{T}, numProfiles::Int64, numSamplingPerProfile::Int64; kargs...) where T
  nodes = zeros(T, 3, numSamplingPerProfile, numProfiles)
  pos = collect((0:numSamplingPerProfile-1)/numSamplingPerProfile .- 0.5)
  for l = 1:numProfiles
    for k = 1:numSamplingPerProfile
      nodes[1,k,l] = (-1)^l * pos[k]*xax(l,2*numProfiles)
      nodes[2,k,l] = (-1)^l * pos[k]*yax(l,2*numProfiles)
      nodes[3,k,l] = (-1)^l * pos[k]*zax(l,2*numProfiles)
    end
  end
  return reshape(nodes, 3, numSamplingPerProfile*numProfiles)
end

function kooshballTimes(::Type{T}, numProfiles::Int64, numSamplingPerProfile::Int64; TE=0.0, AQ=1.e-3, kargs...) where T
  times = zeros(T, numSamplingPerProfile, numProfiles)
  for l = 1:numProfiles
    for k = 1:numSamplingPerProfile
      times[k,l] = (TE + AQ*(k-1)/numSamplingPerProfile )
    end
  end
  return return vec(times)
end

function kooshballDensity(numSamplingPerProfile::Int64, numProfiles::Int64, numSlices::Int64)
  density = zeros(numSamplingPerProfile, numProfiles)
  pos = collect((0:numSamplingPerProfile-1)/numSamplingPerProfile .- 0.5)
  for l = 1:numProfiles
    for k = 1:numSamplingPerProfile
      density[k,l] = abs(pos[k])
    end
  end
  return vec(density)
end


text:
module DistributionsX

using Random: AbstractRNG, GLOBAL_RNG, rand!, randn!
import Requires, StatsFuns, NNlib

### Utility

Base.eps(x::AbstractArray) = eps(eltype(x))
Base.one(x::AbstractArray) = one(eltype(x))

"""
    rsimilar(rng, f!, x, n)

Generate random numbers in a container similar to `x`.
The returned variable will have the same shape of `x` if `n == 1` and
will extend an extra dimension with size `n` if `n != 1`.
"""
function _rsimilar(rng::AbstractRNG, f!::Function, x::AbstractArray, dims::Int...)
    u = similar(x, size(x)..., dims...)
    f!(rng, u)
    return u
end

rsimilar(rng, f!, x::AbstractArray, dims::Int...) = _rsimilar(rng, f!, x, dims...)

randsimilar(rng::AbstractRNG, x::AbstractArray, dims::Int...) = rsimilar(rng, rand!, x, dims...)
randnsimilar(rng::AbstractRNG, x::AbstractArray, dims::Int...) = rsimilar(rng, randn!, x, dims...)

### Use below if we don't want to use in-place rand cores

#  function _rsimilar(rng::AbstractRNG, f::Function, x::AbstractArray, dims::Int...)
#      u = f(rng, eltype(x), size(x)..., dims...)
#      return u
#  end

#  rsimilar(rng, f, x::AbstractArray, dims::Int...) = _rsimilar(rng, f, x, dims...)

#  randsimilar(rng::AbstractRNG, x::AbstractArray, dims::Int...) = rsimilar(rng, rand, x, dims...)
#  randnsimilar(rng::AbstractRNG, x::AbstractArray, dims::Int...) = rsimilar(rng, randn, x, dims...)

### Distributions

using Distributions: Distributions, VariateForm, ValueSupport, Discrete, Continuous, Distribution, ContinuousMultivariateDistribution
import Distributions: logpdf, pdf, cdf, invlogcdf, ccdf, rand, mean, std, var, mode, minimum, maximum
import Distributions: Bernoulli, Normal

struct Batch <: VariateForm end
const BatchDistribution{S<:ValueSupport} = Distribution{Batch,S}
const DiscreteBatchDistribution = Distribution{Batch,Discrete}
const ContinuousBatchDistribution = Distribution{Batch,Continuous}
rand(bd::BatchDistribution, dims::Int...; kwargs...) = rand(GLOBAL_RNG, bd, dims...; kwargs...)

include("noise.jl")
export UniformNoise, GaussianNoise
include("normal.jl")
export Normal, NormalStd, BroadcastedNormalStd, NormalVar, BroadcastedNormalVar, NormalLogStd, BroadcastedNormalLogStd, NormalLogVar, BroadcastedNormalLogVar
include("gumbel.jl")
export GumbelSoftmax, GumbelSoftmax2D, GumbelBernoulli, GumbelBernoulliLogit
include("bernoulli.jl")
export BatchBernoulli, BatchBernoulliLogit
export Bernoulli, BernoulliLogit

export logpdf, pdf, cdf, invlogcdf, ccdf, rand, mean, std, var, mode, minimum, maximum
export logpdflogit, logpdfCoV, logrand, logitrand, logvar, kldiv

### X

include("ad.jl")

function __init__()
    Requires.@require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" begin
        import Flux

        Flux.adapt(::Type{Array}, x::UniformNoise{T,<:Any,S}) where {T,S} = UniformNoise(T, :cpu, x.size)
        Flux.adapt(::Type{Array}, x::GaussianNoise{T,<:Any,S}) where {T,S} = GaussianNoise(T, :cpu, x.size)

        @eval Flux.@functor GumbelSoftmax
        @eval Flux.@functor GumbelBernoulli
        @eval Flux.@functor GumbelBernoulliLogit
        @eval Flux.@functor BatchBernoulli
        @eval Flux.@functor BatchBernoulliLogit
    end

    Requires.@require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
        import .CuArrays

        # NOTE: The two functions below achive the following behaviour
        # - Pass `rng` if it's of type `CuArrays.CURAND.RNG`;
        # - Ignore `rng` and use the global of `CuArrays.CURAND` otherwise.
        # The motivation is to avoid scalar operations on GPUs, which is the case when
        # a CPU's RNG is used for inplace random number generation on GPUs.
        rsimilar(rng::CuArrays.CURAND.RNG, f!, x::CuArrays.CuArray, dims::Int...) = _rsimilar(rng, f!, x, dims...)
        rsimilar(::AbstractRNG, f!, x::CuArrays.CuArray, dims::Int...) = _rsimilar(CuArrays.CURAND.generator(), f!, x, dims...)

        ### Use below if we want rsimilar to be reproducible

        # TODO: add a global switch
        #  rsimilar(rng::AbstractRNG, f, x::CuArrays.CuArray, dims::Int...) = _rsimilar(rng, f, x, dims...) |> cu

        CuArrays.cu(x::UniformNoise) = UniformNoise(Float32, :gpu, x.size)
        CuArrays.cu(x::GaussianNoise) = GaussianNoise(Float32, :gpu, x.size)
    end
end

### Test

using Test: @test

"""
    test_stat(fstat, dist, n_samples, atol; samples)

Test the empirical statistic using samples against the true one.
"""
function test_stat(
    fstat::Function,
    dist::BatchDistribution,
    samples::AbstractArray,
    atol::AbstractFloat
)
    dim = length(size(samples))
    stat_est = dropdims(fstat(samples; dims=dim); dims=dim)
    stat_exact = fstat(dist)
    @test stat_est ≈ stat_exact atol=atol
    return samples
end

test_stat(
    fstat::Function,
    dist::BatchDistribution,
    n_samples::Int,
    atol::AbstractFloat
) = test_stat(fstat, dist, rand(dist, n_samples), atol)

end # module


text:
function _init(a, b::Curve{D, T}) where {D, T}
    leaves = PriorityQueue{Interval{T}, T}()
    glb, gub = getbounds(a, b, b.limits)
    soln = b.limits
    leaves[soln] = glb
    return leaves, glb, gub, soln
end

function _init(a::Curve{D, T}, b::Curve{D, T}) where {D, T}
    leaves = PriorityQueue{IntervalBox{2, T}, T}()
    glb, gub = getbounds(a, b, a.limits, b.limits)
    soln = IntervalBox(a.limits, b.limits)
    leaves[soln] = glb
    return leaves, glb, gub, soln
end

function _init(b::Curve, a)
    return _init(a, b)
end

function _loop(a, b, leaves, glb, gub)
    leaf = dequeue!(leaves)

    left, right = bisect(leaf)

    lb, ub = getbounds(a, b, left...)
    lb = lb > glb ? lb : glb
    gub = gub > ub ? ub : gub
    leaves[left] = lb

    lb, ub = getbounds(a, b, right...)
    lb = lb > glb ? lb : glb
    gub = gub > ub ? ub : gub
    leaves[right] = lb

    soln, glb = peek(leaves)
    return glb, gub, soln
end

function _loop(b::Curve, a, leaves, glb, gub)
    return _loop(a, b, leaves, glb, gub)
end

function _loop(a::Curve, b::Curve, leaves, glb, gub)
    leaf = dequeue!(leaves)

    left, right = bisect(leaf)

    lb, ub = getbounds(a, b, left...)
    lb = lb > glb ? lb : glb
    gub = gub > ub ? ub : gub
    leaves[left] = lb

    lb, ub = getbounds(a, b, right...)
    lb = lb > glb ? lb : glb
    gub = gub > ub ? ub : gub
    leaves[right] = lb

    soln, glb = peek(leaves)
    return glb, gub, soln
end

function closest_points(a, b; atol=1e-8)
    leaves, glb, gub, soln = _init(a, b)
    while (gub - glb) > atol
        glb, gub, soln = _loop(a, b, leaves, glb, gub)
    end
    return closest_points(a, b, soln)
end

function closest_points(a, b::Curve{D}, soln) where {D}
    evalb = b(mid(soln))
    return closest_points(a, evalb, @SVector(ones(D)))
end

function closest_points(b::Curve{D}, a, soln) where {D}
    pts = closest_points(a, b, soln)
    reverse(pts)
end

function closest_points(a::Curve{D}, b::Curve{D}, soln) where {D}
    evala = a(mid(soln[1]))
    evalb = b(mid(soln[2]))
    return closest_points(evala, evalb, @SVector(ones(D)))
end

function minimum_distance(a, b; atol=1e-8)
    leaves, glb, gub, soln = _init(a, b)
    while (gub - glb) > atol
        glb, gub, soln = _loop(a, b, leaves, glb, gub)
    end
    return glb
end

function tolerance_verification(a, b, Δ::Real; atol=1e-8)
    leaves, glb, gub, soln = _init(a, b)
    while (gub - Δ) > atol
        glb, gub, soln = _loop(a, b, leaves, glb, gub)
        if glb > Δ
            return true
        end
    end
    return false
end

function collision_detection(a, b; atol=1e-8)
    leaves, glb, gub, soln = _init(a, b)
    while gub > atol
        glb, gub, soln = _loop(a, b, leaves, glb, gub)
        if glb > 0
            return false
        end
    end
    return true
end


text:
using RandomMatrices
using Base.Test

if isdefined(:_HAVE_GSL)
N=5
A=randn(N,N)
B=randn(N,N)
Q=UniformHaar(2, N)

#Test case where there are no symbolic Haar matrices
@test_approx_eq eval(expectation(N, :Q, :(A*B))) A*B
#Test case where there is one pair of symbolic Haar matrices
@test_approx_eq eval(expectedtrace(N, :Q, :(A*Q*B*Q'))) trace(A)*trace(B)/N

println("Case 3")
println("E(A*Q*B*Q'*A*Q*B*Q') = ", eval(expectation(N, :Q, :(A*Q*B*Q'*A*Q*B*Q'))))

for elty in (Float64, Complex128)
	A = Stewart(elty, N)
	@test_approx_eq A'A eye(N)
end

end #_HAVE_GSL


text:
using OrdinaryDiffEq: DiscreteProblem, FunctionMap, solve

change_tspan(prob, tspan::Tuple) =
    DiscreteProblem(prob.f, prob.u0, tspan, prob.p)

change_tspan(prob, tchunk::Integer) = change_tspan(prob, (0, tchunk))

"""Convert `prob.tspan` to `Tuple{Float64, Float64}`."""
float_tspan(prob) = change_tspan(prob, Tuple{Float64, Float64}(prob.tspan))

"""
    solve_discrete(prob::DiscreteProblem; <keyword arguments>)

Solve prob using DifferentialEquations.jl with a bit of hack to make
it work.
See: https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/issues/251
"""
solve_discrete(prob; kwargs...) =
    solve(
        # Convert tspan to Tuple{Float64, Float64}:
        float_tspan(prob),
        # Avoid DifferentialEquations.jl to repeat the initial condition:
        FunctionMap(scale_by_time=false);
        kwargs...)


text:
#=
Author: Stephan Scholz
Year: 2021

This file contains a Simulation of the Heat equation
Boundary conditions: heat transfer and heat radiation (nonlinear)
Dimension: 1D
Spatial Approximation: Finite Differences
Temporal Approximation: Forward-Euler method
=#

using LinearAlgebra, Plots, SparseArrays, DifferentialEquations

# Define the constants for the PDE
const λ = 45.0          # Conductivity
const ρ = 7800.0;       # Density
const cap = 480.0;      # Specific heat capacitivity
const α = λ/(cap * ρ)   # Diffusivity
const L = 0.1           # Length of rod

#Boundary conditions
const h = 10            # Conduction coefficient
const ϵ = 0.6;          # Emissivity
const sb = 5.67*10^(-8) # Stefan-Boltzmann constant
const k = ϵ * sb        # Radiation coefficient

# Discretization
const dx = L/100;                     # Spatial sampling: o--o--o--o--o--o
const N = round(Int, (L / dx + 1));   # Number of elements

const Tf = 300.0    # Final time
const dt = 10^(-2); # Time step width
const t_len = round(Int, (Tf/dt + 1))

#Laplace operator
const M = spdiagm(-1 => ones(N-1), 0 => -2*ones(N), 1 => ones(N-1))
M[1,2] = 2
M[end,end-1] = 2

const Mx = zeros(N)
const C = spzeros(N,2) # Boundary Conditions: Output / Outflow

C[1,1] = 1	  # Left side
C[end,2] = 1  # Right side

const q = α*dt/(dx^2);
if q > 0.5
  error("Numerical stability is NOT guaranteed! STOP program!")
end

# Temperatures
const ϑ₀   = 400.0; # Initial temperature
const ϑamb = 273.0; # Ambient temperature

# Define initial temperature
function init_temperature(ϑoffset, scale, dx, num_points)
	a = scale
	b = ϑoffset

	L = (num_points-1)*dx
	x = 0 : dx : L

	ϑinit = a*sin.(x*2*pi/L) .+ b

	return ϑinit
end

# Define the discretized PDE as an ODE function
function heat_eq(dθ,θ,p,t)

  λ = p[1]
  ρ = p[2]
  cap = p[3]

  mul!(Mx,M,θ[1:end]) # means: Mx = M * Θ[1:end]

  # Now: zero Neumann boundary conditions
  bc = zeros(2)

  # Later: Natural boundary condition: heat conduction + radiation
  bc[1] =  -1*(h*(θ[1] - ϑamb) + k*(θ[1]^4 - ϑamb^4));
  bc[2] =  -1*(h*(θ[end] - ϑamb) + k*(θ[end]^4 - ϑamb^4));

  flux_out =(2*dx/λ)*C*bc                     # Temperature input and output

  dθ .= λ/(ρ*cap) * 1/(dx^2) * (Mx + flux_out)   # Integration of temperatures

end

θinit = init_temperature(ϑ₀, 30.0, dx, N)

tspan = (0.0, Tf) # Time span
param = [λ, ρ, cap] # Parameter that shall be learned

prob = ODEProblem( heat_eq, θinit,tspan, param )
sol = solve(prob,Euler(),dt=dt,progress=true,save_everystep=true,save_start=true) # use 'save_everystep' only for small systems with few timesteps

xspan = 0 : dx : L

# Temperature at the beginning (t=0) vs at the end (t=Tf)
# plot(xspan, [sol[:,1], sol[:,end]], label=["t=0" "t=Tf"], xlabel="Position x", ylabel="Temperature")

# Temperature at left and right side
# plot(sol.t,[sol[1,:], sol[end,:]], label=["Left" "Right"], xlabel="Time t", ylabel="Temperature", legend=:bottomright)

# 2D Heatmap plot: xaxis = time ; yaxis = position
# heatmap(sol.t, xspan, sol[:,:], xlabel="Time t", ylabel="Position x")


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

using ValueShapes
using Test

using Random
using ArraysOfArrays
import TypedTables

@testset "named_tuple_shape" begin
    @testset "functionality" begin
        get_y(x) = x.y

        data = VectorOfSimilarVectors(reshape(collect(1:22), 11, 2))
        ref_table = TypedTables.Table(
            a = [[1 3 5; 2 4 6], [12 14 16; 13 15 17]],
            b = [7, 18],
            c = [4.2, 4.2],
            x = Matrix[[11 21; 12 22], [11 21; 12 22]],
            y = [[8, 9, 10, 11], [19, 20, 21, 22]]
        )

        named_shapes = (
            a = ArrayShape{Real}(2, 3),
            b = ScalarShape{Real}(),
            c = ConstValueShape(4.2),
            x = ConstValueShape([11 21; 12 22]),
            y = ArrayShape{Real}(4)
        )

        shape = @inferred NamedTupleShape(;named_shapes...)
        @test @inferred(NamedTupleShape(named_shapes)) == shape

        @test @inferred(get_y(shape)) === ValueShapes._accessors(shape).y
        @test @inferred(Base.propertynames(shape)) == (:a, :b, :c, :x, :y)
        @test @inferred(totalndof(shape)) == 11

        @test @inferred(shape(data[1])) == ref_table[1]
        @test @inferred(broadcast(shape, data)) == ref_table

        @test @inferred(merge((foo = 42,), shape)) == merge((foo = 42,), named_shapes)
        @test @inferred(NamedTupleShape(;shape...)) == shape

        @test typeof(@inferred(shape(undef))) == NamedTuple{(:a, :b, :c, :x, :y),Tuple{Array{Float64,2},Float64,Float64,Array{Int,2},Array{Float64,1}}}
        @test typeof(@inferred(valshape(shape(undef)))) <: NamedTupleShape
        @test typeof(valshape(shape(undef))(undef)) == NamedTuple{(:a, :b, :c, :x, :y),Tuple{Array{Float64,2},Float64,Float64,Array{Int,2},Array{Float64,1}}}
        @test @inferred(shape(collect(1:11))) == (a = [1 3 5; 2 4 6], b = 7, c = 4.2, x = [11 21; 12 22], y = [8, 9, 10, 11])
        @test_throws ArgumentError shape(collect(1:12))

        @test valshape(shape.(push!(@inferred(VectorOfSimilarVectors(shape)), @inferred(Vector(undef, shape))))[1]) == valshape(shape(undef))
    end

    @testset "examples" begin
        @test begin
            shape = NamedTupleShape(
                a = ArrayShape{Real}(2, 3),
                b = ScalarShape{Real}(),
                c = ConstValueShape(4.2),
                x = ConstValueShape([11 21; 12 22]),
                y = ArrayShape{Real}(4)
            )
            data = VectorOfSimilarVectors{Int}(shape)
            resize!(data, 10)
            rand!(flatview(data), 0:99)
            table = shape.(data)
            fill!(table.b, 42)
            all(x -> x == 42, view(flatview(data), 7, :))
        end
    end
end


text:
mutable struct Shape
    inner::Ptr{primitivShape_t}

    function Shape(inner::Ptr{primitivShape_t})
        @assert inner != C_NULL
        shape = new(inner)
        finalizer(shape, obj -> @call(DeleteShape, (Ptr{primitivShape_t},), obj.inner))
        shape
    end

    function Shape()
        shape_ptr_ptr = Ptr{primitivShape_t}[0]
        @call(CreateShape, (Ptr{Ptr{primitivShape_t}},), shape_ptr_ptr)
        Shape(shape_ptr_ptr[1])
    end

    function Shape(dims::AbstractArray{T}, batch::Integer = 1) where T<:Integer
        shape_ptr_ptr = Ptr{primitivShape_t}[0]
        @call(CreateShapeWithDims, (Ptr{Cuint}, Cint, Cuint, Ptr{Ptr{primitivShape_t}}), convert(Array{UInt32}, dims), size(dims)[1], batch, shape_ptr_ptr)
        Shape(shape_ptr_ptr[1])
    end
end

"""
    depth(shape)

Returns the depth (length of non-1 dimensions) of the shape.

# Arguments
- `shape::Shape` shape object
"""
function depth(shape::Shape)
    retval_ptr = Ref{UInt32}(0)
    @call(GetShapeDepth, (Ptr{primitivShape_t}, Ptr{Cuint}), shape.inner, retval_ptr)
    retval_ptr[]
end

"""
    batch(shape)

Returns the batch size.

# Arguments
- `shape::Shape` shape object
"""
function batch(shape::Shape)
    retval_ptr = Ref{UInt32}(0)
    @call(GetShapeBatchSize, (Ptr{primitivShape_t}, Ptr{Cuint}), shape.inner, retval_ptr)
    retval_ptr[]
end

"""
    volume(shape)

Returns the number of elements in each sample.

# Arguments
- `shape::Shape` shape object
"""
function volume(shape::Shape)
    retval_ptr = Ref{UInt32}(0)
    @call(GetShapeVolume, (Ptr{primitivShape_t}, Ptr{Cuint}), shape.inner, retval_ptr)
    retval_ptr[]
end

"""
    lower_volume(shape)

Returns the number of elements in 1 to specified dim.

# Arguments
- `shape::Shape` shape object
- `dim::Integer` Upper bound of the dimension.
"""
function lower_volume(shape::Shape, dim::Integer)
    retval_ptr = Ref{UInt32}(0)
    @call(GetShapeLowerVolume, (Ptr{primitivShape_t}, Cuint, Ptr{Cuint}), shape.inner, dim, retval_ptr)
    retval_ptr[]
end


text:
using ProgressBars, Plots, StatsBase, Statistics, LaTeXStrings, JLD

function RandomWalker(P,TimeSteps)
    X = 0
    for i in 1:TimeSteps
        X += sample([1,-1], Weights([P,1-P]))
    end
    return X
end

function Linear_fit(X, Y)
    A = [hcat(X) reshape(ones(length(X)), length(X), 1)]
    b = reshape(hcat(Y), length(Y), 1)
    line = (A \ b)
    return line
end

PList = [0.1,0.3,0.5,0.7,0.9]
N = 100000
TimeSteps = 10000
RES = [[] for i = 1:5]

for n in 1:5
    P = PList[n]
    for i in ProgressBar(1:N)
        push!(RES[n],RandomWalker(P,TimeSteps))
    end
    push!(Avg,mean(RES[n]))
end

save("../../Data/Q4/Q4-RW.jld", "Data", RES)

Line = round.(Linear_fit(PList,Avg))

p1 = histogram(RES[1],bins = 50, label = L"P = %$(PList[1]), \ l = \tau = 1", xlabel=L"\langle X_{final} \rangle_{(P)}", ylabel = L"Number", c = :steelblue)
p2 = histogram(RES[2],bins = 50, label = L"P = %$(PList[2]), \ l = \tau = 1", xlabel=L"\langle X_{final} \rangle_{(P)}", ylabel = L"Number", c = :purple)
p3 = histogram(RES[3],bins = 50, label = L"P = %$(PList[3]), \ l = \tau = 1", xlabel=L"\langle X_{final} \rangle_{(P)}", ylabel = L"Number", c = :red)
p4 = histogram(RES[4],bins = 50, label = L"P = %$(PList[4]), \ l = \tau = 1", xlabel=L"\langle X_{final} \rangle_{(P)}", ylabel = L"Number", c = :gold)
p5 = histogram(RES[5],bins = 50, label = L"P = %$(PList[5]), \ l = \tau = 1", xlabel=L"\langle X_{final} \rangle_{(P)}", ylabel = L"Number", c = :green)
p6 = begin
    plot(PList,Avg,label = L"Y = %$(round(Line[1],digits=3)) X %$(round(Line[2],digits=3))", c = :black)
    scatter!(PList,Avg,label = L"Data\ Points",legend = 140,xlabel = L"P",ylabel = L"\langle \overline{X_{final}} \rangle_{(P)}",c = :steelblue)
end

Pmain1 = begin
    plot(plot(plot(p1,p2,layout=(1,2)),plot(p3,p4,layout=(1,2)),plot(p5,p6,layout=(1,2)),layout=(3,1)))
    plot!(plot_title = L"Distribution\ of\ X_{final}\ (N = %$N\ , Time\ Steps=%$TimeSteps)",size = (1000,1000))
end
savefig(Pmain1,"../../Figs/Q4/Q4-tot.pdf")


text:
"""
Package: ForecastPlots

    fplot(x0, x1, px1, t=1:size(x0,1)+size(x1,1), args...;
          plt = palette([:blue,:red],max(2,size(x0,2))), kw...)

Plot time series with its prediction intervals

# Arguments
- `x0`: Vector or Matrix with time series
- `x1`: Vector of Matrix forecast for x0
- `px1`: Prediction intervals for x1 with size(px1) == (size(x1,1),2,N) where `N` is the number of intervals provided with the format [low_x1, high_x1].
- `t`: Values for the temporal axis, it defaults to intengers starting at 1.
- `args...`: `plot` arguments
- `plt`: palette of colors for the time series, it defaults to palette([:blue,:red],max(2,size(x0,2)))
- `kw...`: `plot` keyword arguments

# Returns
Forecast plots

# Examples
```julia-repl
x0 = rand(20);
x1 = rand(10);
px1 = zeros(10,2,2);
px1[:,1,1] .= -0.2; px1[:,2,1] .= 0.2;
px1[:,1,2] .= -0.4; px1[:,2,2] .= 0.4;

fplot(x0,x1,px1)

xb0 = rand(20) .+ 2;
xb1 = rand(10) .+ 2;
px1 = zeros(10,2,2,2)
px1[:,1,1,1] .= -0.2; px1[:,2,1,1] .= 0.2;
px1[:,1,2,1] .= -0.4; px1[:,2,2,1] .= 0.4;
px1[:,1,1,2] .= -0.2; px1[:,2,1,2] .= 0.2;
px1[:,1,2,2] .= -0.4; px1[:,2,2,2] .= 0.4;

xx0 = [x0 xb0]
xx1 = [x1 xb1]
fplot(xx0, xx1, px1)

using Dates
fplot(x0,x1,px1,now()-Day(20+10-1):Day(1):now())
```
"""
function fplot(x0, x1, px1, t=1:size(x0,1)+size(x1,1), args...;
               plt = palette([:blue,:red],max(2,size(x0,2))), kw...)

    sx0 = size(x0)
    sx1 = size(x1)
    spx = size(px1)
    spx = (sx0[1],spx[2:end]...)

    x0 = vcat(x0,repeat([missing],sx1...))
    x1 = vcat(repeat([missing],sx0...),x1)
    px1 = vcat(repeat([missing],spx...), px1)

    plot(t,x0,args...;kw...)

    for i in 1:size(x0,2)
        for j in size(px1,3):-1:1
            plot!(t,x1[:,i]; ribbon = px1[:,:,j,i],
                  linestyle = :dot, label = :none,
                  color = plt[i], fillalpha = 0.2)
        end
        plot!(t,x1[:,i]; linestyle = :dot, label = :none, color = plt[i])
    end
    current()

end


text:
"""
Vshale(args...)

Calculate the Vshale from logs

`` Vshale=\\frac { { GR }_{ log }-{ GR }_{ Sand } }{ { GR }_{ Sh }-{ GR }_{ Sand } } ``

## Input Example

```julia
Vshale(Gr,GrSand,GrShale)
```
"""
function Vshale(Gr,GrSand,GrShale)
    Vsh=(Gr.-GrSand)./(GrShale.-GrSand)
end


text:
using Pkg
Pkg.activate(".")
using Revise, CCPH, Statistics, Distributions, Plots, JLD
import Dates, CSV, Interpolations, AdaptiveMCMC, BlackBoxOptim, Optim

include("./mcmc/metropolis_log_P.jl")
include("./mcmc/analysis.jl")
include("./Weather_Rosinedal_Struct/create_weather_struct_RO.jl")
include("./Script/Fit_Model_2_RO_Data/model_weather_ts_tuning_C_F_Auxiliary.jl")
include("./Transpiration/transpiration.jl")
#include("Script/Fit_Model_2_RO_Data/find_guess.jl")
#include("Script/Fit_Model_2_RO_Data/model_weather_ts_tuning_C_F_GPP_only_Few_para_NoRoot.jl")
#include("Script/Fit_Model_2_RO_Data/model_tuning_C_F_GPP_Ec_New_logP.jl")
#include("Script/Fit_Model_2_RO_Data/model_tuning_C_F_GPP_Ec_New_logP_Opt.jl")
#include("Script/Fit_Model_2_RO_Data/model_tuning_C_F_GPP_Ec_New_logP_Opt_No_Error.jl")
#include("Script/test_model_data.jl")
#include("Script/run_validation.jl")
include("Script/Fit_Model_2_RO_Data/model_tuning_C_F_GPP_Ec_New_logP_Opt_temp.jl")


text:
export multifilter, multifilter!

"""
    multifilter!(x; samples = nothing, experiments = nothing, colnames = nothing)

Filters the `MultiAssayExperiment` `x` in place so that it only contains the specified samples, experiments or column names.
This returns a reference to the modified `x`.

See [`filtersamplemap`](@ref) for the accepted values of `samples`, `experiments` and `colnames`.
The behavior of this function is equivalent to calling [`filtersamplemap!`](@ref) followed by [`dropunused!`](@ref).
The aspects of `x` that are dropped depend on the arguments:

- Unused samples are dropped from the sample data iff `samples != nothing`.
- Unused experiments are dropped from `experiments(x)` iff `experiments != nothing`.
This may result in some experiments with zero columns after filtering,
but is generally more consistent behavior than experiments disappearing without notice.
- Unused columns are dropped from each experiment if `colnames != nothing`, or if `dropcolnames = true` and `samples != nothing`.
The latter condition provides the expected default behavior, where filtering on the samples is expected to propagate to the corresponding columns of each experiment;
this can be disabled by setting `dropcolnames = false`.

# Examples
```jldoctest
julia> using MultiAssayExperiments

julia> x = exampleobject();

julia> multifilter!(x; samples = ["Patient2", "Patient3"], experiments = "foo")
MultiAssayExperiment object
  experiments(1): foo
  sampledata(2): name disease
  metadata(1): version

julia> experiment(x)
100x6 SummarizedExperiments.SummarizedExperiment
  assays(3): foo bar whee
  rownames: Gene1 Gene2 ... Gene99 Gene100
  rowdata(2): name Type
  colnames: foo4 foo5 ... foo8 foo9
  coldata(3): name Treatment Response
  metadata(1): version
```
"""
function multifilter!(x::MultiAssayExperiment; samples = nothing, experiments = nothing, colnames = nothing, dropcolnames = true)
    filtersamplemap!(x; samples = samples, experiments = experiments, colnames = colnames)

    dropsam = (samples != nothing)
    dropexp = (experiments != nothing)
    dropcol = ((dropsam && dropcolnames) || colnames != nothing)

    dropunused!(x; samples = dropsam, experiments = dropexp, colnames = dropcol)
    return x
end

"""
    multifilter!(x; samples = nothing, experiments = nothing, colnames = nothing)

Return a new `MultiAssayExperiment` that has been filtered to only the specified samples, experiments or column names.
This makes a copy of `x` and passes it (and any keyword arguments in `kwargs`) to [`multifilter!`](@ref);
see the latter function for more details.

# Examples
```jldoctest
julia> using MultiAssayExperiments

julia> x = exampleobject();

julia> multifilter(x; samples = ["Patient2", "Patient3"], experiments = "foo")
MultiAssayExperiment object
  experiments(1): foo
  sampledata(2): name disease
  metadata(1): version
```
"""
function multifilter(x::MultiAssayExperiment; kwargs...)
    y = copy(x)
    multifilter!(y; kwargs...)
    return y
end


text:
# http://mathworld.wolfram.com/SphericalSpiral.html

```julia
using CatmullRom, Plots

x(t) = cos(t) / sqrt(1 + t*t)
y(t) = sin(t) / sqrt(1 + t*t)
z(t) = - t/sqrt(1+t*t)

ts = collect(range(-4pi,4pi, length=24));
xs = x.(ts); ys = y.(ts); zs = z.(ts);

czs, cys, cxs = catmullrom(collect(zip(zs,ys,xs)), 36);

plot(zs, ys, xs, linecolor=:darkred, linewidth=2, legend=nothing, size=(600,600))
plot!(czs, cys, cxs, linecolor=:navy, linewidth=2, legend=nothing, size=(600,600))

```


text:
function initialize_tg!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{M, S, A, PSY.TGFixed, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, P <: PSY.PSS}
    tg = PSY.get_prime_mover(dynamic_device)
    τm0 = get_inner_vars(dynamic_device)[τm_var]
    eff = PSY.get_efficiency(tg)
    P_ref = τm0 / eff
    PSY.set_P_ref!(tg, P_ref)
    #Update Control Refs
    PSY.get_ext(dynamic_device)[CONTROL_REFS][P_ref_index] = P_ref
end

function initialize_tg!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{M, S, A, PSY.TGTypeI, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, P <: PSY.PSS}

    #Get mechanical torque to SyncMach
    τm0 = get_inner_vars(dynamic_device)[τm_var]
    #Get Parameters
    tg = PSY.get_prime_mover(dynamic_device)
    R = PSY.get_R(tg)
    Tc = PSY.get_Tc(tg)
    T3 = PSY.get_T3(tg)
    T4 = PSY.get_T4(tg)
    T5 = PSY.get_T5(tg)

    #Get References
    ω_ref = PSY.get_ext(dynamic_device)[CONTROL_REFS][ω_ref_index]
    ω0 = 1.0

    function f!(out, x)
        P_ref = x[1]
        x_g1 = x[2]
        x_g2 = x[3]
        x_g3 = x[4]

        #Compute auxiliary parameters
        P_in = P_ref + (1.0 / R) * (ω_ref - ω0)

        out[1] = P_in - x_g1
        out[2] = (1.0 - T3 / Tc) * x_g1 - x_g2
        out[3] = ((1.0 - T4 / T5) * (x_g2 + (T3 / Tc) * x_g1) - x_g3)
        out[4] = x_g3 + (T4 / T5) * (x_g2 + (T3 / Tc) * x_g1) - τm0
    end
    x0 = [τm0, τm0, (1.0 - T3 / Tc) * τm0, τm0]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Synch. Machine failed")
    else
        sol_x0 = sol.zero
        #Update Control Refs
        PSY.set_P_ref!(tg, sol_x0[1])
        PSY.get_ext(dynamic_device)[CONTROL_REFS][P_ref_index] = sol_x0[1]
        #Update states
        tg_ix = get_local_state_ix(dynamic_device, PSY.TGTypeI)
        tg_states = @view device_states[tg_ix]
        tg_states[1] = sol_x0[2]
        tg_states[2] = sol_x0[3]
        tg_states[3] = sol_x0[4]
    end
end

function initialize_tg!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{M, S, A, PSY.TGTypeII, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, P <: PSY.PSS}

    #Get mechanical torque to SyncMach
    τm0 = get_inner_vars(dynamic_device)[τm_var]
    #Get parameters
    tg = PSY.get_prime_mover(dynamic_device)
    R = PSY.get_R(tg)
    T1 = PSY.get_T1(tg)
    T2 = PSY.get_T2(tg)
    ω_ref = PSY.get_ext(dynamic_device)[CONTROL_REFS][ω_ref_index]
    ω0 = ω_ref

    function f!(out, x)
        P_ref = x[1]
        xg = x[2]

        out[1] = (1.0 / R) * (T1 / T2) * (ω_ref - ω0) + P_ref / 1.0 + xg - τm0
        out[2] = (1.0 / T2) * ((1.0 / R) * (1 - T2 / T1) * (ω_ref - ω0) - xg)
    end
    x0 = [τm0, 0.0]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Synch. Machine failed")
    else
        sol_x0 = sol.zero
        #Update Control Refs
        PSY.set_P_ref!(tg, sol_x0[1])
        PSY.get_ext(dynamic_device)[CONTROL_REFS][P_ref_index] = sol_x0[1]
        #Update states
        tg_ix = get_local_state_ix(dynamic_device, PSY.TGTypeII)
        tg_states = @view device_states[tg_ix]
        tg_states[1] = sol_x0[2]
    end
end


text:
#!/usr/bin/env julia
###### AC-OPF using Nonconvex ######
#
# implementation reference: https://julianonconvex.github.io/Nonconvex.jl/stable/problem/
# currently does not converge due to an upstream issue with the AD backend Zygote: https://github.com/JuliaNonconvex/Nonconvex.jl/issues/130
#

import PowerModels
import Nonconvex
Nonconvex.@load Ipopt

function solve_opf(file_name)
    time_data_start = time()

    data = PowerModels.parse_file(file_name)
    PowerModels.standardize_cost_terms!(data, order=2)
    PowerModels.calc_thermal_limits!(data)
    ref = PowerModels.build_ref(data)[:it][:pm][:nw][0]

    data_load_time = time() - time_data_start

    time_model_start = time()

    model = Nonconvex.DictModel()

    bus_pd = Dict(i => 0.0 for (i,bus) in ref[:bus])
    bus_qd = Dict(i => 0.0 for (i,bus) in ref[:bus])

    bus_gs = Dict(i => 0.0 for (i,bus) in ref[:bus])
    bus_bs = Dict(i => 0.0 for (i,bus) in ref[:bus])

    for (i,bus) in ref[:bus]
        if length(ref[:bus_loads][i]) > 0
            bus_pd[i] = sum(ref[:load][l]["pd"] for l in ref[:bus_loads][i])
            bus_qd[i] = sum(ref[:load][l]["qd"] for l in ref[:bus_loads][i])
        end

        if length(ref[:bus_shunts][i]) > 0
            bus_gs[i] = sum(ref[:shunt][s]["gs"] for s in ref[:bus_shunts][i])
            bus_bs[i] = sum(ref[:shunt][s]["bs"] for s in ref[:bus_shunts][i])
        end
    end

    br_g = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_b = Dict(i => 0.0 for (i,branch) in ref[:branch])

    br_tr = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_ti = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_tm = Dict(i => 0.0 for (i,branch) in ref[:branch])

    br_g_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_b_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_g_to = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_b_to = Dict(i => 0.0 for (i,branch) in ref[:branch])

    for (i,branch) in ref[:branch]
        g, b = PowerModels.calc_branch_y(branch)
        tr, ti = PowerModels.calc_branch_t(branch)

        br_g[i] = g
        br_b[i] = b

        br_tr[i] = tr
        br_ti[i] = ti
        br_tm[i] = tr^2 + ti^2

        br_g_fr[i] = branch["g_fr"]
        br_b_fr[i] = branch["b_fr"]
        br_g_to[i] = branch["g_to"]
        br_b_to[i] = branch["b_to"]
    end

    for (i,bus) in ref[:bus]
        addvar!(model, "va_$(i)", -Inf, Inf) #va
        addvar!(model, "vm_$(i)", bus["vmin"], bus["vmax"], init=1.0) #vm
    end

    for (i,gen) in ref[:gen]
        addvar!(model, "pg_$(i)", gen["pmin"], gen["pmax"]) #pg
        addvar!(model, "qg_$(i)", gen["qmin"], gen["qmax"]) #qg
    end

    for (l,i,j) in ref[:arcs]
        branch = ref[:branch][l]
        addvar!(model, "p_$(l)_$(i)_$(j)", -branch["rate_a"], branch["rate_a"]) #p
        addvar!(model, "q_$(l)_$(i)_$(j)", -branch["rate_a"], branch["rate_a"]) #q
    end

    #total_callback_time = 0.0
    function opf_objective(x::OrderedDict)
        #start = time()
        cost = 0.0
        for (i,gen) in ref[:gen]
            pg = x["pg_$(i)"]
            cost += gen["cost"][1]*pg^2 + gen["cost"][2]*pg + gen["cost"][3]
        end
        #total_callback_time += time() - start
        return cost
    end

    # TBD once AD is resolved
    #add_ineq_constraint!(model, f)
    #add_eq_constraint!(model, f)

    Nonconvex.set_objective!(model, opf_objective)

    model_variables = -1
    model_constraints = -1

    model_build_time = time() - time_model_start

    time_solve_start = time()

    x0 = NonconvexCore.getinit(model)
    #r = optimize(model, IpoptAlg(), x0, options = IpoptOptions(print_level = 0))
    r = Nonconvex.optimize(model, IpoptAlg(), x0)
    cost = 0.0 # TBD
    feasible = false # TBD

    solve_time = time() - time_solve_start
    total_time = time() - time_data_start

    println("")
    println("\033[1mSummary\033[0m")
    println("   case........: $(file_name)")
    println("   variables...: $(model_variables)")
    println("   constraints.: $(model_constraints)")
    println("   feasible....: $(feasible)")
    println("   cost........: $(round(Int, cost))")
    println("   total time..: $(total_time)")
    println("     data time.: $(data_load_time)")
    println("     build time: $(model_build_time)")
    println("     solve time: $(solve_time)")
    # println("      callbacks: $(total_callback_time)")
    println("")

    return Dict(
        "case" => file_name,
        "variables" => model_variables,
        "constraints" => model_constraints,
        "feasible" => feasible,
        "cost" => cost,
        "time_total" => total_time,
        "time_data" => data_load_time,
        "time_build" => model_build_time,
        "time_solve" => solve_time,
        #"time_callbacks" => TBD,
    )
end

if isinteractive() == false
    solve_opf("$(@__DIR__)/data/pglib_opf_case5_pjm.m")
end

max_stars_repo_name: lanl-ansi/rosetta-opf

max_issues_repo_name: lanl-ansi/rosetta-opf

max_forks_repo_name: lanl-ansi/rosetta-opf


text:
using DDD
using Test, StaticArrays, SparseArrays, LinearAlgebra
cd(@__DIR__)

@testset "nodeFE type" begin
    @test nodeTypeFE(0) == 0
    @test 1 == nodeTypeFE(1)
    @test nodeTypeFE(3) == 3.0
    @test zero(nodeTypeFE) == 0
    @test convert(nodeTypeFE, 3.0) == nodeTypeFE(3)
end

@testset "Shape functions" begin
    points = Float64[
        0 0 0
        1 0 0
        1 1 0
        0 1 0
        0 0 1
        1 0 1
        1 1 1
        0 1 1
        0.5 0 0
        0.5 0.5 0
        0 0.5 0
        0 0 0.5
        0.5 0 0.5
        0.5 0.5 0.5
        0 0.5 0.5
    ]
    Nall = shapeFunction(LinearQuadrangle3D(), points[:, 1], points[:, 2], points[:, 3])
    @test all(isapprox.(sum.(Nall), 1))

    dNdSall =
        shapeFunctionDeriv(LinearQuadrangle3D(), points[:, 1], points[:, 2], points[:, 3])
    checkSum = sum.(dNdSall, dims = 2)
    for i in 1:size(points, 1)
        @test vec(checkSum[i]) == SVector(0.0, 0.0, 0.0)

        N = shapeFunction(LinearQuadrangle3D(), points[i, 1], points[i, 2], points[i, 3])
        @test isapprox(Nall[i], N)

        dNdS = shapeFunctionDeriv(
            LinearQuadrangle3D(),
            points[i, 1],
            points[i, 2],
            points[i, 3],
        )

        @test isapprox(dNdSall[i], dNdS)
    end
end

@testset "Checking arbitrary points of the mesh and connectivity" begin
    DictMaterialParameters = loadJSON("./testData/sampleMaterialParameters.json")
    matParams = loadMaterialParameters(DictMaterialParameters)
    femParams = FEMParameters(;
        type = DispatchRegularCuboidMesh(),
        order = LinearElement(),
        model = CantileverLoad(),
        dx = Float64(1009),
        dy = Float64(1013),
        dz = Float64(1019),
        mx = 11,
        my = 13,
        mz = 17,
    )
    regularCuboidMesh = buildMesh(matParams, femParams)

    @test typeof(regularCuboidMesh) <: AbstractRegularCuboidMesh &&
          typeof(femParams.type) <: AbstractRegularCuboidMesh
    @test regularCuboidMesh.order == femParams.order

    testVertices = [
        0 0 0
        1009 0 0
        1009 1013 0
        0 1013 0
        0 0 1019
        1009 0 1019
        1009 1013 1019
        0 1013 1019
    ]
    testC = [
        3.272727272727273 1.272727272727273 1.272727272727273 0 0 0
        1.272727272727273 3.272727272727273 1.272727272727273 0 0 0
        1.272727272727273 1.272727272727273 3.272727272727273 0 0 0
        0 0 0 1.000000000000000 0 0
        0 0 0 0 1.000000000000000 0
        0 0 0 0 0 1.000000000000000
    ]
    vertices =
        reshape(collect(Iterators.flatten(regularCuboidMesh.vertices.vertices)), 3, 8)
    faces = regularCuboidMesh.faces
    normals = regularCuboidMesh.faceNorm

    @test isapprox(vertices', testVertices)

    faceCoord = vertices[:, faces]
    p = faceCoord[:, 2, :] - faceCoord[:, 1, :]
    q = faceCoord[:, 4, :] - faceCoord[:, 1, :]
    n = reshape(
        collect(Iterators.flatten([normalize(p[:, i] × q[:, i]) for i in 1:size(p, 2)])),
        3,
        6,
    )
    @test isapprox(n, normals)

    mx = regularCuboidMesh.mx
    my = regularCuboidMesh.my
    mz = regularCuboidMesh.mz
    faceNorm = regularCuboidMesh.faceNorm
    lbl = regularCuboidMesh.surfElemNode
    coord = regularCuboidMesh.coord
    for (i, val) in enumerate([
        1,
        1 + mx * mz,
        1 + mx * mz + my * mz,
        1 + 2 * mx * mz + my * mz,
        1 + 2 * mx * mz + 2 * my * mz,
        1 + 2 * mx * mz + 2 * my * mz + mx * my,
    ])
        s = coord[:, lbl[:, val]]
        p = SVector{3, eltype(s)}(s[:, 2] - s[:, 1])
        q = SVector{3, eltype(s)}(s[:, 4] - s[:, 1])
        @test normalize(p × q) ≈ faceNorm[:, i]
    end

    @test isapprox(regularCuboidMesh.C, testC)

    @test regularCuboidMesh.dx == femParams.dx
    @test regularCuboidMesh.dy == femParams.dy
    @test regularCuboidMesh.dz == femParams.dz
    @test regularCuboidMesh.mx == femParams.mx
    @test regularCuboidMesh.my == femParams.my
    @test regularCuboidMesh.mz == femParams.mz
    @test regularCuboidMesh.numElem ==
          regularCuboidMesh.mx * regularCuboidMesh.my * regularCuboidMesh.mz
    @test regularCuboidMesh.numNode ==
          (regularCuboidMesh.mx + 1) *
          (regularCuboidMesh.my + 1) *
          (regularCuboidMesh.mz + 1)

    @test regularCuboidMesh.dx / regularCuboidMesh.mx == regularCuboidMesh.w
    @test regularCuboidMesh.dy / regularCuboidMesh.my == regularCuboidMesh.h
    @test regularCuboidMesh.dz / regularCuboidMesh.mz == regularCuboidMesh.d

    idx = [
        9 3
        9 18
        6 5
        5 20
        10 9
    ]

    testCoord1 = 1.0e+02 * [
        7.338181818181819 0 0
        1.834545454545455 0 0
    ]
    testCoord2 = [733.8181818181819 0.0 0.0; 458.6363636363637 77.92307692307692 0.0]
    testCoord3 = 1.0e+02 * [
        4.586363636363637 0 0
        3.669090909090909 0 0
    ]
    testCoord4 = [366.90909090909093 0.0 0.0; 642.0909090909091 77.92307692307692 0.0]
    testCoord5 = 1.0e+02 * [
        8.255454545454546 0 0
        7.338181818181819 0 0
    ]
    coord = regularCuboidMesh.coord
    @test isapprox(coord[:, idx[1, :]]', testCoord1)
    @test isapprox(coord[:, idx[2, :]]', testCoord2)
    @test isapprox(coord[:, idx[3, :]]', testCoord3)
    @test isapprox(coord[:, idx[4, :]]', testCoord4)
    @test isapprox(coord[:, idx[5, :]]', testCoord5)

    testCon = [
        2508 2519 2352 2339 2507
        1518 1529 1362 1349 1517
        2750 2761 2594 2581 2749
        2750 2761 2594 2581 2749
        2703 2714 2547 2534 2702
    ]
    connectivity = regularCuboidMesh.connectivity

    idxCon = [6, 8, 3, 1, 5]
    idxNode = [2002, 1149, 2201, 2201, 2158]
    @test connectivity[idxCon, idxNode]' == testCon

    KTest = [
        -27.845980957491278
        -45.099576129125296
        -64.81800121385578
        -7.379079254079251
        3.5419580419580443
        29.51631701631701
        5.676247771836005
        -9.901489001393466
        -22.704991087344027
        7.379079254079253
    ]
    K = regularCuboidMesh.K
    droptol!(K, 1e-14)

    idxK = [
        CartesianIndex(2194, 2695)
        CartesianIndex(8235, 7734)
        CartesianIndex(5636, 5672)
        CartesianIndex(7785, 8326)
        CartesianIndex(7018, 7524)
        CartesianIndex(3367, 2868)
        CartesianIndex(6289, 5753)
        CartesianIndex(8368, 7903)
        CartesianIndex(724, 764)
        CartesianIndex(5121, 4654)
    ]
    @test isapprox(K[idxK], KTest)
end

@testset "Boundary conditions" begin
    DictMaterialParameters = loadJSON("./testData/sampleMaterialParameters.json")
    matParams = loadMaterialParameters(DictMaterialParameters)
    femParams = FEMParameters(;
        type = DispatchRegularCuboidMesh(),
        order = LinearElement(),
        model = CantileverLoad(),
        dx = Float64(1009),
        dy = Float64(1013),
        dz = Float64(1019),
        mx = 11,
        my = 13,
        mz = 17,
    )

    regularCuboidMesh = buildMesh(matParams, femParams)
    dx = regularCuboidMesh.dx
    dz = regularCuboidMesh.dz
    faces = regularCuboidMesh.faces
    coord = regularCuboidMesh.coord
    connectivity = regularCuboidMesh.connectivity
    surfNode = regularCuboidMesh.surfNode
    coord = regularCuboidMesh.coord

    cantileverBC, forceDisplacement = Boundaries(femParams, regularCuboidMesh)
    uGamma = cantileverBC.uGamma.node
    mGamma = cantileverBC.mGamma.node
    left = findall(x -> x == 0, coord[1, :])

    loadEdge1 = findall(x -> x ≈ dx, coord[1, :])
    loadEdge2 = findall(x -> x ≈ dz, coord[3, :])
    loadEdge = intersect(loadEdge1, loadEdge2)
    @test norm(coord[:, left]) ≈ norm(coord[:, uGamma])
    @test norm(coord[:, loadEdge]) ≈ norm(coord[:, mGamma])
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
# ---
# title: 964. Least Operators to Express Number
# id: problem964
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, Dynamic Programming
# link: <https://leetcode.com/problems/least-operators-to-express-number/description/>
# hidden: true
# ---
#
# Given a single positive integer `x`, we will write an expression of the form
# `x (op1) x (op2) x (op3) x ...` where each operator `op1`, `op2`, etc. is
# either addition, subtraction, multiplication, or division (`+`, `-`, `*`, or
# `/)`.  For example, with `x = 3`, we might write `3 * 3 / 3 + 3 - 3` which is
# a value of 3.
#
# When writing such an expression, we adhere to the following conventions:
#
#   1. The division operator (`/`) returns rational numbers.
#   2. There are no parentheses placed anywhere.
#   3. We use the usual order of operations: multiplication and division happens before addition and subtraction.
#   4. It's not allowed to use the unary negation operator (`-`).  For example, "`x - x`" is a valid expression as it only uses subtraction, but "`-x + x`" is not because it uses negation.
#
# We would like to write an expression with the least number of operators such
# that the expression equals the given `target`.  Return the least number of
# operators used.
#
#
#
# **Example 1:**
#
#
#
#     Input: x = 3, target = 19
#     Output: 5
#     Explanation: 3 * 3 + 3 * 3 + 3 / 3.  The expression contains 5 operations.
#
#
# **Example 2:**
#
#
#
#     Input: x = 5, target = 501
#     Output: 8
#     Explanation: 5 * 5 * 5 * 5 - 5 * 5 * 5 + 5 / 5.  The expression contains 8 operations.
#
#
# **Example 3:**
#
#
#
#     Input: x = 100, target = 100000000
#     Output: 3
#     Explanation: 100 * 100 * 100 * 100.  The expression contains 3 operations.
#
#
#
# **Note:**
#
#   * `2 <= x <= 100`
#   * `1 <= target <= 2 * 10^8`
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    rewind(angles, period=2π, ref_angle=0)

Given angles or other cyclic coordinates with the specified period, such that
the angle + period has the same meaning as the original angle, this function adds or subtracts
multiples of the period such that the result has the same meaning, but now all angles lie in
an interval of length the specified period, centered on the reference angle `ref_angle`.
"""
function rewind(angles; period=2π, ref_angle=0)
    half_period = period / 2
    return ref_angle .+ mod.(angles .- ref_angle .+ half_period, period) .- half_period
end

function rewind!(angles; period=2π, ref_angle=0)
    half_period = period / 2
    angles .= ref_angle .+ mod.(angles .- ref_angle .+ half_period, period) .- half_period
    return angles
end

function unwind(angles; dims=nothing, period=2π, ref_angle=0)
    wound_angles = rewind(angles; period=period, ref_angle=ref_angle)
    return unwrap(wound_angles .- ref_angle; dims=dims, range=period) .+ ref_angle
end

function unwind!(angles; dims=nothing, period=2π, ref_angle=0)
    rewind!(angles; period=period, ref_angle=ref_angle)
    angles .-= ref_angle
    unwrap!(angles; dims=dims, range=period)
    angles .+= ref_angle
    return angles
end

"""
    pix2sky(m::Enmap, pixcoords)

Convert 1-indexed pixels to sky coordinates. The output sky coordinates are determined by WCS,
but usually are in units of degrees.

# Arguments:
- `m::Enmap`: the map that provides a coordinate system
- `pixcoords`: pixcoords should be a 2-d array where "pixcoords[:, i]" is the i-th set of coordinates,
    or a 1-d array representing a single set of coordinates.

# Returns:
- `Array`: same shape as pixcoords

# Examples
```julia-repl
julia> shape, wcs = fullsky_geometry(deg2rad(1))
julia> pix2sky(shape, wcs [1.0, 1.0])
2-element Vector{Float64}:
 180.0
 -90.0
```
"""
function pix2sky end

pix2sky(m::Enmap, p1; safe=true) = pix2sky(size(m), getwcs(m), p1; safe=safe)
pix2sky(m::Enmap, p1, p2; safe=true) = pix2sky(size(m), getwcs(m), p1, p2; safe=safe)
sky2pix(m::Enmap, p1; safe=true) = sky2pix(size(m), getwcs(m), p1; safe=safe)
sky2pix(m::Enmap, p1, p2; safe=true) = sky2pix(size(m), getwcs(m), p1, p2; safe=safe)

pix2sky!(m::Enmap, p1, p2; safe=true) = pix2sky!(size(m), getwcs(m), p1, p2; safe=safe)
sky2pix!(m::Enmap, p1, p2; safe=true) = sky2pix!(size(m), getwcs(m), p1, p2; safe=safe)

## The default fallback (no projection specified) is to call WCS, which calls the C library.
## If you wanted to replicate Pixell behavior, add 1 to x and y of pixcoords.
## We implement custom routines for CAR (Clenshaw-Curtis variant) instead of using these.
function pix2sky(shape, wcs::WCSTransform, pixcoords; safe=true)
    angle_unit = getunit(wcs)
    skycoords = pix_to_world(wcs, pixcoords) .* angle_unit
    if safe
        return unwind(skycoords; dims=2)
    end
    return skycoords
end
function pix2sky!(shape, wcs::WCSTransform, pixcoords, skycoords; safe=true)
    angle_unit = getunit(wcs)
    pix_to_world!(wcs, pixcoords, skycoords)
    skycoords .*= angle_unit
    if safe
        unwind!(skycoords; dims=2)
    end
    return skycoords
end
function sky2pix(shape, wcs::WCSTransform, skycoords; safe=true)
    angle_unit = getunit(wcs)
    inverse_angle_unit = 1 / angle_unit
    pixcoords = world_to_pix(wcs, skycoords .* inverse_angle_unit)
    if safe
        center_pix = shape[begin:begin+1] ./ 2 .+ 1
        pix_periods = abs.(2π ./ (getcdelt(wcs) .* angle_unit))
        rewind!(view(pixcoords, 1, :); period=pix_periods[1], ref_angle=center_pix[1])
        rewind!(view(pixcoords, 2, :); period=pix_periods[2], ref_angle=center_pix[2])
    end
    return pixcoords
end
function sky2pix!(shape, wcs::WCSTransform, skycoords, pixcoords; safe=true)
    angle_unit = getunit(wcs)
    inverse_angle_unit = 1 / angle_unit
    world_to_pix!(wcs, skycoords .* inverse_angle_unit, pixcoords)
    if safe
        center_pix = shape[begin:begin+1] ./ 2 .+ 1
        pix_periods = abs.(2π ./ (getcdelt(wcs) .* angle_unit))
        rewind!(view(pixcoords, 1, :); period=pix_periods[1], ref_angle=center_pix[1])
        rewind!(view(pixcoords, 2, :); period=pix_periods[2], ref_angle=center_pix[2])
    end
    return pixcoords
end
pix2sky(shape, wcs::WCSTransform,ra::Number, dec::Number; safe=true) =
    pix2sky(shape, wcs, Float64[ra, dec]; safe=safe)
sky2pix(shape, wcs::WCSTransform,ra::Number, dec::Number; safe=true) =
    sky2pix(shape, wcs, Float64[ra, dec]; safe=safe)

function pix2sky(shape, wcs::WCSTransform, ra::AbstractVector, dec::AbstractVector; safe=true)
    result = pix2sky(shape, wcs, convert(Array{Float64,2}, [ra'; dec']); safe=safe)
    return result[1,:], result[2,:]
end
function sky2pix(shape, wcs::WCSTransform, ra::AbstractVector, dec::AbstractVector; safe=true)
    result = sky2pix(shape, wcs, convert(Array{Float64,2}, [ra'; dec']); safe=safe)
    return result[1,:], result[2,:]
end
"""
    pix2sky!(m::Enmap, pixcoords, skycoords)

Convert 1-indexed pixels to sky coordinates, in-place. The output sky coordinates are
determined by WCS, but usually are in units of degrees.

# Arguments:
- `m::Enmap`: the map that provides a coordinate system
- `pixcoords`: pixel coordinates should be a 2-d array where "pixcoords[:, i]" is the i-th
    set of coordinates, or a 1-d array representing a single set of coordinates.
- `skycoords`: output array for sky coordinates, must be same same as pixcoords

# Returns:
- `Array`: skycoords

# Examples
```julia-repl
julia> shape, wcs = fullsky_geometry(deg2rad(1))
       pixcoords =  100 .* rand(2,4096 * 2)
       skycoords =  similar(pixcoords)

julia> pix2sky!(shape, wcs, pixcoords, skycoords)
```
"""
function pix2sky!(shape, wcs::CarClenshawCurtis, pixcoords::AbstractArray{TP,2},
                  skycoords::AbstractArray{TS,2}; safe=true) where {TP,TS}

    angle_unit = getunit(wcs)
    α₀, δ₀ = getcrval(wcs) .* angle_unit
    Δα, Δδ = getcdelt(wcs) .* angle_unit
    iα₀, iδ₀ = getcrpix(wcs)

    # compute RA (α) and DEC (δ)
    @inbounds for ipix ∈ axes(pixcoords, 2)
        iα = pixcoords[1, ipix]
        iδ = pixcoords[2, ipix]
        α = α₀ + (iα - iα₀) * Δα
        δ = δ₀ + (iδ - iδ₀) * Δδ
        skycoords[1, ipix] = α
        skycoords[2, ipix] = δ
    end

    if safe
        unwind!(skycoords; dims=2)
    end

    return skycoords
end

# the not-in-place version just creates an output array and calls the in-place one above
function pix2sky(shape, wcs::CarClenshawCurtis,
                 pixcoords::AbstractArray{TP,2}; safe=true) where TP
    skycoords = similar(pixcoords)
    return pix2sky!(shape, wcs, pixcoords, skycoords; safe=safe)
end

"""
    pix2sky(m::Enmap{T,N,AA,<:CarClenshawCurtis}, ra_pixel, dec_pixel)

Compute the sky position of a single position on the sky.

Only implemented for CAR (Clenshaw-Curtis variant) projections, so
the input map is of type `Enmap{T,N,AA,<:CarClenshawCurtis}`.
This takes pixel indices for RA and DEC, and returns a tuple containing
the corresponding RA and DEC.

# Examples
```julia-repl
julia> shape, wcs = fullsky_geometry(deg2rad(1))
julia> pix2sky(shape, wcs, 30.0, 80.0)
(151.0, -11.0)
```
"""
function pix2sky(shape, wcs::CarClenshawCurtis, ra_pixel, dec_pixel; safe=true)
    angle_unit = getunit(wcs)
    α₀, δ₀ = getcrval(wcs) .* angle_unit
    Δα, Δδ = getcdelt(wcs) .* angle_unit
    iα₀, iδ₀ = getcrpix(wcs)
    α = α₀ .+ (ra_pixel .- iα₀) .* Δα
    δ = δ₀ .+ (dec_pixel .- iδ₀) .* Δδ
    if safe
        return rewind(α), rewind(δ)
    end
    return α, δ
end

# when passing a length-2 vector [ra, dec], return a vector. wraps the pix2sky(shape, wcs, ra_pix, dec_pix)
function pix2sky(shape, wcs::CarClenshawCurtis, pixcoords::AbstractVector; safe=true)
    @assert length(pixcoords) == 2
    skycoords = collect(pix2sky(shape, wcs, first(pixcoords), last(pixcoords)))
    if safe
        unwind!(skycoords; dims=2)
    end
    return skycoords
end

"""
    sky2pix(m::Enmap, skycoords)

Convert sky coordinates to 1-indexed pixels. The input sky coordinates are determined by WCS,
but usually are in units of degrees.

# Arguments:
- `m::Enmap`: the map to obtain the coordinates from
- `skycoords`: skycoords should be a 2-d array where "skycoords[:, i]" is the i-th set of coordinates,
    or a 1-d array representing a single set of coordinates.

# Returns:
- `Array`: same shape as skycoords

# Examples
```julia-repl
julia> shape, wcs = fullsky_geometry(deg2rad(1))
julia> sky2pix(shape, wcs, [30.0, 50.0])
2-element Vector{Float64}:
 151.0
 141.0
```
"""
function sky2pix end

"""
    sky2pix!(m::Enmap, skycoords, pixcoords)

Convert sky coordinates to 1-indexed pixels, in-place. The input sky coordinates are
determined by WCS, but usually are in units of degrees.

# Arguments:
- `m::Enmap`: the map that provides a coordinate system
- `skycoords`: sky coordinates should be a 2-d array where "skycoords[:, i]" is the i-th
    set of coordinates, or a 1-d array representing a single set of coordinates.
- `pixcoords`: output array for pixel coordinates, must be same same as pixcoords

# Returns:
- `Array`: pixcoords

# Examples
```julia-repl
julia> shape, wcs = fullsky_geometry(deg2rad(1))
       skycoords =  similar(pixcoords)
       pixcoords =  100 .* rand(2,4096 * 2)
julia> sky2pix!(shape, wcs, skycoords, pixcoords)
```
"""
function sky2pix!(shape, wcs::CarClenshawCurtis, skycoords::AbstractArray{TS,2},
                  pixcoords::AbstractArray{TP,2}; safe=true) where {TS,TP}

    # retrieve WCS info
    angle_unit = getunit(wcs)
    α₀, δ₀ = getcrval(wcs) .* angle_unit
    Δα, Δδ = getcdelt(wcs) .* angle_unit
    iα₀, iδ₀ = getcrpix(wcs)
    Δα⁻¹, Δδ⁻¹ = 1 / Δα, 1 / Δδ

    # compute RA (α) index and DEC (δ) index
    @inbounds for ipix ∈ axes(pixcoords, 2)
        α = skycoords[1, ipix]
        δ = skycoords[2, ipix]
        iα = iα₀ + (α - α₀) * Δα⁻¹
        iδ = iδ₀ + (δ - δ₀) * Δδ⁻¹
        pixcoords[1, ipix] = iα
        pixcoords[2, ipix] = iδ
    end

    if safe
        center_pix = shape[begin:begin+1] ./ 2 .+ 1
        pix_periods = abs.(2π ./ (Δα, Δδ))
        rewind!(view(pixcoords, 1, :); period=pix_periods[1], ref_angle=center_pix[1])
        rewind!(view(pixcoords, 2, :); period=pix_periods[2], ref_angle=center_pix[2])
    end

    return pixcoords
end

# the not-in-place version just creates an output array and calls the in-place one above
function sky2pix(shape, wcs::CarClenshawCurtis,
                 skycoords::AbstractArray{TS,2}; safe=true) where {TS}
    pixcoords = similar(skycoords)
    return sky2pix!(shape, wcs, skycoords, pixcoords; safe=safe)
end

"""
    sky2pix(m::Enmap{T,N,AA,<:CarClenshawCurtis}, ra, dec)

Compute 1-indexed pixels into sky coordinates.

Only implemented for CAR (Clenshaw-Curtis variant) projections. Takes
RA and DEC and returns a tuple containing the corresponding pixel
indices. If vectors of RA and DEC are given, then vectors of
pixel indices will be returned.

# Examples
```julia-repl
julia> shape, wcs = fullsky_geometry(deg2rad(1))
julia> sky2pix(shape, wcs, 30.0, 80.0)
(151.0, 171.0)
```
"""
function sky2pix(shape, wcs::CarClenshawCurtis, ra::Number, dec::Number; safe=true)
    angle_unit = getunit(wcs)
    α₀, δ₀ = getcrval(wcs) .* angle_unit
    Δα, Δδ = getcdelt(wcs) .* angle_unit
    iα₀, iδ₀ = getcrpix(wcs)
    pix_ra = iα₀ + (ra - α₀) / Δα
    pix_dec = iδ₀ + (dec - δ₀) / Δδ

    if safe
        center_pix = shape[begin:begin+1] ./ 2 .+ 1
        pix_ra = rewind(pix_ra; period=abs(2π / Δα), ref_angle=center_pix[1])
        pix_dec = rewind(pix_dec; period=abs(2π / Δδ), ref_angle=center_pix[2])
    end
    return pix_ra, pix_dec
end
function sky2pix(shape, wcs::CarClenshawCurtis,
                 ra::AV, dec::AV; safe=true) where {AV<:AbstractVector}
    angle_unit = getunit(wcs)
    α₀, δ₀ = getcrval(wcs) .* angle_unit
    Δα, Δδ = getcdelt(wcs) .* angle_unit
    iα₀, iδ₀ = getcrpix(wcs)
    Δα⁻¹, Δδ⁻¹ = 1 / Δα, 1 / Δδ
    pix_ra = iα₀ .+ (ra .- α₀) .* Δα⁻¹
    pix_dec = iδ₀ .+ (dec .- δ₀) .* Δδ⁻¹

    if safe
        center_pix = shape[begin:begin+1] ./ 2 .+ 1
        rewind!(pix_ra; period=abs(2π * Δα⁻¹), ref_angle=center_pix[1])
        rewind!(pix_dec; period=abs(2π * Δδ⁻¹), ref_angle=center_pix[2])
    end

    return pix_ra, pix_dec
end

# when passing a vector [ra, dec], return a vector. wraps the sky2pix(shape, wcs, ra, dec).
function sky2pix(shape, wcs::CarClenshawCurtis, skycoords::AbstractVector; safe=true)
    @assert length(skycoords) == 2
    return collect(sky2pix(shape, wcs::CarClenshawCurtis,
        first(skycoords), last(skycoords); safe=safe))
end

"""Check if a WCS is a cylindrical pixelization"""
function iscyl(wcs::WCSTransform)
    # right now we only support CAR
    ctype = wcs.ctype
    if ctype[1] == "RA---CAR" && ctype[2] == "DEC--CAR"
        return true
    end
    return false
end
iscyl(wcs::CarClenshawCurtis) = true

"""Area of a patch given shape and WCS, in steradians."""
function skyarea(shape, wcs::WCSTransform)
    if iscyl(wcs)
        return skyarea_cyl(shape, wcs)
    end
    error("Area for this WCS is not implemented.")
end
skyarea(shape, wcs::CarClenshawCurtis) = skyarea_cyl(shape, wcs)

# generic cylindrical area calculation
function skyarea_cyl(shape, wcs::AbstractWCSTransform)
    N_α, N_δ = shape
    α, δ = pix2sky(shape, wcs, SVector(0., 0.), SVector(0.5, N_δ+0.5); safe=false)
    δ₁, δ₂ = sort(δ)
    δ₁, δ₂ = max(-π/2, δ₁), min(π/2, δ₂)
    Δα, Δδ = getcdelt(wcs) .* getunit(wcs)
    return (sin(δ₂) - sin(δ₁)) * abs(Δα) * N_α
end

"""Generate a similar Enmap whose pixel values are the areas of the pixels in steradians."""
function pixareamap(m::Enmap{T,N,AA,W}) where {T,N,AA,W<:CarClenshawCurtis}
    pixareas = similar(m)
    pixareamap!(pixareas)
end

function pixareamap(shape, wcs::CarClenshawCurtis)
    pixareas = Enmap(Array{Float64}(undef, shape), wcs)
    pixareamap!(pixareas)
end

"""In-place write to pixels the areas of those pixels in steradians."""
function pixareamap!(pixareas::Enmap)
    shape = size(pixareas)
    wcs = getwcs(pixareas)
    ncols, nrows = shape
    Δα, Δδ = abs.(getcdelt(wcs) .* getunit(wcs))

    for i in 1:nrows
        α, δ = pix2sky(shape, wcs, SVector(1., 1.), SVector(i-0.5, i+0.5); safe=false)
        δ₁, δ₂ = sort(δ)
        δ₁, δ₂ = max(-π/2, δ₁), min(π/2, δ₂)
        pixareas[:,i] .= (sin(δ₂) - sin(δ₁)) * Δα
    end

    return pixareas
end

# this set of slices is necessary because colons are somehow not expanded
slice_geometry(shape::Tuple, wcs, sel_x::Colon, s...) = slice_geometry(shape, wcs, 1:shape[1], s...)
slice_geometry(shape::Tuple, wcs, sel_x::Colon, sel_y::Colon, s...) = slice_geometry(shape, wcs, 1:shape[1], 1:shape[2], s...)
slice_geometry(shape::Tuple, wcs, sel_x, sel_y::Colon, s...) = slice_geometry(shape, wcs, sel_x, 1:shape[2], s...)

# this set of slices is neccesary because slices that remove one of the first two dimensions are not allowed, because
# our WCS structure ALWAYS has two axes. However, it can be desirable to write code like "enmap[1,:] .= 0.0" instead
# even if the object on the left is no longer an Enmap
slice_geometry(shape::Tuple, wcs, sel_x::Integer, sel_y::UnitRange, s...) = slice_geometry(shape, wcs, sel_x:sel_x, sel_y, s...)
slice_geometry(shape::Tuple, wcs, sel_x::Integer, sel_y::Integer, s...) = slice_geometry(shape, wcs, sel_x:sel_x, sel_y:sel_y, s...)
slice_geometry(shape::Tuple, wcs, sel_x::UnitRange, sel_y::Integer, s...) = slice_geometry(shape, wcs, sel_x, sel_y:sel_y, s...)

# does not implement "nowrap" like pixell
function slice_geometry(shape_all::Tuple, wcs, sel_x, sel_y, sel_others...)
    other_dims = shape_all[3:end]
    sel = (sel_x, sel_y)

    starts = map(s -> (step(s) > 0) ? first(s) - 1 : first(s), sel)  # handle backward ranges too
    steps = step.(sel)
    sel_sizes = last.(sel) .- first.(sel) .+ steps
    crpix′ = (getcrpix(wcs) .- (starts .+ 0.5)) ./ steps .+ 0.5
    cdelt′ = getcdelt(wcs) .* steps
    shape = sel_sizes .÷ steps
    wcs′ = sliced_wcs(wcs, cdelt′, crpix′)

    return (shape..., other_dims...), wcs′
end
slice_geometry(m::Enmap, sel_all::Vararg) = slice_geometry(size(m), getwcs(m), sel_all...)

function sliced_wcs(wcs::WCSTransform, cdelt′, crpix′)
    new_wcs = deepcopy(wcs)
    new_wcs.cdelt = collect(cdelt′)
    new_wcs.crpix = collect(crpix′)
    return new_wcs
end

function sliced_wcs(wcs::CarClenshawCurtis{T}, cdelt′, crpix′) where T
    new_wcs = CarClenshawCurtis{T}(cdelt′, crpix′, wcs.crval, wcs.unit)
    return new_wcs
end


text:
using PatModules
@mainmodule subpackage ("one", "two") begin

    say_hi = one.say_hi
    say_bye = two.say_bye
end


text:
const IdProp = Tuple{Symbol, Symbol}

struct CallbackId
    state ::Vector{IdProp}
    input ::Vector{IdProp}
    output ::Vector{IdProp}
end

CallbackId(;input ::Union{Vector{IdProp}, IdProp},
            output ::Union{Vector{IdProp}, IdProp},
            state ::Union{Vector{IdProp}, IdProp} = Vector{IdProp}()
            ) = CallbackId(state, input, output)

Base.convert(::Type{Vector{IdProp}}, v::IdProp) = [v]

struct Callback
    func ::Function
    id ::CallbackId
    pass_changed_props ::Bool
end

struct PreventUpdate <: Exception

end

struct NoUpdate
end

no_update() = NoUpdate()

mutable struct Layout
    component::Union{Nothing, Component}
end

"""
    struct DashApp <: Any

Representation of Dash application
"""
struct DashApp
    name ::String
    layout ::Layout
    callbacks ::Dict{Symbol, Callback}
    external_stylesheets ::Vector{String}
    external_scripts ::Vector{String}
    url_base_pathname ::String
    assets_folder ::String
    callable_components ::Dict{Symbol, Component}
    function DashApp(name::String;
        external_stylesheets ::Vector{String} = Vector{String}(),
        external_scripts ::Vector{String} = Vector{String}(),
        url_base_pathname="/",
        assets_folder::String = "assets")

        new(name, Layout(nothing), Dict{Symbol, Callback}(),
            external_stylesheets, external_scripts,
            url_base_pathname, assets_folder,
            Dict{Symbol, Component}()

            )
    end
end

function layout!(app::DashApp, component::Component)
    Base.getfield(app, :layout).component = component
    Components.collect_with_ids!(app.layout, app.callable_components)
end

getlayout(app::DashApp) = Base.getfield(app, :layout).component

function Base.setproperty!(app::DashApp, name::Symbol, value)
    name == :layout ? layout!(app, value) : Base.setfield!(app, name, value)
end

function Base.getproperty(app::DashApp, name::Symbol)
    name == :layout ? getlayout(app) : Base.getfield(app, name)
end

"""
    dash(name::String; external_stylesheets ::Vector{String} = Vector{String}(), url_base_pathname::String="/")
    dash(layout_maker::Function, name::String; external_stylesheets ::Vector{String} = Vector{String}(), url_base_pathname::String="/")

Construct a dash app using callback for layout creation

# Arguments
- `layout_maker::Function` - function for layout creation. Must has signature ()::Component
- `name::String` - Dashboard name
- `external_stylesheets::Vector{String} = Vector{String}()` - vector of external css urls
- `external_scripts::Vector{String} = Vector{String}()` - vector of external js scripts urls
- `url_base_pathname::String="/"` - base url path for dashboard, default "/"
- `assets_folder::String` - a path, relative to the current working directory,
for extra files to be used in the browser. Default `"assets"`

# Examples
```jldoctest
julia> app = dash("Test") do
    html_div() do
        html_h1("Test Dashboard")
    end
end
```
"""

function dash(name::String;
    external_stylesheets ::Vector{String} = Vector{String}(),
    external_scripts ::Vector{String} = Vector{String}(),
    url_base_pathname="/",
    assets_folder::String = "assets")
        result = DashApp(name,
            external_stylesheets=external_stylesheets,
            external_scripts=external_scripts,
            url_base_pathname=url_base_pathname,
            assets_folder = assets_folder
            )
    return result
end

function dash(layout_maker::Function, name::String;
      external_stylesheets ::Vector{String} = Vector{String}(),
      external_scripts ::Vector{String} = Vector{String}(),
      url_base_pathname="/",
      assets_folder::String = "assets")
    result = dash(name,
        external_stylesheets=external_stylesheets,
        external_scripts=external_scripts,
        url_base_pathname=url_base_pathname,
        assets_folder = assets_folder
        )
    layout!(result, layout_maker())
    return result
end

idprop_string(idprop::IdProp) = "$(idprop[1]).$(idprop[2])"

function check_idprop(app::DashApp, id::IdProp)
    if !haskey(app.callable_components, id[1])
        error("The layout havn't component with id `$(id[1])]`")
    end
    if !is_prop_available(app.callable_components[id[1]], id[2])
        error("The component with id `$(id[1])` havn't property `$(id[2])``")
    end
end

function output_string(id::CallbackId)
    if length(id.output) == 1
        return idprop_string(id.output[1])
    end
    return ".." *
    join(map(idprop_string, id.output), "...") *
    ".."
end

"""
    callback!(func::Function, app::Dash, id::CallbackId; pass_changed_props = false)

Create a callback that updates the output by calling function `func`.

If `pass_changed_props` is true then the first argument of callback is an array of changed properties

# Examples

```julia
app = dash("Test") do
    html_div() do
        dcc_input(id="graphTitle", value="Let's Dance!", type = "text"),
        dcc_input(id="graphTitle2", value="Let's Dance!", type = "text"),
        html_div(id="outputID"),
        html_div(id="outputID2")

    end
end
callback!(app, CallbackId(
    state = [(:graphTitle, :type)],
    input = [(:graphTitle, :value)],
    output = [(:outputID, :children), (:outputID2, :children)]
    )
    ) do stateType, inputValue
    return (stateType * "..." * inputValue, inputValue)
end
```

You can use macro `callid` string macro for make CallbackId :

```julia
callback!(app, callid"{graphTitle.type} graphTitle.value => outputID.children, outputID2.children") do stateType, inputValue

    return (stateType * "..." * inputValue, inputValue)
end
```

Using `changed_props`

```julia
callback!(app, callid"graphTitle.value, graphTitle2.value => outputID.children", pass_changed_props = true) do changed, input1, input2
    if "graphTitle.value" in changed
        return input1
    else
        return input2
    end
end
```

"""
function callback!(func::Function, app::DashApp, id::CallbackId; pass_changed_props = false)

    check_callback(func, app, id, pass_changed_props)

    out_symbol = Symbol(output_string(id))

    push!(app.callbacks, out_symbol => Callback(func, id, pass_changed_props))
end

function check_callback(func::Function, app::DashApp, id::CallbackId, pass_changed_props)

    isempty(id.input) && error("The callback method requires that one or more properly formatted inputs are passed.")

    length(id.output) != length(unique(id.output)) && error("One or more callback outputs have been duplicated; please confirm that all outputs are unique.")

    for out in id.output
        if any(x->out in x.id.output, values(app.callbacks))
            error("output \"$(out)\" already registered")
        end
    end

    foreach(x->check_idprop(app,x), id.state)
    foreach(x->check_idprop(app,x), id.input)
    foreach(x->check_idprop(app,x), id.output)

    args_count = length(id.state) + length(id.input)
    pass_changed_props && (args_count+=1)

    !hasmethod(func, NTuple{args_count, Any}) && error("Callback function don't have method with proper arguments")

    for id_prop in id.input
        id_prop in id.output && error("Circular input and output arguments were found. Please verify that callback outputs are not also input arguments.")
    end
end


text:
using SimpleSolvers
using Test

include("optimizers_problems.jl")

n = 1
x = ones(n)
nl = QuasiNewtonOptimizer(x, F)

@test config(nl) == nl.config
@test status(nl) == nl.status

solve!(x, nl)
println(status(nl))


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 = "Qt6ShaderTools"
version = v"6.3.0"

# Set this to true first when updating the version. It will build only for the host (linux musl).
# After that JLL is in the registyry, set this to false to build for the other platforms, using
# this same package as host build dependency.
const host_build = false

# Collection of sources required to build qt6
sources = [
    ArchiveSource("https://download.qt.io/official_releases/qt/$(version.major).$(version.minor)/$version/submodules/qtshadertools-everywhere-src-$version.tar.xz",
                  "5498959b9d37c254bd126ab1320cf86aeb8a31b3ea7e51db666a8f9698afbd6c"),
    ArchiveSource("https://sourceforge.net/projects/mingw-w64/files/mingw-w64/mingw-w64-release/mingw-w64-v10.0.0.tar.bz2",
                  "ba6b430aed72c63a3768531f6a3ffc2b0fde2c57a3b251450dcf489a894f0894")
]

script = raw"""
cd $WORKSPACE/srcdir

mkdir build
cd build/
qtsrcdir=`ls -d ../qtshadertools-*`

case "$target" in

    x86_64-linux-musl*)
        cmake -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_FIND_ROOT_PATH=$prefix -DCMAKE_BUILD_TYPE=Release $qtsrcdir
    ;;

    *mingw*)
        cd $WORKSPACE/srcdir/mingw*/mingw-w64-headers
        ./configure --prefix=/opt/$target/$target/sys-root --enable-sdk=all --host=$target
        make install

        cd ../mingw-w64-crt/
        if [ ${target} == "i686-w64-mingw32" ]; then
            _crt_configure_args="--disable-lib64 --enable-lib32"
        elif [ ${target} == "x86_64-w64-mingw32" ]; then
            _crt_configure_args="--disable-lib32 --enable-lib64"
        fi
        ./configure --prefix=/opt/$target/$target/sys-root --enable-sdk=all --host=$target --enable-wildcard ${_crt_configure_args}
        make -j${nproc}
        make install

        cd ../mingw-w64-libraries/winpthreads
        ./configure --prefix=/opt/$target/$target/sys-root --host=$target --enable-static --enable-shared
        make -j${nproc}
        make install

        cd $WORKSPACE/srcdir/build
        cmake -DQT_HOST_PATH=$host_prefix -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_FIND_ROOT_PATH=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release $qtsrcdir
    ;;

    *)
        cmake -DQT_HOST_PATH=$host_prefix -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_FIND_ROOT_PATH=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release $qtsrcdir
    ;;

esac

cmake --build . --parallel ${nproc}
cmake --install .
install_license $WORKSPACE/srcdir/qt*-src-*/LICENSE.LGPL3
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
if host_build
    platforms = [Platform("x86_64", "linux",cxxstring_abi=:cxx11,libc="musl")]
else
    platforms = expand_cxxstring_abis(filter(!Sys.isapple, supported_platforms()))
    filter!(p -> arch(p) != "armv6l", platforms) # No OpenGL on armv6
    platforms_macos = [ Platform("x86_64", "macos"), Platform("aarch64", "macos") ]
end

# The products that we will ensure are always built
products = [
    LibraryProduct(["Qt6ShaderTools", "libQt6ShaderTools", "QtShaderTools"], :libqt6shadertools),
]

products_macos = [
    FrameworkProduct("QtShaderTools", :libqt6shadertools),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    HostBuildDependency("Qt6Base_jll"),
    Dependency("Qt6Base_jll"),
]

if !host_build
    push!(dependencies, HostBuildDependency("Qt6ShaderTools_jll"))
end

include("../../fancy_toys.jl")

@static if !host_build
    if any(should_build_platform.(triplet.(platforms_macos)))
        build_tarballs(ARGS, name, version, sources, script, platforms_macos, products_macos, dependencies; preferred_gcc_version = v"9", julia_compat="1.6")
    end
end
if any(should_build_platform.(triplet.(platforms)))
    build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"9", julia_compat="1.6")
end