text:
using Penrose
using Test

# Helper functions
# Returns length of an edge between 2 pts
edgelen_fcn(p) = sqrt((p[2][1]-p[1][1])*(p[2][1]-p[1][1])+(p[2][2]-p[1][2])*(p[2][2]-p[1][2]))
# Returns angle (in degrees) of an edge between 2 pts
angle_fcn(p) = (360 + atan(p[2][2]-p[1][2], p[2][1]-p[1][1]) * 180 / pi) % 360
# Returns the mid-point of an edge between 2 pts
midpt_fcn(p) = [(p[2][1]+p[1][1])/2, (p[2][2]+p[1][2])/2]
# Returns the distance between 2 pts
dist_fcn(p1,p2) = sqrt((p2[1]-p1[1])*(p2[1]-p1[1]) + (p2[2]-p1[2])*(p2[2]-p1[2]))

@testset "Dart angle" begin

    # Angles outside [0-360] should get fixed in ctor
    @testset "angle range" begin
        @test Dart(0,0,-40).angle == 320
        @test Dart(0,0,400).angle ==  40
    end

    @testset "angle internal" begin
        internal_angles = [36, 72, 36, 216]
        d = Dart(0,0,0)
        @test (edge_angle(d,1) + (540 - internal_angles[1]))%360 == edge_angle(d,2)
        @test (edge_angle(d,2) + (540 - internal_angles[2]))%360 == edge_angle(d,3)
        @test (edge_angle(d,3) + (540 - internal_angles[3]))%360 == edge_angle(d,4)
        @test (edge_angle(d,4) + (540 - internal_angles[4]))%360 == edge_angle(d,1)
    end

    @testset "angle 0" begin
        d0 = Dart(0,0,0)
        @test edge_angle(d0, 1) == 252
        @test edge_angle(d0, 2) ==  36
        @test edge_angle(d0, 3) == 144
        @test edge_angle(d0, 4) == 288
    end

    @testset "angle +90" begin
        d90 = Dart(0,0,90)
        @test edge_angle(d90, 1) == 342
        @test edge_angle(d90, 2) == 126
        @test edge_angle(d90, 3) == 234
        @test edge_angle(d90, 4) ==  18
    end

    # -270 should be the same as +90
    @testset "angle -270" begin
        m270 = Dart(0,0,-270)
        @test edge_angle(m270, 1) == 342
        @test edge_angle(m270, 2) == 126
        @test edge_angle(m270, 3) == 234
        @test edge_angle(m270, 4) ==  18
    end

    @testset "angle +180" begin
        d180 = Dart(0,0,180)
        @test edge_angle(d180, 1) ==  72
        @test edge_angle(d180, 2) == 216
        @test edge_angle(d180, 3) == 324
        @test edge_angle(d180, 4) == 108
    end

    # -180 should be the same as +180
    @testset "angle -180" begin
        m180 = Dart(0,0,-180)
        @test edge_angle(m180, 1) ==  72
        @test edge_angle(m180, 2) == 216
        @test edge_angle(m180, 3) == 324
        @test edge_angle(m180, 4) == 108
    end

end

@testset "Kite angle" begin

    @testset "angle internal" begin
        internal_angles = [72, 144, 72, 72]
        k = Kite(0,0,0)
        @test (edge_angle(k,1) + (540 - internal_angles[1]))%360 == edge_angle(k,2)
        @test (edge_angle(k,2) + (540 - internal_angles[2]))%360 == edge_angle(k,3)
        @test (edge_angle(k,3) + (540 - internal_angles[3]))%360 == edge_angle(k,4)
        @test (edge_angle(k,4) + (540 - internal_angles[4]))%360 == edge_angle(k,1)
    end

    @testset "angle 0" begin
        k0 = Kite(0,0,0)
        @test edge_angle(k0, 1) == 324
        @test edge_angle(k0, 2) ==  72
        @test edge_angle(k0, 3) == 108
        @test edge_angle(k0, 4) == 216
    end

    @testset "angle 90" begin
        k90 = Kite(0,0,90)
        @test edge_angle(k90, 1) ==  54
        @test edge_angle(k90, 2) == 162
        @test edge_angle(k90, 3) == 198
        @test edge_angle(k90, 4) == 306
    end

    @testset "angle -270" begin
        m0 = Kite(0,0,-270)
        @test edge_angle(m0, 1) ==  54
        @test edge_angle(m0, 2) == 162
        @test edge_angle(m0, 3) == 198
        @test edge_angle(m0, 4) == 306
    end

    @testset "angle 180" begin
        k180 = Kite(0,0,180)
        @test edge_angle(k180, 1) == 144
        @test edge_angle(k180, 2) == 252
        @test edge_angle(k180, 3) == 288
        @test edge_angle(k180, 4) ==  36
    end

    @testset "angle -180" begin
        m180 = Kite(0,0,-180)
        @test edge_angle(m180, 1) == 144
        @test edge_angle(m180, 2) == 252
        @test edge_angle(m180, 3) == 288
        @test edge_angle(m180, 4) ==  36
    end

end

@testset "Dart edge slopes" begin
    d = Dart(0,0,0)

    @test angle_fcn(edge_points(d,1)) == edge_angle(d,1)
    @test angle_fcn(edge_points(d,2)) == edge_angle(d,2)
    @test angle_fcn(edge_points(d,3)) == edge_angle(d,3)
    @test angle_fcn(edge_points(d,4)) == edge_angle(d,4)
end

@testset "Kite edge slopes" begin
    k = Kite(0,0,0)

    @test angle_fcn(edge_points(k,1)) == edge_angle(k,1)
    @test angle_fcn(edge_points(k,2)) == edge_angle(k,2)
    @test angle_fcn(edge_points(k,3)) == edge_angle(k,3)
    @test angle_fcn(edge_points(k,4)) == edge_angle(k,4)
end

@testset "Dart edge lengths" begin
    d = Dart(0,0,0)

    l1 = (1+sqrt(5))/2
    l2 = 1 + l1

    @test edgelen_fcn(edge_points(d,1)) == l1
    @test edgelen_fcn(edge_points(d,2)) == l2
    @test edgelen_fcn(edge_points(d,3)) == l2
    @test edgelen_fcn(edge_points(d,4)) == l1
end

@testset "Kite edge lengths" begin
    k = Kite(0,0,0)

    l1 = (1+sqrt(5))/2
    l2 = 1 + l1

    @test edgelen_fcn(edge_points(k,1)) == l2
    @test edgelen_fcn(edge_points(k,2)) == l1
    @test edgelen_fcn(edge_points(k,3)) == l1
    @test edgelen_fcn(edge_points(k,4)) == l2
end

@testset "Dart edge center" begin
    d = Dart(0,0,0)

    @test dist_fcn(midpt_fcn(edge_points(d,1)), edge_center(d,1)) == 0
    @test dist_fcn(midpt_fcn(edge_points(d,2)), edge_center(d,2)) == 0
    @test dist_fcn(midpt_fcn(edge_points(d,3)), edge_center(d,3)) == 0
    @test dist_fcn(midpt_fcn(edge_points(d,4)), edge_center(d,4)) == 0
end

@testset "Kite edge center" begin
    k = Kite(0,0,0)

    @test dist_fcn(midpt_fcn(edge_points(k,1)), edge_center(k,1)) == 0
    @test dist_fcn(midpt_fcn(edge_points(k,2)), edge_center(k,2)) == 0
    @test dist_fcn(midpt_fcn(edge_points(k,3)), edge_center(k,3)) == 0
    @test dist_fcn(midpt_fcn(edge_points(k,4)), edge_center(k,4)) == 0
end

@testset "Dart geometry" begin
    phi = (1+sqrt(5))/2
    h = sqrt(5+2*sqrt(5))/2

    @testset "Compare geometry to edge_points" begin
        d = Dart(0,0,0)

        pts = geometry(d)
        @test size(pts) == (4,)
        @test dist_fcn(pts[1], edge_points(d,1)[1]) < 5e-8
        @test dist_fcn(pts[2], edge_points(d,1)[2]) < 5e-8
        @test dist_fcn(pts[2], edge_points(d,2)[1]) < 5e-8
        @test dist_fcn(pts[3], edge_points(d,2)[2]) < 5e-8
        @test dist_fcn(pts[3], edge_points(d,3)[1]) < 5e-8
        @test dist_fcn(pts[4], edge_points(d,3)[2]) < 5e-8
        @test dist_fcn(pts[4], edge_points(d,4)[1]) < 5e-8
        @test dist_fcn(pts[1], edge_points(d,4)[2]) < 5e-8
    end

    @testset "Angle = 0" begin
        d = Dart(0,0,0)

        pts = geometry(d)

        @test size(pts) == (4,)
        @test dist_fcn(pts[1], [0, 0]) < 5e-8
        @test dist_fcn(pts[2], [-1/2, -h]) < 5e-8
        @test dist_fcn(pts[3], [phi, 0]) < 5e-8
        @test dist_fcn(pts[4], [-1/2, h]) < 5e-8
    end

    @testset "Angle = 90" begin
        d = Dart(0,0,90)

        pts = geometry(d)

        @test size(pts) == (4,)
        @test dist_fcn(pts[1], [0, 0]) < 5e-8
        @test dist_fcn(pts[2], [h, -1/2]) < 5e-8
        @test dist_fcn(pts[3], [0, phi]) < 5e-8
        @test dist_fcn(pts[4], [-h,-1/2]) < 5e-8
    end

    @testset "origin = 10,10" begin
        d = Dart(10,10,0)

        pts = geometry(d)
        @test size(pts) == (4,)
        # Why's the tolerance worse here?
        @test dist_fcn(pts[1], [10, 10]) < 5e-7
        @test dist_fcn(pts[2], [10-(1/2), 10-h]) < 5e-7
        @test dist_fcn(pts[3], [10+phi, 10]) < 5e-7
        @test dist_fcn(pts[4], [10-(1/2), 10+h]) < 5e-7
    end
end

@testset "Kite geometry" begin
    k = Kite(0,0,0)

    pts = geometry(k)
    @test size(pts) == (4,)
    @test dist_fcn(pts[1], edge_points(k,1)[1]) < 5e-8
    @test dist_fcn(pts[2], edge_points(k,1)[2]) < 5e-8
    @test dist_fcn(pts[2], edge_points(k,2)[1]) < 5e-8
    @test dist_fcn(pts[3], edge_points(k,2)[2]) < 5e-8
    @test dist_fcn(pts[3], edge_points(k,3)[1]) < 5e-8
    @test dist_fcn(pts[4], edge_points(k,3)[2]) < 5e-8
    @test dist_fcn(pts[4], edge_points(k,4)[1]) < 5e-8
    @test dist_fcn(pts[1], edge_points(k,4)[2]) < 5e-8
end

@testset "Vertex 1" begin
    phi = (1+sqrt(5))/2

    d1 = Dart(-phi, 0, 0)
    d2 = placeDartEdge(3, edge_center(d1, 2), edge_angle(d1, 2))
    d3 = placeDartEdge(3, edge_center(d2, 2), edge_angle(d2, 2))
    d4 = placeDartEdge(3, edge_center(d3, 2), edge_angle(d3, 2))
    d5 = placeDartEdge(3, edge_center(d4, 2), edge_angle(d4, 2))

    h = sqrt(5+2*sqrt(5))/2
    k = (2+sqrt(5)) / (1+sqrt(5))
    p = sqrt(10 + sqrt(20))/4

    @test dist_fcn([d1.cx, d1.cy], [-phi, 0]) < 5e-8
    @test d1.angle == 0

    @test dist_fcn([d2.cx, d2.cy], [-1/2,-h]) < 5e-8
    @test d2.angle == 72

    @test dist_fcn([d3.cx, d3.cy], [k,-p]) < 5e-8
    @test d3.angle == 144

    @test dist_fcn([d4.cx, d4.cy], [k,p]) < 5e-8
    @test d4.angle == 216

    @test dist_fcn([d5.cx, d5.cy], [-1/2, h]) < 5e-8
    @test d5.angle == 288

end

@testset "Vertex 2" begin
    d1 = Dart(0, 0, 0)
    k1 = placeKiteEdge(2, edge_center(d1, 4), edge_angle(d1, 4))
    k2 = placeKiteEdge(4, edge_center(k1, 1), edge_angle(k1, 1))

    phi = (1+sqrt(5))/2
    h = sqrt(5+2*sqrt(5))/2
    k = (2+sqrt(5)) / (1+sqrt(5))
    p = sqrt(10 + sqrt(20))/4

    @test dist_fcn([d1.cx, d1.cy], [0, 0]) < 5e-8
    @test d1.angle == 0

    @test dist_fcn([k1.cx, k1.cy], [-k, p]) < 5e-8
    @test k1.angle == 36

    @test dist_fcn([k2.cx, k2.cy], [-k,-p]) < 5e-8
    @test k2.angle == 324
end

@testset "Vertex 3" begin
    phi = (1+sqrt(5))/2

    k1 = Kite(phi, 0, 0)
    k2 = placeKiteEdge(1, edge_center(k1, 4), edge_angle(k1, 4))
    k3 = placeKiteEdge(1, edge_center(k2, 4), edge_angle(k2, 4))
    k4 = placeKiteEdge(1, edge_center(k3, 4), edge_angle(k3, 4))
    k5 = placeKiteEdge(1, edge_center(k4, 4), edge_angle(k4, 4))

    h = sqrt(5+2*sqrt(5))/2
    k = (2+sqrt(5)) / (1+sqrt(5))
    p = sqrt(10 + sqrt(20))/4

    @test dist_fcn([k1.cx, k1.cy], [phi, 0]) < 5e-8
    @test k1.angle == 0

    @test dist_fcn([k2.cx, k2.cy], [1/2, h]) < 5e-8
    @test k2.angle == 72

    @test dist_fcn([k3.cx, k3.cy], [-k, p]) < 5e-8
    @test k3.angle == 144

    @test dist_fcn([k4.cx, k4.cy], [-k,-p]) < 5e-8
    @test k4.angle == 216

    @test dist_fcn([k5.cx, k5.cy], [1/2,-h]) < 5e-8
    @test k5.angle == 288
end

@testset "Vertex 4" begin
    phi = (1+sqrt(5))/2

    d1 = Dart(-phi, 0, 0)
    d2 = placeDartEdge(3, edge_center(d1, 2), edge_angle(d1, 2))
    k1 = placeKiteEdge(1, edge_center(d2, 2), edge_angle(d2, 2))
    k2 = placeKiteEdge(1, edge_center(k1, 4), edge_angle(k1, 4))
    d3 = placeDartEdge(3, edge_center(k2, 4), edge_angle(k2, 4))

    h = sqrt(5+2*sqrt(5))/2
    k = (2+sqrt(5)) / (1+sqrt(5))
    p = sqrt(10 + sqrt(20))/4

    @test dist_fcn([d1.cx, d1.cy], [-phi, 0]) < 5e-8
    @test d1.angle == 0

    @test dist_fcn([d2.cx, d2.cy], [-1/2,-h]) < 5e-8
    @test d2.angle == 72

    @test dist_fcn([k1.cx, k1.cy], [k,-p]) < 5e-8
    @test k1.angle == 324

    @test dist_fcn([k2.cx, k2.cy], [k, p]) < 5e-8
    @test k2.angle == 36

    @test dist_fcn([d3.cx, d3.cy], [-1/2, h]) < 5e-8
    @test d3.angle == 288
end

@testset "Vertex 5" begin
    k1 = Kite(-1, 0, 0)
    d1 = placeDartEdge(4, edge_center(k1, 2), edge_angle(k1, 2))
    k2 = placeKiteEdge(1, edge_center(d1, 3), edge_angle(d1, 3))
    k3 = placeKiteEdge(1, edge_center(k2, 4), edge_angle(k2, 4))
    d2 = placeDartEdge(2, edge_center(k3, 4), edge_angle(k3, 4))

    h = sqrt(5+2*sqrt(5))/2
    k = (2+sqrt(5)) / (1+sqrt(5))
    p = sqrt(10 + sqrt(20))/4

    @test dist_fcn([k1.cx, k1.cy], [-1, 0]) < 5e-8
    @test k1.angle == 0

    @test dist_fcn([d1.cx, d1.cy], [-1/2, -h]) < 5e-8
    @test d1.angle == 324

    @test dist_fcn([k2.cx, k2.cy], [ k, -p]) < 5e-8
    @test k2.angle == 324

    @test dist_fcn([k3.cx, k3.cy], [ k, p]) < 5e-8
    @test k3.angle == 36

    @test dist_fcn([d2.cx, d2.cy], [-1/2, h]) < 5e-8
    @test d2.angle == 36
end

@testset "Vertex 6" begin
    phi = (1+sqrt(5))/2

    d1 = Dart(-phi, 0, 0)
    k1 = placeKiteEdge(4, edge_center(d1, 2), edge_angle(d1, 2))
    k2 = placeKiteEdge(2, edge_center(k1, 3), edge_angle(k1, 3))
    k3 = placeKiteEdge(4, edge_center(k2, 1), edge_angle(k2, 1))
    k4 = placeKiteEdge(2, edge_center(k3, 3), edge_angle(k3, 3))

    h = sqrt(5+2*sqrt(5))/2
    k = (2+sqrt(5)) / (1+sqrt(5))
    p = sqrt(10 + sqrt(20))/4

    @test dist_fcn([d1.cx, d1.cy], [-phi, 0]) < 5e-8
    @test d1.angle == 0

    @test dist_fcn([k1.cx, k1.cy], [-1/2,-h]) < 5e-8
    @test k1.angle == 0

    @test dist_fcn([k2.cx, k2.cy], [ k,-p]) < 5e-8
    @test k2.angle == 216

    @test dist_fcn([k3.cx, k3.cy], [ k, p]) < 5e-8
    @test k3.angle == 144

    @test dist_fcn([k4.cx, k4.cy], [-1/2, h]) < 5e-8
    @test k4.angle == 0
end

@testset "Vertex 7" begin
    k = (2+sqrt(5)) / (1+sqrt(5))
    p = sqrt(10 + sqrt(20))/4

    k1 = Kite(k-1,-p, 108)
    k2 = placeKiteEdge(3, edge_center(k1, 2), edge_angle(k1, 2))
    d1 = placeDartEdge(4, edge_center(k2, 2), edge_angle(k2, 2))
    d2 = placeDartEdge(2, edge_center(d1, 3), edge_angle(d1, 3))

    @test dist_fcn([k1.cx, k1.cy], [k-1, -p]) < 5e-8
    @test k1.angle == 108

    @test dist_fcn([k2.cx, k2.cy], [k-1, p]) < 5e-8
    @test k2.angle == 252

    @test dist_fcn([d1.cx, d1.cy], [-k, p]) < 5e-8
    @test d1.angle == 216

    @test dist_fcn([d2.cx, d2.cy], [-k,-p]) < 5e-8
    @test d2.angle == 144
end


text:
function test_DAVIS_Conv(w; batchsize=32,saveImages=true)
    @time begin
        printlnFlush("Loading weights from $(w)")
        w = loadNetwork(w)
        colors = [Float32[0.0,0.0,0.0], Float32[1.0,1.0,1.0]]
        xtrain,ytrain,x,y = loadData_DAVIS()
        totalloss = totalcorrect = totalcount = imcounter = 0
        batchCount = Int64(ceil(size(y,4)/batchsize))
        for batch=1:batchCount
            printUpdate(@sprintf("Batch: %5d/%-5d",batch,batchCount))
            #Calculate current indices
            startIndex = 1+(batch-1)*batchsize
            endIndex = min(batch*batchsize,size(y,4))
            #Minibatch
            batchx = KnetArray(x[:,:,:,startIndex:endIndex])
            batchy = y[:,:,:,startIndex:endIndex]
            #Accuracy
            ypred = Array(DAVIS_Conv(w,batchx))
            totalloss += pixelwiseSoftloss(ypred,batchy)
            testAcc = pixelwiseAccuracy(ypred,batchy)
            totalcorrect += testAcc[1]
            totalcount += testAcc[2]
            #Image conversion
            if saveImages
                for im=1:size(ypred,4)
                    savePrediction(ypred[:,:,:,im],colors,outputPath(@sprintf("DAVIS/Conv/%05d.png",imcounter)))
                    imcounter += 1
                end
            end
        end
        printlnFlush(@sprintf("Total Loss: %0.10f\tTotal Accuracy: %0.4f%%",totalloss/batchCount,totalcorrect/totalcount*100))
    end
end

function test_DAVIS_LSTM(w; saveImages=true)
    @time begin
        printlnFlush("Loading weights from $(w)")
        @time w = loadNetwork(w)
        colors = [Float32[0.0,0.0,0.0], Float32[1.0,1.0,1.0]]
        totalloss = totalcorrect = totalcount = imcounter = 0
        xtrain,ytrain,x,y = loadData_DAVIS()
        batchCount = size(y,4)
        h = KnetArray(zeros(Float32,28,28,512,1))
        c = KnetArray(zeros(Float32,28,28,512,1))
        for t=1:batchCount
            printUpdate(@sprintf("Time: %5d/%-5d",t,batchCount))
            #Accuracy
            ypred = Array(DAVIS_LSTM(w,h,c,KnetArray(x[:,:,:,t:t])))
            totalloss += pixelwiseSoftloss(ypred,y[:,:,:,t:t])
            testAcc = pixelwiseAccuracy(ypred,y[:,:,:,t:t])
            totalcorrect += testAcc[1]
            totalcount += testAcc[2]
            #Image conversion
            if saveImages
                savePrediction(ypred[:,:,:,1],colors,outputPath(@sprintf("DAVIS/LSTM/%05d.png",imcounter)))
                imcounter += 1
            end
        end
        printlnFlush(@sprintf("Total Loss: %0.10f\tTotal Accuracy: %0.4f%%",totalloss/batchCount,totalcorrect/totalcount*100))
    end
end

function run_DAVIS_Conv_test()
    test_DAVIS_Conv("DAVIS/Conv.jld"; batchsize=32,saveImages=true)
end

function run_DAVIS_LSTM_test()
    test_DAVIS_LSTM("DAVIS/LSTM.jld"; saveImages=true)
end

export run_DAVIS_Conv_test,run_DAVIS_LSTM_test

ext: jl

max_stars_repo_name: eeyrdkl/MScThesis

max_issues_repo_name: eeyrdkl/MScThesis

max_forks_repo_name: eeyrdkl/MScThesis


text:
# CodePoints iterator
#
# Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones
# Licensed under MIT License, see LICENSE.md

struct CodePoints{T<:AbstractString}
    xs::T
end

"""
    codepoints(str)

An iterator that generates the code points of a string

# Examples
```jldoctest
julia> a = str("abc\U1f596")

julia> collect(a)

julia> collect(codepoints(a))
```
"""
codepoints(xs) = CodePoints(xs)
eltype(::Type{<:CodePoints{S}}) where {S} = eltype(S)
length(it::CodePoints) = length(it.xs)
@static if NEW_ITERATE
iterate(it::CodePoints, state=1) = iterate(it.xs, state)
else
start(it::CodePoints) = 1
done(it::CodePoints, state) = state > ncodeunits(it.xs)
next(it::CodePoints, state) = next(it.xs, state)
end


text:
using JuAFEM
using Tensors
using TimerOutputs
using UnicodePlots

const to = TimerOutput();
const dim = 3

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  hex = false
  geoshape = hex ? Hexahedron : Tetrahedron
  refshape = hex ? RefCube    : RefTetrahedron
  order = hex ? 2 : 1;

  corner1 = Vec{dim}((0.0, 0.0, 0.0))
  corner2 = Vec{dim}((10.0, 1.0, 1.0))
  grid = generate_grid(geoshape, (60, 6, 6), corner1, corner2);
  # Extract the left boundary
  addnodeset!(grid, "clamped", x -> norm(x[1]) ≈ 0.0);

  # Interpolations and values
  interpolation_space = Lagrange{dim, refshape, 1}()
  quadrature_rule = QuadratureRule{dim, refshape}(order)
  cellvalues = CellVectorValues(quadrature_rule, interpolation_space);
  facevalues = FaceVectorValues(QuadratureRule{dim-1, refshape}(order), interpolation_space);

  # DofHandler
  dh = DofHandler(grid)
  push!(dh, :u, dim) # Add a displacement field
  close!(dh)

  @time K = create_symmetric_sparsity_pattern(dh); # assemble only upper half since it is symmetric
  fill!(K.data.nzval, 1.0);
  spy(K.data)

  # Boundaryconditions
  dbc = DirichletBoundaryConditions(dh)
  # Add a homogenoush boundary condition on the "clamped" edge
  add!(dbc, :u, getnodeset(grid, "clamped"), (x,t) -> [0.0, 0.0, 0.0], collect(1:dim))
  close!(dbc)
  t = 0.0
  update!(dbc, t)

  # Create the stiffness tensor
  E = 200e9
  ν = 0.3
  λ = E*ν / ((1+ν) * (1 - 2ν))
  μ = E / (2(1+ν))
  δ(i,j) = i == j ? 1.0 : 0.0
  g(i,j,k,l) = λ*δ(i,j)*δ(k,l) + μ*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k))
  C = SymmetricTensor{4, dim}(g);

  function doassemble(cellvalues::CellVectorValues{dim}, facevalues::FaceVectorValues{dim},
                      K::Symmetric, grid::Grid, dh::DofHandler, C::SymmetricTensor{4, dim}) where dim

      f = zeros(ndofs(dh))
      assembler = start_assemble(K, f)

      n_basefuncs = getnbasefunctions(cellvalues)

      fe = zeros(n_basefuncs) # Local force vector
      Ke = Symmetric(zeros(n_basefuncs, n_basefuncs), :U) # Local stiffness mastrix

      t = Vec{3}((0.0, 1e8, 0.0)) # Traction vector
      b = Vec{3}((0.0, 0.0, 0.0)) # Body force
      ɛ = [zero(SymmetricTensor{2, dim}) for i in 1:n_basefuncs]
      @inbounds for (cellcount, cell) in enumerate(CellIterator(dh))
          @timeit to "assem" begin
          fill!(Ke.data, 0)
          fill!(fe, 0)

          reinit!(cellvalues, cell)
          for q_point in 1:getnquadpoints(cellvalues)
              for i in 1:n_basefuncs
                  ɛ[i] = symmetric(shape_gradient(cellvalues, q_point, i))
              end
              dΩ = getdetJdV(cellvalues, q_point)
              for i in 1:n_basefuncs
                  δu = shape_value(cellvalues, q_point, i)
                  fe[i] += (δu ⋅ b) * dΩ
                  ɛC = ɛ[i] ⊡ C
                  for j in i:n_basefuncs # assemble only upper half
                      Ke.data[i, j] += (ɛC ⊡ ɛ[j]) * dΩ # can only assign to parent of the Symmetric wrapper
                  end
              end
          end

          for face in 1:nfaces(cell)
              if onboundary(cell, face) && (cellcount, face) ∈ getfaceset(grid, "right")
                  reinit!(facevalues, cell, face)
                  for q_point in 1:getnquadpoints(facevalues)
                      dΓ = getdetJdV(facevalues, q_point)
                      for i in 1:n_basefuncs
                          δu = shape_value(facevalues, q_point, i)
                          fe[i] += (δu ⋅ t) * dΓ
                      end
                  end
              end
          end
          global_dofs = celldofs(cell)
          assemble!(assembler, global_dofs, fe, Ke)
          end # timer
      end
      return K, f
  end;

  global K, f, u
  reset_timer!(to)
  K, f = doassemble(cellvalues, facevalues, K, grid, dh, C);
  print_timer(to; linechars = :ascii)

  # Modify K and f such that K \ f gives correct boundary conditions
  @time apply!(K, f, dbc)

  @time u = cholesky(K) \ f;

  # Save file
  vtkfile = vtk_grid("cantilever", dh, u)
  vtk_save(vtkfile)

  Base.Test.@test maximum(u)  ≈ 1.919600482922295
  println("Cantilever successful")

end


text:
using Test
using ScottishTaxBenefitModel.Uprating
using DataFrames
using ScottishTaxBenefitModel
using .RunSettings: Settings
using .ModelHousehold
using .ExampleHouseholdGetter
using .FRSHouseholdGetter
using .Definitions
using .ExampleHelpers

prfr = Uprating.load_prices( Settings() )

print( prfr )

@time thesenames = ExampleHouseholdGetter.initialise( Settings() )

## NOTE this test has the 2019 OBR data and 2019Q4 as a target jammed on - will need
## changing with update versions

@testset "uprating tests" begin
    hh = ExampleHouseholdGetter.get_household( "mel_c2_scot" )
    hh.quarter = 1
    hh.interview_year = 2008
    pers = hh.people[SCOT_HEAD]
    # average index 2008 q1=100; 2019 Q4 = 125.9812039916
    pers.income[wages] = 100.0
    uprate!( hh )
    # FIXME hardly a test & needs changed every time the index changes
    @test pers.income[wages] ≈ 100*20.17/14.87 # 2022q2 av wages index
end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
# Only SOCP Constraints

using SparseArrays

include("../other_utils/feasibleCheck.jl")
include("projections.jl")

# -----------------------
# Second-Order Cone Constraints
# -----------------------
@doc raw"""
    AL_coneSlack(A, b, c, d)

Creates a Second Order Cone Constraint (SOCP) that uses slack variables.

We start with
```math
||Ax - b|| ≤ c^{\top} x - d
```

Using Slack variables, this becomes:
```math
||s|| ≤ t
```
```math
Ax - b = s
```
```math
c^{\top} x - d = t
```

So, in more standard form, this is becomes:
```math
||s|| - t ≤ 0
```
```math
Ax - b - s = 0
```
```math
c^{\top}x - d - t = 0
```

To check constraint satisfaction, use:
[`satisified`](@ref), [`whichSatisfied`](@ref), and [`coneSatisfied`](ref)

To evaluate the constraint, use:
- [`getRaw`](@ref) to evaluate each of `||s|| - t`, etc.
- [`getNormToProjVals`](@ref) to evaluate the projection to the affine
equalities and (more critically) the second order cone.

If the original `||Ax - b|| ≤ c'x - d` violation is needed, use
[`coneValOriginal`](@ref)

---

This ONLY works for second order cones. It will NOT work for any other p-cone,
i.e.
```math
||s||_2
```
works, but
```math
||s||_{10}
```
is not covered under this framework.
"""
struct AL_coneSlack <: constraint
    A
    b
    c
    d
end

# Evaluate the constraint.
"""
    getRaw(r::AL_coneSlack, x, s, t)

Evaluate the constraint without projection. See [`AL_coneSlack`](@ref)
"""
getRaw(r::AL_coneSlack, x, s, t) = [norm(s, 2) - t;
                                    r.A * x - r.b - s;
                                    r.c'x - r.d - t]

"""
    coneSatisfied(r::AL_coneSlack, s, t)

Checks if `||s|| - t ≤ 0`. See [`AL_coneSlack`](@ref)
"""
function coneSatisfied(r::AL_coneSlack, s, t)
    return (norm(s, 2) - t ≤ 0)
end

@doc raw"""
    coneValOriginal(r::AL_coneSlack, x)

Evaluates the raw cone value without slack variables, namely:
```math
||Ax - b|| ≤ c^{\top} x - d
```

See [`AL_coneSlack`](@ref)
"""
function coneValOriginal(r::AL_coneSlack, x)
    sInt = norm(r.A * x - r.b, 2)
    tInt = r.c'x - r.d
    return sInt - tInt
end

"""
    whichSatisfied(r::AL_coneSlack, x, s, t)

Check constraint satisfaction. See also [`satisfied`](@ref)
"""
function whichSatisfied(r::AL_coneSlack, x, s, t)
    raw = getRaw(r, x, s, t)
    return [raw[1] ≤ 0; raw[2:end] .== 0]
end

"""
    satisfied(r::AL_coneSlack, x, s, t)

Check constraint satisfaction. See also [`whichSatisfied`](@ref)
"""
function satisfied(r::AL_coneSlack, x, s, t)
    wSat = whichSatisfied(r, x, s, t)
    for iw in wSat
        if !iw
            # A single entry is false
            return false
        end
    end
    # All entries are true
    return true
end

"""
    coneActive(s, t, λ)

Checks if the cone constraint is "active." This prevents the solver from
getting stuck due to an inability to accurately update the dual λ.

We want to "activate" this constraint when `||s|| - t > 0` OR `λ > 0`.
"""
function coneActive(s, t, λ)
    #=
    returns "ACITVE" and "FILLED"

    "ACTIVE" = should return something nonzero
    "FILLED" = should project onto a solid cone, not the boundary

    We want to "activate" this constraint when ||s|| - t > 0 OR λ > 0.
    This yields 4 cases:

    if ||s|| - t > 0:
        OUTSIDE CONE
        if λ > 0:
            ACTIVE
            The solver is compensating for the constraint
        if λ = 0:
            ACTIVE
            The λ has not been initialized, but the constraint is active
    if ||s|| - t ≤ 0:
        INSIDE CONE
        if λ > 0:
            ACTIVE
            Must project to the boundary of the cone
        if λ = 0:
            INACTIVE
            Constraint is fully satisfied
    =#

    cVal = norm(s, 2) - t

    if cVal > 0
        # OUTSIDE CONE and ACTIVe
        active = true
        filled = true
    else
        # INSIDE CONE
        if λ > 0
            # ACTIVE
            # This is the "Solver might be stuck case."
            println("ACTIVE & NOT FILLED")
            active = true
            filled = false
        else
            # INACTIVE
            active = false
            filled = true
        end
    end

    return active, filled

end

"""
    getProjVecs(r::AL_coneSlack, x, s, t, filled = true, verbose = false)

Gets the projection vectors for each of the constraints that make up
[`AL_coneSlack`](@ref). In particular, we project onto the second order cone
and project each of the two affine equalities.

For the exact projections, see [`projSecondOrderCone`](@ref) and
[`projAffineEq`](@ref).

---

The parameter "filled" adjusts the projections onto the solid second order
cone (`||s|| ≤ t`) if `true` or the boundary of the second order cone
(`||s|| = t`) if `false`.

This is toggled based on [`coneActive`](@ref)
"""
function getProjVecs(r::AL_coneSlack, x, s, t, filled = true, verbose = false)
    #=
    We want to get each of the projection vectors.
    For the cone, we have the (s, t) cone given by ||s|| ≤ t.
    We activate this based on the function above

    For the equality constraints we have the affine equality projections
    Ax = b + s
    But this is equivalent to
    a1'x = b1 + s1, a2'x = b2 + s2, ..., am'x = bm + sm
    Where ak is the row k in A

    For the last equality constraint, we simply have c'x = d + t
    =#
    projVecs = []
    signs = []

    # Cone Constraints
    coneproj = projSecondOrderCone(s, t, filled)

    if verbose
        println("Cone Projection: $coneproj")
    end

    push!(projVecs, coneproj)

    # Set of Equality Constraints
    for ind in 1:size(r.b, 1)
        aI = r.A[ind, :]
        bI = r.b[ind]
        sI = s[ind]
        push!(projVecs, projAffineEq(aI, bI + sI, x))
        if (aI'x - bI - sI) ≥ 0
            push!(signs, 1)
        else
            push!(signs, -1)
        end
    end

    # Last Equality Constraint
    push!(projVecs, projAffineEq(r.c, r.d + t, x))
    if (r.c'x - r.d - t) ≥ 0
        push!(signs, 1)
    else
        push!(signs, -1)
    end

    return projVecs, signs
end

"""
    getNormToProjVals(r::AL_coneSlack, x, s, t, λ = 0)

Get the row by row constraint violation based on projections for the affine
constraints. This is better than [`getRaw`](@ref) since it uses the minimum
distance to the second-order cone rather than the 1D evaluation of `||s|| - t`.
It is also better than [`getRaw`](@ref) since it doesn't depend on
the magnitude of the elements in `r.A` and `r.b`.

See also [`AL_coneSlack`](@ref), [`getProjVecs`](@ref),
[`projSecondOrderCone`](@ref), [`projAffineEq`](@ref)

---

The parameter `λ` describes the dual of `||s|| - t` and is used to determine
if the constraint is active. See [`coneActive`](@ref) for more information.
"""
function getNormToProjVals(r::AL_coneSlack, x, s, t, λ = 0)
    #=
    We want to get the projected vector and calculate the distance between the
    original point and the constraint.
    =#
    active, filled = coneActive(s, t, λ)
    projVec, signs = getProjVecs(r, x, s, t, filled)#true)#filled)
    coneDiff = norm([s; t] - projVec[1], 2)
    if !filled#false#!filled
        coneDiff *= -1
    end
    projDiff = [sign * norm(pv - x, 2) for (sign, pv) in
                                            zip(signs, projVec[2:end])]
    return [coneDiff; projDiff]
end

# Lastly, we do some calculus
"""
    getGradC(r::AL_coneSlack, x, s, t, verbose = false)

Calculates the gradient of a constraint.

For `r::AL_coneSlack`, this is actually the Jacobian with respect to the
variables `x`, `s`, and `t` in that order.
"""
function getGradC(r::AL_coneSlack, x, s, t, verbose = false)
    #=
    We want the gradient of the constraints. This is piecewise due to the
    inequality case. This assumes a 2-norm (for now).

    First, let's consider the inequality constraints.
    → [ρ I_λ c(y) + λ]'∇c(y)

    where I_λ is 1 if [λ > 0 OR g(y) > 0] and 0 otherwise.
    ∇c(y) = [0; s/||s||; -1]

    Second, let's consider the equality constraints.
    → J(h(y))'(ρ h(y) + κ)
            x       s         t
    J = [ A      -1 * I       0  ]  h1
        [ c'       0         -1  ]  h2

    If we put these together, we get
            x       s         t
        [0       s/||s||     -1  ]  c1
    J = [ A      -1 * I       0  ]  h1
        [ c'       0         -1  ]  h2
    =#

    sizeJcols = size(x, 1) + size(s, 1) + size(t, 1)
    sizeJrows = size(t, 2) + size(s, 1) + size(t, 2)

    # if satisfied(r, x, s, t)
    #     if verbose
    #         println("Satisfied")
    #     end
    #     return zeros(sizeJrows, sizeJcols)

    # if coneSatisfied(r, s, t)
    #     jacobRow1 = zeros(1, sizeJcols)
    #     if verbose
    #         println("Cone Satisfied")
    #     end
    # else
    #     jacobRow1 = [zeros(size(x')) s'/norm(s,2) -1]
    # end

    jacobRow1 = [zeros(size(x')) s'/norm(s,2) -1]

    if verbose
        println("Row 1")
        display(jacobRow1)
    end
    jacobRow2 = [r.A (Diagonal(-1*ones(size(r.A, 1)))) zeros(size(r.A, 1))]
    if verbose
        println("Row 2")
        display(reshape(jacobRow2, size(r.A, 1), sizeJcols))
    end
    jacobRow3 = [r.c' zeros(1, size(r.A, 1)) -1]
    if verbose
        println("Row 3")
        display(jacobRow3)
    end
    jacob = [jacobRow1; jacobRow2; jacobRow3]
    if verbose
        println("Size expected = $sizeJrows x $sizeJcols ?= $(size(jacob))")
        display(reshape(jacob, size(jacob)))
        # display(reshape(jacob, sizeJrows, sizeJcols))
    end

    return jacob
end

"""
    getHessC(r::AL_coneSlack, x, s, t, λCone = 0)

Calculates the hessian of a constraint.

For `r::AL_coneSlack`, this is the hessian of the *constraint term* as it
appears in the augmented lagrangian. Namely, `H(ρc(x)'c(x) + λ c(x))`

returns H, but user must multiply by ρ.

See the source code for more technical details.
"""
function getHessC(r::AL_coneSlack, x, s, t, λCone = 0)
    #=
    We want the hessian of the combined constraint terms. Namely:
    ρc(x)'c(x) + λ c(x)

    Where
            [ ||s|| - t ]
    c(x) =  [ Ax - b - s]
            [c'x - d - t]

    Again, this assumes a 2-norm for now.

    We have three primal variables (x, s, t). Thus the Hessian is symmetric
    with this (n+m+1)x(n+m+1) dimensionality.

    We separate the Hessian into 6 submatrices
        [A  B  C]
    H = [B' D  E]
        [C' E' F]

    Focusing first on the cone inequality constraint
            n               m                    1
        [0              0                     0        ]   n
    H = [0              ss'/||s||)            -s/||s|| ]   m
        [0              -s'/||s||             1        ]   1

    When [λ > 0 OR g(y) > 0], otherwise it is uniformly zero.

    Then, we focus on the equality constraints
            n               m           1
        [A'A + cc'      -A'         -c      ]   n
    H = [-A             I_m         0       ]   m
        [-c'            0           1       ]   1

    =#

    sizeH = size(x, 1) + size(s, 1) + size(t, 1)

    # Equality Constraints
    A = r.A' * r.A + r.c * r.c'
    B = - r.A'
    C = - r.c
    D = Diagonal(ones(size(s, 1)))
    E = zeros(size(s))
    F = 1

    # if !coneSatisfied(r, s, t) | (λCone > 0)
    #     ns = norm(s, 2)
    #
    #     D += s * s' / (ns^2)
    #     E += s / ns
    #     F += 1
    # end

    # if !coneSatisfied(r, s, t) | (λCone > 0)
    #     ns = norm(s, 2)
    #
    #     # D += s * s' / (ns^2)
    #     E += - s / ns
    #     F += 1
    # end

    # THIS SHOULD BE THE CORRECT ONE
    if !coneSatisfied(r, s, t) | (λCone > 0)
        ns = norm(s, 2)

        D += s * s' / (ns^2)
        E += - s / ns
        F += 1
    end

    # if !coneSatisfied(r, s, t) | (λCone > 0)
    #     ns = norm(s, 2)
    #
    #     D += Diagonal(sign.(s))
    #     E += - s / ns
    #     F += 1
    # end

    # if !coneSatisfied(r, s, t) | (λCone > 0)
    #     ns = norm(s, 2)
    #
    #     D += Diagonal(s)
    #     E += - s / ns
    #     F += 1
    # end

    hess = [A B C; B' D E; C' E' F]
    return Symmetric(hess)
end

# -----------------------
# Second-Order Cone Constraints
# -----------------------
# Specifically has the form ||Ax - b||_p ≤ c'x - d
# Where p denotes which norm (usually p = 2)

"""
Don't use this version. It is not yet functional. See [`AL_coneSlack`](@ref)
instead
"""
struct AL_pCone <: constraint
    A
    b
    c
    d
    p
end

# Evaluate the constraint.
# Raw = Evaluate the function without projection
getRaw(r::AL_pCone, x) = norm(r.A * x - r.b, r.p) - r.c' * x + r.d

# Check that a constraint is satisfied (||Ax - b|| ≤ c'x - d)
satisfied(r::AL_pCone, x) = (getRaw(r, x) ≤ 0)
whichSatisfied(r::AL_pCone, x) = (getRaw(r, x) .≤ 0)

function getProjVecs(r::AL_pCone, x, verbose = false)
    #=
    We want to project onto the cone where
    v = ||Ax - b||
    s = cx - d
    =#
    v = r.A * x - r.b
    s = r.c' * x - r.d

    proj = projSecondOrderCone(v, s, r.p)

    if verbose
        println("x = $x -> Inside? $(satisfied(r, x))")
        println("v = $v, s = $s")
        println("proj = $proj")
    end

    vproj = proj[1:end - 1]
    xproj = r.A \ (r.b + vproj)

    if verbose
        println("vproj = $vproj -> $(norm(vproj, r.p))")
        print("xproj = $xproj -> $(r.A * xproj - r.b) -> ")
        print("$(norm(r.A * xproj - r.b, r.p))")
        println(" vs $(r.c' * xproj - r.d) vs $(proj[end])")
        println("Satisfied? $(satisfied(r, xproj))")
        println()
    end
    return xproj, proj
end

function getNormToProjVals(r::AL_pCone, x)
    #=
    We want to get the projected vector and calculate the distance between the
    original point and the constraint.
    =#
    projVec = getProjVecs(r, x)
    projDiff = projVec[1] - x
    return norm(projDiff, 2)
end

# Lastly, we do some calculus
function getGradC(r::AL_pCone, x, verbose = false)
    #=
    We want the gradient of the constraints. This is piecewise due to the
    inequality case. This assumes a 2-norm (for now)
    =#
    if satisfied(r, x)
        if verbose
            println("Satisfied")
        end
        return zeros(size(r.A, 2))
    else
        if verbose
            println("NOT Satisfied")
        end
        num = (r.A)' * (r.A * x - r.b)
        den = norm(r.A * x - r.b, 2) # Formula is not for a general norm yet
        return (num/den) - r.c
    end
end

function getHessC(r::AL_pCone, x)
    #=
    We want the hessian of the constraints. This is just zero for affine
    constraints, but nonzero for second order cone constraints.

    Again, this assumes a 2-norm for now.
    =#
    if satisfied(r, x)
        return zeros(size(r.A, 2), size(r.A, 2))
    end

    normVal = norm(r.A * x - r.b, 2)
    numGrad = (r.A)' * (r.A * x - r.b)

    term1 = ((r.A)' * r.A) / normVal  # Must be nxn
    term2 = (numGrad * numGrad') / (normVal^2) # Must be nxn
    return term1 + term2
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
module RichVehicleRoutingProblem

import JSON2
import Scanner

# Default values
const MAXNUMBER = 10^9

"""
    solve(data::RvrpInstance, solver)

Starts the solution procedure.
"""
function solve end

"""
    supported_features(solver_type)::Vector{BitSet}

Checks if the date is compatible with the solver provided.
"""
function supported_features end

include("RVRP_data.jl")
include("RVRP_computed_data.jl")
include("RVRP_solution.jl")
include("utils.jl")
include("parser_json.jl")
include("parser_cvrplib.jl")
include("instance_check.jl")
include("solution_check.jl")
include("helpers.jl")

# Solvers
include("vroom/interface.jl")
include("jsprit/interface.jl")
include("mock_solver/mock_solver.jl")

end


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb04xd(datfile, io=stdout)
    ZERO = 0.0e0
    NIN = 5
    NOUT = 6
    MMAX = 20
    NMAX = 20
    LDA = MMAX
    LDU = MMAX
    LDV = NMAX
    MAXMN = max( MMAX, NMAX )
    MNMIN = min( MMAX, NMAX )
    LENGQ = 2*MNMIN-1
    LDWORK = max( 2*NMAX, NMAX*( NMAX+1 )÷2 ) + max( 2*MNMIN + MAXMN, 8*MNMIN - 5 )
    A = Array{Float64,2}(undef, LDA,NMAX)
    U = Array{Float64,2}(undef, LDU,MMAX)
    V = Array{Float64,2}(undef, LDV,NMAX)
    Q = Array{Float64,1}(undef, LENGQ)
    INUL = Array{BlasBool,1}(undef, MAXMN)
    DWORK = Array{Float64,1}(undef, LDWORK)
f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    M = parse(BlasInt, vs[1])
    N = parse(BlasInt, vs[2])
    RANK = parse(BlasInt, vs[3])
    THETA = parse(Float64, replace(vs[4],'D'=>'E'))
    TOL = parse(Float64, replace(vs[5],'D'=>'E'))
    RELTOL = parse(Float64, replace(vs[6],'D'=>'E'))
    JOBU =  vs[7][1]
    JOBV =  vs[8][1]
if ( M<0 || M>MMAX )
elseif ( N<0 || N>NMAX )
elseif ( RANK>MNMIN )
elseif ( RANK<0 && THETA<ZERO )
else

    vs = String[]
    _isz,_jsz = (M,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
    RANK1 = RANK
    THETA1 = THETA
# interp call 1
    INUL .= 0
    RANK, THETA, INFO, IWARN = SLICOT.mb04xd!(JOBU, JOBV, M, N, RANK, THETA, A, U, V, Q, INUL, TOL, RELTOL)
    println(io, "RANK = $RANK")
    println(io, "THETA = $THETA")
    @test INFO == 0
    INFO == 0 || return
    println(io, "IWARN = $IWARN")
    inul = convert(Vector{Bool},INUL)

if ( INFO!=0 )
else
if ( IWARN!=0 )
else
end # if
    LJOBUA = LSAME( JOBU, 'A' )
    LJOBUS = LSAME( JOBU, 'S' )
    LJOBVA = LSAME( JOBV, 'A' )
    LJOBVS = LSAME( JOBV, 'S' )
    WANTU = LJOBUA||LJOBUS
    WANTV = LJOBVA||LJOBVS
    MINMN = min( M, N )
    LOOP = MINMN - 1
    println(io, "Q:")
    show(io, "text/plain", Bidiagonal(Q[1:MINMN],Q[MINMN+1:MINMN+LOOP],'U'))
    println(io)
if ( WANTU )
    NCOLU = LJOBUS ? MINMN : M
# interp output 1
    println(io, "U:")
    _nc = NCOLU
    _nr = M
    show(io, "text/plain", U[1:_nr,1:_nc])
    println(io)
    println(io, "INUL:")
    show(io, "text/plain", inul[1:NCOLU])
    println(io)

end # if
if ( WANTV )
    NCOLV = LJOBVS ? MINMN : N
# interp output 2
    println(io, "V:")
    _nc = NCOLV
    _nr = N
    show(io, "text/plain", V[1:_nr,1:_nc])
    println(io)
    println(io, "INUL:")
    show(io, "text/plain", inul[1:NCOLV])
    println(io)

end # if
end # if
end # if
    close(f)
end # run_mb04xd()


text:
"""
    SpatiotemporalMetrics

Functions for calculating spatiotemporal gait characteristics.

FS = Footstrike
FO = Foot liftoff

Stride time = time between FS of the same foot
Step time = time since previous opposite FS (e.g. left step time is the time between the current left step and the most recent right step)

Swing time = time from FO to FS of the same leg
Stance time = time from FS to FO of the same leg

Single support = time when only one foot is in contact with the ground; equivalent to the
sum of left and right swing when walking, and the sum of the left and right stance when
running

Double support = time when both feet are in contact with the ground; equivalent to the stride time minus the sum of the left and right swing when walking, and zero when running

Walking:

|RFS|     |LFO|     |LFS|     |RFO|     |RFS|     |LFO|     |LFS|
-----------------------------------------------------------------
  |         Right stance        | R swing |
  ---------------------------------------------------
            | L swing |            Left stance      |
  ---------------------------------------------------
  |  L DS   |         |  R DS   |         |
  -----------------------------------------
  |  Left step time   |  Right step time  |
  -----------------------------------------
  |             Stride time               |
  -----------------------------------------

Running:

|RFS|      |RFO|     |LFS|      |LFO|     |RFS|     |RFO|     |LFS|
-------------------------------------------------------------------
  | R stance |         Right swing          |  Float  |
  ---------------------------------------------------------------
             |  Float  | L stance |           Left swing        |
  ---------------------------------------------------------------
  |  Left step time    |  Right step time   |
  -------------------------------------------
  |             Stride time                 |
  -------------------------------------------
"""
module SpatiotemporalMetrics

using ..Biomechanics, LinearAlgebra
using ..Biomechanics: _rotating_diff

export beginwithevent, stridetimes, steptimes, swingstance, swing, stance, singlesupport,
    doublesupport, steplength, stepwidth

function beginwithevent(fs, events...)
    ffs = something(findlast(<(maximum(first, events)), fs), 1)
    pred = >(fs[ffs])
    fi = findfirst.(pred, events)

    return (fs[ffs:end], getindex.(events, range.(fi, lastindex.(events); step=1))...)
end

"""
    steplength(;lfs, rfs, lftpos, rftpos, [AP=1, VT=3, fs=1, ref=:fixed]) -> (;lsteps, rsteps)

Calculate left and right step lengths.

If the eltype's of `lfs` or `rfs` are Integers, they are assumed to be already in units of samples/indices.

# Keyword Arguments:
- `AP=1`: The column for the anteroposterior axis
- `VT=1`: The column for the vertical axis
- `fs=1`: The sampling frequency
"""
function steplength(;lfs, lfo, rfs, rfo, lftpos, rftpos, AP=1, VT=3, fs=1, requiredsteps=0.9)
    size(lftpos, 1) == size(rftpos, 1) ||
        throw(ArgumentError("left and right foot position data have unequal lengths"))

    if eltype(lfs) <: AbstractFloat
        lfs = toindices(lfs, fs)
    end
    if eltype(lfo) <: AbstractFloat
        lfo = toindices(lfo, fs)
    end
    if eltype(rfs) <: AbstractFloat
        rfs = toindices(rfs, fs)
    end
    if eltype(rfo) <: AbstractFloat
        rfo = toindices(rfo, fs)
    end

    max(last(lfs), last(rfs)) ≤ size(lftpos, 1) ||
        throw(ArgumentError("foot strike after the end of foot position data"))

    steplength(lfs, lfo, rfs, rfo; lftpos, rftpos, AP, VT, requiredsteps)
end

function findfloatingsteps(;lfs, lfo, rfs, rfo)
    events = sort!([tuple.(:rfs, rfs, axes(rfs,1)); tuple.(:lfs, lfs, axes(lfs,1));
        tuple.(:lfo, lfo, axes(lfo,1)); tuple.(:rfo, rfo, axes(rfo,1))], by=x->x[2])

    lfloat, rfloat = Int[], Int[]
    for i in 1:lastindex(events)
        if events[i][1] ∈ (:rfs, :lfs)
            opp = events[i][1] === :rfs ? :lfs : :rfs
            fo = events[i][1] === :rfs ? :lfo : :rfo
            j = findprev(==(opp)∘first, events, i)
            if !isnothing(j)
                if fo ∈ first.(@view(events[j:i]))
                    if events[i][1] === :rfs
                        # @show events[i][1], opp, fo, events[j:i]
                        push!(rfloat, events[i][3])
                    else
                        # @show events[i][1], opp, fo, events[j:i]
                        push!(lfloat, events[i][3])
                    end
                end
            else
                if events[i][1] === :rfs
                    push!(rfloat, events[i][3])
                else
                    push!(lfloat, events[i][3])
                end
            end
        end
    end

    lfloat, rfloat
end

function steplength(
    lfs::AbstractVector{Int},
    lfo::AbstractVector{Int},
    rfs::AbstractVector{Int},
    rfo::AbstractVector{Int};
    lftpos, rftpos, AP=1, VT=3, requiredsteps=0.9
)
    size(lftpos, 1) == size(rftpos, 1) ||
        throw(ArgumentError("left and right foot position data have unequal lengths"))
    max(last(lfs), last(rfs)) ≤ size(lftpos, 1) ||
        throw(ArgumentError("foot strike after the end of foot position data"))
    lfloat, rfloat = findfloatingsteps(;lfs, lfo, rfs, rfo)
    minlsteps = requiredsteps isa Int ? requiredsteps : round(Int, requiredsteps*length(lfs))
    minrsteps = requiredsteps isa Int ? requiredsteps : round(Int, requiredsteps*length(rfs))
    # @show minlsteps minrsteps lfloat rfloat length(lfs) length(rfs)
    if (length(lfs) - length(lfloat) < minlsteps) | (length(rfs) - length(rfloat) < minrsteps)
        throw(ArgumentError("number of floating steps exceed limits; try lowering the number of required steps?"))
    end

    if VT === nothing
        slc = [AP]
    else
        slc = sort([AP,VT])
    end

    goodrfs = rfs[setdiff(axes(rfs, 1), rfloat)]
    goodlfs = lfs[setdiff(axes(lfs, 1), lfloat)]

    rsteps = vec(mapslices(norm, rftpos[goodrfs, slc] - lftpos[goodrfs, slc]; dims=2))
    lsteps = vec(mapslices(norm, lftpos[goodlfs, slc] - rftpos[goodlfs, slc]; dims=2))

    return (;lsteps, rsteps)
end

function stepwidth(;lfs, rfs, lftpos, rftpos, ML=2, fs=1)
    size(lftpos, 1) == size(rftpos, 1) ||
        throw(ArgumentError("left and right foot position data have unequal lengths"))

    lfs = toindices(lfs, fs)
    rfs = toindices(rfs, fs)

    max(last(lfs), last(rfs)) ≤ size(lftpos, 1) ||
        throw(ArgumentError("foot strike after the end of foot position data"))

    stepwidth(lfs, rfs; lftpos, rftpos, ML)
end

function stepwidth(lfs::AbstractVector{Int}, rfs::AbstractVector{Int};
    lftpos, rftpos, ML=2
)
    size(lftpos, 1) == size(rftpos, 1) ||
        throw(ArgumentError("left and right foot position data have unequal lengths"))
    max(last(lfs), last(rfs)) ≤ size(lftpos, 1) ||
        throw(ArgumentError("foot strike after the end of foot position data"))

    if first(lfs) > first(rfs)
        lsteps, rsteps = _rotating_diff(rftpos[rfs, ML], lftpos[lfs, ML]; i=1)
    else
        rsteps, lsteps = _rotating_diff(lftpos[lfs, ML], rftpos[rfs, ML]; i=1)
    end

    return (;lsteps, rsteps)
end

stridetimes(fs) = diff(fs)

"""
    steptimes(;lfs, rfs) -> (;lsteps, rsteps)

Calculate the left and right step times, where a left step is the left foot-strike
following a right foot-strike, and vice versa.
"""
function steptimes(;lfs, rfs, fs=1)
    if first(lfs) > first(rfs)
        lsteps, rsteps = _rotating_diff(rfs, lfs)
    else
        rsteps, lsteps = _rotating_diff(lfs, rfs)
    end

    lsteps ./= fs
    rsteps ./= fs

    return (;lsteps, rsteps)
end

"""
    swingstance(fs, fo; normalize=true) -> (;swing, stance)

Calculate the swing and stance times beginning with the first fs.

If `normalize` is true, swing and stance times will be normalized to the stride time given by `fs`. This will strip any `fo` preceding the first `fs` and any `fo` following the last `fs`. (This would remove an initial swing and/or final stance.)
"""
function swingstance(fs, fo; normalize=true)
    if normalize
        fs, fo = beginwithevent(fs, fo)
    end

    if first(fo) < first(fs)
        sw, st = _rotating_diff(fo, fs)
    else
        st, sw = _rotating_diff(fs, fo)
    end

    if normalize
        l = min(length(st), length(sw))
        resize!(st, l)
        resize!(sw, l)

        strides = stridetimes(fs)
        sw ./= strides
        st ./= strides
    end

    return (;swing=sw, stance=st)
end

function swing(fs, fo; normalize=true)
    sw, _ = swingstance(fs, fo; normalize)
    return sw
end

function stance(fs, fo; normalize=true)
    _, st = swingstance(fs, fo; normalize)
    return st
end

"""
    singlesupport(lfs, lfo, rfs, rfo)

Calculate the single support times based on right strides and beginning with the first rfs.

See also: [`swingstance`](@ref), [`swing`](@ref), [`stance`](@ref)
"""
function singlesupport(;lfs, lfo, rfs, rfo, normalize=true)
    rfs, lfo, lfs, rfo = beginwithevent(rfs, lfo, lfs, rfo)

    lsteps, rsteps = steptimes(;rfs,lfs)
    lswings, lstances = swingstance(lfs, lfo; normalize=false)
    rswings, rstances = swingstance(rfs, rfo; normalize=false)
    l = minimum(lastindex, (lsteps, rsteps, lswings, lstances, rswings, rstances))

    ss = @. @views ifelse(rstances[1:l] > lsteps[1:l], rswings[1:l], rstances[1:l]) +
        ifelse(lstances[1:l] > rsteps[1:l], lswings[1:l], lstances[1:l])

    any(≤(0), ss) && throw(error("Single support cannot be ≤ 0. Check your gait events"))

    if normalize
        strides = stridetimes(rfs)
        ss ./= strides[1:l]
    end

    return ss
end

"""
    doublesupport(lfs, lfo, rfs, rfo)

Calculate the double support times based on right strides and beginning with the first rfs.

See also: [`swingstance`](@ref), [`swing`](@ref), [`stance`](@ref)
"""
function doublesupport(;lfs, lfo, rfs, rfo, normalize=true)
    rfs, lfo, lfs, rfo = beginwithevent(rfs, lfo, lfs, rfo)
    events = sort!([tuple.(:rfs, rfs); tuple.(:lfs, lfs); tuple.(:lfo, lfo);
                    tuple.(:rfo, rfo)], by=x->x[2])
    ds = Float64[]
    lds = false
    for i in 1:lastindex(events)-1
        if events[i][1] === :rfs && events[i+1][1] === :lfo # left double-support (DS)
            push!(ds, events[i+1][2] - events[i][2])
            lds = true
        elseif !isempty(ds) && events[i][1] === :lfs && events[i+1][1] === :rfo # right DS
            if lds
                ds[end] += (events[i+1][2] - events[i][2])
                lds = false
            else
                push!(ds, events[i+1][2] - events[i][2])
                lds = true
            end
        elseif events[i+1] === :rfs
            lds = false
        end
    end

    if normalize
        strides = stridetimes(rfs)
        l = min(lastindex(strides), lastindex(ds))
        ds = ds[1:l] ./ strides[1:l]
    end

    return ds
end

function swingintervals(fs::AbstractVector{T}, fo::AbstractVector{T}; step=one(T)) where T
    intervals(fo, fs; step)
end

function stanceintervals(fs::AbstractVector{T}, fo::AbstractVector{T}; step=one(T)) where T
    intervals(fs, fo; step)
end

end


text:
function SeisPSTM(m::Array{Float32,2},d::Array{Float32,2},param::Dict{Any,Any})

	v = get(param,"v",fill(1500.0f0,size(m)))
	adj = get(param,"adj",true)
	nsinc = int32(get(param,"nsinc",5))
	nsinc = isodd(nsinc) ? int32(nsinc) : int32(nsinc - 1)
	aperture = float32(get(param,"aperture",1000.0))
	sx = get(param,"sx",fill(0.0f0,size(d,2)))
	sy = get(param,"sy",fill(0.0f0,size(d,2)))
	gx = get(param,"gx",fill(0.0f0,size(d,2)))
	gy = get(param,"gy",fill(0.0f0,size(d,2)))
	nt = int32(get(param,"nt",size(m,1)))
	ot = float32(get(param,"ot",0.0))
	dt = float32(get(param,"dt",0.001))
	nx = int32(get(param,"nx",size(m,2)))
	ox = float32(get(param,"ox",0.0))
	dx = float32(get(param,"dx",10.0))
	ny = int32(get(param,"ny",1))
	oy = float32(get(param,"oy",0.0))
	dy = float32(get(param,"dy",10.0))

	for itrace = 1 : size(d,2)
		if (adj == true)
			trace = d[:,itrace]
		end
		pstm_op(m,trace,v,sx[itrace],sy[itrace],gx[itrace],gy[itrace],nx,ox,dx,ny,oy,dy,nt,ot,dt,aperture,nsinc,adj)
		if (adj == false)
			d[:,itrace] = trace
		end
	end

end

function pstm_op(m::Array{Float32,2},trace::Array{Float32,1},v::Array{Float32,2},sx::Float32,sy::Float32,gx::Float32,gy::Float32,nx::Int32,ox::Float32,dx::Float32,ny::Int32,oy::Float32,dy::Float32,nt::Int32,ot::Float32,dt::Float32,aperture::Float32,nsinc::Int32,adj::Bool)

	trace = PhaseShift(trace,pi/2)
	trace1 = trace.*0
	trace2 = trace.*0
	Seismic.integrate(length(trace),trace,trace1,false)
	Seismic.integrate(length(trace),trace2,trace1,true)
	cmpx = (sx + gx)/2.0f0
	cmpy = (sy + gy)/2.0f0

	ixi = max(int32(1),int32((cmpx-aperture-ox)/dx))
	ixf = min(int32((cmpx+aperture-ox)/dx),nx)
	iyi = max(int32(1),int32((cmpy-aperture-oy)/dy))
	iyf = min(int32((cmpy+aperture-oy)/dy),ny)

	if (adj)
		for ix = ixi : ixf
			x = ox + ix*dx
			distx = abs(x - cmpx)
			for iy = iyi : iyf
				y = oy + iy*dy
				disty = abs(y - cmpy)
				dist = sqrt(distx^2 + disty^2)
				dists = sqrt((sx-x)^2 + (sy-y)^2)
				distg = sqrt((gx-x)^2 + (gy-y)^2)
				dists2 = dists*dists
				distg2 = distg*distg
				for it = 1 : nt
					t0 = 0.5f0*(ot + it*dt)
					t02 = t0*t0
					v2 = v[it,(ix-1)*ny + iy]*v[it,(ix-1)*ny + iy]
					ts = sqrt(t02 + dists2/v2)
					tg = sqrt(t02 + distg2/v2)
					t = ts + tg
					geoms = sqrt(1.0f0/(t*v[it,(ix-1)*ny + iy]));
					obliq = sqrt(0.5f0*(1.0f0 + (t0*t0/(4.0f0*ts*tg)) - (1.0f0/(ts*tg))*sqrt(ts*ts - t0*t0/4.0f0)*sqrt(tg*tg - t0*t0/4.0f0)));
					fmax = 0.5f0*v2/dx/(dists/ts + distg/tg);
					jt_float = (t-ot)/dt
					jt = int(jt_float)
					k = max(round(1/fmax/dt-1),0)
					it1 = jt-int(k)-1
					it2 = jt+int(k)+1
					if nsinc < it1 && it2 < nt - nsinc
						for isinc = - (nsinc - 1)/2 : (nsinc - 1)/2
							m[it,(ix-1)*ny + iy] += -geoms*obliq*sinc(float32( (isinc)/nsinc ))*trace2[it1 + isinc]/((k+1)^2)
							m[it,(ix-1)*ny + iy] += 2*geoms*obliq*sinc(float32( (isinc)/nsinc ))*trace2[jt + isinc]/((k+1)^2)
							m[it,(ix-1)*ny + iy] += -geoms*obliq*sinc(float32( (isinc)/nsinc ))*trace2[it2 + isinc]/((k+1)^2)
						end
						#m[it,(ix-1)*ny + iy] += 2*geoms*obliq*trace[jt]
					end
				end
			end
		end
	else
		error("write it now")
	end

end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
using Distributions, BulkAndTails

# Some fake data, in this case a simple Cauchy distribution.
const data = rand(TDist(1.0), 5_000)

# you can set print_level=0 to have no Ipopt output. But it is nice to have the
# output, because if something goes wrong you can google "interpreting Ipopt
# output" and then look at the abundant printed information and probably figure
# out what's going wrong with your model.
#
# Note that Ipopt is pretty sensitive, and so there is a chance that this will
# kind of get stuck or something strange seeming and hit the maximum number of
# iterations or something. This isn't necessarily a problem with the model or
# code. It's a hard optimization problem. We suggest trying again with different
# initial values (an optional second arg here).
const (mle, observed_information_matrix) = fit_mle(BulkAndTailsDist, data, print_level=5)

# A vague check for fit quality, which seems pretty good.
for _x in (0.1, 1.0, 10.0, 100.0, 1000.0)
  println("x = $_x:")
  println("Cauchy PDF:      $(pdf(TDist(1.0), _x))")
  println("Fitted BATs PDF: $(pdf(mle,        _x))")
end


text:
using StaticArrays

@testset "StaticArrays" begin
    X = @SMatrix rand(10,4);
    S = @SMatrix(randn(4,7)) |> x -> Symmetric(x * x');
    PDS = PDMat(S)
    @test isbits(PDS)
    @test isbits(X / PDS)
    @test isbits(PDS \ X')
    @test (X / PDS) ≈ Matrix(X) / Matrix(S)
    @test (PDS \ X') ≈ Matrix(S) \ Matrix(X')

    D = PDiagMat(@SVector(rand(4)))
    @test isbits(D)
    @test isbits(X / D)
    @test isbits(D \ X')
    @test (X / D) ≈ Matrix(X) / Matrix(D)
    @test (D \ X') ≈ Matrix(D) \ Matrix(X')
end


text:
function wrap!(x,y,side)
  for i in 1:2
    dx = (y[i] - x[i])%side
    if dx > side/2
      dx = dx - side
    elseif dx < -side/2
      dx = dx + side
    end
    y[i] = x[i] + dx
  end
end

export wrap!


text:
function writeNodeFluxData( fluxReport::Dict{String,DataFrame},
    nodes::NTuple{N,String}, fluxType::Symbol, mpSim::MPsim,
    timeFactor::Float64, reportGenerationTime::Float64, filename::String,
    overwrite::Bool ) where N

    tStart = now()
    filename = string( filename, endswith( filename, ".xlsx") ? "" : ".xlsx" )

    XLSX.openxlsx( filename, mode = overwrite ? "w" : "rw" ) do xf
        dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
        ws = xf[1]
        nReports = 1

        # Prep the new sheet.
        if overwrite
            XLSX.rename!( ws, string( "Report 1  ", dateStr ) )
        else
            nReports = length( filter( name -> startswith( name, "Report" ),
                XLSX.sheetnames( xf ) ) ) + 1
            ws = XLSX.addsheet!( xf, string( "Report ", nReports, "  ",
                dateStr ) )
        end  # if overwrite

        # Sheet header
        ws["A1"] = "Simulation length"
        ws["B1"] = mpSim.simLength / timeFactor
        ws["C1"] = "years"
        ws["A2"] = "Data generation time"
        ws["B2"] = reportGenerationTime
        ws["C2"] = "seconds"
        ws["A3"] = "Excel generation time"
        ws["C3"] = "seconds"
        ws["A4"] = "Flux type"
        ws["B4"] = string( fluxType )
        ws["A5"] = "Nodes"
        ii = 4

        for node in filter( node -> haskey( fluxReport, node ),
            vcat( nodes... ) )
            ii += 1
            ws[ii, 2] = node

            rws = XLSX.addsheet!( xf, string( node, " ", fluxType, "  Report ",
                nReports ) )
            report = fluxReport[node]
            report[:, :timeStart] /= timeFactor
            report[:, :timePoint] /= timeFactor
            nrows, ncols = size( report )
            headers = names( report )

            for jj in 1:ncols
                rws[1, jj] = jj == 1 ? "Period start" : headers[jj]

                for kk in 1:nrows
                    rws[kk + 1, jj] = report[kk, jj]
                end  # for kk in 1:nrows
            end  # for jj in 1:ncols
        end  # for node in filter( ... )

        excelGenerationTime = ( now() - tStart ).value / 1000.0
        ws["B3"] = excelGenerationTime
    end  # XLSX.openxlsx( ... ) do xf

end  # writeNodeFluxData( fluxReport, nodes, fluxType, mpSim, timeFactor,
     #   reportGenerationTime, filename, overwrite )


text:
using VoronoiDelaunay
using VoronoiCells
using Winston

N = 15
x = 1.0 + rand(N)
y = 1.0 + rand(N)
tess1 = DelaunayTessellation2D{IndexablePoint2D}()
pts1 = [IndexablePoint2D(x[n], y[n], n) for n=1:N]
push!(tess1, pts1)

# ----------------------------------------------------------------------
# Plot tesselations

Vx1, Vy1 = getplotxy(VoronoiDelaunay.voronoiedges(tess1))
p1 = plot(Vx1, Vy1)

for p in pts1
	text( getx(p), gety(p), string(getindex(p)) )
end

xlim(1.0, 2.0)
ylim(1.0, 2.0)

#display(p1)
savefig("tess_corner.png")

# ----------------------------------------------------------------------
# Scaled pints

tess2 = DelaunayTessellation2D{IndexablePoint2D}()
pts2 = large2small(pts1)
push!(tess2, pts2)

# ----------------------------------------------------------------------
# Plot scaled tesselation

Vx2, Vy2 = getplotxy(VoronoiDelaunay.voronoiedges(tess2))
p2 = plot(2.0*Vx2-1.5, 2.0*Vy2-1.5, "r")

for p in pts1
	text( getx(p), gety(p), string(getindex(p)) )
end

xlim(1.0, 2.0)
ylim(1.0, 2.0)

#display(p2)
savefig("tess_nocorner.png")

# ----------------------------------------------------------------------
# Both tesselations

#oplot(Vx1, Vy1)

#display(p2)
#savefig("tess_combined.png")


text:
# greedy pairs heuristic by A. Catalan and M. Fisher (2012) doi:10. 2139/ssrn.2166687
function GREEDYPAIRS(Q::Matrix{Int64},
                     capacity::Array{Int64,1})
    if CHECKCAPACITY(Q::Array{Int64,2}, capacity::Array{Int64,1}) == 1
        ## Calculate the coapperance matrix
        capacity_left = copy(capacity)
        ## Sort the pairs of distinct SKUs by decreasing coappearance
        pairs = SORTPAIRS(Q::Array{Int64,2})
        X = Array{Int64,2}(undef,size(Q,2),size(capacity,1)) .= 0
        GREEDYPAIRSMAIN!(pairs::Array{Int64,2}, X::Array{Int64,2}, capacity_left::Array{Int64,1})
        ## While there is unallocated space in the DCs do
        X = FILLUP(X::Array{Int64,2}, Q::Array{Int64,2}, capacity_left::Array{Int64,1})
        return X::Array{Int64,2}
    end
end

max_stars_repo_name: beyondatlas/splitdeliveries

max_issues_repo_name: beyondatlas/splitdeliveries

max_forks_repo_name: beyondatlas/splitdeliveries


text:
module Localize

using Compat, MacroTools

#######################
#    Language Strings

# string literals
export @zh_str, @fr_str, @lang

# additional language string interface
export language, locale

include("strings.jl")

########################
#     Doc System Utils

export get_multidoc, get_sigs, submodules, path, paths, allpaths, dump_str, replace!, generate

include("docsys.jl")
include("generated.jl")

export Translation, translate
include("text.jl")

export @i18n

include("syntax.jl")

end # module


text:
# fld-write.jl

using MIRT: fld_write, fld_read
#import MIRT: fld_write_data_fix

#using MIRT: ir_test_dir, ir_test_dir!
using Test: @test, @testset, @test_throws, @inferred

#=

"""
`fld_write_test1(file, data, ...)`
"""
function fld_write_test1(file, data ;
	raw::Bool=false,
	warn::Bool=false,
	chat::Bool=false, kwarg...
)
	@inferred fld_write(file, data ; raw, warn, kwarg...)
	tmp = fld_read(file ; chat)
	tmp != data && throw("test failed for file = $file")

	rm(file)
	if raw
		tmp = file[1:(end-4)] * ".raw"
		rm(tmp)
	end
	true
end

# test writing and reading AVS `.fld` files of all key types

dir = mktempdir()
file = joinpath(dir, "fld-write-test.fld")
chat = isinteractive()

@testset "write std" begin
	for dtype in (UInt8, Int16, Int32, Float32, Float64)
		for endian in (:le, :be)
			for raw in (false, true)
				chat && @info "dtype=$dtype endian=$endian raw=$raw"
				data = convert.(dtype, [5:8; 1:4])
				@test fld_write_test1(file, data ; endian, raw,
					head = ["dtype=$dtype endian=$endian raw=$raw",],
					check=true, chat=false)
			end
		end
	end
end

@testset "write odd" begin
	for T in (BigFloat, Float16, BigInt, Int64, Bool)
		data = T.([1 0 1])
		@test fld_write_test1(file, data ; check=true, chat=false)
	end
end

@testset "fld convert" begin
#=
	pairs = (
		(BigFloat, Float64),
		(Float16, Float32),
		(BigInt, Int32),
		(Int64, Int32),
		(Bool, UInt8),
	)
	for (in, out) in pairs
		data = in.([1 0 1])
		@test eltype(fld_write_data_fix(data)) == out
	end
=#
	@test_throws ArgumentError fld_write_data_fix(Rational.([1 0 1]))
end

# for future testing of reading "multi" external files
@testset "write multi" begin
#	file = "../tmp2.fld"
#	file = "tmp.fld"
	data = rand(Float32, 7, 9, 2)
	file1 = "$(basename(file)).part1"
	file2 = "$(basename(file)).part2"
	head = [
		"# AVS field file"
		"ndim=3"
		"dim1=7"
		"dim2=9"
		"dim3=2"
		"nspace=3"
		"veclen=1"
		"data=float"
		"field=uniform"
		"variable 1 file=2 filetype=multi skip=0"
		file1
		file2
	]

	open(file, "w") do fid
		for line in head
			println(fid, line)
		end
	end

	file1 = joinpath(dirname(file), file1)
	file2 = joinpath(dirname(file), file2)
	write(file1, data[:,:,1])
	write(file2, data[:,:,2])

#	run(`cat $file`)
#	run(`ls -l $file1`)
#	run(`op range $file`)

	@test_throws String tmp = fld_read(file) # todo some day...
end

=#

@testset "fld io" begin
	x = 1:5
	file = tempname() * ".fld"
	fld_write(file, x)
	y = fld_read(file)
	@test x == y
end


text:
export newtoncg

"""
newtoncg(f,df,H,x)

Newton-CG method for solving min_x f(x)

"""
function newtoncg(f::Function,df::Function,H::Function,x::Vector;maxIter=20,atol=1e-8,out::Int=0,storeInterm::Bool=false,
	lineSearch::Function=(f,df,fk,dfk,xk,pk)->armijo(f,fk,dfk,xk,pk,maxIter=20),tolCG=1e-2,maxIterCG=30,P=d2f->identity)

    his = zeros(maxIter,6)
    X = (storeInterm) ? zeros(length(x),maxIter) : []
    i = 1; flag = -1; LL = []
    while i<=maxIter
        fk  = f(x)
        dfk = df(x)
        his[i,1:2] = [norm(fk) norm(dfk)]
        if storeInterm; X[:,i] = x; end;

        if(norm(dfk)<atol)
            flag = 0
            his = his[1:i,:]
            break
        end

        # get search direction
        d2f = H(x)
        PC  = P(d2f) # preconditioner
        pk,his[i,6],his[i,5],his[i,4] = cg(d2f,-dfk,out=-1,tol=tolCG,maxIter=maxIterCG,M=PC)
        if his[i,6]==-2 && norm(pk)==0
            pk = -dfk
        end
      # line search
        ak,his[i,3] = lineSearch(f,df,fk,dfk,x,pk)
        if his[i,3]==-1
            flag = -3
            his  = his[1:i,:]
            break
        end
      if out>0
            @printf "iter=%04d\t|f|=%1.2e\t|df|=%1.2e\tcgIter=%d\tLS=%1.4f\n" i his[i,1] his[i,2] his[i,4] ak
        end
        # update
        x  += ak*pk
        i+=1
    end
    i = min(maxIter,i)

    if out>=0
        if flag==-1
            println(@sprintf("newtoncg iterated maxIter (=%d) times but reached only atol of %1.2e instead of tol=%1.2e",i,his[i,2],atol))
        elseif flag==-3
            println(@sprintf("newtoncg stopped at iteration %d due to a line search fail.",i))
        elseif out>1
            println(@sprintf("newtoncg achieved desired atol of %1.2e at iteration %d.",atol,i))
        end
    end
    if storeInterm; X = X[:,1:i]; end
    return x,flag,his,X
end


text:
@generated function addscalar(v::Vec{W,T}, s::T) where {W, T <: IntegerTypesHW}
    typ = "i$(8sizeof(T))"
    vtyp = "<$W x $typ>"
    instrs = String[]
    push!(instrs, "%ie = insertelement $vtyp zeroinitializer, $typ %1, i32 0")
    push!(instrs, "%v = add $vtyp %0, %ie")
    push!(instrs, "ret $vtyp %v")
    quote
        $(Expr(:meta,:inline))
        Vec($LLVMCALL( $(join(instrs,"\n")), NTuple{$W,Core.VecElement{$T}}, Tuple{NTuple{$W,Core.VecElement{$T}},$T}, data(v), s ))
    end
end
@generated function addscalar(v::Vec{W,T}, s::T) where {W, T <: FloatingTypes}
    typ = LLVM_TYPES[T]
    vtyp = "<$W x $typ>"
    instrs = String[]
    push!(instrs, "%ie = insertelement $vtyp zeroinitializer, $typ %1, i32 0")
    push!(instrs, "%v = fadd nsz arcp contract afn reassoc $vtyp %0, %ie")
    push!(instrs, "ret $vtyp %v")
    quote
        $(Expr(:meta,:inline))
        Vec($LLVMCALL( $(join(instrs,"\n")), NTuple{$W,Core.VecElement{$T}}, Tuple{NTuple{$W,Core.VecElement{$T}},$T}, data(v), s ))
    end
end

@generated function mulscalar(v::Vec{W,T}, s::T) where {W, T <: IntegerTypesHW}
    typ = "i$(8sizeof(T))"
    vtyp = "<$W x $typ>"
    instrs = String[]
    push!(instrs, "%ie = insertelement $vtyp $(llvmconst(W, T, 1)), $typ %1, i32 0")
    push!(instrs, "%v = mul $vtyp %0, %ie")
    push!(instrs, "ret $vtyp %v")
    quote
        $(Expr(:meta,:inline))
        Vec($LLVMCALL( $(join(instrs,"\n")), NTuple{$W,Core.VecElement{$T}}, Tuple{NTuple{$W,Core.VecElement{$T}},$T}, data(v), s ))
    end
end
@generated function mulscalar(v::Vec{W,T}, s::T) where {W, T <: FloatingTypes}
    typ = LLVM_TYPES[T]
    vtyp = "<$W x $typ>"
    instrs = String[]
    push!(instrs, "%ie = insertelement $vtyp $(llvmconst(W, T, 1.0)), $typ %1, i32 0")
    push!(instrs, "%v = fmul nsz arcp contract afn reassoc $vtyp %0, %ie")
    push!(instrs, "ret $vtyp %v")
    quote
        $(Expr(:meta,:inline))
        Vec($LLVMCALL( $(join(instrs,"\n")), NTuple{$W,Core.VecElement{$T}}, Tuple{NTuple{$W,Core.VecElement{$T}},$T}, data(v), s ))
    end
end

function scalar_maxmin(W, ::Type{T}, ismax) where {T}
    comp = if T <: Signed
        ismax ? "icmp sgt" : "icmp slt"
    elseif T <: Unsigned
        ismax ? "icmp ugt" : "icmp ult"
    elseif ismax
        "fcmp ogt"
    else
        "fcmp olt"
    end
    if T <: Integer
        typ = "i$(8sizeof(T))"
        basevalue = llvmconst(W, T, ismax ? typemin(T) : typemax(T))
    else
        opzero = ismax ? -Inf : Inf
        basevalue = llvmconst(W, T, repr(reinterpret(UInt64, opzero)))
        if T === Float64
            typ = "double"
            # basevalue = llvmconst(W, T, repr(reinterpret(UInt64, opzero)))
        elseif T === Float32
            typ = "float"
            # basevalue = llvmconst(W, T, repr(reinterpret(UInt32, Float32(opzero))))
        # elseif T === Float16
            # typ = "half"
            # basevalue = llvmconst(W, T, repr(reinterpret(UInt16, Float16(opzero))))
        else
            throw("T === $T not currently supported.")
        end
    end
    vtyp = "<$W x $typ>"
    instrs = String[]
    push!(instrs, "%ie = insertelement $vtyp $(basevalue), $typ %1, i32 0")
    push!(instrs, "%selection = $comp $vtyp %0, %ie")
    push!(instrs, "%v = select <$W x i1> %selection, $vtyp %0, $vtyp %ie")
    push!(instrs, "ret $vtyp %v")
    instrs
end
@generated function maxscalar(v::Vec{W,T}, s::T) where {W, T <: NativeTypes}
    instrs = scalar_maxmin(W, T, true)
    quote
        $(Expr(:meta,:inline))
        Vec($LLVMCALL( $(join(instrs,"\n")), NTuple{$W,Core.VecElement{$T}}, Tuple{NTuple{$W,Core.VecElement{$T}},$T}, data(v), s ))
    end
end
@generated function minscalar(v::Vec{W,T}, s::T) where {W, T <: NativeTypes}
    instrs = scalar_maxmin(W, T, false)
    quote
        $(Expr(:meta,:inline))
        Vec($LLVMCALL( $(join(instrs,"\n")), NTuple{$W,Core.VecElement{$T}}, Tuple{NTuple{$W,Core.VecElement{$T}},$T}, data(v), s ))
    end
end
for (f,op) ∈ [(:addscalar,:(+)), (:mulscalar,:(*)), (:maxscalar,:max), (:minscalar,:min)]
    @eval begin
        @inline $f(v::VecUnroll, s) = VecUnroll(($f(first(getfield(v, :data)), s), Base.tail(getfield(v, :data))...))
        @inline $f(v::Vec{W,T}, s::NativeTypes) where {W, T <: NativeTypes} = $f(v, vconvert(T, s))
        @inline $f(s::NativeTypes, v::AbstractSIMD{W,T}) where {W,T <: NativeTypes} = $f(v, s)
        @inline $f(a, b) = $op(a, b)
    end
end


text:
#=
fft:
- Julia version: 0.6.2
- Author: jerrymei
- Date: 2018-05-05
=#
using FFTW

mutable struct rFFTWrapper{T}
    src::Array{T,2}
    src_hat::Array{Complex{T},2}
    dest_hat::Array{Complex{T},2}
    dest::Array{T,2}
    plan::FFTW.Plan
    inv_plan::FFTW.Plan
    src_hat_calculated::Bool
end

function rFFTWrapper{T}(input::AbstractArray{T,2})
    mx,my=size(input)
    nx,ny=div(mx,2) + 1,my
    src=input
    src_hat=Array{Complex{T}}(nx,ny)
    dest_hat=similar(src_hat)
    dest=similar(input)
    plan=plan_rfft(src)
    inv_plan=plan_irfft(src_hat,mx)
    return rFFTWrapper(src,src_hat,dest_hat,dest,plan,inv_plan,false)
end

function set_src!{T}(fft::rFFTWrapper, input::AbstractArray{T,2})
    if !(input === fft.src)
        assert(size(input) == size(fft.src))
        fft.src=input
        fft.src_hat_calculated=false
    end
end

function set_dest!{T}(fft::rFFTWrapper, output::AbstractArray{T,2})
    assert(size(output) == size(fft.dest))
    fft.dest=output
end

function set_not_calculated(fft::rFFTWrapper)
    fft.src_hat_calculated=false
end

function derivx!(fft::rFFTWrapper)
    mx,my=size(fft.src)
    nx,ny=size(fft.src_hat)
    if !fft.src_hat_calculated
        A_mul_B!(fft.src_hat,fft.plan,fft.src)
        fft.src_hat_calculated=true
    end
    @inbounds for x in 0:nx - 1
        for y in 0:ny - 1
            fft.dest_hat[x + 1,y + 1]=fft.src_hat[x + 1,y + 1] * x * 2pi * im
        end
    end
    A_mul_B!(fft.dest,fft.inv_plan,fft.dest_hat)
end

function derivy!(fft::rFFTWrapper)
    mx,my=size(fft.src)
    nx,ny=size(fft.src_hat)
    if !fft.src_hat_calculated
        A_mul_B!(fft.src_hat,fft.plan,fft.src)
        fft.src_hat_calculated=true
    end
    @inbounds for x in 0:nx - 1
        for y in 0:ny - 1
            fft.dest_hat[x + 1,y + 1]= fft.src_hat[x + 1,y + 1] * (2 * y < my?y:y - my) * 2pi * im
        end
    end
    A_mul_B!(fft.dest,fft.inv_plan,fft.dest_hat)
end

function laplacian!(fft::rFFTWrapper)
    mx,my=size(fft.src)
    nx,ny=size(fft.src_hat)
    if !fft.src_hat_calculated
        A_mul_B!(fft.src_hat,fft.plan,fft.src)
        fft.src_hat_calculated=true
    end
    @inbounds for x in 0:nx - 1
        for y in 0:ny - 1
            fft.dest_hat[x + 1,y + 1]= fft.src_hat[x + 1,y + 1] * (x^2 + (2 * y < my?y:y - my)^2) * (2pi * im)^2
        end
    end
    A_mul_B!(fft.dest,fft.inv_plan,fft.dest_hat)
end

function solve_laplacian!(fft::rFFTWrapper)
    mx,my=size(fft.src)
    nx,ny=size(fft.src_hat)
    if !fft.src_hat_calculated
        A_mul_B!(fft.src_hat,fft.plan,fft.src)
        fft.src_hat_calculated=true
    end
    @inbounds for x in 0:nx - 1
        for y in 1:ny - 1
            fft.dest_hat[x + 1,y + 1]= fft.src_hat[x + 1,y + 1] / (x^2 + (2 * y < my?y:y - my)^2) / (2pi * im)^2
        end
    end
    @inbounds for x in 1:nx - 1
        fft.dest_hat[x + 1,1]=fft.src_hat[x + 1,1] / x^2 / (2pi * im)^2
    end
    A_mul_B!(fft.dest,fft.inv_plan,fft.dest_hat)
end

function derivx!{T}(input::AbstractArray{T,2}, output::AbstractArray{T,2}, fft::rFFTWrapper)
"""Compute y derivative"""
    set_src!(fft,input)
    set_dest!(fft,output)
    derivx!(fft)
end

function derivy!{T}(input::AbstractArray{T,2}, output::AbstractArray{T,2}, fft::rFFTWrapper)
"""Compute y derivative"""
    set_src!(fft,input)
    set_dest!(fft,output)
    derivy!(fft)
end

function laplacian!{T}(input::AbstractArray{T,2}, output::AbstractArray{T,2}, fft::rFFTWrapper)
"""Compute the Laplacian"""
    set_src!(fft,input)
    set_dest!(fft,output)
    laplacian!(fft)
end

function solve_laplacian!{T}(input::AbstractArray{T,2}, output::AbstractArray{T,2}, fft::rFFTWrapper)
"""Solve Poisson equation"""
    set_src!(fft,input)
    set_dest!(fft,output)
    solve_laplacian!(fft)
end

max_stars_repo_name: JerryLingjieMei/NonlinearPDE

max_issues_repo_name: JerryLingjieMei/NonlinearPDE

max_forks_repo_name: JerryLingjieMei/NonlinearPDE


text:
module Internals

## Internal request utilities ##
function canvas_headers end
function process_params end
function process_key end
function process_val end
function request end
function paged_request end

## Internal CanvasObject utilities ##
function id end
function pid end
function json2canvas end
function adjust_fieldtype end
function convert′ end

## Endpoint internals ##
function upload_file end

end # module


text:
function mountaincar_rtaa_main(; params = true_params)::Vector{Int64}
    rtaa_steps = []
    model = MountainCar(0.0)
    planner = MountainCarRTAAPlanner(model, 1000, params)
    generateHeuristic!(planner)
    for rock_c in range_of_values
        println("Rock_c is ", rock_c)
        mountaincar = MountainCar(rock_c)
        agent = MountainCarRTAAAgent(mountaincar, planner)
        n_steps = run(agent, max_steps = 1e4)
        push!(rtaa_steps, n_steps)
        clearResiduals!(planner)
    end
    rtaa_steps
end

function mountaincar_true_main()
    true_steps = []
    for rock_c in range_of_values
        println("Rock_c is ", rock_c)
        mountaincar = MountainCar(rock_c)
        model = MountainCar(rock_c)
        planner = MountainCarRTAAPlanner(model, 1000, true_params)
        generateHeuristic!(planner)
        agent = MountainCarRTAAAgent(mountaincar, planner)
        n_steps = run(agent, max_steps = 1e4)
        push!(true_steps, n_steps)
    end
    true_steps
end

function mountaincar_rtaa_grid_search()
    data = fill(0, (length(theta1_values), length(theta2_values), length(range_of_values)))
    for idx1 = 1:length(theta1_values)
        theta1 = theta1_values[idx1]
        for idx2 = 1:length(theta2_values)
            theta2 = theta2_values[idx2]
            params = MountainCarParameters(theta1, theta2)
            rtaa_steps = mountaincar_rtaa_main(params = params)
            data[idx1, idx2, :] = rtaa_steps
        end
    end
    save(
        gridsearch_data_path,
        "data",
        data,
        "theta1",
        collect(theta1_values),
        "theta2",
        collect(theta2_values),
    )
end


text:
losses = [[10, 0],
          [1, 50],
          [0, 200]]

num_actions = length(losses)

function expected_loss_of_action(prob_spam, action)
    #TODO: Return expected loss over a Bernoulli random variable
    #      with mean prob_spam.
    #      Losses are given by the table above.
    losses = [[10, 0],
              [1, 50],
              [0, 200]]
    loss = losses[action]
    length_prob_spam = length(prob_spam)

    non_prob = 1 .- prob_spam
    expect_loss = [prob_spam[i]*loss[1]+non_prob[i]*loss[2] for i in 1:length_prob_spam]
    return expect_loss
end

prob_range = range(0., stop=1., length=500)

using Plots
for action in 1:num_actions
  display(plot!(prob_range, expected_loss_of_action(prob_range, action)))
end

function optimal_action(prob_spam)
    #TODO: return best action given the probability of spam.
    # Hint: Julia's findmin function might be helpful.

    losses = [[10, 0],[1, 50],[0, 200]]
    non_prob_spam = 1-prob_spam
    expected_loss = []
    for i in 1:3
        loss = losses[i]
        push!(expected_loss, loss[1]*prob_spam+loss[2]*non_prob_spam)
    end
    return findmin(expected_loss)
end

prob_range = range(0, stop=1., length=500)
optimal_losses = []
optimal_actions = []
for p in prob_range
    # TODO:  Compute the optimal action and its expected loss for
    # probability of spam given by p.
    res = optimal_action(p)
    push!(optimal_losses,res[1])
    push!(optimal_actions,res[2])

end
optimal_actions
plot(prob_range, optimal_losses, linecolor=optimal_actions)


text:
using LinearAlgebra

function MPM(Y,K, m=0)
    # Matrix pencil method for line spectral estimation from measurements.
    # Can handle multiple snapshots by constructing a low rank Hankel matrix
    # from the data.
    #
    # Inputs:
    # Y: Data matrix
    # K: number of sinusoids
    # m: Hankel parameter of the Hankel matrix of the data
    #
    # Parameters:
    # N = size(Y,1): number of measurements per snapshot
    # L = size(Y,2): number of snapshots
    #
    # Output:
    # f_hat: estimate of the frequencies of the K sinusoids in [rad]

(N,L) = size(Y)

# If not given from the input, use standard pencil parameter m
if m==0
m = floor(Int,N/2)
end

# Construct horizontally stacked Toeplitz matrix
Y_toep = Array{Complex{Float64},2}(undef,N-m,(m+1)*L)

for j = 1:L, i = 1:m+1
    Y_toep[:,i+(j-1)*(m+1)] = Y[i:N-m-1+i,j]
end

clmnsY0 = setdiff(1:(m+1)*L,(1:L).*(m+1))
clmnsY1 = clmnsY0.+1

Y0 = Y_toep[:,clmnsY0]
Y1 = Y_toep[:,clmnsY1]

F = svd(Y0)
p = sortperm(F.S,rev=true)[1:K]

f = eigvals(Diagonal(ones(K)./F.S[p])*F.U[:,p]'*Y1*F.V[:,p])
f = (angle.(f)/(2*pi) .+ 1).%1
end