text:
using TimeSeries, MarketData
FactCheck.setstyle(:compact)
FactCheck.onlystats(true)

facts("time series methods") do

    context("lag takes previous day and timestamps it to next day") do
        @fact lag(cl).values[1]    --> roughly(111.94, atol=.01)
        @fact lag(cl).timestamp[1] --> Date(2000,1,4)
    end

    context("lag accepts other offset values") do
        @fact lag(cl, 9).timestamp[1] --> Date(2000,1,14)
    end

    context("lag operates on 2d arrays") do
        @fact lag(ohlc, 9).timestamp[1] --> Date(2000,1,14)
    end

    context("lag returns 1d from 1d time arrays") do
        @fact ndims(lag(cl).values) --> 1
    end

    context("lag returns 2d from 2d time arrays") do
        @fact ndims(lag(ohlc).values) --> 2
    end

    context("lead takes next day and timestamps it to current day") do
        @fact lead(cl).values[1]    --> roughly(102.5, atol=.1)
        @fact lead(cl).timestamp[1] --> Date(2000,1,3)
    end

    context("lead accepts other offset values") do
        @fact lead(cl, 9).values[1]    --> 100.44
        @fact lead(cl, 9).timestamp[1] --> Date(2000,1,3)
    end

    context("lead operates on 2d arrays") do
        @fact lead(ohlc, 9).timestamp[1] --> Date(2000,1,3)
    end

    context("lead returns 1d from 1d time arrays") do
        @fact ndims(lead(cl).values) --> 1
    end

    context("lead returns 2d from 2d time arrays") do
        @fact ndims(lead(ohlc).values) --> 2
    end

    context("diff calculates 1st-order differences") do
        @fact diff(op).timestamp                              --> diff(op, padding=false).timestamp
        @fact diff(op).values                                 --> diff(op, padding=false).values
        @fact diff(op, padding=false).values[1]               --> op[2].values[1] .- op[1].values[1]
        @fact isequal(diff(op, padding=true).values[1], NaN)  --> true
        @fact diff(op, padding=true).values[2]                --> diff(op).values[1]
        @fact diff(op, padding=true).values[2]                --> op[2].values[1] .- op[1].values[1]
    end

    context("correct simple return value") do
        @fact percentchange(cl, :simple).values              --> percentchange(cl).values
        @fact percentchange(cl).values                       --> percentchange(cl, padding=false).values
        @fact percentchange(cl).values[1]                    --> roughly((102.5-111.94)/111.94, atol=.01)
        @fact percentchange(ohlc).values[1, :]               --> roughly((ohlc.values[2,:] - ohlc.values[1,:]) ./ ohlc.values[1,:])
        @fact percentchange(cl, padding=true).values[1]      --> isnan
        @fact percentchange(cl, padding=true).values[2]      --> roughly((102.5-111.94)/111.94, atol=.01)
        @fact percentchange(ohlc, padding=true).values[2, :] --> roughly((ohlc.values[2,:] - ohlc.values[1,:]) ./ ohlc.values[1,:])
    end

    context("correct log return value") do
        @fact percentchange(cl, :log).values                       --> percentchange(cl, :log, padding=false).values
        @fact percentchange(cl, :log).values[1]                    --> roughly(log(102.5) - log(111.94), atol=.01)
        @fact percentchange(ohlc, :log).values[1, :]               --> roughly(log(ohlc.values[2,:]) - log(ohlc.values[1,:]), atol=.01)
        @fact percentchange(cl, :log, padding=true).values[1]      --> isnan
        @fact percentchange(cl, :log, padding=true).values[2]      --> roughly(log(102.5) - log(111.94))
        @fact percentchange(ohlc, :log, padding=true).values[2, :] --> roughly(log(ohlc.values[2,:]) - log(ohlc.values[1,:]))
    end

    context("moving supplies correct window length") do
        @fact moving(cl, mean, 10).values                                                      --> moving(cl, mean, 10, padding=false).values
        @fact moving(cl, mean, 10).timestamp[1]                                                --> Date(2000,1,14)
        @fact moving(cl, mean, 10).values[1]                                                   --> roughly(mean(cl.values[1:10]))
        @fact moving(cl, mean, 10, padding=true).timestamp[1]                                  --> Date(2000,1, 3)
        @fact moving(cl, mean, 10, padding=true).timestamp[10]                                 --> Date(2000,1,14)
        @fact isequal(moving(cl, mean, 10, padding=true).values[1], NaN)                       --> true
        @fact moving(cl, mean, 10, padding=true).values[10]                                    --> moving(cl, mean, 10).values[1]
        @fact moving(ohlc, mean, 10).values                                                    --> moving(ohlc, mean, 10, padding=false).values
        @fact moving(ohlc, mean, 10).values[1, :]'                                             --> roughly(mean(ohlc.values[1:10, :], 1))
        @fact isequal(moving(ohlc, mean, 10, padding=true).values[1, :], [NaN, NaN, NaN, NaN]) --> true
        @fact moving(ohlc, mean, 10, padding=true).values[10, :]                               --> moving(ohlc, mean, 10).values[1, :]
    end

    context("upto method accumulates") do
        @fact upto(cl, sum).values[10]        --> roughly(sum(cl.values[1:10]))
        @fact upto(cl, mean).values[10]       --> roughly(mean(cl.values[1:10]))
        @fact upto(cl, sum).timestamp[10]     --> Date(2000,1,14)
        # transpose the upto value output from column to row vector but values are identical
        @fact upto(ohlc, sum).values[10, :]'  --> roughly(sum(ohlc.values[1:10, :], 1))
        @fact upto(ohlc, mean).values[10, :]' --> roughly(mean(ohlc.values[1:10, :], 1))
    end
end

facts("base element-wise operators on TimeArray values") do

    context("only values on intersecting Dates computed") do
        @fact (cl[1:2] ./ op[2:3]).values[1] --> roughly(0.94688222)
        @fact (cl[1:4] .+ op[4:7]).values[1] --> roughly(201.12, atol=.01)
        @fact length(cl[1:2] ./ op[2:3])     --> 1
        @fact length(cl[1:4] .+ op[4:7])     --> 1
    end

    context("correct unary operation on TimeArray values") do
        @fact (+cl).values[1]                --> cl.values[1]
        @fact (-cl).values[1]                --> -cl.values[1]
        @fact (!(cl .== op)).values[1]       --> true
        @fact log(cl).values[1]              --> log(cl.values[1])
        @fact sqrt(cl).values[1]             --> sqrt(cl.values[1])
        @fact (+ohlc).values[1,:]            --> ohlc.values[1,:]
        @fact (-ohlc).values[1,:]            --> -(ohlc.values[1,:])
        @fact (!(ohlc .== ohlc)).values[1,1] --> false
        @fact log(ohlc).values[1, :]         --> log(ohlc.values[1, :])
        @fact sqrt(ohlc).values[1, :]        --> sqrt(ohlc.values[1, :])
        @fact_throws sqrt(-ohlc)
    end

    context("correct dot operation between TimeArray values and Int/Float and viceversa") do
        @fact (cl .- 100).values[1]     --> roughly(11.94, atol=.01)
        @fact (cl .+ 100).values[1]     --> roughly(211.94, atol=.01)
        @fact (cl .* 100).values[1]     --> roughly(11194, atol=1)
        @fact (cl ./ 100).values[1]     --> roughly(1.1194, atol=.0001)
        @fact (cl .^ 2).values[1]       --> roughly(12530.5636, atol=.0001)
        @fact (cl .% 2).values[1]       --> cl.values[1] % 2
        @fact (100 .- cl).values[1]     --> roughly(-11.94, atol=.01)
        @fact (100 .+ cl).values[1]     --> roughly(211.94, atol=.01)
        @fact (100 .* cl).values[1]     --> roughly(11194, atol=.01)
        @fact (100 ./ cl).values[1]     --> roughly(0.8933357155619082)
        @fact (2 .^ cl).values[1]       --> 4980784073277740581384811358191616
        @fact (2 .% cl).values[1]       --> 2
        @fact (ohlc .- 100).values[1,:] --> ohlc.values[1,:] .- 100
        @fact (ohlc .+ 100).values[1,:] --> ohlc.values[1,:] .+ 100
        @fact (ohlc .* 100).values[1,:] --> ohlc.values[1,:] .* 100
        @fact (ohlc ./ 100).values[1,:] --> ohlc.values[1,:] ./ 100
        @fact (ohlc .^ 2).values[1,:]   --> ohlc.values[1,:] .^ 2
        @fact (ohlc .% 2).values[1,:]   --> ohlc.values[1,:] .% 2
        @fact (100 .- ohlc).values[1,:] --> 100 .- ohlc.values[1,:]
        @fact (100 .+ ohlc).values[1,:] --> 100 .+ ohlc.values[1,:]
        @fact (100 .* ohlc).values[1,:] --> 100 .* ohlc.values[1,:]
        @fact (100 ./ ohlc).values[1,:] --> 100 ./ ohlc.values[1,:]
        @fact (2 .^ ohlc).values[1,:]   --> 2 .^ ohlc.values[1,:]
        @fact (2 .% ohlc).values[1,:]   --> 2 .% ohlc.values[1,:]
    end

    context("correct non-dot operation between TimeArray values and Int/Float and viceversa") do
        @fact (cl - 100).values[1] --> roughly(11.94, atol=0.1)
        @fact (cl + 100).values[1] --> roughly(211.94, atol=0.1)
        @fact (cl * 100).values[1] --> roughly(11194, atol=0.1)
        @fact (cl / 100).values[1] --> roughly(1.1194, atol=0.001)
        @fact (cl % 2).values[1]   -->  cl.values[1] % 2
        # not supported by Base - reserved for square matrix multiplication
        @fact_throws (cl ^ 2).values[1]
        @fact (100 - cl).values[1] --> roughly(-11.94, atol=0.1)
        @fact (100 + cl).values[1] --> roughly(211.94, atol=0.1)
        @fact (100 * cl).values[1] --> roughly(11194, atol=0.1)
        # not supported by Base - confusion with matrix inverse
        @fact_throws (100 / cl).values
        # not supported by Base
        @fact_throws (2 ^ cl).values
        @fact (2 % cl).values[1]       --> 2 % cl.values[1]
        @fact (ohlc - 100).values[1,:] --> ohlc.values[1,:] - 100
        @fact (ohlc + 100).values[1,:] --> ohlc.values[1,:] + 100
        @fact (ohlc * 100).values[1,:] --> ohlc.values[1,:] * 100
        @fact (ohlc / 100).values[1,:] --> ohlc.values[1,:] / 100
        # not supported by Base - reserved for square matrix multiplication
        @fact_throws (ohlc ^ 2).values[1,:]
        @fact (ohlc % 2).values[1,:]   --> ohlc.values[1,:] % 2
        @fact (100 - ohlc).values[1,:] --> 100 - ohlc.values[1,:]
        @fact (100 + ohlc).values[1,:] --> 100 + ohlc.values[1,:]
        @fact (100 * ohlc).values[1,:] --> 100 * ohlc.values[1,:]
        # not supported by Base - confusion with matrix inverse
        @fact_throws (100 / ohlc).values[1,:]
        # not supported by Base
        @fact_throws (2 ^ ohlc).values[1,:]
        @fact (2 % ohlc).values[1,:] --> 2 % ohlc.values[1,:]
    end

    context("correct mathematical operations between two same-column-count TimeArrays") do
        @fact (cl .+ op).values[1]                              --> roughly(216.82, atol=.01)
        @fact (cl .- op).values[1]                              --> roughly(7.06, atol=.01)
        @fact (cl .* op).values[1]                              --> roughly(11740.2672, atol=0.0001)
        @fact (cl ./ op).values[1]                              --> roughly(1.067315, atol=0.001)
        @fact (cl .% op).values                                 --> cl.values .% op.values
        @fact (cl .^ op).values                                 --> cl.values .^ op.values
        @fact (cl .* (cl.> 200)).values                         --> cl.values .* (cl.values .> 200)
        @fact (basecall(cl, x->round(Int, x)) .* cl).values     --> round(Int, cl.values) .* cl.values
        @fact (ohlc .+ ohlc).values                             --> ohlc.values .+ ohlc.values
        @fact (ohlc .- ohlc).values                             --> ohlc.values .- ohlc.values
        @fact (ohlc .* ohlc).values                             --> ohlc.values .* ohlc.values
        @fact (ohlc ./ ohlc).values                             --> ohlc.values ./ ohlc.values
        @fact (ohlc .% ohlc).values                             --> ohlc.values .% ohlc.values
        @fact (ohlc .^ ohlc).values                             --> ohlc.values .^ ohlc.values
        @fact (ohlc .* (ohlc .> 200)).values                    --> ohlc.values .* (ohlc.values .> 200)
        @fact (basecall(ohlc, x->round(Int, x)) .* ohlc).values --> round(Int, ohlc.values) .* ohlc.values
    end

    context("correct broadcasted mathematical operations between different-column-count TimeArrays") do
        @fact (ohlc .+ cl).values                             --> ohlc.values .+ cl.values
        @fact (ohlc .- cl).values                             --> ohlc.values .- cl.values
        @fact (ohlc .* cl).values                             --> ohlc.values .* cl.values
        @fact (ohlc ./ cl).values                             --> ohlc.values ./ cl.values
        @fact (ohlc .% cl).values                             --> ohlc.values .% cl.values
        @fact (ohlc .^ cl).values                             --> ohlc.values .^ cl.values
        @fact (ohlc .* (cl.> 200)).values                     --> ohlc.values .* (cl.values .> 200)
        @fact (basecall(ohlc, x->round(Int, x)) .* cl).values --> round(Int, ohlc.values) .* cl.values
        @fact (cl .+ ohlc).values                             --> cl.values .+ ohlc.values
        @fact (cl .- ohlc).values                             --> cl.values .- ohlc.values
        @fact (cl .* ohlc).values                             --> cl.values .* ohlc.values
        @fact (cl ./ ohlc).values                             --> cl.values ./ ohlc.values
        @fact (cl .% ohlc).values                             --> cl.values .% ohlc.values
        @fact (cl .^ ohlc).values                             --> cl.values .^ ohlc.values
        @fact (cl .* (ohlc .> 200)).values                    --> cl.values .* (ohlc.values .> 200)
        @fact (basecall(cl, x->round(Int, x)) .* ohlc).values --> round(Int, cl.values) .* ohlc.values
        @fact_throws (ohlc["Open", "Close"] .+ ohlc) # One array must have a single column
    end

    context("correct comparison operations between TimeArray values and Int/Float (and viceversa)") do
        @fact (cl .> 111.94).values[1]      --> false
        @fact (cl .< 111.94).values[1]      --> false
        @fact (cl .>= 111.94).values[1]     --> true
        @fact (cl .<= 111.94).values[1]     --> true
        @fact (cl .== 111.94).values[1]     --> true
        @fact (cl .!= 111.94).values[1]     --> false
        @fact (111.94 .> cl).values[1]      --> false
        @fact (111.94 .< cl).values[1]      --> false
        @fact (111.94 .>= cl).values[1]     --> true
        @fact (111.94 .<= cl).values[1]     --> true
        @fact (111.94 .== cl).values[1]     --> true
        @fact (111.94 .!= cl).values[1]     --> false
        @fact (ohlc .> 111.94).values[1,:]  --> [false, true, false, false]
        @fact (ohlc .< 111.94).values[1,:]  --> [true, false, true, false]
        @fact (ohlc .>= 111.94).values[1,:] --> [false, true, false, true]
        @fact (ohlc .<= 111.94).values[1,:] --> [true, false, true, true]
        @fact (ohlc .== 111.94).values[1,:] --> [false, false, false, true]
        @fact (ohlc .!= 111.94).values[1,:] --> [true, true, true, false]
        @fact (111.94 .> ohlc).values[1,:]  --> [true, false, true, false]
        @fact (111.94 .< ohlc).values[1,:]  --> [false, true, false, false]
        @fact (111.94 .>= ohlc).values[1,:] --> [true, false, true, true]
        @fact (111.94 .<= ohlc).values[1,:] --> [false, true, false, true]
        @fact (111.94 .== ohlc).values[1,:] --> [false, false, false, true]
        @fact (111.94 .!= ohlc).values[1,:] --> [true, true, true, false]
    end

    context("correct comparison operations between TimeArray values and Bool (and viceversa)") do
        @fact ((cl .> 111.94) .== true).values[1]     --> false
        @fact ((cl .> 111.94) .!= true).values[1]     --> true
        @fact (true .== (cl .> 111.94)).values[1]     --> false
        @fact (true .!= (cl .> 111.94)).values[1]     --> true
        @fact ((ohlc .> 111.94).== true).values[1,:]  --> [false, true, false, false]
        @fact ((ohlc .> 111.94).!= true).values[1,:]  --> [true, false, true, true]
        @fact (true .== (ohlc .> 111.94)).values[1,:] --> [false, true, false, false]
        @fact (true .!= (ohlc .> 111.94)).values[1,:] --> [true, false, true, true]
    end

    context("correct comparison operations between same-column-count TimeArrays") do
        @fact (cl .> op).values[1]        --> true
        @fact (cl .< op).values[1]        --> false
        @fact (cl .<= op).values[1]       --> false
        @fact (cl .>= op).values[1]       --> true
        @fact (cl .== op).values[1]       --> false
        @fact (cl .!= op).values[1]       --> true
        @fact (ohlc .> ohlc).values[1,:]  --> [false, false, false, false]
        @fact (ohlc .< ohlc).values[1,:]  --> [false, false, false, false]
        @fact (ohlc .<= ohlc).values[1,:] --> [true, true, true, true]
        @fact (ohlc .>= ohlc).values[1,:] --> [true, true, true, true]
        @fact (ohlc .== ohlc).values[1,:] --> [true, true, true, true]
        @fact (ohlc .!= ohlc).values[1,:] --> [false, false, false, false]
    end

    context("correct comparison operations between different-column-count TimeArrays") do
        @fact (ohlc .> cl).values  --> ohlc.values .> cl.values
        @fact (ohlc .< cl).values  --> ohlc.values .< cl.values
        @fact (ohlc .>= cl).values --> ohlc.values .>= cl.values
        @fact (ohlc .<= cl).values --> ohlc.values .<= cl.values
        @fact (ohlc .== cl).values --> ohlc.values .== cl.values
        @fact (ohlc .!= cl).values --> ohlc.values .!= cl.values
        @fact (cl .> ohlc).values  --> cl.values .> ohlc.values
        @fact (cl .< ohlc).values  --> cl.values .< ohlc.values
        @fact (cl .>= ohlc).values --> cl.values .>= ohlc.values
        @fact (cl .<= ohlc).values --> cl.values .<= ohlc.values
        @fact (cl .== ohlc).values --> cl.values .== ohlc.values
        @fact (cl .!= ohlc).values --> cl.values .!= ohlc.values
        @fact_throws (ohlc["Open", "Close"] .== ohlc) # One array must have a single column
    end

    context("correct bitwise elementwise operations between bool and TimeArrays' values") do
        @fact ((cl .> 100) & true).values[1]      --> true
        @fact ((cl .> 100) | true).values[1]      --> true
        @fact ((cl .> 100) $ true).values[1]      --> false
        @fact (false & (cl .> 100)).values[1]     --> false
        @fact (false | (cl .> 100)).values[1]     --> true
        @fact (false $ (cl .> 100)).values[1]     --> true
        @fact ((ohlc .> 100) & true).values[4,:]  --> [true, true, false, false]
        @fact ((ohlc .> 100) | true).values[4,:]  --> [true, true, true, true]
        @fact ((ohlc .> 100) $ true).values[4,:]  --> [false, false, true, true]
        @fact (false & (ohlc .> 100)).values[4,:] --> [false, false, false, false]
        @fact (false | (ohlc .> 100)).values[4,:] --> [true, true, false, false]
        @fact (false $ (ohlc .> 100)).values[4,:] --> [true, true, false, false]
      end

    context("correct bitwise elementwise operations between same-column-count TimeArrays' boolean values") do
        @fact ((cl .> 100) & (cl .< 120)).values[1]       --> true
        @fact ((cl .> 100) | (cl .< 120)).values[1]       --> true
        @fact ((cl .> 100) $ (cl .< 120)).values[1]       --> false
        @fact ((ohlc .> 100) & (ohlc .< 120)).values[4,:] --> [true, true, false, false]
        @fact ((ohlc .> 100) | (ohlc .< 120)).values[4,:] --> [true, true, true, true]
        @fact ((ohlc .> 100) $ (ohlc .< 120)).values[4,:] --> [false, false, true, true]
        @fact_throws ((ohlc .> 100) $ (cl.< 120)) # Bitwise broadcasting not supported by Base
    end

end

facts("basecall works with Base methods") do

    context("cumsum works") do
        @fact basecall(cl, cumsum).values[2] --> cl.values[1] + cl.values[2]
    end

    context("log works") do
        @fact basecall(cl, log).values[2] --> log(cl.values[2])
    end
end

facts("adding/removing missing rows works") do

    uohlc = uniformspace(ohlc[1:8])

    context("uniform spacing detection works") do
        @fact datetime1 --> uniformspaced
        @fact uohlc     --> uniformspaced
        @fact cl        --> not(uniformspaced)
        @fact ohlc      --> not(uniformspaced)
    end

    context("forcing uniform spacing works") do
        @fact length(uohlc)                               --> 10
        @fact uohlc[5].values                             --> ohlc[5].values
        @fact isequal(uohlc[6].values, [NaN NaN NaN NaN]) --> true
        @fact isequal(uohlc[7].values, [NaN NaN NaN NaN]) --> true
        @fact uohlc[8].values                             --> ohlc[6].values
    end

    context("dropnan works") do
        nohlc = TimeArray(ohlc.timestamp, copy(ohlc.values), ohlc.colnames, ohlc.meta)
        nohlc.values[7:12, 2] = NaN

        @fact dropnan(uohlc).timestamp       --> dropnan(uohlc, :all).timestamp
        @fact dropnan(uohlc).values          --> dropnan(uohlc, :all).values

        @fact dropnan(ohlc, :all).values     --> ohlc.values
        @fact dropnan(nohlc, :all).timestamp --> ohlc.timestamp
        @fact dropnan(uohlc, :all).timestamp --> ohlc[1:8].timestamp
        @fact dropnan(uohlc, :all).values    --> ohlc[1:8].values

        @fact dropnan(ohlc, :any).values     --> ohlc.values
        @fact dropnan(nohlc, :any).timestamp --> ohlc.timestamp[[1:6;13:end]]
        @fact dropnan(nohlc, :any).values    --> ohlc.values[[1:6;13:end], :]
        @fact dropnan(uohlc, :any).timestamp --> ohlc[1:8].timestamp
        @fact dropnan(uohlc, :any).values    --> ohlc[1:8].values
    end
end


text:
@testset "Low Level Utils" begin
    # clear bit
    for i in 0:63
        @test CounterTools.clearbit(typemax(UInt64), i) == typemax(UInt64) - UInt64(2)^i
    end

    # Set bit
    for i in 0:63
        @test CounterTools.setbit(UInt64(0), i) == UInt64(2)^(i)
    end

    # clear and set bits
    @test CounterTools.clearbits(typemax(UInt64), 0:63) == zero(UInt64)
    @test CounterTools.setbits(zero(UInt64), 0:63) == typemax(UInt64)

    @test CounterTools.mask(2,4) == (1 << 2) + (1 << 3) + (1 << 4)
    m = CounterTools.mask(2, 4)
    @test CounterTools.isbitset(m, 2) == true
    @test CounterTools.isbitset(m, 3) == true
    @test CounterTools.isbitset(m, 4) == true
    @test CounterTools.isbitset(m, 1) == false
    @test CounterTools.isbitset(m, 5) == false
end

# Test writing and reading MSRs
@testset "MSRs" begin
    # Get the current state of counters on CPU 2
    #
    # Avoid CPU 1 (IndexZero(0)) because we're messing with the counter settings.
    # CPU 1 often has a watchdog timer ...
    cpu = CounterTools.IndexZero(1)
    ctrl = CounterTools.IA32_PERF_GLOBAL_CTRL_MSR

    initial_state = CounterTools.readmsr(cpu, ctrl)

    CounterTools.disablecounters(cpu)
    @test CounterTools.readmsr(cpu, ctrl) == 0

    # Enable the counters
    CounterTools.enablecounters(cpu)
    config = 0x0000000700000000 | (CounterTools.mask(CounterTools.numcounters()) - 1)
    @test CounterTools.readmsr(cpu, ctrl) == config

    CounterTools.writemsr(cpu, ctrl, initial_state)
end

@testset "Counter Values" begin
    # a = CounterTools.CounterValue(10)
    # b = CounterTools.CounterValue(20)
    # @test b - a == 10
    # @test a - b == 2^47 - 10
end

# @testset "Testing MCFG" begin
#     @test sizeof(CounterTools.MCFGRecord) == 16
#     @test sizeof(CounterTools.MCFGHeader) == 44
# end

ext: jl


text:
using LinearAlgebra
using TensorOperations

using Test
using puMPS
using puMPS.MPS

include("mps.jl")

max_stars_repo_name: FuTen/puMPS

max_issues_repo_name: FuTen/puMPS

max_forks_repo_name: FuTen/puMPS


text:
export BrickletIO4V2Configuration
struct BrickletIO4V2Configuration
    direction::Char
    value::Bool
end

export BrickletIO4V2InputValueCallbackConfiguration
struct BrickletIO4V2InputValueCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletIO4V2AllInputValueCallbackConfiguration
struct BrickletIO4V2AllInputValueCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletIO4V2Monoflop
struct BrickletIO4V2Monoflop
    value::Bool
    time::Integer
    time_remaining::Integer
end

export BrickletIO4V2EdgeCountConfiguration
struct BrickletIO4V2EdgeCountConfiguration
    edge_type::Integer
    debounce::Integer
end

export BrickletIO4V2PWMConfiguration
struct BrickletIO4V2PWMConfiguration
    frequency::Integer
    duty_cycle::Integer
end

export BrickletIO4V2SPITFPErrorCount
struct BrickletIO4V2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletIO4V2Identity
struct BrickletIO4V2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletIO4V2
"""
4-channel digital input/output
"""
mutable struct BrickletIO4V2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletIO4V2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_io4_v2")
        deviceInternal = package.BrickletIO4V2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_value
"""
    $(SIGNATURES)

Sets the output value of all four channels. A value of *true* or *false* outputs
logic 1 or logic 0 respectively on the corresponding channel.

Use :func:`Set Selected Value` to change only one output channel state.

For example: (True, True, False, False) will turn the channels 0-1 high and the
channels 2-3 low.

All running monoflop timers and PWMs will be aborted if this function is called.

.. note::
 This function does nothing for channels that are configured as input. Pull-up
 resistors can be switched on with :func:`Set Configuration`.
"""
function set_value(device::BrickletIO4V2, value)
    device.deviceInternal.set_value(value)
end

export get_value
"""
    $(SIGNATURES)

Returns the logic levels that are currently measured on the channels.
This function works if the channel is configured as input as well as if it is
configured as output.
"""
function get_value(device::BrickletIO4V2)
    return device.deviceInternal.get_value()
end

export set_selected_value
"""
    $(SIGNATURES)

Sets the output value of a specific channel without affecting the other channels.

A running monoflop timer or PWM for the specific channel will be aborted if this
function is called.

.. note::
 This function does nothing for channels that are configured as input. Pull-up
 resistors can be switched on with :func:`Set Configuration`.
"""
function set_selected_value(device::BrickletIO4V2, channel, value)
    device.deviceInternal.set_selected_value(channel, value)
end

export set_configuration
"""
    $(SIGNATURES)

Configures the value and direction of a specific channel. Possible directions
are 'i' and 'o' for input and output.

If the direction is configured as output, the value is either high or low
(set as *true* or *false*).

If the direction is configured as input, the value is either pull-up or
default (set as *true* or *false*).

For example:

* (0, 'i', true) will set channel 0 as input pull-up.
* (1, 'i', false) will set channel 1 as input default (floating if nothing is connected).
* (2, 'o', true) will set channel 2 as output high.
* (3, 'o', false) will set channel 3 as output low.

A running monoflop timer or PWM for the specific channel will be aborted if this
function is called.
"""
function set_configuration(device::BrickletIO4V2, channel, direction, value)
    device.deviceInternal.set_configuration(channel, direction, value)
end

export get_configuration
"""
    $(SIGNATURES)

Returns the channel configuration as set by :func:`Set Configuration`.
"""
function get_configuration(device::BrickletIO4V2, channel)
    return device.deviceInternal.get_configuration(channel)
end

export set_input_value_callback_configuration
"""
    $(SIGNATURES)

This callback can be configured per channel.

The period is the period with which the :cb:`Input Value`
callback is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change within the
period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_input_value_callback_configuration(device::BrickletIO4V2, channel, period, value_has_to_change)
    device.deviceInternal.set_input_value_callback_configuration(channel, period, value_has_to_change)
end

export get_input_value_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration for the given channel as set by
:func:`Set Input Value Callback Configuration`.
"""
function get_input_value_callback_configuration(device::BrickletIO4V2, channel)
    return device.deviceInternal.get_input_value_callback_configuration(channel)
end

export set_all_input_value_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`All Input Value`
callback is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change within the
period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_all_input_value_callback_configuration(device::BrickletIO4V2, period, value_has_to_change)
    device.deviceInternal.set_all_input_value_callback_configuration(period, value_has_to_change)
end

export get_all_input_value_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by
:func:`Set All Input Value Callback Configuration`.
"""
function get_all_input_value_callback_configuration(device::BrickletIO4V2)
    return device.deviceInternal.get_all_input_value_callback_configuration()
end

export set_monoflop
"""
    $(SIGNATURES)

The first parameter is the desired state of the channel (*true* means output *high*
and *false* means output *low*). The second parameter indicates the time that
the channel should hold the state.

If this function is called with the parameters (true, 1500):
The channel will turn on and in 1.5s it will turn off again.

A PWM for the selected channel will be aborted if this function is called.

A monoflop can be used as a failsafe mechanism. For example: Lets assume you
have a RS485 bus and a IO-4 Bricklet 2.0 is connected to one of the slave
stacks. You can now call this function every second, with a time parameter
of two seconds. The channel will be *high* all the time. If now the RS485
connection is lost, the channel will turn *low* in at most two seconds.
"""
function set_monoflop(device::BrickletIO4V2, channel, value, time)
    device.deviceInternal.set_monoflop(channel, value, time)
end

export get_monoflop
"""
    $(SIGNATURES)

Returns (for the given channel) the current value and the time as set by
:func:`Set Monoflop` as well as the remaining time until the value flips.

If the timer is not running currently, the remaining time will be returned
as 0.
"""
function get_monoflop(device::BrickletIO4V2, channel)
    return device.deviceInternal.get_monoflop(channel)
end

export get_edge_count
"""
    $(SIGNATURES)

Returns the current value of the edge counter for the selected channel. You can
configure the edges that are counted with :func:`Set Edge Count Configuration`.

If you set the reset counter to *true*, the count is set back to 0
directly after it is read.

.. note::
 Calling this function is only allowed for channels configured as input.
"""
function get_edge_count(device::BrickletIO4V2, channel, reset_counter)
    return device.deviceInternal.get_edge_count(channel, reset_counter)
end

export set_edge_count_configuration
"""
    $(SIGNATURES)

Configures the edge counter for a specific channel.

The edge type parameter configures if rising edges, falling edges or
both are counted if the channel is configured for input. Possible edge types are:

* 0 = rising
* 1 = falling
* 2 = both

Configuring an edge counter resets its value to 0.

If you don't know what any of this means, just leave it at default. The
default configuration is very likely OK for you.

.. note::
 Calling this function is only allowed for channels configured as input.
"""
function set_edge_count_configuration(device::BrickletIO4V2, channel, edge_type, debounce)
    device.deviceInternal.set_edge_count_configuration(channel, edge_type, debounce)
end

export get_edge_count_configuration
"""
    $(SIGNATURES)

Returns the edge type and debounce time for the selected channel as set by
:func:`Set Edge Count Configuration`.

.. note::
 Calling this function is only allowed for channels configured as input.
"""
function get_edge_count_configuration(device::BrickletIO4V2, channel)
    return device.deviceInternal.get_edge_count_configuration(channel)
end

export set_pwm_configuration
"""
    $(SIGNATURES)

Activates a PWM for the given channel.

You need to set the channel to output before you call this function, otherwise it will
report an invalid parameter error. To turn the PWM off again, you can set the frequency to 0 or any other
function that changes a value of the channel (e.g. :func:`Set Selected Value`).

A running monoflop timer for the given channel will be aborted if this function
is called.
"""
function set_pwm_configuration(device::BrickletIO4V2, channel, frequency, duty_cycle)
    device.deviceInternal.set_pwm_configuration(channel, frequency, duty_cycle)
end

export get_pwm_configuration
"""
    $(SIGNATURES)

Returns the PWM configuration as set by :func:`Set PWM Configuration`.
"""
function get_pwm_configuration(device::BrickletIO4V2, channel)
    return device.deviceInternal.get_pwm_configuration(channel)
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletIO4V2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletIO4V2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletIO4V2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletIO4V2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletIO4V2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletIO4V2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletIO4V2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletIO4V2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletIO4V2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletIO4V2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletIO4V2)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletIO4V2)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletIO4V2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
"""
This is a terrible function. But is there any other way to reliable get the
abstract supertype of an arbitrary type hierarchy without loosing performance?
"""
@generated function go_abstract{T<:AbstractGeometry}(::Type{T})
    ff = T
    while ff.name.name != :AbstractGeometry
       ff = supertype(ff)
    end
    :($ff)
end

@generated function eltype_or{T<:GeometryPrimitive}(::Type{T}, OR)
    ET = go_abstract(T).parameters[2]
    if isa(ET, TypeVar)
        return :(OR)
    else
        return :($ET)
    end
end
@generated function ndims_or{T<:GeometryPrimitive}(::Type{T}, OR)
    N = go_abstract(T).parameters[1]
    if isa(N, TypeVar)
        return :(OR)
    else
        return :($N)
    end
end

function Base.eltype{T<:AbstractGeometry}(x::Type{T})
    eltype_or(T, Any)
end
@generated function Base.ndims{T<:AbstractGeometry}(x::Type{T})
    N = go_abstract(T).parameters[1]
    isa(N, TypeVar) && error("Ndims not given for type $T")
    :($N)
end

Base.eltype{N, T}(x::AbstractGeometry{N, T}) = T
Base.ndims{N, T}(x::AbstractGeometry{N, T}) = N

Base.eltype{T}(::Type{AFG{T}}) = T
Base.eltype{FG <: AFG}(::Type{FG}) = eltype(supertype(FG))


text:
# ExpSmoothing.jl - Tipo para computar suavizamiento exponencial de medidas de inflación

"""
    InflationExpSmoothing{F <: InflationFunction} <: InflationFunction
Función de inflación que computa la inflación con decaimiento exponencial con el
parámetro ```\\alpha`` (el parámetro del decaimiento ) de la trayectoria
interanual dada por la medida de inflación `inflfn` almacenada.

## Ejemplo

Para computar la inflación con decaimiento exponencial de la variación
interanual del IPC:
```julia-repl
julia> inflfn = InflationExpSmoothing(InflationTotalCPI(), 0.8)
(::InflationExpSmoothing{InflationTotalCPI}) (generic function with 5 methods)
```
"""
struct InflationExpSmoothing{F <: InflationFunction} <: InflationFunction
    inflfn::F
    alpha::Float64
end

# 2. Extender el método de nombre
measure_name(esfn::InflationExpSmoothing) = "Suavizamiento exponencial con parámetro $(round(esfn.alpha, digits=4)) parámetros de " * measure_name(esfn.inflfn)

# Función de parámetros
CPIDataBase.params(esfn::InflationExpSmoothing) = (esfn.alpha, )

# Método de conveniencia para definir sobre parámetro de suavizamiento Enteros
InflationExpSmoothing(inflfn::InflationFunction,alpha::Int) =
    InflationExpSmoothing(inflfn::InflationFunction, convert(Float64,alpha))

# Método que opera sobre CountryStructure: computa la trayectoria de inflación
# con la función inflfn y luego computa el promedio con suavizamiento exponencial con el parámetro
# alpha (el parámetro del decaimiento )
function (esfn::InflationExpSmoothing)(cs::CountryStructure)

    # Cómputo usual de inflación
    tray_infl = esfn.inflfn(cs)

    # Algoritmo de promedio con suavizamiento exponencial
    k = esfn.alpha
    k == 1 && return tray_infl
    k == 0 && return ones(length(tray_infl))*tray_infl[1]

    # Computar el promedio con suavizamiento exponencial (in-place)
    es_tray_infl = smoothing_exponential(tray_infl, k)

    # Devolver la trayectoria de inflación
    es_tray_infl
end

# Actualmente, se define el resumen intermensual por el obtenido con la función
# de inflación interna
# function (mafn::InflationExpSmoothing)(base::VarCPIBase)
#     # Obtener resumen intermensual y aplicarle el suavizamiento exponencial
#     vm = mafn.inflfn(base)
#     vm
# end

# Función para computar el promedio con suavizamiento exponencial de la serie `v` con `k`
# como valor de suavizamiento que pondera la historia 0<=k<=1.
# Esta función modifica los elementos de `v` con el promedio de suavizmiento exponencial
# con alpha dada.
# El tipo de suavizamiento que se utiliza es SES (Simple Exponential Smoothing, o
# Suavizamiento Exponencial Simple en español)

function smoothing_exponential(v, k)
    es = similar(v)
    n = length(es)
    es[1] = v[1]
    for j = 2:(n)
        es[j] = es[j-1] + (v[j] - es[(j-1)])*k
    end
    es
end


text:
isinteractive() && include("common.jl")

#####
##### test diagnostics
#####

@testset "summarize tree statistics" begin
    N = 1000
    directions = Directions(UInt32(0))
    function rand_invalidtree()
        if rand() < 0.1
            REACHED_MAX_DEPTH
        else
            left = rand(-5:5)
            right = left + rand(0:5)
            InvalidTree(left, right)
        end
    end
    tree_statistics = [TreeStatisticsNUTS(randn(), rand(0:5), rand_invalidtree(), rand(),
                                          rand(1:30), directions) for _ in 1:N]
    stats = summarize_tree_statistics(tree_statistics)
    # acceptance rates
    @test stats.N == N
    @test stats.a_mean ≈ mean(x -> x.acceptance_rate, tree_statistics)
    @test stats.a_quantiles ==
        quantile((x -> x.acceptance_rate).(tree_statistics), ACCEPTANCE_QUANTILES)
    # termination counts
    @test stats.termination_counts.divergence ==
        count(x -> is_divergent(x.termination), tree_statistics)
    @test stats.termination_counts.max_depth ==
        count(x -> x.termination == REACHED_MAX_DEPTH, tree_statistics)
    @test stats.termination_counts.turning ==
        (N - stats.termination_counts.max_depth - stats.termination_counts.divergence)
    # depth counts
    for (i, c) in enumerate(stats.depth_counts)
        @test count(x -> x.depth == i - 1, tree_statistics) == c
    end
    @test sum(stats.depth_counts) == N
    # misc
    @test 1.8 ≤ EBFMI(tree_statistics) ≤ 2.2 # nonsensical value, just checking calculation
    @test repr(stats) isa AbstractString # just test that it prints w/o error
end

@testset "log acceptance ratios" begin
    ℓ = multivariate_normal(ones(5))
    log2ϵs = -5:5
    N = 13
    logA = explore_log_acceptance_ratios(ℓ, zeros(5), log2ϵs; N = N)
    @test all(isfinite.(logA))
    @test size(logA) == (length(log2ϵs), N)
end

@testset "leapfrog trajectory" begin
    # problem setup
    K = 2
    ℓ = multivariate_normal(ones(K))
    κ = GaussianKineticEnergy(K)
    q = zeros(K)
    Q = evaluate_ℓ(ℓ, q)
    p = ones(K) .* 0.98
    H = Hamiltonian(κ, ℓ)
    ϵ = 0.1
    ixs = 1:15
    ix0 = 5

    # calculate trajectory manually
    zs1 = let z = PhasePoint(Q, p)
        [(z1 = z; z = leapfrog(H, z, ϵ); z1) for _ in ixs]
    end
    πs1 = logdensity.(Ref(H), zs1)
    Δs1 = πs1 .- πs1[ix0]

    # calculate using function
    traj = leapfrog_trajectory(ℓ, zs1[ix0].Q.q, ϵ, ixs .- ix0; κ = κ, p = zs1[ix0].p)
    @test all(isapprox.(map(t -> t.Δ, traj), Δs1; atol = 1e-5))
    @test all(map((t, y) -> t.z.Q.q ≈ y.Q.q, traj, zs1))
    @test all(map((t, y) -> t.z.p ≈ y.p, traj, zs1))
end


text:
###
# LCX Oracle
###
# At the moment only supports cutting planes
import MathProgBase.AbstractMathProgSolver

export LCXOracle

#returns bool, violated ab_cut
#Solves 4 LPs
function ab_cut_(us, vs, z, data, eps_, solver, CASE)
    const TOL = 1e-6
    m = solver==nothing ? Model() : Model(solver=solver)
    @defVar(m, a[1:size(data, 2)])
    @defVar(m, b)

    #define the RHS
    const N = size(data, 1); const d = size(data, 2)
    @defVar(m, t[1:N] >= 0)
    for i = 1:N
        @addConstraint(m, t[i] >= z*(dot(a, vec(data[i, :])) - b)/N)
    end

    #a,b in a box
    @defVar(m, abs_a[1:d])
    for i = 1:d
        @addConstraint(m, abs_a[i] >=  a[i])
        @addConstraint(m, abs_a[i] >= -a[i])
    end
    #Altered from original paper to match the new sampling technique
    @addConstraint(m, sum{abs_a[i], i=1:d} <= 1.)

    data_max = maximum(abs(data))

    if CASE == :1
        @addConstraint(m, dot(a, vs) - (z-1) * b >= 0)
        @addConstraint(m, dot(a, us) - b >= 0)

        #implied bounds on b to help maintain numerical stability
        @addConstraint(m, b >= -data_max)

        @setObjective(m, :Max, dot(a, vs) - (z-1)*b + dot(a, us) - b - sum{t[i], i=1:N})
    elseif CASE == :2
        @addConstraint(m, dot(a, vs) - (z-1) * b >= 0)
        @addConstraint(m, dot(a, us) - b <= 0)

        #implied bounds on b to help maintain numerical stability
        @addConstraint(m, (z-1) * b <= norm(vs, Inf))
        if abs(z-1) <= TOL
            @addConstraint(m, b <= data_max)
        end

        @setObjective(m, :Max, dot(a, vs) - (z-1)*b - sum{t[i], i=1:N})
    elseif CASE == :3
        @addConstraint(m, dot(a, vs) - (z-1) * b <= 0)
        @addConstraint(m, dot(a, us) - b >= 0)

        #implied bounds on b to help maintain numerical stability
        @addConstraint(m, (z-1) * b >= -norm(vs, Inf) )
        if abs(z-1) <= TOL
            @addConstraint(m, b >= -data_max)
        end

        @setObjective(m, :Max, dot(a, us) - b - sum{t[i], i=1:N})
    end
    status = solve(m)
    status != :Optimal && error("Case $CASE optimization failed with status $status")
    getObjectiveValue(m), getValue(a), getValue(b)
end

#returns the first violation
function ab_cut(us, vs, z, data, eps_, Gamma, TOL, solver = nothing)
    obj, astar, bstar = ab_cut_(us, vs, z, data, eps_, solver, :1)
    if obj > Gamma + TOL
        return obj, astar[:], bstar
    end

    obj, astar, bstar = ab_cut_(us, vs, z, data, eps_, solver, :2)
    if obj > Gamma + TOL
        return obj, astar[:], bstar
    end

    obj, astar, bstar = ab_cut_(us, vs, z, data, eps_, solver, :3)
    if obj > Gamma + TOL
        return obj, astar[:], bstar
    end
    return Gamma, Float64[], 0.0
end

#returns zstar,ustar
#Rebuilds the model each time bc no way to remove cnsts from JuMP
function suppFcnLCX(xs, data, eps_, Gamma, cut_sense;
                    lbounds=Float64[], ubounds=Float64[],
                    TOL=1e-6, MAXITER=1000, trace=false, outer_solver=nothing, abcut_solver=nothing)
    #assign bounds if the user didn't give em
    const d = size(data, 2)
    data_max = maximum(abs(data))

    m = outer_solver==nothing ? Model() : Model(solver=outer_solver)
    @defVar(m, us[1:d])
    @defVar(m, vs[1:d])
    @defVar(m, 1 <= z <= 1/eps_)
    @setObjective(m, cut_sense, dot(xs, us))

    #any user bounds
    !isempty(lbounds) && @addConstraint(m, lbounds[i] <= us[i=1:d])
    !isempty(ubounds) && @addConstraint(m, ubounds[i] >= us[i=1:d])

    #Implied bounds on u,v, added for stability
    for ix = 1:d
        @addConstraint(m,  vs[ix] <= data_max * z + Gamma)
        @addConstraint(m, -vs[ix] <= data_max * z + Gamma)
        @addConstraint(m,  us[ix] <= data_max * z + Gamma)
        @addConstraint(m, -us[ix] <= data_max * z + Gamma)
    end

    status = solve(m)
    status != :Optimal && error("LCX support init solve failed: $status")

    uvals = getValue(us)
    vvals = getValue(vs)
    zval  = getValue(z)

    obj, astar, bstar = ab_cut(uvals[:], vvals[:], zval, data, eps_, Gamma, TOL, abcut_solver)

    iter::Int64 = 0
    while obj > Gamma + TOL
        iter > MAXITER && error("LCX Supp: Max Iterations exceeded, $(obj - Gamma)")
        iter += 1

        @defVar(m, t[1:2] >= 0)

        @addConstraint(m, t[1] >= dot(astar[:], vs) - bstar * z + bstar)
        @addConstraint(m, t[2] >= dot(astar[:], us) - bstar)
        @addConstraint(m, t[1] + t[2] <= mean(max(data * astar[:] - bstar, 0))*z + Gamma)

        status = solve(m)
        status != :Optimal && error("LCX support Inner Solver Failed: $status")
        uvals = getValue(us)
        vvals = getValue(vs)
        zval  = getValue(z)

        obj, astar, bstar = ab_cut(uvals[:], vvals[:], zval, data, eps_, Gamma, TOL, abcut_solver)

        trace && println("Iter: $iter \t $(obj-Gamma) \t $(astar') \t $(bstar)") # \t $(uvals[:]) \t $(vvals[:]) \t $zval
    end
    getObjectiveValue(m), getValue(us[:])
end

########################
type LCXOracle <: AbstractOracle
    data::Array{Float64, 2}
    eps_::Float64
    Gamma::Float64
    cut_tol::Float64         # 1e-6
    max_iter::Int64          # 100
    ab_cut_tol::Float64      # 1e-6
    abcut_solver::Union{AbstractMathProgSolver, Void}
    outer_solver::Union{AbstractMathProgSolver, Void}

    lbounds::Vector{Float64}
    ubounds::Vector{Float64}

    # Other options
    debug_printcut::Bool
    trace::Bool

end

#Preferred Interface
function LCXOracle(data, eps_, alpha;
                    cut_tol=1e-6, max_iter=100, ab_cut_tol=1e-6,
                    lbounds=Float64[], ubounds=Float64[],
                    debug_printcut=false, trace=false,
                    numSamples=10_000, numBoots=10_000,
                    Gamma =-1.0, abcut_solver=nothing, outer_solver=nothing)
    @assert (0 < eps_ < 1) "Epsilon invalid: $eps_"
    Gamma_ = Gamma < 0 ? calc_ab_thresh(data, alpha, numBoots, numSamples) : Gamma
    LCXOracle(data, eps_, Gamma_, cut_tol, max_iter,
                ab_cut_tol, abcut_solver, outer_solver, lbounds, ubounds, debug_printcut, trace)
end

suppFcn(xs, w::LCXOracle, cut_sense) =
    suppFcnLCX(xs, w.data, w.eps_, w.Gamma, cut_sense,
        lbounds=w.lbounds, ubounds=w.ubounds, TOL=w.ab_cut_tol, MAXITER=w.max_iter,
        trace=w.trace, outer_solver=w.outer_solver, abcut_solver=w.abcut_solver)

function setup(w::LCXOracle, rm::Model, prefs)
    rd = JuMPeR.getRobust(rm)
    @assert (rd.numUncs == size(w.data, 2)) "Num Uncertainties $(rd.numUncs) doesn't match columns in data $(size(w.data, 2))"
    @assert (length(rd.uncertaintyset) == 0) "Auxiliary constraints on uncertainties not yet supported"

    #set the bounds on the object now.
    #Shallow copy since they should be "const"
    #only really want them if they're not uniformly Infinite
    if ! all(rd.uncLower .== -Inf)
        w.lbounds = rd.uncLower
    end
    if ! all(rd.uncUpper .== Inf)
        w.ubounds = rd.uncUpper
    end

    # Extract preferences we care about
    w.debug_printcut = get(prefs, :debug_printcut, false)
    w.cut_tol        = get(prefs, :cut_tol, w.cut_tol)
end


text:
"""
MOI constraints
"""

"""
Linear constraints
"""
MOI.supports_constraint(::Optimizer, ::Type{SAF}, ::Type{MOI.EqualTo{Float64}}) = true
# currently only a <= b is supported
MOI.supports_constraint(::Optimizer, ::Type{SAF}, ::Type{MOI.LessThan{Float64}}) = true

function check_inbounds(model::Optimizer, aff::SAF)
	for term in aff.terms
	    check_inbounds(model, term.variable_index)
	end
	return
end

function MOI.add_constraint(model::Optimizer, func::SAF, set::MOI.EqualTo{Float64})
    check_inbounds(model, func)

    if length(func.terms) == 1
        fix!(model.inner, model.variable_info[func.terms[1].variable_index.value], convert(Int64, set.value/func.terms[1].coefficient))
        return MOI.ConstraintIndex{SAF, MOI.EqualTo{Float64}}(0)
    end

    lc = LinearConstraint()
    indices = [v.variable_index.value for v in func.terms]
    coeffs = [v.coefficient for v in func.terms]
    lc.fct = eq_sum
    lc.indices = indices
    lc.coeffs = coeffs
    lc.operator = :(==)
    lc.rhs = set.value
    lc.maxs = zeros(Int, length(indices))
    lc.mins = zeros(Int, length(indices))
    lc.pre_maxs = zeros(Int, length(indices))
    lc.pre_mins = zeros(Int, length(indices))
    # this can be changed later in `set_in_all_different!` but needs to be initialized with false
    lc.in_all_different = false
    lc.idx = length(model.inner.constraints)+1

    push!(model.inner.constraints, lc)

    for (i,ind) in enumerate(lc.indices)
        push!(model.inner.subscription[ind], lc.idx)
    end

    return MOI.ConstraintIndex{SAF, MOI.EqualTo{Float64}}(length(model.inner.constraints))
end

# support for a <= b which is written as a-b <= 0
function MOI.add_constraint(model::Optimizer, func::SAF, set::MOI.LessThan{Float64})
    check_inbounds(model, func)

    if set.upper != 0.0
        error("Only constraints of the type `a <= b` are supported but not `a <= b-2`")
    end
    if length(func.terms) != 2
        error("Only constraints of the type `a <= b` are supported but not `a+b <= c` or something with more terms")
    end
    reverse_order = false
    if func.terms[1].coefficient != 1.0 || func.terms[2].coefficient != -1.0
        if func.terms[1].coefficient == -1.0 && func.terms[2].coefficient == 1.0
            # rhs is lhs and other way around
            reverse_order = true
        else
            error("Only constraints of the type `a <= b` are supported but not `2a <= b`. You used coefficients: $(func.terms[1].coefficient) and $(func.terms[2].coefficient) instead of `1.0` and `-1.0`")
        end
    end

    com = model.inner

    svc = SingleVariableConstraint()
    svc.fct = less_than
    if reverse_order
        svc.lhs = func.terms[2].variable_index.value
        svc.rhs = func.terms[1].variable_index.value
    else
        svc.lhs = func.terms[1].variable_index.value
        svc.rhs = func.terms[2].variable_index.value
    end
    svc.indices = [svc.lhs, svc.rhs]
    svc.idx = length(model.inner.constraints)+1

    push!(model.inner.constraints, svc)

    push!(model.inner.subscription[svc.lhs], svc.idx)
    push!(model.inner.subscription[svc.rhs], svc.idx)

    return MOI.ConstraintIndex{SAF, MOI.LessThan{Float64}}(length(model.inner.constraints))
end

MOI.supports_constraint(::Optimizer, ::Type{MOI.VectorOfVariables}, ::Type{AllDifferentSet}) = true

function MOI.add_constraint(model::Optimizer, vars::MOI.VectorOfVariables, set::AllDifferentSet)
    com = model.inner

    constraint = BasicConstraint()
    constraint.fct = all_different
    constraint.indices = Int[vi.value for vi in vars.variables]
    constraint.idx = length(com.constraints)+1

    push!(com.constraints, constraint)
    for (i,ind) in enumerate(constraint.indices)
        push!(com.subscription[ind], constraint.idx)
    end

    return MOI.ConstraintIndex{MOI.VectorOfVariables, AllDifferentSet}(length(com.constraints))
end

MOI.supports_constraint(::Optimizer, ::Type{SAF}, ::Type{NotEqualSet{Float64}}) = true

function MOI.add_constraint(model::Optimizer, aff::SAF, set::NotEqualSet{Float64})
    if set.value != 0.0
        error("Only constraints of the type `a != b` are supported but not `a != b-2`")
    end
    if length(aff.terms) != 2
        error("Only constraints of the type `a != b` are supported but not `[a,b,...] in NotEqualSet` => only two variables. Otherwise use an AllDifferentSet constraint.")
    end
    if aff.terms[1].coefficient != 1.0 || aff.terms[2].coefficient != -1.0
        error("Only constraints of the type `a != b` are supported but not `2a != b`. You used coefficients: $(aff.terms[1].coefficient) and $(aff.terms[2].coefficient) instead of `1.0` and `-1.0`")
    end

    com = model.inner

    constraint = BasicConstraint()
    constraint.fct = not_equal
    constraint.indices = Int[vi.variable_index.value for vi in aff.terms]
    constraint.idx = length(com.constraints)+1

    push!(com.constraints, constraint)
    for (i,ind) in enumerate(constraint.indices)
        push!(com.subscription[ind], constraint.idx)
    end

    return MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, NotEqualSet{Float64}}(length(com.constraints))
end

function set_pvals!(model::CS.Optimizer)
    com = model.inner
    for constraint in com.constraints
        set_pvals!(com, constraint)
    end
end

function set_constraint_hashes!(model::CS.Optimizer)
    com = model.inner
    for ci=1:length(com.constraints)
        com.constraints[ci].hash = constraint_hash(com.constraints[ci])
    end
end

### !=

sense_to_set(::Function, ::Val{:!=}) = NotEqualSet(0.0)

MOIU.shift_constant(set::NotEqualSet, value) = NotEqualSet(set.value + value)


text:
"""
RuleNode

Type for representing nodes in an expression tree.
"""
mutable struct RuleNode
	ind::Int # index in grammar
	_val::Any  #value of _() evals
	children::Vector{RuleNode}
end

RuleNode(ind::Int) = RuleNode(ind, nothing, RuleNode[])
RuleNode(ind::Int, children::Vector{RuleNode}) = RuleNode(ind, nothing, children)
RuleNode(ind::Int, _val::Any) = RuleNode(ind, _val, RuleNode[])

include("recycler.jl")

function Base.:(==)(A::RuleNode, B::RuleNode)
	(A.ind == B.ind) &&
	    (A._val == B._val) &&
	    (length(A.children) == length(B.children)) && #required because zip doesn't check lengths
	    all(isequal(a,b) for (a,b) in zip(A.children, B.children))
    end

function Base.hash(node::RuleNode, h::UInt=zero(UInt))
	retval = hash(node.ind, h)
	for child in node.children
		retval = hash(child, retval)
	end

	return retval
end

function Base.show(io::IO, node::RuleNode; separator=",", last_child::Bool=false)
	print(io, node.ind)
	if !isempty(node.children)
	    print(io, "{")
	    for (i,c) in enumerate(node.children)
		show(io, c, separator=separator, last_child=(i == length(node.children)))
	    end
	    print(io, "}")
	elseif !last_child
	    print(io, separator)
	end
end

"""
Return the number of vertices in the tree rooted at root.
"""
function Base.length(root::RuleNode)
	retval = 1
	for c in root.children
	    retval += length(c)
	end
	return retval
end

"""
	depth(root::RuleNode)

Return the depth of the expression tree rooted at root.
"""
function depth(root::RuleNode)
	retval = 1
	for c in root.children
	    retval = max(retval, depth(c)+1)
	end
	return retval
end

"""
	node_depth(root::RuleNode, node::RuleNode)

Return the depth of node for an expression tree rooted at root.  Depth is 1 when root == node.
"""
function node_depth(root::RuleNode, node::RuleNode)
	root === node && return 1
	for c in root.children
	    d = node_depth(c, node)
	    d > 0 && (return d+1)
	end
	return 0
end


text:
"""
API:
- `Scheduler`: scheduler object
- `Event`: event object
- `run!`: run simulation
- `remove_event!`: remove event
- `stop!`: stop simulation
- `register!`: add event to scheduler
- `reset!`: reset scheduler
- `replay_events`: print out completed events
- `when types:` now, at, after, every

Use help for more information, e.i.
````julia
] register!
````
"""
module DiscreteEventsLite
    using DataStructures, Printf
    export Scheduler, Event, run!, remove_event!, stop!
    export register!, after, at, now, every, remove_events!, reset!
    export replay_events
    include("structs.jl")
    include("functions.jl")
end


text:
"""
    @ctx <function definition> [label]

Define a function that is context-aware i.e. save the current context
before executing the function and restore the context right before
returning to the caller. So, if the function modifies the context
using (see [`@memo`](@ref)), then the change is not visible the caller.
```
"""
macro ctx(ex, label = nothing)
    def = splitdef(ex)
    name = QuoteNode(label !== nothing ? Symbol(label) : def[:name])
    ctx = gensym()
    def[:body] = quote
        $ctx = ContextTracking.context()
        try
            ContextTracking.save($ctx)
            push!($ctx.path, $name)
            $(def[:body])
        finally
            pop!($ctx.path)
            ContextTracking.restore($ctx)
        end
    end
    return esc(combinedef(def))
end

"""
    @memo var = expr
    @memo var

Store the variable/value from the assigment statement in the current
context.
"""
macro memo(ex)
    # capture the variable
    if ex isa Symbol          # @memo var
        x = ex
    elseif ex.head === :(=)   # @memo var = expression
        x = ex.args[1]
    else
        error("@memo must be followed by an assignment or a variable name.")
    end
    sym = QuoteNode(x)
    val = gensym()
    return esc(quote
        $val = $ex
        push!(ContextTracking.context(), $sym => $val)
    end)
end


text:
module LearningAlgorithms
  files = ["KNeighbors","LinearRegression","KMeans","MultinomialBayes"]
  for f in files
    include("$f.jl")
  end
end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
#Basic Julia in Jupyter

using RDataSets
using Dataframes
using Gadfly
set_default_plot_size(5inch, 5inch/golden);
plot(dataset("datasets","iris"), x="SepalWidth", y="SepalLength", color="Species")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Usrsctp"
version = v"0.9.5"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/sctplab/usrsctp.git", "07f871bda23943c43c9e74cc54f25130459de830")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/usrsctp/
mkdir build && cd build
meson --cross-file="${MESON_TARGET_TOOLCHAIN}" --buildtype=release ..
ninja -j${nproc}
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libusrsctp", :libusrsctp)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# # Machine Learning in Julia (continued)

# An introduction to the
# [MLJ](https://alan-turing-institute.github.io/MLJ.jl/stable/)
# toolbox.

# ### Set-up

# Inspect Julia version:

VERSION

# The following instantiates a package environment.

# The package environment has been created using **Julia 1.6** and may not
# instantiate properly for other Julia versions.

using Pkg
Pkg.activate("env")
Pkg.instantiate()

# ## General resources

# - [MLJ Cheatsheet](https://alan-turing-institute.github.io/MLJ.jl/dev/mlj_cheatsheet/)
# - [Common MLJ Workflows](https://alan-turing-institute.github.io/MLJ.jl/dev/common_mlj_workflows/)
# - [MLJ manual](https://alan-turing-institute.github.io/MLJ.jl/dev/)
# - [Data Science Tutorials in Julia](https://juliaai.github.io/DataScienceTutorials.jl/)

# ## Part 3 - Transformers and Pipelines

# ### Transformers

# Unsupervised models, which receive no target `y` during training,
# always have a `transform` operation. They sometimes also support an
# `inverse_transform` operation, with obvious meaning, and sometimes
# support a `predict` operation (see the clustering example discussed
# [here](https://alan-turing-institute.github.io/MLJ.jl/dev/transformers/#Transformers-that-also-predict-1)).
# Otherwise, they are handled much like supervised models.

# Here's a simple standardization example:

using MLJ

x = rand(100);
@show mean(x) std(x);

#-

model = Standardizer() # a built-in model
mach = machine(model, x)
fit!(mach)
xhat = transform(mach, x);
@show mean(xhat) std(xhat);

# This particular model has an `inverse_transform`:

inverse_transform(mach, xhat) ≈ x

# ### Re-encoding the King County House data as continuous

# For further illustrations of transformers, let's re-encode *all* of the
# King County House input features (see [Ex
# 3](#exercise-3-fixing-scitypes-in-a-table)) into a set of `Continuous`
# features. We do this with the `ContinuousEncoder` model, which, by
# default, will:

# - one-hot encode all `Multiclass` features
# - coerce all `OrderedFactor` features to `Continuous` ones
# - coerce all `Count` features to `Continuous` ones (there aren't any)
# - drop any remaining non-Continuous features (none of these either)

# First, we reload the data and fix the scitypes (Exercise 3):

using UrlDownload, CSV
import DataFrames
house_csv = urldownload("https://raw.githubusercontent.com/ablaom/"*
                        "MachineLearningInJulia2020/for-MLJ-version-0.16/"*
                        "data/house.csv");
house = DataFrames.DataFrame(house_csv)
coerce!(house, autotype(house_csv));
coerce!(house, Count => Continuous, :zipcode => Multiclass);
schema(house)

#-

y, X = unpack(house, ==(:price), rng=123);

# Instantiate the unsupervised model (transformer):

encoder = ContinuousEncoder() # a built-in model; no need to @load it

# Bind the model to the data and fit!

mach = machine(encoder, X) |> fit!;

# Transform and inspect the result:

Xcont = transform(mach, X);
schema(Xcont)

# ### More transformers

# Here's how to list all of MLJ's unsupervised models:

models(m->!m.is_supervised)

# Some commonly used ones are built-in (do not require `@load`ing):

# model type                  | does what?
# ----------------------------|----------------------------------------------
# ContinuousEncoder | transform input table to a table of `Continuous` features (see above)
# FeatureSelector | retain or dump selected features
# FillImputer | impute missing values
# OneHotEncoder | one-hot encoder `Multiclass` (and optionally `OrderedFactor`) features
# Standardizer | standardize (whiten) a vector or all `Continuous` features of a table
# UnivariateBoxCoxTransformer | apply a learned Box-Cox transformation to a vector
# UnivariateDiscretizer | discretize a `Continuous` vector, and hence render its elscitypw `OrderedFactor`

# In addition to "dynamic" transformers (ones that learn something
# from the data and must be `fit!`) users can wrap ordinary functions
# as transformers, and such *static* transformers can depend on
# parameters, like the dynamic ones. See
# [here](https://alan-turing-institute.github.io/MLJ.jl/dev/transformers/#Static-transformers-1)
# for how to define your own static transformers.

# ### Pipelines

length(schema(Xcont).names)

# Let's suppose that additionally we'd like to reduce the dimension of
# our data.  A model that will do this is `PCA` from
# `MultivariateStats.jl`:

PCA = @load PCA
reducer = PCA()

# Now, rather simply repeating the work-flow above, applying the new
# transformation to `Xcont`, we can combine both the encoding and the
# dimension-reducing models into a single model, known as a
# *pipeline*. While MLJ offers a powerful interface for composing
# models in a variety of ways, we'll stick to these simplest class of
# composite models for now. The simplest way to construct a pipeline
# is using the Julia's `|>` syntax:

pipe = encoder |> reducer

# Notice that the model `pipe` has other models as hyperparameters
# (with names automatically generated based on the mode type
# name). The hyperparameters of the component models are are now
# *nested*, but we can still access them:

@show pipe.pca.pratio
pipe.pca.pratio = 0.85

# The pipeline model behaves like any other transformer:

mach = machine(pipe, X)
fit!(mach)
Xsmall = transform(mach, X)
schema(Xsmall)

# Want to combine this pre-processing with ridge regression?

RidgeRegressor = @load RidgeRegressor pkg=MLJLinearModels
rgs = RidgeRegressor()
pipe2 = pipe |> rgs

# Now our pipeline is a supervised model, instead of a transformer,
# whose performance we can evaluate:

mach = machine(pipe2, X, y)
evaluate!(mach, measure=mae, resampling=Holdout()) # CV(nfolds=6) is default

# ### Training of composite models is "smart"

# Now notice what happens if we train on all the data, then change a
# regressor hyper-parameter and retrain:

fit!(mach)

#-

pipe2.ridge_regressor.lambda = 0.1
fit!(mach)

# Second time only the ridge regressor is retrained!

# Mutate a hyper-parameter of the `PCA` model and every model except
# the `ContinuousEncoder` (which comes before it will be retrained):

pipe2.pca.pratio = 0.9999
fit!(mach)

# ### Inspecting composite models

# The dot syntax used above to change the values of *nested*
# hyper-parameters is also useful when inspecting the learned
# parameters and report generated when training a composite model:

fitted_params(mach).ridge_regressor

#-

report(mach).pca

# ### Incorporating target transformations

# Next, suppose that instead of using the raw `:price` as the training
# target, we want to use the log-price (a common practice in dealing
# with house price data). However, suppose that we still want to
# report final *predictions* on the original linear scale (and use
# these for evaluation purposes). Then we wrap our supervised model
# using `TransformedTargetModel`, which has to key-word arguments
# `target` and `inverse`.

# First we'll overload `log` and `exp` for broadcasting:
Base.log(v::AbstractArray) = log.(v)
Base.exp(v::AbstractArray) = exp.(v)

# Now for the new pipeline:

rgs_log = TransformedTargetModel(rgs, target=log, inverse=exp)

pipe3 = pipe |> rgs_log
mach = machine(pipe3, X, y)
evaluate!(mach, measure=mae)

# MLJ will also allow you to insert *learned* target
# transformations. For example, we might want to apply
# `Standardizer()` to the target, to standardize it, or
# `UnivariateBoxCoxTransformer()` to make it look Gaussian. Then
# instead of specifying a *function* for `target`, we specify a
# unsupervised *model* (or model type). One does not specify `inverse`
# because only models implementing `inverse_transform` are
# allowed.

# Let's see which of these two options results in a better outcome:

box = UnivariateBoxCoxTransformer(n=20)
stand = Standardizer()

rgs_box = TransformedTargetModel(rgs, target=box)
pipe4 = pipe |> rgs_box
mach = machine(pipe4, X, y)
evaluate!(mach, measure=mae)

#-

pipe4.transformed_target_model_deterministic.target = stand
evaluate!(mach, measure=mae)

# ### Resources for Part 3

# - From the MLJ manual:
#     - [Transformers and other unsupervised models](https://alan-turing-institute.github.io/MLJ.jl/dev/transformers/)
#     - [Linear pipelines](https://alan-turing-institute.github.io/MLJ.jl/dev/linear_pipelines/#Linear-Pipelines)
# - From Data Science Tutorials:
#     - [Composing models](https://juliaai.github.io/DataScienceTutorials.jl/getting-started/composing-models/)

# ### Exercises for Part 3

# #### Exercise 7

# Consider again the Horse Colic classification problem considered in
# Exercise 6, but with all features, `Finite` and `Infinite`:

csv_file = urldownload("https://raw.githubusercontent.com/ablaom/"*
                   "MachineLearningInJulia2020/"*
                   "for-MLJ-version-0.16/data/horse.csv");
horse = DataFrames.DataFrame(csv_file); # convert to data frame
coerce!(horse, autotype(horse));
coerce!(horse, Count => Continuous);
coerce!(horse,
        :surgery               => Multiclass,
        :age                   => Multiclass,
        :mucous_membranes      => Multiclass,
        :capillary_refill_time => Multiclass,
        :outcome               => Multiclass,
        :cp_data               => Multiclass);

y, X = unpack(horse, ==(:outcome));
schema(X)

# (a) Define a pipeline that:
# - uses `Standardizer` to ensure that features that are already
#   continuous are centered at zero and have unit variance
# - re-encodes the full set of features as `Continuous`, using
#   `ContinuousEncoder`
# - uses the `KMeans` clustering model from `Clustering.jl`
#   to reduce the dimension of the feature space to `k=10`.
# - trains a `EvoTreeClassifier` (a gradient tree boosting
#   algorithm in `EvoTrees.jl`) on the reduced data, using
#   `nrounds=50` and default values for the other
#    hyper-parameters

# (b) Evaluate the pipeline on all data, using 6-fold cross-validation
# and `cross_entropy` loss.

# ☆(c) Plot a learning curve which examines the effect on this loss
# as the tree booster parameter `max_depth` varies from 2 to 10.

#


text:
"""
Return the gradient of the learnable parameters w.r.t. some objective
"""
function grad end
function grad! end

"""
Proximal operator of a function (https://en.wikipedia.org/wiki/Proximal_operator)
"""
function prox end
function prox! end

"""
Anything that takes an input and performs some kind
of function to produce an output. For example a linear
prediction function.
"""
abstract type Transformation end
abstract type StochasticTransformation <: Transformation end
abstract type Learnable <: Transformation end

"""
Do a forward pass, and return the output
"""
function transform end
function transform! end

"""
Baseclass for any prediction model that can be minimized.
This means that an object of a subclass contains all the
information needed to compute its own current loss.
"""
abstract type Minimizable <: Learnable end

function update end
function update! end
function learn end
function learn! end

# --------------------------------------------------------------------

import Base: AbstractSet

"A continuous range (inclusive) between a lo and a hi"
struct IntervalSet{T} <: AbstractSet{T}
    lo::T
    hi::T
end
function IntervalSet(lo::A, hi::B) where {A,B}
    T = promote_type(A,B)
    IntervalSet{T}(convert(T,lo), convert(T,hi))
end

# numeric interval
randtype(s::IntervalSet{T}) where T <: Number = Float64
Base.rand(s::IntervalSet{T}, dims::Integer...) where T <: Number = rand(dims...) .* (s.hi - s.lo) .+ s.lo
Base.in(x::Number, s::IntervalSet{T}) where T <: Number = s.lo <= x <= s.hi
Base.length(s::IntervalSet{T}) where T <: Number = 1
Base.:(==)(s1::IntervalSet{T}, s2::IntervalSet{T}) where T = s1.lo == s2.lo && s1.hi == s2.hi

# vector of intervals
randtype(s::IntervalSet{T}) where T <: AbstractVector = Vector{Float64}
Base.rand(s::IntervalSet{T}) where T <: AbstractVector = Float64[rand() * (s.hi[i] - s.lo[i]) + s.lo[i] for i=1:length(s)]
Base.in(x::AbstractVector, s::IntervalSet{T}) where T <: AbstractVector = all(i -> s.lo[i] <= x[i] <= s.hi[i], 1:length(s))
Base.length(s::IntervalSet{T}) where T <: AbstractVector = length(s.lo)

"Set of discrete items"
struct DiscreteSet{T<:AbstractArray} <: AbstractSet{T}
    items::T
end
randtype(s::DiscreteSet) = eltype(s.items)
Base.rand(s::DiscreteSet, dims::Integer...) = rand(s.items, dims...)
Base.in(x, s::DiscreteSet) = x in s.items
Base.length(s::DiscreteSet) = length(s.items)
Base.getindex(s::DiscreteSet, i::Int) = s.items[i]
Base.:(==)(s1::DiscreteSet, s2::DiscreteSet) = s1.items == s2.items

# operations on arrays of sets
randtype(sets::AbstractArray{S,N}) where {S <: AbstractSet, N} = Array{promote_type(map(randtype, sets)...), N}
Base.rand(sets::AbstractArray{S}) where S <: AbstractSet = eltype(randtype(sets))[rand(s) for s in sets]
function Base.rand(sets::AbstractArray{S}, dim1::Integer, dims::Integer...) where S <: AbstractSet
    A = Array{randtype(sets)}(undef, dim1, dims...)
    for i in eachindex(A)
        A[i] = rand(sets)
    end
    A
end
function Base.in(xs::AbstractArray, sets::AbstractArray{S}) where S <: AbstractSet
    size(xs) == size(sets) && all(map(in, xs, sets))
end

"Groups several heterogenous sets. Used mainly for proper dispatch."
struct TupleSet{T<:Tuple} <: AbstractSet{T}
    sets::T
end
TupleSet(sets::AbstractSet...) = TupleSet(sets)

# rand can return arrays or tuples, but defaults to arrays
randtype(sets::TupleSet, ::Type{Vector}) = Vector{promote_type(map(randtype, sets.sets)...)}
Base.rand(sets::TupleSet, ::Type{Vector}) = eltype(randtype(sets, Vector))[rand(s) for s in sets.sets]
randtype(sets::TupleSet, ::Type{Tuple}) = Tuple{map(randtype, sets.sets)...}
Base.rand(sets::TupleSet, ::Type{Tuple}) = map(rand, sets.sets)
function Base.rand(sets::TupleSet, ::Type{OT}, dim1::Integer, dims::Integer...) where OT
    A = Array{randtype(sets, OT)}(undef, dim1, dims...)
    for i in eachindex(A)
        A[i] = rand(sets, OT)
    end
    A
end
Base.length(sets::TupleSet) = sum(length(s) for s in sets.sets)
Base.iterate(sets::TupleSet) = iterate(sets.sets)
Base.iterate(sets::TupleSet, i) = iterate(sets.sets, i)

randtype(sets::TupleSet) = randtype(sets, Vector)
Base.rand(sets::TupleSet, dims::Integer...) = rand(sets, Vector, dims...)
Base.in(x, sets::TupleSet) = all(map(in, x, sets.sets))

"Returns an AbstractSet representing valid input values"
function inputdomain end

"Returns an AbstractSet representing valid output/target values"
function targetdomain end


text:
using Pkg
Pkg.update()

ipks = Pkg.installed()
pkgs = ["PyCall","RCall","RDatasets","LibCURL","JavaCall","DistributedArrays","SymPy"]
for p in pkgs
    !haskey(ipks,p) && Pkg.add(p)
end

using PyCall, PyPlot, Plots, RCall, RDatasets, JavaCall
using DistributedArrays, SymPy, LibCURL
using Printf, Random, Distributed

# Simple call to a C-routine
# Can be used to set the random number seed using the system clock

systime()   = ccall((:clock,"libc"),Int32,())
randomize() = Random.seed!(systime())

systime()

# A 'mad' example (a*b + c)
# The routine is called fma in library: libc
#
mad(a,b,c) = ccall((:fma,"libc"),Float64,(Float64,Float64,Float64),a,b,c)
mad(3.1,5.2,7.4)

# Use the system library to generate some random numbers
[ccall((:rand, "libc"), Int32, ()) for i = 1:6]

# Call a FORTRAN routine from LAPACK to compute the dot product
# between two arrays.
# FORTRAN passes scalar arguments by reference
# Note: the vectors arepassed by reference already.

function compute_dot(DX::Vector{Float64},DY::Vector{Float64})
    @assert length(DX) == length(DY)
    n = length(DX)
    incx = incy = 1
    dotprod = ccall((:ddot, "libLAPACK"),
                   Float64,
                   (Ptr{Int64}, Ptr{Float64}, Ptr{Int64},
                   Ptr{Float64}, Ptr{Int64}),
                   Ref(n), DX, Ref(incx), DY, Ref(incy))
    return dotprod
end

# Test it out
aa = [rand() for i = 1:1000]
bb = [rand() for i = 1:1000]
compute_dot(aa,bb)    # => Should be about 250

# Get the user's home directory ENV["HOME"]
ptr = ccall((:getenv, "libc"), Ptr{UInt8},(Ptr{UInt8},),"HOME")

# My modules are in Julia/MyModules - put this on the LOAD PATH
#
myHome = unsafe_string(ptr)
push!(LOAD_PATH,string(myHome,"/Julia/MyModules"))

#=
The system library call has a third parameter, which when set to zero will create a new variable but not overwrite an existing one.
=#
# Define the function to replace existing variable
evset(var::String, val::String) =
     ccall((:setenv,"libc"),Clong,
           (Cstring,Cstring,Clong),var,val,1);

# The unset routine is quite simple
evunset(evvar::String) =
     ccall((:unsetenv,"libc"),Clong,(Cstring,),evvar);

#=
The system library call has a third parameter, which when set to zero will create a new variable but not overwrite an existing one.
=#
# Define the function to replace existing variable
evset(var::String, val::String) =
     ccall((:setenv,"libc"),Clong,
           (Cstring,Cstring,Clong),var,val,1);

# The unset routine is quite simple
evunset(evvar::String) =
     ccall((:unsetenv,"libc"),Clong,(Cstring,),evvar);

# Set an environment variable PACKT_HOME ...
# ... and check it
packt = evset("PACKT_HOME", string(myHome,"/PacktPub"));

ENV["PACKT_HOME"]

# Now unset it, verify it is so.
evunset("PACKT_HOME");
ENV["PACKT_HOME"]

#=

// Basel function in C

#include<stdio.h>
#include<stdlib.h>

double basel(int N) {

  double s = 0.0L;
  int i;
  double x;

  if (N < 1) return s;
  for (i = 1; i <= N; i++) {
    x = 1.0L/((double) i);
    s += x*x;
  }
  return s;
}

// Horner's method

#include<math.h>

double horner(double x, double aa[], long n) {
  long i;
  double s = aa[n-1];
  if (n > 1) {
    for (i = n-2; i >= 0; i--) {
      s = s*x + aa[i];
    }
  }
  return s;
}

// Build a dynamic library (on OSX) as:
//
// clang -c basel.c horner.c
// libtool -dynamic basel.o horner.o -o libmyfuns.dylib
//         -lSystem -macosx_version_min 10.13
//
// sudo cp libmyfuns.dylib /usr/local/lib

=#

run(`nm -g libmyfuns.dylib`)

basel = ccall((:basel,"libmyfuns"),Float64,(Int64,),10000000)

# Time the function for 10^7 loops
@elapsed ccall((:basel,"libmyfuns"),Float64,(Int64,),10000000)

# C version of Horners method
# Note use of Cdouble, Clong rather than Float64, Int64 etc.
#
x = 2.1;
aa = [1.0, 2.0, 3.0, 4.0, 5.0];

ccall((:horner,"libmyfuns.dylib"),Cdouble,
              (Cdouble,Ptr{Cdouble},Clong),x,aa,length(aa))

# Compute PI in 'C', usual Monte Carlo method
# We ca put the C code in a multiline string

C_code = """
#include <stddef.h>
#include <stdlib.h>

double c_pi(long n) {
    long k = 0L;
    float rmax = (float) RAND_MAX;
    for (long i = 0L; i < n; i++) {
        float x = ((float) rand())/rmax;
        float y = ((float) rand())/rmax;
        if ((x*x + y*y) < 1.0) {
          k++;
        }
    }
    return 4.0*((double) k)/((double) n);
}
"""

# Get a temporary name to create a library
const Clib = tempname()   # ... make a temporary file

# compile to a shared library by piping C_code to gcc
# (works only if you have gcc installed):
using Libdl
tmplib = string(Clib,".",dlext)
open(`gcc -fPIC -O3 -msse3 -xc -shared -o $tmplib -`, "w") do f
    print(f, C_code)
end

# define a Julia function that calls the C function:
c_pi(N::Int64) = ccall(("c_pi", Clib), Float64, (Clong,), N)

using Random
randomize();
c_pi(1000000)

C_Code = """
include <julia.h>
#include <stdio.h>
#include <math.h>

// Only define this once if in an executable ...
// (i.e. not in a shared library) ...
// ... and if we want the fastest code.
JULIA_DEFINE_FAST_TLS()

int main(int argc, char *argv[])
{
/* required: setup the Julia context */
  jl_init();

/* run Julia commands */
  jl_function_t *fnc1 = jl_get_function(jl_base_module, "exp");
  jl_function_t *fnc2 = jl_get_function(jl_base_module, "sin");
  jl_value_t* arg1 = jl_box_float64(-0.3);
  jl_value_t* arg2 = jl_box_float64(3.0);
  jl_value_t* ret1 = jl_call1(fnc1, arg1);
  jl_value_t* ret2 = jl_call1(fnc2, arg2);

/* unbox and setup final result */
  double retD1 = jl_unbox_float64(ret1);
  double retD2 = jl_unbox_float64(ret2);
  double retD3 = retD1*retD2;
  printf("sin(3.0)*exp(-0.3) from Julia API: %e\n", retD3);
  fflush(stdout);

/* Allow Julia time to cleanup pending write requests etc. */
  jl_atexit_hook(0);
  return 0;
}
""";

# Can use a nice command script in 'share/julia/julia-config.jl
./julia-config.jl
usage: julia-config [--cflags | --ldflags |
                     --ldlibs | --allflags]

JULIA_HOME = \
/Applications/Julia-1.1.app/Contents/Resources/julia; export JULIA_HOME

ls $JULIA_HOME
LICENSE.md bin etc include lib share

cc jltest.c -o jltest -std=gnu99 \
      -I$JULIA_HOME/include/julia \
      -DJULIA_ENABLE_THREADING=1 -fPIC \
      -L$JULIA_HOME/lib \
      -Wl,-rpath,$JULIA_HOME/lib \
      -Wl,-rpath,$JULIA_HOME/lib/julia \
      -ljulia

./jltest
sin(3.0)*exp(-0.3) from Julia API: 1.045443e-01

using LibCURL
url = "http://LondonJulia.org/mastering-julia.html"

# init a curl handle
curl = curl_easy_init();

# set the URL and request to follow redirects
curl_easy_setopt(curl, CURLOPT_URL, url);
curl_easy_setopt(curl, CURLOPT_FOLLOWLOCATION, 1);

# setup the callback function to recv data
function curl_write_cb(curlbuf::Ptr{Nothing}, s::Csize_t, n::Csize_t, p_ctxt::Ptr{Nothing})
  sz = s * n
  data = Array{UInt8}(undef,sz)
  ccall(:memcpy, Ptr{Nothing}, (Ptr{Nothing}, Ptr{Nothing}, UInt64), data, curlbuf, sz)
  println(String(data))
  sz::Csize_t
end

c_curl_write_cb =
  @cfunction(curl_write_cb, Csize_t,
             (Ptr{Nothing}, Csize_t, Csize_t, Ptr{Nothing}));
curl_easy_setopt(curl, CURLOPT_WRITEFUNCTION, c_curl_write_cb);

# execute the query
res = curl_easy_perform(curl)

println("curl url exec response : ", res)

# retrieve HTTP code
http_code = Array{Clong}(undef,1)
curl_easy_getinfo(curl, CURLINFO_RESPONSE_CODE, http_code)
println("httpcode : ", http_code)

# release handle
curl_easy_cleanup(curl)

using PyCall

@pyimport numpy.random as nr
aa = nr.rand(4,5)      # aa is a Julia array generated with Numpy

# So an array slice has no overhead(i.e. in Julia)
aa[2:3,2:3]

# Using the SciPy package
@pyimport scipy.optimize as so
so.ridder(x -> x*cos(x),1,π)

@pyimport scipy.integrate as si
si.quad(x -> x*sin(x),1,π)

# ... and plotting with matplotlib
# @pyimport matplotlib.pyplot as plt
x = range(0,stop=10,length=1000);
y = sin.(3*x + 4*cos.(2*x));

plt.plot(x, y, color="red", linewidth=2.0, linestyle="--")

plt.show()

using SymPy

# We need to work with a special type: a symbol (in the SymPy sense)
u = symbols("u")
x = symbols("x", real=true)
y1, y2 = symbols("y1, y2", positive=true)
alpha = symbols("alpha", integer=true, positive=true)

typeof(x)

# We can solve some algebraic equations
solve(u^2 + 1)

# Perform functional expansions
p = expand(prod([sin(x^(-i)) for i in 1.0:1.0:5.0]))

# SymPy can use the Plot API to display graphs
x = symbols("x")
using Plots
pyplot()

Plots.plot(p, 0.0, 1.0)

# For a more complex example, which we have seen before, ...
# ... the sum of the Basel series
#
i, n = symbols("i, n")
sn = sum(1/i^2, (i, 1, n))

limit(sn, n, oo)

# An initial valued problem
#
y = SymFunction("y")
a, x = symbols("a,x")
eqn = y'(x) - 3*x*y(x) - 1

using RCall

# Perform a t-test on a set of 1000 normally distributed variates
#
x = randn(1000);
R"t.test($x)"

# Optimise function using BFGS (Broyden–Fletcher–Goldfarb–Shannon§§§§§≠±±) algorithm
#
f(x) = 10*sin(0.3*x)*sin(1.3*x^2) + 0.00001*x^4 + 0.2*x + 40
R"optim(0, $(x -> f(x)), method='BFGS')"

@rlibrary MASS

 geyser = rcopy(R"MASS::geyser")

round(sum(geyser[:waiting])/sum(geyser[:duration]),digits=5)

rcall(:summary, geyser[:waiting])

rcall(:summary, geyser[:duration])

# Use R to grab some financial data
# and plot it.
using RCall

R"""
library(data.table)
library(scales)
library(ggplot2)

link <- "Files/fin_data.csv"
dt <- data.table(read.csv(link))
dt[, date := as.Date(date)]

# create indexed values
dt[, idx_price := price/price[1], by = ticker]
"""

R"""
ggplot(dt, aes(x = date, y = idx_price, color = ticker)) +
       geom_line() + theme_bw() +
       xlab("Date") + ylab("Pricen(Indexed 2000 = 1)") +
       scale_color_discrete(name = "Company")
"""

using JavaCall
JavaCall.init(["-Xmx128M"])     # Initialise and ask for additional memory

#=
 May set the current directory as the classpath in the init() call
 viz. JavaCall.init(["-Xmx128M -Djava.class.path=$(@__DIR__)"])

 Ignore the (possible) segmentation fault
 See: http://juliainterop.github.io/JavaCall.jl/faq.html

 This works inside in a Notebook, despite the segmentation fault (on OSX)
 In the REPL, julia may nned to be started with a --handle-signals=no
 option to disable Julia's signal handler.
 [This may cause issues with handling ^C in Julia programs.]
=#

#=

 Ignore the Segmentation fault on OSX, as explained above

=#

# JavaCall needs to modify Julia variables before using them in functions
a = JString("Hello, Blue Eyes")

# The JString as a single field (:ptr)
fieldnames(typeof(a))

a.ptr

# The pointer can now be used in a Java function

b = ccall(JavaCall.jnifunc.GetStringUTFChars, Ptr{UInt8},
            (Ptr{JavaCall.JNIEnv}, Ptr{Nothing}, Ptr{Nothing}),
                                 JavaCall.penv, a.ptr, C_NULL)

unsafe_string(b)

# Similar to PyCall, we can use JavaCall to perform some mathematics

jlm = @jimport "java.lang.Math"
jcall(jlm, "exp", jdouble, (jdouble,), pi*pi/6.0)

# ... and to access Internet libraries
jnu = @jimport java.net.URL
jurl = jnu((JString,), "http://LondonJulia.org/mastering-julia.html")
www = jcall(jurl, "getHost", JString,())

# Finally let's import the ArrayList class, ...
# ... define an instance and add some items to it
JArrayList = @jimport(java.util.ArrayList)
a = JArrayList(())
jcall(a, "add", jboolean, (JObject,), "'Twas ")
jcall(a, "add", jboolean, (JObject,), "brillig, ")
jcall(a, "add", jboolean, (JObject,), "and ")
jcall(a, "add", jboolean, (JObject,), "the ")
jcall(a, "add", jboolean, (JObject,), "slithy ")
jcall(a, "add", jboolean, (JObject,), "toves,")

# Now iterate thru' the array and push it on a Julia array ...
# ... converting the bit type to an (unsafe) string.
#
t = Array{Any, 1}()
for i in JavaCall.iterator(a)
  push!(t, unsafe_string(i))
end

# Evaluate the  result.
join(t)

# Get a webpage using curl or wget
# Notice the backticks
# First check that curl (or wget) is available, ...
# ... otherwise you will need to install it.
#
cmd = `which curl`
typeof(cmd)

#=
 Commands now need to be run, i.e. they are not executed immediately
 Since commands run as separate tasks, it is usually preferable to
 suppress the output of the run() function

 The task output will goto STDOUT (by default) but can be captured
 and used by the Julia program
=#
run(cmd)

# We can use curl to get the webpage ...
# .. and can do it in one step
proc = run(`curl "http://LondonJulia.org/mastering-julia.html"`);

# Look at Base:process.jl for definition of Julia process structure

#=
mutable struct Process <: AbstractPipe
    cmd::Cmd
    handle::Ptr{Cvoid}
    in::IO
    out::IO
    err::IO
    exitcode::Int64
    termsignal::Int32
    exitnotify::Condition
    closenotify::Condition
    . . .
    . . .
end
=#

proc.exitcode   # Zero is normal exit status

cd(string(ENV["HOME"],"/PacktPub/Datasets/Alice"))
pwd()

# Needs to be where the chapter 5 files are and may be changed depending
# on the user's setup

# Define a function to use the wordcount ('wc') utility
wc(f) = isfile(f) && run(`wc $f`)

# Count the number of occurences of 'beaver" in the Hunting of the Snark,
# piping the output to a file

txtfile = "hunting-the-snark.txt";
logfile = "hunting-the-snark.log";

run(pipeline(`grep -i beaver $txtfile`,stdout=logfile));

wc(logfile);

# Do this again for the bellman but append the output to the log file.

run(pipeline(`grep -i bellman $txtfile`,stdout=logfile,append=true));
wc(logfile);

#=
  Do we have any lines with both beaver and bellman
  Note in the 'grep' we ignored case (-i) but will stiil find
  plurals and possibly punctuation, viz.: Beavers.
  I'll attend to this later.
=#
run(pipeline(`grep -bellman $txtfile`,`grep -i beaver`));

function ftop(dir=".", nf=20)
  @assert nf > 0
  efind =  `find $dir -type f -iname "*" -exec du -sh "{}" + `
  try
    run(pipeline(efind,`sort -rh`,`head -$nf`))
    println("Done.")
  catch
    println("No files found.")
  end
end

dd = ENV["HOME"]*"/PacktPub"
ftop(dd,10)

# Code to run from the command line
# Source: ftop-main.jl

#=
# See if any arguments have been passed
# If so process them otherwise set the defaults
#
nargs = size(ARGS)[1]
if nargs > 2
  println("usage: ftop [dir [nfiles]]")
  exit(-1)
else
  dir = (nargs > 0) ? ARGS[1] : "."
  nf  = (nargs > 1) ? ARGS[2] : 20
end

# Main 'find' command
efind =  `find $dir -type f -iname "*" -exec du -sh "{}" + `

# Get and print the directory to be searched
cwd = string(pwd(),"/",dir)
println("Directory: $cwd")

# Then run the pipeline
try
  run(pipeline(efind,`sort -rh`,`head -$nf`))
  println("Done.")
catch
  println("No files found.")
end
=#

#=
 Perl has falled out of fashion with the rise of Python but still
 but still remains one of the best methods for data munging.

  Unix distros and OSX (normally) have Perl avaiable but in Windows
  it needs to installed and on the executable path.

  Julia performance in handling string is not one of its greatest strenghs
  so munging large files can successfully done using Perl.

  Note:
  Julia has introduced an analytical engine (JuliaDB) to tackle the
  processing of large datasets, employing some clever memory management
  techniques and we will discuss this in the next chapter

=#

# On OSX and Linux, there is a word list in /usr/share/dict
# This command outputs an palindromes of 6 or more letters
# (I am not try to teach Perl but most of the syntax is easily followed)
cmd = `perl -nle 'print if $_ eq reverse && length > 5' /usr/share/dict/words`
run(cmd);

# Here is a work-around using temporary files
# Capture the output of palindrome command and

tmpfile = mktemp(tempdir());   # Alternative to tempfile() call
cmd = `perl -nle 'print if $_ eq reverse && length > 5' /usr/share/dict/words`
run(pipeline(cmd; stdout=tmpfile));

dmp = read(tmpfile);
run(`rm -f $tmpfile`);  # Good idea to remove the temporary file

# Dump is a byte array
(eltype(dmp),length(dmp))

# Convert the byte array to a string
# Note the carriage returns
ss = String(dmp)

# Remove the trailing \n and then split into words
sa = split(chomp(ss),"\n")
typeof(sa)

# PERL is good (and quick) for processing strings
# A large pipeline (probably not on Windows)
# Get top 10 words in the Hunting of the Snark, ignoring blank lines
# I'll leave it the the reader ti work out what the command do!
#
cd(ENV["HOME"]*"/PacktPub/Datasets/Alice")

fl = "hunting-the-snark.txt";
c1 = `perl -pne 'tr/[A-Z]/[a-z]/' $(fl)`;
c2 = `perl -ne 'print join("\n", split(/\s+/,$_));print("\n")'`;

run(pipeline(c1,c2,`sort`,`grep -ve '^$'`,`uniq -c`,`sort -rn`,`head -10`));

# Look at occurences of Bellman
#
run(pipeline(c1,c2,`sort`,`grep -ve '^$'`,`uniq -c`,`sort -rn`,`grep -i Bellman`));

# Suspect there is a problem with punctuation marks
# Here are the bottom 10 words.
#
run(pipeline(c1,c2,`sort`,`grep -ve '^$'`,`uniq -c`,`sort -rn`,`tail -10`));

# So we need to add an extra task in the pipe
# If we do this through a function we can pass any poem
#
function munge(fl)
  c1 = `perl -pne 'tr/[A-Z]/[a-z]/' $(fl)`;
  c2 = `perl -pne 's/()[[:punct:]]//g'`;
  c3 = `perl -ne 'print join("\n", split(/\s+/,$_));print("\n")'`
  c4 = `grep -ve '^$'`
  read(pipeline(c1,c2,c3,`sort`,c4,`uniq -c`,`sort -rn` ))
end

text = munge("hunting-the-snark.txt");
(eltype(text),length(text))

lines = split(String(text),"\n")
n = length(lines)

s2 = [split(lines[i]) for i = n-9:n]

pwd()

# We can use the process in/out channels to capture the I/O
# (rev.pl is in the parent folder)
jabber = "jabberwocky.txt";
proc = open(`rev.pl $jabber`,"r+");
close(proc.in);

proc.out

poem = readlines(proc.out);
close(proc.out)

poem

# We are not limited to just Perl(5)
# Get Perl6 from https://rakudo.org/files/
# Put it on the path OR setup a symbolic link to the binary:
# viz: perl6 -> /Applications/Rakudo/bin/perl6

run(`which perl6`);

# Use Perl6 to find the line of the greatest length
#
run(`perl6 -e 'my $max="";
      for (lines) {$max = $_ if .chars > $max.chars};
      END { $max.say }' hunting-the-snark.txt`);

# These techniques can be used with any command processor
# This turns all characters after 40th spot to RED
#
run(`ruby -e 'w = $*.shift; $<.each { |l| puts "#{l}\e[31m#{l.chop!.slice!(w.to_i..-1)}\e[0m" }' 40 hunting-the-snark.txt`);

# And also with Python
# Encode the jabberwocky as Base 64

cd(ENV["HOME"]*"/PacktPub/Datasets/Alice")
f1="jabberwocky.txt"
f2="jabberwocky.b64"
b64encode = `python -c 'import base64,sys; base64.encode(open(sys.argv[1],"rb"),open(sys.argv[2],"wb"))' $f1 $f2`
run(b64encode);

# Check it
run(`cat $f2`);

# Now reverse the process
#
run(`base64 --decode $f2`);

cd(string(ENV["HOME"],"/PacktPub/Chp05")); # You may have to change this

# I got an error after exiting the following process
# So have wrapped it in a try/catch block
try run(`wc $(readdir())`) catch end

# Define a macro to trap any run(CMD) errors
macro traprun(c)
  quote
    if typeof($(esc(c))) == Cmd
      try
        run($(esc(c)))
     catch
     end
    end
  end
end

# Filter on a regular expression
# This will get rid of the directory warnings above
function filter(pat::Regex, dir=".")
  a = Any[]
  for f in readdir(dir)
    occursin(pat,f) && push!(a, f)
  end
  return a
end

# ... and find all the files ending in '.txt'

@traprun `wc $(filter(r"\.txt$"))`;

# Find all the Jupyter notebooks
cd("..")
@traprun `find "/Users/malcolm/PacktPub" -name Chp\*.ipynb`;

# Location of an Apache access log

cd(homedir()*"/PacktPub/Chp05")
logf = "Files/access_log";
run(`/usr/bin/wc -l $logf`);

# First 5 lines - using head command
# Can use the '$' prefix inside the backticks
run(`head -5 $logf`)

# Using PERL this is not a problem
# Do a "tail -5" as a one-liner to show last 5 lines

cmd = `perl -ne 'push @a, $_; @a = @a[@a-5..$#a]; END { print @a }' $logf`
run(cmd);

;cat Code/hcount.pl

# Show the top 10 IP addresses in the weblog

println("Top Ten Hitters ($logf)");
@time run(`Code/hcount.pl $logf`);


text:
genepop(normpath(joinpath(@__DIR__, "..", "precompile", "precompile")) * ".gen", silent = true);
structure(normpath(joinpath(@__DIR__, "..", "precompile", "precompile")) * ".str", silent = true);
delimited(normpath(joinpath(@__DIR__, "..", "precompile", "precompile")) * ".csv", silent = true);
vcf(normpath(joinpath(@__DIR__, "..", "precompile", "precompile")) * ".vcf.gz", silent = true) ;