text:
module MatrixDirichletTest

using Test
using ReactiveMP
using Distributions
using Random

@testset "MatrixDirichlet" begin

    @testset "common" begin
        @test MatrixDirichlet <: Distribution
        @test MatrixDirichlet <: ContinuousDistribution
        @test MatrixDirichlet <: MatrixDistribution

        @test value_support(MatrixDirichlet) === Continuous
        @test variate_form(MatrixDirichlet)  === Matrixvariate
    end

    @testset "vague" begin
        @test_throws MethodError vague(MatrixDirichlet)

        d1 = vague(MatrixDirichlet, 3)

        @test typeof(d1) <: MatrixDirichlet
        @test mean(d1) == ones(3, 3) ./ sum(ones(3, 3), dims = 1)

        d2 = vague(MatrixDirichlet, 4)

        @test typeof(d2) <: MatrixDirichlet
        @test mean(d2) == ones(4, 4) ./ sum(ones(4, 4), dims = 1)

        @test vague(MatrixDirichlet, 3, 3) == vague(MatrixDirichlet, (3, 3))
        @test vague(MatrixDirichlet, 4, 4) == vague(MatrixDirichlet, (4, 4))
        @test vague(MatrixDirichlet, 3, 4) == vague(MatrixDirichlet, (3, 4))
        @test vague(MatrixDirichlet, 4, 3) == vague(MatrixDirichlet, (4, 3))

        d3 = vague(MatrixDirichlet, 3, 4)

        @test typeof(d3) <: MatrixDirichlet
        @test mean(d3) == ones(3, 4) ./ sum(ones(3, 4), dims = 1)
    end

    @testset "entropy" begin
        @test entropy(MatrixDirichlet([ 1.0 1.0; 1.0 1.0; 1.0 1.0 ])) ≈ -1.3862943611198906
        @test entropy(MatrixDirichlet([ 1.2 3.3; 4.0 5.0; 2.0 1.1 ])) ≈ -3.1139933152617787
        @test entropy(MatrixDirichlet([ 0.2 3.4; 5.0 11.0; 0.2 0.6 ])) ≈ -11.444984495104693
    end

    @testset "logmean" begin
        @test logmean(MatrixDirichlet([ 1.0 1.0; 1.0 1.0; 1.0 1.0 ])) ≈ [-1.5000000000000002 -1.5000000000000002; -1.5000000000000002 -1.5000000000000002; -1.5000000000000002 -1.5000000000000002]
        @test logmean(MatrixDirichlet([ 1.2 3.3; 4.0 5.0; 2.0 1.1 ])) ≈ [-2.1920720408623637 -1.1517536610071326; -0.646914475838374 -0.680458481634953; -1.480247809171707 -2.6103310904778305]
        @test logmean(MatrixDirichlet([ 0.2 3.4; 5.0 11.0; 0.2 0.6 ])) ≈ [-6.879998107291004 -1.604778825293528; -0.08484054226701443 -0.32259407259407213; -6.879998107291004 -4.214965875553984]
    end

    @testset "prod" begin
        d1 = MatrixDirichlet([ 0.2 3.4; 5.0 11.0; 0.2 0.6 ])
        d2 = MatrixDirichlet([ 1.2 3.3; 4.0 5.0; 2.0 1.1 ])
        d3 = MatrixDirichlet([ 1.0 1.0; 1.0 1.0; 1.0 1.0 ])

        @test prod(ProdAnalytical(), d1, d2) == MatrixDirichlet([0.3999999999999999 5.699999999999999; 8.0 15.0; 1.2000000000000002 0.7000000000000002])
        @test prod(ProdAnalytical(), d1, d3) == MatrixDirichlet([0.19999999999999996 3.4000000000000004; 5.0 11.0; 0.19999999999999996 0.6000000000000001])
        @test prod(ProdAnalytical(), d2, d3) == MatrixDirichlet([1.2000000000000002 3.3; 4.0 5.0; 2.0 1.1])
    end

    @testset "promote_variate_type" begin
        @test_throws MethodError promote_variate_type(Univariate, MatrixDirichlet)

        @test promote_variate_type(Multivariate, Dirichlet)  === Dirichlet
        @test promote_variate_type(Matrixvariate, Dirichlet) === MatrixDirichlet

        @test promote_variate_type(Multivariate, MatrixDirichlet)  === Dirichlet
        @test promote_variate_type(Matrixvariate, MatrixDirichlet) === MatrixDirichlet
    end

end

end


text:
using Test
using ClimaCoupler: Clock, tick!, stop_time_exceeded

@testset "Clock" begin
    time_info = (start = 0.0, dt = 0.5, stop = 2.0)
    clock = Clock(time_info...)

    tick!(clock)
    @test clock.time == time_info.dt
    while !stop_time_exceeded(clock)
        tick!(clock)
    end
    @test clock.time == time_info.stop
end

ext: jl

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
using Pkg, BinaryBuilder, SHA, Dates
if !isdefined(Pkg, :Artifacts)
    error("This must be run with Julia 1.3+!")
end
using Pkg.Artifacts, Pkg.PlatformEngines

include("../common.jl")

# Metadata
name = "Rootfs"
version = VersionNumber("$(year(today())).$(month(today())).$(day(today()))")
verbose = "--verbose" in ARGS

# We begin by downloading the alpine rootfs and using THAT as a bootstrap rootfs.
rootfs_url = "https://github.com/alpinelinux/docker-alpine/raw/v3.12/x86_64/alpine-minirootfs-3.12.0-x86_64.tar.gz"
rootfs_hash = "0beb54cf9bf69d085f9fcd291ff28b3335184d08b706d535f425e8180851edc9"
mkpath(joinpath(@__DIR__, "build"))
mkpath(joinpath(@__DIR__, "products"))
rootfs_targz_path = joinpath(@__DIR__, "build", "rootfs.tar.gz")
download_verify(rootfs_url, rootfs_hash, rootfs_targz_path; verbose=verbose, force=true)

# Unpack the rootfs (using `tar` on the local machine), then pack it up again (again using tools on the local machine) and squashify it:
rootfs_extracted = joinpath(@__DIR__, "build", "rootfs_extracted")
rm(rootfs_extracted; recursive=true, force=true)
mkpath(rootfs_extracted)
success(`tar -C $(rootfs_extracted) -zxf $(rootfs_targz_path)`)

# In order to launch our rootfs, we need at bare minimum, a sandbox/docker entrypoint.  Ensure those exist.
# We check in `sandbox`, but if it gets modified, we would like it to be updated.  In general, we hope
# to build it inside of BB, but in the worst case, where we've made an incompatible change and the last
# checked-in version of `sandbox` won't do, we need to be able to build it and build it now.  So if
# we're on a compatible platform, and `sandbox` is outdated, build it.
sandbox_path = joinpath(@__DIR__, "bundled", "utils", "sandbox")
if is_outdated(sandbox_path, "$(sandbox_path).c")
    try
        build_platform = platform_key_abi(String(read(`gcc -dumpmachine`)))
        @assert isa(build_platform, Linux)
        @assert arch(build_platform) == :x86_64
        verbose && @info("Rebuilding sandbox for initial bootstrap...")
        success(`gcc -O2 -static -static-libgcc -o $(sandbox_path) $(sandbox_path).c`)
    catch
        if isfile(sandbox_path)
            @warn("Sandbox outdated and we can't build it, continuing, but be warned, initial bootstrap might fail!")
        else
            error("Sandbox missing and we can't build it!  Build it somewhere else via `gcc -O2 -static -static-libgcc -o sandbox sandbox.c`")
        end
    end
end
libc_paths = [joinpath(@__DIR__, "bundled", "libs", "libc", "$(libc)-$(arch)") for libc in ("musl", "glibc") for arch in ("i686", "x86_64")]
download_libcs = joinpath(@__DIR__, "bundled", "libs", "libc", "download_libcs.sh")
if any(is_outdated.(libc_paths, download_libcs))
    success(`$download_libcs`)
end

# Copy in `sandbox` and `docker_entrypoint.sh` since we need those just to get up in the morning.
# Also set up a DNS resolver, since that's important too.  Yes, this work is repeated below, but
# we need a tiny world to set up our slightly larger world inside of.
verbose && @info("Constructing barebones bootstrap RootFS shard...")
mkpath(joinpath(rootfs_extracted, "overlay_workdir"))
cp(sandbox_path, joinpath(rootfs_extracted, "sandbox"); force=true)
cp(joinpath(@__DIR__, "bundled", "utils", "docker_entrypoint.sh"), joinpath(rootfs_extracted, "docker_entrypoint.sh"); force=true)
open(joinpath(rootfs_extracted, "etc", "resolv.conf"), "w") do io
    for resolver in ["8.8.8.8", "8.8.4.4", "4.4.4.4", "1.1.1.1"]
        println(io, "nameserver $(resolver)")
    end
end

# we really like bash, and it's annoying to have to polymorphise, so just lie for the stage 1 bootstrap
cp(joinpath(@__DIR__, "bundled", "utils", "fake_bash.sh"), joinpath(rootfs_extracted, "bin", "bash"); force=true)
cp(joinpath(@__DIR__, "bundled", "conf", "profile"), joinpath(rootfs_extracted, "etc", "profile"); force=true)
cp(joinpath(@__DIR__, "bundled", "conf", "profile.d"), joinpath(rootfs_extracted, "etc", "profile.d"); force=true)
rootfs_unpacked_hash, rootfs_squashfs_hash = generate_artifacts(rootfs_extracted, name, version)

# Slip these barebones rootfs images into our BB install location, overwriting whatever Rootfs shard would be chosen:
verbose && @info("Binding barebones bootstrap RootFS shards...")

insert_compiler_shard(name, version, rootfs_unpacked_hash, :unpacked)
insert_compiler_shard(name, version, rootfs_squashfs_hash, :squashfs)
Core.eval(BinaryBuilder.BinaryBuilderBase, :(bootstrap_list = Symbol[:rootfs]))

# PHWEW.  Okay.  Now, we do some of the same steps over again, but within BinaryBuilder, where
# we can actulaly run tools inside of the rootfs (e.g. if we're building on OSX through docker)

# Sources we build from
sources = [
    ArchiveSource(rootfs_url, rootfs_hash),
    # Objconv is very useful
    ArchiveSource("https://github.com/staticfloat/objconv/archive/v2.49.tar.gz",
                  "5fcdf0eda828fbaf4b3d31ba89b5011f649df3a7ef0cc7520d08fe481cac4e9f"),
    # As is patchelf
    GitSource("https://github.com/NixOS/patchelf.git",
              "e1e39f3639e39360ceebb2f7ed533cede4623070"),
    # We need a very recent version of meson to build gtk stuffs, so let's just grab the latest
    ArchiveSource("https://github.com/mesonbuild/meson/releases/download/0.52.0/meson-0.52.0.tar.gz",
                  "d60f75f0dedcc4fd249dbc7519d6f3ce6df490033d276ef1cf27453ef4938d32"),
    # And also our own local patches, utilities, etc...
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
set -x
# Get build tools ready. Note that they do not pollute the eventual Rootfs image;
# they are only within this currently-running, ephemeral, pocket universe
apk add build-base curl autoconf automake linux-headers gawk python3 py3-setuptools bison git

# $prefix is our chroot under construction
mv bin dev etc home lib media mnt proc root run sbin srv sys tmp usr var $prefix/
cd $prefix

# Setup DNS resolution
printf '%s\n' \
    "nameserver 8.8.8.8" \
    "nameserver 8.8.4.4" \
    "nameserver 4.4.4.4" \
    > etc/resolv.conf

# Insert system mountpoints
touch ./dev/null
touch ./dev/ptmx
touch ./dev/urandom
mkdir ./dev/pts
mkdir ./dev/shm

## Install foundational packages within the chroot
NET_TOOLS="curl wget git openssl ca-certificates"
MISC_TOOLS="python2 python3 py3-pip sudo file libintl patchutils grep zlib"
FILE_TOOLS="tar zip unzip xz findutils squashfs-tools unrar rsync"
INTERACTIVE_TOOLS="bash gdb vim nano tmux strace"
BUILD_TOOLS="make patch gawk autoconf automake libtool bison flex pkgconfig cmake ninja ccache"
apk add --update --root $prefix ${NET_TOOLS} ${MISC_TOOLS} ${FILE_TOOLS} ${INTERACTIVE_TOOLS} ${BUILD_TOOLS}

# chgrp and chown should be no-ops since we run in a single-user mode
rm -f ./bin/chown ./bin/chgrp
touch ./bin/chown ./bin/chgrp
chmod +x ./bin/chown ./bin/chgrp

# usr/libexec/git-core takes up a LOT of space because it uses hardlinks; convert to symlinks:
echo "Replacing hardlink versions of git helper commands with symlinks..."
sha() { sha256sum "$1" | awk '{ print $1 }'; }
GIT_SHA=$(sha ./usr/libexec/git-core/git)
for f in ./usr/libexec/git-core/*; do
    if [[ -f "${f}" ]] && [[ "$(basename ${f})" != "git" ]] && [[ $(sha "${f}") == ${GIT_SHA} ]]; then
        ln -svf git "${f}"
    fi
done

# Install utilities we'll use.  Many of these are compatibility shims, look at
# the files themselves to discover why we use them.
mkdir -p ./usr/local/bin ./usr/local/share/configure_scripts
cp -vd ${WORKSPACE}/srcdir/utils/tar_wrapper.sh ./usr/local/bin/tar
cp -vd ${WORKSPACE}/srcdir/utils/update_configure_scripts.sh ./usr/local/bin/update_configure_scripts
cp -vd ${WORKSPACE}/srcdir/utils/flagon ./usr/local/bin/flagon
cp -vd ${WORKSPACE}/srcdir/utils/fake_uname.sh ./usr/bin/uname
mv ./sbin/sysctl ./sbin/_sysctl
cp -vd ${WORKSPACE}/srcdir/utils/fake_sysctl.sh ./sbin/sysctl
cp -vd ${WORKSPACE}/srcdir/utils/fake_sha512sum.sh ./usr/local/bin/sha512sum
cp -vd ${WORKSPACE}/srcdir/utils/dual_libc_ldd.sh ./usr/bin/ldd
cp -vd ${WORKSPACE}/srcdir/utils/atomic_patch.sh ./usr/local/bin/atomic_patch
cp -vd ${WORKSPACE}/srcdir/utils/install_license.sh ./usr/local/bin/install_license
cp -vd ${WORKSPACE}/srcdir/utils/config.* ./usr/local/share/configure_scripts/
chmod +x ./usr/local/bin/*

# Deploy configuration
cp -vd ${WORKSPACE}/srcdir/conf/nsswitch.conf ./etc/nsswitch.conf

# Clear out all previous profile stuff first
rm -f ${prefix}/etc/profile.d/*
cp -vd ${WORKSPACE}/srcdir/conf/profile ${prefix}/etc/
cp -vd ${WORKSPACE}/srcdir/conf/profile.d/* ${prefix}/etc/profile.d/

# Install vim configuration/.vimrc
cp -vd ${WORKSPACE}/srcdir/conf/vimrc ${prefix}/etc/vim/vimrc
mkdir -p ${prefix}/etc/vim/
git clone https://github.com/VundleVim/Vundle.vim ${prefix}/etc/vim/bundle/Vundle.vim
chroot ${prefix} vim -E -u /etc/vim/vimrc -c PluginInstall -c qall

# Put sandbox and docker entrypoint into the root, to be used as `init` replacements.
gcc -g -O2 -static -static-libgcc -o ${prefix}/sandbox $WORKSPACE/srcdir/utils/sandbox.c
cp -vd ${WORKSPACE}/srcdir/utils/docker_entrypoint.sh ${prefix}/docker_entrypoint.sh

# Build the BB service client
gcc -O2 -std=c99 -static -static-libgcc -o ${prefix}/bin/bb $WORKSPACE/srcdir/utils/bb.c

# Move over libc loaders from our bundled directory.  These have very strict location requirements,
# so they MUST exist in /lib and /lib64.  These were generated from the `Glibc` and `Musl` builders
# in this Rootfs directory.  Note that `musl` always goes to `/lib`, not `/lib64`
mkdir -p ${prefix}/lib64 ${prefix}/lib
cp -vd ${WORKSPACE}/srcdir/libs/libc/glibc-x86_64/* ${prefix}/lib64
cp -vd ${WORKSPACE}/srcdir/libs/libc/glibc-i686/* ${prefix}/lib
cp -vd ${WORKSPACE}/srcdir/libs/libc/musl-x86_64/* ${prefix}/lib
cp -vd ${WORKSPACE}/srcdir/libs/libc/musl-i686/* ${prefix}/lib

# Build/install meson
cd ${WORKSPACE}/srcdir/meson-*/
python3 setup.py build
python3 setup.py install --prefix=/usr --root="${prefix}"

# Build/install objconv
cd ${WORKSPACE}/srcdir/objconv*/
g++ -O2 -o ${prefix}/usr/bin/objconv src/*.cpp

# Build/install patchelf
cd ${WORKSPACE}/srcdir/patchelf*/
./bootstrap.sh
./configure --prefix=${prefix}/usr
make -j${nproc}
make install

# Create mount points for future bindmounts
mkdir -p ${prefix}/overlay_workdir ${prefix}/meta

# Sneak some test suites into the build environment, so we can do sanity checks easily
cp -vdR ${WORKSPACE}/srcdir/testsuite ${prefix}/usr/share/

# Add some licenses for quick-and-dirty licensing
cp -vdR ${WORKSPACE}/srcdir/licenses ${prefix}/usr/share/

## Cleanup
# We can never extract these files, because they are too fancy  :(
rm -rf ${prefix}/usr/share/terminfo

# These cause case sensitivity errors
rm -rf ${prefix}/usr/share/perl5/core_perl/pod
rm -rf ${prefix}/etc/vim/bundle/vim-colorschemes/colors/darkBlue.vim

# Cleanup .pyc/.pyo files as they're not redistributable
find ${prefix}/usr -type f -name "*.py[co]" -delete -or -type d -name "__pycache__" -delete
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [host_platform]

# The products that we will ensure are always built
products = Product[
    ExecutableProduct("sandbox", :sandbox, "./"),
    ExecutableProduct("bash", :bash),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarball
verbose && @info("Building full RootfS shard...")
ndARGS = filter(x -> !occursin("--deploy", x), ARGS)
build_info = build_tarballs(ndARGS, name, version, sources, script, platforms, products, dependencies; skip_audit=true)

# Upload the shards
if any(occursin.("--deploy", ARGS))
    upload_and_insert_shards("JuliaPackaging/Yggdrasil", name, version, build_info)
end

max_stars_repo_name: Micket/Yggdrasil

max_issues_repo_name: Micket/Yggdrasil

max_forks_repo_name: Micket/Yggdrasil


text:
# See documentation for `crawl_call` for documentation on what we're doing here.
function crawl(::Val{:call}, f, args...)
    return :(JSAST(
        :call,
        $(crawl(f)),
        $(crawl.(args)...),
    ))
end
function crawl(::Val{:call}, f::Symbol, args...)
    return crawl_call(f, args...)
end

"""
Crawl a _simple_ function call.

_Simple_ is defined to be a function call where the function name is a symbol.
In general, the function name may be any expression (e.g., in
`@js console.log("foo")`, the function `console.log` is not a symbol and has its
own JSAST since it's using the dot-operator). However, for simple cases, when
the function is just a symbol, we can dispatch to the `crawl_call` function to
perform further transformation if necessary.

This is useful for transforming syntactic constructs which are represented in
Julia as function calls (e.g. we transform `+(1, 2)` into `1 + 2`), as well as
for translating some Julia-isms.
"""
function crawl_call(f::Symbol, args...)
    return crawl_call(Val(f), args...)
end

# Default for generic methods
crawl_call(::Val{M}, args...) where {M} = :(
    JSAST(
        :call,
        $(crawl(M)),
        $(crawl.(args)...),
    )
)

function deparse(::Val{:call}, f, args...)
    jsstring(deparse(f), "(", join(deparse.(args), ", "), ")")
end


text:
using PyCall
import VersionParsing: vparse

const min_sparse_ir_version = "0.9.2"

function check_version_sparse_ir()
    sparse_ir = pyimport("sparse_ir")
    if vparse(sparse_ir.__version__) < vparse(min_sparse_ir_version)
        error("sparse-ir v$(sparse_ir.__version__) is found, but required version >=$(min_sparse_ir_version)! " *
              "PyCall was built with $(PyCall.python). " *
              "If you are using the Conda backend for PyCall, try \"using Conda;Conda.update()\" on the Julia REPL. " *
              "If not, try \"pip3 install -U sparse-ir\". ")
    end
end


text:
export PiecewiseUniformMassDistribution

"Initial 1D piecewise uniform mass distribution."
struct PiecewiseUniformMassDistribution <: InitialMassDistribution{1}
    number_of_particles::UInt
    boundary_positions::Vector{Float}
    mass_densities::Vector{Float}

    """
        PiecewiseUniformMassDistribution(
            exterior_boundary_positions::Bounds{1,<:Number},
            interior_boundary_positions::AbstractVector{Float},
            mass_densities::AbstractVector{Float},
            number_of_particles::Integer,
        )

    Create a new piecewise uniform mass distribution with the given mass densities
    between the given boundary positions, with the given number of particles.
    """
    function PiecewiseUniformMassDistribution(
        exterior_boundary_positions::Bounds{1,<:Number},
        interior_boundary_positions::AbstractVector{Float},
        mass_densities::AbstractVector{Float},
        number_of_particles::Integer,
    )
        @check number_of_particles > 0
        @checkeq(
            length(mass_densities),
            length(interior_boundary_positions) + 1
        )
        boundary_positions = [
            lower(exterior_boundary_positions)[]
            interior_boundary_positions
            upper(exterior_boundary_positions)[]
        ]
        @check all(diff(boundary_positions) .> 0)
        @check all(mass_densities .> 0)
        new(number_of_particles, boundary_positions, mass_densities)
    end
end

function (distribution::PiecewiseUniformMassDistribution)()
    m =
        sum(
            distribution.mass_densities .*
            diff(distribution.boundary_positions),
        ) / distribution.number_of_particles

    positions::Vector{SVectorF{1}} = []
    mass_densities::Vector{Float} = []
    for (ρ, lower_bound, upper_bound) in zip(
        distribution.mass_densities,
        distribution.boundary_positions[1:end-1],
        distribution.boundary_positions[2:end],
    )
        extent = upper_bound - lower_bound
        local_number_of_particles = UInt(round(ρ * extent / m))
        cell_extent = extent / local_number_of_particles
        start = lower_bound + 0.5 * cell_extent
        append!(
            positions,
            map(
                x -> SA_F64[x],
                range(
                    start,
                    step = cell_extent,
                    length = local_number_of_particles,
                ),
            ),
        )
        append!(mass_densities, fill(ρ, local_number_of_particles))
    end

    MassDistributionData(m, positions, mass_densities)
end

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
# test/runtests.jl

using LinearMapsAA
using Test: @test, detect_ambiguities

@test LinearMapAA(:test)

@test length(detect_ambiguities(LinearMapsAA)) == 0


text:
using Manopt, Manifolds, ManifoldsBase, Random, LinearAlgebra, BenchmarkTools
Random.seed!(42)

function run_rayleigh_experiment(n::Int)
    A = randn(n, n)
    A = (A + A') / 2
    F(X::Array{Float64,1}) = X' * A * X
    ∇F(X::Array{Float64,1}) = 2 * (A * X - X * X' * A * X)
    M = Sphere(n - 1)
    x = random_point(M)
    return quasi_Newton(
        M,
        F,
        ∇F,
        x;
        #memory_size=-1,
        stopping_criterion=StopWhenAny(
            StopAfterIteration(max(1000)), StopWhenGradientNormLess(10^(-6))
        ),
        debug=[:Iteration, " ", :Cost, "\n", 1, :Stop],
    )
end
io = IOBuffer()

for n in [100]
    b = @benchmark run_rayleigh_experiment($n) samples = 30
    show(io, "text/plain", b)
    s = String(take!(io))
    println("Benchmarking $(n):\n", s, "\n\n")
end


text:
# # [How to use `@?`](@id tutorial)

using Maybe
using Maybe: @something

using LiterateTest                                                     #src
using Test                                                             #src

# ## Introduction

# Julia `Base` provides functions `findfirst` and `findlast` that
# returns an integer when an element is found:

findfirst(x -> gcd(x, 42) == 21, 50:200)

# and `nothing` if not:

@assert findlast(x -> gcd(x, 42) == 23, 50:200) === nothing

# It is rather tedious to combine such functions:

function find_some_random_range_1(data)
    f(x) = gcd(x, 42) == 21
    i = findfirst(f, data)
    i === nothing && return nothing
    j = findlast(f, data)
    j === nothing && return nothing
    return data[i:j]
end

@assert find_some_random_range_1(50:200) === 63:189
@assert find_some_random_range_1(30:50) === nothing

# To solve this issue, Maybe.jl provides a macro [`@?`](@ref) that
# lets you write

function find_some_random_range_2(data)
    f(x) = gcd(x, 42) == 21
    @? begin
        i = findfirst(f, data)
        j = findlast(f, data)
        return data[i:j]
    end
end

@assert find_some_random_range_2(50:200) === Some(63:189)
@assert find_some_random_range_2(30:50) === nothing

# Similarly, `@?` is also useful for indexing into arrays,
# dictionaries, etc. that may fail.  For example

dict = Dict(:a => 1, :b => nothing, :c => 2)
@assert (@? dict[:a] + dict[:c]) == Some(3)
@assert (@? dict[:a] + dict[:non_existing_key]) === nothing

# This is explained in more details in [Indexing section](@ref
# tutorials-lift-macro-indexing) below.

# ## How it works

# The above example `find_some_random_range_2` is roughly equivalent
# to

function find_some_random_range_3(data)
    f(x) = gcd(x, 42) == 21
    mi = findfirst(f, data)
    mi === nothing && return nothing  # (1)
    i = something(mi)                 # (2)
    mj = findlast(f, data)
    mj === nothing && return nothing  # (1)
    j = something(mj)                 # (2)
    md = Maybe.getindex(data, i:j)    # (3)
    md === nothing && return nothing  # (1)
    d = something(md)                 # (2′)
    return Some(d)                    # (4)
end

@assert find_some_random_range_3(50:200) === Some(63:189)
@assert find_some_random_range_3(30:50) === nothing

# Observe that:
#
# (1) If a function returns `nothing`, the whole evaluation
# short-circuits and evaluates to `nothing`.  (Side notes:
# short-circuiting is not actually implemented using `return` in `@?`
# so that it can be used outside functions.)
#
# (2) The returned value is always unwrapped by `something`.  Thus, it
# works with "ordinary" functions like `+` as well as a function
# returning `Some` (like `Maybe.getindex`); see (2′).
#
# (3) Indexing dispatches to [`Maybe.getindex`](@ref).
#
# (4) Finally, the returned result is always ensured to be wrapped by
# `Some`.

# ## More examples

# Consider a function that returns `nothing` on "failure":

maybe_positive(x) = x > 0 ? x : nothing;

# When a function call in `@?` is evaluated to a non-`nothing`, the
# returned value is wrapped in `Some`:

@test begin
    @? maybe_positive(1)
end === Some(1)

# When `Some` appears in the argument positions, they are
# automatically un-wrapped:

@test begin
    @? maybe_positive(1) + 1
end === Some(2)

# `@?` is evaluated to `nothing` when the first sub-expression is
# evaluated to `nothing`:

@test begin
    global r = @? maybe_positive(-1) + 1
end === nothing
@assert r === nothing

# Literal `nothing` is transformed to `Some(nothing)`:

@test begin
    @? nothing
end === Some(nothing)

# `@?` terminates the call chain immediately when it sees `nothing` as
# the return value.

ARG_HISTORY = []
demo(label, x) = (push!(ARG_HISTORY, label => x); x)

@evaltest "r = @? demo(2, identity(demo(1, nothing)))" begin
    @test ans === nothing
    @test ARG_HISTORY == [1 => nothing]
end
@assert r === nothing

# Note that `demo(2, ...)` is not called:

ARG_HISTORY

# This can be avoided by prefixing the function name by `$`:

empty!(ARG_HISTORY)

@evaltest raw"@? $demo(2, $identity($demo(1, nothing)))" begin
    @test ans === Some(nothing)
    @test ARG_HISTORY == [1 => nothing, 2 => nothing]
end

# Now `demo(2, ...)` is called:

ARG_HISTORY

# This is because `@?` automatically unwraps `Some` in the call chain.
# It means that `Some` acts like the identity function inside `@?`:

@test begin
    @? Some(Some(Some(nothing)))
end === Some(nothing)

# On the other hand, `identity` acts like "unwrap" function:

@test begin
    global r = @? identity(identity(Some(nothing)))
end === nothing
@assert r === nothing

# Finally, with `$Some`:

@test begin
    @? $Some($Some($Some(nothing)))
end === Some(Some(Some(Some(nothing))))

# `$(expression)` can be used to evaluate the whole `expression` in
# the normal context

@test begin
    @? $(Some(Some(Some(nothing))))
end == Some(Some(Some(Some(nothing))))

# ## [Indexing](@id tutorials-lift-macro-indexing)

# Index access is also lifted.  Consider a dictionary

dict = Dict(:a => Dict(:b => nothing, :c => 2));

# Any value stored in the container (here, a `Dict`) is returned as a
# `Some`:

@test begin
    @? dict[:a][:c]
end === Some(2)

# Thus, `nothing` stored in the container is returned as
# `Some(nothing)`

@test begin
    @? dict[:a][:b]
end === Some(nothing)

# On the other hand, accessing non-existing index returns `nothing`:

@assert (@? dict[:a][:d]) === nothing
@assert (@? dict[:d]) === nothing

# Index access can be "fused" with other function calls:

@assert (@? dict[:a][:c] + 1) === Some(3)
@assert (@? dict[:a][:d] + 1) === nothing

# Existing `nothing` value can be normalized using `something` as
# usual inside `@?`.  This is because `@?` automatically unwraps
# `Some` once.

@assert (@? something(dict[:a][:b], 0) + 1) === Some(1)
@assert (@? something(dict[:a][:c], 0) + 1) === Some(3)
@assert (@? something(dict[:a][:d], 0) + 1) === nothing

# Since `identity` acts like unwrapping operation inside `@?`, it can
# be used for normalizing non-existing key and existing `nothing` to
# the same value:

@assert something((@? identity(dict[:a][:b])), 0) + 1 === 1
@assert something((@? identity(dict[:a][:c])), 0) + 1 === 3
@assert something((@? identity(dict[:a][:d])), 0) + 1 === 1

# Indexing also works with arrays

vectors = [[1, 2], [3, 4, 5]]

@assert (@? vectors[1][2]) === Some(2)
@assert (@? vectors[3][4]) === nothing
@assert (@? vectors[1][3]) === nothing

# ## `@? return`

# `return` in `@?` is a powerful pattern for returning a non-`nothing`
# value.

function first_something(dict)
    @? return dict[:a]
    @? return dict[:b]
    @? return dict[:c]
    return nothing
end

@assert first_something(Dict()) == nothing
@assert first_something(Dict(:c => 3)) == Some(3)
@assert first_something(Dict(:a => 1, :c => 3)) == Some(1)

# Use `$return(x)` to avoid the returned value to be automatically
# wrapped by `Some`.  Note that the parentheses are required.

function first_something2(dict)
    @? $return(dict[:a])
    @? $return(dict[:b])
    @? $return(dict[:c])
    return nothing
end

@assert first_something2(Dict()) == nothing
@assert first_something2(Dict(:c => 3)) == 3
@assert first_something2(Dict(:a => 1, :c => 3)) == 1

# Alternatively, combined with [`@something`](@ref):

function first_something3(dict)
    return @something {
        @? dict[:a];
        @? dict[:b];
        @? dict[:c];
        0;  # fallback
    }
end

@assert first_something3(Dict()) === 0
@assert first_something3(Dict(:c => 3)) === 3
@assert first_something3(Dict(:a => 1, :c => 3)) === 1

# Note that `@?` can have multiple statements.  Any of the
# sub-expression evaluating to `nothing` short-circuits to the end of
# `@?` block.

function a_plus_b(dict)
    @? begin
        a = dict[:a]
        b = dict[:b]
        return a + b
    end
    return 0
end

@assert a_plus_b(Dict(:a => 1)) == 0
@assert a_plus_b(Dict(:a => 1, :b => 2)) == Some(3)

# Note that `@?` works with other control flows like `break`:

function something_positive_add3(xs, idx)
    found = nothing
    for i in idx
        @? if xs[i] > 0  # out-of-bound access is ignored
            found = i
            break
        end
    end
    return @? xs[found] + 3
end

@assert something_positive_add3([-1, 0, 1], [0, 1, 3]) === Some(4)
@assert something_positive_add3([-1, 0, 1], [0, 1, 2]) === nothing

# ## Combining `@?` and `@something`

# As explained, using `@?` in `@something` is well supported.
# However, further nesting is not supported (as the same expression
# would be processed by `@?` twice).  Assigning to an intermediate
# variable is a safe way to use the result of `@something` in `@?`.
# The pattern `@something(..., return)` is useful inside functions for
# this.

function extract_a_and_bc(x)
    c = @something {
        @? x[:b][:c];
        @? x[:b][:ccc];
        return;  # filter out if none of them exist
    }
    return @? (a = x[:a], c = c)
end

@assert extract_a_and_bc(Dict(:a => 1, :b => Dict(:c => 2))) === Some((a = 1, c = 2))
@assert extract_a_and_bc(Dict(:a => 1)) === nothing
@assert extract_a_and_bc(Dict(:a => 1, :b => Dict())) === nothing
@assert extract_a_and_bc(Dict(:b => Dict(:c => 2))) === nothing
@assert extract_a_and_bc(Dict(:a => 10, :b => Dict(:ccc => 20))) === Some((a = 10, c = 20))

# In a rare situation, it may be useful to use `$(...)` to nest
# `@something`-of-`@?`s in `@?`:

function a_plus_b_plus_c_or_d_times_2(x)
    @? begin
        p = x[:a] + x[:b]
        q = $(@something {
            @? p + x[:c];
            @? p + x[:d];
            Some(nothing);  # `return` works as well
        })
        return 2q
    end
end

@assert a_plus_b_plus_c_or_d_times_2(Dict(:a => 1, :b => 2, :c => 3)) == Some(12)
@assert a_plus_b_plus_c_or_d_times_2(Dict(:a => 1, :b => 2, :d => 3)) == Some(12)
@assert a_plus_b_plus_c_or_d_times_2(Dict(:a => 1, :b => 2)) === nothing

# ## Functions

# When `@?` sees functions (including closures and `do` blocks ), it
# converts them recursively to operate on `Union{Some{T},Nothing}`.

@? get_a_plus_b(dict) =
    get(dict, :a_plus_b) do
        dict[:a] + dict[:b]
    end

@assert get_a_plus_b(Dict(:a_plus_b => 1)) == Some(1)
@assert get_a_plus_b(Dict(:a => 1, :b => 2)) == Some(3)
@assert get_a_plus_b(Dict(:b => 1)) == nothing

# Thus, functions created with `@?` work nicely in `@?`:

@test begin
    @? begin
        x = get_a_plus_b(Dict(:a_plus_b => 1))
        y = get_a_plus_b(Dict(:a => 1, :b => 2))
        x + y
    end
end === Some(4)

@testset "get_a_plus_b" begin
    @test begin
        @? begin
            x = get_a_plus_b(Dict(:a_plus_b => 1))
            y = get_a_plus_b(Dict(:a => 1))
            x + y
        end
    end == nothing
    @test begin
        @? begin
            x = get_a_plus_b(Dict())
            y = get_a_plus_b(Dict(:a => 1, :b => 2))
            x + y
        end
    end == nothing
end

# ## Debugging

# You can include debugging information to `@? expr` with `:debug`
# flag as in `@? :debug expr`.  This inserts `@debug` logging
# statements for every possible short-circuit points.  When using the
# standard Julia logger, this information can be printed by setting
# the environment variable:
#
# ```julia
# ENV["JULIA_DEBUG"] = "all"  # or narrower scope like "MyPackage"
# ```
#
# This logging statement prints the expression that is evaluated to
# `nothing` and all the local variables using `Base.@locals`.

#! format: off
using Logging                                                        # hide
OUTPUT = let stdout = IOBuffer()                                       #src
withenv("JULIA_DEBUG" => "all") do                                   # hide
with_logger(SimpleLogger(stdout, Logging.Debug)) do                  # hide
dict = Dict(:a => Dict(:b => nothing, :c => 2))
@? :debug dict[:a][:d]
end                                                                  # hide
end                                                                  # hide
    String(take!(stdout))                                              #src
end                                                                    #src
#! format: on

@testset begin
    @test occursin("evaluating = (Maybe.getindex)(dict[:a], :d)", OUTPUT)
end

# `@? :debug expr` also inserts a call to a no-op function
# `Maybe._break`.  When using Debugger.jl, the state just before exit
# can be examined by adding it to the breakpoint with `bp add
# Maybe._break`.
#
# ```julia
# julia> using Maybe, Debugger
#
# julia> dict = Dict(:a => Dict(:b => nothing, :c => 2));
#
# julia> f(dict) = @? :debug dict[:a][:d];
#
# julia> @enter f(dict)
# In f(dict) at REPL[4]:1
# >1  f(dict) = @? :debug dict[:a][:d]
#
# About to run: return Maybe.Implementations.nothing
# 1|debug> bp add Maybe._break
# [ Info: added breakpoint for function _break
# 1] _break
#
# 1|debug> c
# Hit breakpoint:
# failed to lookup source code, showing lowered code:
# In _break() at /home/takafumi/.julia/dev/Maybe/src/lift.jl:51
# >1  1 ─     return Maybe.Implementations.nothing
#
# About to run: $(Expr(:meta, :noinline))
# 1|debug> up
# In f(dict) at REPL[4]:1
# >1  f(dict) = @? :debug dict[:a][:d]
#
# About to run: (Maybe._break)()
# 2|debug>
# ```


text:
#=
Doing dot product with a layer of neurons and multiple inputs
Associated YT NNFS tutorial: https://www.youtube.com/watch?v=tMrbN67U9d4
=#

inputs = [1.0, 2.0, 3.0, 2.5]
weights = [ 0.2   0.8  -0.5   1.0;
            0.5  -0.91  0.26 -0.5;
           -0.26 -0.27  0.17  0.87]

biases = [2.0, 3.0, 0.5]

output = weights * inputs + biases
println(output)

max_stars_repo_name: yarikbratashchuk/NNfSiX

max_issues_repo_name: yarikbratashchuk/NNfSiX

max_forks_repo_name: yarikbratashchuk/NNfSiX


text:
using Klimakoffer

NT = 48

mesh = Mesh()

co2_concentration = 315.0

model = Model(mesh, NT; co2_concentration = co2_concentration,compute_albedo=true)

discretization1 = Discretization(mesh, model, NT)

GlobTemp = compute_equilibrium!(discretization1)


text:
function system!(f::Function, k::Kosmo)
    type_tup = func_types(f)
    v = create_combinations(type_tup, k._component_pool)
    s = System(f, v, type_tup)
    push!(k._system_pool, s)
end

"""
Returns the Function's argument types in order.
"""
func_types(f::Function) = Tuple(first(methods(f)).sig.types[2:end])


text:
export add_propagator, inject!, extract

abstract type AbstractPropagator end

struct CompositePropagator <: AbstractPropagator
    propagators::Vector
end

const GLOBAL_PROPAGATOR = CompositePropagator([])

"""
    add_propagator(p::AbstractPropagator)

Set a new propagator as the global composite propagator.
"""
add_propagator(p::AbstractPropagator) = add_propagator(GLOBAL_PROPAGATOR, p)
add_propagator(cp::CompositePropagator, p::AbstractPropagator) = push!(cp.propagators, p)

Base.keys(p::CompositePropagator) = (k for x in p.propagators for k in keys(x))

"""
    inject(carrier, [global_propagator], [current_context])

Injects the value into a carrier. For example, into the headers of an HTTP request.
"""
function inject!(carrier, propagator::AbstractPropagator=GLOBAL_PROPAGATOR, ctx::Context=current_context()) end

function inject!(
    carrier,
    propagator::CompositePropagator,
    ctx::Context
)
    for p in propagator.propagators
        inject!(carrier, p, ctx)
    end
end

"""
    extract(carrier, [global_propagator], [current_context])

Extracts the value from an incoming request. For example, from the headers of an HTTP request.
"""
function extract(carrier, propagator::AbstractPropagator=global_propagator(), ctx::Context=current_context()) end

function extract(
    carrier,
    propagator::CompositePropagator,
    ctx::Context
)
    for p in propagator.propagators
        ctx = extract(carrier, p, ctx)
    end
    ctx
end


text:
#runs the basic panel form regressions
function olsretcontributionL13(df::AbstractDataFrame)
  specs::FMSpecs = FMSpecs(6)

  #set up the results
  push!(specs,
    xspec=Meta.parse("Llreturnwin"),
    xnames = [:intercept, :LagReturn],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)
  push!(specs,
    xspec=Meta.parse("Llreturn3yrwin"),
    xnames = [:intercept, :LagReturn3yr],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)
  push!(specs,
    xspec=Meta.parse("Llreturnnetbmsp5001yrwin"),
    xnames = [:intercept, :LagReturn],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)

  push!(specs,
    xspec=Meta.parse("Llreturnnetbmsp5003yrwin"),
    xnames = [:intercept, :LagReturn3yr],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)

  push!(specs,
    xspec=Meta.parse("Llreturnnetbmff31yrwin"),
    xnames = [:intercept, :LagReturn],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)

  push!(specs,
    xspec=Meta.parse("Llreturnnetbmff33yrwin"),
    xnames = [:intercept, :LagReturn3yr],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)

  #run the regressions
  computeFMLMresults!(df, specs)

  return specs
end

#runs the basic panel form regressions
function olscontributionregretcon(df::AbstractDataFrame)
  specs::FMSpecs = FMSpecs(6)

  #set up the results
  push!(specs,
    xspec=Meta.parse("Lpcontributions + Lpprogramexpense"),
    xnames = [:intercept, :LagContrib, :LagExp],
    clusterspec = :ein,
    withinspec = :ein,
    yspec=:lreturnwin)

  push!(specs,
    xspec=Meta.parse("Lpcontributions + Lpprogramexpense"),
    xnames = [:intercept, :LagContrib, :LagExp],
    clusterspec = :categoryfiltered,
    withinspec = :categoryfiltered,
    yspec=:lreturnwin)

  push!(specs,
  xspec=Meta.parse("Lpcontributions + Lpprogramexpense"),
  xnames = [:intercept, :LagContrib, :LagExp],
    clusterspec = :ein,
    withinspec = :ein,
    yspec=:lreturnnetbmsp5001yrwin)

  push!(specs,
  xspec=Meta.parse("Lpcontributions + Lpprogramexpense"),
  xnames = [:intercept, :LagContrib, :LagExp],
    clusterspec = :categoryfiltered,
    withinspec = :categoryfiltered,
    yspec=:lreturnnetbmsp5001yrwin)

  push!(specs,
  xspec=Meta.parse("Lpcontributions + Lpprogramexpense"),
  xnames = [:intercept, :LagContrib, :LagExp],
    clusterspec = :ein,
    withinspec = :ein,
    yspec=:lreturnnetbmff31yrwin)

  push!(specs,
  xspec=Meta.parse("Lpcontributions + Lpprogramexpense"),
  xnames = [:intercept, :LagContrib, :LagExp],
    clusterspec = :categoryfiltered,
    withinspec = :categoryfiltered,
    yspec=:lreturnnetbmff31yrwin)

  #run the regressions
  computeFMLMresults!(df, specs)

  return specs
end

#runs the basic panel form regressions
function olscontributionregconret(df::AbstractDataFrame)
  specs::FMSpecs = FMSpecs(6)

  #set up the results
  push!(specs,
    xspec=Meta.parse("Llreturnwin + Lpprogramexpense"),
    xnames = [:intercept, :LagReturn, :LagExp],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)
  push!(specs,
    xspec=Meta.parse("Llreturn3yrwin + Lpprogramexpense"),
    xnames = [:intercept, :LagReturn3yr, :LagExp],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)
  push!(specs,
    xspec=Meta.parse("Llreturnnetbmsp5001yrwin + Lpprogramexpense"),
    xnames = [:intercept, :LagReturn, :LagExp],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)

  push!(specs,
    xspec=Meta.parse("Llreturnnetbmsp5003yrwin + Lpprogramexpense"),
    xnames = [:intercept, :LagReturn3yr, :LagExp],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)

  push!(specs,
    xspec=Meta.parse("Llreturnnetbmff31yrwin + Lpprogramexpense"),
    xnames = [:intercept, :LagReturn, :LagExp],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)

  push!(specs,
    xspec=Meta.parse("Llreturnnetbmff33yrwin + Lpprogramexpense"),
    xnames = [:intercept, :LagReturn3yr, :LagExp],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)

  #run the regressions
  computeFMLMresults!(df, specs)

  return specs
end

function olsretcontributionLL13(df::AbstractDataFrame)
  specs::FMSpecs = FMSpecs(6)

  #set up the results
  push!(specs,
    xspec=Meta.parse("Llreturnwin + LLlreturnwin"),
    xnames = [:intercept, :LagReturn, :Lag2xReturn],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)

  push!(specs,
    xspec=Meta.parse("Llreturn3yrwin + L3Llreturn3yrwin"),
    xnames = [:intercept, :LagReturn3yr, :Lag2xReturn3yr],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)

  push!(specs,
    xspec=Meta.parse("Llreturnnetbmsp5001yrwin + LLlreturnnetbmsp5001yrwin"),
    xnames = [:intercept, :LagReturn, :Lag2xReturn],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)

  push!(specs,
    xspec=Meta.parse("Llreturnnetbmsp5003yrwin + L3Llreturnnetbmsp5003yrwin"),
    xnames = [:intercept, :LagReturn3yr, :Lag2xReturn3yr],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)

  push!(specs,
    xspec=Meta.parse("Llreturnnetbmff31yrwin + LLlreturnnetbmff31yrwin"),
    xnames = [:intercept, :LagReturn, :Lag2xReturn],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)

  push!(specs,
    xspec=Meta.parse("Llreturnnetbmff33yrwin + L3Llreturnnetbmff33yrwin"),
    xnames = [:intercept, :LagReturn3yr, :Lag2xReturn3yr],
    yspec=:pcontributions,
    clusterspec = :ein,
    withinspec = :ein)

  #run the regressions
  computeFMLMresults!(df, specs)

  return specs
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# Model
include_model("acrobot")

# Horizon
T = 101

# Time step
tf = 5.0
h = tf / (T - 1)

# ul <= u <= uu
ul, uu = control_bounds(model, T, -10.0, 10.0)

# Initial and final states
x1 = [0.0; 0.0; 0.0; 0.0]
xT = [π; 0.0; 0.0; 0.0]
xl, xu = state_bounds(model, T, x1 = x1, xT = xT)

# Objective
obj = quadratic_tracking_objective(
        [t < T ? Diagonal(1.0 * ones(model.n)) : Diagonal(ones(model.n)) for t = 1:T],
        [Diagonal(1.0e-1 * ones(model.m)) for t = 1:T-1],
        [xT for t = 1:T],
        [zeros(model.m) for t = 1:T])

# Problem
prob = trajectory_optimization_problem(model,
           obj,
           T,
           h = h,
           xl = xl,
           xu = xu,
           ul = ul,
           uu = uu)

# Trajectory initialization
x0 = linear_interpolation(x1, xT, T) # linear interpolation on state
u0 = random_controls(model, T, 0.001) # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

# Solve
@time z, info = solve(prob, copy(z0), tol = 1.0e-6, c_tol = 1.0e-6)

# Visualize
using Plots
x, u = unpack(z, prob)
plot(hcat(x...)', width = 2.0)
plot(hcat(u...)', width = 2.0, linetype = :steppost)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model, x, Δt = h)

# Model-predictive control
include_snopt()

function run_mpc()
    T_mpc = 25
    x_ref = [x...,[x[end] for i = 1:T_mpc]...]
    u_ref = [u...,[u[end] for i = 1:T_mpc]...]
    Q_mpc = Diagonal(100.0 * ones(model.n))
    R_mpc = Diagonal(0.1 * ones(model.m))
    K, P = tvlqr(model, x, u, h, [Q_mpc for t = 1:T], [R_mpc for t = 1:T-1])
    P_mpc = [P..., [P[end] for i = 1:T_mpc]...]

    model_sim = Acrobot{RK3, FixedTime}(4, 1, 4,
        1.0 + 0.05, 0.33, 1.0, 0.5, 1.0 - 0.05, 0.33, 1.0, 0.5, 9.81, 0.1, 0.1)

    shift = 1

    x_hist = [x_ref[shift] + 1.0e-1 * randn(model.n)]
    u_hist = []

    for t = 1:T-1
    	println("T = $t")
    	ul_mpc, uu_mpc = control_bounds(model, T_mpc, -10.0, 10.0)
    	xl_mpc, xu_mpc = state_bounds(model, T_mpc, x1 = x_hist[end])

    	# Objective
    	obj_mpc = quadratic_tracking_objective(
    	        [t < T_mpc ? Q_mpc : P_mpc[T_mpc + shift - 1] for t = 1:T_mpc],
    	        [R_mpc for t = 1:T_mpc-1],
    	        [x_ref[shift + t - 1] for t = 1:T_mpc],
    		[u_ref[shift + t - 1] for t = 1:T_mpc-1])

    	# Problem
    	prob_mpc = trajectory_optimization_problem(model,
    			obj_mpc,
    			T_mpc,
    			h = h,
    			ul = ul_mpc,
    			uu = uu_mpc,
    			xl = xl_mpc,
    			xu = xu_mpc)

    	# Pack trajectories into vector
    	z0_mpc = pack(x_ref[shift .+ (1:T_mpc)], u_ref[shift .+ (1:T_mpc-1)], prob_mpc)

    	# Solve
    	@time z_mpc, info = solve(prob_mpc, copy(z0_mpc),
            nlp = :SNOPT7, mapl = 0)
    	x_mpc, u_mpc = unpack(z_mpc, prob_mpc)

    	push!(x_hist,
    		fd(model_sim, x_hist[end], u_mpc[1], 1.0e-3 * randn(model_sim.d), h, nothing))
    	push!(u_hist, u_mpc[1])
    	shift += 1
    end
    return x_hist, u_hist
end

x_hist_mpc, u_hist_mpc = run_mpc()

plot(hcat(x...)',
	label = "", color = :black, width = 2.0)
plot!(hcat(x_hist_mpc...)',
	label = "", color = :magenta, width = 1.0)

plot(hcat(u...)',
	label = "u (ref)", color = :black, width = 2.0)
plot!(hcat(u_hist_mpc...)',
	label = "u (mpc)", color = :magenta, width = 1.0)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model, x_hist_mpc, Δt = h)


text:
using IntroAlgoCh2
using Plots, BenchmarkTools, Random

amounts = 1:300

insertionsort_times = zeros(length(amounts))
mergesort_times = zeros(length(amounts))
mergesort_optimized_times = zeros(length(amounts))
sort_times = zeros(length(amounts))

for N=amounts
    GC.gc()
    l = shuffle!(collect(1:N))

    a = copy(l)
    ti = @elapsed insertionsort!(a);

    GC.gc()
    b = copy(l)
    tm = @elapsed mergesort_iterative!(b);

    GC.gc()
    c = copy(l)
    tmo = @elapsed mergesort_iterative_optimized!(c);

    GC.gc()
    d = copy(l)
    ts = @elapsed sort!(d);

    insertionsort_times[N] = ti#mean(ti).time
    mergesort_times[N] = tm#mean(tm).time
    mergesort_optimized_times[N] = tmo
    sort_times[N] = ts
end

plot(amounts[10:end], insertionsort_times[10:end], label="insert");
plot!(amounts[10:end], mergesort_times[10:end], label="merge")
plot!(amounts[10:end], mergesort_optimized_times[10:end], label="merge_optimized")
plot!(amounts[10:end], sort_times[10:end], label="std sort")

if false
s = 1:100_000 |> collect |> shuffle!
begin
    GC.gc()
    t = copy(s)
    @btime mergesort_iterative!($t)
end
begin
    GC.gc()
    t = copy(s)
    @btime mergesort_iterative_optimized!($t)
end
begin
    GC.gc()
    t = copy(s)
    @btime sort!($t)
end
end;


text:
# General operator stuff

get_elementgrid(ao::AbstractOperator) = ao.elementgrid

Base.size(ao::AbstractOperator) = size(get_elementgrid(ao))

function Base.show(io::IO, ao::AbstractOperator)
    s = size(ao)
    print(io, "Operator $(s[end])^3 elements, $(s[1])^3 gauspoints per element")
end

Base.length(::AbstractOperator{N}) where N = N
Base.firstindex(::AbstractOperator) = 1
Base.lastindex(ao::AbstractOperator) = length(ao)

function Base.getindex(so::AbstractOperator{1}, i::Int)
    @assert i == 1
    return so
end

function Base.iterate(v::AbstractOperator, i=0)
    if i < length(v)
        return (v[i+1], i+1)
    else
        return nothing
    end
end

Unitful.unit(ao::AbstractOperator) = ao.unit
Unitful.dimension(ao::AbstractOperator) = dimension(unit(ao))

function Unitful.uconvert(u::Unitful.Units, op::AbstractOperator)
    @assert dimension(u) == dimension(op)
    unit(op) == u && return op
    conv = ustrip(uconvert(u, 1*unit(op))) * u / unit(op)
    return conv*op
end

Base.:(+)(ao::AbstractOperator) = ao
Base.:(-)(ao::AbstractOperator) = -1*ao

dot(v1::AbstractOperator{N}, v2::AbstractOperator{N}) where N = sum(v1.*v2)

function cross(v1::AbstractOperator{3}, v2::AbstractOperator{3})
    r1 = v1[2]*v2[3] - v1[3]*v2[2]
    r2 = v1[1]*v2[3] - v1[3]*v2[1]
    r3 = v1[1]*v2[2] - v1[2]*v2[1]
    return VectorOperator(r1,r2,r3)
end

Base.:(*)(op::AbstractOperator{1}, qs::QuantumState) = op(qs)
Base.:(*)(op::AbstractOperator, qs::QuantumState) = pmap(f->f(qs), op)

## Scalar Operator

"""
    ScalarOperator{T}

Type for operators that are simple multiplications. This means that the
operator is simply an array.

# Fields
- `elementgrid::CubicElementGrid`  :  grid information
- `vals::T`                        :  operator values
- `unit::Unitful.FreeUnits`        :  unit for operator

# Examples
```jldoctest
julia> ceg = CubicElementGrid(5,2,32)
Cubic elements grid with 2^3 elements with 32^3 Gauss points

julia> op = ScalarOperator(ceg, ones(size(ceg)))
Operator 2^3 elements, 32^3 gauspoints per element

julia> sin_op = sin(op)
Operator 2^3 elements, 32^3 gauspoints per element

julia> sin_op + 3op
Operator 2^3 elements, 32^3 gauspoints per element

julia> ψ = QuantumState(ceg, ones(size(ceg)))
Quantum state

julia> -2op * ψ
Quantum state
```
"""
struct ScalarOperator{T} <: AbstractScalarOperator
    elementgrid::CubicElementGrid
    vals::T
    unit::Unitful.FreeUnits
    function ScalarOperator(ceg::CubicElementGrid, ψ::AbstractArray{<:Any,6}; unit::Unitful.FreeUnits=NoUnits)
        new{typeof(ψ)}(ceg, ψ, unit)
    end
end

get_values(so::ScalarOperator) = so.vals

(so::ScalarOperator)(ψ::AbstractArray{<:Any,6}) = QuantumState(so.elementgrid, so.vals.*ψ, unit(so))
function (so::ScalarOperator)(qs::QuantumState)
    @assert size(so) == size(qs)
    return QuantumState(so.elementgrid, so.vals.*qs.psi, unit(so)*unit(qs))
end

for OP in (:(+), :(-))
    @eval begin
        function Base.$OP(a::ScalarOperator,b::ScalarOperator)
            @assert dimension(a) == dimension(b) "Dimension missmatch"
            c = uconvert(unit(a), b)
            ScalarOperator(a.elementgrid, $OP.(a.vals,c.vals); unit=a.unit)
        end
        function Base.$OP(a::ScalarOperator, x::Number)
            @assert dimension(a) == dimension(x) "Dimension missmatch"
            y = ustrip(uconvert(unit(a), x))
            return ScalarOperator(a.elementgrid, $OP.(a.vals,y); unit=a.unit)
        end
        function Base.$OP(x::Number, a::ScalarOperator)
            @assert dimension(a) == dimension(x) "Dimension missmatch"
            y = ustrip(uconvert(unit(a), x))
            return ScalarOperator(a.elementgrid, $OP.(y,a.vals); unit=a.unit)
        end
    end
end

for OP in (:(/), :(*))
    @eval begin
        function Base.$OP(a::ScalarOperator,b::ScalarOperator)
            ScalarOperator(a.elementgrid, $OP.(a.vals,b.vals); unit=$OP(unit(a), unit(b)))
        end
        function Base.$OP(a::ScalarOperator, x::Number)
            return ScalarOperator(a.elementgrid, $OP.(a.vals,ustrip(x)); unit=$OP(unit(a), unit(x)))
        end
        function Base.$OP(x::Number, a::ScalarOperator)
            return ScalarOperator(a.elementgrid, $OP.(ustrip(x),a.vals); unit=$OP(unit(x), unit(a)))
        end
    end
end

Base.:(-)(a::ScalarOperator) = -1*a
Base.:(^)(a::ScalarOperator, x::Number) = ScalarOperator(a.elementgrid, (^).(a.vals,x); unit=unit(a)^x)

for op in (:sqrt, :cbrt, :real, :imag)
    @eval Base.$op(so::ScalarOperator) = ScalarOperator(so.elementgrid, $op.(so.vals); unit=$op(so.unit))
end

for op in (:sin, :cos, :tan, :exp, :log)
    @eval function Base.$op(so::ScalarOperator)
        @assert dimension(so) == dimension(NoUnits) "Operator needs to be dimensionless"
        so1 = auconvert(so)
        return ScalarOperator(so1.elementgrid, $op.(so1.vals))
     end
end

function Unitful.uconvert(a::Unitful.Units, so::ScalarOperator)
    @assert dimension(a) == dimension(so) "Dimension missmatch"
    a == unit(so) && return so
    conv = ustrip(uconvert(a, 1*unit(so)))
    return ScalarOperator(so.elementgrid, so.vals.*conv; unit=a)
end

## Vector Operator

"""
    VectorOperator{N}

Operetor that can be considered to be a vector. This type is build out of several
scalar operators.

# Fields
- `elementgrid::CubicElementGrid`  : grid information
- `operators::Vector{AbstractOperator{1}}`    : scalar operator from which vector
    is build

# Construction
    VectorOperator(op::AbstractOperator{1}...)
"""
struct VectorOperator{N} <: AbstractOperator{N}
    elementgrid::CubicElementGrid
    operators::Vector{AbstractOperator{1}}
    function VectorOperator(op::AbstractOperator{1}...)
        si = size.(op)
        @assert all(map(x-> x==si[begin], si)) "Scalar operators have different sizes"
        new{length(op)}(op[1].elementgrid, [op...])
    end
    function VectorOperator(op::AbstractVector{<:AbstractOperator{1}})
        si = size.(op)
        @assert all(map(x-> x==si[begin], si)) "Scalar operators have different sizes"
        new{length(op)}(op[1].elementgrid, op)
    end
end

Base.firstindex(::VectorOperator) = 1
Base.getindex(v::VectorOperator, i...) = v.operators[i...]
Base.lastindex(v::VectorOperator) = length(v)

Unitful.unit(v::VectorOperator) = unit(v[begin])

for OP in (:(+), :(-))
    @eval begin
        function Base.$OP(a::VectorOperator{N},b::VectorOperator{N}) where N
            VectorOperator($OP(a.operators,b.operators))
        end
        function Base.$OP(a::VectorOperator{N},b::AbstractVector) where N
            @assert length(b) == N
            VectorOperator($OP(a.operators,b))
        end
        function Base.$OP(a::AbstractVector,b::VectorOperator{N}) where N
            @assert length(a) == N
            VectorOperator($OP(a,b.operators))
        end
    end
end

for OP in (:(+), :(-), :(*))
    @eval begin
        Base.$OP(a::VectorOperator, x::Number) = VectorOperator($OP.(a.operators,x))
        Base.$OP(x::Number, a::VectorOperator) = VectorOperator($OP.(x,a.operators))

        function Base.$OP(v::VectorOperator,s::ScalarOperator)
            @assert size(v) == size(s)
            VectorOperator( [$OP(x,s) for x in v.operators] )
        end
        function Base.$OP(s::ScalarOperator,v::VectorOperator)
            @assert size(v) == size(s)
            VectorOperator( [$OP(s,x) for x in v.operators] )
        end
    end
end

for OP in [:(/)]
    @eval begin
        function Base.$OP(v::VectorOperator,s::ScalarOperator)
            @assert size(v) == size(s)
            VectorOperator( [$OP(x,s) for x in v.operators] )
        end
        Base.$OP(a::VectorOperator, x::Number) = VectorOperator($OP.(a.operators,x))
    end
end

Base.:(-)(a::VectorOperator) = VectorOperator( [-x for x in a]  )

## Position operator

"""
    position_operator(ceg::CubicElementGrid)

Returns position operator as a [`VectorOperator`](@ref).
"""
function position_operator(ceg::CubicElementGrid)
    x(v) = v[1]
    y(v) = v[2]
    z(v) = v[3]
    return VectorOperator(
                          ScalarOperator(ceg, x.(ceg); unit=u"bohr"),
                          ScalarOperator(ceg, y.(ceg); unit=u"bohr"),
                          ScalarOperator(ceg, z.(ceg); unit=u"bohr")
                         )
end

## Derivative Operator

"""
    DerivativeOperator{NG,N}

Calculates derivative when operating. `N` is coordinate that is derived.

# Fields
- `elementgrid::CubicElementGrid`   :  grid information
- `dt::DerivativeTensor{NG}`        :  derivative in tensor form

# Creation
    DerivativeOperator(ceg::CubicElementGrid, coordinatate::Int=1)
    DerivativeOperator(ceg::CubicElementGrid, dt::DerivativeTensor, coordinatate::Int=1)
"""
struct DerivativeOperator{NG,N} <: AbstractOperator{1}
    elementgrid::CubicElementGrid
    dt::DerivativeTensor{NG}
    function DerivativeOperator(ceg::CubicElementGrid, coordinatate::Int=1)
        @assert coordinatate ∈ 1:3 "coordinate needs to be 1, 2 or 3"
        dt = DerivativeTensor(ceg)
        new{size(ceg)[1], coordinatate}(ceg, dt)
    end
    function DerivativeOperator(ceg::CubicElementGrid, dt::DerivativeTensor, coordinatate::Int=1)
        @assert coordinatate ∈ 1:3 "coordinate needs to be 1, 2 or 3"
        @assert size(ceg)[1] == size(dt)[1]
        new{size(ceg)[1], coordinatate}(ceg, dt)
    end
end

function (d::DerivativeOperator{<:Any,1})(qs::QuantumState)
    return QuantumState(d.elementgrid, operate_x(d.dt, qs.psi), unit(qs)*unit(d) )
end

function (d::DerivativeOperator{<:Any,2})(qs::QuantumState)
    return QuantumState(d.elementgrid, operate_y(d.dt, qs.psi), unit(qs)*unit(d) )
end

function (d::DerivativeOperator{<:Any,3})(qs::QuantumState)
    return QuantumState(d.elementgrid, operate_z(d.dt, qs.psi), unit(qs)*unit(d) )
end

Unitful.unit(::DerivativeOperator) = u"bohr"^-1

## Gradient Operator

"""
    GradientOperator{NG}

Calculates gradient when operating.

# Fields
- `elementgrid::CubicElementGrid`   :  grid information
- `dt::DerivativeTensor{NG}`        :  derivative in tensor form

# Creation
    GradientOperator(ceg::CubicElementGrid)
"""
struct GradientOperator{NG} <: AbstractOperator{3}
    elementgrid::CubicElementGrid
    dt::DerivativeTensor{NG}
    function GradientOperator(ceg::CubicElementGrid)
        dt = DerivativeTensor(ceg)
        new{size(ceg)[1]}(ceg, dt)
    end
end

function (go::GradientOperator)(ψ::QuantumState)
    return pmap(f->f(ψ), go)
end

function (go::GradientOperator)(a::AbstractVector{<:AbstractQuantumState})
    @assert length(a) == 3
    @assert all(map(x->size(x)==size(go), a))
    out = operate_x(go.dt, a[1].ψ)
    out .+= operate_y(go.dt, a[2].ψ)
    out .+= operate_z(go.dt, a[3].ψ)
    return QuantumState(go.ceg, out, unit(a)*unit(g))
end

#(go::GradientOperator)(ψ::QuantumState) = go(ψ.ψ)
Base.length(::GradientOperator) = 3
Base.show(io::IO, go::GradientOperator) = print(io::IO, "Gradient operator size=$(size(go))")

function Base.getindex(go::GradientOperator, i::Int)
    @assert 1<=i<=3
    return DerivativeOperator(go.elementgrid, go.dt, i)
end

Unitful.unit(::GradientOperator) = u"bohr"^-1

## derivative operations
# these are low level stuff to use only in special cases

function operate_x!(dψ::AbstractArray{<:Any,6}, dt::DerivativeTensor, ψ::AbstractArray{<:Any,6})
    tmp = similar(ψ, size(dt)...)
    tmp .= dt
    @tensor dψ[i,j,k,I,J,K] = tmp[i,l] * ψ[l,j,k,I,J,K]
    return dψ
end

function operate_y!(dψ::AbstractArray{<:Any,6}, dt::DerivativeTensor, ψ::AbstractArray{<:Any,6})
    tmp = similar(ψ, size(dt)...)
    tmp .= dt
    @tensor dψ[i,j,k,I,J,K] = tmp[j,l] * ψ[i,l,k,I,J,K]
    return dψ
end

function operate_z!(dψ::AbstractArray{<:Any,6}, dt::DerivativeTensor, ψ::AbstractArray{<:Any,6})
    tmp = similar(ψ, size(dt)...)
    tmp .= dt
    @tensor dψ[i,j,k,I,J,K] = tmp[k,l] * ψ[i,j,l,I,J,K]
    return dψ
end

function operate_x(dt::DerivativeTensor, ψ::AbstractArray{<:Any,6})
    tmp = similar(ψ)
    return operate_x!(tmp, dt, ψ)
end

function operate_y(dt::DerivativeTensor, ψ::AbstractArray{<:Any,6})
    tmp = similar(ψ)
    return operate_y!(tmp, dt, ψ)
end

function operate_z(dt::DerivativeTensor, ψ::AbstractArray{<:Any,6})
    tmp = similar(ψ)
    return operate_z!(tmp, dt, ψ)
end

## Constant operator

"""
    ConstantTimesOperator{TO,T,N}

Type to represent constant times operator. It is recommended not to create this
directly. This type is created when there is no simple way to implement operator
times constant operation.

# Fields
- `elementgrid::CubicElementGrid`  : grid information
- `op::TO`                         : operator
- `a::T`                           : constant
"""
struct ConstantTimesOperator{TO,T,N} <: AbstractOperator{N}
    elementgrid::CubicElementGrid
    op::TO
    a::T
    function ConstantTimesOperator(op::AbstractOperator{N}, a=1) where N
        if ustrip(imag(a)) ≈ 0
            new{typeof(op),typeof(real(a)), N}(op.elementgrid, op, real(a))
        else
            new{typeof(op),typeof(a), N}(op.elementgrid, op, a)
        end
    end
end

Base.:(*)(op::AbstractOperator, a::Number) = ConstantTimesOperator(op,a)
Base.:(*)(a::Number, op::AbstractOperator) = ConstantTimesOperator(op,a)
Base.:(/)(op::AbstractOperator, a::Number) = ConstantTimesOperator(op,1/a)

Base.getindex(cto::ConstantTimesOperator,i::Int) = cto.a*cto.op[i]

Unitful.unit(op::ConstantTimesOperator) = unit(op.op) * unit(op.a)

function (cto::ConstantTimesOperator{TO,T,1})(qs::QuantumState) where {TO,T}
    return cto.a * cto.op(qs)
end

function dot(cto1::ConstantTimesOperator{TO,T,N}, cto2::ConstantTimesOperator{TO,T,N}) where {TO,T,N}
    return ConstantTimesOperator(dot(cto1.op,cto2.op), cto1.a*cto2.a)
end

function dot(cto::ConstantTimesOperator{TO,T,N}, op::AbstractOperator{N}) where {TO,T,N}
    return ConstantTimesOperator(dot(cto.op,op), cto.a)
end

function dot(op::AbstractOperator{N}, cto::ConstantTimesOperator{TO,T,N}) where {TO,T,N}
    return ConstantTimesOperator(dot(cto.op,op), cto.a)
end

## Operator sum

"""
    OperatorSum{TO1, TO2, N}

Sum of two operator that cannot directly summed. It is recommended not to create
this directly.

# Fields
- `elementgrid::CubicElementGrid`   :  grid inforamation
- `op1::TO1`                        :  operator 1
- `op2::TO2`                        :  operator 2
"""
struct OperatorSum{TO1, TO2, N} <: AbstractCompositeOperator{N}
    elementgrid::CubicElementGrid
    op1::TO1
    op2::TO2
    function OperatorSum(op1::AbstractOperator{N}, op2::AbstractOperator{N}) where N
        @assert size(op1) == size(op2)
        @assert dimension(op1) == dimension(op2) "Operators have different unit dimensions"
        if unit(op1) == unit(op2)
            new{typeof(op1),typeof(op2),N}(op1.elementgrid, op1, op2)
        else
            new{typeof(op1),typeof(op2),N}(get_elementgrid(op1), op1, uconvert(unit(op1), op2))
        end
    end
end

function Base.:(+)(op1::AbstractOperator{N}, op2::AbstractOperator{N}) where N
    return OperatorSum(op1,op2)
end

function Base.:(-)(op1::AbstractOperator{N}, op2::AbstractOperator{N}) where N
    return OperatorSum(op1,-op2)
end

function (os::OperatorSum{TO1,TO2, 1})(qs::QuantumState) where {TO1,TO2}
    return os.op1(qs) + os.op2(qs)
end

Unitful.unit(sos::OperatorSum) = unit(sos.op1)

Base.getindex(op::OperatorSum, i::Int) = op.op1[i] + op.op2[i]

## Scalar operator product

"""
    OperatorProduct{TO1, TO2}

Producto of two operators that cannot be them selves multiplied. These to operator
are operated after each other. Operator 2 is operated first and then operator 1.
It is recommended not to create this directly.

# Fields
- `elementgrid::CubicElementGrid`   :  grid information
- `op1::TO1`                        :  operator 1
- `op2::TO2`                        :  operator 2
"""
struct OperatorProduct{TO1, TO2} <: AbstractCompositeOperator{1}
    elementgrid::CubicElementGrid
    op1::TO1
    op2::TO2
    function OperatorProduct(op1::AbstractOperator{1}, op2::AbstractOperator{1})
        @assert size(op1) == size(op2)
        @assert length(op1) == length(op2)  "Operators have different length"
        new{typeof(op1),typeof(op2)}(get_elementgrid(op1), op1, op2)
    end
end

function Base.:(*)(op1::AbstractOperator{1}, op2::AbstractOperator{1})
    return OperatorProduct(op1,op2)
end

(op::OperatorProduct)(qs::QuantumState) = op.op1(op.op2(qs))

Unitful.unit(op::OperatorProduct) = unit(op.op1) * unit(op.op2)

get_elementgrid(op::OperatorProduct) = get_elementgrid(op.op1)

## Laplace Operator

"""
    LaplaceOperator

∇² operator

# Fields
- `elementgrid::CubicElementGrid`   :  grid information
- `g::GradientOperator`             :  gradient operator

# Creation
    LaplaceOperator(g::GradientOperator)
"""
struct LaplaceOperator <: AbstractOperator{1}
    elementgrid::CubicElementGrid
    g::GradientOperator
    function LaplaceOperator(g::GradientOperator)
        new(g.elementgrid, g)
    end
    function LaplaceOperator(ceg::CubicElementGrid)
        new(ceg, GradientOperator(ceg))
    end
end

Unitful.unit(lo::LaplaceOperator) = unit(lo.g)^2

Base.:(*)(g1::GradientOperator, g2::GradientOperator) = LaplaceOperator(g1)

dot(g1::GradientOperator, g2::GradientOperator) = LaplaceOperator(g1)

function (lo::LaplaceOperator)(qs::QuantumState)
    tmp = similar(qs.psi, size(lo.g.dt))
    tmp .= lo.g.dt
    @tensor w[i,j]:=tmp[i,k]*tmp[k,j]
    @tensor ϕ[i,j,k,I,J,K]:= w[i,ii]*qs.psi[ii,j,k,I,J,K] + w[j,jj]*qs.psi[i,jj,k,I,J,K] + w[k,kk]*qs.psi[i,j,kk,I,J,K]
    return QuantumState(qs.elementgrid, ϕ, unit(qs)*unit(lo.g)^2)
end

## Special operators

"""
    momentum_operator(ceg::CubicElementGrid)
    momentum_operator(qs::QuantumState)
    momentum_operator(H::HamiltonOperator)
    momentum_operator(H::HamiltonOperatorFreeParticle)
    momentum_operator(H::HamiltonOperatorMagneticField)

Return momentum operator. For Hamiltonian operator input this returns
canonical momentum.
"""
function momentum_operator(ceg::CubicElementGrid)
    c = -im * 1u"ħ_au"
    g = GradientOperator(ceg)
    return c*g
end

momentum_operator(qs::QuantumState) = momentum_operator(qs.elementgrid)

function kinetic_energy_operator(ceg::CubicElementGrid, m=1u"me_au")
    @assert dimension(m) == dimension(u"kg")
    return HamiltonOperatorFreeParticle(ceg::CubicElementGrid; m=m)
end

## Hamilton operator

"""
    HamiltonOperatorFreeParticle{T}

Free particle Hamiltonian

# Fields
- `elementgrid::CubicElementGrid`   :  grid information
- `∇²::LaplaceOperator`             :  Laplace operator
- `m::T`                            :  mass of the particle

# Construction
    HamiltonOperatorFreeParticle(ceg::CubicElementGrid;m=1u"me_au")
"""
struct HamiltonOperatorFreeParticle{T} <: AbstractHamiltonOperator
    elementgrid::CubicElementGrid
    ∇²::LaplaceOperator
    m::T
    function HamiltonOperatorFreeParticle(ceg::CubicElementGrid;m=1u"me_au")
        @assert dimension(m) == dimension(u"kg")
        new{typeof(m)}(ceg, LaplaceOperator(ceg), m)
    end
end

Unitful.unit(H::HamiltonOperatorFreeParticle) =u"hartree*bohr^2"*unit(H.∇²)

function (H::HamiltonOperatorFreeParticle)(qs::QuantumState)
    return uconvert(u"hartree*bohr^2", u"ħ"^2/(-2H.m))*(H.∇²*qs)
end

"""
    HamiltonOperator{TF,TV}

Hamiltonian without magnetic field.

# Fields
- `elementgrid::CubicElementGrid`         :  grid information
- `T::HamiltonOperatorFreeParticle{TF}`   :  kinetic energy operator
- `V::TV`                                 :  potential energy operator

# Creation
    HamiltonOperator(V::AbstractOperator{1}; m=1u"me_au")
"""
struct HamiltonOperator{TF,TV} <: AbstractHamiltonOperator
    elementgrid::CubicElementGrid
    T::HamiltonOperatorFreeParticle{TF}
    V::TV
    function HamiltonOperator(V::AbstractOperator{1}; m=1u"me_au")
        @assert dimension(V) == dimension(u"J")
        @assert dimension(m) == dimension(u"kg")
        T = HamiltonOperatorFreeParticle(V.elementgrid; m=m)
        new{typeof(T.m),typeof(V)}(V.elementgrid, T, V)
    end
end

Unitful.unit(H::HamiltonOperator) = unit(H.T)

(H::HamiltonOperator)(qs::QuantumState) = H.T*qs + H.V*qs

"""
    HamiltonOperatorMagneticField{TF,TV,TA,Tq}

Hamiltonian with magnetic field.

# Fields
- `elementgrid::CubicElementGrid`         :  grid information
- `T::HamiltonOperatorFreeParticle{TF}`   :  kinetic energy operator
- `V::TV`                                 :  potential enegy operator
- `A::TA`                                 :  vector potential
- `q::Tq`                                 :  electric charge

# Construction
    HamiltonOperatorMagneticField(Args...; Kwargs...)

## Arguments for construction
- `V::AbstractOperator{1}`    :  potential energy
- `A::AbstractOperator{3}`    :  vector potential

## Keywords for construction
- `m=1u"me_au"`               :  mass of the particle
- `q=-1u"e_au"`               :  charge for the particle
"""
struct HamiltonOperatorMagneticField{TF,TV,TA,Tq} <: AbstractHamiltonOperator
    elementgrid::CubicElementGrid
    T::HamiltonOperatorFreeParticle{TF}
    V::TV
    A::TA
    q::Tq
    function HamiltonOperatorMagneticField(
            V::AbstractOperator{1},
            A::AbstractOperator{3};
            m=1u"me_au",
            q=-1u"e_au"
        )
        @assert dimension(V) == dimension(u"J")
        @assert dimension(A) == dimension(u"T*m")
        @assert dimension(m) == dimension(u"kg")
        @assert dimension(q) == dimension(u"e_au")
        T = HamiltonOperatorFreeParticle(V.elementgrid; m=m)
        new{typeof(T.m),typeof(V),typeof(A),typeof(q)}(V.elementgrid, T, V, A, q)
    end
end

Unitful.unit(H::HamiltonOperatorMagneticField) = unit(H.T)

function (H::HamiltonOperatorMagneticField)(ψ::QuantumState)
    #TODO make this better
    p = momentum_operator(H.T)
    ϕ = H.q^2*(H.A⋅H.A)*ψ
    ϕ += H.q*(p⋅H.A)*ψ
    ϕ += H.q*(H.A⋅p)*ψ
    ϕ /= 2H.T.m
    ϕ += H.V*ψ
    return H.T*ψ + ϕ
end

"""
    vector_potential(ceg, Bx, By, Bz)

Gives vector potential for constant magnetic field.
"""
function vector_potential(ceg, Bx, By, Bz)<
    @assert dimension(Bx) == dimension(By) == dimension(Bz) == dimension(u"T")
    r = position_operator(ceg)
    # A = r×B/2
    Ax = 0.5*(r[2]*Bz-r[3]*By)
    Ay = 0.5*(r[3]*Bx-r[1]*Bz)
    Az = 0.5*(r[1]*By-r[2]*Bx)
    return VectorOperator(Ax, Ay, Az)
end

function momentum_operator(H::HamiltonOperatorMagneticField)
    p = momentum_operator(H.elementgrid)
    return p + H.q*H.A
end

momentum_operator(H::HamiltonOperator) = momentum_operator(H.elementgrid)
momentum_operator(H::HamiltonOperatorFreeParticle) = momentum_operator(H.elementgrid)

## Projection operator

struct ProjectionOperator <: AbstractOperator{1}
    state::QuantumState
end

function (po::ProjectionOperator)(qs::QuantumState)
    s = bracket(po.state,qs)
    return s * po.state
end

get_elementgrid(po::ProjectionOperator) = get_elementgrid(po.state)
Unitful.unit(po::ProjectionOperator) = unit(po.state)

## Density operator

density_operator(qs::QuantumState) = ScalarOperator(qs.elementgrid, ketbra(qs, qs))
density_operator(sd::SlaterDeterminant) = sum( x -> 2*density_operator(x), sd.orbitals )

function density_operator(sd::SlaterDeterminant, occupations::AbstractVector)
    @assert length(sd.orbitals) == length(occupations)
    sum( x -> x[2]*density_operator(x[1]), zip(sd.orbitals, occupations) )
end


text:
# ---
# title: 491. Increasing Subsequences
# id: problem491
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search
# link: <https://leetcode.com/problems/increasing-subsequences/description/>
# hidden: true
# ---
#
# Given an integer array, your task is to find all the different possible
# increasing subsequences of the given array, and the length of an increasing
# subsequence should be at least 2.
#
#
#
# **Example:**
#
#
#
#     Input: [4, 6, 7, 7]
#     Output: [[4, 6], [4, 7], [4, 6, 7], [4, 6, 7, 7], [6, 7], [6, 7, 7], [7,7], [4,7,7]]
#
#
#
#
# **Constraints:**
#
#   * The length of the given array will not exceed 15.
#   * The range of integer in the given array is [-100,100].
#   * The given array may contain duplicates, and two equal integers should also be considered as a special case of increasing sequence.
#
#
## @lc code=start
using LeetCode

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


text:
#---------------------------------------------------------------------#
#This function plots the solution
#Written by F.X. Giraldo on April 19 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function plot_solution(q,coord,space_method,time)

    #Plot Solution
    text=string(space_method,", time= ","$time")
    p1=plot(coord,q[1,:],xlabel="x",ylabel="h(x,t)",legend=true,lw=3,label=[space_method],title=text,ylims=(-0.5,0.5))
    p2=plot(coord,q[2,:],xlabel="x",ylabel="U(x,t)",legend=true,lw=3,label=[space_method],title=text,ylims=(-0.5,0.5))
    display(plot(p1,p2))
    #display(p1)

end


text:
export ⊙

@concrete terse struct PointwiseProductMeasure{T} <: AbstractMeasure
    data::T
end

Base.size(μ::PointwiseProductMeasure) = size(μ.data)

function Base.show(io::IO, μ::PointwiseProductMeasure)
    io = IOContext(io, :compact => true)
    print(io, join(string.(μ.data), " ⊙ "))
end

function Base.show_unquoted(io::IO, μ::PointwiseProductMeasure, indent::Int, prec::Int)
    io = IOContext(io, :compact => true)
    if Base.operator_precedence(:*) ≤ prec
        print(io, "(")
        show(io, μ)
        print(io, ")")
    else
        show(io, μ)
    end
    return nothing
end

Base.length(m::PointwiseProductMeasure{T}) where {T} = length(m.data)

⊙(args...) = pointwiseproduct(args...)

@inline function logdensity(d::PointwiseProductMeasure, x)
    sum((logdensity(dⱼ, x) for dⱼ in d.data))
end

function sampletype(d::PointwiseProductMeasure)
    @inbounds sampletype(first(d.data))
end

basemeasure(d::PointwiseProductMeasure) = @inbounds basemeasure(first(d.data))


text:
export AbstractInputSampler, ContinousInputSampler, ConstantInputSampler, resample!, reset!

"""
    abstract type AbstractInputSampler end

Abstract type for all `PBTSamplers`. Each concret type musst implement

    `Base.getindex(sampler::ConcreteSamplerType, i)`

to get the i-th input sample.
"""
abstract type AbstractInputSampler end

"""
    resample!(s::AbstractInputSampler)

Clear the cache of the sampler. This will result in new samples.
"""
resample!(s::AbstractInputSampler) = empty!(s._cache)

"""
    reset!(s::AbstractInputSampler; seed=1)

Reset the RNG of the input sampler and clear cache.
"""
function reset!(s::AbstractInputSampler; seed=1)
    resample!(s)
    s.rng = MersenneTwister(seed)
end

"""
    ContinousInputSampler{SF} <: AbstractInputSampler

Input sampler to create time continous inputs of the form `(t) -> input`.

    ContinousInputSampler(f, t_span; seed=1)

Create an input sampler based on sample function f.
- `f`: function with signature `(rng, tspan, i)` which returns a anonymous function
  `(t) -> input` which represents the input.
- `tspan`: can be used inside `f` to determine the input
- `seed`: optional seed the RNG
"""
mutable struct ContinousInputSampler{SF} <: AbstractInputSampler
    rng::MersenneTwister
    samplef::SF
    tspan::Tuple{Float64, Float64}
    _cache::Vector{Any}

    ContinousInputSampler(f, t_span; seed=1) = new{typeof(f)}(MersenneTwister(seed), f, t_span, Any[])
end

function Base.getindex(s::ContinousInputSampler, idx::Int)
    clength = lastindex(s._cache)
    if idx > clength
        resize!(s._cache, idx)
        for i in clength+1:idx
            s._cache[i] = s.samplef(s.rng, s.tspan, i)
        end
    end
    return s._cache[idx]
end

"""
    ConstantInputSampler{SF} <: AbstractInputSampler

Input sampler to create constant inputs.

    ConstantInputSampler(f; seed=1)

Create an input sampler based on sample function f.
- `f`: function with signature `(rng, i)` which returns the input
- `seed`: optional seed the RNG
"""
mutable struct ConstantInputSampler{SF} <: AbstractInputSampler
    rng::MersenneTwister
    samplef::SF
    _cache::Vector{Any}

    ConstantInputSampler(f; seed=1) = new{typeof(f)}(MersenneTwister(seed), f, Any[])
end

function Base.getindex(s::ConstantInputSampler, idx::Int)
    clength = lastindex(s._cache)
    if idx > clength
        resize!(s._cache, idx)
        for i in clength+1:idx
            s._cache[i] = s.samplef(s.rng, i)
        end
    end
    return s._cache[idx]
end

Base.getindex(s::AbstractInputSampler, idxs::UnitRange{Int}) = [s[i] for i in idxs]


text:
export Result,
    nresults,
    nsolutions,
    nfinite,
    nsingular,
    nat_infinity,
    nfailed,
    nnonsingular,
    nreal,
    ntracked,
    finite,
    results,
    mapresults,
    failed,
    at_infinity,
    singular,
    nonsingular,
    seed,
    solutions,
    real_solutions,
    multiplicities,
    statistics,
    multiplicities!

######################
## MultiplicityInfo ##
######################

"""
    MultiplicityInfo

This contains informations about the multiplicities of the solutions.
"""
struct MultiplicityInfo
    multiplicities::Dict{Int,Vector{Vector{Int}}}
    # This set stores all path numbers which we want to ignore
    # if we don't show multiple results
    multiple_indicator::Set{Int32}
end

function MultiplicityInfo(pathresults::Vector{<:PathResult}; tol = 1e-6)
    multiple_indicator = Set{Int32}()
    multiplicities = compute_multiplicities(pathresults; tol = float(tol))
    for clusters in values(multiplicities), cluster in clusters
        for i = 2:length(cluster)
            push!(multiple_indicator, path_number(pathresults[cluster[i]]))
        end
    end
    MultiplicityInfo(multiplicities, multiple_indicator)
end

function compute_multiplicities(result::Vector{<:PathResult}; kwargs...) where {T}
    D = Dict{Int,Vector{Vector{Int}}}()
    compute_multiplicities!(D, result; kwargs...)
end
function compute_multiplicities!(
    D::Dict,
    result::Vector{<:PathResult};
    tol::Float64 = 1e-6,
) where {T}
    for m in multiplicities(solution, result, tol = tol)
        if haskey(D, length(m))
            push!(D[length(m)], m)
        else
            D[length(m)] = [m]
        end
    end
    D
end

"""
    assign_multiplicities!(path_results::Vector{<:PathResult}, I::MultiplicityInfo)

Assign the path results their multiplicity according to the multiplicity info `I`.
"""
function assign_multiplicities!(path_results::Vector{<:PathResult}, I::MultiplicityInfo)
    #assign multiplicities
    for r in path_results
        r.multiplicity[] = 1
    end
    for (k, clusters) in I.multiplicities, cluster in clusters, vᵢ in cluster
        path_results[vᵢ].multiplicity[] = k
    end
    path_results
end

is_multiple_result(r::PathResult, I::MultiplicityInfo) =
    path_number(r) ∈ I.multiple_indicator

"""
    Result{V<:AbstractVector}

The result of `solve`. This is a wrapper around the results of each single path
([`PathResult`](@ref)) and it contains some additional informations like a random seed to
replicate the result.
"""
struct Result{V}
    pathresults::Vector{PathResult{V}}
    tracked_paths::Int
    seed::Int
    retracked_paths::Int
    multiplicity_info::Base.RefValue{MultiplicityInfo}
end

function Result(
    pathresults::Vector{PathResult{V}},
    tracked_paths::Int,
    seed::Int;
    retracked_paths::Int = 0,
    multiplicity_tol::Float64 = 1e-6,
) where {V}
    multiplicity_info = MultiplicityInfo(pathresults; tol = multiplicity_tol)
    assign_multiplicities!(pathresults, multiplicity_info)
    Result(pathresults, tracked_paths, seed, retracked_paths, Ref(multiplicity_info))
end

Base.length(r::Result) = length(r.pathresults)
Base.getindex(r::Result, I) = getindex(r.pathresults, I)

Base.iterate(r::Result) = iterate(r.pathresults)
Base.iterate(r::Result, state) = iterate(r.pathresults, state)
Base.lastindex(r::Result) = lastindex(r.pathresults)
Base.eltype(r::Type{Result{V}}) where {V} = PathResult{V}
solution_type(::Result{V}) where {V} = V

is_multiple_result(r::PathResult, R::Result) = is_multiple_result(r, R.multiplicity_info[])
is_multiple_result(r::PathResult, R::Vector{<:PathResult}) = false

"""
    multiplicities!(result::Result; tol=1e-6)

Compute the multiplicities of the solutions in `result` with respect to the given tolerance.
"""
function multiplicities!(result::Result; tol = 1e-6)
    multiplicity_info = MultiplicityInfo(result.pathresults; tol = tol)
    result.multiplicity_info[] = multiplicity_info
    assign_multiplicities!(result.pathresults, multiplicity_info)
    result
end

const Results = Union{Result,Vector{<:PathResult}}
const ProjectiveResult = Result{<:PVector}

"""
    nresults(
        result;
        only_real = false,
        real_tol = 1e-6,
        only_nonsingular = false,
        singular_tol = 1e10,
        onlyfinite = true,
    )

The number of solutions which satisfy the corresponding predicates.

## Example
```julia
result = solve(F)
# Get all non-singular results where all imaginary parts are smaller than 1e-8
nresults(result, only_real=true, real_tol=1e-8, only_nonsingular=true)
```
"""
function nresults(
    R::Results;
    only_real = false,
    real_tol = 1e-6,
    only_nonsingular = false,
    only_singular = false,
    singular_tol = 1e10,
    onlyfinite = true,
    multiple_results = false,
)
    count(R) do r
        (
            !only_real || is_real(r, real_tol)
        ) && (
            !only_nonsingular || is_nonsingular(r, singular_tol)
        ) && (
            !only_singular || is_singular(r, singular_tol)
        ) && (
            !onlyfinite || isfinite(r) || is_projective(r)
        ) && (multiple_results || !is_multiple_result(r, R))
    end
end

"""
    statistics(
        R::Result;
        only_real = false,
        real_tol = 1e-6,
        only_nonsingular = false,
        only_singular = false,
        singular_tol = 1e10,
    )

Statistic about the number of (real) singular and non-singular solutions etc.
Returns a named tuple with the statistics.

## Example
```julia
julia> statistics(solve([x^2+y^2-2, 2x+3y-1]))
    (
     nonsingular = 2,
     singular = 0,
     real_nonsingular = 2,
     real_singular = 0,
     real = 2,
     at_infinity = 0,
     failed = 0,
     total = 2,
    )
"""
function statistics(
    R::Results,
    only_real = false,
    real_tol = 1e-6,
    only_nonsingular = false,
    only_singular = false,
    singular_tol = 1e10,
)

    failed = at_infinity = nonsingular = singular = real_nonsingular = real_singular = 0
    singular_with_multiplicity = real_singular_with_multiplicity = 0
    for r in R
        is_multiple_result(r, R) && continue

        if is_failed(r)
            failed += 1
        elseif is_singular(r, singular_tol)
            if is_real(r, real_tol)
                real_singular += 1
                real_singular_with_multiplicity += unpack(multiplicity(r), 1)
            end
            singular += 1
            singular_with_multiplicity += unpack(multiplicity(r), 1)
        elseif !is_projective(r) && !isfinite(r)
            at_infinity += 1
        else # finite, nonsingular
            if is_real(r, real_tol)
                real_nonsingular += 1
            end
            nonsingular += 1
        end
    end
    (
        nonsingular = nonsingular,
        singular = singular,
        singular_with_multiplicity = singular_with_multiplicity,
        real_nonsingular = real_nonsingular,
        real_singular = real_singular,
        real_singular_with_multiplicity = real_singular_with_multiplicity,
        real = real_nonsingular + real_singular,
        at_infinity = at_infinity,
        failed = failed,
        total = R.tracked_paths,
    )
end

"""
    nfinite(result)

The number of finite solutions.
"""
nfinite(R::Results) = count(isfinite, R)

"""
    nsolutions(result)

The number of solutions.
"""
nsolutions(R::Results) = nresults(R)

"""
    nsingular(
        result;
        singular_tol = 1e10,
        multiplicitytol = 1e-5,
        counting_multiplicities = false,
        kwargs...,
    )

The number of singular solutions. A solution is considered singular if its windingnumber is
larger than 1 or the condition number is larger than `tol`.
If `counting_multiplicities=true` the number of singular solutions times their
multiplicities is returned.
"""
function nsingular(
    R::Results;
    singular_tol = 1e10,
    counting_multiplicities = false,
    kwargs...,
)
    S = results(
        R;
        only_singular = true,
        multiple_results = false,
        singular_tol = singular_tol,
        kwargs...,
    )
    isempty(S) && return 0
    counting_multiplicities && return sum(multiplicity, S)
    length(S)
end

"""
    nat_infinity(result)

The number of solutions at infinity.
"""
nat_infinity(R::Results) = count(is_at_infinity, R)

"""
    nafailed(result)

The number of failed paths.
"""
nfailed(R::Results) = count(is_failed, R)

"""
    nnonsingular(result; tol=1e-10)

The number of non-singular solutions.
"""
nnonsingular(R::Result; tol = 1e10) = count(r -> is_nonsingular(r, tol), R)

"""
    nreal(result; tol=1e-6)

The number of real solutions where all imaginary parts of each solution
are smaller than `tol`.
"""
nreal(R::Results; tol = 1e-6) = count(r -> is_real(r, tol), R)

"""
    ntracked(R::Result)

Returns the total number of paths tracked.
"""
ntracked(R::Result) = R.tracked_paths

"""
    seed(result)

The random seed used in the computation.
"""
seed(result::Result) = result.seed

# Filtering
"""
    results(result; only_real=false, real_tol=1e-6, only_nonsingular=false,
                onlysigular=false, singular_tol=1e10, onlyfinite=true, multiple_results=false)

Return all `PathResult`s for which the given conditions apply.

## Example

```julia
R = solve(F)

# This gives us all PathResults considered non-singular and real (but still as a complex vector).
real_solutions = results(R, only_real=true, only_nonsingular=true)
```
"""
results(R::Results; kwargs...) = mapresults(identity, R; kwargs...)
results(f::Function, R::Results; kwargs...) = mapresults(f, R; kwargs...)

"""
    mapresults(f::Function, result; conditions...)

Apply the function `f` to all `PathResult`s for which the given conditions apply. For the
possible conditions see [`results`](@ref).

## Example
```julia
# This gives us all solutions considered real (but still as a complex vector).
real_solutions = mapresults(solution, R, only_real=true)
```
"""
function mapresults(
    f::Function,
    R::Results;
    only_real = false,
    real_tol = 1e-6,
    only_nonsingular = false,
    only_singular = false,
    singular_tol = 1e10,
    onlyfinite = true,
    multiple_results = false,
)
    [
        f(r)
        for r in R if (!only_real || is_real(
            r,
            real_tol,
        )) && (!only_nonsingular || is_nonsingular(
            r,
            singular_tol,
        )) && (!only_singular || is_singular(
            r,
            singular_tol,
        )) && (!onlyfinite || isfinite(
            r,
        ) || is_projective(r)) && (multiple_results || !is_multiple_result(r, R))
    ]
end

"""
    solutions(result; conditions...)

Return all solution (as `Vector`s) for which the given conditions apply.
For the possible `conditions` see [`results`](@ref).

## Example
```julia
julia> @polyvar x y
julia> result = solve([(x-2)y, y+x+3]);
julia> solutions(result)
[[2.0+0.0im, -5.0+0.0im], [-3.0+0.0im, 0.0+0.0im]]
```
"""
function solutions(result::Results; kwargs...)
    mapresults(solution, result; kwargs...)
end

"""
    real_solutions(result; tol=1e-6, conditions...)

Return all real solution (as `Vector`s of reals) for which the given conditions apply.
For the possible `conditions` see [`results`](@ref). Note that `only_real` is always `true`
and `real_tol` is now `tol`.

## Example
```julia
julia> @polyvar x y
julia> result = solve([(x-2)y, y+x+3]);
julia> real_solutions(result)
[[2.0, -5.0], [-3.0, 0.0]]
```
"""
function real_solutions(result::Results; only_real = true, tol = 1e-6, kwargs...)
    mapresults(real_vector ∘ solution, result; only_real = true, real_tol = tol, kwargs...)
end
@deprecate realsolutions(result; kwargs...) real_solutions(result; kwargs...)

"""
    nonsingular(result::Results; conditions...)

Return all `PathResult`s for which the solution is non-singular. This is just a shorthand
for `results(R; only_nonsingular=true, conditions...)`. For the possible `conditions` see
[`results`](@ref).
"""
nonsingular(R::Results; kwargs...) = results(R; only_nonsingular = true, kwargs...)

"""
    singular(R::Results; tol=1e10, multiple_results=false, kwargs...)

Return all `PathResult`s for which the solution is singular. A solution is labeled singular
if the condition number is greater than `singular_tol`, or if the winding number is > 1.
If `multiple_results=false` only one point from each cluster of multiple solutions is returned.
If If `multiple_results=true` all singular solutions in `R` are returned.
For the possible `kwargs` see [`results`](@ref).
"""
function singular(R::Results; tol = 1e10, kwargs...)
    results(R; only_singular = true, singular_tol = tol, kwargs...)
end

"""
    finite(result::AffineResults; conditions...)

Return all `PathResult`s for which the solution is finite. This is just a shorthand
for `results(R; onlyfinite=true, conditions...)`. For the possible `conditions` see
[`results`](@ref).
"""
finite(R::Results; kwargs...) = results(R; onlyfinite = true, kwargs...)

"""
    real(result, tol=1e-6)

Get all results where the solutions are real with the given tolerance `tol`.
See [`is_real`](@ref) for details regarding the determination of 'realness'.
"""
Base.real(R::Results; tol = 1e-6) = [r for r in R if is_real(r, tol)]

"""
    failed(result)

Get all results where the path tracking failed.
"""
failed(R::Results) = [r for r in R if is_failed(r)]

"""
    at_infinity(result::AffineResult)

Get all results where the solutions is at infinity.
"""
at_infinity(R::Results) = [r for r in R if is_at_infinity(r)]

"""
    multiplicities(V::Results; tol=1e-6)

Returns a `Vector` of `Vector{PathResult}`s grouping the `PathResult`s whose solutions
appear with multiplicities *greater* 1 in 'V'.
Two solutions are regarded as equal, when their pairwise distance is less than 'tol'.
"""
function multiplicities(results::Results; tol = 1e-6)
    map(i -> results[i], multiplicities(solution, results; tol = tol))
end

####Show function

function Base.show(io::IO, x::Result)
    s = statistics(x)
    header = "Result{$(solution_type(x))} with $(s.nonsingular + s.singular) solutions"
    println(io, header)
    println(io, "="^(length(header)))
    println(
        io,
        "• $(s.nonsingular) non-singular $(plural("solution", s.nonsingular)) ($(s.real_nonsingular) real)",
    )
    println(
        io,
        "• $(s.singular) singular $(plural("solution", s.singular)) ($(s.real_singular) real)",
    )
    s.at_infinity > 0 && println(
        io,
        "• $(s.at_infinity) $(plural("solution", s.at_infinity)) at infinity",
    )
    s.failed > 0 && println(io, "• $(s.failed) failed $(plural("path", s.failed))")
    println(io, "• $(ntracked(x)) paths tracked")
    println(io, "• random seed: $(seed(x))")
    if s.singular > 0
        println(io, "• multiplicity table of singular solutions:")
        singular_multiplicities_table(io, x, s)
    end
end

function Base.show(io::IO, x::ProjectiveResult)
    s = statistics(x)
    header = "Result{$(solution_type(x))} with $(s.nonsingular + s.singular) solutions"
    println(io, header)
    println(io, "="^(length(header)))
    println(
        io,
        "• $(s.nonsingular) non-singular $(plural("solution", s.nonsingular)) ($(s.real_nonsingular) real)",
    )
    println(
        io,
        "• $(s.singular) singular $(plural("solution", s.singular)) ($(s.real_singular) real)",
    )
    s.failed > 0 && println(io, "• $(s.failed) failed $(plural("path", s.failed))")
    println(io, "• $(ntracked(x)) paths tracked")
    println(io, "• random seed: $(seed(x))")
    if s.singular > 0
        println(io, "• multiplicity table of singular solutions:")
        singular_multiplicities_table(io, x, s)
    end
end

TreeViews.hastreeview(::Result) = true
TreeViews.hastreeview(::ProjectiveResult) = true
TreeViews.numberofnodes(::Result) = 8
TreeViews.numberofnodes(::ProjectiveResult) = 7
TreeViews.treelabel(io::IO, x::Result, ::MIME"application/prs.juno.inline") = print(
    io,
    "" *
        "Result{$(solution_type(x))}" *
        "",
)

function TreeViews.nodelabel(io::IO, x::Result, i::Int, ::MIME"application/prs.juno.inline")
    s = statistics(x)
    if i == 1
        print(io, "Paths tracked")
    elseif i == 2 && s.nonsingular > 0
        print(io, "$(s.nonsingular) non-singular ($(s.real_nonsingular) real)")
    elseif i == 3 && s.singular > 0
        print(io, "$(s.singular) singular ($(s.real_singular) real)")
    elseif i == 4 && (s.real_nonsingular + s.real_singular) > 0
        print(io, "$(s.real_nonsingular+s.real_singular) real")
    elseif i == 5 && s.at_infinity > 0
        print(io, "$(s.at_infinity) at_infinity")
    elseif i == 6 && s.failed > 0
        print(io, "$(s.failed) failed")
    elseif i == 7
        print(io, "Random seed used")
    elseif i == 8 && s.singular > 0
        print(io, "  multiplicity table of singular solutions: \n")
        singular_multiplicities_table(io, x, s)
    end
end

function TreeViews.treenode(r::Result, i::Integer)
    s = statistics(r)
    if i == 1
        return ntracked(r)
    elseif i == 2 && s.nonsingular > 0
        return finite(r, only_nonsingular = true)
    elseif i == 3 && s.singular > 0
        return finite(r, only_singular = true)
    elseif i == 4 && (s.real_nonsingular + s.real_singular) > 0
        return finite(r, only_real = true)
    elseif i == 5 && s.at_infinity > 0
        return at_infinity(r)
    elseif i == 6 && s.failed > 0
        return failed(r)
    elseif i == 7
        return seed(r)
    elseif i == 8 && s.singular > 0
        return missing
    end
    missing
end

function TreeViews.nodelabel(
    io::IO,
    x::ProjectiveResult,
    i::Int,
    ::MIME"application/prs.juno.inline",
)
    s = statistics(x)
    if i == 1
        print(io, "Paths tracked")
    elseif i == 2 && s.nonsingular > 0
        print(io, "$(s.nonsingular) non-singular ($(s.real_nonsingular) real)")
    elseif i == 3 && s.singular > 0
        print(io, "$(s.singular) singular ($(s.real_singular) real)")
    elseif i == 4 && (s.real_nonsingular + s.real_singular) > 0
        print(io, "$(s.real_nonsingular + s.real_singular) real solutions")
    elseif i == 5 && s.failed > 0
        print(io, "$(s.failed) failed")
    elseif i == 6
        print(io, "Random seed used")
    elseif i == 7 && s.singular > 0
        print(io, "  multiplicity table of singular solutions: \n")
        singular_multiplicities_table(io, x, s)
    end
end

function TreeViews.treenode(r::ProjectiveResult, i::Integer)
    s = statistics(r)
    if i == 1
        return length(r)
    elseif i == 2 && s.nonsingular > 0
        return finite(r, only_nonsingular = true)
    elseif i == 3 && s.singular > 0
        return finite(r, only_singular = true)
    elseif i == 4 && (s.real_nonsingular + s.real_singular) > 0
        return finite(r, only_real = true)
    elseif i == 5 && s.failed > 0
        return failed(r)
    elseif i == 6
        return seed(r)
    elseif i == 7 && s.singular > 0
        return missing
    end
    missing
end

plural(singularstr, n) = n == 1 ? singularstr : singularstr * "s"

function singular_multiplicities_table(io, result::Result, stats = statistics(result))
    M = result.multiplicity_info[].multiplicities
    if isempty(M)
        n_higher_mults_total = 0
    else
        n_higher_mults_total = sum(length, values(M))
    end

    headers = ["mult.", "total", "# real", "# non-real"]
    mult_1_exists = n_higher_mults_total < stats.singular
    data = Matrix{Int}(undef, mult_1_exists + length(M), 4)

    n_higher_real = 0
    i = mult_1_exists + 1
    for k in sort(collect(keys(M)))
        data[i, 1] = k
        n_real_solsᵢ = count(Mᵢ -> is_real(result.pathresults[Mᵢ[1]]), M[k])
        n_solsᵢ = length(M[k])
        data[i, 2] = n_solsᵢ
        data[i, 3] = n_real_solsᵢ
        data[i, 4] = n_solsᵢ - n_real_solsᵢ
        n_higher_real += n_real_solsᵢ
        i += 1
    end

    if mult_1_exists
        data[1, 1] = 1
        data[1, 2] = stats.singular - n_higher_mults_total
        data[1, 3] = stats.real_singular - n_higher_real
        data[1, 4] = data[1, 2] - data[1, 3]
    end

    PrettyTables.pretty_table(
        io,
        data,
        headers,
        PrettyTables.unicode;
        alignment = :c,
        header_crayon = PrettyTables.Crayon(bold = false),
        border_crayon = PrettyTables.Crayon(faint = true),
    )
end