text:
#!/usr/bin/env julia

# Given some model parameters, synthesize and plot the channel maps and
# integrated spectrum

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--norad"
    help = "Use the image already here."
    action = :store_true
    "--config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskJockey.constants
using DiskJockey.image
using DiskJockey.hmodel
# using constants
# using image
# using model
using HDF5

import PyPlot.plt
using LaTeXStrings
import Images

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]

# First contour is at 3 sigma, and then contours go up (or down) in multiples of spacing
function get_levels(rms::Float64, vmax::Float64, spacing=3)
    levels = Float64[]

    val = 3 * rms
    while (val < vmax)
        append!(levels, [val])
        val += rms * spacing
    end

    return levels
end

# function plot_beam(ax, BMAJ, BMIN, xy=(1,-1))
#     BMAJ = 3600. * header["BMAJ"] # [arcsec]
#     BMIN = 3600. * header["BMIN"] # [arcsec]
#     BPA =  header["BPA"] # degrees East of North
#     # from matplotlib.patches import Ellipse
#     ax[:add_artist](PyPlot.matplotlib[:patches][:Ellipse](xy=xy, width=BMIN, height=BMAJ, angle=BPA, facecolor="0.8", linewidth=0.2))
# end

# Plot the raw channel maps directly from RADMC
function plot_chmaps(img::image.RawImage)

    (im_ny, im_nx) = size(img.data)[1:2] # y and x dimensions of the image

    # CO 2-1 rest frame
    # lam0 = cc/230.538e9 * 1e4 # [microns]
    # nlam = length(img.lams)

    # convert wavelengths to velocities
    # vels = c_kms * (img.lams .- lam0)/lam0

    fig, ax = plt.subplots(nrows=2, ncols=12, figsize=(12, 2.8))

    ext = (1, im_nx, 1, im_ny) # Python array convention for Matplotlib

    for row=1:2
        for col=1:12
            iframe = col + (row - 1) * 12

            if iframe > nlam
                # Stop if we run out of channels
                break
            end

            frame = img.data[:,:,iframe]
            frame += 1e-99 #Add a tiny bit so that we don't have log10(0)
            max = maximum(log10(frame))
            ax[row, col][:imshow](log10(frame), vmin=max - 6, vmax=max, interpolation="none", origin="lower", cmap=plt.get_cmap("PuBu"), extent=ext)
            levels = linspace(max - 0.8, max, 5)
            ax[row, col][:contour](log10(frame), origin="lower", colors="k", levels=levels, linestyles="solid", linewidths=0.2, extent=ext)

            if col != 1 || row != 2
                ax[row, col][:xaxis][:set_ticklabels]([])
                ax[row, col][:yaxis][:set_ticklabels]([])
            else
                ax[row, col][:set_xlabel](L"$X$")
                ax[row, col][:set_ylabel](L"$Y$")
            end

            ax[row, col][:annotate](@sprintf("%.1f", vels[iframe]), (0.1, 0.8), xycoords="axes fraction", size=8)

        end
    end

    fig[:subplots_adjust](hspace=0.01, wspace=0.05, top=0.9, bottom=0.1, left=0.05, right=0.95)

    plt.savefig("channel_maps_raw_image.png")

end

# Plot the channel maps using sky convention
function plot_chmaps(img::image.SkyImage, fname="channel_maps_sky.png")

    (im_ny, im_nx) = size(img.data)[1:2] # y and x dimensions of the image

    # Image needs to be flipped along RA dimension
    ext = (img.ra[end], img.ra[1], img.dec[1], img.dec[end])

    # Figure out how many plots we'll have.
    ncols = 8
    nrows = iceil(nlam/ncols)

    fig, ax = plt.subplots(nrows=nrows, ncols=ncols, figsize=(12, 1.5 * nrows))

    for row=1:nrows
        for col=1:ncols
            iframe = col + (row - 1) * ncols

            if col != 1 || row != nrows
                ax[row, col][:xaxis][:set_ticklabels]([])
                ax[row, col][:yaxis][:set_ticklabels]([])
            else
                ax[row, col][:set_xlabel](L"$\Delta \alpha$ ('')")
                ax[row, col][:set_ylabel](L"$\Delta \delta$ ('')")
            end

            if iframe > nlam
                # Stop if we run out of channels
                # Plot a blank square
                ax[row, col][:imshow](zeros((im_ny, im_nx)), cmap=plt.get_cmap("PuBu"), vmin=0, vmax=20, extent=ext, origin="lower")

            else
                #Flip the frame for Sky convention
                frame = fliplr(img.data[:,:,iframe])
                frame += 1e-99 #Add a tiny bit so that we don't have log10(0)
                lframe = log10(frame)
                max = maximum(lframe)
                ix,iy = ind2sub(size(lframe), indmax(lframe))

                # ax[row, col][:imshow](lframe, extent=ext, interpolation="none", vmin=max - 6, vmax=max, origin="lower", cmap=plt.get_cmap("PuBu"))
                ax[row, col][:imshow](frame, extent=ext, interpolation="none", origin="lower", cmap=plt.get_cmap("PuBu"), norm=norm)

                # println("frame ", iframe)
                # println("levels ", levels)
                # println("extrema ", extrema(frame))
                if length(levels) > 0
                    ax[row, col][:contour](frame, origin="lower", colors="k", levels=levels, extent=ext, linestyles="solid", linewidths=0.2)
                end

                ax[row, col][:annotate](@sprintf("%.1f", vels[iframe]), (0.1, 0.8), xycoords="axes fraction", size=8)
            end

        end
    end

    fig[:subplots_adjust](hspace=0.06, wspace=0.01, top=0.9, bottom=0.1, left=0.05, right=0.95)

    plt.savefig(fname)

end

# Plot the raw array for the Sky image
function plot_chmaps_data(img::image.SkyImage)

    (im_ny, im_nx) = size(img.data)[1:2] # y and x dimensions of the image

    ll = sin(img.ra .* arcsec)
    mm = sin(img.dec .* arcsec)
    ext = (ll[1], ll[end], mm[1], mm[end])

    nlam = length(img.lams)

    fig, ax = plt.subplots(nrows=2, ncols=12, figsize=(12, 2.8))

    for row=1:2
        for col=1:12
            iframe = col + (row - 1) * 12

            if iframe > nlam
                # Stop if we run out of channels
                break
            end

            frame = img.data[:,:,iframe]
            frame += 1e-99 #Add a tiny bit so that we don't have log10(0)
            max = maximum(log10(frame))
            ax[row, col][:imshow](log10(frame), extent=ext, vmin=max - 6, vmax=max, interpolation="none", origin="lower", cmap=plt.get_cmap("PuBu"))
            levels = linspace(max - 0.8, max, 5)
            ax[row, col][:contour](log10(frame), origin="lower", colors="k", levels=levels, extent=ext, linestyles="solid", linewidths=0.2)

            if col != 1 || row != 2
                ax[row, col][:xaxis][:set_ticklabels]([])
                ax[row, col][:yaxis][:set_ticklabels]([])
            else
                ax[row, col][:set_xlabel](L"$ll$")
                ax[row, col][:set_ylabel](L"$mm$")
            end

            ax[row, col][:annotate](@sprintf("%.1f", vels[iframe]), (0.1, 0.8), xycoords="axes fraction", size=8)

        end
    end

    fig[:subplots_adjust](hspace=0.01, wspace=0.05, top=0.9, bottom=0.1, left=0.05, right=0.95)

    plt.savefig("channel_maps_sky_raw.png")

end

# Plot the spatially-integrated spectrum
function plot_spectrum(img::image.SkyImage)

    fig = plt.figure()
    ax = fig[:add_subplot](111)

    spec = imToSpec(img)

    ax[:plot](vels, spec[:,2], ls="steps-mid")
    # ax[:plot](vels, reverse(spec[:,2]), ls="steps-mid")
    ax[:set_ylabel](L"$f_\nu$ [Jy]")
    ax[:set_xlabel](L"$v$ [km/s]")

    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt.savefig("spectrum.png")
end

pp = config["parameters"]
params = ["M_star", "r_c", "T_10m", "q_m", "T_10a", "q_a", "T_freeze", "X_freeze", "gamma", "h", "delta", "logM_gas", "ksi", "dpc", "incl", "PA", "vel", "mu_RA", "mu_DEC"]
nparam = length(params)
starting_param = Array(Float64, nparam)

for i=1:nparam
    starting_param[i] = pp[params[i]][1]
end

# Convert logM_gas to M_gas
starting_param[12] = 10^starting_param[12]

pars = Parameters(starting_param...)
#
# vel = pars.vel # [km/s]
# # RADMC conventions for inclination and PA
incl = pars.incl # [deg]
PA = pars.PA # [deg] Position angle runs counter clockwise
npix = config["npix"] # number of pixels
#
# # Doppler shift the dataset wavelength to rest-frame wavelength
# beta = vel/c_kms # relativistic Doppler formula
# shift_lams =  lams .* sqrt((1. - beta) / (1. + beta)) # [microns]
#
# grd = config["grid"]
# grid = Grid(grd["nr"], grd["ntheta"], grd["r_in"], grd["r_out"], true)
#
# # Comment out afterward
# # nchan = 251
# # global const vels = linspace(-10., 10, nchan) # [km/s]
# # # CO 2-1 rest frame
# # lam0 = cc/230.538e9 * 1e4 # [microns]
# # # convert velocities to wavelengths
# # shift_lams = lam0 * (vels/c_kms + 1)

if !parsed_args["norad"]
    run(`radmc3d image incl $incl posang $PA npix $npix loadlambda`)
end

im = imread()

# plot_chmaps(im)

skim = imToSky(im, pars.dpc)

# Do the velocity conversion here
global nlam = length(skim.lams)
# convert wavelengths to velocities
global vels = c_kms * (skim.lams .- lam0)/lam0

# get the colorscale corresponding to the velocities
vmin, vmax = extrema(skim.data)
vvmax = maxabs(skim.data)

ldata = log10(skim.data + 1e-99)
vlmax = maximum(ldata)

# if in(config, "beam")
beam = config["beam"]
rms = beam["rms"] # Jy/beam
BMAJ = beam["BMAJ"]/2 # semi-major axis [arcsec]
BMIN = beam["BMIN"]/2 # semi-minor axis [arcsec]
BAVG = (BMAJ + BMIN)/2
BPA = beam["BPA"] # position angle East of North [degrees]

println("Beam sigma ", BAVG, " [arcsec]")

arcsec_ster = (4.25e10)
# Convert beam from arcsec^2 to Steradians
rms = rms/(pi * BMAJ * BMIN) * arcsec_ster
# println("RMS: ", rms, " vmax: ", vmax)
global levels = get_levels(rms, vmax)

norm = PyPlot.matplotlib[:colors][:Normalize](vlmax - 8, vlmax)

# println("Plotting hires maps")
# plot_chmaps(skim)

println("bluring maps")
sk_blur = blur(skim, [BAVG, BAVG])

# Now redo all this for the blurred
vmin, vmax = extrema(sk_blur.data)
global levels = get_levels(rms, vmax)
norm = PyPlot.matplotlib[:colors][:Normalize](0, vmax)

println("Plotting blured maps")
plot_chmaps(sk_blur, "channel_maps_blur.png")

plot_spectrum(skim)

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
# Copyright (c) 2019 Mathieu Besançon, Oscar Dowson, and contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestCHighs

using HiGHS
using Test

function small_test_model()
    cc = [1.0, -2.0]
    cl = [0.0, 0.0]
    cu = [10.0, 10.0]
    ru = [2.0, 1.0]
    rl = [0.0, 0.0]
    astart = [0, 2, 4]
    aindex = [0, 1, 0, 1]
    avalue = [1.0, 2.0, 1.0, 3.0]
    return (cc, cl, cu, rl, ru, astart, aindex, avalue)
end

function test_Direct_C_call()
    (colcost, collower, colupper, rowlower, rowupper, astart, aindex, avalue) =
        small_test_model()
    n_col = convert(Cint, size(colcost, 1))
    n_row = convert(Cint, size(rowlower, 1))
    n_nz = convert(Cint, size(aindex, 1))

    colcost = convert(Array{Cdouble}, colcost)
    collower = convert(Array{Cdouble}, collower)
    colupper = convert(Array{Cdouble}, colupper)

    rowlower = convert(Array{Cdouble}, rowlower)
    rowupper = convert(Array{Cdouble}, rowupper)
    matstart = convert(Array{Cint}, astart)
    matindex = convert(Array{Cint}, aindex)
    matvalue = convert(Array{Cdouble}, avalue)

    colvalue, coldual =
        (Array{Cdouble,1}(undef, n_col), Array{Cdouble,1}(undef, n_col))
    rowvalue, rowdual =
        (Array{Cdouble,1}(undef, n_row), Array{Cdouble,1}(undef, n_row))
    colbasisstatus, rowbasisstatus =
        (Array{Cint,1}(undef, n_col), Array{Cint,1}(undef, n_row))

    modelstatus = Ref{Cint}(42)
    status = HiGHS.Highs_lpCall(
        n_col,
        n_row,
        n_nz,
        HiGHS.kHighsMatrixFormatColwise,
        HiGHS.kHighsObjSenseMaximize,
        0.0,
        colcost,
        collower,
        colupper,
        rowlower,
        rowupper,
        matstart,
        matindex,
        matvalue,
        colvalue,
        coldual,
        rowvalue,
        rowdual,
        colbasisstatus,
        rowbasisstatus,
        modelstatus,
    )
    @test status == 0
    @test modelstatus[] == HiGHS.kHighsModelStatusOptimal
end

function test_create()
    ptr = Highs_create()
    @test ptr != C_NULL

    objective = Highs_getObjectiveValue(ptr)
    @test objective == 0.0

    # model status, the second parameter is
    # bool scaled_model
    modelstatus = Highs_getModelStatus(ptr)
    @test modelstatus[] == 0 # uninitialized LP

    return Highs_destroy(ptr)
end

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith(string(name), "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
end

end

TestCHighs.runtests()

ext: jl


text:
using Revise
using LinearAlgebra, StaticArrays, GeometryPrimitives, ArrayInterface, ChainRules, ForwardDiff, Zygote, FiniteDifferences, BenchmarkTools, LoopVectorization, Tullio
using Test, ChainRulesTestUtils
using Zygote: @ignore, dropgrad, Context, _pullback, @adjoint
using ForwardDiff: Dual, Tag
using ChainRules: @non_differentiable

ChainRulesCore.rrule(T::Type{<:SArray}, xs::Number...) = ( T(xs...), dv -> (nothing, dv...) )
ChainRulesCore.rrule(T::Type{<:SArray}, x::AbstractArray) = ( T(x), dv -> (nothing, dv) )
ChainRulesCore.rrule(T::Type{<:SMatrix}, xs::Number...) = ( T(xs...), dv -> (nothing, dv...) )
ChainRulesCore.rrule(T::Type{<:SMatrix}, x::AbstractMatrix) = ( T(x), dv -> (nothing, dv) )
ChainRulesCore.rrule(T::Type{<:SVector}, xs::Number...) = ( T(xs...), dv -> (nothing, dv...) )
ChainRulesCore.rrule(T::Type{<:SVector}, x::AbstractVector) = ( T(x), dv -> (nothing, dv) )

function sum2(op,arr)
    return sum(op,arr)
end
function sum2adj( Δ, op, arr )
    n = length(arr)
    g = x->Δ*Zygote.gradient(op,x)[1]
    return ( nothing, map(g,arr))
end
@adjoint function sum2(op,arr)
    return sum2(op,arr),Δ->sum2adj(Δ,op,arr)
end

# AD rules for fast norms of types SVector{2,T} and SVector{2,3}

function _norm2_back_SV2r(x::SVector{2,T}, y, Δy) where T<:Real
    ∂x = Vector{T}(undef,2)
    ∂x .= x .* (real(Δy) * pinv(y))
    return reinterpret(SVector{2,T},∂x)[1]
end

function _norm2_back_SV3r(x::SVector{3,T}, y, Δy) where T<:Real
    ∂x = Vector{T}(undef,3)
    ∂x .= x .* (real(Δy) * pinv(y))
    return reinterpret(SVector{3,T},∂x)[1]
end

function _norm2_back_SV2c(x::SVector{2,T}, y, Δy) where T<:Complex
    ∂x = Vector{T}(undef,2)
    ∂x .= conj.(x) .* (real(Δy) * pinv(y))
    return reinterpret(SVector{2,T},∂x)[1]
end

function _norm2_back_SV3c(x::SVector{3,T}, y, Δy) where T<:Complex
    ∂x = Vector{T}(undef,3)
    ∂x .= conj.(x) .* (real(Δy) * pinv(y))
    return reinterpret(SVector{3,T},∂x)[1]
end

function ChainRulesCore.rrule(::typeof(norm), x::SVector{3,T}) where T<:Real
	y = LinearAlgebra.norm(x)
	function norm_pb(Δy)
		∂x = Thunk() do
			_norm2_back_SV3r(x, y, Δy)
		end
		return ( NO_FIELDS, ∂x )
	end
	norm_pb(::Zero) = (NO_FIELDS, Zero())
    return y, norm_pb
end

function ChainRulesCore.rrule(::typeof(norm), x::SVector{2,T}) where T<:Real
	y = LinearAlgebra.norm(x)
	function norm_pb(Δy)
		∂x = Thunk() do
			_norm2_back_SV2r(x, y, Δy)
		end
		return ( NO_FIELDS, ∂x )
	end
	norm_pb(::Zero) = (NO_FIELDS, Zero())
    return y, norm_pb
end

function ChainRulesCore.rrule(::typeof(norm), x::SVector{3,T}) where T<:Complex
	y = LinearAlgebra.norm(x)
	function norm_pb(Δy)
		∂x = Thunk() do
			_norm2_back_SV3c(x, y, Δy)
		end
		return ( NO_FIELDS, ∂x )
	end
	norm_pb(::Zero) = (NO_FIELDS, Zero())
    return y, norm_pb
end

function ChainRulesCore.rrule(::typeof(norm), x::SVector{2,T}) where T<:Complex
	y = LinearAlgebra.norm(x)
	function norm_pb(Δy)
		∂x = Thunk() do
			_norm2_back_SV2c(x, y, Δy)
		end
		return ( NO_FIELDS, ∂x )
	end
	norm_pb(::Zero) = (NO_FIELDS, Zero())
    return y, norm_pb
end

# #### Polygon fn defs
# function v_pgn(wₜₒₚ::T1,t_core::T1,θ::T1,edge_gap::T1,n_core::T1,n_subs::T1,Δx::T2,Δy::T2) where {T1<:Real,T2<:Real}
# 	t_subs = (Δy -t_core - edge_gap )/2.
# 	c_subs_y = -Δy/2. + edge_gap/2. + t_subs/2.
# 	ε_core = ε_tensor(n_core)
# 	ε_subs = ε_tensor(n_subs)
# 	wt_half = wₜₒₚ / 2
# 	wb_half = wt_half + ( t_core * tan(θ) )
# 	@show tc_half = t_core / 2
# 	verts =   [   	wt_half     -wt_half     -wb_half    wb_half
# 					tc_half     tc_half    -tc_half      -tc_half    ]'
#     return SMatrix{4,2}(verts)
# end
# v_pgn(p::AbstractVector) = v_pgn(p[1],p[2],p[3],p[6],p[4],p[5],Δx,Δy)
# f_onbnd_poly(pin,abs∆xe) =  (p = Zygote.@ignore pin; onbnd = Zygote.@ignore ( abs∆xe .≤ Base.rtoldefault(Float64) * maximum(abs.((-)(bounds(p)...))) ); onbnd)
# f_isout_poly(p,∆xe) =  (isout = Zygote.@ignore ( (∆xe.>0) .| f_onbnd_poly(p,abs.(∆xe))); isout)
# function l_bnds_poly(v::SMatrix{K,2,T})::SVector{2,T} where {K,T<:Real}
# 	@tullio (min) res[b] :=  (v)[a,b]
# end
# function u_bnds_poly(v::SMatrix{K,2,T})::SVector{2,T} where {K,T<:Real}
# 	@tullio (max) res[b] :=  (v)[a,b]
# end
# bounds2(s::Polygon) = l_bnds_poly(s.v), u_bnds_poly(s.v)
# function _Δxe_poly(x::SVector{2,T},v::SMatrix{K,2,T},n::SMatrix{K,2,T})::SVector{K,T}  where {K,T<:Real}
# 	@tullio out[k] := n[k,a] * ( x[a] - v[k,a] ) # edge line eqns for a K-point Polygon{K} `s`
# end
# function _Δxe_poly2(x::SVector{2},v::SMatrix{K,2},n::SMatrix{K,2})  where {K}   #
# 	@tullio out[k] := n[k,a] * ( x[a] - v[k,a] ) # edge line eqns for a K-point Polygon{K} `s`
# end
# Δxe_poly(x::AbstractVector{<:Real},s::Polygon) = _Δxe_poly2(x,s.v,s.n)
# function f_onbnd2(Δxe::SVector{K},s::Polygon{K})::SVector{K} where K
# 	map(x->abs(x)≤rmin,Δxe)
# end
# function f_onbnd2(x::SVector{2},s::Polygon{K})::SVector{K} where K
# 	map(x->abs(x)≤rmin,Δxe_poly(x,s))
# end
# function _Δx_poly(x::SVector{2},v::SMatrix{K,2},n::SMatrix{K,2},kmax::Int)::SVector{2}  where {K} #,T<:Real}
# 	# @tullio Δx[i] := n[$kmax,a] * ( v[$kmax,a] - x[a] * n[$kmax,i]	# works but gradient doesn't vectorize
# 	( @tullio Δx_factor := n[$kmax,a] * ( v[$kmax,a] - x[a] ) )  * SVector(n[kmax,1],n[kmax,2])
# end
# function _Δx_poly3(x::SVector{2},v::SMatrix{K,2},n::SMatrix{K,2},kmax::Int)  where {K} #,T<:Real}
# 	@tullio Δx[i] := n[$kmax,a] * ( v[$kmax,a] - x[a] ) * n[$kmax,i]	# works but gradient doesn't vectorize
# 	# ( @tullio Δx_factor := n[$kmax,a] * ( v[$kmax,a] - x[a] ) )  * SVector(n[kmax,1],n[kmax,2])
# end
# Δx_poly(x::AbstractVector{<:Real},s::Polygon)::SVector{2} = _Δx_poly(x,s.v,s.n,argmax(_Δxe_poly(x,s.v,s.n)))
#
#
# function f_Δxe(n::AbstractArray{<:Real},x::AbstractVector{<:Real},v::AbstractArray{<:Real})::AbstractVector{<:Real}
# 	# edge line eqns for a K-point Polygon{K} `s`:
# 	# ∆xe::SVector{K} = sum(s.n .* (x' .- s.v), dims=Val(2))[:,1]
# 	@tullio out[k] := n[k,a] * ( x[a] - v[k,a] )
# end
# function f_sz(n::AbstractArray,x::AbstractVector,v::AbstractArray)::AbstractVector
# 	# Determine if x is outside of edges, inclusive.
# 	# sz = abs.((-)(bounds(s)...))  # SVector{2}
# 	# onbnd = abs∆xe .≤ Base.rtoldefault(Float64) * max(sz.data...)  # SVector{K}
# 	@tullio sz (-) := n[k,a] * ( x[a] - v[k,a] )
# end
# @non_differentiable f_onbnd_poly(::Any,::Any)
# @non_differentiable f_isout_poly(::Any,::Any)
#
# function _∆x_poly(x::SVector{2},v::SMatrix{K,2},n::SMatrix{K,2},kmax::Int)::SVector{2}  where {K} #,T<:Real}
# 	# @tullio Δx[i] := n[$kmax,a] * ( v[$kmax,a] - x[a] * n[$kmax,i]	# works but gradient doesn't vectorize
# 	( @tullio Δx_factor := n[$kmax,a] * ( v[$kmax,a] - x[a] ) )  * SVector(n[kmax,1],n[kmax,2])
# end
# ∆x_poly(x::AbstractVector{<:Real},s::Polygon)::SVector{2} = _∆x_poly(x,s.v,s.n,argmax(_∆xe_poly(x,s.v,s.n)))
#
# function surfpt_nearby2(x::SVector{2,T1}, s::Polygon{K,K2,D,T2})::Tuple{SVector{2, <:Real}, SVector{2, <:Real}} where {K,K2,D,T1<:Real,T2<:Real}
#     #∆xe::SVector{K,T} = @inbounds sum(s.n .* (x' .- s.v), dims=2)[1:K,1]  # Calculate the signed distances from x to edge lines.
# 	∆xe = Δxe_poly(x,s) ##∆xe_poly(x,s)
# 	# abs∆xe::SVector{K,T} = abs.(∆xe)
#     onbnd = Zygote.@ignore( f_onbnd_poly(s, abs.(∆xe) ) )# abs∆xe .≤ Base.rtoldefault(Float64) * max(sz.data...)  # SVector{K}
#     # isout::BitVector = Zygote.@ignore(  )#(∆xe.>0) .| onbnd  # SVector{K}
#     cout = Zygote.@ignore( count(f_isout_poly(s,∆xe)) )
#     if cout == 2  # x is outside two edges
# 		# println("case 1")
# 		∆xv = ones(4)*x' - s.v #x' .- s.v::SMatrix{K,2,T,K2}
#         l∆xv = map(sqrt,mapreduce(abs2,+,∆xv,dims=2))# sum(abs2,Δxv;dims=2) # hypot.(∆xv[1:K,1], ∆xv[1:K,2])
#         imin = 2 #argmin(l∆xv)
#         surf = @inbounds SVector(s.v[imin,1],s.v[imin,2])
#         imin₋₁ =1 # mod1(imin-1,K)
#         if onbnd[imin] && onbnd[imin₋₁]  # x is very close to vertex imin
# 			# println("case 1.1")
# 			nout = SVector( normalize( [ s.n[imin,1]+s.n[imin₋₁,1], s.v[imin,2]+s.n[imin₋₁,2] ]) )
# 			# nout = reinterpret(SVector{2,T}, normalize( [ s.n[imin,1]+s.n[imin₋₁,1],s.v[imin,2]+s.n[imin₋₁,2] ]))  #  s.n[imin,:] + s.n[imin₋₁,:]
#         else
# 			# println("case 1.2")
# 			nout = SVector( normalize( [ x[1] - s.v[imin,1] , x[2] - s.v[imin,2] ] ) )
#             # nout = reinterpret(SVector{2,T}, normalize( [ x[1] - s.v[imin,1] , x[2] - s.v[imin,2] ] ) )[1]
#         end
#         # nout = normalize(nout)
#     else  # cout ≤ 1 or cout ≥ 3
# 		# println("case 2")
#         imax::Int = Zygote.@ignore(argmax(∆xe))
#         # vmax = SVector{2,T}(s.v[imax,1], s.v[imax,2])
# 		# nout  = SVector{2,T}(s.n[imax,1], s.n[imax,2])
#         # ∆x = (nout⋅(vmax-x)) .* nout
#         # surf = x + ∆x
# 		surf = x + _∆x_poly(x,s.v,s.n,imax) # ∆x = (nmax⋅(vmax-x)) .* nmax
#         nout = @inbounds SVector(s.n[imax,1],s.n[imax,2] )
# 		# vmax = SVector( 1.009770432073105, -0.35 )  #::SVector{2,T} = SVector{2,T}(s.v[imax,1], s.v[imax,2])
# 		# nout = SVector( 0.9749279121818236, 0.2225209339563144 )# = SVector{2,T}(s.n[imax,1], s.n[imax,2])
#         # ∆x = SVector( 0.006711245692134241, 0.0015317980342586456 ) #::SVector{2,T} = (nout⋅(vmax-x)) .* nout
#         # surf = SVector( 0.9067112456921342, 0.10153179803425866 ) # = x + ∆x
#         # nout = nmax
#     end
#     return surf, nout::SVector
# end
#
#
#
# f_onbnd(bin,absdin) =  (b = Zygote.@ignore bin; absd = Zygote.@ignore absdin; abs.(b.r.-absd) .≤ Base.rtoldefault(Float64) .* b.r)  # basically b.r .≈ absd but faster
#
# f_isout(b,absd) =  (isout = Zygote.@ignore ((b.r.<absd) .| f_onbnd(b,absd) ); isout)
# # isout(bin) =  (b = Zygote.@ignore bin; (b.r.<absd) .| (abs.(b.r.-absd) .≤ Base.rtoldefault(Float64) .* b.r))
# f_signs(d) =  (signs = Zygote.@ignore (copysign.(1.0,d)'); signs)
#
#
# function surfpt_nearby2(x::AbstractVector{<:Real}, b::GeometryPrimitives.Box{2,4,D,T}) where {D,T<:Real}
#     ax = inv(b.p)
#     n0 = b.p ./  [ sqrt(b.p[1,1]^2 + b.p[1,2]^2) sqrt(b.p[2,1]^2 + b.p[2,2]^2)  ]
#     d = Array(b.p * (x - b.c))
#     cosθ = diag(n0*ax)
#
#     n = n0 .* copysign.(1.0,d)' #f_signs(d)
#     absd = abs.(d)
#     ∆ = (b.r .- absd) .* cosθ
#     # onbnd = abs.(dropgrad(b.r).-dropgrad(absd)) .≤ Base.rtoldefault(Float64) .* dropgrad(b.r)  # basically b.r .≈ absd but faster
#     # isout = (dropgrad(b.r).<dropgrad(absd)) .| dropgrad(onbnd)
#     onbnd = f_onbnd(b,absd)
#     isout = f_isout(b,absd)
#     projbnd =  all(.!isout .| onbnd)
#     # onbnd_float =  reinterpret(T,onbnd) #Float64.(onbnd)
#     # isout_float =  reinterpret(T,isout) #Float64.(isout)
#     # if ( ( abs(b.r[1]-absd[1]) ≤ Base.rtoldefault(Float64) * b.r[1] ) | (b.r[1]<absd[1]) ) | ( ( abs(b.r[2]-absd[2]) ≤ Base.rtoldefault(Float64) * b.r[2] ) | ( b.r[2] < absd[2] ) )
#     if count(isout) == 0
#         # = not(isout[1]) &  not(isout[2]) case, point is inside box
#         l∆x, i = findmin(∆)  # find closest face
#         nout = n[i,:]
#         ∆x = l∆x * nout
#     else
#         ∆x = n' * (∆ .* isout)
#         # = isout[1] | isout[2] case 1: at least one dimension outside or on boundry
#         if all(.!isout .| onbnd)
#             nout0 = n' * onbnd
#         else
#             nout0 = -∆x
#         end
#         nout = nout0 / norm(nout0)
#     end
#
#     return SVector{2}(x+∆x), SVector{2}(nout)
# end

@non_differentiable KDTree(::Any)
@non_differentiable copysign(::Any,::Any)
ChainRules.refresh_rules()
Zygote.refresh()

const εᵥ = I
function ε_tensor(n::T)::SMatrix{3,3,T,9} where T<:Real
    n² = n^2
    ε =     [	n²      0. 	    0.
                0. 	    n² 	    0.
                0. 	    0. 	    n²  ]
end
function ridge_wg(wₜₒₚ::T1,t_core::T1,θ::T1,edge_gap::T1,n_core::T1,n_subs::T1,Δx::T2,Δy::T2)::Vector{<:GeometryPrimitives.Shape} where {T1<:Real,T2<:Real}
    t_subs = (Δy -t_core - edge_gap )/2.
    c_subs_y = -Δy/2. + edge_gap/2. + t_subs/2.
    ε_core = ε_tensor(n_core)
    ε_subs = ε_tensor(n_subs)
    wt_half = wₜₒₚ / 2
    wb_half = wt_half + ( t_core * tan(θ) )
    tc_half = t_core / 2
    verts =   [   wt_half     -wt_half     -wb_half    wb_half
                    tc_half     tc_half    -tc_half      -tc_half    ]'
    core = GeometryPrimitives.Polygon(                                                                # Instantiate 2D polygon, here a trapazoid
                    SMatrix{4,2}(verts),                                                                    # v: polygon vertices in counter-clockwise order
                    ε_core,                                                                            # data: any type, data associated with box shape
                )
    ax = [      1.     0.
                0.     1.      ]
    b_subs = GeometryPrimitives.Box(                                                        # Instantiate N-D box, here N=2 (rectangle)
                    [0. , c_subs_y],                    # c: center
                    [Δx - edge_gap, t_subs ],        # r: "radii" (half span of each axis)
                    ax,                                    # axes: box axes
                    ε_subs,                                                # data: any type, data associated with box shape
                )
    return [core,b_subs]
end
p = [
    1.7,                #   top ridge width         `w_top`         [μm]
    0.7,                #   ridge thickness         `t_core`        [μm]
    π / 14.0,           #   ridge sidewall angle    `θ`             [radian]
    2.4,                #   core index              `n_core`        [1]
    1.4,                #   substrate index         `n_subs`        [1]
    0.5,                #   vacuum gap at boundaries `edge_gap`     [μm]
]
Δx,Δy,Δz,Nx,Ny,Nz = 6., 4., 1., 128, 128, 1
rwg(p) = ridge_wg(p[1],p[2],p[3],p[6],p[4],p[5],Δx,Δy)
shapes = ridge_wg(p[1],p[2],p[3],p[6],p[4],p[5],Δx,Δy)
s = shapes[1]

function τ_trans(ε::AbstractMatrix{T}) where T<:Real
    return @inbounds SMatrix{3,3,T,9}(
        -1/ε[1,1],      ε[2,1]/ε[1,1],                  ε[3,1]/ε[1,1],
        ε[1,2]/ε[1,1],  ε[2,2] - ε[2,1]*ε[1,2]/ε[1,1],  ε[3,2] - ε[3,1]*ε[1,2]/ε[1,1],
        ε[1,3]/ε[1,1],  ε[2,3] - ε[2,1]*ε[1,3]/ε[1,1],  ε[3,3] - ε[3,1]*ε[1,3]/ε[1,1]
    )
end

function τ⁻¹_trans(τ::AbstractMatrix{T}) where T<:Real
    return @inbounds SMatrix{3,3,T,9}(
        -1/τ[1,1],          -τ[2,1]/τ[1,1],                 -τ[3,1]/τ[1,1],
        -τ[1,2]/τ[1,1],     τ[2,2] - τ[2,1]*τ[1,2]/τ[1,1],  τ[3,2] - τ[3,1]*τ[1,2]/τ[1,1],
        -τ[1,3]/τ[1,1],     τ[2,3] - τ[2,1]*τ[1,3]/τ[1,1],  τ[3,3]- τ[3,1]*τ[1,3]/τ[1,1]
    )
end

import Base: findfirst
function Base.findfirst(p::SVector{N}, s::Vector{S},backup::Shape{N})::Shape{N} where {N,S<:Shape{N}}
    for i in eachindex(s)
        b = bounds(s[i])
        if all(b[1] .< p .< b[2]) && p ∈ s[i]  # check if p is within bounding box is faster
            return s[i]
        end
    end
    return backup
end

function Base.findfirst(p::SVector{N}, kd::KDTree{N},backup::Shape{N})::Shape{N} where {N}
    if isempty(kd.s)
        if p[kd.ix] ≤ kd.x
            return findfirst(p, kd.left, backup)
        else
            return findfirst(p, kd.right, backup)
        end
    else
        return findfirst(p, kd.s, backup)
    end
end

Base.findfirst(p::Vector{<:Real}, kd::KDTree{N}, backup) where {N} = findfirst(SVector{N}(p), kd, backup)
Base.findfirst(p::Vector{<:Real}, s::Vector{<:Shape{N}}, backup) where {N} = findfirst(SVector{N}(p), s, backup)

function findfirst2(p::SVector{N}, s::Vector{S},bu) where {N,S<:Shape{N}}
    for i in eachindex(s)
        b = @inbounds bounds(s[i])
        if @inbounds all(b[1] .< p .< b[2]) && p ∈ s[i]  # check if p is within bounding box is faster
            return @inbounds s[i]
        end
    end
    return bu
end

function findfirst2(p::SVector{N}, kd::KDTree{N},bu) where {N}
    if isempty(kd.s)
        if p[kd.ix] ≤ kd.x
            return findfirst2(p, kd.left,bu)
        else
            return findfirst2(p, kd.right,bu)
        end
    else
        return findfirst2(p, kd.s,bu)
    end
end

function avg_param2(param1, param2, n12, rvol1)
    n = normalize(n12) #n12 / norm(n12) #sqrt(sum2(abs2,n12))
    # Pick a vector that is not along n.
    h = any(iszero.(n)) ? n × normalize(iszero.(n)) :  n × SVector(1., 0. , 0.)
	v = n × h
    # Create a local Cartesian coordinate system.
    S = [n h v]  # unitary
    τ1 = τ_trans(transpose(S) * param1 * S)  # express param1 in S coordinates, and apply τ transform
    τ2 = τ_trans(transpose(S) * param2 * S)  # express param2 in S coordinates, and apply τ transform
    τavg = τ1 .* rvol1 + τ2 .* (1-rvol1)  # volume-weighted average
    return S * τ⁻¹_trans(τavg) * transpose(S)  # apply τ⁻¹ and transform back to global coordinates
end

function make_ei(shapes::Vector{<:Shape};Δx::T,Δy::T,Δz::T,Nx::Int,Ny::Int,Nz::Int,
	 	δx::T=Δx/Nx, δy::T=Δy/Ny, δz::T=Δz/Nz, x=( ( δx .* (0:(Nx-1))) .- Δx/2. ),
		y=( ( δy .* (0:(Ny-1))) .- Δy/2. ), z=( ( δz .* (0:(Nz-1))) .- Δz/2. ) ) where T<:Real
    tree = Zygote.@ignore(KDTree(shapes))
    eibuf = Zygote.Buffer(randn(T,2),3,3,Nx,Ny,Nz)
	# map((a,b)->εₛ(shapes,a,b;tree,δx,δy),x,y)
	for i ∈ eachindex(x), j ∈ eachindex(y)
	# 	# eps = εₛ(shapes,Zygote.dropgrad(tree),Zygote.dropgrad(g.x[i]),Zygote.dropgrad(g.y[j]),Zygote.dropgrad(g.δx),Zygote.dropgrad(g.δy))
		eps = εₛ(shapes,dropgrad(x)[i],dropgrad(y)[j];tree=dropgrad(tree),δx=dropgrad(δx),δy=dropgrad(δy))
		epsi = inv(eps) # inv( (eps' + eps) / 2) # Hermitian(inv(eps))  # inv(Hermitian(eps)) #   # inv(eps)
        eibuf[:,:,i,j,1] = epsi # inv(εₛ(shapes,x[i],y[j];tree,δx,δy)) # epsi #(epsi' + epsi) / 2
    end
    return real(copy(eibuf)) #HybridArray{Tuple{3,3,Dynamic(),Dynamic(),Dynamic()},T,5,5,Array{T,5}}( real(copy(eibuf)) )
end

function make_ei2(shapes::Vector{<:Shape},x,y)
    tree = KDTree(shapes)
    # eibuf = Zygote.Buffer(x,3,3,Nx,Ny,Nz)
	eibuf = Zygote.Buffer(zeros(SMatrix{3,3,Float64,9},2),Nx,Ny,Nz)
	# map((a,b)->εₛ(shapes,a,b;tree,δx,δy),x,y)
	for i ∈ eachindex(x), j ∈ eachindex(y)
	# 	# eps = εₛ(shapes,Zygote.dropgrad(tree),Zygote.dropgrad(g.x[i]),Zygote.dropgrad(g.y[j]),Zygote.dropgrad(g.δx),Zygote.dropgrad(g.δy))
		# eps = εₛ(shapes,x[i],y[j];tree,δx,δy)
		# epsi = inv(eps) # inv( (eps' + eps) / 2) # Hermitian(inv(eps))  # inv(Hermitian(eps)) #   # inv(eps)
        eibuf[i,j,1] = inv(εₛ(shapes,x[i],y[j];tree,δx,δy)::SMatrix{3,3,Float64,9}) # epsi #(epsi' + epsi) / 2
    end
    return real(copy(eibuf)) #HybridArray{Tuple{3,3,Dynamic(),Dynamic(),Dynamic()},T,5,5,Array{T,5}}( real(copy(eibuf)) )
end

function make_ei3(eibuf,shapes::Vector{<:Shape},x,y)
    tree = KDTree(shapes)
    # eibuf = Zygote.Buffer(x,3,3,Nx,Ny,Nz)
	# eibuf = Zygote.Buffer(zeros(SMatrix{3,3,Float64,9},2),Nx,Ny,Nz)
	# map((a,b)->inv(εₛ(shapes,a,b;tree,δx,δy)),x,y)
	for i ∈ eachindex(x), j ∈ eachindex(y)
	# 	# eps = εₛ(shapes,Zygote.dropgrad(tree),Zygote.dropgrad(g.x[i]),Zygote.dropgrad(g.y[j]),Zygote.dropgrad(g.δx),Zygote.dropgrad(g.δy))
		# eps = εₛ(shapes,x[i],y[j];tree,δx,δy)
		# epsi = inv(eps) # inv( (eps' + eps) / 2) # Hermitian(inv(eps))  # inv(Hermitian(eps)) #   # inv(eps)
        @inbounds eibuf[i,j,1] = inv(εₛ(shapes,x[i],y[j];tree,δx,δy)::SMatrix{3,3}) # epsi #(epsi' + epsi) / 2
    end
    return eibuf #HybridArray{Tuple{3,3,Dynamic(),Dynamic(),Dynamic()},T,5,5,Array{T,5}}( real(copy(eibuf)) )
end

make_ei2_fwd(shapes,x,y) = Zygote.forwarddiff(shapes) do shapes
	make_ei2(shapes,x,y)
end

## Polygon testing

function foo1(p)
	s = ridge_wg(p[1],p[2],p[3],p[6],p[4],p[5],6.,4.)[1]
	s.n[2,1] * s.v[3,2]
end
function foo2(p)
	s = ridge_wg(p[1],p[2],p[3],p[6],p[4],p[5],6.,4.)[1]
	xy = SVector(0.9,0.1)
	spn = surfpt_nearby(xy,s)
	spn[1][1] * spn[2][2]
end
function foo3(p::Vector)
	s = ridge_wg(p[1],p[2],p[3],p[6],p[4],p[5],6.,4.)[1]::Polygon{4,8}
	xy = SVector(0.9,0.1)
	spn = surfpt_nearby2(xy,s)
	spn[1][1] * spn[2][2]
end

function foo4(p::Vector)
	s = ridge_wg(p[1],p[2],p[3],p[6],p[4],p[5],6.,4.)[2]
	xy = SVector(0.9,0.1)
	spn = surfpt_nearby(xy,s)
	spn[1][1] * spn[2][2]
end

function foo5(p::Vector)
	s = ridge_wg(p[1],p[2],p[3],p[6],p[4],p[5],6.,4.)[2]::Box{2,4}
	xy = SVector(0.9,0.1)
	spn = surfpt_nearby2(xy,s)
	spn[1][1] * spn[2][2]
end

using OptiMode: εₛ, make_εₛ⁻¹

function goo1(p)
	s::Vector{<:Shape2} = rwg(p)
	εₛ(s,0.9,0.1;tree=KDTree(s),δx=Δx/Nx,δy=Δy/Ny)
end

function goo2(p)
	s::Vector{<:Shape2} = rwg(p)
	εₛ(s,0.9,0.1;tree=KDTree(s),δx=Δx/Nx,δy=Δy/Ny)
end

function eip(shapes::AbstractVector{<:GeometryPrimitives.Shape{2,4,D}},x::Real,y::Real;tree::KDTree,δx::Real,δy::Real,npix_sm::Int=1) where D
    s1 = @ignore(findfirst(SVector(x+δx/2.,y+δy/2.),tree))
    s2 = @ignore(findfirst(SVector(x+δx/2.,y-δy/2.),tree))
    s3 = @ignore(findfirst(SVector(x-δx/2.,y-δy/2.),tree))
    s4 = @ignore(findfirst(SVector(x-δx/2.,y+δy/2.),tree))

    ε1 = isnothing(s1) ? εᵥ : s1.data
    ε2 = isnothing(s2) ? εᵥ : s2.data
    ε3 = isnothing(s3) ? εᵥ : s3.data
    ε4 = isnothing(s4) ? εᵥ : s4.data
	return (ε1,ε2,ε3,ε4)
end

function eip2(shapes::AbstractVector{<:GeometryPrimitives.Shape{2,4,D}},x::Real,y::Real;tree::KDTree,δx::Real,δy::Real,npix_sm::Int=1) where D
	mapres1 = map((op1,op2)->findfirst([op1(x,δx),op2(y,δy)],tree),(+,+,-,-),(+,-,-,+))
	mapres2 = map(x->( isnothing(x) ? εᵥ : x.data ), mapres1 )
	return mapres2
end

##
println("@btime foo1, ∇foo1")
@show foo1(p)
@btime foo1($p)
@show Zygote.gradient(foo1,p)
@btime Zygote.gradient(foo1,$p)

println("@btime foo2, ∇foo2")
@show foo2(p)
@btime foo2($p)
@show Zygote.gradient(foo2,p)
@btime Zygote.gradient(foo2,$p)

println("@btime foo3, ∇foo3")
@show foo3(p)
@btime foo3($p)
@show Zygote.gradient(foo3,p)
@btime Zygote.gradient(foo3,$p)

## polygon step by step
K=4
@show v = v_pgn(p[1],p[2],p[3],p[6],p[4],p[5],Δx,Δy)
@show w = v .- mean(v, dims=1)  # v in center-of-mass coordinates
@show ϕ = mod.(atan.(w[:,2], w[:,1]), 2π)  # SVector{K}: angle of vertices between 0 and 2π; `%` does not work for negative angle
if !issorted(ϕ)	# TODO: make sort_verts shuffling fn with AD rules, currently unsorted verts would break differentiability
    # Do this only when ϕ is not sorted, because the following uses allocations.
    ind = MVector{K}(sortperm(ϕ))  # sortperm(::SVector) currently returns Vector, not MVector
    v = v[ind,:]  # SVector{K}: sorted v
end
# Calculate the increases in angle between neighboring edges.
# ∆v = vcat(diff(v, dims=1), SMatrix{1,2}(v[1,:]-v[end,:]))  # SMatrix{K,2}: edge directions
@show ∆v = vcat(diff(v, dims=1), transpose(v[1,:]-v[end,:]))
@show ∆z = ∆v[:,1] + im * ∆v[:,2]  # SVector{K}: edge directions as complex numbers
@show icurr = ntuple(identity, Val(K-1))
@show inext = ntuple(x->x+1, Val(K-1))
@show ∆ϕ = angle.(∆z[SVector(inext)] ./ ∆z[SVector(icurr)])  # angle returns value between -π and π
@show  n0 = [∆v[:,2] -∆v[:,1]]  # outward normal directions to edges
@show  n = n0 ./ hypot.(n0[:,1],n0[:,2])  # normalize

function foo(v)
	@show w = v .- mean(v, dims=1)  # v in center-of-mass coordinates
	@show ϕ = mod.(atan.(w[:,2], w[:,1]), 2π)  # SVector{K}: angle of vertices between 0 and 2π; `%` does not work for negative angle
	if !issorted(ϕ)	# TODO: make sort_verts shuffling fn with AD rules, currently unsorted verts would break differentiability
	    # Do this only when ϕ is not sorted, because the following uses allocations.
	    ind = MVector{K}(sortperm(ϕ))  # sortperm(::SVector) currently returns Vector, not MVector
	    v = v[ind,:]  # SVector{K}: sorted v
	end

	# Calculate the increases in angle between neighboring edges.
	# ∆v = vcat(diff(v, dims=1), SMatrix{1,2}(v[1,:]-v[end,:]))  # SMatrix{K,2}: edge directions
	@show ∆v = vcat(diff(v, dims=1), transpose(v[1,:]-v[end,:]))
	@show ∆z = ∆v[:,1] + im * ∆v[:,2]  # SVector{K}: edge directions as complex numbers
	@show icurr = ntuple(identity, Val(K-1))
	@show inext = ntuple(x->x+1, Val(K-1))
	@show ∆ϕ = angle.(∆z[SVector(inext)] ./ ∆z[SVector(icurr)])  # angle returns value between -π and π
	@show  n0 = [∆v[:,2] -∆v[:,1]]  # outward normal directions to edges
	@show  n = n0 ./ hypot.(n0[:,1],n0[:,2])  # normalize
end
## box step by step
# spnby
T = Float64
@show x = SVector(2.8,-1.1)
@show b = rwg(p)[2]
@show ax = inv(b.p)  # axes: columns are unit vectors
@show n = (b.p ./ sqrt.(sum(abs2,b.p,dims=2)[:,1]))  # b.p normalized in row direction
@show cosθ = sum(ax.*n', dims=1)[1,:]  # equivalent to diag(n*ax)
# cosθ = diag(n*ax)  # faster than SVector(ntuple(i -> ax[:,i]⋅n[i,:], Val(N)))
# @assert all(cosθ .≥ 0)

@show d = b.p * (x - b.c)
@show n = n .* copysign.(1.0,d)  # operation returns SMatrix (reason for leaving n untransposed)
@show absd = abs.(d)
@show onbnd = abs.(b.r.-absd) .≤ Base.rtoldefault(T) .* b.r  # basically b.r .≈ absd but faster
@show isout = (b.r.<absd) .| onbnd
@show ∆ = (b.r .- absd) .* cosθ  # entries can be negative
if count(isout) == 0  # x strictly inside box; ∆ all positive
	@show l∆x, i = findmin(∆)  # find closest face
	@show nout = n[i,:]
	@show ∆x = l∆x * nout
else  # x outside box or on boundary in one or multiple directions
	@show ∆x = n' * (∆ .* isout)  # project out .!isout directions
	@show nout = all(.!isout .| onbnd) ? n'*onbnd : -∆x
	@show nout = normalize(nout)
end

# bounds
signmatrix(b::Box{1}) = SMatrix{1,1}(1)
signmatrix(b::Box{2}) = SMatrix{2,2}(1,1, -1,1)
signmatrix(b::Box{3}) = SMatrix{3,4}(1,1,1, -1,1,1, 1,-1,1, 1,1,-1)
@show A = inv(b.p) .* b.r'
@show m = maximum(abs.(A * signmatrix(b)), dims=2)[:,1]
## box step by step
K=4
@show v = v_pgn(p[1],p[2],p[3],p[6],p[4],p[5],Δx,Δy)
@show w = v .- mean(v, dims=1)  # v in center-of-mass coordinates
@show ϕ = mod.(atan.(w[:,2], w[:,1]), 2π)  # SVector{K}: angle of vertices between 0 and 2π; `%` does not work for negative angle
if !issorted(ϕ)	# TODO: make sort_verts shuffling fn with AD rules, currently unsorted verts would break differentiability
    # Do this only when ϕ is not sorted, because the following uses allocations.
    ind = MVector{K}(sortperm(ϕ))  # sortperm(::SVector) currently returns Vector, not MVector
    v = v[ind,:]  # SVector{K}: sorted v
end
# Calculate the increases in angle between neighboring edges.
# ∆v = vcat(diff(v, dims=1), SMatrix{1,2}(v[1,:]-v[end,:]))  # SMatrix{K,2}: edge directions
@show ∆v = vcat(diff(v, dims=1), transpose(v[1,:]-v[end,:]))
@show ∆z = ∆v[:,1] + im * ∆v[:,2]  # SVector{K}: edge directions as complex numbers
@show icurr = ntuple(identity, Val(K-1))
@show inext = ntuple(x->x+1, Val(K-1))
@show ∆ϕ = angle.(∆z[SVector(inext)] ./ ∆z[SVector(icurr)])  # angle returns value between -π and π
@show  n0 = [∆v[:,2] -∆v[:,1]]  # outward normal directions to edges
@show  n = n0 ./ hypot.(n0[:,1],n0[:,2])  # normalize

##

@tullio ∆z2[i] := ∆v[i,j] + im * ∆v[i,j]

∆v *  SMatrix{2,1}(1. + 0im, 0 + 1im)

∆v2 = vcat(diff(v, dims=1), SMatrix{1,2}(v[1,1]-v[end,1],v[1,2]-v[end,2]))

@tullio n[i,j] = Δv

∆v2 ≈ ∆v
∆z2 ≈ ∆z
n0 = [∆v[:,2] -∆v[:,1]]  # outward normal directions to edges
n = n0 ./ hypot.(n0[:,1],n0[:,2])  # normalize

function mydiff(v::SMatrix{K,2,T}) where {K,T}
	SMatrix{K,2,T}(vcat(diff(v, dims=1), [v[1,1]-v[K,1]   v[1,2]-v[K,2]]))
end

circdiff0(v) = vcat(diff(v, dims=1), transpose(v[1,:]-v[end,:]))
dcircdiff0(v) = Zygote.gradient(x->circdiff0(x)[1,1],v)
circdiff0(v)
@btime circdiff0($v)
dcircdiff0(v)
@btime dcircdiff0($v)

circdiff1(v) = vcat(diff(v, dims=1), SMatrix{1,2}(v[1,1]-v[end,1],v[1,2]-v[end,2]))
dcircdiff1(v) = Zygote.gradient(x->circdiff1(x)[1,1],v)
circdiff1(v)
@btime circdiff1($v)
dcircdiff1(v)
@btime dcircdiff1($v)
Zygote.gradient(x->circdiff1(x[1,1]),v)

circdiff2(v) = v - circshift(v,1)
dcircdiff2(v) = Zygote.gradient(x->circdiff2(x)[1,1],v)
circdiff2(v)
@btime circdiff2($v)
dcircdiff2(v)
@btime dcircdiff2($v)

diff(vcat(v,v),dims=1)[1:K,1:2] ≈ vcat(diff(v, dims=Val(1)), SMatrix{1,2}(v[1,:]-v[end,:]))
diff(vcat(v,v),dims=1)[1:K,1:2] ≈ vcat(diff(v, dims=1), SMatrix{1,2}(v[1,1]-v[end,1],v[1,2]-v[end,2]))
f1(v) = ( K=size(v)[1]; sum(abs2,diff(vcat(v,v),dims=1)[1:K,1:2]) )
f2(v) = sum(abs2,vcat(diff(v, dims=1), SMatrix{1,2}(v[1,1]-v[end,1],v[1,2]-v[end,2])))
f3(v) = sum(abs2,vcat(diff(v, dims=Val(1)), SMatrix{1,2}(v[1,:]-v[end,:])))
f1(v)
f2(v)
f3(v)
Zygote.gradient(f1,v)
Zygote.gradient(f2,v)
Zygote.gradient(f3,v)

##
function surfpt_nearby2(x::AbstractVector, s::GeometryPrimitives.Sphere{2})
    nout = x==s.c ? SVector(1.0,0.0) : # nout = e₁ for x == s.c
                    normalize(x-s.c)
    return s.c+s.r*nout, nout
end

@btime bounds($s) # 771.861 ns (11 allocations: 848 bytes)
@btime bounds2($s) # 35.546 ns (2 allocations: 64 bytes)

# match original surfpt_nearby(::Polygon) code with vectorized, type-stable, differentiable fns
x = SVector(0.9,0.1)
x = SVector(1.9,1.1)
s = rwg(p)[1]	# point near edge and polygon to use for testing
# original
∆xe = sum(s.n .* (x' .- s.v), dims=Val(2))[:,1]  # SVector{K}: values of equations of edge lines
abs∆xe = abs.(∆xe)  # SVector{K}
# Determine if x is outside of edges, inclusive.
sz = abs.((-)(bounds(s)...))  # SVector{2}
onbnd = abs∆xe .≤ Base.rtoldefault(Float64) * max(sz.data...)  # SVector{K}
isout = (∆xe.>0) .| onbnd  # SVector{K}
cout = count(isout)

∆xv = x' .- s.v
l∆xv = hypot.(∆xv[:,1], ∆xv[:,2])
imin = argmin(l∆xv)
surf = s.v[imin,:]
imin₋₁ = mod1(imin-1,length(l∆xv))

imax = argmax(∆xe)
vmax, nmax = s.v[imax,:], s.n[imax,:]

∆x = (nmax⋅(vmax-x)) .* nmax

surf = x + ∆x
nout = nmax

∆xv = ones(4)*x' - s.v #x' .- s.v::SMatrix{K,2,T,K2}
l∆xv = map(sqrt,mapreduce(abs2,+,∆xv,dims=2))
sum(abs2,Δxv;dims=2) # hypot.(∆xv[1:K,1], ∆xv[1:K,2])
imin = argmin(l∆xv)
surf = @inbounds SVector{2,Float64}(s.v[imin,1],s.v[imin,2])
imin₋₁ = mod1(imin-1,K)

# new
∆xe = f_Δxe(s.n,x,s.v)  # SVector{K}: values of equations of edge lines
# abs∆xe = abs.(∆xe)  # SVector{K}
# Determine if x is outside of edges, inclusive.
# sz = abs.((-)(bounds(s)...))  # SVector{2}
# onbnd = abs∆xe .≤ Base.rtoldefault(Float64) * max(sz.data...)  # SVector{K}
# isout = (∆xe.>0) .| onbnd  # SVector{K}

function τ_trans(ε::AbstractMatrix{T}) where T<:Real
    return @inbounds SMatrix{3,3,T,9}(
        -1/ε[1,1],      ε[2,1]/ε[1,1],                  ε[3,1]/ε[1,1],
        ε[1,2]/ε[1,1],  ε[2,2] - ε[2,1]*ε[1,2]/ε[1,1],  ε[3,2] - ε[3,1]*ε[1,2]/ε[1,1],
        ε[1,3]/ε[1,1],  ε[2,3] - ε[2,1]*ε[1,3]/ε[1,1],  ε[3,3] - ε[3,1]*ε[1,3]/ε[1,1]
    )
end

function τ⁻¹_trans(τ::AbstractMatrix{T}) where T<:Real
    return @inbounds SMatrix{3,3,T,9}(
        -1/τ[1,1],          -τ[2,1]/τ[1,1],                 -τ[3,1]/τ[1,1],
        -τ[1,2]/τ[1,1],     τ[2,2] - τ[2,1]*τ[1,2]/τ[1,1],  τ[3,2] - τ[3,1]*τ[1,2]/τ[1,1],
        -τ[1,3]/τ[1,1],     τ[2,3] - τ[2,1]*τ[1,3]/τ[1,1],  τ[3,3]- τ[3,1]*τ[1,3]/τ[1,1]
    )
end

function avg_param(param1, param2, n12, rvol1)
    n = n12 / norm(n12) #sqrt(sum2(abs2,n12))

    # Pick a vector that is not along n.
    if any(n .== 0)
    	htemp1 = (n .== 0)
    else
    	htemp1 = SVector(1., 0. , 0.)
    end

    # Create two vectors that are normal to n and normal to each other.
    htemp2 = n × htemp1
    h = htemp2 / norm(htemp2) #sqrt(sum2(abs2,htemp2))
    vtemp = n × h
    v = vtemp / norm(vtemp) #sqrt(sum2(abs2,vtemp))
    # Create a local Cartesian coordinate system.
    S = [n h v]  # unitary

    τ1 = τ_trans(transpose(S) * param1 * S)  # express param1 in S coordinates, and apply τ transform
    τ2 = τ_trans(transpose(S) * param2 * S)  # express param2 in S coordinates, and apply τ transform

    τavg = τ1 .* rvol1 + τ2 .* (1-rvol1)  # volume-weighted average

    return S * τ⁻¹_trans(τavg) * transpose(S)  # apply τ⁻¹ and transform back to global coordinates
end

function εₛ(shapes::AbstractVector{<:GeometryPrimitives.Shape{2,4,D}},x::Real,y::Real;tree::KDTree,δx::Real,δy::Real,npix_sm::Int=1)::D where D
    # x1,y1 = x+npix_sm*δx/2.,y+npix_sm*δy/2.
    # x2,y2 = x+npix_sm*δx/2.,y-npix_sm*δy/2.
    # x3,y3 = x-npix_sm*δx/2.,y-npix_sm*δy/2.
    # x4,y4 = x-npix_sm*δx/2.,y+npix_sm*δy/2.
    # x1 = x+δx/2.
    # y1 = y+δy/2.
    # x2 = x+δx/2.
    # y2 = y-δy/2.
    # x3 = x-δx/2.
    # y3 = y-δy/2.
    # x4 = x-δx/2.
    # y4 = y+δy/2.

    s1 = @ignore(findfirst([x+δx/2.,y+δy/2.],tree))
    s2 = @ignore(findfirst([x+δx/2.,y-δy/2.],tree))
    s3 = @ignore(findfirst([x-δx/2.,y-δy/2.],tree))
    s4 = @ignore(findfirst([x-δx/2.,y+δy/2.],tree))

    ε1 = isnothing(s1) ? εᵥ : s1.data
    ε2 = isnothing(s2) ? εᵥ : s2.data
    ε3 = isnothing(s3) ? εᵥ : s3.data
    ε4 = isnothing(s4) ? εᵥ : s4.data

    if (ε1==ε2==ε3==ε4)
        return ε1
    else
        sinds = @ignore ( [ isnothing(ss) ? length(shapes)+1 : findfirst(isequal(ss),shapes) for ss in [s1,s2,s3,s4]] )
        n_unique = @ignore( length(unique(sinds)) )
        if n_unique==2
            s_fg = @ignore(shapes[minimum(dropgrad(sinds))])
            r₀,nout = surfpt_nearby2([x; y], s_fg)
            # bndry_pxl[i,j] = 1
            # nouts[i,j,:] = nout
            # vxl = (SVector{2}(x-δx/2.,y-δy/2.), SVector{2}(x+δx/2.,y+δy/2.))
            rvol = volfrac((SVector{2}(x-δx/2.,y-δy/2.), SVector{2}(x+δx/2.,y+δy/2.)),nout,r₀)
            sind_bg = @ignore(maximum(dropgrad(sinds))) #max(sinds...)
            ε_bg = sind_bg > length(shapes) ? εᵥ : shapes[sind_bg].data
            # return avg_param(
            #         s_fg.data,
            #         ε_bg,
            #         [nout[1];nout[2];0],
            #         rvol,)
            return s_fg.data.^2 + 8*ε_bg
        else
            return (ε1+ε2+ε3+ε4)/4.
        end
    end
end

@adjoint function Broadcast.broadcasted(::typeof(εₛ), x::Union{T,AbstractArray{<:T}}) where {T<:Real}
    y, back = Zygote.broadcast_forward(εₛ, x)
    y, ȳ -> (nothing, back(ȳ)...)
end

function foo(p::Vector{<:Real},x::Real,y::Real;Δx::Real,Δy::Real,Nx::Int,Ny::Int)
    s::Vector{<:Shape2} = rwg(p)
    εₛ(s,x,y;tree=KDTree(s),δx=Δx/Nx,δy=Δy/Ny)
end
foop(x) = foo(x,0.9,0.1;Δx=6.,Δy=4.,Nx=128,Ny=128)[1,1]

function goo1(p::Vector{<:Real},x::Real,y::Real)
    s::Shape2 = rwg(p)[1]
    r₀,nout = surfpt_nearby2([x; y], s)
	return r₀[1] * r₀[2] * nout[1] * nout[2]
end
function goo2(p::Vector{<:Real},x::Real,y::Real)::Real
    s::Shape2 = rwg(p)[2]
    r₀,nout = surfpt_nearby2([x; y], s)
	return r₀[1] * r₀[2]
end
goop1(p) = goo1(p,0.9,0.1)
goop2(p::Vector{<:Real})::Real = goo2(p,0.9,0.1)

function hoo(p::Vector{<:Real},x::Real,y::Real)
    s::Vector{<:Shape2} = rwg(p)[1]
    r₀,nout = surfpt_nearby2([x; y], s)
	volfrac((SVector{2}(x-δx/2.,y-δy/2.), SVector{2}(x+δx/2.,y+δy/2.)),nout,r₀)
end
hoop(p) = goo1(p,0.9,0.1)

pgon = rwg(p)[1] # polygon

function goo(p)
	local xy = SVector(0.9,0.1)
	pgon::Polygon = rwg(p)[1]
	(r,n) = surfpt_nearby(xy,pgon)
	return r[1]
end

dgoo(p) = Zygote.gradient(goo,p)

goo(p)
dgoo(p)
@btime goo($p)

_,spn_pb = Zygote._pullback(Context(),(x,y)->( (r,n) = surfpt_nearby(x,y) ; return (r[1] * r[2] * n[1] * n[2])), SVector(0.9,0.1),pgon)

es_pb(1)
es_pb2(1)

foop(p)
goop1(p)
goop2(p)
hoop(p)

ForwardDiff.gradient(goop2,p)

p̄_FD3 = FiniteDifferences.grad(central_fdm(3,1),foop,p)[1]
@test Zygote.gradient(foop,p)[1] ≈ p̄_FD3
@test ForwardDiff.gradient(foop,p) ≈ p̄_FD3

@btime Zygote.gradient(foop,$p)[1]
@btime ForwardDiff.gradient(foop,$p)

@btime Zygote.gradient(goop1,$p)[1]
@btime ForwardDiff.gradient(goop1,$p)
@btime Zygote.gradient(goop2,$p)[1]
@btime ForwardDiff.gradient(goop2,$p)
@btime Zygote.gradient(hoop,$p)[1]
@btime ForwardDiff.gradient(hoop,$p)

##
M = randn(Float64,3,3)
Ms = SMatrix{3,3,Float64,9}(M)

##
a,b,c = randn(),randn(),randn()
v = [a,b,c]
vs = SVector(a,b,c)
norm(v)
norm(vs)
Zygote.gradient(norm,v)[1]
Zygote.gradient(norm,vs)[1]
@btime Zygote.gradient(norm,$v)[1]
@btime Zygote.gradient(norm,$vs)[1]

vsr = reinterpret(reshape,vs)

_, norm_pb = Zygote._pullback(Context(),norm,v)
_, norm_pb = Zygote._pullback(Context(),norm,vs)

normalize(v)
normalize(vs)
Zygote.gradient(x->normalize(x)[1],v)[1]
Zygote.gradient(x->normalize(x)[1],vs)[1]

@btime normalize($v)
@btime normalize($vs)
@btime Zygote.gradient(x->normalize(x)[1],$v)[1]
@btime Zygote.gradient(x->normalize(x)[1],$vs)[1]


text:
function adjoint_solve(prob, sol_init, u₀, tstops, p, rtol, atol)
    de = ModelingToolkit.modelingtoolkitize(prob)
    sym_jac = eval(ModelingToolkit.generate_jacobian(de)[1])
    function adjoint!(du, u, p, t)
        du = -Base.invokelatest(sym_jac, sol_init(t), p, t)' * u
    end
    rev_tstops = reverse(tstops) # reverse tstops for λ(T) = u₀?
    rev_tspan = (rev_tstops[1], rev_tstops[end])
    adj_prob = ODEProblem(adjoint!, u₀, rev_tspan, p)
    λ = solve(adj_prob, sol_init.alg, reltol=rtol, abstol=atol, tstops=rev_tstops)
    return λ
end

function condition_num(λₜ)
    K = LinearAlgebra.norm(λₜ, 1) + LinearAlgebra.norm(λₜ[1])
    return K
end

function adjoint_step_ctrl(λ₁ₜ, λ₂ₜ, E₂, Eₙ, C, gtol)
    K₁ = condition_num(λ₁ₜ)
    K₂ = condition_num(λ₂ₜ)
    K = (E₂ / Eₙ) * sqrt(K₁^2 + K₂^2)
    h = sqrt((gtol / (K * C))) # T in denominator?
    return h
end


text:
abstract type Beverage end

# template
function prepare(b::Beverage)
    boil_water()
    brew!(b)
    pour_in_cup(b)
    add_condiments!(b)
end

boil_water() = println("Boil water...")

function brew! end

function pour_in_cup end

function add_condiments! end

mutable struct Coffee <: Beverage
    content::String
    condiment::String
end

function brew!(c::Coffee)
    c.content = "coffee"
    println("Brew $(c.content)...")
end

pour_in_cup(::Coffee) = println("Pour coffee in cup...")

function add_condiments!(c::Coffee)
    c.condiment = "milk and sugar"
    println("Add $(c.condiment)...")
end

mutable struct Tea <: Beverage
    content::String
    condiment::String
end

function brew!(t::Tea)
    t.content = "black tea"
    println("Brew $(t.content)...")
end

pour_in_cup(::Tea) = println("Pour tea in cup.")

function add_condiments!(t::Tea)
    t.condiment = "lemon"
    println("Add $(t.condiment)...")
end


text:
### A Pluto.jl notebook ###
# v0.14.3

using Markdown
using InteractiveUtils

# ╔═╡ 40650588-5eb8-450d-824b-ac509bf7894f
using Pkg, DrWatson

# ╔═╡ 9f979b0d-bf6f-4d6f-8def-5be5926711a5
begin
	@quickactivate "NumericalMethodsforEngineers"
	using NumericalMethodsforEngineers
	using Test
end

# ╔═╡ fadae5bc-da5d-49c2-9e1b-fa239bf6577b
md" ### f-02-02.jl"

# ╔═╡ 32ddb193-70a2-4650-8c67-d5501461cf58
begin
	a = [2. -3. 1.; 1. -1. -2.; 3. 1. -1.]
	b = [7., -2., 0.]
end

# ╔═╡ 9fa8c5ec-2c7a-49ac-9cf9-7bfdb7650d01
md" ##### Julia version:"

# ╔═╡ bdc5eb8f-07a6-40df-b2ba-08903f896cb0
c = a\b

# ╔═╡ 44e675da-00a6-45c9-a47f-ac8195aa87f8
f = lu(a)

# ╔═╡ 25800554-5b76-48bc-96bc-0084a3a9177e
md" ##### Fortran version:"

# ╔═╡ af7dba26-b167-4d17-8d94-c8d4accb5d86
lower, upper = lufac(a)

# ╔═╡ 60b977ad-8e1a-49f1-b0d1-51c4f61037dc
subfor(lower, b)

# ╔═╡ f0b8a692-8ad2-4ea9-9dad-9bd5d2cd612b
begin
	subfor!(lower, b)
	b
end

# ╔═╡ 96f0435c-9b8d-48d5-97eb-1c434208c184
subbac(upper, b)

# ╔═╡ c29fed2a-ad84-4873-a550-9d65bd187c83
begin
	subbac!(upper, b)
	b
end

# ╔═╡ a30f098e-4ea3-4ec4-9fa3-373f8c17eb5c
md" ##### Solution vector:"

# ╔═╡ 3c89317c-1154-49bd-93ce-09c631014bbb
round.(a * c; digits=14) == b

# ╔═╡ f9921b74-c79a-4e2a-8375-03e983462875
md" ### End of f-02-02.jl"

# ╔═╡ Cell order:
# ╟─fadae5bc-da5d-49c2-9e1b-fa239bf6577b
# ╠═40650588-5eb8-450d-824b-ac509bf7894f
# ╠═9f979b0d-bf6f-4d6f-8def-5be5926711a5
# ╠═32ddb193-70a2-4650-8c67-d5501461cf58
# ╟─9fa8c5ec-2c7a-49ac-9cf9-7bfdb7650d01
# ╠═bdc5eb8f-07a6-40df-b2ba-08903f896cb0
# ╠═44e675da-00a6-45c9-a47f-ac8195aa87f8
# ╟─25800554-5b76-48bc-96bc-0084a3a9177e
# ╠═af7dba26-b167-4d17-8d94-c8d4accb5d86
# ╠═60b977ad-8e1a-49f1-b0d1-51c4f61037dc
# ╠═f0b8a692-8ad2-4ea9-9dad-9bd5d2cd612b
# ╟─a30f098e-4ea3-4ec4-9fa3-373f8c17eb5c
# ╠═96f0435c-9b8d-48d5-97eb-1c434208c184
# ╠═c29fed2a-ad84-4873-a550-9d65bd187c83
# ╠═3c89317c-1154-49bd-93ce-09c631014bbb
# ╟─f9921b74-c79a-4e2a-8375-03e983462875


text:
using Distributed
addprocs(16)
@everywhere include("get_objective_lorenz.jl")
get_Javg_vs_s()

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
# ---
# title: Faceting
# cover: assets/faceting.png
# description: Generating a grid of plots
# author: "[Pietro Vertechi](https://github.com/piever)"
# ---

using AlgebraOfGraphics, CairoMakie
set_aog_theme!() #src

# ## Facet grid

df = (x=rand(100), y=rand(100), i=rand(["a", "b", "c"], 100), j=rand(["d", "e", "f"], 100))
plt = data(df) * mapping(:x, :y, col=:i, row=:j)
draw(plt)

# ## Facet wrap

df = (x=rand(100), y=rand(100), l=rand(["a", "b", "c", "d", "e"], 100))
plt = data(df) * mapping(:x, :y, layout=:l)
draw(plt)

# ## Adding traces to only some subplots

df1 = (x=rand(100), y=rand(100), i=rand(["a", "b", "c"], 100), j=rand(["d", "e", "f"], 100))
df2 = (x=[0, 1], y=[0.5, 0.5], i=fill("a", 2), j=fill("e", 2))
layers = data(df1) * visual(Scatter) + data(df2) * visual(Lines)
fg = draw(layers * mapping(:x, :y, col=:i, row=:j))

# save cover image #src
mkpath("assets") #src
save("assets/faceting.png", fg) #src


text:
#### n-stream kernels
#  - delta-scaling
#  - adding two sets of properties
#  - extracting subsets
#  - validity checking

# ---------------------------------
# increment 1scalar by nstream
function increment_by_gpoint!(op_1::OneScalar{FT}, op_2::NStream{FT}) where FT
#    integer,                             intent(in   ) :: ncol, nlay, ngpt
#    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau
#    real(FT), dimension(ncol,nlay,ngpt), intent(in   ) :: op_2.tau, op_2.ssa

#    integer  :: icol, ilay, igpt

  op_1.tau = op_1.tau .+ op_2.tau * ( FT(1) - op_2.ssa)
end

# ---------------------------------
# increment 2stream by nstream
function increment_by_gpoint!(op_1::TwoStream{FT}, op_2::NStream{FT}) where FT
#    integer,                              intent(in   ) :: ncol, nlay, ngpt, nmom2
#    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau, op_1.ssa, op_1.g
#    real(FT), dimension(ncol,nlay,ngpt), intent(in   ) :: op_2.tau, op_2.ssa
#    real(FT), dimension(nmom2,
#                        ncol,nlay,ngpt), intent(in   ) :: op_2.p

#    integer  :: icol, ilay, igpt
#    real(FT) :: tau12, tauscat12

  ngpt = get_ngpt(op_1)
  nlay = get_nlay(op_1)
  ncol = get_ncol(op_1)

  for igpt = 1:ngpt
    for ilay = 1:nlay
      for icol = 1:ncol
        # t=op_1.tau + op_2.tau
        tau12 = op_1.tau[icol,ilay,igpt] + op_2.tau[icol,ilay,igpt]
        # w=(op_1.tau*op_1.ssa + op_2.tau*op_2.ssa) / t
        tauscat12 =
           op_1.tau[icol,ilay,igpt] * op_1.ssa[icol,ilay,igpt] +
           op_2.tau[icol,ilay,igpt] * op_2.ssa[icol,ilay,igpt]
        op_1.g[icol,ilay,igpt] =
          (op_1.tau[icol,ilay,igpt] * op_1.ssa[icol,ilay,igpt] * op_1.g[   icol,ilay,igpt]+
           op_2.tau[icol,ilay,igpt] * op_2.ssa[icol,ilay,igpt] * op_2.p[1, icol,ilay,igpt]) / max(ε_machine(FT),tauscat12)
        op_1.ssa[icol,ilay,igpt] = tauscat12 / max(ε_machine(FT),tau12)
        op_1.tau[icol,ilay,igpt] = tau12
      end
    end
  end
end

# increment nstream by 1scalar
function increment_by_gpoint!(op_1::NStream{FT}, op_2::OneScalar{FT}) where FT
#    integer,                              intent(in   ) :: ncol, nlay, ngpt
#    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau, op_1.ssa
#    real(FT), dimension(ncol,nlay,ngpt), intent(in   ) :: op_2.tau

#    integer  :: icol, ilay, igpt
#    real(FT) :: tau12

  ngpt = get_ngpt(op_1)
  nlay = get_nlay(op_1)
  ncol = get_ncol(op_1)

  for igpt = 1:ngpt
    for ilay = 1:nlay
      for icol = 1:ncol
        tau12 = op_1.tau[icol,ilay,igpt] + op_2.tau[icol,ilay,igpt]
        op_1.ssa[icol,ilay,igpt] = op_1.tau[icol,ilay,igpt] * op_1.ssa[icol,ilay,igpt] / max(ε_machine(FT),tau12)
        op_1.tau[icol,ilay,igpt] = tau12
        # p is unchanged
      end
    end
  end
end

# ---------------------------------
# increment nstream by 2stream
function increment_by_gpoint!(op_1::NStream{FT}, op_2::TwoStream{FT}) where FT
#    integer,                              intent(in   ) :: ncol, nlay, ngpt, nmom1
#    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau, op_1.ssa
#    real(FT), dimension(nmom1,
#                        ncol,nlay,ngpt), intent(inout) :: op_1.p
#    real(FT), dimension(ncol,nlay,ngpt), intent(in   ) :: op_2.tau, op_2.ssa, op_2.g

#    integer  :: icol, ilay, igpt
#    real(FT) :: tau12, tauscat12
#    real(FT), dimension(nmom1) :: temp_moms # TK
#    integer  :: imom  #TK

  ngpt = get_ngpt(op_1)
  nlay = get_nlay(op_1)
  ncol = get_ncol(op_1)

  for igpt = 1:ngpt
    for ilay = 1:nlay
      for icol = 1:ncol
        tau12 = op_1.tau[icol,ilay,igpt] + op_2.tau[icol,ilay,igpt]
        tauscat12 =
           op_1.tau[icol,ilay,igpt] * op_1.ssa[icol,ilay,igpt] +
           op_2.tau[icol,ilay,igpt] * op_2.ssa[icol,ilay,igpt]
        #
        # Here assume Henyey-Greenstein
        #
        temp_moms[1] = op_2.g[icol,ilay,igpt]
        for imom = 2:nmom1
          temp_moms[imom] = temp_moms[imom-1] * op_2.g[icol,ilay,igpt]
        end
        op_1.p[1:nmom1, icol,ilay,igpt] =
            (op_1.tau[icol,ilay,igpt] * op_1.ssa[icol,ilay,igpt] * op_1.p[1:nmom1, icol,ilay,igpt] +
             op_2.tau[icol,ilay,igpt] * op_2.ssa[icol,ilay,igpt] * temp_moms[1:nmom1]  ) / max(ε_machine(FT),tauscat12)
        op_1.ssa[icol,ilay,igpt] = tauscat12 / max(ε_machine(FT),tau12)
        op_1.tau[icol,ilay,igpt] = tau12
      end
    end
  end
end

# ---------------------------------
# increment nstream by nstream
function increment_by_gpoint!(op_1::NStream{FT}, op_2::NStream{FT}) where FT
#    integer,                              intent(in   ) :: ncol, nlay, ngpt, nmom1, nmom2
#    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau, op_1.ssa
#    real(FT), dimension(nmom1,
#                        ncol,nlay,ngpt), intent(inout) :: op_1.p
#    real(FT), dimension(ncol,nlay,ngpt), intent(in   ) :: op_2.tau, op_2.ssa
#    real(FT), dimension(nmom2,
#                        ncol,nlay,ngpt), intent(in   ) :: op_2.p

#    integer  :: icol, ilay, igpt, mom_lim
#    real(FT) :: tau12, tauscat12

  ngpt = get_ngpt(op_1)
  nlay = get_nlay(op_1)
  ncol = get_ncol(op_1)

  mom_lim = min(nmom1, nmom2)
  for igpt = 1:ngpt
    for ilay = 1:nlay
      for icol = 1:ncol
        tau12 = op_1.tau[icol,ilay,igpt] + op_2.tau[icol,ilay,igpt]
        tauscat12 =
           op_1.tau[icol,ilay,igpt] * op_1.ssa[icol,ilay,igpt] +
           op_2.tau[icol,ilay,igpt] * op_2.ssa[icol,ilay,igpt]
        #
        # If op2 has more moments than op1 these are ignored;
        #   if it has fewer moments the higher orders are assumed to be 0
        #
        op_1.p[1:mom_lim, icol,ilay,igpt] =
            (op_1.tau[icol,ilay,igpt] * op_1.ssa[icol,ilay,igpt] * op_1.p[1:mom_lim, icol,ilay,igpt] +
             op_2.tau[icol,ilay,igpt] * op_2.ssa[icol,ilay,igpt] * op_2.p[1:mom_lim, icol,ilay,igpt]) / max( ε_machine(FT),tauscat12)
        op_1.ssa[icol,ilay,igpt] = tauscat12 / max(eps,tau12)
        op_1.tau[icol,ilay,igpt] = tau12
      end
    end
  end
end

# increment 1scalar by nstream
function increment_bybnd!(op_1::OneScalar{FT}, op_2::NStream{FT}) where FT
#    integer,                             intent(in   ) :: ncol, nlay, ngpt, nbnd
#    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: tau1
#    real(FT), dimension(ncol,nlay,nbnd), intent(in   ) :: tau2, ssa2
#    integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims # Starting and ending gpoint for each band

#    integer :: ibnd, igpt

  ngpt = get_ngpt(op_1)
  nlay = get_nlay(op_1)
  ncol = get_ncol(op_1)
  gpt_lims = get_band_lims_gpoint(op_1)
  nbnd = get_nband(op_1)

  for ibnd = 1:nbnd
    for igpt = gpt_lims[1, ibnd]:gpt_lims[2, ibnd]
      op_1.tau[:,:,igpt] .= op_1.tau[:,:,igpt] .+ op_2.tau[:,:,ibnd] .* ( FT(1) .- op_2.ssa[:,:,ibnd])
    end
  end
end

"""
    increment_bybnd!(op_1::TwoStream{FT}, op_2::NStream{FT}) where FT

integer,                             intent(in   ) :: ncol, nlay, ngpt, nmom2, nbnd
real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau, op_1.ssa, op_1.g
real(FT), dimension(ncol,nlay,nbnd), intent(in   ) :: op_2.tau, op_2.ssa
real(FT), dimension(nmom2,
                    ncol,nlay,nbnd), intent(in   ) :: op_2.p
integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims # Starting and ending gpoint for each band

integer  :: icol, ilay, igpt, ibnd
real(FT) :: tau12, tauscat12
"""
function increment_bybnd!(op_1::TwoStream{FT}, op_2::NStream{FT}) where FT

  ngpt = get_ngpt(op_1)
  nlay = get_nlay(op_1)
  ncol = get_ncol(op_1)
  gpt_lims = get_band_lims_gpoint(op_1)
  nbnd = get_nband(op_1)

  for ibnd = 1:nbnd
    for igpt = gpt_lims[1, ibnd]:gpt_lims[2, ibnd]
      for ilay = 1:nlay
        for icol = 1:ncol
          # t=op_1.tau + op_2.tau
          tau12 = op_1.tau[icol,ilay,igpt] + op_2.tau[icol,ilay,ibnd]
          # w=(op_1.tau*op_1.ssa + op_2.tau*op_2.ssa) / t
          tauscat12 =
             op_1.tau[icol,ilay,igpt] * op_1.ssa[icol,ilay,igpt] +
             op_2.tau[icol,ilay,ibnd] * op_2.ssa[icol,ilay,ibnd]
          op_1.g[icol,ilay,igpt] =
            (op_1.tau[icol,ilay,igpt] * op_1.ssa[icol,ilay,igpt] * op_1.g[   icol,ilay,igpt]+
             op_2.tau[icol,ilay,ibnd] * op_2.ssa[icol,ilay,ibnd] * op_2.p[1, icol,ilay,ibnd]) / max(ε_machine(FT),tauscat12)
          op_1.ssa[icol,ilay,igpt] = tauscat12 / max(ε_machine(FT),tau12)
          op_1.tau[icol,ilay,igpt] = tau12
        end
      end
    end
  end
end

# increment nstream by 1scalar
function increment_bybnd!(op_1::NStream{FT}, op_2::OneScalar{FT}) where FT
#    integer,                             intent(in   ) :: ncol, nlay, ngpt, nbnd
#    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau, op_1.ssa
#    real(FT), dimension(ncol,nlay,nbnd), intent(in   ) :: op_2.tau
#    integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims # Starting and ending gpoint for each band

#    integer  :: icol, ilay, igpt, ibnd
#    real(FT) :: tau12

  ngpt = get_ngpt(op_1)
  nlay = get_nlay(op_1)
  ncol = get_ncol(op_1)
  gpt_lims = get_band_lims_gpoint(op_1)
  nbnd = get_nband(op_1)

  for ibnd = 1:nbnd
    for igpt = gpt_lims[1, ibnd]:gpt_lims[2, ibnd]
      for ilay = 1:nlay
        for icol = 1:ncol
          tau12 = op_1.tau[icol,ilay,igpt] + op_2.tau[icol,ilay,ibnd]
          op_1.ssa[icol,ilay,igpt] = op_1.tau[icol,ilay,igpt] * op_1.ssa[icol,ilay,igpt] / max(ε_machine(FT),tau12)
          op_1.tau[icol,ilay,igpt] = tau12
          # p is unchanged
        end
      end
    end
  end
end

# ---------------------------------
# increment nstream by 2stream
function increment_bybnd!(op_1::NStream{FT}, op_2::TwoStream{FT}) where FT
#    integer,                             intent(in   ) :: ncol, nlay, ngpt, nmom1, nbnd
#    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau, op_1.ssa
#    real(FT), dimension(nmom1,
#                        ncol,nlay,ngpt), intent(inout) :: op_1.p
#    real(FT), dimension(ncol,nlay,nbnd), intent(in   ) :: op_2.tau, op_2.ssa, op_2.g
#    integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims # Starting and ending gpoint for each band

#    integer  :: icol, ilay, igpt, ibnd
#    real(FT) :: tau12, tauscat12
#    real(FT), dimension(nmom1) :: temp_moms # TK
#    integer  :: imom  #TK

  ngpt = get_ngpt(op_1)
  nlay = get_nlay(op_1)
  ncol = get_ncol(op_1)
  gpt_lims = get_band_lims_gpoint(op_1)
  nbnd = get_nband(op_1)

  for ibnd = 1:nbnd
    for igpt = gpt_lims[1, ibnd]:gpt_lims[2, ibnd]
      for ilay = 1:nlay
        for icol = 1:ncol
          tau12 = op_1.tau[icol,ilay,igpt] + op_2.tau[icol,ilay,ibnd]
          tauscat12 =
             op_1.tau[icol,ilay,igpt] * op_1.ssa[icol,ilay,igpt] +
             op_2.tau[icol,ilay,ibnd] * op_2.ssa[icol,ilay,ibnd]
          #
          # Here assume Henyey-Greenstein
          #
          temp_moms[1] = op_2.g[icol,ilay,ibnd]
          for imom = 2:nmom1
            temp_moms[imom] = temp_moms[imom-1] * op_2.g[icol,ilay,ibnd]
          end
          op_1.p[1:nmom1, icol,ilay,igpt] =
              (op_1.tau[icol,ilay,igpt] * op_1.ssa[icol,ilay,igpt] * op_1.p[1:nmom1, icol,ilay,igpt] +
               op_2.tau[icol,ilay,ibnd] * op_2.ssa[icol,ilay,ibnd] * temp_moms[1:nmom1]  ) / max(ε_machine(FT),tauscat12)
          op_1.ssa[icol,ilay,igpt] = tauscat12 / max(ε_machine(FT),tau12)
          op_1.tau[icol,ilay,igpt] = tau12
        end
      end
    end
  end
end

# ---------------------------------
# increment nstream by nstream
function increment_bybnd!(op_1::NStream{FT}, op_2::NStream{FT}) where FT
#    integer,                             intent(in   ) :: ncol, nlay, ngpt, nmom1, nmom2, nbnd
#    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau, op_1.ssa
#    real(FT), dimension(nmom1,
#                        ncol,nlay,ngpt), intent(inout) :: op_1.p
#    real(FT), dimension(ncol,nlay,nbnd), intent(in   ) :: op_2.tau, op_2.ssa
#    real(FT), dimension(nmom2,
#                        ncol,nlay,nbnd), intent(in   ) :: op_2.p
#    integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims # Starting and ending gpoint for each band

#    integer  :: icol, ilay, igpt, ibnd, mom_lim
#    real(FT) :: tau12, tauscat12

  ngpt = get_ngpt(op_1)
  nlay = get_nlay(op_1)
  ncol = get_ncol(op_1)
  gpt_lims = get_band_lims_gpoint(op_1)
  nbnd = get_nband(op_1)

  mom_lim = min(nmom1, nmom2)
  for ibnd = 1:nbnd
    for igpt = gpt_lims[1, ibnd]:gpt_lims[2, ibnd]
      for ilay = 1:nlay
        for icol = 1:ncol
          tau12 = op_1.tau[icol,ilay,igpt] + op_2.tau[icol,ilay,ibnd]
          tauscat12 =
             op_1.tau[icol,ilay,igpt] * op_1.ssa[icol,ilay,igpt] +
             op_2.tau[icol,ilay,ibnd] * op_2.ssa[icol,ilay,ibnd]
          #
          # If op2 has more moments than op1 these are ignored;
          #   if it has fewer moments the higher orders are assumed to be 0
          #
          op_1.p[1:mom_lim, icol,ilay,igpt] =
              (op_1.tau[icol,ilay,igpt] * op_1.ssa[icol,ilay,igpt] * op_1.p[1:mom_lim, icol,ilay,igpt] +
               op_2.tau[icol,ilay,ibnd] * op_2.ssa[icol,ilay,ibnd] * op_2.p[1:mom_lim, icol,ilay,ibnd]) / max( ε_machine(FT),tauscat12)
          op_1.ssa[icol,ilay,igpt] = tauscat12 / max(eps,tau12)
          op_1.tau[icol,ilay,igpt] = tau12
        end
      end
    end
  end
end


text:
module RocketRerunOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: rerun()" begin

    println("Testing: operator rerun()")

    run_proxyshowcheck("Rerun", rerun(1), args = (check_subscription = true, ))

    run_testset([
        (
            source = from(1:5) |> safe() |> map(Int, (d) -> d == 4 ? throw(4) : d) |> rerun(3),
            values = @ts([ 1:3, 1:3, 1:3, 1:3, e(4) ])
        ),
        (
            source = from(1:5) |> safe() |> map(Int, (d) -> d == 4 ? throw(4) : d) |> rerun(3) |> take(6),
            values = @ts([ 1:3, 1:3, c ])
        ),
        (
            source = timer(100, 30) |> safe() |> map(Int, (d) -> d > 2 ? throw("d") : d) |> rerun(1),
            values = @ts(100 ~ [0] ~ 30 ~ [1] ~ 30 ~ [2] ~ 100 ~ [0] ~ 30 ~ [1] ~ 30 ~ [2] ~ e("d"))
        ),
        (
            source = timer(100, 30) |> safe() |> switch_map(Int, (d) -> d > 1 ? faulted(Int, "$d") : of(d)) |> rerun(2),
            values = @ts(100 ~ [0] ~ 30 ~ [1] ~ 100 ~ [0] ~ 30 ~ [1] ~ 100 ~ [0] ~ 30 ~ [1] ~ e("2"))
        ),
        (
            source = completed() |> rerun(2),
            values = @ts(c)
        ),
        (
            source = faulted("e") |> rerun(2),
            values = @ts(e("e")),
        ),
        (
            source = never() |> rerun(2),
            values = @ts()
        )
    ])

end

end


text:
export get_outputs, set_inputs, process

function get_outputs(ind::CGPInd)::Array{Float64}
    # doesn't re-process, just gives outputs
    outputs = Array{Float64}(undef, length(ind.outputs))
    for i in eachindex(outputs)
        outputs[i] = ind.buffer[ind.outputs[i]]
    end
    outputs
end

function set_inputs(ind::CGPInd, inputs::Array{Float64})::Nothing
    for i in eachindex(inputs)
        ind.buffer[i] = inputs[i]
    end
end

function process(ind::CGPInd)::Array{Float64}
    for i in eachindex(ind.nodes)
        n = ind.nodes[i]
        if n.active
            ind.buffer[i] = n.f(ind.buffer[n.x], ind.buffer[n.y])
        end
    end
    get_outputs(ind)
end

function process(ind::CGPInd, inputs::Array{Float64})::Array{Float64}
    set_inputs(ind, inputs)
    process(ind)
end


text:
@testset "AuthTokenTestCase" begin

    @testset "test_should_generate_kerberos_auth_token_correctly" begin
        auth = kerberos_auth("I am a base64 service ticket")
        @test auth["scheme"] == "kerberos"
        @test auth["principal"] == ""
        @test auth["credentials"] == "I am a base64 service ticket"
        @test !haskey(auth, "realm")
        @test !haskey(auth, "parameters")
    end

    @testset "test_should_generate_basic_auth_without_realm_correctly" begin
        auth = basic_auth("molly", "meoooow")
        @test auth["scheme"] == "basic"
        @test auth["principal"] == "molly"
        @test auth["credentials"] == "meoooow"
        @test !haskey(auth, "realm")
        @test !haskey(auth, "parameters")
    end

    @testset "test_should_generate_base_auth_with_realm_correctly" begin
        auth = basic_auth("molly", "meoooow", "cat_cafe")
        @test auth["scheme"] == "basic"
        @test auth["principal"] == "molly"
        @test auth["credentials"] == "meoooow"
        @test auth["realm"] == "cat_cafe"
        @test !haskey(auth, "parameters")
    end

    @testset "test_should_generate_custom_auth_correctly" begin
        auth = custom_auth("cat", "molly", "meoooow", "cat_cafe", age="1", color="white")
        @test auth["scheme"] == "cat"
        @test auth["principal"] == "molly"
        @test auth["credentials"] == "meoooow"
        @test auth["realm"] == "cat_cafe"
        @test auth["parameters"] == Dict("age"=>"1", "color"=>"white")
    end
end


text:
function dbern(t, order)
out = [];
for o = 0:1:order
   if o ==0
      out = db(t,order,o)
   else
      out = [out db(t,order,o)];
   end

end
return out
end


text:
using HTTP
using Gumbo
using Cascadia
using URIs
using Dates

# Get the titles and URLs of Zenn's Julia tag
# ZennのJuliaタグの文書のタイトルとURLを取得
function get_zenn()
    zenn_url = "https://zenn.dev/topics/julia?order=latest"
    # get a html page : HTMLを取得
    response = HTTP.request("GET", zenn_url)
    if response.status != 200
        return nothing
    end

    code = String(response.body)
    doc = Gumbo.parsehtml(code)

    # get count of children : 子要素の数を取得
    len = length(doc.root.children)

    # get contents : ... and ...
    head = nothing
    body = nothing
    for i in 1:len
        if tag(doc.root[i]) == :head
            head = doc.root[i]
        elseif tag(doc.root[i]) == :body
            body = doc.root[i]
        end
    end

    # get main contents
    # 本文を取得
    list = []
    if !isnothing(body)
        # Get by specifying css selector -> get an array of matching elements
        # セレクターを指定して取得　→　一致する要素の配列を返す
        qs = eachmatch(Selector(".ArticleList_link__vf_6E"), body)
        for s in qs
            url = getattr(s, "href", nothing)
            if !isnothing(url)
                u = resolvereference(zenn_url, url)
                url_str = "$u"
                title = nodeText(s)
                pair = (title, url_str)
                push!(list, pair)
            end
        end
    end
    list
end

function main()
    # 日付付きのファイル名を作成
    filename = "zenn" * Dates.format(now(), "yyyymmdd") * ".txt"
    open(filename, "w") do f
        list = get_zenn()
        if isnothing(list)
            println(f, "Failed to get HTML.")
        else
            for (title, url) in list
                println(f, title * "\t" * url)
            end
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

max_stars_repo_name: julialangjp/ScrapingZennTopicsByJulia

max_issues_repo_name: julialangjp/ScrapingZennTopicsByJulia

max_forks_repo_name: julialangjp/ScrapingZennTopicsByJulia


text:
@testset "$(rpad("Dataseries",80))" begin

    dt = Float64
    nt = 10

    nd = 1
    ni = 1
    ds = SDataSeries(dt, nt)
    @test typeof(ds) != typeof(SDataSeries(dt, nd, nt))
    @test typeof(ds) != typeof(SDataSeries(dt, nd, nt, ni))
    @test typeof(ds) <: AbstractArray{dt,1}
    @test firstindex(ds)   == 0
    @test firstindex(ds,1) == firstindex(ds.d,1) - 1
    @test firstindex(ds,2) == 1
    @test lastindex(ds)    == nt
    @test lastindex(ds,1)  == lastindex(ds.d,1) - 1
    @test lastindex(ds,2)  == 1
    @test axes(ds)   == (0:nt,)
    @test axes(ds,1) == 0:nt
    @test axes(ds,2) == 1:1
    @test size(ds.d) == (nt+1,)
    @test size(ds.d) == size(ds)
    @test ndims(ds)  == 1
    @test eltype(ds) == dt
    @test parent(ds) == ds.d

    for i in 0:nt
        ds[i] = i
    end

    @test ds.d[1] == ds[0]
    @test ds.d[end] == ds[nt]
    @test ds.d[end] == ds[end]

    reset!(ds)
    @test ds[0] == ds[end]

    nd = 2
    ni = 1
    ds = SDataSeries(dt, nd, nt, ni)
    @test typeof(ds) == typeof(SDataSeries(dt, nd, nt))
    @test typeof(ds) <: AbstractArray{dt,2}
    @test firstindex(ds)   == firstindex(ds.d)
    @test firstindex(ds,1) == firstindex(ds.d,1)
    @test firstindex(ds,2) == firstindex(ds.d,2) - 1
    @test firstindex(ds,3) == 1
    @test lastindex(ds)    == lastindex(ds.d)
    @test lastindex(ds,1)  == lastindex(ds.d,1)
    @test lastindex(ds,2)  == lastindex(ds.d,2) - 1
    @test lastindex(ds,3)  == 1
    @test axes(ds)   == (1:nd, 0:nt)
    @test axes(ds,1) == 1:nd
    @test axes(ds,2) == 0:nt
    @test axes(ds,3) == 1:1
    @test size(ds.d) == (nd, nt+1)
    @test size(ds.d) == size(ds)
    @test ndims(ds)  == 2
    @test eltype(ds) == dt
    @test parent(ds) == ds.d

    for i in 0:nt
        ds[1,i] = i
        ds[2,i] = 2^i
    end

    @test ds.d[1,1] == ds[1,0]
    @test ds.d[1:ds.nd,1] == ds[1:ds.nd,0]

    @test ds[1,end] == ds.d[1,end]
    @test ds[1,end] == ds[1,nt]
    @test ds[1:ds.nd,nt] == ds[1:ds.nd,end]

    @test ds.d[:,1] == ds[:,0]

    reset!(ds)
    @test ds[1,0] == ds[1,end]

    tx = rand(nd)
    ds[:,0] .= tx
    @test ds.d[:,1] == tx

    nd = 2
    ni = 2
    ds = SDataSeries(dt, nd, nt, ni)
    @test typeof(ds) <: AbstractArray{dt,3}
    @test firstindex(ds)   == firstindex(ds.d)
    @test firstindex(ds,1) == firstindex(ds.d,1)
    @test firstindex(ds,2) == firstindex(ds.d,2) - 1
    @test firstindex(ds,3) == firstindex(ds.d,3)
    @test firstindex(ds,4) == 1
    @test lastindex(ds)    == lastindex(ds.d)
    @test lastindex(ds,1)  == lastindex(ds.d,1)
    @test lastindex(ds,2)  == lastindex(ds.d,2) - 1
    @test lastindex(ds,3)  == lastindex(ds.d,3)
    @test lastindex(ds,4)  == 1
    @test axes(ds)   == (1:nd, 0:nt, 1:ni)
    @test axes(ds,1) == 1:nd
    @test axes(ds,2) == 0:nt
    @test axes(ds,3) == 1:ni
    @test axes(ds,4) == 1:1
    @test size(ds.d) == (nd, nt+1, ni)
    @test size(ds.d) == size(ds)
    @test ndims(ds)  == 3
    @test eltype(ds) == dt
    @test parent(ds) == ds.d

    for j in 1:ni
        for i in 1:nt
            ds[1,i,j] = j*i
        end
    end

    @test ds.d[1,1,1] == ds[1,0,1]
    @test ds.d[1:ds.nd,1,1] == ds[1:ds.nd,0,1]

    @test ds[1,end,1] == ds.d[1,end,1]
    @test ds[1,nt,1] == ds[1,end,1]
    @test ds[1:ds.nd,nt,1] == ds[1:ds.nd,end,1]

    @test ds.d[:,1,1] == ds[:,0,1]
    @test ds.d[:,1,:] == ds[:,0,:]

    reset!(ds)
    @test ds[1,0,1] == ds[1,end,1]

    tx = rand(nd)
    ds[:,0,1] .= tx
    @test ds.d[:,1,1] == tx

end


text:
class Node {
  int elem;
  Node next;

  void setElem (int e)  { elem = e; }
  void setNext (Node n) { next = n; }

  int  getElem () { return elem; }
  Node getNext () { return next; }

}

class IntQueue {
  Node front;
  Node rear;

  boolean isEmpty () { return front == (Node)null; }

  void insert (int x) {
    Node last = new Node;
    last.setElem(x);
    if (self.isEmpty())
      front = last;
    else
      rear.setNext(last);
    rear = last;
  }

  int first () { return front.getElem(); }

  void rmFirst () {
    front = front.getNext();
  }

  int size () {
      Node n = front;
      int res = 0;
      while (n != (Node)null) {
        n = n.getNext();
        res++;
      }
     return res;
  }
}

int f (int x) {
  return x*x + 3;
}

int main () {
  IntQueue q = new IntQueue;
  q.insert(f(3));
  q.insert(5);
  q.insert(7);
  printInt(q.first());
  q.rmFirst();
  printInt(q.size());
  return 0;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule METIS4_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("METIS4")
JLLWrappers.@generate_main_file("METIS4", UUID("40b5814e-7855-5c9f-99f7-a735ce3fdf8b"))
end  # module METIS4_jll


text:
using DataFrames, DataArrays

"""
    tabulate(data)

Constructs a `Table` from the input data, attempting to make a structure with
This may or may not make copies of data, depending if the ingester decides it would
be beneficial (e.g. accelerate `getindex` from *O(N)* to *O(1)* or *O(log N)*).

Some examples include

 * `DataFrame` becomes a `Table` with column indices as `Label`s (`:x` becomes `l"x"`)
    and type-stable vectors as columns. Makes copies of columns which are `Nullable{T}`
    if the column contains `NA`, or else simply `T` if not.
"""
function tabulate(df::DataFrame)
    data = ()
    for n ∈ names(df)
        data = (data..., Label{n}() => make_faster(df[n]))
    end
    return Table(data)
end

function tabulate(v::AbstractVector{T}) where T
    Table(make_faster(make_faster.(v)))
end

function tabulate(d::Associative{String, <:Any})
    data = ()
    for n ∈ keys(d)
        data = (data..., Label{Symbol(n)}() => make_faster(d[n]))
    end
    return Table(data)
end

function tabulate(d::Associative{Symbol, <:Any})
    data = ()
    for n ∈ keys(d)
        data = (data..., Label{n}() => make_faster(d[n]))
    end
    return Table(data)
end

function make_faster(d::Associative{String, Any})
    data = ()
    for n ∈ keys(d)
        data = (data..., Label{Symbol(n)}() => d[n])
    end
    return data
end

function make_faster(d::Associative{Symbol, Any})
    data = ()
    for n ∈ keys(d)
        data = (data..., Label{n}() => d[n])
    end
    return data
end

function make_faster(vec::DataVector{T}) where {T}
    has_nulls = false
    for x ∈ vec
        if x === NA
            has_nulls = true
            break
        end
    end
    if has_nulls
        out = Vector{Nullable{T}}(length(vec))
        for i ∈ 1:length(vec)
            if vec[i] === NA
                out[i] = Nullable{T}()
            else
                out[i] = Nullable{T}(vec[i])
            end
        end
    else
        out = Vector{T}(length(vec))
        out[:] = vec[:]
    end
    return out
end

function make_faster(vec::AbstractVector{T}) where {T}
    if isleaftype(T)
        return vec
    end

    T2 = Union{}
    for x ∈ vec
        T2 = promote_type(T2, typeof(x))
    end

    out = similar(vec, T2)
    out[:] = vec[:]

    return out
end


text: @enum SizeUnit::Int8 SIZE_ABSOLUTE SIZE_SWAPCHAIN_RELATIVE SIZE_VIEWPORT_RELATIVE


text:
import Base: isequal, ==, show, *, length, >
using SymPy

"A plane geometric object."
abstract type GeoObject end
abstract type GeoShape <: GeoObject end

show(io::IO, objs::Vector{T}) where T <: GeoObject = print(io, "[", join(objs, ", "), "]")

"""
    Point(x, y)

A point at coordinate `(x, y)`.

# Examples
```julia
A = Point(0,0); B = Point(1, 3); C = Point(4,2)
```
"""
struct Point <: GeoObject
    "x cooridnate"
    x::Sym
    "y cooridnate"
    y::Sym

    Point(x, y) = new(Sym(x), Sym(y))
end

"""
# Function

    Point()

A point at coordinate `(0, 0)`.
"""
Point() = Point(0, 0)

"Check if two points are at the same coordinate."
(==)(p1::Point, p2::Point) = p1.x==p2.x && p1.y==p2.y

show(io::IO, pt::Point) = print(io, "Point($(pt.x), $(pt.y))")

"""
    Triangle(A, B, C)

A triangle with vertices `A, B, C`.

# Examples
```julia
A = Point(0,0); B = Point(1, 3); C = Point(4,2)
tri = Triangle(A, B, C)
```
"""
struct Triangle <: GeoShape
    "Vertex A"
    A::Point
    "Vertex B"
    B::Point
    "Vertex C"
    C::Point
end
show(io::IO, tri::Triangle) = print(io, "Triangle($(tri.A), $(tri.B), $(tri.C))")

"""
# Function

    Triangle()

A triangle with vertices `(0, 0), (0, 1), (1,1)`.
"""
Triangle() = Triangle(Point(0, 0), Point(0, 1), Point(1,1))

"""
    Triangle(ax, ay, bx, by, cx, cy)

A triangle with vertices `(ax, ay), (bx, by), (cx, cy)`.

# Examples
```julia
tri = Triangle(0, 0, 1, 3, 4, 2)
```
"""
Triangle(ax, ay, bx, by, cx, cy) = Triangle(Point(ax, ay), Point(bx, by), Point(cx, cy))

"Check if two triangles are the same."
(==)(t1::Triangle, t2::Triangle) = vertices(t1) == vertices(t2)

"Get the list of vertices of the triangle `tri`."
function vertices(tri::Triangle)
    [tri.A, tri.B, tri.C]
end

"""
    Edge(src, dst)

An edge from point `src` to point `dst`.

# Examples
```julia
e = Edge(Point(0, 0), Point(1, 3))
```
"""
struct Edge <: GeoShape
    "Starting point of an edge."
    src::Point
    "Ending point of an edge."
    dst::Point
end

"Check if two edges are at the same."
(==)(e1::Edge, e2::Edge) = e1.src==e2.src && e1.dst==e2.dst

"Get the list of edges of the triangle `tri`."
function edges(tri::Triangle)
    elist = Edge[]
    pts = vertices(tri)
    for i in 1:length(pts)-1
        push!(elist, Edge(pts[i], pts[i+1]))
    end
    push!(elist, Edge(pts[length(pts)], pts[1]))
    elist
end

"Mulitpliy the lengths of two edges."
(*)(e1::Edge, e2::Edge) = length(e1) * length(e2)

"Mulitpliy a number and an edges length."
(*)(s1::Sym, e2::Edge) = s1 * length(e2)

"Compare lengths of two edges."
(>)(e1::Edge, e2::Edge) = length(e1) > length(e2)

"Find the length of an edge."
length(e::Edge) = distance(e.src, e.dst)

"""
    Circle(c, r)

A circle with center at `c` and radius `r`.

# Examples
```julia
circle = Circle(Point(0, 0), 1)
```
"""
struct Circle <: GeoShape
    center::Point
    radius::Number
end
show(io::IO, c::Circle) = print(io, "Circle($(c.center), $(c.radius))")

"Check if two circles are the same."
(==)(c1::Circle, c2::Circle) = c1.center == c2.center && c1.radius == c2.radius

"Get the center of a a circle."
center(c::Circle) = c.center


text:
"""
    abstract type DiscreteModel{Dc,Dp} <: Grid

Abstract type holding information about a physical grid,
the underlying grid topology, and a labeling of
the grid faces. This is the information that typically provides a mesh
generator, and it is what one needs to perform a simulation.

The `DiscreteModel` interface is defined by overloading the methods:

- [`get_grid(model::DiscreteModel)`](@ref)
- [`get_grid_topology(model::DiscreteModel)`](@ref)
- [`get_face_labeling(g::DiscreteModel)`](@ref)

The interface is tested with this function:

- [`test_discrete_model`](@ref)

"""
abstract type DiscreteModel{Dc,Dp} <: Grid{Dc,Dp} end

"""
    get_grid(model::DiscreteModel)
"""
function get_grid(model::DiscreteModel)
  @abstractmethod
end

"""
    get_grid_topology(model::DiscreteModel)
"""
function get_grid_topology(model::DiscreteModel)
  @abstractmethod
end

"""
    get_face_labeling(g::DiscreteModel)
"""
function get_face_labeling(g::DiscreteModel)
  @abstractmethod
end

# Testers

"""
    test_discrete_model(model::DiscreteModel)
"""
function test_discrete_model(model::DiscreteModel{Dc,Dp}) where {Dc,Dp}
  D = Dc
  grid = get_grid(model)
  test_grid(grid)
  test_grid(model)
  topo = get_grid_topology(model)
  test_grid_topology(topo)
  @test num_cell_dims(topo) == Dc
  @test num_point_dims(topo) == Dp
  labels = get_face_labeling(model)
  @test num_faces(labels) == num_faces(topo)
  for d in 0:D
    @test num_faces(labels,d) == num_faces(topo,d)
  end
end

# Delegators to the underlying grid

get_cell_node_ids(g::DiscreteModel) = get_cell_node_ids(get_grid(g))

get_node_coordinates(g::DiscreteModel) = get_node_coordinates(get_grid(g))

get_cell_type(g::DiscreteModel) = get_cell_type(get_grid(g))

get_reffes(g::DiscreteModel) = get_reffes(get_grid(g))

# Default API

"""
    num_dims(model::DiscreteModel)
"""
num_dims(model::DiscreteModel) = num_dims(get_grid_topology(model))

"""
    num_cell_dims(model::DiscreteModel)
"""
num_cell_dims(model::DiscreteModel) = num_cell_dims(get_grid_topology(model))

"""
    num_point_dims(model::DiscreteModel)
"""
num_point_dims(model::DiscreteModel) = num_point_dims(get_grid_topology(model))

"""
    num_faces(g::DiscreteModel,d::Integer)
    num_faces(g::DiscreteModel)
"""
num_faces(g::DiscreteModel,d::Integer) = num_faces(get_grid_topology(g),d)
num_faces(g::DiscreteModel) = num_faces(get_grid_topology(g))

"""
    num_cells(g::DiscreteModel)
"""
num_cells(g::DiscreteModel) = num_cells(get_grid_topology(g))

"""
    num_facets(g::DiscreteModel)
"""
num_facets(g::DiscreteModel) = num_facets(get_grid_topology(g))

"""
    num_edges(g::DiscreteModel)
"""
num_edges(g::DiscreteModel) = num_edges(get_grid_topology(g))

"""
    num_vertices(g::DiscreteModel)
"""
num_vertices(g::DiscreteModel) = num_vertices(get_grid_topology(g))

"""
    num_nodes(g::DiscreteModel)
"""
num_nodes(g::DiscreteModel) = num_nodes(get_grid(g))

"""
    get_polytopes(model::DiscreteModel)
"""
function get_polytopes(model::DiscreteModel)
  topo = get_grid_topology(model)
  get_polytopes(topo)
end

"""
    get_face_nodes(g::DiscreteModel,d::Integer)
"""
function get_face_nodes(g::DiscreteModel,d::Integer)
  compute_face_nodes(g,d)
end

"""
    get_face_nodes(g::DiscreteModel)
"""
function get_face_nodes(g::DiscreteModel)
  compute_face_nodes(g)
end

"""
    compute_face_nodes(model::DiscreteModel,d::Integer)
"""
function compute_face_nodes(model::DiscreteModel,d::Integer)

  if d == num_cell_dims(model)
    return get_cell_node_ids(model)
  end

  topo = get_grid_topology(model)
  D = num_cell_dims(topo)
  cell_to_nodes = Table(get_cell_node_ids(model))
  cell_to_faces = Table(get_faces(topo,D,d))
  cell_to_ctype = get_cell_type(model)
  reffes = get_reffes(model)
  ctype_to_lface_to_lnodes = map( (reffe)-> get_face_nodes(reffe,d) , reffes )
  nfaces = num_faces(topo,d)

  face_to_nodes = generate_face_to_vertices(
    cell_to_nodes,
    cell_to_faces,
    cell_to_ctype,
    ctype_to_lface_to_lnodes,
    nfaces)

  face_to_nodes

end

"""
    compute_face_nodes(model::DiscreteModel)
"""
function compute_face_nodes(model::DiscreteModel)
  D = num_cell_dims(model)
  data = [compute_face_nodes(model,d) for d in 0:D]
  append_tables_globally(data...)
end

"""
    get_face_own_nodes(g::DiscreteModel,d::Integer)
"""
function get_face_own_nodes(g::DiscreteModel,d::Integer)
  compute_face_own_nodes(g,d)
end

"""
    get_face_own_nodes(g::DiscreteModel)
"""
function get_face_own_nodes(g::DiscreteModel)
  compute_face_own_nodes(g)
end

"""
    compute_face_own_nodes(model::DiscreteModel,d::Integer)
"""
function compute_face_own_nodes(model::DiscreteModel,d::Integer)

  topo = get_grid_topology(model)
  D = num_cell_dims(topo)
  cell_to_nodes = Table(get_cell_node_ids(model))
  cell_to_faces = Table(get_faces(topo,D,d))
  cell_to_ctype = get_cell_type(model)
  reffes = get_reffes(model)
  ctype_to_lface_to_lnodes = map( (reffe)-> get_face_own_nodes(reffe,d) , reffes )
  nfaces = num_faces(topo,d)

  face_to_own_nodes = generate_face_to_vertices(
    cell_to_nodes,
    cell_to_faces,
    cell_to_ctype,
    ctype_to_lface_to_lnodes,
    nfaces)

  face_to_own_nodes

end

"""
    compute_face_own_nodes(model::DiscreteModel)
"""
function compute_face_own_nodes(model::DiscreteModel)
  D = num_cell_dims(model)
  data = [compute_face_own_nodes(model,d) for d in 0:D]
  append_tables_globally(data...)
end

"""
    get_vertex_node(g::DiscreteModel)
"""
function get_vertex_node(g::DiscreteModel)
  compute_vertex_node(g)
end

"""
    compute_vertex_node(g::DiscreteModel)
"""
function compute_vertex_node(g::DiscreteModel)
  d=0
  vertex_to_nodes = Table(get_face_own_nodes(g,d))
  vertex_to_nodes.data
end

"""
    get_node_face_owner(g::DiscreteModel)
"""
function get_node_face_owner(g::DiscreteModel)
  compute_node_face_owner(g)
end

"""
    compute_node_face_owner(g::DiscreteModel)
"""
function compute_node_face_owner(g::DiscreteModel)
  face_to_own_nodes = Table(get_face_own_nodes(g))
  node_to_face_owner = zeros(Int32,num_nodes(g))
  _compute_node_face_owner!(node_to_face_owner,face_to_own_nodes)
  node_to_face_owner
end

function  _compute_node_face_owner!(node_to_face_owner,face_to_own_nodes)
  for face in 1:length(face_to_own_nodes)
    pini = face_to_own_nodes.ptrs[face]
    pend = face_to_own_nodes.ptrs[face+1]-1
    for p in pini:pend
      node = face_to_own_nodes.data[p]
      node_to_face_owner[node] = face
    end
  end
end

"""
    get_reffaces(::Type{ReferenceFE{d}},model::DiscreteModel) where d
"""
function get_reffaces(::Type{ReferenceFE{d}},model::DiscreteModel) where d
  reffaces,_ = compute_reffaces(ReferenceFE{d},model)
  reffaces
end

"""
    get_face_type(g::DiscreteModel,d::Integer)

Index to the vector `get_reffaces(ReferenceFE{d},g)`
"""
function get_face_type(g::DiscreteModel,d::Integer)
  _, face_to_ftype = compute_reffaces(ReferenceFE{d},g)
  face_to_ftype
end

"""
    compute_reffaces(::Type{ReferenceFE{d}}, g::DiscreteModel) where d
"""
function compute_reffaces(::Type{ReferenceFE{d}}, g::DiscreteModel) where d
  D = num_cell_dims(g)
  topo = get_grid_topology(g)
  ctype_to_reffe = get_reffes(g)
  ctype_to_lftype_to_refface = [ get_reffaces(ReferenceFE{d},reffe) for reffe in ctype_to_reffe]
  ctype_to_lface_to_lftype = [ get_face_type(reffe,d) for reffe in ctype_to_reffe]
  t = _generate_ftype_to_refface(Val{d}(),ctype_to_lftype_to_refface,ctype_to_lface_to_lftype)
  ftype_to_refface, ctype_to_lface_to_ftype = t
  cell_to_faces = Table(get_faces(topo,D,d))
  cell_to_ctype = get_cell_type(g)
  nfaces = num_faces(g,d)
  face_to_ftype = generate_face_to_face_type(
    cell_to_faces, cell_to_ctype, ctype_to_lface_to_ftype, nfaces)
  (collect1d(ftype_to_refface), face_to_ftype)
end

function compute_reffaces(::Type{ReferenceFE{D}}, g::DiscreteModel{D}) where D
  (get_reffes(g), get_cell_type(g))
end

"""
    get_reffaces(model::DiscreteModel)
"""
function get_reffaces(model::DiscreteModel)
  reffaces, _ , _ = compute_reffaces(model)
  reffaces
end

"""
    get_face_type(model::DiscreteModel)
"""
function get_face_type(model::DiscreteModel)
  _, face_to_ftype, _ = compute_reffaces(model)
  face_to_ftype
end

"""
    get_reffaces_offsets(model::DiscreteModel)
"""
function get_reffaces_offsets(model::DiscreteModel)
  _, _, offsets = compute_reffaces(model)
  offsets
end

"""
    compute_reffaces(g::DiscreteModel)
"""
function compute_reffaces(g::DiscreteModel)
  D = num_cell_dims(g)
  d_to_refdfaces = Vector{LagrangianRefFE}[]
  d_to_dface_to_ftype = Vector{Int8}[]
  for d in 0:D
    push!(d_to_refdfaces,get_reffaces(ReferenceFE{d},g))
    push!(d_to_dface_to_ftype,get_face_type(g,d))
  end
  d_to_offset = zeros(Int,D+1)
  for d in 1:D
    d_to_offset[d+1] = d_to_offset[d] + length(d_to_refdfaces[d])
    d_to_dface_to_ftype[d+1] .+= d_to_offset[d+1]
  end
  (vcat(d_to_refdfaces...), vcat(d_to_dface_to_ftype...), d_to_offset)
end

"""
    Grid(::Type{ReferenceFE{d}},model::DiscreteModel) where d
"""
function Grid(::Type{ReferenceFE{d}},model::DiscreteModel) where d
  node_coordinates = collect1d(get_node_coordinates(model))
  cell_to_nodes = Table(get_face_nodes(model,d))
  cell_to_type = collect1d(get_face_type(model,d))
  reffes = get_reffaces(ReferenceFE{d},model)
  UnstructuredGrid(node_coordinates, cell_to_nodes, reffes, cell_to_type)
end

function Grid(::Type{ReferenceFE{d}},model::DiscreteModel{d}) where d
  get_grid(model)
end

"""
    simplexify(model::DiscreteModel)
"""
function simplexify(model::DiscreteModel)
  umodel = UnstructuredDiscreteModel(model)
  simplexify(umodel)
end

function ReferenceFE(model::DiscreteModel,args...;kwargs...)
  ctype_to_polytope = get_polytopes(model)
  cell_to_ctype = get_cell_type(model)
  ctype_to_reffe = map(p->ReferenceFE(p,args...;kwargs...),ctype_to_polytope)
  cell_to_reffe = expand_cell_data(ctype_to_reffe,cell_to_ctype)
  cell_to_reffe
end

function ReferenceFE(model::DiscreteModel,basis::ModalC0,args...;kwargs...)
  ctype_to_polytope = get_polytopes(model)
  @assert length(ctype_to_polytope) == 1 "Only one polytope expected"
  compute_cell_to_modalC0_reffe(ctype_to_polytope[1],num_cells(model),args...;kwargs...)
end

# IO

function to_dict(model::DiscreteModel)
  umodel = UnstructuredDiscreteModel(model)
  to_dict(umodel)
end

function from_dict(::Type{DiscreteModel},dict::Dict{Symbol,Any})
  from_dict(UnstructuredDiscreteModel,dict)
end

"""
    DiscreteModelFromFile(filename::AbstractString)
"""
function DiscreteModelFromFile(filename::AbstractString)
  base, extension = splitext(filename)
  s = Symbol(extension[2:end])
  DiscreteModelFromFile(filename,Val(s))
end

function DiscreteModelFromFile(filename::AbstractString,::Any)
  @notimplemented
end

function DiscreteModelFromFile(filename::AbstractString,::Val{:json})
  model = from_json_file(DiscreteModel,filename)
  model
end

# GenericDiscreteModel

struct GenericDiscreteModel{Dc,Dp} <: DiscreteModel{Dc,Dp}
  grid::Grid
  grid_topology::GridTopology
  labels::FaceLabeling
  function GenericDiscreteModel(grid::Grid{Dc,Dp},grid_topology::GridTopology,labels::FaceLabeling) where {Dc,Dp}
    new{Dc,Dp}(grid,grid_topology,labels)
  end
end

"""
"""
function DiscreteModel(grid::Grid,grid_topology::GridTopology,labels::FaceLabeling)
  GenericDiscreteModel(grid,grid_topology,labels)
end

get_grid(model::GenericDiscreteModel) = model.grid

get_grid_topology(model::GenericDiscreteModel) = model.grid_topology

get_face_labeling(model::GenericDiscreteModel) = model.labels

function DiscreteModel(::Type{<:Polytope{D}},model::DiscreteModel{D}) where D
  model
end

"""
"""
function DiscreteModel(::Type{<:Polytope{D}},model::DiscreteModel) where D
  grid = Grid(ReferenceFE{D},model)
  topo = GridTopology(Polytope{D},model)
  labeling = FaceLabeling(get_face_labeling(model),D)
  DiscreteModel(grid,topo,labeling)
end

"""
"""
function GridTopology(::Type{<:Polytope{D}},model::DiscreteModel) where D
  GridTopology(Polytope{D},get_grid_topology(model))
end


text:
#
# Copyright (c) Chernobylll

# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMI

# our simulation setup
dt = 0.01
t_start = 0.0
t_stop = 8.0
saveat = t_start:dt:t_stop

# this FMU runs under Windows/Linux
pathToFMU = joinpath(dirname(@__FILE__), "../model/OpenModelica/v1.17.0/SpringFrictionPendulum1D.fmu")

# this FMU runs only under Windows
if Sys.iswindows()
    pathToFMU = joinpath(dirname(@__FILE__), "../model/Dymola/2020x/SpringFrictionPendulum1D.fmu")
end

# load the FMU container
myFMU = fmiLoad(pathToFMU)

# print some useful FMU-information into the REPL
fmiInfo(myFMU)

# make an instance from the FMU
fmiInstantiate!(myFMU; loggingOn=true)

# setup the experiment, start time = 0.0 (optional for setup=true)
#fmiSetupExperiment(myFMU, t_start)

# enter and exit initialization (optional for setup=true)
#fmiEnterInitializationMode(myFMU)
#fmiExitInitializationMode(myFMU)

# run the FMU in mode Co-Simulation (CS) with adaptive step size (CVODE) but fixed save points dt,
# result values are stored in `savedValues`
recordValues = ["mass.s", "mass.v"]
success, savedValues = fmiSimulateCS(myFMU, t_start, t_stop; recordValues=recordValues, setup=true)

# plot the results
fmiPlot(myFMU, recordValues, savedValues)

# unload the FMU, remove unpacked data on disc ("clean up")
fmiUnload(myFMU)

max_stars_repo_name: Chernobylll/fmi-codes-and-variants

max_issues_repo_name: Chernobylll/fmi-codes-and-variants

max_forks_repo_name: Chernobylll/fmi-codes-and-variants


text:
function insertrng!(m::Module)
    id = 0
    # replace all SBitstream constants and re-route parameters
    for v in vertices(m.dfg)
        inputs = getinputs(m.dfg, v)
        for (i, input) in enumerate(inputs)
            # replace constant
            if isconstant(input) && (jltypeof(input) <: SBitstream)
                val = Net(float(value(input));
                          name = string(float(value(input))),
                          class = :constant)
                op = (name = Symbol(""), type = typeof(SBitstream), broadcasted = false)
                input = setclass(input, :internal)
                input = setname(input, "net_rng_$id")
                id += 1
                inputs[i] = input
                set_prop!(m.dfg, v, :inputs, inputs)
                addnode!(m, [val], [input], op)
            elseif isconstant(input) && (jltypeof(input) <: AbstractArray{<:SBitstream})
                val = Net(float.(value(input));
                          name = string.(float.(value(input))),
                          class = :constant)
                op = (name = Symbol(""), type = typeof(SBitstream), broadcasted = true)
                input = setclass(input, :internal)
                input = setname(input, "net_rng_$id")
                id += 1
                inputs[i] = input
                set_prop!(m.dfg, v, :inputs, inputs)
                addnode!(m, [val], [input], op)
            end

            # re-route parameter
            # the name is fixed in the second loop below
            if isparameter(input) &&
               (jltypeof(input) <: Union{SBitstream, AbstractArray{<:SBitstream}})
                newnet = Net(value(input);
                             name = "$(name(input))_rng", class = :internal, signed = true)
                inputs[i] = newnet
                set_prop!(m.dfg, v, :inputs, inputs)
            end
        end
    end

    # replace all SBitstream parameters
    for (name, value) in m.parameters
        if typeof(value) <: SBitstream
            rngin = Net(float(value); name = name, class = :parameter)
            rngout = Net(value; name = "$(name)_rng", class = :internal, signed = true)
            op = (name = Symbol(""), type = typeof(SBitstream), broadcasted = false)
            m.parameters[name] = float(value)
            addnode!(m, [rngin], [rngout], op)
        elseif typeof(value) <: AbstractArray{<:SBitstream}
            rngin = Net(float.(value); name = name, class = :parameter)
            rngout = Net(value; name = "$(name)_rng", class = :internal)
            op = (name = Symbol(""), type = typeof(SBitstream), broadcasted = true)
            m.parameters[name] = float.(value)
            addnode!(m, [rngin], [rngout], op)
        end
    end

    return m
end