text:
using Compat

function advance_trace(s1, s2, U1interp, U2interp, dt)
    u1 = U1interp.(s1,s2)
    u2 = U2interp.(s1,s2)
    @. s1 += dt*u1
    @. s2 += dt*u2
end

function streakplot(;X=[], Y=[], interior=[], dataname=nothing, ilist=[],
                    startx=[],
                    starty=[])
    x = X[1:end,1]
    y = Y[1,1:end]
    traces = length(startx)
    s1 = Array{Array{Float64}}(undef, traces)
    s2 = Array{Array{Float64}}(undef, traces)
    for k=1:traces
        s1[k] = [startx[k]]
        s2[k] = [starty[k]]
    end
    # Iterate
    t = load(dataname(ilist[1]))["ti"] # Get initial time
    stepnum = 0
    for i=ilist
        stepnum += 1
        datafile = dataname(i)
        info("Loading ", datafile)
        U1,U2 = [],[]
        dt = 0.0
        @time begin
            U1 = h5read(datafile, "U1")
            U2 = h5read(datafile, "U2")
            dt = h5read(datafile, "ti") - t
        end
        # Set free stream speed to get traces coming in and out
        @. U1[ (U1==0.0) & (U2==0.0) ] = 1.0
        U1interp = LinearInterpolation((x,y), U1, extrapolation_bc=0)
        U2interp = LinearInterpolation((x,y), U2, extrapolation_bc=0)
        t += dt
        for k=1:traces
            advance_trace(s1[k], s2[k], U1interp, U2interp, dt)
        end
        # Add points, fill in gaps and trim
        if mod(stepnum, 10) == 0
            for k=1:traces
                # Add points to beginning
                unshift!(s1[k], startx[k])
                unshift!(s2[k], starty[k])
                # Interpolate into line
                kx, ky = s1[k], s2[k]
                m = 1
                while m <= length(kx)-1
                    dx = kx[m+1]-kx[m]
                    dy = ky[m+1]-ky[m]
                    if sqrt(dx^2+dy^2) > 0.1
                        newx = kx[m]+dx/2
                        newy = ky[m]+dy/2
                        splice!(kx, m+1:m, newx)
                        splice!(ky, m+1:m, newy)
                    else
                        m += 1
                    end
                end
                xmax = maximum(x)-0.1
                inside = @. kx < xmax
                s1[k] = kx[inside]
                s2[k] = ky[inside]
            end
        end
    end
    # Plot
    begin
        marker = zeros(size(X))
        marker[interior] = 1.0
        markerInterp = LinearInterpolation((x,y), marker, extrapolation_bc=0)
        for k=1:traces
            kx, ky = s1[k], s2[k]
            #inside = @. kx < 5
            inside = markerInterp.(kx, ky) .> 0
            plot(kx[inside], ky[inside], "-k")
        end
    end
end


text:
import AmplNLWriter

@testset "[nl_convert] check special conversion cases" begin
    special_cases = [:cbrt, :abs2, :inv, :log2, :log1p, :exp2, :expm1, :sec,
                     :csc, :cot, :sind, :cosd, :tand, :asind, :acosd, :atand,
                     :secd, :cscd, :cotd, :sech, :csch, :coth, :asech, :acsch]
    for func in special_cases
        x = rand()
        expr = Expr(:call, func, x)
        @test isapprox(eval(AmplNLWriter.convert_formula(expr)), eval(expr), atol=1e-6)
    end
    # These functions need input >1
    for func in [:acoth, :asec, :acsc, :acot, :asecd, :acscd, :acotd]
        x = rand() + 1
        expr = Expr(:call, func, x)
        @test isapprox(eval(AmplNLWriter.convert_formula(expr)),
               eval(expr), atol=1e-6)
    end
end

@testset "[nl_convert] check numeric values" begin
    x = rand()
    @test AmplNLWriter.convert_formula(:($x)) == :($x)
    x = -rand()
    @test AmplNLWriter.convert_formula(:($x)) == :($x)
end

@testset "[nl_convert] check unary, binary and n-ary plus" begin
    expr = :(+(1))
    @test AmplNLWriter.convert_formula(expr) == :(1)
    expr = :(1 + 2)
    @test AmplNLWriter.convert_formula(expr) == :(1 + 2)
    expr = :(1 + 2 + 3)
    @test AmplNLWriter.convert_formula(expr) == :(sum(1, 2, 3))
end

@testset "[nl_convert] check unary, binary and n-ary minus" begin
    expr = :(- x)
    @test AmplNLWriter.convert_formula(expr) == :(neg(x))
    expr = :(x - y)
    @test AmplNLWriter.convert_formula(expr) == :(x - y)
    expr = :(x - y - z)
    @test AmplNLWriter.convert_formula(expr) == :((x - y) - z)
end

@testset "[nl_convert] check n-ary multiplication" begin
    expr = :(x * y * z)
    @test AmplNLWriter.convert_formula(expr) == :(x * (y * z))
end

@testset "[nl_convert] check comparison expansion" begin
    expr = :(1 < 2 < 3)
    @test AmplNLWriter.convert_formula(expr) == :(1 < 2 && 2 < 3)
    expr = :(1 < 2 < 3 < 4)
    @test AmplNLWriter.convert_formula(expr) == :((1 < 2 && 2 < 3) && 3 < 4)
end

@testset "[nl_convert] check user defined functions error" begin
    m = Model(solver=AmplNLSolver("any_solver"))

    myf(x,y) = (x-1)^2+(y-2)^2
    JuMP.register(m, :myf, 2, myf, autodiff=true)

    @variable(m, x[1:2] >= 0.5)
    @NLobjective(m, Min, myf(x[1], x[2]))

    @test_throws Exception solve(m)

    # clean up temp file
    for file in readdir(AmplNLWriter.solverdata_dir)
        if file != ".gitignore"
            rm(joinpath(AmplNLWriter.solverdata_dir, file))
        end
    end
end

ext: jl


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

name = "CouenneBuilder"
version = v"0.1.0"

# Collection of sources required to build CouenneBuilder
sources = [
    "https://github.com/ampl/coin/archive/v20170826.tar.gz" =>
    "b3c7395c22552614107161e1672c232a4fce4effca409863814a7fe623997651",

    # BLAS for Bonmin build
    "https://www.coin-or.org/download/source/BuildTools/BuildTools-Blas-1.0.7.tgz" =>
    "4fe4eaa46f82ef3da125401e312dd5a2bc18b38f5b09e547775e8edadd1826ba",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/coin-20170826/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DBUILD_SHARED_LIBS=ON .
make
make install

# cd Bonmin/
# ./configure --prefix=$prefix --host=$target

mkdir -p $prefix/bin
install bin/* $prefix/bin/

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    FreeBSD(:x86_64),
    Linux(:x86_64, :glibc),
    MacOS(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libCbcSolver", :libCbcSolver),
    LibraryProduct(prefix, "libipopt", :libipopt),
    LibraryProduct(prefix, "libOsi", :libOsi),
    ExecutableProduct(prefix, "cbc", :cbc),
    LibraryProduct(prefix, "libCoinUtils", :libCoinUtils),
    LibraryProduct(prefix, "libCgl", :libCgl),
    LibraryProduct(prefix, "libOsiCbc", :libOsiCbc),
    LibraryProduct(prefix, "libOsiCommonTests", :libOsiCommonTests),
    ExecutableProduct(prefix, "couenne", :couenne),
    LibraryProduct(prefix, "libCouenne", :libCouenne),
    LibraryProduct(prefix, "libOsiClp", :libOsiClp),
    LibraryProduct(prefix, "libbonmin", :libbonmin),
    LibraryProduct(prefix, "libCbc", :libCbc),
    LibraryProduct(prefix, "libClp", :libClp)
]

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

]

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

max_stars_repo_name: jgoldfar/CouenneBuilder

max_issues_repo_name: jgoldfar/CouenneBuilder

max_forks_repo_name: jgoldfar/CouenneBuilder


text:
using JSON

function get_precomputed()
    inlayer = "gs://neuroglancer/pinky40_v11/watershed_mst_trimmed_sem_remap"
    outlayer = "gs://neuroglancer/pinky40_v11/analysis/contact_dicts/"
    p = ContactAnalysis.Precomputed.PrecomputedWrapper(inlayer, 1)
    return inlayer, outlayer, p
end

function schedule_tasks(ranges, inlayer, outlayer)
    for r in ranges
        payload_info = ContactAnalysis.CountEdgesTask.CountEdgesPayloadInfo(inlayer,
                ContactAnalysis.CountEdgesTask.slices_to_str([r...]), 1, outlayer)
        basic_info = SimpleTasks.Tasks.BasicTask.Info(0, "CountEdgesTask", inlayer, [""])
        ContactAnalysis.AWSScheduler.schedule_count_edges(basic_info, payload_info, queue_name="task-queue-pinky")
    end
end

function schedule_pinky40()
    inlayer, outlayer, p = get_precomputed()
    ranges = make_ranges(p)
    schedule_tasks(ranges, inlayer, outlayer)
end

function map_chunks(sz, o, chunk, overlap)
    """
    Create ranges for 3D chunks

    Args:
        sz: dataset size (collection of 3 ints)
        o: origin voxel (collection of 3 ints)
        chunk: size of the chunk (collection of 3 ints)
        overlap: number of voxels in positive overlap (collection of 3 ints)

    Returns:
        list of ranges describing chunks in larger volume
    """
    # pre = PrecomputedWrapper(segmentation_storage, scale_idx)
    # sz = pre.val[:_scale]["size"]
    # o = pre.val[:_scale]["voxel_offset"]
    ranges = []
    for x_start in o[1]:max(chunk[1],1):o[1]+sz[1]-1
        x_end = min(x_start+chunk[1]+overlap[1]-1, o[1]+sz[1]-1)
        for y_start in o[2]:max(chunk[2],1):o[2]+sz[2]-1
            y_end = min(y_start+chunk[2]+overlap[2]-1, o[2]+sz[2]-1)
            for z_start in o[3]:max(chunk[3],1):o[3]+sz[3]-1
                z_end = min(z_start+chunk[3]+overlap[3]-1, o[3]+sz[3]-1)
                r = x_start:x_end, y_start:y_end, z_start:z_end
                push!(ranges, r)
            end
        end
    end
    return ranges
end

"""
Given a precomputed PyCall object, chunk size, and overlap, give list of ranges
"""
function make_ranges(p, chunk_size=[256,256,64], overlap=[1,1,1])
    return map_chunks(p.val[:_scale]["size"], p.val[:_scale]["voxel_offset"],
                                                        chunk_size, overlap)
end

"""
Pull all contact analysis files from gcloud
"""
function copy_results_locally()
    Base.run(`gsutil -m cp gs://neuroglancer/pinky40_v11/analysis/contact_dicts/* ~/data/contact_dicts/`)
end

"""
For a directory of contact analysis json dicts, merge them all together
"""
function merge_json_dicts(dir)
    contacts = Dict{String,Float64}()
    for fn in readdir(dir)
        d = JSON.parsefile(joinpath(dir, fn), dicttype=Dict{String,Float64}, use_mmap=true)
        merge!(contacts, d)
    end
    int_contacts = Dict{Tuple{Int64,Int64},Float64}()
    for (k, v) in contacts
        ki = (map(Int, map(parse, split(k[2:end-1], r",")))...)
        int_contacts[ki] = v
    end
    return contacts
end

function str_to_slice_start_stop(s)
    m=match(r"(\d+)-(\d+)_(\d+)-(\d+)_(\d+)-(\d+)",s)
    bounds = map(x->parse(Int,x), m.captures)
    return [bounds[1], bounds[3], bounds[5], bounds[2], bounds[4], bounds[6]]
end

function str_to_slices(s)
    m=match(r"(\d+)-(\d+)_(\d+)-(\d+)_(\d+)-(\d+)",s)
    bounds = map(x->parse(Int,x), m.captures)
    return (bounds[1]:bounds[2], bounds[3]:bounds[4], bounds[5]:bounds[6])
end

function slices_to_str(s)
    t=map(slice_to_str,s)
    return "$(t[1])_$(t[2])_$(t[3])"
end

function slice_to_str(s)
    return "$(s.start)-$(s.stop)"
end

function get_folder(name)
    return joinpath(homedir(), "data", name)
end

"""
Compute center on output of str_to_slice_start_stop
"""
 function get_center(s)
    ctr = round(Int, [(s[1]+s[4])/2, (s[2]+s[5])/2, (s[3]+s[6])/2])
end

"""
Horizontally concatenate center and bbox info to an edge list
"""
 function append_locations(fn, edge_list)
    s = str_to_slice_start_stop(fn)
    ctr = repmat(get_center(s), 1, size(edge_list,1))'
    bbox = repmat(s, 1, size(edge_list,1))'
    return hcat(edge_list[:,1:2], ctr, edge_list[:,3], bbox)
end

"""
Create contact edge list
"""
 function create_contact_edge_list(fn)
    d = JSON.parsefile(fn, dicttype=Dict{AbstractString,Float64}, use_mmap=true)
    contacts = Array{Any,2}(length(d), 3)
    segs = Set{Int64}()
    for (i, (k, v)) in enumerate(d)
        ki = (map(Int, map(parse, split(k[2:end-1], r",")))...)
        contacts[i,:] = [ki..., v]
        push!(segs, ki[1])
        push!(segs, ki[2])
    end
    return contacts, segs
end

"""
Convert contact dicts into contact edge list and list of segments for proximity
"""
function resave_contact_dict(inputs)
    i, n, dir, fn = inputs
    print("\r$i / $n")
    contacts, segs = create_contact_edge_list(joinpath(dir, fn))
    contacts = append_locations(fn, contacts)
    contacts_dst_fn = joinpath(get_folder("touches"), fn[1:end-17])
    segs_dst_fn = joinpath(get_folder("segs"), fn[1:end-17])
    writedlm(contacts_dst_fn, contacts, ',')
    writedlm(segs_dst_fn, segs, ',')
end

"""
Write out contact lists & sets of segments
"""
function dir_to_contact_lists()
    dir = get_folder("contact_dicts")
    files = readdir(dir)
    inputs = [(i, length(files), dir, fn) for (i,fn) in enumerate(files)]
    pmap(resave_contact_dict, inputs)
end

"""
From chunk with chunk_range, find range of neighbor in relative direction
"""
function get_neighbor(chunk_range, chunk_size, direction=[1,0,0])
    return (chunk_range[1] + chunk_size[1]*direction[1],
            chunk_range[2] + chunk_size[2]*direction[2],
            chunk_range[3] + chunk_size[3]*direction[3])
end

function get_start(chunk_range)
    return [chunk_range[1].start, chunk_range[2].start, chunk_range[3].start]
end

function get_end(chunk_range)
    return [chunk_range[1].stop, chunk_range[2].stop, chunk_range[3].stop]
end

"""
Given two grided coords, order them
"""
function order_coords(coord1, coord2)
    return reduce(|, coord1 .<= coord2) ? (coord1, coord2) : (coord2, coord1);
end

function write_proximities(inputs)
    i, n, directions, chunk, chunk_r = inputs
    print("\r$i / $n")
    segs_fn = joinpath(get_folder("segs"), slices_to_str(chunk_r))
    if isfile(segs_fn)
        chunk_start = get_start(chunk_r)
        chunk_stop = get_end(chunk_r)
        for d in directions
            neigh_r = get_neighbor(chunk_r, chunk, d)
            neigh_fn = joinpath(get_folder("segs"), slices_to_str(neigh_r))
            if isfile(neigh_fn)
                neigh_start = get_start(neigh_r)
                neigh_stop = get_end(neigh_r)
                cube_start, _ = order_coords(chunk_start, neigh_start)
                _, cube_stop = order_coords(chunk_stop, neigh_stop)
                ctr = get_center([cube_start..., cube_stop...])
                try
                    segs = readdlm(segs_fn, Int64)
                    combined_segs = vcat(segs, readdlm(neigh_fn, Int64))
                    unique_segs = sort(unique(segs))
                    if unique_segs[1] == 0
                        unique_segs = unique_segs[2:end]
                    end
                    n = length(unique_segs)
                    k = 1
                    proximity = Array{Any,2}(Int((n^2-n)/2), 9)
                    for (i, seg1) in enumerate(unique_segs[1:end-1])
                        for seg2 in unique_segs[i+1:end]
                            proximity[k,:] = [seg1, seg2,
                                        round(sum(chunk.*[8,8,40]/1000.*d),2),
                                        cube_start..., cube_stop...]
                            k += 1
                        end
                    end
                    cube_range = (cube_start[1]:cube_stop[1],
                                    cube_start[2]:cube_stop[2],
                                    cube_start[3]:cube_stop[3])
                    r_name = slices_to_str(cube_range)
                    proximity_fn = joinpath(get_folder("proximity"), r_name)
                    # proximity = append_locations(proximity_fn, proximity)
                    writedlm(proximity_fn, proximity, ',')
                catch e
                    println((e, chunk_start, neigh_start))
                end
            end
        end
    end
end

"""
Compile list of all segment pairs in nearest-neighbor chunks
"""
function compute_all_proximities(sz, origin, chunk=[256,256,64], overlap=[1,1,1])
    ranges = map_chunks(sz, origin, chunk, overlap)
    # directions = ([0,0,0])
    directions = ([ 1, 1, 1],
                  [ 1, 0, 1],
                  [ 0,-1, 1],
                  [ 0, 1, 1],
                  [ 0, 0, 1],
                  [ 0,-1, 1],
                  [-1, 1, 1],
                  [-1, 0, 1],
                  [-1,-1, 1],
                  [ 1, 1, 0],
                  [ 1, 0, 0],
                  [ 0, 1, 0],
                  [ 0, 0, 0])
    inputs = [(i, length(ranges), directions, chunk, chunk_r)
                            for (i,chunk_r) in enumerate(ranges)]
    pmap(write_proximities, inputs)
end

"""
Add index in first column of the edge list
"""
function index_edges(edge_list)
    id = 1:size(edge_list,1)
    return hcat(id, edge_list)
end

"""
Combine directory of edge lists into one large edge list
"""
function compile_edge_lists(name="touches")
    dir = get_folder(name)
    files = readdir(dir)
    edge_list = append_locations(files[1], readdlm(joinpath(dir,files[1]), Any))
    for (i, fn) in enumerate(files[2:end])
        print("\r$i / $(length(files))")
        try
            next_edge_list = append_locations(fn, readdlm(joinpath(dir,fn), Any))
            edge_list = vcat(edge_list, next_edge_list)
        end
    end
    edge_fn = joinpath(homedir(), "data", string(name, ".csv"))
    writedlm(edge_fn, index_edges(edge_list), ',')
end


text:
using Documenter, ComputationalBiologyCourse

makedocs(
    sitename = "Computational Biology Course",
    pages = [
        "Home" => "index.md",
        "Syllabus" => "Syllabus.md",
        "Lesson 1 - Getting Started" => [
            "Lesson1/index.md",
            "Windows Users" => "Lesson1/0_windows.md",
            "Using the Terminal" => "Lesson1/1_terminal.md",
            "Git for Version Control" => "Lesson1/2_git.md"
            ]
    ],
    authors = "Kevin Bonham, PhD",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true")
)

const repo="github.com/wellesley-bisc195/computing_bio_course.git"
const PR = get(ENV, "TRAVIS_PULL_REQUEST", "false")
if PR == "false"
    # Normal case, only deply docs if merging to master or release tagged
    deploydocs(repo=repo)
else
    @info "Deploying review docs for PR #$PR"
    # TODO: remove most of this once https://github.com/JuliaDocs/Documenter.jl/issues/1131 is resolved

    # Overwrite Documenter's function for generating the versions.js file
    foreach(Base.delete_method, methods(Documenter.Writers.HTMLWriter.generate_version_file))
    Documenter.Writers.HTMLWriter.generate_version_file(_, _) = nothing
    # Overwrite necessary environment variables to trick Documenter to deploy
    ENV["TRAVIS_PULL_REQUEST"] = "false"
    ENV["TRAVIS_BRANCH"] = "master"

    deploydocs(
        devurl="preview-PR$(PR)",
        repo=repo,
    )
end


text:
"""
    generate_tokenizer_from_sed(sed_script_path)

This returns Julia code, that is the translation of a simple sed script
for tokenizing.
This doesn't fully cover all the functionality of sed,
but it covers enough for many purposes.
"""
function generate_tokenizer_from_sed(sed_script, extended=false)::Expr
    code = quote
        ss = input
    end
    for src_line in eachline(sed_script)
        src_line=strip(src_line)
        length(src_line)==0 && continue # skip blanks
        src_line[1]=='#' && continue # skip comments

        #sed lines are `<op><sep><pattern><sep><replacement><sep>flags`
        seperator = src_line[2]
        op, pattern, replacement, flags = split(src_line, seperator)
        @assert(op=="s") #substitute
        @assert(flags=="g" || flags=="", "Unsupported flags: $flags") #substitute

        if extended
        else
            # Normal sed uses `\(` instead of `(` for grouping
            pattern=replace(pattern, raw"\(" => "(")
            pattern=replace(pattern, raw"\)" => ")")

            #sed accepts `&` as whole match
            replacement=replace(replacement, "&" => raw"\0")
        end

        push!(code.args, :(
            ss=replace(ss,
                       $(Regex(pattern)) =>
                       $(Base.SubstitutionString(replacement)))
        ))
    end
    push!(code.args, :(split(ss)))
    code
end

let script = joinpath(@__DIR__, "penn.sed")
    @eval function penn_tokenize(input::AbstractString)
        $(generate_tokenizer_from_sed(script, false))
    end
end

"""
    penn_tokenize(input::AbstractString)

"... to produce Penn Treebank tokenization on arbitrary raw text.
Yeah, sure" quote Robert MacIntyre

Tokenization does a number of things like seperate out contractions:
"shouldn't" becomes ["should", "n't"]
Most other punctuation becomes &'s.
Exception is periods which are not touched.
The input should be a single sentence;
but it will likely be relatively fine if it isn't.
Depends exactly what you want it for.

This is a direct (automatic) translation of the original sed script.

If you want to mess with exactly what it does it is actually really easy.
copy the penn.sed file from this repo, modify it to your hearts content.
There are some lines you can uncomment out.
You can generate a new tokenizer using:

```
@generated function custom_tokenizer(input::AbstractString)
    generate_tokenizer_from_sed(joinpath(@__DIR__, "custom.sed"))
end
```
"""
penn_tokenize

let script = joinpath(@__DIR__, "improved_penn.sed")
    @eval function improved_penn_tokenize(input::AbstractString)
        $(generate_tokenizer_from_sed(script, true))
    end
end

"""
    improved_penn_tokenize(input::AbstractString)

The Improved Penn Treebank tokenizer.
This is a port of NLTK's modified Penn Tokenizer.
It has a bundle of minor changes,
that I don't think are actually documented anywhere.
But things like `cannot cannot` become `can not can not`
where as the original would produce `can not cannot`.

The tokenizer still seperates out contractions:
"shouldn't" becomes ["should", "n't"]

The input should be a single sentence;
but again it will likely be relatively fine if it isn't.
Depends exactly what you want it for.

This matches NLTK's `nltk.tokenize.TreeBankWordTokenizer.tokenize`
"""
improved_penn_tokenize


text:
"""
    read_drifters(pth,lst;chnk=Inf,rng=(missing,missing))

Read near-surface [drifter data](https://doi.org/10.1002/2016JC011716) from the
[Global Drifter Program](https://doi.org/10.25921/7ntx-z961) into a DataFrame.

Note: need to use NetCDF.jl as NCDatasets.jl errors when TIME = Inf
```
pth,list=drifter_files()
df=read_drifters( pth*lst[end],chnk=1000,rng=(2014.1,2014.2) )
#sort!(df, [:t, :lat])
#CSV.write(pth*"Drifter_hourly_2005_2019.csv", df)
```
"""
function read_drifters(fil::String;chnk=Inf,rng=(missing,missing))
   t=ncread(fil,"TIME")
   t_u=ncgetatt(fil,"TIME","units")
   lo=ncread(fil,"LON")
   la=ncread(fil,"LAT")
   ID=ncread(fil,"ID")

   ii=findall(isfinite.(lo.*la.*t))
   t_ii=t[ii]
   t_ii=timedecode.(t_ii, t_u)
   tmp=dayofyear.(t_ii)+(hour.(t_ii) + minute.(t_ii)/60 ) /24
   t_ii=year.(t_ii)+tmp./daysinyear.(t_ii)
   isfinite(rng[1]) ? jj=findall( (t_ii.>rng[1]).&(t_ii.<=rng[2])) : jj=1:length(ii)

   ii=ii[jj]
   t=t_ii[jj]
   lo=lo[ii]
   la=la[ii]
   ID=ID[ii]

   df = DataFrame(ID=Int[], lon=Float64[], lat=Float64[], t=Float64[])

   !isinf(chnk) ? nn=Int(ceil(length(ii)/chnk)) : nn=1
   for jj=1:nn
      #println([jj nn])
      !isinf(chnk) ? i=(jj-1)*chnk.+(1:chnk) : i=(1:length(ii))
      i=i[findall(i.<length(ii))]
      append!(df,DataFrame(lon=lo[i], lat=la[i], t=t[i], ID=Int.(ID[i])))
   end

   return df
end

"""
    read_drifters( pth, lst )

Read near-surface [hourly drifter data](https://doi.org/10.1002/2016JC011716) from the
[Global Drifter Program](https://doi.org/10.25921/7ntx-z961) into a DataFrame.

Note: need to use NetCDF.jl as NCDatasets.jl errors when TIME = Inf

```
pth,list=drifter_files()
df=read_drifters( pth, lst)
```
"""
function read_drifters( pth, lst )
   df = DataFrame([fill(Int, 1) ; fill(Float64, 3)], [:ID, :lon, :lat, :t])
   for fil in lst
      println(fil)
      append!(df,read_drifters( pth*fil,chnk=10000,rng=(2005.0,2020.0) ))
      #println(size(df))
   end
   return df
end

function drifter_files()
   pth="Drifter_hourly_v013/"
   lst=["driftertrajGPS_1.03.nc","driftertrajWMLE_1.02_block1.nc","driftertrajWMLE_1.02_block2.nc",
      "driftertrajWMLE_1.02_block3.nc","driftertrajWMLE_1.02_block4.nc","driftertrajWMLE_1.02_block5.nc",
      "driftertrajWMLE_1.02_block6.nc","driftertrajWMLE_1.03_block7.nc"]
   return pth,list
end


text:
module PhaseReduction

export AbstractDiffSys, SpecialDiffSys, DiffSys, CoupledDiffSys, integrate,
    PeriodicOrbit, PhaseSensitivity, PhaseCoupling,
    vec_field, jacobian,
    vanderPol, FitzHughNagumo, HodgkinHuxley, FastSpikingNeuron

using LinearAlgebra

include("utils.jl")
include("reduction.jl")
include("diffeqs/vanderPol.jl")
include("diffeqs/FitzHughNagumo.jl")
include("diffeqs/HodgkinHuxley.jl")
include("diffeqs/FastSpikingNeuron.jl")

end


text:
module Makeitso

#export Target
#export make
export @target
export @make
#export @update!

using JLD2
using FileIO

STORE_DIR = "store"
setstore(s) = (global STORE_DIR=s)

mutable struct Target
    deps::Vector{Target}
    recipe
    timestamp
    cache
    name
    hash
end

function Target(name, recipe, deps::Vector{Target}, hash)
    t = Target(deps, recipe, 0.0, nothing, name, hash)
end

function make(target, level=0)
    for t in target.deps
        make(t, level+1)
    end

    varname = String(target.name)
    filename = String(target.name) * ".jld2"
    filename = joinpath(STORE_DIR, filename)

    # No file means this is the first run ever
    if !isfile(filename)
        update!(target)
        @info "level $level dep $(target.name): computed from dependencies [initial computation]."
        return target.cache
    end

    # Target was made in previous session. Is it up-to-date?
    deps_stamp = reduce(max, getfield.(target.deps, :timestamp), init=0.0)
    if target.cache == nothing
        d = load(filename)
        if d["hash"] != target.hash
            update!(target)
            @info "level $level dep $(target.name): computed from dependencies [recipe modified]."
        elseif d["timestamp"] < deps_stamp
            update!(target)
            @info "level $level dep $(target.name): computed from dependencies [store out-of-date]."
        else
            target.timestamp = d["timestamp"]
            target.cache = d[varname]
            @info "level $level dep $(target.name): restored from disk."
        end
        return target.cache
    end

    # Target was computed in this session. Is it up-to-date?
    if target.timestamp < deps_stamp
        update!(target)
        @info "level $level dep $(target.name): computed from dependencies [memory cache out-of-date]."
        return target.cache
    else
        @info "level $level dep $(target.name): retrieved from memory cache."
        return target.cache
    end
end

function update!(target::Target)

    varname = String(target.name)
    filename = String(target.name) * ".jld2"
    filename = joinpath(STORE_DIR, filename)
    mkpath(STORE_DIR)

    target.cache = target.recipe(getfield.(target.deps, :cache)...)
    target.timestamp = time()
    save(filename, Dict(
        varname => target.cache,
        "timestamp" => target.timestamp,
        "hash" => target.hash))
end

# function update!(target::Target, val)
#
#     target.cache = val
#     target.timestamp = time()
#
#     varname = String(target.name)
#     filename = String(target.name) * ".jld2"
#
#     save(filename, Dict(
#         varname => target.cache,
#         "timestamp" => target.timestamp,
#         "hash" => target.hash))
# end

# Position independent hashing of expressions
pihash(x) = pihash(x, zero(UInt))
pihash(x::Expr,h) = pihash(x.args, pihash(x.head, h))
pihash(x::LineNumberNode,h) = h
function pihash(x::Array,h)
    for y in x
        h = pihash(y, h)
    end
    h
end
pihash(x::Any,h) = hash(x,h)

# macro update!(var)
#     :(update!($(esc(Symbol("target_",var)))))
# end
#
# macro update!(var, val)
#     :(update!($(esc(Symbol("target_",var))), $(esc(val))))
# end

macro target(out, recipe)

    @assert out isa Symbol
    @assert recipe.head == :->

    out_target = Symbol("target_", out)
    file_name = String(out) * ".jld2"
    file_name = joinpath(STORE_DIR, file_name)

    vnames = [] # A, B, C
    tnames = [] # target, target_B, target_C

    tp = recipe.args[1]
    if tp isa Symbol
        push!(vnames, tp)
        push!(tnames, esc(Symbol("target_", tp)))
    else
        @assert tp.head == :tuple
        for arg in tp.args
            push!(vnames, arg)
            push!(tnames, esc(Symbol("target_", arg)))
        end
    end

    exists = isdefined(__module__, out_target)
    recipe_hash = pihash(recipe)
    if exists
        xp = quote
            if $recipe_hash != $(esc(out_target)).hash
                $(esc(out_target)).deps = Target[$(tnames...)]
                $(esc(out_target)).recipe = $(esc(recipe))
                $(esc(out_target)).timestamp = 0.0
                $(esc(out_target)).cache = nothing
                $(esc(out_target)).hash = $recipe_hash
                isfile($file_name) && rm($file_name)
            end
        end
    else
        xp = :($(esc(out_target)) = Target($(QuoteNode(out)), $(esc(recipe)), Target[$(tnames...)], $recipe_hash))
    end
    return xp
end

macro make(vname)
    tname = Symbol("target_", vname)
    :(make($(esc(tname))))
end

end # module


text:
function solve_a(instructions)
    steps = 0

    index = 1
    len = length(instructions)
    while index < len + 1
        # println("before: $(index)")
        offset = instructions[index]
        instructions[index] += 1
        index += offset
        steps += 1
        # println("after: $(index), $(offset), ")
    end
    steps
end

function solve_b(instructions)
    steps = 0

    index = 1
    len = length(instructions)
    while index < len + 1
        # println("before: $(index)")
        offset = instructions[index]

        if offset >= 3
            instructions[index] -= 1
        else
            instructions[index] += 1
        end

        index += offset
        steps += 1
        # println("after: $(index), $(offset), ")
    end
    steps
end

data = parse.(Int, readlines("input.txt"));
sol_a = solve_a(copy(data))
println("sol a: $(sol_a)")

sol_b = solve_b(copy(data))
println("sol b: $(sol_b)")


text:
module BattleSim

include("warriors.jl")

end

max_stars_repo_name: ordovician/BattleSim

max_issues_repo_name: ordovician/BattleSim

max_forks_repo_name: ordovician/BattleSim


text:
using .Threads

abstract type Message end
struct Query <: Message end
struct Stop  <: Message end

function actor(ch::Channel)
    put!(ch, threadid())
    while true
        token = take!(ch)
        if token isa Query
            put!(ch, threadid())
        elseif token isa Stop
            break
        else
        end
    end
end

function startup()
    @threads for i in 1:nthreads()
        @async actor(CH[threadid()])
    end
end

const CH = Vector{Channel}(undef, nthreads())
const cnt = zeros(Int, nthreads())

function init_all()
    for i in eachindex(CH)
        CH[i] = Channel()
        cnt[i] = 0
    end
end

function results()
    for i in eachindex(CH)
        if isready(CH[i])
            cnt[i] += take!(CH[i]) ÷ i
        end
    end
end

function test_thrd(n::Int)
    init_all()
    for i in 1:n
        startup()
        results()
        foreach(c->put!(c,Stop()), CH)
    end
    cnt
end


text:
function test_NfOrdCls()
  Qx, x = PolynomialRing(FlintQQ, "x")

  K1, a1 = NumberField(x^3 - 2, "a")
  O1 = EquationOrder(K1)

  K2, a2 = NumberField(x - 2, "a")
  O2 = EquationOrder(K2)

  f3 = x^64 - 64*x^62 + 1952*x^60 - 37760*x^58 + 520144*x^56 - 5430656*x^54 + 44662464*x^52 - 296854272*x^50 + 1623421800*x^48 - 7398867840*x^46 + 28362326720*x^44 - 92043777280*x^42 + 254005423840*x^40 - 597659820800*x^38 + 1200442440064*x^36 - 2057901325824*x^34 + 3006465218196*x^32 - 3732682723968*x^30 + 3922021702720*x^28 - 3467892873984*x^26 + 2561511781920*x^24 - 1565841089280*x^22 + 782920544640*x^20 - 315492902400*x^18 + 100563362640*x^16 - 24754058496*x^14 + 4559958144*x^12 - 602516992*x^10 + 53796160*x^8 - 2968064*x^6 + 87296*x^4 - 1024*x^2 + 2

  K3, a3 = NumberField(f3, "a")
  O3 = EquationOrder(K3)

  @test_and_infer(signature, (O1, ), (1, 1))
  @test_and_infer(signature, O2, (1, 0))
  @test_and_infer(signature, (O3, ), (64, 0))

end


text:
"""

    random_graph(n,p=0.5, has_self_loops=false)

generate a random directed graph `g` with `n` nodes and the probability to include edge is `p`. If `has_self_loops`, `g` will include self loops.
"""
function random_graph(n::Integer, p::Real = 0.5; has_self_loops = false)
    g = DiGraph()
    for i = 1:n
        ind = 1
        for j = ind:n
            add_node!(g, j)
            if rand() <= p && (i != j || has_self_loops)
                add_edge!(g, i, j)
            end
        end
    end
    return g
end

"""
    random_evolving_graph(nv, nt, p=0.5, is_directed=true, has_self_loops=false)

Generate a random evolving graph `g` with `nv` nodes, `nt` timestamps. The probability to include each edge is equal to `p`. If `has_self_loops`,
we allow `g` to have self loops.
"""
function random_evolving_graph(g::IntAdjacencyList, nt::Int, p::Real = 0.5;
                               has_self_loops = false)
    nn = g.nnodes
    for t in 1:nt
        for i = 1:nn
            g.is_directed? ind =1 : ind = i
            for j = ind:nn
                    if (rand() <= p && (i != j || has_self_loops))
                        add_edge!(g, i, j, t)
                    end
            end
        end
    end
    g
end
function random_evolving_graph(nv::Int, nt::Int,p::Real = 0.5;
                               is_directed = true,has_self_loops  = false)
    g = IntAdjacencyList(nv, nt, is_directed = is_directed)
    random_evolving_graph(g, nt, p, has_self_loops = has_self_loops)
end


text:
@testset "Knapsack" begin
    mass    = [1, 5, 3, 7, 2, 10, 5]
    utility = [1, 3, 5, 2, 5,  8, 3]

    function fitness(n::AbstractVector)
        total_mass = sum(mass .* n)
        return (total_mass <= 20) ? sum(utility .* n) : 0
    end

    initpop = collect(rand(Bool,length(mass)))

    best, invbestfit, generations, tolerance, history = ga(
        x -> 1 / fitness(x),                    # Function to MINIMISE
        length(initpop),                        # Length of chromosome
        initPopulation = initpop,
        selection = roulette,                   # Options: sus
        mutation = inversion,                   # Options:
        crossover = singlepoint,                # Options:
        mutationRate = 0.2,
        crossoverRate = 0.5,
        ɛ = 0.1,                                # Elitism
        iterations = 20,
        populationSize = 50,
        interim = true);

    @test fitness(best) == 21.
    @test 1. /invbestfit == 21.
    @test sum(mass .* best) <= 20
end


text:
"""
    ot_plan(c, μ, ν; kwargs...)

Compute the optimal transport plan for the Monge-Kantorovich problem with source and target
marginals `μ` and `ν` and cost `c`.

The optimal transport plan solves
```math
\\inf_{\\gamma \\in \\Pi(\\mu, \\nu)} \\int c(x, y) \\, \\mathrm{d}\\gamma(x, y)
```
where ``\\Pi(\\mu, \\nu)`` denotes the couplings of ``\\mu`` and ``\\nu``.

See also: [`ot_cost`](@ref)
"""
function ot_plan end

"""
    ot_cost(c, μ, ν; kwargs...)

Compute the optimal transport cost for the Monge-Kantorovich problem with source and target
marginals `μ` and `ν` and cost `c`.

The optimal transport cost is the scalar value
```math
\\inf_{\\gamma \\in \\Pi(\\mu, \\nu)} \\int c(x, y) \\, \\mathrm{d}\\gamma(x, y)
```
where ``\\Pi(\\mu, \\nu)`` denotes the couplings of ``\\mu`` and ``\\nu``.

See also: [`ot_plan`](@ref)
"""
function ot_cost end

#############
# Discrete OT
#############

"""
    emd(μ, ν, C, optimizer)

Compute the optimal transport plan `γ` for the Monge-Kantorovich problem with source
histogram `μ`, target histogram `ν`, and cost matrix `C` of size `(length(μ), length(ν))`
which solves
```math
\\inf_{γ ∈ Π(μ, ν)} \\langle γ, C \\rangle.
```

The corresponding linear programming problem is solved with the user-provided `optimizer`.
Possible choices are `Tulip.Optimizer()` and `Clp.Optimizer()` in the `Tulip` and `Clp`
packages, respectively.
"""
function emd(μ, ν, C, model::MOI.ModelLike)
    # check size of cost matrix
    nμ = length(μ)
    nν = length(ν)
    size(C) == (nμ, nν) || error("cost matrix `C` must be of size `(length(μ), length(ν))`")
    nC = length(C)

    # define variables
    x = MOI.add_variables(model, nC)
    xmat = reshape(x, nμ, nν)

    # define objective function
    T = float(eltype(C))
    zero_T = zero(T)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(float.(vec(C)), x), zero_T),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    # add non-negativity constraints
    for xi in x
        MOI.add_constraint(model, MOI.SingleVariable(xi), MOI.GreaterThan(zero_T))
    end

    # add constraints for source
    for (i, μi) in zip(axes(xmat, 1), μ) # eachrow(xmat) is not available on Julia 1.0
        f = MOI.ScalarAffineFunction(
            [MOI.ScalarAffineTerm(one(μi), xi) for xi in view(xmat, i, :)], zero(μi)
        )
        MOI.add_constraint(model, f, MOI.EqualTo(μi))
    end

    # add constraints for target
    for (i, νi) in zip(axes(xmat, 2), ν) # eachcol(xmat) is not available on Julia 1.0
        f = MOI.ScalarAffineFunction(
            [MOI.ScalarAffineTerm(one(νi), xi) for xi in view(xmat, :, i)], zero(νi)
        )
        MOI.add_constraint(model, f, MOI.EqualTo(νi))
    end

    # compute optimal solution
    MOI.optimize!(model)
    status = MOI.get(model, MOI.TerminationStatus())
    status === MOI.OPTIMAL || error("failed to compute optimal transport plan: ", status)
    p = MOI.get(model, MOI.VariablePrimal(), x)
    γ = reshape(p, nμ, nν)

    return γ
end

"""
    emd2(μ, ν, C, optimizer; plan=nothing)

Compute the optimal transport cost (a scalar) for the Monge-Kantorovich problem with source
histogram `μ`, target histogram `ν`, and cost matrix `C` of size `(length(μ), length(ν))`
which is given by
```math
\\inf_{γ ∈ Π(μ, ν)} \\langle γ, C \\rangle.
```

The corresponding linear programming problem is solved with the user-provided `optimizer`.
Possible choices are `Tulip.Optimizer()` and `Clp.Optimizer()` in the `Tulip` and `Clp`
packages, respectively.

A pre-computed optimal transport `plan` may be provided.
"""
function emd2(μ, ν, C, optimizer; plan=nothing)
    γ = if plan === nothing
        # compute optimal transport plan
        emd(μ, ν, C, optimizer)
    else
        # check dimensions
        size(C) == (length(μ), length(ν)) ||
            error("cost matrix `C` must be of size `(length(μ), length(ν))`")
        size(plan) == size(C) || error(
            "optimal transport plan `plan` and cost matrix `C` must be of the same size",
        )
        plan
    end
    return dot(γ, C)
end

###################################
# Semidiscrete and continuous 1D OT
###################################

"""
    ot_plan(c, μ::ContinuousUnivariateDistribution, ν::UnivariateDistribution)

Compute the optimal transport plan for the Monge-Kantorovich problem with univariate
distributions `μ` and `ν` as source and target marginals and cost function `c` of
the form ``c(x, y) = h(|x - y|)`` where ``h`` is a convex function.

In this setting, the optimal transport plan is the Monge map
```math
T = F_\\nu^{-1} \\circ F_\\mu
```
where ``F_\\mu`` is the cumulative distribution function of `μ` and ``F_\\nu^{-1}`` is the
quantile function of `ν`.

See also: [`ot_cost`](@ref), [`emd`](@ref)
"""
function ot_plan(c, μ::ContinuousUnivariateDistribution, ν::UnivariateDistribution)
    # Use T instead of γ to indicate that this is a Monge map.
    T(x) = quantile(ν, cdf(μ, x))
    return T
end

"""
    ot_cost(
        c, μ::ContinuousUnivariateDistribution, ν::UnivariateDistribution; plan=nothing
    )

Compute the optimal transport cost for the Monge-Kantorovich problem with univariate
distributions `μ` and `ν` as source and target marginals and cost function `c` of
the form ``c(x, y) = h(|x - y|)`` where ``h`` is a convex function.

In this setting, the optimal transport cost can be computed as
```math
\\int_0^1 c(F_\\mu^{-1}(x), F_\\nu^{-1}(x)) \\mathrm{d}x
```
where ``F_\\mu^{-1}`` and ``F_\\nu^{-1}`` are the quantile functions of `μ` and `ν`,
respectively.

A pre-computed optimal transport `plan` may be provided.

See also: [`ot_plan`](@ref), [`emd2`](@ref)
"""
function ot_cost(
    c, μ::ContinuousUnivariateDistribution, ν::UnivariateDistribution; plan=nothing
)
    cost, _ = if plan === nothing
        quadgk(0, 1) do q
            return c(quantile(μ, q), quantile(ν, q))
        end
    else
        quadgk(0, 1) do q
            x = quantile(μ, q)
            return c(x, plan(x))
        end
    end
    return cost
end

################
# Discrete 1D OT
################

# internal iterator for discrete one-dimensional OT problems
# it returns tuples that consist of the indices of the source and target histograms
# and the optimal flow between the corresponding points
struct Discrete1DOTIterator{T,M,N}
    mu::M
    nu::N
end

# histograms `μ` and `ν` are expected to be iterators of the histograms where the
# corresponding support is sorted
function Discrete1DOTIterator(μ, ν)
    T = Base.promote_eltype(μ, ν)
    return Discrete1DOTIterator{T,typeof(μ),typeof(ν)}(μ, ν)
end

Base.IteratorEltype(::Type{<:Discrete1DOTIterator}) = Base.HasEltype()
Base.eltype(::Type{<:Discrete1DOTIterator{T}}) where {T} = Tuple{Int,Int,T}

Base.length(d::Discrete1DOTIterator) = length(d.mu) + length(d.nu) - 1

# we iterate through the source and target histograms
function Base.iterate(
    d::Discrete1DOTIterator{T}, (i, j, μnext, νnext)=(1, 1, iterate(d.mu), iterate(d.nu))
) where {T}
    # if we are done with iterating through the source and/or target histogram,
    # iteration is stopped
    if μnext === nothing || νnext === nothing
        return nothing
    end

    # unpack next values and states of the source and target histograms
    μiter, μstate = μnext
    νiter, νstate = νnext

    # compute next value of the iterator: indices of source and target histograms
    # and optimal flow between the corresponding points
    min_iter, max_iter = minmax(μiter, νiter)
    iter = (i, j, min_iter)

    # compute next state of the iterator
    diff = max_iter - min_iter
    state = if μiter < max_iter
        # move forward in the source histogram
        (i + 1, j, iterate(d.mu, μstate), (diff, νstate))
    else
        # move forward in the target histogram
        (i, j + 1, (diff, μstate), iterate(d.nu, νstate))
    end

    return iter, state
end

"""
    ot_plan(c, μ::DiscreteNonParametric, ν::DiscreteNonParametric)

Compute the optimal transport cost for the Monge-Kantorovich problem with univariate
discrete distributions `μ` and `ν` as source and target marginals and cost function `c`
of the form ``c(x, y) = h(|x - y|)`` where ``h`` is a convex function.

In this setting, the optimal transport plan can be computed analytically. It is returned as
a sparse matrix.

See also: [`ot_cost`](@ref), [`emd`](@ref)
"""
function ot_plan(_, μ::DiscreteNonParametric, ν::DiscreteNonParametric)
    # unpack the probabilities of the two distributions
    μprobs = probs(μ)
    νprobs = probs(ν)

    # create the iterator
    # note: support of `DiscreteNonParametric` is sorted
    iter = Discrete1DOTIterator(μprobs, νprobs)

    # create arrays for the indices of the two histograms and the optimal flow between the
    # corresponding points
    n = length(iter)
    I = Vector{Int}(undef, n)
    J = Vector{Int}(undef, n)
    W = Vector{Base.promote_eltype(μprobs, νprobs)}(undef, n)

    # compute the sparse optimal transport plan
    @inbounds for (idx, (i, j, w)) in enumerate(iter)
        I[idx] = i
        J[idx] = j
        W[idx] = w
    end
    γ = sparse(I, J, W, length(μprobs), length(νprobs))

    return γ
end

"""
    ot_cost(
        c, μ::DiscreteNonParametric, ν::DiscreteNonParametric; plan=nothing
    )

Compute the optimal transport cost for the Monge-Kantorovich problem with discrete
univariate distributions `μ` and `ν` as source and target marginals and cost function `c`
of the form ``c(x, y) = h(|x - y|)`` where ``h`` is a convex function.

In this setting, the optimal transport cost can be computed analytically.

A pre-computed optimal transport `plan` may be provided.

See also: [`ot_plan`](@ref), [`emd2`](@ref)
"""
function ot_cost(c, μ::DiscreteNonParametric, ν::DiscreteNonParametric; plan=nothing)
    return _ot_cost(c, μ, ν, plan)
end

# compute cost from scratch if no plan is provided
function _ot_cost(c, μ::DiscreteNonParametric, ν::DiscreteNonParametric, ::Nothing)
    # unpack the probabilities of the two distributions
    μprobs = probs(μ)
    νprobs = probs(ν)

    # create the iterator
    # note: support of `DiscreteNonParametric` is sorted
    iter = Discrete1DOTIterator(μprobs, νprobs)

    # compute the cost
    μsupport = support(μ)
    νsupport = support(ν)
    cost = sum(w * c(μsupport[i], νsupport[j]) for (i, j, w) in iter)

    return cost
end

# if a sparse plan is provided, we just iterate through the non-zero entries
function _ot_cost(
    c, μ::DiscreteNonParametric, ν::DiscreteNonParametric, plan::SparseMatrixCSC
)
    # extract non-zero flows
    I, J, W = findnz(plan)

    # compute the cost
    μsupport = support(μ)
    νsupport = support(ν)
    cost = sum(w * c(μsupport[i], νsupport[j]) for (i, j, w) in zip(I, J, W))

    return cost
end

# fallback: compute cost matrix (probably often faster to compute cost from scratch)
function _ot_cost(c, μ::DiscreteNonParametric, ν::DiscreteNonParametric, plan)
    return dot(plan, StatsBase.pairwise(c, support(μ), support(ν)))
end

################
# OT Gaussians
################

"""
    ot_cost(::SqEuclidean, μ::MvNormal, ν::MvNormal)

Compute the squared 2-Wasserstein distance between normal distributions `μ` and `ν` as
source and target marginals.

In this setting, the optimal transport cost can be computed as
```math
W_2^2(\\mu, \\nu) = \\|m_\\mu - m_\\nu \\|^2 + \\mathcal{B}(\\Sigma_\\mu, \\Sigma_\\nu)^2,
```
where ``\\mu = \\mathcal{N}(m_\\mu, \\Sigma_\\mu)``,
``\\nu = \\mathcal{N}(m_\\nu, \\Sigma_\\nu)``, and ``\\mathcal{B}`` is the Bures metric.

See also: [`ot_plan`](@ref), [`emd2`](@ref)
"""
function ot_cost(::SqEuclidean, μ::MvNormal, ν::MvNormal)
    return sqeuclidean(μ.μ, ν.μ) + sqbures(μ.Σ, ν.Σ)
end

"""
    ot_cost(::SqEuclidean, μ::Normal, ν::Normal)

Compute the squared 2-Wasserstein distance between univariate normal distributions `μ` and
`ν` as source and target marginals.

See also: [`ot_plan`](@ref), [`emd2`](@ref)
"""
function ot_cost(::SqEuclidean, μ::Normal, ν::Normal)
    return (μ.μ - ν.μ)^2 + (μ.σ - ν.σ)^2
end

"""
    ot_plan(::SqEuclidean, μ::MvNormal, ν::MvNormal)

Compute the optimal transport plan for the Monge-Kantorovich problem with multivariate
normal distributions `μ` and `ν` as source and target marginals and cost function
``c(x, y) = \\|x - y\\|_2^2``.

In this setting, for ``\\mu = \\mathcal{N}(m_\\mu, \\Sigma_\\mu)`` and
``\\nu = \\mathcal{N}(m_\\nu, \\Sigma_\\nu)``, the optimal transport plan is the Monge
map
```math
T \\colon x \\mapsto m_\\nu
+ \\Sigma_\\mu^{-1/2}
{\\big(\\Sigma_\\mu^{1/2} \\Sigma_\\nu \\Sigma_\\mu^{1/2}\\big)}^{1/2}\\Sigma_\\mu^{-1/2}
(x - m_\\mu).
```

See also: [`ot_cost`](@ref), [`emd`](@ref)
"""
function ot_plan(::SqEuclidean, μ::MvNormal, ν::MvNormal)
    Σμsqrt = μ.Σ^(-1 / 2)
    A = Σμsqrt * sqrt(_gaussian_ot_A(μ.Σ, ν.Σ)) * Σμsqrt
    mμ = μ.μ
    mν = ν.μ
    T(x) = mν + A * (x - mμ)
    return T
end

"""
    ot_plan(::SqEuclidean, μ::Normal, ν::Normal)

Compute the optimal transport plan for the Monge-Kantorovich problem with
normal distributions `μ` and `ν` as source and target marginals and cost function
``c(x, y) = \\|x - y\\|_2^2``.

See also: [`ot_cost`](@ref), [`emd`](@ref)
"""
function ot_plan(::SqEuclidean, μ::Normal, ν::Normal)
    mμ = μ.μ
    mν = ν.μ
    a = ν.σ / μ.σ
    T(x) = mν + a * (x - mμ)
    return T
end


text:
using BSON: @load
using Flux
using StatsPlots
# get the first layer parameters for influence analysis
@load "neuralmodel.bson" nnmodel
beta = nnmodel.layers[1].W # get first layer betas
z = maximum(abs.(beta),dims=1);
heatmap(z, xlabel="statistic", title="Importance of inputs, bright=high, dark=low")
savefig("ImportanceOfStatistics.png")

max_stars_repo_name: mcreel/SimulatedNeuralMoments

max_issues_repo_name: mcreel/SimulatedNeuralMoments


text:
# Before running this, please make sure to activate and instantiate the
# tutorial-specific package environment, using this
# [`Project.toml`](https://raw.githubusercontent.com/juliaai/DataScienceTutorials.jl/gh-pages/__generated/A-composing-models/Project.toml) and
# [this `Manifest.toml`](https://raw.githubusercontent.com/juliaai/DataScienceTutorials.jl/gh-pages/__generated/A-composing-models/Manifest.toml), or by following
# [these](https://juliaai.github.io/DataScienceTutorials.jl/#learning_by_doing) detailed instructions.

# ## Generating dummy data
# Let's start by generating some dummy data with both numerical values and categorical values:

using MLJ
using PrettyPrinting

KNNRegressor = @load KNNRegressor

# input

X = (age    = [23, 45, 34, 25, 67],
     gender = categorical(['m', 'm', 'f', 'm', 'f']))

# target

height = [178, 194, 165, 173, 168];

# Note that the scientific type of `age` is `Count` here:

scitype(X.age)

# We will want to coerce that to `Continuous` so that it can be given to a regressor that expects such values.

# ## Declaring a pipeline

# A typical workflow for such data is to one-hot-encode the categorical data and then apply some regression model on the data.
# Let's say that we want to apply the following steps:
# 1. One hot encode the categorical features in `X`
# 1. Standardize the target variable (`:height`)
# 1. Train a KNN regression model on the one hot encoded data and the Standardized target.

# The `Pipeline` constructor helps you define such a simple (non-branching) pipeline of steps to be applied in order:

pipe = Pipeline(
    coercer = X -> coerce(X, :age=>Continuous),
    one_hot_encoder = OneHotEncoder(),
    transformed_target_model = TransformedTargetModel(
        model = KNNRegressor(K=3);
        target=UnivariateStandardizer()
    )
)

# Note the coercion of the `:age` variable to Continuous since `KNNRegressor` expects `Continuous` input.
# Note also the `TransformedTargetModel` which allows one to learn a transformation (in this case Standardization) of the
# target variable to be passed to the `KNNRegressor`.

# Hyperparameters of this pipeline can be accessed (and set) using dot syntax:

pipe.transformed_target_model.model.K = 2
pipe.one_hot_encoder.drop_last = true;

# Evaluation for a pipe can be done with the `evaluate!` method; implicitly it will construct machines that will contain the fitted parameters etc:

evaluate(
    pipe,
    X,
    height,
    resampling=Holdout(),
    measure=rms
) |> pprint

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
##################
# `AbstractPass` #
##################

abstract type AbstractPass end

#########
# `Tag` #
#########

# this @pure annotation has official vtjnash approval :p
Base.@pure _pure_objectid(x) = objectid(x)

abstract type AbstractContextName end

struct Tag{N<:AbstractContextName,X,E#=<:Union{Nothing,Tag}=#} end

Tag(::Type{N}, ::Type{X}) where {N,X} = Tag(N, X, Nothing)

Tag(::Type{N}, ::Type{X}, ::Type{E}) where {N,X,E} = Tag{N,_pure_objectid(X),E}()

#################
# `BindingMeta` #
#################
# We define these here because we need them to define `Context`,
# but most code that works with these types is in src/tagged.jl

mutable struct BindingMeta
    data::Any
    BindingMeta() = new()
end

const BindingMetaDict = Dict{Symbol,BindingMeta}
const BindingMetaDictCache = IdDict{Module,BindingMetaDict}

#############
# `Context` #
#############

struct DisableHooks end

"""
```
Context{N<:Cassette.AbstractContextName,
        M<:Any,
        P<:Cassette.AbstractPass,
        T<:Union{Nothing,Cassette.Tag},
        B<:Union{Nothing,Cassette.BindingMetaDictCache},
        H<:Union{Nothing,Cassette.DisableHooks}}
```

A type representing a Cassette execution context. This type is normally interacted with
through type aliases constructed via `Cassette.@context`:

```
julia> Cassette.@context MyCtx
Cassette.Context{nametype(MyCtx),M,P,T,B,H} where H<:Union{Nothing,DisableHooks}
                                            where B<:Union{Nothing,IdDict{Module,Dict{Symbol,BindingMeta}}}
                                            where P<:AbstractPass
                                            where T<:Union{Nothing,Tag}
                                            where M
```

# Constructors

Given a context type alias named e.g. `MyCtx`, an instance of the type can be constructed via:

```
MyCtx(; metadata = nothing, pass = Cassette.NoPass())
```

To construct a new context instance using an existing context instance as a template, see
the `similarcontext` function.

To enable contextual tagging for a given context instance, see the [`enabletagging`](@ref) function.

# Fields

- `name::N<:Cassette.AbstractContextName`: a parameter used to disambiguate different
    contexts for overloading purposes (e.g. distinguishes `MyCtx` from other `Context` type
    aliases).

- `metadata::M<:Any`: trace-local metadata as provided to the context constructor

- `tag::T<:Union{Nothing,Tag}`: the tag object that is attached to values when they are
    tagged w.r.t. the context instance

- `pass::P<:Cassette.AbstractPass`: the Cassette pass that will be applied to all method
    bodies encountered during contextual execution (see the [`@pass`](@ref) macro for details).

- `bindingscache::B<:Union{Nothing,BindingMetaDictCache}}`: storage for metadata associated
    with tagged module bindings

- `hooktoggle::H<:Union{Nothing,DisableHooks}`: configuration toggle for disabling
    the `overdub` pass's `prehook`/`posthook` injection (see [`disablehooks`](@ref)
    for details)
"""
struct Context{N<:AbstractContextName,
               M<:Any,
               T<:Union{Nothing,Tag},
               P<:AbstractPass,
               B<:Union{Nothing,BindingMetaDictCache},
               H<:Union{Nothing,DisableHooks}}
    name::N
    metadata::M
    tag::T
    pass::P
    bindingscache::B
    hooktoggle::H
end

const ContextUntagged{N<:AbstractContextName} = Context{N,<:Any,Nothing}
const ContextTagged{T<:Tag,N<:AbstractContextName} = Context{N,<:Any,T}
const ContextWithPass{P<:AbstractPass,N<:AbstractContextName} = Context{N,<:Any,<:Union{Nothing,Tag},P}
const ContextWithHookToggle{H<:Union{Nothing,DisableHooks},N<:AbstractContextName} = Context{N,<:Any,<:Union{Nothing,Tag},<:AbstractPass,<:Union{Nothing,BindingMetaDictCache},H}

function Context(name::AbstractContextName; metadata = nothing, pass::AbstractPass = NO_PASS)
    return Context(name, metadata, nothing, pass, nothing, nothing)
end

"""
```
similarcontext(context::Context;
               metadata = context.metadata,
               pass = context.pass)
```

Return a copy of the given `context`, where the copy's `metadata` and/or `pass`
fields are replaced with those provided via the corresponding keyword arguments.
"""
function similarcontext(context::Context;
                        metadata = context.metadata,
                        pass = context.pass)
    return Context(context.name, metadata, context.tag, pass,
                   context.bindingscache, context.hooktoggle)
end

"""
```
disablehooks(context::Cassette.Context)
```

Return of copy of the given `context` with `prehook`/`posthook` injection
disabled for the context. Disabling hook injection can reduce IR bloat in
scenarios where these hooks are not being utilized.
"""
function disablehooks(context::Context)
    return Context(context.name, context.metadata, context.tag, context.pass,
                   context.bindingscache, DisableHooks())
end

"""
```
enabletagging(context::Cassette.Context, f)
```

Return a copy of the given `context` with the tagging system enabled for the contextual
execution of `f`.

Cassette uses the type of `f` to generate the `tag` field of the returned instance.

Note that it is generally unsafe to use the returned instance to contextually execute
functions other than `f`. Specifically, in cases of nested contextual execution where
both inner and outer contexts employ the tagging system, improper application of the
tagged system could cause (for example) separate contexts to erroneously interfere with
each other's metadata propagation.

See also: [`hastagging`](@ref)
"""
function enabletagging(context::Context, f)
    return Context(context.name, context.metadata,
                   Tag(typeof(context.name), typeof(f)),
                   context.pass, BindingMetaDictCache(), context.hooktoggle)
end

"""
```
hastagging(::Type{<:Cassette.Context})
```

Returns `true` if the given type indicates that the contextual tagging system is enabled
for context instances of the type, returns `false` otherwise.

# Example

```
julia> Cassette.@context MyCtx;

julia> ctx = MyCtx();

julia> Cassette.hastagging(typeof(ctx))
false

julia> ctx = Cassette.enabletagging(ctx, sum);

julia> Cassette.hastagging(typeof(ctx))
true
```

See also: [`enabletagging`](@ref)
"""
hastagging(::Type{<:ContextTagged}) = true
hastagging(::Type{<:ContextUntagged}) = false

hashooks(::Type{<:ContextWithHookToggle{Nothing}}) = true
hashooks(::Type{<:ContextWithHookToggle{DisableHooks}}) = false

tagtype(::C) where {C<:Context} = tagtype(C)
tagtype(::Type{<:ContextTagged{T}}) where {T} = T

nametype(::Type{<:Context{N}}) where {N} = N

passtype(::Type{<:ContextWithPass{P}}) where {P} = P

############
# @context #
############

"""
```
Cassette.@context Ctx
```

Define a new Cassette context type with the name `Ctx`. In reality, `Ctx` is simply a type
alias for `Cassette.Context{Cassette.nametype(Ctx)}`.

Note that `Cassette.overdub` is automatically overloaded w.r.t. `Ctx` to define several
primitives by default. A full list of these default primitives can be obtained by running:

```
methods(Cassette.overdub, (Ctx, Vararg{Any}))
```

Note also that many of the default primitives' signatures only match when contextual tagging
is enabled.

See also: [`Context`](@ref)
"""
macro context(_Ctx)
    @assert isa(_Ctx, Symbol) "context name must be a Symbol"
    # These are not necessarily hygienic, but allow for `@context C` to be
    # repeatedly declared (as a no-op for redundant declarations) without
    # error, similarly to repeated `struct` declarations.
    CtxName = esc(Symbol("##$(_Ctx)#Name"))
    CtxTagged = esc(Symbol("##$(_Ctx)#Tagged"))
    Ctx = esc(_Ctx)
    N, M, T, P, H = esc(:N), esc(:M), esc(:T), esc(:P), esc(:H)
    return quote
        struct $CtxName <: AbstractContextName end

        Base.show(io::IO, ::Type{$CtxName}) = print(io, "nametype(", $(string(_Ctx)), ")")

        const $Ctx{$M,$T<:Union{Nothing,Tag}} = Context{$CtxName,$M,$T}
        const $CtxTagged = ContextTagged{$T,$CtxName} where {$T<:Tag}

        $Ctx(; kwargs...) = Context($CtxName(); kwargs...)

        @inline Cassette.overdub(::C, ::Typeof(Tag), ::Type{N}, ::Type{X}) where {C<:$Ctx,N,X} = Tag(N, X, tagtype(C))

        @inline Cassette.overdub(ctx::$Ctx, ::typeof(Core._apply), f, args...) = Core._apply(overdub, (ctx, f), args...)

        # TODO: There are certain non-`Core.Builtin` functions which the compiler often
        # relies upon constant propagation/tfuncs to infer, instead of specializing on
        # them, such as `isdispatchtuple`. Such functions should generally be contextual
        # primitives by default for the sake of performance, and we should add more of
        # them here as we encounter them.
        @inline Cassette.overdub(ctx::$Ctx, f::Typeof(Base.isdispatchtuple), T::Type) = fallback(ctx, f, T)
        @inline Cassette.overdub(ctx::$Ctx, f::Typeof(Base.eltype), T::Type) = fallback(ctx, f, T)
        @inline Cassette.overdub(ctx::$Ctx, f::Typeof(Base.convert), T::Type, t::Tuple) = fallback(ctx, f, T, t)
        @inline Cassette.overdub(ctx::$Ctx{<:Any,Nothing}, f::Typeof(Core.kwfunc), x) = fallback(ctx, f, x)
        @inline Cassette.overdub(ctx::$Ctx{<:Any,Nothing}, f::Typeof(Base.getproperty), x::Any, s::Symbol) = fallback(ctx, f, x, s)

        # the below primitives are only active when the tagging system is enabled (`typeof(ctx) <: CtxTagged`)

        @inline Cassette.overdub(ctx::C, f::Typeof(tag), value, ::C, metadata) where {C<:$CtxTagged} = fallback(ctx, f, value, ctx, metadata)
        @inline Cassette.overdub(ctx::$CtxTagged, ::Typeof(Array{T,N}), undef::UndefInitializer, args...) where {T,N} = tagged_new_array(ctx, Array{T,N}, undef, args...)
        @inline Cassette.overdub(ctx::$CtxTagged, ::Typeof(Core.Module), args...) = tagged_new_module(ctx, args...)
        @inline Cassette.overdub(ctx::$CtxTagged, ::Typeof(Core.tuple), args...) = tagged_new_tuple(ctx, args...)
        @inline Cassette.overdub(ctx::$CtxTagged, ::Typeof(Base.nameof), args...) = tagged_nameof(ctx, m)
        @inline Cassette.overdub(ctx::$CtxTagged, ::Typeof(Core.getfield), args...) = tagged_getfield(ctx, args...)
        @inline Cassette.overdub(ctx::$CtxTagged, ::Typeof(Core.setfield!), args...) = tagged_setfield!(ctx, args...)
        @inline Cassette.overdub(ctx::$CtxTagged, ::Typeof(Core.arrayref), args...) = tagged_arrayref(ctx, args...)
        @inline Cassette.overdub(ctx::$CtxTagged, ::Typeof(Core.arrayset), args...) = tagged_arrayset(ctx, args...)
        @inline Cassette.overdub(ctx::$CtxTagged, ::Typeof(Base._growbeg!), args...) = tagged_growbeg!(ctx, args...)
        @inline Cassette.overdub(ctx::$CtxTagged, ::Typeof(Base._growend!), args...) = tagged_growend!(ctx, args...)
        @inline Cassette.overdub(ctx::$CtxTagged, ::Typeof(Base._growat!), args...) = tagged_growat!(ctx, args...)
        @inline Cassette.overdub(ctx::$CtxTagged, ::Typeof(Base._deletebeg!), args...) = tagged_deletebeg!(ctx, args...)
        @inline Cassette.overdub(ctx::$CtxTagged, ::Typeof(Base._deleteend!), args...) = tagged_deleteend!(ctx, args...)
        @inline Cassette.overdub(ctx::$CtxTagged, ::Typeof(Base._deleteat!), args...) = tagged_deleteat!(ctx, args...)
        @inline Cassette.overdub(ctx::$CtxTagged, ::Typeof(Core.typeassert), args...) = tagged_typeassert(ctx, args...)

        @inline function Cassette.overdub(ctx::$CtxTagged, f::Core.IntrinsicFunction, args...)
            if f === Base.sitofp
                return tagged_sitofp(ctx, args...)
            elseif f === Base.sle_int
                return tagged_sle_int(ctx, args...)
            else # TODO: add more cases
                return fallback(ctx, f, args...)
            end
        end

        $Ctx
    end
end

###############################
# contextual dispatch methods #
###############################

"""
```
prehook(context::Context, f, args...)
```

Overload this Cassette method w.r.t. a given context in order to define a new contextual
prehook for that context.

To understand when/how this method is called, see the documentation for [`overdub`](@ref).

Invoking `prehook` is a no-op by default (it immediately returns `nothing`).

See also: [`overdub`](@ref), [`posthook`](@ref), [`recurse`](@ref), [`fallback`](@ref)

# Examples

Simple trace logging:

```
julia> Cassette.@context PrintCtx;

julia> Cassette.prehook(::PrintCtx, f, args...) = println(f, args)

julia> Cassette.overdub(PrintCtx(), /, 1, 2)
float(1,)
AbstractFloat(1,)
Float64(1,)
sitofp(Float64, 1)
float(2,)
AbstractFloat(2,)
Float64(2,)
sitofp(Float64, 2)
/(1.0, 2.0)
div_float(1.0, 2.0)
0.5
```

Counting the number of method invocations with one or more arguments of a given type:

```
julia> mutable struct Count{T}
           count::Int
       end

julia> Cassette.@context CountCtx;

julia> Cassette.prehook(ctx::CountCtx{Count{T}}, f, arg::T, args::T...) where {T} = (ctx.metadata.count += 1)

# count the number of calls of the form `f(::Float64, ::Float64...)`
julia> ctx = CountCtx(metadata = Count{Float64}(0));

julia> Cassette.overdub(ctx, /, 1, 2)
0.5

julia> ctx.metadata.count
2
```
"""
@inline prehook(::Context, ::Vararg{Any}) = nothing

"""
```
posthook(context::Context, output, f, args...)
```

Overload this Cassette method w.r.t. a given context in order to define a new contextual
posthook for that context.

To understand when/how this method is called, see the documentation for [`overdub`](@ref).

Invoking `posthook` is a no-op by default (it immediately returns `nothing`).

See also: [`overdub`](@ref), [`prehook`](@ref), [`recurse`](@ref), [`fallback`](@ref)

# Examples

Simple trace logging:

```
julia> Cassette.@context PrintCtx;

julia> Cassette.posthook(::PrintCtx, output, f, args...) = println(output, " = ", f, args)

julia> Cassette.overdub(PrintCtx(), /, 1, 2)
1.0 = sitofp(Float64, 1)
1.0 = Float64(1,)
1.0 = AbstractFloat(1,)
1.0 = float(1,)
2.0 = sitofp(Float64, 2)
2.0 = Float64(2,)
2.0 = AbstractFloat(2,)
2.0 = float(2,)
0.5 = div_float(1.0, 2.0)
0.5 = /(1.0, 2.0)
0.5
```

Accumulate the sum of all numeric scalar outputs encountered in the trace:

```
julia> mutable struct Accum
           x::Number
       end

julia> Cassette.@context AccumCtx;

julia> Cassette.posthook(ctx::AccumCtx{Accum}, out::Number, f, args...) = (ctx.metadata.x += out)

julia> ctx = AccumCtx(metadata = Accum(0));

julia> Cassette.overdub(ctx, /, 1, 2)
0.5

julia> ctx.metadata.x
13.0
```
"""
@inline posthook(::Context, ::Vararg{Any}) = nothing

"""
```
fallback(context::Context, f, args...)
```

Overload this Cassette method w.r.t. a given context in order to define a new contextual
execution fallback for that context.

To understand when/how this method is called, see the documentation for [`overdub`](@ref) and
[`canrecurse`](@ref).

By default, invoking `fallback(context, f, args...)` will simply call `f(args...)` (with all
arguments automatically untagged, if `hastagging(typeof(context))`).

See also:  [`canrecurse`](@ref), [`overdub`](@ref), [`recurse`](@ref), [`prehook`](@ref), [`posthook`](@ref)
"""
@inline fallback(ctx::Context, args...) = call(ctx, args...)

@inline call(::ContextUntagged, f, args...) = f(args...)
@inline call(context::ContextTagged, f, args...) = untag(f, context)(ntuple(i -> untag(args[i], context), Val(nfields(args)))...)

# TODO: This is currently needed to force the compiler to specialize on the type arguments
# to `Core.apply_type`. In the future, it would be best for Julia's compiler to better handle
# varargs calls to such functions with type arguments, or at least provide a better way to
# force specialization on the type arguments.
@inline call(::ContextUntagged, f::typeof(Core.apply_type), ::Type{A}, ::Type{B}) where {A,B} = f(A, B)
@inline call(::ContextTagged, f::typeof(Core.apply_type), ::Type{A}, ::Type{B}) where {A,B} = f(A, B)

"""
```
canrecurse(context::Context, f, args...)
```

Return `true` if `f(args...)` has a lowered IR representation that Cassette can overdub,
return `false` otherwise.

Alternatively, but equivalently:

Return `false` if `recurse(context, f, args...)` directly translates to
`fallback(context, f, args...)`, return `true` otherwise.

Note that unlike `overdub`, `fallback`, etc., this function is not intended to be overloaded.

See also:  [`overdub`](@ref), [`fallback`](@ref), [`recurse`](@ref)
"""
@inline canrecurse(ctx::Context, f, @nospecialize(args...)) = !(isa(untag(f, ctx), Core.Builtin) || _iscompilerfunc(untag(f, ctx)))
@inline canrecurse(ctx::Context, ::typeof(Core._apply), f, @nospecialize(args...)) = Core._apply(canrecurse, (ctx, f), args...)
@inline canrecurse(ctx::Context, ::typeof(Core.invoke), f, @nospecialize(_), @nospecialize(args...)) = canrecurse(ctx, f, args...)

_iscompilerfunc(::F) where {F} = Core.Compiler.typename(F).module === Core.Compiler


text:
# ==============================================================================
# Max-plus system linear
# ==============================================================================

# ==============================================================================
# Creation of max-plus system linear

### First system
S1 = mpsyslin(mpsparse([1 2 3; 4 5 6; 7 8 9]), mpsparse([0;0;0]), mpsparse([0 0 0]), mpsparse(mpeye(Int64, 3,3)), mpzeros(Int64, 3,1))

@test S1.A == [1 2 3; 4 5 6; 7 8 9]
@test S1.B == reshape([0;0;0], length([0;0;0]), 1)
@test S1.C == [0 0 0]
@test S1.D == mpeye(Int64, 3,3)
@test S1.x0 == mpzeros(Int64, 3,1)

S2 = mpsyslin(mpsparse([1 2 3; 4 5 6; 7 8 9]), mpsparse([0;0;0]), mpsparse([0 0 0]), mpsparse(mpeye(Int64, 3,3)))
S3 = mpsyslin(mpsparse([1 2 3; 4 5 6; 7 8 9]), mpsparse([0;0;0]), mpsparse([0 0 0]))
S4 = mpsyslin(MP([1 2 3; 4 5 6; 7 8 9]), MP([0;0;0]), MP([0 0 0]), mpeye(Int64, 3,3), MP([0; 0; 0]))
S5 = mpsyslin(MP([1 2 3; 4 5 6; 7 8 9]), MP([0;0;0]), MP([0 0 0]), mpeye(Int64, 3,3))
S6 = mpsyslin(MP([1 2 3; 4 5 6; 7 8 9]), MP([0;0;0]), MP([0 0 0]))

@test S1.A == S2.A == S3.A == S4.A == S5.A == S6.A
@test S1.B == S2.B == S3.B == S4.B == S5.B == S6.B
@test S1.C == S2.C == S3.C == S4.C == S5.C == S6.C
@test S1.D == S2.D == S3.D == S4.D == S5.D == S6.D
@test S1.x0 == S2.x0 == S3.x0 == S5.x0 == S6.x0
@test S4.x0 == mpones(Int64, 3,1)

### Second system
S7 = mpsyslin(mpsparse([1 2; 3 4]), mpsparse([0; 0]), mpsparse([0 0]), mpsparse(mpeye(Int64, 2,2)), mpzeros(Int64, 2,1))

@test S7.A == [1 2; 3 4]
@test S7.B == reshape([0;0], length([0;0]), 1)
@test S7.C == [0 0]
@test S7.D == mpeye(Int64, 2,2)
@test S7.x0 == mpzeros(Int64, 2,1)

S8 = mpsyslin(mpsparse([1 2; 3 4]), mpsparse([0; 0]), mpsparse([0 0]), mpsparse(mpeye(Int64, 2,2)))
S9 = mpsyslin(mpsparse([1 2; 3 4]), mpsparse([0; 0]), mpsparse([0 0]))
S10 = mpsyslin(MP([1 2; 3 4]), MP([0; 0]), MP([0 0]), mpeye(Int64, 2,2), MP([0; 0]))
S11 = mpsyslin(MP([1 2; 3 4]), MP([0; 0]), MP([0 0]), mpeye(Int64, 2,2))
S12 = mpsyslin(MP([1 2; 3 4]), MP([0; 0]), MP([0 0]))

@test S7.A == S8.A == S9.A == S10.A == S11.A == S12.A
@test S7.B == S8.B == S9.B == S10.B == S11.B == S12.B
@test S7.C == S8.C == S9.C == S10.C == S11.C == S12.C
@test S7.D == S8.D == S9.D == S10.D == S10.D == S11.D == S12.D
@test S7.x0 == S8.x0 == S9.x0 == S11.x0 == S12.x0
@test S10.x0 == mpones(Int64, 2,1)

### Third system
A = MP([1 2; 3 4]); B = MP([5 6; 7 8]); C = MP([9 10; 11 12]); D = MP([13 14; 15 16])
x0 = MP([17; 18])

S13 = mpsyslin(A,B,C,D,x0)
S14 = mpsyslin(A,B,C,D)
S15 = mpsyslin(A,B,C)
@test S13.A == S14.A == S15.A == A
@test S13.B == S14.B == S15.B == B
@test S13.C == S14.C == S15.C == C
@test S13.D == S14.D == D
@test S13.x0 == reshape(x0, length(x0), 1)
@test S14.x0 == S15.x0 == mpzeros(Int64, 2, 1)

# ==============================================================================
# Simulation. Note here we should call the explicit systems but here
# they are the same than implicit.

res = mpsimul(S7, MP(1:10), false)
@test res == MP([37])

res = mpsimul(S7, MP(1:10), true)
@test res == MP([1 5 9 13 17 21 25 29 33 37])
@test res == mpsimul(S8, MP(1:10), true)
@test res == mpsimul(S9, MP(1:10), true)
@test res == mpsimul(S11, MP(1:10), true)
@test res == mpsimul(S12, MP(1:10), true)

@test mpsimul(S10, MP(1:10), true) == MP([4  8  12  16  20  24  28  32  36  40])

# ==============================================================================
# Composition of systems

S1 = mpsyslin(MP([1.0 2; 3 4]), MP([0.0; 0]), MP([0.0 0]), mpeye(Float64, 2,2))
S2 = mpsyslin(MP([1.0 2 3; 4 5 6; 7 8 9]), MP([0.0;0;0]), MP([0.0 0 0]), mpeye(Float64, 3,3))

# Diagonal composition
S = S1 | S2

@test S.A == [1.0 2 mp0 mp0 mp0; 3 4 mp0 mp0 mp0; mp0 mp0 1 2 3; mp0 mp0 4 5 6; mp0 mp0 7 8 9]
@test S.B == [mp1 mp0; mp1 mp0; mp0 mp1; mp0 mp1; mp0 mp1]
@test S.C == [mp1 mp1 mp0 mp0 mp0; mp0 mp0 mp1 mp1 mp1]
@test S.D == mpeye(Float64, 5,5)
@test S.x0 == mpzeros(Float64, 5,1)

# Parallel composition
S = S1 + S2

@test S.A == [1.0 2 mp0 mp0 mp0; 3 4 mp0 mp0 mp0; mp0 mp0 1 2 3; mp0 mp0 4 5 6; mp0 mp0 7 8 9]
@test S.B == mpones(Float64, 5,1)
@test S.C == mpones(Float64, 1,5)
@test S.D == mpeye(Float64, 5,5)
@test S.x0 == mpzeros(Float64, 5,1)

# Series composition
S = S1 * S2

@test S.A == [1.0 2 3 mp0 mp0; 4 5 6 mp0 mp0; 7 8 9 mp0 mp0; mp0 mp0 mp0 1 2; mp0 mp0 mp0 3 4]
@test S.B == [mpones(Float64, 3,1); dense(mpzeros(Float64, 2,1))]
@test S.C == [dense(mpzeros(Float64, 1,3)) mpones(Float64, 1,2)]
D = mpeye(Float64, 5,5); D[4:5, 1:3] .= mp1;
@test S.D == D
@test S.x0 == mpzeros(Float64, 5,1)

# Input in common
S = [S1; S2]

@test S.A == [1.0 2 mp0 mp0 mp0; 3 4 mp0 mp0 mp0; mp0 mp0 1 2 3; mp0 mp0 4 5 6; mp0 mp0 7 8 9]
@test S.B == mpones(Float64, 5,1)
@test S.C == [mp1 mp1 mp0 mp0 mp0; mp0 mp0 mp1 mp1 mp1]
@test S.D == mpeye(Float64, 5,5)
@test S.x0 == mpzeros(Float64, 5,1)

# Output addition FIXME syntax S = [S1, S2] is not working
S = [S1 S2]

@test S.A == [1.0 2 mp0 mp0 mp0; 3 4 mp0 mp0 mp0; mp0 mp0 1 2 3; mp0 mp0 4 5 6; mp0 mp0 7 8 9]
@test S.B == [mp1 mp0; mp1 mp0; mp0 mp1; mp0 mp1; mp0 mp1]
@test S.C == mpones(Float64, 1,5)
@test S.D == mpeye(Float64, 5,5)
@test S.x0 == mpzeros(Float64, 5,1)

# Feedback composition FIXME syntax /.
S = S1 / S2

@test S.A == [1.0 2 mp0 mp0 mp0; 3 4 mp0 mp0 mp0; mp0 mp0 1 2 3; mp0 mp0 4 5 6; mp0 mp0 7 8 9]
@test S.B == [mpones(Float64, 2,1); dense(mpzeros(Float64, 3,1))]
@test S.C == [mpones(Float64, 1,2) dense(mpzeros(Float64, 1,3))]
D = mpeye(Float64, 5,5); D[3:5, 1:2] .= mp1; D[1:2, 3:5] .= mp1;
@test S.D == D
@test S.x0 == mpzeros(Float64, 5,1)

# ==============================================================================
#

S1 = mpsyslin(MP([1.0 2; 3 4]), MP([0.0; 0]), MP([0.0 0]), mpeye(Float64, 2,2))

S = 5 * S1
@test S.A == [1.0 2; 3 4]
@test S.B == reshape([0.0; 0], length([0;0]), 1)
@test S.C == [5.0 5]
@test S.D == mpeye(Float64, 2,2)
@test S.x0 == mpzeros(Float64, 2,1)
@test S == MP(5.0) * S1

S = S1 * 5
@test S.A == [1.0 2; 3 4]
@test S.B == reshape([5.0; 5], length([0;0]), 1)
@test S.C == [0.0 0]
@test S.D == mpeye(Float64, 2,2)
@test S.x0 == mpzeros(Float64, 2,1)
@test S1 * 5 == S1 * MP(5.0)

# ==============================================================================
# Implicit to Explicit

S1 = mpsyslin(MP([1.0 2; 3 4]), MP([0.0; 0]), MP([0.0 0]), mpeye(Float64, 2,2))
S = mpexplicit(S1)

@test S.A == S1.A
@test S.B == S1.B
@test S.C == S1.C
@test S.D == mpeye(Float64, 2,2)
@test S.x0 == S1.x0

S2 = mpsyslin(MP([1.0 2 3; 4 5 6; 7 8 9]), MP([0.0;0;0]), MP([0.0 0 0]), mpeye(Float64, 3,3))
S = mpexplicit(S2)

@test S.A == S2.A
@test S.B == S2.B
@test S.C == S2.C
@test S.D == mpeye(Float64, 3,3)
@test S.x0 == S2.x0


text:
"Hello, world"
mystringexample1 = "Hello, world"
println(mystringexample1)

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
@testset "cg" begin

#test cg.jl

#positive semidefinite matrix
A=randn(200,100);  A=A'*A;
xt=randn(100)
b=A*xt

#inexact
(x,flag,relres,iter1) = cg(A,b,tol=1e-5,maxIter=1000,out=0);
@test norm(A*x-b)/norm(b) <= 1e-5

#test if very accurate solution can be achieved
(x,flag,relres,iter) = cg(A,b,tol=1e-14,maxIter=1000,out=0);
@test norm(A*x-b)/norm(b) <= 1.01e-14

#test with very good initial guess
x=deepcopy(xt) .+ eps();
(x,flag,relres,iter2) = cg(A,b,tol=1e-5,maxIter=1000,x=x,out=0);
@test norm(A*x-b)/norm(b) <= 1e-5
@test iter2<iter1 #test if the very good initial guess reduced the number of iterations

#test if the output is the input when the initial guess is the true solution
x=deepcopy(xt);
(x,flag,relres,iter2) = cg(A,b,tol=1e-14,maxIter=1000,x=x,out=0);
@test norm(A*x-b)/norm(b) <= 1e-14
@test iter2==1
@test x==xt

#test if simple preconditioner reduces the number of iterations
A=diagm(0 => (1:1:100))
(x1,flag,relres1,iter1) = cg(A,b,tol=1e-12,maxIter=1000,out=0); #unpreconditioned cg
Prec(input) = input./diag(A); #Jacobi preconditioner
(x2,flag,relres2,iter2) = cg(A,b,tol=1e-12,maxIter=1000,M=Prec,out=0);
@test iter2<iter1

end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
function hello()
    # https://docs.julialang.org/en/v1/manual/functions/#The-return-Keyword
    return "Hello, World!"
end

max_stars_repo_name: joaopalmeiro/learning-julia

max_issues_repo_name: joaopalmeiro/learning-julia

max_forks_repo_name: joaopalmeiro/learning-julia


text:
using AnalogDiscovery

function go(fs::Float64)
  cDev = enumDevices()
  @show ("Available",cDev)
  myDev = Int32(0)
  if cDev > 0
    if enumDeviceIsOpened(myDev)
      error("Error opening device, already in use.")
    end
    hdwf = deviceOpen(myDev)
  else
    error("No devices available.")
  end

  analogInFrequencySet(hdwf,fs)
  bmin,bmax = analogInBufferSizeInfo(hdwf)
  bsize = min(Int32(round(fs)),bmax)
  analogInBufferSizeSet(hdwf,bsize)
  analogInAcquisitionModeSet(hdwf,ACQMODERECORD)

  chVoltMin,chVoltMax,chSteps = analogInChannelRangeInfo(hdwf)
  chCount = analogInChannelCount(hdwf)
  for ch in Int32(0):Int32(chCount-1)
    analogInChannelEnableSet(hdwf,ch,true)
    analogInChannelFilterSet(hdwf,ch,FILTERAVERAGE)
    analogInChannelRangeSet(hdwf,ch,chVoltMax)
  end

  # sleep(2.0)
  # Collect 20-sec worth of data from all channels
  blks = Int(ceil(20*fs/bsize))
  chData = [Vector{Float64}(blks*bsize) for i in 1:chCount]
  analogInRecordLengthSet(hdwf,0.0)
  analogInConfigure(hdwf,false,true) # start acquistion

  # Similar to AnalogIn_Record.py example in the Waveforms SDK
  cSamples = 0
  nSamples = blks*bsize
  chDataBuf = Vector{Float64}(bsize)
  while cSamples < nSamples
    status = analogInStatus(hdwf,true)
    if (cSamples == 0) && ((status == DWFSTATECONFIG) || (status == DWFSTATEPREFILL) || (status == DWFSTATEARMED))
      # Acquisition not yet started
      continue
    end
    cAvailable,cLost,cCorrupted = analogInStatusRecord(hdwf)
    # @show (cSamples,cAvailable,cLost,cCorrupted)
    cSamples += cLost
    if cAvailable==0
      continue
    end
    if cSamples+cAvailable > nSamples
      cAvailable = nSamples-cSamples
    end
    # get samples
    for ch in Int32(0):Int32(chCount-1)
      analogInStatusData!(chDataBuf,hdwf,ch,Int32(cAvailable))
      chData[ch+1][cSamples+(1:cAvailable)] = chDataBuf[1:cAvailable]
    end
    cSamples += cAvailable
  end

  deviceClose(hdwf)
  return chData
end

const fs = 512.0
chData = go(fs)


text:
# generate guassian data from a precision matrix
function generateData(Ω, n)
  p = size(Ω, 1)
  d, U = eig(Ω)
  for i=1:p
    d[i] = 1. / sqrt(d[i])
  end
  X = randn(n, p) * (U * diagm(d) * U')
end

############################

function generateSparseSPD(p, α=0.95, smallest_coef=.1, largest_coef=.9)

  aux = rand(p, p)
  d = Uniform(smallest_coef, largest_coef)
  for c=1:p, r=c+1:p
    aux[r,c] = aux[r,c] < α ? 0. : rand(d)
  end
  tril!(aux, -1)

  # permutation = randperm(p)
  # aux = aux[permutation, :]
  # aux = aux[:, permutation]
  # @show aux
  chol = aux - eye(p)
  prec = chol * chol'

  return prec
end

# takes as an input a precision matrix,
# returns a new precision matrix
# off diagonal elements have some probability of having flipped sign
function generateRandomDifference(Ω::StridedMatrix, probability_flip)

  p = size(Ω, 1)
  Ωy = copy(Ω)

  for c=1:p, r=c+1:p
    if rand() < probability_flip
      v = -Ω[r,c]
      Ωy[r,c] = Ωy[c,r] = v
    end
  end

  return Ωy
end

# ρ -- strength of an edge
# generate
function genLatentModel(Ω::StridedMatrix, r, relStrength)
  p = size(Ω, 1)
  evmax = eigmax(Ω)

  U = randn(p, r)
  L = Symmetric( U*U' )
  L = L / norm(L, 2.) * relStrength * evmax

  (Ω, L)
end


text:
using Sigma
using AbstractDomains

a = flip(1)
omegatron = Sigma.abstract_omega(1)
omegatron = HyperBox([Interval(0.501,.51)])
b = convert(Sigma.DRealRandVar{Bool},a)
@show call(b,omegatron)


text:
using Oceananigans
using Oceananigans.Units
using Oceananigans.TurbulenceClosures: Horizontal
using GLMakie

grid = RectilinearGrid(size = (60, 60, 1),
                       x = (0, 1200kilometers),
                       y = (0, 1200kilometers),
                       z = (-4kilometers, 0),
                       halo = (3, 3, 3),
                       topology = (Bounded, Bounded, Bounded))

# Specify coriolis, viscosity / turbulence closure, and boundary condition
coriolis = BetaPlane(f₀=1e-4, β=1e-11)
closure = ScalarDiffusivity(ν=100, isotropy=Horizontal())

no_slip = ValueBoundaryCondition(0)
wind_stress(x, y, t) = - 1e-4 * cos(2π * x / 1200kilometers)
boundary_conditions = (;
    u = FieldBoundaryConditions(top=FluxBoundaryCondition(wind_stress), south=no_slip, north=no_slip),
    v = FieldBoundaryConditions(east=no_slip, west=no_slip)
)

model = HydrostaticFreeSurfaceModel(; grid, coriolis, closure, boundary_conditions,
                                      momentum_advection = UpwindBiasedFifthOrder(),
                                      tracers = (), buoyancy = nothing)

# Time-stepping
simulation = Simulation(model, Δt=20seconds, stop_time=5days)

progress(sim) = @info "Iter: $(iteration(sim)), time: $(prettytime(sim))"
simulation.callbacks[:progress] = Callback(progress, IterationInterval(100))

run!(simulation)

# Visualize the x-velocity field
fig, ax, pl = heatmap(interior(model.velocities.u)[:, :, 1])
display(fig)


text:
module CUDAapi

using Libdl
using Logging

# FIXME: replace with an additional log level when we depend on 0.7+
macro trace(ex...)
    esc(:(@debug $(ex...)))
end

include("util.jl")

include("discovery.jl")
include("availability.jl")

include("complex.jl")
include("library_types.jl")

include("call.jl")

end


text:
module Curls

export Curl

using ...Common
using SparseArrays

struct Curl{N,T}
	curlB::SparseMatrixCSC{T,Int}
	curlE::SparseMatrixCSC{T,Int}
end

function Curl(sim::Simulation{1})

	∂ₓE = E_derivative(length(sim),sim.dx)
	posE = map(p->p.x,sim.x)
	PE = invperm(sortperm(posE))

	∂ₓB = B_derivative(length(sim),sim.dx)
	posB = expand_B_pos(posE,sim.dx)
	PB = invperm(sortperm(posB))

	permute!(∂ₓE,PB,PE) # not sure why this is the correct permutation
	permute!(∂ₓB,PE,PB) # but it works

	I3 = sparse([3,2],[2,3],[1,-1],3,3)
	curlB = kron(I3,∂ₓB)
	curlE = kron(I3,∂ₓE)

	return Curl{1,Float64}(curlB,curlE)
end

function Curl(sim::Simulation{2}) end

function Curl(sim::Simulation{3}) end

function B_derivative(n::Int,dx::Float64)
	rowsm = 1:n
	rowsp = 1:n
	colsm = 1:n
	colsp = 2:n+1
	vals = ones(n)
	return sparse(vcat(rowsm,rowsp),vcat(colsm,colsp),vcat(-vals,vals)/dx,n,n+1)
end

function E_derivative(n::Int,dx::Float64)
	rowsm = 2:n
	rowsp = 2:n
	colsm = 1:n-1
	colsp = 2:n
	vals = ones(n-1)
	return sparse(vcat([1,n+1],rowsm,rowsp),vcat([1,n],colsm,colsp),vcat([2,-2],-vals,vals)/dx,n+1,n)
end

expand_B_pos(sim::Simulation{1}) = expand_B_pos(map(p->p.x,sim.x),sim.dx)
expand_B_pos(posE::Vector,dx::Float64) = vcat(posE.-dx/2,maximum(posE)+dx/2)

end # module


text:
include("solvers/azsolver.jl")
include("solvers/generic_azsolver.jl")
include("solvers/lowranksolver.jl")
include("solvers/randomized.jl")

# solver needs also a problemstyle and a solverstyle
solver(ap::ApproximationProblem;
            problemstyle = ProblemStyle(ap),
            samplingstyle = SamplingStyle(ap),
            solverstyle = SolverStyle(samplingstyle, ap), options...) =
        solver(problemstyle, solverstyle, ap; samplingstyle=samplingstyle, options...)

# solver needs a solverstyle ap and operator
function solver(pstyle::DictionaryOperatorStyle, solverstyle::SolverStyle, ap::ApproximationProblem;
            samplingstyle = SamplingStyle(ap), options...)
    S = samplingoperator(samplingstyle, ap; options...)
    A = discretization(samplingstyle, ap, S; options...)
    solver(solverstyle, ap, A; S = S, options...)
end

# samplingoperator is no TensorProductOperator
function solver(pstyle::DictionaryOperatorStyle, solverstyle::ProductSolverStyle, ap::ApproximationProblem;
            samplingstyle = SamplingStyle(ap), options...)
    S = samplingoperator(samplingstyle, ap; options...)
    A = discretization(samplingstyle, ap; options...)
    solver(solverstyle, ap, A; S = S, options...)
end

function solver(::InverseStyle, ap::ApproximationProblem, A::AbstractOperator; weightedAZ=false, options...)
    if weightedAZ
        AZ_Cweight = haskey(options,:AZ_Cweight) ? options[:AZ_Cweight] : error("No options `AZ_Cweight`")
        AZ_Cweight*inv(A*AZ_Cweight)
    else
        inv(A)
    end
end
function solver(::DirectStyle, ap::ApproximationProblem, A::AbstractOperator; weightedAZ=false, options...)
    if weightedAZ
        AZ_Cweight = haskey(options,:AZ_Cweight) ? options[:AZ_Cweight] : error("No options `AZ_Cweight`")
        AZ_Cweight*directsolver(A*AZ_Cweight; options...)
    else
        directsolver(A; options...)
    end
end
function solver(::IterativeStyle, ap::ApproximationProblem, A::AbstractOperator; weightedAZ=false, options...)
    if weightedAZ
        AZ_Cweight = haskey(options,:AZ_Cweight) ? options[:AZ_Cweight] : error("No options `AZ_Cweight`")
        AZ_Cweight*iterativesolver(A*AZ_Cweight; options...)
    else
        iterativesolver(A; options...)
    end
end

solver(style::AZStyle, ap::ApproximationProblem, A::AbstractOperator; problemstyle=ProblemStyle(ap), options...) =
    solver(problemstyle, style, ap, A; options...)
solver(pstyle::GenericOperatorStyle, solverstyle::AZStyle, ap::ApproximationProblem, A; options...) =
    GenericAZSolver(ap, A; options...)
solver(pstyle::DictionaryOperatorStyle, style::AZStyle, ap::ApproximationProblem, A::AbstractOperator; options...) =
    solver(style, ap, A, AZ_Zt(pstyle, ap; (options)...); options...)

function solver(::AZStyle, ap::ApproximationProblem, A::AbstractOperator, Zt::AbstractOperator;
            B=nothing, smallcoefficients=false, smallcoefficients_atol=NaN, smallcoefficients_rtol=NaN, verbose=false, options...)
    if smallcoefficients
        w = BasisFunctions.quadweights(sampling_grid(ap; options...), measure(ap; options...))
        normF = abs(sqrt(sum(w .* B.^2)))
        if !isnan(smallcoefficients_rtol)
            verbose && println("Change smallcoefficients relative tolerance to absolute tolerance rtol*||f||")
            smallcoefficients_atol = smallcoefficients_rtol*normF
            smallcoefficients_rtol = NaN
        end

        AZSolver(A, Zt; smallcoefficients=smallcoefficients, smallcoefficients_rtol=smallcoefficients_rtol,
                smallcoefficients_atol=smallcoefficients_atol, verbose=verbose, options...)
    else
        AZSolver(A, Zt; B=B, verbose=verbose, options...)
    end
end

solver(::DualStyle, ap::ApproximationProblem, A::AbstractOperator; options...) = dualdiscretization(ap; options...)'

solver(solverstyle::ProductSolverStyle, ap::ApproximationProblem, A::AbstractOperator; samplingstyle=SamplingStyle(ap), options...) =
    solver(solverstyle, samplingstyle, ap, A; options...)
function solver(solverstyle::ProductSolverStyle, samplingstyle::ProductSamplingStyle, ap::ApproximationProblem, A::AbstractOperator; S, options...)
    solvere, sse, ape, Ae, Se = elements(solverstyle), elements(samplingstyle), elements(ap), elements(A), BasisFunctions.productelements(S)
    @assert length(solvere) == length(sse) == length(ape) == length(Ae) == length(Se)
    TensorProductOperator(
		map( (solversi, ssi, api, Ai, Si)->solver(solversi, api, Ai; S=Si, samplingstyle=ssi, options...), solvere, sse, ape, Ae, Se)...
	)
end

solver(style::AZSmoothStyle, ap::ApproximationProblem, A::AbstractOperator; options...) = solver(style, ap, A, AZ_Zt(ap; options...); options...)
include("solvers/smoothsolver.jl")
solver(::AZSmoothStyle, ap::ApproximationProblem, A::AbstractOperator, Zt; options...) =
    solver(AZStyle(), ap, A, Zt; options..., weightedAZ=true, AZ_Cweight=sobolevAZweight(A; options...))


text:
using Test
using ForwardDiff;
using FinancialToolbox
Dual_=ForwardDiff.Dual
print_colored("Starting Forward Diff Dual Numbers Test\n",:green)

#Test Parameters
testToll=1e-14;
spot=10;K=10;r=0.02;T=2.0;sigma=0.2;d=0.01;
assert_(value,toll)=@test abs(value)<toll
#EuropeanCall Option
PriceCall=blsprice(spot,K,r,T,sigma,d);
DeltaCall=blsdelta(spot,K,r,T,sigma,d);
ThetaCall=blstheta(spot,K,r,T,sigma,d);
RhoCall=blsrho(spot,K,r,T,sigma,d);
LambdaCall=blslambda(spot,K,r,T,sigma,d);
VannaCall=blsvanna(spot,K,r,T,sigma,d);

#EuropeanPut Option
PricePut=blsprice(spot,K,r,T,sigma,d,false);
DeltaPut=blsdelta(spot,K,r,T,sigma,d,false);
ThetaPut=blstheta(spot,K,r,T,sigma,d,false);
RhoPut=blsrho(spot,K,r,T,sigma,d,false);
LambdaPut=blslambda(spot,K,r,T,sigma,d,false);
VannaPut=blsvanna(spot,K,r,T,sigma,d,false);

#Equals for both Options
Gamma=blsgamma(spot,K,r,T,sigma,d);
Vega=blsvega(spot,K,r,T,sigma,d);

########DUAL NUMBERS
DerToll=1e-13;
#Function definition
#Call
FcallDual(spot)=blsprice(spot,K,r,T,sigma,d);
GcallDual(r)=blsprice(spot,K,r,T,sigma,d);
HcallDual(T)=blsprice(spot,K,r,T,sigma,d);
LcallDual(sigma)=blsprice(spot,K,r,T,sigma,d);
PcallDual(spot)=blsprice(spot,K,r,T,sigma,d)*spot.value/blsprice(spot.value,K,r,T,sigma,d);
VcallDual(sigma)=blsdelta(spot,K,r,T,sigma,d);
#Put
FputDual(spot)=blsprice(spot,K,r,T,sigma,d,false);
GputDual(r)=blsprice(spot,K,r,T,sigma,d,false);
HputDual(T)=blsprice(spot,K,r,T,sigma,d,false);
PputDual(spot)=blsprice(spot,K,r,T,sigma,d,false)*spot.value/blsprice(spot.value,K,r,T,sigma,d,false);
VPutDual(sigma)=blsdelta(spot,K,r,T,sigma,d,false);
#Input
SpotDual=Dual_(spot,1.0);
rDual=Dual_(r,1.0);
TDual=Dual_(T,1.0);
SigmaDual=Dual_(sigma,1.0);

#Automatic Differentiation Test
#TEST
print_colored("--- European Call Sensitivities: DualNumbers\n",:yellow)
print_colored("-----Testing Delta\n",:blue);
assert_(FcallDual(SpotDual).partials[1]-DeltaCall,DerToll)
print_colored("-----Testing Rho\n",:blue);
assert_(GcallDual(rDual).partials[1]-RhoCall,DerToll)
print_colored("-----Testing Theta\n",:blue);
assert_(-HcallDual(TDual).partials[1]-ThetaCall,DerToll)
print_colored("-----Testing Lambda\n",:blue);
assert_(PcallDual(SpotDual).partials[1]-LambdaCall,DerToll)
print_colored("-----Testing Vanna\n",:blue);
assert_(VcallDual(SigmaDual).partials[1]-VannaCall,DerToll)
#TEST
print_colored("--- European Put Sensitivities: DualNumbers\n",:yellow)
print_colored("-----Testing Delta\n",:blue);
assert_(FputDual(SpotDual).partials[1]-DeltaPut,DerToll)
print_colored("-----Testing Rho\n",:blue);
assert_(GputDual(rDual).partials[1]-RhoPut,DerToll)
print_colored("-----Testing Theta\n",:blue);
assert_(-HputDual(TDual).partials[1]-ThetaPut,DerToll)
print_colored("-----Testing Lambda\n",:blue);
assert_(PputDual(SpotDual).partials[1]-LambdaPut,DerToll)
print_colored("-----Testing Vanna\n",:blue);
assert_(VPutDual(SigmaDual).partials[1]-VannaPut,DerToll)

print_colored("-----Testing Vega\n",:blue);
assert_(LcallDual(SigmaDual).partials[1]-Vega,DerToll)
print_colored("Dual Numbers Test Passed\n",:green)
println("")

#TEST OF INPUT VALIDATION
print_colored("Starting Input Validation Test Dual\n",:magenta)
print_colored("----Testing Negative  Spot Price \n",:cyan)
@test_throws(ErrorException, blsprice(-SpotDual,K,r,T,sigma,d))
@test_throws(ErrorException, blkprice(-SpotDual,K,r,T,sigma))
@test_throws(ErrorException, blsdelta(-SpotDual,K,r,T,sigma,d))
@test_throws(ErrorException, blsgamma(-SpotDual,K,r,T,sigma,d))
@test_throws(ErrorException, blstheta(-SpotDual,K,r,T,sigma,d))
@test_throws(ErrorException, blsrho(-SpotDual,K,r,T,sigma,d))
@test_throws(ErrorException, blspsi(-SpotDual,K,r,T,sigma,d))
@test_throws(ErrorException, blslambda(-SpotDual,K,r,T,sigma,d))
@test_throws(ErrorException, blsvanna(-SpotDual,K,r,T,sigma,d))
@test_throws(ErrorException, blsvega(-SpotDual,K,r,T,sigma,d))

print_colored("----Testing Negative  Strike Price \n",:cyan)
KK=Dual_(K,0);
@test_throws(ErrorException, blsprice(spot,-KK,r,T,sigma,d))
@test_throws(ErrorException, blkprice(spot,-KK,r,T,sigma))
@test_throws(ErrorException, blsdelta(spot,-KK,r,T,sigma,d))
@test_throws(ErrorException, blsgamma(spot,-KK,r,T,sigma,d))
@test_throws(ErrorException, blstheta(spot,-KK,r,T,sigma,d))
@test_throws(ErrorException, blsrho(spot,-KK,r,T,sigma,d))
@test_throws(ErrorException, blspsi(spot,-KK,r,T,sigma,d))
@test_throws(ErrorException, blslambda(spot,-KK,r,T,sigma,d))
@test_throws(ErrorException, blsvanna(spot,-KK,r,T,sigma,d))
@test_throws(ErrorException, blsvega(spot,-KK,r,T,sigma,d))

print_colored("----Testing Negative  Time to Maturity \n",:cyan)
@test_throws(ErrorException, blsprice(spot,K,r,-TDual,sigma,d))
@test_throws(ErrorException, blkprice(spot,K,r,-TDual,sigma))
@test_throws(ErrorException, blsdelta(spot,K,r,-TDual,sigma,d))
@test_throws(ErrorException, blsgamma(spot,K,r,-TDual,sigma,d))
@test_throws(ErrorException, blstheta(spot,K,r,-TDual,sigma,d))
@test_throws(ErrorException, blsrho(spot,K,r,-TDual,sigma,d))
@test_throws(ErrorException, blspsi(spot,K,r,-TDual,sigma,d))
@test_throws(ErrorException, blslambda(spot,K,r,-TDual,sigma,d))
@test_throws(ErrorException, blsvanna(spot,K,r,-TDual,sigma,d))
@test_throws(ErrorException, blsvega(spot,K,r,-TDual,sigma,d))

print_colored("----Testing Negative  Volatility \n",:cyan)
@test_throws(ErrorException, blsprice(spot,K,r,T,-SigmaDual,d))
@test_throws(ErrorException, blkprice(spot,K,r,T,-SigmaDual))
@test_throws(ErrorException, blsdelta(spot,K,r,T,-SigmaDual,d))
@test_throws(ErrorException, blsgamma(spot,K,r,T,-SigmaDual,d))
@test_throws(ErrorException, blstheta(spot,K,r,T,-SigmaDual,d))
@test_throws(ErrorException, blsrho(spot,K,r,T,-SigmaDual,d))
@test_throws(ErrorException, blspsi(spot,K,r,T,-SigmaDual,d))
@test_throws(ErrorException, blslambda(spot,K,r,T,-SigmaDual,d))
@test_throws(ErrorException, blsvanna(spot,K,r,T,-SigmaDual,d))
@test_throws(ErrorException, blsvega(spot,K,r,T,-SigmaDual,d))

print_colored("Dual Input Validation Test Passed\n",:magenta)


text:
using DistributionVectors
using Test, Distributions, LinearAlgebra, Missings, StatsBase, MissingStrategies

function boot_dvsum(dv, nboot = 100_000)
    nboot = 100_000
    x = rand(dv, nboot);
    sums = vec(sum(x, dims = 2));
    params(fit(LogNormal,sums))
end

function boot_dvsums_acf(dv, acf, nboot = 10_000)
    # test sum formula by bootstrap sample
    μ, σ = params(dv)
    Sigma = Diagonal(σ) * cormatrix_for_acf(length(dv), acf) * Diagonal(σ);
    dn = MvNormal(disallowmissing(μ), Symmetric(Sigma));
    x = rand(dn, nboot) .|> exp
    sums = vec(sum(x, dims = 1))
    #density(sums)
    drsum = fit(LogNormal, sums)
end

@testset "sumlognormals" begin

@testset "two vars uncorrelated" begin
    # generate nSample values of two lognormal random variables
    d1 = LogNormal(log(110), 0.25)
    d2 = LogNormal(log(100), 0.15)
    dv = SimpleDistributionVector(d1, d2);
    dsum = @inferred sum(dv)
    @testset "no missings" begin
        @test params(dsum)[1] ≈ log(210) rtol = 0.02
        # checked with random numbers
        #boot_dvsum(dv)
        @test params(dsum)[2] ≈ 0.15 rtol = 0.02
    end;
    @testset "with missings" begin
        dv = SimpleDistributionVector(d1, missing);
        @test isequal(collect(dv),[d1, missing])
        @test_throws Exception dsum2 = sum(dv) # without skipmissing
        dsum2 = @inferred sum(skipmissing(dv))
        @test dsum2 == d1
        #a,b = sum(dv,SkipMissing())
        dsum3 = @inferred sum(dv,SkipMissing())
        @test dsum3 == d1
        #@btime sum(skipmissing($dv)) # does not allocate
    end;
    @testset "with gapfilling flag" begin
      dv = SimpleDistributionVector(d1, d2, d1);
      isgapfilled = [true, false, false]
      dsum4 = @inferred sum(dv)
      dsum = @inferred sum(dv, isgapfilled=isgapfilled)
      #@code_warntype sum(dv, isgapfilled=isgapfilled)
      @test mean(dsum) == mean(dsum4)
      @test std(dsum) > std(dsum4)
      # mean function
      dmean = @inferred mean(dv; isgapfilled=isgapfilled)
      @test dmean.μ ≈ dsum.μ / 3
      # relative error like non-gapfilled
      dmeann = mean(dv[.!isgapfilled])
      @test std(dmean)/mean(dmean) ≈ std(dmeann)/mean(dmeann)
    end;
    @testset "with missings and gapfilling flag" begin
      dv = SimpleDistributionVector(d1, d2, d1, missing);
      isgapfilled = [true, true, false, false]
      dsum5 = @inferred sum(dv[1:3], isgapfilled=isgapfilled[1:3])
      dsum = @inferred sum(dv, isgapfilled=isgapfilled, SkipMissing())
      #@code_warntype sum(dv, isgapfilled=isgapfilled, SkipMissing())
      @test dsum == dsum5
    end;
end;

@testset "few correlated vars" begin
  mu = log.([110,100,80,120,160.0])
  sigma = log.([1.2,1.5,1.1,1.3,1.1])
  acf1 = @inferred AutoCorrelationFunction([1, 0.4,0.1])
  n = length(mu)
  corrM = @inferred cormatrix_for_acf(n, acf1)
  dv = SimpleDistributionVector(LogNormal{eltype(mu)}, mu, sigma);
  mum = allowmissing(mu); mum[1] = missing
  dvm = SimpleDistributionVector(LogNormal{eltype(mu)}, mum, sigma)
  #
  @testset "matrix without missing" begin
    dsum = @inferred sum(dv, Symmetric(corrM))
    # checked with random numbers
    #boot_dvsums_acf(dv, acf1)
    @test params(dsum)[2] ≈ 0.128 rtol = 0.02
    @test mean(dsum) ≈ sum(mean.(dv))
    # acf variant
    dsum3 = @inferred sum(dv, acf1)
    @test all(params(dsum3) .≈ params(dsum))
    # mean function
    dmean = @inferred mean(dv, Symmetric(corrM))
    @test dmean.μ ≈ dsum.μ / length(dv)
    @test dmean.σ ≈ dsum.σ
    dmean_acf = @inferred mean(dv, acf1)
    @test dmean_acf == dmean
  end;
  @testset "matrix with missing" begin
    @test_throws ErrorException dsumm = sum(dvm, Symmetric(corrM))
    #S = similar(mum);
    dsumm = @inferred sum(dvm, Symmetric(corrM), SkipMissing())
    # checked from sample
    # boot_dvsums_acf(dvm[2:end], acf1)
    #@test σstar(dsumm) ≈ 1.15 rtol = 0.02
    @test params(dsumm)[2] ≈ 0.141 rtol = 0.02
    @test_throws ErrorException dsumm2 = sum(dvm, acf1)
    dsumm2 = @inferred sum(dvm, acf1, SkipMissing())
    @test all(params(dsumm2) .≈ params(dsumm))
  end;
  @testset "with gapfilling flag" begin
    isgapfilled = fill(false, length(dv)); isgapfilled[4:end] .= true
    dsum4 = sum(dv, Symmetric(corrM))
    dsum = @inferred sum(dv, Symmetric(corrM), isgapfilled=isgapfilled)
    #@code_warntype sum(dv, isgapfilled=isgapfilled)
    @test mean(dsum) == mean(dsum4)
    @test std(dsum) > std(dsum4)
  end;
  @testset "with missings and gapfilling flag" begin
    isgapfilled = fill(false, length(dvm)); isgapfilled[4:end] .= true
    dsum4 = sum(dvm, Symmetric(corrM), SkipMissing())
    dsum4b = @inferred sum(
      dvm, Symmetric(corrM), SkipMissing(), isgapfilled=isgapfilled)
    #@code_warntype sum(dvm, isgapfilled=isgapfilled)
    @test mean(dsum4b) == mean(dsum4)
    @test std(dsum4b) > std(dsum4)
    # acf variant
    dsum5b = @inferred sum(
      dvm, acf1, SkipMissing(), isgapfilled=isgapfilled)
    @test dsum5b == dsum4b
    dsum5c = @inferred sum(
      dvm, acf1, SkipMissing(), isgapfilled=isgapfilled,
      method = Val(:bandedmatrix))
    @test dsum5c == dsum4b
    @test_throws ErrorException sum(
      dvm, acf1, SkipMissing(), isgapfilled=isgapfilled,
      method = Val(:badmethodarg))
  end;
end;

end; # testset "sumlognormals"
# see also conde in sumlognormals_benchmark.jl
# could not stay here because of btime macro