text:
import JUDI.get_computational_nt, JUDI.time_modeling

function get_computational_nt(srcGeometry::H5GeometryOOC, recGeometry::H5GeometryOOC, model::Model; dt=nothing)
  # Determine number of computational time steps
  nsrc = length(srcGeometry.container)
  nt = Array{Integer}(undef, nsrc)
  dtComp = calculate_dt(model; dt=dt)
  for j=1:nsrc
    dt_rec = recGeometry.dt[j]
    nt_rec = recGeometry.nt[j]
    dt_src = srcGeometry.dt[j]
    nt_src = srcGeometry.nt[j]
    ntRec = Int(dt_rec*(nt_rec-1) ÷ dtComp) + 1
    ntSrc = Int(dt_src*(nt_src-1) ÷ dtComp) + 1
    nt[j] = max(ntRec, ntSrc)
  end
  return nt
end

# Setup time-domain linear or nonlinear foward and adjoint modeling and interface to devito
function time_modeling(model_full::Model, srcGeometry::H5GeometryOOC, srcData::JUDI.ArrayOrNot,
                      recGeometry::H5GeometryOOC, recData::JUDI.ArrayOrNot,
                      dm::JUDI.PhysOrNot, op::Symbol, options::JUDIOptions)
  # Load full geometry for out-of-core geometry containers
  recGeometryIC = Geometry(recGeometry)
  srcGeometryIC = Geometry(srcGeometry)

  # Reutrn directly for J*0
  if op==:born
    if norm(dm) == 0 && options.return_array == false
      return judiVector(recGeometryIC, zeros(Float32, recGeometryIC.nt[1], length(recGeometryIC.xloc[1])))
    elseif norm(dm) == 0 && options.return_array == true
      return vec(zeros(Float32, recGeometryIC.nt[1], length(recGeometryIC.xloc[1])))
    end
  end

  # limit model to area with sources/receivers
  if options.limit_m == true
    model = deepcopy(model_full)
    model, dm = limit_model_to_receiver_area(srcGeometryIC, recGeometryIC, model, options.buffer_size; pert=dm)
  else
    model = model_full
  end

  # Set up Python model structure
  modelPy = devito_model(model, options, dm)

  # Remove receivers outside the modeling domain (otherwise leads to segmentation faults)
  recGeometryIC, recData = remove_out_of_bounds_receivers(recGeometryIC, recData, model)

  # Devito interface
  argout = devito_interface(modelPy, srcGeometryIC, srcData, recGeometryIC, recData, dm, options)
  # Extend gradient back to original model size
  if op==:adjoint_born && options.limit_m==true
    argout = extend_gradient(model_full, model, argout)
  end

  if options.save_data_to_disk
    nSamp = size(argout.data[1])[1]
    nTrc = size(argout.data[1])[2]

    (length(srcGeometry.container[1].ind) != nTrc) && (@error "Number of indexes doesn't match the number of traces\n")
    if @time !srcGeometry.container.seis.writeTrace(argout.data[1],srcGeometry.container[1].ind,0)
      @warn "Unable to write traces for pkey: $(srcGeometry.container[1].pkey)\n
            pkeyval: $(srcGeometry.container[1].pkeyvals[1])\n"
    end
    srcGeometry.container.seis.getH5File().flush()
  end

  return argout
end

# Saving to disk utilities
# save_to_disk(shot, srcGeometry, srcData, options, ::Val) = shot
# save_to_disk(shot::judiVector, srcGeometry, srcData, options, ::Val{false}) = shot

# function save_to_disk(data::)
#   shot.geometry[1]

#   container = write_shot_record(srcGeometry, srcData, shot.geometry[1], shot.data[1], options)
#   return judiVector(container)
# end

# function save_to_disk(shot::judiVector)
#   shot.geometry[1]

#   container = write_shot_record(srcGeometry, srcData, shot.geometry[1], shot.data[1], options)
#   return judiVector(container)
# end

# function write_shot_record(srcGeometry::GeometryIC, srcData, recGeometry::GeometryIC, recData, options)
#   q = judiVector(srcGeometry, srcData)
#   d = judiVector(recGeometry, recData)
#   pos = [srcGeometry.xloc[1][1], srcGeometry.yloc[1][1],  srcGeometry.zloc[1][1]]
#   pos = join(["_"*string(trunc(p; digits=2)) for p in pos])
#   file = join([string(options.file_name), pos,".segy"])
#   block_out = judiVector_to_SeisBlock(d, q)
#   segy_write(join([options.file_path,"/",file]), block_out)
#   container = scan_file(join([options.file_path,"/",file]),
#                         ["GroupX", "GroupY", "dt", "SourceSurfaceElevation", "RecGroupElevation"];
#                         chunksize=256)
#   return container
# end


text:
module Ntl

include("src/utils.jl")
include("src/models.jl")
include("src/samplers.jl")
include("src/fit.jl")
include("src/evaluate.jl")
include("src/plot.jl")
include("src/synthetic.jl")

end

ext: jl

max_stars_repo_name: realseanla/ntl-stick-breaking-julia

max_issues_repo_name: realseanla/ntl-stick-breaking-julia

max_forks_repo_name: realseanla/ntl-stick-breaking-julia


text:
using Documenter, FinEtools, FinEtoolsAcoustics

makedocs(
	modules = [FinEtoolsAcoustics],
	doctest = false, clean = true,
	format = Documenter.HTML(prettyurls = false),
	authors = "Petr Krysl",
	sitename = "FinEtoolsAcoustics.jl",
	pages = Any[
	"Home" => "index.md",
	"How to guide" => "guide/guide.md",
	"Types and Functions" => Any[
		"man/types.md",
		"man/functions.md"]
		]
	)

deploydocs(
    repo = "github.com/PetrKryslUCSD/FinEtoolsAcoustics.jl.git",
)


text:
@inline function add_dbfp_db(x::DoubleFloat{T}, y::T) where {T<:AbstractFloat}
    res = DoubleFloat{T}(add_ddfp_dd(HILO(x), y))
    return !isnan(HI(res)) ? res : add_dbfp_db_nonfinite(x,y)
end

@inline function sub_dbfp_db(x::DoubleFloat{T}, y::T) where {T<:AbstractFloat}
    res = DoubleFloat{T}(sub_ddfp_dd(HILO(x), y))
    return !isnan(HI(res))  ? res : sub_dbfp_db_nonfinite(x,y)
end

@inline function mul_dbfp_db(x::DoubleFloat{T}, y::T) where {T<:AbstractFloat}
    res = DoubleFloat{T}(mul_ddfp_dd(HILO(x), y))
    return !isnan(HI(res))  ? res : mul_dbfp_db_nonfinite(x,y)
end

@inline function dvi_dbfp_db(x::DoubleFloat{T}, y::T) where {T<:AbstractFloat}
    res = DoubleFloat{T}(dvi_ddfp_dd(HILO(x), y))
    return !isnan(HI(res)) ? res : dvi_dbfp_db_nonfinite(x,y)
end

@inline function add_dbfp_db_nonfinite(x::DoubleFloat{T}, y::T) where {T<:AbstractFloat}
    z = HI(x) + y
    return DoubleFloat{T}(z, zero(T))
end

@inline function sub_dbfp_db_nonfinite(x::DoubleFloat{T}, y::T) where {T<:AbstractFloat}
    z = HI(x) - y
    return DoubleFloat{T}(z, zero(T))
end

@inline function mul_dbfp_db_nonfinite(x::DoubleFloat{T}, y::T) where {T<:AbstractFloat}
    z = HI(x) * y
    return DoubleFloat{T}(z, zero(T))
end

@inline function dvi_dbfp_db_nonfinite(x::DoubleFloat{T}, y::T) where {T<:AbstractFloat}
    z = div(HI(x), y)
    return DoubleFloat{T}(z, zero(T))
end


text:
import Pkg
Pkg.rm(["Grids", "GridGames"])
Pkg.develop(["Grids", "GridGames"])
using MuZero
using Base.Threads: @spawn
include("game.jl")
include("muzero_config.jl")
config = getconfig(Tac)
actor_pool_size = Threads.nthreads() - 1
config.num_actors = actor_pool_size
@show Threads.nthreads(), config.num_actors
learner = Learner(
    config;
    uuid="b234fc40-a07b-4418-9fbc-adb0278dd4e1",
    dir="~/models/MuZeroModels.jl/tictactoe/",
    play_only=true,
)
play!(learner)


text:
#! /usr/bin/julia

# Rosetta Code, Terminal control/Ringing the terminal bell

println("This should ring a bell.\a")

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using WordTokenizers

@testset "POS" begin
    pos = PoSTagger()

    @testset "Basic" begin
        str = "The very first major corpus of English for computer analysis was the Brown Corpus."
        @test pos(str) ==  ["DT", "RB", "JJ", "JJ", "NN", "IN", "JJ", "IN", "NN", "NN", "VBD", "DT", "NNP", "NNP", "."]

        str = "If the Irish win the World Cup this year, it will be their 3rd time in a row."
        @test pos(str) == ["IN", "DT", "NNP", "VBP", "DT", "NNP", "NNP", "DT", "NN", ",", "PRP", "MD", "VB", "PRP\$", "CD", "JJ", "NN", "IN", "DT", "NN", "."]
    end

    @testset "Unknown Unicode characters" begin
        # Making sure that the pos model handles for unknown unicode characters
        str = "आ β⬰ 5¥ "
        @test length(pos(str)) == length(WordTokenizers.tokenize(str))

        str = "You owe John Doe 5¥."
        @test pos(str) ==  ["PRP", "VBP", "NNP", "NNP", "CD", "NNP", "."]
    end

    @testset "Documents and Corpus" begin
        pos = PoSTagger()

        text1 = "A little too small"
        text2 = "Here Foo Bar, please have some chocolate."

        sd = StringDocument(text1)
        td = TokenDocument(text2)

        tags = pos(sd)
        @test length(tags) == length(WordTokenizers.split_sentences(text1))
        @test length(tags[1]) == length(WordTokenizers.tokenize(text1))
        uniq1 = unique(vcat(tags...))
        @test "RB" ∈ uniq1
        @test "JJ" ∈ uniq1

        tags = pos(td)
        @test length(tags) == length(WordTokenizers.split_sentences(text2))
        @test length(tags[1]) == length(WordTokenizers.tokenize(text2))
        @test findall( x -> x == "NNP", vcat(tags...)) == [2, 3]

        crps = Corpus([sd, td])
        tags = pos(crps)

        @test length(tags) == length(crps.documents)
        @test tags[1] == pos(crps.documents[1])
        @test tags[2] == pos(crps.documents[2])
    end
end


text:
module CompilerWrapper

using JavaCall

include("PluginWrapper/PluginWrapper.jl")

end


text:
"""
    Unitcell

Struct containing the positions of basis sites, primitive translations and bonds for a lattice graph.
* `basis   :: Vector{SVector{3, Float64}}`       : position of basis sites in unitcell. basis[1] has to be the origin.
* `vectors :: Vector{SVector{3, Float64}}`       : primitive translations of the lattice
* `bonds   :: Vector{Vector{SVector{4, Int64}}}` : bonds connecting basis sites
Use `get_unitcell` to load the unitcell for a specific lattice and `lattice_avail` to print available lattices.
"""
struct Unitcell
    basis   :: Vector{SVector{3, Float64}}
    vectors :: Vector{SVector{3, Float64}}
    bonds   :: Vector{Vector{SVector{4, Int64}}}
end

# generate unitcell dummy
function get_unitcell_empty()

    basis   = Vector{SVector{3, Float64}}(undef, 1)
    vectors = Vector{SVector{3, Float64}}(undef, 1)
    bonds   = Vector{Vector{SVector{4, Int64}}}(undef, 1)
    uc      = Unitcell(basis, vectors, bonds)

    return uc
end

# load custom 2D unitcells
include("unitcell_lib/square.jl")
include("unitcell_lib/honeycomb.jl")
include("unitcell_lib/kagome.jl")
include("unitcell_lib/triangular.jl")

# load custom 3D unitcells
include("unitcell_lib/cubic.jl")
include("unitcell_lib/fcc.jl")
include("unitcell_lib/bcc.jl")
include("unitcell_lib/hyperhoneycomb.jl")
include("unitcell_lib/hyperkagome.jl")
include("unitcell_lib/pyrochlore.jl")
include("unitcell_lib/diamond.jl")

# print available lattices
function lattice_avail() :: Nothing

    println("###################### 2D Lattices ######################")
    println("square")
    println("honeycomb")
    println("kagome")
    println("triangular")
    println()
    println("###################### 3D Lattices ######################")
    println("cubic")
    println("fcc")
    println("bcc")
    println("hyperhoneycomb")
    println("hyperkagome")
    println("pyrochlore")
    println("diamond")

    return nothing
end

"""
    get_unitcell(
        name :: String
        )    :: Unitcell

Returns unitcell for lattice name. Use `lattice_avail` to print available lattices.
"""
function get_unitcell(
    name :: String
    )    :: Unitcell

    uc = get_unitcell_empty()

    if name == "square"
        uc = get_unitcell_square()
    elseif name == "honeycomb"
        uc = get_unitcell_honeycomb()
    elseif name == "kagome"
        uc = get_unitcell_kagome()
    elseif name == "triangular"
        uc = get_unitcell_triangular()
    elseif name == "cubic"
        uc = get_unitcell_cubic()
    elseif name == "fcc"
        uc = get_unitcell_fcc()
    elseif name == "bcc"
        uc = get_unitcell_bcc()
    elseif name == "hyperhoneycomb"
        uc = get_unitcell_hyperhoneycomb()
    elseif name == "hyperkagome"
        uc = get_unitcell_hyperkagome()
    elseif name == "pyrochlore"
        uc = get_unitcell_pyrochlore()
    elseif name == "diamond"
        uc = get_unitcell_diamond()
    else
        error("Unitcell $(name) unknown.")
    end
end


text:
#module Csv2mat
using MAT
using DataFrames
#export csv2mat
a = 13

"""
convert a csv table file to a .mat file for import to Matlab
csv2mat(filename, destination="log.mat"; startline=0, writeNames = true, subsample = 1, separator = ',')
"""
function csv2mat(filename, destination="log.mat"; startline=0, writeNames = true, subsample = 1, separator = ',')

    # Read the csv file into a dataframe
    println("Reading file $filename of size $(filesize(filename)/1_000_000) MB")
    df = readtable(filename, separator = separator, header=true, skipstart=startline)
    println("Preprocessing")
    (rows, cols) = size(df)
    # replace NA values with 0
    for  j = 1:cols, i = 1:rows
        ismissing(df[i,j]) && (df[i,j] = 0)
    end

    # Remove strange characters from variable names
    n = names(df)
    stripchars = [' ', '_', ':']
    for (i,symbol) in enumerate(n)
        n[i] = Symbol(strip(string(symbol),stripchars))
        writeNames && println(n[i])
    end
    names!(df,n) # Update the column names

    # Write the results to a .mat-file
    written = text = 0
    println("Creating save dict")
    d = Dict{String, Vector{<:Real}}()
    for i in 1:cols
        if !isa(df[1,i],Number)
            text += 1
            continue
        end
        written += 1
        a = replace(convert(Array,df[1:subsample:end,i]), missing=>0)
        d[string(n[i])] = a
    end
    println("Write the results to $destination")
    @time matwrite(destination, d)
    printstyled("Wrote $((round(Int,ceil(size(df,1)/subsample)),written)) entries to log.mat\n", color=:green)
    printstyled("Skipped $(text) columns with text\n", color=:yellow)

    # Perform requested plotting
    # if doplot
    #     colsToPlot = cols
    #     constants = Dict{Int64,Float64}()
    #     for i = 1:cols
    #         if all(df[:,i] .== df[1,i])
    #             typeof(df[1,i]) <: Number && push!(constants,i,df[1,i])
    #             colsToPlot-= 1
    #         end
    #         if !(typeof(df[1,i])<:Number)
    #             colsToPlot-= 1
    #         end
    #     end
    #     numberOfFigures = round(Int64,ceil(colsToPlot/subplots))
    #     print("Plotting all data in $(numberOfFigures) windows ... ")
    #     for window = 0:(numberOfFigures-1)
    #         plotrows::Int64 = ceil(sqrt(minimum([colsToPlot,subplots])))
    #         plotcols::Int64 = ceil(minimum([colsToPlot,subplots])/plotrows)
    #         doplot && plot(layout=(plotrows,plotcols))
    #         # Write the results to a .mat-file
    #         plotrow = plotcol = 0
    #         skip = 0;
    #
    #         for ii in (window*subplots).+collect(1:minimum([colsToPlot,subplots]))
    #             i = ii+skip
    #             while haskey(constants,i) || !(typeof(df[1,i])<:Number)
    #                 skip += 1
    #                 i = ii+skip
    #             end
    #             try
    #                 plot!(df[:,i], subplot=plotrow*plotcols + plotcol+1, title=replace(string(n[i]),"_"," "))
    #             catch ex
    #                 print_with_color(:red,"Could not plot column: $(string(n[i]))\n")
    #             end
    #             plotcol = (plotcol + 1) % plotcols
    #             if plotcol == 0
    #                 plotrow = (plotrow + 1) % plotrows
    #             end
    #         end
    #         #winston && display(t)
    #         colsToPlot -= plotrows*plotcols
    #     end
    #     println("done")
    #     print_with_color(:yellow,"Constants (not plotted):\n")
    #     for i in keys(constants)
    #         println(string(n[i], ":\t", round(constants[i],4)))
    #     end
    # end

end


text:
#=
for (fn) in (:ClCl, :ClOp, :OpCl, :OpOp)
    @eval ($fn){R1<:Real, R2<:Real}(a::R1, b::R2) = ($fn)(promote(a,b)...)
end
=#

# grasp
# it takes values presumed exact
# adjusting as needed for (from closedclosed into Grasp)
# so there is no risk of value exclusion
# or
# it takes values as grasped
# and adjusts as needed (from Grasp into closedclosed)
# so there is no risk of value exclusion
#
# grasped  yields a grasp_ed interval
# released yields a twice grasp_ed interval

grasp{R<:Real}(::Type{ClCl}, a::R, b::R) = (a,b)
grasp{R<:Real}(::Type{ClOp}, a::R, b::R) = (a,nextfloat(b))
grasp{R<:Real}(::Type{OpCl}, a::R, b::R) = (prevfloat(a),b)
grasp{R<:Real}(::Type{OpOp}, a::R, b::R) = widen(a,b)

grasped{G<:Grasp, R<:Real}(::Type{G}, a::R, b::R) = Rvl{G,R}( grasp(::Type{G},a,b)... )
release{G<:Grasp, R<:Real}(::Type{G}, a::R, b::R) = grasp(::Type{G}, a, b)
release{G<:Grasp, R<:Real}(this::Rvl{G,R})  = release(::Type{G}, values(this)... )
released{G<:Grasp, R<:Real}(this::Rvl{G,R}) = Rvl{G,R}( release(::Type{G},a,b)... )

if ::isdefined(Main,:ValidatedNumerics)

  import ValidatedNumerics.Interval, ValidatedNumerics.@interval

  convert{R<:Real}(::Type{Rvl{ClCl}}, x::Interval{R}) = ClCl(x.lo,x.hi)
  convert{R<:Real}(::Type{Rvl}, x::Interval{R}) = convert(Rvl{ClCl}, x)
  function  convert{G<:Grasp,R<:Real}(::Type{Rvl{G,R}, x::Interval{R})
      lo, hi = grasp(G, x.lo, x.hi)
      Rvl{G,R}(lo, hi)
  end
  function  convert{G<:Grasp,R<:Real}(::Type{Interval{R}}, x::Rvl{G,R})
      lo, hi = release(G, x.lo, x.hi)
      Interval(lo, hi)
  end

end


text:
module ParameterSpacePartitions
    using Requires, Distributions, ConcreteStructs, LinearAlgebra
    using ThreadedIterables, SpecialFunctions, ComponentArrays

    export find_partitions,
        adapt!,
        no_adaption!,
        estimate_volume,
        bias_correction,
        psp_slices

    export Options,
        Results

    include("init.jl")
    include("structs.jl")
    include("sampler.jl")
    include("volume.jl")
    include("TestModels.jl")

end

max_stars_repo_name: itsdfish/ParameterSpacePartitions

max_forks_repo_name: itsdfish/ParameterSpacePartitions


text:
#import OpenCL; const cl = OpenCL

function getPotentialKernel{T <: FloatingPoint}(:: Type{T})
        nType = T == Float64 ? "double" : "float"
        nPi = T == Float64 ? "M_PI" : "M_PI_F"
        nPi4 = T == Float64 ? "M_PI_4" : "M_PI_4_F"

        return "
        #if defined(cl_khr_fp64)  // Khronos extension available?
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
        #elif defined(cl_amd_fp64)  // AMD extension available?
        #pragma OPENCL EXTENSION cl_amd_fp64 : enable
        #endif
        #define number $nType
        #define number4 $(nType)4
        #define number8 $(nType)8
        #define PI $nPi
        #define PI_4 $nPi4

        #define Afield(x,y) a[y*D2 + x]
        #define Bfield(x,y) b[y*D2 + x]
        #define Area(x,y)   d[y*D2 + x]

        #define Apot(x,y) aout[y*D2 + x]
        #define Bpot(x,y) bout[y*D2 + x]

        static number sum(const number8 n) {
            return n.s0 + n.s1 +n.s2 + n.s3 + n.s4 + n.s5 + n.s6 + n.s7;
        }

        __kernel void potential(
                      __global const number *a,
                      __global const number *b,
                      __global const number4 *d,
                      __global number *aout,
                      __global number *bout,
                      const int D1,
                      const int D2,
                      const number repulsion) {

            int i = get_global_id(0);
            int j = get_global_id(1);

            int west = j-1;
            int east = j+1;
            int north = i+1;
            int south = i-1;

            if(j == 0)
                west = D2 - 1;
            if(j == D2 - 1)
                east = 0;
            if(i == D1 - 1)
                north = 0;
            if(i == 0)
                south = D1 - 1;

            const number4 a_ij = Area(i, j);

            const number selfRepulsion = 1*repulsion; // making this lower than repulsion allows for neighbors to have relative potential, so increases the chance that the hydrophobe will flow.
            const number r_ij = selfRepulsion * Afield(i,j);

            number8 area;
            area.s0 = Area(north, east).s0;
            area.s1 = Area(north, j   ).s1;
            area.s2 = Area(north, west).s2;
            area.s3 = Area(i    , west).s3;
            area.s4 = Area(south, west).s0;
            area.s5 = Area(south, j   ).s1;
            area.s6 = Area(south, east).s2;
            area.s7 = Area(i    , east).s3;

            number8 r_af;
            r_af.s0 = Afield(north, east);
            r_af.s1 = Afield(north, j   );
            r_af.s2 = Afield(north, west);
            r_af.s3 = Afield(i    , west);
            r_af.s4 = Afield(south, west);
            r_af.s5 = Afield(south, j   );
            r_af.s6 = Afield(south, east);
            r_af.s7 = Afield(i    , east);

            r_af = repulsion * r_af;

            Bpot(i,j) =   sum(r_af * area) + r_ij;

            Apot(i,j) = repulsion *(Bfield(north,east) * a_ij.s0 +
                                    Bfield(north,j   ) * a_ij.s1 +
                                    Bfield(north,west) * a_ij.s2 +
                                    Bfield(i,west    ) * a_ij.s3 +
                                    Bfield(south,west) * a_ij.s0 +
                                    Bfield(south,j   ) * a_ij.s1 +
                                    Bfield(south,east) * a_ij.s2 +
                                    Bfield(i,east    ) * a_ij.s3 + Bfield(i,j));
    }
"
end

function potentialCL!{T <: FloatingPoint}(
    a_buff, b_buff, area_buff,
    aout_buff, bout_buff,
    repulsion :: T, d1 :: Int64, d2 :: Int64,
    ctx, queue, program)

    k_p = Kernel(program, "potential")

    call(queue, k_p, (d1,d2), nothing, a_buff, b_buff, area_buff, aout_buff, bout_buff, int32(d1), int32(d2), repulsion)
end


text:
"""
 Tests for  https://link.springer.com/article/10.1007/s10543-007-0130-3 with test problems as in the paper.
 RS1, RS2
"""

import Statistics # for mean values of trajectories
import LinearAlgebra # for the normn
using StochasticDiffEq
using Test
using Random
#using DiffEqGPU

function generate_weak_solutions(prob, alg, dts, numtraj; ensemblealg=EnsembleThreads())
  sols = []
  for i in 1:length(dts)
    sol = solve(prob,alg;ensemblealg=ensemblealg,dt=dts[i],adaptive=false,save_start=false,save_everystep=false,weak_timeseries_errors=false,weak_dense_errors=false,trajectories=Int(numtraj))
    println(i)
    push!(sols,sol)
  end
  return sols
end

function strat_weak_errors(prob, alg, dts, numtraj; ensemblealg=EnsembleThreads())
  errors = Array{typeof(prob.prob.u0)}(undef, length((dts)), 2)
  for i in 1:length(dts)
    sol = solve(prob, alg, dt=dts[i], #save_start=false,save_everystep=true,weak_timeseries_errors=false,weak_dense_errors=false,
    trajectories=numtraj, adaptive=false)

    computed_exp = Statistics.mean([_sol.u[end] for _sol in sol])
    true_exp1 = prob.prob.u0*exp(1.0*(prob.prob.p[1]+0.5*prob.prob.p[2]^2))
    #true_exp1b = prob.u0*exp(1.0*(prob.p[1][1]))
    true_exp2 = Statistics.mean([_sol.u_analytic[end] for _sol in sol])
    println(i)
    errors[i,1] = LinearAlgebra.norm(computed_exp-true_exp1)
    errors[i,2] = LinearAlgebra.norm(computed_exp-true_exp2)
  end
  return errors
end

function prob_func(prob, i, repeat)
    remake(prob,seed=seeds[i])
end

"""
 Test Scalar SDEs (oop)
"""

@info "Scalar oop noise"

numtraj = Int(1e5) # in the paper they use 1e9
u₀ = 0.1
p = [1.5, 0.1]
f(u,p,t) = p[1]*u
g(u,p,t) = p[2]*u
dts = 1 .//2 .^(5:-1:1)
tspan = (0.0,1.0)

h1(z) = z

seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

prob = SDEProblem(f,g,u₀,tspan, p)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h1(sol[end]),false),
        prob_func = prob_func
        )
_solutions = @time generate_weak_solutions(ensemble_prob, RS1(), dts, numtraj, ensemblealg=EnsembleThreads())

errors = [LinearAlgebra.norm(Statistics.mean(sol.u) - u₀*exp(1.0*(p[1]+0.5*p[2]^2))) for sol in _solutions]
m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
@test -(m-2) < 0.3
println("RS1:", m)

#using Plots; convergence_plot = plot(dts, errors, xaxis=:log, yaxis=:log)
#savefig(convergence_plot, "RS1-final-"*string(numtraj)*".pdf")

# f_strat(u0,p,t,W) = @.(u0*exp(p[1]*t+p[2]*W))
# seed = 100
# Random.seed!(seed)
# seeds = rand(UInt, numtraj)
# prob_test = SDEProblem(SDEFunction(f,g,analytic=f_strat),g,u₀,tspan,p)
# ensemble_test = EnsembleProblem(prob_test;
#         #output_func = (sol,i) -> (h1(sol[end]),false),
#         prob_func = prob_func
#         )
# errors2 = @time strat_weak_errors(ensemble_test, RS1(), dts, numtraj, ensemblealg=EnsembleThreads())
#
# @test m ≈ log(errors2[end,1]/errors2[1,1])/log(dts[end]/dts[1])
# @test -(log(errors2[end,2]/errors2[1,2])/log(dts[end]/dts[1])-2) < 0.3

seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

_solutions = @time generate_weak_solutions(ensemble_prob, RS2(), dts, numtraj, ensemblealg=EnsembleThreads())

errors = [LinearAlgebra.norm(Statistics.mean(sol.u) - u₀*exp(1.0*(p[1]+0.5*p[2]^2))) for sol in _solutions]
m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
@test -(m-2) < 0.3

println("RS2:", m)

"""
 Test Scalar SDEs (iip)
"""

@info "Scalar iip noise"

f!(du,u,p,t) = du[1] = p[1]*u[1]
g!(du,u,p,t) = du[1] = p[2]*u[1]

seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

prob = SDEProblem(f!,g!,[u₀],tspan, p)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h1(sol[end]),false),
        prob_func = prob_func
        )
_solutions = @time generate_weak_solutions(ensemble_prob, RS1(), dts,
  numtraj, ensemblealg=EnsembleThreads())

errors = [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]+0.5*p[2]^2))) for sol in _solutions]
m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
@test -(m-2) < 0.3
println("RS1:", m)

seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

prob = SDEProblem(f!,g!,[u₀],tspan, p)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h1(sol[end]),false),
        prob_func = prob_func
        )
_solutions = @time generate_weak_solutions(ensemble_prob, RS2(), dts,
  numtraj, ensemblealg=EnsembleThreads())

errors = [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]+0.5*p[2]^2))) for sol in _solutions]
m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
@test -(m-2) < 0.3
println("RS1:", m)

"""
 Test non-commutative noise SDEs (iip)
"""

@info "Non-commutative noise"

u₀ = [0.1,0.1]
function f2!(du,u,p,t)
  du[1] = 5//4*u[2]-5//4*u[1]
  du[2] = 1//4*u[1]-1//4*u[2]
end
function g2!(du,u,p,t)
  du[1,1] = sqrt(3)/2*(u[1]-u[2])
  du[1,2] = 1//2*(u[1]+u[2])
  #du[2,1] = 0
  du[2,2] = u[1]
end
dts = 1 .//2 .^(4:-1:1)
tspan = (0.0,1.0)

h2(z) = z^2 # E(x_i) = 1/10 exp(1/2t) or E(x_1* x_2) = 1/100 exp(2t)

prob = SDEProblem(f2!,g2!,u₀,tspan,noise_rate_prototype=zeros(2,2))
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h2(sol[end][1]),false),
        prob_func = prob_func
        )

numtraj = Int(1e6)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

_solutions = @time generate_weak_solutions(ensemble_prob, RS1(), dts, numtraj, ensemblealg=EnsembleThreads())

errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(2)) for sol in _solutions]
m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
@test -(m-2) < 0.3

println("RS1:", m)

numtraj = Int(1e6)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

_solutions = @time generate_weak_solutions(ensemble_prob, RS2(), dts, numtraj, ensemblealg=EnsembleThreads())

errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(2)) for sol in _solutions]
m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
@test -(m-2) < 0.3

println("RS2:", m)

"""
 Test Diagonal noise SDEs (iip)
"""

@info "Diagonal noise"

u₀ = [0.1,0.1]
function f3!(du,u,p,t)
  du[1] = 299//200*u[1]
  du[2] = 299//200*u[2]
end
function g3!(du,u,p,t)
  du[1] = 1//10*u[1]
  du[2] = 1//10*u[1]
end
dts = 1 .//2 .^(5:-1:1)
tspan = (0.0,1.0)

h3(z) = z^2 # == 1//10**exp(3//2*t) if h3(z) = z and  == 1//100**exp(301//100*t) if h3(z) = z^2 )

prob = SDEProblem(f3!,g3!,u₀,tspan)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h3(sol[end][1]),false),
        prob_func = prob_func
        )

numtraj = Int(4e4)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

_solutions = @time generate_weak_solutions(ensemble_prob, RS1(), dts, numtraj, ensemblealg=EnsembleThreads())

errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
@test -(m-2) < 0.3

println("RS1:", m)

numtraj = Int(1e6)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

_solutions = @time generate_weak_solutions(ensemble_prob, RS2(), dts, numtraj, ensemblealg=EnsembleThreads())

errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
@test -(m-2) < 0.3

println("RS2:", m)


text:
#=
  M & M problem.

  From Think Bayes, page 6f.
  """
  M&M’s are small candy-coated chocolates that come in a variety of colors.
  Mars, Inc., which makes M&M’s, changes the mixture of colors from time
  to time.

  In 1995, they introduced blue M&M’s. Before then, the color mix in a bag
  of plain M&M’s was 30% Brown, 20% Yellow, 20% Red, 10% Green, 10%
  Orange, 10% Tan. Afterward it was 24% Blue , 20% Green, 16% Orange,
  14% Yellow, 13% Red, 13% Brown.

  Suppose a friend of mine has two bags of M&M’s, and he tells me that one
  is from 1994 and one from 1996. He won’t tell me which is which, but he
  gives me one M&M from each bag. One is yellow and one is green. What is
  the probability that the yellow one came from the 1994 bag?

  """

  Page 7 (the table): The answer is 20/27: ~0.74074.

   Distributions of variable mix0
   mix1994    =>    7442  (0.744200)
   mix1996    =>    2558  (0.255800)

   Distributions of variable mix1
   mix1996    =>    7442  (0.744200)
   mix1994    =>    2558  (0.255800)

   Distributions of variable color0
   yellow     =>   10000  (1.000000)

   Distributions of variable color1
   green      =>   10000  (1.000000)

  Cf ~/webppl/m_and_m_problem.wppl

=#

using Turing, StatsPlots

include("jl_utils.jl")

@model function m_and_m_problem()
    brown  = 1
    yellow = 2
    red    = 3
    green  = 4
    orange = 5
    tan    = 6
    blue   = 7
    colors = [brown,yellow,red,green,orange,tan,blue]

    mix1994 = 1
    mix1996 = 2
    mixes = [mix1994, mix1996]

    # First pick a bag in mix0 and then pick the other bag in mix1
    mix0 ~ UniformDraw(mixes)
    mix1 ~ Dirac(mix0 === mix1994 ? mix1996 : mix1994)

    # The mixes of each colors for the two years
    mix1994ps = [30,20,20,10,10,10,0]
    mix1996ps = [13,14,13,20,16,0,24]
    color0 ~ mix0 == mix1994 ? DiscreteNonParametric(colors,simplex(mix1994ps)) :
                               DiscreteNonParametric(colors,simplex(mix1996ps))

    color1 ~ mix1 == mix1994 ? DiscreteNonParametric(colors,simplex(mix1994ps)) :
                               DiscreteNonParametric(colors,simplex(mix1996ps))

    true ~ Dirac(color0 == yellow)
    true ~ Dirac(color1 == green)

end

model = m_and_m_problem()

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(15), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns,:mix0,[:mix1994, :mix1996])
show_var_dist_pct(chns,:mix1,[:mix1994, :mix1996])

println()
show_var_dist_pct(chns,:color0,["brown","yellow","red","green","orange","tan","blue"])
show_var_dist_pct(chns,:color1,["brown","yellow","red","green","orange","tan","blue"])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
module NfMatModule

using Hecke
import Nemo
import Hecke: is_zero_entry
"""
In Antic, `nf_elem` is a union of
-  Degree 1 Field:
   - fmpz  numerator
   - fmpz  denominator
-  Degree 2 Field:
   - fmpz [3] coeffs 0, 1, and 2 to have a buffer for products
   - fmpz  denominator
-  Degree > 2 Field:
   - fmpq_poly which is
     - fmpz * numerator coeffs
     - fmpz   denominator
     - slong  alloc
     - slong  length
All in all, at most 4 Ptr/long/Int. This is "modelled" here as
4 Ints.

Important: the denominator has to be initioalised in all cases, a
binary zero is fine for the other fields.
Due to the different position of the denominator in deg 2, we have
2 constructors below.
Further `zero` is constructing a valid 0 in degree 1 and > 2
while `one` is used to build a zero in degree 2.

Given and `nf_elem_raw` are immutable (otherwise the memory layout does not
work, ie. a flat array) and thus cannot have a finalizer, this type can only
be used inside a larger sturcture via `Vector{nf_elem_raw}` as in
the matrices. This would also allow polynomials...
"""
struct nf_elem_raw
  a::Int
  b::Int
  c::Int
  d::Int
  function nf_elem_raw()
    new(0,1,0,0)  #the den is in pos 2 for deg == 1 or > 2
                  #and in pos 4 for deg == 2
  end
  function nf_elem_raw(::Int)
    new(0,0,0,1)  #the den is in pos 2 for deg == 1 or > 2
                  #and in pos 4 for deg == 2
  end
end

function Base.zero(::Type{nf_elem_raw})
  return nf_elem_raw()
end

function Base.one(::Type{nf_elem_raw})
  return nf_elem_raw(1)
end

mutable struct NfMatElem <: MatElem{nf_elem}
  entries:: Vector{nf_elem_raw}
  rows:: Vector{Int}
  nrows::Int
  ncols::Int
  base_ring::AnticNumberField

  function NfMatElem(K::AnticNumberField, r::Int, c::Int)
    if degree(K) == 2
      entries = ones(nf_elem_raw, r*c)
    else
      entries = zeros(nf_elem_raw, r*c)
    end
    rows = zeros(Int, r)
    for i=1:r
      rows[i] = (i-1)*c
    end
    M = new(entries, rows, r, c, K)
    #= to support `view`s, the finalizer is associated
       to the `entries` *not* the matrix.
       Since the finalizer get only array and not the field,
       different finalizer have to be used depending on the degree
    =#
    #=
    for i=1:r
      for j=1:c
        init!(getindex_raw(M, i, j), K)
      end
    end
    =#
    if degree(K) == 1
      finalizer(NfMatElem_clear1, entries)
    elseif degree(K) == 2
      finalizer(NfMatElem_clear2, entries)
    else
      finalizer(NfMatElem_clear3, entries)
    end
    return M
  end
  #creates a view M[r1:r2, c1:c2]
  function NfMatElem(M::NfMatElem, r1::Int, r2::Int, c1::Int, c2::Int)
    rows = M.rows[r1:r2]
    for i=1:length(rows)
      rows[i] += c1-1
    end
    return new(M.entries, rows, r2-r1+1, c2-c1+1, base_ring(M))
  end
  function NfMatElem(M::NfMatElem)
    K = base_ring(M)
    if degree(K) == 2
      entries = ones(nf_elem_raw, length(M.entries))
    else
      entries = zeros(nf_elem_raw, length(M.entries))
    end
    for i=1:length(entries)
      if M.entries[i].a != 0 #wrong for deg2 , here .d needs to be tested
        ccall((:nf_elem_set, Nemo.libantic), Cvoid, (Ptr{nf_elem_raw}, Ptr{nf_elem_raw}, Ref{AnticNumberField}), pointer(entries, i), pointer(M.entries, i), K)
      end
    end
    if degree(K) == 1
      finalizer(NfMatElem_clear1, entries)
    elseif degree(K) == 2
      finalizer(NfMatElem_clear2, entries)
    else
      finalizer(NfMatElem_clear3, entries)
    end
    return new(entries, copy(M.rows), nrows(M), ncols(M), K)
  end
end

function NfMatElem_clear1(en::Vector{nf_elem_raw})
  #degree 1 case...
  for i=1:length(en)
    p = pointer(en, i)
    ccall((:fmpz_clear, Nemo.libflint), Cvoid, (Ptr{nf_elem_raw}, ), p)
    p += sizeof(Int)
    ccall((:fmpz_clear, Nemo.libflint), Cvoid, (Ptr{nf_elem_raw}, ), p)
  end
end

function NfMatElem_clear2(en::Vector{nf_elem_raw})
  for i=1:length(en)
    p = pointer(en, i)
    ccall((:fmpz_clear, Nemo.libflint), Cvoid, (Ptr{nf_elem_raw}, ), p)
    p += sizeof(Int)
    ccall((:fmpz_clear, Nemo.libflint), Cvoid, (Ptr{nf_elem_raw}, ), p)
    p += sizeof(Int)
    ccall((:fmpz_clear, Nemo.libflint), Cvoid, (Ptr{nf_elem_raw}, ), p)
    p += sizeof(Int)
    ccall((:fmpz_clear, Nemo.libflint), Cvoid, (Ptr{nf_elem_raw}, ), p)
  end
end

function NfMatElem_clear3(en::Vector{nf_elem_raw})
  for i=1:length(en)
    ccall((:fmpq_poly_clear, Nemo.libflint), Cvoid, (Ptr{nf_elem_raw},), pointer(en, i))
  end
end

function Base.deepcopy_internal(M::NfMatElem, d::IdDict)
  return NfMatElem(M::NfMatElem)
end

function Base.checkbounds(::Type{Bool}, M::NfMatElem, i::Int, j::Int)
  return 1 <= i <= M.nrows && 1 <= j <= M.ncols
end

function Base.checkbounds(M::NfMatElem, i::Int, j::Int)
  Base.checkbounds(Bool, M, i, j) || throw(BoundsError(M, (i,j)))
end

Hecke.nrows(M::NfMatElem) = M.nrows
Hecke.ncols(M::NfMatElem) = M.ncols
Hecke.base_ring(M::NfMatElem) = M.base_ring

@inline rows(M::NfMatElem) = M.rows
@inline ent(M::NfMatElem) = M.entries

####
#
# view
#
####
function Base.view(M::NfMatElem, r1::Int , r2::Int, c1::Int, c2::Int)
  @assert 0 < r1 <= r2 <= nrows(M)
  @assert 0 < c1 <= c2 <= ncols(M)
  return NfMatElem(M, r1, r2, c1, c2)
end

function Base.view(M::NfMatElem, r::UnitRange{Int}, c::UnitRange{Int})
  return Base.view(M, r.start, r.stop, c.start, c.stop)
end

function Base.view(M::NfMatElem, r::Colon, c::UnitRange{Int})
  return Base.view(M, 1:nrows(M), c)
end

function Base.view(M::NfMatElem, r::UnitRange{Int}, c::Colon)
  return Base.view(M, r, 1:ncols(M))
end

####
#
# Creation/ conversion
#
####
function Hecke.matrix(K::AnticNumberField, r::Int, c::Int, a::Vector{nf_elem})
  M = NfMatElem(K, r, c)
  for i=1:r
    for j=1:c
      M[i,j] = a[(i-1)*c+j]
    end
  end
  return M
end

#convert from "classical" implementation of MatElem{nf_elem}
function Hecke.matrix(M::Generic.MatSpaceElem{nf_elem})
  K = base_ring(M)
  N = NfMatElem(K, nrows(M), ncols(M))
  for i=1:nrows(M)
    for j=1:ncols(M)
      N[i,j] = M[i,j]
    end
  end
  return N
end

function Hecke.zero_matrix(K::AnticNumberField, r::Int, c::Int)
  return NfMatElem(K, r, c)
end

function Hecke.identity_matrix(K::AnticNumberField, n::Int)
  M = NfMatElem(K, n, n)
  for i=1:n
    M[i,i] = one(K)
  end
  return M
end

function Hecke.identity_matrix(M::NfMatElem)
  return identity_matrix(base_ring(M), nrows(M))
end

function Array(M::NfMatElem)
  return collect(M)
end

function Generic.Mat{nf_elem}(M::NfMatElem)
  N =  Generic.MatSpaceElem{nf_elem}(Array(M))
  N.base_ring = M.base_ring
  return N
end

Hecke.dense_matrix_type(::Type{nf_elem}) = NfMatElem

####
#
# some concatenation
#
####
function Base.vcat(M::NfMatElem, N::NfMatElem)
  K = base_ring(M)
  @assert K == base_ring(N)
  @assert ncols(M) == ncols(N)
  MN = NfMatElem(base_ring(M), nrows(M)+nrows(N), ncols(M))
  for i=1:nrows(M)
    for j=1:ncols(M)
      MN[i,j] = getindex_raw(M, i, j)
    end
  end
  for i=1:nrows(N)
    for j=1:ncols(N)
      MN[i+nrows(M), j] = getindex_raw(N, i, j)
    end
  end
  return MN
end

function Hecke.vcat!(M::NfMatElem, N::NfMatElem)
  @assert ncols(M) == ncols(N)
  n = nrows(M)
  M.nrows += nrows(N)
  K = base_ring(M)
  @assert base_ring(N) == base_ring(M)
  if degree(K) == 2
    append!(M.entries, ones(nf_elem_raw,  length(N.entries)))
  else
    append!(M.entries, zeros(nf_elem_raw,  length(N.entries)))
  end
  for i=1:nrows(N)
    push!(M.rows, (i+n-1)*ncols(M))
    for j=1:ncols(M)
      if !is_zero_entry(N, i, j)
        M[i+n, j] = getindex_raw(N, i, j)
      end
    end
  end
end

function Base.hcat(M::NfMatElem, N::NfMatElem)
  K = base_ring(M)
  @assert K == base_ring(N)
  @assert nrows(M) == nrows(N)
  MN = NfMatElem(base_ring(M), nrows(M), ncols(M)+ncols(N))
  for i=1:nrows(M)
    for j=1:ncols(M)
      MN[i,j] = getindex_raw(M, i, j)
    end
    for j=1:ncols(N)
      MN[i, j+ncols(M)] = getindex_raw(N, i, j)
    end
  end
  return MN
end

####
#
# Access - the checkbounds can be switched of for the raw versino
# the generic versino does alloc and copy, thus is slow anyway
#
####
function Base.getindex(M::NfMatElem, r::Int, c::Int)
  checkbounds(M, r, c)
  a = base_ring(M)()
  ccall((:nf_elem_set, Nemo.libantic), Cvoid, (Ref{nf_elem}, Ptr{nf_elem_raw}, Ref{AnticNumberField}), a, pointer(M.entries, M.rows[r]+c), base_ring(M))
  return a
end

function Base.getindex(M::NfMatElem, r::AbstractVector{Int}, c::AbstractVector{Int})
  N = NfMatElem(base_ring(M), length(r), length(c))
  ii = 1
  for i = r
    jj = 1
    for j = c
      N[ii, jj] = M[i,j]
      jj += 1
    end
    ii += 1
  end
  return N
end

_rrange(::NfMatElem, r::Int) = r:r
_crange(::NfMatElem, c::Int) = c:c
_rrange(M::NfMatElem, r::Colon) = 1:nrows(M)
_crange(M::NfMatElem, c::Colon) = 1:ncols(M)
_rrange(::NfMatElem, r::UnitRange{Int}) = r
_crange(::NfMatElem, c::UnitRange{Int}) = c
_rrange(::NfMatElem, r::Vector{Int}) = r
_crange(::NfMatElem, c::Vector{Int}) = c

function Base.getindex(M::NfMatElem, r::Union{Colon, Int64, AbstractVector{Int64}}, c::Union{Colon, Int64, AbstractVector{Int64}})
  return Base.getindex(M, _rrange(M, r), _crange(M, c))
end

@inline function getindex_raw(M::NfMatElem, r::Int, c::Int)
  @boundscheck checkbounds(M, r, c)
  en = ent(M)
  ro = rows(M)
  @assert ro[r]+c <= length(en)
  return @inbounds pointer(en, ro[r]+c)
end

function Base.setindex!(M::NfMatElem, a::nf_elem, r::Int, c::Int)
  checkbounds(M, r, c)
  ccall((:nf_elem_set, Nemo.libantic), Cvoid, (Ptr{nf_elem_raw}, Ref{nf_elem}, Ref{AnticNumberField}), pointer(M.entries, M.rows[r]+c), a, base_ring(M))
end

function Base.setindex!(M::NfMatElem, a::Int, r::Int, c::Int)
  checkbounds(M, r, c)
  ccall((:nf_elem_set_si, Nemo.libantic), Cvoid, (Ptr{nf_elem_raw}, Int, Ref{AnticNumberField}), pointer(M.entries, M.rows[r]+c), a, base_ring(M))
end

function Base.setindex!(M::NfMatElem, a::fmpz, r::Int, c::Int)
  checkbounds(M, r, c)
  ccall((:nf_elem_set_fmpz, Nemo.libantic), Cvoid, (Ptr{nf_elem_raw}, Ref{fmpz}, Ref{AnticNumberField}), pointer(M.entries, M.rows[r]+c), a, base_ring(M))
end

#hopefully not used?
function Base.setindex!(M::NfMatElem, a::nf_elem_raw, r::Int, c::Int)
  @boundscheck checkbounds(M, r, c)
  ccall((:nf_elem_set, Nemo.libantic), Cvoid, (Ptr{nf_elem_raw}, Ref{nf_elem_raw}, Ref{AnticNumberField}), pointer(M.entries, M.rows[r]+c), a, base_ring(M))
end

@inline function Base.setindex!(M::NfMatElem, a::Ptr{nf_elem_raw}, r::Int, c::Int)
  @boundscheck checkbounds(M, r, c)
  ccall((:nf_elem_set, Nemo.libantic), Cvoid, (Ptr{nf_elem_raw}, Ref{nf_elem_raw}, Ref{AnticNumberField}), pointer(M.entries, M.rows[r]+c), a, base_ring(M))
end

function Base.setindex!(M::NfMatElem, N::NfMatElem, r::UnitRange{Int}, c::UnitRange{Int})
  for i = r
    for j = c
      M[i,j] = N[i-r.start+1,j-c.start+1] #getindex_raw(N, i, j)
    end
  end
end

#= yields ambiguities
function Base.setindex!(M::NfMatElem, N::NfMatElem, r::Union{Colon, Int64, AbstractVector{Int64}}, c::Union{Colon, Int64, AbstractVector{Int64}})
  Base.setindex!(M, N, _rrange(M, r), _crange(M, c))
end
=#

function Hecke.swap_rows!(M::NfMatElem, i::Int, j::Int)
  M.rows[i], M.rows[j] = M.rows[j], M.rows[i]
end

function Hecke.swap_rows(M::NfMatElem, i::Int, j::Int)
  N = deepcopy(M)
  swap_rows!(N, i, j)
  return N
end

@inline function swap_entry!(M::NfMatElem, r1::Int, c1::Int, r2::Int, c2::Int)
  @boundscheck checkbounds(M, r1, c1)
  @boundscheck checkbounds(M, r2, c2)
  M.entries[M.rows[r1]+c1], M.entries[M.rows[r2]+c2] =
     M.entries[M.rows[r2]+c2], M.entries[M.rows[r1]+c1]
end

function Hecke.swap_cols!(M::NfMatElem, i::Int, j::Int)
  for r=1:nrows(M)
    swap_entry!(M, r, i, r, j)
  end
end

function Hecke.swap_cols(M::NfMatElem, i::Int, j::Int)
  N = deepcopy(M)
  swap_cols!(N, i, j)
  return N
end
#=
function Base.show(io::IO, M::NfMatElem)
  print(io, "$(M.nrows) x $(M.ncols) matrix over $(M.base_ring)")
end
=#

@inline function Hecke.is_zero_entry(M::NfMatElem, i::Int, j::Int)
  p = getindex_raw(M, i, j)
  reduce!(p, base_ring(M))
  return ccall((:nf_elem_is_zero, Nemo.libantic), Cint, (Ptr{nf_elem_raw}, Ref{AnticNumberField}), p, base_ring(M)) == 1
end

function Hecke.is_zero_row(M::NfMatElem, i::Int)
  return all(x->Hecke.is_zero_entry(M, i, x), 1:ncols(M))
end

function Hecke.is_zero_column(M::NfMatElem, i::Int)
  return all(x->Hecke.is_zero_entry(M, x, i), 1:nrows(M))
end

@inline function isone_entry(M::NfMatElem, i::Int, j::Int)
  p = getindex_raw(M, i, j)
  reduce!(p, base_ring(M))
  return ccall((:nf_elem_is_one, Nemo.libantic), Cint, (Ptr{nf_elem_raw}, Ref{AnticNumberField}), p, base_ring(M)) == 1
end

@inline function mul!(a::Ptr{nf_elem_raw}, b::Ptr{nf_elem_raw}, c::Ptr{nf_elem_raw}, K::AnticNumberField, red::Bool = true)
  ccall((:nf_elem_mul_red, Nemo.libantic), Cvoid, (Ptr{nf_elem_raw}, Ptr{nf_elem_raw}, Ptr{nf_elem_raw}, Ref{AnticNumberField}, Cint), a, b, c, K, red)
end

@inline function mul!(a::Ptr{nf_elem_raw}, b::Ptr{nf_elem_raw}, c::nf_elem, K::AnticNumberField, red::Bool = true)
  ccall((:nf_elem_mul_red, Nemo.libantic), Cvoid, (Ptr{nf_elem_raw}, Ptr{nf_elem_raw}, Ref{nf_elem}, Ref{AnticNumberField}, Cint), a, b, c, K, red)
end

@inline function mul!(a::nf_elem, b::Ptr{nf_elem_raw}, c::Ptr{nf_elem_raw}, K::AnticNumberField, red::Bool = true)
  ccall((:nf_elem_mul_red, Nemo.libantic), Cvoid, (Ref{nf_elem}, Ptr{nf_elem_raw}, Ptr{nf_elem_raw}, Ref{AnticNumberField}, Cint), a, b, c, K, red)
end

@inline function mul!(a::nf_elem, b::nf_elem, c::Ptr{nf_elem_raw}, K::AnticNumberField, red::Bool = true)
  ccall((:nf_elem_mul_red, Nemo.libantic), Cvoid, (Ref{nf_elem}, Ref{nf_elem}, Ptr{nf_elem_raw}, Ref{AnticNumberField}, Cint), a, b, c, K, red)
end

@inline function reduce!(a::Ptr{nf_elem_raw}, K::AnticNumberField)
  ccall((:nf_elem_reduce, Nemo.libantic), Cvoid, (Ptr{nf_elem_raw}, Ref{AnticNumberField}), a, K)
end

@inline function add!(a::Ptr{nf_elem_raw}, b::Ptr{nf_elem_raw}, c::Ptr{nf_elem_raw}, K::AnticNumberField)
  ccall((:nf_elem_add, Nemo.libantic), Cvoid, (Ptr{nf_elem_raw}, Ptr{nf_elem_raw}, Ptr{nf_elem_raw}, Ref{AnticNumberField}), a, b, c, K)
end

@inline function add!(a::Ptr{nf_elem_raw}, b::Ptr{nf_elem_raw}, c::nf_elem, K::AnticNumberField)
  ccall((:nf_elem_add, Nemo.libantic), Cvoid, (Ptr{nf_elem_raw}, Ptr{nf_elem_raw}, Ref{nf_elem}, Ref{AnticNumberField}), a, b, c, K)
end

@inline function sub!(a::Ptr{nf_elem_raw}, b::Ptr{nf_elem_raw}, c::Ptr{nf_elem_raw}, K::AnticNumberField)
  ccall((:nf_elem_sub, Nemo.libantic), Cvoid, (Ptr{nf_elem_raw}, Ptr{nf_elem_raw}, Ptr{nf_elem_raw}, Ref{AnticNumberField}), a, b, c, K)
end

@inline function sub!(a::Ptr{nf_elem_raw}, b::Ptr{nf_elem_raw}, c::nf_elem, K::AnticNumberField)
  ccall((:nf_elem_sub, Nemo.libantic), Cvoid, (Ptr{nf_elem_raw}, Ptr{nf_elem_raw}, Ref{nf_elem}, Ref{AnticNumberField}), a, b, c, K)
end

function Hecke.divide_row!(M::NfMatElem, i::Int, a::nf_elem)
  ai = inv(a)
  K = parent(a)
  @assert base_ring(M) === K
  for j=1:ncols(M)
    mul!(getindex_raw(M, i, j), getindex_raw(M, i, j), a, K)
  end
end

function Hecke.transpose(M::NfMatElem)
  N = zero_matrix(base_ring(M), ncols(M), nrows(M))
  for i=1:nrows(M)
    for j=1:ncols(M)
      N[j,i] = getindex_raw(M, i, j)
    end
  end
  return N
end

function reduce!(M::NfMatElem, q::NfMatElem, piv::Vector{Int})
  K = base_ring(M)
  if !true

    @assert length(piv) == ncols(q)
    for i=1:length(piv)
      if piv[i] != 0
        @assert q[piv[i], i] == 1
        @assert all(j->(piv[i] == j) || (q[j, i] == 0), 1:nrows(q))
      end
    end

  a = zero_matrix(K, nrows(M), nrows(q))
  for l=1:length(piv)
    if piv[l] == 0
      continue
    end
    for i=1:nrows(M)
      a[i, piv[l]] = getindex_raw(M, i, l)
    end
  end
  N = mul_KS(a, q)
  #=
  @assert N == a*q
  for l=1:length(piv)
    if piv[l] == 0
      continue
    end
    for i=1:nrows(M)
      @assert N[i, l] == M[i, l]
    end
  end
  T = M-N
  =#

  Nemo.sub!(M, M, N)
  #=
  @assert M == T

  for l=1:length(piv)
    if piv[l] == 0
      continue
    end
    for i=1:nrows(M)
      if !is_zero_entry(M, i, l)
      end
      @assert is_zero_entry(M, i, l)
    end
  end
  =#

  return
end

  t = K()
  for i=1:nrows(M)
    for l=1:length(piv)
      if piv[l] == 0
        continue
      end
      reduce!(getindex_raw(M, i, l), K)
      if is_zero_entry(M, i, l)
        continue
      end
      #M[i, :] -= q[piv[l], :]*M[i, j]
      for k=1:ncols(M)
        if k == l
          continue
        end
        mul!(t, getindex_raw(q, piv[l], k), getindex_raw(M, i, l), K, false)
        sub!(getindex_raw(M, i, k), getindex_raw(M, i, k), t, K)
      end
      zero!(getindex_raw(M, i, l), K)
    end
  end
  reduce!(M)
end

#standart Gauss, no brain.
# assumes 1:start is already in ref with pivot array in piv
# operates in start:stop
#  start:stop is reduced modulo 1:start-1
#  then Gauss continues
#
# if det == true then the scaling of multiplied together to find the
# determinant.
function _ref!(M::NfMatElem;
   piv  ::Vector{Int} = zeros(Int, ncols(M)),
   start::Int = 1, stop::Int = nrows(M),
   det  ::Bool = false,
   trafo::Union{Nothing, NfMatElem} = nothing)

  K = base_ring(M)
  t = K()
  de = K(1)
  for i=1:stop
    if i < start
      j = findfirst(isequal(i), piv)
      @assert j !== nothing
    else
      j = 1
      while j <= ncols(M) && Hecke.is_zero_entry(M, i, j)
        j += 1
      end
      if j>ncols(M)
        continue
      end
      @inbounds piv[j] = i
      @assert !Hecke.is_zero_entry(M, i, j)
      if det
        Nemo.mul!(de, de, M[i,j])
      end
      if !isone_entry(M, i, j)
        ccall((:nf_elem_inv, Nemo.libantic), Cvoid, (Ref{nf_elem}, Ptr{nf_elem_raw}, Ref{AnticNumberField}), t, getindex_raw(M, i, j), K)
        k = j
        while k<= ncols(M)
          reduce!(getindex_raw(M, i, k), K)
          mul!(getindex_raw(M, i, k), getindex_raw(M, i, k), t, K)
          k += 1
        end
        if trafo !== nothing
          for k=1:ncols(trafo)
            mul!(getindex_raw(trafo, i, k), getindex_raw(trafo, i, k), t, K)
          end
        end
      end
    end
    @assert isone_entry(M, i, j)
    for r = max(start, i+1):nrows(M)
      if Hecke.is_zero_entry(M, r, j)
        continue
      end
      s = getindex_raw(M, r, j)
      reduce!(s, K)
      if trafo !== nothing
        for k=1:ncols(trafo)
          if is_zero_entry(trafo, i, k)
            continue
          end
          Mrk = getindex_raw(trafo, r, k)
          mul!(t, s, getindex_raw(trafo, i, k), K, false)
          sub!(Mrk, Mrk, t, K)
          reduce!(Mrk, K)
        end
      end
      for k=ncols(M):-1:j
        if Hecke.is_zero_entry(M, i, k)
          continue
        end
        Mrk = getindex_raw(M, r, k)
        mul!(t, s, getindex_raw(M, i, k), K, false)
        sub!(Mrk, Mrk, t, K)
      end
    end
  end
  reduce!(M)
  return piv, de
end

function ref!(M::NfMatElem)
  piv, _ = _ref!(M)
  permute_rows!(M, piv)
  return length(piv) - length(findall(iszero, piv))
end

function ref(M::NfMatElem)
  N = deepcopy(M)
  return ref!(N), N
end

Nemo.iszero(a::Ptr{nf_elem_raw}, K::AnticNumberField) =
  ccall((:nf_elem_is_zero, Nemo.libantic), Int, (Ptr{nf_elem_raw}, Ref{AnticNumberField}), a, K) == 1

function reduce_up!(M::NfMatElem, piv::Vector{Int})
  K = base_ring(M)
  t = K()
  for j = 1:ncols(M)
    @inbounds p = piv[j]
    iszero(p) && continue
    for k=1:p-1
      s = getindex_raw(M, k, j)
      if iszero(s, K)
        continue
      end
      for h=ncols(M):-1:j
        mul!(t, getindex_raw(M, p, h), s, K)
        sub!(getindex_raw(M, k, h), getindex_raw(M, k, h), t, K)
      end
    end
  end
end

function Hecke.rref!(M::NfMatElem)
  piv = _ref!(M)[1]
  reduce_up!(M, piv)
  permute_rows!(M, piv)
  return length(piv) - length(findall(iszero, piv))
end

function rref(M::NfMatElem)
  N = deepcopy(M)
  return rref!(N), N
end

function permute_rows!(M::NfMatElem, pi::Vector{Int})
  j = 1
  for k=1:length(pi)
    i = pi[k]
    if i==0
      continue
    end
    if j != i
      @show :swap, i, j
      pi[i], pi[j] = pi[j], pi[i]
      swap_rows!(M, j, i)
    end
    j += 1
  end
  return nothing
end

function inv_new(M::NfMatElem)
  @assert nrows(M) == ncols(M)
  T = hcat(M, identity_matrix(M))
  rref!(T)
  return T[:, ncols(M)+1:end]
end

Base.similar(M::NfMatElem, r::Int, c::Int) = zero_matrix(base_ring(M), r, c)
Base.similar(M::NfMatElem) = zero_matrix(base_ring(M), nrows(M), ncols(M))

#not used.
function init!(A::Ptr{nf_elem_raw}, K::AnticNumberField)
  ccall((:nf_elem_init, Nemo.libantic), Cvoid, (Ptr{nf_elem_raw}, Ref{AnticNumberField}), A, K)
end

function Base.:*(x::NfMatElem, y::NfMatElem)
  ncols(x) != nrows(y) && error("Incompatible matrix dimensions")
  K = base_ring(x)
  @assert K == base_ring(y)

  A = similar(x, nrows(x), ncols(y))

  C = K()
  for i = 1:nrows(x)
    for j = 1:ncols(y)
      for k = 1:ncols(x)
        mul!(C, getindex_raw(x, i, k), getindex_raw(y, k, j), K, false)
        add!(getindex_raw(A, i, j), getindex_raw(A, i, j), C, K)
      end
      reduce!(getindex_raw(A, i, j), K)
    end
  end
  return A
end

function Hecke.zero!(a::Ptr{nf_elem_raw}, K::AnticNumberField)
  ccall((:nf_elem_zero, Nemo.libantic), Cvoid, (Ptr{nf_elem_raw}, Ref{AnticNumberField}), a, K)
end

#add version is used in strassen
#otherwise, classical matrix multiplication
function Nemo.mul!(C::NfMatElem, A::NfMatElem, B::NfMatElem, add::Bool = false)
  K = base_ring(A)
  @assert base_ring(B) == base_ring(C) == K
  T = K()
  S = K()
  for i=1:nrows(A)
    for j=1:ncols(B)
      Hecke.zero!(S)
      for k=1:ncols(A)
        mul!(T, getindex_raw(A, i, k), getindex_raw(B, k, j), K, false)
        Hecke.add!(S, S, T)
      end
      Hecke.reduce!(S)
      if add
        add!(getindex_raw(C, i, j), getindex_raw(C, i, j), S, K)
      else
        C[i,j] = S
      end
    end
  end
end

function Base.:-(A::NfMatElem, B::NfMatElem)
  C = zero_matrix(base_ring(A), nrows(A), ncols(A))
  return Nemo.sub!(C, A, B)
end

function Nemo.sub!(C::NfMatElem, A::NfMatElem, B::NfMatElem)
  K = base_ring(A)
  for i=1:nrows(A)
    for j=1:ncols(A)
      sub!(getindex_raw(C, i, j), getindex_raw(A, i, j), getindex_raw(B, i, j), K)
    end
  end
  return C
end

function Base.:+(A::NfMatElem, B::NfMatElem)
  C = zero_matrix(base_ring(A), nrows(A), ncols(A))
  return Nemo.add!(C, A, B)
end

function Nemo.add!(C::NfMatElem, A::NfMatElem, B::NfMatElem)
  K = base_ring(A)
  for i=1:nrows(A)
    for j=1:ncols(A)
      add!(getindex_raw(C, i, j), getindex_raw(A, i, j), getindex_raw(B, i, j), K)
    end
  end
  return C
end

#############################################################################
#
# matrix multiplication via fmpq_mat:
# write A as sum A_i t^t for A_i fmpq_mat's
# same for B and form the product, using A_i*B_j thus using fast(?)
# fmpq_mat_mul
#
# (ONLY available for degree(K) > 2, defaults to normal otherwise)
#
#############################################################################
function Hecke.coeff(M::NfMatElem, i::Int)
  return map(x->coeff(x, i), M)
end

function getindex_raw(M::fmpq_mat, i::Int, j::Int)
  return ccall((:fmpq_mat_entry, Nemo.libflint), Ptr{fmpq}, (Ref{fmpq_mat}, Int, Int), M, i-1, j-1)
end

function coeff!(m::fmpq_mat, M::NfMatElem, n::Int)
  K = base_ring(M)
  for i=1:nrows(M)
    for j=1:ncols(M)
      ccall((:nf_elem_get_coeff_fmpq, Nemo.libantic), Cvoid, (Ptr{fmpq}, Ptr{nf_elem_raw}, Int, Ref{AnticNumberField}), getindex_raw(m, i, j), getindex_raw(M, i, j), n, K)
    end
  end
end

function Hecke.denominator(a::Ptr{nf_elem_raw}, K::AnticNumberField)
  d = fmpz()
  denominator!(d, a, K)
  return d
end

function denominator!(d::fmpz, a::Ptr{nf_elem_raw}, K::AnticNumberField)
  ccall((:nf_elem_get_den, Nemo.libantic), Nothing, (Ref{fmpz}, Ptr{nf_elem_raw}, Ref{AnticNumberField}), d, a, K)
end

function setcoeff!(M::NfMatElem, n::Int, m::fmpq_mat)
  K = base_ring(M)
  @assert degree(K) > 2 #for now
  for i=1:nrows(M)
    for j=1:ncols(M)
      ccall((:fmpq_poly_set_coeff_fmpq, Nemo.libflint), Cvoid, (Ptr{nf_elem_raw}, Int, Ptr{fmpq}), getindex_raw(M, i, j), n, getindex_raw(m, i, j))
    end
  end
end

function reduce!(M::NfMatElem)
  K = base_ring(M)
  for i=1:nrows(M)
    for j=1:ncols(M)
      reduce!(getindex_raw(M, i, j), K)
    end
  end
end

function mul_new(A::NfMatElem, B::NfMatElem)
  K = base_ring(A)
  degree(K) > 2 || return A*B
  @assert degree(K) > 2
  C = similar(A, nrows(A), ncols(B))
  a = [coeff(A, i) for i=0:degree(K)-1]
  b = [coeff(B, i) for i=0:degree(K)-1]
  d = zero_matrix(QQ, nrows(C), ncols(C))
  for l=0:2*degree(K)-1
    c = zero_matrix(QQ, nrows(C), ncols(C))
    for i=max(0, l-degree(K)+1):min(degree(K)-1, l)
      #i+j = l, j = l-j
      Nemo.mul!(d, a[i+1], b[l-i+1])
      Nemo.add!(c, c, d)
    end
    setcoeff!(C, l, c)
  end
  reduce!(C)
  return C
end

#############################################################################
#
# multiplication via Kronnecker segmentation via fmpz_poly_mat
#
# seems to be fastest.
#
# Think of a NfMatElem A as a matrix of fmpq_poly
# Goal:
# C = A*B
# write A as dA * X for dA diagonal in Z and X integral (fmpz_poly)
# write B as Y * dB for bD diagonal in Z and Y integral
# C = dA * X * Y * dB
# and the middle product is fast via flint.
#
# (only for degree(K) > 2, defaults to normal otherwise)
#############################################################################
struct fmpz_poly_raw
  coeffs::Int #Ptr{fmpz}
  a::Int
  l::Int
  function fmpz_poly_raw()
    return new(0,0,0)
  end
  function fmpz_poly_raw(c::Int, a::Int, l::Int)
    return new(c, a, l)
  end
end

Base.zero(::Type{fmpz_poly_raw}) = fmpz_poly_raw()

struct fmpq_poly_raw
  coeffs::Int #Ptr{fmpz}
  den::Int #fmpz
  a::Int
  l::Int
  function fmpq_poly_raw()
    return new(0,1,0,0)
  end
  function fmpq_poly_raw(c::Int, d::Int, a::Int, l::Int)
    return new(c, d, a, l)
  end
end

Base.zero(::Type{fmpq_poly_raw}) = fmpq_poly_raw()

#minimal support to have printing (for debugging)
mutable struct fmpz_poly_mat <: MatElem{fmpz_poly}
   entries::Int #Ptr{fmpz_poly_raw}
   r::Int
   c::Int
   d::Int #Ptr{Ptr{fmpz_poly_raw}}
   function fmpz_poly_mat(r::Int, c::Int)
     mat = new()
     ccall((:fmpz_poly_mat_init, Nemo.libflint), Cvoid, (Ref{fmpz_poly_mat}, Int, Int), mat, r, c)
     finalizer(fmpz_poly_mat_clear, mat)
     return mat
   end
end

function fmpz_poly_mat_clear(M::fmpz_poly_mat)
  ccall((:fmpz_poly_mat_clear, Nemo.libflint), Cvoid, (Ref{fmpz_poly_mat}, ), M)
end

Hecke.base_ring(::fmpz_poly_mat) = Hecke.Globals.Zx
Hecke.nrows(M::fmpz_poly_mat) = M.r
Hecke.ncols(M::fmpz_poly_mat) = M.c

function Base.getindex(M::fmpz_poly_mat, i::Int, j::Int)
  f = Hecke.Globals.Zx()
  ccall((:fmpz_poly_set, Nemo.libflint), Cvoid, (Ref{fmpz_poly}, Ptr{fmpz_poly_raw}), f, getindex_raw(M, i, j))
  return f
end

function Base.checkbounds(M::fmpz_poly_mat, i::Int, j::Int)
  (1 <= i <= M.r && 1 <= j <= M.c) || throw(BoundsError(M, (i,j)))
end

@inline function getindex_raw(M::fmpz_poly_mat, i::Int, j::Int)
  @boundscheck checkbounds(M, i, j)
  return ccall((:fmpz_poly_mat_entry, Nemo.libflint), Ptr{fmpz_poly_raw}, (Ref{fmpz_poly_mat}, Int, Int), M, i-1, j-1)
end

@inline function Base.setindex!(M::fmpz_poly_mat, f::Ptr{nf_elem_raw}, i::Int, j::Int)
  @boundscheck checkbounds(M, i, j)
  ff = reinterpret(Ptr{Int}, f)
  B = [fmpz_poly_raw(unsafe_load(ff, 1), unsafe_load(ff, 3), unsafe_load(ff, 4))]
  Base.GC.@preserve B begin
    b = pointer(B, 1)
    ccall((:fmpz_poly_set, Nemo.libflint), Cvoid, (Ptr{fmpz_poly_raw}, Ptr{fmpz_poly_raw}), getindex_raw(M, i, j), b)
  end
end

@inline function Base.setindex!(M::NfMatElem, f::Ptr{fmpz_poly_raw}, i::Int, j::Int)
  @boundscheck checkbounds(M, i, j)
  ff = reinterpret(Ptr{Int}, f)
  B = [fmpq_poly_raw(unsafe_load(ff, 1), 1, unsafe_load(ff, 2), unsafe_load(ff, 3))]
  Base.GC.@preserve B begin
    b = pointer(B, 1)
    ccall((:fmpq_poly_set, Nemo.libflint), Cvoid, (Ptr{nf_elem_raw}, Ptr{fmpq_poly_raw}), getindex_raw(M, i, j), b)
  end
end

function mul_KS(A::NfMatElem, B::NfMatElem)
  K = base_ring(A)
  degree(K) > 2 || return A*B

  dA = [reduce(lcm, map(denominator, A[i, :]), init = fmpz(1)) for i=1:nrows(A)]
  tA = fmpz_poly_mat(nrows(A), ncols(A))
  tB = fmpz_poly_mat(nrows(B), ncols(B))
  tC = fmpz_poly_mat(nrows(A), ncols(B))
  Base.GC.@preserve tA tB tC begin
    for i=1:nrows(A)
      for j=1:ncols(A)
        tA[i,j] = getindex_raw(A, i, j)
        s = Nemo.divexact(dA[i], denominator(getindex_raw(A, i, j), K))
        if !isone(s)
          ccall((:fmpz_poly_scalar_mul_fmpz, Nemo.libflint), Cvoid, (Ptr{fmpz_poly_raw}, Ptr{fmpz_poly_raw}, Ref{fmpz}), getindex_raw(tA, i, j), getindex_raw(tA, i, j), s)
        end
      end
    end

    dB = [reduce(lcm, map(denominator, B[:, i]), init = fmpz(1)) for i=1:ncols(B)]
    for i=1:nrows(B)
      for j=1:ncols(B)
        tB[i,j] = getindex_raw(B, i, j)
        s = Nemo.divexact(dB[j], denominator(getindex_raw(B, i, j), K))
        if !isone(s)
          ccall((:fmpz_poly_scalar_mul_fmpz, Nemo.libflint), Cvoid, (Ptr{fmpz_poly_raw}, Ptr{fmpz_poly_raw}, Ref{fmpz}), getindex_raw(tB, i, j), getindex_raw(tB, i, j), s)
        end
      end
    end

    ccall((:fmpz_poly_mat_mul, Nemo.libflint), Cvoid, (Ref{fmpz_poly_mat}, Ref{fmpz_poly_mat}, Ref{fmpz_poly_mat}), tC, tA, tB)

    C = zero_matrix(K, nrows(A), ncols(B))
    for i=1:nrows(C)
      for j=1:ncols(C)
        C[i,j] = getindex_raw(tC, i, j)
        ccall((:nf_elem_set_den, Nemo.libantic), Cvoid, (Ptr{nf_elem_raw}, Ref{fmpz}, Ref{AnticNumberField}), getindex_raw(C, i, j), dA[i]*dB[j], K)
        reduce!(getindex_raw(C, i, j), K)
      end
    end
  end
  return C
end

export NfMatElem

end # NfMatModule

module TestNfMatMul
using Hecke
using Main.NfMatModule

"""
    test_mul(k::AnticNumberField, a::Int, r::AbstractVector{fmpz})

Creates two random `a x a` matrix with coefficients in `k` of size in `r`
and multiplies them in various ways
"""
function test_mul(k::AnticNumberField, a::Int, r::AbstractVector{fmpz})
  return test_mul(k, a, a, a, r)
end

"""
    test_mul(k::AnticNumberField, a::Int, b::Int, c::Int, r::AbstractVector{fmpz})

Create a `a x b` matrix and a `b x c` matrix with coefficients in `k` of size
`r` and mupliplies them in various ways.
"""
function test_mul(k::AnticNumberField, a::Int, b::Int, c::Int, r::AbstractVector{fmpz})
  M = matrix(k, a, b, [rand(k, r) for i=1:a*b])
  N = matrix(k, b, c, [rand(k, r) for i=1:b*c])
  C = zero_matrix(k, a, c)
  m = Generic.Mat{nf_elem}(M)
  n = Generic.Mat{nf_elem}(N)
  Base.GC.gc()
  @time m*n
  Base.GC.gc()
  @time M*N
  if isdefined(Main, :Strassen)
    Base.GC.gc()
    @time Main.Strassen.mul_strassen!(C, M, N)
  end
  Base.GC.gc()
  @time NfMatModule.mul_new(M, N)
  Base.GC.gc()
  @time NfMatModule.mul_KS(M, N)
end

end

module SpinModule
using Hecke
using Main.NfMatModule
import Main.NfMatModule.NfMatElem

function Hecke.mul!(a::NfMatElem, M::SMat{nf_elem}, b::NfMatElem)
  #a = M*b
  K = base_ring(a)
  s = K()
  for i=1:ncols(b)
    t = Main.NfMatModule.getindex_raw(a, 1, i)
    zero!(t, K)
    for (j, v) = M[i]
      Main.NfMatModule.mul!(s, v, Main.NfMatModule.getindex_raw(b, 1, j), K, false)
      Main.NfMatModule.add!(t, t, s, K)
    end
    Main.NfMatModule.reduce!(t, K)
    a[1, i] = t
  end
end

#TODO: use inplace here and add trafo....
# - think about _ref:
#   - options for full?
#   - trafo only for start..stop? as an option?
function spin(A::SMat{nf_elem}, b::NfMatElem)
  b = deepcopy(b)
  c = similar(b)
  s = deepcopy(b)
  k = base_ring(b)
  trafo = identity_matrix(k, 1)
  piv, _ = Main.NfMatModule._ref!(s, trafo = trafo)
  while true
    Hecke.mul!(c, A, b)
    vcat!(s, c)
    if nrows(trafo) < nrows(s)
      tr = identity_matrix(k, 2*nrows(trafo))
      for i=1:nrows(trafo)
        for j=1:nrows(trafo)
          tr[i,j] = Main.NfMatModule.getindex_raw(trafo, i, j)
        end
      end
      trafo = tr
    end
    Main.NfMatModule._ref!(s, start = nrows(s), piv = piv, trafo = trafo)
    if is_zero_row(s, nrows(s))
      return s, piv, trafo#[end, :]
    end
    c, b = b, c
  end
end

function spin(A::SMat{nf_elem}, b::NfMatElem, quo::NfMatElem, qpiv ::Vector{Int})
  b = deepcopy(b)
  k = base_ring(b)
  Main.NfMatModule.reduce!(b, quo, qpiv)
  if is_zero_row(b, 1)
    return b, zeros(Int, ncols(b)), identity_matrix(k, 1)
  end
  c = similar(b)
  s = deepcopy(b)
  trafo = identity_matrix(k, 1)
  piv, _ = Main.NfMatModule._ref!(s, trafo = trafo)
  while true
    Hecke.mul!(c, A, b)
    Main.NfMatModule.reduce!(c, quo, qpiv)
    vcat!(s, c)
    if nrows(trafo) < nrows(s)
      tr = identity_matrix(k, 2*nrows(trafo))
      for i=1:nrows(trafo)
        for j=1:nrows(trafo)
          tr[i,j] = Main.NfMatModule.getindex_raw(trafo, i, j)
        end
      end
      trafo = tr
    end

    Main.NfMatModule._ref!(s, start = nrows(s), piv = piv, trafo = trafo)
    if is_zero_row(s, nrows(s))
      return s, piv, trafo#[end, :]
    end
    c, b = b, c
  end
end

function charpoly_fac_elem(A::SMat{nf_elem})
  k = base_ring(A)

  kx, _ = PolynomialRing(k)
  lf = elem_type(kx)[]

  q = zero_matrix(k, 0, ncols(A))
  pq = zeros(Int, ncols(A))

  e = zero_matrix(k, 1, ncols(A))

  while true
    i = 1
    while i <= ncols(A) && pq[i] != 0
      i += 1
    end
    if i>ncols(A)
      return lf
    end
    for j=1:ncols(A)
      zero!(Main.NfMatModule.getindex_raw(e, 1, j), k)
    end
    e[1,i] = 1

    s, p, t = spin(A, e, q, pq)
    @assert is_zero_row(s, nrows(s))
    @assert !is_zero_row(s, nrows(s)-1)

    push!(lf, kx(vec(collect(view(t, nrows(s):nrows(s), 1:nrows(s))))))
    lb = nrows(q)
#    Main.NfMatModule.reduce_up!(s, p)
#    Main.NfMatModule.reduce!(q, s, p)
    vcat!(q, view(s, 1:nrows(s)-1, :))
    @show size(q)
    for i=1:length(p)
      if p[i] != 0
        @assert pq[i] == 0
        pq[i] = p[i]+lb
      end
    end
    @show MSet(pq)
  end
end

end #SpinModule


text:
# This is to facilitate conversion of MATLAB codes to Julia, if needed.

@doc """
meshgrid is similar to MATLAB's function. It takes in a range, vector, matrix or similar
and makes an array of those values.
In the case of a cartesian grid, it would like an 'X' and 'Y' vectors and make a matrix filled with
the X and Y coordinates:
julia> X = 0:3
julia> Y = 0:3
julia> XX,YY = meshgrid(X,Y)
julia> XX
4×4 Array{Int64,2}:
 0  0  0  0
 1  1  1  1
 2  2  2  2
 3  3  3  3
 julia> YY
 4×4 Array{Int64,2}:
 0  1  2  3
 0  1  2  3
 0  1  2  3
 0  1  2  3

 If a single vector/range is inputted, for example X = 0:3 , it will output interpret it as meshgrid(X,X)

"""-> function meshgrid(xs, ys)
    [xs[i] for i in 1:length(xs), j in 1:length(ys)], [ys[j] for i in 1:length(xs), j in 1:length(ys)]
end

function meshgrid(xs;dims=3)
    if dims==3
        meshgrid(xs, xs,xs)
    elseif dims==2
        meshgrid(xs, xs)
    end

end

function meshgrid(xs,ys,zs)
    return ([xs[i] for i in 1:length(xs), j in 1:length(ys), k in 1:length(zs)],
    [ys[j] for i in 1:length(xs), j in 1:length(ys), k in 1:length(zs)],
    [zs[k] for i in 1:length(xs), j in 1:length(ys), k in 1:length(zs)])

end


text:
### Closed-Form UB for the MAD Case

#Assumes distribution already standardized:  mu = 1, M = 1
delta_l(S) = -log(S) / lambertw(-exp(-1) / S)
delta_m(S) = log(S) / (1 + log(S))
delta_h(S) = (S - 1) / S

#Assumes Standardized:  mu = 1, S = 1
function _vopp_ub_MAD_cf(S, D)
	if 0 <= D <= delta_l(S)
		return -lambertw((D - 1) / exp(1)) / (1 - D)
	elseif D <= delta_m(S)
		return log(S) / D
	elseif D <= delta_h(S)
		return -lambertw( -exp(-1) / S / (1 - D) )
	else
		throw("D $D < $(delta_h(S)) is too large for Scale $S.  Infeasible")
	end
	return -1.
end

function vopp_ub_MAD(S, M, D; method = :ClosedForm)
	if method == :ClosedForm
		Sc = vopp.comp_Sc(S, M)
		return _vopp_ub_MAD_cf(Sc, D / M)
	elseif method == :Optimization
		return vopp_ub_MAD_opt(S, M, D)[1]
	else
		throw("Method must be one of :ClosedForm or :Optimization : $method")
	end
	return -1.
end

#assumes things already standardized
function tight_dist_ub_MAD(x, S, D)
	alpha = 1/vopp_ub_MAD(S, 1., D)
	#match which regime
	if 0 <= D <= delta_l(S)  #low heterogeneity
		if x <= alpha
			return 1.
		elseif x <= 1
			return alpha / x
		elseif x <= S
			return D / log(S) / x
		else
			return 0.
		end
	elseif D <= delta_m(S)  #medium
		if x == 0
			return 1
		elseif x <= exp(1) * S^(1 - 1/D)
			return alpha / exp(1) * S^(1/D - 1)
		elseif x <= S
			return alpha / x
		else
			return 0.
		end
	elseif D <= delta_h(S) #high
		if x == 0
			return 1.
		elseif x <= alpha/(1-D)
			return 1 - D
		elseif x <= S
			return alpha / x
		else
			return 0.
		end
	else
		throw("D not valid; $D not in [0, $(delta_h(S))")
	end
end


text:
using JET, InteractiveUtils

const CC = JET.CC

import .CC:
    widenconst, ⊑, Bottom

import JET:
    AbstractAnalyzer,
    JETAnalyzer,
    AbstractGlobal,
    get_result,
    ToplevelConfig,
    virtual_process,
    virtualize_module_context,
    gen_virtual_module,
    ToplevelErrorReport,
    InferenceErrorReport,
    print_reports,
    get_reports,
    get_cache

function subtypes_recursive!(t, ts)
    push!(ts, t)
    if isabstracttype(t)
        for ct in subtypes(t)
            subtypes_recursive!(ct, ts)
        end
    end
    return ts
end

let
    ts = Type[]
    subtypes_recursive!(ToplevelErrorReport, ts)
    subtypes_recursive!(InferenceErrorReport, ts)
    for t in ts
        canonicalname = Symbol(parentmodule(t), '.', nameof(t))
        canonicalpath = Symbol.(split(string(canonicalname), '.'))

        modpath = Expr(:., canonicalpath[1:end-1]...)
        symname = Expr(:., last(canonicalpath))
        ex = Expr(:import, Expr(:(:), modpath, symname))
        Core.eval(@__MODULE__, ex)
    end
end

import Base.Meta:
    isexpr

"""
    @fixturedef

Creates a virtual module and defines fixtures from toplevel expression `ex`
"""
macro fixturedef(ex)
    @assert isexpr(ex, :block)
    return quote let
        vmod = $gen_virtual_module()
        for x in $(ex.args)
            Core.eval(vmod, x)
        end
        vmod # return virtual module
    end end
end

"""
    @analyze_toplevel [jetconfigs...] ex -> JET.VirtualProcessResult

Enters JET analysis from toplevel expression `ex`, and returns the analysis result.
"""
macro analyze_toplevel(xs...)
    jetconfigs = filter(iskwarg, xs)
    xs′ = filter(!iskwarg, xs)
    @assert length(xs′) == 1
    ex = first(xs′)
    return _analyze_toplevel(ex, __source__, jetconfigs)
end

"""
    @analyze_toplevel2 [jetconfigs...] ex -> (Module, JET.VirtualProcessResult)

Works similarly to `@analyze_toplevel`, but also creates a virtual module beforehand and
returns that as well after the whole analysis for the later inspection.
"""
macro analyze_toplevel2(xs...)
    jetconfigs = filter(iskwarg, xs)
    xs′ = filter(!iskwarg, xs)
    @assert length(xs′) == 1
    ex = first(xs′)

    vmod = gensym(:vmod)
    jetconfigs = (:(context = $vmod), jetconfigs...,)
    jetconfigs = (:(virtualize = false), jetconfigs...,)
    ex2 = _analyze_toplevel(ex, __source__, jetconfigs)
    return :(let
        $(esc(vmod)) = $gen_virtual_module()
        ret2 = $ex2
        $(esc(vmod)), ret2
    end)
end

iskwarg(@nospecialize(x)) = isexpr(x, :(=))

function _analyze_toplevel(ex, lnn, jetconfigs)
    toplevelex = (isexpr(ex, :block) ?
                  Expr(:toplevel, lnn, ex.args...) : # flatten here
                  Expr(:toplevel, lnn, ex)
                  ) |> QuoteNode
    return :(let
        analyzer = $(GlobalRef(JET, :JETAnalyzer))(; $(map(esc, jetconfigs)...))
        config = ToplevelConfig(; $(map(esc, jetconfigs)...))
        $virtual_process($toplevelex,
                         $(string(lnn.file)),
                         analyzer,
                         config,
                         )
    end)
end

# `report_file` with silent top-level logger
report_file2(args...; kwargs...) =
    report_file(args...; toplevel_logger = nothing, kwargs...)

is_concrete(mod, sym) = isdefined(mod, sym) && !isa(getfield(mod, sym), AbstractGlobal)
is_abstract(mod, sym) = isdefined(mod, sym) && isa(getfield(mod, sym), AbstractGlobal)
isa_abstract(x, @nospecialize(typ)) = isa(x, AbstractGlobal) && x.t ⊑ typ

# JET will try to concretize global variable when its type is a constant at analysis time,
# but the starategy is a bit complicated right now and may change in the future
# these utilities allow robust testing to check if a object is successfully analyzed by JET whichever it's concretized or abstracted
is_analyzed(mod, sym) = isdefined(mod, sym) # essentially, `is_concrete(mod, sym) || is_abstract(mod, sym)`
isa_analyzed(x, @nospecialize(typ)) = isa_abstract(x, typ) || isa(x, typ)

# fresh execution/benchmark tools
include(normpath(@__DIR__, "..", "benchmark", "JETBenchmarkUtils.jl"))
using .JETBenchmarkUtils


text:
# Set up the environment to run this example. Make sure you're within the folder that this
# file lives in.
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using Turing, Stheno, RDatasets, Random, Plots, Statistics
rng = MersenneTwister(5)
T = Float64

# Use RDatasets to load a binary classification dataset, we use features from this dataset to predict it's color
# The crabs data frame has 200 rows and 8 columns, describing 5 morphological measurements on 50 crabs
# each of two colour forms and both sexes, of the species Leptograpsus variegatus collected at Fremantle, W. Australia.
crabs = dataset("MASS","crabs")
crabs = crabs[shuffle(1:size(crabs, 1)), :]
train = crabs[1:div(end,2), :]
test = crabs[div(end,2)+1:end, :]

train_y = Array{Bool}(undef,size(train, 1))
train_y[train.Sp.=="B"].=0
train_y[train.Sp.=="O"].=1
train_x = Matrix(transpose(convert(Array,train[:,4:end])))

test_y = Array{Bool}(undef, size(test, 1))
test_y[test.Sp.=="B"].=0
test_y[test.Sp.=="O"].=1
test_x = Matrix(transpose(convert(Array, test[:, 4:end])))

# The probabilistic model for the binary classification problem can be formulated as:
# f ~ GP(μ, K)
# y|f ~ Bernoulli(sigmoid(f))
# This model is build using Turing.jl, please refer to https://turing.ml/dev/docs/using-turing/get-started for details.

σ(x) = one(T) / (one(T)+exp(-x))

function build_gp(logl, σ², X)
    ard_eq_kernel = σ² * stretch(SEKernel(), exp.(-logl))
    gp = GP(ard_eq_kernel, GPC())
    prior = gp(ColVecs(X), T(0.01))
    gp, prior
end

# The Turing model used to estimate the posterior distribution,
# the latent variable is f & the parameter is logl
@model gpc_learn(X, y) = begin
    logl ~ Normal(T(0.0), T(2.0))
    _, prior = build_gp(logl, T(1.0), X)
    f ~ prior
    for i in eachindex(y)
      	y[i] ~ Bernoulli(σ(f[i]))
    end
end

# Function used to infer the label for newly inputs
# NOTE: for simplicity, we use MAP estimation here instead of integrate over the latent and parameters
function gpc_infer(x, logl, Xtrain, fsamples)
    nsamples = size(fsamples, 2)
    fxs = []
    for i in 1:nsamples
        gp, prior = build_gp(logl[i], T(1.0), Xtrain)
        conditioned_gp = gp | Obs(prior, fsamples[:, i])
        posterior = conditioned_gp(ColVecs(x))
        push!(fxs, mean.(marginals(posterior)))
    end
    fx_mean = vec(mean(hcat(fxs...), dims=2))
    p = σ.(fx_mean)
    y = Int.(p .> T(0.5))
    y
end

# Marginalize over non Gaussian likelihood is intractable,
# we use MCMC to compute the approximate posterior
model = gpc_learn(train_x, train_y)
mcmc_samples = sample(model, HMC(0.01, 10), 5000);

# store the sampled parameter & latent variable, descard the first 1000 samples
logl_df = mcmc_samples[:logl]
logl = vec(logl_df.value.data)
logl = Array{T}(logl)
reserve_logl = logl[1001:end]

fsamples_df = mcmc_samples[:f]
fsamples = Matrix(transpose(dropdims(fsamples_df.value.data, dims=3)))
fsamples = convert.(T, fsamples)
reserve_fsamples = fsamples[:, 1001:end]

# prediction
pred_y = gpc_infer(test_x, reserve_logl, train_x, reserve_fsamples)
# determine the accuracy of our prediction
function accuracy(pred_y, y)
    N = length(y)
    N_neq = sum(abs.(pred_y .- y))
    r = T(1.0) - N_neq / N
    r
end

@show accuracy(pred_y, test_y)


text:
function DSSTPropagator(arg0::ODEIntegrator)
    return DSSTPropagator((ODEIntegrator,), arg0)
end

function DSSTPropagator(arg0::ODEIntegrator, arg1::PropagationType)
    return DSSTPropagator((ODEIntegrator, PropagationType), arg0, arg1)
end

function DSSTPropagator(arg0::ODEIntegrator, arg1::PropagationType, arg2::AttitudeProvider)
    return DSSTPropagator((ODEIntegrator, PropagationType, AttitudeProvider), arg0, arg1, arg2)
end

function add_additional_equations(obj::AbstractIntegratedPropagator, arg0::AdditionalEquations)
    return jcall(obj, "addAdditionalEquations", void, (AdditionalEquations,), arg0)
end

function add_additional_state_provider(obj::AbstractPropagator, arg0::AdditionalStateProvider)
    return jcall(obj, "addAdditionalStateProvider", void, (AdditionalStateProvider,), arg0)
end

function add_event_detector(obj::AbstractIntegratedPropagator, arg0::EventDetector)
    return jcall(obj, "addEventDetector", void, (EventDetector,), arg0)
end

function add_force_model(obj::DSSTPropagator, arg0::DSSTForceModel)
    return jcall(obj, "addForceModel", void, (DSSTForceModel,), arg0)
end

function clear_events_detectors(obj::AbstractIntegratedPropagator)
    return jcall(obj, "clearEventsDetectors", void, ())
end

function compute_mean_state(::Type{DSSTPropagator}, arg0::SpacecraftState, arg1::AttitudeProvider, arg2::Collection)
    return jcall(DSSTPropagator, "computeMeanState", SpacecraftState, (SpacecraftState, AttitudeProvider, Collection), arg0, arg1, arg2)
end

function compute_mean_state(::Type{DSSTPropagator}, arg0::SpacecraftState, arg1::AttitudeProvider, arg2::Collection, arg3::jdouble, arg4::jint)
    return jcall(DSSTPropagator, "computeMeanState", SpacecraftState, (SpacecraftState, AttitudeProvider, Collection, jdouble, jint), arg0, arg1, arg2, arg3, arg4)
end

function compute_osculating_state(::Type{DSSTPropagator}, arg0::SpacecraftState, arg1::AttitudeProvider, arg2::Collection)
    return jcall(DSSTPropagator, "computeOsculatingState", SpacecraftState, (SpacecraftState, AttitudeProvider, Collection), arg0, arg1, arg2)
end

function get_additional_state_providers(obj::AbstractPropagator)
    return jcall(obj, "getAdditionalStateProviders", List, ())
end

function get_all_force_models(obj::DSSTPropagator)
    return jcall(obj, "getAllForceModels", List, ())
end

function get_attitude_provider(obj::AbstractPropagator)
    return jcall(obj, "getAttitudeProvider", AttitudeProvider, ())
end

function get_basic_dimension(obj::AbstractIntegratedPropagator)
    return jcall(obj, "getBasicDimension", jint, ())
end

function get_calls(obj::AbstractIntegratedPropagator)
    return jcall(obj, "getCalls", jint, ())
end

function get_events_detectors(obj::AbstractIntegratedPropagator)
    return jcall(obj, "getEventsDetectors", Collection, ())
end

function get_frame(obj::AbstractPropagator)
    return jcall(obj, "getFrame", Frame, ())
end

function get_generated_ephemeris(obj::AbstractIntegratedPropagator)
    return jcall(obj, "getGeneratedEphemeris", BoundedPropagator, ())
end

function get_initial_state(obj::AbstractPropagator)
    return jcall(obj, "getInitialState", SpacecraftState, ())
end

function get_managed_additional_states(obj::AbstractIntegratedPropagator)
    return jcall(obj, "getManagedAdditionalStates", Vector{JString}, ())
end

function get_mode(obj::AbstractPropagator)
    return jcall(obj, "getMode", jint, ())
end

function get_mu(obj::AbstractIntegratedPropagator)
    return jcall(obj, "getMu", jdouble, ())
end

function get_orbit_type(obj::DSSTPropagator)
    return jcall(obj, "getOrbitType", OrbitType, ())
end

function get_position_angle_type(obj::DSSTPropagator)
    return jcall(obj, "getPositionAngleType", PositionAngle, ())
end

function get_pv_coordinates(obj::AbstractPropagator, arg0::AbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedPVCoordinates, (AbsoluteDate, Frame), arg0, arg1)
end

function get_satellite_revolution(obj::DSSTPropagator)
    return jcall(obj, "getSatelliteRevolution", jint, ())
end

function get_selected_coefficients(obj::DSSTPropagator)
    return jcall(obj, "getSelectedCoefficients", Set, ())
end

function initial_is_osculating(obj::DSSTPropagator)
    return jcall(obj, "initialIsOsculating", jboolean, ())
end

function is_additional_state_managed(obj::AbstractIntegratedPropagator, arg0::JString)
    return jcall(obj, "isAdditionalStateManaged", jboolean, (JString,), arg0)
end

function propagate(obj::AbstractIntegratedPropagator, arg0::AbsoluteDate)
    return jcall(obj, "propagate", SpacecraftState, (AbsoluteDate,), arg0)
end

function propagate(obj::AbstractIntegratedPropagator, arg0::AbsoluteDate, arg1::AbsoluteDate)
    return jcall(obj, "propagate", SpacecraftState, (AbsoluteDate, AbsoluteDate), arg0, arg1)
end

function remove_force_models(obj::DSSTPropagator)
    return jcall(obj, "removeForceModels", void, ())
end

function reset_initial_state(obj::DSSTPropagator, arg0::SpacecraftState)
    return jcall(obj, "resetInitialState", void, (SpacecraftState,), arg0)
end

function set_attitude_provider(obj::DSSTPropagator, arg0::AttitudeProvider)
    return jcall(obj, "setAttitudeProvider", void, (AttitudeProvider,), arg0)
end

function set_ephemeris_mode(obj::AbstractIntegratedPropagator)
    return jcall(obj, "setEphemerisMode", void, ())
end

function set_ephemeris_mode(obj::AbstractIntegratedPropagator, arg0::OrekitStepHandler)
    return jcall(obj, "setEphemerisMode", void, (OrekitStepHandler,), arg0)
end

function set_initial_state(obj::DSSTPropagator, arg0::SpacecraftState)
    return jcall(obj, "setInitialState", void, (SpacecraftState,), arg0)
end

function set_initial_state(obj::DSSTPropagator, arg0::SpacecraftState, arg1::PropagationType)
    return jcall(obj, "setInitialState", void, (SpacecraftState, PropagationType), arg0, arg1)
end

function set_interpolation_grid_to_fixed_number_of_points(obj::DSSTPropagator, arg0::jint)
    return jcall(obj, "setInterpolationGridToFixedNumberOfPoints", void, (jint,), arg0)
end

function set_interpolation_grid_to_max_time_gap(obj::DSSTPropagator, arg0::jdouble)
    return jcall(obj, "setInterpolationGridToMaxTimeGap", void, (jdouble,), arg0)
end

function set_master_mode(obj::AbstractIntegratedPropagator, arg0::OrekitStepHandler)
    return jcall(obj, "setMasterMode", void, (OrekitStepHandler,), arg0)
end

function set_master_mode(obj::AbstractPropagator, arg0::jdouble, arg1::OrekitFixedStepHandler)
    return jcall(obj, "setMasterMode", void, (jdouble, OrekitFixedStepHandler), arg0, arg1)
end

function set_mu(obj::DSSTPropagator, arg0::jdouble)
    return jcall(obj, "setMu", void, (jdouble,), arg0)
end

function set_reset_at_end(obj::AbstractIntegratedPropagator, arg0::jboolean)
    return jcall(obj, "setResetAtEnd", void, (jboolean,), arg0)
end

function set_satellite_revolution(obj::DSSTPropagator, arg0::jint)
    return jcall(obj, "setSatelliteRevolution", void, (jint,), arg0)
end

function set_selected_coefficients(obj::DSSTPropagator, arg0::Set)
    return jcall(obj, "setSelectedCoefficients", void, (Set,), arg0)
end

function set_slave_mode(obj::AbstractIntegratedPropagator)
    return jcall(obj, "setSlaveMode", void, ())
end

function tolerances(::Type{DSSTPropagator}, arg0::jdouble, arg1::Orbit)
    return jcall(DSSTPropagator, "tolerances", Vector{Vector{jdouble}}, (jdouble, Orbit), arg0, arg1)
end


text:
@testset "API functions" begin
    test_ex(APIFunction(func_by_name(:vkEnumeratePhysicalDevices), false), :(
        function _enumerate_physical_devices(instance)::ResultTypes.Result{Vector{PhysicalDevice},VulkanError}
            pPhysicalDeviceCount = Ref{UInt32}()
            @repeat_while_incomplete begin
                @check @dispatch instance vkEnumeratePhysicalDevices(instance, pPhysicalDeviceCount, C_NULL)
                pPhysicalDevices = Vector{VkPhysicalDevice}(undef, pPhysicalDeviceCount[])
                @check @dispatch instance vkEnumeratePhysicalDevices(instance, pPhysicalDeviceCount, pPhysicalDevices)
            end
            PhysicalDevice.(pPhysicalDevices, identity, instance)
        end
    ))

    test_ex(APIFunction(func_by_name(:vkGetPhysicalDeviceProperties), false), :(
        function _get_physical_device_properties(physical_device)::_PhysicalDeviceProperties
            pProperties = Ref{VkPhysicalDeviceProperties}()
            @dispatch instance(physical_device) vkGetPhysicalDeviceProperties(physical_device, pProperties)
            from_vk(_PhysicalDeviceProperties, pProperties[])
        end
    ))

    test_ex(APIFunction(func_by_name(:vkEnumerateInstanceVersion), false), :(
        function _enumerate_instance_version()::ResultTypes.Result{VersionNumber,VulkanError}
            pApiVersion = Ref{UInt32}()
            @check @dispatch nothing vkEnumerateInstanceVersion(pApiVersion)
            from_vk(VersionNumber, pApiVersion[])
        end
    ))

    test_ex(APIFunction(func_by_name(:vkEnumerateInstanceExtensionProperties), false), :(
        function _enumerate_instance_extension_properties(; layer_name = C_NULL)::ResultTypes.Result{Vector{_ExtensionProperties},VulkanError}
            pPropertyCount = Ref{UInt32}()
            @repeat_while_incomplete begin
                @check @dispatch nothing vkEnumerateInstanceExtensionProperties(layer_name, pPropertyCount, C_NULL)
                pProperties = Vector{VkExtensionProperties}(undef, pPropertyCount[])
                @check @dispatch nothing vkEnumerateInstanceExtensionProperties(layer_name, pPropertyCount, pProperties)
            end
            from_vk.(_ExtensionProperties, pProperties)
        end
    ))

    test_ex(APIFunction(func_by_name(:vkGetGeneratedCommandsMemoryRequirementsNV), false), :(
        function _get_generated_commands_memory_requirements_nv(device, info::_GeneratedCommandsMemoryRequirementsInfoNV, next_types::Type...)::_MemoryRequirements2
            memory_requirements = initialize(_MemoryRequirements2, next_types...)
            pMemoryRequirements = Ref(Base.unsafe_convert(VkMemoryRequirements2, memory_requirements))
            GC.@preserve memory_requirements begin
                @dispatch device vkGetGeneratedCommandsMemoryRequirementsNV(device, info, pMemoryRequirements)
                _MemoryRequirements2(pMemoryRequirements[], Any[memory_requirements])
            end
        end
    ))

    test_ex(APIFunction(func_by_name(:vkGetInstanceProcAddr), false), :(_get_instance_proc_addr(name::AbstractString; instance = C_NULL)::FunctionPtr = vkGetInstanceProcAddr(instance, name)))
    test_ex(APIFunction(func_by_name(:vkGetInstanceProcAddr), true), :(_get_instance_proc_addr(name::AbstractString, fptr::FunctionPtr; instance = C_NULL)::FunctionPtr = vkGetInstanceProcAddr(instance, name, fptr)))

    test_ex(APIFunction(func_by_name(:vkGetPhysicalDeviceSurfacePresentModesKHR), false), :(
        function _get_physical_device_surface_present_modes_khr(physical_device; surface = C_NULL)::ResultTypes.Result{Vector{PresentModeKHR},VulkanError}
            pPresentModeCount = Ref{UInt32}()
            @repeat_while_incomplete begin
                @check @dispatch instance(physical_device) vkGetPhysicalDeviceSurfacePresentModesKHR(physical_device, surface, pPresentModeCount, C_NULL)
                pPresentModes = Vector{VkPresentModeKHR}(undef, pPresentModeCount[])
                @check @dispatch instance(physical_device) vkGetPhysicalDeviceSurfacePresentModesKHR(physical_device, surface, pPresentModeCount, pPresentModes)
            end
            pPresentModes
        end
    ))

    test_ex(APIFunction(func_by_name(:vkGetRandROutputDisplayEXT), false), :(
        function _get_rand_r_output_display_ext(physical_device, dpy::Ptr{vk.Display}, rr_output::vk.RROutput)::ResultTypes.Result{DisplayKHR,VulkanError}
            pDisplay = Ref{VkDisplayKHR}()
            @check @dispatch instance(physical_device) vkGetRandROutputDisplayEXT(physical_device, dpy, rr_output, pDisplay)
            DisplayKHR(pDisplay[], identity, physical_device)
        end
    ))

    test_ex(APIFunction(func_by_name(:vkRegisterDeviceEventEXT), false), :(
        function _register_device_event_ext(device, device_event_info::_DeviceEventInfoEXT; allocator = C_NULL)::ResultTypes.Result{Fence,VulkanError}
            pFence = Ref{VkFence}()
            @check @dispatch device vkRegisterDeviceEventEXT(device, device_event_info, allocator, pFence)
            Fence(pFence[], (x -> _destroy_fence(device, x; allocator)), device)
        end
    ))

    test_ex(APIFunction(func_by_name(:vkCreateInstance), false), :(
        function _create_instance(create_info::_InstanceCreateInfo; allocator = C_NULL)::ResultTypes.Result{Instance,VulkanError}
            pInstance = Ref{VkInstance}()
            @check @dispatch nothing vkCreateInstance(create_info, allocator, pInstance)
            @fill_dispatch_table Instance(pInstance[], x -> _destroy_instance(x; allocator))
        end
    ))

    test_ex(APIFunction(func_by_name(:vkCreateDebugReportCallbackEXT), true), :(
        function _create_debug_report_callback_ext(instance, create_info::_DebugReportCallbackCreateInfoEXT, fptr_create::FunctionPtr, fptr_destroy::FunctionPtr; allocator = C_NULL)::ResultTypes.Result{DebugReportCallbackEXT,VulkanError}
            pCallback = Ref{VkDebugReportCallbackEXT}()
            @check vkCreateDebugReportCallbackEXT(instance, create_info, allocator, pCallback, fptr_create)
            DebugReportCallbackEXT(pCallback[], (x -> _destroy_debug_report_callback_ext(instance, x, fptr_destroy; allocator)), instance)
        end
    ))

    test_ex(APIFunction(func_by_name(:vkCreateGraphicsPipelines), false), :(
        function _create_graphics_pipelines(device, create_infos::AbstractArray; pipeline_cache = C_NULL, allocator = C_NULL)::ResultTypes.Result{Tuple{Vector{Pipeline},Result},VulkanError}
            pPipelines = Vector{VkPipeline}(undef, pointer_length(create_infos))
            @check @dispatch device vkCreateGraphicsPipelines(device, pipeline_cache, pointer_length(create_infos), create_infos, allocator, pPipelines)
            (Pipeline.(pPipelines, x -> _destroy_pipeline(device, x; allocator), device), _return_code)
        end
    ))

    test_ex(APIFunction(func_by_name(:vkAllocateDescriptorSets), false), :(
        function _allocate_descriptor_sets(device, allocate_info::_DescriptorSetAllocateInfo)::ResultTypes.Result{Vector{DescriptorSet},VulkanError}
            pDescriptorSets = Vector{VkDescriptorSet}(undef, allocate_info.vks.descriptorSetCount)
            @check @dispatch device vkAllocateDescriptorSets(device, allocate_info, pDescriptorSets)
            DescriptorSet.(pDescriptorSets, identity, getproperty(allocate_info, :descriptor_pool))
        end
    ))

    test_ex(APIFunction(func_by_name(:vkMergePipelineCaches), false), :(
        _merge_pipeline_caches(device, dst_cache, src_caches::AbstractArray)::ResultTypes.Result{Result,VulkanError} = @check(@dispatch(device, vkMergePipelineCaches(device, dst_cache, pointer_length(src_caches), src_caches)))
    ))

    test_ex(APIFunction(func_by_name(:vkGetFenceFdKHR), false), :(
        function _get_fence_fd_khr(device, get_fd_info::_FenceGetFdInfoKHR)::ResultTypes.Result{Int,VulkanError}
            pFd = Ref{Int}()
            @check @dispatch device vkGetFenceFdKHR(device, get_fd_info, pFd)
            pFd[]
        end
    ))

    test_ex(APIFunction(func_by_name(:vkGetDeviceGroupPeerMemoryFeatures), false), :(
        function _get_device_group_peer_memory_features(device, heap_index::Integer, local_device_index::Integer, remote_device_index::Integer)::PeerMemoryFeatureFlag
            pPeerMemoryFeatures = Ref{VkPeerMemoryFeatureFlags}()
            @dispatch device vkGetDeviceGroupPeerMemoryFeatures(device, heap_index, local_device_index, remote_device_index, pPeerMemoryFeatures)
            pPeerMemoryFeatures[]
        end
    ))

    test_ex(APIFunction(func_by_name(:vkUpdateDescriptorSets), false), :(
        _update_descriptor_sets(device, descriptor_writes::AbstractArray, descriptor_copies::AbstractArray)::Cvoid = @dispatch device vkUpdateDescriptorSets(device, pointer_length(descriptor_writes), descriptor_writes, pointer_length(descriptor_copies), descriptor_copies)
    ))

    test_ex(APIFunction(func_by_name(:vkCmdSetViewport), false), :(
        _cmd_set_viewport(command_buffer, viewports::AbstractArray)::Cvoid = @dispatch device(command_buffer) vkCmdSetViewport(command_buffer, 0, pointer_length(viewports), viewports)
    ))

    test_ex(APIFunction(func_by_name(:vkCmdSetLineWidth), false), :(
        _cmd_set_line_width(command_buffer, line_width::Real)::Cvoid = @dispatch device(command_buffer) vkCmdSetLineWidth(command_buffer, line_width)
    ))

    test_ex(APIFunction(func_by_name(:vkDestroyInstance), false), :(
        _destroy_instance(instance; allocator = C_NULL)::Cvoid = @dispatch instance vkDestroyInstance(instance, allocator)
    ))

    test_ex(APIFunction(func_by_name(:vkMapMemory), false), :(
        function _map_memory(device, memory, offset::Integer, size::Integer; flags = 0)::ResultTypes.Result{Ptr{Cvoid},VulkanError}
            ppData = Ref{Ptr{Cvoid}}()
            @check @dispatch device vkMapMemory(device, memory, offset, size, flags, ppData)
            ppData[]
        end
    ))

    test_ex(APIFunction(func_by_name(:vkEnumeratePhysicalDeviceQueueFamilyPerformanceQueryCountersKHR), false), :(
        function _enumerate_physical_device_queue_family_performance_query_counters_khr(physical_device, queue_family_index::Integer)::ResultTypes.Result{Tuple{Vector{_PerformanceCounterKHR},Vector{_PerformanceCounterDescriptionKHR}},VulkanError}
            pCounterCount = Ref{UInt32}()
            @repeat_while_incomplete begin
                @check @dispatch instance(physical_device) vkEnumeratePhysicalDeviceQueueFamilyPerformanceQueryCountersKHR(physical_device, queue_family_index, pCounterCount, C_NULL, C_NULL)
                pCounters = Vector{VkPerformanceCounterKHR}(undef, pCounterCount[])
                pCounterDescriptions = Vector{VkPerformanceCounterDescriptionKHR}(undef, pCounterCount[])
                @check @dispatch instance(physical_device) vkEnumeratePhysicalDeviceQueueFamilyPerformanceQueryCountersKHR(physical_device, queue_family_index, pCounterCount, pCounters, pCounterDescriptions)
            end
            (from_vk.(_PerformanceCounterKHR, pCounters), from_vk.(_PerformanceCounterDescriptionKHR, pCounterDescriptions))
        end
    ))

    test_ex(APIFunction(func_by_name(:vkGetPipelineCacheData), false), :(
        function _get_pipeline_cache_data(device, pipeline_cache)::ResultTypes.Result{Tuple{UInt,Ptr{Cvoid}},VulkanError}
            pDataSize = Ref{UInt}()
            @repeat_while_incomplete begin
                @check @dispatch device vkGetPipelineCacheData(device, pipeline_cache, pDataSize, C_NULL)
                pData = Libc.malloc(pDataSize[])
                @check @dispatch device vkGetPipelineCacheData(device, pipeline_cache, pDataSize, pData)
                if _return_code == VK_INCOMPLETE
                    Libc.free(pData)
                end
            end
            (pDataSize[], pData)
        end
    ))

    test_ex(APIFunction(func_by_name(:vkWriteAccelerationStructuresPropertiesKHR), false), :(
        _write_acceleration_structures_properties_khr(device, acceleration_structures::AbstractArray, query_type::QueryType, data_size::Integer, data::Ptr{Cvoid}, stride::Integer)::ResultTypes.Result{Result,VulkanError} = @check @dispatch device vkWriteAccelerationStructuresPropertiesKHR(device, pointer_length(acceleration_structures), acceleration_structures, query_type, data_size, data, stride)
    ))

    test_ex(APIFunction(func_by_name(:vkGetQueryPoolResults), false), :(
        _get_query_pool_results(device, query_pool, first_query::Integer, query_count::Integer, data_size::Integer, data::Ptr{Cvoid}, stride::Integer; flags = 0)::ResultTypes.Result{Result,VulkanError} = @check @dispatch device vkGetQueryPoolResults(device, query_pool, first_query, query_count, data_size, data, stride, flags)
    ))

    test_ex(APIFunction(func_by_name(:vkGetFenceStatus), false), :(
        _get_fence_status(device, fence)::ResultTypes.Result{Result,VulkanError} = @check(@dispatch(device, vkGetFenceStatus(device, fence)))
    ))

    test_ex(APIFunction(func_by_name(:vkGetSwapchainImagesKHR), false), :(
        function _get_swapchain_images_khr(device, swapchain)::ResultTypes.Result{Vector{Image},VulkanError}
            pSwapchainImageCount = Ref{UInt32}()
            @repeat_while_incomplete begin
                @check @dispatch device vkGetSwapchainImagesKHR(device, swapchain, pSwapchainImageCount, C_NULL)
                pSwapchainImages = Vector{VkImage}(undef, pSwapchainImageCount[])
                @check @dispatch device vkGetSwapchainImagesKHR(device, swapchain, pSwapchainImageCount, pSwapchainImages)
            end
            Image.(pSwapchainImages, identity, device)
        end
    ))

    @testset "Automatic reconstruction of create infos" begin
        test_ex(APIFunction(create_func(:vkCreateInstance), false), :(
            _create_instance(enabled_layer_names::AbstractArray, enabled_extension_names::AbstractArray; allocator = C_NULL, next = C_NULL, flags = 0, application_info = C_NULL) = _create_instance(_InstanceCreateInfo(enabled_layer_names, enabled_extension_names; next, flags, application_info); allocator)
        ))

        test_ex(APIFunction(create_func(:vkCreateDebugReportCallbackEXT), true), :(
            _create_debug_report_callback_ext(instance, pfn_callback::FunctionPtr, fptr_create::FunctionPtr, fptr_destroy::FunctionPtr; allocator = C_NULL, next = C_NULL, flags = 0, user_data = C_NULL) = _create_debug_report_callback_ext(instance, _DebugReportCallbackCreateInfoEXT(pfn_callback; next, flags, user_data), fptr_create, fptr_destroy; allocator)
        ),
        )
    end
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction