text:
using DDD
using Test, LinearAlgebra, StaticArrays
cd(@__DIR__)

@testset "Calculating sigma_tilde" begin
    numNode = 2
    numSeg = numNode - 1
    len = numNode + 1
    xrange = range(-300, 300, length = len)
    yrange = range(-300, 300, length = len)

    X = ones(length(yrange)) .* xrange'
    Y = ones(length(xrange))' .* yrange
    Z = zeros(length(xrange))' .* zeros(len)
    points = [X[:]'; Y[:]'; Z[:]']

    l = Float64[0; 0; 1]
    b = Float64[1; 0; 0]
    n = b × l
    a = 5 * norm(b)

    links = zeros(Int, 2, numSeg)
    slipPlane = zeros(3, numSeg)
    bVec = zeros(3, numSeg)
    coord = [zeros(len)'; zeros(len)'; xrange']
    label = zeros(nodeTypeDln, len)
    nodeVel = similar(coord)
    nodeForce = similar(coord)
    for i in 1:numSeg
        links[:, i] .= (i, i + 1)
        slipPlane[:, i] = n
        bVec[:, i] = b
    end

    matParams = MaterialParameters(; crystalStruct = BCC(), μ = 1.0, μMag = 1.0, ν = 0.28)
    dlnParams = DislocationParameters(;
        coreRad = a,
        coreRadMag = 1.0,
        minSegLen = a + 2,
        maxSegLen = a + 3,
        minArea = a + 1,
        maxArea = a + 2,
        dragCoeffs = (edge = 1.0, screw = 1.0, climb = 1.0, line = 1.0),
        mobility = mobBCC(),
    )
    network = DislocationNetwork(;
        links = links,
        slipPlane = slipPlane,
        bVec = bVec,
        coord = coord,
        label = label,
        nodeVel = nodeVel,
        nodeForce = nodeForce,
    )
    makeConnect!(network)
    getSegmentIdx!(network)

    σTilde = calc_σTilde(points, dlnParams, matParams, network)'

    test_σTilde = [
        0.000248130106066 0.000035465440391 -0.000011334371408 0.000035426063059 -0.000045958542661 -0.000059335461657
        0 0 0 -0.000260417766102 0.000030180624794 0
        -0.000248130106066 -0.000035465440391 0.000011334371408 0.000035426063059 -0.000045958542661 0.000059335461657
        0.000260598566573 0.000390644731208 0.000015657608042 0 -0.000207853177578 0
        0 0 0 0 -0.001871493087898 0
        -0.000260598566573 -0.000390644731208 -0.000015657608042 0 -0.000207853177578 0
        0.000248130106066 0.000035465440391 -0.000011334371408 -0.000035426063059 -0.000045958542661 0.000059335461657
        0 0 0 0.000260417766102 0.000030180624794 0
        -0.000248130106066 -0.000035465440391 0.000011334371408 -0.000035426063059 -0.000045958542661 -0.000059335461657
    ]

    @test isapprox(σTilde[:, 1:4], test_σTilde[:, 1:4])
    @test isapprox(σTilde[:, 6], test_σTilde[:, 5])
    @test isapprox(σTilde[:, 5], test_σTilde[:, 6])

    σ = zeros(6, size(points, 2))
    calc_σTilde!(σ, points, dlnParams, matParams, network)
    @test isequal(σ', σTilde)

    ##

    numNode = 2
    numSeg = numNode - 1
    len = numNode + 1
    xrange = range(-300, 300, length = len)
    yrange = range(-300, 300, length = len)

    X = ones(length(yrange)) .* xrange'
    Y = ones(length(xrange))' .* yrange
    Z = zeros(length(xrange))' .* zeros(len)
    points = [X[:]'; Y[:]'; Z[:]']

    l = Float64[0; 0; 1]
    b = Float64[0; 0; 1]
    n = b × l
    a = 5 * norm(b)

    links = zeros(Int, 2, numSeg)
    slipPlane = zeros(3, numSeg)
    bVec = zeros(3, numSeg)
    coord = [zeros(len)'; zeros(len)'; xrange']
    label = zeros(nodeTypeDln, len)
    nodeVel = similar(coord)
    nodeForce = similar(coord)
    for i in 1:numSeg
        links[:, i] .= (i, i + 1)
        slipPlane[:, i] = n
        bVec[:, i] = b
    end

    matParams = MaterialParameters(; crystalStruct = BCC(), μ = 1.0, μMag = 1.0, ν = 0.28)
    dlnParams = DislocationParameters(;
        coreRad = a,
        coreRadMag = 1.0,
        minSegLen = a + 2,
        maxSegLen = a + 3,
        minArea = a + 1,
        maxArea = a + 2,
        dragCoeffs = (edge = 1.0, screw = 1.0, climb = 1.0, line = 1.0),
        mobility = mobBCC(),
    )
    network = DislocationNetwork(;
        links = links,
        slipPlane = slipPlane,
        bVec = bVec,
        coord = coord,
        label = label,
        nodeVel = nodeVel,
        nodeForce = nodeForce,
    )
    makeConnect!(network)
    getSegmentIdx!(network)

    σTilde = calc_σTilde(points, dlnParams, matParams, network)'

    test_σTilde =
        1.0e-03 * [
            0 0 0 0 -0.076573455482457 0.076573455482457
            0 0 0 0 -0.187565879762850 0
            0 0 0 0 -0.076573455482457 -0.076573455482457
            0 0 0 0 0 0.187565879762850
            0 0 0 0 0 0
            0 0 0 0 0 -0.187565879762850
            0 0 0 0 0.076573455482457 0.076573455482457
            0 0 0 0 0.187565879762850 0
            0 0 0 0 0.076573455482457 -0.076573455482457
        ]

    @test isapprox(σTilde[:, 1:4], test_σTilde[:, 1:4])
    @test isapprox(σTilde[:, 6], test_σTilde[:, 5])
    @test isapprox(σTilde[:, 5], test_σTilde[:, 6])

    σ = zeros(6, size(points, 2))
    calc_σTilde!(σ, points, dlnParams, matParams, network)
    @test isequal(σ', σTilde)

    ##

    numNode = 50
    numSeg = numNode - 1
    len = numNode + 1
    xrange = range(-300, 300, length = len)
    yrange = range(-300, 300, length = len)

    X = ones(length(yrange)) .* xrange'
    Y = ones(length(xrange))' .* yrange
    Z = zeros(length(xrange))' .* zeros(len)
    points = [X[:]'; Y[:]'; Z[:]']

    l = Float64[0; 0; 1]
    b = Float64[0; 0; 1]
    n = b × l
    a = 5 * norm(b)

    links = zeros(Int, 2, numSeg)
    slipPlane = zeros(3, numSeg)
    bVec = zeros(3, numSeg)
    coord = [zeros(len)'; zeros(len)'; xrange']
    label = zeros(nodeTypeDln, len)
    nodeVel = similar(coord)
    nodeForce = similar(coord)
    for i in 1:numSeg
        links[:, i] .= (i, i + 1)
        slipPlane[:, i] = n
        bVec[:, i] = b
    end

    matParams = MaterialParameters(; crystalStruct = BCC(), μ = 1.0, μMag = 1.0, ν = 0.28)
    dlnParams = DislocationParameters(;
        coreRad = a,
        coreRadMag = 1.0,
        minSegLen = a + 2,
        maxSegLen = a + 3,
        minArea = a + 1,
        maxArea = a + 2,
        dragCoeffs = (edge = 1.0, screw = 1.0, climb = 1.0, line = 1.0),
        mobility = mobBCC(),
    )
    network = DislocationNetwork(;
        links = links,
        slipPlane = slipPlane,
        bVec = bVec,
        coord = coord,
        label = label,
        nodeVel = nodeVel,
        nodeForce = nodeForce,
    )
    makeConnect!(network)
    getSegmentIdx!(network)

    stress = reshape(calc_σTilde(points, dlnParams, matParams, network), 6, len, :)

    σ = zeros(6, size(points, 2))
    calc_σTilde!(σ, points, dlnParams, matParams, network)
    σ = reshape(σ, 6, len, :)

    @test isequal(σ, stress)

    ##
    l = Float64[0; 0; 1]
    b = Float64[1; 0; 0]
    n = b × l
    a = 5 * norm(b)

    links = zeros(Int, 2, numSeg)
    slipPlane = zeros(3, numSeg)
    bVec = zeros(3, numSeg)
    coord = [zeros(len)'; zeros(len)'; xrange']
    label = zeros(nodeTypeDln, len)
    nodeVel = similar(coord)
    nodeForce = similar(coord)
    for i in 1:numSeg
        links[:, i] .= (i, i + 1)
        slipPlane[:, i] = n
        bVec[:, i] = b
    end

    matParams = MaterialParameters(; crystalStruct = BCC(), μ = 1.0, μMag = 1.0, ν = 0.28)
    dlnParams = DislocationParameters(;
        coreRad = a,
        coreRadMag = 1.0,
        minSegLen = a + 2,
        maxSegLen = a + 3,
        minArea = a + 1,
        maxArea = a + 2,
        dragCoeffs = (edge = 1.0, screw = 1.0, climb = 1.0, line = 1.0),
        mobility = mobBCC(),
    )
    network = DislocationNetwork(;
        links = links,
        slipPlane = slipPlane,
        bVec = bVec,
        coord = coord,
        label = label,
        nodeVel = nodeVel,
        nodeForce = nodeForce,
    )
    makeConnect!(network)
    getSegmentIdx!(network)

    stress = reshape(calc_σTilde(points, dlnParams, matParams, network), 6, len, :)

    σ = zeros(6, size(points, 2))
    calc_σTilde!(σ, points, dlnParams, matParams, network)
    σ = reshape(σ, 6, len, :)

    @test isequal(σ, stress)
end

@testset "Calculate uTilde" begin
    dlnParams = DislocationParameters(; mobility = mobBCC())
    matParams = MaterialParameters(; crystalStruct = BCC(), ν = 0.28)
    femParams = FEMParameters(;
        type = DispatchRegularCuboidMesh(),
        order = LinearElement(),
        model = CantileverLoad(),
        dx = 1013.0,
        dy = 1987.0,
        dz = 2999.0,
        mx = 3,
        my = 5,
        mz = 7,
    )
    intParams = IntegrationParameters(; method = AdaptiveEulerTrapezoid())
    slipSystem = SlipSystem(;
        crystalStruct = BCC(),
        slipPlane = Float64[-1; 1; 0],
        bVec = Float64[1; 1; 1],
    )
    dx, dy, dz = femParams.dx, femParams.dy, femParams.dz
    segLen = (dx + dy + dz) / 30
    prismLoop = DislocationLoop(;
        loopType = loopPrism(),    # Prismatic loop, all segments are edge segments.
        numSides = 8,   # 8-sided loop.
        nodeSide = 1,   # One node per side, if 1 nodes will be in the corners.
        numLoops = 1,   # Number of loops of this type to generate when making a network.
        segLen = segLen * SVector{8}(ones(8)),  # Length of each segment between nodes, equal to the number of nodes.
        slipSystemIdx = 1, # Slip System (assuming slip systems are stored in a file, this is the index).
        slipSystem = slipSystem,  # Slip system of the segments.
        label = SVector{8, nodeTypeDln}(1, 1, 1, 1, 1, 1, 1, 1),    # Node labels, has to be equal to the number of nodes.
        buffer = 0,   # Buffer to increase the dislocation spread.
        range = SMatrix{3, 2, Float64}(dx / 2, dy / 2, dz / 2, dx / 2, dy / 2, dz / 2),  # Distribution range
        dist = Zeros(),  # Loop distribution.
    )
    shearLoop = DislocationLoop(;
        loopType = loopShear(),    # Prismatic loop, all segments are edge segments.
        numSides = 8,   # 8-sided loop.
        nodeSide = 1,   # One node per side, if 1 nodes will be in the corners.
        numLoops = 1,   # Number of loops of this type to generate when making a network.
        segLen = segLen * SVector{8}(ones(8)),  # Length of each segment between nodes, equal to the number of nodes.
        slipSystemIdx = 1, # Slip System (assuming slip systems are stored in a file, this is the index).
        slipSystem = slipSystem,  # Slip plane of the segments.
        label = SVector{8, nodeTypeDln}(1, 1, 1, 1, 1, 1, 1, 1),    # Node labels, has to be equal to the number of nodes.
        buffer = 0,   # Buffer to increase the dislocation spread.
        range = SMatrix{3, 2, Float64}(dx / 2, dy / 2, dz / 2, dx / 2, dy / 2, dz / 2),  # Distribution range
        dist = Zeros(),  # Loop distribution.
    )
    network = DislocationNetwork([prismLoop, shearLoop])
    regularCuboidMesh = buildMesh(matParams, femParams)
    cantileverBC, forceDisplacement = Boundaries(femParams, regularCuboidMesh)

    uGammaTest = [
        1
        5
        9
        13
        17
        21
        25
        29
        33
        37
        41
        45
        49
        53
        57
        61
        65
        69
        73
        77
        81
        85
        89
        93
        97
        101
        105
        109
        113
        117
        121
        125
        129
        133
        137
        141
        145
        149
        153
        157
        161
        165
        169
        173
        177
        181
        185
        189
        32
        64
        96
        128
        160
        192
    ]

    @test [32, 64, 96, 128, 160] == setdiff(uGammaTest, cantileverBC.uGammaDln)

    uTilde = calc_uTilde(regularCuboidMesh, cantileverBC, matParams, network)

    ux = uTilde[1:3:end]
    uy = uTilde[2:3:end]
    uz = uTilde[3:3:end]
    uxTest = [
        0.002902351345488887,
        -0.0008253133519966701,
        -0.00021650827733668473,
        0.0006009438041774093,
        0.0029629460677653143,
        0.0020545681539812758,
        0.0005322552345154338,
        -0.0005404119399346637,
        -0.0005134326512485413,
        -0.0004903899090375116,
        -6.880892240796615e-6,
        0.0004742726937890663,
        0.005397126562750102,
        0.009359641926982321,
        0.01170673090883372,
        0.007291725493998129,
        0.0018495485834692132,
        4.236557903560154e-5,
        -0.0017021045567632155,
        -0.002843394458582297,
        -0.0024107193135325805,
        0.0005391062331057765,
        0.0018200619637690372,
        0.0012347172382416167,
        0.006602657554255273,
        0.00477962525180953,
        0.00038355589475357044,
        -0.00187535023767425,
        0.015312108727395595,
        0.012116022126695285,
        -0.002632912972606305,
        -0.0057063039169111505,
        0.025588846219305697,
        0.016235303324262124,
        -0.015472906634952174,
        -0.010385046987745625,
        0.012677799867875676,
        -0.004826705167431574,
        -0.023681161774827538,
        -0.0029539571317446283,
        -8.098032109580808e-5,
        -0.005475618765255888,
        -0.002317412841145468,
        0.0021198441533540995,
        -0.0009284282909694205,
        -0.0014265217871674203,
        -4.468666199580377e-5,
        0.0012357735270278093,
        0.0008253133519967209,
        -0.0029023513454888464,
        0.0005404119399346326,
        -0.000532255234515522,
        -0.002054568153981224,
        -0.0029629460677652792,
    ]
    uyTest = [
        0.0067363815836566756,
        0.0014278511601506676,
        -0.0005688368686884211,
        -0.004853663242082534,
        0.0054292511205919925,
        0.0032679235999618647,
        0.001853515275120392,
        0.0015749470920581972,
        -0.0009098880504148402,
        -0.0011964949197497293,
        -0.002159488471687042,
        -0.003753398090388903,
        0.012104398173005865,
        0.020456327309446968,
        0.02531347315436703,
        0.015967697461884905,
        0.004258310757917565,
        0.000136935016830809,
        0.0021021084169944824,
        0.0019008757293132049,
        -0.002757046866739879,
        -0.010749158149593692,
        -0.011816314972861354,
        -0.008061677900843734,
        0.011133590727917276,
        0.006594615229743162,
        0.0033694655216900955,
        0.0029191581465615618,
        0.024101063034655975,
        0.01567205085188615,
        0.0076048333077830595,
        0.005486079417588696,
        0.03910246222571273,
        0.03110200505064155,
        0.01434813388039752,
        0.004463311985544097,
        0.0210078769595474,
        0.009726316969711102,
        0.001130961200306315,
        -0.009841537292647575,
        0.0015061696382455902,
        -0.0023875791696053553,
        -0.006296297689713076,
        -0.01224629630206807,
        -0.001057544898685026,
        -0.0017861516064492313,
        -0.0037456489215073716,
        -0.0069507076726078865,
        -0.0014278511601506416,
        -0.006736381583656756,
        -0.0015749470920582325,
        -0.0018535152751203109,
        -0.0032679235999618314,
        -0.005429251120592,
    ]
    uzTest = [
        0.007979193931819731,
        -0.0015951911425738472,
        0.0012361385228138791,
        -0.005673325438189695,
        0.008293744049109469,
        0.005993888724568718,
        0.0020111921045511444,
        -0.0008199479417571705,
        0.0011040570890579328,
        -0.0005229167518047081,
        -0.00342451549514798,
        -0.005491419587598859,
        0.011382337565374375,
        0.013338843818961211,
        0.008721215722078795,
        0.0014455787776552964,
        0.0007183293235734062,
        0.0014301694225654616,
        -0.002556389674620718,
        -0.0030809269872969602,
        -0.002751202465258223,
        -0.0049295376454319825,
        -0.007995531449452647,
        -0.007572525652961138,
        0.014145193432864162,
        0.010859230173295977,
        0.0022227278814157146,
        -0.0025468808611998463,
        0.021927866951300702,
        0.019333802591945715,
        0.0005872825065838683,
        -0.005515269881330632,
        0.018591615331351674,
        0.020285444264391116,
        -0.0017472723394227748,
        -0.006145908832257655,
        0.002491294207535403,
        0.007236926045437346,
        0.0001805379885031523,
        -0.006972815701468629,
        0.0027768377925153826,
        0.00465361358956461,
        -0.005817365214834877,
        -0.011041518738134993,
        0.0022959756942656615,
        0.0005704038917998867,
        -0.00519688720088853,
        -0.008537520866297306,
        0.0015951911425738617,
        -0.007979193931819718,
        0.0008199479417571756,
        -0.0020111921045511297,
        -0.005993888724568653,
        -0.008293744049109446,
    ]
    @test isapprox(ux, uxTest, rtol = 5e-5)
    @test isapprox(uy, uyTest, rtol = 5e-5)
    @test isapprox(uz, uzTest, rtol = 5e-5)

    calc_uTilde!(forceDisplacement, regularCuboidMesh, cantileverBC, matParams, network)

    uGammaDofs = cantileverBC.uDofsDln
    @test isapprox(uTilde, forceDisplacement.uTilde[uGammaDofs])

    uTildeTest = [
        0.002902351345488887,
        -0.0017021045567632155,
        0.012677799867875676,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0014278511601506676,
        0.0019008757293132049,
        0.009726316969711102,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        -0.0015951911425738472,
        -0.0030809269872969602,
        0.007236926045437346,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        -0.00021650827733668473,
        -0.0024107193135325805,
        -0.023681161774827538,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        -0.0005688368686884211,
        -0.002757046866739879,
        0.001130961200306315,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0067363815836566756,
        0.0021021084169944824,
        0.0210078769595474,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0029629460677653143,
        0.0018200619637690372,
        -8.098032109580808e-5,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        -0.0005134326512485413,
        0.00038355589475357044,
        -0.0009284282909694205,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        -0.0009098880504148402,
        0.0033694655216900955,
        -0.001057544898685026,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0011040570890579328,
        0.0022227278814157146,
        0.0022959756942656615,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        -0.0004903899090375116,
        -0.00187535023767425,
        -0.0014265217871674203,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0005322552345154338,
        0.006602657554255273,
        -0.002317412841145468,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0054292511205919925,
        -0.011816314972861354,
        0.0015061696382455902,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        -0.0011964949197497293,
        0.0029191581465615618,
        -0.0017861516064492313,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        -0.0005229167518047081,
        -0.0025468808611998463,
        0.0005704038917998867,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        -6.880892240796615e-6,
        0.015312108727395595,
        -4.468666199580377e-5,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        -0.002159488471687042,
        0.024101063034655975,
        -0.0037456489215073716,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.001853515275120392,
        0.011133590727917276,
        -0.006296297689713076,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.008293744049109469,
        -0.007995531449452647,
        0.0027768377925153826,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        -0.00342451549514798,
        0.021927866951300702,
        -0.00519688720088853,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0004742726937890663,
        0.012116022126695285,
        0.0012357735270278093,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        -0.003753398090388903,
        0.01567205085188615,
        -0.0069507076726078865,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        -0.005491419587598859,
        0.019333802591945715,
        -0.008537520866297306,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0020111921045511444,
        0.014145193432864162,
        -0.005817365214834877,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0020545681539812758,
        0.0012347172382416167,
        -0.005475618765255888,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.005397126562750102,
        -0.002632912972606305,
        0.0008253133519967209,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.012104398173005865,
        0.0076048333077830595,
        -0.0014278511601506416,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.011382337565374375,
        0.0005872825065838683,
        0.0015951911425738617,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.009359641926982321,
        -0.0057063039169111505,
        -0.0029023513454888464,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        -0.0005404119399346637,
        0.00477962525180953,
        0.0021198441533540995,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0032679235999618647,
        -0.008061677900843734,
        -0.0023875791696053553,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.020456327309446968,
        0.005486079417588696,
        -0.006736381583656756,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.013338843818961211,
        -0.005515269881330632,
        -0.007979193931819718,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.01170673090883372,
        0.025588846219305697,
        0.0005404119399346326,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.02531347315436703,
        0.03910246222571273,
        -0.0015749470920582325,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0015749470920581972,
        0.006594615229743162,
        -0.01224629630206807,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.005993888724568718,
        -0.007572525652961138,
        0.00465361358956461,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.008721215722078795,
        0.018591615331351674,
        0.0008199479417571756,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.007291725493998129,
        0.016235303324262124,
        -0.000532255234515522,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.015967697461884905,
        0.03110200505064155,
        -0.0018535152751203109,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0014455787776552964,
        0.020285444264391116,
        -0.0020111921045511297,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        -0.0008199479417571705,
        0.010859230173295977,
        -0.011041518738134993,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.007979193931819731,
        -0.002556389674620718,
        0.002491294207535403,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0018495485834692132,
        -0.015472906634952174,
        -0.002054568153981224,
        0.0012361385228138791,
        -0.002751202465258223,
        0.0001805379885031523,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0007183293235734062,
        -0.0017472723394227748,
        -0.005993888724568653,
        0.0006009438041774093,
        0.0005391062331057765,
        -0.0029539571317446283,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        4.236557903560154e-5,
        -0.010385046987745625,
        -0.0029629460677652792,
        -0.004853663242082534,
        -0.010749158149593692,
        -0.009841537292647575,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.000136935016830809,
        0.004463311985544097,
        -0.005429251120592,
        -0.005673325438189695,
        -0.0049295376454319825,
        -0.006972815701468629,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0014301694225654616,
        -0.006145908832257655,
        -0.008293744049109446,
        -0.0008253133519966701,
        -0.002843394458582297,
        -0.004826705167431574,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.004258310757917565,
        0.01434813388039752,
        -0.0032679235999618314,
    ]
    @test isapprox(uTildeTest, forceDisplacement.uTilde, rtol = 5e-5)
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
using Random, Test
using TupleTools

using Base: tail, front

n = 10

p = randperm(n)
ip = invperm(p)

t = (p...,)

@test @inferred(TupleTools.argtail2(1, 2, 3, 4)) == (3, 4)
@test @inferred(TupleTools.tail2(t)) == t[3:n]
@test @inferred(TupleTools.unsafe_tail(t)) == t[2:n]
@test @inferred(TupleTools.unsafe_front(t)) == t[1:n-1]
@test @inferred(TupleTools.unsafe_tail(())) == ()
@test @inferred(TupleTools.unsafe_front(())) == ()
@test @inferred(TupleTools.vcat()) == ()
@test @inferred(TupleTools.vcat((1, 2))) == (1, 2)
@test @inferred(TupleTools.getindices(t, ())) == ()
@test @inferred(TupleTools.getindices(t, (1,2,3))) == t[1:3]

@test @inferred(TupleTools.deleteat((1, 2), (1, ))) == (2, )
for i = 1:n
    @test @inferred(TupleTools.deleteat(t, i)) == (deleteat!(copy(p), i)...,)
    @test @inferred(TupleTools.insertat(t, i, (1,2,3))) == (vcat(p[1:i-1], [1,2,3], p[i+1:n])...,)
end
@test @inferred(TupleTools.deleteat((1,2,3,4,5,6), (3,1,5))) == (2,4,6)
@test_throws BoundsError TupleTools.deleteat(t, 0)
@test_throws BoundsError TupleTools.deleteat(t, n+1)
@test_throws BoundsError TupleTools.insertat(t, 0, (1,2,3))
@test_throws BoundsError TupleTools.insertat(t, n+1, (1,2,3))
for i = 0:n
    @test @inferred(TupleTools.insertafter(t, i, (1,2,3))) == (vcat(p[1:i], [1,2,3], p[i+1:n])...,)
end
@test_throws BoundsError TupleTools.insertafter(t, -1, (1,2,3))
@test_throws BoundsError TupleTools.insertafter(t, n+1, (1,2,3))
@test @inferred(TupleTools.vcat((1,2,3),4,(5,),(),(6,7,8))) == (1,2,3,4,5,6,7,8)
@test @inferred(TupleTools.flatten((1,2,3),4,(5,),(),(6,7,8))) == (1,2,3,4,5,6,7,8)
@test @inferred(TupleTools.flatten((1,(2,3)),4,(5,),(),((6,),(7,(8,))))) == (1,2,3,4,5,6,7,8)

@test @inferred(TupleTools.sum(())) == 0
@test @inferred(TupleTools.sum((1, ))) == 1
@test @inferred(TupleTools.sum(t)) == sum(t)
@test @inferred(TupleTools.cumsum(())) == ()
@test @inferred(TupleTools.cumsum((1, ))) == (1, )
@test @inferred(TupleTools.cumsum(t)) == (cumsum(p)...,)
@test @inferred(TupleTools.prod(())) == 1
@test @inferred(TupleTools.prod((2, ))) == 2
@test @inferred(TupleTools.prod(t)) == prod(t)
@test @inferred(TupleTools.cumprod(())) == ()
@test @inferred(TupleTools.cumprod((1, ))) == (1, )
@test @inferred(TupleTools.cumprod(t)) == (cumprod(p)...,)

for a in (t, (1, ))
    @test @inferred(TupleTools.findmin(a)) == findmin(a)
    @test @inferred(TupleTools.findmax(a)) == findmax(a)
    @test @inferred(TupleTools.minimum(a)) == minimum(a)
    @test @inferred(TupleTools.maximum(a)) == maximum(a)
    @test @inferred(TupleTools.argmin(a)) == argmin(a)
    @test @inferred(TupleTools.argmax(a)) == argmax(a)
end

@test @inferred(TupleTools.sort((1, ))) == (1, )
@test @inferred(TupleTools.sort(())) == ()
@test @inferred(TupleTools.sort(t; rev = true)) == (sort(p; rev = true)...,)
@test @inferred(TupleTools.sort(t; rev = false)) == (sort(p; rev = false)..., )
@test @inferred(TupleTools.sortperm(t)) == (sortperm(p)...,)
@test @inferred(TupleTools.sortperm(t; rev=true)) == (sortperm(p; rev=true)..., )
@test @inferred(TupleTools.invperm(t)) == (ip...,)
@test @inferred(TupleTools.isperm(t)) == true
@test @inferred(TupleTools.isperm((1,2,1))) == false
@test @inferred(TupleTools.permute(t, t)) == (p[p]...,)

@test @inferred(TupleTools.vcat()) == ()
@test @inferred(TupleTools.diff(())) == ()
@test @inferred(TupleTools.diff((1, ))) == ()
@test @inferred(TupleTools.diff((1, 2, 3))) == (1, 1)

ext: jl

max_stars_repo_name: Jutho/TupleTools

max_issues_repo_name: Jutho/TupleTools

max_forks_repo_name: Jutho/TupleTools


text:
@testset "Random sampling on unit, n-dimensional sphere" begin

    uv1 = BlackBoxOptim.sample_unit_hypersphere(1)
    @test size(uv1, 1) == 1

    uv2 = BlackBoxOptim.sample_unit_hypersphere(2)
    @test size(uv2, 1) == 2

    uv3 = BlackBoxOptim.sample_unit_hypersphere(3, 4)
    @test size(uv3, 1) == 3
    @test size(uv3, 2) == 4

    for i in 1:100
        n = rand(1:100)
        num = rand(2:10)
        u = BlackBoxOptim.sample_unit_hypersphere(n, num)
        @test size(u, 1) == n
        @test size(u, 2) == num
        for j in 1:num
            @test isapprox( norm(u[:,j]), 1.0 )
        end
    end

end

@testset "Random direction generator" begin

    rdg1 = BlackBoxOptim.RandomDirectionGen(2, 3)
    ds1 = BlackBoxOptim.directions_for_k(rdg1, 1)
    @test size(ds1) == (2, 3)

    rdg2 = BlackBoxOptim.RandomDirectionGen(10, 17)
    ds2 = BlackBoxOptim.directions_for_k(rdg2, 1)
    @test size(ds2) == (10, 17)

end

@testset "Mirrored random direction generator" begin

    mrdg1 = BlackBoxOptim.MirroredRandomDirectionGen(2, 4)
    ds1 = BlackBoxOptim.directions_for_k(mrdg1, 1)
    @test size(ds1) == (2, 4)
    @test ds1[:,3] == -ds1[:,1]
    @test ds1[:,4] == -ds1[:,2]

    mrdg2 = BlackBoxOptim.MirroredRandomDirectionGen(10, 6)
    ds2 = BlackBoxOptim.directions_for_k(mrdg2, 1)
    @test size(ds2) == (10, 6)
    @test ds2[:,4] == -ds2[:,1]
    @test ds2[:,5] == -ds2[:,2]
    @test ds2[:,6] == -ds2[:,3]

    # Must be even number of directions
    @test_throws ArgumentError BlackBoxOptim.MirroredRandomDirectionGen(10, 1)
    @test_throws ArgumentError BlackBoxOptim.MirroredRandomDirectionGen(10, 3)
    @test_throws ArgumentError BlackBoxOptim.MirroredRandomDirectionGen(10, 7)

end


text:
using MomentBasedEstimators
using Test

facts("Testing basic interface") do
    context("Test from vignette for the R gmm package. Automatic Differentiation") do
        include("normal_dist.jl")
        cft = [3.843, 2.067]
        cfe = coef(step_1)
        sdt = [0.06527666675890574,0.04672629429325811]
        sde = stderror(step_1, TwoStepGMM())
        for j = 1:2
            @test cfe[j] ≈ cft[j] atol = 0.001
            @test sde[j] ≈ sdt[j] atol = 0.00001
        end
        Je, pe = MomentBasedEstimators.J_test(gmm_iid_mgr)
        ## This is the objective value, which is the J-test
        ## for other softwares
        ov     = 1.4398836656920428
        Jt, pt = (1.4378658264483137,0.23048500853597673)
        @test Je ≈ Jt atol = 0.0000001
        @test pe ≈ pt atol = 0.0000001
        @test objval(gmm_iid_mgr) ≈ ov atol = 0.0001
        @test objval(gmm_iid_mgr) < ov
    end

    context("Example 13.5 from Greene (2012) -- verified with Stata") do
        include("gamma_dist.jl")
        cf_stata = [3.358432, .1244622]
        cfe = coef(two_step)
        for j = 1:2
            @test cfe[j] ≈ cf_stata[j] atol=1e-3
        end
        Je, pe = MomentBasedEstimators.J_test(two_step)
        ## The following test is conditional on version.
        ## What changed is the tol for pinv between 0.3 and 0.4
        ##
        if VERSION < v"0.4-"
            Jt, pt = (1.7433378725860427,0.41825292894063326)
        else
            Jt, pt = (3.0875508893757986,0.21357324340764258)
        end
        ## This is the stata J-test
        ov = 1.97522
        @test Je ≈ Jt atol=1e-6
        @test pe ≈ pt atol=1e-6
        @test objval(two_step) ≈ ov atol = 01e-03
        @test objval(two_step) < ov
    end

    context("Instrumental variables -- verified with Stata") do
        include("instrumental_variables.jl")
        cf_stata = [-.0050209, .0708816, .0593661]
        cfe = coef(gmm2s)
        V = vcov(gmm2s)
        for j = 1:2
            @fact cfe[j] --> roughly(cf_stata[j], atol=1e-3)
        end

        V_stata = [ .00478911  -.00272551   .00163441;
                    -.00272551   .01026473  -.00201291;
                    .00163441   -.00201291   .00983372]

        for j = 1:length(V)
            @fact V[j] --> roughly(V_stata[j], atol = 0.01)
        end

        Je, pe = MomentBasedEstimators.J_test(gmm2s)
        ov = 1.191
        @fact objval(gmm2s) --> roughly(ov, atol = 0.02)
    end
    context("Instrumental variables large --- verified by asymptotics") do
        include("instrumental_variables_large.jl")
        @test status(gmm_base)    == :Optimal
        @test status(gmm_one)     == :Optimal
        @test status(el_base)     == :Optimal
        @test status(kl_base)     == :Optimal
        @test status(cue_base)    == :Optimal
        @test status(cue_uncon)   == :Optimal
        @test status(kl_ana_grad) == :Optimal
        @test status(kl_ana_full) == :Optimal

        @test coef(el_base) ≈ coef(kl_base)  atol = 0.01
        @test coef(kl_base) ≈ coef(kl_ana_full)

        @test objval(kl_base) ≈ objval(kl_ana_full)

        @test vcov(kl_base)  ≈ vcov(kl_ana_full)
        @test vcov(kl_base)  ≈ vcov(kl_ana_full)
        @test vcov(gmm_base) ≈ vcov(el_base) atol = 0.01

        @test vcov(el_base, weighted = false)      ≈ vcov(gmm_base) atol = 0.001
        @test vcov(kl_base, weighted = false)      ≈ vcov(gmm_base) atol = 0.001
        @test vcov(kl_ana_grad, weighted = false)  ≈ vcov(gmm_base) atol = 0.001
        @test vcov(kl_ana_full, weighted = false)  ≈ vcov(gmm_base) atol = 0.001
        @test vcov(cue_base, weighted = false)     ≈ vcov(gmm_base) atol = 0.001

        @test vcov(el_base, weighted = true)       ≈ vcov(gmm_base) atol = 0.001
        @test vcov(kl_base, weighted = true)       ≈ vcov(gmm_base) atol = 0.001
        @test vcov(kl_ana_grad, weighted = true)   ≈ vcov(gmm_base) atol = 0.001
        @test vcov(kl_ana_full, weighted = true)   ≈ vcov(gmm_base) atol = 0.001
        @test vcov(cue_base, weighted = true) ≈ vcov(gmm_base) atol =  0.001

        @test stderr(gmm_base)' --> sqrt(vcov(gmm_base))

        @test stderr(kl_base)' --> sqrt(vcov(kl_base))
        @test stderr(el_base)' --> sqrt(vcov(el_base))

        @test stderr(kl_base, weighted = true) --> stderr(kl_base)
        @test stderr(el_base, weighted = true) --> stderr(el_base)

        @test vcov(kl_base_truncated)  --> roughly(vcov(kl_base), 0.01)
        @test vcov(cue_base_truncated) --> roughly(vcov(cue_base), 0.01)

        @test J_test(kl_base_truncated)[1] --> roughly(J_test(kl_base)[1], 0.05)
        @test J_test(cue_base_truncated)[1]  --> roughly(J_test(cue_base)[1], 0.16)
        @test J_test(kl_base_truncated)[2] --> roughly(J_test(kl_base)[2], 0.02)
        @test J_test(cue_base_truncated)[2]  --> roughly(J_test(cue_base)[2], 0.02)

        tmp = J_test(gmm_base)
        @fact tmp[1] --> roughly(7.937738532483664, 1e-07)
        @fact tmp[2] --> roughly(0.5404326890480416, 1e-07)

        tmp = J_test(el_base)
        @fact tmp[1] --> roughly(7.825790062416562, 1e-07)
        @fact tmp[2] --> roughly(0.5517934008321658, 1e-07)

        @fact coeftable(gmm_base).cols[1] --> coef(gmm_base)
        @fact coeftable(gmm_base).cols[2] --> stderr(gmm_base)

        @fact coeftable(el_base).cols[1] --> coef(el_base)
        @fact coeftable(el_base).cols[2] --> stderr(el_base)

    end
    ## context("Instrumental variables large --- verified by asymptotics") do
    ##     Vgmm = vcov(gmm)
    ##     Vmd  = vcov(md)

    ##     for j = 1:length(V)
    ##         @fact Vgmm[j] --> roughly(Vmd[j], atol = 0.00001)
    ##     end

    ##     Vmd_1 = vcov(md, false, :unweighted)
    ##     Vmd_2 = vcov(md, false, :weighted)

    ##     for j = 1:length(V)
    ##         @fact Vmd_1[j] --> roughly(Vmd_2[j], atol = 0.00001)
    ##     end
    ## end
end

## facts("Test utilities") do
##     context("test row_kron") do
##         h = ["a" "b"; "c" "d"]
##         z = ["1" "2" "3"; "4" "5" "6"]
##         want = ["a1" "a2" "a3" "b1" "b2" "b3"; "c4" "c5" "c6" "d4" "d5" "d6"]
##         @fact MomentBasedEstimators.row_kron(h, z) --> want

##         # now test on some bigger matrices
##         a = randn(400, 3)
##         b = randn(400, 5)
##         out = MomentBasedEstimators.row_kron(a, b)
##         @fact size(out) --> (400, 15)

##         rows_good = true
##         for row=1:400
##             rows_good &= out[row, :] == kron(a[row, :], b[row, :])
##         end
##         @fact rows_good --> true
##     end

##     context("test max_args") do
##         foo(x, y, z) = nothing  # standard 3 args
##         bar(x, z=100) = nothing  # standard 2 args with default value
##         baz = (x, y, z)-> nothing  # anonymous 3 args
##         qux(a; b=100) = nothing  # standard 1 with 1 kwarg (kwarg not counted)

##         @fact MomentBasedEstimators.max_args(foo) --> 3
##         @fact MomentBasedEstimators.max_args(bar) --> 2
##         @fact MomentBasedEstimators.max_args(baz) --> 3
##         @fact MomentBasedEstimators.max_args(qux) --> 1
##     end
## end

FactCheck.exitstatus()


text:
module ComplexPhasePortrait
using Reactive, IntervalSets, RecipesBase
import AbstractPlotting
import Images
import Colors
import Colors: RGB, HSL
import AbstractPlotting: plot!, Scene, AbstractScene, ScenePlot
export portrait,
       PTproper, PTcgrid, PTstepphase, PTstepmod, phaseplot

export (..)

abstract type PortraitType <: Any end
struct PTproper <: PortraitType end
struct PTcgrid <: PortraitType end
struct PTstepmod <: PortraitType end
struct PTstepphase <: PortraitType end

const brighten = 0.1

"""
portrait builds a phase portrait over a give complex grid.

The PortraitType determines what type of portrait is generated:

* PTproper - a proper phase portrait
* PTcgrid - a conformal grid
* PTstepphase - a stepped phase plot
* PTstepmod - a stepped modulus plot

Named arguments:

* ctype - is the colormap type. The default "standard" is a basic HSL map
  with H∈[0,360], S = 1.0, and L = 0.5. Setting `ctype="nist"` gives the
  [NIST color scheme](http://dlmf.nist.gov/help/vrml/aboutcolor).
* pres - is the phase step resolution (number of phase jumps), default
  `pres=20`.
"""

function portrait(fval::Array{Complex{Float64},2}; kwargs...)
    portrait(fval, PTproper; kwargs...)
end

function portrait(fval::Array{Complex{Float64},2}, ::Type{PTproper};
                  kwargs...)
    args = baseArgs(fval; kwargs...)
    phaseToImage!(args[1:3]...)
    return args[1]
end

function portrait(fval::Array{Complex{Float64},2}, ::Type{PTcgrid};
                  pres=20, kwargs...)
    (img, nphase, cm, farg) = baseArgs(fval; kwargs...)
    lowb = sqrt(0.75^2*(1.0 - brighten) + brighten)
    black = (sawfun(farg, 1/pres, lowb, 1.0)
             .*sawfun.(log.(abs.(fval)), 2pi/pres, lowb, 1.0))
    phaseToImage!(img, nphase, black, cm)
    return img
end

function portrait(fval::Array{Complex{Float64},2}, ::Type{PTstepmod};
                  pres=20, kwargs...)
    (img, nphase, cm) = baseArgs(fval; kwargs...)
    phaseToImage!(img, nphase, sawfun(log.(abs.(fval)), 2pi/pres, 0.75, 1.0), cm)
    return img
end

function portrait(fval::Array{Complex{Float64},2}, ::Type{PTstepphase};
                  pres=20, kwargs...)
    (img, nphase, cm, farg) = baseArgs(fval; kwargs...)
    phaseToImage!(img, nphase, sawfun(farg, 1/pres, 0.75, 1.0), cm)
    return img
end

function baseArgs(fval::Array{Complex{Float64},2}; kwargs...)
    (farg, nphase, cm) = setupPhase(fval; kwargs...)
    (m, n) = size(nphase)
    img = Array{RGB{Float64}}(undef, m, n)
    return img, nphase, cm, farg
end

function setupPhase(fval; kwargs...)
    cm = baseColorMap(;kwargs...)
    nc = length(cm)
    farg = (angle.(-fval) .+ pi)./2pi
    nphase = stepfun(farg, nc)

    return (farg, nphase, cm)
end

"Color map for complex phase portrait with 600 elements."
function baseColorMap(;ctype="standard")
    if ctype == "nist"
        nc = 900
        cm = range(HSL(0.0, 1.0, 0.5), stop=HSL(360.0, 1.0, 0.5), length=nc)
        idx = [1:Int(nc/6); Int(nc/6)+1:2:Int(nc/2); Int(nc/2)+1:2*Int(nc/3);
              2*Int(nc/3)+1:2:nc]
        cm = cm[idx]
    else
        cm = range(HSL(0.0, 1.0, 0.5), stop=HSL(360.0, 1.0, 0.5), length=600)
    end
    cm = convert(Array{RGB{Float64}}, cm)
end

"Integer step function with period 1 such that [0,1]⟶[1,nmax]."
function stepfun(x, nmax)
    y = x .- floor.(x)
    y = convert(Array{UInt16}, floor.(nmax*y) .+ 1)
end

"Sawtooth function over reals with period dx onto [a,b]."
function sawfun(x, dx, a, b)
    x = x/dx
    x = x .- floor.(x)
    x = a .+ (b - a)*x
end

function phaseToImage!(img, pidx, cm)
    n = size(pidx, 1)
    for j in 1:size(pidx, 2), i in 1:n
        img[n-i+1,j] = cm[pidx[i,j]]
    end
    nothing
end

function phaseToImage!(img, pidx, black, cm)
    n = size(pidx, 1)
    for j in 1:size(pidx, 2), i in 1:n
        img[n-i+1,j] = cm[pidx[i,j]]*black[i,j]
    end
    nothing
end

function _range(d::ClosedInterval)
    a,b = endpoints(d)
    range(a; stop=b, length=500)
end

_range(d::AbstractVector) = d

portrait(x::AbstractVector, y::AbstractVector, f::Function) =
    portrait(convert(AbstractMatrix{ComplexF64}, f.(x' .+ im .*y)))

portrait(x::ClosedInterval, y::ClosedInterval, f::Function) =
    portrait(_range(x), _range(y), f)

## Recipe for Plots

@userplot PhasePlot

@recipe function f(c::PhasePlot)
    xin, yin, ff = c.args
    xx, yy = _range(xin),  _range(yin)
    Z = ff.(xx' .+ im.*yy)
    yflip := false
    @series xx, yy, portrait(Matrix{ComplexF64}(Z[end:-1:1,:]))
end

## Recipe for Makie

@AbstractPlotting.recipe(Phase) do scene
    default_theme(scene)
end

function plot!(plot::Phase{Tuple{X,Y,F}}) where {X<:AbstractVector,Y<:AbstractVector,F<:AbstractMatrix}
    x, y, f = value.(plot[1:3])::Tuple{X,Y,F}
    image!(plot, x, y, portrait(convert(AbstractMatrix{ComplexF64}, f)))
end

# avoid compile time issues
function phase(x::AbstractVector, y::AbstractVector, f::Function; kwds...)
    z = x' .+ im.*y
    a_x,b_x = first(x),last(x)
    a_y,b_y = first(y),last(y)
    r = (b_y-a_y)/(b_x-a_x)
    N = 400
    s = Scene(resolution=(max(N,floor(Int,N/r)), max(N,floor(Int, r *N))))
    phase!(s, x, y, f.(z); limits = FRect(a_x,a_y,b_x-a_x,b_y-a_y), kwds...)
end
function phase!(plot::Union{AbstractScene,ScenePlot}, x::AbstractVector, y::AbstractVector, f::Function; kwds...)
    z = x' .+ im.*y
    phase!(plot, x, y, f.(z); kwds...)
end

phase(x::ClosedInterval, y::ClosedInterval, f::Function; kwds...) =
    phase(_range(x), _range(y), f; kwds...)

phase!(plot::Union{AbstractScene,AbstractPlotting.ScenePlot}, x::ClosedInterval, y::ClosedInterval, f::Function; kwds...) =
    phase!(plot, _range(x), _range(y), r; kwds...)

end # module


text:
using DynamicPolynomials
using MultivariatePolynomials
using Test
using LinearAlgebra

include("mono.jl")
include("poly.jl")
include("comp.jl")

# TODO move to MultivariatePolynomials.jl
@testset "Subs with no variables" begin
    @polyvar x
    t = convert(termtype(x, Int), 3)
    @test t == @inferred subs(t, x => x + 1)
    @test t == @inferred subs(t, x => x + 1.0)
    @test t == @inferred subs(t, x => 1x)
    @test t == @inferred subs(t, x => 1.0x)
    @test t == @inferred subs(t, x => 1.0)
end

module newmodule
    using Test
    import DynamicPolynomials
    @testset "Polyvar macro hygiene" begin
        # Verify that the @polyvar macro works when the package has been activated
        # with `import` instead of `using`.
        DynamicPolynomials.@polyvar x y
        @test isa(x, DynamicPolynomials.PolyVar)
        @test isa(y, DynamicPolynomials.PolyVar)
    end
end

include("mvp.jl")


text:
#=

      Linear Layer

=#

export LinearLayer

type LinearLayer{D<:Device, F<:Float} <: Layer
  # Device
  device::Device

  # The parameters
  weight::matrixParameters(D,F)
  bias::matrixParameters(D,F)

  # State
  output::denseRealMatrix(D,F)
  input_gradient::denseRealMatrix(D,F)

  function LinearLayer(device, weight, bias)
    new(device, weight, bias, array(device, F, 0, 0), array(device, F, 0, 0))
  end
end

"Initialize a linear module with the input and output size"
function LinearLayer{D<:Device, F<:Float}(d::D, ::Type{F}, inputSize::Int64, outputSize::Int64)
  LinearLayer{D,F}(d, matrixParameters(d, F, inputSize, outputSize), matrixParameters(d, F, 1, outputSize))
end

function forward!{D<:Device, F<:Float}(linear::LinearLayer{D,F}, input::DenseMatrix)
  @ensuresize linear.output, size(input,1), size(linear.weight.values, 2)

  gemm!('N', 'N', 1., input, linear.weight.values, 0., linear.output)
  broadcast!(+, linear.output, linear.output, linear.bias.values)
end

function compute_inputgradient!{D<:Device, F<:Float}(linear::LinearLayer{D,F}, input::DenseMatrix, gradOutput::DenseMatrix)
  @ensuresize linear.input_gradient, size(gradOutput,1), size(linear.weight.values, 1)

  gemm!('N', 'T', 1., gradOutput, linear.weight.values, 0., linear.input_gradient)
end

function update_gradient!{D<:Device, F<:Float}(linear::LinearLayer{D,F}, input::DenseMatrix, gradOutput::DenseMatrix, scale::F=1.)
  if !isnull(linear.weight.gradient)
    gemm!('T', 'N', scale, input, gradOutput, 0., get(linear.weight.gradient))
  end

  if !isnull(linear.bias.gradient)
    axpy!(scale, sum(gradOutput, 1), get(linear.bias.gradient))
  end
end


text:
#=--- TerminalGraphics / libsixel.jl --------------------------------------=#

"""
The sub-module `LibSixel` provides some basic wrappers to _libsixel_.
It requires _libsixel_ version 1.0.6 or above for full functionality.
"""
module LibSixel

using ColorTypes
using FixedPointNumbers: N0f8

struct SixelError <: Exception
    msg::String
end

import Base: showerror
showerror(io::IO, err::SixelError) = print(io, err.msg)

const PtrToSixelAllocator = Ptr{Void}

"""
The `SixelEncoder` converts files (given as filenames) or raw image data
into a sixel encoded escape sequence and dumps the data to STDOUT.

Equivalent to the `struct` defined by `libsixel::sixel_encoder_t`.
Apparently these two definitions must always be in sync.
"""
struct SixelEncoder
    ref                   ::Cuint                    # reference counter
    allocator             ::PtrToSixelAllocator
    reqcolors             ::Cint
    color_option          ::Cint
    mapfile               ::Cstring
    builtin_palette       ::Cint
    method_for_diffuse    ::Cint
    method_for_largest    ::Cint
    method_for_rep        ::Cint
    quality_mode          ::Cint
    method_for_resampling ::Cint
    loop_mode             ::Cint
    palette_type          ::Cint
    f8bit                 ::Cint
    finvert               ::Cint
    fuse_macro            ::Cint
    fignore_delay         ::Cint
    complexion            ::Cint
    fstatic               ::Cint
    pixelwidth            ::Cint
    pixelheight           ::Cint
    percentwidth          ::Cint
    percentheight         ::Cint
    clipx                 ::Cint
    clipy                 ::Cint
    clipwidth             ::Cint
    clipheight            ::Cint
    clipfirst             ::Cint
    macro_number          ::Cint
    penetrate_multiplexer ::Cint
    encode_policy         ::Cint
    pipe_mode             ::Cint
    verbose               ::Cint
    has_gri_arg_limit     ::Cint
    bgcolor               ::Cstring
    outfd                 ::Cint
    finsecure             ::Cint
    cancel_flag           ::Ptr{Cint}
    dither_cache          ::Ptr{Void}
end

# offset value of pixelformat
const SIXEL_FORMATTYPE_COLOR     = UInt8(0)
const SIXEL_FORMATTYPE_GRAYSCALE = UInt8(1 << 6)
const SIXEL_FORMATTYPE_PALETTE   = UInt8(1 << 7)

@enum( PixelFormat::Cint
     , Pixel_RGB555   = (SIXEL_FORMATTYPE_COLOR     | 0x01) # 15bpp
     , Pixel_RGB565   = (SIXEL_FORMATTYPE_COLOR     | 0x02) # 16bpp
     , Pixel_RGB888   = (SIXEL_FORMATTYPE_COLOR     | 0x03) # 24bpp
     , Pixel_BGR555   = (SIXEL_FORMATTYPE_COLOR     | 0x04) # 15bpp
     , Pixel_BGR565   = (SIXEL_FORMATTYPE_COLOR     | 0x05) # 16bpp
     , Pixel_BGR888   = (SIXEL_FORMATTYPE_COLOR     | 0x06) # 24bpp
     , Pixel_ARGB8888 = (SIXEL_FORMATTYPE_COLOR     | 0x10) # 32bpp
     , Pixel_RGBA8888 = (SIXEL_FORMATTYPE_COLOR     | 0x11) # 32bpp
     , Pixel_GRAY1    = (SIXEL_FORMATTYPE_GRAYSCALE | 0x00) #  1bpp grayscale
     , Pixel_GRAY2    = (SIXEL_FORMATTYPE_GRAYSCALE | 0x01) #  2bpp grayscale
     , Pixel_GRAY4    = (SIXEL_FORMATTYPE_GRAYSCALE | 0x02) #  4bpp grayscale
     , Pixel_GRAY8    = (SIXEL_FORMATTYPE_GRAYSCALE | 0x03) #  8bpp grayscale
     , Pixel_AG88     = (SIXEL_FORMATTYPE_GRAYSCALE | 0x13) # 16bpp gray+alpha
     , Pixel_GA88     = (SIXEL_FORMATTYPE_GRAYSCALE | 0x23) # 16bpp gray+alpha
     , Pixel_PAL1     = (SIXEL_FORMATTYPE_PALETTE   | 0x00) #  1bpp palette
     , Pixel_PAL2     = (SIXEL_FORMATTYPE_PALETTE   | 0x01) #  2bpp palette
     , Pixel_PAL4     = (SIXEL_FORMATTYPE_PALETTE   | 0x02) #  4bpp palette
     , Pixel_PAL8     = (SIXEL_FORMATTYPE_PALETTE   | 0x03) #  8bpp palette
)

# Mapping of Colorant to PixelFormat, only lists the supported ones
const supported_colormaps = Dict(
      ColorTypes.RGB{N0f8}   => Pixel_RGB888,
      ColorTypes.ARGB{N0f8}  => Pixel_ARGB8888,
      ColorTypes.RGBA{N0f8}  => Pixel_RGBA8888,
      ColorTypes.BGR{N0f8}   => Pixel_BGR888,
      ColorTypes.Gray{N0f8}  => Pixel_GRAY8,
      ColorTypes.AGray{N0f8} => Pixel_AG88,
      ColorTypes.GrayA{N0f8} => Pixel_GA88,
)

#= currently not required...
@enum( Colorspace::Cint
     , Colorspace_AUTO     = 0   # choose palette type automatically
     , Colorspace_HLS      = 1   # HLS colorspace (hue, lightness, saturation)
     , Colorspace_RGB      = 2   # RGB colorspace
)

@enum( DitherPaletteType::Cint
     , Palette_MONO_DARK   = 0x0  # monochrome terminal with dark background
     , Palette_MONO_LIGHT  = 0x1  # monochrome terminal with dark background
     , Palette_XTERM16     = 0x2  # xterm 16color
     , Palette_XTERM256    = 0x3  # xterm 256color
     , Palette_VT340_MONO  = 0x4  # vt340 monochrome
     , Palette_VT340_COLOR = 0x5  # vt340 color
     , Palette_G1          = 0x6  # 1bit grayscale
     , Palette_G2          = 0x7  # 2bit grayscale
     , Palette_G4          = 0x8  # 4bit grayscale
     , Palette_G8          = 0x9  # 8bit grayscale
)
=#

const encoder_options = Dict{Symbol,Cint}(
    :INPUT            => 'i',   # -i, --input: specify input file name.
    :OUTPUT           => 'o',   # -o, --output: specify output file name.
    :OUTFILE          => 'o',   # -o, --outfile: specify output file name.
    :SEVENBIT_MODE    => '7',   # -7, --7bit-mode: for 7bit terminals or printers (default)
    :EIGHTBIT_MODE    => '8',   # -8, --8bit-mode: for 8bit terminals or printers
    :COLORS           => 'p',   # -p COLORS, --colors=COLORS: specify number of colors
    :MAPFILE          => 'm',   # -m FILE, --mapfile=FILE: specify set of colors
    :MONOCHROME       => 'e',   # -e, --monochrome: output monochrome sixel image
    :INSECURE         => 'k',   # -k, --insecure: allow to connect to SSL sites without certs
    :INVERT           => 'i',   # -i, --invert: assume the terminal background color
    :HIGH_COLOR       => 'I',   # -I, --high-color: output 15bpp sixel image
    :USE_MACRO        => 'u',   # -u, --use-macro: use DECDMAC and DEVINVM sequences
    :MACRO_NUMBER     => 'n',   # -n MACRONO, --macro-number=MACRONO: specify macro register number
    :COMPLEXION_SCORE => 'C',   # -C COMPLEXIONSCORE, --complexion-score=COMPLEXIONSCORE: specify an number argument for the score of complexion correction.
    :IGNORE_DELAY     => 'g',   # -g, --ignore-delay: render GIF animation without delay
    :STATIC           => 'S',   # -S, --static: render animated GIF as a static image
    :DIFFUSION        => 'd',   # -d DIFFUSIONTYPE, --diffusion=DIFFUSIONTYPE: choose diffusion method which used with -p option.
                                #          DIFFUSIONTYPE is one of them:
                                #            auto     -> choose diffusion type automatically (default)
                                #            none     -> do not diffuse
                                #            fs       -> Floyd-Steinberg method
                                #            atkinson -> Bill Atkinson's method
                                #            jajuni   -> Jarvis, Judice & Ninke
                                #            stucki   -> Stucki's method
                                #            burkes   -> Burkes' method
    :FIND_LARGEST     => 'f',   # -f FINDTYPE, --find-largest=FINDTYPE: choose method for finding the largest dimension of median cut boxes for splitting, make sense only when -p option (color reduction) is specified
                                #         FINDTYPE is one of them:
                                #           auto -> choose finding method automatically (default)
                                #           norm -> simply comparing the range in RGB space
                                #           lum  -> transforming into luminosities before the comparison
    :SELECT_COLOR     => 's',   # -s SELECTTYPE, --select-color=SELECTTYPE: choose the method for selecting representative color from each median-cut box, make sense only when -p option (color reduction) is specified
                                #        SELECTTYPE is one of them:
                                #          auto      -> choose selecting method automatically (default)
                                #          center    -> choose the center of the box
                                #          average    -> calculate the color average into the box
                                #          histogram -> similar with average but considers color histogram
    :CROP             => 'c',   # -c REGION, --crop=REGION: crop source image to fit the specified geometry. REGION should be formatted as '%dx%d+%d+%d'
    :WIDTH            => 'w',   # -w WIDTH, --width=WIDTH: resize image to specified width WIDTH is represented by the following syntax
                                #          auto       -> preserving aspect ratio (default)
                                #          <number>%  -> scale width with given percentage
                                #          <number>   -> scale width with pixel counts
                                #          <number>px -> scale width with pixel counts
    :HEIGHT           => 'h',   # -h HEIGHT, --height=HEIGHT: resize image to specified height HEIGHT is represented by the following syntax
                                #           auto       -> preserving aspect ratio (default)
                                #           <number>%  -> scale height with given percentage
                                #           <number>   -> scale height with pixel counts
                                #           <number>px -> scale height with pixel counts
    :RESAMPLING       => 'r',   # -r RESAMPLINGTYPE, --resampling=RESAMPLINGTYPE: choose resampling filter used with -w or -h option (scaling) RESAMPLINGTYPE is one of them:
                                #          nearest  -> Nearest-Neighbor method
                                #          gaussian -> Gaussian filter
                                #          hanning  -> Hanning filter
                                #          hamming  -> Hamming filter
                                #          bilinear -> Bilinear filter (default)
                                #          welsh    -> Welsh filter
                                #          bicubic  -> Bicubic filter
                                #          lanczos2 -> Lanczos-2 filter
                                #          lanczos3 -> Lanczos-3 filter
                                #          lanczos4 -> Lanczos-4 filter
    :QUALITY          => 'q',   # -q QUALITYMODE, --quality=QUALITYMODE: select quality of color quanlization.
                                #          auto -> decide quality mode automatically (default)
                                #          low  -> low quality and high speed mode
                                #          high -> high quality and low speed mode
                                #          full -> full quality and careful speed mode
    :LOOPMODE         => 'l',   # -l LOOPMODE, --loop-control=LOOPMODE: select loop control mode for GIF animation.
                                #          auto    -> honor the setting of GIF header (default)
                                #          force   -> always enable loop
                                #          disable -> always disable loop
    :PALETTE_TYPE     => 't',   # -t PALETTETYPE, --palette-type=PALETTETYPE: select palette color space type
                                #          auto -> choose palette type automatically (default)
                                #          hls  -> use HLS color space
                                #          rgb  -> use RGB color space
    :BUILTIN_PALETTE  => 'b',   # -b BUILTINPALETTE, --builtin-palette=BUILTINPALETTE: select built-in palette type
                                #          xterm16    -> X default 16 color map
                                #          xterm256   -> X default 256 color map
                                #          vt340mono  -> VT340 monochrome map
                                #          vt340color -> VT340 color map
                                #          gray1      -> 1bit grayscale map
                                #          gray2      -> 2bit grayscale map
                                #          gray4      -> 4bit grayscale map
                                #          gray8      -> 8bit grayscale map
    :ENCODE_POLICY    => 'E',   # -E ENCODEPOLICY, --encode-policy=ENCODEPOLICY: select encoding policy
                                #          auto -> choose encoding policy automatically (default)
                                #          fast -> encode as fast as possible
                                #          size -> encode to as small sixel sequence as possible
    :BGCOLOR          => 'B',   # -B BGCOLOR, --bgcolor=BGCOLOR: specify background color, BGCOLOR is represented by the following syntax
                                #          #rgb
                                #          #rrggbb
                                #          #rrrgggbbb
                                #          #rrrrggggbbbb
                                #          rgb:r/g/b
                                #          rgb:rr/gg/bb
                                #          rgb:rrr/ggg/bbb
                                #          rgb:rrrr/gggg/bbbb
    :PENETRATE        => 'P',   # -P, --penetrate: penetrate GNU Screen using DCS pass-through sequence
    :PIPE_MODE        => 'D',   # -D, --pipe-mode: read source images from stdin continuously
    :VERBOSE          => 'v',   # -v, --verbose: show debugging info
    :VERSION          => 'V',   # -V, --version: show version and license info
    :HELP             => 'H',   # -H, --help: show this help
)

#
#
#  Default encoder instance
#
#

#=-------------------------------------------------------------------------=#
#
#   Binary dependencies & initalisation
#
#=-------------------------------------------------------------------------=#

const depsfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error( "'TerminalGraphics' is not properly installed.\n" *
           "Please run Pkg.build(\"TerminalGraphics\"), then restart Julia." )
end

# use this everywhere accessing libsixel the encoder instance is required
const PtrToSixelEncoder = Ptr{SixelEncoder}

"Handle (pointer) to the default 'libsixel' encoder instance."
const encoder = Ref{PtrToSixelEncoder}(C_NULL)

__init__() = (global encoder[] = encoder_new())

#=-------------------------------------------------------------------------=#
#
#   Function wrappers
#
#=-------------------------------------------------------------------------=#

function encoder_new() ::PtrToSixelEncoder
    enc = Ref{PtrToSixelEncoder}()
    ccall( (:sixel_encoder_new, libsixel), Cint
         , (Ref{PtrToSixelEncoder}, PtrToSixelAllocator)
         , enc, C_NULL
         ) |> throw_on_sixelerror
    return enc[]
end

encoder_unref(enc::PtrToSixelEncoder) ::Void =
    ccall( (:sixel_encoder_unref, libsixel), Void, (PtrToSixelEncoder,), enc )

encoder_setopt(opt::Symbol, val=nothing) = encoder_setopt(encoder[], opt, val)
function encoder_setopt(enc::PtrToSixelEncoder, opt::Symbol, val)
    opt_num = encoder_options[opt]
    v = val === nothing ? C_NULL : string(val)
    ccall( (:sixel_encoder_setopt, libsixel), Cint
         , (PtrToSixelEncoder, Cint, Cstring)
         , enc, opt_num, v) |> throw_on_sixelerror
end

encoder_encode(filename::AbstractString) = encoder_encode(encoder[], filename)
encoder_encode(enc::PtrToSixelEncoder, filename::AbstractString) ::Void =
    ccall( (:sixel_encoder_encode, libsixel), Cint, (PtrToSixelEncoder, Cstring)
         , enc, string(filename)) |> throw_on_sixelerror

encoder_encode_bytes(img::Matrix) = encoder_encode_bytes(encoder[], img)
function encoder_encode_bytes(enc::PtrToSixelEncoder, img::Matrix{C}) where {C<:Colorant}
    eltype(C) !== N0f8 && throw(ArgumentError(
            "'libsixel' can only display 8bit per color channel. " *
            "Convert your image!"))
    pixelformat = supported_colormaps[C]
    width, height = size(img)
    ccall( (:sixel_encoder_encode_bytes, libsixel), Cint
         , (PtrToSixelEncoder, Ptr{Void}, Cint, Cint, Cint, Ptr{Void}, Cint)
         ,  enc, img, width, height, pixelformat, C_NULL, 0
         ) |> throw_on_sixelerror
end

function throw_on_sixelerror(x::Cint) ::Void
    x == 0 && return
    s = ccall( (:sixel_helper_format_error, libsixel), Cstring, (Cint,), x)
    throw(SixelError("LibSixel error: " * Base.unsafe_string(s)))
end

end # module TerminalGraphics.LibSixel

#=--- end of file ---------------------------------------------------------=#


text:
##############
# DiffResult #
##############

@compat type DiffResult{O,V,D<:Tuple}
    value::V
    derivs::D # ith element = ith-order derivative
    function (::Type{DiffResult{O,V,D}}){O,V,D}(value::V, derivs::NTuple{O,Any})
        return new{O,V,D}(value, derivs)
    end
end

"""
    DiffResult(value, derivs::Tuple)

Return a `DiffResult` instance where values will be stored in the provided `value` storage
and derivatives will be stored in the provided `derivs` storage.

Note that the arguments can be `Number`s or `AbstractArray`s, depending on the dimensionality
of your target function.
"""
DiffResult{V,O}(value::V, derivs::NTuple{O,Any}) = DiffResult{O,V,typeof(derivs)}(value, derivs)

"""
    DiffResult(value, derivs...)

Equivalent to `DiffResult(value, derivs::Tuple)`, where `derivs...` is the splatted form of `derivs::Tuple`.
"""
DiffResult(value, derivs...) = DiffResult(value, derivs)

"""
    GradientResult(x::AbstractArray)

Construct a `DiffResult` that can be used for gradient calculations where `x` is the
input to the target function.

Note that `GradientResult` allocates its own storage; `x` is only used for type and
shape information. If you want to allocate storage yourself, use the `DiffResult`
constructor instead.
"""
GradientResult(x::AbstractArray) = DiffResult(first(x), similar(x))

"""
    JacobianResult(x::AbstractArray)

Construct a `DiffResult` that can be used for Jacobian calculations where `x` is the
input to the target function. This method assumes that the target function's output
dimension equals its input dimension.

Note that `JacobianResult` allocates its own storage; `x` is only used for type and
shape information. If you want to allocate storage yourself, use the `DiffResult`
constructor instead.
"""
JacobianResult(x::AbstractArray) = DiffResult(similar(x), similar(x, length(x), length(x)))

"""
    JacobianResult(y::AbstractArray, x::AbstractArray)

Construct a `DiffResult` that can be used for Jacobian calculations where `x` is the
input to the target function, and `y` is the output (e.g. when taking the Jacobian
of `f!(y, x)`).

Like the single argument version, `y` and `x` are only used for type and
shape information and are not stored in the returned `DiffResult`.
"""
JacobianResult(y::AbstractArray, x::AbstractArray) = DiffResult(similar(y), similar(y, length(y), length(x)))

"""
    HessianResult(x::AbstractArray)

Construct a `DiffResult` that can be used for Hessian calculations where `x` is the
input to the target function.

Note that `HessianResult` allocates its own storage; `x` is only used for type and
shape information. If you want to allocate storage yourself, use the `DiffResult`
constructor instead.
"""
HessianResult(x::AbstractArray) = DiffResult(first(x), similar(x), similar(x, length(x), length(x)))

Base.eltype(r::DiffResult) = eltype(typeof(r))
Base.eltype{O,V,D}(::Type{DiffResult{O,V,D}}) = eltype(V)

@compat Base.:(==)(a::DiffResult, b::DiffResult) = a.value == b.value && a.derivs == b.derivs

Base.copy(r::DiffResult) = DiffResult(copy(r.value), map(copy, r.derivs))

# value/value! #
#--------------#

"""
    value(r::DiffResult)

Return the primal value stored in `r`.

Note that this method returns a reference, not a copy. Thus, if `value(r)` is mutable,
mutating `value(r)` will mutate `r`.
"""
value(r::DiffResult) = r.value

"""
    value!(r::DiffResult, x)

Copy `x` into `r`'s value storage, such that `value(r) == x`.
"""
value!(r::DiffResult, x::Number) = (r.value = x; return r)
value!(r::DiffResult, x::AbstractArray) = (copy!(value(r), x); return r)

"""
    value!(f, r::DiffResult, x)

Like `value!(r::DiffResult, x)`, but with `f` applied to each element, such that `value(r) == map(f, x)`.
"""
value!(f, r::DiffResult, x::Number) = (r.value = f(x); return r)
value!(f, r::DiffResult, x::AbstractArray) = (map!(f, value(r), x); return r)

# derivative/derivative! #
#------------------------#

"""
    derivative(r::DiffResult, ::Type{Val{i}} = Val{1})

Return the `ith` derivative stored in `r`, defaulting to the first derivative.

Note that this method returns a reference, not a copy. Thus, if `derivative(r)` is mutable,
mutating `derivative(r)` will mutate `r`.
"""
derivative{i}(r::DiffResult, ::Type{Val{i}} = Val{1}) = r.derivs[i]

"""
    derivative!(r::DiffResult, x, ::Type{Val{i}} = Val{1})

Copy `x` into `r`'s `ith` derivative storage, such that `derivative(r, Val{i}) == x`.
"""
@generated function derivative!{O,i}(r::DiffResult{O}, x::Number, ::Type{Val{i}} = Val{1})
    newderivs = Expr(:tuple, [i == n ? :(x) : :(derivative(r, Val{$n})) for n in 1:O]...)
    return quote
        r.derivs = $newderivs
        return r
    end
end

function derivative!{i}(r::DiffResult, x::AbstractArray, ::Type{Val{i}} = Val{1})
    copy!(derivative(r, Val{i}), x)
    return r
end

"""
    derivative!(f, r::DiffResult, x, ::Type{Val{i}} = Val{1})

Like `derivative!(r::DiffResult, x, Val{i})`, but with `f` applied to each element,
such that `derivative(r, Val{i}) == map(f, x)`.
"""
@generated function derivative!{O,i}(f, r::DiffResult{O}, x::Number, ::Type{Val{i}} = Val{1})
    newderivs = Expr(:tuple, [i == n ? :(f(x)) : :(derivative(r, Val{$n})) for n in 1:O]...)
    return quote
        r.derivs = $newderivs
        return r
    end
end

function derivative!{i}(f, r::DiffResult, x::AbstractArray, ::Type{Val{i}} = Val{1})
    map!(f, derivative(r, Val{i}), x)
    return r
end

# special-cased methods #
#-----------------------#

"""
    gradient(r::DiffResult)

Return the gradient stored in `r` (equivalent to `derivative(r)`).

Note that this method returns a reference, not a copy. Thus, if `gradient(r)` is mutable,
mutating `gradient(r)` will mutate `r`.
"""
gradient(r::DiffResult) = derivative(r)

"""
    gradient!(r::DiffResult, x)

Copy `x` into `r`'s gradient storage, such that `gradient(r) == x`.
"""
gradient!(r::DiffResult, x) = derivative!(r, x)

"""
    gradient!(f, r::DiffResult, x)

Like `gradient!(r::DiffResult, x)`, but with `f` applied to each element,
such that `gradient(r) == map(f, x)`.
"""
gradient!(f, r::DiffResult, x) = derivative!(f, r, x)

"""
    jacobian(r::DiffResult)

Return the Jacobian stored in `r` (equivalent to `derivative(r)`).

Note that this method returns a reference, not a copy. Thus, if `jacobian(r)` is mutable,
mutating `jacobian(r)` will mutate `r`.
"""
jacobian(r::DiffResult) = derivative(r)

"""
    jacobian!(r::DiffResult, x)

Copy `x` into `r`'s Jacobian storage, such that `jacobian(r) == x`.
"""
jacobian!(r::DiffResult, x) = derivative!(r, x)

"""
    jacobian!(f, r::DiffResult, x)

Like `jacobian!(r::DiffResult, x)`, but with `f` applied to each element,
such that `jacobian(r) == map(f, x)`.
"""
jacobian!(f, r::DiffResult, x) = derivative!(f, r, x)

"""
    hessian(r::DiffResult)

Return the Hessian stored in `r` (equivalent to `derivative(r, Val{2})`).

Note that this method returns a reference, not a copy. Thus, if `hessian(r)` is mutable,
mutating `hessian(r)` will mutate `r`.
"""
hessian(r::DiffResult) = derivative(r, Val{2})

"""
    hessian!(r::DiffResult, x)

Copy `x` into `r`'s Hessian storage, such that `hessian(r) == x`.
"""
hessian!(r::DiffResult, x) = derivative!(r, x, Val{2})

"""
    hessian!(f, r::DiffResult, x)

Like `hessian!(r::DiffResult, x)`, but with `f` applied to each element,
such that `hessian(r) == map(f, x)`.
"""
hessian!(f, r::DiffResult, x) = derivative!(f, r, x, Val{2})


text:
module SMC

using Compat
# using ExpFamily
using Distributions

import Base.length, Base.mean

const Int   = Int64
const Float = Float64

include("hmm.jl")
include("kalman.jl")

include("particles.jl")
include("resample.jl")
include("proposal.jl")
include("particlefilter.jl")
include("particlesmoother.jl")

end # module


text:
"""
    haskey(prop::Properties, key::String)

Check if the property-dictionary has the given key. See `Base.haskey`.
"""
haskey(prop::Properties, key::String) = haskey(prop.properties, key)

"""
    delete!(prop::Properties, key::String)

Delete the pair identified by the given key.
See `Base.delete!`.
"""
function delete!(prop::Properties, key::String)
    delete!(prop.properties, key)
    prop
end

"""
    keys(prop::Properties)

Return all keys stored in the property-dictionary.
See `Base.keys`.
"""
keys(prop::Properties) = keys(prop.properties)

"""
    values(prop::Properties)

Return all values stored in the property-dictionary.
See `Base.values`.
"""
values(prop::Properties) = values(prop.properties)

"""
    pairs(prop::Properties)

Return an iterator over `key => value` pairs.
See `Base.pairs`.
"""
pairs(prop::Properties) = pairs(prop.properties)

"""
    function sizehint!(prop::Properties, n::Integer)

From the `Base` docs: "*Suggest that collection `prop` reserve capacity for at least `n` elements.
This can improve performance.*"
"""
function sizehint!(prop::Properties, n::Integer)
    sizehint!(prop.properties, n)
    prop
end

"""
    keytype(prop::Properties)

Return type of `prop`'s keys. From it's definition it is always `String`.
"""
keytype(prop::Properties) = String

"""
    valtype(prop::Properties)

Return type of `prop`'s values. From it's definition it is always `String`.
"""
valtype(prop::Properties) = String


text:
export add_exist_abstract, add_univ_abstract, add_or_abstract

function add_exist_abstract(mgr::Ptr{Manager}, f::Ptr{Node}, cube::Ptr{Node})
    result = Cudd_addExistAbstract(mgr, f, cube)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end

function add_univ_abstract(mgr::Ptr{Manager}, f::Ptr{Node}, cube::Ptr{Node})
    result = Cudd_addUnivAbstract(mgr, f, cube)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end

function add_or_abstract(mgr::Ptr{Manager}, f::Ptr{Node}, cube::Ptr{Node})
    result = Cudd_addOrAbstract(mgr, f, cube)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end


text:
module PCG

include("types.jl")
import .Types

include("arrays_caches.jl")
import .ArraysCaches

include("geometry.jl")
import .Geometry

include("topologies/Topologies.jl")
import .Topologies

include("storages/Storages.jl")
import .Storages

include("universes.jl")
import .Universes

include("neighborhoods.jl")
import .Neighborhoods

include("operations.jl")
import .Operations

include("recorders/Recorders.jl")
import .Recorders

include("spaces.jl")
import .Spaces

end

max_stars_repo_name: Tiendil/pcg-julia

max_issues_repo_name: Tiendil/pcg-julia

max_forks_repo_name: Tiendil/pcg-julia


text:
module GLPKInterfaceLP

using Compat
using Compat.SparseArrays
using Compat.LinearAlgebra

import GLPK
import MathProgBase
const MPB = MathProgBase
using ..GLPKInterfaceBase

export GLPKSolverLP

mutable struct GLPKMathProgModelLP <: GLPKMathProgModel
    inner::GLPK.Prob
    method::Symbol
    param::Union{GLPK.SimplexParam, GLPK.InteriorParam}
    infeasible_bounds::Bool
end

mutable struct GLPKSolverLP <: MPB.AbstractMathProgSolver
    presolve::Bool
    method::Symbol
    opts
    function GLPKSolverLP(;presolve::Bool=false, method::Symbol=:Simplex, opts...)
        method in [:Simplex, :Exact, :InteriorPoint] ||
            error("""
                  Unknown method for GLPK LP solver: $method
                         Allowed methods:
                           :Simplex
                           :Exact
                           :InteriorPoint""")
        new(presolve, method, opts)
    end
end

function Base.copy(m::GLPKMathProgModelLP)

    m2inner = GLPK.Prob()

    GLPK.copy_prob(m2inner, m.inner, GLPK.ON)

    return GLPKMathProgModelLP(m2inner, m.method, deepcopy(m.param), m.infeasible_bounds)
end

function MPB.LinearQuadraticModel(s::GLPKSolverLP)
    if s.method == :Simplex || s.method == :Exact
        param = GLPK.SimplexParam()
        if s.presolve
            param.presolve = GLPK.ON
        end
    elseif s.method == :InteriorPoint
        param = GLPK.InteriorParam()
        if s.presolve
            Compat.@warn("Ignored option: presolve")
        end
    else
        error("This is a bug")
    end
    param.msg_lev = GLPK.MSG_ERR
    for (k,v) in s.opts
        i = findfirst(x->x==k, fieldnames(typeof(param)))
        if (VERSION < v"0.7-" && i > 0) || (VERSION >= v"0.7-" && i !== nothing)
            t = typeof(param).types[i]
            setfield!(param, i, convert(t, v))
        else
            Compat.@warn("Ignored option: $(string(k))")
        end
    end
    lpm = GLPKMathProgModelLP(GLPK.Prob(), s.method, param, false)
    return lpm
end

function MPB.setparameters!(s::GLPKSolverLP; mpboptions...)
    opts = collect(Any, s.opts)
    for (optname, optval) in mpboptions
        if optname == :TimeLimit
            push!(opts, (:tm_lim,round(Int,1000*optval))) # milliseconds
        elseif optname == :Silent
            if optval == true
                push!(opts, (:msg_lev,GLPK.MSG_OFF))
            end
        else
            error("Unrecognized parameter $optname")
        end
    end
    s.opts = opts
    nothing
end

function MPB.setparameters!(m::GLPKMathProgModelLP; mpboptions...)
    for (optname, optval) in mpboptions
        if optname == :TimeLimit
            m.param.tm_lim = round(Int,1000*optval)
        elseif optname == :Silent
            if optval == true
                m.param.msg_lev = GLPK.MSG_OFF
                m.smplxparam.msg_lev = GLPK.MSG_OFF
            end
        else
            error("Unrecognized parameter $optname")
        end
    end
end

function MPB.optimize!(lpm::GLPKMathProgModelLP)
    lpm.infeasible_bounds = false
    lp = lpm.inner
    for c in 1:MPB.numvar(lpm)
        if GLPK.get_col_lb(lp, c) > GLPK.get_col_ub(lp, c)
            lpm.infeasible_bounds = true
            break
        end
    end
    if !lpm.infeasible_bounds
        for r in 1:MPB.numconstr(lpm)
            if GLPK.get_row_lb(lp, r) > GLPK.get_row_ub(lp, r)
                lpm.infeasible_bounds = true
                break
            end
        end
    end
    if !lpm.infeasible_bounds
        if lpm.method == :Simplex
            solve = GLPK.simplex
        elseif lpm.method == :Exact
            solve = GLPK.exact
        elseif lpm.method == :InteriorPoint
            solve = GLPK.interior
        else
            error("bug")
        end
        return solve(lpm.inner, lpm.param)
    end
end

function MPB.status(lpm::GLPKMathProgModelLP)
    if lpm.infeasible_bounds
        return :Infeasible
    end
    if lpm.method == :Simplex || lpm.method == :Exact
        get_status = GLPK.get_status
    elseif lpm.method == :InteriorPoint
        get_status = GLPK.ipt_status
    else
        error("bug")
    end
    s = get_status(lpm.inner)
    if s == GLPK.OPT
        return :Optimal
    elseif s == GLPK.INFEAS
        return :Infeasible
    elseif s == GLPK.UNBND
        return :Unbounded
    elseif s == GLPK.FEAS
        return :Feasible
    elseif s == GLPK.NOFEAS
        return :Infeasible
    elseif s == GLPK.UNDEF
        return :Undefined
    else
        error("Internal library error")
    end
end

function MPB.getobjval(lpm::GLPKMathProgModelLP)
    if lpm.infeasible_bounds
        if GLPK.get_obj_dir(lpm.inner) == GLPK.MAX
            return -Inf
        else
            return Inf
        end
    end
    if lpm.method == :Simplex || lpm.method == :Exact
        get_obj_val = GLPK.get_obj_val
    elseif lpm.method == :InteriorPoint
        get_obj_val = GLPK.ipt_obj_val
    else
        error("bug")
    end
    return get_obj_val(lpm.inner)
end

function check_feasible_bounds(lpm::GLPKMathProgModelLP, name::String)
    if lpm.infeasible_bounds
        error("$name is not available when some constraint bounds are infeasible (lower bound > upper bound)")
    end
end

function MPB.getsolution(lpm::GLPKMathProgModelLP)
    check_feasible_bounds(lpm, "getsolution")
    lp = lpm.inner
    n = GLPK.get_num_cols(lp)

    if lpm.method == :Simplex || lpm.method == :Exact
        get_col_prim = GLPK.get_col_prim
    elseif lpm.method == :InteriorPoint
        get_col_prim = GLPK.ipt_col_prim
    else
        error("bug")
    end

    return [get_col_prim(lp, i) for i in 1:n]
end

function MPB.getconstrsolution(lpm::GLPKMathProgModelLP)
    check_feasible_bounds(lpm, "getconstrsolution")
    lp = lpm.inner
    m = GLPK.get_num_rows(lp)

    if lpm.method == :Simplex || lpm.method == :Exact
        get_row_prim = GLPK.get_row_prim
    elseif lpm.method == :InteriorPoint
        get_row_prim = GLPK.ipt_row_prim
    else
        error("bug")
    end

    return [get_row_prim(lp, i) for i in 1:m]
end

function MPB.getreducedcosts(lpm::GLPKMathProgModelLP)
    check_feasible_bounds(lpm, "getreducedcosts")
    lp = lpm.inner
    n = GLPK.get_num_cols(lp)

    if lpm.method == :Simplex || lpm.method == :Exact
        get_col_dual = GLPK.get_col_dual
    elseif lpm.method == :InteriorPoint
        get_col_dual = GLPK.ipt_col_dual
    else
        error("bug")
    end

    return [get_col_dual(lp, i) for i in 1:n]
end

function MPB.getconstrduals(lpm::GLPKMathProgModelLP)
    check_feasible_bounds(lpm, "getconstrduals")

    lp = lpm.inner
    m = GLPK.get_num_rows(lp)

    if lpm.method == :Simplex || lpm.method == :Exact
        get_row_dual = GLPK.get_row_dual
    elseif lpm.method == :InteriorPoint
        get_row_dual = GLPK.ipt_row_dual
    else
        error("bug")
    end

    return [get_row_dual(lp, i) for i in 1:m]
end

# The functions getinfeasibilityray and getunboundedray are adapted from code
# taken from the LEMON C++ optimization library. This is the copyright notice:
#
### Copyright (C) 2003-2010
### Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
### (Egervary Research Group on Combinatorial Optimization, EGRES).
###
### Permission to use, modify and distribute this software is granted
### provided that this copyright notice appears in all copies. For
### precise terms see the accompanying LICENSE file.
###
### This software is provided "AS IS" with no warranty of any kind,
### express or implied, and with no claim as to its suitability for any
### purpose.

function MPB.getinfeasibilityray(lpm::GLPKMathProgModelLP)
    if lpm.infeasible_bounds
        # See https://github.com/JuliaOpt/GLPKMathProgInterface.jl/pull/34
        return zeros(MPB.numconstr(lpm))
    end

    lp = lpm.inner

    if lpm.method == :Simplex || lpm.method == :Exact
    elseif lpm.method == :InteriorPoint
        error("getinfeasibilityray is not available when using the InteriorPoint method")
    else
        error("bug")
    end

    m = GLPK.get_num_rows(lp)

    ray = zeros(m)

    ur = GLPK.get_unbnd_ray(lp)
    if ur != 0
        if ur <= m
            k = ur
            get_stat = GLPK.get_row_stat
            get_bind = GLPK.get_row_bind
            get_prim = GLPK.get_row_prim
            get_ub = GLPK.get_row_ub
        else
            k = ur - m
            get_stat = GLPK.get_col_stat
            get_bind = GLPK.get_col_bind
            get_prim = GLPK.get_col_prim
            get_ub = GLPK.get_col_ub
        end

        get_stat(lp, k) == GLPK.BS || error("unbounded ray is primal (use getunboundedray)")

        ray[get_bind(lp, k)] = (get_prim(lp, k) > get_ub(lp, k)) ? -1 : 1

        GLPK.btran(lp, ray)
    else
        eps = 1e-7
        for i = 1:m
            idx = GLPK.get_bhead(lp, i)
            if idx <= m
                k = idx
                get_prim = GLPK.get_row_prim
                get_ub = GLPK.get_row_ub
                get_lb = GLPK.get_row_lb
            else
                k = idx - m
                get_prim = GLPK.get_col_prim
                get_ub = GLPK.get_col_ub
                get_lb = GLPK.get_col_lb
            end

            res = get_prim(lp, k)
            if res > get_ub(lp, k) + eps
                ray[i] = -1
            elseif res < get_lb(lp, k) - eps
                ray[i] = 1
            else
                continue # ray[i] == 0
            end

            if idx <= m
                ray[i] *= GLPK.get_rii(lp, k)
            else
                ray[i] /= GLPK.get_sjj(lp, k)
            end
        end

        GLPK.btran(lp, ray)

        for i = 1:m
            ray[i] /= GLPK.get_rii(lp, i)
        end
    end

    return ray
end

function MPB.getunboundedray(lpm::GLPKMathProgModelLP)
    check_feasible_bounds(lpm, "getreducedcosts")

    lp = lpm.inner

    if lpm.method == :Simplex || lpm.method == :Exact
    elseif lpm.method == :InteriorPoint
        error("getunboundedray is not available when using the InteriorPoint method")
    else
        error("bug")
    end

    m = GLPK.get_num_rows(lp)
    n = GLPK.get_num_cols(lp)

    ray = zeros(n)

    ur = GLPK.get_unbnd_ray(lp)
    if ur != 0
        if ur <= m
            k = ur
            get_stat = GLPK.get_row_stat
            get_dual = GLPK.get_row_dual
        else
            k = ur - m
            get_stat = GLPK.get_col_stat
            get_dual = GLPK.get_col_dual
            ray[k] = 1
        end

        get_stat(lp, k) != GLPK.BS || error("unbounded ray is dual (use getinfeasibilityray)")

        for (ri, rv) in zip(GLPK.eval_tab_col(lp, ur)...)
            ri > m && (ray[ri - m] = rv)
        end

        if (GLPK.get_obj_dir(lp) == GLPK.MAX) ⊻ (get_dual(lp, k) > 0)
            ray .*= -1.0
        end
    else
        for i = 1:n
            ray[i] = GLPK.get_col_prim(lp, i)
        end
    end

    return ray
end

end


text:
immutable AffineTransform{T, Tv, N}
    scalefwd::Matrix{T}
    scaleinv::Matrix{T}
    offset::Vector{Tv}
    temp::Vector{Tv}

    function (::Type{AffineTransform{T,Tv,N}}){T,Tv,N}(S::AbstractMatrix, iS::AbstractMatrix, t::AbstractVector)
        if size(iS) != (N,N) || size(S) != (N,N) || size(t) != (N,)
            throw(DimensionMismatch("Inputs must be $N-dimensional"))
        end
        temp = Vector{eltype(t)}(length(t))
        new{T,Tv,N}(S, iS, t, temp)
    end
end
function AffineTransform(S, offset)
    A, iA = promote(S, inv(S))
    AffineTransform{eltype(A), eltype(offset), length(offset)}(A, iA, offset)
end

ndims{T,Tv,N}(tform::AffineTransform{T,Tv,N}) = N
eltype{T,Tv}(tf::AffineTransform{T,Tv}) = Tv

function show(io::IO, tf::AffineTransform)
    println(io, typeof(tf), ":")
    println(io, "matrix: ", tf.scalefwd)
    println(io, "translation: ", tf.offset)
end

scale(tf::AffineTransform, s::Vector) = AffineTransform(tf.scalefwd * Diagonal(s), tf.offset)
scale(s::Vector, tf::AffineTransform) = AffineTransform(Diagonal(s) * tf.scalefwd, tf.offset.*s)

*(a::AffineTransform, v::AbstractVecOrMat) = tformfwd(a, v)
\(a::AffineTransform, x::AbstractVecOrMat) = tforminv(a, x)

if VERSION < v"0.6.0-dev.1996"
    function tformfwd{T,Tv}(a::AffineTransform{T,Tv}, x::AbstractVector)
        tformfwd!(Array(Base.promote_eltype_op(*, a.scalefwd, x), length(x)), a, x)
    end
    tforminv{T,Tv}(a::AffineTransform{T,Tv}, x::AbstractVector) = tforminv!(Array(Base.promote_eltype_op((*), a.scalefwd, x), length(x)), a, x)
    tformfwd{T,Tv}(a::AffineTransform{T,Tv}, x::AbstractMatrix) = tformfwd!(Array(Base.promote_eltype_op((*), a.scalefwd, x), size(x,1), size(x,2)), a, x)
    tforminv{T,Tv}(a::AffineTransform{T,Tv}, x::AbstractMatrix) = tforminv!(Array(Base.promote_eltype_op((*), a.scalefwd, x), size(x,1), size(x,2)), a, x)
else
    function tformfwd{T,Tv}(a::AffineTransform{T,Tv}, x::AbstractVector)
        tformfwd!(Vector{eltype(a.scalefwd*x)}(length(x)), a, x)
    end
    tforminv{T,Tv}(a::AffineTransform{T,Tv}, x::AbstractVector) = tforminv!(Vector{eltype(a.scalefwd*x)}(length(x)), a, x)
    tformfwd{T,Tv}(a::AffineTransform{T,Tv}, x::AbstractMatrix) = tformfwd!(Matrix{eltype(a.scalefwd*x)}(size(x)), a, x)
    tforminv{T,Tv}(a::AffineTransform{T,Tv}, x::AbstractMatrix) = tforminv!(Matrix{eltype(a.scalefwd*x)}(size(x)), a, x)
end

# Implementation of a.scalefwd*X .+ a.offset
function tformfwd!(dest, a::AffineTransform, X::AbstractVecOrMat)
    A_mul_B!(dest, a.scalefwd, X)
    @inbounds for j = 1:size(X,2), i = 1:size(X,1)
        dest[i,j] += a.offset[i]
    end
    dest
end

# Implementation of a.scaleinv*(X .- a.offset)
function tforminv!(dest, a::AffineTransform, X::AbstractVecOrMat)
    size(X,1) == ndims(a) || throw(DimensionMismatch("Dimensionality $(size(X,1)) of inputs does not match $(ndims(a)) for the transformation"))
    A_mul_B!(dest, a.scaleinv, X)
    A_mul_B!(a.temp, a.scaleinv, a.offset)
    @inbounds for j = 1:size(X,2), i = 1:size(X,1)
        dest[i,j] -= a.temp[i]
    end
    dest
end

# Faster versions for particular dimensionalities (with vectors)
for VecType in (AbstractVector, CartesianIndex)
    @eval begin
        function tforminv!{T,Tv}(v, a::AffineTransform{T,Tv,3}, x::$VecType)
            length(x) == 3 || throw(DimensionMismatch("vector must be 3-dimensional"))
            t = a.temp
            o = a.offset
            @inbounds t[1] = x[1]-o[1]
            @inbounds t[2] = x[2]-o[2]
            @inbounds t[3] = x[3]-o[3]
            gemv3!(v, a.scaleinv, t)
        end
        function tformfwd!{T,Tv}(x, a::AffineTransform{T,Tv,3}, v::$VecType)
            length(x) == 3 || throw(DimensionMismatch("vector must be 3-dimensional"))
            gemv3!(x, a.scalefwd, v)
            o = a.offset
            @inbounds x[1] += o[1]
            @inbounds x[2] += o[2]
            @inbounds x[3] += o[3]
            x
        end
    end
end

function tforminv{T,Tv}(a::AffineTransform{T,Tv,3}, x1::Number, x2::Number, x3::Number)
    o = a.offset
    @inbounds t1 = x1 - o[1]
    @inbounds t2 = x2 - o[2]
    @inbounds t3 = x3 - o[3]
    gemv3(a.scaleinv, t1, t2, t3)
end

function tformfwd{T,Tv}(a::AffineTransform{T,Tv,3}, x1::Number, x2::Number, x3::Number)
    t1, t2, t3 = gemv3(a.scalefwd, x1, x2, x3)
    o = a.offset
    @inbounds r1 = t1+o[1]
    @inbounds r2 = t2+o[2]
    @inbounds r3 = t3+o[3]
    r1, r2, r3
end

for VecType in (AbstractVector, CartesianIndex)
    @eval begin
        function tforminv!{T,Tv}(v, a::AffineTransform{T,Tv,2}, x::$VecType)
            length(x) == 2 || throw(DimensionMismatch("vector must be 2-dimensional"))
            t = a.temp
            o = a.offset
            @inbounds t[1] = x[1]-o[1]
            @inbounds t[2] = x[2]-o[2]
            gemv2!(v, a.scaleinv, t)
        end
        function tformfwd!{T,Tv}(x, a::AffineTransform{T,Tv,2}, v::$VecType)
            length(x) == 2 || throw(DimensionMismatch("vector must be 2-dimensional"))
            gemv2!(x, a.scalefwd, v)
            t = a.offset
            @inbounds x[1] += t[1]
            @inbounds x[2] += t[2]
            x
        end
    end
end

function tforminv{T,Tv}(a::AffineTransform{T,Tv,2}, x1::Number, x2::Number)
    o = a.offset
    @inbounds t1 = x1 - o[1]
    @inbounds t2 = x2 - o[2]
    gemv2(a.scaleinv, t1, t2)
end

function tformfwd{T,Tv}(a::AffineTransform{T,Tv,2}, x1::Number, x2::Number)
    t1, t2 = gemv2(a.scalefwd, x1, x2)
    o = a.offset
    @inbounds r1 = t1+o[1]
    @inbounds r2 = t2+o[2]
    r1, r2
end

# For just the 3rd coordinate (we sometimes only need that one)
function tforminv3rd{T,Tv}(a::AffineTransform{T,Tv,3}, x::AbstractVector)
    o = a.offset
    A = a.scaleinv
    A[3]*(x[1]-o[1]) + A[6]*(x[2]-o[2]) + A[9]*(x[3]-o[3])
end

function scaled_add3!(x, s::Number, r)
    @inbounds x[1]+=s*r[1]
    @inbounds x[2]+=s*r[2]
    @inbounds x[3]+=s*r[3]
end

function gemv3!(res, A, x)
    @inbounds res[1] = A[1]*x[1] + A[4]*x[2] + A[7]*x[3]
    @inbounds res[2] = A[2]*x[1] + A[5]*x[2] + A[8]*x[3]
    @inbounds res[3] = A[3]*x[1] + A[6]*x[2] + A[9]*x[3]
    res
end

function gemv3(A, x1::Number, x2::Number, x3::Number)
    @inbounds r1 = A[1]*x1 + A[4]*x2 + A[7]*x3
    @inbounds r2 = A[2]*x1 + A[5]*x2 + A[8]*x3
    @inbounds r3 = A[3]*x1 + A[6]*x2 + A[9]*x3
    r1, r2, r3
end

function scaled_add2!(x, s::Number, r)
    @inbounds x[1]+=s*r[1]
    @inbounds x[2]+=s*r[2]
end

function gemv2!(res, A, x)
    @inbounds res[1] = A[1]*x[1] + A[3]*x[2]
    @inbounds res[2] = A[2]*x[1] + A[4]*x[2]
    res
end

function gemv2(A, x1::Number, x2::Number)
    @inbounds r1 = A[1]*x1 + A[3]*x2
    @inbounds r2 = A[2]*x1 + A[4]*x2
    r1, r2
end

### Identity transformation
tformeye(T, Tv, nd) = AffineTransform(eye(T, nd), zeros(Tv, nd))
tformeye(T, nd) = tformeye(T, T, nd)
tformeye(nd) = AffineTransform(eye(nd), zeros(nd))

### Rotations
# Right-hand rotation around an axis (axis points towards viewer)

rotation2(angle) = [cos(angle) -sin(angle); sin(angle) cos(angle)]

function tformrotate(angle)
    A = rotation2(angle)
    AffineTransform(A, zeros(eltype(A),2))
end

# The following assumes uaxis is normalized
function _rotation3(uaxis::Vector, angle)
    if length(uaxis) != 3
        error("3d rotations only")
    end
    ux, uy, uz = uaxis[1], uaxis[2], uaxis[3]
    c = cos(angle)
    s = sin(angle)
    cm = one(typeof(c)) - c
    R = [c+ux*ux*cm       ux*uy*cm-uz*s      ux*uz*cm+uy*s;
         uy*ux*cm+uz*s      c+uy*uy*cm       uy*uz*cm-ux*s;
         uz*ux*cm-uy*s    uz*uy*cm+ux*s      c+uz*uz*cm]
end

function rotation3{T}(axis::Vector{T}, angle)
    n = norm(axis)
    axisn = n>0 ? axis/n : (tmp = zeros(T,length(axis)); tmp[1] = 1)
    _rotation3(axisn, angle)
end

# Angle/axis representation where the angle is the norm of the vector (so axis is not normalized)
function rotation3{T}(axis::Vector{T})
    n = norm(axis)
    axisn = n>0 ? axis/n : (tmp = zeros(typeof(one(T)/1),length(axis)); tmp[1] = 1; tmp)
    _rotation3(axisn, n)
end

# Convert rotation matrix to axis representation (angle/axis with the norm of the vector coding the angle)
function rotationparameters{T}(R::Matrix{T})
    size(R, 1) == size(R, 2) || error("Matrix must be square")
    if size(R, 1) == 2
        return [atan2(-R[1,2],R[1,1])]
    elseif size(R, 1) == 3
        # See Itzhack Y. Bar-Itzhack.  "New Method for Extracting the Quaternion from a Rotation Matrix",
        # Journal of Guidance, Control, and Dynamics, Vol. 23, No. 6 (2000), pp. 1085-1087.
        # On https://en.wikipedia.org/wiki/Rotation_matrix#Conversions
        # This is said to be robust even if R is not orthogonal
        Q12p, Q13p, Q23p = R[1,2]+R[2,1], R[1,3]+R[3,1], R[2,3]+R[3,2]
        Q12m, Q13m, Q23m = R[1,2]-R[2,1], R[1,3]-R[3,1], R[2,3]-R[3,2]
        K = [R[1,1]-R[2,2]-R[3,3]  Q12p  Q13p  Q23m;
             Q12p  R[2,2]-R[1,1]-R[3,3]  Q23p -Q13m;
             Q13p  Q23p  R[3,3]-R[1,1]-R[2,2]  Q12m;
             Q23m -Q13m  Q12m  R[1,1]+R[2,2]+R[3,3]]
        D, V = eig(K)
        imax = indmax(D)
        q = V[:,imax]
        ax = q[1:3]
        n = norm(ax)
        if n > 0
            ax = ax/n
        end
        th = -2*atan(n/q[4])
        return th*ax
    else
        error("Rotations in $(size(R, 1)) dimensions not supported")
    end
end

tformtranslate(trans::Vector) = AffineTransform(eye(length(trans)), trans)

tformscale(scale::Real, ndims::Int) = AffineTransform(scale*eye(typeof(scale),ndims), zeros(typeof(scale),ndims))

function tformrotate(axis::Vector, angle)
    if length(axis) == 3
        return AffineTransform(rotation3(axis, angle), zeros(eltype(axis),3))
    else
        error("Dimensionality ", length(axis), " not supported")
    end
end

function tformrotate(x::Vector)
    if length(x) == 3
        return AffineTransform(rotation3(x), zeros(eltype(x),3))
    else
        error("Dimensionality ", length(x), " not supported")
    end
end

# The following is useful in optimization routines, which tend to work on parameter vectors
function tformrigid(p::Vector)
    if length(p) == 3
        # 2d: [angle, shiftx, shifty]
        return AffineTransform(rotation2(p[1]), p[2:3])
    elseif length(p) == 6
        # 3d: [axis, shift]
        return AffineTransform(rotation3(p[1:3]), p[4:6])
    else
        error("Rigid transformations with $(length(p)) parameters not supported")
    end
end

function(*)(T2::AffineTransform, T1::AffineTransform)
    A1 = T1.scalefwd
    A2 = T2.scalefwd
    o = T2.offset + A2*T1.offset
    AffineTransform(A2*A1, o)
end

function(\)(T2::AffineTransform, T1::AffineTransform)
    A1 = T1.scalefwd
    A2 = T2.scaleinv
    o = A2*(T1.offset-T2.offset)
    AffineTransform(A2*A1, o)
end