text:
## Day 16: Packet Decoder ##
############################

day16_test_1 = "8A004A801A8002F478"
day16_test_2 = "620080001611562C8802118E34"
day16_test_3 = "C0015000016115A2E0802F182340"
day16_test_4 = "A0016C880162017C3686B18A3D4780"
day16_input = readlines(open("input-16"))[1]

function input_to_binary_list(input)
    binary_list = []
    for c in input
        binary = string(parse(Int, c, base = 16), base = 2)
        for _ = length(binary):3
            push!(binary_list, 0)
        end
        for d in binary
            push!(binary_list, parse(Int, d, base = 2))
        end
    end
    binary_list
end

function parse_packet!(input)
    version = grab_int!(input, 3)
    typeid = grab_int!(input, 3)
    if typeid == 4
        content = parse_content!(input)
    else
        length_type_id = grab_int!(input, 1)
        if length_type_id == 0
            # The next 15 bits are the total length of the packet, in bits.
            bits_to_parse = grab_int!(input, 15)
            remaining_data_size = length(input)
            content = []
            while length(input) + bits_to_parse > remaining_data_size
                push!(content, parse_packet!(input))
            end
        else
            # The next 11 bits is the amount of contained packets.
            packets_to_parse = grab_int!(input, 11)
            content = []
            for _ = 1:packets_to_parse
                push!(content, parse_packet!(input))
            end
        end
    end
    (version, typeid, content)
end

function grab_int!(input, digits)::Int
    result = 0
    for _ = 1:digits
        result = result * 2 + popfirst!(input)
    end
    result
end

function parse_content!(input)
    value = 0
    while true
        next_digit = grab_int!(input, 5)
        if next_digit >= 16
            # We can continue, but need to strip of a leading 1
            value = value * 16 + next_digit - 16
        else
            # We've reached the end of the content
            value = value * 16 + next_digit
            break
        end
    end
    value
end

function version_sum(packet)
    version, typeid, content = packet
    if typeid == 4
        version
    else
        sum(version_sum.(content)) + version
    end
end

@assert 16 == version_sum(parse_packet!(input_to_binary_list(day16_test_1)))
@assert 12 == version_sum(parse_packet!(input_to_binary_list(day16_test_2)))
@assert 23 == version_sum(parse_packet!(input_to_binary_list(day16_test_3)))
@assert 31 == version_sum(parse_packet!(input_to_binary_list(day16_test_4)))
@assert 1038 == @show version_sum(parse_packet!(input_to_binary_list(day16_input)))

function package_value(packet)
    version, typeid, content = packet
    if typeid == 0
        sum(package_value.(content))
    elseif typeid == 1
        prod(package_value.(content))
    elseif typeid == 2
        minimum(package_value.(content))
    elseif typeid == 3
        maximum(package_value.(content))
    elseif typeid == 4
        content
    elseif typeid == 5
        @assert length(content) == 2
        package_value(content[1]) > package_value(content[2]) ? 1 : 0
    elseif typeid == 6
        @assert length(content) == 2
        package_value(content[1]) < package_value(content[2]) ? 1 : 0
    elseif typeid == 7
        @assert length(content) == 2
        package_value(content[1]) == package_value(content[2]) ? 1 : 0
    else
        throw("Unknown typeid: $typeid")
    end
end

@assert 3 == package_value(parse_packet!(input_to_binary_list("C200B40A82")))
@assert 54 == package_value(parse_packet!(input_to_binary_list("04005AC33890")))
@assert 7 == package_value(parse_packet!(input_to_binary_list("880086C3E88112")))
@assert 9 == package_value(parse_packet!(input_to_binary_list("CE00C43D881120")))
@assert 1 == package_value(parse_packet!(input_to_binary_list("D8005AC2A8F0")))
@assert 0 == package_value(parse_packet!(input_to_binary_list("F600BC2D8F")))
@assert 0 == package_value(parse_packet!(input_to_binary_list("9C005AC2F8F0")))
@assert 1 == package_value(parse_packet!(input_to_binary_list("9C0141080250320F1802104A08")))

@assert 246761930504 == @show package_value(parse_packet!(input_to_binary_list(day16_input)))


text:
# __precompile__()

module MotorPower
using Dierckx
export motorVI
export shigleygear
export planetarygear

"""
out1, out2, Po, eta_M, mass_M, Ro, VmaxmotorVI(input1,input2,Kv,I0;numMotors=1,Case=0,mtype = "hobby")

Author: Kevin Moore
Date: 26 June 2017

# Inputs:
- Q: N-m
- omega: rad/s

# Parameters:
- Kv: rad/s/volt
- Ro: ohms
- I0: amps
- numMotors: #
- case::Bool: 1 is Q, omega in 0 is V, I in

# Constants:

# Outputs:
- V or Q: volts or N-m
- Im or omega: amps or rad/s
- Pshaft: watts
- eta_M: efficiency base 1
- mass_M: kg

Method: Mark Drela 1st order model

Assumptions/Limitiations: thermal resistance increase neglected

References:
"""
function motormass(Im, Kv, Vm, numMotors=1, mtype = "hobby") #Kv in is rad/s/V
    if mtype == "hobby" #TODO: I think that the units for the original data for R0 are incorrect, should be ohms and not mOhms
        k_Imax = [136.24,6.15153, -0.186447, 0.00131255, -0.00304932, 6.66054e-5] #TODO: revisit fit with cross validation
        mass_M  = (k_Imax[1] + k_Imax[2] * Im / numMotors + k_Imax[3] * (Kv * 30 / pi) +
                    k_Imax[4] * (Im / numMotors) .^ 2 + k_Imax[5] * Im / numMotors .* (Kv * 30 / pi) +
                    k_Imax[6] * (Im / numMotors) .^ 2)/1000 #This fit is for Kv - RPM/V
        mass_M = mass_M * numMotors
    end

    if mtype == "astroflight"
        mass_M = 2.464 * Im / numMotors ./ (Kv * 30 / pi) + 0.368 #Kg, kv in RPM/V
        mass_M = mass_M * numMotors
    end

    if mtype == "FB"
        specificPower = 0.6 #kW/Kg
        mass_M = Im * Vm / specificPower #Kg
        mass_M = mass_M * numMotors
    end

    return mass_M
end #motormass

function motorVI(input1, input2, Kv, I0;numMotors=1, Case=0, mtype = "hobby")
    Vmax = 40000 / (Kv * 30 / pi)

    if mtype == "hobby"
        k_params = [0.000144198, 7.73731e-5, 2.17626e-7, -3.77547e-5, -5.66064e-7, 3.1259e-12, -1.29173e-7, 1.4195e-7, 9.93959e-11, -1.6806e-14]
        Ro = k_params[1] + k_params[2] * I0 + k_params[3] * (Kv * 30 / pi) + k_params[4] * I0 .^ 2 +
                k_params[5] * I0 .* (Kv * 30 / pi) + k_params[6] * (Kv * 30 / pi) .^ 2 + k_params[7] * I0 .^ 3 +
                k_params[8] * I0 .^ 2 .* (Kv * 30 / pi) + k_params[9] * I0 .* (Kv * 30 / pi) .^ 2 +
                k_params[10] * (Kv * 30 / pi) .^ 3
        tune = Ro * 0.2
        Ro = Ro - tune
        Ro = 1
    end

    if mtype == "astroflight" || mtype == "FB"
        Ro = .0467 * I0 .^ -1.892 #excel fit, TODO: assume FM motor is similar
    end

    if Case == 1
        Q = input1
        omega = input2
        #R2= Ro.*alpha.*delTmax./Imax^2  # include temperature and peak effects
        Ri= Ro #+R2.*I^2

        # Current Model
        Ii = I0 + Kv .* Q #Kv in rad/s/V

        # Voltage Model
        Vm = omega ./ Kv
        V = Vm + Ri .* Ii

        # Derived Relations
        Pshaft = ((V - (omega ./ Kv)) .* (1.0 ./ Ri) - I0) .* (omega ./ Kv)
        eta_M = (1.0 - ((I0 .* Ri) ./ (V - omega ./ Kv))) .* (omega ./ (V .* Kv))
        Pm = Q .* omega ./ eta_M
        Im = Pm ./ V

        # Motor mass
        mass_M = motormass(Im, Kv, V, numMotors, mtype)

        # #Prevent motor from generating power
        # if Vm < 0
        #     Vm = 0
        # end
        # if Im <0
        #     Im = 0
        # end

        return  V, Im, Pshaft, eta_M, mass_M, Ro, Vmax
    else #Adapted from University of Auburn Master's thesis
        Vm = input1
        Im = input2

        #R2= Ro.*alpha.*delTmax./Imax^2  # include temperature and peak effects
        Ri= Ro #+R2.*I^2
        N_Pol = 12
        N_Phase = 3
        omegaMax = Vm * Kv
        L = 2 / (N_Pol * omegaMax) * (Vm - Ri * I0)
        Tem = (N_Phase * N_Pol) / 2 * L * Im
        omega = (Vm / (N_Pol * L / 2) - Ri / (N_Phase * (N_Pol * L / 2)^2) * Tem)
        Kt = 1 / Kv
        Q = Kt * Im
        Pi = Vm * Im
        Po = Q * omega #(Vm - Im * Ri) * (Im - I0) #TODO: this can become very wrong and be exploited by the optimizer: use Q*omega for now
        eta_M = Po / Pi

        mass_M = motormass(Im, Kv, Vm, numMotors, mtype)

        return  Q, omega, Po, eta_M, mass_M, Ro, Vmax
    end
end #motorVI

"""
motorcontroller

Author: Kevin Moore
Date: 26 June 2017

Updates:

Units: SI

Inputs:
Vm: volts
Im: amps

Parameters:

Constants:

Outputs:
eta_E: efficiency base 1
mass_E: kg

Method: Empirical fits, ohms law

Assumptions/Limitiations: accuracy of empirical data

References: http://www.drivecalc.de
"""

function ESC(Vm, Im; esctype="hobby")
    if esctype=="hobby" #From drivecalc, also consistent with astroflight hobby esc's
        mass_E = (7.1261 * exp(.0114 .* Im)) / 1000 #Kg
        R_E = 130.36 * Im ^ (-0.9) / 1000 #ohms
        eta_E = (Vm - (Im * R_E)) / Vm
    elseif esctype=="FB"
        specificPower = 4444.0 #W/Kg
        mass_E = Vm * Im / specificPower #Assume linear scaling
        eta_E = 0.97 #TODO...
    elseif esctype=="AstroFlightHighVoltage"
        specificPower = 22059.0 #W/Kg
        mass_E = Vm * Im / specificPower #Assume linear scaling
        eta_E = 0.97 #TODO...
    # else #TODO... add more
    end

    return mass_E, eta_E
end #esc

"""
motorVI

Author: Kevin Moore
Date: 26 June 2017

Updates:

Units: SI

Inputs:
RPMi: RPM
RPMo: RPM
Pin: watts into pinion
Tprop: N prop thrust
Dp: m pinion diameter
F: m gears' face width
Pd: teeth/meter
PAR: parameters

Parameters:
PAR[:gearPitch] = 20*pi/180
PAR[:Cp] = 2300 sqrt.(psi) #steel # must review entire gearbox model if changed
PAR[:Hbp] = 240 #brindell hardness pinion
PAR[:Hbg] = 200 #brindell hardness gear
PAR[:rhoSteel] = 8050 #kg/m3
PAR[:rhoAl] = 2800 #kg/m3
PAR[:SF] = 1.25 #safety factor
PAR[:bearPerRod] = 2

Constants:

Outputs:
M_b: kg
plateArea: m2 (for gearbox drag)
Qgb: Nm output torque
cpb: pinion bending constraint
cgb: gear bending constraint
cpw: pinion wear constraint
cgw: gear wear constraint

Method: Shigely

Assumptions/Limitiations: changing materials requires changing much of the model

References: Shigley 14-4
"""
function shigleygear(RPMi, RPMo, Pin, Tprop, Dp, F, Pd, gearPitch, Cp, Hbp, Hbg,
                        rhoSteel, rhoAl, SF, bearPerRod, Ndays)

    #--- Unpack parameters ---#
    # gearPitch =PAR[:gearPitch]
    # Cp = PAR[:Cp]#steel # must review entire gearbox model if changed
    # Hbp = PAR[:Hbp]#brindell hardness pinion
    # Hbg = PAR[:Hbg]#brindell hardness gear
    # rhoSteel = PAR[:rhoSteel]#kg/m3
    # rhoAl = PAR[:rhoAl]#kg/m3
    # SF =PAR[:SF] #safety factor
    # bearPerRod = PAR[:bearPerRod]
    # Ndays = PAR[:Ndays]

    # ----- Gear Calculation ----- #
    Dp = Dp * 39.1901 # pinion diameter inches
    F = F * 39.1901 # gear face width inches
    Pd = Pd / 39.1901 # diametrical pitch teeth/inch
    Tprop = Tprop * 0.224809

    Np = RPMi * 60 * 24 * Ndays #full 90 days of running
    Ng = RPMo * 60 * 24 * Ndays #full 90 days of running

    Dg = RPMi / RPMo*Dp

    HP = Pin / 745.7

    V = pi * Dp * RPMi / 12 #ft/min

    Wt = 33000 * HP / V #lbf

    Ko = 1
    Qv = 6

    B =  0.25 * (12 - Qv)^(2/3)
    A = 50 + 56 * (1 - B)

    Kv = ((A + sqrt.(V)) / A)^B

    #lewis form factor table 14-2
    teeth = [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 24, 26, 28, 30,
                34, 38, 43, 50, 60, 75, 100, 150, 300, 400]
    Y = [0.245, 0.261, 0.277, 0.29, 0.296, 0.303, 0.309, 0.312, 0.322,
            0.328, 0.331, 0.337, 0.346, 0.353, 0.359, 0.371, 0.384,
            0.397, 0.409, 0.422, 0.435, 0.447, 0.46, 0.472, 0.48]
    spl_Y = Dierckx.Spline1D(teeth, Y) #TODO: replace with akima spline?
    Tp = pi * Dp * Pd
    Tg = pi * Dg * Pd
    Yp = spl_Y(Tp) #lewis form factor
    Yg = spl_Y(Tg) #lewis form factor

    KSp = 1.192 * (F * sqrt.(Yp) / Pd) ^ .0535 #size factor
    KSg = 1.192 * (F * sqrt.(Yg) / Pd) ^ .0535 #size factor

    Cmc = 1 #uncrowned
    Cpfp = F / 10 / Dp - 0.025
    Cpfg = F / 10 / Dg - 0.025
    Cpm = 1
    Cma = .27
    Ce = 1

    Kmp = 1 + Cmc * (Cpfp * Cpm + Cma * Ce)
    Kmg = 1 + Cmc * (Cpfg * Cpm + Cma * Ce)

    Kb = 1 #const thickness gears

    Ynp = 1.3558 * Np ^ -.0178
    Yng = 1.3558 * Ng ^ -.0178

    #Table 14.10
    #reliability = [0.9999, 0.999, 0.99, 0.9, 0.5]
    #Krel = [1.5, 1.25, 1, 0.85, 0.7]
    Kr = 1 # for R = 0.99
    mg = RPMi / RPMo
    Kt = 1
    Ks = 1
    Cf = 1
    mn = 1
    Jp = 0.35 #TODO may need to fit a curve...
    Jg = 0.4

    I = cos(gearPitch) * sin(gearPitch) / (2 * mn) * mg / (mg + 1)

    Stp = 102 * Hbp + 16400 #bending strength grade 2 hardened steel
    Stg = 102 * Hbp + 16400

    Scp = 349 * Hbp + 34300 #contact fatigue strength grade 2 hardened steel
    Scg = 349 * Hbg + 34300

    Znp = 1.4488 * Np ^ -0.023
    Zng = 1.4488 * Ng ^ -0.023

    Aprime = 8.98e-3 * Hbp / Hbg - 8.29e-3 # for 1.2<Hbp / Hbg<1.7
    Ch = 1 + Aprime * (mg - 1)

    #Pinion Tooth bending
    sigmap = Wt * Ko * Kv * Ks * Pd / F * Kmp * Kb / Jp
    SFbp = Stp * Ynp / (Kt * Kr) / sigmap

    #Gear Tooth Bending
    sigmag = Wt * Ko * Kv * Ks * Pd / F * Kmg * Kb / Jg
    SFbg = Stg * Yng / (Kt * Kr) / sigmag

    #Pinion Tooth Wear
    # println("Wt: $Wt Ko: $Ko Kv: $Kv Ks: $Ks Kmp: $Kmp Dp: $Dp F: $F Cf: $Cf I: $I ")
    sigmaCp = Cp * (Wt * Ko * Kv * Ks * Kmp / Dp / F * Cf / I) ^ 0.5
    SFwp = Scp * Znp / (Kt * Kr) / sigmaCp

    #Gear Tooth Wear
    sigmaCg = Cp * (Wt * Ko * Kv * Ks * Kmg / Dp / F * Cf / I) ^ 0.5
    SFwg = Scg * Zng / (Kt * Kr) / sigmaCg

    #Gear Masses
    gearVolume = pi * F / 39.1901 * ((Dp / 2 / 39.1901) ^ 2 + (Dg / 2 / 39.1901) ^ 2)
    gearMass = gearVolume * rhoSteel #kg

    # ----- Bearing Calculation ----- #
    xDp = Np / 1e6
    xDg = Ng / 1e6

    #C10p = 1.2*Tprop / bearPerRod / 10*(xDp / (.02+4.439(1-.99)^(1 / 1.483)))^(1 / 3) #Pinion supported by motor
    C10g = 1.2 * Tprop * SF / bearPerRod * (xDg / (0.02 + 4.439 * (1 - 0.99) ^ (1 / 1.483))) ^ (1 / 3)

    bearingID = [10, 12, 15, 17, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95] #mm
    bearingOD = [30, 32, 35, 40, 47, 52, 62, 72, 80, 85, 90, 100, 110, 120, 125, 130, 140, 150, 160, 170] #mm
    bearingWidth = [9, 10, 11, 12, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 28, 30, 32] ##
    deepGrooveC10Rating = [5.07, 6.98, 7.8, 9.56, 12.7, 14.0, 19.5, 25.5, 30.7, 33.2, 35.1, 43.6, 47.5, 55.9, 61.8, 66.3, 70.2, 83.2, 95.6, 108.0] #Kn

    spl_RatVWidth = Dierckx.Spline1D(deepGrooveC10Rating, bearingWidth)
    spl_WidthVID = Dierckx.Spline1D(bearingWidth, bearingID)
    spl_WidthVOD = Dierckx.Spline1D(bearingWidth, bearingOD)

    #bear_wp = spl_RatVWidth(C10p) / 1000
    bear_wg = spl_RatVWidth(C10g) / 1000

    #bear_IDp = spl_WidthVID(bear_wp) / 1000
    bear_IDg = spl_WidthVID(bear_wg) / 1000

    #bear_IDg = spl_WidthVID(bear_wg) / 1000
    bear_ODg = spl_WidthVOD(bear_wg) / 1000

    #bearMassp = pi * ((bearODp / 2) ^ 2 - (bearIDp / 2) ^ 2) * bear_wp * rhoSteel
    bearMassg = pi * ((bear_ODg / 2) ^ 2 - (bear_IDg / 2) ^ 2) * bear_wg * rhoSteel
    bearMass = bearPerRod * (bearMassg) #kg, 2 bearings of each

    # ----- Shaft Calculation ----- #
    shaftLen = bear_IDg * 6 # 4D spacing between bearings and 2D to attach prop TODO
    shaftMass = pi * (bear_IDg / 2) ^ 2 * shaftLen * rhoAl #kg

    #Mounting Plate Mass (designed like y shape)
    plateVol = (Dg + Dp) / 39.1901 * bear_ODg * 1.2 * bear_wg + (bear_wg * bear_ODg * 1.2 * Dg / 39.1901)
    plateArea = (Dg + Dp) / 39.1901 * bear_ODg * 1.2 #m^2
    #slightly wider than the bearing, as long as the gears, same depth as bearing
    mountMass = plateVol * rhoAl #Kg TODO estimate stress?

    # ----- Outputs ----- #
    gearMass = gearMass * 0.6 * 0.8 #Assume 40% is drilled out and inner sections are thinner by ~15% TODO:
    M_b = gearMass #+shaftMass#+bearMass+gearMass+mountMass #gearbox mass in Kg

    #constraints
    cpb = (SF - SFbp) #/SFbp
    cgb = (SF - SFbg) #/SFbg
    cpw = (SF - SFwp) #/SFwp
    cgw = (SF - SFwg) #/SFwg

    eta_GB = 0.95 #TODO

    # # set mass and everything else to 0 for gear ratio of 1 (no gearbox)
    # sharpness = 1000.0
    # x_smooth = mg*sharpness-sharpness
    # if  x_smooth < 1.0 && x_smooth > 0.0
    #     correction = 3*x_smooth^2 - 2*x_smooth^3
    #     M_b = M_b*correction
    #     plateArea = plateArea*correction
    #     eta_GB = eta_GB+(1-eta_GB)*(1-correction)
    #     # cpb = cpb*correction
    #     # cgb = cgb*correction
    #     # cpw = cpw*correction
    #     # cgw = cgw*correction
    #
    # end
    # x_smooth2 = 1/mg*sharpness-sharpness
    # if  x_smooth2 < 1.0 && x_smooth2 > 0.0
    #     correction = 3*x_smooth2^2 - 2*x_smooth2^3
    #     M_b = M_b*correction
    #     plateArea = plateArea*correction
    #     eta_GB = eta_GB+(1-eta_GB)*(1-correction)
    #     # cpb = cpb*correction
    #     # cgb = cgb*correction
    #     # cpw = cpw*correction
    #     # cgw = cgw*correction
    #
    # end

    Qgb = Pin / (RPMo * pi / 30) * eta_GB

    return eta_GB, M_b, plateArea, Qgb, cpb, cgb, cpw, cgw
end #function

function planetarygear(Qin, ratio)
    M_b = (30.818 .* Qin + 1.1846) / 1000 #kg
    eta_GB = 0.95 #TODO

    # # set mass and everything else to 0 for gear ratio of 1 (no gearbox)
    # sharpness = 1000.0
    # x_smooth = ratio*sharpness-sharpness
    # if  x_smooth < 1.0 && x_smooth > 0.0
    #     correction = 3*x_smooth^2 - 2*x_smooth^3
    #     M_b = M_b*correction
    #     eta_GB = eta_GB+(1-eta_GB)*(1-correction)
    # end
    #
    # x_smooth2 = 1/ratio*sharpness-sharpness
    # if  x_smooth2 < 1.0 && x_smooth2 > 0.0
    #     correction = 3*x_smooth2^2 - 2*x_smooth2^3
    #     M_b = M_b*correction
    #     eta_GB = eta_GB+(1-eta_GB)*(1-correction)
    # end

    return eta_GB, M_b
end

function wire_1(volts, amps, distance, coreType, wireDensity)

    if coreType==1 #ie a single wire

        # line fits from data from http://www.engineeringtoolbox.com/wire-gauges-d_419.html
        diameter = (0.0359 .* amps + 0.6777) / 1000 #for diameter in mm, convert to meters  #R² = 0.99407
        wireVolume = pi * (diameter ./ 2) .^ 2 .* distance
        wireMass = wireVolume * wireDensity
        #copperWeightN = wireMass .* g
        resistancePerKm = 22.447 * (diameter * 1000) .^ -2.007 #for diameter in mm R² = 0.99948
        resistance = resistancePerKm / 1000 .* distance #convert to per meter and multiply by meter

        powerTransEfficiency = 1 - amps .* resistance ./ volts
    else
        println("TODO: inclue different core types")
    end

    return wireMass, powerTransEfficiency
end #wire_1

end # module

ext: jl


text:
using DelimitedFiles

function ReadData(fname::AbstractString)
    dataStr = readdlm(fname, String)

    data = zeros(Int, length(dataStr), length(dataStr[1]))
    for i = 1:length(dataStr)
        data[i,:] = map(x -> parse(Int, x), split(dataStr[i],""))
    end

    return data
end

function FindMinima(data::Matrix)
    vOffsetU = circshift(data, (-1,0)) # data shifted up one
    hOffsetL = circshift(data, (0,-1)) # data shifted to the left
    vOffsetD = circshift(data, (1,0))  # data shifted up one
    hOffsetR = circshift(data, (0,1))  # data shifted to the left

    # set the edges to be 9 (the maximum possible value)
    vOffsetU[size(data)[1],:] = fill(9, size(data)[2])
    hOffsetL[:,size(data)[2]] = fill(9, size(data)[1])
    vOffsetD[1,:] = fill(9, size(data)[2])
    hOffsetR[:,1] = fill(9, size(data)[1])

    return intersect([findall(x -> x<0, data-offset) for offset in [vOffsetU, hOffsetL, vOffsetD, hOffsetR]]...)
end

function Part1(data::Matrix)
    minimaPositions = FindMinima(data)

    minima = map(x -> getindex(data, x), minimaPositions)

    return sum(minima)+length(minima)
end

function Part2(data::Matrix)
    vOffsetU = circshift(data, (-1,0)) # data shifted up one
    hOffsetL = circshift(data, (0,-1)) # data shifted to the left
    vOffsetD = circshift(data, (1,0))  # data shifted up one
    hOffsetR = circshift(data, (0,1))  # data shifted to the left

    # set the edges to be 9 (the maximum possible value)
    vOffsetU[size(data)[1],:] = fill(9, size(data)[2])
    hOffsetL[:,size(data)[2]] = fill(9, size(data)[1])
    vOffsetD[1,:] = fill(9, size(data)[2])
    hOffsetR[:,1] = fill(9, size(data)[1])

    minimaPositions = intersect([findall(x -> x<0,  data-offset) for offset in [vOffsetU, hOffsetL, vOffsetD, hOffsetR]]...)
    maximaPositions = intersect([findall(x -> x>=0, data-offset) for offset in [vOffsetU, hOffsetL, vOffsetD, hOffsetR]]...) # this will catch saddle points

    basinSizes = []

    for x in minimaPositions
        # start building basin
        currentBasin = [x]

        basinSize = 0

        Q = [x]
        explored = [x]
        while ! isempty(Q)
            v = pop!(Q)
            if getindex(data, v) == 9
                continue
            else
                basinSize += 1
            end
            i, j = v[1], v[2]
            for k=-1:2:1
                if CartesianIndex(i+k, j) ∉ explored && 1 <= i+k <= size(data)[1]
                    push!(explored, CartesianIndex(i+k, j))
                    push!(Q, CartesianIndex(i+k, j))
                end

                if CartesianIndex(i, j+k) ∉ explored && 1 <= j+k <= size(data)[2]
                    push!(explored, CartesianIndex(i, j+k))
                    push!(Q, CartesianIndex(i, j+k))
                end
            end
        end

        push!(basinSizes, basinSize)
    end

    sort!(basinSizes, rev=true)

    return basinSizes[1]*basinSizes[2]*basinSizes[3]

end

#data = ReadData("../data/test9.txt")
data = ReadData("../data/day9.txt")

println(Part2(data))


text:
using ConstrainedDynamics
using ConstrainedDynamics: ∂g∂ʳpos, ∂g∂ʳvel, discretizestate!, LVᵀmat, vrotate
using ForwardDiff
using Rotations
using Rotations: params
using LinearAlgebra

function transtest3()
    Δt = 0.01
    xa1 = rand(3)
    qa1 = rand(UnitQuaternion)
    xb1 = rand(3)
    qb1 = rand(UnitQuaternion)

    va1 = rand(3)
    ωa1 = rand(3)
    vb1 = rand(3)
    ωb1 = rand(3)

    va2 = rand(3)
    ωa2 = rand(3)
    vb2 = rand(3)
    ωb2 = rand(3)

    pa = rand(3)
    pb = rand(3)

    origin = Origin{Float64}()
    body1 = Body(Box(1., 1., 1., 1.))
    body2 = Body(Box(1., 1., 1., 1.))

    oc1 = EqualityConstraint(Floating(origin, body1))
    oc2 = EqualityConstraint(Floating(origin, body2))
    joint1 = EqualityConstraint(ConstrainedDynamics.Translational3{Float64}(body1, body2, p1=pa, p2=pb))

    mech = Mechanism(origin, [body1;body2], [oc1;oc2;joint1])
    discretizestate!(body1,xa1,qa1,va1,va2,ωa1,ωa2,Δt)
    discretizestate!(body2,xb1,qb1,vb1,vb2,ωb1,ωb2,Δt)

    res = ForwardDiff.jacobian(transfunc3pos, [getxqkvector(body1.state);getxqkvector(body2.state);pa;pb])
    X1 = res[1:3,1:3]
    Q1 = res[1:3,4:7] * LVᵀmat(getqk(body1.state))
    X2 = res[1:3,8:10]
    Q2 = res[1:3,11:14] * LVᵀmat(getqk(body2.state))

    n11 = norm(X1 - ∂g∂ʳpos(mech, joint1, 1)[1:3,1:3])
    n12 = norm(Q1 - ∂g∂ʳpos(mech, joint1, 1)[1:3,4:6])
    n21 = norm(X2 - ∂g∂ʳpos(mech, joint1, 2)[1:3,1:3])
    n22 = norm(Q2 - ∂g∂ʳpos(mech, joint1, 2)[1:3,4:6])

    res = ForwardDiff.jacobian(transfunc3vel, [getsolestimate(body1.state);getsolestimate(body2.state);pa;pb])
    V1 = res[1:3,4:6]
    W1 = res[1:3,11:13]
    V2 = res[1:3,17:19]
    W2 = res[1:3,24:26]

    n31 = norm(V1 - ∂g∂ʳvel(mech, joint1, 1)[1:3,1:3])
    n32 = norm(W1 - ∂g∂ʳvel(mech, joint1, 1)[1:3,4:6])
    n41 = norm(V2 - ∂g∂ʳvel(mech, joint1, 2)[1:3,1:3])
    n42 = norm(W2 - ∂g∂ʳvel(mech, joint1, 2)[1:3,4:6])

    # display((n11, n12, n21, n22))
    # display((n31, n32, n41, n42))
    return n11+n12+n21+n22+n31+n32+n41+n42
end

function transtest2()
    Δt = 0.01
    xa1 = rand(3)
    qa1 = rand(UnitQuaternion)
    xb1 = rand(3)
    qb1 = rand(UnitQuaternion)

    va1 = rand(3)
    ωa1 = rand(3)
    vb1 = rand(3)
    ωb1 = rand(3)

    va2 = rand(3)
    ωa2 = rand(3)
    vb2 = rand(3)
    ωb2 = rand(3)

    pa = rand(3)
    pb = rand(3)

    v = rand(3)

    origin = Origin{Float64}()
    body1 = Body(Box(1., 1., 1., 1.))
    body2 = Body(Box(1., 1., 1., 1.))

    oc1 = EqualityConstraint(Floating(origin, body1))
    oc2 = EqualityConstraint(Floating(origin, body2))
    joint1 = EqualityConstraint(ConstrainedDynamics.Translational2{Float64}(body1, body2, p1=pa, p2=pb, axis=v))
    V12 = joint1.constraints[1].V12

    mech = Mechanism(origin, [body1;body2], [oc1;oc2;joint1])
    discretizestate!(body1,xa1,qa1,va1,va2,ωa1,ωa2,Δt)
    discretizestate!(body2,xb1,qb1,vb1,vb2,ωb1,ωb2,Δt)

    res = ForwardDiff.jacobian(transfunc2pos, [getxqkvector(body1.state);getxqkvector(body2.state);pa;pb;V12[1,:];V12[2,:]])
    X1 = res[1:2,1:3]
    Q1 = res[1:2,4:7] * LVᵀmat(getqk(body1.state))
    X2 = res[1:2,8:10]
    Q2 = res[1:2,11:14] * LVᵀmat(getqk(body2.state))

    n11 = norm(X1 - ∂g∂ʳpos(mech, joint1, 1)[1:2,1:3])
    n12 = norm(Q1 - ∂g∂ʳpos(mech, joint1, 1)[1:2,4:6])
    n21 = norm(X2 - ∂g∂ʳpos(mech, joint1, 2)[1:2,1:3])
    n22 = norm(Q2 - ∂g∂ʳpos(mech, joint1, 2)[1:2,4:6])

    res = ForwardDiff.jacobian(transfunc2vel, [getsolestimate(body1.state);getsolestimate(body2.state);pa;pb;V12[1,:];V12[2,:]])
    V1 = res[1:2,4:6]
    W1 = res[1:2,11:13]
    V2 = res[1:2,17:19]
    W2 = res[1:2,24:26]

    n31 = norm(V1 - ∂g∂ʳvel(mech, joint1, 1)[1:2,1:3])
    n32 = norm(W1 - ∂g∂ʳvel(mech, joint1, 1)[1:2,4:6])
    n41 = norm(V2 - ∂g∂ʳvel(mech, joint1, 2)[1:2,1:3])
    n42 = norm(W2 - ∂g∂ʳvel(mech, joint1, 2)[1:2,4:6])

    # display((n11, n12, n21, n22))
    # display((n31, n32, n41, n42))
    return n11+n12+n21+n22+n31+n32+n41+n42
end

function transtest1()
    Δt = 0.01
    xa1 = rand(3)
    qa1 = rand(UnitQuaternion)
    xb1 = rand(3)
    qb1 = rand(UnitQuaternion)

    va1 = rand(3)
    ωa1 = rand(3)
    vb1 = rand(3)
    ωb1 = rand(3)

    va2 = rand(3)
    ωa2 = rand(3)
    vb2 = rand(3)
    ωb2 = rand(3)

    pa = rand(3)
    pb = rand(3)

    v = rand(3)

    origin = Origin{Float64}()
    body1 = Body(Box(1., 1., 1., 1.))
    body2 = Body(Box(1., 1., 1., 1.))

    oc1 = EqualityConstraint(Floating(origin, body1))
    oc2 = EqualityConstraint(Floating(origin, body2))
    joint1 = EqualityConstraint(ConstrainedDynamics.Translational1{Float64}(body1, body2, p1=pa, p2=pb, axis=v))
    v = joint1.constraints[1].V3'

    mech = Mechanism(origin, [body1;body2], [oc1;oc2;joint1])
    discretizestate!(body1,xa1,qa1,va1,va2,ωa1,ωa2,Δt)
    discretizestate!(body2,xb1,qb1,vb1,vb2,ωb1,ωb2,Δt)

    res = ForwardDiff.gradient(transfunc1pos, [getxqkvector(body1.state);getxqkvector(body2.state);pa;pb;v])
    X1 = res[1:3]'
    Q1 = res[4:7]' * LVᵀmat(getqk(body1.state))
    X2 = res[8:10]'
    Q2 = res[11:14]' * LVᵀmat(getqk(body2.state))

    n11 = norm(X1 - ∂g∂ʳpos(mech, joint1, 1)[1:3]')
    n12 = norm(Q1 - ∂g∂ʳpos(mech, joint1, 1)[4:6]')
    n21 = norm(X2 - ∂g∂ʳpos(mech, joint1, 2)[1:3]')
    n22 = norm(Q2 - ∂g∂ʳpos(mech, joint1, 2)[4:6]')

    res = ForwardDiff.gradient(transfunc1vel, [getsolestimate(body1.state);getsolestimate(body2.state);pa;pb;v])
    V1 = res[4:6]'
    W1 = res[11:13]'
    V2 = res[17:19]'
    W2 = res[24:26]'

    n31 = norm(V1 - ∂g∂ʳvel(mech, joint1, 1)[1:3]')
    n32 = norm(W1 - ∂g∂ʳvel(mech, joint1, 1)[4:6]')
    n41 = norm(V2 - ∂g∂ʳvel(mech, joint1, 2)[1:3]')
    n42 = norm(W2 - ∂g∂ʳvel(mech, joint1, 2)[4:6]')

    # display((n11, n12, n21, n22))
    # display((n31, n32, n41, n42))
    return n11+n12+n21+n22+n31+n32+n41+n42
end

function rottest3()
    Δt = 0.01
    xa1 = rand(3)
    qa1 = rand(UnitQuaternion)
    xb1 = rand(3)
    qb1 = rand(UnitQuaternion)

    va1 = rand(3)
    ωa1 = rand(3)
    vb1 = rand(3)
    ωb1 = rand(3)

    va2 = rand(3)
    ωa2 = rand(3)
    vb2 = rand(3)
    ωb2 = rand(3)

    qoffset = rand(UnitQuaternion)

    origin = Origin{Float64}()
    body1 = Body(Box(1., 1., 1., 1.))
    body2 = Body(Box(1., 1., 1., 1.))

    oc1 = EqualityConstraint(Floating(origin, body1))
    oc2 = EqualityConstraint(Floating(origin, body2))
    joint1 = EqualityConstraint(ConstrainedDynamics.Rotational3{Float64}(body1, body2, qoffset = qoffset))

    mech = Mechanism(origin, [body1;body2], [oc1;oc2;joint1])
    discretizestate!(body1,xa1,qa1,va1,va2,ωa1,ωa2,Δt)
    discretizestate!(body2,xb1,qb1,vb1,vb2,ωb1,ωb2,Δt)

    res = ForwardDiff.jacobian(rotfunc3pos, [getxqkvector(body1.state);getxqkvector(body2.state);params(qoffset)])
    X1 = res[1:3,1:3]
    Q1 = res[1:3,4:7] * LVᵀmat(getqk(body1.state))
    X2 = res[1:3,8:10]
    Q2 = res[1:3,11:14] * LVᵀmat(getqk(body2.state))

    n11 = norm(X1 - ∂g∂ʳpos(mech, joint1, 1)[1:3,1:3])
    n12 = norm(Q1 - ∂g∂ʳpos(mech, joint1, 1)[1:3,4:6])
    n21 = norm(X2 - ∂g∂ʳpos(mech, joint1, 2)[1:3,1:3])
    n22 = norm(Q2 - ∂g∂ʳpos(mech, joint1, 2)[1:3,4:6])

    res = ForwardDiff.jacobian(rotfunc3vel, [getsolestimate(body1.state);getsolestimate(body2.state);params(qoffset)])
    V1 = res[1:3,4:6]
    W1 = res[1:3,11:13]
    V2 = res[1:3,17:19]
    W2 = res[1:3,24:26]

    n31 = norm(V1 - ∂g∂ʳvel(mech, joint1, 1)[1:3,1:3])
    n32 = norm(W1 - ∂g∂ʳvel(mech, joint1, 1)[1:3,4:6])
    n41 = norm(V2 - ∂g∂ʳvel(mech, joint1, 2)[1:3,1:3])
    n42 = norm(W2 - ∂g∂ʳvel(mech, joint1, 2)[1:3,4:6])

    # display((n11, n12, n21, n22))
    # display((n31, n32, n41, n42))
    return n11+n12+n21+n22+n31+n32+n41+n42
end

function rottest2()
    Δt = 0.01
    xa1 = rand(3)
    qa1 = rand(UnitQuaternion)
    xb1 = rand(3)
    qb1 = rand(UnitQuaternion)

    va1 = rand(3)
    ωa1 = rand(3)
    vb1 = rand(3)
    ωb1 = rand(3)

    va2 = rand(3)
    ωa2 = rand(3)
    vb2 = rand(3)
    ωb2 = rand(3)

    qoffset = rand(UnitQuaternion)

    v = rand(3)

    origin = Origin{Float64}()
    body1 = Body(Box(1., 1., 1., 1.))
    body2 = Body(Box(1., 1., 1., 1.))

    oc1 = EqualityConstraint(Floating(origin, body1))
    oc2 = EqualityConstraint(Floating(origin, body2))
    joint1 = EqualityConstraint(ConstrainedDynamics.Rotational2{Float64}(body1, body2, axis = v, qoffset = qoffset))
    V12 = joint1.constraints[1].V12

    mech = Mechanism(origin, [body1;body2], [oc1;oc2;joint1])
    discretizestate!(body1,xa1,qa1,va1,va2,ωa1,ωa2,Δt)
    discretizestate!(body2,xb1,qb1,vb1,vb2,ωb1,ωb2,Δt)

    res = ForwardDiff.jacobian(rotfunc2pos, [getxqkvector(body1.state);getxqkvector(body2.state);params(qoffset);V12[1,:];V12[2,:]])
    X1 = res[1:2,1:3]
    Q1 = res[1:2,4:7] * LVᵀmat(getqk(body1.state))
    X2 = res[1:2,8:10]
    Q2 = res[1:2,11:14] * LVᵀmat(getqk(body2.state))

    n11 = norm(X1 - ∂g∂ʳpos(mech, joint1, 1)[1:2,1:3])
    n12 = norm(Q1 - ∂g∂ʳpos(mech, joint1, 1)[1:2,4:6])
    n21 = norm(X2 - ∂g∂ʳpos(mech, joint1, 2)[1:2,1:3])
    n22 = norm(Q2 - ∂g∂ʳpos(mech, joint1, 2)[1:2,4:6])

    res = ForwardDiff.jacobian(rotfunc2vel, [getsolestimate(body1.state);getsolestimate(body2.state);params(qoffset);V12[1,:];V12[2,:]])
    V1 = res[1:2,4:6]
    W1 = res[1:2,11:13]
    V2 = res[1:2,17:19]
    W2 = res[1:2,24:26]

    n31 = norm(V1 - ∂g∂ʳvel(mech, joint1, 1)[1:2,1:3])
    n32 = norm(W1 - ∂g∂ʳvel(mech, joint1, 1)[1:2,4:6])
    n41 = norm(V2 - ∂g∂ʳvel(mech, joint1, 2)[1:2,1:3])
    n42 = norm(W2 - ∂g∂ʳvel(mech, joint1, 2)[1:2,4:6])

    # display((n11, n12, n21, n22))
    # display((n31, n32, n41, n42))
    return n11+n12+n21+n22+n31+n32+n41+n42
end

function rottest1()
    Δt = 0.01
    xa1 = rand(3)
    qa1 = rand(UnitQuaternion)
    xb1 = rand(3)
    qb1 = rand(UnitQuaternion)

    va1 = rand(3)
    ωa1 = rand(3)
    vb1 = rand(3)
    ωb1 = rand(3)

    va2 = rand(3)
    ωa2 = rand(3)
    vb2 = rand(3)
    ωb2 = rand(3)

    qoffset = rand(UnitQuaternion)

    v = rand(3)

    origin = Origin{Float64}()
    body1 = Body(Box(1., 1., 1., 1.))
    body2 = Body(Box(1., 1., 1., 1.))

    oc1 = EqualityConstraint(Floating(origin, body1))
    oc2 = EqualityConstraint(Floating(origin, body2))
    joint1 = EqualityConstraint(ConstrainedDynamics.Rotational1{Float64}(body1, body2, axis = v, qoffset = qoffset))
    V3 = joint1.constraints[1].V3'

    mech = Mechanism(origin, [body1;body2], [oc1;oc2;joint1])
    discretizestate!(body1,xa1,qa1,va1,va2,ωa1,ωa2,Δt)
    discretizestate!(body2,xb1,qb1,vb1,vb2,ωb1,ωb2,Δt)

    res = ForwardDiff.gradient(rotfunc1pos, [getxqkvector(body1.state);getxqkvector(body2.state);params(qoffset);V3])
    X1 = res[1:3]'
    Q1 = res[4:7]' * LVᵀmat(getqk(body1.state))
    X2 = res[8:10]'
    Q2 = res[11:14]' * LVᵀmat(getqk(body2.state))

    n11 = norm(X1 - ∂g∂ʳpos(mech, joint1, 1)[1:3]')
    n12 = norm(Q1 - ∂g∂ʳpos(mech, joint1, 1)[4:6]')
    n21 = norm(X2 - ∂g∂ʳpos(mech, joint1, 2)[1:3]')
    n22 = norm(Q2 - ∂g∂ʳpos(mech, joint1, 2)[4:6]')

    res = ForwardDiff.gradient(rotfunc1vel, [getsolestimate(body1.state);getsolestimate(body2.state);params(qoffset);V3])
    V1 = res[4:6]'
    W1 = res[11:13]'
    V2 = res[17:19]'
    W2 = res[24:26]'

    n31 = norm(V1 - ∂g∂ʳvel(mech, joint1, 1)[1:3]')
    n32 = norm(W1 - ∂g∂ʳvel(mech, joint1, 1)[4:6]')
    n41 = norm(V2 - ∂g∂ʳvel(mech, joint1, 2)[1:3]')
    n42 = norm(W2 - ∂g∂ʳvel(mech, joint1, 2)[4:6]')

    # display((n11, n12, n21, n22))
    # display((n31, n32, n41, n42))
    return n11+n12+n21+n22+n31+n32+n41+n42
end

for i=1:10
    @test isapprox(transtest3(), 0.0; atol = 1e-8)
    @test isapprox(transtest2(), 0.0; atol = 1e-8)
    @test isapprox(transtest1(), 0.0; atol = 1e-8)
    @test isapprox(rottest3(), 0.0; atol = 1e-8)
    @test isapprox(rottest2(), 0.0; atol = 1e-8)
    @test isapprox(rottest1(), 0.0; atol = 1e-8)
end


text:
using Docile
using Base.Test

include("macro-spec.jl")
include("method-docs.jl")
include("loop-generated-docs.jl")
include("old-tests.jl")

# fix reversed order in query macro
let results = @query doctest(Docile)
    @test length(results) == 1
    @test isa(results[1][1], Method)
    @test isa(results[1][2], Docile.Entry{:method})
end


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = page_model()
    include("../src/components/LGcycle.jl")

    add_comp!(m, LGcycle)

    set_param!(m, :LGcycle, :e_globalLGemissions, readpagedata(m,"test/validationdata/$valdir/e_globalLGemissions.csv"))
    set_param!(m, :LGcycle, :y_year, [2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.,2250.,2300.]) #real value
    set_param!(m, :LGcycle, :y_year_0, 2015.) #real value
    set_param!(m, :LGcycle, :rtl_g_landtemperature, readpagedata(m,"test/validationdata/$valdir/rtl_g_landtemperature.csv"))

    p=load_parameters(m)
    set_leftover_params!(m,p) #important for setting left over component values

    # run Model
    run(m)

    conc=m[:LGcycle,  :c_LGconcentration]
    conc_compare=readpagedata(m,"test/validationdata/$valdir/c_LGconcentration.csv")

    @test conc ≈ conc_compare rtol=1e-4
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
using Jecco.AdS5_3_1

grid = SpecCartGrid3D(
    x_min            = -5.0,
    x_max            =  5.0,
    x_nodes          =  32,
    y_min            = -5.0,
    y_max            =  5.0,
    y_nodes          =  32,
    u_outer_min      =  0.1,
    u_outer_max      =  1.2,
    u_outer_domains  =  2,
    u_outer_nodes    =  32,
    u_inner_nodes    =  12,
    fd_order         =  4,
    sigma_diss       =  0.2,
)

id = BlackBrane_xi1(
    a40    = -1.0,
    AH_pos = 1.0,
    xi_0   = 0.0,
    xi_nx  = 1,
    xi_Ax  = 0.05,
    xi_ny  = 1,
    xi_Ay  = 0.05,
    xmax = grid.x_max,
    xmin = grid.x_min,
    ymax = grid.y_max,
    ymin = grid.y_min,
)

evoleq = AffineNull(
    phi0           = 0.0,
    potential      = ZeroPotential(),
    gaugecondition = Advect_xi(xi_vx = 0.1, xi_vy = 0.1),
)

diag = DiagAH(
    find_AH_every_t    = 0.01,
)

io = InOut(
    out_boundary_every  = 10,
    out_bulk_every      = 200,
    out_gauge_every     = 10,
)

integration = Integration(
    dt              = 0.002,
    tmax            = 10.0,
    ODE_method      = AdS5_3_1.AB3(),
    filter_poststep = true,
)

run_model(grid, id, evoleq, diag, integration, io)


text:
### extirpolation.jl
#
# Copyright (C) 2016 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: periodogram, lomb scargle, extirpolation
#
# This file is a part of LombScargle.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# Utility functions used for implementing the fast algorithm proposed here:
# * Press, W. H., Rybicki, G. B. 1989, ApJ, 338, 277 (URL:
#   http://dx.doi.org/10.1086/167197,
#   Bibcode: http://adsabs.harvard.edu/abs/1989ApJ...338..277P)
# The actual implementation is adapted from Astropy project.  See
# stats/lombscargle/implementations/utils.py in Astropy source code.
#
### Code:

function add_at!{N1<:Number,N2<:Number,N3<:Number}(arr::AbstractVector{N1},
                                                   ind::AbstractVector{N2},
                                                   vals::AbstractVector{N3})
    for i in eachindex(ind)
        arr[mod(ind[i], length(arr)) + 1] += vals[i]
    end
end

function extirpolate{RE<:Real,NU<:Number}(X::AbstractVector{RE},
                                          Y::AbstractVector{NU},
                                          N::Integer=0, M::Integer=4)
    @assert length(X) == length(Y)
    x, y = collect(X), collect(Y)
    if N <= 0
        # Get the maximum of "X", `maximum' has a faster method for ranges.
        N = trunc(Int, maximum(X) + 0.5*M + 1)
    end
    # Now use legendre polynomial weights to populate the results array; This is
    # an efficient recursive implementation (See Press et al. 1989)
    result = zeros(NU, N)
    # first take care of the easy cases where x is an integer
    integers = find(isinteger, x)
    add_at!(result, trunc(Int, x[integers]), y[integers])
    deleteat!(x, integers)
    deleteat!(y, integers)
    # For each remaining x, find the index describing the extirpolation range.
    # i.e. ilo[i] < x[i] < ilo[i] + M with x[i] in the center, adjusted so that
    # the limits are within the range 0...N
    ilo = clamp(trunc(Int, x - div(M, 2)), 0, N - M)
    numerator = y .* [prod(x[j] - ilo[j] - (0:(M - 1))) for j in eachindex(x)]
    denominator = factorial(M - 1)
    @inbounds for j in 0:(M - 1)
        if j > 0
            denominator *= j/(j - M)
        end
        ind = ilo + (M - 1 - j)
        add_at!(result, ind, numerator ./ (denominator * (x .- ind)))
    end
    return result
end

function trig_sum{R1<:Real,R2<:Real}(t::AbstractVector{R1},
                                     h::AbstractVector{R2}, df::Real,
                                     N::Integer, f0::Real=0.0,
                                     freq_factor::Integer=1,
                                     oversampling::Integer=5, Mfft::Integer=4)
    @assert Mfft > 0
    df *= freq_factor
    f0 *= freq_factor
    @assert df > 0
    # `ifft' can take arrays of any length in input, but it's faster when the
    # length is exactly a power of 2.
    Nfft = nextpow2(N * oversampling)
    t0 = minimum(t)
    if f0 > 0
        h = h .* exp(2im * pi * f0 * (t - t0))
    end
    tnorm = mod(((t - t0) * Nfft * df), Nfft)
    grid = extirpolate(tnorm, h, Nfft, Mfft)
    fftgrid = Nfft * ifft(grid)[1:N]
    if t0 != 0
        f = f0 + df * (0:N - 1)
        fftgrid .*= exp(2im * pi * t0 * f)
    end
    C = real(fftgrid)
    S = imag(fftgrid)
    return C, S
end


text:
println("2^62=", 2^62)
println("2^63=", 2^63)


text:
using GeneralizedSampling
using Base.Test

println("Testing miscellaneous functions...")

#=
Testing the miscellaneous functions:
- had!, hadc!, yphad!, conj! should the same result as the standard not in-place functions.
- isuniform for matrices should return true for output from grid and false otherwise.
- In 2D the weights/Voronoi areas are known beforehand for some point configurations.
- For some point configurations, the density is known beforehand.
- Testing the relevant wavelet names for isdaubechies and ishaar.
- Assembling the spilt parts should yield the original array.
=#

# ------------------------------------------------------------
# had!, hadc!, yphad!, conj!

M = rand(2:10)
N = rand(2:10)

# had!
A = rand(M, N)
B = rand(M, N)
hadtmp = A.*B
had!(A, B)
@test hadtmp == A

# hadc!
A = rand(M, N) + im*rand(M,N)
B = rand(M, N) + im*rand(M,N)
hadctmp = A.*conj(B)
hadc!(A, B)
@test hadctmp == A

# yphad!
a = rand(M)
b = rand(M)
y = rand(M)
yphadtmp = y + a.*b
yphad!(y, a, b)
@test yphadtmp == y

# conj!
conj!(A, B)
@test A == conj(B)

# ------------------------------------------------------------
# isuniform and grid

for M in 15:16, scaling in rand(2)
	x = grid(M, scaling)
	@test isuniform(x)

	x[1] += 1
	@test !isuniform(x)

	@test !isuniform(rand(M))

	P = grid( (M,M), scaling )
	@test isuniform(P)

	P[1] += 1
	@test !isuniform(P)

	@test !isuniform(rand(M,2))
end

# ------------------------------------------------------------
# weights

# A corner in the "middle" of each quadrant
xi = [0.5 0.5 ; -0.5 0.5 ; -0.5 -0.5 ; 0.5 -0.5]
@test_approx_eq weights(xi, 1.0) ones(4)

# ------------------------------------------------------------
# Density

# 1D: Equidistant points
for N in 3:10
	x = linspace(-1,1,N)
	@test_approx_eq density(x,1) 2/(N-1)
end

@test density( [-1.0; 1.0], 1 ) == 2.0

# A corner in the "middle" of each quadrant
xi = [0.5 0.5 ; -0.5 0.5 ; -0.5 -0.5 ; 0.5 -0.5]
@test_approx_eq density(xi, 1.0) sqrt(2)/2

# ------------------------------------------------------------
# ishaar and isdaubechies

for N in 2:10, prefix in ["db", "DB", "Db", "dB"]
	name = string(prefix, N)
	@test isdaubechies(name)
	@test !ishaar(name)
end

@test ishaar("haar")
@test ishaar("db1")

@test !isdaubechies("ddb2")

# ------------------------------------------------------------
# Split

# Vector
x = rand( rand(5:10) )
L, I, R = split(x,2)
@test x == vcat(L, I, R)

# Matrix
A = rand(rand(5:10), rand(5:10))
S = split(A,2)
@test A == [ S.LL S.LI S.LR ; S.IL S.II S.IR ; S.RL S.RI S.RR ]


text:
export Sha1, Oid, raw, iszero, @sha1_str

const OID_RAWSZ = 20
const OID_HEXSZ = OID_RAWSZ * 2
const OID_MINPREFIXLEN = 4

# immutable Oid
#   id1::Uint8
#   id2::Uint8
#   ...
#   id20::Uint8
# end

@eval begin
    $(Expr(:type, false, :Oid,
        Expr(:block,
            [Expr(:(::), symbol("id$i"), :Uint8) for i=1:OID_RAWSZ]...)))
end

# default Oid constructor (all zeros)
Oid() = @eval begin
    $(Expr(:call, :Oid, [:(zero(Uint8)) for _=1:OID_RAWSZ]...))
end

Oid(ptr::Ptr{Oid}) = unsafe_load(ptr)::Oid

Oid(ptr::Ptr{Uint8}) = begin
    if ptr == C_NULL
        throw(ArgumentError("NULL pointer passed to Oid() constructor"))
    end
    oid_ptr = Oid[Oid()]
    ccall((:git_oid_fromraw, libgit2), Void, (Ptr{Oid}, Ptr{Uint8}), oid_ptr, ptr)
    return oid_ptr[1]
end

Oid(id::Array{Uint8,1}) = begin
    if length(id) != OID_RAWSZ
        throw(ArgumentError("invalid raw buffer size"))
    end
    return Oid(pointer(id))
end

Oid(id::String) = begin
    bstr = bytestring(id)
    if sizeof(bstr) != OID_HEXSZ
        throw(ArgumentError("invalid hex size"))
    end
    oid_ptr = Oid[Oid()]
    @check ccall((:git_oid_fromstrp, libgit2), Cint,
                 (Ptr{Oid}, Ptr{Cchar}), oid_ptr, bstr)
    return oid_ptr[1]
end

Oid(id::Oid) = id

Base.copy!(arr::Array{Uint8, 1}, id::Oid) = begin
    unsafe_copy!(convert(Ptr{Uint8}, arr), convert(Ptr{Uint8},  hex2bytes(hex(id))), OID_RAWSZ)
    return arr
end

raw(id::Oid) = copy!(Array(Uint8, OID_RAWSZ), id)

const _hexstr = Array(Uint8, OID_HEXSZ)
Base.hex(id::Oid) = begin
    ccall((:git_oid_nfmt, libgit2), Void,
          (Ptr{Uint8}, Csize_t, Ptr{Oid}), _hexstr, OID_HEXSZ, &id)
    return bytestring(convert(Ptr{Uint8}, _hexstr), OID_HEXSZ)
end

Base.string(id::Oid) = hex(id)

Base.show(io::IO, id::Oid) = print(io, "Oid($(string(id)))")

Base.hash(id::Oid) = hash(hex(id))

Base.cmp(id1::Oid, id2::Oid) = int(ccall((:git_oid_cmp, libgit2), Cint,
                                         (Ptr{Oid}, Ptr{Oid}), &id1, &id2))

Base.(:(==))(id1::Oid, id2::Oid) = cmp(id1, id2) == 0
Base.isequal(id1::Oid, id2::Oid) = cmp(id1, id2) == 0
Base.isless(id1::Oid, id2::Oid)  = cmp(id1, id2) < 0

iszero(id::Oid) = begin
    bytes = raw(id)
    for i=1:OID_RAWSZ
        if bytes[i] != zero(Uint8)
            return false
        end
    end
    return true
end

immutable Sha1
    data::ASCIIString

    Sha1(s::String) = begin
        bstr = bytestring(s)::ASCIIString
        if sizeof(bstr) != OID_HEXSZ
            throw(ArgumentError("invalid sha1 string length $(length(bstr))"))
        end
        for c in bstr
            if !('0' <= c <= '9' || 'a' <= c <= 'f' || 'A' <= c <= 'F')
                throw(ArgumentError("not a hexadecimal string: $(repr(s))"))
            end
        end
        return new(bstr)
    end
end

macro sha1_str(s)
    Sha1(s)
end

Oid(sha::Sha1) = Oid(hex2bytes(sha.data))


text:
function eosshow(io::IO, ::MIME"text/plain", Base.@nospecialize(model::EoSModel))
    print(io, typeof(model))
    length(model) == 1 && println(io, " with 1 component:")
    length(model) > 1 && println(io, " with ", length(model), " components:")
    for i in 1:length(model)
        print(io, " \"", model.components[i], "\"")
        println(io)
    end
    if hasfield(typeof(model),:params)
        paramnames = fieldnames(typeof(model.params))
        len_params = length(paramnames)
        !iszero(len_params) && print(io,"Contains parameters: ")
        firstloop = true
        for fieldname in paramnames
            firstloop == false && print(io, ", ")
            print(io, fieldname)
            firstloop = false
        end
    end
end

function eosshow(io::IO, Base.@nospecialize(model::EoSModel))
    print(io, typeof(model))
    firstloop = true
    print(io, "(")
    for i in 1:length(model.components)
        firstloop == false && print(io, ", ")
        print(io, "\"", model.components[i], "\"")
        firstloop = false
    end
    print(io, ")")
end

function gc_eosshow(io::IO, ::MIME"text/plain", Base.@nospecialize(model::EoSModel))
    print(io, typeof(model))
    length(model) == 1 && println(io, " with 1 component:")
    length(model) > 1 && println(io, " with ", length(model), " components:")
    param = model.groups
    for i in 1:length(param.components)

        print(io, " \"", param.components[i], "\": ")
        firstloop = true
        for j in 1:length(param.n_groups[i])
            firstloop == false && print(io, ", ")
            print(io, "\"", param.groups[i][j], "\" => ", param.n_groups[i][j])
            firstloop = false
        end
        i != length(param.components) && println(io)
    end
    println(io)
    print(io, "Contains parameters: ")
    firstloop = true
    for fieldname in fieldnames(typeof(model.params))
        firstloop == false && print(io, ", ")
        print(io, fieldname)
        firstloop = false
    end
end

export eosshow


text:
module R_constraints
    include(dirname(@__FILE__)*"\\Utils_cons.jl")
    include(dirname(dirname(@__FILE__))*"\\Moments\\Moments.jl")
    using .Utils_cons
    using .Moments
    const uc = Utils_cons
    const mom = Moments

    export make_mon_expo_keys,
           make_PSD_con,
           make_ord4_con,
           make_loc_cons_S_inf,
           make_loc_cons_S₂,
           make_loc_cons_S₂₁,
           make_G_con

    make_mon_expo_keys(n,t::Int) = Moments.make_mon_expo(n,t*2)

    """L([x,y]ₜ[x,y]ₜᵀ) ⪰ 0"""
    function make_PSD_con(d,t,Lx)
        MB = mom.get_ℝ_block_diag(d,t)
        Bs = [keys(MB)...]
        PSD_con = Dict(zip(Bs, map(b -> uc.idx2var(Lx, MB[b]),Bs)))
        return PSD_con
    end

    # """ρ = aaᵀ⊗bbᵀ , ρᵢⱼᵢⱼ = 0 ⟹ aᵢbⱼ = 0 ⟹ L(xᵢyⱼ) = 0 """
    #
    # function zeroprop(ρ,d,t,Lx)
    #     nair((i,j)) = vcat(mom.eᵢ(d[1],i),mom.eᵢ(d[2],j))
    #     zdiags =[(i,j) for i in 1:d[1], j in 1:d[2] if ρ[i+(j-1)*d[1],i+(j-1)*d[1]] == 0 ]
    #     return uc.idx2var(Lx, [k .* m  for m ∈  map(nair, zdiags) for k in vcat(1,3:t[1])])
    # end

    """L(xxᵀ⊗yyᵀ)"""
    make_ord4_con(d,Lx) = uc.idx2var(Lx,mom.make_xxᵀ_tens_yyᵀ(d))

    """ L(g⋅η) ⪰ 0 for η ∈ even-degree-principle-submatrices of [x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ
                        g = √(ρₘₐₓ) - xᵢ² , √(ρₘₐₓ) - yᵢ² and i ∈ [d] """
    function make_loc_con(eᵢ,sqr_ρ,B,Lx)
        Lxᵢ = uc.idx2var(Lx, B + repeat([2*eᵢ], size(B)...)) #  L(xᵢ²⋅η)
        Lρ = sqr_ρ*uc.idx2var(Lx, B)  # L(√ρₘₐₓ⋅η)
        return  Lρ - Lxᵢ
    end

    function make_loc_cons_S_inf(ρ,d,t,Lx)
        n     = sum(d)
        MB    = mom.get_ℝ_block_diag(d,t.- 1)
        sqr_ρ = sqrt(maximum(ρ))

        loc_con = Dict()
        for b in keys(MB)
            for k in 1:d[1]
                loc_con[(b,"x²_$k")] = make_loc_con(uc.eᵢ(n,k),sqr_ρ,MB[b],Lx) # Constraint:  L((√(ρₘₐₓ) - xᵢ²))⋅η) ⪰ 0
            end
            for k in 1:d[2]
                loc_con[(b,"y²_$k")] = make_loc_con(uc.eᵢ(n,k+d[1]),sqr_ρ,MB[b],Lx) # Constraint:  L((√(ρₘₐₓ) - yᵢ²))⋅η) ⪰ 0
            end
        end
        return loc_con
    end

    """L(g⋅η) ⪰ 0 ∀ η ∈ even-degree-principle-submatrices of [x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ
                       g = √(Tr(ρ)) - ∑xᵢ² , √(Tr(ρ)) - ∑yᵢ² """
    function make_loc_cons_S₂(ρ,d,t,Lx)
        d₁,d₂ = d
        n     = sum(d)
        MB    = mom.get_ℝ_block_diag(d,t.- 1)
        loc_con = Dict() ; g₂ = Dict()
        stρ   = sqrt(sum([ρ[i,i] for i in 1:n]))
        for b in keys(MB)
            xRterm = sum([uc.idx2var(Lx, MB[b] .+ 2* [uc.eᵢ(n,k)])    for k in 1:d₁ ] ) # L((∑xᵢ²)⋅[x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ )
            yRterm = sum([uc.idx2var(Lx, MB[b] .+ 2* [uc.eᵢ(n,k+d₁)]) for k in 1:d₂ ]) # L((∑yᵢ²)⋅[x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ )
            loc_con[b,"x"] = stρ*uc.idx2var(Lx, MB[b]) - xRterm # √(Tr(ρ))⋅L([x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ) - L((∑xᵢ²)⋅[x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ )
            loc_con[b,"y"] = stρ*uc.idx2var(Lx, MB[b]) - yRterm # √(Tr(ρ))⋅L([x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ) - L((∑yᵢ²)⋅[x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ )
            g₂[b] = xRterm - yRterm
        end
        return loc_con, g₂
    end

    """ L((Tr(ρ) - ∑xᵢ²)⋅η) ⪰ 0 ∀ η ∈ "even-degree"-principle-submatrices of [x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ
        L((∑yᵢ² - 1)⋅η) = 0 """
    function make_loc_cons_S₂₁(ρ,d,t,Lx)
        d₁,d₂ = d
        n     = sum(d)
        MB    = mom.get_ℝ_block_diag(d,t.- 1)
        trρ   = sum([ρ[i,i] for i in 1:size(ρ)[1]])
        loc_con = Dict()
        g₂ = Dict()
        for b in keys(MB)
            MBLx   = uc.idx2var(Lx, MB[b]) # L([x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ)
            xRterm = sum([uc.idx2var(Lx, MB[b] .+ [2*uc.eᵢ(n,k)]) for k in 1:d₁ ]) # L((∑xᵢ²)⋅[x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ )
            yRterm = sum([uc.idx2var(Lx, MB[b] .+ [2*uc.eᵢ(n,k+d₁)]) for k in 1:d₂ ]) # L((∑yᵢ²)⋅[x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ )
            loc_con[b] = MBLx - xRterm #  tr(ρ)L([x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ) - L((∑xᵢ²)⋅[x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ ) ⪰ 0
            g₂[b]      = MBLx - yRterm # tr(ρ)L([x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ) - L((∑yᵢ²)⋅[x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ ) = 0
        end
        return loc_con, g₂
    end

    """ ρ⊗L(η) - L( (xxᵀ⊗yyᵀ) ⊗ η ) ⪰ 0 ∀ η ∈ even-degree-principle-submatrices of ([x,y]₌ₜ₋₂[x,y]₌ₜ₋₂ᵀ) """
    function make_G_con(ρ,d,t,Lx)
        MB          = mom.get_ℝ_block_diag(d,t.- 2)  # exponents of [x, y]ₜ₋₂[x, y]ᵀₜ₋₂
        xxᵀtensyyᵀ  = mom.make_xxᵀ_tens_yyᵀ(d)       # exponents of xxᵀ⊗yyᵀ

        G_con = Dict()
        for b in keys(MB)
            Lterm       = kron(ρ, uc.idx2var(Lx,MB[b]))  # ρ⊗L([x, y]ₜ₋₂[x, y]ᵀₜ₋₂)
            Rterm       = uc.idx2var(Lx, mom.var_kron(xxᵀtensyyᵀ,MB[b]))  # L( (xxᵀ⊗yyᵀ) ⊗ ([x, y]ₜ₋₂[x, y]ᵀₜ₋₂) )

            G_con[b]    = Lterm - Rterm # ρ⊗L([x,y]ₜ₋₂[x,y]ᵀₜ₋₂) - L( (xxᵀ⊗yyᵀ) ⊗ ([x,y]ₜ₋₂[x,y]ᵀₜ₋₂) )
        end
        return G_con
    end
end

max_stars_repo_name: JAndriesJ/sep-rank-bounding

max_issues_repo_name: JAndriesJ/sep-rank-bounding

max_forks_repo_name: JAndriesJ/sep-rank-bounding


text:
using MittagLefflerFunctions
using Polynomials
using PyPlot
using SpecialFunctions: erfcx

α = 1/2
β = 1.0
N = 5
r = 2N + 1
m = r + 1
n = r

#Eref = MLQuad(α, β, 2N, :hyperbola)
Eref(z) = erfcx(-z)

EP = MLPade(α, β, m, n, :svd)
p, q = Polynomial(EP.p), Polynomial(EP.q)
χ = roots(q)
EPade(z) = p(-z) / q(-z)

EQ = MLQuad(α, β, N, :hyperbola)
w, C, A = EQ.qs.w, EQ.qs.C, EQ.qs.A
wα = w.^α
P = [ conj.(wα[N:-1:1]); wα[0:N] ]
R1 = -A * C .* w.^(α-β)
R = [ conj.(R1[N:-1:1]); R1[0:N] ]

function mlf_quad1(z)
    s = 0.0
    for j = 1:length(P)
        s += R[j] / ( z - P[j] )
    end
    return s
end

figure(1)
plot(real.(-χ), imag.(-χ), "o",
     real.(P), imag.(P), "o")
legend((L"$-\chi_j$", L"$P_j$"))
grid(true)

Nx = 301
Ny = 301
xmin, xmax = -5.0, 3.0
ymin, ymax = -4.0, 4.0
x = range(xmin, xmax, length=Nx)
y = range(ymin, ymax, length=Ny)

z = Complex{Float64}[ complex(x[j], y[i]) for i=1:Ny, j=1:Nx ]

err_quad = Float64[ abs(Eref(z[i,j])-mlf_quad1(z[i,j])) for i = 1:Ny, j=1:Nx ]
err_pade = Float64[ abs(Eref(z[i,j])-EPade(z[i,j])) for i = 1:Ny, j=1:Nx ]

figure(2)
contour(x, y, log10.(err_quad), levels=collect(-9:4))
plot(real.(P), imag.(P), "ro")
colorbar()
xlabel(L"$\Re z$")
ylabel(L"$\Im z$")
grid(true)
savefig("figure11L.pdf")

figure(3)
contour(x, y, log10.(err_pade), levels=collect(-14:4))
plot(real.(-χ), imag.(-χ), "ro")
colorbar()
xlabel(L"$\Re z$")
ylabel(L"$\Im z$")
grid(true)
savefig("figure11R.pdf")


text:
using RunTests
using Base.Test

@testmodule SeventhSetOfTests begin
  function test_hello()
    @test true
  end
end


text:
export Model
export write_to_file!, predict, loss_function, update_model!, update_model!_2

import .ArcEager: transitions_number, set_labels!, execute_transition, zero_cost_transitions, transition_costs, gold_scores, is_transition_valid
import .ArcEager: GoldState, FORBIDDEN_COST

using Flux, CUDA
using Zygote
using JLD2
using DependencyParserTest

#=
  Input layer dimension: batch_size * embeddings_size
  hidden layer weights dimension(h): hidden_size * (batch_size * embeddings_size)
  bias dineansion: hidden_size
  output layer weights dimensions: labels_num * hidden_size
=#

# CUDA.allowscalar(false)

mutable struct Model
  embeddings
  hidden_layer::Dense
  output_layer::Dense
  word_ids::Dict{String, Integer}
  tag_ids::Dict{String, Integer}
  label_ids::Dict{String, Integer}
  gpu_available
  hidden_accumulator

  function Model(
    embeddings,
    hidden_layer::Dense,
    output_layer::Dense,
    word_ids::Dict{String, Integer},
    tag_ids::Dict{String, Integer},
    label_ids::Dict{String, Integer},
  )
    model = new(embeddings, hidden_layer, output_layer, word_ids, tag_ids, label_ids)
    model.gpu_available = CUDA.functional()
    model
  end

  function Model(settings::Settings, system::ParsingSystem, embeddings_file::String, connlu_sentences::Vector{ConnluSentence})
    model = new()

    loaded_embeddings, embedding_ids = read_embeddings_file(embeddings_file)
    embeddings_size = length(loaded_embeddings[begin, :])
    if embeddings_size != settings.embeddings_size
      ArgumentError("Incorrect embeddings dimensions. Given: $(embeddings_size). In settings: $(settings.embeddings_size)") |> throw
    end

    set_corpus_data!(model, connlu_sentences)
    set_labels!(system, model.label_ids |> keys |> collect)

    model.embeddings = rand(Float64, length(model.word_ids) + length(model.tag_ids) + length(model.label_ids), embeddings_size)
    match_embeddings!(model.embeddings, loaded_embeddings, embedding_ids, model.word_ids |> keys |> collect)

    model.hidden_layer = Dense(settings.batch_size * embeddings_size, settings.hidden_size)
    model.output_layer = Dense(settings.hidden_size, transitions_number(system), bias=false, sigmoid)
    model.hidden_accumulator = zeros(Float32, settings.hidden_size)

    model.gpu_available = CUDA.functional()

    model
  end

  function Model(model_file::String)
    lines = readlines(model_file)
    info_line = lines[begin]
    words_count, tags_count, labels_count, embeddings_size = split(info_line, " ") |> info_values -> map(value -> parse(Int32, value), info_values)
    total_ids_count =  words_count + tags_count + labels_count

    word_ids = Dict{String, Integer}()
    tag_ids = Dict{String, Integer}()
    label_ids = Dict{String, Integer}()
    embeddings = Matrix{Float64}(undef, total_ids_count, embeddings_size)

    # read all embeddings for words, tags and labels
    for i = 1:total_ids_count
      line = lines[i + 1]
      entity, embedding = split(line) |> values -> [values[begin], map(value -> parse(Float64, value), values[begin+1:end])]
      embeddings[i, :] = embedding

      # Depend on index define entity as word or as tag or as label
      if 1 <= i <= words_count
        word_ids[entity] = i
      elseif words_count < i <= words_count + tags_count
        tag_ids[entity] = i
      else
        label_ids[entity] = i
      end
    end

    info_line = lines[total_ids_count + 2]
    batch_size, hidden_size, labels_num = split(info_line) |> info_values -> map(value -> parse(Int32, value), info_values)

    hidden_layer = Matrix{Float64}(undef, hidden_size, batch_size * embeddings_size)
    output_layer = Matrix{Float64}(undef, labels_num, hidden_size)

    # read hidden layer weights and bias
    for i = 1:hidden_size
      line = lines[i + total_ids_count + 2]
      hidden_layer[i, :] = split(line) |> values -> map(value -> parse(Float64, value), values)
    end
    hidden_bias_line = lines[begin + total_ids_count + hidden_size + 2]
    hidden_bias = split(hidden_bias_line) |> values -> map(value -> parse(Float64, value), values)

    # read softmax layer weights
    for i = 1:labels_num
      line = lines[i + total_ids_count + hidden_size + 3]
      output_layer[i, :] = split(line) |> values -> map(value -> parse(Float64, value), values)
    end

    hidden_accumulator = zeros(Float32, hidden_size)

    model = new(embeddings, Dense(hidden_layer, hidden_bias), Dense(output_layer, false, sigmoid), word_ids, tag_ids, label_ids, undef, hidden_accumulator)
    model.gpu_available = CUDA.functional()

    enable_cuda(model)

    model
  end
end

Flux.trainable(m::Model) = (m.embeddings, m.hidden_layer.weight, m.hidden_layer.bias, m.output_layer.weight,)

function enable_cuda(model::Model)
  if model.gpu_available
    @info "enable cuda"
    model.embeddings = cu(model.embeddings)
    model.hidden_layer = fmap(cu, model.hidden_layer)
    model.output_layer = fmap(cu, model.output_layer)
    model.hidden_accumulator = cu(model.hidden_accumulator)
  end
end

function calculate_hidden(model, input)
  fill!(model.hidden_accumulator, zero(Float32))

  embeddings_size = length(model.embeddings[begin, :])
  batch_size = length(input)
  hidden_weight = model.hidden_layer.weight

  for i in 1:batch_size
    offset = (i - 1) * embeddings_size
    hidden_slice = view(hidden_weight, :, (offset + 1) : (offset + embeddings_size))

    model.hidden_accumulator .+= hidden_slice * input[i]
  end

  model.hidden_accumulator .+= model.hidden_layer.bias

  model.hidden_accumulator
end

# some of operators in calculate_hidden not compilable for Zygote
function calculate_hidden_train(model, input)
  result = zeros(Float32, length(model.hidden_layer.weight[:, begin]))
  if model.gpu_available
    result = cu(result)
  end

  embeddings_size = length(model.embeddings[begin, :])
  batch_size = length(input)
  hidden_weight = model.hidden_layer.weight

  for i in 1:batch_size
    offset = (i - 1) * embeddings_size
    hidden_slice = view(hidden_weight, :, (offset + 1) : (offset + embeddings_size))

    result += hidden_slice * input[i]
  end

  result += model.hidden_layer.bias

  Flux.dropout(result, 0.5, dims = 1)
end

function predict(model::Model, batch)
  (calculate_hidden(model, batch) .^ 3) |>
    model.output_layer |>
    result -> reshape(result, 1, :)
end

function train_predict_tree(model::Model, sentence::Sentence, context::TrainingContext)
  config = Configuration(sentence)
  transitions_number = 0

  while !is_terminal(config)
    transition = predict_transition(model, context.settings, context.system, config)
    transition === nothing && break
    execute_transition(config, transition, context.system)

    transitions_number += 1
    transitions_number > LIMIT_TRANSITIONS_NUMBER && break
  end

  config.tree
end

function predict_train(model::Model, batch)
  take_batch_embeddings(model, batch) |>
    batch_emb ->(calculate_hidden_train(model, batch_emb) .^ 3) |>
    Flux.normalise |>
    model.output_layer |>
    result -> reshape(result, 1, :)
end

function predict_transition(model::Model, settings::Settings, system::ParsingSystem, config::Configuration)
  scores = form_batch(model, settings, config) |> batch -> take_batch_embeddings(model, batch) |>
    batch -> predict(model, batch) |> collect |>
    result -> filter(score -> score[end] > 0 && is_transition_valid(config, system.transitions[score[begin]], system), collect(enumerate(result)))

  length(scores) == 0 && return

  score_index = findmax(map(score -> score[end], scores))[end]

  system.transitions[scores[score_index][begin]]
end

function update_model!(model::Model, dataset, training_context::TrainingContext)
  ps = Flux.params(model)

  loss = (sample) -> begin
    sum(sample) do (batch, gold)
      predict_train(model, batch) |> scores -> transition_loss(scores, gold)
    end + L2_norm(ps, training_context.settings)
  end

  test_sample = first(dataset)

  evalcb = () -> begin
   @show test_loss(model, test_sample, training_context)
   CUDA.reclaim()
  end

  throttle_cb = Flux.throttle(evalcb, 10)

  Flux.train!(loss, ps, dataset, training_context.optimizer, cb=throttle_cb)
end

function loss_function(entropy_sum, ps, settings::Settings)
  entropy_sum + L2_norm(ps, settings)
end

function L2_norm(ps, settings::Settings)
  sqnorm(x) = sum(abs2, x)

  sum(sqnorm, ps) * (settings.reg_weight / 2)
end

function transition_loss(scores, gold)
  Flux.Losses.binary_focal_loss(scores, gold, γ=2)
end

function test_transition_loss(scores, gold)
  Flux.Losses.binarycrossentropy(scores, gold)
end

function train!(model::Model, training_context::TrainingContext)
  training_context.optimizer = ADAM()

  training_context.test_dataset = Flux.DataLoader(build_dataset(model, training_context.test_connlu_sentences[begin:begin+800], training_context), batchsize=training_context.settings.sample_size, shuffle=true)

  raw_dataset = build_dataset(model, training_context.connlu_sentences, training_context)
  train_datasets = Flux.DataLoader(raw_dataset, batchsize=Int64(round(length(raw_dataset) / 10)), shuffle=true)

  train_epoch = () -> begin
    for dataset in train_datasets
      samples = Flux.DataLoader(dataset, batchsize=training_context.settings.sample_size, shuffle=true)

      update_model!(model, samples, training_context)

      GC.gc()
      CUDA.reclaim()

      @info "Dataset ends, start test"
      test_training_scores(model, training_context)

      GC.gc()
      CUDA.reclaim()
    end
  end

  @info "Start training"
  Flux.@epochs 500 train_epoch()
end

function build_dataset(model, sentences_batch, training_context)
  train_samples = []
  frequency = Dict()
  mutex = ReentrantLock()
  sent_mutex = ReentrantLock()
  freq_mutex = ReentrantLock()
  sentences_number = 0

  parse_gold_tree = (connlu_sentence) -> begin
    batch = []

    sentence = zip(connlu_sentence.token_doc.tokens, connlu_sentence.pos_tags) |> collect |> Sentence

    config = Configuration(sentence)
    process_config(model, config, connlu_sentence.gold_tree, training_context, 0, mutex, batch)

    for (_, gold) in batch
      gold_correct = findall(e -> e == 1, gold)

      lock(freq_mutex) do
        if haskey(frequency, gold_correct)
          frequency[gold_correct] += 1
        else
          frequency[gold_correct] = 1
        end
      end
    end

    lock(sent_mutex) do
      train_samples = vcat(train_samples, batch)
      sentences_number += 1
      if sentences_number % 250 == 0
        @info "Sentences processed: $sentences_number"
      end
    end
  end

  try
    Threads.@sync begin
      @info "Build dataset..."
      for connlu_sentence in sentences_batch
        Threads.@spawn parse_gold_tree(connlu_sentence)
      end
    end
  catch err
    println(err.task.exception)
    return err
  end

  @info "Do some upsample...."

  train_samples = map(train_samples) do sample
    gold = sample[end]
    gold_correct = findall(e -> e != 0, gold)

    repeat_number = if frequency[gold_correct] < 300
      7
    elseif frequency[gold_correct] < 700
      5
    elseif frequency[gold_correct] < 1000
      3
    else
      1
    end

    frequency[gold_correct] += repeat_number - 1

    if model.gpu_available
      gold = cu(gold)
    end

    [(sample[begin], gold) for i in 1:repeat_number]
  end |> Iterators.flatten |> collect

  @info frequency
  @info "Dataset builded..."

  train_samples
end

function process_config(model, config, gold_tree, context, transition_number, mutex, train_samples)
  is_terminal(config) && return
  transition_number >= LIMIT_TRANSITIONS_NUMBER && return

  gold_state = GoldState(gold_tree, config, context.system)
  gold  = transition_costs(gold_state) |> gold_scores

  zero_transitions = zero_cost_transitions(gold_state)
  length(zero_transitions) == 0 && return

  batch = form_batch(model, context.settings, config)

  lock(mutex) do
    push!(train_samples, (batch, gold))
  end

  if rand(1:100) in 1:(context.beam_coef * 100)
    foreach(zero_transitions) do transition
      next_config = deepcopy(config)
      execute_transition(next_config, transition, context.system)

      process_config(model, next_config, gold_tree, context, transition_number + 1, mutex, train_samples)
    end
  else
    transition = rand(zero_transitions)
    execute_transition(config, transition, context.system)

    process_config(model, config, gold_tree, context, transition_number + 1, mutex, train_samples)
  end
end

function test_training_scores(model::Model, context::TrainingContext)
  losses = []
  parsed_trees_file = "tmp/parsed_trees.txt"
  mutex = ReentrantLock()
  losses_mutex = ReentrantLock()
  parse_tree = (connlu_sentence, file) -> begin
    sentence = zip(connlu_sentence.token_doc.tokens, connlu_sentence.pos_tags) |> collect |> Sentence

    tree = train_predict_tree(model, sentence, context)
    tree_text = convert_to_string(tree)

    lock(mutex) do
      write(file, tree_text)
      write(file, "\n")
      write(file, "\n")
    end
  end

  parse_sample = (sample) -> begin
    loss = test_loss(model, sample, context)

    lock(losses_mutex) do
      push!(losses, loss)
    end
  end

  open(parsed_trees_file, "w") do file
    for i in 1:length(context.test_connlu_sentences)
      connlu_sentence = context.test_connlu_sentences[i]

      parse_tree(connlu_sentence, file)
    end
  end

  test_dataset = context.test_dataset |> collect

  for i in 1:length(test_dataset)
    sample = test_dataset[i]

    parse_sample(sample)
  end

  conllu_source = DependencyParserTest.Sources.ConlluSource(context.test_connlu_file)
  parser_source = DependencyParserTest.Sources.CoreNLPSource(parsed_trees_file)

  uas, las = DependencyParserTest.Benchmark.test_accuracy(conllu_source, parser_source)
  best_loss = min(losses...)
  worst_loss = max(losses...)
  avg_loss = sum(losses) / length(losses)

  open(context.training_results_file, "a") do file
    result_line = "UAS=$(uas), LAS=$(las), best_loss=$(best_loss), worst_loss=$(worst_loss), avg_loss=$(avg_loss)\n"
    write(file, result_line)
  end

  write_to_file!(model, context.model_file * "_last.txt")

  if uas > context.best_uas || (uas == context.best_uas && las > context.best_las)
    context.best_uas = uas
    context.best_las = las
    context.best_loss = avg_loss

    write_to_file!(model, context.model_file * "_best.txt")
  end
end

function test_loss(model::Model, sample, context)
  ps = Flux.params(model)

  sum(sample) do (batch, gold)
    take_batch_embeddings(model, batch) |>
    batch_emb -> predict(model, batch_emb) |>
      scores -> test_transition_loss(scores, gold)
  end + L2_norm(ps, context.settings)
end

#=
File structure
known_words_count known_pos_tags_count known_labels_count embeddings_size
words embeddings
pos_tags embeddings
labels embeddings
batch_size hidden_size label_nums
hidden weights
hidden bias
softmax weights
=#
function write_to_file!(model::Model, filename::String)
  embeddings = model.embeddings |> collect
  hidden_weight = model.hidden_layer.weight |> collect
  hidden_bias = model.hidden_layer.bias |> collect
  output_weight = model.output_layer.weight |> collect

  open(filename, "w") do file
    embeddings_size = length(embeddings[begin, :])

    sort_by_value(pair) = pair[end]
    write(file, "$(length(model.word_ids)) $(length(model.tag_ids)) $(length(model.label_ids)) $(embeddings_size)\n")
    known_entities = vcat(
      sort(collect(model.word_ids), by = sort_by_value),
      sort(collect(model.tag_ids), by = sort_by_value),
      sort(collect(model.label_ids), by = sort_by_value)
    )

    foreach(known_entities) do (entity, entity_id)
      embedding = embeddings[entity_id, :]
      write(file, "$(entity) $(join(embedding, " "))\n")
    end

    hidden_size = length(hidden_weight[:, begin])
    labels_num = length(output_weight[:, begin])
    batch_size = Int32(length(hidden_weight[begin, :]) / embeddings_size)

    write(file, "$(batch_size) $(hidden_size) $(labels_num)\n")

    for i = 1:hidden_size
      write(file, join(hidden_weight[i, :], " "))
      write(file, "\n")
    end

    write(file, join(hidden_bias, " "))
    write(file, "\n")

    for i = 1:labels_num
      write(file, join(output_weight[i, :], " "))
      write(file, "\n")
    end
  end
end

function read_embeddings_file(filename::String)
  lines = readlines(filename)
  words_count, dimension = split(lines[begin]) |> (numbers -> map(number -> parse(Int64, number), numbers))
  deleteat!(lines, 1)

  embeddings = zeros(Float64, words_count, dimension)
  embedding_ids = Dict{String, Integer}()

  foreach(enumerate(lines)) do (index, line)
    splitted = split(line, " ")
    word = line[begin]
    embedding_ids[string(word)] = index

    for i = 1:dimension
      embeddings[index, i] = parse(Float64, splitted[i + 1])
    end
  end

  [embeddings, embedding_ids]
end

function match_embeddings!(embeddings::Matrix{Float64}, loaded_embeddings::Matrix{Float64}, embedding_ids::Dict{String, Integer}, known_words::Vector{String})
  embeddings_size = size(embeddings)[end]

  foreach(enumerate(known_words)) do (index, word)
    embedding_id = if haskey(embedding_ids, word)
      embedding_ids[word]
    elseif haskey(embedding_ids, lowercase(word))
      embedding_ids[lowercase(word)]
    else
      0
    end

    if embedding_id > 0
      embeddings[index, :] = loaded_embeddings[embedding_id, :]
    else
      embeddings[index, :] = [rand(-100:100) / 10000.0 for i in 1:embeddings_size]
    end
  end
end

function set_corpus_data!(model::Model, connlu_sentences::Vector{ConnluSentence})
  seq_id = Sequence()

  model.word_ids = Dict{String, Integer}()
  model.tag_ids = Dict{String, Integer}()
  model.label_ids = Dict{String, Integer}()

  corpus = map(conllu_sentence -> conllu_sentence.token_doc, connlu_sentences) |> Corpus
  update_lexicon!(corpus)

  lexicon(corpus) |> keys |> collect |> words -> foreach(word -> model.word_ids[word] = next!(seq_id), words)
  model.word_ids[UNKNOWN_TOKEN]= next!(seq_id)
  model.word_ids[NULL_TOKEN]= next!(seq_id)
  model.word_ids[ROOT_TOKEN]= next!(seq_id)

  map(conllu_sentence -> conllu_sentence.pos_tags, connlu_sentences) |>
    Iterators.flatten |>
    collect |>
    unique |>
    tags -> foreach(tag -> model.tag_ids[tag] = next!(seq_id), tags)
  model.tag_ids[UNKNOWN_TOKEN]= next!(seq_id)
  model.tag_ids[NULL_TOKEN]= next!(seq_id)
  model.tag_ids[ROOT_TOKEN]= next!(seq_id)

  map(conllu_sentence -> conllu_sentence.gold_tree.nodes, connlu_sentences) |>
    Iterators.flatten |>
    collect |>
    nodes -> map(node -> node.label, nodes) |>
    unique |>
    labels -> foreach(label -> model.label_ids[label] = next!(seq_id), labels)
  model.label_ids[NULL_TOKEN]= next!(seq_id)
end

#=
while batch_size only is 48 there is structure of batch
1-12 - word_ids
13-24 - tag_ids
25-32 - label_ids
=#
const POS_OFFSET = 12
const LABEL_OFFSET = 24
const STACK_OFFSET = 4
const STACK_NUMBER = 4

function form_batch(model::Model, settings::Settings, config::Configuration)
  batch = zeros(Integer, settings.batch_size)

  word_id_by_word_index(word_index::Integer) = get_token(config, word_index) |> token -> get_word_id(model, token)
  tag_id_by_word_index(word_index::Integer) = get_tag(config, word_index) |> tag -> get_tag_id(model, tag)
  label_id_by_word_index(word_index::Integer) = get_label(config, word_index) |> label -> get_label_id(model, label)

  # add top two stack elements and top two buffers elems with their's tags
  for i = 1:2
    stack_word_index = get_stack_element(config, i)
    buffer_word_index = get_buffer_element(config, i)

    batch[i] = word_id_by_word_index(stack_word_index)
    batch[i + POS_OFFSET] = tag_id_by_word_index(stack_word_index)
    batch[i + 2] = word_id_by_word_index(buffer_word_index)
    batch[i + POS_OFFSET + 2] = tag_id_by_word_index(buffer_word_index)
  end

  #=
    Add:
    The first and second leftmost / rightmost children of the first stack and buffer element
  =#

  set_word_data_by_index_with_offset = function (word_index::Integer, additional_offset)
    batch[STACK_OFFSET + additional_offset] = word_id_by_word_index(word_index)
    batch[STACK_OFFSET + POS_OFFSET + additional_offset] = tag_id_by_word_index(word_index)
    batch[LABEL_OFFSET + additional_offset] = label_id_by_word_index(word_index)
  end

  set_children_data = (config_word_index, offset_number) -> begin
    get_left_child(config.tree, config_word_index) |> word_index -> set_word_data_by_index_with_offset(word_index, offset_number * STACK_NUMBER + 1)
    get_right_child(config.tree, config_word_index) |> word_index -> set_word_data_by_index_with_offset(word_index, offset_number * STACK_NUMBER + 2)
    get_left_child(config.tree, config_word_index, child_number=2) |> word_index -> set_word_data_by_index_with_offset(word_index, offset_number * STACK_NUMBER + 3)
    get_right_child(config.tree, config_word_index, child_number=2) |> word_index -> set_word_data_by_index_with_offset(word_index, offset_number * STACK_NUMBER + 4)
  end

  set_children_data(get_stack_element(config, 1), 0)
  set_children_data(get_buffer_element(config, 1), 1)

  batch
end

function take_batch_embeddings(model, batch)
  take_embedding(word_index::Integer) = view(model.embeddings, word_index, :)
  batch = map(take_embedding, batch)

  batch
end

function get_word_id(model::Model, word::String)
  haskey(model.word_ids, word) ? model.word_ids[word] : model.word_ids[UNKNOWN_TOKEN]
end

function get_tag_id(model::Model, tag::String)
  haskey(model.tag_ids, tag) ? model.tag_ids[tag] : model.tag_ids[UNKNOWN_TOKEN]
end

function get_label_id(model::Model, label::String)
  label == EMPTY_LABEL ? model.label_ids[NULL_TOKEN] : model.label_ids[label]
end

function try_unsafe_free(x)
  hasmethod(CUDA.unsafe_free!, Tuple{typeof(x)}) && CUDA.unsafe_free!(x)
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
abstract type AbstractPacket end
abstract type AbstractRegister end
abstract type AbstractValueHolder end
mutable struct Packet <: AbstractPacket
    val::Float64
end

function __repr__(self::AbstractPacket)::String
    return AbstractPacket(self.val)
end

function __eq__(self::AbstractPacket, other::AbstractPacket)::Bool
    return __key(self) == __key(other)
end

function __lt__(self::AbstractPacket, other::AbstractPacket)::Bool
    return __key(self) < __key(other)
end

function __le__(self::AbstractPacket, other::AbstractPacket)::Bool
    return __key(self) <= __key(other)
end

function __gt__(self::AbstractPacket, other::AbstractPacket)::Bool
    return __key(self) > __key(other)
end

function __ge__(self::AbstractPacket, other::AbstractPacket)::Bool
    return __key(self) >= __key(other)
end

function __key(self::AbstractPacket)
    (self.val)
end

mutable struct Register <: AbstractRegister
    PACKET::AbstractPacket
    VALUE::Int64
end

function __repr__(self::AbstractRegister)::String
    return AbstractRegister(self.PACKET, self.VALUE)
end

function __eq__(self::AbstractRegister, other::AbstractRegister)::Bool
    return __key(self) == __key(other)
end

function __lt__(self::AbstractRegister, other::AbstractRegister)::Bool
    return __key(self) < __key(other)
end

function __le__(self::AbstractRegister, other::AbstractRegister)::Bool
    return __key(self) <= __key(other)
end

function __gt__(self::AbstractRegister, other::AbstractRegister)::Bool
    return __key(self) > __key(other)
end

function __ge__(self::AbstractRegister, other::AbstractRegister)::Bool
    return __key(self) >= __key(other)
end

function __key(self::AbstractRegister)
    (__key(self.PACKET), self.VALUE)
end

mutable struct ValueHolder <: AbstractValueHolder
    val::Int64
    strVal::String
end

function __repr__(self::AbstractValueHolder)::String
    return AbstractValueHolder(self.val, self.strVal)
end

function __eq__(self::AbstractValueHolder, other::AbstractValueHolder)::Bool
    return __key(self) == __key(other)
end

function __lt__(self::AbstractValueHolder, other::AbstractValueHolder)::Bool
    return __key(self) < __key(other)
end

function __le__(self::AbstractValueHolder, other::AbstractValueHolder)::Bool
    return __key(self) <= __key(other)
end

function __gt__(self::AbstractValueHolder, other::AbstractValueHolder)::Bool
    return __key(self) > __key(other)
end

function __ge__(self::AbstractValueHolder, other::AbstractValueHolder)::Bool
    return __key(self) >= __key(other)
end

function __key(self::AbstractValueHolder)
    (self.val, self.strVal)
end

if abspath(PROGRAM_FILE) == @__FILE__
    c1 = ValueHolder(10, "10")
    @assert(__eq__(c1, ValueHolder(10, "10")))
    c2 = ValueHolder(10, "10")
    @assert(__le__(c1, c2))
    @assert(__ge__(c1, c2))
    c1 = ValueHolder(11, "10")
    @assert(__lt__(c2, c1))
    @assert(__gt__(c1, c2))
    c3 = Packet(2.4)
    @assert(__eq__(c3, Packet(2.4)))
    c4 = Packet(2.4)
    @assert(__le__(c3, c4))
    @assert(__ge__(c3, c4))
    c3 = Packet(2.5)
    @assert(__lt__(c4, c3))
    @assert(__gt__(c3, c4))
    c5 = Register(Packet(2.2), 10)
    @assert(__eq__(c5, Register(Packet(2.2), 10)))
    c6 = Register(Packet(2.2), 10)
    @assert(__le__(c5, c6))
    @assert(__ge__(c5, c6))
    c5 = Register(Packet(2.3), 10)
    @assert(__lt__(c6, c5))
    @assert(__gt__(c5, c6))
end


text:
### Utilities used by @match macro
# author: Kevin Squire (@kmsquire)

#
# ismatch
#

ismatch(r::AbstractRange, s) = s in r
ismatch(c::Char, s::Number) = false
ismatch(s::Number, c::Char) = false
ismatch(r::Regex, s::AbstractString) = occursin(r, s)
ismatch(r, s) = (r == s)

#
# slicedim
#
# "view" version of slicedim

function _slicedim(A::AbstractArray, d::Integer, i::Integer)
    if (d < 1) | (d > ndims(A))
        throw(BoundsError())
    end
    sz = size(A)
    # Force 1x..x1 slices to extract the value
    # Note that this is no longer a reference.
    otherdims = [sz...]
    splice!(otherdims, d)
    if all(otherdims .== 1)
        A[[ n == d ? i : 1 for n in 1:ndims(A) ]...]
    else
        view(A, [ n == d ? i : (1:sz[n]) for n in 1:ndims(A) ]...)
    end
end

function _slicedim(A::AbstractArray, d::Integer, i)
    if (d < 1) | (d > ndims(A))
        throw(BoundsError())
    end
    sz = size(A)
    view(A, [ n == d ? i : (1:sz[n]) for n in 1:ndims(A) ]...)
end

_slicedim(A::AbstractVector, d::Integer, i::Integer) =
    (if (d < 0) | (d > 1);  throw(BoundsError()) end;  A[i])

_slicedim(A::AbstractVector, d::Integer, i) =
    (if (d < 0) | (d > 1);  throw(BoundsError()) end;  view(A, i))

function slicedim(A::AbstractArray, s::Integer, i::Integer, from_end::Bool=false)
    d = s + max(ndims(A) - 2, 0)
    from_end && (i = size(A, d) - i)
    _slicedim(A, d, i)
end

function slicedim(A::AbstractArray, s::Integer, i::Integer, j::Integer)
    d = s + max(ndims(A) - 2, 0)
    j = size(A, d) - j # this is the distance from the end of the dim size
    _slicedim(A, d, i:j)
end

#
# getvars
#
# get all symbols in an expression

getvars(e)         = Symbol[]
getvars(e::Symbol) = startswith(string(e), '@') ? Symbol[] : Symbol[e]

function getvars(e::Expr)
    if isexpr(e, :call)
        getvars(e.args[2:end])
    else
        getvars(e.args)
    end
end
getvars(es::AbstractArray) = unique!(vcat([getvars(e) for e in es]...))

#
# getvar
#
# get the symbol from a :(::) expression

getvar(x::Expr) = isexpr(x, :(::)) ? x.args[1] : x
getvar(x::Symbol) = x

#
# check_dim_size_expr
#
# generate an expression to check the size of a variable dimension against an array of expressions

function check_dim_size_expr(val, dim, ex::Expr)
    if length(ex.args) == 0 || !any([isexpr(e, :(...)) for e in ex.args])
        :(Match.checkdims($val, $dim, $(length(ex.args))))
    else
        :(Match.checkdims2($val, $dim, $(length(ex.args))))
    end
end

function checkdims(val::AbstractArray, dim, dimsize)
    dim = dim + max(ndims(val) - 2, 0)
    dim <= ndims(val) && size(val, dim) == dimsize
end

checkdims(val, dim, dimsize) = false

function checkdims2(val::AbstractArray, dim, dimsize)
    dim = dim + max(ndims(val) - 2, 0)
    dim <= ndims(val) && size(val, dim) >= dimsize - 1
end

checkdims2(val, dim, dimsize) = false

#
# check_tuple_len_expr
#
# generate an expression to check the length of a tuple variable against a tuple expression

function check_tuple_len_expr(val, ex::Expr)
    if length(ex.args) == 0 || !any([isexpr(e, :(...)) for e in ex.args])
        :(length($val) == $(length(ex.args)))
    else
        :(length($val) >= $(length(ex.args) - 1))
    end
end

function check_tuple_len(val::Expr, ex::Expr)
    if !isexpr(val, :tuple) || !isexpr(ex, :tuple)
        false
    elseif length(ex.args) == 0 || !any([isexpr(e, :(...)) for e in ex.args])
        length(val.args) == length(ex.args)
    else
        length(val.args) >= length(ex.args) - 1
    end
end

#
# joinexprs
#
# join an array of (e.g., true/false) expressions with an operator

function joinexprs(exprs::AbstractArray, oper::Symbol, default=:nothing)
    len = length(exprs)

    len == 0 ? default :
    len == 1 ? exprs[1] :
               Expr(oper, joinexprs(view(exprs, 1:(len - 1)), oper, default), exprs[end])
end

#
# let_expr
#
# generate an optional let expression

function let_expr(expr, assignments::AbstractArray)
    length(assignments) == 0 && return expr
    assignment_block = Expr(:block, assignments...)
    return Expr(:let, assignment_block, expr)
end

#
# array_type_of
#
# modify x::Type => x::AbstractArray{Type}

function array_type_of(ex::Expr)
    if isexpr(ex, :(::))
        :($(ex.args[1])::AbstractArray{$(ex.args[2])})
    else
        ex
    end
end

array_type_of(sym::Symbol) = :($sym::AbstractArray)


text:
module CodeCompressor

using MLStyle

export compresscode, compresssource, @compresscode_str

@data CodeState begin
    InitCode
    NormalCode(instring::Bool, bracketnest::Int, sind::Int)
    NormalString(instring::Bool, istriple::Bool, sind::Int)
    CommandString(instring::Bool, istriple::Bool, sind::Int)
    Comment(instring::Bool, islines::Bool, sind::Int)
end
#strings = Set([normalstring, triplestring, commandstring])

"""
    compresscode(code; deletespace=true, deletebreak=true, deletecomment=true)

compresscode returns compressed code.

If you want to compress jl file, you use `compresssource`.

# Examples
```jldoctest
julia> compresscode(\"""
          a =   1 +1 #hogemi
       #= memomemo =#
       println("a $("b" * "c")?", a, true ? 10 : 1);
       \""")
" a = 1 +1 ;;println(\"a bc?\", a, true ? 10 : 1);;"

```
"""
function compresscode(code; deletespace=true, deletebreak=true, deletecomment=true)
    begin
        eval(code)
    end

    io = IOBuffer()
    state = CodeState[InitCode]
    codelength = length(code)
    ind = 1
    bracketcount = 0
    wio(s) = write(io, s)
    while ind <= codelength
        @match state[end] begin
            InitCode => begin
                while ind <= codelength
                    c = code[ind]
                    if !deletespace && c in " \t"
                        break
                    elseif !deletebreak && c == '\n'
                        break
                    else
                        break
                    end
                    ind += 1
                end
                pop!(state)
                push!(state, NormalCode(false, 0, ind))
            end
            NormalCode(instring, bracketnest, sind) => begin
                c = code[ind]
                writeup = true
                adds = ""
                bind = ind-1
                # for commen
                if c == '#'
                    if code[ind+1] == '='
                        push!(state, Comment(instring, true, ind))
                    else
                        eind = ind
                        while eind <= codelength && code[eind] != '\n'
                            eind += 1
                        end
                        ind = eind
                        if deletecomment
                            writeup = false
                        else
                            bind = eind-1
                        end
                        pop!(state)
                        push!(state, NormalCode(instring, bracketcount, ind))
                    end
                # for string
                elseif c == '"'
                    if code[ind+1:ind+2] == "\"\""
                        push!(state, NormalString(instring, true, ind))
                    else
                        push!(state, NormalString(instring, false, ind))
                    end
                # for command
                elseif c == '`'
                    if code[ind+1:ind+2] == "``"
                        push!(state, CommandString(instring, true, ind))
                    else
                        push!(state, CommandString(instring, false, ind))
                    end
                elseif c == '('
                    bracketcount += 1
                    writeup = false
                    ind += 1
                elseif c == ')'
                    # for code in string
                    bracketcount -= 1
                    ind += 1
                    writeup = false
                    if instring && bracketnest == bracketcount
                        pop!(state)
                    end
                # for break
                elseif deletebreak && c == '\n'
                    adds = ";"
                    ind += 1
                    pop!(state)
                    push!(state, NormalCode(instring, bracketcount, ind))
                # for space
                elseif deletespace && c in " \t"
                    while c in " \t"
                        ind += 1
                        c = code[ind]
                    end
                    adds = " "
                    pop!(state)
                    push!(state, NormalCode(instring, bracketcount, ind))
                else
                    writeup = false
                    ind += 1
                end
                if writeup && !instring
                    wio(code[sind:bind])
                    if adds != ""
                        wio(adds)
                    end
                end
            end
            NormalString(instring, istriple, sind) => begin
                bind = ind
                while ind <= codelength
                    if !istriple && code[ind] == '"'
                        bind = ind
                        ind += 1
                        pop!(state)
                        push!(state, NormalCode(instring, bracketcount, ind))
                        break
                    elseif istriple && code[ind:ind+2] == "\"\"\""
                        bind = ind+2
                        ind += 3
                        pop!(state)
                        push!(state, NormalCode(instring, bracketcount, ind))
                        break
                    else code[ind:ind+1] == "\$("
                        bind = ind
                        ind += 1
                        push!(NormalCode(true, bracketcount, ind))
                        break
                    end
                    ind += 1
                end
                if !instring
                    wio(code[sind:bind])
                end
            end
            CommandString(instring, istriple, sind) => begin
                c = code[ind]
            end
            Comment(instring, islines, sind) && if islines end => begin
                eind = ind+2
                while eind+1 <= codelength && code[eind:eind+1] != "=#"
                    eind += 1
                end
                if !deletecomment && !instring
                    wio(code[sind:eind+1])
                end
                ind = eind+2
                pop!(state)
                push!(state, NormalCode(instring, bracketcount, ind))
            end
        end
    end

    compcode = String(take!(io))
    close(io)

    return compcode
end

"""
    compresssource(filename, deletespace=true, deletebreak=true, deletecomment=false)

compresssource reads the jl source file and returns the compressed code.

"""
function compresssource(filename, deletespace=true, deletebreak=true, deletecomment=true)
    code = read(filename) |> String
    return compresscode(code; deletespace=deletespace, deletebreak=deletebreak, deletecomment=deletecomment)
end

"""
@compresscode_str is a macro to return compressed code.

"""
macro compresscode_str(code)
    compresscode(code)
end

end # module


text:
module Graphics

using Reexport
using Requires

include("Graphviz.jl")
include("GraphvizWiringDiagrams.jl")
include("TikZ.jl")
include("TikZWiringDiagrams.jl")

@reexport using .GraphvizWiringDiagrams
@reexport using .TikZWiringDiagrams

function __init__()
  @require TikzPictures="37f6aa50-8035-52d0-81c2-5a1d08754b2d" include("TikZExternal.jl")
end

end


text:
#!/usr/bin/env julia

using Unitful, UnitfulAstro

table = (("Erg", :erg, :J),
         ("Dyne", :dyn, :N),
         ("Arcminute", :arcminute, :°),
         ("Arcsecond", :arcsecond, :°),
         ("Julian year", :yr, :d),
         ("Astronomical unit", :AU, :cm),
         ("Light year", :ly, :cm),
         ("Parsec", :pc, :cm),
         ("Jansky", :Jy, u"erg/s/cm^2/Hz"),
         ("Angstrom", :angstrom, :nm),
         ("Solar radius", :Rsun, :cm),
         ("Solar irradiance", :Ssun, u"erg/s/cm^2"),
         ("Solar luminosity", :Lsun, u"erg/s"),
         ("G × Solar mass", :GMsun, u"m^3/s^2"),
         ("Solar mass", :Msun, :g),
         ("Earth radius (equatorial)", :Rearth_e, :cm),
         ("Earth radius (polar)", :Rearth_p, :cm),
         ("Earth radius", :Rearth, :cm),
         ("G × Earth mass", :GMearth, u"m^3/s^2"),
         ("Earth mass", :Mearth, :g),
         ("Jupiter radius (equatorial)", :Rjup_e, :cm),
         ("Jupiter radius (polar)", :Rjup_p, :cm),
         ("Jupiter radius", :Rjup, :cm),
         ("G × Jupiter mass", :GMjup, u"m^3/s^2"),
         ("Jupiter mass", :Mjup, :g),
         ("Fifty-one ergs", :foe, :erg),
         ("Solar flux unit", :SFU, :kJy),
         ("Total electron content unit", :TECU, u"cm^-2"))

column1 = String[]
column2 = String[]
column3 = String[]

for (name, unit, base) in table
    if isa(base, Symbol)
        value = uconvert(getfield(UnitfulAstro, base), 1*getfield(UnitfulAstro, unit))
    else
        value = uconvert(base, 1*getfield(UnitfulAstro, unit))
    end
    push!(column1, name)
    push!(column2, "`UnitfulAstro.$(string(unit))`")
    push!(column3, "`$(repr(value))`")
end

A = maximum(length.(column1))
B = maximum(length.(column2))
C = maximum(length.(column3))

println("| $(lpad("Name", A)) | $(lpad("Binding", B)) | $(lpad("Value", C)) |")
println("|", "-"^(A+2), "|", "-"^(B+2), "|", "-"^(C+2), "|")
for (a, b, c) in zip(column1, column2, column3)
    println("| $(lpad(a, A)) | $(lpad(b, B)) | $(lpad(c, C)) |")
end


text:
struct DirectMCEstimate
    num_samples::Int
end

name(x::DirectMCEstimate) = "Direct MC"

struct DirectMCResult{Samples} <: SimulationResult
    samples::Samples
end

log_marginal(result::DirectMCResult{<:AbstractMatrix}) = vec(logmeanexp(result.samples, dims=2))
log_marginal(result::DirectMCResult{<:AbstractVector}) = logmeanexp(result.samples)

function simulate(algorithm::DirectMCEstimate, initial, system; kwargs...)
    samples = zeros(Float64, algorithm.num_samples)
    for i in 1:algorithm.num_samples
        signal = sample(initial, system)
        samples[i] = -energy_difference(signal, system)
    end
    DirectMCResult(samples)
end


text:
export FQE

mutable struct FQE{
    P<:GaussianNetwork,
    C<:NeuralNetworkApproximator,
    C_T<:NeuralNetworkApproximator,
    R<:AbstractRNG,
 } <: AbstractLearner
    policy::P
    q_network::C
    target_q_network::C_T
    n_evals::Int
    γ::Float32
    batch_size::Int
    update_freq::Int
    update_step::Int
    tar_update_freq::Int
    rng::R
    #logging
    loss::Float32
end

"""
    FQE(;kwargs...)

See [Hyperparameter Selection for Offline Reinforcement Learning](https://arxiv.org/abs/2007.09055)

# Keyword arguments
- `policy`, the policy for which FQE should be performed.
- `q_network`, critic to evaluate the Q value of `state`, `action` pair.
- `target_q_network`, target critic used for evaluating target Q values.
- `n_evals::Int`, number of evaluations to perform to return the performance of the policy.
- `γ::Float32 = 0.99f0`, discount factor.
- `batch_size::Int = 32`.
- `update_freq::Int = 50`, frequency of updating the `target_q_network`.
- `update_step::Int = 0`.
- `tar_update_freq::Int = 50`
- `rng::AbstractRNG = Range.GLOBAL_RNG`.

`policy` is expected to be a pre-trained [`GaussianNetwork`](@ref) with a particular choice of hyperparameters
preferrably trained using the same `dataset`.
"""
function FQE(;
    policy,
    q_network,
    target_q_network,
    n_evals=20,
    γ=0.99f0,
    batch_size=32,
    update_freq=1,
    update_step=0,
    tar_update_freq = 50,
    rng=Random.GLOBAL_RNG,
)
    copyto!(q_network, target_q_network) #force sync
    FQE(
        policy,
        q_network,
        target_q_network,
        n_evals,
        γ,
        batch_size,
        update_freq,
        update_step,
        tar_update_freq,
        rng,
        0.0f0,
    )
end

Flux.functor(x::FQE) = (Q = x.q_network, Qₜ = x.target_q_network),
y -> begin
    x = @set x.q_network = y.Q
    x = @set x.target_q_network = y.Qₜ
    x
end

function (l::FQE)(env)
    s = send_to_device(device(l.policy), state(env))
    s = Flux.unsqueeze(s, ndims(s) + 1)
    action = dropdims(l.policy(l.rng, s; is_sampling=true), dims=2)
end

function (l::FQE)(env, ::Val{:Eval})
    results = []
    D = device(l.policy)
    for _ in 1:l.n_evals
        reset!(env)
        s = send_to_device(D, state(env))
        s = Flux.unsqueeze(s, ndims(s)+1)
        a = dropdims(l.policy(l.rng, s; is_sampling=true), dims=2)
        input = vcat(s, a)
        result = l.q_network(input)
        push!(results, result[])
    end
    mean(results), results
end

function (l::FQE)(state::AbstractArray, action::AbstractArray)
    D = device(l.q_network)
    s = send_to_device(D, state)
    a = send_to_device(D, reshape(action, :, l.batch_size))
    input = vcat(s, a)
    value = l.q_network(input)
end

function RLBase.update!(l::FQE, batch::NamedTuple{SARTS})
    policy = l.policy
    Q, Qₜ = l.q_network, l.target_q_network

    D = device(Q)
    s, a, r, t, s′ = (send_to_device(D, batch[x]) for x in SARTS)
    γ = l.γ
    batch_size = l.batch_size

    loss_func = Flux.Losses.mse

    q′ = Qₜ(vcat(s′, policy(s′)[1])) |> vec

    target = r .+ γ .* (1 .- t) .* q′

    gs = gradient(params(Q)) do
        q = Q(vcat(s, reshape(a, :, batch_size))) |> vec
        loss = loss_func(q, target)
        Zygote.ignore() do
            l.loss = loss
        end
        loss
    end
    Flux.Optimise.update!(Q.optimizer, params(Q), gs)

    if l.update_step % l.tar_update_freq == 0
        Qₜ = deepcopy(Q)
    end
end


text:
struct Literal
    name::Int
    pos::Bool
end
Base.:!(l::Literal) = Literal(l.name,!l.pos)
Base.show(io::IO, l::Literal) = print(io, l.pos ? '+' : '-', l.name)

struct Clause
    lits::Array{Literal,1}
end
Base.length(c::Clause) = length(c.lits)
Base.iterate(c::Clause) = (c.lits[begin], firstindex(c.lits))
Base.iterate(c::Clause, i) = i >= lastindex(c.lits) ? nothing : (c.lits[i+1], i+1)
Base.getindex(c::Clause, i...) = c.lits[i...]
Base.show(io::IO, c::Clause) = begin
    print(io, '(')

    for (i,l) in enumerate(c)
        print(io, l)
        i != length(c) && print(io, " ∨ ")
    end

    print(io, ')')
end

Base.in(l::Literal, c::Clause) = l ∈ c.lits

struct Formula
    clauses::Array{Clause,1}
end
Base.length(f::Formula) = length(f.clauses)
Base.iterate(f::Formula) = (f.clauses[begin], firstindex(f.clauses))
Base.iterate(f::Formula, i) = i >= lastindex(f.clauses) ? nothing : (f.clauses[i+1], i+1)
Base.getindex(f::Formula, i...) = f.clauses[i...]
Base.show(io::IO, f::Formula) = begin
    print(io, '(')

    for (i,l) in enumerate(f)
        print(io, l)
        i != length(f) && print(io, " ∧\n ")
    end

    print(io, ')')
end

struct Trail
    decs::Vector{Tuple{Literal, Bool}}
end
Base.isempty(t::Trail) = isempty(t.decs)
Base.iterate(t::Trail) = iterate(t.decs)
Base.iterate(t::Trail, i) = iterate(t.decs, i)
Base.length(t::Trail) = length(t.decs)
Base.getindex(t::Trail, i...) = t.decs[i...]
Base.push!(t::Trail, x) = push!(t.decs, x)
Base.filter(f, t::Trail) = filter(f, t.decs)

Base.in(l::Literal, t::Trail) = any(l == x for (x,_) in t)

Base.show(io::IO, t::Trail) = begin
    print(io, '[')

    for (i,(l,d)) in enumerate(t)
        d && print(io, '|')
        print(io, l)
        i != length(t) && print(io, ", ")
    end

    print(io ,']')
end

# convenience constructors
Literal(literal::Int) = Literal(abs(literal), literal > 0)
Clause(v::Vector{Int}) = Clause(Literal.(v))
Formula(v::Vector{Vector{Int}}) = Formula(Clause.(v))
solve(v::Vector{Vector{Int}}) = solve(Formula(v))


text:
"""
CircStats.Datasets includes a number of example circular statistical datasets.
"""
module Datasets

"""Azimuths at which 714 mallards disappeared from sight.

#### Reference
- Matthews, G.V.T. (1961).  "Nonsense" orientation in mallard anas platyrhynchos
  and its relation to experiments on bird navigation.  Ibis, 103a, 211-230.
"""
const Matthews_1961_mallards = (0:20:340, [40, 22, 20, 9, 6, 3, 3, 1, 6, 3, 11, 22, 24, 58, 136, 138, 143, 69])

end # module


text: ispangram(input::AbstractString) = length(Set(matchall(r"[a-z]", lowercase(input)))) == 26