text:
# C- entrypoints for all lowlevel callback functions

export fuse_session_new, fuse_session_destroy, fuse_session_mount, fuse_session_unmount, fuse_session_loop
export fuse_reply_attr, fuse_reply_bmap, fuse_reply_buf, fuse_reply_create, fuse_reply_data, fuse_reply_entry, fuse_reply_err
export fuse_reply_ioctl, fuse_reply_ioctl_iov, fuse_reply_ioctl_retry, fuse_reply_iov, fuse_reply_lock, fuse_reply_lseek, fuse_reply_none
export fuse_reply_attrfuse_reply_none, fuse_reply_open, fuse_reply_poll, fuse_reply_readlink, fuse_reply_statfs, fuse_reply_write, fuse_reply_xattr
export fuse_add_direntry, fuse_add_direntry_plus, fuse_req_ctx

export fuse_req_ctx, fuse_req_getgroups, fuse_req_interrupt_func, fuse_req_interrupted, fuse_req_userdata

# F_INIT = 1 (26)
function Ginit(init::Function, fs::Any)
    Cinit = let init = init, fs = fs
        function (::Ptr{Nothing}, conn::Ptr{FuseConnInfo})
            docall() do
                init(fs, CStruct{FuseConnInfo}(conn))
            end
        end
    end
    @cfunction $Cinit Cvoid (Ptr{Nothing}, Ptr{FuseConnInfo})
end
# F_DESTROY = 2
function Gdestroy(destroy::Function, fs::Any)
    function Cdestroy(::Ptr{Nothing})
        docall() do
            destroy(fs)
        end
    end
    (@cfunction $Cdestroy Cvoid (Ptr{Nothing},))
end
# F_LOOKUP = 3 (1)
function Glookup(lookup::Function, fs::Any)
    function Clookup(req::FuseReq, parent::FuseIno, name::Cstring)
        docall(req) do
            name = unsafe_string(name)
            error = lookup(fs, req, parent, name)
            error
        end
        nothing
    end
    (@cfunction $Clookup  Cvoid (FuseReq, FuseIno, Cstring))
end
# F_FORGET = 4
function Gforget(forget::Function, fs::Any)
    function Cforget(req::FuseReq, ino::FuseIno, lookup::Cuint64_t)
        docall(req) do
            forget(fs, req, ino, lookup)
        end
    end
    (@cfunction $Cforget Cvoid (FuseReq, FuseIno, Cuint64_t))
end
# F_GETATTR = 5
function Ggetattr(getattr::Function, fs::Any)
    function Cgetattr(req::FuseReq, ino::FuseIno, fi::Ptr{FuseFileInfo})
        docall(req) do
            getattr(fs, req, ino, CStruct(fi))
        end
        nothing
    end
    (@cfunction $Cgetattr Cvoid (FuseReq, FuseIno, Ptr{FuseFileInfo}))
end
# F_SETATTR = 6
function Gsetattr(setattr::Function, fs::Any)
    function Csetattr(req::FuseReq, ino::FuseIno, attr::Ptr{Cstat}, to_set::Cint, fi::Ptr{FuseFileInfo})
        docall(req) do
             setattr(fs, req, ino, CStruct(attr), to_set, CStruct(fi))
        end
    end
    (@cfunction $Csetattr Cvoid (FuseReq, FuseIno, Ptr{Cstat}, Cint, Ptr{FuseFileInfo}))
end
# F_READLINK = 7
function Greadlink(readlink::Function, fs::Any)
    function Creadlink(req::FuseReq, ino::FuseIno)
        docall(req) do
            readlink(fs, req, ino)
        end
    end
    (@cfunction $Creadlink Cvoid (FuseReq, FuseIno))
end
# F_MKNOD = 8
function Gmknod(mknod::Function, fs::Any)
    function Cmknod(req::FuseReq, parent::FuseIno, name::Cstring, mode::FuseMode, rdev::FuseDev)
        docall(req) do
            name = unsafe_string(name)
            mknod(fs, req, parent, name, mode, rdev)
        end
    end
    (@cfunction $Cmknod Cvoid (FuseReq, FuseIno, Cstring, FuseMode, FuseDev))
end
# F_MKDIR = 9
function Gmkdir(mkdir::Function, fs::Any)
    function Cmkdir(req::FuseReq, parent::FuseIno, name::Cstring, mode::FuseMode)
        docall(req) do
            name = unsafe_string(name)
            mkdir(fs, req, parent, name, mode)
        end
    end
    (@cfunction $Cmkdir Cvoid (FuseReq, FuseIno, Cstring, FuseMode))
end
# F_UNLINK = 10
function Gunlink(unlink::Function, fs::Any)
    function Cunlink(req::FuseReq, parent::FuseIno, name::Cstring)
        docall(req) do
            name = unsafe_string(name)
            unlink(fs, req, parent, name)
        end
    end
    (@cfunction $Cunlink Cvoid (FuseReq, FuseIno, Cstring))
end
# F_RMDIR = 11
function Grmdir(rmdir::Function, fs::Any)
    function Crmdir(req::FuseReq, parent::FuseIno, name::Cstring)
        docall(req) do
            name = unsafe_string(name)
            rmdir(fs, req, parent, name)
        end
    end
    (@cfunction $Crmdir Cvoid (FuseReq, FuseIno, Cstring))
end
# F_SYMLINK = 12
function Gsymlink(symlink::Function, fs::Any)
    function Csymlink(req::FuseReq, link::Cstring, parent::FuseIno, name::Cstring)
        docall(req) do
            name = unsafe_string(name)
            link = unsafe_string(link)
            symlink(fs, req, link, parent, name)
        end
    end
    (@cfunction $Csymlink Cvoid (FuseReq, Cstring, FuseIno, Cstring))
end
# F_RENAME = 13
function Grename(rename::Function, fs::Any)
    function Crename(req::FuseReq, parent::FuseIno, name::Cstring, newparent::FuseIno, newname::Cstring, flags::Cuint)
        docall(req) do
            name = unsafe_string(name)
            newname = unsafe_string(newname)
            rename(fs, req, parent, name, newparent, newname, flags)
        end
    end
    (@cfunction $Crename Cvoid (FuseReq, FuseIno, Cstring, FuseIno, Cstring, Cuint))
end
# F_LINK = 14
function Glink(link::Function, fs::Any)
    function Clink(req::FuseReq, ino::FuseIno, newparent::FuseIno, newname::Cstring)
        docall(req) do
            newname = unsafe_string(newname)
            link(fs, req, ino, newparent, newname)
        end
    end
    (@cfunction $Clink Cvoid (FuseReq, FuseIno, FuseIno, Cstring))
end
# F_OPEN = 15
function Gopen(open::Function, fs::Any)
    function Copen(req::FuseReq, ino::FuseIno, fi::Ptr{FuseFileInfo})
        docall(req) do
            open(fs, req, ino, CStruct(fi))
        end
    end
    (@cfunction $Copen Cvoid (FuseReq, FuseIno, Ptr{FuseFileInfo}))
end
# F_READ = 16
function Gread(read::Function, fs::Any)
    function Cread(req::FuseReq, ino::FuseIno, size::Csize_t, off::Coff_t, fi::Ptr{FuseFileInfo})
        docall(req) do
            read(fs, req, ino, size, off, CStruct(fi))
        end
    end
    (@cfunction $Cread Cvoid (FuseReq, FuseIno, Csize_t, Coff_t, Ptr{FuseFileInfo}))
end
# F_WRITE = 17
function Gwrite(write::Function, fs::Any)
    function Cwrite(req::FuseReq, ino::FuseIno, cbuf::Ptr{UInt8}, size::Csize_t, off::Csize_t, fi::Ptr{FuseFileInfo})
        docall(req) do
            write(fs, req, ino, CVector{UInt8}(cbuf), size, off, CStruct(fi))
        end
    end
    (@cfunction $Cwrite Cvoid (FuseReq, FuseIno, Ptr{UInt8}, Csize_t, Csize_t, Ptr{FuseFileInfo}))
end
# F_FLUSH = 18
function Gflush(flush::Function, fs::Any)
    function Cflush(req::FuseReq, ino::FuseIno, fi::Ptr{FuseFileInfo})
        docall(req) do
            flush(fs, req, ino, CStruct(fi))
        end
    end
    (@cfunction $Cflush Cvoid (FuseReq, FuseIno, Ptr{FuseFileInfo}))
end
# F_RELEASE = 19
function Grelease(release::Function, fs::Any)
    function Crelease(req::FuseReq, ino::FuseIno, fi::Ptr{FuseFileInfo})
        docall(req) do
            release(fs, req, ino, CStruct(fi))
        end
    end
    (@cfunction $Crelease Cvoid (FuseReq, FuseIno, Ptr{FuseFileInfo}))
end
# F_FSYNC = 20
function Gfsync(fsync::Function, fs::Any)
    function Cfsync(req::FuseReq, ino::FuseIno, datasync::Cint, fi::Ptr{FuseFileInfo})
        docall(req) do
            fsync(fs, req, ino, datasync, CStruct(fi))
        end
    end
    (@cfunction $Cfsync Cvoid (FuseReq, FuseIno, Cint, Ptr{FuseFileInfo}))
end
# F_OPENDIR = 21 (37)
function Gopendir(opendir::Function, fs::Any)
    function Copendir(req::FuseReq, ino::FuseIno, fi::Ptr{FuseFileInfo})
        docall(req) do
            opendir(fs, req, ino, CStruct(fi))
        end
    end
    (@cfunction $Copendir Cvoid (FuseReq, FuseIno, Ptr{FuseFileInfo}))
end
# F_READDIR = 22
function Greaddir(readdir::Function, fs::Any)
    function Creaddir(req::FuseReq, ino::FuseIno, size::Csize_t, off::Csize_t, fi::Ptr{FuseFileInfo})
        docall(req) do
            readdir(fs, req, ino, size, off, CStruct(fi))
        end
    end
    (@cfunction $Creaddir Cvoid (FuseReq, FuseIno, Csize_t, Csize_t, Ptr{FuseFileInfo}))
end
# F_RELEASEDIR = 23
function Greleasedir(releasedir::Function, fs::Any)
    function Creleasedir(req::FuseReq, ino::FuseIno, fi::Ptr{FuseFileInfo})
        docall(req) do
            releasedir(fs, req, ino, CStruct(fi))
        end
    end
    (@cfunction $Creleasedir Cvoid (FuseReq, FuseIno, Ptr{FuseFileInfo}))
end
# F_FSYNCDIR = 24
function Gfsyncdir(fsyncdir::Function, fs::Any)
    function Cfsyncdir(req::FuseReq, ino::FuseIno, datasync::Cint, fi::Ptr{FuseFileInfo})
        docall(req) do
            fsyncdir(fs, req, ino, datasync, CStruct(fi))
        end
    end
    (@cfunction $Cfsyncdir Cvoid (FuseReq, FuseIno, Cint, Ptr{FuseFileInfo}))
end
# F_STATFS = 25
function Gstatfs(statfs::Function, fs::Any)
    function Cstatfs(req::FuseReq, ino::FuseIno)
        docall(req) do
            statfs(fs, req, ino)
        end
    end
    (@cfunction $Cstatfs Cvoid (FuseReq, FuseIno))
end
# F_SETXATTR = 26
function Gsetxattr(setxattr::Function, fs::Any)
    function Csetxattr(req::FuseReq, ino::FuseIno, name::Cstring, value::Ptr{UInt8}, size::Csize_t, flags::Cint)
        docall(req) do
            name = unsafe_string(name)
            value = copy(unsafe_wrap(Array, value, size))
            setxattr(fs, req, ino, name, value, flags)
        end
    end
    (@cfunction $Csetxattr Cvoid (FuseReq, FuseIno, Cstring, Ptr{UInt8}, Csize_t, Cint))
end
# F_GETXATTR = 27
function Ggetxattr(getxattr::Function, fs::Any)
    function Cgetxattr(req::FuseReq, ino::FuseIno, name::Cstring, size::Csize_t)
        docall(req) do
            name = unsafe_string(name)
            getxattr(fs, req, ino, name, size)
        end
    end
    (@cfunction $Cgetxattr Cvoid (FuseReq, FuseIno, Cstring, Csize_t))
end
# F_LISTXATTR = 28
function Glistxattr(listxattr::Function, fs::Any)
    function Clistxattr(req::FuseReq, ino::FuseIno, size::Csize_t)
        docall(req) do
            listxattr(fs, req, ino, size)
        end
    end
    (@cfunction $Clistxattr Cvoid (FuseReq, FuseIno, Csize_t))
end
# F_REMOVEXATTR = 29
function Gremovexattr(removexattr::Function, fs::Any)
    function Cremovexattr(req::FuseReq, ino::FuseIno, name::Cstring)
        docall(req) do
            name = unsafe_string(name)
            removexattr(fs, req, ino, name)
        end
    end
    (@cfunction $Cremovexattr Cvoid (FuseReq, FuseIno, Cstring))
end
# F_ACCESS = 30 (34)
function Gaccess(access::Function, fs::Any)
    function Caccess(req::FuseReq, ino::FuseIno, mask::Cint)
        docall(req) do
            access(fs, req, ino, mask)
        end
    end
    (@cfunction $Caccess Cvoid (FuseReq, FuseIno, Cint))
end
# F_CREATE = 31
function Gcreate(create::Function, fs::Any)
    function Ccreate(req::FuseReq, parent::FuseIno, name::Cstring, mode::FuseMode, fi::Ptr{FuseFileInfo})
        docall(req) do
            name = unsafe_string(name)
            create(fs, req, parent, name, mode, CStruct(fi))
        end
    end
    (@cfunction $Ccreate Cvoid (FuseReq, FuseIno, Cstring, FuseMode, Ptr{FuseFileInfo}))
end
# F_GETLK = 32
function Ggetlk(getlk::Function, fs::Any)
    function Cgetlk(req::FuseReq, ino::FuseIno, fi::Ptr{FuseFileInfo}, lock::Ptr{Cflock})
        docall(req) do
            getlk(fs, req, ino, CStruct(fi), CStruct(lock))
        end
    end
    (@cfunction $Cgetlk Cvoid (FuseReq, FuseIno, Ptr{FuseFileInfo}, Ptr{Cflock}))
end
# F_SETLK = 33
function Gsetlk(setlk::Function, fs::Any)
    function Csetlk(req::FuseReq, ino::FuseIno, fi::Ptr{FuseFileInfo}, lock::Ptr{Cflock}, sleep::Cint)
        docall(req) do
            setlk(fs, req, ino, CStruct(fi), CStruct(lock), sleep)
        end
    end
    (@cfunction $Csetlk Cvoid (FuseReq, FuseIno, Ptr{FuseFileInfo}, Ptr{Cflock}, Cint))
# F_BMAP = 34
function Gbmap(bmap::Function, fs::Any)
    function Cbmap(req::FuseReq, ino::FuseIno, blocksize::Csize_t, idx::Cuint64_t)
        docall(req) do
            bmap(fs, req, ino, blocksize, idx)
        end
    end
    (@cfunction $Cbmap Cvoid (FuseReq, FuseIno, Csize_t, Cuint64_t))
end
end
# F_IOCTL = 35
function Gioctl(ioctl::Function, fs::Any)
    function Cioctl(req::FuseReq, ino::FuseIno, cmd::Cuint, arg::Ptr{Nothing}, fi::Ptr{FuseFileInfo}, flags::Cuint, in_buf::Ptr{Ptr{Cvoid}}, in_bufsz::Csize_t, out_bufsz::Csize_t)
        docall(req) do
            ioctl(fs, req, ino, cmd, arg, CStruct(fi), flags, in_buf, in_bufsz, out_bufsz)
        end
    end
    (@cfunction $Cioctl Cvoid (FuseReq, FuseIno, Cuint, Ptr{Cvoid}, Ptr{FuseFileInfo}, Cuint, Ptr{Ptr{Cvoid}}, Csize_t, Csize_t))
end
# F_POLL = 36
function Gpoll(poll::Function, fs::Any)
    function Cpoll(req::FuseReq, ino::FuseIno, fi::Ptr{FuseFileInfo}, ph::Ptr{FusePollHandle})
        docall(req) do
            poll(fs, req, ino, CStruct(fi), CStruct(ph))
        end
    end
    (@cfunction $Cpoll Cvoid (FuseReq, FuseIno, Ptr{FuseFileInfo}, Ptr{FusePollHandle}))
end
# F_WRITE_BUF = 37
function Gwrite_buf(write_buf::Function, fs::Any)
    function Cwrite_buf(req::FuseReq, ino::FuseIno, bufv::Ptr{FuseBufvec}, off::Coff_t, fi::Ptr{FuseFileInfo})
        docall(req) do
            write_buf(fs, req, ino, CStruct(bufv), off, CStruct(fi))
        end
    end
    (@cfunction $Cwrite_buf Cvoid (FuseReq, FuseIno, Ptr{FuseBufvec}, Coff_t, Ptr{FuseFileInfo}))
end
# F_RETRIEVE_REPLY = 38
function Gretrieve_reply(retrieve_reply::Function, fs::Any)
    function Cretrieve_reply(req::FuseReq, cookie::Ptr{Cvoid}, ino::FuseIno, off::Coff_t, bufv::Ptr{FuseBufvec})
        docall(req) do
            retrieve_reply(fs, req, cookie, ino, off, CStruct(bufv))
        end
    end
    (@cfunction $Cretrieve_reply Cvoid (FuseReq, Ptr{Cvoid}, FuseIno, Coff_t, Ptr{FuseBufvec}))
end
# F_FORGET_MULTI = 39
function Gforget_multi(forget_multi::Function, fs::Any)
    function Cforget_multi(req::FuseReq, count::Csize_t, forgets::Ptr{FuseForgetData})
        docall(req) do
            forget_multi(fs, req, count, CStruct{forgets})
        end
    end
    (@cfunction $Cforget_multi Cvoid (FuseReq, Csize_t, Ptr{FuseForgetData}))
end
# F_FLOCK = 40
function Gflock(flock::Function, fs::Any)
    function Cflock(req::FuseReq, ino::FuseIno, fi::Ptr{FuseFileInfo}, op::Cint)
        docall(req) do
            flock(fs, req, ino, CStruct(fi), op)
        end
    end
    (@cfunction $Cflock Cvoid (FuseReq, FuseIno, Ptr{FuseFileInfo}, Cint))
end
# F_FALLOCATE = 41
function Gfallocate(fallocate::Function, fs::Any)
    function Cfallocate(req::FuseReq, ino::FuseIno, mode::Cint, off::Coff_t, length::Coff_t, fi::Ptr{FuseFileInfo})
        docall(req) do
            fallocate(fs, req, ino, mode, off, length, CStruct(fi))
        end
    end
    (@cfunction $Cfallocate Cvoid (FuseReq, FuseIno, Cint, Coff_t, Coff_t, Ptr{FuseFileInfo}))
end
# F_READDIRPLUS = 42
function Greaddirplus(readdirplus::Function, fs::Any)
    #, size_t size, off_t off, struct fuse_file_info *fi)
    function Creaddirplus(req::FuseReq, ino::FuseIno, size::Csize_t, off::Coff_t, fi::Ptr{FuseFileInfo})
        docall(req) do
            readdirplus(fs, req, ino, size, off, CStruct(fi))
        end
    end
    (@cfunction $Creaddirplus Cvoid (FuseReq, FuseIno, Csize_t, Coff_t, Ptr{FuseFileInfo}))
end
# F_COPY_FILE_RANGE = 43
function Gcopy_file_range(copy_file_range::Function, fs::Any)
    function Ccopy_file_range(req::FuseReq, ino_in::FuseIno, off_in::Coff_t, fi_in::Ptr{FuseFileInfo}, ino_out::FuseIno, off_out::Coff_t, fi_out::Ptr{FuseFileInfo}, len::Csize_t, flags::Cint)
        docall(req) do
            copy_file_range(fs, req, ino_in, off_in, CStruct(fi_in), ino_out, off_out, CStruct(fi_out), len, flags)
        end
    end
    (@cfunction $Ccopy_file_range Cvoid (FuseReq, FuseIno, Coff_t, Ptr{FuseFileInfo}, FuseIno, Coff_t, Ptr{FuseFileInfo}, Csize_t, Cint))
end
# F_LSEEK = 44
function Glseek(lseek::Function, fs::Any)
    function Clseek(req::FuseReq, ino::FuseIno, off::Coff_t, whence::Cint, fi::Ptr{FuseFileInfo})
        docall(req) do
            lseek(fs, req, ino, off, whence, CStruct(fi))
        end
    end
    (@cfunction $Clseek Cvoid (FuseReq, FuseIno, Coff_t, Cint, Ptr{FuseFileInfo}))
end

# reply functions to be called inside callback functions

function fuse_reply_attr(req::FuseReq, attr::CStructAccess{Cstat}, attr_timeout::Real)
    ccall((:fuse_reply_attr, :libfuse3), Cint, (FuseReq, Ptr{Cstat}, Cdouble), req, attr, attr_timeout)
end
function fuse_reply_bmap(req::FuseReq, idx::Integer)
    ccall((:fuse_reply_bmap, :libfuse3), Cint, (FuseReq, UInt64), req, idx)
end
function fuse_reply_buf(req::FuseReq, buf::Vector{UInt8}, size::Integer)
    ccall((:fuse_reply_buf, :libfuse3), Cint, (FuseReq, Ptr{UInt8}, Csize_t), req, pointer(buf), size)
end
function fuse_reply_create(req::FuseReq, e::CStructAccess{FuseEntryParam}, fi::CStruct{FuseFileInfo})
    ccall((:fuse_reply_create, :libfuse3), Cint, (FuseReq, Ptr{FuseEntryParam}, Ptr{FuseFileInfo}), req, e, fi)
end
function fuse_reply_data(req::FuseReq, bufv::CStructAccess{FuseBufvec}, flags::FuseBufCopyFlags)
    ccall((:fuse_reply_data, :libfuse3), Cint, (FuseReq, Ptr{FuseBufvec}, Cint), req, bufv, flags.flag)
end
function fuse_reply_entry(req::FuseReq, entry::CStructAccess{FuseEntryParam})
    ccall((:fuse_reply_entry, :libfuse3), Cint, (FuseReq, Ptr{FuseEntryParam}), req, entry)
end
function fuse_reply_none(req::FuseReq)
    ccall((:fuse_reply_none, :libfuse3), Cint, (FuseReq,), req)
end
function fuse_reply_err(req::FuseReq, err::Integer)
    ccall((:fuse_reply_err, :libfuse3), Cint, (FuseReq, Cint), req, err)
end
function fuse_reply_ioctl(req::FuseReq, result::Integer, buf::CStruct, size::Integer)
    ccall((:fuse_reply_ioctl, :libfuse3), Cint, (FuseReq, Cint, Ptr{Nothing}, Csize_t), req, result, buf, size)
end
function fuse_reply_ioctl_iov(req::FuseReq, result::Integer, iov::CStructAccess{Ciovec}, count::Integer)
    ccall((:fuse_reply_ioctl_iov, :libfuse3), Cint, (FuseReq, Cint, Ptr{Ciovec}, Cint), req, result, iov, count)
end
function fuse_reply_ioctl_retry(req::FuseReq, in_iov::CStructAccess{Ciovec}, in_count::Integer, out_iov::CStruct{Ciovec}, out_count::Integer)
    ccall((:fuse_reply_ioctl_retry, :libfuse3), Cint, (FuseReq, Ptr{Ciovec}, Csize_t, Ptr{Ciovec}, Csize_t), req, in_iov, in_count, out_iov, out_count)
end
function fuse_reply_iov(req::FuseReq, iov::CStructAccess{Ciovec}, count::Integer)
    ccall((:fuse_reply_iov, :libfuse3), Cint, (FuseReq, Ptr{Ciovec}, Cint), req, iov, count)
end
function fuse_reply_lock(req::FuseReq, lock::CStructAccess{Cflock})
    ccall((:fuse_reply_lock, :libfuse3), Cint, (FuseReq, Ptr{Cflock}), req, lock)
end
function fuse_reply_lseek(req::FuseReq, off::Integer)
    ccall((:fuse_reply_lseek, :libfuse3), Cint, (FuseReq, Coff_t), req, off)
end
function fuse_reply_open(req::FuseReq, fi::CStructAccess{FuseFileInfo})
    ccall((:fuse_reply_open, :libfuse3), Cint, (FuseReq, Ptr{FuseFileInfo}), req, fi)
end
function fuse_reply_poll(req::FuseReq, revents::Integer )
    ccall((:fuse_reply_entry, :libfuse3), Cint, (FuseReq, Cuint), req, revents)
end
function fuse_reply_readlink(req::FuseReq, link::Vector{UInt8})
    ccall((:fuse_reply_readlink, :libfuse3), Cint, (FuseReq, Ptr{UInt8}), req, link)
end
function fuse_reply_statfs(req::FuseReq, stbuf::CStructAccess{Cstatvfs})
    ccall((:fuse_reply_statfs, :libfuse3), Cint, (FuseReq, Ptr{Cstatvfs}), req, stbuf)
end
function fuse_reply_write(req::FuseReq, count::Integer)
    ccall((:fuse_reply_write, :libfuse3), Cint, (FuseReq, Csize_t), req, count)
end
function fuse_reply_xattr(req::FuseReq, count::Integer)
    ccall((:fuse_reply_xattr, :libfuse3), Cint, (FuseReq, Csize_t), req, count)
end

# accessors for req
function fuse_req_getgroups(req::FuseReq, list::Vector{Cgid_t})
    ccall((:fuse_req_getgroups, :libfuse3), Cint, (FuseReq, Cint, Ptr{Cgid_t}), req, length(list), pointer_from_vector(list))
end
function fuse_req_interrupt_func(req::FuseReq, func::Ptr{Nothing}, data::Any)
    ccall((:fuse_req_ctx, :libfuse3), Cvoid, (FuseReq, Ptr{Nothing}, Ptr{Nothing}), req, func, data)
end
function fuse_req_interrupted(req::FuseReq)
    ccall((:fuse_req_interrupted, :libfuse3), Cint, (FuseReq,), req)
end
function fuse_req_userdata(req::FuseReq)
    ptr_to_userdata(ccall((:fuse_req_userdata, :libfuse3), Ptr{Any}, (FuseReq,), req))
end

# helper function to support readdir
function fuse_add_direntry(req::FuseReq, buf::AbstractVector{UInt8}, name::String, st::CStructAccess{Cstat}, off::Integer)
    bufsz = length(buf)
    ccall((:fuse_add_direntry, :libfuse3), Csize_t, (FuseReq, Ptr{UInt8}, Csize_t, Cstring, Ptr{Cstat}, Coff_t), req, buf, bufsz, name, st, off)
end

function fuse_add_direntry_plus(req::FuseReq, buf::AbstractVector{UInt8}, name::String, e::CStructAccess{FuseEntryParam}, off::Integer)
    bufsz = length(buf)
    ccall((:fuse_add_direntry_plus, :libfuse3), Csize_t, (FuseReq, Ptr{UInt8}, Csize_t, Cstring, Ptr{FuseEntryParam}, Coff_t), req, buf, bufsz, name, e, off)
end

function fuse_req_ctx(req::FuseReq)
    cdata = ccall((:fuse_req_ctx, :libfuse3), Ptr{FuseCtx}, (FuseReq,), req)
    CStruct(cdata)
end

# session functions

function fuse_session_new(fargs::CStructAccess{FuseCmdlineArgs}, callbacks::FuseLowlevelOps)
    ccall((:fuse_session_new, :libfuse3), Ptr{Cvoid},
        (Ptr{FuseCmdlineArgs}, Ref{FuseLowlevelOps}, Cint, Ref{Nothing}),
        fargs, callbacks, sizeof(callbacks), C_NULL)
end

function fuse_session_mount(se, mountpoint)
    ccall((:fuse_session_mount, :libfuse3), Cint, (Ptr{Nothing}, Cstring), se, mountpoint)
end

function fuse_session_loop(se)
    ccall((:fuse_session_loop, :libfuse3), Cint, (Ptr{Nothing},), se)
end

function fuse_session_loop_mt(se, cfg::FuseLoopConfig)
    ccall((:fuse_session_loop_mt, :libfuse3), Cint, (Ptr{Nothing}, Ref{FuseLoopConfig}), se, cfg)
end

function fuse_session_unmount(se)
    ccall((:fuse_session_unmount, :libfuse3), Cvoid, (Ptr{Nothing},), se)
end

function fuse_session_destroy(se)
    ccall((:fuse_session_destroy, :libfuse3), Cvoid, (Ptr{Nothing},), se)
end

function ptr_to_userdata(u::Ptr)
    println("ptr_to_userdata($u)")
    ud = Base.unsafe_pointer_to_objref(u)
    ud isa Ref ? getindex(ud) : ud
end


text:
#__precompile__()
module StaticNetsW

using Utilities,Distributions, StatsBase,Optim, LineSearches, StatsFuns,Roots

export splitVec, distributeAinB, distributeAinB!, distributeAinVecN,maxLargeVal

## STATIC NETWORK MODEL
abstract type Ergm end
abstract type ErgmW <: Ergm end
abstract type ErgmBinW <: Ergm end
abstract type ErgmWcount <: ErgmW end
abstract type ErgmBin <: Ergm end

#constants
targetErrValStaticNets = 0.005
targetErrValStaticNetsW = 1e-5
bigConstVal = 10^6
maxLargeVal =  1e40# 1e-10 *sqrt(prevfloat(Inf))
minSmallVal = 1e2*eps()

#-------------------------Count valued Weighted Directed Networks---------------------------------------

struct  ErgmDirWCount1 <: ErgmW #Bin stands for (Binary) Adjacency matrix
    "ERGM for directed networks with geometrically(NB??) distributed count weights,
    in and out strenghts as statistics and possibility to have groups of nodes
    associated with a single pair of in out parameters. "
    obs:: Array{<:Int,1} #[InStr;OutStr] degrees of each node (even if groups are present)
    Par::Array{Array{<:Real,1},1} # One parameter per group of nodes  UNRESTRICTED VERSION  [[InPar;OutPar], scalePar ]
    groupsInds::Array{<:Real,1} #The index of the groups each node belongs to
    # two parameters per node (In and out)
 end
fooErgmDirWCount1 =  ErgmDirWCount1(ones(Int,6),[zeros(6),zeros(1)],ones(3))
ErgmDirWCount1(StrIO:: Array{<:Int,1}) =  ErgmDirWCount1(StrIO,zeros(length(StrIO[:])),[ones(round(Int,length(StrIO[:])/2)),ones(1)])

function DegSeq2graphDegSeq(Model::ErgmDirWCount1, StrIO::Array{<:Real,1}; Amat::BitArray{2} = trues(round((Int,length(StrIO)/2)),round(Int,(length(StrIO)/2)) ))
    #From a degree sequence find another one that is certanly graphicable
    #The matrix associated with the starting degree sequence is very hubs oriented

    # One stub for each unit weight
    N, StrI,StrO = splitVec(StrIO)
    sum(StrI) == sum(StrO)?(): (println(sum(StrI) - sum(StrO));error())
    inds = Vector(1:N)
    StubsI_inds = Int.(sortrows([StrI inds],rev=true))#decreasing sort
    StubsO_inds = Int.(sortrows([StrO inds],rev=true))

    mat = zeros(Int64,N,N)
    for i in 1:N
        j=1
        while j < N
            ind_i = StubsI_inds[i,2]
            ind_j = StubsO_inds[j,2]
            while (StubsI_inds[i,1]>0)& (StubsO_inds[j,1]>0)&(ind_i!=ind_j)&Amat[i,j]
                mat[ind_i,ind_j] += 1
                StubsI_inds[i,1]  -= 1
                StubsO_inds[j,1]  -= 1
                #println(i,j, StubsI_inds[i,1],StubsO_inds[j,1]  )
            end
            j += 1
        end
    end
    #display(StubsI_inds[:,1])
    ((sum(mat[Array{Bool}(eye(mat))])) != 0) ?error():()
    graphSeq= [squeeze(sum(mat,2),2); squeeze(sum(mat,1),1)]
    return graphSeq
 end

#-------------------------Countinuosly Weighted   Directed Networks Fixed A---------------------------------------

struct  ErgmDirW1Afixed <: ErgmW #Bin stands for (Binary) Adjacency matrix
    "ERGM for directed networks with Bin A fixed and Gamma distributed
     Continuous  Weights and finite probability of observing zeros.
     Thes statistics are in and out strenghts, and a binary adjacency matrix.  "
    obs:: Array{<:Real,1}  # [InStr;OutStr]
    Amat:: BitArray{2} #Binary Adjacency matrix Amat
    Par::Array{Array{<:Real,1},1} # One parameter per group of nodes  UNRESTRICTED VERSION  [ [InPar;OutPar], ScalePar]
 end
fooErgmDirW1Afixed =  ErgmDirW1Afixed(ones(6),trues(3,3),[zeros(6),zeros(1)])
ErgmDirW1Afixed(StrIO::Array{<:Real,1},Amat::BitArray{2}) =  ErgmDirW1Afixed(StrIO,Amat,[zeros(Float64,length(StrIO[:])), zeros(1)]  )

function DegSeq2graphDegSeq(Model::ErgmDirW1Afixed,StrIO::Array{<:Real,1};Amat::BitArray{2}=Model.Amat)
    #From a strengths sequence find another one that is certanly graphicable

    # pe deg2graph deg posso importare le funzioni in c++ scritte da veraart oppure fare
    # una discretizzazione (ad esempio di modo che la strength piu' piccola sia divisa in 10N
    # unità) e poi passarla in pasto alla versione Discreta pesata
   N,StrI,StrO = splitVec(StrIO)
    sum(StrI) == sum(StrO)?():error()

    minContinuousVal = minimum(StrIO[StrIO.!=0])
    # choose the minimum unit weight such that the smallest node has  10N stubs
    unitW = minContinuousVal/(10N)
    #display(minContinuousVal)
    StrCountIO = ceil(Int,StrIO./(unitW))

    graphSeq = unitW .* DegSeq2graphDegSeq(fooErgmDirWCount1,StrCountIO;Amat = Amat)
    return graphSeq
 end

bndPar2uBndPar(Model::ErgmDirW1Afixed, bndParNodesIO::Array{<:Real,1} ;indMin::Int = -1 ) =
    bndPar2uBndPar(Model, bndParNodesIO ;indMin = indMin)

uBndPar2bndPar(Model::ErgmDirW1Afixed, uBndParNodesIO::Array{<:Real,1} ;indMin::Int = -1 ) =
    uBndPar2bndPar(Model, uBndParNodesIO ;indMin = indMin)

samplSingMatCan(Model::ErgmDirW1Afixed,expMat::Array{<:Real,2},Amat::BitArray{2};α=1.0::Real) =
    (expMat[.!Amat] = 0; samplSingMatCan(Model::ErgmDirW1,expMat;α=α) )

function expMatrix2(Model::ErgmDirW1Afixed, parNodesIO::Array{<:Real,1},Amat::BitArray{2} ;logalpha = 0.0 )
    "Given the vector of model parameters (groups), return the product matrix(often useful
    in likelihood computation) and the expected matrix"
    α = exp( logalpha )
    N,parI,parO = splitVec(parNodesIO)
    parMat = exp.(parI) .+ exp.(parO)'
    expMat = putZeroDiag(α./( parMat) )
    expMat[.!Amat] = 0
    infInd = isfinite.(expMat)
    all(infInd) ?() :( display([squeeze(prod(infInd,2),2) squeeze(prod(infInd,1),1)]); error())
    return parMat,expMat
 end
expMatrix(Model::ErgmDirW1Afixed, parNodesIO::Array{<:Real,1} ,Amat::BitArray{2};logalpha = 0.0  ) =
        expMatrix2(Model, parNodesIO,Amat;logalpha = logalpha )[2]# return only the expected matrix given the parameters
expValStatsFromMat(Model::ErgmDirW1Afixed,expMat::Array{<:Real,2}  ) =
                        expValStatsFromMat(fooErgmDirBin1,expMat)

expValStats(Model::ErgmDirW1Afixed, parNodesIO::Array{<:Real,1},Amat::BitArray{2} ) = expValStatsFromMat(Model, expMatrix(Model,parNodesIO ,Amat)) # if only the paramters are given as input
function firstOrderCond(Model::ErgmDirW1Afixed;strIO::Array{<:Real,1} = Model.obs,parNodesIO::Array{<:Real,1} = Model.Par, Amat::BitArray{2}=Model.Amat )
    "Given the model, the strengths sequence, the parameters and Adjacency matrx ,
    return the First order conditions. Gradient of the loglikelihood, or system
    of differences between expected and observed values of the statistics (either
    degrees(In and Out) of single nodes or of a group). "
    foc = strIO .-expValStats(Model,parNodesIO,Amat)
    return foc
 end
function identify!(Model::ErgmDirW1Afixed,parIO::Array{<:Real,1})
    "Given a vector of parameters, return the transformed vector that verifies
    an identification condition. Do not do anything if the model is identified."
    #set the first of the in parameters equal to one (Restricted version)
    # N,parI,parO = splitVec(parIO)
    # parIO = log.([ exp.(parI) - exp.(parI[1]) ;  exp.(parO) + exp.(parI[1]) ])
    return parIO
  end

function estimateSlow(Model::ErgmDirW1Afixed; strIO::Array{<:Real,1} = Model.obs,
                Amat::BitArray{2}=falses(10,10) , targetErr  =targetErrValStaticNetsW,bigConst   = bigConstVal,identIter = false,startVals::Array{Float64,1} =zeros(10) )
    "Given model type, observations of the statistics and groups assignments,
    estimate the parameters."
    N,strI,strO = splitVec(strIO)
    (sum(strI) - sum(strO) )< 10^Base.eps()? S =sum(strO) :error()
    SperLink  = S/(N*N-N)
    # The parameters inside this function are Restricted to be in [0,1] hence
    # they need to be converted (log(par)) in order to be used outside

    # if no Amat is given assume the one implied by the strengths sequence
    sum(Amat ) == 0? Amat = (strI.>0).*(strO'.>0):()
    #uniform parameter as starting value
    unifstartval = 0.5/SperLink;
    bigConstVal = 1e4
    if sum(startVals) == 0
        local parI = unifstartval*ones(Float64,N)
        local parO = unifstartval*ones(Float64,N)
    else
        N,parI,parO = splitVec(startVals)
    end

    #bigNumber = 1e3./targetErr
    eps =  targetErr*1e-12 #10*Base.eps()
    i=0
    maxIt = 500
    relErr = 1
    tol_x = 10*Base.eps()

    function tmpFindZero(objFun::Function,brak::Array{<:Real,1},start::Real)
        #    wrapper find_zero options
        println((brak[1],brak[2]))
        println((objFun(brak[1]),objFun(brak[2])))
        #println()
        out = find_zero(objFun,(brak[1],brak[2]),Bisection())
        println((out,objFun(out)))
        return out
    end
    degsI = sumSq(Amat,2)
    degsO = sumSq(Amat,1)
    zerIndI = (strI .== 0 )
    (zerIndI == (degsI.==0) )? () : error()
    zerIndO = (strO .== 0 )
    @show (zerIndO)
    @show  (degsO.==0)
    (zerIndO == (degsO.==0) )? () : error()

    while (i<maxIt)&(relErr>targetErr)
        for n=1:N
            indGzO= (parO.>0).& (.!zerIndO); indGzO[n] = false
            indLzO= (parO.<0).& (.!zerIndO); indLzO[n] = false
            lBndParI = maximum( [(sum(indGzO)>0? -minimum(parO[indGzO]): [] ); (sum(indLzO)>0 ? -minimum(parO[indLzO]) : []) ] )

            println((1,n))
            #put to zero each component of the gradient1
            function objfunI_n(parI_n::Real)
                parVecOthGroups = parI_n .+ parO[Amat[:,n]]
                expStrI =  sum(  1./( parVecOthGroups) )   ##parO_n = parO[n]   + (NgroupsMemb[ng,2]-1)* ( 1/(parI_ng + parO_ng)) )
                return  strI[n] - expStrI

            end
            #find zero of given component

            if strI[n] <= 1e-3
                #error("Remove zero col and rows")
                parI[n] = bigConstVal
                #println(0)
            else
                tmpBnd = 1/(strI[n]/(N*(N-1)) )
                #println((lBndParI))
                brakets =[ lBndParI + 1/S ;  bigConstVal]
                #ERRORE PROBABILMENTE IN SCAMIO TRA MATRICE A E A'
                #DEVO ANCHE CORREGGERE IL LOWER BOUND PER TENERE CONTO DGLI
                #ELEMENTI NULLI DELLA AMAT
                parI[n] = tmpFindZero(objfunI_n,brakets,parI[n])
            end
        end
        parUnIO = ([parI;parO])
        identIter? parUnIO = identify(fooErgmDirW1Afixed,parUnIO):()
        #println(parI)
        err = firstOrderCond(Model; strIO = strIO,parGroupsIO = parUnIO , groupsInds = groupsInds)
        #display(err)

        for n=1:N

            indGzI=(parI.>0).& (.!zerIndI); indGzI[n] = false
            indLzI=(parI.<0).& (.!zerIndI); indLzI[n] = false
            lBndParO = maximum( [(sum(indGzI)>0? -minimum(parI[indGzI]): [] ); (sum(indLzI)>0 ? -minimum(parI[indLzI]) : []) ] )

            #println((2,nInd))
            #put to zero each component of the gradient1
            function objfunO_n(parO_n::Real)
                parVecOthGroups = parO_n .+ parI[1:end ][Amat[n,1:end]]
                expStrO = sum(1./(parVecOthGroups))#NgroupsMemb[ng,2]* sum( NgroupsMemb[1:end .!=ng ,2].*(1./( parVecOthGroups)) )    + (NgroupsMemb[ng,2]-1)* (1/(parI_ng + parO_ng)) )
                return  strO[n] - expStrO
            end
            #find zero of given component
            if strO[n]<= 1e-3
                #error("Remove zero col and rows")
                parO[n] = bigConstVal
                #println(bigConstVal)
                #println(0)
            else
                tmpBnd = 1/(strGroupsO[n]/(NgroupsMemb[n,2]*(N-1)) )
                #println((lBndParO))
                brakets =[lBndParO + 1/S ;  bigConstVal ]
                parO[n] = tmpFindZero(objfunO_n,brakets,parO[n])
            end

        end
        parUnIO = ([parI;parO])
        #println(parO)
        err = firstOrderCond(Model; strIO = strIO,parGroupsIO = parUnIO , groupsInds = groupsInds)
        #display(err)
    #    tmpMat = par.*par';expMat = putZeroDiag!(tmpMat./(1+tmpMat)) ;expDeg = squeeze(sum(expMat,2),2)
        relErrVec = abs.(err./strGroups)[strGroups .>= 1e-3]
        #println(relErrVec)
        #take into account also groups with zero strength. multiply by target error
        # implies that the expected strength for groups with zero strength
        # has to be lower than 1
        absErrRescVec = abs.(err[strGroups.==0]) .* targetErr
        relErr = maximum([relErrVec;absErrRescVec])
        isnan(relErr)?error():()
        # display(err')
         #println((mean(relErrVec[deg.!=0]),relErr))
        i+=1
        #i==1?error():()
    end
    #println([parI parO])
    i<maxIt?():(println(i);error())
    # the parameter of the first index has to be unique because I want it to be the only infinite
    strIO[1]==0 ? parI[1] = bigConstVal*50:()
    parIO = [parI;parO]
    #println(unique(parIO))
    α = exp(Model.Par[2])
    parIOout = parIO #.* α

    outPar = identify!(Model,parIO)
    outMod = ErgmDirW1(strIO,[outPar,Model.Par[2]],groupsInds)

    return outPar, i , outMod
 end
function sampl(Model::ErgmDirW1Afixed,Nsample::Int;  parNodesIO::Array{Array{<:Real,1}}=Model.Par,Amat::BitArray{2}=Model.Amat  )
    parNodesIO = [parI[groupsInds];parO[groupsInds]]
    expMat = expMatrix(Model, parNodesIO,Amat)
    N = length(expMat[:,1])
    SampleMats =zeros(Float64,N,N,Nsample)
    logα = parArray[2][1]
    α = exp(logα)
    for s=1:Nsample
        SampleMats[:,:,s] = samplSingMatCan(Model,expMat,Amat,α)
    end
    return SampleMats
 end
function estimateIPFMat(strIO::Array{<:Real,1},A::Array{<:Real,2} ; targetErr  =100 * targetErrValStaticNetsW)
    "Given In and Out strength sequences and a binary matrix A, return the matrix obtained with IPF "
    N,sI,sO = splitVec(strIO)
    nnzInds = strIO.!=0
    ~,nnzIndsI,nnzIndsO = splitVec(nnzInds)
    W = Float64.(A)
    it = 0;maxAbsRelErr = 10
    while maxAbsRelErr > targetErr
        sIprev = sum(W,2)

        W[nnzIndsI,:] = W[nnzIndsI,:].*(sI[nnzIndsI]./sIprev[nnzIndsI])
        sOprev = sum(W,1)
        W[:,nnzIndsO] = (sO[nnzIndsO]./sOprev[nnzIndsO])'.*W[:,nnzIndsO]
        strIO_IPF = expValStatsFromMat(Model,W)
        err = strIO_IPF - strIO
        relErr = err[nnzInds]./strIO[nnzInds]
        maxAbsRelErr = maximum(abs.(relErr))
        it +=1
    end
    return W,it
 end

 #-------------------------Countinuosly Weighted Directed Networks---------------------------------------

 struct  ErgmDirW1 <: ErgmW #Bin stands for (Binary) Adjacency matrix
     "ERGM for directed networks with Gamma distributed Continuous Weights,
     in and out strenghts as statistics and possibility to have groups of nodes
     associated with a single pair of in out parameters. "
     obs:: Array{<:Real,1} #[InStr;OutStr] degrees of each node (even if groups are present)
     Par::Array{Array{<:Real,1},1} # One parameter per group of nodes  UNRESTRICTED VERSION  [ [InPar;OutPar], ScalePar]
     groupsInds::Array{Int,1} #The index of the groups each node belongs to
     # two parameters per node (In and out)
  end
 fooErgmDirW1 =  ErgmDirW1(ones(6),[zeros(6),zeros(1)],ones(3))
 ErgmDirW1(StrIO:: Array{<:Real,1}) =  ErgmDirW1(StrIO,[zeros(Float64,length(StrIO[:])), zeros(1)] ,Int.(1:round(Int,length(StrIO[:])/2))  )

 function DegSeq2graphDegSeq(Model::ErgmDirW1,StrIO::Array{<:Real,1})
     #From a degree sequence find another one that is certanly graphicable

     # pe deg2graph deg posso importare le funzioni in c++ scritte da veraart oppure fare
     # una discretizzazione (ad esempio di modo che la strength piu' piccola sia divisa in 10N
     # unità) e poi passarla in pasto alla versione Discreta pesata
     N,StrI,StrO = splitVec(StrIO)
     sum(StrI) == sum(StrO)?():error()

     minContinuousVal = minimum(StrIO[StrIO.!=0])
     # choose the minimum unit weight such that the smallest node has  10N stubs
     unitW = minContinuousVal/(10N)
     #display(minContinuousVal)
     StrCountIO = ceil(Int,StrIO./(unitW))
     graphSeq = unitW .* DegSeq2graphDegSeq(fooErgmDirWCount1,StrCountIO)
     return graphSeq
  end

 function bndPar2uBndPar(Model::ErgmDirW1, bndParNodesIO::Array{<:Real,1} ;indMin::Int = -1 )
     N,bndParI, bndParO = splitVec(bndParNodesIO)
     bndParI[1] == 0 ? () :error()
     uBndParO = log.(bndParO)
     indMin == -1 ? indMin = indmin(bndParO):()
     #println(indMin)
     # Define the unbounded version of the In par s.t. it can really takeall allowed values
     # without violating λ_i +λ_j > 0
     uBndParI = log.(bndParO[indMin] .+ bndParI )
     uBndParI[1] = -Inf
     uBndParIO = [uBndParI;uBndParO]
     return  uBndParIO
  end
 function uBndPar2bndPar(Model::ErgmDirW1, uBndParNodesIO::Array{<:Real,1} ;indMin::Int = -1 )
     N,uBndParI, uBndParO = splitVec(uBndParNodesIO)
     uBndParI[1] == -Inf ? () :error()
     expLambdaO = exp.(uBndParO)
     indMin == -1 ? indMin = indmin(expLambdaO):()
     # Define the unbounded version of the In par s.t. it can really takeall allowed values
     # without violating λ_i +λ_j > 0
     bndParI = - expLambdaO[indMin] .+ exp.(uBndParI)
     #identification
     bndParI[1] = 0
     bndParO =   expLambdaO + minSmallVal

     #bndParO[indzeroO] = minSmallVal
     bndParIO = [bndParI;bndParO]
     indInf = bndParIO.>= maxLargeVal

     bndParIO[indInf] = maxLargeVal
     return  bndParIO
  end
 function linSpacedPar(Model::ErgmDirW1,Nnodes::Int;Ngroups = Nnodes,deltaN::Int=3,graphConstr = true)
     error()
     "for a given static model return a widely spaced sequence of groups parameters"
     N=Nnodes
     #generate syntetic degree distributions, equally spaced, and estimate the model
     tmpdegsIO = round.([ Vector(linspace(deltaN,N-1-deltaN,N));Vector(linspace(deltaN,N-1-deltaN,N))])
     if graphConstr
         tmpdegsIO  =   DegSeq2graphDegSeq(fooErgmDirBin1,tmpdegsIO)
     end

     if Ngroups>N#if more groups than nodes are required throw an error
         error()
     else
         groupsInds = distributeAinVecN(Vector(1:Ngroups),Nnodes)
     end
     out,it,estMod = estimate(ErgmDirW1(tmpdegsIO);groupsInds = groupsInds)
     return out,tmpdegsIO
  end
 function samplSingMatCan(Model::ErgmDirW1,expMat::Array{<:Real,2};α=1.0::Real)
     """
     given the expected matrix sample one ranUtilities matrix from the corresponding pdf
     """
     N1,N2 = size(expMat)
     N1!=N2?error():N=N1
     out = zeros(Float64,N,N)
     for c = 1:N
         for r=1:N
             if r!=c
                 #The Gamma in Distributions takes as input α and θ = expVal/α
                 #println(expMat[r,c]/α)
                 expMat[r,c]/α > 0 ? () :  (println(expMat[r,c]);  error() )
                 out[r,c] = rand(Gamma(α,expMat[r,c]/α))
             end
         end
     end
     return out
  end
 function expMatrix2(Model::ErgmDirW1, parNodesIO::Array{<:Real,1} ;logalpha = 0.0 )
     "Given the vector of model parameters (groups), return the product matrix(often useful
     in likelihood computation) and the expected matrix"
     α = exp( logalpha )
     N,parI,parO = splitVec(parNodesIO)
     parMat = parI .+ parO'
     expMat = putZeroDiag(α./( parMat) )
     infInd = isfinite.(expMat)
     all(infInd) ?() :( display([squeeze(prod(infInd,2),2) squeeze(prod(infInd,1),1)]); error())
     return parMat,expMat
  end
 expMatrix(Model::ErgmDirW1, parNodesIO::Array{<:Real,1} ;logalpha = 0.0  ) =
         expMatrix2(Model, parNodesIO;logalpha = logalpha )[2]# return only the expected matrix given the parameters
 expValStatsFromMat(Model::ErgmDirW1,expMat::Array{<:Real,2}  ) =
                         expValStatsFromMat(fooErgmDirBin1,expMat)

 expValStats(Model::ErgmDirW1, parNodesIO::Array{<:Real,1} ) = expValStatsFromMat(Model, expMatrix(Model,parNodesIO )) # if only the paramters are given as input
 function firstOrderCond(Model::ErgmDirW1;strIO::Array{<:Real,1} = Model.obs,
                     strGroups::Array{<:Real,1} = zeros(4), parGroupsIO::Array{<:Real,1}=Model.Par, groupsInds::Array{<:Real,1} = Model.groupsInds )
     "Given the model, the degree sequence, the parameters and groups assignements,
     return the First order conditions. Gradient of the loglikelihood, or system
     of differences between expected and observed values of the statistics (either
     degrees(In and Out) of single nodes or of a group). "
     N,strI,strO = splitVec(strIO)
     groupsNames = unique(groupsInds)
     NG = length(groupsNames)
     NGcheck,parI,parO = splitVec(parGroupsIO)
     NG==NGcheck? ():error()
     if sum(strGroups)==0
         strGroups = zeros(2NG) ; for i =1:N strGroups[groupsInds[i]] += strI[i];strGroups[NG + groupsInds[i]] += strO[i]; end
     end
     parNodesIO = [parI[groupsInds];parO[groupsInds]]
     N,expStrI,expStrO = splitVec(expValStats(Model,parNodesIO))
     expGroupsStr = zeros(2NG) ; for i =1:N expGroupsStr[groupsInds[i]] += expStrI[i];expGroupsStr[NG + groupsInds[i]] += expStrO[i]; end
     # println(groupsInds)
     # println(size(strGroups))
     # println(size(expGroupsStr))
     foc = strGroups .- expGroupsStr
     return foc
  end
 function identify!(Model::ErgmDirW1,parGroupsIO::Array{<:Real,1})
     "Given a vector of parameters, return the transformed vector that verifies
     an identification condition. Do not do anything if the model is identified."
     #set the first of the in parameters equal to one (Restricted version)
     N,parI,parO = splitVec(parGroupsIO)
     parIO = [ parI - parI[1] ;  parO + parI[1] ]
     return parIO
   end

 function estimate(Model::ErgmDirW1; strIO::Array{<:Real,1} = Model.obs,
                     groupsInds = Model.groupsInds, targetErr  =targetErrValStaticNetsW,bigConst   = bigConstVal)
     "Given model type, observations of the statistics and groups assignments,
     estimate the parameters."
     N,strI,strO = splitVec(strIO)
     groupsNames = unique(groupsInds)
     NG = length(groupsNames)
     tmpDic = countmap(groupsInds)
     #number of nodes per group
     NgroupsMemb = sortrows([[i for i in keys(tmpDic)] [  tmpDic[i] for i in keys(tmpDic) ] ])
     # agrregate degree of each group
     strGroups = zeros(Float64,2NG) ; for i =1:N strGroups[groupsInds[i]] += strI[i]; strGroups[NG + groupsInds[i]] += strO[i]; end
     ~,strGroupsI , strGroupsO = splitVec(strGroups)
     zerIndI = strGroupsI.==0
     zerIndO = strGroupsO.==0
     #sort the groups for degrees
     indsStrGroupsOrderI = Int.(sortrows([strGroupsI Vector(1:NG)], rev = true)[:,2])
     indsStrGroupsOrderO = Int.(sortrows([strGroupsO Vector(1:NG)], rev = true)[:,2])

     (sum(strGroupsI) - sum(strGroupsO) )< 10^Base.eps()? S =sum(strGroupsO) :error()
     SperLink  = S/(N*N-N)
     # The parameters inside this function are Restricted to be in [0,1] hence
     # they need to be converted (log(par)) in order to be used outside

     #uniform parameter as starting value
     unifstartval = 0.5/SperLink;
     bigConstVal = 1e4
     local parI = unifstartval*ones(Float64,NG)
     local parO = unifstartval*ones(Float64,NG)
     if NG==1
         error()
         parI[1] = 0.0
         parO[1] = log(unifstartval)
         outParUN = [ parI; parO ]
         outMod = ErgmDirW1(Model.obs,outParUN,groupsInds)
         i = 0
     else
         #bigNumber = 1e3./targetErr
         eps =  targetErr*1e-12 #10*Base.eps()

         i=0
         maxIt = 500
         relErr = 1
         tol_x = 10*Base.eps()

         function tmpFindZero(objFun::Function,brak::Array{<:Real,1},start::Real)
             #    wrapper find_zero options
             #println((brak[1],brak[2]))
             #println((objFun(brak[1]),objFun(brak[2])))
             #println()
             out = find_zero(objFun,(brak[1],brak[2]),Bisection())
             #println((out,objFun(out)))
             return out
         end
         while (i<maxIt)&(relErr>targetErr)

             for nInd=1:NG
                 ng = indsStrGroupsOrderI[nInd]

                 indGzO= (parO.>0).& (.!zerIndO); indGzO[ng] = false
                 indLzO= (parO.<0).& (.!zerIndO); indLzO[ng] = false
                 lBndParI = maximum( [(sum(indGzO)>0? -minimum(parO[indGzO]): [] ); (sum(indLzO)>0 ? -minimum(parO[indLzO]) : []) ] )

                 #println((1,nInd))
                 #put to zero each component of the gradient1
                 function objfunI_ng(parI_ng::Real)
                     parVecOthGroups = parI_ng .+ parO[1:end .!=ng]
                     parO_ng = parO[ng]
                     expGroupStrI = NgroupsMemb[ng,2]*( sum( NgroupsMemb[1:end .!=ng ,2].*( 1./( parVecOthGroups)) )    + (NgroupsMemb[ng,2]-1)* ( 1/(parI_ng + parO_ng)) )
                     return  strGroupsI[ng] - expGroupStrI

                 end
                 #find zero of given component

                 if strGroupsI[ng] <= 1e-3
                     #error("Remove zero col and rows")
                     parI[ng] = bigConstVal
                     #println(0)
                 else
                     tmpBnd = 1/(strGroupsI[ng]/(NgroupsMemb[ng,2]*(N-1)) )
                     #println((lBndParI))
                     brakets =[ lBndParI + 1/S ;  bigConstVal]
                     parI[ng] = tmpFindZero(objfunI_ng,brakets,parI[ng])
                 end
             end
             parUnIO = ([parI;parO])
             #println(parI)
             err = firstOrderCond(Model; strIO = strIO,parGroupsIO = parUnIO , groupsInds = groupsInds)
             #display(err)

             for nInd=1:NG
                 ng = indsStrGroupsOrderO[nInd]

                 indGzI=(parI.>0).& (.!zerIndI); indGzI[ng] = false
                 indLzI=(parI.<0).& (.!zerIndI); indLzI[ng] = false
                 lBndParO = maximum( [(sum(indGzI)>0? -minimum(parI[indGzI]): [] ); (sum(indLzI)>0 ? -minimum(parI[indLzI]) : []) ] )

                 #println((2,nInd))
                 #put to zero each component of the gradient1
                 function objfunO_ng(parO_ng::Real)
                     parVecOthGroups = parO_ng .+ parI[1:end .!=ng]
                     parI_ng = parI[ng]
                     expGroupStrO = NgroupsMemb[ng,2]*( sum( NgroupsMemb[1:end .!=ng ,2].*(1./( parVecOthGroups)) )    + (NgroupsMemb[ng,2]-1)* (1/(parI_ng + parO_ng)) )
                     return  strGroupsO[ng] - expGroupStrO
                 end
                 #find zero of given component
                 if strGroupsO[ng]<= 1e-3
                     #error("Remove zero col and rows")
                     parO[ng] = bigConstVal
                     #println(bigConstVal)
                     #println(0)
                 else
                     tmpBnd = 1/(strGroupsO[ng]/(NgroupsMemb[ng,2]*(N-1)) )
                     #println((lBndParO))
                     brakets =[lBndParO + 1/S ;  bigConstVal ]
                     parO[ng] = tmpFindZero(objfunO_ng,brakets,parO[ng])
                 end

             end
             parUnIO = ([parI;parO])
             #println(parO)
             err = firstOrderCond(Model; strIO = strIO,parGroupsIO = parUnIO , groupsInds = groupsInds)
             #display(err)
         #    tmpMat = par.*par';expMat = putZeroDiag!(tmpMat./(1+tmpMat)) ;expDeg = squeeze(sum(expMat,2),2)
             relErrVec = abs.(err./strGroups)[strGroups .>= 1e-3]
             #println(relErrVec)
             #take into account also groups with zero strength. multiply by target error
             # implies that the expected strength for groups with zero strength
             # has to be lower than 1
             absErrRescVec = abs.(err[strGroups.==0]) .* targetErr
             relErr = maximum([relErrVec;absErrRescVec])
             isnan(relErr)?error():()
             # display(err')
              #println((mean(relErrVec[deg.!=0]),relErr))
             i+=1
             #i==1?error():()
         end
         #println([parI parO])
         i<maxIt?():(println(i);error())
         # the parameter of the first index has to be unique because I want it to be the only infinite
         strIO[1]==0 ? parI[1] = bigConstVal*50:()
         parIO = [parI;parO]
         #println(unique(parIO))
         α = exp(Model.Par[2])
         parIOout = parIO #.* α

         outPar = identify!(Model,parIO)
         outMod = ErgmDirW1(strIO,[outPar,Model.Par[2]],groupsInds)
     end
     return outPar, i , outMod
  end
 function sampl(Model::ErgmDirW1,Nsample::Int;  parArray::Array{Array{<:Real,1}}=Model.Par, groupsInds = Model.groupsInds )
     parGroupsIO = parArray[1]
     NGcheck,parI,parO = splitVec(parGroupsIO)
     parNodesIO = [parI[groupsInds];parO[groupsInds]]
     expMat = expMatrix(Model, parNodesIO)
     N = length(expMat[:,1])
     SampleMats =zeros(Float64,N,N,Nsample)
     logα = parArray[2][1]
     α = exp(logα)
     for s=1:Nsample
         SampleMats[:,:,s] = samplSingMatCan(Model,expMat,α)
     end
     return SampleMats
  end

#--------------------------- Snapshot sequence functions

function zeroStrParFun(Model::ErgmDirW1Afixed,N::Int;method="SMALL" ,smallProb = 1e-1)
    #define a number big enough to play the role of Inf for
    # purposes of sampling N(N-1) bernoulli rvs with prob 1/(1+exp(  bigNumb))
    # prob of sampling degree > 0 (N-1)/(1+exp(bigNumb)) < 1e-6
    if uppercase(method) == "INF"
        # zero degrees have -Inf fitnesses
        zeroDegPar = - Inf #log((N-1)/smallProb - 1)
    elseif uppercase(method) == "10AS21"
        # such that the difference  fit(2) - fit(1) == fit(1)-fit(0) , where
        # fit(n) means the fitness associated with a node of deg n.
        # when more than one choice is possible (the fitness depends on both degrees)
        # than choose fit(2) and fit(1) that maximize the difference

    elseif uppercase(method) == "SMALL"
        # such that the probability of observing a link is very small
        zeroDegPar = 100
    end

    return zeroDegPar
 end

function estimate(model::ErgmDirW1Afixed; strIO::Array{<:Real,1} = Model.obs,
                Amat::BitArray{2}=falses(10,10) , targetErr  =0.001,identIter = false)
    ## Write and test a new way to estimate DirW1 extending the iteration of
    # Chatterjee, Diaconis and  Sly to weighted networks
    N,strI,strO = splitVec(strIO)
    sum(Amat) == 0 ?  Amat = (strI.>0).*(strO'.>0) :()
    lstrI = log.(strI);lstrO = log.(strO)
    (sum(strI) - sum(strO) )< 10^Base.eps()? S =sum(strO) :error()
    SperLink  = S/(N*N-N)
    unifstartval =  -0.5*log(SperLink)
    nnzInds = strIO.!=0
    #inds = (ldegI.==0)
    #function that defines the iteration
    function phiFunc(vParUn::Array{<:Real,1})
        ~,vParUnI,vParUnO = splitVec((vParUn))
        matI = exp.( -vParUnI .+ vParUnO' )
        matO =  1./matI
        mat2sumI =  putZeroDiag(1./(1+matI ))

        mat2sumI[.!Amat ] =0
        mat2sumO =  putZeroDiag(1./(1+matO ))
        mat2sumO[.!Amat ] =0
    #    println(putZeroDiag(1./(1+matI )))
        outPhiI = ( - lstrI + log.( sumSq( mat2sumI ,2) ) )
        outPhiO = ( - lstrO + log.( sumSq( mat2sumO,1) ) )
        #println(vParReI[ldegI.==0])
        #println(log.(sumSq(matI,2))[inds])
        #println(matI[inds,:][1,1:20])
        #println(matI[inds,:][2,1:20])
        return [outPhiI;outPhiO]
     end
    function checkFit(vParUn::Vector{<:Real})
        #println(vParUn)
        expMat = expMatrix(model,vParUn,Amat)

        expStrIO =  [sumSq(expMat,2) ; sumSq(expMat,1) ]
        errsIO = (expStrIO - strIO)[nnzInds]
        relErrsIO = abs.(errsIO)./(strIO[nnzInds])
        # mInd = indmax(relErrsIO)
        # println((mInd,relErrsIO[mInd],strIO[nnzInds][mInd],expStrIO[nnzInds][mInd],vParUn[nnzInds][mInd] ))
        # println(sum(expMat[mInd,:]) )
        # tmp = Vector(1:2N)
        # println(tmp[nnzInds][mInd])

        return prod(relErrsIO .< targetErr)
     end

    maxIt = 200
    it=0
    it2 = 0
    vParUn = [-lstrI;-lstrO]# 0*unifstartval.*ones(2N)
    #starting from an identified set of par. and then identify at each step
    #amounts to doing an identified step. i.e. moving summing an identified
    #vector
    identIter?vParUn =  identify!(model,vParUn):()
    while (!checkFit(vParUn)) & (it<maxIt)
        println(it)
        vParUn = phiFunc(vParUn); it+=1
        identIter?vParUn =  identify!(model,vParUn):()
    end
    #if did not obtain the required precision try with the slower method
    vParUn[.!nnzInds] = zeroStrParFun(model,N)
    if it==maxIt

        startSlowVals = exp.(vParUn)
        outSlow = estimateSlow(model; strIO = strIO,targetErr = targetErr,
                            startVals =startSlowVals)

        vParUn = identify!(model,outSlow[1])
        it2 = outSlow[2]
    end
    #zeroInds = degIO.==0
    #vParUn[zeroInds] = bigNegParFun(model,N)
    return vParUn,it,it2
 end
 #

struct SnapSeqNetDirW1
    obsT:: Array{<:Real,2}  # In and Out strengths for each t . The last
                            # dimension is always time. 2NxT
    Par::Array{<:Real,2} # In and Out Parameters [InW_T;outW_T] for each t
 end
SnapSeqNetDirW1(obs::Array{<:Real,2}) = SnapSeqNetDirW1(obs,zeros(obs))
fooSnapSeqNetDirW1 = SnapSeqNetDirW1(zeros(2,3))
function estimate(Model::SnapSeqNetDirW1; strIO_T::Array{<:Real,2}=Model.obsT,
                        targetErr::Real=1e-3,identPost=false,identIter=false,zeroDegFit="small")
   N2,T = size(strIO_T);N = round(Int,N2/2)
   UnParT = zeros(2N,T)
   prog = 0
   for t = 1:T
       str_t = strIO_T[:,t]
       UnPar_t , Niter1,Niter2  = estimate(fooErgmDirW1Afixed; strIO = str_t ,Amat = falses(10,10), targetErr =  targetErr,identIter=identIter)
       identPost?UnParT[:,t] = identify!(fooErgmDirW1Afixed,UnPar_t): UnParT[:,t] = UnPar_t
        round(t/T,1)>prog? (prog = round(t/T,1);println(prog)):()
        #println((t,Niter1,Niter2))
   end
   UnParT[strIO_T.==0] = zeroDegParFun(fooErgmDirW1Afixed,N;method=zeroDegFit )
 return UnParT
 end

 function sampl(Model::SnapSeqNetDirW1, vParUnIO_T::Array{<:Real,2}; Nsample::Int=1  )
    N2,T = size(vParUnIO_T);N = round(Int,N2/2)
    SampleMats =zeros(Float64,Nsample,N,N,T)
    for s=1:Nsample
        for t = 1:T
            SampleMats[s,:,:,t] = samplSingMatCan(fooErgmDirW1Afixed,expMatrix(fooErgmDirW1Afixed,vParUnIO_T[:,t]))
        end
    end
    Nsample == 1 ? (return squeeze(SampleMats,1)):(return SampleMats )
  end

 #---------------------------------Binary DIRECTED Networks with predetermined COVARIATES

end

##

ext: jl


text:
using NCDatasets
using Formatting

include("./SSM_z_res.jl")

fn_o = ARGS[1]

Nz = length(zs_SSM) - 1
N_zs = length(zs_SSM)

Dataset(fn_o, "c") do ds

    defDim(ds, "Nz", Nz)
    defDim(ds, "zs", N_zs)

    defVar(ds, "zs", zs_SSM, ("zs",))

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
@inbounds @fastmath function boundary_force_snap_p!(issp::Int64,pac::Paramc,pass::Vector{Paramss},pap::Paramp)
	p=pap.p
	boundary=pass[issp].boundary[5]
	ps=view(p[1],:,:,1)
	bs=view(boundary,:,:,1)
	copyto!(ps,bs)
end
@inbounds @fastmath function boundary_force_snap_vxvz!(issp::Int64,pac::Paramc,pass::Vector{Paramss},pap::Paramp)
	# initial conditions from boundary for first propagating field only
	p=pap.p
	pw=p[1]
	boundary=pass[issp].boundary[5]
	ps=view(pw,:,:,2)
	bs=view(boundary,:,:,2)
	copyto!(ps,bs)
	ps=view(pw,:,:,3)
	bs=view(boundary,:,:,3)
	copyto!(ps,bs)
end
@fastmath @inbounds function boundary_force!(it::Int64,issp::Int64,pac::Paramc,pass::Vector{Paramss},pap::Paramp)
	boundary=pass[issp].boundary
	p=pap.p
	ibx0=pac.ibx0; ibz0=pac.ibz0; ibx1=pac.ibx1; ibz1=pac.ibz1
	boundaryf_l!(it,ibx0,ibx1,ibz0,ibz1,p,boundary)
	boundaryf_r!(it,ibx0,ibx1,ibz0,ibz1,p,boundary)
	boundaryf_t!(it,ibx0,ibx1,ibz0,ibz1,p,boundary)
	boundaryf_b!(it,ibx0,ibx1,ibz0,ibz1,p,boundary)
end
@fastmath @inbounds function boundaryf_l!(it,ibx0,ibx1,ibz0,ibz1,p,boundary)
	pw=p[1]
	for ix=1:3
		@simd for iz=1:ibz1-ibz0+1
			pw[ibz0+iz-1,ibx0+ix-1,1] = boundary[4][iz,ix,it]
		end
	end
end
@fastmath @inbounds function boundaryf_r!(it,ibx0,ibx1,ibz0,ibz1,p,boundary)
	pw=p[1]
	for ix=1:3
		@simd for iz=1:ibz1-ibz0+1
			pw[ibz0+iz-1,ibx1-ix+1,1] = boundary[2][iz,ix,it]
		end
	end
end
@fastmath @inbounds function boundaryf_t!(it,ibx0,ibx1,ibz0,ibz1,p,boundary)
	pw=p[1]
	@simd for ix=1:ibx1-ibx0+1
		for iz=1:3
			pw[ibz0+iz-1,ibx0+ix-1,1] = boundary[1][iz,ix,it]
		end
	end
end
@fastmath @inbounds function boundaryf_b!(it,ibx0,ibx1,ibz0,ibz1,p,boundary)
	pw=p[1]
	@simd for ix=1:ibx1-ibx0+1
		for iz=1:3
			pw[ibz1-iz+1,ibx0+ix-1,1] = boundary[3][iz,ix,it]
		end
	end
end

@inbounds @fastmath function boundary_save_snap_p!(issp::Int64,pac::Paramc,pass::Vector{Paramss},pap::Paramp)
	p=pap.p
	pw=p[1]
	boundary=pass[issp].boundary[5]
	ps=view(pw,:,:,1)
	bs=view(boundary,:,:,1)
	copyto!(bs, ps)
end
@inbounds @fastmath function boundary_save_snap_vxvz!(issp::Int64,pac::Paramc,pass::Vector{Paramss},pap::Paramp)
	p=pap.p
	pw=p[1]
	boundary=pass[issp].boundary[5]
	#vx
	bs=view(boundary,:,:,2)
	ps=view(pw,:,:,2)
	copyto!(bs, ps)
	rmul!(bs,-1.)
	# vz
	bs=view(boundary,:,:,3)
	ps=view(pw,:,:,3)
	copyto!(bs, ps)
	rmul!(bs,-1.)
end
@fastmath @inbounds function boundary_save!(it::Int64,issp::Int64,pac::Paramc,pass::Vector{Paramss},pap::Paramp)
	boundary=pass[issp].boundary
	p=pap.p
	ibx0=pac.ibx0; ibz0=pac.ibz0; ibx1=pac.ibx1; ibz1=pac.ibz1
	boundarys_l!(it,ibx0,ibx1,ibz0,ibz1,p,boundary)
	boundarys_r!(it,ibx0,ibx1,ibz0,ibz1,p,boundary)
	boundarys_t!(it,ibx0,ibx1,ibz0,ibz1,p,boundary)
	boundarys_b!(it,ibx0,ibx1,ibz0,ibz1,p,boundary)
end
@fastmath @inbounds function boundarys_l!(it,ibx0,ibx1,ibz0,ibz1,p,boundary)
	pw=p[1]
	for ix=1:3
		@simd for iz=1:ibz1-ibz0+1
			boundary[4][iz,ix,it] = pw[ibz0+iz-1,ibx0+ix-1,1]
		end
	end
end
@fastmath @inbounds function boundarys_r!(it,ibx0,ibx1,ibz0,ibz1,p,boundary)
	pw=p[1]
	for ix=1:3
		@simd for iz=1:ibz1-ibz0+1
			boundary[2][iz,ix,it] = pw[ibz0+iz-1,ibx1-ix+1,1]
		end
	end
end
@fastmath @inbounds function boundarys_t!(it,ibx0,ibx1,ibz0,ibz1,p,boundary)
	pw=p[1]
	@simd for ix=1:ibx1-ibx0+1
		for iz=1:3
			boundary[1][iz,ix,it] = pw[ibz0+iz-1,ibx0+ix-1,1]
		end
	end
end
@fastmath @inbounds function boundarys_b!(it,ibx0,ibx1,ibz0,ibz1,p,boundary)
	pw=p[1]
	@simd for ix=1:ibx1-ibx0+1
		for iz=1:3
			boundary[3][iz,ix,it] = pw[ibz1-iz+1,ibx0+ix-1,1]
		end
	end
end


text:
using AutomotiveDrivingModels, AutomotivePOMDPs, POMDPs, ProgressMeter

params = UrbanParams(nlanes_main=1,
                     crosswalk_pos =  [VecSE2(6, 0., pi/2), VecSE2(-6, 0., pi/2), VecSE2(0., -5., 0.)],
                     crosswalk_length =  [14.0, 14., 14.0],
                     crosswalk_width = [4.0, 4.0, 3.1],
                     stop_line = 22.0)
env = UrbanEnv(params=params);

mdp = PedCarMDP(env = env, pos_res=3., vel_res=2., ped_birth=0.7, ped_type=VehicleDef(AgentClass.PEDESTRIAN, 1.0, 3.0))

egos = get_ego_states(mdp.env, mdp.pos_res, mdp.vel_res)

for (i, e) in enumerate(egos)
    @assert ind2ego(mdp.env, i, mdp.pos_res, mdp.vel_res) == e
end

routes = AutomotivePOMDPs.get_car_routes(mdp.env)

for route in routes
    cars = get_car_states(mdp.env, route, mdp.pos_res, mdp.vel_res)
    for (i, c) in enumerate(cars)
        @assert ind2car(mdp.env, i, route, mdp.pos_res, mdp.vel_res) == c
    end
end

peds = get_ped_states(mdp.env, mdp.pos_res, mdp.vel_ped_res)

for (i, p) in enumerate(peds)
    if ind2ped(mdp.env, i, mdp.pos_res, mdp.vel_ped_res) != p
        println(i)
        break
    end
end

state_space = states(mdp)

@showprogress for (i, s) in enumerate(state_space)
    @assert ind2state(mdp, i) == s
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
module InverseOptimalControl

import ..DynamicsModelInterface
import ..JuMPUtils
import ..CostUtils
import Ipopt
import JuMP

using ..ForwardOptimalControl: linear_dynamics_constraints, forward_quadratic_objective
using JuMP: @variable, @constraint, @objective
using UnPack: @unpack

export solve_inverse_optimal_control

#============================================= LQ case =============================================#

"The lagrangian of the forward LQR problem."
function lqr_lagrangian(x, u, λ; Q, R, A, B)
    c_dyn = linear_dynamics_constraints(x, u; A, B)
    @views forward_quadratic_objective(x, u; Q, R) +
           sum(λ[:, t]' * c_dyn[:, t] for t in axes(x)[2][1:(end - 1)])
end

"The hand-written gradient of the forward LQR problem in x."
function lqr_lagrangian_grad_x(x, u, λ; Q, R, A, B)
    T = size(x, 2)
    hcat(
        2 * Q * x[:, 1] - A[1]' * λ[:, 1], # special handling of first time step
        reduce(hcat, (2 * Q * x[:, t] + λ[:, t - 1] - A[t]' * λ[:, t] for t in 2:(T - 1))),
        2 * Q * x[:, T] + λ[:, T - 1],
    )
end

"The hand-written gradient of the forward LQR problem in u."
function lqr_lagrangian_grad_u(x, u, λ; Q, R, A, B)
    T = size(x, 2)
    hcat(reduce(hcat, (2 * R * u[:, t] - B[t]' * λ[:, t] for t in 1:(T - 1))), 2 * R * u[:, T])
end

"""
Solves aninverse LQR problem using JuMP.

`x̂` is the (observed) state trajectory for which we wish to find the cost parametrization.

`Q̃` and `R̃` are iterables of quadatic state and control cost matrices for which the weight
vectors`q` and `r` are to be estimated.
"""
function solve_inverse_lqr(
    x̂,
    Q̃,
    R̃;
    A,
    B,
    r_sqr_min = 1e-5,
    solver = Ipopt.Optimizer,
    solver_attributes = (; print_level = 3),
)
    T = size(x̂)[2]
    n_states, n_controls = size(only(unique(B)))
    opt_model = JuMP.Model(solver)
    JuMPUtils.set_solver_attributes!(opt_model; solver_attributes...)

    # decision variable
    q = @variable(opt_model, [1:length(Q̃)], lower_bound = 0)
    r = @variable(opt_model, [1:length(R̃)], lower_bound = 0)
    x = @variable(opt_model, [1:n_states, 1:T])
    u = @variable(opt_model, [1:n_controls, 1:T])
    λ = @variable(opt_model, [1:n_states, 1:(T - 1)]) # multipliers of the forward LQR condition

    # initial condition
    @constraint(opt_model, x[:, 1] .== x̂[:, 1])
    # dynamics
    @constraint(opt_model, linear_dynamics_constraints(x, u; A, B) .== 0)
    # Optimality conditions (KKT) of forward LQR show up as a constraints
    Q = sum(q .* Q̃)
    R = sum(r .* R̃)
    ∇ₓL = @constraint(opt_model, lqr_lagrangian_grad_x(x, u, λ; Q, R, A, B)[:, 2:end] .== 0)
    ∇ᵤL = @constraint(opt_model, lqr_lagrangian_grad_u(x, u, λ; Q, R, A, B) .== 0)
    # regularization
    @constraint(opt_model, r' * r >= r_sqr_min)
    @constraint(opt_model, r' * r + q' * q == 1)

    @objective(opt_model, Min, sum((x .- x̂) .^ 2))
    @time JuMP.optimize!(opt_model)

    JuMPUtils.isconverged(opt_model),
    JuMPUtils.get_values(; q, r, x, u, λ, ∇ₓL, ∇ᵤL),
    opt_model,
    JuMP.value.(Q),
    JuMP.value.(R)
end

#=========================================== Non-LQ-Case ===========================================#

function solve_inverse_optimal_control(
    y;
    control_system,
    cost_model,
    observation_model,
    fixed_inputs = (),
    init = (),
    solver = Ipopt.Optimizer,
    solver_attributes = (; print_level = 3),
    cmin = 1e-5,
    max_observation_error = nothing,
    verbose = false,
    init_with_observation = true,
)
    T = size(y)[2]
    @unpack n_states, n_controls = control_system

    opt_model = JuMP.Model(solver)
    JuMPUtils.set_solver_attributes!(opt_model; solver_attributes...)

    # decision variable
    weights = @variable(opt_model, [keys(cost_model.weights)],)
    x = @variable(opt_model, [1:n_states, 1:T])
    u = @variable(opt_model, [1:n_controls, 1:T])
    λ = @variable(opt_model, [1:n_states, 1:(T - 1)]) # multipliers of the forward optimality condition

    # initialization
    if hasproperty(init, :weights) && !isnothing(init.weights)
        verbose && @info "Using weight guess: $(init.weights)"
        for k in keys(init.weights)
            JuMP.set_start_value(weights[k], init.weights[k])
        end
    else
        verbose && @info "Using weight default"
        JuMP.set_start_value.(weights, 1 / length(weights))
    end

    # Initialization
    if init_with_observation
        # Note: This initialization is only works if `observation_model.expected_observation`
        # effectively creates an array view into x (extracting sub-components of the state).
        JuMP.set_start_value.(observation_model.expected_observation(x), y)
    end
    JuMPUtils.init_if_hasproperty!(u, init, :u)

    # constraints
    if iszero(observation_model.σ)
        @constraint(opt_model, observation_model.expected_observation(x[:, 1]) .== y[:, 1])
    end
    DynamicsModelInterface.add_dynamics_constraints!(control_system, opt_model, x, u)
    # figure out which inputs we control and fix all others
    controlled_inputs = filter(i -> i ∉ fixed_inputs, 1:n_controls)
    for i in fixed_inputs
        JuMP.fix.(u[i, :], init.u[i, :])
    end
    # Require forward-optimality for all *controlled* inputs.
    df = DynamicsModelInterface.add_dynamics_jacobians!(control_system, opt_model, x, u)
    dJ = cost_model.add_objective_gradients!(opt_model, x, u; weights)
    @constraint(
        opt_model,
        [t = 2:(T - 1)],
        dJ.dx[:, t] + λ[:, t - 1] - (λ[:, t]' * df.dx[:, :, t])' .== 0
    )
    @constraint(opt_model, dJ.dx[:, T] + λ[:, T - 1] .== 0)
    @constraint(
        opt_model,
        [t = 1:(T - 1)],
        dJ.du[controlled_inputs, t] - (λ[:, t]' * df.du[:, controlled_inputs, t])' .== 0
    )
    @constraint(opt_model, dJ.du[controlled_inputs, T] .== 0)
    # regularization
    @constraint(opt_model, weights .>= cmin)
    @constraint(opt_model, sum(weights) == 1)

    y_expected = observation_model.expected_observation(x)
    # Sometimes useful for debugging: Only search in a local neighborhood of the demonstration if we
    # have an error-bound on the noise.
    if !isnothing(max_observation_error)
        @constraint(opt_model, (y_expected - y) .^ 2 .<= max_observation_error^2)
    end

    # The inverse objective: match the observed demonstration
    @objective(opt_model, Min, sum(el -> el^2, y_expected .- y))

    time = @elapsed JuMP.optimize!(opt_model)
    verbose && @info time

    solution = merge(
        CostUtils.namedtuple(JuMPUtils.get_values(; weights)),
        JuMPUtils.get_values(; x, u, λ),
    )

    JuMPUtils.isconverged(opt_model), solution, opt_model
end

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -65.76627584345755 <= q <= 59.26265199263469)

                     add_NL_constraint(m, :(gelu(-0.4079491356153899 + -0.18134844314194787*gelu(0.2817697329190594 + 0.989386077401103*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + 0.6068435500303653*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + 0.8530798012021474*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + 0.8376755183348754*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + 0.6594644406874322*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + -0.2976334456396539*gelu(-0.4984715860832294 + -0.6460783499544172*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + 0.0952691245496573*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + -0.5434919689873383*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + -0.2400778955734424*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + 0.904210133814042*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + -0.11940205422195982*gelu(-0.9402956128299329 + 0.19522910845028774*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + -0.5236489579886725*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + -0.09802794343243137*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + -0.5359271533353085*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + 0.5362141752531211*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + -0.6350558501595733*gelu(-0.842982128122256 + 0.6833459994096582*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + 0.15560340828874075*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + 0.5368387122163134*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + -0.6086259049043772*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + -0.5403730914454505*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + 0.4683577877900529*gelu(0.6067063734003204 + 0.8661978259504042*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + -0.06088102390312189*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + 0.9836506901320718*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + 0.7952476421961925*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + -0.6703006410599177*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3])))) + gelu(0.9379679817528435 + -0.4660824895766815*gelu(0.2817697329190594 + 0.989386077401103*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + 0.6068435500303653*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + 0.8530798012021474*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + 0.8376755183348754*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + 0.6594644406874322*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + 0.32761179552569075*gelu(-0.4984715860832294 + -0.6460783499544172*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + 0.0952691245496573*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + -0.5434919689873383*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + -0.2400778955734424*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + 0.904210133814042*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + 0.6880683972685175*gelu(-0.9402956128299329 + 0.19522910845028774*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + -0.5236489579886725*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + -0.09802794343243137*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + -0.5359271533353085*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + 0.5362141752531211*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + 0.5888723964207281*gelu(-0.842982128122256 + 0.6833459994096582*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + 0.15560340828874075*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + 0.5368387122163134*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + -0.6086259049043772*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + -0.5403730914454505*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + -0.05595383404471299*gelu(0.6067063734003204 + 0.8661978259504042*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + -0.06088102390312189*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + 0.9836506901320718*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + 0.7952476421961925*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + -0.6703006410599177*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3])))) + gelu(-0.3048668135843893 + -0.3680078347752045*gelu(0.2817697329190594 + 0.989386077401103*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + 0.6068435500303653*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + 0.8530798012021474*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + 0.8376755183348754*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + 0.6594644406874322*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + -0.09989706667818732*gelu(-0.4984715860832294 + -0.6460783499544172*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + 0.0952691245496573*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + -0.5434919689873383*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + -0.2400778955734424*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + 0.904210133814042*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + -0.6733365176188579*gelu(-0.9402956128299329 + 0.19522910845028774*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + -0.5236489579886725*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + -0.09802794343243137*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + -0.5359271533353085*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + 0.5362141752531211*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + -0.7967060490111804*gelu(-0.842982128122256 + 0.6833459994096582*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + 0.15560340828874075*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + 0.5368387122163134*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + -0.6086259049043772*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + -0.5403730914454505*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + 0.8483184805894783*gelu(0.6067063734003204 + 0.8661978259504042*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + -0.06088102390312189*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + 0.9836506901320718*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + 0.7952476421961925*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + -0.6703006410599177*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3])))) + gelu(-0.8260325513973243 + 0.9282935832138137*gelu(0.2817697329190594 + 0.989386077401103*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + 0.6068435500303653*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + 0.8530798012021474*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + 0.8376755183348754*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + 0.6594644406874322*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + -0.5968467249691596*gelu(-0.4984715860832294 + -0.6460783499544172*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + 0.0952691245496573*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + -0.5434919689873383*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + -0.2400778955734424*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + 0.904210133814042*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + -0.003449906971975558*gelu(-0.9402956128299329 + 0.19522910845028774*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + -0.5236489579886725*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + -0.09802794343243137*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + -0.5359271533353085*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + 0.5362141752531211*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + 0.5427745947324003*gelu(-0.842982128122256 + 0.6833459994096582*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + 0.15560340828874075*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + 0.5368387122163134*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + -0.6086259049043772*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + -0.5403730914454505*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + 0.4794292841893828*gelu(0.6067063734003204 + 0.8661978259504042*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + -0.06088102390312189*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + 0.9836506901320718*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + 0.7952476421961925*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + -0.6703006410599177*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3])))) + gelu(-0.3442450542855955 + 0.3679171232252214*gelu(0.2817697329190594 + 0.989386077401103*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + 0.6068435500303653*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + 0.8530798012021474*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + 0.8376755183348754*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + 0.6594644406874322*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + 0.4969759600474797*gelu(-0.4984715860832294 + -0.6460783499544172*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + 0.0952691245496573*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + -0.5434919689873383*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + -0.2400778955734424*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + 0.904210133814042*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + -0.17345395101703964*gelu(-0.9402956128299329 + 0.19522910845028774*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + -0.5236489579886725*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + -0.09802794343243137*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + -0.5359271533353085*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + 0.5362141752531211*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + -0.24637355767613478*gelu(-0.842982128122256 + 0.6833459994096582*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + 0.15560340828874075*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + 0.5368387122163134*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + -0.6086259049043772*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + -0.5403730914454505*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))) + 0.4183638306902213*gelu(0.6067063734003204 + 0.8661978259504042*gelu(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])) + -0.06088102390312189*gelu(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])) + 0.9836506901320718*gelu(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])) + 0.7952476421961925*gelu(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])) + -0.6703006410599177*gelu(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Functions for AR(k) model

# Generate data from autoregression model with specified noise
# n = sample size (number of time points)
# a = autoregression coefficients
# f_noise(t) = noise at time t
function generate(n,a,f_noise)
    k = length(a)
    x = zeros(n+k)
    for t = k+1:n+k
        x[t] = dot(x[t-1:-1:t-k],a) + f_noise(t)
    end
    x = x[k+1:n+k]
    return x
end

# Log of robustified marginal likelihood for AR(k) model
# x = data
# s = stddev of data, about the mean function
# k = order of model (# of coefs)
# sa = stddev of Normal priors on coefs a_i
# zeta = raise lik to this power
function log_marginal(x,k,s,sa,zeta)
    n = length(x)

    M = zeros(k,k)
    for i = 1:k, j = 1:k
        r = 0.0
        for t = (max(i,j)+1):n
            r += x[t-i]*x[t-j]
        end
        M[i,j] = r/s^2
    end
    v = zeros(k)
    for i = 1:k
        r = 0.0
        for t = (i+1):n
            r += x[t]*x[t-i]
        end
        v[i] = r/s^2
    end

    Lambda = zeta*M + eye(k)/sa^2

    log_m = 0.5*zeta^2*dot(v,Lambda\v) - k*log(sa) - 0.5*logdet(Lambda) - zeta*0.5*n*log(2*pi*s^2) - zeta*0.5*dot(x,x)/s^2

    return log_m
end

max_stars_repo_name: jwmi/CoarsenedPosterior

max_issues_repo_name: jwmi/CoarsenedPosterior

max_forks_repo_name: jwmi/CoarsenedPosterior


text:
using BrainFlow

BrainFlow.enable_dev_logger(BrainFlow.BOARD_CONTROLLER)

params = BrainFlowInputParams()
board_shim = BrainFlow.BoardShim(BrainFlow.SYNTHETIC_BOARD, params)

BrainFlow.prepare_session(board_shim)
BrainFlow.start_stream(board_shim, 45000, "file://data.csv:w")
sleep(1)
BrainFlow.insert_marker(1.0, board_shim)
sleep(1)
BrainFlow.stop_stream(board_shim)
data = BrainFlow.get_current_board_data(256, board_shim)
BrainFlow.release_session(board_shim)


text:
@testset "Shuffle" begin
    for s in ["A", "C", "G", "T"]
        kmer = DNAKmer(s)
        @test kmer === shuffle(kmer)
    end

    function count(kmer)
        a = c = g = t = 0
        for x in kmer
            a += x == DNA_A
            c += x == DNA_C
            g += x == DNA_G
            t += x == DNA_T
        end
        return a, c, g, t
    end

    for k in 1:64, _ in 1:10
        kmer = rand(DNAKmer{k})
        @test count(kmer) == count(shuffle(kmer))
        if k ≥ 30
            @test kmer != shuffle(kmer)
        end
    end
end


text:
RecipesBase.@recipe function f(mtg::MultiScaleTreeGraph.Node;
    mode = "2d",
    node_color = :black,
    colormap = :viridis,
    edge_color = nothing,
    color_missing = RGBA(0, 0, 0, 0.3)
)

    if edge_color === nothing
        edge_color = node_color
    end

    colormap = get_colormap(colormap)
    df_coordinates, node_color_in, edge_color, text_color =
        mtg_XYZ_color(mtg, node_color, edge_color, colormap, color_missing = color_missing)

    x = df_coordinates.XX
    y = df_coordinates.YY
    z = df_coordinates.ZZ

    for i in 2:size(df_coordinates, 1)
        x2 = [df_coordinates.XX_from[i], df_coordinates.XX[i]]
        y2 = [df_coordinates.YY_from[i], df_coordinates.YY[i]]
        z2 = [df_coordinates.ZZ_from[i], df_coordinates.ZZ[i]]
        edge_col = edge_color[i][1] # Plot is not compatible with line gradients

        RecipesBase.@series begin
            label := ""
            seriescolor := edge_col
            if mode == "2d"
                seriestype := :line
                x2, y2
            else
                seriestype := :line3d
                x2, z2, y2
            end
        end
    end

    hover_arg = string.(
        "name: `node_", df_coordinates.id,
        "`, link: `", df_coordinates.link,
        "`, symbol: `", df_coordinates.symbol,
        "`, index: `", df_coordinates.index, "`",
        if Symbol(node_color) in get_attributes(mtg)
            string.(", $node_color: ", df_coordinates[:, node_color])
        else
            ""
        end
    )

    RecipesBase.@series begin
        label := ""
        palette := colormap
        color := node_color_in
        # marker_z := node_color
        colorbar_entry := false
        hover := hover_arg
        if mode == "2d"
            seriestype := :scatter
            x, y
        else
            seriestype := :scatter3d
            x, z, y
        end

    end
end

"""
    plot(opf::MultiScaleTreeGraph.Node; kwargs...)
    plot!(opf::MultiScaleTreeGraph.Node; kwargs...)

Make a diagram of the MTG tree, paired with a `Plots.jl` backend.

See also [`diagram`](@ref) for the same plot with a `Makie.jl` backend.

## Attributes

- `mode = "2d"`: The mode for plotting, either "2d" or "3d"
- `node_color = :black`: the node color, can be a color or any MTG attribute
- `edge_color = node_color`: same as `node_color`, but for the edges
- `colormap = :viridis`: the colormap used for coloring
- `color_missing = RGBA(0, 0, 0, 0.3)`: The color used for missing values

# Examples

```julia
# import Pkg; Pkg.add("PlotlyJS")
using Plots, PlantGeom
plotlyjs()

file = joinpath(dirname(dirname(pathof(PlantGeom))),"test","files","simple_plant.opf")
# file = joinpath(dirname(dirname(pathof(PlantGeom))),"test","files","coffee.opf")

opf = read_opf(file)

plot(opf, node_color = :Length)
```
"""
RecipesBase.plot!, RecipesBase.plot


text:
################################################################################
# ProblemSize
################################################################################

mutable struct ProblemSize{SVu,SVx,SVz}
    N::Int
    n::Int
    m::Int
    p::Int
    ni::Vector{Int}
    mi::Vector{Int}
    pu::SVu
    px::SVx
    pz::SVz
    S::Int
end

function ProblemSize(N::Int, model::AbstractGameModel)
    n = model.n
    m = model.m
    p = model.p
    S = n*p*(N-1) + m*(N-1) + n*(N-1)
    return ProblemSize(
        N,
        model.n,
        model.m,
        model.p,
        model.ni,
        model.mi,
        model.pu,
        model.px,
        model.pz,
        S
        )
end

import Base.==
function (==)(p1::ProblemSize, p2::ProblemSize)
    out = true
    for name in fieldnames(ProblemSize)
        out &= getfield(p1, name) == getfield(p2, name)
    end
    return out
end


text:
export
    ispentagonal,
    nthpentagonal, npentagonal,
    allpentagonal, somepentagonal, exactpentagonal

ispentagonal(n::Integer) = begin
    t = 24*n + 1
    r = isqrt(t)
    (r^2 == t) && (r % 6 == 5)
end

nthpentagonal(n::Integer) = (n * (3n-1)) ÷ 2
npentagonal(n::Int, T::Type = Int) = collect(exactpentagonal(n, T))

allpentagonal(T::Type = Int) = PentagonalIterator{T}()
somepentagonal(xmax::T) where {T<:Integer} = @pipe allpentagonal(T) |> takewhile(x -> x <= xmax, _)
exactpentagonal(n::Int, T::Type = Int) = @pipe allpentagonal(T) |> take(_, n)

struct PentagonalIterator{T<:Integer}
end

Base.start(::PentagonalIterator{T}) where {T<:Integer} = zero(T), one(T)
Base.next(::PentagonalIterator{T}, state) where {T<:Integer} = begin
    s, n = state
    s += 3*n - 2
    s, (s, n + 1)
end
Base.done(::PentagonalIterator, _) = false

Base.eltype(it::PentagonalIterator) = Base.eltype(typeof(it))
Base.eltype(::Type{PentagonalIterator{T}}) where {T<:Integer} = T

Base.iteratorsize(::PentagonalIterator) = Base.IsInfinite()