text:
#########################################################
#                                                       #
#  marginal effects of exogenous determinants on E(u)   #
#                                                       #
#########################################################

#? ----------- truncated normal, marginal effect function -------

function marg_trun( # PorC::Int64,
                    pos::NamedTuple,
                    coef::Array{Float64, 1},
                    Zmarg, Wmarg)

         #* No need `modelid` because each has to be customized (cannot accept
         #*   qvar())

         #* Zmarg is the marginal effect, which acts like "coeff" in a usual
         #*   MLE model.
         #*   (1) This is obs-by-obs, so for example,
         #*       use X'*coef[...] rather than X*coef[....].
         #*   (2) Marg eff is on E(u), not BC or JLMS.
         #*   (3) The line `return ....` cannot be used.

       z_pre = Zmarg'*coef[pos.begz : pos.endz]  # mu
       w_pre = Wmarg'*coef[pos.begw : pos.endw] # log_σᵤ²

       μ   = z_pre
       σᵤ  = exp(0.5 * w_pre)

       Λ = μ/σᵤ

       uncondU = σᵤ* (Λ + normpdf(Λ) / normcdf(Λ)) # kx1
end

#? ----------- truncated normal, get marginal effect -------

function get_marg(::Type{Trun}, # PorC::Int64,
                  pos::NamedTuple, num::NamedTuple, coef::Array{Float64, 1},
                  Z::Matrix, q, W::Matrix)

    mm_z = Array{Float64}(undef, num.nofz, num.nofobs)
    mm_w = Array{Float64}(undef, num.nofw, num.nofobs)

    @inbounds for i in 1:num.nofobs
        @views marg = ForwardDiff.gradient(marg -> marg_trun(pos, coef,
                                                      marg[1 : num.nofz],
                                                      marg[num.nofz+1 : num.nofmarg]),
                                    vcat(Z[i,:], W[i,:]) );

        mm_z[:,i] = marg[1 : num.nofz]
        mm_w[:,i] = marg[num.nofz+1 : end]

     end  # for i in 1:num.nofobs

     margeff = DataFrame(mm_z', _dicM[:μ]) # the base set
        mm_w = DataFrame(mm_w', _dicM[:σᵤ²])

      #* purge off the constant var's marginal effect from the DataFrame
         margeff = nonConsDataFrame(margeff, Z)
            mm_w = nonConsDataFrame(mm_w, W)

      #* if same var in different equations, add up the marg eff
         margeff = addDataFrame(margeff, mm_w)

      #* prepare info for printing
         margMean = (; zip(Symbol.(names(margeff)) , round.(mean.(eachcol(margeff)); digits=5))...)

      #* modify variable names to indicate marginal effects
         newname = Symbol.(fill("marg_", (size(margeff,2), 1)) .* names(margeff))
         margeff = rename!(margeff, vec(newname))

     return  margeff, margMean

end

#? ---------------------------------------------------------
#? ----------- half normal, marginal effect function -------
#? ---------------------------------------------------------

function marg_half( # PorC::Int64,
     pos::NamedTuple, coef::Array{Float64, 1},
     Wmarg)

     w_pre = Wmarg'*coef[pos.begw : pos.endw] # log_σᵤ²
     σᵤ  = exp(0.5 * w_pre)
     uncondU = sqrt(2/π) * σᵤ # kx1

end

#? ----------- half normal, get marginal effect -------

function get_marg(::Type{Half}, # PorC::Int64,
     pos::NamedTuple, num::NamedTuple, coef::Array{Float64, 1},
     z, q, W::Matrix)

     mm_w = Array{Float64}(undef, num.nofw, num.nofobs)

     @inbounds for i in 1:num.nofobs
     @views marg = ForwardDiff.gradient(marg -> marg_half(pos, coef,
                                             marg[1 : num.nofw]),
                         vcat(W[i,:]) );

     mm_w[:,i] = marg[1 : end]

     end  # for i in 1:num.nofobs

     margeff = DataFrame(mm_w', _dicM[:σᵤ²]) # the base set

     #* purge off the constant var's marginal effect from the DataFrame
         margeff = nonConsDataFrame(margeff, W)

     #* if same var in different equations, add up the marg eff
         #margeff = addDataFrame(margeff, mm_w)

      #* prepare info for printing
         margMean = (; zip(Symbol.(names(margeff)) , round.(mean.(eachcol(margeff)); digits=5))...)

      #* modify variable names to indicate marginal effects
         newname = Symbol.(fill("marg_", (size(margeff,2), 1)) .* names(margeff))
         margeff = rename!(margeff, vec(newname))

     return  margeff, margMean
end

#? ---------------------------------------------------------
#? ---------- Exponential normal, marginal effect function -
#? ---------------------------------------------------------

function marg_expo( pos::NamedTuple, coef::Array{Float64, 1},
     Wmarg)

     w_pre = Wmarg'*coef[pos.begw : pos.endw] # log_σᵤ²
     λ   = exp(0.5*w_pre)
     uncondU = λ # kx1

end

#? ----------- Exponential, get marginal effect -------

function get_marg(::Type{Expo}, # PorC::Int64,
     pos::NamedTuple, num::NamedTuple, coef::Array{Float64, 1},
     z, q, W::Matrix)

     mm_w = Array{Float64}(undef, num.nofw, num.nofobs)

     @inbounds for i in 1:num.nofobs
     @views       marg = ForwardDiff.gradient(marg -> marg_expo( pos, coef,
                                           marg[1 : num.nofw]),
                                 vcat(W[i,:]) );

     mm_w[:,i] = marg[1 : end]

     end  # for i in 1:num.nofobs

     margeff = DataFrame(mm_w', _dicM[:σᵤ²]) # the base set

     #* purge off the constant var's marginal effect from the DataFrame
         margeff = nonConsDataFrame(margeff, W)

     #* if same var in different equations, add up the marg eff
         #margeff = addDataFrame(margeff, mm_w)

      #* prepare info for printing
         margMean = (; zip(Symbol.(names(margeff)) , round.(mean.(eachcol(margeff)); digits=5))...)

      #* modify variable names to indicate marginal effects
         newname = Symbol.(fill("marg_", (size(margeff,2), 1)) .* names(margeff))
         margeff = rename!(margeff, vec(newname))

     return  margeff, margMean
end

#? ---------------------------------------------------------
#? ------ Scaling Property Model, marginal effect function -
#? ---------------------------------------------------------

function marg_scal( # PorC::Int64,
     pos::NamedTuple, coef::Array{Float64, 1},
     Zmarg, Qmarg, Wmarg)

     z_pre = Zmarg'*coef[pos.begz : pos.endz]  # mu
     q_pre = Qmarg'*coef[pos.begq : pos.endq]
     w_pre = Wmarg'*coef[pos.begw : pos.endw] # log_σᵤ²

     μ   = (z_pre) * exp(q_pre)             # Q is _cons; μ here is the after-mutiplied-by-scaling function
     σᵤ  = exp(0.5 * w_pre + q_pre)
     Λ   = μ/σᵤ

     uncondU = σᵤ* (Λ + normpdf(Λ) / normcdf(Λ))

end

#? ----------- scaling property model, get marginal effect -------

function get_marg(::Type{Trun_Scale}, # PorC::Int64,
     pos::NamedTuple, num::NamedTuple, coef::Array{Float64, 1},
     Z::Matrix, Q::Matrix, W::Matrix)

     mm_z = Array{Float64}(undef, num.nofz, num.nofobs)
     mm_q = Array{Float64}(undef, num.nofq, num.nofobs)
     mm_w = Array{Float64}(undef, num.nofw, num.nofobs)

     @inbounds for i in 1:num.nofobs
     @views marg = ForwardDiff.gradient(marg -> marg_scal( pos, coef,
                                             marg[1 : num.nofz],
                                             marg[num.nofz+1 : num.nofz+num.nofq],
                                             marg[num.nofz+num.nofq+1 : num.nofmarg]),
                         vcat(Z[i,:], Q[i,:], W[i,:]) );

     mm_z[:,i] = marg[1 : num.nofz]
     mm_q[:,i] = marg[num.nofz+1 : num.nofz+num.nofq]
     mm_w[:,i] = marg[num.nofz+num.nofq+1 : end]

     end  # for i in 1:num.nofobs

     margeff = DataFrame(mm_z', _dicM[:μ]) # the base set
     mm_q = DataFrame(mm_q', _dicM[:hscale])
     mm_w = DataFrame(mm_w', _dicM[:σᵤ²])

     #* purge off the constant var's marginal effect from the DataFrame
     margeff = nonConsDataFrame(margeff, Z)
     mm_q = nonConsDataFrame(mm_q, Q)
     mm_w = nonConsDataFrame(mm_w, W)

     #* if same var in different equations, add up the marg eff
     margeff = addDataFrame(margeff, mm_q)
     margeff = addDataFrame(margeff, mm_w)

     #* prepare info for printing
     margMean = (; zip(Symbol.(names(margeff)) , round.(mean.(eachcol(margeff)); digits=5))...)

     #* modify variable names to indicate marginal effects
     newname = Symbol.(fill("marg_", (size(margeff,2), 1)) .* names(margeff))
     margeff = rename!(margeff, vec(newname))

     return  margeff, margMean

end

#? ---------------------------------------------------------
#? - panel FE Wang and Ho, Half normal, marginal effect function -
#? ---------------------------------------------------------

function marg_fewhh( # PorC::Int64,
     pos::NamedTuple, coef::Array{Float64, 1},
     Qmarg, Wmarg)

     q_pre = Qmarg'*coef[pos.begq : pos.endq]
     w_pre = Wmarg'*coef[pos.begw : pos.endw] # log_σᵤ²

     σᵤ  = exp(0.5*w_pre + q_pre)
     uncondU = sqrt(2/π) * σᵤ
end

#? -- panel FE Wang and Ho, Half normal, , get marginal effect ----

function get_marg(::Type{PFEWHH}, # PorC::Int64,
     pos::NamedTuple, num::NamedTuple, coef::Array{Float64, 1},
     z, Q::Matrix, W::Matrix)

        #* Note that Y and X are within-transformed by `getvar`,
        #* but Q, W, V are still in the original level.

     mm_q = Array{Float64}(undef, num.nofq, num.nofobs)
     mm_w = Array{Float64}(undef, num.nofw, num.nofobs)

     @inbounds for i in 1:num.nofobs
          @views marg = ForwardDiff.gradient(marg -> marg_fewhh( pos, coef,
                                             marg[1 : num.nofq],
                                             marg[num.nofq+1 : num.nofmarg]),
                         vcat( Q[i,:], W[i,:]) );

     mm_q[:,i] = marg[1 : num.nofq]
     mm_w[:,i] = marg[num.nofq+1 : end]

     end  # for i in 1:num.nofobs

     margeff = DataFrame(mm_q', _dicM[:hscale])
     mm_w    = DataFrame(mm_w', _dicM[:σᵤ²])

     #* purge off the constant var's marginal effect from the DataFrame
     margeff = nonConsDataFrame(margeff, Q)
     mm_w    = nonConsDataFrame(mm_w, W)

     #* if same var in different equations, add up the marg eff
     margeff = addDataFrame(margeff, mm_w)

     #* prepare info for printing
     margMean = (; zip(Symbol.(names(margeff)) , round.(mean.(eachcol(margeff)); digits=5))...)

     #* modify variable names to indicate marginal effects
     newname = Symbol.(fill("marg_", (size(margeff,2), 1)) .* names(margeff))
     margeff = rename!(margeff, vec(newname))

     return  margeff, margMean
end

#? --------------------------------------------------------------------
#? - panel FE Wang and Ho, truncated normal, marginal effect function -
#? --------------------------------------------------------------------

function marg_fewht( # PorC::Int64,
     pos::NamedTuple, coef::Array{Float64, 1},
     Zmarg, Qmarg, Wmarg)

       z_pre = Zmarg'*coef[pos.begz : pos.endz]  # mu, a scalar
       q_pre = Qmarg'*coef[pos.begq : pos.endq]
       w_pre = Wmarg'*coef[pos.begw : pos.endw] # log_σᵤ²

      μ  = exp(q_pre)*z_pre
      σᵤ = exp(q_pre + 0.5*w_pre)
      Λ  = μ/σᵤ

    uncondU = σᵤ* (Λ + normpdf(Λ) / normcdf(Λ)) # kx1

end

#? -- panel FE Wang and Ho, truncated normal, , get marginal effect ----

function get_marg(::Type{PFEWHT}, # PorC::Int64,
     pos::NamedTuple, num::NamedTuple, coef::Array{Float64, 1},
     Z::Matrix, Q::Matrix, W::Matrix)

        #* Note that Y and X are within-transformed by `getvar`,
        #* but Q, W, V are still in the original level.

     mm_z = Array{Float64}(undef, num.nofz, num.nofobs)
     mm_q = Array{Float64}(undef, num.nofq, num.nofobs)
     mm_w = Array{Float64}(undef, num.nofw, num.nofobs)

     @inbounds for i in 1:num.nofobs
     @views marg = ForwardDiff.gradient(marg -> marg_fewht( pos, coef,
                                             marg[1 : num.nofz],
                                             marg[num.nofz+1 : num.nofz+num.nofq],
                                             marg[num.nofz+num.nofq+1 : num.nofmarg]),
                         vcat( Z[i,:], Q[i,:], W[i,:]) );

     mm_z[:,i] = marg[1 : num.nofz]
     mm_q[:,i] = marg[num.nofz+1 : num.nofz+num.nofq]
     mm_w[:,i] = marg[num.nofz+num.nofq+1 : end]

     end  # for i in 1:num.nofobs

     margeff = DataFrame(mm_z', _dicM[:μ]) # the base set
     mm_q = DataFrame(mm_q', _dicM[:hscale])
     mm_w = DataFrame(mm_w', _dicM[:σᵤ²])

     #* purge off the constant var's marginal effect from the DataFrame
     margeff = nonConsDataFrame(margeff, Z)
     mm_q = nonConsDataFrame(mm_q, Q)
     mm_w = nonConsDataFrame(mm_w, W)

     #* if same var in different equations, add up the marg eff
     margeff = addDataFrame(margeff, mm_q)
     margeff = addDataFrame(margeff, mm_w)

      #* prepare info for printing
      margMean = (; zip(Symbol.(names(margeff)) , round.(mean.(eachcol(margeff)); digits=5))...)

      #* modify variable names to indicate marginal effects
      newname = Symbol.(fill("marg_", (size(margeff,2), 1)) .* names(margeff))
      margeff = rename!(margeff, vec(newname))

     return  margeff, margMean
end

#? --------------------------------------------------------------------
#? - Panel Time Decay Model, marginal effect function -
#? --------------------------------------------------------------------

function marg_pdecay( # PorC::Int64,
     pos::NamedTuple, coef::Array{Float64, 1},
     Zmarg, Qmarg, Wmarg)

     #* same as marg_fewht

     z_pre = Zmarg'*coef[pos.begz : pos.endz]  # mu, a scalar
     q_pre = Qmarg'*coef[pos.begq : pos.endq]
     w_pre = Wmarg'*coef[pos.begw : pos.endw] # log_σᵤ²

     μ  = exp(q_pre)*z_pre
     σᵤ = exp(q_pre + 0.5*w_pre)
     Λ  = μ/σᵤ

     uncondU = σᵤ* (Λ + normpdf(Λ) / normcdf(Λ)) # kx1

end

#? -- panel time decay model, get marginal effect ----

function get_marg(::Type{PanDecay}, # PorC::Int64,
     pos::NamedTuple, num::NamedTuple, coef::Array{Float64, 1},
     Z::Matrix, Q::Matrix, W::Matrix)

        #* same as fewht

     mm_z = Array{Float64}(undef, num.nofz, num.nofobs)
     mm_q = Array{Float64}(undef, num.nofq, num.nofobs)
     mm_w = Array{Float64}(undef, num.nofw, num.nofobs)

     @inbounds for i in 1:num.nofobs
     @views marg = ForwardDiff.gradient(marg -> marg_pdecay( pos, coef,
                                             marg[1 : num.nofz],
                                             marg[num.nofz+1 : num.nofz+num.nofq],
                                             marg[num.nofz+num.nofq+1 : num.nofmarg]),
                         vcat( Z[i,:], Q[i,:], W[i,:]) );

     mm_z[:,i] = marg[1 : num.nofz]
     mm_q[:,i] = marg[num.nofz+1 : num.nofz+num.nofq]
     mm_w[:,i] = marg[num.nofz+num.nofq+1 : end]

     end  # for i in 1:num.nofobs

     margeff = DataFrame(mm_z', _dicM[:μ]) # the base set
     mm_q = DataFrame(mm_q', _dicM[:gamma])
     mm_w = DataFrame(mm_w', _dicM[:σᵤ²])

     #* purge off the constant var's marginal effect from the DataFrame
     margeff = nonConsDataFrame(margeff, Z)
     mm_q = nonConsDataFrame(mm_q, Q)
     mm_w = nonConsDataFrame(mm_w, W)

     #* if same var in different equations, add up the marg eff
     margeff = addDataFrame(margeff, mm_q)
     margeff = addDataFrame(margeff, mm_w)

     #* prepare info for printing
     margMean = (; zip(Symbol.(names(margeff)) , round.(mean.(eachcol(margeff)); digits=5))...)

     #* modify var name to indicate marginal effect
     newname = Symbol.(fill("marg_", (size(margeff,2), 1)) .* names(margeff))
     margeff = rename!(margeff, vec(newname))

     return  margeff, margMean

end

#? --------------------------------------------------------------------
#? - Panel Kumbhakar 1990 Model, marginal effect function -
#? --------------------------------------------------------------------

function marg_pkumb90( # PorC::Int64,
     pos::NamedTuple, coef::Array{Float64, 1},
     Zmarg, Qmarg, Wmarg)

     #* same as marg_fewht

     z_pre = Zmarg'*coef[pos.begz : pos.endz]  # mu, a scalar
     q_pre = Qmarg'*coef[pos.begq : pos.endq]
     w_pre = Wmarg'*coef[pos.begw : pos.endw] # log_σᵤ²

     μ  = (2 ./(1+exp(q_pre)))*z_pre # exp(q_pre)*z_pre
     σᵤ = (2 ./(1+exp(q_pre)))*exp(0.5*w_pre)    # exp(q_pre + 0.5*w_pre)
     Λ  = μ/σᵤ

     uncondU = σᵤ* (Λ + normpdf(Λ) / normcdf(Λ)) # kx1

end

#? -- panel Kumbhakar 1990 model, get marginal effect ----

function get_marg(::Type{PanKumb90}, # PorC::Int64,
     pos::NamedTuple, num::NamedTuple, coef::Array{Float64, 1},
     Z::Matrix, Q::Matrix, W::Matrix)

        #* same as fewht

     mm_z = Array{Float64}(undef, num.nofz, num.nofobs)
     mm_q = Array{Float64}(undef, num.nofq, num.nofobs)
     mm_w = Array{Float64}(undef, num.nofw, num.nofobs)

     @inbounds for i in 1:num.nofobs
     @views marg = ForwardDiff.gradient(marg -> marg_pkumb90( pos, coef,
                                             marg[1 : num.nofz],
                                             marg[num.nofz+1 : num.nofz+num.nofq],
                                             marg[num.nofz+num.nofq+1 : num.nofmarg]),
                         vcat( Z[i,:], Q[i,:], W[i,:]) );

     mm_z[:,i] = marg[1 : num.nofz]
     mm_q[:,i] = marg[num.nofz+1 : num.nofz+num.nofq]
     mm_w[:,i] = marg[num.nofz+num.nofq+1 : end]

     end  # for i in 1:num.nofobs

     margeff = DataFrame(mm_z', _dicM[:μ]) # the base set
     mm_q = DataFrame(mm_q', _dicM[:gamma])
     mm_w = DataFrame(mm_w', _dicM[:σᵤ²])

     #* purge off the constant var's marginal effect from the DataFrame
     margeff = nonConsDataFrame(margeff, Z)
     mm_q = nonConsDataFrame(mm_q, Q)
     mm_w = nonConsDataFrame(mm_w, W)

     #* if same var in different equations, add up the marg eff
     margeff = addDataFrame(margeff, mm_q)
     margeff = addDataFrame(margeff, mm_w)

     #* prepare info for printing
     margMean = (; zip(Symbol.(names(margeff)) , round.(mean.(eachcol(margeff)); digits=5))...)

     #* modify var name to indicate marginal effect
     newname = Symbol.(fill("marg_", (size(margeff,2), 1)) .* names(margeff))
     margeff = rename!(margeff, vec(newname))

     return  margeff, margMean

end

#? -----------------------------------------------------------------------------
#? - panel FE CSW (JoE 2014), CSN model; Half normal, marginal effect function -
#? -----------------------------------------------------------------------------

#? -- panel FE CSW (JoE 2014), CSN model; Half normal, get marginal effect ----

function get_marg(::Type{PFECSWH},
     pos::NamedTuple, num::NamedTuple, coef::Array{Float64, 1},
     z, q, w)

     # This model does not allow exogenous determinants of inefficiency
     # thus no marginal effect.

     return  (), ()
end

#? -----------------------------------------------------------------------------
#? - panel true random effect model, half normal and truncated normal
#? -----------------------------------------------------------------------------

#? --  Half normal ---------------

function get_marg(::Type{PTREH},
     pos::NamedTuple, num::NamedTuple, coef::Array{Float64, 1},
     z, q, w)

     # This model does not allow exogenous determinants of inefficiency
     # thus no marginal effect.

     return  (), ()
end

#? --  truncated normal ---------------

function get_marg(::Type{PTRET},
     pos::NamedTuple, num::NamedTuple, coef::Array{Float64, 1},
     z, q, w)

     # This model does not allow exogenous determinants of inefficiency
     # thus no marginal effect.

     return  (), ()
end


text:
function generateNodeFluxPlot( fluxReport::Dict{String,DataFrame},
    nodes::NTuple{N,String}, fluxType::Symbol, timeFactor::Float64,
    showPlot::Bool, savePlot::Bool, filename::String,
    plotFunction::Function ) where N

    fileroot = nothing
    extension = nothing

    # Get file name ready.
    if savePlot
        # If the user doesn't specify a filename, a random one is generated.
        if filename == ""
            dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
            fileroot = string( string( "Node flux plot (", dateStr, ")",
                labelInfix[:fluxType] ) )
            extension = "svg"
        else
            extension = split( filename, "." )
            fileroot = length( extension ) == 1 ? extension[1] :
                join( extension[1:(end - 1)], "." )
            extension = length( extension ) == 1 ? "" : extension[end]
            extension = extension ∈ extensions ? extension : "svg"
        end  # if filename == ""

        fileroot = string( fileroot )

        if !ispath( dirname( fileroot ) )
            mkpath( dirname( fileroot ) )
        end  # if !ispath( dirname( fileroot ) )
    end  # if savePlot

    # Generate plots.
    for node in unique( nodes )
        if !haskey( fluxReport, node )
            continue
        end  # if !haskey( fluxReport, node )

        plotfilename = string( fileroot, " ", node, ".", extension )
        plt = plotFunction( fluxReport[node], string( node ), fluxType,
            timeFactor )

        # Show plot if needed.
        if showPlot
            gui( plt )
        end  # if showPlot

        # Save plot if needed.
        if savePlot
            savefig( plt, plotfilename )
        end  # if savePlot
    end  # for node in unique( nodes )

end  # generateNodeFluxPlot( fluxReport, nodes, fluxType, timeFactor, showPlot,
     #   savePlot, fileroot, extension, plotFunction )

function normalFluxPlot( fluxReport::DataFrame, node::String, fluxType::Symbol,
    timeFactor::Float64 )

    timePoints = fluxReport[:, :timePoint]
    labels = generateFluxLabels( names( fluxReport )[4:end], fluxType )
    labels = vcat( string( "Flux", labelInfix[fluxType][1], node ),
        labels )
    plotTitle = string( "Breakdown of ", labels[1] )
    timePoints = fluxReport[:, :timePoint] / timeFactor
    dataPoints = Matrix( fluxReport[:, 3:end] )

    plt = nothing

    if fluxType === :within
        ymax = size( dataPoints, 2 ) == 1 ? 20 :
            max( maximum( dataPoints[:, 2:end] ), 20 )
        plt = plot( show=false, size=(960, 540), title=plotTitle,
            label=labels[1], lw=3, color=:black, ylim=[0, ymax] )
    else
        ymax = max( maximum( dataPoints[:, 1] ), 20 )
        plt = plot( timePoints, dataPoints[:, 1], show=false, size=(960, 540),
            title=plotTitle, label=labels[1], lw=3, color=:black,
            ylim=[0, ymax] )
    end  # if fluxType === :within

    plot!( plt, timePoints, dataPoints[:, 2:end],
        label=hcat( labels[2:end]... ), lw=2 )

    return plt

end  # normalFluxPlot( fluxReport, node, fluxType, timeFactor )

function stackedFluxPlot( fluxReport::DataFrame, node::String, fluxType::Symbol,
    timeFactor::Float64 )

    timePoints = fluxReport[:, :timePoint]
    labels = generateFluxLabels( names( fluxReport )[4:end], fluxType )
    plotTitle = string( "Breakdown of ", labels[1] )
    timePoints = fluxReport[:, :timePoint] / timeFactor
    dataPoints = Matrix( fluxReport[:, 4:end] )
    nPoints, nSeries = size( dataPoints )
    dataPoints = hcat( fill( 0, nPoints ), cumsum( dataPoints, dims=2 ) )
    ymax = max( maximum( fluxReport[:, 3] ), 20 )

    plt = plot( show=false, size=(960,540), title=plotTitle, ylim=[0, ymax] )

    for ii in 1:nSeries
        plot!( plt, timePoints, dataPoints[:, ii + 1], lw=2, lalpha=1.0, fillto=dataPoints[:, ii], falpha=0.5, label=labels[ii] )
    end  # for ii in 1:nSeries

    return plt

end  # stackedFluxPlot( fluxReport, node, fluxType, timeFactor )

function generateFluxLabels( labels::Vector{String}, fluxType::Symbol )

    # Split the labels into transition, source node, and target node.
    labelBreakdown = split.( labels, ": " )
    transitions = getindex.( labelBreakdown, 1 )
    labelBreakdown = split.( getindex.( labelBreakdown, 2 ), " => " )
    sourceNodes = getindex.( labelBreakdown, 1 )
    targetNodes = getindex.( labelBreakdown, 2 )

    if fluxType === :within
        return string.( sourceNodes, labelInfix[fluxType][2], targetNodes )
    end  # if fluxType === :within

    return string.( transitions, labelInfix[fluxType][2],
        fluxType === :in ? sourceNodes : targetNodes )

end  # generateFluxLabels( labels, plotType )

ext: jl


text:
using TuringModels, StatsFuns

Turing.setadbackend(:reverse_diff);
#Turing.turnprogress(false);

d = DataFrame(CSV.read(joinpath(@__DIR__, "..", "..", "data", "chimpanzees.csv"),
  delim=';'));
size(d) # Should be 504x8

# pulled_left, condition, prosoc_left
@model m10_3(y, x₁, x₂) = begin
    α ~ Normal(0, 10)
    βp ~ Normal(0, 10)
    βpC ~ Normal(0, 10)

    for i ∈ 1:length(y)
        p = logistic(α + (βp + βpC * x₁[i]) * x₂[i])
        y[i] ~ Binomial(1, p)
    end
end;

chns = sample(m10_3(d[:,:pulled_left], d[:,:condition], d[:,:prosoc_left]),
  Turing.NUTS(0.95), 2000);

# Rethinking result

m_10_03t_result = "
      Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
 a    0.05   0.13      -0.15       0.25  3284    1
 bp   0.62   0.22       0.28       0.98  3032    1
 bpC -0.11   0.26      -0.53       0.29  3184    1
";

# Describe the draws

describe(chns)

# End of m10.03t.jl


text:
using Distributions
using DifferentialEquations
using DataFrames
using CSV
using PyPlot

include("SICRmodelsnew.jl")
include("SICUmodels.jl")

"""
	function mcmc(data,cols,p,N,total::Int,sicrfn!)

Metropolis-Hastings MCMC of the SIR ODE model assuming a Poisson likelihood
for the cumulative number of cases and deaths, and a noninformative prior
data = DataFrame of
cols = [active cases, deaths, recovered cases]
p[1:end-1] = array of parameters used in the model
p[end] = number infected on day 0 of data (Jan 22, 2020)
N = population of data region
total = number of MCMC steps
output = (loglikelihoods, parameters, maximum loglikelihood, max loglikelihood parameters)
sicrfn! = ODE file to use

"""
function mcmc(data::DataFrame,cols::Array,p::Array,N::Float64,total::Int,sicrfn!,temp::Float64 = 1.)

	chi = ll(data,cols,p,N,sicrfn!)
	accept = 0
	pt = copy(p)
	chiml = chi
	pml = copy(p)

	chiout = Array{Float64,1}(undef,total)
	pout = Array{Float64,2}(undef,total,length(p))

	for step in 1:total
		sampleparams!(pt,p)
		chit = ll(data,cols,pt,N,sicrfn!)
		if rand() < exp((chi - chit)/temp)
			chi = chit
			p = copy(pt)
			accept += 1
			if chi < chiml
				chiml = chi
				pml = copy(p)
			end
		end
		chiout[step] = chi
		pout[step,:] = p
	end
	return chiout, pout, accept/total, chiml, pml
end

# generate MCMC step guess from current parameter value r using logNormal distirbution for all
# parameters except start day with scale parameter manually adjusted so approximately a third of guesses are accepted
function sampleparams!(pt,p)
	for i in eachindex(p)
		d = Distributions.LogNormal(log(p[i]),.002)
		pt[i] = rand(d)
	end
end

# loglikelihood of data given parameters at r using a Poisson likelihood function
function ll(data,cols,p,N,sicrfn!)
	firstday = findfirst(data[:,cols[1]] .>0)
	lastday = length(data[:,cols[1]])
	tspan = (firstday,lastday)
	uin = uinitial(data,cols,firstday,N)
	prediction,days = model(p,uin,float.(tspan),N,sicrfn!)
	chi = 0
	for set in eachindex(cols)
		for i in eachindex(days)
			d = Distributions.Poisson(max(prediction[set,i],2.))
			chi -= loglikelihood(d,[data[days[i],cols[set]]])
		end
	end
	return chi
end

function uinitial(data,cols,firstday,N)
	uin = Array{Float64,1}(undef,length(cols))
	for i in eachindex(cols)
		uin[i] = data[firstday,cols[i]]/N
	end
	return uin
end

# produce arrays of cases and deaths
function model(pin,uin,tspan,N,sicrfn!)
	sol = modelsol(pin,uin,tspan,N,sicrfn!)
	return sol[1:3,:]*N, round.(Int,sol.t)
end

function modelsol(pin,uin,tspan,N,sicrfn!)
	p = pin[1:end-1]
	I0 = (p[1]- p[3])/p[2] * uin[1] + pin[end]/N
	Z0 = p[1]/p[2]*uin[1]
	u0 = vcat(uin,[I0,Z0])
	prob = ODEProblem(sicrfn!,u0,tspan,p)
	solve(prob,saveat=1.)
end

function makehistograms(out)
	for i in 1:size(out)[2]
		subplot(2,3,i)
		hist(out[:,i],bins=100,density=true)
		title("param$(i)",fontsize=15)
	end
	tight_layout()
end

# returns predictions for C,R,D, matching data, and array of days
# solution array is in same row col order as data
# data columns have been ordered to match solutions.
function modelprediction(data,cols,p,N,sicrfn!)
	firstday = findfirst(data[:,cols[1]] .>0)
	lastday = length(data[:,cols[1]])
	sol = modelprediction(data,cols,p,N,sicrfn!,lastday)
	return sol[1:length(cols),:]'*N, data[firstday:end,cols], Int.(sol.t)
end

# returns solutions to model from day of first case to requested last day
function modelprediction(data,cols,p,N,sicrfn!,lastday)
	firstday = findfirst(data[:,cols[1]] .>0)
	tspan = (firstday,lastday)
	uin = uinitial(data,cols,firstday,N)
	modelsol(p,float.(uin),float.(tspan),N,sicrfn!)
end

# returns statistics for posteriors and other measures
function measures(psamples,pml)
	nparams = size(psamples,2)
	pout = Array{Tuple,1}(undef,nparams)
	for i in 1:nparams
		p = psamples[:,i]
		pout[i] = pml[i],mean(p),std(p),quantile(p,[.025,.5,.975])
	end
	R0 = pml[1]/(pml[2]+pml[3])
	return pout,R0
end

max_stars_repo_name: cannin/covid-sicr-test

max_issues_repo_name: cannin/covid-sicr-test

max_forks_repo_name: cannin/covid-sicr-test


text:
mutable struct MuffinTinRadialVars
    # scale factor for number of muffin-tin points
    nrmtscf::Float64
    # number of muffin-tin radial points for each species
    nrmt::Vector{Int64}
    # maximum nrmt over all the species
    nrmtmax::Int64
    # optional default muffin-tin radius for all atoms
    rmtall::Float64
    # minimum allowed distance between muffin-tin surfaces
    rmtdelta::Float64
    # muffin-tin radii
    rmt::Vector{Float64}
    # total muffin-tin volume
    omegamt::Float64
    # radial step length for coarse mesh
    lradstp::Int64
    # number of coarse radial mesh points
    nrcmt::Vector{Int64}
    # maximum nrcmt over all the species
    nrcmtmax::Int64
    # coarse muffin-tin radial mesh
    rcmt::Array{Float64,2}
    # r^l on fine radial mesh
    rlmt::Array{Float64,3}
    # r^l on coarse radial mesh
    rlcmt::Array{Float64,3}
    # weights for spline integration on fine radial mesh
    wrmt::Array{Float64,2}
    # weights for spline partial integration on fine radial mesh
    wprmt::Array{Float64,3}
    # weights for spline integration on coarse radial mesh
    wrcmt::Array{Float64,2}
    # weights for spline partial integration on coarse radial mesh
    wprcmt::Array{Float64,3}
    # maximum allowable angular momentum for augmented plane waves
    maxlapw::Int64 # parameter =50
    # maximum angular momentum for augmented plane waves
    lmaxapw::Int64
    # (lmaxapw+1)^2
    lmmaxapw::Int64
    # maximum angular momentum on the outer part of the muffin-tin
    lmaxo::Int64
    # (lmaxo+1)^2
    lmmaxo::Int64
    # maximum angular momentum on the inner part of the muffin-tin
    lmaxi::Int64
    # (lmaxi+1)^2
    lmmaxi::Int64
    # fraction of muffin-tin radius which constitutes the inner part
    fracinr::Float64
    # number of fine/coarse radial points on the inner part of the muffin-tin
    nrmti::Vector{Int64}
    nrcmti::Vector{Int64}
    # index to (l,m) pairs
    idxlm::Array{Int64,2}
    # inverse index to (l,m) pairs
    idxil::Vector{Int64}
    idxim::Vector{Int64}
    # number of fine/coarse points in packed muffin-tins
    npmti::Vector{Int64}
    npmt::Vector{Int64}
    npcmti::Vector{Int64}
    npcmt::Vector{Int64}
    # maximum number of points over all packed muffin-tins
    npmtmax::Int64
    npcmtmax::Int64
end

function MuffinTinRadialVars(Nspecies)

    nrmtscf = 0
    nrmt = zeros(Int64,Nspecies)
    nrmtmax = 0
    rmtall = 0
    rmtdelta = 0.05 # default
    rmt = zeros(Nspecies)
    omegamt = 0.0
    lradstp = 4
    nrcmt = zeros(Int64,Nspecies)
    nrcmtmax = 0
    rcmt = zeros(Float64,1,1)
    rlmt = zeros(Float64,1,1,1)
    rlcmt = zeros(Float64,1,1,1)
    wrmt = zeros(Float64,1,1)
    wprmt = zeros(Float64,1,1,1)
    wrcmt = zeros(Float64,1,1)
    wprcmt = zeros(Float64,1,1,1)
    maxlapw = 50
    lmaxapw = 8 # default
    lmmaxapw = (lmaxapw+1)^2
    lmaxo  = 6
    lmmaxo = (lmaxo+1)^2
    lmaxi  = 1 # min(lmaxi,lmaxo)
    lmmaxi = (lmaxi+1)^2
    fracinr = 0.01

    nrmti = zeros(Int64,Nspecies)
    nrcmti = zeros(Int64,Nspecies)

    idxlm = zeros(Int64,1,1)
    idxil = zeros(Int64,1)
    idxim = zeros(Int64,1)

    npmti = zeros(Int64,Nspecies)
    npmt = zeros(Int64,Nspecies)
    npcmti = zeros(Int64,Nspecies)
    npcmt = zeros(Int64,Nspecies)

    npmtmax = 0
    npcmtmax = 0

    return MuffinTinRadialVars(
        nrmtscf, nrmt, nrmtmax, rmtall, rmtdelta, rmt, omegamt, lradstp,
        nrcmt, nrcmtmax, rcmt, rlmt, rlcmt, wrmt, wprmt, wrcmt, wprcmt,
        maxlapw, lmaxapw, lmmaxapw, lmaxo, lmmaxo, lmaxi, lmmaxi, fracinr,
        nrmti, nrcmti, idxlm, idxil, idxim, npmti, npmt, npcmti, npcmt,
        npmtmax, npcmtmax,
    )

end

function init_zero!( mtr_vars::MuffinTinRadialVars )

    nrmt = mtr_vars.nrmt
    nspecies = size(nrmt)[1]

    nrcmt = mtr_vars.nrcmt
    lradstp = mtr_vars.lradstp

    # make the muffin-tin mesh commensurate with lradstp
    for is in 1:nspecies
        nrmt[is] = nrmt[is] - (nrmt[is]-1)%lradstp
        nrcmt[is] =( nrmt[is] - 1)/lradstp + 1
    end
    mtr_vars.nrmtmax = maximum(nrmt)
    mtr_vars.nrcmtmax = maximum(nrcmt)

    @assert mtr_vars.lmaxo <= mtr_vars.lmaxapw

    return
end

# Should be called after genrmesh
function init_packed_mtr!( mtr_vars::MuffinTinRadialVars )
    #
    nrmti = mtr_vars.nrmti
    nspecies = size(nrmti,1)
    nrmt = mtr_vars.nrmt
    #
    nrcmt = mtr_vars.nrcmt
    nrcmti = mtr_vars.nrcmti
    #
    lmmaxo = mtr_vars.lmmaxo
    lmmaxi = mtr_vars.lmmaxi
    #
    mtr_vars.npmtmax  = 1
    mtr_vars.npcmtmax = 1
    #

    for is in 1:nspecies
        #
        mtr_vars.npmti[is] = lmmaxi*nrmti[is]
        mtr_vars.npmt[is] = mtr_vars.npmti[is] + lmmaxo*(nrmt[is] - nrmti[is])
        mtr_vars.npmtmax = max(mtr_vars.npmtmax, mtr_vars.npmt[is])
        #
        mtr_vars.npcmti[is] = lmmaxi*nrcmti[is]
        mtr_vars.npcmt[is] = mtr_vars.npcmti[is] + lmmaxo*(nrcmt[is] - nrcmti[is])
        mtr_vars.npcmtmax = max(mtr_vars.npcmtmax, mtr_vars.npcmt[is])
    end
    return
end


text:
struct niw_stable_hyperparams <: distribution_hyper_params
    κ::Float64
    m::AbstractArray{Float64}
    ν::Float64
    ψ::AbstractArray{Float64}
end

mutable struct niw_stable_sufficient_statistics <: sufficient_statistics
    N::Float64
    points_sum::AbstractArray{Float64,1}
    S::AbstractArray{Float64,2}
end

function calc_posterior(prior:: niw_stable_hyperparams, suff_statistics::niw_stable_sufficient_statistics)
    if suff_statistics.N == 0
        return prior
    end
    κ = prior.κ + suff_statistics.N
    ν = prior.ν + suff_statistics.N

    m = (prior.m.*prior.κ + suff_statistics.points_sum) / κ
    ψ = (prior.ν * prior.ψ + prior.κ*prior.m*prior.m' -κ*m*m'+ suff_statistics.S) / ν
    ψ = Matrix(Hermitian(ψ))
    if  isposdef(ψ) == false
        println(ψ)
        println(m)
        println(κ*m*m')
        println(suff_statistics)
        println(prior)
    end
    return niw_stable_hyperparams(κ,m,ν,ψ)
end

#function calc_posterior(prior:: niw_hyperparams, suff_statistics::niw_sufficient_statistics)
#    κ = prior.κ + suff_statistics.N
#    m = (κ * prior.m + suff_statistics.points_sum) / κ
#    ν = prior.ν + suff_statistics.N
#    ψ = prior.ψ + suff_statistics.S + prior.κ*prior.m*prior.m'-κ*m*m'
#    return niw_hyperparams(κ,m,ν,ψ)
#end

function sample_distribution(hyperparams::niw_stable_hyperparams)
    Σ = rand(Distributions.InverseWishart(hyperparams.ν, hyperparams.ν* hyperparams.ψ))
    μ = rand(Distributions.MvNormal(hyperparams.m, Σ/hyperparams.κ))
    Σ = Matrix{Float64}(I, length(μ), length(μ)) * 1.0
    return mv_gaussian(μ,Σ,inv(Σ),logdet(Σ))
end

function create_sufficient_statistics(hyper::niw_stable_hyperparams,posterior::niw_stable_hyperparams,points::AbstractArray{Float64,2}, pts_to_group = 0)
    if size(points,2) == 0
        return niw_stable_sufficient_statistics(size(points,2),zeros(length(hyper.m)),zeros(length(hyper.m),length(hyper.m)))
    end
    pts = Array(points)
    points_sum = @view sum(pts, dims = 2)[:]
    S = pts * pts'
    return niw_stable_sufficient_statistics(size(points,2),points_sum,S)
end

function log_marginal_likelihood(hyper::niw_stable_hyperparams, posterior_hyper::niw_stable_hyperparams, suff_stats::niw_stable_sufficient_statistics)
    D = length(suff_stats.points_sum)
    logpi = log(pi)
    return -suff_stats.N*D/2*logpi +
        log_multivariate_gamma(posterior_hyper.ν/2, D)-
        log_multivariate_gamma(hyper.ν/2, D) +
         (hyper.ν/2)*logdet(hyper.ψ*hyper.ν)-
         (posterior_hyper.ν/2)*logdet(posterior_hyper.ψ*posterior_hyper.ν) +
         (D/2)*(log(hyper.κ)-(D/2)*log(posterior_hyper.κ))
end

function aggregate_suff_stats(suff_l::niw_stable_sufficient_statistics, suff_r::niw_stable_sufficient_statistics)
    return niw_stable_sufficient_statistics(suff_l.N+suff_r.N, suff_l.points_sum + suff_r.points_sum, suff_l.S+suff_r.S)
end


text:
"""
Helper type for rendering arrows in the visualizer.

Usage:

Given a `Visualizer` named `vis`, we can create an arrow with:

```
ArrowVisualizer arrow(vis)
setobject!(arrow)
```

We can then set the base and direction of the arrow with:

```
settransform!(arrow, Point(0.0, 0.0, 0.0), Vec(0.5, 0.5, 0.0))
```

Note that this `settransform!` should work correctly inside an `Animation`, as long as
you do `setobject!` before starting the animation.
"""
struct ArrowVisualizer{V<:AbstractVisualizer}
    vis::V
end

function setobject!(vis::ArrowVisualizer, material::AbstractMaterial=defaultmaterial();
        shaft_material::AbstractMaterial=material,
        head_material::AbstractMaterial=material)
    settransform!(vis, zero(Point{3, Float64}), zero(Vec{3, Float64}))
    shaft = Cylinder(zero(Point{3, Float64}), Point(0.0, 0.0, 1.0), 1.0)
    setobject!(vis.vis[:shaft], shaft, shaft_material)
    head = Cone(zero(Point{3, Float64}), Point(0.0, 0.0, 1.0), 1.0)
    setobject!(vis.vis[:head], head, head_material)
    vis
end

function Base.show(io::IO, v::ArrowVisualizer)
    print(io, "MeshCat ArrowVisualizer with path $(v.vis.path)/arrow")
end

function settransform!(vis::ArrowVisualizer, base::Point{3}, vec::Vec{3};
        shaft_radius=0.01,
        max_head_radius=2*shaft_radius,
        max_head_length=max_head_radius)
    vec_length = norm(vec)
    T = typeof(vec_length)
    rotation = if vec_length > eps(T)
        rotation_between(SVector(0., 0., 1.), vec)
    else
        one(UnitQuaternion{Float64})
    end |> LinearMap

    vis_tform = Translation(base) ∘ rotation
    settransform!(vis.vis, vis_tform)

    shaft_length = max(vec_length - max_head_length, 0)
    shaft_scaling_diag = SVector(shaft_radius, shaft_radius, shaft_length)
    if iszero(shaft_length)
        # This case is necessary to ensure that the shaft
        # completely disappears in animations.
        shaft_scaling_diag = zero(shaft_scaling_diag)
    end
    shaft_scaling = LinearMap(Diagonal(shaft_scaling_diag))
    settransform!(vis.vis[:shaft], shaft_scaling)

    head_length = vec_length - shaft_length
    head_radius = max_head_radius * head_length / max_head_length
    head_scaling = LinearMap(Diagonal(SVector(head_radius, head_radius, head_length)))
    head_tform = Translation(shaft_length * Vec(0, 0, 1)) ∘ head_scaling
    # head_tform = Translation(base) ∘ rotation ∘ Translation(shaft_length * Vec(0, 0, 1)) ∘ head_scaling
    settransform!(vis.vis[:head], head_tform)

    vis
end

function settransform!(vis::ArrowVisualizer, base::Point{3}, tip::Point{3}; kwargs...)
    settransform!(vis, base, Vec(tip - base); kwargs...)
end


text:
using VersionBenchmarks
using Test


text:
# playEkemono.jl

##文法
include("grammer.jl")

#けもの世界
const サーバル=:サーバル
const かばん=:カバン
const ラッキービースト=:ラッキービースト
const 博士=:博士
const 助手=:助手
const 図書館=:図書館
const サバンナ地方=:サバンナ地方
const ジャングル地方=:ジャングル地方
const セルリアン=:セルリアン
const アライさん=:アライさん
const ロイヤルペンギン=:ロイヤルペンギン
const コウテイペンギン=:コウテイペンギン
const ジェンツーペンギン=:ジェンツーペンギン
const イワトビペンギン=:イワトビペンギン
const フンボルトペンギン=:フンボルトペンギン
const ペパプ=:ペパプ

けもの会議=[[サーバル は かばん が 好き],
	[ラッキービースト は かばん に 話す],
	[博士 は 図書館 に 住む],
	[助手 は 図書館 に 住む]
]

qk1=[[誰 は 図書館 に 住む]]

qk2=[[サーバル は かばん を 食べる]]
qk3=[[サーバル は 誰 を 食べる]]

m1 = solveQ([誰], [[誰 は 図書館 に 住む]],けもの会議)
st1 = statements([誰], [[誰 は 図書館 に 住む]],m1)

print(stringT(st1))

findFact([誰], [サーバル は かばん を 食べる],けもの会議)
findFact([誰], [サーバル は 誰 を 食べる],けもの会議)

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
using Revise
using POMDPModels
using POMDPs
using POMDPModelTools
using StaticArrays
using Random
using DiscreteValueIteration
using GlobalApproximationValueIteration
using Flux
using Statistics
using Test

include("test/test_with_nonlinear_gfa.jl")
include("test/test_with_linear_gfa.jl")


text:
import ModiaMath

include("../../../src/language/Instantiation.jl")
include("../../../src/language/Execution.jl")
include("../../../src/language/ModiaLogging.jl")

include("../../../src/symbolic/BLTandPantelides/BLTandPantelidesUtilities.jl")
include("../../../src/symbolic/BLTandPantelides/BLTandPantelides.jl")

include("../../../src/symbolic/Utilities.jl")
include("../../../src/symbolic/ExactlyRemoveSingularities.jl")
include("../../../src/symbolic/StateSelection.jl")

include("../../../src/symbolic/DAEquations/Synchronous.jl")
include("../../../src/symbolic/DAEquations/SymbolicTransform.jl")
include("../../../src/symbolic/DAEquations/BasicStructuralTransform.jl")


text:
import Base.Iterators: product

abstract type DesignType end

struct Opt <: DesignType
end
struct Rand <: DesignType
end

const Optimize = Opt()
const Randomize = Rand()

abstract type ModelType end

struct Stat <: DesignType
end
struct Dyn <: DesignType
end

const Static = Stat()
const Dynamic = Dyn()

"""
**Model**

Creates a model object containing a log likelihood function and prior distributions

- `prior`: a vector of distribution objects for each parameter. A uniform prior is used if no
prior is passed.
- `loglike`: a function that computes the log likelihood

Constructor

````julia
Model(args...; prior=nothing, loglike, kwargs...)
````
*Examples*

```julia
# Default uninform prior
loglike(θ, design, data) = ...

model = Model(;loglike)
```

```julia
# Custom priors
prior = [Beta(5,5),...]
loglike(θ, design, data) = ...

model = Model(;loglike, prior)
```
"""
mutable struct Model{F<:Function,P}
    prior::P
    loglike::F
end

function Model(args...; prior=nothing, loglike, kwargs...)
    return Model(prior, (x...)->loglike(x..., args...; kwargs...))
end

"""
    Optimizer(;task, model, grid_design, grid_parms, grid_response)

`Optimizer` constructs a model object for adaptive design optimization

- `model`: a model object
- `grid_design`:a grid of design parameters
- `grid_parms`: a grid of model parameters
- `grid_data`: a grid of data
- `log_like`: a three dimensional array of precomputed log likelihoods
- `marg_log_like`: a two dimensional array containing marginal log likelihoods for design and data
- `priors`: a multidimensional array of prior probabilities for parameters
- `log_post`: a one dimensional array of log posterior probabilities for parameters
- `entropy`: a two dimensional array of entropy values for parameter and design combinations
-
"""
mutable struct Optimizer{A,MT,M<:Model,T1,T2,T3,T4,T5,T6,T7,T8,T9,
    T10,T11}
    design_type::A
    model_type::MT
    model::M
    design_grid::T1
    parm_grid::T2
    data_grid::T3
    log_like::Array{Float64,3}
    marg_log_like::T4
    priors::T5
    log_post::Vector{Float64}
    entropy::T6
    marg_entropy::Vector{Float64}
    cond_entropy::Vector{Float64}
    mutual_info::Vector{Float64}
    best_design::T7
    parm_names::T8
    design_names::T9
    model_state::T10
    update_state!::T11
end

function Optimizer(args...; design_type=Optimize, model, design_list, parm_list, data_list,
    state_type=nothing, model_type=Static, update_state! = nothing, kwargs...)
    design_names,design_grid = to_grid(design_list)
    parm_names,parm_grid = to_grid(parm_list)
    _,data_grid = to_grid(data_list)
    dims = map(length, (parm_grid,design_grid,data_grid))
    model_state = create_state(model_type, state_type, dims, args...; kwargs...)
    log_like = loglikelihood(model, design_grid, parm_grid, data_grid, model_type, model_state)
    priors = prior_probs(model, parm_grid)
    post = priors[:]
    log_post = log.(post)
    entropy = compute_entropy(log_like)
    marg_log_like = marginal_log_like(log_post, log_like)
    marg_entropy = marginal_entropy(marg_log_like)
    cond_entropy = conditional_entropy(entropy, post)
    mutual_info = mutual_information(marg_entropy, cond_entropy)
    best_design = find_best_design(mutual_info, design_grid, design_names)
    return Optimizer(design_type, model_type, model, design_grid, parm_grid,
        data_grid, log_like,marg_log_like, priors, log_post, entropy, marg_entropy,
        cond_entropy, mutual_info, best_design, parm_names, design_names,
        model_state, update_state!)
end


text:
using DataFrames, RDatasets

abstract NLregMod{T<:FloatingPoint}     # nonlinear regression model

type MicMen{T<:FloatingPoint} <: NLregMod{T}
    conc::Vector{T}
end
MicMen{T<:FloatingPoint}(c::DataArray{T,1}) = MicMen(vector(c))

type fMicMen <: Base.Func{2} end
Base.evaluate{T<:FloatingPoint}(::fMicMen,x::T,K::T) = x/(K+x)
#result_type{T<:FloatingPoint}(::fMicMen,::Type{T},::Type{T}) = T

StatsBase.nobs(m::MicMen) = length(m.conc)
pnames(m::MicMen) = ["Vm", "K"]
nlinpars(m::MicMen) = 1
function modelmat!{T<:FloatingPoint}(mm::Matrix{T},m::MicMen{T},nlpars::StridedVector{T})
    (size(mm) == (nobs(m),1) && length(nlpars) == 1) || error("Dimension mismatch")
    map!(fMicMen(), mm, m.conc, nlpars[1])
end
function modelmat{T<:FloatingPoint}(m::NLregMod{T},nlpars::StridedVector{T})
    modelmat!(Array(T,(nobs(m),nlinear(m))), m, nlpars)
end
function mmder!{T<:FloatingPoint}(mmd::Array{T,3},mm::Matrix{T},m::MicMen{T},nlpars::StridedVector{T})
    n = nobs(m)
    (size(mmd) == (n,1,1) && size(mm) == (n,1) && length(nlpars) == 1) || error("Dimension mismatch")
    map!(fMicMen(), mm, m.conc, nlpars[1])
end
function expctd!{T<:FloatingPoint}(mu::Vector{T},m::NLregMod{T},pars::Vector{T})
    nl = nlinear(m)
    BLAS.gemv!('N', one(T),
               modelmat(m,sub(pars,(nl+1):length(pars))),
               sub(pars,1:nl), 0., mu)
end
expctd{T<:FloatingPoint}(m::NLregMod{T},pars::Vector{T}) = expctd!(Array(T,(nobs(m),)),m,pars)
npars{T<:FloatingPoint}(m::NLregMod{T}) = length(pnames(m))
nnlpars{T<:FloatingPoint}(m::NLregMod{T}) = length(pnames(m)) - nlinpars(m)

type plinear{T<:FloatingPoint}
    m::NLregMod{T}
    y::Vector{T}
    mu::Vector{T}
    pars::Vector{T}
    mm::Matrix{T}
    tr::Matrix{T}
end
function plinear{T<:FloatingPoint}(m::NLregMod{T},y::Vector{T})
    n = length(y); nobs(m) == n || error("Dimension mismatch")
    nl = nlinpars(m); mm = Array(T, n, nl)
    plinear(m, y, Array(T,n), Array(T,length(pnames(m))), mm, similar(mm), zeros(T,nl,nl))
end
plinear{T<:FloatingPoint}(m::NLregMod{T}, y::DataArray{T,1}) = plinear(m,vector(y))

function StatsBase.deviance{T<:FloatingPoint}(pl::plinear{T},nlpars::Vector{T})
    m = pl.m; nl = size(pl.mm,2); # number of conditionally linear parameters
    lin = 1:nl; nonlin = nl + (1:(length(pnames(m)) - nl))
    copy!(sub(pl.pars, nlin), nlpars)   # record the current values
    _,pl.tr = qrfact!(modelmat!(pl.vs,m,nlpars)) # decompose model matrix for linear pars
    LAPACK.gemqrt!('L','T',pl.vs,pl.tr,copy!(pl.mu,pl.y)) # create Q'y in pl.mu
    BLAS.trsv!('U','N','N',sub(pl.vs,lin,lin),copy!(sub(pl.pars,lin),sub(pl.mu,lin))) # solve
    fill!(sub(pl.mu, (nl+1):nobs(m)), zero(T))
    LAPACK.gemqrt!('L','N',pl.vs,pl.tr,pl.mu)
    sqdiffsum(pl.y,pl.mu)
end

abstract Ptransform

type VClka2Vkka <: Ptransform end
function pars!{T<:FloatingPoint}(t::VClka2Vkka,out::Matrix{T},in::Matrix{T})
    copy!(out, in)
    n,p = size(out); p == 3 || error("3 columns required in 'in' and 'out'")
    for i in 1:n
        @inbounds out[i,2] = in[i,2]/in[i,1]
    end
    out
end
pars{T<:FloatingPoint}(t::VClka2Vkka,in::Matrix{T})=pars!(t,similar(in),in)
function parsjac!{T<:FloatingPoint}(::VClka2Vkka,out::Matrix{T},jac::Array{T,3},in::Matrix{T})
    copy!(out, in); n,p = size(out)
    p == 3 || error("3 columns required in 'in' and 'out'")
    (n,p,p) == size(jac) || error("dimension mismatch")
    for i in 1:n
        V = in[i,1]
        Cl = in[i,2]
        out[i,2] = Cl/V
        jac[i,2,1] = -Cl/(V * V)
        jac[i,2,2] = 1. / V
    end
    out,jac
end
pnmsin(::Type{VClka2Vkka}) = ["V","Cl","ka"]
pnmsout(::Type{VClka2Vkka}) = ["V","k","ka"]

# Oral, single dose, 1 compartment using k and ka as parameters.
type OralSd1Vkka{T<:FloatingPoint} <: NLRegFD
    time::Vector{T}
    dose::Float64
end
StatsBase.nobs(m::OralSd1Vkka) = length(m.time)
pnames(m::OralSd1Vkka) = ["V","k","ka"]

# Evaluation functor.  Result should be multiplied by dose/V
type fOralSd1kka <: Base.Func{3} end
function Base.evaluate{T<:FloatingPoint}(::fOralSd1kka,t::T,k::T,ka::T)
    ka*(exp(-k*t)-exp(-ka*t))/(ka-k)
end
#result_type{T<:FloatingPoint}(::fOralSd1kka,::Type{T},::Type{T},::Type{T}) = T

function expctd!{T<:FloatingPoint}(mu::Vector{T},m::OralSd1Vkka,x::Array{T})
    map1!(Multiply(), map!(fOralSd1kka(), mu, m.time, elcol(x,2), elcol(x,3), m.dose ./ elcol(x,1)))
end
# PK model for single oral dose, 1 compartment with parameters V, Cl and ka
type OralSd1VClka <: NLRegFD
    time::Vector{Float64}
    dose::Float64
end
OralSd1VClka(time::Vector{Float64}) = OralSd1VClka(time,1.) # default is unit dose
StatsBase.nobs(m::OralSd1VClka) = length(m.time)
pnames(m::OralSd1VClka) = ["V","Cl","ka"]
npar(m::NonlinearRegModel) = length(pnames(m))

function expctd!{T<:Float64}(mu::Vector{T}, m::OralSd1VClka, x::Vector{T})
    length(x) == 3 || error("trailing dimension of x must be 3")
    map1!(Multiply(), map!(fOralSd1kka(), mu, m.time, x[2]/x[1], x[3]), m.dose ./ x[1])
end

function expctd!{T<:Float64}(eta::Vector{T}, m::OralSd1VClka, x::Matrix{T})
    n,p = size(x); p == 3 || error("trailing dimension of x must be 3")
#    V = x[1,:]; Cl = x[2,:]; ka = x[3,:]; k = Cl./V;
    map1!(Multiply(), map!(fOralSd1kka(),m.time,x[2,:]./x[1,:],x[3,:]), m.dose./x[1,:])
end
expctd{T<:Float64}(m::NonlinearRegModel,x::Vector{T}) = expctd!(Array(T,(nobs(m),)),m,x)
expctd(m::NonlinearRegModel, x::Matrix{Float64}) = expctd!(Array(T,(nobs(m),)),m,x)

function expctdjac!{T<:Float64}(mu::Vector{T},jac::Matrix{T},m::OralSd1VClka,x::Vector{T})
    length(x) == 3 || error("length(x) = $(length(x)), should be 3")
    t = m.time; n = length(t); d = m.dose;
    V = x[1]
    Cl = x[2]
    ka = x[3]
    k = Cl/V                            # e2
    e1 = d/V
    e3 = ka - k
    e4 = ka / e3
    e5 = e1 * e4
    e14 = V * V
    e15 = Cl/e14
    e20 = e3 * e3
    e28 = 1./V
    for i in 1:n
        ti = t[i]
        e8 = exp(-k*ti)
        e11 = exp(-ka*ti)
        e12 = e8 - e11
        expctd[i] = e5*e12
        jac[i,1] = e5*(e8*(e15*ti)) - (e1*(ka*e15/e20) + d/e14*e4)*e12
        jac[i,2] = e1 * (ka * e28/e20) * e12 - e5 * (e8 * (e28 * ti))
        jac[i,3] = e1 * (1/e3 - ka/e20) * e12 + e5 * (e11 * ti)
    end
    mu, jac
end
const step = sqrt(eps())
const steps = [-step, step]
const mults = 1. + steps
function expctdjac(m::NLRegFD, x::Vector{Float64})
    p = length(x); pred = expctd(m, x); n = length(pred)
    jac = zeros(n,p)
    for j in 1:p
        par = copy(x); pjs = x[j] == 0. ? steps : x[j] * mults
        par[j] = pjs[2]
        rj = expctd(m, par)
        par[j] = pjs[1]
        jac[:,j] = (rj - expctd(m, par))/diff(pjs)
    end
    pred, jac
end
m = OralSd1VClka([0.25, 0.57, 1.12, 2.02, 3.82, 5.1, 7.03, 9.05, 12.12, 24.37], 4.02)

type NonlinearLS                    # nonlinear least squares problems
    pars::Vector{Float64}
    incr::Vector{Float64}
    obs::Vector{Float64}
    eta::Vector{Float64}                # expected response
    qtr::Vector{Float64}
    jacob::Matrix{Float64}              # Jacobian matrix
    qr::QR{Float64}
    m::NonlinearRegModel
    rss::Float64 # residual sum of squares at last successful iteration
    tolsqr::Float64    # squared tolerance for orthogonality criterion
    minfac::Float64
    mxiter::Int
    fit::Bool
end

function NonlinearLS(m::NonlinearRegModel, obs::Vector{Float64}, init::Vector{Float64})
    n = length(obs); p = length(init);
    eta, jacob = expctdjac(m, init)
    resid = obs - eta
    NonlinearLS(init, zero(init), obs, eta, resid, zero(resid), jacob,
                qrfact(jacob), m, sum(resid.^2), 1e-8, 0.5^10, 1000, false)
end
nlm = NonlinearLS(m, [2.84, 6.57, 10.5, 9.66, 8.58, 8.36, 7.47, 6.89, 5.94, 3.28], exp([-1,-4,0.55]))
# evaluate expected response and residual at m.pars + fac * m.incr
# return residual sum of squares
function updtres(nl::NonlinearLS, fac::Float64)
    expctd!(nl.m, nl.eta, nl.pars + fac * nl.incr)
    sqdiffsum(nl.obs,nl.eta)
end

# Create the QR factorization, qtr = Q'resid, solve for the increment and
# return the numerator of the squared convergence criterion
function qtr(nl::NonlinearLS)
    vs = nl.qr.vs; inc = nl.incr; qt = nl.qtr
    copy!(vs, nl.jacob); copy!(qt, nl.resid)
    _, T = LAPACK.geqrt3!(vs)
    copy!(nl.qr.T,T)
    LAPACK.gemqrt!('L','T',vs,T,qt)
    s = 0.; p = size(vs,2)
    for i in 1:p qti = qt[i]; s += qti * qti; inc[i] = qti  end
    BLAS.trsv!('U','N','N',sub(vs,1:p,1:p),inc)
    s
end

function gnfit(nl::NonlinearLS)          # Gauss-Newton nonlinear least squares
    if !nl.fit
        converged = false; rss = nl.rss
        for i in 1:nl.mxiter
            crit = qtr(nl)/nl.rss # evaluate increment and orthogonality cvg. crit.
            converged = crit < nl.tolsqr
            f = 1.
            while f >= nl.minfac
                rss = updtres(nl,f)
                if rss < nl.rss break end
                f *= 0.5                    # step-halving
            end
            if f < nl.minfac
                error("Failure to reduce rss at $(nl.pars) with incr = $(nl.incr) and minfac = $(nl.minfact)")
            end
            nl.rss = rss
            nl.pars += f * nl.incr
            if converged break end
            nl.eta, nl.jacob = expctdjac(nl.m, nl.pars)  # evaluate Jacobian
        end
        if !converged error("failure to converge in $(nl.mxiter) iterations") end
        nl.fit = true
    end
    nl
end

function Base.show(io::IO, nl::NonlinearLS)
    gnfit(nl)
    n,p = size(nl.jacob)
    s2 = nl.rss/float(n-p)
    varcov = s2 * symmetrize!(potri!('U', nl.qr.vs[1:p,:])[1],'U')
    stderr = sqrt(diag(varcov))
    t_vals = nl.pars./stderr
    println(io, "Model fit by nonlinear least squares to $n observations\n")
    println(io, DataFrame(parameter=nl.pars,stderr=stderr,t_value=nl.pars./stderr))
    println("Residual sum of squares at estimates = $(nl.rss)")
    println("Residual standard error = $(sqrt(s2)) on $(n-p) degrees of freedom")
end


text:
## A small driver file to run the 4
# synthetic data tests in parallel
# data calibrated to the RossmanKaggleData
#Run like this
#julia -p 4 -L syntheticDataHarness.jl test_SyntheticData.jl numRuns outPathStub d usePoisson
#passed arguments
	#ARGS[1] numRuns
	#ARGs[2] partial path for output.
	#ARGS[3] d
	#ARGS[4] bool for using the poisson amount of data assumption.

const numRuns = parse(Int, ARGS[1])
const spath = ARGS[2]
const param_path = "../RossmanKaggleData/CleanedData/"
const d = parse(Int, ARGS[3])
const usePoisson = parse(Bool, ARGS[4])

const s = .95
K_grid = vcat(round.(Int, 2 .^(3:.25:10)), 1115)
N_grid = [10]

outPath = "$(spath)_syntheticRossman_$(s)_$(usePoisson)_$(4*numRuns)"
#old naming specification.  keep for a bit.
#outPath = "$(spath)_Ross_$(maximum(K_grid))_$(d)_$(N)_$(s)_$(usePoisson)_$(4*numRuns)"

#First read in the data and parse it appropriately
#Do this on a single processor bc it should be fast.
ps_full = readdlm("../RossmanKaggleData/CleanedData/ps_full$(d).csv", ',')
supp_full = readdlm("../RossmanKaggleData/CleanedData/support$(d).csv", ',')

start_time = time_ns()
file_a = @spawn convInKtest(numRuns, K_grid, supp_full, ps_full, "$(outPath)_a_", N_grid, s, seed=8675309, usePoisson=usePoisson)
file_b = @spawn convInKtest(numRuns, K_grid, supp_full, ps_full, "$(outPath)_b_", N_grid, s, seed=5167462266, usePoisson=usePoisson)
file_c = @spawn convInKtest(numRuns, K_grid, supp_full, ps_full, "$(outPath)_c_", N_grid, s, seed=5164174290, usePoisson=usePoisson)
file_d = @spawn convInKtest(numRuns, K_grid, supp_full, ps_full, "$(outPath)_d_", N_grid, s, seed=112456059, usePoisson=usePoisson)

# ######
file_a = fetch(file_a)
file_b = fetch(file_b)
file_c = fetch(file_c)
file_d = fetch(file_d)

time_stamp = (time_ns() - start_time) * 1e-9

##read everyone in, throw away a line
data, header = readdlm(file_a, ',', header=true)

data_t = readdlm(file_b, ',', skipstart=1)
data_t[:, 1] .+= numRuns
data = vcat(data, data_t)

data_t = readdlm(file_c, ',', skipstart=1)
data_t[:, 1] .+= 2numRuns
data = vcat(data, data_t)

data_t = readdlm(file_d, ',', skipstart=1)
data_t[:, 1] .+= 3numRuns
data = vcat(data, data_t)

f = open("$(outPath).csv", "w")
writedlm(f, header, ',')
writedlm(f, data, ',')
close(f)

println("Num Paths: \t $(4*numRuns) \t Time:", time_stamp )

max_stars_repo_name: CausalML/DataPooling

max_issues_repo_name: CausalML/DataPooling

max_forks_repo_name: CausalML/DataPooling


text:
@testset "conic linear problems" begin
    @testset "conic_linear1_test" begin
        #=
        primal
            min -3x - 2y - 4z
        s.t.
            x +  y +  z == 3 :w_1
            y +  z == 2      :w_2
            x>=0 y>=0 z>=0

        dual
            max 3w_1 + 2w_2
        s.t.
            [-w_1 - 3.0,
             -w_1 - w_2 - 2.0,
             -w_1 - w_2 - 4.0] in Nonnegatives
        =#
        primal_model = conic_linear1_test()
        dual_model, primal_dual_map = dual_model_and_map(primal_model)

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 2
        list_of_cons = MOI.get(dual_model, MOI.ListOfConstraints())
        @test list_of_cons == [(VAF{Float64}, MOI.Nonnegatives)]
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VAF{Float64},MOI.Nonnegatives}(),
        ) == 1
        obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
        @test obj_type == SAF{Float64}
        @test MOI.get(dual_model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
        @test MOI.constant(obj) == 0.0
        @test Set(MOI.coefficient.(obj.terms)) == Set([3.0; 2.0])

        ci = first(
            MOI.get(
                dual_model,
                MOI.ListOfConstraintIndices{VAF{Float64},MOI.Nonnegatives}(),
            ),
        )
        eq_con1_fun = MOI.get(dual_model, MOI.ConstraintFunction(), ci)
        @test MOI.coefficient.(eq_con1_fun.terms) == -ones(5)
        @test MOI.constant(eq_con1_fun) == [-3.0, -2.0, -4.0]
        @test MOI.get(dual_model, MOI.ConstraintSet(), ci) ==
              MOI.Nonnegatives(3)

        primal_con_dual_var = primal_dual_map.primal_con_dual_var
        ci_zero = first(
            MOI.get(
                primal_model,
                MOI.ListOfConstraintIndices{VAF{Float64},MOI.Zeros}(),
            ),
        )
        @test primal_con_dual_var[ci_zero] == VI.(1:2)
        ci_nneg = first(
            MOI.get(
                primal_model,
                MOI.ListOfConstraintIndices{VVF,MOI.Nonnegatives}(),
            ),
        )
        @test !haskey(primal_con_dual_var, ci_nneg)
        @test primal_dual_map.constrained_var_dual[ci_nneg] == ci

        for i in 1:3
            vi = VI(i)
            @test !haskey(primal_dual_map.primal_var_dual_con, vi)
            @test primal_dual_map.constrained_var_idx[vi] == (ci_nneg, i)
        end
    end

    @testset "conic_linear3_test" begin
        #=
        primal
            min  3x + 2y - 4z + 0s
        s.t.
            x           -  s == -4    :w_1
                y            == -3    :w_2
            x      +  z      == 12    :w_3
            y <= 0
            z >= 0
            s zero

        dual
            max -4w_4 - 3w_5 + 12w_6
        s.t
            w_1 + w_3 == 3
            [-w_2 + 2] in Nonpositives
            [-w_3 - 4] in Nonnegatives
        =#
        primal_model = conic_linear3_test()
        dual_model, primal_dual_map = dual_model_and_map(primal_model)

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 3
        list_of_cons = MOI.get(dual_model, MOI.ListOfConstraints())
        @test Set(list_of_cons) == Set(
            [
                (VAF{Float64}, MOI.Nonpositives)
                (VAF{Float64}, MOI.Nonnegatives)
                (SAF{Float64}, MOI.EqualTo{Float64})
            ],
        )
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.EqualTo{Float64}}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VAF{Float64},MOI.Nonnegatives}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VAF{Float64},MOI.Nonpositives}(),
        ) == 1
        obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
        @test obj_type == SAF{Float64}
        @test MOI.get(dual_model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
        @test MOI.constant(obj) == 0.0
        @test Set(MOI.coefficient.(obj.terms)) == Set([-4.0; 12.0; -3.0])

        ci_eq = first(
            MOI.get(
                dual_model,
                MOI.ListOfConstraintIndices{SAF{Float64},MOI.EqualTo{Float64}}(),
            ),
        )
        eq_con1_fun = MOI.get(dual_model, MOI.ConstraintFunction(), ci_eq)
        eq_con1_set = MOI.get(dual_model, MOI.ConstraintSet(), ci_eq)
        @test MOI.coefficient.(eq_con1_fun.terms) == [1.0; 1.0]
        @test MOI.constant.(eq_con1_fun) == 0.0
        @test MOI.constant(eq_con1_set) == 3.0
        ci_np = first(
            MOI.get(
                dual_model,
                MOI.ListOfConstraintIndices{VAF{Float64},MOI.Nonpositives}(),
            ),
        )
        eq_con2_fun = MOI.get(dual_model, MOI.ConstraintFunction(), ci_np)
        @test MOI.get(dual_model, MOI.ConstraintSet(), ci_np) ==
              MOI.Nonpositives(1)
        @test MOI.coefficient.(eq_con2_fun.terms) == -[1.0]
        @test MOI.constant(eq_con2_fun) == [2.0]
        ci_nn = first(
            MOI.get(
                dual_model,
                MOI.ListOfConstraintIndices{VAF{Float64},MOI.Nonnegatives}(),
            ),
        )
        eq_con3_fun = MOI.get(dual_model, MOI.ConstraintFunction(), ci_nn)
        @test MOI.get(dual_model, MOI.ConstraintSet(), ci_nn) ==
              MOI.Nonnegatives(1)
        @test MOI.coefficient.(eq_con3_fun.terms) == -[1.0]
        @test MOI.constant(eq_con3_fun) == [-4.0]

        primal_con_dual_var = primal_dual_map.primal_con_dual_var
        ci_zero = first(
            MOI.get(primal_model, MOI.ListOfConstraintIndices{VVF,MOI.Zeros}()),
        )
        @test !haskey(primal_con_dual_var, ci_zero)
        ci_nneg = first(
            MOI.get(
                primal_model,
                MOI.ListOfConstraintIndices{VVF,MOI.Nonnegatives}(),
            ),
        )
        @test !haskey(primal_con_dual_var, ci_nneg)
        @test primal_dual_map.constrained_var_dual[ci_nneg] == ci_nn
        ci_npos = first(
            MOI.get(
                primal_model,
                MOI.ListOfConstraintIndices{VVF,MOI.Nonpositives}(),
            ),
        )
        @test !haskey(primal_con_dual_var, ci_npos)
        @test primal_dual_map.constrained_var_dual[ci_npos] == ci_np
        ci_aff_zero = first(
            MOI.get(
                primal_model,
                MOI.ListOfConstraintIndices{VAF{Float64},MOI.Zeros}(),
            ),
        )
        @test primal_con_dual_var[ci_aff_zero] == VI.(1:3)

        primal_var_dual_con = primal_dual_map.primal_var_dual_con
        @test primal_var_dual_con[VI(1)] == ci_eq
        for i in 2:4
            vi = VI(i)
            @test !haskey(primal_var_dual_con, vi)
        end
        @test primal_dual_map.constrained_var_idx[VI(2)] == (ci_npos, 1)
        @test primal_dual_map.constrained_var_idx[VI(3)] == (ci_nneg, 1)
        @test primal_dual_map.constrained_var_idx[VI(4)] == (ci_zero, 1)
    end
end


text:
# array methods
linspace(a,b,n) = LinRange(a,b,n) |> collect

"""
    x = xvec(λ,N)

Create `x` values with correct periodicity for box specified by length `λ`, using `N` points.
"""
xvec(L,N) = LinRange(-L/2,L/2,N+1)[2:end] |> collect

"""
    k = kvec(λ,N)

Create `k` values with correct periodicity for box specified by length `λ` for number of points `N`.
"""
function kvec(λ,N)
    # @assert iseven(N)
    # nk = 0:Int(N/2)
    # k = [nk[1:end-1];-reverse(nk[2:end])]*2*π/λ
    k = fftfreq(N)*N*2*π/λ    # N not restricted
    return k
end

"""
    X = xvecs(L,N)

Create a tuple containing the spatial coordinate array for each spatial dimension.
"""
function xvecs(L,N)
    X = []
    for (λ,ν) in zip(L,N)
        x = xvec(λ,ν)
        push!(X,x)
    end
    return X |> Tuple
end

"""
    K = kvecs(L,N)

Create a tuple containing the spatial coordinate array for each spatial dimension.
"""
function kvecs(L,N)
    K = []
    for (λ,ν) in zip(L,N)
        k = kvec(λ,ν)
        push!(K,k)
    end
    return K |> Tuple
end

"""
    k² = k2(K)

Create the kinetic energy array `k²` on the `k`-grid defined by the tuple `K`.
"""
function k2(K)
    kind = Iterators.product(K...)
    return map(k-> sum(abs2.(k)),kind)
end

"""
    X,K,dX,dK = makearrays(L,N)

Create all `x` and `k` arrays for box specified by tuples `L=(Lx,...)` and `N=(Nx,...)`.
For convenience, differentials `dX`, `dK` are also reaturned. `L` and `N` must be tuples of equal length.
"""
function makearrays(L,N)
    @assert length(L) == length(N)
    X = xvecs(L,N)
    K = kvecs(L,N)
    dX = Float64[]; dK = Float64[]
    for j ∈ eachindex(X)
        x = X[j]; k = K[j]
        dx = x[2]-x[1]; dk = k[2]-k[1]
        push!(dX,dx)
        push!(dK,dk)
    end
    dX = dX |> Tuple
    dK = dK |> Tuple
    return X,K,dX,dK
end

"""
    A = crandn_array(M)

Make placeholder `2x2x...` complex `randn()` array of `M` dimensions."""
function crandn_array(M)
    a = Int.(ones(M)).+1
    return randn(a...) |> complex
end

"""
    A = crandnpartition(D,M)

Make placeholder ArrayPartition vector of length `M`, containing `2x2x...` rank D complex matrices.
"""
function crandnpartition(D,N)
    a = crandn_array(D)
    args = []
    for j = 1:N
        push!(args,a)
    end
    return ArrayPartition(args...)
end


text:
abstract type AbstractLoadFormulation <: AbstractDeviceFormulation end
abstract type AbstractControllablePowerLoadFormulation <: AbstractLoadFormulation end
struct StaticPowerLoad <: AbstractLoadFormulation end
struct InterruptiblePowerLoad <: AbstractControllablePowerLoadFormulation end
struct DispatchablePowerLoad <: AbstractControllablePowerLoadFormulation end

########################### dispatchable load variables ####################################
function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    ::PSIContainer,
) where {T <: ActivePowerVariable, U <: PSY.ElectricLoad}
    return AddVariableSpec(;
        variable_name = make_variable_name(T, U),
        binary = false,
        expression_name = :nodal_balance_active,
        sign = -1.0,
        lb_value_func = x -> 0.0,
        ub_value_func = x -> PSY.get_max_active_power(x),
    )
end

function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    ::PSIContainer,
) where {T <: ReactivePowerVariable, U <: PSY.ElectricLoad}
    return AddVariableSpec(;
        variable_name = make_variable_name(T, U),
        binary = false,
        expression_name = :nodal_balance_reactive,
        sign = -1.0,
        lb_value_func = x -> 0.0,
        ub_value_func = x -> PSY.get_max_reactive_power(x),
    )
end

function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    ::PSIContainer,
) where {T <: OnVariable, U <: PSY.ElectricLoad}
    return AddVariableSpec(; variable_name = make_variable_name(T, U), binary = true)
end

####################################### Reactive Power Constraints #########################
"""
Reactive Power Constraints on Controllable Loads Assume Constant power_factor
"""
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ReactivePowerVariable},
    ::Type{<:PSY.ElectricLoad},
    ::Type{<:AbstractControllablePowerLoadFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
)
    return DeviceRangeConstraintSpec(;
        custom_psi_container_func = custom_reactive_power_constraints!,
    )
end

function custom_reactive_power_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::Type{<:AbstractControllablePowerLoadFormulation},
) where {T <: PSY.ElectricLoad}
    time_steps = model_time_steps(psi_container)
    constraint = JuMPConstraintArray(undef, [PSY.get_name(d) for d in devices], time_steps)
    assign_constraint!(psi_container, REACTIVE, T, constraint)

    for t in time_steps, d in devices
        name = PSY.get_name(d)
        pf = sin(atan((PSY.get_max_reactive_power(d) / PSY.get_max_active_power(d))))
        reactive = get_variable(psi_container, REACTIVE_POWER, T)[name, t]
        real = get_variable(psi_container, ACTIVE_POWER, T)[name, t] * pf
        constraint[name, t] = JuMP.@constraint(psi_container.JuMPmodel, reactive == real)
    end
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerVariable},
    ::Type{T},
    ::Type{<:DispatchablePowerLoad},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.ElectricLoad}
    if (!use_parameters && !use_forecasts)
        return DeviceRangeConstraintSpec(;
            range_constraint_spec = RangeConstraintSpec(;
                constraint_name = make_constraint_name(
                    RangeConstraint,
                    ActivePowerVariable,
                    T,
                ),
                variable_name = make_variable_name(ActivePowerVariable, T),
                limits_func = x -> (min = 0.0, max = PSY.get_active_power(x)),
                constraint_func = device_range!,
                constraint_struct = DeviceRangeConstraintInfo,
            ),
        )
    end

    return DeviceRangeConstraintSpec(;
        timeseries_range_constraint_spec = TimeSeriesConstraintSpec(
            constraint_name = make_constraint_name(RangeConstraint, ActivePowerVariable, T),
            variable_name = make_variable_name(ActivePowerVariable, T),
            parameter_name = use_parameters ? ACTIVE_POWER : nothing,
            forecast_label = "max_active_power",
            multiplier_func = x -> PSY.get_max_active_power(x),
            constraint_func = use_parameters ? device_timeseries_param_ub! :
                              device_timeseries_ub!,
        ),
    )
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerVariable},
    ::Type{T},
    ::Type{<:InterruptiblePowerLoad},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.ElectricLoad}
    if (!use_parameters && !use_forecasts)
        return DeviceRangeConstraintSpec(;
            range_constraint_spec = RangeConstraintSpec(;
                constraint_name = make_constraint_name(
                    RangeConstraint,
                    ActivePowerVariable,
                    T,
                ),
                variable_name = make_variable_name(ActivePowerVariable, T),
                bin_variable_names = [make_variable_name(OnVariable, T)],
                limits_func = x -> (min = 0.0, max = PSY.get_active_power(x)),
                constraint_func = device_semicontinuousrange!,
                constraint_struct = DeviceRangeConstraintInfo,
            ),
        )
    end

    return DeviceRangeConstraintSpec(;
        timeseries_range_constraint_spec = TimeSeriesConstraintSpec(
            constraint_name = make_constraint_name(RangeConstraint, ActivePowerVariable, T),
            variable_name = make_variable_name(ActivePowerVariable, T),
            bin_variable_name = make_variable_name(OnVariable, T),
            parameter_name = use_parameters ? ON : nothing,
            forecast_label = "max_active_power",
            multiplier_func = x -> PSY.get_max_active_power(x),
            constraint_func = use_parameters ? device_timeseries_ub_bigM! :
                              device_timeseries_ub_bin!,
        ),
    )
end

########################## Addition to the nodal balances ##################################

function NodalExpressionSpec(
    ::Type{T},
    ::Type{<:PM.AbstractPowerModel},
    use_forecasts::Bool,
) where {T <: PSY.ElectricLoad}
    return NodalExpressionSpec(
        "max_active_power",
        REACTIVE_POWER,
        use_forecasts ? x -> PSY.get_max_reactive_power(x) : x -> PSY.get_reactive_power(x),
        -1.0,
        T,
    )
end

function NodalExpressionSpec(
    ::Type{T},
    ::Type{<:PM.AbstractActivePowerModel},
    use_forecasts::Bool,
) where {T <: PSY.ElectricLoad}
    return NodalExpressionSpec(
        "max_active_power",
        ACTIVE_POWER,
        use_forecasts ? x -> PSY.get_max_active_power(x) : x -> PSY.get_active_power(x),
        -1.0,
        T,
    )
end

############################## FormulationControllable Load Cost ###########################
function AddCostSpec(
    ::Type{T},
    ::Type{DispatchablePowerLoad},
    ::PSIContainer,
) where {T <: PSY.ControllableLoad}
    cost_function = x -> isnothing(x) ? 1.0 : PSY.get_variable(x)
    return AddCostSpec(;
        variable_type = ActivePowerVariable,
        component_type = T,
        variable_cost = cost_function,
        multiplier = OBJECTIVE_FUNCTION_NEGATIVE,
    )
end

function AddCostSpec(
    ::Type{T},
    ::Type{InterruptiblePowerLoad},
    ::PSIContainer,
) where {T <: PSY.ControllableLoad}
    cost_function = x -> isnothing(x) ? 1.0 : PSY.get_fixed(x)
    return AddCostSpec(;
        variable_type = OnVariable,
        component_type = T,
        fixed_cost = cost_function,
        multiplier = OBJECTIVE_FUNCTION_NEGATIVE,
    )
end


text:
# SPDX-License-Identifier: X11
# 2020-12-22
# Day 6, Part 2

function main()
  input = split.(
    split(readchomp(stdin), "\n\n"; keepempty=false),
    "\n";
    keepempty=false)

  (println ∘ sum ∘ map)(countanswers, input)
end

function countanswers(group)::Integer
  dict = (Dict ∘ zip)('a':'z', trues(26))

  for person ∈ group
    curr = (Dict ∘ zip)('a':'z', falses(26))
    for c ∈ person
      curr[c] = true
    end

    for (c, b) ∈ curr
      !b && (dict[c] = false)
    end
  end

  return (count ∘ values)(dict)
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie
CairoMakie.activate!() # HIDE
let
    x = y = -10:0.11:10
    y1d =  sin.(x) ./ x
    lower = y1d .- 0.1
    upper = y1d .+ 0.1
    fig = Figure(resolution = (600, 400))
    ax = Axis(fig, xlabel = "x", ylabel = "y")
    lines!(x, y1d, color = :black)
    band!(x, lower, upper; color = ("#E69F00", 0.2))
    fig[1,1] = ax
    fig
    #save("ConfidenceRegion.png"), fig, px_per_unit = 2.0)
    save(joinpath(@__DIR__, "output", "ConfidenceRegion.png"), fig, px_per_unit = 2.0) # HIDE
end

using Pkg # HIDE
Pkg.status("CairoMakie") # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
function _read(valfn)
    try; return parse(Int, read(valfn, String))
    catch ingored end
    return 0
end

let
    @info("Running race tests")
    # must match the values on race_script.jl
    lkfn = joinpath(@__DIR__, "lock")
    valfn = joinpath(@__DIR__, "sum.txt")
    logfn = joinpath(@__DIR__, "log.txt")
    try
        script = joinpath(@__DIR__, "race_script.jl")
        @assert isfile(script)

        rm(valfn; force = true)
        rm(lkfn; force = true)
        rm(logfn; force = true)
        write(valfn, "0")

        currptoj = Base.current_project(@__DIR__)

        nprocs = 20
        @info("Spawning $(nprocs) competing processes")
        for t in 1:nprocs

            # run this for debug
            # julia_cmd = strip(string(Base.julia_cmd()), ['`'])
            # plogfn = joinpath(@__DIR__, "log$(t).txt")
            # jlsrc = "$(julia_cmd) --project=$(currptoj) --startup-file=no $(script) 2>&1 > $(plogfn)"
            # run(`bash -c $(jlsrc)`; wait = false)

            julia_cmd = Base.julia_cmd()
            run(`$(julia_cmd) --project=$(currptoj) --startup-file=no $(script)`; wait = false)

            sleep(0.1)
        end
        # wait
        mt0 = -1
        val = 0
        N = nprocs * 10 # must match the values on race_script.jl

        @info("Reading")
        print(val, "/", N, "\r")
        for _ in 1:nprocs
            while true
                sleep(5.0)

                val = _read(valfn)
                mt = mtime(valfn)
                iszero(mt) && continue
                val != 0 && (mt == mt0) && break
                mt0 = mt

                print(val, "/", N, "\r")
            end
            (val == N) && break

            println(val, "/", N)
            println("waiting...")
        end
        println(val, "/", N)
        ok_res = val >= N * 0.99 # > 99% success
        @test ok_res

        # deb info
        if !ok_res && isfile(logfn)
            println("\n\n", "-"^60)
            println("LOG", "\n")
            println(read(logfn, String))
            println("\n\n", "-"^60)
        end

    finally
        rm(lkfn; force = true)
        rm(valfn; force = true)
        rm(logfn; force = true)
    end
end


text:
#
#      Rn - The manifold of the n-dimensional (real valued) Euclidean space
#
# Manopt.jl, R. Bergmann, 2019
import LinearAlgebra: I, norm
import Base: exp, log, show
export Euclidean, RnPoint, RnTVector
export distance, exp, log, norm, dot, manifoldDimension, show, getValue
export zeroTVector, tangentONB, randomMPoint, randomTVector
export validateMPoint, validateTVector, typeofMPoint, typeofTVector
# Types
# ---

@doc doc"""
    Euclidean <: Manifold

The manifold $\mathcal M = \mathbb R^n$ of the $n$-dimensional Euclidean vector
space. We employ the notation $\langle\cdot,\cdot\rangle$ for the inner product
and $\lVert\cdot\rVert_2$ for its induced norm.

# Abbreviation

`Rn`

# Constructor

    Euclidean(n)

construct the n-dimensional Euclidean space $\mathbb R^n$.
"""
struct Euclidean <: Manifold
  name::String
  dimension::Int
  abbreviation::String
  Euclidean(dimension::Int) = new("$dimension-dimensional Euclidean space",dimension,"R$dimension")
end

@doc doc"""
    RnPoint <: MPoint

the point $x\in\mathcal M$ for $\mathcal M=\mathbb R^n$ represented by an
$n$-dimensional `Vector{T}`, where `T <: AbstractFloat`.
"""
struct RnPoint{T<:AbstractFloat} <: MPoint
  value::Vector{T}
  RnPoint{T}( value::Vector{T} ) where T<:AbstractFloat = new(value)
  RnPoint{T}( value::T ) where T<:AbstractFloat = new([value])
end
RnPoint(value::T) where {T <: AbstractFloat} = RnPoint{T}(value)
RnPoint(value::Vector{T}) where {T <: AbstractFloat} = RnPoint{T}(value)
getValue(x::RnPoint) = length(x.value)==1 ? x.value[1] : x.value

@doc doc"""
    RnTVector <: TVector

the point $\xi\in\mathcal M$ for $\mathcal M=\mathbb R^n$ represented by an
$n$-dimensional `Vector{T}`, where `T <: AbstractFloat`.
"""
struct RnTVector{T <: AbstractFloat}  <: TVector
  value::Vector{T}
  RnTVector{T}(value::Vector{T})  where {T <: AbstractFloat}  = new(value)
  RnTVector{T}(value::T) where {T <: AbstractFloat}  = new([value])
end
RnTVector(value::T) where {T <: AbstractFloat} = RnTVector{T}(value)
RnTVector(value::Vector{T})  where {T <: AbstractFloat}  = RnTVector{T}(value)

getValue(ξ::RnTVector) = length(ξ.value)==1 ? ξ.value[1] : ξ.value

# Traits
# ---
# (a) Rn is a MatrixManifold
@traitimpl IsMatrixM{Euclidean}
@traitimpl IsMatrixP{RnPoint}
@traitimpl IsMatrixTV{RnTVector}

# Functions
# ---
@doc doc"""
    distance(M,x,y)

compute the Euclidean distance $\lVert x - y\rVert$
"""
function distance(M::Euclidean,x::RnPoint{T},y::RnPoint{T})::T where {T <: AbstractFloat}
    if length(getValue(x)) > 1
        return norm( getValue(x) - getValue(y) )
    else
        return abs( getValue(x) - getValue(y) )
    end
end
@doc doc"""
    dot(M,x,ξ,ν)
Computes the Euclidean inner product of `ξ` and `ν`, i.e.
$\langle\xi,\nu\rangle = \displaystyle\sum_{k=1}^n \xi_k\nu_k$.
"""
dot(M::Euclidean,x::RnPoint{T},ξ::RnTVector{T}, ν::RnTVector{T}) where {T <: AbstractFloat} = dot( getValue(ξ) , getValue(ν) )
@doc doc"""
    exp(M,x,ξ[, t=1.0])

compute the exponential map on the [`Euclidean`](@ref) manifold `M`, i.e.
$x+t*\xi$, where the scaling parameter `t` is optional.
"""
exp(M::Euclidean,x::RnPoint{T},ξ::RnTVector{T},t::Float64=1.0) where {T <: AbstractFloat} = RnPoint(getValue(x) + t*getValue(ξ) )
@doc doc"""
    log(M,x,y)

computes the logarithmic map on the [`Euclidean`](@ref) manifold `M`, i.e. $y-x$.
"""
log(M::Euclidean,x::RnPoint{T},y::RnPoint{T})  where {T <: AbstractFloat} = RnTVector( getValue(y) - getValue(x) )
@doc doc"""
    manifoldDimension(x)

return the manifold dimension of the [`RnPoint`](@ref) `x`, i.e. $n$.
"""
manifoldDimension(x::RnPoint) = length( getValue(x) )
@doc doc"""
    manifoldDimension(M)

return the manifold dimension of the [`Euclidean`](@ref) manifold `M`,
i.e. the length of the vectors stored in `M.dimension`, i.e. $n$.
"""
manifoldDimension(M::Euclidean) = M.dimension
@doc doc"""
    parallelTransport(M,x,y,ξ)

compute the parallel transport  the [`Euclidean`](@ref) manifold `M`, which is
the identity.
"""
parallelTransport(M::Euclidean, x::RnPoint{T}, y::RnPoint{T}, ξ::RnTVector{T})  where {T <: AbstractFloat} = ξ
@doc doc"""
    randomMPoint(M[,T=Float64])

generate a random point on the [`Euclidean`](@ref) manifold `M`, where the
optional parameter determines the type of the entries of the
resulting [`RnPoint`](@ref).
"""
randomMPoint(M::Euclidean, T::DataType=Float64) = RnPoint( randn(T,M.dimension) )

doc"""
    randomTVector(M,x,:Gaussian[,σ=1.0])

generate a Gaussian random vector on the [`Euclidean`](@ref) manifold `M` with
standard deviation `σ`.
"""
randomTVector(M::Euclidean, x::RnPoint{T}, ::Val{:Gaussian}, σ::Float64=1.0) where {T} = RnTVector( σ * randn(T,M.dimension) )

@doc doc"""
    (Ξ,κ) = tangentONB(M,x,y)

compute an ONB within the tangent space $T_x\mathcal M$ at the [`MPoint`](@ref)
on the [`Euclidean`](@ref) manifold `M`, such that $\xi=\log_xy$
is the first vector and compute the eigenvalues of the curvature tensor
$R(\Xi,\dot g)\dot g$, where $g=g_{x,\xi}$ is the geodesic with $g(0)=x$,
$\dot g(0) = \xi$, i.e. $\kappa_1$ corresponding to $\Xi_1=\xi$ is zero.

# See also
 [`jacobiField`](@ref), [`adjointJacobiField`](@ref).
"""
tangentONB(M::Euclidean, x::RnPoint{T}, y::RnPoint{T}) where {T <: AbstractFloat}  = tangentONB(M,x,log(M,x,y))
tangentONB(M::Euclidean, x::RnPoint{T}, ξ::RnTVector{T}) where {T <: AbstractFloat}  =
  [ RnTVector( Matrix{T}(I,manifoldDimension(x),manifoldDimension(x))[:,i] )
    for i in 1:manifoldDimension(x)], zeros(manifoldDimension(x))

typeofTVector(::Type{RnPoint{T}}) where T = RnTVector{T}
typeofMPoint(::Type{RnTVector{T}}) where T = RnPoint{T}

@doc doc"""
    typicalDistance(M)

returns the typical distance on the [`Euclidean`](@ref) manifold `M`: $\sqrt{n}$.
"""
typicalDistance(M::Euclidean) = sqrt(M.dimension)

@doc doc"""
    validateMPoint(M,x)

Checks that a [`RnPoint`](@ref) `x` has a valid value for a point on the
[`Euclidean`](@ref) manifold `M`$=\mathbb R^n$, which is the case if the
dimensions fit.
"""
validateMPoint(M::Euclidean, x::RnPoint) = manifoldDimension(M) == manifoldDimension(x)

@doc doc"""
    validateTVector(M,x,ξ)

Checks, that the [`RnTVector`](@ref) `ξ` is a valid tangent vector in the
tangent space of the [`RnPoint`](@ref) `x` ont the [`Euclidean`](@ref)
manifold `M`, which is always the case as long as their vector dimensions agree.
"""
validateTVector(M::Euclidean,x::RnPoint,ξ::RnTVector) = length(getValue(x) )== length(getValue(ξ))

@doc doc"""
    ξ = zeroTVector(M,x)

returns a zero vector in the tangent space $T_x\mathcal M$ of the
[`RnPoint`](@ref) $x\in\mathbb R^n$ on the [`Euclidean`](@ref) manifold `M`.
"""
function zeroTVector(M::Euclidean, x::RnPoint{T}) where {T <: AbstractFloat}
    return RnTVector(  zero( getValue(x) )  )
end
#
#
# --- Display functions for the objects/types
show(io::IO, M::Euclidean) = print(io, "The $(M.name)");
show(io::IO, x::RnPoint) = print(io, "Rn($( getValue(x) ))");
show(io::IO, ξ::RnTVector) = print(io, "RnT($( getValue(ξ) ))");


text:
using COSMO, Test, LinearAlgebra, Random, SparseArrays
rng = Random.MersenneTwister(41)

# This test is precision agnostic and will be run with T precision
if !@isdefined(UnitTestFloats)
    UnitTestFloats = [Float64] #if not run in full test setup, just do it for one float type
end

@testset "Interface" begin
for T in UnitTestFloats

    @testset "Interface (T = $(T))" begin
        # set!() -  to directly set the problem data without building constraints first

        P = T[4. 1;1 2]
        q = T[1; 1.]
        A = T[1. 1;1 0; 0 1]
        l = T[1.;0;0]
        u = T[1.;0.7;0.7]
        Aa = [A; -A]
        b = [u; -l]
        convex_sets = [COSMO.Nonnegatives{T}(3); COSMO.Nonnegatives{T}(3)]

        model = COSMO.Model{T}()
        COSMO.set!(model, P, q, Aa, b, convex_sets )
        res = COSMO.optimize!(model)
        @test isapprox(res.obj_val, 1.88, atol = 1e-3)

        qf = T[1; 1; 1]
        Pf = zeros(T, 1, 1)
        Af = T[1 2; 1 2]
        bf = T[1; 2]
        model = COSMO.Model{T}()
        @test_throws DimensionMismatch COSMO.set!(model, P, qf, Aa, b, convex_sets)
        @test_throws DimensionMismatch COSMO.set!(model, Pf, q, Aa, b, convex_sets)
        @test_throws DimensionMismatch COSMO.set!(model, P, q, Af, b, convex_sets)
        @test_throws DimensionMismatch COSMO.set!(model, P, q, Aa, bf, convex_sets)

        # assemble! function test different combinations of P an q

        #P not sparse, q dense
        P = rand(rng, T, 2, 2)
        q = rand(rng, T, 2)
        A = rand(rng, T, 1, 2)
        b = one(T)
        constraint = COSMO.Constraint(A, b, COSMO.Nonnegatives)
        model = COSMO.Model{T}();
        COSMO.assemble!(model, P, q, constraint)
        @test model.p.P == P
        @test model.p.q == q
        @test model.p.A == -A
        @test model.p.b == [b]

        # sparse input data
        COSMO.empty_model!(model)
        COSMO.assemble!(model, sparse(P), sparse(q), constraint)
        @test model.p.P == P
        @test model.p.q == q

        A = rand(rng, T, 5, 1)
        b = rand(rng, T, 5)
        constraint = COSMO.Constraint(A, b, COSMO.Nonnegatives)

        # P number, q vector
        COSMO.empty_model!(model)
        COSMO.assemble!(model, one(T), [one(T)], constraint)
        @test model.p.P == hcat([one(T)])
        @test model.p.q == [one(T)]

        # P vector, q number
        COSMO.empty_model!(model)
        COSMO.assemble!(model, [one(T)], one(T), constraint)
        @test model.p.P == hcat([one(T)])
        @test model.p.q == [one(T)]

        # P number, q number
        COSMO.empty_model!(model)
        COSMO.assemble!(model, one(T), one(T), constraint, settings = COSMO.Settings{T}(verbose = true))
        @test model.p.P == hcat([one(T)])
        @test model.p.q == [one(T)]

        # P vector, q matrix
        COSMO.empty_model!(model)
        COSMO.assemble!(model, [one(T)], hcat([one(T)]), constraint)
        @test model.p.P == hcat([one(T)])
        @test model.p.q == [one(T)]

        # Test case where dimension of an A in one constraint is inconsistent with P
        A = one(T)
        b = zero(T)
        P = spdiagm(0=> ones(T, 2))
        q = rand(rng, T, 2)
        constraint = COSMO.Constraint(A, b, COSMO.Nonnegatives)
        model = COSMO.Model{T}();
        @test_throws DimensionMismatch COSMO.assemble!(model, P, q, [constraint])

        # Disallow models with PSD constraints of type BigFloat
        if T == BigFloat
            model = COSMO.Model{T}();
            P = rand(rng, T, 4, 4)
            q = rand(rng, T, 4)
            A = spdiagm(0 => ones(T, 4))
            b = zeros(T, 4)
            constraint = COSMO.Constraint(A, b, COSMO.PsdCone)
            @test_throws ArgumentError COSMO.assemble!(model, P, q, [constraint])
        end
    end
end
end


text:
using CellFreeModelGenerationKit

# path to test metabolism file -
path_to_test_metabolism_file = "./test/data/Test.vff"

# load the file -
file_buffer = read_model_document(path_to_test_metabolism_file)

# parse the metabolism section of the file_buffer -
result = parse_vff_metabolic_section(file_buffer)
if (isa(result.value,Exception) == true)
    throw(result.value)
end
df = result.value

max_stars_repo_name: varnerlab/CellFreeModelGenerationKit

max_issues_repo_name: varnerlab/CellFreeModelGenerationKit

max_forks_repo_name: varnerlab/CellFreeModelGenerationKit


text:
module GreedImages

using Images, Reel

using ...PCG.Geometry
using ...PCG.Types
using ...PCG.Universes

export Sprite, GreedImageRecorder, node_position, canvas_size, choose_sprite

const Image = Array{RGB{Normed{UInt8,8}}, 2}

struct Sprite
    color::RGBA
    _image::Image
end

struct GreedImageRecorder <: Recorder
    cell_size::Size{Int64}
    duration::Int64
    filename::String

    _frames::Vector{Image}
end

GreedImageRecorder(cell_size, duration, filename) = GreedImageRecorder(cell_size,
                                                                       duration,
                                                                       filename,
                                                                       [])

function node_position end
function canvas_size end
function choose_sprite end

function Universes.record_state!(universe::Universe, recorder::GreedImageRecorder)

    image_size = canvas_size(universe.topology, recorder)

    canvas = fill(RGBA(0, 0, 0, 1), yx(image_size))

    universe(complete_turn=false) do element
        sprite = choose_sprite(recorder, element)

        position = node_position(recorder, element.topology_index, image_size)

        image = sprite._image

        sprite_size = size(image)

        x = position.x
        y = position.y

        for i=y:(y+sprite_size[1]-1), j=x:(x+sprite_size[2]-1)
            canvas[i, j] += image[i-y+1, j-x+1]
        end

    end

    push!(recorder._frames, canvas)
end

function Universes.finish_recording!(recorder::GreedImageRecorder)
    # TODO: specify fps directly?
    frames = Frames(MIME("image/png"), fps=1.0 / (recorder.duration / 1000))

    # TODO: rewrite to dot syntax
    for frame in recorder._frames
        push!(frames, frame)
    end

    write(recorder.filename, frames)

end

end

max_stars_repo_name: Tiendil/pcg-julia

max_issues_repo_name: Tiendil/pcg-julia

max_forks_repo_name: Tiendil/pcg-julia


text:
"""
    find_most_massive_halo(filebase::String [, nfiles::Int=1])

Reads the selected file and its subfiles and returns position, virial radius
and a HaloID object that contains the subfile which contains the most massive
halo and the position in the block.

"""
function find_most_massive_halo(filebase::String, nfiles::Int=1)

    Mmax     = 0.0
    POS      = zeros(Float32, 3)
    RVIR     = 0.0
    max_file = 0
    max_id   = 0

    sub_input = select_file(filebase, 0)
    # read the info block
    info = read_info(sub_input)

    # check if the massblock is called MVIR or MTOP
    mass_block = "MVIR"
    rvir_block = "RVIR"
    if size(info[getfield.(info, :block_name) .== mass_block],1) == 0
        mass_block = "MTOP"
        rvir_block = "RTOP"
    end

    @showprogress "Reading files..." for i = 0:nfiles-1

        sub_input = select_file(filebase, i)

        M = read_subfind(sub_input, mass_block)
        max_test = findmax(M)

        if max_test[1] > Mmax

            max_file = i
            max_id   = max_test[2][1]

            # store position and virial radius of most massive halo
            POS  = read_subfind(sub_input, "GPOS")[:,max_id]
            RVIR = read_subfind(sub_input, rvir_block)[max_id]

            # store new maximum mass
            Mmax = max_test[1]
        end

    end # for

    return POS, RVIR, HaloID(max_file, max_id)
end

"""
    filter_subfind(filebase::String, blockname::String, filter_function::Function [, nfiles::Integer=1])

Selects entries in subfind block that fulfill the 'filter_funcion' requirements and
returns a 'SubfindFilter' object.

# Examples
```julia
julia> find_mass_gt_1e15(M) = ( (M > 1.e15) ? true : false )
find_mass_gt_1e15 (generic function with 1 method)
julia> filtered_subfind = filter_subfind(filebase, "MVIR", find_mass_gt_1e15)
[...]
```

"""
function filter_subfind(filebase::String, blockname::String, filter_function::Function, nfiles::Integer=1)

    # allocate empty array for SubfindFilter objects
    A = Vector{HaloID}(undef, 0)

    # loop over all files in parallel
    @threads for i = 0:nfiles-1

        # find correct input file
        sub_input = select_file(filebase, i)

        # read block
        block = read_subfind(sub_input, blockname)

        # apply filter function
        selection = filter_function.(block)

        # select matches
        correct_selection = findall(selection)

        if size(selection[selection],1) > 0

            # create array of integers for easy storing
            id_array = collect(1:size(selection,1))[correct_selection]

            # create HaloID object and push it to the storage array
            for j = 1:size(id_array,1)
                push!(A, HaloID(i, id_array[j]))
            end

        end # if entries > 0

    end # for loop

    return A
end

"""
    filter_subfind(filebase::String, blockname::String, filter_function::Function [, nfiles::Integer=1])

Selects entries in subfind block that fulfill the 'filter_funcion' requirements and
returns a 'SubfindFilter' object.

# Examples
```julia
julia> find_mass_gt_1e15(M) = ( (M > 1.e15) ? true : false )
find_mass_gt_1e15 (generic function with 1 method)
julia> filtered_subfind = filter_subfind(filebase, "MVIR", find_mass_gt_1e15)
[...]
```

"""
function filter_subfind(filebase::String, filter_function::Function, nfiles::Integer=1)

    # allocate empty array for SubfindFilter objects
    A = Vector{HaloID}(undef, 0)

    # loop over all files in parallel
    @threads for i = 0:nfiles-1

        # find correct input file
        sub_input = select_file(filebase, i)

        # apply filter function
        id_array = filter_function(sub_input)

        if size(id_array,1) > 0
            # create HaloID object and push it to the storage array
            for j = 1:size(id_array,1)
                push!(A, HaloID(i, id_array[j]))
            end

        end # if entries > 0

    end # for loop

    return A
end


text:
"""
Dykstra

"""
mutable struct Dykstra <: FOSAlgorithm
    direct::Bool
    options
end
Dykstra(;direct=false, kwargs...) = Dykstra(direct, kwargs)

struct DykstraData{T1,T2} <: FOSSolverData
    p::Array{Float64,1}
    q::Array{Float64,1}
    y::Array{Float64,1}
    S1::T1
    S2::T2
end

function init_algorithm!(alg::Dykstra, model::AbstractFOSModel)
    S1, S2, n, status_generator = get_sets_and_status(alg, model)
    data = DykstraData(zeros(n), zeros(n), Array{Float64,1}(undef, n), S1, S2)
    return data, status_generator
end

function Base.step(::Dykstra, data::DykstraData, x, i, status::AbstractStatus, longstep=nothing)
    p,q,y,S1,S2 = data.p,data.q,data.y,data.S1,data.S2

    prox!(y, S1, x .+ p)
    addprojeq(longstep, y, x .+ p)
    p .= x .+ p .- y
    prox!(x, S2, y .+ q)
    checkstatus(status, x)
    addprojineq(longstep, x, y .+ q)
    q .= y .+ q .- x
    return
end

function getsol(::Dykstra, data::DykstraData, x)
    y,S1,S2 = data.y,data.S1,data.S2
    tmp1 = similar(x) #Ok to allocate here, could overwrite p, q, if safe warmstart?
    prox!(tmp1, S1, x)
    prox!(y, S2, tmp1)
    return y #We reuse y here
end

support_longstep(::Dykstra) = true
projections_per_step(::Dykstra) = (1,1)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayTrustedRootCertificate <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id

    function ApplicationGatewayTrustedRootCertificate(;id=nothing)
        o = new()
        validate_property(ApplicationGatewayTrustedRootCertificate, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        o
    end
end # type ApplicationGatewayTrustedRootCertificate

const _property_map_ApplicationGatewayTrustedRootCertificate = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"))
const _property_types_ApplicationGatewayTrustedRootCertificate = Dict{Symbol,String}(Symbol("id")=>"String")
Base.propertynames(::Type{ ApplicationGatewayTrustedRootCertificate }) = collect(keys(_property_map_ApplicationGatewayTrustedRootCertificate))
Swagger.property_type(::Type{ ApplicationGatewayTrustedRootCertificate }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayTrustedRootCertificate[name]))}
Swagger.field_name(::Type{ ApplicationGatewayTrustedRootCertificate }, property_name::Symbol) =  _property_map_ApplicationGatewayTrustedRootCertificate[property_name]

function check_required(o::ApplicationGatewayTrustedRootCertificate)
    true
end

function validate_property(::Type{ ApplicationGatewayTrustedRootCertificate }, name::Symbol, val)
end