text:
# calculate crc checksums for data.  this process is closely related
# to the routines in IntModN (particularly GF2Poly), but the
# traditional implementations are optimized to the point where
# implementing them from "first principles" makes little sense
# (although we can use IntModN to check the results here).

# the basic idea is to do a polynomial division along the byte stream.
# it's all explained very nicely at
# http://en.wikipedia.org/wiki/Computation_of_cyclic_redundancy_checks
# and http://www.zlib.net/crc_v3.txt

module CRC

export crc, make_tables, spec, CHECK, NoTables, Single, Multiple,
       CRC_3_ROHC, CRC_4_ITU, CRC_5_EPC, CRC_5_ITU, CRC_5_USB,
       CRC_6_CDMA2000_A, CRC_6_CDMA2000_B, CRC_6_DARC, CRC_6_ITU,
       CRC_7, CRC_7_ROHC, CRC_8, CRC_8_CDMA2000, CRC_8_DARC,
       CRC_8_DVB_S2, CRC_8_EBU, CRC_8_I_CODE, CRC_8_ITU, CRC_8_MAXIM,
       CRC_8_ROHC, CRC_8_WCDMA, CRC_10, CRC_10_CDMA2000, CRC_11,
       CRC_12_3GPP, CRC_12_CDMA2000, CRC_12_DECT, CRC_13_BBC,
       CRC_14_DARC, CRC_15, CRC_15_MPT1327, CRC_16_ARC,
       CRC_16_AUG_CCITT, CRC_16_BUYPASS, CRC_16_CCITT_FALSE,
       CRC_16_CDMA2000, CRC_16_DDS_110, CRC_16_DECT_R, CRC_16_DECT_X,
       CRC_16_DNP, CRC_16_EN_13757, CRC_16_GENIBUS, CRC_16_MAXIM,
       CRC_16_RIELLO, CRC_16_TELEDISK, CRC_16_USB, CRC_16_CRC_A,
       CRC_16_KERMIT, CRC_16_MODBUS, CRC_16_X_25, CRC_16_XMODEM,
       CRC_24, CRC_24_FLEXRAY_A, CRC_24_FLEXRAY_B, CRC_31_PHILIPS,
       CRC_32, CRC_32_BZIP2, CRC_32_C, CRC_32_D, CRC_32_MPEG_2,
       CRC_32_POSIX, CRC_32_Q, CRC_32_JAMCRC, CRC_32_XFER, CRC_40_GSM,
       CRC_64, CRC_64_WE, CRC_64_XZ, CRC_82_DARC, ALL, main

import Base.Cartesian: @nexprs
import Base: ==, isless, print
using ArgParse

typealias U Unsigned

# ---- CRC specifications from http://www.zlib.net/crc_v3.txt

CHECK = b"123456789"   # universal check vector

immutable Spec{P<:U}
    width::Int    # polynomial degree
    poly::P       # generating poly with msb missing
    init::P       # initial remainder
    refin::Bool   # reflect input
    refout::Bool  # reflect output
    xorout::P     # xored with final remainder
    check::P      # checksum for CHECK
    function Spec(width::Int, poly::P, init::P, refin::Bool, refout::Bool, xorout::P, check::P)
        @assert width <= 8 * sizeof(P)
        new(width, poly, init, refin, refout, xorout, check)
    end
end

spec{P<:U}(poly::P, init::P, refin::Bool, refout::Bool, xorout::P, check::P) =
    Spec{P}(8*sizeof(P), poly, init, refin, refout, xorout, check)

spec{P<:U}(width::Int, poly::P, init::P, refin::Bool, refout::Bool, xorout::P, check::P) =
    Spec{P}(width, poly, init, refin, refout, xorout, check)

==(a::Spec, b::Spec) = a.width == b.width && a.poly == b.poly && a.init == b.init && a.refin == b.refin && a.refout == b.refout && a.xorout == b.xorout

isless(a::Spec, b::Spec) = a.width < b.width || (a.width == b.width && (a.poly < b.poly || a.poly == b.poly && (a.init < b.init || (a.init == b.init && !a.refin && (a.refin != b.refin || (a.refin == b.refin && !a.refout && (a.refout != b.refout || a.xorout < b.xorout)))))))

function print{P<:U}(io::IO, s::Spec{P})
    n = 1 + div(s.width - 1, 4)
    h = x -> "0x" * hex(x, n)
    print(io, "width=$(s.width) poly=$(h(s.poly)) init=$(h(s.init)) refin=$(s.refin) refout=$(s.refout) xorout=$(h(s.xorout)) check=$(h(s.check))")
end

# http://reveng.sourceforge.net/crc-catalogue/1-15.htm
CRC_3_ROHC =         spec(3, 0x03, 0x07, true,  true,  0x00, 0x06)
CRC_4_ITU =          spec(4, 0x03, 0x00, true,  true,  0x00, 0x07)
CRC_5_EPC =          spec(5, 0x09, 0x09, false, false, 0x00, 0x00)
CRC_5_ITU =          spec(5, 0x15, 0x00, true,  true,  0x00, 0x07)
CRC_5_USB =          spec(5, 0x05, 0x1f, true,  true,  0x1f, 0x19)
CRC_6_CDMA2000_A =   spec(6, 0x27, 0x3f, false, false, 0x00, 0x0d)
CRC_6_CDMA2000_B =   spec(6, 0x07, 0x3f, false, false, 0x00, 0x3b)
CRC_6_DARC =         spec(6, 0x19, 0x00, true,  true,  0x00, 0x26)
CRC_6_ITU =          spec(6, 0x03, 0x00, true,  true,  0x00, 0x06)
CRC_7 =              spec(7, 0x09, 0x00, false, false, 0x00, 0x75)
CRC_7_ROHC =         spec(7, 0x4f, 0x7f, true,  true,  0x00, 0x53)
CRC_8 =              spec(0x07, 0x00, false, false, 0x00, 0xf4)
CRC_8_CDMA2000 =     spec(0x9b, 0xff, false, false, 0x00, 0xda)
CRC_8_DARC =         spec(0x39, 0x00, true,  true,  0x00, 0x15)
CRC_8_DVB_S2 =       spec(0xd5, 0x00, false, false, 0x00, 0xbc)
CRC_8_EBU =          spec(0x1d, 0xff, true,  true,  0x00, 0x97)
CRC_8_I_CODE =       spec(0x1d, 0xfd, false, false, 0x00, 0x7e)
CRC_8_ITU =          spec(0x07, 0x00, false, false, 0x55, 0xa1)
CRC_8_MAXIM =        spec(0x31, 0x00, true,  true,  0x00, 0xa1)
CRC_8_ROHC =         spec(0x07, 0xff, true,  true,  0x00, 0xd0)
CRC_8_WCDMA =        spec(0x9b, 0x00, true,  true,  0x00, 0x25)
CRC_10 =             spec(10, 0x0233, 0x0000, false, false, 0x0000, 0x0199)
CRC_10_CDMA2000 =    spec(10, 0x03d9, 0x03ff, false, false, 0x0000, 0x0233)
CRC_11 =             spec(11, 0x0385, 0x001a, false, false, 0x0000, 0x05a3)
CRC_12_3GPP =        spec(12, 0x080f, 0x0000, false, true,  0x0000, 0x0daf)
CRC_12_CDMA2000 =    spec(12, 0x0f13, 0x0fff, false, false, 0x0000, 0x0d4d)
CRC_12_DECT =        spec(12, 0x080f, 0x0000, false, false, 0x0000, 0x0f5b)
CRC_13_BBC =         spec(13, 0x1cf5, 0x0000, false, false, 0x0000, 0x04fa)
CRC_14_DARC =        spec(14, 0x0805, 0x0000, true,  true,  0x0000, 0x082d)
CRC_15 =             spec(15, 0x4599, 0x0000, false, false, 0x0000, 0x059e)
CRC_15_MPT1327 =     spec(15, 0x6815, 0x0000, false, false, 0x0001, 0x2566)

# http://reveng.sourceforge.net/crc-catalogue/16.htm#crc.cat.crc-16-ccitt-false
CRC_16_ARC =         spec(0x8005, 0x0000, true,  true,  0x0000, 0xbb3d)
CRC_16_AUG_CCITT =   spec(0x1021, 0x1d0f, false, false, 0x0000, 0xe5cc)
CRC_16_BUYPASS   =   spec(0x8005, 0x0000, false, false, 0x0000, 0xfee8)
CRC_16_CCITT_FALSE = spec(0x1021, 0xffff, false, false, 0x0000, 0x29b1)
CRC_16_CDMA2000 =    spec(0xc867, 0xffff, false, false, 0x0000, 0x4c06)
CRC_16_DDS_110 =     spec(0x8005, 0x800d, false, false, 0x0000, 0x9ecf)
CRC_16_DECT_R =      spec(0x0589, 0x0000, false, false, 0x0001, 0x007e)
CRC_16_DECT_X =      spec(0x0589, 0x0000, false, false, 0x0000, 0x007f)
CRC_16_DNP =         spec(0x3d65, 0x0000, true,  true,  0xffff, 0xea82)
CRC_16_EN_13757 =    spec(0x3d65, 0x0000, false, false, 0xffff, 0xc2b7)
CRC_16_GENIBUS =     spec(0x1021, 0xffff, false, false, 0xffff, 0xd64e)
CRC_16_MAXIM =       spec(0x8005, 0x0000, true,  true,  0xffff, 0x44c2)
CRC_16_RIELLO =      spec(0x8bb7, 0x0000, false, false, 0x0000, 0xd0db)
CRC_16_TELEDISK =    spec(0x1021, 0x89ec, true,  true,  0x0000, 0x26b1)
CRC_16_USB =         spec(0x8005, 0xffff, true,  true,  0xffff, 0xb4c8)
CRC_16_CRC_A =       spec(0x1021, 0xc6c6, true,  true,  0x0000, 0xbf05)
CRC_16_KERMIT =      spec(0x1021, 0x0000, true,  true,  0x0000, 0x2189)
CRC_16_MODBUS =      spec(0x8005, 0xffff, true,  true,  0x0000, 0x4b37)
CRC_16_X_25 =        spec(0x1021, 0xffff, true,  true,  0xffff, 0x906e)
CRC_16_XMODEM =      spec(0x1021, 0x0000, false, false, 0x0000, 0x31c3)

# http://reveng.sourceforge.net/crc-catalogue/17plus.htm
CRC_24 =             spec(24, 0x00864cfb, 0x00b704ce, false, false, 0x00000000, 0x0021cf02)
CRC_24_FLEXRAY_A =   spec(24, 0x005d6dcb, 0x00fedcba, false, false, 0x00000000, 0x007979bd)
CRC_24_FLEXRAY_B =   spec(24, 0x005d6dcb, 0x00abcdef, false, false, 0x00000000, 0x001f23b8)
CRC_31_PHILIPS =     spec(31, 0x04c11db7, 0x7fffffff, false, false, 0x7fffffff, 0x0ce9e46c)
CRC_32 =             spec(0x04c11db7, 0xffffffff, true,  true,  0xffffffff, 0xcbf43926)
CRC_32_BZIP2 =       spec(0x04c11db7, 0xffffffff, false, false, 0xffffffff, 0xfc891918)
CRC_32_C =           spec(0x1edc6f41, 0xffffffff, true,  true,  0xffffffff, 0xe3069283)
CRC_32_D =           spec(0xa833982b, 0xffffffff, true,  true,  0xffffffff, 0x87315576)
CRC_32_MPEG_2 =      spec(0x04c11db7, 0xffffffff, false, false, 0x00000000, 0x0376e6e7)
CRC_32_POSIX =       spec(0x04c11db7, 0x00000000, false, false, 0xffffffff, 0x765e7680)
CRC_32_Q =           spec(0x814141ab, 0x00000000, false, false, 0x00000000, 0x3010bf7f)
CRC_32_JAMCRC =      spec(0x04c11db7, 0xffffffff, true,  true,  0x00000000, 0x340bc6d9)
CRC_32_XFER =        spec(0x000000af, 0x00000000, false, false, 0x00000000, 0xbd0be338)
CRC_40_GSM =         spec(40, 0x0000000004820009, 0x0000000000000000, false, false, 0x000000ffffffffff, 0x000000d4164fc646)
CRC_64 =             spec(0x42f0e1eba9ea3693, 0x0000000000000000, false, false, 0x0000000000000000, 0x6c40df5f0b497347)
CRC_64_WE =          spec(0x42f0e1eba9ea3693, 0xffffffffffffffff, false, false, 0xffffffffffffffff, 0x62ec59e3f1a4f00a)
CRC_64_XZ =          spec(0x42f0e1eba9ea3693, 0xffffffffffffffff, true,  true,  0xffffffffffffffff, 0x995dc9bbdf1939fa)
CRC_82_DARC =        spec(82, 0x0308c0111011401440411, 0x000000000000000000000, true,  true,   0x000000000000000000000, 0x09ea83f625023801fd612)

ALL = Dict{Any,Any}(:CRC_3_ROHC=>CRC_3_ROHC,
                            :CRC_4_ITU=>CRC_4_ITU,
                            :CRC_5_EPC=>CRC_5_EPC,
                            :CRC_5_ITU=>CRC_5_ITU,
                            :CRC_5_USB=>CRC_5_USB,
                            :CRC_6_CDMA2000_A=>CRC_6_CDMA2000_A,
                            :CRC_6_CDMA2000_B=>CRC_6_CDMA2000_B,
                            :CRC_6_DARC=>CRC_6_DARC,
                            :CRC_6_ITU=>CRC_6_ITU, :CRC_7=>CRC_7,
                            :CRC_7_ROHC=>CRC_7_ROHC, :CRC_8=>CRC_8,
                            :CRC_8_CDMA2000=>CRC_8_CDMA2000,
                            :CRC_8_DARC=>CRC_8_DARC,
                            :CRC_8_DVB_S2=>CRC_8_DVB_S2,
                            :CRC_8_EBU=>CRC_8_EBU,
                            :CRC_8_I_CODE=>CRC_8_I_CODE,
                            :CRC_8_ITU=>CRC_8_ITU,
                            :CRC_8_MAXIM=>CRC_8_MAXIM,
                            :CRC_8_ROHC=>CRC_8_ROHC,
                            :CRC_8_WCDMA=>CRC_8_WCDMA,
                            :CRC_10=>CRC_10,
                            :CRC_10_CDMA2000=>CRC_10_CDMA2000,
                            :CRC_11=>CRC_11,
                            :CRC_12_3GPP=>CRC_12_3GPP,
                            :CRC_12_CDMA2000=>CRC_12_CDMA2000,
                            :CRC_12_DECT=>CRC_12_DECT,
                            :CRC_13_BBC=>CRC_13_BBC,
                            :CRC_14_DARC=>CRC_14_DARC,
                            :CRC_15=>CRC_15,
                            :CRC_15_MPT1327=>CRC_15_MPT1327,
                            :CRC_16_ARC=>CRC_16_ARC,
                            :CRC_16_AUG_CCITT=>CRC_16_AUG_CCITT,
                            :CRC_16_BUYPASS=>CRC_16_BUYPASS,
                            :CRC_16_CCITT_FALSE=>CRC_16_CCITT_FALSE,
                            :CRC_16_CDMA2000=>CRC_16_CDMA2000,
                            :CRC_16_DDS_110=>CRC_16_DDS_110,
                            :CRC_16_DECT_R=>CRC_16_DECT_R,
                            :CRC_16_DECT_X=>CRC_16_DECT_X,
                            :CRC_16_DNP=>CRC_16_DNP,
                            :CRC_16_EN_13757=>CRC_16_EN_13757,
                            :CRC_16_GENIBUS=>CRC_16_GENIBUS,
                            :CRC_16_MAXIM=>CRC_16_MAXIM,
                            :CRC_16_RIELLO=>CRC_16_RIELLO,
                            :CRC_16_TELEDISK=>CRC_16_TELEDISK,
                            :CRC_16_USB=>CRC_16_USB,
                            :CRC_16_CRC_A=>CRC_16_CRC_A,
                            :CRC_16_KERMIT=>CRC_16_KERMIT,
                            :CRC_16_MODBUS=>CRC_16_MODBUS,
                            :CRC_16_X_25=>CRC_16_X_25,
                            :CRC_16_XMODEM=>CRC_16_XMODEM,
                            :CRC_24=>CRC_24,
                            :CRC_24_FLEXRAY_A=>CRC_24_FLEXRAY_A,
                            :CRC_24_FLEXRAY_B=>CRC_24_FLEXRAY_B,
                            :CRC_31_PHILIPS=>CRC_31_PHILIPS,
                            :CRC_32=>CRC_32,
                            :CRC_32_BZIP2=>CRC_32_BZIP2,
                            :CRC_32_C=>CRC_32_C, :CRC_32_D=>CRC_32_D,
                            :CRC_32_MPEG_2=>CRC_32_MPEG_2,
                            :CRC_32_POSIX=>CRC_32_POSIX,
                            :CRC_32_Q=>CRC_32_Q,
                            :CRC_32_JAMCRC=>CRC_32_JAMCRC,
                            :CRC_32_XFER=>CRC_32_XFER,
                            :CRC_40_GSM=>CRC_40_GSM, :CRC_64=>CRC_64,
                            :CRC_64_WE=>CRC_64_WE,
                            :CRC_64_XZ=>CRC_64_XZ,
                            :CRC_82_DARC=>CRC_82_DARC)

# --- Utilities

# http://stackoverflow.com/questions/2602823/in-c-c-whats-the-simplest-way-to-reverse-the-order-of-bits-in-a-byte
function reflect_bits(n::UInt8)
    n = (n & 0xf0) >>> 4 | (n & 0x0f) << 4;
    n = (n & 0xcc) >>> 2 | (n & 0x33) << 2;
    (n & 0xaa) >>> 1 | (n & 0x55) << 1;
end

const REFLECT_8 = UInt8[reflect_bits(i) for i in 0x00:0xff]

reflect(u::UInt8) = REFLECT_8[u+1]

for (T,S) in ((UInt16, UInt8), (UInt32, UInt16), (UInt64, UInt32), (UInt128, UInt64))
    n = 8 * sizeof(S)
    mask::S = -one(S)
    @eval reflect(u::$T) = (convert($T, reflect(convert($S, u & $mask))) << $n) | reflect(convert($S, (u >>> $n) & $mask))
end

function reflect{T<:U}(size, u::T)
    width = 8 * sizeof(T)
    @assert size <= width "cannot reflect a value larger than the representation"
    reflect(u) >>> (width - size)
end

function largest(T, TS...)
    for t in TS
        if sizeof(t) > sizeof(T)
            T = t
        end
    end
    T
end

function fastest(T, TS...)
    L = largest(T, TS...)
    if L == UInt32 && UInt32 != UInt
        UInt64  # round up for speed
    else
        L
    end
end

function pad{A<:U}(::Type{A}, width)
    8 * sizeof(A) - width
end

# --- CRC calculations

# a calculation may be a "natural" long division (padded to the size
# of the accumulator) or "reverse the rest of the world".  this
# depends on whether the input input bytes are taken as given
# (Forwards) or reflected (Backwards).

abstract Direction{A<:U}

immutable Backwards{A<:U}<:Direction{A}
    poly::A
    init::A
    function Backwards(spec::Spec)
        poly = reflect(spec.width, convert(A, spec.poly))
        init = reflect(spec.width, convert(A, spec.init))
        new(poly, init)
    end
end

immutable Forwards{A<:U}<:Direction{A}
    pad_p::Int
    poly::A
    carry::A
    pad_8::Int
    init::A
    function Forwards(spec::Spec)
        pad_p = pad(A, spec.width)
        poly = convert(A, spec.poly) << pad_p
        carry = one(A) << pad(A, 1)
        pad_8 = pad(A, 8)
        init = convert(A, spec.init) << pad_p
        new(pad_p, poly, carry, pad_8, init)
    end
end

# a calculation may use no, one, or many tables...

abstract Tables{A<:U}

immutable NoTables{A<:U}<:Tables{A}
    NoTables(direcn::Direction{A}) = new()
end

immutable Multiple{A<:U}<:Tables{A}
    tables::Vector{Vector{A}}
    function Multiple(direcn::Direction{A})
        tables = Vector{A}[Array(A, 256) for _ in 1:sizeof(A)]
        fill_table(direcn, tables[1])
        for index in zero(UInt8):convert(UInt8, 255)
            remainder = tables[1][index + 1]
            for t in 2:sizeof(A)
                # chaining gives the contribution of a byte to the
                # remainder after being shifted through t other zero bytes
                # (the different bytes are then combined with xor).
                remainder = chain(direcn, tables[1], remainder)
                tables[t][index + 1] = remainder
            end
        end
        new(tables)
    end
end

immutable Single{A<:U}<:Tables{A}
    table::Vector{A}
    # when using multiple tables we need a related single table for the
    # "last few bytes"
    Single(tables::Multiple{A}) = new(tables.tables[1])
    Single(direcn::Direction{A}) = new(fill_table(direcn, Array(A, 256)))
end

# the basic lookup table is just the remainder for each value
function fill_table{A<:U}(direcn, table::Vector{A})
    for index in zero(UInt8):convert(UInt8, 255)
        table[index + 1] = extend(direcn, NoTables{A}(direcn), [index], zero(A))
    end
    table
end

function chain{A<:U}(direcn::Forwards{A}, table::Vector{A}, remainder::A)
    (remainder << 8) $ table[1 + ((remainder >>> direcn.pad_8) & 0xff)]
end

function chain{A<:U}(direcn::Backwards{A}, table::Vector{A}, remainder::A)
    (remainder >>> 8) $ table[1 + (remainder & 0xff)]
end

# main entry point.
function crc{P<:U, T<:Tables}(spec::Spec{P}; tables::Type{T}=Multiple)
    A = fastest(P, tables == Multiple ? UInt : UInt8)
    direcn = spec.refin ? Backwards{A}(spec) : Forwards{A}(spec)
    remainder::A = direcn.init
    tables = tables{A}(direcn)
    function handler(data::Vector{UInt8}; append=false)
        remainder = append ? remainder : direcn.init
        remainder = extend(direcn, tables, data, remainder)
        finalize(spec, direcn, remainder)
    end
    function handler(io::IO; append=false, buflen=1000000)
        buffer = Array(UInt8, buflen)
        remainder = append ? remainder : direcn.init
        while (nb = readbytes!(io, buffer)) > 0
            remainder = extend(direcn, tables, buffer[1:nb], remainder)
        end
        finalize(spec, direcn, remainder)
    end
    function handler(data::AbstractString; append=false)
        handler(convert(Vector{UInt8}, data), append=append)
    end
    handler
end

# calculations assume that the remainder is either reflected or
# padded; these functions correct for this to return the final result.

function finalize{P<:U, A<:U}(spec::Spec{P}, direcn::Forwards{A}, remainder::A)
    remainder = convert(P, remainder >> direcn.pad_p)
    remainder = spec.refout ? reflect(spec.width, remainder) : remainder
    remainder $ spec.xorout
end

function finalize{P<:U, A<:U}(spec::Spec{P}, direcn::Backwards{A}, remainder::A)
    remainder = spec.refout ? remainder : reflect(spec.width, remainder)
    convert(P, remainder) $ spec.xorout
end

# all the routines below calculate a new value of remainder for the
# data given (along with the algorithm type, available lookup tables,
# etc).

const UNROLL = 16  # only get small improvements past this

function extend{A<:U}(direcn::Forwards{A}, tables::NoTables, data::Vector{UInt8}, remainder::A)
    for word::UInt8 in data
        remainder::A = remainder $ (convert(A, word) << direcn.pad_8)
        for _ in 1:8
            if remainder & direcn.carry == direcn.carry
                remainder = (remainder << 1) $ direcn.poly
            else
                remainder <<= 1
            end
        end
    end
    remainder
end

function extend{A<:U}(direcn::Forwards{A}, tables::Single{A}, data::Vector{UInt8}, remainder::A)
    for i in 1:length(data)
        @inbounds word::UInt8 = data[i]
        remainder::A = remainder $ (convert(A, word) << direcn.pad_8)
        remainder = (remainder << 8) $ tables.table[1 + remainder >>> direcn.pad_8]
    end
    remainder
end

for A in (UInt16, UInt32, UInt64, UInt128)
    n_tables = sizeof(A)
    @eval begin
        function extend(direcn::Forwards{$A}, tables::Multiple{$A}, data::Vector{$A}, remainder::$A)
            word::$A, tmp::$A, remainder::$A = zero($A), zero($A), remainder
            i = 1
            while true
                @nexprs $UNROLL _ -> begin  # unroll inner loop
                    if i > length(data)
                        break
                    end
                    @inbounds word = data[i]
                    tmp, remainder = remainder, zero($A)
                    @nexprs $n_tables t -> begin  # unroll table
                        # access similar to reflected but we need to
                        # handle little-endian bytes which turn out
                        # wrong for this case
                        remainder $= tables.tables[t][((word >>> ($n_tables - t)*8) $ (tmp >>> (t-1)*8)) & 0xff + 1]
                    end
                    i += 1
                end
            end
            remainder
        end
    end
end

function extend{A<:U}(direcn::Backwards{A}, tables::NoTables, data::Vector{UInt8}, remainder::A)
    for word::UInt8 in data
        remainder::A = remainder $ convert(A, word)
        for _ in 1:8
            if remainder & one(A) == one(A)
                remainder = (remainder >>> 1) $ direcn.poly
            else
                remainder >>>= 1
            end
        end
    end
    remainder
end

function extend{A<:U}(direcn::Backwards{A}, tables::Single{A}, data::Vector{UInt8}, remainder::A)
    for i in 1:length(data)
        @inbounds word::UInt8 = data[i]
        remainder::A = remainder $ (convert(A, word))
        remainder = (remainder >>> 8) $ tables.table[1 + remainder & 0xff]
    end
    remainder
end

# generate code for processing a bunch of bytes in a single machine
# word, using multiple tables.  stolen from libz.
for A in (UInt16, UInt32, UInt64, UInt128)
    n_tables = sizeof(A)
    @eval begin
        function extend(direcn::Backwards{$A}, tables::Multiple{$A}, data::Vector{$A}, remainder::$A)
            word::$A, tmp::$A, remainder::$A = zero($A), zero($A), remainder
            i = 1
            while true
                @nexprs $UNROLL _ -> begin  # unroll inner loop
                    if i > length(data)
                        break
                    end
                    @inbounds word = data[i]
                    tmp, remainder = remainder $ word, zero($A)
                    @nexprs $n_tables t -> begin  # unroll table access
                        remainder $= tables.tables[t][(tmp >>> ($n_tables - t)*8) & 0xff + 1]
                    end
                    i += 1
                end
            end
            remainder
        end
    end
end

# short-circuit "all UInt8" avoiding a self-recursive loop below
function extend(direcn::Direction{UInt8}, tables::Multiple{UInt8}, data::Vector{UInt8}, remainder::UInt8)
    extend(direcn, Single{UInt8}(tables), data, remainder)
end

function extend{A<:U}(direcn::Direction{A}, tables::Multiple{A}, data::Vector{UInt8}, remainder::A)
    # this is "clever" - we alias the array of bytes to native machine
    # words and then process those, which typically (64 bits) loads 8
    # bytes at a time.
    tail = length(data) % sizeof(A)
    remainder = extend(direcn, tables, reinterpret(A, data), remainder)
    # slurp up the final bytes that didn't fill a complete machine word.
    extend(direcn, Single{A}(tables), data[end-tail+1:end], remainder)
end

# ---- Command line

# julia -e "using CRC; main(ARGS)" ...

function main(args)

    s = ArgParseSettings("Calculate the CRC for the given files")

    @add_arg_table s begin
        "--list", "-l"
        help = "list available CRC algorithms"
        action = :store_true
        "--decimal", "-d"
        help = "show checksums as decimal values (default is hex)"
        action = :store_true
        "--crc", "-c"
        help = "name the CRC to use"
        default = "CRC_32"
        "--append", "-a"
        help = "combine the data from all files"
        action = :store_true
        "files"
        help = "the files to read (- for stdin)"
        nargs = '*'
    end

    parsed_args = parse_args(args, s)

    if parsed_args["list"]
        names = sort(collect(keys(ALL)), by=n->ALL[n])
        for name in names
            @printf("%s %s\n", name, ALL[name])
        end
    end

    name = symbol(parsed_args["crc"])
    if !haskey(ALL, name)
        error("CRC $name is not defined")
    end

    append = parsed_args["append"]
    spec = ALL[name]
    c = crc(spec)
    if parsed_args["decimal"]
        fmt = x -> dec(x)
    else
        fmt = x -> "0x" * hex(x, 1 + div(spec.width-1, 4))
    end
    sum = 0
    for file in parsed_args["files"]
        if file == "-"
            sum = c(STDIN, append=append)
        else
            open(file, "r") do f
                sum = c(f, append=append)
            end
        end
        if !append
            println("$(fmt(sum)) $file")
        end
    end
    if append
        println("$(fmt(sum))")
    end

end

end


text:
module SSLTestConstruction
using SimpleSDMLayers
using Test

# Correct type
M = rand(Bool, (5,5))
S = SimpleSDMPredictor(M, 0.0, 1.0, 0.0, 1.0)
@test typeof(S) == SimpleSDMPredictor{Bool}

# Correct type
M = rand(Bool, (5,5))
R = SimpleSDMResponse(M, 0.0, 1.0, 0.0, 1.0)
@test typeof(R) == SimpleSDMResponse{Bool}

# Construction only from a matrix
R = SimpleSDMResponse(M)
@test typeof(R) == SimpleSDMResponse{Bool}
@test R.left == -180.
@test R.right == 180.
@test R.top == 90.
@test R.bottom == -90.

R = SimpleSDMPredictor(M)
@test typeof(R) == SimpleSDMPredictor{Bool}
@test R.left == -180.
@test R.right == 180.
@test R.top == 90.
@test R.bottom == -90.

# Construction for a matrix of nothing and values
M = Matrix{Union{Nothing,Float64}}(nothing, (3,5))
R = SimpleSDMResponse(M)
@test R.left == -180.
@test R.right == 180.
@test R.top == 90.
@test R.bottom == -90.

end

ext: jl


text:
#
# This script runs all tests in current folder, which are
# in format: "test-<name>.jl". Other files will be skipped.
#
# @author DeadbraiN
#
using FactCheck

#
# All tests should have format: "test-<name>.jl"
#
test_regex = r"^test-.*\.jl$"
#
# This line is needed after migrating to julia v0.4, because
# include() doesn't search files(modules) in current folder.
#
include("../src/global/ImportFolders.jl")
#
# All test files should be in "test" folder
#
test_files = filter(n -> ismatch(test_regex, n), readdir("test"))
if length(test_files) == 0
    error("No test files found. Make sure you're running from the root directory and test files are in 'test' folder.")
end
#
# Include tests means run them
#
for test_file in test_files
    println("\n\nRunning $test_file")
    include(test_file)
end
#
# return the overall exit status
#
exitstatus()

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
using WebIO
using Test
using Blink, JSExpr

if !(@isdefined open_window)
    include("./test-utils.jl")
end

@testset "Test observable updates are sent to all connections" begin
    w1, w2 = open_window(), open_window()
    s = Scope(
        dom=node(
            :p,
            "initial",
            attributes=Dict(
                "data-testid" => "foobar",
            ),
        ),
    )
    obs = Observable{Any}(s, "obs", "")
    onjs(obs, @js function (value)
        _webIOScope.dom.querySelector("[data-testid=foobar]").innerText = value
    end)

    body!(w1, s)
    body!(w2, s)
    # Sleep to allow the JS to execute.
    sleep(0.25)

    value = (w) -> @js w document.querySelector("[data-testid=foobar]").innerText
    @test value(w1) == "initial"
    @test value(w2) == "initial"

    obs[] = "updated"
    # Sleep to allow the JS to execute.
    sleep(0.25)
    @test value(w1) == "updated"
    @test value(w2) == "updated"
end


text:
immutable Skellam <: DiscreteUnivariateDistribution
    mu1::Float64
    mu2::Float64
    Skellam(m1::Real, m2::Real) = new(float64(m1), float64(m2))
end

isupperbounded(::Union(Skellam, Type{Skellam})) = false
islowerbounded(::Union(Skellam, Type{Skellam})) = false
isbounded(::Union(Skellam, Type{Skellam})) = false

minimum(::Union(Skellam, Type{Skellam})) = -Inf
maximum(::Union(Skellam, Type{Skellam})) = Inf

insupport(::Skellam, x::Real) = isinteger(x) && zero(x) <= x
insupport(::Type{Skellam}, x::Real) = isinteger(x) && zero(x) <= x

kurtosis(d::Skellam) = 1.0 / (d.mu1 + d.mu2)

mean(d::Skellam) = d.mu1 - d.mu2

function mgf(d::Skellam, t::Real)
    exp(-(d.mu1 + d.mu2) + d.mu1 * exp(t) + d.mu2 * exp(-t))
end

function cf(d::Skellam, t::Real)
    exp(-(d.mu1 + d.mu2) + d.mu1 * exp(im * t) + d.mu2 * exp(im * -t))
end

function pdf(d::Skellam, k::Real)
    isinteger(k) || return 0.0
    exp(-(d.mu1 + d.mu2)) * (d.mu1 / d.mu2)^(k / 2.0) *
        real(besseli(k, 2.0 * sqrt(d.mu1 * d.mu2)))
end

function logpdf(d::Skellam, k::Real)
    isinteger(k) || return 0.0
     -(d.mu1 + d.mu2) + (k / 2.0) * log(d.mu1) - log(d.mu2) +
           log(real(besseli(k, 2.0 * sqrt(d.mu1 * d.mu2))))
end

rand(d::Skellam) = rand(Poisson(d.mu1)) - rand(Poisson(d.mu2))

skewness(d::Skellam) = (d.mu1 - d.mu2) / (d.mu1 + d.mu2)^1.5

var(d::Skellam) = d.mu1 + d.mu2


text:
module Utils

field_name_pairs(value, field_names, args...) = NamedTuple(name => value for name in field_names)
field_name_pairs(t::Union{Tuple, AbstractArray}, field_names, args...) =
    NamedTuple(name => t[i] for (i, name) in enumerate(field_names))

function field_name_pairs(nt::Union{NamedTuple, Dict}, field_names, nt_name="")
    nt_field_names = keys(nt)
    nt_summary = summary(nt)

    # Validate user-supplied NamedTuple
    all(name ∈ nt_field_names for name in field_names) ||
        throw(ArgumentError("$nt_name $nt_summary must have values for every field in $field_names " *
                            " but only has values for $nt_field_names"))

    return nt
end

function prettyvector(v::AbstractVector, bookends=3)
    separator = " … "
    N = length(v)

    if N < 2bookends + 4
        content = string.(v) .* ", "
        content[end] = content[end][1:end-2]
        return string("[", content..., "]")
    else
        separator = " … "
        beginning = [string(v[i]) for i=1:bookends]
        ending = [string(v[end+1-i]) for i=1:bookends]

        for i = 1:bookends-1
            beginning[i] *= ", "
            ending[i] *= ", "
        end

        N = length(v)

        return string("[", beginning..., separator, ending..., "] ($N elements)")
    end
end

end # module


text:
using Measurements
using Printf
include("/home/jefter66/Projects/Lablib/src/Stats.jl")
include("/home/jefter66/Projects/Lablib/src/LeastSquares.jl")
γ(h1,h3) = return (h1) / (h1 - h3);

Δh = 0.1;

h1 = [18.1 ± Δh, 16.0 ± Δh, 17.6 ± Δh];
h3 = [5.7, 4,5.1];

ɣs = [γ(h1[1],h3[1]), γ(h1[2],h3[2]), γ(h1[3],h3[3])];
deviation(ɣs);

A(d) = return π * (d^2 / 4);
# ΔA = 0.00001;

P0 = 691 * 133.322 ± 0.01; # ΔP = 0.01;
g = 9.81;

d = 0.016 ± 0.001; # Δd = 0.001;
m = 0.0167 ± 0.0001; # Δm = 0.0001;
V = 0.0104  ± 0.0001; # ΔV = 0.0001;
P = (m * g)/A(d) + P0;
T = 1.2  ± 0.01; # ΔT = 0.01;

𝛄 = 4 * π^2 * ((m * V)/(P * A(d)^2 * T^2));
#x, Δx = m * V, Δm * V + ΔV * m;
#y, Δy = P*A(d)^2 * T^2, 2*T*ΔT*(P*A(d)^2)  +  (ΔP * A(d)^2 + 2 * A(d) * ΔA * P);
#Δ𝛄 = 4 * π^2 * (Δx * y + Δy * x) / (y^2);

@printf("𝛄 = %.6f ± %.6f", 𝛄.val, 𝛄.err);

x = [24.7,1.0,-196.0,97.0];
y = [69.2,62.8,17.7,82.3];

leastSquares(x,y);

max_stars_repo_path: Fisica experimental/02/06/p6-calcs.jl

max_issues_repo_path: Fisica experimental/02/06/p6-calcs.jl

max_forks_repo_path: Fisica experimental/02/06/p6-calcs.jl


text:
using PropertyGraphs
using LightGraphs

mutable struct VertexProperties
    color::String
    size::Float64
end

# -------- String vertex labels. --------

pg = SimplePropertyGraph{Int, String, VertexProperties}()

add_vertex!(pg, "a")
add_vertex!(pg, "b")
add_vertex!(pg, "c")

"b" in pg
"asdf" in pg
"b" ∉ pg
"asdf" ∉ pg

add_edge!(pg, "a", "b")
add_edge!(pg, "b", "c")
add_edge!(pg, pg("a", "c"))

pg["a"] = VertexProperties("red", 3.4)
pg["b"] = VertexProperties("blue", 1.9)
pg["c"] = VertexProperties("yellow", 7.0)

pg["a"]
pg["b"]
pg["c"]

pg("a")
pg("b")
pg("c")
pg("a", "b")

label(pg, 1)
label(pg, 2)
label(pg, 3)
label.(pg, [1, 2, 3])

pg["a"].color
pg["a"].color = "pink"
pg["a"].color

pg["b"].size
pg["b"].size = 100.2
pg["b"].size

paths = dijkstra_shortest_paths(pg, pg("a"))
enumerate_paths(paths, pg("c"))

# -------- Integer vertex labels. --------

pg2 = SimplePropertyGraph{Int, Int, VertexProperties}()

add_vertex!(pg2, 101)
add_vertex!(pg2, 102)
add_vertex!(pg2, 103)

add_edge!(pg2, 101, 102)
add_edge!(pg2, 102, 103)
add_edge!(pg2, 101, 103)

pg2[101] = VertexProperties("red", 3.4)
pg2[102] = VertexProperties("blue", 1.9)
pg2[103] = VertexProperties("yellow", 7.0)

pg2[101]
pg2[102]
pg2[103]

# -------- Directed graph. --------

pg = SimplePropertyDiGraph{Int, String, VertexProperties}()

add_vertex!(pg, "a")
add_vertex!(pg, "b")
add_vertex!(pg, "c")

add_edge!(pg, "a", "b")
add_edge!(pg, "b", "c")
add_edge!(pg, pg("a", "c"))

pg["a"] = VertexProperties("red", 3.4)
pg["b"] = VertexProperties("blue", 1.9)
pg["c"] = VertexProperties("yellow", 7.0)

paths = dijkstra_shortest_paths(pg, pg("a"))
enumerate_paths(paths, pg("c"))

# -------- Remove vertex. --------

pg = SimplePropertyDiGraph{Int, String, VertexProperties}()

add_vertex!(pg, "a")
add_vertex!(pg, "b")
add_vertex!(pg, "c")
add_vertex!(pg, "d")
add_vertex!(pg, "e")

pg["a"] = VertexProperties("red", 3.4)
pg["b"] = VertexProperties("blue", 1.9)
pg["c"] = VertexProperties("yellow", 7.0)
pg["d"] = VertexProperties("pink", 42)
pg["e"] = VertexProperties("orange", 100)

add_edge!(pg, "a", "b")
add_edge!(pg, "b", "c")
add_edge!(pg, "c", "d")
add_edge!(pg, "d", "e")
add_edge!(pg, "e", "a")

write_dot("/Users/bieganek/test.dot", pg)

rem_vertex!(pg, "b")


text:
###
### Stochastic Gradient Hamiltonian Samplers
###

####
#### Stochastic Gradient Hamiltonian Monte Carlo sampler.
####

"""
    SGHMC(n_iters::Int, learning_rate::Float64, momentum_decay::Float64)

Stochastic Gradient Hamiltonian Monte Carlo sampler.

Usage:

```julia
SGHMC(1000, 0.01, 0.1)
```

Arguments:

- `n_iters::Int` : Number of samples to pull.
- `learning_rate::Float64` : The learning rate.
- `momentum_decay::Float64` : Momentum decay variable.

"""
mutable struct SGHMC{AD, space, metricT <: AHMC.AbstractMetric} <: StaticHamiltonian{AD}
    n_iters::Int       # number of samples
    learning_rate::Float64   # learning rate
    momentum_decay::Float64   # momentum decay
end
SGHMC(args...) = SGHMC{ADBackend()}(args...)
function SGHMC{AD}(n_iters::Int, learning_rate::Float64, momentum_decay::Float64, ::Type{metricT}, space::Tuple) where {AD, metricT <: AHMC.AbstractMetric}
    return SGHMC{AD, space, metricT}(n_iters, learning_rate, momentum_decay)
end

function SGHMC{AD}(
    n_iters,
    learning_rate,
    momentum_decay,
    ::Tuple{};
    kwargs...
) where AD
    return SGHMC{AD}(n_iters, learning_rate, momentum_decay; kwargs...)
end
function SGHMC{AD}(
    n_iters,
    learning_rate,
    momentum_decay,
    space::Symbol...;
    metricT=AHMC.UnitEuclideanMetric
) where AD
    return SGHMC{AD}(n_iters, learning_rate, momentum_decay, metricT, space)
end

function step(
    model,
    spl::Sampler{<:SGHMC},
    vi::VarInfo,
    is_first::Val{true};
    kwargs...
)
    spl.selector.tag != :default && link!(vi, spl)

    # Initialize velocity
    v = zeros(Float64, size(vi[spl]))
    spl.info[:v] = v

    spl.selector.tag != :default && invlink!(vi, spl)
    return vi, true
end

function step(
    model,
    spl::Sampler{<:SGHMC},
    vi::VarInfo,
    is_first::Val{false};
    kwargs...
)
    # Set parameters
    η, α = spl.alg.learning_rate, spl.alg.momentum_decay
    spl.info[:eval_num] = 0

    Turing.DEBUG && @debug "X-> R..."
    if spl.selector.tag != :default
        link!(vi, spl)
        runmodel!(model, vi, spl)
    end

    Turing.DEBUG && @debug "recording old variables..."
    θ, v = vi[spl], spl.info[:v]
    _, grad = gradient_logp(θ, vi, model, spl)
    verifygrad(grad)

    # Implements the update equations from (15) of Chen et al. (2014).
    Turing.DEBUG && @debug "update latent variables and velocity..."
    θ .+= v
    v .= (1 - α) .* v .+ η .* grad .+ rand.(Normal.(zeros(length(θ)), sqrt(2 * η * α)))

    Turing.DEBUG && @debug "saving new latent variables..."
    vi[spl] = θ

    Turing.DEBUG && @debug "R -> X..."
    spl.selector.tag != :default && invlink!(vi, spl)

    Turing.DEBUG && @debug "always accept..."
    return vi, true
end

####
#### Stochastic Gradient Langevin Dynamics sampler.
####

"""
    SGLD(n_iters::Int, ϵ::Float64)

 Stochastic Gradient Langevin Dynamics sampler.

Usage:

```julia
SGLD(1000, 0.5)
```

Arguments:

- `n_iters::Int` : Number of samples to pull.
- `ϵ::Float64` : The scaling factor for the learing rate.

Reference:

Welling, M., & Teh, Y. W. (2011).  Bayesian learning via stochastic gradient Langevin dynamics.
In Proceedings of the 28th international conference on machine learning (ICML-11) (pp. 681-688).
"""
mutable struct SGLD{AD, space, metricT <: AHMC.AbstractMetric} <: StaticHamiltonian{AD}
    n_iters :: Int       # number of samples
    ϵ :: Float64   # constant scale factor of learning rate
end
function SGLD{AD}(n_iters::Int, ϵ::Float64, ::Type{metricT}, space::Tuple) where {AD, metricT <: AHMC.AbstractMetric}
    return SGLD{AD, space, metricT}(n_iters, ϵ)
end

SGLD(args...; kwargs...) = SGLD{ADBackend()}(args...; kwargs...)

function SGLD{AD}(
    n_iters,
    ϵ,
    ::Tuple{};
    kwargs...
) where AD
    return SGLD{AD}(n_iters, ϵ; kwargs...)
end

function SGLD{AD}(
    n_iters,
    ϵ,
    space::Symbol...;
    metricT=AHMC.UnitEuclideanMetric
) where AD
    return SGLD{AD}(n_iters, ϵ, metricT, space)
end

function step(
    model,
    spl::Sampler{<:SGLD},
    vi::VarInfo,
    is_first::Val{true};
    kwargs...
)
    spl.selector.tag != :default && link!(vi, spl)

    mssa = AHMC.Adaptation.ManualSSAdaptor(AHMC.Adaptation.MSSState(spl.alg.ϵ))
    spl.info[:adaptor] = AHMC.NaiveHMCAdaptor(AHMC.UnitPreconditioner(), mssa)

    spl.selector.tag != :default && invlink!(vi, spl)
    return vi, true
end

function step(
    model,
    spl::Sampler{<:SGLD},
    vi::VarInfo,
    is_first::Val{false};
    kwargs...
)
    # Update iteration counter
    spl.info[:i] += 1
    spl.info[:eval_num] = 0

    Turing.DEBUG && @debug "compute current step size..."
    γ = .35
    ϵ_t = spl.alg.ϵ / spl.info[:i]^γ # NOTE: Choose γ=.55 in paper
    mssa = spl.info[:adaptor].ssa
    mssa.state.ϵ = ϵ_t

    Turing.DEBUG && @debug "X-> R..."
    if spl.selector.tag != :default
        link!(vi, spl)
        runmodel!(model, vi, spl)
    end

    Turing.DEBUG && @debug "recording old variables..."
    θ = vi[spl]
    _, grad = gradient_logp(θ, vi, model, spl)
    verifygrad(grad)

    Turing.DEBUG && @debug "update latent variables..."
    θ .+= ϵ_t .* grad ./ 2 .- rand.(Normal.(zeros(length(θ)), sqrt(ϵ_t)))

    Turing.DEBUG && @debug "always accept..."
    vi[spl] = θ

    Turing.DEBUG && @debug "R -> X..."
    spl.selector.tag != :default && invlink!(vi, spl)

    return vi, true
end


text:
const fp64max = floatmax(Float64)
const maxhyp_fp64max = 710.4758600739439
const minhyp_fp64max = 5.562684646268003e-309

function sinh(x::DoubleFloat{T}) where {T<:IEEEFloat}
    abs(HI(x)) > maxhyp_fp64max && return DoubleFloat{T}(copysign(Inf, HI(x)))
    !isfinite(HI(x)) && return x
    iszero(x) && return zero(x)
    signbit(x) && return -sinh(-x)

    y = exp(x) - exp(-x)
    z = DoubleFloat{T}(y.hi*0.5, y.lo*0.5)
    return z
end

function cosh(x::DoubleFloat{T}) where {T<:IEEEFloat}
    abs(HI(x)) > maxhyp_fp64max && return DoubleFloat{T}(Inf)
    !isfinite(x) && return abs(x)
    iszero(x) && return one(x)
    signbit(x) && return cosh(-x)

    y = exp(x) + exp(-x)
    z = DoubleFloat{T}(y.hi*0.5, y.lo*0.5)
    return z
end

function tanh(x::DoubleFloat{T}) where {T<:IEEEFloat}
    abs(HI(x)) > maxhyp_fp64max && return copysign(one(DoubleFloat{T}), x)
    return sinh(x) / cosh(x)
end
#=
function tanh(x::DoubleFloat{T}) where {T<:IEEEFloat}
    x < 0 && return -tanh(-x)
    iszero(x) && return zero(x)
    !isfinite(x) && return nan(typeof(x))
    epos = exp(x)
    eneg = exp(-x)
    z = (epos - eneg) / (epos + eneg)
    return z
end
=#

function csch(x::DoubleFloat{T}) where {T<:IEEEFloat}
    abs(HI(x)) > maxhyp_fp64max && return copysign(zero(DoubleFloat{T}), x)
    return inv(sinh(x))
end

function sech(x::DoubleFloat{T}) where {T<:IEEEFloat}
    abs(HI(x)) > maxhyp_fp64max && return zero(DoubleFloat{T})
    return inv(cosh(x))
end

function coth(x::DoubleFloat{T}) where {T<:IEEEFloat}
    abs(HI(x)) > maxhyp_fp64max && return copysign(one(DoubleFloat{T}), x)
    return cosh(x) / sinh(x)
end


text:
#=

tribe wrapper.

Ignacio Quintero Mächler

t(-_-t)

April 27 2017

=#

"""
    tribe(tree_file::String,
          data_file::String,
          out_file ::String;
          min_dt   ::Float64           = 0.01,
          niter    ::Int64             = 50_000,
          nburn    ::Int64             = 50_000,
          nthin    ::Int64             = 500,
          saveXY   ::Tuple{Bool,Int64} = (false, 1_000),
          saveDM   ::Tuple{Bool,Int64} = (false, 1_000),
          ωxprior  ::NTuple{2,Float64} = (0.,10.),
          ω1prior  ::NTuple{2,Float64} = (0.,10.),
          ω0prior  ::NTuple{2,Float64} = (0.,10.),
          σ²prior  ::Float64           = 1e-1,
          λprior   ::Float64           = 1e-1,
          weight   ::NTuple{5,Float64} = (0.15,0.05,0.02,0.02,5e-3),
          λ1i      ::Float64           = 1.0,
          λ0i      ::Float64           = 0.5,
          ωxi      ::Float64           = 0.0,
          ω1i      ::Float64           = 0.0,
          ω0i      ::Float64           = 0.0,
          fix_ωx   ::Bool              = false,
          fix_ω1   ::Bool              = false,
          fix_ω0   ::Bool              = false,
          delim    ::Char              = '\t',
          eol      ::Char              = '\r')

Run tribe. Wrapper for all functions.
"""
function tribe(tree_file::String,
               data_file::String,
               out_file ::String;
               min_dt   ::Float64           = 0.01,
               niter    ::Int64             = 50_000,
               nburn    ::Int64             = 50_000,
               nthin    ::Int64             = 500,
               saveXY   ::Tuple{Bool,Int64} = (false, 1_000),
               saveDM   ::Tuple{Bool,Int64} = (false, 1_000),
               ωxprior  ::NTuple{2,Float64} = (0.,10.),
               ω1prior  ::NTuple{2,Float64} = (0.,10.),
               ω0prior  ::NTuple{2,Float64} = (0.,10.),
               σ²prior  ::Float64           = 1e-1,
               λprior   ::Float64           = 1e-1,
               weight   ::NTuple{5,Float64} = (0.15,0.05,0.02,0.02,5e-3),
               λ1i      ::Float64           = 1.0,
               λ0i      ::Float64           = 0.5,
               ωxi      ::Float64           = 0.0,
               ω1i      ::Float64           = 0.0,
               ω0i      ::Float64           = 0.0,
               fix_ωx   ::Bool              = false,
               fix_ω1   ::Bool              = false,
               fix_ω0   ::Bool              = false,
               delim    ::Char              = '\t',
               eol      ::Char              = '\r')

  tip_values, tip_areas, tree, bts =
    read_data_tribe(tree_file, data_file)

  X, Y, B, ncoup, δt, tree, si =
    initialize_data(tip_values, tip_areas, min_dt, tree, bts)

  R = tribe_mcmc(X, Y, ncoup, δt,
                 deepcopy(tree.ed),
                 deepcopy(tree.el),
                 B,
                 niter    = niter,
                 nthin    = nthin,
                 nburn    = nburn,
                 saveXY   = saveXY,
                 saveDM   = saveDM,
                 ωxprior  = ωxprior,
                 ω1prior  = ω1prior,
                 ω0prior  = ω0prior,
                 σ²prior  = σ²prior,
                 λprior   = λprior,
                 out_file = out_file,
                 weight   = weight,
                 λ1i      = λ1i,
                 λ0i      = λ0i,
                 ωxi      = ωxi,
                 ω1i      = ω1i,
                 ω0i      = ω0i,
                 σ²i      = si,
                 fix_ωx   = fix_ωx,
                 fix_ω1   = fix_ω1,
                 fix_ω0   = fix_ω0)

  return R

end

"""
    tribe(tip_values::Dict{Int64,Float64},
          tip_areas ::Dict{Int64,Array{Int64,1}},
          tree      ::rtree,
          bts       ::Array{Float64,1},
          out_file  ::String;
          min_dt    ::Float64           = 0.01,
          niter     ::Int64             = 500_000,
          nburn     ::Int64             = 500_000,
          nthin     ::Int64             = 1_000,
          saveXY    ::Tuple{Bool,Int64} = (false, 1_000),
          saveDM    ::Tuple{Bool,Int64} = (false, 1_000),
          ωxprior   ::NTuple{2,Float64} = (0.,10.),
          ω1prior   ::NTuple{2,Float64} = (0.,10.),
          ω0prior   ::NTuple{2,Float64} = (0.,10.),
          σ²prior   ::Float64           = 1e-1,
          λprior    ::Float64           = 1e-1,
          weight    ::NTuple{5,Float64} = (0.15,0.05,0.02,0.02,5e-3),
          λ1i       ::Float64           = 1.0,
          λ0i       ::Float64           = 0.4,
          ωxi       ::Float64           = 0.,
          ω1i       ::Float64           = 0.,
          ω0i       ::Float64           = 0.,
          fix_ωx    ::Bool              = false,
          fix_ω1    ::Bool              = false,
          fix_ω0    ::Bool              = false,
          delim     ::Char              = '\t',
          eol       ::Char              = '\r')

Run tribe for simulations. Wrapper for all functions.
"""
function tribe(tip_values::Dict{Int64,Float64},
               tip_areas ::Dict{Int64,Array{Int64,1}},
               tree      ::rtree,
               bts       ::Array{Float64,1},
               out_file  ::String;
               min_dt    ::Float64           = 0.01,
               niter     ::Int64             = 500_000,
               nburn     ::Int64             = 500_000,
               nthin     ::Int64             = 1_000,
               saveXY    ::Tuple{Bool,Int64} = (false, 1_000),
               saveDM    ::Tuple{Bool,Int64} = (false, 1_000),
               ωxprior   ::NTuple{2,Float64} = (0.,10.),
               ω1prior   ::NTuple{2,Float64} = (0.,10.),
               ω0prior   ::NTuple{2,Float64} = (0.,10.),
               σ²prior   ::Float64           = 1e-1,
               λprior    ::Float64           = 1e-1,
               weight    ::NTuple{5,Float64} = (0.15,0.05,0.02,0.02,5e-3),
               λ1i       ::Float64           = 1.0,
               λ0i       ::Float64           = 0.4,
               ωxi       ::Float64           = 0.,
               ω1i       ::Float64           = 0.,
               ω0i       ::Float64           = 0.,
               fix_ωx    ::Bool              = false,
               fix_ω1    ::Bool              = false,
               fix_ω0    ::Bool              = false,
               delim     ::Char              = '\t',
               eol       ::Char              = '\r')

  X, Y, B, ncoup, δt, tree, si =
    initialize_data(tip_values, tip_areas, min_dt, tree, bts)

  R = tribe_mcmc(X, Y, ncoup, δt,
                 deepcopy(tree.ed),
                 deepcopy(tree.el),
                 B,
                 niter    = niter,
                 nthin    = nthin,
                 nburn    = nburn,
                 saveXY   = saveXY,
                 saveDM   = saveDM,
                 ωxprior  = ωxprior,
                 ω1prior  = ω1prior,
                 ω0prior  = ω0prior,
                 σ²prior  = σ²prior,
                 λprior   = λprior,
                 out_file = out_file,
                 weight   = weight,
                 λ1i      = λ1i,
                 λ0i      = λ0i,
                 ωxi      = ωxi,
                 ω1i      = ω1i,
                 ω0i      = ω0i,
                 σ²i      = si,
                 fix_ωx   = fix_ωx,
                 fix_ω1   = fix_ω1,
                 fix_ω0   = fix_ω0)

  return R

end

"""
    tribe(out_file::String;
          niter   ::Int64             = 500_000,
          nburn   ::Int64             = 500_000,
          nthin   ::Int64             = 1_000,
          ωxprior ::NTuple{2,Float64} = (0.,10.),
          ω1prior ::NTuple{2,Float64} = (0.,10.),
          ω0prior ::NTuple{2,Float64} = (0.,10.),
          σ²prior ::Float64           = 1e-1,
          λprior  ::Float64           = 1e-1,
          weight  ::NTuple{4,Float64} = (0.15,0.05,0.02,0.02),
          σ²i     ::Float64           = 1.,
          ωxi     ::Float64           = 0.,
          ω1i     ::Float64           = 0.01,
          ω0i     ::Float64           = 0.01,
          λ1i     ::Float64           = 1.0,
          λ0i     ::Float64           = 0.2,
          fix_ωx  ::Bool              = false,
          fix_ω1  ::Bool              = false,
          fix_ω0  ::Bool              = false)

Run tribe **under the prior**. Wrapper for all functions.
"""
function tribe(out_file::String;
               niter   ::Int64             = 500_000,
               nburn   ::Int64             = 500_000,
               nthin   ::Int64             = 1_000,
               ωxprior ::NTuple{2,Float64} = (0.,10.),
               ω1prior ::NTuple{2,Float64} = (0.,10.),
               ω0prior ::NTuple{2,Float64} = (0.,10.),
               σ²prior ::Float64           = 1e-1,
               λprior  ::Float64           = 1e-1,
               weight  ::NTuple{4,Float64} = (0.15,0.05,0.02,0.02),
               σ²i     ::Float64           = 1.,
               ωxi     ::Float64           = 0.,
               ω1i     ::Float64           = 0.01,
               ω0i     ::Float64           = 0.01,
               λ1i     ::Float64           = 1.0,
               λ0i     ::Float64           = 0.2,
               fix_ωx  ::Bool              = false,
               fix_ω1  ::Bool              = false,
               fix_ω0  ::Bool              = false)

  R = tribe_mcmc(out_file,
                 niter    = niter,
                 nthin    = nthin,
                 nburn    = nburn,
                 ωxprior  = ωxprior,
                 ω1prior  = ω1prior,
                 ω0prior  = ω0prior,
                 σ²prior  = σ²prior,
                 λprior   = λprior,
                 weight   = weight,
                 λ1i      = λ1i,
                 λ0i      = λ0i,
                 ωxi      = ωxi,
                 ω1i      = ω1i,
                 ω0i      = ω0i,
                 σ²i      = σ²i,
                 fix_ωx   = fix_ωx,
                 fix_ω1   = fix_ω1,
                 fix_ω0   = fix_ω0)

  return R

end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
using ForwardDiff
using LinearAlgebra

let
	f(x) = x+x
	g(x) = 5*x
	h(x) = dot(x, x)

	# Once we get ForwardDiff working, this should work
	# h(x::VecE) = dot(x, x)

	@test ForwardDiff.jacobian(f, VecE2(1.0, 2.0)) == [2.0 0.0; 0.0 2.0]
	@test ForwardDiff.jacobian(g, VecE2(1.0, 2.0)) == [5.0 0.0; 0.0 5.0]
	@test ForwardDiff.gradient(h, VecE2(1.0, 2.0)) == [2.0, 4.0]
	@test ForwardDiff.hessian(h, VecE2(1.0, 2.0)) == [2.0 0.0; 0.0 2.0]

	@test ForwardDiff.jacobian(f, VecE3(1.0, 2.0, 3.0)) == 2.0*Matrix{Float64}(I, 3, 3)
	@test ForwardDiff.jacobian(g, VecE3(1.0, 2.0, 3.0)) == 5.0*Matrix{Float64}(I, 3, 3)
	@test ForwardDiff.gradient(h, VecE3(1.0, 2.0, 3.0)) == [2.0, 4.0, 6.0]
	@test ForwardDiff.hessian(h, VecE3(1.0, 2.0, 3.0)) == 2.0*Matrix{Float64}(I, 3, 3)

	@test ForwardDiff.jacobian(f, VecSE2(1.0, 2.0, 3.0)) == 2.0*Matrix{Float64}(I, 3, 3)

	# Once we get ForwardDiff working, this should work...
	# h(x::VecSE2) = normsquared(VecE2(x))
	# @test ForwardDiff.jacobian(g, VecSE2(1.0, 2.0, 3.0)) == diagm([5.0, 5.0, 1.0])
	# @test ForwardDiff.gradient(h, VecSE2(1.0, 2.0, 3.0)) == [2.0, 4.0, 0.0]
	# @test ForwardDiff.hessian(h, VecSE2(1.0, 2.0, 3.0)) == [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 0.0]
end


text:
@testset "Testing ecdf.jl" begin
    # set up

    obs = rand(40)
    e = ECDF(obs)

    # 1.1 is larger than all obs, so ecdf should be 1
    @test ecdf(e, 1.1) ≈ 1.0

    # -1.0 is small than all obs, so ecdf should be 0
    @test ecdf(e, -1.0) ≈ 0.0

    # larger values should have larger values on ecdf
    let x = rand()
        F_1 = ecdf(e, x)
        F_2 = ecdf(e, x*1.1)
        @test F_1 <= F_2
    end

end  # testset


text:
using Flux, DiffEqFlux, DiffEqSensitivity, StochasticDiffEq, Zygote, Test

function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end
function lotka_volterra(u,p,t)
  x, y = u
  α, β, δ, γ = p
  dx = α*x - β*x*y
  dy = -δ*y + γ*x*y
  [dx,dy]
end
function lotka_volterra_noise(du,u,p,t)
  du[1] = 0.01u[1]
  du[2] = 0.01u[2]
end
function lotka_volterra_noise(u,p,t)
  [0.01u[1],0.01u[2]]
end
prob = SDEProblem(lotka_volterra,lotka_volterra_noise,[1.0,1.0],(0.0,10.0))
p = [2.2, 1.0, 2.0, 0.4]
function predict_fd_sde(p)
  solve(prob,SOSRI(),p=p,saveat=0.0:0.1:0.5,sensealg=ForwardDiffSensitivity())[1,:]
end
loss_fd_sde(p) = sum(abs2,x-1 for x in predict_fd_sde(p))
loss_fd_sde(p)

prob = SDEProblem{false}(lotka_volterra,lotka_volterra_noise,[1.0,1.0],(0.0,10.0))
p = [2.2, 1.0, 2.0, 0.4]
function predict_fd_sde(p)
  solve(prob,SOSRI(),p=p,saveat=0.0:0.1:0.5,sensealg=ForwardDiffSensitivity())[1,:]
end
loss_fd_sde(p) = sum(abs2,x-1 for x in predict_fd_sde(p))
loss_fd_sde(p)

@test !iszero(Zygote.gradient(loss_fd_sde,p)[1])

prob = SDEProblem(lotka_volterra,lotka_volterra_noise,[1.0,1.0],(0.0,0.5))
function predict_rd_sde(p)
  solve(prob,SOSRI(),p=p,saveat=0.0:0.1:0.5,sensealg=TrackerAdjoint())[1,:]
end
loss_rd_sde(p) = sum(abs2,x-1 for x in predict_rd_sde(p))
@test !iszero(Zygote.gradient(loss_rd_sde,p)[1])

prob = SDEProblem{false}(lotka_volterra,lotka_volterra_noise,[1.0,1.0],(0.0,0.5))
function predict_rd_sde(p)
  solve(prob,SOSRI(),p=p,saveat=0.0:0.1:0.5,sensealg=TrackerAdjoint())[1,:]
end
loss_rd_sde(p) = sum(abs2,x-1 for x in predict_rd_sde(p))
@test !iszero(Zygote.gradient(loss_rd_sde,p)[1])


text:
function MOI.add_variable(opt::Optimizer)
    add_variable(opt.form)
    return VI(num_variables(opt.form))
end

function MOI.add_variables(opt::Optimizer, N::Int)
    return [MOI.add_variable(opt) for i in 1:N]
end

function MOI.is_valid(opt::Optimizer, v::VI)
    return 1 <= v.value <= num_variables(opt.form)
end

function MOI.get(opt::Optimizer, ::MOI.VariablePrimal, vi::VI)
    MOI.throw_if_not_valid(opt, vi)
    return opt.result.best_solution[vi.value]
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{SV},
    ::Type{<:_V_BOUND_SETS},
)
    return true
end

# TODO: This might "lie" to you: For example, if you set
# both a GT and LT constraint, this will report it as an
# IN constraint.
function _get_scalar_set(form::DMIPFormulation, cvi::CVI)
    l, u = get_bounds(form, cvi)
    if -Inf < l == u < Inf
        return ET
    elseif -Inf < l && u == Inf
        return GT
    elseif l == -Inf && u < Inf
        return LT
    elseif -Inf < l <= u < Inf
        return IN
    else
        @assert l == -Inf && u == Inf
        return nothing
    end
end

function MOI.is_valid(opt::Optimizer, c::CI{SV,S}) where {S<:_V_BOUND_SETS}
    MOI.is_valid(opt, VI(c.value)) || return false
    cvi = CVI(c.value)
    return S == _get_scalar_set(opt.form, cvi)
end

function MOI.get(
    opt::Optimizer,
    ::MOI.ConstraintFunction,
    c::CI{SV,<:_V_BOUND_SETS},
)
    MOI.throw_if_not_valid(opt, c)
    return SV(VI(c.value))
end

# TODO: Do we need to throw if bounds are already set?
function MOI.add_constraint(opt::Optimizer, f::SV, s::ET)
    MOI.throw_if_not_valid(opt, f.variable)
    cvi = convert(CVI, f.variable)
    set_bounds!(opt.form, cvi, IN(s.value, s.value))
    return CI{SV,ET}(index(cvi))
end
function MOI.add_constraint(opt::Optimizer, f::SV, s::GT)
    MOI.throw_if_not_valid(opt, f.variable)
    cvi = convert(CVI, f.variable)
    l, u = get_bounds(opt.form, cvi)
    set_bounds!(opt.form, cvi, IN(s.lower, u))
    return CI{SV,GT}(index(cvi))
end
function MOI.add_constraint(opt::Optimizer, f::SV, s::LT)
    MOI.throw_if_not_valid(opt, f.variable)
    cvi = convert(CVI, f.variable)
    l, u = get_bounds(opt.form, cvi)
    set_bounds!(opt.form, cvi, IN(l, s.upper))
    return CI{SV,LT}(index(cvi))
end
function MOI.add_constraint(opt::Optimizer, f::SV, s::IN)
    MOI.throw_if_not_valid(opt, f.variable)
    cvi = convert(CVI, f.variable)
    l, u = get_bounds(opt.form, cvi)
    set_bounds!(opt.form, cvi, IN(s.lower, s.upper))
    return CI{SV,IN}(index(cvi))
end

MOI.supports_constraint(::Optimizer, ::Type{SV}, ::Type{<:_V_INT_SETS}) = true
function MOI.add_constraint(
    opt::Optimizer,
    f::SV,
    set::S,
) where {S<:_V_INT_SETS}
    vi = f.variable
    MOI.throw_if_not_valid(opt, vi)
    cvi = convert(CVI, vi)
    if get_variable_kind(opt.form, cvi) !== nothing
        error(
            "Already set variable integrality of $(get_variable_kind(opt.form, cvi)) for $(vi); cannot overwrite to $set.",
        )
    end
    set_variable_kind!(opt.form, cvi, set)
    return CI{SV,S}(index(cvi))
end

function MOI.get(opt::Optimizer, ::MOI.NumberOfVariables)
    return num_variables(opt.form)
end

function MOI.get(opt::Optimizer, ::MOI.ListOfVariableIndices)
    return [VI(i) for i in 1:num_variables(opt.form)]
end

function MOI.get(
    opt::Optimizer,
    ::MOI.NumberOfConstraints{SV,S},
) where {S<:_V_BOUND_SETS}
    cnt = 0
    for cvi in all_variables(opt.form)
        if S == _get_scalar_set(opt.form, cvi)
            cnt += 1
        end
    end
    return cnt
end

function MOI.get(
    opt::Optimizer,
    ::MOI.ListOfConstraintIndices{SV,S},
) where {S<:_V_BOUND_SETS}
    indices = CI{SV,S}[]
    for cvi in all_variables(opt.form)
        if S == _get_scalar_set(opt.form, cvi)
            push!(indices, CI{SV,S}(i))
        end
    end
    return indices
end

function MOI.supports(
    ::Optimizer,
    ::MOI.ConstraintPrimal,
    ::CI{SV,<:_V_BOUND_SETS},
)
    return true
end
function MOI.get(
    opt::Optimizer,
    ::MOI.ConstraintPrimal,
    ci::CI{SV,S},
) where {S<:_V_BOUND_SETS}
    vi = VI(ci.value)
    MOI.throw_if_not_valid(opt, vi)
    return MOI.get(opt, MOI.VariablePrimal(), vi)
end


text:
module ShapeModels

using FunctionalData
using FunctionalDataUtils
using HDF5
using LinearAlgebra
using MultivariateStats
using Pkg
using Plots
using Statistics

export PCAShapeModel, shape, coeffs, clamp, meanshape, modeshapes, nmodes, vec, modesstd
export axisij, plotshape, plotshapes
export maxcoeffvec, mincoeffvec

struct PCAShapeModelCoeffs
	modes
	rot
	scale
	translation
    ndims
end

struct PCAShapeModel
	aligned
	pca
    ndims
    nlandmarks
    center
    maxtranslation
    buf::Array{Float64,2}
end

modesstd(a::PCAShapeModel) = sqrt.(principalvars(a.pca))

function PCAShapeModelCoeffs(a::PCAShapeModel, x)
    if a.ndims == 2
        PCAShapeModelCoeffs(x[1:end-4], x[end-3:end-3], x[end-2], x[end-1:end], a.ndims)
    else
        PCAShapeModelCoeffs(x[1:end-6], x[end-5:end-4], x[end-3], x[end-2:end], a.ndims)
    end
end
import Base.vec
vec(a::PCAShapeModelCoeffs) = [a.modes; a.rot; a.scale; a.translation]

eye(a) = diagm(0 => ones(a))
function PCAShapeModel(landmarks::Array{T,3};
                       percentage = 0.98,
                       center = zeros(size(landmarks,1)),
                       maxtranslation = typemax(Float32)*ones(size(landmarks,1)),
                       optcoeffs = false) where T<:Real

    ndims, nlandmarks, nshapes = size(landmarks)

    aligned = copy(landmarks)
    scales = zeros(nshapes)
    aligned = @p map aligned x->x .- mean(x, dims = 2)
    scales = @p map aligned x->@p distance zeros(sizem(aligned),1) x | mean
    relativescales = scales ./ mean(scales)

    rotmatrices = Array{Any}(undef, nshapes)
    rotmatrices[1] = eye(ndims)
    aligned[:,:,1] = aligned[:,:,1] ./= relativescales[1]
    for i = 2:nshapes
        aligned[:,:,i] ./= relativescales[i]
        # https://en.wikipedia.org/wiki/Kabsch_algorithm
        U,E,V = svd(aligned[:,:,1] * aligned[:,:,i]')
        Eprime = eye(length(E))
        Eprime[end] = sign(det(U*V))
        rotmatrices[i] = V' * Eprime' * U
        aligned[:,:,i] = rotmatrices[i] * aligned[:,:,i]
    end

	alignedshapes = aligned
    aligned = reshape(aligned, (ndims*nlandmarks, nshapes))
	pca = fit(PCA, aligned, pratio = percentage)

    model = PCAShapeModel(alignedshapes, pca, ndims, nlandmarks, center, maxtranslation, zeros(ndims, nlandmarks))

    ## compute optimal coeffs for training data
    if optcoeffs
        invrot = @p map rotmatrices inv
        angles = @p map invrot x->asin(x[2,1])
        modes  = @p transform pca aligned | unstack
        translations = @p map landmarks (x->mean(x, dims = 2)) | map minus center | unstack
        coeffs = @p zip modes angles relativescales-1 translations| collect | map x->PCAShapeModelCoeffs(x[1], x[2], x[3], x[4], ndims)
        model, coeffs
    else
        model
    end
end

import Base.reshape
col(a) = reshape(a, length(a),1)
row(a) = reshape(a, 1, length(a))
reshape(a::PCAShapeModel, b) = reshape(b, a.ndims, a.nlandmarks)
meanshape(a::PCAShapeModel) = shape(a, zeros(nmodes(a)))

maxcoeffvec(a::PCAShapeModel) = vcat(
    2.5*modesstd(a),
    a.ndims == 2 ? 0.3 : [0.3, 0.3, 0.3],
    0.2,
    a.maxtranslation)
mincoeffvec(a::PCAShapeModel) = -maxcoeffvec(a)

nmodes(a::PCAShapeModel) = outdim(a.pca) + (a.ndims==2 ? 4 : 6)

function modeshapes(a::PCAShapeModel, ind, at = range(mincoeffvec(a)[ind], stop = maxcoeffvec(a)[ind], length = 10))
    @assert ind>0 && ind <= nmodes(a)
	r = zeros(a.ndims, a.nlandmarks, length(at))
	for i = 1:length(at)
    	v = zeros(nmodes(a))
		v[ind] = at[i]
		r[:,:,i] = shape(a, v)
	end
	r
end

rotmatrix(coeffs::PCAShapeModelCoeffs) = rotmatrix(coeffs.scale, coeffs.rot...)
rotmatrix(scale, alpha) = (scale+1)*[cos(alpha) -sin(alpha); sin(alpha) cos(alpha)]
function rotmatrix(scale, rotm, rotn, roto)
    Rm = [1 0 0; 0 cos(rotm) -sin(rotm); 0 sin(rotm) cos(rotm)]
    Rn = [cos(rotn) 0 sin(rotn); 0 1 0; -sin(rotn) 0 cos(rotn)]
    Ro = [cos(roto) -sin(roto) 0; sin(roto) cos(roto) 0; 0 0 1]
    (scale+1)*Rm*Rn*Ro
end

function shape(a::PCAShapeModel, coeffs)
    if sizem(coeffs) != nmodes(a)
        error("Cant create shape from coeffs, dimension mismatch. size(coeffs) == $(size(coeffs)), should be $(nmodes(a))")
    end
    shape(a, PCAShapeModelCoeffs(a, coeffs))
end
shape(a::PCAShapeModel, coeffs::PCAShapeModelCoeffs) = shape!(a.buf, a, coeffs)
function shape!(buf, a::PCAShapeModel, coeffs::PCAShapeModelCoeffs)
    r = reconstruct(a.pca, coeffs.modes)
    r = reshape(a, r)
    buf[:] = rotmatrix(coeffs)*r .+ coeffs.translation .+ a.center
end

(a::PCAShapeModel)(coeffs) = shape(a, coeffs)

function coeffs(a::PCAShapeModel, coords::Array{T,2}) where T<:Real
    # TODO
end

import Base.clamp
clamp(a::PCAShapeModel, coeffs::Vector) = clamp(a, col(coeffs))
function clamp(a::PCAShapeModel, coeffs::Array{T,2}) where T<:Real
    # TODO
    error("write me")
end

examplelandmarks() = examplelandmarks(:hands2d)
function examplelandmarks(a::Symbol)
    if a==:hands2d
        return h5read(joinpath(dirname(@__FILE__()),"../data/2Dlandmarks.hdf5"),"landmarks")
    elseif a==:lungs
        return dictread(joinpath(dirname(@__FILE__()),"../data/lungs.dictfile"))[:landmarks]
    else
        error("unknown dataset '$a'")
    end
end

function exampleimages(a::Symbol)
    if a == :hands2d
        shapes = examplelandmarks(a)
        return [inpolygon(at(shapes,i).-mean(at(shapes,i)).+[200 250]',zeros(500,400)) for i in 1:14]
    else
        error("unknown dataset '$a'")
    end
end

function axisij()
    a = axis()
    axis((a[1],a[2],a[4],a[3]))
end

function plotshape(a, args...; kargs...)
    if size(a,1)==2
        plot(vec(a[2,:]), vec(a[1,:]), args...; kargs...)
        # axisij()
        # xlabel("N")
        # ylabel("M")
    elseif size(a,1)==3
        plot3d(vec(a[1,:]), vec(a[2,:]), vec(a[3,:]), args...; kargs...)
        # xlabel("M")
        # ylabel("N")
        # zlabel("O")
    else
        error("Can't plot data of size $(size(a))")
    end
    # axis("equal")
end

function plotshapes(a)
    if size(a,1)==3
        for i = 1:last(size(a))
            plotshape(a[:,:,i])
            hold(true)
        end
        hold(false)
    else
        gridplot(a)
    end
end

at(a,i) = slicedim(a,ndims(a),i)

function gridplot(a)
    N = last(size(a))
    sm = floor(sqrt(N))
    sn = ceil(N/sm)
    for m = 1:sm, n = 1:sn
        datai = (n-1)*sm+m
        ploti = (m-1)*sn+n
        if datai <= N
            subplot(sm,sn,ploti)
            plotshape(at(a,round(Int,datai)))
        end
    end
end

end # module


text:
type Garray
    ahandle::Array{Ptr{Void}}
    atyp::DataType
    elem_size::Int64
    access_iob::IOBuffer
    access_arr::Array
end

const GarrayMemoryHandle = IOBuffer

function Garray(T::DataType, elem_size::Int64, num_elems::Int64)
    a = Garray([C_NULL], T, elem_size, IOBuffer(), [])
    r = ccall((:garray_create, libgasp), Cint, (Ptr{Void}, Int64, Int64,
              Ptr{Int64}, Ptr{Void}), ghandle[1], num_elems, a.elem_size,
              C_NULL, pointer(a.ahandle, 1))
    if r != 0
        error("construction failure")
    end
    global num_garrays
    num_garrays = num_garrays+1
    finalizer(a, (function(a)
                    ccall((:garray_destroy, libgasp),
                            Void, (Ptr{Void},), a.ahandle[1])
                    global num_garrays
                    num_garrays = num_garrays-1
                    exiting && num_garrays == 0 && __shutdown__()
                  end))
    return a
end

function length(ga::Garray)
    ccall((:garray_length, libgasp), Int64, (Ptr{Void},), ga.ahandle[1])
end

function elemsize(ga::Garray)
    ccall((:garray_elemsize, libgasp), Int64, (Ptr{Void},), ga.ahandle[1])
end

function get(ga::Garray, lo::Int64, hi::Int64)
    adjlo = lo - 1
    adjhi = hi - 1
    getlen = hi - lo + 1
    cbuflen = getlen * ga.elem_size
    cbuf = Array{UInt8}(cbuflen)
    r = ccall((:garray_get, libgasp), Cint, (Ptr{Void}, Int64, Int64,
              Ptr{Void}), ga.ahandle[1], adjlo, adjhi, cbuf)
    if r != 0
        error("Garray get failed")
    end
    iob = IOBuffer(cbuf)
    buf = Array{ga.atyp}(getlen)
    for i = 1:length(buf)
        try
            buf[i] = deserialize(iob)
        catch exc
            if !isa(exc, UndefRefError) && !isa(exc, BoundsError)
                rethrow()
            end
            # this is expected when an array element is uninitialized
        end
        seek(iob, i * ga.elem_size)
    end
    return buf, iob
end

function put!(ga::Garray, lo::Int64, hi::Int64, buf::Array)
    adjlo = lo - 1
    adjhi = hi - 1
    putlen = hi - lo + 1
    cbuflen = putlen * ga.elem_size
    cbuf = Array{UInt8}(cbuflen)
    iob = IOBuffer(cbuf, true, true)
    for i = 1:length(buf)
        serialize(iob, buf[i])
        seek(iob, i * ga.elem_size)
    end
    r = ccall((:garray_put, libgasp), Cint, (Ptr{Void}, Int64, Int64,
              Ptr{Void}), ga.ahandle[1], adjlo, adjhi, cbuf)
    if r != 0
        error("Garray put failed")
    end
end

function distribution(ga::Garray, rank::Int64)
    lo = Ref{Int64}(0)
    hi = Ref{Int64}(0)
    r = ccall((:garray_distribution, libgasp), Cint, (Ptr{Void}, Int64,
            Ptr{Int64}, Ptr{Int64}), ga.ahandle[1], rank-1, lo, hi)
    if r != 0
        error("could not get distribution")
    end
    llo = lo[] + 1
    lhi = hi[] + 1
    return llo, lhi
end

function access(ga::Garray, lo::Int64, hi::Int64)
    p = [C_NULL]
    r = ccall((:garray_access, libgasp), Cint, (Ptr{Void}, Int64, Int64,
              Ptr{Ptr{Void}}), ga.ahandle[1], lo-1, hi-1, pointer(p, 1))
    if r != 0
        error("could not get access")
    end
    acclen = hi - lo + 1
    buf = Array{ga.atyp}(acclen)
    if length(buf) == 0
        return buf
    end
    cbuflen = acclen * ga.elem_size
    iob = IOBuffer(unsafe_wrap(Array, convert(Ptr{UInt8}, p[1]), cbuflen),
                   true, true)

    for i = 1:length(buf)
        try
            buf[i] = deserialize(iob)
        catch exc
            if !isa(exc, UndefRefError) && !isa(exc, BoundsError)
                rethrow()
            end
            # this is expected when an array element is uninitialized
        end
        seek(iob, i * ga.elem_size)
    end
    seek(iob, 0)
    ga.access_iob = iob
    ga.access_arr = buf
    return buf
end

function flush(ga::Garray)
    if ga.access_arr != []
        for i = 1:length(ga.access_arr)
            try serialize(ga.access_iob, ga.access_arr[i])
            catch exc
                if !isa(exc, UndefRefError)
                    rethrow()
                end
            end
            seek(ga.access_iob, i * ga.elem_size)
        end
        ga.access_arr = []
    end
    ccall((:garray_flush, libgasp), Void, (Ptr{Void},), ga.ahandle[1])
end


text:
using Documenter, AbstractTrees

makedocs(
    modules = [AbstractTrees],
    sitename = "AbstractTrees.jl",
    authors = "Keno Fischer",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages = [
        "Home" => "index.md",
        "Iteration" => "iteration.md",
        "Internals" => "internals.md",
        "FAQ" => "faq.md",
    ],
)

deploydocs(
    repo = "github.com/JuliaCollections/AbstractTrees.jl.git",
)


text:
################################################################################
# Roadway
################################################################################

abstract type RoadwayOptions
end

mutable struct Roadway{T}
    l::Int
    lane::Vector{Lane{T}}
	opts::RoadwayOptions
end

function Roadway(lanes::Vector{Lane{T}}, opts::RoadwayOptions) where {T}
    l = length(lanes)
    sort!(lanes, by = x -> x.id)
    if [lanes[i].id for i=1:l] != (1:l)
        # @show "Reindexing the lanes"
        for i = 1:l
            lanes[i].id = i
        end
    end
    return Roadway{T}(l,lanes,opts)
end

function add_lane!(roadway::Roadway{T}, lane::Lane{T}) where {T}
    roadway.l += 1
    lane.id = roadway.l
    push!(roadway.lane, lane)
    return nothing
end

function test_centerline_continuity(f; v::T=1.0, s_min::T=-10.0, s_max::T=10.0,
	sample::Int=2000) where T
	# plt = scatter(legend=false)
	out = true
	step = (s_max - s_min)/sample
	vs_ = f(s_min,v)
	vs = f(s_min,v)
	for k = 1:sample
		s = s_min + k*step
		vs = f(s,v)
		out &= norm([vs.x, vs.y] - [vs_.x, vs_.y]) <= (step + 1e-5)
		vs_ = vs
		# scatter!([vs.x], [vs.y])
	end
	# display(plt)
	return out
end


text:
module MixtureModel_2D

using SPPL

n = @sppl begin
    X ~ SPPL.Normal(0, 2)
    Y ~ 0.6 * SPPL.Normal(8, 1) | 0.4 * SPPL.Gamma(-3, 3)
end

e0 = ((-4 < n.X) < 4) & ((1 < n.Y^2) < 4)
e1 = ((-1 < n.X) < 1) & ((-1.5 < n.Y) < 6)
modelc = condition(model, e0 | e1)
e = dnf_to_disjoint_union(e0 | e1)
println([n.model.prob(x) for x in e.subexprs])

end # module


text:
function vertex_input_attribute_descriptions(::Type{T}, binding, formats=Format.(fieldtypes(T))) where {T}
    VertexInputAttributeDescription.(
        0:fieldcount(T)-1,
        binding,
        formats,
        fieldoffset.(T, 1:fieldcount(T)),
    )
end

Vulkan.VertexInputBindingDescription(::Type{T}, binding; input_rate = VERTEX_INPUT_RATE_VERTEX) where {T} =
    VertexInputBindingDescription(binding, sizeof(T), input_rate)


text:
using .JSort

using .JSort
import Base: read, read!, fill!, write, close

mutable struct QuickReader <: AbstractReader
    parameters::Parameters
    input::CBuffer
    labr_size::Int
    initialtime::Int
    eγ::Vector{Vector{Float32}}
    source::String
end

function QuickReader(parameters::Parameters)
    j(x) = joinpath(parameters.savepath, x)
    source = joinpath(parameters.readpath, "raw.bin")
    if ~isfile(source)
        error("The file $source does not exist")
    end
    buffer = CBuffer(source)
    egs = [Float32[] for i in 1:33]

    QuickReader(parameters, buffer, 0, -1, egs, source)
end

read(reader::QuickReader, gate::Gate) = read(reader, gate=gate)
read(reader::QuickReader, gate::Gate2D) = read(reader, pgate=gate)

function read(reader::QuickReader; gate::Gate=Gate(), pgate=nothing)
    close(reader.input)
    reader.input = CBuffer(reader.source)

    qkinzpoly = reader.parameters.qkinz
    timecalib = reader.parameters.calibrator[:t]
    timecorrection = reader.parameters.calibrator[:tc]
    gammacalib = reader.parameters.calibrator[:γ]
    ecalib = reader.parameters.calibrator[:e]
    Δecalib = reader.parameters.calibrator[:Δe]
    prompt = reader.parameters.gates[:prompt]
    background = reader.parameters.gates[:background]
    event = Event()
    read_events = front = back = 0
    # TODO Same type as calibrators
    # TODO Buffer2D are not overwritten
    e = Δe = eΔe = Float32(0.0)
    eₓᵗ = eᵧ = Float32(0.0)
    t = tc = Float32(0.0)
    reader.eγ = [Float32[] for i in 1:33]
    for eγ in reader.eγ
        sizehint!(eγ, 2^13)
    end

    while isgood(reader)
        read_events += 1
        fill!(event, reader)
        front, back = event.front, event.back
        e = ecalib[front, back](event.e)
        Δe = Δecalib[front, back](event.de)
        if !isnothing(pgate)
            !pgate(e, Δe) && continue
        end
        # Theoretical excitation energy from kinematics
        eΔe = e + Δe
        # This is a `constant` poly. No calibration needed
        eₓᵗ = qkinzpoly[front](eΔe)
        !gate(eₓᵗ) && continue

        @views for labr in event.labr[1:reader.labr_size]
            labr.channel ∈ reader.parameters.badchannels && continue
            # Gamma detectors
            eᵧ = gammacalib[labr.channel](labr.adc)

            # Time
            if reader.initialtime == -1
                reader.initialtime = event.labr[1].tdc
            end
            t = timecalib[labr.channel](labr.tdc)
            !prompt(t) && continue
            #tc = timecorrection[labr.channel](t, eΔe)
            # Also correct using eᵧ?
            push!(reader.eγ[labr.channel], eᵧ)
        end
    end
end

function close(reader::QuickReader)
    close(reader.input)
end

function ChannelDetector(reader::QuickReader)
    data = GammaDetector[GammaDetector(reader.eγ[chn], chn) for chn in 1:33]
    ChannelDetector(data)
end

plot(r::QuickReader; kwargs...) = plot(ChannelDetector(r); kwargs...)

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


max_stars_repo_name: saurabhcommand/Hello-world

max_issues_repo_name: saurabhcommand/Hello-world

max_forks_repo_name: saurabhcommand/Hello-world


text:
## This file implements Example H071 using the MathProgBase Interface.
## This is the example code from the MathProgBase test file from https://github.com/JuliaOpt/MathProgBase.jl/blob/master/test/nlp.jl , slightly modified to remove the tests.
## This interface is solver independent, but will use Ipopt as the solver here.

using MathProgBase
import MathProgBase.MathProgSolverInterface   ##In Julia, importing means that the prefix 'MathProgSolverInterface.' must be included before calling any functions within this module
using Ipopt

tic()  #Start the clock

# Here the type represents the complete instance, but it
# could also store instance data.
# The abstract type AbstractNLPEvaluator is used by solvers for accessing the objective function f and constraints g.
# Solvers may query the value, gradients, Hessian-vector products, and the Hessian of the Lagrangian.
type HS071 <: MathProgSolverInterface.AbstractNLPEvaluator
end

# hs071
# min x1 * x4 * (x1 + x2 + x3) + x3
# st  x1 * x2 * x3 * x4 >= 25
#     x1^2 + x2^2 + x3^2 + x4^2 = 40
#     1 <= x1, x2, x3, x4 <= 5
# Start at (1,5,5,1)
# End at (1.000..., 4.743..., 3.821..., 1.379...)

## Initialize must be called before any other methods. The vector requested_features lists features requested by the solver.
function MathProgSolverInterface.initialize(d::HS071, requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in [:Grad, :Jac, :Hess])
            error("Unsupported feature $feat")
        end
    end
end

## Returns which features are available for the model HS071
MathProgSolverInterface.features_available(d::HS071) = [:Grad, :Jac, :Hess]

################################
# Define the Objective Function
################################
MathProgSolverInterface.eval_f(d::HS071, x) = x[1] * x[4] * (x[1] + x[2] + x[3]) + x[3]

###################
# Constraint set
##################
function MathProgSolverInterface.eval_g(d::HS071, g, x)
    g[1] = x[1]   * x[2]   * x[3]   * x[4]
    g[2] = x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2
end

#####################################
# Gradient of the Objective Function
#####################################
function MathProgSolverInterface.eval_grad_f(d::HS071, grad_f, x)
    grad_f[1] = x[1] * x[4] + x[4] * (x[1] + x[2] + x[3])
    grad_f[2] = x[1] * x[4]
    grad_f[3] = x[1] * x[4] + 1
    grad_f[4] = x[1] * (x[1] + x[2] + x[3])
end

###################################################################################################################################
##These functions collect the sparsity structures of the Jacobian matrix and the Hessian of the Lagrangian (row first then columns)
MathProgSolverInterface.jac_structure(d::HS071) = [1,1,1,1,2,2,2,2],[1,2,3,4,1,2,3,4]
# lower triangle only
MathProgSolverInterface.hesslag_structure(d::HS071) = [1,2,2,3,3,3,4,4,4,4],[1,1,2,1,2,3,1,2,3,4]

#########################################################################################################################
##OR you could define the values of the non-zero elements. The function below does this for the Jacobian matrix and saves
## them to the matrix J. There is no need to create the matrix J.
####################################################################
function MathProgSolverInterface.eval_jac_g(d::HS071, J, x)
    # Constraint (row) 1
    J[1] = x[2]*x[3]*x[4]  # 1,1
    J[2] = x[1]*x[3]*x[4]  # 1,2
    J[3] = x[1]*x[2]*x[4]  # 1,3
    J[4] = x[1]*x[2]*x[3]  # 1,4
    # Constraint (row) 2
    J[5] = 2*x[1]  # 2,1
    J[6] = 2*x[2]  # 2,2
    J[7] = 2*x[3]  # 2,3
    J[8] = 2*x[4]  # 2,4
end

####################################################################################################
# Similarly for the Hessian -  σ = weight on objective function, μ = multipliers on the constraints
####################################################################################################
function MathProgSolverInterface.eval_hesslag(d::HS071, H, x, σ, μ)
    # Again, only lower left triangle (since it is symmetric)
    # Objective
    H[1] = σ * (2*x[4])               # 1,1
    H[2] = σ * (  x[4])               # 2,1
    H[3] = 0                          # 2,2
    H[4] = σ * (  x[4])               # 3,1
    H[5] = 0                          # 3,2
    H[6] = 0                          # 3,3
    H[7] = σ* (2*x[1] + x[2] + x[3])  # 4,1
    H[8] = σ * (  x[1])               # 4,2
    H[9] = σ * (  x[1])               # 4,3
    H[10] = 0                         # 4,4

    # First constraint
    H[2] += μ[1] * (x[3] * x[4])  # 2,1
    H[4] += μ[1] * (x[2] * x[4])  # 3,1
    H[5] += μ[1] * (x[1] * x[4])  # 3,2
    H[7] += μ[1] * (x[2] * x[3])  # 4,1
    H[8] += μ[1] * (x[1] * x[3])  # 4,2
    H[9] += μ[1] * (x[1] * x[2])  # 4,3

    # Second constraint
    H[1]  += μ[2] * 2  # 1,1
    H[3]  += μ[2] * 2  # 2,2
    H[6]  += μ[2] * 2  # 3,3
    H[10] += μ[2] * 2  # 4,4

end

# Now put this all together to define and solve the model.
# First, define a model object and the solver to use
    m = MathProgSolverInterface.model(MathProgBase.defaultNLPsolver )
# The vector of lower bounds on the variables
    l = [1,1,1,1]
# The vector of upper bounds on the variables
    u = [5,5,5,5]
# Vector of lower bounds on the constraints
    lb = [25, 40]
# Vector of constraint upper bounds
    ub = [Inf, 40]
# Load the nonlinear programming problem into the model - syntax for loadnonlinearproblem!() is
# loadnonlinearproblem!(m::AbstractMathProgModel, numVar, numConstr, l, u, lb, ub, sense, d::AbstractNLPEvaluator)
    MathProgSolverInterface.loadnonlinearproblem!(m, 4, 2, l, u, lb, ub, :Min, HS071())
## Set the initial values for the optimisation
    MathProgSolverInterface.setwarmstart!(m,[1,5,5,1])

## Solve the optimisation problem
    MathProgSolverInterface.optimize!(m)
## Let's see the termination status after solving
    stat = MathProgSolverInterface.status(m)

# Now let's see the solution.
    x = MathProgSolverInterface.getsolution(m)

toc()

println(x)

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
if !isdefined(Main, :FCFP)
    include(joinpath(@__DIR__, "FCFP.jl"))
end

"""
The ``LocalOptimizer`` module define methods to optimize locally a Fixed
Cost Flow Problem (FCFP) from an intial solution. It consists mainly by
rerouting main flow arteries. A flow artery is defined as a sequence of
edges where the flow is constant (for now strictly more than one edge).
The optimizer tries to find a path replacing the artery but with a smallest
cost. To do that it computes the shortest path between some nodes that are
already visited by the flow (e.g. parent nodes of the artery) and the nodes
that the artery is sustaining. If such a path is found, the artery will be
removed from the flow graph and replaced by the found path.

For now, the hypothesis of the optimizer are the following:
    * The flow graph is a set of DAG (directed acyclic graphs)
    * Edges have no capacity
    * The initial solution is feasible

"""
module LocalOptimizer

    #include("DSSP.jl")
    include(joinpath(@__DIR__ , "graph_utils.jl"))
    import Main.FCFP.FCFPInstance
    using DataStructures

    # TODO: catch errors (is not a DAG etc), and except with returning the initial flows

    mutable struct LocalOptimizerInstance
        init_flows::FlowDictType
        working_flows::FlowDictType
        flows::FlowDictType
        optimization_graph::MultiDiGraph
        fcfp::FCFPInstance
        costs::Dict{EdgeType, WeightType}
        roots::Array{VertexType, 1}
        treated_arteries::Set
        logger
    end

    function _init_loi(flows::Dict, fcfp::FCFPInstance, logger)
        optimization_graph = MultiDiGraph(keys(flows))
        i_flows = FlowDictType(flows)
        flow_dag = get_flow_dag(optimization_graph, i_flows)
        roots = [n for n in get_nodes(flow_dag) if has_no_parent(flow_dag, n)]
        treated_arteries = Set()
        costs = Dict{EdgeType, WeightType}()
        sizehint!(costs, length(fcfp.fixed_costs))
        for (u, v, _) in keys(fcfp.fixed_costs)
            costs[(u, v)] = minimum([fcfp.fixed_costs[(u, v, k)] for k=optimization_graph.succ[u][v]])
        end
        return LocalOptimizerInstance(i_flows, FlowDictType(flows), FlowDictType(flows),
                                        optimization_graph, fcfp, costs,
                                        roots, treated_arteries, logger)
    end

    LocalOptimizerInstance(flows::Dict, fcfp::FCFPInstance, logger) = _init_loi(flows, fcfp, logger)

    function validate_move(loi::LocalOptimizerInstance)
        loi.flows = copy(loi.working_flows)
    end

    function optimize_locally(flows::Dict, fcfp::FCFPInstance, logger=println) :: Dict
        start = time()
        loi = LocalOptimizerInstance(flows, fcfp, logger)
        initial_cost = sum(loi.fcfp.fixed_costs[e] for (e, f) in loi.flows if f > 0)
        clean_flow_two_cycles(loi.working_flows)
        break_cycles(loi, loi.working_flows)
        validate_move(loi)
        while replace_arteries(loi)
            gain = initial_cost - sum(loi.fcfp.fixed_costs[e] for (e, f) in loi.flows if f > 0)
            loi.logger("Cost gain: $(gain)")
        end
        gain = initial_cost - sum(loi.fcfp.fixed_costs[e] for (e, f) in loi.flows if f > 0)
        optimization_time = time() - start
        loi.logger("Arteries flows rerouting postprocess: gain=$(gain),  time: $(optimization_time)")
        if check_solution(loi, loi.flows)
            return loi.flows
        else
            loi.logger("Local optimization failed")
            return loi.init_flows
        end
    end

    function replace_arteries(loi::LocalOptimizerInstance) :: Bool
        flow_dag = get_flow_dag(loi.optimization_graph, loi.flows)
        arteries = get_flow_arteries(loi, flow_dag)
        roots_descendants = Array{Set{VertexType}, 1}([descendants(flow_dag, r, false) for r in loi.roots])
        roots_marginal_capacity = [loi.fcfp.production[r] - sum(get(loi.fcfp.heat_demand, n, 0.0)
                                                                for n in roots_descendants[r_index])
                                   for (r_index, r) in enumerate(loi.roots)]
        while !isempty(arteries)
            a = dequeue!(arteries)
            # determine sub tree
            end_node = last(a)[2]
            sub_tree_nodes = descendants(flow_dag, end_node, false)
            if isempty(sub_tree_nodes)
                continue
            end
            sub_tree_demand = sum(get(loi.fcfp.heat_demand, n, 0.0) for n in sub_tree_nodes)
            # determine parent sub tree
            start_node = first(a)[1]
            current_roots_indexes = Set{Int64}()
            other_roots_indexes = Set{Int64}()
            for (d_index, d) in enumerate(roots_descendants)
                if (start_node in d)
                    push!(current_roots_indexes, d_index)
                elseif (roots_marginal_capacity[d_index] >= sub_tree_demand)
                    push!(other_roots_indexes, d_index)
                end
            end
            concerned_sub_trees_indexes = union(current_roots_indexes, other_roots_indexes)
            parent_dag_trees_nodes = cascade_union([roots_descendants[d_index] for d_index in concerned_sub_trees_indexes])
            parent_dag_trees_nodes = setdiff(parent_dag_trees_nodes, union(sub_tree_nodes,
                                                                           Set(e[2] for e in a)))
            if isempty(parent_dag_trees_nodes)
                continue
            end
            # Add 0-weight edges in the sub-graph
            modified_costs = add_zero_weight_edges_from_nodes(sub_tree_nodes, loi.costs)
            path, cost = dijkstra_multisource(loi.optimization_graph,
                                              parent_dag_trees_nodes,
                                              end_node,
                                              modified_costs)
            edge_like_path = get_edge_like_path(path, loi.optimization_graph, loi.fcfp.fixed_costs)
            if edge_like_path != a
                potential_cycle_creation = false
                if any((n != end_node) && (n in sub_tree_nodes) for n in path)
                    potential_cycle_creation = true
                end
                new_start = first(path)
                if any(new_start in roots_descendants[d_index] for d_index in current_roots_indexes)
                    loi.logger("Moving flow inside current dag")
                    reroute_artery(loi, a, edge_like_path, flow_dag, potential_cycle_creation)
                else
                    loi.logger("Moving flow to another dag")
                    new_roots_candidates = Set{VertexType}([loi.roots[r] for r in other_roots_indexes if new_start in roots_descendants[r]])
                    current_roots = Set{VertexType}([loi.roots[r] for r in current_roots_indexes])
                    move_artery_to_other_dag(loi, a, edge_like_path, flow_dag, current_roots, new_roots_candidates, potential_cycle_creation)
                end
                validate_move(loi)
                return true
            else
                push!(loi.treated_arteries, a)
            end
        end
        return false
    end

    function reroute_artery(loi::LocalOptimizerInstance, a::Array{MultiEdgeType, 1},
                            new_a::Array{MultiEdgeType, 1}, flow_dag::MultiDiGraph, potential_cycle_creation=false)
        artery_flow = loi.working_flows[first(a)]
        for e in a
            loi.working_flows[e] = zero(WeightType)
        end
        new_start = first(new_a)[1]
        start_node = first(a)[1]
        if new_start != start_node
            lca = lowest_common_ancestor(flow_dag, start_node, new_start)
            if lca === nothing
                throw(ErrorException("Unable to find common ancestor"))
            end
            # 1: Browse back the flow to not show flow in this direction
            if lca != start_node
                back_path, _ = dijkstra_multisource(flow_dag, Set{VertexType}([lca]), start_node, loi.costs)
                edge_back_path = get_edge_like_path(back_path, flow_dag, loi.fcfp.fixed_costs)
                for e in edge_back_path
                    if artery_flow <= loi.working_flows[e]
                        loi.working_flows[e] -= artery_flow
                    else
                        potential_cycle_creation = true
                        loi.working_flows[(e[2], e[1], e[3])] = artery_flow - loi.working_flows[e]
                        loi.working_flows[e] = zero(WeightType)
                    end
                end
            end
            # 2: Browse down the flow to the start of the new artery
            if lca != new_start
                (down_path, _) = dijkstra_multisource(flow_dag, Set{VertexType}([lca]), new_start, loi.costs)
                edge_down_path = get_edge_like_path(down_path, flow_dag, loi.fcfp.fixed_costs)
                for e in edge_down_path
                    loi.working_flows[e] += artery_flow
                end
            end
        end
        for e in new_a
            loi.working_flows[e] += artery_flow  # Adding flow because their may be already a flow on it
        end
        clean_flow_two_cycles(loi.working_flows)
        if potential_cycle_creation
            break_cycles(loi, loi.working_flows, Set{VertexType}([new_start]))
        end
    end

    function move_artery_to_other_dag(loi::LocalOptimizerInstance,
                                      a::Array{MultiEdgeType, 1},
                                      new_a::Array{MultiEdgeType, 1},
                                      flow_dag::MultiDiGraph,
                                      current_roots:: Set{VertexType},
                                      new_roots_candidates::Set{VertexType},
                                      potential_cycle_creation=false)
        artery_flow = loi.working_flows[first(a)]
        for e in a
            loi.working_flows[e] = zero(WeightType)
        end
        # Find path from old root to old start and remove flow
        old_start = first(a)[1]
        residual_artery_flow = artery_flow
        for r in current_roots  # Removing artery flow taking into account multiple roots
            root_flow = sum(loi.working_flows[e] for e in out_edges(flow_dag, r))
            flow_to_remove = min(root_flow, residual_artery_flow)
            (back_path, _) = dijkstra_multisource(flow_dag, Set{VertexType}([r]), old_start, loi.costs)
            edge_back_path = get_edge_like_path(back_path, flow_dag, loi.fcfp.fixed_costs)
            for e in edge_back_path
                loi.working_flows[e] -= flow_to_remove
            end
            residual_artery_flow -= flow_to_remove
            if residual_artery_flow <= 0.0
                break
            end
        end
        # Find path from new root to new start and add flow
        new_start = first(new_a)[1]
        (down_path, _) = dijkstra_multisource(flow_dag, new_roots_candidates, new_start, loi.costs)
        edge_down_path = get_edge_like_path(down_path, flow_dag, loi.fcfp.fixed_costs)
        for e in edge_down_path
            loi.working_flows[e] += artery_flow
        end
        # Set new artery flow
        for e in new_a
            loi.working_flows[e] += artery_flow  # Adding flow because their may be already a flow on it
        end
        # Clean flows
        clean_flow_two_cycles(loi.working_flows)
        if potential_cycle_creation
            break_cycles(loi, loi.working_flows, Set{VertexType}([new_start]))
        end
    end

    function get_edge_like_path(path::Array, graph::MultiDiGraph, costs::Dict) :: Array{MultiEdgeType, 1}
        edge_like_path = Array{MultiEdgeType, 1}()
        for i in 1:(length(path) - 1)
            u, v = path[i], path[i + 1]
            min_edge = (u, v, first(graph.succ[u][v]))
            min_weight = costs[min_edge]
            for k in graph.succ[u][v]
                if costs[(u, v, k)] < min_weight
                    min_edge = (u, v, k)
                    min_weight = costs[min_edge]
                end
            end
            push!(edge_like_path, min_edge)
        end
        return edge_like_path
    end

    function add_zero_weight_edges_from_nodes(nodes::Set{VertexType}, costs::Dict) :: Dict
        modified_costs = copy(costs)
        for (u, v) in keys(costs)
            if (u in nodes) && (v in nodes)
                modified_costs[(u, v)] = zero(WeightType)
                modified_costs[(v, u)] = zero(WeightType)
            end
        end
        return modified_costs
    end

    function remove_zero_weight_edges(graph, added_edges, costs)
        remove_edges!(graph, added_edges)
        for e in added_edges
            delete!(costs, e)
        end
    end

    function cascade_union(sets::Array{Set{T}, 1}) :: Set{T} where {T<:Any}
        u = Set{T}()
        for s in sets
            u = union(u, s)
        end
        return u
    end

    function get_flow_arteries(loi::LocalOptimizerInstance, flow_dag::MultiDiGraph) :: PriorityQueue{Array{MultiEdgeType, 1}, WeightType}
        arteries = PriorityQueue{Array{MultiEdgeType, 1}, WeightType}()
        explored_nodes = Set{VertexType}()
        for r in loi.roots
            nodes_stack = Array{VertexType, 1}([r])
            while !isempty(nodes_stack)
                n = pop!(nodes_stack)
                if n in explored_nodes
                    continue
                end
                push!(explored_nodes, n)
                out_flow_edges = collect(out_edges(flow_dag, n))
                if (length(out_flow_edges) != 1) || (! has_one_parent(flow_dag, n))
                    for e in out_flow_edges
                        push!(nodes_stack, e[2])
                    end
                    continue
                end
                # A new artery is found, explore it in depth
                current_artery = [first(in_edges(flow_dag, n)), out_flow_edges[1]]
                sub_nodes_stack = [out_flow_edges[1][2]]
                while !isempty(sub_nodes_stack)
                    sub_node = pop!(sub_nodes_stack)
                    push!(explored_nodes, sub_node)
                    out_flow_edges = collect(out_edges(flow_dag, sub_node))
                    if (length(out_flow_edges) != 1) || (! has_one_parent(flow_dag, sub_node))
                        # End of the artery
                        artery_cost = - sum(loi.fcfp.fixed_costs[e] for e in current_artery)
                        if current_artery in loi.treated_arteries
                            artery_cost += 1e5
                        end
                        arteries[current_artery] = artery_cost
                        for e in out_flow_edges
                            push!(nodes_stack, e[2])
                        end
                        break
                    else
                        push!(current_artery, out_flow_edges[1])
                        push!(sub_nodes_stack, out_flow_edges[1][2])
                    end
                end
            end
        end
        return arteries
    end

    function sort_arteries(loi::LocalOptimizerInstance, arteries) :: Array
        failed_arteries = []
        new_arteries = []
        for a in arteries
            if a in loi.treated_arteries
                push!(failed_arteries, a)
            else
                push!(new_arteries, a)
            end
        end
        failed_lengths = [sum(loi.fcfp.fixed_costs[e] for e in a) for a in failed_arteries]
        sorted_failed_arteries = [x for (_, x) in sort(collect(zip(failed_lengths, failed_arteries)), rev=true)]
        new_lengths = [sum(loi.fcfp.fixed_costs[e] for e in a) for a in new_arteries]
        sorted_new_arteries = [x for (_, x) in sort(collect(zip(new_lengths, new_arteries)), rev=true)]
        return vcat(sorted_new_arteries, sorted_failed_arteries)
    end

    function get_flow_dag(optimization_graph::MultiDiGraph, flows::FlowDictType) :: MultiDiGraph
        edges_to_keep = [e for (e, f) in flows if f > 0]
        return get_edge_sub_graph(optimization_graph, edges_to_keep)
    end

    """
    Remove the sub-optimal flow 2-cycles allowed with the flow conservation.
    Flows dictionary is modified inplace.
    """
    function clean_flow_two_cycles(flows::FlowDictType)
        for (e, flow) in flows
            if flow > 0
                reversed_e = (e[2], e[1], e[3])
                if (reversed_e in keys(flows)) && (flows[reversed_e] > 0)
                    reversed_flow = flows[reversed_e]
                    cycle_flow = min(flow, reversed_flow)
                    flows[e] -= cycle_flow
                    flows[reversed_e] -= cycle_flow
                end
            end
        end
    end

    """
    Break cycles in a given flow dict by subtracting the minimal flow.
    """
    function break_cycles(loi::LocalOptimizerInstance, flows::FlowDictType,
                          start_nodes=nothing::Union{Nothing, Set{VertexType}})
        while true
            flow_dag = get_flow_dag(loi.optimization_graph, flows)
            cycle = find_cycle(flow_dag, start_nodes)
            if cycle !== nothing
                min_flow, min_flow_edge = minimum([(flows[e], e) for e=cycle])
                loi.logger("Breaking cycle :$(min_flow)")
                for e in cycle
                    flows[e] -= min_flow
                end
                flows[min_flow_edge] = zero(FlowType)
            else
                break
            end
        end
    end

    """
    Check that solution is correct with respect to the constraint of the FCFP.
    """
    function check_solution(loi::LocalOptimizerInstance, flows::FlowDictType) :: Bool
        # Check that flow graph is a DAG
        flow_dag = get_flow_dag(loi.optimization_graph, flows)
        c = find_cycle(flow_dag)
        if c !== nothing
            return false
        end
        # Check flow conservation
        tol = 1e-5
        nodes = get_nodes(flow_dag)
        for n in nodes
            out_flow = sum(Array{WeightType}([flows[e] for e in out_edges(flow_dag, n)]))
            in_flow = sum(Array{WeightType}([flows[e] for e in in_edges(flow_dag, n)]))
            if n in loi.fcfp.building_nodes
                continue
            elseif n in loi.fcfp.production_nodes
                if out_flow > (loi.fcfp.production[n] + tol)
                    return false
                end
            else
                if abs(out_flow - in_flow) > tol
                    return false
                end
            end
        end
        return true
    end

    export optimize_locally
end

max_stars_repo_name: Planheat/Planheat-Tool

max_issues_repo_name: Planheat/Planheat-Tool

max_forks_repo_name: Planheat/Planheat-Tool