text:
# This file covers everything related to recruitment.

# The functions of the Recruitment type require SimJulia and ResumableFunctions.

# The functions of the Retirement type require the Personnel, PersonnelDatabase,
#   ManpowerSimulation, Attrition, and Retirement types.
requiredTypes = [ "manpowerSimulation",
    "attrition", "retirement" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

# This function sets the recruitment schedule.
export setRecruitmentSchedule
function setRecruitmentSchedule( recScheme::Recruitment, freq::T1,
    offset::T2 = 0.0 ) where T1 <: Real where T2 <: Real

    if freq <= 0.0
        error( "Recruitment frequency must be > 0.0." )
    end  # if freq <= 0.0

    recScheme.recruitFreq = freq
    recScheme.recruitOffset = offset % freq + ( offset < 0.0 ? freq : 0.0 )

end  # setRecruitmentSchedule( recScheme, freq, offset )

export setRecruitmentLimits
"""
```
setRecruitmentLimits( recScheme::Recruitment,
                      minRec::T1,
                      maxRec::T2 )
    where T1 <: Real
    where T2 <: Real
```
This function sets the minimum and maximum number of people to recruit using
recruitment scheme `recScheme` to `minRec` and `maxRec` respectively. If the
values don't make sense (maxRec <= 0 or minRec > maxRec), no change is made.

This function returns `nothing`. If the arguments don't make sense, the function
issues a warning.
"""
function setRecruitmentLimits( recScheme::Recruitment, minRec::T1,
    maxRec::T2 )::Void where T1 <: Real where T2 <: Real

    if maxRec < 0
        warn( "Max recruitment must be > 0. Not making any changes." )
        return
    end  # if maxRec <= 0

    if minRec > maxRec
        warn( "Min recruitment must be <= max recruitment. Not making any changes." )
        return
    end  # if minRec > maxRec

    recScheme.isAdaptive = true
    recScheme.minRecruit = max( 0, minRec )
    recScheme.maxRecruit = maxRec

    return

end  # setRecruitmentLimits( recScheme, minRec, maxRec )

# This function sets the state to recruit into. Note that it doesn't check if
#   the state is defined in the simulation. If it isn't, the function will
#   perform general recruitment instead.
export setRecruitState
function setRecruitState( recScheme::Recruitment, stateName::String )::Void

    recScheme.recState = stateName
    return

end  # setRecruitState( recScheme, stateName )

# This function sets the fixed number of people to recruit during every
#   recruitment cycle.
export setRecruitmentFixed
function setRecruitmentFixed( recScheme::Recruitment, amount::T ) where T <: Integer

    if amount < 0
        warn( "Negative number of people to recruit entered. Not making changes to recruitment scheme." )
        return
    end  # if amount <= 0

    recScheme.isAdaptive = false
    recScheme.recDistType = :disc
    recScheme.recDistNodes = Dict( Float64( amount ) => 1.0 )
    recScheme.recDist = function() return amount end

end  # setRecruitmentFixed( recScheme, amount )

# This function sets the recruitment distribution of the recruitment scheme to
#   the given distribution with the given nodes.
export setRecruitmentDistribution
function setRecruitmentDistribution( recScheme::Recruitment,
    recDist::Dict{Int, Float64}, recDistType::Symbol )

    # Check if the distribution type is known.
    if recDistType ∉ [ :disc, :pUnif, :pLin ]
        warn( "Unknown distribution type. Recruitment distribution not set." )
        return
    end  # if recDistType ∉ [ :disc, :pUnif, :pLin ]

    # Get the list of proper nodes (recruitment >= 0 and p >= 0).
    tmpNodes = collect( keys( recDist ) )
    tmpNodes = tmpNodes[ map( node -> ( node >= 0.0 ) &&
        ( recDist[ node ] >= 0.0 ), tmpNodes ) ]
    sort!( tmpNodes )

    # Check if there are sufficient nodes.
    if length( tmpNodes ) < ( recDistType === :disc ? 1 : 2 )
        warn( "Not enough distribution nodes. Recruitment distribution not set." )
        return
    end  # if length( tmpNodes ) < ( recDistType === :disc ? 1 : 2 )

    # Check if the total probability mass is non-zero.
    pMass = 0.0
    foreach( node -> pMass += recDist[ node ], tmpNodes )

    if ( pMass == 0.0 ) || ( ( recDistType === :pUnif ) &&
        ( recDist[ tmpNodes[ end ] ] == pMass ) )
        warn( "Proposed distribution has 0 probability mass. Recruitment distribution not set." )
        return
    end  # if pMass == 0.0

    # Set the distribution
    recScheme.isAdaptive = false
    recScheme.recDistType = recDistType
    recScheme.recDistNodes = recDist
    distFuncs = Dict{Symbol, Function}(
        :disc => setDiscRecDist,
        :pUnif => setPunifRecDist,
        :pLin => setPlinRecDist )
    ( distFuncs[ recDistType ] )( recScheme, recDist, tmpNodes )

end  # setRecruitmentDistribution( recScheme, recDist, recDistType )

# This function sets the recruitment distribution of the recruitment scheme to a
#   discrete distribution with the given nodes.
function setDiscRecDist( recScheme::Recruitment, recDist::Dict{Int, Float64},
    nodes::Vector{Int} )

    # Get the point probabilities of the nodes.
    pNodes = map( node -> recDist[ node ], nodes )
    pNodes /= sum( pNodes )

    recScheme.recDist = function()
        return nodes[ rand( Categorical( pNodes ) ) ]
    end

end  # setDiscAgeDist( recScheme, ageDist, nodes )

# This function sets the recruitment distribution of the recruitment scheme to a
#   piecewise uniform distribution with the given nodes.
function setPunifRecDist( recScheme::Recruitment, recDist::Dict{Int, Float64},
    nodes::Vector{Int} )

    # Get the point probabilities of the intervals.
    pInts = map( node -> recDist[ node ], nodes[1:(end-1)] )
    pInts /= sum( pInts )

    recScheme.recDist = function()
        intI = rand( Categorical( pInts ) )
        return rand( DiscreteUniform( nodes[ intI ], nodes[ intI + 1 ] - 1 ) )
    end

end  # setPunifRecDist( recScheme, recDist, nodes )

# This function sets the recruitment distribution of the recruitment scheme to a
#   piecewise linear distribution with the given nodes.
function setPlinRecDist( recScheme::Recruitment, recDist::Dict{Int, Float64},
    nodes::Vector{Int} )

    # Get the point probabilities of the intervals.
    pointWeights = map( node -> recDist[ node ], nodes )
    nodeDiff = nodes[ 2:end ] - nodes[ 1:(end-1) ]
    weightDiff = pointWeights[ 2:end ] - pointWeights[ 1:(end-1) ]
    bracketWeights = map(
        ii -> 0.5 * ( pointWeights[ ii + 1 ] * ( nodeDiff[ ii ] + 1.0 ) +
        pointWeights[ ii ] * ( nodeDiff[ ii ] - 1.0 ) ),
        1:length( weightDiff ) )
    cumulWeights = cumsum( vcat( pointWeights[ 1 ], bracketWeights ) )

    recScheme.recDist = function()
        weight = rand( Uniform( 0, cumulWeights[ end ] ) )
        intI = findlast( weight .>= cumulWeights )

        if intI == 0
            return nodes[ 1 ]
        end  # if intI == 0

        wTilde = weight - cumulWeights[ intI ]
        weightNodeRatio = weightDiff[ intI ] / nodeDiff[ intI ]
        polyEq = Poly( [ -wTilde, pointWeights[ intI ] + weightNodeRatio / 2,
            weightNodeRatio / 2 ] )
        pos = filter( root -> 0 <= root < nodeDiff[ intI ],
            roots( polyEq ) )[ 1 ]

        return nodes[ intI ] + ceil( Int, pos )
    end

end  # setPunifRecDist( recScheme, recDist, nodes )

# This function sets the recruitment age to a single, fixed number.
export setRecruitmentAge
function setRecruitmentAge( recScheme::Recruitment, age::T ) where T <: Real

    if age < 0.0
        warn( "Age at recruitment must be ⩾ 0.0. Not making changes to recruitment scheme." )
        return
    end

    recScheme.ageDistType = :disc
    recScheme.ageDistNodes = Dict( Float64( age ) => 1.0 )
    recScheme.ageDist = function() return age end

end  # setRecruitmentAge( recScheme, age )

# This function sets the age distribution of the recruitment scheme to the
#   given distribution with the given nodes.
export setAgeDistribution
function setAgeDistribution( recScheme::Recruitment,
    ageDist::Dict{Float64, Float64}, ageDistType::Symbol )

    # Check if the distribution type is known.
    if ageDistType ∉ [ :disc, :pUnif, :pLin ]
        warn( "Unknown distribution type. Age distribution not set." )
        return
    end  # if recDistType ∉ [ :disc, :pUnif, :pLin ]

    # Get the list of proper nodes (age >= 0 and p >= 0).
    tmpNodes = collect( keys( ageDist ) )
    tmpNodes = tmpNodes[ map( node -> ( node >= 0.0 ) &&
        ( ageDist[ node ] >= 0.0 ), tmpNodes ) ]
    sort!( tmpNodes )

    # Check if there are sufficient nodes.
    if length( tmpNodes ) < ( ageDistType === :disc ? 1 : 2 )
        warn( "Not enough distribution nodes. Age distribution not set." )
        return
    end  # if length( tmpNodes ) < ( recDistType === :disc ? 1 : 2 )

    # Check if the total probability mass is non-zero.
    pMass = 0.0
    foreach( node -> pMass += ageDist[ node ], tmpNodes )

    if ( pMass == 0.0 ) || ( ( ageDistType === :pUnif ) &&
        ( ageDist[ tmpNodes[ end ] ] == pMass ) )
        warn( "Proposed distribution has 0 probability mass. Age distribution not set." )
        return
    end  # if pMass == 0.0

    # Set the distribution
    recScheme.ageDistType = ageDistType
    recScheme.ageDistNodes = ageDist
    distFuncs = Dict{Symbol, Function}(
        :disc => setDiscAgeDist,
        :pUnif => setPUnifAgeDist,
        :pLin => setPLinAgeDist )
    ( distFuncs[ ageDistType ] )( recScheme, ageDist, tmpNodes )

end  # setAgeDistribution( recScheme, ageDist, recDistType )

# This function sets the age distribution of the recruitment scheme to a
#   discrete distribution with the given nodes.
function setDiscAgeDist( recScheme::Recruitment,
    ageDist::Dict{Float64, Float64}, nodes::Vector{Float64} )

    # Get the point probabilities of the nodes.
    pNodes = map( node -> ageDist[ node ], nodes )
    pNodes /= sum( pNodes )

    recScheme.ageDist = function()
        return nodes[ rand( Categorical( pNodes ) ) ]
    end

end  # setDiscAgeDist( recScheme, ageDist, nodes )

# This function sets the age distribution of the recruitment scheme to a
#   piecewise uniform distribution with the given nodes.
# The probability at node t_i is the (unweighted) probability that the result
#   lies in the interval [ t_i, t_{i+1} ]. Hence, the probabilty at the last
#   node is ignored.
function setPUnifAgeDist( recScheme::Recruitment,
    ageDist::Dict{Float64, Float64}, nodes::Vector{Float64} )

    # Get the point probabilities of the intervals.
    pInts = map( node -> ageDist[ node ], nodes[1:(end-1)] )
    pInts /= sum( pInts )

    recScheme.ageDist = function()
        intI = rand( Categorical( pInts ) )
        return rand( Uniform( nodes[ intI ], nodes[ intI + 1 ] ) )
    end

end  # setPUnifAgeDist( recScheme, ageDist, nodes )

# This function sets the age distribution of the recruitment scheme to a
#   piecewise linear distribution with the given nodes.
function setPLinAgeDist( recScheme::Recruitment,
    ageDist::Dict{Float64, Float64}, nodes::Vector{Float64} )

    # Get the point probabilities of the intervals.
    pNodes = map( node -> ageDist[ node ], nodes )
    pMasses = map( ii -> ( nodes[ ii + 1 ] - nodes[ ii ] ) *
        ( pNodes[ ii ] + pNodes[ ii + 1 ] ) / 2.0, eachindex( nodes[ 2:end ] ) )
    pMasses /= sum( pMasses )

    recScheme.ageDist = function()
        intI = rand( Categorical( pMasses ) )
        pDiff = pNodes[ intI + 1 ] - pNodes[ intI ]

        if pDiff == 0
            dist = Uniform( nodes[ intI ], nodes[ intI + 1 ] )
        else
            mu = pDiff > 0 ? nodes[ intI + 1 ] : nodes[ intI ]
            sigma = max( pNodes[ intI ], pNodes[ intI + 1 ] ) *
                ( nodes[ intI + 1 ] - nodes[ intI ] ) /
                abs( pNodes[ intI + 1 ] - pNodes[ intI ] )
            dist = Truncated( SymTriangularDist( mu, sigma ), nodes[ intI ],
                nodes[ intI + 1 ] )
        end  # if pDiff == 0

        return rand( dist )
    end

end  # setPLinAgeDist( recScheme, ageDist, nodes )

# This function generates and returns the number of personnel members to recruit
#   in a recruitment cycle.
function generatePoolSize( mpSim::ManpowerSimulation, recScheme::Recruitment )

    nrToRecruit = recScheme.maxRecruit

    if recScheme.isAdaptive && ( mpSim.personnelTarget > 0 )
        personnelNeeded = mpSim.personnelTarget - mpSim.personnelSize
        nrToRecruit = max( min( personnelNeeded, nrToRecruit ),
            recScheme.minRecruit )
    elseif !recScheme.isAdaptive
        nrToRecruit = recScheme.recDist()
    end  # if ( mpSim.personnelTarget > 0 ) && ...

    return nrToRecruit

end  # generatePoolSize( mpSim, recScheme )

function generatePoolSize( mpSim::ManpowerSimulation, recScheme::Recruitment,
    recState::State )

    nrToRecruit = recScheme.maxRecruit

    if recScheme.isAdaptive
        persToOrgTarget = mpSim.personnelTarget > 0 ? mpSim.personnelTarget -
            mpSim.personnelSize : typemax( Int )
        persToStateTarget = recState.stateTarget >= 0 ? recState.stateTarget -
            length( recState.inStateSince ) : typemax( Int )
        nrToRecruit = max( min( persToOrgTarget, persToStateTarget,
            nrToRecruit ), recScheme.minRecruit )
    elseif !recScheme.isAdaptive
        nrToRecruit = recScheme.recDist()
    end  # if ( mpSim.personnelTarget > 0 ) && ...

    return nrToRecruit

end  # generatePoolSize( mpSim, recScheme )

# This function generates a single personnel member using the information in
#   the recruitment scheme.
function createPerson( mpSim::ManpowerSimulation, recScheme::Recruitment,
    recState::T = nothing ) where T <: Union{Void, State}

    # Create the person in the database.
    id = "Sim" * string( mpSim.resultSize + 1 )
    ageAtRecruitment = recScheme.ageDist()

    # Generate initial values.
    initVals = Dict{String, Any}()
    initValsFixed = Dict{String, Bool}()

    for ii in eachindex( mpSim.initAttrList )
        attr = mpSim.initAttrList[ ii ]
        initVals[ attr.name ] = generateAttrValue( attr )
        initValsFixed[ attr.name ] = attr.isFixed
    end  # for attr in mpSim.initAttrList

    initPersStates = []

    if isa( recState, State )
        # When recruiting for a specific state, set the required attributes, and
        #   put entity in that state.
        for attr in keys( recState.requirements )
            initVals[ attr ] = recState.requirements[ attr ][ 1 ]
            initValsFixed[ attr ] = false
        end  # for attr in keys( recState.requirements )

        initPersStates = [ recState ]
    else
        # Otherwise, identify all initial states the person belongs to and add
        #   entry info to each of those states.
        initPersStates = collect( Iterators.filter(
            state -> isPersonnelOfState( initVals, state ),
            keys( mpSim.initStateList ) ) )
            # XXX Iterators.filter is needed to avoid deprecation warnings.
    end  # if isa( recState, State )

    # Check if each person can be assigned to exactly one state.
    if mpSim.isWellDefined
        if isempty( initPersStates )
            mpSim.isWellDefined = false
            warn( "Person created that couldn't be assigned to any state. Please check system configuration for consistency." )
        elseif length( initPersStates ) > 1
            mpSim.isWellDefined = false
            warn( "Person created that can be assigned to multiple states. Please check system configuration for consistency." )
        end  # if isempty( initPersStates )
    end  # if mpSim.isWellDefined

    for state in initPersStates
        state.inStateSince[ id ] = now( mpSim )
        state.isLockedForTransition[ id ] = false
    end  # for state in initPersStates

    stateRetAge = isempty( initPersStates ) ? 0 :
        initPersStates[ 1 ].stateRetAge

    stateNames = map( state -> state.name, initPersStates )
    timeOfRetirement = computeExpectedRetirementTime( mpSim,
        mpSim.retirementScheme, ageAtRecruitment, stateRetAge, now( mpSim ) )
    attrScheme = determineAttritionScheme( initPersStates, mpSim )
    timeOfAttr = generateTimeOfAttrition( attrScheme, now( mpSim ) )

    # Add person to the personnel database.
    command = "INSERT INTO $(mpSim.personnelDBname)
        ($(mpSim.idKey), status, timeEntered, ageAtRecruitment,
            expectedRetirementTime, expectedAttritionTime, attritionScheme"

    if !isempty( initVals )
        command *= ", '$(join( keys( initVals ), "', '" ))'"
    end  # if !isempty( initVals )

    command *= ") VALUES
        ('$id', 'active', $(now( mpSim )), $ageAtRecruitment, $timeOfRetirement,
            $timeOfAttr, '$(attrScheme.name)'"

    if !isempty( initVals )
        command *= ", '$(join( map( attrName -> initVals[ attrName ], keys( initVals ) ), "', '" ))'"
    end # if !isempty( initVals )

    command *= ")"
    SQLite.execute!( mpSim.simDB, command )

    # Add entry of person to the history database.
    command = "INSERT INTO $(mpSim.historyDBname)
        ($(mpSim.idKey), attribute, timeIndex, strValue) VALUES
        ('$id', 'status', $(now( mpSim )), 'active')"

    # Additional variable attributes.
    for attrName in keys( initVals )
        if !initValsFixed[ attrName ]
            command *= ", ('$id', '$attrName', $(now( mpSim )), '$(initVals[ attrName ])')"
        end  # !initValsFixed[ attrName ]
    end  # for attrName in keys( initVals )

    SQLite.execute!( mpSim.simDB, command )

    # Add recruitment event and all applicable initial states to transition
    #   database.
    command = "INSERT INTO $(mpSim.transitionDBname)
        ($(mpSim.idKey), timeIndex, transition, endState) VALUES
        ('$id', $(now( mpSim )), '$(recScheme.name)', 'active')"

    # Initiate all applicable transition processes.
    for state in initPersStates
        command *= ", ('$id', $(now( mpSim )), '$(recScheme.name)', '$(state.name)')"

        # for trans in mpSim.initStateList[ state ]
        #     @process transitionProcess( mpSim.sim, trans, id, mpSim )
        # end  # for trans in mpSim.initStateList[ state ]
    end  # for state in initPersStates

    SQLite.execute!( mpSim.simDB, command )

    # Adjust the size of the personnel database.
    mpSim.personnelSize += 1
    mpSim.resultSize += 1

end  # createPerson( mpSim, recScheme )

# This function performs the recruitment part of a single recruitment cycle.
function recruitmentCycle( mpSim::ManpowerSimulation, recScheme::Recruitment )

    nrToRecruit = generatePoolSize( mpSim, recScheme )
    recState = recScheme.recState == "" ? nothing :
        mpSim.stateList[ recScheme.recState ]

    if !isa( recState, Void )
        nrToRecruit = generatePoolSize( mpSim, recScheme, recState )
    end  # if !isa( recState, Void )

    # Stop here if no recruitments happen in this period.
    if nrToRecruit == 0
        return
    end  # if nrToRecruit == 0

    foreach( ii -> createPerson( mpSim, recScheme, recState ), 1:nrToRecruit )

    return

end  # recruitmentCycle( mpSim, recScheme )

# This is the process in the simulation for a single recruitment scheme.
# Version for JuliaBox
@resumable function recruitProcess( sim::Simulation,
    schemeNr::Integer, mpSim::ManpowerSimulation )

    processTime = Dates.Millisecond( 0 )
    tStart = now()

    recScheme = mpSim.recruitmentSchemes[ schemeNr ]
    timeToWait = recScheme.recruitOffset
    priority = mpSim.phasePriorities[ :recruitment ]
    priority += recScheme.isAdaptive ? 0 : 1

    while now( sim ) + timeToWait <= mpSim.simLength
        processTime += now() - tStart
        @yield timeout( sim, timeToWait, priority = priority )
        tStart = now()
        timeToWait = recScheme.recruitFreq
        recruitmentCycle( mpSim, recScheme )
    end  # while now( sim ) + timeToWait <= mpSim.simLength

    processTime += now() - tStart
    println( "Recruitment process for '$(recScheme.name)' took $(processTime.value / 1000) seconds." )

end

function Base.show( io::IO, recScheme::Recruitment )

    print( io, "Recruitment schedule: $(recScheme.recruitFreq) (+ $(recScheme.recruitOffset))\n" )
    print( io, recScheme.isAdaptive ? "A" : "Non-a" )
    print( io, "daptive recruitment scheme" )

    if recScheme.isAdaptive
        print( io, "\nRecruitment per cycle: $(recScheme.minRecruit) - $(recScheme.maxRecruit)" )
    end  # if recScheme.isAdaptive || ...

end  # show( io, recScheme )

# ==============================================================================
# Non-exported methods
# ==============================================================================

function readRecruitmentScheme( sheet::XLSX.Worksheet,
    ii::T )::Recruitment where T <: Integer

    dataColNr = ii * 5 - 3
    name = sheet[ XLSX.CellRef( 5, dataColNr ) ]
    recScheme = Recruitment( name, sheet[ XLSX.CellRef( 6, dataColNr ) ],
        sheet[ XLSX.CellRef( 7, dataColNr ) ] )
    recState = sheet[ XLSX.CellRef( 8, dataColNr ) ]
    recState = isa( recState, Missings.Missing ) ||
        lowercase( recState ) == "active" ? "" : recState
    setRecruitState( recScheme, recState )
    isAdaptive = sheet[ XLSX.CellRef( 11, dataColNr ) ] == "YES"
    isRandom = sheet[ XLSX.CellRef( 12, dataColNr ) ] == "YES"
    nRow = 16
    numNodes = sheet[ XLSX.CellRef( nRow + 1, dataColNr ) ]

    if isAdaptive
        minRec = sheet[ XLSX.CellRef( 9, dataColNr ) ] === nothing ? 0 :
            sheet[ XLSX.CellRef( 9, dataColNr ) ]
        setRecruitmentLimits( recScheme, minRec,
            sheet[ XLSX.CellRef( 10, dataColNr ) ] )
    elseif isRandom
        distType = distTypes[ sheet[ XLSX.CellRef( nRow, dataColNr ) ] ]
        minNodes = distType == "Pointwise" ? 1 : 2
        recDist = Dict{Int, Float64}()

        for jj in (1:numNodes) + 2
            node = sheet[ XLSX.CellRef( nRow + jj, dataColNr ) ]
            weight = sheet[ XLSX.CellRef( nRow + jj, dataColNr + 1 ) ]

            if isa( node, Real ) && !haskey( recDist, node ) && ( node >= 0 ) &&
                ( weight >= 0 )
                recDist[ floor( Int, node ) ] = weight
            end  # if isa( node, Real ) && ...
        end  # for ii in (1:numNodes) + 2

        if length( recDist ) < minNodes
            error( "Recruitment type $name has an insufficient number of valid nodes defined for its population size distribution." )
        end  # if numNodes < minNodes

        setRecruitmentDistribution( recScheme, recDist, distType )
    else
        setRecruitmentFixed( recScheme, sheet[ XLSX.CellRef( 10, dataColNr ) ] )
    end  # if isAdaptive

    isFixedAge = sheet[ XLSX.CellRef( 13, dataColNr ) ] == "YES"

    # Add the age distribution.
    if isFixedAge
        setRecruitmentAge( recScheme,
            sheet[ XLSX.CellRef( 14, dataColNr ) ] * 12.0 )
    else
        # Get to the start of the age distribution
        nRow += numNodes + 4
        distType = distTypes[ sheet[ XLSX.CellRef( nRow, dataColNr ) ] ]
        numNodes = sheet[ XLSX.CellRef( nRow + 1, dataColNr ) ]
        minNodes = distType == "Pointwise" ? 1 : 2
        ageDist = Dict{Float64, Float64}()

        for jj in (1:numNodes) + 2
            age = sheet[ XLSX.CellRef( nRow + jj, dataColNr ) ]
            pMass = sheet[ XLSX.CellRef( nRow + jj, dataColNr + 1 ) ]

            # Only add the entry if it makes sense.
            if isa( age, Real ) && !haskey( ageDist, age ) && ( age >= 0 ) &&
                ( pMass >= 0 )
                ageDist[ age * 12.0 ] = pMass
            end  # if isa( age, Real ) && ...
        end  # for jj in (1:numNodes) + 2

        if length( ageDist ) < minNodes
            error( "Recruitment type $name has an insufficient number of valid nodes defined for its recruitment age distribution." )
        end  # if numNodes < ( distType == "Pointwise" ? 1 : 2 )

        setAgeDistribution( recScheme, ageDist, distType )
    end  # if isFixedAge

    return recScheme

end  # function generateRecruitmentScheme( s, ii )


text:
using BenchmarkTools, Filtering

SUITE = BenchmarkGroup()

SUITE["ar_linear"] = BenchmarkGroup(["models","linear",])

β = [0.0;1.0]
δ = 1.0
φ = 0.95
ση = 1.0
σε = 1.0

θ0 = [δ;atanh(φ);ση;σε]

β = [0.0;1.0]

N = 1000

arh = ARH(β)

SUITE["ar_linear"]["simulate"] = @benchmarkable Filtering.simulate(arh,θ0,1000)

X,Y = Filtering.simulate(arh,θ0,1000)

SUITE["ar_linear"]["kalman"] = @benchmarkable kalman_filter(arh,θ0,Y,fill([1.0],length(Y)))
SUITE["ar_linear"]["particle"] = @benchmarkable particle_filter(arh,θ0,Y,N,fill([1.0],length(Y)))

ext: jl


text:
module DoctorDocstrings

using TerminalMenus
using Printf
using TypedTables
using REPL
import InteractiveUtils: clipboard, @edit

export list_no_docs, diagnosedocs, fixdocs, pickandcopy

getdocs(f) = Docs.Text(Docs.doc(f))
hasdocs(f) = !occursin(r"^No documentation found.", string(getdocs(f)))

function hasexamples(f)
    s = string(getdocs(f))
    return occursin("Examples", s) && occursin("jldoctest", s) && occursin("julia>", s)
end

function make_data(mod)
    length(names(mod)) < 1 && Error("Your $mod does not export any names")
    lennames = length(names(mod))
    col1 = [getfield(mod, i) for i in names(mod) if isdefined(mod, i)]
    col2 = [hasdocs(j) for j in col1]
    col3 = [hasexamples(k) for k in col1]

    f1 = "Total %"
    f2 = @sprintf("%.2f", count(hasdocs(j) for j in col1)/lennames)
    f3 = @sprintf("%.2f", count(hasexamples(k) for k in col1)/lennames)
    footer = [f1 f2 f3]

    data = vcat(hcat(col1, col2, col3), footer)
    return data
end

#function format_data(data, fun = false)
#    if !fun # Boo! :(
#        nope = '-'
#        yup = "Yes"
#    else    # Yay! :)
#        nope = '💩'
#        yup = '🎉'
#    end
#    f_data = copy(data)
#    len = size(data)[2]
#    # God this is ugly
#    for i in len-1:len
#        f_data[i,2] = f_data[i, 2] ? yup : nope
#        f_data[i,3] = f_data[i, 3] ? yup : nope
#    end
#    f_data
#end

function list_no_docs(mod, footer = true, fun = false)
    data = make_data(mod)
#    f_data = format_data(data, fun)
    table = Table(Functions = data[:,1], Docs = data[:,2], Examples = data[:,3])
#    return f_data, table
    return footer ? table : table[1:end-1, :]
end

# TODO Fix buggy shell and pkg starting lines
function history_parser(n = 50, maxlen = 1000)
    count = 0
    lines = reverse(readlines(REPL.find_hist_file()))
    len = length(lines)
    res = String[]
    i = 2
    while count < n && i < min(len, maxlen)
        temp = ""
        start = ""
        if endswith(lines[i], "julia")
            start = "julia> "
        elseif endswith(lines[i], "shell")
                start = ";"
        elseif endswith(lines[i], "pkg")
                start = "]"
        end
        # skip over the # time
        i += 2
        while !startswith(lines[i], "# mode: ") && i < min(len, maxlen)
            temp = lines[i][2:end] * '\n' * temp
            i += 1
        end
        push!(res, start * temp)
        count += 1
    end
    reverse!(res)
end

function pickandcopy(n=25, keeprepl = false)
    options = history_parser(n)
    menu = MultiSelectMenu(options, pagesize = n)
    choiceset = request("Pick your history to copy:", menu)
    picks = [options[i] for i in sort(collect(choiceset))]
    latest_entry = picks[end]
    res = join(picks, '\n')
    if length(choiceset) > 0
        @info "Copied $(length(choiceset)) items to clipboard"
    else
        @info "Something got borked. Try again."
    end
    clipboard(res)
    # TODO Fix this horrible hack to not get the "julia> "
    # latest_entry is a string, so this should "work"
    res, latest_entry[7:end]
end

findfixables(table) = table.Functions[.!(table.Docs .| table.Examples)]

function buildpastestring(template_str, copypicks)
    template_str * """
Examples
====
```jldoctest
""" * '\n' * copypicks * "\n```\n\"\"\""
end

#1. Get REPL history DONE
#2. find first fixable doc string
#2.5 Choose a template
#3. append Example header * REPL History to template,
#3.5 put that into the clipboard
#4. open editor
function latesteditprompt(latest_entry)
    @info "Templated docstrings with example has been pasted into your system clipboard"
    @info "Would you like to @edit $(latest_entry[1:end-1]) ? Y/N"
    if any(occursin.(["yes", "Y", "YES", "y", "Yes"], readline()))
        expr = Meta.parse(latest_entry)
        @eval @edit $(expr)
    else
        @info "Didn't get a Yes/YES/y/Y/Yes answer. Appointment canceled."
    end
end

function instructprompt(doctarget = "your functions")
    @info "You should now write a few examples of a docstring for $doctarget in the REPL"
    @info "You will be able to copy/paste them by calling 'fixdocs()'"
end

function diagnosedocs(mod, verbose = true)
    table = list_no_docs(mod, false, false)
    fixables = findfixables(table)
    if isempty(fixables)
        return println("Yay! All your exports have docs and examples!")
    end
    if verbose
        instructprompt()
    end
    return Table(NeedFixes = fixables)
end

function fixdocs(template = picktemplate(docstr_templates), nhistory = 25)
    picks, latest_entry = pickandcopy(nhistory)
    pastestring = buildpastestring(template, picks)
    clipboard(pastestring)
    latesteditprompt(latest_entry)
end

function picktemplate(docstr_templates)
    menu = RadioMenu(["Basic", "BlueStyle", "Type", "Package quickstart"])
    choice = request("Choose your template", menu)
    if choice != -1
        nothing
    else
        println("Menu canceled")
    end
    docstr_templates[choice]
end

doc_example_template = """
\"\"\"
    foo(T) -> S

This function calls `foo`.

"""
BlueStyle_func_template = """
    mysearch(array::MyArray{T}, val::T; verbose=true) where {T} -> Int

Searches the `array` for the `val`. For some reason we don't want to use Julia's
builtin search :)

# Arguments
- `array::MyArray{T}`: the array to search
- `val::T`: the value to search for
#
# Keywords
- `verbose::Bool=true`: print out progress details

# Returns
- `Int`: the index where `val` is located in the `array`

# Throws
- `NotFoundError`: I guess we could throw an error if `val` isn't found.

Examples
====
 """

type_template = """
     MyArray{T,N}

My super awesome array wrapper!

# Fields
- `data::AbstractArray{T,N}`: stores the array being wrapped
- `metadata::Dict`: stores metadata about the array
"""

quickstart_template = """

Welcome to MY_PACKAGE.jl!

This packages is useful for a) ... b) ... c) ...

Here's the top 5 workflows in this package

Examples
====
"""

docstr_templates = [doc_example_template, BlueStyle_func_template, type_template, quickstart_template]

function docstring_wizard(mod, quiet = true)
    if quiet
        println(banner)
    end
    t = list_no_docs(mod)

    todo_docs = t.Functions[findall(t.Docs)]
    todo_examples = t.Functions[findall(t.Examples)]
    todos = string.(Set([todo_docs, todo_examples]))
    docstr_templates = ["none", "example", "verbose"]
    menu = RadioMenu(docstr_templates)
    choice = request("Choose your template format", menu)

    todos_menu = RadioMenu(todos)
    choice = request("Choose a function to work on", todos_menu)

    if choice != -1
        println("")
    end
end

banner = raw"""
   ..____.......____................._........_...................._._.
   .|.._.\._.__|.._.\..___...___.___|.|_._.__(_)_.__...__._.___...(_).|
   .|.|.|.|.'__|.|.|.|/._.\./.__/.__|.__|.'__|.|.'_.\./._ ./.__|..|.|.|
   .|.|_|.|.|_.|.|_|.|.(_).|.(__\__.\.|_|.|..|.|.|.|.|.(_|.\__.\_.|.|.|
   .|____/|_(_)|____/.\___/.\___|___/\__|_|..|_|_|.|_|\__,.|___(_)/.|_|
   ...................................................|___/.....|__/...
   For the Julia community with <3
""";

end


text:
using LinearAlgebra: Diagonal, Hermitian, ishermitian
using QuantumLattices: contentnames, kind, dimension, azimuth, rcoord, update!, matrix
using QuantumLattices: PID, CPID, Point, Lattice, FID, Fock, NID, Phonon, Index, Hilbert, Metric, OIDToTuple, Parameters
using QuantumLattices: Hopping, Onsite, Pairing, PhononKinetic, PhononPotential
using QuantumLattices: ReciprocalPath, @rectangle_str
using TightBindingApproximation

@time @testset "prerequisite" begin
    @test promote_type(TBAKind{:TBA}, TBAKind{:TBA}) == TBAKind{:TBA}
    @test promote_type(TBAKind{:BdG}, TBAKind{:BdG}) == TBAKind{:BdG}
    @test promote_type(TBAKind{:TBA}, TBAKind{:BdG}) == promote_type(TBAKind{:BdG}, TBAKind{:TBA}) == TBAKind{:BdG}

    @test TBAKind(Hopping) == TBAKind(Onsite) == TBAKind(:TBA)
    @test TBAKind(Pairing) == TBAKind(PhononKinetic) == TBAKind(PhononPotential) == TBAKind(:BdG)
    @test TBAKind(Tuple{Hopping, Onsite}) == TBAKind(:TBA)
    @test TBAKind(Tuple{Hopping, Onsite, Pairing}) == TBAKind(:BdG)

    @test Metric(TBAKind(:TBA), Hilbert(PID(1)=>Fock{:b}(1, 1, 1))) == OIDToTuple(:site, :orbital, :spin)
    @test Metric(TBAKind(:TBA), Hilbert(CPID(1, 1)=>Fock{:f}(1, 1, 1))) == OIDToTuple(:scope, :site, :orbital, :spin)
    @test Metric(TBAKind(:BdG), Hilbert(PID(1)=>Fock{:b}(1, 1, 1))) == OIDToTuple(:nambu, :site, :orbital, :spin)
    @test Metric(TBAKind(:BdG), Hilbert(CPID(1, 1)=>Fock{:f}(1, 1, 1))) == OIDToTuple(:nambu, :scope, :site, :orbital, :spin)
    @test Metric(TBAKind(:BdG), Hilbert(PID(1)=>Phonon(2))) == OIDToTuple(:tag, :site, :dir)
    @test Metric(TBAKind(:BdG), Hilbert(CPID(1, 1)=>Phonon(2))) == OIDToTuple(:tag, :scope, :site, :dir)

    @test commutator(TBAKind(:TBA), Hilbert(PID(1)=>Fock{:f}(1, 2, 2))) == nothing
    @test commutator(TBAKind(:TBA), Hilbert(PID(1)=>Fock{:b}(1, 2, 2))) == nothing
    @test commutator(TBAKind(:BdG), Hilbert(PID(1)=>Fock{:f}(1, 2, 2))) == nothing
    @test commutator(TBAKind(:BdG), Hilbert(PID(1)=>Fock{:b}(1, 2, 2))) == Diagonal([1, 1, -1, -1])
    @test commutator(TBAKind(:BdG), Hilbert(PID(1)=>Phonon(2))) == Hermitian([0 0 -im 0; 0 0 0 -im; im 0 0 0; 0 im 0 0])

    contentnames(AbstractTBA) == (:H, :commutator)
    contentnames(TBA) == (:lattice, :H, :commutator)
end

@time @testset "TBA" begin
    lattice = Lattice(:Square,
        [Point(PID(1), [0.0, 0.0])],
        vectors=[[1.0, 0.0], [0.0, 1.0]],
        neighbors=1
        )
    hilbert = Hilbert(PID(1)=>Fock{:f}(1, 1, 2))
    t = Hopping(:t, 1.0, 1)
    μ = Onsite(:μ, 0.0, modulate=true)

    tba = TBA(lattice, hilbert, (t, μ))
    @test kind(tba) == kind(typeof(tba)) == TBAKind(:TBA)
    @test valtype(tba) == valtype(typeof(tba)) == Float64
    @test dimension(tba) == 1
    @test Parameters(tba) == (t=1.0, μ=0.0)

    m = matrix(tba)
    @test size(m) == (1, 1)
    @test m[1, 1] == m.H[1, 1]
    @test ishermitian(m) == ishermitian(typeof(m)) == true

    A(t, μ; k) = hcat(2t*cos(k[1])+2t*cos(k[2])+μ)
    tbaₐ = TBA(lattice, A, (t=1.0, μ=0.0))
    path = ReciprocalPath(lattice.reciprocals, rectangle"Γ-X-M-Γ", length=8)
    for kv in pairs(path)
        m = matrix(tba; kv...)
        mₐ = matrix(tbaₐ; kv...)
        @test m.H ≈ mₐ.H ≈ Hermitian(A(1.0, 0.0; kv...))
    end
    update!(tba, μ=0.5)
    update!(tbaₐ, μ=0.5)
    for kv in pairs(path)
        m = matrix(tba; kv...)
        mₐ = matrix(tbaₐ; kv...)
        @test m.H ≈ mₐ.H ≈ Hermitian(A(1.0, 0.5; kv...))
    end

    Δ = Pairing(:Δ, Complex(0.1), 1, amplitude=bond->exp(im*azimuth(rcoord(bond))))
    bdg = TBA(lattice, hilbert, (t, μ, Δ))
    @test kind(bdg) == kind(typeof(bdg)) == TBAKind(:BdG)
    @test valtype(bdg) == valtype(typeof(bdg)) == Complex{Float64}
    @test Parameters(bdg) == (t=1.0, μ=0.5, Δ=Complex(0.1))

    A(t, μ, Δ; k) = [2t*cos(k[1])+2t*cos(k[2])+μ 2im*Δ*sin(k[1])+2Δ*sin(k[2]); -2im*Δ*sin(k[1])+2Δ*sin(k[2]) -2t*cos(k[1])-2t*cos(k[2])-μ]
    bdgₐ = TBA(lattice, A, (t=1.0, μ=0.5, Δ=0.1))
    path = ReciprocalPath(lattice.reciprocals, rectangle"Γ-X-M-Γ", length=8)
    for kv in pairs(path)
        m = matrix(bdg; kv...)
        mₐ = matrix(bdgₐ; kv...)
        @test m.H ≈ mₐ.H ≈ Hermitian(A(1.0, 0.5, 0.1; kv...))
    end

    samplesets = [
        SampleNode(lattice.reciprocals, [0.0, 0.0], [1], [-4.0]),
        SampleNode(lattice.reciprocals, [0.5, 0.5], [1], [+4.0])
    ]
    op = optimize!(tba, samplesets)[2]
    @test isapprox(op.minimizer, [-1.0, 0.0], atol=10^-10)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 60b53bf6-dfa2-4fc5-b725-9e4239f7132c
using Colors, SpectralVis

# ╔═╡ 23835930-c0e0-11eb-21a7-49aedc5bc4b0
md"CIE1960 UCS"

# ╔═╡ 148c4957-ed02-4023-841f-9c5e82a6f796
plotly()

# ╔═╡ 144388e1-7642-4aba-b402-04a2114d9268
function cie60(c::XYZ)
	U = 2 * c.x / 3
	V = c.y
	W = 0.5 * (-c.x + 3c.y + c.z)
	u = U / (U+V+W)
	v = V / (U+V+W)
	u,v
end

# ╔═╡ f55e9217-f9b6-43e5-b29e-aa532f4adb5c
cie60(XYZ(.1, .5, .1))

# ╔═╡ 50a68b35-6e4d-4cd2-bf77-7aff52afeb83
begin
	temp = 6500
	meg = 1e7
	inv = meg / temp

	range = map(x -> 10000000/x, collect(10:500:100000))
	vals = collect(range)
	u = zeros(length(vals))
	v = zeros(length(vals))

	for k in 1:length(vals)
		t = vals[k]
		wp = blackbody_whitepoint(SPECENV,t)
		a,b = cie60(wp)
		u[k] = a
		v[k] = b
	end

	plot(xlimits=(0.15,0.25))
	scatter!(u,v,ratio = 1.0, color = :orange, markersize = 2, label = "knots")

		spl = cubicspline(u,v,:natural)
	sin = interp(spl,u[1]:0.001:u[end])
	plot!(sin[1],sin[2], color = :red, label = "spline")

	fine = map(x -> x, collect(100:10:20000))
	f = zeros(length(fine))
	g = zeros(length(fine))

	for l in 1:length(fine)
		te = fine[l]
		wpo = blackbody_whitepoint(SPECENV,te)
		c,d = cie60(wpo)
		f[l] = c
		g[l] = d
	end

	scatter!(f,g, color = :green, markersize = 1, markerstyle = :cross, label="truth")
end

# ╔═╡ ef6cd9f0-8416-472a-a4bb-43d5265d5757

# ╔═╡ acf3ed94-54cc-408f-bdae-c2eb16c9bbb7
range

# ╔═╡ b09c644b-449f-4756-aa82-c9aeb6e2265c
length(range)

# ╔═╡ Cell order:
# ╟─23835930-c0e0-11eb-21a7-49aedc5bc4b0
# ╠═60b53bf6-dfa2-4fc5-b725-9e4239f7132c
# ╠═148c4957-ed02-4023-841f-9c5e82a6f796
# ╠═144388e1-7642-4aba-b402-04a2114d9268
# ╠═f55e9217-f9b6-43e5-b29e-aa532f4adb5c
# ╠═50a68b35-6e4d-4cd2-bf77-7aff52afeb83
# ╠═ef6cd9f0-8416-472a-a4bb-43d5265d5757
# ╠═acf3ed94-54cc-408f-bdae-c2eb16c9bbb7
# ╠═b09c644b-449f-4756-aa82-c9aeb6e2265c


text:
# Discrete operators types and methods.
# Based on the LinearMaps.jl package.

############
# AbstractDiscreteOp
############

"""
    abstract type AbstractDiscreteOp

A linear operator that can be evaluated using its forward map. It also
supports additions and compositions with others `::AbstractDiscreteOp`.
"""
abstract type AbstractDiscreteOp end

const DiscreteOpTuple = Tuple{Vararg{AbstractDiscreteOp}}

Base.:+(d::AbstractDiscreteOp) = d
Base.:+(d::AbstractDiscreteOp,_) = abstractmethod(d)
Base.:+(_,d::AbstractDiscreteOp) = abstractmethod(d)
Base.:-(d::AbstractDiscreteOp,_) = abstractmethod(d)
Base.:-(_,d::AbstractDiscreteOp) = abstractmethod(d)
Base.:*(d::AbstractDiscreteOp,_) = abstractmethod(d)
Base.:*(_,d::AbstractDiscreteOp) = abstractmethod(d)
Base.size(d::AbstractDiscreteOp) = abstractmethod(d)
Base.size(d::AbstractDiscreteOp,i::Integer) = size(d)[i]

"""
    materialize(d::AbstractDiscreteOp)

Returns the explicit matrix representation of the underlying operator in
`d::AbstractDiscreteOp`. This method can be expensive due to
matrix additions and multiplications.
"""
materialize(d::AbstractDiscreteOp) = abstractmethod(d)
_materialize(d::AbstractDiscreteOp,x) = materialize(d)*x

_check_dim_mul(A,B) = @assert(size(A,2) == size(B,1))
_check_dim_sum(A,B) = @assert(size(A) == size(B))

############
# UniformScalingDiscreteOp
############
struct UniformScalingDiscreteOp{T<:Number} <: AbstractDiscreteOp
    λ::T
    s::Int64  # size
end
Base.size(d::UniformScalingDiscreteOp) = (d.s,d.s)

function Base.:*(d::UniformScalingDiscreteOp,x::AbstractVecOrMat)
    _check_dim_mul(d,x)
    return d.λ*x
end
function Base.:*(d1::UniformScalingDiscreteOp,d2::UniformScalingDiscreteOp)
    _check_dim_mul(d1,d2)
    return UniformScalingDiscreteOp(d1.λ*d2.λ, size(d1,1))
end
function Base.:*(d::UniformScalingDiscreteOp,α::Number)
    return UniformScalingDiscreteOp(d.λ*α, size(d,2))
end
function Base.:*(α::Number,d::UniformScalingDiscreteOp)
    return UniformScalingDiscreteOp(α*d.λ, size(d,1))
end

materialize(d::UniformScalingDiscreteOp) = UniformScaling(d.λ)

############
# DiscreteOp
############
struct DiscreteOp{D} <: AbstractDiscreteOp
    op::D
end
DiscreteOp(d::AbstractDiscreteOp) = d
Base.size(d::DiscreteOp) = size(d.op)
Base.:*(d::DiscreteOp,x::AbstractVecOrMat) = d.op*x
materialize(d::DiscreteOp) = d.op

############
# CompositeDiscreteOp
############
struct CompositeDiscreteOp{D<:DiscreteOpTuple} <: AbstractDiscreteOp
    maps::D
    function CompositeDiscreteOp(maps::D) where D<:DiscreteOpTuple
        @assert length(maps) > 1
        return new{D}(maps)
    end
end
CompositeDiscreteOp(d1::AbstractDiscreteOp, d2::AbstractDiscreteOp) = CompositeDiscreteOp((d1,d2))
CompositeDiscreteOp(d1::AbstractDiscreteOp, d2::CompositeDiscreteOp) = CompositeDiscreteOp((d1,d2.maps...))
CompositeDiscreteOp(d1::CompositeDiscreteOp, d2::AbstractDiscreteOp) = CompositeDiscreteOp((d1.maps...,d2))
CompositeDiscreteOp(d1::CompositeDiscreteOp, d2::CompositeDiscreteOp) = CompositeDiscreteOp((d1.maps...,d2.maps...))

Base.size(d::CompositeDiscreteOp) = (size(first(d.maps),1), size(last(d.maps),2))

function Base.:*(d::CompositeDiscreteOp,x::AbstractVecOrMat)
    _check_dim_mul(d,x)
    # evaluate product from right to left
    y = last(d.maps)*x
    length(d.maps) == 2 && return first(d.maps)*y
    return CompositeDiscreteOp(d.maps[1:end-1])*y
end
function Base.:*(d1::AbstractDiscreteOp, d2::AbstractDiscreteOp)
    _check_dim_mul(d1,d2)
    return CompositeDiscreteOp(d1, d2)
end
function Base.:*(α::Number, d::AbstractDiscreteOp)
    return CompositeDiscreteOp(UniformScalingDiscreteOp(α,size(d,1)), d)
end
function Base.:*(d::AbstractDiscreteOp, α::Number)
    return CompositeDiscreteOp(d, UniformScalingDiscreteOp(α,size(d,2)))
end
Base.:*(u::UniformScaling, d::AbstractDiscreteOp) = u.λ*d
Base.:*(d::AbstractDiscreteOp, u::UniformScaling) = d*u.λ
function Base.:-(d::AbstractDiscreteOp)
    return CompositeDiscreteOp(UniformScalingDiscreteOp(-1,size(d,1)), d)
end

function materialize(d::CompositeDiscreteOp)
    return _materialize(d, true)
end
function _materialize(d::CompositeDiscreteOp, x)
    # evaluate operators from right to left
    m1 = materialize(d.maps[end])
    y = _mymul(m1,x)   # m1*x
    if length(d.maps) == 2
        m2 = materialize(d.maps[1])
        return _mymul(m2,y)  # m2*y
    end
    return _materialize(CompositeDiscreteOp(d.maps[1:end-1]), y)
end
function _mymul(A, B)
    # this function patches a bug that
    # occurs when a Diagonal{<:SMatrix}
    # is multiplied by a Matrix{<:SMatrix}
    # with SMatrices of different sizes
    if A isa Diagonal{<:SMatrix} && B isa Matrix{<:SMatrix}
        return A.diag.*B
    elseif A isa Matrix{<:SMatrix} && B isa Diagonal{<:SMatrix}
        return A.*permutedims(B.diag)
    else
        return A*B
    end
end

############
# LinearCombinationDiscreteOp
############
struct LinearCombinationDiscreteOp{D<:DiscreteOpTuple} <: AbstractDiscreteOp
    maps::D
    function LinearCombinationDiscreteOp(maps::D) where D<:DiscreteOpTuple
        @assert length(maps) > 1
        return new{D}(maps)
    end
end
LinearCombinationDiscreteOp(d1::AbstractDiscreteOp, d2::AbstractDiscreteOp) = LinearCombinationDiscreteOp((d1,d2))
LinearCombinationDiscreteOp(d1::AbstractDiscreteOp, d2::LinearCombinationDiscreteOp) = LinearCombinationDiscreteOp((d1,d2.maps...))
LinearCombinationDiscreteOp(d1::LinearCombinationDiscreteOp, d2::AbstractDiscreteOp) = LinearCombinationDiscreteOp((d1.maps...,d2))
LinearCombinationDiscreteOp(d1::LinearCombinationDiscreteOp, d2::LinearCombinationDiscreteOp) = LinearCombinationDiscreteOp((d1.maps...,d2.maps...))

Base.size(d::LinearCombinationDiscreteOp) = size(first(d.maps))

function Base.:*(d::LinearCombinationDiscreteOp,x::AbstractVecOrMat)
    _check_dim_mul(d,x)
    # evaluate product from left to right
    y = first(d.maps)*x
    for n in 2:length(d.maps)
        # TODO: use mul! instead
        y += d.maps[n]*x
    end
    return y
end
function Base.:+(d1::AbstractDiscreteOp, d2::AbstractDiscreteOp)
    _check_dim_sum(d1,d2)
    return LinearCombinationDiscreteOp(d1, d2)
end
function Base.:+(u::UniformScaling, d::AbstractDiscreteOp)
    @assert size(d,1) == size(d,2)
    return UniformScalingDiscreteOp(u.λ,size(d,1)) + d
end
function Base.:+(d::AbstractDiscreteOp, u::UniformScaling)
    @assert size(d,1) == size(d,2)
    return d + UniformScalingDiscreteOp(u.λ,size(d,2))
end
Base.:-(u::UniformScaling, d::AbstractDiscreteOp) = u+(-d)
Base.:-(d::AbstractDiscreteOp, u::UniformScaling) = d+(-u)
Base.:-(d1::AbstractDiscreteOp, d2::AbstractDiscreteOp) = d1+(-d2)

function materialize(d::LinearCombinationDiscreteOp)
    # find a map that is not an UniformScalingDiscreteOp
    map_index = 1
    for _ in eachindex(d.maps)
        !isa(d.maps[map_index], UniformScalingDiscreteOp) && break
        map_index += 1
    end
    msg = """Cannot materialize a LinearCombinationDiscreteOp
    whose maps are all UniformScalingDiscreteOp. Not implemented."""
    @assert (map_index ≤ length(d.maps)) msg
    # materialize map and
    # evaluate the sum from left to right
    y = materialize(d.maps[map_index])
    for n in eachindex(d.maps)
        n == map_index && continue
        y += materialize(d.maps[n])
    end
    return y
end

############
# GMRES and solvers
############
"""
    struct DiscreteOpGMRES{D<:AbstractDiscreteOp,T<:Number,V}

A wrapper around 'op::AbstractDiscreteOp' to use in conjunction with
the `IterativeSolvers.gmres!` method.
"""
struct DiscreteOpGMRES{D<:AbstractDiscreteOp,T<:Number,V}
    # T: scalar matrix element type
    # V: element type of the vectors that `op::AbstractDiscreteOp`
    # receives as input (e.g. SVector{3,Float64} or ComplexF64)
    op::D
    s::Int64  # scalar matrix size
end
function DiscreteOpGMRES(op::AbstractDiscreteOp, σ::AbstractVector{V}) where V
    @assert size(op,1) == size(op,2)
    @assert size(op,2) == length(σ)
    if V <: Number
        s = length(σ)  # size
        return DiscreteOpGMRES{typeof(op),V,V}(op, s)
    elseif V <: SVector
        T = eltype(V)
        s = length(σ) * length(V)  # size
        return DiscreteOpGMRES{typeof(op),T,V}(op, s)
    else
        notimplemented()
    end
end
Base.size(g::DiscreteOpGMRES) = (g.s, g.s)
Base.size(g::DiscreteOpGMRES, i) = size(g)[i]
Base.eltype(::DiscreteOpGMRES{D,T}) where {D,T} = T
function LinearAlgebra.mul!(yvec::AbstractVector{T}, g::DiscreteOpGMRES{D,T,V}, xvec::AbstractVector{T}) where {D,T,V}
    if T === V
        yvec .= g.op*xvec
    else
        x = reinterpret(V, xvec)
        yvec .= reinterpret(T, g.op*x)
    end
    return yvec
end

function IterativeSolvers.gmres!(σ::Density{V},A::AbstractDiscreteOp,μ::Density{V},args...;kwargs...) where V
    log = haskey(kwargs,:log) ? kwargs[:log] : false
    Aop = DiscreteOpGMRES(A, μ)
    if V <: Number
        if log
            vals,hist = gmres!(σ.vals,Aop,μ.vals,args...;kwargs...)
            return σ,hist
        else
            vals = gmres!(σ.vals,Aop,μ.vals,args...;kwargs...)
            return σ
        end
    elseif V <: SVector
        σ_vec = reinterpret(eltype(V),σ.vals)
        μ_vec = reinterpret(eltype(V),μ.vals)
        if log
            vals,hist = gmres!(σ_vec,Aop,μ_vec,args...;kwargs...)
            return σ,hist
        else
            vals = gmres!(σ_vec,Aop,μ_vec,args...;kwargs...)
            return σ
        end
    else
        notimplemented()
    end
end

function Base.:\(A::AbstractDiscreteOp,σ::Density{V}) where V
    Amat = materialize(A)  # assemble the full matrix
    if V <: Number
        @assert eltype(Amat) === V
        vals = Amat \ σ.vals
        return Density(vals, σ.mesh)
    elseif V <: SVector
        T = eltype(V)
        N = length(V)
        @assert eltype(Amat) <: SMatrix{N,N,T}
        Ascalar = blockmatrix_to_matrix(Amat)
        σ_vec = reinterpret(T, σ.vals)
        vals_vec = Ascalar \ σ_vec
        vals = reinterpret(V, vals_vec) |> collect
        return Density(vals, σ.mesh)
    else
        notimplemented()
    end
end

############
# show
############
function _show_size(D)
    return "$(size(D,1))×$(size(D,2))"
end

function _show_typeof(D)
    return _show_size(D) * ' ' * string(typeof(D))
end
function _show_typeof(D::Number)
    return string(typeof(D))
end
function _show_typeof(D::Nystrom.AbstractDiscreteOp)
    return _show_size(D) * ' ' * split(string(typeof(D)),'{')[1]
end

function _show_content(D::Nystrom.AbstractDiscreteOp, i)
    abstractmethod(D)
end
function _show_content(D::Nystrom.UniformScalingDiscreteOp, i)
    content = _show_typeof(D) * " with scale: $(D.λ)" * '\n' * ' '^i * _show_typeof(D.λ)
    return content
end
function _show_content(D::Nystrom.DiscreteOp, i)
    content = _show_typeof(D) * " of" *'\n' * ' '^i * _show_typeof(D.op)
    return content
end
function _show_content(D::Union{Nystrom.CompositeDiscreteOp,Nystrom.LinearCombinationDiscreteOp}, i)
    content = _show_typeof(D) * " with $(length(D.maps)) maps:"
    for m in D.maps
        content *= '\n' * ' '^i * _show_content(m, i+2)
    end
    return content
end

Base.show(io::IO, D::Nystrom.AbstractDiscreteOp) = print(io, _show_content(D, 2))

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -37.14009556780016 <= q <= 48.365016418887166)

                     add_NL_constraint(m, :(sigmoid(0.24902435571631587 + -0.02312729681609138*sigmoid(-0.46870412880299783 + 0.1700038462852933*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.44738483682757924*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.5112184633122538*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.17714072857386487*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9529511507705823*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.5418072450548097*sigmoid(0.5019192867350393 + 0.2997113620022893*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.5721484689942011*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.4073624372378828*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + -0.07630883358086304*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.3103910161241421*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.9789829058919253*sigmoid(-0.652369829242545 + -0.998242239199818*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.041041478696029365*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.4348093907280841*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.7773011179238121*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7537599741964742*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.13446361488468206*sigmoid(-0.40141742757451926 + 0.7289092837988314*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + 0.16077592157977172*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.6378484913808045*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9326240829024273*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9638370498324353*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.3069132295129271*sigmoid(0.17680227033398044 + 0.8293506390949057*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.09356004658531791*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.16358082594109513*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9029174038261001*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7318887196209252*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3])))) + sigmoid(-0.017552134571214406 + -0.41004797399293524*sigmoid(-0.46870412880299783 + 0.1700038462852933*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.44738483682757924*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.5112184633122538*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.17714072857386487*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9529511507705823*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.47019448283071075*sigmoid(0.5019192867350393 + 0.2997113620022893*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.5721484689942011*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.4073624372378828*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + -0.07630883358086304*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.3103910161241421*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.35795867979675133*sigmoid(-0.652369829242545 + -0.998242239199818*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.041041478696029365*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.4348093907280841*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.7773011179238121*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7537599741964742*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.44715886126593407*sigmoid(-0.40141742757451926 + 0.7289092837988314*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + 0.16077592157977172*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.6378484913808045*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9326240829024273*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9638370498324353*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.22638221307132733*sigmoid(0.17680227033398044 + 0.8293506390949057*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.09356004658531791*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.16358082594109513*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9029174038261001*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7318887196209252*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3])))) + sigmoid(-0.4029592795849051 + -0.6197606146665708*sigmoid(-0.46870412880299783 + 0.1700038462852933*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.44738483682757924*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.5112184633122538*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.17714072857386487*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9529511507705823*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.8464388391219866*sigmoid(0.5019192867350393 + 0.2997113620022893*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.5721484689942011*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.4073624372378828*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + -0.07630883358086304*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.3103910161241421*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.6262663939568953*sigmoid(-0.652369829242545 + -0.998242239199818*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.041041478696029365*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.4348093907280841*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.7773011179238121*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7537599741964742*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.06090941414894013*sigmoid(-0.40141742757451926 + 0.7289092837988314*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + 0.16077592157977172*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.6378484913808045*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9326240829024273*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9638370498324353*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.6699649674430836*sigmoid(0.17680227033398044 + 0.8293506390949057*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.09356004658531791*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.16358082594109513*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9029174038261001*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7318887196209252*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3])))) + sigmoid(0.6579858782021573 + -0.21911632564812766*sigmoid(-0.46870412880299783 + 0.1700038462852933*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.44738483682757924*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.5112184633122538*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.17714072857386487*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9529511507705823*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.7292607357990395*sigmoid(0.5019192867350393 + 0.2997113620022893*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.5721484689942011*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.4073624372378828*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + -0.07630883358086304*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.3103910161241421*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.23882038501163416*sigmoid(-0.652369829242545 + -0.998242239199818*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.041041478696029365*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.4348093907280841*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.7773011179238121*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7537599741964742*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.6855544316424598*sigmoid(-0.40141742757451926 + 0.7289092837988314*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + 0.16077592157977172*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.6378484913808045*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9326240829024273*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9638370498324353*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.2868272308560216*sigmoid(0.17680227033398044 + 0.8293506390949057*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.09356004658531791*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.16358082594109513*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9029174038261001*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7318887196209252*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3])))) + sigmoid(-0.10147744083820065 + -0.7175110849623545*sigmoid(-0.46870412880299783 + 0.1700038462852933*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.44738483682757924*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.5112184633122538*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.17714072857386487*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9529511507705823*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.9583278917345623*sigmoid(0.5019192867350393 + 0.2997113620022893*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.5721484689942011*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.4073624372378828*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + -0.07630883358086304*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.3103910161241421*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.7313133539296679*sigmoid(-0.652369829242545 + -0.998242239199818*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.041041478696029365*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.4348093907280841*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.7773011179238121*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7537599741964742*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.008349755673621395*sigmoid(-0.40141742757451926 + 0.7289092837988314*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + 0.16077592157977172*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.6378484913808045*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9326240829024273*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9638370498324353*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.13663535205769417*sigmoid(0.17680227033398044 + 0.8293506390949057*sigmoid(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.09356004658531791*sigmoid(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.16358082594109513*sigmoid(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9029174038261001*sigmoid(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7318887196209252*sigmoid(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using BenchmarkTools
using InteractiveUtils

struct MyTypeAny
    x
    y
end

struct MyTypeConcrete
    x::Vector{Int64}
    y::Vector{Float64}
end

struct MyTypeAbstract
    x::Vector{Integer}
    y::Vector{Real}
end

struct MyTypeParametric{T1<:Integer,T2<:Real}
    x::Vector{T1}
    y::Vector{T2}
end

N=100
@benchmark mt = MyTypeAny(rand(Int64,N), rand(Float64, N))
@benchmark mtc = MyTypeConcrete(rand(Int64,N), rand(Float64, N))
@benchmark mta = MyTypeAbstract(rand(Int64,N), rand(Float64, N))
@benchmark mtp = MyTypeParametric(rand(Int64,N),rand(Float64, N))

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using LaTeXStrings

@cell "heatmap_with_labels" begin
    fig = Figure(resolution = (600, 600))
    ax = fig[1, 1] = Axis(fig)
    values = RNG.rand(10, 10)

    heatmap!(ax, values)

    text!(ax,
        string.(round.(vec(values'), digits = 2)),
        position = [Point2f(x, y) for x in 1:10 for y in 1:10],
        align = (:center, :center),
        color = ifelse.(vec(values') .< 0.3, :white, :black),
        textsize = 12)
    fig
end

@cell "data space" begin
    pos = [Point2f(0, 0), Point2f(10, 10)]
    fig = text(
        ["0 is the ORIGIN of this", "10 says hi"],
        position = pos,
        axis = (aspect = DataAspect(),),
        markerspace = :data,
        align = (:center, :center),
        textsize = 2)
    scatter!(pos)
    fig
end

@cell "single_strings_single_positions" begin
    scene = Scene(camera = campixel!, resolution = (800, 800))

    points = [Point(x, y) .* 200 for x in 1:3 for y in 1:3]
    scatter!(scene, points, marker = :circle, markersize = 10px)

    i = 1
    for halign in (:right, :center, :left), valign in (:top, :center, :bottom)

        for rotation in (-pi/6, 0.0, pi/6)
            text!(scene, string(halign) * "/" * string(valign) *
                    " " * string(round(rad2deg(rotation), digits = 0)) * "°",
                color = (:black, 0.5),
                position = points[i],
                align = (halign, valign),
                rotation = rotation)
        end
        i += 1
    end
    scene
end

@cell "multi_strings_multi_positions" begin
    scene = Scene(camera = campixel!, resolution = (800, 800))

    angles = (-pi/6, 0.0, pi/6)
    points = [Point(x, y) .* 200 for x in 1:3 for y in 1:3 for angle in angles]
    aligns = [(halign, valign) for halign in
        (:right, :center, :left) for valign in (:top, :center, :bottom) for rotation in angles]
    rotations = [rotation for _ in
        (:right, :center, :left) for _ in (:top, :center, :bottom) for rotation in angles]

    strings = [string(halign) * "/" * string(valign) *
        " " * string(round(rad2deg(rotation), digits = 0)) * "°"
            for halign in (:right, :center, :left)
            for valign in (:top, :center, :bottom)
            for rotation in angles]

    scatter!(scene, points, marker = :circle, markersize = 10px)

    text!(scene, strings, position = points, align = aligns, rotation = rotations,
        color = [(:black, alpha) for alpha in LinRange(0.3, 0.7, length(points))])

    scene
end

@cell "single_strings_single_positions_justification" begin
    scene = Scene(camera = campixel!, resolution = (800, 800))

    points = [Point(x, y) .* 200 for x in 1:3 for y in 1:3]
    scatter!(scene, points, marker = :circle, markersize = 10px)

    symbols = (:left, :center, :right)

    for ((justification, halign), point) in zip(Iterators.product(symbols, symbols), points)

        t = text!(scene, "a\nshort\nparagraph",
            color = (:black, 0.5),
            position = point,
            align = (halign, :center),
            justification = justification)

        bb = boundingbox(t)
        wireframe!(scene, bb, color = (:red, 0.2))
    end

    for (p, al) in zip(points[3:3:end], (:left, :center, :right))
        text!(scene, "align :" * string(al), position = p .+ (0, 80),
            align = (:center, :baseline))
    end

    for (p, al) in zip(points[7:9], (:left, :center, :right))
        text!(scene, "justification\n:" * string(al), position = p .+ (80, 0),
            align = (:center, :top), rotation = pi/2)
    end

    scene
end

@cell "multi_boundingboxes" begin
    scene = Scene(camera = campixel!, resolution = (800, 800))

    t1 = text!(scene,
        fill("makie", 4),
        position = [(200, 200) .+ 60 * Point2f(cos(a), sin(a)) for a in pi/4:pi/2:7pi/4],
        rotation = pi/4:pi/2:7pi/4,
        align = (:left, :center),
        textsize = 30,
        markerspace = :data
    )

    wireframe!(scene, boundingbox(t1), color = (:blue, 0.3))

    t2 = text!(scene,
        fill("makie", 4),
        position = [(200, 600) .+ 60 * Point2f(cos(a), sin(a)) for a in pi/4:pi/2:7pi/4],
        rotation = pi/4:pi/2:7pi/4,
        align = (:left, :center),
        textsize = 30,
        markerspace = :pixel
    )

    wireframe!(scene, boundingbox(t2), color = (:red, 0.3))

    scene
end

@cell "single_boundingboxes" begin
    scene = Scene(camera = campixel!, resolution = (800, 800))

    for a in pi/4:pi/2:7pi/4

        t = text!(scene,
            "makie",
            position = (200, 200) .+ 60 * Point2f(cos(a), sin(a)),
            rotation = a,
            align = (:left, :center),
            textsize = 30,
            markerspace = :data
        )

        wireframe!(scene, boundingbox(t), color = (:blue, 0.3))

        t2 = text!(scene,
            "makie",
            position = (200, 600) .+ 60 * Point2f(cos(a), sin(a)),
            rotation = a,
            align = (:left, :center),
            textsize = 30,
            markerspace = :pixel
        )

        # these boundingboxes should be invisible because they only enclose the anchor
        wireframe!(scene, boundingbox(t2), color = (:red, 0.3))

    end
    scene
end

@cell "text_in_3d_axis" begin
    text(
        fill("Makie", 7),
        rotation = [i / 7 * 1.5pi for i in 1:7],
        position = [Point3f(0, 0, i/2) for i in 1:7],
        color = [cgrad(:viridis)[x] for x in LinRange(0, 1, 7)],
        align = (:left, :baseline),
        textsize = 1,
        markerspace = :data
    )
end

@cell "empty_lines" begin
    scene = Scene(camera = campixel!, resolution = (800, 800))

    t1 = text!(scene, "Line1\nLine 2\n\nLine4",
        position = (200, 400), align = (:center, :center), markerspace = :data)

    wireframe!(scene, boundingbox(t1), color = (:red, 0.3))

    t2 = text!(scene, "\nLine 2\nLine 3\n\n\nLine6\n\n",
        position = (400, 400), align = (:center, :center), markerspace = :data)

    wireframe!(scene, boundingbox(t2), color = (:blue, 0.3))

    scene
end

@cell "3D screenspace annotations" begin
    positions = RNG.rand(Point3f, 10)
    fig, ax, p = meshscatter(positions, color=:white)
    text!(
        fill("Annotation", 10),
        position = positions,
        align = (:center, :center),
        textsize = 20,
        markerspace = :pixel,
        overdraw=false)
    fig
end

@cell "Text offset" begin
    f = Figure(resolution = (1000, 1000))
    barplot(f[1, 1], 3:5)
    text!("bar 1", position = (1, 3), offset = (0, 10), align = (:center, :baseline))
    text!(["bar 2", "bar 3"], position = [(2, 4), (3, 5)],
        offset = [(0, -10), (0, -20)],
        align = (:center, :top), color = :white)

    scatter(f[1, 2], Point2f(0, 0))
    text!("hello", position = (0, 0), offset = (40, 0), align = (:left, :center))
    text!("hello", position = (0, 0), offset = (40, 0), align = (:left, :center),
        rotation = -pi/4)
    text!("hello", position = (0, 0), offset = (40, 0), align = (:left, :center),
        rotation = pi/4)

    scatter(f[2, 1], Point2f[(0, 0), (10, 0), (20, 10)])
    text!("ABC", markerspace = :data, offset = (0, 0), color = (:red, 0.3), align = (:left, :baseline))
    text!("ABC", markerspace = :data, offset = (10, 0), color = (:green, 0.3), align = (:left, :baseline))
    text!("ABC", markerspace = :data, offset = (20, 10), color = (:blue, 0.3), align = (:left, :baseline))

    LScene(f[2, 2], show_axis=false)
    scatter!(Point3f[(0, 0, 0), (2, 2, 2)])
    text!("hello", position = Point3f(1, 1, 1), offset = (10, 10))

    f
end

@cell "Log10 text" begin
    barplot([1, 10, 100], fillto = 0.1, axis = (yscale = log10,))
    text!(["bar 1", "bar 2", "bar 3"], position = [(1, 1), (2, 10), (3, 100)],
        offset = (0, -10), color = :white, align = (:center, :top))
    tightlimits!(current_axis(), Bottom())
    current_figure()
end

@cell "latex strings" begin
    f, ax , l = lines(cumsum(RNG.randn(1000)),
        axis = (
                title = L"\sum_k{x y_k}",
                xlabel = L"\lim_{x →\infty} A^j v_{(a + b)_k}^i \sqrt{23.5} x!= \sqrt{\frac{1+6}{4+a+g}}\int_{0}^{2π} \sin(x) dx",
                ylabel = L"x + y - sin(x) × tan(y) + \sqrt{2}",
            ),
        figure = (fontsize = 18,)
    )
    text!(L"\int_{0}^{2π} \sin(x) dx", position = (500, 0))
    Legend(f[1, 2], [l, l, l], [L"\sum{xy}", L"a\int_0^5x^2+2ab", L"||x-y||^2"])
    f
end

@cell "latex simple" begin
    s = Scene(camera = campixel!)
    t = text!(s,
        L"\sqrt{2}",
        position = (50, 50),
        rotation = pi/2,
        markerspace = :data)
    s
end

@cell "latex bb" begin
    s = Scene(camera = campixel!)
    t = text!(s,
        L"\int_0^5x^2+2ab",
        position = Point2f(50, 50),
        rotation = 0.0,
        markerspace = :data)
    wireframe!(s, boundingbox(t))
    s
end

@cell "latex updates" begin
    s = Scene(camera = campixel!)
    st = Stepper(s)
    textnode = Observable([L"\int_0^5x^2+2ab", L"\int_0^5x^2+2ab"])
    posnode = Observable(Point2f[(50, 50), (100, 100)])

    t = text!(s,
        textnode,
        position = posnode,
        rotation = 0.0,
        markerspace = :data)

    Makie.step!(st)
    ## change lengths
    textnode.val = push!(textnode[], L"\int_0^5x^2+2ab")
    posnode[] = push!(posnode[], Point2f(150, 150))
    Makie.step!(st)
    st
end

@cell "update annotation style" begin
    s = Scene(camera = campixel!)
    st = Stepper(s)
    textposnode = Observable([
        (L"\int_0^5x^2+2ab", Point2f(50, 50)),
        (L"\int_0^5x^2+2ab", Point2f(100, 100)),
    ])

    t = text!(s,
        textposnode,
        markerspace = :data)

    Makie.step!(st)
    ## change lengths
    textposnode[] = push!(textposnode[], (L"\int_0^5x^2+2ab", Point2f(150, 150)))
    Makie.step!(st)
    st
end

@cell "latex ticks" begin
    lines(0..25, x -> 4 * sin(x) / (cos(3x) + 4), figure = (fontsize = 25,),
        axis = (
            xticks = (0:10:20, [L"10^{-3.5}", L"10^{-4.5}", L"10^{-5.5}"]),
            yticks = ([-1, 0, 1], [L"\sum_%$i{xy}" for i in 1:3]),
            yticklabelrotation = pi/8,
            title = L"\int_0^1{x^2}",
            xlabel = L"\sum_k{x_k ⋅ y_k}",
            ylabel = L"\int_a^b{\sqrt{abx}}"
        ),
    )
end

@cell "dynamic latex ticks" begin
    lines(0..25, x -> 4 * sin(x) / (cos(3x) + 4),
        figure = (fontsize = 16,),
        axis = (xtickformat = (xs -> [L"e^{\sqrt{%$x}}+\sum" for x in xs]), ))
end


text:
export backslash_solve!, build_kkt_system

function reset_value_and_add_vars!(x::AffineOrConstant, unique_vars_map::Dict{UInt64, AffineExpr})
  if x.head != :constant
    x.value = nothing
    if x.head == :variable
      unique_vars_map[object_id(x)] = x
    else
      for child in x.children
        reset_value_and_add_vars!(child, unique_vars_map)
      end
    end
  end
end

function reset_values_and_get_vars!(p::Problem)
  p.status = "not yet solved"
  p.optval = nothing
  unique_vars_map = Dict{UInt64, AffineExpr}()
  p.objective.value = nothing
  for affine in p.objective.affines
    reset_value_and_add_vars!(affine, unique_vars_map)
  end
  for constraint in p.constraints
    reset_value_and_add_vars!(constraint.canon_form, unique_vars_map)
  end
  return unique_vars_map
end

function get_var_ranges_and_num_vars(unique_vars_map::Dict{UInt64, AffineExpr})
  vars_to_ranges_map = Dict{UInt64, Tuple{Int64, Int64}}()
  num_vars = 0
  for (id, var) in unique_vars_map
    sz = var.size[1] * var.size[2]
    vars_to_ranges_map[id] = (num_vars + 1, num_vars + sz)
    num_vars += sz
  end
  return vars_to_ranges_map, num_vars
end

function get_num_rows(affines::Array{AffineExpr})
  num_rows = 0
  for affine in affines
    num_rows += affine.size[1] * affine.size[2]
  end
  return num_rows
end

function coalesce_affine_exprs(vars_to_ranges_map::Dict{UInt64, Tuple{Int64, Int64}}, num_vars::Int64, affines::Array{AffineExpr})
  num_rows = get_num_rows(affines)
  coefficient = spzeros(num_rows, num_vars)
  constant = spzeros(num_rows, 1)
  row = 0
  for affine in affines
    row_sz = affine.size[1] * affine.size[2]
    for (var, coeff) in affine.vars_to_coeffs_map
      col_range = vars_to_ranges_map[var]
      coefficient[row + 1 : row + row_sz, col_range[1] : col_range[2]] = coeff.value
    end
    constant[row + 1 : row + row_sz] = affine.constant.value
    row += row_sz
  end
  return coefficient, constant
end

function build_kkt_system(A, b, C, d)
  if size(A, 1) == 0
    num_rows = size(C, 1)
    num_cols = size(C, 2)
  else
    num_rows = num_cols = size(A, 2) + size(C, 1)
  end

  coefficient = spzeros(num_rows, num_cols)
  constant = spzeros(num_rows, 1)

  if size(A, 1) != 0
    coefficient[1 : size(A, 2), 1 : size(A, 2)] = 2 * A' * A
    constant[1 : size(A, 2)] = -2 * A' * b
    if size(C, 1) != 0
      coefficient[1 : size(C, 2), size(A, 2) + 1 : end] = C'
      coefficient[size(A, 2) + 1 : end, 1 : size(C, 2)] = C
      constant[size(A, 2) + 1 : end] = -d
    end
  elseif size(C, 1) != 0
    coefficient[:, 1 : size(C, 2)] = C
    constant = -d
  else
    error("KKT system should not be empty")
  end

  return coefficient, constant
end

function populate_vars!(unique_vars_map::Dict{UInt64, AffineExpr}, vars_to_ranges_map::Dict{UInt64, Tuple{Int64, Int64}}, solution)
  for (id, var) in unique_vars_map
    var.value = spzeros(var.size...)
    var.value[:] = solution[vars_to_ranges_map[id][1] : vars_to_ranges_map[id][2]]
  end
end

function backslash_solve!(p::Problem)
  unique_vars_map = reset_values_and_get_vars!(p)
  vars_to_ranges_map, num_vars = get_var_ranges_and_num_vars(unique_vars_map)
  num_objectives = length(p.objective.affines)
  scaled_obj_affines = Array(AffineExpr, num_objectives)
  for i in 1:num_objectives
    scaled_obj_affines[i] = sqrt(p.objective.multipliers[i]) * p.objective.affines[i]
  end
  A, b = coalesce_affine_exprs(vars_to_ranges_map, num_vars, scaled_obj_affines)
  canon_forms = AffineExpr[]
  for constraint in p.constraints
    push!(canon_forms, constraint.canon_form)
  end
  C, d = coalesce_affine_exprs(vars_to_ranges_map, num_vars, canon_forms)
  if size(A, 1) == 0
    if size(C, 1) < size(C, 2)
      error("Underdetermined system of equations.")
    elseif size(C, 2) < size(C, 1)
      error("Overdetermined system of equaitons.")
    end
  end
  solution = nothing
  # return A, b, C, d
  if size(A, 1) > 0 || size(C, 1) > 0
    coefficient, constant = build_kkt_system(A, b, C, d)
    # Julia 0.4 sparse backslash uses Cholesky decomposition
    # and requires coefficient to be positive definite
    if !isposdef(coefficient)
      coefficient = full(coefficient)
    end
    try
      solution = coefficient\full(constant)
    end
  end
  if solution == nothing
    p.status = "KKT singular"
  else
    p.status = "solved"
    p.optval = p.objective.scalar
    if size(A, 1) > 0
      p.optval += sum((A*solution[1:num_vars] + b).^2)
    end
    populate_vars!(unique_vars_map, vars_to_ranges_map, solution)
  end
  return p.optval
end


text:
# build hamiltonian as a function of Bloch momentum
"""
    hamiltoniank_full!(mat::AbstractMatrix, sys::TBHamiltonian{D, T1}, k::SVector{D, T2}) where {D, T1<:Number, T2<:Number}

Returns the Bloch Hamiltonian (as a full matrix) of a periodic system, `sys` as a function of the momentum `k`.
The hoppings of the Hamiltonian must be sorted
"""
function hamiltoniank_full!(mat::AbstractMatrix, sys::TBHamiltonian{D, T1}, k::SVector{D, T2}) where {D, T1<:Number, T2<:Number}

    @assert size(mat) == (norbs(sys), norbs(sys))

    i, j = sys.hoppings[1].to, sys.hoppings[1].from
    accu = 0.0im

    for hop in sys.hoppings

        if (i, j) == (hop.to, hop.from)
            # while from and to are the same accumulate sum_n h^{i,j}_{n} e^{ikn}
            phase = cdot(k, sys.latbasis*hop.dir) + cdot(k, sys.orbitals[hop.to].pos - sys.orbitals[hop.from].pos)

            accu += hop.val*exp(-1im*phase)

        else
            # from and to have changed: store from, to and value of sum_n h^{i,j}_{n} e^{ikn}
            mat[i, j] = accu

            # reset accumulator
            phase = cdot(k, sys.latbasis*hop.dir) + cdot(k, sys.orbitals[hop.to].pos - sys.orbitals[hop.from].pos)

            accu = hop.val*exp(-1im*phase)

            # update i, j
            i, j = hop.to, hop.from

        end

    end

    # write last entry of matrix
    mat[i, j] = accu

    return mat

end

"""
    hamiltoniank_full(sys::TBHamiltonian{D, T1}, k::SVector{D, T2}) where {D, T1<:Number, T2<:Number}

Returns the Bloch Hamiltonian (as a full matrix) of a periodic system, `sys` as a function of the momentum `k`.
The hoppings of the Hamiltonian must be sorted
"""
function hamiltoniank_full(sys::TBHamiltonian{D, T1}, k::SVector{D, T2}) where {D, T1<:Number, T2<:Number}

    dim = norbs(sys)
    mat = zeros(Complex{Float64}, dim, dim)

    hamiltoniank_full!(mat, sys, k)

    return mat
end

"""
    hamiltoniank_sparse(sys::TBHamiltonian{D, T1}, k::SVector{D, T2}) where {D, T1<:Number, T2<:Number}

Returns the Bloch Hamiltonian (as a sparse matrix) of a periodic system, `sys` as a function of the momentum `k`.
The hoppings of the Hamiltonian must be sorted
"""
function hamiltoniank_sparse(sys::TBHamiltonian{D, T1}, k::SVector{D, T2}) where {D, T1<:Number, T2<:Number}

    iA = Int[]
    jA = Int[]
    V = Complex{Float64}[]

    i, j = sys.hoppings[1].to, sys.hoppings[1].from
    accu = 0.0im

    for hop in sys.hoppings

        if (i, j) == (hop.to, hop.from)
            # while from and to are the same accumulate sum_n h^{i,j}_{n} e^{ikn}
            phase = cdot(k, sys.latbasis*hop.dir) + cdot(k, sys.orbitals[hop.to].pos - sys.orbitals[hop.from].pos)

            accu += hop.val*exp(-1im*phase)

        else
            # from and to have changed: store from, to and value of sum_n h^{i,j}_{n} e^{ikn}
            push!(iA, i)
            push!(jA, j)
            push!(V, accu)

            # reset accumulator
            phase = cdot(k, sys.latbasis*hop.dir) + cdot(k, sys.orbitals[hop.to].pos - sys.orbitals[hop.from].pos)

            accu = hop.val*exp(-1im*phase)

            # update i, j
            i, j = hop.to, hop.from

        end

    end

    # write last entry of matrix
    push!(iA, i)
    push!(jA, j)
    push!(V, accu)

    return sparse(iA, jA, V)

end

"""
    hamiltoniank(sys::TBHamiltonian{D, T1}, k::SVector{D, T2}; sparse=false) where {D, T1<:Number, T2<:Number}

Returns the Bloch Hamiltonian (as a sparse or full matrix) of a periodic system, `sys` as a function of the momentum `k`.
The hoppings of the Hamiltonian must be sorted.
"""
function hamiltoniank(sys::TBHamiltonian{D, T1}, k::SVector{D, T2}; sparse=false) where {D, T1<:Number, T2<:Number}

    if sparse
        return hamiltoniank_sparse(sys, k)
    else
        return hamiltoniank_full(sys, k)
    end
end

"""
    hamiltoniank(sys::TBHamiltonian{D, T}, k::AbstractVector; sparse=false) where {D, T<:Number}

Returns the Bloch Hamiltonian (as a sparse or full matrix) of a periodic system, `sys` as a function of the momentum `k`.
The hoppings of the Hamiltonian must be sorted.
It is checked if the hoppings of the hamiltonian are sorted and if the length of the `k` vectors matches the dimensionality of the system.
"""
function hamiltoniank(sys::TBHamiltonian{D, T}, k::AbstractVector; sparse=false) where {D, T<:Number}
    @assert length(k) == D "length of k vector much be the same as the dimensionality of the system"
    @assert sys.sorted[1] "hoppings of the Hamiltonian must be sorted"

    return hamiltoniank(sys, toSVector(k), sparse=sparse)
end

"""
    hamiltoniank!(mat::AbstractMatrix, sys::TBHamiltonian{D, T}, k::AbstractVector; sparse=false) where {D, T<:Number}

In-place version of `hamiltoniank`. Only works for full matrices.
"""
function hamiltoniank!(mat::AbstractMatrix, sys::TBHamiltonian{D, T}, k::AbstractVector) where {D, T<:Number}
    @assert length(k) == D "length of k vector much be the same as the dimensionality of the system"
    @assert sys.sorted[1] "hoppings of the Hamiltonian must be sorted"

    hamiltoniank_full!(mat, sys, toSVector(k))

end


text:
"""
    solve(
        equation::DiffusionEquation{ClosedIntervalDomain{T},F},
        u,
        tlist,
        M,
        N;
        method = "filterfirst",
        boundary_conditions = "exact",
        solver = QNDF(),
        abstol = 1e-4,
        reltol = 1e-3,
        degmax = 10,
        λ = 1e-6
    ) where {T,F}

Solve `equation` from `tlist[1]` to `tlist[2]` with initial conditions `u` and a
discretization of `M` and `N` points. If `method` is `"filterfirst"`, the equation is filtered then
discretized. If `method` is `"discretizefirst"`, the equation is discretized then filtered.
"""
function solve(
    equation::DiffusionEquation{ClosedIntervalDomain{T},F},
    u,
    tlist,
    M,
    N;
    method = "filterfirst",
    boundary_conditions = "exact",
    solver = QNDF(),
    abstol = 1e-4,
    reltol = 1e-3,
    degmax = 10,
    λ = 1e-6,
) where {T,F}

    @unpack domain, filter, f, g_a, g_b = equation
    x = discretize(domain, M)
    ξ = discretize(domain, N)

    # Get matrices
    D = diffusion_matrix(domain, N)
    W = filter_matrix(filter, domain, M, N; degmax, λ)
    R = reconstruction_matrix(filter, domain, M, N; degmax, λ)

    W₀ = W[:, 1]
    Wₙ = W[:, end]

    # Initial conditions
    ū = apply_filter(u, filter, domain)
    uₕ = u.(ξ)
    # ūₕ = ū.(x)
    ūₕ = filter_matrix(filter, domain, M, N) * uₕ

    function f!(fₕ, t)
        fₕ[1] = 0
        fₕ[2:end-1] .= f.(ξ[2:end-1], t)
        fₕ[end] = 0
    end

    if boundary_conditions == "exact"
        Mass = [
            spzeros(1, N + 1)
            spzeros(N - 1, 1) sparse(I, N - 1, N - 1) spzeros(N - 1, 1)
            spzeros(1, N + 1)
        ]
        J = [
            -1 spzeros(1, N)
            D[2:end-1, :]
            spzeros(1, N) -1
        ]
        γ_a = g_a
        γ_b = g_b
    elseif boundary_conditions == "derivative"
        Mass = sparse(I, N + 1, N + 1)
        # Zero out boundary u
        J = D
        J[[1, end], :] .= 0
        dropzeros!(J)
        γ_a = t -> derivative(g_a, t)
        γ_b = t -> derivative(g_b, t)
    else
        error("Unknown boundary conditions")
    end

    if method == "filterfirst"
        error("Not implemented")
    elseif method == "discretizefirst"
        # Solve discretized-then-filtered problem
        p = (; Ju = zero(x), J = W * J * R, fₕ = zero(ξ), Wf = zero(x))
        function Mdu!(du, u, p, t)
            # @show t
            @unpack Ju, J, fₕ, Wf = p
            f!(fₕ, t)
            mul!(Wf, W, fₕ)
            mul!(Ju, J, u)
            @. du = Ju + Wf + γ_a(t) * W₀ + γ_b(t) * Wₙ
        end
        odefunction = ODEFunction(
            Mdu!,
            jac = (J, u, p, t) -> (J .= p.J),
            jac_prototype = p.J,
            mass_matrix = W * Mass * R,
        )
        problem = ODEProblem(odefunction, ūₕ, tlist, p)
        solution = OrdinaryDiffEq.solve(problem, solver; abstol, reltol)
    else
        error("Unknown method")
    end

    solution
end

"""
    solve(
        equation::DiffusionEquation{PeriodicIntervalDomain{T},F},
        u,
        tlist,
        M,
        N;
        method = "filterfirst",
        solver = QNDF(),
        abstol = 1e-4,
        reltol = 1e-3,
        degmax = 10,
        λ = 1e-6
    ) where {T,F}

Solve `equation` from `tlist[1]` to `tlist[2]` with initial conditions `u` and a
discretization of `M` and `N` points. If `method` is `"filterfirst"`, the equation is filtered then
discretized. If `method` is `"discretizefirst"`, the equation is discretized then filtered.
"""
function solve(
    equation::DiffusionEquation{PeriodicIntervalDomain{T},F},
    u,
    tlist,
    M,
    N;
    method = "filterfirst",
    solver = QNDF(),
    abstol = 1e-4,
    reltol = 1e-3,
    degmax = 10,
    λ = 1e-6,
) where {T,F}
    @unpack domain, filter = equation

    # Domain
    x = discretize(domain, M)
    ξ = discretize(domain, N)

    # Get matrices
    D = diffusion_matrix(domain, N)
    W = filter_matrix(filter, domain, M, N; degmax, λ)
    R = reconstruction_matrix(filter, domain, M, N; degmax, λ)
    A = spdiagm(α.(ξ))

    # Initial conditions
    uₕ = u.(ξ)
    ūₕ = W * uₕ

    if method == "filterfirst"
        error("Not implemented")
        # Filter
        h = filter.width
        dh(ξ) = 0.0
        α(ξ) = 1 / 3 * dh(ξ) * h(ξ)
    elseif method == "discretizefirst"
        p = (; J = W * D * R)
        du!(du, u, p, t) = mul!(du, p.J, u)
        odefunction = ODEFunction(
            du!,
            jac = (J, u, p, t) -> (J .= p.J),
            jac_prototype = p.J,
            mass_matrix = W * R,
        )
        problem = ODEProblem(odefunction, ūₕ, tlist, p)
        solution = OrdinaryDiffEq.solve(problem, solver; abstol, reltol)
    else
        error("Unknown method")
    end

    solution
end