text:
"""
    AlgoDeforLinearModule

Module for algorithms used in linear deformation models.
"""
module AlgoDeforLinearModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.AlgoBaseModule: dcheck!
import Arpack: eigs
import SparseArrays: spzeros
import LinearAlgebra: mul!
my_A_mul_B!(C, A, B) = mul!(C, A, B)
import FinEtools.FieldModule: Field, ndofs, setebc!, numberdofs!, applyebc!, scattersysvec!
import FinEtools.NodalFieldModule: NodalField, nnodes
import FinEtools.FEMMBaseModule: associategeometry!, distribloads, fieldfromintegpoints, elemfieldfromintegpoints
import FinEtools.FEMMDeforLinearBaseModule: stiffness, mass, nzebcloadsstiffness, thermalstrainloads, inspectintegpoints
import FinEtools.FEMMDeforLinearMSModule: stiffness, mass, nzebcloadsstiffness, thermalstrainloads, inspectintegpoints
import FinEtools.DeforModelRedModule: stresscomponentmap
import FinEtools.ForceIntensityModule: ForceIntensity
import FinEtools.MeshModificationModule: meshboundary
import FinEtools.MeshExportModule: vtkexportmesh
import LinearAlgebra: eigen, qr, dot, cholesky

"""
    AlgoDeforLinearModule.linearstatics(modeldata::FDataDict)

Algorithm for static linear deformation (stress) analysis.

`modeldata` = dictionary with values for keys

* "fens"  = finite element node set
* "regions"  = array of region dictionaries
* "essential_bcs" = array of essential boundary condition dictionaries
* "traction_bcs" = array of traction boundary condition dictionaries
* "temperature_change" = dictionary of data for temperature change

For each region (connected piece of the domain made of a particular material),
mandatory, the  region dictionary  contains values for keys:
* "femm" = finite element model machine (mandatory);

For essential boundary conditions (optional) each dictionary
would hold
  + "displacement" = fixed (prescribed) displacement (scalar),  or
            a function with signature
                function w = f(x)
            If not given, zero displacement assumed.
  + "component" = which component is prescribed  (1, 2, 3)?
  + "node_list" = list of nodes on the boundary to which the condition applies
            (mandatory)

For traction boundary conditions (optional) each dictionary
would hold
  + "femm" = finite element model machine (mandatory);
  + "traction_vector" = traction vector,  either  a constant or  a function
        Positive  when outgoing.

Output:
modeldata= the dictionary on input is augmented with
* "geom" = the nodal field that is the geometry
* "u" = the nodal field that is the computed displacement
* "dT" = the nodal field that is the temperature increment
"""
function linearstatics(modeldata::FDataDict)

    # For traction boundary conditions (optional):
    # model_data.boundary_conditions.traction = cell array of struct,
    #           each piece of surface with traction boundary condition gets one
    #           element of the array with a struct with the attributes
    #     traction=traction (vector), supply a zero for component in which
    #           the boundary condition is inactive
    #     fes = finite element set on the boundary to which
    #                       the condition applies
    #     integration_rule= integration rule
    #
    # For body loads (optional):
    # model_data.body_load = cell array of struct,
    #          each piece of the domain can have each its own body load
    #     force  = force density vector
    #     fes = finite element set to which the load applies
    #     integration_rule= integration rule
    #
    # For multi point constraints (MPC) (optional):
    # model_data.mpc= cell array of structs, each for one MPC.
    #      node_list = list of node numbers involved in the MPC,
    #      dof_list= numbers of degrees of freedom for the nodes above,
    #      umultipliers=multipliers for the nodes above,
    #      penfact=the penalty factor to multiply  the constraint matrix,
    #          The MPC looks like this: sum_i m_i u_{dof(i),node(i)} =0
    #          where m_i is the multiplier.

    # Lists of recognized keys for the data dictionaries:
    modeldata_recognized_keys = ["fens", "regions",
    "essential_bcs",  "traction_bcs", "temperature_change",
    "factorize"]
    essential_bcs_recognized_keys = ["displacement", "node_list", "component"]
    traction_bcs_recognized_keys = ["femm", "traction_vector"]
    regions_recognized_keys = ["femm", "body_load"]
    temperature_change_recognized_keys = ["temperature"]

    # Extract the nodes
    fens=get(()->error("Must get fens!"), modeldata, "fens")

    # Construct the geometry field
    geom = NodalField(fens.xyz)

    # Construct the displacement field
    u = NodalField(zeros(nnodes(geom), ndofs(geom)))

    # Construct the temperature field
    temp = NodalField(zeros(nnodes(geom), 1))

    modeldata["timing"] = FDataDict()

    tstart = time()
    # Apply the essential boundary conditions on the displacement field
    essential_bcs = get(modeldata, "essential_bcs", nothing);
    if (essential_bcs != nothing)
        for j = 1:length(essential_bcs)
            ebc = essential_bcs[j]
            dcheck!(ebc, essential_bcs_recognized_keys)
            fenids = get(()->error("Must get node list!"), ebc, "node_list");
            displacement = get(ebc, "displacement", nothing);
            u_fixed = zeros(FFlt, length(fenids)); # default is  zero displacement
            if (displacement != nothing) # if it is nonzero,
                if (typeof(displacement) <: Function) # it could be a function
                    for k = 1:length(fenids)
                        u_fixed[k] = displacement(geom.values[fenids[k],:])[1];
                    end
                else # or it could be a constant
                    fill!(u_fixed, displacement);
                end
            end
            component = get(ebc, "component", 0); # which component?
            setebc!(u, fenids[:], true, component, u_fixed);
        end
        applyebc!(u);
    end

    # Number the equations
    numberdofs!(u)           #,Renumbering_options); # NOT DONE <<<<<<<<<<<<<<<<<
    modeldata["timing"]["essential_bcs"] = time() - tstart

    tstart = time()
    # Initialize the heat loads vector
    F = zeros(FFlt,u.nfreedofs);
    # Construct the system stiffness matrix
    K = spzeros(u.nfreedofs,u.nfreedofs); # (all zeros, for the moment)
    regions = get(()->error("Must get region list!"), modeldata, "regions")
    for i = 1:length(regions)
        region = regions[i]
        dcheck!(region, regions_recognized_keys)
        femm = region["femm"];
        # Give the  FEMM a chance  to precompute  geometry-related quantities
        femm = associategeometry!(femm, geom);
        # Add up all the conductivity matrices for all the regions
        K = K + stiffness(femm, geom, u);
        # Loads due to the essential boundary conditions on the displacement field
        essential_bcs = get(modeldata, "essential_bcs", nothing);
        if (essential_bcs != nothing) # there was at least one EBC applied
            F = F + nzebcloadsstiffness(femm, geom, u);
        end
    end
    modeldata["timing"]["stiffness"] = time() - tstart

    # # Process the body load
    # body_load = get(modeldata, "body_load", nothing);
    # if (body_load  !=nothing)
    #     for j=1:length(model_data.body_load)
    #         body_load =model_data.body_load{j};
    #         femm = femm_deformation_linear (struct ('material',[],...
    #             'fes',body_load.fes,...
    #             'integration_rule',body_load.integration_rule));
    #         fi= force_intensity(struct('magn',body_load.force));
    #         F = F + distrib_loads(femm, sysvec_assembler, geom, u, fi, 3);
    #     end
    #     clear body_load fi  femm
    # end

    tstart = time()
    # Process the traction boundary condition
    traction_bcs = get(modeldata, "traction_bcs", nothing);
    if (traction_bcs != nothing)
        for j=1:length(traction_bcs)
            tractionbc = traction_bcs[j]
            dcheck!(tractionbc, traction_bcs_recognized_keys)
            traction_vector = tractionbc["traction_vector"];
            if (typeof(traction_vector) <: Function)
                fi = ForceIntensity(FFlt, ndofs(geom), traction_vector);
            else
                fi = ForceIntensity(traction_vector);
            end
            femm = tractionbc["femm"]
            F = F + distribloads(femm, geom, u, fi, 2);
        end
    end
    modeldata["timing"]["traction_bcs"] = time() - tstart

    tstart = time()
    # Process the thermal strain  loading
    temperature_change = get(modeldata, "temperature_change", nothing);
    if (temperature_change != nothing)
        dcheck!(temperature_change, temperature_change_recognized_keys)
        # zero temperature change is a reasonable default
        temp = NodalField(zeros(size(fens.xyz,1),1))
        temperature = get(temperature_change, "temperature", nothing);
        if (temperature != nothing) # if it is nonzero,
            if (typeof(temperature) <: Function) # it could be a function
                for k = 1:count(fens)
                    temp.values[k] = temperature(geom.values[k,:])[1];
                end
            else # or it could be a constant
                fill!(temp.values, temperature);
            end
        end
        for i = 1:length(regions)
            region = regions[i]
            femm = region["femm"];
            F = F + thermalstrainloads(femm, geom, u, temp)
        end
    end
    modeldata["timing"]["temperature_change"] = time() - tstart

    # # Process the nodal force boundary condition
    # if (isfield(model_data.boundary_conditions, 'nodal_force' ))
    #     for j=1:length(model_data.boundary_conditions.nodal_force)
    #         nodal_force =model_data.boundary_conditions.nodal_force{j};
    #         femm = femm_deformation_linear (struct ('material',[],...
    #             'fes',fe_set_P1(struct('conn',reshape(nodal_force.node_list,[],1))),...
    #             'integration_rule',point_rule));
    #         fi= force_intensity(struct('magn',nodal_force.force));
    #         F = F + distrib_loads(femm, sysvec_assembler, geom, u, fi, 0);
    #     end
    #     clear nodal_force fi femm
    # end

    # # Apply multi point constraints
    # if isfield(model_data,'mpc')
    #     for i=1:length(model_data.mpc)
    #         mpc =model_data.mpc{i};
    #         dofnums=0*mpc.umultipliers;# Construct an array of the degree of freedom numbers
    #         for kx=1:length(mpc.node_list)
    #             dofnums(kx)=u.dofnums(mpc.node_list(kx),mpc.dof_list(kx));
    #         end
    #         # Now call the utility function to calculate the constraint matrix
    #         [Kmpc,Fmpc]=apply_penalty_mpc(u.nfreedofs,dofnums,mpc.umultipliers,0.0,mpc.penfact);
    #         K = K + Kmpc;
    #         F = F + Fmpc;
    #     end
    #     clear Kmpc Fmpc
    # end

    tstart = time()
    # Solve the system of linear algebraic equations
    K = cholesky(K);
    U = K\F;
    scattersysvec!(u, U[:])
    modeldata["timing"]["solution"] = time() - tstart

    # Update the model data
    setindex!(modeldata, geom, "geom");
    setindex!(modeldata, u, "u");
    setindex!(modeldata, temp, "temp");
    setindex!(modeldata, temp, "dT");
    setindex!(modeldata, dot(F,U)/2, "work");
    return modeldata            # ... And return the updated model data
end

"""
    AlgoDeforLinearModule.exportdeformation(modeldata::FDataDict)

Algorithm for exporting of the deformation for visualization in Paraview.

`modeldata` = dictionary with values for keys

* "fens"  = finite element node set
* "regions"  = array of region dictionaries
* "geom" = geometry field
* "u" = displacement field, or
* "us" = array of  tuples (name, displacement field)
* "postprocessing" = dictionary  with values for keys
  + "boundary_only" = should only the boundary of the  regions be rendered?
                      Default is render the interior.
  + "file" = name of the  postprocessing file

For each region (connected piece of the domain made of a particular material),
mandatory, the  region dictionary  contains values for keys:
* "femm" = finite element mmodel machine (mandatory);

Output: modeldata updated with
* modeldata["postprocessing"]["exported"] = array of data dictionaries, one for
        each exported file. The data is stored with the keys:
    - "file" - names of exported file
    - "field" - nodal or elemental field
"""
function exportdeformation(modeldata::FDataDict)
    modeldata_recognized_keys = ["fens", "regions", "geom", "u", "postprocessing"]
    postprocessing_recognized_keys = ["boundary_only", "file"]
    # Defaults
    boundary_only = false;
    ffile = "deformation"
    dcheck!(modeldata, modeldata_recognized_keys)

    # Let's have a look at what's been specified
    postprocessing = get(modeldata, "postprocessing", nothing);
    if (postprocessing != nothing)
        dcheck!(postprocessing, postprocessing_recognized_keys)
        boundary_only =  get(postprocessing, "boundary_only", boundary_only);
        ffile =  get(postprocessing, "file", ffile);
    end

    fens = get(()->error("Must get fens!"), modeldata, "fens")
    geom = get(()->error("Must get geometry field!"), modeldata, "geom");
    u = get(modeldata, "u", nothing);
    us = get(modeldata, "us", nothing);
    if us == nothing
        us = [("u", u)]
    end

    # Export one file for each region
    modeldata["postprocessing"]["exported"] = Array{FDataDict, 1}()
    regions = get(()->error("Must get region!"), modeldata, "regions")
    for i = 1:length(regions)
        region = regions[i]
        femm = region["femm"]
        rfile = ffile * "$i" * ".vtk";
        vectors = Tuple{String, FFltMat}[]
        for ixxxx = 1:length(us)
            push!(vectors, (us[ixxxx][1], us[ixxxx][2].values))
        end
        if boundary_only
            bfes = meshboundary(femm.integdomain.fes);
            vtkexportmesh(rfile, fens, bfes;  vectors=vectors)
        else
            vtkexportmesh(rfile, fens, femm.integdomain.fes; vectors=vectors)
        end
        ed = FDataDict("file"=>rfile, "field"=>u, "region"=>i,
            "type"=>"displacement")
        push!(modeldata["postprocessing"]["exported"], ed)
    end

    return modeldata
end

"""
    AlgoDeforLinearModule.exportstress(modeldata::FDataDict)

Algorithm for exporting of the stress for visualization in Paraview.

`modeldata` = dictionary with values for keys

* "fens"  = finite element node set
* "regions"  = array of region dictionaries
* "geom" = geometry field
* "u" = displacement field
* "postprocessing" = dictionary  with values for keys
  + "boundary_only" = should only the boundary of the  regions be rendered?
                      Default is render the interior.
  + "file" = name of the  postprocessing file
  + "quantity" = quantity to be exported (default :Cauchy)
  + "component" = which component of the quantity?
  + "outputcsys" = output coordinate system
  + "inspectormeth" = inspector method to pass to `inspectintegpoints()`
  + "extrap" = method for extrapolating from the quadrature points to the nodes
    within one element

For each region (connected piece of the domain made of a particular material),
mandatory, the  region dictionary  contains values for keys:
* "femm" = finite element mmodel machine (mandatory);

Output: modeldata updated with
* modeldata["postprocessing"]["exported"] = array of data dictionaries, one for
        each exported file. The data is stored with the keys:
    - "file" - name of exported file
    - "field" - nodal field
"""
function exportstress(modeldata::FDataDict)
    modeldata_recognized_keys = ["fens", "regions", "geom", "u",
    "dT", "postprocessing"]
    postprocessing_recognized_keys = ["boundary_only", "file", "quantity",
    "component", "outputcsys", "nodevalmethod", "reportat"]
    # Defaults
    boundary_only = false;
    ffile = "stress"
    dcheck!(modeldata, modeldata_recognized_keys)
    quantity = :Cauchy
    component = 1
    outputcsys = nothing
    reportat = :default
    nodevalmethod = :invdistance
    # Let's have a look at what's been specified
    postprocessing = get(modeldata, "postprocessing", nothing);
    if (postprocessing != nothing)
        dcheck!(postprocessing, postprocessing_recognized_keys)
        boundary_only = get(postprocessing, "boundary_only", boundary_only);
        ffile = get(postprocessing, "file", ffile);
        quantity = get(postprocessing, "quantity", quantity);
        component = get(postprocessing, "component", component);
        outputcsys = get(postprocessing, "outputcsys", outputcsys);
        nodevalmethod = get(postprocessing, "nodevalmethod", nodevalmethod);
        reportat = get(postprocessing, "reportat", reportat);
    end

    fens = get(()->error("Must get fens!"), modeldata, "fens")
    geom = get(()->error("Must get geometry field!"), modeldata, "geom");
    u = get(()->error("Must get displacement field!"), modeldata, "u");
    dT = get(modeldata, "dT", nothing);

    context = []
    if (outputcsys != nothing)
        push!(context, (:outputcsys, outputcsys))
    end
    if (nodevalmethod != nothing)
        push!(context, (:nodevalmethod, nodevalmethod))
    end
    if (reportat != nothing)
        push!(context, (:reportat, reportat))
    end

    # Export a file for each region
    modeldata["postprocessing"]["exported"] = Array{FDataDict, 1}()
    regions = get(()->error("Must get region!"), modeldata, "regions")
    for i = 1:length(regions)
        region = regions[i]
        femm = region["femm"]
        rfile = ffile * "-" * string(quantity) * string(component) * "-region $i" * ".vtk";
        if (typeof(component) == Symbol)
            componentnum = stresscomponentmap(femm.mr)[component]
        else
            componentnum = component
        end
        componentname = length(componentnum) > 1 ? "" : "$(componentnum)"
        # Note that we are creating a field  separately for each region.  This is
        # important  for the following reason: if the regions were of different
        # materials, or if they were of the same material but with different material
        # axes orientation, averaging across the material interface  would not make
        # sense.
        if (dT != nothing)
            fld = fieldfromintegpoints(femm, geom, u, dT, quantity, componentnum;
            context...)
        else
            fld = fieldfromintegpoints(femm, geom, u, quantity, componentnum;
            context...)
        end
        if boundary_only
            bfes = meshboundary(femm.integdomain.fes);
            vtkexportmesh(rfile, fens, bfes;
            scalars=[(string(quantity)*componentname, fld.values)],
            vectors=[("u", u.values)])
        else
            vtkexportmesh(rfile, fens, femm.integdomain.fes;
            scalars=[(string(quantity)*componentname, fld.values)],
            vectors=[("u", u.values)])
        end
        ed = FDataDict("file"=>rfile, "field"=>fld, "region"=>i,
            "type"=>"nodal stress",
            "quantity"=>quantity, "component"=>component, "outputcsys"=>outputcsys,
            "nodevalmethod"=>nodevalmethod, "reportat"=>reportat)
        push!(modeldata["postprocessing"]["exported"], ed)
    end

    return modeldata
end

"""
    AlgoDeforLinearModule.exportstresselementwise(modeldata::FDataDict)

Algorithm for exporting of the elementwise stress for visualization in Paraview.

`modeldata` = dictionary with values for keys

* "fens"  = finite element node set
* "regions"  = array of region dictionaries
* "geom" = geometry field
* "u" = displacement field
* "postprocessing" = dictionary  with values for keys
  + "boundary_only" = should only the boundary of the  regions be rendered?
                      Default is render the interior.
  + "file" = name of the  postprocessing file
  + "quantity" = quantity to be exported (default :Cauchy)
  + "component" = which component of the quantity?
  + "outputcsys" = output coordinate system

For each region (connected piece of the domain made of a particular material),
mandatory, the  region dictionary  contains values for keys:
* "femm" = finite element mmodel machine (mandatory);

Output: modeldata updated with
* modeldata["postprocessing"]["exported"] = array of data dictionaries, one for
        each exported file. The data is stored with the keys:
    - "file" - name of exported file
    - "field" - elemental field
"""
function exportstresselementwise(modeldata::FDataDict)
    modeldata_recognized_keys = ["fens", "regions", "geom", "u",
    "dT", "postprocessing"]
    postprocessing_recognized_keys = ["boundary_only", "file", "quantity",
    "component", "outputcsys" ]
    # Defaults
    boundary_only = false;
    ffile = "stress"
    dcheck!(modeldata, modeldata_recognized_keys)
    quantity = :Cauchy
    component = 1
    outputcsys = nothing
    # Let's have a look at what's been specified
    postprocessing = get(modeldata, "postprocessing", nothing);
    if (postprocessing != nothing)
        dcheck!(postprocessing, postprocessing_recognized_keys)
        boundary_only = get(postprocessing, "boundary_only", boundary_only);
        ffile = get(postprocessing, "file", ffile);
        quantity = get(postprocessing, "quantity", quantity);
        component = get(postprocessing, "component", component);
        outputcsys = get(postprocessing, "outputcsys", outputcsys);
    end

    fens = get(()->error("Must get fens!"), modeldata, "fens")
    geom = get(()->error("Must get geometry field!"), modeldata, "geom");
    u = get(()->error("Must get displacement field!"), modeldata, "u");
    dT = get(modeldata, "dT", nothing);

    context = []
    if (outputcsys != nothing)
        push!(context, (:outputcsys, outputcsys))
    end

    # Export a file for each region
    modeldata["postprocessing"]["exported"] = Array{FDataDict, 1}()
    regions = get(()->error("Must get region!"), modeldata, "regions")
    for i = 1:length(regions)
        region = regions[i]
        femm = region["femm"]
        rfile = ffile * "-" * string(quantity) * string(component) * "-region $i" * ".vtk";
        if (typeof(component) == Symbol)
            componentnum = stresscomponentmap(femm.mr)[component]
        else
            componentnum = component
        end
        componentname = length(componentnum) > 1 ? "" : "$(componentnum)"
        # Note that we are creating a field  separately for each region.  This is
        # important  for the following reason: if the regions were of different
        # materials, or if they were of the same material but with different material
        # axes orientation, averaging across the material interface  would not make
        # sense.
        if (dT != nothing)
            fld = elemfieldfromintegpoints(femm, geom, u, dT, quantity, componentnum;
            context...)
        else
            fld = elemfieldfromintegpoints(femm, geom, u, quantity, componentnum;
            context...)
        end
        if boundary_only
            bfes = meshboundary(femm.integdomain.fes);
            vtkexportmesh(rfile, fens, bfes;
            scalars=[(string(quantity)*componentname, fld.values)],
            vectors=[("u", u.values)])
        else
            vtkexportmesh(rfile, fens, femm.integdomain.fes;
            scalars=[(string(quantity)*componentname, fld.values)],
            vectors=[("u", u.values)])
        end
        ed = FDataDict("file"=>rfile, "field"=>fld, "region"=>i,
            "type"=>"elemental stress",
            "quantity"=>quantity, "component"=>component, "outputcsys"=>outputcsys)
        push!(modeldata["postprocessing"]["exported"], ed)
    end

    return modeldata
end

"""
    AlgoDeforLinearModule.modal(modeldata::FDataDict)

Modal (free-vibration) analysis solver.

`modeldata` = dictionary with values for keys

* "fens"  = finite element node set
* "regions"  = array of region dictionaries
* "essential_bcs" = array of essential boundary condition dictionaries

For each region (connected piece of the domain made of a particular material),
mandatory, the  region dictionary  contains values for keys:
* "femm" = finite element mmodel machine (mandatory);

For essential boundary conditions (optional) each dictionary
would hold
  + "displacement" = fixed (prescribed) displacement (scalar): only zero
        displacement is  allowed for modal analysis.
  + "component" = which component is prescribed  (1, 2, 3)?
  + "node_list" = list of nodes on the boundary to which the condition applies
            (mandatory)

Control parameters:
* "neigvs" = number of eigenvalues/eigenvectors to compute
* "omega_shift"= angular frequency shift for mass shifting
* "use_lumped_mass" = true or false?  (Default is false: consistent mass)

Output:
modeldata= the dictionary on input is augmented with
* "geom" = the nodal field that is the geometry
* "u" = the nodal field that is the computed displacement
* "neigvs" = Number of computed eigenvectors
* "W" = Computed eigenvectors, neigvs columns
* "omega" =  Computed angular frequencies, array of length neigvs
* "raw_eigenvalues" = Raw computed eigenvalues
"""
function modal(modeldata::FDataDict)

    # For multi point constraints (MPC) (optional):
    # model_data.mpc= cell array of structs, each for one MPC.
    #      mpc.node_list = list of node numbers involved in the MPC,
    #      mpc.dof_list= numbers of degrees of freedom for the nodes above,
    #      mpc.umultipliers=multipliers for the nodes above,
    #      mpc.penfact=the penalty factor to multiply  the constraint matrix,
    #          The MPC looks like this: sum_i m_i u_{dof(i),node(i)} =0
    #          where m_i is the multiplier.

    # Lists of recognized keys for the data dictionaries:
    modeldata_recognized_keys = ["fens", "regions",
    "essential_bcs", "neigvs", "omega_shift", "use_lumped_mass"]
    essential_bcs_recognized_keys = ["displacement", "node_list", "component"]
    regions_recognized_keys = ["femm", "femm_stiffness", "femm_mass", "body_load"]

    neigvs = get(modeldata, "neigvs", 7); # Number of eigenvalues

    omega_shift = get(modeldata, "omega_shift", 0.0); # Mass shifting

    use_factorization = get(modeldata, "use_factorization", false); # Factorization?

    use_lumped_mass = get(modeldata, "use_lumped_mass", false); # Lumped mass?

    # Extract the nodes
    fens = get(()->error("Must get fens!"), modeldata, "fens")

    # Construct the geometry field
    geom = NodalField(fens.xyz)

    # Construct the displacement field
    u = NodalField(zeros(nnodes(geom),ndofs(geom)))

    # Apply the essential boundary conditions on the displacement field
    essential_bcs = get(modeldata, "essential_bcs", nothing);
    if (essential_bcs != nothing)
        for j = 1:length(essential_bcs)
            ebc = essential_bcs[j]
            dcheck!(ebc, essential_bcs_recognized_keys)
            fenids = get(()->error("Must get node list!"), ebc, "node_list");
            displacement = get(ebc, "displacement", nothing);
            u_fixed = zeros(FFlt, length(fenids)); # only zero displacement accepted
            component = get(ebc, "component", 0); # which component?
            setebc!(u, fenids[:], true, component, u_fixed);
        end
        applyebc!(u);
    end

    # Number the equations
    numberdofs!(u)           #,Renumbering_options); # NOT DONE <<<<<<<<<<<<<<<<<

    # Construct the system stiffness matrix
    K = spzeros(u.nfreedofs,u.nfreedofs); # (all zeros, for the moment)
    regions = get(()->error("Must get region list!"), modeldata, "regions")
    for i = 1:length(regions)
        region = regions[i]
        dcheck!(region, regions_recognized_keys)
        if "femm_stiffness"  in keys(region)
            femm = region["femm_stiffness"];
        else
            femm = get(()->error("Must get femm or femm_stiffness!"), region, "femm")
        end
        # Give the  FEMM a chance  to precompute  geometry-related quantities
        femm = associategeometry!(femm, geom);
        # Add up all the stiffness matrices for all the regions
        K = K + stiffness(femm, geom, u);
    end

    # Construct the system mass matrix
    M = spzeros(u.nfreedofs,u.nfreedofs); # (all zeros, for the moment)
    regions = get(()->error("Must get region list!"), modeldata, "regions")
    for i = 1:length(regions)
        region = regions[i]
        dcheck!(region, regions_recognized_keys)
        if "femm_mass"  in keys(region)
            femm = region["femm_mass"];
        else
            femm = get(()->error("Must get femm or femm_mass!"), region, "femm")
        end
        # Give the  FEMM a chance  to precompute  geometry-related quantities
        femm = associategeometry!(femm, geom);
        # Add up all the mass matrices for all the regions
        M = M + mass(femm, geom, u);
    end

    # Options for the eigenproblem solution

    # Solve
    # if (~ use_factorization )
    #     # This is one way of solving the eigenvalue problem, just pass the matrices
    #     [W,Omega]= eigs(K+omega_shift*M, M, neigvs, 'SM', evopts);
    # else
    # This form uses the factorized matrix and has the potential of being much faster
    # Factorize the left-hand side matrix for efficiency (Choleski)
    # [mA,status] = chol(K+omega_shift*M,'lower');#,'vector',prm
    # if ( status ~= 0 ) error('Choleski factorization failed'), end
    # clear K; # Not needed anymore
    # mAt= mA';
    # [W,Omega]= eigs(@(bv)mAt\(mA\bv), u.nfreedofs, M, neigvs, 'SM', evopts);
    #          [W,Omega]= eigen(full(K+omega_shift*M), full(M));

    d,v,nev,nconv = eigs(K+omega_shift*M, M; nev=neigvs, which=:SM)
    broadcast!(+, d, d, -omega_shift);

    modeldata["raw_eigenvalues"] = d;
    #    Subtract the mass-shifting Angular frequency
    if any(imag(d) .!= 0.0)
        d=real(d);
    end
    if any(real(d) .< 0.0)
        d = abs.(d);
    end
    #    Sort  the angular frequencies by magnitude.  Make sure all
    #    imaginary parts of the eigenvalues are removed.
    ix =sortperm(d);

    # Update the model data: store geometry
    modeldata["geom"] = geom;
    # Store the displacement field
    modeldata["u"] = u;
    # Number of computed eigenvectors
    modeldata["neigvs"] = nev;
    #  Computed eigenvectors: we are ignoring the imaginary part here
    #  because the modal analysis is presumed to have been performed for
    #  an undamped structure
    modeldata["W"] = real(v[:,ix]);
    #  Computed angular frequencies
    modeldata["omega"] = sqrt.(d[ix]);
    return modeldata
end

"""
    AlgoDeforLinearModule.exportmode(modeldata::FDataDict)

Algorithm for exporting of the mmode shape for visualization in Paraview.

`modeldata` = dictionary with values for keys

* "fens"  = finite element node set
* "regions"  = array of region dictionaries
* "geom" = geometry field
* "u" = displacement field
* "W" = Computed free-vibration eigenvectors, neigvs columns
* "omega" =  Computed free-vibration angular frequencies, array of length neigvs
* "postprocessing" = dictionary  with values for keys
  + "boundary_only" = should only the boundary of the  regions be rendered?
                      Default is render the interior.
  + "file" = name of the  postprocessing file
  + "mode" = which mode should be visualized?
  + "component" = which component of the quantity?
  + "outputcsys" = output coordinate system

For each region (connected piece of the domain made of a particular material),
mandatory, the  region dictionary  contains values for keys:
* "femm" = finite element mmodel machine (mandatory);

Output: modeldata updated with
* modeldata["postprocessing"]["exported"] = see `exportdeformation()`
"""
function exportmode(modeldata::FDataDict)
    modeldata_recognized_keys = ["fens", "regions", "geom", "u",
    "omega", "W",
    "postprocessing"]
    postprocessing_recognized_keys = ["boundary_only", "file", "mode"]
    mode = 1;
    dcheck!(modeldata, modeldata_recognized_keys)

    # Let's have a look at what's been specified
    postprocessing = get(modeldata, "postprocessing", nothing);
    if (postprocessing != nothing)
        dcheck!(postprocessing, postprocessing_recognized_keys)
        mode =  get(postprocessing, "mode", mode);
    end

    omega = modeldata["omega"]

    # Scatter the desired mode
    W = modeldata["W"]
    if typeof(mode)<:Int
        @assert 0 < mode <= length(omega) "Invalid mode number $mode"
        scattersysvec!(modeldata["u"], W[:,mode])
    else
        us = Tuple{String, Field}[]
        u = modeldata["u"]
        for ixxxx in mode
            @assert 0 < ixxxx <= length(omega) "Invalid mode number $ixxxx"
            scattersysvec!(u, W[:,ixxxx])
            push!(us, ("mode_$(ixxxx)", deepcopy(u)))
        end
        modeldata["us"] = us
    end

    return exportdeformation(modeldata)
end

"""
    ssit(K, M; nev::Int=6, evshift::FFlt = 0.0,
        v0::FFltMat = Array{FFlt}(0, 0),
        tol::FFlt = 1.0e-3, maxiter::Int = 300, verbose::Bool=false)

Subspace  Iteration (block inverse power) method.

Block inverse power method for k smallest eigenvalues of the generalized
eigenvalue problem
           `K*v= lambda*M*v`

# Arguments
* `K` =  square symmetric stiffness matrix (if necessary mass-shifted),
* `M` =  square symmetric mass matrix,

# Keyword arguments
* `v0` =  initial guess of the eigenvectors (for instance random),
* `nev` = the number of eigenvalues sought
* `tol` = relative tolerance on the eigenvalue, expressed in terms of norms of the
      change of the eigenvalue estimates from iteration to iteration.
* `maxiter` =  maximum number of allowed iterations
* `withrr` = with Rayleigh-Ritz problem solved to improve the subspace?  (default
    is false)
* `verbose` = verbose? (default is false)

#  Return
* `labm` = computed eigenvalues,
* `v` = computed eigenvectors,
* `nconv` = number of converged eigenvalues
* `niter` = number of iterations taken
* `nmult` = ignore this output
* `lamberr` = eigenvalue errors, defined as  normalized  differences  of
    successive  estimates of the eigenvalues
"""
function ssit(K, M; nev::Int=6, evshift::FFlt = 0.0,
    v0::FFltMat = Array{FFlt}(0, 0),
    tol::FFlt = 1.0e-3, maxiter::Int = 300, withrr::Bool=false,
    verbose::Bool=false)
    @assert nev >= 1
    v = deepcopy(v0)
    if isempty(v0)
        v = rand(size(K, 1), nev)
    end
    @assert nev <= size(v, 2)
    nvecs = size(v, 2)  # How many eigenvalues are iterated?
    plamb = zeros(nvecs)  # previous eigenvalue
    lamb = zeros(nvecs)
    lamberr = zeros(nvecs)
    converged = falses(nvecs)  # not yet
    niter = 0
    nconv = 0
    nmult = 0
    factor = cholesky(K+evshift*M)
    Kv = zeros(size(K, 1), size(v, 2))
    Mv = zeros(size(M, 1), size(v, 2))
    for i = 1:maxiter
        u = factor\(M*v)
        factorization = qr(u)  # ; full=falseeconomy factorization
        v = Array(factorization.Q)
        my_A_mul_B!(Kv, K, v)
        my_A_mul_B!(Mv, M, v)
        for j = 1:nvecs
            lamb[j] = dot(v[:, j], Kv[:, j]) / dot(v[:, j], Mv[:, j])
            lamberr[j] = abs(lamb[j] - plamb[j])/abs(lamb[j])
            converged[j] = lamberr[j] <= tol
        end
        nconv = length(findall(converged[1:nev]))
        verbose && println("nconv = $(nconv)")
        if nconv >= nev # converged on all requested eigenvalues
            break
        end
        if withrr
            decomp = eigen(transpose(v)*Kv, transpose(v)*Mv)
            ix = sortperm(abs.(decomp.values))
            rrd = decomp.values[ix]
            rrv = decomp.vectors[:, ix]
            v = v*rrv
        end
        plamb, lamb = lamb, plamb # swap the eigenvalue arrays
        niter = niter + 1
    end
    return lamb, v, nconv, niter, nmult, lamberr
end
# (d,[v,],nconv,niter,nmult,resid)
# eigs returns the nev requested eigenvalues in d, the corresponding Ritz vectors
# v (only if ritzvec=true), the number of converged eigenvalues nconv, the number
# of iterations niter and the number of matrix vector multiplications nmult, as
# well as the final residual vector resid.

end


text:
####################
# Continuous Lagrange
####################

struct ContinuousLagrange{dim,shape,qorder,gorder} <: FiniteElement{dim,shape,qorder,gorder}
    nodal_base_coefs::Matrix{Float64}
    geom_base_coefs::Matrix{Float64}
    topology::Dict{Int,Int}
    geom_topology::Dict{Int,Int}
    M::Matrix{Float64}
end

function Base.show(io::IO, fe::ContinuousLagrange{dim,shape,qorder,gorder}) where {dim, shape, qorder, gorder}
  println(io, "$dim D Continuous Lagrange Finite Element")
  println(io, "order: $qorder")
  println(io, "Reference Shape: ",split("$shape",".")[end], " with geometric order: $gorder")
end

function ContinuousLagrange(shape::Symbol,qorder)
  shape = map_shape_symbols(shape)
  ContinuousLagrange(getdim(shape), typeof(shape), qorder,1)
end
ContinuousLagrange(shape::Type{s},qorder) where {s <: Shape} = ContinuousLagrange(getdim(shape()), shape, qorder,1)

function ContinuousLagrange(dim,shape,qorder,gorder)
    nodal_base_coefs, topology = _nodal_data(dim, shape, qorder)
    if qorder == gorder
        geom_base_coefs = nodal_base_coefs
        geom_topology = topology
    else
        geom_base_coefs, geom_topology = _nodal_data(dim, shape, gorder)
    end
    if gorder == 1
        M = one(geom_base_coefs)
    else
        M = _nodal_geom_data(dim, shape, gorder, size(geom_base_coefs,1), geom_base_coefs)
    end
    ContinuousLagrange{dim, shape, qorder, gorder}(nodal_base_coefs, geom_base_coefs, topology, geom_topology, M)
end

function _nodal_data(dim::Int,shape, order::Int)
    nodal_points, topology = get_nodal_points(shape, order)
    ip_prime = getdefaultdualbasis(shape,order)
    nbasefuncs = getnbasefunctions(ip_prime)
    V = [value(ip_prime, j, nodal_points[i]) for i=1:nbasefuncs,j = 1:nbasefuncs]   #l_i(ϕ) = ϕ(x_i)
    return inv(V), topology
end

function _nodal_geom_data(dim::Int,shape::Type{Shape{N}}, order::Int, n_geom_basefuncs::Int, geom_base_coefs::Matrix) where {N}
    # Matrix to get spacial coordinates
    nodal_points, topology = get_nodal_points(shape, order)
    T = eltype(nodal_points[1])
    qrs = QuadratureRule{shape,N,T}(fill(T(NaN), length(nodal_points)), nodal_points) # weights will not be used
    n_qpoints = length(getweights(qrs))
    M =    fill(zero(T)           * T(NaN), n_geom_basefuncs, n_qpoints)
    for (qp, ξ) in enumerate(qrs.points)
        for i in 1:n_geom_basefuncs
            M[i, qp] = Tensors.dot(geom_base_coefs[:,i], value(getdefaultdualbasis(shape,order), ξ))
        end
    end
    M
end

@inline getnbasefunctions(::ContinuousLagrange{1,Segment,order}) where {order} = order + 1
@inline getnbasefunctions(::ContinuousLagrange{2,Triangle,order}) where {order} = Int((order+1)*(order+2)/2)
@inline getnbasefunctions(::ContinuousLagrange{3,Tetrahedron,order}) where {order} = Int((order+1)*(order+2)*(order+3)/6)
@inline getnbasefunctions(::ContinuousLagrange{dim,HyperCube{dim},order}) where {dim,order} = (order+1)^dim

@inline getngeombasefunctions(::ContinuousLagrange{1,Segment,order,gorder}) where {order,gorder} = gorder + 1
@inline getngeombasefunctions(::ContinuousLagrange{2,Triangle,order,gorder}) where {order,gorder} = Int((gorder+1)*(gorder+2)/2)
@inline getngeombasefunctions(::ContinuousLagrange{3,Tetrahedron,order,gorder}) where {order,gorder} = Int((gorder+1)*(gorder+2)*(gorder+3)/6)
@inline getngeombasefunctions(::ContinuousLagrange{dim,HyperCube{dim},order,gorder}) where {dim,order,gorder} = (gorder+1)^dim
@inline gettopology(ip::ContinuousLagrange) = ip.topology
gettopology(mesh::AbstractPolytopalMesh,cell, ip::ContinuousLagrange) = ip.topology
@inline getgeomtopology(ip::ContinuousLagrange) = ip.geom_topology

"""
value(ip::ContinuousLagrange{dim,shape,order}, k::Int, ξ::Vec{dim,T}) where {dim,shape,order, T}
Compute value of Continuous Lagrange Finite Element basis `j` at point ξ on shape `shape`
"""
function value(ip::ContinuousLagrange{dim,shape,order}, k::Int, ξ::Tensors.Vec{dim,T}) where {dim,shape,order, T}
    Tensors.dot(ip.nodal_base_coefs[:,k], value(getdefaultdualbasis(shape,order), ξ))
end

"""
geom_value(ip::ContinuousLagrange{dim,shape,order}, k::Int, ξ::Vec{dim,T}) where {dim,shape,order, T}
Compute value of Continuous Lagrange Finite Element basis `j` at point ξ on shape `shape`
"""
function geom_value(ip::ContinuousLagrange{dim,shape,order,gorder}, k::Int, ξ::Tensors.Vec{dim,T}) where {dim,shape,order,gorder, T}
    Tensors.dot(ip.geom_base_coefs[:,k], value(getdefaultdualbasis(shape,gorder), ξ))
end

ext: jl

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
module Typeclass
export @class, @instance, @instance!, Typeclass_stub, Class

#for clarity
import Base.Meta.quot
isSymbol(x::Symbol)=true
isSymbol(x)=false
isQuote(x::Expr)=x.head==:quote
isQuote(x)=false
isDoubleColon(x::Expr)=x.head==:(::)
isDoubleColon(x)=false
isCall(x::Expr)=x.head==:call
isCall(x)=false
isFunction(x::Expr)=x.head==:function
isFunction(x)=false
isCurly(x::Expr)=x.head==:curly
isCurly(x)=false
isAssign(x::Expr)=x.head==:(=)
isAssign(x)=false

#Parse the type signature of a given function call
function parse_sig(call_expr)
    if call_expr.head != :call
        error("Expected call expression: $call_expr")
    else
        sym=call_expr.args[1]

        typesig=[]
        for typex in call_expr.args[2:end]
            if isDoubleColon(typex)
                #parse correctly for foo(::bar) and untyped
                s=length(typex.args)==2 ? typex.args[2] : typex.args[1]
                push!(typesig,s)
            else
                push!(typesig,:Any)
            end
        end
        (sym::Symbol,typesig)
    end
end

#clean up input; the supplied return values are not used right now, so not really important
unblock(b::Expr) = b.head==:block && length(b.args)==2 && b.args[1].head==:line? b.args[2] :error("Unblock is blocked")

#Parse a function declaration
function parse_declaration(dec::Expr)
    if isCall(dec) #foo(x::Int,y)
        return_type = Any
        call_expr = dec
        default_impl = nothing
    elseif isDoubleColon(dec) #foo(x::Int,y)::Int
        return_type = dec.args[2]
        call_expr = dec.args[1]
        default_impl = nothing
    #foo(x::Int,y) = (x+y) [::Int]
    elseif isAssign(dec)
        call_expr = dec.args[1]
        RHS = unblock(dec.args[2]) #was mishandling blocks, which are generally generated.
        if isDoubleColon(RHS) #foo(x::Int,y) = (x+y)::Int
            return_type = RHS.args[2]
            default_impl = RHS.args[1]
        else #foo(x::Int,y) = (x+y)
            return_type = Any
            default_impl = RHS
        end
    elseif isFunction(dec) #currently ignores return type info
        call_expr = dec.args[1]
        RHS = dec.args[2]
        return_type = Any
        default_impl = RHS
    else
        error("Unsupported function declaration $(dec), head must be :(::), :(=), :call, or :function")
    end
    call_expr,return_type,default_impl
end

immutable Typeclass_stub{T} end

function _class_code(args...)
    declarations=filter!(x->x.head!=:line,args[end].args)
    class_name=args[1]
    type_parameters=args[2:end-1]
    parsed_declarations = map(parse_declaration,declarations)
    class_decl=Expr(:(=),class_name,Class(class_name,type_parameters,parsed_declarations)) #was final line

    #stubs provide reflection info and prevent problems with derefing later
    funsyms = [ce.args[1] for (ce,j,k) in parsed_declarations]
    stub_block=Expr(:block)
    for sym in funsyms
        stub_decl=Expr(:(=),Expr(:call,sym,Expr(:(::),Expr(:curly,:Typeclass_stub,quot(class_name)))),true)
        push!(stub_block.args,stub_decl)
    end
    Expr(:block,class_decl,stub_block,nothing)
end

type Class
    class_name::Symbol
    type_parameters #symbols or curlies
    declarations #triple call_expr,return_type, default_implementation
end

import Base.show
function show(io::IO,C::Class)
    println(io,"Typeclass ",C.class_name)
    println(io,"Type parameters ",C.type_parameters)
    println(io,"Declarations ",C.declarations)
end

macro class(args...)
    esc(_class_code(args...))
end

typesub(binding,s::Symbol) = haskey(binding,s)? binding[s] :s
typesub(binding,e::Expr)   = Expr(e.head,[typesub(binding,a) for a in e.args]...)
typesub(binding,x)         = x

function _instance_code(overwrite,typeclass,implementing_type,args)
    instance_declarations = length(args)>2? filter!(x->x.head!=:line,args[end].args) : Any[] #allows bare "@instance Foo Bar" with no user declarations
    parsed_instance_declarations = map(parse_declaration,instance_declarations)
    instance_declared_sigs = [parse_sig(pid[1]) for pid in parsed_instance_declarations]
    class_name=args[1]
    implementing_type_name=args[2]

    if length(args)>2
        bindings = Dict(zip(typeclass.type_parameters,args[2:end-1]))
        implementing_type_names=args[2:end-1]
    else
        bindings = Dict(zip(typeclass.type_parameters,args[2:end])) #end=2
        implementing_type_names=args[2:end]
    end

    instance_block=Expr(:block)
    for d in instance_declarations
        push!(instance_block.args,d)     #run raw user declarations
    end

    class_block=Expr(:block)
    for parsed_dec in typeclass.declarations
        call_expr=parsed_dec[1]
        return_type=parsed_dec[2]
        default_impl=parsed_dec[3]

        LHS=typesub(bindings,call_expr)
        fsym,fsig=parse_sig(LHS)
        if !(parse_sig(LHS) in instance_declared_sigs) #if specialized, parsed_dec was not declared by user
            default_impl==nothing? error("Class method $LHS required by $class_name not implemented in instance declaration for $implementing_type_name") :nothing
            RHS=typesub(bindings,default_impl)
            decl=Expr(:(=),LHS,RHS)
            check=Expr(:call,isempty,Expr(:call,:methods,fsym,Expr(:tuple,fsig...))) #will error if fsym is undef, and try blocks cause hyg probs, hence stub
            warning=Expr(:call,:info,"Instance delcaration omitted method $(LHS) is already defined, skipping specializing $class_name class method")
            checked_decl = Expr(:if, check, decl, warning)
            push!(class_block.args, overwrite? decl : checked_decl)
        end
    end

    #allow adding Typeclass constraints to functions and types, as well as runtime checks
    registration_block=Expr(:block) #uses Graphs.jl method
    reg_LHS=Expr(:call,Symbol(string("implements_",string(class_name))), [Expr(:(::),implementing_type_name) for implementing_type_name in implementing_type_names]...)
    push!(registration_block.args,Expr(:(=),reg_LHS,:true))

    Expr(:block,registration_block,instance_block,class_block,nothing)
end

macro instance(args...)
    typeclass=args[1]
    implementing_type=args[2]
    e=Expr(:call,_instance_code,false,typeclass,implementing_type,args) #stopping here just returns the code
    esc(Expr(:call,:eval,e)) #like calling eval( @instance ...) and stopping at previous line; added esc when made a module
end

macro instance!(args...) #version that overwrites existing methods
    typeclass=args[1]
    implementing_type=args[2]
    e=Expr(:call,_instance_code,true,typeclass,implementing_type,args)
    esc(Expr(:call,:eval,e))
end

macro ins(args...)
    global typeclass=args[1]
    println(typeclass)
    implementing_type=args[2]
    println(implementing_type)
    esc(Expr(:call,_instance_code,typeclass,implementing_type,args)) #stopping here just returns the code; added esc when made a module
end
#args=@id <decl>
#_instance_code(args[1],args[2],args) #to debug

end #module


text:
# Julia program to check if a number is a power of two.

function power_of_two(n)
    # 0 is not a power of 2
    if(n != 0)
        if((n & (n-1)) == false)
            return true
        end
    end
    return false
end

print("Enter a number: ")
n = readline()
n = parse(Int, n)
res = power_of_two(n)
if(res)
    print("The given number is a power of 2.")
else
    print("The given number is not a power of 2.")
end

"""
Time Complexity - O(1)
Space Complexity - O(1)

SAMPLE INPUT AND OUTPUT

SAMPLE I
Enter a number: 4
The given number is a power of 2.

SAMPLE II
Enter a number: 0
The given number is not a power of 2.
"""

max_issues_repo_name: adarshnjena/NeoAlgo

max_forks_repo_name: adarshnjena/NeoAlgo


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

module Versions

export VersionBound, VersionRange, VersionSpec, semver_spec, isjoinable

################
# VersionBound #
################
struct VersionBound
    t::NTuple{3,UInt32}
    n::Int
    function VersionBound(tin::NTuple{n,Integer}) where n
        n <= 3 || throw(ArgumentError("VersionBound: you can only specify major, minor and patch versions"))
        n == 0 && return new((0,           0,      0), n)
        n == 1 && return new((tin[1],      0,      0), n)
        n == 2 && return new((tin[1], tin[2],      0), n)
        n == 3 && return new((tin[1], tin[2], tin[3]), n)
        error("invalid $n")
    end
end
VersionBound(t::Integer...) = VersionBound(t)
VersionBound(v::VersionNumber) = VersionBound(v.major, v.minor, v.patch)

Base.getindex(b::VersionBound, i::Int) = b.t[i]

function ≲(v::VersionNumber, b::VersionBound)
    b.n == 0 && return true
    b.n == 1 && return v.major <= b[1]
    b.n == 2 && return (v.major, v.minor) <= (b[1], b[2])
    return (v.major, v.minor, v.patch) <= (b[1], b[2], b[3])
end

function ≲(b::VersionBound, v::VersionNumber)
    b.n == 0 && return true
    b.n == 1 && return v.major >= b[1]
    b.n == 2 && return (v.major, v.minor) >= (b[1], b[2])
    return (v.major, v.minor, v.patch) >= (b[1], b[2], b[3])
end

function isless_ll(a::VersionBound, b::VersionBound)
    m, n = a.n, b.n
    for i = 1:min(m, n)
        a[i] < b[i] && return true
        a[i] > b[i] && return false
    end
    return m < n
end

stricterlower(a::VersionBound, b::VersionBound) = isless_ll(a, b) ? b : a

# Comparison between two upper bounds
function isless_uu(a::VersionBound, b::VersionBound)
    m, n = a.n, b.n
    for i = 1:min(m, n)
        a[i] < b[i] && return true
        a[i] > b[i] && return false
    end
    return m > n
end

stricterupper(a::VersionBound, b::VersionBound) = isless_uu(a, b) ? a : b

# `isjoinable` compares an upper bound of a range with the lower bound of the next range
# to determine if they can be joined, as in [1.5-2.8, 2.5-3] -> [1.5-3]. Used by `union!`.
# The equal-length-bounds case is special since e.g. `1.5` can be joined with `1.6`,
# `2.3.4` can be joined with `2.3.5` etc.

function isjoinable(up::VersionBound, lo::VersionBound)
    up.n == 0 && lo.n == 0 && return true
    if up.n == lo.n
        n = up.n
        for i = 1:(n - 1)
            up[i] > lo[i] && return true
            up[i] < lo[i] && return false
        end
        up[n] < lo[n] - 1 && return false
        return true
    else
        l = min(up.n, lo.n)
        for i = 1:l
            up[i] > lo[i] && return true
            up[i] < lo[i] && return false
        end
    end
    return true
end

Base.hash(r::VersionBound, h::UInt) = hash(hash(r.t, h), r.n)

# Hot code
function VersionBound(s::AbstractString)
    s = strip(s)
    s == "*" && return VersionBound()
    first(s) == 'v' && (s = SubString(s, 2))
    l = lastindex(s)

    p = findnext('.', s, 1)
    b = p === nothing ? l : (p-1)
    i = parse(Int64, SubString(s, 1, b))
    p === nothing && return VersionBound(i)

    a = p+1
    p = findnext('.', s, a)
    b = p === nothing ? l : (p-1)
    j = parse(Int64, SubString(s, a, b))
    p === nothing && return VersionBound(i, j)

    a = p+1
    p = findnext('.', s, a)
    b = p === nothing ? l : (p-1)
    k = parse(Int64, SubString(s, a, b))
    p === nothing && return VersionBound(i, j, k)

    error("invalid VersionBound string $(repr(s))")
end

################
# VersionRange #
################
struct VersionRange
    lower::VersionBound
    upper::VersionBound
    # NOTE: ranges are allowed to be empty; they are ignored by VersionSpec anyway
    function VersionRange(lo::VersionBound, hi::VersionBound)
        # lo.t == hi.t implies that digits past min(lo.n, hi.n) are zero
        # lo.n < hi.n example: 1.2-1.2.0 => 1.2.0
        # lo.n > hi.n example: 1.2.0-1.2 => 1.2
        lo.t == hi.t && (lo = hi)
        return new(lo, hi)
    end
end
VersionRange(b::VersionBound=VersionBound()) = VersionRange(b, b)
VersionRange(t::Integer...)                  = VersionRange(VersionBound(t...))
VersionRange(v::VersionNumber)               = VersionRange(VersionBound(v))
VersionRange(lo::VersionNumber, hi::VersionNumber) = VersionRange(VersionBound(lo), VersionBound(hi))

# The vast majority of VersionRanges are in practice equal to "1"
const VersionRange_1 = VersionRange(VersionBound("1"), VersionBound("1"))
function VersionRange(s::AbstractString)
    s == "1" && return VersionRange_1
    p = split(s, "-")
    if length(p) != 1 && length(p) != 2
        throw(ArgumentError("invalid version range: $(repr(s))"))
    end
    lower = VersionBound(p[1])
    upper = length(p) == 1 ? lower : VersionBound(p[2])
    return VersionRange(lower, upper)
end

function Base.isempty(r::VersionRange)
    for i = 1:min(r.lower.n, r.upper.n)
        r.lower[i] > r.upper[i] && return true
        r.lower[i] < r.upper[i] && return false
    end
    return false
end

function Base.print(io::IO, r::VersionRange)
    m, n = r.lower.n, r.upper.n
    if (m, n) == (0, 0)
        print(io, '*')
    elseif m == 0
        print(io, "0-")
        join(io, r.upper.t, '.')
    elseif n == 0
        join(io, r.lower.t, '.')
        print(io, "-*")
    else
        join(io, r.lower.t[1:m], '.')
        if r.lower != r.upper
            print(io, '-')
            join(io, r.upper.t[1:n], '.')
        end
    end
end
Base.show(io::IO, r::VersionRange) = print(io, "VersionRange(\"", r, "\")")

Base.in(v::VersionNumber, r::VersionRange) = r.lower ≲ v ≲ r.upper

Base.intersect(a::VersionRange, b::VersionRange) = VersionRange(stricterlower(a.lower, b.lower), stricterupper(a.upper, b.upper))

function Base.union!(ranges::Vector{<:VersionRange})
    l = length(ranges)
    l == 0 && return ranges

    sort!(ranges, lt=(a, b) -> (isless_ll(a.lower, b.lower) || (a.lower == b.lower && isless_uu(a.upper, b.upper))))

    k0 = 1
    ks = findfirst(!isempty, ranges)
    ks === nothing && return empty!(ranges)

    lo, up, k0 = ranges[ks].lower, ranges[ks].upper, 1
    for k = (ks + 1):l
        isempty(ranges[k]) && continue
        lo1, up1 = ranges[k].lower, ranges[k].upper
        if isjoinable(up, lo1)
            isless_uu(up, up1) && (up = up1)
            continue
        end
        vr = VersionRange(lo, up)
        @assert !isempty(vr)
        ranges[k0] = vr
        k0 += 1
        lo, up = lo1, up1
    end
    vr = VersionRange(lo, up)
    if !isempty(vr)
        ranges[k0] = vr
        k0 += 1
    end
    resize!(ranges, k0 - 1)
    return ranges
end

###############
# VersionSpec #
###############
struct VersionSpec
    ranges::Vector{VersionRange}
    VersionSpec(r::Vector{<:VersionRange}) = new(union!(r))
    VersionSpec(vs::VersionSpec) = vs
end

VersionSpec(r::VersionRange) = VersionSpec(VersionRange[r])
VersionSpec(v::VersionNumber) = VersionSpec(VersionRange(v))
const _all_versionsspec = VersionSpec(VersionRange())
VersionSpec() = _all_versionsspec
VersionSpec(s::AbstractString) = VersionSpec(VersionRange(s))
VersionSpec(v::AbstractVector) = VersionSpec(map(VersionRange, v))

# Hot code
function Base.in(v::VersionNumber, s::VersionSpec)
    for r in s.ranges
        v in r && return true
    end
    return false
end

Base.copy(vs::VersionSpec) = VersionSpec(vs)

const empty_versionspec = VersionSpec(VersionRange[])
const _empty_symbol = "∅"

Base.isempty(s::VersionSpec) = all(isempty, s.ranges)
@assert isempty(empty_versionspec)
# Hot code, measure performance before changing
function Base.intersect(A::VersionSpec, B::VersionSpec)
    (isempty(A) || isempty(B)) && return copy(empty_versionspec)
    ranges = Vector{VersionRange}(undef, length(A.ranges) * length(B.ranges))
    i = 1
    @inbounds for a in A.ranges, b in B.ranges
        ranges[i] = intersect(a, b)
        i += 1
    end
    VersionSpec(ranges)
end
Base.intersect(a::VersionNumber, B::VersionSpec) = a in B ? VersionSpec(a) : empty_versionspec
Base.intersect(A::VersionSpec, b::VersionNumber) = intersect(b, A)

function Base.union(A::VersionSpec, B::VersionSpec)
    A == B && return A
    Ar = copy(A.ranges)
    append!(Ar, B.ranges)
    union!(Ar)
    return VersionSpec(Ar)
end

Base.:(==)(A::VersionSpec, B::VersionSpec) = A.ranges == B.ranges
Base.hash(s::VersionSpec, h::UInt) = hash(s.ranges, h + (0x2fd2ca6efa023f44 % UInt))

function Base.print(io::IO, s::VersionSpec)
    isempty(s) && return print(io, _empty_symbol)
    length(s.ranges) == 1 && return print(io, s.ranges[1])
    print(io, '[')
    for i = 1:length(s.ranges)
        1 < i && print(io, ", ")
        print(io, s.ranges[i])
    end
    print(io, ']')
end
Base.show(io::IO, s::VersionSpec) = print(io, "VersionSpec(\"", s, "\")")

###################
# Semver notation #
###################

function semver_spec(s::String; throw = true)
    ranges = VersionRange[]
    for ver in strip.(split(strip(s), ','))
        range = nothing
        found_match = false
        for (ver_reg, f) in ver_regs
            if occursin(ver_reg, ver)
                range = f(match(ver_reg, ver))
                found_match = true
                break
            end
        end
        if !found_match
            if throw
                error("invalid version specifier: \"$s\"")
            else
                return nothing
            end
        end
        push!(ranges, range)
    end
    return VersionSpec(ranges)
end

function semver_interval(m::RegexMatch)
    @assert length(m.captures) == 4
    n_significant = count(x -> x !== nothing, m.captures) - 1
    typ, _major, _minor, _patch = m.captures
    major =                           parse(Int, _major)
    minor = (n_significant < 2) ? 0 : parse(Int, _minor)
    patch = (n_significant < 3) ? 0 : parse(Int, _patch)
    if n_significant == 3 && major == 0 && minor == 0 && patch == 0
        error("invalid version: \"0.0.0\"")
    end
    # Default type is :caret
    vertyp = (typ == "" || typ == "^") ? :caret : :tilde
    v0 = VersionBound((major, minor, patch))
    if vertyp === :caret
        if major != 0
            return VersionRange(v0, VersionBound((v0[1],)))
        elseif minor != 0
            return VersionRange(v0, VersionBound((v0[1], v0[2])))
        else
            if n_significant == 1
                return VersionRange(v0, VersionBound((0,)))
            elseif n_significant == 2
                return VersionRange(v0, VersionBound((0, 0,)))
            else
                return VersionRange(v0, VersionBound((0, 0, v0[3])))
            end
        end
    else
        if n_significant == 3 || n_significant == 2
            return VersionRange(v0, VersionBound((v0[1], v0[2],)))
        else
            return VersionRange(v0, VersionBound((v0[1],)))
        end
    end
end

const _inf = VersionBound("*")
function inequality_interval(m::RegexMatch)
    @assert length(m.captures) == 4
    typ, _major, _minor, _patch = m.captures
    n_significant = count(x -> x !== nothing, m.captures) - 1
    major =                           parse(Int, _major)
    minor = (n_significant < 2) ? 0 : parse(Int, _minor)
    patch = (n_significant < 3) ? 0 : parse(Int, _patch)
    if n_significant == 3 && major == 0 && minor == 0 && patch == 0
        error("invalid version: 0.0.0")
    end
    v = VersionBound(major, minor, patch)
    if occursin(r"^<\s*$", typ)
        nil = VersionBound(0, 0, 0)
        if v[3] == 0
            if v[2] == 0
                v1 = VersionBound(v[1]-1)
            else
                v1 = VersionBound(v[1], v[2]-1)
            end
        else
            v1 = VersionBound(v[1], v[2], v[3]-1)
        end
        return VersionRange(nil, v1)
    elseif occursin(r"^=\s*$", typ)
        return VersionRange(v)
    elseif occursin(r"^>=\s*$", typ) || occursin(r"^≥\s*$", typ)
           return VersionRange(v, _inf)
    else
        error("invalid prefix $typ")
    end
end

function hyphen_interval(m::RegexMatch)
    @assert length(m.captures) == 6
    _lower_major, _lower_minor, _lower_patch, _upper_major, _upper_minor, _upper_patch = m.captures
    if isnothing(_lower_minor)
        lower_bound = VersionBound(parse(Int, _lower_major))
    elseif isnothing(_lower_patch)
        lower_bound = VersionBound(parse(Int, _lower_major),
                                   parse(Int, _lower_minor))
    else
        lower_bound = VersionBound(parse(Int, _lower_major),
                                   parse(Int, _lower_minor),
                                   parse(Int, _lower_patch))
    end
    if isnothing(_upper_minor)
        upper_bound = VersionBound(parse(Int, _upper_major))
    elseif isnothing(_upper_patch)
        upper_bound = VersionBound(parse(Int, _upper_major),
                                   parse(Int, _upper_minor))
    else
        upper_bound = VersionBound(parse(Int, _upper_major),
                                   parse(Int, _upper_minor),
                                   parse(Int, _upper_patch))
    end
    return VersionRange(lower_bound, upper_bound)
end

const version = "v?([0-9]+?)(?:\\.([0-9]+?))?(?:\\.([0-9]+?))?"
const ver_regs =
Pair{Regex,Any}[
                Regex("^([~^]?)?$version\$") => semver_interval, # 0.5 ^0.4 ~0.3.2
                Regex("^((?:≥\\s*)|(?:>=\\s*)|(?:=\\s*)|(?:<\\s*)|(?:=\\s*))v?$version\$")  => inequality_interval,# < 0.2 >= 0.5,2
                Regex("^[\\s]*$version[\\s]*?\\s-\\s[\\s]*?$version[\\s]*\$") => hyphen_interval, # 0.7 - 1.3
]

end


text:
module BinaryWrappers

export @generate_wrappers

using Scratch
using JLLWrappers

const wrapper_key = "binarywrappers_v$(VERSION.major).$(VERSION.minor)"

# use @generate_wrappers instead to automatically deduce the calling module
function generate_wrappers(m::Module, caller::Union{Module, Base.UUID, Nothing})
    # we generate wrappers per minor julia version
    # the scratch will belong to the jll which the wrappers are generated for
    # and the usage is tied to the module calling the `@generate` macro.
    target = get_scratch!(m, wrapper_key, caller)

    binpath(name) = joinpath(target, "bin", name)
    mkpath(binpath(""))

    bindir = joinpath(getproperty(m, :artifact_dir), "bin")
    sourcebinary = joinpath(bindir,"\$(basename \$0)")
    libpath = getproperty(m, :LIBPATH)[]

    wrapper = """
        #!/bin/sh
        # Since we cannot run these binaries through the usual julia commands we need
        # this wrapper that sets up the correct library paths.
        export $(JLLWrappers.LIBPATH_env)="$(libpath)"
        exec $(sourcebinary) "\$@"
        """

    # For shell scripts we use `source` (.) instead of exec to avoid macOS stripping
    # the LIBPATH. This means \$0 will point to our wrapper.
    # For 4ti2 this means that its own wrapper scripts will call our wrapper for the
    # binaries, so strictly speaking adjusting the LIBPATH here is not necessary.
    scriptwrapper = """
        # We cannot run shell scripts directly as macOS will remove
        # DYLD_FALLBACK_LIBRARY_PATH for any subshells.
        # So we source the original script instead.
        export $(JLLWrappers.LIBPATH_env)="$(libpath)"
        . $(sourcebinary) "\$@"
        """

    # POSIX compatible shells
    shellre = r"^#!/bin/(ba|da|z|k)?sh"

    for bin in readdir(bindir)
        if isfile(joinpath(bindir, bin))
            (tmpfile, tmpio) = mktemp(binpath(""),cleanup=false)
            # shell scripts use a different wrapper because macOS...
            m = match(shellre, String(read(joinpath(bindir, bin), 11)))
            if m != nothing
                write(tmpio, m.match * "\n" * scriptwrapper)
            else
                write(tmpio, wrapper)
            end
            close(tmpio)
            chmod(tmpfile, 0o755)
            # using mv would introduce some race conditions due to concurrent deletes
            Base.Filesystem.rename(tmpfile, binpath(bin))
        end
    end
    return binpath("")
end

# this behaves slightly different than the @get_scratch! function:
# it will associate the scratch to the module passed as argument
# and use the calling module for the scratch usage (for gc)
macro generate_wrappers(m::Symbol)
    uuid = Base.PkgId(__module__).uuid
    return quote
        generate_wrappers($(esc(m)), $(esc(uuid)))
    end
end

end


text:
## Defines functions to be used for beta scaling studies

# create a ferromagnet on a d random regular graph with or without weak impurity
function ferro_random_regular_model(N::Integer,d::Integer,α::Real,β::Real,FLAG_weak_impurity=false)
    #=
    # choose with equal probability those to be changed to -β
    # should select random numbers uniformly without replacement
    # also need to make sure resulting matrix is symmetric
    =#

    rand_reg_gm = random_regular_graph(N,d)
    temp_adj_matrix = adjacency_matrix(rand_reg_gm)

    U_temp_adj_matrix = triu(temp_adj_matrix)

    # adj_matrix is a sparse matrix
    R, C, V = findnz(U_temp_adj_matrix)
    V = β*V

    # one random index to change to -α
    ind_α = 4
    if FLAG_weak_impurity
        V[ind_α] = -α
    else
        V[ind_α] = α
    end

    U_adj_matrix = sparse(R, C, V, N, N)
    adj_matrix = U_adj_matrix + U_adj_matrix'
    adj_matrix = Array(adj_matrix)

    # Needed for plotting
    struct_adj_matrix = Array(temp_adj_matrix)

    return adj_matrix, struct_adj_matrix
end

# create a spin glass model on a d random regular graph
function spin_glass_random_regular_model(N::Integer,d::Integer,α::Real,β::Real)
    #=
    # choose with equal probability those to be changed to -β
    # should select random numbers uniformly without replacement
    # also need to make sure resulting matrix is symmetric
    =#

    rand_reg_gm = random_regular_graph(N,d)
    temp_adj_matrix = adjacency_matrix(rand_reg_gm)

    U_temp_adj_matrix = triu(temp_adj_matrix)

    # adj_matrix is a sparse matrix
    R, C, V = findnz(U_temp_adj_matrix)
    V = β*V

    # indices to be changed to negative β (randperm or randsubseq)
    # Ref: https://discourse.julialang.org/t/sampling-without-replacement/1073
    #ind_neg_β = randsubseq(1:length(V),0.5)
    ind_neg_β = randperm(length(V))[1:Int(N/2)]

    # Also find two random indices to change to +/- α
    ind_α = randperm(length(V))[1:2]

    V[ind_neg_β] = -β
    V[ind_α[1]] = α
    V[ind_α[2]] = -α

    U_adj_matrix = sparse(R, C, V, N, N)
    adj_matrix = U_adj_matrix + U_adj_matrix'
    adj_matrix = Array(adj_matrix)

    # Needed for plotting
    struct_adj_matrix = Array(temp_adj_matrix)

    return adj_matrix, struct_adj_matrix
end

# create a spin glass model on a double periodic lattice
function spin_glass_lattice(N::Integer,α::Real,β::Real)
    #=
    # choose with equal probability those to be changed to -β
    # should select random numbers uniformly without replacement
    # also need to make sure resulting matrix is symmetric
    =#
    temp_adj_matrix = double_periodic_lattice(Int(sqrt(N)),Int(sqrt(N)))
    temp_adj_matrix = sparse(temp_adj_matrix)

    U_temp_adj_matrix = triu(temp_adj_matrix)

    # adj_matrix is a sparse matrix
    R, C, V = findnz(U_temp_adj_matrix)
    V = β*V

    # indices to be changed to negative β (randperm or randsubseq)
    # Ref: https://discourse.julialang.org/t/sampling-without-replacement/1073
    #ind_neg_β = randsubseq(1:length(V),0.5)
    ind_neg_β = randperm(length(V))[1:Int(N/2)]

    # Also find two random indices to change to +/- α
    ind_α = randperm(length(V))[1:2]
    # HACK FOR NOW
    ind_α[1] = 4
    ind_α[2] = 28

    V[ind_neg_β] = -β
    V[ind_α[1]] = α
    V[ind_α[2]] = -α

    U_adj_matrix = sparse(R, C, V, N, N)
    adj_matrix = U_adj_matrix + U_adj_matrix'
    adj_matrix = Array(adj_matrix)

    # Needed for plotting
    struct_adj_matrix = Array(temp_adj_matrix)

    return adj_matrix, struct_adj_matrix
end

# replace values of beta in an existing spin glasss graphical model
function replace_beta_gm(adj_matrix::Array, α::Real, β::Real)
    mod_adj_matrix = deepcopy(adj_matrix)

    # Replace values greater than α
    mod_adj_matrix[adj_matrix.>α] .= β

    # Replace values smaller than -α
    mod_adj_matrix[adj_matrix.<-α] .= -β

    return mod_adj_matrix
end

# create a double periodic lattice of size (m,n)
function double_periodic_lattice(m::Integer,n::Integer)
    n_spins = m*n
    adj_matrix = zeros(n_spins,n_spins)

    for i = 1:(m*n)
        # top
        adj_matrix[i,( ((i-1) % m) == 0 ? (i+m-1) : i-1 )] = 1

        # bottom
        adj_matrix[i,( ((i+1) % m) == 1 ? (i+1-m) : i+1 )] = 1

        # left
        adj_matrix[i,( (i+m) > n_spins ? (i+m) % n_spins : i+m )] = 1

        # right
        adj_matrix[i,( (i-m) <= 0 ? (n_spins-m+i) : i-m )] = 1
    end

    return adj_matrix
end

# create a periodic ferromagnetic lattice with or without weak impurity
function ferro_lattice(N::Integer,α::Real,β::Real,FLAG_weak_impurity=false)
    temp_adj_matrix = double_periodic_lattice(Int(sqrt(N)),Int(sqrt(N)))
    temp_adj_matrix = sparse(temp_adj_matrix)

    U_temp_adj_matrix = triu(temp_adj_matrix)

    # adj_matrix is a sparse matrix
    R, C, V = findnz(U_temp_adj_matrix)
    V = β*V

    # one random index to change to -α
    #ind_α = rand(1:length(V))
    ind_α = 4
    if FLAG_weak_impurity
        V[ind_α] = -α
    else
        V[ind_α] = α
    end

    U_adj_matrix = sparse(R, C, V, N, N)
    adj_matrix = U_adj_matrix + U_adj_matrix'
    adj_matrix = Array(adj_matrix)

    # Needed for plotting
    struct_adj_matrix = Array(temp_adj_matrix)

    return adj_matrix, struct_adj_matrix
end

# create a non periodic lattice of size (m,n) for plotting purposes
function non_periodic_lattice(adj_matrix_gm::Array, m::Integer,n::Integer)
    # Map of the nodes in the plotted gm
    node_map_plot = Array{Int64,2}(undef,m,n)
    n_spins = m*n
    adj_matrix_plot = zeros(n_spins,n_spins)
    node_map_plot[1:n_spins] = 1:n_spins

    # Map of the indices in the plotted map to the original map
    values_map_plot = Array{Int64,2}(undef,m,n)
    values_map_plot[2:m-1,2:n-1] = 1:((m-2)*(n-2))

    # Specify edges between nodes at the boundary to be removed
    boundary_nodes = vcat(node_map_plot[:,1],node_map_plot[:,n],node_map_plot[1,2:n-1],node_map_plot[m,2:n-1])
    corner_nodes = [1,m,m*(n-1)+1,m*n]

    # Apply periodicity
    values_map_plot[2:(m-1),1] = copy(values_map_plot[2:(m-1),(n-1)])
    values_map_plot[2:(m-1),n] = copy(values_map_plot[2:(m-1),2])
    values_map_plot[1,2:(n-1)] = copy(values_map_plot[(m-1),2:(n-1)])
    values_map_plot[m,2:(n-1)] = copy(values_map_plot[2,2:(n-1)])
    values_map_plot[corner_nodes] = zeros(length(corner_nodes))

    # Create edges
    for i = 1:(m*n)
        if !(i in corner_nodes)
            # top
            if ((i-1) % m) != 0 && !((i-1) in corner_nodes)
                adj_matrix_plot[i,i-1] = adj_matrix[values_map_plot[i],values_map_plot[i-1]]
            end
            # bottom
            if ((i+1) % m) != 1 && !((i+1) in corner_nodes)
                adj_matrix_plot[i,i+1] = adj_matrix[values_map_plot[i],values_map_plot[i+1]]
            end
            # right
            if (i+m) <= n_spins && !((i+m) in corner_nodes)
                adj_matrix_plot[i,i+m] = adj_matrix[values_map_plot[i],values_map_plot[i+m]]
            end
            # left
            if (i-m) > 0 && !((i-m) in corner_nodes)
                adj_matrix_plot[i,i-m] = adj_matrix[values_map_plot[i],values_map_plot[i-m]]
            end
        end
    end
    for i in boundary_nodes
        for j in boundary_nodes
            adj_matrix_plot[i,j] = 0
        end
    end

    struct_adj_matrix_plot = sign.(abs.(adj_matrix_plot))

    return adj_matrix_plot, struct_adj_matrix_plot
end

# Function to draw a periodic lattice using Junet (no longer works in Julia 1.0)
function draw_periodic_lattice(adj_matrix::Array, m::Integer, n::Integer, α::Real, β::Real, gm_file::String)
    # Get the graph that Junet can plot
    adj_matrix_junet, struct_adj_matrix_junet = non_periodic_lattice(adj_matrix,m+2,n+2)

    # Create Junet Graph
    gm = Junet.Graph(struct_adj_matrix_junet,directed=false)

    # Layout
    N = nodecount(gm)
    x = Array{Float64,1}(N)
    y = Array{Float64,1}(N)

    node_map_gm = zeros(m+2,n+2)
    node_map_gm[1:N] = 1:N

    for i=1:N
        y[i], x[i] = ind2sub(node_map_gm,find(isequal(i),node_map_gm)[1])
        if i in node_map_gm[:,1]
            x[i] += 0.3
        elseif i in node_map_gm[:,n+2]
            x[i] += -0.3
        elseif i in node_map_gm[1,:]
            y[i] += +0.3
        elseif i in node_map_gm[m+2,:]
            y[i] += -0.3
        end
    end

    # Decide node border colors
    node_border_colors_gm = Array{String,1}(N)
    for i=1:N
        if (i % (m+2) == 1) || (i % (m+2) == 0) || (i <= m+2) || (i >= (m+2)*(n+1))
            node_border_colors_gm[i] = "white"
        else
            node_border_colors_gm[i] = "black"
        end
    end

    h = Junet.plot(gm, size=(300,300), layout=(x,y), node_color="white", node_border_color=node_border_colors_gm,
        node_border_width=2, edge_color=[edge_color_gm(i, adj_matrix_junet, α, β) for i = 1:edgecount(gm)],
        edge_width=3, format=:eps, file=gm_file)
end

# Function to draw a periodic lattice using GraphPlot.jl
function draw_periodic_lattice2(adj_matrix::Array, m::Integer, n::Integer, α::Real, β::Real)
    # Get the graph that Junet can plot
    adj_matrix_plot, struct_adj_matrix_plot = non_periodic_lattice(adj_matrix,m+2,n+2)

    # Layout
    N = (m+2)*(n+2)
    x = Array{Float64,1}(undef,N)
    y = Array{Float64,1}(undef,N)

    node_map_gm = zeros(m+2,n+2)
    node_map_gm[1:N] = 1:N

    for i=1:N
        y[i], x[i] = Tuple(findfirst(isequal(i),node_map_gm))
        if i in node_map_gm[:,1]
            x[i] += 0.3
        elseif i in node_map_gm[:,n+2]
            x[i] += -0.3
        elseif i in node_map_gm[1,:]
            y[i] += +0.3
        elseif i in node_map_gm[m+2,:]
            y[i] += -0.3
        end
    end

    # Decide node border colors
    node_border_colors_gm = Array{String,1}(undef,N)
    for i=1:N
        if (i % (m+2) == 1) || (i % (m+2) == 0) || (i <= m+2) || (i >= (m+2)*(n+1))
            node_border_colors_gm[i] = "white"
        else
            node_border_colors_gm[i] = "black"
        end
    end

    gm = SimpleGraph(struct_adj_matrix_plot)
    gplot(gm, x, y)
end

# decide the edge color to be used when plotting a graphical model
function edge_color_gm(ind_edge::Integer, adj_matrix::Array, α::Real, β::Real)
    A = sparse(adj_matrix)
    R, C, V = findnz(A)
    v = V[ind_edge]

    if v ≈ α
        edge_color="orange"
    elseif v ≈ β
        edge_color="red"
    elseif v ≈ -α
        edge_color="cyan"
    elseif v ≈ -β
        edge_color="blue"
    end

    return edge_color
end

# Assert if the correct gm has been obtained
function assert_correct_gm_old(learned_gm::Array, true_gm::Array, τ::Real)
    learned_struct_gm = zeros(size(learned_gm))

    for i = 1:size(learned_gm,1)
        for j = 1:size(learned_gm,2)
            if abs(learned_gm[i,j]) <= τ
                learned_struct_gm[i,j] = 0
            else
                learned_struct_gm[i,j] = sign(learned_gm[i,j])
            end
        end
    end

    true_struct_gm = sign.(true_gm)

    return all(learned_struct_gm .≈ true_struct_gm), learned_struct_gm, true_struct_gm
    #return all(learned_struct_gm .≈ true_struct_gm)
end

# Assert if the correct gm has been obtained
function assert_correct_gm(learned_gm::Array, true_gm::Array, τ::Real)
    learned_struct_gm = Array{Int8,2}(undef,size(learned_gm))
    true_struct_gm = Array{Int8,2}(undef,size(true_gm))

    for i = 1:size(learned_gm,1)
        for j = 1:size(learned_gm,2)
            if abs(learned_gm[i,j]) < τ
                learned_struct_gm[i,j] = 0
            else
                learned_struct_gm[i,j] = sign(learned_gm[i,j])
            end
        end
    end

    true_struct_gm = sign.(true_gm)

    #return all(learned_struct_gm .== true_struct_gm), learned_struct_gm, true_struct_gm
    return all(learned_struct_gm .== true_struct_gm)
end

# Assert if the correct gm has been obtained
function assert_correct_gm2(learned_gm::Array, true_gm::Array, τ::Real)
    learned_struct_gm = Array{Int8,2}(size(learned_gm))
    true_struct_gm = Array{Int8,2}(size(true_gm))

    for i = 1:size(learned_gm,1)
        for j = 1:size(learned_gm,2)
            if abs(learned_gm[i,j]) < τ
                learned_struct_gm[i,j] = 0
            else
                learned_struct_gm[i,j] = sign(learned_gm[i,j])
            end
        end
    end

    true_struct_gm = sign.(true_gm)

    return all(learned_struct_gm .== true_struct_gm), learned_struct_gm, true_struct_gm
end

# get the optimal number of samples for a given graphical model -- Glauber Dynamics
function get_M_opt(true_adj_matrix::Array, learning_method::GMLFormulation, learning_algo::GMLMethod, τ=0.2, L_success=45, M_guess=1000, M_factor=0.1)
    # Graphical model
    true_gm = FactorGraph(true_adj_matrix)

    # Initialize
    N_trials = 0    # Number of attempts
    L_trials = 0    # Number of successful trials so far
    M = copy(M_guess)   # Number of samples
    FLAG_correct_gm = false # FLAG if correct gm is recovered
    samples = sample(true_gm, M)

    while L_success > L_trials
        # Learn the graphical model
        learned_adj_matrix = learn(samples, learning_method, learning_algo)

        # Assert if correct GM
        FLAG_correct_gm = assert_correct_gm(learned_adj_matrix, true_adj_matrix, τ)

        if FLAG_correct_gm
            N_trials += 1
            L_trials += 1
            samples = sample(true_gm, M)
        else
            N_trials = 0
            L_trials = 0
            M = Int(floor(((1 + M_factor)*M)/100))*100
            samples = sample(true_gm, M)
        end
        @printf("FLAG=%d, M=%d, n_trials=%d, P=%d/%d\n", FLAG_correct_gm, M, N_trials, L_trials, L_success)
    end

    return M
end

# get the optimal number of samples for a given graphical model -- Glauber Dynamics
function get_M_opt_glauber_dynamics(true_adj_matrix::Array, learning_method::GMLFormulation, learning_algo::GMLMethod, sampling_regime::SamplingRegime, τ=0.2, L_success=45, M_guess=1000, M_factor=0.1)
    # Graphical model
    true_gm = FactorGraph(true_adj_matrix)

    # Initialize
    N_trials = 0    # Number of attempts
    L_trials = 0    # Number of successful trials so far
    M = copy(M_guess)   # Number of samples
    FLAG_correct_gm = false # FLAG if correct gm is recovered
    samples_T, samples_mixed = gibbs_sampling2(true_gm, M, sampling_regime)

    while L_success > L_trials
        # Learn the graphical model
        learned_adj_matrix = learn_glauber_dynamics(samples_T, learning_method, learning_algo)

        # Assert if correct GM
        FLAG_correct_gm = assert_correct_gm(learned_adj_matrix, true_adj_matrix, τ)

        if FLAG_correct_gm
            N_trials += 1
            L_trials += 1
            samples_T, samples_mixed = gibbs_sampling2(true_gm, M, sampling_regime)
        else
            N_trials = 0
            L_trials = 0
            M = Int(floor(((1 + M_factor)*M)/50))*50
            samples_T, samples_mixed = gibbs_sampling2(true_gm, M, sampling_regime)
        end
        @printf("FLAG=%d, M=%d, n_trials=%d, P=%d/%d\n", FLAG_correct_gm, M, N_trials, L_trials, L_success)
    end

    return M
end

# get the optimal number of samples for a given graphical model -- Glauber Dynamics
function get_M_opt_glauber_dynamics_regularization(true_adj_matrix::Array, learning_method::GMLFormulation, learning_algo::GMLMethod, sampling_regime::SamplingRegime, τ=0.2, L_success=45, M_guess=1000, M_factor=0.1)
    # Graphical model
    true_gm = FactorGraph(true_adj_matrix)

    # Initialize
    N_trials = 0    # Number of attempts
    L_trials = 0    # Number of successful trials so far
    M = copy(M_guess)   # Number of samples
    FLAG_correct_gm = false # FLAG if correct gm is recovered
    samples_T, samples_mixed = gibbs_sampling2(true_gm, M, sampling_regime)

    while L_success > L_trials
        # Learn the graphical model
        learned_adj_matrix = GraphicalModelLearning.learn_glauber_dynamics_regularization(samples_T, learning_method, learning_algo)

        # Assert if correct GM
        FLAG_correct_gm = assert_correct_gm(learned_adj_matrix, true_adj_matrix, τ)

        if FLAG_correct_gm
            N_trials += 1
            L_trials += 1
            samples_T, samples_mixed = gibbs_sampling2(true_gm, M, sampling_regime)
        else
            N_trials = 0
            L_trials = 0
            M = Int(floor(((1 + M_factor)*M)/50))*50
            samples_T, samples_mixed = gibbs_sampling2(true_gm, M, sampling_regime)
        end
        @printf("FLAG=%d, M=%d, n_trials=%d, P=%d/%d\n", FLAG_correct_gm, M, N_trials, L_trials, L_success)
    end

    return M
end

# Function to carry out one active learning return
function active_learning_run(m0::Integer, mtot::Integer, mbatch::Integer, num_spins::Integer, true_adj_matrix::Array, learning_method::GMLFormulation, learning_algo::GMLMethod, τ=0.2, FLAG_verbose=true)
    #=
    Inputs:
    m0 -  size of the first batch
    mbatch - preferred size of samples in each batch
    mtot - total number of samples
    true_adj_matrix is basically my oracle here

    Assuming the sampling regime is M-regime and hence not mentioned here
    =#
    # Mon oracle
    true_gm = FactorGraph(true_adj_matrix)

    # Create alphabet
    config_number = 2^num_spins
    X_U = [config for config=0:(config_number -1)]
    p_U = (1/config_number)*ones(config_number)

    # Get the initial set of samples
    samples_M, samples_mixed = gibbs_sampling_query(true_gm, m0, X_U, p_U, M_regime())
    m_remaining = mtot-m0
    M = copy(m0)

    # learn the graphical model based on the samples
    learned_adj_matrix = learn_glauber_dynamics(samples_M, RISE(), NLP())

    if FLAG_verbose
        FLAG_correct_gm = assert_correct_gm(learned_adj_matrix, true_adj_matrix, τ)
        @printf("Active: FLAG=%d, M=%d \n", FLAG_correct_gm, M)
    end

    while m_remaining > 0
        # create the entropy distribution
        s_configs_learned = entropy_configs(learned_adj_matrix,X_U,n)
        q_s_learned = s_configs_learned/sum(s_configs_learned)

        # modify the distribution slightly by adding a small amount of uniform distribution
        lambda = 1 - 1/((M)^(1/6))
        q_s = lambda*q_s_learned + (1-lambda)*p_U

        # get a batch of samples and update remaining number of queries
        mtemp = min(m_remaining,mbatch)
        samples_M_temp, samples_mixed_temp = gibbs_sampling_query(true_gm, mtemp, X_U, q_s, M_regime(), false)

        m_remaining = m_remaining - mtemp
        M += mtemp

        # add to the existing set of samples
        samples_M = GraphicalModelLearning.add_histograms(samples_M,samples_M_temp)

        # learn the graphical model based on the samples
        learned_adj_matrix = learn_glauber_dynamics(samples_M, RISE(), NLP())

        if FLAG_verbose
            FLAG_correct_gm = assert_correct_gm(learned_adj_matrix, true_adj_matrix, τ)
            @printf("Active: FLAG=%d, M=%d \n", FLAG_correct_gm, M)
        end
    end

    return learned_adj_matrix
end

# get the optimal number of samples for a given graphical model -- Glauber Dynamics
function get_M_opt_glauber_dynamics_AL(true_adj_matrix::Array, τ=0.2, L_success=45, M_guess=1000, M_factor=0.1)
    # Graphical model
    true_gm = FactorGraph(true_adj_matrix)
    n = true_gm.varible_count

    # Initialize
    N_trials = 0    # Number of attempts
    L_trials = 0    # Number of successful trials so far
    M = copy(M_guess)   # Number of samples
    FLAG_correct_gm = false # FLAG if correct gm is recovered

    learned_adj_matrix = Array{Float64,2}(undef,size(true_adj_matrix))

    # Parameters of run
    #m0 = min(Int(floor(M/2)), 5000)
    m0 = Int(floor(M/3))
    mbatch = Int(floor((M-m0)/15))

    while L_success > L_trials
        # Carry out an independent active learning run
        learned_adj_matrix = active_learning_run(m0, M, mbatch, n, true_adj_matrix, RISE(), NLP(), τ, false)

        # Assert if correct GM
        FLAG_correct_gm = assert_correct_gm(learned_adj_matrix, true_adj_matrix, τ)

        if FLAG_correct_gm
            N_trials += 1
            L_trials += 1
        else
            N_trials = 0
            L_trials = 0
            M = Int(floor(((1 + M_factor)*M)/50))*50

            m0 = Int(floor(M/3))
            mbatch = Int(floor((M-m0)/15))
        end
        @printf("FLAG=%d, M=%d, n_trials=%d, P=%d/%d\n", FLAG_correct_gm, M, N_trials, L_trials, L_success)
    end

    return M
end

# Evaluate the gradient of the objective function afterwards
function get_gradient_value(samples::Array, learned_gm::Array, learning_method::GMLFormulation, learning_method_name::String)
    num_conf, num_spins, num_samples = data_info(samples)

    grad_value = 0
    grad_value_i_array = Array{Float64,1}(num_spins)

    for current_spin = 1:num_spins
        nodal_stat  = [ samples[k, 1 + num_spins + current_spin] * (i == current_spin ? 1 : samples[k, 1 + i]) for k=1:num_conf_i , i=1:num_spins]

        x = learned_gm[current_spin,:]

        if learning_method_name == "RISE"
            iso_obj_deriv(x,k) = (1/num_samples)*( sum(samples[:,1].*(exp.(-1*(nodal_stat*x))).*(-1*nodal_stat[:,k])) )
        elseif learning_method_name == "RPLE"
            iso_obj_deriv(x,k) = (1/num_samples)*( sum(samples[:,1].*(1 ./ (1 + exp.(2*(nodal_stat*x)))).*(-2*nodal_stat[:,k])) )
        end

        for ind_spin = 1:num_spins
            grad_value_i_array[ind_spin] = iso_obj_deriv(x,ind_spin)
        end
        grad_value += norm(grad_value_i_array,Inf)
    end
    return grad_value/num_spins
end

# Evaluate the gradient of the objective function afterwards
function get_gradient_value_glauber_dynamics(samples::Array, learned_gm::Array, learning_method::GMLFormulation, learning_method_name::String)
    num_conf, num_spins, num_samples = data_info_glauber_dynamics(samples)

    grad_value = 0
    grad_value_i_array = Array{Float64,1}(num_spins)

    for current_spin = 1:num_spins
        samples_i = samples[find(isequal(current_spin),samples[:,2]),:]
        num_conf_i = size(samples_i,1)
        num_samples_i = sum(samples_i[:,1])

        nodal_stat  = [ samples_i[k, 2 + num_spins + current_spin] * (i == current_spin ? 1 : samples_i[k, 2 + i]) for k=1:num_conf_i , i=1:num_spins]

        x = learned_gm[current_spin,:]

        if learning_method_name == "RISE"
            iso_obj_deriv(x,k) = (1/num_samples_i)*( sum(samples_i[:,1].*(exp.(-1*(nodal_stat*x))).*(-1*nodal_stat[:,k])) )
        elseif learning_method_name == "RPLE"
            iso_obj_deriv(x,k) = (1/num_samples_i)*( sum(samples_i[:,1].*(1 ./ (1 + exp.(2*(nodal_stat*x)))).*(-2*nodal_stat[:,k])) )
        end

        for ind_spin = 1:num_spins
            grad_value_i_array[ind_spin] = iso_obj_deriv(x,ind_spin)
        end
        grad_value += norm(grad_value_i_array,Inf)
    end
    return grad_value/num_spins
end

# Evaluate the gradient of the objective function afterwards
function get_hessian_value_glauber_dynamics(samples::Array, learned_gm::Array, learning_method::GMLFormulation, learning_method_name::String)
    num_conf, num_spins, num_samples = data_info_glauber_dynamics(samples)

    grad_value = 0
    hessian_matrix = Array{Float64,1}(num_spins-1,num_spins-1)

    current_spin = rand(1:num_spins)
    samples_i = samples[find(isequal(current_spin),samples[:,2]),:]
    num_conf_i = size(samples_i,1)
    num_samples_i = sum(samples_i[:,1])

    nodal_stat  = [ samples_i[k, 2 + num_spins + current_spin] * (i == current_spin ? 1 : samples_i[k, 2 + i]) for k=1:num_conf_i , i=1:num_spins]
    cov_stat = [ samples_i[k, 2 + num_spins + current_spin] * (i == current_spin ? 1 : samples_i[k, 2 + i]) for k=1:num_conf_i , i=1:num_spins]

    x = learned_gm[current_spin,:]

    if learning_method_name == "RISE"
        iso_obj_deriv(x,k,l) = (1/num_samples_i)*( sum(samples_i[:,1].*(exp.(-1*(nodal_stat*x))).*(-1*nodal_stat[:,k])) )
    elseif learning_method_name == "RPLE"
        iso_obj_deriv(x,k,l) = (1/num_samples_i)*( sum(samples_i[:,1].*(1 ./ (1 + exp.(2*(nodal_stat*x)))).*(-2*nodal_stat[:,k])) )
    end

    for ind_spin = 1:num_spins
        grad_value_i_array[ind_spin] = iso_obj_deriv(x,ind_spin)
    end

    return grad_value/num_spins
end


text:
struct TransDateTimeVariate <: AbstractDateTimeVariate
    name::String
    basevariate::AbstractDateTimeVariate
    transform::Function
end

Base.length(var::TransDateTimeVariate) = length(var.basevariate)

function slice(dtvariate::TransDateTimeVariate, fromobs::Integer, toobs::Integer, slicelength::Integer)
    basevar = dtvariate.basevariate
    f = dtvariate.transform
    g = (ms::Int64) -> begin
        dt = Dates.epochms2datetime(ms)
        Dates.datetime2epochms(f(dt))
    end
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slices = slice(basevar, fromobs, toobs, slicelength)
    mapslice(g, slices, slicelength, Int64)
end

function Base.map(dtvariate::TransDateTimeVariate, dataframe::AbstractDataFrame)
    TransDateTimeVariate(dtvariate.name, map(dtvariate.basevariate, dataframe), dtvariate.transform)
end

function datetimevariate(f, basevariate::AbstractDateTimeVariate)
    TransDateTimeVariate("$(f)($(getname(basevariate)))", basevariate, f)
end

# function Base.broadcast(f, basevariate::AbstractDateTimeVariate)
#     if typeof(f(zero(Int64))) <: AbstractFloat
#         TransCovariate("$(f)($(getname(basecovariate)))", basecovariate, f)
#     else
#         TransCovBoolVariate{S}("$(f)($(getname(basecovariate)))", basecovariate, f)
#     end
# end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
module Swagger

using HTTP
using JSON
using MbedTLS
using Dates

import Base: convert, show, summary, getindex, keys, length, getproperty, setproperty!, propertynames
import JSON: lower
import Base: iterate

abstract type SwaggerModel end
abstract type SwaggerApi end

include("client.jl")
include("json.jl")
include("val.jl")

end # module Swagger


text:
using Printf

function subroutine(number)
        @printf("Αυτή είναι μια υπορουτίνα με αριθμό %d.\n", number)
end

println("Γεια σου, κόσμε!")

if (1 == 1)
    println("Εκτέλεση if: Αληθές και όχι ψευδές!")
end

print("Εκτέλεση βρόγχου: ")

for i = 0:9
    @printf("%d ", i)
end

println("")

subroutine(5)

max_stars_repo_name: ubuntu-gr/OlesOiGlosses

max_issues_repo_name: ubuntu-gr/OlesOiGlosses

max_forks_repo_name: ubuntu-gr/OlesOiGlosses