text:
module thick_pipe_examples
using FinEtools
using PyCall

struct MyIData
    c::FInt
    r::FFltVec
    s::FFltVec
end

function thick_pipe_axi()
    println("Thick pipe with internal pressure: axially symmetric model")
    #=
    This is a simple modification of the full three-dimensional simulation of
    the tutorial pub_thick_pipe that implements the axially-symmetric model
    reduction procedure.

    An infinitely long thick walled cylindrical pipe
    with inner boundary radius of 3 mm and outer boundary radius of 9 mm is
    subjected to an internal pressure of 1.0 MPa. A wedge   with thickness of
    2 mm and a 90-degree angle sector is considered for the finite element
    analysis. The material properties are taken as  isotropic linear elastic
    with $E=1000$ MPa and $\nu=0.4999$ to represent nearly incompressible
    behavior. This problem has been proposed to by MacNeal and Harder as a
    test of an element's ability to represent the  response of a nearly
    incompressible material. The plane-strain condition is assumed in the
    axial direction of the pipe which together with the radial symmetry
    confines the material in all but the radial direction and therefore
    amplifies the numerical difficulties associated with the confinement of
    the nearly incompressible material.

    There is an analytical solution to this problem. Timoshenko and Goodier
    presented the original solution of Lame in their textbook. We are going
    to compare with  both the stress distribution (radial and hoop stresses)
    and the displacement of the inner  cylindrical surface.

    References:
    - Macneal RH, Harder RL (1985) A proposed standard set of problems to test
    finite element accuracy. Finite Elements in Analysis and Design 1: 3-20.
    - Timoshenko S. and Goodier J. N., Theory of Elasticity, McGraw-Hill, 2nd ed., 1951.

    =#

    # Internal radius of the pipe.
    a = 3*phun("MM");
    ##
    # External radius of the pipe.
    b = 9*phun("MM");
    ##
    # Thickness of the slice.
    t = 2*phun("MM");

    ##
    # Geometrical tolerance.
    tolerance   = a/10000.;
    ##
    # Young's modulus and Poisson's ratio.
    E = 1000*phun("MEGA*PA");
    nu = 0.499;
    ##
    # Applied pressure on the internal surface.
    press =   1.0*phun("MEGA*PA");

    ##
    # Analytical solutions.   Radial stress:
    radial_stress(r)  = press*a.^2/(b^2-a^2).*(1-b^2./r.^2);
    ##
    # Circumferential (hoop) stress:
    hoop_stress(r) = press*a.^2/(b^2-a^2).*(1+b^2./r.^2);

    ##
    # Radial displacement:
    radial_displacement(r) = press*a^2*(1+nu)*(b^2+r.^2*(1-2*nu))/(E*(b^2-a^2).*r);;

    ##
    # Therefore the radial displacement of the loaded surface will be:
    urex  =  radial_displacement(a);

    ##
    # The mesh parameters: The numbers of element edges axially,
    # and through the thickness of the pipe wall (radially).

    na = 1; nt = 10;

    ##
    # Note that the material object needs to be created with the proper
    # model-dimension reduction in effect.  In this case that is the axial symmetry
    # assumption.
    MR = DeforModelRed2DAxisymm
    axisymmetric = true

    # Create the mesh and initialize the geometry.  First we are going
    # to construct the block of elements with the first coordinate
    # corresponding to the thickness in the radial direction, and the second
    # coordinate is the thickness in the axial direction.
    fens,fes  =   Q8block(b-a, t, nt, na);

    # Extract the boundary  and mark the finite elements on the
    # interior surface.
    bdryfes = meshboundary(fes);

    bcl = selectelem(fens, bdryfes, box=[0.,0.,-Inf,Inf], inflate=tolerance);
    internal_fenids= connectednodes(subset(bdryfes,bcl));
    # Now  shape the block  into  the actual wedge piece of the pipe.
    for i=1:count(fens)
        fens.xyz[i,:] = fens.xyz[i,:] + [a; 0.0];
    end

    # now we create the geometry and displacement fields
    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),2)) # displacement field

    # The plane-strain condition in the axial direction  is specified by selecting nodes
    # on the plane y=0 and y=t.
    l1 = selectnode(fens; box=[-Inf Inf 0.0 0.0], inflate = tolerance)
    setebc!(u, l1, true, 2, 0.0)
    l1 = selectnode(fens; box=[-Inf Inf t t], inflate = tolerance)
    setebc!(u, l1, true, 2, 0.0)

    applyebc!(u)
    numberdofs!(u)

    # The traction boundary condition is applied in the radial
    # direction.

    el1femm =  FEMMBase(IntegData(subset(bdryfes,bcl), GaussRule(1, 3), axisymmetric))
    fi = ForceIntensity([press; 0.0]);
    F2= distribloads(el1femm, geom, u, fi, 2);

    # Property and material
    material = MatDeforElastIso(MR,  E, nu)

    femm = FEMMDeforLinear(MR, IntegData(fes, GaussRule(2, 2), axisymmetric), material)

    K = stiffness(femm, geom, u)
    #K=cholesky(K)
    U =  K\(F2)
    scattersysvec!(u,U[:])

    # Transfer the solution of the displacement to the nodes on the
    # internal cylindrical surface and convert to
    # cylindrical-coordinate displacements there.
    uv=u.values[internal_fenids,:]
    # Report the  relative displacement on the internal surface:
    println("(Approximate/true displacement) at the internal surface: $( mean(uv[:,1])/urex*100  ) %")

    # Produce a plot of the radial stress component in the cylindrical
    # coordinate system. Note that this is the usual representation of
    # stress using nodal stress field.

    fld = fieldfromintegpoints(femm, geom, u, :Cauchy, 1)

    File =  "thick_pipe_sigmax.vtk"
    vtkexportmesh(File, fens, fes; scalars=[("sigmax", fld.values)])

    # Produce a plot of the solution components in the cylindrical
    # coordinate system.

    function inspector(idat::MyIData, elnum, conn, xe,  out,  xq)
        push!(idat.r, xq[1])
        push!(idat.s, out[idat.c])
        return idat
    end

    idat = MyIData(1, FFltVec[], FFltVec[])
    idat = inspectintegpoints(femm, geom, u, collect(1:count(fes)), inspector, idat, :Cauchy)

    @pyimport matplotlib.pyplot as plt
    plt.style[:use]("seaborn-whitegrid")
    fig = plt.figure()
    ax = plt.axes()
    # Plot the analytical solution.
    r = linspace(a,b,100);
    ax[:plot](r, radial_stress(r), linestyle="solid", marker=:o, label="Analytical")
    # Plot the computed  integration-point data
    ax[:plot](idat.r, idat.s, linestyle="none", marker=:o, label="FEA")
    # plt.xlim(0.0, 1.0)
    # plt.axis("equal")
    plt.legend()
    ax[:set_xlabel]("Radial distance")
    ax[:set_ylabel]("Radial stress")
    ax[:grid](linestyle="--", linewidth=0.5, color=".25", zorder=-10)
    plt.show()

    # ##
    # # *Regular quadratic triangle*
    # ##
    # # We start with the workhorse of most commonly used finite element
    # # packages, the quadratic triangle. Similarly to the quadratic
    # # tetrahedron in the 3-D version of this tutorial (pub_thick_pipe), the
    # # stress is polluted with oscillations. Definitely not as bad as in the
    # # plane-strain simulations when the triangles were distorted into
    # # shapes with curved edges, but the disturbances are there.
    # description ='T6';# tetrahedron
    # mf =@T6_block;
    # femmf =@(fes)femm_deformation_linear(struct('fes',fes,...
    #     'material',mater,'integration_rule',tri_rule(struct('npts',3))));
    # surface_integration_rule=gauss_rule(struct('dim',1, 'order', 3));
    # execute_simulation (description, mf, femmf, surface_integration_rule);

    # ##
    # # The same remedy of selective reduced integration as in full 3-D models
    # # will also work here.  Is demonstrated by the simulation with the
    # # selective reduced integration quadratic triangle.
    # ##
    # # *Selective reduced
    # # integration quadratic triangle*
    # description ='T6-SRI';
    # mf =@T6_block;
    # femmf =@(fes)femm_deformation_linear_sri(struct('fes',fes,...
    #     'material',mater,...
    #     'integration_rule_volumetric',tri_rule(struct('npts',1)),...
    #     'integration_rule_deviatoric',tri_rule(struct('npts',3))));
    # surface_integration_rule=gauss_rule(struct('dim',1, 'order', 3));
    # execute_simulation (description, mf, femmf, surface_integration_rule);

    # ##
    # # The selective reduced integration works very well with the T6 triangle.
    # ##
    # # An element that is often used in these situations is the uniformly
    # # under integrated serendipity (8-node) quadrilateral.
    # ##
    # # *Reduced integration serendipity quadrilateral*
    # ##
    # #  The same finite
    # # element model machine as above is used, and the integration is the 2
    # # x 2 Gauss rule (one order lower than that required for full
    # # integration which would be 3 x 3).
    # description ='Q8R';
    # mf =@Q8_block;
    # femmf =@(fes)femm_deformation_linear(struct('fes',fes,...
    #     'material',mater,...
    #     'integration_rule',gauss_rule(struct('dim',2, 'order',2))));
    # surface_integration_rule=gauss_rule(struct('dim',1, 'order',3));
    # execute_simulation (description, mf, femmf, surface_integration_rule);

    # ##
    # # *Full integration serendipity quadrilateral*
    # ##
    # # Using the full-integration Gauss rule of 3 x 3 points clearly leads to disaster.
    # description ='Q8';
    # mf =@Q8_block;
    # femmf =@(fes)femm_deformation_linear(struct('fes',fes,...
    #     'material',mater,...
    #     'integration_rule',gauss_rule(struct('dim',2, 'order',3))));
    # surface_integration_rule=gauss_rule(struct('dim',1, 'order',3));
    # execute_simulation (description, mf, femmf, surface_integration_rule);

    ##
    # The stress is now totally unacceptable.

    ## Discussion
    #
    ##
    # The axially symmetric model is clearly very effective
    # computationally, as the size is much reduced compared to the 3-D
    # model.  In conjunction with  uniform or selective reduced integration
    # it can be very accurate as well.

    #pub_thick_pipe_axi()
    # end

end # thick_pipe_axi

function thick_pipe_ps()

    ## Thick pipe with internal pressure: plane strain
    #

    ##
    # Link to the  <matlab:edit('pub_thick_pipe_ps') m-file>.

    ## Description
    ##
    # This is a simple modification of the full three-dimensional simulation of
    # the tutorial pub_thick_pipe takes advantage of the plane-strain model
    # reduction procedure.
    ##
    # An infinitely long thick walled cylindrical pipe
    # with inner boundary radius of 3 mm and outer boundary radius of 9 mm is
    # subjected to an internal pressure of 1.0 MPa. A wedge   with thickness of
    # 2 mm and a 90-degree angle sector is considered for the finite element
    # analysis. The material properties are taken as  isotropic linear elastic
    # with $E=1000$ MPa and $\nu=0.4999$ to represent nearly incompressible
    # behavior. This problem has been proposed to by MacNeal and Harder as a
    # test of an element's ability to represent the  response of a nearly
    # incompressible material. The plane-strain condition is assumed in the
    # axial direction of the pipe which together with the radial symmetry
    # confines the material in all but the radial direction and therefore
    # amplifies the numerical difficulties associated with the confinement of
    # the nearly incompressible material.
    ##
    # There is an analytical solution to this problem. Timoshenko and Goodier
    # presented the original solution of Lame in their textbook. We are going
    # to compare with  both the stress distribution (radial and hoop stresses)
    # and the displacement of the inner  cylindrical surface.

    ##
    #
    #
    #

    #
    #

    #
Figure 1. Definition of the geometry of the internally pressurized thick pipe

    #

    #

    ##
    # References:
    #
    # # Macneal RH, Harder RL (1985) A proposed standard set of problems to test
    # finite element accuracy. Finite Elements in Analysis and Design 1: 3-20.
    #
    # # Timoshenko S. and Goodier J. N., Theory of Elasticity, McGraw-Hill, 2nd ed., 1951.

    ## Solution
    #

    ##
    # Internal radius of the pipe.
    a = 3*phun("MM");
    ##
    # External radius of the pipe.
    b = 9*phun("MM");
    ##
    # Thickness of the slice.
    t = 2*phun("MM");

    ##
    # Geometrical tolerance.
    tolerance  =a/10000.;
    ##
    # Young's modulus and Poisson's ratio.
    E = 1000*phun("MEGA*PA");
    nu = 0.499;
    ##
    # Applied pressure on the internal surface.
    press = 1.0*phun("MEGA*PA");

    ##
    # Analytical solutions.   Radial stress:
    radial_stress(r) =press*a.^2/(b^2-a^2).*(1-(b^2)./r.^2);
    ##
    # Circumferential (hoop) stress:
    hoop_stress(r)=press*a.^2/(b^2-a^2).*(1+(b^2)./r.^2);

    ##
    # Radial displacement:
    radial_displacement(r)=press*a^2*(1+nu)*(b^2+r.^2*(1-2*nu))/(E*(b^2-a^2).*r);;

    ##
    # Therefore the radial displacement of the loaded surface will be:
    urex = radial_displacement(a);

    ##
    # The mesh parameters: The numbers of element edges axially,
    # and through the thickness of the pipe wall (radially).

    nc=3; nt=3;

    ##
    # Note that the material object needs to be created with the proper
    # model-dimension reduction in mind.  In this case that is the axial symmetry
    # assumption.
    MR = DeforModelRed2DStrain

    # Create the mesh and initialize the geometry.  First we are going
    # to construct the block of elements with the first coordinate
    # corresponding to the angle, and the second
    # coordinate is the thickness in the radial direction.
    anglrrange = 90.0/180*pi;
    fens,fes =  Q8block(anglrrange, b-a, nc, nt);

    # Extract the boundary  and mark the finite elements on the
    # interior surface.
    bdryfes = meshboundary(fes);
    bcl = selectelem(fens, bdryfes, box=[-Inf,Inf,0.,0.], inflate=tolerance);
    internal_fenids = connectednodes(subset(bdryfes,bcl));
    # Now  shape the block  into  the actual wedge piece of the pipe.
    ayr=fens.xyz;
    for i=1:count(fens)
        angl=ayr[i,1]; r=a+ayr[i,2];
        fens.xyz[i,:] = [r*sin(angl),(r*cos(angl))];
    end

    # now we create the geometry and displacement fields
    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),2)) # displacement field

    # The symmetry boundary condition  is specified by selecting nodes
    # on the plane x=0.
    l1 = selectnode(fens; box=[0.0 0.0 -Inf Inf], inflate = tolerance)
    setebc!(u, l1, true, 1, 0.0)
    # The second symmetry boundary condition is specified by selecting
    # nodes on the plane y=0.
    l1 = selectnode(fens; box=[-Inf Inf 0.0 0.0], inflate = tolerance)
    setebc!(u, l1, true, 2, 0.0)

    applyebc!(u)
    numberdofs!(u)

    # The traction boundary condition is applied in the radial
    # direction.

    el1femm =  FEMMBase(IntegData(subset(bdryfes,bcl), GaussRule(1, 3)))
    function pressureloading!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        copy!(forceout, XYZ/norm(XYZ)*press)
        return forceout
    end
    fi = ForceIntensity(FFlt, 2, pressureloading!); # pressure normal to the internal cylindrical surface
    F2 = distribloads(el1femm, geom, u, fi, 2);

    # Property and material
    material = MatDeforElastIso(MR, E, nu)

    femm = FEMMDeforLinear(MR, IntegData(fes, GaussRule(2, 2)), material)

    K =stiffness(femm, geom, u)
    #K=cholesky(K)
    U=  K\(F2)
    scattersysvec!(u, U[:])

    # Transfer the solution of the displacement to the nodes on the
    # internal cylindrical surface and convert to
    # cylindrical-coordinate displacements there.
    uv = u.values[internal_fenids,:];
    ur = zeros(FFlt,length(internal_fenids));

    for j = 1:length(internal_fenids)
        n = fens.xyz[internal_fenids[j],:];
        n = n'/norm(n);# normal to the cylindrical internal surface
        ur[j] = dot(vec(uv[j,:]),vec(n));
    end

    # Report the  relative displacement on the internal surface:
    println("(Approximate/true displacement) at the internal surface: $( mean(ur)/urex*100  ) %")

    # Produce a plot of the radial stress component in the cylindrical
    # coordinate system. Note that this is the usual representation of
    # stress using nodal stress field.

    fld = fieldfromintegpoints(femm, geom, u, :Cauchy, 1)

    File =  "thick_pipe_sigmax.vtk"
    vtkexportmesh(File, fens, fes; scalars=[("sigmax", fld.values)])

    # Produce a plot of the solution components in the cylindrical
    # coordinate system.
    # Plot the analytical solution.

    function inspector(idat::MyIData, elnum, conn, xe,  out,  xq)
        function outputRm(c)
            theNormal=c;
            r=norm(theNormal);# distance from the axis of symmetry
            theNormal =theNormal/r;# compute the unit normal vector
            e1p=[theNormal';0.];# local cylind. coordinate  basis vectors
            e3p=[0.,0.,1.]';# this one points along the axis of the cylinder
            e2p=cross(vec(e3p),vec(e1p));# this one is along the hoop direction
            R= [vec(e1p) vec(e2p) vec(e3p)];# transformation matrix for the stress
            return R
        end
        Rm=outputRm(xq)
        tm=zeros(FFlt,3,3)
        stress4vto3x3t!(tm, out);# stress in global XYZ
        tpm = Rm'*tm*Rm;#  stress matrix in cylindrical coordinates
        sp=zeros(FFlt,6)
        stress3x3tto6v!(sp, tpm);# stress vector in cylindr. coord.
        push!(idat.r,norm(xq))
        push!(idat.s,sp[idat.c])
        return idat
    end

    idat = MyIData(1, FFltVec[], FFltVec[])
    idat = inspectintegpoints(femm, geom, u, collect(1:count(fes)),
    inspector, idat, :Cauchy)
    # show(idat)

    @pyimport matplotlib.pyplot as plt
    plt.style[:use]("seaborn-whitegrid")
    fig = plt.figure()
    ax = plt.axes()
    # Plot the analytical solution.
    r = linspace(a,b,100);
    ax[:plot](r, radial_stress(r), linestyle="solid", marker=:o, label="Analytical")
    # Plot the computed  integration-point data
    ax[:plot](idat.r, idat.s, linestyle="none", marker=:o, label="FEA")
    # plt.xlim(0.0, 1.0)
    # plt.axis("equal")
    plt.legend()
    ax[:set_xlabel]("Radial distance")
    ax[:set_ylabel]("Radial stress")
    ax[:grid](linestyle="--", linewidth=0.5, color=".25", zorder=-10)
    plt.show()

end # thick_pipe_ps

function thick_pipe_ps_T6()

    ## Thick pipe with internal pressure: plane strain
    #

    ##
    # Link to the  <matlab:edit('pub_thick_pipe_ps') m-file>.

    ## Description
    ##
    # This is a simple modification of the full three-dimensional simulation of
    # the tutorial pub_thick_pipe takes advantage of the plane-strain model
    # reduction procedure.
    ##
    # An infinitely long thick walled cylindrical pipe
    # with inner boundary radius of 3 mm and outer boundary radius of 9 mm is
    # subjected to an internal pressure of 1.0 MPa. A wedge   with thickness of
    # 2 mm and a 90-degree angle sector is considered for the finite element
    # analysis. The material properties are taken as  isotropic linear elastic
    # with $E=1000$ MPa and $\nu=0.4999$ to represent nearly incompressible
    # behavior. This problem has been proposed to by MacNeal and Harder as a
    # test of an element's ability to represent the  response of a nearly
    # incompressible material. The plane-strain condition is assumed in the
    # axial direction of the pipe which together with the radial symmetry
    # confines the material in all but the radial direction and therefore
    # amplifies the numerical difficulties associated with the confinement of
    # the nearly incompressible material.
    ##
    # There is an analytical solution to this problem. Timoshenko and Goodier
    # presented the original solution of Lame in their textbook. We are going
    # to compare with  both the stress distribution (radial and hoop stresses)
    # and the displacement of the inner  cylindrical surface.

    ##
    #
    #
    #

    #
    #

    #
Figure 1. Definition of the geometry of the internally pressurized thick pipe

    #

    #

    ##
    # References:
    #
    # # Macneal RH, Harder RL (1985) A proposed standard set of problems to test
    # finite element accuracy. Finite Elements in Analysis and Design 1: 3-20.
    #
    # # Timoshenko S. and Goodier J. N., Theory of Elasticity, McGraw-Hill, 2nd ed., 1951.

    ## Solution
    #

    ##
    # Internal radius of the pipe.
    a = 3*phun("MM");
    ##
    # External radius of the pipe.
    b = 9*phun("MM");
    ##
    # Thickness of the slice.
    t = 2*phun("MM");

    ##
    # Geometrical tolerance.
    tolerance  =a/10000.;
    ##
    # Young's modulus and Poisson's ratio.
    E = 1000*phun("MEGA*PA");
    nu = 0.499;
    ##
    # Applied pressure on the internal surface.
    press = 1.0*phun("MEGA*PA");

    ##
    # Analytical solutions.   Radial stress:
    radial_stress(r) =press*a.^2/(b^2-a^2).*(1-(b^2)./r.^2);
    ##
    # Circumferential (hoop) stress:
    hoop_stress(r)=press*a.^2/(b^2-a^2).*(1+(b^2)./r.^2);

    ##
    # Radial displacement:
    radial_displacement(r)=press*a^2*(1+nu)*(b^2+r.^2*(1-2*nu))/(E*(b^2-a^2).*r);;

    ##
    # Therefore the radial displacement of the loaded surface will be:
    urex = radial_displacement(a);

    ##
    # The mesh parameters: The numbers of element edges axially,
    # and through the thickness of the pipe wall (radially).

    nc=3; nt=3;

    ##
    # Note that the material object needs to be created with the proper
    # model-dimension reduction in mind.  In this case that is the axial symmetry
    # assumption.
    MR = DeforModelRed2DStrain

    # Create the mesh and initialize the geometry.  First we are going
    # to construct the block of elements with the first coordinate
    # corresponding to the angle, and the second
    # coordinate is the thickness in the radial direction.
    anglrrange = 90.0/180*pi;
    fens,fes =  T6block(anglrrange, b-a, nc, nt);

    # Extract the boundary  and mark the finite elements on the
    # interior surface.
    bdryfes = meshboundary(fes);
    bcl = selectelem(fens, bdryfes, box=[-Inf,Inf,0.,0.], inflate=tolerance);
    internal_fenids = connectednodes(subset(bdryfes,bcl));
    # Now  shape the block  into  the actual wedge piece of the pipe.
    ayr=fens.xyz;
    for i=1:count(fens)
        angl=ayr[i,1]; r=a+ayr[i,2];
        fens.xyz[i,:] = [r*sin(angl),(r*cos(angl))];
    end

    # now we create the geometry and displacement fields
    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),2)) # displacement field

    # The symmetry boundary condition  is specified by selecting nodes
    # on the plane x=0.
    l1 = selectnode(fens; box=[0.0 0.0 -Inf Inf], inflate = tolerance)
    setebc!(u, l1, true, 1, 0.0)
    # The second symmetry boundary condition is specified by selecting
    # nodes on the plane y=0.
    l1 = selectnode(fens; box=[-Inf Inf 0.0 0.0], inflate = tolerance)
    setebc!(u, l1, true, 2, 0.0)

    applyebc!(u)
    numberdofs!(u)

    # The traction boundary condition is applied in the radial
    # direction.

    el1femm =  FEMMBase(IntegData(subset(bdryfes,bcl), GaussRule(1, 3)))
    function pressureloading!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        copy!(forceout, XYZ/norm(XYZ)*press)
        return forceout
    end
    fi = ForceIntensity(FFlt, 2, pressureloading!); # pressure normal to the internal cylindrical surface
    F2 = distribloads(el1femm, geom, u, fi, 2);

    # Property and material
    material = MatDeforElastIso(MR, E, nu)

    femm = FEMMDeforLinear(MR, IntegData(fes, TriRule(3)), material)

    K =stiffness(femm, geom, u)
    #K=cholesky(K)
    U=  K\(F2)
    scattersysvec!(u, U[:])

    # Transfer the solution of the displacement to the nodes on the
    # internal cylindrical surface and convert to
    # cylindrical-coordinate displacements there.
    uv = u.values[internal_fenids,:];
    ur = zeros(FFlt,length(internal_fenids));

    for j = 1:length(internal_fenids)
        n = fens.xyz[internal_fenids[j],:];
        n = n'/norm(n);# normal to the cylindrical internal surface
        ur[j] = dot(vec(uv[j,:]),vec(n));
    end

    # Report the  relative displacement on the internal surface:
    println("(Approximate/true displacement) at the internal surface: $( mean(ur)/urex*100  ) %")

    # Produce a plot of the radial stress component in the cylindrical
    # coordinate system. Note that this is the usual representation of
    # stress using nodal stress field.

    fld = fieldfromintegpoints(femm, geom, u, :Cauchy, 1)

    File =  "thick_pipe_sigmax.vtk"
    vtkexportmesh(File, fens, fes; scalars=[("sigmax", fld.values)])

    # Produce a plot of the solution components in the cylindrical
    # coordinate system.
    # Plot the analytical solution.

    function inspector(idat::MyIData, elnum, conn, xe,  out,  xq)
        function outputRm(c)
            theNormal=c;
            r=norm(theNormal);# distance from the axis of symmetry
            theNormal =theNormal/r;# compute the unit normal vector
            e1p=[theNormal';0.];# local cylind. coordinate  basis vectors
            e3p=[0.,0.,1.]';# this one points along the axis of the cylinder
            e2p=cross(vec(e3p),vec(e1p));# this one is along the hoop direction
            R= [vec(e1p) vec(e2p) vec(e3p)];# transformation matrix for the stress
            return R
        end
        Rm=outputRm(xq)
        tm=zeros(FFlt,3,3)
        stress4vto3x3t!(tm, out);# stress in global XYZ
        tpm = Rm'*tm*Rm;#  stress matrix in cylindrical coordinates
        sp=zeros(FFlt,6)
        stress3x3tto6v!(sp, tpm);# stress vector in cylindr. coord.
        push!(idat.r,norm(xq))
        push!(idat.s,sp[idat.c])
        return idat
    end

    idat = MyIData(1, FFltVec[], FFltVec[])
    idat = inspectintegpoints(femm, geom, u, collect(1:count(fes)),
    inspector, idat, :Cauchy)
    # show(idat)

    @pyimport matplotlib.pyplot as plt
    plt.style[:use]("seaborn-whitegrid")
    fig = plt.figure()
    ax = plt.axes()
    # Plot the analytical solution.
    r = linspace(a,b,100);
    ax[:plot](r, radial_stress(r), linestyle="solid", marker=:o, label="Analytical")
    # Plot the computed  integration-point data
    ax[:plot](idat.r, idat.s, linestyle="none", marker=:o, label="FEA")
    # plt.xlim(0.0, 1.0)
    # plt.axis("equal")
    plt.legend()
    ax[:set_xlabel]("Radial distance")
    ax[:set_ylabel]("Radial stress")
    ax[:grid](linestyle="--", linewidth=0.5, color=".25", zorder=-10)
    plt.show()

    #pub_thick_pipe_ps()

end # thick_pipe_ps_T6

function allrun()
    println("#####################################################")
    println("# thick_pipe_axi ")
    thick_pipe_axi()
    println("#####################################################")
    println("# thick_pipe_ps ")
    thick_pipe_ps()
    println("#####################################################")
    println("# thick_pipe_ps_T6 ")
    thick_pipe_ps_T6()
    return true
end # function allrun

end # module thick_pipe_examples


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule sgtsnepi_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("sgtsnepi")
JLLWrappers.@generate_main_file("sgtsnepi", UUID("c2c51ba6-9464-585c-93d5-ba434ab08fad"))
end  # module sgtsnepi_jll

ext: jl


text:
module SymArrays

export SymArray, contract, contract!, SymArr_ifsym

include("symarray.jl")

include("contractions.jl")

end # module


text:
### Working with data

# usage of read with Char
julia> read(STDIN, Char)
j
'j': ASCII/Unicode U+006a (category Ll: Letter, lowercase)

# usage of read with Bool
julia> read(STDIN, Bool)
true
true

# usage of read with Int8
julia> read(STDIN, Int8)
23
50

# method signature for Bool Type
julia> @which read(STDIN,Bool)
read(s::IO, ::Type{Bool}) in Base at io.jl:370

# method signature for Char Type
julia> @which read(STDIN,Char)
read(s::IO, ::Type{Char}) in Base at io.jl:403

# method signature for Int8 Type
julia> @which read(STDIN,Int8)
read(s::IO, ::Type{Int8}) in Base at io.jl:365

# Or Alternatively
julia> for val in [:Bool, :Char, :Int8]
          @eval println(@which read(STDIN, $val))
       end
read(s::IO, ::Type{Bool}) in Base at io.jl:370
read(s::IO, ::Type{Char}) in Base at io.jl:403
read(s::IO, ::Type{Int8}) in Base at io.jl:365

julia> statement = readline()
julia is fast!
"julia is fast!"

julia> readline(chomp = false)
julia
"julia\n"

julia> methods(readline)
# 4 methods for generic function "readline":
readline(s::IOStream; chomp) in Base at iostream.jl:234
readline() in Base at io.jl:190
readline(filename::AbstractString; chomp) in Base at io.jl:184
readline(s::IO; chomp) in Base at io.jl:190

julia> number = parse(Int64, readline())
23

julia> println(number)
23

# taking input from user and performing actions on it
# take input from user
julia> numbers = readline()
23 45 67 89
"23 45 67 89"

# use the split function to convert string to array
julia> split(numbers)
4-element Array{SubString{String},1}:
 "23"
 "45"
 "67"
 "89"

# the easy way
julia> for item in split(numbers)
           println(parse(Int64, item))
       end
23
45
67
89

# the array comprehension way
julia> [parse(Int64, item) for item in split(numbers)]
4-element Array{Int64,1}:
 23
 45
 67
 89

### Working with Text Files
# Mode Description
# –––– –––––––––––––––––––––––––––––
# r    read
# r+   read, write
# w    write, create, truncate
# w+   read, write, create, truncate
# a    write, create, append
# a+   read, write, create, append

julia> methods(open)
# 8 methods for generic function "open":
open(fname::AbstractString) in Base at iostream.jl:113
open(fname::AbstractString, rd::Bool, wr::Bool, cr::Bool, tr::Bool, ff::Bool) in Base at iostream.jl:103
open(fname::AbstractString, mode::AbstractString) in Base at iostream.jl:132
open(f::Function, cmds::Base.AbstractCmd, args...) in Base at process.jl:599
open(f::Function, args...) in Base at iostream.jl:150
open(cmds::Base.AbstractCmd) in Base at process.jl:575
open(cmds::Base.AbstractCmd, mode::AbstractString) in Base at process.jl:575
open(cmds::Base.AbstractCmd, mode::AbstractString, other::Union{Base.FileRedirect, IO, RawFD}) in Base at process.jl:575

# sample file
shell> cat sample.txt
Julia is a great and powerful language
It is homoiconic.
It supports static as well as dynamic typing.
Julia walks like python runs like C.
Julia uses multiple dispatch.

# inside julia,
# read file from the device folder in the current path
julia> file = open("sample.txt")
IOStream(<file sample.txt>)

# oread the contents of file
julia> file_data = readlines(file)
5-element Array{String,1}:
 "Julia is a great and powerful language"
 "It is homoiconic."
 "It supports static as well as dynamic typing."
 "Julia walks like python runs like C."
 "Julia uses multiple dispatch."

# running an enumerate over the file provides us with a counter too!
julia> enumerate(file_data)
Enumerate{Array{String,1}}(String["Julia is a great and powerful language", "It is homoiconic.", "It supports static as well as dynamic typing.", "Julia walks like python runs like C.", "Julia uses multiple dispatch."])

# using enumerate to our advantage to get the length of each line
julia> for lines in enumerate(file_data)
           println(lines[1],"-> ", lines[2])
       end
1-> Julia is a great and powerful language
2-> It is homoiconic.
3-> It supports static as well as dynamic typing.
4-> Julia walks like python runs like C.
5-> Julia uses multiple dispatch.

# convert everything to uppercase
julia> for line in file_data
           println(uppercase(line))
       end
JULIA IS A GREAT AND POWERFUL LANGUAGE
IT IS HOMOICONIC.
IT SUPPORTS STATIC AS WELL AS DYNAMIC TYPING.
JULIA WALKS LIKE PYTHON RUNS LIKE C.
JULIA USES MULTIPLE DISPATCH.

# reverse each line!!
julia> for line in file_data
           println(reverse(line))
       end
egaugnal lufrewop dna taerg a si ailuJ
.cinociomoh si tI
.gnipyt cimanyd sa llew sa citats stroppus tI
.C ekil snur nohtyp ekil sklaw ailuJ
.hctapsid elpitlum sesu ailuJ

# to simply count the number of lines in a file
julia> countlines("sample.txt")
5

# print the first line of the file
julia> first(file_data)
"Julia is a great and powerful language"

# print last line of the file
julia> last(file_data)
"Julia uses multiple dispatch."

###Working with CSV and Delimited File Formats

julia> csvfile = readcsv("sample.csv")
5×3 Array{Any,2}:
 1  "James"       "UK"
 2  "Lora"        "UK"
 3  "Raj"         "India"
 4  "Rangnatham"  "Sri lanka"
 5  "Azhar"       "Bangladesh"

# getting only the names of the people
julia> csvfile[:,2]
5-element Array{Any,1}:
 "James"
 "Lora"
 "Raj"
 "Rangnatham"
 "Azhar"

# getting only the top 3 data
julia> csvfile[1:3,:]
3×3 Array{Any,2}:
 1  "James"  "UK"
 2  "Lora"   "UK"
 3  "Raj"    "India"

# reverse sorting the rows
julia> sortrows(csvfile, rev=true)
5×3 Array{Any,2}:
 5  "Azhar"       "Bangladesh"
 4  "Rangnatham"  "Sri lanka"
 3  "Raj"         "India"
 2  "Lora"        "UK"
 1  "James"       "UK"

# A pipe seperated file with the same data as above
shell> cat sample.psv
1|"James"|"UK"
2|"Lora"|"UK"
3|"Raj"|"India"
4|"Rangnatham"|"Sri lanka"
5|"Azhar"|"Bangladesh"

# using '|' as delimiter
julia> readdlm("sample.psv",'|')
5×3 Array{Any,2}:
 1  "James"       "UK"
 2  "Lora"        "UK"
 3  "Raj"         "India"
 4  "Rangnatham"  "Sri lanka"
 5  "Azhar"       "Bangladesh"

### Working with DataFrames
# install the package
julia> Pkg.add("DataFrames")

# first time load takes times due to precompilation
julia> using DataFrames

# check the functions provided by DataFrames
julia> names(DataFrames)
254-element Array{Symbol,1}:
 :&
 Symbol("@csv2_str")
 Symbol("@csv_str")
 Symbol("@data")
 Symbol("@formula")
 Symbol("@pdata")
 Symbol("@tsv_str")
 Symbol("@wsv_str")
 Symbol("@~")
 :AIC
 :AICc
 :AbstractContrasts
 :AbstractDataArray
 :AbstractDataFrame
 :AbstractDataMatrix
 :AbstractDataVector
 ⋮
 :winsor
 :winsor!
 :wmean
 :wmean!
 :wmedian
 :wquantile
 :writetable
 :wsample
 :wsample!
 :wsum
 :wsum!
 :xor
 :zscore
 :zscore!
 :|

julia> DataMatrix{Any} == DataArray{Any, 2}
true

julia> DataVector{Any} == DataArray{Any, 1}
true

# NA values
julia> a = [1,2,3,nothing,5,6]
6-element Array{Any,1}:
 1
 2
 3
  nothing
 5
 6

# try to access the element
julia> a[4]

# check the data type
julia> typeof(a[4])
Void

julia> using DataArrays

julia> a = DataArray([1.1, 2.2, 3.3, 4.4, 5.5, 6.6])
6-element DataArrays.DataArray{Float64,1}:
 1.1
 2.2
 3.3
 4.4
 5.5
 6.6

julia> a[1] = NA
NA

julia> a
6-element DataArrays.DataArray{Float64,1}:
  NA
 2.2
 3.3
 4.4
 5.5
 6.6

julia> typeof(a[1])
DataArrays.NAtype

julia> true || a
true

julia> true && a
6-element DataArrays.DataArray{Float64,1}:
  NA
 2.2
 3.3
 4.4
 5.5
 6.6

julia> mean(a)
NA

julia> mean(a[2:length(a)])
4.4

# DataArrays
# creating a 2D data matrix
julia> data_matrix = DataArray([1 2 3;4 5 6])
2×3 DataArrays.DataArray{Int64,2}:
 1  2  3
 4  5  6

# Creating a 1D vector
julia> data_vector = DataArray([1 2 3 4 5 6 7 8])
1×8 DataArrays.DataArray{Int64,2}:
 1  2  3  4  5  6  7  8

# Creating a column major 1D array
julia> data_vector_column_major = DataArray([1,2,3,4])
4-element DataArrays.DataArray{Int64,1}:
 1
 2
 3
 4

# assigning the first value as NA
julia> data_vector_column_major[1] = NA
NA

# dropping NA values
julia> data_vector_column_major
4-element DataArrays.DataArray{Int64,1}:
  NA
 2
 3
 4

# DataFrames
julia> dframe = DataFrame(Names = ["John","Ajay"], Age = [27,28])
2×2 DataFrames.DataFrame
│ Row │ Names  │ Age │
├─────┼────────┼─────┤
│ 1   │ "John" │ 27  │
│ 2   │ "Ajay" │ 28  │

julia> dframe.columns[1]
2-element DataArrays.DataArray{String,1}:
 "John"
 "Ajay"

julia> typeof(dframe.columns[1])
DataArrays.DataArray{String,1}

### Linear Algebra and Differential Calculus

julia> A = rand(3,3)
3×3 Array{Float64,2}:
 0.821807   0.828687  0.974031
 0.996824   0.805663  0.274284
 0.0341033  0.224237  0.39982

julia> ones(5)
5-element Array{Float64,1}:
 1.0
 1.0
 1.0
 1.0
 1.0

# Vector
julia> Vector{Float64} == Array{Float64,1}
true

# Matrix
julia> Matrix{Float64} == Array{Float64,2}
true

# Multiplication
# creat a random matrix of 3x3 dimension
julia> A
3×3 Array{Float64,2}:
 0.673465  0.880229  0.100458
 0.752117  0.545464  0.0180286
 0.531316  0.221628  0.179626

# Easily, we could do
julia> b = 2 * A
3×3 Array{Float64,2}:
 1.34693  1.76046   0.200917
 1.50423  1.09093   0.0360573
 1.06263  0.443255  0.359252

# Julia's advantage! Omit the '*' while multiplication
julia> b = 2A
3×3 Array{Float64,2}:
 1.34693  1.76046   0.200917
 1.50423  1.09093   0.0360573
 1.06263  0.443255  0.359252

# Matrix operations
# Matrix Transpose
julia> transpose_of_A = A'
3×3 Array{Float64,2}:
 0.673465  0.752117   0.531316
 0.880229  0.545464   0.221628
 0.100458  0.0180286  0.179626

# Matrix inversion
julia> X = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> inv(X)
2×2 Array{Float64,2}:
 -2.0   1.0
  1.5  -0.5

# Matrix Determinant
julia> M = [10 10; 10 10]
2×2 Array{Int64,2}:
 10  10
 10  10

julia> det(M)
0.0

# Matrix Multiplication
julia> M * M'
2×2 Array{Int64,2}:
 200  200
 200  200

# Matrix Division
julia> div.(M, 10)
2×2 Array{Int64,2}:
 1  1
 1  1

# Eigenvalues of Matrix
julia> eigvals(M)
2-element Array{Float64,1}:
  0.0
 20.0

# linear equations
# define the value of x
julia> x = 5
5

# solve the equation
julia> equation = 3x^2 + 4x + 3
98

# Differential calculus
julia> Pkg.add("Calculus")
julia> using Calculus

# list of functions supported by Calculus
julia> names(Calculus)
22-element Array{Symbol,1}:
 Symbol("@sexpr")
 :AbstractVariable
 :BasicVariable
 :Calculus
 :SymbolParameter
 :Symbolic
 :SymbolicVariable
 :check_derivative
 :check_gradient
 :check_hessian
 :check_second_derivative
 :deparse
 :derivative
 :differentiate
 :hessian
 :integrate
 :jacobian
 :processExpr
 :second_derivative
 :simplify
 :symbolic_derivative_bessel_list
 :symbolic_derivatives_1arg

julia> f(x) = sin(x)
f (generic function with 1 method)

julia> f'(1.0) - cos(1.0)
-5.036193684304635e-12

julia> f''(1.0) - (-sin(1.0))
-6.647716624952338e-7

julia> f'''(1.0) - (-cos(1.0))
0.11809095011119602

julia> differentiate("cos(x) + sin(x) + exp(-x) * cos(x)", :x)
:(-(sin(x)) + cos(x) + (-(exp(-x)) * cos(x) + exp(-x) * -(sin(x))))

julia> differentiate("cos(x) + sin(y) + exp(-x) * cos(y)", [:x, :y])
2-element Array{Any,1}:
 :(-(sin(x)) + -(exp(-x)) * cos(y))
 :(cos(y) + exp(-x) * -(sin(y)))

### Statistics

julia> x = [10,20,30,40,50]
5-element Array{Int64,1}:
 10
 20
 30
 40
 50

# computing the mean
julia> mean(x)
30.0

# computing the median
julia> median(x)
30.0

# computing the sum
julia> sum(x)
150

# computing the standard deviation
julia> std(x)
15.811388300841896

# computing the variance
julia> var(x)
250.0

# cummax() : To find the cumulative maximum.
# cummin() : To find the cumulative minimum.
# cumsum() : To find the cumulative maximum.
# cumprod() : To find the cumulative product.

# old implementations
julia> cummax(x)
5-element Array{Int64,1}:
 10
 20
 30
 40
 50

julia> cummin(x)
5-element Array{Int64,1}:
 10
 10
 10
 10
 10

julia> cumsum(x)
5-element Array{Int64,1}:
  10
  30
  60
 100
 150

julia> cumprod(x)
5-element Array{Int64,1}:
       10
      200
     6000
   240000
 12000000

# using accumulate function
julia> accumulate(+,x)
5-element Array{Int64,1}:
  10
  30
  60
 100
 150

julia> accumulate(*,x)
5-element Array{Int64,1}:
       10
      200
     6000
   240000
 12000000

julia> accumulate(max,x)
5-element Array{Int64,1}:
 10
 20
 30
 40
 50

julia> accumulate(min,x)
5-element Array{Int64,1}:
 10
 10
 10
 10
 10

# Basic Statistics using DataFrames

julia> dframe = DataFrame(Subjects = ["Maths","Physics","Chemistry"],Marks = [90,85,95])
3×2 DataFrames.DataFrame
│ Row │ Subjects    │ Marks │
├─────┼─────────────┼───────┤
│ 1   │ "Maths"     │ 90    │
│ 2   │ "Physics"   │ 85    │
│ 3   │ "Chemistry" │ 95    │

julia> describe(dframe)
Subjects
Summary Stats:
Length:         3
Type:           String
Number Unique:  3
Number Missing: 0
% Missing:      0.000000

Marks
Summary Stats:
Mean:           90.000000
Minimum:        85.000000
1st Quartile:   87.500000
Median:         90.000000
3rd Quartile:   92.500000
Maximum:        95.000000
Length:         3
Type:           Int64
Number Missing: 0
% Missing:      0.000000

# using pandas!
julia> Pkg.add("Pandas")
julia> using Pandas

julia> pandasDataframe = Pandas.DataFrame(Dict(:Subjects => ["Maths","Physics","Chemistry"],:Marks => [90,85,95]))
   Marks   Subjects
0     90      Maths
1     85    Physics
2     95  Chemistry

julia> Pandas.describe(pandasDataframe)
       Marks
count    3.0
mean    90.0
std      5.0
min     85.0
25%     87.5
50%     90.0
75%     92.5
max     95.0

julia> pandasDataframe[:Subjects]
0        Maths
1      Physics
2    Chemistry
Name: Subjects, dtype: object

julia> pandasDataframe[:Marks]
0    90
1    85
2    95
Name: Marks, dtype: int64

# simple query operation
julia> query(pandasDataframe,:(Marks>90))
   Marks   Subjects
2     95  Chemistry

### Advanced Statistics Topics
# Distributions

julia> Pkg.add("Distributions")
julia> using Distributions

# setting the seed for normal distribution
julia> srand(123)
MersenneTwister(UInt32[0x0000007b], Base.dSFMT.DSFMT_state(Int32[1464307935, 1073116007, 222134151, 1073120226, -290652630, 1072956456, -580276323, 1073476387, 1332671753, 1073438661  …  138346874, 1073030449, 1049893279, 1073166535, -1999907543, 1597138926, -775229811, 32947490, 382, 0]), [1.23253, 1.95067, 1.54183, 1.85035, 1.28927, 1.24474, 1.40729, 1.95055, 1.39839, 1.92576  …  1.79227, 1.83391, 1.89061, 1.74502, 1.57469, 1.24833, 1.69181, 1.48955, 1.40392, 1.75348], 382)

# initialize the distribution
julia> distribution = Normal()
Distributions.Normal{Float64}(μ=0.0, σ=1.0)

# create the random array of 10 element following Normal distribution
julia> x = rand(distribution, 10)
10-element Array{Float64,1}:
  1.19027
  2.04818
  1.14265
  0.459416
 -0.396679
 -0.664713
  0.980968
 -0.0754831
  0.273815
 -0.194229

# binomial
julia> Binomial()
Distributions.Binomial{Float64}(n=1, p=0.5)

# cauchy
julia> Cauchy()
Distributions.Cauchy{Float64}(μ=0.0, σ=1.0)

# poissons
julia> Poisson()
Distributions.Poisson{Float64}(λ=1.0)

# TimeSeries

julia> Pkg.add("TimeSeries")
julia> Pkg.add("MarketData")

julia> using TimeSeries
julia> dates  = collect(Date(2017,8,1):Date(2017,8,5))
5-element Array{Date,1}:
 2017-08-01
 2017-08-02
 2017-08-03
 2017-08-04
 2017-08-05

julia> sample_time = TimeArray(dates, rand(length(dates)))
5x1 TimeSeries.TimeArray{Float64,1,Date,Array{Float64,1}} 2017-08-01 to 2017-08-05
2017-08-01 | 0.7059
2017-08-02 | 0.292
2017-08-03 | 0.2811
2017-08-04 | 0.7929
2017-08-05 | 0.2092

julia> fieldnames(sample_time)
4-element Array{Symbol,1}:
 :timestamp
 :values
 :colnames
 :meta

julia> sample_time.timestamp
5-element Array{Date,1}:
 2017-08-01
 2017-08-02
 2017-08-03
 2017-08-04
 2017-08-05

julia> sample_time.values
5-element Array{Float64,1}:
 0.70586
 0.291978
 0.281066
 0.792931
 0.20923

julia> sample_time.colnames
1-element Array{String,1}:
 ""

julia> sample_time.meta

# retrieve the first element
julia> head(sample_time)
1x1 TimeSeries.TimeArray{Float64,2,Date,Array{Float64,2}} 2017-08-01 to 2017-08-01

2017-08-01 | 0.7179

# retrieve the last element
julia> tail(sample_time)
1x1 TimeSeries.TimeArray{Float64,2,Date,Array{Float64,2}} 2017-08-05 to 2017-08-05

2017-08-05 | 0.4142

### Hypothesis Testing
julia> Pkg.add('HypothesisTests")
julia> using HypothesisTests

julia> using Distributions

julia> srand(123)
MersenneTwister(UInt32[0x0000007b], Base.dSFMT.DSFMT_state(Int32[1464307935, 1073116007, 222134151, 1073120226, -290652630, 1072956456, -580276323, 1073476387, 1332671753, 1073438661  …  138346874, 1073030449, 1049893279, 1073166535, -1999907543, 1597138926, -775229811, 32947490, 382, 0]), [1.04643, 1.18883, 1.92848, 1.52435, 1.80384, 1.20354, 1.40414, 1.04937, 1.25594, 1.66531  …  1.86043, 1.53826, 1.54179, 1.83724, 1.85687, 1.14965, 1.14176, 1.03677, 1.17827, 1.21738], 382)

julia> sampleOne = rand(Normal(), 10)
10-element Array{Float64,1}:
  1.19027
  2.04818
  1.14265
  0.459416
 -0.396679
 -0.664713
  0.980968
 -0.0754831
  0.273815
 -0.194229

julia> testOne = OneSampleTTest(sampleOne)
One sample t-test
-----------------
Population details:
    parameter of interest:   Mean
    value under h_0:         0
    point estimate:          0.47641935520300993
    95% confidence interval: (-0.13332094295432084, 1.0861596533603408)

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.11093746407653728

Details:
    number of observations:   10
    t-statistic:              1.7675319478229796
    degrees of freedom:       9
    empirical standard error: 0.2695393176852065

julia> @which OneSampleTTest(sampleOne)
HypothesisTests.OneSampleTTest(v::AbstractArray{T,1}) where T<:Real

julia> pvalue(testOne)
0.11093746407653728

# the most significant
julia> pvalue(testOne, tail=:right)
0.05546873203826864

# the least significant
julia> pvalue(testOne, tail=:left)
0.9445312679617314

julia> BinomialTest(25, 1000, 0.50)
Binomial test
-------------
Population details:
    parameter of interest:   Probability of success
    value under h_0:         0.5
    point estimate:          0.025
    95% confidence interval: (0.01624253569688223, 0.036684823051923)

Test summary:
    outcome with 95% confidence: reject h_0
    two-sided p-value:           9.125992175283306e-252

Details:
    number of observations: 1000
    number of successes:    25

### optimisation

# JuMP

julia> Pkg.add("JuMP")
julia> Pkg.add("Clp")
julia> using JuMP
julia> using Clp

# creating a model without a solver
julia> m = JuMP.Model()
Feasibility problem with:
 * 0 linear constraints
 * 0 variables
Solver is default solver

# but we want to use Clp as solver, we have
julia> m = JuMP.Model(solver = Clp.ClpSolver())
Feasibility problem with:
 * 0 linear constraints
 * 0 variables
Solver is ClpMathProg

# create optimiser.jl
using JuMP
using Clp

m = Model(solver = ClpSolver())
@variable(m, 0 <= a <= 2 )
@variable(m, 0 <= b <= 10 )

@objective(m, Max, 5a + 3*b )
@constraint(m, 1a + 5b <= 3.0 )

print(m)

status = solve(m)

println("Objective value: ", getobjectivevalue(m))
println("a = ", getvalue(a))
println("b = ", getvalue(b))

# run optimiser.jl
mymachine:user$ julia optimiser.jl
Max 5 a + 3 b
Subject to
 a + 5 b ≤ 3
 0 ≤ a ≤ 2
 0 ≤ b ≤ 10
Objective value: 10.6
a = 2.0
b = 0.2

# Convex.jl

julia> Pkg.add("Convex")
julia> Pkg.add("SCS")
julia> X = Variable(2, 2)
Variable of
size: (2, 2)
sign: Convex.NoSign()
vexity: Convex.AffineVexity()

julia> y = Variable()
Variable of
size: (1, 1)
sign: Convex.NoSign()
vexity: Convex.AffineVexity()

julia> p = minimize(vecnorm(X) + y, 2 * X <= 1, X' + y >= 1, X >= 0, y >= 0)
Problem:
minimize AbstractExpr with
head: +
size: (1, 1)
sign: Convex.NoSign()
vexity: Convex.ConvexVexity()

subject to
Constraint:
<= constraint
lhs: AbstractExpr with
head: *
size: (2, 2)
sign: Convex.NoSign()
vexity: Convex.AffineVexity()

rhs: 1
vexity: Convex.AffineVexity()
		Constraint:
>= constraint
lhs: AbstractExpr with
head: +
size: (2, 2)
sign: Convex.NoSign()
vexity: Convex.AffineVexity()

rhs: 1
vexity: Convex.AffineVexity()
		Constraint:
>= constraint
lhs: Variable of
size: (2, 2)
sign: Convex.NoSign()
vexity: Convex.AffineVexity()
rhs: 0
vexity: Convex.AffineVexity()
		Constraint:
>= constraint
lhs: Variable of
size: (1, 1)
sign: Convex.NoSign()
vexity: Convex.AffineVexity()
rhs: 0
vexity: Convex.AffineVexity()
current status: not yet solved

julia> solve!(p)
----------------------------------------------------------------------------
	SCS v1.2.6 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012-2016
----------------------------------------------------------------------------
Lin-sys: sparse-direct, nnz in A = 25
eps = 1.00e-04, alpha = 1.80, max_iters = 20000, normalize = 1, scale = 5.00
Variables n = 7, constraints m = 19
Cones:	primal zero / dual free vars: 1
	linear vars: 13
	soc vars: 5, soc blks: 1
Setup time: 1.09e-04s
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0|      inf       inf       nan       inf       inf       inf  3.97e-05
    60| 6.94e-06  8.04e-05  4.05e-05  1.00e+00  1.00e+00  2.20e-16  1.69e-04
----------------------------------------------------------------------------
Status: Solved
Timing: Solve time: 1.74e-04s
	Lin-sys: nnz in L factor: 51, avg solve time: 5.07e-07s
	Cones: avg projection time: 9.80e-08s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 5.8351e-19, dist(y, K*) = 0.0000e+00, s'y/|s||y| = 3.3410e-18
|Ax + s - b|_2 / (1 + |b|_2) = 6.9388e-06
|A'y + c|_2 / (1 + |c|_2) = 8.0432e-05
|c'x + b'y| / (1 + |c'x| + |b'y|) = 4.0465e-05
----------------------------------------------------------------------------
c'x = 1.0000, -b'y = 0.9999
============================================================================

julia> println(round(X.value, 2))
[0.0 0.0; 0.0 0.0]

julia> println(y.value)
1.0000106473513

julia> p.optval
1.0000111259613518

max_stars_repo_name: PacktPublishing/Learning-Julia

max_issues_repo_name: PacktPublishing/Learning-Julia

max_forks_repo_name: PacktPublishing/Learning-Julia


text:
# Simulation using native rk45 implementation
# Tentatively to simulate two rigid bodies constrained by one revolute joint

include("RigidBody.jl")
include("OrientationConversion.jl")
include("Joint.jl")
include("Force.jl")
include("extF.jl")
include("rk45.jl")
include("rk4.jl")

using LinearAlgebra
using Revise
using DifferentialEquations
using StaticArrays

global GravityInInertial = MVector{3}([0.0,0.0,-9.806])
# global extFList = Vector{extForces}(undef,2)

struct solSim
    # Structure to store results of Simulation for each rigid body
    r::Matrix{Float64}
    v::Matrix{Float64}
    β::Matrix{Float64}
    βdot::Matrix{Float64}
end

## RK45
function simulateRK45(tEnd::Float64,tSpan::Float64,j::Joint...;
                  g::MArray{Tuple{3},Float64,1,3}=[0.0;0.0;-9.806])#,extFVec::Vector{extForces}=Vector{extForces}(undef,1))
    # Ellipsis (...) to facilitate supplying variable number of arguments to the function
    # Initial condition (all bodies connected)
    global GravityInInertial = g
    # global extFList = extFVec
    X0 = j[1].RB1.x
    for k = 1:length(j)
        append!(X0,j[k].RB2.x)
    end

    # Using native RK implementation
    tSim, sol = rk45(tEnd,0.0,X0,j,mainDynRK)#,eps = 1e-6)
    nBodies = length(j) + 1
    ntStops = length(tSim)
    sol = reshape(sol,(14*nBodies,ntStops))

    solFinal = Vector{solSim}(undef,length(j))
    for i=1:length(j)
        rSol = Matrix{Float64}(undef,ntStops,3)
        βSol = Matrix{Float64}(undef,ntStops,4)
        vSol = Matrix{Float64}(undef,ntStops,3)
        βdotSol = Matrix{Float64}(undef,ntStops,4)
        for t=1: ntStops
            rSol[t,:] = sol[14*i+1:14*i+3,t]
            βSol[t,:] = sol[14*i+4:14*i+7,t]
            vSol[t,:] = sol[14*i+8:14*i+10,t]
            βdotSol[t,:] = sol[14*i+11:14*i+14,t]
        end
        solFinal[i] = solSim(rSol,vSol,βSol,βdotSol)
    end
    return tSim, solFinal
end

## RK4
function simulateRK4(tEnd::Float64,tSpan::Float64,j::Joint...;
                  g::MArray{Tuple{3},Float64,1,3}=[0.0;0.0;-9.806])#,extFVec::Vector{extForces}=Vector{extForces}(undef,1))
    # Ellipsis (...) to facilitate supplying variable number of arguments to the function
    # Initial condition (all bodies connected)
    global GravityInInertial = g
    # global extFList = extFVec
    X0 = j[1].RB1.x
    for k = 1:length(j)
        append!(X0,j[k].RB2.x)
    end

    # Using native RK implementation
    tSim, sol = rk4(tEnd,0.0,X0,j,mainDynRK, h = tSpan)
    nBodies = length(j) + 1
    ntStops = length(tSim)
    sol = reshape(sol,(14*nBodies,ntStops))

    solFinal = Vector{solSim}(undef,length(j))
    for i=1:length(j)
        rSol = Matrix{Float64}(undef,ntStops,3)
        vSol = Matrix{Float64}(undef,ntStops,3)
        βSol = Matrix{Float64}(undef,ntStops,4)
        βdotSol = Matrix{Float64}(undef,ntStops,4)
        for t=1: ntStops
            rSol[t,:] = sol[14*i+1:14*i+3,t]
            vSol[t,:] = sol[14*i+8:14*i+10,t]
            βSol[t,:] = sol[14*i+4:14*i+7,t]
            βdotSol[t,:] = sol[14*i+11:14*i+14,t]
        end
        solFinal[i] = solSim(rSol,vSol,βSol,βdotSol)
    end
    return tSim, solFinal
end

function mainDynRK(t::Float64, x::Vector{Float64}, p::Tuple{Vararg{Joint}})
    dx = Vector{Float64}(undef,length(x))
    return mainDynODE!(dx,x,p,t)
end

function mainDynODE!(dX::Vector{Float64}, X::Vector{Float64}, j::Tuple{Vararg{Joint}}, t::Float64)
    # ODE function to be used as per DifferentialEquations convention
    # Create extForcesList storing extForces for each rigid body
    # Create ForceConstraints Array storing constraint forces acting on each rigid body
    global GravityInInertial
    # Update RigidBodies
    updateRigidBody!(j[1].RB1,X[1:14])
    for k=1:length(j)
        updateRigidBody!(j[k].RB2,X[14*k+1:14*(k+1)])
    end

    extFListCopy = externalFTotal(t,j) # Generates all the external forces explicitly specified, through joint actions or directly.

    unconstrF, constrF = Constraint(j, extFListCopy, GravityInInertial)
    # unconstrF is the matrix of the unconstrained Forces acting on the robot.
    # constrF is the matrix of the constraint Forces acting on the system developed using UK formulation

    dX[1:14] = zeros(14)
    for k=1:length(j)
        dX[14*k+1:14*(k+1)] = rbDynQuat(j[k].RB2, unconstrF[:,j[k].RB2.bodyID], constrF[:,j[k].RB2.bodyID])
    end
    # mainDyn!(dX, X, j, unconstrF, constrF)
    return dX
end

function externalFTotal(t::Float64, j::Tuple{Vararg{Joint}})
    # Returns externally applied forces in total, not including gravity. Includes Forces from joints and other explicitly specified external forces.
    extFList = extF(t,j)

    # Generate forces from actuated joints on each body
    ForceJoints = Matrix{Float64}(undef,6,2*length(j))
    for k=1:length(j)
        ForceJoints[:,2*k-1], ForceJoints[:,2*k] = genJointF(t,j[k])
    end

    # Add ForceJoints to extFList
    extFListCopy = deepcopy(extFList)
    for k=1:length(j)
        extFListCopy[j[k].RB1.bodyID].Forces =
        vcat(extFListCopy[j[k].RB1.bodyID].Forces, reshape(ForceJoints[1:3,2*k-1],(1,3)))
        extFListCopy[j[k].RB1.bodyID].Positions =
        vcat(extFListCopy[j[k].RB1.bodyID].Positions, reshape(j[k].pos1,(1,3)))
        extFListCopy[j[k].RB1.bodyID].Torques =
        vcat(extFListCopy[j[k].RB1.bodyID].Torques, reshape(ForceJoints[4:6,2*k-1],(1,3)))

        extFListCopy[j[k].RB2.bodyID].Forces =
        vcat(extFListCopy[j[k].RB2.bodyID].Forces, reshape(ForceJoints[1:3,2*k],(1,3)))
        extFListCopy[j[k].RB2.bodyID].Positions =
        vcat(extFListCopy[j[k].RB2.bodyID].Positions, reshape(j[k].pos2,(1,3)))
        extFListCopy[j[k].RB2.bodyID].Torques =
        vcat(extFListCopy[j[k].RB2.bodyID].Torques, reshape(ForceJoints[4:6,2*k],(1,3)))
    end
    return extFListCopy
end

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
module Discord

using Dates
using Distributed
using HTTP
using JSON
using OpenTrick
using Setfield
using TimeToLive

const DISCORD_JL_VERSION = v"0.1.0"
const API_VERSION = 6
const DISCORD_API = "https://discordapp.com/api"

const TTLDict = Dict{DataType, Union{Period, Nothing}}

function locked(f::Function, x::Threads.AbstractLock)
    lock(x)
    try f() finally unlock(x) end
end

function catchmsg(e::Exception)
    return sprint(showerror, e) * sprint(Base.show_backtrace, catch_backtrace())
end

include("types.jl")
include("events.jl")
include("state.jl")
include("limiter.jl")
include("client.jl")
include("gateway.jl")
include("rest.jl")
include("crud.jl")
include("commands.jl")
include("helpers.jl")
include("Defaults.jl")

end


text:
# This file contains the functions needed to define the gans
using LinearAlgebra

function crossEntTrue(x)
  return max(x, 0) - x + log(1 + exp(-abs(x)))
end

function crossEntFalse(x)
  return max(x, 0) + log(1 + exp(-abs(x)))
end

# creates a random orthogonal matrix for initialization
function rand_orth(m, n)
  @assert(m>=n)
  A = randn(m,n)
  return Matrix(qr(A).Q)
end

function ReLU(x)
  return max(x,zero(x))
end

# This function applies a dense net of specified depth, size, and nonlinearity
# with weights given by the first argument to the input data.
function denseNet(weights, data, inDim, outDim, depth, n_hidden, nonlin)
  offs = 1
  res = nonlin.(reshape(weights[offs : (offs + inDim * n_hidden - 1)],
                        (n_hidden, inDim)) * data)
  offs += inDim * n_hidden
  for k = 1 : depth
    res = nonlin.(reshape(weights[offs : (offs + n_hidden^2 - 1)],
                          (n_hidden, n_hidden)) * res)
    offs += n_hidden^2
  end
  return reshape(weights[offs : (offs + outDim * n_hidden - 1)], (outDim, n_hidden)) * res
end

# Creates a vector of suitable size to be the weights of a dense net and
# initializes each layer as a random orthonormal matrix.
function init_denseNet(inDim, outDim, depth, n_hidden, scale=0.8)
  if n_hidden >= inDim
    out = (scale * rand_orth(n_hidden, inDim))[:]
  else
    out = (scale * rand_orth(inDim, n_hidden)')[:]
  end

  for k = 1 : depth
    append!(out, scale * rand_orth(n_hidden, n_hidden)[:])
  end
  append!(out, (scale * rand_orth(n_hidden,outDim)')[:])
  return out
end

max_stars_repo_name: f-t-s/CGD


text:
"""
**ShotProfileEWEM**

*Shot Profile Elastic Wave Equation Migration and Demigration of 3D isotropic data.*

**IN**

* m : vector of filenames of image (comprised of [mpp;mps1;mps2])
* d : vector of filenames of data (comprised of [ux;uy;uz])
* adj : flag for adjoint (migration), or forward (demigration) (default=true)
* damping = 1000. : damping for deconvolution imaging condition
* vp = "vp" : seis file containing the p-wave velocity (should have same x and z dimensions as the desired image)
* vs = "vs" : seis file containing the s-wave velocity (should have same x and z dimensions as the desired image)
* angx = "angx" : seis file containing incidence angles in the x direction for each shot
* angy = "angy" : seis file containing incidence angles in the y direction for each shot
* wav = "wav" : seis file containing the source wavelet (in time domain)
* sz = 0. : source depth (Dev: read this from source wavelet file for variable source depth)
* gz = 0. : receiver depth (Dev: read this from data file for variable source depth (but then what to do in fwd op?))
* nangx = 1 : number of angle bins in x direction
* oangx = 0. : min angle in x direction (angle between source incidence angle and reflector normal in Degrees)
* dangx = 1. : angle increment in x direction
* nangy = 1 : number of angle bins in y direction
* oangy = 0. : min angle in y direction (angle between source incidence angle and reflector normal in Degrees)
* dangy = 1. : angle increment in y direction
* fmin = 0. : min frequency to process (Hz)
* fmax = 80. : max frequency to process (Hz)
* padt = 1 : pad factor for the time axis
* padx = 1 : pad factor for the spatial axes
* verbose = false : flag for error / debugging messages
* sx = [0.] : array of source X positions (meters)
* sy = [0.] : array of source Y positions (meters)

**OUT**

*Credits: AS, 2015*

"""

function ShotProfileEWEM(m::Array{AbstractString,1},d::Array{AbstractString,1},adj=true;damping=1000.,vp="vp.seis",vs="vs.seis",angx="angx.seis",angy="angy.seis",wav="wav.seis",sz=0.,gz=0.,nangx=1,oangx=0,dangx=1,nangy=1,oangy=0,dangy=1,fmin=0,fmax=80,padt=1,padx=1,verbose=false,sx=[0],sy=[0])

	nshot = length(sx)
	vel,h,extent = SeisRead(vp)
        min_imx = h[1].imx
        max_imx = h[end].imx
        nx = max_imx - min_imx + 1
        ox = h[1].mx
        dx = h[2].mx - ox
        min_imy = h[1].imy
        max_imy = h[end].imy
        ny = max_imy - min_imy + 1
        oy = h[1].my
        dy = ny > 1 ? h[nx+1].my - oy : dx
        nz = h[1].n1
        dz = h[1].d1
        oz = h[1].o1

	shot_list = Array(Shot3C,nshot)
	for ishot = 1 : nshot
		shot_list[ishot] = Shot3C(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
		shot_list[ishot].ux = join([d[1] "_shot_" Int(floor(sx[ishot])) "_" Int(floor(sy[ishot]))])
		shot_list[ishot].uy = join([d[2] "_shot_" Int(floor(sx[ishot])) "_" Int(floor(sy[ishot]))])
		shot_list[ishot].uz = join([d[3] "_shot_" Int(floor(sx[ishot])) "_" Int(floor(sy[ishot]))])
		shot_list[ishot].mpp = join([m[1] "_shot_" Int(floor(sx[ishot])) "_" Int(floor(sy[ishot]))])
		shot_list[ishot].mps1 = join([m[2] "_shot_" Int(floor(sx[ishot])) "_" Int(floor(sy[ishot]))])
		shot_list[ishot].mps2 = join([m[3] "_shot_" Int(floor(sx[ishot])) "_" Int(floor(sy[ishot]))])
		shot_list[ishot].angx = join(["angx_shot_" Int(floor(sx[ishot])) "_" Int(floor(sy[ishot]))])
		shot_list[ishot].angy = join(["angy_shot_" Int(floor(sx[ishot])) "_" Int(floor(sy[ishot]))])
		shot_list[ishot].vp = vp
		shot_list[ishot].vs = vs
		shot_list[ishot].wav = wav
		shot_list[ishot].damping = damping
		shot_list[ishot].sx = sx[ishot]
		shot_list[ishot].sy = sy[ishot]
		shot_list[ishot].sz = sz
		shot_list[ishot].gz = gz
		shot_list[ishot].fmin = fmin
		shot_list[ishot].fmax = fmax
		shot_list[ishot].padt = padt
		shot_list[ishot].padx = padx
		shot_list[ishot].adj = (adj == true) ? "y" : "n"
		shot_list[ishot].verbose = (verbose == true) ? "y" : "n"
	end

	if (nangx != 1 || nangy != 1)
		for ishot = 1 : nshot
			SeisWindow(angx,shot_list[ishot].angx,key=["sx","sy"],minval=[sx[ishot],sy[ishot]],maxval=[sx[ishot],sy[ishot]])
			SeisWindow(angy,shot_list[ishot].angy,key=["sx","sy"],minval=[sx[ishot],sy[ishot]],maxval=[sx[ishot],sy[ishot]])
		end
	end

	if (adj == true)
		for ishot = 1 : nshot
			SeisWindow(d[1],shot_list[ishot].ux,key=["sx","sy"],minval=[sx[ishot],sy[ishot]],maxval=[sx[ishot],sy[ishot]])
			SeisWindow(d[2],shot_list[ishot].uy,key=["sx","sy"],minval=[sx[ishot],sy[ishot]],maxval=[sx[ishot],sy[ishot]])
			SeisWindow(d[3],shot_list[ishot].uz,key=["sx","sy"],minval=[sx[ishot],sy[ishot]],maxval=[sx[ishot],sy[ishot]])
		end
		@sync @parallel for ishot = 1 : nshot
			a = shotewem(shot_list[ishot])
		end
		j = 1
		gather = zeros(Float32,nz,nangx*nangy)
		extent = Seismic.Extent(nz,nx,ny,nangx,nangy,
				oz,ox,oy,oangx,oangy,
				dz,dx,dy,dangx,dangy,
				"Depth","mx","my","angx","angy",
				"s","m","m","degrees","degrees",
				"")
		for imx = 1 : nx
			for imy = 1 : ny
				h = Header[]
				for iangx = 1 : nangx
					for iangy = 1 : nangy
						push!(h,Seismic.InitSeisHeader())
						h[(iangx-1)*nangy + iangy].tracenum = convert(typeof(h[1].tracenum),j + iangx*(nangy-1) + iangy)
						h[(iangx-1)*nangy + iangy].o1 = convert(typeof(h[1].o1),0)
						h[(iangx-1)*nangy + iangy].n1 = convert(typeof(h[1].n1),nz)
						h[(iangx-1)*nangy + iangy].d1 = convert(typeof(h[1].d1),dz)
						h[(iangx-1)*nangy + iangy].imx = convert(typeof(h[1].imx),imx-1 + min_imx)
						h[(iangx-1)*nangy + iangy].imy = convert(typeof(h[1].imy),imy-1 + min_imy)
						h[(iangx-1)*nangy + iangy].mx = convert(typeof(h[1].mx),dx*(imx-1 + min_imx))
						h[(iangx-1)*nangy + iangy].my = convert(typeof(h[1].my),dy*(imy-1 + min_imy))
						h[(iangx-1)*nangy + iangy].iang = convert(typeof(h[1].iang),iangx-1)
						h[(iangx-1)*nangy + iangy].ang = convert(typeof(h[1].ang),(iangx-1)*dangx + oangx)
						h[(iangx-1)*nangy + iangy].iaz = convert(typeof(h[1].iaz),iangy-1)
						h[(iangx-1)*nangy + iangy].az = convert(typeof(h[1].az),(iangy-1)*dangy + oangy)
					end
				end
				SeisWrite(m[1],gather,h,extent,itrace=j)
				SeisWrite(m[2],gather,h,extent,itrace=j)
				SeisWrite(m[3],gather,h,extent,itrace=j)
				j += nangx*nangy
			end
		end

		mpp_m = ParseDataName(m[1])
		mpp_h = ParseHeaderName(m[1])
		mps1_m = ParseDataName(m[2])
		mps1_h = ParseHeaderName(m[2])
		mps2_m = ParseDataName(m[3])
		mps2_h = ParseHeaderName(m[3])
		stream_mpp = open(mpp_m,"a+")
		stream_hpp = open(mpp_h,"a+")
		stream_mps1 = open(mps1_m,"a+")
		stream_hps1 = open(mps1_h,"a+")
		stream_mps2 = open(mps2_m,"a+")
		stream_hps2 = open(mps2_h,"a+")

		for ishot = 1 : nshot
			mpp_shot,hpp_shot,extent  = SeisRead(shot_list[ishot].mpp)
			mps1_shot,hps1_shot,extent  = SeisRead(shot_list[ishot].mps1)
			mps2_shot,hps2_shot,extent  = SeisRead(shot_list[ishot].mps2)
			if (nangx != 1 || nangy != 1)
				angx_shot,h_ang,extent = SeisRead(shot_list[ishot].angx)
				angy_shot,h_ang,extent = SeisRead(shot_list[ishot].angy)
			else
				angx_shot = 0.*mpp_shot
				angy_shot = 0.*mpp_shot
			end

			nx_shot = size(mpp_shot,2)
			for ix = 1 : nx_shot
				itrace = (hpp_shot[ix].imx)*ny*nangx*nangy + (hpp_shot[ix].imy)*nangx*nangy
				position_m = 4*nz*itrace
				seek(stream_mpp,position_m)
				m1 = read(stream_mpp,Float32,nz*nangx*nangy)
				m1 = reshape(m1,convert(Int64,nz),convert(Int64,nangx*nangy))
				seek(stream_mps1,position_m)
				m2 = read(stream_mps1,Float32,nz*nangx*nangy)
				m2 = reshape(m2,convert(Int64,nz),convert(Int64,nangx*nangy))
				seek(stream_mps2,position_m)
				m3 = read(stream_mps2,Float32,nz*nangx*nangy)
				m3 = reshape(m3,convert(Int64,nz),convert(Int64,nangx*nangy))
				for iz = 1 : nz
					# bilinear interpolation of angx and angy onto grid points
					angx = angx_shot[iz,ix]
					iangx = Int(floor((angx - oangx)/dangx)) + 1
					angx1 = (iangx-1)*dangx + oangx
					angx2 = angx1 + dangx
					angy = angy_shot[iz,ix]
					iangy = Int(floor((angy - oangy)/dangy)) + 1
					angy1 = (iangy-1)*dangy + oangy
					angy2 = angy1 + dangy
					w1 = (angx2-angx)*(angy2-angy)/((angx2-angx1)*(angy2-angy1))
					w2 = (angx-angx1)*(angy2-angy)/((angx2-angx1)*(angy2-angy1))
					w3 = (angx2-angx)*(angy-angy1)/((angx2-angx1)*(angy2-angy1))
					w4 = (angx-angx1)*(angy-angy1)/((angx2-angx1)*(angy2-angy1))
					if (iangx >= 1 && iangx <= nangx && iangy >= 1 && iangy <= nangy)
						m1[iz,(iangx-1)*nangy + iangy]   += w1*mpp_shot[iz,ix]
						if (iangx < nangx)
							m1[iz,(iangx)*nangy   + iangy]   += w2*mpp_shot[iz,ix]
						end
						if (nangy > 1)
							m1[iz,(iangx-1)*nangy + iangy+1] += w3*mpp_shot[iz,ix]
							if (iangx < nangx)
								m1[iz,(iangx)*nangy   + iangy+1] += w4*mpp_shot[iz,ix]
							end
						end
						m2[iz,(iangx-1)*nangy + iangy]   += w1*signf(angx)*signf(angy)*mps1_shot[iz,ix]
						if (iangx < nangx)
							m2[iz,(iangx)*nangy   + iangy]   += w2*signf(angx)*signf(angy)*mps1_shot[iz,ix]
						end
						if (nangy > 1)
							m2[iz,(iangx-1)*nangy + iangy+1] += w3*signf(angx)*signf(angy)*mps1_shot[iz,ix]
							if (iangx < nangx)
								m2[iz,(iangx)*nangy   + iangy+1] += w4*signf(angx)*signf(angy)*mps1_shot[iz,ix]
							end
						end
						m3[iz,(iangx-1)*nangy + iangy]   += w1*signf(angx)*signf(angy)*mps2_shot[iz,ix]
						if (iangx < nangx)
							m3[iz,(iangx)*nangy   + iangy]   += w2*signf(angx)*signf(angy)*mps2_shot[iz,ix]
						end
						if (nangy > 1)
							m3[iz,(iangx-1)*nangy + iangy+1] += w3*signf(angx)*signf(angy)*mps2_shot[iz,ix]
							if (iangx < nangx)
								m3[iz,(iangx)*nangy   + iangy+1] += w4*signf(angx)*signf(angy)*mps2_shot[iz,ix]
							end
						end
					end
				end
				seek(stream_mpp,position_m)
				write(stream_mpp,convert(Array{Float32,1},m1[:]))
				seek(stream_mps1,position_m)
				write(stream_mps1,convert(Array{Float32,1},m2[:]))
				seek(stream_mps2,position_m)
				write(stream_mps2,convert(Array{Float32,1},m3[:]))
			end
			SeisRemove(shot_list[ishot].ux)
			SeisRemove(shot_list[ishot].uy)
			SeisRemove(shot_list[ishot].uz)
			SeisRemove(shot_list[ishot].mpp)
			SeisRemove(shot_list[ishot].mps1)
			SeisRemove(shot_list[ishot].mps2)
			if (nangx != 1 || nangy != 1)
				SeisRemove(shot_list[ishot].angx)
				SeisRemove(shot_list[ishot].angy)
			end
		end
		close(stream_mpp)
		close(stream_hpp)
		close(stream_mps1)
		close(stream_hps1)
		close(stream_mps2)
		close(stream_hps2)
	else

		mpp_m = ParseDataName(m[1])
		mpp_h = ParseHeaderName(m[1])
		mps1_m = ParseDataName(m[2])
		mps1_h = ParseHeaderName(m[2])
		mps2_m = ParseDataName(m[3])
		mps2_h = ParseHeaderName(m[3])
		stream_mpp = open(mpp_m,"r")
		stream_hpp = open(mpp_h,"r")
		stream_mps1 = open(mps1_m,"r")
		stream_hps1 = open(mps1_h,"r")
		stream_mps2 = open(mps2_m,"r")
		stream_hps2 = open(mps2_h,"r")

		for ishot = 1 : nshot

			sx = shot_list[ishot].sx
			sy = shot_list[ishot].sy
			h_shot = Array(Header,nx*ny)
			for imx = 1 : nx
				for imy = 1 : ny
					h_shot[j] = Seismic.InitSeisHeader()
					h_shot[j].tracenum = convert(typeof(h_shot[1].tracenum),j)
					h_shot[j].o1 = convert(typeof(h_shot[1].o1),0)
					h_shot[j].n1 = convert(typeof(h_shot[1].n1),nz)
					h_shot[j].d1 = convert(typeof(h_shot[1].d1),dz)
					h_shot[j].imx = convert(typeof(h_shot[1].imx),imx-1 + min_imx)
					h_shot[j].imy = convert(typeof(h_shot[1].imy),imy-1 + min_imy)
					h_shot[j].mx = convert(typeof(h_shot[1].mx),dx*h_shot[j].imx)
					h_shot[j].my = convert(typeof(h_shot[1].my),dy*h_shot[j].imy)
					j += 1
				end
			end

			mpp_shot = zeros(nz,nx*ny)
			mps1_shot = zeros(nz,nx*ny)
			mps2_shot = zeros(nz,nx*ny)
			if (nangx != 1 || nangy != 1)
				angx_shot,h_ang,extent = SeisRead(shot_list[ishot].angx)
				angy_shot,h_ang,extent = SeisRead(shot_list[ishot].angy)
			else
				angx_shot = 0.*mpp_shot
				angy_shot = 0.*mpp_shot
			end
			for imx = 1 : nx
				for imy = 1 : ny
					ix = (imx - 1)*ny + imy
					h_shot[ix].sx = convert(typeof(h[1].sx),shot_list[ishot].sx)
					h_shot[ix].sy = convert(typeof(h[1].sy),shot_list[ishot].sy)
					itrace = (imx-1 + min_imx)*ny*nangx*nangy + (imy-1 + min_imy)*nangx*nangy
					position_m = 4*nz*itrace
					seek(stream_mpp,position_m)
					m1 = read(stream_mpp,Float32,nz*nangx*nangy)
					m1 = reshape(m1,convert(Int64,nz),convert(Int64,nangx*nangy))
					seek(stream_mps1,position_m)
					m2 = read(stream_mps1,Float32,nz*nangx*nangy)
					m2 = reshape(m2,convert(Int64,nz),convert(Int64,nangx*nangy))
					seek(stream_mps2,position_m)
					m3 = read(stream_mps2,Float32,nz*nangx*nangy)
					m3 = reshape(m3,convert(Int64,nz),convert(Int64,nangx*nangy))
					for iz = 1 : nz
						# bilinear interpolation of ang and az onto grid points
						angx = angx_shot[iz,ix]
						iangx = Int(floor((angx - oangx)/dangx)) + 1
						angx1 = (iangx-1)*dangx + oangx
						angx2 = angx1 + dangx
						angy = angy_shot[iz,ix]
						iangy = Int(floor((angy - oangy)/dangy)) + 1
						angy1 = (iangy-1)*dangy + oangy
						angy2 = angy1 + dangy
						w1 = (angx2-angx)*(angy2-angy)/((angx2-angx1)*(angy2-angy1))
						w2 = (angx-angx1)*(angy2-angy)/((angx2-angx1)*(angy2-angy1))
						w3 = (angx2-angx)*(angy-angy1)/((angx2-angx1)*(angy2-angy1))
						w4 = (angx-angx1)*(angy-angy1)/((angx2-angx1)*(angy2-angy1))
						if (iangx >= 1 && iangx <= nangx && iangy >= 1 && iangy <= nangy)
							mpp_shot[iz,ix] += w1*m1[iz,(iangx-1)*nangy + iangy]
							if (iangx < nangx)
								mpp_shot[iz,ix] += w2*m1[iz,(iangx)*nangy   + iangy]
							end
							if (nangy > 1)
								mpp_shot[iz,ix] += w3*m1[iz,(iangx-1)*nangy + iangy+1]
								if (iangx < nangx)
									mpp_shot[iz,ix] += w4*m1[iz,(iangx)*nangy   + iangy+1]
								end
							end
							mps1_shot[iz,ix] += w1*signf(angx)*signf(angy)*m2[iz,(iangx-1)*nangy + iangy]
							if (iangx < nangx)
								mps1_shot[iz,ix] += w2*signf(angx)*signf(angy)*m2[iz,(iangx)*nangy   + iangy]
							end
							if (nangy > 1)
								mps1_shot[iz,ix] += w3*signf(angx)*signf(angy)*m2[iz,(iangx-1)*nangy + iangy+1]
								if (iangx < nangx)
									mps1_shot[iz,ix] += w4*signf(angx)*signf(angy)*m2[iz,(iangx)*nangy   + iangy+1]
								end
							end
							mps2_shot[iz,ix] += w1*signf(angx)*signf(angy)*m3[iz,(iangx-1)*nangy + iangy]
							if (iangx < nangx)
								mps2_shot[iz,ix] += w2*signf(angx)*signf(angy)*m3[iz,(iangx)*nangy   + iangy]
							end
							if (nangy > 1)
								mps2_shot[iz,ix] += w3*signf(angx)*signf(angy)*m3[iz,(iangx-1)*nangy + iangy+1]
								if (iangx < nangx)
									mps2_shot[iz,ix] += w4*signf(angx)*signf(angy)*m3[iz,(iangx)*nangy   + iangy+1]
								end
							end
						end
					end
				end
			end
			SeisWrite(shot_list[ishot].mpp,mpp_shot,h_shot)
			SeisWrite(shot_list[ishot].mps1,mps1_shot,h_shot)
			SeisWrite(shot_list[ishot].mps2,mps2_shot,h_shot)
		end
		close(stream_mpp)
		close(stream_hpp)
		close(stream_mps1)
		close(stream_hps1)
		close(stream_mps2)
		close(stream_hps2)
		@sync @parallel for ishot = 1 : nshot
			a = shotewem(shot_list[ishot])
		end
		j = 1
		for ishot = 1 : nshot
			ux_shot,h_shot = SeisRead(shot_list[ishot].ux)
			SeisWrite(d[1],ux_shot,h_shot,itrace=j)
			uy_shot,h_shot = SeisRead(shot_list[ishot].uy)
			SeisWrite(d[2],uy_shot,h_shot,itrace=j)
			uz_shot,h_shot = SeisRead(shot_list[ishot].uz)
			SeisWrite(d[3],uz_shot,h_shot,itrace=j)
			j += size(ux_shot,2)
			SeisRemove(shot_list[ishot].ux)
			SeisRemove(shot_list[ishot].uy)
			SeisRemove(shot_list[ishot].uz)
			SeisRemove(shot_list[ishot].mpp)
			SeisRemove(shot_list[ishot].mps1)
			SeisRemove(shot_list[ishot].mps2)
			if (nangx != 1 || nangy != 1)
				SeisRemove(shot_list[ishot].angx)
				SeisRemove(shot_list[ishot].angy)
			end
		end
	end
end

type Shot3C
	ux
	uy
	uz
	mpp
	mps1
	mps2
	angx
	angy
	vp
	vs
	wav
	damping
	sx
	sy
	sz
	gz
	fmin
	fmax
	padt
	padx
	adj
	verbose
end

function shotewem(shot)

	@compat if (Libdl.find_library(["shotewem"]) == "")
		error("Couldn't find shared library shotewem.so, make sure it is installed and check LD_LIBRARY_PATH environmental variable.")
	end
	argv = ["0",
	join(["adj=",shot.adj]),
	join(["ux=",shot.ux]), join(["uy=",shot.uy]), join(["uz=",shot.uz]),
	join(["mpp=",shot.mpp]), join(["mps1=",shot.mps1]), join(["mps2=",shot.mps2]),
	join(["vp=",shot.vp]), join(["vs=",shot.vs]), join(["wav=",shot.wav]),
	join(["damping=",shot.damping]), join(["sx=",shot.sx]), join(["sy=",shot.sy]),  join(["sz=",shot.sz]),  join(["gz=",shot.gz]),
	join(["fmin=",shot.fmin]),  join(["fmax=",shot.fmax]),
	join(["padt=",shot.padt]),  join(["padx=",shot.padx]),
	join(["verbose=",shot.verbose]) ]
	@compat a = ccall((:main, "shotewem"), Int32, (Int32, Ptr{Ptr{UInt8}}), length(argv), argv)
	return(a)

end

function signf(a)

	b = a < 0.0 ? -1.0 : 1.0
	return(b)

end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ 2d49dec2-688f-11eb-1528-216d9f93be24
using Pkg

# ╔═╡ 5cc0237e-688f-11eb-2b9a-255add42261f
using MechanicalSketch

# ╔═╡ 501c0e00-688f-11eb-2ef6-010b120c3e6d
Pkg.status()

# ╔═╡ 7e85c4c0-688f-11eb-3b19-25c83b4a3df3
import MechanicalSketch.placeimage

# ╔═╡ 8b8abef0-688f-11eb-2d65-4303ac879bf1
import MechanicalSketch.@imagematrix

# ╔═╡ a61295e0-688f-11eb-2382-836970fb11d9

# ╔═╡ Cell order:
# ╠═2d49dec2-688f-11eb-1528-216d9f93be24
# ╠═501c0e00-688f-11eb-2ef6-010b120c3e6d
# ╠═5cc0237e-688f-11eb-2b9a-255add42261f
# ╠═7e85c4c0-688f-11eb-3b19-25c83b4a3df3
# ╠═8b8abef0-688f-11eb-2d65-4303ac879bf1
# ╠═a61295e0-688f-11eb-2382-836970fb11d9


text:
using DynamicalSystemsBase
using Statistics, Test

ds = Systems.lorenz()
ALG = SimpleATsit5()

println("\nTesting internal norms...")

# %%
@testset begin "internal norm: tangent"
i1 = tangent_integrator(ds, 3; diffeq = (alg = ALG,))
i2 = tangent_integrator(ds, 3; diffeq = (alg = ALG, internalnorm = DynamicalSystemsBase._standardnorm))

step!(i1)
step!(i2)
N = 1000
dts = zeros(2, N)
for i in 1:N
    step!(i1); step!(i2)
    dts[:, i] .= (i1.dt, i2.dt)
end

dt1 = mean(dts[1, :])
dt2 = mean(dts[2, :])

@test dt1 > dt2
end
# %%
@testset begin "internal norm: parallel"
s = [get_state(ds), get_state(ds) .+ rand(3)]
i1 = parallel_integrator(ds, deepcopy(s);
diffeq = (alg = ALG, internalnorm = DynamicalSystemsBase._standardnorm))
i2 = parallel_integrator(ds, deepcopy(s); diffeq = (alg = ALG,))

step!(i1)
step!(i2)

N = 2000
dts = zeros(2, N)
for i in 1:N
    step!(i1); step!(i2)
    dts[:, i] .= (i1.dt, i2.dt)
end

dt1 = mean(dts[1, :])
dt2 = mean(dts[2, :])

@test dt2 > dt1
end


text:
using NeuralGRN, LightGraphs, Plots
using Test

@testset "Problem Construction Test" begin
    g = complete_graph(3)
    @test_throws AssertionError SteadyStateGRNProblem{Int64, Float32}(g, [1, 2])
    @test_throws AssertionError TemporalGRNProblem{Int64, Float32}(g,  [1 2 3; 4 5 6])

    @test TemporalGRNProblem{Int64, Float32}(g, [1 2; 3 4; 5 6]).names == string.(1:3)
end

@testset "Plotting Tests" begin
    g = complete_graph(3)

    @test typeof(plot_problem_graph(TemporalGRNProblem{Int64, Float32}(g,
        [1 2; 3 4; 5 6]))) == Plots.Plot{Plots.GRBackend}
end


text:
@testset "LU solver" begin
    n = 20
    m = 10
    A = rand(n, n)
    X = rand(n, m)
    B = rand(n, m)
    x = rand(n)
    b = rand(n)

    sol = ContactImplicitMPC.lu_solver(A)
    ContactImplicitMPC.linear_solve!(sol, X, A, B)
    @test norm(A * X - B, Inf) < 1e-10

    sol = ContactImplicitMPC.lu_solver(A)
    ContactImplicitMPC.linear_solve!(sol, x, A, b)
    @test norm(A * x - b, Inf) < 1e-10
end

# n = 20
# m = 10
# A = rand(n, n)
# X = rand(n, m)
# B = rand(n, m)
# x = rand(n)
# b = rand(n)
#
# sol = ContactImplicitMPC.lu_solver(A)
# @benchmark ContactImplicitMPC.linear_solve!($sol, $X, $A, $B)
# @test norm(A * X - B, Inf) < 1e-10
#
# sol = ContactImplicitMPC.lu_solver(A)
# @benchmark ContactImplicitMPC.linear_solve!($sol, $x, $A, $b)
# @test norm(A * x - b, Inf) < 1e-10
# @benchmark factorize!($sol, $A)
# @benchmark ldiv!($x, $(sol.lu), $b)


text:
#!/usr/bin/julia

using ProgressBars
using Profile

const OP_MAP = Dict(
    "add" => "+",
    "mul" => "*",
    "div" => "div",
    "mod" => "rem",
    "eql" => "==",
)

function compile(instlist)
    linelist = String[]
    push!(linelist, "function (inputs::AbstractVector{T}) where T <: Integer")

    # push!(linelist, "rev_inputs = Int[]")
    # push!(linelist, "bigend, digit = divrem(n, 10)")
    # push!(linelist, "while bigend > 0")
    # push!(linelist, "append!(rev_inputs, digit)")
    # push!(linelist, "bigend, digit = divrem(bigend, 10)")
    # push!(linelist, "end")
    # push!(linelist, "inputs = reverse(rev_inputs)")

    # push!(linelist, "regs = zeros(T, 4)")
    push!(linelist, "w = 0")
    push!(linelist, "x = 0")
    push!(linelist, "y = 0")
    push!(linelist, "z = 0")
    # push!(linelist, "next = iterate(inputs)")
    push!(linelist, "input_i = 1")
    push!(linelist, "input_max = length(inputs)")
    for inst in instlist
        tokens = split(inst, " ")
        if length(tokens) == 2
            # push!(linelist, "if input_i <= input_max")
            # push!(linelist, "(item, state) = next")
            push!(linelist, "if inputs[input_i] <= 0 || inputs[input_i] > 9")
            push!(linelist, "return 0, 0, 0, 0")
            push!(linelist, "end")
            push!(linelist, "$(tokens[2]) = inputs[input_i]")
            push!(linelist, "input_i += 1")
            # push!(linelist, "next = iterate(inputs, state)")
            # push!(linelist, "end")
        else
            op, a, b = tokens
            opfunc = OP_MAP[op]
            push!(linelist, "$(a) = $(opfunc)($(a), $(b))")
        end
    end
    push!(linelist, "return w, x, y, z")
    push!(linelist, "end")
    expr_str = join(linelist, "\n")
    return Meta.parse(expr_str)
end

function run_func(monad_func, num_iter=9^14)
    works = 0
    digits = repeat([1], 14)
    # for _ in 1:num_iter #ProgressBar(1:num_iter)
    for _ in ProgressBar(1:cld(num_iter, 9^8))
        finished = true
        for _ in 1:9^8
            finished = true
            # @info [reverse(digits)...]
            # @info digits
            w, x, y, z = monad_func(digits)
            if z == 0
                @info digits
                works = sum(d * 10^(14-i) for (i, d) in enumerate(digits))
                finished = true
                break
            end
            for di in 14:-1:1
                next_d = digits[di] + 1
                if next_d <= 9
                    digits[di] = next_d
                    finished = false
                    break
                end
                digits[di] = 1
            end
            if finished
                break
            end
            # regs = monad_func([digits...])
            # @info regs
        end
        if finished
            break
        end
    end
    return works
end

function make_func()
    lines = open("input/24.txt", "r") do io
        map(strip, readlines(io))
    end
    monad_expr = compile(lines)
    monad_func = eval(monad_expr)
    return monad_expr, monad_func

    # @time count_trees(maze, 1, 3)
    # @timev count_trees(maze, 1, 3)
    # @info "Found answer" count_trees(maze, 1, 3)
    # @timev prod(count_trees(maze, l...) for l in check_list)
    # @info "Found answer" prod(count_trees(maze, l...) for l in check_list)
end

if abspath(PROGRAM_FILE) == @__FILE__
    _, monad_func = make_func()
    @time run_func(monad_func, 10^6)
    # Profile.clear_malloc_data()
    @timev run_func(monad_func, 9^14)
    # Profile.print()
end


text:
function InterSatellitesRange(arg0::ObservableSatellite, arg1::ObservableSatellite, arg2::jboolean, arg3::AbsoluteDate, arg4::jdouble, arg5::jdouble, arg6::jdouble)
    return InterSatellitesRange((ObservableSatellite, ObservableSatellite, jboolean, AbsoluteDate, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function is_two_way(obj::InterSatellitesRange)
    return jcall(obj, "isTwoWay", jboolean, ())
end


text:
type SharedMemoryInputs <: Inputs
    ratings_file::FileSpec
    nusers::Int
    nitems::Int
    R::Nullable{InputRatings}
    RT::Nullable{InputRatings}
    item_idmap::Nullable{InputIdMap}
    user_idmap::Nullable{InputIdMap}

    function SharedMemoryInputs(file::FileSpec)
        new(file, 0, 0, nothing, nothing, nothing, nothing)
    end
end

function clear(inp::SharedMemoryInputs)
    inp.R = nothing
    inp.RT = nothing
    inp.item_idmap = nothing
    inp.user_idmap = nothing
end

function localize!(inp::SharedMemoryInputs)
    if !isnull(inp.R)
        R = get(inp.R)
        isa(R, SharedRatingMatrix) && (inp.R = copy(R))
    end

    if !isnull(inp.RT)
        RT = get(inp.RT)
        isa(RT, SharedRatingMatrix) && (inp.RT = copy(RT))
    end

    if !isnull(inp.item_idmap)
        item_idmap = get(inp.item_idmap)
        isa(item_idmap, SharedVector) && (inp.item_idmap = copy(item_idmap))
    end

    if !isnull(inp.user_idmap)
        user_idmap = get(inp.user_idmap)
        isa(user_idmap, SharedVector) && (inp.user_idmap = copy(user_idmap))
    end
    nothing
end

function share!(inp::SharedMemoryInputs)
    if !isnull(inp.R)
        R = get(inp.R)
        isa(R, SharedRatingMatrix) || (inp.R = share(R))
    end

    if !isnull(inp.RT)
        RT = get(inp.RT)
        isa(RT, SharedRatingMatrix) || (inp.RT = share(RT))
    end

    if !isnull(inp.item_idmap)
        item_idmap = get(inp.item_idmap)
        isa(item_idmap, SharedVector) || (inp.item_idmap = share(item_idmap))
    end

    if !isnull(inp.user_idmap)
        user_idmap = get(inp.user_idmap)
        isa(user_idmap, SharedVector) || (inp.user_idmap = share(user_idmap))
    end

    nothing
end

# Note:
# Filtering out causes the item and user ids to change.
# We need to keep a mapping to be able to match in the recommend step.
function filter_empty(R::RatingMatrix; only_items=Int64[])
    if !isempty(only_items)
        max_num_items = maximum(only_items)
        if size(R, 2) < max_num_items
            RI, RJ, RV = findnz(R)
            R = sparse(RI, RJ, RV, size(R, 1), max_num_items)
        end
        R = R'
        R = R[only_items, :]
        R = R'
        non_empty_items = only_items
    end

    U = sum(R, 2)
    non_empty_users = find(U)
    R = R[non_empty_users, :]

    if isempty(only_items)
        R = R'
        P = sum(R, 2)
        non_empty_items = find(P)
        R = R[non_empty_items, :]
        R = R'
    end

    R, non_empty_items, non_empty_users
end

function ensure_loaded(inp::SharedMemoryInputs; only_items=Int64[])
    if isnull(inp.R)
        @logmsg("loading inputs...")
        t1 = time()
        A = read_input(inp.ratings_file)

        if isa(A, SparseMatrixCSC)
            R = convert(SparseMatrixCSC{Float64,Int64}, A)
        else
            # separate the columns and make them of appropriate types
            users   = convert(Vector{Int64},   A[:,1])
            items   = convert(Vector{Int64},   A[:,2])
            ratings = convert(Vector{Float64}, A[:,3])

            # create a sparse matrix
            R = sparse(users, items, ratings)
        end

        R, item_idmap, user_idmap = filter_empty(R; only_items=only_items)
        inp.R = R
        inp.nusers, inp.nitems = size(R)
        inp.item_idmap = (extrema(item_idmap) == (1,length(item_idmap))) ? nothing : item_idmap
        inp.user_idmap = (extrema(user_idmap) == (1,length(user_idmap))) ? nothing : user_idmap
        inp.RT = R'
        t2 = time()
        isnull(inp.item_idmap) && @logmsg("no need to map item_ids")
        isnull(inp.user_idmap) && @logmsg("no need to map user_ids")
        @logmsg("time to load inputs: $(t2-t1) secs")
    end
    nothing
end

item_idmap(inp::SharedMemoryInputs) = isnull(inp.item_idmap) ? Int64[] : get(inp.item_idmap)
user_idmap(inp::SharedMemoryInputs) = isnull(inp.user_idmap) ? Int64[] : get(inp.user_idmap)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Markdown

import Base: show, ==
import Core: @doc_str

include(joinpath("parse", "config.jl"))
include(joinpath("parse", "util.jl"))
include(joinpath("parse", "parse.jl"))

include(joinpath("Common", "Common.jl"))
include(joinpath("GitHub", "GitHub.jl"))
include(joinpath("IPython", "IPython.jl"))
include(joinpath("Julia", "Julia.jl"))

include(joinpath("render", "plain.jl"))
include(joinpath("render", "html.jl"))
include(joinpath("render", "latex.jl"))
include(joinpath("render", "rst.jl"))

include(joinpath("render", "terminal", "render.jl"))

export readme, license, @md_str, @doc_str

parse(markdown::AbstractString; flavor = julia) = parse(IOBuffer(markdown), flavor = flavor)
parse_file(file::AbstractString; flavor = julia) = parse(readstring(file), flavor = flavor)

readme(pkg::AbstractString; flavor = github) = parse_file(Pkg.dir(pkg, "README.md"), flavor = flavor)
readme(pkg::Module; flavor = github) = readme(string(pkg), flavor = flavor)

license(pkg::AbstractString; flavor = github) = parse_file(Pkg.dir(pkg, "LICENSE.md"), flavor = flavor)
license(pkg::Module; flavor = github) = license(string(pkg), flavor = flavor)

function mdexpr(s, flavor = :julia)
    md = parse(s, flavor = Symbol(flavor))
    esc(toexpr(md))
end

function docexpr(s, flavor = :julia)
    quote
        let md = $(mdexpr(s, flavor))
            md.meta[:path] = @__FILE__
            md.meta[:module] = current_module()
            md
        end
    end
end

macro md_str(s, t...)
    mdexpr(s, t...)
end

doc_str(md, file, mod) = (md.meta[:path] = file; md.meta[:module] = mod; md)
doc_str(md::AbstractString, file, mod) = doc_str(parse(md), file, mod)

macro doc_str(s::AbstractString, t...)
    :($(doc_str)($(mdexpr(s, t...)), $(Base).@__FILE__, $(current_module)()))
end

function Base.display(d::Base.REPL.REPLDisplay, md::Vector{MD})
    for md in md
        display(d, md)
    end
end

end

max_stars_repo_name: odow/julia

max_issues_repo_name: odow/julia

max_forks_repo_name: odow/julia


text:
using NOVAS
using Documenter

DocMeta.setdocmeta!(NOVAS, :DocTestSetup, :(using NOVAS); recursive=true)

makedocs(; modules=[NOVAS], authors="Kiran Shila <me@kiranshila.com> and contributors",
         repo="https://github.com/kiranshila/NOVAS.jl/blob/{commit}{path}#{line}",
         sitename="NOVAS.jl",
         format=Documenter.HTML(; prettyurls=get(ENV, "CI", "false") == "true",
                                canonical="https://kiranshila.github.io/NOVAS.jl",
                                assets=[asset("https://cdn.plot.ly/plotly-2.6.3.min.js")]),
         pages=["Home" => "index.md", "Examples" => "examples.md",
                "Benchmarks" => "benchmarks.md", "API" => "api.md"])

deploydocs(; repo="github.com/kiranshila/NOVAS.jl", devbranch="main")