text:
const INPUTS = Dict{Symbol, Any}(
    :auth                    => ("BUMPSTDLIBS_TOKEN",                   nothing),
    :close_old_pull_requests => ("BUMPSTDLIBS_CLOSE_OLD_PULL_REQUESTS", "true"),
    :push_if_no_changes      => ("BUMPSTDLIBS_PUSH_IF_NO_CHANGES",      "false"),
    :stdlibs_to_include      => ("BUMPSTDLIBS_STDLIBS_TO_INCLUDE",      "all"),
)

function post_process_input(::Val, value)
    return value
end

function post_process_input(::Val{:auth}, value)
    return GitHub.authenticate(value)
end

function post_process_input(::Val{:close_old_pull_requests}, value)
    return parse(Bool, value)
end

function post_process_input(::Val{:push_if_no_changes}, value)
    return parse(Bool, value)
end

function get_input_from_environment(input::Symbol,
                                    env_var_name::Union{AbstractString, Nothing} = nothing)
    if haskey(INPUTS, input)
        default_env_var_name, default_value = INPUTS[input]
        if env_var_name isa Nothing
            env_var_name_to_use = default_env_var_name
        else
            env_var_name_to_use = env_var_name
        end
        env_var_contents = strip(get(ENV, env_var_name_to_use, ""))
        if isempty(env_var_contents)
            if default_value isa Nothing
                throw(ArgumentError("Either the `$(_env_var_name)` environment variable is undefined, or it is defined but empty."))
            else
                return post_process_input(Val(input), default_value)
            end
        else
            return post_process_input(Val(input), env_var_contents)
        end
    end
    throw(ArgumentError("$(input) is not a valid input"))
end


text:
module test_core

include("../../src/HijriConverter.jl")
using Test
using Dates
using .HijriConverter

hijri_obj_1 = HijriConverter.Hijri(1410, 8, 13)
hijri_obj_2 = HijriConverter.Hijri(1460, 4, 21)
hijri_obj_3 = HijriConverter.Hijri(1343, 1, 1)
hijri_obj_4 = HijriConverter.Hijri(1500, 12, 30)

function TestHijri()
#        @test hijri_obj.__repr__() == "Hijri(1410, 8, 13)"

#       @test hijri_obj.__str__() == "1410-08-13"

#        @test self.hijri_obj.__hash__() == hash(("Hijri", 1410, 8, 13))

#    @pytest.mark.parametrize("test_input", ["__gt__", "__ge__", "__lt__", "__le__"])
 #   def test_comparison_notimplemented(self, test_input):
  #      assert getattr(self.hijri_obj, test_input)("1410-08-13") == NotImplemented
end

function test_equality()
        @test hijri_obj_1 == HijriConverter.Hijri(1410, 8, 13)
        @test hijri_obj_1 != HijriConverter.Hijri(1410, 8, 14)
        @test hijri_obj_1 != "1410-08-13"
end

function test_ordering() #TODO: implement comparison
        @test hijri_obj_1 > HijriConverter.Hijri(1410, 8, 12)
        @test hijri_obj_1 >= HijriConverter.Hijri(1410, 8, 13)
        @test hijri_obj_1 < HijriConverter.Hijri(1410, 8, 14)
        @test hijri_obj_1 <= HijriConverter.Hijri(1410, 8, 13)
end

function test_fromisoformat()
        @test HijriConverter.fromisoformat("1410-08-13") == hijri_obj_1
        @test HijriConverter.fromisoformat("1460-04-21") == hijri_obj_2
        @test HijriConverter.fromisoformat("1343-01-01") == hijri_obj_3
        @test HijriConverter.fromisoformat("1500-12-30") == hijri_obj_4
end

function test_today()
        @test HijriConverter.to_gregorian(HijriConverter.hijri_today()) == Dates.today()
end

function test_year()
        @test hijri_obj_1.year == 1410
        @test hijri_obj_2.year == 1460
        @test hijri_obj_3.year == 1343
        @test hijri_obj_4.year == 1500
end

function test_month()
        @test hijri_obj_1.month == 8
        @test hijri_obj_2.month == 4
        @test hijri_obj_3.month == 1
        @test hijri_obj_4.month == 12
end

function test_day()
        @test hijri_obj_1.day == 13
        @test hijri_obj_2.day == 21
        @test hijri_obj_3.day == 1
        @test hijri_obj_4.day == 30
end

function test_datetuple()
        @test HijriConverter.datetuple(hijri_obj_1) == (1410, 8, 13)
end

function test_isoformat()
        @test HijriConverter.isoformat(hijri_obj_1) == "1410-08-13"
end

function test_dmyformat()
        @test HijriConverter.dmyformat(hijri_obj_1) == "13/08/1410"
        @test HijriConverter.dmyformat(hijri_obj_1; padding=false) == "13/8/1410"
        @test HijriConverter.dmyformat(hijri_obj_1; separator=".") == "13.08.1410"
end

function test_month_length()
        @test HijriConverter.month_length(hijri_obj_1) == 29
end

function test_month_name()
        @test HijriConverter.month_name(hijri_obj_1) == "Sha’ban"
        @test HijriConverter.month_name(hijri_obj_1, "en") == "Sha’ban"
        @test HijriConverter.month_name(hijri_obj_1, "en-US") == "Sha’ban"

        @test HijriConverter.month_name(hijri_obj_1, "ar") == "شعبان"
        @test HijriConverter.month_name(hijri_obj_1, "az") == "şaban"
        @test HijriConverter.month_name(hijri_obj_1, "bn") == "শাবান"
        @test HijriConverter.month_name(hijri_obj_1, "tr") == "Şaban"
end

function test_weekday()
        @test HijriConverter.weekday(hijri_obj_1) == 5
end

function test_iso_weekday()
        @test HijriConverter.isoweekday(hijri_obj_1) == 6
end

function test_day_name()
        @test HijriConverter.day_name(hijri_obj_1) == "Saturday"
        @test HijriConverter.day_name(hijri_obj_1, "en") == "Saturday"
        @test HijriConverter.day_name(hijri_obj_1, "en-US") == "Saturday"
end

function test_notation()
        @test HijriConverter.notation(hijri_obj_1) == "AH"
        @test HijriConverter.notation(hijri_obj_1, "en") == "AH"
        @test HijriConverter.notation(hijri_obj_1, "en-US") == "AH"
end

function test_to_julian()
        @test HijriConverter.to_julian(hijri_obj_1) == 2447961
end

function test_to_gregorian()
        @test HijriConverter.datetuple(HijriConverter.to_gregorian(hijri_obj_1)) == (1990, 3, 10)
        @test HijriConverter.datetuple(HijriConverter.to_gregorian(hijri_obj_2)) == (2038, 5, 26)
        @test HijriConverter.datetuple(HijriConverter.to_gregorian(hijri_obj_3)) == (1924, 8, 2)
        @test HijriConverter.datetuple(HijriConverter.to_gregorian(hijri_obj_4)) == (2077, 11, 16)
end

function test_month_index()
        @test HijriConverter._month_index(hijri_obj_1) == 811
end
end

ext: jl

max_stars_repo_name: ceferisbarov/hijri-converter

max_issues_repo_name: ceferisbarov/hijri-converter


text:
function jarnik(g::AbstractGraph, start::Vertex)
    # clear vertex marks
    g.mark = g.mark * 0

    # init minimum spanning tree, distance and vertex vectors
    mst = typeof(g)(countvertices(g))
    vertices = fill(start, countvertices(g))
    distance = fill(Inf, countvertices(g))
    distance[start] = 0
    for i = 1:countvertices(g)
        v = minimumvertex!(g, distance)
        if (v != start)
            setedge!(mst, vertices[v], v, getweight(g, vertices[v], v))
            setedge!(mst, v, vertices[v], getweight(g, v, vertices[v]))
            # println((vertices[v], v), " ", getweight(g, vertices[v], v))
        end
        (distance[v] == Inf) && throw(ErrorException("Unreachable vertices."))
        setmark!(g, v, 1)
        let w = first(g, v)
            while w <= countvertices(g)
                if distance[w] > getweight(g, v, w)
                    distance[w] = getweight(g, v, w)
                    vertices[w] = v
                end
                w = next(g, v, w)
            end
        end
    end
    # clear vertex marks
    g.mark = g.mark * 0
    return mst
end


text:
#############################################################################
# models.jl
# Handles the models of problems
#
# a PrimalModel has the form:
#   minimize g(x) + f(x)
#   s.t.     x \in s, s in ContiSet
#            x \in C, C is convex,
# where g(x) is a convex function, and f(x) is the Lovasz extension of a submodular function
#
# a DualModel has the form:
#  minimize g(x)
#  s.t.      x \in P, P is a polyhedron associated with a submodular function
##            x \in C, C is convex
#############################################################################

export PrimalModel
export ConvexProblem, ConvexLovasz, ConvexLovaszAbs
export DualModel
export AssocPolyConstrained, LPoverAssocPoly
export get_model

### models of the primal problems
abstract type PrimalModel <: SCOPEModel end

### convex problem
mutable struct ConvexProblem <: PrimalModel end
### Convex function + Lovasz extension
mutable struct ConvexLovasz <: PrimalModel
  convex_part::Problem
  lovasz::LovaszExtAtom
end
### Convex function + Lovasz extension on absolute values
mutable struct ConvexLovaszAbs <: PrimalModel
  convex_part::Problem
  lovaszabs::LovaszExtAbsAtom
end

### models of the dual problems
abstract type DualModel <: SCOPEModel end

### convex optimization over a polyhedron associated with a submodular function
mutable struct AssocPolyConstrained <: DualModel
  prob::Problem
  poly::AssocPoly
end
### linear programming over associated polyhedra
mutable struct LPoverAssocPoly <: DualModel end

function get_model(objective::AbstractExpr, constraints::AbstractArray=Constraint[])
  objective_cv = get_cv(objective)
  objective_sv = get_sv(objective)
  constraints_cv = []
  constraints_sv = []
  for i = 1:length(constraints)
    constraints_cv = vcat(constraints_cv, get_cv(constraints[i].lhs), get_cv(constraints[i].rhs))
    constraints_sv = vcat(constraints_sv, get_sv(constraints[i].lhs), get_sv(constraints[i].rhs))
  end
  if objective.head == :+
    lonum = 0                                         # number of lovasz extensions
    loind = 1                                         # index of lovasz extensions
    lanum = 0                                         # number of lovasz + abs
    laind = 1                                         # index of lovasz + abs
    childnum = length(objective.children)
    for i = 1:childnum
      if objective.children[i].head == :lovasz
        lonum += 1
        loind = i
      elseif objective.children[i].head == :lovaszabs
        lanum += 1
        laind = i
      end
    end
    if (lonum + lanum) > 1
      error("Cannot solve objective functions with multiple Lovasz extensions!")
    elseif lonum == 1
      lovasz = objective.children[loind]
      if loind == 1
        objective1 = objective.children[2]
        if childnum > 2
          for i = 3 : childnum
            objective1 += objective.children[i]
          end
        end
      else
        objective1 = objective.children[1]
        for i = 2:(loind - 1)
          objective1 += objective.children[i]
        end
        for i = (loind + 1) : childnum
          objective1 += objective.children[i]
        end
      end
      # convex + Lovasz extension
      if vexity(objective1) == ConvexVexity() || vexity(objective1) == AffineVexity()
        q = Problem(:minimize, objective1, constraints)
        return ConvexLovasz(q, lovasz)
      else
        warn("The objective function is the sum of the Lovasz extension of a submodular function and a non-DCP-convex funcion; solution may be inaccurate.")
        q = Problem(:minimize, objective1, constraints)
        return ConvexLovasz(q, lovasz)
      end
    elseif lanum == 1
      lovasz = objective.children[laind]
      if laind == 1
        objective1 = objective.children[2]
        if childnum > 2
          for i = 3 : childnum
            objective1 += objective.children[i]
          end
        end
      else
        objective1 = objective.children[1]
        for i = 2:(laind - 1)
          objective1 += objective.children[i]
        end
        for i = (laind + 1) : childnum
          objective1 += objective.children[i]
        end
      end
      # convex + Lovasz extension on absolute values
      if vexity(objective1) == ConvexVexity() || vexity(objective1) == AffineVexity()
        q = Problem(:minimize, objective1, constraints)
        return ConvexLovaszAbs(q, lovasz)
      else
        warn("The objective function is the sum of the Lovasz extension of a submodular function and a non-DCP-convex funcion; solution may be inaccurate.")
        q = Problem(:minimize, objective1, constraints)
        return ConvexLovaszAbs(q, lovasz)
      end
    else
      return ConvexProblem()
    end
  elseif length(constraints) == 1
    if isa(constraints[1], SetConstraint)             # set-constrained optimization problems
      if isa(constraints[1].rhs, AssocPoly) && length(objective_sv) == 0 && length(objective_cv) == 1 && objective_cv[1].id_hash == constraints_cv[1].id_hash
        if vexity(objective) == AffineVexity()
          return LPoverAssocPoly()
        elseif vexity(objective) == ConvexVexity()
          return AssocPolyConstrained(Problem(:minimize, objective), constraints[1].rhs)
        end
      end
    else
      return ConvexProblem()
    end
  else
    return ConvexProblem()                            # fail-safe option, TODO: add other one
  end
end


text:
using JuMP, Gurobi

## routing of restored wavelength
function WaveRerouting(OpticalTopo, failed_IPedge, failed_fibers_index, rerouting_K)
    nfailedge = size(failed_IPedge, 1)
    rehoused_IProutingEdge = []  # routing of rehoused IP link in terms of fiber link index
    rehoused_IProuting = []  # routing of rehoused IP link in terms of fiber node sequence
    failedIPbranckIndexAll = []  # index of rehoused IP link branches
    failedIPbrachIndexGroup = []  # index group of rehoused IP link branches

    ## an optical topology after fiber cut (an optical scenario)
    num_nodes = length(OpticalTopo["nodes"])
    optical_graph = LightGraphs.SimpleDiGraph(num_nodes)
    distances = Inf*ones(num_nodes, num_nodes)
    num_edges = length(OpticalTopo["links"])

    for i in 1:num_edges
        if i in failed_fibers_index  # delete the failed fiber
            continue
        else
            LightGraphs.add_edge!(optical_graph, OpticalTopo["links"][i][1], OpticalTopo["links"][i][2])
            distances[OpticalTopo["links"][i][1], OpticalTopo["links"][i][2]] = OpticalTopo["fiber_length"][i]
            distances[OpticalTopo["links"][i][2], OpticalTopo["links"][i][1]] = OpticalTopo["fiber_length"][i]
        end
    end

    # find rerouting_K paths for each failed IP links, pay attention: IP links are bidirectional!
    global_IPbranch = 1
    for w in 1:nfailedge
        # println("failed IP src-dst: ", failed_IPedge[w][1], "-", failed_IPedge[w][2])
        state = LightGraphs.yen_k_shortest_paths(optical_graph, failed_IPedge[w][1], failed_IPedge[w][2], distances, rerouting_K)
        paths = state.paths
        # println("paths: ", paths)
        if length(paths) <= rerouting_K
            path_edges = []
            for p in 1:length(paths)
                k_path_edges = []
                for i in 1:length(paths[p])-1
                    e = findfirst(x -> x == (paths[p][i], paths[p][i+1]), OpticalTopo["links"])
                    append!(k_path_edges, e)
                end
                push!(path_edges, k_path_edges)
            end

            append!(rehoused_IProutingEdge, path_edges)
            append!(rehoused_IProuting, paths)
            append!(failedIPbranckIndexAll, range(global_IPbranch, length=length(paths), step=1))
            push!(failedIPbrachIndexGroup, range(global_IPbranch, length=length(paths), step=1))
            global_IPbranch += length(paths)
        end
    end

    return rehoused_IProutingEdge, rehoused_IProuting, failedIPbranckIndexAll, failedIPbrachIndexGroup
end

## wavelength assignment of restored wavelength, generating restoration link options that maximize restored capacity, this is an ILP
function RestoreILP(GRB_ENV, Fibers, FibercapacityCode, failedIPedges, failedIPBranchRoutingFiber, failedIPbranckIndexAll, failedIPbrachIndexGroup, failed_IP_initialbw, rerouting_K)
    println("solving restoration wavelength assignment ILP considering wavelength continuity")

    nFibers = length(Fibers)
    nfailedIPedges = length(failedIPedges)
    nfailedIPedgeBranchAll = length(failedIPbranckIndexAll)  # this number can be small than nfailedIPedges * nfailedIPedgeBranchPerLink
    nfailedIPedgeBranchPerLink = rerouting_K
    nwavelength = size(FibercapacityCode, 2)
    uni_failedIPedges = []
    reverse_failedIPedges = []
    for edge_index in 1:nfailedIPedges
        e = findfirst(x -> x == (failedIPedges[edge_index][2], failedIPedges[edge_index][1], failedIPedges[edge_index][3]), failedIPedges)
        if edge_index < e
            push!(uni_failedIPedges, edge_index)
            push!(reverse_failedIPedges, e)
        end
    end

    # how each IP branch is routed on fibers
    L = zeros(nfailedIPedgeBranchAll, nFibers)
    for t in 1:nfailedIPedgeBranchAll  # nfailedIPedgeBranchAll is global indexed
        for f in 1:nFibers
            if in(f, failedIPBranchRoutingFiber[t])
                L[t,f] = 1
            end
        end
    end

    if length(failedIPBranchRoutingFiber) > 0  # if this scenario has failures
        model = Model(() -> Gurobi.Optimizer(GRB_ENV))
        set_optimizer_attribute(model, "OutputFlag", 0)
        set_optimizer_attribute(model, "Threads", 32)

        @variable(model, restored_bw[1:nfailedIPedges] >= 0, Int)
        @variable(model, IPBranch_bw[1:nfailedIPedges, 1:nfailedIPedgeBranchPerLink] >= 0, Int)  # bandwidth allocation for all IP branches
        @variable(model, lambda[1:nfailedIPedgeBranchAll, 1:nFibers, 1:nwavelength] >=0, Bin)  # if IP link's branch use fiber and wavelength

        # Equation 14, wavelength resource used only once if the resource is usable
        for w in 1:nwavelength
            for f in 1:nFibers
                @constraint(model, sum(lambda[t,f,w] for t in 1:nfailedIPedgeBranchAll) <= FibercapacityCode[f,w])
            end
        end

        # Equation 15, translate wavelength usage to IPBranch_bw
        for l in 1:nfailedIPedges
            for t in 1:nfailedIPedgeBranchPerLink  # t is the index for branches of the failIP link, not global branch index
                if t <= length(failedIPbrachIndexGroup[l])
                    for f in 1:nFibers
                        @constraint(model, IPBranch_bw[l,t]*L[failedIPbrachIndexGroup[l][t],f] == sum(lambda[failedIPbrachIndexGroup[l][t],f,w] for w in 1:nwavelength))
                    end
                else
                    @constraint(model, IPBranch_bw[l,t] == 0)
                end
            end
        end

        # Equation 16, wavelength continuity
        for t in 1:nfailedIPedgeBranchAll
            for f in failedIPBranchRoutingFiber[t]
                for ff in failedIPBranchRoutingFiber[t]
                    for w in 1:nwavelength
                        @constraint(model, lambda[t,f,w]*L[t,f] == lambda[t,ff,w]*L[t,ff])
                    end
                end
            end
        end

        # Equation 17, restored bw should no larger than initial bw, 100 is per wavelength gbps
        for l in 1:nfailedIPedges
            @constraint(model, restored_bw[l]*100 <= failed_IP_initialbw[l])
            @constraint(model, restored_bw[l] == sum(IPBranch_bw[l,t] for t in 1:nfailedIPedgeBranchPerLink))
        end

        # Auxiliary: bidirectional link bandwidth equal
        for e in 1:length(uni_failedIPedges)
            @constraint(model, restored_bw[uni_failedIPedges[e]] == restored_bw[reverse_failedIPedges[e]])
        end

        @objective(model, Max, sum(restored_bw[l] for l in 1:nfailedIPedges))  # maximizing total restorable bandwidth capacity
        optimize!(model)

        return value.(restored_bw), objective_value(model), value.(IPBranch_bw)

    else
        restored_bw = zeros(nfailedIPedges)
        IPBranch_bw = zeros(nfailedIPedges, nfailedIPedgeBranchPerLink)

        return restored_bw, 0, IPBranch_bw
    end
end

## wavelength assignment of restored wavelength, generating restoration link options that maximize restored capacity, this is a LP relaxed from previous ILP
function RestoreLP(GRB_ENV, Fibers, FibercapacityCode, failedIPedges, failedIPBranchRoutingFiber, failedIPbranckIndexAll, failedIPbrachIndexGroup, failed_IP_initialbw, rerouting_K)
    println("solving restoration wavelength assignment relaxed LP considering wavelength continuity")

    nFibers = length(Fibers)
    nfailedIPedges = length(failedIPedges)
    nfailedIPedgeBranchAll = length(failedIPbranckIndexAll)  # this number can be small than nfailedIPedges * nfailedIPedgeBranchPerLink
    nfailedIPedgeBranchPerLink = rerouting_K
    nwavelength = size(FibercapacityCode, 2)
    # println(nwavelength)
    uni_failedIPedges = []
    reverse_failedIPedges = []
    for edge_index in 1:nfailedIPedges
        e = findfirst(x -> x == (failedIPedges[edge_index][2], failedIPedges[edge_index][1], failedIPedges[edge_index][3]), failedIPedges)
        if edge_index < e
            push!(uni_failedIPedges, edge_index)
            push!(reverse_failedIPedges, e)
        end
    end

    # how each IP branch is routed on fibers
    L = zeros(nfailedIPedgeBranchAll, nFibers)
    for t in 1:nfailedIPedgeBranchAll  # nfailedIPedgeBranchAll is global indexed
        for f in 1:nFibers
            if in(f, failedIPBranchRoutingFiber[t])
                L[t,f] = 1
            end
        end
    end

    if length(failedIPBranchRoutingFiber) > 0  # if this scenario has failures
        model = Model(() -> Gurobi.Optimizer(GRB_ENV))
        set_optimizer_attribute(model, "OutputFlag", 0)
        set_optimizer_attribute(model, "Threads", 32)

        @variable(model, restored_bw[1:nfailedIPedges] >= 0)  # integer constraint relaxed
        @variable(model, IPBranch_bw[1:nfailedIPedges, 1:nfailedIPedgeBranchPerLink] >= 0)  # integer constraint relaxed, bandwidth allocation for all IP branches
        @variable(model, 0 <= lambda[1:nfailedIPedgeBranchAll, 1:nFibers, 1:nwavelength] <=1)  # integer constraint relaxed, if IP link's branch use fiber and wavelength

        # Equation 14, wavelength resource used only once if the resource is usable
        for w in 1:nwavelength
            for f in 1:nFibers
                @constraint(model, sum(lambda[t,f,w] for t in 1:nfailedIPedgeBranchAll) <= FibercapacityCode[f,w])
            end
        end

        # Equation 15, translate wavelength usage to IPBranch_bw
        for l in 1:nfailedIPedges
            for t in 1:nfailedIPedgeBranchPerLink  # t is the index for branches of the failIP link, not global branch index
                if t <= length(failedIPbrachIndexGroup[l])
                    for f in 1:nFibers
                        @constraint(model, IPBranch_bw[l,t]*L[failedIPbrachIndexGroup[l][t],f] == sum(lambda[failedIPbrachIndexGroup[l][t],f,w] for w in 1:nwavelength))
                    end
                else
                    @constraint(model, IPBranch_bw[l,t] == 0)
                end
            end
        end

        # Equation 16, wavelength continuity
        for t in 1:nfailedIPedgeBranchAll
            for f in failedIPBranchRoutingFiber[t]
                for ff in failedIPBranchRoutingFiber[t]
                    for w in 1:nwavelength
                        @constraint(model, lambda[t,f,w]*L[t,f] == lambda[t,ff,w]*L[t,ff])
                    end
                end
            end
        end

        # Equation 17, restored bw should no larger than initial bw, 100 is per wavelength gbps
        for l in 1:nfailedIPedges
            @constraint(model, restored_bw[l]*100 <= failed_IP_initialbw[l])
            @constraint(model, restored_bw[l] == sum(IPBranch_bw[l,t] for t in 1:nfailedIPedgeBranchPerLink))
        end

        # Auxiliary: bidirectional link bandwidth equal
        for e in 1:length(uni_failedIPedges)
            @constraint(model, restored_bw[uni_failedIPedges[e]] == restored_bw[reverse_failedIPedges[e]])
        end

        @objective(model, Max, sum(restored_bw[l] for l in 1:nfailedIPedges))  # maximizing total restorable bandwidth capacity
        optimize!(model)

        return value.(restored_bw), objective_value(model), value.(IPBranch_bw)

    else
        restored_bw = zeros(nfailedIPedges)
        IPBranch_bw = zeros(nfailedIPedges, nfailedIPedgeBranchPerLink)

        return restored_bw, 0, IPBranch_bw
    end
end

## check if a particular ticket satisfy the RWA constraints
function RestoreRWACheck(GRB_ENV, ticket_restored_bw, Fibers, FibercapacityCode, failedIPedges, failedIPBranchRoutingFiber, failedIPbranckIndexAll, failedIPbrachIndexGroup, failed_IP_initialbw, rerouting_K, verbose)
    if verbose println("[Checking] randomized rounding solution feasibility with optical-layer restoration ILP") end
    nFibers = length(Fibers)
    nfailedIPedges = length(failedIPedges)
    nfailedIPedgeBranchAll = length(failedIPbranckIndexAll)  # this number can be small than nfailedIPedges * nfailedIPedgeBranchPerLink
    nfailedIPedgeBranchPerLink = rerouting_K
    nwavelength = size(FibercapacityCode, 2)
    uni_failedIPedges = []
    reverse_failedIPedges = []
    for edge_index in 1:nfailedIPedges
        e = findfirst(x -> x == (failedIPedges[edge_index][2], failedIPedges[edge_index][1], failedIPedges[edge_index][3]), failedIPedges)
        if edge_index < e
            push!(uni_failedIPedges, edge_index)
            push!(reverse_failedIPedges, e)
        end
    end

    # how each IP branch is routed on fibers
    L = zeros(nfailedIPedgeBranchAll, nFibers)
    for t in 1:nfailedIPedgeBranchAll  # nfailedIPedgeBranchAll is global indexed
        for f in 1:nFibers
            if in(f, failedIPBranchRoutingFiber[t])
                L[t,f] = 1
            end
        end
    end

    if length(failedIPBranchRoutingFiber) > 0  # if this scenario has failures
        model = Model(() -> Gurobi.Optimizer(GRB_ENV))
        set_optimizer_attribute(model, "OutputFlag", 0)
        set_optimizer_attribute(model, "Threads", 32)

        @variable(model, restored_bw[1:nfailedIPedges] >= 0, Int)
        @variable(model, IPBranch_bw[1:nfailedIPedges, 1:nfailedIPedgeBranchPerLink] >= 0, Int)  # bandwidth allocation for all IP branches
        @variable(model, lambda[1:nfailedIPedgeBranchAll, 1:nFibers, 1:nwavelength] <= 0, Bin)  # if IP link's branch use fiber and wavelength

        # Equation 14, wavelength resource used only once if the resource is usable
        for w in 1:nwavelength
            for f in 1:nFibers
                @constraint(model, sum(lambda[t,f,w] for t in 1:nfailedIPedgeBranchAll) <= FibercapacityCode[f,w])
            end
        end

        # Equation 15, translate wavelength usage to IPBranch_bw
        for l in 1:nfailedIPedges
            for t in 1:nfailedIPedgeBranchPerLink  # t is the index for branches of the failIP link, not global branch index
                if t <= length(failedIPbrachIndexGroup[l])
                    for f in 1:nFibers
                        @constraint(model, IPBranch_bw[l,t]*L[failedIPbrachIndexGroup[l][t],f] == sum(lambda[failedIPbrachIndexGroup[l][t],f,w] for w in 1:nwavelength))
                    end
                else
                    @constraint(model, IPBranch_bw[l,t] == 0)
                end
            end
        end

        # Equation 16, wavelength continuity
        for t in 1:nfailedIPedgeBranchAll
            for f in failedIPBranchRoutingFiber[t]
                for ff in failedIPBranchRoutingFiber[t]
                    for w in 1:nwavelength
                        @constraint(model, lambda[t,f,w]*L[t,f] == lambda[t,ff,w]*L[t,ff])
                    end
                end
            end
        end

        # Equation 17, restored bw should no larger than initial bw, 100 is per wavelength gbps
        for l in 1:nfailedIPedges
            @constraint(model, restored_bw[l]*100 <= failed_IP_initialbw[l])
            @constraint(model, restored_bw[l] == sum(IPBranch_bw[l,t] for t in 1:nfailedIPedgeBranchPerLink))
        end

        # Auxiliary: bidirectional link bandwidth equal
        for e in 1:length(uni_failedIPedges)
            @constraint(model, restored_bw[uni_failedIPedges[e]] == restored_bw[reverse_failedIPedges[e]])
        end

        # do not solve the model, but just check the solution feasibility
        result_dict = Dict(restored_bw[1] => ticket_restored_bw[1])
        for l in 2:nfailedIPedges
            result_dict = merge!(result_dict, Dict(restored_bw[l] => ticket_restored_bw[l]))
        end
        # println("result_dict: ", result_dict)
        feasibility = primal_feasibility_report(model, result_dict, skip_missing = true)
        # println("feasibility: ", feasibility)

        if length(feasibility) == 0
            if verbose
                printstyled(" - $(ticket_restored_bw) feasibility true - ", color=:green)
                println(feasibility)
            end
            return true
        else
            if verbose
                printstyled(" - $(ticket_restored_bw) feasibility false - ", color=:yellow)
                println(feasibility)
            end
            return false
        end

    else
        return true  # non failure scenario
    end
end

## randomized rounding based on relaxed LP solution
function RandomRounding(GRB_ENV, LP_restored_bw, restored_bw_rwa, failedIPedges, failed_IP_initialbw, ticket_set_size, option_gap, OpticalTopo, rehoused_IProutingEdge, failedIPbranckindex, failedIPbrachGroup, optical_rerouting_K, verbose)
    # println("Randomized rounding gap: ", option_gap)
    ilp_check = 1
    initial_wavelength_num = failed_IP_initialbw ./ 100  # 1 wavelength = 100 Gbps
    probabilities_up = []
    nrestored = size(LP_restored_bw, 1)
    nfailedIPedges = length(failedIPedges)
    # identify bidirectional IP relationship
    uni_failedIPedges = []
    reverse_failedIPedges = []
    for edge_index in 1:nfailedIPedges
        e = findfirst(x -> x == (failedIPedges[edge_index][2], failedIPedges[edge_index][1], failedIPedges[edge_index][3]), failedIPedges)
        if edge_index < e
            push!(uni_failedIPedges, edge_index)
            push!(reverse_failedIPedges, e)
        end
    end

    for i in 1:length(uni_failedIPedges)
        push!(probabilities_up, LP_restored_bw[uni_failedIPedges[i]] - floor(Int, LP_restored_bw[uni_failedIPedges[i]]))
    end

    outer_iteration = min(ticket_set_size*nrestored*10) # a large number for generating options
    theoretical_max_ticket = floor(Int128, sqrt(prod(initial_wavelength_num)))
    if sqrt(prod(initial_wavelength_num)) - floor(sqrt(prod(initial_wavelength_num))) > 0.001 && verbose
        println("initial_wavelength_num $(initial_wavelength_num)")
    end
    if verbose println("The max number of possible ticket in this failure scenario: $(theoretical_max_ticket)") end
    restored_bw = zeros(outer_iteration, nrestored)
    real_restored_bw = zeros(ticket_set_size, nrestored)  # generated tickets for this scenario, restorable bw for each failed IP link

    # have the first randoized rounding result to be the planning (arrow naive) result
    rounded_index = 1
    for b in 1:length(uni_failedIPedges)
        real_restored_bw[rounded_index,uni_failedIPedges[b]] = restored_bw_rwa[uni_failedIPedges[b]]
        real_restored_bw[rounded_index,reverse_failedIPedges[b]] = restored_bw_rwa[reverse_failedIPedges[b]]
    end

    # have the second randomized rouding result to be the generic rounding results from LP
    if ticket_set_size >= 2
        rounded_index = 2
        for b in 1:length(uni_failedIPedges)
            real_restored_bw[rounded_index,uni_failedIPedges[b]] = floor(LP_restored_bw[uni_failedIPedges[b]])
            real_restored_bw[rounded_index,reverse_failedIPedges[b]] = floor(LP_restored_bw[reverse_failedIPedges[b]])
        end
    end

    # starting from the third ticket, randomized generation
    progress = ProgressMeter.Progress(min(ticket_set_size,theoretical_max_ticket), .1, "Running randomized rounding $(ticket_set_size) tickets (this scenario of $(nrestored/2) failed links has max $(outer_iteration) trials)...\n", 50)
    if ticket_set_size > 2
        for m in 1:outer_iteration
            # m'th rounding attempt to get restored_bw[outer_iteration, nrestored]
            for b in 1:length(uni_failedIPedges)  # for each failed link to restore
                rounding_range = max(initial_wavelength_num[uni_failedIPedges[b]], floor(LP_restored_bw[uni_failedIPedges[b]]))  # rouding stride
                # rounding_range = initial_wavelength_num[uni_failedIPedges[b]] - floor(LP_restored_bw[uni_failedIPedges[b]])  # rouding stride
                rd = rand(1) # generate a random number between 0 and 1
                # probability of the b'th failed link
                if probabilities_up[b] == 0  ## handling integer LP solution
                    if rd[1] < 0.3  # round up
                        restored_bw[m,uni_failedIPedges[b]] = min(initial_wavelength_num[uni_failedIPedges[b]], ceil(Int, LP_restored_bw[uni_failedIPedges[b]]) + rand(1:rounding_range))
                        restored_bw[m,reverse_failedIPedges[b]] = restored_bw[m,uni_failedIPedges[b]]
                    elseif rd[1] < 0.7  # stay still
                        restored_bw[m,uni_failedIPedges[b]] = LP_restored_bw[uni_failedIPedges[b]]
                        restored_bw[m,reverse_failedIPedges[b]] = restored_bw[m,uni_failedIPedges[b]]
                    else  # round down
                        restored_bw[m,uni_failedIPedges[b]] = max(0, floor(LP_restored_bw[uni_failedIPedges[b]]) - rand(1:rounding_range))
                        restored_bw[m,reverse_failedIPedges[b]] = restored_bw[m,uni_failedIPedges[b]]
                    end
                else
                    if rd[1] < probabilities_up[b]  # round up
                        restored_bw[m,uni_failedIPedges[b]] = min(initial_wavelength_num[uni_failedIPedges[b]], ceil(Int, LP_restored_bw[uni_failedIPedges[b]]) + rand(1:rounding_range) - 1)
                        restored_bw[m,reverse_failedIPedges[b]] = restored_bw[m,uni_failedIPedges[b]]
                    else  # round down
                        restored_bw[m,uni_failedIPedges[b]] = max(0, floor(Int, LP_restored_bw[uni_failedIPedges[b]]) - rand(1:rounding_range) + 1)
                        restored_bw[m,reverse_failedIPedges[b]] = restored_bw[m,uni_failedIPedges[b]]
                    end
                end
            end

            # filter out bad options
            if round(Int, sum(restored_bw[m,:])) <= round(Int, sum(LP_restored_bw)) && round(Int, sum(restored_bw[m,:])) > option_gap*round(Int, sum(LP_restored_bw))  # make sure the tickets are within reasonable range of total restorable capacity
                # distill randomized tickets by removing duplicates
                duplicate_sign = 0
                for r in 1:rounded_index
                    if real_restored_bw[r,:] == restored_bw[m,:]
                        duplicate_sign = 1
                        # println("duplicate_sign")
                        break
                    end
                end
                if duplicate_sign == 0  # this is an unique ticket generated from randomized rounding
                    ## if all the restorable links happens to round down from LP, then it is feasible no need to check (this is probably not a good ticket)
                    weak_ticket = 1
                    for b in 1:nrestored
                        if restored_bw[m,b] > real_restored_bw[2,b]  # second ticket is LP result
                            weak_ticket = 0
                            break
                        end
                    end
                    if ilp_check == 1 && weak_ticket == 0  # not a weak ticket and need to run feasibility check
                        Fibers = OpticalTopo["links"]
                        FibercapacityCode = OpticalTopo["capacityCode"]
                        failedIPBranchRoutingFiber = rehoused_IProutingEdge
                        failedIPbranckIndexAll = failedIPbranckindex
                        failedIPbrachIndexGroup = failedIPbrachGroup
                        rerouting_K = optical_rerouting_K
                        check_time = @timed check_sign = RestoreRWACheck(GRB_ENV, restored_bw[m,:], Fibers, FibercapacityCode, failedIPedges, failedIPBranchRoutingFiber, failedIPbranckIndexAll, failedIPbrachIndexGroup, failed_IP_initialbw, rerouting_K, verbose)
                        if verbose println("check time $(check_time[2])") end
                        if check_sign  # feasibility check passed
                            rounded_index += 1  # append this ticket to real_restored_bw as feasible ticket
                            if verbose println("Get one feasible ticket $(restored_bw[m,:])") end
                            for b in 1:nrestored
                                real_restored_bw[rounded_index,b] = restored_bw[m,b]
                            end
                            ProgressMeter.next!(progress, showvalues = [])
                        end
                    else
                        rounded_index += 1  # append this ticket to real_restored_bw as feasible ticket
                        if verbose println("Get one feasible ticket $(restored_bw[m,:])") end
                        for b in 1:nrestored
                            real_restored_bw[rounded_index,b] = restored_bw[m,b]
                        end
                        ProgressMeter.next!(progress, showvalues = [])
                    end
                end
            end

            # we get desired number of rounded results or get the theoretical max number of tickets
            if rounded_index >= min(ticket_set_size, theoretical_max_ticket)
                if verbose println("Get desired number $(rounded_index) of tickets") end
                break
            end
        end
    end

    return real_restored_bw
end

## convert failure scenarios to IP link failure information
function ReadFailureScenario(scenario, IPedges, IPcapacity)
    nedges = size(IPedges,1)
    failed_IPedge = []
    failed_IP_initialindex = []
    failed_IP_initialbw = []
    for x in 1:nedges
        if scenario[x] < 1   # link failure
            push!(failed_IPedge, IPedges[x])
            push!(failed_IP_initialindex, x)
            push!(failed_IP_initialbw, IPcapacity[x])
        end
    end

    return failed_IPedge, failed_IP_initialindex, failed_IP_initialbw
end

## routing of restored wavelength for all IP links as a database only for OPTIMAL cross-layer formulation
function AllIPWaveRerouting(OpticalTopo, IPedges, rerouting_K)
    nIPedges = size(IPedges, 1)
    rehoused_IProutingEdge = []
    rehoused_IProutingPaths = []
    IPbranckIndexAll = []
    IPbrachIndexGroup = []

    ## an optical topology without fiber cut
    num_nodes = length(OpticalTopo["nodes"])
    optical_graph = LightGraphs.SimpleDiGraph(num_nodes)
    distances = Inf*ones(num_nodes, num_nodes)

    num_edges = length(OpticalTopo["links"])
    for i in 1:num_edges
        LightGraphs.add_edge!(optical_graph, OpticalTopo["links"][i][1], OpticalTopo["links"][i][2])
        distances[OpticalTopo["links"][i][1], OpticalTopo["links"][i][2]] = OpticalTopo["fiber_length"][i]
        distances[OpticalTopo["links"][i][2], OpticalTopo["links"][i][1]] = OpticalTopo["fiber_length"][i]
    end

    # find rerouting_K paths for each IP links, pay attention: IP links are bidirectional!
    global_IPbranch = 1
    for w in 1:nIPedges
        state = LightGraphs.yen_k_shortest_paths(optical_graph, IPedges[w][1], IPedges[w][2], distances, rerouting_K)
        paths = state.paths
        # println("paths: ", paths)
        if length(paths) <= rerouting_K
            path_edges = []
            for p in 1:length(paths)
                k_path_edges = []
                for i in 1:length(paths[p])-1
                    e = findfirst(x -> x == (paths[p][i], paths[p][i+1]), OpticalTopo["links"])
                    append!(k_path_edges, e)
                end
                push!(path_edges, k_path_edges)
            end

            append!(rehoused_IProutingEdge, path_edges)
            append!(rehoused_IProutingPaths, paths)
            append!(IPbranckIndexAll, range(global_IPbranch, length=length(paths), step=1))
            push!(IPbrachIndexGroup, range(global_IPbranch, length=length(paths), step=1))
            global_IPbranch += length(paths)
        end
    end

    return rehoused_IProutingEdge, rehoused_IProutingPaths, IPbranckIndexAll, IPbrachIndexGroup
end

max_stars_repo_name: hipersys-team/arrow

max_issues_repo_name: arrow-bow/arrow

max_forks_repo_name: arrow-bow/arrow


text:
#=
testODE:
- Julia version:
- Author: ymocquar
- Date: 2019-11-13
=#
using HOODESolver
using DifferentialEquations
using LinearAlgebra
using Plots
using Random
function twoscales_solve( par_u0::PrepareU0, order, t, nb)
end
function fctMain(n_tau)
    seed=123321
    Random.seed!(seed)
    u0 = rand(BigFloat, 4)
    B = 2rand(BigFloat, 4, 4) - ones(BigFloat, 4, 4)
    println("seed = $seed B=$B")
    nbeps=48
    nbmaxtest=3
    paseps = big"10.0"^0.125
    ordmin=3
    ordmax=14
    ordpas = 1
    t_max = big"1.0"
    A = [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
    for order = ordmin:ordpas:ordmax
        y = ones(Float64, nbeps, nbmaxtest)
        x=zeros(Float64,nbeps)
        for i_eps=1:nbeps
            epsilon = big"1.0"/paseps^i_eps
            x[i_eps] = epsilon
            fct = u -> B*u
            parphi = PreparePhi(epsilon, n_tau, A, fct, B)
            println("prepareU0 eps=$epsilon n_tau=$n_tau")
            @time par_u0 = PrepareU0(parphi, order+2, u0)
            solref = getexactsol(parphi, u0, t_max)
            eps_v = convert(Float32,epsilon)
            println("epsilon = $eps_v solref=$solref")
            for indc=1:nbmaxtest
                nb = 10*10^indc
                @time parGen = PrepareTwoScalePureAB(nb, t_max, order, par_u0)
                @time sol = twoscales_pure_ab(parGen, only_end=true)
                println("solref=$solref")
                println("nb=$nb sol=$sol")
                diff=solref-sol
                println("nb=$nb dt=$(1.0/nb) normInf=$(norm(diff,Inf)) norm2=$(norm(diff))")
                y[i_eps, indc] = norm(diff,Inf)
                println("epsilon=$epsilon result=$y")
            end
        end
        labels=Array{String,2}(undef, 1, nbmaxtest)
        for i=1:nbmaxtest
            nb = 10*10^i
            labels[1,i] = " delta t,order=$(1/nb), $order "
            mindiff = minimum(y[:,i])
            bornediff = min(mindiff*1e20,1)
            for j=1:nbeps
                if y[j,i] >= bornediff
                    y[j,i] = NaN
                end
            end
        end
        gr()
        p=Plots.plot(
    x,
    y,
    xlabel="epsilon",
    xaxis=:log,
    ylabel="error",
    yaxis=:log,
    legend=:topleft,
    label=labels,
    marker=2
)
        prec_v = precision(BigFloat)
        Plots.savefig(p,"out/r_$(prec_v)_$(order)_$(n_tau)_epsilon_lgn_2.pdf")
    end
end
fctMain(32)


text:
#! format: off
using LinearElasticity
using Documenter

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

makedocs(;
    modules=[LinearElasticity],
    authors="Qi Zhang <singularitti@outlook.com>",
    repo="https://github.com/MineralsCloud/LinearElasticity.jl/blob/{commit}{path}#{line}",
    sitename="LinearElasticity.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://MineralsCloud.github.io/LinearElasticity.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/MineralsCloud/LinearElasticity.jl",
)


text:
using BoundaryValueDiffEq
using DiffEqBase, OrdinaryDiffEq, DiffEqDevTools
using Test
using NLsolve

@time @testset "Shooting Method Tests" begin
include("shooting_tests.jl")
include("orbital.jl")
end

@time @testset "Collocation Method (MIRK) Tests" begin include("mirk_convergence_tests.jl") end


text:
using GLFW
using VulkanCore
using VulkanCore.LibVulkan

@assert GLFW.VulkanSupported()

include(joinpath(@__DIR__, "vkhelper.jl"))

const WIDTH = 800
const HEIGHT = 600

## init GLFW window
GLFW.WindowHint(GLFW.CLIENT_API, GLFW.NO_API)    # not to create an OpenGL context
GLFW.WindowHint(GLFW.RESIZABLE, 0)
window = GLFW.CreateWindow(WIDTH, HEIGHT, "Vulkan")

## init Vulkan
## creating instance
appInfoRef = VkApplicationInfo(
    "Application Name: Create Instance",
    v"1.0.0",
    "No Engine Name",
    v"1.0.0",
    VK_API_VERSION_1_2,
) |> Ref

layers = String[]
extensions = GLFW.GetRequiredInstanceExtensions()
@assert check_extensions(extensions)
createInfoRef = VkInstanceCreateInfo(appInfoRef, layers, extensions) |> Ref

instanceRef = Ref(VkInstance(C_NULL))
result = GC.@preserve appInfoRef layers extensions vkCreateInstance(createInfoRef, C_NULL, instanceRef)
@assert result == VK_SUCCESS "failed to create instance!"

## main loop
while !GLFW.WindowShouldClose(window)
    GLFW.PollEvents()
end

## cleaning up
vkDestroyInstance(instanceRef[], C_NULL)
GLFW.DestroyWindow(window)

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
# NLopt constraints ------------------------------------------------------------------------

# 1.5*x1 == x2 (aka 1.5*x1 - x2 == 0)
#= function eq_constraint(x, grad)
    if length(grad) > 0
        grad .= 0.0
        grad[1] = 1.5
        grad[2] = -1.0
    end
    1.5*x[1] - x[2]
end =#

# x30*x31 ≥ 0.6 (aka 0.6 - x30*x31 ≤ 0)
function ineq_constraint(x, grad)
    if length(grad) > 0
        grad .= 0.0
        grad[30] = -x[31]
        grad[31] = -x[30]
    end
    0.6 - x[30]*x[31]
end

constrained_optimizer = SemOptimizerNLopt(;
    algorithm = :AUGLAG,
    local_algorithm = :LD_LBFGS,
    local_options = Dict(
        :ftol_rel => 1e-6
    ),
    # equality_constraints = NLoptConstraint(;f = eq_constraint, tol = 1e-14),
    inequality_constraints = NLoptConstraint(;f = ineq_constraint, tol = 0.0),
)

model_ml_constrained = Sem(
    specification = spec,
    data = dat,
    optimizer = constrained_optimizer
)

solution_constrained = sem_fit(model_ml_constrained)

# NLopt option setting ---------------------------------------------------------------------

model_ml_maxeval = Sem(
    specification = spec,
    data = dat,
    optimizer = SemOptimizerNLopt,
    options = Dict(:maxeval => 10)
)

############################################################################################
### test solution
############################################################################################

@testset "ml_solution_maxeval" begin
    solution_maxeval = sem_fit(model_ml_maxeval)
    @test solution_maxeval.optimization_result.problem.numevals == 10
    @test solution_maxeval.optimization_result.result[3] == :MAXEVAL_REACHED
end

@testset "ml_solution_constrained" begin
    solution_constrained = sem_fit(model_ml_constrained)
    @test solution_constrained.solution[31]*solution_constrained.solution[30] >= 0.6
    @test all(abs.(solution_constrained.solution) .< 10)
    @test_skip solution_constrained.optimization_result.result[3] == :FTOL_REACHED
    @test abs(solution_constrained.minimum - 21.21) < 0.01
end


text:
#
# T = Dual{Nothing,Float64,6}
#
# function partial_length(T::Type)
#     a = T(0.0)
#     return length(a.partials)
# end
#
# rand_partials(N::Int64) = Partials(Tuple(rand(N)))
#
# function rand_sv3_zero(T::Type)
#     (T == Float64) && (return zeros(SVector{3,Float64}))
#     N = partial_length(T)
#     v1 = T(0.0, rand_partials(N))
#     v2 = T(0.0, rand_partials(N))
#     v3 = T(0.0, rand_partials(N))
#     return SVector{3,T}(v1, v2, v3)
# end
#
# function rand_sv3(T::Type)
#     (T == Float64) && (return randn(SVector{3,Float64}))
#     N = partial_length(T)
#     v1 = T(randn(), rand_partials(N))
#     v2 = T(randn(), rand_partials(N))
#     v3 = T(randn(), rand_partials(N))
#     return SVector{3,T}(v1, v2, v3)
# end
#
# is_nan_in_partials(b::Float64) = false
# function is_nan_in_partials(b::Dual{Tag_Type,Float64,N}) where {Tag_Type,N}
#     return any(isnan.(b.partials))
# end
# function is_nan_in_partials(b::SVector{3,T}) where {T}
#     return any(is_nan_in_partials.(b))
# end
#
# @testset "div_by_zero" begin
#     for d_type = [Float64, T]
#         sv3_zero = rand_sv3_zero(d_type)
#         safe_norm_zero = safe_norm(sv3_zero)
#         @test !is_nan_in_partials(safe_norm_zero)
#         @test !is_nan_in_partials(safe_inv_norm(sv3_zero))
#         @test sv3_zero == safe_normalize(sv3_zero)
#         @test !is_nan_in_partials(safe_normalize(sv3_zero))
#         @test !is_nan_in_partials(safe_inv_norm²(sv3_zero))
#         for k = 1:128
#             rv = rand_sv3(d_type)
#             rv_norm = norm(rv)
#             inv_rv_norm = 1 ./ rv_norm
#             rv_normalize = rv .* inv_rv_norm
#             @test norm²(rv) ≈ rv_norm^2
#             @test unsafe_norm(rv) ≈ rv_norm
#             @test safe_norm(rv) ≈ rv_norm
#             @test unsafe_inv_norm(rv) ≈ inv_rv_norm
#             @test safe_inv_norm(rv) ≈ inv_rv_norm
#             @test safe_normalize(rv) ≈ rv_normalize
#             @test unsafe_normalize(rv) ≈ rv_normalize
#             @test safe_inv_norm²(rv) ≈ inv_rv_norm^2
#         end
#         @test zero(T) == safe_scalar_divide(zero(T), rand(T) + 0.5)
#         @test zero(T) == safe_scalar_divide(zero(T), zero(T))
#         @test_throws ErrorException safe_scalar_divide(one(T), zero(T))
#     end
# end
#
# @testset "div by zero type stability test" begin
#     for d_type = [Float64, T]
#         rv = rand_sv3(d_type)
#         rv_0 = rand_sv3_zero(d_type)
#         for fun_name = [norm², unsafe_norm, safe_norm, unsafe_inv_norm, safe_inv_norm, safe_normalize,
#             unsafe_normalize, safe_inv_norm²]
#
#             ans_1 = fun_name(rv)
#             ans_2 = fun_name(rv_0)
#             is_type_stable = typeof(ans_1) == typeof(ans_2)
#             if !is_type_stable
#                 println("for DataType $d_type function $fun_name is not type stable")
#             end
#             @test is_type_stable
#         end
#     end
# end


text:
include("LQSolver.jl")
include("subsolver.jl")
include("executors.jl")
include("penaltyterm.jl")


text:
module MOF

import ..FileFormats
import OrderedCollections
import JSON
import VecMathOptInterface

const MOI = VecMathOptInterface

const SCHEMA_PATH = joinpath(@__DIR__, "mof.1.0.schema.json")
const VERSION = v"1.0"
const SUPPORTED_VERSIONS = (v"1.0", v"0.6", v"0.5", v"0.4")

const OrderedObject = OrderedCollections.OrderedDict{String,Any}
const UnorderedObject = Dict{String,Any}
const Object = Union{OrderedObject,UnorderedObject}

function _parse_mof_version(version)
    return VersionNumber(version["major"], version["minor"])
end

struct Nonlinear <: MOI.AbstractScalarFunction
    expr::Expr
end
Base.copy(nonlinear::Nonlinear) = Nonlinear(copy(nonlinear.expr))
MOI.Utilities.canonicalize!(nonlinear::Nonlinear) = nonlinear

MOI.Utilities.@model(
    InnerModel,
    (MOI.ZeroOne, MOI.Integer),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
        MOI.Semicontinuous,
        MOI.Semiinteger,
    ),
    (
        MOI.Reals,
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.ExponentialCone,
        MOI.DualExponentialCone,
        MOI.NormOneCone,
        MOI.NormInfinityCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.RelativeEntropyCone,
        MOI.RootDetConeTriangle,
        MOI.RootDetConeSquare,
        MOI.LogDetConeTriangle,
        MOI.LogDetConeSquare,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.PositiveSemidefiniteConeSquare,
    ),
    (MOI.PowerCone, MOI.DualPowerCone, MOI.SOS1, MOI.SOS2),
    (Nonlinear,),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

# Indicator is handled by UniversalFallback.

const Model = MOI.Utilities.UniversalFallback{InnerModel{Float64}}

struct Options
    print_compact::Bool
    warn::Bool
end

function get_options(m::Model)
    return get(m.model.ext, :MOF_OPTIONS, Options(false, false))
end

"""
    Model(; kwargs...)

Create an empty instance of FileFormats.MOF.Model.

Keyword arguments are:

 - `print_compact::Bool=false`: print the JSON file in a compact format without
   spaces or newlines.

 - `warn::Bool=false`: print a warning when variables or constraints are renamed
"""
function Model(; print_compact::Bool = false, warn::Bool = false)
    model = MOI.Utilities.UniversalFallback(InnerModel{Float64}())
    model.model.ext[:MOF_OPTIONS] = Options(print_compact, warn)
    return model
end

function Base.show(io::IO, ::Model)
    print(io, "A MathOptFormat Model")
    return
end

include("nonlinear.jl")

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

end


text:
@testset "basis.jl" begin
    @testset "ZBasis" begin
        @testset "FullZBasis" begin
            @test_throws ArgumentError zbasis(-1)
            @test zbasis(5) isa FullZBasis
            @test SpinSymmetry._indices(zbasis(5)) == 1:2^5

            @test zbasis(5) == FullZBasis(5)
            @test zbasis(6) != FullZBasis(5)
            @test hash(zbasis(5)) == hash(FullZBasis(5))
            @test hash(zbasis(6)) != hash(FullZBasis(5))
        end

        @testset "ZBlockBasis" begin
            @test_throws ArgumentError zbasis(-1, 0)
            @test_throws ArgumentError zbasis(2,-1)
            @test_throws ArgumentError zbasis(2,3)

            @test zbasis(2,1) isa ZBlockBasis

            @test zbasis(5,2) == ZBlockBasis(5,2)
            @test zbasis(6,2) != ZBlockBasis(5,2)
            @test zbasis(6,3) != ZBlockBasis(6,2)
            @test hash(zbasis(5,2)) == hash(ZBlockBasis(5,2))
            @test hash(zbasis(6,2)) != hash(ZBlockBasis(5,2))
            @test hash(zbasis(6,3)) != hash(ZBlockBasis(6,2))

            # small correctness check
            @test SpinSymmetry._indices(zbasis(2,0)) == [1]
            @test SpinSymmetry._indices(zbasis(2,1)) == [2,3]
            @test SpinSymmetry._indices(zbasis(2,2)) == [4]

            digitsum(k) = sum(parse.(Int, [string(k-1; base=2)...]))
            for k in 0:10
                zblock = zbasis(10, k)
                @test all(digitsum.(SpinSymmetry._indices(zblock)) .== k)
            end
        end
    end

    @testset "symmetrized_basis" begin
        # different argument forms
        @test symmetrized_basis(6) == symmetrized_basis(zbasis(6))
        @test symmetrized_basis(6, 2) == symmetrized_basis(zbasis(6, 2))
        @test symmetrized_basis(6) isa SymmetrizedBasis
        @test symmetrized_basis(zbasis(6), Flip(6), 0) == symmetrized_basis(6, Flip(6), 0)
        @test symmetrized_basis(zbasis(6, 2), Flip(6), 0) == symmetrized_basis(6, 2, Flip(6), 0)

        # equality
        @test symmetrized_basis(6, Flip(6), 0, Shift(6), 1) == symmetrized_basis(6, Shift(6), 1, Flip(6), 0)
        @test symmetrized_basis(6, Flip(6), 0, Shift(6, 5), 1) == symmetrized_basis(6, Shift(6, -1), 1, Flip(6), 0)
        @test symmetrized_basis(6, Flip(6), 0, Shift(6), 0) != symmetrized_basis(6, Shift(6), 1, Flip(6), 0)
        @test symmetrized_basis(6, Flip(6), 1, Shift(6), 0) != symmetrized_basis(6, Shift(6), 1, Flip(6), 0)

        let b1 = symmetrized_basis(6, Flip(6), 0, Shift(6), 1),
            b2 = symmetrized_basis(6, Shift(6), 1, Flip(6), 0)
            @test SpinSymmetry._phase_factors(1:2^6, b1.symmetries, b1.sectors)  ==
                SpinSymmetry._phase_factors(1:2^6, b2.symmetries, b2.sectors)
        end
    end

    @testset "basissize" begin
        let basis = zbasis(10)
            @test basissize(basis) == length(SpinSymmetry._indices(basis))
        end

        for k in 0:10
            let basis = zbasis(10, k)
                @test basissize(basis) == length(SpinSymmetry._indices(basis))
            end
        end

        @test basissize(symmetrized_basis(10, Flip(10), 0)) == 2^9
        @test basissize(symmetrized_basis(10, 5,  Flip(10), 0)) == binomial(10, 5)/2
        @test basissize(symmetrized_basis(10, 4,  Flip(10), 0)) == binomial(10, 4)
    end

    @testset "symmetrize stuff" begin
        using LinearAlgebra, Random, SparseArrays
        Random.seed!(5)

        @testset "return types" begin
            N = 11

            real_basis1 = zbasis(N, div(N-1, 2))
            real_basis2 = symmetrized_basis(N, Shift(N), 0)
            complex_basis = symmetrized_basis(N, Shift(N), 1)
            σz = spdiagm([1,-1])
            make_op(op, N) = sum(kron(kron(I(2^(N-i)), op), I(2^(i-1))) for i in 1:N)
            op = make_op(σz, N)

            for (op, should_be_sparse) in ((op, true), (Matrix(op), false))
                let symm = symmetrize_operator(op, real_basis1)
                    @test issparse(symm) == should_be_sparse
                    @test eltype(symm) <: Real
                end
                let symm = symmetrize_operator(op, real_basis2)
                    @test issparse(symm) == should_be_sparse
                    @test eltype(symm) <: Real
                end
                let symm = symmetrize_operator(op, complex_basis)
                    @test issparse(symm) == should_be_sparse
                    @test eltype(symm) <: Complex
                end
            end
        end

        N = 5

        state = normalize!(rand(2^N))
        operator = rand(2^N, 2^N)
        base = symmetrized_basis(N, Flip(N), 0)
        mid = 2^(N-1)
        front = 1:mid
        back = 2^N:-1:mid+1

        @test symmetrize_state(state, symmetrized_basis(N, 0)) ≈ [state[1]]
        @test symmetrize_state(state, zbasis(N, N)) ≈ [state[end]]
        @test symmetrize_operator(operator, symmetrized_basis(N, 0)) ≈ [operator[1,1]]
        @test symmetrize_operator(operator, zbasis(N, N)) ≈ [operator[end, end]]

        symm_state = symmetrize_state(state, base)
        @test symm_state ≈ symmetrize_state(state, N, Flip(N), 0)
        @test symm_state ≈ 1/√2 .* (state[front] .+ state[back])

        asymm_state = symmetrize_state(state, N, Flip(N), 1)
        @test asymm_state ≈ 1/√2 .* (state[front] .- state[back])

        @test_throws ArgumentError symmetrize_state(zeros(2^N+1), base)

        symm_operator = symmetrize_operator(operator, base)
        @test symm_operator ≈ symmetrize_operator(operator, N, Flip(N), 0)
        @test symm_operator ≈ 1/2 .* (operator[front, front] .+ operator[back, back] .+ operator[front, back] .+ operator[back, front])

        asymm_operator = symmetrize_operator(operator, N, Flip(N), 1)
        @test asymm_operator ≈ 1/2 .* (operator[front, front] .+ operator[back, back] .- operator[front, back] .- operator[back, front])

        @test_throws ArgumentError symmetrize_operator(zeros(2^N+1, 2^N+1), base)
        @test_throws ArgumentError symmetrize_operator(zeros(2^(N-1), 2^(N+1)), base)
    end
end


text:
# ---
# title: 221. Maximal Square
# id: problem221
# author: Indigo
# date: 2021-05-29
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/maximal-square/description/>
# hidden: true
# ---
#
# Given an `m x n` binary `matrix` filled with `0`'s and `1`'s, _find the
# largest square containing only_ `1`'s _and return its area_.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/26/max1grid.jpg)
#
#
#
#     Input: matrix = [["1","0","1","0","0"],["1","0","1","1","1"],["1","1","1","1","1"],["1","0","0","1","0"]]
#     Output: 4
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/26/max2grid.jpg)
#
#
#
#     Input: matrix = [["0","1"],["1","0"]]
#     Output: 1
#
#
# **Example 3:**
#
#
#
#     Input: matrix = [["0"]]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `m == matrix.length`
#   * `n == matrix[i].length`
#   * `1 <= m, n <= 300`
#   * `matrix[i][j]` is `'0'` or `'1'`.
#
#
## @lc code=start
using LeetCode

function maximal_square(matrix::Matrix{Char})
    length(matrix) == 0 && return 0
    maxs = any(==('1'), matrix)
    maxs == 0 && return 0
    m, n = size(matrix)
    dp = fill(0, size(matrix))
    for i in 1:m
        dp[i, 1] = matrix[i, 1] == '1'
    end
    for j in 1:n
        dp[1, j] = matrix[1, j] == '1'
    end
    for i in 2:m, j in 2:n
        if matrix[i, j] == '1'
            dp[i, j] = min(dp[i, j - 1], dp[i - 1, j], dp[i - 1, j - 1]) + 1
            maxs = max(maxs, dp[i, j])
        end
    end
    return maxs * maxs
end
## @lc code=end


text:
# Read in file of helper functions
# This file incluse: cleantext, countwords
include("../word_count_helpers.jl")

# Load the file names
dataLoc = "../../../data/word_count/";
fnames = dataLoc.*readdir(dataLoc)

# Iterate through file names
for fname in fnames

    # Read in file and clean the text
    f = open(fname)
    text = readlines(f)
    text = cleantext(text)

    # Count number of times each word appears
    wordCounts = countwords(text)

    # Sort and print the top 5 words with their counts
    rankedwords = sort(collect(wordCounts), by=x->x[2], rev=true)
    println(rankedwords[1:5,:])

end

max_stars_repo_name: llsc-supercloud/teaching-examples

max_issues_repo_name: llsc-supercloud/teaching-examples

max_forks_repo_name: llsc-supercloud/teaching-examples


text:
import DPFEHM
import GaussianRandomFields
import Optim
import PyPlot
import Random
import Zygote

Random.seed!(0)

#set up the grid
mins = [0, 0, 0]; maxs = [50, 50, 1]#size of the domain, in meters
ns = [30, 30, 3]#number of nodes on the grid
coords, neighbors, areasoverlengths, _ = DPFEHM.regulargrid3d(mins, maxs, ns)#build the grid

#set up the boundary conditions
Qs = zeros(size(coords, 2))
injectionnode = 1#inject in the lower left corner
Qs[injectionnode] = 1e-6#m^3/s
dirichletnodes = Int[size(coords, 2)]#fix the pressure in the upper right corner
dirichleths = zeros(size(coords, 2))
dirichleths[size(coords, 2)] = 0.0

#set up the conductivity field
lambda = 50.0#meters -- correlation length of log-conductivity
sigma = 1.0#standard deviation of log-conductivity
mu = -9.0#mean of log conductivity -- ~1e-4 m/s, like clean sand here https://en.wikipedia.org/wiki/Hydraulic_conductivity#/media/File:Groundwater_Freeze_and_Cherry_1979_Table_2-2.png
cov = GaussianRandomFields.CovarianceFunction(2, GaussianRandomFields.Matern(lambda, 1; σ=sigma))

alphas = fill(0.375, length(neighbors))
Ns = fill(1.25, length(neighbors))

x_pts = range(mins[1], maxs[1]; length=ns[1])
y_pts = range(mins[2], maxs[2]; length=ns[2])
num_eigenvectors = 200
grf = GaussianRandomFields.GaussianRandomField(cov, GaussianRandomFields.KarhunenLoeve(num_eigenvectors), x_pts, y_pts)
logKs = zeros(reverse(ns)...)
logKs2d = mu .+ GaussianRandomFields.sample(grf)'#generate a random realization of the log-conductivity field
for i = 1:ns[3]#copy the 2d field to each of the 3d layers
    v = view(logKs, i, :, :)
    v .= logKs2d
end

#plot the log-conductivity
fig, ax = PyPlot.subplots()
img = ax.imshow(logKs[1, :, :], origin="lower")
ax.title.set_text("Conductivity Field")
fig.colorbar(img)
display(fig)
println()
PyPlot.close(fig)

logKs2Ks_neighbors(Ks) = exp.(0.5 * (Ks[map(p->p[1], neighbors)] .+ Ks[map(p->p[2], neighbors)]))
function solveforh(logKs, dirichleths)
    @assert length(logKs) == length(Qs)
    if maximum(logKs) - minimum(logKs) > 25
        return fill(NaN, length(Qs))#this is needed to prevent the solver from blowing up if the line search takes us somewhere crazy
    else
        Ks_neighbors = logKs2Ks_neighbors(logKs)
        psi0 = zeros(size(coords, 2))
        return reshape(DPFEHM.richards_steadystate(psi0, Ks_neighbors, neighbors, areasoverlengths, dirichletnodes, dirichleths, coords, alphas, Ns, Qs), reverse(ns)...)
    end
end

#now do a forward solve for the head
print("forward solve time")
@time h = solveforh(logKs, dirichleths)#solve for the head
#plot the head at the bottom of the domain
fig, ax = PyPlot.subplots()
img = ax.imshow(h[1, :, :], origin="lower")
ax.title.set_text("Head")
fig.colorbar(img)
display(fig)
println()
PyPlot.close(fig)

#now compute the gradient of a function involving solveforh
gradient_node = div(size(coords, 2), 2) + 500
gradient_node_x = coords[1, gradient_node]
gradient_node_y = coords[2, gradient_node]
print("forward and gradient time")
@time grad = Zygote.gradient((x, y)->solveforh(x, y)[gradient_node], logKs, dirichleths)#calculate the gradient (which involves a redundant calculation of the forward pass)
function_evaluation, back = Zygote.pullback((x, y)->solveforh(x, y)[gradient_node], logKs, dirichleths)#this pullback thing lets us not redo the forward pass
print("gradient time")
@time grad = back(1.0)#compute the gradient of a function involving solveforh
#plot the gradient of the function w.r.t. the logK at the bottom of the domain
fig, ax = PyPlot.subplots()
img = ax.imshow(grad[1][1, :, :], origin="lower", extent=[mins[1], maxs[1], mins[2], maxs[2]])
ax.plot([gradient_node_x], [gradient_node_y], "r.", ms=10, alpha=0.5)
ax.title.set_text("Gradient of head at dot w.r.t. logK at bottom of domain")
fig.colorbar(img)
display(fig)
println()
PyPlot.close(fig)


text:
using Serialization
using LightGraphs
using SimpleWeightedGraphs
using MetaGraphs
using Bijections
using Statistics
import LightGraphs.Parallel
import Base.Threads.@threads

mutable struct GGraphData
    genotype_vertex_mapping::Bijection
    genotype_graph::SimpleWeightedDiGraph

    phenotype_graph::SimpleWeightedDiGraph

    neutral_networks::Vector
    neutral_network_graph::SimpleWeightedDiGraph

    fitness_values::Dict

    GGraphData() = new(Bijection(), SimpleWeightedDiGraph(0), SimpleWeightedDiGraph(0), [], SimpleWeightedDiGraph(0), Dict())
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
# This code demonstrates the use of SolveDSGE for solving rational expectations models
# to second-order accuracy.  The demonstration models are the first and second examples
# from Schmitt-Grohe and Uribe (2004)

#-------------------------------------------------------------------------------------------#

# Bring the SolveDSGE module into global scope

using SolveDSGE

# Bring the NLsolve module into global scope

# Example one from Schmitt-Grohe and Uribe (2004)

# Specify the model in terms of a function whose fix-point in the model's nonstochastic steady state

function model_eqns(x)

  f = zeros(length(x))

  f[1] = x[4] - rho*x[1]
  f[2] = (1.0-delta)*exp(x[2]) + exp(x[1])*exp(x[2])^alpha - exp(x[3]) - exp(x[5])
  f[3] = exp(x[3])^(-gama) - betta*exp(x[6])^(-gama)*(1.0-delta+alpha*exp(x[4])*exp(x[5])^(alpha-1.0))
  f[4] = x[4] - x[1]
  f[5] = x[5] - x[2]
  f[6] = x[6] - x[3]

#  println(x)

  return f

end

# Specify the model in terms of a function to compute its first derivatives
# The function above (model_eqns) can also be used for this purpose

function model_eqns_2(x)

  f = zeros(3)

  f[1] = x[4] - rho*x[1]
  f[2] = (1.0-delta)*exp(x[2]) + exp(x[1])*exp(x[2])^alpha - exp(x[3]) - exp(x[5])
  f[3] = exp(x[3])^(-gama) - betta*exp(x[6])^(-gama)*(1.0-delta+alpha*exp(x[4])*exp(x[5])^(alpha-1.0))

  return f

end

# Specify the three model equations separately in order to compute the second derivatives

# Technology shock equation

function model_eqn_1(x)

  f = [0.0]

  f[1] = x[4] - rho*x[1]

  return f

end

# Capital accumulation equation

function model_eqn_2(x)

  f = [0.0]

  f[1] = (1.0-delta)*exp(x[2]) + exp(x[1])*exp(x[2])^alpha - exp(x[3]) - exp(x[5])

  return f

end

# Consumption Euler equation

function model_eqn_3(x)

  f = [0.0]

  f[1] = exp(x[3])^(-gama) - betta*exp(x[6])^(-gama)*(1.0-delta+alpha*exp(x[4])*exp(x[5])^(alpha-1.0))

  return f

end

# Assign values to model parameters

betta  = 0.95
alpha = 0.3
gama  = 2.0
delta = 1.0
rho   = 0.0
etta   = [1.0, 0.0]
sigma = [1.0]

# Initial guess for computing steady state

x = [0.01, -1.8, -0.9]
x = [x; x]

# Solve for the nonstochastic steady state

(z, f_zero, iters) = newton(model_eqns, x, 1e-14, 1000)

ss = copy(z)

# Compute the model's first derivatives at the nonstochastic steady state

deriv1 = derivative(model_eqns_2,ss)

# Compute the model's second derivatives at the nonstochastic steady state

deriv21 = hessian(model_eqn_1,ss)
deriv22 = hessian(model_eqn_2,ss)
deriv23 = hessian(model_eqn_3,ss)

deriv2 = [deriv21; deriv22; deriv23]

nx = 2
ny = 1
cutoff = 1.0

# Put the model into Gomme_Klein_Form type

m_gk = Gomme_Klein_Form(nx,ny,deriv1,deriv2,etta,sigma)

# Do some basic checking to determine whether the model matrices are conformable

assessment = check_model_form(m_gk)
println(assessment)

# Solve the model using Gomme and Klein (2011)

soln_gk = solve_re(m_gk,cutoff)
responses_gk = impulses(soln_gk,5,1)

# Put the model into Lombardo_Sutherland_Form type

m_ls = Lombardo_Sutherland_Form(nx,ny,deriv1,deriv2,etta,sigma)

# Solve the model using Lombardo and Sutherland (2007)

soln_ls = solve_re(m_ls,cutoff)
responses_ls_1 = impulses(soln_ls,5,1)

# Convert Lombardo_Sutherland_Soln to Gomme_Klein_Soln

new_soln_ls = convert_second_order(soln_ls)
responses_ls_2 = impulses(new_soln_ls,5,1)

# Check to see that the two solutions are the same

println(maximum(abs,new_soln_ls.hxx-soln_gk.hxx))
println(maximum(abs,new_soln_ls.gxx-soln_gk.gxx))

println(maximum(abs,responses_gk-responses_ls_1))

# Example two from Schmitt-Grohe and Uribe (2004)

# Specify the model in terms of a function whose fix-point in the model's nonstochastic steady state

function model2_eqns(x)

  f = zeros(length(x))

  f[1]  = x[7] - rho*x[1]
  f[2]  = x[8] - rho*x[2]
  f[3]  = exp(x[5]) + exp(x[6]) + exp(x[9]) - (1.0-delta)*exp(x[3]) + exp(x[10]) - (1.0-delta)*exp(x[4]) - exp(x[1])*exp(x[3])^alpha - exp(x[2])*exp(x[4])^alpha
  f[4]  = exp(x[5]) - exp(x[6])
  f[5]  = exp(x[5])^(-gama) - betta*exp(x[11])^(-gama)*(1.0-delta+alpha*exp(x[7])*exp(x[9])^(alpha-1.0))
  f[6]  = exp(x[6])^(-gama) - betta*exp(x[12])^(-gama)*(1.0-delta+alpha*exp(x[8])*exp(x[10])^(alpha-1.0))
  f[7]  = x[7] - x[1]
  f[8]  = x[8] - x[2]
  f[9]  = x[9] - x[3]
  f[10] = x[10] - x[4]
  f[11] = x[11] - x[5]
  f[12] = x[12] - x[6]

  return f

end

# Specify the model in terms of a function to compute its first derivatives
# The function above (model2_eqns) can also be used for this purpose

function model2_eqns_2(x)

  f = zeros(length(x))

  f[1]  = x[7] - rho*x[1]
  f[2]  = x[8] - rho*x[2]
  f[3]  = exp(x[5]) + exp(x[6]) + exp(x[9]) - (1.0-delta)*exp(x[3]) + exp(x[10]) - (1.0-delta)*exp(x[4]) - exp(x[1])*exp(x[3])^alpha - exp(x[2])*exp(x[4])^alpha
  f[4]  = exp(x[5]) - exp(x[6])
  f[5]  = exp(x[5])^(-gama) - betta*exp(x[11])^(-gama)*(1.0-delta+alpha*exp(x[7])*exp(x[9])^(alpha-1.0))
  f[6]  = exp(x[6])^(-gama) - betta*exp(x[12])^(-gama)*(1.0-delta+alpha*exp(x[8])*exp(x[10])^(alpha-1.0))

  return f

end

# Specify the three model equations separately in order to compute the second derivatives

# First technology shock equation

function model2_eqn_1(x)

  f = [0.0]

  f[1]  = x[7] - rho*x[1]

  return f

end

# Second technology shock equation

function model2_eqn_2(x)

  f = [0.0]

  f[1]  = x[8] - rho*x[2]

  return f

end

# Capital accumulation equation

function model2_eqn_3(x)

  f = [0.0]

  f[1]  = exp(x[5]) + exp(x[6]) + exp(x[9]) - (1.0-delta)*exp(x[3]) + exp(x[10]) - (1.0-delta)*exp(x[4]) - exp(x[1])*exp(x[3])^alpha - exp(x[2])*exp(x[4])^alpha

  return f

end

# Simple way of specifying second capital accumulation equation

function model2_eqn_4(x)

  f = [0.0]

  f[1]  = exp(x[5]) - exp(x[6])

  return f

end

# Euler equation for first consumption good

function model2_eqn_5(x)

  f = [0.0]

  f[1]  = exp(x[5])^(-gama) - betta*exp(x[11])^(-gama)*(1.0-delta+alpha*exp(x[7])*exp(x[9])^(alpha-1.0))

  return f

end

# Euler equation for second consumption good

function model2_eqn_6(x)

  f = [0.0]

  f[1]  = exp(x[6])^(-gama) - betta*exp(x[12])^(-gama)*(1.0-delta+alpha*exp(x[8])*exp(x[10])^(alpha-1.0))

  return f

end

# Assign values to model parameters

betta = 0.95
alpha = 0.3
delta = 0.1
gama  = 2.0
rho   = 0.0
etta   = [1.0 0.0; 0.0 1.0; 0.0 0.0; 0.0 0.0]
sigma = [1.0 0.0; 0.0 1.0]

# Initial guess for computing steady state

x = [0.1, 0.1, 1.0, 1.0, 0.1, 0.1]
x = [x; x]

# Solve for the nonstochastic steady state
(z, f_zero, iters) = newton(model2_eqns, x, 1e-14, 1000)

ss = copy(z)

# Compute model's first derivatives at steady state

deriv1 = derivative(model2_eqns_2,ss)

# Compute model's second derivatives at steady state

deriv21 = hessian(model2_eqn_1,ss)
deriv22 = hessian(model2_eqn_2,ss)
deriv23 = hessian(model2_eqn_3,ss)
deriv24 = hessian(model2_eqn_4,ss)
deriv25 = hessian(model2_eqn_5,ss)
deriv26 = hessian(model2_eqn_6,ss)

deriv2 = [deriv21; deriv22; deriv23; deriv24; deriv25; deriv26]

nx = 4
ny = 2
cutoff = 1.0

# Put the model into Gomme_Klein_Form type

m_gk = Gomme_Klein_Form(nx,ny,deriv1,deriv2,etta,sigma)

# Do some basic checking to determine whether the model matrices are conformable

assessment = check_model_form(m_gk)
println(assessment)

# Solve the model using Gomme and Klein (2011)

soln_gk = solve_re(m_gk,cutoff)
responses_gk = impulses(soln_gk,5,1)

# Put the model into Lombard_Sutherland_Form type

m_ls = Lombardo_Sutherland_Form(nx,ny,deriv1,deriv2,etta,sigma)

# Solve the model using Lombardo and Sutherland (2007)

soln_ls = solve_re(m_ls,cutoff)
responses_ls_1 = impulses(soln_ls,5,1)

# Convert Lombardo_Sutherland_Soln to Gomme_Klein_Soln

new_soln_ls = convert_second_order(soln_ls)
responses_ls_2 = impulses(new_soln_ls,5,1)

# Check to see that the two solutions are the same

println(maximum(abs,new_soln_ls.hxx-soln_gk.hxx))
println(maximum(abs,new_soln_ls.gxx-soln_gk.gxx))

println(maximum(abs,responses_gk-responses_ls_1))


text:
let
    n = 50
    p = 2
    t = (1:n)/n
    Ut,Vt = egrss.generators(t,p);
    K = egrss.full(Ut,Vt)

    Wt,c = egrss.ldl(Ut,Vt)
    Lref = cholesky(K).L
    L = egrss.full_tril(Ut,Wt,1.0)
    @test( isapprox(L,Lref./diag(Lref)') )
    @test( isapprox(map(sqrt,c),diag(Lref)) )

    alpha = 1e-6
    Wt,c = egrss.ldl(Ut,Vt,alpha)
    Lref = cholesky(egrss.full(Ut,Vt,alpha)).L
    L = egrss.full_tril(Ut,Wt,1.0)
    @test( isapprox(L,Lref./diag(Lref)') )
    @test( isapprox(map(sqrt,c),diag(Lref)) )

    d = collect((1:n)/n) .+ 1e-3
    Wt,c = egrss.ldl(Ut,Vt,d)
    Lref = cholesky(K + Diagonal(d)).L
    L = egrss.full_tril(Ut,Wt,1.0)
    @test( isapprox(L,Lref./diag(Lref)') )
    @test( isapprox(map(sqrt,c),diag(Lref)) )

end

max_stars_repo_name: martinandersen/EGRSS

max_issues_repo_name: martinandersen/EGRSS

max_forks_repo_name: martinandersen/EGRSS


text:
mutable struct SingleLayer{F, FP, A}
    σ::F
    σ′::FP
    W::A
end

SingleLayer(in::Integer, out::Integer, σ, σ′; init=(dims...)->zeros(Float32, dims...)) =
    SingleLayer(σ, σ′, init(out, in))

(layer::SingleLayer)(x) = layer.σ.(layer.W*x)
deriv(layer::SingleLayer, x) = layer.σ′.(layer.W*x)

# Sparse Updates

Linear(in::Integer, out::Integer) =
    SingleLayer(in, out, identity, (x)->1.0)

# SparseLayer(num_weights, args...; kwargs...) =
#     SingleLayer(num_weights, out, args...; kwargs...)

# sigmoid = Flux.sigmoid
sigmoid′(x) = sigmoid(x)*(1.0-sigmoid(x))
identity′(x) = 1.0
tanh′(x) = 1 - (tanh(x))^2

mutable struct TabularLayer{A}
    W::A
end
Tabular(dims::Integer...) = TabularLayer(zeros(dims...))

@forward TabularLayer.W Base.getindex, Base.setindex!

(layer::TabularLayer)(x) = layer[x]
(layer::TabularLayer)(x::Array{T, 1}) where {T<:Integer} = layer[x...]

mutable struct SparseLayer{F, FP, A, P}
    σ::F
    σ′::FP
    W::A
    ϕ::P
end

SparseLayer(num_weights::Integer, out::Integer, σ=identity, σ′=identity′; init=(dims...)->zeros(Float32, dims...)) =
    SparseLayer(σ, σ′, init(out, num_weights), zeros(Float32, out, num_weights))

(layer::SparseLayer)(x::Array{Int64, 1}) = layer.σ.(sum(layer.W[x]))
deriv(layer::SparseLayer, x::Array{Int64, 1}) = layer.σ′.(sum(layer.W[x]))

(layer::SparseLayer)(x::Array{CartesianIndex{1}, 1}) = layer.σ.(sum(layer.W[x]))
deriv(layer::SparseLayer, x::Array{CartesianIndex{1}, 1}) = layer.σ′.(sum(layer.W[x]))


text:
#Distribute a optigraph among workers.  Each worker should have the same first stage model.  Each worker will be allocated some of the nodes in the original modelgraph
function distribute_optigraph(graph::OptiGraph,to_workers::Vector{Int64};remote_name = :graph)
    #NOTE: Does not yet support subgraphs.  Need to aggregate first.
    to_workers = sort(to_workers)
    n_nodes = num_nodes(graph)
    n_workers = length(to_workers)
    nodes_per_worker = Int64(floor(n_nodes/n_workers))
    nodes = all_nodes(graph)
    node_indices = [getindex(graph,node) for node in nodes]
    user_pips_data = PipsNLP._setup_pips_nlp_data!(graph)

    #Allocate optinodes onto workers.  Split it up evenly
    allocations = []
    node_indices = []
    j = 1
    while  j <= n_nodes
        if j + nodes_per_worker > n_nodes
            push!(allocations,nodes[j:end])
            push!(node_indices,[getindex(graph,node) for node in nodes[j:end]])
        else
            push!(allocations,nodes[j:j+nodes_per_worker - 1])
            push!(node_indices,[getindex(graph,node) for node in nodes[j:j+nodes_per_worker - 1]])
        end
        j += nodes_per_worker
    end

    println("Distributing optigraph among workers: $to_workers")
    channel = RemoteChannel(1)
    remote_references = []
    #Fill channel with sets of nodes to send
    #TODO: Make this runs in parallel with an @async?
    @sync begin
        for (i,worker) in enumerate(to_workers)
            @spawnat(1, put!(channel, allocations[i]))
            ref1 = @spawnat worker begin
                Core.eval(Main, Expr(:(=), :nodes, take!(channel)))
            end
            wait(ref1)
            ref2 = @spawnat worker Core.eval(Main, Expr(:(=), remote_name, PipsNLP._create_worker_optigraph(getfield(Main,:nodes),
            node_indices[i],
            user_pips_data,
            n_nodes)))
            push!(remote_references,ref2)
        end
    end
    #graph.ext[:allocations] = allocations
    return remote_references
end

function fill_solution!(graph::OptiGraph,remote_name::Symbol,worker::Int64)
    for i = 1:num_all_nodes(graph)
        x_sol = @fetchfrom worker getnode(getfield(Main,remote_name),i).ext[:pips_data].x_sol
        status = @fetchfrom worker termination_status(getnode(getfield(Main,remote_name),i))

        node = getnode(graph,i)
        Plasmo.set_node_primals(node,all_variables(node),x_sol)
        Plasmo.set_node_status(node,status)
    end
    return nothing
end

function _create_worker_optigraph(optinodes::Vector{OptiNode},
    node_indices::Vector{Int64},
    user_pips_data,
    n_nodes)

    graph = OptiGraph()
    graph.node_idx_map = Dict{OptiNode,Int64}()
    graph.ext[:user_pips_data] = user_pips_data

    #Add nodes to worker's graph.  Each worker should have the same number of nodes, but some will be empty.
    for i = 1:n_nodes
        node = add_node!(graph)
        node.ext[:pips_data] = PIPSNLPData()  #all nodes need this
    end
    for (i,index) in enumerate(node_indices)
        graph.optinodes[index] = optinodes[i]
    end

    return graph
end


text:
using Clang.Generators
using Clang.LibClang.Clang_jll

cd(@__DIR__)

# Location of the extracted wasmer.tar.gz
wasmer_location = joinpath(get(ENV, "WASMER_LOCATION", "../wasmer"), "lib/c-api/")

options = load_options(joinpath(@__DIR__, "generator.toml"))

# add compiler flags
args = get_default_args()
push!(args, "-I" * wasmer_location)

headers = [
    joinpath(wasmer_location, "../../target/release/build/wasmer-c-api-56a33427fae3bc8b/out/wasmer.h"),
]

# create context
ctx = create_context(headers, args, options)

# build without printing so we can do custom rewriting
build!(ctx, BUILDSTAGE_NO_PRINTING)

# custom rewriter
function rewrite!(e::Expr)
    # if Meta.isexpr(e, :function)
    #     pushfirst!(e.args[2].args, Expr(:macrocall, Symbol("@show"), LineNumberNode(1), :libwasmtime))
    #     return e
    # end

    # if Meta.isexpr(e, :struct, 3) && e.args[2] == :wasmtime_extern
    #     e.args = [
    #         e.args[1],
    #         e.args[2],
    #         [:($(Symbol("_pad", i))::$(t)) for (i, t) in enumerate((UInt8, UInt16, UInt32))]...,
    #         e.args[end],
    #         # :(wasmtime_extern(kind, of) = new(kind, 0, 0, 0, of)),
    #     ]
    #     return e
    # end

    Meta.isexpr(e, :const) || return e

    eq = e.args[1]
    if eq.head === :(=) && eq.args[1] === :WASM_EMPTY_VEC
        e.args[1].args[2] = nothing
    elseif eq.head === :(=) && eq.args[1] === :wasm_name && eq.args[2] === :wasm_byte_vec
        e.args[1].args[2] = :wasm_byte_vec_t
    elseif eq.head === :(=) && eq.args[1] === :wasm_byte_t
        e.args[1].args[2] = :UInt8
    end

    return e
end

function rewrite!(dag::ExprDAG)
    for node in get_nodes(dag)
        for expr in get_exprs(node)
            rewrite!(expr)
        end
    end
end

rewrite!(ctx.dag)

# print
build!(ctx, BUILDSTAGE_PRINTING_ONLY)


text:
module MBOSHI

import ..LOCATIONS
import ..LABELS
import ..NON_SPEECH_UNITS

using StringEncodings

const URL = "https://github.com/besacier/mboshi-french-parallel-corpus.git"

function filter_word(w)
    if w == "</s>"
        return LABELS[:sil]
    elseif w == "<UNK>"
        return LABELS[:unk]
    elseif w == "[silence]"
        return LABELS[:sil]
    else
        return w
    end
end

function filter_pronun(p)
    if p == "."
        return [LABELS[:sil]]
    elseif p == "!!"
        return [LABELS[:unk]]
    else
        return ["$c" for c in p]
    end
end

function filter_phone(p)
    if p == '!' return LABELS[:unk] end
    if p == '.' return LABELS[:sil] end
    p
end

function prepare(datadir)
    mkpath(joinpath(datadir, LOCATIONS[:local]))

    # Download the github directory
    repopath = joinpath(datadir, LOCATIONS[:local], "mboshi-github")
    if ! ispath(repopath)
        @info "Downloading the data, this may take a while..."
        run(`git clone $URL $repopath`)
    else
        @info "Data aleady downloaded"
    end

    # Get the list of WAV files
    allwavdir = joinpath(repopath, "full_corpus_newsplit", "all")
    allwavs = filter(x -> endswith(x, "wav"), readdir(allwavdir))

    # Some WAV files have incorrect time stamp. In this step we create
    # a new "fixed" set of WAV files. This step needs `sox` installed.
    fixedwavdir = mkpath(joinpath(datadir, LOCATIONS[:local], "fixed_wavs"))
    if ! ispath(joinpath(fixedwavdir, ".done"))
        @info "Fixing the time stamp of the WAV files"
        for fname in allwavs
            inpath = joinpath(allwavdir, fname)
            outpath = joinpath(fixedwavdir, fname)
            run(`sox --ignore-length $inpath $outpath`)
        end

        open(joinpath(fixedwavdir, ".done"), "w") do f println(f, today()) end
    else
        @info "WAV files time stamp already fixed"
    end

    # Extract the forced alignments of the data + lexicon
    lexicon = Set()
    phones = Set()
    alifile = joinpath(repopath, "forced_alignments_supervised_spkr",
                       "limsi-align", "trn0.seg")
    full_ali_file = joinpath(datadir, LOCATIONS[:local], "full_ali")
    open(full_ali_file, "w") do f_ali
        open(alifile, enc"latin1", "r") do f
            uttid = ""
            phn_stack = []
            for line in eachline(f)
                if startswith(line, "#@ sid")
                    if uttid ≠ "" println(f_ali) end
                    uttid = split(split(line, "=")[2], ".")[1]
                    print(f_ali, uttid)
                elseif startswith(line, "#@ word")
                    tokens = split(line)
                    phn_stack = [c for c in tokens[4]]
                    word = filter_word(split(tokens[2], "=")[2])
                    pronun = filter_pronun(tokens[4])
                    push!(phones, pronun...)
                    push!(lexicon, (word, pronun))
                elseif ! startswith(line, "#")
                    label = filter_phone(popfirst!(phn_stack))
                    tokens = split(line)
                    s1, s2, s3 = tokens[5:7]
                    duration = parse(Int, s1) + parse(Int, s2) + parse(Int, s3)
                    for i in 1:duration
                        print(f_ali, " ", label)
                    end
                end
            end
        end
    end

    full_ali = Dict()
    open(full_ali_file, "r") do f
        for line in eachline(f)
            tokens = split(line)
            uttid = tokens[1]
            ali = join(tokens[2:end], " ")
            full_ali[uttid] = ali
        end
    end

    ###################################################################
    # Prepare each set

    for set in ["dev", "train"]
        @info "Preparing $set set..."

        setdir = mkpath(joinpath(datadir, set))

        wavdir = joinpath(repopath, "full_corpus_newsplit", set)
        wavs = filter(x -> endswith(x, "wav"), readdir(wavdir))

        # Extract uttids and speakers labels
        uttids = []
        speakers = Set{String}()
        utt2spk = Dict()
        utt2wav = Dict()
        for fname in wavs
            tokens = split(replace(fname, ".wav" => ""), "_")
            speaker = tokens[1]
            uttid = join(tokens[2:end], "_")
            push!(uttids, uttid)
            push!(speakers, speaker)
            utt2spk[uttid] = speaker
            utt2wav[uttid] = abspath(joinpath(fixedwavdir, fname))
        end

        uttids = sort(uttids)
        speakers = sort([s for s in speakers])

        # uttids
        open(joinpath(setdir, LOCATIONS[:uttids]), "w") do f
            for uttid in uttids println(f, uttid) end
        end

        # wav.scp
        open(joinpath(setdir, LOCATIONS[:wavs]), "w") do f
            for uttid in uttids println(f, uttid, "\t", utt2wav[uttid]) end
        end

        # uttids_speakers
        open(joinpath(setdir, LOCATIONS[:utt2spk]), "w") do f
            for uttid in uttids println(f, uttid, "\t", utt2spk[uttid]) end
        end

        # speakers
        open(joinpath(setdir, LOCATIONS[:speakers]), "w") do f
            for s in speakers println(f, s) end
        end

        # ali
        open(joinpath(setdir, LOCATIONS[:ali]), "w") do f
            for uttid in uttids println(f, uttid, "\t", full_ali[uttid]) end
        end

        # trans.wrd
        open(joinpath(setdir, LOCATIONS[:trans] * ".wrd"), "w") do f
            for uttid in uttids
                speaker = utt2spk[uttid]
                fname = "$(join([speaker, uttid], "_")).mb.cleaned"
                open(joinpath(wavdir, fname), "r") do f2
                    println(f, uttid, "\t", readline(f2))
                end
            end
        end

        # trans.wrd
        open(joinpath(setdir, LOCATIONS[:trans] * ".wrd.fr"), "w") do f
            for uttid in uttids
                speaker = utt2spk[uttid]
                fname = "$(join([speaker, uttid], "_")).fr.cleaned.noPunct"
                open(joinpath(wavdir, fname), "r") do f2
                    println(f, uttid, "\t", readline(f2))
                end
            end
        end
    end

    ###################################################################
    # Create the full split
    #
    fulldir = mkpath(joinpath(datadir, "full"))

    tgts = [LOCATIONS[:ali], LOCATIONS[:speakers], LOCATIONS[:trans]*".wrd",
            LOCATIONS[:trans]*".wrd.fr", LOCATIONS[:uttids],
            LOCATIONS[:uttids], LOCATIONS[:utt2spk], LOCATIONS[:wavs]]
    for fname in tgts
        traindir = mkpath(joinpath(datadir, "train"))
        devdir = mkpath(joinpath(datadir, "dev"))
        open(joinpath(fulldir, fname), "w") do f
            run(pipeline(pipeline(`cat $(traindir)/$(fname) $(devdir)/$(fname)`, `sort`, `uniq`), stdout=f))
        end
    end

    ###################################################################
    # Prepare the dictionary and the lexicon
    langdir = mkpath(joinpath(datadir, LOCATIONS[:lang]))

    words = Set()
    for (w, p) in lexicon push!(words, w) end
    open(joinpath(langdir, LOCATIONS[:words]), "w") do f
        for w in sort([w for w in words])
            println(f, w)
        end
    end

    isspeechunit(u) = u ∉ NON_SPEECH_UNITS
    open(joinpath(langdir, LOCATIONS[:units]), "w") do f
        for p in sort([p for p in phones])
            type = isspeechunit(p) ? LABELS[:speechunit] : LABELS[:nonspeeechunit]
            println(f, p, "\t", type)
        end
    end

    lex_entries = []
    for (w,p) in lexicon
        push!(lex_entries, "$w\t$(join(p, " "))")
    end

    open(joinpath(langdir, LOCATIONS[:lexicon]), "w") do f
        for p in sort([p for p in lex_entries])
            println(f, p)
        end
    end
end

end # module

max_issues_repo_name: lucasondel/SpeechDatasets

max_forks_repo_name: lucasondel/SpeechDatasets


text:
module TestLeafOptions

using ComoniconTypes
using ComoniconTargetExpr
using ComoniconTargetExpr: emit_expr, emit_body, emit_norm_body, emit_dash_body
using Test

const test_args = Ref{Vector{Any}}()
const test_kwargs = Ref{Vector{Any}}()

function foo(;kwargs...)
    test_kwargs[] = [kwargs...]
end

cmd = Entry(;
    version=v"1.1.0",
    root=LeafCommand(;
        fn=foo,
        name="leaf",
        options=Dict(
            "option-a" => Option(;sym=:option_a, hint="int", type=Int, short=true),
            "option-b" => Option(;sym=:option_b, hint="float64", type=Float64),
        ),
        flags=Dict(
            "flag-a" => Flag(;sym=:flag_a, short=true),
            "flag-b" => Flag(;sym=:flag_b),
        )
    )
)

eval(emit_expr(cmd))

@testset "test leaf options" begin
    @test command_main(["--option-a=3", "--option-b", "1.2", "-f", "--flag-b"]) == 0
    @test test_kwargs[] == [:option_a=>3, :option_b=>1.2, :flag_a=>true, :flag_b=>true]
    @test command_main(["-o=3", "--option-b", "1.2", "-f", "--flag-b"]) == 0
    @test test_kwargs[] == [:option_a=>3, :option_b=>1.2, :flag_a=>true, :flag_b=>true]
    @test command_main(["-o3", "--option-b", "1.2", "-f", "--flag-b"]) == 0
    @test test_kwargs[] == [:option_a=>3, :option_b=>1.2, :flag_a=>true, :flag_b=>true]
    @test command_main(["--option-a", "--option-b", "1.2", "-f", "--flag-b"]) == 1
    @test command_main(["-o", "--option-b", "1.2", "-f", "--flag-b"]) == 1
end

end


text:
function gen_lattice_hexagonal( a::Float64; coa=8.0/3.0 )
    LL = zeros(3,3)
    LL[:,1] = [1.0, 0.0, 0.0]
    LL[:,2] = [cos(pi/3.0), sin(pi/3.0), 0.0]
    LL[:,3] = [0.0, 0.0, coa]
    return a*LL
end

function gen_lattice_sc( a::Float64 )
    LL = zeros(3,3)
    LL[:,1] = [ a, 0.0, 0.0]
    LL[:,2] = [ 0.0, a, 0.0]
    LL[:,3] = [ 0.0, 0.0, a]
    return LL
end

function gen_lattice_fcc( a::Float64 )
    LL = zeros(3,3)
    LL[:,1] = [-1.0, 0.0, 1.0]
    LL[:,2] = [ 0.0, 1.0, 1.0]
    LL[:,3] = [-1.0, 1.0, 0.0]
    return 0.5*a*LL
end

function gen_lattice_bcc(a)
    LL = zeros(3,3)
    LL[:,1] = [ 1.0,  1.0, 1.0]
    LL[:,2] = [-1.0,  1.0, 1.0]
    LL[:,3] = [-1.0, -1.0, 1.0]
    return 0.5*a*LL
end