text:
print("Hello world. walls")

include("PRM.jl")
#using Distributions
glvisualize()

pfont = Plots.Font("sans-serif", 8, :left, :left, 0.0, RGB{U8}(0.0, 0.0, 0.0))
glvisualize(titlefont = pfont, guidefont = pfont)

####################################
# Randomly Generate Clutter (obstacles)
####################################
# implement function to create clutter
# Let's say on average we want to create 4 obstacles... (otherwise, *most* of
# the time we will product one large rectangle that is falling out of the room

#targetNumObs = 2
function genClutter(roomperimeter, param, targetNumObs, clutterPercentage, roomWidth, roomHeight)
    obstacles = Vector{HyperRectangle}()
    roomArea = roomWidth * roomHeight
    targetSumObsArea= roomArea * clutterPercentage
    sumObsArea = 0
    while sumObsArea < targetSumObsArea
        x,y = rand(Uniform(0, roomWidth),2)

        randWidth, randHeight = rand(Uniform(1, roomWidth/targetNumObs),2)
        #randWidth, randHeight = rand(1:roomWidth/targetNumObs,2)
        protoObstacle = HyperRectangle(Vec(x, y), Vec(randWidth, randHeight)) #Todo
        if contains(roomperimeter, protoObstacle)
            push!(obstacles, protoObstacle)
            sumObsArea += randWidth*randHeight
        end
    end
    actualSumObsArea = sumObsArea
    #print("obstacles generated")
    return actualSumObsArea, obstacles
end
####

function genTrapWalls()
    trapWalls = Vector{LineSegment}()
    lineL = LineSegment( Point(5,5), Point(5,10))
    lineR = LineSegment( Point(10,5), Point(10,10))
    #lineBot = LineSegment( Point(5,5), Point(10,5))
    lineTop = LineSegment( Point(5,10), Point(10,10))
    lineBotRight = LineSegment( Point(8,5), (10,5))
    lineBotLeft = LineSegment( Point(5,5), (7,5))
    push!(trapWalls, lineL, lineR, lineTop , lineBotLeft, lineBotRight) #opening at top
    startstate = (7,7) #start state inside the trap
    return startstate, trapWalls
end

function genTrapWalls(doorSize)
    trapWalls = Vector{LineSegment}()
    lineL = LineSegment( Point(5,5), Point(5,10))
    lineR = LineSegment( Point(10,5), Point(10,10))
    #lineBot = LineSegment( Point(5,5), Point(10,5))
    lineTop = LineSegment( Point(5,10), Point(10,10))
    lstop = 7.5 - (doorSize/2.0)
    rstart = lstop + doorSize

    lineBotRight = LineSegment( Point(rstart,5), (10,5))
    lineBotLeft = LineSegment( Point(5,5), (lstop,5))

    push!(trapWalls, lineL, lineR, lineTop , lineBotLeft, lineBotRight) #opening at top
    startstate = (7,7) #start state inside the trap
    return startstate, trapWalls
end

function testGenWalls()

    ####################################
    ## PARAMETERS
    ####################################
    numPts = 100
    connectRadius = 10
    param = algT.AlgParameters(numPts, connectRadius)

    targetNumObs = 3
    clutterPercentage = 0.1
    roomWidth,roomHeight  = 20,20

    startstate = Point(7.,7)
    goalstate = Point(19.,19)

    ####################################
    ## INIT
    ####################################
    obstacles = Vector{HyperRectangle}()
    perimeter = HyperRectangle(Vec(0.,0), Vec(roomWidth,roomHeight))
    wallsPerimeter = algfxn.decompRect(perimeter)

    walls = Vector{LineSegment}()
    for l in wallsPerimeter
        push!(walls, l)
    end

    #lineLeftTop = LineSegment( Point(5,10), Point(7,10))
    #lineRightTop = LineSegment( Point(8,10), Point(10,10))

    #push!(walls, lineL, lineR, lineBot, lineLeftTop, lineRightTop) #opening at bottom
    push!(walls, lineL, lineR, lineTop , lineBotLeft, lineBotRight) #opening at top

    #plot!(walls, color =:black)
    #plot!(obstacles, fillalpha=0.5)

    obsArea, obstacles = genClutter(perimeter, param, targetNumObs, clutterPercentage, roomWidth, roomHeight)

####################################
# Run PRM once and display config space plot
####################################

    #obsArea, obstacles = genClutter(perimeter, param, targetNumObs, clutterPercentage, roomWidth, roomHeight)
    r = algT.Room(roomWidth,roomHeight,walls,obstacles)
    roomPlot = plotfxn.plotRoom(r)

    ## Run preprocessing
    nodeslist, edgeslist = preprocessPRM(r, param)
    print("\n ---- pre-processed --- \n")

    ## Query created RM

    pathcost, isPathFound, solPath = queryPRM(startstate, goalstate, nodeslist, edgeslist, obstacles)
    print("\n ---- queried ---- \n")

    ## Plot path found
    timestamp = Base.Dates.now()
    title = "PRM with numPts =$numPts, maxDist=$connectRadius, \npathcost = $pathcost,
            timestamp=$timestamp)\n\n"
    roadmap = algT.roadmap(startstate, goalstate, nodeslist, edgeslist)

    #prmPlot= plotfxn.plotPRM(roomPlot, roadmap, solPath, title::String, afont)
    prmPlot= plotfxn.plotPRM(roomPlot, roadmap, solPath, title::String)
    plot!(prmPlot, size=(900,600))

    print("\n --- Time --- \n")
    @show timestamp
    #print("\n --- Obstacles --- \n")
    #@show obstacles
    #print("\n --- Nodes --- \n")
    #@show nodeslist
    #print("\n --- Edges --- \n")
    #@show edgeslist
    # print("\n -------- \n")

    gui(prmPlot)
        print("\n --- Press ENTER to continue --- \n")
        #readline(STDIN)
        # todo: save one "representative" figure from each run
        # print("\n --- continue --- \n")
    @show  isPathFound

end

function walls_vs_nSamples()
    ####################################
    # Run multiple trials and scatterplot all cost vs area for all runs
    # due to varying absolute area
    ####################################
    ####################################
    ## PARAMETERS - RRT
    ####################################
    connectRadius = 8.
    numPts = 0.
    param = algT.AlgParameters(numPts, connectRadius)
    targetNumObs = 5
    percentClutter = 0.15
    roomWidth,roomHeight  = 20,20

    startstate, trapwalls = genTrapWalls() #bugtrap walls
    goalstate = Point(19.,19)

    ####################################
    ## PARAMETERS - plot clutter
    ####################################
    N = 30

    #clutterPercentages = [0 0.01 0.05 0.1 0.12 0.15 0.2 0.25 0.3 0.4 0.5]
    #clutterPercentages = [0 0.1 0.2 0.3 0.5]

    numPtsList = [10 20 30 80 100]

    ####################################
    ## Init
    ####################################
    listCosts = Vector{Float32}()
    listpSucc= Vector{Float32}()

    obstacles = Vector{HyperRectangle}()
    perimeter = HyperRectangle(Vec(0.,0), Vec(roomWidth,roomHeight))
    wallsPerimeter = algfxn.decompRect(perimeter)

    walls = Vector{LineSegment}()

    for l in wallsPerimeter
        push!(walls, l)
    end
    for l in trapwalls
        push!(walls, l)
    end

    ####################################
    # Run Trials
    ####################################

    plot()
    tic()

    stddevs = Vector{Float64}()
    avgCosts = Vector{Float32}() #one per %clutter

    print("\n ------ \n")
    for numPts in numPtsList
        idx, nSuccess, pathcost = 0,0,0.
        pathcosts = Vector{Float32}() #one per %clutter
        stddev_n = 0.

        param = algT.AlgParameters(numPts, connectRadius)

        while idx < N
            idx += 1
            obstacles = Vector{HyperRectangle}()
            obsArea, obstacles = genClutter(perimeter, param, targetNumObs, percentClutter, roomWidth, roomHeight)
            r = algT.Room(roomWidth,roomHeight,walls,obstacles)
            nodeslist, edgeslist = preprocessPRM(r, param)
            pathcost, isPathFound, solPath = queryPRM(startstate, goalstate, nodeslist, edgeslist, obstacles)

            if isPathFound
                nSuccess += 1
            end

            if !(pathcost == Void)
                push!(pathcosts, pathcost)
            end
        end

        totalcost = sum(pathcosts)
        avgCost = totalcost /  nSuccess # average cost across the nTrials
        #@show nSuccess
        pSucc = nSuccess / N

        sumSqdError = sum([cost - avgCost for cost in pathcosts].^2)
        #@show sumSqdError
        #@show pathcosts
        #@show avgCost
        stddev_n = sqrt( sumSqdError / (nSuccess -1))
        #@show stddev_n

        #@printf("For the iter of %d the avg cost was %d across %d trials", maxIter, avgCost, nTrials)
        push!(avgCosts, avgCost) #todo combine pt1
        push!(stddevs, stddev_n) #todo combine pt2
        push!(listpSucc, pSucc)
        print("\n ------ \n")
    end

    timeExperiment = toc();
    print("Time to run experiment: $timeExperiment\n")
    timestamp = Base.Dates.now()

    ####################################
    # Plot Results
    ####################################
    # clutter (area) on X
    # pathcost on Y
    sizeplot = (800,800)

    supTitle="TRAP prm with maxdist=$connectRadius, "*  "Avgd across #samples = $N" *
            "\n(time to run:$(ceil(timeExperiment, 2)) "*
            "\ntimestamp=$timestamp)"
    costTitle= "\nnumPts% vs pathcost\n"

    stddevs = stddevs' #why is this needed for the errors to plot correctly (not constant error)?

    pPRMcost = scatter(numPtsList, avgCosts',
        markercolor = :black,
        title = supTitle * costTitle, ylabel = "euclidean path cost", xlabel = "numPts",
        yaxis=((0,50), 0:5:40), yerr = stddevs) #this automatically centers the 1 stddev, I think. So each side is +- 0.5 stddev?
         # yerr = yerrCost)

    successTitle = "numPts vs pSuccess"
    pPRMsuccess = scatter(numPtsList, listpSucc',
        color = :orange, markersize= 6,
        title = successTitle, title_location=:center, ylabel = ("P(success)=numSucc/$N trials"), xlabel = "numPts",
        yaxis=((0, 1.2), 0:0.1:1))

    plot(pPRMcost, pPRMsuccess, layout=(2,1), legend=false,
         xaxis=((0,320), 0:50:300),
        size = sizeplot)
end

function doorSize_vs_pathCost()
    ####################################
    # Run multiple trials and scatterplot all cost vs area for all runs
    # due to varying absolute area
    ####################################
    ####################################
    ## PARAMETERS - RRT
    ####################################
    connectRadius = 4.
    numPts = 100.
    param = algT.AlgParameters(numPts, connectRadius)
    targetNumObs = 4
    percentClutter = 0.1
    roomWidth,roomHeight  = 20,20

    goalstate = Point(19.,19)

    ####################################
    ## PARAMETERS - plot clutter
    ####################################
    N = 80

    #clutterPercentages = [0 0.01 0.05 0.1 0.12 0.15 0.2 0.25 0.3 0.4 0.5]
    #clutterPercentages = [0 0.1 0.2 0.3 0.5]

    #numPtsList = [10 20 30 80 100]
    doorSizes = [ 0.2 0.5 0.7 1 1.5 2 2.5 3 3.5 4 4.5 4.8]
    #doorSizes = [0.5 1 2 3 4 4.5]
    #doorSizes = [1]

    ####################################
    ## Init
    ####################################
    listCosts = Vector{Float32}()
    listpSucc= Vector{Float32}()

    obstacles = Vector{HyperRectangle}()
    perimeter = HyperRectangle(Vec(0.,0), Vec(roomWidth,roomHeight))
    wallsPerimeter = algfxn.decompRect(perimeter)

    walls = Vector{LineSegment}()

    for l in wallsPerimeter
        push!(walls, l)
    end

    ####################################
    # Run Trials
    ####################################

    plot()
    tic()

    stddevs = Vector{Float64}()
    avgCosts = Vector{Float32}() #one per %clutter

    print("\n ------ \n")
    for doorSize in doorSizes

        obstacles = Vector{HyperRectangle}()
        perimeter = HyperRectangle(Vec(0.,0), Vec(roomWidth,roomHeight))
        wallsPerimeter = algfxn.decompRect(perimeter)

        walls = Vector{LineSegment}()
        startstate,trapwalls = genTrapWalls(doorSize)

        for l in wallsPerimeter
            push!(walls, l)
        end
        for l in trapwalls
            push!(walls, l)
        end

        idx, nSuccess, pathcost = 0,0,0.
        pathcosts = Vector{Float32}() #one per %clutter
        stddev_n = 0.

        param = algT.AlgParameters(numPts, connectRadius)

        while idx < N
            idx += 1
            obstacles = Vector{HyperRectangle}()
            obsArea, obstacles = genClutter(perimeter, param, targetNumObs, percentClutter, roomWidth, roomHeight)
            r = algT.Room(roomWidth,roomHeight,walls,obstacles)
            nodeslist, edgeslist = preprocessPRM(r, param)
            pathcost, isPathFound, solPath = queryPRM(startstate, goalstate, nodeslist, edgeslist, obstacles)

            if isPathFound
                nSuccess += 1
            end

            if !(pathcost == Void)
                push!(pathcosts, pathcost)
            end
        end

        totalcost = sum(pathcosts)
        avgCost = totalcost /  nSuccess # average cost across the nTrials
        #@show nSuccess
        pSucc = nSuccess / N

        sumSqdError = sum([cost - avgCost for cost in pathcosts].^2)
        #@show sumSqdError
        #@show pathcosts
        #@show avgCost
        stddev_n = sqrt( sumSqdError / (nSuccess -1))
        #@show stddev_n

        #@printf("For the iter of %d the avg cost was %d across %d trials", maxIter, avgCost, nTrials)
        push!(avgCosts, avgCost) #todo combine pt1
        push!(stddevs, stddev_n) #todo combine pt2
        push!(listpSucc, pSucc)
        print("\n ------ \n")
    end

    timeExperiment = toc();
    print("Time to run experiment: $timeExperiment\n")
    timestamp = Base.Dates.now()

    ####################################
    # Plot Results
    ####################################
    # clutter (area) on X
    # pathcost on Y
    sizeplot = (600,800)

    supTitle="TRAP prm with numPts=$numPts, maxdist=$connectRadius,
    targetNumObs=$targetNumObs, percentClutter=$percentClutter,
    Avgd across #samples = $N, \n(time to run:$(ceil(timeExperiment, 2)) secs "*
            "timestamp=$timestamp)"

    costTitle= "\ndoorSizes vs pathcost\n"

    stddevs = stddevs' #why is this needed for the errors to plot correctly (not constant error)?

    pPRMcost = scatter(doorSizes, avgCosts',
        markercolor = :black,
        title = supTitle * costTitle, ylabel = "euclidean path cost", xlabel = "doorSizes",
        yaxis=((20,50), 20:5:40), yerr = stddevs) #this automatically centers the 1 stddev, I think. So each side is +- 0.5 stddev?
         # yerr = yerrCost)

    successTitle = "doorSizes vs pSuccess"
    pPRMsuccess = scatter(doorSizes, listpSucc',
        color = :orange, markersize= 6,
        title = successTitle, title_location=:center, ylabel = ("P(success)=numSucc/$N trials"), xlabel = "doorSizes",
        yaxis=((0, 1.2), 0:0.1:1))

    plot(pPRMcost, pPRMsuccess, layout=(2,1), legend=false,
         xaxis=((0,5), 0:0.5:5),
        size = sizeplot)
end

function plotDoorSize()
    ####################################
    ## PARAMETERS - RRT
    ####################################
    connectRadius = 8.
    numPts = 100.
    param = algT.AlgParameters(numPts, connectRadius)
    targetNumObs = 3
    percentClutter = 0.1
    roomWidth,roomHeight  = 20,20

    goalstate = Point(19.,19)
    doorSize = 4
    ####################################
    ## PARAMETERS -INIT
    ####################################

        obstacles = Vector{HyperRectangle}()
        perimeter = HyperRectangle(Vec(0.,0), Vec(roomWidth,roomHeight))
        wallsPerimeter = algfxn.decompRect(perimeter)

        walls = Vector{LineSegment}()
        startstate,trapwalls = genTrapWalls(doorSize)

        for l in wallsPerimeter
            push!(walls, l)
        end
        for l in trapwalls
            push!(walls, l)
        end

        idx, nSuccess, pathcost = 0,0,0.
        pathcosts = Vector{Float32}() #one per %clutter
        stddev_n = 0.

        param = algT.AlgParameters(numPts, connectRadius)
        obstacles = Vector{HyperRectangle}()
        obsArea, obstacles = genClutter(perimeter, param, targetNumObs, percentClutter, roomWidth, roomHeight)
        r = algT.Room(roomWidth,roomHeight,walls,obstacles)
    roomPlot = plotfxn.plotRoom(r)

        nodeslist, edgeslist = preprocessPRM(r, param)
        pathcost, isPathFound, solPath = queryPRM(startstate, goalstate, nodeslist, edgeslist, obstacles)

    ## Plot path found
    timestamp = Base.Dates.now()
    title = "PRM with #doorsize=$doorSize, samples=$numPts, maxDist=$connectRadius,
    pathcost = $(round(pathcost,2)), targetNumObs=$targetNumObs, percentClutter=$percentClutter, \ntimestamp=$timestamp)\n"
    targetNumObs = 3
    percentClutter = 0.1
    roadmap = algT.roadmap(startstate, goalstate, nodeslist, edgeslist)

    #prmPlot= plotfxn.plotPRM(roomPlot, roadmap, solPath, title::String, afont)
    prmPlot= plotfxn.plotPRM(roomPlot, roadmap, solPath, title::String)

    print("\n --- Time --- \n")
    @show timestamp
    #print("\n --- Obstacles --- \n")
    #@show obstacles
    #print("\n --- Nodes --- \n")
    #@show nodeslist
    #print("\n --- Edges --- \n")
    #@show edgeslist
    # print("\n -------- \n")

    #f = Plots.Font(05,"Liberation Serif")
    #plot!(prmPlot,titlefont = f)

    gui(prmPlot)
        print("\n --- Press ENTER to continue --- \n")
        #readline(STDIN)
        # todo: save one "representative" figure from each run
        # print("\n --- continue --- \n")
    @show  isPathFound
end

#walls_vs_nSamples()
doorSize_vs_pathCost()

#plotDoorSize()


text:
module GetDigit_n_xc
    using Printf
    using Roots

    const INITX = BigFloat(5)
    const MAX_PRECISION = 100
    const π_str = "3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798"

    function getxc(func, dh::BigFloat, digit, half_presicion)
        setprecision(MAX_PRECISION)
        if half_presicion
            f = x -> func(x, dh) - BigFloat(π) + BigFloat(5) * BigFloat(10)^(BigFloat(-digit - 1))
        else
            f = x -> func(x, dh) - BigFloat(π) + BigFloat(10)^(BigFloat(-digit - 1))
        end

        return find_zero(f, INITX, Order8())
    end

    function getxc(func, maxn::Int64, digit, half_presicion)
        setprecision(MAX_PRECISION)
        if half_presicion
            f = x -> func(x, maxn) - BigFloat(π) + BigFloat(5) * BigFloat(10)^(BigFloat(-digit - 1))
        else
            f = x -> func(x, maxn) - BigFloat(π) + BigFloat(10)^(BigFloat(-digit - 1))
        end

        return find_zero(f, INITX, Order8())
    end

    function getn_linear(func, xc, maxn, precision, half_precision)
        for n = 1:maxn + 1
            if n == maxn + 1
                return nothing
            end

            setprecision(MAX_PRECISION)
            res = abs(func(xc, n) - BigFloat(π))
            @printf("n = %d, res = %.50f\n", n, res)

            if half_precision && res < BigFloat(5) * BigFloat(10)^(BigFloat(-precision - 1))
                return n
            elseif !half_precision && res < BigFloat(10)^(BigFloat(-precision - 1))
                return n
            end
        end

        return nhi != MAXN ? nhi : nothing
    end

    function getn_even_linear(func, xc, maxn, precision, half_precision)
        maxn = maxn & 1 == 0 ? maxn : maxn + 1
        for n = 2:2:maxn + 2
            if n == maxn + 2
                return nothing
            end

            setprecision(MAX_PRECISION)
            res = abs(func(xc, n) - BigFloat(π))
            @printf("n = %d, res = %.50f\n", n, res)

            if half_precision && res < BigFloat(5) * BigFloat(10)^(BigFloat(-precision - 1))
                return n
            elseif !half_precision && res < BigFloat(10)^(BigFloat(-precision - 1))
                return n
            end
        end

        return nhi != MAXN ? nhi : nothing
    end

    function getdigit(func, xc, n, precision, half_precision)
        plo = precision
        phi = MAX_PRECISION

        # 正しい桁を二分探索で求める
        while phi - plo > 1
            p = (phi + plo) >> 1

            setprecision(p)
            res = abs(func(xc, n) - BigFloat(π))
            if half_precision && res < BigFloat(5) * BigFloat(10)^(BigFloat(-precision - 1))
                phi = p
            elseif !half_precision && res < BigFloat(10)^(BigFloat(-precision - 1))
                phi = p
            else
                plo = p
            end

            @printf("p = %d, plo = %d, phi = %d, res = %.100f\n", p, plo, phi, res)
        end

        return phi != MAX_PRECISION ? phi : nothing
    end

    function getdigit_2(func, xc, n, digit, precision)
        for p = precision:-1:digit
            setprecision(p)

            s = roundoff(func(xc, n), digit)
            @printf("p = %d, res = %s\n", p, s)
            if s != roundoff(π_str, digit)
                setprecision(p + 1)

                return roundoff(func(xc, n), digit) == roundoff(π_str, digit) ? p + 1 : getdigit_3(func, xc, n, digit, p + 1)
            end
        end

        return nothing
    end

    function getdigit_3(func, xc, n, digit, precision)
        p = precision

        while true
            setprecision(p)

            s = roundoff(func(xc, n), digit)
            @printf("p = %d, res = %s\n", p, s)
            if s == roundoff(π_str, digit)
                break
            end

            p += 1
        end

        return p
    end

    function roundoff(value::BigFloat, digit)
        fmt = "%." * string(digit) * "f"

        return @eval @sprintf($fmt, $value)
    end

    function roundoff(value::String, digit)
        res = collect(value)[1:digit + 3]
        if res[digit + 3] >= '5' && res[digit + 2] == '9'
            res[digit + 2] = '0'
            res[digit + 1] = Char(res[digit + 1] + 1)
        elseif res[digit + 3] >= '5'
            res[digit + 2] = Char(res[digit + 2] + 1)
        end

        return join(res[1:digit + 2])
    end
end

ext: jl


text:
"""
    mean_overlap(a, y0, x0, x1, y1)

Mean overlap of an interval of length `a` and interval `[x0,x1]`,
when moving inside `[y0,y1]`.

Assumptions `a <= y1 - y0` and ` y0 <= x0 <= x1 <= y1`.

Return `result` with `min(x1 - x0, a) <= result <= a`.
"""
function mean_overlap(a::T, y0::T, x0::T, x1::T, y1::T) where T <: Real
    x0 = max(x0, y0)
    x1 = min(x1, y1)
    y0 <= x0 <= x1 <= y1 || throw(ArgumentError("violating ordering of boundaries"))
    a = min(a, y1 - y0)
    y10 = y1 - y0 - a
    if y10 <= (eps())*a
        return min(x1 - x0, a)
    end
    m11 = max(x0 - y0 - a, 0)
    m12 = max(min(x0+a, x1) - max(y0+a, x0), 0)
    d12 = (min(x0+a, x1) + max(y0+a, x0)) / 2 - x0
    m13 = max(min(x0+a, y1) - max(y0+a, x1), 0)
    m22 = max(x1 - x0 - a, 0)
    m23 = max(min(x1+a, y1) - max(x0+a, x1), 0)
    d23 = -(min(x1+a, y1) + max(x0+a, x1)) /2 + x1 + a
    m33 = max(y1 - x1 - a, 0)

    #println(m11 + m12 + m13 + m22 + m23 + m33 - y1 + y0 + a)
    #display([m11 m12 m13 m22 m23 m33;0 d12 x1 - x0 a d23 0])
    (m12*d12 + m13*(x1-x0) + m22*a + m23*d23) / y10
end

function sumeff(points, vmin, vmax, poly, pmin, pmax, xy, tol)
    pmi = pmin[xy] - tol
    pma = pmax[xy] + tol
    vmi = vmin[xy] - tol
    vma = vmax[xy] + tol
    sum(mean_overlap.(all_boxlengths(poly, xy) .+ 2tol, pmi, vmi, vma, pma))
end

function countmid(points, pmin, pmax, xy, tol)
    npoints = length(points)
    v = (vertex(points, k, xy) for k = 1:npoints)
    pmi = pmin[xy] - tol
    pma = pmax[xy] + tol
    low, mid, hi = sum([x < pmi, pmi <= x <= pma, x > pma] for x in v)
    mid
end

function estimate_effort(points, vmin, vmax, poly, pmin, pmax, tol)
    sumx = sumeff(points, vmin, vmax, poly, pmin, pmax, 1, tol)
    sumy = sumeff(points, vmin, vmax, poly, pmin, pmax, 2, tol)
    sumx0 = sumeff(points, vmin, vmax, poly, pmin, pmax, 1, 0.0)
    sumy0 = sumeff(points, vmin, vmax, poly, pmin, pmax, 2, 0.0)

    xmid = countmid(points, pmin, pmax, 1, tol)
    ymid = countmid(points, pmin, pmax, 2, tol)
    xmid0 = countmid(points, pmin, pmax, 1, 0.0)
    ymid0 = countmid(points, pmin, pmax, 2, 0.0)

    xmid * sumx, ymid * sumy, xmid0 * sumx0, ymid0 * sumy0
end


text:
@safetestset "PaduaTransforms" begin include("test_PaduaTransforms.jl") end

using PoincareInvariants.SecondChebyshevPlans

@safetestset "Differentiation" begin
    using ..SecondChebyshevPlans: DiffPlan, differentiate!

    @test DiffPlan{Float64}(5).D == [0  1  0  3  0   5;
                                     0  0  4  0  8   0;
                                     0  0  0  6  0  10;
                                     0  0  0  0  8   0;
                                     0  0  0  0  0  10;
                                     0  0  0  0  0   0]

    @test DiffPlan{Float64}(4).D == [0  1  0  3  0;
                                     0  0  4  0  8;
                                     0  0  0  6  0;
                                     0  0  0  0  8;
                                     0  0  0  0  0]

    coeffs = [1  2  3  0;
              4  5  0  0;
              6  0  0  0;
              0  0  0  0]

    P = DiffPlan{Float64}(3)

    ∂x = zeros(4, 4)
    ∂y = zeros(4, 4)

    differentiate!(∂x, ∂y, P, coeffs)

    @test ∂x == [2  12  0  0;
                 5   0  0  0;
                 0   0  0  0;
                 0   0  0  0]

    @test ∂y == [ 4  5  0  0;
                 24  0  0  0;
                  0  0  0  0;
                  0  0  0  0]

    ndcoeffs = ntuple(_ -> rand(4, 4), 6)
    nd∂x = ntuple(_ -> zeros(4, 4), 6)
    nd∂y = ntuple(_ -> zeros(4, 4), 6)

    differentiate!(nd∂x, nd∂y, P, ndcoeffs)

    for i in 1:6
        differentiate!(∂x, ∂y, P, ndcoeffs[i])
        @test nd∂x[i] == ∂x
        @test nd∂y[i] == ∂y
    end
end

@safetestset "Integration" begin
    using ..SecondChebyshevPlans: getintweights, integrate

    # integrating odd polynomials over symmetric boundary conditions gives 0
    @test all(getintweights(100)[2:2:end] .== 0)
    @test all(getintweights(999)[2:2:end] .== 0)

    coeffs = [1 4 7 0;
              2 5 8 0;
              3 6 9 0;
		      0 0 0 0]

    @test integrate(coeffs, getintweights(3)) ≈ 1 * 4 + 3 * -4/3 + 7 * -4/3 + 9 * 4/9 atol=5eps()
end

@safetestset "getintegrand! with CallIntPlan" begin
    using ..SecondChebyshevPlans.PaduaTransforms
    using ..SecondChebyshevPlans: DiffPlan, differentiate!, CallIntPlan, getintegrand!
    using PoincareInvariants

    @testset "6 Points in 2 Dimensions" begin
        degree = 2
        D = 2
        phasecoeffs = (
            Float64[1 2 3;
                    4 5 0;
                    6 0 0],
            Float64[1 4 6;
                    2 5 0;
                    3 0 0]
        )

        ∂xcoeffs = ntuple(_ -> zeros(degree+1, degree+1), D)
        ∂ycoeffs = ntuple(_ -> zeros(degree+1, degree+1), D)
        differentiate!(∂xcoeffs, ∂ycoeffs, DiffPlan{Float64}(degree), phasecoeffs)

        phasepoints = Matrix{Float64}(undef, getpaduanum(degree), D)
        iplan = InvPaduaTransformPlan{Float64}(degree)
        invpaduatransform!(phasepoints, iplan, phasecoeffs)

        ω(v, t, p) = [0 -1; 1  0]

        plan = CallIntPlan{Float64, D}(degree)

        intcoeffs = zeros(degree+1, degree+1)

        getintegrand!(intcoeffs, plan, ω, phasepoints, 0, nothing, ∂xcoeffs, ∂ycoeffs)

        @test intcoeffs ≈ [-72   -82  -30;
                           -82  -432    0;
                           -30     0    0] atol=5eps()
    end

    @testset "50 Points in 6 Dimensions" begin
        D = 6
        N = 50
        pointnum = nextpaduanum(N)
        degree = getdegree(pointnum)

        phasecoeffs = ntuple(_ -> zeros(degree+1, degree+1), D)
        phasecoeffs[1][1, 1] = 0.5  # 0.5
        phasecoeffs[1][1, 2] = 0.5  # 0.5 x
        phasecoeffs[D ÷ 2 + 1][1, 1] = 0.5  # 0.5
        phasecoeffs[D ÷ 2 + 1][2, 1] = 0.5  # 0.5 y

        ∂xcoeffs = ntuple(_ -> zeros(degree+1, degree+1), D)
        ∂ycoeffs = ntuple(_ -> zeros(degree+1, degree+1), D)
        differentiate!(∂xcoeffs, ∂ycoeffs, DiffPlan{Float64}(degree), phasecoeffs)

        phasepoints = Matrix{Float64}(undef, getpaduanum(degree), D)
        iplan = InvPaduaTransformPlan{Float64}(degree)
        invpaduatransform!(phasepoints, iplan, phasecoeffs)

        ω = canonical_two_form

        plan = CallIntPlan{Float64, D}(degree)

        intcoeffs = zeros(degree+1, degree+1)

        getintegrand!(intcoeffs, plan, ω, phasepoints, 0, nothing, ∂xcoeffs, ∂ycoeffs)

        testcoeffs = zeros(degree+1, degree+1)
        testcoeffs[1, 1] = 0.25

        @test testcoeffs ≈ intcoeffs atol=1eps()
    end
end

@safetestset "compute! with SecondChebyshevPlan (OOP)" begin
    using ..SecondChebyshevPlans.PaduaTransforms
    using ..SecondChebyshevPlans: SecondChebyshevPlan, compute!, DiffPlan
    using PoincareInvariants

    D = 12
    N = 200
    ω = canonical_two_form
    pointnum = nextpaduanum(N)
    degree = getdegree(pointnum)
    plan = SecondChebyshevPlan{Float64, D}(ω, pointnum)

    testphasecoeffs = [zeros(degree+1, degree+1) for _ in 1:D]
    test∂x = [zeros(degree+1, degree+1) for _ in 1:D]
    test∂y = [zeros(degree+1, degree+1) for _ in 1:D]
    testintcoeffs = zeros(degree+1, degree+1)

    phasepoints = getpaduapoints(degree) do x, y
        ntuple(D) do i
            i == 1 && return x
            i == D ÷ 2 + 1 && return y
            return 0
        end
    end

    pinv = SecondPoincareInvariant{Float64, D}(ω, pointnum, plan)
    @test compute!(pinv, phasepoints, 0, nothing) ≈ 4 atol=20eps()

    testphasecoeffs[1][1, 2] = 1  # 1 x
    testphasecoeffs[D ÷ 2 + 1][2, 1] = 1 # 1 y

    test∂x[1][1, 1] = 1
    test∂y[D ÷ 2 + 1][1, 1] = 1

    for d in 1:D
        @test maximum(abs, plan.phasecoeffs[d] .- testphasecoeffs[d]) / eps() < 10
        @test maximum(abs, plan.∂x[d] .- test∂x[d]) / eps() < 50
        @test maximum(abs, plan.∂y[d] .- test∂y[d]) / eps() < 50
    end

    testintcoeffs[1, 1] = 1
    @test maximum(abs, plan.intcoeffs .- testintcoeffs) / eps() < 50
end

@safetestset "getpoints, getpointspec and getpointnum" begin
    using ..SecondChebyshevPlans.PaduaTransforms
    using ..SecondChebyshevPlans: SecondChebyshevPlan, getpointspec, getpoints

    for N in [10, 321, 2178], T in [Float32, Float64]
        @test getpointspec(N, SecondChebyshevPlan) == nextpaduanum(N)

        @inferred Vector{T} getpoints((x, y) -> x, T, N, SecondChebyshevPlan)
        @inferred Matrix{T} getpoints((x, y) -> (x, y), T, N, SecondChebyshevPlan)

        f(x, y) = 5 * x, cos(x*y), x+y
        rescale(x, y) = ((x, y) .+ 1) ./ 2

        pnts = getpoints(f, T, N, SecondChebyshevPlan)
        testpnts = getpaduapoints(T, nextdegree(N)) do x, y
            px, py = rescale(x, y)
            return f(px, py)
        end

        @test pnts isa Matrix{T}
        @test size(pnts) == (nextpaduanum(N), 3)
        @test pnts ≈ testpnts
    end
end


text:
gf() = begin
	CSV.read("c:/data/AirOnTimeCSV/airOT198710.csv")
end

gf2() = begin
	p = "c:/data/AirOnTimeCSV/"
	f = joinpath.(p, readdir(p))
	sort!(f, by = x->filesize(x), rev=true)
	reduce(vcat, CSV.read.(f[1:100]))
end

iow() = begin
	open("c:/data/bin.bin", "w")
end

ior() = begin
	open("c:/data/bin.bin", "r")
end


text:
@testset "Minimum assignment problem tests" begin
    begin
        C = [1 4 5; 5 7 6; 5 8 8]
        x =  minimumAssignmentProblem(C)
        @test vecdot(C, x) == 15
    end

    begin
        C = [9 2 7 8; 6 4 3 7; 5 8 1 8; 7 6 9 4]
        x = minimumAssignmentProblem(C)
        @test vecdot(C, x) == 13
    end
end


text:
module PaperFigures
    using   OceanWaveSpectralFitting,
            RecipesBase,
            LaTeXStrings,
            DSP,
            Statistics,
            JLD2,
            Dates
    import OceanWaveSpectralFitting.WhittleLikelihoodInference: HermitianPlot
    include("../freqdirspectra/FrequencyDirectionSpectra.jl")
    include("../exampleanalysis/PaperAnalysis.jl")
    include("../simulation/PaperSimulation.jl")
    include("generic.jl")
    for i in 1:11
        include("Figure$(i).jl")
    end
end


text:
# @inline ldof(S) = Int(2S+1)
function creator(M::Integer)
    ret = Float64[sqrt(m) for m in M:-1:1]
    return diagm(1=>ret)
end
function creator(M::Integer, L::Integer, i::Integer)
    ld = M+1
    leftN = ld^(i-1)
    rightN = ld^(L-i)
    kron(eye(leftN), kron(creator(M), eye(rightN)))
end

function annihilator(M::Integer)
    ret = Float64[sqrt(m+1) for m in M-1:-1:0]
    return diagm(-1=>ret)
end
function annihilator(M::Integer, L::Integer, i::Integer)
    ld = M+1
    leftN = ld^(i-1)
    rightN = ld^(L-i)
    kron(eye(leftN), kron(annihilator(M), eye(rightN)))
end

function number(M::Integer)
    ret = Float64[m for m in M:-1:0]
    return diagm(0=>ret)
end
function number(M::Integer, L::Integer, i::Integer)
    ld = M+1
    leftN = ld^(i-1)
    rightN = ld^(L-i)
    kron(eye(leftN), kron(number(M), eye(rightN)))
end

function numberdensity(M::Integer,L::Integer,
                       k::Integer=0)
    ld = M+1
    N = ld^L
    n = number(M)
    res = zeros(N)
    invV = 1.0/L
    for i in 1:N
        for (j,m) in enumerate(digits(i-1,base=ld,pad=L))
            res[i] += n[m+1,m+1] * cospi(k*invV*(j-1))
        end
    end
    res .*= 1.0/L
    return diagm(0=>res)
end

function bosonchain(M::Integer, L::Integer,
                    t::Real, V::Real, U::Real,
                    mu::Real,
                    Guni::Real, Gstag::Real)
    ld = M+1
    c = creator(M)
    a = annihilator(M)
    x = c+a
    n = number(M)
    i = eye(ld)

    G  = [-(Guni+Gstag), -(Guni-Gstag)]
    onsite = [-mu.*n .+ U.*n.*(n.-1.0) .+ G[i].*x for i in 1:2]
    bondH = V.*kron(n,n) .- t.*(kron(c,a) .+ kron(a,c))

    N = ld^L
    H = zeros(N,N)
    for i in 1:L
        H .+= kron(eye(ld^(i-1)), kron(onsite[mod1(i,2)], eye(ld^(L-i))))
    end
    H .+= V.*kron(n, kron(eye(div(N,ld*ld)), n))
    H .-= t .* kron(c, kron(eye(div(N,ld*ld)), a))
    H .-= t .* kron(a, kron(eye(div(N,ld*ld)), c))

    leftN = 1
    rightN = ld^(L-2)
    for i in 1:L-1
        H .+= kron(eye(leftN), kron(bondH, eye(rightN)))
        leftN *= ld
        rightN = div(rightN,ld)
    end
    return H
end

doc"""
\mathcal{H} = -t\sum_i(a_i c_{i+1} + c_i a_{i+1}) + V\sum_i n_i n_{i+1} + U \sum_i n_i(n_i-1) - mu \sum_i n_i - G \sum_i (c_i + a_i)
"""
struct BosonChainSolver <: Solver
    ef :: Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
    M :: Int
    L :: Int
    function BosonChainSolver(M::Integer, L::Integer
                              ;
                              t::Real=1.0, V::Real=0.0, U::Real=0.0,
                              mu::Real=0.0, Guni::Real=0.0, Gstag::Real=0.0,
                             )
        new(eigen(bosonchain(M, L, t, V, U, mu, Guni, Gstag)), M, L)
    end
end

creator(solver::BosonChainSolver) = creator(solver.M)
creator(solver::BosonChainSolver, i::Integer) = creator(solver.M, solver.L, i)
annihilator(solver::BosonChainSolver) = annihilator(solver.M)
annihilator(solver::BosonChainSolver, i::Integer) = annihilator(solver.M, solver.L, i)
basis(solver::BosonChainSolver) = number(solver.M)
basis(solver::BosonChainSolver, i::Integer) = number(solver.M, solver.L, i)
orderparameter(solver::BosonChainSolver, k::Integer=0) = numberdensity(solver.M, solver.L, k)


text:
# ------------------------------------------------------------------------------------------
# # Álgebra lineal básica en Julia
# Autor: Andreas Noack Jensen (MIT) (http://www.econ.ku.dk/phdstudent/noack/)
# (con edición de Jane Herriman)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Julia's syntax is very similar to other languages but there are some important
# differences. Define a matrix of random normal variates
# ------------------------------------------------------------------------------------------

A = rand(1:4,3,3)

# ------------------------------------------------------------------------------------------
# Definir un vector de unos
# ------------------------------------------------------------------------------------------

x = fill(1.0, (3))

# ------------------------------------------------------------------------------------------
# Notamos que $A$ tiene el tipo Array{Int64,2} pero $x$ tiene tipo Array{Int64,1}. Julia
# define los alias a Vector{Type}=Array{Type,1} y Matrix{Type}=Array{Type,2}.
#
# Muchas de las operaciones básicas son idénticas a otros lenguajes
#
# #### Multiplicación
# ------------------------------------------------------------------------------------------

b = A*x

# ------------------------------------------------------------------------------------------
# #### Traspuestas
# Como en otros lenguajes, `A'` es la transpuesta conjugada mientras que `A.'` es sólo la
# traspuesta
# ------------------------------------------------------------------------------------------

Asym = A + A'

# ------------------------------------------------------------------------------------------
# #### Multiplicación traspuesta
# Julia nos permite escribir esto sin *
# ------------------------------------------------------------------------------------------

Apd = A'A

# ------------------------------------------------------------------------------------------
# #### Resolviendo sistemas lineales
# El problema $Ax=b$ para $A$ cuadrada se resulve con la función \.
# ------------------------------------------------------------------------------------------

A\b

# ------------------------------------------------------------------------------------------
# #### Sistemas sobredeterminados
# Cuando nuestra matriz es alta (número de renglones mayores al número de columnas), tenemos
# un sistema lineal sobredeterminado.
#
#
# En este caso \ calcula la de mínimos cuadrados
# ------------------------------------------------------------------------------------------

Atall = rand(3, 2)
display(Atall)
Atall\b

# ------------------------------------------------------------------------------------------
# La función \ también sirve ocn problemas deficientes de rango de mínimos cuadrados. En
# este caso, la solución no es única y Julia regresa el valor con la menor norma.
#
# Para crear un problema de rango deficiente de mínimos cuadrados, vamos a crear una matriz
# deficiente en rango con columnas linealmente dependientes
# ------------------------------------------------------------------------------------------

v = randn(3)
rankdef = [v v]

rankdef\b

# ------------------------------------------------------------------------------------------
# #### Sistemas indeterminados
# cuando A es corta (número de columnas mayor al número de renglones), tenemos un sistema
# indeterminado
#
# En este caso \ regresa la solución con la norma mínima
# ------------------------------------------------------------------------------------------

Ashort = rand(2, 3)
display(Ashort)
Ashort\b[1:2]

# ------------------------------------------------------------------------------------------
# ### Ejercicios
#
#
#
# ```
# A = [
#  1  2  3  4  5  6  7  8  9  10
#  1  2  3  4  5  6  7  8  9  10
#  1  2  3  4  5  6  7  8  9  10
#  1  2  3  4  5  6  7  8  9  10
#  1  2  3  4  5  6  7  8  9  10
#  1  2  3  4  5  6  7  8  9  10
#  1  2  3  4  5  6  7  8  9  10
#  1  2  3  4  5  6  7  8  9  10
#  1  2  3  4  5  6  7  8  9  10
#  1  2  3  4  5  6  7  8  9  10
#  ]
# ```
#
# Quieres obtener
#
# ```
# A = [
#  7  8  9  10  1  2  3  4  5  6
#  7  8  9  10  1  2  3  4  5  6
#  7  8  9  10  1  2  3  4  5  6
#  7  8  9  10  1  2  3  4  5  6
#  7  8  9  10  1  2  3  4  5  6
#  7  8  9  10  1  2  3  4  5  6
#  7  8  9  10  1  2  3  4  5  6
#  7  8  9  10  1  2  3  4  5  6
#  7  8  9  10  1  2  3  4  5  6
#  7  8  9  10  1  2  3  4  5  6
#  ]
# ```
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 10.2 Toma el producto de un vector `v` con sí mismo.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 10.3 Toma el producto de un vector `v` con sí mismo.
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/intro-to-julia-ES/10. Algebra lineal basica.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using SimpleIntegrals
using SimpleIntegrals: trapezoid
using Test
using InteractiveUtils: subtypes
using BenchmarkTools: @benchmark

using Unitful: Unitful

@testset "perf" begin
    N = 10000
    @info "Benchmarking on array of $N elements"
    for (xs, ys) in [
                     (sort!(randn(N)), randn(N)),
                     (1:N, randn(N)),
                     (sort!(randn(N))*Unitful.cm, randn(N)*Unitful.V),
                    ]
        b = @benchmark integral($xs, $ys) evals=2 samples=2
        @test b.allocs < 50
        display(b)

        window = (first(xs), last(xs))
        b = @benchmark integral($(xs), $ys, window=$window) evals=2 samples=2
        @test b.allocs < 50
        display(b)
    end
end

@testset "Unitful" begin
    xs = [0.0, 1.0] * Unitful.s
    ys = [1, 2] * Unitful.m
    @inferred integral(xs, ys)
    @test integral(xs, ys) ≈ 1.5*Unitful.m * Unitful.s
end

@testset "spike" begin
    xs = [0,1,2]
    ys = [0,1,0]
    @test trapezoid(xs,ys) ≈ 1
    @test trapezoid(xs,ys, 0, 1) ≈ 0.5
    @test trapezoid(xs,ys, 1, 2) ≈ 0.5
    @test trapezoid(xs,ys, 0.1, 1) ≈ 0.5 - 0.1^2/2
    @test trapezoid(xs,ys, 0., 0.9) ≈ 0.9^2/2
    a,b = sort!(rand(2))
    a = 0.1; b =0.9
    @test trapezoid(xs,ys, a,b) ≈ (b-a) * (b+a) /2
end

@testset "constant function" begin
    for _ in 1:100
        n = rand(1:10)
        xs = sort!(randn(n))
        c = randn()
        ys = fill(c, size(xs))
        @test trapezoid(xs, ys) ≈ c * (last(xs) - first(xs))

        a,b = extrema(xs)
        @test trapezoid(xs,ys) ≈ trapezoid(xs,ys,a,b)

        w1,w2 = sort!(rand(2))
        a = minimum(xs) + w1 * (maximum(xs) - minimum(xs))
        b = minimum(xs) + w2 * (maximum(xs) - minimum(xs))
        @test trapezoid(xs,ys,a,b) ≈ c * (b-a)
    end
end

@testset "edgecases" begin
    @test trapezoid(Float64[], Float64[]) === 0.
    @test trapezoid(Float64[1], Float64[1],1,1) === 0.
end

function random_xs_ys(T,n)
    random(T,n) = T.(randn(n))
    xs = sort!(random(T,n))
    ys = random(T,n)
    xs, ys
end

@testset "pasting" begin
    for _ in 1:100
        n1 = rand(1:10)
        n2 = rand(1:10)
        xs1, ys1 = random_xs_ys(Float64,n1)
        xs2, ys2 = random_xs_ys(Float64,n2)
        xs2 .= xs2 .- first(xs2) .+ last(xs1)
        @assert first(xs2) == last(xs1)
        @test trapezoid(xs1,ys1) + trapezoid(xs2, ys2) ≈
            trapezoid([xs1;xs2], [ys1;ys2])
    end
    for _ in 1:100
        n = rand(1:10)
        xs,ys = random_xs_ys(Float64, n)
        w1,w2,w3 = sort!(rand(3))
        a = minimum(xs) + w1 * (maximum(xs) - minimum(xs))
        b = minimum(xs) + w2 * (maximum(xs) - minimum(xs))
        c = minimum(xs) + w3 * (maximum(xs) - minimum(xs))
        @test trapezoid(xs,ys,a,b) + trapezoid(xs,ys,b,c) ≈
            trapezoid(xs,ys,a,c)
    end
end

@testset "typestability" begin
    for T in subtypes(AbstractFloat)
        xs, ys = random_xs_ys(T,0)
        @inferred trapezoid(xs, ys)
        @test T == typeof(trapezoid(xs,ys))

        xs, ys = random_xs_ys(T,2)
        a = b = first(xs)
        @inferred trapezoid(xs, ys, a, b)
        @test T == typeof(trapezoid(xs, ys, a, b))
    end
end

@testset "api" begin
    xs, ys = random_xs_ys(Float64,10)
    @test trapezoid(xs,ys) == integral(xs,ys)
    a,b = extrema(xs)
    @test trapezoid(xs,ys,a,b) == integral(xs,ys, window=(a,b))
end

@testset "Range and Vector consistent" begin
    for _ in 1:100
        xs = range(randn(), length=rand(2:30), step=10rand())
        ys = randn(length(xs))
        @test integral(xs, ys) ≈ integral([xs;],ys)
        w1,w2= sort!(rand(2))
        a = minimum(xs) + w1 * (maximum(xs) - minimum(xs))
        b = minimum(xs) + w2 * (maximum(xs) - minimum(xs))
        @test integral(xs, ys, window=(a,b)) ≈ integral([xs;],ys, window=(a,b))
    end
end


text:
using LinearAlgebraicRepresentation, ViewerGL
Lar = LinearAlgebraicRepresentation
GL = ViewerGL

#filename = "/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/Lar.svg"
#V,EV = Lar.svg2lar(filename)

filename = "/home/lomilaria/space/LinearAlgebraicRepresentation.jl/test/svg/Lar.svg"
V,EV = Lar.svg2lar(filename)

GL.VIEW([
	GL.GLLines(V,EV,GL.COLORS[1]),
	GL.GLFrame2
]);

function pointsinout(V,EV, n=2_000_000)
	point_in = []
	point_out = []
	point_on = []
	classify = Lar.pointInPolygonClassification(V,EV)
	for k=1:n
		queryPoint = rand(1,2)
		inOut = classify(queryPoint)
		# println("k = $k, queryPoint = $queryPoint, inOut = $inOut")
		if inOut=="p_in"
			push!(point_in, queryPoint);
		elseif inOut=="p_out"
			push!(point_out, queryPoint);
		elseif inOut=="p_on"
			push!(point_on, queryPoint);
		end
	end
	#GL.GLPoints(queryPoint,GL.COLORS[2])
	return point_in,point_out,point_on
end

#funzione parallelizzata
function pointsinout_multithreading(V,EV, n=2_000_000,t=Threads.nthreads())

  point_in = [[] for i in 1:t]
  point_out = [[] for i in 1:t]
  point_on = [[] for i in 1:t]

  classify = Lar.pointInPolygonClassification(V,EV)

  Threads.@threads for k=1:n

    #queryPoint = rand(rnglist[Threads.threadid()],1,2)
    queryPoint = rand(1,2)
    inOut = classify(queryPoint)
    # println("k = $k, queryPoint = $queryPoint, inOut = $inOut")
    if inOut=="p_in"
      push!(point_in[Threads.threadid()], queryPoint);
    elseif inOut=="p_out"
      push!(point_out[Threads.threadid()], queryPoint);
    elseif inOut=="p_on"
      push!(point_on[Threads.threadid()], queryPoint);
    end
  end
  #GL.GLPoints(queryPoint,GL.COLORS[2])
  #print(point_in,point_out,point_on)
  return vcat(point_in...),vcat(point_out...),vcat(point_on...)
end

function pointsinout_simd(V,EV, n=2_000_000)
	point_in = []
	point_out = []
	point_on = []
	classify = Lar.pointInPolygonClassification(V,EV)
	@inbounds @simd for k=1:n
		pushQueryPoint(point_in,point_out,point_on,classify)
	end
	#GL.GLPoints(queryPoint,GL.COLORS[2])
	return point_in,point_out,point_on
end

function pointsinout_distributed(V,EV, n=2_000_000)
	point_in =  SharedArray{Array{Any,1}}
	point_out =  SharedArray{Array{Any,1}}
	point_on =  SharedArray{Array{Any,1}}
	classify = Lar.pointInPolygonClassification(V,EV)
	@sync @distributed for k=1:n
		queryPoint = rand(1,2)
		inOut = classify(queryPoint)
		if inOut=="p_in"
			push!(point_in, queryPoint);
		elseif inOut=="p_out"
			push!(point_out, queryPoint);
		elseif inOut=="p_on"
			push!(point_on, queryPoint);
		end
	end
	#GL.GLPoints(queryPoint,GL.COLORS[2])
	return point_in,point_out,point_on
end

function pointsinout_tasks(V,EV, n=2_000_000)
	point_in = []
	point_out = []
	point_on = []
	classify = Lar.pointInPolygonClassification(V,EV)
	@sync for k=1:n
		@async pushQueryPoint(point_in,point_out,point_on,classify)
	end
	#GL.GLPoints(queryPoint,GL.COLORS[2])
	return point_in,point_out,point_on
end

function pushQueryPoint(point_in,point_out,point_on,classify)
	queryPoint = rand(1,2)
	inOut = classify(queryPoint)
	if inOut=="p_in"
		push!(point_in, queryPoint);
	elseif inOut=="p_out"
		push!(point_out, queryPoint);
	elseif inOut=="p_on"
		push!(point_on, queryPoint);
	end
end
points_in, points_out, points_on = pointsinout(V,EV);
pointsin = [vcat(points_in...) zeros(length(points_in),1)]
pointsout = [vcat(points_out...) zeros(length(points_out),1)]

polygon = [GL.GLLines(V,EV,GL.COLORS[1])];
in_mesh = [GL.GLPoints(pointsin, GL.COLORS[2])]
out_mesh = [GL.GLPoints(pointsout, GL.COLORS[3])]

result = cat([polygon,in_mesh,out_mesh])
GL.VIEW(result);


text:
using Plots
using GraphRecipes
using CSV

tickets_df = CSV.read("833829401_T_DB1B_MARKET.csv")

coords_df = CSV.read("833829401_T_MASTER_CORD.csv")

US_airports = coords_df[coords_df.AIRPORT_COUNTRY_NAME .== "United States",:]

y = Array{Float64}(collect(US_airports[:LATITUDE])[1:(end-1)])

x = Array{Float64}(collect(US_airports[:LONGITUDE])[1:(end-1)])

tmp = findall(x.>0)
x[tmp] = x[tmp] .- 360

US_tickets = tickets_df[(tickets_df.ORIGIN_COUNTRY .== "US") .& (tickets_df.DEST_COUNTRY .== "US"),:]

US_airports_ids = US_airports[:AIRPORT_ID][1:(end-1)]
US_airports_ids_mapping = Dict(US_airports_ids[i]=>i for i in 1:length(US_airports_ids))
US_tickets_src = US_tickets[:ORIGIN_AIRPORT_ID]
US_tickets_dst = US_tickets[:DEST_AIRPORT_ID]

src,dst = Int64[],Int64[]
for i in 1:length(US_tickets_src)
    push!(src,US_airports_ids_mapping[US_tickets_src[i]])
    push!(dst,US_airports_ids_mapping[US_tickets_dst[i]])
    # convert each edge to undirected
    push!(dst,US_airports_ids_mapping[US_tickets_src[i]])
    push!(src,US_airports_ids_mapping[US_tickets_dst[i]])
end

using SparseArrays
function writeSMAT(filename::AbstractString, A::SparseMatrixCSC{T,Int}; values::Bool=true) where T
    open(filename, "w") do outfile
        write(outfile, join((size(A,1), size(A,2), nnz(A)), " "), "\n")

        rows = rowvals(A)
        vals = nonzeros(A)
        m, n = size(A)
        for j = 1:n
           for nzi in nzrange(A, j)
              row = rows[nzi]
              val = vals[nzi]
              if values
                write(outfile, join((row-1, j-1, val), " "), "\n")
              else
                write(outfile, join((row-1, j-1, 1), " "), "\n")
              end
           end
        end
    end
end

A = sparse(src,dst,ones(length(src)),length(x),length(x))
writeSMAT("flight-tickets.smat",A)

xy = hcat(x,y)

using DelimitedFiles
writedlm("us-airports.xy", xy)
writedlm("us-airports.labels", US_airports[:DISPLAY_AIRPORT_NAME][1:(end-1)])
writedlm("us-airports-metadata.csv",
  zip(US_airports[:AIRPORT_ID][1:(end-1)], US_airports[:AIRPORT_STATE_NAME][1:(end-1)], US_airports[:CITY_MARKET_ID][1:(end-1)]))

max_stars_repo_name: empet/more-graphs

max_issues_repo_name: empet/more-graphs

max_forks_repo_name: empet/more-graphs


text:
module PosteriorPlots

using DataFrames, Statistics, Plots, StatsPlots, StatsBase, Random, KernelDensity, StatsBase, MCMCChains, SampleChains

include("PlotHelpers.jl")
include("Inference.jl")
include("PosteriorPredictiveCheck.jl")

# Exports

export plot_posterior_intervals
export plot_posterior_hist
export plot_posterior_density
export plot_posterior_check
export plot_hist_check
export plot_ecdf_check

end


text:
using SafeTestsets
@safetestset "my_f_tests" begin include("my_f_tests.jl") end


text:
include("1.jl")
function stepgood(r::Robot)
    putmarker!(r)
    i = 1
    while isborder(r,Nord)
        backing(r,HorizonSide(mod(i, 4)))
        putmarker!(r)
        i = mod(i+2,4)
    end
end

max_stars_repo_name: bubasooo/proggggggggg

max_issues_repo_name: bubasooo/proggggggggg

max_forks_repo_name: bubasooo/proggggggggg


text:
const Qobj = Dict{String,Any}

"""
    TargetQobjQASM

Qobj (of type QASM) target. See https://arxiv.org/pdf/1809.03452.pdf
for specifications.
"""
Base.@kwdef struct TargetQobjQASM
    nshots::Int = 1024
    seed::Int = 1
    max_credits::Int = 3
end

# NOTE: Qobj can only generated after running at least :julia optimization pass
#       this is because we need to inline all quantum routines
#       result program can only have constant intrinsic routine

mutable struct CodeGenQobjState
    nqubits::Int

    pc::Int
    stmt::Any

    # NOTE:
    # we alloate memory_slots for all measurement
    # at level of codegen since measure cannot be
    # easily removed it can affect the whole state
    memory_offset::Int
    memory_slots::Int
    # alloate for used measure result
    register_map::Dict{Int,Vector{Int}}
    # allocate qubits on contiguous address
    qubits_map::Dict{Int,Int}
end

function allocate_qobj_qubits!(qubits, locs::Locations, curr::Ref{Int})
    for k in locs.storage
        if !haskey(qubits, k)
            qubits[k] = curr[]
            curr[] += 1
        end
    end
    return qubits
end

function CodeGenQobjState(ci::CodeInfo)
    memory_slots = 0
    register_head = 0
    register_map = Dict{Int,Vector{Int}}()
    qubits_map = Dict{Int,Int}()
    measure_used = measure_ssa_uses!(Set{Int64}(), ci)
    qmem_ptr = Ref{Int}(0)

    for (v, stmt) in eachstmt(ci)
        is_quantum_statement(stmt) || continue
        qt = quantum_stmt_type(stmt)

        if qt === :measure
            _, locs = obtain_const_measure_stmt(stmt, ci)
            allocate_qobj_qubits!(qubits_map, locs, qmem_ptr)
            memory_slots += length(locs)

            if v in measure_used
                register_map[v] = collect(register_head:register_head+length(locs)-1)
                register_head += length(locs)
            end
        elseif qt === :gate
            _, _, locs = obtain_const_gate_stmt(stmt, ci)
            allocate_qobj_qubits!(qubits_map, locs, qmem_ptr)
        elseif qt === :ctrl
            _, _, locs, ctrl = obtain_const_ctrl_stmt(stmt, ci)
            allocate_qobj_qubits!(qubits_map, locs, qmem_ptr)
            allocate_qobj_qubits!(qubits_map, ctrl.storage, qmem_ptr)
        end
    end

    return CodeGenQobjState(qmem_ptr[], 0, nothing, 0, memory_slots, register_map, qubits_map)
end

function get_qubit_addrs(ctrl::CtrlLocations, st::CodeGenQobjState)
    return get_qubit_addrs(ctrl.storage, st)
end

function get_qubit_addrs(storage, st::CodeGenQobjState)
    qubits = Vector{Int}(undef, length(storage))
    for (k, l) in enumerate(storage)
        qubits[k] = st.qubits_map[l]
    end
    return qubits
end

function get_qubit_addrs(locs::Locations, st::CodeGenQobjState)
    return get_qubit_addrs(locs.storage, st)
end

"scan used measurement result"
function measure_ssa_uses!(used, ci::CodeInfo)
    measure_ssa = Int[]
    for (v, stmt) in eachstmt(ci)
        is_quantum_statement(stmt) || continue
        qt = quantum_stmt_type(stmt)

        if qt === :measure
            push!(measure_ssa, v)
        end
    end

    for (v, stmt) in eachstmt(ci)
        stmt isa GotoIfNot || continue
        cond = stmt.cond
        cond_stmt = ci.code[stmt.cond.id]
        for useref in userefs(cond_stmt)
            val = useref[]
            if isa(val, SSAValue) && val.id in measure_ssa
                push!(used, val.id)
            end
        end
    end
    return used
end

# NOTE: we only generate experiment Qobj
# the complete Qobj should be handled in
# client before sending to IBM

function codegen(target::TargetQobjQASM, ci::CodeInfo)
    st = CodeGenQobjState(ci)
    spec = ci.parent.specTypes.parameters[2]
    prog = Qobj(
        "header" => Qobj(
            "name" => routine_name(spec),
            "n_qubits" => st.nqubits,
            "memory_slots" => st.memory_slots,
        ),
        "config" => codegen_config(target, ci, st),
        "instructions" => codegen_inst(target, ci, st),
    )
end

function codegen_config(target::TargetQobjQASM, ci::CodeInfo, st::CodeGenQobjState)
    return Qobj(
        "n_qubits" => st.nqubits,
        "shots" => target.nshots,
        "memory_slots" => st.memory_slots, # TODO: generate this from ci
        "seed" => target.seed,
        "max_credits" => target.max_credits,
    )
end

function codegen_inst(target::TargetQobjQASM, ci::CodeInfo, st::CodeGenQobjState)
    prog = Any[]
    st.pc = 1
    while st.pc <= length(ci.code)
        inst = nothing
        st.stmt = ci.code[st.pc]
        if st.stmt isa Expr && is_quantum_statement(st.stmt)
            inst = codegen_stmt(target, ci, st)
            st.pc += 1
        elseif st.stmt isa GotoIfNot
            inst = codegen_ifnot(target, ci, st)
        else
            st.pc += 1
        end

        if !isnothing(inst)
            if inst isa Vector
                append!(prog, inst)
            else
                push!(prog, inst)
            end
        end
    end
    return prog
end

function codegen_stmt(target::TargetQobjQASM, ci::CodeInfo, st::CodeGenQobjState)
    qt = quantum_stmt_type(st.stmt)
    if qt === :gate
        inst = codegen_gate(target, ci, st)
    elseif qt === :ctrl
        inst = codegen_ctrl(target, ci, st)
    elseif qt === :measure
        inst = codegen_measure(target, ci, st)
    elseif qt === :barrier
        inst = codegen_barrier(target, ci, st)
    else
        error("incompatible statement for Qobj: $(st.stmt)")
    end
    return inst
end

function codegen_gate(::TargetQobjQASM, ci::CodeInfo, st::CodeGenQobjState)
    gate, gt, locs = obtain_const_gate_stmt(st.stmt, ci)

    gt <: RoutineSpec &&
        error("non-intrinsic routine used: $gt, turn on the optimizer or revise your program")
    gt <: IntrinsicRoutine || error("invalid gate type: $gate::$gt")
    gate isa IntrinsicRoutine || error("gate $gate is not constant")
    # TODO: check if the intrinsic is available in target hardware
    #       might need a new field in the target to indicate hardware
    #       basis gate.

    qobj = @match gate begin
        ::Intrinsics.XGate => Qobj("name" => "x")

        ::Intrinsics.YGate => Qobj("name" => "y")

        ::Intrinsics.ZGate => Qobj("name" => "z")

        ::Intrinsics.HGate => Qobj("name" => "h")

        ::Intrinsics.SGate => Qobj("name" => "s")

        ::Intrinsics.TGate => Qobj("name" => "t")

        # TODO: sdg, tdg

        Intrinsics.shift(θ) => Qobj("name" => "u1", "params" => Any[θ])

        Intrinsics.Rx(θ) => Qobj("name" => "rx", "params" => Any[θ])

        Intrinsics.Ry(θ) => Qobj("name" => "ry", "params" => Any[θ])

        Intrinsics.Rz(θ) => Qobj("name" => "rz", "params" => Any[θ])

        _ => error("invalid gate for Qobj backend: $gate")
    end

    # NOTE: intrinsic instructions only supports single qubit
    # operation, we need to expand the program for our syntax
    # sugar here, given we don't have type information in macros
    if length(locs) == 1
        qobj["qubits"] = get_qubit_addrs(locs, st)
        return Any[qobj]
    else
        qobjs = Any[]
        for l in locs.storage
            obj = copy(qobj)
            obj["qubits"] = get_qubit_addrs(l, st)
            push!(qobjs, obj)
        end
        return qobjs
    end
end

function codegen_ctrl(::TargetQobjQASM, ci::CodeInfo, st::CodeGenQobjState)
    gate, gt, locs, ctrl = obtain_const_ctrl_stmt(st.stmt, ci)
    gt <: RoutineSpec &&
        error("non-intrinsic routine used: $gt, turn on the optimizer or revise your program")
    gt <: IntrinsicRoutine || error("invalid gate type: $gate::$gt")

    all(ctrl.flags) || error("inverse ctrl is not supported in Qobj backend")
    gate isa IntrinsicRoutine || error("gate $gate is not constant")

    # control gates support in qelib1.inc
    qobj = @match (gate, length(ctrl), length(locs)) begin
        (::Intrinsics.XGate, 1, 1) => Qobj("name" => "cx")

        (::Intrinsics.YGate, 1, 1) => Qobj("name" => "cy")

        (::Intrinsics.ZGate, 1, 1) => Qobj("name" => "cz")

        (::Intrinsics.HGate, 1, 1) => Qobj("name" => "ch")

        (::Intrinsics.XGate, 2, 1) => Qobj("name" => "ccx")

        (Intrinsics.Rz(θ), 1, 1) => Qobj("name" => "crz", "params" => Any[θ])

        (Intrinsics.shift(θ), 1, 1) => Qobj("name" => "cu1", "params" => Any[θ])

        (Intrinsics.UGate(α, β, γ), 1, 1) => Qobj("name" => "cu3", "params" => Any[α, β, γ])

        _ => error("invalid control statement for Qobj backend, got: $(st.stmt)")
    end

    qobj["qubits"] = vcat(get_qubit_addrs(ctrl, st), get_qubit_addrs(locs, st))
    return qobj
end

function codegen_measure(::TargetQobjQASM, ci::CodeInfo, st::CodeGenQobjState)
    _, locs = obtain_const_measure_stmt(st.stmt, ci)

    qobj = Qobj(
        "name" => "measure",
        "qubits" => get_qubit_addrs(locs, st),
        "memory" => collect(st.memory_offset:st.memory_offset+length(locs)-1),
    )

    if haskey(st.register_map, st.pc)
        qobj["register"] = st.register_map[st.pc]
    end

    return qobj
end

function codegen_barrier(::TargetQobjQASM, ci::CodeInfo, st::CodeGenQobjState)
    locs = obtain_const(st.stmt.args[2], ci)::Locations

    return Qobj("name" => "barrier", "qubits" => get_qubit_addrs(locs, st))
end

function codegen_ifnot(::TargetQobjQASM, ci::CodeInfo, st::CodeGenQobjState)
    cond = st.stmt.cond::SSAValue
    creg, right = obtain_qasm_ifnot_cond(cond, ci)

    if creg isa SSAValue
        haskey(st.register_map, creg.id) ||
            error("classical post processing measurement result in Qobj is not supported")
        condition = st.register_map[creg.id]
    else
        error("invalid value for measurement result: $creg, possible fix: run :julia optimization pass before codegen")
    end

    length(condition) == 1 || error("Qobj only allows single bit condition for if statement")

    # find the first quantum stmt after goto
    move_to_first_quantum_stmt(ci, st)

    # check if this is actually a reset statement:
    #
    # c = @measure loc
    # if c == 1
    #   loc => X
    # end

    if right == 1
        qt = quantum_stmt_type(st.stmt)
        if qt === :gate
            gate, gt, locs = obtain_const_gate_stmt(st.stmt, ci)
            if gate === Intrinsics.X
                return Qobj("name" => "reset", "qubits" => get_qubit_addrs(locs, st))
            end
        end
    end

    qobj = codegen_stmt(target, ci, st)
    if qobj isa Vector # return from codegen_gate
        for each in qobj
            each["conditional"] = condition[]
        end
    else
        qobj["conditional"] = condition[]
    end
    return qobj
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    EmpiricalVariogram(partition, var₁, var₂=var₁; [parameters])

Compute the empirical (cross-)variogram of the spatial `partition` for
variables `var₁` and `var₂` as described in Hoffimann & Zadrozny 2019.

## References

* Hoffimann, J and Zadrozny, B. 2019. [Efficient variography with partition variograms]
  (https://www.sciencedirect.com/science/article/pii/S0098300419302936)
"""
function EmpiricalVariogram(partition::Partition, var₁::Symbol, var₂::Symbol=var₁; kwargs...)
  # retain spatial data with at least 2 points
  filtered = Iterators.filter(d -> nelements(d) > 1, partition)

  @assert !isempty(filtered) "invalid partition of spatial data"

  γ(d) = EmpiricalVariogram(d, var₁, var₂; kwargs...)
  foldxt(merge, Map(γ), collect(filtered))
end

"""
    DirectionalVariogram(direction, data, var₁, var₂=var₁; dtol=1e-6, [parameters])

Computes the empirical (cross-)variogram for the variables `var₁` and `var₂` stored in
geospatial `data` along a given `direction` with band tolerance `dtol`.

Optional parameters include the parameters for [`EmpiricalVariogram`](@ref) and the
parameters for [`DirectionPartition`](@ref).
"""
function DirectionalVariogram(dir, data, var₁, var₂=var₁; dtol=1e-6, kwargs...)
  p = partition(data, DirectionPartition(dir; tol=dtol))
  EmpiricalVariogram(p, var₁, var₂; kwargs...)
end

"""
    PlanarVariogram(normal, data, var₁, var₂=var₁; ntol=1e-6, [parameters])

Computes the empirical (cross-)variogram for the variables `var₁` and `var₂` stored in
geospatial `data` along a plane perpendicular to a `normal` direction with plane
tolerance `ntol`.

Optional parameters include the parameters for [`EmpiricalVariogram`](@ref) and the
parameters for [`PlanePartition`](@ref).
"""
function PlanarVariogram(normal, data, var₁, var₂=var₁; ntol=1e-6, kwargs...)
  p = partition(data, PlanePartition(normal; tol=ntol))
  EmpiricalVariogram(p, var₁, var₂; kwargs...)
end


text:
module ValueFunctionSolver
    export gss, it_solve, index_vcattensor, tensor_grid

    const invphi = (sqrt(5) - 1) / 2
    const invphi2 = (3 - sqrt(5)) / 2
    """
    Golden section search.

    Given a function f with a single local minumum in
    the interval [a,b], gss returns a subset interval
    [c,d] that contains the minimum with d-c <= tol.

    example:
    f(x) = (x-2)^2
    a = 1
    b = 5
    tol = 1e-5
    bracket = gss(f, a, b, tol = tol)
    result:
    (1.9999959837979107, 2.0000050911830893)
    """
    function gss(f,a::Float64,b::Float64; tol::Float64=1e-5)
        (a,b)=(min(a,b),max(a,b))
        h = b - a
        if h <= tol
            return (a,b)
        end
        # required steps to achieve tolerance
        n = Int(ceil(log(tol/h)/log(invphi)))
        c = a + invphi2 * h
        d = a + invphi * h
        yc = f(c)
        yd = f(d)

        for k in 1:n
            if yc < yd
                b = d
                d = c
                yd = yc
                h = invphi*h
                c = a + invphi2 * h
                yc = f(c)
            else
                a = c
                c = d
                yc = yd
                h = invphi*h
                d = a + invphi * h
                yd = f(d)
            end
        end
        if yc < yd
            return (a,d)
        else
            return (c,b)
        end
    end

    function it_solve(Vguess, update_V!; tol = 1e-5, maxit = 1000, verbose = true)
        V = copy(Vguess)
        dist = tol*1e2
        for it = 1:maxit
            if verbose
                println("Iteration $it")
            end
            update_V!(V)
            dist = maximum(abs, V-Vguess)
            if verbose
                println("Distance $dist")
            end
            if isnan(dist)
                throw("Distance NaN")
            end
            if dist < tol
                return V, dist
            end
            Vguess[:] = V[:]
        end
        println("Did not reach specified tolerance within maximum number of iterations. Current dist: $dist")
        return V, dist
    end

    """
    index_vcattensor(iDim::Array{Int}, Idim::Array{Int})

    Row Index into vcat of a matrix. Subtensor position in each dimension. First dimension moves first.

    Example:
    I1 = 2
    I2 = 3
    X = [i1*i2 for i1=1:I1, i2=1:I2]

    2×3 Array{Int64,2}:
    1  2  3
    2  4  6

    Xlong = X[:]

    6-element Array{Int64,1}:
    1
    2
    2
    4
    3
    6

    Then:
    X[i1, i2, i3] == Xlong[index_vcattensor([i1, i2, i3], [I1, I2, I3])]
    """
    function index_vcattensor(idim::Array{Int}, Idim::Array{Int})
       myind = 0
       # iterate up to second to last index
       for i in length(idim):-1:2
           myind += (idim[i] - 1)*prod(Idim[1:i-1])
       end
       myind += idim[1]
       return myind
   end

   """
   tensor_grid(xs)

   Vcat of tensor grid. First dimension moves first.
   where xs is a collection of unidimensional grids
   """
   function tensor_grid(xs)
       N = [length(xs[i]) for i = 1:length(xs)]
       totN = prod(N)
       d = length(N)
       xmat = Array{eltype(xs[1])}(totN,d)
       for i in 1:d
           repout = Int(totN/prod(N[1:i]))
           repin = Int(totN/N[i]/repout)
           xmat[:, i] = repeat(xs[i], inner=repin, outer=repout)
       end
       return xmat
    end

    """
    tensor_product(xs)

    product along row of Vcat of tensor grid. First dimension moves first.
    where xs is a collection of unidimensional grids
    """
   function tensor_product(xs)
       N = [length(xs[i]) for i = 1:length(xs)]
       totN = prod(N)
       d = length(N)
       xmat = ones(eltype(xs[1]),totN)
       for i in 1:d
           repout = Int(totN/prod(N[1:i]))
           repin = Int(totN/N[i]/repout)
           xmat .= xmat.*repeat(xs[i], inner=repin, outer=repout)
       end
       return xmat
    end

    """
    tensor_grid(xs, ws)
    """
    function tensor_grid(xs,ws)
        N = [length(xs[i]) for i = 1:length(xs)]
        Nw = [length(ws[i]) for i = 1:length(xs)]
        @assert all(N.==Nw)

        totN = prod(N)
        d = length(N)
        xmat = Array{eltype(xs[1])}(totN,d)
        w = ones(eltype(ws[1]),totN)
        for i in 1:d
            repout = Int(totN/prod(N[1:i]))
            repin = Int(totN/N[i]/repout)
            xmat[:, i] = repeat(xs[i], inner=repin, outer=repout)
            w .= w.*repeat(ws[i], inner=repin, outer=repout)
        end
        return xmat, w
     end

end


text:
#defines variables found in DIMEN.INC for use in Julia

const IBX = 200
const IIX = 400
const NBX = 40
const NJX = 30
const NGX = 8
const NCX = 1200
const NPX = 400
const NNX = 250
const NAJX = 40
const NPNTX = 201
const NTIMX = 201
const NSENX = 36
const NRITZX = 800
const NFLPX = 20
const NENGX = 12
const NXGUSX = 10
const NYGUSX = 10
const NSETX = 3
const KPX = 21
const JBX = 102
const NGVX = 12*NJX + NNX
const NGPX = 27 + 12*NSENX
const NGFX  = NFLPX + NENGX
const NGUSX = NXGUSX*NYGUSX
const NFRPX = NFLPX + NENGX +   NXGUSX*(2*NYGUSX-1)
const NFRLX = NFLPX + NENGX + 2*NXGUSX*(2*NYGUSX-1)
const NGVQX = NNX*(IIX-NGX-2*NJX-NENGX) + 27*NJX
const NGPQX = IIX*18
const NGFQX = IIX*NGFX
const NGLX  = NGVX + NGPX + NGFX + 1
const NBRKX = NPX + NJX + NGX + NCX


text:
using Flux
using Flux: @functor
using Zygote
using MAT
#using Plots
using Statistics
using Distributions
using LinearAlgebra
using Random
using BSON
using JLD
using PyPlot
using CUDA
#using Base.Iterators: partition

Xtrain = JLD.load("Data/XYtrain_XYvalid_processed_T1w.jld","Xtrain")
Ytrain = JLD.load("Data/XYtrain_XYvalid_processed_T1w.jld","Ytrain")
XValid = JLD.load("Data/XYtrain_XYvalid_processed_T1w.jld","XValid")
YValid = JLD.load("Data/XYtrain_XYvalid_processed_T1w.jld","YValid")

idx_shuffle = shuffle(collect(1:size(Xtrain,4)))
Xtrain = Xtrain[:,:,1,idx_shuffle]
Ytrain = Ytrain[:,:,:,idx_shuffle]

Xtrain = reshape(Xtrain,size(Xtrain,1),size(Xtrain,1),1,size(Xtrain,3))

Xtrain = gpu(Xtrain)
Ytrain = gpu(Ytrain)

batch_size = 10

train_loader = Flux.Data.DataLoader(Xtrain,Ytrain, batchsize=batch_size,shuffle=true)

XValid = XValid[:,:,1,:]
XValid = reshape(XValid,size(XValid,1),size(XValid,1),1,size(XValid,3))

XValid = gpu(XValid)
YValid = gpu(YValid)

channels_in=1
channels_out=1

# Unet define
function _random_normal(shape...)
  return Float32.(rand(Normal(0.f0,0.02f0),shape...))
end

expand_dims(x,n::Int) = reshape(x,ones(Int64,n)...,size(x)...)
function squeeze(x)
    if size(x)[end] != 1
        return dropdims(x, dims = tuple(findall(size(x) .== 1)...))
    else
        # For the case BATCH_SIZE = 1
        int_val = dropdims(x, dims = tuple(findall(size(x) .== 1)...))
        return reshape(int_val,size(int_val)...,1)
    end
end

BatchNormWrap(out_ch) = BatchNorm(out_ch)

UNetConvBlock(in_chs, out_chs, kernel = (3, 3)) =
    Chain(Conv(kernel, in_chs=>out_chs,pad = (1, 1);init=_random_normal),
	BatchNormWrap(out_chs),
	x->leakyrelu.(x,0.2f0))

ConvDown(in_chs,out_chs,kernel = (4,4)) =
  Chain(Conv(kernel,in_chs=>out_chs,pad=(1,1),stride=(2,2);init=_random_normal),
	BatchNormWrap(out_chs),
	x->leakyrelu.(x,0.2f0))

struct UNetUpBlock
  upsample
end

@functor UNetUpBlock

UNetUpBlock(in_chs::Int, out_chs::Int; kernel = (3, 3), p = 0.5f0) =
    UNetUpBlock(Chain(x->leakyrelu.(x,0.2f0),
       		ConvTranspose((2, 2), in_chs=>out_chs,
			stride=(2, 2);init=_random_normal),
		BatchNormWrap(out_chs),
		Dropout(p)))

function (u::UNetUpBlock)(x, bridge)
  x = u.upsample(x)
  return cat(x, bridge, dims = 3)
end

struct Unet
  conv_down_blocks
  conv_blocks
  conv_blocks_2nd_down
  up_blocks
  conv_blocks_1st_up
  conv_blocks_2nd_up
end

@functor Unet

function Unet(channels::Int = 1)
  conv_down_blocks = Chain(ConvDown(64,64),
		      ConvDown(128,128),
		      ConvDown(256,256),
		      ConvDown(512,512))

  conv_blocks = Chain(UNetConvBlock(channels, 64),
		 UNetConvBlock(64, 128),
		 UNetConvBlock(128, 256),
		 UNetConvBlock(256, 512),
		 UNetConvBlock(512, 1024),
		 UNetConvBlock(1024, 1024))

  # add by Xiaodong Ma: 2nd conv layer in the first half
  conv_blocks_2nd_down = Chain(UNetConvBlock(64, 64),
		 UNetConvBlock(128, 128),
		 UNetConvBlock(256, 256),
		 UNetConvBlock(512, 512))

  up_blocks = Chain(UNetUpBlock(1024, 512),
		UNetUpBlock(512, 256),
		UNetUpBlock(256, 128),
		UNetUpBlock(128, 64,p = 0.0f0),
		Chain(x->leakyrelu.(x,0.2f0),
    Conv((1, 1), 64=>channels;init=_random_normal)))

  # add by Xiaodong Ma: 1st conv layer in the second half
  conv_blocks_1st_up = Chain(UNetConvBlock(1024, 512),
		 UNetConvBlock(512, 256),
		 UNetConvBlock(256, 128),
     UNetConvBlock(128, 64))

  # add by Xiaodong Ma: 1st conv layer in the second half
  conv_blocks_2nd_up = Chain(UNetConvBlock(512, 512),
		 UNetConvBlock(256, 256),
		 UNetConvBlock(128, 128),
		 UNetConvBlock(64, 64))

  Unet(conv_down_blocks,
       conv_blocks,
       conv_blocks_2nd_down,
       up_blocks,
       conv_blocks_1st_up,
       conv_blocks_2nd_up)
end

function (u::Unet)(x::AbstractArray)
  op = u.conv_blocks_2nd_down[1](u.conv_blocks[1](x)) #out_chs:64

  x1 = u.conv_blocks_2nd_down[2](u.conv_blocks[2](u.conv_down_blocks[1](op))) #out_chs:128
  x2 = u.conv_blocks_2nd_down[3](u.conv_blocks[3](u.conv_down_blocks[2](x1))) #out_chs:256
  x3 = u.conv_blocks_2nd_down[4](u.conv_blocks[4](u.conv_down_blocks[3](x2))) #out_chs:512
  x4 = u.conv_blocks[5](u.conv_down_blocks[4](x3)) #out_chs:1024

  up_x4 = u.conv_blocks[6](x4) #out_chs:1024

  up_x1 = u.conv_blocks_2nd_up[1](u.conv_blocks_1st_up[1](u.up_blocks[1](up_x4, x3)))  #out_chs:512
  up_x2 = u.conv_blocks_2nd_up[2](u.conv_blocks_1st_up[2](u.up_blocks[2](up_x1, x2))) #out_chs:256
  up_x3 = u.conv_blocks_2nd_up[3](u.conv_blocks_1st_up[3](u.up_blocks[3](up_x2, x1))) #out_chs:128
  up_x5 = u.conv_blocks_2nd_up[4](u.conv_blocks_1st_up[4](u.up_blocks[4](up_x3, op))) #out_chs:64
  tanh.(u.up_blocks[end](up_x5)) #out_chs:channels
end

u = Unet(channels_in)
u = gpu(u)

#loss(x,y) = mean(abs,(u(x).-y))
loss(x,y) = mean((u(x).-y).^2)
validationloss() = mean(Flux.mse(u(XValid[:,:,:,(idx-1)*100+1:idx*100]),YValid[:,:,:,(idx-1)*100+1:idx*100]) for idx=1:round(Int,size(XValid,4)/100))
#validationloss() = mean(mean(abs,(u(XValid[:,:,:,(idx-1)*100+1:idx*100]).-YValid[:,:,:,(idx-1)*100+1:idx*100])) for idx=1:round(Int,size(XValid,4)/100))

opt = ADAM(0.001)

function my_custom_train!(loss, ps, data, opt)
    # training_loss is declared local so it will be available for logging outside the gradient calculation.
    local training_loss
    ps = Params(ps)
    for d in data
      if d isa AbstractArray{<:Number}
        gs = Flux.gradient(ps) do
          training_loss = loss(d)
          return training_loss
        end
      else
        gs = Flux.gradient(ps) do
          training_loss = loss(d...)
          return training_loss
        end
      end
      # Insert whatever code you want here that needs training_loss, e.g. logging.
      # logging_callback(training_loss)
      # Insert what ever code you want here that needs gradient.
      # E.g. logging with TensorBoardLogger.jl as histogram so you can see if it is becoming huge.
      #println("training_loss = ",training_loss,"; gradient[1] = ",gs[ps[1]],"; W1 = ",ps.order.data[1][1],"; b1 = ",ps.order.data[2][1])
      println("training_loss = ",training_loss)
      Flux.update!(opt, ps, gs)
      # Here you might like to check validation set accuracy, and break out to do early stopping.
    end
  end

loss_valid = zeros(10)
for epoch in 1:10
  #@time Flux.train!(loss, Flux.params(u), train_loader, opt, cb = cb)
  @time my_custom_train!(loss, Flux.params(u), train_loader, opt)
  #@time Flux.train!(loss, Flux.params(u), train_set, opt)
  #model = cpu(u)
  #BSON.@save string("Results/model_unet_ptx_processed_7dir_gpu_epoch",string(epoch),".bson") model
  #u = gpu(model)
  loss_valid[epoch] = validationloss()
end
model = cpu(u)
BSON.@save string("Results/model_unet_ptx_processed_7dir_gpu_orig_epoch10.bson") model

#test
#=
XValid = JLD.load("XYtrain_XYvalid_processed.jld","XValid")
YValid = JLD.load("XYtrain_XYvalid_processed.jld","YValid")
=#

Xtest = zeros(Float32,size(XValid,1),size(XValid,2),2,1)
Ytest = zeros(Float32,size(XValid,1),size(XValid,2),1,1)

idx_dir_test = 2
idx_slc_test = 15
idx_test = idx_slc_test + (idx_dir_test-1)*100
Xtest[:,:,:,1] = XValid[:,:,:,idx_test]
Ytest[:,:,:,1] = YValid[:,:,:,idx_test]
Xtest_output = u(gpu(Xtest))

fig, ax = plt.subplots()
PyPlot.gray()
ax.imshow(abs.(Xtest[:,:,1,1]))

fig2, ax2 = plt.subplots()
PyPlot.gray()
#savefig("Xttrain_40.png",fig)
ax2.imshow(abs.(Xtest_output[:,:,1,1]))

fig3, ax3 = plt.subplots()
PyPlot.gray()
ax3.imshow(abs.(Ytest[:,:,1,1]))

#BSON.@load string("Results/model_unet_ptx_processed_7dir_gpu_epoch",string(10),".bson") model
model = cpu(u)
slc_list = collect(1:100)
dir_list = collect(1:7)
XtestOutput4write = zeros(size(XValid,1),size(XValid,2),length(slc_list),length(dir_list))

Xtest = zeros(Float32,size(XValid,1),size(XValid,2),1,1)
Ytest = zeros(Float32,size(XValid,1),size(XValid,2),1,1)
for idx_dir_test in 1:length(dir_list)
  for idx_slc_test in 1:length(slc_list)
    idx_test = idx_slc_test + (idx_dir_test-1)*100
    Xtest[:,:,:,1] = XValid[:,:,:,idx_test]
    #Ytest[:,:,1,1] = YValid[:,:,idx_slc_test,idx_dir_test]
    XtestOutput4write[:,:,idx_slc_test,idx_dir_test] = model(Xtest)
  end
end

#=
slc_list = [10 25]
dir_list = [1 2]
XtestOutput4write_slc10_25_dir1_2 = zeros(Float32,size(XValid,1),size(XValid,2),length(slc_list),length(dir_list),10)
for idx_epoch in 1:10
  BSON.@load string("Results/model_unet_ptx_processed_7dir_gpu_epoch",string(idx_epoch),".bson") model
  Xtest = zeros(Float32,size(XValid,1),size(XValid,2),1,1)
  Ytest = zeros(Float32,size(XValid,1),size(XValid,2),1,1)
  for idx_dir_test in 1:length(dir_list)
    for idx_slc_test in 1:length(slc_list)
      idx_test = idx_slc_test + (idx_dir_test-1)*100
      Xtest[:,:,:,1] = XValid[:,:,:,idx_test]
      XtestOutput4write_slc10_25_dir1_2[:,:,idx_slc_test,idx_dir_test,idx_epoch] = model(Xtest)
    end
  end
end
=#

file = matopen("Results/XValidOutput_subj5_processed_7dir_gpu_orig_epoch10.mat", "w")
write(file, "XtestOutput4write", XtestOutput4write)
#write(file, "XtestOutput4write_slc10_25_dir1_2", XtestOutput4write_slc10_25_dir1_2)
close(file)

# testing on training data
BSON.@load string("model_unet_ptx_processed_7dir_epoch",string(10),".bson") u
XValid = Xtrain[:,:,:,1:700]
XValid = reshape(XValid,size(XValid,1),size(XValid,2),100,7)
XtestOutput4write = zeros(size(XValid))

Xtest = zeros(Float32,size(XValid,1),size(XValid,2),1,1)
Ytest = zeros(Float32,size(XValid,1),size(XValid,2),1,1)
for idx_dir_test in 1:size(XValid,4)
  for idx_slc_test in 1:size(XValid,3)
    Xtest[:,:,1,1] = XValid[:,:,idx_slc_test,idx_dir_test]
    #Ytest[:,:,1,1] = YValid[:,:,idx_slc_test,idx_dir_test]
    XtestOutput4write[:,:,idx_slc_test,idx_dir_test] = u(Xtest)
  end
end

file = matopen("XValidOutput_traingsubj1_processed_7dir_epoch10.mat", "w")
write(file, "XtestOutput4write", XtestOutput4write)
close(file)


text:
using Gadfly, DataArrays, Datetime, RDatasets

economics = dataset("HistData", "Prostitutes")
dates = Date[date(d) for d in economics[:Date]]

try
    economics[:Date] = dates
catch
    economics["Date"] = dates
end

p = plot(economics, x=:Date, y=:Count, Geom.line)


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

"""
    ConstValueDist <: Distributions.Distribution

Represents a delta distribution for a constant value of arbritrary type.

Calling `varshape` on a `ConstValueDist` will yield a
[`ConstValueShape`](@ref).
"""
struct ConstValueDist{VF<:VariateForm,T} <: Distribution{VF,Continuous}
    value::T
end

export ConstValueDist

ConstValueDist(x::T) where {T<:Real} = ConstValueDist{Univariate,T}(x)
ConstValueDist(x::T) where {T<:AbstractVector{<:Real}} = ConstValueDist{Multivariate,T}(x)
ConstValueDist(x::T) where {T<:AbstractMatrix{<:Real}} = ConstValueDist{Matrixvariate,T}(x)
ConstValueDist(x::NamedTuple{names}) where names = ConstValueDist{NamedTupleVariate{names},typeof(x)}(x)

_pdf_impl(d::ConstValueDist, x) = d.value == x ? float(eltype(d))(Inf) : float(eltype(d))(0)
_logpdf_impl(d::ConstValueDist, x) = d.value == x ? float(eltype(d))(Inf) : float(eltype(d))(-Inf)

Distributions.pdf(d::ConstValueDist{Univariate}, x::Real) = _pdf_impl(d, x)
Distributions.logpdf(d::ConstValueDist{Univariate}, x::Real) = _logpdf_impl(d, x)

Distributions._pdf(d::ConstValueDist{Multivariate}, x::AbstractVector{<:Real}) = _pdf_impl(d, x)
Distributions._logpdf(d::ConstValueDist{Multivariate}, x::AbstractVector{<:Real}) = log(pdf(d, x))

Distributions._pdf(d::ConstValueDist{Matrixvariate}, x::AbstractMatrix{<:Real}) = _pdf_impl(d, x)
Distributions._logpdf(d::ConstValueDist{Matrixvariate}, x::AbstractMatrix{<:Real}) = log(pdf(d, x))

Distributions.pdf(d::ConstValueDist{<:NamedTupleVariate{names}}, x::NamedTuple{names}) where names = _pdf_impl(d, x)
Distributions.logpdf(d::ConstValueDist{<:NamedTupleVariate{names}}, x::NamedTuple{names}) where names = log(pdf(d, x))

Distributions.insupport(d::ConstValueDist{Univariate}, x::Real) = x == d.value
Distributions.insupport(d::ConstValueDist{Multivariate}, x::AbstractVector{<:Real}) = x == d.value
Distributions.insupport(d::ConstValueDist{Matrixvariate}, x::AbstractMatrix{<:Real}) = x == d.value
Distributions.insupport(d::ConstValueDist{<:NamedTupleVariate{names}}, x::NamedTuple{names}) where names = x == d.value

Distributions.cdf(d::ConstValueDist{Univariate}, x::Real) = d.value <= x ? Float32(1) : Float32(0)
Distributions.quantile(d::ConstValueDist{Univariate}, q::Real) = d.value # Sensible?
Distributions.minimum(d::ConstValueDist{Univariate}) = d.value
Distributions.maximum(d::ConstValueDist{Univariate}) = d.value

StatsBase.mean(d::ConstValueDist) = d.value

StatsBase.mode(d::ConstValueDist) = d.value

Base.size(d::ConstValueDist{<:PlainVariate}) = size(d.value)
Base.length(d::ConstValueDist{<:PlainVariate}) = prod(size(d))
Base.eltype(d::ConstValueDist{<:PlainVariate}) = eltype(d.value)

Random.rand(rng::AbstractRNG, d::ConstValueDist) = d.value

Distributions._rand!(rng::AbstractRNG, d::ConstValueDist{<:Multivariate}, x::AbstractVector{<:Real}) = copyto!(x, d.value)
Distributions._rand!(rng::AbstractRNG, d::ConstValueDist{<:Matrixvariate}, x::AbstractMatrix{<:Real}) = copyto!(x, d.value)

Random.rand(rng::AbstractRNG, d::ConstValueDist{<:StructVariate}, dims::Dims) = Fill(d.value, dims)
Random.rand!(rng::AbstractRNG, d::ConstValueDist{<:StructVariate}, A::AbstractArray) = fill!(A, d.value)

ValueShapes.varshape(d::ConstValueDist) = ConstValueShape(d.value)

Statistics.var(d::ConstValueDist) = zero(d.value)


text:
using Base.Test

numfail = 0

# function mpi_test_printer(hdl::Task)
#     global numfail
#     for t = hdl
#         println("handler:", t)
#         if (t.succeed)
#             print(".")
#         else
#             numfail += 1
#             println("")
#             dump(t)
#             println("")
#         end
#     end
#     println("")
# end

function main()
    for x in ARGS
        include(x)
    end

    exit(numfail)
end

main()


text:
# Autogenerated wrapper script for SDL2_net_jll for x86_64-w64-mingw32
export libSDL2_net

using SDL2_jll
JLLWrappers.@generate_wrapper_header("SDL2_net")
JLLWrappers.@declare_library_product(libSDL2_net, "SDL2_net.dll")
function __init__()
    JLLWrappers.@generate_init_header(SDL2_jll)
    JLLWrappers.@init_library_product(
        libSDL2_net,
        "bin\\SDL2_net.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# This file holds the definition of the functions pertaining to the Attrition
#   type.

export  setAttritionName!,
        setAttritionRate!,
        setAttritionPeriod!,
        setAttritionCurve!

"""
```
setAttritionName!(
    attrition::Attrition,
    name::String )
```
This function sets the name of the attrition scheme `attrition` to `name`.

This function returns `true`, indicating the name has been set successfully.
"""
function setAttritionName!( attrition::Attrition, name::String )::Bool

    attrition.name = name
    return true

end  # setAttritionName!( attrition, name )

"""
```
setAttritionRate!(
    attrition::Attrition,
    rate::Real )
```
This function sets the attrition curve of the attrition scheme `attrition` to a constant value of `rate` per period for all terms. If the given rate is not in the interval [0.0, 1.0), the function doesn't make a change and gives a warning.

This function returns `true` if the rate was set successfully, `false` otherwise.
"""
function setAttritionRate!( attrition::Attrition, rate::Real )::Bool

    if !( 0.0 <= rate < 1.0 )
        @warn "Attrition rate must be in the range [0.0, 1.0), not making any changes."
        return false
    end  # if !( 0.0 <= rate < 1.0 )

    attrition.curvePoints = [0.0]
    attrition.rates = [rate]
    computeDistPars( attrition )
    return true

end  # setAttritionRate!( attrition, rate )

"""
```
setAttritionPeriod!(
    attrition::Attrition,
    period::Real )
```
This function sets the attrition period of the attrition scheme `attrition` to `period`. If the given period is ⩽ 0, the function doesn't make a change and gives a warning.

This function returns `true` if the period was set successfully, `false` otherwise.
"""
function setAttritionPeriod!( attrition::Attrition, period::Real )::Bool

    if period <= 0.0
        @warn "Attrition period must be > 0.0, not making any changes."
        return false
    end  # if period <= 0.0

    attrition.period = period
    computeDistPars( attrition )
    return true

end  # setAttritionPeriod!( attrition, period )

"""
```
setAttritionCurve!(
    attrition::Attrition,
    curve::Dict{Float64,Float64} )
```
This function sets the attrition curve of attrition scheme `attrition` to `curve`. The attrition curve is passed as a `Dict{Float64,Float64}` with the keys as the time a person exists in the simulation, and the values as the attrition rates per period from the time specified by that key to the next.

This function ignores non-sensical attrition rates, that is, rates outside the interval [0.0, 1.0). This function takes the rate at the negative time point closest to zero as the zero rate, and ignores the other negative time points.

If the first time point is non-zero positive, a zero attrition rate is set from time point 0 to that time point. If there are no eligible curve points, this function makes no changes.

This function returns `true` if the attrition rate curve is set successfully, `false otherwise`.
"""
function setAttritionCurve!( attrition::Attrition,
    curve::Dict{T1,T2} )::Bool where {T1 <: Real, T2 <: Real}

    terms = Float64.( sort( collect( keys( curve ) ) ) )
    rates = Float64.( map( term -> curve[term], terms ) )

    # Filter unsuitable attrition rates.
    isRateOkay = 0 .<= rates .< 1
    terms = terms[isRateOkay]
    rates = rates[isRateOkay]

    # If no suitable ones remain, make no changes.
    if isempty( terms )
        @warn "All attrition rates in the given curve are invalid, not making any changes."
        return false
    end  # if isempty( terms )

    # Clip the vectors to the last negative term point. If there is none, add a
    #   0 rate at the start.
    zeroInd = findlast( terms .<= 0 )

    if zeroInd isa Nothing
        terms = vcat( 0.0, terms )
        rates = vcat( 0.0, rates )
    else
        terms = terms[zeroInd:end]
        terms[1] = 0.0
        rates = rates[zeroInd:end]
    end  # if zeroInd isa Nothing

    # Filter out successive attrition rates that are the same.
    isDifferentFromPrevious = vcat( true,
        rates[2:end] .!= rates[1:(end - 1)] )
    attrition.curvePoints = terms[isDifferentFromPrevious]
    attrition.rates = rates[isDifferentFromPrevious]
    computeDistPars( attrition )
    return true

end  # setAttritionCurve!( attrition, curve )

"""
```
setAttritionCurve!(
    attrition::Attrition,
    curve::Array{Float64, 2} )
```
This function sets the attrition curve of attrition scheme `attrition` to
`curve`. The attrition curve is passed as a 2-dimensional `Array{Float64, 2}`
with 2 columns. The first column has the time a person exists in the simulation,
and the second column has the attrition rates per period from the time specified
by that key to the next.

If the dimensions of the curve array are incorrect, or if the curve array contains duplicate time points, this function makes no changes.

All remarks for the version of this function with a `Dict` are valid here as well.

This function returns `true` if the attrition rate curve is set successfully, `false otherwise`.
"""
function setAttritionCurve!( attrition::Attrition,
    curve::Array{T, 2} )::Bool where T <: Real

    # Check for correct dimensions.
    if size( curve, 2 ) != 2
        @warn "Invalid array given to define attrition curve, not making any changes."
        return false
    end  # if size( curve, 2 ) != 2

    # Check for duplicate terms.
    if length( curve[:, 1] ) != length( unique( curve[:, 1] ) )
        @warn "Duplicate entries in the terms of the attrition curve, not making any changes."
        return false
    end  # if length( curve[:, 1] ) != length( unique( curve[:, 1] ) )

    # Generate curve dictionary.
    curveDict = Dict{Float64,Float64}()

    for ii in eachindex( curve[:, 1] )
        curveDict[curve[ii, 1]] = curve[ii, 2]
    end  # for ii in eachindex( curve[:, 1] )

    return setAttritionCurve!( attrition, curveDict )

end  # setAttritionCurve!( attrition, curve )

function Base.show( io::IO, attrition::Attrition )::Nothing

    print( io, "Attrition scheme '", attrition.name, "'" )

    if all( attrition.rates .== 0.0 )
        print( io, "\n  No attrition." )
        return
    end  # if all( attrition.rates .== 0.0 )

    print( io, "\n  Attrition period: ", attrition.period )

    if length( attrition.rates ) == 1
        print( io, "\n  Attrition rate: ", attrition.rates[1] * 100.0, "%" )
    else
        print( io, "\n  Attrition curve" )

        for ii in eachindex( attrition.rates )
            print( io, "\n    term: ", attrition.curvePoints[ii], ";" )
            print( io, " rate: ", attrition.rates[ii] * 100, "%" )
        end  # for ii in eachindex( attrition.rates )
    end  # if length( attrition.rates ) == 1

    return

end  # Base.show( io, attrition )

include( joinpath( privPath, "attrition.jl" ) )


text:
# helper.jl
# utilities used in multiple test files

macro NOTinferred(ex) # flag where @inferred fails
    :($(esc(ex)))
end


text: include("discResonance.jl")