text:
module TestExamples3DMHD

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_3d_dgsem")

@testset "MHD" begin
  @trixi_testset "elixir_mhd_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
      l2   = [0.017590099293094203, 0.017695875823827714, 0.017695875823827686, 0.017698038279620777, 0.07495006099352074, 0.010391801950005755, 0.010391801950005759, 0.010393502246627087, 2.524766553484067e-16],
      linf = [0.28173002819718196, 0.3297583616136297, 0.32975836161363004, 0.356862935505337, 1.2893514981209626, 0.10950981489747313, 0.10950981489747136, 0.11517234329681891, 2.0816911067714202e-15])
  end

  @trixi_testset "elixir_mhd_ec.jl with initial_condition=initial_condition_constant" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
      l2   = [4.270231310667203e-16, 2.4381208042014784e-15, 5.345107673575357e-15, 3.00313882171883e-15, 1.7772703118758417e-14, 1.0340110783830874e-15, 1.1779095371939702e-15, 9.961878521814573e-16, 8.1201730630719145e-16],
      linf = [2.4424906541753444e-15, 2.881028748902281e-14, 2.4646951146678475e-14, 2.3092638912203256e-14, 2.3447910280083306e-13, 1.7763568394002505e-14, 1.0436096431476471e-14, 2.042810365310288e-14, 7.057203733035201e-15],
      atol = 1000*eps(),
      initial_condition=initial_condition_constant)
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [0.0032217291057246157, 0.009511644936958913, 0.004217358459420256, 0.011591709179125335, 0.009456218722393708, 0.00916500047763897, 0.005069863732625444, 0.011503011541926135, 0.003988175543749985],
      linf = [0.01188593784273051, 0.03638015998373141, 0.01568200398945724, 0.04666974730787579, 0.031235294705421968, 0.03316343064943483, 0.011539436992528018, 0.04896687646520839, 0.018714054039927555])
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl with flux_derigs_etal" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [0.003755235939722358, 0.009062519246840721, 0.004096299856228109, 0.011429935838448906, 0.006897420817511043, 0.00900886245212482, 0.004926537542780259, 0.01153285554590683, 0.0037842060148666886],
      linf = [0.012982853115883541, 0.0320228076558316, 0.011575276754611022, 0.04425778643430531, 0.02478109022285846, 0.03198699034954189, 0.009761077061886558, 0.04433669321441455, 0.01618905441148782],
      volume_flux = (flux_derigs_etal, flux_nonconservative_powell))
  end

  @trixi_testset "elixir_mhd_alfven_wave_mortar.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave_mortar.jl"),
      l2   = [0.001879021634926363, 0.007032724521848316, 0.0032793932234187325, 0.009056594733320348, 0.007514150120617965, 0.007328739509868727, 0.00309794018112387, 0.009026356949274878, 0.0035732583778049776],
      linf = [0.013734346970999622, 0.06173467158736011, 0.02183946452704291, 0.06258216169457917, 0.03672304497348122, 0.055120532123884625, 0.018202716205672487, 0.06133688282205586, 0.019888161885935608],
      tspan = (0.0, 0.25))
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl with Orszag-Tang setup + flux_hll" begin
    # OBS! This setup does not make much sense and is only used to exercise all components of the
    # flux_hll implementation
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [0.004391143689111404, 0.04144737547475548, 0.041501307637678286, 0.04150353006408862, 0.03693135855995625, 0.021125605214031118, 0.03295607553556973, 0.03296235755245784, 7.16035229384135e-6],
      linf = [0.017894703320895378, 0.08486850681397005, 0.0891044523165206, 0.08492024792056754, 0.10448301878352373, 0.05381260695579509, 0.0884774018719996, 0.07784546966765199, 7.71609149516089e-5],
      initial_condition = function initial_condition_orszag_tang(x, t, equations::IdealGlmMhdEquations3D)
        # The classical Orszag-Tang vortex test case adapted to 3D. Setup is taken from
        # Table 4 of the paper
        # - M. Bohm, A. R. Winters, G. J. Gassner, D. Derigs, F. Hindenlang, & J. Saur (2020)
        #   An entropy stable nodal discontinuous Galerkin method for the resistive MHD
        #   equations. Part I: Theory and numerical verification
        #   [doi: 10.1016/j.jcp.2018.06.027](https://doi.org/10.1016/j.jcp.2018.06.027)
        # Domain must be [0, 1]^3 , γ = 5/3
        rho = 25.0 / (36.0 * pi)
        v1 = -sin(2.0*pi*x[3])
        v2 =  sin(2.0*pi*x[1])
        v3 =  sin(2.0*pi*x[2])
        p = 5.0 / (12.0 * pi)
        B1 = -sin(2.0*pi*x[3]) / (4.0*pi)
        B2 =  sin(4.0*pi*x[1]) / (4.0*pi)
        B3 =  sin(4.0*pi*x[2]) / (4.0*pi)
        psi = 0.0
        return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations)
      end,
      surface_flux = (flux_hll, flux_nonconservative_powell),
      volume_flux  = (flux_central, flux_nonconservative_powell),
      coordinates_min = (0.0, 0.0, 0.0),
      coordinates_max = (1.0, 1.0, 1.0),
      initial_refinement_level=3,
      cfl = 1.1,
      tspan = (0.0, 0.06))
  end

  @trixi_testset "elixir_mhd_ec_shockcapturing.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec_shockcapturing.jl"),
      l2   = [0.0186712969755079, 0.01620736832264799, 0.01620736832264803, 0.016207474382769683, 0.07306422729650594, 0.007355137041002365, 0.0073551370410023425, 0.00735520932001833, 0.000506140942330923],
      linf = [0.28040713666979633, 0.27212885844703694, 0.2721288584470349, 0.2837380205051839, 0.7915852408267114, 0.08770240288089526, 0.08770240288089792, 0.08773409387876674, 0.050221095224119834])
  end
end

end # module


text:
using Manifolds
using BenchmarkTools
using StaticArrays
using LinearAlgebra
using Random

# Define a parent BenchmarkGroup to contain our suite
SUITE = BenchmarkGroup()

Random.seed!(12334)

function add_manifold(M::Manifold, pts, name;
    test_tangent_vector_broadcasting = true,
    retraction_methods = [],
    inverse_retraction_methods = [],
    point_distributions = [],
    tvector_distributions = [])

    SUITE["manifolds"][name] = BenchmarkGroup()
    tv = log(M, pts[1], pts[2])
    tv1 = log(M, pts[1], pts[2])
    tv2 = log(M, pts[1], pts[3])
    p = similar(pts[1])
    SUITE["manifolds"][name]["similar"] = @benchmarkable similar($(pts[1]))
    SUITE["manifolds"][name]["log"] = @benchmarkable log($M, $(pts[1]), $(pts[2]))
    SUITE["manifolds"][name]["log!"] = @benchmarkable log!($M, $tv, $(pts[1]), $(pts[2]))
    for iretr ∈ inverse_retraction_methods
        SUITE["manifolds"][name]["inverse_retract: "*string(iretr)] = @benchmarkable inverse_retract($M, $(pts[1]), $(pts[2]), $iretr)
        SUITE["manifolds"][name]["inverse_retract!: "*string(iretr)] = @benchmarkable inverse_retract!($M, $tv, $(pts[1]), $(pts[2]), $iretr)
    end
    SUITE["manifolds"][name]["exp"] = @benchmarkable exp($M, $(pts[1]), $tv1)
    SUITE["manifolds"][name]["exp!"] = @benchmarkable exp!($M, $p, $(pts[1]), $tv1)
    for retr ∈ retraction_methods
        SUITE["manifolds"][name]["retract: "*string(retr)] = @benchmarkable retract($M, $(pts[1]), $tv1, $retr)
        SUITE["manifolds"][name]["retract!: "*string(retr)] = @benchmarkable retract!($M, $p, $(pts[1]), $tv1, $retr)
    end
    SUITE["manifolds"][name]["norm"] = @benchmarkable norm($M, $(pts[1]), $tv1)
    SUITE["manifolds"][name]["inner"] = @benchmarkable inner($M, $(pts[1]), $tv1, $tv2)
    SUITE["manifolds"][name]["distance"] = @benchmarkable distance($M, $(pts[1]), $(pts[2]))
    SUITE["manifolds"][name]["isapprox (pt)"] = @benchmarkable isapprox($M, $(pts[1]), $(pts[2]))
    SUITE["manifolds"][name]["isapprox (tv)"] = @benchmarkable isapprox($M, $(pts[1]), $tv1, $tv2)
    SUITE["manifolds"][name]["2 * tv1 + 3 * tv2"] = @benchmarkable 2 * $tv1 + 3 * $tv2
    SUITE["manifolds"][name]["tv = 2 * tv1 + 3 * tv2"] = @benchmarkable $tv = 2 * $tv1 + 3 * $tv2
    if test_tangent_vector_broadcasting
        SUITE["manifolds"][name]["tv = 2 .* tv1 .+ 3 .* tv2"] = @benchmarkable $tv = 2 .* $tv1 .+ 3 .* $tv2
        SUITE["manifolds"][name]["tv .= 2 .* tv1 .+ 3 .* tv2"] = @benchmarkable $tv .= 2 .* $tv1 .+ 3 .* $tv2
    end
    for pd ∈ point_distributions
        distr_name = string(pd)
        distr_name = distr_name[1:min(length(distr_name), 50)]
        SUITE["manifolds"][name]["point distribution "*distr_name] = @benchmarkable rand($pd)
    end
    for tvd ∈ point_distributions
        distr_name = string(tvd)
        distr_name = distr_name[1:min(length(distr_name), 50)]
        SUITE["manifolds"][name]["tangent vector distribution "*distr_name] = @benchmarkable rand($tvd)
    end
end

# General manifold benchmarks
function add_manifold_benchmarks()

    SUITE["manifolds"] = BenchmarkGroup()

    s2 = Manifolds.Sphere(2)
    array_s2 = ArrayManifold(s2)
    r2 = Manifolds.Euclidean(2)

    pts_r2 = [Size(2)([1.0, 1.0]),
              Size(2)([-2.0, 3.0]),
              Size(2)([3.0, -2.0])]

    add_manifold(r2,
                 pts_r2,
                 "Euclidean{2} -- SizedArray")

    add_manifold(r2,
                 [MVector{2,Float64}([1.0, 1.0]),
                  MVector{2,Float64}([-2.0, 3.0]),
                  MVector{2,Float64}([3.0, -2.0])],
                  "Euclidean{2} -- MVector")

    ud_sphere = Manifolds.uniform_distribution(s2, Size(3)([1.0, 0.0, 0.0]))
    gtd_sphere = Manifolds.normal_tvector_distribution(s2, Size(3)([1.0, 0.0, 0.0]), 1.0)

    pts_s2 = [Size(3)([1.0, 0.0, 0.0]),
              Size(3)([0.0, 1.0, 0.0]),
              Size(3)([0.0, 0.0, 1.0])]

    add_manifold(s2,
                 pts_s2,
                 "Sphere{2} -- SizedArray";
                 point_distributions = [ud_sphere],
                 tvector_distributions = [gtd_sphere])

    add_manifold(array_s2,
                 [Size(3)([1.0, 0.0, 0.0]),
                  Size(3)([0.0, 1.0, 0.0]),
                  Size(3)([0.0, 0.0, 1.0])],
                 "ArrayManifold{Sphere{2}} -- SizedArray";
                 test_tangent_vector_broadcasting = false)

    retraction_methods_rot = [Manifolds.PolarRetraction(),
                              Manifolds.QRRetraction()]

    inverse_retraction_methods_rot = [Manifolds.PolarInverseRetraction(),
                                      Manifolds.QRInverseRetraction()]

    so2 = Manifolds.Rotations(2)
    angles = (0.0, π/2, 2π/3)
    add_manifold(so2,
                 [Size(2, 2)([cos(ϕ) -sin(ϕ); sin(ϕ) cos(ϕ)]) for ϕ in angles],
                 "Rotations(2) -- SizedArray",
                 retraction_methods = retraction_methods_rot,
                 inverse_retraction_methods = inverse_retraction_methods_rot)

    so3 = Manifolds.Rotations(3)
    angles = (0.0, π/2, 2π/3)
    add_manifold(so3,
                 [Size(3, 3)([cos(ϕ) -sin(ϕ) 0; sin(ϕ) cos(ϕ) 0; 0 0 1]) for ϕ in angles],
                 "Rotations(3) -- SizedArray",
                 retraction_methods = retraction_methods_rot,
                 inverse_retraction_methods = inverse_retraction_methods_rot)

    m_prod = Manifolds.ProductManifold(s2, r2)
    shape_s2r2_array = Manifolds.ShapeSpecification(Manifolds.ArrayReshaper(), s2, r2)
    shape_s2r2_static = Manifolds.ShapeSpecification(Manifolds.StaticReshaper(), s2, r2)

    pts_prd_base = [[1.0, 0.0, 0.0, 0.0, 0.0],
                    [0.0, 1.0, 0.0, 1.0, 0.0],
                    [0.0, 0.0, 1.0, 0.0, 0.1]]
    pts_prod_static = map(p -> Manifolds.ProductArray(shape_s2r2_static, p), pts_prd_base)
    add_manifold(m_prod, pts_prod_static, "ProductManifold with ProductArray (static)")

    pts_prod_array = map(p -> Manifolds.ProductArray(shape_s2r2_array, p), pts_prd_base)
    add_manifold(m_prod, pts_prod_array, "ProductManifold with ProductArray (array)")

    pts_prod_mpoints = [Manifolds.ProductMPoint(p[1], p[2]) for p in zip(pts_s2, pts_r2)]
    add_manifold(m_prod, pts_prod_mpoints, "ProductManifold with MPoint";
        test_tangent_vector_broadcasting = false)
end

add_manifold_benchmarks()

ext: jl


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

name = "libportaudio"
version = v"19.6.0"

# Collection of sources required to build libportaudio. Not all of these
# are used for all platforms.
sources = [
    "http://portaudio.com/archives/pa_stable_v190600_20161030.tgz" =>
        "f5a21d7dcd6ee84397446fa1fa1a0675bb2e8a4a6dceb4305a8404698d8d1513",

    "https://github.com/jackaudio/jack2/releases/download/v1.9.12/jack2-1.9.12.tar.gz" =>
        "deefe2f936dc32f59ad3cef7e37276c2035ef8a024ca92118f35c9a292272e33",

    "ftp://ftp.alsa-project.org/pub/lib/alsa-lib-1.1.6.tar.bz2" =>
        "5f2cd274b272cae0d0d111e8a9e363f08783329157e8dd68b3de0c096de6d724",

    "http://www.steinberg.net/sdk_downloads/ASIOSDK2.3.1.zip" =>
        "31074764475059448a9b7a56f103f4723ed60465e0e9d1a9446ca03dcf840f04"
        ]

# Bash recipe for building across all platforms
script = raw"""
set -e
cd $WORKSPACE/srcdir
# set up our platform-specific switches
case "$target" in
    *linux*)
        jack_platform=linux
        build_alsa=true
        build_jack=true
    ;;
    *freebsd*)
        jack_platform=linux
        build_alsa=true
        build_jack=true
    ;;
    *darwin*)
        jack_platform=darwin
        build_alsa=false
        build_jack=true
    ;;
    *mingw32*)
        jack_platform=msys
        build_alsa=false
        build_jack=true
    ;;
esac
if [[ "$build_alsa" == "true" ]]; then
    cd alsa-lib-1.1.6/
    ./configure --prefix=$prefix --host=$target
    make
    make install
    cd ..
fi
if [[ "$build_jack" == "true" ]]; then
    cd jack2-1.9.12/
    ./waf configure --prefix=$prefix --platform=$jack_platform
    ./waf build
    ./waf install
    cd ..
fi

# move the ASIO SDK to where CMake can find it
mv "asiosdk2.3.1 svnrev312937/ASIOSDK2.3.1" asiosdk2.3.1
rmdir "asiosdk2.3.1 svnrev312937"

mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DCMAKE_DISABLE_FIND_PACKAGE_PkgConfig=ON ../portaudio/
make
make install
"""

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

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libportaudio", :libportaudio)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://s3.us-east-2.amazonaws.com/ringbuffersbuilder/build_pa_ringbuffer.v19.6.0.jl"
]

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

max_stars_repo_name: ssfrr/PortAudioBuilder

max_issues_repo_name: ssfrr/PortAudioBuilder

max_forks_repo_name: ssfrr/PortAudioBuilder


text:
using PyCall
using DPMMSubClusters.TorchWrapper: torch_model_2d, torch_model_10d, torch_model_128d

function create_first_local_cluster(group::local_group)
    suff =
        create_sufficient_statistics(group.model_hyperparams.distribution_hyper_params, [])
    post = group.model_hyperparams.distribution_hyper_params
    dist = sample_distribution(post)
    cp = cluster_parameters(
        group.model_hyperparams.distribution_hyper_params,
        dist,
        suff,
        post,
    )
    cpl = deepcopy(cp)
    cpr = deepcopy(cp)
    splittable = splittable_cluster_params(
        cp,
        cpl,
        cpr,
        [0.5, 0.5],
        false,
        ones(burnout_period + 5) * -Inf,
    )
    cp.suff_statistics.N = size(group.points, 2)
    cpl.suff_statistics.N = sum(group.labels_subcluster .== 1)
    cpl.suff_statistics.N = sum(group.labels_subcluster .== 2)
    cluster = local_cluster(
        splittable,
        group.model_hyperparams.total_dim,
        cp.suff_statistics.N,
        cpl.suff_statistics.N,
        cpl.suff_statistics.N,
    )
    @sync for i in (nworkers() == 0 ? procs() : workers())
        @spawnat i split_first_cluster_worker!(group)
    end
    return cluster
end

function create_outlier_local_cluster(group::local_group, outlier_params)
    suff = create_sufficient_statistics(outlier_params, outlier_params, Array(group.points))
    post = calc_posterior(outlier_params, suff)
    dist = sample_distribution(post)
    cp = cluster_parameters(outlier_params, dist, suff, post)
    cpl = deepcopy(cp)
    cpr = deepcopy(cp)
    splittable = splittable_cluster_params(
        cp,
        cpl,
        cpr,
        [0.5, 0.5],
        false,
        ones(burnout_period + 5) * -Inf,
    )
    cp.suff_statistics.N = size(group.points, 2)
    cpl.suff_statistics.N = sum(group.labels_subcluster .== 1)
    cpl.suff_statistics.N = sum(group.labels_subcluster .== 2)
    cluster = local_cluster(
        splittable,
        group.model_hyperparams.total_dim,
        cp.suff_statistics.N,
        cpl.suff_statistics.N,
        cpl.suff_statistics.N,
    )
    # @sync for i in (nworkers()== 0 ? procs() : workers())
    #     @spawnat i split_first_cluster_worker!(group)
    # end
    return cluster
end

function sample_sub_clusters!(group::local_group)
    for i in (nworkers() == 0 ? procs() : workers())
        @spawnat i sample_sub_clusters_worker!(
            group.points,
            group.labels,
            group.labels_subcluster,
        )
    end
end

function sample_sub_clusters_worker!(group_points, group_labels, group_labels_subcluster)

    pts = localpart(group_points)
    labels = localpart(group_labels)
    sub_labels = localpart(group_labels_subcluster)

    for (i, v) in enumerate(clusters_vector)
        create_subclusters_labels!(
            (@view sub_labels[labels.==i]),
            (@view pts[:, labels.==i]),
            v,
        )
    end

end

function create_subclusters_labels!(
    labels::AbstractArray{Int64,1},
    points::AbstractArray{Float32,2},
    cluster_params::thin_cluster_params,
)
    if size(labels, 1) == 0
        return
    end
    parr = zeros(Float32, length(labels), 2)
    log_likelihood!((@view parr[:, 1]), points, cluster_params.l_dist)
    log_likelihood!((@view parr[:, 2]), points, cluster_params.r_dist)
    parr[:, 1] .+= log(cluster_params.lr_weights[1])
    parr[:, 2] .+= log(cluster_params.lr_weights[2])
    sample_log_cat_array!(labels, parr)
end

function sample_labels!(group::local_group, final::Bool, no_more_splits::Bool)
    sample_labels!(group.labels, group.points, final, no_more_splits)
end

function sample_labels!(
    labels::AbstractArray{Int64,1},
    points::AbstractArray{Float32,2},
    final::Bool,
    no_more_splits::Bool,
)
    for i in (nworkers() == 0 ? procs() : workers())
        @spawnat i sample_labels_worker!(labels, points, final, no_more_splits)
    end
end

function sample_labels_worker!(
    labels::AbstractArray{Int64,1},
    points::AbstractArray{Float32,2},
    final::Bool,
    no_more_splits::Bool,
)
    indices = localindices(points)[2]
    lbls = localpart(labels)
    pts = localpart(points)
    log_weights = log.(clusters_weights)
    parr = zeros(Float32, length(indices), length(clusters_vector))
    tic = time()
    for (k, cluster) in enumerate(clusters_vector)
        log_likelihood!(reshape((@view parr[:, k]), :, 1), pts, cluster.cluster_dist)
    end
    for (k, v) in enumerate(clusters_weights)
        parr[:, k] .+= log(v)
    end

    if final
        lbls .= mapslices(argmax, parr, dims = [2])[:]
    else
        sample_log_cat_array!(lbls, parr)
    end
end

function update_splittable_cluster_params!(splittable_cluser::splittable_cluster_params)
    cpl = splittable_cluser.cluster_params_l
    cpr = splittable_cluser.cluster_params_r
    cp = splittable_cluser.cluster_params

    cp.posterior_hyperparams = calc_posterior(cp.hyperparams, cp.suff_statistics)
    cpl.posterior_hyperparams = calc_posterior(cpl.hyperparams, cpl.suff_statistics)
    cpr.posterior_hyperparams = calc_posterior(cpr.hyperparams, cpr.suff_statistics)

end

function create_suff_stats_dict_worker(
    group_pts,
    group_labels,
    group_sublabels,
    hyper_params,
    indices,
)
    suff_stats_dict = Dict()
    if indices == nothing
        indices = collect(1:length(clusters_vector))
    end

    points = localpart(group_pts)
    labels = localpart(group_labels)
    sublabels = localpart(group_sublabels)
    for index in indices
        pts = @view points[:, labels.==index]
        sub_labels = @view sublabels[labels.==index]

        cpl_suff = create_sufficient_statistics(
            hyper_params,
            hyper_params,
            @view pts[:, sub_labels.==1]
        )
        cpr_suff = create_sufficient_statistics(
            hyper_params,
            hyper_params,
            @view pts[:, sub_labels.==2]
        )
        cp_suff = create_sufficient_statistics(hyper_params, hyper_params, pts)
        suff_stats_dict[index] = thin_suff_stats(cp_suff, cpl_suff, cpr_suff)
    end
    return suff_stats_dict
end

function create_suff_stats_dict_node_leader(
    group_pts,
    group_labels,
    group_sublabels,
    hyper_params,
    proc_ids,
    indices,
)
    leader_suff_dict = Dict()
    workers_suff_dict = Dict()
    if indices == nothing
        indices = collect(1:length(clusters_vector))
    end
    for i in proc_ids
        workers_suff_dict[i] = remotecall(
            create_suff_stats_dict_worker,
            i,
            group_pts,
            group_labels,
            group_sublabels,
            hyper_params,
            indices,
        )
    end
    suff_stats_vectors = [[] for i = 1:length(indices)]
    cluster_to_index = Dict([indices[i] => i for i = 1:length(indices)])
    workers_suff_dict_fetched = Dict([k => fetch(v) for (k, v) in workers_suff_dict])
    for (k, v) in workers_suff_dict_fetched
        for (cluster, suff) in v
            push!(suff_stats_vectors[cluster_to_index[cluster]], suff)
        end
    end
    for (k, v) in enumerate(suff_stats_vectors)
        if length(v) > 0
            cp_suff = reduce(aggregate_suff_stats, [x.cluster_suff for x in v])
            cpl_suff = reduce(aggregate_suff_stats, [x.l_suff for x in v])
            cpr_suff = reduce(aggregate_suff_stats, [x.r_suff for x in v])
            cluster_index = indices[k]
            leader_suff_dict[cluster_index] = thin_suff_stats(cp_suff, cpl_suff, cpr_suff)
        end
    end

    return leader_suff_dict
end

function update_suff_stats_posterior!(
    group::local_group,
    indices = nothing,
    use_leader::Bool = true,
)
    workers_suff_dict = Dict()
    if indices == nothing
        indices = collect(1:length(group.local_clusters))
    end
    if use_leader
        for i in collect(keys(leader_dict))
            workers_suff_dict[i] = remotecall(
                create_suff_stats_dict_node_leader,
                i,
                group.points,
                group.labels,
                group.labels_subcluster,
                group.model_hyperparams.distribution_hyper_params,
                leader_dict[i],
                indices,
            )
        end
    else
        for i in (nworkers() == 0 ? procs() : workers())
            workers_suff_dict[i] = @spawnat i create_suff_stats_dict_worker(
                group.points,
                group.labels,
                group.labels_subcluster,
                group.model_hyperparams.distribution_hyper_params,
                indices,
            )
        end
    end
    suff_stats_vectors = [[] for i = 1:length(indices)]
    cluster_to_index = Dict([indices[i] => i for i = 1:length(indices)])
    workers_suff_dict_fetched = Dict([k => fetch(v) for (k, v) in workers_suff_dict])
    for (k, v) in workers_suff_dict_fetched
        for (cluster, suff) in v
            push!(suff_stats_vectors[cluster_to_index[cluster]], suff)
        end
    end
    for (index, v) in enumerate(indices)
        # if outlier_mod > 0 && v == 1
        #     continue
        # end
        if length(suff_stats_vectors[index]) == 0
            continue
        end
        cluster = group.local_clusters[v]
        cp = cluster.cluster_params
        cp.cluster_params.suff_statistics = reduce(
            aggregate_suff_stats,
            [x.cluster_suff for x in suff_stats_vectors[index]],
        )
        cp.cluster_params_l.suff_statistics =
            reduce(aggregate_suff_stats, [x.l_suff for x in suff_stats_vectors[index]])
        cp.cluster_params_r.suff_statistics =
            reduce(aggregate_suff_stats, [x.r_suff for x in suff_stats_vectors[index]])
        cluster.points_count = cp.cluster_params.suff_statistics.N
        update_splittable_cluster_params!(cluster.cluster_params)
    end

end

function split_first_cluster_worker!(group::local_group)
    sub_labels = localpart(group.labels_subcluster)
    pts = localpart(group.points)
    sub_labels .= rand(1:2, length(sub_labels))
end

function split_cluster_local_worker!(
    labels,
    sub_labels,
    points,
    indices::Vector{Int64},
    new_indices::Vector{Int64},
)
    labels = localpart(labels)
    sub_labels = localpart(sub_labels)
    pts = localpart(points)
    for (i, index) in enumerate(indices)
        cluster_sub_labels = @view sub_labels[labels.==index]
        cluster_labels = @view labels[labels.==index]
        cluster_points = @view pts[:, labels.==index]
        cluster_labels[cluster_sub_labels.==2] .= new_indices[i]

        cluster_sub_labels .= rand(1:2, length(cluster_sub_labels))

    end
end

function split_cluster_local!(
    group::local_group,
    cluster::local_cluster,
    index::Int64,
    new_index::Int64,
)

    l_split = copy_local_cluster(cluster)
    l_split.cluster_params = create_splittable_from_params(
        cluster.cluster_params.cluster_params_r,
        group.model_hyperparams.α,
    )
    cluster.cluster_params = create_splittable_from_params(
        cluster.cluster_params.cluster_params_l,
        group.model_hyperparams.α,
    )

    l_split.points_count = l_split.cluster_params.cluster_params.suff_statistics.N
    cluster.points_count = cluster.cluster_params.cluster_params.suff_statistics.N

    group.local_clusters[new_index] = l_split

end

function merge_clusters_worker!(
    group::local_group,
    indices::Vector{Int64},
    new_indices::Vector{Int64},
)
    labels = localpart(group.labels)
    sub_labels = localpart(group.labels_subcluster)

    for (i, index) in enumerate(indices)
        cluster_sub_labels = @view localpart(group.labels_subcluster)[labels.==index]
        cluster_sub_labels .= 1
        cluster_sub_labels =
            @view localpart(group.labels_subcluster)[labels.==new_indices[i]]
        cluster_sub_labels .= 2
        labels[labels.==new_indices[i]] .= index
    end
end

function merge_clusters!(group::local_group, index_l::Int64, index_r::Int64)
    new_splittable_cluster = merge_clusters_to_splittable(
        group.local_clusters[index_l].cluster_params.cluster_params,
        group.local_clusters[index_r].cluster_params.cluster_params,
        group.model_hyperparams.α,
    )
    group.local_clusters[index_l].cluster_params = new_splittable_cluster
    group.local_clusters[index_l].points_count += group.local_clusters[index_r].points_count
    group.local_clusters[index_r].points_count = 0
    group.local_clusters[index_r].cluster_params.cluster_params.suff_statistics.N = 0
    group.local_clusters[index_r].cluster_params.splittable = false
end

function should_split_local!(
    should_split::AbstractArray{Float32,1},
    cluster_params::splittable_cluster_params,
    α::Float32,
    final::Bool,
)
    cpl = cluster_params.cluster_params_l
    cpr = cluster_params.cluster_params_r
    cp = cluster_params.cluster_params
    if final || cpl.suff_statistics.N == 0 || cpr.suff_statistics.N == 0
        should_split .= 0
        return
    end
    post = calc_posterior(cp.hyperparams, cp.suff_statistics)
    lpost = calc_posterior(cp.hyperparams, cpl.suff_statistics)
    rpost = calc_posterior(cp.hyperparams, cpr.suff_statistics)

    log_likihood_l = log_marginal_likelihood(cpl.hyperparams, lpost, cpl.suff_statistics)
    log_likihood_r = log_marginal_likelihood(cpr.hyperparams, rpost, cpr.suff_statistics)
    log_likihood = log_marginal_likelihood(cp.hyperparams, post, cp.suff_statistics)
    log_HR =
        log(α) +
        logabsgamma(cpl.suff_statistics.N)[1] +
        log_likihood_l +
        logabsgamma(cpr.suff_statistics.N)[1] +
        log_likihood_r - (logabsgamma(cp.suff_statistics.N)[1] + log_likihood)
    if log_HR > log(rand())
        should_split .= 1
    end
    # print("LOG HR IS: ")
    # println(log_HR)
end

function check_and_split!(group::local_group, final::Bool, init_type::AbstractString)
    split_arr = zeros(Float32, length(group.local_clusters))
    for (index, cluster) in enumerate(group.local_clusters)
        if outlier_mod > 0 && index == 1
            continue
        end
        if cluster.cluster_params.splittable == true &&
           cluster.cluster_params.cluster_params.suff_statistics.N > 1
            should_split_local!(
                (@view split_arr[index, :]),
                cluster.cluster_params,
                group.model_hyperparams.α,
                final,
            )
        end
    end
    new_index = length(group.local_clusters) + 1
    indices = Vector{Int64}()
    new_indices = Vector{Int64}()
    resize!(group.local_clusters, Int64(length(group.local_clusters) + sum(split_arr)))
    for i = 1:length(split_arr)
        if split_arr[i] == 1

            push!(indices, i)
            push!(new_indices, new_index)
            split_cluster_local!(group, group.local_clusters[i], i, new_index)
            new_index += 1
        end
    end
    all_indices = vcat(indices, new_indices)
    if length(indices) > 0
        for i in (nworkers() == 0 ? procs() : workers())
            @spawnat i split_cluster_local_worker!(
                group.labels,
                group.labels_subcluster,
                group.points,
                indices,
                new_indices,
            )
        end

        # INSERT HERE:
        if init_type == "splitnet_2d"
            for i in all_indices
                println("doing splitnet 2d init")
                splitnet_2d_init!(group,i)
            end
        elseif init_type == "splitnet_10d"
            for i in all_indices
                println("doing splitnet 10d init")
                splitnet_10d_init!(group,i)
            end
        elseif init_type == "splitnet_128d"
            for i in all_indices
                println("doing splitnet 128d init")
                splitnet_128d_init!(group,i)
            end
        elseif init_type == "eigensplit"
            for i in all_indices
                println("doing eigensplit init")
                eigensplit_init!(group,i)
            end
        elseif init_type == "kmeans"
            for i in all_indices
                println("doing kmeans init")
                kmeans_init!(group,i)
            end
        end

    end
    return all_indices
end

function check_and_merge!(group::local_group, final::Bool)
    mergable = zeros(Float32, 1)
    indices = Vector{Int64}()
    new_indices = Vector{Int64}()
    for i = 1:length(group.local_clusters)
        if outlier_mod > 0 && i == 1
            continue
        end
        for j = i+1:length(group.local_clusters)
            if (
                group.local_clusters[i].cluster_params.splittable == true &&
                group.local_clusters[j].cluster_params.splittable == true &&
                group.local_clusters[i].cluster_params.cluster_params.suff_statistics.N >
                0 &&
                group.local_clusters[j].cluster_params.cluster_params.suff_statistics.N > 0
            )
                should_merge!(
                    mergable,
                    group.local_clusters[i].cluster_params.cluster_params,
                    group.local_clusters[j].cluster_params.cluster_params,
                    group.model_hyperparams.α,
                    final,
                )
            end
            if mergable[1] == 1
                merge_clusters!(group, i, j)
                push!(indices, i)
                push!(new_indices, j)
            end
            mergable[1] = 0
        end
    end
    for i in (nworkers() == 0 ? procs() : workers())
        @spawnat i merge_clusters_worker!(group, indices, new_indices)
    end
    return indices
end

function sample_clusters!(group::local_group, first::Bool)
    points_count = Vector{Float32}()
    local_workers = procs(1)[2:end]
    cluster_params_futures = Dict()
    for (i, cluster) in enumerate(group.local_clusters)
        cluster_params_futures[i] =
            sample_cluster_params(cluster.cluster_params, group.model_hyperparams.α, first)
    end
    for (i, cluster) in enumerate(group.local_clusters)
        if outlier_mod > 0 && i == 1
            continue
        end
        cluster.cluster_params = fetch(cluster_params_futures[i])
        cluster.points_count = cluster.cluster_params.cluster_params.suff_statistics.N
        push!(points_count, cluster.points_count)
    end
    push!(points_count, group.model_hyperparams.α)
    group.weights = rand(Dirichlet(Float64.(points_count)))[1:end-1] .* (1 - outlier_mod)
    if outlier_mod > 0
        group.weights = vcat([outlier_mod], group.weights)
    end
end

function create_thin_cluster_params(cluster::local_cluster)
    return thin_cluster_params(
        cluster.cluster_params.cluster_params.distribution,
        cluster.cluster_params.cluster_params_l.distribution,
        cluster.cluster_params.cluster_params_r.distribution,
        cluster.cluster_params.lr_weights,
    )
end

function remove_empty_clusters_worker!(labels, pts_count)
    labels = localpart(labels)
    removed = 0
    for (index, count) in enumerate(pts_count)
        if count == 0
            labels[labels.>index-removed] .-= 1
            removed += 1
        end
    end
end

function remove_empty_clusters!(group::local_group)
    new_vec = Vector{local_cluster}()
    removed = 0
    pts_count = Vector{Int64}()
    for (index, cluster) in enumerate(group.local_clusters)
        push!(pts_count, cluster.points_count)
        if cluster.points_count > 0 ||
           (outlier_mod > 0 && index == 1) ||
           (outlier_mod > 0 && index == 2 && length(group.local_clusters) == 2)
            push!(new_vec, cluster)
        end
    end
    @sync for i in (nworkers() == 0 ? procs() : workers())
        @spawnat i remove_empty_clusters_worker!(group.labels, pts_count)
    end
    group.local_clusters = new_vec
end

function rand_subclusters_labels!(labels::AbstractArray{Int64,1})    #lr_arr = create_array(zeros(Float32,length(labels), 2))
    if size(labels, 1) == 0
        return
    end
    labels .= rand(1:2, size(labels, 1))
end

function reset_bad_clusters_worker!(
    indices::Vector{Int64},
    group_points::AbstractArray{Float32,2},
    group_labels::AbstractArray{Int64,1},
    group_labels_subcluster::AbstractArray{Int64,1},
)
    labels = localpart(group_labels)
    sub_labels = localpart(group_labels_subcluster)
    pts = localpart(group_points)
    for i in indices
        rand_subclusters_labels!((@view sub_labels[labels.==i]))
    end
end

function reset_splitted_clusters!(group::local_group, bad_clusters::Vector{Int64})

    for i in bad_clusters
        group.local_clusters[i].cluster_params.logsublikelihood_hist =
            ones(burnout_period + 5) * -Inf
    end
    for i in (nworkers() == 0 ? procs() : workers())
        @spawnat i reset_bad_clusters_worker!(
            bad_clusters,
            group.points,
            group.labels,
            group.labels_subcluster,
        )
    end
    update_suff_stats_posterior!(group, bad_clusters)
end

function reset_bad_clusters!(group::local_group)
    bad_clusters = Vector{Int64}()
    for (i, c) in enumerate(group.local_clusters)
        cl = c.cluster_params.cluster_params_l
        cr = c.cluster_params.cluster_params_r
        if cl.suff_statistics.N == 0 || cr.suff_statistics.N == 0
            push!(bad_clusters, i)
            c.cluster_params.logsublikelihood_hist = ones(burnout_period + 5) * -Inf
            c.cluster_params.splittable = false
        end
    end
    @sync for i in (nworkers() == 0 ? procs() : workers())
        @spawnat i reset_bad_clusters_worker!(
            bad_clusters,
            group.points,
            group.labels,
            group.labels_subcluster,
        )
    end
    update_suff_stats_posterior!(group, bad_clusters)
end

function broadcast_cluster_params(params_vector, weights_vector)
    @sync for (k, v) in leader_dict
        @spawnat k broadcast_to_node(params_vector, weights_vector, v)
    end
end

#Below code is to prevent some bug with deseraliztion
function broadcast_to_node(params_vector, weights_vector, proc_ids)
    refs = Dict()
    @sync for i in proc_ids
        refs[i] = @spawnat i set_global_data(params_vector, weights_vector)
    end
    alltrue = false
    responses = Dict([k => fetch(v) for (k, v) in refs])
    for (k, v) in responses
        while v == false
            v = remotecall_fetch(set_global_data, k, params_vector, weights_vector)
        end
    end
end

function set_global_data(params_vector, weights_vector)
    succ = true
    try
        global clusters_vector = params_vector
        global clusters_weights = weights_vector
    catch
        succ = false
    end
    return succ
end

function group_step(group::local_group, no_more_splits::Bool, final::Bool, first::Bool, init_type)
    sample_clusters!(group, false)
    broadcast_cluster_params(
        [create_thin_cluster_params(x) for x in group.local_clusters],
        group.weights,
    )
    sample_labels!(group, (hard_clustering ? true : final), no_more_splits)
    sample_sub_clusters!(group)
    update_suff_stats_posterior!(group)
    reset_bad_clusters!(group)
    if no_more_splits == false
        indices = []
        indices = check_and_split!(group, final, init_type)
        update_suff_stats_posterior!(group, indices)
        check_and_merge!(group, final)
    end
    remove_empty_clusters!(group)
    return
end

# PyTORCH MODEL:

function splitnet_2d_init!(group::local_group,cluster_num::Int64)
    pts = Array(group.points)
    labels = Array(group.labels)
    cluster_pts = pts[:,labels .==cluster_num]

    y = torch_model_2d.infer(cluster_pts)
    # preds = argmax(y, dims=2)
    preds = round.(y) .+ 1
    # split_labels =  dropdims([ i for i in preds ], dims=2)
    split_labels =  dropdims(preds, dims=2)

    sub_labels = localpart(group.labels_subcluster)
    sub_labels[labels .== cluster_num] .= split_labels

    # println("SPLITNET COUNT")
    # print(count(split_labels.==1))
    # print(count(split_labels.==2))
end

function splitnet_10d_init!(group::local_group,cluster_num::Int64)
    pts = Array(group.points)
    labels = Array(group.labels)
    cluster_pts = pts[:,labels .==cluster_num]

    y = torch_model_10d.infer(cluster_pts)
    # preds = argmax(y, dims=2)
    preds = round.(y) .+ 1
    # split_labels =  dropdims([ i for i in preds ], dims=2)
    split_labels =  dropdims(preds, dims=2)

    sub_labels = localpart(group.labels_subcluster)
    sub_labels[labels .== cluster_num] .= split_labels
end

function splitnet_128d_init!(group::local_group,cluster_num::Int64)
    pts = Array(group.points)
    labels = Array(group.labels)
    cluster_pts = pts[:,labels .==cluster_num]

    y = torch_model_128d.infer(cluster_pts)
    # preds = argmax(y, dims=2)
    preds = round.(y) .+ 1
    # split_labels =  dropdims([ i for i in preds ], dims=2)
    split_labels =  dropdims(preds, dims=2)

    sub_labels = localpart(group.labels_subcluster)
    sub_labels[labels .== cluster_num] .= split_labels
end

# KMEANS:

function kmeans_init!(group::local_group,cluster_num::Int64)
    pts = Array(group.points)
    labels = Array(group.labels)
    cluster_pts = pts[:,labels .==cluster_num]

    if size(cluster_pts,2) == 1
        return [1]
    end

    split_labels = kmeans(cluster_pts, 2).assignments
    sub_labels = localpart(group.labels_subcluster)
    sub_labels[labels .== cluster_num] .= split_labels

    # println("KMEANS COUNT")
    # print(count(split_labels.==1))
    # print(count(split_labels.==2))
end

# EIGENSPLIT

function eigensplit_init!(group::local_group,cluster_num::Int64)
    N = group.local_clusters[cluster_num].cluster_params.cluster_params.suff_statistics.N
    XXT = group.local_clusters[cluster_num].cluster_params.cluster_params.suff_statistics.S / N
    μ = group.local_clusters[cluster_num].cluster_params.cluster_params.suff_statistics.points_sum / N
    M = XXT - μ*μ'

    F = eigen(M)

    vals1 = copy(F.values)
    vals2 = copy(F.values)

    vals1[end] = 1e-5
    vals2[end-1] = 1e-5

    vals1, vals2

    subcov1 = F.vectors * Diagonal(vals1) * F.vectors'
    subcov2 = F.vectors * Diagonal(vals2) * F.vectors'
    invcov1 = inv(subcov1)
    invcov2 = inv(subcov2)

    @sync for i in workers()
        @spawnat i cov_split_worker!(group, cluster_num, invcov1,invcov2)
    end
end

function cov_split_worker!(group::local_group,cluster_num::Int64, invcov1, invcov2)
    malhan_dist1 = x -> x' * invcov1 * x
    malhan_dist2 = x -> x' * invcov2 * x
    one_or_two = x -> malhan_dist1(x) > malhan_dist2(x) ? 1 : 2
    labels = localpart(group.labels)
    sublabels = @view localpart(group.labels_subcluster)[labels .== cluster_num]
    pts = @view localpart(group.points)[:,labels .== cluster_num]
    for i=1:length(sublabels)
        sublabels[i] = one_or_two(pts[:,i])
    end
end

# function eigensplit_init!(group::local_group,cluster_num::Int64)
#     N = group.local_clusters[cluster_num].cluster_params.cluster_params.suff_statistics.N
#     XXT = group.local_clusters[cluster_num].cluster_params.cluster_params.suff_statistics.S / N
#     μ = group.local_clusters[cluster_num].cluster_params.cluster_params.suff_statistics.points_sum / N
#     M = XXT - μ*μ'
#     # println("vecs:")
#     F = eigen(M)
#     vecs = F.vectors
#     vals = F.values
#     # println("eigsolve:")
#     # v1 = eigsolve(M, 1, maxiter = 10)[2][1]
#     mxindx = findmax(vals)[2]
#     v1=vecs[mxindx,:]
#     means_futures = Vector{Future}(undef,nworkers())
#     for i in workers()
#         means_futures[i-1] =  @spawnat i tranform_points_worker!(group.points,group.labels,(@eval $cluster_num),(@eval $v1),(@eval $μ))
#     end
#     min_pts = [fetch(x) != nothing ? fetch(x)[1] : nothing for x in means_futures if fetch(x) != nothing]
#     max_pts = [fetch(x) != nothing ? fetch(x)[2] : nothing for x in means_futures if fetch(x) != nothing]
#     if length(min_pts) == 0
#         return
#     end

#     min_mean = minimum(min_pts)
#     max_mean = maximum(max_pts)
#     # min_mean = -100
#     # max_mean = +100
#     #Distributed K-Means
#     iter = 0
#     converged = false
#     while(iter < max_split_iter && !converged)
#         means_futures = Vector{Future}(undef,nworkers())
#         @sync for i in workers()
#             means_futures[i-1] = @spawnat i kmeans_iter_worker!(min_mean,max_mean)
#         end
#         min_means = [fetch(x) != nothing ? fetch(x)[1] : nothing for x in means_futures if fetch(x) != nothing]
#         max_means = [fetch(x) != nothing ? fetch(x)[2] : nothing for x in means_futures if fetch(x) != nothing]
#         min_sum = 0
#         min_count = 0
#         for x in min_means
#             min_sum += x[1]
#             min_count += x[2]
#         end
#         max_sum = 0
#         max_count = 0
#         for x in max_means
#             max_sum += x[1]
#             max_count += x[2]
#         end
#         new_min = min_sum / min_count
#         new_max = max_sum / max_count
#         if new_min == min_mean && new_max == max_mean
#             converged = true
#         else
#             min_mean = new_min
#             max_mean = new_max
#         end
#         iter += 1
#     end
#     @sync for i in workers()
#         @spawnat i set_smart_labels_in_worker!(cluster_num,group.labels, group.labels_subcluster,group.points)
#     end
# end

function set_smart_labels_in_worker!(cluster_num,labels,sub_labels,points)
    labels = localpart(labels)
    sub_labels = localpart(sub_labels)
    sub_labels[labels .== cluster_num] .= split_labels
end

function tranform_points_worker!(all_pts,all_labels,cluster_num::Int64,v,μ)

    pts = localpart(all_pts)
    labels = localpart(all_labels)
    rel_points = pts[:,labels .== cluster_num] .- μ
    global transformed_pts = v' * rel_points

    if length(transformed_pts) > 1
        return percentile((@view transformed_pts[:]),0.10), percentile((@view transformed_pts[:]),0.90)
    end
    return nothing
end

function kmeans_iter_worker!(min_mean,max_mean)
    global split_labels = [(abs(x - min_mean) < abs(x-max_mean) ? 1 : 2) for x in (@view transformed_pts[:])]
    if length(split_labels) > 0
        return (sum(transformed_pts[split_labels .== 1]), sum(split_labels .== 1)), (sum(transformed_pts[split_labels .== 2]), sum(split_labels .== 2))
    end
    return nothing
end

function eigens_init!(group::local_group,cluster_num::Int64)
    pts = Array(group.points)
    labels = Array(group.labels)
    cluster_pts = pts[:,labels .==cluster_num]

    if size(points,2) == 1
        return [1]
    end

    m = MultivariateStats.fit(PCA,points; maxoutdim=1)
    v1=projection(m)
    new_pts = transform(m,points)
    a,minpt = findmin(new_pts)
    a,maxpt = findmax(new_pts)
    minpt = minpt[2]
    maxpt = maxpt[2]
    result = kmeans(new_pts, 2, init = [minpt,maxpt])
    return result.assignments
end


text:
using GLMakie
using AbstractPlotting
AbstractPlotting.inline!(true)

function volume_video!(sim,data_func;name="file.mp4",duration=1,step=0.1)
    # Set up viz data and figure
    dat = Node(data_func(sim))
    fig = volume(dat,colorrange=(π,4π),algorithm=:absorption,
                    backgroundcolor = RGBf0(0.9, 0.9, 0.9))

    # Run simulation and update figure data
    t₀ = round(sim_time(sim))
    t = range(t₀,t₀+duration;step)
    record(fig,name,t) do tᵢ
        sim_step!(sim,tᵢ)
        dat[] = data_func(sim)
        println("simulation ",round(Int,(tᵢ-t₀)/duration*100),"% complete")
    end
    return sim,fig
end

function contour_video!(sim,data_func,geom_func;name="file.mp4",duration=1,step=0.1)
    # Set up viz data and figure
    fig = contour(geom_func(sim),levels=[0.5])
    dat = Node(data_func(sim))
    contour!(dat,levels=[-7,7],colormap=:balance,alpha=0.2,colorrange=[-7,7])

    # Run simulation and update figure data
    t₀ = round(sim_time(sim))
    t = range(t₀,t₀+duration;step)
    record(fig,name,t) do tᵢ
        sim_step!(sim,tᵢ)
        dat[] = data_func(sim)
        println("simulation ",round(Int,(tᵢ-t₀)/duration*100),"% complete")
    end
    return sim,fig
end


text:
# define the generative model for the multiple pedestrian crosswalk environment

### REWARD MODEL ##################################################################################

function POMDPs.reward(pomdp::OCPOMDP, s::OCState, a::OCAction, sp::OCState)
    r = 0.
    ego = sp[findfirst(EGO_ID, sp)]
    if is_crash(sp)
        r += pomdp.collision_cost
    end
    if ego.state.posG.x >= pomdp.ego_goal
        r += pomdp.goal_reward
    elseif a.acc > 0.
        r += pomdp.action_cost
    else
        r += pomdp.action_cost
    end
    return r
end

### TERMINAL STATES ###############################################################################

function POMDPs.isterminal(pomdp::OCPOMDP, s::OCState)
    ego = s[findfirst(EGO_ID, s)]
    return is_crash(s) || ego.state.posG.x >= pomdp.ego_goal
end

### TRANSITION MODEL ##############################################################################

function POMDPs.gen(::DDNNode{:sp}, pomdp::OCPOMDP, s::OCState, a::OCAction, rng::AbstractRNG)
    actions = Array{Any}(undef, length(s))
    pomdp.models[1].a = a
    sp = deepcopy(s) #XXX bad
    # max_id = 0
    # for veh in sp
    #     if veh.id > max_id
    #         max_id = veh.id
    #     end
    # end
    # if rand(rng) < pomdp.p_birth && max_id < pomdp.max_peds+1 && !pomdp.no_ped
    # if rand(rng) < 0.01 && max_id < pomdp.max_peds+1 && !pomdp.no_ped
    #     # println("Spawning new pedestrians")
    #     new_ped = initial_pedestrian(pomdp, sp, rng)
    #     pomdp.models[new_ped.id] = ConstantPedestrian(dt = pomdp.ΔT)#TODO parameterized
    #     push!(sp, new_ped)
    #     actions = Array{Any}(undef, length(sp))
    # end
    get_actions!(actions, sp, pomdp.env.roadway, pomdp.models)
    tick!(sp, pomdp.env.roadway, actions, pomdp.ΔT)
    clean_scene!(pomdp.env, sp)
    return sp
end

# function clean_scene!(env::CrosswalkEnv, scene::Scene)
#     for veh in scene
#         if veh.id == 1
#             continue
#         end
#         y = veh.state.posG.y
#         if y >= env.params.crosswalk_length/4
#             deleteat!(scene, findfirst(veh.id, scene))
#         end
#     end
# end

function AutomotiveDrivingModels.propagate(veh::Vehicle, action::OCAction, roadway::Roadway, Δt::Float64)
    x_ = veh.state.posG.x + veh.state.v*Δt + action.acc*Δt^2/2
    if x_ <= veh.state.posG.x
        x_ = veh.state.posG.x
    end
    v_ = veh.state.v + action.acc*Δt
    v_ = clamp(v_, 0., 8.) # parameterized
    return VehicleState(VecSE2(x_, veh.state.posG.y, veh.state.posG.θ), roadway, v_)
end

### INITIAL STATES ################################################################################

function POMDPs.initialstate(pomdp::OCPOMDP, rng::AbstractRNG)
    pomdp.no_ped = rand(rng) < pomdp.no_ped_prob
    scene = Scene()
    if !pomdp.no_ped
        for i=1:pomdp.max_peds
            if rand(rng) < pomdp.p_birth # pedestrian appear
                new_ped = initial_pedestrian(pomdp, scene, rng, true)
                pomdp.models[new_ped.id] = ConstantPedestrian(dt = pomdp.ΔT)
                push!(scene, new_ped)
            end
        end
    end
    ego = initial_ego(pomdp, rng)
    push!(scene, ego)
    return scene
end

"""
    initial_ego(pomdp::OCPOMDP, rng::AbstractRNG)
generate an initial state for the ego car, returns a Vehicle object
"""
function initial_ego(pomdp::OCPOMDP, rng::AbstractRNG)
    env = pomdp.env
    x0 = pomdp.ego_start
    y0 = 0. # parameterize
    v0 = rand(rng, Distributions.Uniform(6., pomdp.env.params.speed_limit)) #TODO parameterize
    return Vehicle(VehicleState(VecSE2(x0, y0, 0.0), env.roadway.segments[1].lanes[1], env.roadway, v0),
                   VehicleDef(), 1)
end

"""
    initial_pedestrian(pomdp::OCPOMDP, scene::Scene, rng::AbstractRNG)
Create a new pedestrian entity at its initial state
"""
function initial_pedestrian(pomdp::OCPOMDP, scene::Scene, rng::AbstractRNG, first_scene::Bool = false)
    env = pomdp.env
    # lateral position
    d_lat = Distributions.Uniform(env.params.roadway_length - env.params.crosswalk_width + 2, env.params.roadway_length + env.params.crosswalk_width-1)
    x0 = 0.5*rand(rng, d_lat)

    # vertical position
    y0 = pomdp.ped_start
    if first_scene
        y0 = rand(rng, pomdp.ped_start:pomdp.ped_goal)
    end

    # orientation
    θ = π/2

    #velocity
    v0 = rand(rng, Distributions.Uniform(0., env.params.ped_max_speed))
    cw_roadway = Roadway([RoadSegment(2, [env.crosswalk])]);
    ped_initialstate = VehicleState(VecSE2(x0, y0, θ), env.crosswalk, cw_roadway, v0);

    # new id, increment last id
    max_id = 0
    for veh in scene
        if veh.id > max_id
            max_id = veh.id
        end
    end
    id = max_id + 1
    if max_id == 0
        id = 2
    end

    return Vehicle(ped_initialstate, PEDESTRIAN_DEF, id)
end

function POMDPs.initialstate_distribution(pomdp::OCPOMDP)
    return GenerativeDist(pomdp)
end

### Observations ##################################################################################

# uncomment for LIDAR observation
# function POMDPs.generate_o(pomdp::OCPOMDP, s::OCState, a::OCAction, sp::OCState, rng::AbstractRNG)
#     observe!(pomdp.sensor, sp, pomdp.env.roadway, EGO_ID)
#     return pomdp.sensor
# end
#
# function POMDPs.generate_o(pomdp::OCPOMDP, s::OCState, rng::AbstractRNG)
#     observe!(pomdp.sensor, s, pomdp.env.roadway, EGO_ID)
#     return pomdp.sensor
# end
#
# function POMDPs.convert_o(::Type{Vector{Float64}}, o::OCObs, pomdp::OCPOMDP)
#     return [o.ranges/pomdp.sensor.max_range, o.range_rates/pomdp.env.params.speed_limit]
# end

# function POMDPs.convert_o(::Type{Vector{Float64}}, o::OCObs, pomdp::OCPOMDP)
#     return [o.ranges/pomdp.sensor.max_range, o.range_rates/pomdp.env.params.speed_limit]
# end

# uncomment for vector representation
function POMDPs.gen(::DDNNode{:o}, pomdp::OCPOMDP, s::OCState, a::OCAction, sp::OCState, rng::AbstractRNG)
    pos_noise = pomdp.pos_obs_noise
    vel_noise = pomdp.vel_obs_noise
    o = zeros(2 + 2*pomdp.max_peds)
    ego = sp[findfirst(EGO_ID, sp)].state
    o[1] = ego.posG.x
    o[2] = ego.v
    ped_off = get_off_the_grid(pomdp)
    for i=2:pomdp.max_peds+1
        o[2*i - 1] = ped_off.posG.y
        o[2*i] = 0.
    end
    for veh in sp
        if veh.id == EGO_ID
            continue
        end
        @assert veh.id <= pomdp.max_peds+1
        ped = veh.state
        if is_observable_fixed(ped, ego, pomdp.env)
            o[2*veh.id - 1] = ped.posG.y + pos_noise*randn(rng)
            o[2*veh.id] = ped.v + vel_noise*randn(rng)
        end
    end
    # rescale
    o[1] /= pomdp.ego_goal
    o[2] /= pomdp.env.params.speed_limit
    for i=2:pomdp.max_peds+1
        o[2*i - 1] /= pomdp.ped_goal
        o[2*i] /= 2. # XXX parameterized
    end
    return o
end

function POMDPs.initialobs(pomdp::OCPOMDP, s::OCState, rng::AbstractRNG)
    return gen(DDNNode(:o), pomdp, s, OCAction(0.), s, rng::AbstractRNG)
end

# function POMDPs.convert_o(::Type{Vector{Float64}}, o::OCObs, pomdp::OCPOMDP)
#     o[1] /= pomdp.ego_goal
#     o[2] /= pomdp.env.params.speed_limit
#     for i=2:pomdp.max_peds+1
#         o[2*i - 1] /= pomdp.ped_goal
#         o[2*i] /= 2. # XXX parameterized
#     end
#     return o
# end

### All together ##################################################################################

function POMDPs.gen(::DDNOut{(:sp,:o,:r)}, pomdp::OCPOMDP, s::OCState, a::OCAction, rng::AbstractRNG)
    sp = gen(DDNOut(:sp), pomdp, s, a, rng)
    o = gen(DDNNode(:o), pomdp, s, a, sp, rng)
    r = reward(pomdp, s, a, sp)
    return sp, o, r
end

#### Helpers

"""
    get_off_the_grid(pomdp::OCPOMDP)
return the off the grid pedestrian state
"""
function get_off_the_grid(pomdp::OCPOMDP)
    env = pomdp.env
    x_ped = 0.5*(env.params.roadway_length - env.params.crosswalk_width + 1)
    cl = pomdp.env.params.crosswalk_length
    pos = VecSE2(x_ped, -cl/2 - 1, pi/2)
    return VehicleState(pos, pomdp.env.roadway, Inf)
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
function part1(input)
    bits = parse_input(input)
    version_sum = Ref(0)
    decode_packet(bits, version_sum)
    return version_sum[]
end

function part2(input)
    bits = parse_input(input)
    return decode_packet(bits, Ref(0))
end

mutable struct Bits
    data::Vector{UInt8}
    pos::Int
    bits::UInt64
    bits_left::Int
end

Base.length(bits::Bits) = 4 * (length(bits.data) - bits.pos + 1) + bits.bits_left

function get_bits!(bits, n)
    if bits.bits_left < n
        while bits.bits_left <= 60 && bits.pos <= length(bits.data)
            bits.bits <<= 4
            c = bits.data[bits.pos]
            bits.bits |= c > UInt8('9') ? (c - 0x37) : (c - 0x30)
            bits.pos += 1
            bits.bits_left += 4
        end
    end
    bits.bits_left -= n
    out = bits.bits >> bits.bits_left
    bits.bits &= 1 << bits.bits_left - 1
    return out
end

function parse_input(input)
    return Bits(input.data, 1, UInt64(0), 0)
end

function decode_packet(bits, version_sum)
    version, id = get_bits!(bits, 3), get_bits!(bits, 3)
    version_sum[] += version
    if id == 4
        value = 0
        while true
            more = get_bits!(bits, 1) == 1
            value = (value << 4) | get_bits!(bits, 4)
            more || return value
        end
    else
        values = Int[]
        if get_bits!(bits, 1) == 0
            num_bits = get_bits!(bits, 15)
            target = length(bits) - num_bits
            while length(bits) > target
                push!(values, decode_packet(bits, version_sum))
            end
        else
            num_packages = get_bits!(bits, 11)
            for i = 1:num_packages
                push!(values, decode_packet(bits, version_sum))
            end
        end
        id < 4 && return (sum, prod, minimum, maximum)[id + 1](values)
        return Int((>, <, ==)[id - 4](values...))
    end
end


text:
using GraphMatFun, LsqFit, LinearAlgebra,Random
include("exp_reset_all.jl");
#include("simulationtools.jl");
include("simulationtools.jl");

f=exp
m=7
target_n=200 # Measure error wrt this discr
n=50;        # Optimize wrt this discr
rho=3.59;  # SID rho m=7

its=5;
droptol0=1e-10;
opt_kwargs=Dict(:logger=> 0,:γ0 => 0.5,:droptol => droptol0,
                :linlsqr => :real_svd,:maxit => its)
base = init_state_mult(f,rho,n,m; eltype=Complex{BigFloat})
base.params[:n]=n;
base.params[:target_n]=target_n;
base.params[:opt_kwargs]=opt_kwargs;
base.params[:kickit_mode]=1;

# First the optimization free
sid_org=init_state_mult!(deepcopy(base),:sid,m,showmeta=true)

#bbc_org=init_state_mult!(deepcopy(base),:bbc,m,showmeta=true)

datadir=joinpath("..","..","data","exp");
#filename=joinpath(datadir,"exp_sastre_m6_opt_rho2_7.cgr");
filename=joinpath(datadir,"exp_sastre_m6.cgr");
sastre_org=init_state_file!(deepcopy(base),:sastre,filename,showmeta=true,
                            scale_and_square=true)
#filename="simulations/newgraphs/exp_m4_sastre+GN_0_69.cgr"
#sastre_org=init_state_mult!(deepcopy(base),:sastre,m,showmeta=true)

#sastre_org.cref=get_degopt_crefs(sastre_org.graph);

ps_org=init_state_mult!(deepcopy(base),:ps,m,showmeta=true)

mono_org=init_state_mult!(deepcopy(base),:mono,m,showmeta=true)

(mono,simlist,commandlist)=
     run_sequence(mono_org,"sssGGGssdddsssddssdddskKssddssdddddssssddddssdddsssddddssssdddssssssddsddsssddssddsdsgggggggdssssssGGGGGGGssszq");

(ps,simlist,commandlist)=
     run_sequence(ps_org, "ksssGGGssdddddddddsddddddsdddddddddddsssddddddddddssgggggsssGGGssssGGsddsssddssddssskssddssddsssdddsssdddsggggggsssssssssssggsggssssssskssssgGGGGGGssssssssssssssGGGGGssszq");

filename=joinpath(datadir,"exp_sastre_m6_opt_rho2_7.cgr");
sastre0=init_state_file!(deepcopy(base),:sastre,filename,showmeta=true,
                            scale_and_square=true)
(sastre,simlist,commandlist)=
     run_sequence(sastre0,
                  "sdddddddddddddddddddddddddddddddddsdssGGGssddssSzq");

sid0=init_state_mult!(deepcopy(base),:sid,m,showmeta=true)
(sid,simlist,commandlist)=
     run_sequence(sid0,
     "ddddddSSSzkSSSSGGGSSSdddddSddSddSSddSddSSdddSdddSSSdddSddSzdddSdddSSSzSdddSSSNSddSdsSS2Sz2zq");
# "ddddddSSSzkSSSSGGGSSSdddddSddSddSSddSddSSdddSdddSSSdddSddSzdddSdddSSSzSdddSSSNSddSds2SSSq");

#(bbc,simlist,commandlist)=
#     run_sequence(bbc_org,
#     "sssssddssssddsssssssssssdddssdddssssdddssdddssssssssssdggdq")
#

include("exp_print_all.jl");
println("Run include(exp_save_all.jl) if you want to save");

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
# uses file produced by `11.6_produce.jl`.
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("colorscheme.jl"))
using Random, OrdinaryDiffEq
import GLMakie

data = wload(datadir("FitzHugh", "spiralwave.jld2"))
@unpack uout, tvec, params = data

heatobs = GLMakie.Observable(uout[1][1,:,:])
tobs = GLMakie.Observable(0.0)
titobs = GLMakie.lift(t -> "t = $(t)", tobs)

fig = GLMakie.Figure(resolution = (600, 550))
ax = fig[1,1] = GLMakie.Axis(fig; title = titobs)
hmap = GLMakie.heatmap!(ax, heatobs; colormap = :tokyo, colorrange = (-0.2, 1))
cb = GLMakie.Colorbar(fig[1, 2], hmap; width = 20)
display(fig)

GLMakie.record(
        fig, projectdir("animations", "11", "fitzhugh_spiralwave.mp4"),
        1:length(tvec); framerate = 15
    ) do i

    tobs[] = tvec[i]
    heatobs[] = uout[i][1,:,:]
end

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using Revise
using Antenna

elem_pattern = include("read_element_re.jl")
ploting_pattern(elem_pattern,min=-35)
set_param(f=2.7e9)
point = point_rectangle(Nx = 8,Ny=8, dx=0.111/2, dy=0.111/2,pattern=elem_pattern)
# point = point_rectangle(Nx = 4,Ny=4, dx=0.111/2, dy=0.111/2,pattern=elem_pattern)
res_pattern = cal_pattern(point, point, 0,0)
@show radiated_power(res_pattern)
ploting_pattern(res_pattern, min=-30)
Antenna.ploting_pattern_2D(res_pattern)


text:
# # Batched Generalized Minimal Residual
# In this tutorial we describe the basics of using the batched gmres iterative solver.
# At the end you should be able to
# 1. Use BatchedGeneralizedMinimalResidual to solve batches of linear systems
# 2. Construct a columnwise linear solver with BatchedGeneralizedMinimalResidual

# ## What is the Generalized Minimal Residual Method?
# The  Generalized Minimal Residual Method (GMRES) is a [Krylov subspace](https://en.wikipedia.org/wiki/Krylov_subspace) method for solving linear systems:
# ```math
#  Ax = b
# ```
# See the [wikipedia](https://en.wikipedia.org/wiki/Generalized_minimal_residual_method) for more details.

# ## What is the Batched Generalized Minimal Residual Method?
# As the name suggests it solves a whole bunch of independent GMRES problems

# ## Basic Example
# First we must load a few things
using ClimateMachine,
    ClimateMachine.LinearSolvers,
    ClimateMachine.BatchedGeneralizedMinimalResidualSolver
using LinearAlgebra, Random, Plots

# Next we define two linear systems that we would like to solve simultaneously.
# The matrix for the first linear system is
A1 = [
    2.0 -1.0 0.0
    -1.0 2.0 -1.0
    0.0 -1.0 2.0
];
# And the right hand side is
b1 = ones(typeof(1.0), 3);
# The exact solution to the first linear system is
x1_exact = [1.5, 2.0, 1.5];
# The matrix for the first linear system is
A2 = [
    2.0 -1.0 0.0
    0.0 2.0 -1.0
    0.0 0.0 2.0
];
# And the right hand side is
b2 = ones(typeof(1.0), 3);
# The exact solution to second linear system is
x2_exact = [0.875, 0.75, 0.5];

# We now define a function that performs the action of each linear operator independently.
function closure_linear_operator(A1, A2)
    function linear_operator!(x, y)
        mul!(view(x, :, 1), A1, view(y, :, 1))
        mul!(view(x, :, 2), A2, view(y, :, 2))
        return nothing
    end
    return linear_operator!
end;

# To understand how this works let us construct an instance
# of the linear operator and apply it to a vector
linear_operator! = closure_linear_operator(A1, A2);
# Let us see what the action of this linear operator is
y1 = ones(typeof(1.0), 3);
y2 = ones(typeof(1.0), 3) * 2.0;
y = [y1 y2];
x = copy(y);
linear_operator!(x, y);
x
# We see that the first column is `A1 * [1 1 1]'`
# and the second column is `A2 * [2 2 2]'`
# that is,
[A1 * y1 A2 * y2]

# We are now ready to set up our Batched Generalized Minimal Residual solver
# We must now set up the right hand side of the linear system
b = [b1 b2];
# as well as the exact solution, (to verify convergence)
x_exact = [x1_exact x2_exact];
# !!! warning
#     For BatchedGeneralizedMinimalResidual the assumption is that each column of b is independent and corresponds to a batch. This will come back later.

# We now use an instance of the solver
linearsolver = BatchedGeneralizedMinimalResidual(b);
# As well as an initial guess, denoted by the variable x
x1 = ones(typeof(1.0), 3);
x2 = ones(typeof(1.0), 3);
x = [x1 x2];
# To solve the linear system, we just need to pass to the linearsolve! function
iters = linearsolve!(linear_operator!, linearsolver, x, b)
# which is guaranteed to converge in 3 iterations since `length(b1)=length(b2)=3`
# We can now check that the solution that we computed, x
x
# has converged to the exact solution
x_exact
# Which indeed it has.
# ## Advanced Example

# We now go through a more advanced application of the Batched Generalized Minimal Residual solver
# !!! warning
#     Iterative methods should be used with preconditioners!
# The first thing we do is define a linear operator that mimics
# the behavior of a columnwise operator in ClimateMachine
function closure_linear_operator!(A, tup)
    function linear_operator!(y, x)
        alias_x = reshape(x, tup)
        alias_y = reshape(y, tup)
        for i6 in 1:tup[6]
            for i4 in 1:tup[4]
                for i2 in 1:tup[2]
                    for i1 in 1:tup[1]
                        tmp = alias_x[i1, i2, :, i4, :, i6][:]
                        tmp2 = A[i1, i2, i4, i6] * tmp
                        alias_y[i1, i2, :, i4, :, i6] .=
                            reshape(tmp2, (tup[3], tup[5]))
                    end
                end
            end
        end
    end
end;
# Next we define the array structure of an MPIStateArray
# in its true high dimensional form
tup = (2, 2, 5, 2, 10, 2);
# We define our linear operator as a random matrix
Random.seed!(1234);
B = [
    randn(tup[3] * tup[5], tup[3] * tup[5])
    for i1 in 1:tup[1], i2 in 1:tup[2], i4 in 1:tup[4], i6 in 1:tup[6]
];
columnwise_A = [
    B[i1, i2, i4, i6] + 3 * (i1 + i2 + i4 + i6) * I
    for i1 in 1:tup[1], i2 in 1:tup[2], i4 in 1:tup[4], i6 in 1:tup[6]
];
# as well as its inverse
columnwise_inv_A = [
    inv(columnwise_A[i1, i2, i4, i6])
    for i1 in 1:tup[1], i2 in 1:tup[2], i4 in 1:tup[4], i6 in 1:tup[6]
];
columnwise_linear_operator! = closure_linear_operator!(columnwise_A, tup);
columnwise_inverse_linear_operator! =
    closure_linear_operator!(columnwise_inv_A, tup);
# The structure of an MPIStateArray is related to its true
# higher dimensional form as follows:
mpi_tup = (tup[1] * tup[2] * tup[3], tup[4], tup[5] * tup[6]);
# We now define the right hand side of our Linear system
b = randn(mpi_tup);
# As well as the initial guess
x = copy(b);
x += randn(mpi_tup) * 0.1;
# In the previous tutorial we mentioned that it is assumed that
# the right hand side is an array whose column vectors all independent linear
# systems. But right now the array structure of ``x`` and ``b`` do not follow
# this requirement.
# To handle this case we must pass in additional arguments that tell the
# linear solver how to reconcile these differences.
# The first thing that the linear solver must know of is the higher tensor
# form of the MPIStateArray, which is just the `tup` from before
reshape_tuple_f = tup;
# The second thing it needs to know is which indices correspond to a column
# and we want to make sure that these are the first set of indices that appear
# in the permutation tuple (which can be thought of as enacting
# a Tensor Transpose).
permute_tuple_f = (5, 3, 4, 6, 1, 2);
# It has this format since the 3 and 5 index slots
# are the ones associated with traversing a column. And the 4 index
# slot corresponds to a state.
# We also need to tell our solver which kind of Array struct to use
ArrayType = Array;

# We are now ready to finally define our linear solver, which uses a number
# of keyword arguments
gmres = BatchedGeneralizedMinimalResidual(
    b,
    ArrayType = ArrayType,
    m = tup[3] * tup[5] * tup[4],
    n = tup[1] * tup[2] * tup[6],
    reshape_tuple_f = reshape_tuple_f,
    permute_tuple_f = permute_tuple_f,
    atol = eps(Float64) * 10^2,
    rtol = eps(Float64) * 10^2,
);
# ```m``` is the number of gridpoints along a column. As mentioned previously,
# this is `tup[3]*tup[5]*tup[4]`. The ```n``` term corresponds to the batch size
# or the number of columns in this case. ```atol``` and ```rtol``` are relative and
# absolute tolerances

# All the hard work is done, now we just call our linear solver
iters = linearsolve!(
    columnwise_linear_operator!,
    gmres,
    x,
    b,
    max_iters = tup[3] * tup[5] * tup[4],
)
# We see that it converged in less than `tup[3]*tup[5] = 50` iterations.
# Let us verify that it is indeed correct by computing the exact answer
# numerically and comparing it against the iterative solver.
x_exact = copy(x);
columnwise_inverse_linear_operator!(x_exact, b);
# Now we can compare with some norms
norm(x - x_exact) / norm(x_exact)
columnwise_linear_operator!(x_exact, x);
norm(x_exact - b) / norm(b)
# Which we see are small, given our choice of ```atol``` and ```rtol```.
# The struct also keeps a record of its convergence rate
# in the residual member. The convergence rate of each column
# can be visualized via
plot(log.(gmres.residual[1:iters, :]) / log(10));
plot!(legend = false, xlims = (1, iters), ylims = (-15, 2));
plot!(ylabel = "log10 residual", xlabel = "iterations")


text:
using Documenter, SpikingNN

makedocs(;
    modules=[SpikingNN],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Usage" => [
            "Neurons" => "neuron.md",
            "Populations" => "population.md",
            "Inputs" => "input.md",
            "Networks" => "network.md"
        ],
        "Models" => [
            "Neuron Models" => "neuron-models.md",
            "Synapse Models" => "synapse-models.md",
            "Threshold Models" => "threshold-models.md"
        ],
    ],
    repo="https://github.com/darsnack/SpikingNN.jl/blob/{commit}{path}#L{line}",
    sitename="SpikingNN.jl",
    authors="Kyle Daruwalla",
    assets=String[],
)

deploydocs(;
    repo="github.com/darsnack/SpikingNN.jl",
)


text:
function instantiate_eltype(ex::Expr)
    if istensor(ex)
        return Expr(:call, :eltype, gettensorobject(ex))
    elseif ex.head == :call && (ex.args[1] == :+ || ex.args[1] == :- || ex.args[1] == :* || ex.args[1] == :/)
        if length(ex.args) > 2
            return Expr(:call, :promote_type, map(instantiate_eltype, ex.args[2:end])...)
        else
            return instantiate_eltype(ex.args[2])
        end
    elseif ex.head == :call && ex.args[1] == :conj
        return instantiate_eltype(ex.args[2])
    elseif isscalarexpr(ex)
        return :(typeof($ex))
    else
        # return :(eltype($ex)) # would probably lead to doing the same operation twice
        throw(ArgumentError("unable to determine eltype"))
    end
end
instantiate_eltype(ex) = Expr(:call, :typeof, ex)

function instantiate_scalar(ex::Expr)
    if ex.head == :call && ex.args[1] == :scalar
        @assert length(ex.args) == 2 && istensorexpr(ex.args[2])
        return :(scalar($(instantiate(nothing, 0, ex.args[2], 1, [], [], true))))
    elseif ex.head == :call
        return Expr(ex.head, ex.args[1], map(instantiate_scalar, ex.args[2:end])...)
    else
        return Expr(ex.head, map(instantiate_scalar, ex.args)...)
    end
end
instantiate_scalar(ex::Symbol) = ex
instantiate_scalar(ex) = ex

function instantiate(dst, β, ex::Expr, α, leftind::Vector{Any}, rightind::Vector{Any}, istemporary = false)
    if isgeneraltensor(ex)
        return instantiate_generaltensor(dst, β, ex, α, leftind, rightind, istemporary)
    elseif ex.head == :call && (ex.args[1] == :+ || ex.args[1] == :-) # linear combination
        return instantiate_linearcombination(dst, β, ex, α, leftind, rightind, istemporary)
    elseif ex.head == :call && ex.args[1] == :* && length(ex.args) == 3 # multiplication: should be pairwise by now
        if isscalarexpr(ex.args[2])
            return instantiate(dst, β, ex.args[3], Expr(:call, :*, instantiate_scalar(ex.args[2]), α), leftind, rightind, istemporary)
        elseif isscalarexpr(ex.args[3])
            return instantiate(dst, β, ex.args[2], Expr(:call, :*, α, instantiate_scalar(ex.args[3])), leftind, rightind, istemporary)
        else
            return instantiate_contraction(dst, β, ex, α, leftind, rightind, istemporary)
        end
    elseif ex.head == :call && ex.args[1] == :/ && length(ex.args) == 3
        return instantiate(dst, β, ex.args[2], Expr(:call, :/, α, instantiate_scalar(ex.args[3])), leftind, rightind, istemporary)
    elseif ex.head == :call && ex.args[1] == :\ && length(ex.args) == 3
        return instantiate(dst, β, ex.args[3], Expr(:call, :\, instantiate_scalar(ex.args[2]), α), leftind, rightind, istemporary)
    end
    throw(ArgumentError("problem with parsing $ex"))
end

function instantiate_generaltensor(dst, β, ex::Expr, α, leftind::Vector{Any}, rightind::Vector{Any}, istemporary = false)
    src, srcleftind, srcrightind, α2, conj = decomposegeneraltensor(ex)
    srcind = vcat(srcleftind, srcrightind)
    conjarg = conj ? :(:C) : :(:N)

    p1 = (map(l->findfirst(isequal(l), srcind), leftind)...,)
    p2 = (map(l->findfirst(isequal(l), srcind), rightind)...,)

    αsym = gensym()
    if dst === nothing
        dst = gensym()
        if istemporary
            initex = quote
                $αsym = $α*$α2
                $dst = cached_similar_from_indices($(QuoteNode(dst)), promote_type(eltype($src), typeof($αsym)), $p1, $p2, $src, $conjarg)
            end
        else
            initex = quote
                $αsym = $α*$α2
                $dst = similar_from_indices(promote_type(eltype($src),typeof($αsym)), $p1, $p2, $src, $conjarg)
            end
        end
    else
        initex = :($αsym = $α*$α2)
    end

    if hastraceindices(ex)
        traceind = unique(setdiff(setdiff(srcind,leftind), rightind))
        q1 = (map(l->findfirst(isequal(l), srcind), traceind)...,)
        q2 = (map(l->findlast(isequal(l), srcind), traceind)...,)
        if any(x->(x===nothing), (p1...,p2...,q1...,q2...)) ||
            !isperm((p1...,p2...,q1...,q2...)) ||
            length(srcind) != length(leftind) + length(rightind) + 2*length(traceind)
            err = "trace: $(tuple(srcleftind..., srcrightind...)) to $(tuple(leftind..., rightind...)))"
            return :(throw(IndexError($err)))
        end
        return quote
            $initex
            trace!($α * $α2, $src, $conjarg, $β, $dst, $p1, $p2, $q1, $q2)
            # $dst
        end
    else
        if any(x->(x===nothing), (p1..., p2...)) || !isperm((p1...,p2...)) ||
            length(srcind) != length(leftind) + length(rightind)
            err = "add: $(tuple(srcleftind..., srcrightind...)) to $(tuple(leftind..., rightind...)))"
            return :(throw(IndexError($err)))
        end
        return quote
            $initex
            add!($α * $α2, $src, $conjarg, $β, $dst, $p1, $p2)
            # $dst
        end
    end
end
function instantiate_linearcombination(dst, β, ex::Expr, α, leftind::Vector{Any}, rightind::Vector{Any}, istemporary = false)
    if ex.head == :call && (ex.args[1] == :+ || ex.args[1] == :-) # addition: add one by one
        if dst === nothing
            αnew = Expr(:call, :*, α, Expr(:call, :one, instantiate_eltype(ex)))
            ex1 = instantiate(dst, β, ex.args[2], αnew, leftind, rightind, istemporary)
            dst = gensym()
            returnex = :($dst = $ex1)
        else
            returnex = instantiate(dst, β, ex.args[2], α, leftind, rightind, istemporary)
        end
        αnew = (ex.args[1] == :-) ? Expr(:call, :-, α) : α
        for k = 3:length(ex.args)
            ex1 = instantiate(dst, true, ex.args[k], αnew, leftind, rightind)
            returnex = quote
                $returnex
                $ex1
            end
        end
        return quote
            $returnex
            # $dst
        end
    else
        throw(ArgumentError("unable to instantiate linear combination: $ex"))
    end
end
function instantiate_contraction(dst, β, ex::Expr, α, leftind::Vector{Any}, rightind::Vector{Any}, istemporary = false)
    @assert ex.head == :call && ex.args[1] == :* && length(ex.args) == 3 &&
        istensorexpr(ex.args[2]) && istensorexpr(ex.args[3])
    exA = ex.args[2]
    exB = ex.args[3]

    indA = getindices(exA)
    indB = getindices(exB)
    cind = intersect(indA, indB)
    indC = vcat(leftind, rightind)
    oindA = intersect(indA, indC) # in the order they appear in A
    oindB = intersect(indB, indC) # in the order they appear in B

    symA = gensym()
    symB = gensym()
    symC = gensym()
    symTC = gensym()

    # prepare tensors or tensor expressions
    if dst === nothing
        TA = instantiate_eltype(exA)
        TB = instantiate_eltype(exB)
        TC = Expr(:call, :promote_type, TA, TB, :(typeof($α)))
    else
        TC = Expr(:call, :eltype, dst)
    end

    if !isgeneraltensor(exA) || hastraceindices(exA)
        initA = instantiate(nothing, false, exA, true, oindA, cind, true)
        poA = ((1:length(oindA))...,)
        pcA = length(oindA) .+ ((1:length(cind))...,)
        conjA = :(:N)
        initA = Expr(:(=), symA, initA)
        αA = 1
    else
        A, indlA, indrA, αA, conj = decomposegeneraltensor(exA)
        indA = vcat(indlA, indrA)
        poA = (map(l->findfirst(isequal(l), indA), oindA)...,)
        pcA = (map(l->findfirst(isequal(l), indA), cind)...,)
        TA = dst === nothing ? :(float(eltype($A))) : :(eltype($dst))
        conjA = conj ? :(:C) : :(:N)
        initA = Expr(:(=), symA, A)
    end

    if !isgeneraltensor(exB) || hastraceindices(exB)
        initB = instantiate(nothing, false, exB, true, cind, oindB, true)
        poB = length(cind) .+ ((1:length(oindB))...,)
        pcB = ((1:length(cind))...,)
        conjB = :(:N)
        initB = Expr(:(=), symB, initB)
        αB = 1
    else
        B, indlB, indrB, αB, conj = decomposegeneraltensor(exB)
        indB = vcat(indlB, indrB)
        poB = (map(l->findfirst(isequal(l), indB), oindB)...,)
        pcB = (map(l->findfirst(isequal(l), indB), cind)...,)
        conjB = conj ? :(:C) : :(:N)
        initB = Expr(:(=), symB, B)
    end

    oindAB = vcat(oindA, oindB)
    p1 = (map(l->findfirst(isequal(l), oindAB), leftind)...,)
    p2 = (map(l->findfirst(isequal(l), oindAB), rightind)...,)
    if any(x->(x===nothing), (poA..., pcA..., poB..., pcB..., p1..., p2...)) ||
        !(isperm((poA...,pcA...)) && length(indA) == length(poA)+length(pcA)) ||
        !(isperm((pcB...,poB...)) && length(indB) == length(poB)+length(pcB)) ||
        !(isperm((p1...,p2...)) && length(oindAB) == length(p1)+length(p2))
        err = "contraction: $(tuple(leftind..., rightind...)) from $(tuple(indA...,)) and $(tuple(indB...,)))"
        return :(throw(IndexError($err)))
    end
    if dst === nothing
        if istemporary
            initC = :($symC = cached_similar_from_indices($(QuoteNode(symC)), $symTC, $poA, $poB, $p1, $p2, $symA, $symB, $conjA, $conjB))
        else
            initC = :($symC = similar_from_indices($symTC, $poA, $poB, $p1, $p2, $symA, $symB, $conjA, $conjB))
        end
    else
        initC = :($symC = $dst)
    end

    return quote
        $symTC = $TC
        $initA
        $initB
        $initC
        contract!($α*$αA*$αB, $symA, $conjA, $symB, $conjB, $β, $symC,
                    $poA, $pcA, $poB, $pcB, $p1, $p2,
                    $((gensym(),gensym(),gensym())))
        # $symC
    end
end


text:
using Random
using LinearAlgebra: tril, triu, diagind
import ..LMAX, ..Norms

@testset "Equality of legendre[!] ($(typeof(N)) with argument type $T)" for N in Norms, T in NumTypes
    x = T(0.5)
    # On the first pass, use out-of-place interface
    #
    # Fill the full matrix of values just once.
    λ2 = legendre(N, 0:LMAX, 0:LMAX, x)
    for m in 0:LMAX
        # Check that vector filling function returns same answers
        λ1 = legendre(N, 0:LMAX, m, x)
        @test @view(λ1[m+1:end]) == @view(λ2[m+1:end,m+1])
        for ℓ in m:LMAX
            # Then check that scalar output matches as well.
            @test legendre(N, ℓ, m, x) == λ1[ℓ+1]
        end
    end

    # Now on the second pass, use the in-place interface
    λ2 = fill(T(NaN), LMAX+1, LMAX+1)
    λ1 = fill(T(NaN), LMAX+1)
    λ0 = fill(T(NaN))

    legendre!(N, λ2, LMAX, LMAX, x)
    # Test for equality in legendre! and legendre; wrap with tril to exlude NaNs used later
    @test tril(λ2) == legendre(N, 0:LMAX, 0:LMAX, x)
    for m in 0:LMAX
        # Check that vector filling function returns same answers
        legendre!(N, λ1, LMAX, m, x)
        @test @view(λ1[m+1:end]) == @view(λ2[m+1:end,m+1])
        for ℓ in m:LMAX
            # Then check that scalar output matches as well.
            legendre!(N, λ0, ℓ, m, x)
            @test λ0[] == λ1[ℓ+1]
        end
    end

    # As an implementation detail, the outputs arrays are not cleared outside of the valid
    # domain when operated up by the in-place functions, so initial or previous values
    # should remain. Check these invariants as a way to notice an implementation change.
    #
    # For the matrix output, the upper triangle (excluding diagonal) should be all NaN.
    @test isequal(triu(λ2, 1), triu(fill(T(NaN), LMAX+1, LMAX+1), 1))
    # The vector output should accumulate the (m,m) term at each step, i.e. the diagonal
    # of the matrix output.
    @test λ1 == λ2[diagind(λ2)]
end

# Test failed with world age problems when generated functions were incorrectly used;
# numerical derivative via newer world's Dual numbers tests for this error.
@testset "Derivatives via dual numbers" begin
    using ForwardDiff: derivative
    ctab = LegendreUnitCoeff{Float64}(LMAX)

    # Numerical derivative for Plm via recurrence relations:
    function num_deriv1(l, m, x)
        d = -(l+1) * x * Plm(l, m, x) + (l - m + 1) * Plm(l+1, m, x)
        return d / (x^2 - 1)
    end
    dual_deriv1(l, m, x) = derivative(z -> Plm(l, m, z), x)
    dual_deriv1_tab(l, m, x) = derivative(z -> ctab(l, m, z), x)

    for l in 0:LMAX, m in 0:l
        x = 2rand() - 1
        @test num_deriv1(l, m, x) ≈ dual_deriv1(l, m, x)
        @test num_deriv1(l, m, x) ≈ dual_deriv1_tab(l, m, x)
    end
end

@testset "Return inference ($(typeof(N)) with argument type $T)" for N in Norms, T in NumTypes
    Λ₀ = fill(0.0)
    Λ₁ = zeros(LMAX+1)
    Λ₂ = zeros(LMAX+1, LMAX+1)
    # 4 cases in recursion to handle:
    #   1. Initial condition
    #   2. Boost along diagonal
    #   3. Boost in ℓ from (m,m) -> (m+1, m)
    #   4. Boost again in ℓ from (m+1,m) -> (m+2, m) [and further]
    cases = ((0,0), (1,1), (2,1), (3,1))
    @testset "(ℓ,m) == ($ℓ, $m)" for (ℓ, m) in cases
        # Return type matches that of the argument for allocating interface, even though
        # internal calculation will promote.
        #
        # single (ℓ,m)
        @test @inferred(legendre(N, ℓ, m, one(T))) isa T
        @test @inferred(legendre!(N, Λ₀, LMAX, LMAX, one(T))) isa typeof(Λ₀)
        # all ℓ for fixed m
        @test @inferred(legendre(N, 0:ℓ, m, one(T))) isa Vector{T}
        @test @inferred(legendre!(N, Λ₁, LMAX, LMAX, one(T))) isa typeof(Λ₁)
        # all ℓ and m
        @test @inferred(legendre(N, 0:ℓ, 0:m, one(T))) isa Matrix{T}
        @test @inferred(legendre!(N, Λ₂, LMAX, LMAX, one(T))) isa typeof(Λ₂)
    end
end

@testset "Internal type promotion" begin
    # The implementation should promote all internal calculations, so given arguments
    # 0.5f0 (32-bit) and 0.5e0 (64-bit) which can both be exactly represented, calculate
    # that the former is calculated at 64-bit precision when the 64-bit output encourages
    # it.
    @test λlm!(fill(0.0), LMAX, LMAX, 0.5f0) == λlm!(fill(0.0), LMAX, LMAX, 0.5e0)
    # Do the same for a normalization which contains its own inherent element type.
    # Now take advantage of the fact that 0.2 is not exactly representable and have
    # different approximations at different accuracy to induce numerical differences
    # based on whether the values have been promoted or not.
    bigΛ! = LegendreSphereCoeff{BigFloat}(LMAX)
    λ1 = bigΛ!(fill(0f0), LMAX, LMAX, 0.2f0)
    λ2 = bigΛ!(fill(big(0f0)), LMAX, LMAX, big(0.2f0))
    @test λ1[] == Float32(λ2[])
end

@testset "Broadcasting" begin
    x = 0.5
    z = range(-1, 1, length=10)
    Λ = zeros(length(z), LMAX + 1, LMAX + 1)

    # Single (l,m)
    @test λlm.(LMAX, LMAX, x) == λlm(LMAX, LMAX, x)
    @test λlm.(LMAX, LMAX, z) == λlm(LMAX, LMAX, z)
    # All l for fixed m
    @test λlm.(0:LMAX, LMAX, x) == λlm(0:LMAX, LMAX, x)
    @test λlm.(0:LMAX, LMAX, z) == λlm(0:LMAX, LMAX, z)
    # All l and m over multiple z
    @test λlm.(0:LMAX, 0:LMAX, x) == λlm(0:LMAX, 0:LMAX, x)
    @test λlm.(0:LMAX, 0:LMAX, z) == λlm(0:LMAX, 0:LMAX, z)

    # Test that the single (l,m) for scalar and vector arguments inferrs correctly, since
    # the return is branched on whether the argument is 0-dimensional or not.
    # Dot syntax is not a call expression, so have to test the broadcasted() call directly.
    @test @inferred(Base.broadcasted(legendre, λlm, LMAX, LMAX, x)) isa Float64
    @test @inferred(Base.broadcasted(legendre, λlm, LMAX, LMAX, z)) isa Vector{Float64}
end

@testset "Shape preservation ($(ndims(x))-dimensional array)" for x in (fill(1.0), ones(1), ones(1,1))
    @test ndims(λlm(LMAX, LMAX, x)) == ndims(x)
    @test ndims(λlm(0:LMAX, LMAX, x)) == ndims(x) + 1
    @test ndims(λlm(0:LMAX, 0:LMAX, x)) == ndims(x) + 2

    @test ndims(λlm.(LMAX, LMAX, x)) == ndims(x)
    @test ndims(λlm.(0:LMAX, LMAX, x)) == ndims(x) + 1
    @test ndims(λlm.(0:LMAX, 0:LMAX, x)) == ndims(x) + 2
end

@testset "Axes preservation (offset axes $axs)" for axs in ((-5:5,), (-5:5, -5:5))
    using OffsetArrays
    using Base: OneTo

    X = reshape(collect(range(-1, 1, length=prod(length.(axs)))), axs...)

    @test axes(λlm(LMAX, LMAX, X)) == axs
    @test axes(λlm(0:LMAX, LMAX, X)) == (axs..., OneTo(LMAX+1))
    @test axes(λlm(0:LMAX, 0:LMAX, X)) == (axs..., OneTo(LMAX+1), OneTo(LMAX+1))

    @test axes(λlm.(LMAX, LMAX, X)) == axs
    @test axes(λlm.(0:LMAX, LMAX, X)) == (axs..., OneTo(LMAX+1))
    @test axes(λlm.(0:LMAX, 0:LMAX, X)) == (axs..., OneTo(LMAX+1), OneTo(LMAX+1))

    Λ  = λlm(0:LMAX, 0:LMAX, parent(X)) # Normal axes
    Λ′ = λlm(0:LMAX, 0:LMAX, X)         # Offset axes
    @test parent(Λ′) == Λ # Equality of values
    @test_throws DimensionMismatch λlm!(Λ, LMAX, LMAX, X) # Mismatched axes
end

@testset "Preallocated work space" begin
    using AssociatedLegendrePolynomials: unsafe_legendre!
    x1 = 0.5
    xN = range(-1, 1, length=10)
    Λ1 = zeros(LMAX+1, LMAX+1)
    ΛN = zeros(length(xN), LMAX+1, LMAX+1)
    norm = LegendreUnitNorm()
    work1 = AssociatedLegendrePolynomials.Work(norm, Λ1, x1)
    workN = AssociatedLegendrePolynomials.Work(norm, ΛN, xN)
    # Check equality before allocations to ensure the methods have been compiled.
    @test unsafe_legendre!(norm, copy(Λ1), LMAX, LMAX, x1) ==
            unsafe_legendre!(work1, Λ1, LMAX, LMAX, x1)
    @test unsafe_legendre!(norm, copy(ΛN), LMAX, LMAX, xN) ==
            unsafe_legendre!(workN, ΛN, LMAX, LMAX, xN)
    # Normal version should allocate at least as much space as 8 working buffers
    @test 8sizeof(work1.z) ≤ @allocated unsafe_legendre!(norm, Λ1, LMAX, LMAX, x1)
    @test 8sizeof(workN.z) ≤ @allocated unsafe_legendre!(norm, ΛN, LMAX, LMAX, xN)
    # The pre-allocated buffer space should then permit zero allocation calls
    @test 0 == @allocated unsafe_legendre!(work1, Λ1, LMAX, LMAX, x1)
    @test 0 == @allocated unsafe_legendre!(workN, ΛN, LMAX, LMAX, xN)
end


text:
# sino_geom.jl

using MIRT: sino_geom
using MIRT: sino_geom_help, sino_geom_plot_grids
using MIRT: image_geom, prompt
import MIRT: sino_geom_gamma_dfs

using Plots: plot!, plot, scatter, gui
using Test: @test, @test_throws, @inferred

"""
    sino_geom_show()
show an example of each of the 4 main geometries
"""
function sino_geom_show( ; kwarg...)
	down = 4
	ig = image_geom(nx=512, fov=500)
	ig = ig.down(down)

	arg = (nb = 888, na = down*8)
	nb_moj = round(Int, ig.nx*down*sqrt(2)) # match rfov
	sg_list = (
		sino_geom(:par ; d=ig.fovs[1]/arg.nb, arg..., kwarg...),
		sino_geom(:ge1 ; dfs=0, arg..., kwarg...), # arc
		sino_geom(:ge1 ; dfs=Inf, arg..., nb=960, kwarg...), # flat
		sino_geom(:moj ; d=1, arg..., nb=nb_moj, kwarg...),
		)

	nsg = length(sg_list)
	pl = Array{Any}(undef, nsg)
	for ii = 1:nsg
		sg = sg_list[ii].down(down)
		pl[ii] = plot(); sino_geom_plot!(sg, plot!)
		plot!(pl[ii], xlim=[-1,1]*550, ylim=[-1,1]*550)
		plot!(pl[ii], xtick=[-1,0,1]*250)
	end
	return pl
end

down = 8
ig = image_geom(nx=512, fov=500).down(down)
ig = image_geom(nx=ig.nx, dx=ig.dx, mask=ig.circ())

sg_list = (
	sino_geom(:par ; down=down, d=4),
	sino_geom(:ge1 ; down=down, orbit_start=20, dfs=0), # arc
	sino_geom(:ge1 ; down=down, orbit_start=20, dfs=Inf), # flat
	sino_geom(:moj ; down=down, d=4*sqrt(2)),
	sino_geom(:fan ; down=down, orbit=:short),
)

sg_list[1].help

nsg = length(sg_list)
pl = Array{Any}(undef, nsg)
for ii = 1:nsg
	sg = sg_list[ii]

	sg.ad[2]
	sg.rfov
	#@inferred
	sg.down(2)
	#@inferred
	sg.over(2)
	sg.dim
	sg.w
	sg.ones
	sg.zeros
	sg.dr
	sg.ds
	sg.r
	sg.s
	sg.ad
	sg.ar
	sg.xds
	sg.yds

	if sg isa SinoFan
		sg.gamma
		sg.gamma_max
		sg.dfs
		sg.dso
		sg.orbit_short
	end

	sg.grid
	sg.plot_grid(plot)
#	gui(); prompt()

	if sg isa SinoMoj
		sg.d_moj(0)
		sg.d_ang # angular dependent d for :moj
	end

	@test sg.shape(vec(sg.ones)) == sg.ones
	sg.taufun(ig.x, 0*ig.x)
	sg.unitv()
	sg.unitv(ib=1, ia=2)

	pl[ii] = plot(); sg.plot!(plot! ; ig=ig)
#	gui(); prompt()
end

#@inferred todo
@test sino_geom(:ge1 ; units=:cm) isa SinoFanArc

@test_throws String sino_geom(:badhow)
@test_throws String sino_geom(:ge1 ; dfs=-1)
@test_throws String sino_geom(:ge1 ; units=:bad)

sg = sino_geom(:ge1 ; orbit=:short)
@test sino_geom_gamma_dfs(sg) isa Vector

show(isinteractive() ? stdout : devnull, sg)
show(isinteractive() ? stdout : devnull, MIME("text/plain"), sg)

@inferred sino_geom_help()

pg = sino_geom_plot_grids(plot)
ps = sino_geom_show()

plot(pg..., ps..., layout=(2,4))
prompt()
plot(pl[1:4]...)