text:
using Test, RDKitMinimalLib

molblockv2000 = """

     RDKit          2D

 13 13  0  0  0  0  0  0  0  0999 V2000
    8.8810   -2.1206    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    8.8798   -2.9479    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    9.5946   -3.3607    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   10.3110   -2.9474    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   10.3081   -2.1170    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    9.5928   -1.7078    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   11.0210   -1.7018    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   11.7369   -2.1116    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
   11.0260   -3.3588    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
   11.0273   -4.1837    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   11.7423   -4.5949    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   10.3136   -4.5972    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
   11.0178   -0.8769    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  2  0
  5  7  1  0
  3  4  2  0
  7  8  2  0
  4  9  1  0
  4  5  1  0
  9 10  1  0
  2  3  1  0
 10 11  1  0
  5  6  2  0
 10 12  2  0
  6  1  1  0
  7 13  1  0
M  END
"""

@testset "io" begin
    mol = get_mol(molblockv2000)
    @test mol.mol_size[] == 0x000000000000042a
    qmol = get_qmol("c1ccccc1")
    @test qmol.mol_size[] == 0x00000000000001eb
    @test get_smiles(mol) == "CC(=O)Oc1ccccc1C(=O)O"
    @test get_smarts(qmol) == "c1ccccc1"
    @test get_cxsmiles(mol) == "CC(=O)Oc1ccccc1C(=O)O |(11.7423,-4.5949,;11.0273,-4.1837,;10.3136,-4.5972,;11.026,-3.3588,;10.311,-2.9474,;9.5946,-3.3607,;8.8798,-2.9479,;8.881,-2.1206,;9.5928,-1.7078,;10.3081,-2.117,;11.021,-1.7018,;11.7369,-2.1116,;11.0178,-0.8769,)|"
    @test occursin("V2000", get_molblock(mol))
    @test occursin("V3000", get_v3kmolblock(mol))
    @test occursin("commonchem", get_json(mol))
    @test get_inchi(mol) == "InChI=1S/C9H8O4/c1-6(10)13-8-5-3-2-4-7(8)9(11)12/h2-5H,1H3,(H,11,12)"
    @test get_inchi_for_molblock(molblockv2000) == "InChI=1S/C9H8O4/c1-6(10)13-8-5-3-2-4-7(8)9(11)12/h2-5H,1H3,(H,11,12)"
    @test get_inchikey_for_inchi("InChI=1S/C9H8O4/c1-6(10)13-8-5-3-2-4-7(8)9(11)12/h2-5H,1H3,(H,11,12)") == "BSYNRYMUTXBXSQ-UHFFFAOYSA-N"
end

@testset "drawing" begin
    mol = get_mol(molblockv2000)
    @test occursin("width='350px'", get_svg(mol, Dict{String,Any}("height" => 300, "width" => 350)))

    qmol = get_qmol("c1ccccc1")
    smatch = get_substruct_match(mol, qmol)
    @test occursin("fill:#FF7F7F", get_svg(mol, smatch))
end

@testset "calculators" begin
    mol = get_mol(molblockv2000)
    fp_details = Dict{String,Any}("nBits" => 64, "radius" => 2)

    @test get_morgan_fp(mol, fp_details) == "1110010110110100110000110000000001001001000000010010000101001100"
    @test count(i -> (i == '1'), get_morgan_fp(mol)) == 24
    tb = get_morgan_fp_as_bytes(mol)
    @test sum([count_ones(b) for b in tb]) == 24

    @test get_rdkit_fp(mol, fp_details) == "1111111111111111111111111111111111111111111111111111111111111111"
    @test count(i -> (i == '1'), get_rdkit_fp(mol)) == 354
    tb = get_rdkit_fp_as_bytes(mol)
    @test sum([count_ones(b) for b in tb]) == 354

    @test get_pattern_fp(mol, fp_details) == "1111111111111101111111111111111011011111111111111111111111111111"
    @test count(i -> (i == '1'), get_pattern_fp(mol)) == 173
    tb = get_pattern_fp_as_bytes(mol)
    @test sum([count_ones(b) for b in tb]) == 173

    @test get_descriptors(mol) == Dict{String, Any}("CrippenMR" => 44.7103, "lipinskiHBD" => 1.0, "kappa1" => 9.2496, "kappa3" => 2.29741, "NumHBD" => 1.0, "NumSpiroAtoms" => 0.0, "chi1n" => 3.61745, "FractionCSP3" => 0.11111, "chi0n" => 6.98135, "NumUnspecifiedAtomStereoCenters" => 0.0, "NumAmideBonds" => 0.0, "NumAromaticHeterocycles" => 0.0, "NumAtoms" => 21.0, "tpsa" => 63.6, "NumHeteroatoms" => 4.0, "Phi" => 2.63916, "exactmw" => 180.04225, "hallKierAlpha" => -1.83999, "CrippenClogP" => 1.31009, "chi3v" => 1.37115, "NumSaturatedHeterocycles" => 0.0, "NumAromaticRings" => 1.0, "labuteASA" => 74.75705, "amw" => 180.15899, "NumAliphaticHeterocycles" => 0.0, "NumAtomStereoCenters" => 0.0, "chi4v" => 0.88717, "NumRotatableBonds" => 2.0, "lipinskiHBA" => 4.0, "NumRings" => 1.0, "chi0v" => 6.98135, "NumHeterocycles" => 0.0, "NumHeavyAtoms" => 13.0, "chi3n" => 1.37115, "chi1v" => 3.61745, "NumBridgeheadAtoms" => 0.0, "kappa2" => 3.70925, "chi2n" => 1.37115, "NumHBA" => 3.0, "chi4n" => 0.88717, "NumSaturatedRings" => 0.0, "NumAliphaticRings" => 0.0, "chi2v" => 1.37115)
end

@testset "standardization" begin
    # normalize (sulfoxide)
    mol = get_mol("CS(C)=O")
    normalize(mol)
    smiles = get_smiles(mol)
    @test smiles == "C[S+](C)[O-]"

    # cleanup
    mol = get_mol("[Pt]CCN(=O)=O", Dict{String,Any}("sanitize" => false))
    smiles = get_smiles(mol)
    @test smiles == "O=N(=O)CC[Pt]"
    cleanup(mol)
    smiles = get_smiles(mol)
    @test smiles == "[CH2-]C[N+](=O)[O-].[Pt+]"

    # fragment_parent
    fragment_parent(mol)
    smiles = get_smiles(mol)
    @test smiles == "[CH2-]C[N+](=O)[O-]"

    # charge_parent
    charge_parent(mol)
    smiles = get_smiles(mol)
    @test smiles == "CC[N+](=O)[O-]"

    mol = get_mol("[Pt]CCN(=O)=O", Dict{String,Any}("sanitize" => false))
    charge_parent(mol)
    smiles = get_smiles(mol)
    @test smiles == "CC[N+](=O)[O-]"

    mol = get_mol("[Pt]CCN(=O)=O", Dict{String,Any}("sanitize" => false))
    charge_parent(mol, Dict{String,Any}("skipStandardize" => true))
    smiles = get_smiles(mol)
    @test smiles == "[CH2-]C[N+](=O)[O-].[Pt+]"

    # neutralize
    mol = get_mol("[CH2-]CN(=O)=O", Dict{String,Any}("sanitize" => false))
    neutralize(mol)
    smiles = get_smiles(mol)
    @test smiles == "CCN(=O)=O"

    # reionize
    mol = get_mol("[O-]c1cc(C(=O)O)ccc1")
    reionize(mol)
    smiles = get_smiles(mol)
    @test smiles == "O=C([O-])c1cccc(O)c1"

    # canonical_tautomer
    mol = get_mol("OC(O)C(=N)CO")
    canonical_tautomer(mol)
    smiles = get_smiles(mol)
    @test smiles == "NC(CO)C(=O)O"
end

@testset "coordinates" begin
    mol = get_mol(molblockv2000)
    set_3d_coords(mol)
    @test occursin("RDKit          3D", get_molblock(mol))
    set_2d_coords(mol)
    @test occursin("RDKit          2D", get_molblock(mol))

    mol = get_mol("CC(=O)Oc1ccccc1C(=O)O")
    set_3d_coords(mol)
    template = get_mol("CC(=O)Nc1ccc(O)cc1")
    set_2d_coords(template)
    set_2d_coords_aligned(mol, template)
    @test occursin("RDKit          2D", get_molblock(mol))
end

@testset "modification" begin
    mol = get_mol(molblockv2000)
    hcount = occursin("H", get_molblock(mol))
    add_hs(mol)
    @test occursin("H", get_molblock(mol)) > hcount

    remove_all_hs(mol)
    @test occursin("H", get_molblock(mol)) == hcount
end

@testset "substructure" begin
    mol = get_mol("CC(=O)Oc1ccccc1C(=O)O")
    qmol = get_qmol("c1ccccc1")
    @test get_substruct_match(mol, qmol) == Dict{String,Any}("bonds" => Any[4, 5, 6, 7, 8, 12], "atoms" => Any[4, 5, 6, 7, 8, 9])
    @test get_substruct_matches(mol, qmol) == Any[Dict{String,Any}("bonds" => Any[4, 5, 6, 7, 8, 12], "atoms" => Any[4, 5, 6, 7, 8, 9])]
end


text:
import Base: show

"""
    UserNSRunner

A `UserNSRunner` represents an "execution context", an object that bundles all
necessary information to run commands within the container that contains
our crossbuild environment.  Use `run()` to actually run commands within the
`UserNSRunner`, and `runshell()` as a quick way to get an interactive shell
within the crossbuild environment.
"""
mutable struct UserNSRunner <: Runner
    sandbox_cmd::Cmd
    env::Dict{String, String}
    platform::Platform
end

function platform_def_mapping(platform)
    tp = triplet(platform)
    mapping = Pair{String,String}[
        shards_dir(tp) => joinpath("/opt", tp)
    ]

    # We might also need the x86_64-linux-gnu platform for bootstrapping,
    # so make sure that's always included
    if platform != Linux(:x86_64)
        ltp = triplet(Linux(:x86_64))
        push!(mapping, shards_dir(ltp) => joinpath("/opt", ltp))
    end

    # If we're trying to run macOS and we have an SDK directory, mount that!
    if platform == MacOS()
        sdk_version = "MacOSX10.10.sdk"
        sdk_shard_path = shards_dir(sdk_version)
        push!(mapping, sdk_shard_path => joinpath("/opt", tp, sdk_version))
    end

    # Reverse mapping order, because `sandbox` reads them backwards
    reverse!(mapping)
    return mapping
end

function UserNSRunner(workspace_root::String; cwd = nothing,
                      platform::Platform = platform_key(),
                      extra_env=Dict{String, String}(),
                      verbose::Bool = false,
                      mappings::Vector = platform_def_mapping(platform),
                      workspaces::Vector = Pair[])
    global use_ccache

    # Ensure the rootfs for this platform is downloaded and up to date.
    # Also, since we require the Linux(:x86_64) shard for HOST_CC....
    update_rootfs(triplet.([platform, Linux(:x86_64)]); verbose=verbose)

    # Check to make sure we're not going to try and bindmount within an
    # encrypted directory, as that triggers kernel bugs
    check_encryption(workspace_root; verbose=verbose)

    # Construct environment variables we'll use from here on out
    envs = merge(target_envs(triplet(platform)), extra_env)

    # the workspace_root is always a workspace, and we always mount it first
    insert!(workspaces, 1, workspace_root => "/workspace")

    # If we're enabling ccache, then mount in a read-writeable volume at /root/.ccache
    if use_ccache
        if !isdir(ccache_dir())
            mkpath(ccache_dir())
        end
        push!(workspaces, ccache_dir() => "/root/.ccache")
    end

    # Construct sandbox command
    sandbox_cmd = `$(sandbox_path())`
    if verbose
        sandbox_cmd = `$sandbox_cmd --verbose`
    end
    sandbox_cmd = `$sandbox_cmd --rootfs $(rootfs_dir())`
    if cwd != nothing
        sandbox_cmd = `$sandbox_cmd --cd $cwd`
    end

    # Add in read-only mappings and read-write workspaces
    for (outside, inside) in mappings
        sandbox_cmd = `$sandbox_cmd --map $outside:$inside`
    end
    for (outside, inside) in workspaces
        sandbox_cmd = `$sandbox_cmd --workspace $outside:$inside`
    end

    # Check to see if we need to run privileged containers.
    if runner_override == "privileged"
        # Next, prefer `sudo`, but allow fallback to `su`. Also, force-set
        # $LD_LIBRARY_PATH with these commands, because it is typically
        # lost and forgotten.  :(
        if sudo_cmd() == `sudo`
            sandbox_cmd = `$(sudo_cmd()) -E LD_LIBRARY_PATH=$(envs["LD_LIBRARY_PATH"]) $sandbox_cmd`
        else
            sandbox_cmd = `$(sudo_cmd()) "$sandbox_cmd"`
        end
    end

    # Finally, return the UserNSRunner in all its glory
    return UserNSRunner(sandbox_cmd, envs, platform)
end

function show(io::IO, x::UserNSRunner)
    p = x.platform
    # Displays as, e.g., Linux x86_64 (glibc) UserNSRunner
    write(io, "$(typeof(p).name.name)", " ", arch(p), " ",
          Compat.Sys.islinux(p) ? "($(p.libc)) " : "",
          "UserNSRunner")
end

function Base.run(ur::UserNSRunner, cmd, logpath::AbstractString; verbose::Bool = false, tee_stream=stdout)
    if runner_override == "privileged"
        msg = "Running privileged container via `sudo`, may ask for your password:"
        BinaryProvider.info_onchange(msg, "privileged", "userns_run_privileged")
    end

    did_succeed = true
    cd(rootfs_dir()) do
        oc = OutputCollector(setenv(`$(ur.sandbox_cmd) -- $(cmd)`, ur.env); verbose=verbose, tee_stream=tee_stream)

        did_succeed = wait(oc)

        if !isempty(logpath)
            # Write out the logfile, regardless of whether it was successful
            mkpath(dirname(logpath))
            open(logpath, "w") do f
                # First write out the actual command, then the command output
                println(f, cmd)
                print(f, merge(oc))
            end
        end
    end

    # Return whether we succeeded or not
    return did_succeed
end

const AnyRedirectable = Union{Base.AbstractCmd, Base.TTY, IOStream}
function run_interactive(ur::UserNSRunner, cmd::Cmd; stdin = nothing, stdout = nothing, stderr = nothing)
    if runner_override == "privileged"
        msg = "Running privileged container via `sudo`, may ask for your password:"
        BinaryProvider.info_onchange(msg, "privileged", "userns_run_privileged")
    end

    cd(rootfs_dir()) do
        cmd = setenv(`$(ur.sandbox_cmd) -- $(cmd)`, ur.env)
        if stdin isa AnyRedirectable
            cmd = pipeline(cmd, stdin=stdin)
        end
        if stdout isa AnyRedirectable
            cmd = pipeline(cmd, stdout=stdout)
        end
        if stderr isa AnyRedirectable
            cmd = pipeline(cmd, stderr=stderr)
        end

        if stdout isa IOBuffer
            if !(stdin isa IOBuffer)
                stdin = devnull
            end
            out, process = open(cmd, "r", stdin)
            @schedule begin
                while !eof(out)
                    write(stdout, read(out))
                end
            end
            wait(process)
        else
            run(cmd)
        end
    end
end

function runshell(ur::UserNSRunner, args...; kwargs...)
    run_interactive(ur, `/bin/bash`, args...; kwargs...)
end

function runshell(::Type{UserNSRunner}, platform::Platform = platform_key(); verbose::Bool = false)
    ur = UserNSRunner(
        pwd();
        cwd="/workspace/",
        platform=platform,
        verbose=verbose
    )
    return runshell(ur)
end

function probe_unprivileged_containers(;verbose::Bool=false)
    # Ensure the base rootfs is available
    update_rootfs(String[]; verbose=false)

    # Ensure we're not about to make fools of ourselves by trying to mount an
    # encrypted directory, which triggers kernel bugs.  :(
    check_encryption(pwd())

    # Construct an extremely simple sandbox command
    sandbox_cmd = `$(sandbox_path()) --rootfs $(rootfs_dir())`
    cmd = `$(sandbox_cmd) -- /bin/bash -c "echo hello julia"`

    if verbose
        Compat.@info("Probing for unprivileged container capability...")
    end
    oc = OutputCollector(cmd; verbose=verbose, tail_error=false)
    return wait(oc) && merge(oc) == "hello julia\n"
end

"""
    is_ecryptfs(path::AbstractString; verbose::Bool=false)

Checks to see if the given `path` (or any parent directory) is placed upon an
`ecryptfs` mount.  This is known not to work on current kernels, see this bug
for more details: https://bugzilla.kernel.org/show_bug.cgi?id=197603

This method returns whether it is encrypted or not, and what mountpoint it used
to make that decision.
"""
function is_ecryptfs(path::AbstractString; verbose::Bool=false)
    # Canonicalize `path` immediately, and if it's a directory, add a "/" so
    # as to be consistent with the rest of this function
    path = abspath(path)
    if isdir(path)
        path = abspath(path * "/")
    end

    if verbose
        Compat.@info("Checking to see if $path is encrypted...")
    end

    # Get a listing of the current mounts.  If we can't do this, just give up
    if !isfile("/proc/mounts")
        if verbose
            Compat.@info("Couldn't open /proc/mounts, returning...")
        end
        return false
    end
    mounts = String(read("/proc/mounts"))

    # Grab the fstype and the mountpoints
    mounts = [split(m)[2:3] for m in split(mounts, "\n") if !isempty(m)]

    # Canonicalize mountpoints now so as to dodge symlink difficulties
    mounts = [(abspath(m[1]*"/"), m[2]) for m in mounts]

    # Fast-path asking for a mountpoint directly (e.g. not a subdirectory)
    direct_path = [m[1] == path for m in mounts]
    parent = if any(direct_path)
        mounts[findfirst(direct_path)]
    else
        # Find the longest prefix mount:
        parent_mounts = [m for m in mounts if startswith(path, m[1])]
        parent_mounts[indmax(map(m->length(m[1]), parent_mounts))]
    end

    # Return true if this mountpoint is an ecryptfs mount
    return parent[2] == "ecryptfs", parent[1]
end

function check_encryption(workspace_root::AbstractString;
                          verbose::Bool = false)
    # If we've explicitly allowed ecryptfs, just quit out immediately
    global allow_ecryptfs
    if allow_ecryptfs
        return
    end
    msg = []

    is_encrypted, mountpoint = is_ecryptfs(workspace_root; verbose=verbose)
    if is_encrypted
        push!(msg, replace(strip("""
        Will not launch a user namespace runner within $(workspace_root), it
        has been encrypted!  Change your working directory to one outside of
        $(mountpoint) and try again.
        """), "\n" => " "))
    end

    is_encrypted, mountpoint = is_ecryptfs(rootfs_dir(); verbose=verbose)
    if is_encrypted
        push!(msg, replace(strip("""
        Cannot mount rootfs at $(rootfs_dir()), it has been encrypted!  Change
        your rootfs cache directory to one outside of $(mountpoint) by setting
        the BINARYBUILDER_ROOTFS_DIR environment variable and try again.
        """), "\n" => " "))
    end

    is_encrypted, mountpoint = is_ecryptfs(shards_dir(); verbose=verbose)
    if is_encrypted
        push!(msg, replace(strip("""
        Cannot mount rootfs shards within $(shards_dir()), it has been
        encrypted!  Change your shard cache directory to one outside of
        $(mountpoint) by setting the BINARYBUILDER_SHARDS_DIR environment
        variable and try again.
        """), "\n" => " "))
    end

    if !isempty(msg)
        error(join(msg, "\n"))
    end
end

ext: jl


text:
module FixtureLearners

importall Orchestra.Types
importall Orchestra.Util

export MLProblem,
       Classification,
       FeatureClassification,
       NumericFeatureClassification,
       PerfectScoreLearner,
       AlwaysSameLabelLearner,
       fit_and_transform!,
       fit!,
       transform!

abstract MLProblem
abstract Classification <: MLProblem

# NOTE(svs14): Currently hardcoded example.
#              Consider turning into rule-based generator.
train_dataset = [
  1.0        1 "b"  2 "c" "a";
  2.0        2 "b"  3 "c" "a";
  nan(3.0)   3 "b"  4 "c" "a";
  -1.0      -1 "d" -2 "c" "b";
  -2.0      -2 "d" -3 "c" "b";
  nan(-3.0) -3 "d" -4 "c" "b";
  1.0        1 "a"  1 "a" "c";
  2.0        2 "b"  2 "b" "c";
  nan(3.0)   3 "c"  3 "c" "c";
  0.0        0 "e"  1 "a" "d";
  0.0        0 "e"  2 "b" "d";
  nan(0.0)   0 "e"  3 "c" "d";
]
test_dataset = [
  4.0        4 "b"  5 "c" "a";
  nan(5.0)   5 "b"  6 "c" "a";
  -4.0      -4 "d" -5 "c" "b";
  nan(-5.0) -5 "d" -6 "c" "b";
  4.0        4 "d"  4 "d" "c";
  nan(5.0)   5 "e"  5 "e" "c";
  0.0        0 "e"  4 "d" "d";
  nan(0.0)   0 "e"  5 "e" "d";
]

type FeatureClassification <: Classification
  train_instances::Matrix
  test_instances::Matrix
  train_labels::Vector
  test_labels::Vector

  function FeatureClassification()
    train_instances = train_dataset[:, 1:end-1]
    test_instances = test_dataset[:, 1:end-1]
    train_labels = train_dataset[:, end]
    test_labels = test_dataset[:, end]
    new(
      train_instances,
      test_instances,
      train_labels,
      test_labels
    )
  end
end

type NumericFeatureClassification <: Classification
  train_instances::Matrix
  test_instances::Matrix
  train_labels::Vector
  test_labels::Vector

  function NumericFeatureClassification()
    train_instances = convert(Array{Real, 2}, train_dataset[:, [2,4]])
    test_instances = convert(Array{Real, 2}, test_dataset[:, [2,4]])
    train_labels = convert(Array{String, 1}, train_dataset[:, end])
    test_labels = convert(Array{String, 1}, test_dataset[:, end])
    new(
      train_instances,
      test_instances,
      train_labels,
      test_labels
    )
  end
end

function fit_and_transform!(transformer::Transformer, problem::MLProblem, seed=1)
    srand(seed)
    fit!(transformer, problem.train_instances, problem.train_labels)
    return transform!(transformer, problem.test_instances)
end

type PerfectScoreLearner <: TestLearner
  model
  options

  function PerfectScoreLearner(options=Dict())
    default_options = {
      :output => :class,
      :problem => NumericFeatureClassification()
    }
    new(nothing, nested_dict_merge(default_options, options))
  end
end

function fit!(
  psl::PerfectScoreLearner, instances::Matrix, labels::Vector)

  problem = psl.options[:problem]

  dataset = [
    problem.train_instances problem.train_labels;
    problem.test_instances problem.test_labels
  ]
  instance_label_map = [
    dataset[i,1:2] => dataset[i,3] for i=1:size(dataset, 1)
  ]

  psl.model = {
    :map => instance_label_map
  }
end

function transform!(
  psl::PerfectScoreLearner, instances::Matrix)

  num_instances = size(instances, 1)
  predictions = Array(String, num_instances)
  for i in 1:num_instances
    predictions[i] = psl.model[:map][instances[i,:]]
  end
  return predictions
end

type AlwaysSameLabelLearner <: TestLearner
  model
  options

  function AlwaysSameLabelLearner(options=Dict())
    default_options = {
      :output => :class,
      :label => nothing
    }
    new(nothing, nested_dict_merge(default_options, options))
  end
end

function fit!(awsl::AlwaysSameLabelLearner, instances::Matrix, labels::Vector)
  if awsl.options[:label] == nothing
    awsl.model = {
      :label => first(labels)
    }
  else
    awsl.model = {
      :label => awsl.options[:label]
    }
  end
end

function transform!(awsl::AlwaysSameLabelLearner, instances::Matrix)
  return fill(awsl.model[:label], size(instances, 1))
end

end # module


text:
#=
mask.jl
Methods related to an image support mask.
2019-06-22 Jeff Fessler, University of Michigan
2020-06-30 embed! getindex!
=#

export embed, embed!, mask_or, maskit
#export mask_outline
export getindex!

#using ImageFiltering: imfilter, centered
#using SparseArrays: sparse, findnz, AbstractSparseVector

"""
    getindex!(y::AbstractVector, x::AbstractArray{T,D}, mask::AbstractArray{Bool,D})
Equivalent to the in-place `y .= x[mask]` but is non-allocating.

For non-Boolean indexing, just use `@views y .= A[index]`, per
https://discourse.julialang.org/t/efficient-non-allocating-in-place-getindex-for-bitarray/42268
"""
@inline function getindex!(
    y::AbstractVector,
    x::AbstractArray{T,D},
    mask::AbstractArray{Bool,D},
) where {T,D}
    axes(y) == (Base.OneTo(count(mask)),) || throw("y axes $(axes(y))")
    axes(mask) == axes(x) || throw(DimensionMismatch("x vs mask"))
    accum = 1
    @inbounds for i in 1:LinearIndices(mask)[findlast(mask)]
        y[accum] = x[i]
        accum += mask[i]
    end
    return y
end

"""
    mask_or(mask)
compress 3D mask to 2D by logical `or` along `z` direction
"""
@inline mask_or(mask::AbstractArray{Bool,2}) = mask
mask_or(mask::AbstractArray{Bool,3}) =
    dropdims(sum(mask, dims=3) .> 0, dims=3)
mask_or(mask::AbstractArray{Bool}) = throw("ndims(mask) = $(ndims(mask))")

#=
"""
    mask_outline(mask)
return outer boundary of 2D mask (or mask_or for 3D)
"""
function mask_outline(mask::AbstractMatrix{Bool})
    tmp = imfilter(mask, centered(ones(Int32,3,3))) # dilate
#    tmp = tmp[2:end-1,2:end-1] # 'same'
    return (tmp .> 1) .& (.! mask)
end
mask_outline(mask::AbstractArray{Bool,3}) = mask_outline(mask_or(mask))
=#

"""
    embed!(array, v, mask ; filler=0)
Embed vector `v` of length `count(mask)`
into elements of `array` where `mask` is `true`,
setting the remaining elements to `filler` (default 0).
"""
embed!(array::AbstractArray{T,D}, v::AbstractVector{<:Number},
    mask::AbstractArray{Bool,D} ; filler::T = zero(T)
) where {T, D} =
    setindex!(fill!(array, filler), v, mask)

"""
    array = embed(v, mask ; filler=0)

Embed vector `v` of length `count(mask)`
into elements of an array where `mask` is `true`; see `embed!`.
"""
embed(v::AbstractVector{T}, mask::AbstractArray{Bool} ;
    filler::Number = zero(T)) where {T <: Number} =
    embed!(fill(filler, size(mask)), v, mask)

"""
    array = embed(matrix::AbstractMatrix{<:Number}, mask::AbstractArray{Bool})

Embed each column of `matrix` into `mask` then `cat` along next dimension
In:
* `matrix [count(mask) L]`
* `mask [(N)]`

Out:
* `array [(N) L]`
"""
function embed(x::AbstractMatrix{<:Number}, mask::AbstractArray{Bool})
    L = size(x,2)
    out = zeros(eltype(x), prod(size(mask)), L)
    for l=1:L
        out[:,l] = vec(embed(x[:,l], mask))
    end
    reshape(out, size(mask)..., L)
end

#=
this is needed only for "unpacking" sparse system matrices so ignore for now
# image_geom_embed_sparse
# called by image_geom_embed
# function _embed_sparse(x::Array{T}) where {T <: Number}
function embed(x::AbstractSparseVector{<:Number},
                mask::AbstractArray{Bool,N} where N)
    i, v = findnz(x)
    ind = findall(mask)
    j = ind(j)
    return sparsevec(i, j, a, size(x,1), length(mask))
end
=#

"""
    maskit(x::AbstractArray{<:Number}, mask)
opposite of embed
"""
function maskit(x::AbstractArray, mask::AbstractArray{Bool})
    dim = size(x)
    ((ndims(x) >= ndims(mask)) && (size(x)[1:ndims(mask)] == size(mask))) ||
        throw(DimensionMismatch("size x $(size(x)) vs mask $(size(mask))"))
    x = reshape(x, length(mask), :)
    x = x[vec(mask),:] # reshape(mask, prod(_dim()))]
    return length(dim) == ndims(mask) ?
        dropdims(x, dims=2) : # squeeze
        reshape(x, :, dim[(1+ndims(mask)):end]...)
end


text:
using Test
using QuantumOptics
using LinearAlgebra

@testset "stochastic_semiclassical" begin

b_spin = SpinBasis(1//2)
sz = sigmaz(b_spin)
sm = sigmam(b_spin)
sp = sigmap(b_spin)
zero_op = 0*sz
γ = 0.1
noise_op = 0.5γ*sz

H = γ*(sp + sm)
Hs = [noise_op]

ψ0 = spindown(b_spin)

rates = [0.1]
J = [sm]
C = [sm]
Jdagger = dagger.(J)
C .*= rates
Cdagger = dagger.(C)

u0 = ComplexF64[0.1, 0.5]
ψ_sc = semiclassical.State(ψ0, u0)
ρ_sc = dm(ψ_sc)

rates_mat = [0.1 0.05; 0.05 0.1]

dt = 1/100.0
T = [0:0.1:10;]
T_short = [0:dt:dt;]

# Function definitions for schroedinger_semiclassical
function fquantum(t, psi, u)
    return H
end
function fclassical(du, u, psi, t)
    du[1] = u[2]
    du[2] = sin(u[1])*cos(u[1])
end
function fquantum_stoch(t, psi, u)
    Hs
end
function fclassical_stoch(du, u, psi, t)
    du[1] = 0.2*u[1]
    du[2] = 0.2*u[2]
end
function fclassical_stoch2(du, u, psi, t)
    du[1,1] = 0.2*u[1]
    du[2,2] = 0.2*u[2]
end
function fclassical_stoch_ndiag(du, u, psi, t)
    du[1,1] = 0.2*u[1]
    du[1,2] = 0.1*u[1]
    du[2,3] = -0.1u[2]
end

# Function definitions for master_semiclassical
function fquantum_master(t, rho, u)
    H, J, Jdagger
end
function fstoch_q_master(t, rho, u)
    10 .*C, 10 .*Cdagger
end

# Test semiclassical schroedinger
tout, ψt_sc = stochastic.schroedinger_semiclassical(T, ψ_sc, fquantum, fclassical;
            fstoch_quantum=fquantum_stoch, dt=dt, normalize_state=true)
for ψ=ψt_sc
    @test norm(ψ.quantum) ≈ 1.0
end
tout, ψt_sc = stochastic.schroedinger_semiclassical(T_short, ψ_sc, fquantum, fclassical;
            fstoch_classical=fclassical_stoch, dt=dt)
tout, ψt_sc = stochastic.schroedinger_semiclassical(T_short, ψ_sc, fquantum, fclassical;
            fstoch_quantum=fquantum_stoch, fstoch_classical=fclassical_stoch2,
            noise_processes=1,
            noise_prototype_classical=zeros(ComplexF64, 2,2), dt=dt)
tout, ψt_sc = stochastic.schroedinger_semiclassical(T_short, ψ_sc, fquantum, fclassical;
            fstoch_classical=fclassical_stoch_ndiag,
            noise_prototype_classical=zeros(ComplexF64, 2, 3), dt=dt)

# Semiclassical master
tout, ρt = stochastic.master_semiclassical(T, ρ_sc, fquantum_master, fclassical;
            fstoch_quantum=fstoch_q_master, dt=dt)
tout, ρt = stochastic.master_semiclassical(T, ρ_sc, fquantum_master, fclassical;
            fstoch_classical=fclassical_stoch, dt=dt)
tout, ρt = stochastic.master_semiclassical(T, ψ_sc, fquantum_master, fclassical;
            fstoch_quantum=fstoch_q_master, fstoch_classical=fclassical_stoch2,
            noise_prototype_classical=zeros(ComplexF64, 2, 2), dt=dt)

# Test error messages
@test_throws ArgumentError stochastic.schroedinger_semiclassical(T, ψ_sc, fquantum, fclassical)
@test_throws ArgumentError stochastic.schroedinger_semiclassical(T, ψ_sc, fquantum, fclassical;
        fstoch_quantum=fquantum_stoch, fstoch_classical=fclassical_stoch)
@test_throws ArgumentError stochastic.master_semiclassical(T, ρ_sc, fquantum_master, fclassical)

end # testset


text:
# module model_fit

using DifferentialEquations, Flux, Optim, DiffEqFlux, DiffEqSensitivity, Plots
using BlackBoxOptim
using BenchmarkTools
using RecursiveArrayTools # for VectorOfArray
# -------------------------------------------------------------------------

function generate_data(sol)
  t = collect(range(0,stop=10,length=101))

  solution = sum(sol(t),dims=1)
  # randomized = VectorOfArray([(sol(t[i]) + .02randn(size(sol,1))) for i in 1:length(t)])
  # dataset = convert(Array,randomized)
  randomized = solution .+ .01randn(size(solution))
  dataset = randomized
  display(scatter!(t,dataset'))
  return dataset
end

function generate_data(sol,t)

  solution = sum(sol(t),dims=1)
  # randomized = VectorOfArray([(sol(t[i]) + .02randn(size(sol,1))) for i in 1:length(t)])
  # dataset = convert(Array,randomized)
  randomized = solution .+ .01randn(size(solution))
  dataset = randomized
  display(scatter!(t,dataset'))
  return dataset
end

function run_model_2(model,p,data,densities,constants,t) # version of run function with multiple models
    f = (du,u,p,t) -> model(du,u,10 .^ p,t,densities,constants)
    tmp_prob = ODEProblem(f,u₀,tspan,p)
    # tmp_sol = solve(tmp_prob,AutoVern7(KenCarp4()),saveat=t, abstol=1e-8,reltol=1e-8)
    tmp_sol = solve(tmp_prob,Vern7(),saveat=t, abstol=1e-5,reltol=1e-5)
    # stack Σ of solution across n species
    tmp_sol
end

function loss_2(p,data,model,dens,cons,t,unique_t)
    # for i in 1:siz
    tmp_sol = run_model_2(model,p,data,dens,cons,unique_t)

    # @show Σ_sol
    # @show solution
    # if any((s.retcode != :Success for s in solution))
    # @show typeof(solution.retcode)
    if tmp_sol.retcode != :Success
      # @show solution
      L = Inf
    else
      Σ_sol = sum(Array(tmp_sol),dims=1)
      Δy = 0.
      for (i, tᵢ) in enumerate(t)
          ind = unique_t .== tᵢ
          yᵢ_exp = Σ_sol[ind][1]
          yᵢ_obs = data[i]
          Δyᵢ = yᵢ_obs .- yᵢ_exp
          # setindex!(Δy, Δyᵢ, i)
          # Δy[i] = Δyᵢ
          Δy = vcat(Δy,Δyᵢ)
          # @show Δy
          # @show tᵢ, yᵢ_obs, Δyᵢ
      end
      L = sum(abs2, Δy[2:end]) #, Σ_sol
    end
    return L
end

# function callback(p,l,tmp_sol,data,model,dens,cons,t)
function callback(p,l,data,model,dens,cons,t,unique_t)
  @show l, p

  # Σ_sol = run_model_2(model,p,data,dens,cons,unique_t)

  # fig = plot(unique_t,Σ_sol')
  # scatter!(fig,t,data,legend=:bottomright)
  # display(fig)

  false
end

function scan_param_space(p,data,model,dens,cons,t,unique_t,lb,ub)
      #create rand mesh
      n = 10000
      ps = rand(n,size(p,1)) .* 10 .- 8
      errors = zeros(n)
      for i in 1:n
          # error, a = loss(ps[i,:],data,model,dens,cons,t)
          error = loss_2(ps[i,:],data,model,dens,cons,t,unique_t)
          errors[i] = error
          # @show error
      end
      min_ind = argmin(errors)
      min_error = minimum(errors)
      @show findmin(errors)
      p_best = ps[min_ind,:]
      @show p_best
      if size(p,1) == 2
          fig = scatter(ps[:,1],ps[:,2],log10.(errors),
          zcolor = log10.(errors),camera=[0,90])
          scatter!([p_best[1]],[p_best[2]],[log10(min_error)],markercolor=:green,markersize=5)
          display(fig)
      end
      return p_best
end

function optimization(model,data,dens,cons,options,t)
# fit function - given model, data, densities, options get best parameters
  unique_t = unique(t)
  loss_in = (p) -> loss_2(p,data,model,dens,cons,t,unique_t)
  # callback_in = (p,l,tmp_sol) -> callback(p,l,tmp_sol,data,model,dens,cons,t)
  callback_in = (p,l) -> callback(p,l,data,model,dens,cons,t,unique_t)

  # @unpack lb,ub = options
  lb,ub = options

  start = scan_param_space(pinit,data,model,dens,cons,t,unique_t,lb,ub) # !! change behavior so that it doesnt need pinit

  res = DiffEqFlux.sciml_train(loss_in,start,BFGS(),cb=callback_in)

  rates = res.minimizer
  l = res.minimum

  # bootstrap_rates = bootstrapping(rates,data,model,dens,cons,t)
  bootstrap_rates = []

  f = (du,u,p,t) -> model(du,u,p,t,dens,cons) # optimize in log space
  prob = ODEProblem(f,u₀,tspan,10 .^ rates)
  sol = solve(prob,Vern7())
  h = plot!(sol)
  display(h)
  @show model, rates, start

  hes_for = ForwardDiff.hessian(loss_in,rates)
  hes_zyg = Zygote.hessian(loss_in,rates)
  @show hes_for ≈ hes_zyg

  return rates, hes_zyg, l, bootstrap_rates
end

function bootstrapping(rates,data,model,dens,cons,t)
  #create sample
  data = n
  B_size = 10
  # rand_mat = rand(1:length(data),(length(data),B_size))

  rand_mat = [rand(1:length(data),length(data)) for i in 1:B_size]

  bootstrap_rates = zeros(size(rates))
  # bootstrap_rates = []
  for i in 1:B_size
    resampled_data = data[rand_mat[i]]
    resampled_t = t[rand_mat[i]]
    sort_ind = sortperm(resampled_t)
    resampled_data = resampled_data[sort_ind]
    resampled_t = resampled_t[sort_ind]
    unique_t = unique(resampled_t)

    # loss_in = (p) -> loss_2(p,resampled_data,model,dens,cons,resampled_t,unique_t)
    # callback_in = (p,l) -> callback(p,l,resampled_data,model,dens,cons,resampled_t,unique_t)
    loss_in = (p) -> loss_2(p,data,model,dens,cons,t,unique(t))
    callback_in = (p,l) -> callback(p,l,data,model,dens,cons,t,unique(t))

    # @show loss_in(rates)
    res = optimize(loss_in,-9, 2, rates,Fminbox(NelderMead()))
    # @time res = DiffEqFlux.sciml_train(loss_in,rates,BFGS(),
    #                                 lower_bounds = -9,
    #                                 upper_bounds = 2,
    #                                 cb = callback_in,
    #                                 maxiters=50,
    #                                 f_tol = 1e-4)

    # @show res.minimizer
    bootstrap_rates= hcat(bootstrap_rates,res.minimizer)
    # push!(bootstrap_rates,res.minimizer)
    # Σ_sol = run_model_2(model,p,resampled_data,dens,cons,unique_t)
  end

  return bootstrap_rates[:,2:end]
end


text:
module AppendedTriangulationsTests

using Test
using Gridap.ReferenceFEs
using Gridap.Arrays
using Gridap.Geometry
using Gridap.Fields
using Gridap.Integration
using LinearAlgebra: ⋅
using Gridap.CellData
using Gridap.FESpaces
using Gridap.ReferenceFEs
using FillArrays

domain = (0,1,0,1)
partition = (10,10)
model = CartesianDiscreteModel(domain,partition)

ncells = num_cells(model)
nin = ceil(Int,2*ncells/3)
cell_to_mask = fill(false,ncells)
cell_to_mask[1:nin] .= true

grid = get_grid(model)

Ω_in = RestrictedTriangulation(grid,cell_to_mask)
Ω_out = RestrictedTriangulation(grid,.! cell_to_mask)
Ω = lazy_append(Ω_out,Ω_in)
test_triangulation(Ω)

order = 1
degree = 2*order
quad_in = CellQuadrature(Ω_in,degree)
quad_out = CellQuadrature(Ω_out,degree)
quad = CellQuadrature(Ω,degree)
@test isa(quad.trian,AppendedTriangulation)
@test isa(quad.cell_quad,AppendedArray)
@test isa(quad.cell_point,AppendedArray)
@test isa(quad.cell_weight,AppendedArray)

V = TestFESpace(model,ReferenceFE(lagrangian,Float64,order),conformity=:H1)

v(x) = x[1]+x[2]

vh = interpolate(v,V)

e = v - vh
el2 = sqrt(sum(integrate(e*e,quad)))
@test el2 < 1.0e-8

x = get_cell_points(quad)

dv = get_cell_shapefuns(V)
du = get_cell_shapefuns_trial(V)
@test isa(dv(x),AppendedArray)
@test isa(∇(dv)(x),AppendedArray)
@test isa(du(x),AppendedArray)
@test isa(∇(du)(x),AppendedArray)

cellmat = integrate( ∇(dv)⋅∇(du), quad )

@test isa(cellmat,AppendedArray)
@test isa(cellmat.a,Fill)
@test isa(cellmat.b,Fill)

end # module


text:
## using TimeData

## include("/home/chris/research/julia/TimeData/src/TimeData.jl")
include("/home/chris/research/julia/Econometrics/src/Econometrics.jl")
using TimeData
using Winston
using Distributions

tm = readTimedata("/home/chris/research/asset_mgmt/data/datesLogRet.csv")

## get 5 lowest minima stocks
mapFunc(x) = minimum(x)
critFunc(x) = islowest(x, 5)
mostNegativeReturns = TimeData.getVars(tm, mapFunc, critFunc)

#########################
## get moving averages ##
#########################

import Winston.plot
function plot(tm::Timematr)
    plot(core(tm))
end

movAvgs = movAvg(tm, 150)
plot(movAvgs)

marketMeans = rowmeans(tm)
plot(marketMeans)
var(core(marketMeans))

marketMovAvg = movAvg(marketMeans, 400)
plot(marketMovAvg)
Winston.ylim(-0.6, 0.6)

######################################
## mean-shifting induced dependency ##
######################################

marketMovAvg = movAvg(marketMeans, 20)
plot(marketMovAvg)
Winston.ylim(-0.6, 0.6)

## for each day, simulate independent normally distributed values
nDays = size(marketMovAvg, 1)

using Distributions
d = Normal(0, 1.6)
simInnov = rand(d, nDays, 2)
simVals = simInnov + repmat(core(marketMovAvg), 1, 2)
cor(simVals[:, 1], simVals[:,2])

## plot autocor
p = []
for ii=1:20
    if(ii==1)
        p = autocorr(core(tm[1])[:])
    else
        p = autocorr!(core(tm[ii])[:])
    end
end
p
autocorr_red!(simVals[:, 1])

## look at ranks!!

######################################
## correlations for all frequencies ##
######################################

maxAggregation = 50

## include("/home/chris/research/julia/TimeData/src/TimeData.jl")
## tm = readTimedata("/home/chris/research/asset_mgmt/data/datesLogRet.csv")
nObs = size(tm, 2)

frequ = 50
for ii=1:20
    eqVal = true
    x = 0
    y = 0
    while(eqVal)
        (x, y) = [rand(1:nObs) rand(1:nObs)]
        if x !== y
            eqVal = false
        end
    end
    println("first: ", x, " , second: ", y)
    a = Econometrics.frequCorrPlot(tm, x, y, frequ)
    plot(a[:, 1], a[:, 2], "+b")
    ylim(-0.2, 1)
    ## p = FramedPlot(
    ##                xrange=(0,frequ),
    ##                yrange=(-0.3,1))
    ## ## pts = Points(a[:, 1], a[:, 2], kind="filled circle")
    ## pts = Points(a[:, 1], a[:, 2], kind="plus")
    ## add(p, pts)
    fname = string("./pics/corrFrequ", ii, ".png")
    file(fname)
end

    ## create simulated data with no autocorrelation

nSim = 15000
using Distributions
## create two-dimensional normal distribution
rho = 0.5
sigma1 = 1.4
sigma2 = 1.4
covMatr = [sigma1^2 sigma1*sigma2*rho; sigma1*sigma2*rho sigma2^2]
d = MvNormal(covMatr)

simInnov = rand(d, nSim)
tmSim = Timematr(simInnov')
a = TimeData.frequCorrPlot(tmSim, 1, 2, frequ)
plot(a[:, 1], a[:, 2], "+b")
ylim(-0.2, 1)

###########
## GARCH ##
###########

testData = tm[1]

## load econometrics library
include("/home/chris/research/julia/Econometrics/src/Econometrics.jl")

(sigmaHat, paramsHat) = Econometrics.garchFit(testData)

plot(sigmaHat)

sigma = sqrt(var(core(testData)))
initVal = [0, sigma^2, 0, 0, sigma]

params1 =  [0.0338456,
            0.463194,
            0.616413,
            0.193743,
            1.39494]

params2 = [ 0.056642,
           0.0467392,
           0.932111,
           0.0485675,
           1.01403]

Econometrics.garchLLH(paramsHat, core(testData))
Econometrics.garchLLH(params2, core(testData))

d = Normal(0, sigma)
nllh2 = -sum(log(pdf(d, core(testData))))

kk = x -> Econometrics.garchLLH(x, core(testData))
kk(initVal)

garchObj = garchFit(core(tm[1])[:])
sigs = garchObj.sigma
plot(sigs)


text:
module StaleB

# StaleB does not know about StaleC when it is being built.
# However, if StaleC is loaded first, we get `"insert_backedges"`
# invalidations.
using StaleA

# This will be invalidated if StaleC is loaded
useA() = StaleA.stale("hello")

# force precompilation
useA()

end


text:
"""
    DataSubset(data, [indices], [obsdim])

Description
============

Used to represent a subset of some `data` of arbitrary type by
storing which observation-indices the subset spans. Furthermore,
subsequent subsettings are accumulated without needing to access
actual data.

The main purpose for the existence of `DataSubset` is to delay
data access and movement until an actual batch of data (or single
observation) is needed for some computation. This is particularily
useful when the data is not located in memory, but on the hard
drive or some remote location. In such a scenario one wants to
load the required data only when needed.

This type is usually not constructed manually, but instead
instantiated by calling [`datasubset`](@ref),
[`shuffleobs`](@ref), or [`splitobs`](@ref)

In case `data` is some `Tuple`, the constructor will be mapped
over its elements. That means that the constructor returns a
`Tuple` of `DataSubset` instead of a `DataSubset` of `Tuple`.

Arguments
==========

- **`data`** : The object describing the dataset. Can be of any
    type as long as it implements [`getobs`](@ref) and
    [`nobs`](@ref) (see Details for more information).

- **`indices`** : Optional. The index or indices of the
    observation(s) in `data` that the subset should represent.
    Can be of type `Int` or some subtype of `AbstractVector`.

- **`obsdim`** : Optional. If it makes sense for the type of `data`,
    `obsdim` can be used to specify which dimension of `data`
    denotes the observations. It can be specified in a type-stable
    manner as a positional argument (see `?LearnBase.ObsDim`), or
    more conveniently as a smart keyword argument.

Methods
========

- **`getindex`** : Returns the observation(s) of the given
    index/indices as a new `DataSubset`. No data is copied aside
    from the required indices.

- **`nobs`** : Returns the total number observations in the subset
    (**not** the whole data set underneath).

- **`getobs`** : Returns the underlying data that the
    `DataSubset` represents at the given relative indices. Note
    that these indices are in "subset space", and in general will
    not directly correspond to the same indices in the underlying
    data set.

Details
========

For `DataSubset` to work on some data structure, the desired type
`MyType` must implement the following interface:

- `LearnBase.getobs(data::MyType, idx, [obsdim::ObsDimension])` :
    Should return the observation(s) indexed by `idx`.
    In what form is up to the user.
    Note that `idx` can be of type `Int` or `AbstractVector`.

- `StatsBase.nobs(data::MyType, [obsdim::ObsDimension])` :
    Should return the total number of observations in `data`

The following methods can also be provided and are optional:

- `LearnBase.getobs(data::MyType)` :
    By default this function is the identity function.
    If that is not the behaviour that you want for your type,
    you need to provide this method as well.

- `LearnBase.datasubset(data::MyType, idx, obsdim::ObsDimension)` :
    If your custom type has its own kind of subset type, you can
    return it here. An example for such a case are `SubArray` for
    representing a subset of some `AbstractArray`.
    Note: If your type has no use for `obsdim` then dispatch on
    `::ObsDim.Undefined` in the signature.

- `LearnBase.getobs!(buffer, data::MyType, [idx], [obsdim::ObsDimension])` :
    Inplace version of `getobs(data, idx, obsdim)`. If this method
    is provided for `MyType`, then `eachobs` and `eachbatch`
    (among others) can preallocate a buffer that is then reused
    every iteration. Note: `buffer` should be equivalent to the
    return value of `getobs(::MyType, ...)`, since this is how
    `buffer` is preallocated by default.

- `LearnBase.gettargets(data::MyType, idx, [obsdim::ObsDimension])` :
    If `MyType` has a special way to query targets without
    needing to invoke `getobs`, then you can provide your own
    logic here. This can be useful when the targets of your are
    always loaded as metadata, while the data itself remains on
    the hard disk until actually needed.

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)
- Tom Breloff (Github: https://github.com/tbreloff)

Examples
=========

```julia
X, y = MLDataUtils.load_iris()

# The iris set has 150 observations and 4 features
@assert size(X) == (4,150)

# Represents the 80 observations as a DataSubset
subset = DataSubset(X, 21:100)
@assert nobs(subset) == 80
@assert typeof(subset) <: DataSubset
# getobs indexes into the subset
@assert getobs(subset, 1:10) == X[:, 21:30]

# You can also work with data that uses some other dimension
# to denote the observations.
@assert size(X') == (150,4)
subset = DataSubset(X', 21:100, obsdim = :first) # or "obsdim = 1"
@assert nobs(subset) == 80

# To specify the obsdim in a type-stable way, use positional arguments
# provided by the submodule `ObsDim`.
@inferred DataSubset(X', 21:100, ObsDim.First())

# Subsets also works for tuple of data. (useful for labeled data)
subset = DataSubset((X,y), 1:100)
@assert nobs(subset) == 100
@assert typeof(subset) <: Tuple # Tuple of DataSubset

# The lowercase version tries to avoid boxing into DataSubset
# for types that provide a custom "subset", such as arrays.
# Here it instead creates a native SubArray.
subset = datasubset(X, 1:100)
@assert nobs(subset) == 100
@assert typeof(subset) <: SubArray

# Also works for tuples of arbitrary length
subset = datasubset((X,y), 1:100)
@assert nobs(subset) == 100
@assert typeof(subset) <: Tuple # tuple of SubArray

# Split dataset into training and test split
train, test = splitobs(shuffleobs((X,y)), at = 0.7)
@assert typeof(train) <: Tuple # of SubArray
@assert typeof(test)  <: Tuple # of SubArray
@assert nobs(train) == 105
@assert nobs(test) == 45
```

see also
=========

[`datasubset`](@ref),  [`getobs`](@ref), [`nobs`](@ref),
[`splitobs`](@ref), [`shuffleobs`](@ref),
[`KFolds`](@ref), [`BatchView`](@ref), [`ObsView`](@ref),
"""
struct DataSubset{T, I<:Union{Int,AbstractVector}, O<:ObsDimension}
    data::T
    indices::I
    obsdim::O

    function DataSubset{T,I,O}(data::T, indices::I, obsdim::O) where {T,I,O}
        if T <: Tuple
            error("inner constructor should not be called using a Tuple")
        end
        1 <= minimum(indices) || throw(BoundsError(data, indices))
        maximum(indices) <= nobs(data, obsdim) || throw(BoundsError(data, indices))
        new{T,I,O}(data, indices, obsdim)
    end
end

DataSubset(data::T, indices::I, obsdim::O) where {T,I,O} =
    DataSubset{T,I,O}(data, indices, obsdim)

# don't nest subsets
function DataSubset(subset::DataSubset, indices, obsdim)
    @assert subset.obsdim == obsdim
    DataSubset(subset.data, _view(subset.indices, indices), obsdim)
end

function Base.show(io::IO, subset::DataSubset)
    if get(io, :compact, false)
        print(io, "DataSubset{", typeof(subset.data), "} with " , nobs(subset), " observations")
    else
        print(io, summary(subset), "\n ", nobs(subset), " observations")
    end
end

function Base.summary(subset::DataSubset)
    io = IOBuffer()
    print(io, typeof(subset).name.name, "(")
    Base.showarg(io, subset.data, false)
    print(io, ", ")
    Base.showarg(io, subset.indices, false)
    print(io, ", ")
    print(io, replace(string(subset.obsdim), "LearnBase." => ""))
    print(io, ')')
    first(readlines(seek(io,0)))
end

# compare if both subsets cover the same observations of the same data
# we don't care how the indices are stored, just that they match
# in order and values
function Base.:(==)(s1::DataSubset,s2::DataSubset)
    s1.data == s2.data &&
        all(i1==i2 for (i1,i2) in zip(s1.indices,s2.indices)) &&
        s1.obsdim == s2.obsdim
end

Base.length(subset::DataSubset) = length(subset.indices)

Base.lastindex(subset::DataSubset) = length(subset)

Base.getindex(subset::DataSubset, idx) =
    DataSubset(subset.data, _view(subset.indices, idx), subset.obsdim)

nobs(subset::DataSubset) = length(subset)

getobs(subset::DataSubset) =
    getobs(subset.data, subset.indices, subset.obsdim)

function getobs(subset::DataSubset, obsdim::ObsDimension)
    @assert obsdim === subset.obsdim
    getobs(subset)
end

getobs(subset::DataSubset, idx) =
    getobs(subset.data, _view(subset.indices, idx), subset.obsdim)

getobs!(buffer, subset::DataSubset) =
    getobs!(buffer, subset.data, subset.indices, subset.obsdim)

getobs!(buffer, subset::DataSubset, idx) =
    getobs!(buffer, subset.data, _view(subset.indices, idx), subset.obsdim)

# compatibility with nested functions
default_obsdim(subset::DataSubset) = subset.obsdim

nobs(subset::DataSubset, ::ObsDim.Undefined) = nobs(subset)

function nobs(subset::DataSubset, obsdim::ObsDimension)
    @assert obsdim === subset.obsdim
    nobs(subset)
end

function getobs(subset::DataSubset, idx, obsdim::ObsDimension)
    @assert obsdim === subset.obsdim
    getobs(subset, idx)
end

function getobs!(buffer, subset::DataSubset, idx, obsdim::ObsDimension)
    @assert obsdim === subset.obsdim
    getobs!(buffer, subset, idx)
end

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

"""
    datasubset(data, [indices], [obsdim])

Returns a lazy subset of the observations in `data` that
correspond to the given `indices`. No data will be copied except
of the indices. It is similar to calling `DataSubset(data,
[indices], [obsdim])`, but returns a `SubArray` if the type of
`data` is `Array` or `SubArray`. Furthermore, this function may
be extended for custom types of `data` that also want to provide
their own subset-type.

If instead you want to get the subset of observations
corresponding to the given `indices` in their native type, use
`getobs`.

If it makes sense for the type of `data`, `obsdim` can be used to
specify which dimension of `data` denotes the observations. It
can be specified in a type-stable manner as a positional argument
(see `?LearnBase.ObsDim`), or more conveniently as a smart
keyword argument.

see `DataSubset` for more information.
"""
datasubset(data, indices, obsdim) =
    DataSubset(data, indices, obsdim)

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

for fun in (:DataSubset, :datasubset)
    @eval begin
        ($fun)(data, indices; obsdim = default_obsdim(data)) =
            ($fun)(data, indices, convert(LearnBase.ObsDimension,obsdim))

        function ($fun)(data; obsdim = default_obsdim(data))
            nobsdim = convert(LearnBase.ObsDimension,obsdim)
            ($fun)(data, 1:nobs(data, nobsdim), nobsdim)
        end

        # No-op
        ($fun)(subset::DataSubset) = subset

        # allow type-stable way to just provide the obsdim
        ($fun)(data, obsdim::ObsDimension) =
            ($fun)(data, 1:nobs(data, obsdim), obsdim)

        ($fun)(data::Tuple, obsdim::ObsDimension) =
            ($fun)(data, 1:nobs(data, obsdim), obsdim)

        ($fun)(data::Tuple, obsdim::Tuple) =
            ($fun)(data, 1:nobs(data, obsdim), obsdim)

        # map DataSubset over the tuple
        function ($fun)(tup::Tuple)
            _check_nobs(tup)
            map(data -> ($fun)(data), tup)
        end

        function ($fun)(tup::Tuple, indices)
            _check_nobs(tup)
            map(data -> ($fun)(data, indices), tup)
        end

        function ($fun)(tup::Tuple, indices, obsdim::ObsDimension)
            _check_nobs(tup, obsdim)
            map(data -> ($fun)(data, indices, obsdim), tup)
        end

        @generated function ($fun)(tup::Tuple, indices, obsdims::Tuple)
            N = length(obsdims.types)
            quote
                _check_nobs(tup, obsdims)
                # This line generates a tuple of N elements:
                # (datasubset(tup[1], indices, obsdims[1]), datasu...
                $(Expr(:tuple, (:(($($fun))(tup[$i], indices, obsdims[$i])) for i in 1:N)...))
            end
        end
    end
end

# --------------------------------------------------------------------
# Arrays

datasubset(A::SubArray; kw...) = A

# catch the undefined setting for consistency.
# should never happen by accident
datasubset(A::AbstractArray, idx, obsdim::ObsDim.Undefined) =
    throw(MethodError(datasubset, (A, idx, obsdim)))

datasubset(A::AbstractSparseArray, idx, obsdim::ObsDimension) =
    DataSubset(A, idx, obsdim)

@generated function datasubset(A::AbstractArray{T,N}, idx, obsdim::ObsDimension) where {T,N}
    @assert N > 0
    if N == 1 && idx <: Integer
        :(view(A, idx))
    elseif obsdim <: ObsDim.First
        :(view(A, idx, $(fill(:(:),N-1)...)))
    elseif obsdim <: ObsDim.Last || (obsdim <: ObsDim.Constant && obsdim.parameters[1] == N)
        :(view(A, $(fill(:(:),N-1)...), idx))
    else # obsdim <: ObsDim.Constant
        DIM = obsdim.parameters[1]
        DIM > N && throw(DimensionMismatch("the given obsdim=$DIM is greater than the number of available dimensions N=$N"))
        :(view(A, $(fill(:(:),DIM-1)...), idx, $(fill(:(:),N-DIM)...)))
    end
end


text:
"""
    getdata_arlett_animation([seed]) -> DataFrame

Start and finish one Arlett walk and return every tenth coordinate in a DataFrame.

Other simulation controls include:
    set_FLOW_BIAS()
    set_CATALASE_ON_WALLS()
"""
function getdata_arlett_animation(seed::Int=randseed())
    init_N = NUMBER_OF_WALKS
    set_NUMBER_OF_WALKS(1)
    arlett_print(seed)
    println("NOTE: # of steps (max) is ignored")
    Random.seed!(seed)
    x_arr, y_arr = ([], [])
    steps_sofar = 0
    every_tenth_frame = 10
    # Record every tenth xy position instead of run_sim!()
    peroxidexy = spawnrandompoint()
    while peroxidexy != undef
        if steps_sofar % every_tenth_frame == 0
            x, y = peroxidexy
            compressedX = convert(Float16, x)
            compressedY = convert(Float16, y)
            push!(x_arr, compressedX)
            push!(y_arr, compressedY)
        end
        peroxidexy, collisiontype = one_step!(peroxidexy)
        steps_sofar += 1
    end
    df = DataFrame(x_coordinate = x_arr,
                    y_coordinate = y_arr)
    current_df(df)
    current_seed(seed)
    current_path("out/animations/")
    set_NUMBER_OF_WALKS(init_N)
    return df
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
#=
Construct a 100 x 100 grid and run discreteVI on it. Also run global approx VI where any integral point in the grid
can be sampled. Compare the values at the integer grid points. Confirm that the max difference reduces with increasing
number of samples.
=#

# Feature vector conversion
# For point (x,y), fv is (1 x y xy x^2 y^2)
function GlobalApproximationValueIteration.convert_featurevector(::Type{SVector{10, Float64}}, s::GWPos, mdp::SimpleGridWorld)

    x = s[1]
    y = s[2]

    v = SVector{10, Float64}(1, x, y, x*y, x^2, y^2, x^3, x^2*y, x*y^2, y^3)

    return v
end

# Sample a specific integral point in the grid
function GlobalApproximationValueIteration.sample_state(mdp::SimpleGridWorld, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    x = rand(rng, 1:mdp.size[1])
    y = rand(rng, 1:mdp.size[2])

    return GWPos(x, y)
end

function test_absolute_error()

    rng = MersenneTwister(1234)

    # Set solver and grid world parameters
    MAX_ITERS = 500
    NUM_SAMPLES = 1000

    SIZE_X = 5
    SIZE_Y = 5
    REWARD_COV_PROB = 0.4

    rewards = Dict{GWPos, Float64}()
    for x = 1:SIZE_X
        for y = 1:SIZE_Y
            if rand(rng) < REWARD_COV_PROB
                rewards[GWPos(x, y)] = 1
            end
        end
    end

    # Create the MDP for a typical grid world
    mdp = SimpleGridWorld(size=(SIZE_X, SIZE_Y), rewards=rewards)

    # Create the linear function approximation with 10 weight parameters, initialized to zero
    lin_gfa = LinearGlobalFunctionApproximator(zeros(10))

    # Initialize the global approximation solver with the linear approximator and solve the MDP to obtain the policy
    gfa_solver = GlobalApproximationValueIterationSolver(lin_gfa, num_samples=NUM_SAMPLES, max_iterations=MAX_ITERS, verbose=true, fv_type=SVector{10, Float64})
    gfa_policy = solve(gfa_solver, mdp)

    error_arr = Vector{Float64}(undef, 0)

    for state in states(mdp)

        full_val = value(policy, state)
        approx_val = value(gfa_policy, state)
        abs_diff = abs(full_val - approx_val)

        push!(error_arr, abs_diff)
    end

    @show mean(error_arr)
    @show maximum(error_arr)

    return (mean(error_arr) < 0.04 && maximum(error_arr) < 0.35)
end

function test_relative_error()

    rng = MersenneTwister(2378)

    # Attempt to approximate globally with N samples and M iterations. As N increases, the average error should decrease
    MAX_ITERS = 500
    NUM_SAMPLES_LOW = 30
    NUM_SAMPLES_HI = 1000

    # Grid probabilities
    SIZE_X = 5
    SIZE_Y = 5
    REWARD_COV_PROB = 0.4

    # Generate a block of reward states from 40,40 to 60,60
    rewards = Dict{GWPos, Float64}()
    for x = 1:SIZE_X
        for y = 1:SIZE_Y
            if rand(rng) < REWARD_COV_PROB
                rewards[GWPos(x, y)] = 1
            end
        end
    end

    # Create MDP
    mdp = SimpleGridWorld(size=(SIZE_X, SIZE_Y), rewards=rewards)

    lin_gfa_1 = LinearGlobalFunctionApproximator(zeros(10))
    lin_gfa_2 = LinearGlobalFunctionApproximator(zeros(10))

    solver_low = GlobalApproximationValueIterationSolver(lin_gfa_1; num_samples=NUM_SAMPLES_LOW, max_iterations=MAX_ITERS, verbose=true, fv_type=SVector{10, Float64})
    solver_hi = GlobalApproximationValueIterationSolver(lin_gfa_2; num_samples=NUM_SAMPLES_HI, max_iterations=MAX_ITERS, verbose=true, fv_type=SVector{10, Float64})

    policy_low = solve(solver_low, mdp)
    policy_hi = solve(solver_hi, mdp)

    # Now solve with dicrete VI
    solver = ValueIterationSolver(max_iterations=1000, verbose=true)
    policy = solve(solver, mdp)

    err_arr_low = Vector{Float64}(undef, 0)
    err_arr_hi = Vector{Float64}(undef, 0)

    for state in states(mdp)

        full_val = value(policy, state)

        approx_val_low = value(policy_low, state)
        approx_val_hi = value(policy_hi, state)

        push!(err_arr_low, abs(full_val - approx_val_low))
        push!(err_arr_hi, abs(full_val - approx_val_hi))

    end

    @show mean(err_arr_low), mean(err_arr_hi)
    @show maximum(err_arr_low), maximum(err_arr_hi)

    return (mean(err_arr_low) > mean(err_arr_hi) && maximum(err_arr_low) > maximum(err_arr_hi))
end

# @test test_absolute_error() == true
@test test_relative_error() == true


text:
using TropicalNumbers, TropicalYao.Reversible.LogLikeNumbers
using Test, NiLang.AD, NiLang, Random
using ForwardDiff
using NiLang.NiLangCore: default_constructor

@testset "GVar" begin
    x = Tropical(0.4)
	@test (~Tropical)(x) == 0.4
    @test GVar(x) isa Tropical{<:GVar}
    @test grad(x) isa Tropical{Float64}
    @test (~GVar)(GVar(Tropical(0.4))) == x

    x = GVar(Tropical(0))
    @test zero(x) isa Tropical{GVar{Int, Int}}

    x = ULogarithmic(0.4)
    @test GVar(x) isa ULogarithmic{<:GVar}
    @test grad(x) isa ULogarithmic{Float64}
    @test (~GVar)(GVar(ULogarithmic(0.4))) == x
end

@testset "basic instructions" begin
	x, y = Tropical(1), Tropical(2)
	@instr x *= y
	@test x == Tropical(3)
	@test y == Tropical(2)
	x, y, bk = Tropical(3), Tropical(2), false
	@instr i_unsafe_addto(y, x, bk)
	@test x == Tropical(2)
	@test y == Tropical(3)
	@test bk == true
end

@testset "basic instructions, ULogarithmic" begin
    x, y = default_constructor(ULogarithmic{Int}, 1),
    default_constructor(ULogarithmic{Int}, 2)
	@instr x *= y
    @test x == default_constructor(ULogarithmic{Int}, 3)
    @test y == default_constructor(ULogarithmic{Int}, 2)

	@test PlusEq(gaussian_log)(1.0, 2.0) == (1.0+log(1+exp(2.0)), 2.0)
	@test check_grad(PlusEq(gaussian_log), (1.0, 2.0); iloss=1)

    x, y,z = default_constructor(ULogarithmic{Float64}, 7.0),
    default_constructor(ULogarithmic{Float64}, 2.0),
    default_constructor(ULogarithmic{Float64}, 3.0)
	@instr x *= y + z
	@test check_inv(MulEq(+), (x, y, z))
	@test x.log ≈ log(exp(7.0) * (exp(2.0) + exp(3.0)))
    x, y,z = default_constructor(ULogarithmic{Float64}, 7.0),
    default_constructor(ULogarithmic{Float64}, 5.0),
    default_constructor(ULogarithmic{Float64}, 3.0)
	@instr x *= y - z
	@test x.log ≈ log(exp(7.0) * (exp(5.0) - exp(3.0)))
	function muleq(f, x::T, y, z) where T
        x = default_constructor(ULogarithmic{T}, x)
        y = default_constructor(ULogarithmic{T}, y)
        z = default_constructor(ULogarithmic{T}, z)
		x *= f(y, z)
		x.log
	end
	g1 = ForwardDiff.gradient(arr->muleq(+, arr...), [7.0, 5.0, 3.0])
    x, y,z = default_constructor(ULogarithmic{Float64}, 7.0),
    default_constructor(ULogarithmic{Float64}, 5.0),
    default_constructor(ULogarithmic{Float64}, 3.0)
	@instr (MulEq(+))(x, y, z)
	@instr GVar(x)
	@instr GVar(y)
	@instr GVar(z)
	@instr x.log.g += 1
	@instr (~MulEq(+))(x, y, z)
	@test grad(x.log) ≈ g1[1]
	@test grad(y.log) ≈ g1[2]
	@test grad(z.log) ≈ g1[3]

	g2 = ForwardDiff.gradient(arr->muleq(-, arr...), [7.0, 5.0, 3.0])
    x, y,z = default_constructor(ULogarithmic{Float64}, 2.0),
    default_constructor(ULogarithmic{Float64}, 5.0),
    default_constructor(ULogarithmic{Float64}, 3.0)
	@instr (MulEq(-))(x, y, z)
	@instr GVar(x)
	@instr GVar(y)
	@instr GVar(z)
	@instr x.log.g += 1
	@instr (~MulEq(-))(x, y, z)
	@test grad(x.log) ≈ g2[1]
	@test grad(y.log) ≈ g2[2]
	@test grad(z.log) ≈ g2[3]
end


text:
function interp1extra(Xgrid,Ygrid,val,interpObj)
# I just don't know why the Grid guys didn't add linear extrapolation!
# interpObj: You should call the interpolation function first, this
#            program just perform the linear extrapolation if you are
#            out of the top boundary of Xgrid
# Example:
#   intfV  = InterpIrregular(Agrid1,EV1, BCnil, InterpLinear);
#   VA1    = interp1extra(Agrid1,EV1,A1,intfV); # instead of intfV[A1]

    if val<Xgrid[1] # Below the minimum point
    	Vval=Ygrid[1];
    elseif val>Xgrid[1] && val<Xgrid[length(Xgrid)]  # Interior Solution
		#interpObj  = InterpIrregular(Xgrid,Ygrid, BCnil, InterpLinear); #Already done outside this function
        Vval=interpObj[val]
    else  # Above the maximum point (linear extrapolation: li=x*b+c)
        x2=Xgrid[length(Xgrid)]
        x1 =Xgrid[length(Xgrid)-1]

        l2 =Ygrid[length(Xgrid)]
        l1 =Ygrid[length(Xgrid)-1]

        b=(l2-l1)/(x2-x1);
        c=l2-b*x2;

    	Vval=val*b+c;
    end

    return Vval

end

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
using Buffers
using Test
using Zygote

@testset "Buffers" begin
    @testset "basic" begin
        w = rand(2, 3)
        b = rand(2)
        c = rand(2)

        buff = Buffer([])
        @test length(buff) == 0

        push!(buff, w)
        push!(buff, b)
        push!(buff, c)
        @test length(buff) == 3
        @test buff[1] === w
        @test buff[2] === b
        @test buff[3] === c

        deleteat!(buff, 2)
        @test length(buff) == 2
        @test buff[1] === w
        @test buff[2] === c

        deleteat!(buff, 1)
        @test length(buff) == 1
        @test buff[1] === c
    end

    @testset "gradcheck" begin
        @test gradient([1, 2, 3]) do x
            b = Buffer(x)
            b[:] = x
            return sum(copy(b))
        end == ([1, 1, 1],)

        function vstack(xs)
            buf = Buffer(xs, length(xs), 5)
            for i = 1:5
                buf[:, i] = xs
            end
            return copy(buf)
        end

        function _scaledpairedcumsum(xs)
            n = length(xs)
            n > 0 || return copy(xs)
            buf = Buffer(xs)
            buf[1:2] = xs[1:2]
            for i in 3:2:(n - 1)
                buf[i:(i + 1)] = buf[(i - 2):(i - 1)] + i * xs[i:(i + 1)]
            end
            return copy(buf)
        end

        @test gradient(x -> sum(vstack(x)), [1, 2, 3]) == ([5, 5, 5],)

        @test Zygote.pullback(_scaledpairedcumsum, 1:4)[2]([0, 0, 0, 1]) == ([0, 1, 0, 3],)
        @test Zygote.pullback(_scaledpairedcumsum, 1:4)[2]([0, 0, 1, 0]) == ([1, 0, 3, 0],)
        @test Zygote.pullback(_scaledpairedcumsum, 1:4)[2]([1, 1, 1, 1]) == ([2, 2, 3, 3],)

        buf = Buffer([1, 2, 3])
        buf[1] = 1
        copy(buf)
        @test_throws ErrorException buf[1] = 1
        @test eltype(buf) === Int
        @test length(buf) === 3
        @test ndims(buf) === 1
        @test size(buf) === (3,)
        @test size(buf, 2) === 1
        @test axes(buf) == (1:3,)
        @test axes(buf, 2) == 1:1
        @test eachindex(buf) == 1:3
        @test stride(buf, 2) === 3
        @test strides(buf) === (1,)
        @test collect(buf) == collect(copy(buf))
        @test Base.IteratorSize(typeof(buf)) == Base.IteratorSize(typeof([1, 2, 3]))

        @test gradient([1, 2, 3]) do xs
            b = Buffer(xs)
            b .= xs .* 2
            return sum(copy(b))
        end == ([2, 2, 2],)

        @test gradient(2) do x
            b = Buffer([])
            push!(b, x)
            push!(b, 3)
            prod(copy(b))
        end == (3,)
    end

    VERSION > v"1.3-" && include("threads.jl")
end


text:
abstract type QState end

function probability(state::QState, config::Vector{Int})
	if numQubits(state) != length(config)
		error("Given configuration dimension does not match the state")
	end
	if any(config .> 1) || any(config .<0)
		error("Invalid configuration input.")
	end
	ρ_config = ρ(state, config)
	return Real(ρ_config*conj(ρ_config))
end


text:
module MagnusNudecker

using LinearAlgebra

function selection_matrix(m::Int)
    cd = round(Int, m*(m+1)/2)
    rd = m*m

    D = zeros(Int, rd, cd)
    @inbounds for j = 1:m
        for i = 1:j
            r_ij = round(Int, (j*j-j)/2 + i)
            h_ij = round(Int, m*(j-1) + i)
            h_ji = round(Int, m*(i-1) + j)
            D[h_ij, r_ij] = 1
            D[h_ji, r_ij] = 1
        end
    end
    D
end

function Hmat(i, j, m, n)
    H = zeros(m,n)
    H[i,j] = 1
    H
end

function commutation_matrix(m::Int, n::Int = m)
    K = zeros(m*n, m*n)
    for i = 1:m
        for j = 1:n
            Hij = Hmat(i, j, m, n)
            K += kron(Hij, Hij')
        end
    end
    K
end

function Nm(m::Int)
    Km = commutation_matrix(m,m)
    for j = 1:m^2
        Km[j,j] = Km[j,j] + 1.0
    end
    Km./2
end

function Emat(i, j, m)
    E = zeros(m,m)
    E[i,j] = 1
    E
end

function Tmat(i, j, m)
    if i ≠ j
        Emat(i, j, m) + Emat(j, i, m)
    else
        Emat(i,i,m)
    end
end

function uvec(i, j, m)
    r = round(Int, m*(m+1)/2)
    u = zeros(r, 1)
    idx = round(Int, (j-1)*m + i - j*(j-1)/2)
    u[idx,1] = 1
    u
end

function Dm(m::Int)
    D = zeros(m^2, round(Int, m*(m+1)/2))
    for j = 1:m
        for i = j:m
            D += vec(Tmat(i, j, m))*uvec(i,j,m)'
        end
    end
    D
end

function Lmt(m::Int)
    D = Dm(m)
    for j = 1:(m-1)
        D[j*m+(1:j),:] = 0.0
    end
    D
end

function Lm(m::Int)
    Lmt(m)'
end

export Lm, Dm, Nm, selection_matrix, commutation_matrix

end # module


text:
@testset "1.two-sum.jl" begin
@test two_sum([2, 7, 11, 15], 9) == (1, 2)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct ComputeOperationsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Gets a list of compute operations.
Param: api_version::String (required)
Return: ComputeOperationListResult
"""
function _swaggerinternal_operationsList(_api::ComputeOperationsApi, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ComputeOperationListResult, "/providers/Microsoft.Compute/operations", ["azure_auth"])
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function operationsList(_api::ComputeOperationsApi, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_operationsList(_api, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function operationsList(_api::ComputeOperationsApi, response_stream::Channel, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_operationsList(_api, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export operationsList


text:
# Julia1_3.1.jl  Conditional Control Structures
#
# Boolean values are
true, false # are Boolean values (type Bool)
#   (one byte values that are 1 or 0).
# Flow control - Conditionals if-block
# if-elseif-else-end can be nested to any level
x=2; y=1
if x<y
    println("x is less than y")
elseif x>y
    println("x is greater than y")
else
    println("x is equal to y")
end

# you can have a chain of elseifs for
# handling different cases

# watch out for x=1//3; y=1/3, won't be equal
# however rationalize(1/3)==1//3
# 1 ≈ 0.999999999 where ≈ means approximately equal
# 3.825==153/40 this is false
# 3.825≈153/40  this is true

# one liner conditional borrowed from C
# cond? true-expr:false-expr
x=-2
y = x<0 ? -x : x  # gets absolute value of x
# same as y=abs(x)

# other one liner conditional executions
# if <cond>; <expression>; end can be reduced to:
#  where && read as "and then"
x=-3
x<0 && (x=-x) # if x<0 and then x set to -x
# if !<cond>; <expression>; end   can be reduced to:
#  where || read as "or else"
x=-4
!(x<0) || (x=-x) # if not x<0 bypass else x set to -x
# in this case both perform x=abs(x)

# the above operations use a side effect of eval
# when used in this fashion
#  It can obsfuscate code

# Now as clarity intended.
# logical operator && is "and"
# logical operator || is non-exclusice "or"
# a && b is true, if and only if, both a and b are true.
# a || b is true, if either a or b or both are true.
 3<5 && 7>6  # is true
 3>5 || 2<5  # is true, note had to continue to eval
 3>5 && 7>6  # 3>5 false so needn't eval rest

# numeric(and alphanumeric)comparison operators
# < less than, > greater than, == equal to
# <= or ≤  less than or equal,     ≤ is Char(0x2264)
# >= or ≥  greater than or equal,  ≥ is Char(0x2265)
# != or ≠  not equal,              ≠ is Char(0x2260)
# != or ≠  same as !(x==y) which is often preferred
# many other operators for checking types, etc.

# operators can be chained:
1<2≤2<3==3>2≥1==1<3≠5
# results in true

# logical Operators
# && = AND, || = OR,
# xor(x,y) XOR true if x or y true but not both!
xor(true,true)

x=1; y=2
x<5 && y>2
# false
0.0≤x≤1.0 && 1.0≤y≤2.0
# true


text:
using IndirectTrajOpt
using Test
using SafeTestsets

@time begin
@time @safetestset "CR3BP Jacobian" begin include("testCr3bpEomJac.jl") end
end


text:
#
# separable.jl --
#
# Implementation of efficient separables filters by means of the van
# Herk-Gil-Werman algorithm.
#
#------------------------------------------------------------------------------
#
# This file is part of the `LocalFilters.jl` package licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2018, Éric Thiébaut.
#

module Separable

using Compat

using ..LocalFilters

using ..LocalFilters:
    IndexInterval,
    axes,
    cartesianregion

import ..LocalFilters:
    dilate!,
    dilate,
    erode!,
    erode,
    localfilter!,
    localfilter

"""
# Local filter by the van Herk-Gil-Werman algorithm

```julia
localfilter!([dst = A,] A, dims, op, rngs [, w])
```

overwrites the contents of `dst` with the result of applying van
Herk-Gil-Werman algorithm to filter array `A` along dimension(s) `dims` with
(associative) binary operation `op` and contiguous structuring element(s)
defined by the interval(s) `rngs`.  Optional argument `w` is a workspace array
which is automatically allocated if not provided; otherwise, it must be a
vector of same element type as `A` which is resized (with [`resize!`](@ref)) as
needed.  The destination `dst` must have the same indices as the source `A`
(see [`axes`](@ref)).  Operation can be done in-place; that is, `dst` and `A`
can be the same.

Argument `dims` specifies along which dimension(s) of `A` the filter is to be
applied, it can be a single integer, several integers or a colon `:` to specify
all dimensions.  Dimensions are processed in the order given by `dims` (the
same dimension may appear several times) and there must be a matching interval
in `rngs` to specify the structuring element (except that if `rngs` is a single
interval, it is used for every dimension in `dims`).  An interval is either an
integer or an integer unit range in the form `kmin:kmax` (an interval specified
as a single integer, say `k`, is the same as specifying `k:k`).

Assuming mono-dimensional arrays `A` and `dst`, the single filtering pass:

```julia
localfilter!(dst, A, :, op, rng)
```

yields:

```
dst[j] = A[j-kmax] ⋄ A[j-kmax+1] ⋄ A[j-kmax+2] ⋄ ... ⋄ A[j-kmin]
```

for all `j ∈ [first(axes(A,1)):last(axes(A,1))]`, with `x ⋄ y = op(x, y)`,
`kmin = first(rng)` and `kmax = last(rng)`.  Note that if `kmin = kmax = k`
(which occurs if `rng` is a simple integer), the result of the filter is to
operate a simple shift by `k` along the corresponding dimension and has no
effects if `k = 0`.  This can be exploited to not filter some dimension(s).

The out-place version, allocates the destination array and is called as:

```julia
localfilter([T,] A, dims, op, rngs [, w])
```

with `T` the element type of the result (by default `T = eltype(A)`).

## Examples

The in-place *morphological erosion* (local minimum) of the array `A` on a
centered structuring element of width 7 in every dimension can be applied by:

```julia
localfilter!(A, :, min, -3:3)
```

Index interval `0` may be specified to do nothing along the corresponding
dimension.  For instance, assuming `A` is a three-dimensional array:

```julia
localfilter!(A, :, max, (-3:3, 0, -4:4))
```

overwrites `A` its *morphological dilation* (*i.e.* local maximum) in a
centered local neighborhood of size `7×1×9` (nothing is done along the second
dimension).  The same result may be obtained with:

```julia
localfilter!(A, (1,3), max, (-3:3, -4:4))
```

where the second dimension is omitted from the list of dimensions.

The *local average* of the two-dimensional array `A` on a centered
structuring element of size 11×11 can be computed as:

```julia
localfilter(A, :, +, (-5:5, -5:5))*(1/11)
```

## Efficiency and restrictions

The van Herk-Gil-Werman algorithm is very fast for rectangular structuring
elements.  It takes at most 3 operations to filter an element along a given
dimension whatever the width `p` of the considered local neighborhood.  For
`N`-dimensional arrays, the algorithm requires only `3N` operations per element
instead of `p^N - 1` operations for a naive implementation.  This however
requires to make a pass along each dimension so memory page faults may reduce
the performances.  This is somewhat attenuated by the fact that the algorithm
can be applied in-place.  For efficient multi-dimensional out-of-place
filtering, it is recommended to make the first pass with a fresh destination
array and then all other passes in-place on the destination array.

To apply the van Herk-Gil-Werman algorithm, the structuring element must be
separable along the dimensions and its components must be contiguous.  In other
words, the algorithm is only applicable for `N`-dimensional rectangular
neighborhoods.  The structuring element may however be off-centered by
arbitrary offsets along each dimension.

To take into account boundary conditions (for now only nearest neighborhood is
implemented) and allow for in-place operation, the algorithm allocates a
workspace array.

"""
function localfilter!(dst::AbstractArray{T,N},
                      A::AbstractArray{<:Any,N}, d::Int,
                      op::Function, kmin::Int, kmax::Int,
                      w::Vector{T}) where {T,N}
    #
    # A monodimensional local filter involving binary operation ⋄ on a
    # neighborhood [kmin:kmax] yields the following array B when filtering the
    # array A:
    #
    #     B[j] = A[j-kmax] ⋄ A[j-kmax+1] ⋄ ... ⋄ A[j-kmin]
    #
    # Assuming ⋄ is an associative binary operation, van Herk-Gil-Werman (HGW)
    # algorithm consists in building 2 auxilliary arrays, R and S, as follows:
    #
    #     R[1] = A[j-kmin-1]
    #     R[i] = A[j-kmin-i] ⋄ R[i-1]   for i ∈ [2:p-1]
    #     S[1] = A[j-kmin]
    #     S[i+1] = S[i] ⋄ A[j-kmin+i]   for i ∈ [1:p-1]
    #
    # with p = kmax - kmin + 1 the length of the structuring element.
    # Graphically (for p = 7):
    #
    #             j-kmax      j-kmin
    #              ↓           ↓
    #     R[p-1] ← □⋄□⋄□⋄□⋄□⋄□ □             → S[1]
    #     R[p-2] ←   □⋄□⋄□⋄□⋄□ □⋄□           → S[2]
    #     R[p-3] ←     □⋄□⋄□⋄□ □⋄□⋄□         → S[3]
    #      ...           □⋄□⋄□ □⋄□⋄□⋄□          ...
    #     R[2]   ←         □⋄□ □⋄□⋄□⋄□⋄□     → S[p-2]
    #     R[1]   ←           □ □⋄□⋄□⋄□⋄□⋄□   → S[p-1]
    #                          □⋄□⋄□⋄□⋄□⋄□⋄□ → S[p]
    #
    # where each □ represents a cell in A.  Thanks to the recursion, this only
    # takes 2p - 3 operations ⋄. Then, assuming the result is stored in array
    # B, we have:
    #
    #     B[j]     = R[p-1]⋄S[1]
    #     B[j+1]   = R[p-2]⋄S[2]
    #     ...
    #     B[j+p-2] = R[1]⋄S[p-1]
    #     B[j+p-1] = S[p]
    #
    # which requires p - 1 more operations ⋄.  So the total number of
    # operations ⋄ is 3p - 4 for computing p values in B, hence 3 - 4/p
    # operations per cell on average.
    #
    # The following implementation exploits a single workspace array to
    # temporarily store A[j] to allow for in-place appllication of the filter
    # and with suplementary values to account for boundary conditions.  The
    # workspace array is also used to store R[i] while the values of S[i] are
    # computed on the fly.
    #

    # Check arguments and get dimensions.
    1 ≤ d ≤ N || throw(ArgumentError("out of bounds dimension index"))
    kmin ≤ kmax || throw(ArgumentError("invalid structuring element interval"))
    inds = axes(A)
    axes(dst) == inds || throw(DimensionMismatch("source and destination must have the same indices"))
    jmin, jmax = first(inds[d]), last(inds[d])
    if kmin == kmax == 0 || jmin > jmax
        # Nothing to do!
        return copyto!(dst, A)
    end

    # Get index bounds on other parts of input array A.
    R1 = cartesianregion(inds[1:d-1])
    R2 = cartesianregion(inds[d+1:N])

    if kmin == kmax
        # Perform a simple shift: `dst[j] = A[j-k]`.  To allow for in-place
        # operation without using any temporaries, we walk along the dimension
        # in the forward direction if the shift is negative and in the
        # backward direction otherwise.
        k = kmin
        if k < 0
            _shiftarray!(dst, A, R1, jmin:jmax, R2, k)
        else
            _shiftarray!(dst, A, R1, jmax:-1:jmin, R2, k)
        end
    else
        # Apply van Herk-Gil-Werman algorithm.
        _localfilter!(dst, A, R1, jmin, jmax, R2, op, kmin, kmax, w)
    end
    return dst
end

# Private methods to break type uncertainty.

function _shiftarray!(dst::AbstractArray{<:Any,N}, A::AbstractArray{<:Any,N},
                      R1, rng::AbstractRange{Int}, R2, k::Int) where {N}
    jmin, jmax = minimum(rng), maximum(rng)
    @inbounds for J2 in R2, J1 in R1
        @simd for j in rng
            jp = clamp(j - k, jmin, jmax)
            dst[J1,j,J2] = A[J1,jp,J2]
        end
    end
    return nothing
end

function _localfilter!(dst::AbstractArray{T,N},
                       A::AbstractArray{<:Any,N},
                       R1, jmin::Int, jmax::Int, R2,
                       op::Function, kmin::Int, kmax::Int,
                       w::Vector{T}) where {T,N}
    n = jmax - jmin + 1 # length of dimension in A
    p = kmax - kmin + 1 # length of neighborhood
    imin, imax = p, workspacelength(n, p) # range for storing A in W
    length(w) ≥ imax || resize!(w, imax)
    off = imin - jmin + kmax # offset such that W[j-k+off] ≡ A[j-k]
    m = off - kmin # W[j+m] ≡ A[j-kmin]
    pm1 = p - 1

    @inbounds for J2 in R2, J1 in R1
        # Fill the workspace W[imin:imax] with A[jmin-kmax:jmax-kmin] taking
        # care of boundary conditions (here we assume nearest neighbor
        # conditions).
        @simd for i in imin:imax
            j = clamp(i - off, jmin, jmax)
            w[i] = A[J1,j,J2]
        end

        # Process the input by blocks of as much as p elements (less at the end
        # of the range).
        for j = jmin:p:jmax
            # Compute auxilliary array W[1:p-1] ≡ R[1:p-1].
            jpm = j + m
            w[1] = w[jpm-1] # R[1] = A[j-kmin-1]
            @simd for i in 2:pm1
                # R[i] = A[j-kmin-i] ⋄ R[i-1]   for i ∈ [2:p-1]
                w[i] = op(w[jpm-i], w[i-1])
            end

            # Apply the recursion to compute at least 1 and at most p resulting
            # values.
            s = w[jpm] # S[1] = A[j-kmin]
            dst[J1,j,J2] = op(w[pm1], s) # B[j] = R[p-1]⋄S[1]
            @simd for i in 1:min(p-2, jmax-j)
                s = op(s, w[jpm+i]) # S[i+1] = S[i]⋄A[j-kmin+i] for i ∈ [p-1:1]
                dst[J1,j+i,J2] = op(w[pm1-i], s) # B[j+i] = R[p-1-i]⋄S[i+1]
            end
            if j + pm1 ≤ jmax
                # B[j+p-1] = S[p-1]⋄A[j-kmin+p-1]
                dst[J1,j+pm1,J2] = op(s,w[jpm+pm1])
            end
        end
    end
    return nothing
end

# Provide destination.

function localfilter(A::AbstractArray{T,N},
                     dims::Union{Colon, Integer, Tuple{Vararg{Integer}},
                                 AbstractVector{<:Integer}},
                     op::Function, args...) where {T,N}
    return localfilter(T, A, dims, op, args...)
end

function localfilter(::Type{T}, A::AbstractArray{<:Any,N},
                     dims::Union{Colon, Integer, Tuple{Vararg{Integer}},
                                 AbstractVector{<:Integer}},
                     op::Function, args...) where {T,N}
    return localfilter!(similar(Array{T,N}, axes(A)), A, dims, op, args...)
end

@doc @doc(localfilter!) localfilter

# In-place operation.

function localfilter!(A::AbstractArray{T,N},
                      dims::Union{Colon, Integer, Tuple{Vararg{Integer}},
                                  AbstractVector{<:Integer}},
                      op::Function,
                      rngs::Union{IndexInterval, Integer,
                                  AbstractVector{<:IndexInterval},
                                  Tuple{Vararg{IndexInterval}}},
                      args...) where {T,N}
    return localfilter!(A, A, dims, op, rngs, args...)
end

# Wrapper methods when destination is specified.
#
# The filter is applied along all chosen dimensions.  To reduce page memory
# faults, the operation is performed out-of-place (unless dst and A are the
# same) for the first chosen dimension and the operation is performed in-place
# for the other chosen dimensions.

function localfilter!(dst::AbstractArray{T,N},
                      A::AbstractArray{<:Any,N},
                      d::Integer,
                      op::Function,
                      rng::IndexInterval,
                      w::Vector{T} = workspace(T, A, d, rng)) where {T,N}
    return localfilter!(dst, A, Int(d), op, Int(first(rng)), Int(last(rng)), w)
end

function localfilter!(dst::AbstractArray{T,N},
                      A::AbstractArray{<:Any,N},
                      ::Colon,
                      op::Function,
                      rng::IndexInterval,
                      w::Vector{T} = workspace(T, A, :, rng)) where {T,N}
    kmin, kmax = Int(first(rng)), Int(last(rng))
    N ≥ 1 || return copyto!(dst, A)
    localfilter!(dst, A, 1, op, kmin, kmax, w)
    for d in 2:N
        localfilter!(dst, d, op, kmin, kmax, w)
    end
    return dst
end

function localfilter!(dst::AbstractArray{T,N},
                      A::AbstractArray{<:Any,N},
                      ::Colon,
                      op::Function,
                      rngs::Union{AbstractVector{<:IndexInterval},
                                  Tuple{Vararg{IndexInterval}}},
                      w::Vector{T} = workspace(T, A, :, rngs)) where {T,N}
    length(rngs) == N || throw(DimensionMismatch("there must be as many intervals as dimensions"))
    N ≥ 1 || return copyto!(dst, A)
    localfilter!(dst, A, 1, op, rngs[1], w)
    for d in 2:N
        localfilter!(dst, d, op, rngs[d], w)
    end
    return dst
end

function localfilter!(dst::AbstractArray{T,N},
                      A::AbstractArray{<:Any,N},
                      dims::Union{AbstractVector{<:Integer},
                                  Tuple{Vararg{Integer}}},
                      op::Function,
                      rng::IndexInterval,
                      w::Vector{T} = workspace(T, A, dims, rng)) where {T,N}
    m = length(dims)
    m ≥ 1 || return copyto!(dst, A)
    localfilter!(dst, A, dims[1], op, rng, w)
    for d in 2:m
        localfilter!(dst, dims[d], op, rng, w)
    end
    return dst
end

function localfilter!(dst::AbstractArray{T,N},
                      A::AbstractArray{<:Any,N},
                      dims::Union{AbstractVector{<:Integer},
                                  Tuple{Vararg{Integer}}},
                      op::Function,
                      rngs::Union{AbstractVector{<:IndexInterval},
                                  Tuple{Vararg{IndexInterval}}},
                      w::Vector{T} = workspace(T, A, dims, rngs)) where {T,N}
    m = length(dims)
    length(rngs) == m || throw(DimensionMismatch("list of dimensions and list of intervals must have the same length"))
    m ≥ 1 || return copyto!(dst, A)
    localfilter!(dst, A, dims[1], op, rngs[1], w)
    for d in 2:m
        localfilter!(dst, dims[d], op, rngs[d], w)
    end
    return dst
end

# Basic morphological operations.

for (f, op) in ((:erode, min), (:dilate, max))
    fp = Symbol(f, "!")
    @eval begin

        function $f(A::AbstractArray{T,N},
                    dims::Union{Colon, Integer, Tuple{Vararg{Integer}},
                                AbstractVector{<:Integer}},
                    rngs::Union{IndexInterval, Tuple{Vararg{IndexInterval}},
                                AbstractVector{<:IndexInterval}},
                    args...) where {T,N}
            return localfilter(T, A, dims, $op, rngs, args...)
        end

        function $fp(A::AbstractArray{T,N},
                     dims::Union{Colon, Integer, Tuple{Vararg{Integer}},
                                 AbstractVector{<:Integer}},
                     rngs::Union{IndexInterval, Tuple{Vararg{IndexInterval}},
                                 AbstractVector{<:IndexInterval}},
                     args...) where {T,N}
            return localfilter!(A, dims, $op, rngs, args...)
        end

        function $fp(dst::AbstractArray{T,N},
                     A::AbstractArray{T,N},
                     dims::Union{Colon, Integer, Tuple{Vararg{Integer}},
                                 AbstractVector{<:Integer}},
                     rngs::Union{IndexInterval, Tuple{Vararg{IndexInterval}},
                                 AbstractVector{<:IndexInterval}},
                     args...) where {T,N}
            return localfilter!(dst, A, dims, $op, rngs, args...)
        end

    end
end

"""

```julia
workspacelength(n, p)
```

yields the minimal length of the workspace array for applying the van
Herk-Gil-Werman algorithm along a dimension of length `n` with a structuring
element of width `p`.  If `n < 1` or `p ≤ 1`, zero is returned because there is
no needs for a workspace array.

```julia
workspacelength(A, dims, rngs)
```

yields the minimal length of the workspace array for applying the van
Herk-Gil-Werman algorithm along the dimension(s) `dims` of array `A` with a
structuring element defined by the interval(s) `rngs`.  If arguments are not
compatible, ero is returned because there is no needs for a workspace array.

"""
workspacelength(n::Int, p::Int) = (n < 1 || p ≤ 1 ? 0 : n + 2*(p - 1))

function workspacelength(A::AbstractArray{<:Any,N},
                         d::Integer,
                         rng::IndexInterval)::Int where {N}
    1 ≤ d ≤ N || return 0
    n = length(axes(A, d))
    p = length(rng)
    return workspacelength(n, p)
end

function workspacelength(A::AbstractArray{<:Any,N},
                         ::Colon,
                         rng::IndexInterval)::Int where {N}
    N ≥ 1 || return 0
    p = length(rng)
    n = reduce(max, map(length, axes(A)))
    return workspacelength(n, p)
end

function workspacelength(A::AbstractArray{<:Any,N},
                         ::Colon,
                         rngs::Union{AbstractVector{<:IndexInterval},
                                     Tuple{Vararg{IndexInterval}}})::Int where {N}
    len = 0
    if length(rngs) == N
        for d in 1:N
            n = length(axes(A, d))
            p = length(rngs[d])
            len = max(len, workspacelength(n, p))
        end
    end
    return len
end

function workspacelength(A::AbstractArray{<:Any,N},
                         dims::Union{AbstractVector{<:Integer},
                                     Tuple{Vararg{Integer}}},
                         rng::IndexInterval)::Int where {N}
    len = 0
    p = length(rng)
    for dim in dims
        n = length(axes(A, dim))
        len = max(len, workspacelength(n, p))
    end
    return len
end

function workspacelength(A::AbstractArray{<:Any,N},
                         dims::Union{AbstractVector{<:Integer},
                                     Tuple{Vararg{Integer}}},
                         rngs::Union{AbstractVector{<:IndexInterval},
                                     Tuple{Vararg{IndexInterval}}})::Int where {N}
    len = 0
    if (m = length(dims)) == length(rngs)
        for d in 1:m
            n = length(axes(A, dims[d]))
            p = length(rngs[d])
            len = max(len, workspacelength(n, p))
        end
    end
    return len
end

"""

```julia
workspace(T, A, dims, rngs)
```

yields a workspace array for applying the van Herk-Gil-Werman algorithm along
the dimension(s) `dims` of array `A` with a structuring element defined by the
interval(s) `rngs`.  The element type of the workspace is `T` which is that of
`A` by default.

"""
function workspace(::Type{T},
                   A::AbstractArray{<:Any,N},
                   dims::Union{Colon, Integer, AbstractVector{<:Integer},
                               Tuple{Vararg{Integer}}},
                   rngs::Union{IndexInterval, AbstractVector{<:IndexInterval},
                               Tuple{Vararg{IndexInterval}}}) where {T,N}
    return Array{T}(undef, workspacelength(A, dims, rngs))
end

end # module


text:
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 4., Zadania 1.–4.
#

module Interpolation

using PyPlot

FL = Float64

"""
### Checks if the provided vectors are of equal length.

### Arguments
- `x` — first vector
- `y` — second vector

### Return value:
- `Nothing` — the function either throws an error if there is a mismatch in the vector lengths or does nothing.

"""
function _assertDimensionError(x, y)
    if length(x) != length(y)
        throw(DimensionMismatch("provided vectors are of not same length"))
    end
end

"""
### Calculates the difference quotients for a given set of nodes and a set of values of a given function in those nodes.

### Arguments
- `x` — vector of nodes
- `f` — vector of values of a given function in those nodes

### Return value:
- a vector of difference quotients

"""
function ilorazyRoznicowe(x::Vector{FL}, f::Vector{FL})

    # assert vector lengths
    _assertDimensionError(x, f)

    len = length(x)

    # copy the values (the first iteration)
    tmp = copy(f)

    # the output vector that contains all the difference quotients
    output = Vector{FL}()
    push!(output, tmp[1])

    # perform the rest of the iterations
    for p in (len-1):-1:1
        # m = len-p

        # calc all the difference quotients in this column
        for k in 1:p
            a = (tmp[k+1] - tmp[k])
            # m = len-p
            b = (x[k + (len-p)] - x[k])
            tmp[k] = a/b
        end

        push!(output, tmp[1])

    end

    return output

end

"""
### Calculates the value of the interpolation polynomial in a given point.

### Arguments
- `x` — vector of nodes
- `fx` — vector of difference quotients
- `t` — the point in which the function will calculate the value of the interpolation polynomial

### Return value:
- value of the interpolation polynomial in the provided point `t`

"""
function warNewton(x::Vector{FL}, fx::Vector{FL}, t::FL)

    # assert vector length
    _assertDimensionError(x, fx)

    # here `n` is actually equal to $n+1$ relative to the algorithm explanation in the attached document
    n = length(x)

    # implement the recursive function
    # special case wₙ(z) = f[x₀,…,xₙ] if k == n
    # otherwise return normal recursive value
    w = (z,k) -> (k == n ?
                  fx[n]
                : fx[k] + (z - x[k]) * w(z,k+1))

    # in Julia we count starting at `1`, so relative to the theory behind this algorithm it’s w₀(t) = Nₙ(t)
    return w(t, 1)

end

"""
### Calculates all the coefficients of the Newton’s interpolation polynomial in its natural polynomial form.

### Arguments
- `x` — vector of nodes
- `fx` — vector of difference quotients

### Return value:
- vector of coefficients of the Newton’s interpolation polynomial

"""
function naturalna(x::Vector{FL}, fx::Vector{FL})

    # assert vector length
    _assertDimensionError(x, fx)

    n = length(x)

    coefficients = [fx[n]]

    for k in (n-1):-1:1

        xₖ = x[k]

        new_b₀ = -coefficients[1] * xₖ + fx[k]

        # when calculating the new coefficients the `m`-th coefficient is copied over to the new `m+1`-th coefficient,
        # whilst the last coefficient is calculated using the previous last coefficient and the difference quotient
        new_coefficients = vcat([new_b₀], coefficients)

        # we iterate over all coefficients except the first one and the last one
        for m in n-k:-1:2

            # b^{(k)}_m = b^{(k+1)}_{m-1} - b^{(k+1)}_m ⋅ xₖ
            new_coefficients[m] = coefficients[m-1] - coefficients[m] * xₖ

        end

        coefficients = new_coefficients

    end

    return coefficients

end

"""
### Draws a graph of a given function and its interpolation.

### Arguments
- `f` — function to examine
- `a` — beginning of the interval to examine
- `b` — ending of the interval to examine
- `n` — the “resolution” (the degree of the interpolation polynomial) of the interpolation

### Return value:
- `Nothing` — the function invokes PyPlot plotting function.

"""
function rysujNnfx(f, a::FL, b::FL, n::Int)

    h = (b-a)/n
    res = h/10
    x = vec([i for i in a:h:b])

    fx = ilorazyRoznicowe(x, vec([f(i) for i in x]))

    XA = a:res:b

    YA_f = [f(i) for i in XA]
    YA_N = [warNewton(x, fx, i) for i in XA]

    plot(XA, YA_N)
    plot(XA, YA_f)

end

end

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
function test_gaussian_kl_loss_layer(backend::Backend, T, eps)
  println("-- Testing GaussianKLLossLayer on $(typeof(backend)){$T}...")

  ############################################################
  # Prepare Data for Testing
  ############################################################
  tensor_dim = rand(2:5)
  dims = tuple(rand(6:11, tensor_dim)...)
  println("    > $dims")
  mus = rand(T, dims)
  sigmas = sqrt.(rand(T, dims).^2)

  ############################################################
  # Setup
  ############################################################
  weight = T(1.1)
  layer  = GaussianKLLossLayer(; bottoms=[:predictions, :labels], weight=weight)
  mu_blob  = make_blob(backend, T, dims)
  sigma_blob = make_blob(backend, T, dims)

  mu_diff_blob  = make_blob(backend, T, dims)
  sigma_diff_blob  = make_blob(backend, T, dims)

  copy!(mu_blob, mus)
  copy!(sigma_blob, sigmas)
  inputs = Blob[mu_blob, sigma_blob]
  diffs = Blob[mu_diff_blob, sigma_diff_blob]

  state  = setup(backend, layer, inputs, diffs)

  forward(backend, state, inputs)

  n = length(mu_blob)
  loss = 0.5(sum(mus.^2 + sigmas.^2 - 2log.(sigmas)) - n)
  loss *= weight/get_num(mu_blob)
  @test -eps < loss-state.loss < eps

  backward(backend, state, inputs, diffs)
  grad = mus
  grad *= weight/get_num(mu_blob)
  diff = similar(grad)
  copy!(diff, diffs[1])
  @test all(-eps .< grad - diff .< eps)

  grad = sigmas - 1./sigmas
  grad *= weight/get_num(mu_blob)
  diff = similar(grad)
  copy!(diff, diffs[2])
  @test all(-eps .< grad - diff .< eps)

  shutdown(backend, state)
end

function test_gaussian_kl_loss_layer(backend::Backend)
  test_gaussian_kl_loss_layer(backend, Float32, 5e-2)
  test_gaussian_kl_loss_layer(backend, Float64, 1e-8)
end

if test_gpu
  test_gaussian_kl_loss_layer(backend_gpu)
end

if test_cpu
  test_gaussian_kl_loss_layer(backend_cpu)
end


text:
using AntsArticle
using Test

@testset "AntsArticle.jl" begin
    @test true
end


text:
module QuAlgorithmZoo

using LuxurySparse, LinearAlgebra
using Yao
using Yao.Intrinsics
using Yao.Registers
using Yao.Blocks
import Yao.Blocks: mat, dispatch!, niparameters, iparameters, setiparameters!, cache_key, print_block, _make_rot_mat, apply!, PrimitiveBlock
import Base: ==, copy, hash
import Yao.Intrinsics: ishermitian, isreflexive, isunitary

export openbox

"""
    openbox(block::AbstractBlock) -> AbstractBlock

For a black box, like QFTBlock, you can get its white box (loyal simulation) using this function.
"""
function openbox end

include("Miscellaneous.jl")
include("Adam.jl")
include("QFT.jl")
include("CircuitBuild.jl")
include("QCOptProblem.jl")
include("RotBasis.jl")
include("Grover.jl")
include("PhaseEstimation.jl")
include("HHL.jl")
include("hamiltonian_solvers.jl")
include("HadamardTest.jl")

end # module