text:
using Base: has_fast_linear_indexing, first_step_last_ascending
include("mydct.jl")

using MAT
using LinearAlgebra

matrix = [231 32 233 161 24 71 140 245
247 40 248 245 124 204 36 107
234 202 245 167 9 217 239 173
193 190 100 167 43 180 8 70
11 24 210 177 81 243 8 112
97 195 203 47 125 114 165 181
193 70 174 167 41 30 127 245
87 149 57 192 65 129 178 228]

MHmatrix = [
1.11e+03    4.40e+01    7.59e+01    -1.38e+02   3.50e+00    1.22e+02   1.95e+02    -1.01e+02
7.71e+01    1.14e+02    -2.18e+01   4.13e+01    8.77e+00    9.90e+01   1.38e+02    1.09e+01
4.48e+01    -6.27e+01   1.11e+02    -7.63e+01   1.24e+02    9.55e+01   -3.98e+01   5.85e+01
-6.99e+0    -4.02e+01   -2.34e+01   -7.67e+01   2.66e+01    -3.68e+01   6.61e+01    1.25e+02
-1.09e+0    -4.33e+01   -5.55e+01   8.17e+00    3.02e+01    -2.86e+01   2.44e+00    -9.41e+01
-5.38e+0    5.66e+01    1.73e+02    -3.54e+01   3.23e+01    3.34e+01   -5.81e+01   1.90e+01
7.88e+01    -6.45e+01   1.18e+02    -1.50e+01   -1.37e+02   -3.06e+01   -1.05e+02   3.98e+01
1.97e+01    -7.81e+01   9.72e-01    -7.23e+01   -2.15e+01   8.13e+01   6.37e+01    5.90e+00
]

cmatrix = [1118.7500000000002 44.02219262301787 75.91905032609426 -138.57241099707306 3.5000000000000284 122.07805518803062 195.04386762362958 -101.60490592795372; 77.1900790187901 114.86820590697009 -21.80144211456541 41.364135057036215 8.77720597607919 99.08296204231638 138.17151573211936 10.909279534077866; 44.83515365361005 -62.75244638273558 111.61411421683839 -76.37896579814382 124.42215968070717 95.59841936381827 -39.82879694884997 58.52376699413914; -69.98366474768355 -40.240894480279415 -23.497050834570445 -76.73205936982573 26.645774950452864 -36.83282895356625 66.18914845485334 125.42973061799121; -109.00000000000003 -43.34308566057378 -55.54369078865144 8.173470828327527 30.25000000000001 -28.660243733699012 2.4414982233614957 -94.14370254649862; -5.387835905170836 56.63450089125683 173.0215190411498 -35.42344938713219 32.38782492363563 33.45767275218619 -58.1167863718725 19.022561487011163; 78.84396931190855 -64.592409552383 118.67120305115003 -15.09048397585811 -137.31692787267204 -30.619666281198647 -105.11411421683843 39.813049706882936; 19.78824382848893 -78.18134089948808 0.9723118598351412 -72.34641800960696 -21.578163250360877 81.29990354778818 63.71037820527631 5.90618071066943]

v = [231 32 233 161 24 71 140 245]
c = [401.99 6.60002 109.167 -112.786 65.4074 121.831 116.656 28.8004]

function equalsMat(mat, emat)
    if round.(mat, digits=2)==round.(emat, digits=2)
        return true
    else
        return false
    end

end

function testDCT2(mat, cmat)
    if equalsMat(mDCT2n3(mat), cmat)
        return true
    else
        return false
    end
end

function testDCT2(mat)
    if equalsMat(mDCT2n3(mat), dct(mat))
        return true
    else
        return false
    end
end

function testiDCT2(cmat, mat)
    if equalsMat(miDCT2n3(cmat), mat)
        return true
    else
        return false
    end
end

function testiDCT2(cmat)
    if equalsMat(miDCT2n3(cmat), idct(cmat))
        return true
    else
        return false
    end
end

function testDCT(v, c)
    if equalsMat(mDCT(v), c)
        if equalsMat(v, miDCT(c))
            return true
        end
    else
        return false
    end
end

## VENGONO TESTATE TUTTE LE DCT

function testALL(mat, cmat)
    if (testDCT2(mat, cmat) && (testDCT2(mat))
        && testiDCT2(cmat, mat) && testiDCT2(cmat))
        return true
    else
        return false
    end
end

function testALL(mat)
    cmat = dct(mat)
    imat = idct(mat)
    return testALL(mat, cmat) && testALL(imat, mat)
end

if testALL(v,c)
    println("All vector test PASSED")
else
    println("All vector test NOT Passed")
end
if testALL(matrix, cmatrix)
    println("All matrix test PASSED")
else
    println("All matrix test NOT Passed")
end

max_stars_repo_name: illoxian/bitmap-compressor-tool

max_issues_repo_name: illoxian/bitmap-compressor-tool

max_forks_repo_name: illoxian/bitmap-compressor-tool


text:
using DifferentialEquations, Plots, CellMLToolkit, Sundials, TimeseriesPrediction, DynamicalSystems

model_root = joinpath(splitdir(pathof(CellMLToolkit))[1], "..", "models")
ml = CellModel(joinpath(model_root, "ohara_rudy_cipa_v1_2017.cellml.xml"));

tspan = (0, 15000.0)
prob = ODEProblem(ml, tspan);

p = prob.p; p[201] = 175.0; prob = remake(prob, p=p)
tt = collect(tspan[begin]:1.0:tspan[end])

sol = solve(prob, CVODE_BDF(linear_solver=:GMRES), dtmax=0.5, saveat=tt)
plot(sol, vars=49, linewidth=3)

data = sol[49,:]

train_inds = 1:10000
testy_inds = 10001:length(tt)

trainData = data[train_inds] .+ 1e-0 .* randn(Float64, length(train_inds))

#theiler = estimate_delay(trainData, "mi_min") # estimate a Theiler window
#Tmax = 200 # maximum possible delay
#Y, τ_vals, ts_vals, Ls, εs = pecuzal_embedding(trainData; τs = 0:1:Tmax , w = theiler, econ = true)

τ_vals = 0:1:10
d = length(τ_vals)
ntype = FixedMassNeighborhood(3)
pred = localmodel_tsp(trainData, d, τ_vals, length(testy_inds); ntype=ntype)

plot!(tt[train_inds], trainData, linewidth=2, label = "training (trunc.)")
plot!(tt[testy_inds], data[testy_inds], linewidth=2, label = "actual signal")
plot!(tt[testy_inds], pred[2:end], color=:red, linestyle=:dash, label="predicted")
plot!(title = "Training points: $(length(train_inds)), predicted points: $(length(testy_inds))", xlabel="\$ t \$", ylabel = "\$ V(t) \$")
plot!(xlims=[tt[train_inds[end-1000]], tt[end]])
savefig("./testing.pdf")

ext: jl

max_stars_repo_name: cmarcotte/timeseriesprediction

max_issues_repo_name: cmarcotte/timeseriesprediction

max_forks_repo_name: cmarcotte/timeseriesprediction


text:
"""
    Contact{T,N}

    Abstract type contain contact information associated with a parent Body object.
"""
abstract type Contact{T,N} end

# constraints
constraint(model::Contact, body::Body, λ, timestep) = constraint(model, next_configuration(body.state, timestep)..., λ)

# constraint Jacobians
constraint_jacobian_configuration(model::Contact, body::Body, id, λ, timestep) = constraint_jacobian_configuration(model, next_configuration(body.state, timestep)..., current_configuration_velocity(body.state)..., λ, timestep)
constraint_jacobian_velocity(model::Contact, body::Body, id, λ, timestep) = constraint_jacobian_velocity(model, next_configuration(body.state, timestep)..., current_configuration_velocity(body.state)..., λ, timestep)

# impulses
impulse_map(model::Contact, body::Body, id, λ, timestep) = impulse_map(model, next_configuration(body.state, timestep)..., λ)

function impulse_map(model::Contact, x::AbstractVector, q::Quaternion, λ)
    X = force_mapping(model, x, q)
    # Q = - X * q * skew(model.contact_point - vector_rotate(model.offset, inv(q)))
    Q = - X * rotation_matrix(q) * skew(model.contact_point - vector_rotate(model.offset, inv(q)))
    return [X'; Q']
end

# force mapping
function force_mapping(model::Contact, x::AbstractVector, q::Quaternion)
    X = [model.surface_normal_projector;
         szeros(1,3);
         model.surface_projector]
    return X
end

# utilities
Base.length(model::Contact{T,N}) where {T,N} = N
neutral_vector(model::Contact{T,N}) where {T,N} = sones(T, Int(N/2))
cone_degree(model::Contact{T,N}) where {T,N} = Int(N/2)


text:
using IndirectArrays, MappedArrays, OrderedCollections
using Test, FixedPointNumbers, Colors

@testset "values::AbstractVector" begin
    colors = [RGB(1,0,0) RGB(0,1,0);
              RGB(0,0,1) RGB(1,0,0)]
    index0 = [1 3;
              2 1]
    for indexT in (Int8, Int16, UInt8, UInt16)
        A = IndirectArray{indexT}(colors)
        @test eltype(A) == RGB{N0f8}
        @test size(A) == (2,2)
        @test ndims(A) == 2
        @test A[1,1] === A[1] === RGB(1,0,0)
        @test A[2,1] === A[2] === RGB(0,0,1)
        @test A[1,2] === A[3] === RGB(0,1,0)
        @test A[2,2] === A[4] === RGB(1,0,0)
        @test isa(eachindex(A), AbstractUnitRange)
        @test A.index == index0
    end
    x = IndirectArray(colors[:])
    @test x == IndirectArray{UInt8}(colors[:])
    xc = copy(x)
    @test x == xc
    xc[2], xc[3] = RGB(0,1,0), RGB(0,1,1)
    @test xc == IndirectArray([RGB(1,0,0), RGB(0,1,0), RGB(0,1,1), RGB(1,0,0)])
    @test append!(x, x) == IndirectArray([colors[:]; colors[:]])
    @test append!(x, IndirectArray([RGB(1,0,0), RGB(1,1,0)])) ==
        IndirectArray([colors[:]; colors[:]; [RGB(1,0,0), RGB(1,1,0)]])
    # Append with non-IndirectArray
    @test append!(IndirectArray(colors[:]), IndirectArray(colors[:])) ==
          append!(IndirectArray(colors[:]), colors[:])

    # Bounds checking upon construction
    index_ob = copy(index0)
    index_ob[1] = 5   # out-of-bounds
    unsafe_ia(idx, vals) = (@inbounds ret = IndirectArray(idx, vals); ret)
      safe_ia(idx, vals) = (ret = IndirectArray(idx, vals); ret)
    @test_throws BoundsError safe_ia(index_ob, colors[1:3])
    # This requires inlining, which means it fails on Travis since we turn
    # off inlining for better coverage stats
    # B = unsafe_ia(index_ob, colors)
    # @test_throws BoundsError B[1]
    # @test B[2] == RGB(0,0,1)

    # Non-Arrays
    a = [0.1  0.4;
        0.33 1.0]
    f(x) = round(Int, 99*x) + 1   # maps 0-1 to 1-100
    m = mappedarray(f, a)
    cmap = colormap("RdBu", 100)
    img = IndirectArray(m, cmap)
    @test img == [cmap[11] cmap[41];
                  cmap[34] cmap[100]]
end

@testset "values::AbstractDict" begin
    # With Dicts
    a = ['a' 'b';
        'q' 'j']
    v = freeze(Dict('a' => "apple", 'b' => "book", 'q' => "quit", 'j' => "jolly"))
    A = IndirectArray(a, v)
    @test A[1,1] == "apple"
    @test A[2,1] == "quit"
    @test A[1,2] == "book"
    @test A[2,2] == "jolly"
    @test_throws BoundsError IndirectArray(a, Dict('a' => "apple", 'b' => "book"))
end


text:
import RIOPA
import RIOPA: Ctrl, DataGen
import Test: @testset, @test, @test_throws

struct TestDataGenTag <: RIOPA.DataGen.DataGenTag end

@testset "DataGen" begin
    @test DataGen.get_tag(nothing) == DataGen.DefaultDataGenTag()
    @test_throws KeyError DataGen.get_tag("test")
    DataGen.add("test", TestDataGenTag())
    @test DataGen.get_tag("test") == TestDataGenTag()

    config = RIOPA.default_config()
    ds = Ctrl.configure_dataset(config[:datasets][1])
    str_cfg = ds.cfg.streams[1]
    @test RIOPA.get_payload_group_id(0, 4, str_cfg) == 1
    @test RIOPA.get_payload_group_id(1, 4, str_cfg) == 1
    @test RIOPA.get_payload_group_id(2, 4, str_cfg) == 2
    @test RIOPA.get_payload_group_id(3, 4, str_cfg) == 2

    @test RIOPA.get_payload_group_id(0, 4, ds.cfg.streams[2]) == 1
    @test RIOPA.get_payload_group_id(1, 4, ds.cfg.streams[2]) == 2
    @test RIOPA.get_payload_group_id(2, 4, ds.cfg.streams[2]) == 2
    @test RIOPA.get_payload_group_id(3, 4, ds.cfg.streams[2]) == 2

    @test RIOPA.get_payload_group_id(0, 5, ds.cfg.streams[1]) == 1
    @test RIOPA.get_payload_group_id(1, 5, ds.cfg.streams[1]) == 1
    @test RIOPA.get_payload_group_id(2, 5, ds.cfg.streams[1]) == 2
    @test RIOPA.get_payload_group_id(3, 5, ds.cfg.streams[1]) == 2
    @test RIOPA.get_payload_group_id(4, 5, ds.cfg.streams[1]) == 2

    @test RIOPA.get_payload_group_id(0, 1, str_cfg) == 2

    DataGen.generate!(DataGen.DefaultDataGenTag(), ds)
    @test 2000 <= length(ds.streams[1].data.vec) <= 2400
    @test 4000 <= length(ds.streams[2].data.vec) <= 4800
    ds.curr_step += 1
    DataGen.generate!(DataGen.DefaultDataGenTag(), ds)
    @test 4000 <= length(ds.streams[1].data.vec) <= 4800
    @test 6000 <= length(ds.streams[2].data.vec) <= 7200

    config[:datasets][1][:data_streams][1][:proc_payloads][1][:ratio] = 1.0
    @test_throws DataGen.ProcessPayloadRatioError Ctrl.configure_dataset(
        config[:datasets][1],
    )
end


text:
include("../julia/ProgramTunnel_fs_new.jl")

using .ProgramTunnel_fs
using Formatting
using JSON

PTI = ProgramTunnelInfo(reverse_role=true, rotate=30)

sarrs = [ collect(Float64, 1:5),  collect(Float64, 101:105) / 2.0 ]
rarrs = [ collect(Float64, 1:5),  collect(Float64, 1:5)           ]

for t=1:100
    msg = format("[SSM] This is the {:d} time", t)
    println("MSG: ", msg)

    for i = 1:2
        for k = 1:5
            sarrs[i][k] = t*100 + i*10 + k
        end
        rarrs[i] .= 0.0
    end

    sendData(PTI, msg, sarrs)
    msg = recvData!(PTI, rarrs)

    println("Received msg:", msg)
    println("Recved arrs: ")
    print(json(rarrs, 4))
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using OpenImageIO
using Base.Test

import OpenImageIO: channeltype

@testset "TypeDesc" begin
    scalardesc(t) = TypeDesc(t, OpenImageIO.TypeDesc_SCALAR,
                             OpenImageIO.TypeDesc_NOXFORM, 0)
    tupledesc(t, n)  = TypeDesc(t, OpenImageIO.TypeDesc_SCALAR,
                                OpenImageIO.TypeDesc_NOXFORM, n)

    @test TypeDesc(UInt8)   == scalardesc(OpenImageIO.TypeDesc_UINT8)
    @test TypeDesc(Float64) == scalardesc(OpenImageIO.TypeDesc_DOUBLE)
    @test TypeDesc(Void)    == scalardesc(OpenImageIO.TypeDesc_PTR)

    @test TypeDesc(NTuple{2,Int16})   == tupledesc(OpenImageIO.TypeDesc_INT16,2)
    @test TypeDesc(NTuple{3,Float16}) == tupledesc(OpenImageIO.TypeDesc_HALF,3)

    @test channeltype(TypeDesc(OpenImageIO.TypeDesc_UINT32)) == UInt32
    @test channeltype(TypeDesc(OpenImageIO.TypeDesc_FLOAT))  == Float32

    @test string(TypeDesc(Int64)) == "TypeDesc(\"int64\")"
end

@testset "ImageInput" begin
    @test read_image(ImageInput_open("rgbw.png")) == [(0xff,0x00,0x00) (0x00,0xff,0x00);
                                                      (0x00,0x00,0xff) (0xff,0xff,0xff)]
    # No image present
    @test_throws ErrorException ImageInput_open("not_present.png")
    # Present, but not an image
    @test_throws ErrorException ImageInput_open("runtests.jl")
end


text:
@testset "6.2.4 (d+e x)^m cosh(a+b x+c x^2)^n" begin
    (a, b, c, d, e, x, ) = @variables a b c d e x

    #= ::Package:: =#

    #= ::Section:: =#
    #=Integrands*of*the*form*(d+e*x)^m*cosh(a+b*x+c*x^2)^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(d+e*x)^m*cosh(a+b*x+c*x^2)^n=#

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [x^2*cosh(a + b*x + c*x^2), x, 12, (b^2*ℯ^(-a + b^2/(4*c))*sqrt(pi)*erf((b + 2*c*x)/(2*sqrt(c))))/(16*c^(5/2)) + (ℯ^(-a + b^2/(4*c))*sqrt(pi)*erf((b + 2*c*x)/(2*sqrt(c))))/(8*c^(3/2)) + (b^2*ℯ^(a - b^2/(4*c))*sqrt(pi)*erfi((b + 2*c*x)/(2*sqrt(c))))/(16*c^(5/2)) - (ℯ^(a - b^2/(4*c))*sqrt(pi)*erfi((b + 2*c*x)/(2*sqrt(c))))/(8*c^(3/2)) - (b*sinh(a + b*x + c*x^2))/(4*c^2) + (x*sinh(a + b*x + c*x^2))/(2*c)]
    @test_int [x*cosh(a + b*x + c*x^2), x, 6, -((b*ℯ^(-a + b^2/(4*c))*sqrt(pi)*erf((b + 2*c*x)/(2*sqrt(c))))/(8*c^(3/2))) - (b*ℯ^(a - b^2/(4*c))*sqrt(pi)*erfi((b + 2*c*x)/(2*sqrt(c))))/(8*c^(3/2)) + sinh(a + b*x + c*x^2)/(2*c)]
    @test_int [cosh(a + b*x + c*x^2), x, 5, (ℯ^(-a + b^2/(4*c))*sqrt(pi)*erf((b + 2*c*x)/(2*sqrt(c))))/(4*sqrt(c)) + (ℯ^(a - b^2/(4*c))*sqrt(pi)*erfi((b + 2*c*x)/(2*sqrt(c))))/(4*sqrt(c))]
    @test_int [cosh(a + b*x + c*x^2)/x, x, 0, Unintegrable(cosh(a + b*x + c*x^2)/x, x)]
    @test_int [cosh(a + b*x + c*x^2)/x^2 - b*sinh(a + b*x + c*x^2)/x, x, 7, -(cosh(a + b*x + c*x^2)/x) - (1/2)*sqrt(c)*ℯ^(-a + b^2/(4*c))*sqrt(pi)*erf((b + 2*c*x)/(2*sqrt(c))) + (1/2)*sqrt(c)*ℯ^(a - b^2/(4*c))*sqrt(pi)*erfi((b + 2*c*x)/(2*sqrt(c)))]

    @test_int [x^2*cosh(a + b*x - c*x^2), x, 12, -((b^2*ℯ^(a + b^2/(4*c))*sqrt(pi)*erf((b - 2*c*x)/(2*sqrt(c))))/(16*c^(5/2))) - (ℯ^(a + b^2/(4*c))*sqrt(pi)*erf((b - 2*c*x)/(2*sqrt(c))))/(8*c^(3/2)) - (b^2*ℯ^(-a - b^2/(4*c))*sqrt(pi)*erfi((b - 2*c*x)/(2*sqrt(c))))/(16*c^(5/2)) + (ℯ^(-a - b^2/(4*c))*sqrt(pi)*erfi((b - 2*c*x)/(2*sqrt(c))))/(8*c^(3/2)) - (b*sinh(a + b*x - c*x^2))/(4*c^2) - (x*sinh(a + b*x - c*x^2))/(2*c)]
    @test_int [x*cosh(a + b*x - c*x^2), x, 6, -((b*ℯ^(a + b^2/(4*c))*sqrt(pi)*erf((b - 2*c*x)/(2*sqrt(c))))/(8*c^(3/2))) - (b*ℯ^(-a - b^2/(4*c))*sqrt(pi)*erfi((b - 2*c*x)/(2*sqrt(c))))/(8*c^(3/2)) - sinh(a + b*x - c*x^2)/(2*c)]
    @test_int [cosh(a + b*x - c*x^2), x, 5, -((ℯ^(a + b^2/(4*c))*sqrt(pi)*erf((b - 2*c*x)/(2*sqrt(c))))/(4*sqrt(c))) - (ℯ^(-a - b^2/(4*c))*sqrt(pi)*erfi((b - 2*c*x)/(2*sqrt(c))))/(4*sqrt(c))]
    @test_int [cosh(a + b*x - c*x^2)/x, x, 0, Unintegrable(cosh(a + b*x - c*x^2)/x, x)]
    @test_int [cosh(a + b*x - c*x^2)/x^2 - b*sinh(a + b*x - c*x^2)/x, x, 7, -(cosh(a + b*x - c*x^2)/x) + (1/2)*sqrt(c)*ℯ^(a + b^2/(4*c))*sqrt(pi)*erf((b - 2*c*x)/(2*sqrt(c))) - (1/2)*sqrt(c)*ℯ^(-a - b^2/(4*c))*sqrt(pi)*erfi((b - 2*c*x)/(2*sqrt(c)))]

    @test_int [x^2*cosh(1/4 + x + x^2), x, 12, (-(3/16))*sqrt(pi)*erf((1/2)*(-1 - 2*x)) - (1/16)*sqrt(pi)*erfi((1/2)*(1 + 2*x)) - (1/4)*sinh(1/4 + x + x^2) + (1/2)*x*sinh(1/4 + x + x^2)]
    @test_int [x*cosh(1/4 + x + x^2), x, 6, (1/8)*sqrt(pi)*erf((1/2)*(-1 - 2*x)) - (1/8)*sqrt(pi)*erfi((1/2)*(1 + 2*x)) + (1/2)*sinh(1/4 + x + x^2)]
    @test_int [cosh(1/4 + x + x^2), x, 5, (-(1/4))*sqrt(pi)*erf((1/2)*(-1 - 2*x)) + (1/4)*sqrt(pi)*erfi((1/2)*(1 + 2*x))]
    @test_int [cosh(1/4 + x + x^2)/x, x, 0, Unintegrable(cosh(1/4 + x + x^2)/x, x)]
    @test_int [cosh(1/4 + x + x^2)/x^2, x, 6, -(cosh(1/4 + x + x^2)/x) + (1/2)*sqrt(pi)*erf((1/2)*(-1 - 2*x)) + (1/2)*sqrt(pi)*erfi((1/2)*(1 + 2*x)) + Unintegrable(sinh(1/4 + x + x^2)/x, x)]

    @test_int [x^2*cosh(a + b*x + c*x^2)^2, x, 14, x^3/6 + (b^2*ℯ^(-2*a + b^2/(2*c))*sqrt(pi/2)*erf((b + 2*c*x)/(sqrt(2)*sqrt(c))))/(32*c^(5/2)) + (ℯ^(-2*a + b^2/(2*c))*sqrt(pi/2)*erf((b + 2*c*x)/(sqrt(2)*sqrt(c))))/(32*c^(3/2)) + (b^2*ℯ^(2*a - b^2/(2*c))*sqrt(pi/2)*erfi((b + 2*c*x)/(sqrt(2)*sqrt(c))))/(32*c^(5/2)) - (ℯ^(2*a - b^2/(2*c))*sqrt(pi/2)*erfi((b + 2*c*x)/(sqrt(2)*sqrt(c))))/(32*c^(3/2)) - (b*sinh(2*a + 2*b*x + 2*c*x^2))/(16*c^2) + (x*sinh(2*a + 2*b*x + 2*c*x^2))/(8*c)]
    @test_int [x*cosh(a + b*x + c*x^2)^2, x, 8, x^2/4 - (b*ℯ^(-2*a + b^2/(2*c))*sqrt(pi/2)*erf((b + 2*c*x)/(sqrt(2)*sqrt(c))))/(16*c^(3/2)) - (b*ℯ^(2*a - b^2/(2*c))*sqrt(pi/2)*erfi((b + 2*c*x)/(sqrt(2)*sqrt(c))))/(16*c^(3/2)) + sinh(2*a + 2*b*x + 2*c*x^2)/(8*c)]
    @test_int [cosh(a + b*x + c*x^2)^2, x, 7, x/2 + (ℯ^(-2*a + b^2/(2*c))*sqrt(pi/2)*erf((b + 2*c*x)/(sqrt(2)*sqrt(c))))/(8*sqrt(c)) + (ℯ^(2*a - b^2/(2*c))*sqrt(pi/2)*erfi((b + 2*c*x)/(sqrt(2)*sqrt(c))))/(8*sqrt(c))]
    @test_int [cosh(a + b*x + c*x^2)^2/x, x, 2, (1/2)*Unintegrable(cosh(2*a + 2*b*x + 2*c*x^2)/x, x) + log(x)/2]

    @test_int [x^2*cosh(a + b*x - c*x^2)^2, x, 14, x^3/6 - (b^2*ℯ^(2*a + b^2/(2*c))*sqrt(pi/2)*erf((b - 2*c*x)/(sqrt(2)*sqrt(c))))/(32*c^(5/2)) - (ℯ^(2*a + b^2/(2*c))*sqrt(pi/2)*erf((b - 2*c*x)/(sqrt(2)*sqrt(c))))/(32*c^(3/2)) - (b^2*ℯ^(-2*a - b^2/(2*c))*sqrt(pi/2)*erfi((b - 2*c*x)/(sqrt(2)*sqrt(c))))/(32*c^(5/2)) + (ℯ^(-2*a - b^2/(2*c))*sqrt(pi/2)*erfi((b - 2*c*x)/(sqrt(2)*sqrt(c))))/(32*c^(3/2)) - (b*sinh(2*a + 2*b*x - 2*c*x^2))/(16*c^2) - (x*sinh(2*a + 2*b*x - 2*c*x^2))/(8*c)]
    @test_int [x*cosh(a + b*x - c*x^2)^2, x, 8, x^2/4 - (b*ℯ^(2*a + b^2/(2*c))*sqrt(pi/2)*erf((b - 2*c*x)/(sqrt(2)*sqrt(c))))/(16*c^(3/2)) - (b*ℯ^(-2*a - b^2/(2*c))*sqrt(pi/2)*erfi((b - 2*c*x)/(sqrt(2)*sqrt(c))))/(16*c^(3/2)) - sinh(2*a + 2*b*x - 2*c*x^2)/(8*c)]
    @test_int [cosh(a + b*x - c*x^2)^2, x, 7, x/2 - (ℯ^(2*a + b^2/(2*c))*sqrt(pi/2)*erf((b - 2*c*x)/(sqrt(2)*sqrt(c))))/(8*sqrt(c)) - (ℯ^(-2*a - b^2/(2*c))*sqrt(pi/2)*erfi((b - 2*c*x)/(sqrt(2)*sqrt(c))))/(8*sqrt(c))]
    @test_int [cosh(a + b*x - c*x^2)^2/x, x, 2, (1/2)*Unintegrable(cosh(2*a + 2*b*x - 2*c*x^2)/x, x) + log(x)/2]

    @test_int [x^2*cosh(1/4 + x + x^2)^2, x, 14, x^3/6 + (1/16)*sqrt(pi/2)*erf((1 + 2*x)/sqrt(2)) - (1/16)*sinh(1/2 + 2*x + 2*x^2) + (1/8)*x*sinh(1/2 + 2*x + 2*x^2)]
    @test_int [x*cosh(1/4 + x + x^2)^2, x, 8, x^2/4 - (1/16)*sqrt(pi/2)*erf((1 + 2*x)/sqrt(2)) - (1/16)*sqrt(pi/2)*erfi((1 + 2*x)/sqrt(2)) + (1/8)*sinh(1/2 + 2*x + 2*x^2)]
    @test_int [cosh(1/4 + x + x^2)^2, x, 7, x/2 + (1/8)*sqrt(pi/2)*erf((1 + 2*x)/sqrt(2)) + (1/8)*sqrt(pi/2)*erfi((1 + 2*x)/sqrt(2))]
    @test_int [cosh(1/4 + x + x^2)^2/x, x, 2, (1/2)*Unintegrable(cosh(1/2 + 2*x + 2*x^2)/x, x) + log(x)/2]

    @test_int [(d + e*x)^2*cosh(a + b*x + c*x^2), x, 12, (e^2*ℯ^(-a + b^2/(4*c))*sqrt(pi)*erf((b + 2*c*x)/(2*sqrt(c))))/(8*c^(3/2)) + ((2*c*d - b*e)^2*ℯ^(-a + b^2/(4*c))*sqrt(pi)*erf((b + 2*c*x)/(2*sqrt(c))))/(16*c^(5/2)) - (e^2*ℯ^(a - b^2/(4*c))*sqrt(pi)*erfi((b + 2*c*x)/(2*sqrt(c))))/(8*c^(3/2)) + ((2*c*d - b*e)^2*ℯ^(a - b^2/(4*c))*sqrt(pi)*erfi((b + 2*c*x)/(2*sqrt(c))))/(16*c^(5/2)) + (e*(2*c*d - b*e)*sinh(a + b*x + c*x^2))/(4*c^2) + (e*(d + e*x)*sinh(a + b*x + c*x^2))/(2*c)]
    @test_int [(d + e*x)*cosh(a + b*x + c*x^2), x, 6, ((2*c*d - b*e)*ℯ^(-a + b^2/(4*c))*sqrt(pi)*erf((b + 2*c*x)/(2*sqrt(c))))/(8*c^(3/2)) + ((2*c*d - b*e)*ℯ^(a - b^2/(4*c))*sqrt(pi)*erfi((b + 2*c*x)/(2*sqrt(c))))/(8*c^(3/2)) + (e*sinh(a + b*x + c*x^2))/(2*c)]
    @test_int [cosh(a + b*x + c*x^2)/(d + e*x), x, 0, Unintegrable(cosh(a + b*x + c*x^2)/(d + e*x), x)]

    @test_int [(d + e*x)^2*cosh(a + b*x + c*x^2)^2, x, 14, (d + e*x)^3/(6*e) + (e^2*ℯ^(-2*a + b^2/(2*c))*sqrt(pi/2)*erf((b + 2*c*x)/(sqrt(2)*sqrt(c))))/(32*c^(3/2)) + ((2*c*d - b*e)^2*ℯ^(-2*a + b^2/(2*c))*sqrt(pi/2)*erf((b + 2*c*x)/(sqrt(2)*sqrt(c))))/(32*c^(5/2)) - (e^2*ℯ^(2*a - b^2/(2*c))*sqrt(pi/2)*erfi((b + 2*c*x)/(sqrt(2)*sqrt(c))))/(32*c^(3/2)) + ((2*c*d - b*e)^2*ℯ^(2*a - b^2/(2*c))*sqrt(pi/2)*erfi((b + 2*c*x)/(sqrt(2)*sqrt(c))))/(32*c^(5/2)) + (e*(2*c*d - b*e)*sinh(2*a + 2*b*x + 2*c*x^2))/(16*c^2) + (e*(d + e*x)*sinh(2*a + 2*b*x + 2*c*x^2))/(8*c)]
    @test_int [(d + e*x)*cosh(a + b*x + c*x^2)^2, x, 8, (d + e*x)^2/(4*e) + ((2*c*d - b*e)*ℯ^(-2*a + b^2/(2*c))*sqrt(pi/2)*erf((b + 2*c*x)/(sqrt(2)*sqrt(c))))/(16*c^(3/2)) + ((2*c*d - b*e)*ℯ^(2*a - b^2/(2*c))*sqrt(pi/2)*erfi((b + 2*c*x)/(sqrt(2)*sqrt(c))))/(16*c^(3/2)) + (e*sinh(2*a + 2*b*x + 2*c*x^2))/(8*c)]
    @test_int [cosh(a + b*x + c*x^2)^2/(d + e*x), x, 2, (1/2)*Unintegrable(cosh(2*a + 2*b*x + 2*c*x^2)/(d + e*x), x) + log(d + e*x)/(2*e)]

    #= ::Subsubsection:: =#
    #=n<0=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*(d+e*x)^m*cosh(a+b*x+c*x^2)^(n/2)=#
    end

max_stars_repo_path: test/6 Hyperbolic functions/6.2 Hyperbolic cosine/6.2.4 (d+e x)^m cosh(a+b x+c x^2)^n.jl


text:
@mlj_model mutable struct APIx0 <: Supervised
    f0::Int
end
@mlj_model mutable struct APIx0b <: Supervised
    f0::Int
end

mutable struct APIx1 <: Static end

@testset "fit-x" begin
    m0 = APIx0(f0=1)
    m1 = APIx0b(f0=3)
    # no weight support: fallback
    M.fit(m::APIx0, v::Int, X, y) = (5, nothing, nothing)
    @test fit(m0, 1, randn(2), randn(2), 5) == (5, nothing, nothing)
    # with weight support: use
    M.fit(m::APIx0b, v::Int, X, y, w) = (7, nothing, nothing)
    @test fit(m1, 1, randn(2), randn(2), 5) == (7, nothing, nothing)
    # default fitted params
    @test M.fitted_params(m1, 7) == (fitresult=7,)
    # static
    s1 = APIx1()
    @test fit(s1, 1, 0) == (nothing, nothing, nothing)

    #update fallback = fit
    @test update(m0, 1, 5, nothing, randn(2), 5) == (5, nothing, nothing)
end


text:
function Base.:+(x::Variable, y::Variable)
    return LinearVariables([x.idx,y.idx],[1,1])
end

function Base.:-(x::Variable, y::Variable)
    return LinearVariables([x.idx,y.idx],[1,-1])
end

function Base.:+(x::LinearVariables, y::Variable)
    lv = LinearVariables(x.indices, x.coeffs)
    push!(lv.indices, y.idx)
    push!(lv.coeffs, 1)
    return lv
end

function Base.:-(x::LinearVariables, y::Variable)
    lv = LinearVariables(x.indices, x.coeffs)
    push!(lv.indices, y.idx)
    push!(lv.coeffs, -1)
    return lv
end

function Base.:+(x::Variable, y::Int)
    return LinearVariables([x.idx], [1])+y
end

function Base.:+(x::LinearVariables, y::Int)
    lv = LinearVariables(x.indices, x.coeffs)
    push!(lv.indices, 0)
    push!(lv.coeffs, y)
    return lv
end

function Base.:-(x::Variable, y::Int)
    return LinearVariables([x.idx], [1])-y
end

function Base.:-(x::LinearVariables, y::Int)
    return x+(-y)
end

function Base.:+(x::Variable, y::LinearVariables)
    return y+x # commutative
end

function Base.:+(x::LinearVariables, y::LinearVariables)
    return LinearVariables(vcat(x.indices, y.indices), vcat(x.coeffs, y.coeffs))
end

function Base.:-(x::LinearVariables, y::LinearVariables)
    return LinearVariables(vcat(x.indices, y.indices), vcat(x.coeffs, -y.coeffs))
end

function Base.:*(x::Int, y::Variable)
    return LinearVariables([y.idx],[x])
end

function Base.:*(y::Variable, x::Int)
    return LinearVariables([y.idx],[x])
end

function simplify(x::LinearVariables)
    set_indices = Set(x.indices)
    # if unique
    if length(set_indices) == length(x.indices)
        if 0 in set_indices
            indices = Int[]
            coeffs = Int[]
            lhs = 0
            for i=1:length(x.indices)
                if x.indices[i] == 0
                    lhs = x.coeffs[i]
                else
                    push!(indices, x.indices[i])
                    push!(coeffs, x.coeffs[i])
                end
            end
            return indices, coeffs, lhs
        else
            return x.indices, x.coeffs, 0
        end
    end

    perm = sortperm(x.indices)
    x.indices = x.indices[perm]
    x.coeffs = x.coeffs[perm]
    indices = Int[x.indices[1]]
    coeffs = Int[x.coeffs[1]]

    last_idx = x.indices[1]

    for i=2:length(x.indices)
        if x.indices[i] == last_idx
            coeffs[end] += x.coeffs[i]
        else
            last_idx = x.indices[i]
            push!(indices, x.indices[i])
            push!(coeffs, x.coeffs[i])
        end
    end

    lhs = 0
    if indices[1] == 0
        lhs = coeffs[1]
        indices = indices[2:end]
        coeffs = coeffs[2:end]
    end

    return indices, coeffs, lhs
end


text:
module ImageQualityIndexes

using OffsetArrays
using ImageCore
# Where possible we avoid a direct dependency to reduce the number of [compat] bounds
using ImageCore.MappedArrays
using ImageCore: NumberLike, Pixel, GenericImage, GenericGrayImage
using ImageDistances, ImageFiltering
using Statistics: mean, std
using Base.Iterators: repeated, flatten

include("generic.jl")
include("psnr.jl")
include("ssim.jl")
include("msssim.jl")
include("colorfulness.jl")
include("deprecations.jl")

export
    # generic
    assess,

    # Peak signal-to-noise ratio
    PSNR,
    assess_psnr,

    # Structral Similarity
    SSIM,
    assess_ssim,

    # Multi Scale Structural Similarity
    MSSSIM,
    assess_msssim,

    # Colorfulness
    HASLER_AND_SUSSTRUNK_M3,
    hasler_and_susstrunk_m3,
    colorfulness

end # module


text:
"""
    metareader(md::MatrixDescriptor, key::AbstractString)

Return specific data files (matrix, rhs, solution, or other metadata).
The `key` must be contained in data.metadata or an `DataError` is thrown.
Hint: `md.A, md.b, md.x` also deliver the matrix, rhs, and solution.
So this is needed specifically for other metadata.
"""
function metareader(mdesc::MatrixDescriptor{<:RemoteMatrixData}, name::AbstractString)
    name = metastring_reverse(mdesc.data, name)
    get!(mdesc.cache, name) do
        _metareader(mdesc.data, name)
    end
end

function metareader(mdesc::MatrixDescriptor{<:GeneratedMatrixData}, name::AbstractString)
    fillcache!(mdesc)
    dao = mdesc.cache[]
    if name == "A" && dao isa AbstractArray
        dao
    else
        try
            getproperty(dao, Symbol(name))
        catch
            daterr("this instance of '$(typeof(dao))' has no property '$name'")
        end
    end
end

function _metareader(data::RemoteMatrixData, name::AbstractString)
    if name in data.metadata
        path = joinpath(dirname(matrixfile(data)), name)
        endswith(name, ".mtx") ? mmread(path) : read(path, String)
    else
        daterr("$(data.name) has no metadata $name")
    end
end

function fillcache!(mdesc::MatrixDescriptor{<:GeneratedMatrixData})
    dao = mdesc.cache[]
    if dao === nothing
        dao = mdesc.cache[] = mdesc.data.func(mdesc.args...)
        dao !== nothing || daterr("function $(mdesc.data.func) returned `nothing`")
    end
    dao
end

# This a the preferred API to access metadata.
import Base: getproperty, propertynames, getindex

function getproperty(mdesc::MatrixDescriptor{T}, s::Symbol) where T
    s in (:data, :args, :cache) && return getfield(mdesc, s)
    s in metasymbols(mdesc) && return metareader(mdesc, string(s))
    if T <: RemoteMatrixData
        s in (:m, :n, :nnz, :dnz) && return getfield(mdesc.data.header, s)
    else
        s == :m && return size(mdesc.A, 1)
        s == :n && return size(mdesc.A, 2)
        s == :nnz && return count(mdesc.A .!= 0)
        s == :dnz && return 0
    end
    metareader(mdesc, string(s))
end

function propertynames(mdesc::MatrixDescriptor{T}; private=false) where T
    props = Symbol[]
    append!(props, metasymbols(mdesc))
    append!(props, [:m, :n, :nnz, :dnz])
    private && append!(props, fieldnames(MatrixDescriptor))
    props
end

function getproperty(data::RemoteMatrixData, s::Symbol)
    s in (:name, :id, :header, :properties, :metadata) && return getfield(data, s)
    s in fieldnames(MetaInfo) && return getfield(data.header, s)
    getfield(data, s)
end

function propertynames(data::RemoteMatrixData; private=false)
    props = [PROPS...]
    private && append!(props, fieldnames(RemoteMatrixData))
    props
end

function propertynames(data::GeneratedMatrixData; private=false)
    private ? fieldnames(GeneratedMatrixData) : [:name, :id]
end

"""
    metasymbols(md::MatrixDescriptor)

Return a vector of symbols, which point to metadata of the problem.
These symbols `s` may be used to access the objects by `getproperty(md, s)`
or by `getindex(md, s)`. The syntax `md.S` is preferred, if `S` is a constant
`Julia` word. In any case `md[s]` is possible.

Example:
`md = mdopen("*/bfly"); A = md.A; co = md.coord; txt10 = md["Gname_10.txt"]`
"""
metasymbols(md::MatrixDescriptor{<:RemoteMatrixData}) = metasymbols(md.data)
function metasymbols(data::RemoteMatrixData)
    Symbol.(metastring.(data.name, metadata(data)))
end
function metasymbols(md::MatrixDescriptor{<:GeneratedMatrixData})
    mdc = md.cache[]
    mdc isa Array || mdc === nothing ? [:A] : propertynames(mdc)
end
metasymbols(data::MatrixData) = [:A]

function Base.getindex(mdesc::MatrixDescriptor, name::Union{Symbol,AbstractString})
    metareader(mdesc, string(name))
end

function (mdesc::MatrixDescriptor)(name::Union{Symbol,AbstractString})
    metareader(mdesc, string(name))
end

#internal helper to select special metadata (matrix, rhs, or solution)
function metaname(data::RemoteMatrixData, exli::AbstractString...)
    base = rsplit(data.name, '/', limit=2)[end]
    meda = metadata(data)
    for ext in exli
        f = metastring_reverse(data, ext)
        if f in meda
            return f
        end
    end
    daterr("unknown metadata extensions: $exli - available $(String.(data.metadata))")
end

"""
    metastring(name, metaname)

In the standard cases convert full meta name to abbreviated form.

+ given `name = "abc"`, then
+ `"abc.mtx" => "A"`
+ `"abc_def.mtx" => "def"`
+ `"abc_def.txt" => "def.txt"`
+ `"abc-def.mtx" => "abc-def.mtx"`
"""
function metastring(name::AbstractString, metaname::AbstractString)
    MTX = ".mtx"
    base = split(name, '/')[end]
    meta = metaname
    if meta == string(base, MTX)
        meta = "A"
    elseif startswith(meta, string(base, '_'))
        meta = meta[sizeof(base)+2:end]
    end
    es = rsplit(meta, '.', limit=2)
    if endswith(meta, MTX) && meta != metaname
        meta = meta[1:end-sizeof(MTX)]
    end
    meta
end

"""
    metastring_reverse(data::MatrixData, abbreviation::String)

Given a `MatrixData` object and an abbreviation, return a matching full metadata name
for the abbreviation. If no full name is found, return the abbreviation unchanged.
"""
function metastring_reverse(data::RemoteMatrixData, metaabbr::AbstractString)
    MTX = ".mtx"
    base = split(data.name, '/')[end]
    metaabbr == "A" && return string(base, MTX)
    mdata = metadata(data)
    meta = string(base, '_', metaabbr)
    meta in mdata && return meta
    meta = string(meta, MTX)
    meta in mdata && return meta
    metaabbr
end
metastring_reverse(::MatrixData, metaabbr::AbstractString) = metaabbr


text:
#=
caller_name.jl
Use Stacktrace to show information about the
file, line number and calling function

Useful for messages and debugging.

Jeff Fessler, 2019-06-27, University of Michigan
=#

export caller_name

using Test: @test

"""
`caller_name() or caller_name(;level=4)`

Return "filename line fun():" as a string
to describe where this function was called.

Stack levels:
- 1: #caller_name
- 2: caller_name()
- 3: function that invoked caller()
- 4: the calling function we want to return
Hence the default level is 4,
but we increment it by one in case user says `@show caller_name()`
in which case stack[3] is a macro expansion.
"""
function caller_name(; level::Integer=4)
	stack = stacktrace()
	level += (stack[3].func == Symbol("macro expansion")) # trick
	(level < 1 || level > length(stack)) && throw("bad level $level")
	frame = stack[level]
	file = string(frame.file)
	file = basename(file)
	line = frame.line
	func = frame.func
#	display("$which, $file $line $func(): ")

	return "$file $line $func(): "
end

"""
`caller_name(:test)`
"""
function caller_name(test::Symbol)
	test != :test && throw(ArgumentError("test $test"))

	function f2()
		caller_name()
	end

	line = 2 + @__LINE__
	function f1()
		f2() # this is two lines below @__LINE__ above
	end

	@test isa(f1(), String)
	@test f1()[end-12:end] == ".jl $line f1(): "

	true
end

# caller_name(:test)


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# Maximum Expected Coverage Location Problem (MEXCLP)
# From paper: "A maximum expected covering location model: formulation, properties and heuristic solution"

"""
	solveMexclp!(sim::Simulation;
		numAmbs::Int = sim.numAmbs,
		busyFraction::Float = 0.5,
		demandWeights::Dict{Priority,Float} = Dict([p => 1.0 for p in priorities]),
		stationCapacities::Vector{Int} = [station.capacity for station in sim.stations],
		results::Dict = Dict())
Solves the Maximum Expected Coverage Location Problem (MEXCLP) for `sim` and returns the number of ambulances to assign to each station, also the converse is returned - a station index for each ambulance.
The problem assumes that all ambulances are equivalent.

# Keyword arguments
- `numAmbs` is the number of ambulances to solve for, must be >= 1
- `busyFraction` is the fraction of time that ambulances are busy; should be within [0,1] though this is not enforced
- `demandWeights` is the weight to apply to each demand priority for the objective function
- `stationCapacities` is the maximum number of ambulances that each station can hold
- `results` will store results of mexclp such as the objective value and decision variable values
"""
function solveMexclp!(sim::Simulation;
	numAmbs::Int = sim.numAmbs,
	busyFraction::Float = 0.5,
	demandWeights::Dict{Priority,Float} = Dict([p => 1.0 for p in priorities]),
	stationCapacities::Vector{Int} = [station.capacity for station in sim.stations],
	results::Dict = Dict())

	@assert(numAmbs >= 1, "need at least 1 ambulance for mexclp")
	@assert(sim.travel.numSets == 1) # otherwise need to solve mexclp for each travel set?

	@assert(length(stationCapacities) == sim.numStations)
	@assert(all(stationCapacities .>= 0), "station capacities must be non-negative")
	stationCapacities = min.(stationCapacities, numAmbs) # reduce values where capacity > numAmbs
	@assert(sum(stationCapacities) >= numAmbs, "the total capacity for ambulances at stations is less than the number of ambulances")

	# initialise demand and demand coverage data if not already initialised
	sim.demand.initialised || initDemand!(sim)
	@assert(sim.demand.numSets == 1) # otherwise need to solve mexclp for each demand set?
	sim.demandCoverage.initialised || initDemandCoverage!(sim)

	# shorthand
	numStations = sim.numStations
	currentTime = sim.startTime

	# get demand point coverage data
	pointData = Dict{Vector{Int},Float}() # pointData[pointStations] = pointDemand
	for demandPriority in priorities
		if !haskey(demandWeights, demandPriority) || demandWeights[demandPriority] == 0
			continue
		end

		pointsCoverageMode = getPointsCoverageMode!(sim, demandPriority, currentTime)
		demandMode = getDemandMode!(sim.demand, demandPriority, currentTime)
		pointSetsDemands = getPointSetsDemands!(sim, demandPriority, currentTime; pointsCoverageMode = pointsCoverageMode) * demandMode.rasterMultiplier * demandWeights[demandPriority]

		for (i,stationSet) in enumerate(pointsCoverageMode.stationSets)
			get!(pointData, stationSet, 0.0)
			pointData[stationSet] += pointSetsDemands[i]
		end
	end
	pointStations = collect(keys(pointData))
	pointDemands = collect(values(pointData))
	# point j has demand pointDemands[j] and is covered by stations pointStations[j]
	numPoints = length(pointDemands) # shorthand

	# calculate benefit of covering each point with a kth ambulance
	marginalBenefit = (busyFraction.^[0:numAmbs-1;])*(1-busyFraction) # point cover benefit values, for single demand
	pointCoverCountBenefit = [pointDemands[j] * marginalBenefit for j = 1:numPoints] # pointCoverCountBenefit[j][k] is the benefit of adding a kth ambulance to cover point j

	#######
	# IP

	# shorthand:
	a = numAmbs
	s = numStations
	p = numPoints

	model = Model()

	if !all(v -> issorted(v, rev=true), pointCoverCountBenefit)
		set_optimizer(model, with_optimizer(GLPK.Optimizer)) # solve speed not tested

		@variables(model, begin
			(x[i=1:s] >= 0, Int) # x[i] = number of ambulances assigned to station i
			(y[j=1:p,k=1:a], Bin) # y[j,k] = true if point j is covered by at least k ambulances
		end)

		@constraints(model, begin
			(useAllAmbs, sum(x) == a)
			(pointCoverCount[j=1:p], sum(y[j,:]) == sum(x[pointStations[j]]))
			(pointCoverOrder[j=1:p, k=1:(a-1)], y[j,k] >= y[j,k+1])
			(stationMaxNumAmbs[i=1:s], x[i] <= stationCapacities[i])
		end)
	else
		# pointCoverCountBenefit[j] is non-increasing (because busyFraction >= 0), so:
		# - can have y variables as non-binary (but still constrained to be between 0 and 1)
		# - can leave out the constraint 'pointCoverOrder'
		# - could have x variables as non-integer, though if there are multiple solutions then x values might not be naturally integer, so will leave the integer constraint

		set_optimizer(model, with_optimizer(Cbc.Optimizer, logLevel=0)) # solve speed not tested

		@variables(model, begin
			(x[i=1:s] >= 0, Int) # x[i] = number of ambulances assigned to station i
			(0 <= y[j=1:p,k=1:a] <= 1) # y[j,k] = true if point j is covered by at least k ambulances
		end)

		@constraints(model, begin
			(useAllAmbs, sum(x) == a)
			(pointCoverCount[j=1:p], sum(y[j,:]) == sum(x[pointStations[j]]))
			(stationMaxNumAmbs[i=1:s], x[i] <= stationCapacities[i])
		end)
	end

	@expressions(model, begin
		# (expectedPointCoverage[j=1:p], sum(y[j,k] * pointCoverCountBenefit[j][k] for k=1:a))
		(expectedCoverage, sum(y[j,k] * pointCoverCountBenefit[j][k] for j=1:p, k=1:a))
	end)

	# solve
	@objective(model, Max, expectedCoverage)
	optimize!(model)
	@assert(termination_status(model) == MOI.OPTIMAL)
	xValue = JuMP.value.(x); yValue = JuMP.value.(y) # JuMP and LightXML both export value()
	objVal = JuMP.value(expectedCoverage)

	results[:x] = xValue
	results[:y] = yValue
	results[:objVal] = objVal

	# extract solution
	stationsNumAmbs = convert(Vector{Int}, round.(xValue)) # solution; stationsNumAmbs[i] gives number of ambulances to allocate to station i
	pointSlotCover = convert(Array{Bool,2}, round.(yValue)) # pointSlotCover[j,k] = true if point j is covered by at least k ambulances

	if checkMode
		# check constraints
		@assert(sum(stationsNumAmbs) == a) # useAllAmbs constraint
		@assert(all(j -> sum(pointSlotCover[j,:]) == sum(stationsNumAmbs[pointStations[j]]), 1:p)) # pointCoverCount constraint
		@assert(all(i -> 0 <= stationsNumAmbs[i] <= stationCapacities[i], 1:s)) # stationMaxNumAmbs constraint
		@assert(all(j -> issorted(pointSlotCover[j], rev=true), 1:p)) # pointCoverOrder constraint
	end

	# convert to a deployment
	deployment = stationsNumAmbsToDeployment(stationsNumAmbs)

	return stationsNumAmbs, deployment # also return model?
end


text:
# This file contains main function to generate network ans solve differential equations
# and to run the series of repetitions of given model.
# Apart from that, functions also make some initial and final states analysis and save the results to Matlab file.
# Code is optimized

# module MeasureHeiderBalance

using DrWatson
using DifferentialEquations
using Plots
using LinearAlgebra
using Dates

# function to measure HB in the case of a single-layered system with attributes
# n - number of agents
# attr - attribute type; contains all information about attributes
# gamma - coupling strength
# maxtime - maximal time of calculating differential equations
# ode_fun - function calculating derivatives
# solver - function solving different equations
# show_plot - parameter whether the solution will be plotted. If not transition values are not saved.
# input is optional extra parameters with initial conditions of RL and AL (u0, xy_attr).
# If input is empty, random initial conditions are generated.
# If not, the first value (possibly empty) should contain initial conditions for RL.
# The 2nd value of input may contain initial conditions for AL.
#
# Returns the tuple of:
# ishb_sim_par - array of Boolean values [is it balanced state, is AL in balanced state,
# similarity between RL and AL, balanced ratio of RL, is it paradise state,
# is it hell state, counts od different triads (Deltas: Delta0, Delta1, Delta2, Delta3),
# is it weakly balanced state, local polarization measure, is the system globally polarized],
# time of simulation to reach stable solution (or maxtime if not reached),
# values of RL weights at the end of simulation,
# initial conditions of RL,
# initial conditions of AL,
# whole sol (if transition values are not saved, see show_plot, then they are not here).
function calc_heider_attr(
    n::Int,
    attr::AbstractAttributes,
    gamma::Float64,
    maxtime::Float64,
    ode_fun::Function,
    solver,
    show_plot::Bool,
    input...,
)
    #initial conditions
    need_init_u0 = true
    need_init_xy = true
    if length(input) > 0
        if !isempty(input[1])
            need_init_u0 = false
            u0 = input[1]
        end
        if length(input) == 2
            if !isempty(input[2])
                need_init_xy = false
                xy_attr = input[2]
            end
        end
    end

    if need_init_u0
        u0 = triu((rand(n, n) * 2) .- 1, 1)
    end
    if need_init_xy
        val0_attr = get_attributes(attr, n)
        xy_attr = get_attribute_layer_weights(attr, val0_attr)
    end

    # help variable
    mask = triu(trues(size(u0)), 1)

    condition_here(u, t, integrator) = condition2(u, t, integrator, mask)
    affect!(integrator) = terminate!(integrator)
    cb = DiscreteCallback(condition_here, affect!)

    #ode and model parameters
    tspan = (0.0, maxtime)

    lay1mul = zeros(n, n)
    x_sim = zeros(n, n)

    p = (n, gamma .* xy_attr, lay1mul, x_sim, mask)

    #solving
    prob = ODEProblem(ode_fun, u0, tspan, p)
    sol = solve(
        prob,
        solver,
        reltol = 1e-6,
        abstol = 1e-12,
        callback = cb,
        isoutofdomain = (u, p, t) -> any(x -> abs.(x) >= 1, u),
        save_everystep = show_plot,
    )

    #estimating output
    paradise = is_paradise(sol.u[end], n)
    Deltas = get_triad_counts(sol.u[end], n)
    weak_balance_in_complete_graph = Deltas[1+1] == 0
    local_polarization = get_local_polarization(Deltas)
    global_polarization = weak_balance_in_complete_graph && (!paradise)
    ishb_sim_par = [
        is_hb(sol.u[end], n),
        is_hb(xy_attr, n),
        get_similarity(sol.u[end], xy_attr, n),
        get_balanced_ratio(sol.u[end], n),
        paradise,
        is_hell(sol.u[end], n),
        Deltas,
        weak_balance_in_complete_graph,
        local_polarization,
        global_polarization,
    ]

    if show_plot

        h = plot(
            sol,
            linewidth = 5,
            title = "Solution to the linear ODE with a thick line",
            xaxis = "Time (t)",
            yaxis = "weights(t)",
            ylim = [-1, +1],
            vars = reshape(1:n^2, n, n)[mask],
        )
        legend = false
        display(h)
        #gui()
    end

    return (ishb_sim_par, sol.t[end], sol.u[end], u0, xy_attr, sol)
end
export calc_heider_attr

# Function creating results file. It parses parameters creating a unique file name and saves initially this file.
# For parameters description see `using_heider_attr`.
# Returns a Tuple with:
#   object of `Result` DataType,
#   `filename`
#   filename without extension
function initialize_file(
    n::Int,
    attr::AbstractAttributes,
    gammas::Vector{Float64},
    maxtime::Float64,
    ode_fun_name::String,
    files_folder::Vector{String},
    filename_prefix::String,
)

    r = Result(n, attr, gammas, maxtime, ode_fun_name)

    prefix = filename_prefix * Dates.format(now(), "yyyy-mm-ddTHH:MM:SS")
    file_params = savename(prefix, r, "mat", sort = false)
    file_params = replace(file_params, "attr_degeneracy" => "v")

    filename = projectdir(files_folder..., file_params)
    save_result(r, filename) #''allocating'' place
    return r, filename
end

# Function using `calc_heider_attr` function to simulate a number of repetitions
# of solving the system having given parameters and random initial conditions:
# Thus, this simulates the influence of attributes on preventing from forming
# a polarized state.
#
# Parameters:
# n - number of agents
# attr - attribute type; contains all information about attributes
# gammas - a vector of coupling strengths to simulate
# zmax - number of repetitions
# maxtime - maximal time of calculating differential equations
# ode_fun_name - function name calculating derivatives (`string`)
#
# Additional named arguments are:
# disp_each - display status every how many ratio of realizations (if 0 then do not display) (default 0.5)
# disp_more_every - display status every certain number of seconds (if 0, then do not display) (default 600)
# save_each - how often (in seconds) should results be saved (if 0, then save at the end) (default 600)
# files_folder - folder name inside the project the results should be saved (default "data").
#       This can be also an array of folder names in the correct folder structure.
# filename_prefix - string that should start each simulation results' file (default "")
#
# `disp_more_every` and `save_each` work like that, that if the specified time has past
# then sth is displayed/saved. But it doesn't mean exact time of action.
# When the specified time has past means when the current function
# (heider function with parameters) finishes.
#
# Returns `Result` object.
function using_heider_attr(
    n::Int,
    attr::AbstractAttributes,
    gammas::Vector{Float64},
    zmax::Int,
    maxtime::Float64,
    ode_fun_name::String;
    disp_each = 0.5,
    disp_more_every = 600,
    save_each = 600,
    files_folder::Vector{String} = ["data"],
    filename_prefix::String = "",
)

    ode_fun = getfield(PolarizationFramework, Symbol(ode_fun_name))
    solver = AutoTsit5(Rodas5(autodiff = false))

    r, filename = initialize_file(
        n,
        attr,
        gammas,
        maxtime,
        ode_fun_name,
        files_folder,
        filename_prefix,
    )

    # time preparation
    if disp_more_every != 0
        time_disp = time()
    end

    if save_each != 0
        time_save = time()
    end

    firstline = 1
    realization_counter = 0
    for i = 1:length(gammas)
        gamma1 = gammas[i]

        #zeroing data arrays
        HB = zeros(zmax)
        HB_x = zeros(zmax)
        BR = zeros(zmax) #ratio of balanced triads
        HB_attr = zeros(zmax)
        HB_only_weights = zeros(zmax)
        paradise = zeros(zmax)
        hell = zeros(zmax)
        Deltas = zeros(4, zmax)
        weak_balance_in_complete_graph = zeros(zmax)
        local_polarization = zeros(zmax)
        global_polarization = zeros(zmax)
        initial_neg_links_count = zeros(zmax)
        links_destab_changed = zeros(4, zmax)

        sim = zeros(zmax)
        x_attr_sim = zeros(zmax)
        stab = zeros(zmax)
        times = zeros(zmax)
        for rep = 1:zmax
            #simulation
            (ishb_sim_par, t, u, u0, xy_attr, sol) =
                calc_heider_attr(n, attr, gamma1, maxtime, ode_fun, solver, false)
            realization_counter += 1

            #work on results
            HB_x[rep],
            HB_attr[rep],
            x_attr_sim[rep],
            BR[rep],
            paradise[rep],
            hell[rep],
            Deltas[:, rep],
            weak_balance_in_complete_graph[rep],
            local_polarization[rep],
            global_polarization[rep] = ishb_sim_par

            initial_neg_links_count[rep] = sum(u0 .< 0)
            links_destab_changed[3, rep] = sum(u[u0.>0] .< 0) #number of initial pos links that changed to negative
            links_destab_changed[4, rep] = sum(u[u0.<0] .> 0) #number of initial neg links that changed to positive

            if t < maxtime #we have stability
                HB[rep] = ishb_sim_par[1]
                stab[rep] = 1
            end

            times[rep] = t

            #checking time/repetitions and eventually informing about progress and saving
            if disp_more_every != 0
                if disp_more_every < time() - time_disp
                    time_disp = time()
                    # displaying
                    g = attr.g
                    v = r.attr_degeneracy
                    a = r.attr_name
                    # display_res(@ntuple(rep, a, n, g, v, gamma1))
                    display_res(; rep, a, n, g, v, gamma1)
                end
            end

            if save_each != 0
                if save_each < time() - time_save
                    time_save = time()
                    # partial saving
                    fields = (
                        HB,
                        HB_x,
                        HB_attr,
                        sim,
                        x_attr_sim,
                        BR,
                        paradise,
                        hell,
                        initial_neg_links_count,
                        links_destab_changed,
                        Deltas,
                        weak_balance_in_complete_graph,
                        local_polarization,
                        global_polarization,
                        stab,
                        times,
                        i,
                        rep,
                        firstline,
                    )
                    update_result!(r, fields)
                    save_result(r, filename)
                    save_result(r, filename, ext = "jld2")
                end
            end

            if disp_each != 0
                if disp_each <= realization_counter / zmax
                    realization_counter = 0
                    # displaying
                    g = attr.g
                    v = r.attr_degeneracy
                    a = r.attr_name
                    # display_res(@ntuple(rep, a, n, g, v, gamma1))
                    display_res(; rep, a, n, g, v, gamma1)
                end
            end
        end
        fields = (
            HB,
            HB_x,
            HB_attr,
            sim,
            x_attr_sim,
            BR,
            paradise,
            hell,
            initial_neg_links_count,
            links_destab_changed,
            Deltas,
            weak_balance_in_complete_graph,
            local_polarization,
            global_polarization,
            stab,
            times,
            i,
            zmax,
            firstline,
        )
        update_result!(r, fields)
        firstline += 2
    end

    save_result(r, filename)
    save_result(r, filename, ext = "jld2")

    return r
end
export using_heider_attr

# Function using `calc_heider_attr` function to simulate a number of repetitions
# of solving the system having given parameters and almost a balanced RL.
# Agents are divided into two groups and form positive links (+0.99) inside those groups
# and negative (-0.99) to agents outside their group.
# Thus, this simulates destabilization of the polarized, initial state.
#
# Parameters:
# n - number of agents
# attr - attribute type; contains all information about attributes
# gammas - a vector of coupling strengths to simulate
# zmax - number of repetitions
# maxtime - maximal time of calculating differential equations
# ode_fun_name - function name calculating derivatives (`string`)
#
# Additional named arguments are:
# disp_each - display status every how many ratio of realizations (if 0 then do not display) (default 0.5)
# disp_more_every - display status every certain number of seconds (if 0, then do not display) (default 600)
# save_each - how often (in seconds) should results be saved (if 0, then save at the end) (default 600)
# files_folder - folder name inside the project the results should be saved (default "data")
#       This can be also an array of folder names in the correct folder structure.
# filename_prefix - string that should start each simulation results' file (default "")
#
# `disp_more_every` and `save_each` work like that, that if the specified time has past
# then sth is displayed/saved. But it doesn't mean exact time of action.
# When the specified time has past means when the current function
# (heider function with parameters) finishes.
#
# Returns `Result` object.
function using_heider_attr_destab(
    n::Int,
    attr::AbstractAttributes,
    gammas::Vector{Float64},
    larger_size::Int,
    zmax::Int,
    maxtime::Float64,
    ode_fun_name::String;
    disp_each = 0.5,
    disp_more_every = 600,
    save_each = 600,
    files_folder::Vector{String} = ["data"],
    filename_prefix::String = "",
)

    ode_fun = getfield(PolarizationFramework, Symbol(ode_fun_name))
    solver = AutoTsit5(Rodas5(autodiff = false))

    rl_weights = init_random_balanced_relations(n, larger_size)

    r, filename = initialize_file(
        n,
        attr,
        gammas,
        maxtime,
        ode_fun_name,
        files_folder,
        filename_prefix,
    )

    # time prepariation
    if disp_more_every != 0
        time_disp = time()
    end

    if save_each != 0
        time_save = time()
    end

    firstline = 1
    realization_counter = 0
    for i = 1:length(gammas)
        gamma1 = gammas[i]

        #zeroing data arrays
        HB = zeros(zmax)
        HB_x = zeros(zmax)
        BR = zeros(zmax) #ratio of balanced triads
        HB_attr = zeros(zmax)
        # HB_only_weights = zeros(zmax);
        paradise = zeros(zmax)
        hell = zeros(zmax)
        Deltas = zeros(4, zmax)
        weak_balance_in_complete_graph = zeros(zmax)
        local_polarization = zeros(zmax)
        global_polarization = zeros(zmax)
        initial_neg_links_count = ones(zmax) * larger_size * (n - larger_size)
        links_destab_changed = zeros(4, zmax)

        sim = zeros(zmax)
        x_attr_sim = zeros(zmax)
        stab = zeros(zmax)
        times = zeros(zmax)

        for rep = 1:zmax
            val0_attr = get_attributes(attr, n)
            al_weights = get_attribute_layer_weights(attr, val0_attr)

            (pos_destab, neg_destab) =
                get_destabilized_links_count(rl_weights, al_weights, gamma1)
            links_destab_changed[1, :] .= pos_destab
            links_destab_changed[2, :] .= neg_destab
            if pos_destab + neg_destab == 0 #no destabilization => no sense of further calculations
                HB[rep] = 1
                HB_x[rep] = 1
                BR[rep] = 1
                HB_attr[rep] = is_hb(al_weights, n)
                # HB_only_weights = zeros(zmax);
                paradise[rep] = is_paradise(rl_weights, n)
                hell[rep] = 0
                Deltas[:, rep] = get_triad_counts(rl_weights, n)
                local_polarization[rep] = get_local_polarization(Deltas[:, rep])
                global_polarization[rep] = 1 - paradise[rep]
                weak_balance_in_complete_graph[rep] = 1

                x_attr_sim[rep] = get_similarity(rl_weights, al_weights, n)
                stab[rep] = 1
                times[rep] = 0

                links_destab_changed[3, rep] = 0
                links_destab_changed[4, rep] = 0
            else
                #simulation
                (ishb_sim_par, t, u, u0, xy_attr, sol) = calc_heider_attr(
                    n,
                    attr,
                    gamma1,
                    maxtime,
                    ode_fun,
                    solver,
                    false,
                    rl_weights,
                    al_weights,
                )

                #work on results
                HB_x[rep],
                HB_attr[rep],
                x_attr_sim[rep],
                BR[rep],
                paradise[rep],
                hell[rep],
                Deltas[:, rep],
                weak_balance_in_complete_graph[rep],
                local_polarization[rep],
                global_polarization[rep] = ishb_sim_par

                links_destab_changed[3, rep] = sum(u[u0.>0] .< 0) #number of initial pos links that changed to negative
                links_destab_changed[4, rep] = sum(u[u0.<0] .> 0) #number of initial neg links that changed to positive

                if t < maxtime #we have stability
                    HB[rep] = ishb_sim_par[1]
                    stab[rep] = 1
                end

                times[rep] = t
            end
            realization_counter += 1

            #checking time/repetitions and eventually informing about progress and saving
            if disp_more_every != 0
                if disp_more_every < time() - time_disp
                    time_disp = time()
                    # displaying
                    # display_res(gamma1, rep)
                    g = attr.g
                    v = r.attr_degeneracy
                    a = r.attr_name
                    # display_res(@ntuple(rep, a, n, g, v, gamma1, larger_size))
                    display_res(; rep, a, n, g, v, gamma1, larger_size)
                end
            end

            if save_each != 0
                if save_each < time() - time_save
                    time_save = time()
                    # partial saving
                    fields = (
                        HB,
                        HB_x,
                        HB_attr,
                        sim,
                        x_attr_sim,
                        BR,
                        paradise,
                        hell,
                        initial_neg_links_count,
                        links_destab_changed,
                        Deltas,
                        weak_balance_in_complete_graph,
                        local_polarization,
                        global_polarization,
                        stab,
                        times,
                        i,
                        rep,
                        firstline,
                    )
                    update_result!(r, fields)
                    save_result(r, filename)
                    save_result(r, filename, ext = "jld2")
                end
            end

            if disp_each != 0
                if disp_each <= realization_counter / zmax
                    realization_counter = 0
                    # displaying
                    g = attr.g
                    v = r.attr_degeneracy
                    a = r.attr_name
                    # display_res(@ntuple(rep, a, n, g, v, gamma1, larger_size))
                    display_res(; rep, a, n, g, v, gamma1, larger_size)
                end
            end
        end
        fields = (
            HB,
            HB_x,
            HB_attr,
            sim,
            x_attr_sim,
            BR,
            paradise,
            hell,
            initial_neg_links_count,
            links_destab_changed,
            Deltas,
            weak_balance_in_complete_graph,
            local_polarization,
            global_polarization,
            stab,
            times,
            i,
            zmax,
            firstline,
        )
        update_result!(r, fields)
        firstline += 2
    end

    save_result(r, filename)
    save_result(r, filename, ext = "jld2")

    return r
end
export using_heider_attr_destab

# end

max_stars_repo_name: pjgorski/PolarizationFramework

max_issues_repo_name: pjgorski/PolarizationFramework

max_forks_repo_name: pjgorski/PolarizationFramework


text:
function cudnnPoolingForward(x, y; window=2, padding=0, stride=1, mode=0)
    @assert padding==0 && stride==1 && mode==0
    # x: (W,H,C,N)
    y2 = zeros(y)
    Wx,Hx,C,Nx = size(x);
    Wy,Hy,K,Ny = size(y);
    @assert Nx == Ny && C==K
    for n in 1:Nx, c in 1:C, j in 1:Hx-psize+1, i in 1:Wx-psize+1
        y2[i,j,c,n] = maximum(x[i:i+psize-1,j:j+psize-1,c,n])
    end
    return y2
end

function cudnnPoolingBackward(y, dy, x, dx; window=2, padding=0, stride=1, mode=0)
    @assert padding==0 && stride==1 && mode==0
    # x: (W,H,C,N)
    dx2 = zeros(x)
    Wx,Hx,C,Nx = size(x);
    Wy,Hy,K,Ny = size(y);
    @assert Nx == Ny && C==K
    for n in 1:Nx, c in 1:C, j in 1:Hx-window+1, i in 1:Wx-window+1
        # y2[i,j,c,n] = maximum()
        a = x[i:i+window-1,j:j+window-1,c,n]
        mi,mj = ind2sub(a,indmax(a))
        dx2[i+mi-1,j+mj-1,c,n] = dy[i,j,c,n]
    end
    return dx2
end


text:
using ThreadingUtilities
using Test
using VectorizationBase

import Aqua
import InteractiveUtils


text:
#=
Online plotting of new spikes in non-paused display
=#
function draw_spike(rhd::RHD2000,han::Gui_Handles)

    spike_num=han.sc.spike
    s=han.sc.s
    o=han.sc.o
    reads=han.draws

    ctx=copy(han.sc.ctx2s)
    paint_with_alpha(ctx,0.0)

    Cairo.translate(ctx,0.0,han.sc.h2/2)
    Cairo.scale(ctx,han.sc.w2/han.wave_points,s)

    for i=1:rhd.nums[spike_num]

        kk=rhd.buf[i,spike_num].inds.start
        k_end=rhd.buf[i,spike_num].inds.stop
        if (kk>0) & (k_end<size(rhd.v,1))

            move_to(ctx,1,rhd.v[kk,spike_num]-o)

            for k=2:han.wave_points
                y = rhd.v[kk+k-1,spike_num]-o
                line_to(ctx,k,y)
            end

	    set_line_width(ctx,0.5);
	    @inbounds select_color(ctx,rhd.buf[i,spike_num].id)
	    stroke(ctx)

            #Add spike to buffer
            if han.buf.count > 0
                mycount=1
                for k=1:han.wave_points
                    han.buf.spikes[mycount,han.buf.ind] = rhd.v[kk+k-1,spike_num]
                    mycount+=1
                end
                han.buf.clus[han.buf.ind] = rhd.buf[i,spike_num].id-1
                han.buf.count+=1
                han.buf.ind+=1
                if han.buf.count>500
                    han.buf.count=500
                end
                if han.buf.ind>500
                    han.buf.ind=1
                end
            end

            update_isi(rhd,han,i)
        end
    end

    identity_matrix(ctx)

    set_source(han.sc.ctx2s,ctx)
    SpikeSorting.mask_surface(han.sc.ctx2s,ctx,0.0,0.0)
    fill(han.sc.ctx2s)

    set_source(han.sc.ctx2,ctx)
    SpikeSorting.mask_surface(han.sc.ctx2,ctx,0.0,0.0)
    fill(han.sc.ctx2)

    nothing
end


text:
# Imports
using Random
using LinearAlgebra
using SparseArrays
using Arpack: eigs

# ESN struct
mutable struct ESN
    inSize::Integer # Tamanho da entrada
    outSize::Integer # Tamanho da saída
    resSize::Integer # Tamanho do reservatório
    leaking_rate::AbstractFloat # leaking rate
    sparsity::AbstractFloat
    spectral_radius::AbstractFloat
    W
    Win
    Wout
    states
    x
    Yt

    function ESN(; inSize::Integer=1,
        outSize::Integer=1,
        resSize::Integer=500,
        leaking_rate::AbstractFloat=0.5,
        sparsity::AbstractFloat=0.6,
        spectral_radius::AbstractFloat=0.9)

        Random.seed!(42)
        # pesos da entrada
        Win = rand(Float64, (resSize, 1+inSize)) .- 0.5

        # pesos do reservatório
        W = sprand(Float64,resSize,resSize, 1-sparsity)
        W = W - 0.5*(W .!= 0)

        rhoW = 0

        try
            rhoW = maximum(abs.(eigs(W)[1])) # mais rápido
        catch e
            rhoW = maximum(abs.(eigvals(Matrix(W)))) #funciona sempre
        end

        W = W * spectral_radius / rhoW

        esn = new()
        esn.inSize = inSize
        esn.outSize = outSize
        esn.resSize = resSize
        esn.leaking_rate = leaking_rate
        esn.sparsity = sparsity
        esn.spectral_radius = spectral_radius
        esn.Win = Win
        esn.W = W

        return esn
    end
end

function fit!(esn::ESN, data)
    # matriz de estados
    states = zeros(Float64, (1+esn.inSize+esn.resSize, length(data)-1))

    # matriz "target"
    Yt = (data[(2):(length(data))])'

    # rodando o reservatório e armazenando valores na matriz de estados
    x = zeros(Float64, (esn.resSize, 1))

    for t = 1:length(data)-1
        u = data[t]
        x = (1-esn.leaking_rate)*x + esn.leaking_rate*tanh.( esn.Win*[1;u] + esn.W*x )
        states[:, t] = [1;u;x]
    end

    reg = 1e-8  # regularization

    esn.Wout = Yt*states' * pinv(states*states' + reg*I(1+esn.inSize+esn.resSize))

    esn.states = states
    esn.x = x
    esn.Yt = Yt

end

function predict!(esn::ESN, steps, v1)
    Y = zeros(esn.outSize, steps)
    u = esn.Yt[length(esn.Yt)]
    x = esn.x
    for t = 1:steps
        x = (1-esn.leaking_rate)*x + esn.leaking_rate*tanh.( esn.Win*[1;u] + esn.W*x )
        y = esn.Wout*[1;u;x]
        Y[:,t] = y
        # generative mode:
        u = y
    end

    return Y
end


text:
"""
An abstract type for which all concrete types expose a database-like
interface.

**Common methods**

An AbstractDataTable is a two-dimensional table with Symbols for
column names. An AbstractDataTable is also similar to an Associative
type in that it allows indexing by a key (the columns).

The following are normally implemented for AbstractDataTables:

* [`describe`](@ref) : summarize columns
* [`dump`](@ref) : show structure
* `hcat` : horizontal concatenation
* `vcat` : vertical concatenation
* `names` : columns names
* [`names!`](@ref) : set columns names
* [`rename!`](@ref) : rename columns names based on keyword arguments
* [`eltypes`](@ref) : `eltype` of each column
* `length` : number of columns
* `size` : (nrows, ncols)
* [`head`](@ref) : first `n` rows
* [`tail`](@ref) : last `n` rows
* `convert` : convert to an array
* `NullableArray` : convert to a NullableArray
* [`completecases`](@ref) : boolean vector of complete cases (rows with no nulls)
* [`dropnull`](@ref) : remove rows with null values
* [`dropnull!`](@ref) : remove rows with null values in-place
* [`nonunique`](@ref) : indexes of duplicate rows
* [`unique!`](@ref) : remove duplicate rows
* `similar` : a DataTable with similar columns as `d`

**Indexing**

Table columns are accessed (`getindex`) by a single index that can be
a symbol identifier, an integer, or a vector of each. If a single
column is selected, just the column object is returned. If multiple
columns are selected, some AbstractDataTable is returned.

```julia
d[:colA]
d[3]
d[[:colA, :colB]]
d[[1:3; 5]]
```

Rows and columns can be indexed like a `Matrix` with the added feature
of indexing columns by name.

```julia
d[1:3, :colA]
d[3,3]
d[3,:]
d[3,[:colA, :colB]]
d[:, [:colA, :colB]]
d[[1:3; 5], :]
```

`setindex` works similarly.
"""
@compat abstract type AbstractDataTable end

##############################################################################
##
## Interface (not final)
##
##############################################################################

# index(dt) => AbstractIndex
# nrow(dt) => Int
# ncol(dt) => Int
# getindex(...)
# setindex!(...) exclusive of methods that add new columns

##############################################################################
##
## Basic properties of a DataTable
##
##############################################################################

immutable Cols{T <: AbstractDataTable} <: AbstractVector{Any}
    dt::T
end
Base.start(::Cols) = 1
Base.done(itr::Cols, st) = st > length(itr.dt)
Base.next(itr::Cols, st) = (itr.dt[st], st + 1)
Base.length(itr::Cols) = length(itr.dt)
Base.size(itr::Cols, ix) = ix==1 ? length(itr) : throw(ArgumentError("Incorrect dimension"))
Base.size(itr::Cols) = (length(itr.dt),)
@compat Base.IndexStyle(::Type{<:Cols}) = IndexLinear()
Base.getindex(itr::Cols, inds...) = getindex(itr.dt, inds...)

# N.B. where stored as a vector, 'columns(x) = x.vector' is a bit cheaper
columns{T <: AbstractDataTable}(dt::T) = Cols{T}(dt)

Base.names(dt::AbstractDataTable) = names(index(dt))
_names(dt::AbstractDataTable) = _names(index(dt))

"""
Set column names

```julia
names!(dt::AbstractDataTable, vals)
```

**Arguments**

* `dt` : the AbstractDataTable
* `vals` : column names, normally a Vector{Symbol} the same length as
  the number of columns in `dt`
* `allow_duplicates` : if `false` (the default), an error will be raised
  if duplicate names are found; if `true`, duplicate names will be suffixed
  with `_i` (`i` starting at 1 for the first duplicate).

**Result**

* `::AbstractDataTable` : the updated result

**Examples**

```julia
dt = DataTable(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
names!(dt, [:a, :b, :c])
names!(dt, [:a, :b, :a])  # throws ArgumentError
names!(dt, [:a, :b, :a], allow_duplicates=true)  # renames second :a to :a_1
```

"""
function names!(dt::AbstractDataTable, vals; allow_duplicates=false)
    names!(index(dt), vals; allow_duplicates=allow_duplicates)
    return dt
end

function rename!(dt::AbstractDataTable, args...)
    rename!(index(dt), args...)
    return dt
end
rename!(f::Function, dt::AbstractDataTable) = rename!(dt, f)

rename(dt::AbstractDataTable, args...) = rename!(copy(dt), args...)
rename(f::Function, dt::AbstractDataTable) = rename(dt, f)

"""
Rename columns

```julia
rename!(dt::AbstractDataTable, from::Symbol, to::Symbol)
rename!(dt::AbstractDataTable, d::Associative)
rename!(f::Function, dt::AbstractDataTable)
rename(dt::AbstractDataTable, from::Symbol, to::Symbol)
rename(f::Function, dt::AbstractDataTable)
```

**Arguments**

* `dt` : the AbstractDataTable
* `d` : an Associative type that maps the original name to a new name
* `f` : a function that has the old column name (a symbol) as input
  and new column name (a symbol) as output

**Result**

* `::AbstractDataTable` : the updated result

**Examples**

```julia
dt = DataTable(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
rename(x -> @compat(Symbol)(uppercase(string(x))), dt)
rename(dt, @compat(Dict(:i=>:A, :x=>:X)))
rename(dt, :y, :Y)
rename!(dt, @compat(Dict(:i=>:A, :x=>:X)))
```

"""
(rename!, rename)

"""
Return element types of columns

```julia
eltypes(dt::AbstractDataTable)
```

**Arguments**

* `dt` : the AbstractDataTable

**Result**

* `::Vector{Type}` : the element type of each column

**Examples**

```julia
dt = DataTable(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
eltypes(dt)
```

"""
eltypes(dt::AbstractDataTable) = map!(eltype, Vector{Type}(size(dt,2)), columns(dt))

Base.size(dt::AbstractDataTable) = (nrow(dt), ncol(dt))
function Base.size(dt::AbstractDataTable, i::Integer)
    if i == 1
        nrow(dt)
    elseif i == 2
        ncol(dt)
    else
        throw(ArgumentError("DataTables only have two dimensions"))
    end
end

Base.length(dt::AbstractDataTable) = ncol(dt)
Base.endof(dt::AbstractDataTable) = ncol(dt)

Base.ndims(::AbstractDataTable) = 2

##############################################################################
##
## Similar
##
##############################################################################

Base.similar(dt::AbstractDataTable, dims::Int) =
    DataTable(Any[similar(x, dims) for x in columns(dt)], copy(index(dt)))

##############################################################################
##
## Equality
##
##############################################################################

# Imported in DataTables.jl for compatibility across Julia 0.4 and 0.5
@compat(Base.:(==))(dt1::AbstractDataTable, dt2::AbstractDataTable) = isequal(dt1, dt2)

function Base.isequal(dt1::AbstractDataTable, dt2::AbstractDataTable)
    size(dt1, 2) == size(dt2, 2) || return false
    isequal(index(dt1), index(dt2)) || return false
    for idx in 1:size(dt1, 2)
        isequal(dt1[idx], dt2[idx]) || return false
    end
    return true
end

##############################################################################
##
## Associative methods
##
##############################################################################

Base.haskey(dt::AbstractDataTable, key::Any) = haskey(index(dt), key)
Base.get(dt::AbstractDataTable, key::Any, default::Any) = haskey(dt, key) ? dt[key] : default
Base.isempty(dt::AbstractDataTable) = ncol(dt) == 0

##############################################################################
##
## Description
##
##############################################################################

head(dt::AbstractDataTable, r::Int) = dt[1:min(r,nrow(dt)), :]
head(dt::AbstractDataTable) = head(dt, 6)
tail(dt::AbstractDataTable, r::Int) = dt[max(1,nrow(dt)-r+1):nrow(dt), :]
tail(dt::AbstractDataTable) = tail(dt, 6)

"""
Show the first or last part of an AbstractDataTable

```julia
head(dt::AbstractDataTable, r::Int = 6)
tail(dt::AbstractDataTable, r::Int = 6)
```

**Arguments**

* `dt` : the AbstractDataTable
* `r` : the number of rows to show

**Result**

* `::AbstractDataTable` : the first or last part of `dt`

**Examples**

```julia
dt = DataTable(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
head(dt)
tail(dt)
```

"""
(head, tail)

# get the structure of a DT
"""
Show the structure of an AbstractDataTable, in a tree-like format

```julia
dump(dt::AbstractDataTable, n::Int = 5)
dump(io::IO, dt::AbstractDataTable, n::Int = 5)
```

**Arguments**

* `dt` : the AbstractDataTable
* `n` : the number of levels to show
* `io` : optional output descriptor

**Result**

* nothing

**Examples**

```julia
dt = DataTable(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
dump(dt)
```

"""
function Base.dump(io::IO, dt::AbstractDataTable, n::Int, indent)
    println(io, typeof(dt), "  $(nrow(dt)) observations of $(ncol(dt)) variables")
    if n > 0
        for (name, col) in eachcol(dt)
            print(io, indent, "  ", name, ": ")
            dump(io, col, n - 1, string(indent, "  "))
        end
    end
end

# summarize the columns of a DT
# TODO: clever layout in rows
"""
Summarize the columns of an AbstractDataTable

```julia
describe(dt::AbstractDataTable)
describe(io, dt::AbstractDataTable)
```

**Arguments**

* `dt` : the AbstractDataTable
* `io` : optional output descriptor

**Result**

* nothing

**Details**

If the column's base type derives from Number, compute the minimum, first
quantile, median, mean, third quantile, and maximum. Nulls are filtered and
reported separately.

For boolean columns, report trues, falses, and nulls.

For other types, show column characteristics and number of nulls.

**Examples**

```julia
dt = DataTable(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
describe(dt)
```

"""
StatsBase.describe(dt::AbstractDataTable) = describe(STDOUT, dt)
function StatsBase.describe(io, dt::AbstractDataTable)
    for (name, col) in eachcol(dt)
        println(io, name)
        describe(io, col)
        println(io, )
    end
end

##############################################################################
##
## Miscellaneous
##
##############################################################################

function _nonnull!(res, col)
    for (i, el) in enumerate(col)
        res[i] &= !_isnull(el)
    end
end

function _nonnull!(res, col::NullableArray)
    for (i, el) in enumerate(col.isnull)
        res[i] &= !el
    end
end

function _nonnull!(res, col::NullableCategoricalArray)
    for (i, el) in enumerate(col.refs)
        res[i] &= el > 0
    end
end

"""
Indexes of complete cases (rows without null values)

```julia
completecases(dt::AbstractDataTable)
```

**Arguments**

* `dt` : the AbstractDataTable

**Result**

* `::Vector{Bool}` : indexes of complete cases

See also [`dropnull`](@ref) and [`dropnull!`](@ref).

**Examples**

```julia
dt = DataTable(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
dt[[1,4,5], :x] = Nullable()
dt[[9,10], :y] = Nullable()
completecases(dt)
```

"""
function completecases(dt::AbstractDataTable)
    res = trues(size(dt, 1))
    for i in 1:size(dt, 2)
        _nonnull!(res, dt[i])
    end
    res
end

"""
Remove rows with null values.

```julia
dropnull(dt::AbstractDataTable)
```

**Arguments**

* `dt` : the AbstractDataTable

**Result**

* `::AbstractDataTable` : the updated copy

See also [`completecases`](@ref) and [`dropnull!`](@ref).

**Examples**

```julia
dt = DataTable(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
dt[[1,4,5], :x] = Nullable()
dt[[9,10], :y] = Nullable()
dropnull(dt)
```

"""
dropnull(dt::AbstractDataTable) = deleterows!(copy(dt), find(!, completecases(dt)))

"""
Remove rows with null values in-place.

```julia
dropnull!(dt::AbstractDataTable)
```

**Arguments**

* `dt` : the AbstractDataTable

**Result**

* `::AbstractDataTable` : the updated version

See also [`dropnull`](@ref) and [`completecases`](@ref).

**Examples**

```julia
dt = DataTable(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
dt[[1,4,5], :x] = Nullable()
dt[[9,10], :y] = Nullable()
dropnull!(dt)
```

"""
dropnull!(dt::AbstractDataTable) = deleterows!(dt, find(!, completecases(dt)))

function Base.convert(::Type{Array}, dt::AbstractDataTable)
    convert(Matrix, dt)
end
function Base.convert(::Type{Matrix}, dt::AbstractDataTable)
    T = reduce(promote_type, eltypes(dt))
    T <: Nullable && (T = eltype(T))
    convert(Matrix{T}, dt)
end
function Base.convert{T}(::Type{Array{T}}, dt::AbstractDataTable)
    convert(Matrix{T}, dt)
end
function Base.convert{T}(::Type{Matrix{T}}, dt::AbstractDataTable)
    n, p = size(dt)
    res = Matrix{T}(n, p)
    idx = 1
    for (name, col) in zip(names(dt), columns(dt))
        any(isnull, col) && error("cannot convert a DataTable containing null values to array (found for column $name)")
        copy!(res, idx, convert(Vector{T}, col))
        idx += n
    end
    return res
end

function Base.convert(::Type{NullableArray}, dt::AbstractDataTable)
    convert(NullableMatrix, dt)
end
function Base.convert(::Type{NullableMatrix}, dt::AbstractDataTable)
    T = reduce(promote_type, eltypes(dt))
    T <: Nullable && (T = eltype(T))
    convert(NullableMatrix{T}, dt)
end
function Base.convert{T}(::Type{NullableArray{T}}, dt::AbstractDataTable)
    convert(NullableMatrix{T}, dt)
end
function Base.convert{T}(::Type{NullableMatrix{T}}, dt::AbstractDataTable)
    n, p = size(dt)
    res = NullableArray(T, n, p)
    idx = 1
    for col in columns(dt)
        copy!(res, idx, col)
        idx += n
    end
    return res
end

"""
Indexes of duplicate rows (a row that is a duplicate of a prior row)

```julia
nonunique(dt::AbstractDataTable)
nonunique(dt::AbstractDataTable, cols)
```

**Arguments**

* `dt` : the AbstractDataTable
* `cols` : a column indicator (Symbol, Int, Vector{Symbol}, etc.) specifying the column(s) to compare

**Result**

* `::Vector{Bool}` : indicates whether the row is a duplicate of some
  prior row

See also [`unique`](@ref) and [`unique!`](@ref).

**Examples**

```julia
dt = DataTable(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
dt = vcat(dt, dt)
nonunique(dt)
nonunique(dt, 1)
```

"""
function nonunique(dt::AbstractDataTable)
    gslots = row_group_slots(dt)[3]
    # unique rows are the first encountered group representatives,
    # nonunique are everything else
    res = fill(true, nrow(dt))
    @inbounds for g_row in gslots
        (g_row > 0) && (res[g_row] = false)
    end
    return res
end

nonunique(dt::AbstractDataTable, cols::Union{Real, Symbol}) = nonunique(dt[[cols]])
nonunique(dt::AbstractDataTable, cols::Any) = nonunique(dt[cols])

if isdefined(Base, :unique!) # Julia >= 0.7
    import Base.unique!
end

unique!(dt::AbstractDataTable) = deleterows!(dt, find(nonunique(dt)))
unique!(dt::AbstractDataTable, cols::Any) = deleterows!(dt, find(nonunique(dt, cols)))

# Unique rows of an AbstractDataTable.
Base.unique(dt::AbstractDataTable) = dt[(!).(nonunique(dt)), :]
Base.unique(dt::AbstractDataTable, cols::Any) = dt[(!).(nonunique(dt, cols)), :]

"""
Delete duplicate rows

```julia
unique(dt::AbstractDataTable)
unique(dt::AbstractDataTable, cols)
unique!(dt::AbstractDataTable)
unique!(dt::AbstractDataTable, cols)
```

**Arguments**

* `dt` : the AbstractDataTable
* `cols` :  column indicator (Symbol, Int, Vector{Symbol}, etc.)
specifying the column(s) to compare.

**Result**

* `::AbstractDataTable` : the updated version of `dt` with unique rows.
When `cols` is specified, the return DataTable contains complete rows,
retaining in each case the first instance for which `dt[cols]` is unique.

See also [`nonunique`](@ref).

**Examples**

```julia
dt = DataTable(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
dt = vcat(dt, dt)
unique(dt)   # doesn't modify dt
unique(dt, 1)
unique!(dt)  # modifies dt
```

"""
(unique, unique!)

# Count the number of missing values in every column of an AbstractDataTable.
function colmissing(dt::AbstractDataTable) # -> Vector{Int}
    nrows, ncols = size(dt)
    missing = zeros(Int, ncols)
    for j in 1:ncols
        missing[j] = countnull(dt[j])
    end
    return missing
end

function without(dt::AbstractDataTable, icols::Vector{Int})
    newcols = setdiff(1:ncol(dt), icols)
    dt[newcols]
end
without(dt::AbstractDataTable, i::Int) = without(dt, [i])
without(dt::AbstractDataTable, c::Any) = without(dt, index(dt)[c])

##############################################################################
##
## Hcat / vcat
##
##############################################################################

# hcat's first argument must be an AbstractDataTable
# Trailing arguments (currently) may also be NullableVectors, Vectors, or scalars.

# hcat! is defined in datatables/datatables.jl
# Its first argument (currently) must be a DataTable.

# catch-all to cover cases where indexing returns a DataTable and copy doesn't
Base.hcat(dt::AbstractDataTable, x) = hcat!(dt[:, :], x)
Base.hcat(dt1::AbstractDataTable, dt2::AbstractDataTable) = hcat!(dt1[:, :], dt2)

Base.hcat(dt::AbstractDataTable, x, y...) = hcat!(hcat(dt, x), y...)
Base.hcat(dt1::AbstractDataTable, dt2::AbstractDataTable, dtn::AbstractDataTable...) = hcat!(hcat(dt1, dt2), dtn...)

@generated function promote_col_type(cols::AbstractVector...)
    elty = Base.promote_eltype(cols...)
    if elty <: Nullable
        elty = eltype(elty)
    end
    if elty <: CategoricalValue
        elty = elty.parameters[1]
    end
    if any(col -> eltype(col) <: Nullable, cols)
        if any(col -> col <: Union{AbstractCategoricalArray, AbstractNullableCategoricalArray}, cols)
            return :(NullableCategoricalVector{$elty})
        else
            return :(NullableVector{$elty})
        end
    else
        if any(col -> col <: Union{AbstractCategoricalArray, AbstractNullableCategoricalArray}, cols)
            return :(CategoricalVector{$elty})
        else
            return :(Vector{$elty})
        end
    end
end

"""
    vcat(dts::AbstractDataTable...)

Vertically concatenate `AbstractDataTables` that have the same column names in
the same order.

# Example
```jldoctest
julia> dt1 = DataTable(A=1:3, B=1:3);
julia> dt2 = DataTable(A=4:6, B=4:6);
julia> vcat(dt1, dt2)
6×2 DataTables.DataTable
│ Row │ A │ B │
├─────┼───┼───┤
│ 1   │ 1 │ 1 │
│ 2   │ 2 │ 2 │
│ 3   │ 3 │ 3 │
│ 4   │ 4 │ 4 │
│ 5   │ 5 │ 5 │
│ 6   │ 6 │ 6 │
```
"""
Base.vcat(dt::AbstractDataTable) = dt
function Base.vcat(dts::AbstractDataTable...)
    isempty(dts) && return DataTable()
    allheaders = map(names, dts)
    if all(h -> length(h) == 0, allheaders)
        return DataTable()
    end
    uniqueheaders = unique(allheaders)
    if length(uniqueheaders) > 1
        unionunique = union(uniqueheaders...)
        coldiff = setdiff(unionunique, intersect(uniqueheaders...))
        if !isempty(coldiff)
            # if any datatables are a full superset of names, skip them
            filter!(u -> Set(u) != Set(unionunique), uniqueheaders)
            estrings = Vector{String}(length(uniqueheaders))
            for (i, u) in enumerate(uniqueheaders)
                matching = find(h -> u == h, allheaders)
                headerdiff = setdiff(coldiff, u)
                cols = join(headerdiff, ", ", " and ")
                args = join(matching, ", ", " and ")
                estrings[i] = "column(s) $cols are missing from argument(s) $args"
            end
            throw(ArgumentError(join(estrings, ", ", ", and ")))
        else
            estrings = Vector{String}(length(uniqueheaders))
            for (i, u) in enumerate(uniqueheaders)
                indices = find(a -> a == u, allheaders)
                estrings[i] = "column order of argument(s) $(join(indices, ", ", " and "))"
            end
            throw(ArgumentError(join(estrings, " != ")))
        end
    else
        header = uniqueheaders[1]
        cols = Vector{Any}(length(header))
        for i in 1:length(cols)
            data = [dt[i] for dt in dts]
            lens = map(length, data)
            cols[i] = promote_col_type(data...)(sum(lens))
            offset = 1
            for j in 1:length(data)
                copy!(cols[i], offset, data[j])
                offset += lens[j]
            end
        end
        return DataTable(cols, header)
    end
end

##############################################################################
##
## Hashing
##
## Make sure this agrees with isequals()
##
##############################################################################

function Base.hash(dt::AbstractDataTable)
    h = hash(size(dt)) + 1
    for i in 1:size(dt, 2)
        h = hash(dt[i], h)
    end
    return @compat UInt(h)
end

## Documentation for methods defined elsewhere

"""
Number of rows or columns in an AbstractDataTable

```julia
nrow(dt::AbstractDataTable)
ncol(dt::AbstractDataTable)
```

**Arguments**

* `dt` : the AbstractDataTable

**Result**

* `::AbstractDataTable` : the updated version

See also [`size`](@ref).

NOTE: these functions may be depreciated for `size`.

**Examples**

```julia
dt = DataTable(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
size(dt)
nrow(dt)
ncol(dt)
```

"""
# nrow, ncol


text:
module JMcDM

using DataFrames
using LinearAlgebra
using JuMP
using Cbc
using StatsBase
using Statistics

# includes
include("types.jl")
include("utilities.jl")
include("topsis.jl")
include("vikor.jl")
include("electre.jl")
include("moora.jl")
include("dematel.jl")
include("ahp.jl")
include("nds.jl")
include("singlecriterion.jl")
include("game.jl")
include("dataenvelop.jl")
include("grey.jl")
include("saw.jl")
include("aras.jl")
include("wpm.jl")
include("waspas.jl")
include("edas.jl")
include("marcos.jl")
include("mabac.jl")
include("mairca.jl")
include("copras.jl")
include("promethee.jl")
include("cocoso.jl")
include("critic.jl")
include("entropy.jl")
include("codas.jl")

include("summary.jl")

# export MCDM types
export MCDMResult
export TopsisResult
export VikorResult
export ElectreResult
export MooraResult
export AHPConsistencyResult
export AHPResult
export NDSResult
export GreyResult
export SawResult
export ARASResult
export WPMResult
export WASPASResult
export EDASResult
export MARCOSResult
export MABACResult
export MAIRCAResult
export COPRASResult
export PrometheeResult
export CoCoSoResult
export CRITICResult
export EntropyResult
export CODASResult

# export game type
export GameResult

#  export SCDM types
export SCDMResult
export LaplaceResult
export MaximinResult
export MaximaxResult
export MinimaxResult
export MiniminResult
export SavageResult
export HurwiczResult
export MLEResult
export ExpectedRegretResult
export DataEnvelopResult

# export utility functions
export euclidean
export normalize
export colmaxs
export colmins
export unitize
export makeDecisionMatrix
export reverseminmax
export makeminmax

#  export MCDM tools
export topsis
export vikor
export electre
export moora
export dematel
export ahp_RI, ahp_consistency, ahp
export nds
export grey
export saw
export aras
export wpm
export waspas
export edas
export marcos
export mabac
export mairca
export copras
export promethee, prometLinear, prometVShape, prometUsual, prometQuasi, prometLevel
export cocoso
export critic
export entropy
export codas

#  export SCDM tools
export laplace
export maximin
export maximax
export minimax
export minimin
export savage
export hurwicz
export mle
export expectedregret

# export game solver
export game

# export data envelop
export dataenvelop

#  export summary function
export summary

end # module

max_stars_repo_name: bahadirfyildirim/JMcDM

max_issues_repo_name: bahadirfyildirim/JMcDM

max_forks_repo_name: bahadirfyildirim/JMcDM


text:
# Abstract representation of demand.

#=
Demands are parameterized in terms of how they are registered in time and where
they are located.  Functions will be provided to convert `Demand` into the
appropriate type `StaticLoad`, `InterruptibleLoad`, etc. that is properly
located at a `Bus`.
=#

# Import packages.

using Dates, TimeSeries

# Timeseries for allowable ranges of demand.

"""
An "envelope" of minimum and maximum allowable demands at each time point.

# Type parameters
- `T <: TimeType`: timestamp

# Example
```
example = TimeArray(
    [Time(0) , Time(8)   , Time(9) , Time(17)  , Time(18)],
    [(5., 8.), (10., 12.), (8., 8.), (11., 13.), (3., 4.)]
) :: Envelope{Time}
```
"""
const Envelope{T <: TimeType} = TimeArray{Tuple{Float64,Float64},1,T,Array{Tuple{Float64,Float64},1}}

"""
An "envelope" of minimum and maximum allowable demands at each time point, with a location for the demand.

# Type parameters
- `T <: TimeType`: timestamp
- `L`            : network location

# Example
```
example = TimeArray(
    [Time(0)             , Time(8)             , Time(9)                 , Time(17)           , Time(18)            ],
    [("Home #23", 5., 8.), ("Road #14", 0., 0.), ("Workplace #3", 8., 8.), ("Road #9", 0., 0.), ("Home #23", 3., 4.)]
) :: LocatedEnvelope{Time}
```
"""
const LocatedEnvelope{T <: TimeType, L} = TimeArray{Tuple{L,Float64,Float64},1,T,Array{Tuple{L,Float64,Float64},1}}

# Abstract interface for demand.

"""
The most abstract type of demand.

# Type parameters
- `T <: TimeType`: timestamp
- `L`            : network location
"""
abstract type Demand{T <: TimeType, L} end

"""
The "envelope" of minimum and maximum allowable demands at each time pont, with a location for the demand.

This must be implemented by subtypes of `Demand`.

# Arguments
- `demand :: Demand{T,L}`: the demand
"""
function envelope(demand :: Demand{T,L}) :: LocatedEnvelope{T, L} where L where T <: TimeType
end


text:
function _extract_tsfresh_feats(column, shift; ka...)
    @from tsfresh imports extract_relevant_features
    @from tsfresh.utilities.dataframe_functions imports roll_time_series
    dst = @sprintf("tsfresh/%s_%s.h5", column, shift)
    isfile(dst) && return dst
    df = pd.read_parquet("df.parquet", columns = ["代码", column])
    if df[column].dtype == "uint8"
        df[column] = df[column] / 128 - 1
    end
    df[column] = df[column].replace([Inf, -Inf], NaN).fillna(0)
    y = pd.read_pickle("y.pkl")["y"]
    if length(y) == length(df)
        df = roll_time_series(df, "代码", nothing, nothing, 1, shift)
        df.drop(columns = "代码", inplace = true)
        y.index = df.groupby("id")["sort"].last().sort_values().index
        df_erf = extract_relevant_features(df, y; column_id = "id", column_sort = "sort", ka...)
    else
        df_erf = extract_relevant_features(df, y; column_id = "代码", ka...)
        df_erf.drop(columns = "代码", inplace = true)
    end
    df_erf = df_erf.astype("float32")
    df_erf.columns = ["$c-$shift" for c in df_erf.columns]
    to_hdf5(df_erf, dst)
    @eval GC.gc()
    return dst
end

function extract_tsfresh_feats(df, y = nothing; shifts = ["20T"], horizon = "20T", remove = true, ka...)
    Δt = df["时间戳"].groupby(df["代码"]).diff().median()
    shifts = @. ceil(Int, parsefreq(shifts) / Δt)
    if isnothing(y)
        horizon = ceil(Int, parsefreq(horizon) / Δt)
        y = df["涨幅"].groupby(df["代码"]).rolling(horizon).sum().groupby("代码").shift(-horizon).fillna(0)
        y = Series(y.reset_index(level = "代码", drop = true).sort_index().astype("float32"), name = "y")
    end
    df_meta, df_fea = split_metafeat(df)
    df_ts = pdhcat(df_fea, df_meta[["代码"]])
    df_ts.to_parquet("df.parquet")
    y.to_frame().to_pickle("y.pkl")
    !isdir("tsfresh") && mkdir("tsfresh")
    h5s = pmap(Iterators.product(df_fea.columns, shifts)) do (column, shift)
        _extract_tsfresh_feats(column, shift; ka...)
    end
    if length(df) != length(y)
        df = df.drop_duplicates(subset = "代码", keep = "last")
        df.reset_index(inplace = true)
        @assert length(df) == length(y)
    end
    df = pdhcat(df, read_hdf5.(h5s)...)
    remove && foreach(rm, ["df.parquet", "y.pkl", h5s...])
    return df
end

function extract_talib_feats(df, bycol = "code"; keep_columns = [])
    ohlc_cols = ["open", "high", "low", "close", "volume"]
    df_ohlc = df[ohlc_cols].astype("double")
    df_talib = df[df.columns.drop(ohlc_cols)]
    @from talib imports get_function_groups
    normclose = [get_function_groups()["Overlap Studies"];
                get_function_groups()["Price Transform"];
                ["APO", "MACD", "MACDEXT", "MACDFIX", "MINUS_DM", "MOM", "PLUS_DM"];
                ["TRANGE", "ATR", "LINEARREG", "LINEARREG_INTERCEPT", "STDDEV", "TSF"]]
    groupfuncs = [(g, f) for (g, fs) in get_function_groups() for f in fs]
    dfs = @showprogress 10 pmap(groupfuncs) do (g, f)
        g = split(g, ' ')[1]
        g == "Math" && return
        f == "MAVP" && return
        f == "VAR" && return
        c = "talib:" * g * "_" * f
        !isempty(keep_columns) && c ∉ keep_columns && return
        println(c)
        py"""
import gc
import pandas as pd
from talib.abstract import Function
def apptalib(x):
    df = Function($f)(x)
    if len(df.shape) == 1:
        df = df.to_frame()
    return df
df = $df_ohlc.groupby($df_talib[$bycol]).apply(apptalib)
if $f in $$normclose:
    for c in df.columns:
        df[c] = df[c] / $df_ohlc["close"]
df = df.astype('float32')
df.columns = [$c if c == 0 else $c + '_' + c  for c in df]
gc.collect()
"""
        return py"df"
    end
    return pdhcat(df_talib, filter(!isnothing, dfs)...)
end

function TradingFCParameters(tscal = 1)
    Dict(
        "abs_energy" => nothing,
        "count_above_mean" => nothing,
        "first_location_of_maximum" => nothing,
        "first_location_of_minimum" => nothing,
        "kurtosis" => nothing,
        "last_location_of_maximum" => nothing,
        "last_location_of_minimum" => nothing,
        "longest_strike_above_mean" => nothing,
        "longest_strike_below_mean" => nothing,
        "maximum" => nothing,
        "mean" => nothing,
        "mean_abs_change" => nothing,
        "mean_change" => nothing,
        "mean_second_derivative_central" => nothing,
        "median" => nothing,
        "minimum" => nothing,
        # "sample_entropy" => nothing,
        "skewness" => nothing,
        "standard_deviation" => nothing,
        "variance" => nothing,
        # "time_reversal_asymmetry_statistic" => [Dict("lag" => tscal * lag) for lag in [1, 5]],
        # "c3" => [Dict("lag" => tscal * lag) for lag in [1, 5]],
        # "symmetry_looking" => [Dict("r" => r * 0.05) for r in [1, 5, 10]],
        # "large_standard_deviation" => [Dict("r" => r * 0.05) for r in [1, 5, 10]],
        # "quantile" => [Dict("q" => q) for q in [0.25, 0.75]],
        # "autocorrelation" => [Dict("lag" => tscal * lag) for lag in [1, 5, 10]],
        # "agg_autocorrelation" => [Dict("f_agg" => s, "maxlag" => 40 * tscal) for s in ["mean", "median", "var"]],
        # "partial_autocorrelation" => [Dict("lag" => tscal * lag) for lag in [1, 5, 10]],
        # !"number_peaks" => [Dict("n" => tscal * n) for n in [1, 5, 10]],
        "binned_entropy" => [Dict("max_bins" => 10)],
        # !"fft_aggregated" => [Dict("aggtype" => s) for s in ["centroid", "variance", "skew", "kurtosis"]],
        # "max_langevin_fixed_point" => [Dict("m" => 3, "r" => 30)],
        # !"linear_trend" => [Dict("attr" => "pvalue"), Dict("attr" => "rvalue"), Dict("attr" => "slope"), Dict("attr" => "stderr")],
        # "agg_linear_trend" => [Dict("attr" => attr, "chunk_len" => tscal * i, "f_agg" => f)
        #                      for attr in ["rvalue", "slope", "stderr"]
        #                      for i in [5, 10, 50]
        #                      for f in ["max", "min", "mean", "var"]],
        # "augmented_dickey_fuller" => [Dict("attr" => "teststat"), Dict("attr" => "pvalue"), Dict("attr" => "usedlag")],
        # "ratio_beyond_r_sigma" => [Dict("r" => x) for x in [1, 5]],
        # "ar_coefficient" => [Dict("coeff" => coeff, "k" => k) for coeff in 0:4 for k in [10]],
    )
end