text:
using Random, Tar
import Tar_jll

Base.@kwdef struct DockerExecutor <: SandboxExecutor
    label::String = Random.randstring(10)
    privileges::Symbol = :privileged
end

function cleanup(exe::DockerExecutor)
    success(`docker system prune --force --filter=label=$(docker_image_label(exe))`)
end

Base.show(io::IO, exe::DockerExecutor) = write(io, "Docker Executor")

function executor_available(::Type{DockerExecutor}; verbose::Bool = false)
    # don't even try to exec if it doesn't exist
    if Sys.which("docker") === nothing
        if verbose
            @info("No `docker` command found; docker unavailable")
        end
        return false
    end

    # Return true if we can `docker ps`; if we can't, then there's probably a permissions issue
    if !success(`docker ps`)
        if verbose
            @warn("Unable to run `docker ps`; perhaps you're not in the `docker` group?")
        end
        return false
    end
    return with_executor(DockerExecutor) do exe
        return probe_executor(exe; test_read_only_map=true, test_read_write_map=true, verbose)
    end
end

timestamps_path() = joinpath(@get_scratch!("docker_timestamp_hashes"), "path_timestamps.toml")
function load_timestamps()
    path = timestamps_path()
    if !isfile(path)
        return Dict()
    end
    try
        return TOML.parsefile(path)
    catch e
        @error("couldn't load $(path)", exception=e)
        return Dict()
    end
end

function save_timestamp(image_name::String, timestamp::Float64)
    timestamps_toml = entry = sprint() do io
        timestamps = load_timestamps()
        timestamps[image_name] = timestamp
        TOML.print(io, timestamps)
    end
    open(timestamps_path(), write=true) do io
        write(io, timestamps_toml)
    end
end

docker_image_name(root_path::String, uid::Cint, gid::Cint) = "sandbox_rootfs:$(string(Base._crc32c(root_path), base=16))-$(uid)-$(gid)"
docker_image_label(exe::DockerExecutor) = string("org.julialang.sandbox.jl=", exe.label)
function should_build_docker_image(root_path::String, uid::Cint, gid::Cint)
    # If the image doesn't exist at all, always return true
    image_name = docker_image_name(root_path, uid, gid)
    if !success(`docker image inspect $(image_name)`)
        return true
    end

    # If this image has been built before, compare its historical timestamp to the current one
    curr_ctime = max_directory_ctime(root_path)
    prev_ctime = get(load_timestamps(), image_name, 0.0)
    return curr_ctime != prev_ctime
end

"""
    build_docker_image(root_path::String)

Docker doesn't like volume mounts within volume mounts, like we do with `sandbox`.
So we do things "the docker way", where we construct a rootfs docker image, then mount
things on top of that, with no recursive mounting.  We cut down on unnecessary work
somewhat by quick-scanning the directory for changes and only rebuilding if changes
are detected.
"""
function build_docker_image(root_path::String, uid::Cint, gid::Cint; verbose::Bool = false)
    image_name = docker_image_name(root_path, uid, gid)
    if should_build_docker_image(root_path, uid, gid)
        max_ctime = max_directory_ctime(root_path)
        if verbose
            @info("Building docker image $(image_name) with max timestamp $(max_ctime)")
        end

        # Build the docker image
        open(`docker import - $(image_name)`, "w", verbose ? stdout : devnull) do io
            # We need to record permissions, and therefore we cannot use Tar.jl.
            # Some systems (e.g. macOS) ship with a BSD tar that does not support the
            # `--owner` and `--group` command-line options. Therefore, if Tar_jll is
            # available, we use the GNU tar provided by Tar_jll. If Tar_jll is not available,
            # we fall back to the system tar.
            cd(root_path) do
            # tar = Tar_jll.is_available() ? Tar_jll.tar() : `tar`
            # run(pipeline(`$(tar) -c --owner=$(uid) --group=$(gid) .`, stdout=io))
            Tar_jll.tar() do tar
                run(pipeline(`$(tar) -c --owner=$(uid) --group=$(gid) .`, stdout=io))
            end
            end
        end

        # Record that we built it
        save_timestamp(image_name, max_ctime)
    end

    return image_name
end

function commit_previous_run(exe::DockerExecutor, image_name::String)
    ids = split(readchomp(`docker ps -a --filter label=$(docker_image_label(exe)) --format "{{.ID}}"`))
    if isempty(ids)
        return image_name
    end

    # We'll take the first docker container ID that we get, as its the most recent, and commit it.
    image_name = "sandbox_rootfs_persist:$(first(ids))"
    run(`docker commit $(first(ids)) $(image_name)`)
    return image_name
end

function build_executor_command(exe::DockerExecutor, config::SandboxConfig, user_cmd::Cmd)
    # Build the docker image that corresponds to this rootfs
    image_name = build_docker_image(config.read_only_maps["/"], config.uid, config.gid; verbose=config.verbose)

    if config.persist
        # If this is a persistent run, check to see if any previous runs have happened from
        # this executor, and if they have, we'll commit that previous run as a new image and
        # use it instead of the "base" image.
        image_name = commit_previous_run(exe, image_name)
    end

    # Begin building `docker` args
    if exe.privileges === :privileged  # this is the default
        # pros: allows you to do nested execution. e.g. the ability to run `Sandbox` inside `Sandbox`
        # cons: may allow processes inside the Docker container to access secure environment variables of processes outside the container
        privilege_args = String["--privileged"]
    elseif exe.privileges === :no_new_privileges
        # pros: may prevent privilege escalation attempts
        # cons: you won't be able to do nested execution
        privilege_args = String["--security-opt", "no-new-privileges"]
    elseif exe.privileges === :unprivileged
        # cons: you won't be able to do nested execution; privilege escalation may still work
        privilege_args = String[]
    else
        throw(ArgumentError("invalid value for exe.privileges: $(exe.privileges)"))
    end
    cmd_string = String["docker", "run", privilege_args..., "-i", "--label", docker_image_label(exe)]

    # If we're doing a fully-interactive session, tell it to allocate a psuedo-TTY
    if all(isa.((config.stdin, config.stdout, config.stderr), Base.TTY))
        push!(cmd_string, "-t")
    end

    # Start in the right directory
    append!(cmd_string, ["-w", config.pwd])

    # Add in read-only mappings (skipping the rootfs)
    for (dst, src) in config.read_only_maps
        if dst == "/"
            continue
        end
        append!(cmd_string, ["-v", "$(src):$(dst):ro"])
    end

    # Add in read-write mappings
    for (dst, src) in config.read_write_maps
        append!(cmd_string, ["-v", "$(src):$(dst)"])
    end

    # Apply environment mappings, first from `config`, next from `user_cmd`.
    for (k, v) in config.env
        append!(cmd_string, ["-e", "$(k)=$(v)"])
    end
    if user_cmd.env !== nothing
        for pair in user_cmd.env
            append!(cmd_string, ["-e", pair])
        end
    end

    # Add in entrypoint, if it is set
    if config.entrypoint !== nothing
        append!(cmd_string, ["--entrypoint", config.entrypoint])
    end

    # For each platform requested by `multiarch`, ensure its matching interpreter is registered,
    # but only if we're on Linux.  If we're on some other platform, like macOS where Docker is
    # implemented with a virtual machine, we just trust the docker folks to have set up the
    # relevant `binfmt_misc` mappings properly.
    if Sys.islinux()
        register_requested_formats!(config.multiarch_formats; verbose=config.verbose)
    end

    # Set the user and group
    append!(cmd_string, ["--user", "$(config.uid):$(config.gid)"])

    # Finally, append the docker image name user-requested command string
    push!(cmd_string, image_name)
    append!(cmd_string, user_cmd.exec)

    docker_cmd = Cmd(cmd_string)

    # If the user has asked that this command be allowed to fail silently, pass that on
    if user_cmd.ignorestatus
        docker_cmd = ignorestatus(docker_cmd)
    end

    return docker_cmd
end

sanitize_key(name) = replace(name, ':' => '-')

"""
    export_docker_image(image::String,
                        output_dir::String = <default scratch location>;
                        verbose::Bool = false,
                        force::Bool = false)

Exports the given docker image name to the requested output directory.  Useful
for pulling down a known good rootfs image from Docker Hub, for future use by
Sandbox executors.  If `force` is set to true, will overwrite a pre-existing
directory, otherwise will silently return.
"""
function export_docker_image(image_name::String,
                             output_dir::String = @get_scratch!("docker-$(sanitize_key(image_name))");
                             force::Bool = false,
                             verbose::Bool = false)
    if ispath(output_dir) && !isempty(readdir(output_dir))
        if force
            rmdir(output_dir; force=true, recursive=true)
        else
            if verbose
                @warn("Will not overwrite pre-existing directory $(output_dir)")
            end
            return output_dir
        end
    end

    # Get a container ID ready to be passed to `docker export`
    container_id = readchomp(`docker create $(image_name) /bin/true`)

    # Get the ID of that container (since we can't export by label, sadly)
    if isempty(container_id)
        if verbose
            @warn("Unable to create conatiner based on $(image_name)")
        end
        return nothing
    end

    # Export the container filesystem to a directory
    try
        mkpath(output_dir)
        open(`docker export $(container_id)`) do tar_io
            Tar.extract(tar_io, output_dir) do hdr
                # Skip known troublesome files
                return hdr.type ∉ (:chardev,)
            end
        end
    finally
        run(`docker rm -f $(container_id)`)
    end
    return output_dir
end

"""
    pull_docker_image(image::String,
                      output_dir::String = <default scratch location>;
                      verbose::Bool = false,
                      force::Bool = false)

Pulls and saves the given docker image name to the requested output directory.
Useful for pulling down a known good rootfs image from Docker Hub, for future use
by Sandbox executors.  If `force` is set to true, will overwrite a pre-existing
directory, otherwise will silently return.
"""
function pull_docker_image(image_name::String,
                           output_dir::String = @get_scratch!("docker-$(sanitize_key(image_name))");
                           force::Bool = false,
                           verbose::Bool = false)
    if ispath(output_dir) && !isempty(readdir(output_dir))
        if force
            rmdir(output_dir; force=true, recursive=true)
        else
            if verbose
                @warn("Will not overwrite pre-existing directory $(output_dir)")
            end
            return output_dir
        end
    end

    # Pull the latest version of the image
    try
        run(`docker pull $(image_name)`)
    catch
        if verbose
            @warn("Cannot pull", image_name)
        end
        return nothing
    end

    # Once the image is pulled, export it to given output directory
    return export_docker_image(image_name, output_dir; force, verbose)
end


text:
photo(; kwargs...) = requestimage("place/photo"; kwargs...)

"""
Returns a photo corresponding to `photoref`, with specified `maxwidth`.

When you get place information using a Place Details request, photo references
will be returned for relevant photographic content. The Nearby Search and Text
Search requests also return a single photo reference per place, when relevant.
Using the Photo service you can then access the referenced photos and resize the
image to the optimal size for your application.

https://developers.google.com/places/web-service/photos

## Parameters
* `key`: your google maps API key
* `photoref`: A string identifier that uniquely identifies a photo. Photo
    references are returned from either a Place Search or Place Details request.
* `maxwidth`: Specifies the maximum desired width (between 1 and 1600 pixels) of
    the image. If the image is smaller than `maxwidth`, the original image will
    be returned. If the image is larger, it will be scaled to match, restricted
    to its original aspect ratio.
"""
photo(key::String, photoref::String, maxwidth::Int=400; kwargs...) =
    photo(key=key, photoreference=photoref, maxwidth=maxwidth; kwargs...)

ext: jl


text:
module LevelKRacing

using AutomotiveDrivingModels
using AutoViz
using NearestNeighbors

import AutomotiveDrivingModels: get_actions!, observe!, get_name
import Base.rand
import PyPlot

export
    # HierarchicalRecedingHorizonController
    HRHC,
    curveDist,
    wrap_to_π,
    kdProject,
    generateObstacleMap,
    updateObstacleMap!,
    generateMotionMap,
    screenCollision,
    tailgateAvoidance,
    getSuccessorStates,
    loopProjectionKD,
    computeTrajectory,
    screenTrajectory,
    checkCollision,
    calculateObjective,
    plot_stϕ,
    plotHRHCInfo,
    plotObjectiveHorizon,
    plotSplineRoadway,

    # LevelKVisualizations
    DemoState,
    DemoMotionPrimitives,
    DemoIncreasedHorizonBehavior,
    DemoBuildingBlocks,
    DemoObjective,
    DemoTailgateAvoidance,
    DemoObserveActObjective,
    DemoObstacleMap,

    # SplineRaceWay
    Raceway,

    # SplineUtils
    ClosedB_Spline,
    B_SplineDerivative,
    ResampleSplineEven,
    GenSplineRoadway,
    PlotSplineRoadway

include(Pkg.dir("LevelKRacing","src","controllers","HierarchicalRecedingHorizonController.jl"))
include(Pkg.dir("LevelKRacing","src","visualization","LevelKVisualizations.jl"))
include(Pkg.dir("LevelKRacing","src","utils","SplineRaceWay.jl"))
include(Pkg.dir("LevelKRacing","src","utils","SplineUtils.jl"))

end


text:
using ProgressMeter, Plots, StatsBase, Statistics, LaTeXStrings, JLD

function RandomWalker(P,L,X0)
    X = X0
    Tₛ= 0
    while abs(X) != round(Int,L/2)
        X += sample([1,-1], Weights([P,1-P]))
        Tₛ+= 1
    end
    return Tₛ
end

PList = hcat(0:0.01:1)
X0List = hcat(-20:20)
runnum = 10000

Data = zeros(length(X0List),length(PList),2)
progress = Progress(length(X0List)*length(PList)*runnum)

for i in 1:length(X0List)
    X0 = X0List[i]
    for j in 1:length(PList)
        P = PList[j]
        stepruns = []
        for run in 1:runnum
            push!(stepruns,RandomWalker(P,40,X0))
            next!(progress)
            update!(progress)
        end
        Data[i,j,1] = mean(stepruns)
        Data[i,j,2] = std(stepruns)
    end
end

save("../../Data/Q5/Q5-LT.jld", "Data", Data)
Data=load("../../Data/Q5/Q5-LT.jld")["Data"]

P1 = heatmap(0:0.01:1,-20:20,Data[:,:,1], xaxis=nothing,legend = nothing, yaxis=nothing,title = L"Average\ Life\ time")
P2 = contour(0:0.01:1,-20:20,Data[:,:,1], xaxis=nothing, fill = true ,legend = nothing,ylabel = L"X_{init}",title = L"Average\ Life\ time")
P3 = heatmap(0:0.01:1,-20:20,Data[:,:,2], yaxis=nothing, legend = nothing, xlabel = L"P",title = L"STD\ of\ Life\ time")
P4 = contour(0:0.01:1,-20:20,Data[:,:,2], fill = true ,legend = nothing, xlabel = L"P",ylabel = L"X_{init}",title = L"STD\ of\ Life\ time")

Pmain1 = begin
    P_tot = [P2 P4;P1 P3]
    CBar = plot(heatmap((0:400).*ones(401,1),xlabel ="",
        legend=:none, xticks=:none, yticks=(1:40:401, string.(0:40:400))),
            heatmap((0:300).*ones(301,1),xlabel = L"range",
                legend=:none, xticks=:none, yticks=(1:30:301, string.(0:30:300))),
                    layout = (2,1))
    plot(P_tot...,CBar,
        layout = @layout[grid(2,2) a{0.035w}],
        plot_title = L"Life\ time\ of\ a\ Random\ Walker\ (%$runnum \ runs)",size = (800,600))
end
savefig(Pmain1,"../../Figs/Q5/Q5-cont.pdf")

X0Plot = []
for X0 in [11,21,31]
    PLT = begin
        plot(PList,Data[X0,:,1], linestyle = :dot, c = :black)
        scatter!(PList,Data[X0,:,1], yerr = Data[X0,:,2], c = :black, markersize = 1.75)
        scatter!(PList,Data[X0,:,1], ribbon = Data[X0,:,2], markersize = 1.5, c = :purple, markerstrokewidth=0.2)
        plot!(title = L"X_{init} = %$(X0List[X0])",
            xlabel = L"P", ylabel = L"Avg\ Life\ Time", legend = :none)
    end
    push!(X0Plot, PLT)
end

X0P = plot(X0Plot...,layout = (1,3))

PPlot = []
for P in [26,51,76]
    PLT = begin
        plot(X0List,Data[:,P,1], linestyle = :dot, c = :black)
        scatter!(X0List,Data[:,P,1], yerr = Data[:,P,2], c = :black, markersize = 1.75)
        scatter!(X0List,Data[:,P,1], ribbon = Data[:,P,2], markersize = 1.5, c = :steelblue, markerstrokewidth=0.2)
        plot!(title = L"P = %$(PList[P])",
            xlabel = L"X0", ylabel = L"Avg\ Life\ Time", legend = :none)
    end
    push!(PPlot, PLT)
end
PPlot[1]
PP = plot(PPlot...,layout = (1,3))

Pmain2 = plot(PP,X0P,layout = (2,1),size = (850,690), plot_title = L"Life\ time\ of\ a\ Random\ Walker\ (%$runnum \ runs)")
savefig(Pmain2,"../../Figs/Q5/Q5-scat.pdf")


text:
__precompile__(true)
module WinPath

import WinReg
import JSON
import Base

export in, out, getConfig

const BASE = WinReg.HKEY_LOCAL_MACHINE
const PATH = "System\\ControlSet001\\Control\\Session Manager\\Environment"
const VALUENAME = "Path"
const CONFIG_PATH = "C:\\WinPath\\config.json"

struct Config
    root::AbstractString
    outPath::AbstractString
    inPath::AbstractString
    function Config(root::AbstractString, outPath::AbstractString, inPath::AbstractString)
        new(root, outPath, inPath)
    end
    function Config(jsonDict::Dict{String, Any})
        root = jsonDict["root"]
        outPath = jsonDict["outPath"]
        inPath = jsonDict["inPath"]
        new(root, outPath, inPath)
    end
end

function Base.show(io::IO, config::Config)
    print(io, string("root: ", config.root, "\n",
                     "outPath: ", config.outPath, "\n",
                     "inPath: ", config.inPath))
end

function config()
    jsonDict = JSON.parsefile(CONFIG_PATH)
    config = Config(jsonDict)
    if !isdir(config.root)
        mkdir(config.root)
    end
    return config
end

function getPath()
    return WinReg.querykey(BASE, PATH, VALUENAME)
end

function out()
    path = getPath()
    list = split(path, ";")
    list = filter!(x -> x != "", list)
    list = sort(list)
    outData = join(list, "\n")
    open(CONFIG.outPath, "w+") do file
        write(file, outData)
    end
    return nothing
end

function in()
    list = String[]
    open(CONFIG.outPath, "r") do outfile
        for ln in readlines(outfile)
            push!(list, ln)
        end
        list = sort(list)
        list = join(list, ";")
        open(CONFIG.inPath, "w+") do infile
            write(infile, list)
        end
    end
    return nothing
end

function getConfig()
    println(CONFIG)
end

const CONFIG = config()

end


text:
## Task切换不会使用调用栈
## Task可以以任意次序切换，而不像函数必须返回上层调用

## 这种流程控制使得解决特定问题更简便
## 例如生产-消费问题

function producer(c::Channel)
    put!(c, "start")
    for n=1:4
        put!(c, 2n)
    end
    put!(c, "stop")
end;
#  不需显式关闭Channel, Channel 对 Task 的绑定同时也意味着 Channel 的生命周期与绑定的 Task 一致。当 Task 结束时，Channel 对象会自动关闭。多个 Channel 可以绑定到一个 Task，反之亦然。

chan = Channel(producer);

println(take!(chan))
println(take!(chan))
println(take!(chan))

# 用for-loop读取剩余的值
for x in chan
    println(x)
end

# 多数 Task 切换是在等待如 I/O 请求的事件，由 Julia Base 里的调度器执行。调度器维持一个可运行 Task 的队列，并执行一个事件循环，来根据例如收到消息等外部事件来重启 Task

max_stars_repo_name: SenZhangAI/example

max_issues_repo_name: SenZhangAI/example

max_forks_repo_name: SenZhangAI/example


text:
# include("../src/plotSol.jl")
# include("../src/OrientationConversion.jl")
# include("../src/simulate.jl")
#
# using JLD
# clearconsole()
#
# savedSol = load("ModRob1dSecSimRK45NoTol.jld") # DifferentialEquations
# # sol = load("ModRob1dSecRK4.jld") # Native RK
# solFinal = savedSol["solFinal"]
# tSim = savedSol["tSim"]

## Assigning solutions
solBody1 = solFinal[1]
solGimbal1 = solFinal[2]
solProp11 = solFinal[3]
solProp12 = solFinal[4]
solBox = solFinal[5]
solBody2 = solFinal[6]
solGimbal2 = solFinal[7]
solProp21 = solFinal[8]
solProp22 = solFinal[9]

## Plotting
jointLoc2 = Matrix{Float64}(undef,length(tSim),3)
jointLoc3 = Matrix{Float64}(undef,length(tSim),3)
jointLoc4 = Matrix{Float64}(undef,length(tSim),3)
jointLoc5 = Matrix{Float64}(undef,length(tSim),3)
jointLoc6 = Matrix{Float64}(undef,length(tSim),3)
jointLoc7 = Matrix{Float64}(undef,length(tSim),3)
jointLoc8 = Matrix{Float64}(undef,length(tSim),3)
jointLoc9 = Matrix{Float64}(undef,length(tSim),3)

ωBody1 = Matrix{Float64}(undef,length(tSim),3)
ωGimbal1 = Matrix{Float64}(undef,length(tSim),3)
ωProp11 = Matrix{Float64}(undef,length(tSim),3)
ωProp12 = Matrix{Float64}(undef,length(tSim),3)
ωGimbal1InBody1 = Matrix{Float64}(undef,length(tSim),3)
ωRel1 = Matrix{Float64}(undef,length(tSim),3)
ωRel2 = Matrix{Float64}(undef,length(tSim),3)
ωRel3 = Matrix{Float64}(undef,length(tSim),3)

ωBox = Matrix{Float64}(undef,length(tSim),3)

ωBody2 = Matrix{Float64}(undef,length(tSim),3)
ωGimbal2 = Matrix{Float64}(undef,length(tSim),3)
ωProp21 = Matrix{Float64}(undef,length(tSim),3)
ωProp22 = Matrix{Float64}(undef,length(tSim),3)
ωGimbal2InBody2 = Matrix{Float64}(undef,length(tSim),3)
ωRel4 = Matrix{Float64}(undef,length(tSim),3)
ωRel5 = Matrix{Float64}(undef,length(tSim),3)
ωRel6 = Matrix{Float64}(undef,length(tSim),3)
for i=1:length(tSim)
    Body1.dcm = quat2dcm(solBody1.β[i,:])
    Gimbal1.dcm = quat2dcm(solGimbal1.β[i,:])
    Prop11.dcm = quat2dcm(solProp11.β[i,:])
    Prop12.dcm = quat2dcm(solProp12.β[i,:])
    Box.dcm = quat2dcm(solBox.β[i,:])

    ωBox[i,:] = angVel(solBox.β[i,:],solBox.βdot[i,:])
    ωBody1[i,:] = angVel(solBody1.β[i,:],solBody1.βdot[i,:])
    ωGimbal1[i,:] = angVel(solGimbal1.β[i,:],solGimbal1.βdot[i,:])
    ωProp11[i,:] = angVel(solProp11.β[i,:],solProp11.βdot[i,:])
    ωProp12[i,:] = angVel(solProp12.β[i,:],solProp12.βdot[i,:])

    ωGimbal1InBody1[i,:] = quat2dcm(solBody1.β[i,:])*
                         transpose(quat2dcm(solGimbal1.β[i,:]))*
                         ωGimbal1[i,:]
    ωRel1[i,:] = ωGimbal1InBody1[i,:] - ωBody1[i,:]
    ωRel2[i,:] = quat2dcm(solGimbal1.β[i,:])*
                 transpose(quat2dcm(solProp11.β[i,:]))*
                 ωProp11[i,:] - ωGimbal1[i,:]
    ωRel3[i,:] = quat2dcm(solGimbal1.β[i,:])*
                 transpose(quat2dcm(solProp12.β[i,:]))*
                 ωProp12[i,:] - ωGimbal1[i,:]

    jointLoc2[i,:] = solBody1.r[i,:] + transpose(Body1.dcm)*rjBody  -                    solGimbal1.r[i,:] -                     transpose(Gimbal1.dcm)*rjGimbal1
    jointLoc3[i,:] = solGimbal1.r[i,:] + transpose(Gimbal1.dcm)*rjGimbal2 - solProp11.r[i,:] - transpose(Prop11.dcm)*rjProp1
    jointLoc4[i,:] = solGimbal1.r[i,:] +
                     transpose(Gimbal1.dcm)*rjGimbal3 -
                     solProp12.r[i,:] - transpose(Prop12.dcm)*rjProp2
    jointLoc5[i,:] = solBody1.r[i,:] + transpose(Body1.dcm)*rjCac1 - solBox.r[i,:] - transpose(Box.dcm)*rjBox1
end
for i=1:length(tSim)
     Body2.dcm = quat2dcm(solBody2.β[i,:])
     Gimbal2.dcm = quat2dcm(solGimbal2.β[i,:])
     Prop21.dcm = quat2dcm(solProp21.β[i,:])
     Prop22.dcm = quat2dcm(solProp22.β[i,:])
     Box.dcm = quat2dcm(solBox.β[i,:])

     ωBody2[i,:] = angVel(solBody2.β[i,:],solBody2.βdot[i,:])
     ωGimbal2[i,:] = angVel(solGimbal2.β[i,:],solGimbal2.βdot[i,:])
     ωProp21[i,:] = angVel(solProp21.β[i,:],solProp21.βdot[i,:])
     ωProp22[i,:] = angVel(solProp22.β[i,:],solProp22.βdot[i,:])

     ωGimbal2InBody2[i,:] = quat2dcm(solBody2.β[i,:])*
                          transpose(quat2dcm(solGimbal2.β[i,:]))*
                          ωGimbal2[i,:]
     ωRel4[i,:] = ωGimbal2InBody2[i,:] - ωBody2[i,:]
     ωRel5[i,:] = quat2dcm(solGimbal2.β[i,:])*
                  transpose(quat2dcm(solProp21.β[i,:]))*
                  ωProp21[i,:] - ωGimbal2[i,:]
     ωRel6[i,:] = quat2dcm(solGimbal2.β[i,:])*
                  transpose(quat2dcm(solProp22.β[i,:]))*
                  ωProp22[i,:] - ωGimbal2[i,:]

     jointLoc7[i,:] = solBody2.r[i,:] + transpose(Body2.dcm)*rjBody  -                    solGimbal2.r[i,:] -                     transpose(Gimbal2.dcm)*rjGimbal1
     jointLoc8[i,:] = solGimbal2.r[i,:] + transpose(Gimbal2.dcm)*rjGimbal2 - solProp21.r[i,:] - transpose(Prop21.dcm)*rjProp1
     jointLoc9[i,:] = solGimbal2.r[i,:] +
                      transpose(Gimbal2.dcm)*rjGimbal3 -
                      solProp22.r[i,:] - transpose(Prop22.dcm)*rjProp2
     jointLoc6[i,:] = solBody2.r[i,:] + transpose(Body2.dcm)*rjCac2 - solBox.r[i,:] - transpose(Box.dcm)*rjBox2
end

## Plotting
# Joint Locations
plotPos(tSim,jointLoc2)
plotPos(tSim,jointLoc3)
plotPos(tSim,jointLoc4)
plotPos(tSim,jointLoc5)
plotPos(tSim,jointLoc6)
plotPos(tSim,jointLoc7)
plotPos(tSim,jointLoc8)
plotPos(tSim,jointLoc9)

plotErrNorm(tSim,solBody1.β)
plotErrNorm(tSim,solGimbal1.β)
plotErrNorm(tSim,solProp11.β)
plotErrNorm(tSim,solProp12.β)
plotErrNorm(tSim,solBody2.β)
plotErrNorm(tSim,solGimbal2.β)
plotErrNorm(tSim,solProp21.β)
plotErrNorm(tSim,solProp22.β)

plotAngVel(tSim,ωBody1)
plotAngVel(tSim,ωGimbal1)
plotAngVel(tSim,ωProp11)
plotAngVel(tSim,ωProp12)
plotAngVel(tSim,ωRel1)
plotAngVel(tSim,ωRel2)
plotAngVel(tSim,ωRel3)

plotAngVel(tSim,ωBox)

plotAngVel(tSim,ωBody2)
plotAngVel(tSim,ωGimbal2)
plotAngVel(tSim,ωProp21)
plotAngVel(tSim,ωProp22)
plotAngVel(tSim,ωRel4)
plotAngVel(tSim,ωRel5)
plotAngVel(tSim,ωRel6)

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
module MessagePassingAssemblersTest

using Test
using ForneyLab
using ForneyLab: assembleBreaker!, assembleClamp!, assembleInferenceAlgorithm!, assemblePosteriorFactor!, assembleSchedule!, assembleInitialization!, assembleMarginalTable!, condense, flatten, setTargets!, messagePassingSchedule

@testset "assembleClamp!" begin
    g = FactorGraph()
    nd = Clamp(Variable(), 1.0)
    assembleClamp!(nd, ProbabilityDistribution)
    @test nd.dist_or_msg == ProbabilityDistribution
end

@testset "assembleBreaker!" begin
    breaker_entry = ScheduleEntry()
    assembleBreaker!(breaker_entry, GaussianMeanPrecision, ())
    @test breaker_entry.family == GaussianMeanPrecision
    @test breaker_entry.initialize == true
    @test breaker_entry.dimensionality == ()

    breaker_entry = ScheduleEntry()
    assembleBreaker!(breaker_entry, Union{Gamma, Wishart}, (1,1))
    @test breaker_entry.family == Union{Gamma, Wishart}
    @test breaker_entry.initialize == true
    @test breaker_entry.dimensionality == (1,1)
end

@testset "assembleSchedule!" begin
    fg = FactorGraph()
    @RV x ~ GaussianMeanPrecision(0.0, 1.0)
    GaussianMeanPrecision(x, 0.0, 1.0)
    pfz = PosteriorFactorization()
    pf = PosteriorFactor(fg)
    setTargets!(pf, pfz, target_variables=Set{Variable}([x]))
    pf.schedule = messagePassingSchedule(pf)
    algo = InferenceAlgorithm(pfz)
    algo.target_to_marginal_entry = Dict()
    algo.interface_to_schedule_entry = ForneyLab.interfaceToScheduleEntry(algo)
    assembleSchedule!(pf)
    @test pf.schedule[3].message_update_rule == ForneyLab.SPGaussianMeanPrecisionOutNPP
    @test pf.schedule[6].message_update_rule == ForneyLab.SPGaussianMeanPrecisionOutNPP
end

@testset "assembleInitialization!" begin
    # Expectation propagation
    fg = FactorGraph()
    @RV x ~ GaussianMeanPrecision(0.0, 1.0)
    @RV y ~ Probit(x)
    placeholder(y, :y)
    pfz = PosteriorFactorization()
    pf = PosteriorFactor(fg)
    pf.ep_sites = Set{Interface}([fg.nodes[:probit_1].i[:in1]])
    setTargets!(pf, pfz, target_variables=Set{Variable}([x]))
    pf.schedule = messagePassingSchedule(pf)
    algo = InferenceAlgorithm(pfz)
    algo.target_to_marginal_entry = Dict()
    algo.interface_to_schedule_entry = ForneyLab.interfaceToScheduleEntry(algo)
    assembleSchedule!(pf)
    assembleInitialization!(pf)
    @test pf.schedule[5].message_update_rule == ForneyLab.EPProbitIn1PG
    @test pf.schedule[3].initialize

    # Nonlinear
    f(z) = z
    fg = FactorGraph()
    @RV x ~ GaussianMeanPrecision(0.0, 1.0)
    @RV y ~ Nonlinear{Unscented}(x, g=f)
    GaussianMeanPrecision(y, 0.0, 1.0)

    pfz = PosteriorFactorization()
    pf = PosteriorFactor(fg)
    setTargets!(pf, pfz, target_variables=Set{Variable}([x]))
    pf.schedule = messagePassingSchedule(pf)
    algo = InferenceAlgorithm(pfz)
    algo.target_to_marginal_entry = Dict()
    algo.interface_to_schedule_entry = ForneyLab.interfaceToScheduleEntry(algo)
    assembleSchedule!(pf)
    assembleInitialization!(pf)
    @test pf.schedule[7].message_update_rule == ForneyLab.SPNonlinearUTIn1GG
    @test pf.schedule[3].initialize

    # Optimize
    fg = FactorGraph()
    @RV x ~ GaussianMeanPrecision(0.0, 1.0)
    placeholder(x, :x)
    pfz = PosteriorFactorization()
    pf = PosteriorFactor(fg)
    setTargets!(pf, pfz, target_variables=Set{Variable}([x]))
    pf.schedule = messagePassingSchedule(pf)
    algo = InferenceAlgorithm(pfz)
    algo.target_to_marginal_entry = Dict()
    algo.interface_to_schedule_entry = ForneyLab.interfaceToScheduleEntry(algo)
    assembleSchedule!(pf)
    assembleInitialization!(pf)
    @test pf.schedule[3].initialize
end

@testset "assembleMarginalTable!" begin
    # Nothing rule
    fg = FactorGraph()
    @RV x ~ GaussianMeanPrecision(0.0, 1.0)
    pfz = PosteriorFactorization()
    pf = PosteriorFactor(fg)
    setTargets!(pf, pfz, target_variables=Set{Variable}([x]))
    pf.schedule = messagePassingSchedule(pf)
    pf.marginal_table = marginalTable(x)
    algo = InferenceAlgorithm(pfz)
    algo.interface_to_schedule_entry = ForneyLab.interfaceToScheduleEntry(algo)
    algo.target_to_marginal_entry = ForneyLab.targetToMarginalEntry(algo)
    assembleMarginalTable!(pf)
    @test pf.marginal_table[1].marginal_update_rule == Nothing
    @test pf.marginal_table[1].marginal_id == :x
    @test pf.marginal_table[1].inbounds == [pf.schedule[3]]

    # Product rule
    fg = FactorGraph()
    @RV x ~ GaussianMeanPrecision(0.0, 1.0)
    GaussianMeanPrecision(x, 0.0, 1.0)
    pfz = PosteriorFactorization()
    pf = PosteriorFactor(fg)
    setTargets!(pf, pfz, target_variables=Set{Variable}([x]))
    pf.schedule = messagePassingSchedule(pf)
    pf.marginal_table = marginalTable(x)
    algo = InferenceAlgorithm(pfz)
    algo.interface_to_schedule_entry = ForneyLab.interfaceToScheduleEntry(algo)
    algo.target_to_marginal_entry = ForneyLab.targetToMarginalEntry(algo)
    assembleMarginalTable!(pf)
    @test pf.marginal_table[1].marginal_update_rule == ForneyLab.Product
    @test pf.marginal_table[1].marginal_id == :x
    @test pf.marginal_table[1].inbounds == [pf.schedule[3], pf.schedule[6]]

    # Marginal rule
    fg = FactorGraph()
    @RV x ~ GaussianMeanPrecision(0.0, 1.0)
    @RV y ~ GaussianMeanPrecision(x, 1.0)
    GaussianMeanPrecision(y, 0.0, 1.0)
    pfz = PosteriorFactorization([x,y], ids=[:XY])
    pf = pfz.posterior_factors[:XY]
    setTargets!(pf, pfz, external_targets=true)
    pf.schedule = messagePassingSchedule(pf)
    pf.marginal_table = marginalTable(pf)
    algo = InferenceAlgorithm(pfz)
    algo.interface_to_schedule_entry = ForneyLab.interfaceToScheduleEntry(algo)
    algo.target_to_marginal_entry = ForneyLab.targetToMarginalEntry(algo)
    assembleMarginalTable!(pf)
    @test pf.marginal_table[3].marginal_update_rule == ForneyLab.MGaussianMeanPrecisionGGD
    @test pf.marginal_table[3].marginal_id == :y_x
    @test length(pf.marginal_table[3].inbounds) == 3
end

@testset "assemblePosteriorFactor!" begin
    fg = FactorGraph()
    @RV x ~ GaussianMeanPrecision(0.0, 1.0)
    GaussianMeanPrecision(x, 0.0, 1.0)
    pfz = PosteriorFactorization()
    pf = PosteriorFactor(fg)
    setTargets!(pf, pfz, target_variables=Set{Variable}([x]))
    schedule = messagePassingSchedule(pf)
    pf.schedule = condense(flatten(schedule))
    pf.marginal_table = marginalTable(x)
    algo = InferenceAlgorithm()
    algo.interface_to_schedule_entry = ForneyLab.interfaceToScheduleEntry(algo)
    algo.target_to_marginal_entry = ForneyLab.targetToMarginalEntry(algo)
    assemblePosteriorFactor!(pf)
    @test pf.schedule[1].schedule_index == 1
    @test pf.schedule[1].message_update_rule == ForneyLab.SPGaussianMeanPrecisionOutNPP
    @test pf.schedule[1].inbounds == [nothing, fg.nodes[:clamp_1], fg.nodes[:clamp_2]]
    @test pf.schedule[2].schedule_index == 2
    @test pf.schedule[2].message_update_rule == ForneyLab.SPGaussianMeanPrecisionOutNPP
    @test pf.schedule[2].inbounds == [nothing, fg.nodes[:clamp_3], fg.nodes[:clamp_4]]
    @test pf.marginal_table[1].marginal_id == :x
    @test pf.marginal_table[1].marginal_update_rule == ForneyLab.Product
    @test pf.marginal_table[1].inbounds == [schedule[3], schedule[6]]
end

@testset "assembleInferenceAlgorithm!" begin
    fg = FactorGraph()
    @RV x ~ GaussianMeanPrecision(0.0, 1.0)
    GaussianMeanPrecision(x, 0.0, 1.0)
    pfz = PosteriorFactorization()
    pf = PosteriorFactor(fg)
    setTargets!(pf, pfz, target_variables=Set{Variable}([x]))
    schedule = messagePassingSchedule(pf)
    pf.schedule = condense(flatten(schedule))
    pf.marginal_table = marginalTable(x)
    algo = InferenceAlgorithm(pfz)
    assembleInferenceAlgorithm!(algo)

    @test length(algo.interface_to_schedule_entry) == 2
    @test length(algo.target_to_marginal_entry) == 1
end

end # module


text:
module HelloJulia

export go

using IJulia

const ROOT = joinpath(@__DIR__, "..")
const NOTEBOOKS = joinpath(ROOT, "notebooks")

go() = begin
    notebook(dir=NOTEBOOKS)
end

end # module


text:
export ray_class_group, narrow_class_group, MapRayClassGrp

add_verbose_scope(:RayFacElem)
add_assert_scope(:RayFacElem)

###############################################################################
#
#  Map Type
#
###############################################################################

mutable struct MapRayClassGrp <: Map{GrpAbFinGen, FacElemMon{Hecke.NfOrdIdlSet}, HeckeMap, MapRayClassGrp}
  header::Hecke.MapHeader{GrpAbFinGen, FacElemMon{Hecke.NfOrdIdlSet}}
  defining_modulus::Tuple{NfOrdIdl, Vector{InfPlc}}
  fact_mod::Dict{NfOrdIdl, Int} #The factorization of the finite part of the defining modulus

  gens::Tuple{Vector{NfOrdIdl}, Vector{GrpAbFinGenElem}}

  #Dictionaries to cache preimages. Used in the action on the ray class group
  prime_ideal_preimage_cache::Dict{NfOrdIdl, GrpAbFinGenElem}
  prime_ideal_cache::Vector{NfOrdIdl}

  clgrpmap::MapClassGrp
  powers::Vector{Tuple{NfOrdIdl, NfOrdIdl}}
  groups_and_maps::Vector{Tuple{GrpAbFinGen, GrpAbFinGenToAbsOrdQuoRingMultMap}}
  disc_log_inf_plc::Dict{InfPlc, GrpAbFinGenElem} #The infinite places and the corresponding discrete logarithm.
  gens_mult_grp_disc_log::Vector{Tuple{NfOrdElem, GrpAbFinGenElem}}

  function MapRayClassGrp()
    z = new()
    z.prime_ideal_preimage_cache = Dict{NfOrdIdl, GrpAbFinGenElem}()
    return z
  end
end

defining_modulus(mR) = mR.defining_modulus

################################################################################
#
#  Function that stores the principal generators element of the powers
#  of the generators in the class group map
#
################################################################################

function __assure_princ_gen(c::Hecke.ClassGrpCtx{SMat{fmpz}}, nquo::Int)
  module_trafo_assure(c.M)
  C = c.dl_data[3]
  OK = order(c)
  s = c.dl_data[1]
  if length(c.dl_data) == 4
    T = c.dl_data[4]
  else
    T = inv(c.dl_data[2])
    c.dl_data = (s, c.dl_data[2], C, T)
  end
  RelHnf = c.M.basis
  gens = c.FB.ideals
  rels = vcat(c.R_gen, c.R_rel)
  trafo = c.M.trafo
  res = Tuple{FacElem{NfOrdIdl, NfOrdIdlSet}, FacElem{nf_elem, AnticNumberField}}[]
  diff = ppio(Int(C.snf[end]), nquo)[2]
  diff_gens = ncols(T) - ngens(C)
  for i = 1:ngens(C)
    if nquo != -1
      if iscoprime(C.snf[i], nquo)
        continue
      end
      el = diff*sub(T, i+diff_gens:i+diff_gens, 1:ncols(T))
      ex = Int(ppio(Int(C.snf[i]), nquo)[1])
    else
      el = sub(T, i+diff_gens:i+diff_gens, 1:ncols(T))
      ex = Int(C.snf[i])
    end
    els_r = Tuple{Int, fmpz}[]
    DI = Dict{NfOrdIdl, fmpz}()
    for j = 1:ncols(el)
      if !iszero(el[1, j])
        add_to_key!(DI, gens[j+s-1], el[1, j])
        push!(els_r, (j+s-1, ex*el[1, j]))
      end
    end
    r = sparse_row(FlintZZ, els_r, sort = false)
    sol, d = solve_ut(RelHnf, r)
    @assert isone(d)
    rs = zeros(fmpz, c.M.bas_gens.r + c.M.rel_gens.r)

    for (p,v) in sol
      rs[p] = v
    end

    for i in length(trafo):-1:1
      apply_right!(rs, trafo[i])
    end

    e = FacElem(rels, rs)
    e = reduce_mod_units([e], _get_UnitGrpCtx_of_order(OK))[1]
    dd = e.fac
    for i = dd.idxfloor:length(dd.vals)
      if dd.slots[i] == 0x1 && iszero(dd.vals[i])
        dd.count -= 1
        dd.slots[i] = 0x0
      end
    end
    I = FacElem(DI)
    J, a = reduce_ideal(I)
    inv!(a)
    pow!(a, ex)
    mul!(e, e, a)
    push!(res, (FacElem(Dict(J => 1)), e))
  end
  return res
end

function _assure_princ_gen(mC::MapClassGrp)
  if isdefined(mC, :princ_gens)
    return nothing
  end
  C = domain(mC)
  OK = order(codomain(mC))
  K = nf(OK)
  if order(domain(mC)) == 1
    res1 = Vector{Tuple{FacElem{NfOrdIdl, NfOrdIdlSet}, FacElem{nf_elem, AnticNumberField}}}()
    if ngens(domain(mC)) == 1
      push!(res1, (FacElem(Dict(ideal(OK, 1) => 1)), FacElem(Dict(K(1) => 1))))
    end
    mC.princ_gens = res1
    return nothing
  end
  try
    c = _get_ClassGrpCtx_of_order(OK)
    if !isdefined(c, :dl_data)
      throw(AccessorNotSetError())
    end
    res = __assure_princ_gen(c, mC.quo)
    @hassert :RayFacElem 1 isconsistent(mC, res)
    mC.princ_gens = res
    return nothing
  catch e
    if !(isa(e, AccessorNotSetError))
      rethrow(e)
    end
    c = _get_ClassGrpCtx_of_order(OK.lllO)
    reslll = __assure_princ_gen(c, mC.quo)
    res = Vector{Tuple{FacElem{NfOrdIdl, NfOrdIdlSet}, FacElem{nf_elem, AnticNumberField}}}(undef, length(reslll))
    for i = 1:length(res)
      fe = Dict{NfOrdIdl, fmpz}()
      for (k, v) in reslll[i][1]
        fe[IdealSet(OK)(k)] = v
      end
      res[i] = (FacElem(fe), reslll[i][2])
    end
    @hassert :RayFacElem 1 isconsistent(mC, res)
    mC.princ_gens = res
    return nothing
  end
end

function isconsistent(mC, res)
  C = domain(mC)
  OK = order(codomain(mC))
  for i = 1:length(res)
    I = evaluate(res[i][1]).num
    if mC\I != C[i]
      return false
    end
    e = evaluate(res[i][2])
    J = ideal(OK, OK(e))
    if I^Int(C.snf[i]) != J
      return false
    end
  end
  return true
end

################################################################################
#
#  Class group as a ray class group
#
################################################################################

function class_as_ray_class(C::GrpAbFinGen, mC::MapClassGrp, exp_class::Function,  m::NfOrdIdl, expo::Int)

  O = order(m)
  X = abelian_group(rels(C))

  if expo != -1
    Q, mQ = quo(C, expo, false)
    local disclog1
    let Q = Q, mC = mC, mQ = mQ, X = X
      function disclog1(J::NfOrdIdl)
        return mQ(mC\(J))
      end

      function disclog1(J::FacElem{NfOrdIdl, NfOrdIdlSet})
        a = X[0]
        for (f, k) in J.fac
          a += k*disclog(f)
        end
        return a
      end
    end

    local expo_map
    let mQ = mQ, exp_class = exp_class
      function expo_map(el::GrpAbFinGenElem)
        return exp_class(mQ\el)
      end
    end

    mp1 = Hecke.MapRayClassGrp()
    mp1.header = Hecke.MapHeader(Q, FacElemMon(parent(m)), expo_map, disclog1)
    mp1.fact_mod = Dict{NfOrdIdl, Int}()
    mp1.defining_modulus = (m, InfPlc[])
    mp1.clgrpmap = mC
    return Q, mp1
  end

  local disclog
  let X = X, mC = mC
    function disclog(J::NfOrdIdl)
      return X((mC\J).coeff)
    end

    function disclog(J::FacElem{NfOrdIdl, NfOrdIdlSet})
      a = X[0]
      for (f, k) in J.fac
        a += k*disclog(f)
      end
      return a
    end
  end

  mp = Hecke.MapRayClassGrp()
  mp.header = Hecke.MapHeader(X, FacElemMon(parent(m)), exp_class, disclog)
  mp.fact_mod = Dict{NfOrdIdl, Int}()
  mp.defining_modulus = (m, InfPlc[])
  mp.clgrpmap = mC
  return X, mp
end

function empty_ray_class(m::NfOrdIdl)
  O = order(parent(m))
  X = abelian_group(Int[])

  local exp
  let O = O
    function exp(a::GrpAbFinGenElem)
      return FacElem(Dict(ideal(O,1) => fmpz(1)))
    end
  end

  local disclog
  let X = X
    function disclog(J::Union{NfOrdIdl, FacElem{NfOrdIdl}})
      return id(X)
    end
  end

  mp = Hecke.MapRayClassGrp()
  mp.header = Hecke.MapHeader(X, FacElemMon(parent(m)) , exp, disclog)
  mp.defining_modulus = (m, InfPlc[])
  return X,mp

end

##############################################################################
#
#  Functions for the evaluation of factored elements
#
###############################################################################

#
#  Multiple elements evaluation
#
function fac_elems_eval(p::NfOrdIdl, q::NfOrdIdl, elems::Vector{FacElem{nf_elem, AnticNumberField}}, exponent::fmpz)
  return _eval_quo(elems, p, q, exponent)
end

function _preproc(el::FacElem{nf_elem, AnticNumberField}, exponent::fmpz)
  K = base_ring(el)
  OK = maximal_order(K)
  Qx = parent(K.pol)
  x = Dict{nf_elem, fmpz}()
  for (f, k) in el
    l = mod(k,exponent)
    if iszero(l)
      continue
    end
    if f in OK
      add_to_key!(x, f, l)
    else
      d = denominator(f, OK)
      add_to_key!(x, K(d), exponent-l)
      n = d*f
      c = numerator(content(Qx(n)))
      if isone(c)
        add_to_key!(x, n, l)
      else
        add_to_key!(x, divexact(n, c), l)
        add_to_key!(x, K(c), l)
      end
    end
  end
  if !isempty(x)
    return FacElem(x)
  else
    return FacElem(Dict(K(1)=> 1))
  end
end

function _preproc(p::NfOrdIdl, el::FacElem{nf_elem, AnticNumberField}, exponent::fmpz)
  O = order(p)
  K = nf(O)
  Qx = parent(K.pol)
  x = Dict{nf_elem, fmpz}()
  P = minimum(p, copy = false)
  for (f, k) in el
    l = mod(k,exponent)
    if iszero(l)
      continue
    end
    n = numerator(f)
    d = denominator(f)
    if !isone(d)
      com, uncom = ppio(d, P)
      if !isone(uncom)
        add_to_key!(x, K(mod(uncom, P)), exponent-l)
      end
      if !isone(com)
        e, b = ispower(com)
        add_to_key!(x, K(b), -e*l)
      end
    end
    c = numerator(content(Qx(n)))
    if isone(c)
      add_to_key!(x, n, l)
    else
      add_to_key!(x, divexact(n, c), l)
      com, uncom = ppio(c, P)
      if !isone(uncom)
        add_to_key!(x, K(mod(uncom, P)), l)
      end
      if !isone(com)
        e, b = ispower(com)
        add_to_key!(x, K(b), e*l)
      end
    end
  end
  if !isempty(x)
    return FacElem(x)
  else
    return FacElem(Dict(K(1)=> 1))
  end
end

function _preproc(p::NfOrdIdl, elems::Vector{FacElem{nf_elem, AnticNumberField}}, exponent::fmpz)
  newelems = FacElem{nf_elem, AnticNumberField}[_preproc(p, x, exponent) for x in elems]
  return newelems
end

function _powermod(a::nf_elem, i::Int, p::fmpz)
  if iszero(i)
    return one(parent(a))
  elseif isone(i)
    b = mod(a, p)
    return b
  else
    bit = ~((~UInt(0)) >> 1)
    while (UInt(bit) & i) == 0
      bit >>= 1
    end
    z = deepcopy(a)
    bit >>= 1
    while bit != 0
      mul!(z, z, z)
      z = mod(z, p)
      if (UInt(bit) & i) != 0
        mul!(z, z, a)
        z = mod(z, p)
      end
      bit >>= 1
    end
    return z
  end
end

function _ev_quo(Q, mQ, elems, p, exponent, multiplicity::Int)
  el = elem_type(Q)[one(Q) for i = 1:length(elems)]
  anti_uni = anti_uniformizer(p)
  powers = Dict{Int, nf_elem}()
  powers[1] = anti_uni
  O = order(p)
  F, mF = ResidueField(O, p)
  for i = 1:length(elems)
    J = elems[i]
    vp = fmpz(0)
    for (f, k1) in J
      k = mod(k1, exponent)
      if iszero(k)
        continue
      end
      if isinteger(f)
        inte = numerator(coeff(f, 0))
        vpp, np = remove(inte, minimum(p, copy = false))
        mul!(el[i], el[i], Q(np)^k)
        vp += vpp*k
        continue
      end
      el1 = O(f, false)
      if !iszero(mF(el1))
        if !isone(k)
          mul!(el[i], el[i], mQ(el1)^k)
        else
          mul!(el[i], el[i], mQ(el1))
        end
        continue
      end
      val = valuation(f, p)
      if haskey(powers, val)
        act_el = O(powers[val]*f, false)
      else
        exp_av = div(multiplicity*val, ramification_index(p))
        anti_val = _powermod(anti_uni, val, minimum(p)^(exp_av+1))
        powers[val] = anti_val
        act_el = O(anti_val*f, false)
      end
      if !isone(k)
        mul!(el[i], el[i], mQ(act_el)^k)
      else
        mul!(el[i], el[i], mQ(act_el))
      end
    end
    vp = mod(vp, exponent)
    if !iszero(vp)
      if haskey(powers, ramification_index(p))
        eli = minimum(p, copy = false)*powers[ramification_index(p)]
      else
        powers[ramification_index(p)] = _powermod(anti_uni, ramification_index(p), minimum(p)^(multiplicity+1))
        eli = minimum(p, copy = false)*powers[ramification_index(p)]
      end
      if isone(vp)
        mul!(el[i], el[i], mQ(O(eli, false)))
      else
        mul!(el[i], el[i], mQ(O(eli, false))^vp)
      end
    end
  end
  return NfOrdElem[mQ\el[i] for i=1:length(el)]
end

function _eval_quo(elems::Vector{FacElem{nf_elem, AnticNumberField}}, p::NfOrdIdl, q::NfOrdIdl, exponent::fmpz)
  O = order(p)
  if p == q
    if fits(Int, p.minimum)
      Q, mQ = ResidueFieldSmall(O, p)
      return _ev_quo(Q, mQ, elems, p, exponent, 1)
    else
      Q, mQ = ResidueField(O, p)
      return _ev_quo(Q, mQ, elems, p, exponent, 1)
    end
  else
    Q, mQ = quo(O, q)
    mult = Int(clog(norm(q), norm(p)))
    return _ev_quo(Q, mQ, elems, p, exponent, mult)
  end
end

################################################################################
#
#  n-part of the class group
#
################################################################################

function iscoprime(a, b)
  return isone(gcd(a, b))
end

function n_part_class_group(mC::Hecke.MapClassGrp, n::Integer)
  O = order(codomain(mC))
  C = domain(mC)
  @assert issnf(C)
  K = nf(O)
  if iscoprime(exponent(C), n)
    G = abelian_group(fmpz[])
    local exp1
    let O = O
      function exp1(a::GrpAbFinGenElem)
        return ideal(O, one(O))
      end
    end

    local disclog1
    let G = G
      function disclog1(I::NfOrdIdl)
        return G[0]
      end
    end

    mp=Hecke.MapClassGrp()
    mp.quo = n
    mp.header=Hecke.MapHeader(G, mC.header.codomain, exp1, disclog1)
    mp.princ_gens = Tuple{FacElem{NfOrdIdl, NfOrdIdlSet}, FacElem{nf_elem, AnticNumberField}}[]
    return G, mp
  end

  ind = findfirst(x -> !iscoprime(x, n), C.snf)
  diff = ppio(exponent(C), fmpz(n))[2]

  invariants = fmpz[ppio(x, fmpz(n))[1] for x in C.snf[ind:end]]
  #filter!(!isone, invariants)
  G = abelian_group(invariants)
  local exp2
  let O = O, G = G
    function exp2(a::GrpAbFinGenElem)
      new_coeff = zero_matrix(FlintZZ, 1, ngens(C))
      for i = 1:ngens(G)
        new_coeff[1, i+ind-1] = a[i]*diff
      end
      return mC(C(new_coeff))
    end
  end

  local disclog2
  let G = G, mC = mC, C = C, diff = diff
    idiff = invmod(diff, exponent(G))
    function disclog2(I::NfOrdIdl)
      if I.is_principal == 1
        return id(G)
      end
      x=idiff*(mC\I)
      y = zero_matrix(FlintZZ, 1, ngens(G))
      for i=ind:ngens(C)
        y[1,i-ind+1]=x.coeff[1,i]
      end
      return GrpAbFinGenElem(G, y)
    end
  end

  mp = Hecke.MapClassGrp()
  mp.header = Hecke.MapHeader(G, mC.header.codomain, exp2, disclog2)
  mp.quo = Int(exponent(G))
  if isdefined(mC, :princ_gens)
    princ_gens = Vector{Tuple{FacElem{NfOrdIdl, NfOrdIdlSet}, FacElem{nf_elem, AnticNumberField}}}(undef, length(vect))
    for i = 1:length(princ_gens)
      princ_gens[i] = (mC.princ_gens[ind+i-1][1]^diff, mC.princ_gens[ind+i-1][2])
    end
    mp.princ_gens = princ_gens
  end
  return G, mp
end

################################################################################
#
#  Make positive
#
################################################################################

#makes the element x positive at all the embeddings adding a multiple of a
#TODO: Do this properly!
function make_positive(x::NfOrdElem, a::fmpz)
  els = conjugates_real(elem_in_nf(x))
  m = fmpz(0)
  for i=1:length(els)
    if ispositive(els[i])
      continue
    end
    y = abs(lower_bound(els[i], fmpz))
    mu = div(y, a)
    m = max(m, mu+1)
  end
  #@hassert :RayFacElem 1 iscoprime(ideal(parent(x),x), ideal(parent(x), a))
  #@hassert :RayFacElem 1 iscoprime(ideal(parent(x),x+fmpz(m)*a), ideal(parent(x), a))
  @hassert :RayFacElem 1 istotally_positive(x+m*a)
  el_to_add = m*a
  return x+el_to_add
end

###################################################################################
#
#  Narrow Class Group
#
###################################################################################

@doc Markdown.doc"""
    narrow_class_group(O::NfOrd) -> GrpAbFinGen, Map

Computes the narrow (or strict) class group of $O$, ie. the group of invertable
ideals modulo principal ideals generated by elements that are
positive at all real places.
"""
function narrow_class_group(O::NfOrd)
  @assert ismaximal_known_and_maximal(O)
  K = nf(O)
  plc = real_places(K)
  return ray_class_group(ideal(O, 1), real_places(K))
end

###################################################################################
#
#  Ray Class Group
#
###################################################################################

function ray_class_group(O::NfOrd, D::Dict{NfOrdIdl, Int}, inf_plc::Vector{InfPlc} = InfPlc[]; n_quo::Int = -1, GRH::Bool = true)
  I = ideal(O, 1)
  minI = fmpz(1)
  for (p, v) in D
    q = p^v
    I *= q
    minI = lcm(minI, minimum(q))
  end
  I.minimum = minI
  return ray_class_group(I, inf_plc, GRH = GRH, n_quo = n_quo, lp = D)
end

#
# We compute the group using the sequence U -> (O/m)^* _> Cl^m -> Cl -> 1
#
@doc Markdown.doc"""
    ray_class_group(m::NfOrdIdl, inf_plc::Vector{InfPlc}; n_quo::Int, lp::Dict{NfOrdIdl, Int}) -> GrpAbFinGen, MapRayClassGrp

Given an ideal $m$ and a set of infinite places of $K$,
this function returns the corresponding ray class group as an abstract group $\mathcal {Cl}_m$ and a map going
from the group into the group of ideals of $K$ that are coprime to $m$.
If `n_quo` is set, it will return the group modulo `n_quo`. The factorization of $m$ can be given with the keyword argument `lp`.
"""
function ray_class_group(m::NfOrdIdl, inf_plc::Vector{InfPlc} = Vector{InfPlc}(); GRH::Bool = true, n_quo::Int = -1, lp::Dict{NfOrdIdl, Int} = factor(m))

  O = order(m)
  K = nf(O)

  C, mC = class_group(O, GRH = GRH)
  expC = exponent(C)
  diffC = fmpz(1)
  if n_quo != -1
    C, mC = n_part_class_group(mC, n_quo)
    diffC = divexact(expC, exponent(C))
    expC = exponent(C)
  end
  if isone(m) && isempty(inf_plc)
    local exp_c
    let mC = mC
      function exp_c(a::GrpAbFinGenElem)
        return FacElem(Dict(mC(a) => 1))
      end
    end
    return class_as_ray_class(C, mC, exp_c, m, n_quo)
  end
  _assure_princ_gen(mC)

  exp_class, Kel = find_coprime_representatives(mC, m, lp)

  if n_quo != -1
    powers = Vector{Tuple{NfOrdIdl, NfOrdIdl}}()
    quo_rings = Tuple{NfOrdQuoRing, Hecke.AbsOrdQuoMap{NfAbsOrd{AnticNumberField,nf_elem},NfAbsOrdIdl{AnticNumberField,nf_elem},NfAbsOrdElem{AnticNumberField,nf_elem}}}[]
    groups_and_maps = Tuple{GrpAbFinGen, Hecke.GrpAbFinGenToAbsOrdQuoRingMultMap{NfAbsOrd{AnticNumberField,nf_elem},NfAbsOrdIdl{AnticNumberField,nf_elem},NfAbsOrdElem{AnticNumberField,nf_elem}}}[]
    for (pp, vv) in lp
      dtame = Dict{NfOrdIdl, Int}()
      dwild = Dict{NfOrdIdl, Int}()
      npp = norm(pp)
      qq = ideal(O, 1)
      if !iscoprime(npp-1, n_quo)
        dtame[pp] = 1
        qq = pp
      end
      if vv > 1 && !iscoprime(npp, n_quo)
        dwild[pp] = vv
        qq = pp^vv
      end
      if isempty(dtame) && isempty(dwild)
        continue
      end
      push!(powers, (pp, qq))
      Q, mQ = quo(O, qq)
      if pp == qq
        Q.factor = dtame
      else
        Q.factor = dwild
      end
      push!(quo_rings, (Q, mQ))
      push!(groups_and_maps, _mult_grp_mod_n(quo_rings[end][1], dtame, dwild, n_quo))
    end
  else
    powers = Tuple{NfOrdIdl, NfOrdIdl}[(p, p^v) for (p, v) in lp]
    quo_rings = Tuple{NfOrdQuoRing, Hecke.AbsOrdQuoMap{NfAbsOrd{AnticNumberField,nf_elem},NfAbsOrdIdl{AnticNumberField,nf_elem},NfAbsOrdElem{AnticNumberField,nf_elem}}}[quo(O, q) for (p, q) in powers]
    groups_and_maps = Tuple{GrpAbFinGen, Hecke.GrpAbFinGenToAbsOrdQuoRingMultMap{NfAbsOrd{AnticNumberField,nf_elem},NfAbsOrdIdl{AnticNumberField,nf_elem},NfAbsOrdElem{AnticNumberField,nf_elem}}}[_multgrp(x[1], true) for x in quo_rings]
  end
  if isempty(groups_and_maps)
    nG = 0
    expon = fmpz(1)
  else
    nG = sum(ngens(x[1]) for x in groups_and_maps)
    expon = lcm([exponent(x[1]) for x in groups_and_maps])
  end

  if n_quo == -1 || iseven(n_quo)
    p = InfPlc[ x for x in inf_plc if isreal(x) ]
  else
    p = InfPlc[]
  end
  H, eH, lH = infinite_primes_map(O, p, m)
  expon = lcm(expon, exponent(H))

  U, mU = unit_group_fac_elem(O, GRH = GRH)

  # We construct the relation matrix and evaluate units and relations with the class group in the quotient by m
  # Then we compute the discrete logarithms

  if n_quo == -1
    R = zero_matrix(FlintZZ, ngens(C)+nG+ngens(H)+ngens(U), ngens(H)+ngens(C)+nG)
  else
    R = zero_matrix(FlintZZ, 2*(ngens(C)+nG+ngens(H))+ngens(U), ngens(C)+ngens(H)+nG)
    for i = 1:ncols(R)
      R[i+ngens(C)+nG+ngens(H)+ngens(U), i] = n_quo
    end
  end
  for i=1:ngens(C)
    R[i,i] = C.snf[i]
  end
  ind = 1
  for s = 1:length(quo_rings)
    G = groups_and_maps[s][1]
    @assert issnf(G)
    for i = 1:ngens(G)
      R[i+ngens(C)+ind-1, i+ngens(C)+ind-1] = G.snf[i]
    end
    ind += ngens(G)
  end
  for i = 1:ngens(H)
    R[ngens(C)+nG+i, ngens(C)+nG+i] = 2
  end

  @vprint :RayFacElem 1 "Collecting elements to be evaluated; first, units \n"
  tobeeval1 = Vector{FacElem{nf_elem, AnticNumberField}}(undef, ngens(U)+ngens(C))
  for i = 1:ngens(U)
    tobeeval1[i] = mU(U[i])
  end
  for i = 1:ngens(C)
    tobeeval1[i+ngens(U)] = mC.princ_gens[i][2]*(FacElem(Dict(Kel[i] => C.snf[i])))
  end
  tobeeval = _preproc(m, tobeeval1, expon)

  ind = 1
  for i = 1:length(groups_and_maps)
    exp_q = gcd(expon, norm(powers[i][2])- divexact(norm(powers[i][2]), norm(powers[i][1])))
    @vtime :RayFacElem 3 evals = fac_elems_eval(powers[i][1], powers[i][2], tobeeval, exp_q)
    Q = quo_rings[i][1]
    mG = groups_and_maps[i][2]
    for j = 1:ngens(U)
      a = (mG\Q(evals[j])).coeff
      for s = 1:ncols(a)
        R[j+nG+ngens(H)+ngens(C), ngens(C)+s+ind-1] = a[1, s]
      end
    end

    for j = 1:ngens(C)
      a = (mG\Q(evals[j+ngens(U)])).coeff
      for s = 1:ncols(a)
        R[j, ngens(C)+ind+s-1] = -a[1, s]
      end
    end
    ind += ngens(groups_and_maps[i][1])
  end

  if !isempty(p)
    for j = 1:ngens(U)
      a = lH(tobeeval1[j]).coeff
      for s = 1:ncols(a)
        R[j+nG+ngens(C)+ngens(H), ngens(C)+ind-1+s] = a[1, s]
      end
    end
    for j = 1:ngens(C)
      a = lH(tobeeval1[j+ngens(U)]).coeff
      for s = 1:ncols(a)
        R[j, ngens(C)+ind-1+s] = -a[1, s]
      end
    end
  end

  X = abelian_group(R)
  if n_quo != -1
    X.exponent = n_quo
  end

  local disclog
  let X = X, mC = mC, C = C, exp_class = exp_class, powers = powers, groups_and_maps = groups_and_maps, quo_rings = quo_rings, lH = lH, diffC = diffC, n_quo = n_quo, m = m, expon = expon
    invd = invmod(fmpz(diffC), expon)
    # Discrete logarithm
    function disclog(J::FacElem{NfOrdIdl, NfOrdIdlSet})
      @vprint :RayFacElem 1 "Disc log of element $J \n"
      a = id(X)
      for (f, k) in J
        a += k*disclog(f)
      end
      return a
    end

    function disclog(J::NfOrdIdl)
      @hassert :RayFacElem 1 iscoprime(J, m)
      if isone(J)
        @vprint :RayFacElem 1 "J is one \n"
        return id(X)
      end
      coeffs = zero_matrix(FlintZZ, 1, ngens(X))
      if J.is_principal == 1 && isdefined(J, :princ_gen)
        z = FacElem(Dict(J.princ_gen.elem_in_nf => diffC))
      else
        L = mC\J
        for i = 1:ngens(C)
          coeffs[1, i] = L[i]
        end
        @vprint :RayFacElem 1 "Disc log of element J in the Class Group: $(L.coeff) \n"
        s = exp_class(L)
        inv!(s)
        add_to_key!(s.fac, J, 1)
        pow!(s, diffC)
        @vprint :RayFacElem 1 "This ideal is principal: $s \n"
        z = principal_generator_fac_elem(s)
      end
      ii = 1
      z1 = _preproc(m, z, expon)
      for i = 1:length(powers)
        P, Q = powers[i]
        exponq = gcd(expon, norm(Q)-divexact(norm(Q), norm(P)))
        el = fac_elems_eval(P, Q, FacElem{nf_elem, AnticNumberField}[z1], exponq)
        y = (invd*(groups_and_maps[i][2]\quo_rings[i][1](el[1]))).coeff
        for s = 1:ncols(y)
          coeffs[1, ii-1+ngens(C)+s] = y[1, s]
        end
        ii += ngens(groups_and_maps[i][1])
      end
      if !isempty(p)
        b = lH(z).coeff
        for s = 1:ncols(b)
          coeffs[1, ii-1+s+ngens(C)] = b[1, s]
        end
      end
      return GrpAbFinGenElem(X, coeffs)
    end
  end

  Dgens = Tuple{NfOrdElem, GrpAbFinGenElem}[]
  ind = 1
  #For the exponential map and other purposes, we need generators of the full multiplicative group
  #In particular, we need the idempotents...
  for i = 1:length(powers)
    P, Q = powers[i]
    mG = groups_and_maps[i][2]
    J = ideal(O, 1)
    minJ = fmpz(1)
    mins = fmpz(1)
    for (PP, vPP) in lp
      if minimum(PP, copy = false) != minimum(P, copy = false)
        mins = lcm(mins, minimum(PP, copy = false)^vPP)
        continue
      end
      if PP != P
        Jm = PP^vPP
        J *= Jm
        minJ = lcm(minJ, minimum(Jm))
      end
    end
    J.minimum = minJ
    i1, i2 = idempotents(Q, J)
    if !isone(mins)
      d, u1, v1 = gcdx(minimum(Q, copy = false), mins)
      i1 = i1*(u1*minimum(Q, copy = false) + v1*mins) + u1*minimum(Q, copy = false) *i2
      i2 = v1*mins*i2
    end
    gens = mG.generators
    if isempty(p)
      if haskey(mG.tame, P)
        gens_tame = mG.tame[P].generators
        for s = 1:length(gens_tame)
          gens_tame[s] = gens_tame[s]*i2 + i1
        end
        mG.tame[P].generators = gens_tame
      end
      if haskey(mG.wild, P)
        gens_wild = mG.wild[P].generators
        for s = 1:length(gens_wild)
          gens_wild[s] = gens_wild[s]*i2 + i1
        end
        mG.wild[P].generators = gens_wild
      end
      for s = 1:length(gens)
        push!(Dgens, (gens[s].elem*i2+i1, X[ngens(C)+ind-1+s]))
      end

    else
      if haskey(mG.tame, P)
        gens_tame = mG.tame[P].generators
        for s = 1:length(gens_tame)
          gens_tame[s] = make_positive(gens_tame[s]*i2 + i1, minimum(m, copy = false))
        end
        mG.tame[P].generators = gens_tame
      end
      if haskey(mG.wild, P)
        gens_wild = mG.wild[P].generators
        for s = 1:length(gens_wild)
          gens_wild[s] = make_positive(gens_wild[s]*i2 + i1, minimum(m, copy = false))
        end
        mG.wild[P].generators = gens_wild
      end
      for s = 1:length(gens)
        elgen = make_positive(gens[s].elem*i2 + i1, minimum(m, copy = false))
        push!(Dgens, (elgen, X[ngens(C)+ind-1+s]))
      end
    end
    ind += length(gens)
  end

  local expo
  let C = C, O = O, groups_and_maps = groups_and_maps, exp_class = exp_class, eH = eH, H = H, K = K, Dgens = Dgens, X = X, p = p
    function expo(a::GrpAbFinGenElem)
      b = GrpAbFinGenElem(C, sub(a.coeff, 1:1, 1:ngens(C)))
      res = exp_class(b)
      for i = 1:nG
        if !iszero(a.coeff[1, ngens(C)+i])
          add_to_key!(res.fac, ideal(O, Dgens[i][1]), a.coeff[1, ngens(C)+i])
        end
      end
      for i = 1:length(p)
        if !iszero(a.coeff[i+nG+ngens(C)])
          add_to_key!(res.fac, ideal(O, O(1+eH(H[i]))), 1)
        end
      end
      return res
    end
  end

  ind = 1
  for i = 1:length(powers)
    mG = groups_and_maps[i][2]
    for (prim, mprim) in mG.tame
      dis = zero_matrix(FlintZZ, 1, ngens(X))
      to_be_c = mprim.disc_log.coeff
      for i = 1:length(to_be_c)
        dis[1, ind-1+i+ngens(C)] = to_be_c[1, i]
      end
      mprim.disc_log = GrpAbFinGenElem(X, dis)
    end
    ind += ngens(domain(mG))
  end

  disc_log_inf = Dict{InfPlc, GrpAbFinGenElem}()
  for i = 1:length(p)
    eldi = zero_matrix(FlintZZ, 1,  ngens(X))
    eldi[1, ngens(X) - length(p) + i] = 1
    disc_log_inf[p[i]] = GrpAbFinGenElem(X, eldi)
  end

  mp = MapRayClassGrp()
  mp.header = Hecke.MapHeader(X, FacElemMon(parent(m)), expo, disclog)
  mp.fact_mod = lp
  mp.defining_modulus = (m, inf_plc)
  mp.powers = powers
  mp.groups_and_maps = groups_and_maps
  mp.disc_log_inf_plc = disc_log_inf
  mp.gens_mult_grp_disc_log = Dgens
  mp.clgrpmap = mC
  return X, mp

end

##################################################################################
#
#  Ray Class Group over QQ
#
##################################################################################

function ray_class_groupQQ(O::NfOrd, modulus::Int, inf_plc::Bool, n_quo::Int)

  R=ResidueRing(FlintZZ, modulus, cached=false)
  U, mU = unit_group_mod(R, n_quo)
  U.exponent = n_quo
  if inf_plc
    function disc_log1(I::NfOrdIdl)
      @assert gcd(minimum(I),modulus)==1
      i = Int(mod(I.minimum, modulus))
      return mU\(R(i))
    end

    function expon1(a::GrpAbFinGenElem)
      x=mU(a)
      return FacElem(Dict{NfOrdIdl, fmpz}(ideal(O,lift(x)) => 1))
    end

    mp=Hecke.MapRayClassGrp()
    mp.header = Hecke.MapHeader(U, FacElemMon(parent(ideal(O,1))) , expon1, disc_log1)
    mp.defining_modulus = (ideal(O, modulus), real_places(nf(O)))
    return U, mp

  elseif isodd(n_quo)

    function disc_log2(I::NfOrdIdl)
      @assert gcd(minimum(I),modulus)==1
      i=Int(mod(I.minimum, modulus))
      return mU\(R(i))
    end

    function expon2(a::GrpAbFinGenElem)
      x=mU(a)
      return FacElem(Dict{NfOrdIdl, fmpz}(ideal(O,lift(x)) => 1))
    end

    mp = Hecke.MapRayClassGrp()
    mp.header = Hecke.MapHeader(U, FacElemMon(parent(ideal(O,1))) , expon2, disc_log2)
    mp.defining_modulus = (ideal(O, modulus), InfPlc[])
    return U,mp

  else

    Q, mQ = quo(U, GrpAbFinGenElem[mU\(R(-1))], false)

    function disc_log(I::NfOrdIdl)
      i=Int(mod(minimum(I), modulus))
      return mQ(mU\(R(i)))
    end

    function expon(a::GrpAbFinGenElem)
      x=mU(mQ\a)
      return FacElem(Dict{NfOrdIdl, fmpz}(ideal(O,x) => 1))
    end

    mp=Hecke.MapRayClassGrp()
    mp.header = Hecke.MapHeader(Q, FacElemMon(parent(ideal(O,1))) , expon, disc_log)
    mp.defining_modulus = (ideal(O, modulus), InfPlc[])
    return Q,mp

  end

end

##################################################################################
#
#  Action of the Galois Group on the Ray Class Group
#
##################################################################################

#changes the generators of the multiplicative group so that they are coprime to a
function change_into_coprime(mR::MapRayClassGrp, a::fmpz)
  m = minimum(mR.defining_modulus[1])
  com, uncom = ppio(a, m)
  if uncom == 1
    return nothing
  end
  _, s, t = gcdx(uncom, m)
  gens = mR.gens_mult_grp_disc_log
  for i = 1:length(gens)
    gens[i] = (gens[i][1]*s*uncom+m*t, gens[i][2])
  end
  mR.gens_mult_grp_disc_log = gens
  return nothing
end

#  Find convenient ideals that generate the ray class groups
function find_gens(mR::MapRayClassGrp; coprime_to::fmpz = fmpz(-1))

  O = order(codomain(mR))
  R = domain(mR)
  m = mR.defining_modulus[1]
  mm = minimum(m)
  if coprime_to != -1
    mm = lcm(mm, coprime_to)
  end
  mm = lcm(mm, discriminant(EquationOrder(nf(O))))
  if isdefined(mR, :gens)
    if coprime_to == -1
      return mR.gens[1], mR.gens[2]
    else
      found = false
      sR = GrpAbFinGenElem[]
      lp = NfOrdIdl[]
      for i = 1:length(mR.gens[1])
        if iscoprime(mR.gens[1][i], coprime_to)
          push!(lp, mR.gens[1][i])
          push!(sR, mR.gens[2][i])
        else
          found = true
        end
      end
      if !found
        return lp, sR
      end
    end
  else
    sR = GrpAbFinGenElem[]
    lp = NfOrdIdl[]
  end

  q, mq = quo(R, sR, false)
  s, ms = snf(q)

  #  First, generators of the multiplicative group.
  #  If the class group is trivial, they are enough
  if isdefined(mR, :fact_mod) && !isempty(mR.fact_mod)
    if coprime_to != -1
      # First, I change them in order to be coprime to coprime_to
      change_into_coprime(mR, coprime_to)
    end
    gens_m = mR.gens_mult_grp_disc_log
    for i = 1:length(gens_m)
      el_q = mq(gens_m[i][2])
      if iszero(el_q)
        continue
      end
      if isprime_power(order(s))
        el_in_s = ms\el_q
        found = false
        for j = 1:ngens(s)
          if iscoprime(el_in_s[j], s.snf[j])
            found = true
            break
          end
        end
        if found
          push!(sR, gens_m[i][2])
          push!(lp, ideal(O, gens_m[i][1]))
          q, mq = quo(R, sR, false)
          s, ms = snf(q)
          if order(s) == 1
            if !isdefined(mR, :gens)
              mR.gens = (lp, sR)
            end
            return lp, sR
          end
        end
      else
        push!(sR, gens_m[i][2])
        push!(lp, ideal(O, gens_m[i][1]))
        q, mq = quo(R, sR, false)
        s, ms = snf(q)
        if order(s) == 1
          if !isdefined(mR, :gens)
            mR.gens = (lp, sR)
          end
          return lp, sR
        end
      end
    end
  end

  if !isempty(mR.defining_modulus[2])
    S, ex, lo = infinite_primes_map(O, mR.defining_modulus[2], m)
    for i=1:length(mR.defining_modulus[2])
      pl = mR.defining_modulus[2][i]
      if iscomplex(pl)
        continue
      end
      f = mR.disc_log_inf_plc[pl]
      if iszero(mq(f))
        continue
      else
        el_in_snf = ms\mq(f)
        found = false
        for j = 1:ngens(s)
          if iscoprime(s.snf[j], el_in_snf[j])
            found = true
            break
          end
        end
        if found
          I = ideal(O, 1+ex(S[i])) # Careful! This way I am assuming an ordering!
          push!(sR, f)
          push!(lp, I)
          q, mq = quo(R, sR, false)
          s, ms = snf(q)
          if order(s)==1
	          if !isdefined(mR, :gens)
              mR.gens = (lp, sR)
            end
            return lp, sR
          end
        end
      end
    end
  end

  #This means that the class group is non trivial. I need primes generating the class group
  mC = mR.clgrpmap
  C = domain(mC)
  primes_class_group = Vector{NfOrdIdl}()
  disc_log_primes_class_grp = Vector{GrpAbFinGenElem}()
  if isdefined(mC, :small_gens)
    for x in mC.small_gens
      if !iscoprime(mm, minimum(x, copy = false))
        continue
      end
      push!(primes_class_group, x)
      push!(disc_log_primes_class_grp, mC\x)
    end
  end
  Q1, mQ1 = quo(C, disc_log_primes_class_grp, false)
  S1, mS1 = snf(Q1)
  p1 = fmpz(2)
  while order(S1) != 1
    p1 = next_prime(p1)
    if !iscoprime(p1, mm)
      continue
    end
    lP = prime_decomposition(O, p1)
    for i = 1:length(lP)-1
      x = lP[i][1]
      f = mC\x
      el_in_snf = mS1\(mQ1(f))
      found = false
      for j = 1:ngens(S1)
        if iscoprime(S1.snf[j], el_in_snf[j])
          found = true
          break
        end
      end
      if found
        push!(primes_class_group, x)
        push!(disc_log_primes_class_grp, f)
        Q1, mQ1 = quo(C, disc_log_primes_class_grp, false)
        S1, mS1 = snf(Q1)
        if order(S1) == 1
          break
        end
      end
    end
  end
  for i = 1:length(primes_class_group)
    push!(lp, primes_class_group[i])
    push!(sR, mR\primes_class_group[i])
  end
  if !isdefined(mR, :gens)
    mR.gens = (lp, sR)
  end
  return lp, sR
end

function induce_action(mR::Union{MapRayClassGrp, MapClassGrp}, Aut::Vector{Hecke.NfToNfMor}, mp::GrpAbFinGenMap = id_hom(domain(mR)))
  R = domain(mR)
  G = Vector{GrpAbFinGenMap}(undef, length(Aut))
  if isempty(Aut)
    return G
  end
  K = domain(Aut[1])
  O = order(codomain(mR))

  #  Instead of applying the automorphisms to the elements given by mR, I choose small primes
  #  generating the group and study the action on them. In this way, I take advantage of the cache of the
  #  class group map

  Igens, IPgens, subs, IPsubs = find_gens_for_action(mR)
  genstot = vcat(subs, IPsubs)
  for k = 1:length(Aut)
    images = Vector{GrpAbFinGenElem}(undef, length(Igens)+length(IPgens))
    for i=1:length(Igens)
      J = induce_image(Aut[k], Igens[i])
      images[i] = mR\J
    end
    for i = 1:length(IPgens)
      Pn = induce_image(Aut[k],IPgens[i])
      images[i+length(Igens)] = mR.disc_log_inf_plc[Pn]
    end

    if mp != id_hom(R)
      G[k] = hom(genstot, images, check = true)
    else
      G[k] = hom(GrpAbFinGenElem[mp(x) for x = genstot], GrpAbFinGenElem[mp(x) for x = images], check = true)
    end
    @hassert :RayFacElem 1 isbijective(G[k])
  end
  return G

end

function find_gens_for_action(mR::MapClassGrp)

	lp, sR = find_gens(pseudo_inv(mR), PrimesSet(2, -1))
  ip = InfPlc[]
  sR1 = GrpAbFinGenElem[]
	return lp, ip, sR, sR1

end

#  Find places that generate the ray class group
function find_gens_for_action(mR::MapRayClassGrp)
  O = order(codomain(mR))
  R = domain(mR)
  m = mR.defining_modulus[1]
  mm = minimum(m)
  lp = NfOrdIdl[]
  sR = GrpAbFinGenElem[]
  ip = InfPlc[]
  sR1 = GrpAbFinGenElem[]
  q, mq = quo(R, sR, false)

  #  First, generators of the multiplicative group.
  #  If the class group is trivial, they are enough
  if isdefined(mR, :fact_mod) && !isempty(mR.fact_mod)
    gens_and_logs = mR.gens_mult_grp_disc_log
    for j = 1:length(gens_and_logs)
    	if iszero(mq(gens_and_logs[j][2]))
    	  continue
    	end
    	push!(sR, gens_and_logs[j][2])
    	push!(lp, ideal(O, gens_and_logs[j][1]))
      q, mq = quo(R, sR, false)
      if order(q) == 1
        return lp, ip, sR, sR1
      end
    end
  end

  if !isempty(defining_modulus(mR)[2])
    dlog_dict = mR.disc_log_inf_plc
    for (p, v) in dlog_dict
      if iszero(mq(v))
        continue
      end
      push!(ip, p)
      push!(sR1, v)
      q, mq = quo(R, vcat(sR, sR1), false)
      if order(q) == 1
        return lp, ip, sR, sR1
      end
    end
  end
  #Now, gens of class group. Those are cached in the class group map
  mC = mR.clgrpmap
  if isdefined(mC, :small_gens)
    for P in mC.small_gens
      if iscoprime(minimum(P, copy = false), mm)
        push!(lp, P)
        push!(sR, mR\P)
      end
    end
  end
  q, mq = quo(R, vcat(sR, sR1), false)
  if order(q) != 1
    lpmC, srmC = find_gens(pseudo_inv(mC), PrimesSet(2, -1), mm)
    for P in lpmC
      push!(lp, P)
      push!(sR, mR\P)
    end
  end
  return lp, ip, sR, sR1
end

################################################################################
#
#  Generator 1 mod m
#
################################################################################

@doc Markdown.doc"""
    has_principal_generator_1_mod_m(I::NfOrdIdl, m::NfOrdIdl, inf_plc::Vector{InfPlc} = InfPlc[]) -> Bool, NfOrdElem

Given an ideal $I$, this function checks if the ideal is trivial in the ray class group mod ($m$, inf_plc).
If this is the case, we also return a generator which is 1 mod $m$. If not, the second return value is wrong.
"""
function has_principal_generator_1_mod_m(I::Union{NfOrdIdl, FacElem{NfOrdIdl, NfOrdIdlSet}}, m::NfOrdIdl, inf_plc::Vector{InfPlc} = InfPlc[]; GRH::Bool = true)

  # This function could be optimized if I cache some stuff from the construction
  # of the ray class group, but only in the case of the full ray_class_group
  # and not in the quotient.

  @assert all(x -> isreal(x), inf_plc)

  O = order(m)
  C, mC = class_group(O, GRH = GRH)
  fl, gen = isprincipal_fac_elem(I)
  if !fl
    return fl, gen
  end
  U, mU = unit_group_fac_elem(O, GRH = GRH)

  lp = factor(m)
  powers = Tuple{NfOrdIdl, NfOrdIdl}[(x, x^v) for (x, v) in lp]
  quo_rings = Tuple{NfOrdQuoRing, NfOrdQuoMap}[quo(O, q) for (x, q) in powers]
  groups_and_maps = Tuple{GrpAbFinGen, Hecke.GrpAbFinGenToAbsOrdQuoRingMultMap{NfAbsOrd{AnticNumberField,nf_elem},NfAbsOrdIdl{AnticNumberField,nf_elem},NfAbsOrdElem{AnticNumberField,nf_elem}}}[multiplicative_group(Q[1]) for Q in quo_rings]
  invariants = Vector{fmpz}()
  for x in groups_and_maps
    append!(invariants, x[1].snf)
  end
  for x in inf_plc
    push!(invariants, fmpz(2))
  end
  G = abelian_group(invariants)

  expo = exponent(G)
  tobeeval1 = FacElem{nf_elem, AnticNumberField}[mU(x) for x in gens(U)]
  push!(tobeeval1, gen)
  tobeeval = _preproc(m, tobeeval1, expo)
  coeffs = Vector{fmpz_mat}(undef, length(tobeeval))
  for i = 1:length(coeffs)
    coeffs[i] = zero_matrix(FlintZZ, 1, ngens(G))
  end
  ii = 1
  for i = 1:length(powers)
    P, Q = powers[i]
    exponq = gcd(expo, norm(Q)-divexact(norm(Q), norm(P)))
    els = fac_elems_eval(P, Q, tobeeval, exponq)
    for t = 1:length(els)
      el = els[t]
      y = (groups_and_maps[i][2]\quo_rings[i][1](el)).coeff
      for s = 1:ncols(y)
        coeffs[t][1, ii-1+s] = y[1, s]
      end
    end
    ii += ngens(groups_and_maps[i][1])
  end
  if !isempty(inf_plc)
    H, eH, lH = infinite_primes_map(O, inf_plc, m)
    for t = 1:length(tobeeval)
      el = lH(tobeeval[i])
      for s = 1:ngens(H)
        coeffs[t][1, ii-1+s] = el[s]
      end
    end
  end
  els_G = GrpAbFinGenElem[G(x) for x in coeffs]
  S, mS = sub(G, els_G[1:end-1])
  fl1, coord = haspreimage(mS, els_G[end])
  if !fl1
    return false, gen
  end
  @assert ngens(S) == ngens(U)
  for i = 1:ngens(U)
    if coord[i] != 0
      mul!(gen, gen, tobeeval1[i]^-Int(coord[i]))
    end
  end
  return true, gen
end

function disc_log_generalized_ray_class_grp(I::Union{NfOrdIdl, FacElem{NfOrdIdl, NfOrdIdlSet}}, mr::MapRayClassGrp)
  el = mr\I
  J = mr(el)
  I1 = I * inv(J)
  fl1, gen1 = has_principal_generator_1_mod_m(I1, mr.defining_modulus[1], mr.defining_modulus[2])
  @assert fl1
  return gen1, el
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

module test_withoutMacros

import ..ModiaMath
using Test
using LinearAlgebra
using StaticArrays

include(joinpath(ModiaMath.path, "examples", "withoutMacros", "Simulate_Pendulum.jl"))
import .Simulate_Pendulum
pendulumWithoutMacro = ModiaMath.SimulationModel(Simulate_Pendulum.PendulumWithoutMacro(L=0.5, m=1.0, d=0.1, phi0_deg=60.0))

include(joinpath(ModiaMath.path, "examples", "withoutMacros", "Simulate_FreeBodyRotation.jl"))
import .Simulate_FreeBodyRotation
freeBodyRotationWithoutMacro = ModiaMath.SimulationModel(Simulate_FreeBodyRotation.FreeBodyRotationWithoutMacro(), stopTime=5.0, tolerance=1e-8)

@testset "ModiaMath: withoutMacro/Simulate_Pendulum.jl" begin
    result = ModiaMath.simulate!(pendulumWithoutMacro, stopTime=10.0, tolerance=1e-8, interval=0.1, log=true)
    phi = result.series["phi"]
    w   = result.series["w"]

    @test isapprox(phi[end], 0.120289; atol=1e-3 )
    @test isapprox(w[end], 0.273965; atol=1e-2 )
end

@testset "ModiaMath: withoutMacro/Simulate_FreeBodyRotation.jl" begin
    result = ModiaMath.simulate!(freeBodyRotationWithoutMacro, log=true)
    ModiaMath.plot(result, [:q, :w, :derw, :tau])
end

end


text:
import LinearInterpolations as LI

using ArgCheck
using CoordinateTransformations
using ConstructionBase: setproperties
using StaticArrays
using AxisArrayConversion: to
using LinearAlgebra

export create_interpolate
export pullback, pullback!

"""
    struct ThreadsEx end

Executes the computation using multiple CPU threads.
"""
struct ThreadsEx end

function create_interpolate(data; kw...)
    obj = AC.to(NamedTuple, data)
    axes = Tuple(obj.axes)
    itp = LI.Interpolate(axes, obj.values; kw...)
    return itp
end

function pullback(f, axes, data; executor=ThreadsEx(), kw...)
    AC.check_consistency(data)
    data_nt    = AC.to(NamedTuple, data)
    axes_nt    = AC.name_axes(axes)
    out_nt = _pullback_namedtuple(f, axes_nt, data_nt, executor; kw...)
    T      = AC.roottype(typeof(data))
    out    = AC.to(T, out_nt)
    return out
end
function pullback!(f, out, data; executor=ThreadsEx(), kw...)
    AC.check_consistency(out)
    AC.check_consistency(data)
    out_nt  = AC.to(NamedTuple, out)
    data_nt = AC.to(NamedTuple, data)
    _pullback_namedtuple!(f, out_nt, data_nt, executor; kw...)
    T = AC.roottype(typeof(out))
    return AC.to(T, out_nt)
end

function _pullback_namedtuple(f, axes::NamedTuple, data::NamedTuple, ex::ThreadsEx; kw...)
    dims = map(length, Tuple(axes))
    T = float(eltype(data.values))
    vals = similar(data.values, T, dims)
    out  = (axes=axes, values=vals)
    return _pullback_namedtuple!(f, out, data, ex; kw...)
end

@noinline function _pullback_namedtuple!(f, out::NamedTuple, data::NamedTuple, ex::ThreadsEx;
                                         kw...)
    itp = create_interpolate(data; kw...)
    let out=out, f=f, itp=itp
        Threads.@threads for ci in CartesianIndices(out.values)
            pt0 = SVector(map(getindex, Tuple(out.axes), Tuple(ci)))
            pt = f(pt0)
            out.values[ci] = itp(pt)
        end
    end
    return out
end

function pullback_axes(f, axes)
    pullback_axes_tuple(f, Tuple(axes))
end

function pullback_axes_tuple(trans::Translation, axes::NTuple{N,Any}) where {N}
    @argcheck length(trans.translation) == length(axes)
    map(axes, NTuple{N}(trans.translation)) do ax, x
        ax .- x
    end
end

function norm0(arr; atol=0, rtol=sqrt(eps(float(eltype(arr)))))
    m = norm(arr, Inf)
    ret = 0
    thresh = max(atol, rtol*m)
    for x in arr
        if norm(x) >= thresh
            ret += 1
        end
    end
    ret
end

function pullback_axes_tuple(am::AffineMap, axes)
    axes = pullback_axes_tuple(Translation(am.translation), axes)
    axes = pullback_axes_tuple(LinearMap(am.linear), axes)
end

function pullback_axes_tuple(am::LinearMap, axes)
    sp = splitperm(am.linear)
    isp = inv_splitperm(sp)
    forward_axes_splitperm(isp, axes)
end

function splitperm(A::AbstractMatrix)
    @argcheck size(A,1) == size(A,2)
    N = size(A,1)
    @argcheck norm0(A) == N
    colmaxima = NTuple{N}(findmax(abs.(A), dims=(1,))[2])
    perm = map(colmaxima) do ci
        Tuple(ci)[1]
    end

    scalings = map(colmaxima) do ci
        A[ci]
    end
    (;perm, scalings)
end

function forward_axes_splitperm(sp, axes)
    axes = axes .* sp.scalings
    axes = getindex.(Ref(axes), sp.perm)
    map(axes) do ax
        first(ax) > last(ax) ? reverse(ax) : ax
    end
end

function matrix_from_splitperm(sp)
    N = length(sp.perm)
    T = eltype(sp.scalings)
    ret = zeros(T, N,N)
    for (i,j, s) in zip(1:N,sp.perm, sp.scalings)
        ret[j,i] = s
    end
    ret
end

function inv_splitperm(sp)
    perm = invperm(sp.perm)
    scalings = map(sp.scalings) do s
        1/s
    end
    scalings = getindex.(Ref(scalings), sp.perm)
    (;perm, scalings)
end

################################################################################
##### common_ground
################################################################################
function common_grid(grids)
    n = length(first(grids))
    axes_values = map(ntuple(identity, n)) do i
        axes_i = map(grids) do grid
            @argcheck length(grid) == n
            grid[i]
        end
        return common_axis(axes_i)
    end
    grid1 = first(grids)
    if grid1 isa Tuple
        axes_values
    else
        NamedTuple{propertynames(grid1)}(axes_values)
    end
end

function common_axis(axes)
    x_lo, x_hi = common_axis_bounds(axes)
    ret = filter(first(axes)) do x
        x_lo < x < x_hi
    end
    pushfirst!(ret, x_lo)
    push!(ret, x_hi)
    return ret
end

function common_axis_bounds(axs)
    x_lo = maximum(minimum, axs)
    x_hi = minimum(maximum, axs)
    return x_lo, x_hi
end

"""
    new_axis_arrays = common_ground(axis_arrays)

Take a collection of axis array like objects, with possibly distinct but overlapping
axes and shrink them such that their axes coincide.
"""
function common_ground(arrs)
    nts = map(arr -> AxisArrayConversion.to(NamedTuple, arr), arrs)
    axs = common_grid(map(nt -> nt.axes, nts))

    map(arrs, nts) do orig, nt
        nt_common = pullback(identity, axs, nt)
        T = AxisArrayConversion.roottype(typeof(orig))
        AxisArrayConversion.to(T, nt_common)
    end
end

function flip_decreasing_axes(nt::NamedTuple)
    flip_needs = map(nt.axes) do ax
        last(ax) < first(ax)
    end
    if any(flip_needs)
        axes_new = map(nt.axes, flip_needs) do ax, needs_flip
            needs_flip ? reverse(ax) : ax
        end
        inds_new = map(Base.axes(nt.values), flip_needs) do inds, needs_flip
            needs_flip ? reverse(inds) : inds
        end
        vals_new = nt.values[inds_new...]
        return (axes=axes_new, values=vals_new)
    else
        return nt
    end
end

function flip_decreasing_axes(arr)
    nt = AxisArrayConversion.to(NamedTuple, arr)
    nt_out = flip_decreasing_axes(nt)
    T = AxisArrayConversion.roottype(typeof(arr))
    AxisArrayConversion.to(T, nt_out)
end


text:
module FuNN

using Random

export

    # Classifier
    NeuralNetwork,

    # Main methods
    train!,
    predict,
    one_hot,
    one_cold,

    # Activation and loss functions
    relu,
    logistic,
    softmax,
    cross_entropy_binary,
    cross_entropy_multi,
    mean_squared_error

include("./NeuralNetwork.jl")
include("./back_propagation.jl")
include("./forward_propagation.jl")
include("./main_methods.jl")
include("./update_methods.jl")

end # module


text:
export iauPnm06a
"""
Form the matrix of precession-nutation for a given date (including
frame bias), IAU 2006 precession and IAU 2000A nutation models.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1,date2 double       TT as a 2-part Julian Date (Note 1)

Returned:
   rnpb        double[3][3] bias-precession-nutation matrix (Note 2)

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The matrix operates in the sense V(date) = rnpb * V(GCRS), where
   the p-vector V(date) is with respect to the true equatorial triad
   of date date1+date2 and the p-vector V(GCRS) is with respect to
   the Geocentric Celestial Reference System (IAU, 2000).

Called:
   iauPfw06     bias-precession F-W angles, IAU 2006
   iauNut06a    nutation, IAU 2006/2000A
   iauFw2m      F-W angles to r-matrix

Reference:

   Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855.

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPnm06a(double date1, double date2, double rnpb[3][3])

function iauPnm06a(date1::Real, date2::Real)
   rbpn = zeros(Float64, 3, 3)

   ccall((:iauPnm06a, libsofa_c), Cvoid,
       (Cdouble, Cdouble, Ptr{Cdouble}),
       convert(Float64, date1),
       convert(Float64, date2),
       rbpn)

   return SMatrix{3,3}(rbpn')
end


text:
__precompile__()

module DeepQLearning

push!(LOAD_PATH, "./")
push!(LOAD_PATH, "../../src/")

# using Simulation_Env_module_nnpolicy
using Simulation_Env_module

# package code goes here
using Distributions, StatsBase, Parameters
using TensorFlow
using POMDPs, POMDPToolbox, DeepRL
using JLD

const tf = TensorFlow

export DeepQLearningSolver,
       QNetworkArchitecture,
       DeepRecurrentQLearningSolver,
       RecurrentQNetworkArchitecture,
       # tf helpers
       flatten,
       dense,
       conv2d,
       mlp,
       cnn_to_mlp,
       get_train_vars_by_name,
       init_session,
       logg_scalar,

       # replay buffer
       DQExperience,
       ReplayBuffer,
       PrioritizedReplayBuffer,
       is_full,
       max_size,
       add_exp!,
       populate_replay_buffer!,
       update_priorities!,

       # training
       TrainGraph,
       DQNPolicy,
       LSTMPolicy,
       reset_hidden_state!,
       get_action,
       get_fused_q_vals,
       get_action!,
       eval_q,
       eval_lstm,
       action,
       solve,
       save,
       restore

include("tf_helpers.jl")

"""
    QNetworkArchitecture
Specify the Q network architecture as convolutional layers followed by a multi layer perceptron
- `fc::Vector{Int64}` the number of nodes and hidden layers [8, 32] => 2 layers of 8 and 32 nodes
- `conv::Vector{Tuple{Int64, Vector{Int64}, Int64}}` [(#filters, (kernel_size,kernel_size), stride)]
"""
@with_kw mutable struct QNetworkArchitecture
   fc::Vector{Int64} = Vector{Int64}[]
   conv::Vector{Tuple{Int64, Vector{Int64}, Int64}} = Vector{Tuple{Int64, Tuple{Int64, Int64}, Int64}}[]
end

"""
    DeepQLearningSolver

## Fields:
- `arch::QNetworkArchitecture` Specify the architecture of the Q network default = QNetworkArchitecture(conv=[], fc=[])
- `lr::Float64` learning rate default = 0.005
- `max_steps::Int64` total number of training step default = 1000
- `target_update_freq::Int64` frequency at which the target network is updated default = 500
- `batch_size::Int64` batch size sampled from the replay buffer default = 32
- `train_freq::Int64` frequency at which the active network is updated default  = 4
- `log_freq::Int64` frequency at which to logg info default = 100
- `eval_freq::Int64` frequency at which to eval the network default = 100
- `num_ep_eval::Int64` number of episodes to evaluate the policy default = 100
- `eps_fraction::Float64` fraction of the training set used to explore default = 0.5
- `eps_end::Float64` value of epsilon at the end of the exploration phase default = 0.01
- `double_q::Bool` double q learning udpate default = true
- `dueling::Bool` dueling structure for the q network default = true
- `prioritized_replay::Bool` enable prioritized experience replay default = true
- `prioritized_replay_alpha::Float64` default = 0.6
- `prioritized_replay_epsilon::Float64` default = 1e-6
- `prioritized_replay_beta::Float64` default = 0.4
- `buffer_size::Int64` size of the experience replay buffer default = 1000
- `max_episode_length::Int64` maximum length of a training episode default = 100
- `train_start::Int64` number of steps used to fill in the replay buffer initially default = 200
- `grad_clip::Bool` enables gradient clipping default = true
- `clip_val::Float64` maximum value for the grad norm default = 10.0
- `rng::AbstractRNG` random number generator default = MersenneTwister(0)
- `verbose::Bool` default = true
"""
@with_kw mutable struct DeepQLearningSolver
    arch::QNetworkArchitecture = QNetworkArchitecture(conv=[], fc=[])
    lr::Float64 = 0.005
    max_steps::Int64 = 1000
    target_update_freq::Int64 = 500
    batch_size::Int64 = 32
    train_freq::Int64  = 4
    log_freq::Int64 = 100
    eval_freq::Int64 = 100
    num_ep_eval::Int64 = 100
    eps_fraction::Float64 = 0.5
    eps_start::Float64 = 1.0
    eps_end::Float64 = 0.01
    double_q::Bool = true
    dueling::Bool = true
    prioritized_replay::Bool = true
    prioritized_replay_alpha::Float64 = 0.6
    prioritized_replay_epsilon::Float64 = 1e-6
    prioritized_replay_beta::Float64 = 0.4
    buffer_size::Int64 = 1000
    max_episode_length::Int64 = 100
    train_start::Int64 = 200
    grad_clip::Bool = true
    clip_val::Float64 = 10.0
    rng::AbstractRNG = MersenneTwister(0)
    if_log::Bool = true
    logdir::String = "log"
    save_freq::Int64 = 10000
    verbose::Bool = true
    weighted_correction_sum::Bool = false
    train::Bool = true
end

include("recurrent_deep_q_learning.jl")
include("policy.jl")
include("graph.jl")
include("experience_replay.jl")
include("episode_replay.jl")
include("prioritized_experience_replay.jl")
include("q_network.jl")
include("solver.jl")
include("saver.jl")

end # module

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
# run a Verlet tarjectory from a double well in 1D
using Plots
using Printf
using LinearAlgebra
using ForwardDiff
using QuadGK
using BasicMD

include("potentials.jl")

M = 1.0;
q₀ = [-0.5];
p₀ = [1.0];
Δt = 1e-1;
n_iters = 10^3;

V = x->DoubleWell(x);
cfg = ForwardDiff.GradientConfig(V, q₀);
gradV! = (gradV, x)-> ForwardDiff.gradient!(gradV, V, x, cfg);

sampler = Verlet(gradV!, M, Δt);

Q₀ = copy(q₀);
P₀ = copy(p₀);
sample_trajectory!([Q₀, P₀], sampler, options=MDOptions(n_iters=n_iters));
@printf("In Place (Q,P) after %d iterations: (%g,%g)\n",n_iters, Q₀[1], P₀[1]);

Xvals = sample_trajectory([q₀, p₀], sampler,
    options=MDOptions(n_iters=n_iters,n_save_iters=n_iters));
Q = Xvals[end][1][1]
P = Xvals[end][2][1]
@printf("(Q,P) after %d iterations: (%g,%g)\n",n_iters, Q[1], P[1]);

Xvals = sample_trajectory([q₀, p₀], sampler, options=MDOptions(n_iters=n_iters));
Q = Xvals[end][1][1]
P = Xvals[end][2][1]
@printf("(Q,P) after %d iterations: (%g,%g)\n",n_iters, Q[1], P[1]);
plot(Δt *(1:n_iters), [X[1][1] for X in Xvals],label="q")
plot!(Δt *(1:n_iters), [X[2][1] for X in Xvals],label="p")
xlabel!("t")

max_stars_repo_name: gideonsimpson/JuBasicMD


text:
  function setvalue!(value,thisseq,h,x,lossIDX,target,Tmax,W,H1)
      H=length(value[h[1]])
                 for tind=1:Tmax
               value[h[tind]]=ones(H,1)/H1
           end

            T=length(thisseq)
            for tind=1:Tmax-length(thisseq)
                value[x[tind]]=0
            end
            counter=0

            for tind=2:Tmax
                value[lossIDX[tind]]=0 # initialise all to zero
                value[target[tind]]=BitArray(W,1)
                #fill!(value[target[tind]],false)
            end

            counter=0
            for tind=Tmax-length(thisseq)+2:Tmax
                counter+=1
                fill!(value[target[tind]],false)
                value[target[tind]][thisseq[counter+1]]=true
                value[lossIDX[tind]]=1
                value[x[tind-1]]=thisseq[counter]
            end
      end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
# Local compute version

# @show ARGS
include("parsecommands.jl")

## Load all required packages
using Distributed

# check_procs(4) # make sure there are 4 processes waiting before loading packages

using Dates, Statistics
@everywhere using Caesar
@everywhere using JLD2
using CoordinateTransformations, Rotations, StaticArrays

using AprilTags
using Images, ImageDraw
# using ImageView

@everywhere begin
using Fontconfig
using Cairo
using Compose
using Gadfly
using RoMEPlotting
# using FileIO
# using GeometryTypes # using MeshCat
end

# setup configuration
using YAML

include(joinpath(dirname(@__FILE__),"configParameters.jl") )

# Figure export folder
global currdirtime = now()
if parsed_args["previous"] != ""
    # currdirtime = "2018-10-28T23:17:30.067"
    # currdirtime = "2018-11-03T22:48:51.924"
    # currdirtime = "2018-11-07T01:36:52.274"
    # currdirtime = "2018-11-08T21:56:16.991"
  # currdirtime = "2018-10-28T23:17:30.067"
  # currdirtime = "2018-11-03T22:48:51.924"
  # currdirtime = "2018-11-07T01:36:52.274"
  currdirtime = parsed_args["previous"]
end

resultsparentdir = joinpath(datadir, "results")
resultsdir = joinpath(resultsparentdir, "$(currdirtime)")

if parsed_args["previous"] == ""
  # When running fresh from new data
  include(joinpath(dirname(@__FILE__),"createResultsDirs.jl"))
end

# Utils required for this processing script
include(joinpath(dirname(@__FILE__),"racecarUtils.jl") )
include(joinpath(dirname(@__FILE__),"cameraUtils.jl") )
# include(joinpath(dirname(@__FILE__),"visualizationUtils.jl") )

@everywhere function plotRacecarInterm(fgl::G, resultsdirl, psyml::Symbol)::Nothing where G <: AbstractDFG
  @show ls(fgl)
  pl = drawPosesLandms(fgl, spscale=0.1, drawhist=false, meanmax=:mean) #,xmin=-3,xmax=6,ymin=-5,ymax=2);
  Gadfly.draw(PNG(joinpath(resultsdirl, "images", "$(psyml).png"),15cm, 10cm),pl)
  pl = drawPosesLandms(fgl, spscale=0.1, meanmax=:mean) # ,xmin=-3,xmax=3,ymin=-2,ymax=2);
  Gadfly.draw(PNG(joinpath(resultsdirl, "images", "hist_$(psyml).png"),15cm, 10cm),pl)
  pl = plotPose2Vels(fgl, Symbol("$(psyml)"), coord=Coord.Cartesian(xmin=-1.0, xmax=1.0))
  Gadfly.draw(PNG(joinpath(resultsdirl, "images", "vels_$(psyml).png"),15cm, 10cm),pl)
  # save combined image with tags
  nothing
end

global tag_bag = Dict()
if parsed_args["previous"] == ""
  tag_bag = detectTagsInImgs(datafolder, imgfolder, resultsdir, camidxs)
  # save the tag bag file for future use
  @save resultsdir*"/tag_det_per_pose.jld2" tag_bag
else
  @load resultsdir*"/tag_det_per_pose.jld2" tag_bag
end

fg = main(resultsdir, camidxs, tag_bag, jldfile=parsed_args["jldfile"], failsafe=parsed_args["failsafe"], show=parsed_args["show"]  )

0

# fails
# key 1 not found
# julia101 -p 4 apriltag_and_zed_slam.jl --previous "2018-11-09T01:42:33.279" --jldfile "racecar_fg_x299.jld2" --folder_name "labrun7" --failsafe

#


text:
type LMMCurveState <: CurveState
  numberOfRates::Int
  rateTimes::Vector{Float64}
  rateTaus::Vector{Float64}
  firstIdx::Int
  discRatios::Vector{Float64}
  forwardRates::Vector{Float64}
  cmSwapRates::Vector{Float64}
  cmSwapAnnuties::Vector{Float64}
  cotSwapRates::Vector{Float64}
  cotAnnuities::Vector{Float64}
  firstCotAnnuityComped::Int
end

function LMMCurveState(rateTimes::Vector{Float64})
  numberOfRates = isempty(rateTimes) ? 0 : length(rateTimes) - 1
  rateTaus = check_increasing_times_and_calculate_taus(rateTimes)

  discRatios = ones(numberOfRates + 1)
  forwardRates = Vector{Float64}(numberOfRates)
  cmSwapRates = Vector{Float64}(numberOfRates)
  cmSwapAnnuties = fill(rateTaus[numberOfRates], numberOfRates)
  cotSwapRates = Vector{Float64}(numberOfRates)
  cotAnnuities = fill(rateTaus[numberOfRates], numberOfRates)

  return LMMCurveState(numberOfRates, rateTimes, rateTaus, numberOfRates+1, discRatios, forwardRates, cmSwapRates, cmSwapAnnuties, cotSwapRates,
                      cotAnnuities, numberOfRates+1)
end

function set_on_forward_rates!(lmm::LMMCurveState, rates::Vector{Float64}, firstValidIndex::Int=1)
  length(rates) == lmm.numberOfRates || error("rates mismatch")
  firstValidIndex <= lmm.numberOfRates || error("first valid index must be less than or equal to numberOfRates")

  # first copy input
  lmm.firstIdx = firstValidIndex
  lmm.forwardRates[lmm.firstIdx:end] = rates[lmm.firstIdx:end]

  # then calculate the discount ratios
  # taken care at constructor time discRatios[numberOfRates] = 1.0
  for i = lmm.firstIdx:lmm.numberOfRates
    lmm.discRatios[i+1] = lmm.discRatios[i] / (1.0 + lmm.forwardRates[i] * lmm.rateTaus[i])
  end

  # lazy evaluation of :
  # - coterminal swap rates/annuities
  # - constant maturity swap rates/annuities

  lmm.firstCotAnnuityComped = lmm.numberOfRates+1
end

function discount_ratio(lmm::LMMCurveState, i::Int, j::Int)
  lmm.firstIdx <= lmm.numberOfRates || error("curve state not initialized yet")
  min(i, j) >= lmm.firstIdx || error("invalid index")
  max(i, j) <= lmm.numberOfRates + 1 || error("invalid index")

  return lmm.discRatios[i] / lmm.discRatios[j]
end

function coterminal_from_discount_ratios!(lmm::LMMCurveState, firstValidIndex::Int, discountFactors::Vector{Float64}, taus::Vector{Float64},
                                        cotSwapRates::Vector{Float64}, cotSwapAnnuities::Vector{Float64})
  nCotSwapRates = length(cotSwapRates)
  length(taus) == nCotSwapRates || error("taus size != cotSwapRate size")
  length(cotSwapAnnuities) == nCotSwapRates || error("cotSwapAnnuities size != cotSwapRate size")
  length(discountFactors) == nCotSwapRates + 1 || error("discountFactors size != cotSwapRate size + 1")

  cotSwapAnnuities[nCotSwapRates] = taus[nCotSwapRates] * discountFactors[nCotSwapRates + 1]
  cotSwapRates[nCotSwapRates] = (discountFactors[nCotSwapRates] - discountFactors[nCotSwapRates + 1]) / cotSwapAnnuities[nCotSwapRates]

  for i = nCotSwapRates:-1:firstValidIndex + 1
    cotSwapAnnuities[i-1] = cotSwapAnnuities[i] + taus[i-1] * discountFactors[i]
    cotSwapRates[i-1] = (discountFactors[i-1] - discountFactors[nCotSwapRates+1]) / cotSwapAnnuities[i-1]
  end

  return cotSwapAnnuities, cotSwapRates
end

function forward_rate(lmm::LMMCurveState, i::Int)
  lmm.firstIdx <= lmm.numberOfRates || error("curve state not initialized yet")
  (i >= lmm.firstIdx && i <= lmm.numberOfRates) || error("invalid index")

  return lmm.forwardRates[i]
end

function coterminal_swap_annuity(lmm::LMMCurveState, numeraire::Int, i::Int)
  lmm.firstIdx <= lmm.numberOfRates || error("curve state not initialized yet")
  (numeraire >= lmm.firstIdx && numeraire <= lmm.numberOfRates + 1) || error("invalid numeraire")
  (i >= lmm.firstIdx && i <= lmm.numberOfRates) || error("invalid index")

  if lmm.firstCotAnnuityComped <= i
    return lmm.cotAnnuities[i] / lmm.discRatios[numeraire]
  end

  if lmm.firstCotAnnuityComped == lmm.numberOfRates+1
    lmm.cotAnnuities[lmm.numberOfRates] = lmm.rateTaus[lmm.numberOfRates] * lmm.discRatios[lmm.numberOfRates + 1]
    lmm.firstCotAnnuityComped -= 1
  end

  for j = lmm.firstCotAnnuityComped-1:-1:i
    lmm.cotAnnuities[j] = lmm.cotAnnuities[j+1] + lmm.rateTaus[j] * lmm.discRatios[j+1]
  end

  lmm.firstCotAnnuityComped = i

  return lmm.cotAnnuities[i] / lmm.discRatios[numeraire]
end

function coterminal_swap_rate(lmm::LMMCurveState, i::Int)
  lmm.firstIdx <= lmm.numberOfRates || error("curve state not initialized yet")
  (i >= lmm.firstIdx && i <= lmm.numberOfRates) || error("invalid index")

  res = (lmm.discRatios[i] / lmm.discRatios[lmm.numberOfRates+1] - 1.0) / coterminal_swap_annuity(lmm, lmm.numberOfRates+1, i)

  return res
end

function coterminal_swap_rates!(lmm::LMMCurveState)
  lmm.firstIdx <= lmm.numberOfRates || error("curve state not initialized yet")
  coterminal_from_discount_ratios!(lmm, lmm.firstIdx, lmm.discRatios, lmm.rateTaus, lmm.cotSwapRates, lmm.cotAnnuities)

  return lmm.cotSwapRates
end

clone(lmm::LMMCurveState) = LMMCurveState(lmm.numberOfRates, copy(lmm.rateTimes), copy(lmm.rateTaus), lmm.firstIdx, copy(lmm.discRatios), copy(lmm.forwardRates),
                                          copy(lmm.cmSwapRates), copy(lmm.cmSwapAnnuties), copy(lmm.cotSwapRates), copy(lmm.cotAnnuities), lmm.firstCotAnnuityComped)