text:
# This file is derived from julia iobuffer.jl. License is MIT: https://julialang.org/license

## work with String as UInt8 buffer via I/O primitives and share with tokens/substrings ##

# Stateful string

const MaybeIO = Union{Nothing, IO}

"a constant for a default allocation limit in IOShared and other resizable array structures"
const DEFAULTLIMIT = 1024%UInt32

"value indicating no mark ist set"
const NOMARK =typemax(UInt32)

"""
IO buffer which can share content with SubString and Token instances.

Reading of tokens or substrings can be done without heap allocations -
such returned items simply share the internal buffer, IOShare instance
will "copy on write" to protect shared data.

An IOShared instance operates in one of the following, mutually exclusive, modes:

1) *IO queue:* like IOBuffer, readable and writable. with some additional
text processing methods giving a "mutable string" behavior.

2) *buffered input stream:* a data source (another IO instance) is attached.
 When a read operation requests more bytes than available in the internal
 buffer, data from attached source is appended, and already read content
 is removed if not blocked by a mark. IOShared is usually only read.
 Writing is possible, but unusual: if will effectively insert data into the
 input stream at the current write offset position, which is under control
 of IOShare. More realistic is insertion at the read position, e.g. to
 replace some include notation by the content to include.

3) *buffered output stream:* a data sink (another IO instance) is attached.
 When a write operation requests more bytes than free space is available
 in the internal buffer, all content is written to attached data sink and
 then removed if not blocked by a mark. IOShared is usually only written.
 Read access is permitted, but will effectively delete data from the
 output stream. Similar to 2), insert/delete operations are more relevant
 to manipulate the stream befor it gets finally written to the sink.
 Use mark to ensure that content behind mark is not written to data sink,
 but remove mark as soon as possible.

 All operating flavours maintain two internal offsets, marking the currently
 used part of the internal buffer.
 One is the current read position, the other the current write position.

 IOShared queues allow (intermittent reads and writes).
 Current content is only changed by explicit content operations (read,
 write, resize and some more). Offsets into content, however, can change, due
 to reorganization. Already read content (before readofs) can be removed, reducing
 offsets of remaining content.

 IOshared buffered streams have only a possibly small window of the whole
 data in its internal buffer. Reading from a buffered output stream is
 tecnically possible, but will effectively delete data  the stream: read data is no more written
 to the data sink. Similarly, writing to a buffered input stream will insert
 data into the stream, which is usually not intended. it removes data before  data in data sink,
 operations can cause flushing of current content to data sink. With attached
 data source, read operations may trigger appending data from source to
 internal buffer. Time and size of those implicit content changes are under
 control of the IOBuffer implementation.

 You can govern implicit content changes by use of the mark mechanism: setting the
 mark to a content position, effectively prevents implicit content changes
 between mark position and write position.  However content before the mark
 position may get removed by implicit changes, moving all offsets.
 To also prevent this, set mark position to 0. This and only this guarantees that
 content offsets will be changed by explicit operations, only.

In addition to the operating mode, every IOShared can act as a *heap for
flyweight tokens:* vectors of flyweight tokens can be registered to the
IOShared, using it for all content which needs a buffer.
[`Tokenvector`](@ref) is based on this feature. Compared to Vector{Token},
the overhead per token is halved, and memory locality is
improved for operations on the vector. Heap space is allocated from the end
of the internal buffer, reducing the *limit* offset. This way, it does not
conflict with the operating modes for reading/writing.

IOShared has an internal garbage collection for its heap: if it runs out
of free memory in its buffer, it allocates a new buffer, and copies
living heap entries to a contiguous memory block. If no other objects
keep a reference to the old buffer, it gets (julia) garbage collected.
Sharing content with SubString, HToken and BToken will prevent garbage
collection of the old buffer. This can have
effects similar to memory leaks in classical C code: a small SubString/Token
blocks a possibly large buffer. For excellent memory efficiency, avoid
sharing only small parts of an IOShared with SubString/HToken/BToken
structures. With TokenVector{VToken} and VToken elements, you can leave all
content in one IOShared buffer, and avoid most content copying.
A typical use case is a parse process, which keeps only a small percentage of
source text in tokens.

"""
mutable struct IOShared <: IO
    buffer :: String # PRIVATE! memory with byte data.
    readofs :: UInt32 # read position offset / number of consumed bytes
    writeofs :: UInt32 # write (append) position offset / end of defined data
    shared :: UInt32 # offset in buffer behind last shared byte (0: nothing shared)
    mark :: UInt32 # offset of data not to flush/forget. typemax: no active mark
    limit :: UInt32 # buffer limit for read/write operations, memory above is used by put function
    preferredlimit :: UInt32 # limit in case of reallocation due to sharing
    registered :: Vector{HybridFly} # true: mark is always 0 (offsets are guaranteed not to change)
    io::MaybeIO # nothing or if ioread: source for fillup, else target for flush
    eolRemoved :: Int32 # <0: undefined. else: no. of eol bytes which were removed from buffer by reorg-s
    eol :: UInt8 #  last byte of an end-of-line-sequence (default: 10, change to 13 for Mac text files)
    ioread :: Bool # true: io is source for fillup. false: io is target for flush or nothing
    compressOnPut :: Bool # true: put does lookup content, share instead of append if found
    function IOShared(io::MaybeIO,ioread::Bool,limit::UInt32)
        ioread && io===nothing && error("ioread must be false if io is nothing")
        z = zero(UInt32)
        s = _string_n(limit)
        new(s,z,z,z,NOMARK,limit,limit,Vector{HybridFly}(),io,-1,10,ioread,false)
    end
    function IOShared(offset::UInt32, size::UInt64, s::String)
        @boundscheck checkrange(offset,size,s)
        size32= UInt32(size) # check on overflow necessary: cannot handle too long strings
        preferredlimit = max(size32,DEFAULTLIMIT) # well ... just a very simple guess
        limit = UInt32(ncodeunits(s)) # check on overflow necessary: cannot handle too long strings
        new(s,offset,offset+size32,offset+size32,NOMARK,limit,preferredlimit,Vector{HybridFly}(),nothing,-1,10,ioread,false)
    end
end

IOShared(io::MaybeIO,ioread::Bool) = IOShared(io,ioread,DEFAULTLIMIT)

IOShared(limit::UInt32=DEFAULTLIMIT) = IOShared(nothing,false,limit)

IOShared(s::SubString{String}) = IOShared(UInt32(s.offset),s.ncodeunits%UInt64,s.string)

IOShared(s::String) = IOShared(zero(UInt32),usize(s),s)

IOShared(t::BToken) = IOShared(t.offset,t.len,t.buffer)

IOShared(t::AbstractString) = IOShared(BToken(t))

## internal API

#"offset boundscheck, including "
#Base.checkbounds(::Type{Bool}, io::IOShared, ofs::UInt32) = io.readofs <= ofs <= io.writeofs

"""
    _share(io::IOShared,share::UInt32) :: String

Return a buffer reference. Like take!, this function gives you efficient access to
the whole content of io, but unlike take!, io continues using buffer also for write
operations.

That buffer is typed as String, but somehow violates String immutability:
only the first share bytes in the returned buffer are guaranteed not to change.
If io needs to change its buffer in the first share bytes, it allocates a
new buffer and copies needed content to it.

Of course, the returned String must be treated immutable - writing to it
via pointer operations may result in system crashes.

_share is regarded a private function. To obtain the text content of an IOShared,
use read operations or  [`get`](@ref) .
"""
function _share(io::IOShared,share::UInt32)
    if (io.shared<share)
        io.shared = share
    end
    return io.buffer
end

"share the whole internal buffer. ATTENTION: may contain uninitialized data"
_share(io::IOShared) = _share(io,usize32(io.buffer))

"currently loaded content size in bytes. Changes e.g. on fill, flush, read, write"
usize(io::IOShared) = usize32(io) % UInt64

usize32(io::IOShared) = io.writeofs-io.readofs

"""
    unsafe_put(heap::IOShared, t::Token{F})

put contents of t into the heap of an IOShared and return a [`FlyToken`](@ref) referencing
this IOShared with the same content.

Is marked unsafe, because the following preconditions are not checked (task of caller):

 * enough free space in heap
 * token is not direct
 * token size > 0

Method is used internally when storing in a registered [`TokenVector`](@ref),
and by realloc.

Only module-intern use, from outside use [`put`](@ref), which ensures preconditions.

If heap.compressOnPut is true, a search for the contents in heap is performed,
found contents is referenced. This can reduce memory consumption
but increases CPU use. Default is false, change it with compressOnPut(heap,true)
"""
function unsafe_put(heap::IOShared, t::Token{F}) where F
    size = usize32(t)
    if heap.compressOnPut
        # try to find in in heap space
        ofs = locate(t,heap.limit,usize(heap.buffer)-heap.limit,heap.buffer)
        if ofs<typemax(UInt32)
            return (t.fly & !OFFSET_BITS) + ofs
        end
    end
    # not found: append data in heap
    # ensure_writeable(heap,size) # already done by caller
    heap.limit -= size
    ofs = heap.limit
    b = heap.buffer
    bt = t.buffer
    GC.@preserve b bt unsafe_copyto!(pointer(b)+ofs, pointer(bt)+offset(t), size)
    return (t.fly & !OFFSET_BITS) + ofs
end

"""
    put(heap::IOShared, s, ::Token{F})

return a [`Token`](@ref) with a string representation of s. Returned token is
direct or uses heap as its backing buffer for its content.

s can be anything for which a Token{F} constructor exists,
in particular AbstractString, AbstractToken, Real, Integer, Bool.
If s isa AbstractToken, its category is preserved,
for any other type of s, a category derived from the type of s is used,
e. g. T_INT for integers or T_TEXT for strings.

put stores contents in a way not interfering with any read and write access on heap.
The IOShared parameter is named heap to indicate how it works: put uses a "heap area"
in the internal buffer at the physical end of the buffer, growing towards
lower adresses. heap area is share-protected, not available for read and write access,
and preserved on a buffer reallocation.

If s has its content already in heap (s is a token or substring using heap-s buffer)
or returned token is direct, heap is not changed.

If heap.compressOnPut is true, a search for the content is performed,
found content is referenced, avoiding content move and memory consumption.
The  price is increased CPU use: put is not O(1) but O(total used buffer size),
in this case. Default is false, change it with compressOnPut(heap,true)
"""
put(heap::IOShared, s, ::Type{Token{F}}) where F = put(heap,Token{F}(s))

function put(heap::IOShared, t::Token{F}) where F
    size = usize32(t)
    if !isdirect(t) && size>0 && t.buffer !== heap.buffer
        if F <: HybridFly && size<=MAX_DIRECT_SIZE
            return HToken(DirectFly(t))
        end
        #= search here makes no sense, is performed by unsafe_put iff heap.compressOnPut is true
        # already in heap?
        ofs = locate(t,heap.limit,usize(heap.buffer),heap.buffer)
        if ofs<typemax(UInt32)
            return Token{F}((t.fly & !OFFSET_BITS) + ofs,heap.buffer)
        end
        # or in content area?
        ofs = locate(t,0,heap.writeofs,heap.buffer)
        if ofs<typemax(UInt32)
            return Token{F}((t.fly & !OFFSET_BITS) + ofs,_share(heap,ofs+size))
        end
        # not found
        =#
        ensure_writeable(heap,size)
        f = unsafe_put(heap, t)
        t = Token{F}(f,heap.buffer) # we must not use share: all above offset heap.limit is implicitly shared
    end
    return t
end

"set heap.compressOnPut for subsequent put calls"
function compressOnPut(heap::IOShared, compress::Bool)
    heap.compressOnPut = compress
    return nothing
end

"compute tuple () no. of lines before offset ofs, offset of begin of line containing position ofs)"
function _lineofs(io::IOShared, ofs::UInt32) :: Tuple{UInt32,UInt32}
    lines = zero(UInt32)
    ofs = 0%UInt32 # position of begin of line
    while (oeol=locate(io.eol,ofs,(io.readofs-ofs)%UInt64,io.buffer))<io.readofs
        ofs = oeol+one(UInt32)
        lines += one(UInt32)
    end
    return (lines, ofs)
end

"""
    _eolAdjust(io::IOShared, preservedOfs::UInt32) ::UInt32

called if a reorg function wants to remove data section 0..preservedOfs,
before content destruction.

if eol processing is not activated (default), return preservedOfs

Otherwise, counts all EOL markers in buffer up to offset preservedOfs and
add io to eolRemoved. Return the offset of the begin of the last text line
to be removed. This will enforce that buffer contains the whole line, e.g.
to generate line listings in error situations.
"""
function _eolAdjust(io::IOShared, preservedOfs::UInt32)
    if io.eolRemoved>=0

    end
    return preservedOfs
end

"""
    linepos(io::IOShared) :: Pair{Int,Int}

return -1 => -1 (undefined) or lineNo =>index of the current read position in line.
lineNo and index are 1-based (first byte of a file has (1,1) as its linepos)

text line structure is undefinded by default.
To enable line number tracking, call [`tracklines`](@ref).

linepos is not cheap: it scans for end-of-line markers throughout the
internal buffer, which is quite expensive for large buffers. It is
intended for error reporting in parsers.

Consider reading complete lines and counting lines in application code,
if you frequently need line number information.

"""
function linepos(io::IOShared) :: Pair{Int,Int}
    io.eolRemoved<0 && return -1 => -1
    lines, ofs = _linepos(io,io.readofs)
    return lines+io.eolRemoved+1 => io.readofs - ofs + 1
end

"""
    tracklines(io::IOShared, eol::UInt8, linesRemoved::Int)

activate line number tracking for *io*, using *eol* as marker byte for an end-of-line sequence
and initialize the number of already removed lines with *linesRemoved*. Call it immediately
after constructing io with *linesRemoved*==0.

Line number tracking has its price: it requires scanning for end-of-line markers in
possibly huge text buffers. Therefore, it is disabled by default.
"""
function tracklines(io::IOShared, eol::UInt8, linesRemoved::Int)
    io.eol = eol
    io.linesRemoved = linesRemoved
    return nothing
end

"""
move memory within an IOShared buffer: insert/delete at ofs

Method is regarded private and intended to be used within IOShared, only.
External API is function resize.

ofs is internal offset into io-s buffer.
PRECOND: limit large enough, no share at ofs, no delete beyond io.writeofs
ofs <= io.writeofs, if resize<0: ofs-resize <= io.writeofs

readofs, writeofs and mark are adjusted in io.
No reallocation - if resize would violate sharing, an error is thrown.

Resize of IOShared having registered tokens can fail and is not recommended!

Resize returs true on a sucessful operation, else false.
"""
function _resize(io::IOShared, ofs::UInt32, resize::Int)
    @boundscheck begin
        (io.shared <= ofs # never chance shared contents
        && ofs <= io.writeofs # not beyond end of content
        && ofs <= io.writeofs + resize  # no delete beyond end of content (resize<0)
        && io.writeofs+resize <= io.limit # enough space (relevant if resize>0)
        ) || error("preconditions violated for resize by $resize at $ofs in $io")
    end
    #= NO!! unnecessary!! all registered tokens are share-protected!!
    if ofs<io.writeofs
        # we need to adjust registered tokens referencing nonheap buffer memory
        for v in io.registered
            for i in 1:length(v)
                t = v[i]
                s = t.len
                if !isdirect(t) && s>0
                    # we have a heap based content, test for need to relocate it
                    o = offset(t)
                    if o >= ofs
                        # adjust offset if necessary
                        if o < io.writeofs
                            v[i] += resize
                        end
                    else
                        if o+s>ofs
                            # not nice - token would be destroyed and must be relocated
                            if io.writeofs+s>io.limit
                                return false  # we give up - not enough space to safely relocate heap token.
                            end
                            v[i] = hf(unsafe_put(io,BToken(bf(t),io.buffer)))
                        end
                    end
                end
            end
        end
    end
    =#
    s = io.buffer
    ptr = pointer(s)
    GC.@preserve s begin
        if resize>0 #insert is the most frequently used case
            ccall(:memmove, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
                  ptr+ofs+resize, ptr+ofs, io.writeofs-ofs)

        else #delete/noop
            resize == 0 && return nothing
            ccall(:memmove, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
                  ptr+ofs, ptr+ofs-resize, io.writeofs-ofs+resize)
        end
    end
    io.writeofs += resize
    ofs<io.readofs && (io.readofs += resize)
    io.mark != NOMARK && ofs<io.mark && (io.mark += resize)
    return nothing
end

"""
allocate a new buffer and preserve content, resizing at offset ofs.

Method tries to flush first, to reduce allocation as well as
the amount of bytes to copy.

io is changed, adjusted ofs is returned.

Inserted bytes (resize>0) are uninitialized.

ofs precondition: io.readOfs <= ofs <= io.writeOfs

realloc is used internally, external use is
not recommended.

"""
function _realloc!(io::IOShared, ofs::UInt32, resize::Int) ::UInt32
    flush(io) # write to sink to decrease used area
    preservedOfs = min(io.readofs, io.mark) # begin of protected content which must be copied. end is io.writeofs
    if io.eolRemoved >= 0
        # we have EOL tracking. reduce preservedOfs such that is the begin of a text line.
        # update io.eolRemoved with number of lines below preservedOfs
        (lines, preservedOfs) = _lineofs(io,preservedOfs)
        io.eolRemoved += lines
    end
    need = io.writeofs - preservedOfs + resize # minimum size of new buffer (without heap).
    @boundscheck begin
        (io.readofs <= ofs
        && ofs <= io.writeofs
        && need >= 0 && need <= typemax(UInt32)
        && (ofs-resize <= io.writeofs)
        ) || boundserror(io,ofs,resize)
    end
    # compute size of all tokens which must be put into heap because content is not preserved
    relocSize = 0 # size of all tokens not already in heap
    heapUsed = 0 # size of tokens already in heap
    oldLimit = io.limit
    for v in io.registered
        for t in v
            s = t.len
            o = offset(t)
            if !isdirect(t)
                if o<oldLimit
                    # we will relocate all registered tokens not already in heap
                    # reason: we want shared==0 on return.
                    relocSize += s
                else
                    heapUsed += s
                end
            end
        end
    end
    heapSize = usize(io.buffer)-oldLimit
    need += relocSize + min(heapSize,heapUsed) # keep old heap if already compressed
    #
    # new limit: several considerations.
    # (1) at least, it must be >= need
    # (2) should not be smaller than preferredlimit
    # (3) ensure that write(UInt8) is O(1) ==> enlarge need such that free space is O(need)
    newBuffersize = max(io.preferredlimit,UInt32((need*3)>>1)) # 50% more than need
    oldbuffer = io.buffer
    newbuffer = _string_n(newBuffersize)
    io.buffer = newbuffer
    io.limit = newBuffersize
    keptHeapOffset = typemax(UInt32) # all offsets above are in old heap and only to adjust
    keptHeapDelta = 0
    GC.@preserve oldbuffer newbuffer begin # for safety. might be unnecessary
        if heapSize<heapUsed
            # keep old heap as contiguous memory block.
            # move it now, so that relocation can search in new heap.
            ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
              pointer(newBuffer)+newBuffersize-heapSize, pointer(oldBuffer)+oldLimit, heapSize)
            io.limit -= heapSize
            keptHeapOffset = oldLimit
            keptHeapDelta = io.limit-oldLimit
        end
        # treatment of registered tokens: adjust offsets, maybe move content
        for v in io.registered
            for i in 1:length(v)
                t = v[i]
                if !isdirect(t) && t.len>0
                    o = offset(t)
                    if o >= keptHeapOffset
                        v[i] = t+keptHeapDelta
                    else
                        v[i] = hf(unsafe_put(io,BToken(t,oldbuffer)))
                    end
                end
            end
        end
        src = pointer(oldbuffer)
        newbuffer = io.buffer
        dst = pointer(newbuffer)
        count = (ofs - preservedOfs) %UInt # size of reserved content before ofs
        ## copy first part: preservedOfs..ofs
        ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
            dst, src+preservedOfs, count)
        # copy 2nd part: ofs..io.writeofs
        if resize<0
            count2 =  (io.writeofs-ofs+resize) %UInt
            if count2>0
                ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
                    dst+count, src+ofs-resize, count2)
            end
        else
            count2 =  (io.writeofs-ofs) %UInt
            if count2>0
                ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
                    dst+count+resize, src+ofs, count2)
            end
        end
        ofs -= preservedOfs
        io.readofs -= preservedOfs
        io.mark != NOMARK && (io.mark -= preservedOfs)
        io.writeofs += resize-preservedOfs
    end
    io.shared = 0
    return ofs
end

"""
    shrink!(io:IOShared, reserve::UInt32, unshare::Bool = true)

if unused memory of io exceeds reserve, reallocate internal buffer of io such that exactly reserve bytes are unused by content or heap.

If unshare is true, all registered flyweight tokens are moved to heap, which can increase used memory. But it ensures
that content area is unshared on return.

No memory optimization is performed - that should be done before.
"""
function shrink!(io::IOShared, reserve::UInt32, unshare::Bool = true)
    if io.limit - io.writeofs > reserved
        if unshare
            oldLimit = io.limit
            for v in io.registered
                for t in v
                    s = t.len
                    o = offset(t)
                    if !isdirect(t)
                        if o<oldLimit
                            # we will relocate all registered tokens not already in heap
                            # reason: we want shared==0 on return.
                            put(io,BToken(t,io.buffer))
                        end
                    end
                end
            end
            io.shared = 0
        end
        # compute size of all tokens which must be put into heap because content is not preserved
        if io.limit - io.writeofs > reserved
            toshrink = io.limit - io.writeofs - reserved
            heapsize = usize32(io.buffer)-io.limit
            newBuffersize = usize32(io.buffer)-toshrink
            oldbuffer = io.buffer
            newbuffer = _string_n(newBuffersize)
            io.buffer = newbuffer
            io.limit -= toshrink
            GC.@preserve oldbuffer newbuffer begin # for safety. might be unnecessary
                # keep old heap as contiguous memory block.
                ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
                pointer(newBuffer)+newBuffersize-heapSize, pointer(oldBuffer)+oldLimit, heapSize)
                # keep content area as is
                ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
                    newbuffer, oldbuffer, io.writeofs)
            end
        end
    end
    return nothing
end

"""

    modify!(io::IOShared, ofs::UInt32, resize::Int)

Ensure that io.buffer is writeable starting at ofs, and resize current content at ofs.

If resize is negative, deletion is restricted to the current content.
If io.shared>ofs, or if io.limit < ofs+resize, a new buffer is allocated.

Buffer allocation may change offsets, i.e. io.readofs, io.writeofs and ofs.
Instance variables are adjusted, the adjusted value of ofs is returned by modify.

In case of an insert/append, content at ofs..ofs+resize is uninitialized.

Method is used package-internally, but also exported as central low level resize method.

Common use is appending, i.e. ofs==io.writeofs, resize>0.
Any other use makes things complicated if a data sink or data source is assigned to io.
"""
function modify!(io::IOShared, ofs::UInt32, resize::Int)
    if io.shared<=ofs
        # good chances for a noop
        if io.writeofs+resize<=io.limit
            # insert/delete
            if ofs==io.writeofs && resize >=0
                # normal append: expected to be the most frequent case
                io.writeOfs += resize
                return ofs
            end
            _resize(io,ofs,resize)
            return ofs
        end
        # not enough space at end - can we reuse space at beginning?
        if io.eolRemoved<0
            # we require no eol tracking, because it does not allow resize work with 0 < resize pos < readofs
            flush(io)
            preservedOfs = min(io.readofs,io.mark)
            delta = io.shared%Int - preservedOfs # bytes deleteable at shared if <0
            if io.writeofs+resize+delta <= io.limit
                # removing obsolete content after shared is enough
                _resize(io,io.shared,delta) # remove free space at beginning
                ofs += delta
                _resize(io,ofs,resize)
                return ofs
            end
        end
    end
    # nothing but realloc helps ...
    return _realloc(io,ofs,resize)
end

function ensure_writeable(io::IOShared, count::UInt32)
    if io.shared>io.writeofs || io.writeofs+count>io.limit
        _realloc(io,io.writeofs,count%Int)
        io.writeofs -= count # we only ensured space for writing. write operation will change offset
    end
    return nothing
end

function Base.ensureroom(io::IOShared, nshort::Int)
    ensure_writable(io, UInt32(nshort))
    return io
end

"""
    ensure_readable(io::IOShared, count::UInt64)::Bool

ensure that count bytes can be read from io or throw an error.
"""
function ensure_readable(io::IOShared, count::UInt32)
    fillup(io,count)
    usize32(io) >= count && return nothing
    error("read request for $count bytes fails - end-of-data")
end

"""
    fillup(io::IO,count::UInt32)

try to ensure that count bytes are readable from io.
In contrast to  [`ensure_readable`](@ref), no error
is thrown if not successful.

if io has an attached data source, and usize(io)<count,
fillup will try to put additional count-usize(io) bytes
from data source into buffer of io.

n, the effectively number of bytes which are transferred from data source,
may vary a lot:

 * n==0 (because nothing is available),

 * 0 < n < count-usize(io) (because only n bytes are available)

 * n == count-usize(io) (very unlikely)

 * n > count-usize(io) (data sink has more than n bytes,
   and we can put more than count bytes into io without reallocation)

fillup is called internally on read requests to ensure enough
bytes are available in io-s buffer.

Search operations on IOBuffer never do fillup internally. To get
valid search results, you must call fillup before searching with
a count so large that match results beyond count bytes are
irrelevant.

NOOP if io.ioread==false
"""
function fillup(io::IOShared, count::UInt32)
    if io.writeofs-io.readofs < count && io.ioread
        count -=  io.writeofs-io.readofs # what is additionally requested
        avail = min(bytesavailable(io.io), typemax(UInt32)) % UInt32
        count > avail && (count = avail)
        # count is now fill request, reduced if less is available
        ensure_writable(io,count)
        # can we fill more than requested?
        if io.limit-io.writeofs > count
            count = min(avail, io.limit-io.writeofs)
        end
        # finally, do fill
        s = io.buffer
        ptr = pointer(s)
        GC.@preserve s unsafe_read(io.io,ptr+io.writeofs,count%UInt)
        io.writeofs += count
        # TODO: check if we have reached EOF and close io in that case?
    end
    nothing
end

locate(what::Union{AbstractString,AbstractChar},io::IOShared) = locate(BToken(what), io.readofs,usize(io),io.buffer)

locate(what::UInt8,io::IOShared) = locate(what, io.readofs,usize(io),io.buffer)

#TODO an best alg anpassen.
"""
    locate(pool::IOShared,ofs::UInt32,size::UInt64,s::String)

Search in pool for a byte sequence of size bytes, beginning at
offset ofs (0-based!) in String s.

Return the offset (0-based) in pool.buffer or typemax(UInt32) (not found).

TODO use boyer-moore search if size and usize(pool) are 'large'.
TODO search for last 32 bytes if size>32 - that maximizes skip distance
"""
function locate(pool::IOShared,ofs::UInt32,size::UInt64,s::String)
    if size == 0
        return pool.readofs
    end
    if size == 1
        return locate(pool, codeunit(s,ofs+1))
    end
    s_ptr = pointer(s)
    s_end = s_ptr +size
    buffer = pool.buffer
    p_ptr = pointer(buffer)+pool.readofs
    p_size = usize(pool)
    GC.@preserve buffer s begin
        cs_last = unsafe_load(s_ptr+size-1)
        if size == 1
            # single byte search
            return locate(io, cs_last)
        end
        # TODO: classical brute force search for very small size eg <=3 or 7
        # TODO: search with boyer-moore for larger sizes eg>15
        # search adopted from in _searchindex(s::ByteArray, t::ByteArray, i::Integer)
        # build minihash union
        bloom = 0%UInt64 # ORed minihashes s[1..size-1]
        skip = size # skip if last byte matches, some other byte not
        j = size-1
        # build bloom mask and maximal skip
        while (j-=1) >= 0
            cs = unsafe_load(s_ptr+j)
            bloom |= minihash(cs)
            if cs==cs_last
                skip = min(skip,size-j-1)
            end
        end
        # now bloom is ORed hash of all bytes in s except the last one
        # minihash(byte) & bloom == 0 implies: byte does not occur in s,
        # except last byte.
        # loop for all possible matches of last byte: i=size-1..p_size-1
        i = size-1
        while i < p_size
            cp = unsafe_load(p_ptr+i) # last byte in current match candidate
            #skip = 0 # we have not determined skip distance
            if cp==cs_last
                ## last byte matches. test successive backwards up to i-size+1
                j = size-1
                p0 = p_ptr+i-j # pointer to offset 0 in current match
                while (j-=1) >= 0
                    if unsafe_load(p0+j) != unsafe_load(s_ptr+j)
                        break
                    end
                end
                if j<0
                    #match
                    return i-(size-1) # offset of 1st byte
                end
                # no match.
                i += skip
            elseif bloom & minihash(cp) == 0
                # cp does not occur in s
                i += size
            else
                # default skip: 1 byte
                i += 1
            end
        end
    end
end

"fast search in buffer (no fillup!)"
function Base.occursin(delim::UInt8, io::IOShared)
    buf = io.buffer
    q = GC.@preserve buf ccall(:memchr,Ptr{UInt8},(Ptr{UInt8},Int32,Csize_t),pointer(buf)+io.readofs,delim,usize(io))
    return q != C_NULL
end

Base.match(r::Regex, io::IOShared) = match(r,substr(io))

#= unfertig!!!
function Base.readuntil(io::IOShared, delim::UInt8; keep::Bool=false)

    rel = 0 # relative offset
    while true
        fillup(io,rel+1)
        if io.readofs+rel >= writeofs
            error("Delimiter not found in data")
        end
        if codeunit(io.buffer)

    end
    lb = 70
    A = StringVector(lb)
    nread = 0
    nout = 0
    data = io.data
    for i = io.readofs : io.writeofs
        @inbounds b = data[i]
        nread += 1
        if keep || b != delim
            nout += 1
            if nout > lb
                lb = nout*2
                resize!(A, lb)
            end
            @inbounds A[nout] = b
        end
        if b == delim
            break
        end
    end
    io.readofs += nread
    if lb != nout
        resize!(A, nout)
    end
    A
end
=#

## Base API

"""
if not io.ioread, write all unmarked current contents to data sink.

written content is removed from this.
no buffer reorganization, but flushing can increase
re-use of an unshared buffer
"""
function Base.flush(io::IOShared)
    if io !== nothing && !io.ioread
        endofs = min(io.mark,io.writeofs)
        if  endofs > io.readofs
            sswrite(io.io,io.readofs,endofs-io.readofs, io.buffer)
            io.readofs = endofs
        end
    end
    nothing
end

"return a new IOShared with same (shared) contents, but no mark, default preferredlimit"
function Base.copy(b::IOShared)
    b.io === nothing || error("cannot copy IOShared with attached source/sink")
    return IOShared(io.readofs,usize(io),io.buffer)
end

function Base.show(io::IO, b::IOShared)
    print(io,   "IOShared(  readofs=",b.readofs,
                ", writeofs=", b.writeofs,
                ", shared=", b.shared,
                ", limit=", b.limit,
                ", '")
    len = io.writeofs-io.readofs
    if len<=11
        print(io,SubString(buffer,readofs+1,readofs+len))
    else
        print(io,SubString(buffer,readofs+1,readofs+5),
        '~',SubString(buffer,size-4,size))
    end
    print(io,"')")
end

Base.skip(io::IOShared, delta::Integer) = seek(io,io.readofs+delta)

## read functions

function Base.unsafe_read(from::IOShared, p::Ptr{UInt8}, nb::UInt)
    #from.readable || throw(ArgumentError("read failed, IOBuffer is not readable"))
    ensure_readable(from,nb)
    b = from.buffer
    GC.@preserve b unsafe_copyto!(p, pointer(b)+from.readofs, nb)
    from.readofs += nb
    nothing
end

"""
    Base.read(io::IOShared, ::Type{BToken})

Read a Token in serialized format, sharing the buffer (no content copying).
This is a main benefit of IOShared over IOBuffer:
IOShared is still mutable, it keeps track of shared portions
and reallocates if shared portions have to be changed.

"""
function Base.read(io::IOShared, ::Type{Token{T}}) where T
    p = read(io,Packed31)
    size = t.len
    ensure_readable(io,size)
    if T <: HToken && size <= MAX_DIRECT_SIZE
        return HToken(read(io,p,DirectFly))
    end
    t = BufferFly(p) + io.readofs
    io.readofs += size
    @inbounds return Token{T}(t,_share(io,io.readofs))
end

"read bytes as a shared substring. Similar available as a peek method"
function Base.read(io::IOShared, ::Type{SubString}, size::UInt64)
    ensure_readable(io,size)

    t = BufferFly(p) + io.readofs
    # reference and skip string data in IOShared instance
    ofs = io.readofs
    io.readofs += size
    @inbounds return substring(ofs,size,_share(io,io.readofs))
end

"read all content in buffer as a SubString (no fillUp!)"
Base.read(io::IOShared, ::Type{SubString}) = read(io,usize(io))

function Base.read(from::IOShared, T::Union{Type{Int16},Type{UInt16},Type{Int32},Type{UInt32},Type{Int64},Type{UInt64},Type{Int128},Type{UInt128},Type{Float16},Type{Float32},Type{Float64}})
    nb = sizeof(T)
    ensure_readable(from,nb)
    b = from.buffer
    GC.@preserve b begin
        ptr::Ptr{T} = pointer(b)+from.readofs
        x = unsafe_load(ptr)
    end
    from.readofs += nb
    return x
end

# not efficient for IOShared - consider read(SubString)
"return up to requested count bytes as Vector"
function Base.read(src::IOShared, count::Integer = typemax(Int))
    size = min(count, bytesavailable(src))
    b = Vector{UInt8}(undef, size)
    readbytes!(src, b, size)
    return b
end

Base.read(from::IOShared, ::Type{Ptr{T}}) where {T} = convert(Ptr{T}, read(from, UInt))

function Base.read(from::IOShared, ::Type{UInt8})
    ensure_readable(from,1)
    b = from.buffer
    GC.@preserve b ret = unsafe_load(pointer(b)+from.readofs)
    from.readofs += 1
    return ret
end

#= inherit it from IO!
function Base.readbytes!(src::IOShared, b::Vector{UInt8}, count=length(b))
    # resize if necessary
    size = min(count, bytesavailable(src))
    if (size > length(b))
        resize!(b, size)
    end
    unsafe_read(src,pointer(b),size)
    return size
end
=#

## write functions

function Base.unsafe_write(to::IOShared, p::Ptr{UInt8}, nb::UInt)
    ensure_writeable(to,nb)
    b = to.buffer
    GC.@preserve b unsafe_copyto!(pointer(b)+to.writeofs, p, nb)
    to.writeofs += nb
    return nothing
end

function Base.write(to::IOShared, byte::UInt8)
    ensure_writeable(to,1)
    b = to.buffer
    GC.@preserve b unsafe_store!(pointer(b)+to.writeofs, byte)
    to.writeofs += 1
    return nothing
end

"write all available bytes (no waiting for async streams)"
function Base.write(to::IOShared, from::IO)
    # would work, but allocates intermediate Vector #  = write(to,read(from))
    b = to.buffer
    size = UInt32(bytesavailable(from))
    ensure_writeable(to,size)
    GC.@preserve b unsafe_read(from,pointer(b)+to.writeofs,size)
    to.writeofs += size
    return nothing
end

"write all available bytes in from's buffer (with fillup)"
function Base.write(to::IOShared, from::IOShared)
    to===from && error("write(io,io) not allowed for io::IOShared ")
    sswrite(to,from.readofs,available(from),from.buffer)
end

# size related operations

function Base.eof(io::IOShared)
    fillup(io,1%Uint32)
    return io.readofs == io.writeofs
end

function Base.bytesavailable(io::IOShared)
    ret = usize(io)%Int
    if io.ioread
        ret += bytesavailable(io.io)
    end
end

## offset related operations

Base.position(io::IOShared) = io.readofs %Int

"""
    seek(io::IOShared, n::Integer)

Seek moves read position given by position(io) to a new offset.

Please take into account that read and write operations may move
content within the internal buffer, offsets are usually not stable.

Use seek within currently defined content, given by the offsets
from position(io) up to position(io)+usize(io) directly before seek operation,
or use mark to "freeze" content above mark position, and seek
relative to mark position.

If mark is 0, all offsets are preserved except for explicit resize operations.
"""
function Base.seek(io::IOShared, n::Integer)
    # different from IOBuffer which projects on valid offsets, due to a REPL requirement
    (n < 0 || n > io.writeofs) && throw(ArgumentError("Attempt to seek outside IOShared boundaries."))
    io.readofs = n%UInt32
    return io
end

Base.seekstart(io::IOShared) = seek(io,0)

function Base.seekend(io::IOShared)
    io.readofs = io.writeofs
    return io
end

function Base.mark(io::IOShared)
    io.mark = io.readofs
end

Base.ismarked(io::IOShared) = io.mark != NOMARK

function Base.unmark(io::IOShared)
    ret = ismarked(io)
    io.mark = NOMARK
    return ret
end

function Base.reset(io::IOShared)
    io.mark == NOMARK && error("IOShared is not marked")
    io.readofs = io.mark
    io.mark = NOMARK
    return io.readofs
end

"""
    truncate(io::IOShared, n::Integer)

Seek operation for write position. truncate can move write position behind
current write position, this is equivalent to writing 0x00 bytes such that
write position becomes n.

If the new write position is smaller than the mark position,
mark position is removed.

n is a 0-based offset. It is declared Integer to be compatible with
the method signature in Base, value must be unsigned.

Warning: take care if io.io is not nothing
"""
function Base.truncate(io::IOShared, n::Integer)
    n < 0 && error("truncate failed, n must be ≥ 0, got $n")
    n<io.writeofs && length(io.registered)>0 && error("truncate not implemented for IOShared used as heap. ")
    pos = UInt32(n) # does check for 32 bit overflow
    if pos >io.writeofs
        ensure_writeable(io,pos-io.writeofs)
        for i in io.writeofs:pos
            write(io,0x00)
        end
    end
    io.writeofs = pos
    io.readofs>io.writeofs && (io.readofs = io.writeofs)
    io.writeofs < io.mark && io.mark != NOMARK && unmark(io)
    return io
end

"peek a byte at offset (aka 0-based index) in data of this"
function byte(s::IOShared, ofs::UInt32)
    ensure_readable(s,ofs+1)
    b = s.buffer
    ret = GC.@preserve b unsafe_load(pointer(b)+s.readofs+ofs)
    return ret
end

"lookahead 1 byte"
Base.@propagate_inbounds function Base.peek(from::IOShared)
    ensure_readable(from,1)
    b = from.buffer
    GC.@preserve b begin
        ptr = pointer(b)+from.readofs
        x = unsafe_load(ptr)
    end
    return x
end

"lookahead n bytes, returned as a BToken"
function Base.peek(io::IOShared, size::UInt64)
    ensure_readable(io,size)
    # reference string data in IOShared instance
    @inbounds return BToken(io.readofs,size,_share(io,io.readofs+size))
end

## other IO function overloads for IOShared

Base.isreadable(io::IOShared) = true
Base.iswritable(io::IOShared) = true

"similar to take!(IOBuffer) but returns a substring instead of a mutable byte array"
function Base.take!(io::IOShared)
    ss = SubString(io.readofs,usize(io), shared(io,io.writeofs))
    io.readofs = io.writeofs = 0
    return ss
end

@noinline function Base.close(io::IOShared)
    if io.io !== nothing
        flush(io)
        close(io.io)
        io.io = nothing
    end
    io.ioread = false
    io.writeofs = 0
    io.readofs = 0
    io.mark = typemax(UInt32)
    io.compressOnPut = false
    io.preserveOffsets = false
    io.limit = 0 # enforces reallocation on any attempt to use again, used also as flag for isopen
    io.buffer = EMPTYSTRING
    nothing
end

Base.isopen(io::IOShared) = io.limit !=0

"equivalent to read everything"
function Base.empty!(io::IOShared)
    flush(io)
    io.readofs = io.writeofs
    if io.ioread
        close(io.io)
        io.io = nothing
        io.ioread = false
    end
    return io
end


text:
#unum-helpers.jl
#helper functions for the unum constructor, and functions that will frequently
#be used with the constructor.
###########################################################
#Utility functions

################################################################################
# frac_trim - actively trims the fraction to size, setting a ubit as appropriate.

function trim(frac::UInt64, fsize::UInt16)
  frac & mask_top(fsize)
end

function trim!{FSS}(frac::ArrayNum{FSS}, fsize::UInt16)
  middle_cell = div(fsize, 0x0040) + 1
  frac.a[middle_cell] &= mask_top(fsize % 0x0040)
  for idx = (middle_cell + 1):__cell_length(FSS)
    frac.a[idx] = 0
  end
end
#generates the frac_trim!() functions
doc"""
  `Unums.frac_trim!(x::Unum, fsize::UInt16)` strictly trims a fraction to a
  certain size.  This doesn't set the ubit flag if digits are cut off, for that,
  use `Unums.trim_and_set_ubit!`
"""
@fracproc trim fsize

doc"""
  `Unums.exact_trim!(x::Unum)` detects exact unums and trims them to the shortest
  fsize.
"""
@universal function exact_trim!(x::Unum)
  is_ulp(x) && return x  #kick out if we're an ulp.
  tz = ctz(x.fraction)
  if FSS < 7
    x.fsize = (tz == 64) ? 0 : (0x003f - tz)
  else
    x.fsize = (tz == 2 ^ FSS) ? 0 : (2 ^ FSS - 1 - tz)
  end
  return x
end

################################################################################
# trim_and_set_ubit - sets the ubit of the fraction if there is trailing data past the
# desired fsize.

needs_ubit(frac::UInt64, fsize::UInt16) = (frac & mask_bot(fsize)) != 0

function needs_ubit{FSS}(frac::ArrayNum{FSS}, fsize::UInt16)
  accum::UInt64 = zero(UInt64)
  middle_cell = div(fsize, 0x0040) + 1
  accum = frac.a[middle_cell] & mask_bot(fsize % 0x0040)
  for idx = (middle_cell + 1):__cell_length(FSS)
    accum |= frac.a[idx]
  end
  return (accum != 0)
end

doc"""
  `Unums.trim_and_set_ubit!(x::Unum, fsize::UInt16)` is a utility to adjust the
  fsize of a unum.  It checks to see if setting the fsize would need require the
  ubit flag to be thrown, throws it, and then resizes it.
"""
@universal function trim_and_set_ubit!(x::Unum, fsize::UInt16)
  x.flags |= (needs_ubit(x.fraction, fsize) * UNUM_UBIT_MASK)
  frac_trim!(x, fsize)
  x.fsize = min(fsize, max_fsize(FSS))
  return x
end

trim_and_set_ubit!{ESS,FSS}(x::UnumSmall{ESS,FSS}) = trim_and_set_ubit!(x, max_fsize(FSS))
trim_and_set_ubit!{ESS}(x::UnumSmall{ESS,6}) = nothing
trim_and_set_ubit!{ESS,FSS}(x::UnumLarge{ESS,FSS}) = nothing

doc"""
  `Unums.rsh_and_set_ubit!(x::Unum, shift::UInt16)` is a utility that shifts the
  unum fraction to the right and sets the ubit if the fraction gets shifted past
  the right side of the unum.
"""
@universal function rsh_and_set_ubit!(x::Unum, shift::UInt16, needs_guard::Bool = false)
  mfsize::UInt16 = max_fsize(FSS)
  x.flags |= (needs_ubit(x.fraction, mfsize - shift + (needs_guard * o16)) * UNUM_UBIT_MASK)
  x.fsize = min(x.fsize + shift, mfsize)
  frac_rsh!(x, shift)
  return x
end

################################################################################
# EXPONENT ENCODING AND DECODING
doc"""
`Unums.encode_exp(::Int64)` returns a duple `(esize, exponent)` which returns the
most compact Unum representation for an unbiased exponent.  This does not
account for subnormal representations.
"""
function encode_exp(unbiasedexp::Int64)
  #make sure our unbiased exponent is a signed integer
  unbiasedexp = Int64(unbiasedexp)
  esize = to16(64 - leading_zeros(UInt64(abs(unbiasedexp - 1))))
  (esize, UInt64(unbiasedexp + 1 << esize - 1))
end

doc"""
`Unums.decode_exp(::UInt16, ::UInt64)` takes a duple `(esize, exponent)` and
returns the unbiased exponent which is represented by this value.  This does
not account for subnormal representations.
"""
#the inverse operation is finding the unbiased exponent of an Unum.
decode_exp(esize::UInt16, exponent::UInt64) = Int64(exponent) - (1 << esize) + 1

doc"""
`Unums.max_esize(::Int64)` retrieves the maximum possible esize value based on
the ESS value.
"""
max_esize(ESS::Int64) = to16((1 << ESS) - 1)

doc"""
`Unums.max_biased_exponent(::Int64)` retrieves the maximum possible biased exponent
for a given ESS.  Passing a UInt16 instead to `max_biased_exponent(::UInt16)`
signals that your are passing an esize value, instead of a ESS value.  Note that
there is no correspoding min_biased_exponent, because that is always 0.
"""
max_biased_exponent(ESS::Int64) = (o64 << (1 << ESS)) - o64
max_biased_exponent(esize::UInt16) = (o64 << (esize + 1)) - o64

doc"""
`Unums.max_exponent(::Int64)` retrieves the maximum possible unbiased exponent for a given ESS.
"""
max_exponent(ESS::Int64) = Int64(1 << (1 << ESS - 1))

doc"""
`Unums.min_exponent(::Int64)` retrieves the minimum possible unbiased exponent
for a given ESS, not counting subnormal representations.
`min_exponent(ESS::Int64, FSS::Int64)` accounts for this.
"""
min_exponent(ESS::Int64) = Int64(-(1 << (1 << ESS - 1)) + 2)
#and then a minimum exponent that takes into account subnormality.
min_exponent(ESS::Int64, FSS::Int64) = min_exponent(ESS) - (max_fsize(FSS) + 1)

doc"""
`Unums.zero_ulp_index(esize, fsize)` generates an index, based on esize and fsize
that corresponds to how many places to the left of 2 the corresponding "zero ulp"
will be.
"""
#calculates how many decimal places to the right of two a given zero ulp represents.
zero_ulp_index(esize::UInt16, fsize::UInt16) = (o16 << esize) - o16 + fsize

doc"""
`Unums.zero_ulp_params(::Type{Unum{ESS,FSS}}, index)` finds the best parameters
for the zero ulp, given the index of the desired ulp and the constraints set by
ESS and FSS.
"""
@universal function zero_ulp_params(T::Type{Unum}, index::UInt16)
  #first check if we need max_esize.
  mesize = max_esize(ESS)
  mfsize = max_fsize(FSS)

  #are we exceeding the max_esize limit?
  max_eindex = o16 << max_esize(ESS) - o16
  (index > max_eindex) && return (mesize, min(mfsize, index - max_eindex))

  #next, take the logarithm of the index to figure out how far we go.
  esize::UInt16 = 0x000F - clz(index + o16)
  #the remainder becomes the fraction
  fsize::UInt16 = index - o16 << esize + o16
  #unless we have to shortchange it.
  (esize, min(fsize, mfsize))
end

ext: jl

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
# Interactive plot arrays: random walks as examples

using PlotlyJS
using Distributions

# Simulate some random walks
labs = ["J","U","L","I","A"];
N = 10000;
x = 1:N;
y = cumsum(rand(Normal(0,1), N, length(labs)),dims=1);
randomwalks = [scatter(;x=x,y=y[:,i],mode="lines",name=labs[i]) for i in 1:length(labs)];

layout = Layout(
  title = "Intracellular population dynamics",
  xaxis_title="time (d)",
  yaxis_title="Change in # molecules",
  font_family = "Arial",
  font_size = 14
)

# Plot all walks together
plot(randomwalks,layout)

# Plot walks in an array
layout["title"]=""
plots = [plot(randomwalks[i],layout) for i in 1:length(labs)];
[plots[1] plots[2]]
[plots[1] plots[2]  plots[3]
plots[4] plots[5] plot()]


text:
using Random
using Distributions
using MultiAgentPathFinding
using MultiAgentAllocationTransit
using BenchmarkTools
using Statistics
using JSON
using Logging
global_logger(SimpleLogger(stderr, Logging.Error))

rng = MersenneTwister(6789)

# Script arguments relating to transit files etc.
const city_params_file = ARGS[1]
const stop_coords_file = ARGS[2]
const trips_file = ARGS[3]
const drone_params_file = ARGS[4]
const bb_params_file = ARGS[5]
const out_file_pref = ARGS[6]
const N_DEPOTS = parse(Int64, ARGS[7])
const N_AGENTS = parse(Int64, ARGS[8])
const N_TRIALS = parse(Int64, ARGS[9])

# MAPF-TN params
const TRANSIT_CAP_RANGE = (3, 5)
const ECBS_WEIGHT = 1.1

function main()

    # Load the various parameters
    city_params = parse_city_params(city_params_file)
    bb_params = parse_city_params(bb_params_file)
    drone_params = parse_drone_params(drone_params_file)

    lat_dist = Uniform(bb_params.lat_start, bb_params.lat_end)
    lon_dist = Uniform(bb_params.lon_start, bb_params.lon_end)

    # Transit Graph Preprocessing
    tg = load_transit_graph_latlong(stop_coords_file, trips_file, TRANSIT_CAP_RANGE, rng)
    tg, stop_idx_to_trips = transit_graph_preprocessing(tg, MultiAgentAllocationTransit.distance_lat_lon_euclidean, drone_params)

    mapf_results = Dict("cap_range"=>TRANSIT_CAP_RANGE, "weight"=>ECBS_WEIGHT,
                        "depots"=>N_DEPOTS, "agents"=>N_AGENTS)

    search_times = Float64[]
    conflicts = Int64[]
    max_transit_options = Float64[]
    max_path_dists = Float64[]
    makespans = Float64[]

    # Always ignore the first trial
    TRIAL = 1
    while TRIAL <= N_TRIALS+1

        # Inner loop - number of trials
        N_SITES = 2 * N_AGENTS

        println("Trial $(TRIAL)")

        # Generate depots and sites
        depots = [LatLonCoords((lat = rand(rng, lat_dist), lon = rand(rng, lon_dist))) for i = 1:N_DEPOTS]
        sites = [LatLonCoords((lat = rand(rng, lat_dist), lon = rand(rng, lon_dist))) for i = 1:N_SITES]
        depot_sites = vcat(depots, sites)

        # Off transit preprocessing
        otg = OffTransitGraph(depots = depots, sites = sites)

        aug_trips_fws_dists = augmented_trip_meta_graph_fws_dists(tg, MultiAgentAllocationTransit.distance_lat_lon_euclidean,
                                                              length(depots), length(sites),
                                                              vcat(depots, sites),
                                                              drone_params)

        state_graph, depot_sites_to_vtx, trip_to_vtx_range = setup_state_graph(tg, otg)

        # Set up the shadow env to do task allocation with
        env = MAPFTransitEnv(off_transit_graph = otg, transit_graph = tg, state_graph = state_graph,
                             agent_states = AgentState[], depot_sites_to_vtx = depot_sites_to_vtx, trip_to_vtx_range = trip_to_vtx_range,
                             stop_idx_to_trips = stop_idx_to_trips, aug_trips_fws_dists = aug_trips_fws_dists,
                             drone_params = drone_params, dist_fn = MultiAgentAllocationTransit.distance_lat_lon_euclidean,
                             curr_site_points = [], threshold_global_conflicts = 4)

        # run the task allocation, obtain the agent tasks and true number of agents
        cost_fn(i, j) = MultiAgentAllocationTransit.distance_lat_lon_euclidean(depot_sites[i], depot_sites[j])

        agent_tours = task_allocation(N_DEPOTS, N_SITES, N_AGENTS,
                                      depot_sites, cost_fn)

        println("Task Allocation Done!")

        agent_tasks = get_agent_task_set(agent_tours, N_DEPOTS, N_SITES)
        true_n_agents = length(agent_tasks)

        # Now set n_agent dependent env params
        env.agent_states = [AgentState(task=agt_task) for agt_task in agent_tasks]
        env.curr_site_points = zeros(Int64, true_n_agents)
        env.plan_ref_times = zeros(true_n_agents)

        # Setup ECBS Solver; Run search
        initial_states = Vector{MAPFTransitVertexState}(undef, true_n_agents)
        for i = 1:true_n_agents
            initial_states[i] = env.state_graph.vertices[env.depot_sites_to_vtx[string("d-", agent_tasks[i].origin)]]
        end

        solver = ECBSSolver{MAPFTransitVertexState,MAPFTransitAction,Float64,Makespan,MAPFTransitConflict,MAPFTransitConstraints,MAPFTransitEnv}(env = env, weight = ECBS_WEIGHT)

        # Try catch for conflicts
        t = 0.0
        solution = PlanResult{MAPFTransitVertexState,MAPFTransitAction,Float64}[]
        try
            t = @elapsed solution = search!(solver, initial_states)
        catch e
            if isa(e, DomainError)
                println("More than $(env.threshold_global_conflicts) conflicts; exiting!")
            end
            continue
        end

        println("$t seconds; $(solver.num_global_conflicts) conflicts")

        # Compute env diagnostics
        set_solution_diagnostics!(env, solution)

        n_valid_firstpaths = length(env.valid_path_dists)

        if n_valid_firstpaths > 0.67 * true_n_agents

            msp = maximum([s.cost for s in solution])

            push!(search_times, t)
            push!(conflicts, solver.num_global_conflicts)
            push!(max_transit_options, maximum(env.valid_transit_options))
            push!(max_path_dists, maximum(env.valid_path_dists))
            push!(makespans, msp)

            TRIAL = TRIAL + 1
        end
    end

    mapf_results["results"] =   Dict("times" => search_times[2:end],
                                    "conflicts" => conflicts[2:end],
                                    "max_transit_options" => max_transit_options[2:end],
                                    "max_path_dists" => max_path_dists[2:end],
                                    "makespans"=>makespans)

    out_file_name = string(out_file_pref,"_",N_DEPOTS,"deps_",N_AGENTS,"_agts.json")
    open(out_file_name, "w") do f
        JSON.print(f, mapf_results, 2)
    end
end

main()


text:
T = 1.5
dt = 1/1000

tt = 0.:dt:T

x0 = ℝ{2}(2.0, 1.0)

L = @SMatrix [1. 0.]
Σ = @SMatrix [0.1]
v = ℝ{1}(2.5)

γ = 0.7

ϵ = 0.00001


text:
export nnlasso_kernel, lasso_kernel, nnlasso
using SparseArrays
using LinearAlgebra

function nnlasso(A,y,tau; x :: Vector{Float64} = tau*rand(size(A,2))/(size(A,2)*2), l = ones(size(A,2)+1)/(2*size(A,2)))
  K = Array(A'*A);
  b = Array(A'*y);
  return nnlasso_kernel(K,b,tau)
end

function nnlasso_kernel(K,b,tau; x :: Vector{Float64} = tau*rand(length(b))/(length(b)*2), l = ones(length(b)+1)/(2*length(b)))
  K = Array(K)
  b = Array(b)
  n = size(K,2)
  if (n == 0)
    return zeros(0)
  end

  function f_g_h(x)
    0.5*dot(x,K*x) - dot(x,b) ,K*x -b,K
  end
  p = OptimizationProblem(n,n+1,f_g_h,[-I; ones(1,n)], [zeros(n);tau])
  x,l = primalDualSolve(p :: OptimizationProblem, x, l, 10.0, 1e-12 , 1e-12 ,1000, 0.1 , 0.4)
  x[x.<1e-9] .= 0.0
  return x
end

function lasso_kernel(K,b,tau; x :: Vector{Float64} = tau*rand(length(b))/(length(b)*2), l = ones(2*length(b)+1)/(2*length(b)))
  K = Array(K)
  b = Array(b)
  n = size(K,2)
  if ( n == 0)
    return zeros(0)
  end
  function f_g_h(x)
    x = x[1:n]
    f,g,h = 0.5*dot(x,K*x) - dot(x,b) ,K*x -b,K
    f,[g; zeros(n)], [h zeros(n,n); zeros(n,n) zeros(n,n)]
  end
  consrt1 = [eye(n) -eye(n)]
  consrt2 = [-eye(n) -eye(n)]
  p = OptimizationProblem(2*n,2*n + 1,f_g_h,[consrt1;consrt2; [zeros(1,n) ones(1,n)]], [zeros(2*n);tau])
  x,l = primalDualSolve(p :: OptimizationProblem, [x; 1.1*abs(x)], l, 10.0, 1e-10 , 1e-10 ,1000, 0.1 , 0.4)
  x = x[1:n]
  x[abs(x).<1e-9] = 0.0
  return x
end


text:
# This file calcules the ellipticintegral
#
#    ∫_0^z dz/(sqrt(1-z^2)*sqrt(a^2-z^2))
#
# by writing the integrand as a Cauchy transform
# cauchyintegral is an integral of the Cauchy transform

using ApproxFun, SingularIntegralEquations

function ellipticintegral(a)
    xp=Fun(identity,[a,1.])
    xm=Fun(identity,[-a,-1.])
    up=2.0im./(sqrt(1+xp).*sqrt(a+xp).*sqrt(xp-1.).*sqrt(a-xp))
    um=2.0im./(sqrt(-1-xm).*sqrt(a+xm).*sqrt(1.-xm).*sqrt(a-xm))
    g(z)=cauchyintegral(up,z)+cauchyintegral(um,z)
    g
end


text:
# Autogenerated wrapper script for FFTW_jll for powerpc64le-linux-gnu-libgfortran5
export libfftw3, libfftw3_mpi, libfftw3f, libfftw3f_mpi

using MPItrampoline_jll
using MicrosoftMPI_jll
JLLWrappers.@generate_wrapper_header("FFTW")
JLLWrappers.@declare_library_product(libfftw3, "libfftw3.so.3")
JLLWrappers.@declare_library_product(libfftw3_mpi, "libfftw3_mpi.so.3")
JLLWrappers.@declare_library_product(libfftw3f, "libfftw3f.so.3")
JLLWrappers.@declare_library_product(libfftw3f_mpi, "libfftw3f_mpi.so.3")
function __init__()
    JLLWrappers.@generate_init_header(MPItrampoline_jll, MicrosoftMPI_jll)
    JLLWrappers.@init_library_product(
        libfftw3,
        "lib/libfftw3.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libfftw3_mpi,
        "lib/libfftw3_mpi.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libfftw3f,
        "lib/libfftw3f.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libfftw3f_mpi,
        "lib/libfftw3f_mpi.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
abstract type AbstractProjector end

struct L2Projector{CV<:CellValues} <: AbstractProjector
    fe_values::CV
    M_cholesky #::SuiteSparse.CHOLMOD.Factor{Float64}
    dh::MixedDofHandler
    set::Vector{Int}
    node2dof_map::Dict{Int64, Array{Int64,N} where N}
end

function L2Projector(fe_values::Ferrite.Values, interp::Interpolation,
    grid::Ferrite.AbstractGrid, set=1:getncells(grid), fe_values_mass::Ferrite.Values=fe_values)

    dim, T, shape = typeof(fe_values).parameters

    # Create an internal scalar valued field. This is enough since the projection is done on a component basis, hence a scalar field.
    dh = MixedDofHandler(grid)
    field = Field(:_, interp, 1)
    fh = FieldHandler([field], Set(set))
    push!(dh, fh)
    _, vertex_dict, _, _ = __close!(dh)

    M = _assemble_L2_matrix(fe_values_mass, set, dh)  # the "mass" matrix
    M_cholesky = cholesky(M)  # TODO maybe have a lazy eval instead of precomputing? / JB
    return L2Projector(fe_values, M_cholesky, dh, collect(set), vertex_dict[1])
end

function L2Projector(qr::QuadratureRule, func_ip::Interpolation,
    grid::Ferrite.AbstractGrid, set=1:getncells(grid), qr_mass::QuadratureRule=qr,
    geom_ip::Interpolation = default_interpolation(typeof(grid.cells[first(set)])))

    _check_same_celltype(grid, collect(set)) # TODO this does the right thing, but gives the wrong error message if it fails

    fe_values = CellScalarValues(qr, func_ip, geom_ip)
    fe_values_mass = CellScalarValues(qr_mass, func_ip, geom_ip)

    # Create an internal scalar valued field. This is enough since the projection is done on a component basis, hence a scalar field.
    dh = MixedDofHandler(grid)
    field = Field(:_, func_ip, 1) # we need to create the field, but the interpolation is not used here
    fh = FieldHandler([field], Set(set))
    push!(dh, fh)
    _, vertex_dict, _, _ = __close!(dh)

    M = _assemble_L2_matrix(fe_values_mass, set, dh)  # the "mass" matrix
    M_cholesky = cholesky(M)  # TODO maybe have a lazy eval instead of precomputing? / JB
    return L2Projector(fe_values, M_cholesky, dh, collect(set), vertex_dict[1])
end

function _assemble_L2_matrix(fe_values, set, dh)

    n = Ferrite.getn_scalarbasefunctions(fe_values)
    M = create_symmetric_sparsity_pattern(dh)
    assembler = start_assemble(M)

    Me = zeros(n, n)
    cell_dofs = zeros(Int, n)

    function symmetrize_to_lower!(K)
       for i in 1:size(K, 1)
           for j in i+1:size(K, 1)
               K[j, i] = K[i, j]
           end
       end
    end

    ## Assemble contributions from each cell
    for cellnum in set
        celldofs!(cell_dofs, dh, cellnum)

        fill!(Me, 0)
        Xe = getcoordinates(dh.grid, cellnum)
        reinit!(fe_values, Xe)

        ## ∭( v ⋅ u )dΩ
        for q_point = 1:getnquadpoints(fe_values)
            dΩ = getdetJdV(fe_values, q_point)
            for j = 1:n
                v = shape_value(fe_values, q_point, j)
                for i = 1:j
                    u = shape_value(fe_values, q_point, i)
                    Me[i, j] += v ⋅ u * dΩ
                end
            end
        end
        symmetrize_to_lower!(Me)
        assemble!(assembler, cell_dofs, Me)
    end
    return M
end

"""
    project(Vector{Vector{<:Tensor}}}, L2Projector, project_to_nodes=true)

Makes a L2 projection of tensor values to the nodes of the grid. This is commonly used for turning values, computed at integration points, to nodal values, so they can be visualized easier. It is also useful for error estimation and recovering more accurate estimations of secondary unknowns.

If the parameter `project_to_nodes` is true, then the projection returns the values in the order of the mesh nodes. If false, it returns the values corresponding to the degrees of freedom for a scalar field over the domain, which is useful if one wants to interpolate the projected values.
"""
function project(
    vars::Array{Array{T,1},1},
    proj::L2Projector;
    project_to_nodes=true) where {T<:AbstractTensor}

    M = length(vars[1][1].data)
    projected_vals = _project(vars, proj, M)
    if project_to_nodes
        # NOTE we may have more projected values than verticies in the mesh => not all values are returned
        nnodes = getnnodes(proj.dh.grid)
        reordered_vals = fill(NaN, nnodes, size(projected_vals, 2))
        for node = 1:nnodes
            if haskey(proj.node2dof_map, node)
                reordered_vals[node, :] = projected_vals[proj.node2dof_map[node], :]
            end
        end
        return [T(Tuple(reordered_vals[i,:])) for i=1:size(reordered_vals, 1)]
    else
        # convert back to the original tensor type
        return [T(Tuple(reordered_vals[i,:])) for i=1:size(reordered_vals, 1)]
    end
end

function _project(vars, proj::L2Projector, M::Integer)
    # Assemble the multi-column rhs, f = ∭( v ⋅ x̂ )dΩ
    # The number of columns corresponds to the length of the data-tuple in the tensor x̂.

    f = zeros(ndofs(proj.dh), M)
    fe_values = proj.fe_values
    n = getnbasefunctions(proj.fe_values)
    fe = zeros(n, M)

    cell_dofs = zeros(Int, n)
    nqp = getnquadpoints(fe_values)

    ## Assemble contributions from each cell
    for cellnum in proj.set
        celldofs!(cell_dofs, proj.dh, cellnum)
        fill!(fe, 0)
        Xe = getcoordinates(proj.dh.grid, cellnum)
        cell_vars = vars[cellnum]
        reinit!(fe_values, Xe)

        for q_point = 1:nqp
            dΩ = getdetJdV(fe_values, q_point)
            qp_vars = cell_vars[q_point]
            for i = 1:n
                v = shape_value(fe_values, q_point, i)
                fe[i, :] += v * [qp_vars.data[i] for i=1:M] * dΩ
            end
        end

        # Assemble cell contribution
        for (num, dof) in enumerate(cell_dofs)
            f[dof, :] += fe[num, :]
        end
    end

    # solve for the projected nodal values
    return proj.M_cholesky\f

end


text:
include(joinpath(ROOT,"..","src_PowSysMod", "PowSysMod_body.jl"))

function add_to_dict!(dict, key, val)
  if !haskey(dict, key)
    dict[key] = 0
  end
  dict[key] += val
end

"""
pb = read_QCQP(filepath::String)

Read a QCQP file (e.g. WB2.dat) and builds a corresponding Problem type object.
"""
function read_QCQP(filepath::String)
  pb = Problem()

  data = readdlm(filepath)
  line_ind = 1
  while data[line_ind, 1] == "VAR_TYPE"
    add_variable!(pb, Variable(data[line_ind, 3], Complex))
    line_ind += 1
  end

  obj = Polynomial()
  while data[line_ind, 2] == "OBJ"
    obj = obj + build_monomial_from_line(data[line_ind, :], pb)
    line_ind += 1
  end
  set_objective!(pb, obj)

  ## Add constraints
  while line_ind < size(data, 1)
    line = data[line_ind, :]
    current_cstrName = String(line[2])
    current_p = Polynomial()
    current_lb = 0
    current_ub = 0
    while string(line[2]) == current_cstrName && line_ind <= size(data, 1)
      if line[1] == "UB"
        current_ub = line[5] + im*line[6]
      elseif line[1] == "LB"
        current_lb = line[5] + im*line[6]
      else
        current_p = current_p + build_monomial_from_line(line, pb)
      end
      line_ind += 1
      if line_ind <= size(data, 1)
        line = data[line_ind, :]
      end
    end
    add_constraint!(pb, current_cstrName, current_lb << current_p << current_ub)
  end
  pb
end

function build_monomial_from_line(line::Vector{Any}, pb::Problem)
  var1 = read_var(String(line[3]), pb)
  var2 = read_var(String(line[4]), pb)
  (line[5] + im * line[6]) * conj(var1) * var2
end

"""
var = read_var(varName, pb)

Return the Problem variable corresponding to varName if it exists in pb, or else
warn, add it to pb and return the variable.
"""
function read_var(varName::String, pb::Problem)
  var = Variable()
  if varName == "NONE"
    var = Variable("1", Complex)
  elseif !has_variable(pb, varName)
    var = Variable(varName, Complex)
    warn("Variable ", varName, " is not known in the current Problem. Adding ", var, ".")
    add_variable!(pb, var)
  else
    var = Variable(varName, pb.variables[varName])
  end
  var
end

extractWord(l, x) = String(matchall(r"\S+", l)[x])

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
#Additional constants for repeated structs
const TGOV1DU = SteamTurbineGov1

"""
Parse .dyr file into a dictionary indexed by bus number.
Each bus number key has a dictionary indexed by component type and id.

Comments in .dyr files are not supported (beginning of lines with //).

"""
function _parse_dyr_file(file::AbstractString)
    dyr_text = read(file, String)
    start = 1
    parsed_values = Dict{Int, Dict}()
    models = Vector{String}()
    while start < length(dyr_text)
        val = findnext('/', dyr_text, start)
        if isnothing(val)
            break
        end
        text = strip(dyr_text[start:(val - 1)])
        if !isempty(text)
            line = replace(text, "\'" => "")
            line = replace(line, "," => " ")
            val_array = strip.(split(line))
            bus = parse(Int, val_array[1])
            model = string(val_array[2])
            id = string(val_array[3])
            component_dict = get!(parsed_values, bus, Dict{Tuple{String, String}, Array}())
            component_dict[(model, id)] = parse.(Float64, val_array[4:end])
        end
        start = val + 1
    end
    return parsed_values
end

function _parse_input_types(_v, val)
    #If the parameter is a Float64, then use the value directly as the argument.
    #Typically uses for the resistance that is not available in .dyr files.
    if isa(_v, Float64)
        if isnan(_v)
            error("nan for $(_v)")
        end
        return _v
        #If the parameter is an Int, then use the integer as the key of the value in the dictionary.
    elseif isa(_v, Int)
        return val[_v]
    elseif _v == "NaN"
        return NaN
        #If the parameter is a tuple (as a string), then construct the tuple directly.
    elseif isa(_v, String)
        m = match(r"^\((\d+)\s*,\s*(\d+)\)$", _v)
        if m !== nothing
            _tuple_ix = parse.(Int, m.captures)
            return Tuple(val[_tuple_ix])
        else
            error("String $(_v) not recognized for parsing")
        end
    else
        error("invalid input value $val")
    end
    return
end

"""
Populate arguments in a vector for each dynamic component (except Shafts).
Returns a vector with the parameter values of the argument of each component.

"""
function _populate_args(param_map::Vector, val)
    struct_args = Vector{Any}(undef, length(param_map))
    _populate_args!(struct_args, param_map, val, "")
    return struct_args
end

function _populate_args!(struct_args, param_map::Vector, val, ::String)
    for (ix, _v) in enumerate(param_map)
        parameter_value = _parse_input_types(_v, val)
        if !all(isnan.(parameter_value))
            struct_args[ix] = parameter_value
        end
    end
    return
end

function _populate_args!(struct_args, param_dic::Dict, val, id::String)
    param_map = param_dic[id]
    _populate_args!(struct_args, param_map, val, id)
end

"""
Create a SingleMass shaft struct directly using the parameter mapping.

"""
function _make_shaft(param_map, val)
    constructor_shaft =
        (args...) -> InteractiveUtils.getfield(PowerSystems, :SingleMass)(args...)
    shaft_args = _populate_args(param_map, val)
    return constructor_shaft(shaft_args...)
end

function _make_source(g::StaticInjection, r::Float64, x::Float64)
    return Source(
        name = get_name(g),
        available = true,
        bus = get_bus(g),
        active_power = get_active_power(g),
        reactive_power = get_reactive_power(g),
        R_th = r,
        X_th = x,
    )
end

function _assign_missing_components!(
    bus_dict,
    component_table,
    ::Type{T},
) where {T <: DynamicGenerator}
    for va in values(bus_dict)
        if ismissing(va[component_table["AVR"]])
            va[component_table["AVR"]] = AVRFixed(1.0)
        end
        if ismissing(va[component_table["TurbineGov"]])
            va[component_table["TurbineGov"]] = TGFixed(1.0)
        end
        if ismissing(va[component_table["PSS"]])
            va[component_table["PSS"]] = PSSFixed(0.0)
        end
    end
    return
end

function _assign_missing_components!(
    va,
    component_table,
    ::Type{T},
) where {T <: DynamicInverter}
    if ismissing(va[component_table["Converter"]])
        va[component_table["Converter"]] = AverageConverter(750.0, 2.75)
    end
    if ismissing(va[component_table["Filter"]])
        va[component_table["Filter"]] = RLFilter(0.0, 0.0)
    end
    if ismissing(va[component_table["FrequencyEstimator"]])
        va[component_table["FrequencyEstimator"]] = FixedFrequency()
    end
    if ismissing(va[component_table["DCSource"]])
        va[component_table["DCSource"]] = FixedDCSource(750.0)
    end
    return
end

"""
Parse a .dyr file directly from its name by constructing its dictionary of dictionaries.

"""
function _parse_dyr_components(dyr_file::AbstractString)
    ext = splitext(dyr_file)[2]
    if lowercase(ext) in [".dyr"]
        data = _parse_dyr_file(dyr_file)
        return _parse_dyr_components(data)
    else
        throw(DataFormatError("$dyr_file is not a .dyr file type"))
    end
end

"""
Parse dictionary of dictionaries of data (from `_parse_dyr_file`) into a dictionary of struct components.
The function receives the parsed dictionary and constructs a dictionary indexed by bus, that contains a
dictionary with each dynamic generator and dynamic inverter components (indexed via its id).

For Generators, each dictionary indexed by id contains a vector with 5 of its components:
* Machine
* Shaft
* AVR
* TurbineGov
* PSS

For Inverters, each dictionary indexed by id contains a vector with 7 of its components:
* Converter
* ActivePowerControl
* ReactivePowerControl
* InnerControl
* DCSource
* FrequencyEstimator
* Filter

"""
function _parse_dyr_components(data::Dict)
    yaml_mapping = open(PSSE_DYR_MAPPING_FILE) do io
        aux = YAML.load(io)
        return aux
    end
    #param_map contains the mapping between struct params and dyr file.
    param_map = yaml_mapping["parameter_mapping"][1]
    #gen_map contains all the supported structs for generators
    gen_map = yaml_mapping["generator_mapping"][1]
    #inv_map contains al the supported structs for inverters
    inv_map = yaml_mapping["inverter_mapping"][1]
    #dic will contain the dictionary index by bus.
    #Each entry will be a dictionary, with id as keys, that contains the vector of components
    dic = Dict{Int, Any}()
    _parse_dyr_generator_components!(dic, data, gen_map, param_map)
    _parse_dyr_inverter_components!(dic, data, inv_map, param_map)
    return dic
end

"""
Parse dictionary of data (from `_parse_dyr_file`) into a dictionary of struct components.
The function receives the parsed dictionary and constructs a dictionary indexed by bus, that contains a
dictionary with each dynamic generator indexed by its id.

"""
function _parse_dyr_generator_components!(
    dic::Dict,
    data::Dict,
    gen_map::Dict,
    param_map::Dict,
)
    component_table =
        Dict("Machine" => 1, "Shaft" => 2, "AVR" => 3, "TurbineGov" => 4, "PSS" => 5)
    for (bus_num, bus_data) in data #bus_data is a dictionary with values per component (key)
        bus_dict = Dict{String, Any}()
        for (componentID, componentValues) in bus_data
            #Fill array of 5 components per generator
            temp = get!(bus_dict, componentID[2], fill!(Vector{Any}(undef, 5), missing))
            #Only create if name is in the supported keys in the mapping
            if componentID[1] in keys(gen_map)
                #Get the component dictionary
                components_dict = gen_map[componentID[1]]
                for (gen_field, struct_as_str) in components_dict
                    #Get the parameter indices w/r to componentValues
                    params_ix = param_map[struct_as_str]
                    #Construct Shaft if it appears
                    if occursin("Shaft", struct_as_str)
                        temp[2] = _make_shaft(params_ix, componentValues)
                        continue
                    end
                    #Construct Components
                    component_constructor =
                        (args...) ->
                            InteractiveUtils.getfield(PowerSystems, Symbol(struct_as_str))(
                                args...,
                            )
                    struct_args = _populate_args(params_ix, componentValues)
                    if struct_as_str == "BaseMachine"
                        struct_args = [
                            0.0,    #R
                            1e-5,   #X
                            1.0,    #eq_p
                        ]
                    end
                    temp[component_table[gen_field]] = component_constructor(struct_args...)
                end
            else
                @warn "$(componentID[1]) at bus $bus_num, id $(componentID[2]), not supported in PowerSystems.jl. Skipping data."
            end
        end
        #Assign generic components if there were not provided in Dynamic Data
        _assign_missing_components!(bus_dict, component_table, DynamicGenerator)
        #Store dictionary of components in a dictionary indexed by bus
        dic[bus_num] = bus_dict
    end
end

"""
Parse dictionary of data (from `_parse_dyr_file`) into a dictionary of struct components.
The function receives the parsed dictionary and constructs a dictionary indexed by bus, that contains a
dictionary with each dynamic inverter indexed by its id.

"""
function _parse_dyr_inverter_components!(
    dic::Dict,
    data::Dict,
    inv_map::Dict,
    param_map::Dict,
)
    component_table = Dict(
        "Converter" => 1,
        "ActivePowerControl" => 2,
        "ReactivePowerControl" => 3,
        "InnerControl" => 4,
        "DCSource" => 5,
        "FrequencyEstimator" => 6,
        "Filter" => 7,
    )
    for (bus_num, bus_data) in data #bus_data is a dictionary with values per component (key is a tuple of PSSE name and number ID)
        #bus_dict will add a vector of components structs for each generator ID on the bus_num key of the master dictionary 'dic'
        #bus_dict_values will only contain the vector of parameters that will be used to construct the structs
        bus_dict_values = Dict{String, Any}()
        for (componentID, componentValues) in bus_data
            #ComponentID is a tuple:
            #ComponentID[1] is the PSSE name
            #ComponentID[2] is the number ID of the generator/inverter
            #Fill array of 7 components per inverter
            #Only create if name is in the supported keys in the mapping
            if componentID[1] in keys(inv_map)
                #Get the component dictionary
                temp = get!(bus_dict_values, componentID[2], Dict{Any, Vector{Any}}())
                components_dict = inv_map[componentID[1]]
                for (inv_field, struct_as_str) in components_dict
                    param_vec = get!(
                        temp,
                        (inv_field, struct_as_str),
                        _instantiate_param_vector_size(struct_as_str, param_map),
                    )
                    params_ix = param_map[struct_as_str]
                    _populate_args!(param_vec, params_ix, componentValues, componentID[1])
                end
            end
        end
        #Store vector of components in a dictionary indexed by bus
        for (device_id, device_info) in bus_dict_values
            temp_vec = fill!(Vector{Any}(undef, 7), missing)
            for (k, struct_args) in device_info
                component_type, struct_as_str = k
                _convert_argument_types!(struct_as_str, struct_args)
                component_constructor =
                    (args...) ->
                        InteractiveUtils.getfield(PowerSystems, Symbol(struct_as_str))(
                            args...,
                        )
                temp_vec[component_table[component_type]] =
                    component_constructor(struct_args...)
            end
            _assign_missing_components!(temp_vec, component_table, DynamicInverter)
            dic[bus_num][device_id] = temp_vec
        end
    end
end

"""
Construct appropiate vector size for components that collect parameters from
more than 2 PSS/E components

"""
function _instantiate_param_vector_size(str::AbstractString, param_map::Dict)
    if str == "ActiveRenewableControllerAB"
        return fill!(Vector{Any}(undef, 17), missing)
    elseif str == "ReactiveRenewableControllerAB"
        return fill!(Vector{Any}(undef, 25), missing)
    elseif str == "RECurrentControlB"
        return fill!(Vector{Any}(undef, 12), missing)
    elseif str in keys(param_map)
        return fill!(Vector{Any}(undef, length(param_map[str])), missing)
    else
        error("String $(str) not supported in the parser")
    end
end

"""
Convert specific parameters to types that are not Float64 for
specific inverter components

"""
function _convert_argument_types!(str::AbstractString, struct_args::Vector)
    if str == "ActiveRenewableControllerAB"
        struct_args[1:5] .= Int.(struct_args[1:5])
        struct_args[4] = string(struct_args[4])
    elseif str == "ReactiveRenewableControllerAB"
        struct_args[1:8] .= Int.(struct_args[1:8])
        struct_args[4] = string(struct_args[4])
    elseif str == "RECurrentControlB"
        struct_args[1:2] .= Int.(struct_args[1:2])
    elseif str == "RenewableEnergyConverterTypeA"
        #No changes to struct_args
    else
        error("$str not defined for dynamic component arguments")
    end
end

"""
Add to a system already created the dynamic components.
The system should already be parsed from a .raw file.

## Examples:
```julia
dyr_file = "Example.dyr"
add_dyn_injectors!(sys, dyr_file)
```

"""
function add_dyn_injectors!(sys::System, dyr_file::AbstractString)
    bus_dict_gen = _parse_dyr_components(dyr_file)
    add_dyn_injectors!(sys, bus_dict_gen)
end

function add_dyn_injectors!(sys::System, bus_dict_gen::Dict)
    @info "Generators provided in .dyr, without a generator in .raw file will be skipped."
    for g in collect(get_components(ThermalStandard, sys))
        _num = get_number(get_bus(g))
        _name = get_name(g)
        _id = split(_name, "-")[end]
        temp_dict = get(bus_dict_gen, _num, nothing)
        if temp_dict === nothing
            @warn "Generator at bus $(_num), id $(_id), not found in Dynamic Data.\nVoltage Source will be used to model it."
            r, x = get_ext(g)["z_source"]
            if x == 0.0
                @warn "No series reactance found. Setting it to 1e-6"
                x = 1e-6
            end
            s = _make_source(g, r, x)
            remove_component!(typeof(g), sys, _name)
            add_component!(sys, s)
        elseif all(.!ismissing.(temp_dict[_id]))
            if length(temp_dict[_id]) == 5 #Generator has 5 components
                #Obtain Machine from Dictionary
                machine = temp_dict[_id][1]
                #Update R,X from RSORCE and XSORCE from RAW file
                r, x = get_ext(g)["z_source"]
                set_R!(machine, r)
                #Obtain Shaft from dictionary
                shaft = temp_dict[_id][2]
                if typeof(machine) == BaseMachine
                    if x == 0.0
                        @warn "No series reactance found. Setting it to 1e-6"
                        x = 1e-6
                    end
                    set_Xd_p!(machine, x)
                    if get_H(shaft) == 0.0
                        @info "Machine at bus $(_num), id $(_id) has zero inertia. Modeling it as Voltage Source"
                        s = _make_source(g, r, x)
                        remove_component!(typeof(g), sys, _name)
                        add_component!(sys, s)
                        #Don't add DynamicComponent in case of adding Source
                        continue
                    end
                end
                #Add Dynamic Generator
                dyn_gen = DynamicGenerator(get_name(g), 1.0, temp_dict[_id]...)
                add_component!(sys, dyn_gen, g)
            elseif length(temp_dict[_id]) == 7 #Inverter has 7 components (6 if Outer is put together)
                dyn_inv = DynamicInverter(get_name(g), 1.0, temp_dict[_id]...)
                add_component!(sys, dyn_inv, g)
            else
                @assert false
            end
        else
            error("Generator at bus $(_num), id $(_id), name $(_name), not supported")
        end
    end
end


text:
export ExplicitSolverType

"""
# Description
    ExplicitSolverType(;
        solver_method = LSRK54CarpenterKennedy,
    )

This solver type constructs an ODE solver using an explicit
Runge-Kutta method.

# Arguments
- `solver_method` (Function): Function defining the explicit
    Runge-Kutta solver.
    Default: `LSRK54CarpenterKennedy`
"""
struct ExplicitSolverType <: AbstractSolverType
    # Function for an explicit Runge-Kutta method
    solver_method::Function

    function ExplicitSolverType(; solver_method = LSRK54CarpenterKennedy)

        return new(solver_method)
    end
end

"""
    getdtmodel(ode_solver::AbstractSolverType, bl)

A function which returns a model representing the dynamics
with the most restrictive time-stepping requirements.
"""
function getdtmodel(::ExplicitSolverType, bl)
    # For explicit methods, the entire model itself
    # contributes to the total stability of the time-integrator
    return bl
end

"""
# Description
    function solversetup(
        ode_solver::ExplicitSolverType,
        dg,
        Q,
        dt,
        t0,
        diffusion_direction,
    )

Creates an explicit ODE solver.
"""
function solversetup(
    ode_solver::ExplicitSolverType,
    dg,
    Q,
    dt,
    t0,
    diffusion_direction,
)

    solver = ode_solver.solver_method(dg, Q; dt = dt, t0 = t0)

    return solver
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using FirstOrderSolvers: KKTMatrix, AffinePlusLinear
using Convex

Random.seed!(10)

#Test KKTMatrix
A = randn(10,20)
M1 = [I A';A -I]
M2 = KKTMatrix(A)
x = randn(30)
y1 = randn(30)
y2 = randn(30)
mul!(y1, M1, x)
mul!(y2, M2, x)
@test y1 ≈ y2

mul!(y1, transpose(M1), x)
mul!(y2, transpose(M2), x)
@test y1 ≈ y2

#Test AffinePlusLinear

x0 = randn(20)
z0 = randn(10)
q = randn(20)
b = randn(10)

# β = 1
β = 1
S2 = AffinePlusLinear(A, b, q, β)
y2 = Array{Float64,1}(undef, 30)
FirstOrderSolvers.prox!(y2, S2, [x0;z0])
x2 = view(y2,1:20)
z2 = view(y2,21:30)

# #This test is equivalent to y3
# x = Variable(20)
# z = Variable(10)
# p = minimize(1/2*norm(x-x0)^2+1/2*norm(z-z0)^2+dot(q,x), [A*x-z==b])
# solve!(p, FirstOrderSolvers.DR(eps=1e-9))
# x1 = x.value
# z1 = z.value
# y1 = [x1; z1]
# @test y1 ≈ y2

y3 = [I A'; A -I]\[x0-q+A'z0; b]
@test y3 ≈ y2

β = -1
S2 = AffinePlusLinear(A, b, q, β)
y2 = Array{Float64,1}(undef, 30)
FirstOrderSolvers.prox!(y2, S2, [x0;z0])
x2 = view(y2,1:20)
z2 = view(y2,21:30)

# #This test is equivalent to y3
# x = Variable(20)
# z = Variable(10)
# p = minimize(1/2*norm(x-x0)^2+1/2*norm(z-z0)^2+dot(q,x), [A*x-β*z==b])
# solve!(p, FirstOrderSolvers.DR(eps=1e-9))
# x1 = x.value
# z1 = z.value
# y1 = [x1; z1]
# @test y1 ≈ y2

y3 = [I -A'; A I]\[x0-q-A'z0; b]
@test y3 ≈ y2


text:
using DataFrames
using StatPlots
using DataArrays
#Importing in dateset
cpi = readtable("cpi.csv")
#setting as date format
cpidate = Date(cpi[1])
#seperating dataset
cpidata = cpi[2]
#normalizing
cpinormal = cpi[2]/sum(cpi[2])
#reading in higher ed dataset
highered = readtable("highered.csv")
#seperating out date
highereddate = DataArray(highered[1])
#making as Date
highereddate = Date(highereddate)
#normalizing
higherednormal = highered[2]/sum(highered[2])
#plot unnormalized cpi vs highered
plot(cpidate, cpidata, label = "Urban CPI")
plot!(highereddate, highered[2], label = "Higher Education CPI")
title!("Urban CPI vs. Higher Education")
xaxis!("Date")
yaxis!("Price")

#plot normalized Plots
plot(cpidate, cpinormal)
plot!(highereddate, higherednormal)

max_stars_repo_name: raagagrawal/highered

max_issues_repo_name: raagagrawal/highered

max_forks_repo_name: raagagrawal/highered


text:
using Printf

const tcode = (Regex=>Char)[r"A|B|C|Ä|Å|Á|Â|Ç" => '2',
                            r"D|E|F|È|Ê|É" => '3',
                            r"G|H|I|Í" => '4',
                            r"J|K|L" => '5',
                            r"M|N|O|Ó|Ö|Ô|Ñ" => '6',
                            r"P|Q|R|S" => '7',
                            r"T|U|V|Û|Ü" => '8',
                            r"W|X|Y|Z" => '9']

function tpad(str::IOStream)
    tnym = (String=>Array{String,1})[]
    for w in eachline(str)
        w = chomp(w)
        t = uppercase(w)
        for (k,v) in tcode
            t = replace(t, k, v)
        end
        t = replace(t, r"\D", '1')
        tnym[t] = [get(tnym, t, String[]), w]
    end
    return tnym
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
   set_inactive!(state, agent_idx, post_list)

Sets an agent's state to inactive and deletes all his connections to others and posts

# Arguments
- `state`: a tuple of the current graph and agent_list
- `agent_idx`: agent index
- `post_list`: List of all published posts in network

See also: [`Agent`](@ref), [`tick!`](@ref)
"""
function set_inactive!(
    state::Tuple{AbstractGraph, AbstractArray},
    agent_idx::Integer,
    post_list::AbstractArray
)

    graph, agent_list = state
    this_agent = agent_list[agent_idx]

    this_agent.active = false
    agent_edges = [
        e for e in edges(graph) if (src(e) == agent_idx || dst(e) == agent_idx)
    ]

    for e in agent_edges
        rem_edge!(graph, e)
    end

    empty!(this_agent.feed)

    for p in post_list
        if p.source_agent == agent_idx
            p.weight = -1
        end
    end

    return state
end

# suppress output of include()
;


text:
module PEPSAD

# peps and models
export PEPS, randomizePEPS!, inner_network, mpo, localham, checklocalham, Model, prime
# optimizations
export gradient_descent, optimize, generate_inner_network, extract_data

include("peps.jl")
include("models.jl")
include("optimizations.jl")

end


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

# ToDo: Replace with custom UnshapedUvDist:
unshaped(d::UnivariateDistribution) = Distributions.Product(Fill(d, 1))


text:
using HorizonSideRobots

mutable struct Coord
    x::Int
    y::Int
    Coord() = new(0,0)
    Coord(x::Int,y::Int) = new(x,y)
end

function move!(coord::Coord, side::HorizonSide)
    if side==Nord
        coord.y += 1
    elseif side==Sud
        coord.y -= 1
    elseif side==Ost
        coord.x += 1
    else
        coord.x -= 1
    end
end

get_coords(coord::Coord) = (coord.x, coord.y)

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
# using Revise, Test, Plots
using BifurcationKit, LinearAlgebra, Setfield
const BK = BifurcationKit

# display internal information about the branch
function displayBr(contRes)
	println("#"^50)
	for ii in 1:length(contRes.branch)
		println("- $ii --------------")
		println("step = ", contRes[ii][end])
		println("eiv = "); display(contRes.eig[ii].eigenvals)
		println("stab = ", contRes[ii].stable)
		println("n_uns = ", contRes[ii].n_unstable)
	end
end

function testBranch(br)
	# test if stability works
	# test that stability corresponds
	out = true
	for ii in eachindex(br.branch)
		@test br.eig[ii].eigenvals == eigenvals(br, br.eig[ii].step, false)
		# compute number of unstable eigenvalues
		isstable, n_u, n_i = BK.isStable(br.contparams, br.eig[ii].eigenvals)
		# test that the stability matches the one in eig
		@test br[ii].n_unstable == n_u
		@test br[ii].stable  == isstable
		# we test that step = ii-1
		@test br.branch[ii][end] == ii-1
		# test that the field `step` match in the structure
		@test br.branch[ii][end] == br.eig[ii].step
	end
	# test about bifurcation points
	for bp in br.specialpoint
		id = bp.idx
		# test that the states marked as bifurcation points are always after true bifurcation points
		@test abs(br[id].n_unstable - br[id-1].n_unstable) > 0
		# test that the bifurcation point belongs to the interval
		@test bp.interval[1] <= bp.param <= bp.interval[2]
		# test that bp.param = br[id].param
		@test bp.param == br[id].param
	end
end

NL(x) = -x^3
dNL(x) = -3x^2
Ftb(x,p) = -x .+ (p.L * x) .* p.λ .+ NL.(x)

function Jtb(x, p)
	J = copy(p.L .* p.λ)
	for i in eachindex(x)
		J[i,i] += dNL(x[i]) - 1
	end
	return J
end

par = (L = Diagonal([1.0/ii for ii in 1:5 for jj in 1:ii]), λ = .0)
append!(par.L.diag, [1/6. 1/6.5 1/6.75 1/6.875])

# ensemble of bifurcation points
specialpoints = unique(1 ./par.L.diag);
dimBif = [ii for ii in 1:5]; append!(dimBif, [1 1 1 1])

x0 = zeros(size(par.L, 1))

optc = ContinuationPar(pMin = -1., pMax = 10., ds = 0.1, maxSteps = 150, detectBifurcation = 2, saveEigenvectors = false)
br1, = continuation(Ftb, Jtb, x0, par, (@lens _.λ), optc; verbosity = 0)
testBranch(br1)

br2, = continuation(Ftb, Jtb, x0, par, (@lens _.λ), setproperties(optc; detectBifurcation = 3, pMax = 10.3, nInversion = 4, tolBisectionEigenvalue = 1e-7); plot=false, verbosity = 0)
testBranch(br2)
for bp in br2.specialpoint
	@test bp.interval[1] <= bp.param <= bp.interval[2]
end

specialpoint2 = [bp.param for bp in br2.specialpoint]
@test specialpoint2 > specialpoints
@test norm(specialpoints - specialpoint2, Inf) < 3e-3
dimBif2 = [abs(bp.δ[1]) for bp in br2.specialpoint]
@test dimBif2 == dimBif

# case where bisection "fails". Test whether the bifurcation point belongs to the specified interval
br3, = continuation(Ftb, Jtb, x0, par, (@lens _.λ), setproperties(optc; detectBifurcation = 3, pMax = 10.3, nInversion = 8, tolBisectionEigenvalue = 1e-7); verbosity = 0)
testBranch(br3)

# case where bisection "fails". Test whether the bifurcation point belongs to the specified interval
# in this case, we test if coming from above, and having no inversion, still leads to correct result
br4, = continuation(Ftb, Jtb, x0, (@set par.λ = 0.95), (@lens _.λ), setproperties(optc; detectBifurcation = 3, pMax = 1.95, nInversion = 8, ds = 0.7, dsmax = 1.5, maxBisectionSteps = 1); verbosity = 0)
testBranch(br4)
####################################################################################################
# this example is to test failures in Newton annd how it affects the bifurcation points labels
using ForwardDiff
F = (x, p; k = 3) -> (@. p * x -  x^k/k)
Jac_m = (x, p; k = 2) -> diagm(0 => p .- x.^k)

opts = ContinuationPar(dsmax = 0.1, dsmin = 1e-5, ds = 0.001, maxSteps = 130, pMin = -3., pMax = 0.1, saveSolEveryStep = 0, newtonOptions = NewtonPar(tol = 1e-8, verbose = false, maxIter = 4), detectBifurcation=3, nInversion=4)

br4, = continuation(F, Jac_m, zeros(1), -0.1, (@lens _), opts; verbosity = 0, plot=false,
	printSolution = (x,p)->x[1])
testBranch(br4)
####################################################################################################
using ForwardDiff
function Ftb(X, p)
	p1, p2, k = p
	x, y = X
	out = similar(X)
	out[1] = p1 + x - y - x^k/k
	out[2] = p1 + y + x - 2y^k/k
	out
end

Jtb2 = (X, p) -> ForwardDiff.jacobian(z -> Ftb(z,p), X)
par = (p1 = -3., p2=-3., k=3)

opts = ContinuationPar(dsmax = 0.1, ds = 0.001, maxSteps = 135, pMin = -3., pMax = 4.0, newtonOptions = NewtonPar(maxIter = 5), detectBifurcation = 3, nInversion = 6, dsminBisection = 1e-9, maxBisectionSteps = 15, nev = 2)

br, = continuation(Ftb, Jtb2, -2ones(2), par, (@lens _.p1), (@set opts.detectBifurcation = 3);
	plot = false, verbosity = 0,
	printSolution = (x, p) -> x[1],)
	show(br)
testBranch(br)

br, = continuation(Ftb, Jtb2, -2ones(2), par, (@lens _.p1), opts)
testBranch(br)


text:
"""
    function d2f_PgPg_nonzeroslope(args...; kwargs...)

Calculates the second derivative of the total generation cost function for the generators
w.r.t. the active power output (Pg). Since the cost function for each generator is assumed
to be polinomial of the form: Gencost = αPg² + βPg + γ, the derivative is only calculated
for the generators with non zero slopes. Their second derivative is: 2α, for all the other
generators the second derivative is zero. Also, only self elements have a second derivative.
Note: a non zero slope generator will have the parameter α ≠ 0.

# Arguments
- `sys::System`:                            Power system in p.u. (from PowerSystems.jl)
- `res::PowerSimulations.
    OperationsProblemResults`:              Results of the solved OPF for the system
                                             (from PowerSimulations.jl)

# Keywords
- `dual_gen_tol::Float64 = 1e-1`:           Tolerance to identify any binding generators

"""
function d2f_PgPg_nonzeroslope(
    sys::System,
    res::PowerSimulations.OperationsProblemResults;
    dual_gen_tol::Float64 = 1e-1
    )
    # Define elements
    gen_loc = 0
    gen_nonzero = 0
    d2f_PgPgnonzeroslope_array  = Array{Float64}(undef, 0)
    gens_thermal_nonzeroslope_names = Array{String}(undef, 0)
    gens_thermal_nonzeroslope_busnumbers = Array{Int64}(undef, 0)

    # Get all thermal Generators
    gens_thermal = get_components(ThermalStandard, sys)
    all_PGenThermal = get_variables(res)[:P__ThermalStandard] #Optimised PGen

    # Get dual variables from the OPF results
    dualPgub = get_duals(res)[:P_ub__ThermalStandard__RangeConstraint]
    dualPglb = get_duals(res)[:P_lb__ThermalStandard__RangeConstraint]

    # Identify Non zero Slope Generators and calculate their derivative
    for gen_thermal in gens_thermal
        gen_loc = gen_loc +1
        if get_bustype(gen_thermal.bus) ≠ BusTypes.REF  #skips slack
            bindub = !isapprox(dualPgub[1, gen_loc], 0.0; atol = dual_gen_tol) #is binding?
            bindlb = !isapprox(dualPglb[1, gen_loc], 0.0; atol = dual_gen_tol) #is binding?
            if !bindub && !bindlb #skips binding Pg gens Dual = 0
                α, β = get_cost(get_variable(get_operation_cost(gen_thermal)))
                if α ≠ 0 # Only gens with α(Pg^2)
                    #Only self elements have second derivative
                    gen_nonzero = gen_nonzero + 1
                    resize!(d2f_PgPgnonzeroslope_array, gen_nonzero)
                    resize!(gens_thermal_nonzeroslope_names, gen_nonzero)
                    resize!(gens_thermal_nonzeroslope_busnumbers, gen_nonzero)
                    gen_thermal_nonzeroslope_name = gen_thermal.name
                    d2f_PgPgnonzeroslope_array[gen_nonzero] = 2*α
                    gens_thermal_nonzeroslope_names[gen_nonzero] =
                        gen_thermal_nonzeroslope_name
                    gens_thermal_nonzeroslope_busnumbers[gen_nonzero] =
                        gen_thermal.bus.number
                end
            end
        end
    end
    # Matrix of information
    d2f_PgPgnonzeroslope_MatrixInfo = Array{Any}(undef, (gen_nonzero, 3))
    d2f_PgPgnonzeroslope  = Array{Float64}(undef, (gen_nonzero, 1))
    d2f_PgPgnonzeroslope_MatrixInfo[:, 1] = gens_thermal_nonzeroslope_names
    d2f_PgPgnonzeroslope_MatrixInfo[:, 2] = d2f_PgPgnonzeroslope_array
    d2f_PgPgnonzeroslope_MatrixInfo[:, 3] = gens_thermal_nonzeroslope_busnumbers
    # Sort by bus to match with PTDF structure
    d2f_PgPgnonzeroslope_MatrixInfo =
        d2f_PgPgnonzeroslope_MatrixInfo[sortperm(d2f_PgPgnonzeroslope_MatrixInfo[:, 3]), :]

    # Second Derivative Matrix
    d2f_PgPgnonzeroslope[:,1]  = d2f_PgPgnonzeroslope_MatrixInfo[:, 2]
    d2f_PgPgnonzeroslope_matrix = Diagonal(d2f_PgPgnonzeroslope[:, 1])

    return (d2f_PgPgnonzeroslope_matrix, d2f_PgPgnonzeroslope_MatrixInfo)
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD