text:
using Distributed
using JLD2

include("../BioClim/src/required.jl")

## Get & prepare data
@time @everywhere begin
    # Load data from CSV files
    #df = CSV.read("/Users/kiristern/Documents/GitHub/SDM/data/pred_prey/Puma_concolor.csv", header=true)
    # Prepare data (select columns, arrange values)
    # df = prepare_gbif_data(pred_prey)
    df = pred_prey
    # Separate species
    taxa_occ = [df[df.species .== u,:] for u in unique(df.species)]
end

## Load predictions
@load "../SDM/data/predictions_procyon_lotor.jld2" predictions
pred = predictions[1]
p1 = predictions[2]
p2 = predictions[3]
p3 = predictions[4]
p4 = predictions[5]
p5 = predictions[6]
p6 = predictions[7]
p7 = predictions[8]

# Combine layers
layers = [pred, p1, p2, p3,p4,p5,p6,p7]

# Create function
function expand_layers(layers::Array{SDMLayer{Float64},1})
    # Get minimum coordinates
    min_lon = min(map(x -> x.left, layers)...)
    max_lon = max(map(x -> x.right, layers)...)
    min_lat = min(map(x -> x.bottom, layers)...)
    max_lat = max(map(x -> x.top, layers)...)

    # Get grid size (rounding should prevent problems with last decimal)
    grid_size_lons = 1/round(1/stride(layers[1],1))
    grid_size_lats = 1/round(1/stride(layers[1],2))

    # Get coordinate range of newlayer -> original layers must have same stride
    lons_newlayers = min_lon+grid_size_lons:2*grid_size_lons:max_lon-grid_size_lons
    lats_newlayers = min_lat+grid_size_lats:2*grid_size_lats:max_lat-grid_size_lats

    # Create expanded layers
    newlayers = []
    for i in 1:length(layers)
        # Get position of original layers in new layer
        ## find where the point = 0 (= min value; all others are the absolute value)
        m_lon = findmin(abs.(layers[i].left+grid_size_lons .- lons_newlayers))[2]
        M_lon = findmin(abs.(layers[i].right-grid_size_lons .- lons_newlayers))[2]
        m_lat = findmin(abs.(layers[i].bottom+grid_size_lats .- lats_newlayers))[2]
        M_lat = findmin(abs.(layers[i].top-grid_size_lats .- lats_newlayers))[2]
        # Create new grid
        newgrid = fill(NaN, length(lats_newlayers), length(lons_newlayers))
        # Fill in original values
        newgrid[(m_lat:M_lat), (m_lon:M_lon)] .= layers[i].grid
        # Convert to SDMLayer
        newlayer = SDMLayer(newgrid, min_lon, max_lon, min_lat, max_lat)
        # Export result
        push!(newlayers, newlayer)
    end
    return newlayers
end

@time newlayers = expand_layers(layers)

#find probability of species occuring together
#create array, probs, using pred data
probs = copy(newlayers[1].grid)
#for all prey data
@time for j in 2:length(newlayers)
    for k in 1:length(newlayers[1].grid)
        #if prey occurence is not NaN value AND if the probability of overlap is not NaN (ie. no overlap)
        if !isnan(probs[k]) && !isnan(newlayers[j][k])
            #then multiply the probability of sp 1,2,3... by the next sp
            probs[k] = probs[k] .* newlayers[j][k]
        end
    end
end
probs

#check if probability of occurence changed (thus overlap in species occurence) -- if layer[1] value changes or not
presence = probs .== newlayers[1].grid
#if no change, change to NaN
probs[presence] .= NaN
#filter only occurence data (ie. not NaN values)
filter(!isnan, probs)

# Convert to SDMLayer
map_layer = SDMLayer(probs, newlayers[1].left, newlayers[1].right, newlayers[1].bottom, newlayers[1].top)

plot_layer = plotSDM(map_layer)
plot_layer = plot!(title="Probability of Occurence: Procyon lotor")

savefig("../SDM/figures/sdm_map-Procyon_lotor_predictions_title.pdf")

##save individual maps
#map preds
worldmap = clip(worldshape(50), pred)

sdm_plot = plot([0.0], lab = " ", msw = 0.0, ms = 0.0, size = (900, 450), frame=:box,
                title = "Procyon lotor")
xaxis!(sdm_plot, (pred.left, pred.right), "Longitude")
yaxis!(sdm_plot, (pred.bottom, pred.top), "Latitude")

for p in worldmap
    sh = Shape([pp.x for pp in p.points], [pp.y for pp in p.points])
    plot!(sdm_plot, sh, c = :lightgrey, lab = " ")
end

heatmap!(
    sdm_plot,
    longitudes(pred), latitudes(pred), pred.grid,
    aspectratio = 92.60/60.75, c = :BuPu,
    clim = (0.0, maximum(filter(!isnan, pred.grid)))
    )

for p in worldmap
    xy = map(x -> (x.x, x.y), p.points)
    plot!(sdm_plot, xy, c=:grey, lab = " ", leg=false)
end
sdm_plot

savefig("../SDM/figures/sdm_map-Procyon_lotor.pdf")

#map prey[x]
worldmap = clip(worldshape(50), p7)

sdm_plot = plot([0.0], lab = " ", msw = 0.0, ms = 0.0, size = (900, 450), frame=:box,
                title = "Anas platyrhynchos")
xaxis!(sdm_plot, (p7.left, p7.right), "Longitude")
yaxis!(sdm_plot, (p7.bottom, p7.top), "Latitude")

for p in worldmap
    sh = Shape([pp.x for pp in p.points], [pp.y for pp in p.points])
    plot!(sdm_plot, sh, c = :lightgrey, lab = " ")
end

heatmap!(
    sdm_plot,
    longitudes(p7), latitudes(p7), p7.grid,
    aspectratio = 92.60/60.75, c = :BuPu,
    clim = (0.0, maximum(filter(!isnan, p7.grid)))
    )

for p in worldmap
    xy = map(x -> (x.x, x.y), p.points)
    plot!(sdm_plot, xy, c=:grey, lab = " ", leg=false)
end
sdm_plot

savefig("../SDM/figures/sdm_map-Anas_platyrhynchos.pdf")

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
"""
    Module with an implementation of CORK (CORK version 0.2)
     http://twr.cs.kuleuven.be/research/software/nleps/cork.php
"""
module NEPSolver_CORK # Will change later

    using LinearAlgebra

    export cork
    export TargetType
    export CORKInputType

    type TargetType
        Sigma::Complex64
        shift::Complex64
    end

    type CORKInputType
        A # Length d array with n x n matrices
        B # Length d array with n x n matrices
        M # (d-1) x d matrix
        N # (d-1) x d matrix

    end
    """

"""
#    function cork(L)
#
#    end
#
#    function cork(L,k::Integer)
#
#    end
#
#    function cork(L,k::Integer,sigma::AbstractFloat)
#
#    end

    function cork(L::CORKInputType,k::Integer,target::TargetType)

        #
        n=size(L.A[1],1)

        ##  Check inputs
        #nb = ceil( (m + maxrest*(m-p) + 1) / floor(length(shifts)) );

        nb=10; # hard coded
        shifts=repmat([target.shift],nb)

        ## Initialize

        # Default values:

        m = max(3k,20);
        maxrest = 50;
        p = max(2k,10);

        d=size(L.A,1)

        ## Initialization
        v0=ones(n)

        # Compute the Q-matrix
        Q,U0 = qr(reshape(v0/opnorm(v0),n,1));
        Q = Matrix(Q)
        r = size(Q,2);

        # Init tensor U
        #U = zeros(r,m+1,d);
        U = zeros(r+1,m+1,d); #Don't use size(U) since it is different from CORK
        U[1:r,1,1:size(U0,1)] = U0;
        j = 1;

        # Init Hessenberg matrices H and K
        H = zeros(m+1,m);
        K = zeros(m+1,m);

        # Storage of eigenpairs and residuals
        X = NaN*zeros(n,m);
        lambda = NaN*zeros(m);
        res = NaN*zeros(m);

        # History storing
        Lam = NaN*zeros(m,m + maxrest*(m-p));
        Res = NaN*zeros(m,m + maxrest*(m-p));

        # What's this for ??
        J = [1;zeros(m + maxrest*(m-p),1)];
        R = [r;zeros(m + maxrest*(m-p),1)];

        # For LU
        LU = Array{Any}(1);
        SHIFTS = [];

        # Don't know how to do cellfun in julia

        # Which matrices in L.A and L.B have non-zero elements
        NNZ_A=zeros(Bool,d)
        NNZ_B=zeros(Bool,d)
        for i=1:d
            NNZ_A[i]=countnz(L.A[i])>0
            NNZ_B[i]=countnz(L.B[i])>0
        end

        i=1
        nbrest = 0;
        l = 0;
        count = 0;
        flag = true;

        println("shifts:",shifts)
        while i <= m + maxrest*(m-p)

            println("size(Q):",size(Q))
            println("r:",r," i=",i," j=",j, " Q=",Q)
            println("U=",U)
            r = corkstep!(shifts[i],r,j,i,Q,U,L,NNZ_A,NNZ_B)
            println("r=",r)
            return
            # r = corkstep(shifts(i),r,j,i);

            #            flag = ritzpairs(r,j,i,l);
            #            if j==m
            #                if nbrest < maxrest
            #                    nbrest = nbrest + 1;
            #                    [r,j,l] = implicitrestart(r,nbrest);
            #                else
            #                    break
            #                end
            #            end

            i+=1
            j+=1
        end

    end

    function corkstep!(shift,r,j,i,Q,U,L,NNZ_A,NNZ_B)
        println("Iteration ",i)
        eta=1/sqrt(2);
        d=size(L.A,1)

        # Next vector
        v=Q[:,1:r]*reshape(U[1:r,j,1:d],r,d);

        (v1,MsN1inv0,MsN1invN)=backslash(shift,v,L,NNZ_A,NNZ_B)
        # Not yet implemented

        # level 1 orthogonalization
        q = copy(v1);
        delta = Inf; nborth = 0;
        maxreorth=2
        while (opnorm(q) < eta*delta) && (nborth <= maxreorth)
            delta = opnorm(q);
            if nborth == 0
                u1 = Q[:,1:r]'*v1;
                q = q - Q[:,1:r]*u1;
            else
                q = q - Q[:,1:r]*(Q[:,1:r]'*q);
            end
            nborth = nborth + 1;
            println("Reorth Q ",nborth)

        end
        delta = opnorm(q)+eps();
        #println("XXXXXXXXXXX")

        println("delta=",delta)
        # Update Q
        if delta > eps()
            rnew = r + 1;
            q = q/delta;
            #Q[:,rnew] = q;
            Q=hcat(Q,q)
            println("size(Q)=",size(Q))
        else
            rnew = r;
        end
        println("rnew=",rnew)
        # level 2 orthogonalization
        if rnew > r
            ## How to expand tensors? (Workaround: init with larger)
            U[rnew,:,:] = 0;
            u1 = [u1;q'*v1];
        end
        tmp=broadcast(*,MsN1inv0.',u1)

        println("size(u1)=",size(u1))
        u2=reshape(U[1:rnew,j,1:d],rnew,d)*MsN1invN.'-tmp

        println("*************** tmp=", tmp)
        println(" u2 =",u2)
        ##
        u = [u1;reshape(u2,size(u2,1)*size(u2,2),1)];
        println(" opnorm(u) =",opnorm(u))
        while (opnorm(u) < eta*delta) && (nborth <= maxreorth)
            delta = opnorm(u);
            h = reshape(conj(permute(U(1:rnew,1:j,1:d),[2,1,3])),[],d*rnew)*u;

        end
        rnew=1
        return rnew,Q

    end
    function backslash(shift,y,L,NNZ_A,NNZ_B)

        reuselu=false
        newlu=true

        # initialize
        d=size(L.A,1)
        m0 = L.M[1:d-1,1];
        M1 = L.M[1:d-1,2:d];
        n0 = L.N[1:d-1,1];
        N1 = L.N[1:d-1,2:d];

        # Kronecker coefficients
        MsN1 = M1 - shift*N1;
        msn0 = m0 - shift*n0;
        MsN1inv0 = MsN1\msn0;
        MsN1invN = MsN1\L.N[1:d-1,1:d];

        ## Build At
        At = L.A[1] - shift*L.B[1];
        for ii = 2:d
            At = At - MsN1inv0[ii-1]*(L.A[ii] - shift*L.B[ii]);
            println("MsN1inv0[ii-1]=",MsN1inv0[ii-1])
        end

        println("opnorm(At)=",opnorm(At))

        println(d)

        ## intermediate x used in yt
        x = y*MsN1invN.';

        ## build yt
        yt = L.B[1]*y[:,1];
        for ii = 2:d
            if NNZ_A[ii] && NNZ_B[ii]
                yt = yt - L.A[ii]*x[:,ii-1] + L.B[ii]*(shift*x[:,ii-1]+y[:,ii]);
            elseif NNZ.A[ii]
                yt = yt - L.A[ii]*x[:,ii-1];
            else
                yt = yt + L.B[ii]*(shift*x[:,ii-1]+y[:,ii]);
            end
        end

        ## build x1
        println("At",size(At))
        x1 = At\yt;

        return x1,MsN1inv0,MsN1invN

    end
end

ext: jl


text:
using myJuliaUtils
using Test
using LinearAlgebra
using Statistics

m = 0
reps = 5
for i = 1:reps
    K = round(Int, 1 + rand() * 9);
    uV = randn(K, K);
    uV = uV'uV
    epsilon = cholesky(uV).L * randn(K, 50000000);
    global m += cov(epsilon') ≈ uV
end
@test m == reps

max_stars_repo_name: zymbuzz/myJuliaUtils

max_issues_repo_name: zymbuzz/myJuliaUtils

max_forks_repo_name: zymbuzz/myJuliaUtils


text:
@testset "`validate_date`" begin
    @test validate_date("6/1/2017") == "6/1/2017"
    @test validate_date("06/01/2017") == "6/1/2017"
    @test validate_date("6-1-2017") == "6/1/2017" # i guess we'll handle dashes
    @test validate_date(Date(2017, 6, 1)) == "6/1/2017"

    @test_throws ArgumentError validate_date("5/32/2017") # invalid date rounds to next month
    @test_throws ArgumentError validate_date("6/1/17") # full year required
    @test_throws ArgumentError validate_date("June 1 2017") # no spelled-out parsing handled
    @test_throws ArgumentError validate_date("1 June 2017") # no spelled-out parsing handled
end

@testset "`request_meetings`" begin
    meetings = request_meetings("6/1/2020", "6/1/2020")
    @test nrow(meetings) == 2
    m = first(meetings)
    r = Meeting(m) # Just test that this is possible
    @test m.name == "Public Health and Public Safety Committee"
    @test m.date == DateTime("2020-06-01T18:00:00")
    @test isa(m.date, DateTime)
    @test isa(m.id, Int)
    @test m.link == "http://somervillecityma.iqm2.com/Citizens/Detail_Meeting.aspx?ID=3163"

    # Test that current meeting serialization hasn't changed
    # If this test _fails_, need to (a) bump `meeting_version` and (b) resave current
    # version of cached meeting asset:
    test_meeting_path = joinpath(TEST_ASSETS, "v$(SomervilleCouncilParser.meeting_version)",
                                 "meeting.arrow")
    # Legolas.write(test_meeting_path, DataFrame(m), SomervilleCouncilParser.meeting_schema)  # Uncomment to save new version
    @test isfile(test_meeting_path)
    previous_mtg = DataFrame(Legolas.read(test_meeting_path))
    @test previous_mtg == DataFrame(m)
end


text:
const JOB_QUEUE = Ref(Dict{String, Tuple{Int, Jobs.JobState}}())

function main_loop(s::Server)
    job_dirs_procs = Dict{String,Task}()
    queue!(JOB_QUEUE[], s, true)
    for (j, info) in JOB_QUEUE[]
        if info[2] == Jobs.Running
            job_dirs_procs[j] = spawn_worker(j)
        end
    end
    while true
        queue!(JOB_QUEUE[], s)
        handle_workflow_runners!(job_dirs_procs)
        if length(job_dirs_procs) < s.max_concurrent_jobs
            handle_job_submission!(s, job_dirs_procs)
        end
        sleep(SLEEP_TIME)
    end
end

function handle_workflow_runners!(job_dirs_procs)
    to_rm = String[]
    for (k, t) in job_dirs_procs
        if istaskdone(t)
            if istaskfailed(t)
                @error "Task failed with $(t.result)"
            else
                t_ = fetch(t)
                @info """Workflow for job directory $(k) done."""
            end
            push!(to_rm, k)
        end
    end
    if !isempty(to_rm)
        pop!.((job_dirs_procs,), to_rm)
    end
end

function spawn_worker(jobdir)
    # TODO: implement workflows again
    # if ispath(joinpath(jobdir, ".workflow/environment.jld2"))
    #     env_dat = JLD2.load(joinpath(jobdir, ".workflow/environment.jld2"))
    #     proc = addprocs(1; exeflags = "--project=$(env_dat["project"])")[1]
    #     using_expr = Base.Meta.parse("""using $(join(env_dat["modules"], ", "))""")
    #     ctx_path = joinpath(jobdir, ".workflow/ctx.jld2")

    #     Distributed.remotecall_eval(Distributed.Main, proc, using_expr)
        # f = remotecall(DFControl.run_queue, proc, job, JLD2.load(ctx_path)["ctx"];
                       # sleep_time = SLEEP_TIME)
        # return proc, f
    # else
    id = next_jobid()
    return Threads.@spawn begin
        with_logger(job_logger(id)) do
            @info (timestamp = Dates.now(), jobdir = jobdir)
            sleep(3*SLEEP_TIME)
            info = state(jobdir)
            while info == Jobs.Pending || info == Jobs.Running || info == Jobs.Submitted
                sleep(SLEEP_TIME)
                info = state(jobdir)
                @info (timestamp = Dates.now(), jobdir = jobdir, state = info)
            end
            outputdata(jobdir, map(x -> splitext(splitpath(x.infile)[end])[1], FileIO.read_job_script(joinpath(jobdir, "job.tt"))[2]))
            @info (timestamp = Dates.now(), jobdir = jobdir, state = info)
        end
    end
end

# Jobs are submitted by the daemon, using supplied job jld2 from the caller (i.e. another machine)
# Additional files are packaged with the job
function handle_job_submission!(s::Server, job_dirs_procs)
    s = DFC.Server("localhost")
    lines = filter(!isempty, readlines(PENDING_JOBS_FILE))
    write(PENDING_JOBS_FILE, "")
    if length(lines) + length(job_dirs_procs) > s.max_concurrent_jobs

        to_submit = lines[1:s.max_concurrent_jobs - length(job_dirs_procs)]
        open(PENDING_JOBS_FILE, "a") do f
            for l in lines[s.max_concurrent_jobs - length(job_dirs_procs) + 1:end]
                write(f, l * "\n")
            end
        end
    else
        to_submit = lines
    end

    if !isempty(to_submit)
        curdir = pwd()
        while !isempty(to_submit)
            j = pop!(to_submit)
            try
                info = submit(s, j)
                JOB_QUEUE[][j] = info
                job_dirs_procs[j] = spawn_worker(j)
            catch
                sleep(SLEEP_TIME)
                push!(to_submit, j)
            end
        end
        cd(curdir)
    end
end

function run_queue(job::Job, ctx::Dict; sleep_time = 10)
    logger = workflow_logger(job)
    with_logger(logger) do
        qd = queued_dir(job)
        queued_steps = readdir(qd)
        order = sortperm(parse.(Int, getindex.(splitext.(queued_steps), 1)))
        fd = finished_dir(job)
        if !ispath(qd) || isempty(queued_steps)
            return
        end
        if !ispath(fd)
            mkpath(fd)
        else
            for f in readdir(fd)
                rm(joinpath(fd, f))
            end
        end

        for f in queued_steps[order]
            step_file = joinpath(qd, f)
            t = include(step_file)
            @eval $(t)($(job), $(ctx))
            JLD2.jldsave(joinpath(job, ".workflow/ctx.jld2"); ctx=ctx)
            while isrunning(job)
                sleep(sleep_time)
            end
            Service.outputdata(job)
            mv(step_file, joinpath(fd, f); force = true)
        end
    end
    return true
end


text:
export black!, red!, green!, yellow!, blue!, magenta!, cyan!, white!
export color

COLORS = Dict(
"black"         => 0,
"red"           => 1,
"green"         => 2,
"yellow"        => 3,
"blue"          => 4,
"magenta"       => 5,
"cyan"          => 6,
"white"         => 7,
"default"       => 9
)

MODES = Dict(
"default"       => 0,
"bold"          => 1,
"underline"     => 4,
"blink"         => 5,
"swap"          => 7,
"hide"          => 8
)

function color(color::String, str::String; background::String="default", mode::String="default")
    FColor = "\e[$(30 + COLORS[color])m"
    BColor = "\e[$(40 + COLORS[background])m"
    Action = "\e[$(MODES[mode])m"
    return Action * BColor * FColor * str * "\e[0m"
end

black!(str; mode::String="default")   = color("black",   string(str), mode=mode)
red!(str; mode::String="default")     = color("red",     string(str), mode=mode)
green!(str; mode::String="default")   = color("green",   string(str), mode=mode)
yellow!(str; mode::String="default")  = color("yellow",  string(str), mode=mode)
blue!(str; mode::String="default")    = color("blue",    string(str), mode=mode)
magenta!(str; mode::String="default") = color("magenta", string(str), mode=mode)
cyan!(str; mode::String="default")    = color("cyan",    string(str), mode=mode)
white!(str; mode::String="default")   = color("white",   string(str), mode=mode)

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
#!/usr/bin/env julia

using Yao
using Yao.ConstGate # needed for P1 = 0.5*(I - sigma_z) block
using Plots
using BitBasis
using StatsBase
using GraphRecipes
using GRUtils

#=
H(t) = Ω(t) ∑_i σ_i^x - δ(t) ∑_i n_i + u ∑_ij n_i n_j
=#

const u = 1.35
const Ω_max = 1.89
const δ_0 = -1.0
const δ_max = 1.0

function get_edges(graph::Vector{NTuple{2, Float64}})
    Nv = size(graph)[1]
    edges = falses(Nv, Nv)
    for i in 1:(Nv-1)
        xi, yi = graph[i]
        for j in (i+1):Nv
            xj, yj = graph[j]

            dij = sqrt((xi - xj)^2. + (yi - yj)^2.)
            if dij <= 1.0
                edges[i,j] = true
            end
        end
    end
    return findall(edges)
end

function Ω(t::Float64)
    if 0 <= t <= 0.25
        return (Ω_max / 0.25) * t
    elseif 0.25 < t <= 0.69
        return Ω_max
    elseif 0.69 < t <= 1
        return - Ω_max * t / 0.31 + Ω_max * (1 + 0.69/0.31)
    end
end

function δ(t::Float64)
    slope = (δ_0 - δ_max)/(0.25 - 0.69)
    if 0 <= t <= 0.25
        return δ_0
    elseif 0.25 < t <= 0.69
        return t * slope + (δ_max - slope * 0.69)
    elseif 0.69 < t <= 1
        return δ_max
    end
end

function hamiltonian(graph::Vector{NTuple{2, Float64}}, edges::Vector{CartesianIndex{2}}, t::Float64)
    # the UD-MIS Hamiltonian
    Nv = size(graph)[1] # number of vertices

    interaction_term = map(1:size(edges)[1]) do i
        l,m = edges[i][1], edges[i][2]
        repeat(Nv,u*P1,(l,m))
    end |> sum
    interaction_term - δ(t)*sum(map(i->put(Nv,i=>P1), 1:Nv)) + Ω(t)*sum(map(i->put(Nv,i=>X), 1:Nv))
end

function run_annealing(graph::Vector{NTuple{2, Float64}}, edges::Vector{CartesianIndex{2}}, dt::Float64)
    psi_t = zero_state(size(graph)[1])
    for t in 0:dt:1.0
        h = hamiltonian(graph, edges, t)
        psi_t = psi_t |> TimeEvolution(h, dt * 100)
    end
    return psi_t
end

graph = [(0.3461717838632017, 1.4984640297338632),
         (0.6316400411846113, 2.5754677320579895),
         (1.3906262250927481, 2.164978861396621),
         (0.66436005100802, 0.6717919819739032),
         (0.8663329771713457, 3.3876341010035995),
         (1.1643107343501296, 1.0823066243402013)
        ]
edges = get_edges(graph)
dt = 0.001

psi = run_annealing(graph, edges, dt)
open("task2_data.dat","w") do io
    for sample in measure(psi; nshots=10_000)
        println(io, sample)
    end
end

# added histogram to show degeneration of the solutions
nshots = 1000
samples = measure(psi; nshots=nshots)
# @show samples


text:
module AMQPTestThroughput

using AMQPClient, Test, Random

const JULIA_HOME = Sys.BINDIR

const EXCG_DIRECT = "amq.direct"
const QUEUE1 = "queue1"
const ROUTE1 = "key1"
const MSG_SIZE = 1024
const NMSGS = 10^5
const no_ack = true

const M = Message(rand(UInt8, 1024), content_type="application/octet-stream", delivery_mode=PERSISTENT)

testlog(msg) = println(msg)

function setup(;virtualhost="/", host="localhost", port=AMQPClient.AMQP_DEFAULT_PORT, auth_params=AMQPClient.DEFAULT_AUTH_PARAMS)
    # open a connection
    testlog("opening connection...")
    conn = connection(;virtualhost=virtualhost, host=host, port=port, auth_params=auth_params)

    # open a channel
    testlog("opening channel...")
    chan1 = channel(conn, AMQPClient.UNUSED_CHANNEL, true)
    @test chan1.id == 1

    # create and bind queues
    testlog("creating queues...")
    success, name, message_count, consumer_count = queue_declare(chan1, QUEUE1)
    @test success
    @test message_count == 0

    @test queue_bind(chan1, QUEUE1, EXCG_DIRECT, ROUTE1)

    conn, chan1
end

function teardown(conn, chan1, delete=false)
    testlog("closing down...")
    if delete
        success, message_count = queue_purge(chan1, QUEUE1)
        @test success
        @test message_count == 0

        @test queue_unbind(chan1, QUEUE1, EXCG_DIRECT, ROUTE1)
        success, message_count = queue_delete(chan1, QUEUE1)
        @test success
        @test message_count == 0
    end

    # close channels and connection
    close(chan1)
    AMQPClient.wait_for_state(chan1, AMQPClient.CONN_STATE_CLOSED)
    @test !isopen(chan1)

    close(conn)
    AMQPClient.wait_for_state(conn, AMQPClient.CONN_STATE_CLOSED)
    @test !isopen(conn)
end

function publish(conn, chan1)
    testlog("starting basic publisher...")
    # publish N messages
    for idx in 1:NMSGS
        basic_publish(chan1, M; exchange=EXCG_DIRECT, routing_key=ROUTE1)
        if (idx % 10000) == 0
            println("publishing $idx ...")
            sleep(1)
        end
    end
end

function consume(conn, chan1)
    testlog("starting basic consumer...")
    # start a consumer task
    global msg_count = 0
    global start_time = time()
    global end_time = 0
    consumer_fn = (rcvd_msg) -> begin
        global msg_count
        global end_time
        msg_count += 1
        if ((msg_count % 10000) == 0) || (msg_count == NMSGS)
            #basic_ack(chan1, 0; all_upto=true)
            println("ack sent $msg_count ...")
        end
        no_ack || basic_ack(chan1, rcvd_msg.delivery_tag)
        if msg_count == NMSGS
            end_time = time()
        end
    end
    success, consumer_tag = basic_consume(chan1, QUEUE1, consumer_fn; no_ack=no_ack)
    @test success

    # wait to receive all messages
    while msg_count < NMSGS
        println("$msg_count of $NMSGS messages processed")
        sleep(2)
    end

    # cancel the consumer task
    @test basic_cancel(chan1, consumer_tag)

    # time to send and receive
    total_time = max(end_time - start_time, 1)
    println("time to send and receive $NMSGS messages: $(end_time - start_time) secs @ $(NMSGS/total_time) msgs per second")
end

function run_publisher()
    conn, chan1 = AMQPTestThroughput.setup()
    AMQPTestThroughput.publish(conn, chan1)
    AMQPTestThroughput.teardown(conn, chan1, false) # exit without destroying queue
    nothing
end

function run_consumer()
    conn, chan1 = AMQPTestThroughput.setup()
    AMQPTestThroughput.consume(conn, chan1)
    println("waiting for publisher to exit gracefully...")
    sleep(10)  # wait for publisher to exit gracefully
    AMQPTestThroughput.teardown(conn, chan1, true)
    nothing
end

function spawn_test(script, flags)
    opts = Base.JLOptions()
    inline_flag = opts.can_inline == 1 ? `` : `--inline=no`
    cov_flag = (opts.code_coverage == 1) ? `--code-coverage=user` :
                 (opts.code_coverage == 2) ? `--code-coverage=all` :
                 ``
    srvrscript = joinpath(dirname(@__FILE__), script)
    srvrcmd = `$(joinpath(JULIA_HOME, "julia")) $cov_flag $inline_flag $srvrscript $flags`
    println("Running tests from ", script, "\n", "="^60)
    ret = run(srvrcmd)
    println("Finished ", script, "\n", "="^60)
    nothing
end

function runtests()
    println("starting consumer")
    consumer = @async spawn_test("test_throughput.jl", "--runconsumer")
    sleep(5)
    println("starting publisher")
    publisher = @async spawn_test("test_throughput.jl", "--runpublisher")
    wait(consumer)
    wait(publisher)
    nothing
end

end # module AMQPTestThroughput

!isempty(ARGS) && (ARGS[1] == "--runpublisher") && AMQPTestThroughput.run_publisher()
!isempty(ARGS) && (ARGS[1] == "--runconsumer")  && AMQPTestThroughput.run_consumer()


text:
using Flux
using WindMixing
using JLD2
using FileIO
using OceanParameterizations
using OrdinaryDiffEq
using BenchmarkTools

train_files = ["-1e-3"]

𝒟train = data(train_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)

PATH = pwd()

OUTPUT_PATH = "D:\\University Matters\\Massachusetts Institute of Technology\\CLiMA Project\\OceanParameterizations.jl\\training_output"
FILE_PATH = joinpath(OUTPUT_PATH, "test_nonmutating.jld2")

FILE_PATH_NN = joinpath(PATH, "extracted_training_output", "NDE_training_convective_adjustment_1sim_-1e-3_2_extracted.jld2")

file = jldopen(FILE_PATH_NN, "r")
uw_NN = file["neural_network/uw"]
vw_NN = file["neural_network/vw"]
wT_NN = file["neural_network/wT"]

train_epochs = [10]
train_tranges = [1:40:1153]
train_optimizers = [[ADAM(0.01)] for i in 1:length(train_epochs)]
timestepper = ROCK4()

benchmark_mutating = []
benchmark_nonmutating = []

tranges = [1:10:100, 1:20:200, 1:10:200, 1:20:400, 1:40:400]

for trange in tranges
    @info "$(trange), mutating"
    t = @benchmark train_NDE_convective_adjustment(uw_NN, vw_NN, wT_NN, 𝒟train, $trange, ROCK4(), [ADAM(0.01)], 1, FILE_PATH, 1, 1, 10f0, 1)
    push!(benchmark_mutating, t)
    @info "$(trange), non-mutating"
    t = @benchmark train_NDE_convective_adjustment_nonmutating(uw_NN, vw_NN, wT_NN, 𝒟train, $trange, ROCK4(), [ADAM(0.01)], 1, FILE_PATH, 1, 1, 10f0, 1)
    push!(benchmark_nonmutating, t)
end

@info minimum.(benchmark_mutating)
@info minimum.(benchmark_nonmutating)


text:
using Base.Test
using LJH
# Write file 1 as an LJH 2.1 file and files 2 as LJH 2.2 with identical data.
name21, f1 = mktemp()
name22, f2 = mktemp()

dt = 9.6e-6
npre = 200
nsamp = 1000

LJH.writeljhheader(f1, dt, npre, nsamp; version="2.1.0")
LJH.writeljhheader(f2, dt, npre, nsamp; version="2.2.0")

rowcount = collect(5000:1000:10000)
timestamps = [Int64(round(r*(dt/30)*1e6)) for r in rowcount]
N = length(rowcount)

data = rand(UInt16, nsamp, N)
data[1,:] = 0xffff

LJH.writeljhdata(f1, data, rowcount)
LJH.writeljhdata(f2, data, rowcount, timestamps)

close(f1)
close(f2)
# Now check that the data are readable
ljh21 = LJHGroup(name21)
ljh22 = LJHGroup(name22)
# Test that the header info is correct
for ljh in (ljh21, ljh22)
    @test record_nsamples(ljh) == nsamp
    @test pretrig_nsamples(ljh) == npre
    @test frametime(ljh) == dt
    @test length(ljh) == N

    # Test indexed access
    ranges = (1:N, N:-1:1, 1:3:N, N:-2:1)
    for r in ranges
        for i = r
            record = ljh[i]
            @test record.data==data[:,i]
            @test record.timestamp_usec==timestamps[i]
        end
    end
    row(ljh)
    column(ljh)
    num_rows(ljh)
    num_columns(ljh)

    # test slices
    for (r1,r2) in zip(collect(ljh[2:4]), collect(ljh)[2:4])
        @test r1.data==r2.data
    end
end

# Now a group corresponding to 3 files (actually, same one 3 times)
grp = LJHGroup([name22, name21, name22])
for j=1:3N
    record = grp[j]
    d,r,t = record.data, record.rowcount, record.timestamp_usec
    i = mod(j-1, N)+1
    @test d==data[:,i]
    @test t==timestamps[i]
end

for (r1,r2) in zip(collect(grp[3:10]), collect(grp)[3:10])
    @test r1.data == r2.data
    @test r1.rowcount == r2.rowcount
    @test r1.timestamp_usec == r2.timestamp_usec
end
@test lengths(grp) == [N,N,N]

grp.ljhfiles[1].record_nsamples=0
@test_throws AssertionError record_nsamples(grp)

data_r, rowcount_r, timestamp_usec_r = get_data_rowcount_timestamp(ljh22)
for i = 1:N
  @assert data[:,1]==data_r[1]
end
@assert rowcount_r==rowcount
@assert timestamp_usec_r == timestamps


text:
@testset "Basic CSV.Source" begin

#test on non-existent file
@test_throws ArgumentError CSV.Source("");

#test where datarow > headerrow
@test_throws ArgumentError CSV.Source(joinpath(dir, "test_no_header.csv");datarow=1,header=2);

#test various encodings
f = CSV.Source(joinpath(dir, "test_utf8_with_BOM.csv"))
@test Data.header(Data.schema(f)) == ["col1","col2","col3"]

f = CSV.Source(joinpath(dir, "test_utf8.csv"), allowmissing=:auto)
sch = Data.schema(f)
@test f.options.delim == UInt8(',')
@test size(sch, 2) == 3
@test size(sch, 1) == 3
@test Data.header(sch) == ["col1","col2","col3"]
@test Data.types(sch) == (Float64,Float64,Float64)
ds = CSV.read(f)
@test ds[1][1] == 1.0
@test ds[1][2] == 4.0
@test ds[1][3] == 7.0
@test ds[2][1] == 2.0
sch2 = Data.schema(ds)
@test Data.header(sch2) == Data.header(sch)
@test Data.types(sch) == (Float64,Float64,Float64)

f = CSV.Source(joinpath(dir, "test_utf8.csv"))
si = CSV.write(joinpath(dir, "new_test_utf8.csv"), f)
so = CSV.Source(si)
@test so.options.delim == UInt8(',')
@test size(Data.schema(so), 2) == 3
@test size(Data.schema(so), 1) == 3
@test Data.header(Data.schema(so)) == ["col1","col2","col3"]
@test Data.types(Data.schema(so)) == (Union{Float64,Missing},Union{Float64,Missing},Union{Float64,Missing})
ds = CSV.read(so)
@test ds[1][1] == 1.0
@test ds[1][2] == 4.0
@test ds[1][3] == 7.0
@test ds[2][1] == 2.0
@test Data.types(Data.schema(f)) == Data.types(Data.schema(so)) == Data.types(Data.schema(ds))
f = si = so = ds = nothing; GC.gc(); GC.gc()
try
rm(joinpath(dir, "new_test_utf8.csv"))
catch e
end

# f = CSV.Source(joinpath(dir, "test_utf16_be.csv"))
# f = CSV.Source(joinpath(dir, "test_utf16_le.csv"))
# f = CSV.Source(joinpath(dir, "test_utf16.csv"))
f = CSV.Source(joinpath(dir, "test_windows.csv"))
@test size(Data.schema(f), 2) == 3
@test size(Data.schema(f), 1) == 3

#test one column file
f = CSV.Source(joinpath(dir, "test_single_column.csv"), allowmissing=:auto)
@test f.options.delim == UInt8(',')
@test f.options.quotechar == UInt8('"')
@test f.options.escapechar == UInt8('\\')
@test size(Data.schema(f), 2) == 1
@test size(Data.schema(f), 1) == 3
@test Data.header(Data.schema(f)) == ["col1"]
@test Data.types(Data.schema(f)) == (Int64,)
ds = CSV.read(f)
@test ds[1][1] == 1
@test ds[1][2] == 2
@test ds[1][3] == 3

#test empty file
f = CSV.Source(joinpath(dir, "test_empty_file.csv"))
@test ismissing(size(Data.schema(f), 1))

#test file with just newlines
f = CSV.Source(joinpath(dir, "test_empty_file_newlines.csv"))
@test size(Data.schema(f), 2) == 1
@test size(Data.schema(f), 1) == 9
@test Data.header(Data.schema(f)) == [""]
@test Data.types(Data.schema(f)) == (Missing,)

#test with various quotechars, escapechars
f = CSV.Source(joinpath(dir, "test_simple_quoted.csv"))
@test size(Data.schema(f), 2) == 2
@test size(Data.schema(f), 1) == 1
ds = CSV.read(f)
@test String(ds[1][1]) == "quoted field 1"
@test String(ds[2][1]) == "quoted field 2"
f = CSV.Source(joinpath(dir, "test_quoted_delim_and_newline.csv"))
@test size(Data.schema(f), 2) == 2
@test size(Data.schema(f), 1) == 1

f = CSV.Source(joinpath(dir, "test_quoted_numbers.csv"); categorical=false, strings=:weakref, allowmissing=:auto)
@test size(Data.schema(f), 2) == 3
@test size(Data.schema(f), 1) == 3
ds = CSV.read(f)
@test Data.types(Data.schema(f)) == (WeakRefString{UInt8}, Int64, Int64)

#test various newlines
f = CSV.Source(joinpath(dir, "test_crlf_line_endings.csv"), allowmissing=:auto)
@test Data.header(Data.schema(f)) == ["col1","col2","col3"]
@test size(Data.schema(f), 2) == 3
@test Data.types(Data.schema(f)) == (Int64,Int64,Int64)
ds = CSV.read(f)
@test ds[1][1] == 1
f = CSV.Source(joinpath(dir, "test_newline_line_endings.csv"), allowmissing=:auto)
@test Data.header(Data.schema(f)) == ["col1","col2","col3"]
@test size(Data.schema(f), 2) == 3
@test Data.types(Data.schema(f)) == (Int64,Int64,Int64)
f = CSV.Source(joinpath(dir, "test_mac_line_endings.csv"), allowmissing=:auto)
@test Data.header(Data.schema(f)) == ["col1","col2","col3"]
@test size(Data.schema(f), 2) == 3
@test Data.types(Data.schema(f)) == (Int64,Int64,Int64)

end # testset

@testset "CSV.Source keyword arguments" begin

#test headerrow, datarow, footerskips
f = CSV.Source(joinpath(dir, "test_no_header.csv"); header=0, datarow=1, allowmissing=:auto)
@test Data.header(Data.schema(f)) == ["Column1","Column2","Column3"]
@test Data.types(Data.schema(f)) == (Float64,Float64,Float64)
@test size(Data.schema(f), 2) == 3
@test size(Data.schema(f), 1) == 3
f = CSV.Source(joinpath(dir, "test_2_footer_rows.csv"); header=4, datarow=5, footerskip=2, allowmissing=:auto)
@test Data.header(Data.schema(f)) == ["col1","col2","col3"]
@test Data.types(Data.schema(f)) == (Int64,Int64,Int64)
@test size(Data.schema(f), 2) == 3
@test size(Data.schema(f), 1) == 3

#test dates, dateformats
f = CSV.Source(joinpath(dir, "test_dates.csv"); types=[Date], dateformat="yyyy-mm-dd", allowmissing=:auto)
@test size(Data.schema(f), 2) == 1
@test size(Data.schema(f), 1) == 3
@test Data.types(Data.schema(f)) == (Date,)
ds = CSV.read(f)
@test ds[1][1] == Date(2015,1,1)
@test ds[1][2] == Date(2015,1,2)
@test ds[1][3] == Date(2015,1,3)
f = CSV.Source(joinpath(dir, "test_excel_date_formats.csv"); dateformat="mm/dd/yy", allowmissing=:auto)
@test size(Data.schema(f), 2) == 1
@test size(Data.schema(f), 1) == 3
@test Data.types(Data.schema(f)) == (Date,)
f = CSV.Source(joinpath(dir, "test_excel_date_formats.csv"); types=[Date], dateformat="mm/dd/yy", allowmissing=:auto)
@test size(Data.schema(f), 2) == 1
@test size(Data.schema(f), 1) == 3
@test Data.types(Data.schema(f)) == (Date,)
f = CSV.Source(joinpath(dir, "test_datetimes.csv"); dateformat="yyyy-mm-dd HH:MM:SS.s", allowmissing=:auto)
@test size(Data.schema(f), 2) == 1
@test size(Data.schema(f), 1) == 3
@test Data.types(Data.schema(f)) == (DateTime,)
ds = CSV.read(f)
@test ds[1][1] == DateTime(2015,1,1)
@test ds[1][2] == DateTime(2015,1,2,0,0,1)
@test ds[1][3] == DateTime(2015,1,3,0,12,0,1)

#test bad types
f = CSV.Source(joinpath(dir, "test_float_in_int_column.csv"); types=[Int,Int,Int])
@test_throws CSV.ParsingException CSV.read(f)

#test missing values
f = CSV.Source(joinpath(dir, "test_missing_value_NULL.csv"); categorical=false, allowmissing=:auto)
@test size(Data.schema(f), 2) == 3
@test size(Data.schema(f), 1) == 3
@test Data.types(Data.schema(f)) == (Float64,String,Float64)
ds = CSV.read(f)
@test ds[1][1] == 1.0
@test string(ds[2][1]) == "2.0"
@test string(ds[2][2]) == "NULL"
f = CSV.Source(joinpath(dir, "test_missing_value_NULL.csv"); missingstring="NULL", allowmissing=:auto)
@test size(Data.schema(f), 2) == 3
@test size(Data.schema(f), 1) == 3
@test f.options.missingstring == codeunits("NULL")
@test Data.types(Data.schema(f)) == (Float64,Union{Float64, Missing},Float64)
ds = CSV.read(f)
@test ds[1][1] == 1.0
@test ds[2][1] == 2.0
@test ismissing(ds[2][2])

# uses default missing value ""
f = CSV.Source(joinpath(dir, "test_missing_value.csv"), allowmissing=:auto)
@test size(Data.schema(f), 2) == 3
@test size(Data.schema(f), 1) == 3
@test Data.types(Data.schema(f)) == (Float64,Union{Float64, Missing},Float64)

f = CSV.Source(joinpath(dir, "test_header_range.csv");header=1:3,allowmissing=:auto)
@test size(Data.schema(f), 2) == 3
@test size(Data.schema(f), 1) == 3
@test Data.header(Data.schema(f)) == ["col1_sub1_part1","col2_sub2_part2","col3_sub3_part3"]
ds = CSV.read(f)

f = CSV.Source(joinpath(dir, "test_header_range.csv");header=["col1_sub1_part1","col2_sub2_part2","col3_sub3_part3"],datarow=4,allowmissing=:auto)
@test size(Data.schema(f), 2) == 3
@test size(Data.schema(f), 1) == 3
@test Data.header(Data.schema(f)) == ["col1_sub1_part1","col2_sub2_part2","col3_sub3_part3"]
ds = CSV.read(f)

f = CSV.Source(joinpath(dir, "test_basic.csv");types=Dict(2=>Float64),allowmissing=:auto)
@test size(Data.schema(f), 2) == 3
@test size(Data.schema(f), 1) == 3
@test Data.types(Data.schema(f)) == (Int64,Float64,Int64)
ds = CSV.read(f)

f = CSV.Source(joinpath(dir, "test_basic_pipe.csv");delim='|',allowmissing=:auto)
@test size(Data.schema(f), 2) == 3
@test size(Data.schema(f), 1) == 3
@test Data.types(Data.schema(f)) == (Int64,Int64,Int64)
@test f.options.delim == UInt8('|')
ds = CSV.read(f)

f = CSV.Source(joinpath(dir, "test_basic_pipe.csv");delim='|',footerskip=1,allowmissing=:auto)
@test size(Data.schema(f), 2) == 3
@test size(Data.schema(f), 1) == 2
@test Data.types(Data.schema(f)) == (Int64,Int64,Int64)
@test f.options.delim == UInt8('|')
ds = CSV.read(f)

t = tempname()
f = CSV.Sink(t)

f = CSV.Source(joinpath(dir, "test_missing_value_NULL.csv"))
types = Data.types(Data.schema(f))
@test CSV.parsefield(f, types[1]) == 1.0
@test CSV.parsefield(f, types[2]) == "2.0"
@test CSV.parsefield(f, types[3]) == 3.0

t = tempname()
f = open(t, "w+")
Base.write(f, String(read(joinpath(dir, "test_missing_value_NULL.csv"))))
seekstart(f)
source = CSV.Source(f; header=[], datarow=2, allowmissing=:none)
df = CSV.read(source)
@test Data.header(Data.schema(df)) == ["Column1", "Column2", "Column3"]

Data.reset!(source)
df2 = CSV.read(source)
@test df == df2

@test_throws ArgumentError CSV.Source(f; types = [Int, Int, Int, Int])
close(f)
f = source = nothing; GC.gc(); GC.gc()
try
rm(t)
catch e
end

# test tab-delimited missing values
d = CSV.read(joinpath(dir, "test_tab_null_empty.txt"); delim='\t')
@test ismissing(d[2][2])

d = CSV.read(joinpath(dir, "test_tab_null_string.txt"); delim='\t', missingstring="NULL")
@test ismissing(d[2][2])

# read a write protected file
let fn = tempname()
    open(fn, "w") do f
        write(f, "Julia")
    end
    chmod(fn, 0o444)
    CSV.read(fn)
    GC.gc(); GC.gc()
    try
    rm(fn)
    catch e
    end
end

# CSV with header and no data is treated the same as an empty buffer with header supplied
df1 = CSV.read(IOBuffer("a,b,c"))
df2 = CSV.read(IOBuffer(""); header=["a", "b", "c"])
@test size(Data.schema(df1)) == (0, 3)
@test size(Data.schema(df2)) == (0, 3)
@test df1 == df2

# Adding transforms to CSV with header but no data returns empty frame as expected
# (previously the lack of a ::String dispatch in the transform function caused an error)
transforms = Dict{Int, Function}(2 => x::Integer -> "b$x")
df1 = CSV.read(IOBuffer("a,b,c\n1,2,3\n4,5,6"); allowmissing=:none, transforms=transforms)
df2 = CSV.read(IOBuffer("a,b,c\n1,b2,3\n4,b5,6"); allowmissing=:none)
@test size(Data.schema(df1)) == (2, 3)
@test size(Data.schema(df2)) == (2, 3)
@test df1 == df2
df3 = CSV.read(IOBuffer("a,b,c"); allowmissing=:none, transforms=transforms)
df4 = CSV.read(IOBuffer("a,b,c"); allowmissing=:none)
@test size(Data.schema(df3)) == (0, 3)
@test size(Data.schema(df4)) == (0, 3)
@test df3 == df4

let fn = tempname()
    df = CSV.read(IOBuffer("a,b,c\n1,2,3\n4,5,6"), CSV.Sink(fn); allowmissing=:none, transforms=transforms)
    @test String(read(fn)) == "a,b,c\n1,b2,3\n4,b5,6\n"
    try
    rm(fn)
    catch e
    end
end

let fn = tempname()
    df = CSV.read(IOBuffer("a,b,c"), CSV.Sink(fn); allowmissing=:none, transforms=transforms)
    @test String(read(fn)) == "a,b,c\n"
    try
    rm(fn)
    catch e
    end
end

source = IOBuffer("col1,col2,col3") # empty dataset
df = CSV.read(source; transforms=Dict(2 => floor))
@test size(Data.schema(df)) == (0, 3)
@test Data.types(Data.schema(df)) == (Any, Any, Any)

# Integer overflow; #100
@test_throws OverflowError CSV.read(joinpath(dir, "int8_overflow.csv"); types=[Int8])

# dash as missingstring; #92
df = CSV.read(joinpath(dir, "dash_as_null.csv"); missingstring="-")
@test ismissing(df[1][2])

df = CSV.read(joinpath(dir, "plus_as_null.csv"); missingstring="+")
@test ismissing(df[1][2])

# #83
df = CSV.read(joinpath(dir, "comma_decimal.csv"); delim=';', decimal=',')
@test df[1][1] === 3.14
@test df[1][2] === 1.0
@test df[2][1] === Int64(1)
@test df[2][2] === Int64(1)

# #86
df = CSV.read(joinpath(dir, "double_quote_quotechar_and_escapechar.csv"); escapechar='"')
@test size(df) == (24, 5)
@test df[5][24] == "NORTH DAKOTA STATE \"\"A\"\" #1"

# #84
df = CSV.read(joinpath(dir, "census.txt"); delim='\t', allowmissing=:auto)
@test eltype(df[9]) == Float64
@test size(df) == (3, 9)

# #79
df = CSV.read(joinpath(dir, "bools.csv"), allowmissing=:auto)
@test eltype(df[1]) == Bool
@test df[1] == [true, false, true, false]
@test df[2] == [false, true, true, false]
@test df[3] == [1, 2, 3, 4]

# #64
df = CSV.read(joinpath(dir, "attenu.csv"), missingstring="NA", types=Dict(3=>Union{Missing, String}))
@test size(df) == (182, 5)

f = CSV.Source(joinpath(dir, "test_null_only_column.csv"), categorical=false, missingstring="NA", allowmissing=:auto)
@test size(Data.schema(f)) == (3, 2)
ds = CSV.read(f)
@test Data.types(Data.schema(f)) == (String, Missing)
@test all(ismissing, ds[2])

# #107
df = CSV.read(IOBuffer("1,a,i\n2,b,ii\n3,c,iii"); datarow=1)
@test size(df) == (3, 3)

# #115 (Int64 -> Union{Int64, Missing} -> Union{String, Missing} promotion)
df = CSV.read(joinpath(dir, "attenu.csv"), missingstring="NA", rows_for_type_detect=200, allowmissing=:auto)
@test size(df) == (182, 5)
@test Data.types(Data.schema(df)) == (Int64, Float64, Union{Missings.Missing, String}, Float64, Float64)

# #137
tbl = DataFrame(a=[11,22], dt=[Date(2017,12,7), Date(2017,12,14)])
tbl[:dttm] = DateTime.(tbl[:dt])
CSV.write("test.tsv", tbl; delim='\t')
df = CSV.read("test.tsv"; delim='\t', allowmissing=:auto)
@test Data.types(Data.schema(df)) == (Int64, Date, DateTime)
df = nothing; GC.gc(); GC.gc()
try
rm("test.tsv")
catch e
end

end # testset

@testset "CSV.Source various files" begin

#other various files found around the internet
f = CSV.Source(joinpath(dir, "baseball.csv"); rows_for_type_detect=35, strings=:weakref)
@test size(Data.schema(f), 2) == 15
@test size(Data.schema(f), 1) == 35
@test Data.header(Data.schema(f)) == ["Rk","Year","Age","Tm","Lg","","W","L","W-L%","G","Finish","Wpost","Lpost","W-L%post",""]
@test Data.types(Data.schema(f)) == (Union{Int64, Missing},Union{Int64, Missing},Union{Int64, Missing},Union{CategoricalString{UInt32}, Missing},Union{CategoricalString{UInt32}, Missing},Union{WeakRefString{UInt8}, Missing},Union{Int64, Missing},Union{Int64, Missing},Union{Float64, Missing},Union{Int64, Missing},Union{Float64, Missing},Union{Int64, Missing},Union{Int64, Missing},Union{Float64, Missing},Union{CategoricalString{UInt32}, Missing})
ds = CSV.read(f)

f = CSV.Source(joinpath(dir, "FL_insurance_sample.csv"); types=Dict(10=>Float64,12=>Float64), allowmissing=:auto)
@test size(Data.schema(f), 2) == 18
@test size(Data.schema(f), 1) == 36634
@test Data.header(Data.schema(f)) == ["policyID","statecode","county","eq_site_limit","hu_site_limit","fl_site_limit","fr_site_limit","tiv_2011","tiv_2012","eq_site_deductible","hu_site_deductible","fl_site_deductible","fr_site_deductible","point_latitude","point_longitude","line","construction","point_granularity"]
@test Data.types(Data.schema(f)) == (Int64,CategoricalString{UInt32},CategoricalString{UInt32},Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Int64,Float64,Float64,CategoricalString{UInt32},CategoricalString{UInt32},Int64)
ds = CSV.read(f)

f = CSV.Source(joinpath(dir, "FL_insurance_sample.csv");types=Dict("eq_site_deductible"=>Float64,"fl_site_deductible"=>Float64),allowmissing=:auto)
@test size(Data.schema(f), 2) == 18
@test size(Data.schema(f), 1) == 36634
@test Data.header(Data.schema(f)) == ["policyID","statecode","county","eq_site_limit","hu_site_limit","fl_site_limit","fr_site_limit","tiv_2011","tiv_2012","eq_site_deductible","hu_site_deductible","fl_site_deductible","fr_site_deductible","point_latitude","point_longitude","line","construction","point_granularity"]
@test Data.types(Data.schema(f)) == (Int64,CategoricalString{UInt32},CategoricalString{UInt32},Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Int64,Float64,Float64,CategoricalString{UInt32},CategoricalString{UInt32},Int64)
ds = CSV.read(f)

f = CSV.Source(joinpath(dir, "SacramentocrimeJanuary2006.csv"), allowmissing=:auto)
@test size(Data.schema(f), 2) == 9
@test size(Data.schema(f), 1) == 7584
@test Data.header(Data.schema(f)) == ["cdatetime","address","district","beat","grid","crimedescr","ucr_ncic_code","latitude","longitude"]
@test Data.types(Data.schema(f)) == (CategoricalString{UInt32},String,Int64,CategoricalString{UInt32},Int64,CategoricalString{UInt32},Int64,Float64,Float64)
ds = CSV.read(f)

f = CSV.Source(joinpath(dir, "Sacramentorealestatetransactions.csv"), allowmissing=:auto)
@test size(Data.schema(f), 2) == 12
@test size(Data.schema(f), 1) == 985
@test Data.header(Data.schema(f)) == ["street","city","zip","state","beds","baths","sq__ft","type","sale_date","price","latitude","longitude"]
@test Data.types(Data.schema(f)) == (String,CategoricalString{UInt32},Int64,CategoricalString{UInt32},Int64,Int64,Int64,CategoricalString{UInt32},CategoricalString{UInt32},Int64,Float64,Float64)
ds = CSV.read(f)

f = CSV.Source(joinpath(dir, "SalesJan2009.csv"); types=Dict(3=>WeakRefString{UInt8},7=>Union{WeakRefString{UInt8}, Missing}), strings=:weakref, allowmissing=:auto)
@test size(Data.schema(f), 2) == 12
@test size(Data.schema(f), 1) == 998
@test Data.header(Data.schema(f)) == ["Transaction_date","Product","Price","Payment_Type","Name","City","State","Country","Account_Created","Last_Login","Latitude","Longitude"]
@test Data.types(Data.schema(f)) == (WeakRefString{UInt8},CategoricalString{UInt32},WeakRefString{UInt8},CategoricalString{UInt32},WeakRefString{UInt8},WeakRefString{UInt8},Union{WeakRefString{UInt8}, Missing},CategoricalString{UInt32},WeakRefString{UInt8},WeakRefString{UInt8},Float64,Float64)
ds = CSV.read(f)

f = CSV.Source(joinpath(dir, "stocks.csv"), allowmissing=:auto)
@test size(Data.schema(f), 2) == 2
@test size(Data.schema(f), 1) == 30
@test Data.header(Data.schema(f)) == ["Stock Name","Company Name"]
@test Data.types(Data.schema(f)) == (String,String)
ds = CSV.read(f)

f = CSV.Source(joinpath(dir, "TechCrunchcontinentalUSA.csv"); types=Dict(4=>Union{String, Missing},5=>Union{String, Missing}), allowmissing=:auto)
@test size(Data.schema(f), 2) == 10
@test size(Data.schema(f), 1) == 1460
@test Data.header(Data.schema(f)) == ["permalink","company","numEmps","category","city","state","fundedDate","raisedAmt","raisedCurrency","round"]
@test Data.types(Data.schema(f)) == (CategoricalString{UInt32},CategoricalString{UInt32},Union{Int64, Missing},Union{String, Missing},Union{String, Missing},CategoricalString{UInt32},CategoricalString{UInt32},Int64,CategoricalString{UInt32},CategoricalString{UInt32})
ds = CSV.read(f)

f = CSV.Source(joinpath(dir, "Fielding.csv"); allowmissing=:none, types=Dict("GS"=>Union{Int, Missing},"PO"=>Union{Int, Missing},"A"=>Union{Int, Missing},"E"=>Union{Int, Missing},"DP"=>Union{Int, Missing},"PB"=>Union{Int, Missing},"InnOuts"=>Union{Int, Missing},"WP"=>Union{Int, Missing},"SB"=>Union{Int, Missing},"CS"=>Union{Int, Missing},"ZR"=>Union{Int, Missing}))
@test size(Data.schema(f), 2) == 18
@test size(Data.schema(f), 1) == 167938
@test Data.header(Data.schema(f)) == ["playerID","yearID","stint","teamID","lgID","POS","G","GS","InnOuts","PO","A","E","DP","PB","WP","SB","CS","ZR"]
@test Data.types(Data.schema(f)) == (CategoricalArrays.CategoricalString{UInt32}, Int64, Int64, CategoricalArrays.CategoricalString{UInt32}, CategoricalArrays.CategoricalString{UInt32}, CategoricalArrays.CategoricalString{UInt32}, Int64, Union{Int, Missings.Missing}, Union{Int, Missings.Missing}, Union{Int, Missings.Missing}, Union{Int, Missings.Missing}, Union{Int, Missings.Missing}, Union{Int, Missings.Missing}, Union{Int, Missings.Missing}, Union{Int, Missings.Missing}, Union{Int, Missings.Missing}, Union{Int, Missings.Missing}, Union{Int, Missings.Missing})
ds = CSV.read(f)

f = CSV.Source(joinpath(dir, "latest (1).csv"); header=0, missingstring="\\N", types=Dict(13=>Union{Float64, Missing},17=>Union{Int, Missing},18=>Union{Float64, Missing},20=>Union{Float64, Missing}), allowmissing=:auto)
@test size(Data.schema(f), 2) == 25
@test size(Data.schema(f), 1) == 1000
@test Data.header(Data.schema(f)) == ["Column$i" for i = 1:size(Data.schema(f), 2)]
@test Data.types(Data.schema(f)) == (CategoricalArrays.CategoricalString{UInt32}, CategoricalArrays.CategoricalString{UInt32}, Int64, Int64, CategoricalArrays.CategoricalString{UInt32}, Int64, CategoricalArrays.CategoricalString{UInt32}, Int64, Date, Date, Int64, CategoricalArrays.CategoricalString{UInt32}, Union{Float64, Missings.Missing}, Union{Float64, Missings.Missing}, Union{Float64, Missings.Missing}, Union{Float64, Missings.Missing}, Union{Int, Missings.Missing}, Union{Float64, Missings.Missing}, Float64, Union{Float64, Missings.Missing}, Union{Float64, Missings.Missing}, Union{Int64, Missings.Missing}, Float64, Union{Float64, Missings.Missing}, Union{Float64, Missings.Missing})
ds = CSV.read(f)

f = CSV.Source(joinpath(dir, "pandas_zeros.csv"), allowmissing=:none)
@test size(Data.schema(f), 2) == 50
@test size(Data.schema(f), 1) == 100000
@test Data.header(Data.schema(f)) == [string(i) for i = 0:49]
@test Data.types(Data.schema(f)) == (fill(Int64,50)...,)
@time ds = CSV.read(f)

end # testset

@testset "CSV.TransposedSource" begin

# CSV.TransposedSource
df = CSV.read(joinpath(dir, "transposed.csv"); transpose=true)
@test size(df) == (3, 3)
@test Data.header(Data.schema(df)) == ["col1", "col2", "col3"]
@test df[1][1] == 1
@test df[1][2] == 2
@test df[1][3] == 3

df = CSV.read(joinpath(dir, "transposed_1row.csv"); transpose=true)
@test size(df) == (1, 1)

df = CSV.read(joinpath(dir, "transposed_emtpy.csv"); transpose=true)
@test size(df) == (0, 1)

df = CSV.read(joinpath(dir, "transposed_extra_newline.csv"); transpose=true)
@test size(df) == (2, 2)

df = CSV.read(joinpath(dir, "transposed_noheader.csv"); transpose=true, header=0)
@test size(df) == (2, 3)
@test Data.header(Data.schema(df)) == ["Column1", "Column2", "Column3"]

df = CSV.read(joinpath(dir, "transposed_noheader.csv"); transpose=true, header=["c1", "c2", "c3"])
@test size(df) == (2, 3)
@test Data.header(Data.schema(df)) == ["c1", "c2", "c3"]

end # testset


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

 ## Basic functions ##

isreal(x::AbstractArray) = all(isreal,x)
iszero(x::AbstractArray) = all(iszero,x)
isreal(x::AbstractArray{<:Real}) = true

## Constructors ##

"""
    vec(a::AbstractArray) -> AbstractVector

Reshape the array `a` as a one-dimensional column vector. Return `a` if it is
already an `AbstractVector`. The resulting array
shares the same underlying data as `a`, so it will only be mutable if `a` is
mutable, in which case modifying one will also modify the other.

# Examples
```jldoctest
julia> a = [1 2 3; 4 5 6]
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

julia> vec(a)
6-element Vector{Int64}:
 1
 4
 2
 5
 3
 6

julia> vec(1:3)
1:3
```

See also [`reshape`](@ref), [`dropdims`](@ref).
"""
vec(a::AbstractArray) = reshape(a,length(a))
vec(a::AbstractVector) = a

_sub(::Tuple{}, ::Tuple{}) = ()
_sub(t::Tuple, ::Tuple{}) = t
_sub(t::Tuple, s::Tuple) = _sub(tail(t), tail(s))

"""
    dropdims(A; dims)

Return an array with the same data as `A`, but with the dimensions specified by
`dims` removed. `size(A,d)` must equal 1 for every `d` in `dims`,
and repeated dimensions or numbers outside `1:ndims(A)` are forbidden.

The result shares the same underlying data as `A`, such that the
result is mutable if and only if `A` is mutable, and setting elements of one
alters the values of the other.

See also: [`reshape`](@ref), [`vec`](@ref).

# Examples
```jldoctest
julia> a = reshape(Vector(1:4),(2,2,1,1))
2×2×1×1 Array{Int64, 4}:
[:, :, 1, 1] =
 1  3
 2  4

julia> b = dropdims(a; dims=3)
2×2×1 Array{Int64, 3}:
[:, :, 1] =
 1  3
 2  4

julia> b[1,1,1] = 5; a
2×2×1×1 Array{Int64, 4}:
[:, :, 1, 1] =
 5  3
 2  4
```
"""
dropdims(A; dims) = _dropdims(A, dims)
function _dropdims(A::AbstractArray, dims::Dims)
    for i in eachindex(dims)
        1 <= dims[i] <= ndims(A) || throw(ArgumentError("dropped dims must be in range 1:ndims(A)"))
        length(axes(A, dims[i])) == 1 || throw(ArgumentError("dropped dims must all be size 1"))
        for j = 1:i-1
            dims[j] == dims[i] && throw(ArgumentError("dropped dims must be unique"))
        end
    end
    d = ()
    for i = 1:ndims(A)
        if !in(i, dims)
            d = tuple(d..., axes(A, i))
        end
    end
    reshape(A, d::typeof(_sub(axes(A), dims)))
end
_dropdims(A::AbstractArray, dim::Integer) = _dropdims(A, (Int(dim),))

## Unary operators ##

conj(x::AbstractArray{<:Real}) = x
conj!(x::AbstractArray{<:Real}) = x

real(x::AbstractArray{<:Real}) = x
imag(x::AbstractArray{<:Real}) = zero(x)

+(x::AbstractArray{<:Number}) = x
*(x::AbstractArray{<:Number,2}) = x

# index A[:,:,...,i,:,:,...] where "i" is in dimension "d"

"""
    selectdim(A, d::Integer, i)

Return a view of all the data of `A` where the index for dimension `d` equals `i`.

Equivalent to `view(A,:,:,...,i,:,:,...)` where `i` is in position `d`.

See also: [`eachslice`](@ref).

# Examples
```jldoctest
julia> A = [1 2 3 4; 5 6 7 8]
2×4 Matrix{Int64}:
 1  2  3  4
 5  6  7  8

julia> selectdim(A, 2, 3)
2-element view(::Matrix{Int64}, :, 3) with eltype Int64:
 3
 7

julia> selectdim(A, 2, 3:4)
2×2 view(::Matrix{Int64}, :, 3:4) with eltype Int64:
 3  4
 7  8
```
"""
@inline selectdim(A::AbstractArray, d::Integer, i) = _selectdim(A, d, i, _setindex(i, d, map(Slice, axes(A))...))
@noinline function _selectdim(A, d, i, idxs)
    d >= 1 || throw(ArgumentError("dimension must be ≥ 1, got $d"))
    nd = ndims(A)
    d > nd && (i == 1 || throw(BoundsError(A, (ntuple(Returns(Colon()),d-1)..., i))))
    return view(A, idxs...)
end

function circshift(a::AbstractArray, shiftamt::Real)
    circshift!(similar(a), a, (Integer(shiftamt),))
end
circshift(a::AbstractArray, shiftamt::DimsInteger) = circshift!(similar(a), a, shiftamt)
"""
    circshift(A, shifts)

Circularly shift, i.e. rotate, the data in an array. The second argument is a tuple or
vector giving the amount to shift in each dimension, or an integer to shift only in the
first dimension.

See also: [`circshift!`](@ref), [`circcopy!`](@ref), [`bitrotate`](@ref), [`<<`](@ref).

# Examples
```jldoctest
julia> b = reshape(Vector(1:16), (4,4))
4×4 Matrix{Int64}:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

julia> circshift(b, (0,2))
4×4 Matrix{Int64}:
  9  13  1  5
 10  14  2  6
 11  15  3  7
 12  16  4  8

julia> circshift(b, (-1,0))
4×4 Matrix{Int64}:
 2  6  10  14
 3  7  11  15
 4  8  12  16
 1  5   9  13

julia> a = BitArray([true, true, false, false, true])
5-element BitVector:
 1
 1
 0
 0
 1

julia> circshift(a, 1)
5-element BitVector:
 1
 1
 1
 0
 0

julia> circshift(a, -1)
5-element BitVector:
 1
 0
 0
 1
 1
```
"""
function circshift(a::AbstractArray, shiftamt)
    circshift!(similar(a), a, map(Integer, (shiftamt...,)))
end

## Other array functions ##

"""
    repeat(A::AbstractArray, counts::Integer...)

Construct an array by repeating array `A` a given number of times in each dimension, specified by `counts`.

See also: [`fill`](@ref), [`Iterators.repeated`](@ref), [`Iterators.cycle`](@ref).

# Examples
```jldoctest
julia> repeat([1, 2, 3], 2)
6-element Vector{Int64}:
 1
 2
 3
 1
 2
 3

julia> repeat([1, 2, 3], 2, 3)
6×3 Matrix{Int64}:
 1  1  1
 2  2  2
 3  3  3
 1  1  1
 2  2  2
 3  3  3
```
"""
function repeat(A::AbstractArray, counts...)
    return _RepeatInnerOuter.repeat(A, outer=counts)
end

"""
    repeat(A::AbstractArray; inner=ntuple(Returns(1), ndims(A)), outer=ntuple(Returns(1), ndims(A)))

Construct an array by repeating the entries of `A`. The i-th element of `inner` specifies
the number of times that the individual entries of the i-th dimension of `A` should be
repeated. The i-th element of `outer` specifies the number of times that a slice along the
i-th dimension of `A` should be repeated. If `inner` or `outer` are omitted, no repetition
is performed.

# Examples
```jldoctest
julia> repeat(1:2, inner=2)
4-element Vector{Int64}:
 1
 1
 2
 2

julia> repeat(1:2, outer=2)
4-element Vector{Int64}:
 1
 2
 1
 2

julia> repeat([1 2; 3 4], inner=(2, 1), outer=(1, 3))
4×6 Matrix{Int64}:
 1  2  1  2  1  2
 1  2  1  2  1  2
 3  4  3  4  3  4
 3  4  3  4  3  4
```
"""
function repeat(A::AbstractArray; inner = nothing, outer = nothing)
    return _RepeatInnerOuter.repeat(A, inner=inner, outer=outer)
end

module _RepeatInnerOuter

function repeat(arr; inner=nothing, outer=nothing)
    check(arr, inner, outer)
    arr, inner, outer = resolve(arr, inner, outer)
    repeat_inner_outer(arr, inner, outer)
end

to_tuple(t::Tuple) = t
to_tuple(x::Integer) = (x,)
to_tuple(itr) = tuple(itr...)

function pad(a, b)
    N = max(length(a), length(b))
    Base.fill_to_length(a, 1, Val(N)), Base.fill_to_length(b, 1, Val(N))
end
function pad(a, b, c)
    N = max(max(length(a), length(b)), length(c))
    Base.fill_to_length(a, 1, Val(N)), Base.fill_to_length(b, 1, Val(N)), Base.fill_to_length(c, 1, Val(N))
end

function resolve(arr::AbstractArray{<:Any, N}, inner::NTuple{N, Any}, outer::NTuple{N,Any}) where {N}
    arr, inner, outer
end
function resolve(arr, inner, outer)
    dims, inner, outer = pad(size(arr), to_tuple(inner), to_tuple(outer))
    reshape(arr, dims), inner, outer
end
function resolve(arr, inner::Nothing, outer::Nothing)
    return arr, inner, outer
end
function resolve(arr, inner::Nothing, outer)
    dims, outer = pad(size(arr), to_tuple(outer))
    reshape(arr, dims), inner, outer
end
function resolve(arr, inner, outer::Nothing)
    dims, inner = pad(size(arr), to_tuple(inner))
    reshape(arr, dims), inner, outer
end

function check(arr, inner, outer)
    if inner !== nothing
        # TODO: Currently one based indexing is demanded for inner !== nothing,
        # but not for outer !== nothing. Decide for something consistent.
        Base.require_one_based_indexing(arr)
        if any(<(0), inner)
            throw(ArgumentError("no inner repetition count may be negative; got $inner"))
        end
        if length(inner) < ndims(arr)
            throw(ArgumentError("number of inner repetitions ($(length(inner))) cannot be less than number of dimensions of input array ($(ndims(arr)))"))
        end
    end
    if outer !== nothing
        if any(<(0), outer)
            throw(ArgumentError("no outer repetition count may be negative; got $outer"))
        end
        if (length(outer) < ndims(arr)) && (inner !== nothing)
            throw(ArgumentError("number of outer repetitions ($(length(outer))) cannot be less than number of dimensions of input array ($(ndims(arr)))"))
        end
    end
end

repeat_inner_outer(arr, inner::Nothing, outer::Nothing) = arr
repeat_inner_outer(arr, ::Nothing, outer) = repeat_outer(arr, outer)
repeat_inner_outer(arr, inner, ::Nothing) = repeat_inner(arr, inner)
repeat_inner_outer(arr, inner, outer) = repeat_outer(repeat_inner(arr, inner), outer)

function repeat_outer(a::AbstractMatrix, (m,n)::NTuple{2, Any})
    o, p = size(a,1), size(a,2)
    b = similar(a, o*m, p*n)
    for j=1:n
        d = (j-1)*p+1
        R = d:d+p-1
        for i=1:m
            c = (i-1)*o+1
            @inbounds b[c:c+o-1, R] = a
        end
    end
    return b
end

function repeat_outer(a::AbstractVector, (m,)::Tuple{Any})
    o = length(a)
    b = similar(a, o*m)
    for i=1:m
        c = (i-1)*o+1
        @inbounds b[c:c+o-1] = a
    end
    return b
end

function repeat_outer(arr::AbstractArray{<:Any,N}, dims::NTuple{N,Any}) where {N}
    insize  = size(arr)
    outsize = map(*, insize, dims)
    out = similar(arr, outsize)
    for I in CartesianIndices(arr)
        for J in CartesianIndices(dims)
            TIJ = map(Tuple(I), Tuple(J), insize) do i, j, d
                i + d * (j-1)
            end
            IJ = CartesianIndex(TIJ)
            @inbounds out[IJ] = arr[I]
        end
    end
    return out
end

function repeat_inner(arr, inner)
    basedims = size(arr)
    outsize = map(*, size(arr), inner)
    out = similar(arr, outsize)
    for I in CartesianIndices(arr)
        for J in CartesianIndices(inner)
            TIJ = map(Tuple(I), Tuple(J), inner) do i, j, d
                (i-1) * d + j
            end
            IJ = CartesianIndex(TIJ)
            @inbounds out[IJ] = arr[I]
        end
    end
    return out
end

end#module

"""
    eachrow(A::AbstractVecOrMat)

Create a generator that iterates over the first dimension of vector or matrix `A`,
returning the rows as `AbstractVector` views.

See also [`eachcol`](@ref), [`eachslice`](@ref), [`mapslices`](@ref).

!!! compat "Julia 1.1"
     This function requires at least Julia 1.1.

# Example

```jldoctest
julia> a = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> first(eachrow(a))
2-element view(::Matrix{Int64}, 1, :) with eltype Int64:
 1
 2

julia> collect(eachrow(a))
2-element Vector{SubArray{Int64, 1, Matrix{Int64}, Tuple{Int64, Base.Slice{Base.OneTo{Int64}}}, true}}:
 [1, 2]
 [3, 4]
```
"""
eachrow(A::AbstractVecOrMat) = (view(A, i, :) for i in axes(A, 1))

"""
    eachcol(A::AbstractVecOrMat)

Create a generator that iterates over the second dimension of matrix `A`, returning the
columns as `AbstractVector` views.

See also [`eachrow`](@ref) and [`eachslice`](@ref).

!!! compat "Julia 1.1"
     This function requires at least Julia 1.1.

# Example

```jldoctest
julia> a = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> first(eachcol(a))
2-element view(::Matrix{Int64}, :, 1) with eltype Int64:
 1
 3

julia> collect(eachcol(a))
2-element Vector{SubArray{Int64, 1, Matrix{Int64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}}:
 [1, 3]
 [2, 4]
```
"""
eachcol(A::AbstractVecOrMat) = (view(A, :, i) for i in axes(A, 2))

"""
    eachslice(A::AbstractArray; dims)

Create a generator that iterates over dimensions `dims` of `A`, returning views that select all
the data from the other dimensions in `A`.

Only a single dimension in `dims` is currently supported. Equivalent to `(view(A,:,:,...,i,:,:
...)) for i in axes(A, dims))`, where `i` is in position `dims`.

See also [`eachrow`](@ref), [`eachcol`](@ref), [`mapslices`](@ref), and [`selectdim`](@ref).

!!! compat "Julia 1.1"
     This function requires at least Julia 1.1.

# Example

```jldoctest
julia> M = [1 2 3; 4 5 6; 7 8 9]
3×3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9

julia> first(eachslice(M, dims=1))
3-element view(::Matrix{Int64}, 1, :) with eltype Int64:
 1
 2
 3

julia> collect(eachslice(M, dims=2))
3-element Vector{SubArray{Int64, 1, Matrix{Int64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}}:
 [1, 4, 7]
 [2, 5, 8]
 [3, 6, 9]
```
"""
@inline function eachslice(A::AbstractArray; dims)
    length(dims) == 1 || throw(ArgumentError("only single dimensions are supported"))
    dim = first(dims)
    dim <= ndims(A) || throw(DimensionMismatch("A doesn't have $dim dimensions"))
    inds_before = ntuple(Returns(:), dim-1)
    inds_after = ntuple(Returns(:), ndims(A)-dim)
    return (view(A, inds_before..., i, inds_after...) for i in axes(A, dim))
end

max_issues_repo_name: simsurace/julia

max_forks_repo_name: simsurace/julia


text:
# Compute petri net adjacency matrix

using SparseArrays

# The adjacency matrix that is returned has dimensions (n_species + n_reactions)^2
function petri_net(rn::ReactionNetwork)
	n_species = rn.n_species
	reactions = rn.reactions

	n_reactions = length(reactions)

	n = n_species + n_reactions

	# start with empty adjacency matrix
	adjm = spzeros(Int,n,n)

	for j = 1:length(reactions)
		r = reactions[j]
		nr = length(r.reactants)
		for i=1:nr
			if r.kf != 0
				adjm[r.reactants[i],n_species+j] = r.stoichr[i]
			end
			if r.kr != 0
				adjm[n_species+j,r.reactants[i]] = r.stoichr[i]
			end
		end
		nr = length(r.products)
		for i=1:nr
			if r.kf != 0
				adjm[n_species+j,r.products[i]] = r.stoichp[i]
			end
			if r.kr != 0
				adjm[r.products[i],n_species+j] = r.stoichp[i]
			end
		end

	end

	return adjm
end


text:
module HSARuntime

# copied from CUDAdrv/src/CUDAdrv.jl
const ext = joinpath(dirname(@__DIR__), "deps", "ext.jl")
isfile(ext) || error("HSARuntime.jl has not been built, please run Pkg.build(\"HSARuntime\").")
include(ext)
if !configured
    # default (non-functional) values for critical variables,
    # making it possible to _load_ the package at all times.
    const libhsaruntime_version = v"0.0"
    const libhsaruntime_vendor = "none"
    const libhsaruntime_path = nothing
end

atexit() do
    configured && HSA.shut_down()
end

### Imports ###

using CEnum
using Setfield
using Libdl

### Exports ###

export HSAAgent, HSAQueue, HSAExecutable, HSAKernelInstance, HSAArray, HSASignal
export get_agents, profile, get_first_isa, launch!
export get_default_agent, get_default_queue

### HSA Runtime Wrapper ###

include(joinpath(@__DIR__, "HSA", "HSA.jl"))
import .HSA: Agent, Queue, Executable, Status, Signal

include("extras.jl")
include("error.jl")
include("agent.jl")
include("queue.jl")
include("signal.jl")
include("executable.jl")
include("kernel.jl")
include("memory.jl")
include("array.jl")

function __init__()
    # We want to always be able to load the package, regardless of
    # whether HSA libraries are found
    configured || return

    # Make sure we load the library found by the last `] build`
    push!(Libdl.DL_LOAD_PATH, dirname(libhsaruntime_path))
    # TODO: Do the same (if possible) for the debug library

    HSA.init() |> check

    agents = get_agents(:gpu)
    if length(agents) > 0
        DEFAULT_AGENT[] = first(agents)
    end
end

end


text:
import LinearAlgebra: dot

"""
  shape2D(s::Real,t::Real,order::Int)

Computes shape functions of linear or quadratic order on a quadrilateral element

INPUT:
  s:      [Real]
              horizontal location in computation domain [-1,1]^2
  t:      [Real]
              vertical location in computation domain [-1,1]^2
  order:  [Int]
              Order of basis function (1 -> linear, 2-> quadratic)

OUTPUT:
  phi:    [Vector{Real}]
              Nodal basis functions evaluated at (s,t)
  dphids  [Vector{Real}]
              Derivative in s-direction of nodal basis function
                evaluated at (s,t)
  dphidt:  [Vector{Real}]
              Derivative in t-direction of nodal basis function
                evaluated at (s,t)
"""
function shape2D(s::T,t::T,order::Int) where T<:Real
  ein = one(s)
  tt = typeof(s)

  # linear element
  if order == 1
    phi    = Array{tt}(undef,4)
    dphids = Array{tt}(undef,4)
    dphidt = Array{tt}(undef,4)

    phi[1]    = 0.25*(ein-s)*(ein-t)
    phi[2]    = 0.25*(ein+s)*(ein-t)
    phi[3]    = 0.25*(ein+s)*(ein+t)
    phi[4]    = 0.25*(ein-s)*(ein+t)
    dphids[1] = -0.25*(ein-t)
    dphids[2] =  0.25*(ein-t)
    dphids[3] =  0.25*(ein+t)
    dphids[4] = -0.25*(ein+t)
    dphidt[1] = -0.25*(ein-s)
    dphidt[2] = -0.25*(ein+s)
    dphidt[3] =  0.25*(ein+s)
    dphidt[4] =  0.25*(ein-s)

  # quadratic element
  elseif order == 2
    phi    = Array{tt}(undef,9)
    dphids = Array{tt}(undef,9)
    dphidt = Array{tt}(undef,9)

    ells1  = 0.5*s*(s-ein); ellt1 = 0.5*t*(t-ein);
    ells2  = ein-(s*s);     ellt2 = ein-(t*t);
    ells3  = 0.5*s*(s+ein); ellt3 = 0.5*t*(t+ein);
    dells1 = s-0.5;         dellt1 = t-0.5;
    dells2 = -2*s;          dellt2 = -2*t;
    dells3 = s+0.5;         dellt3 = t+0.5;

    phi[1] = ells1*ellt1
    phi[2] = ells3*ellt1
    phi[3] = ells3*ellt3
    phi[4] = ells1*ellt3
    phi[5] = ells2*ellt1
    phi[6] = ells3*ellt2
    phi[7] = ells2*ellt3
    phi[8] = ells1*ellt2
    phi[9] = ells2*ellt2
    dphids[1] = dells1*ellt1
    dphids[2] = dells3*ellt1
    dphids[3] = dells3*ellt3
    dphids[4] = dells1*ellt3
    dphids[5] = dells2*ellt1
    dphids[6] = dells3*ellt2
    dphids[7] = dells2*ellt3
    dphids[8] = dells1*ellt2
    dphids[9] = dells2*ellt2
    dphidt[1] = ells1*dellt1
    dphidt[2] = ells3*dellt1
    dphidt[3] = ells3*dellt3
    dphidt[4] = ells1*dellt3
    dphidt[5] = ells2*dellt1
    dphidt[6] = ells3*dellt2
    dphidt[7] = ells2*dellt3
    dphidt[8] = ells1*dellt2
    dphidt[9] = ells2*dellt2
  end

  return phi,dphids,dphidt
end

"""
  shape2D!(phi::Vector{T},dphids::Vector{T},dphidt::Vector{T},
           s::T,t::T,order::Int) where T<:Real

In-place version of shape2D; see shape2D function info.
"""
function shape2D!(phi,dphids,dphidt,s::T,t::T,order::Int) where T<:Real
  ein = one(s)

  # linear element
  if order == 1
    phi[1]    = 0.25*(ein-s)*(ein-t)
    phi[2]    = 0.25*(ein+s)*(ein-t)
    phi[3]    = 0.25*(ein+s)*(ein+t)
    phi[4]    = 0.25*(ein-s)*(ein+t)
    dphids[1] = -0.25*(ein-t)
    dphids[2] =  0.25*(ein-t)
    dphids[3] =  0.25*(ein+t)
    dphids[4] = -0.25*(ein+t)
    dphidt[1] = -0.25*(ein-s)
    dphidt[2] = -0.25*(ein+s)
    dphidt[3] =  0.25*(ein+s)
    dphidt[4] =  0.25*(ein-s)

  # quadratic element
  elseif order == 2
    ells1  = 0.5*s*(s-ein);  ellt1 = 0.5*t*(t-ein);
    ells2  = ein-(s*s);      ellt2 = ein-(t*t);
    ells3  = 0.5*s*(s+ein);  ellt3 = 0.5*t*(t+ein);
    dells1 = s-0.5;          dellt1 = t-0.5;
    dells2 = -2*s;           dellt2 = -2*t;
    dells3 = s+0.5;          dellt3 = t+0.5;

    phi[1] = ells1*ellt1
    phi[2] = ells3*ellt1
    phi[3] = ells3*ellt3
    phi[4] = ells1*ellt3
    phi[5] = ells2*ellt1
    phi[6] = ells3*ellt2
    phi[7] = ells2*ellt3
    phi[8] = ells1*ellt2
    phi[9] = ells2*ellt2
    dphids[1] = dells1*ellt1
    dphids[2] = dells3*ellt1
    dphids[3] = dells3*ellt3
    dphids[4] = dells1*ellt3
    dphids[5] = dells2*ellt1
    dphids[6] = dells3*ellt2
    dphids[7] = dells2*ellt3
    dphids[8] = dells1*ellt2
    dphids[9] = dells2*ellt2
    dphidt[1] = ells1*dellt1
    dphidt[2] = ells3*dellt1
    dphidt[3] = ells3*dellt3
    dphidt[4] = ells1*dellt3
    dphidt[5] = ells2*dellt1
    dphidt[6] = ells3*dellt2
    dphidt[7] = ells2*dellt3
    dphidt[8] = ells1*dellt2
    dphidt[9] = ells2*dellt2
  end

  nothing
end

"""
  shapeEval(coef::Vector{Real},basis::Vector{Real})

Evaluates the sum given by the basis values and coefficients. Most useful for
evaluating functions interior to finite elements when (s,t) are known

INPUT:
  coef:    [Vector{Real}]
              Coefficients for some quantity evaluated at some point
  basis:   [Vector{Real}]
              Nodal Basis functions evaluated at some point

OUTPUT:
  -        [Real]
              Interpolated value
"""
function shapeEval(coef::Vector{T},basis::Vector{T}) where T<:Real
  return dot(coef,basis)
end

"""
  derivShape2D(s,t,xNodes,yNodes,order)

Evaluates derivatives of shape functions

INPUT:
  s:       [Real]
              Horizontal component of point in computational domain [-1,1]^2
  t:       [Real]
              Vertical component of point in computational domain [-1,1]^2
  xNodes:  [Vector{Real}]
              Horizontal component of element nodes
  yNodes:  [Vector{Real}]
              Vertical component of element nodes
  order:   [Int]
              Order of basis polynomials on element

OUTPUT:
  phi:     [Vector{Real}]
              Nodal basis function evaluated at node (s,t)
  dphidx:  [Vector{Real}]
              Derivative in x-direction of nodal basis function
                evaluated at (s,t)
  dphidy:  [Vector{Real}]
              Derivative in y-direction of nodal basis function
                evaluated at (s,t)
  jac:     [Real]
              Jacobian of transformation from computation to physical
                domain (or is it the other way around...?)
"""
function derivShape2D(s::T,t::T,xNodes::Vector{T},yNodes::Vector{T},order::Int) where T<:Real
  if order==1
    numNodes=4
  elseif order==2
    numNodes=9
  else
    error("invalid order. Must be 1 or 2")
  end

  # evaluate shape functions
  phi,dphids,dphidt = shape2D(s,t,order)

  dxds = shapeEval(xNodes,dphids)
  dxdt = shapeEval(xNodes,dphidt)
  dyds = shapeEval(yNodes,dphids)
  dydt = shapeEval(yNodes,dphidt)

  # evaluates jacobian
  jac = dxds*dydt - dxdt*dyds

  # check element jacobian
  if jac < 1e-9
    println("Bad element warning ...\n")

    if jac <= 0.0
      error("singular Jacobian ... Aborted ...")
    end
  end

  # evaluate derivatives
  dphidx = @.  dphids*dydt - dphidt*dyds
  dphidy = @. -dphids*dxdt + dphidt*dxds

  return phi,dphidx,dphidy,jac
end
"""
  derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s,t,xNodes,yNodes,order)

Evaluates derivatives of shape functions in-place
"""
function derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s,t,xNodes,yNodes,order::Int)
  (order==1 ? numNodes=4 : numNodes=9)

  dxds = shapeEval(xNodes,dphids)
  dxdt = shapeEval(xNodes,dphidt)
  dyds = shapeEval(yNodes,dphids)
  dydt = shapeEval(yNodes,dphidt)

  # evaluate derivatives
  for i=1:numNodes
    dphidx[i] =  dphids[i]*dydt - dphidt[i]*dyds
    dphidy[i] = -dphids[i]*dxdt + dphidt[i]*dxds
  end

  nothing
end

"""
  jacCalc(xNodes,yNodes,dphids,dphidt)

Computes jacobian of element
"""
function jacCalc(xNodes,yNodes,dphids,dphidt)
  dxds = shapeEval(xNodes,dphids)
  dxdt = shapeEval(xNodes,dphidt)
  dyds = shapeEval(yNodes,dphids)
  dydt = shapeEval(yNodes,dphidt)

  # evaluates jacobian
  jac = dxds*dydt - dxdt*dyds

  # check element jacobian
  if jac < 1e-9
    println("Bad element warning ...\n")

    if jac <= 0.0
      error("singular Jacobian ... Aborted ...")
    end
  end

  return jac
end

"""
  GaussQuadPoints2D(order::Int)

Produces Gauss Quadrature weights and points, normalized to (-1,1), for a given order

INPUT:
  order:  [Int]
             Order of gauss quadrature

OUTPUT:
  w:      [Vector{Real}]
             Weights of Guass Quadrature values
  s:      [Vector{Real}]
             Horizontal components of Gauss Quadrature nodes
  t:      [Vector{Real}]
             Horizontal components of Gauss Quadrature nodes
"""
function GaussQuadPoints2D(order::Int;tt=Float64::DataType)
  if order==2        # 2x2 Gauss points
    w = Array{tt,1}(undef,4)
    s = Array{tt,1}(undef,4)
    t = Array{tt,1}(undef,4)

    gpt=1.0e0/sqrt(3.0e0);
    s[1] = -gpt; t[1] = -gpt; w[1]=1.0;
    s[2] =  gpt; t[2] = -gpt; w[2]=1.0;
    s[3] =  gpt; t[3] =  gpt; w[3]=1.0;
    s[4] = -gpt; t[4] =  gpt; w[4]=1.0;
  elseif order==3    # 3x3 Gauss points
    w = zeros(tt,9)
    s = zeros(tt,9)
    t = zeros(tt,9)

    gpt=sqrt(0.6);
    s[1] = -gpt; t[1] = -gpt; w[1]=25/81
    s[2] =  gpt; t[2] = -gpt; w[2]=25/81
    s[3] =  gpt; t[3] =  gpt; w[3]=25/81
    s[4] = -gpt; t[4] =  gpt; w[4]=25/81
    s[5] =  0.0; t[5] = -gpt; w[5]=40/81
    s[6] =  gpt; t[6] =  0.0; w[6]=40/81
    s[7] =  0.0; t[7] =  gpt; w[7]=40/81
    s[8] = -gpt; t[8] =  0.0; w[8]=40/81
    s[9] =  0.0; t[9] =  0.0; w[9]=64/81
  elseif order==4      # 4x4 Gauss points
    w = zeros(tt,16)
    s = zeros(tt,16)
    t = zeros(tt,16)

    tmp  = 2/7*sqrt(6/5)
    gpt1 = sqrt(3/7-tmp)
    gpt2 = sqrt(3/7+tmp)
    wp   = (18+sqrt(30))/36
    wm   = (18-sqrt(30))/36
    s[1]  =  gpt1; t[1]  =  gpt1; w[1]  = wp*wp
    s[2]  =  gpt1; t[2]  = -gpt1; w[2]  = wp*wp
    s[3]  =  gpt1; t[3]  =  gpt2; w[3]  = wp*wm
    s[4]  =  gpt1; t[4]  = -gpt2; w[4]  = wp*wm
    s[5]  = -gpt1; t[5]  =  gpt1; w[5]  = wp*wp
    s[6]  = -gpt1; t[6]  = -gpt1; w[6]  = wp*wp
    s[7]  = -gpt1; t[7]  =  gpt2; w[7]  = wp*wm
    s[8]  = -gpt1; t[8]  = -gpt2; w[8]  = wp*wm

    s[9]  =  gpt2; t[9]  =  gpt1; w[9]  = wm*wp
    s[10] =  gpt2; t[10] = -gpt1; w[10] = wm*wp
    s[11] =  gpt2; t[11] =  gpt2; w[11] = wm*wm
    s[12] =  gpt2; t[12] = -gpt2; w[12] = wm*wm
    s[13] = -gpt2; t[13] =  gpt1; w[13] = wm*wp
    s[14] = -gpt2; t[14] = -gpt1; w[14] = wm*wp
    s[15] = -gpt2; t[15] =  gpt2; w[15] = wm*wm
    s[16] = -gpt2; t[16] = -gpt2; w[16] = wm*wm
  else
    error("Check Gauss point integration specification")
  end

  return w,s,t
end

"""
  GaussEdge(edgeNumber,order)

Returns weights and point values for edges of gauss points
edgeNumber = 1 -> between points 1 and 2
edgeNumber = 2 -> between points 2 and 3
edgeNumber = 3 -> between points 3 and 4
edgeNumber = 4 -> between points 4 and 1

INPUT:
  edgeNumber:  [Int]
                  Index to determine where edge of interest lies
  order:       [Int]
                  Order of basis interpolating function

OUTPUT:
  w:           [Vector{Real}]
                  Weights of Guass Quadrature values
  s:           [Vector{Real}]
                  Horizontal components of Gauss Quadrature nodes
  t:           [Vector{Real}]
                  Horizontal components of Gauss Quadrature nodes
"""
function GaussEdge(edgeNumber::Int,order::Int;tt=Float64::DataType)
  w,_ = GaussQuadPoints1D(order;tt=tt)

  if edgeNumber == 1
    _,s = GaussQuadPoints1D(order;tt=tt)
    t   = -ones(order+1)
  elseif edgeNumber==2
    s   = ones(order+1)
    _,t = GaussQuadPoints1D(order;tt=tt)
  elseif edgeNumber==3
    _,s = GaussQuadPoints1D(order;tt=tt)
    t   = ones(order+1)
  elseif edgeNumber==4
    s   = -ones(order+1)
    _,t = GaussQuadPoints1D(order;tt=tt)
  else
    error("Invalid Edge Number entered. Must be 1-4, but $(edgeNumber) was entered")
  end

  return w::Vector{tt},s::Vector{tt},t::Vector{tt}
end

"""
  FEMstats(xy,cm,dNodes,nNodes,F)

Prints to the terminal useful information about the system being solved

INPUT:
  mesh:   [AbstractMesh]
             Mesh object containing nodes and connectivity
  prob:   [AbstractProblem]
             Problem object containing boundary information

OUTPUT:
  none
"""
function FEMstats(mesh,prob)
  numTnodes = length(mesh.xy)
  Nel       = length(mesh.cm)
  numDnodes = length(prob.bcid[:dNodes])
  numNnodes = length(prob.bcid[:nNodes])

  println("# of total nodes:        ",numTnodes)
  println("# of elements:           ",Nel)
  println("# of Dirichlet BC nodes: ",numDnodes)
  println("# of Neumann BC nodes:   ",numNnodes)
  println("% boundary nodes:        ",
          round(numTnodes/(numDnodes+numNnodes),digits=2),"%")
end
function FEMstats(mesh)
  numTnodes = length(mesh.xy)
  Nel       = length(mesh.cm)

  println("# of total nodes:        ",numTnodes)
  println("# of elements:           ",Nel)
end

"""
	GaussQuadPoints1D(order)

Produces Gauss Quadrature weights and points, normalized to (-1,1), for a given order
"""
function GaussQuadPoints1D(order::Int;tt=Float64::DataType)
	w = zeros(tt,order)
	s = zeros(tt,order)

	if order == 2
		w[1] = one(tt)
		w[2] = one(tt)
		s[1] = sqrt(1/3)
		s[2] = -s[1]
  elseif order == 3
    w[1] = 8/9
    w[2] = 5/9
    w[3] = w[2]
    s[1] = 0.0
    s[2] =  sqrt(3/5)
    s[3] = -s[2]
  end

  return w, s
end


text:
export concat

import Base: show

"""
    concat(sources...)
    concat(sources::S) where { S <: Tuple }

Combines multiple Observables to create an Observable which sequentially emits all values from given Observable and then moves on to the next.
All values of each passed Observable merged into a single Observable, in order, in serial fashion.

# Arguments
- `sources`: input sources

# Examples
```jldoctest
using Rocket

source1 = of(1)
source2 = of(2)

subscribe!(concat(source1, source2), logger())
;

# output
[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Completed
```

```
using Rocket

source1 = of(1) |> async()
source2 = of(2)

subscribe!(concat(source1, source2), logger())
;

# output
[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Completed
```
See also: [`Subscribable`](@ref), [`subscribe!`](@ref)
"""
concat()                                = error("concat operator expects at least one inner observable on input")
concat(args...)                         = concat(args)
concat(sources::S) where { S <: Tuple } = ConcatObservable{union_type(sources), S}(sources)

##

@subscribable struct ConcatObservable{D, S} <: Subscribable{D}
    sources :: S
end

function on_subscribe!(observable::ConcatObservable{D, S}, actor::A) where { D, S, A }
    inner = ConcatInnerActor{D, S, A}(observable.sources, actor)
    subscription = subscribe!(observable.sources[1], inner)
    if get_current_index(inner) === 1
        set_subscription!(inner, subscription)
    end
    return ConcatSubscription(inner)
end

##

mutable struct ConcatInnerActor{D, S, A} <: Actor{D}
    sources       :: S
    actor         :: A
    current_index :: Int
    subscription  :: Teardown

    ConcatInnerActor{D, S, A}(sources::S, actor::A) where { D, S, A } = begin
        return new(sources, actor, 1, voidTeardown)
    end
end

get_current_index(actor::ConcatInnerActor)         = actor.current_index
set_current_index!(actor::ConcatInnerActor, index) = actor.current_index = index

get_subscription(actor::ConcatInnerActor)                = actor.subscription
set_subscription!(actor::ConcatInnerActor, subscription) = actor.subscription = subscription

function on_next!(actor::ConcatInnerActor{D}, data::D) where D
    next!(actor.actor, data)
end

function on_error!(actor::ConcatInnerActor, err)
    error!(actor.actor, err)
end

function on_complete!(actor::ConcatInnerActor)
    cindex = get_current_index(actor)
    if cindex === length(actor.sources)
        complete!(actor.actor)
    else
        set_current_index!(actor, cindex + 1)
        subscription = subscribe!(actor.sources[cindex + 1], actor)
        if get_current_index(actor) === cindex
            set_subscription!(actor, subscription)
        end
    end
end

##

struct ConcatSubscription{A} <: Teardown
    inner :: A
end

as_teardown(::Type{ <: ConcatSubscription}) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::ConcatSubscription)
    return unsubscribe!(get_subscription(subscription.inner))
end

Base.show(io::IO, ::ConcatObservable{D}) where D  = print(io, "ConcatObservable($D)")
Base.show(io::IO, ::ConcatSubscription)           = print(io, "ConcatSubscription()")


text:
"""
	unfold(ref, to_unf, samps=nothing; isomap=:default, optim=:default)

Unfold the input points `to_unf` based on the reference points informed in `ref`.
Returns a coordinate matrix with the unfolded points. Or a tuple of matrices if
`to_unf` is a list of points.

## Parameters:

* `to_unf`  - coordinate matrix with the points that will be unfolded
* `ref`     - coordinate matrix with the reference points for unfolding
* `isomap`  - additional parameters for isomap step; either :default or a
  NamedTuple with the keys shown below
* `optim`   - additional parameters for optimization step; either :default or a
  NamedTuple with the keys shown below

### `isomap` parameters:

*Default:* isomap = (search=:knn, neigh=16, anchors=1500, reftol=0.01)
* `search`  - search type to build neighbors graph for Isomap (:knn for
  k-nearest neighbor or :radius for radius search).
* `neigh`   - number of neighbors (for search=:knn) or radius distance
  (for search=:radius) to build neighbors graph for Isomap.
* `anchors` - number of anchors for landmark isomap
* `reftol`  - distance to which reference points are considered duplicates

### `optim` parameters:

*Default:* optim = (search=:knn, neigh=50, maxerr=5, nchunks=4)
* `search`  - search type to optimize locally (:knn for k-nearest neighbor or
  :radius for radius search).
* `neigh`   - number of neighbors (for search=:knn) or radius distance
  (for search=:radius) to local optimization.
* `maxerr`  - the maximum accepted distortion for unfolding optimization.
* `nchunks` - number of rounds of optimization.
"""
function unfold(to_unf::AbstractMatrix, ref::AbstractMatrix;
	            isomap=:default, optim=:default)

	# read isomap and optim parameters or assign the defaults below
	ipars = (search=:knn, neigh=16, anchors=1500, reftol=0.01)
	opars = (search=:knn, neigh=50, maxerr=5, nchunks=4)

	ipars = updatepars(ipars, isomap)
	opars = updatepars(opars, optim)

	# conversions if necessary
	!(ref[1] isa Float64) && (ref = Float64.(ref))
	!(to_unf[1] isa Float64) && (to_unf = Float64.(to_unf))

	# random seed
	Random.seed!(1234567890)

	# pre-process reference points
	ref = remove_duplicates(ref, tol=ipars.reftol)
	resol   = get_resolution(ref)

	# do landmark isomap at reference points
	unf_ref = landmark_isomap(ref, ipars.search, ipars.neigh, ipars.anchors)

	# get initial guess
	normals, good = getnormals(ref, ipars.search, ipars.neigh)
	unf = firstguess(to_unf, view(ref,:,good), view(unf_ref,:,good), view(normals,:,good))

	# unfold points in random chunks
	ndom = size(to_unf,2)
	shuffled_ids = shuffle(1:ndom)
	nchunks = ceil(Int, ndom/opars.nchunks)
	chunks  = collect(Iterators.partition(shuffled_ids, nchunks))

	for (i, ids) in enumerate(chunks)
		# only reference surface as conditioner
		known_org = ref[:,good]
		known_unf = unf_ref[:,good]
		ids_to_unf = ids

		# add unfolded from previous chunks
		if i > 1
			extra_ids = vcat(chunks[1:(i-1)]...)
			extra_org = view(to_unf,:,extra_ids)
			extra_unf = view(unf,:,extra_ids)
			good_, bad_ = error_ids(extra_org, extra_unf, opars.search,
			                        opars.neigh, opars.maxerr)

			known_org = hcat(known_org, view(extra_org, :, good_))
			known_unf = hcat(known_unf, view(extra_unf, :, good_))
			ids_to_unf = Int.(union(ids, view(extra_ids, bad_)))
		end

		# unfold points
		initguess    = view(unf, :, ids_to_unf)
		chunk_to_unf = view(to_unf, :, ids_to_unf)
		unf[:,ids_to_unf] .= opt(known_org, known_unf, chunk_to_unf,
		                             opars.search, opars.neigh, initguess)
	end

	unf
end

function unfold(to_unf::AbstractVector, ref::AbstractMatrix;
	            isomap=:default, optim=:default)

	merged_to_unf = reduce(hcat, to_unf)
	unf = unfold(merged_to_unf, ref, isomap=isomap, optim=optim)
	j_ = cumsum(size.(to_unf, 2))
	i_ = [i == 1 ? 1 : j_[i-1]+1 for i in 1:length(j_)]
	[unf[:,i:j] for (i,j) in zip(i_,j_)]
end

"""
	unfold(to_unf, ref, unf_ref; optim=:default)

Unfold the input points `to_unf` based on points already unfolded `unf_ref` and
their correspondent in the original space `ref`. Useful to unfold samples after
the domain is already unfolded, for example. Returns a coordinate matrix with
the unfolded points. Or a tuple of matrices if `to_unf` is a list of points.

## Parameters:

* `to_unf`  - coordinate matrix with the points that will be unfolded
* `ref`     - coordinate matrix with the reference points for unfolding
* `unf_ref` - coordinate matrix with the unfolded reference points
* `optim`   - additional parameters for optimization step; either :default or a
  NamedTuple with the keys shown below

### `optim` parameters:

*Default:* optim = (search=:knn, neigh=50, maxerr=5, nchunks=1)
* `search`  - search type to optimize locally (:knn for k-nearest neighbor or
  :radius for radius search).
* `neigh`   - number of neighbors (for search=:knn) or radius distance
  (for search=:radius) to local optimization.
* `maxerr`  - the maximum accepted distortion for unfolding optimization.
* `nchunks` - number of rounds of optimization.
"""
function unfold(to_unf::AbstractMatrix, ref::AbstractMatrix, unf_ref::AbstractMatrix;
	optim=:default)

	# read optim parameters or assign the defaults below
	opars = (search=:knn, neigh=50, maxerr=5, nchunks=2)
	opars = updatepars(opars, optim)

	# conversions if necessary
	!(to_unf[1] isa Float64) && (to_unf = Float64.(to_unf))

	# unfold points in random chunks
	ndom = size(to_unf,2)
	shuffled_ids = shuffle(1:ndom)
	nchunks = ceil(Int, ndom/opars.nchunks)
	chunks  = collect(Iterators.partition(shuffled_ids, nchunks))
	unf = firstguess(to_unf, ref, unf_ref)

	for (i, ids) in enumerate(chunks)
		# only reference surface as conditioner
		known_org = ref
		known_unf = unf_ref
		ids_to_unf = ids

		# add unfolded from previous chunks
		if i > 1
			extra_ids = vcat(chunks[1:(i-1)]...)
			extra_org = view(to_unf,:,extra_ids)
			extra_unf = view(unf,:,extra_ids)
			good_, bad_ = error_ids(extra_org, extra_unf, opars.search,
									opars.neigh, opars.maxerr)

			known_org = hcat(known_org, view(extra_org, :, good_))
			known_unf = hcat(known_unf, view(extra_unf, :, good_))
			ids_to_unf = Int.(union(ids, view(extra_ids, bad_)))
		end

		# unfold points
		initguess    = view(unf, :, ids_to_unf)
		chunk_to_unf = view(to_unf, :, ids_to_unf)
		unf[:,ids_to_unf] .= opt(known_org, known_unf, chunk_to_unf,
									 opars.search, opars.neigh, initguess)
	end
	unf
end

function unfold(to_unf::AbstractVector, ref::AbstractMatrix,
	            unf_ref::AbstractMatrix; optim=:default)

	merged_to_unf = reduce(hcat, to_unf)
	unf = unfold(merged_to_unf, ref, unf_ref, optim=optim)
	j_ = cumsum(size.(to_unf, 2))
	i_ = [i == 1 ? 1 : j_[i-1]+1 for i in 1:length(j_)]
	[unf[:,i:j] for (i,j) in zip(i_,j_)]
end

function updatepars(default, newpars)
	if newpars != :default
		for (k,v) in zip(keys(newpars), values(newpars))
			default = Setfield.setindex(default, v, k)
		end
	end
	default
end