text:
module momap2para3_2
using FinEtools
using FinEtools.MeshSelectionModule: vselect
using FinEtools.MeshExportModule
using FinEtools.MeshExportModule: VTKWrite
using Test
import LinearAlgebra: norm
function test()
    A = 50.0*phun("m") # length  of loaded rectangle
    B = 200.0*phun("m") # length  of loaded rectangle
    C = 100.0*phun("m") # span of the plate

    # Select how find the mesh should be
    Refinement = 2
    nA, nB, nC = Refinement * 1, Refinement * 2, Refinement * 4;
    xs = reshape(collect(linearspace(0.0, A, nA + 1)), nA + 1, 1)
    ys = reshape(collect(linearspace(0.0, B, nB + 1)), nB + 1, 1)
    zs = reshape(collect(linearspace(0.0, C, nC + 1)), nC + 1, 1)
    fensc,fesc = T10blockx(xs, ys, zs, :b)
    fc = NodalField(zeros(count(fensc), 1))
    for i = 1:count(fensc)
        x, y, z = fensc.xyz[i, :]
        fc.values[i, :] .= sin(2*x/A) * cos(6.5*y/B) * sin(3*z/C-1.0)
    end
    File = "momap2para3_2-coarse.vtu"
    VTKWrite.vtkwrite(File, fensc, fesc; scalars = [("fc", fc.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    Refinement = Refinement + 1
    nA, nB, nC = Refinement * 1, Refinement * 2, Refinement * 4;
    xs = reshape(collect(linearspace(0.0, A, nA + 1)), nA + 1, 1)
    ys = reshape(collect(linearspace(0.0, B, nB + 1)), nB + 1, 1)
    zs = reshape(collect(linearspace(0.0, C, nC + 1)), nC + 1, 1)
    fensf,fesf = T10blockx(xs, ys, zs, :b)
    ff = NodalField(zeros(count(fensf), 1))
    tolerance = min(A/nA, B/nB, C/nC)/1000.0

    referenceff = NodalField(zeros(count(fensf), 1))
    for i = 1:count(fensf)
        x, y, z = fensf.xyz[i, :]
        referenceff.values[i, :] .= sin(2*x/A) * cos(6.5*y/B) * sin(3*z/C-1.0)
    end
    File = "momap2para3_2-reference.vtu"
    VTKWrite.vtkwrite(File, fensf, fesf; scalars = [("referenceff", referenceff.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    ff = transferfield!(ff, fensf, fesf, fc, fensc, fesc, tolerance)
    File = "momap2para3_2-fine.vtu"
    VTKWrite.vtkwrite(File, fensf, fesf; scalars = [("ff", ff.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    diffff = NodalField(referenceff.values - ff.values)
    femm  = FEMMBase(IntegDomain(fesf, SimplexRule(3, 4)))
    geom = NodalField(fensf.xyz)
    error = integratefieldfunction(femm, geom, diffff, (x, v) -> norm(v), 0.0)
    ref = integratefieldfunction(femm, geom, referenceff, (x, v) -> norm(v), 0.0)
    # println("error/ref = $(error/ref)")
    @test abs(error/ref - 0.02541369940759616) < 1.0e-4
end
end
using .momap2para3_2
 momap2para3_2.test()

module mmmQ4blockneous2mm_2
using FinEtools
using FinEtools.MeshExportModule: VTKWrite
using Test
function test()
    rho=1.21*1e-9;# mass density
    c =345.0*1000;# millimeters per second
    bulk= c^2*rho;
    Lx=1900.0;# length of the box, millimeters
    Ly=800.0; # length of the box, millimeters

    fens,fes = Q8block(Lx,Ly,3,2); # Mesh
    # show(fes.conn)
    fens.xyz = xyz3(fens)
    fens.xyz[:, 3] .= (fens.xyz[:, 1].^2 + fens.xyz[:, 2].^2)/1.e3
    File = "mesh.vtu"
    VTKWrite.vtkwrite(File, fens, fes)
    # rm(File)
    # @async run(`"paraview.exe" $File`)
end
end
using .mmmQ4blockneous2mm_2
mmmQ4blockneous2mm_2.test()

module miimportexportm_2
using FinEtools
using FinEtools.MeshImportModule
using FinEtools.MeshExportModule: VTKWrite
using Test
function test()
    output = MeshImportModule.import_NASTRAN(dirname(@__FILE__) * "/" * "Slot-coarser.nas";
    allocationchunk = 13)
    # show(fes.conn[count(fes), :])
    File = "Slot-coarser.vtu"
    VTKWrite.vtkwrite(File, output["fens"], output["fesets"][1])
    rm(File)
    @test output["fesets"][1].conn[count(output["fesets"][1]), :] == NTuple{10,Int64}[(143, 140, 144, 138, 361, 363, 176, 519, 781, 520)]
    # @async run(`"paraview.exe" $File`)
end
end
using .miimportexportm_2
miimportexportm_2.test()

module mmmLshapemmm_2
using FinEtools
using FinEtools.MeshExportModule: VTKWrite
using Test
function test()
    # println("""
    # Meshing  for an L-shaped membrane.
    # """)

    t0 = time()

    W= 100. # strip width
    L = 200. # length of the strip
    nL = 5
    nW = 5
    tolerance = W/nW/1.0e5
    Meshes = Array{Tuple{FENodeSet, AbstractFESet},1}()
    push!(Meshes, Q4quadrilateral([0.0 0.0; W W], nW, nW))
    push!(Meshes, Q4quadrilateral([-L 0.0; 0.0 W], nL, nW))
    push!(Meshes, Q4quadrilateral([0.0 -L; W 0.0], nW, nL))
    fens, outputfes = mergenmeshes(Meshes, tolerance);
    fes = cat(outputfes[1], cat(outputfes[2], outputfes[3]))
    @test count(fens) == 96
    @test fes.conn == NTuple{4,Int64}[(1, 2, 8, 7), (7, 8, 14, 13), (13, 14, 20, 19), (19, 20, 26, 25), (25, 26, 32, 31), (2, 3, 9, 8), (8, 9, 15, 14), (14, 15, 21, 20), (20, 21, 27, 26), (26, 27, 33, 32), (3, 4, 10, 9), (9, 10, 16, 15), (15, 16, 22, 21), (21, 22, 28, 27), (27, 28, 34, 33), (4, 5, 11, 10), (10, 11, 17, 16), (16, 17, 23, 22), (22, 23, 29, 28), (28, 29, 35, 34), (5, 6, 12, 11), (11, 12, 18, 17), (17, 18, 24, 23), (23, 24, 30, 29), (29, 30, 36, 35), (37, 38, 43, 42), (42, 43, 48, 47), (47, 48, 53, 52), (52, 53, 58, 57), (57, 58, 63, 62), (38, 39, 44, 43), (43, 44, 49, 48), (48, 49, 54, 53), (53, 54, 59, 58), (58, 59, 64, 63), (39, 40, 45, 44), (44, 45, 50, 49), (49, 50, 55, 54), (54, 55, 60, 59), (59, 60, 65, 64), (40, 41, 46, 45), (45, 46, 51, 50), (50, 51, 56, 55), (55, 56, 61, 60), (60, 61, 66, 65), (41, 1, 7, 46), (46, 7, 13, 51), (51, 13, 19, 56), (56, 19, 25, 61), (61, 25, 31, 66), (67, 68, 74, 73), (73, 74, 80, 79), (79, 80, 86, 85), (85, 86, 92, 91), (91, 92, 2, 1), (68, 69, 75, 74), (74, 75, 81, 80), (80, 81, 87, 86), (86, 87, 93, 92), (92, 93, 3, 2), (69, 70, 76, 75), (75, 76, 82, 81), (81, 82, 88, 87), (87, 88, 94, 93), (93, 94, 4, 3), (70, 71, 77, 76), (76, 77, 83, 82), (82, 83, 89, 88), (88, 89, 95, 94), (94, 95, 5, 4), (71, 72, 78, 77), (77, 78, 84, 83), (83, 84, 90, 89), (89, 90, 96, 95), (95, 96, 6, 5)]
    geom = NodalField(fens.xyz)

    File =  "L_shape.vtu"
    VTKWrite.vtkwrite(File, connasarray(fes), geom.values, VTKWrite.Q4);
    # @async run(`"paraview.exe" $File`)
    rm(File)
    VTKWrite.vtkwrite(File, fes.conn, geom.values, VTKWrite.Q4);
    # @async run(`"paraview.exe" $File`)
    rm(File)

    # println("Done")
    true
end
end
using .mmmLshapemmm_2
0mmmLshapemmm_2.test()

module mmsmoothingm1_2
using FinEtools
using FinEtools.MeshExportModule: VTKWrite
using Test
import LinearAlgebra: norm
function test()

    # println("""
    # Meshing, deforming  and smoothing
    # """)

    A= 100. # strip width
    N = 16
    tolerance = A/N/1.0e5
    fens,fes = T3block(A, A, N, N)

    bnl = connectednodes(meshboundary(fes))
    for ixxxx = 1:length(bnl)
        x, y = fens.xyz[bnl[ixxxx], :]
        fens.xyz[bnl[ixxxx], 1] += A/N*sin(2*pi*y/A)
        fens.xyz[bnl[ixxxx], 2] += -A/N*sin(2*pi*x/A)
    end

    File =  "mesh_smoothing_before.vtu"
    VTKWrite.vtkwrite(File, fens, fes);
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end
    before = [50.0, 43.75]
    # println("$(fens.xyz[Int(N^2/2), :] )")
    @test norm(before-fens.xyz[Int(N^2/2), :]) < 1.e-4

    fixedv = falses(count(fens))
    fixedv[bnl] .= true
    fens = meshsmoothing(fens, fes; fixedv = fixedv, method = :laplace, npass = 100)

    after = [50.0438, 44.0315]
    # println("$(fens.xyz[Int(N^2/2), :] )")
    @test norm(fens.xyz[Int(N^2/2), :]-after) < 1.e-4

    geom = NodalField(fens.xyz)

    File =  "mesh_smoothing_after.vtu"
    VTKWrite.vtkwrite(File, connasarray(fes), geom.values, VTKWrite.T3);
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    # println("Done")
    true
end
end
using .mmsmoothingm1_2
mmsmoothingm1_2.test()

module mmAbaqusexport3_2
using FinEtools
using FinEtools.MeshImportModule
using FinEtools.MeshExportModule: VTKWrite
using Test
function test()
    ## Solid cylinder/taper/sphere—-temperature loading; quadratic brick mesh

    # The mesh  will be created in a very coarse representation from the
    # key points in the drawing. The first coordinate is radial, the second coordinate is axial.
    rz=[1.     0.;#A
    1.4    0.;#B
    0.995184726672197   0.098017140329561;
    1.393258617341076 0.137223996461385;
    0.980785  0.195090;#
    1.37309939 0.27312645;
    0.956940335732209   0.290284677254462
    1.339716470025092 0.406398548156247
    0.9238795  0.38268;#C
    1.2124  0.7;#D
    0.7071  0.7071;#E
    1.1062  1.045;#F
    0.7071  (0.7071+1.79)/2;#(E+H)/2
    1.      1.39;#G
    0.7071  1.79;#H
    1.      1.79;#I
    ]*phun("M")
    tolerance =1.e-6*phun("M")

    # This is the quadrilateral mesh of the cross-section.   It will be modified and
    # refined as  we go.
    fens = FENodeSet(rz);
    fes = FESetQ4([1 2 4 3; 3 4 6 5; 5 6 8 7; 7 8 10 9; 9 10 12 11; 11 12 14 13; 13 14 16 15]);

    nref = 1;
    for ref = 1:nref
        fens,fes = Q4refine(fens,fes);
        list = selectnode(fens, distance=1.0+0.1/2^nref, from=[0. 0.], inflate=tolerance);
        fens.xyz[list,:] = FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);
    end

    ##
    # The mesh is extruded by sweeping around the axis of symmetry.
    # Only a single layer of elements is generated of internal angle
    # |angslice|.
    nLayers = 7;
    angslice  = 5*pi/16;

    ##
    # First the mesh is extruded to a block whose third dimension
    # represents the angular coordinate.
    fens,fes = H8extrudeQ4(fens, fes, nLayers,
    (rz,k)->[rz[1],rz[2],0.0]-(k)/nLayers*[0.,0.,angslice]);
    # The block is now converted  to the axially symmetric geometry by using the
    # third (angular) coordinate  to sweep out  an axially symmetric domain. The
    # ccoordinates of the nodes at this point are |rza|,  radial distance,
    # Z-coordinate, angle.
    sweep(rza) = [-rza[1]*sin(rza[3]+angslice/2.0), rza[1]*cos(rza[3]+angslice/2.0), rza[2]]
    for j=1:size(fens.xyz,1)
        fens.xyz[j,:] = sweep(fens.xyz[j,:])
    end

    AE = AbaqusExporter("LE11NAFEMS_H8_B.inp");
    HEADING(AE, "LE11NAFEMS: Linear bricks.");
    PART(AE, "part1");
    END_PART(AE);
    ASSEMBLY(AE, "ASSEM1");
    INSTANCE(AE, "INSTNC1", "PART1");
    NODE(AE, fens.xyz);
    ELEMENT(AE, "c3d8rh", "AllElements", 1, connasarray(fes))
    END_INSTANCE(AE);
    END_ASSEMBLY(AE);
    STEP_PERTURBATION_BUCKLE(AE, 1)
    CLOAD(AE, 1, 2, 10.0)
    EL_PRINT(AE, "AllElements", "S")
    ENERGY_PRINT(AE);
    END_STEP(AE);
    close(AE)

    output = MeshImportModule.import_ABAQUS("./LE11NAFEMS_H8_B.inp";
    allocationchunk = 11)
    fens, fes = output["fens"], output["fesets"][1]

    File = "LE11NAFEMS_H8.vtu"
    VTKWrite.vtkwrite(File, fens, fes)
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    try rm(AE.filename) catch end
end
end
using .mmAbaqusexport3_2
mmAbaqusexport3_2.test()

module momap2para3_2a
using FinEtools
using FinEtools.MeshSelectionModule: vselect
using FinEtools.MeshExportModule
using FinEtools.MeshExportModule: VTKWrite
using Test
import LinearAlgebra: norm
function test()
    A = 50.0*phun("m") # length  of loaded rectangle
    B = 200.0*phun("m") # length  of loaded rectangle
    C = 100.0*phun("m") # span of the plate

    # Select how find the mesh should be
    Refinement = 2
    nA, nB, nC = Refinement * 1, Refinement * 2, Refinement * 4;
    xs = reshape(collect(linearspace(0.0, A, nA + 1)), nA + 1, 1)
    ys = reshape(collect(linearspace(0.0, B, nB + 1)), nB + 1, 1)
    zs = reshape(collect(linearspace(0.0, C, nC + 1)), nC + 1, 1)
    fensc,fesc = T10blockx(xs, ys, zs, :b)
    fc = NodalField(zeros(count(fensc), 1))
    for i = 1:count(fensc)
        x, y, z = fensc.xyz[i, :]
        fc.values[i, :] .= sin(2*x/A) * cos(6.5*y/B) * sin(3*z/C-1.0)
    end
    File = "momap2para3_2a-coarse.vtu"
    VTKWrite.vtkwrite(File, fensc, fesc; scalars = [("fc", fc.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    Refinement = Refinement + 1
    nA, nB, nC = Refinement * 1, Refinement * 2, Refinement * 4;
    xs = reshape(collect(linearspace(0.0, A, nA + 1)), nA + 1, 1)
    ys = reshape(collect(linearspace(0.0, B, nB + 1)), nB + 1, 1)
    zs = reshape(collect(linearspace(0.0, C, nC + 1)), nC + 1, 1)
    fensf,fesf = T10blockx(xs, ys, zs, :b)
    ff = NodalField(zeros(count(fensf), 1))
    tolerance = min(A/nA, B/nB, C/nC)/1000.0

    referenceff = NodalField(zeros(count(fensf), 1))
    for i = 1:count(fensf)
        x, y, z = fensf.xyz[i, :]
        referenceff.values[i, :] .= sin(2*x/A) * cos(6.5*y/B) * sin(3*z/C-1.0)
    end
    File = "momap2para3_2a-reference.vtu"
    VTKWrite. vtkwrite(File, fensf, fesf; scalars = [("referenceff", referenceff.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    ff = transferfield!(ff, fensf, fesf, fc, fensc, fesc, tolerance)
    File = "momap2para3_2a-fine.vtu"
    VTKWrite.vtkwrite(File, fensf, fesf; scalars = [("ff", ff.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    diffff = NodalField(referenceff.values - ff.values)
    femm  = FEMMBase(IntegDomain(fesf, SimplexRule(3, 4)))
    geom = NodalField(fensf.xyz)
    error = integratefieldfunction(femm, geom, diffff, (x, v) -> norm(v), 0.0)
    ref = integratefieldfunction(femm, geom, referenceff, (x, v) -> norm(v), 0.0)
    # println("error/ref = $(error/ref)")
    @test abs(error/ref - 0.02541369940759616) < 1.0e-4
end
end
using .momap2para3_2a
momap2para3_2a.test()

module momap2para4_2
using FinEtools
using FinEtools.MeshSelectionModule: vselect
using FinEtools.MeshExportModule
using FinEtools.MeshExportModule: VTKWrite
using Test
import LinearAlgebra: norm
function test()
    A = 50.0*phun("m") # length  of loaded rectangle
    B = 200.0*phun("m") # length  of loaded rectangle
    C = 100.0*phun("m") # span of the plate

    # Select how find the mesh should be
    Refinement = 2
    nA, nB, nC = Refinement * 1, Refinement * 6, Refinement * 4;
    xs = reshape(collect(linearspace(0.0, A, nA + 1)), nA + 1, 1)
    ys = reshape(collect(linearspace(0.0, B, nB + 1)), nB + 1, 1)
    zs = reshape(collect(linearspace(0.0, C, nC + 1)), nC + 1, 1)
    fensc,fesc = T10blockx(xs, ys, zs, :b)

    centroidpc = centroidparametric(fesc)
    N = bfun(fesc, centroidpc)
    NT = transpose(N)

    fc = ElementalField(zeros(count(fesc), 1))
    for i = 1:count(fesc)
        c = [k for k in fesc.conn[i]]
        centroid = NT * fensc.xyz[c, :]
        x, y, z = centroid
        fc.values[i, :] .= sin(2*x/A) * cos(6.5*y/B) * sin(3*z/C-1.0)
    end
    File = "momap2para3-coarse.vtu"
    VTKWrite.vtkwrite(File, fensc, fesc; scalars = [("fc", fc.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    Refinement = Refinement + 1
    nA, nB, nC = Refinement * 1, Refinement * 6, Refinement * 4;
    xs = reshape(collect(linearspace(0.0, A, nA + 1)), nA + 1, 1)
    ys = reshape(collect(linearspace(0.0, B, nB + 1)), nB + 1, 1)
    zs = reshape(collect(linearspace(0.0, C, nC + 1)), nC + 1, 1)
    fensf,fesf = T10blockx(xs, ys, zs, :b)
    tolerance = min(A/nA, B/nB, C/nC)/1000.0

    ff = ElementalField(zeros(count(fesf), 1))
    referenceff = ElementalField(zeros(count(fesf), 1))
    for i = 1:count(fesf)
        c = [k for k in fesf.conn[i]]
        centroid = NT * fensf.xyz[c, :]
        x, y, z = centroid
        referenceff.values[i, :] .= sin(2*x/A) * cos(6.5*y/B) * sin(3*z/C-1.0)
    end
    File = "momap2para3-reference.vtu"
    VTKWrite.vtkwrite(File, fensf, fesf; scalars = [("referenceff", referenceff.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    ff = transferfield!(ff, fensf, fesf, fc, fensc, fesc, tolerance)
    File = "momap2para3-fine.vtu"
    VTKWrite.vtkwrite(File, fensf, fesf; scalars = [("ff", ff.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    diffff = ElementalField(referenceff.values - ff.values)
    femm  = FEMMBase(IntegDomain(fesf, SimplexRule(3, 4)))
    geom = NodalField(fensf.xyz)
    error = integratefieldfunction(femm, geom, diffff, (x, v) -> norm(v), 0.0)
    ref = integratefieldfunction(femm, geom, referenceff, (x, v) -> norm(v), 0.0)
    # println("error/ref = $(error/ref)")
    @test abs(error/ref - 0.19808425992688541) < 1.0e-4
end
end
using .momap2para4_2
momap2para4_2.test()

module momap2para5_2
using FinEtools
using FinEtools.MeshSelectionModule: vselect
using FinEtools.MeshExportModule
using FinEtools.MeshExportModule: VTKWrite
using Test
import LinearAlgebra: norm
function test()
    A = 50.0*phun("m") # length  of loaded rectangle
    B = 200.0*phun("m") # length  of loaded rectangle
    C = 100.0*phun("m") # span of the plate

    Meshing = T4blockx
    # Select how find the mesh should be
    Refinement = 2
    nA, nB, nC = Refinement * 1, Refinement * 6, Refinement * 4;
    xs = reshape(collect(linearspace(0.0, A, nA + 1)), nA + 1, 1)
    ys = reshape(collect(linearspace(0.0, B, nB + 1)), nB + 1, 1)
    zs = reshape(collect(linearspace(0.0, C, nC + 1)), nC + 1, 1)
    fensc,fesc = Meshing(xs, ys, zs, :b)

    centroidpc = centroidparametric(fesc)
    N = bfun(fesc, centroidpc)
    NT = transpose(N)

    fc = ElementalField(zeros(count(fesc), 1))
    for i = 1:count(fesc)
        c = [k for k in fesc.conn[i]]
        centroid = NT * fensc.xyz[c, :]
        x, y, z = centroid
        fc.values[i, :] .= sin(2*x/A) * cos(6.5*y/B) * sin(3*z/C-1.0)
    end
    File = "momap2para3-coarse.vtu"
    VTKWrite.vtkwrite(File, fensc, fesc; scalars = [("fc", fc.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    Refinement = Refinement + 1
    nA, nB, nC = Refinement * 1, Refinement * 6, Refinement * 4;
    xs = reshape(collect(linearspace(0.0, A, nA + 1)), nA + 1, 1)
    ys = reshape(collect(linearspace(0.0, B, nB + 1)), nB + 1, 1)
    zs = reshape(collect(linearspace(0.0, C, nC + 1)), nC + 1, 1)
    fensf,fesf = Meshing(xs, ys, zs, :b)
    tolerance = min(A/nA, B/nB, C/nC)/1000.0

    ff = ElementalField(zeros(count(fesf), 1))
    referenceff = ElementalField(zeros(count(fesf), 1))
    for i = 1:count(fesf)
        c = [k for k in fesf.conn[i]]
        centroid = NT * fensf.xyz[c, :]
        x, y, z = centroid
        referenceff.values[i, :] .= sin(2*x/A) * cos(6.5*y/B) * sin(3*z/C-1.0)
    end
    File = "momap2para3-reference.vtu"
    VTKWrite.vtkwrite(File, fensf, fesf; scalars = [("referenceff", referenceff.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    ff = transferfield!(ff, fensf, fesf, fc, fensc, fesc, tolerance)
    File = "momap2para3-fine.vtu"
    VTKWrite.vtkwrite(File, fensf, fesf; scalars = [("ff", ff.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    diffff = ElementalField(referenceff.values - ff.values)
    femm  = FEMMBase(IntegDomain(fesf, SimplexRule(3, 4)))
    geom = NodalField(fensf.xyz)
    error = integratefieldfunction(femm, geom, diffff, (x, v) -> norm(v), 0.0)
    ref = integratefieldfunction(femm, geom, referenceff, (x, v) -> norm(v), 0.0)
    # println("error/ref = $(error/ref)")
    @test abs(error/ref - 0.19808425992688541) < 1.0e-4
end
end
using .momap2para5_2
momap2para5_2.test()

module momap2para6_2
using FinEtools
using FinEtools.MeshSelectionModule: vselect
using FinEtools.MeshExportModule
using FinEtools.MeshExportModule: VTKWrite
using Test
import LinearAlgebra: norm
function test()
    A = 50.0*phun("m") # length  of loaded rectangle
    B = 200.0*phun("m") # length  of loaded rectangle
    C = 100.0*phun("m") # span of the plate

    Meshing = H20blockx
    # Select how find the mesh should be
    Refinement = 2
    nA, nB, nC = Refinement * 1, Refinement * 6, Refinement * 4;
    xs = collect(linearspace(0.0, A, nA + 1))
    ys = collect(linearspace(0.0, B, nB + 1))
    zs = collect(linearspace(0.0, C, nC + 1))
    fensc,fesc = Meshing(xs, ys, zs)

    centroidpc = centroidparametric(fesc)
    N = bfun(fesc, centroidpc)
    NT = transpose(N)

    fc = ElementalField(zeros(count(fesc), 1))
    for i = 1:count(fesc)
        c = [k for k in fesc.conn[i]]
        centroid = NT * fensc.xyz[c, :]
        x, y, z = centroid
        fc.values[i, :] .= sin(2*x/A) * cos(6.5*y/B) * sin(3*z/C-1.0)
    end
    File = "momap2para3-coarse.vtu"
    VTKWrite.vtkwrite(File, fensc, fesc; scalars = [("fc", fc.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    Refinement = Refinement + 1
    nA, nB, nC = Refinement * 1, Refinement * 6, Refinement * 4;
    xs = collect(linearspace(0.0, A, nA + 1))
    ys = collect(linearspace(0.0, B, nB + 1))
    zs = collect(linearspace(0.0, C, nC + 1))
    fensf,fesf = Meshing(xs, ys, zs)
    tolerance = min(A/nA, B/nB, C/nC)/1000.0

    ff = ElementalField(zeros(count(fesf), 1))
    referenceff = ElementalField(zeros(count(fesf), 1))
    for i = 1:count(fesf)
        c = [k for k in fesf.conn[i]]
        centroid = NT * fensf.xyz[c, :]
        x, y, z = centroid
        referenceff.values[i, :] .= sin(2*x/A) * cos(6.5*y/B) * sin(3*z/C-1.0)
    end
    File = "momap2para3-reference.vtu"
    VTKWrite.vtkwrite(File, fensf, fesf; scalars = [("referenceff", referenceff.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    ff = transferfield!(ff, fensf, fesf, fc, fensc, fesc, tolerance)
    File = "momap2para3-fine.vtu"
    VTKWrite.vtkwrite(File, fensf, fesf; scalars = [("ff", ff.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    diffff = ElementalField(referenceff.values - ff.values)
    femm  = FEMMBase(IntegDomain(fesf, GaussRule(3, 3)))
    geom = NodalField(fensf.xyz)
    error = integratefieldfunction(femm, geom, diffff, (x, v) -> norm(v), 0.0)
    ref = integratefieldfunction(femm, geom, referenceff, (x, v) -> norm(v), 0.0)
    # println("error/ref = $(error/ref)")
    @test abs(error/ref - 0.3602708839379691) < 1.0e-4
end
end
using .momap2para6_2
momap2para6_2.test()

module momap2para12_2
using FinEtools
using FinEtools.MeshSelectionModule: vselect
using FinEtools.MeshExportModule
using FinEtools.MeshExportModule: VTKWrite
using Test
import LinearAlgebra: norm
function test()
    A = 50.0*phun("m") # length  of loaded rectangle
    B = 200.0*phun("m") # length  of loaded rectangle
    C = 100.0*phun("m") # span of the plate

    Meshing = T6blockx
    # Select how find the mesh should be
    Refinement = 2
    nA, nB, nC = Refinement * 1, Refinement * 6, Refinement * 4;
    xs = collect(linearspace(0.0, A, nA + 1))
    ys = collect(linearspace(0.0, B, nB + 1))
    fensc,fesc = Meshing(xs, ys)

    centroidpc = centroidparametric(fesc)
    N = bfun(fesc, centroidpc)
    NT = transpose(N)

    fc = ElementalField(zeros(count(fesc), 1))
    for i = 1:count(fesc)
        c = [k for k in fesc.conn[i]]
        centroid = NT * fensc.xyz[c, :]
        x, y = centroid
        fc.values[i, :] .= sin(2*x/A) * cos(6.5*y/B)
    end
    File = "momap2para12_2-coarse.vtu"
    VTKWrite.vtkwrite(File, fensc, fesc; scalars = [("fc", fc.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    Refinement = Refinement + 1
    nA, nB, nC = Refinement * 1, Refinement * 6, Refinement * 4;
    xs = collect(linearspace(0.0, A, nA + 1))
    ys = collect(linearspace(0.0, B, nB + 1))
    zs = collect(linearspace(0.0, C, nC + 1))
    fensf,fesf = Meshing(xs, ys)
    tolerance = min(A/nA, B/nB, C/nC)/1000.0

    ff = ElementalField(zeros(count(fesf), 1))
    referenceff = ElementalField(zeros(count(fesf), 1))
    for i = 1:count(fesf)
        c = [k for k in fesf.conn[i]]
        centroid = NT * fensf.xyz[c, :]
        x, y = centroid
        referenceff.values[i, :] .= sin(2*x/A) * cos(6.5*y/B)
    end
    File = "momap2para12_2-reference.vtu"
    VTKWrite.vtkwrite(File, fensf, fesf; scalars = [("referenceff", referenceff.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    ff = transferfield!(ff, fensf, fesf, fc, fensc, fesc, tolerance)
    File = "momap2para12_2-fine.vtu"
    VTKWrite.vtkwrite(File, fensf, fesf; scalars = [("ff", ff.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    diffff = ElementalField(referenceff.values - ff.values)
    femm  = FEMMBase(IntegDomain(fesf, SimplexRule(2, 3)))
    geom = NodalField(fensf.xyz)
    error = integratefieldfunction(femm, geom, diffff, (x, v) -> norm(v), 0.0)
    ref = integratefieldfunction(femm, geom, referenceff, (x, v) -> norm(v), 0.0)
    # println("error/ref = $(error/ref)")
    @test abs(error/ref -  0.19654730310399787) < 1.0e-4
end
end
using .momap2para12_2
momap2para12_2.test()

module mt4refine1_2
using FinEtools
using FinEtools.MeshExportModule
using FinEtools.MeshExportModule: VTKWrite
using Test
function test()
    xs = collect(linearspace(0.0, pi / 2, 5))
    ys = collect(linearspace(0.0, 1.0, 6).^2)
    zs = collect(linearspace(0.0, 1.0, 3))
    fens, fes = T4blockx(xs, ys, zs, :a)
    # for i = 1:count(fens)
    #     a, y, z = fens.xyz[i,:]
    #     fens.xyz[i,1] = sin(a) * (y + 0.5)
    #     fens.xyz[i,2] = cos(a) * (y + 0.5)
    #     fens.xyz[i,3] = z
    # end
    @test count(fes) == 240
    bfes = meshboundary(fes)
    @test count(bfes) == 2*2*(4*5 + 5*2 + 4*2)
    fens, fes = T4refine(fens, fes)
    for i = 1:count(fens)
        a, y, z = fens.xyz[i,:]
        fens.xyz[i,1] = sin(a) * (y + 0.5)
        fens.xyz[i,2] = cos(a) * (y + 0.5)
        fens.xyz[i,3] = z
    end
    @test count(fes) == 240*8
    bfes = meshboundary(fes)
    @test count(bfes) == 4*2*2*(4*5 + 5*2 + 4*2)

    geom  =  NodalField(fens.xyz)
    femm  =  FEMMBase(IntegDomain(fes, SimplexRule(3, 4)))
    V = integratefunction(femm, geom, (x) ->  1.0)

    File = "Refine-T4-a.vtu"
    VTKWrite.vtkwrite(File, fens, bfes)
    rm(File)
    # @async run(`"paraview.exe" $File`)
end
end
using .mt4refine1_2
mt4refine1_2.test()

module mimportexportm1_2
using FinEtools
using FinEtools.MeshImportModule
using FinEtools.MeshExportModule
using FinEtools.MeshExportModule: VTKWrite
using Test
function test()
    output = MeshImportModule.import_NASTRAN(dirname(@__FILE__) * "/" * "cylinder.nas";    allocationchunk = 13, expectfixedformat = true)
    @test count(output["fens"]) == 1406
    @test count(output["fesets"][1]) == 829
    # show(fes.conn[count(fes), :])
    File = "cylinder.vtu"
    VTKWrite.vtkwrite(File, output["fens"], output["fesets"][1])
    rm(File)
    #   @test output["fesets"][1].conn[count(output["fesets"][1]), :] == NTuple{10,Int64}[(143, 140, 144, 138, 361, 363, 176, 519, 781, 520)]
    # @async run(`"paraview.exe" $File`)
end
end
using .mimportexportm1_2
mimportexportm1_2.test()

module momap2para61_2
using FinEtools
using FinEtools.MeshSelectionModule: vselect
using FinEtools.MeshExportModule
using FinEtools.MeshExportModule: VTKWrite
using Test
import LinearAlgebra: norm
function test()
    A = 50.0*phun("m") # length  of loaded rectangle
    B = 200.0*phun("m") # length  of loaded rectangle
    C = 100.0*phun("m") # span of the plate

    # Select how find the mesh should be
    Refinement = 2
    nA, nB, nC = Refinement * 1, Refinement * 2, Refinement * 4;
    xs = reshape(collect(linearspace(0.0, A, nA + 1)), nA + 1, 1)
    ys = reshape(collect(linearspace(0.0, B, nB + 1)), nB + 1, 1)
    zs = reshape(collect(linearspace(0.0, C, nC + 1)), nC + 1, 1)
    fensc,fesc = T10blockx(xs, ys, zs, :b)
    fc = NodalField(zeros(count(fensc), 1))
    for i = 1:count(fensc)
        x, y, z = fensc.xyz[i, :]
        fc.values[i, :] .= sin(2*x/A) * cos(6.5*y/B) * sin(3*z/C-1.0)
    end
    File = "momap2para61_2-coarse.vtu"
    VTKWrite.vtkwrite(File, fensc, fesc; scalars = [("fc", fc.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    Refinement = Refinement + 1
    nA, nB, nC = Refinement * 1, Refinement * 2, Refinement * 4;
    xs = reshape(collect(linearspace(0.0, A, nA + 1)), nA + 1, 1)
    ys = reshape(collect(linearspace(0.0, B, nB + 1)), nB + 1, 1)
    zs = reshape(collect(linearspace(0.0, C, nC + 1)), nC + 1, 1)
    fensf,fesf = T10blockx(xs, ys, zs, :b)
    ff = NodalField(zeros(count(fensf), 1))
    tolerance = min(A/nA, B/nB, C/nC)/1000.0

    referenceff = NodalField(zeros(count(fensf), 1))
    for i = 1:count(fensf)
        x, y, z = fensf.xyz[i, :]
        referenceff.values[i, :] .= sin(2*x/A) * cos(6.5*y/B) * sin(3*z/C-1.0)
    end
    File = "momap2para61_2-reference.vtu"
    VTKWrite.vtkwrite(File, fensf, fesf; scalars = [("referenceff", referenceff.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    ff = transferfield!(ff, fensf, fesf, fc, fensc, fesc, tolerance)
    File = "momap2para61_2-fine.vtu"
    VTKWrite.vtkwrite(File, fensf, fesf; scalars = [("ff", ff.values), ("ffcopy", ff.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    diffff = NodalField(referenceff.values - ff.values)
    femm  = FEMMBase(IntegDomain(fesf, SimplexRule(3, 4)))
    geom = NodalField(fensf.xyz)
    error = integratefieldfunction(femm, geom, diffff, (x, v) -> norm(v), 0.0)
    ref = integratefieldfunction(femm, geom, referenceff, (x, v) -> norm(v), 0.0)
    # println("error/ref = $(error/ref)")
    @test abs(error/ref - 0.02541369940759616) < 1.0e-4
end
end
using .momap2para61_2
momap2para61_2.test()

module momap2para6378_2
using FinEtools
using FinEtools.MeshSelectionModule: vselect
using FinEtools.MeshExportModule
using FinEtools.MeshExportModule: VTKWrite
using Test
import LinearAlgebra: norm
function test()
    A = 50.0*phun("m") # length  of loaded rectangle
    B = 200.0*phun("m") # length  of loaded rectangle
    C = 100.0*phun("m") # span of the plate

    # Select how find the mesh should be
    Refinement = 2
    nA, nB, nC = Refinement * 1, Refinement * 6, Refinement * 4;
    xs = reshape(collect(linearspace(0.0, A, nA + 1)), nA + 1, 1)
    ys = reshape(collect(linearspace(0.0, B, nB + 1)), nB + 1, 1)
    zs = reshape(collect(linearspace(0.0, C, nC + 1)), nC + 1, 1)
    fensc,fesc = T10blockx(xs, ys, zs, :b)

    centroidpc = centroidparametric(fesc)
    N = bfun(fesc, centroidpc)
    NT = transpose(N)

    fc = ElementalField(zeros(count(fesc), 1))
    for i = 1:count(fesc)
        c = [k for k in fesc.conn[i]]
        centroid = NT * fensc.xyz[c, :]
        x, y, z = centroid
        fc.values[i, :] .= sin(2*x/A) * cos(6.5*y/B) * sin(3*z/C-1.0)
    end
    File = "momap2para3-coarse.vtu"
    VTKWrite.vtkwrite(File, fensc, fesc; scalars = [("fc", fc.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    Refinement = Refinement + 1
    nA, nB, nC = Refinement * 1, Refinement * 6, Refinement * 4;
    xs = reshape(collect(linearspace(0.0, A, nA + 1)), nA + 1, 1)
    ys = reshape(collect(linearspace(0.0, B, nB + 1)), nB + 1, 1)
    zs = reshape(collect(linearspace(0.0, C, nC + 1)), nC + 1, 1)
    fensf,fesf = T10blockx(xs, ys, zs, :b)
    tolerance = min(A/nA, B/nB, C/nC)/1000.0

    ff = ElementalField(zeros(count(fesf), 1))
    referenceff = ElementalField(zeros(count(fesf), 1))
    for i = 1:count(fesf)
        c = [k for k in fesf.conn[i]]
        centroid = NT * fensf.xyz[c, :]
        x, y, z = centroid
        referenceff.values[i, :] .= sin(2*x/A) * cos(6.5*y/B) * sin(3*z/C-1.0)
    end
    File = "momap2para3-reference.vtu"
    VTKWrite.vtkwrite(File, fensf, fesf; scalars = [("referenceff", referenceff.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    ff = transferfield!(ff, fensf, fesf, fc, fensc, fesc, tolerance)
    File = "momap2para3-fine.vtu"
    VTKWrite.vtkwrite(File, fensf, fesf; scalars = [("ff", ff.values), ("ffcopy", ff.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    diffff = ElementalField(referenceff.values - ff.values)
    femm  = FEMMBase(IntegDomain(fesf, SimplexRule(3, 4)))
    geom = NodalField(fensf.xyz)
    error = integratefieldfunction(femm, geom, diffff, (x, v) -> norm(v), 0.0)
    ref = integratefieldfunction(femm, geom, referenceff, (x, v) -> norm(v), 0.0)
    # println("error/ref = $(error/ref)")
    @test abs(error/ref - 0.19808425992688541) < 1.0e-4
end
end
using .momap2para6378_2
momap2para6378_2.test()

module mcircle3
using FinEtools
using FinEtools.MeshExportModule: VTK
using Test
function test()
   fens,fes = T3circlen(1.3, 5); # Mesh
   # File = "mesh.vtk"
   #   VTK.vtkexportmesh(File, fens, fes)
   l = selectnode(fens, box = [0.0, 1.3, 0.0, 0.0], inflate = 0.01)
   @test l == [1, 3, 7, 12, 16, 28, 32, 44, 49]
   true
end
end
using .mcircle3
mcircle3.test()

module mcircle3a
using FinEtools
using FinEtools.MeshExportModule: VTK
using Test
function test()
    r = 1.3
    m, n = 8, 4
    fens,fes = T3block(pi/2, r, m, n); # Mesh
    for i in 1:count(fens)
        a = pi/2 - fens.xyz[i, 1]
        r = fens.xyz[i, 2]
        fens.xyz[i, :] .= (r*cos(a), r*sin(a))
    end
    fens, newn = fusenodes(fens, fens, r/1000)
    updateconn!(fes, newn)
    connected = findunconnnodes(fens, fes);
    fens, newn = compactnodes(fens, connected);
    @test count(fens) == (m+1)*(n+1) - (m+1) + 1
    fes = renumberconn!(fes, newn);
    keep = fill(true, count(fes))
    ca = connasarray(fes)
    for j in 1:count(fes)
        keep[j] = (length(unique(ca[j, :])) == 3)
    end
    l =  findall(x -> x == true, keep)
    fes = subset(fes, l)
    @test count(fes) == 2*m*n - m
    # File = "mesh.vtk"
    # VTK.vtkexportmesh(File, fens, fes)
    bfes = meshboundary(fes)
    # File = "mesh-boundary.vtk"
    # VTK.vtkexportmesh(File, fens, bfes)
    @test count(bfes) == m+2*n
    true
end
end
using .mcircle3a
mcircle3a.test()

module mcircle3b
using FinEtools
using FinEtools.MeshExportModule: VTK
using Test
function test()
    r = 1.3
    nc, nr = 8, 4
    fens,fes = T3circleseg(pi/2, r, nc, nr); # Mesh
    # File = "mesh.vtk"
    # VTK.vtkexportmesh(File, fens, fes)
    bfes = meshboundary(fes)
    # File = "mesh-boundary.vtk"
    # VTK.vtkexportmesh(File, fens, bfes)
    @test count(fens) == (nc+1)*(nr+1) - (nc+1) + 1
    @test count(fes) == 2*nc*nr - nc
    @test count(bfes) == nc+2*nr
    true
end
end
using .mcircle3b
mcircle3b.test()

module ml2extrude1
using FinEtools
using FinEtools.MeshExportModule: VTK
using Test
function test()
    X = [0.0 0.0;
        0.0 -0.2;
        0.3 -0.2;
        0.3 -0.4
    ]
    l2fens = FENodeSet(X);
    conn = [1 2; 2 3; 3 4]
    l2fes = FESetL2(conn);
    # fens,fes = Q4extrudeL2(l2fens, l2fes, 5, (x, k) -> vcat(vec(x), (k-1)*1)); # Mesh
    ex = function (x, k)
        x .+= k*0.1
        x
    end
    fens,fes = Q4extrudeL2(l2fens, l2fes, 5, ex); # Mesh
    # File = "mesh.vtk"
    # VTK.vtkexportmesh(File, fens, fes)
    # bfes = meshboundary(fes)
    # File = "mesh-boundary.vtk"
    # VTK.vtkexportmesh(File, fens, bfes)
    @test count(fens) == 24
    @test count(fes) == 15
    # @test count(bfes) == nc+2*nr
    true
end
end
using .ml2extrude1
ml2extrude1.test()

module ml2extrude2
using FinEtools
using FinEtools.MeshExportModule: VTK
using Test
function test()
    X = [0.0 0.0 0.0;
        0.0 -0.2 0.0;
        0.3 -0.2 0.0;
        0.3 -0.4 0.0
    ]
    l2fens = FENodeSet(X);
    conn = [1 2; 2 3; 3 4]
    l2fes = FESetL2(conn);
    # fens,fes = Q4extrudeL2(l2fens, l2fes, 5, (x, k) -> vcat(vec(x), (k-1)*1)); # Mesh
    ex = function (x, k)
        x[3] += (k-1)*0.1
        x
    end
    fens,fes = Q4extrudeL2(l2fens, l2fes, 5, ex); # Mesh
    # File = "mesh.vtk"
    # VTK.vtkexportmesh(File, fens, fes)
    bfes = meshboundary(fes)
    # File = "mesh-boundary.vtk"
    # VTK.vtkexportmesh(File, fens, bfes)
    @test count(fens) == 24
    @test count(fes) == 15
    @test count(bfes) == 5+5+3+3
    true
end
end
using .ml2extrude2
ml2extrude2.test()


text:
module AllanDeviations

#
# Exports/
#
export AllanTauDescriptor, AllTaus, QuarterOctave, HalfOctave, Octave, HalfDecade, Decade

export allandev
export mallandev
export hadamarddev
export timedev
export totaldev
export mtie
#
# /Exports
#

#
# Types/
#
abstract type AllanTauDescriptor end
struct AllTaus <: AllanTauDescriptor end
struct QuarterOctave <: AllanTauDescriptor end
struct HalfOctave <: AllanTauDescriptor end
struct Octave <: AllanTauDescriptor end
struct HalfDecade <: AllanTauDescriptor end
struct Decade <: AllanTauDescriptor end
#
# /Types
#

#
# Helper Functions/
#
function frequencytophase(data::Array{T, 1}, rate::AbstractFloat) where T
	dt = 1 / rate
	n = length(data) + 1
	dataPrime = zeros(T, n)
	walkingSum = zero(T)
	@inbounds for i in 2:n #spare the first element so that the phase begins with zero
		walkingSum += data[i - 1]
		dataPrime[i] = walkingSum * dt
	end
	dataPrime
end

#tau-descriptor to m
function taudescription_to_m(::Type{AllTaus}, rate::AbstractFloat, n::Int)
	1:(n - 2)
end
function taudescription_to_m(::Type{Decade}, rate::AbstractFloat, n::Int)
	10 .^(0:Int(floor(log10(n))))
end
function taudescription_to_m(::Type{HalfDecade}, rate::AbstractFloat, n::Int)
	5 .^(0:Int(floor(log(5.0, n))))
end
function taudescription_to_m(::Type{Octave}, rate::AbstractFloat, n::Int)
	2 .^(0:Int(floor(log2(n))))
end
function taudescription_to_m(::Type{HalfOctave}, rate::AbstractFloat, n::Int)
	unique(Int.(floor.(
		1.5 .^(1:Int(floor(log(1.5, n))))
		)))
end
function taudescription_to_m(::Type{QuarterOctave}, rate::AbstractFloat, n::Int)
	unique(Int.(floor.(
		1.25 .^(1:Int(floor(log(1.25, n))))
		)))
end
#tau with custom log base value to m
function taudescription_to_m(taus::AbstractFloat, rate::AbstractFloat, n::Int)
	if taus <= 1.0
		error("Custom `taus`-log scale must be greater than 1.0")
	end
	unique(Int.(floor.(
		taus .^(0:Int(floor(log(taus, n))))
		)))
end
#tau with custom array to m
function taudescription_to_m(taus::Array{Float64}, rate::AbstractFloat, n::Int)
	m = unique(Int.(floor.(rate .* taus)))
	m[m .>= 1]
end
#tau with custom tau-length count
function taudescription_to_m(count::Integer, rate::AbstractFloat, n::Int)
	unique(Int.(floor.(
		1.125 .^(1:(log(1.125, n) / count):log(1.125, n))
		)))
end
#
# /Helper Functions
#

#
# Exported functions/
#

include("dev_allan.jl")
include("dev_mallan.jl")
include("dev_hadamard.jl")
include("dev_time.jl")
include("dev_total.jl")
include("dev_mtie.jl")

#
# /Exported functions
#

end # AllanDeviations

ext: jl


text:
using LabDataSources
using Test

@testset "LabDataSources.jl" begin
    # Write your tests here.
end


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

# -------------
# COMPOSITIONS
# -------------

"""
    alr(c)

Additive log-ratio transformation of composition `c`.
"""
function alr(c::Composition{D}) where {D}
  w = components(c) .+ eps()
  SVector(ntuple(i -> log(w[i] / w[D]), D-1))
end

"""
    alrinv(x)

Inverse alr transformation of coordinates `x`.
"""
alrinv(x::SVector{D,T}) where {D,T<:Real} =
  Composition(𝒞([exp.(x); SVector(one(T))]))

alrinv(x::AbstractVector) = alrinv(SVector{length(x)}(x))

# -------
# TABLES
# -------

"""
    ALR([refvar])

Additive log-ratio transform following the
[TableTransforms.jl](https://github.com/JuliaML/TableTransforms.jl)
interface.

Optionally, specify the reference variable `refvar` for the ratios.
Default to the last column of the input table.
"""
struct ALR <: LogRatio
  refvar::Union{Symbol,Nothing}
end

ALR() = ALR(nothing)

refvar(transform::ALR, vars) =
  isnothing(transform.refvar) ? last(vars) : transform.refvar

newvars(::ALR, n) = collect(n)[begin:end-1]

oldvars(::ALR, vars, rvar) = [collect(vars); rvar]

function applymatrix(::ALR, X)
  L = log.(X .+ eps())
  L[:,begin:end-1] .- L[:,end]
end

function revertmatrix(::ALR, Y)
  E = [exp.(Y) ones(size(Y,1))]
  mapslices(𝒞, E, dims=2)
end


text:
# A demonstration how to use the new SparsePSHDegreeM degree type
# see also ?ACE.RPI.SparsePSHDegreeM

using ACE, JuLIP
using ACE: z2i, i2z, order
#---

zTi = AtomicNumber(:Ti)
zAl = AtomicNumber(:Al)

# the weights for the n
Dn = Dict( "default" => 1.0,
            (zTi, zAl) => 2.0,    # weak interaction between species
            (zAl, zTi) => 2.0,    # weak interaction between species
            (zAl, zAl) => 1.2,    # slightly smaller basis for Al than for Ti
            (2, zTi, zAl) => 0.8  # but for 3-body we override this ...
         )

# the weights for the l
Dl = Dict( "default" => 1.5, )    # let's do nothing special here...
                                  # but same options are available for for n

# the degrees
Dd = Dict( "default" => 10,
           1 => 20,     # N = 1
           2 => 20,     # ...
           3 => 15,
           (2, zTi) => 25   # an extra push for the 3-body Ti basis
        )                   # (probably a dumb idea, just for illustration!)

Deg = ACE.RPI.SparsePSHDegreeM(Dn, Dl, Dd)

#---

# generate basis
# - note that degree is already incorporated into Deg
#   but we can still enlarge it e.g. by using maxdeg = 1.2, 1.5, 2.0, ...
basis = ACE.Utils.rpi_basis(species = [:Ti, :Al],
                              N = 5,
                              r0 = rnn(:Ti),
                              D = Deg,
                              maxdeg = 1)

#---
# analyse the basis a bit to see the effect this had

iAl = z2i(basis, AtomicNumber(:Al))
iTi = z2i(basis, AtomicNumber(:Ti))

println("The Ti site has more basis functions than the Al site")
@show length(basis.pibasis.inner[iAl])
@show length(basis.pibasis.inner[iTi])

specTi = collect(keys(basis.pibasis.inner[iTi].b2iAA))
specTi2 = specTi[ order.(specTi) .== 2 ]
specTi3 = specTi[ order.(specTi) .== 3 ]

println("The Ti-2N interaction has more basis functions than Ti-3N interaction")
@show length(specTi2)
@show length(specTi3)

specTi2_Ti = specTi2[ [ all(b.z == zTi for b in B.oneps) for B in specTi2 ] ]
specTi2_Al = specTi2[ [ all(b.z == zAl for b in B.oneps) for B in specTi2 ] ]

println("The Ti-Ti interaction has more basis functions than Ti-Al interaction")
@show length(specTi2_Ti)
@show length(specTi2_Al)


text:
"""
# ComplexVisual

A module to visualize holomorphic/meromorphic functions.

This is done by showing how sets are changed/deformed by functions
(holomorphic/meromorphic functions interpreted as transformations)
and complex phase portraits.

Often the sets (that are transformed) depend on one degree of freedom (which
may be changed by the mouse in the `ComplexVisualGtk` frontend).
"""
module ComplexVisual

using ComplexPortraits
using Colors
using Cairo
using Printf

include("monkeypatch.jl")

include("basic.jl")
include("rectangle.jl")
include("canvas.jl")
include("context.jl")
include("contextstyle.jl")
include("painter.jl")
include("eximage.jl")
include("layout.jl")
include("text.jl")
include("axis.jl")
include("decoration.jl")
include("scene.jl")
include("lrdomains.jl")
include("slider.jl")
include("winding.jl")
include("compcodomains.jl")

"""macro for importing the *huge* set of symbols."""
macro import_huge()
    quote
        @ComplexVisual.import_basic_huge
        @ComplexVisual.import_rectangle_huge
        @ComplexVisual.import_canvas_huge
        @ComplexVisual.import_context_huge
        @ComplexVisual.import_contextstyle_huge
        @ComplexVisual.import_painter_huge
        @ComplexVisual.import_eximage_huge
        @ComplexVisual.import_layout_huge
        @ComplexVisual.import_text_huge
        @ComplexVisual.import_axis_huge
        @ComplexVisual.import_decoration_huge
        @ComplexVisual.import_scene_huge
        @ComplexVisual.import_lrdomains_huge
        @ComplexVisual.import_slider_huge
        @ComplexVisual.import_winding_huge
        @ComplexVisual.import_compcodomains_huge
    end
end

end

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
function record(alg, cf, lconfig)
    logs = synclogger.SyncLogger(cf, lconfig)
    logs.connect()
    try
        alg(logs)
    catch e
        @show e
        # Catch keyboard interrupts if desired
    end
    logs.disconnect()
end

function log_posn(uri=_first_available())
    play(uri) do cf
        lconfig = logger.LogConfig(name="Position", period_in_ms=100)
        lconfig.add_variable("stateEstimate.x", "float")
        lconfig.add_variable("stateEstimate.y", "float")
        lconfig.add_variable("stateEstimate.z", "float")
        record(cf, lconfig) do logs
            for entry in logs
                data = entry[2]
                print("x = ", data["stateEstimate.x"], ", ")
                print("y = ", data["stateEstimate.y"], ", ")
                print("z = ", data["stateEstimate.z"], "\n")
            end
        end
    end
end

function log_quat(uri=_first_available())
    play(uri) do cf
        lconfig = logger.LogConfig(name="Quaternion", period_in_ms=100)
        lconfig.add_variable("stateEstimate.qx", "float")
        lconfig.add_variable("stateEstimate.qy", "float")
        lconfig.add_variable("stateEstimate.qz", "float")
        lconfig.add_variable("stateEstimate.qw", "float")
        record(cf, lconfig) do logs
            for entry in logs
                data = entry[2]
                print("qx = ", data["stateEstimate.qx"], ", ")
                print("qy = ", data["stateEstimate.qy"], ", ")
                print("qz = ", data["stateEstimate.qz"], ", ")
                print("qw = ", data["stateEstimate.qw"], "\n")
            end
        end
    end
end

function log_config(uri=_first_available())
    play(uri) do cf
        lconfig = logger.LogConfig(name="Configuration", period_in_ms=100)
        lconfig.add_variable("stateEstimateZ.x", "int16_t")
        lconfig.add_variable("stateEstimateZ.y", "int16_t")
        lconfig.add_variable("stateEstimateZ.z", "int16_t")
        lconfig.add_variable("stateEstimate.qx", "float")
        lconfig.add_variable("stateEstimate.qy", "float")
        lconfig.add_variable("stateEstimate.qz", "float")
        lconfig.add_variable("stateEstimate.qw", "float")
        record(cf, lconfig) do logs
            for entry in logs
                data = entry[2]
                print("x = ", data["stateEstimateZ.x"]/100.0, ", ")
                print("y = ", data["stateEstimateZ.y"]/100.0, ", ")
                print("z = ", data["stateEstimateZ.z"]/100.0, ", ")
                print("qx = ", data["stateEstimate.qx"], ", ")
                print("qy = ", data["stateEstimate.qy"], ", ")
                print("qz = ", data["stateEstimate.qz"], ", ")
                print("qw = ", data["stateEstimate.qw"], "\n")
            end
        end
    end
end

#  function plot_state()
#  end


text:
Base.summary(df::AbstractDataFrame) =
    @sprintf("%d×%d %s", size(df)..., nameof(typeof(df)))
Base.summary(io::IO, df::AbstractDataFrame) = print(io, summary(df))

"""
    DataFrames.ourstrwidth(io::IO, x::Any, buffer::IOBuffer, truncstring::Int)

Determine the number of characters that would be used to print a value.
"""
function ourstrwidth(io::IO, x::Any, buffer::IOBuffer, truncstring::Int)
    truncate(buffer, 0)
    ourshow(IOContext(buffer, :compact=>get(io, :compact, true)), x, truncstring)
    return textwidth(String(take!(buffer)))
end

function truncatestring(s::AbstractString, truncstring::Int)
    truncstring <= 0 && return s
    totalwidth = 0
    for (i, c) in enumerate(s)
        totalwidth += textwidth(c)
        if totalwidth > truncstring
            return first(s, i-1) * '…'
        end
    end
    return s
end

"""
    DataFrames.ourshow(io::IO, x::Any, truncstring::Int)

Render a value to an `IO` object compactly using print.
`truncstring` indicates the approximate number of text characters width to truncate
the output (if it is a non-positive value then no truncation is applied).
"""
function ourshow(io::IO, x::Any, truncstring::Int; styled::Bool=false)
    io_ctx = IOContext(io, :compact=>get(io, :compact, true), :typeinfo=>typeof(x))
    sx = sprint(print, x, context=io_ctx)
    sx = escape_string(sx, ()) # do not escape "
    sx = truncatestring(sx, truncstring)
    styled ? printstyled(io_ctx, sx, color=:light_black) : print(io_ctx, sx)
end

const SHOW_TABULAR_TYPES = Union{AbstractDataFrame, DataFrameRow, DataFrameRows,
                                 DataFrameColumns, GroupedDataFrame}

# workaround Julia 1.0 for Char
ourshow(io::IO, x::Char, truncstring::Int; styled::Bool=false) =
    ourshow(io, string(x), styled=styled, truncstring)

ourshow(io::IO, x::Nothing, truncstring::Int; styled::Bool=false) =
    ourshow(io, "", styled=styled, truncstring)
ourshow(io::IO, x::SHOW_TABULAR_TYPES, truncstring::Int; styled::Bool=false) =
    ourshow(io, summary(x), truncstring, styled=styled)

function ourshow(io::IO, x::Markdown.MD, truncstring::Int)
    r = repr(x)
    truncstring <= 0 && return chomp(truncstring)
    len = min(length(r, 1, something(findfirst(==('\n'), r), lastindex(r)+1)-1), truncstring)
    return print(io, len < length(r) - 1 ? first(r, len)*'…' : first(r, len))
end

# AbstractChar: https://github.com/JuliaLang/julia/pull/34730 (1.5.0-DEV.261)
# Irrational: https://github.com/JuliaLang/julia/pull/34741 (1.5.0-DEV.266)
if VERSION < v"1.5.0-DEV.261" || VERSION < v"1.5.0-DEV.266"
    function ourshow(io::IO, x::T, truncstring::Int) where T <: Union{AbstractChar, Irrational}
        io = IOContext(io, :compact=>get(io, :compact, true), :typeinfo=>typeof(x))
        show(io, x)
    end
end

# For most data frames, especially wide, columns having the same element type
# occur multiple times. batch_compacttype ensures that we compute string
# representation of a specific column element type only once and then reuse it.

function batch_compacttype(types::Vector{Any}, maxwidths::Vector{Int})
    @assert length(types) == length(maxwidths)
    cache = Dict{Any, String}()
    return map(types, maxwidths) do T, maxwidth
        get!(cache, T) do
            compacttype(T, maxwidth)
        end
    end
end

function batch_compacttype(types::Vector{Any}, maxwidth::Int)
    cache = Dict{Type, String}()
    return map(types) do T
        get!(cache, T) do
            compacttype(T, maxwidth)
        end
    end
end

"""
    compacttype(T::Type, maxwidth::Int=8, initial::Bool=true)

Return compact string representation of type `T`.

For displaying data frame we do not want string representation of type to be
longer than `maxwidth`. This function implements rules how type names are
cropped if they are longer than `maxwidth`.
"""
function compacttype(T::Type, maxwidth::Int)
    maxwidth = max(8, maxwidth)

    T === Any && return "Any"
    T === Missing && return "Missing"

    sT = string(T)
    textwidth(sT) ≤ maxwidth && return sT

    if T >: Missing
        T = nonmissingtype(T)
        sT = string(T)
        suffix = "?"
        textwidth(sT) ≤ maxwidth && return sT * suffix
    else
        suffix = ""
    end

    maxwidth -= 1 # we will add "…" at the end

    # This is only type display shortening so we
    # are OK with any T whose name starts with CategoricalValue here
    if startswith(sT, "CategoricalValue") || startswith(sT, "CategoricalArrays.CategoricalValue")
        sT = string(nameof(T))
        if textwidth(sT) ≤ maxwidth
            return sT * "…" * suffix
        else
            return (maxwidth ≥ 11 ? "Categorical…" : "Cat…") * suffix
        end
    elseif T isa Union
        return "Union…" * suffix
    else
        sTfull = sT
        sT = string(nameof(T))
    end

    # handle the case when the type printed is not parametric but string(T)
    # prefixed it with the module name which caused it to be overlong
    textwidth(sT) ≤ maxwidth + 1 && endswith(sTfull, sT) && return sT

    cumwidth = 0
    stop = 0
    for (i, c) in enumerate(sT)
        cumwidth += textwidth(c)
        if cumwidth ≤ maxwidth
            stop = i
        else
            break
        end
    end
    return first(sT, stop) * "…" * suffix
end

function _show(io::IO,
               df::AbstractDataFrame;
               allrows::Bool = !get(io, :limit, false),
               allcols::Bool = !get(io, :limit, false),
               rowlabel::Symbol = :Row,
               summary::Bool = true,
               eltypes::Bool = true,
               rowid = nothing,
               truncate::Int = 32,
               kwargs...)

    _check_consistency(df)

    names_str = names(df)
    names_len = Int[textwidth(n) for n in names_str]
    maxwidth = Int[max(9, nl) for nl in names_len]
    types = Any[eltype(c) for c in eachcol(df)]
    types_str = batch_compacttype(types, maxwidth)

    if allcols && allrows
        crop = :none
    elseif allcols
        crop = :vertical
    elseif allrows
        crop = :horizontal
    else
        crop = :both
    end

    # For consistency, if `kwargs` has `compact_printng`, we must use it.
    compact_printing::Bool = get(kwargs, :compact_printing, get(io, :compact, true))

    num_rows, num_cols = size(df)

    # By default, we align the columns to the left unless they are numbers,
    # which is checked in the following.
    alignment = fill(:l, num_cols)

    # Create the dictionary with the anchor regex that is used to align the
    # floating points.
    alignment_anchor_regex = Dict{Int, Vector{Regex}}()

    # Regex to align real numbers.
    alignment_regex_real = [r"\."]

    # Regex for columns with complex numbers.
    #
    # Here we are matching `+` or `-` unless it is not at the beginning of the
    # string or an `e` precedes it.
    alignment_regex_complex = [r"(?<!^)(?<!e)[+-]"]

    for i = 1:num_cols
        type_i = nonmissingtype(types[i])

        if type_i <: Complex
            alignment_anchor_regex[i] = alignment_regex_complex
            alignment[i] = :r
        elseif type_i <: Real
            alignment_anchor_regex[i] = alignment_regex_real
            alignment[i] = :r
        elseif type_i <: Number
            alignment[i] = :r
        end
    end

    # Make sure that `truncate` does not hide the type and the column name.
    maximum_columns_width = Int[truncate == 0 ? 0 : max(truncate + 1, l, textwidth(t))
                                for (l, t) in zip(names_len, types_str)]

    # Check if the user wants to display a summary about the DataFrame that is
    # being printed. This will be shown using the `title` option of
    # `pretty_table`.
    title = summary ? Base.summary(df) : ""

    # If `rowid` is not `nothing`, then we are printing a data row. In this
    # case, we will add this information using the row name column of
    # PrettyTables.jl. Otherwise, we can just use the row number column.
    if (rowid === nothing) || (ncol(df) == 0)
        show_row_number::Bool = get(kwargs, :show_row_number, true)
        row_names = nothing

        # If the columns with row numbers is not shown, then we should not
        # display a vertical line after the first column.
        vlines = fill(1, show_row_number)
    else
        nrow(df) != 1 &&
            throw(ArgumentError("rowid may be passed only with a single row data frame"))

        # In this case, if the user does not want to show the row number, then
        # we must hide the row name column, which is used to display the
        # `rowid`.
        if !get(kwargs, :show_row_number, true)
            row_names = nothing
            vlines = Int[]
        else
            row_names = [string(rowid)]
            vlines = Int[1]
        end

        show_row_number = false
    end

    # Print the table with the selected options.
    pretty_table(io, df;
                 alignment                   = alignment,
                 alignment_anchor_fallback   = :r,
                 alignment_anchor_regex      = alignment_anchor_regex,
                 compact_printing            = compact_printing,
                 crop                        = crop,
                 crop_num_lines_at_beginning = 2,
                 ellipsis_line_skip          = 3,
                 formatters                  = (_pretty_tables_general_formatter,),
                 header                      = (names_str, types_str),
                 header_alignment            = :l,
                 hlines                      = [:header],
                 highlighters                = (_PRETTY_TABLES_HIGHLIGHTER,),
                 maximum_columns_width       = maximum_columns_width,
                 newline_at_end              = false,
                 nosubheader                 = !eltypes,
                 row_name_alignment          = :r,
                 row_name_crayon             = Crayon(),
                 row_name_column_title       = string(rowlabel),
                 row_names                   = row_names,
                 row_number_alignment        = :r,
                 row_number_column_title     = string(rowlabel),
                 show_row_number             = show_row_number,
                 title                       = title,
                 vcrop_mode                  = :middle,
                 vlines                      = vlines,
                 kwargs...)

    return nothing
end

"""
    show([io::IO, ]df::AbstractDataFrame;
         allrows::Bool = !get(io, :limit, false),
         allcols::Bool = !get(io, :limit, false),
         allgroups::Bool = !get(io, :limit, false),
         rowlabel::Symbol = :Row,
         summary::Bool = true,
         eltypes::Bool = true,
         truncate::Int = 32,
         kwargs...)

Render a data frame to an I/O stream. The specific visual
representation chosen depends on the width of the display.

If `io` is omitted, the result is printed to `stdout`,
and `allrows`, `allcols` and `allgroups` default to `false`.

# Arguments
- `io::IO`: The I/O stream to which `df` will be printed.
- `df::AbstractDataFrame`: The data frame to print.
- `allrows::Bool `: Whether to print all rows, rather than
  a subset that fits the device height. By default this is the case only if
  `io` does not have the `IOContext` property `limit` set.
- `allcols::Bool`: Whether to print all columns, rather than
  a subset that fits the device width. By default this is the case only if
  `io` does not have the `IOContext` property `limit` set.
- `allgroups::Bool`: Whether to print all groups rather than
  the first and last, when `df` is a `GroupedDataFrame`.
  By default this is the case only if `io` does not have the `IOContext` property
  `limit` set.
- `rowlabel::Symbol = :Row`: The label to use for the column containing row numbers.
- `summary::Bool = true`: Whether to print a brief string summary of the data frame.
- `eltypes::Bool = true`: Whether to print the column types under column names.
- `truncate::Int = 32`: the maximal display width the output can use before
  being truncated (in the `textwidth` sense, excluding `…`).
  If `truncate` is 0 or less, no truncation is applied.
- `kwargs...`: Any keyword argument supported by the function `pretty_table` of
  PrettyTables.jl can be passed here to customize the output.

# Examples
```jldoctest
julia> using DataFrames

julia> df = DataFrame(A=1:3, B=["x", "y", "z"]);

julia> show(df, show_row_number=false)
3×2 DataFrame
 A      B
 Int64  String
───────────────
     1  x
     2  y
     3  z
```
"""
Base.show(io::IO,
          df::AbstractDataFrame;
          allrows::Bool = !get(io, :limit, false),
          allcols::Bool = !get(io, :limit, false),
          rowlabel::Symbol = :Row,
          summary::Bool = true,
          eltypes::Bool = true,
          truncate::Int = 32,
          kwargs...) =
    _show(io, df; allrows=allrows, allcols=allcols, rowlabel=rowlabel,
          summary=summary, eltypes=eltypes, truncate=truncate, kwargs...)

Base.show(df::AbstractDataFrame;
          allrows::Bool = !get(stdout, :limit, true),
          allcols::Bool = !get(stdout, :limit, true),
          rowlabel::Symbol = :Row,
          summary::Bool = true,
          eltypes::Bool = true,
          truncate::Int = 32,
          kwargs...) =
    show(stdout, df;
         allrows=allrows, allcols=allcols, rowlabel=rowlabel, summary=summary,
         eltypes=eltypes, truncate=truncate, kwargs...)


text:
#
# ======================================================================
#   Import some chemistry and materials science data tables from
#   ASE. Without this data, JuLIP can do very little!
# ======================================================================
#
using PyCall, JSON
@pyimport ase.data as ase_data
asedata = Dict(
     :symbols => ase_data.chemical_symbols,
      :masses => ase_data.atomic_masses,
   :refstates => ase_data.reference_states
   )

write(@__DIR__() * "/asedata.json", JSON.json(asedata, 0))

# ======================================================================

# NOTE:
# some other data that we could consider adding
# asedata.atomic_numbers
# ase_data.atomic_names
# ase_data.covalent_radii
# ase_data.ground_state_magnetic_moments
# ase_data.vdw_radii
#  can we get some more stuff like electron affinity somewhere?

# function rnn_old(species::Symbol)
#    X = positions(bulk(species) * 2)
#    return minimum( norm(X[n]-X[m]) for n = 1:length(X) for m = n+1:length(X) )
# end
#
# _rnn = fill(-1.0, length(_symbols))
# for n = 2:length(_symbols)
#    z = n-1
#    try
#       _rnn[n] = rnn_old(JuLIP.Chemistry.chemical_symbol(z))
#    catch
#    end
# end

# ======================================================================

# get access to the atomic numbers
ase_emt = pyimport("ase.calculators.emt")
ase_data = pyimport("ase.data")

function emt_default_parameters()
   # import everything we need from ASE
   Bohr = ase_emt.Bohr
   # get_atomic_numbers(id::AbstractString) = ase_data.atomic_numbers[id]
   ase_parameters = ase_emt.parameters
   beta = ase_emt.beta
   # convert ASE style Dict to a new Dict where parameters
   # are stored in a type instead of tuple.
   params = Dict{String, Dict{String,Float64}}()
   maxseq = maximum([par[2] for par in values(ase_parameters)]) * Bohr  # ✓
   rc = beta * maxseq * 0.5 * (sqrt(3) + sqrt(4))   # ✓
   rr = rc * 2 * sqrt(4) / (sqrt(3) + sqrt(4))   # ✓
   acut = log(9999.0) / (rr - rc)      # ✓

   for (key, p) in ase_parameters   # p is a tuple of parameters, key a species
      s0 = p[2] * Bohr
      eta2 = p[4] / Bohr
      kappa = p[5] / Bohr
      x = eta2 * beta * s0
      gamma1 = 0.0
      gamma2 = 0.0
      for (i, n) in enumerate([12; 6; 24])
         r = s0 * beta * sqrt(i)
         x = n / (12.0 * (1.0 + exp(acut * (r - rc))))
         gamma1 += x * exp(-eta2 * (r - beta * s0))
         gamma2 += x * exp(-kappa / beta * (r - beta * s0))
      end
      params[key] = Dict{String, Float64}(
                     "E0" => p[1], "s0" => s0, "V0" => p[3], "η2" => eta2,
                     "κ" => kappa, "λ" => p[6] / Bohr, "n0" => p[7] / Bohr^3,
                     "γ1" => gamma1, "γ2" => gamma2 )
      par = params[key]
      par["Cpair"] = -0.5 * par["V0"] / par["γ2"] / par["n0"]
   end

   return Dict("params" => params, "acut" => acut, "rc" => rc, "beta" => beta)
end

emt_data = emt_default_parameters()
write(@__DIR__() * "/emt.json", JSON.json(emt_data, 0))


text:
"""
    Newton_Raphson(f, fp, x; tol = 1e-14, MITR = 100)
---
Given function `f` and its derivative `fp`, this function finds a root for `f = 0` with
a starter `x`
"""
function Newton_Raphson(f, fp, x; tol = 1e-14, MITR = 100)
    xnew, xold = x, Inf
    fn, fo = f(xnew), Inf

    ctr = 1

    while (ctr < 100) && (abs(xnew - xold) > tol) && ( abs(fn - fo) > tol )
	x = xnew - f(xnew)/fp(xnew) # update step
	xnew, xold = x, xnew
        fn, fo = f(xnew), fn
	ctr = ctr + 1
    end

    if ctr == MITR
	error("Did not converge in $MITR steps")
    else
	xnew, ctr
    end
end

"""
    Newton_Raphson(f, x; tol = 1e-14, MITR = 100, inc = 1e-8)
---
Given function `f`, this function numerically finds a root for `f = 0` with a starter `x` using `secant` method.
"""
function Newton_Raphson(f, x; tol = 1e-14, MITR = 100, inc = 1e-8)
    fp(x1, x2) = (f(x1) - f(x2))/(x1 - x2)

    fn, fo = f(x1), Inf

    ctr = 1

    x2, x1 = x + inc, x
    while (ctr < 100) && (abs(x2 - x1) > tol) && ( abs(fn - fo) > tol )
        x2, x1 = x2 - f(x2)/fp(x1, x2), x2
        fn, fo = f(x2), fn
	ctr = ctr + 1
    end

    if ctr == 100
	error("Did not converge in 100 steps")
    else
	x2, ctr
    end
end


text:
# FIXME:: grooming

# --[ (a*b+c*d) -> muladd ]
# LJR

ex = :(a1*b1 + a2*b2 + a3*b3 + a4*b4)
# :(a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4)

co = [Symbol.(ex.args[i].args[2:3]) for i in 2:length(ex.args)]
#4-element Array{Array{Symbol,1},1}:
# Symbol[:a1, :b1]
# Symbol[:a2, :b2]
# Symbol[:a3, :b3]
# Symbol[:a4, :b4]

res = Expr(:call, :*, co[end]...)
# :(a4 * b4)

for i in (length(co)-1):-1:1
    res = :(muladd( $(co[i]...), $res ))
end

res
#:(muladd(a1, b1, muladd(a2, b2, muladd(a3, b3, a4 * b4))))

function _muladd_transform(ex)
    if length(ex.args) == 2
        return ex.args[2]
    else
    a, b = ex.args[2].args[2:end]
    rest = _muladd_transform(Expr(ex.head, :+, ex.args[3:end]...))
    return :($(Base.muladd)($a, $b, $rest))
    end
end

function _muladd_rewrite(ex)
    co = [Symbol.(ex.args[i].args[2:3]) for i in 2:length(ex.args)]
    res = Expr(:call, :*, co[end]...)
    for i in (length(co)-1):-1:1
        res = :(muladd( $(co[i]...), $res ))
    end
    res
end

# Symbolic diff
# LJR

struct SymbolFun{T} end
SymbolFun(S::Symbol) = SymbolFun{S}()

D(ex::Symbol, x::Symbol) = ex == x ? 1 :0
D(ex::Number, x::Symbol) = 0

function D(ex::Expr, x)
    @assert ex.head == :call
    D(SymbolFun(ex.args[1]), ex.args[2:end], x)
end

D(fun::SymbolFun{T}, args, x) where T = error("string(T) function is not supported!")

# Power Rule
function D(::SymbolFun{:^}, args, x)
    a, b = args[1], args[2]
    da, db = D(a, x), D(b, x)
    if da == 0 && db == 0    # if a^b is constant: a^b -> 0
        return 0
    elseif db == 0           # if b is constant: a^b -> b*da*a^(b-1)
        return :( $b* $da * ($a ^ ($b - 1)) )
    else                     # a(x)^b(x) -> http://www.wolframalpha.com/input/?source=frontpage-immediate-access&i=d%2Fdx+a(x)%5Eb(x)
                             # I don't know, just ask Wolfram Alpha!
        :( $a ^ ($b - 1) * ($da * $b + $a * log($a) * $db) )
    end
end

# Product Rule
# d/dx (f * g) = (d/dx f) * g + f * (d/dx g)
# d/dx (f * g * h) = (d/dx f) * g * h + f * (d/dx g) * h + ...
function D(::SymbolFun{:*}, args, x)
    N = length(args)
    outer_args = Vector{Any}(N)
    for i in 1:N
        inner_args = Vector{Any}(N)
        for j in 1:N
            i == j ? inner_args[j] = D(args[j], x) : inner_args[j] = args[j]
        end
        outer_args[i] = Expr(:call, :*, inner_args...)
    end
    Expr(:call, :+, outer_args...)
end

# The Quotient Rule
# d/dx (f / g) = ((d/dx f) * g - f * (d/dx g)) / g^2
function D(::SymbolFun{:/}, args, x)
    f, g = args[1], args[2]
    df, dg = D(f, x), D(g, x)
    if df == 0 && dg == 0
        return 0
    elseif df == 0
        return :( -$dg * $f / $g^2 )
    elseif dg == 0
        return :( $df / $g )
    else
        return :( ($df * $g - $f * $dg) / $g^2 )
    end
end

function D(::SymbolFun{:+}, args, x)
    terms = Any[:+]
    for y in args
        dx = D(y, x)
        if dx != 0
            push!(terms, dx)
        end
    end
    if (length(terms) == 1)
        return 0
    elseif (length(terms) == 2)
        return terms[2]
    else
        return Expr(:call, terms...)
    end
end

function D(::SymbolFun{:-}, args, x)
    terms = Any[:-]
    term1 = D(args[1], x)
    push!(terms, term1)
    for y in args[2:end]
        dx = D(y, x)
        if dx != 0
            push!(terms, dx)
        end
    end
    if term1 != 0 && length(terms) == 2 && length(args) >= 2
        return term1
    elseif (term1 == 0 && length(terms) == 2)
        return 0
    else
        return Expr(:call, terms...)
    end
end

D(fun::SymbolFun{T}, args, x) where T = error("$(string(T)) function is not supported!")

[1:57]
 ```julia> ex = :(11x^3 + 12 + 13x^5 + 14x^9)
:(11 * x ^ 3 + 12 + 13 * x ^ 5 + 14 * x ^ 9)

julia> Meta.show_sexpr(ex)
(:call, :+, (:call, :*, 11, (:call, :^, :x, 3)), 12, (:call, :*, 13, (:call, :^,
:x, 5)), (:call, :*, 14, (:call, :^, :x, 9)))```

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
# constructor for the vecchia factorization approxiamtion

"""
Vecchia Factorization approximation:
```
vecchia(Σ, blk_sizes) -> R⁻¹ M R⁻ᴴ
```
Pivoted Vecchia Factorization approximation:
```
vecchia(Σ, blk_sizes, perm) -> Pᵀ R⁻¹ M R⁻ᴴ P
```
where the argument where `Σ` is an `AbstractMatrix` or a `Function` taking indices `(i,j)` and
returning `Σ[i,j]`.
"""
function vecchia end

"""
Construct individual factors in the Pivoted Vecchia approximation `Pᵀ R⁻¹ M R⁻ᴴ P`.
```
R_M_P(Σ, blk_sizes, perm) -> R, M, P
```
where the argument where `Σ` is an `AbstractMatrix` or a `Function` taking indices `(i,j)` and
returning `Σ[i,j]`. The return values have the following types:
```
R <: Ridiagonal
M <: Midiagonal
P <: Piv
```
"""
function R_M_P end

function vecchia(Σ::Union{AbstractMatrix, Function}, blk_sizes::AbstractVector{<:Integer}, perm::AbstractVector{<:Integer})
	R, M, P = R_M_P(Σ, blk_sizes, perm)
    # println("Vecchia Factorization approximation: Pᵀ R⁻¹ M R⁻ᴴ P")
	return P' * inv(R) * M * inv(R)' * P
end

function vecchia(Σ::Union{AbstractMatrix, Function}, blk_sizes::AbstractVector{<:Integer})
	R, M, P = R_M_P(Σ, blk_sizes)
    # println("Vecchia Factorization approximation: R⁻¹ M R⁻ᴴ")
	return inv(R) * M * inv(R)'
end

function R_M_P(Σ::AbstractMatrix{T}, blk_sizes::AbstractVector{<:Integer}, perm::AbstractVector{<:Integer}=1:sum(blk_sizes)) where T
	LinearAlgebra.checksquare(Σ)
	@assert isperm(perm)

	blk_indices = blocks(PseudoBlockArray(perm, blk_sizes))
	N = length(blk_sizes)
	M = Vector{Matrix{T}}(undef, N)
	R = Vector{Matrix{T}}(undef, N-1)
	for ic in 1:N # loops over the column block index
		if ic == 1
			M[ic] = Σ[blk_indices[ic], blk_indices[ic]]
		else
			U 		= cholesky(Σ[blk_indices[ic-1], blk_indices[ic-1]]).U
			C 		= Σ[blk_indices[ic], blk_indices[ic-1]] / U
			R[ic-1] = - C / U'
			M[ic]   = Σ[blk_indices[ic], blk_indices[ic]] - C*C'
		end
	end

    return Ridiagonal(R), Midiagonal(M), Piv(perm)
end

# QUESTION: it seems redundant to have two of these R_M_P methods. Perhaps we should just have this one then use dispatch
# on vecchia(...) above to pass the correct argument
function R_M_P(Σfun::Function, blk_sizes::AbstractVector{<:Integer}, perm::AbstractVector{<:Integer}=1:sum(blk_sizes))
	@assert isperm(perm)

	blk_indices = blocks(PseudoBlockArray(perm, blk_sizes))
	N = length(blk_sizes)
	T = typeof(Σfun(1,1)) # This seems brittle. To do it right, check how `map` does it
	M = Vector{Matrix{T}}(undef, N)
	R = Vector{Matrix{T}}(undef, N-1)
	for ic in 1:N # loops over the column block index
		if ic == 1
			M[ic] = Σfun.(blk_indices[ic], blk_indices[ic]')
		else
			U 		= cholesky(Σfun.(blk_indices[ic-1], blk_indices[ic-1]')).U
			C 		= Σfun.(blk_indices[ic], blk_indices[ic-1]') / U
			R[ic-1] = - C / U'
			M[ic]   = Σfun.(blk_indices[ic], blk_indices[ic]') - C*C'
		end
	end

    return Ridiagonal(R), Midiagonal(M), Piv(perm)
end


text:
"""
    AVSfldIO
module for AVS .fld file IO
"""
module AVSfldIO

using FileIO: File, @format_str

include("fld-read.jl")
include("fld-write.jl")

# the two key methods that will be used by FileIO:

"""
    data = load(ff::File{format"AVSfld"} ; kwargs...)
AVSfld load method
"""
load(ff::File{format"AVSfld"} ; kwargs...) =
   fld_read(ff.filename ; kwargs...)

"""
    save(ff::File{format"AVSfld"}, data::AbstractArray{<:Real} ; kwargs...)
AVSfld save method
"""
save(ff::File{format"AVSfld"}, data::AbstractArray{<:Real} ; kwargs...) =
   fld_write(ff.filename, data ; kwargs...)

save(ff::File{format"AVSfld"}, data ; kwargs...) =
    throw(ArgumentError("data must be an array of reals for AVSfld"))

end # module

max_stars_repo_name: JeffFessler/AVSfldIO

max_issues_repo_name: JeffFessler/AVSfldIO


text:
function genInterpolants(data)
    # Get rows of data
    rows    = size(data, 1)

    # Allocate matrix for satelite data
    satData = zeros(Int(floor(rows / 30)), 4)

    # Generate interpolants for each GPS Satelite
    local df
    @showprogress "Generating spline interpolants..." for i in 1:32
        # Cound epochs of data for sat i and put in satData
        count = 0
        idxs  = [1,3,4,5]
        for j in 1:rows
            if data[j, 2] == i
                count += 1
                @views satData[count, :] .= data[j, idxs]
            end
        end

        if i == 1
            x  = CubicSpline(satData[1:count,2], satData[1:count,1])
            y  = CubicSpline(satData[1:count,3], satData[1:count,1])
            z  = CubicSpline(satData[1:count,4], satData[1:count,1])

            df = DataFrame(SAT = i, X = [x], Y = [y], Z = [z])
        else
            x  = CubicSpline(satData[1:count,2], satData[1:count,1])
            y  = CubicSpline(satData[1:count,3], satData[1:count,1])
            z  = CubicSpline(satData[1:count,4], satData[1:count,1])
            push!(df, Dict(
                "SAT" => i, "X" => x, "Y" => y, "Z" => z
            ))
        end
    end
    return df
end

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
@testset "Weighted" begin
    @testset "Known V, Unknown W" begin
        Y = [[0.597, 1.299, 3.27, -1.27, 7.59, -6.89],
             [-1.0, -4.67, 1.88, -1.39, 5.09, 2.4],
             [-5.63, -11.4, -8.56, -4.54, 22.5, 15.2],
             [-19.1, -13.6, -12.0, -17.7, 21.0, 28.5],
             [-25.5, -17.8, -21.8, -20.1, 44.4, 42.6],
             [-34.1, -23.5, -25.6, -22.8, 58.4, 53.8]]
        Y = to_matrix(Y)
        η = [0.45, 0.45, 0.1]

        F = [1. 0. 0.; 0. 1. 0.]
        G = [1. 0. 1.; 0. 1. 1.; 0. 0. 1.]
        V = Symmetric(diagm([12., 12.]))

        # Filter

        a, R, m, C = kfilter(Y, F, G, V, η, 0.8)

        @test a[end,:] ≈ [-18.794803420834505, -18.62873963779488, -4.093136957816489]
        @test R[end][3,3] ≈ 1.0514508656899235
        @test C[end][3,3] ≈ 0.5565127730469597

        # Smoother

        s, S = ksmoother(G, a, R, m, C)

        @test s[1,:] ≈ [1.6363856970098998, 0.6806573714011395, -2.4041282160885644]
        @test diag(S[1]) ≈ [5.401963801663216, 5.401963801663269, 6.64771001788035]

        # Fitted

        f, Q = fitted(F, V, s, S)

        @test f[1,:] ≈ [1.6363856970098998, 0.6806573714011395]
        @test diag(Q[end]) ≈ [17.27516049014339, 17.27516049014339]

        # Forecast

        fh, Qh = forecast(F, G, V, 0.8, m[end,:], C[end], 10)

        @test fh[end,:] ≈ [-59.3151969404975, -57.705730042143614]
        @test diag(Qh[end]) ≈ [479.334143445396, 479.334143445396]
    end

    @testset "Unknown W, Unknown V" begin
        Y = [[0.597, 1.299, 3.27, -1.27, 7.59, -6.89],
             [-1.0, -4.67, 1.88, -1.39, 5.09, 2.4],
             [-5.63, -11.4, -8.56, -4.54, 22.5, 15.2],
             [-19.1, -13.6, -12.0, -17.7, 21.0, 28.5],
             [-25.5, -17.8, -21.8, -20.1, 44.4, 42.6],
             [-34.1, -23.5, -25.6, -22.8, 58.4, 53.8]]
        Y = to_matrix(Y)
        η = [0.7, 0.2, 0.1]

        F = [1. 0. 0.; 0. 1. 0.]
        G = [1. 0. 1.; 0. 1. 1.; 0. 0. 1.]

        # Estimate
        θ, Σ, _ = estimate(Y, F, G, η, 0.8)

        @test θ[end,:] ≈ [-19.05578625907453, -16.79719558628691, -3.266323626397268]
        @test Σ ≈ [245.88141984197944 0; 0 186.51152750423392]

        # Estimate with exclusion
        θ, Σ, _ = estimate(Y, F, G, η, 0.8, exclude_under=0.15)

        @test θ[end,:] ≈ [-27.303258240406077, -23.792801359688326, -4.886906081116004]
        @test Σ ≈ [18.275634284242617 0; 0 2.814105974385866]

        # Evolutional Covariances
        Σₑ = evolutional_covariances(Y, F, G, Σ, η, 0.8)

        @test Σₑ[end][3,3] ≈ 0.05964878602083959

        # Fitted
        f, Q = fitted(F, Σ, θ, Σₑ)

        @test f[end,:] ≈ [-27.303258240406077, -23.792801359688326]
        @test Q[end][1,2] ≈ 0.508856381882239

        # Forecast
        fh, Qh = forecast(F, G, Σ, 0.8, θ[end,:], Σₑ[end], 10)

        @test fh[end,:] ≈ [-76.17231905156612, -72.66186217084837]
        @test Qh[end][1,2] ≈ 57.17048592729376
    end
end


text:
module DisplacementWrapper

using JavaCall

const BodiesElements = @jimport org.orekit.data.BodiesElements
const Class = @jimport java.lang.Class
const DataProvidersManager = @jimport org.orekit.data.DataProvidersManager
const Frame = @jimport org.orekit.frames.Frame
const GeodeticPoint = @jimport org.orekit.bodies.GeodeticPoint
const IERSConventions = @jimport org.orekit.utils.IERSConventions
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const Object = @jimport java.lang.Object
const OceanLoading = @jimport org.orekit.models.earth.displacement.OceanLoading
const OceanLoadingCoefficients = @jimport org.orekit.models.earth.displacement.OceanLoadingCoefficients
const OceanLoadingCoefficientsBLQFactory = @jimport org.orekit.models.earth.displacement.OceanLoadingCoefficientsBLQFactory
const OneAxisEllipsoid = @jimport org.orekit.bodies.OneAxisEllipsoid
const PVCoordinatesProvider = @jimport org.orekit.utils.PVCoordinatesProvider
const StationDisplacement = @jimport org.orekit.models.earth.displacement.StationDisplacement
const TidalDisplacement = @jimport org.orekit.models.earth.displacement.TidalDisplacement
const Tide = @jimport org.orekit.models.earth.displacement.Tide
const Vector3D = @jimport org.hipparchus.geometry.euclidean.threed.Vector3D

export OceanLoading
export OceanLoadingCoefficients
export OceanLoadingCoefficientsBLQFactory
export StationDisplacement
export TidalDisplacement
export Tide
export displacement
export equals
export get_class
export get_coefficients
export get_delaunay_multipliers
export get_doodson_multipliers
export get_doodson_number
export get_nb_species
export get_nb_tides
export get_phase
export get_rate
export get_site_location
export get_site_name
export get_sites
export get_south_amplitude
export get_south_phase
export get_tau_multiplier
export get_tide
export get_west_amplitude
export get_west_phase
export get_zenith_amplitude
export get_zenith_phase
export hash_code
export notify
export notify_all
export to_string
export wait

include("ocean_loading.jl")
include("ocean_loading_coefficients.jl")
include("ocean_loading_coefficients_blq_factory.jl")
include("station_displacement.jl")
include("tidal_displacement.jl")
include("tide.jl")

end


text:
# =================================================================
# Laub-Loomis model
# See https://easychair.org/publications/paper/gjfh
# =================================================================

using Reachability, MathematicalSystems, LazySets, TaylorIntegration
using Reachability: solve

# Equations of motion
# We write the function such that the operations are either unary or binary:
@taylorize function laubloomis!(t, x, dx)
    dx[1] = 1.4*x[3] - 0.9*x[1]
    dx[2] = 2.5*x[5] - 1.5*x[2]
    dx[3] = 0.6*x[7] - 0.8*(x[2]*x[3])
    dx[4] = 2 - 1.3*(x[3]*x[4])
    dx[5] = 0.7*x[1] - (x[4]*x[5])
    dx[6] = 0.3*x[1] - 3.1*x[6]
    dx[7] = 1.8*x[6] - 1.6*(x[2]*x[7])
    return dx
end

function laubloomis(; T=20.0, W=0.01, plot_vars=[0, 4],
                      property=(t,x)->x[4] < 4.5,
                      project_reachset=true)

    # equations, x' = f(x(t))
    𝐹 = BlackBoxContinuousSystem(laubloomis!, 7)

    X0c = [1.2, 1.05, 1.5, 2.4, 1.0, 0.1, 0.45]
    X0 = Hyperrectangle(X0c, fill(W, 7))

    # instantiate the IVP
    𝑃 = InitialValueProblem(𝐹, X0)

    # general options
    𝑂 = Options(:T=>T, :plot_vars=>plot_vars, :property=>property,
                :project_reachset=>project_reachset, :mode=>"check")

    return (𝑃, 𝑂)
end

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
using CSO
using Test

using DataFrames, CSV, Dates, StatsBase, Plots, MLBase, GLM, Combinatorics, DecisionTree

@testset "CSO.jl" begin
    include("dataAnalysis_test.jl")
    include("dataProcessing_test.jl")
end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
using OptimMPB
using Base.Test

type EX1 <: MathProgBase.AbstractNLPEvaluator end

MathProgBase.features_available(d::EX1) = [:Grad]
MathProgBase.eval_f(d::EX1, x) = (309.0 - 5.0 * x[1])^2 + (17.0 - x[2])^2

function MathProgBase.eval_grad_f(d::EX1, gr, x)
  gr[1] = -10.0 * (309.0 - 5.0 * x[1])
  gr[2] = -2.0 * (17.0 - x[2])
end

function grad_f(gr, x)
  gr[1] = -10.0 * (309.0 - 5.0 * x[1])
  gr[2] = -2.0 * (17.0 - x[2])
end

function MathProgBase.jac_structure(d::EX1)
    Int[],Int[]
end

function MathProgBase.eval_jac_g(d::EX1, J, x)
    nothing
end

for method in (AcceleratedGradientDescent, BFGS, ConjugateGradient,
               GradientDescent, LBFGS, MomentumGradientDescent, NelderMead)
  s = OptimSolver(method)
  m = MathProgBase.NonlinearModel(s)
  MathProgBase.loadproblem!(m, 2, [-100.,-100.], [200.,200.], :Min, EX1())
  MathProgBase.setwarmstart!(m, [0.1,-0.1])
  MathProgBase.optimize!(m)

  @test Optim.minimizer(m.inner.out) == MathProgBase.getsolution(m)
  @test Optim.minimum(m.inner.out) == MathProgBase.getobjval(m)
  @test MathProgBase.getvarLB(m) == [-100., -100]
  @test MathProgBase.getvarUB(m) == [200., 200]
  @test MathProgBase.status(m) == Optim.converged(m.inner.out)
  @test norm(Optim.minimizer(m.inner.out) .- [309.0 / 5.0, 17.0]) < 0.01
  if method != NelderMead
      @test Optim.g_converged(m.inner.out)
  end
end

for method in (AcceleratedGradientDescent, BFGS, ConjugateGradient,
               GradientDescent, LBFGS, MomentumGradientDescent, NelderMead)
  s = OptimSolver(method)
  m = MathProgBase.NonlinearModel(s)
  MathProgBase.loadproblem!(m, 2, [-Inf,-Inf], [+Inf,+Inf], :Min, EX1())
  MathProgBase.setwarmstart!(m, [0.1,-0.1])
  MathProgBase.optimize!(m)

  @test Optim.minimizer(m.inner.out) == MathProgBase.getsolution(m)
  @test Optim.minimum(m.inner.out) == MathProgBase.getobjval(m)
  @test MathProgBase.getvarLB(m) == [-Inf, -Inf]
  @test MathProgBase.getvarUB(m) == [+Inf, +Inf]
  @test MathProgBase.status(m) == Optim.converged(m.inner.out)
  @test norm(Optim.minimizer(m.inner.out) .- [309.0 / 5.0, 17.0]) < 0.01
  if method != NelderMead
      @test Optim.g_converged(m.inner.out)
  end
end

for method in (BFGS, LBFGS, GradientDescent)
  s = OptimSolver(method, show_trace = true, store_trace = true, extended_trace = true)
  m = MathProgBase.NonlinearModel(s)
  MathProgBase.loadproblem!(m, 2, [-100.,-100.], [200.,200.], :Min, EX1())
  MathProgBase.setwarmstart!(m, [0.1,-0.1])
  MathProgBase.optimize!(m)

  @test Optim.minimizer(m.inner.out) == MathProgBase.getsolution(m)
  @test Optim.minimum(m.inner.out) == MathProgBase.getobjval(m)
  @test MathProgBase.getvarLB(m) == [-100., -100]
  @test MathProgBase.getvarUB(m) == [200., 200]
  @test MathProgBase.status(m) == Optim.converged(m.inner.out)
  @test norm(Optim.minimizer(m.inner.out) - [309.0 / 5.0, 17.0]) < 0.01
  @test Optim.g_converged(m.inner.out)
end

type EX2 <: MathProgBase.AbstractNLPEvaluator end

MathProgBase.features_available(d::EX2) = []
MathProgBase.eval_f(d::EX2, x) = (309.0 - 5.0 * x[1])^2 + (17.0 - x[2])^2

# function MathProgBase.eval_grad_f(d::EX2, gr, x)
#   gr[1] = -10.0 * (309.0 - 5.0 * x[1])
#   gr[2] = -2.0 * (17.0 - x[2])
# end

for method in (BFGS, LBFGS, GradientDescent)
  s = OptimSolver(method)
  m = MathProgBase.NonlinearModel(s)
  MathProgBase.loadproblem!(m, 2, [-100.,-100.], [200.,200.], :Min, EX2())
  MathProgBase.setwarmstart!(m, [0.1,-0.1])
  MathProgBase.optimize!(m)

  @test Optim.minimizer(m.inner.out) == MathProgBase.getsolution(m)
  @test Optim.minimum(m.inner.out) == MathProgBase.getobjval(m)
  @test MathProgBase.getvarLB(m) == [-100., -100]
  @test MathProgBase.getvarUB(m) == [200., 200]
  @test MathProgBase.status(m) == Optim.converged(m.inner.out)
  @test norm(Optim.minimizer(m.inner.out) .- [309.0 / 5.0, 17.0]) < 0.01
  @test Optim.g_converged(m.inner.out)
end

## CSMINWEL
# using CsminWel
# s = OptimSolver(Csminwel, show_trace = true, store_trace = true, extended_trace = true)
# m = MathProgBase.NonlinearModel(s)
# MathProgBase.loadproblem!(m, 2, [-Inf,-Inf], [+Inf,+Inf], :Min, EX1())
# MathProgBase.setwarmstart!(m, [0.1,-0.1])
# MathProgBase.optimize!(m)
# @test Optim.minimizer(m.inner.out) == MathProgBase.getsolution(m)
# @test Optim.minimum(m.inner.out) == MathProgBase.getobjval(m)
# @test MathProgBase.getvarLB(m) == [-Inf, -Inf]
# @test MathProgBase.getvarUB(m) == [+Inf, +Inf]
# @test MathProgBase.status(m) == Optim.converged(m.inner.out)
# @test norm(Optim.minimizer(m.inner.out) - [309.0 / 5.0, 17.0]) < 0.01
# @test Optim.g_converged(m.inner.out)
# @test MathProgBase.getsense(m) == 1
# @test MathProgBase.numvar(m) == 2
# @test MathProgBase.numconstr(m) == 0


text:
#=
    EodData internal utility functions
=#

# ==================
# Internal constants
const WS = "http://ws.eoddata.com/data.asmx"
const HOST_ADDRESS = "ws.eoddata.com"
const CONTENT_TYPE = "application/x-www-form-urlencoded"
const REQUEST_TIMEOUT = 600.0

# =========
# Functions

# get_response()
# --------------
# Returns an xml tree.
# INPUT: Web Service Call, Web Service Call Parameters
# OUTPUT: Xml tree of the type ::ETree
# http://ws.eoddata.com/data.asmx?op=CountryList
function get_response(call::String, params::Dict{ASCIIString, ASCIIString})
    # Get response
    resp = HTTPC.post("$WS$call",
                      params,
                      RequestOptions(headers=[("Host",HOST_ADDRESS)],
                                    content_type=CONTENT_TYPE,
                                    request_timeout=REQUEST_TIMEOUT))
    validate_http_response(resp) && return xp_parse(bytestring(resp.body))
end

# set_countries()
# --------------
# Returns a collection of countries
# INPUT: XML tree
# OUTPUT:  Dict() of countries of type ::OrderedDict{String, Country}
function set_countries(xml_tree::ETree)
    countries = OrderedDict{String, Country}()
    for c_xml in find(xml_tree, "/RESPONSE/COUNTRIES/CountryBase")
        countries[c.code] = c = Country(c_xml)
    end
    return countries
end

# set_data_formats()
# --------------
# Returns a collection of data formats
# INPUT: XML tree
# OUTPUT:  Dict() of data formats of type ::OrderedDict{String, Exchange}
function set_data_formats(xml_tree::ETree)
    formats = OrderedDict{String, DataFormat}()
    for df_xml in find(xml_tree, "/RESPONSE/DATAFORMATS/DATAFORMAT")
        formats[df.code] = df = DataFormat(df_xml)
    end
    return formats
end

# set_exchanges()
# --------------
# Returns a collection of exchanges
# INPUT: XML tree
# OUTPUT:  Dict() of exchanges of type ::OrderedDict{String, Exchange}
function set_exchanges(xml_tree::ETree)
    exchanges = OrderedDict{String, Exchange}()
    for ex_xml in find(xml_tree, "/RESPONSE/EXCHANGES/EXCHANGE")
        exchanges[ex.code] = ex = Exchange(ex_xml)
    end
    return exchanges
end

# set_fundamentals()
# --------------
# Returns a collection of fundamentals
# INPUT: XML tree
# OUTPUT:  Dict() of fundamentals of type ::OrderedDict{String, Fundamental}
function set_fundamentals(xml_tree::ETree)
    fundamentals = OrderedDict{String, Fundamental}()
    for fl_xml in find(xml_tree, "/RESPONSE/FUNDAMENTALS/FUNDAMENTAL")
        fundamentals[fl.ticker_code] = fl = Fundamental(fl_xml)
    end
    return fundamentals
end

# set_quotes()
# --------------
# Returns a collection of quotes
# INPUT: XML tree
# OUTPUT:  Dict() of quotes of type ::OrderedDict{String, Quote}
function set_quotes(xml_tree::ETree)
    quotes = OrderedDict{String, Quote}()
    for qt_xml in find(xml_tree, "/RESPONSE/QUOTES/QUOTE")
        quotes[qt.ticker_code * "_" * string(qt.date_time)] = qt = Quote(qt_xml)
    end
    return quotes
end

# set_quotes_2()
# --------------
# Returns a collection of quotes
# INPUT: XML tree
# OUTPUT:  Dict() of quotes of type ::OrderedDict{String, Quote_2}
function set_quotes_2(xml_tree::ETree)
    quotes = OrderedDict{String, Quote_2}()
    for qt_2_xml in find(xml_tree, "/RESPONSE/QUOTES2/QUOTE2")
        quotes[qt_2.ticker_code * "_" * string(qt_2.date_time)] = qt_2 = Quote_2(qt_2_xml)
    end
    return quotes
end

# set_splits()
# --------------
# Returns a collection of splits
# INPUT: XML tree
# OUTPUT:  Dict() of splits of type ::OrderedDict{String, Split}
function set_splits(xml_tree::ETree)
    splits = OrderedDict{String, Split}()
    for sp_xml in find(xml_tree, "/RESPONSE/SPLITS/SPLIT")
        splits[sp.ticker_code * "_" * string(sp.date_time)] = sp = Split(sp_xml)
    end
    return splits
end

# set_technicals()
# --------------
# Returns a collection of technicals
# INPUT: XML tree
# OUTPUT: Dict() of technicals of type ::OrderedDict{String, Technical}
function set_technicals(xml_tree::ETree)
    # Shred xml_tree into a OrderedDict{String, Technical}
    technicals = OrderedDict{String, Technical}()
    for tl_xml in find(xml_tree, "/RESPONSE/TECHNICALS/TECHNICAL")
        technicals[tl.ticker_code] = tl = Technical(tl_xml)
    end
    return  technicals
end

# set_tickers()
# --------------
# Returns a collection of tickers
# INPUT: XML tree
# OUTPUT: Dict() of tickers of type ::OrderedDict{String, Ticker}
function set_tickers(xml_tree::ETree)
    tickers = OrderedDict{String, Ticker}()
    for tk_xml in find(xml_tree, "/RESPONSE/SYMBOLS/SYMBOL")
        tickers[tk.code] = tk = Ticker(tk_xml)
    end
    return tickers
end

# set_tickers_2()
# --------------
# Returns a collection of tickers
# INPUT: XML tree
# OUTPUT: Dict() of tickers of type ::OrderedDict{String, Ticker_2}
function set_tickers_2(xml_tree::ETree)
    tickers = OrderedDict{String, Ticker_2}()
    for tk_2_xml in find(xml_tree, "/RESPONSE/SYMBOLS2/SYMBOL2")
        tickers[tk_2.code] = tk_2 = Ticker_2(tk_2_xml)
    end
    return tickers
end

# set_ticker_changes()
# --------------
# Returns a collection of ticker changes
# INPUT: XML tree
# OUTPUT: Dict() of ticker changes of type ::OrderedDict{String, TickerChange}
function set_ticker_changes(xml_tree::ETree)
    ticker_changes = OrderedDict{String, TickerChange}()
    for tc_xml in find(xml_tree, "/RESPONSE/SYMBOLCHANGES/SYMBOLCHANGE")
        ticker_changes[tc.old_ticker_code * "_" * string(tc.date_time)] = tc = TickerChange(tc_xml)
    end
    return ticker_changes
end

# validate_http_response()
# --------------
# Returns true if an http code = 200 is received, error otherwise
# INPUT: Http response
# OUTPUT: Bool, or error
function validate_http_response(resp::Response)
    if resp.http_code != 200
        error("validate_http_response() failed with the following information: $(resp.headers)")
    else
        return true
    end
end

# validate_xml()
# --------------
# Returns true if the web service call returned a message of success, error otherwise
# INPUT: XML tree
# OUTPUT: Bool, or error
function validate_xml(xml_tree::ETree)
    message = lowercase(strip(find(xml_tree, "/RESPONSE[1]{Message}")))
    if message != "success"
        error("validate_xml() failed with message returned of: $message")
    else
        return true
    end
end

# validate_xml_login()
# --------------
# Returns true if the web service call to login returned a message of success, error otherwise
# INPUT: XML tree
# OUTPUT: Bool, or error
function validate_xml_login(xml_tree::ETree)
    message = lowercase(strip(find(xml_tree, "/LOGINRESPONSE[1]{Message}")))
    if message != "login successful"
        error("validate_xml_login() failed with message returned of: $message")
    else
        return true
    end
end


text: infinity() = Inf

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _mag2geo{T<:AbstractFloat}(lat::T, long::T, pole_lat::T, pole_long::T)
    r       = 1.0 # Distance from planet center.  Value unimportant -- just need
                  # a length for conversion to rectangular coordinates
    lat  = deg2rad(lat)
    long = deg2rad(long)
    alt  = lat*0.0 + r

    # convert to rectangular coordinates
    #   x-axis: defined by the vector going from Earth's center towards
    #        the intersection of the equator and Greenwich's meridian.
    #   z-axis: axis of the geographic poles
    #   y-axis: defined by y=z^x
    x = alt*cos(lat)*cos(long)
    y = alt*cos(lat)*sin(long)
    z = alt*sin(lat)

    # First rotation: in the plane of the current meridian from magnetic pole to
    # geographic pole.
    togeolat = Array(T, 3, 3)
    togeolat[:,1] = [cos(pi/2 - pole_lat), 0.0, -sin(pi/2 - pole_lat)]
    togeolat[:,2] = [                 0.0, 1.0,                   0.0]
    togeolat[:,3] = [sin(pi/2 - pole_lat), 0.0,  cos(pi/2 - pole_lat)]
    out = togeolat * [x, y, z]

    # Second rotation matrix: rotation around plane of the equator, from the
    # meridian containing the magnetic poles to the Greenwich meridian.
    maglong2geolong = Array(T, 3, 3)
    maglong2geolong[:,1] = [ cos(pole_long), sin(pole_long), 0.0]
    maglong2geolong[:,2] = [-sin(pole_long), cos(pole_long), 0.0]
    maglong2geolong[:,3] = [            0.0,            0.0, 1.0]
    out = maglong2geolong * out

    geolat  = rad2deg(atan2(out[3], hypot(out[1], out[2])))
    geolong = rad2deg(atan2(out[2], out[1]))
    # I don't care about that one...just put it there for completeness' sake
    # magalt  = vecnorm(out) - r
    return geolat, geolong
end

"""
    mag2geo(latitude, longitude[, year]) -> geographic_latitude, geographic_longitude

### Purpose ###

Convert from geomagnetic to geographic coordinates.

### Explanation ###

Converts from geomagnetic (latitude, longitude) to geographic (latitude,
longitude).  Altitude is not involved in this function.

### Arguments ###

* `latitude`: geomagnetic latitude (North), in degrees.
* `longitude`: geomagnetic longitude (East), in degrees.
* `year` (optional numerical argument): the year in which to perform conversion.
  If omitted, defaults to current year.

The coordinates can be passed as arrays of the same length.

### Output ###

The 2-tuple of geographic (latitude, longitude) coordinates, in degrees.

If geomagnetic coordinates are given as arrays, a 2-tuple of arrays of the same
length is returned.

### Example ###

Find position of North Magnetic Pole in 2016

``` julia
mag2geo(90, 0, 2016)
# => (86.395,-166.29000000000002)
```

### Notes ###

This function uses list of North Magnetic Pole positions provided by World
Magnetic Model (https://www.ngdc.noaa.gov/geomag/data/poles/NP.xy).

`geo2mag` converts geographic coordinates to geomagnetic coordinates.

Code of this function is based on IDL Astronomy User's Library.
"""
mag2geo(lat::Real, long::Real, year::Real=Dates.year(Dates.now())) =
    _mag2geo(promote(float(lat), float(long),
                     deg2rad(POLELATLONG[year][1]::Float64),
                     deg2rad(POLELATLONG[year][2]::Float64))...)

function mag2geo{LA<:Real, LO<:Real}(lat::AbstractArray{LA},
                                     long::AbstractArray{LO},
                                     year::Real=Dates.year(Dates.now()))
    @assert length(lat) == length(long)
    typela   = typeof(float(one(LA)))
    geolat   = similar(lat, typela)
    geolong  = similar(lat, typela)
    polelat  = deg2rad(POLELATLONG[year][1])
    polelong = deg2rad(POLELATLONG[year][2])
    for i in eachindex(lat)
        geolat[i], geolong[i] =
            _mag2geo(promote(float(lat[i]), float(long[i]),
                             polelat, polelong)...)
    end
    return geolat, geolong
end


text:
"""
    module TestLinearEquations - Test LinearEquations(..) in ModiaBase/src/EquationAndStateInfo.jl
"""
module TestLinearEquations

using Test
using Modia
using Modia.TimerOutputs

println("... Test TestLinearEquations")

# Solve A*x = b as residual equation
const A = [1.0 2.0 3.0 4.0 5.0;
           0.0 2.0 0.0 0.0 0.0;
           0.0 0.0 3.0 0.0 0.0;
           0.0 0.0 0.0 4.0 0.0;
           0.0 0.0 0.0 0.0 5.0]
const b = [-1.0, -2.0, -3.0, -4.0, -5.0]
const isInitial = true       # the first call must be performed with isInitial=true
const time  = 0.0            # only for warning/error messages
const timer = TimerOutput()  # timer to measure the time to solve the linear equation system
const leq   = Modia.LinearEquations(["x2", "x1"], [:x1], [1,5], 1, false)

computeResidual1(x1)    = A*x1-b
computeResidual2(x1,x2) = x1[1] + 2*x1[2] + 3*x1[3] + 4*x2

# Solution
const x1_sol = A\b
const x2_sol = -(x1_sol[1] + 2*x1_sol[2] + 3*x1_sol[3])/4.0

# Test function for ODE mode
function solveLinearEquation(leq,isInitial,time,timer)
    local x1::Vector{Float64}, x2
    leq.mode = -3
    while Modia.LinearEquationsIteration!(leq, isInitial, time, timer)
        x1 = leq.x_vec[1]
        x2 = leq.x[1]
        append!(leq.residuals, computeResidual1(x1))
        append!(leq.residuals, computeResidual2(x1,x2))
    end
    return (x1,x2)
end

# Test function for DAE mode
function solveLinearEquation(leq,isInitial,solve_leq,isStoreResult,time,timer)
    local x1, x2
    leq.mode = -3
    while Modia.LinearEquationsIteration!(leq, isInitial, solve_leq, isStoreResult, time, timer)
        x1 = leq.x_vec[1]
        x2 = leq.x[1]
        append!(leq.residuals, computeResidual1(x1))
        append!(leq.residuals, computeResidual2(x1,x2))
    end
    return (x1,x2)
end

@testset "Test TestLinearEquations" begin
    # Test ODE mode (x is always explicitly computed by solving the linear equation system)
    leq.odeMode = true
    (x1,x2) = solveLinearEquation(leq,isInitial,time,timer)

    d1   = x1_sol-x1
    err1 = sqrt(d1'*d1) + abs(x2_sol - x2)
    @test isapprox(err1, 0.0, atol=1e-13)

    # Test DAE mode
    leq.odeMode = false

    # x is computed explicitly at initialization or at an event
    (x1,x2) = solveLinearEquation(leq,false,true,false,time,timer)
    d2   = x1_sol-x1
    err2 = sqrt(d2'*d2) + abs(x2_sol - x2)
    @test isapprox(err2, 0.0, atol=1e-13)

    # x is computed by DAE solver
    leq.x = [6.0, 1.0, 2.0, 3.0, 4.0, 5.0]  # Provide x from DAE solver
    (x1,x2) = solveLinearEquation(leq,false,false,false,time,timer)
    residuals = leq.residuals      # Store residuals in DAE solver
    d3 = leq.x - vcat(x2,x1)
    d4 = vcat(A*leq.x[2:6] - b, leq.x[2] + 2*leq.x[3] + 3*leq.x[4] + 4*leq.x[1]) - residuals
    err3 = sqrt(d3'*d3)
    err4 = sqrt(d4'*d4)
    @test isapprox(err3, 0.0, atol=1e-15)
    @test isapprox(err4, 0.0, atol=1e-15)

    # x is computed by DAE solver at a communication point (leq.residuals is not provided)
    leq.x = [6.0, 1.0, 2.0, 3.0, 4.0, 5.0]  # Provide x from DAE solver
    (x1,x2) = solveLinearEquation(leq,false,false,true,time,timer)
    d5 = leq.x - vcat(x2,x1)
    err5 = sqrt(d5'*d5)
    @test isapprox(err5, 0.0, atol=1e-15)
end

end


text:
module DocumentFormat

import FilePathsBase
import JuliaFormatter

export format, isformatted

function default_juliaformatter_config()
    return (
        indent=4,
        annotate_untyped_fields_with_any=false,
        join_lines_based_on_source=true,
        trailing_comma=nothing,
        margin=10_000,
        always_for_in=nothing,
        whitespace_in_kwargs=false
    )
end

function format(path::FilePathsBase.AbstractPath)
    if isdir(path)
        JuliaFormatter.format(string(path); default_juliaformatter_config()...)
    else
        error("Invalid path.")
    end

    return nothing
end

function isformatted(path::FilePathsBase.AbstractPath)
    if isdir(path)
        return JuliaFormatter.format(string(path); overwrite=false, default_juliaformatter_config()...)
    else
        error("Invalid path.")
    end
end

end


text:
using MPGenbb
using Test

@testset "MPGenbb.jl" begin
    # Write your tests here.
end