text:
module TestManager

using ModelBuilds, HelperFunctions, TestSpecDef, JuMP, BenchmarkTools, Profile, JLD

export Test

function Test(testSpec::TestSpec)
    dataTypeID = testSpec.dataTypeID
    problemNbr = testSpec.problemNbr

	# set default time and number of samples for benchmarking process
	BenchmarkTools.DEFAULT_PARAMETERS.seconds = testSpec.benchMarkSeconds
	BenchmarkTools.DEFAULT_PARAMETERS.samples = testSpec.benchMarkSamples

    dataDict = load(testSpec.dataDict)
    n = dataDict["n"]
    coeff = dataDict["arrayCoeff"]
    bound = dataDict["arrayBound"]

    if dataTypeID in [DAAS20 DS]
        str_con = dataDict["str_con"]
        str_obj = dataDict["str_obj"]
    end

    directMode = testSpec.directMode
    optimization = testSpec.optimization
    writeToLp = testSpec.writeToLpFile
    LpFile = string("Results/LpFiles/" , string(dataTypeID) ,"/", string(dataTypeID) ,"_", problemNbr ,".lp")
    writeToExcel = testSpec.writeToExcel
    excelFile = testSpec.excelFile
    excelLocation = testSpec.excelLocation
    sheetNbr = testSpec.sheetNbr
    writeProfileLog = testSpec.writeProfileLog
    objValue = 0.0
    optimizer = testSpec.optimizer

    println("Start Test")

    if writeToLp == true
        if dataTypeID == A
            model = arrayModel(coeff, bound, n, directMode, optimization)
        elseif dataTypeID == AL
            model = arrayModelLongName(coeff, bound, n, directMode, optimization)
        elseif dataTypeID == ALB
            model = arrayModelLongBasename(coeff, bound, n, directMode, optimization)
        elseif dataTypeID == DAA
            model = denseAxisArrayModel(coeff, bound, n, directMode, optimization)
        elseif dataTypeID == DAAL
            model = denseAxisArrayModelLongName(coeff, bound, n, directMode, optimization)
        elseif dataTypeID == DAALS
            model = denseAxisArrayModelLongNameShortBasename(coeff, bound, n, directMode, optimization)
        elseif dataTypeID == DAALB
            model = denseAxisArrayModelLongBasename(coeff, bound, n, directMode, optimization)
        elseif dataTypeID == DAAS20
            model = denseAxisArrayModelStr20(coeff, bound, str_obj, str_con, directMode, optimization)
        elseif dataTypeID == D
            model = dictModel(coeff, bound, n, directMode, optimization)
        elseif dataTypeID == DS
            model = dictModelStr(coeff, bound, str_obj, str_con, directMode, optimization)
        elseif dataTypeID == DD
            model = dictionaryModel(coeff, bound, n, directMode, optimization)
        elseif dataTypeID == ND
            model = namedDimsModel(coeff, bound, n, directMode, optimization)
        end

        JuMP.write_to_file(model, LpFile)

        if optimization == true
            objValue = JuMP.objective_value(model)
        end

        empty!(model)
        println(".lp file was saved succesfully")
    end

    if writeToExcel == true

        println("benchmarking started")
        #first without saving values

        #save in every case only the second benchmark
        if dataTypeID == A
            @benchmark arrayModel($coeff, $bound, $n, $directMode, $optimization, $optimizer)
            benchmark = @benchmark arrayModel($coeff, $bound, $n, $directMode, $optimization, $optimizer)
        elseif dataTypeID == AL
            @benchmark arrayModelLongName($coeff, $bound, $n, $directMode, $optimization, $optimizer)
            benchmark = @benchmark arrayModelLongName($coeff, $bound, $n, $directMode, $optimization, $optimizer)
        elseif dataTypeID == ALB
            @benchmark arrayModelLongBasename($coeff, $bound, $n, $directMode, $optimization, $optimizer)
            benchmark = @benchmark arrayModelLongBasename($coeff, $bound, $n, $directMode, $optimization, $optimizer)
        elseif dataTypeID == DAA
            @benchmark denseAxisArrayModel($coeff, $bound, $n, $directMode, $optimization, $optimizer)
            benchmark = @benchmark denseAxisArrayModel($coeff, $bound, $n, $directMode, $optimization, $optimizer)
        elseif dataTypeID == DAAL
            @benchmark denseAxisArrayModelLongName($coeff, $bound, $n, $directMode, $optimization, $optimizer)
            benchmark = @benchmark denseAxisArrayModelLongName($coeff, $bound, $n, $directMode, $optimization, $optimizer)
        elseif dataTypeID == DAALS
            @benchmark denseAxisArrayModelLongNameShortBasename($coeff, $bound, $n, $directMode, $optimization, $optimizer)
            benchmark = @benchmark denseAxisArrayModelLongNameShortBasename($coeff, $bound, $n, $directMode, $optimization, $optimizer)
        elseif dataTypeID == DAALB
            @benchmark denseAxisArrayModelLongBasename($coeff, $bound, $n, $directMode, $optimization, $optimizer)
            benchmark = @benchmark denseAxisArrayModelLongBasename($coeff, $bound, $n, $directMode, $optimization, $optimizer)
        elseif dataTypeID == DAAS20
            @benchmark denseAxisArrayModelStr20($coeff, $bound, $str_obj, $str_con, $directMode, $optimization, $optimizer)
            benchmark = @benchmark denseAxisArrayModelStr20($coeff, $bound, $str_obj, $str_con, $directMode, $optimization, $optimizer)
        elseif dataTypeID == D
            @benchmark dictModel($coeff, $bound, $n, $directMode, $optimization, $optimizer)
            benchmark = @benchmark dictModel($coeff, $bound, $n, $directMode, $optimization, $optimizer)
        elseif dataTypeID == DS
            @benchmark dictModelStr($coeff, $bound, $str_obj, $str_con, $directMode, $optimization, $optimizer)
            benchmark = @benchmark dictModelStr($coeff, $bound, $str_obj, $str_con, $directMode, $optimization, $optimizer)
        elseif dataTypeID == DD
            @benchmark dictionaryModel($coeff, $bound, $n, $directMode, $optimization, $optimizer)
            benchmark = @benchmark dictionaryModel($coeff, $bound, $n, $directMode, $optimization, $optimizer)
        elseif dataTypeID == ND
            @benchmark namedDimsModel($coeff, $bound, $n, $directMode, $optimization, $optimizer)
            benchmark = @benchmark namedDimsModel($coeff, $bound, $n, $directMode, $optimization, $optimizer)
        end

        println("benchmarking finished")
		#location specifices where to save the benchmark in the excel file
		#save the benchmark in the excel file
		benchmarkToExcel(excelFile, sheetNbr, benchmark, excelLocation, objValue)
        println("excel file was saved succesfully")
    end
end

end

max_stars_repo_name: SonntagK/JuMPContainerBenchmarks

max_issues_repo_name: SonntagK/JuMPContainerBenchmarks

max_forks_repo_name: SonntagK/JuMPContainerBenchmarks


text:
module TestPlotRecipes

using Test
using Plots
import ThreadPools: LoggedQueuePool, LoggedStaticPool
using ThreadPools

include("util.jl")

@testset "plot recipes" begin

    @testset "log" begin
        plot(readlog("testlog.txt"))
        png("_test.png")
        @test isfile("_test.png")
        rm("_test.png")
    end

    @testset "LoggedQueuePool" begin
        pool = LoggedQueuePool()
        pool.log = readlog("testlog.txt")
        plot(pool)
        close(pool)
        png("_test.png")
        @test isfile("_test.png")
        rm("_test.png")
    end

    @testset "LoggedStaticPool" begin
        pool = LoggedStaticPool(1:Threads.nthreads(), ReentrantLock(), [], time_ns(), readlog("testlog.txt"))
        plot(pool)
        close(pool)
        png("_test.png")
        @test isfile("_test.png")
        rm("_test.png")
    end
end

end # module

ext: jl


text:
# Regular functions for evaluating derivatives and gradients of Chebyshev polynomials

# Regular functions for derivatives

function chebyshev_derivative(weights::AbstractArray{T,N},x::AbstractArray{R,1},pos::S,order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,R<:Number,N,S<:Integer}

  poly = Array{Array{R,2},1}(undef,N)
  @inbounds for i = 1:N
    if i === pos
      poly[i] = derivative_of_chebyshev_polynomial(order[i],normalize_node(x[i],domain[:,i]))
    else
      poly[i] = chebyshev_polynomial(order[i],normalize_node(x[i],domain[:,i]))
    end
  end

  derivative = zero(T)
  @inbounds for i in CartesianIndices(weights)
    poly_product = poly[1][i[1]]
    @inbounds for j = 2:N
      poly_product *= poly[j][i[j]]
    end
    derivative += weights[i]*poly_product
  end

  return derivative*(2.0/(domain[1,pos]-domain[2,pos]))

end

function chebyshev_derivative(weights::AbstractArray{T,N},x::AbstractArray{R,1},pos::S,order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,R<:Number,N,S<:Integer}

  poly = Array{Array{R,2},1}(undef,N)
  @inbounds for i = 1:N
    if i === pos
      poly[i] = derivative_of_chebyshev_polynomial(order,normalize_node(x[i],domain[:,i]))
    else
      poly[i] = chebyshev_polynomial(order,normalize_node(x[i],domain[:,i]))
    end
  end

  derivative = zero(T)
  @inbounds for i in CartesianIndices(weights)
    if sum(Tuple(i)) <= order+N
      poly_product = poly[1][i[1]]
      @inbounds for j = 2:N
        poly_product *= poly[j][i[j]]
      end
      derivative += weights[i]*poly_product
    end
  end

  return derivative*(2.0/(domain[1,pos]-domain[2,pos]))

end

function chebyshev_derivative(cheb_poly::ChebPoly,x::AbstractArray{R,1},pos::S) where {T<:AbstractFloat,R<:Number,S<:Integer}

  derivative = chebyshev_derivative(cheb_poly.weights,x,pos,cheb_poly.order,cheb_poly.domain)

  return derivative

end

function chebyshev_derivative(cheb::ChebInterpRoots,x::AbstractArray{R,1},pos::S) where {T<:AbstractFloat,R<:Number,S<:Integer}

  weights    = chebyshev_weights_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)
  derivative = chebyshev_derivative(weights,x,pos,cheb.order,cheb.domain)

  return derivative

end

function chebyshev_derivative(cheb::ChebInterpExtrema,x::AbstractArray{R,1},pos::S) where {T<:AbstractFloat,R<:Number,S<:Integer}

  weights    = chebyshev_weights_extrema_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)
  derivative = chebyshev_derivative(weights,x,pos,cheb.order,cheb.domain)

  return derivative

end

function chebyshev_derivative(cheb::ChebInterpExtended,x::AbstractArray{R,1},pos::S) where {T<:AbstractFloat,R<:Number,S<:Integer}

  weights    = chebyshev_weights_extended_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)
  derivative = chebyshev_derivative(weights,x,pos,cheb.order,cheb.domain)

  return derivative

end

# Regular functions for gradients

function chebyshev_gradient(weights::AbstractArray{T,N},x::AbstractArray{R,1},order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,R<:Number,N,S<:Integer}

  gradient = Array{R,2}(undef,1,N)

  @inbounds for i = 1:N
    gradient[i] = chebyshev_derivative(weights,x,i,order,domain)
  end

  return gradient

end

function chebyshev_gradient(weights::AbstractArray{T,N},x::AbstractArray{R,1},order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,R<:Number,N,S<:Integer}

  gradient = Array{R,2}(undef,1,N)

  @inbounds for i = 1:N
    gradient[i] = chebyshev_derivative(weights,x,i,order,domain)
  end

  return gradient

end

function chebyshev_gradient(cheb_poly::ChebPoly,x::AbstractArray{R,1}) where {R<:Number}

  gradient = chebyshev_gradient(cheb_poly.weights,x,cheb_poly.order,cheb_poly.domain)

  return gradient

end

function chebyshev_gradient(cheb::ChebInterpRoots,x::AbstractArray{R,1}) where {R<:Number}

  weights  = chebyshev_weights_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)
  gradient = chebyshev_gradient(weights,x,cheb.order,cheb.domain)

  return gradient

end

function chebyshev_gradient(cheb::ChebInterpExtrema,x::AbstractArray{R,1}) where {R<:Number}

  weights  = chebyshev_weights_extrema_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)
  gradient = chebyshev_gradient(weights,x,cheb.order,cheb.domain)

  return gradient

end

function chebyshev_gradient(cheb::ChebInterpExtended,x::AbstractArray{R,1}) where {R<:Number}

  weights  = chebyshev_weights_extended_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)
  gradient = chebyshev_gradient(weights,x,cheb.order,cheb.domain)

  return gradient

end

function chebyshev_gradient(cheb::ChebInterpVertesi,x::AbstractArray{R,1}) where {R<:Number}

  weights  = chebyshev_weights_vertesi_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)
  gradient = chebyshev_gradient(weights,x,cheb.order,cheb.domain)

  return gradient

end

function chebyshev_derivative(weights::AbstractArray{T,N},pos::S,order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T <: AbstractFloat,N,S <: Integer}

  function chebderiv(x::AbstractArray{R,1}) where {R<:Number}

    return chebyshev_derivative(weights,x,pos,order,domain)

  end

  return chebderiv

end

function chebyshev_derivative(weights::AbstractArray{T,N},pos::S,order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T <: AbstractFloat,N,S <: Integer}

  function chebderiv(x::AbstractArray{R,1}) where {R<:Number}

    return chebyshev_derivative(weights,x,pos,order,domain)

  end

  return chebderiv

end

function chebyshev_derivative(cheb_poly::ChebPoly,pos::S) where {S <: Integer}

  function chebderiv(x::AbstractArray{R,1}) where {R<:Number}

    return chebyshev_derivative(cheb_poly,x,pos)

  end

  return chebderiv

end

function chebyshev_derivative(cheb::ChebInterpRoots,pos::S) where {S <: Integer}

  weights = chebyshev_weights_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)

  function chebderiv(x::AbstractArray{R,1}) where {R<:Number}

    return chebyshev_derivative(weights,x,pos,cheb.order,cheb.domain)

  end

  return chebderiv

end

function chebyshev_derivative(cheb::ChebInterpExtrema,pos::S) where {S <: Integer}

  weights = chebyshev_weights_extrema_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)

  function chebderiv(x::AbstractArray{R,1}) where {R<:Number}

    return chebyshev_derivative(weights,x,pos,cheb.order,cheb.domain)

  end

  return chebderiv

end

function chebyshev_derivative(cheb::ChebInterpExtended,pos::S) where {S <: Integer}

  weights = chebyshev_weights_extended_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)

  function chebderiv(x::AbstractArray{R,1}) where {R<:Number}

    return chebyshev_derivative(weights,x,pos,cheb.order,cheb.domain)

  end

  return chebderiv

end

function chebyshev_derivative(cheb::ChebInterpVertesi,pos::S) where {S <: Integer}

  weights = chebyshev_weights_vertesi_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)

  function chebderiv(x::AbstractArray{R,1}) where {R<:Number}

    return chebyshev_derivative(weights,x,pos,cheb.order,cheb.domain)

  end

  return chebderiv

end

function chebyshev_gradient(weights::AbstractArray{T,N},order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T <: AbstractFloat,N,S <: Integer}

  function chebgrad(x::AbstractArray{R,1}) where {R<:Number}

    return chebyshev_gradient(weights,x,pos,order,domain)

  end

  return chebgrad

end

function chebyshev_gradient(weights::AbstractArray{T,N},order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T <: AbstractFloat,N,S <: Integer}

  function chebgrad(x::AbstractArray{R,1}) where {R<:Number}

    return chebyshev_gradient(weights,x,pos,order,domain)

  end

  return chebgrad

end

function chebyshev_gradient(cheb_poly::ChebPoly) where {S <: Integer}

  function chebgrad(x::AbstractArray{R,1}) where {R<:Number}

    return chebyshev_gradient(cheb_poly,x)

  end

  return chebgrad

end

function chebyshev_gradient(cheb::ChebInterpRoots) where {S <: Integer}

  weights = chebyshev_weights_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)

  function chebgrad(x::AbstractArray{R,1}) where {R<:Number}

    return chebyshev_gradient(weights,x,cheb.order,cheb.domain)

  end

  return chebgrad

end

function chebyshev_gradient(cheb::ChebInterpExtrema) where {S <: Integer}

  weights = chebyshev_weights_extrema_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)

  function chebgrad(x::AbstractArray{R,1}) where {R<:Number}

    return chebyshev_gradient(weights,x,cheb.order,cheb.domain)

  end

  return chebgrad

end

function chebyshev_gradient(cheb::ChebInterpExtended) where {S <: Integer}

  weights = chebyshev_weights_extended_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)

  function chebgrad(x::AbstractArray{R,1}) where {R<:Number}

    return chebyshev_gradient(weights,x,cheb.order,cheb.domain)

  end

  return chebgrad

end

function chebyshev_gradient(cheb::ChebInterpVertesi) where {S <: Integer}

  weights = chebyshev_weights_vertesi_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)

  function chebgrad(x::AbstractArray{R,1}) where {R<:Number}

    return chebyshev_gradient(weights,x,cheb.order,cheb.domain)

  end

  return chebgrad

end

max_issues_repo_name: RJDennis/ChebyshevApprox

max_forks_repo_name: RJDennis/ChebyshevApprox


text:
#=##############################################################################
# DESCRIPTION
    Utilities.
# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Sep 2018
  * License   : MIT License
=###############################################################################

"""
`simplewing(b, ar, tr, twist_root, twist_tip, lambda, gamma;
bodytype=RigidWakeBody,
span_NDIVS="automatic", rfl_NDIVS="automatic",
airfoil_root="naca6412.dat", airfoil_tip="naca6412.dat",
airfoil_path=def_rfl_path)`

Generates a symmetric single-section wing.

**ARGUMENTS**
  * `b::Real`         : Span.
  * `ar::Real`        : Aspect ratio defined as b/c_tip.
  * `tr::Real`        : Taper ratio defined as c_tip/c_root.
  * `twist_root::Real`: (deg) twist of the root.
  * `twist_tip::Real` : (deg) twist of the tip.
  * `lambda::Real`    : (deg) sweep.
  * `gamma::Real`     : (deg) dihedral.

**OPTIONAL ARGUMENTS**
  * `bodytype::Type{LBodyTypes}`: Type of lifting body to generate.
  * `span_NDIVS::ndivstype`     : Spanwise divisions.
  * `rfl_NDIVS::ndivstype`    : Chordwise divisions.
  * `airfoil_root::String`      : File to root airfoil contour.
  * `airfoil_tip::String`       : File to tip airfoil contour.
  * `airfoil_path::String`      : Path to airfoil files.

NOTE: See gt.multidscretize for a description of arguments of type `ndivstype`.
NOTE2: In the current implementation, sweep and dihedral are done about the LE.
"""
function simplewing(b::RType, ar::RType, tr::RType, twist_root::RType,
                      twist_tip::RType, lambda::RType, gamma::RType;
                      bodytype::Type{LBodyTypes}=RigidWakeBody,
                      span_NDIVS::ndivstype=nothing,
                      rfl_NDIVS::ndivstype=nothing,
                      airfoil_root::String="naca6412.dat",
                      airfoil_tip::String="naca6412.dat",
                      airfoil_path::String=def_rfl_path,
                      spl_s::Real=0.0000001,
                      rflspl_s::Real=0.00000001,
                      verify_spline::Bool=true,
                      verify_rflspline::Bool=true,
                      opt_args...
                      )

  # ----------------- GEOMETRY DESCRIPTION -------------------------------------
  c_tip = b/ar                        # Tip chord
  c_root = c_tip/tr                   # Root chord
  semispan = b/2                      # (m) semi-span length

  y_tip = b/2
  x_tip = y_tip*tan(lambda*pi/180)
  z_tip = y_tip*tan(gamma*pi/180)

  chords = [0.00 c_root/semispan;     # (semi-span position, chord c/semib)
            1.00 c_tip/semispan]

  twists = [0.0 twist_root;           # (semi-span position, twist (deg))
            1.0 twist_tip]

  x_pos = [0.00 0;                    # (semi-span position, LE x-position x/semib)
           1.00 x_tip/semispan]

  z_pos = [0.00 0;                    # (semi-span position, LE z-position x/semib)
           1.00 z_tip/semispan]

  airfoil_files = [(0.0, airfoil_root), # (semi-span position, airfoil file)
                   (1.0, airfoil_tip)]

 # ----------------- DISCRETIZATION 0000----------------------------------------

  # Defines divisions
  if span_NDIVS==nothing
    b_NDIVS = [(1.0, 35, 20.0, true)]  # Span cell sections
  else
    b_NDIVS = span_NDIVS
  end

  if rfl_NDIVS==nothing
    urfl_NDIVS = [(0.25, 7,   10.0, false),   # Cells on upper side of airfoils
                  (0.50,  5,    1.0, true),
                  (0.25,  6, 1/10.0, false)]
  else
    urfl_NDIVS = rfl_NDIVS
  end

  lrfl_NDIVS = urfl_NDIVS             # Cells on lower side of airfoils

  # ----------------- LOFTING PARAMETERS ---------------------------------------
  b_low = -1.0                        # Lower bound of span lofting
  b_up = 1.0                          # Upper bound of span lofting
  symmetric = true                    # Lofting symmetric about b=0
  spl_k = 1                           # Spline order of distributions along span
  # spl_s = 0.0000001                 # Spline smoothing of distribution along span
  # rflspl_s = 0.00000001             # Spline smoothing of airfoil cross sections.
  # verify_spline = false             # Plots the splined distributions
  # verify_rflspline = true           # Plots the splined airfoil cross sections

  return generate_loft_liftbody(bodytype, airfoil_files, airfoil_path,
                                        urfl_NDIVS, lrfl_NDIVS,
                                        semispan, b_low, b_up, b_NDIVS,
                                        chords, twists, x_pos, z_pos;
                                        dimsplit=1,
                                        symmetric=symmetric,
                                        spl_k=spl_k, spl_s=spl_s,
                                        verify_spline=verify_spline,
                                        verify_rflspline=verify_rflspline,
                                        rflspl_s=rflspl_s,
                                        opt_args...
                                    )
end


text:
using ContextualBandits, AdaptivePartition
using Test
using Random: seed!
using Distributions: Normal

@testset "Contexts" begin
    include("contexts/test_UniformContext.jl")
end

@testset "Policies" begin
    include("policies/test_FixedPolicy.jl")
    include("policies/test_PartitionTS.jl")
    include("policies/test_KNNUCB.jl")
end

@testset "Rewards" begin
    include("rewards/test_GasussianRewards.jl")
end

@testset "Simulate" begin
    include("test_simulate.jl")
end


text:
export DirectMaxwellSolver

# To do:
#
# We can create a similar type that reduces the number of factorization.  During the Newton
# iteration for solving the nonlinear equation, the matrix changes only diagonally.
# Furthermore, this diagonal change is parametrized by ω.  Therefore, we can perform
# factorization only once at a single ω, and use the idea of fast frequency sweep for other
# ω's.  Because γ(ω) has a pole, we may need to use the Padé approximation instead of the
# simple power series expansion, for which γ(ω) will generate infinitely many terms because
# γ(ω) is infinitely differentiable.
#
# A simiar idea cannot be applied to solving the lasing equation, where the diagonal change
# involves a change in the hole-burning term and is not parametrized by the scalar ω.
# However, we may construct a new parametrized matrix for each hole-burning term, by fixing
# the hole burning term and changing only ω.  The downside of this is that the coefficients
# of the series expansion of the solution needs to be calculated for every ω.  (For the
# nonlasing equation, this needed to be calculated only once for varying ω during the Newton
# iteration...  Is this right?  This is correct when both A and b are parametrized by ω, but
# in our case b is probably not, so this may not be correct.)  However, if the major
# performance improvement comes from the reduced number of factorizations rather than the
# reduced number of linear solves, this scheme should still have a huge benefit.
#
# However, all these ideas become irrelevant for 3D problems, where linear systems will be
# solved iteratively and therefore no factorization will be performed.  Hence, implementing
# the above ideas is not a priority.

mutable struct DirectMaxwellSolver{VR<:AbsVecFloat,VC<:AbsVecComplexF,MC<:AbsMatComplexF,FC<:Factorization{ComplexF},Kₑ₋₁} <: AbstractMaxwellSolver
    CC::MC  # double curl operator ∇ × μ⁻¹ ∇ ×
    A::MC  # Maxwell operator
    Me2c::MC  # operator interpolating edge Ew for w = 1, ..., Kₑ at voxel corners
    Mc2e::NTuple{Kₑ₋₁,MC}  # operator interpolating voxel-corner Ew for w = 1, ..., Kₑ at Ev for v that is circularly shifted 1, ..., Kₑ
    vtemp::VC  # temporary storage for complex vector
    wtemp::VC  # temporary storage for complex vector
    fact::FC  # matrix factors of A

    function DirectMaxwellSolver{VR,VC,MC,FC,Kₑ₋₁}(CC::AbsMatNumber, A::AbsMatNumber,
                                                 Me2c::AbsMatNumber, Mc2e::NTuple{Kₑ₋₁,AbsMatNumber},
                                                 vtemp::AbsVecNumber, wtemp::AbsVecNumber
                                                 ) where {VR<:AbsVecFloat,VC<:AbsVecComplexF,MC<:AbsMatComplexF,FC<:Factorization{ComplexF},Kₑ₋₁}
        isequal(size(CC)...) || @error "CC should be square, but size(CC) = $(size(CC))."

        size(CC)==size(A)==size(Me2c) || @error "size.((CC,A,Me2c)) = $(size.((CC,A,Me2c))) should be all same."

        for k = 1:Kₑ₋₁
            size(Mc2e[k])==size(CC) || @error "size(Mc2e[$k]) = $(size(Mc2e[k])) and size(CC) = $(size(CC)) should be the same."
        end

        size(CC,1)==length(vtemp)==length(wtemp) ||
            @error "size.((CC,vtemp,wtemp),1) = $(size.((CC,vtemp,wtemp),1)) should be all same."

        return new(CC, A, Me2c, Mc2e, vtemp, wtemp)
    end
end

# Below, CC is shared between modes, but A is not (it is overwritten for new modes).
# Therefore we need to keep a separate copy of A.
function DirectMaxwellSolver(CC::MC, Me2c::MC, Mc2e::Tuple{Vararg{MC}}; Vtype::Type=Vector) where {MC<:AbsMatComplexF}
    Kₑ₋₁ = length(Mc2e)
    FC = Base.return_types(lu, tuple(MC))[1];
    n = size(CC, 2)

    return DirectMaxwellSolver{Vtype{Float},Vtype{ComplexF},MC,FC,Kₑ₋₁}(CC, similar(CC), Me2c, Mc2e, VecComplexF(undef,n), VecComplexF(undef,n))
end

Base.similar(mxsolver::DirectMaxwellSolver) = DirectMaxwellSolver(mxsolver.CC, mxsolver.Me2c, mxsolver.Mc2e)
Base.size(mxsolver::DirectMaxwellSolver) = size(mxsolver.CC)

# Factorize the matrix A and cache it.
function SALTBase.init_mxsolver!(mxsolver::DirectMaxwellSolver, ω::Number, ε::AbsVecNumber)
    mA = size(mxsolver.A, 1)
    nε = length(ε)
    mA==nε || @error "size(mxsolver.A,1) = $mA and length(ε) = $nε should be same."

    mxsolver.A .= mxsolver.CC  # initialize; works for sparse matrices with same nonzero entry pattern
    for i = 1:nε
        mxsolver.A[i,i] -= ω^2 * ε[i]
    end

    t = @elapsed mxsolver.fact = lu(mxsolver.A)
    # @info "time factorization: $t"

    return nothing
end

function SALTBase.mxsolve!(x::VC, mxsolver::DirectMaxwellSolver{VR,VC}, b::VC) where {VR,VC}
    t = @elapsed ldiv!(x, mxsolver.fact, b)
    # @info "time mxsolve!: $t"

    return nothing
end

function SALTBase.mxsolve_transpose!(x::VC, mxsolver::DirectMaxwellSolver{VR,VC}, b::VC) where {VR,VC}
    t = @elapsed ldiv!(x, transpose(mxsolver.fact), b)
    # @info "time mxsolve_transpose!: $t"

    return nothing
end

function SALTBase.mxapply!(b::VC, mxsolver::DirectMaxwellSolver{VR,VC}, x::VC) where {VR,VC}
    mul!(b, mxsolver.A, x)

    return nothing
end

function SALTBase.abs2_interp!(abs2ψ::VR, mxsolver::DirectMaxwellSolver{VR,VC}, ψ::VC) where {VR,VC}
    # Calculate the amplitude squares of Ex, Ey, Ez at their Yee grid locations.
    abs2ψ .= abs2.(ψ)

    # Interpolate Ex, Ey, Ez at the Yee grid voxel corner locations.
    ψc = mxsolver.vtemp
    mul!(ψc, mxsolver.Me2c, ψ)

    # Calculate the amplitude squares of Ex, Ey, Ez at the Yee grid's Ez, Ex, Ey locations.
    ψe = mxsolver.wtemp
    for M = mxsolver.Mc2e
        mul!(ψe, M, ψc)
        abs2ψ .+= abs2.(ψe)
    end
end

# In IterativeMaxwellData, the only things I need to store would be the things that I need
# for applying the Maxwell operation or averaging operation matrix-free, such as ε, ω, and
# grid.  The functions themselves that apply the matrix-free Maxwell operation or
# matrix-free averaging operation to a column vector should be implemented in MaxwellFDFD.
# mxsolve!, etc that take IterativeMaxwellData simply need to call those functions with the
# arguments it stores.


text:
#=
function Quantum2dSquareSimulator{T}(lattice::Quantum2dSquareLattice{T}, dimM::Int, wholeiteration::Int)
function (simulator::Quantum2dSquareSimulator)()
=#
type Quantum2dSquareSimulator{T} <: Quantum2dSimulator{T}
	dimM::Int
	wholeiteration::Int
	normalizationFactor::Array{T}
	data4energy::Array{T}
	countiteration::Int
	trottercount::Int
	lattice::Quantum2dSquareLattice{T}

	function Quantum2dSquareSimulator{T}(lattice::Quantum2dSquareLattice{T}, dimM::Int)
		this = new{T}()
		this.dimM = dimM
		this.wholeiteration = getTrotteriteration(lattice) * 3
		this.lattice = lattice
		this.normalizationFactor = Array{T}(this.wholeiteration)
		this.data4energy = Array{T}(this.wholeiteration)
		initializeCount!(this)
		return this
	end
end

Quantum2dSquareSimulator{T}(lattice::Quantum2dSquareLattice{T}, dimM::Int) =
Quantum2dSquareSimulator{T}(lattice,dimM)

#=
<info>
tensorT[l,r,f,b,u,v]

<The steps in simulator>
renormalizex(tensorT1, tensorT2; isTensorT::Int = 1)
renormalizey and renormalizez uses renormalizeX but with permutedims on input and output tensors

	getTensorU(tensorT, dimM)
	getTensorV is from getTensorU with permutedims
		getTenmatMMd(tensorT)  -> tenmatM
		svd(tenmatM)
		truncate
		Ul, lambdaVector, Uld = svd(matrixMMd)
		trunUl = truncMatrixU(Ul,dimM)
		tensorUl = matU2tenU(trunUl)

	getNewTensorT(tensorT1, tensorT2, tensorU, tensorV)

	setTensorT

	return renoramalizationFactor

=#

function (simulator::Quantum2dSquareSimulator)()
	initializeCount!(simulator)
	while true
		whichIsTensorT = 1
		countUp!(simulator)
		renormalizeX!(simulator; whichIsTensorT = whichIsTensorT)
		countUp!(simulator)
		renormalizeY!(simulator; whichIsTensorT = whichIsTensorT)
		countUp!(simulator,"trotter")
		renormalizeZ!(simulator; whichIsTensorT = whichIsTensorT)
		whichIsTensorT = 2
		countUp!(simulator)
		renormalizeX!(simulator; whichIsTensorT = whichIsTensorT)
		countUp!(simulator)
		renormalizeY!(simulator; whichIsTensorT = whichIsTensorT)
		countUp!(simulator,"trotter")
		renormalizeZ!(simulator; whichIsTensorT = whichIsTensorT)
		if isDone(simulator)
			break
		end
	end
	freeenergy = getFreeEnergy(simulator)
	magnetization = getExpectationValue(simulator)
	if isPottsModel(simulator)
		NOS = convert(T,getStates(simulator))
		magnetization = (NOS * magnetization - one(T))/(NOS - one(T))
	end

	return freeenergy, magnetization
end

#---
# functions in simulator
#=
function renormalize(tensorT1, tensorT2, dimM; whichIsTensorT = 1)
return newTensorT, newTensorTtilde, normalizationFactor
	in simulate_2dQ_renormalize.jl file
=#
function renormalizeZ!{T}(simulator::Quantum2dSquareSimulator{T}; whichIsTensorT::Int = 1)
	# time axis
	tensorT1, tensorT2 = getTensorT(simulator, whichIsTensorT)
	newTensorT, newTensorTtilde, normalizationFactor = renormalize(tensorT1, tensorT2, getDimM(simulator); whichIsTensorT = whichIsTensorT)

	setTensorT!(simulator.lattice, newTensorT)
	setTensorT!(simulator.lattice, newTensorTtilde; tilde = true)
	setNormalizationfactor!(simulator, normalizationFactor)
end

function renormalizeX!{T}(simulator::Quantum2dSquareSimulator{T}; whichIsTensorT::Int = 1)
	tensorT1, tensorT2 = getTensorT(simulator, whichIsTensorT)
	newTensorT, newTensorTtilde, normalizationFactor = renormalize(permutedims(tensorT1,[3,4,5,6,1,2]), permutedims(tensorT2,[3,4,5,6,1,2]), getDimM(simulator); whichIsTensorT = whichIsTensorT)

	setTensorT!(simulator.lattice, permutedims(newTensorT,[5,6,1,2,3,4]))
	setTensorT!(simulator.lattice, permutedims(newTensorTtilde,[5,6,1,2,3,4]); tilde = true)
	setNormalizationfactor!(simulator, normalizationFactor)
end

function renormalizeY!{T}(simulator::Quantum2dSquareSimulator{T}; whichIsTensorT::Int = 1)
	tensorT1, tensorT2 = getTensorT(simulator, whichIsTensorT)
	newTensorT, newTensorTtilde, normalizationFactor = renormalize(permutedims(tensorT1,[5,6,1,2,3,4]), permutedims(tensorT2,[5,6,1,2,3,4]), getDimM(simulator); whichIsTensorT = whichIsTensorT)

	setTensorT!(simulator.lattice, permutedims(newTensorT,[3,4,5,6,1,2]))
	setTensorT!(simulator.lattice, permutedims(newTensorTtilde,[3,4,5,6,1,2]); tilde = true)
	setNormalizationfactor!(simulator, normalizationFactor)
end

#---
# test Function
function testQuantum2dSquareSimulator()
	externalfield = 1.0e-13
	temperature = 0.0
	spininfo = SpinInfo("quantum_ising_2",externalfield, temperature)
	trotterparameter = 0.01
	trotteriteration = 2
	trotterinfo = TrotterInfo(trotterparameter, trotteriteration)
	testModel = QuantumIsingModel(spininfo, trotterinfo)
	lattice = buildLattice("quantum_2d_square", testModel)
	dimM = 3
	simulator = Quantum2dSquareSimulator(lattice, dimM)
	simulator()
	return simulator
end


text:
"""
Maximize sharp coefficient allocation.
"""
function max_sharpe(Σ, r̄, rf)
    one = ones(size(r̄, 1))
    invΣ = pinv(Σ, 1E-25)
    v = invΣ * (r̄ .- one * rf)
    return v
end

"""
Equal weights.
"""
function equal_weights(numA)
    return fill(1 / numA, numA)
end

max_stars_repo_name: andrewrosemberg/PortfolioOpt

max_issues_repo_name: andrewrosemberg/PortfolioOpt

max_forks_repo_name: andrewrosemberg/PortfolioOpt


text:
_isbound(var) = !isnothing(var._op)

_getid(var::Ghost.Variable) = var.id
_getid(call::Ghost.Call) = call.id
_getid(x) = x

_gettapeop(var::Ghost.Variable) =
    _isbound(var) ? var._op :
                    error("Cannot get operation on tape for unbound variable.")

_gettapeval(input::Ghost.Input) = input.val
_gettapeval(call::Ghost.Call) = call.val
_gettapeval(var::Ghost.Variable) =
    _isbound(var) ? _gettapeval(var._op) :
                    error("Cannot get value on tape for unbound variable.")
_gettapeval(x) = x

_issimulatable(op) = false
_issimulatable(var::Ghost.Variable) =
    (_gettapeval(var) isa Base.Broadcast.Broadcasted) ? _issimulatable(_gettapeop(var)) :
                                                        _gettapeval(var) isa SBitstreamLike
_issimulatable(input::Ghost.Input) = input.val isa SBitstreamLike
_issimulatable(call::Ghost.Call) = any(_issimulatable, call.args)

_isbcast(::typeof(Base.broadcasted)) = true
_isbcast(x) = false

_isinput(::Ghost.Input) = true
_isinput(x::Ghost.Variable) = _isinput(_gettapeop(x))
_isinput(x) = false

_isvariable(x::Ghost.AbstractOp) = true
_isvariable(x::Ghost.Variable) = true
_isvariable(x) = false

_isstruct(f::T) where T = isstructtype(T) && !(f isa Function)


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

"""
Traverse the entries in a tree and its subtrees in post or pre order.

Function parameter should have following signature:

    (Cstring, Ptr{Void}, Ptr{Void}) -> Cint
"""
function treewalk(f::Function, tree::GitTree, payload=Any[], post::Bool = false)
    cbf = cfunction(f, Cint, (Cstring, Ptr{Void}, Ptr{Void}))
    cbf_payload = Ref{typeof(payload)}(payload)
    @check ccall((:git_tree_walk, :libgit2), Cint,
                  (Ptr{Void}, Cint, Ptr{Void}, Ptr{Void}),
                   tree.ptr, post, cbf, cbf_payload)
    return cbf_payload
end

function filename(te::GitTreeEntry)
    str = ccall((:git_tree_entry_name, :libgit2), Cstring, (Ptr{Void},), te.ptr)
    str != C_NULL && return unsafe_string(str)
    return nothing
end

function filemode(te::GitTreeEntry)
    return ccall((:git_tree_entry_filemode, :libgit2), Cint, (Ptr{Void},), te.ptr)
end

function object(repo::GitRepo, te::GitTreeEntry)
    obj_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_tree_entry_to_object, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Ref{Void}),
                   obj_ptr_ptr, repo.ptr, te.ptr)
    return GitAnyObject(obj_ptr_ptr[])
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using ParallelDifferentialEvolution
using Test
import Distributed: pmap

@testset "ParallelDifferentialEvolution.jl" begin
    function bounds(x)
        minx = -10.0
        maxx = 10.0
        x .* (maxx .- minx) .+ minx
    end
    f(x) = sum(bounds(x).^2) / length(x)
    res = diffevo(f, 1)
    @test isapprox(bounds(res.m)[1], 0.0, atol=1e-6)
    @test isapprox(res.fm, 0.0, atol=1e-13)
    res = diffevo(f, 1, cb=logcb)
    @test isapprox(bounds(res.m)[1], 0.0, atol=1e-6)
    @test isapprox(res.fm, 0.0, atol=1e-13)
    res = diffevo(f, 3)
    @test all(isapprox.(bounds(res.m), 0.0, atol=1e-6))
    @test isapprox(res.fm, 0.0, atol=1e-13)
    res = diffevo(f, 10, maxiter=2000)
    @test all(isapprox.(bounds(res.m), 0.0, atol=1e-6))
    @test isapprox(res.fm, 0.0, atol=1e-13)
    res = diffevo(f, 10, maxiter=2000, fmap=pmap)
    @test all(isapprox.(bounds(res.m), 0.0, atol=1e-6))
    @test isapprox(res.fm, 0.0, atol=1e-13)
    res = diffevo(f, 3, seeds=[[0.2,0.1,0.9], [0.7,0.2,0.8]])
    @test all(isapprox.(bounds(res.m), 0.0, atol=1e-6))
    @test isapprox(res.fm, 0.0, atol=1e-13)
end


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module CommReLU

using Random
using ExtractMacro
using ..Interface
using ..Common

using ..DeltaE.ArraySets

export GraphCommReLU

import ..Interface: energy, delta_energy, update_cache!, neighbors

function gen_ξ(N::Integer, P::Integer)
    ξ = BitArray(undef, P, N)
    ξv = Vector{BitVector}(undef, P)
    for a = 1:P
        ξv[a] = bitrand(N)
        ξ[a,:] = ξv[a]
    end
    # the ξv representation is for energy
    # the ξ representation is for delta_energy and update_cache!
    y = bitrand(P)
    return ξ, ξv, y
end

mutable struct Stabilities
    K2::Int
    p1::Vector{ArraySet}
    m1::Vector{ArraySet}
    p2::ArraySet
    m2::ArraySet
    Δ1s::Vector{Vector{Int}}
    Δ2s::Vector{Int}
    ξsi::BitVector
    last_move::Int
    function Stabilities(K2::Integer, P::Integer)
        p1 = [ArraySet(P) for k in 1:K2]
        m1 = [ArraySet(P) for k in 1:K2]
        p2 = ArraySet(P)
        m2 = ArraySet(P)
        Δ1s = [zeros(P) for k in 1:K2]
        Δ2s = zeros(P)
        ξsi = BitArray(undef, P)
        return new(K2, p1, m1, p2, m2, Δ1s, Δ2s, ξsi, 0)
    end
end

mutable struct GraphCommReLU <: SimpleGraph{Int}
    N::Int
    K1::Int
    K2::Int
    P::Int
    ξ::BitMatrix
    ξv::Vector{BitVector}
    y::BitVector
    stab::Stabilities
    tmps::BitVector
    function GraphCommReLU(K2::Int, ξ::BitMatrix, ξv::Vector{BitVector}, y::BitVector)
        P, N = size(ξ)
        @assert N % K2 == 0
        K1 = N ÷ K2
        @assert length(ξv) == P
        @assert all(length(ξ1) == N for ξ1 in ξv)
        @assert length(y) == P
        iseven(K1) || throw(ArgumentError("K1 must be even, given: $K1"))
        iseven(K2) || throw(ArgumentError("K2 must be even, given: $K2"))
        stab = Stabilities(K2, P)
        tmps = BitVector(undef, N)
        return new(N, K1, K2, P, ξ, ξv, y, stab, tmps)
    end
end

"""
    GraphCommReLU(K1::Integer, K2::Integer, P::Integer; fc = false) <: SimpleGraph{Int}

A `SimpleGraph` implementing a two-layer binary committee machine with `K2` hidden units,
where each hidden unit has `K1` binary (\$±1\$) synapses, trained on `P` random i.i.d. \$±1\$
patterns. The output of each hidden unit is computed through a ReLU function \$max(0,x)\$,
and then half of the `K2` units have a negative weight associated to their output.
Both `K1` and `K2` must be even.
The default is to generate a tree-like machine, but it can be made fully-connected with the
argument `fc=true`.

The energy of the model is computed as the number of misclassified patterns.
"""
function GraphCommReLU(K1::Integer, K2::Integer, P::Integer; fc::Bool = false)
    Kin = fc ? K1 : K1 * K2
    ξ, ξv, y = gen_ξ(Kin, P)
    if fc
        ξ = repeat(ξ, outer=(1,K2))
        ξv = [repeat(ξ1, K2) for ξ1 in ξv]
    end
    GraphCommReLU(K2, ξ, ξv, y)
end

function Base.empty!(stab::Stabilities)
    @extract stab : p1 m1 p2 m2 Δ1s Δ2s ξsi
    map!(empty!, p1, p1)
    map!(empty!, m1, m1)
    empty!(p2)
    empty!(m2)
    map!(v->fill!(v, 0), Δ1s, Δ1s)
    fill!(Δ2s, 0)
    stab.last_move = 0
    return stab
end

function energy(X::GraphCommReLU, C::Config)
    @assert X.N == C.N
    @extract C : s
    @extract X : N K1 K2 P ξv y stab tmps
    @extract stab : p1 m1 p2 m2 Δ1s Δ2s

    E = 0
    empty!(stab)
    tmps = BitArray(undef, K1)

    for a = 1:P
        Δ2 = 0
        for k = 1:K2
            c = 2 * (2k ≤ K2) - 1
            unsafe_copyto!(tmps, 1, s, (k-1)*K1 + 1, K1)
            map!(⊻, tmps, tmps, ξv[a][(k-1)*K1 .+ (1:K1)])
            Δ1 = K1 - 2 * sum(tmps)
            Δ1s[k][a] = Δ1
            Δ2 += c * max(Δ1, 0)
            if Δ1 > 0
                push!(p1[k], a)
            elseif Δ1 == 0
                push!(m1[k], a)
            end
        end
        Δ2 *= 2y[a] - 1
        Δ2s[a] = Δ2
        if Δ2 == 2
            push!(p2, a)
        elseif Δ2 ≤ 0
            Δ2 == 0 && push!(m2, a)
            E += 1
        end
    end

    return E
end

function update_cache!(X::GraphCommReLU, C::Config, move::Int)
    @assert X.N == C.N
    @extract C : s
    @extract X : N K1 K2 P ξ y stab
    @extract stab : p1 m1 p2 m2 Δ1s Δ2s ξsi last_move
    @assert N == length(s)

    si = s[move]
    k2 = (move - 1) ÷ K1 + 1
    # k1 = (move - 1) % K1 + 1
    last_move ≠ move && unsafe_copyto!(ξsi, 1, ξ, (move-1)*P + 1, P)
    stab.last_move = move
    si && flipbits!(ξsi)
    Δ1k = Δ1s[k2]
    p1k = p1[k2]
    m1k = m1[k2]
    c = 2 * (2k2 ≤ K2) - 1
    # @assert ξsi == ξ[:,move] ⊻ si
    @inbounds for a = 1:P
        o = 2y[a] - 1
        oldΔ2 = Δ2s[a] * o
        oldΔ1 = Δ1k[a]
        newΔ1 = oldΔ1 + (2 - 4 * ξsi[a])

        newΔ2 = oldΔ2
        if oldΔ1 ≤ 0 && newΔ1 > 0
            @assert oldΔ1 == 0 && newΔ1 == 2
            delete!(m1k, a)
            push!(p1k, a)
            newΔ2 += 2c
        elseif oldΔ1 > 0 && newΔ1 ≤ 0
            @assert oldΔ1 == 2 && newΔ1 == 0
            delete!(p1k, a)
            push!(m1k, a)
            newΔ2 -= 2c
        elseif oldΔ1 > 0 && newΔ1 > 0
            newΔ2 += c * (newΔ1 - oldΔ1)
        elseif oldΔ1 < 0 && newΔ1 == 0
            push!(m1k, a)
        elseif oldΔ1 == 0
            @assert newΔ1 < 0
            delete!(m1k, a)
        end
        Δ1k[a] = newΔ1

        newΔ2 == oldΔ2 && continue
        oldΔ2 *= o
        newΔ2 *= o

        if oldΔ2 > 2 && newΔ2 == 2
            push!(p2, a)
        elseif oldΔ2 == 2
            delete!(p2, a)
            newΔ2 == 0 && push!(m2, a)
        elseif oldΔ2 == 0
            delete!(m2, a)
            newΔ2 == 2 && push!(p2, a)
        elseif oldΔ2 < 0 && newΔ2 == 0
            push!(m2, a)
        end
        Δ2s[a] = newΔ2
    end
    si && flipbits!(ξsi)

    return
end

function delta_energy_naive(X::GraphCommReLU, C::Config, move::Int)
    @assert C.N == X.N
    @extract C : s
    @extract X : N P ξ stab
    @assert 1 ≤ move ≤ N

    spinflip!(C, move)
    E1 = energy(X, C)
    spinflip!(C, move)
    E0 = energy(X, C)
    return E1 - E0
end

function delta_energy(X::GraphCommReLU, C::Config, move::Int)
    @assert C.N == X.N
    @extract C : s
    @extract X : N K1 K2 P ξ y stab
    @extract stab : p1 m1 p2 m2 ξsi last_move
    @assert 1 ≤ move ≤ N

    si = s[move]
    k2 = (move - 1) ÷ K1 + 1
    # k1 = (move - 1) % K1 + 1
    c = 2 * (2k2 ≤ K2) - 1

    last_move ≠ move && unsafe_copyto!(ξsi, 1, ξ, (move-1)*P + 1, P)
    stab.last_move = move
    si && flipbits!(ξsi)
    ΔE = 0
    @inbounds for a in p2
        o = 2y[a] - 1
        if c == o && a ∈ p1[k2]
            ΔE += ~ξsi[a]
        elseif c ≠ o && (a ∈ m1[k2] || a ∈ p1[k2])
            ΔE += ξsi[a]
        end
    end
    @inbounds for a in m2
        o = 2y[a] - 1
        if c == o && (a ∈ m1[k2] || a ∈ p1[k2])
            ΔE -= ξsi[a]
        elseif c ≠ o && a ∈ p1[k2]
            ΔE -= ~ξsi[a]
        end
    end
    si && flipbits!(ξsi)
    # td = delta_energy_naive(X, C, move)
    # ΔE == td || @show ΔE, td
    # @assert ΔE == td
    return ΔE
end

function check_delta(X::GraphCommReLU, C::Config, move::Int)
    e0 = energy(X, C)
    delta = delta_energy(X, C, move)
    spinflip!(C, move)
    e1 = energy(X, C)
    spinflip!(C, move)

    (e1-e0) == delta || (@show e1,e0,delta,e1-e0; error())
end

neighbors(X::GraphCommReLU, i::Int) = AllButOne(X.N, i)

end


text:
@testset "$(format("Sampling: UniformModel[VecOrMat]"))" begin
    #==================== Vector Case ====================#
    let
        tc1 = rand(10000)
        output = Sampler(UniformModel(0.5)).draw(tc1)
        @test count(x->ismissing(x),output) > 1000
        @test count(x->!ismissing(x),output) > 1000
        # stronger
        tc2 = rand(10000)
        output = Sampler(UniformModel(0.1)).draw(tc2)
        @test count(x->ismissing(x),output) >1000
        @test count(x->!ismissing(x),output) <1000
    end

    #==================== Matrix Case ====================#
    let
        @test 4 <= count(x -> !ismissing(x),Sampler(UniformModel(0.1)).draw(ones(10,10))) <= 10
        @test 10 < count(x -> ismissing(x), Sampler(UniformModel(0.1)).draw(ones(10,10)))
    end
    #==================== Factory Mode ====================#
    let
        sampler = provide(Sampler{UniformModel}(),rate = 0.5)
        tc1 = rand(10000)
        output = sampler.draw(tc1)
        @test count(x->ismissing(x),output) > 1000
        @test count(x->!ismissing(x),output) > 1000
        # stronger
        tc2 = rand(10000)
        sampler2 = provide(Sampler{UniformModel}(),rate = 0.1)
        output = sampler2.draw(tc2)
        @test count(x->ismissing(x),output) >1000
        @test count(x->!ismissing(x),output) <=1000
        @test 4 <= count(x -> !ismissing(x),sampler2.draw(ones(10,10))) <= 10
        @test 10 < count(x -> ismissing(x),sampler2.draw(ones(10,10)))
    end
end


text:
#PS4 Dayuan Xie

# using Pkg
# Pkg.add("Distributions")

using Optim
using HTTP
using GLM
using LinearAlgebra
using Random
using Statistics
using DataFrames
using CSV
using FreqTables
using Distributions

function q6()
# Question 1

url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS4-mixture/nlsw88t.csv"
df = CSV.read(HTTP.get(url).body)
X = [df.age df.white df.collgrad]
Z = hcat(df.elnwage1, df.elnwage2, df.elnwage3, df.elnwage4,
df.elnwage5, df.elnwage6, df.elnwage7, df.elnwage8)
y = df.occ_code

function mlogit_with_Z(theta, X, Z, y)

    alpha = theta[1:end-1]
    gamma = theta[end]
    K = size(X,2)
    J = length(unique(y))
    N = length(y)
    bigY = zeros(N,J)
    for j=1:J
        bigY[:,j] = y.==j
    end
    bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

    T = promote_type(eltype(X),eltype(theta))
    num   = zeros(T,N,J)
    dem   = zeros(T,N)
    for j=1:J
        num[:,j] = exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)
        dem .+= num[:,j]
    end

    P = num./repeat(dem,1,J)

    loglike = -sum( bigY.*log.(P) )

    return loglike
end
startvals = [2*rand(7*size(X,2)).-1; .1]
truevals = [0.040374454953787886, 0.24399409314825987, -1.5713213321281576,
            0.04332547647972951, 0.14685584442757782, -2.959104052831784,
            0.10205742855969692, 0.7473089603290985, -4.120051101278431,
            0.037562869245200775, 0.6884901936101743, -3.6557712276056433,
            0.020454321679879015, -0.3584005467776664, -4.376931619035733,
            0.10746368376194072, -0.5263737776237152, -6.199194024108684,
            0.11688248621330447, -0.28705521539111295, -5.322249141580821, 1.307479115935751]
#  startvals = truevals.*rand(size(truevals))
td = TwiceDifferentiable(theta -> mlogit_with_Z(theta, X, Z, y), startvals; autodiff = :forward)
theta_hat_optim_ad = optimize(td, truevals, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
theta_hat_mle_ad = theta_hat_optim_ad.minimizer
H  = Optim.hessian!(td, theta_hat_mle_ad)
theta_hat_mle_ad_se = sqrt.(diag(inv(H)))
println(theta_hat_mle_ad)
println(theta_hat_mle_ad_se)
println([theta_hat_mle_ad theta_hat_mle_ad_se])

# Question 2
println("Yes.The Estimated coefficient γ is 1.307479115935751.")

# Question 3
function lgwt(N0::Integer,a::Real=-1,b::Real=1)
    N  = N0-1
    N1 = N+1
    N2 = N+2
    xu = range(-1,stop=1,length=N1)
    y = cos.((2*(0:N) .+ 1)*pi/(2*N .+ 2))  .+  ( 0.27/N1 ) .* sin.( pi .* xu .* N/N2 )
    L  = zeros(N1,N2)
    Lp = zeros(N1,N2)
    y0 = 2
    vareps = 2e-52
    i = 0
    while norm(y.-y0,Inf)>vareps
        L[:,1]  .= 1
        Lp[:,1] .= 0
        L[:,2] .= y
        for k=2:N1
            L[:,k+1] = ( (2*k-1)*y .* L[:,k] .- (k-1)*L[:,k-1] )/k
        end
        Lp = (N2)*( L[:,N1] .- y .* L[:,N2] )./(1 .- y.^2)
        y0 = y
        y  = y0 - L[:,N2]./Lp
        i+=1
    end

    x = (a.*(1 .- y) .+ b .* (1 .+ y))./2
    w=(b-a)./((1 .- y.^2).*Lp.^2)*(N2/N1)^2
    return x,w
end

# Question 3a
d = Normal(0,1)
nodes, weights = lgwt(7,-4,4)
println(sum(weights.*pdf.(d,nodes)))
println(sum(weights.*nodes.*pdf.(d,nodes)))

# Question 3b
d = Normal(0,2)
nodes, weights = lgwt(7,-5*sqrt(2),5*sqrt(2))
println(sum(weights.*nodes.*nodes.*pdf.(d,nodes)))

nodes, weights = lgwt(10,-5*sqrt(2),5*sqrt(2))
println(sum(weights.*nodes.*nodes.*pdf.(d,nodes)))

# Estimated value, 3.925601931584661 and 3.9769753072954845
# True value σ^2=4

# Question 3c
D=1000000
μ=0
σ=sqrt(2)
d = Normal(μ,σ^2)
x=rand(Uniform(-5*σ,5*σ),D)
println(sum(((5*σ-(-5*σ))/D).*x.*x.*pdf.(d,x)))
println(sum(((5*σ-(-5*σ))/D).*x.*pdf.(d,x)))
println(sum(((5*σ-(-5*σ))/D).*pdf.(d,x)))

D=1000
x=rand(Uniform(-5*σ,5*σ),D)
println(sum(((5*σ-(-5*σ))/D).*x.*x.*pdf.(d,x)))
println(sum(((5*σ-(-5*σ))/D).*x.*pdf.(d,x)))
println(sum(((5*σ-(-5*σ))/D).*pdf.(d,x)))

# Question 4
mu = theta_hat_mle_ad[end]
sigma = theta_hat_mle_ad_se[end]
function mixed_mlogit_with_Z(theta, X, Z, y)

    alpha = theta[1:end-1]
    gamma = theta[end]
    K = size(X,2)
    J = length(unique(y))
    N = length(y)
    bigY = zeros(N,J)
    for j=1:J
        bigY[:,j] = y.==j
    end
    bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

    T = promote_type(eltype(X),eltype(theta))
    num   = zeros(T,N,J)
    dem   = zeros(T,N)
    for j=1:J
        num[:,j] = exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)
        dem .+= num[:,j]
    end

    P = num./repeat(dem,1,J)

    d = Normal(mu,sigma^2)
    nodes, weights = lgwt(7,-4*sigma,4*sigma)
    mixed_p=sum(weights.*P.*pdf.(d,nodes))
    mixed_loglike = -sum( bigY.*log.(mixed_p) )

    return mixed_loglike
end

startvals = [2*rand(7*size(X,2)).-1; .1]
td = TwiceDifferentiable(theta ->mixed_loglike, startvals; autodiff = :forward)
mixed_theta_hat_optim_ad = optimize(dt,startvals, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
mixed_theta_hat_mle_ad = mixed_theta_hat_optim_ad.minimizer
H  = Optim.hessian!(td, mixed_theta_hat_mle_ad)
mixed_theta_hat_mle_ad_se = sqrt.(diag(inv(H)))
println(mixed_theta_hat_mle_ad)
println([mixed_theta_hat_mle_ad mixed_theta_hat_mle_ad_se])

# Question 5
mu = theta_hat_mle_ad[end]
sigma = theta_hat_mle_ad_se[end]
function mc_mixed_mlogit_with_Z(theta, X, Z, y)

    alpha = theta[1:end-1]
    gamma = theta[end]
    K = size(X,2)
    J = length(unique(y))
    N = length(y)
    bigY = zeros(N,J)
    for j=1:J
        bigY[:,j] = y.==j
    end
    bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

    T = promote_type(eltype(X),eltype(theta))
    num   = zeros(T,N,J)
    dem   = zeros(T,N)
    for j=1:J
        num[:,j] = exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)
        dem .+= num[:,j]
    end

    P = num./repeat(dem,1,J)

    d = Normal(mu,sigma^2)
    weights = (4*sigma-(-4*sigma))/size(P)
    mc_mixed_p=sum(weights.*P.*pdf.(d,P))
    mc_mixed_loglike = -sum( bigY.*log.(mc_mixed_p) )

    return mc_mixed_loglike
end

startvals = [2*rand(7*size(X,2)).-1; .1]
td = TwiceDifferentiable(theta ->mc_mixed_loglike, startvals; autodiff = :forward)
mc_mixed_theta_hat_optim_ad = optimize(dt,startvals, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
mc_mixed_theta_hat_mle_ad = mc_mixed_theta_hat_optim_ad.minimizer
H  = Optim.hessian!(td, mc_mixed_theta_hat_mle_ad)
mc_mixed_theta_hat_mle_ad_se = sqrt.(diag(inv(H)))
println(mc_mixed_theta_hat_mle_ad)
println([mc_mixed_theta_hat_mle_ad mc_mixed_theta_hat_mle_ad_se])

end

# Question 6
q6()


text:
using Test

include("_dummy_model.jl")

@testset "machines" begin
    include("machines.jl")
end

@testset "controls" begin
    include("controls.jl")
end


text:
# Distilled from JuMP's test/JuMPExtension.jl

"""
    ModelView{T,S}

A wrapper around any JuMP `model` to prefix variables and constraints with
`<component.class>::` and store them in an `component.objects`-Dictionary instead.
Every other operation is just handed through to the inner `model`.
"""
struct ModelView{T,S} <: JuMP.AbstractModel
    jumpmodel::T
    component::S
end

prefix(view::ModelView, name::String) = issimple(view.component) ? name : string(naming(view.component), "::", name)

Base.broadcastable(view::ModelView) = Ref(view)
JuMP.object_dictionary(view::ModelView) = get!(model(view.component).jumpobjects, naming(view.component), Dict{Symbol}{Any}())

JuMP.variable_type(view::ModelView) = variable_type(view.jumpmodel)
JuMP.add_variable(view::ModelView, v::JuMP.AbstractVariable, name::String="") =
    model(view.component).jumpnames ? add_variable(view.jumpmodel, v, prefix(view, name)) : add_variable(view.jumpmodel, v)
JuMP.delete(view::ModelView, variable_ref::VariableRef) = delete(view.jumpmodel, variable_ref)
JuMP.is_valid(view::ModelView, variable_ref::VariableRef) = is_valid(view.jumpmodel, variable_ref)
JuMP.num_variables(view::ModelView) = num_variables(view.jumpmodel)

JuMP.constraint_type(view::ModelView) = constraint_type(view.jumpmodel)
JuMP.add_constraint(view::ModelView, c::JuMP.AbstractConstraint, name::String="") =
    model(view.component).jumpnames ? add_constraint(view.jumpmodel, c, prefix(view, name)) : add_constraint(view.jumpmodel, c)
JuMP.delete(view::ModelView, constraint_ref::ConstraintRef) = delete(view.jumpmodel, constraint_ref)
JuMP.is_valid(view::ModelView, constraint_ref::ConstraintRef) = is_valid(view.jumpmodel, constraint_ref)
JuMP.num_constraints(view::ModelView) = num_constraints(view.jumpmodel)

# Objective
JuMP.set_objective(view::ModelView, sense::MOI.OptimizationSense, f::JuMP.AbstractJuMPScalar) = set_objective(view.jumpmodel, sense, f)
JuMP.set_objective_sense(view::ModelView, sense) = set_objective_sense(view.jumpmodel, sense)
JuMP.objective_sense(view::ModelView) = objective_sense(view.jumpmodel)

JuMP.objective_function_type(view::ModelView) = objective_function_type(view.jumpmodel)
JuMP.objective_function(view::ModelView) = objective_function(view.jumpmodel)

# Names
JuMP.variable_by_name(view::ModelView, name::String) =
    model(view.component).jumpnames ? variable_by_name(view.jumpmodel, prefix(view, name)) : nothing
JuMP.constraint_by_name(view::ModelView, name::String) =
    model(view.component).jumpnames ? constraint_by_name(view.jumpmodel, prefix(view, name)) : nothing

# Show
JuMP.show_backend_summary(io::IO, view::ModelView) = show_backend_summary(io, view.jumpmodel)
JuMP.show_objective_function_summary(io::IO, view::ModelView) = show_objective_function_sumary(io, view.jumpmodel)
JuMP.objective_function_string(print_mode, view::ModelView) = objective_function_string(print_mode, view.jumpmodel)
JuMP.show_constraints_summary(io::IO, view::ModelView) = show_constraints_summary(io, view.jumpmodel)
JuMP.constraints_string(print_mode, view::ModelView) = constraints_string(print_mode, view.jumpmodel)

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
@testset "VecRep" begin
    @testset "Intersection" begin
        h = HalfSpace([1, 0], 1) ∩ HyperPlane([0, 1], 1)
        @test nhyperplanes(h) == 1
        @test nhalfspaces(h) == 1
        hs = h.halfspaces
        hp = h.hyperplanes
        _hp = h.hyperplanes.hyperplanes
        intersect!(h, HalfSpace([1, 1], 2))
        @test h.halfspaces === hs
        @test h.hyperplanes === hp
        @test h.hyperplanes.hyperplanes === _hp
        @test nhyperplanes(h) == 1
        @test nhalfspaces(h) == 2
        intersect!(h, HyperPlane([1, 1], 3))
        @test h.halfspaces === hs
        @test h.hyperplanes === hp
        @test h.hyperplanes.hyperplanes === _hp
        @test nhyperplanes(h) == 2
        @test nhalfspaces(h) == 2
    end
end


text:
module LDrawParser

using LightGraphs
using GraphUtils
using GeometryBasics
using Rotations, CoordinateTransformations
using Parameters
using Logging
using StaticArrays

export
    get_part_library_dir,
    set_part_library_dir!,
    find_part_file,

    SubModelPlan,
    DATModel,
    SubFileRef,
    BuildingStep,
    MPDModel,

    model_name,
    has_model,
    get_model,
    has_part,
    get_part,

    parse_ldraw_file!,
    parse_ldraw_file,

    populate_part_geometry!,
    ldraw_base_transform

global PART_LIBRARY_DIR = "/scratch/ldraw_parts_library/ldraw/"
get_part_library_dir() = deepcopy(PART_LIBRARY_DIR)
function set_part_library_dir!(path)
    PART_LIBRARY_DIR = path
end

function find_part_file(name,library=get_part_library_dir())
    if isfile(name)
        return name
    end
    for (root,dirs,_) in walkdir(library)
        for dir in dirs
            d = joinpath(root,dir)
            for p in [joinpath(d,name),joinpath(d,lowercase(name))]
                if isfile(p)
                    return p
                end
            end
        end
    end
    println("Part file ",name," not found in library at ",library)
    return nothing
end

"""
    FILE_TYPE

All LDraw files carry the LDR (default), DAT or MPD extension.

Official Parts
Part | Subpart | Primitive | 8_Primitive | 48_Primitive | Shortcut
Unofficial Parts
Unofficial_Part| Unofficial_Subpart | Unofficial_Primitive |
Unofficial_8_Primitive | Unofficial_48_Primitive | Unofficial_Shortcut

The file type is usually prefaced in one of the following ways
    0 !LDRAW_ORG <type> (qualifier(s)) (update-tag)
    0 LDRAW_ORG <type> update-tag
    0 Official LCAD <type> update-tag
    0 Unofficial <type>
    0 Un-official <type>
"""
@enum FILE_TYPE begin
    MODEL
    PART
    SUBPART
    PRIMITIVE
    SHORTCUT
    NONE_FILE_TYPE
end

const FILE_TYPE_HEADERS = Set{String}([
    "!LDRAW_ORG",
    "LDRAW_ORG",
    "OFFICIAL LCAD",
    "OFFICIAL",
    "UNOFFICIAL",
    "UN-OFFICIAL"
])
const FILE_TYPE_DICT = Dict{String,FILE_TYPE}(
    "MODEL"                     => MODEL,
    "PART"                      => PART,
    "UNOFFICIAL_PART"           => PART,
    "SUBPART"                   => SUBPART,
    "UNOFFICIAL_SUBPART"        => SUBPART,
    "PRIMITIVE"                 => PRIMITIVE,
    "8_PRIMITIVE"               => PRIMITIVE,
    "48_PRIMITIVE"              => PRIMITIVE,
    "UNOFFICIAL_PRIMITIVE"      => PRIMITIVE,
    "UNOFFICIAL_8_PRIMITIVE"    => PRIMITIVE,
    "UNOFFICIAL_48_PRIMITIVE"   => PRIMITIVE,
    "UNOFFICIAL_SHORTCUT"       => PRIMITIVE,
    "SHORTCUT"                  => SHORTCUT,
)
parse_file_type(k::AbstractString) = get(FILE_TYPE_DICT,uppercase(k), NONE_FILE_TYPE)

"""
    COMMAND_CODE

The line type of a line is the first number on the line. The line types are:
0. META          # 0 !<META command> <additional parameters>
1. SUB_FILE_REF  # 1 <colour> x y z a b c d e f g h i <file>
2. LINE          # 2 <colour> x1 y1 z1 x2 y2 z2
3. TRIANGLE      # 3 <colour> x1 y1 z1 x2 y2 z2 x3 y3 z3
4. QUADRILATERAL # 4 <colour> x1 y1 z1 x2 y2 z2 x3 y3 z3 x4 y4 z4
5. OPTIONAL_LINE # 5 <colour> x1 y1 z1 x2 y2 z2 x3 y3 z3 x4 y4 z4
"""
@enum COMMAND_CODE begin
    META          = 0 # 0 !<META command> <additional parameters>
    SUB_FILE_REF  = 1 # 1 <colour> x y z a b c d e f g h i <file>
    LINE          = 2 # 2 <colour> x1 y1 z1 x2 y2 z2
    TRIANGLE      = 3 # 3 <colour> x1 y1 z1 x2 y2 z2 x3 y3 z3
    QUADRILATERAL = 4 # 4 <colour> x1 y1 z1 x2 y2 z2 x3 y3 z3 x4 y4 z4
    OPTIONAL_LINE = 5 # 5 <colour> x1 y1 z1 x2 y2 z2 x3 y3 z3 x4 y4 z4
end
parse_command_code(k::AbstractString) = COMMAND_CODE(parse(Int,k))

"""
    META_COMMAND

0 !<META command> <additional parameters>
- `!` is used to positively identify this as a META command. (Note: A few
    official meta commands do not start with a ! in order to preserve backwards compatibility, however, all new official META commands must start with a ! and it is strongly recommended that new unofficial meta-commands also start with a !)
- `<META command>` is any string in all caps
- `<additional parameters>` is any string. Note that if a META command does not
    require any additional parameter, none should be given.
"""
@enum META_COMMAND begin
    FILE
    NAME
    STEP
    ROTSTEP_BEGIN # ROTSTEP (<x-angle> <y-angle> <z-angle> [(REL| ADD | ABS)] | END)
    ROTSTEP_END
    FILE_TYPE_DECLARATION
    OTHER_META_COMMAND
end

const META_COMMAND_DICT = Dict{String,META_COMMAND}(
    "FILE"          =>FILE,
    "STEP"          =>STEP,
    "ROTSTEP"       =>ROTSTEP_BEGIN,
    "ROTSTEP END"   =>ROTSTEP_END,
    "NAME"          =>NAME,
    "NAME:"         =>NAME,
    [k=>FILE_TYPE_DECLARATION for k in FILE_TYPE_HEADERS]...
)
parse_meta_command(k::AbstractString) = get(META_COMMAND_DICT, uppercase(k), OTHER_META_COMMAND)

@enum ROTATION_MODE begin
    REL
    ADD
    ABS
end
const ROTATION_MODE_DICT = Dict{String,ROTATION_MODE}(
    "REL"=>REL,
    "ADD"=>ADD,
    "ABS"=>ABS,
)
parse_rotation_mode(k::AbstractString) = get(ROTATION_MODE_DICT,uppercase(k),REL)

function parse_line(line)
    if Base.Sys.isunix()
        line = replace(line,"\\"=>"/") # switch from Windows directory delimiters to Unix
    end
    split_line = split(line)
    return split_line
end
const SplitLine = Vector{A} where {A<:AbstractString}
for op in [:parse_file_type, :parse_command_code]
    @eval $op(split_line::SplitLine) = $op(split_line[1])
end

function parse_meta_command(split_line::SplitLine)
    if split_line[1] == "0"
        return parse_meta_command(split_line[2:end])
    end
    cmd = parse_meta_command(split_line[1])
    if cmd == ROTSTEP_BEGIN # NOTE not necessary, as it only affects the VIEW (not the model)
        if parse_meta_command(join(split_line[1:2]," ")) == ROTSTEP_END
            cmd = ROTSTEP_END
        end
    end
    return cmd
end

const Point3D   = Point{3,Float64}

"""
    NgonElement

Represents geometry from an LDraw file
"""
struct NgonElement{E}
    color::Int
    geom::E
end

"""
    OptionalLineElement

Represents optional line geometry from an LDraw file
"""
struct OptionalLineElement
    color::Int
    geom::Line
    control_pts::Line
end

"""
    SubFileRef

Represents a sub-file reference from an LDraw file. Encodes the placement of a
part or submodel.
"""
struct SubFileRef
    color::Int
    pos::Point3D
    rot::Mat{3,3,Float64}
    file::String
end
Base.summary(s::SubFileRef) = string("SubFileRef → ",s.file," : ",s.pos)
model_name(r::SubFileRef) = r.file
function build_transform(ref::SubFileRef)
    # @warn "Need to account for the weird LDraw coordinate system at some point"
    Translation(ref.pos[1],ref.pos[2],ref.pos[3]) ∘ LinearMap(ref.rot)
end
function SubFileRef(ref::SubFileRef,T::AffineMap)
    new_ref = SubFileRef(
        ref.color,
        Point3D(T.translation...),
        T.linear,
        ref.file
    )
end

"""
    BuildingStep

Represents a sequence of part placements that make up a building step in a LDraw
file.
"""
struct BuildingStep
    parent::String # points to the SubModelPlan to which it belongs
    lines::Vector{SubFileRef}
    BuildingStep(p::String) = new(p,Vector{SubFileRef}())
end
Base.push!(step::BuildingStep,ref::SubFileRef) = push!(step.lines,ref)
n_lines(s::BuildingStep) = length(s.lines)

"""
    SubModelPlan

Represents the sequence of building steps that make up a sub model in an LDraw
file
"""
struct SubModelPlan
    name::String
    steps::Vector{BuildingStep}
    SubModelPlan(name::String) = new(
        name,
        Vector{BuildingStep}([BuildingStep(name)])
        )
end
model_name(r::SubModelPlan) = r.name
n_build_steps(m::SubModelPlan) = length(m.steps)
n_components(m::SubModelPlan) = sum(map(n_lines,m.steps))
BuildingStep(p::SubModelPlan) = BuildingStep(model_name(p))
Base.summary(n::SubModelPlan) = string("SubModelPlan: ",model_name(n),": ",
    n_build_steps(n)," building steps, ",n_components(n)," components")

const Quadrilateral{Dim,T} = GeometryBasics.Ngon{Dim,T,4,Point{Dim,T}}
mutable struct Toggle
    status::Bool
end
function set_status!(t::Toggle,val=true)
    t.status = val
end
get_status(t::Toggle) = copy(t.status)

"""
    DATModel

Encodes the raw geometry of a LDraw part stored in a .dat file. It is possible
to avoid populating the geometry fields, which is useful for large models or
models that use parts from the LDraw library.
"""
struct DATModel
    name::String
    line_geometry::Vector{NgonElement{Line{3,Float64}}}
    triangle_geometry::Vector{NgonElement{Triangle{3,Float64}}}
    quadrilateral_geometry::Vector{NgonElement{Quadrilateral{3,Float64}}}
    optional_line_geometry::Vector{OptionalLineElement}
    subfiles::Vector{SubFileRef} # points to other DATModels
    populated::Toggle
    DATModel(name::String) = new(
        name,
        Vector{NgonElement{Line{3,Float64}}}(),
        Vector{NgonElement{Triangle{3,Float64}}}(),
        Vector{NgonElement{Quadrilateral{3,Float64}}}(),
        Vector{OptionalLineElement}(),
        Vector{String}(),
        Toggle(false)
    )
end
model_name(r::DATModel) = r.name
function extract_surface_geometry(m::DATModel)
    elements = Vector{GeometryBasics.Ngon}()
    for vec in (m.triangle_geometry,m.quadrilateral_geometry)
        for e in vec
            push!(elements,e.geom)
        end
    end
    elements
end
function extract_geometry(m::DATModel)
    elements = Vector{GeometryBasics.Ngon}()
    for vec in (
            m.line_geometry,
            m.triangle_geometry,
            m.quadrilateral_geometry,
            m.optional_line_geometry
        )
        for e in vec
            push!(elements,e.geom)
        end
    end
    elements
end
function extract_points(m::DATModel)
    pts = Vector{Point3{Float64}}()
    for e in extract_geometry(m)
        for pt in e.points
            push!(pts,pt)
        end
    end
    pts
end
function incorporate_geometry!(m::DATModel,ref::SubFileRef,child::DATModel,scale=1.0)
    t = build_transform(ref)
    incorporate_geometry!(m,child,t,scale)
end
function incorporate_geometry!(m::DATModel,child::DATModel,t,scale=1.0)
    # @info "incorporating geometry from $(child.name) into $(m.name)"
    for (parent_geometry,child_geometry) in zip(
            (m.line_geometry, m.triangle_geometry, m.quadrilateral_geometry,
                m.optional_line_geometry),
            (child.line_geometry, child.triangle_geometry, child.quadrilateral_geometry,
                child.optional_line_geometry)
        )
        for element in child_geometry
            push!(parent_geometry, t(element)*scale)
        end
    end
    return m
end
function DATModel(part::DATModel,T,scale=1.0)
    new_part = DATModel(part.name)
    incorporate_geometry!(new_part,part,T,scale)
    set_status!(new_part.populated,get_status(part.populated))
    new_part
end

"""
    MPDModel

The MPD model stores the information contained in a .mpd or .ldr file. This
includes a submodel tree (stored implicitly in a dictionary that maps model_name
to SubModelPlan) and a part list. The first model in MPDModel.models is the main
model. All the following are submodels of that model and/or each other.
"""
@with_kw_noshow struct MPDModel #<: AbstractCustomNEDiGraph{CustomNode{Union{SubModelPlan,DATModel},String},CustomEdge{SubModelRef,String},String}
    models      ::Dict{String,SubModelPlan} = Dict{String,SubModelPlan}() # each file is a list of steps
    parts       ::Dict{String,DATModel} = Dict{String,DATModel}()
    sub_parts   ::Dict{String,DATModel} = Dict{String,DATModel}()
end
has_part(m::MPDModel,k) = haskey(m.parts,k) || haskey(m.sub_parts,k)
function get_part(m::MPDModel,k)
    if haskey(m.parts,k)
        return m.parts[k]
    elseif haskey(m.sub_parts,k)
        return m.sub_parts[k]
    end
    return nothing
end
has_model(m::MPDModel,k) = haskey(m.models,k)
get_model(m::MPDModel,k) = get(m.models,k,nothing)
function set_part!(m::MPDModel,part,k=part.name)
    if haskey(m.sub_parts,k)
        m.sub_parts[k] = part
    else
        m.parts[k] = part
    end
    part
end
function add_part!(m::MPDModel,k)
    @assert !has_part(m,k)
    m.parts[k] = DATModel(k)
end
function add_sub_part!(m::MPDModel,k)
    @assert !has_part(m,k)
    m.sub_parts[k] = DATModel(k)
end
part_keys(m::MPDModel) = keys(m.parts)
sub_part_keys(m::MPDModel) = keys(m.sub_parts)
all_part_keys(m::MPDModel) = Base.Iterators.flatten((part_keys(m),sub_part_keys(m)))

points_to_part(m::MPDModel,ref::SubFileRef) = has_part(m,ref.file)
points_to_model(m::MPDModel,ref::SubFileRef) = has_model(m,ref.file)

# struct LDRGeometry
#     lines::Vector{NgonElement{Line{3,Float64}}}
#     triangles::Vector{NgonElement{Triangle{3,Float64}}}
#     quadrilaterals::Vector{NgonElement{Quadrilateral{3,Float64}}}
#     optional_lines::Vector{OptionalLineElement}
#     populated::Toggle
#     LDRGeometry() = new(
#         Vector{NgonElement{Line{3,Float64}}}(),
#         Vector{NgonElement{Triangle{3,Float64}}}(),
#         Vector{NgonElement{Quadrilateral{3,Float64}}}(),
#         Vector{OptionalLineElement}(),
#         Toggle(false)
#     )
# end
#
# struct LDRModel
#     name::String
#     file_type::FILE_TYPE
#     geometry::LDRGeometry
#     building_steps::Vector{BuildingStep}
#     sub_file_refs::Vector{SubFileRef}
# end

@with_kw_noshow mutable struct MPDModelState
    file_type::FILE_TYPE = NONE_FILE_TYPE
    active_model::String = ""
    active_part::String = ""
    rot_stack::Vector{SMatrix{3,3,Float64}} = [one(SMatrix{3,3,Float64})]
end

Base.summary(s::MPDModelState) = string(
"MPDModelState(file_type=$(s.file_type), active_model=$(s.active_model), active_part=$(s.active_part))")

update_state(state::MPDModelState) = MPDModelState(state) # TODO deal with single step macro commands, etc.
function active_building_step(submodel::SubModelPlan,state)
    @assert !isempty(submodel.steps)
    active_step = submodel.steps[end]
end
function active_submodel(model::MPDModel,state)
    @assert !isempty(model.models)
    return model.models[state.active_model]
end
function set_active_model!(model::MPDModel,state,name)
    # @assert !haskey(model.models,name)
    if !haskey(model.models,name)
        model.models[name] = SubModelPlan(name)
    else
        @debug "$name is already in model!"
    end
    return MPDModelState(state,active_model=name)
end
function active_building_step(model::MPDModel,state)
    active_model = active_submodel(model,state)
    return active_building_step(active_model,state)
end
function set_new_active_building_step!(model::SubModelPlan)
    push!(model.steps,BuildingStep(model))
    return model
end
function set_new_active_building_step!(model::MPDModel,state)
    active_model = active_submodel(model,state)
    set_new_active_building_step!(active_model)
    return model
end
function active_part(model::MPDModel,state)
    # @assert !isempty(model.parts)
    @assert has_part(model,state.active_part)
    # return model.parts[state.active_part]
    return get_part(model,state.active_part)
end
function set_active_part!(model::MPDModel,state,name)
    # @assert !has_part(model,name) "$name is already a part in model!"
    if !has_part(model,name)
        add_part!(model,name)
    else
        @debug "$name is already a part in model!"
    end
    @debug "Active part = $name"
    return MPDModelState(state,active_part=name)
end
function add_sub_file_placement!(model::MPDModel,state,ref)
    # TODO figure out how to place a subfile that is not part of a build step,
    # but is rather (presumably) a subfile of a .dat model
    if state.file_type == MODEL
        if !isempty(state.active_model) # != ""
            push!(active_building_step(model,state),ref)
        end
    else
        if !isempty(state.active_part) # != ""
            push!(active_part(model,state).subfiles,ref)
        end
    end
    if !has_part(model,ref.file)
        if isempty(state.active_part)
            add_part!(model,ref.file)
        else
            add_sub_part!(model,ref.file)
        end
    end
    return state
end

"""
    parse_ldraw_file!

Args:
    - model
    - filename or IO
"""
function parse_ldraw_file!(model,io,state = MPDModelState())
    # state = MPDModelState()
    for line in eachline(io)
        # try
            if length(line) == 0
                continue
            end
            split_line = parse_line(line)
            if isempty(split_line[1])
                continue
            end
            code = parse_command_code(split_line)
            @debug "LINE: $line"
            @debug "code: $code"
            if code == META
                state = read_meta_line!(model,state,split_line)
            elseif code == SUB_FILE_REF
                state = read_sub_file_ref!(model,state,split_line)
            # Geometry
            elseif code == LINE
                state = read_line!(model,state,split_line)
            elseif code == TRIANGLE
                state = read_triangle!(model,state,split_line)
            elseif code == QUADRILATERAL
                state = read_quadrilateral!(model,state,split_line)
            elseif code == OPTIONAL_LINE
                state = read_optional_line!(model,state,split_line)
            end
        # catch e
        #     bt = catch_backtrace()
        #     showerror(stdout,e,bt)
        #     rethrow(e)
        # end
    end
    return model
end
function parse_ldraw_file!(model,filename::String,args...)
    open(find_part_file(filename),"r") do io
        parse_ldraw_file!(model,io,args...)
    end
end
parse_ldraw_file(io,args...) = parse_ldraw_file!(MPDModel(),io,args...)
parse_color(c) = parse(Int,c)

"""
    read_meta_line(model,state,line)

Modifies the model and parser_state based on a META command. For example, the
FILE meta command indicates the beginning of a new file, so this creates a new
active model into which subsequent building steps will be placed.
The STEP meta command indicates the end of the current step, which prompts the
parser to close the current build step and begin a new one.
"""
function read_meta_line!(model,state,line)
    @assert parse_command_code(line[1]) == META
    if length(line) < 2
        @debug "Returning because length(line) < 2. Usually this means the end of the file"
        return state
    end
    # cmd = line[2]
    cmd = parse_meta_command(line[2:end])
    @debug "cmd: $cmd"
    if cmd == FILE || cmd == NAME
        filename = join(line[3:end]," ")
        ext = lowercase(splitext(filename)[2])
        if ext == ".dat"
            state = set_active_part!(model,state,filename)
            if state.file_type == NONE_FILE_TYPE
                state.file_type = PART
            end
        elseif ext == ".mpd" || ext == ".ldr"
            state = set_active_model!(model,state,filename)
            if state.file_type == NONE_FILE_TYPE
                state.file_type = MODEL
            end
        end
        @debug "file = $filename"
    elseif cmd == STEP
        set_new_active_building_step!(model,state)
    # elseif cmd == ROTSTEP_BEGIN
    #     rx,ry,rz = [parse(Float64,line[3]),parse(Float64,line[4]),parse(Float64,line[5])]
    #     rot_mode = parse_rotation_mode(line[6])
    #     rmat = RotMatrix(RotXYZ(rx,ry,rz)).mat
    #     if rot_mode == ABS
    #         empty!(state.rot_stack)
    #         push!(state.rot_stack, rmat)
    #     elseif rot_mode == REL
    #         @warn "Need to be sure that RotXYZ does things in the right order"
    #         push!(state.rot_stack, rmat)
    #     elseif rot_mode == ADD
    #         push!(state.rot_stack, rmat)
    #     end
    # elseif cmd == ROTSTEP_END
    #     if length(state.rot_stack) > 1
    #         pop!(state.rot_stack) # remove the last
    #     else
    #         state.rot_stack[1] = one(eltype(state.rot_stack))
    #     end
    elseif cmd == FILE_TYPE_DECLARATION
        state.file_type = parse_file_type(line[3])
        if state.file_type == NONE_FILE_TYPE
            @debug "file type not resolved on line : $line"
        end
        @debug "file_type=$(state.file_type)"
    else
        # TODO Handle other META commands, especially BFC
    end
    return state
end

"""
    read_sub_file_ref

Receives a SUB_FILE_REF line (with the leading SUB_FILE_REF id stripped)
"""
function read_sub_file_ref!(model,state,line)
    @assert parse_command_code(line[1]) == SUB_FILE_REF
    @assert length(line) >= 15 "$line"
    color = parse_color(line[2])
    # coordinate of part
    x,y,z = parse.(Float64,line[3:5])
    # rotation of part
    rot_mat = collect(transpose(reshape(parse.(Float64,line[6:14]),3,3)))
    file = join(line[15:end]," ")
    # TODO add a line struct to the model
    ref = SubFileRef(
        color,
        Point3D(x,y,z),
        Mat{3,3,Float64}(rot_mat),
        file
    )
    add_sub_file_placement!(model,state,ref)
    # push!(model.sub_file_refs,ref)
    return state
end

"""
    read_line!

For reading lines of type LINE
"""
function read_line!(model,state,line)
    @assert parse_command_code(line[1]) == LINE
    @assert length(line) == 8 "$line"
    color = parse_color(line[2])
    p1 = Point3D(parse.(Float64,line[3:5]))
    p2 = Point3D(parse.(Float64,line[6:8]))
    # add to model
    push!(
        active_part(model,state).line_geometry,
        NgonElement(color,Line(p1,p2))
        )
    return state
end

"""
    read_triangle!

For reading lines of type TRIANGLE
"""
function read_triangle!(model,state,line)
    @assert parse_command_code(line[1]) == TRIANGLE
    @assert length(line) == 11 "$line"
    color = parse_color(line[2])
    p1 = Point3D(parse.(Float64,line[3:5]))
    p2 = Point3D(parse.(Float64,line[6:8]))
    p3 = Point3D(parse.(Float64,line[9:11]))
    # add to model
    push!(
        active_part(model,state).triangle_geometry,
        NgonElement(color,Triangle(p1,p2,p3))
        )
    return state
end

"""
    read_quadrilateral!

For reading lines of type QUADRILATERAL
"""
function read_quadrilateral!(model,state,line)
    @assert parse_command_code(line[1]) == QUADRILATERAL
    @assert length(line) == 14 "$line"
    color = parse_color(line[2])
    p1 = Point3D(parse.(Float64,line[3:5]))
    p2 = Point3D(parse.(Float64,line[6:8]))
    p3 = Point3D(parse.(Float64,line[9:11]))
    p4 = Point3D(parse.(Float64,line[12:14]))
    # add to model
    push!(
        active_part(model,state).quadrilateral_geometry,
        NgonElement(color,GeometryBasics.Quadrilateral(p1,p2,p3,p4))
        )
    return state
end

"""
    read_optional_line!

For reading lines of type OPTIONAL_LINE
"""
function read_optional_line!(model,state,line)
    @assert parse_command_code(line[1]) == OPTIONAL_LINE
    @assert length(line) == 14
    color = parse_color(line[2])
    p1 = Point3D(parse.(Float64,line[3:5]))
    p2 = Point3D(parse.(Float64,line[6:8]))
    p3 = Point3D(parse.(Float64,line[9:11]))
    p4 = Point3D(parse.(Float64,line[12:14]))
    # add to model
    push!(
        active_part(model,state).optional_line_geometry,
        OptionalLineElement(
            color,
            Line(p1,p2),
            Line(p3,p4)
        ))
    return state
end

"""
    populate_part_geometry!(model,frontier=Set(collect(part_keys(model))))

Load all geometry into `model.parts`. Loading is recursive, so that geometry
will be loaded through arbitrary levels of nested subparts until finally being
stored in each atomic part that is referenced by the main model(s).
"""
function populate_part_geometry!(model,frontier=Set(collect(part_keys(model))))
    explored = Set{String}()
    while !isempty(frontier)
        subcomponent = pop!(frontier)
        push!(explored,subcomponent)
        partfile = find_part_file(subcomponent)
        if partfile === nothing
            @warn "Can't find file $subcomponent to populate geometry. Skipping..."
            continue
        end
        parse_ldraw_file!(model,partfile,MPDModelState(active_part=subcomponent))
        # for k in keys(model.parts)
        for k in all_part_keys(model)
            if !(k in explored)
                push!(frontier,k)
            end
        end
    end
    # go back through and recursively load geometry from subcomponents
    frontier = explored
    explored = Set{String}()
    while !isempty(frontier)
        name = pop!(frontier)
        part = get_part(model,name)
        recurse_part_geometry!(model,part,explored)
        setdiff!(frontier,explored)
    end
    model
end
function recurse_part_geometry!(model,part::DATModel,explored)
    if get_status(part.populated)
        push!(explored,part.name)
        return part
    end
    for ref in part.subfiles
        if has_part(model,ref.file)
            subpart = get_part(model,ref.file)
            if !(ref.file in explored)
                recurse_part_geometry!(model,subpart,explored)
            end
            incorporate_geometry!(part,ref,subpart)
        end
    end
    set_status!(part.populated,true)
    push!(explored,part.name)
    part
end

const LDRAW_BASE_FRAME = SMatrix{3,3,Float64}(
    1.0,  0.0,  0.0,
    0.0,  0.0,  -1.0,
    0.0,  1.0,  0.0
)
"""
    ldraw_base_frame()

Returns a rotation matrix that defines the base LDraw coordinate system.
"""
ldraw_base_frame() = deepcopy(LDRAW_BASE_FRAME)
ldraw_base_transform() = Translation(0.0,0.0,0.0) ∘ LinearMap(ldraw_base_frame())

function (t::AffineMap)(g::G) where {G<:GeometryBasics.Ngon}
    G(map(t,g.points))
end
(t::AffineMap)(g::G) where {G<:NgonElement} = G(g.color,t(g.geom))
(t::AffineMap)(g::G) where {G<:OptionalLineElement} = G(g.color,t(g.geom),t(g.control_pts))
function Base.:(*)(r::Rotation,g::G) where {G<:GeometryBasics.Ngon}
    G(map(p->r*p,g.points))
end
function Base.:(*)(g::G,x::Float64) where {G<:GeometryBasics.Ngon}
    G(map(p->p*x,g.points))
end
Base.:(*)(g::G,x::Float64) where {G<:NgonElement} = G(g.color,g.geom*x)
Base.:(*)(g::G,x::Float64) where {G<:OptionalLineElement} = G(g.color,g.geom*x,g.control_pts*x)
function scale_translation(T::AffineMap,scale::Float64)
    compose(Translation(scale*T.translation...), LinearMap(T.linear))
end

"""
    change_coordinate_system!(model::MPDModel,T)

Transform the coordinate system of the entire model
"""
function change_coordinate_system!(model::MPDModel,T=ldraw_base_transform(),scale=1.0)
    Tinv = inv(T)
    for (k,m) in model.models
        for step in m.steps
            for i in 1:length(step.lines)
                ref = step.lines[i]
                t = build_transform(ref)
                new_t = scale_translation(T ∘ t ∘ Tinv,scale)
                step.lines[i] = SubFileRef(ref,new_t)
            end
        end
    end
    for k in collect(all_part_keys(model))
        part = get_part(model,k)
        new_part = DATModel(part,T,scale)
        set_part!(model,new_part,k)
    end
    return model
end

# For rendering
function GeometryBasics.decompose(
        ::Type{TriangleFace{Int}},
        n::GeometryBasics.Ngon{3,Float64,N,Point{3,Float64}}) where {N}
    return SVector{N-1,TriangleFace{Int}}(
        TriangleFace{Int}(1,i,i+1) for i in 1:N-1
    )
end
function GeometryBasics.decompose(
        ::Type{Point{3,Float64}},
        n::GeometryBasics.Ngon)
    return n.points
end
GeometryBasics.faces(n::GeometryBasics.Ngon{3,Float64,4,Point{3,Float64}}) = [
    TriangleFace{Int}(1,2,3),TriangleFace{Int}(1,3,4)
]
GeometryBasics.faces(n::GeometryBasics.Ngon{3,Float64,3,Point{3,Float64}}) = [
    TriangleFace{Int}(1,2,3)
]
GeometryBasics.coordinates(n::GeometryBasics.Ngon{3,Float64,N,Point{3,Float64}}) where {N} = Vector(n.points)

function GeometryBasics.coordinates(v::AbstractVector{G}) where {N,G<:GeometryBasics.Ngon{3,Float64,N,Point{3,Float64}}}
    vcat(map(coordinates,v)...)
end
function GeometryBasics.faces(v::AbstractVector{G}) where {N,G<:GeometryBasics.Ngon{3,Float64,N,Point{3,Float64}}}
    vcat(map(i->map(f->f .+ ((i-1)*N),faces(v[i])),1:length(v))...)
end
function GeometryBasics.faces(v::AbstractVector{G}) where {G<:GeometryBasics.Ngon}
    face_vec = Vector{TriangleFace{Int}}()
    i = 0
    for element in v
        append!(face_vec, map(f->f.+i,faces(element)))
        i = face_vec[end].data[3]
    end
    return face_vec
end
function GeometryBasics.coordinates(v::AbstractVector{G}) where {G<:GeometryBasics.Ngon}
    coords = Vector{Point{3,Float64}}()
    for element in v
        append!(coords, coordinates(element))
    end
    return coords
end

end


text:
module mutil001

using Test
using LinearAlgebra
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.AssemblyModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!

function _test()
   SM = AssemblyModule
   a = SM.SysmatAssemblerSparseCSRSymm(0.0)
   SM.startassembly!(a, 5, 5, 3, 7, 7)
   m = [0.24406   0.599773    0.833404  0.0420141
    0.786024  0.00206713  0.995379  0.780298
    0.845816  0.198459    0.355149  0.224996]
    # @show    m'*m
   SM.assemble!(a, m'*m, [5 2 1 4], [5 2 1 4])
   m = [0.146618  0.53471   0.614342    0.737833
     0.479719  0.41354   0.00760941  0.836455
     0.254868  0.476189  0.460794    0.00919633
     0.159064  0.261821  0.317078    0.77646
     0.643538  0.429817  0.59788     0.958909]
     # @show    m'*m
   SM.assemble!(a, m'*m, [2 3 1 5], [2 3 1 5])
   A = SM.makematrix!(a)
   # @show A

   a = SysmatAssemblerSparseSymm(0.0)
   startassembly!(a, 5, 5, 3, 7, 7)
   m = [0.24406   0.599773    0.833404  0.0420141
    0.786024  0.00206713  0.995379  0.780298
    0.845816  0.198459    0.355149  0.224996]
    # @show    m'*m
   assemble!(a, m'*m, [5 2 1 4], [5 2 1 4])
   m = [0.146618  0.53471   0.614342    0.737833
     0.479719  0.41354   0.00760941  0.836455
     0.254868  0.476189  0.460794    0.00919633
     0.159064  0.261821  0.317078    0.77646
     0.643538  0.429817  0.59788     0.958909]
     # @show    m'*m
   assemble!(a, m'*m, [2 3 1 5], [2 3 1 5])
   A1 = makematrix!(a)
   @test norm(A - A1) / norm(A1) < 1.0e-9
   return true
end

_test()

end # module


text:
using DashBootstrapComponents

carousel = html_div([
    dbc_carousel(
        id = "carousel",
        items = [
            Dict("key" => "1", "src" => "/static/images/slide1.svg"),
            Dict("key" => "2", "src" => "/static/images/slide2.svg"),
            Dict("key" => "3", "src" => "/static/images/slide3.svg"),
        ],
        controls = false,
        indicators = false,
        interval = nothing,
    ),
    dbc_radioitems(
        id = "slide-number",
        options = [
            Dict("label" => "Slide 1", "value" => 0),
            Dict("label" => "Slide 2", "value" => 1),
            Dict("label" => "Slide 3", "value" => 2),
        ],
        value = 0,
        inline = true,
    ),
]);

callback!(app, Output("carousel", "active_index"), Input("slide-number", "value")) do idx
    return idx
end;

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components