text:
"""
    compute_2d_boundaries_indices

Computes the boundaries and indices for a regular spatial sampling,
think e.g. 2x2 degree grid boxes.

"""

function compute_2d_boundaries_indices(lons::Vector{<:Real},
                                       lats::Array{<:Real},
                                       SS::RegularGridCells)

    return calculate_regular_2d_histogram(
        lons, lats,
        SS.delta_lon,
        SS.delta_lat
    )
end

"""
    compute_2d_boundaries_indices

Computes the boundaries and indices for a full ROI, i.e. no further
aggregation and only one (maybe big) grid cell and all scenes fall
into the same [1] grid cell.

"""
function compute_2d_boundaries_indices(lons::Vector{<:Real},
                                       lats::Vector{<:Real},
                                       SS::FullROI)

    lon_grid = [minimum(lons), maximum(lons)]
    lat_grid = [minimum(lats), maximum(lats)]
    idx = repeat([1], length(lons))

    return (lon_grid, lat_grid), idx
end

"""
    aggregate_scenes(IS, SS, TS)

This function will aggregate an array of scenes according to a
regular time sampler. Regular time sampling means "stick all scenes
in a given regular time interval into a separate aggregation container".
Classic example: collect in daily/weekly/monthly intervals.

# Arguments
- 'IS::InstrumentSampling': An instrument sampling object (contains the scenes)
- 'SS::SpatialSampling': A spatial sampling object
- 'TS::RegularTemporalSampling': Regular temporal sampling object

# Returns
- 'Array{Aggregate, 1}': aggregate objects

"""
function aggregate_scenes(IS::InstrumentSampling,
                          SS::SpatialSampling,
                          TS::RegularTemporalSampling; minimum_number=1)

    # For empty instrument samplings, just return an empty container
    if length(IS) == 0
        return Aggregate[]
    end

    if minimum_number < 1
        @error "Minimum number must be >= 1 (you wanted $(minimum_number))"
        return Aggregate[]
    end

    # Just in case the scenes are not ordered in time, established ordering index
    # Perform the temporal aggregation according to TS
    time_idx = sortperm(get_time(IS))

    # Extract first and last times
    start_date = get_time(IS)[time_idx[1]]
    end_date = get_time(IS)[time_idx[end]]

    # This is the time boundary grid
    time_boundary = collect(start_date:Dates.Second(TS.period):end_date+Dates.Second(TS.period))

    # Compute spatial aggregation here
    spatial_bound, spatial_idx = compute_2d_boundaries_indices(
        get_lon(IS),
        get_lat(IS),
        SS
    )

    # This is an empty list of Aggregate objects, which
    # we fill up in the next section of the code.
    aggregates = Aggregate[]

    # General idea:
    # -------------
    # Go through every unique spatial index, and perform a time aggregate
    # on the subset of all scenes which fall into this spatial index. We
    # then end up with a spatio-temporal aggregation according to SS and TS

    for ss_idx in unique(spatial_idx)

        # Remember - spatial index can be of any arbitrary dimension, but
        # will most likely be 1- or 2-dim tuple

        # Run through all time bins
        last_start_index = 1
        for i in 1:length(time_boundary) - 1

            agg_idx = Int[]

            for j in last_start_index:length(IS)

                if (IS.scenes[time_idx[j]].loctime.time >= time_boundary[i]) &
                    (IS.scenes[time_idx[j]].loctime.time < time_boundary[i+1]) &
                    (ss_idx == spatial_idx[time_idx[j]])

                    # time_idx[j] belongs to aggregate bin bounded by
                    # time_boundary[i] and time_boundary[i+1]

                    # Thus, scene with index "j" belongs to spatial subset "ss_idx"
                    # and temporal bin between time_boundary[i] and time_boundary[i+1]

                    # Save current index, so that next iteration of the scene loop
                    # can start from here
                    last_start_index = j

                    # Push this index into index list
                    push!(agg_idx, time_idx[j])

                end
            end

            # Continue to next bin if below threshold
            if length(agg_idx) < minimum_number
                continue
            end

            # Otherwise, create a new aggregate and
            # push it into list
            this_aggregate = Aggregate(
                length(agg_idx),
                agg_idx,
                IS.scenes[agg_idx],
                time_boundary[i],
                time_boundary[i+1],
            )

            push!(aggregates, this_aggregate)

        end # This ends the temporal index loop
    end # This ends the spatial index loop

    return aggregates

end

function calculate_light_response_curve(v_agg::Vector{Aggregate}; sif_truth=false)

    N = (x -> x.N).(v_agg)

    # Calculate the PPFDs for each aggregate
    ppfd_means = (x -> mean(get_ppfd.(x.scenes))).(v_agg)
    ppfd_std = (x -> std(get_ppfd.(x.scenes))).(v_agg)

    # Obtain the SIF values

    # 1) SIF aggregate mean
    if sif_truth
        sif_means = (x -> mean(get_sif_truth.(x.scenes))).(v_agg)
        sif_std = (x -> std(get_sif_truth.(x.scenes))).(v_agg)
    else
        sif_means = (x -> mean(get_sif.(x.scenes))).(v_agg)
        sif_std = (x -> std(get_sif.(x.scenes))).(v_agg)
    end
    # 2) SIF aggregate uncertainty is calculated using the
    #    per-retrieval uncertainty values
    sif_ucerts = (x -> sqrt(1.0 / sum(1.0 ./ (get_sif_ucert.(x.scenes) .^2)))).(v_agg)

    # 3) Get NIRv (radiance), which is NIR (reflectance) * NDVI * NIR (radiance)
    nirv_means = (x -> mean(get_nirv_radiance.(x.scenes))).(v_agg)

    # 4) Calculate some uncertainty for NIRv (radiance)

    # Assume an uncertainty of 2% for radiance
    # https://www.spiedigitallibrary.org/journals/journal-of-applied-remote-sensing/volume-12/issue-3/034001/Updates-of-Moderate-Resolution-Imaging-Spectroradiometer-on-orbit-calibration-uncertainty/10.1117/1.JRS.12.034001.full

    # These are lists of lists
    # scenes of aggregates
    nir_agg = (x -> get_nir.(x.scenes)).(v_agg)
    vis_agg = (x -> get_vis.(x.scenes)).(v_agg)

    # Construct the per-scene uncertainties for any NIRv
    nirv_ucerts = []

    for agg in 1:length(v_agg)

        # These are now lists (scenes)
        nir = nir_agg[agg]
        vis = vis_agg[agg]

        scene_ucerts = []
        for s in 1:v_agg[agg].N

            # say .. some percent of its value?
            ucert_nir = 0.02 * nir[s]
            ucert_vis = 0.02 * vis[s]

            x = 1.0 / (nir[s] + vis[s])^2
            x *= sqrt(
                4 * nir[s]^4 * ucert_vis^2
                + (nir[s]^2 + vis[s]^2 - 2 * nir[s] * vis[s])^2
                * ucert_nir^2
            )

            this_sza = get_sza(v_agg[agg].scenes[s])
            x *= calculate_BOA_irradiance(this_sza, 865.0) * 1000 / pi

            push!(scene_ucerts, x)

        end

        push!(nirv_ucerts, sqrt(1.0 / sum(1.0 ./ (scene_ucerts .^ 2))))

    end

    # 5) Construct the SIF / NIRv (radiance) ratio along with
    #    the uncertainties

    ratio = sif_means ./ nirv_means

    ratio_ucert = sqrt.((ratio .^2) .* (
        (sif_ucerts ./ sif_means) .^ 2 +
        (nirv_ucerts ./ nirv_means) .^ 2
    ))

    return N, ppfd_means, ratio, ratio_ucert

end


text:
# adjoint-project.jl
# test adjoint consistency for SPECT projector/back-projector on very small case

using SPECTrecon: SPECTplan
using SPECTrecon: project!, backproject!
using SPECTrecon: project, backproject
using LinearMapsAA: LinearMapAA
using LinearAlgebra: dot
using Test: @test, @testset

@testset "adjoint-project-matrix" begin
    T = Float32
    nx = 8; ny = nx
    nz = 6
    nview = 7

    mumap = rand(T, nx, ny, nz)

    px = 5
    pz = 3
    psfs = rand(T, px, pz, ny, nview)
    psfs = psfs .+ mapslices(reverse, psfs, dims = [1, 2]) # symmetrize
    psfs = psfs ./ mapslices(sum, psfs, dims = [1, 2])

    dy = T(4.7952)
    idim = (nx,ny,nz)
    odim = (nx,nz,nview)

    for interpmeth in (:one, :two)
        for mode in (:fast, :mem)
            plan = SPECTplan(mumap, psfs, dy; interpmeth, mode)
            forw! = (y,x) -> project!(y, x, plan)
            back! = (x,y) -> backproject!(x, y, plan)
            A = LinearMapAA(forw!, back!, (prod(odim),prod(idim)); T, idim, odim)
            @test Matrix(A)' ≈ Matrix(A')
        end
    end

    for interpmeth in (:one, :two)
        for mode in (:fast, :mem)
            forw = x -> project(x, mumap, psfs, dy; interpmeth, mode)
            back = y -> backproject(y, mumap, psfs, dy; interpmeth, mode)
            A = LinearMapAA(forw, back, (prod(odim),prod(idim)); T, idim, odim)
            @test Matrix(A)' ≈ Matrix(A')
        end
    end

end

@testset "proj!-adj-test-dot" begin
    T = Float64
    nx = 32; ny = nx
    nz = 20
    nview = 30

    mumap = rand(T, nx, ny, nz)

    px = 5
    pz = 3
    psfs = rand(T, px, pz, ny, nview)
    psfs = psfs .+ mapslices(reverse, psfs, dims = [1, 2])
    psfs = psfs ./ mapslices(sum, psfs, dims = [1, 2])

    dy = T(4.7952)

    plan = SPECTplan(mumap, psfs, dy)
    show(isinteractive() ? stdout : devnull, "text/plain", plan)

    for interpmeth in (:one, :two)
        for mode in (:fast, :mem)
            plan = SPECTplan(mumap, psfs, dy; interpmeth, mode)
            x = randn(T, nx, ny, nz)
            y = randn(T, nx, nz, nview)
            output_x = similar(y)
            output_y = similar(x)
            project!(output_x, x, plan)
            backproject!(output_y, y, plan)
            @test isapprox(dot(y, output_x), dot(x, output_y); rtol = 1e-5)
        end
    end

    for interpmeth in (:one, :two)
        for mode in (:fast, :mem)
            plan = SPECTplan(mumap, psfs, dy; interpmeth, mode)
            x = randn(T, nx, ny, nz)
            y = randn(T, nx, nz, nview)
            output_x = similar(y)
            output_y = similar(x)
            backproject!(output_y, y, plan)
            project!(output_x, x, plan)
            @test isapprox(dot(y, output_x), dot(x, output_y); rtol = 1e-5)
        end
    end

    for interpmeth in (:one, :two)
        for mode in (:fast, :mem)
            x = randn(T, nx, ny, nz)
            y = randn(T, nx, nz, nview)
            output_x = project(x, mumap, psfs, dy; interpmeth, mode)
            output_y = backproject(y, mumap, psfs, dy; interpmeth, mode)
            @test isapprox(dot(y, output_x), dot(x, output_y); rtol = 1e-5)
        end
    end

end

ext: jl


text:
abstract type WeightUpdateStrategy end

# cf. https://discourse.julialang.org/t/broadcasting-structs-as-scalars/14310
Base.Broadcast.broadcastable(ws::WeightUpdateStrategy) = Ref(ws)

struct NoUpdateStrategy <: WeightUpdateStrategy end

update_v(v, label::Symbol, update_strategy::NoUpdateStrategy) = v

struct FixedWeightStrategy <: WeightUpdateStrategy
    v_Lin::Float64
    v_Lout::Float64
end

function update_v(v, fb_label::Symbol, update_strategy::FixedWeightStrategy)
    if fb_label == :Lin
        return update_strategy.v_Lin
    elseif fb_label == :Lout
        return update_strategy.v_Lout
    else
        return v
    end
end

get_default_v(::FixedWeightStrategy) = one(Float64)


text:
using Test, SpecialFunctions
import TruncatedNormal as TN

@test TN.tnmom2c(0, -Inf, Inf) == 1
@test TN.tnmom2c(0, -Inf, 0) == 1
@test TN.tnmom2c(0, 0, +Inf) == 1

@test TN.tnmom2c(0, +Inf, +Inf) == Inf
@test TN.tnmom2c(0, -Inf, -Inf) == Inf

@test isnan(TN.tnmom2c(0, 1, 0))
@test isnan(TN.tnmom2c(0, +Inf, 0))
@test isnan(TN.tnmom2c(0, 0, -Inf))

@test isnan(TN.tnmom2c(0, NaN, NaN))
@test isnan(TN.tnmom2c(0, 0, NaN))
@test isnan(TN.tnmom2c(0, NaN, 0))

@test TN.tnmom2c(10, 0, 10) ≈ 85.0423087839426928824034533318
@test TN.tnmom2c(1, -3, 3) ≈ 1.97333692466254147658812248690
@test TN.tnmom2c(1, 3, 4) ≈ 5.15893137101968604368598313652
@test TN.tnmom2c(-1, 3, 4) ≈ 18.2007485133797706592985387254
@test TN.tnmom2c(2, 3, 4) ≈ 1.63802279983964373587970534208
@test TN.tnmom2c(-2, 3, 4) ≈ 27.7216570845598129671048165198

for x = -10:10
    @test TN.tnmom2c(0, x, +Inf) ≈ 1 + x * √(2/π) / erfcx(+x/√2)
    @test TN.tnmom2c(0, -Inf, x) ≈ 1 - x * √(2/π) / erfcx(-x/√2)
end

@test TN.tnmom2c(0, -Inf, Inf, 0, 1) == 1

for x = -10:10
    @test TN.tnmom2c(0, x, +Inf, 0, 1) ≈ TN.tnmom2(x, +Inf)
    @test TN.tnmom2c(0, -Inf, x, 0, 1) ≈ TN.tnmom2(-Inf, x)
end

@test TN.tnmom2c(0, -1e200, 1e200) ≈ 1
@test TN.tnmom2c(0, -1e6, -999000) ≈ 9.9800100000199999999999799599399200002e11
@test TN.tnmom2c(0, +1e6, +Inf) ≈ 1.00000000000199999999999800000e12
@test TN.tnmom2c(0, -Inf, -1e6) ≈ 1.00000000000199999999999800000e12

@test TN.tnmom2c(4, -10, 20, 3, 3//2) ≈ 3.24999999999999956053727803062


text:
using ArgParse

s=ArgParseSettings()
@add_arg_table s begin
    "--part2"
        help="Do part 2"
        action=:store_true
    "--test"
        help="Run test data"
        action=:store_true
end
parsed_args = parse_args(ARGS, s)

mutable struct Board
    squares::Dict
    columns::Array
    rows::Array
    done::Bool
end

function mark_board(board, num)
    if !board.done && haskey(board.squares, num)
        (row, col) = pop!(board.squares, num)
        board.rows[row] += 1 # Crossed off one more number in this row
        board.columns[col] += 1 # Likewise colum
        board.done = (board.rows[row] == 5 || board.columns[col]==5)
        return board.done
    end
    false
end

function find_winner(numbers, boards)
    for num in numbers
        for board in boards
            if mark_board(board, num)
               return num, board
            end
        end
    end
end

function find_loser(numbers, boards)
    number_ofboards = length(boards)
    for num in numbers
        for board in boards
            if mark_board(board, num)
                number_ofboards -= 1
                if number_ofboards == 0
                    return num, board
               end
            end
        end
    end
end

test_suffix = parsed_args["test"] ? "_test" : ""
lines = readlines("input/day_4$test_suffix.txt")
numbers = split(lines[1], ",")
boards = []
for boardIdx in 2:6:length(lines)-1
    squares = []
    for row in 1:5
        row_nums = split(lines[boardIdx+row])
        for (col, val) in enumerate(row_nums)
            push!(squares, (val,(row, col)))
        end
    end
    push!(boards,Board(Dict(squares), zeros(5), zeros(5), false))
end
if parsed_args["part2"]
    (last, the_one) = find_loser(numbers, boards)
else
    (last, the_one) = find_winner(numbers, boards)
end
unmarked = sum(s -> parse(Int32, s), keys(the_one.squares))
answer = unmarked * parse(Int32, last)
println("The answer is: $answer")


text:
@testset "proximal maps" begin
  #
  # prox_TV
  p = [1.,0.,0.]
  q = [0.,1.,0.]
  M = Sphere(2)
  N = PowerManifold(M, NestedPowerRepresentation(), 2)
  @test_throws ErrorException prox_distance(M, 1., p, q, 3)
  @test distance(M, prox_distance(M, distance(M,p,q)/2, p,q,1), shortest_geodesic(M, p, q, 0.5)) ≈ 0
  (r,s) = prox_TV(M,π/4,(p,q))
  @test norm( r - s ) < eps(Float64)
  # i.e. they are moved together
  @test distance(M,r,s) < eps(Float64)
  (t,u) = prox_TV(M,π/8,(p,q));
  @test_throws ErrorException prox_TV(M, π, (p,q), 3)
  # they cross correlate
  @test ( abs(t[1]-u[2])< eps(Float64) && abs(t[2]-u[1]) < eps(Float64) && abs(t[3]-u[3])< eps(Float64) )
  @test distance(M,t,u) == π/4 # and have moved half their distance
  #
  (v,w) = prox_TV(M,1.,(p,q),2)
  vC, wC = shortest_geodesic(M, p, q, [1/3, 2/3])
  @test distance(M, v, vC) ≈ 0
  @test distance(M, w, wC) ≈ 0
  # prox_TV on Power
  T = prox_TV(N, π/8, [p,q])
  @test distance(N, T, [t, u]) ≈ 0
  # parallelprox_TV
  N2 = PowerManifold(M, NestedPowerRepresentation(), 3)
  r = geodesic(M, p, q, 0.5)
  s, t = prox_TV(M, π/16, (r, q) )
  u, v = prox_TV(M, π/16, (p, r) )
  y = prox_parallel_TV(N2, π/16, [[p,r,q], [p,r,q]])
  @test distance(N2, y[1], [p,s,t]) ≈ 0 # even indices in first comp
  @test distance(N2, y[2], [u,v,q]) ≈ 0 # odd in second
  # dimensions of x have to fit, here they don't
  @test_throws ErrorException prox_parallel_TV(N2, π/16, [[p,r,q]])
  # prox_TV2
  p2, r2, q2 = prox_TV2(M,1., (p,r,q) )
  sum(distance.(Ref(M), [p,r,q], [p2, r2, q2] )) ≈ 0
  @test_throws ErrorException prox_TV2(M, 1., (p,r,q), 2) # since prox_TV is only defined for p=1
  distance(PowerManifold(M,NestedPowerRepresentation(), 3), [p2,r2,q2], prox_TV2(PowerManifold(M, NestedPowerRepresentation(), 3),1., [p,r,q])) ≈ 0
  # Circle
  M2 = Circle()
  N2 = PowerManifold(M2,3)
  pS, rS, qS = [-0.5, 0.1, 0.5]
  d = dot([pS, rS, qS], [1., -2., 1.])
  m = min(0.3, abs( Manopt.sym_rem(d)/6) )
  s = sign(Manopt.sym_rem(d))
  pSc, rSc, qSc = Manopt.sym_rem.( [pS, rS, qS] .- m .* s .* [1., -2., 1.] )
  pSr, rSr, qSr = prox_TV2(M2, 0.3, (pS, rS, qS) )
  @test sum( distance.(Ref(M2), [pSc, rSc, qSc], [pSr, rSr, qSr]) ) ≈ 0
  # p=2
  t = 0.3*Manopt.sym_rem(d)/(1+0.3*6.)
  @test sum(
    distance.( Ref(M2),
      [prox_TV2(M2,0.3,(pS, rS, qS),2)...],
      [pS, rS, qS] .- t.*[1., -2., 1.]
  )) ≈ 0
  # others fail
  @test_throws ErrorException prox_TV2(M2, 0.3, (pS, rS, qS),3)
  # Rn
  M3 = Euclidean(1)
  pR, rR, qR = [pS, rS, qS]
  m = min.( Ref(0.3), abs.( [pR, rR, qR].*[1., -2., 1.] )/6 )
  s = sign(d)  # we can reuse d
  pRc, rRc, qRc = [pR, rR, qR]  .- m .* s .* [1., -2., 1.]
  pRr, rRr, qRr = prox_TV2(M3, 0.3, (pR, rR, qR) )
  @test sum( distance.(Ref(M3), [pRc, rRc, qRc], [pRr, rRr, qRr]) ) ≈ 0
  # p=2
  t = 0.3*d/(1+0.3*6.)
  @test sum(
    distance.( Ref(M3),
      [prox_TV2(M3, 0.3, (pR, rR, qR), 2)...],
      [pR, rR, qR] .- t.*[1., -2., 1.]
  )) ≈ 0
  # others fail
  @test_throws ErrorException prox_TV2(M3, 0.3, (pR, rR, qR),3)
  #
  # collaborative integer tests
  #
  @test_throws ErrorException prox_TV2(M3, 0.3, (pS, rS, qS),3)
  ξR, ηR, νR = [pS, rS, qS]
  N3 = PowerManifold(M3,3)
  P = [pR rR qR]
  Ξ = [ξR ηR νR]
  @test prox_collaborative_TV(N3,0.,P, Ξ,1, 1) == Ξ
  @test prox_collaborative_TV(N3,0.,P, Ξ,1., 1) == Ξ
  @test prox_collaborative_TV(N3,0.,P, Ξ,1, 1.) == Ξ
  @test prox_collaborative_TV(N3,0.,P, Ξ,1., 1.) == Ξ

  @test prox_collaborative_TV(N3,0.,P, Ξ, 2, 1) == Ξ
  @test norm(N3,P, prox_collaborative_TV(N3,0.,P, Ξ, 2, Inf)) ≈ 0
  @test norm(N3,P, prox_collaborative_TV(N3,0.,P, Ξ, 1, Inf)) ≈ 0
  @test norm(N3,P, prox_collaborative_TV(N3,0.,P, Ξ, Inf, Inf)) ≈ 0
  @test_throws ErrorException prox_collaborative_TV(N3,0.,P, Ξ, 3,3)
  @test_throws ErrorException prox_collaborative_TV(N3,0.,P, Ξ, 3,1)
  @test_throws ErrorException prox_collaborative_TV(N3,0.,P, Ξ, 3,Inf)
end


text:
function _fourier_structure!(
    gofr::AbstractArray,
    rdfpos::AbstractArray,
    sqobj::StructureFactor,
    dr::AbstractFloat,
)
    sum_total = 0.0

    for i = 1:sqobj.nm
        sqobj.sq[i, 1] = (i - 1.0) * sqobj.dq
    end

    for i = 2:sqobj.nm
        integrand =
            @. rdfpos^2 * gofr * sin(sqobj.sq[i, 1] * rdfpos) / (sqobj.sq[i, 1] * rdfpos)
        sum_total = sum(integrand)
        sum_total *= sqobj.norm_const * dr
        sqobj.sq[i, 2] = 1.0 + sum_total
    end
end

function compute_structure!(
    sqobj::StructureFactor,
    gofr::PairDistributionFunction,
    s::SimulationSystem,
)
    rdfpos = @view gofr.gofr[2:end, 1]
    hfunc = gofr.gofr[2:end, 2] .- 1.0
    _fourier_structure!(hfunc, rdfpos, sqobj, gofr.dr)
end


text:
using OceanParameterizations
using OceanTurb, Plots
using Oceananigans.Grids: Center, Face

include("lesbrary_data.jl")
include("data_containers.jl")
include("animate_prediction.jl")

output_gif_directory = "Output1"

files =  ["free_convection", "strong_wind", "strong_wind_no_coriolis", "weak_wind_strong_cooling",
          "strong_wind_weak_cooling", "strong_wind_weak_heating"]

for test_file in files
    𝒟test = data(test_file,
                        override_scalings=𝒟train.scalings, # use the scalings from the training data
                        animate=false,
                        animate_dir="$(output_gif_directory)/Testing")

    ## KPP Parameterization (no training)
    les = read_les_output(test_file)
    parameters = KPP.Parameters() # default parameters
    KPP_model = closure_kpp_full_evolution(parameters, 𝒟test.T.coarse[:,1], Δt, les)
    predictions = KPP_model()
    T_KPP = (predictions, 𝒟test.T.coarse)
    mse(T_KPP)
    animate_prediction(T_KPP, "T", 𝒟test, test_file; legend_labels=["KPP(T)", "truth"], filename="T_KPP_$(test_file)", directory=output_gif_directory)

    ## TKE Parameterization (no training; use default parameters)
    les = read_les_output(test_file)
    parameters = TKEMassFlux.TKEParameters() # default parameters
    TKE_model = closure_tke_full_evolution(parameters, 𝒟test.T.coarse[:,1], Δt, les)
    predictions = TKE_model()
    T_TKE = (predictions, 𝒟test.T.coarse)
    mse(T_TKE)
    animate_prediction(T_TKE, "T", 𝒟test, test_file; legend_labels=["TKE(T)", "truth"], filename="T_TKE_$(test_file)", directory=output_gif_directory)
end


text:
# Skipping MacroDefinition: imcdCanvasPutImage(_canvas,_image,_x,_y,_w,_h,_xmin,_xmax,_ymin,_ymax){if(_image->color_space==IM_RGB){if(_image->has_alpha)cdCanvasPutImageRectRGBA(_canvas,_image->width,_image->height,(unsignedchar*)_image->data[0],(unsignedchar*)_image->data[1],(unsignedchar*)_image->data[2],(unsignedchar*)_image->data[3],_x,_y,_w,_h,_xmin,_xmax,_ymin,_ymax);elsecdCanvasPutImageRectRGB(_canvas,_image->width,_image->height,(unsignedchar*)_image->data[0],(unsignedchar*)_image->data[1],(unsignedchar*)_image->data[2],_x,_y,_w,_h,_xmin,_xmax,_ymin,_ymax);}elsecdCanvasPutImageRectMap(_canvas,_image->width,_image->height,(unsignedchar*)_image->data[0],_image->palette,_x,_y,_w,_h,_xmin,_xmax,_ymin,_ymax);}
# begin enum imComplex2Real
typealias imComplex2Real UInt32
const IM_CPX_REAL = 0
const IM_CPX_IMAG = 1
const IM_CPX_MAG = 2
const IM_CPX_PHASE = 3
# end enum imComplex2Real
# begin enum imGammaFactor
typealias imGammaFactor Cint
const IM_GAMMA_LINEAR = 0
const IM_GAMMA_LOGLITE = -10
const IM_GAMMA_LOGHEAVY = -1000
const IM_GAMMA_EXPLITE = 2
const IM_GAMMA_EXPHEAVY = 7
# end enum imGammaFactor
# begin enum imCastMode
typealias imCastMode UInt32
const IM_CAST_MINMAX = 0
const IM_CAST_FIXED = 1
const IM_CAST_DIRECT = 2
const IM_CAST_USER = 3
# end enum imCastMode


text:
abstract type AbstractModel end
abstract type AbstractSingleModel <: AbstractModel end
# needs to be defined for each model

value(m::AbstractModel, ::DateObs) = startdate(m)
value(m::AbstractModel, x::ConstObs) = x.val
value(m::AbstractModel, ::Zero) = 0*numeraire(m)
value(m::AbstractModel, c::Amount) = value(m, c.o)
value(m::AbstractModel, c::Scale{ConstObs{T},C}) where {T,C} = c.s.val * value(m, c.c)

value(m::AbstractModel, c::Either) = max(value(m, c.c1), value(m, c.c2))
value(m::AbstractModel, c::Both) = value(m, c.c1) + value(m, c.c2)
value(m::AbstractModel, c::Give) = -value(m, c.c)

@inline valueat(m::AbstractModel, c::Either, i...) =
    max(valueat(m,c.c1,i...),valueat(m,c.c2,i...))

@inline valueat(m::AbstractModel, c::Zero, i...) =
    0.0*numeraire(m)

@inline valueat(m::AbstractModel, c::Amount, i...) =
    valueat(m, c.o, i...)

@inline valueat(m::AbstractModel, c::Both, i...) =
    valueat(m,c.c1, i...) + valueat(m,c.c2, i...)

@inline valueat(m::AbstractModel, c::Give, i...) =
    -valueat(m,c.c,i...)

@inline valueat(m::AbstractModel, c::Scale, i...) =
    valueat(m,c.s,i...)*valueat(m,c.c,i...)

@inline valueat(m::AbstractModel, o::ConstObs, i...) =
    o.val

value(m::AbstractModel, c::WhenAt{Zero}) = value(m, Zero())

value(m::AbstractModel, c::WhenAt{Both{C1,C2}}) where {C1,C2} =
    value(m, When(c.p, c.c.c1)) + value(m, When(c.p, c.c.c2))
value(m::AbstractModel, c::WhenAt{Give{C}}) where {C} =
    -value(m, When(c.p, c.c.c))


text:
function solve()
    n = parse(Int, readline())
    ret = 0
    for i=1:n-1
        ret += 1 / (1 - i/n)
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
#=
find maximum volume hypercube with edges parallel to the axes inside a polyhedron
or an ellipsoid defined with l_1, l_infty, or l_2 ball constraints
=#

struct MaxVolumeJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    n::Int
    epinormeucl_constr::Bool # add an L2 ball constraint
    epinorminf_constrs::Bool # add L1 and Linfty ball constraints
end

function build(inst::MaxVolumeJuMP{T}) where {T <: Float64}
    @assert xor(inst.epinormeucl_constr, inst.epinorminf_constrs)
    n = inst.n

    A = randn(n, n)
    # ensure there will be a feasible solution
    x = randn(n)
    A = A * A' + 10I
    gamma = norm(A * x) / sqrt(n)

    model = JuMP.Model()
    JuMP.@variable(model, t)
    JuMP.@variable(model, x[1:n])
    JuMP.@objective(model, Max, t)
    JuMP.@constraint(model, vcat(t, x) in MOI.GeometricMeanCone(n + 1))

    if inst.epinormeucl_constr
        JuMP.@constraint(model, vcat(gamma, A * x) in JuMP.SecondOrderCone())
    end
    if inst.epinorminf_constrs
        JuMP.@constraint(model, vcat(gamma, A * x) in MOI.NormInfinityCone(n + 1))
        JuMP.@constraint(model, vcat(sqrt(n) * gamma, A * x) in
            MOI.NormOneCone(n + 1))
    end

    return model
end


text:
@inline splitview(a::Tuple...) = splitview(Float64,a)
@inline splitview{T}(::Type{T},a::Tuple...) = splitview(T,a)
@inline splitview(x::AbstractVector,a::Tuple...) = splitview(x,a)

@generated function splitview{T,N}(::Type{T},arr::NTuple{N,Tuple})
    quote
      len = 0
      @nexprs $N (n)->(s_n = len+1; e_n = len+prod(arr[n]); len = e_n)
      x = Array(T,len)
      X = @ntuple $N (n)->(reshape(view(x,s_n:e_n),arr[n]))
      start = @ntuple $N (n)->(s_n)
      stop = @ntuple $N (n)->(e_n)
      return x,X,start,stop
    end
end

@generated function splitview{T,N}(x::AbstractVector{T},arr::NTuple{N,Tuple})
    quote
      len = 0
      @nexprs $N (n)->(s_n = len+1; e_n = len+prod(arr[n]); len = e_n)
      X = @ntuple $N (n)->(reshape(view(x,s_n:e_n),arr[n]))
      start = @ntuple $N (n)->(s_n)
      stop = @ntuple $N (n)->(e_n)
      return X,start,stop
    end
end


text:
# Script to generate a regular grid of lines in ply format for testing

f = open("grid.ply", "w")

N = 7

write(f, """
ply
format ascii 1.0
element vertex $(N^3)
property float x
property float y
property float z
element color $(N^3)
property float r
property float g
property float b
element edge $(3*N^2*(N-1))
property list uchar int vertex_index
end_header
""")

# Vertex positions
for i=1:N
    for j=1:N
        for k=1:N
            write(f, "$i $j $k\n")
        end
    end
end

# Vertex colors
for i=1:N
    for j=1:N
        for k=1:N
            r = rand()
            gb = rand()
            if gb > r
                gb = r
            end
            write(f, "$r $gb $gb\n")
        end
    end
end

# Edges
strides = [1,N,N*N]
for axis=0:2
    for i=1:N, j=1:N, k=1:N-1
        i1 = sum(strides .* (circshift([i,j,k], axis) - 1))
        i2 = sum(strides .* (circshift([i,j,k+1], axis) - 1))
        write(f, "2 $i1 $i2\n")
    end
end

close(f)

max_stars_repo_name: Vishwas-Venkatachalapathy/travis-lab

max_issues_repo_name: Vishwas-Venkatachalapathy/travis-lab

max_forks_repo_name: Vishwas-Venkatachalapathy/travis-lab


text:
export Mesh3, get_ngauss, get_volume

@doc raw"""
`Mesh` holds data structures for an unstructured mesh.

- `nodes`: a $n_v \times 2$ coordinates array
- `edges`: a $n_{\text{edge}} \times 2$ integer array for edges
- `elems`: a $n_e \times 3$ connectivity matrix, 1-based.
- `nnode`, `nedge`, `nelem`: number of nodes, edges, and elements
- `ndof`: total number of degrees of freedoms
- `conn`: connectivity matrix, `nelems × 3` or `nelems × 6`, depending on whether a linear element or a quadratic element is used.
- `lorder`: order of quadrature rule for line integrals
- `elem_type`: type of the element (P1, P2 or BDM1)

"""
mutable struct Mesh3
    nodes::Array{Float64, 2}
    edges::Array{Int64, 2}
    faces::Array{Int64, 2}
    elems::Array{Int64, 2}
    nnode::Int64
    nedge::Int64
    nface::Int64
    nelem::Int64
    ndof::Int64
    conn::Array{Int64, 2}
    lorder::Int64
    elem_ndof::Int64
    elem_type::FiniteElementType
end

"""
    Mesh3(coords::Array{Float64, 2}, elems::Array{Int64, 2},
        order::Int64 = -1, degree::Union{FiniteElementType,Int64} = 1, lorder::Int64 = -1)

- `degree`: 1 for P1 element, 2 for P2 element
- `order`: Integration order for elements
- `lorder`: Integration order for faces
"""
function Mesh3(coords::Array{Float64, 2}, elems::Array{Int64, 2},
    order::Int64 = -1, degree::Union{FiniteElementType,Int64} = 1, lorder::Int64 = -1)
    @assert length(size(coords))==2 && size(coords,2)==3
    @assert length(size(elems))==2 && size(elems,2)==4
    if degree==P1
        degree = 1
    elseif degree==P2
        degree = 2
    end

    if order==-1
        if degree == 1
            order = 2
        elseif degree == 2
            order = 4
        end
    end

    nnode = size(coords, 1)
    nelem = size(elems,1)
    nedges = zeros(Int64, 1)
    c = coords'[:]
    e = Int32.(elems'[:].- 1)
    edges_ptr = @eval ccall((:init_nnfem_mesh3, $LIBMFEM3), Ptr{Clonglong}, (Ptr{Cdouble}, Cint,
            Ptr{Cint}, Cint, Cint, Cint, Ptr{Clonglong}), $c, Int32(size($coords, 1)), $e, Int32(size($elems,1)),
            Int32($order), Int32($degree), $nedges)
    nedges = nedges[1]
    edges = unsafe_wrap(Array{Int64,1}, edges_ptr, (2nedges,), own=true)
    edges = reshape(edges, nedges, 2)
    elem_dof = Int64(@eval ccall((:mfem_get_elem_ndof3, $LIBMFEM3), Cint, ()))
    conn = zeros(Int64, elem_dof * size(elems, 1))
    @eval ccall((:mfem_get_connectivity3, $LIBMFEM3), Cvoid, (Ptr{Clonglong}, ), $conn)
    ndof = Int64(@eval ccall((:mfem_get_ndof3, $LIBMFEM3), Cint, ()))
    conn = reshape(conn, elem_dof, size(elems, 1))'|>Array
    @eval ccall((:mfem_get_element_to_vertices3, $LIBMFEM3), Cvoid, (Ptr{Clonglong}, ), $elems)

    elem_type = missing
    if degree==1
        elem_type = P1
    elseif degree==2
        elem_type = P2
    end

    # get faces
    fset = Set{Tuple{Int64, Int64, Int64}}([])
    for i = 1:nelem
        cc = elems[i,:]
        push!(fset, Tuple(sort(cc[[1;2;3]])))
        push!(fset, Tuple(sort(cc[[1;2;4]])))
        push!(fset, Tuple(sort(cc[[1;3;4]])))
        push!(fset, Tuple(sort(cc[[2;3;4]])))
    end

    faces = vcat([[x[1] x[2] x[3]] for x in fset]...)
    nfaces = size(faces, 1)
    Mesh3(coords, edges,  faces, elems, nnode, nedges, nfaces, nelem, ndof, conn, lorder, elem_dof, elem_type)
end

Base.:copy(mesh::Mesh3) = Mesh3(copy(mesh.nodes),
                            copy(mesh.edges),
                            copy(mesh.faces),
                            copy(mesh.elems),
                            copy(mesh.nnode),
                            copy(mesh.nedge),
                            copy(mesh.nelem),
                            copy(mesh.nface),
                            copy(mesh.ndof),
                            copy(mesh.conn),
                            copy(mesh.lorder),
                            copy(mesh.elem_ndof),
                            copy(mesh.elem_type))

@doc raw"""
    Mesh3(m::Int64, n::Int64, l::Int64,
        h::Float64;
        order::Int64 = 2,
        degree::Union{FiniteElementType, Int64} = 1,
        lorder::Int64 = -1)

Constructs a mesh of a rectangular domain. The rectangle is split into $m\times n$ cells, and each cell is further split into two triangles.
`order` specifies the quadrature rule order. `degree` determines the degree for finite element basis functions.
"""
function Mesh3(m::Int64, n::Int64, l::Int64, h::Float64; order::Int64 = -1,
            degree::Union{FiniteElementType, Int64} = 1, lorder::Int64 = -1)
    coords = zeros((m+1)*(n+1)*(l+1), 3)
    elems = zeros(Int64, 5*m*n*l, 4)
    function ID(i, j, k)
        (k-1)*(n+1)*(m+1) + (j-1)*(m+1) + i
    end
    TE1 = [
        [1; 2; 3; 5],
        [2; 3; 4; 8],
        [3; 5; 7; 8],
        [2; 3; 5; 8],
        [2; 5; 6; 8]
    ]
    TE2 = [
        [1; 2; 4; 6],
        [1; 5; 6; 7],
        [4; 6; 7; 8],
        [1; 4; 6; 7],
        [1; 3; 4; 7]
    ]

    s = 1
    for k = 1:l+1
        for j = 1:m+1
            for i = 1:n+1
                x = (i-1)*h
                y = (j-1)*h
                z = (k-1)*h
                coords[s,:] = [x;y;z]
                s += 1
            end
        end
    end
    K = 0
    for i = 1:n
        for j = 1:m
            for k = 1:l
                IDX = [
                    ID(i, j, k)
                    ID(i+1, j, k)
                    ID(i, j+1, k)
                    ID(i+1, j+1, k)
                    ID(i, j, k+1)
                    ID(i+1, j, k+1)
                    ID(i, j+1, k+1)
                    ID(i+1, j+1, k+1)
                ]
                for s = 1:5
                    if (i+j+k)%2==0
                        elems[s + K,:] = IDX[TE1[s]]
                    else
                        elems[s + K,:] = IDX[TE2[s]]
                    end
                end
                K += 5
            end
        end
    end

    Mesh3(coords, elems, order, degree, lorder)
end

"""
    Mesh3(filename::String; file_format::Union{String, Missing} = missing,
    order::Int64 = 2, degree::Union{FiniteElementType, Int64} = 1, lorder::Int64 = 2)
"""
function Mesh3(filename::String; file_format::Union{String, Missing} = missing,
    order::Int64 = 2, degree::Union{FiniteElementType, Int64} = 1, lorder::Int64 = 2)
    if splitext(filename)[2] == ".mat"
        d = matread(filename)
        return Mesh3(Float64.(d["nodes"]), Int64.(d["elems"]), order, degree, lorder)
    end
    meshio = get_meshio()
    if !ismissing(file_format)
        mesh = meshio.read(filename, file_format = file_format)
    else
    mesh = meshio.read(filename)
    end
    elem = []
    for (mkr, dat) in mesh.cells
        if mkr == "tetra"
            push!(elem, dat)
        end
    end
    elem = vcat(elem...)
    if length(elem)==0
        error("No triangles found in the mesh file.")
    end
    Mesh3(Float64.(mesh.points), Int64.(elem) .+ 1, order, degree, lorder)
end

# """
#     Mesh(; order::Int64 = -1,
#     degree::Union{FiniteElementType, Int64} = 1, lorder::Int64 = -1)

# Creates a mesh with a reference triangle.

# ![](https://raw.githubusercontent.com/ADCMEMarket/ADCMEImages/master/AdFem/mapping.png)

# """
# function Mesh(; order::Int64 = -1,
#     degree::Union{FiniteElementType, Int64} = 1, lorder::Int64 = -1)
#     coords = [
#         1.0 0.0
#         0.0 1.0
#         0.0 0.0
#     ]
#     elems = [
#         1 2 3
#     ]
#     Mesh(coords, elems, order, degree, lorder)
# end

"""
    get_ngauss(mesh::Mesh3)

Return the total number of Gauss points.
"""
function get_ngauss(mesh::Mesh3)
    return Int64(@eval ccall((:mfem_get_ngauss3, $LIBMFEM3), Cint, ()))
end

"""
    get_volume(mesh::Mesh3)

Return the areas of triangles as an array.
"""
function get_volume(mesh::Mesh3)
    a = zeros(size(mesh.elems,1))
    @eval ccall((:mfem_get_area3, $LIBMFEM3), Cvoid, (Ptr{Cdouble}, ), $a)
    a
end

"""
    gauss_nodes(mesh::Mesh3)
"""
function gauss_nodes(mesh::Mesh3)
    ngauss = get_ngauss(mesh)
    x = zeros(ngauss)
    y = zeros(ngauss)
    z = zeros(ngauss)
    @eval ccall((:mfem_get_gauss3, $LIBMFEM3), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), $x, $y, $z)
    [x y z]
end

"""
    fem_nodes(mesh::Mesh3)
"""
function fem_nodes(mesh::Mesh3)
    mesh.nodes
end

"""
    fvm_nodes(mesh::Mesh3)
"""
function fvm_nodes(mesh::Mesh3)
    nnode = size(mesh.nodes, 1)
    nelem = size(mesh.elems, 1)
    out = zeros(nelem, 2)
    for i = 1:nelem
        idx = mesh.elems[i, :]
        out[i, 1] = mean(mesh.nodes[idx, 1])
        out[i, 2] = mean(mesh.nodes[idx, 2])
    end
    return out
end

function _edge_dict(mesh::Mesh3)
    D = Dict{Tuple{Int64, Int64}, Int64}()
    for i = 1:mesh.nedge
        D[(mesh.edges[i,1], mesh.edges[i,2])] = i
        D[(mesh.edges[i,2], mesh.edges[i,1])] = i
    end
    D
end


text:
function subplot_defaults()
    textSize = 8;
    return Dict{Symbol, Any}([
        :xticklabelsize => textSize,
        :yticklabelsize => textSize,
        :xlabelsize => textSize,
        :ylabelsize => textSize
    ])
end

"""
	$(SIGNATURES)

Make a plot that contains several subplots.

# Arguments

- `plotFctV`:  Vector of functions that generate each plot.
    Function signature: `ax = plotFct(fig, pos; kwargs...)`.
    If `plotFct` does not return an `Axis`, `axV` remains empty.
- `figTitle`: Title for entire figure. `nothing` is ignored.
"""
function subplots(plotFctV; figTitle = nothing, kwargs...)
    args = merge(subplot_defaults(), kwargs);
    fig = blank_plot(; args...);
    nPlots = length(plotFctV);
    nRows, nCols = subplot_layout(nPlots);
    axV = Vector{Axis}(undef, nPlots);
    for (j, plotFct) in enumerate(plotFctV)
        pos = subplot_pos(j, nRows, nCols);
        ax = plotFct(fig, pos; color = main_color(), args...);
        if ax isa Axis
            axV[j] = ax;
        end
    end
    add_title!(fig, figTitle);
    return fig, axV
end

"""
	$(SIGNATURES)

Makes Axis objects so that subplots can be written into Figure fig.

It is generally better to make axes one by one.
"""
function make_axes!(fig, nPlots :: Integer; kwargs...)
    args = merge(subplot_defaults(), kwargs);
    nr, nc = subplot_layout(nPlots);
    axV = Vector{Any}(undef, nPlots);
    for ir = 1 : nr
        for ic = 1 : nc
            if (ir * ic) <= nPlots
                j = (ir-1) * nc + ic;
                axV[j] = fig[ir, ic] = Axis(fig; args...);
            end
        end
    end
    return axV
end

"""
	$(SIGNATURES)

Make layout for subplots of `nPlots` equally sized plots.

Return number of rows and columns.

# Example

```julia
fig = blank_plot();
nPlots = 5;
nRows, nCols = subplot_layout(nPlots);
for iPlot = 1 : nPlots
    pos = subplot_pos(iPlot, nRows, nCols);
    line_plot(1:5, 1:5; fig = fig, pos = pos);
end
```
"""
function subplot_layout(nPlots)
    if nPlots <= 3
        nRows = 1;
        nCols = nPlots;
    elseif nPlots == 4
        nRows = 2;
        nCols = 2;
    else
        nCols = 3;
        nRows = ceil(Int, nPlots / nCols);
        @assert (nPlots <= (nRows * nCols)) "$nPlots too high: $nRows, $nCols"
    # else
    #     error("Not implemented for $nPlots");
    end
    return nRows, nCols
end

"""
	$(SIGNATURES)

Row and column position of a linear plot index. Goes together with [`subplot_layout`](@ref).
"""
function subplot_pos(j, nr, nc)
    # Columns first so that we go row by row.
    cIdx = CartesianIndices((nc, nr))[j];
    ir = cIdx[2];
    ic = cIdx[1];
    return ir, ic
end

"""
	$(SIGNATURES)

Link the given axes.
"""
link_axes(args...) = linkaxes!(args...);

"""
	$(SIGNATURES)

Add a title to an entire figure with subplots. This shifts all existing layouts down.

Possible bug: This narrows the figure to the width of the title, unless there are multiple figure columns. Don't use for a single plot. There one should use an axis title.
"""
function add_title!(fig, titleStr; kwargs...)
    fig[0, :] = Label(fig, titleStr; kwargs...);
end

add_title!(fig, titleStr :: Nothing; kwargs...) = nothing;

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


text:
# See https://jp.mathworks.com/help/matlab/ref/peaks.html
function peaks(N=49::Int)
    function formula(x::T, y::T) where T<:Float64
        z = 3(1-x)^2*exp(-(x^2)-(y+1)^2)-10(x/5 - x^3 - y^5)*exp(-x^2-y^2)-(1/3)exp(-(x+1)^2-y^2)
    end
    vec = collect(Float64, LinRange(-3.,3.,N))
    x = repeat(vec', outer=(N,1))
    y = repeat(vec, outer=(1,N))
    z = formula.(x,y)
    return (x,y,z)
end


text:
#########################################################################
#
#   ExGraph type definition and related functions
#
#########################################################################

typealias NSMap BiDict{ExNode, Any}   # ExNode - Symbol map

#####  ExGraph type definitions ######
type ExGraph
  nodes::Vector{ExNode}  # nodes in this graph
  exti::NSMap
  seti::NSMap
  exto::NSMap
  seto::NSMap
end

ExGraph()                   = ExGraph( ExNode[] )
ExGraph(vn::Vector{ExNode}) = ExGraph( vn, NSMap(),
                                           NSMap(),
                                           NSMap(),
                                           NSMap() )

hasnode(m::NSMap, n::ExNode) = haskey(m.kv, n)
hassym( m::NSMap, k)         = haskey(m.vk, k)
nodes(m::NSMap)              = keys(m.kv)
syms( m::NSMap)              = keys(m.vk)

getnode(m::NSMap, k)         = m.vk[k]
getsym( m::NSMap, n::ExNode) = m.kv[n]

function show(io::IO, g::ExGraph)
  tn = fill("", length(g.nodes), 8)

  for (i,n) in enumerate(g.nodes)
    tn[i,1] = "$i"

    if hasnode(g.exti, n)
      sym = g.exti[n]
      tn[i,2] = "$sym >>"
      hassym(g.exto, sym) && (tn[i,3] = "+")
    elseif hasnode(g.seti, n)
      sym = g.seti[n]
      tn[i,2] = "$sym <<"
      hassym(g.seto, sym) && (tn[i,3] = "+")
    end

    tn[i,4] = "[$(subtype(n))]"
    main = isa(n, NFor) ? n.main[1] : n.main

    tn[i,5] = join( map( x -> "$x", indexin(n.parents,    Any[g.nodes...])), ", ")
    tn[i,6] = join( map( x -> "$x", indexin(n.precedence, Any[g.nodes...])), ", ")

    tn[i,7] = repr(main)
    tn[i,8] = "($(typeof(n.val))) $(repr(n.val)[1:min(40, end)])"
  end

  tn = vcat(["node" "symbol" "ext ?" "type" "parents" "precedence" "main" "value"],
        tn)
  sz = maximum(map(length, tn), 1)
  tn = vcat(vec(tn[1,:])', map(s->"-"^s, sz), tn[2:end,:])
  for i in 1:size(tn,1)
    for j in 1:size(tn, 2)
      print(io, rpad(tn[i,j], sz[j]), " | ")
    end
    println(io)
  end

end

#####  ExGraph functions  #####

### checks if n is an ancestor of 'exits', optionnally excluding some nodes from the path
##  !! assumes that g has been evalsorted
function isancestor(anc::ExNode, exits::Vector{ExNode}, g::ExGraph, except::Vector{ExNode})
    gm = Set(exits)
    for n in reverse( g.nodes )
        n in gm     || continue
        n in except && continue
        anc in n.parents && return true
        union!(gm, n.parents)
    end

    false
end

isancestor(n::ExNode, exits::Vector{ExNode}, g::ExGraph) = isancestor(n, exits, g, ExNode[])

# copies a graph and its nodes, leaves onodes references intact
function copy(g::ExGraph)
  g2 = ExGraph()
  nmap = Dict()
  evalsort!(g)
  for n in g.nodes
    n2 = addnode!(g2, copy(n))
    n2.parents    = [ nmap[n] for n in n2.parents    ]
    n2.precedence = [ nmap[n] for n in n2.precedence ]
    nmap[n] = n2
  end

  # update onodes of subgraphs (for loops)
  for n in filter(x -> isa(x, NFor), g2.nodes)
    fg = n.main[2]
    no = NSMap()
    for (k,v) in fg.exto
      no[ nmap[k] ] = v
    end
    fg.exto = no

    no = NSMap()
    for (k,v) in fg.seto
      no[ nmap[k] ] = v
    end
    fg.seto = no
  end

  # copy node mapping and translate inner nodes to newly created ones
  g2.exto = NSMap(g.exto.kv)
  g2.seto = NSMap(g.seto.kv)
  for (k,v) in g.exti
    g2.exti[ nmap[k] ] = v
  end
  for (k,v) in g.seti
    g2.seti[ nmap[k] ] = v
  end

  g2
end

# add a single node
addnode!(g::ExGraph, nn::ExNode) = ( push!(g.nodes, nn) ; return g.nodes[end] )

######## transforms n-ary +, *, max, min, sum, etc...  into binary ops  ######
function splitnary!(g::ExGraph)
  for n in g.nodes  # n = g.nodes[5]
      if isa(n, NCall) && ( length(n.parents) > 3 )
        func = n.parents[1].main
        if isa(n.parents[1], NDot) # if getfield, evaluate function (issue #25)
          fsym = isa(n.parents[1].main, Expr) ?
                     n.parents[1].main.args[1] :  # Expression
                     n.parents[1].main.value      # QuoteNode

          func = getfield(n.parents[1].parents[1].main, fsym)
        end

        if func in [+, *, sum, min, max]
          nn = addnode!(g, NCall(:call, [n.parents[1]; n.parents[3:end]] ) )
          n.parents = [n.parents[1:2]; nn]
        end

      elseif isa(n, NFor)
        splitnary!(n.main[2])

      end
  end
  g
end

####### fuses nodes nr and nk, keeps nk ########
# removes node nr and keeps node nk
#  updates all references to nr
function fusenodes(g::ExGraph, nk::ExNode, nr::ExNode)
  if hasnode(g.exti, nr)
    g.exti[nk] = g.exti[nr]
  end

  # test if nr is associated to a setting node
  # if true, we create an NIn on nk, and associate var to it
  if hasnode(g.seti, nr)
    nn = addnode!(g, NIn(g.seti[nr], [nk]))
    nn.val = "fuse #1"
    g.seti[nn] = g.seti[nr]  # nn replaces nr as set_inode

    if hasnode(g.seto, nr)   # change onodes too (if we are in a subgraph)
      g.seto[nn] = g.seto[nr]  # nn replaces nr as set_onode
    end
  end

  # test if nr is in the precedence of some node and nk is a parent of the same node
  ps = filter(x -> (nr in x.precedence) && (nk in x.parents), g.nodes)
  if length(ps) > 0
    nn = addnode!(g, NIn(nothing, [nk]))
    nn.val = "fuse #2"

    for n in g.nodes
      n.parents    = map(x -> x==nr ? nn : x, n.parents)
      n.precedence = map(x -> x==nr ? nn : x, n.precedence)
    end

  end

  # replace references to nr by nk in parents of other nodes
  for n in filter(n -> n != nr && n != nk, g.nodes)
    if isa(n, NFor)
      g2 = n.main[2]

      # this should not happen...
      @assert !hasnode(g2.seto, nr) "[fusenodes (for)] attempt to fuse set_onode $nr"

      if hasnode(g2.exto, nr)
        if hasnode(g2.exto, nk)  # both nr and nk are used by the for loop
          symr = g2.exto[nr]
          symk = g2.exto[nk]
          fusenodes(g2, getnode(g2.exti, symk), getnode(g2.exti, symr))
        end

        g2.exto[nk] = g2.exto[nr]  # nk replaces nr in g2.exto
      end
    end

    for (i, n2) in enumerate(n.parents)
      n2 == nr && (n.parents[i] = nk)
    end
    for (i, n2) in enumerate(n.precedence)
      n2 == nr && (n.precedence[i] = nk)
    end

  end

  # remove node nr in g
  filter!(n -> n != nr, g.nodes)
end

####### trims the graph to necessary nodes for exitnodes to evaluate  ###########
prune!(g::ExGraph) = prune!(g, collect(nodes(g.seti)))

function prune!(g::ExGraph, exitnodes)
  ns2 = copy(exitnodes)
  evalsort!(g)
  for n in reverse(g.nodes)
    # removed NIn nodes should be removed from For loops too
    if isa(n, NIn) && !(n in ns2) && isa(n.parents[1], NFor)
      fg = n.parents[1].main[2]
      sym = fg.seto[n]
      delete!(fg.seto, n)
    end

    n in ns2 || continue

    if isa(n, NFor)
      g2 = n.main[2]

      # list of g2 nodes whose outer node is in ns2
      exitnodes2 = ExNode[]
      for (k, sym) in g2.seto
        k in ns2 || continue
        push!(exitnodes2, getnode(g2.seti, sym))
      end
      # don't forget reentrant variables
      for (n2, sym) in g2.seti
        (n2 in exitnodes2)               && continue
        hassym(g2.exti, sym)             || continue
        on = getnode(g2.exti, sym)
        isancestor(on, exitnodes2, g2)   || continue
        push!(exitnodes2, n2)
      end
      exitnodes2 = unique(exitnodes2)
      prune!(g2, exitnodes2)

      # update parents
      n.parents = [n.parents[1]; intersect( n.parents, collect(nodes(g2.exto)) ) ]
    end

    # ns2 = union(ns2, n.parents)
    ns2 = vcat(ns2, n.parents)
  end

  # remove unused external inodes in map and corresponding onodes (if they exist)
  for (k,v) in g.exti
    k in ns2 && continue
    delete!(g.exti, k)
    hassym(g.exto, v) && delete!(g.exto, getnode(g.exto, v))
  end
  # reduce set inodes to what was specified in initial exitnodes parameter
  for (k,v) in g.seti
    k in exitnodes && continue
    delete!(g.seti, k)
    hassym(g.seto, v) && delete!(g.seto, getnode(g.seto, v))
  end

  # remove precedence nodes that have disappeared
  for n in ns2
    n.precedence = intersect(n.precedence, ns2)
  end

  # g.nodes = intersect(g.nodes, ns2)
  g.nodes = unique(ns2)
  g
end

####### sort graph to an evaluable order ###########
function evalsort!(g::ExGraph)
  nr = Set{ExNode}(g.nodes)
  ns = ExNode[]

  while length(nr) > 0
    canary = length(ns)
    for n in nr
      any(x -> x in nr, n.parents) && continue
      any(x -> x in nr, n.precedence) && continue
      push!(ns,n) ; delete!(nr,n)
    end
    (canary == length(ns)) && error("[evalsort!] cycle in graph")
  end
  g.nodes = ns

  # separate pass on subgraphs
  map( n -> evalsort!(n.main[2]), filter(n->isa(n, NFor), g.nodes))

  g
end

####### calculate the value of each node  ###########
function calc!(g::ExGraph; params=Dict(), emod = Main)

  # evaluate a symbol
  function myeval(thing)
    local ret

    if haskey(params, thing)
      return params[thing]
    else
      try
        ret = emod.eval(thing)
      catch e
        println("[calc!] can't evaluate $thing in \n $g \n with")
        display(params)
        rethrow(e)
      end
      return ret
    end
  end

  function evaluate(n::NCall)
    local ret
    try
      ret = (n.parents[1].main)([ x.val for x in n.parents[2:end]]...)
    catch e
      eex = Expr(:call, n.parents[1].main, [ x.val for x in n.parents[2:end]]...)
      error("$e when calling $eex in \n$g")
    end
    return ret
  end

  function evaluate(n::NComp)
    local ret
    try
      ret = emod.eval( Expr(:call, n.main, Any[ x.val for x in n.parents]...) )
      #TODO: improve speed with smth like ret = (n.main)([ x.val for x in n.parents]...)
    catch e
      eex = Expr(:call, n.main, Any[ x.val for x in n.parents]...)
      error("$e when calling $eex in \n$g")
    end
    return ret
  end

  function evaluate(n::NExt)
    hasnode(g.exti, n) || return myeval(n.main)

    sym = g.exti[n]  # should be equal to n.main but just to be sure..
    hassym(g.exto, sym) || return myeval(n.main)
    return getnode(g.exto, sym).val  # return node val in parent graph
  end

  evaluate(n::NConst) = n.main

  # evaluate(n::NRef)   = myeval( Expr(:ref , Any[ x.val for x in n.parents]...))
  evaluate(n::NRef)   = getindex(n.parents[1].val, [n2.val for n2 in n.parents[2:end]]...)

  # evaluate(n::NDot)   = myeval( Expr(:.   , n.parents[1].val, n.main) )
  function evaluate(n::NDot)
    fsym = isa(n.main, Expr) ? n.main.args[1] : n.main.value
    getfield(n.parents[1].val, fsym)
  end

  function evaluate(n::NSRef)
    if length(n.parents) >= 3   # regular setindex
      setindex!(n.parents[1].val, n.parents[2].val, [n2.val for n2 in n.parents[3:end]]...)
      return n.parents[1].val
    else
      return n.parents[2].val
    end
  end

  # evaluate(n::NSRef)  = n.parents[1].val
  evaluate(n::NSDot)  = n.parents[1].val

  function evaluate(n::NIn)
      isa(n.parents[1], NFor) && return n.parents[1].val[n]
      n.parents[1].val
  end

  function evaluate(n::NFor)
    g2 = n.main[2]
    is = n.main[1]                    # symbol of loop index
    iter = evaluate(n.parents[1])     #  myeval(n.main[1].args[2])

    params2 = copy(params)
    for is0 in iter
      params2[is] = is0
      calc!(g2, params=params2)
    end

    valdict = Dict()
    for (k, sym) in g2.seto
      valdict[k] = getnode(g2.seti, sym).val
    end

    valdict
  end

  evalsort!(g)
  for (i,n) in enumerate(g.nodes)
    n.val = evaluate(n)
  end
  g
end

###### inserts graph src into dest  ######
# TODO : inserted graph may update variables and necessitate a precedence update

function addgraph!(src::Expr, dest::ExGraph, smap::Dict)
  addgraph!(tograph(src), dest, smap)
end

function addgraph!(src::ExGraph, dest::ExGraph, smap::Dict)
  length(src.exto.kv)>0 && warn("[addgraph] adding graph with external onodes")
  length(src.seto.kv)>0 && warn("[addgraph] adding graph with set onodes")
  # TODO : this control should be done at the deriv_rules.jl level

  ig = copy(src) # make a copy, update references, sort
  exitnode = getnode(ig.seti, nothing) # result of added subgraph

  # evalsort!(ig)

  nmap = Dict()
  for n in ig.nodes  #  n = src[1]
    if isa(n, NExt)
      if haskey(smap, n.main)
        nmap[n] = smap[n.main]

        # should the exitnode be a NExt, it has to be updated
        if n == exitnode
          exitnode = nmap[n]
        end

      else
        # error message removed to leave a chance for resolution
        #   in evalmod (π for example)
        # error("unmapped symbol in source graph $(n.main)")
      end

    elseif isa(n, NFor) # update references, including onodes
      n.parents =    [ haskey(nmap, n2) ? nmap[n2] : n2 for n2 in n.parents    ]
      n.precedence = [ haskey(nmap, n2) ? nmap[n2] : n2 for n2 in n.precedence ]

      g2 = n.main[2]
      for (n2, sym) in g2.exto
        haskey(nmap, n2) || continue
        g2.exto[ nmap[n2]] = sym
      end

      push!(dest.nodes, n)

    else  # update references to NExt that have been remapped
      n.parents =    [ haskey(nmap, n2) ? nmap[n2] : n2 for n2 in n.parents    ]
      n.precedence = [ haskey(nmap, n2) ? nmap[n2] : n2 for n2 in n.precedence ]
      push!(dest.nodes, n)

    end
  end

  # return exitnode of subgraph
  exitnode
end


text:
import MCMCChains: Chains

function create_a3d(noofsamples, noofvariables, noofchains)
   a3d = fill(0.0, noofsamples, noofvariables, noofchains)
   a3d
 end

 function insert_chain!(a3d, chain, posterior, trans)
   for i in 1:size(a3d, 1)
     a3d[i,:,chain] = inverse(trans, posterior[i])
   end
 end

 function insert_chain!(a3d, chain, posterior)
   for i in 1:size(a3d, 1)
     a3d[i,:,chain] = posterior[i, :]
   end
 end

 function create_mcmcchains(a3d, cnames;start=1)
   Chains(a3d, cnames; start=start)
 end

 function create_mcmcchains(a3d, cnames, sections::Dict{Symbol, Vector{String}};
   start=1)
   Chains(a3d, cnames, sections; start=start)
 end


text:
# MetidaFreq.jl

struct ConTab{T <: AbstractMatrix{Int}, R <: Union{Vector{String}, Nothing}, C <: Union{Vector{String}, Nothing}, ID <: Dict} <: AbstractIdData
    tab::T
    rown::R
    coln::C
    id::ID
end

"""
    contab(m::AbstractMatrix{Int};
        rownames::Union{Vector{String}, Nothing} = nothing,
        colnames::Union{Vector{String}, Nothing} = nothing,
        id::Dict = Dict())
"""
function contab(m::AbstractMatrix{Int};
    rownames::Union{Vector{String}, Nothing} = nothing,
    colnames::Union{Vector{String}, Nothing} = nothing,
    id::Dict = Dict())
    ConTab{typeof(m), typeof(rownames), typeof(colnames), typeof(id)}(m, rownames, colnames, id)
end
"""
    contab(v::AbstractVector{Int};
        rownames::Union{Vector{String}, Nothing} = nothing,
        colnames::Union{Vector{String}, Nothing} = nothing,
        id::Dict = Dict())
"""
function contab(v::AbstractVector{Int};
    rownames::Union{Vector{String}, Nothing} = nothing,
    colnames::Union{Vector{String}, Nothing} = nothing,
    id::Dict = Dict())
    contab(permutedims(v); rownames = rownames, colnames = colnames, id = id)
end
"""
    contab(data, row::Symbol, col::Symbol; sort::Union{Nothing, Symbol, AbstractVector{Symbol}} = nothing, id = nothing)
"""
function contab(data, row::Symbol, col::Symbol; sort::Union{Nothing, Symbol, AbstractVector{Symbol}} = nothing, id = nothing)
    cols = Tables.columns(data)
    c    = isnothing(sort) ? (row, col) : isa(sort, Symbol) ? (row, col, sort) : Tuple(append!([row, col], sort))
    res  = contab_(Tuple(Tables.getcolumn(cols, y) for y in c))
    s    = size(res[1])
    rdi = res[2][1]
    cdi = res[2][2]
    rowstr = Vector{String}(undef, length(rdi))
    colstr = Vector{String}(undef, length(cdi))
    if length(s) > 2
        dimn  = length(s) - 2
        ci    = CartesianIndices(s[3:end]) #?
        v     = Vector{ConTab}(undef, length(ci))
        rdict = Vector{Dict}(undef, dimn)
        for i = 1:length(s)-2
            rdict[i] = Dict(v => k for (k,v) in res[2][2+i])
        end
        for i = 1:length(ci)
            id = Vector{Pair}(undef, length(ci[i]))
            for j = 1:length(ci[i])
                id[j] = c[2+j] => rdict[j][ci[i][j]]
            end
            for k in keys(rdi)
                rowstr[rdi[k]] = String(k)
            end
            for k in keys(cdi)
                colstr[cdi[k]] = String(k)
            end
            v[i] = ConTab(Matrix(view(res[1], :, :, ci[i])), rowstr,  colstr, Dict(id))
        end
        return DataSet(v)
    else
        for k in keys(rdi)
            rowstr[rdi[k]] = string(k)
        end
        for k in keys(cdi)
            colstr[cdi[k]] = string(k)
        end
        return ConTab(res[1], rowstr,  colstr, (isnothing(id) ? Dict() : id))
    end
end

function contab_(data::Tuple, T::Type = promote_type(eltype.(data)...))
    #d = Dict{Tuple{eltype.(data)...}, Int}()
    d = Dict{Tuple{map(eltype, data)...}, Int}()
    for (i, element) in enumerate(zip(data...))
        ind = ht_keyindex(d, element)
        if ind > 0
            @inbounds d.vals[ind] += one(Int)
        else
            @inbounds d[element] = one(Int)
        end
    end
    k    = collect(keys(d))
    n    = length(data)
    dims = Vector{Dict{T, Int}}(undef, n)
    @inbounds for i in Base.OneTo(n)
        dims[i] = Dict{eltype(data[i]), Int}()
        s = Set{eltype(data[i])}()
        @inbounds for j in Base.OneTo(length(k))
            push!(s, k[j][i])
        end
        us = collect(s)
        @inbounds for v in Base.OneTo(length(us))
            dims[i][us[v]] = v
        end
    end

    m = zeros(Int, length.(dims)...)
    @inbounds for j in Base.OneTo(length(k))
        m[map(getindex, dims, k[j])...] = d[k[j]]
    end

    #=
    diml = Tuple(map(length, dims))
    li   = LinearIndices(diml)
    m = zeros(Int, length(li))
    @inbounds for j in Base.OneTo(length(k))
        m[li[map(getindex, dims, k[j])...]] = d[k[j]]
    end
    reshape(m, diml...), dims
    =#
    m, dims
end

function getinvindex(data)
    data.pool.invindex
end

function contab_(data::NTuple{n, AbstractCategoricalVector}) where n

    levs = map(levels, data)
    dims = Tuple(map(length, levs))

    #ci = CartesianIndices(dims)
    #m  = zeros(Int, dims...)
    li = LinearIndices(dims)
    a = zeros(Int, length(li))

    r = Vector{Int}(undef, n)
    @inbounds for i in 1:length(data[1])
        @inbounds for j = 1:n
            r[j] = Int(data[j].refs[i])
        end
        a[li[r...]] += 1
        #m[r...] += 1
    end
    #m = reshape(a, dims...)
    reshape(a, dims...), getinvindex.(data)
end

function Base.size(contab::ConTab)
    size(contab.tab)
end
function Base.size(contab::ConTab, dim::Int)
    size(contab.tab, dim)
end
function Base.show(io::IO, contab::ConTab)
    println(io, "  Contingency table::")
    PrettyTables.pretty_table(io, contab.tab; header = contab.coln, row_names = contab.rown, tf = PrettyTables.tf_compact)
    if !isnothing(contab.id) && length(contab.id) > 0
        print(io, "  ID: ")
        for (k,v) in contab.id
            print(io, "$k => $v; ")
        end
    end
end


text:
struct UnionArray{T, N, P} <: Abstract.UnionArray{T, N}
    parent::P

    function UnionArray(parent::P) where {T, N, P <: AbstractArray{T, N}}
        # Make sure that `A.parent` is not a UnionArray/Vector
        P <: UnionArrayImpls && throw(MethodError(UnionArray, (parent,)))
        return new{T, N, P}(parent)
    end
end

const UnionArrayImpls = Union{UnionVector, UnionArray}

function Adapt.adapt_structure(to, A::UnionArray)
    B = Adapt.adapt(to, parent(A))
    if B isa UnionArray
        # `Adapt.adapt` might just call `reshape`, which then wraps the array
        # to `UnionArray`. In this case, there is no need to re-wrap it:
        return B
    else
        # Maybe just throw?
        return UnionArray(B)
    end
end

# TODO: handle `reshape` in Transducers.jl so that we only need one
# `parent(...)` call here:
executor_type(A::UnionArray) = executor_type(parent(parent(A)))

# A very minimal dummy array implementation just for implementing `default_reshape`:
struct DummyArray{N} <: AbstractArray{Any,N}
    dims::NTuple{N,Int}
end
Base.size(A::DummyArray) = A.dims

function default_reshape(A::UnionVector, dims::T) where {T}
    dummy = reshape(DummyArray(size(A)), dims)
    return @set dummy.parent = A
end

ua_reshape(A::UnionVector, dims) = UnionArray(default_reshape(A, dims))
ua_reshape(A::UnionArray, dims) = UnionArray(reshape(A.parent, dims))
ua_reshape(A::UnionVector, ::Tuple{Colon}) = A
function ua_reshape(A::UnionArray, ::Tuple{Colon})
    v = reshape(A.parent, :)
    if v isa UnionVector
        return v
    else
        return UnionArray(v)
    end
end

Base.reshape(A::UnionArrayImpls, dims::Dims) =
    ua_reshape(A, dims) :: UnionArrayImpls
Base.reshape(A::UnionArrayImpls, dims::Tuple{Vararg{Union{Int,Colon}}}) =
    ua_reshape(A, dims) :: UnionArrayImpls

# disambiguation:
Base.reshape(A::UnionArrayImpls, dims::Tuple{Colon}) =
    ua_reshape(A, dims) :: UnionArrayImpls

Base.parent(A::UnionArray) = A.parent
Base.size(A::UnionArray) = size(A.parent)
Base.getindex(A::UnionArray, I...) = A.parent[I...]
Base.setindex!(A::UnionArray, v, I...) = setindex!(A.parent, v, I...)

# Base.showarg(io::IO, A::UnionArray, toplevel) =
#     Base.showarg(io, A.parent, toplevel)

Abstract.UnionArray(
    undef::UndefInitializer,
    VectorType::Type{<:AbstractVector},
    ETS::ElTypeSpec,
    dims::Union{Integer,Tuple{Vararg{Integer}}},
) = Abstract.UnionArray(undef, VectorType, VectorType{UInt8}, ETS, dims)

Abstract.UnionArray(
    undef::UndefInitializer,
    DataVectorType::Type{<:AbstractVector},
    TypeTagVectorType::Type{<:AbstractVector{UInt8}},
    ETS::ElTypeSpec,
    n::Union{Integer,Tuple{Integer}},
) = UnionVector(undef, DataVectorType, TypeTagVectorType, ETS, first(n))

Abstract.UnionArray(
    undef::UndefInitializer,
    DataVectorType::Type{<:AbstractVector},
    TypeTagVectorType::Type{<:AbstractVector{UInt8}},
    ETS::ElTypeSpec,
    dims::Tuple{Integer,Integer,Vararg{Integer}},
) = reshape(UnionVector(undef, DataVectorType, TypeTagVectorType, ETS, prod(dims)), dims)

function Base.fill!(A::UnionArray, x)
    fill!(parent(parent(A)), x)
    return A
end