text:
"""
    Barrier(comm::Comm)

Blocks until `comm` is synchronized.

If `comm` is an intracommunicator, then it blocks until all members of the group have called it.

If `comm` is an intercommunicator, then it blocks until all members of the other group have called it.

# External links
$(_doc_external("MPI_Barrier"))
"""
function Barrier(comm::Comm)
    # int MPI_Barrier(MPI_Comm comm)
    @mpichk ccall((:MPI_Barrier, libmpi), Cint, (MPI_Comm,), comm)
    return nothing
end

"""
    Ibarrier(comm::Comm)

Blocks until `comm` is synchronized.

If `comm` is an intracommunicator, then it blocks until all members of the group have called it.

If `comm` is an intercommunicator, then it blocks until all members of the other group have called it.

# External links
$(_doc_external("MPI_Ibarrier"))
"""
function Ibarrier(comm::Comm)
    req = Request()
    # int MPI_Ibarrier(MPI_Comm comm, MPI_Req req)
    @mpichk ccall((:MPI_Ibarrier, libmpi), Cint, (MPI_Comm, Ptr{MPI_Request}), comm, req)
    return req
end

"""
    Bcast!(buf, comm::Comm; root::Integer=0)

Broadcast the buffer `buf` from `root` to all processes in `comm`.

# See also
- [`bcast`](@ref)

# External links
$(_doc_external("MPI_Bcast"))
"""
Bcast!(buf, comm::Comm; root::Integer=Cint(0)) =
    Bcast!(buf, root, comm)

function Bcast!(buf::Buffer, root::Integer, comm::Comm)
    # int MPI_Bcast(void* buffer, int count, MPI_Datatype datatype, int root,
    #               MPI_Comm comm)
    @mpichk ccall((:MPI_Bcast, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, Cint, MPI_Comm),
                  buf.data, buf.count, buf.datatype, root, comm)
    return buf.data
end
function Bcast!(data, root::Integer, comm::Comm)
    Bcast!(Buffer(data), root, comm)
end

"""
    bcast(obj, comm::Comm; root::Integer=0)

Broadcast the object `obj` from rank `root` to all processes on `comm`. This is
able to handle arbitrary data.

# See also

- [`Bcast!`](@ref)
"""
bcast(obj, comm::Comm; root::Integer=Cint(0)) =
    bcast(obj, root, comm)
function bcast(obj, root::Integer, comm::Comm)
    isroot = Comm_rank(comm) == root
    count = Ref{Cint}()
    if isroot
        buf = MPI.serialize(obj)
        count[] = length(buf)
    end
    Bcast!(count, root, comm)
    if !isroot
        buf = Array{UInt8}(undef, count[])
    end
    Bcast!(buf, root, comm)
    if !isroot
        obj = MPI.deserialize(buf)
    end
    return obj
end

"""
    Scatter!(sendbuf::Union{UBuffer,Nothing}, recvbuf, comm::Comm;
        root::Integer=0)

Splits the buffer `sendbuf` in the `root` process into `Comm_size(comm)` chunks,
sending the `j`-th chunk to the process of rank `j-1` into the `recvbuf` buffer.

`sendbuf` on the root process should be a [`UBuffer`](@ref) (an `Array` can also be
passed directly if the sizes can be determined from `recvbuf`). On non-root processes it
is ignored, and `nothing` can be passed instead.

`recvbuf` is a [`Buffer`](@ref) object, or any object for which `Buffer(recvbuf)` is
defined. On the root process, it can also be [`MPI.IN_PLACE`](@ref), in which case it is
unmodified. For example:
```
if root == MPI.Comm_rank(comm)
    MPI.Scatter!(UBuffer(buf, count), MPI.IN_PLACE, comm; root=root)
else
    MPI.Scatter!(nothing, buf, comm; root=root)
end
```

# See also
- [`Scatterv!`](@ref) if the number of elements varies between processes.

# External links
$(_doc_external("MPI_Scatter"))
"""
Scatter!(sendbuf, recvbuf, comm::Comm; root::Integer=Cint(0)) =
    Scatter!(sendbuf, recvbuf, root, comm)
function Scatter!(sendbuf::UBuffer, recvbuf::Buffer, root::Integer, comm::Comm)
    if sendbuf.nchunks !== nothing && Comm_rank(comm) == root
        @assert sendbuf.nchunks >= Comm_size(comm)
    end
    # int MPI_Scatter(const void* sendbuf, int sendcount, MPI_Datatype sendtype,
    #                 void* recvbuf, int recvcount, MPI_Datatype recvtype, int root,
    #                 MPI_Comm comm)
    @mpichk ccall((:MPI_Scatter, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype,
                   MPIPtr, Cint, MPI_Datatype, Cint, MPI_Comm),
                  sendbuf.data, sendbuf.count, sendbuf.datatype,
                  recvbuf.data, recvbuf.count, recvbuf.datatype, root, comm)
    return recvbuf.data
end
Scatter!(sendbuf::UBuffer, recvbuf, root::Integer, comm::Comm) =
    Scatter!(sendbuf, Buffer(recvbuf), root, comm)
Scatter!(sendbuf::Nothing, recvbuf, root::Integer, comm::Comm) =
    Scatter!(UBuffer(nothing), recvbuf, root, comm)

# determine UBuffer count from recvbuf
Scatter!(sendbuf::AbstractArray{T}, recvbuf::Union{Ref{T},AbstractArray{T}}, root::Integer, comm::Comm) where {T} =
    Scatter!(UBuffer(sendbuf,length(recvbuf)), recvbuf, root, comm)

"""
    Scatter(sendbuf, T, comm::Comm; root::Integer=0)

Splits the buffer `sendbuf` in the `root` process into `Comm_size(comm)` chunks,
sending the `j`-th chunk to the process of rank `j-1` as an object of type `T`.

# See also
- [`Scatter!`](@ref)
"""
Scatter(sendbuf, T, comm; root::Integer=Cint(0)) =
    Scatter(sendbuf, T, root, comm)
Scatter(sendbuf, ::Type{T}, root::Integer, comm::Comm) where {T} =
    Scatter!(sendbuf, Ref{T}(), root, comm)[]

"""
    Scatterv!(sendbuf, recvbuf, comm::Comm; root::Integer=0)

Splits the buffer `sendbuf` in the `root` process into `Comm_size(comm)` chunks and sends
the `j`th chunk to the process of rank `j-1` into the `recvbuf` buffer.

`sendbuf` on the root process should be a [`VBuffer`](@ref). On non-root processes it is
ignored, and `nothing` can be passed instead.

`recvbuf` is a [`Buffer`](@ref) object, or any object for which `Buffer(recvbuf)` is
defined. On the root process, it can also be [`MPI.IN_PLACE`](@ref), in which case it is
unmodified. For example:
```
if root == MPI.Comm_rank(comm)
    MPI.Scatterv!(VBuffer(buf, counts), MPI.IN_PLACE, comm; root=root)
else
    MPI.Scatterv!(nothing, buf, comm; root=root)
end
```

# See also
- [`Scatter!`](@ref) if the number of elements are the same for all processes

# External links
$(_doc_external("MPI_Scatterv"))
"""
Scatterv!(sendbuf, recvbuf, comm::Comm; root::Integer=Cint(0)) =
    Scatterv!(sendbuf, recvbuf, root, comm)
function Scatterv!(sendbuf::VBuffer, recvbuf::Buffer, root::Integer, comm::Comm)
    if Comm_rank(comm) == root
        @assert length(sendbuf.counts) >= Comm_size(comm)
    end
    # int MPI_Scatterv(const void* sendbuf, const int sendcounts[],
    #                  const int displs[], MPI_Datatype sendtype, void* recvbuf,
    #                  int recvcount, MPI_Datatype recvtype, int root, MPI_Comm comm)
    @mpichk ccall((:MPI_Scatterv, libmpi), Cint,
                  (MPIPtr, Ptr{Cint}, Ptr{Cint},  MPI_Datatype,
                   MPIPtr, Cint, MPI_Datatype, Cint, MPI_Comm),
                  sendbuf.data, sendbuf.counts, sendbuf.displs, sendbuf.datatype,
                  recvbuf.data, recvbuf.count, recvbuf.datatype,
                  root, comm)
    return recvbuf.data
end
Scatterv!(sendbuf::VBuffer, recvbuf, root::Integer, comm::Comm) =
    Scatterv!(sendbuf, Buffer(recvbuf), root, comm)
Scatterv!(sendbuf::Nothing, recvbuf, root::Integer, comm::Comm) =
    Scatterv!(VBuffer(nothing), recvbuf, root, comm)

"""
    Gather!(sendbuf, recvbuf, comm::Comm; root::Integer=0)

Each process sends the contents of the buffer `sendbuf` to the `root` process. The `root`
process stores elements in rank order in the buffer buffer `recvbuf`.

`sendbuf` should be a [`Buffer`](@ref) object, or any object for which
[`Buffer_send`](@ref) is defined, with the same length on all processes, and should be the
same length on all processes.

On the root process, `sendbuf` can be [`MPI.IN_PLACE`](@ref) on the root process, in which
case the corresponding entries in `recvbuf` are assumed to be already in place (this
corresponds the behaviour of `MPI_IN_PLACE` in `MPI_Gather`). For example:
```
if root == MPI.Comm_rank(comm)
    MPI.Gather!(MPI.IN_PLACE, UBuffer(buf, count), comm; root=root)
else
    MPI.Gather!(buf, nothing, comm; root=root)
end
```

`recvbuf` on the root process should be a [`UBuffer`](@ref), or can be an `AbstractArray`
if the length can be determined from `sendbuf`. On non-root processes it is ignored and
can be `nothing`.

# See also
- [`Gather`](@ref) for the allocating operation.
- [`Gatherv!`](@ref) if the number of elements varies between processes.
- [`Allgather!`](@ref) to send the result to all processes.

# External links
$(_doc_external("MPI_Gather"))
"""
Gather!(sendbuf, recvbuf, comm::Comm; root::Integer=Cint(0)) =
    Gather!(sendbuf, recvbuf, root, comm)
function Gather!(sendbuf::Buffer, recvbuf::UBuffer, root::Integer, comm::Comm)
    if recvbuf.nchunks !== nothing && Comm_rank(comm) == root
        @assert recvbuf.nchunks >= Comm_size(comm)
    end
    # int MPI_Gather(const void* sendbuf, int sendcount, MPI_Datatype sendtype,
    #                void* recvbuf, int recvcount, MPI_Datatype recvtype, int root,
    #                MPI_Comm comm)
    @mpichk ccall((:MPI_Gather, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, MPIPtr, Cint, MPI_Datatype, Cint, MPI_Comm),
                  sendbuf.data, sendbuf.count, sendbuf.datatype,
                  recvbuf.data, recvbuf.count, recvbuf.datatype, root, comm)
    return recvbuf.data
end
Gather!(sendbuf, recvbuf::UBuffer, root::Integer, comm::Comm) =
    Gather!(Buffer_send(sendbuf), recvbuf, root, comm)
Gather!(sendbuf::Union{Ref,AbstractArray}, recvbuf::AbstractArray, root::Integer, comm::Comm) =
    Gather!(sendbuf, UBuffer(recvbuf, length(sendbuf)), root, comm)
Gather!(sendbuf, recvbuf::Nothing, root::Integer, comm::Comm) =
    Gather!(sendbuf, UBuffer(nothing), root, comm)
Gather!(sendbuf::Nothing, recvbuf, root::Integer, comm::Comm) =
    Gather!(IN_PLACE, recvbuf, root, comm)

"""
    Gather(sendbuf, comm::Comm; root=0)

Each process sends the contents of the buffer `sendbuf` to the `root` process. The `root`
allocates the output buffer and stores elements in rank order.

`sendbuf` can be an `AbstractArray` or a scalar, and should be the same length on all
processes.

# See also
- [`Gather!`](@ref) for the mutating operation.
- [`Gatherv!`](@ref) if the number of elements varies between processes.
- [`Allgather!`](@ref)/[`Allgather`](@ref) to send the result to all processes.

# External links
$(_doc_external("MPI_Gather"))
"""
Gather(sendbuf, comm::Comm; root::Integer=Cint(0)) =
    Gather(sendbuf, root, comm)
Gather(sendbuf::AbstractArray, root::Integer, comm::Comm) =
    Gather!(sendbuf, Comm_rank(comm) == root ? similar(sendbuf, Comm_size(comm) * length(sendbuf)) : nothing, root, comm)
Gather(object::T, root::Integer, comm::Comm) where {T} =
    Gather!(Ref(object), Comm_rank(comm) == root ? Array{T}(undef, Comm_size(comm)) : nothing, root, comm)

"""
    Gatherv!(sendbuf, recvbuf, comm::Comm; root::Integer=0)

Each process sends the contents of the buffer `sendbuf` to the `root` process. The `root`
stores elements in rank order in the buffer `recvbuf`.

`sendbuf` should be a [`Buffer`](@ref) object, or any object for which
[`Buffer_send`](@ref) is defined, with the same length on all processes.

On the root process, `sendbuf` can be [`MPI.IN_PLACE`](@ref), in which case the
corresponding entries in `recvbuf` are assumed to be already in place. For example
```
if root == MPI.Comm_rank(comm)
    Gatherv!(MPI.IN_PLACE, VBuffer(buf, counts), comm; root=root)
else
    Gatherv!(buf, nothing, comm; root=root)
end
```

`recvbuf` on the root process should be a [`VBuffer`](@ref), or can be an `AbstractArray`
if the length can be determined from `sendbuf`. On non-root processes it is ignored and
can be `nothing`.

# See also
- [`Gatherv!`](@ref) for the allocating operation
- [`Gather!`](@ref)
- [`Allgatherv!`](@ref) to send the result to all processes

# External links
$(_doc_external("MPI_Gatherv"))
"""
Gatherv!(sendbuf, recvbuf, comm::Comm; root::Integer=Cint(0)) =
    Gatherv!(sendbuf, recvbuf, root, comm)
function Gatherv!(sendbuf::Buffer, recvbuf::VBuffer, root::Integer, comm::Comm)
    if Comm_rank(comm) == root
        @assert length(recvbuf.counts) >= Comm_size(comm)
    end
    # int MPI_Gatherv(const void* sendbuf, int sendcount, MPI_Datatype sendtype,
    #                 void* recvbuf, const int recvcounts[], const int displs[],
    #                 MPI_Datatype recvtype, int root, MPI_Comm comm)
    @mpichk ccall((:MPI_Gatherv, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, MPIPtr, Ptr{Cint}, Ptr{Cint}, MPI_Datatype, Cint, MPI_Comm),
                  sendbuf.data, sendbuf.count, sendbuf.datatype,
                  recvbuf.data, recvbuf.counts, recvbuf.displs, recvbuf.datatype, root, comm)
    return recvbuf.data
end
Gatherv!(sendbuf, recvbuf::VBuffer, root::Integer, comm::Comm) =
    Gatherv!(Buffer_send(sendbuf), recvbuf, root, comm)
Gatherv!(sendbuf, recvbuf::Nothing, root::Integer, comm::Comm) =
    Gatherv!(sendbuf, VBuffer(nothing), root, comm)

"""
    Allgather!(sendbuf, recvbuf::UBuffer, comm::Comm)
    Allgather!(sendrecvbuf::UBuffer, comm::Comm)

Each process sends the contents of `sendbuf` to the other processes, the result of which
is stored in rank order into `recvbuf`.

`sendbuf` can be a [`Buffer`](@ref) object, or any object for which [`Buffer_send`](@ref)
is defined, and should be the same length on all processes.

`recvbuf` can be a [`UBuffer`](@ref), or can be an `AbstractArray` if the length can be
determined from `sendbuf`.

If only one buffer `sendrecvbuf` is provided, then on each process the data to send is
assumed to be in the area where it would receive its own contribution.

# See also
- [`Allgather`](@ref) for the allocating operation
- [`Allgatherv!`](@ref) if the number of elements varies between processes.
- [`Gather!`](@ref) to send only to a single root process

# External links
$(_doc_external("MPI_Allgather"))
"""
function Allgather!(sendbuf::Buffer, recvbuf::UBuffer, comm::Comm)
    if recvbuf.nchunks !== nothing
        @assert recvbuf.nchunks >= Comm_size(comm)
    end

    # int MPI_Allgather(const void* sendbuf, int sendcount,
    #                   MPI_Datatype sendtype, void* recvbuf, int recvcount,
    #                   MPI_Datatype recvtype, MPI_Comm comm)
    @mpichk ccall((:MPI_Allgather, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, MPIPtr, Cint, MPI_Datatype, MPI_Comm),
                  sendbuf.data, sendbuf.count, sendbuf.datatype,
                  recvbuf.data, recvbuf.count, recvbuf.datatype, comm)
    return recvbuf.data
end
Allgather!(sendbuf, recvbuf::UBuffer, comm::Comm) =
    Allgather!(Buffer_send(sendbuf), recvbuf, comm)

Allgather!(sendbuf::Union{Ref,AbstractArray}, recvbuf::AbstractArray, comm::Comm) =
    Allgather!(sendbuf, UBuffer(recvbuf, length(sendbuf)), comm)

function Allgather!(sendrecvbuf::UBuffer, comm::Comm)
    Allgather!(IN_PLACE, sendrecvbuf, comm)
end

"""
    Allgather(sendbuf, comm)

Each process sends the contents of `sendbuf` to the other processes, who store the results
in rank order allocating the output buffer.

`sendbuf` can be an `AbstractArray` or a scalar, and should be the same size on all
processes.

# See also
- [`Allgather!`](@ref) for the mutating operation
- [`Allgatherv!`](@ref) if the number of elements varies between processes.
- [`Gather!`](@ref) to send only to a single root process

# External links
$(_doc_external("MPI_Allgather"))
"""
function Allgather(sendbuf::AbstractArray, comm::Comm)
    Allgather!(sendbuf, similar(sendbuf, Comm_size(comm) * length(sendbuf)), comm)
end
function Allgather(object::T, comm::Comm) where T
    Allgather!(Ref(object), Array{T}(undef, Comm_size(comm)), comm)
end

"""
    Allgatherv!(sendbuf, recvbuf::VBuffer, comm::Comm)
    Allgatherv!(sendrecvbuf::VBuffer, comm::Comm)

Each process sends the contents of `sendbuf` to all other process. Each process stores the
received in the [`VBuffer`](@ref) `recvbuf`.

`sendbuf` can be a [`Buffer`](@ref) object, or any object for which [`Buffer_send`](@ref)
is defined.

If only one buffer `sendrecvbuf` is provided, then for each process, the data to be sent
is taken from the interval of `recvbuf` where it would store its own data.

# See also
- [`Gatherv!`](@ref) to send the result to a single process

# External links
$(_doc_external("MPI_Allgatherv"))
"""
function Allgatherv!(sendbuf::Buffer, recvbuf::VBuffer, comm::Comm)
    @assert length(recvbuf.counts) >= Comm_size(comm)

    # int MPI_Allgatherv(const void* sendbuf, int sendcount,
    #                    MPI_Datatype sendtype, void* recvbuf, const int recvcounts[],
    #                    const int displs[], MPI_Datatype recvtype, MPI_Comm comm)
    @mpichk ccall((:MPI_Allgatherv, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, MPIPtr, Ptr{Cint}, Ptr{Cint}, MPI_Datatype, MPI_Comm),
                  sendbuf.data, sendbuf.count, sendbuf.datatype,
                  recvbuf.data, recvbuf.counts, recvbuf.displs, recvbuf.datatype,
                  comm)
    return recvbuf.data
end
function Allgatherv!(sendbuf, recvbuf::VBuffer, comm::Comm)
    Allgatherv!(Buffer_send(sendbuf), recvbuf, comm)
end
function Allgatherv!(sendrecvbuf::VBuffer, comm::Comm)
    Allgatherv!(IN_PLACE, sendrecvbuf, comm)
end

"""
    Alltoall!(sendbuf::UBuffer, recvbuf::UBuffer, comm::Comm)
    Alltoall!(sendrecvbuf::UBuffer, comm::Comm)

Every process divides the [`UBuffer`](@ref) `sendbuf` into `Comm_size(comm)` chunks of
equal size, sending the `j`-th chunk to the process of rank `j-1`.  Every process stores
the data received from rank `j-1` process in the `j`-th chunk of the buffer `recvbuf`.

```
rank    send buf                        recv buf
----    --------                        --------
 0      a,b,c,d,e,f       Alltoall      a,b,A,B,α,β
 1      A,B,C,D,E,F  ---------------->  c,d,C,D,γ,ψ
 2      α,β,γ,ψ,η,ν                     e,f,E,F,η,ν
```

If only one buffer `sendrecvbuf` is used, then data is overwritten.

# See also
- [`Alltoall`](@ref) for the allocating operation

# External links
$(_doc_external("MPI_Alltoall"))
"""
function Alltoall!(sendbuf::UBuffer, recvbuf::UBuffer, comm::Comm)
    if sendbuf.data !== Consts.MPI_IN_PLACE[] && sendbuf.nchunks !== nothing
        @assert sendbuf.nchunks >= Comm_size(comm)
    end
    if recvbuf.nchunks !== nothing
        @assert recvbuf.nchunks >= Comm_size(comm)
    end
    # int MPI_Alltoall(const void* sendbuf, int sendcount, MPI_Datatype sendtype,
    #                  void* recvbuf, int recvcount, MPI_Datatype recvtype,
    #                  MPI_Comm comm)
    @mpichk ccall((:MPI_Alltoall, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, MPIPtr, Cint, MPI_Datatype, MPI_Comm),
                  sendbuf.data, sendbuf.count, sendbuf.datatype,
                  recvbuf.data, recvbuf.count, recvbuf.datatype,
                  comm)
    return recvbuf.data
end
Alltoall!(sendbuf::InPlace, recvbuf::UBuffer, comm::Comm) =
    Alltoall!(UBuffer(IN_PLACE), recvbuf, comm)
Alltoall!(sendrecvbuf::UBuffer, comm::Comm) =
    Alltoall!(IN_PLACE, sendrecvbuf, comm)

"""
    Alltoall(sendbuf::UBuffer, comm::Comm)

Every process divides the [`UBuffer`](@ref) `sendbuf` into `Comm_size(comm)` chunks of
equal size, sending the `j`-th chunk to the process of rank `j-1`. Every process allocates
the output buffer and stores the data received from the process on rank `j-1` in the
`j`-th chunk.

```
rank    send buf                        recv buf
----    --------                        --------
 0      a,b,c,d,e,f       Alltoall      a,b,A,B,α,β
 1      A,B,C,D,E,F  ---------------->  c,d,C,D,γ,ψ
 2      α,β,γ,ψ,η,ν                     e,f,E,F,η,ν
```

# See also
- [`Alltoall!`](@ref) for the mutating operation

# External links
$(_doc_external("MPI_Alltoall"))
"""
Alltoall(sendbuf::UBuffer,  comm::Comm) =
    Alltoall!(sendbuf, similar(sendbuf), comm)

"""
    Alltoallv!(sendbuf::VBuffer, recvbuf::VBuffer, comm::Comm)

Similar to [`Alltoall!`](@ref), except with different size chunks per process.

# See also
- [`VBuffer`](@ref)

# External links
$(_doc_external("MPI_Alltoallv"))
"""
function Alltoallv!(sendbuf::VBuffer, recvbuf::VBuffer, comm::Comm)
    if sendbuf.data !== Consts.MPI_IN_PLACE[]
        @assert length(sendbuf.counts) >= Comm_size(comm)
    end
    @assert length(recvbuf.counts) >= Comm_size(comm)
    # int MPI_Alltoallv(const void* sendbuf, const int sendcounts[],
    #                   const int sdispls[], MPI_Datatype sendtype, void* recvbuf,
    #                   const int recvcounts[], const int rdispls[],
    #                   MPI_Datatype recvtype, MPI_Comm comm)
    @mpichk ccall((:MPI_Alltoallv, libmpi), Cint,
                  (MPIPtr, Ptr{Cint}, Ptr{Cint}, MPI_Datatype,
                   MPIPtr, Ptr{Cint}, Ptr{Cint}, MPI_Datatype,
                   MPI_Comm),
                  sendbuf.data, sendbuf.counts, sendbuf.displs, sendbuf.datatype,
                  recvbuf.data, recvbuf.counts, recvbuf.displs, recvbuf.datatype,
                  comm)

    return recvbuf.data
end

### Reduce/Scan

## Reduce

# mutating
"""
    Reduce!(sendbuf, recvbuf, op, comm::Comm; root::Integer=0)
    Reduce!(sendrecvbuf, op, comm::Comm; root::Integer=0)

Performs elementwise reduction using the operator `op` on the buffer `sendbuf` and stores
the result in `recvbuf` on the process of rank `root`.

On non-root processes `recvbuf` is ignored, and can be `nothing`.

To perform the reduction in place, provide a single buffer `sendrecvbuf`.

# See also
- [`Reduce`](@ref) to handle allocation of the output buffer.
- [`Allreduce!`](@ref)/[`Allreduce`](@ref) to send reduction to all ranks.
- [`Op`](@ref) for details on reduction operators.

# External links
$(_doc_external("MPI_Reduce"))
"""
Reduce!(sendrecvbuf, op, comm::Comm; root::Integer=Cint(0)) =
    Reduce!(sendrecvbuf, op, root, comm)
Reduce!(sendbuf, recvbuf, op, comm::Comm; root::Integer=Cint(0)) =
    Reduce!(sendbuf, recvbuf, op, root, comm)

function Reduce!(rbuf::RBuffer, op::Union{Op,MPI_Op}, root::Integer, comm::Comm)
    # int MPI_Reduce(const void* sendbuf, void* recvbuf, int count,
    #                MPI_Datatype datatype, MPI_Op op, int root, MPI_Comm comm)
    @mpichk ccall((:MPI_Reduce, libmpi), Cint,
                  (MPIPtr, MPIPtr, Cint, MPI_Datatype, MPI_Op, Cint, MPI_Comm),
                  rbuf.senddata, rbuf.recvdata, rbuf.count, rbuf.datatype, op, root, comm)
    return rbuf.recvdata
end

# Convert user-provided functions to MPI.Op
function Reduce!(rbuf::RBuffer, op, root::Integer, comm::Comm)
    Reduce!(rbuf, Op(op, eltype(rbuf)), root, comm)
end
Reduce!(sendbuf, recvbuf, op, root::Integer, comm::Comm) =
    Reduce!(RBuffer(sendbuf, recvbuf), op, root, comm)

# inplace
function Reduce!(buf, op, root::Integer, comm::Comm)
    if Comm_rank(comm) == root
        Reduce!(IN_PLACE, buf, op, root, comm)
    else
        Reduce!(buf, nothing, op, root, comm)
    end
end

# allocating
"""
    recvbuf = Reduce(sendbuf, op, comm::Comm; root::Integer=0)

Performs elementwise reduction using the operator `op` on the buffer `sendbuf`, returning
the result `recvbuf` on the process of rank `root`, and `nothing` on non-root processes.

`sendbuf` can also be a scalar, in which case `recvbuf` will be a value of the same type.

# See also
- [`Reduce!`](@ref) for mutating and in-place operations
- [`Allreduce!`](@ref)/[`Allreduce`](@ref) to send reduction to all ranks.
- [`Op`](@ref) for details on reduction operators.

# External links
$(_doc_external("MPI_Reduce"))
"""
Reduce(sendbuf, op, comm::Comm; root::Integer=Cint(0)) =
    Reduce(sendbuf, op, root, comm)
function Reduce(sendbuf::AbstractArray, op, root::Integer, comm::Comm)
    if Comm_rank(comm) == root
        Reduce!(sendbuf, similar(sendbuf), op, root, comm)
    else
        Reduce!(sendbuf, nothing, op, root, comm)
    end
end
function Reduce(object::T, op, root::Integer, comm::Comm) where {T}
    source = Ref(object)
    if Comm_rank(comm) == root
        Reduce!(source, Ref{T}(), op, root, comm)[]
    else
        Reduce!(source, nothing, op, root, comm)
    end
end

## Allreduce

# mutating
"""
    Allreduce!(sendbuf, recvbuf, op, comm::Comm)
    Allreduce!(sendrecvbuf, op, comm::Comm)

Performs elementwise reduction using the operator `op` on the buffer `sendbuf`, storing
the result in the `recvbuf` of all processes in the group.

`Allreduce!` is equivalent to a [`Reduce!`](@ref) operation followed by
a [`Bcast!`](@ref), but can lead to better performance.

If only one `sendrecvbuf` buffer is provided, then the operation is performed in-place.

# See also
- [`Allreduce`](@ref), to handle allocation of the output buffer.
- [`Reduce!`](@ref)/[`Reduce`](@ref) to send reduction to a single rank.
- [`Op`](@ref) for details on reduction operators.

# External links
$(_doc_external("MPI_Allreduce"))
"""
function Allreduce!(rbuf::RBuffer, op::Union{Op,MPI_Op}, comm::Comm)
    # int MPI_Allreduce(const void* sendbuf, void* recvbuf, int count,
    #                   MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
    @mpichk ccall((:MPI_Allreduce, libmpi), Cint,
                  (MPIPtr, MPIPtr, Cint, MPI_Datatype, MPI_Op, MPI_Comm),
                  rbuf.senddata, rbuf.recvdata, rbuf.count, rbuf.datatype, op, comm)
    rbuf.recvdata
end
function Allreduce!(rbuf::RBuffer, op, comm::Comm)
    Allreduce!(rbuf, Op(op, eltype(rbuf)), comm)
end
Allreduce!(sendbuf, recvbuf, op, comm::Comm) =
    Allreduce!(RBuffer(sendbuf, recvbuf), op, comm)

# inplace
Allreduce!(buf, op, comm::Comm) =
    Allreduce!(IN_PLACE, buf, op, comm)

# allocating
"""
    recvbuf = Allreduce(sendbuf, op, comm)

Performs elementwise reduction using the operator `op` on the buffer `sendbuf`, returning
the result in the `recvbuf` of all processes in the group.

`sendbuf` can also be a scalar, in which case `recvbuf` will be a value of the same type.

# See also
- [`Allreduce!`](@ref) for mutating or in-place operations.
- [`Reduce!`](@ref)/[`Reduce`](@ref) to send reduction to a single rank.
- [`Op`](@ref) for details on reduction operators.

# External links
$(_doc_external("MPI_Allreduce"))
"""
Allreduce(sendbuf::AbstractArray, op, comm::Comm) =
    Allreduce!(sendbuf, similar(sendbuf), op, comm)
Allreduce(obj::T, op, comm::Comm) where {T} =
    Allreduce!(Ref(obj), Ref{T}(), op, comm)[]

## Scan

# mutating
"""
    Scan!(sendbuf, recvbuf, op, comm::Comm)
    Scan!(sendrecvbuf, op, comm::Comm)

Inclusive prefix reduction (analagous to `accumulate` in Julia): `recvbuf` on rank `i`
will contain the the result of reducing `sendbuf` by `op` from ranks `0:i`.

If only a single buffer `sendrecvbuf` is provided, then operations will be performed
in-place.

# See also
- [`Scan`](@ref) to handle allocation of the output buffer
- [`Exscan!`](@ref)/[`Exscan`](@ref) for exclusive scan
- [`Op`](@ref) for details on reduction operators.

# External links
$(_doc_external("MPI_Scan"))
"""
function Scan!(rbuf::RBuffer, op::Union{Op,MPI_Op}, comm::Comm)
    # int MPI_Scan(const void* sendbuf, void* recvbuf, int count,
    #              MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
    @mpichk ccall((:MPI_Scan, libmpi), Cint,
                  (MPIPtr, MPIPtr, Cint, MPI_Datatype, MPI_Op, MPI_Comm),
                  rbuf.senddata, rbuf.recvdata, rbuf.count, rbuf.datatype, op, comm)
    rbuf.recvdata
end
function Scan!(rbuf::RBuffer, op, comm::Comm)
    Scan!(rbuf, Op(op, eltype(rbuf)), comm)
end
Scan!(sendbuf, recvbuf, op, comm::Comm) =
    Scan!(RBuffer(sendbuf, recvbuf), op, comm)

# inplace
Scan!(buf, op, comm::Comm) =
    Scan!(IN_PLACE, buf, op, comm)

# allocating
"""
    recvbuf = Scan(sendbuf, op, comm::Comm)

Inclusive prefix reduction (analagous to `accumulate` in Julia): `recvbuf` on rank `i`
will contain the the result of reducing `sendbuf` by `op` from ranks `0:i`.

`sendbuf` can also be a scalar, in which case `recvbuf` will also be a scalar of the same
type.

# See also
- [`Scan!`](@ref) for mutating or in-place operations
- [`Exscan!`](@ref)/[`Exscan`](@ref) for exclusive scan
- [`Op`](@ref) for details on reduction operators.

# External links
$(_doc_external("MPI_Scan"))
"""
Scan(sendbuf::AbstractArray, op, comm::Comm) =
    Scan!(sendbuf, similar(sendbuf), op, comm)
Scan(object::T, op, comm::Comm) where {T} =
    Scan!(Ref(object), Ref{T}(), op, comm)[]

## Exscan
# mutating
"""
    Exscan!(sendbuf, recvbuf, op, comm::Comm)
    Exscan!(sendrecvbuf, op, comm::Comm)

Exclusive prefix reduction (analagous to `accumulate` in Julia): `recvbuf` on rank `i`
will contain the the result of reducing `sendbuf` by `op` from ranks `0:i-1`. The
`recvbuf` on rank `0` is ignored, and the `recvbuf` on rank `1` will contain the contents
of `sendbuf` on rank `0`.

If only a single `sendrecvbuf` is provided, then operations are performed in-place, and
`buf` on rank 0 will remain unchanged.

# See also
- [`Exscan`](@ref) to handle allocation of the output buffer
- [`Scan!`](@ref)/[`Scan`](@ref) for inclusive scan
- [`Op`](@ref) for details on reduction operators.

# External links
$(_doc_external("MPI_Exscan"))
"""
function Exscan!(rbuf::RBuffer, op::Union{Op,MPI_Op}, comm::Comm)
    # int MPI_Exscan(const void* sendbuf, void* recvbuf, int count,
    #                MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
    @mpichk ccall((:MPI_Exscan, libmpi), Cint,
          (MPIPtr, MPIPtr, Cint, MPI_Datatype, MPI_Op, MPI_Comm),
          rbuf.senddata, rbuf.recvdata, rbuf.count, rbuf.datatype, op, comm)
    rbuf.recvdata
end
function Exscan!(rbuf::RBuffer, op, comm::Comm)
    Exscan!(rbuf, Op(op, eltype(rbuf)), comm)
end
Exscan!(sendbuf, recvbuf, op, comm::Comm) =
    Exscan!(RBuffer(sendbuf, recvbuf), op, comm)

# inplace
Exscan!(buf, op, comm::Comm) =
    Exscan!(IN_PLACE, buf, op, comm)

# allocating
"""
    recvbuf = Exscan(sendbuf, op, comm::Comm)

Exclusive prefix reduction (analagous to `accumulate` in Julia): `recvbuf` on rank `i`
will contain the the result of reducing `sendbuf` by `op` from ranks `0:i-1`. The
`recvbuf` on rank `0` is undefined, and the `recvbuf` on rank `1` will contain the contents
of `sendbuf` on rank `0`.

# See also
- [`Exscan!`](@ref) for mutating and in-place operations
- [`Scan!`](@ref)/[`Scan`](@ref) for inclusive scan
- [`Op`](@ref) for details on reduction operators.

# External links
$(_doc_external("MPI_Scan"))
"""
Exscan(sendbuf::AbstractArray, op, comm::Comm) =
    Exscan!(sendbuf, similar(sendbuf), op, comm)
Exscan(object::T, op, comm::Comm) where {T} =
    Exscan!(Ref(object), Ref{T}(), op, comm)[]


text:
# provides some helper tools to deal with NA elements (i.e. elements of type Nullable{T}()).

import DataFrames: DataFrame

"""

A thin wrapper around AbstractArray. The reason to introduce this wrapper is to redefine
the dotted operators such as .+, .-. Those operators will be mapped to arrays, elementwise just as before,
but, if each element is null, those operators will be applied to the one inside the Nullable.
For example,

```julia
julia> AbstractArrayWrapper([Nullable(1), Nullable(2)]) .+ AbstractArrayWrapper([Nullable{Int}(), Nullable(3)])
2-element MultidimensionalTables.AbstractArrayWrapper{Nullable{Int64},1,Array{Nullable{Int64},1}}:
 Nullable{Int64}()
 Nullable(5)
```

Note that this means lifting those dotted operators via the list(AbstractArray) and maybe(Nullable) functors.

It is possible to redefine those operators for AbstractArray, but concerning about compatibility, it may be
best to introduce a new wrapper class for that.

"""
immutable AbstractArrayWrapper{T,N,A<:AbstractArray} <: AbstractArray{T,N}
  a::A #AbstractArray{T,N} # this is a thin array wrapper
end

AbstractArrayWrapper(a) = AbstractArrayWrapper{eltype(a),ndims(a),typeof(a)}(a)

# delegate the functionality of a type into its component.
# use it when you simply wrap some other data type, but want to use the funcitonality defiend for that data type
# for the new type.
macro delegate(source, targets...)
  if !(:head in fieldnames(source) && source.head == :.)
    throw(ArgumentError("cannot create a delegate statement from the macro @delegate."))
  else
    typename = source.args[1]
    fieldname = source.args[2].args[1]
    targetexpr = map(targets) do target
      Expr(:(=),
           Expr(:call, target, Expr(:(::), :x, typename), Expr(:..., :args)),
           Expr(:call, target, Expr(:., :x, QuoteNode(fieldname)), Expr(:..., :args)))
    end
    esc(Expr(:block, targetexpr...))
  end
end

# similar to delegate, but wraps the final result again using the new data type.
macro delegate_and_lift(source, targets...)
  if !(:head in fieldnames(source) && source.head == :.)
    throw(ArgumentError("cannot create a delegate statement from the macro @delegate."))
  else
    typename = source.args[1]
    fieldname = source.args[2].args[1]
    targetexpr = map(targets) do target
      Expr(:(=),
           Expr(:call, esc(target), Expr(:(::), :x, typename), Expr(:..., :args)),
           Expr(:call, typename, Expr(:call, esc(target), Expr(:., :x, QuoteNode(fieldname)), Expr(:..., :args))))
    end
    Expr(:block, targetexpr...)
  end
end

# used in naliftexp.
naliftexpriter(expriter) = begin
  if all(x->:head in fieldnames(x) && x.head == :row, expriter)
    begin
      vcat_expr = vcat(map(x->x.args, expriter)...)
      test_expr = filter(x->x!=:NA, vcat_expr)
      type_expr = Expr(:call, :eltype, Expr(:vect, test_expr...))
      newexpriter = map(expriter) do expr
        rowexpr = map(expr.args) do subexpr
          if subexpr == :NA
            Expr(:call, Expr(:curly, :Nullable, esc(type_expr)))
          else
            naliftexp(subexpr)
          end
        end
        Expr(:row, rowexpr...)
      end
      newexpriter
    end
  else
    begin
      test_expr = filter(x -> x!=:NA, expriter)
      type_expr = Expr(:call, :eltype, Expr(:vect, test_expr...))
      newexpriter = map(expriter) do expr
        if expr == :NA
          Expr(:call, Expr(:curly, :Nullable, esc(type_expr)))
        else
          naliftexp(expr)
        end
      end
      newexpriter
    end
  end
end

# reads an array expression and replaces a marker NA with an appropriate Nullable{T}().
# the way it does it now is to create an array without that NA marker, and check the resulting array type (say, T),
# then replace the NA marker with Nullable{T}().
naliftexp(expr) = begin
  new_expr =
    if :head in fieldnames(expr)
      if expr.head == :vect || expr.head == :vcat || expr.head == :hcat || expr.head == :row
        new_args = naliftexpriter(expr.args)
        Expr(expr.head, new_args...)
      elseif (expr.head == :ref && isupper(string(expr.args[1])[1])) ||
             expr.head == :typed_vcat ||
             expr.head == :typed_hcat
        Expr(expr.head, expr.args[1], naliftexpriter(expr.args[2:end])...)
      else
        nalift_expwrap(expr)
      end
    else
      nalift_expwrap(expr)
    end
  new_expr
end

nalift_expwrap(expr) = Expr(:call, :nalift_nested, esc(expr))

# a simple wrapper of a floating point array. For floating points, there is a dedicated NaN value,
# so there is no need to use the Nullable element type such as in Array{Nullable{Float64}}.
immutable FloatNAArray{T<:AbstractFloat,N,A} <: AbstractArray{Nullable{T},N}
  data::A #AbstractArray{T,N}
end

FloatNAArray{T<:AbstractFloat,N}(data::AbstractArray{T,N}) = FloatNAArray{T,N,typeof(data)}(data)
FloatNAArray{T<:AbstractFloat,N}(data::AbstractArray{Nullable{T},N}) = begin
  nulldata = convert(T, NaN)
  projdata = map(x->x.isnull ? nulldata : x.value, data)
  FloatNAArray(projdata)
end
Base.eltype{T<:AbstractFloat,N,A}(::Type{FloatNAArray{T,N,A}}) = Nullable{T}
Base.getindex{T<:AbstractFloat}(arr::FloatNAArray{T}, arg::Int) = (r=getindex(arr.data, arg);isnan(r) ? Nullable{T}() : Nullable(r))
Base.getindex{T<:AbstractFloat}(arr::FloatNAArray{T}, arg::CartesianIndex) = (r=getindex(arr.data, arg);isnan(r) ? Nullable{T}() : Nullable(r))
Base.getindex{T<:AbstractFloat}(arr::FloatNAArray{T}, args::Int...) = (r=getindex(arr.data, args...);isnan(r) ? Nullable{T}() : Nullable(r))
Base.getindex{T<:AbstractFloat}(arr::FloatNAArray{T}, args...) = FloatNAArray(getindex(arr.data, args...))
Base.setindex!{T<:AbstractFloat}(arr::FloatNAArray{T}, v::Nullable, arg::Int) = (nulldata=convert(T,NaN);setindex!(arr.data, if v.isnull;nulldata else v.value end, arg))
Base.setindex!{T<:AbstractFloat}(arr::FloatNAArray{T}, v::Nullable, args::Int...) = (nulldata=convert(T,NaN);setindex!(arr.data, if v.isnull;nulldata else v.value end, args...))
Base.setindex!{T<:AbstractFloat}(arr::FloatNAArray{T}, v::Nullable, args...) = (nulldata=convert(T,NaN);setindex!(arr.data, if v.isnull;nulldata else v.value end, args...))
Base.setindex!{T<:AbstractFloat}(arr::FloatNAArray{T}, v::AbstractFloat, args::Int...) = setindex!(arr.data, v, args...)
Base.setindex!{T<:AbstractFloat}(arr::FloatNAArray{T}, v::AbstractFloat, arg::Int) = setindex!(arr.data, v, arg)
Base.setindex!{T<:AbstractFloat}(arr::FloatNAArray{T}, v::AbstractFloat, arg::Int) = setindex!(arr.data, v, arg)
Base.setindex!{T<:AbstractFloat,V<:AbstractFloat}(arr::FloatNAArray{T}, v::AbstractArray{V}, args...) = setindex!(arr.data, v, args...)
Base.setindex!{T<:AbstractFloat,V<:Nullable}(arr::FloatNAArray{T}, v::AbstractArray{V}, args...) = (nulldata=convert(T,NaN);setindex!(arr.data, map(x->if x.isnull;nulldata else x.value end, v), args...))

Base.reshape(arr::FloatNAArray, args::Tuple{Vararg{Int}}) = FloatNAArray(reshape(arr.data, args))
Base.reshape(arr::FloatNAArray, args::Int...) = FloatNAArray(reshape(arr.data, args...))
Base.repeat(arr::FloatNAArray; kwargs...) = FloatNAArray(repeat(arr.data; kwargs...))

Base.next{T}(arr::FloatNAArray{T}, state) = ((x,ns)=next(arr.data, state);isnan(x) ? (Nullable{T}(),ns) : (Nullable(x),ns))
Base.linearindexing{T<:AbstractFloat,N,A}(::Type{FloatNAArray{T,N,A}}) = Base.linearindexing(A)
Base.sub(arr::FloatNAArray, args::Union{Colon,Int,AbstractVector}...) = FloatNAArray(sub(arr.data, args...))
Base.slice(arr::FloatNAArray, args::Union{Colon,Int,AbstractVector}...) = FloatNAArray(slice(arr.data, args...))
Base.sub(arr::FloatNAArray, args::Tuple{Vararg{Union{Colon,Int,AbstractVector}}})= FloatNAArray(sub(arr.data, args...))
Base.slice(arr::FloatNAArray, args::Tuple{Vararg{Union{Colon,Int,AbstractVector}}}) = FloatNAArray(slice(arr.data, args...))
@delegate(FloatNAArray.data, Base.start, Base.done, Base.size, Base.find)
@delegate_and_lift(FloatNAArray.data, Base.transpose, Base.permutedims, Base.repeat, Base.reshape, Base.sort, Base.sort!, Base.reverse,
                                      Base.sub, Base.slice)
Base.similar{T,N}(arr::FloatNAArray, ::Type{T}, dims::NTuple{N,Int}) = similar(arr.data, T, dims)
Base.similar{T<:AbstractFloat,N}(arr::FloatNAArray, ::Type{Nullable{T}}, dims::NTuple{N,Int}) = FloatNAArray(similar(arr.data, T, dims))
Base.map(f::Function, arr0::FloatNAArray, arrs::AbstractArray...) = begin
  if isempty(arr0)
    return similar(arr0, Nullable{Any})
  end
  firstelem = f(arr0[1], map(x->x[1], arrs)...)
  returntype = typeof(firstelem)
  result = similar(arr0, returntype)
  floatnaarray_map_inner!(result, firstelem, f, arr0, arrs)
  result
end

floatnaarray_map_inner!(result::AbstractArray, firstelem, f::Function, arr0::FloatNAArray, arrs) = begin
  after_first = false
  result[1] = firstelem #f(arr0[1], map(x->x[1], arrs)...)
  for i in eachindex(arr0,arrs...)
    if after_first
      result[i] = f(arr0[i], map(x->x[i], arrs)...)
    end
    after_first = true
  end
end

Base.repeat(arr::FloatNAArray; kwargs...) = FloatNAArray(repeat(arr.data; kwargs...))
Base.sort(arr::FloatNAArray; kwargs...) = FloatNAArray(sort(arr.data; kwargs...))
Base.sort!(arr::FloatNAArray; kwargs...) = FloatNAArray(sort!(arr.data; kwargs...))
Base.cat(dim::Int, arr::FloatNAArray, arrs::FloatNAArray...) = FloatNAArray(cat(dim, arr.data, map(x->x.data, arrs)))
Base.vcat(arrs::FloatNAArray...) = cat(1, arrs...)
Base.hcat(arrs::FloatNAArray...) = cat(2, arrs...)
Base.sub(arr::FloatNAArray, args::Union{Base.Colon,Int,AbstractVector}...) = FloatNAArray(sub(arr.data, args...))
Base.slice(arr::FloatNAArray, args::Union{Base.Colon,Int,AbstractVector}...) = FloatNAArray(slice(arr.data, args...))
simplify_floatarray(arr::FloatNAArray) = arr
simplify_floatarray{T<:AbstractFloat,N,A}(arr::FloatNAArray{T,N,A}) = arr
simplify_floatarray{T<:AbstractFloat,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}) = arr
simplify_floatarray{T<:AbstractFloat}(arr::AbstractArrayWrapper{Nullable{T}}) =
  AbstractArrayWrapper(simplify_floatarray(arr.a))
simplify_floatarray{T<:AbstractFloat,N}(arr::AbstractArray{Nullable{T},N}) = begin
  if isempty(arr)
    # for some reason, if arr is empty, projected has the type Nullable{T}[] !
    # isempty(projected) ? FloatNAArray(Array(T,fill(0,N)...)) : FloatNAArray(projected)
    FloatNAArray(zeros(T, size(arr)))
  else
    FloatNAArray(map(arr) do elem
      if elem.isnull
        convert(T, NaN)::T
      else
        elem.value::T
      end
    end)
  end
end
simplify_floatarray(arr) = arr
# simplify things. Ideally it can perform some other checks,
# currently, it constrains the element type further if possible,
# use FloatNAArray if necessary, and wrap it over AbstractArrayWrapper if necessary.
simplify_array(arr::AbstractArray) = wrap_array(simplify_floatarray(type_array(arr)))
simplify_array(nonarray) = nonarray

"""

`nalift(arr)`

Lift each element in an array `arr` to `Nullable` if it is not already so.
Unlike `@nalift`, it does not perform lifting recursively.
It returns `arr` itself when applied to a `DictArray`/`LabeledArray`.

##### Examples

```julia
julia> nalift(Any[[1,2,3],[4,5]])
2-element MultidimensionalTables.AbstractArrayWrapper{Nullable{Array{Int64,1}},1,Array{Nullable{Array{Int64,1}},1}}:
 Nullable([1,2,3])
 Nullable([4,5])

julia> nalift([1,2,3])
3-element MultidimensionalTables.AbstractArrayWrapper{Nullable{Int64},1,Array{Nullable{Int64},1}}:
 Nullable(1)
 Nullable(2)
 Nullable(3)

julia> nalift(Any[[1,2,3],[4,5]])
2-element MultidimensionalTables.AbstractArrayWrapper{Nullable{Array{Int64,1}},1,Array{Nullable{Array{Int64,1}},1}}:
 Nullable([1,2,3])
 Nullable([4,5])

julia> nalift(darr(a=[1 2;3 4;5 6], b=[:x :y;:z :w;:u :v]))
3 x 2 DictArray

a b |a b
----+----
1 x |2 y
3 z |4 w
5 u |6 v
```

"""
function nalift end

nalift(x::Nullable) = x
nalift{T}(x::AbstractArray{Nullable{T}}) = wrap_array(x)
nalift{F<:AbstractFloat}(x::FloatNAArray{F}) = wrap_array(x)
nalift{F<:AbstractFloat}(x::AbstractArray{F}) = wrap_array(FloatNAArray(x))
nalift(x::AbstractArray) = begin
  if isempty(x)
    wrap_array(similar(x, Nullable{eltype(x)}))
  elseif isa(x[1], Nullable) && all(elem->isa(elem, Nullable), x)
    # something is fishy... Even though x does not claim it is Nullable, its elements
    #@show "something is fishy... Even though x does not claim it is Nullable, its elements"
    #Nullable[x...]
    result = similar(x, Nullable)
    copy!(result, x)
    wrap_array(result)
  else
    wrap_array(map(Nullable, x))
  end
end
nalift(x::DictArray) = x
nalift(x::LabeledArray) = x
nalift(x::DataFrame) = x
nalift{T,N,AXES,TN}(x::AbstractArray{LabeledArray{T,N,AXES,TN}}) = x
nalift{K,N,VS,SV}(x::AbstractArray{DictArray{K,N,VS,SV}}) = x
nalift(x) = Nullable(x)

# nalift nonmacro version. Read the input array and put Nullable wrappers appropriately.
nalift_nested(x::Nullable) = begin
  if !x.isnull && isa(x.value, AbstractArray)
    Nullable(simplify_array(map(elem->nalift(nalift_nested(elem)), x.value)))
  else
    x
  end
end
nalift_nested(x::DictArray) = x
nalift_nested(x::LabeledArray) = x
nalift_nested(x::DataFrame) = x
nalift_nested{T<:AbstractFloat}(x::AbstractArray{T}) = wrap_array(FloatNAArray(x))
nalift_nested{T}(x::AbstractArray{Nullable{T}}) = begin
  if isempty(x)
    x
  elseif isa(x[1].value, AbstractArray)
    simplify_array(map(elem->Nullable(nalift(nalift_nested(elem.value))), x))
  else
    simplify_array(x)
  end
end
nalift_nested(x::AbstractArray{Nullable}) = begin
  if isempty(x)
    x
  elseif isa(x[1].value, AbstractArray)
    simplify_array(map(elem->Nullable(nalift(nalift_nested(elem.value))), x))
  else
    simplify_array(x)
  end
end
#
#nalift_nested(x::AbstractArray{Nullable}) = simplify_array(x)
nalift_nested(x::AbstractArray) = begin
  if isempty(x)
    wrap_array(similar(x, Nullable{eltype(x)}))
  else
    r = map(elem->nalift_nested(elem), x)
    simplify_array(r)
  end
end
nalift_nested(x) = Nullable(x)

"""

`@nalift(arr)`

Lift each element in an array `arr` to `Nullable` if it is not already so.
It is mainly used to translate a manually typed array expression such as `[1,2,3,NA,5]` into a `Nullable` array.
Unlike `nalift`, it performs lifting recursively.
It returns `arr` itself when applied to a `DictArray`/`LabeledArray`.

##### Examples

```julia
julia> @nalift([1,2,3])
3-element MultidimensionalTables.AbstractArrayWrapper{Nullable{Int64},1,Array{Nullable{Int64},1}}:
 Nullable(1)
 Nullable(2)
 Nullable(3)

julia> @nalift([1,2,NA])
3-element MultidimensionalTables.AbstractArrayWrapper{Nullable{Int64},1,Array{Nullable{Int64},1}}:
 Nullable(1)
 Nullable(2)
 Nullable{Int64}()

julia> @nalift(Any[[1,2,3],[NA,5]])
2-element MultidimensionalTables.AbstractArrayWrapper{Nullable{MultidimensionalTables.AbstractArrayWrapper{Nullable{Int64},1,Array{Nullable{Int64},1}}},1,Array{Nullable{MultidimensionalTables.AbstractArrayWrapper{Nullable{Int64},1,Array{Nullable{Int64},1}}},1}}:
 Nullable([Nullable(1),Nullable(2),Nullable(3)])
 Nullable([Nullable{Int64}(),Nullable(5)])

julia> @nalift(larr(a=[1 2;3 4;5 6], b=[:x :y;:z :w;:u :v]))
3 x 2 LabeledArray

  |1   |2
--+----+----
  |a b |a b
--+----+----
1 |1 x |2 y
2 |3 z |4 w
3 |5 u |6 v
```

"""
macro nalift(expr)
  new_expr = naliftexp(expr)
  quote
    nalift(simplify_array(recursive_wrap_array($new_expr)))
  end
end

# had to devise this hackish function...
recursive_wrap_array{T<:AbstractArray}(x::AbstractArray{Nullable{T}}) = begin
  if T<:AbstractArray
    map(x->x.isnull ? x : Nullable(recursive_wrap_array(x.value)), x)
  else
    wrap_array(x)
  end
end
recursive_wrap_array(x::AbstractArray) = begin
  if (eltype(x) <: Nullable) || (!isempty(x) && isa(x[1], Nullable)) #type(typeof(x[1])) <: AbstractArray)
    wrap_array(map(x->x.isnull ? x : Nullable(recursive_wrap_array(x.value)), x))
  else
    wrap_array(map(recursive_wrap_array, x))
  end
end
recursive_wrap_array(x) = x
recursive_wrap_array(x::DictArray) = x
recursive_wrap_array(x::LabeledArray) = x
recursive_wrap_array{T,N,A<:FloatNAArray}(x::AbstractArrayWrapper{T,N,A}) = x

type NAElementException <: Exception end

"""

`setna!(arr, args...)`

Set the element of an array `arr` at `args` to `NA`.
If `args...` is omitted, all elements are set to `NA`.

* If `arr` is an array of element type `Nullable{T}`, `NA` means `Nullable{T}()`.
* If `arr` is a `DictArray`, `NA` means all fields at that position are `NA`.
* If `arr` is a `LabeledArray`, `NA` means the base of `arr` at that position is `NA`.

##### Examples

```julia
julia> setna!(@nalift([1,2,NA,4,5]))
5-element MultidimensionalTables.AbstractArrayWrapper{Nullable{Int64},1,Array{Nullable{Int64},1}}:
 Nullable{Int64}()
 Nullable{Int64}()
 Nullable{Int64}()
 Nullable{Int64}()
 Nullable{Int64}()

julia> setna!(@nalift([1,2,NA,4,5]), 2)
5-element Array{Nullable{Int64},1}:
 Nullable(1)
 Nullable{Int64}()
 Nullable{Int64}()
 Nullable(4)
 Nullable(5)

julia> setna!(@darr(a=[1 2 NA;4 5 6], b=[:x :y :z;:u :v :w]), 1:2, 1)
2 x 3 DictArray

a b |a b |a b
----+----+----
    |2 y |  z
    |5 v |6 w

julia> setna!(larr(a=[1 2 3;4 5 6], b=[:x :y :z;:u :v :w], axis1=[:X,:Y]), 1, 2:3)
2 x 3 LabeledArray

  |1   |2   |3
--+----+----+----
  |a b |a b |a b
--+----+----+----
X |1 x |    |
Y |4 u |5 v |6 w
```

"""
function setna! end

setna!{T}(arr::AbstractArray{Nullable{T}}) = fill!(arr, Nullable{T}())
setna!(arr::AbstractArray{Nullable}) = fill!(arr, Nullable{Any}())
setna!{T}(arr::AbstractArray{Nullable{T}}, args...) = setindex!(arr, Nullable{T}(), args...)
setna!(arr::AbstractArray{Nullable}, args...) = setindex!(arr, Nullable{Any}(), args...)
setna!{T<:AbstractFloat}(arr::FloatNAArray{T}, args...) = setindex!(arr.data, convert(T,NaN), args...)
setna!{T<:AbstractFloat}(arr::FloatNAArray{T}) = fill!(arr.data, convert(T,NaN))
setna!(arr::AbstractArrayWrapper, args...) = setna!(arr.a, args...)
setna!(arr::DictArray, args...) = begin
  map(tgt -> setna!(tgt, args...), arr.data.values)
  arr
end
setna!(arr::LabeledArray, args...) = begin
  setna!(arr.data, args...)
  arr
end

"""

`igna(arr [, nareplace])`

Ignore null elements from `arr`.
Null elements will be replaced by `nareplace`, if provided.
If not, the behavior is implementation specific: depending on the array type, it may give some default value or raise an error.
Most likely, a nullable element in an array of `Nullable{F}` element type for some `AbstractFloat` `F` can be replaced by a version of `NaN`.
But for other types, it may be better to raise an error.

* `igna(arr::AbstractArray{Nullable{T},N} [, na_replace])`: ignores null elements from `arr` and return an `AbstractArray{T,N}`. A null value is replaced by `na_replace` if provided. Otherwise, the result is implementation specific.

* `igna(ldict::LDict [, na_replace])` ignores null values from `ldict` and replace them with `na_replace` if provided. Otherwise, the result is implementation specific.

##### Examples

```julia
julia> igna(@nalift([1,2,NA,4,5]))
ERROR: MultidimensionalTables.NAElementException()
 in anonymous at /Users/changsoonpark/.julia/v0.4/MultidimensionalTables/src/na/na.jl:315
 in map_to! at abstractarray.jl:1289
 in map at abstractarray.jl:1311
 in igna at /Users/changsoonpark/.julia/v0.4/MultidimensionalTables/src/na/na.jl:313

julia> igna(@nalift([1.0,2.0,NA,4.0,5.0]))
5-element MultidimensionalTables.AbstractArrayWrapper{Float64,1,Array{Float64,1}}:
   1.0
   2.0
 NaN
   4.0
   5.0

julia> igna(@nalift([1,2,NA,4,5]), 3)
5-element MultidimensionalTables.AbstractArrayWrapper{Int64,1,Array{Int64,1}}:
 1
 2
 3
 4
 5

julia> igna(LDict(:a=>Nullable(3), :b=>Nullable{Int}()), 1)
MultidimensionalTables.LDict{Symbol,Int64} with 2 entries:
  :a => 3
  :b => 1
```

"""
function igna end

igna{T}(arr::AbstractArray{Nullable{T}}) = if isempty(arr)
  similar(arr, T)
else
  map(arr) do elem
    if elem.isnull
      throw(NAElementException())
    else
      elem.value
    end
  end
end
igna(arr::AbstractArray{Nullable}) = map(arr) do elem
  if elem.isnull
    throw(NAElementException())
  else
    elem.value
  end
end
igna{T}(arr::AbstractArray{Nullable{T}}, na_replace::T) = if isempty(arr)
  similar(arr, T)
else
  map(arr) do elem
    if elem.isnull
      na_replace
    else
      elem.value
    end
  end
end
igna(arr::AbstractArray{Nullable}, na_replace) = map(arr) do elem
  if elem.isnull
    na_replace
  else
    elem.value
  end
end
igna{T<:AbstractFloat}(arr::FloatNAArray{T}) = arr.data
igna{T<:AbstractFloat}(arr::FloatNAArray{T}, na_replace::T) = if isempty(arr)
  similar(arr, T)
else
  map(arr.data) do elem
    isnan(elem) ? na_replace : elem
  end
end
igna{T<:AbstractFloat}(arr::AbstractArray{Nullable{T}}) = igna(arr, convert(T, NaN))
igna{T<:AbstractFloat}(arr::AbstractArray{Nullable{T}}, na_replace::T) = if isempty(arr)
  similar(arr, T)
else
  map(arr) do elem
    if elem.isnull; na_replace else elem.value end
  end
end
igna{T}(elem::Nullable{T}, na_replace::T) = elem.isnull ? na_replace : elem.value

"""

`ignabool(arr)`

Ignore the `Nullable` part of of either a `Nullable` array or a `Nullable` variable.
It is mainly used in the condition statement for @select or @update, where it assumes that only Nullable(true) chooses the element.  Nullable(false) or Nullable{T}() will be regarded as false.

* `ignabool(::AbstractArray{Nullable{Bool}}) returns an `AbstractArray{Bool}` where null and `Nullable(false)` elements are converted into `false` and `Nullable(true)` into `true`.
* `ignabool(::Nullable{Bool})` converts null and `Nullable(false)` elements into `false` and `Nullable(true)` into true.

##### Examples

```julia
julia> ignabool(Nullable{Bool}())
false

julia> ignabool(Nullable(true))
true

julia> ignabool(Nullable(false))
false

julia> ignabool(@nalift([true true NA;false NA true]))
2x3 MultidimensionalTables.AbstractArrayWrapper{Bool,2,Array{Bool,2}}:
  true   true  false
 false  false   true
```

"""
function ignabool end

ignabool(arr::AbstractArrayWrapper{Nullable{Bool}}) = AbstractArrayWrapper(ignabool(arr.a))
# this turns to be definitely faster than `map`.
ignabool(arr::AbstractArray{Nullable{Bool}}) = begin
  result = similar(arr, Bool)
  for i in eachindex(arr)
    @inbounds result[i] = ignabool(arr[i])
  end
  result
end
#ignabool(arr::AbstractArray{Nullable{Bool}}) = map_typed(ignabool, Bool, arr)
ignabool(elem::Nullable{Bool}) = !elem.isnull && elem.value
ignabool(arr::AbstractArray{Bool}) = arr
ignabool(elem::Bool) = elem

"""

`isna(arr [, coords...])`

Checks `NA` for each element and produces an AbstractArray{Bool} of the same shape as `arr`.
If `coords...` are provided, `isna` checks `NA` at that position.

* If an input array is `AbstractArray{Nullable{T}}`, it checkes whether an element is null.
* If an input array is `DictArray`, it tests whether all values of the dictionary values for each element are null.
* If an input array is `LabeledArray`, it applies `isna` to the base of `arr`.

##### Examples

```julia
julia> t = @darr(a=[1 NA 3;4 5 NA], b=[NA NA :z;:u :v :w])
2 x 3 DictArray

a b |a b |a b
----+----+----
1   |    |3 z
4 u |5 v |  w

julia> isna(t)
2x3 MultidimensionalTables.AbstractArrayWrapper{Bool,2,Array{Bool,2}}:
 false   true  false
 false  false  false

julia> isna(t, 2, 2:3)
1x2 MultidimensionalTables.AbstractArrayWrapper{Bool,2,Array{Bool,2}}:
 false  false

julia> isna(@larr(t, axis1[NA,:Y], axis2[NA,NA,"W"]))
2x3 MultidimensionalTables.AbstractArrayWrapper{Bool,2,Array{Bool,2}}:
 false   true  false
 false  false  false

julia> isna(@nalift([1 2 NA;NA 5 6]))
2x3 MultidimensionalTables.AbstractArrayWrapper{Bool,2,Array{Bool,2}}:
 false  false   true
  true  false  false
```

"""
function isna end

isna{T}(arr::AbstractArray{Nullable{T}}) = map(elem->elem.isnull, arr)
isna(arr::AbstractArray{Nullable}) = map(elem->elem.isnull, arr)
isna{T<:AbstractFloat}(arr::FloatNAArray{T}) = map(isnan, arr.data)
isna{T,N,A}(arr::AbstractArrayWrapper{T,N,A}) = AbstractArrayWrapper(isna(arr.a))
isna{T<:AbstractFloat}(arr::AbstractArray{Nullable{T}}) = map(elem->elem.isnull || isnan(elem.value), arr)
isna(x::Nullable) = x.isnull
isna(arr::AbstractArray, coords...) = isna(arr[coords...])

# this is a little bit out of place...
(==)(::AbstractArrayWrapper, ::DefaultAxis) = false
(==)(::DefaultAxis, ::AbstractArrayWrapper) = false

ext: jl


text:
export continuous
function continuous(map_fn, steps::Int, initial_state::Vector{Float32}=Float32[ 1 ])
    continuous(map_fn, UInt64(steps), initial_state)
end
function continuous(map_fn, steps::UInt64, initial_state::Vector{Float32})
    current_state = initial_state
    uns::Float32 = 0
    history = AutomataContinuousHistory()

    push!(history, current_state, uns)

    for i = 2:steps
        next_state = zeros(Float32, length(current_state) + 2)
        for i in 1:length(next_state)
            next_state[i] = map_fn(state_at(current_state, uns, i-2) + state_at(current_state, uns, i-1) + state_at(current_state, uns, i))
        end
        uns = map_fn(uns * 3)
        current_state = next_state

        push!(history, next_state, uns)
    end

    return history
end


text:
module UnstructuredGeometry

# Dependencies of this module

using Gridap
using Gridap.Helpers
using Gridap.CellValuesGallery
using Gridap.GridGraphs: FullGridGraphFromData
using UnstructuredGrids
using UnstructuredGrids.Kernels: rewind_ptrs!
using UnstructuredGrids.Kernels: length_to_ptrs!
using UnstructuredGrids.Kernels: generate_data_and_ptrs

# Functionality provided by this module

export rewind_ptrs!
export length_to_ptrs!
export generate_data_and_ptrs
export UnstructuredGrid
export FlexibleUnstructuredGrid
import Gridap: points, cells, celltypes, cellorders
import Gridap: FullGridGraph
export cellsdata, cellsptrs
export UGrid
import UnstructuredGrids: UGrid

"""
Struct representing an unstructured grid with efficient memory layout
"""
struct UnstructuredGrid{
  D,Z,P<:AbstractVector{Point{D,Float64}},A<:IndexCellValue{NTuple{Z,Int}},B<:IndexCellValue{Int}} <: Grid{D,Z}
  points::P
  cells_data::Vector{Int}
  cells_ptrs::Vector{Int}
  ctypes::A
  corders::B
end

function UnstructuredGrid(
  points::AbstractArray{Float64,2},
  cells_data,
  cells_ptrs,
  ctypes,
  corders)

  dim, npoints = size(points)
  k = reshape(points,(dim*npoints,))
  _points = reinterpret(Point{dim,Float64},k)

  UnstructuredGrid(
    _points,
    cells_data,
    cells_ptrs,
    ctypes,
    corders)
end

function UnstructuredGrid(grid::Grid)
  @notimplemented
end

UnstructuredGrid(grid::UnstructuredGrid) = grid

points(self::UnstructuredGrid) = CellValueFromArray(self.points)

cells(self::UnstructuredGrid) = CellVectorFromDataAndPtrs(self.cells_data,self.cells_ptrs)

cellsdata(self::UnstructuredGrid) = self.cells_data

cellsptrs(self::UnstructuredGrid) = self.cells_ptrs

celltypes(self::UnstructuredGrid) = self.ctypes

cellorders(self::UnstructuredGrid) = self.corders

"""
Create a UGrid from a UnstructuredGrid
"""
UGrid(grid::UnstructuredGrid) = _unstructured_grid_to_ugrid(grid)

function FullGridGraph(grid::UnstructuredGrid)
  ugrid = UGrid(grid)
  g = generate_full_grid_graph(ugrid)
  data = Array{IndexCellArray,2}(undef,size(g))
  for j in 1:size(g,2)
    for i in 1:size(g,1)
      # if i != j || (i==1 && j==1)
      x = g[i,j]
      data[i,j] = CellVectorFromDataAndPtrs(x.list,x.ptrs)
      # end
    end
  end
  FullGridGraphFromData(data)
end

"""
Struct representing an unstructured grid with inefficient memory layout
"""
struct FlexibleUnstructuredGrid{D,Z,A<:IndexCellValue{NTuple{Z,Int}},B<:IndexCellValue{Int}} <: Grid{D,Z}
  points::Vector{Point{D,Float64}}
  cells::Vector{Vector{Int}}
  ctypes::A
  corders::B
end

points(self::FlexibleUnstructuredGrid) = CellValueFromArray(self.points)

cells(self::FlexibleUnstructuredGrid) = CellValueFromArray(self.cells)

celltypes(self::FlexibleUnstructuredGrid) = self.ctypes

cellorders(self::FlexibleUnstructuredGrid) = self.corders

# Helpers

function _unstructured_grid_to_ugrid(grid::UnstructuredGrid{D}) where D

  x = points(grid)
  npoins = length(x)
  coords = reshape(reinterpret(Float64,x),(D,npoins))

  ctypes, refcells = _setup_ctypes_and_refcells(
    celltypes(grid), cellorders(grid))

  UGrid(
    cellsdata(grid),
    cellsptrs(grid),
    ctypes,
    refcells,
    coords)

end

function _setup_ctypes_and_refcells(cell_to_code, cell_to_order)
  @notimplemented
end

function _setup_ctypes_and_refcells(
  cell_to_code::ConstantCellValue,
  cell_to_order::ConstantCellValue)
  order = cell_to_order.value
  @notimplementedif order != 1
  code = cell_to_code.value
  polytope = Polytope(code)
  refcell = RefCell(polytope)
  ncells = length(cell_to_code)
  cell_to_ctype = ConstantCellValue(1,ncells)
  ctype_to_refcell = [refcell]
  (cell_to_ctype, ctype_to_refcell)
end

end # module Unstructured


text:
"""制御器のまとめ
"""
module Controller

using LinearAlgebra
using Parameters  # 構造体にキーワード引数をつけるためのモジュール
using MatrixEquations  #　行列方程式ソルバー
using ControlSystems  # 制御関係
#using Optim  # 最適化ライブラリ
#using JuMP
using ForwardDiff

include("dynamics.jl")
using .Dynamics: M, C, G, K, D, uncertain_K, uncertain_D, invM

export KinematicController
export PDandFBController
export PassivityBasedController
export PassivityBasedAdaptiveController
export SDREController
export MPCController
export calc_torque
export θp_dot

"""運動学制御器"""
@with_kw struct KinematicController{T}
    K_kin::Matrix{T}
    isUncertainty::Bool
end

"""入力を計算"""
function calc_torque(
    p::KinematicController{T},
    q::Vector{T}, q_dot::Vector{T},
    qd::Vector{T}, qd_dot::Vector{T}, qd_dot_dot::Vector{T},
    ) where T
    ek = qd .- q
    return p.K_kin * ek
end

"""PD + フィードバック線形化制御器"""
@with_kw struct PDandFBController{T}
    Kd::Matrix{T}
    Kp::Matrix{T}
    isUncertainty::Bool
end

"""入力を計算"""
function calc_torque(
    p::PDandFBController{T},
    q::Vector{T}, q_dot::Vector{T},
    qd::Vector{T}, qd_dot::Vector{T}, qd_dot_dot::Vector{T},
    ) where T

    α = M(q)
    q̃ = q .- qd
    q̃_dot = q_dot .- qd_dot

    τ_dash = qd_dot_dot .- p.Kd*q̃_dot .- p.Kp*q̃

    if p.isUncertainty  # 不確かさあり
        β = C(q, q_dot)*q_dot .+ uncertain_D*q_dot .+ uncertain_K*q .+ G(q)
    else
        β = C(q, q_dot)*q_dot .+ D*q_dot .+ K*q .+ G(q)
    end

    α*τ_dash .+ β
end

"""受動性に基づく制御器"""
@with_kw struct PassivityBasedController{T}
    Λ::Matrix{T}
    KG::Matrix{T}
    isUncertainty::Bool
end

"""入力を計算"""
function calc_torque(
    p::PassivityBasedController{T},
    q::Vector{T}, q_dot::Vector{T},
    qd::Vector{T}, qd_dot::Vector{T}, qd_dot_dot::Vector{T},
    ) where T
    #print("hoge")
    q̃ = q .- qd
    q̃_dot = q_dot .- qd_dot

    v = qd_dot .- p.Λ*q̃
    a = qd_dot_dot .- p.Λ*q̃_dot
    r = q̃_dot .+ p.Λ*q̃

    if p.isUncertainty  # 不確かさあり
        return M(q)*a .+ C(q, q_dot)*v .+ G(q) .+ uncertain_K*q .+ uncertain_D*v .- p.KG*r
    else
        return M(q)*a .+ C(q, q_dot)*v .+ G(q) .+ K*q .+ D*v .- p.KG*r
    end
end

"""受動性に基づく適応制御"""
@with_kw struct PassivityBasedAdaptiveController{T}
    invΓ::Matrix{T}
    Λ::Matrix{T}
    KG::Matrix{T}
    isUncertainty::Bool
end

"""パラメータの更新式"""
function θp_dot(
    p::PassivityBasedAdaptiveController{T},
    q::Vector{T}, q_dot::Vector{T}, qd::Vector{T}, qd_dot::Vector{T},
    ) where T

    q̃ = q .- qd
    q̃_dot = q_dot .- qd_dot

    v = qd_dot .- p.Λ*q̃
    r = q̃_dot .+ p.Λ*q̃

    Y = [diagm(q) diagm(v)]

    return -p.invΓ * Y' * r
end

"""入力を計算"""
function calc_torque(
    p::PassivityBasedAdaptiveController{T},
    q::Vector{T}, q_dot::Vector{T},
    qd::Vector{T}, qd_dot::Vector{T},
    θp::Vector{T}
    ) where T
    #print("hoge")

    q̃ = q .- qd
    q̃_dot = q_dot .- qd_dot

    v = qd_dot .- p.Λ*q̃

    r = q̃_dot .+ p.Λ*q̃

    Y = [diagm(q) diagm(v)]

    return Y * θp .- p.KG*r

end

"""SDRE法"""
@with_kw struct SDREController{T}
    Q::Matrix{T}
    R::Matrix{T}
    isUncertainty::Bool
end

"""可制御性判定"""
function isControllable(A, B)
    n, _ = size(A)
    _, m = size(B)
    #println(n, m)
    Co = Matrix{T}(undef, n, n*m)
    for i in 1:n
        Co[:, m*(i-1)+1:m*i] = (A^(i-1))*B
    end
    return rank(Co) == n
end

"""SDRE法の入力を計算"""
function calc_torque(
    p::SDREController{T},
    q::Vector{T}, q_dot::Vector{T},
    qd::Vector{T}, qd_dot::Vector{T}, qd_dot_dot::Vector{T},
    ) where T
    #println(uncertain_K)
    x = [
        qd .- q
        qd_dot .- q_dot
    ]
    invM = inv(M(q))

    if p.isUncertainty  # 不確かさありのとき
        A = [
            zeros(T, 3, 3) Matrix{T}(I, 3, 3)
            -invM*uncertain_K -invM*(C(q, q_dot) .+ uncertain_D)
        ]
        B = [
            zeros(T, 3, 3)
            invM
        ]
    else  # 不確かさ無し
        #println(K)
        A = [
            zeros(T, 3, 3) Matrix{T}(I, 3, 3)
            -invM*K -invM*(C(q, q_dot) .+ D)
        ]
        B = [
            zeros(T, 3, 3)
            invM
        ]
    end

    #println("Co = \n", ctrb(A, B))
    #println("rank = ", ctrb(A, B) |> rank)

    if rank(ctrb(A, B)) == size(A, 1)
        println("可制御!")
    else
        println("ランク落ち...")
    end

    P, _, _ = arec(A, B, p.R, p.Q, zero(B))  # リカッチ方程式を解く
    #println(P)
    opt_gain = inv(p.R) * B' * P  # 最適フィードバックゲイン
    println(-opt_gain*x .- G(q))
    return -opt_gain*x .- G(q)
end

"""MPC制御器"""
@with_kw struct MPCController{T}
    Q::Matrix{T}
    R::Matrix{T}
    n::Int64  # 予測ホライズン
    Δt::T  # 予測ホライズンの刻み時価
    isUncertainty::Bool
end

"""線形化したときのA行列"""
function calc_A!(
    q::Vector{Float64}, q_dot::Vector{Float64}, H::Vector{Float64},
    out::Matrix{Float64}
    )
    ccall(
        (:A, "o/soft_robot/derived/mac2/eqs/c_so/A.so"),
        Cvoid,
        (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
        H[1], H[2], H[3], q[1], q_dot[1], q[2], q_dot[2], q[3], q_dot[3], out
    )
end

"""線形化したときのA行列"""
function calc_A(
    q::Vector{Float64}, q_dot::Vector{Float64}, H::Vector{Float64},
    )
    Z = Matrix{Float64}(undef, 9, 9)
    calc_A!(q, q_dot, H, Z)
    Z
end

"""状態方程式のドリフト項"""
function fx!(
    q::Vector{Float64}, q_dot::Vector{Float64}, H::Vector{Float64},
    out::Vector{Float64}
    )
    ccall(
        (:fx, "o/soft_robot/derived/ikko_dake/eqs/c_so/fx.so"),
        Cvoid,
        (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
        H[1], H[2], H[3], q[1], q_dot[1], q[2], q_dot[2], q[3], q_dot[3], out
    )
end

"""状態方程式のドリフト項"""
function fx(
    q::Vector{Float64}, q_dot::Vector{Float64}, H::Vector{Float64},
    )
    Z = Vector{Float64}(undef, 9)
    fx!(q, q_dot, H, Z)
    Z
end

function calc_ℱ(A::Matrix{T}, n::Int64) where T
    Z = Matrix{T}(undef, 9*n, 9)

    for i in 1:n
        Z[(i-1)*9+1:i*9, :] = A^i
    end

    Z
end

function calc_𝒢(A::Matrix{T}, B::Matrix{T}, n::Int64) where T
    Z = Matrix{T}(undef, 9*n, 3*n)

    for i in 1:n
        for j in 1:n
            if j > i
                Z[(i-1)*9+1:i*9, (j-1)*3+1:j*3] = zero(B)
            else
                Z[(i-1)*9+1:i*9, (j-1)*3+1:j*3] = A^(i-1-j) * B
            end
        end
    end

    Z
end

function calc_𝒮(A::Matrix{T}, n::Int64) where T
    Z = Matrix{T}(undef, 9*n, 9*n)

    for i in 1:n
        for j in 1:n
            if j > i
                Z[(i-1)*9+1:i*9, (j-1)*9+1:j*9] = zero(A)
            else
                Z[(i-1)*9+1:i*9, (j-1)*9+1:j*9] = A^(i-1-j)
            end
        end
    end

    Z
end

function calc_ℋ(C::Matrix{T}, n::Int64) where T
    m = size(C, 1)
    Z = zeros(T, m*n, m*n)
    for i in 1:n
        Z[(i-1)*m+1:i*m, (i-1)*m+1:i*m] = C
    end

    Z
end

"""MPCで入力を計算"""
function calc_torque(
    p::MPCController{T},
    q::Vector{T}, q_dot::Vector{T}, H::Vector{T},
    calc_qd, calc_qd_dot,
    t::T
    ) where T

    println(" ")
    println("t = ", t)
    X₀ = [q; q_dot; H]

    # A, B，C行列を計算
    A = calc_A(q, q_dot, H)
    #_fx(q) = fx(q, q_dot, H)
    #A = ForwardDiff.jacobian(_fx, q)
    println(eigvals(A))
    B = [
        zeros(T, 3, 3)
        Matrix{T}(I, 3, 3)
        zeros(T, 3, 3)
    ]
    C = Matrix{T}(I, 9, 9)

    # 目標状態ベクトルと目標入力ベクトル作成
    Yref = Vector{T}(undef, 9*p.n)
    for i in 1:p.n
        Yref[(i-1)*9+1:i*9] = [
            calc_qd(t + i*p.Δt)
            calc_qd_dot(t + i*p.Δt)
            zeros(T, 3)
        ]
    end
    Uref = zeros(T, 3*p.n)

    # 外乱ベクトル作成
    W = Vector{T}(undef, 9*p.n)
    for i in 1:p.n
        W[(i-1)*9+1:i*9] = [
            zeros(T, 6)
            H
        ] .+ fx(q, q_dot, H)
    end

    # F, G, S行列作成
    ℱ = calc_ℱ(A, p.n)
    𝒢 = calc_𝒢(A, B, p.n)
    𝒮 = calc_𝒮(A, p.n)
    ℋ = calc_ℋ(C, p.n)

    # 重み行列を作成
    𝒬 = zeros(T, 9*p.n, 9*p.n)
    for i in 1:p.n
        𝒬[(i-1)*9+1:i*9, (i-1)*9+1:i*9] = p.Q
    end

    ℛ = zeros(T, 3*p.n, 3*p.n)
    for i in 1:p.n
        ℛ[(i-1)*3+1:i*3, (i-1)*3+1:i*3] = p.R
    end

    # println("G ", size(𝒢))
    # println("ℋ ", size(ℋ))
    # println("Q ", size(𝒬))
    # println("R ", size(ℛ))

    #println("G = ", 𝒢)
    println("HGのランク", rank(ℋ * 𝒢))
    ℳ = 𝒢' * ℋ' * 𝒬 * ℋ * 𝒢 .+ ℛ
    𝒩 = (ℋ*(ℱ*X₀ .+ 𝒮*W .- Yref))' * 𝒬 * ℋ * 𝒢 .- Uref'*ℛ
    # println("M, ", size(ℳ))
    # println("N, ", size(𝒩))
    println("Mの固有値 ", eigvals(ℳ))
    Uopt = -inv(ℳ) * 𝒩'  # 最適入力

    # 最適入力を最適トルクに変換
    print("tau = ", inv(M(q)) * Uopt[1:3])
    return inv(M(q)) * Uopt[1:3]
end

end

using LinearAlgebra
using ForwardDiff
using Zygote
using .Controller
A = Controller.calc_A(
    [0.001, 0.002, 0.005],zeros(Float64, 3),zeros(Float64, 3)
)
eigvals(A)
# x = zeros(Float64, 3)
# function _f(q)
#     q_dot = zeros(Float64, 3)
#     H = zeros(Float64, 3)
#     return Controller.fx([q, 0.0, 0.0], q_dot, H)[1]
# end
#_f(q) = [q[1], q[2]^2, q[3]]
#B = _f(zeros(Float64, 9))
#A = ForwardDiff.jacobian(_f, x)
#A = ForwardDiff.derivative(_f, 0.0)

## Gあってるか確認
# using LinearAlgebra
# using .Controller
# A = Controller.calc_A(
#     [0.001, 0.0002, 0.005],zeros(Float64, 3),zeros(Float64, 3)
# )
# B = [
#     zeros(Float64, 3, 3)
#     Matrix{Float64}(I, 3, 3)
#     zeros(Float64, 3, 3)
# ]
# G = Controller.calc_𝒢(A, B, 2)


text:
#=
An XOR linked list is a more memory efficient doubly linked list. Instead of each node holding next and prev fields, it holds a field named both, which is an XOR of the next node and the previous node. Implement an XOR linked list; it has an add(element) which adds the element to the end, and a get(index) which returns the node at index.

If using a language that has no pointers (such as Python), you can assume you have access to get_pointer and dereference_pointer functions that converts between nodes and memory addresses.
=#

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
mutable struct WasmtimeModule <: AbstractWasmModule
    wasmtime_module_ptr::Ptr{wasmtime_module_t}

    WasmtimeModule(ptr::Ptr{wasmtime_module_t}) =
        finalizer(wasmtime_module_delete, new(ptr))
end

function WasmtimeModule(engine, code)
    mod_ptr = Ref{Ptr{wasmtime_module_t}}(Ptr{wasmtime_module_t}())
    @wt_check GC.@preserve mod wasmtime_module_new(
        engine,
        code,
        length(code),
        Base.pointer_from_objref(mod_ptr),
    )
    WasmtimeModule(mod_ptr[])
end

Base.unsafe_convert(::Type{Ptr{wasmtime_module_t}}, mod::WasmtimeModule) =
    mod.wasmtime_module_ptr

Base.show(io::IO, ::WasmtimeModule) = write(io, "WasmtimeModule()")

function imports(mod::WasmtimeModule)
    module_type = wasmtime_module_type(mod)
    mod_imports = WasmPtrVec(wasm_importtype_t)
    wasmtime_moduletype_imports(module_type, mod_imports)
    WasmImports(mod, mod_imports)
end

mutable struct NotInstanciatedWasmExport <: AbstractWasmExport
    wasm_export_ptr::Ptr{wasm_exporttype_t}

    wasm_module::WasmtimeModule
    name::String

    function NotInstanciatedWasmExport(
        wasm_export_ptr::Ptr{wasm_exporttype_t},
        mod::WasmtimeModule,
    )

        owned_wasm_export_ptr = wasm_exporttype_copy(wasm_export_ptr)
        @assert owned_wasm_export_ptr != C_NULL "Failed to copy WASM export"

        name_vec_ptr = wasm_exporttype_name(owned_wasm_export_ptr)
        name_vec = Base.unsafe_load(name_vec_ptr)
        name = unsafe_string(name_vec.data, name_vec.size)
        wasm_name_delete(name_vec_ptr)

        finalizer(new(owned_wasm_export_ptr, mod, name)) do wasm_export
            wasm_exporttype_delete(wasm_export.wasm_export_ptr)
        end
    end
end

function exports(mod::WasmtimeModule)
    module_type = wasmtime_module_type(mod)
    @assert module_type != C_NULL "Failed to get module type from WasmtimeModule"

    wasm_exports = WasmPtrVec(wasm_exporttype_t)
    wasmtime_moduletype_exports(module_type, wasm_exports)

    exports_vector = map(a -> NotInstanciatedWasmExport(a, mod), wasm_exports)

    WasmExports{WasmtimeModule,NotInstanciatedWasmExport}(mod, exports_vector)
end


text:
"""
    NCalculus

Numerical Differentiation and Integration Module.
"""
module NCalculus

const GaussianRoots = [0.5773502692,-0.5773502692,
        0.7745966692,0.0,-0.7745966692,
        0.8611363116,0.3399810436,-0.3399810436,-0.8611363116,
        0.9061798459,0.5384693101,0.0,-0.5384693101,-0.9061798459]
const GaussianCoef = [1.0,1.0,
    0.5555555556,0.8888888889,0.5555555556,
    0.3478548451,0.6521451549,0.6521451549,0.3478548451,
    0.2369268850,0.4786286705,0.5688888889,0.4786286705,0.2369268850]

"""
    ThreePoint(f::Function, x₀::Real, h::Real; method::String="Endpoint")

input a function `f`, x₀ and h.Args is method (default: "Endpoint", if you set method is't "Endpoint", it will
return MidPoint ans.)
"""
@inline function ThreePoint(f::Function, x₀::Real, h::Real; method::String="Endpoint")
    if method == "Endpoint"
        return 1/h*(-3/2*f(x₀) + 2*f(x₀+h)- 1/2*f(x₀+2*h))
    else
        return 1/(2*h) * (f(x₀+h) - f(x₀- h))
    end
end

"""
    FivePoint(f::Function, x₀::Real, h::Real; method::String="Endpoint")

input a function `f`, x₀ and h.Args is method (default: "Endpoint", if you set method is't "Endpoint", it will
return MidPoint ans.)
"""
@inline function FivePoint(f::Function, x₀::Real, h::Real; method::String="Endpoint")
    if method == "Endpoint"
        return 1/(12*h)*(-25*f(x₀)+48f(x₀ +h)-36f(x₀+2h)+16f(x₀+3h)-3f(x₀+4h))
    else
        return 1/(12*h) * (f(x₀ - 2h) - 8f(x₀ -h)+8f(x₀+h)-f(x₀ + 2h))
    end
end

"""
    Trapezoidal(f::function, a::Real, b::Real)

input a function `f`, and ``x in [a, b]`` calculate ∫f(x)dx.
"""
@inline function Trapezoidal(f::Function, a::Real, b::Real)
    return (b-a)/2 * (f(a) + f(b))
end

"""
    Simpson(f::Function, a::Real, b::Real)

input a function `f`, and ``x in [a, b]`` calculate ∫f(x)dx.
"""
@inline function Simpson(f::Function, a::Real, b::Real)
    return (b-a)/6 * (f(a) + 4*f(a+(b-a)/2) + f(b))
end

"""
    Newton_Cotes(f::Function, a::Real, b::Real; n::Int=0, method::String="Midpoint")

input a function `f`, and ``x in [a, b]``` calculate ``int_a^b f(x)dx``. args ``n∈[1,4]``
"""
@inline function Newton_cotes(f::Function, a::Real, b::Real; n::Int=0, method::String="Midpoint")
    if method == "Open"
        o1 = (b-a)/2 * (f(a) + f(b))
        o2 = (b-a)/6 * (f(a) + 4*f(a+(b-a)/2) + f(b))
        o3 = 3*(b-a)/24 * (f(a) + 3*f(a+(b-a)/3) + 3*f(a+2*(b-a)/3) + f(b))
        o4 = 2*(b-a)/180 * (7*f(a) + 32*f(a+(b-a)/4) + 12*f(a+(b-a)/2) + 32*f(a + 3*(b-a)/4) + 7*f(b))
        o = [o1, o2, o3, o4]
        if 0<n≤4
            return o[n]
        else
            return o
        end
    else
        c1 = 2 * (b-a)/2 * f((b+a)/2)
        c2 = 3/2 * (b-a)/3 * (f((b-a)/3) + f(2*(b-a)/3))
        c3 = 4/3 * (b-a)/4 * (2*f((b-a)/4) - f((b-a)/2) + 2*f(3*(b-a)/4))
        c4 = 5/24 * (b-a)/5 * (11*f((b-a)/5)+f(2*(b-a)/5)+f(3*(b-a)/5)+11f(4*(b-a)/5))
        c = [c1, c2, c3, c4]
        if 0<n≤4
            return c[n]
        else
            return c
        end
    end
end

"""
    Romberg(f::Function, a::Real, b::Real, n::Int; seq_tab::Bool=false)
input a function `f`, and ``x in [a, b]``, set `n`.if Args `seq_tab=true`, will return a table, else return a value.
"""
@inline function Romberg(f::Function, a::Real, b::Real, n::Int; seq_tab::Bool=false)
    h = (b-a)
    result = zeros(n,n)
    result[1,1] = (f(a)+f(b))*h/2
    for i = 2:n
        w = f(a)+f(b)
        for j = 1:(2^(i-1)-1)
            w = w + 2*f(a+j*(h/2^(i-1)))
        end
        result[i,1] = w*(h/2^i)
    end
    for i = 2:n
        for j = i:n
            result[j,i] = result[j,i-1]+(1/(4^(i-1)-1))*(result[j,i-1]-result[j-1,i-1])
        end
    end
    if seq_tab
        return result
    else
        return result[n,n]
    end
end

"""
    Gaussian_Quad(f::Function; n::Int=3, a::Real=-1, b::Real=1)

Gaussian Quadrature must set `n<6`,default calculate `int f(x)dx` in `[-1,1]`. you can set interval `[a, b]` to change default interval.
"""
@inline function Gaussian_Quad(f::Function;n::Int=3,a::Real=-1, b::Real=1)
    #Gaussian Args
    #Root = [0.5773502692,-0.5773502692,
    #    0.7745966692,0.0,-0.7745966692,
    #    0.8611363116,0.3399810436,-0.3399810436,-0.8611363116,
    #    0.9061798459,0.5384693101,0.0,-0.5384693101,-0.9061798459]
    #Coef = [1.0,1.0,
    #    0.5555555556,0.8888888889,0.5555555556,
    #    0.3478548451,0.6521451549,0.6521451549,0.3478548451,
    #    0.2369268850,0.4786286705,0.5688888889,0.4786286705,0.2369268850]

    if a == -1 && b == 1
        return GaussianCoef[sum(1:n-1):sum(1:n-1)+n-1]' * f.(GaussianRoots[sum(1:n-1):sum(1:n-1)+n-1])

    else
        return GaussianCoef[sum(1:n-1):sum(1:n-1)+n-1]' * (b-a)/2 * f.((b-a)/2 .* GaussianRoots[sum(1:n-1):sum(1:n-1)+n-1] .+ (b+a)/2)
    end
end

"""
    MutipleIntegralSimple(f::Function, interval_x::Tuple, interval_y::Tuple)

It will apply the **Composite Trapezoidal rule** to approximate the integral.
"""
@inline function MutipleIntegralSimple(f::Function, interval_x::Tuple, interval_y::Tuple)
    a, b = interval_x
    c, d = interval_y

    return (b-a)*(d-c)/16 * (f(a,c) + f(a,d) +f(b,c) + f(b,d)
    + 2*(f((a+b)/2,c) + f((a+b)/2, d) + f(a, (c+d)/2) + f(b, (c+d)/2))
    + 4*f((a+b)/2, (c+d)/2))
end

"""
    SimpsonDoubleIntegral(f::Function, interval_x::Tuple, y_up::Function, y_down::Function; n::Int =6, m::Int=6)

``f  = f(x, y), x ∈ interval_x, y = yup(x) - ydown(x)``, Args m and n is Separation distance.
"""
@inline function SimpsonDoubleIntegral(f::Function, interval_x::Tuple, y_up::Function, y_down::Function; n::Int =6, m::Int=6)
    a, b = interval_x
    h = (b - a) / n
    J₁, J₂, J₃ = zeros(3)
    for i  = 0:n
        x = a + i * h
        HX = (y_up(x) - y_down(x)) / m
        K₁ = f(x, y_down(x)) + f(x, y_up(x))
        K₂, K₃ = 0, 0
        for j = 1:m-1
            y = y_down(x) + j*HX
            Q = f(x, y)
            if iseven(j)
                K₂+= Q
            else
                K₃+= Q
            end
        end
        L = (K₁ + 2*K₂ + 4*K₃)*HX/3
        if i == 0 || i == n
            J₁ = J₁ + L
        elseif iseven(i)
            J₂+= L
        else
            J₃+= L
        end
    end
    return h*(J₁+2J₂+4J₃)/3
end

"""
     GaussianDoubleIntegral(f::Function, interval_x::Tuple, y_up::Function, y_down::Function, m::Int=5, n::Int=5)

``f  = f(x, y), x ∈ interval_x, y = yup(x) - ydown(x)``, Args m and n is Separation distance.
"""
@inline function GaussianDoubleIntegral(f::Function, interval_x::Tuple, y_up::Function, y_down::Function; m::Int=5, n::Int=5)
    m, n = min(5, m), min(5, n)
    a, b = interval_x
    h₁ = (b - a) / 2
    h₂ = (b + a) / 2
    J = 0
    for i in 1:m
        JX = 0
        x = h₁ * GaussianRoots[sum(1:m-1)+i-1] + h₂
        d₁ = y_up(x)
        c₁ = y_down(x)
        k₁ = (d₁ - c₁)/2
        k₂ = (d₁ + c₁)/2
        for j = 1:n
            y = k₁ * GaussianRoots[sum(1:n-1)+j-1] + k₂
            Q = f(x, y)
            JX += GaussianCoef[sum(1:n-1)+j-1] * Q
        end
        J += GaussianCoef[sum(1:m-1)+i-1] * k₁ * JX
    end
    return h₁*J
end

end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 9d996162-d41c-11ea-09b4-994ed99dd014

# ╔═╡ 9d3661fc-d41c-11ea-1997-175b6d9a4095

# ╔═╡ d5eaa180-d41b-11ea-32f6-6ff398c3e4c2
begin
	using LinearAlgebra
	using PlutoUI
	using StatsKit
	import Plots
	import StatsPlots

	md"""
	# Julia 102: Estatística com Pluto.jl

	Vamos explorar um pouco as distribuições e funções estatísticas no Julia. A primeira coisa a notar é que o Julia básico contém o [pacote](https://docs.julialang.org/en/v1/stdlib/Statistics/) `Statistics`, com operações bem simples como `mean()` e `cov()`, e o [pacote](https://docs.julialang.org/en/v1/stdlib/Random/) `Random`, para a geração de números (pseudo)aleatórios. No entanto, há todo um ambiente de computação estatística em desenvolvimento para a linguagem. Vamos utilizar o [meta-pacote](https://github.com/JuliaStats/StatsKit.jl) `StatsKit` que importa alguns dos principais pacotes estatísticos.

	Inicialmente, vamos explorar algumas funções e visualizações de interesse para a Estatística por meio das distribuicões univariadas do pacote [Distributions](https://juliastats.org/Distributions.jl/stable/univariate/).
	"""
end

# ╔═╡ 61be6a20-d41c-11ea-0eee-654e9ec43f39
begin
	md"""
	Vamos brincar um pouco com uma distribuição Normal na qual escolheremos os parâmetros μ e σ², comparando o seu gráfico com o da Normal Padrão.

	**Selecione o parâmetro μ da Normal**:

	$(@bind μ Slider(-50.0:0.1:50.0, 0))

	**Selecione o parâmetro σ² da Normal**:

	$(@bind σ² Slider(1:0.01:50.0))

	Acima podemos conferir o resultado, em que a Normal Padrão está em azul escuro e a nossa Normal está em vermelho. Note que o valor inicial de μ é `0` e o σ² é `1`, logo as curvas estão sobrepostas.
	"""
end

# ╔═╡ 517a652c-d41c-11ea-2118-fdbf45828452
md"""### Distribuição Normal
**Valor de μ**: $μ

**Valor de σ²**: $σ²
"""

# ╔═╡ 57bc858e-d41c-11ea-392b-f1c840845820
#Visualização das distribuições
begin
	Plots.plot(Normal(μ, σ²), color = :red, lw = 2, legend = false, fill = (0, :red), fillalpha = 0.5)
	Plots.plot!(Normal(0,1), color = :royalblue, lw = 2, legend = false, fill = (0, :royalblue), fillalpha = 0.5)
end

# ╔═╡ 664c20be-d41c-11ea-1467-e1c670479af6
begin

	x = rand(Normal(μ, σ²), 100)

	md"""
	Vamos agora gerar uma amostra aleatória de tamanho 100 da distribuição da nossa Normal(μ, σ²) usando o pacote `Random`. O comando está abaixo.

	"""
end

# ╔═╡ acc667f8-d59c-11ea-0805-41b21640540b
begin
	md"""
	Agora vejamos a distribuição LogNormal.

	**Selecione o parâmetro μₗ da LogNormal**:

	$(@bind μₗ Slider(0.0:0.1:25.0))

	**Selecione o parâmetro σ²ₗ da LogNormal**:

	$(@bind σ²ₗ Slider(0.01:0.01:5.0))

	Note que o valor inicial de μₗ é `1.0` e o σ²ₗ é `0.01`.
	"""
end

# ╔═╡ a3dc9324-d59c-11ea-3787-39b6dc720a4b
md"""### Distribuição LogNormal
**Valor de μₗ**: $μₗ

**Valor de σ²ₗ**: $σ²ₗ
"""

# ╔═╡ a651c700-d59c-11ea-029c-1d55b204a41c
#Visualização da distribuição
Plots.plot(LogNormal(μₗ, sqrt(σ²ₗ)), color = :red, lw = 2, legend = false, fill = (0, :red), fillalpha = 0.5)

# ╔═╡ 6ebf2016-d41c-11ea-1a8d-6566e98d1f34
begin
	md"""
	Agora vejamos a distribuição Gamma e o seu formato de acordo com os parâmetros k e θ, onde θ é o parâmetro de escala. Note que nesta parametrização, temos que a função de densidade é dada por

	\begin{equation}
	f(x| k, \theta) = \frac{x^{k-1} e^{-x/θ}}{\Gamma(k) \theta^{k}} \text{ ,}
	\end{equation}

	para $x > 0$, $k > 0$ e $\theta > 0$.

	**Selecione o parâmetro k da Gamma**:

	$(@bind k Slider(1.0:0.1:25.0))

	**Selecione o parâmetro θ da Gamma**:

	$(@bind θ Slider(0.1:0.1:25.0))

	Novamente podemos conferir o resultado. Note que o valor inicial de k é `1.0` e o θ é `0.1`.
	"""
end

# ╔═╡ 6a7b3bf2-d41c-11ea-014e-8b94632190ff
md"""### Distribuição Gamma
**Valor de k**: $k

**Valor de θ**: $θ
"""

# ╔═╡ 6e2923fe-d41c-11ea-03d2-0956b986677c
#Visualização da distribuição
Plots.plot(Gamma(k, θ), color = :red, lw = 2, legend = false, fill = (0, :red), fillalpha = 0.5)

# ╔═╡ 801c622e-d41c-11ea-0441-e300c24904ae
begin
	md"""
	Agora a distribuição mais versátil de todas: a distribuição Beta e o seu formato de acordo com os parâmetros α e β.

	**Selecione o parâmetro α da Beta**:

	$(@bind α Slider(0.1:0.1:10.0))

	**Selecione o parâmetro β da Beta**:

	$(@bind β Slider(0.1:0.1:10.0))

	Os valores iniciais de α e de Β são ambos `0.1`.
	"""
end

# ╔═╡ 6f79fbca-d41c-11ea-0c49-3ffd9d1d5ab2
md"""### Distribuição Beta
**Valor de α**: $α

**Valor de β**: $β
"""

# ╔═╡ 80571608-d41c-11ea-2e72-b7f20192a991
#Visualização da distribuição
Plots.plot(Beta(α, β), color = :red, lw = 2, legend = false, fill = (0, :red), fillalpha = 0.5)

# ╔═╡ 70188128-d41c-11ea-1e8f-5371e257a3bf
begin
	md"""
	Vemos agora a distribuição que ninguém se lembra como é a função de densidade: a distribuição Χ² e o seu formato de acordo com o parâmetro ν de graus de liberdade.

	**Selecione o parâmetro ν da Chi-Quadrado**:

	$(@bind ν Slider(1:50))

	O valor inicial de ν é `1` e só assume valores inteiros positivos.

	"""
end

# ╔═╡ 7fd8c4ce-d41c-11ea-0790-f7afc7a1a84d
#Verificando o valor do parâmetro
md"""
### Distribuição Chi-Quadrado
**Valor de ν**: $ν
"""

# ╔═╡ 7f661988-d41c-11ea-3167-753a433c34e3
#Visualização da distribuição
Plots.plot(Chisq(ν), color = :red, lw = 2, legend = false, fill = (0, :red), fillalpha = 0.5)

# ╔═╡ 9dcde59a-d41c-11ea-36ec-db54e9680172
begin
	md"""
	Como **com certeza** alguém disse na anterior que sabe sim a forma da função de densidade da Χ², então agora sim a distribuição cuja função de densidade que **ninguém** se lembra. Algum desafiante?

	**Selecione o parâmetro d₁ da F-Snedecor**:

	$(@bind d₁ Slider(1:50))

	**Selecione o parâmetro d₂ da F-Snedecor**:

	$(@bind d₂ Slider(1:50))

	Os valores iniciais de d₁ e d₂ são ambos `1`, e só assumem valores inteiros positivos.
	"""
end

# ╔═╡ 9bcc9a70-d41c-11ea-3235-0ffa6ccc0bd5
#Verificando o valor dos parâmetros
md"""
### Distribuição F-Snedecor
**Valor de d₁**: $d₁

**Valor de d₂**: $d₂
"""

# ╔═╡ 9dfe153a-d41c-11ea-2a80-79631819f1ac
#Visualização da distribuição
Plots.plot(FDist(d₁, d₂), color = :red, lw = 2, legend = false, fill = (0, :red), fillalpha = 0.5)

# ╔═╡ Cell order:
# ╠═d5eaa180-d41b-11ea-32f6-6ff398c3e4c2
# ╟─517a652c-d41c-11ea-2118-fdbf45828452
# ╟─57bc858e-d41c-11ea-392b-f1c840845820
# ╟─61be6a20-d41c-11ea-0eee-654e9ec43f39
# ╟─664c20be-d41c-11ea-1467-e1c670479af6
# ╟─a3dc9324-d59c-11ea-3787-39b6dc720a4b
# ╟─a651c700-d59c-11ea-029c-1d55b204a41c
# ╟─acc667f8-d59c-11ea-0805-41b21640540b
# ╟─6a7b3bf2-d41c-11ea-014e-8b94632190ff
# ╟─6e2923fe-d41c-11ea-03d2-0956b986677c
# ╟─6ebf2016-d41c-11ea-1a8d-6566e98d1f34
# ╟─6f79fbca-d41c-11ea-0c49-3ffd9d1d5ab2
# ╟─80571608-d41c-11ea-2e72-b7f20192a991
# ╟─801c622e-d41c-11ea-0441-e300c24904ae
# ╟─7fd8c4ce-d41c-11ea-0790-f7afc7a1a84d
# ╠═7f661988-d41c-11ea-3167-753a433c34e3
# ╟─70188128-d41c-11ea-1e8f-5371e257a3bf
# ╟─9bcc9a70-d41c-11ea-3235-0ffa6ccc0bd5
# ╟─9dfe153a-d41c-11ea-2a80-79631819f1ac
# ╟─9dcde59a-d41c-11ea-36ec-db54e9680172
# ╠═9d996162-d41c-11ea-09b4-994ed99dd014
# ╠═9d3661fc-d41c-11ea-1997-175b6d9a4095

max_stars_repo_name: fausto-mpj/EduPluto

max_issues_repo_name: fausto-mpj/EduPluto

max_forks_repo_name: fausto-mpj/EduPluto


text:
using Hiccup
using Base64

render(e::Editor, ::Nothing) = render(e, icon("check"))

render(::Console, ::Nothing) = nothing

render(::Inline, x::Union{Float16, Float32, Float64}) =
  isnan(x) || isinf(x) ?
    view(span(".syntax--constant.syntax--numeric", string(x))) :
    Dict(:type => :number, :value => float(x), :full => string(x))

@render Inline x::Expr begin
  text = string(x)
  length(split(text, "\n")) == 1 ?
    Model(Dict(:type => :code, :text => text)) :
    Tree(Text("Code"),
         [Model(Dict(:type => :code, :text => text, :attrs => Dict(:block => true)))])
end

render(::Console, x::Expr) =
  Atom.msg("result", Dict(:type => :code, :text => string(x), :attrs => Dict(:block => true)))

@render Inline x::Text begin
  ls = split(chomp(string(x)), "\n")
  length(ls) > 1 ?
    Tree(Model(ls[1]), c(Model(join(ls[2:end], "\n")))) :
    span(ls[1])
end

showmethod(args...) = which(show, (IO, args...))

const inline_mime = "application/prs.juno.inline"

@render Inline x begin
  fields = fieldnames(typeof(x))

  legacy_inline = false
  # Juno-specific display always takes precedence
  if showable("application/juno+inline", x) && !showable(inline_mime, x)
    legacy_inline = true
    m = which(show, (IO, MIME"application/juno+inline", typeof(x)))
    @warn("""
      The \"application/juno+inline\" MIME type is deprecated. Please use \"$(inline_mime)\" instead.
    """, maxlog=1, _id=:juno_inline_legacy, _file=string(m.file), _line=m.line, _module=m.module)
  end
  if showable(inline_mime, x) || legacy_inline
    io = IOBuffer()
    x′ = show(IOContext(io, :limit => true, :color => true), legacy_inline ? MIME"application/juno+inline"() : inline_mime, x)
    if !(x′ isa Nothing)
      defaultrepr(x′, true)
    else
      Text(String(take!(io)))
    end
  elseif showmethod(MIME"text/plain", typeof(x)) ≠ showmethod(MIME"text/plain", Any)
    Text(filter(isvalid, strlimit(sprint(io -> show(IOContext(io, :limit => true, :color => true), MIME"text/plain"(), x)), 2000)))
  elseif showmethod(typeof(x)) ≠ showmethod(Any)
    Text(filter(isvalid, strlimit(sprint(io -> show(IOContext(io, :limit => true, :color => true), x)), 2000)))
  else
    defaultrepr(x, true)
  end
end

function defaultrepr(x, lazy = false)
  fields = fieldnames(typeof(x))
  if isempty(fields)
    span(c(render(Inline(), typeof(x)), "()"))
  else
    lazy ? LazyTree(typeof(x), () -> [SubTree(Text("$f → "), getfield′(x, f, UNDEF)) for f in fields]) :
           Tree(typeof(x), [SubTree(Text("$f → "), getfield′(x, f, UNDEF)) for f in fields])
  end
end

typ(x) = span(".syntax--support.syntax--type", x)

@render Inline x::Type typ(string(x))

for A in :[Vector, Matrix, AbstractVector, AbstractMatrix].args
  @eval begin
    render(i::Inline, ::Type{$A}) =
      render(i, typ($(string(A))))
    (render(i::Inline, ::Type{$A{T}}) where T) =
      render(i, typ(string($(string(A)), "{$T}")))
  end
end

@render Inline x::Module span(".syntax--keyword.syntax--other", string(x))

@render Inline x::Symbol span(".syntax--constant.syntax--other.syntax--symbol", repr(x))

@render Inline x::Char span(".syntax--string.syntax--quoted.syntax--single", escape_string("'$x'"))

@render Inline x::VersionNumber span(".syntax--string.syntax--quoted.syntax--other", sprint(show, x))

@render Inline _::Nothing span(".syntax--constant", "nothing")

@render Inline _::Missing span(".syntax--constant", "missing")

@render Inline _::Undefined span(".fade", "<undefined>")

import Base.Docs: doc

isanon(f) = startswith(string(typeof(f).name.mt.name), "#")

@render Inline f::Function begin
  name = repr(f)
  name_without_module = string(f)
  binding = Docs.Binding(typeof(f).name.module, Symbol(name_without_module))
  if isanon(f)
    span(".syntax--support.syntax--function", "λ")
  else
    LazyTree(
      span(".syntax--support.syntax--function", name),
      ()->[(CodeTools.hasdoc(binding) ? [md_hlines(doc(binding))] : [])..., methods(f)]
    )
  end
end

@render Inline f::Core.IntrinsicFunction begin
  id = Core.Intrinsics.bitcast(Int32, f)
  span(c(span(".syntax--support.syntax--function", string(f)), " (intrinsic function #$id)"))
end

@render Inline f::Core.Builtin begin
  LazyTree(
    span(c(span(".syntax--support.syntax--function", string(f)), " (built-in function)")),
    ()->[(CodeTools.hasdoc(f) ? [md_hlines(doc(f))] : [])..., methods(f)]
  )
end

# TODO: lazy load a recursive tree
function trim(xs, len = 25)
  if length(xs) ≤ len
    undefs(xs)
  else
    [undefs(xs[1:(len÷2 - iseven(len))]); fade("..."); undefs(xs[end-(len÷2 - 1):end])]
  end
end

pluralize(n::Integer, one, more=one) = string(n, " ", n == 1 ? one : more)
pluralize(xs, one, more=one) = pluralize(length(xs), one, more)

@render i::Inline xs::Vector begin
    LazyTree(span(c(render(i, typeof(xs)), Atom.fade(" with $(pluralize(xs, "element", "elements"))"))), () -> trim(xs))
end

@render i::Inline xs::Set begin
    LazyTree(span(c(render(i, typeof(xs)), Atom.fade(" with $(pluralize(xs, "element", "elements"))"))), () -> trim(collect(xs)))
end

@render Inline xs::AbstractArray begin
  Text(sprint(io -> show(IOContext(io, :limit=>true), MIME"text/plain"(), xs)))
end

@render i::Inline d::AbstractDict begin
  cs = () -> begin
    j = 0
    st = Array{Atom.SubTree}(undef, 0)
    for (key, val) in d
      push!(st, SubTree(span(c(render(i, key), " => ")), val))
      j += 1
      j > 25 && (push!(st, SubTree(span("... => "), span("..."))); break)
    end
    return st
  end
  LazyTree(span(c(typ(string(nameof(typeof(d)),
            "{$(eltype(d).parameters[1]), $(eltype(d).parameters[2])}")), Atom.fade(" with $(pluralize(d, "entry", "entries"))"))), cs)
end

@render Inline x::Number span(".syntax--constant.syntax--numeric", sprint(show, x))

@render i::Inline x::Complex begin
  re, ima = reim(x)
  if signbit(ima)
    span(c(render(i, re), " - ", render(i, -ima), "im"))
  else
    span(c(render(i, re), " + ", render(i, ima), "im"))
  end
end

@render Inline p::Ptr begin
  Row(Atom.fade(string(typeof(p))), Text(" @"),
       span(".syntax--constant.syntax--numeric", c("0x$(string(UInt(p), base=16, pad=Sys.WORD_SIZE>>2))")))
end

# TODO: lazy load the rest of the string
@render i::Inline x::AbstractString begin
  x = collect(x)
  length(x) ≤ 500 ?
    span(".syntax--string", c(render(i, Text(stringmime("text/plain", join(x)))))) :
    Row(span(".syntax--string", c("\"", render(i, Text(escape_string(join(x[1:min(length(x),500)])))))),
        Text("..."))
end

render(i::Inline, x::Irrational{sym}) where sym =
  render(i, span(c(string(sym), " = ", render(i, float(x)), "...")))

@render i::Inline xs::Tuple begin
  if length(xs) === 1
    span(c("(", render(i, first(xs)), ",)"))
  else
    span(c("(", interpose(map(x->render(i, x), xs), ", ")..., ")"))
  end
end

include("methods.jl")
include("markdown.jl")


text:
import PyPlot.spy,PyPlot.xticks,PyPlot.yticks
function spy(A::Assoc)
   im = spy(convert(Array,getadj(logical(A))))
   X = getcol(A)
   Y = getrow(A)

   numIdx = 8

   if length(X) > numIdx
       xIdx = 0:convert(Int64,floor((length(X)-1)/numIdx)):length(X)-1
   else
       xIdx = 0:(length(X)-1)
   end
   if length(Y) > numIdx
       yIdx = 0:convert(Int64,floor((length(Y)-1)/numIdx)):length(Y)-1
   else
       yIdx = 0:(length(Y)-1)
   end
   xticks(xIdx, X[collect(xIdx) .+ 1],rotation="vertical")
   yticks(yIdx, Y[collect(yIdx) .+ 1])

   return im

end


text:
# This file is part of the HORIZONS.jl package
# The HORIZONS.jl package is licensed under the MIT "Expat" License
# Copyright (c) 2017: Jorge Pérez.

#The following methods implement some functionality from JPL's Horizons telnet interface
#For more detailed info about JPL HORIZONS system, visit http://ssd.jpl.nasa.gov/?horizons

# The following code is based on the smb_spk script, as was retrieved from:
# ftp://ssd.jpl.nasa.gov/pub/ssd/SCRIPTS/smb_spk
# Date retrieved: Jul 19, 2017
# Credit: Jon D. Giorgini, NASA-JPL
# Jon.D.Giorgini@jpl.nasa.gov

function smb_spk(flag::String, small_body::ObjectName, start_time::StartStopTime,
        stop_time::StartStopTime, email::String="joe@your.domain.name",
        file_name::String=""; kwargs...)

    smb_spk(flag, small_body, DateTime(start_time), DateTime(stop_time), email, file_name; kwargs...)
end

function smb_spk(flag::String, small_body::ObjectName,
        start_time::DateTime, stop_time::DateTime,
        email::String="joe@your.domain.name", file_name::String="";
        ftp_verbose::Bool=false, timeout::Int=60)

    # TODO: handle spk_ID_override (see smb_spk Expect script)

    # Convert start and stop time from `DateTime`s to `String`s
    small_body_str = string(small_body)
    start_time_str = Dates.format(start_time, HORIZONS_DATE_FORMAT)
    stop_time_str = Dates.format(stop_time, HORIZONS_DATE_FORMAT)

    if flag == "t"
        file_type = "yes"
        ftp_type = :ascii
        ftp_sufx = ".xsp"
    elseif flag == "b"
        file_type  = "no"
        ftp_type = :binary
        ftp_sufx = ".bsp"
    elseif flag == "1"
        file_type = "1"
        ftp_type = :binary
        ftp_sufx = ".bsp"
    elseif flag == "2"
        file_type = "21"
        ftp_type = :binary
        ftp_sufx = ".bsp"
    else
        ArgumentError("Unknown file type: $flag")
    end

    local_file = file_name

    # Connect to Horizons
    proc = ExpectProc(`telnet $HORIZONS_MACHINE 6775`, timeout)

    # Wait for main Horizons prompt, set up, proceed
    idx = expect!(proc, ["unknown host", "Horizons> "])
    if idx == 1
        warn("This system cannot find $HORIZONS_MACHINE")
        close(proc)
    elseif idx == 2
        println(proc, "PAGE")
    end

    idx = expect!(proc, ["Horizons> "], timeout=15)
    if idx == 1
        println(proc, "##2")
    end
    idx = expect!(proc, ["Horizons> "], timeout=15)
    if idx == 1
        println(proc, small_body_str)
    end

    # Handle prompt search/select
    idx = expect!(proc, [r".*Continue.*: $", r".*such object record.*", r".*Select.*<cr>: $"])
    if idx == 1
        println(proc, "yes")
        idx = expect!(proc, [r".*PK.*: $", r".*lay.*: $"])
        if idx ==1
            # TODO: handle spk_ID_override (see smb_spk Expect script)
            println(proc, "S")
        elseif idx == 2
            println(proc, "x")
            throw(ArgumentError("Cancelled -- unique object not found: $small_body_str\nObject not matched to database OR multiple matches found."))
        end
    elseif idx == 2
        # currently unable to reproduce this case on HORIZONS v4.10
        println(proc, "x")
        throw(ArgumentError("No such object record found: $small_body_str"))
    elseif idx == 3
        println(proc, "S")
    end

    # Pick out SPK ID sent by Horizons
    if file_name == ""
        idx = expect!(proc, [r" Assigned SPK object ID:  (.*)\r\r\n \r\r\n Enter your"])
        proc_match = match(r" Assigned SPK object ID:  (.*)\r\r\n \r\r\n Enter your", proc.match)
        # @show proc_match
        spkid = proc_match[1]
        local_file = spkid*ftp_sufx
    end

    # Process prompt for e-mail address
    idx = expect!(proc, [r".*address.*: $"])
    if idx == 1
        println(proc, email)
    end

    # Process e-mail confirmation
    idx = expect!(proc, [r".*yes.*: $"])
    if idx == 1
        println(proc, "yes")
    end

    # Process file type
    idx = expect!(proc, [r".*SPK file format.*: $", r".*YES.*: $"])
    if idx == 1
        println(proc, file_type)
    elseif idx == 2
        println(proc, file_type)
    end

    # Set start date
    idx = expect!(proc, [r".*START.*: "s])
    if idx == 1
        println(proc, start_time_str)
    end

    # Handle start date error or STOP date
    idx = expect!(proc, [r".*try.*: "s, r".*STOP.*: "s])
    if idx == 1
        println(proc, "X")
        throw(ArgumentError("START time $start_time_str outside set SPK limits."))
    elseif idx == 2
        println(proc, stop_time_str)
    end

    # Handle stop date error
    # timeout is set to Inf to allow for file production
    idx = expect!(proc, [r".*large.*", r".*try.*"s, r".*time-span too small.*", r".*YES.*"s], timeout=Inf)
    if idx == 1
        println(proc, "X")
        throw(ArgumentError("Error in STOP date: $stop_time_str\nStop date must not be more than 200 years after start."))
    elseif idx == 2
        println(proc, "X")
        throw(ArgumentError("STOP time $stop_time_str outside set SPK limits."))
    elseif idx == 3
        println(proc, "X")
        throw(ArgumentError("Error in requested length: $start_time_str to $stop_time_str\nTime span of file must be >= 32 days."))
    elseif idx == 4
        # Binary SPK file created.
        # Add more objects to file  [ YES, NO, ? ] :
        println(proc, "NO")
    end

    # Pick out FTP file name
    idx = expect!(proc, [r"File name   : (.*)\r\r\n   File type"], timeout=15)
    proc_match = match(r"File name   : (.*)\r\r\n   File type", proc.match)
    ftp_name = strip(proc_match[1])

    # Close telnet connection
    close(proc)

    # Retrieve file by anonymous FTP and save to file
    ftp_init()
    # workaround `@` in email address
    ftp_email = replace(email, "@" => "_at_")
    ftp = FTP(hostname=HORIZONS_MACHINE, username="anonymous", password=ftp_email, verbose=ftp_verbose)
    cd(ftp, HORIZONS_FTP_DIR)
    if file_name == ""
        io = download(ftp, ftp_name, local_file)
        close(ftp)
        ftp_cleanup()
        return ftp_name, local_file
    else
        io = download(ftp, ftp_name, file_name)
        close(ftp)
        ftp_cleanup()
        return ftp_name, file_name
    end
end