text:
"""
Train model through eigendecomposition
"""
function modelTrain(trainSparse::SparseMatrixCSC{Float64, Int}, trainN::Int, maxK::Int, decomp::String)
  # kernel matrix train nodes
  trainDiag = spdiagm(vec(sqrt(sum(trainSparse.^2, 2))), 0, trainN, trainN)
  trainNorm = \(trainDiag, trainSparse)
  trainDiag = []
  kernelTrain = trainNorm * trainNorm'
  kernelTrain = 0.5*(kernelTrain + kernelTrain')
  inverseD = 1./(sum(kernelTrain, 1)')
  inverseDsparse = spdiagm(vec(inverseD), 0)
  # centering matrix
  centerMatrix = eye(trainN) - ones(trainN, 1) * inverseD' / sum(inverseD)
  # eigendecomposition
  if decomp == "svd"
    # svd support for sparse matrices in version 0.4
    # use svds() with economy size
    error("modelTrain: svd not yet available.")
  elseif decomp == "eigs"
    # Lanczos method
    λ,α = eigs(inverseDsparse * centerMatrix * kernelTrain, nev=maxK+4, which=:LM)
  elseif decomp == "eig"
    # full eigendecomposition
    λ,α = eig(inverseDsparse * centerMatrix * kernelTrain)
  else
    error("modelTrain: wrong decomp argument.")
  end
  λ = real(λ)
  λ[isnan(λ)] = 0.0
  α = real(α)
  if length(λ) < maxK
    println("Less eigenvalues than requested: k = $(maxK), λ = $(length(λ))")
  end
  # sort eigenvalues in descending order and only keep top k-1
  indexλ = resize!(sortperm(λ, rev=true), maxK-1)
  λ = λ[indexλ]
  α = α[:, indexλ]
  # bias terms
  β = (1 / (sum(inverseD))) * (inverseD' * kernelTrain * α);
  # score variables
  eTrain = kernelTrain * α + repeat(β, outer = [trainN, 1])
  CB = Array(Array{Int,2}, maxK-1)
  # generate codebooks for each possible number of clusters k
  for k in 2:maxK
    CB[k-1] = codebook(eTrain[:, 1:k-1], k)
  end
  return trainNorm, λ, α, β, CB
end

"""
Validate model through Balanced Angular Fit criterion (BAF)
"""
function modelValid(validSparse::SparseMatrixCSC{Float64, Int}, validN::Int, trainNorm::SparseMatrixCSC{Float64, Int}, maxK::Int, minK::Int, α::Array{Float64, 2}, β::Array{Float64, 2}, CB::Array{Array{Int, 2}})
  # kernel matrix valid nodes
  kernelValid = validSparse*trainNorm';
  eValid = kernelValid * α + repeat(β, outer = [validN, 1])
  # out-of-sample extension
  qValid = Array(Array{Int,2}, maxK-1)
  for numk in 2:maxK
    qValid[numk-1] = membership(eValid[:,1:numk-1], CB[numk-1])
  end
  # estimate optimal number of clusters using Balanced Angular Fit (BAF)
  baf = Array(Float64, maxK-minK+1)
  for k in 1:maxK-minK+1
    baf[k] = mean(balanced_angular_fit(eValid[:, 1:minK-2+k], qValid[minK-2+k]))
  end
  valueBAF,indexBAF = findmax(baf)
  kBAF = indexBAF+minK-1
  return baf, kBAF, indexBAF, valueBAF
end

"""
Assign cluster indices to unweighted network
Blocks of nodes in parallel
"""
function modelTestU(fileCount::Int, trainNorm::SparseMatrixCSC{Float64, Int}, minN::Int, maxN::Int, α::Array{Float64,2}, β::Array{Float64,2}, CB::Array{Array{Int,2}}, kBAF::Int, maxMB::Int; limitN=-1)
  lastNeighbour = maxN
  if minN == 0
    lastNeighbour += 1
  end
  if limitN >= 0 && limitN <= maxN
    endNode = limitN
  else
    endNode = maxN
  end
  avgLength = @parallel (+) for i in 1:fileCount
    mean(fileData[i][2])
  end
  # determine number of nodes tested simultaneously
  stepSize = min(int(floor(maxMB*1000000 / (16*fileCount*nprocs()*avgLength))), int(ceil((endNode-minN+1)/nprocs())))
  println("Step size $(stepSize)")
  countN = @parallel (+) for n in colon(minN,stepSize,endNode+stepSize)
    if n > endNode
      0
    else
      if n+stepSize > endNode
        nodes = [n:endNode]
      else
        nodes = [n:n+stepSize-1]
      end
      testNodes = Array(Int,0)
      testNeighbours = Array(Int,0)
      println("Starting range $(n)")
      @inbounds for fileIndex in 1:fileCount
        foundNodes = 0
        @inbounds for index in 1:length(fileData[fileIndex][1])
          if fileData[fileIndex][1][index] in nodes
            neighbours = fileData[fileIndex][3][:, index]
            foundEnd = findfirst(neighbours, -1)
            if foundEnd != 0
              deleteat!(neighbours, foundEnd:length(neighbours))
            end
            foundNodes += 1
            append!(testNodes, fill(fileData[fileIndex][1][index], length(neighbours)))
            append!(testNeighbours, neighbours)
          end
          if foundNodes == length(nodes)
            break
          end
        end
      end
      println("Range $(n) collecting done")
      permIndices = unique(testNodes)
      testData = fill(-1, length(nodes), 2)
      if length(permIndices) != 0
        # avoid zero index in neighbours
        if minN == 0
          testNeighbours[testNeighbours .== 0] = lastNeighbour
        end
        @inbounds for i in 1:length(testNodes)
          testNodes[i] = findfirst(permIndices, testNodes[i])
        end
        # kernel matrix single test node
        testSparse = sparse(testNodes, testNeighbours, fill(1.0, length(testNodes)), length(permIndices), lastNeighbour, (w1,w2)->w1)
        kernelTest = testSparse*trainNorm';
        eTest = kernelTest * α + repeat(β, outer = [length(permIndices), 1])
        qTest = membership(eTest[:,1:kBAF-1], CB[kBAF-1])
        # store node and its cluster membership
        @inbounds for i in 1:length(permIndices)
          foundPerm = findfirst(nodes, permIndices[i])
          if foundPerm != 0
            testData[foundPerm, :] = [permIndices[i], qTest[i]]
          end
        end
      end
      println("Range $(n) storing qData")
      qData[nodes .- (minN-1), :] = testData
      length(permIndices)
    end
  end
  println("Nodes tested $(countN)")
end

"""
Assign cluster indices to weighted network
Blocks of nodes in parallel
"""
function modelTestW(fileCount::Int, trainNorm::SparseMatrixCSC{Float64, Int}, minN::Int, maxN::Int, α::Array{Float64,2}, β::Array{Float64,2}, CB::Array{Array{Int,2}}, kBAF::Int, maxMB::Int; limitN=-1)
  lastNeighbour = maxN
  if minN == 0
    lastNeighbour += 1
  end
  if limitN >= 0 && limitN <= maxN
    endNode = limitN
  else
    endNode = maxN
  end
  avgLength = @parallel (+) for i in 1:fileCount
    mean(fileData[i][2])
  end
  # determine number of nodes tested simultaneously
  stepSize = min(int(floor(maxMB*1000000*0.5 / (24*nprocs()*avgLength))), int(ceil((endNode-minN+1)/nprocs())))
  println("Step size $(stepSize)")
  countN = @parallel (+) for n in colon(minN,stepSize,endNode+stepSize)
    if n > endNode
      0
    else
      if n+stepSize > endNode
        nodes = [n:endNode]
      else
        nodes = [n:n+stepSize-1]
      end
      testNodes = Array(Int,0)
      testNeighbours = Array(Int,0)
      testWeights = Array(Float64,0)
      @inbounds for fileIndex in 1:fileCount
        foundNodes = 0
        @inbounds for index in 1:length(fileData[fileIndex][1])
          if fileData[fileIndex][1][index] in nodes
            neighbours = fileData[fileIndex][3][:, index]
            weights = fileData[fileIndex][4][:, index]
            foundEnd = findfirst(neighbours, -1)
            if foundEnd != 0
              deleteat!(neighbours, foundEnd:length(neighbours))
              deleteat!(weights, foundEnd:length(weights))
            end
            foundNodes += 1
            append!(testNodes, fill(fileData[fileIndex][1][index], length(neighbours)))
            append!(testNeighbours, neighbours)
            append!(testWeights, weights)
          end
          if foundNodes == length(nodes)
            break
          end
        end
      end
      permIndices = unique(testNodes)
      testData = fill(-1, length(nodes), 2)
      if length(permIndices) != 0
        # avoid zero index in neighbours
        if minN == 0
          testNeighbours[testNeighbours .== 0] = lastNeighbour
        end
        @inbounds for i in 1:length(testNodes)
          testNodes[i] = findfirst(permIndices, testNodes[i])
        end
        # kernel matrix single test node
        testSparse = sparse(testNodes, testNeighbours, testWeights, length(permIndices), lastNeighbour, (w1,w2)->w1)
        kernelTest = testSparse*trainNorm';
        eTest = kernelTest * α + repeat(β, outer = [length(permIndices), 1])
        qTest = membership(eTest[:,1:kBAF-1], CB[kBAF-1])
        # store node and its cluster membership
        @inbounds for i in 1:length(permIndices)
          foundPerm = findfirst(nodes, permIndices[i])
          if foundPerm != 0
            testData[foundPerm, :] = [permIndices[i], qTest[i]]
          end
        end
      end
      qData[nodes .- (minN-1), :] = testData
      length(permIndices)
    end
  end
  println("Nodes tested $(countN)")
end

"""
Assign cluster indices to network
Single node in parallel -> inefficient
"""
function modelTestSingle(fileCount::Int, trainNorm::SparseMatrixCSC{Float64, Int}, minN::Int, maxN::Int, α::Array{Float64,2}, β::Array{Float64,2}, CB::Array{Array{Int,2}}, k::Int)
  lastNeighbour = maxN
  if minN == 0
    lastNeighbour += 1
  end
  countN = @parallel (+) for n in minN:maxN
    testNeighbours = Array(Array{Int},0)
    testWeights = Array(Array{Float64},0)
    for fileIndex in 1:fileCount
      index = findfirst(fileData[fileIndex][1], n)
      if index != 0
        neighbours = fileData[fileIndex][3][:, index]
        weights = fileData[fileIndex][4][:, index]
        foundIndex = findfirst(neighbours, -1)
        if foundIndex != 0
          deleteat!(neighbours, foundIndex:length(neighbours))
          deleteat!(weights, foundIndex:length(weights))
        end
        push!(testNeighbours, neighbours)
        push!(testWeights, weights)
      end
    end
    testNeighbours = vcat(testNeighbours...)
    testWeights = vcat(testWeights...)
    if length(testNeighbours) != 0
      # avoid zero index in neighbours
      if minN == 0
        testNeighbours[testNeighbours .== 0] = lastNeighbour
      end
      # kernel matrix single test node
      testSparse = sparse(fill(1, length(testNeighbours)), testNeighbours, testWeights, 1, lastNeighbour, (w1,w2)->w1)
      kernelTest = testSparse*trainNorm';
      eTest = kernelTest * α + β
      qTest = membership(eTest[:,1:k-1], CB[k-1])
      # store node and its cluster membership
      qData[n-minN+1, :] = [n qTest]
    else
      qData[n-minN+1, :] = [-1 -1]
    end
    1
  end
  println("Nodes tested $(countN)")
end

"""
Helper functions
"""
function codebook(e::Array{Float64,2}, numk::Int)
  eBin = int(sign(e))
  eUniBin = unique(eBin, 1)
  # TODO method unique does not return index vectors (JuliaLang issue #1845)
  # current solution = nested for-loop comparison
  eCW = Array(Int, size(eBin,1))
  @inbounds for i in 1:size(eBin, 1)
    @inbounds for j in 1:size(eUniBin,1)
      if eUniBin[j,:] == eBin[i,:]
        eCW[i] = j
        break;
      end
    end
  end
  # calculate occurence of each codeword
  m = size(eUniBin, 1)
  eSizeCW = Array(Int, m)
  @inbounds for i in 1:m
    eSizeCW[i] = sum(eCW .== i)
  end
  # sort codewords by descending occurence
  indexCW = sortperm(eSizeCW, rev=true)
  if m < numk
    println("Redundant codewords for $(numk) clusters")
    return eUniBin[indexCW[1:m], :]
  end
  return eUniBin[indexCW[1:numk], :]
end

function membership(e::Array{Float64,2}, CB::Array{Int,2})
  eBin = int(sign(e))
  # Hamming distance is equal to squared Euclidean when using binary vectors
  minValues, qVectors = findmin(dist_euclidean(eBin, CB), 2)
  return int(ceil(qVectors ./ size(qVectors,1)))
end

function balanced_angular_fit(e::Array{Float64,2}, q::Array{Int,2})
  eN, eK = size(e)
  # calculate mean of each cluster
  clustermean = zeros(eK+1, eK)
  for i in 1:eK+1
    # TODO mean function along dim, JuliaLang issue #2265
    clustermean[i, :] = [mean(e[q[:] .== i, :][:,j]) for j in 1:eK]
  end
  clustermeanNorm = sqrt(sum(clustermean.^2, 2))
  # angular similarity between eigenspace projections and each cluster mean
  clusterValues = Array(Float64, eN)
  clusterIndices = Array(Int, eN)
  for i in 1:eN
    clusterSim = clustermean * e[i, :]' ./ (clustermeanNorm * norm(e[i, :]));
    maxValue, maxIndex = findmax(clusterSim, 1)
    clusterValues[i] = maxValue[1]
    clusterIndices[i] = maxIndex[1]
  end
  # cosine similarity entire cluster
  clusterFit = Array(Float64, eK+1)
  for i in 1:eK+1
    indices = findin(clusterIndices, i)
    if isempty(indices)
      clusterFit[i] = 0.0
    else
      clusterFit[i] = mean(clusterValues[indices])
    end
  end
  return clusterFit
end

function dist_euclidean(x::Array{Int,2}, y::Array{Int,2})
  xDim = size(x, 1)
  yDim = size(y, 1)
  xSum = sum(x.*x, 2)
  ySum = sum(y'.*y', 1)
  return sqrt(repeat(xSum, outer=[1,yDim]) + repeat(ySum, outer=[xDim,1]) - 2*x*y')
end


text:
function param_values()::Vector{Float64}
    p::Vector{Float64} = zeros(C.NUM)

    p[C.kscyce] = 0.003
    p[C.kdcyce] = 0.001
    p[C.kdcycee] = 0.0001
    p[C.kdcycea] = 0.03
    p[C.kasse] = 1
    p[C.kdise] = 0.02
    ## CYCA SYNTHESISp[C.DEGRADATION AND P27 BINDING/DISSOCIATION:
    p[C.kscyca] = 0.0025
    p[C.kdcyca] = 0.002
    p[C.kdcycac1] = 0.4
    p[C.kassa] = 1
    p[C.kdisa] = 0.02
    ## P27 SYNTHESIS AND DEGRADATION:
    p[C.ks27] = 0.008
    p[C.kd27] = 0.004
    p[C.kd27e] = 2
    p[C.kd27a] = 2
    ## EMI1 SYNTHESIS AND DEGRADATION:
    p[C.ksemi1] = 0.003
    p[C.kdemi1] = 0.001
    ## CDH1 REGULATION:
    p[C.Cdh1T] = 1
    p[C.kacdh1] = 0.02
    p[C.kicdh1e] = 0.07
    p[C.kicdh1a] = 0.2
    p[C.kasec] = 2
    p[C.kdiec] = 0.02
    ## SKP2 SYNTHESIS AND DEGRADATION:
    p[C.ksskp2] = 0.004
    p[C.kdskp2] = 0.002
    p[C.kdskp2c1] = 0.2
    ## CDK INHIBITOR
    p[C.Inhibitor] = 0.0

    p[C.Emi1T] = 0.0
    p[C.CycElevel] = 1.0

    return p
end

ext: jl

max_stars_repo_name: himoto/Bifurcation

max_issues_repo_name: himoto/Bifurcation

max_forks_repo_name: himoto/Bifurcation


text:
module Ex5_13c
using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

# Example 5.13c - Solution of an EMMP using EMMSYN
println("Example 5.13c")

# define s as an improper transfer function
s = rtf('s');
# enter Gu(s), Gf(s) and Mr(s)
Gu = [s/(s^2+3*s+2) 1/(s+2);
     s/(s+1) 0;
      0 1/(s+2)];
Mr = dss(eye(2));                  # enter Mr(s)
p, mu = size(Gu); mf = mu

# build the synthesis model with additive faults
sysf = fdimodset(dss(Gu), c = 1:mu, f = 1:mu);

# enter reference model
Mr = fdimodset(dss(eye(mf)), f = 1:mf)

atol = 1.e-7                # tolerance for rank tests
sdeg = -1                   # set stability degree

# solve an exact model-matching problem using EMMSYN
Q, R, info = emmsyn(sysf, Mr; atol, sdeg, minimal = false); info

# display results
println("Q = $(dss2rm(Q.sys, atol = 1.e-7))")
println("M = $(dss2rm(info.M, atol = 1.e-7))")

# check solution
G = [sysf.sys; eye(mu,mu+mf)]; F = [zeros(mf,mu) info.M*Mr.sys];
@test iszero(Q.sys*G-F; atol)

end # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
using SimpleParser
using Base.Test

# helper function for testing.
token_producer(s::String) = reverse([token for token in Lexer(s)])

plist_tokens = token_producer("{eggs = spam; foo = (bar, \"foo bar\");}")
strnum_tokens = token_producer("one 2 three 3 five")

@test pop!(plist_tokens) == Token(LBRACE, "{")
@test pop!(plist_tokens) == Token(STRING, "eggs")
@test pop!(plist_tokens) == Token(EQUAL, "=")
@test pop!(plist_tokens) == Token(STRING, "spam")
@test pop!(plist_tokens) == Token(SEMICOLON, ";")
@test pop!(plist_tokens) == Token(STRING, "foo")
@test pop!(plist_tokens) == Token(EQUAL, "=")
@test pop!(plist_tokens) == Token(LPAREN, "(")
@test pop!(plist_tokens) == Token(STRING, "bar")
@test pop!(plist_tokens) == Token(COMMA, ",")
@test pop!(plist_tokens) == Token(STRING, "foo bar")
@test pop!(plist_tokens) == Token(RPAREN, ")")
@test pop!(plist_tokens) == Token(SEMICOLON, ";")
@test pop!(plist_tokens) == Token(RBRACE, "}")
@test pop!(strnum_tokens) == Token(STRING, "one")
@test pop!(strnum_tokens) == Token(NUMBER, "2")
@test pop!(strnum_tokens) == Token(STRING, "three")
@test pop!(strnum_tokens) == Token(NUMBER, "3")
@test pop!(strnum_tokens) == Token(STRING, "five")

lexer = Lexer("{eggs = spam; foo = (bar, \"foo bar\");}")
p = Parser(lexer)
@test look_ahead_type(p) == LBRACE
next_token(p)
@test look_ahead_type(p) == STRING
next_token(p)
@test look_ahead_type(p) == EQUAL
next_token(p)
@test look_ahead_type(p) == STRING
next_token(p)
@test look_ahead_type(p) == SEMICOLON
next_token(p)
@test look_ahead_type(p) == STRING
next_token(p)
@test look_ahead_type(p) == EQUAL
next_token(p)
@test look_ahead_type(p) == LPAREN
next_token(p)
@test look_ahead_type(p) == STRING
next_token(p)
@test look_ahead_type(p) == COMMA
next_token(p)
@test look_ahead_type(p) == STRING
next_token(p)
@test look_ahead_type(p) == RPAREN
next_token(p)
@test look_ahead_type(p) == SEMICOLON
next_token(p)
@test look_ahead_type(p) == RBRACE
next_token(p)
@test look_ahead_type(p) == EOF

strange_tokens = token_producer("\"Dot(u Vector2D) float64\" \"()\$0\"")
@test pop!(strange_tokens) == Token(STRING, "Dot(u Vector2D) float64")
@test pop!(strange_tokens) == Token(STRING, "()\$0")


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

mutable struct MaintenanceRedeployStatus <: SwaggerModel
    isCustomerInitiatedMaintenanceAllowed::Any # spec type: Union{ Nothing, Bool } # spec name: isCustomerInitiatedMaintenanceAllowed
    preMaintenanceWindowStartTime::Any # spec type: Union{ Nothing, DateTime } # spec name: preMaintenanceWindowStartTime
    preMaintenanceWindowEndTime::Any # spec type: Union{ Nothing, DateTime } # spec name: preMaintenanceWindowEndTime
    maintenanceWindowStartTime::Any # spec type: Union{ Nothing, DateTime } # spec name: maintenanceWindowStartTime
    maintenanceWindowEndTime::Any # spec type: Union{ Nothing, DateTime } # spec name: maintenanceWindowEndTime
    lastOperationResultCode::Any # spec type: Union{ Nothing, String } # spec name: lastOperationResultCode
    lastOperationMessage::Any # spec type: Union{ Nothing, String } # spec name: lastOperationMessage

    function MaintenanceRedeployStatus(;isCustomerInitiatedMaintenanceAllowed=nothing, preMaintenanceWindowStartTime=nothing, preMaintenanceWindowEndTime=nothing, maintenanceWindowStartTime=nothing, maintenanceWindowEndTime=nothing, lastOperationResultCode=nothing, lastOperationMessage=nothing)
        o = new()
        validate_property(MaintenanceRedeployStatus, Symbol("isCustomerInitiatedMaintenanceAllowed"), isCustomerInitiatedMaintenanceAllowed)
        setfield!(o, Symbol("isCustomerInitiatedMaintenanceAllowed"), isCustomerInitiatedMaintenanceAllowed)
        validate_property(MaintenanceRedeployStatus, Symbol("preMaintenanceWindowStartTime"), preMaintenanceWindowStartTime)
        setfield!(o, Symbol("preMaintenanceWindowStartTime"), preMaintenanceWindowStartTime)
        validate_property(MaintenanceRedeployStatus, Symbol("preMaintenanceWindowEndTime"), preMaintenanceWindowEndTime)
        setfield!(o, Symbol("preMaintenanceWindowEndTime"), preMaintenanceWindowEndTime)
        validate_property(MaintenanceRedeployStatus, Symbol("maintenanceWindowStartTime"), maintenanceWindowStartTime)
        setfield!(o, Symbol("maintenanceWindowStartTime"), maintenanceWindowStartTime)
        validate_property(MaintenanceRedeployStatus, Symbol("maintenanceWindowEndTime"), maintenanceWindowEndTime)
        setfield!(o, Symbol("maintenanceWindowEndTime"), maintenanceWindowEndTime)
        validate_property(MaintenanceRedeployStatus, Symbol("lastOperationResultCode"), lastOperationResultCode)
        setfield!(o, Symbol("lastOperationResultCode"), lastOperationResultCode)
        validate_property(MaintenanceRedeployStatus, Symbol("lastOperationMessage"), lastOperationMessage)
        setfield!(o, Symbol("lastOperationMessage"), lastOperationMessage)
        o
    end
end # type MaintenanceRedeployStatus

const _property_map_MaintenanceRedeployStatus = Dict{Symbol,Symbol}(Symbol("isCustomerInitiatedMaintenanceAllowed")=>Symbol("isCustomerInitiatedMaintenanceAllowed"), Symbol("preMaintenanceWindowStartTime")=>Symbol("preMaintenanceWindowStartTime"), Symbol("preMaintenanceWindowEndTime")=>Symbol("preMaintenanceWindowEndTime"), Symbol("maintenanceWindowStartTime")=>Symbol("maintenanceWindowStartTime"), Symbol("maintenanceWindowEndTime")=>Symbol("maintenanceWindowEndTime"), Symbol("lastOperationResultCode")=>Symbol("lastOperationResultCode"), Symbol("lastOperationMessage")=>Symbol("lastOperationMessage"))
const _property_types_MaintenanceRedeployStatus = Dict{Symbol,String}(Symbol("isCustomerInitiatedMaintenanceAllowed")=>"Bool", Symbol("preMaintenanceWindowStartTime")=>"DateTime", Symbol("preMaintenanceWindowEndTime")=>"DateTime", Symbol("maintenanceWindowStartTime")=>"DateTime", Symbol("maintenanceWindowEndTime")=>"DateTime", Symbol("lastOperationResultCode")=>"String", Symbol("lastOperationMessage")=>"String")
Base.propertynames(::Type{ MaintenanceRedeployStatus }) = collect(keys(_property_map_MaintenanceRedeployStatus))
Swagger.property_type(::Type{ MaintenanceRedeployStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_MaintenanceRedeployStatus[name]))}
Swagger.field_name(::Type{ MaintenanceRedeployStatus }, property_name::Symbol) =  _property_map_MaintenanceRedeployStatus[property_name]

const _allowed_MaintenanceRedeployStatus_lastOperationResultCode = ["None", "RetryLater", "MaintenanceAborted", "MaintenanceCompleted"]

function check_required(o::MaintenanceRedeployStatus)
    true
end

function validate_property(::Type{ MaintenanceRedeployStatus }, name::Symbol, val)
    if name === Symbol("lastOperationResultCode")
        Swagger.validate_param(name, "MaintenanceRedeployStatus", :enum, val, _allowed_MaintenanceRedeployStatus_lastOperationResultCode)
    end
end


text:
using Test, OrderedCollections, HDF5, BenchmarkTools
include(joinpath(dirname(@__FILE__), "../../src/includeall.jl"))

time_results = false

@testset "Pseudo-Transient Relaxation" begin
    rp = joinpath(dirname(@__FILE__), "../reference/solve/pseudo_transient_relaxation.h5")
    η = h5read(rp, "eta")
    vₑ = h5read(rp, "ve")
    vₑ₁ = h5read(rp, "ve1")
    vₕ = h5read(rp, "vh")
    vₕ₁ = h5read(rp, "vh1")
    Mvₑ = h5read(rp, "Mve")
    Mvₕ = h5read(rp, "Mvh")
    MU = h5read(rp, "MU")
    S = h5read(rp, "S")
    Δ = h5read(rp, "Delta")
    G = h5read(rp, "G")
    stategrid = StateGrid(OrderedDict(:η => η))

    # Test copy from Princeton Initiative script
    myvₑ₁ = upwind_parabolic_pde(η, Mvₑ, MU, S.^2, G, vₑ, Δ)
    myvₕ₁ = upwind_parabolic_pde(η, Mvₕ, MU, S.^2, G, vₕ, Δ)
    @test myvₑ₁ ≈ vₑ₁
    @test myvₕ₁ ≈ vₕ₁

    # Test against an implementation where the operators are constructed rather than the raw LHS matrix
    N = length(η)
    DU1 = zeros(N)
    DU2 = zeros(N)
    DD1 = zeros(N)
    DD2 = zeros(N)
    D01 = zeros(N)
    D02 = zeros(N)
    Σ²0 = zeros(N)
    dX = diff(η)
    Σ²0[2:N-1] .= S[2:N-1].^2 ./ (dX[1:N-2] + dX[2:N-1]) # approx Σ² / (2 * dx): this term is the Σ²/2 coefficient
    DU1[2:N] = max.(MU[1:N - 1], 0.) ./ dX
    DU2[2:N] = Σ²0[1:N - 1] ./ dX
    DD1[1:N - 1] = max.(-MU[2:N], 0.) ./ dX
    DD2[1:N - 1] = Σ²0[2:N] ./ dX
    D01[1:N - 1] .= -DU1[2:N]
    D01[2:N] .-= DD1[1:N - 1]
    D02[1:N - 1] .= -DU2[2:N]
    D02[2:N] .-= DD2[1:N - 1]
    L1 = spdiagm(0 => D01, -1 => DD1[1:N - 1], 1 => DU1[2:N])
    L2 = spdiagm(0 => D02, -1 => DD2[1:N - 1], 1 => DU2[2:N])
    A = spdiagm(0 => (1 - Δ) .* ones(N) + Δ .* Mvₑ) - Δ .* (L1 + L2)
    directvₑ₁ = A \ ((1 - Δ) .* vₑ)
    @test myvₑ₁ ≈ directvₑ₁

    # Test new wrapper scripts utilizing DiffEqOperator for more flexible upwinding (e.g. higher-order)
    L₁ = UpwindDifference(1, 1, diff(η)[1], N)
    L₂ = CenteredDifference(2, 2, diff(η)[1], N)
    Qbc = RobinBC((0., 1., 0.), (0., 1., 0.), diff(η)[1])
    myvₑ₁ = similar(vₑ)
    myvₕ₁ = similar(vₕ)
    pseudo_transient_relaxation!(stategrid, (myvₑ₁, myvₕ₁), (vₑ, vₕ), (Mvₑ, Mvₕ),
                                 MU ./ η, (S ./ η).^2, Δ; uniform = true,
                                 L₁ = L₁, L₂ = L₂, Q = Qbc)

    @test myvₑ₁ ≈ vₑ₁
    @test myvₕ₁ ≈ vₕ₁

    L₁ = UpwindDifference(1, 1, diff(η)[1], N)
    L₂ = CenteredDifference(2, 2, diff(η)[1], N)
    Qbc = RobinBC((0., 1., 0.), (0., 1., 0.), diff(η)[1])
    myvₑ₁ = similar(vₑ)
    myvₕ₁ = similar(vₕ)
    pseudo_transient_relaxation!(stategrid, (myvₑ₁, myvₕ₁), (vₑ, vₕ), (Mvₑ, Mvₕ),
                                 MU ./ η, (S ./ η).^2, Δ; uniform = false,
                                 L₁ = L₁, L₂ = L₂, Q = Qbc)

    @test myvₑ₁ ≈ vₑ₁
    @test myvₕ₁ ≈ vₕ₁

    # Test fast default using BandedMatrices
    myvₑ₁ = similar(vₑ)
    myvₕ₁ = similar(vₕ)
    pseudo_transient_relaxation!(stategrid, (myvₑ₁, myvₕ₁), (vₑ, vₕ), (Mvₑ, Mvₕ),
                                 MU ./ η, (S ./ η).^2, Δ; uniform = true)

    @test @test_matrix_approx_eq myvₑ₁[2:end - 1] vₑ₁[2:end - 1]
    @test @test_matrix_approx_eq myvₕ₁ vₕ₁

    myvₑ₁ = similar(vₑ)
    myvₕ₁ = similar(vₕ)
    pseudo_transient_relaxation!(stategrid, (myvₑ₁, myvₕ₁), (vₑ, vₕ), (Mvₑ, Mvₕ),
                                 MU ./ η, (S ./ η).^2, Δ; uniform = false)

    @test @test_matrix_approx_eq myvₑ₁ vₑ₁
    @test @test_matrix_approx_eq myvₕ₁ vₕ₁
end

if time_results
    println("\nTiming different methods for calculating the implicit time step.\n")
    rp = joinpath(dirname(@__FILE__), "../reference/solve/pseudo_transient_relaxation.h5")
    η = h5read(rp, "eta")
    vₑ = h5read(rp, "ve")
    vₑ₁ = h5read(rp, "ve1")
    vₕ = h5read(rp, "vh")
    vₕ₁ = h5read(rp, "vh1")
    Mvₑ = h5read(rp, "Mve")
    Mvₕ = h5read(rp, "Mvh")
    MU = h5read(rp, "MU")
    S = h5read(rp, "S")
    Δ = h5read(rp, "Delta")
    G = h5read(rp, "G")

    # Test copy from Princeton Initiative script
    println("Princeton Initiative script")
    @btime begin
        upwind_parabolic_pde(η, Mvₑ, MU, S.^2, G, vₑ, Δ)
        upwind_parabolic_pde(η, Mvₕ, MU, S.^2, G, vₕ, Δ)
    end

    # Test new wrapper scripts utilizing DiffEqOperator for more flexible upwinding (e.g. higher-order)
    L₁ = UpwindDifference(1, 1, diff(η)[1], N)
    L₂ = CenteredDifference(2, 2, diff(η)[1], N)
    Qbc = RobinBC((0., 1., 0.), (0., 1., 0.), diff(η)[1])
    myvₑ₁ = similar(vₑ)
    myvₕ₁ = similar(vₕ)
    println("Using DiffEqOperator finite differencing, uniform")
    @btime begin
        pseudo_transient_relaxation!(stategrid, (myvₑ₁, myvₕ₁), (vₑ, vₕ), (Mvₑ, Mvₕ),
                                     MU ./ η, (S ./ η).^2, Δ; uniform = true,
                                     L₁ = L₁, L₂ = L₂, Q = Qbc)
    end

    L₁ = UpwindDifference(1, 1, diff(η)[1], N)
    L₂ = CenteredDifference(2, 2, diff(η)[1], N)
    Qbc = RobinBC((0., 1., 0.), (0., 1., 0.), diff(η)[1])
    myvₑ₁ = similar(vₑ)
    myvₕ₁ = similar(vₕ)
    println("Using DiffEqOperator finite differencing, non-uniform")
    @btime begin
        pseudo_transient_relaxation!(stategrid, (myvₑ₁, myvₕ₁), (vₑ, vₕ), (Mvₑ, Mvₕ),
                                     MU ./ η, (S ./ η).^2, Δ; uniform = false,
                                     L₁ = L₁, L₂ = L₂, Q = Qbc)
    end

    # Test fast default using BandedMatrices
    myvₑ₁ = similar(vₑ)
    myvₕ₁ = similar(vₕ)
    println("Using BandedMatrices, uniform")
    @btime begin
        pseudo_transient_relaxation!(stategrid, (myvₑ₁, myvₕ₁), (vₑ, vₕ), (Mvₑ, Mvₕ),
                                     MU ./ η, (S ./ η).^2, Δ; uniform = true)
    end

    myvₑ₁ = similar(vₑ)
    myvₕ₁ = similar(vₕ)
    println("Using BandedMatrices, non-uniform")
    @btime begin
        pseudo_transient_relaxation!(stategrid, (myvₑ₁, myvₕ₁), (vₑ, vₕ), (Mvₑ, Mvₕ),
                                     MU ./ η, (S ./ η).^2, Δ; uniform = false)
    end
end

nothing


text:
"""
    go2lab(rng, subspace, Toy_data ; repeat = 3 ,sigma = 1)

Simulated the evaluation of constructs in the lab using the given Toy data set.
This function wraps the evaluation of multiple constructs. The mean and standard deviation are returned:

[`_eval_con`](@ref)
"""

function go2lab(rng, subspace, Toy_data; repeat=3 ,sigma=1)
    μ = Vector{Float64}(undef,length(subspace))
    σ = Vector{Float64}(undef,length(subspace))
    for (i,construct) in enumerate(subspace)
         results =  _eval_con(rng,construct,Toy_data,sigma,repeat)
         μ[i] = Statistics.mean(results)
         σ[i] = Statistics.std(results)
    end
    return μ,σ
end

"""
    ground_truth(construct, Toy_data)

Calculates the true value of a construct given the Toy Data set.
"""
ground_truth(construct,Toy_data) = sum([Toy_data[mod] for mod in construct])

"""
_eval_con(rng,construct,Toy_data,sigma,repeat)

Internal function to evaluated a single constructed.
Sigma sets the standard deviation of the normal distribution that sets that generates the activities.
Repeats is the number of time the same constructed is evaluated.

[`_eval_con`](@ref)
"""
function _eval_con(rng,construct,Toy_data,sigma,repeat)
  mean = ground_truth(construct,Toy_data)
  distr = Normal(mean,sigma)
  return  [rand(rng,distr) for rep in 1:repeat]
end

max_stars_repo_name: SamDier/BOMoD

max_issues_repo_name: SamDier/BOMoD

max_forks_repo_name: SamDier/BOMoD


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
contract_left.jl defines methods for the contract_left(x, y) function
"""

# --- Exports

export contract_left
import Base.:(<)

# --- Operator aliases

<(M::AbstractMultivector, N::AbstractMultivector) = contract_left(M, N)

<(M::AbstractMultivector, x::Real) = contract_left(M, x)
<(x::Real, M::AbstractMultivector) = contract_left(x, M)

<(M::AbstractMultivector, v::Vector{<:Real}) = contract_left(M, v)
<(v::Vector{<:Real}, M::AbstractMultivector) = contract_left(v, M)

# --- Method definitions

"""
    contract_left(M::AbstractMultivector,
                  N::AbstractMultivector)::AbstractMultivector
    <(M::AbstractMultivector, N::AbstractMultivector)::AbstractMultivector

Compute the left contraction of `M` with `N`.
"""
function contract_left end

# ------ Specializations involving an AbstractMultivector instance

# M::AbstractMultivector, B::AbstractBlade
# B::AbstractBlade, M::AbstractMultivector
contract_left(M::AbstractMultivector, B::AbstractBlade) =
    Multivector(map(C -> contract_left(C, B), blades(M)))

contract_left(B::AbstractBlade, M::AbstractMultivector) =
    Multivector(map(C -> contract_left(B, C), blades(M)))

# M::AbstractMultivector, B::AbstractScalar
# B::AbstractScalar, M::AbstractMultivector
contract_left(M::AbstractMultivector, B::AbstractScalar) =
    length(M[0]) > 0 ?
        Scalar(value(M[0][1]) * B) :
        zero(Scalar{typeof(norm(M))})

contract_left(B::AbstractScalar, M::AbstractMultivector) = B * M

# M::AbstractMultivector, B::Zero
# B::Zero, M::AbstractMultivector
contract_left(B::Zero, M::AbstractMultivector) = zero(M)
contract_left(M::AbstractMultivector, B::Zero) = zero(M)

# M::AbstractMultivector, x::Real
# x::Real, M::AbstractMultivector
contract_left(M::AbstractMultivector, x::Real) =
    length(M[0]) > 0 ?
        Scalar(value(M[0][1]) * x) :
        zero(Scalar{typeof(norm(M))})

contract_left(x::Real, M::AbstractMultivector) = x * M

# ------ Specializations involving an AbstractBlade instance

# B::AbstractBlade, C::AbstractScalar
# B::AbstractScalar, C::AbstractBlade
contract_left(B::AbstractBlade, C::AbstractScalar) = zero(B)

contract_left(B::AbstractScalar, C::AbstractBlade) =
    Blade(C, volume=value(B) * volume(C))

# B::AbstractBlade, C::Zero
# B::Zero, C::AbstractBlade
contract_left(B::AbstractBlade, C::Zero) = zero(B)
contract_left(B::Zero, C::Blade) = zero(B)

# ------ Specializations involving a Blade instance

# B::Blade, C::Blade
function contract_left(B::Blade, C::Blade)
    # --- Check arguments

    assert_dim_equal(B, C)

    # --- Handle edge cases

    # (B ⋅ C) = 0 if grade(B) > grade(C)
    if grade(B) > grade(C)
        return zero(B)
    end

    # --- Compute (B ⋅ C) = project(B, C) * C
    #     = volume(C) dual(project(B, C), C) * I_C
    #     = (-1)^((grade(C) * grade(C) - 1) / 2) volume(C) project(B, C) / I_C
    #     = (-1)^((grade(C) * grade(C) - 1) / 2) volume(C)
    #       dual(project(B, C), C)
    #
    #     where I_C is the unit blade for the subspace represented by blade `C`
    #     that has the same orientation as basis(C).

    # Compute project(B, C) = (B ⋅ C) / C
    projection = project(B, C)

    # Compute
    #   (-1)^((grade(C) * grade(C) - 1) / 2) volume(C) dual(project(B, C), C)
    mod(grade(C), 4) < 2 ?
        volume(C) * dual(projection, C) :
       -volume(C) * dual(projection, C)
end

# B::Blade, C::Pseudoscalar
# B::Pseudoscalar, C::Blade
function contract_left(B::Blade, C::Pseudoscalar)
    assert_dim_equal(B, C)
    mod(grade(C), 4) < 2 ?
        dual(B, C) * volume(C) :
       -dual(B, C) * volume(C)
end

function contract_left(B::Pseudoscalar, C::Blade)
    assert_dim_equal(B, C)
    zero(B)
end

# B::Blade, v::Vector
# v::Vector, B::Blade
function contract_left(B::Blade, v::Vector{<:Real})
    assert_dim_equal(v, B)
    grade(B) > 1 ? zero(B) : Scalar(volume(B) * basis(B) ⋅ v)
end

function contract_left(v::Vector{<:Real}, B::Blade)
    # --- Check arguments

    assert_dim_equal(v, B)

    # --- Compute (v ⋅ B) = project(v, B) * B
    #     = volume(B) dual(project(v, B), B) * I_B
    #     = (-1)^((grade(B) * grade(B) - 1) / 2) volume(B) project(v, B) / I_B
    #     = (-1)^((grade(B) * grade(B) - 1) / 2) volume(B)
    #       dual(project(v, B), B)
    #
    #     where I_B is the unit blade for the subspace represented by blade `B`
    #     that has the same orientation as basis(B).

    # Compute project(v, B) = (v ⋅ B) / B
    projection = project(v, B)

    # Compute
    #   (-1)^((grade(B) * grade(B) - 1) / 2) volume(B) dual(project(v, B), B)
    mod(grade(B), 4) < 2 ?
        volume(B) * dual(projection, B) :
       -volume(B) * dual(projection, B)
end

# B::Blade, x::Real
# x::Real, B::Blade
contract_left(B::Blade, x::Real) = zero(B)
contract_left(x::Real, B::Blade) = x * B

# ------ Specializations involving a Pseudoscalar instance

# B::Pseudoscalar, C::Pseudoscalar
function contract_left(B::Pseudoscalar, C::Pseudoscalar)
    assert_dim_equal(B, C)
    mod(grade(B), 4) < 2 ?
        Scalar(value(B) * value(C)) :
        Scalar(-value(B) * value(C))
end

# B::Pseudoscalar, C::AbstractScalar
# B::AbstractScalar, B::Pseudoscalar
contract_left(B::Pseudoscalar, C::AbstractScalar) = zero(B)
contract_left(B::AbstractScalar, C::Pseudoscalar) =
    Pseudoscalar(C, value=value(B) * value(C))

# B::Pseudoscalar, C::Zero
# B::Zero, C::Pseudoscalar
contract_left(B::Pseudoscalar, C::Zero) = zero(B)
contract_left(B::Zero, C::Pseudoscalar) = zero(B)

# B::Pseudoscalar, v::Vector
# v::Vector, B::Pseudoscalar
function contract_left(B::Pseudoscalar, v::Vector{<:Real})
    assert_dim_equal(B, v)
    zero(B)
end

function contract_left(v::Vector{<:Real}, B::Pseudoscalar)
    assert_dim_equal(B, v)

    mod(grade(B), 4) < 2 ?
        volume(B) * dual(Blade(v)) :
       -volume(B) * dual(Blade(v))
end

# ------ Specializations involving an AbstractScalar instance

# B::AbstractScalar, C::One
# B::One, C::AbstractScalar
contract_left(B::AbstractScalar, C::One) = B
contract_left(B::One, C::AbstractScalar) = C

# B::AbstractScalar, C::Zero
# B::Zero, C::AbstractScalar
contract_left(B::AbstractScalar, C::Zero) = C
contract_left(B::Zero, C::AbstractScalar) = B

# B::AbstractScalar, x::Real
# x::Real, B::AbstractScalar
contract_left(B::AbstractScalar, x::Real) = B * x
contract_left(x::Real, B::AbstractScalar) = x * B

# B::AbstractScalar, v::Vector
# v::Vector, B::AbstractScalar
contract_left(B::AbstractScalar, v::Vector{<:Real}) = Blade(value(B) * v)
contract_left(v::Vector{<:Real}, B::AbstractScalar) = zero(B)

# ------ Specializations involving a Scalar instance

# B::Scalar, C::Scalar
contract_left(B::Scalar, C::Scalar) = B * C

# ------ Specializations involving a One instance

# B::One, C::One
contract_left(B::One, C::One) = B

# B::One, C::Zero
# B::Zero, C::One
contract_left(B::One, C::Zero) = C
contract_left(B::Zero, C::One) = B

# ------ Specializations involving a Zero instance

# B::Zero, C::Zero
contract_left(B::Zero, C::Zero) = B


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

using Markdown
using InteractiveUtils

# ╔═╡ de7a5e8a-1ee4-11eb-2bbf-0f0eb3e489ec
using Pkg, DrWatson

# ╔═╡ e4eeb13a-1ee4-11eb-0113-91f9a9c3b659
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Plots
end

# ╔═╡ ac3ceabe-1ee4-11eb-3935-313e96aafd62
md"## Listing 2.9"

# ╔═╡ c4bef1a4-1ee9-11eb-0b58-8be033f6fbad
occupancyAnalytic(n,r) =  sum([(-1)^k*binomial(n,k)*(1 - k/n)^r for k in 0:n])

# ╔═╡ 8ce90c92-1eee-11eb-21f7-f78acba42f4a
function occupancyMC(n,r,N)
    fullCount = 0
    for _ in 1:N
        envelopes = zeros(Int,n)
        for k in 1:r
            target = rand(1:n)
            envelopes[target] += 1
        end
        numFilled = sum(envelopes .> 0)
        if numFilled == n
            fullCount += 1
        end
    end
    return fullCount/N
end

# ╔═╡ 8ce947c0-1eee-11eb-0f4e-49531a7b01c6
begin
	max_n, N, Kvals = 100, 10^3, [2,3,4]

	analytic = [[occupancyAnalytic(big(n),big(k*n)) for n in 1:max_n] for k in Kvals]
	monteCarlo = [[occupancyMC(n,k*n,N) for n in 1:max_n] for k in Kvals]

	plot(1:max_n, analytic, c=[:blue :red :green],
		label=["K=2" "K=3" "K=4"])
	scatter!(1:max_n, monteCarlo, mc=:black, shape=:+,
		label="", xlims=(0,max_n),ylims=(0,1),
		xlabel="n Envelopes", ylabel="Probability", legend=:topright)
end

# ╔═╡ ea086874-1eeb-11eb-2cb0-4186697d4294
md"## End of listing 2.9"

# ╔═╡ Cell order:
# ╟─ac3ceabe-1ee4-11eb-3935-313e96aafd62
# ╠═de7a5e8a-1ee4-11eb-2bbf-0f0eb3e489ec
# ╠═e4eeb13a-1ee4-11eb-0113-91f9a9c3b659
# ╠═c4bef1a4-1ee9-11eb-0b58-8be033f6fbad
# ╠═8ce90c92-1eee-11eb-21f7-f78acba42f4a
# ╠═8ce947c0-1eee-11eb-0f4e-49531a7b01c6
# ╟─ea086874-1eeb-11eb-2cb0-4186697d4294


text:
#= ParseState =#
# Makes passing everything around a bit easier.
struct GCommand
    func
    args::Vector
end

mutable struct ParseState
    basename::AbstractString
    cframe::Int
    commands::Vector{GCommand}
    coords::Vector{Matrix{Float64}}
    display::IBuffer
    nframes::Int
    symtab::Dict{Symbol, Any}
    transmat::Matrix{Float64}
end
ParseState() = ParseState("", 1, Vector{GCommand}(), [eye(4)], IBuffer(), 1, Dict(), eye(4))

modifycoord!(ps::ParseState, trans::Matrix{Float64}) = push!(ps.coords, pop!(ps.coords) * trans)

withpm!(ps::ParseState, func, args...) = begin
    tmpp = Edges()
    func(tmpp, args...)
    drawpm!(tmpp * ps.coords[end], ps.display, VIEWVEC, CAMBIENT, LIGHTS, REFLECTION)
end

withem!(ps::ParseState, func, args...) = begin
    tmpe = Edges()
    func(tmpe, args...)
    drawem!(tmpe * ps.coords[end], ps.display, [255, 255, 255])
end

display(ps::ParseState, cmd::Cmd) = begin
    try
        open(io::IO -> dump_ppm_p6(ps.display, io), cmd, "w")
    catch e
        e isa Base.UVError || rethrow()
    end
end

#= Knobs =#
struct Knob
    sframe::Int
    eframe::Int
    svalue::Float64
    evalue::Float64
    frdata::Vector{Float64}
end
Knob(sf::Int, ef::Int, sv::Float64, ev::Float64, fr::Int) = Knob(sf, ef, sv, ev, zeros(fr))

getindex(k::Knob, i::Union{Int,UnitRange}) = k.frdata[i]

# extensions to JuliaParser.Lexer
function peekchars(io::IO, count::Int)
    mark(io)
    ret = String(read(io, count))
    reset(io)
    return ret
end
peekchars(ts::TokenStream, count::Int) = peekchars(ts.io, count)

peekcomment(ts::TokenStream, cm::AbstractString) = peekchars(ts, length(cm)) == cm

function skipcomment(ts::TokenStream, cm::AbstractString)
    @assert peekcomment(ts, cm)
    Lexer.skip_to_eol(ts)
end

function skipws_and_comments(ts::TokenStream, cm::AbstractString)
    while ! eof(ts)
        Lexer.skipws(ts, true)
        peekcomment(ts, cm) || break
        skipcomment(ts, cm)
    end
end

function require_token(ts::TokenStream, nextt=Lexer.next_token)
    if ts.putback !== nothing
        t = ts.putback
    elseif ts.lasttoken !== nothing
        t = ts.lasttoken
    else
        t = nextt(ts)
    end
    eof(t) && throw(Incomplete(:other, diag(here(ts), "incomplete: premature end of input")))
    while Lexer.isnewline(t)
        Lexer.take_token(ts)
        t = nextt(ts)
    end
    eof(t) && throw(Incomplete(:other, diag(here(ts), "incomplete: premature end of input")))
    ts.putback === nothing && Lexer.set_token!(ts, t)
    return t
end

# entry method
function renderfile(contents::Vector{String}, lang::AbstractString="mdl")
    submod = getfield(Render, Symbol(lang))
    parser = getfield(submod, Symbol(lang, :_parser))
    execl  = getfield(submod, Symbol(lang, :_execute))

    ps = ParseState()
    for line in contents
        ts = TokenStream{Lexer.SourceLocToken}(line)
        parser(ts, ps)
    end
    execl(ps)
    ps
end
renderfile(io::IO, lang::AbstractString="mdl") = renderfile(readlines(io), lang)
renderfile(file::AbstractString, lang::AbstractString="mdl") = open((io) -> renderfile(io, lang), file)

max_stars_repo_name: AnAverageHuman/jlgraphics

max_issues_repo_name: AnAverageHuman/jlgraphics

max_forks_repo_name: AnAverageHuman/jlgraphics


text:
@testset "59.spiral-matrix-ii.jl" begin
    @test generate_matrix(3) == [[1, 2, 3], [8, 9, 4], [7, 6, 5]]
    @test generate_matrix(1) == [[1]]
    @test generate_matrix(5) == [
        [1, 2, 3, 4, 5],
        [16, 17, 18, 19, 6],
        [15, 24, 25, 20, 7],
        [14, 23, 22, 21, 8],
        [13, 12, 11, 10, 9],
    ]
end


text:
include("./params.jl")

function pprocess_gamma(lambda,k)
    # Example usage:
    # lambda = 20.0; # Spike-rate (Hz).
    # k = 2.3;       # Shape parameter.
    # x, rs = pprocess_gamma(lambda, k)
    #
    # Input argument
    # y:    Spike-rate.
    # k:    The shape parameter of the ISI distribution.
    #       * To use small k << 1, use small sampling time dt.
    #       * Maximum k allowed is ~10.
    #
    # Output argument
    # x:    Spike time.
    # rs:   Intensity function.

    ################# WRITE SOMETHING #################

    return x, rs
end

max_stars_repo_name: ganow/julia-point-process

max_issues_repo_name: ganow/julia-point-process

max_forks_repo_name: ganow/julia-point-process


text:
# This file contains functions for the Profit Modified Directional Distance Function Efficiency DEA model
"""
    ProfitModifiedDDFDEAModel
An data structure representing a profit Modified Directional Distance Function DEA model.
"""
struct ProfitModifiedDDFDEAModel <: AbstractProfitDEAModel
    n::Int64
    m::Int64
    s::Int64
    Gx::Symbol
    Gy::Symbol
    monetary::Bool
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    lambda::SparseMatrixCSC{Float64, Int64}
    techeff::Vector
    alloceff::Vector
    normalization::Vector
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    deaprofitmddf(X, Y, W, P; Gx, Gy)
Compute profit efficiency using Modified DDF data envelopment analysis model for
inputs `X`, outputs `Y`, price of inputs `W`, and price of outputs `P`.

# Direction specification:

The directions `Gx` and `Gy` can be one of the following symbols.
- `:Observed`: use observed values.

Alternatively, a vector or matrix with the desired directions can be supplied.

# Optional Arguments
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [2; 4; 8; 12; 6; 14; 14; 9.412];

julia> Y = [1; 5; 8; 9; 3; 7; 9; 2.353];

julia> W = [1; 1; 1; 1; 1; 1; 1; 1];

julia> P = [2; 2; 2; 2; 2; 2; 2; 2];

julia> deaprofitmddf(X, Y, W, P, Gx = :Observed, Gy = :Observed)
Profit Modified DDF DEA Model
DMUs = 8; Inputs = 1; Outputs = 1
Returns to Scale = VRS
Gx = Observed; Gy = Observed
─────────────────────────────────────
       Profit   Technical  Allocative
─────────────────────────────────────
1  4.0         7.623e-7    4.0
2  0.5         8.28765e-8  0.5
3  5.68126e-8  4.78858e-8  8.92688e-9
4  0.166667    1.08205e-7  0.166667
5  1.33333     1.16667     0.166667
6  0.571429    0.571429    3.2471e-9
7  0.285714    0.142857    0.142857
8  2.69996     2.54994     0.150021
─────────────────────────────────────
```
"""
function deaprofitmddf(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector},
    W::Union{Matrix,Vector}, P::Union{Matrix,Vector};
    Gx::Union{Symbol,Matrix,Vector}, Gy::Union{Symbol,Matrix,Vector},
    monetary::Bool = false,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::ProfitModifiedDDFDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    nw, mw = size(W, 1), size(W, 2)
    np, sp = size(P, 1), size(P, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nw != nx
        throw(DimensionMismatch("number of rows in W and X ($nw, $nx) are not equal"));
    end
    if np != ny
        throw(DimensionMismatch("number of rows in P and Y ($np, $ny) are not equal"));
    end
    if mw != m
        throw(DimensionMismatch("number of columns in W and X ($mw, $m) are not equal"));
    end
    if sp != s
        throw(DimensionMismatch("number of columns in P and Y ($sp, $s) are not equal"));
    end

    # Build or get user directions
    if typeof(Gx) == Symbol
        Gxsym = Gx

        if Gx == :Observed
            Gx = X
        else
            throw(ArgumentError("Invalid `Gx`"));
        end

    else
        Gxsym = :Custom
    end

    if typeof(Gy) == Symbol
        Gysym = Gy

        if Gy == :Observed
            Gy = Y
        else
            throw(ArgumentError("Invalid `Gy`"));
        end

    else
        Gysym = :Custom
    end

    if (size(Gx, 1) != size(X, 1)) | (size(Gx, 2) != size(X, 2))
        throw(DimensionMismatch("size of Gx and X ($(size(Gx)), $(size(X))) are not equal"));
    end
    if (size(Gy, 1) != size(Y, 1)) | (size(Gy, 2) != size(Y, 2))
        throw(DimensionMismatch("size of Gy and Y ($(size(Gy)), $(size(Y))) are not equal"));
    end

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:LP)
    end

    # Get maximum profit targets and lambdas
    n = nx

    Xtarget, Ytarget, plambdaeff = deamaxprofit(X, Y, W, P, optimizer = optimizer)

    # Profit, technical and allocative efficiency
    maxprofit = sum(P .* Ytarget, dims = 2) .- sum(W .* Xtarget, dims = 2)
    obsprofit = sum(P .* Y, dims = 2) .- sum(W .* X, dims = 2)
    pefficiency = vec(maxprofit .- obsprofit)

    normalization = zeros(n)
    for i in 1:n
        if obsprofit[i] >= 0
            normalization[i] = sum(W[i,:] .* Gx[i,:], dims = 2)[1]
        else
            normalization[i] = sum(P[i,:] .* Gy[i,:], dims = 2)[1]
        end
    end
    normalization = vec(normalization)

    techefficiency = efficiency(deamddf(X, Y, Gx = Gx, Gy = Gy, rts = :VRS, slack = false, optimizer = optimizer))

    if monetary
        techefficiency = techefficiency .* normalization
    else
        pefficiency = pefficiency ./ normalization
    end

    allocefficiency = pefficiency - techefficiency

    return ProfitModifiedDDFDEAModel(n, m, s, Gxsym, Gysym, monetary, names, pefficiency, plambdaeff, techefficiency, allocefficiency, normalization, Xtarget, Ytarget)

end

function Base.show(io::IO, x::ProfitModifiedDDFDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    dmunames = names(x)

    eff = efficiency(x)
    techeff = efficiency(x, :Technical)
    alloceff = efficiency(x, :Allocative)

    if !compact
        print(io, "Profit Modified DDF DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Returns to Scale = VRS")
        print(io, "\n")
        print(io, "Gx = ", string(x.Gx), "; Gy = ", string(x.Gy))
        print(io, "\n")
        show(io, CoefTable(hcat(eff, techeff, alloceff), ["Profit", "Technical", "Allocative"], dmunames))
    end

end


text:
import Pkg; Pkg.add("DescriptorSystems")
using Documenter, DescriptorSystems, FaultDetectionTools, DescriptorSystems
DocMeta.setdocmeta!(FaultDetectionTools, :DocTestSetup, :(using FaultDetectionTools); recursive=true)

makedocs(
  modules  = [FaultDetectionTools],
  sitename = "FaultDetectionTools.jl",
  authors  = "Andreas Varga",
  format   = Documenter.HTML(prettyurls = false),
  pages    = [
    "Home"   => "index.md",
    "Tutorials"   => [
      "FDDbasics.md",
      "SynthesisParadigms.md"
      ],
    "Fault Detection" => [
        "FDIObjects.md",
        "FDIanalysis.md",
        "FDIsynthesis.md",
        "FDIperformance.md"],
    "Model Detection" => [
        "MDObjects.md",
        "MDanalysis.md",
        "MDsynthesis.md",
        "MDperformance.md"
        ],
     "Utilities" => [
        "FDIutils.md"
     ],
     "Index" => "makeindex.md"
  ]
)

deploydocs(
  repo = "github.com/andreasvarga/FaultDetectionTools.jl.git",
  target = "build",
  devbranch = "master"
)


text:
immutable Tag
    rep::AbstractString
end

type TaggedValue
    tag::AbstractString
    value::Any
end

import Base.==, Base.hash, Base.string

function ==(tv::TaggedValue, other::TaggedValue)
    tv.value == other.value && tv.tag == other.tag
end

const hashttaggedvalue_seed = UInt === UInt64 ? 0x8ee14727cfcae1ce : 0x1f487ece

function hash(tv::TaggedValue, h::UInt)
    h = hash(hashttaggedvalue_seed, h)
    h $= hash(tv.tag)
    h $= hash(tv.value)
    return h
end

function string(tv::TaggedValue)
    "$(tv.tag): $(tv.value)"
end


text:
# Part 1
diagnostic = map(x -> parse.(Bool, split(x, "")), readlines("input.txt"))
γ = sum(diagnostic) .>= length(diagnostic)/2
ϵ = .!γ
bit2decimal(x) = sum([2^(i-1) for i in eachindex(x) if Bool(x[end-(i-1)])])
@info "Power consumption = $(prod(bit2decimal, [γ, ϵ]))" bit2decimal(γ) bit2decimal(ϵ)

# Part 2
bitfilter(y, b, i) = filter(y) do x
    target = sum(z -> z[i], y) >= length(y)/2
    length(y) > 1 ? b ? x[i] == target : x[i] != target : true
end

function findrate(z, b)
    inds = (reverse ∘ eachindex ∘ first)(z)
    funs = [x -> bitfilter(x, b, i) for i in inds]
    (only ∘ reduce(∘, funs))(z)
end

O₂, CO₂ = map(x -> findrate(diagnostic, x), [true, false])
@info "Life support rating = $(prod(bit2decimal, [O₂, CO₂]))" bit2decimal(O₂) bit2decimal(CO₂)


text:
#=
Simple Bayesian estimation example on 1D terrain-relative localization.
=#

using DataInterpolations
using Distributions  # measurement models
using ImageFiltering # conv (same), see https://discourse.julialang.org/t/convolution-conv-with-same-size-output/38260/10
using Plots          # uhh, let me be clear
gr()

## 0. Terrain gen
function randu(low::Float64, high::Float64)::Float64
    # rand will return an array, so need to index into it
    return low + (high-low)*rand(Float64,1)[1]
end

function mpe(h)
    # midpoint elevation algorithm
    he = zeros(2*length(h)-1)
    he[1:2:end] = h   # original values
    # midpoints
    for i in 2:2:(length(he)-1)
        (low, high) = minmax(he[i-1],he[i+1])
        he[i] = randu(low,high)
    end
    return he
end

terrain = 4*randn(10)
for i in 1:8
  terrain= mpe(terrain)
end

x = range(-20, 20,length=length(terrain))
plot(x,terrain,label="height")

# interpolated version
terrain_ = LinearInterpolation(terrain, x)

## 1. Prior and motion model
function convs(signal, kernel)
    # convenience wrapper to get same length on output signal
    return imfilter(signal, reflect(centered(kernel)), Fill(0) )
end

function propagate(prior, model)
    # propagate the system prior through the motion model
    convs(prior, model)
end

function neval(dist, x)
    d = pdf.(dist, x)
    d./(diff(x)[1]*sum(d))
end

motion_model =  Normal(1,0.1)
motion_kernel = neval(motion_model, x)
plot(x, motion_kernel, label="kernel")
# prior distribution on position
prior = neval(Normal(-5,0.1), x)
dx = diff(x)[1]
plot!(x, prior./(dx*sum(prior)),label="prior")

post = propagate(prior, motion_kernel)
plot!(x, post./(dx*sum(post)), label="post")

for i=1:3
    post = propagate(post, motion_kernel)
    plot!(x, post./(dx*sum(post)), label="post")
end

## 2. measurement model
function measurementLikelihood( xt,
                                terrain::LinearInterpolation)
  #
  y = terrain(xt)
    # y = terrain[argmin(abs(xt .- x))]
  alti = Normal(y, 1) # meas model
  p = pdf.(alti, terrain.u)
  return p ./ sum(p)
end

terrain_(-5.1) # if terrain is DataInterpolations.jl object, you can do terrain(xt)

pxz = measurementLikelihood( -5, terrain_ )
  # # terrain height at the true position
  # y = terrain[argmax(prior)] # true measurement
  # # measurement model
  # alti = Normal(y, 1) # measurement + uncertaint
  # # compute p(x|z) and plot normalized
  # pxz = pdf.(alti, terrain)
  # pxz./(dx*sum(pxz))

plot(x, pxz, label="likelihood")

## 3. sim
post = neval(Normal(-5,0.1),x) # est position vs. xt
xt= -5.0
for i=1:10
    # true val (true position), like in dirac postion point
    xt= xt+1.0

    # propagate
    post = propagate(post, motion_kernel)

    # measurement
    likelihood = measurementLikelihood(xt, terrain_)
      # y = terrain[argmin(abs(xt-x))]
      # alti = Normal(y, 1) # meas model
      # likelihood = pdf.(alti, terrain)

    # update
    post = post.*likelihood
end

## implements a scalar measurement factor (e.g. elevation against known scalar map)
# addFactor!(fg, [:x1, :topography], ScalarFieldFactor(y))
struct ScalarFieldFactor{T<:SamplableBelief} <: IIF.AbstractRelative__
  altimeter::T
end

function IIF.getSample(s::CalcFactor{<:ScalarFieldFactor}, N::Int=1)
    # 1. generate measurement model -
    #   1a. (high-level): we get a single value (e.g. 20m elevation) - Normal(y, sigma)
    #   1b. (low-level):  we get a vector-valued measurement (e.g. back-scatter strength) - Normal(argmax(y),sigma) [classic] or BSS(y) [full measurement]
    #   NOTE: 1b is not the n-dimensional vector field scenario; it is just a more complicated (un-processed) measurement
    # 2. generate likelihood (eval measurement model over field)
    # 3. draw samples from likelihood

  return (reshape(rand(s.altimeter,N),1,N),)
end

# keep track of IIF #1051
function (s::CalcFactor{<:ScalarFieldFactor})(elevMeasurement, pose_i)
  # assume pose_i is from a Point2
  # magically need access to global constant fields already loaded
  # addConstant!(fg, :topography, ScalarField(....,manifoldDomain))

  # NOTE: topoMap must support 2d interpolation to allow for `topoMap(pose)` querying
  topoMap = s.constants[:topography] # still early estimate of API

  # terrain_ = LinearInterpolation(terrain, x)
  # x,y is test point on prior data
  residual = elevMeasurement - topoMap(pose_i) # topoMap == terrain_

  return residual
end

##

#=
# Below: discussion on syntax for different use cases

## for 1a.
# user only identifies the altimeter reading (map magically happens)
# pxz = measurementLikelihood();
# likelihoodSampler = bss(x, pxz)
# addFactor!(fg, :x6, Point1Prior(bss))

# semantic question on addConstant!(fg,...) vs. addField!(fg, ...)
#  con field is overloaded term
#  con constant doesnt at first caputer

fg = initfg()
# addField!()
addConstant!(fg, :m, ScalarField(......????)) # load DEM here # on what manifold?

addVariable!(fg, :x0, ContinuousEuclid{1})
addVariable!(fg, :x1, ContinuousEuclid{1})
addFactor!(fg, [:x0;:x1], LinearRelative(Normal(1, 0.1)))
addVariable!(fg, :x2, ContinuousEuclid{1})
addFactor!(fg, [:x1;:x2], LinearRelative(Normal(1, 0.1)))
addVariable!(fg, :x3, ContinuousEuclid{1})
addFactor!(fg, [:x2;:x3], LinearRelative(Normal(1, 0.1)))
addVariable!(fg, :x4, ContinuousEuclid{1})
addFactor!(fg, [:x3;:x4], LinearRelative(Normal(1, 0.1)))
altiMeasure = 10
addFactor!(fg, [:x3, :m], Scalar(Normal(altiMeasure, 0.1)))

# localizaton
addConstant!(fg,  :m, ScalarField(...))
addFactor!(fg, [:x1, :m], ScalarMeasurement(Normal(), ))

addConstant!(fg, :wmm, VectorField( "nasa_wmm_here"))
# wmm: V = f(x,y,z)

## DF mental blocker on `Variable` aspect here
addVariable!(fg, :wmm, SpatialVectorField( )) # for bar magnet
addFactor!(fg, [:x7, :wmm], VectorMeas(y_mag))
# mapping
addVariable!(fg, :m, BandLimitedScalarField(.....))
addVariable!(fg, :m, HistogramScalarField(.....))
addVariable!(fg, :m, GriddedScalarField(.....))
##

# variation: static vs dynamic
#  - static does not convey as much as addConstant

# addFactor -> addSampling, addSampler, ...

# addVariable!
# addFactor!
# addConstant?
#  - value varies of space and time

# addVariation good and should be in API for e.g. calibration, etc. reasons.
#  - ambiguous time vs. space varying aspect
#  - whos the target audience of the API naming
#  - too close to sensitivity/perturbation analysis lingo
#  - VicPrk eg, addVariation! (... wheelbase ... ), addVariation!( ... wheeldiameter ... )

#
# for 1b
# (x,h) samples on band-limited terrain
# addFactor!(fg, :x6, :m, ScalarBeam(h,b))

## to get to factor graph
#=
1. 1D localization
  a.  Unaries (known map) Keep map in fg node somehow (gaps/duty-cycle measurements)
  b.  mapping -- terrain as a variable to be estimated; constraints are now
      pairwise (map no longer encapsulated in unaries)
2. 1D via gaps -> 2D no map (line intersect)-> 2D point intersects only
  a. beam pattern somewhere in here too.

1
=#
##


text:
using Plots, JLD2, FileIO, Plots.PlotMeasures

plotly()

function CAnucleartest()
    # resultslist = Dict()
    # capacity = Dict()
    # energy = Dict()
    JLD2.@load "nuclearcosts2.jld2" resultslist capacity energy
    for nuccost in 8500:500:9000
        options, hourinfo, sets, params = buildsetsparams(carboncap=0.0, hydroinvestmentsallowed=false)
        println("\n\n\nNew run: nuclear cost=$nuccost.")
        params.investcost[:nuclear,:_] = nuccost
        model = buildvarsmodel(options, hourinfo, sets, params)
        println("\nSolving model...")
        status = solve(model.modelname)
        println("\nSolve status: $status")
        resultslist[nuccost] = sum(getvalue(model.vars.Systemcost)) / sum(params.demand) * 1000
        capacity[nuccost] = sum(getvalue(model.vars.Capacity[r,:nuclear,:_]) for r in sets.REGION)
        energy[nuccost] = sum(getvalue(model.vars.Electricity[r,:nuclear,:_,h]) for r in sets.REGION, h in sets.HOUR) / sum(params.demand)
        JLD2.@save "nuclearcosts2.jld2" resultslist capacity energy
        println("\nReading results...")
        results = readresults(model, status)
        name = autorunname(model.options) * ", nuccost=$nuccost"
        println("\nSaving results to disk...")
        saveresults(results, name, resultsfile="nuclearruns.jld2")
    end
    resultslist, capacity
end

function build_costs_from_runs(runsname, outname)
    path = "D:\\model runs\\"
    runs = listresults(resultsfile="$path$runsname")
    resultslist = Dict()
    allstatus = Dict()
    for (i, runname) in enumerate(runs)
        println("$i/$(length(runs)): $runname")
        r = loadresults(runname, resultsfile="$path$runsname")
        resultslist[runname] = sum(r.Systemcost)
        allstatus[runname] = r.Status
    end
    JLD2.@save "$path$outname" resultslist allstatus
end

function newrun!(namesuffix, options, hourinfo, sets, params, resultslist, allstatus, path, runsuffix)
    name = (namesuffix == "") ? autorunname(options) : autorunname(options) * ", " * namesuffix
    println("\n\n\nNew run: $name.")
    model = buildvarsmodel(options, hourinfo, sets, params)
    println("\nSolving model...")
    status = solve(model.modelname)
    println("\nSolve status: $status")
    resultslist[name] = sum(getvalue(model.vars.Systemcost))
    allstatus[name] = status
    JLD2.@save "$(path)supergridcosts$runsuffix.jld2" resultslist allstatus
    println("\nReading results...")
    results = readresults(model, status)
    println("\nSaving results to disk...")
    saveresults(results, name, resultsfile="$(path)supergridruns$runsuffix.jld2")
end

function supergridruns1()
    resultslist = Dict()
    allstatus = Dict()
    path = "D:\\model runs\\"
    runsuffix = "_dec6"
    for tm in [:none, :islands, :all]
        for cap in [1, 0.2, 0.1, 0.05, 0.02, 0.01, 0.005, 0.002, 0.001, 0]
            for solarwind in [1, 2], nuc in [false]
                options, hourinfo, sets, params = buildsetsparams(regionset=:Eurasia21,
                    islandindexes=[1:8, 9:15, 16:21], hours=1, maxbioenergy=0.05,
                    sspscenario="ssp2-26", sspyear=2050, datayear=2017,
                    carboncap=cap, solarwindarea=solarwind, nuclearallowed=nuc,
                    transmissionallowed=tm, disabletechs=[:csp])
                newrun!("", options, hourinfo, sets, params, resultslist, allstatus, path, runsuffix)
            end
            for solarwind in [1], nuc in [true]
                options, hourinfo, sets, params = buildsetsparams(regionset=:Eurasia21,
                    islandindexes=[1:8, 9:15, 16:21], hours=1, maxbioenergy=0.05,
                    sspscenario="ssp2-26", sspyear=2050, datayear=2017,
                    carboncap=cap, solarwindarea=solarwind, nuclearallowed=nuc,
                    transmissionallowed=tm, disabletechs=[:csp])
                newrun!("", options, hourinfo, sets, params, resultslist, allstatus, path, runsuffix)
            end
        end
    end
    resultslist, allstatus
end

function supergridruns2()
    path = "D:\\model runs\\"
    runsuffix = "_dec6"
    JLD2.@load "$(path)supergridcosts$runsuffix.jld2" resultslist allstatus
    for tm in [:none, :islands, :all]
        for cap in [1, 0.2, 0.1, 0.05, 0.02, 0.01, 0.005, 0.002, 0.001, 0]
            for solarwind in [1, 2], nuc in [false]
                options, hourinfo, sets, params = buildsetsparams(regionset=:Eurasia21,
                    islandindexes=[1:8, 9:15, 16:21], hours=1, maxbioenergy=0.05,
                    sspscenario="ssp2-26", sspyear=2050, datayear=2017,
                    carboncap=cap, solarwindarea=solarwind, nuclearallowed=nuc,
                    transmissionallowed=tm, disabletechs=[:csp])
                params.transmissioninvestcost .*= 0.5
                params.transmissionfixedcost .*= 0.5
                newrun!("half_transmission_cost", options, hourinfo, sets, params, resultslist, allstatus, path, runsuffix)

                options, hourinfo, sets, params = buildsetsparams(regionset=:Eurasia21,
                    islandindexes=[1:8, 9:15, 16:21], hours=1, maxbioenergy=0.05,
                    sspscenario="ssp2-26", sspyear=2050, datayear=2017,
                    carboncap=cap, solarwindarea=solarwind, nuclearallowed=nuc,
                    transmissionallowed=tm, disabletechs=[])
                newrun!("", options, hourinfo, sets, params, resultslist, allstatus, path, runsuffix)
            end
        end
    end
    resultslist, allstatus
end

function supergridruns3()
    path = "D:\\model runs\\"
    runsuffix = "_dec6"
    JLD2.@load "$(path)supergridcosts$runsuffix.jld2" resultslist allstatus
    for combo in [("ssp2-26", 2017), ("ssp2-26", 1998), ("ssp2-26", 2003), ("ssp1-45", 2017)]
    # for combo in [("ssp2-26", 2003), ("ssp1-45", 2017)]
        sspscen, erayear = combo
        for tm in [:islands, :all], cap in [0.001]
            for solarwind in [1, 2], nuc in [false]
                for solar in [:high, :mid, :low], battery in [:high, :mid, :low]
                    options, hourinfo, sets, params = buildsetsparams(regionset=:Eurasia21,
                        islandindexes=[1:8, 9:15, 16:21], hours=1, maxbioenergy=0.05,
                        sspscenario=sspscen, sspyear=2050, datayear=erayear,
                        carboncap=cap, solarwindarea=solarwind, nuclearallowed=nuc,
                        transmissionallowed=tm, disabletechs=[:csp])
                    solarbatterycosts!(sets, params, solar, battery)
                    newrun!("solar=$solar, battery=$battery", options, hourinfo, sets, params, resultslist, allstatus, path, runsuffix)
                end
            end
        end
    end
    resultslist, allstatus
end

function solarbatterycosts!(sets, params, solar, battery)
    pvcost = params.investcost[:pv,:a1]
    pvroofcost = params.investcost[:pvroof,:a1]
    batterycost = params.investcost[:battery,:_]
    for c in sets.CLASS[:pv]
        if solar == :high
            params.investcost[:pv,c] = pvcost + 300
            params.investcost[:pvroof,c] = pvroofcost + 300
        elseif solar == :mid
            params.investcost[:pv,c] = pvcost
            params.investcost[:pvroof,c] = pvroofcost
        elseif solar == :low
            params.investcost[:pv,c] = pvcost - 300
            params.investcost[:pvroof,c] = pvroofcost - 300
        end
    end
    if battery == :high
        params.investcost[:battery,:_] = batterycost * 1.5
    elseif battery == :mid
        params.investcost[:battery,:_] = batterycost
    elseif battery == :low
        params.investcost[:battery,:_] = batterycost * 0.5
    end
    nothing
end

function gispaperruns(runsuffix="_nov14", discountrate=0.07)
    path = "D:\\model runs\\"
    resultsfile = "$(path)gispaper_mixes$runsuffix.jld2"
    for carboncap in [0.025], landarea in [1, 2], region in [:China6, :Europe8]
        println("\n\n\nNew run: carboncap=$carboncap, landarea=$landarea, region=$region.")
        runmodel(regionset=region, carboncap=carboncap, discountrate=discountrate,
                        nuclearallowed=false, solarwindarea=landarea, resultsfile=resultsfile);
        if region == :Europe8
            println("\n\n\nNew run: carboncap=$carboncap, landarea=$landarea, region=$region, no transmission.")
            runmodel(regionset=region, carboncap=carboncap, discountrate=discountrate,
                    nuclearallowed=false, solarwindarea=landarea, transmissionallowed=:none, resultsfile=resultsfile);
        end
    end
end

# GIS paper figure 4
function plot_gispaper_mixes1(runsuffix="_nov14", carboncap=0.025, discountrate=0.07)
    path = "D:\\model runs\\"
    resultsfile = "$(path)gispaper_mixes$runsuffix.jld2"
    scen = ["default", "high land"]
    resultsnames = ["discountrate=$discountrate, nuclearallowed=false, carboncap=$carboncap, resultsfile=$resultsfile",
                    "discountrate=$discountrate, nuclearallowed=false, carboncap=$carboncap, resultsfile=$resultsfile, solarwindarea=2"]
    chart_energymix_scenarios(scen, resultsnames, resultsfile, size=(480,550), bar_width=0.7, xlims=(0.3,2.7), title="Europe")
    scen = ["default", "high land"]
    resultsnames = ["regionset=China6, discountrate=$discountrate, nuclearallowed=false, carboncap=$carboncap, resultsfile=$resultsfile",
                    "regionset=China6, discountrate=$discountrate, nuclearallowed=false, carboncap=$carboncap, resultsfile=$resultsfile, solarwindarea=2"]
    chart_energymix_scenarios(scen, resultsnames, resultsfile, size=(480,550), bar_width=0.7, xlims=(0.3,2.7), title="China")
end

# GIS paper figure 5
function plot_gispaper_mixes2(runsuffix="_nov14", carboncap=0.025, discountrate=0.07)
    path = "D:\\model runs\\"
    resultsfile = "$(path)gispaper_mixes$runsuffix.jld2"
    scen = ["default", "no transmission"]
    resultsnames = ["discountrate=$discountrate, nuclearallowed=false, carboncap=$carboncap, resultsfile=$resultsfile",
                    "transmissionallowed=none, discountrate=$discountrate, nuclearallowed=false, carboncap=$carboncap, resultsfile=$resultsfile"]
    chart_energymix_scenarios(scen, resultsnames, resultsfile, size=(480,550), bar_width=0.7, xlims=(0.3,2.7), title="Europe")
end

# GIS paper figure 6
function plot_gispaper_springmonth(runsuffix="_nov14", carboncap=0.025, discountrate=0.07)
    path = "D:\\model runs\\"
    resultsfile = "$(path)gispaper_mixes$runsuffix.jld2"
    r = loadresults("transmissionallowed=none, discountrate=$discountrate, nuclearallowed=false, carboncap=$carboncap, resultsfile=$resultsfile", resultsfile=resultsfile)
    annualelec, capac, tcapac, chart = analyzeresults(r)
    chart(:FRA, xlims=(1872,1872+722), ylims=(0,255))
    chart(:SPA, xlims=(1872,1872+722), ylims=(0,120))
end

# GIS paper figures 3 and 7
function plot_gispaper_classes1(runsuffix="_nov14", carboncap=0.025, discountrate=0.07)
    path = "D:\\model runs\\"
    resultsfile = "$(path)gispaper_mixes$runsuffix.jld2"
    r = loadresults("transmissionallowed=none, discountrate=$discountrate, nuclearallowed=false, carboncap=$carboncap, resultsfile=$resultsfile", resultsfile=resultsfile)
    annualelec, capac, tcapac, chart = analyzeresults(r)
    chart(:TOT)
    chart(:BARS, ylims=(0,1150))
end

# GIS paper figure 3 and 8
function plot_gispaper_classes2(runsuffix="_nov14", carboncap=0.025, discountrate=0.07)
    path = "D:\\model runs\\"
    resultsfile = "$(path)gispaper_mixes$runsuffix.jld2"
    r = loadresults("discountrate=$discountrate, nuclearallowed=false, carboncap=$carboncap, resultsfile=$resultsfile", resultsfile=resultsfile)
    annualelec, capac, tcapac, chart = analyzeresults(r)
    chart(:TOT)
    chart(:BARS, ylims=(0,1150))
end

function mergeresults()
    # JLD2.@load "iewcosts1_0.jld2" resultslist allstatus
    # res0, st0 = resultslist, allstatus
    # JLD2.@load "iewcosts1.jld2" resultslist allstatus
    # resultslist[false, 1, :none, 1.0] = res0[false, 1, :none, 1.0]
    # allstatus[false, 1, :none, 1.0] = st0[false, 1, :none, 1.0]
    # JLD2.@save "iewcosts1.jld2" resultslist allstatus
    JLD2.@load "iewcosts2 - part1.jld2" resultslist allstatus
    res1, st1 = resultslist, allstatus
    JLD2.@load "iewcosts2 - part2.jld2" resultslist allstatus
    res2, st2 = resultslist, allstatus
    resultslist = merge(res1, res2)
    allstatus = merge(st1, st2)
    JLD2.@save "iewcosts2.jld2" resultslist allstatus
end

function plot_supergridpaper_lines(ssp="ssp2-26", datayear=2017; runsuffix="_dec6")
    path = "D:\\model runs\\"
    resultsfile = "$(path)supergridruns$runsuffix.jld2"
    oldresults = (runsuffix == "_apr14")
    r = loadresults(runname(1, :all, 1.0, ssp, datayear, false, false, false; oldresults=oldresults), resultsfile=resultsfile)
    totaldemand = sum(r.params[:demand])    # 1.975066479481802e7   # (GWh/yr)
    existinghydro = sumdimdrop(r.Electricity[:hydro,:x0], dims=1)
    JLD2.@load "$(path)supergridcosts$runsuffix.jld2" resultslist allstatus
    res = resultslist
    carboncaps = Any[1000; 200; 100; 50; 20; 10; 5; 2; 1; 0]

    function getresults(a,b,c,d,e,f,g,h)
        cost = get(res, runname(a,b,c,d,e,f,g,h; oldresults=oldresults), NaN)    # M€/year
        return cost > 1e7 ? NaN : cost/(totaldemand - sum(existinghydro))*1000 # €/MWh
    end

    resmat2 = [getresults(1, tm, cap/1000, ssp, datayear, false, false, false)
                    for cap in carboncaps, tm in [:islands, :all]]
    resmat4 = [getresults(2, tm, cap/1000, ssp, datayear, false, false, false)
                    for cap in carboncaps, tm in [:islands, :all]]
    carboncaps[1] = "none"

    plotly()
    p = plot(string.(carboncaps), [resmat2 resmat4], color=[1 2 1 2], line=[:solid :solid :dash :dash])
    display(plot(p, size=(850,450), ylim=(0,61),
                    label=["C" "S" "C-Hland" "S-Hland"],
                    line=3, tickfont=14, legendfont=14,
                    titlefont=16, guidefont=14, xlabel="Global CO2 cap [g CO2/kWh]", ylabel="Average system cost [€/MWh]",
                    legend = :outertopright, left_margin=50px, gridlinewidth=1))
    p2 = plot(string.(carboncaps), resmat2, title="Low solar & wind area", label=["" ""])
    p4 = plot(string.(carboncaps), resmat4, title="High solar & wind area", label=["islands" "all"])
    display(plot(p2, p4, layout=2, size=(1000,450), ylim=(0,60), line=3, tickfont=14, legendfont=14,
                    titlefont=16, guidefont=14, xlabel="Global CO2 cap [g CO2/kWh]", ylabel="Average system cost [€/MWh]",
                    legend = :outertopright, left_margin=50px, gridlinewidth=1))

end

# Figure 2 in the supergrid paper (fig 1 is the eurasia map).
# SG.plot_supergridpaper_lines_transmission("ssp2-26", 2017, false, false, false, hidetitle=true)
# Supplementary figures:
# S1: SG.plot_supergridpaper_lines_transmission("ssp2-26", 2017, false, false, true, hidescenario=true)
# S2: SG.plot_supergridpaper_lines_transmission("ssp2-26", 2017, false, true, false, hidescenario=true)
# S3: SG.plot_supergridpaper_lines_transmission("ssp2-26", 2017, true, false, false, hidescenario=true)
function plot_supergridpaper_lines_transmission(ssp="ssp2-26", datayear=2017,
            allowcsp=false, allownuclear=false, halftmcost=false;
            hidetitle=false, hidescenario=false, runsuffix="_dec6")
    path = "D:\\model runs\\"
    resultsfile = "$(path)supergridruns$runsuffix.jld2"
    oldresults = (runsuffix == "_apr14")
    r = loadresults(runname(1, :all, 1.0, ssp, datayear, allowcsp, allownuclear, halftmcost; oldresults=oldresults), resultsfile=resultsfile)
    totaldemand = sum(r.params[:demand])    # 1.975066479481802e7   # (GWh/yr)
    existinghydro = sumdimdrop(r.Electricity[:hydro,:x0], dims=1)
    JLD2.@load "$(path)supergridcosts$runsuffix.jld2" resultslist allstatus
    res = resultslist
    # JLD2.@load "$(path)supergridcosts2$runsuffix.jld2" resultslist allstatus
    # res = merge(res, resultslist)
    carboncaps = Any[1000; 200; 100; 50; 20; 10; 5; 2; 1; 0]
    titletext = "CSP " * (allowcsp ? "" : "not ") * "allowed" *
        (halftmcost ? ", half transmission cost" : "") *
        (allownuclear ? ", nuclear allowed" : "") *
        (hidescenario ? "" : (oldresults ? " [old results]" : " [$ssp, $datayear]"))
    titletext = hidetitle ? "" : titletext

    function getresults(a,b,c,d,e,f,g,h)
        cost = get(res, runname(a,b,c,d,e,f,g,h; oldresults=oldresults), NaN)    # M€/year
        return cost > 1e7 ? NaN : cost/(totaldemand - sum(existinghydro))*1000 # €/MWh
    end

    resmat1 = [getresults(1, tm, cap/1000, ssp, datayear, allowcsp, allownuclear, halftmcost)
                    for cap in carboncaps, tm in [:none, :islands, :all]]
    resmat2 = [getresults(2, tm, cap/1000, ssp, datayear, allowcsp, allownuclear, halftmcost)
                    for cap in carboncaps, tm in [:none, :islands, :all]]
    carboncaps[1] = "none"

    plotly()
    p = plot(string.(carboncaps), [resmat1 resmat2], color=[3 1 2 3 1 2], line=[:solid :solid :solid :dash :dash :dash])
    display(plot(p, size=(850,500), ylim=(0,81), fg_legend = :transparent,
        label=["R           " "C" "S" "R-Hland" "C-Hland" "S-Hland"],
        line=3, tickfont=14, legendfont=14, yticks=0:10:80, legend = :outertopright,
        titlefont=16, guidefont=14, xlabel="Global CO2 cap [g CO2/kWh]", ylabel="Average system cost [€/MWh]",
        left_margin=50px, top_margin=10px, right_margin=0px, bottom_margin=8px,
        gridlinewidth=1, title=titletext))
    return [resmat1 100*round.(1 .- resmat2./resmat1, digits=3) resmat2]
end

# Other Supergrid figures:
# Figure S4:
# using GlobalEnergyGIS; GE = GlobalEnergyGIS
# Figures S5 & S6:
# using GlobalEnergyGIS; GE = GlobalEnergyGIS
# GISsolar(gisregion="France13", plotmasks=:onlymasks)
# GISwind(gisregion="France13", plotmasks=:onlymasks)
# GISsolar(gisregion="KZK", plotmasks=:onlymasks)
# GISwind(gisregion="KZK", plotmasks=:onlymasks)
# GISsolar(gisregion="CH_E", plotmasks=:onlymasks)
# GISwind(gisregion="CH_E", plotmasks=:onlymasks)
# Figure S7:
# using GlobalEnergyGIS; GE = GlobalEnergyGIS
# using Plots, Plots.PlotMeasures
# plotly()
# absindex, index, meanwind, annualwind = annualwindindex(gisregion="Eurasia21", sites_quantile=0.25, aggregateregions=[1:8, 9:15, 16:21]);
# plot(1979:2019, absindex, c=[1 3 2], label=["Europe"  "CAS"  "China"], ylabel="[m/s]", tickfont=14, lw=2, legend=:outertopright, size=(800,550), legendfontsize=12, marker=:circle, markersize=3, markerstrokewidth=0, fg_legend = :transparent, title="Average wind speeds from ERA5 data (top 25% of pixels)")

runname(land::Int, tm::Symbol, cap::Float64, ssp::String, datayear::Int, allowcsp::Bool,
        allownuclear::Bool, halftmcost::Bool; oldresults=false) =
    "regionset=Eurasia21, " * (tm == :all ? "" : "transmissionallowed=$tm, ") *
    (datayear == 2017 ? "" : "datayear=$datayear, ") *
    (allownuclear ? "" : "nuclearallowed=false, ") *
    (cap == 1.0 ? "" : "carboncap=$cap, ") *
    (land == 1 ? "" : "solarwindarea=2, ") *
    (allowcsp ? "" : "disabletechs=$(oldresults ? "Symbol" : "")[:csp], ") *
    (ssp == "ssp2-34" ? "" : "sspscenario=$ssp, ") *
    "islandindexes=UnitRange{Int64}[1:8, 9:15, 16:21]" *
    (halftmcost ? ", half_transmission_cost" : "")

runname(land::Int, tm::Symbol, ssp::String, datayear::Int, solar::Symbol, battery::Symbol;
        oldresults=false) =
    "regionset=Eurasia21, " * (tm == :all ? "" : "transmissionallowed=$tm, ") *
    (datayear == 2017 ? "" : "datayear=$datayear, ") *
    "nuclearallowed=false, carboncap=0.001, " *
    (land == 1 ? "" : "solarwindarea=2, ") *
    "disabletechs=$(oldresults ? "Symbol" : "")[:csp], " *
    (ssp == "ssp2-34" ? "" : "sspscenario=$ssp, ") *
    "islandindexes=UnitRange{Int64}[1:8, 9:15, 16:21], solar=$solar, battery=$battery"

# Figure 4 in the supergrid paper.
# SG.plot_supergridpaper_bubbles("ssp2-26", 2017, hidescenario=true)
# Supplementary figure S4, S5, S6:
# SG.plot_supergridpaper_bubbles("ssp1-45", 2017, hidescenario=false)
# SG.plot_supergridpaper_bubbles("ssp2-26", 1998, hidescenario=false)
# SG.plot_supergridpaper_bubbles("ssp2-26", 2003, hidescenario=false)
function plot_supergridpaper_bubbles(ssp="ssp2-34", datayear=2017;
                hidescenario=false, runsuffix="_dec6")
    path = "D:\\model runs\\"
    JLD2.@load "$(path)supergridcosts$runsuffix.jld2" resultslist allstatus
    oldresults = (runsuffix == "_apr14")
    res = resultslist
    # showall(keys(res))
    rows = [3 3 3 2 2 2 1 1 1]
    cols = [3 2 1 3 2 1 3 2 1]
    r1 = [res[runname(1,:islands,ssp,datayear,solar,battery; oldresults=oldresults)] /
            res[runname(1,:all,ssp,datayear,solar,battery; oldresults=oldresults)]-1
                for solar in [:high, :mid, :low], battery in [:high, :mid, :low]]
    r2 = [res[runname(2,:islands,ssp,datayear,solar,battery; oldresults=oldresults)] /
            res[runname(2,:all,ssp,datayear,solar,battery; oldresults=oldresults)]-1
                for solar in [:high, :mid, :low], battery in [:high, :mid, :low]]
    bs = 0.4    # bubble size
    annotations1 = [(rows[i]-0.65*bs*sqrt(r1[i]*100/pi), cols[i], text("$(round(r1[i]*100, digits=1))%", :right)) for i=1:9]
    annotations2 = [(rows[i]-0.65*bs*sqrt(r2[i]*100/pi), cols[i], text("$(round(r2[i]*100, digits=1))%", :right)) for i=1:9]
    titlesuffix = hidescenario ? "" : (oldresults ? " [old results]" : " [$ssp, $datayear]")
    s1 = scatter(rows, cols, markersize=reshape(sqrt.(r1*100/pi)*75*bs, (1,9)),
        annotations=annotations1, xlim=(0.4,3.4), ylim=(0.5,3.5), legend=false,
        title="Default solar & wind area" * titlesuffix,
        xlabel="battery cost", ylabel="solar PV cost", color=1, tickfont=14, guidefont=14)
    xticks!([1,2,3],["low","mid","high"])
    yticks!([1,2,3],["low","mid","high"])
    s2 = scatter(rows, cols, markersize=reshape(sqrt.(r2*100/pi)*75*bs, (1,9)),
        annotations=annotations2, xlim=(0.4,3.4), ylim=(0.5,3.5), legend=false,
        title="High solar & wind area" * titlesuffix,
        xlabel="battery cost", ylabel="solar PV cost", color=1, tickfont=14, guidefont=14,
        left_margin=20px)
    xticks!([1,2,3],["low","mid","high"])
    yticks!([1,2,3],["low","mid","high"])
    # display(plot(s2, size=(500,450)))
    display(plot(s1, s2, layout=2, size=(1000,450)))
end

function plot_supergridpaper_bubbles_abscost(transmission=:all, ssp="ssp2-34", datayear=2017;
            hidescenario=false, runsuffix="_dec6")
    path = "D:\\model runs\\"
    JLD2.@load "$(path)supergridcosts$runsuffix.jld2" resultslist allstatus
    oldresults = (runsuffix == "_apr14")
    res = resultslist
    # showall(keys(res))
    rows = [3 3 3 2 2 2 1 1 1]
    cols = [3 2 1 3 2 1 3 2 1]
    r1 = [res[runname(1,transmission,ssp,datayear,solar,battery; oldresults=oldresults)]/1e6
        for solar in [:high, :mid, :low], battery in [:high, :mid, :low]]
    r2 = [res[runname(2,transmission,ssp,datayear,solar,battery; oldresults=oldresults)]/1e6
        for solar in [:high, :mid, :low], battery in [:high, :mid, :low]]
    bs = 0.4    # bubble size
    annotations1 = [(rows[i]-0.65*bs*sqrt(r1[i]/pi), cols[i], text("$(round(r1[i], digits=3))", :right)) for i=1:9]
    annotations2 = [(rows[i]-0.65*bs*sqrt(r2[i]/pi), cols[i], text("$(round(r2[i], digits=3))", :right)) for i=1:9]
    titlesuffix = hidescenario ? "" : (oldresults ? " [old results]" : " [$transmission, $ssp, $datayear]")
    s1 = scatter(rows, cols, markersize=reshape(sqrt.(r1/pi)*75*bs, (1,9)),
    annotations=annotations1, xlim=(0.4,3.4), ylim=(0.5,3.5), legend=false,
    title="Default area" * titlesuffix,
    xlabel="battery cost", ylabel="solar PV cost", color=1, tickfont=14, guidefont=14)
    xticks!([1,2,3],["low","mid","high"])
    yticks!([1,2,3],["low","mid","high"])
    s2 = scatter(rows, cols, markersize=reshape(sqrt.(r2/pi)*75*bs, (1,9)),
    annotations=annotations2, xlim=(0.4,3.4), ylim=(0.5,3.5), legend=false,
    title="High area" * titlesuffix,
    xlabel="battery cost", ylabel="solar PV cost", color=1, tickfont=14, guidefont=14,
    left_margin=20px)
    xticks!([1,2,3],["low","mid","high"])
    yticks!([1,2,3],["low","mid","high"])
    # display(plot(s2, size=(500,450)))
    display(plot(s1, s2, layout=2, size=(1000,450)))
end

function plot_supergridpaper_all_energymixes()
    # PAPER VERSIONS (no R scenario, unused legend techs removed)
    # plot_supergridpaper_energymix("ssp2-34", 2018, false, false, false,
    #     indices=[2,3,5,6], deletetechs=[1,2,6,7,12], runsuffix="_apr14")
    # plot_supergridpaper_energymix("ssp2-34", 2018, true, false, false,
    #     indices=[2,3,4,5,6], deletetechs=[1,2,12], runsuffix="_apr14")
    # plot_supergridpaper_energymix("ssp2-34", 2018, false, true, false,
    #     indices=[1,2,3,4,5,6], deletetechs=[2,6,7,12], runsuffix="_apr14")
    # plot_supergridpaper_energymix("ssp2-34", 2018, false, false, true,
    #     indices=[2,3,4,5,6], deletetechs=[1,2,6,7,12], runsuffix="_apr14")

    # DISCUSSION VERSIONS (only infeasible runs removed)
    plot_supergridpaper_energymix("ssp2-26", 2017, false, false, false,
        indices=[2,3,4,5,6], runsuffix="_dec6")
    plot_supergridpaper_energymix("ssp2-26", 2017, true, false, false,
        indices=[2,3,4,5,6], runsuffix="_dec6")
    plot_supergridpaper_energymix("ssp2-26", 2017, false, true, false,
        indices=[1,2,3,4,5,6], runsuffix="_dec6")
    plot_supergridpaper_energymix("ssp2-26", 2017, false, false, true,
        indices=[2,3,4,5,6], runsuffix="_dec6")
    # plot_supergridpaper_energymix("ssp2-26", 2017, false, false, false,
    #     indices=[2,3,5,6], runsuffix="_dec6",
    #     extrasuffix=", solar=mid, battery=mid", ylims=(0,34))
    # plot_supergridpaper_energymix("ssp2-26", 1998, false, false, false,
    #     indices=[2,3,5,6], runsuffix="_dec6",
    #     extrasuffix=", solar=mid, battery=mid", ylims=(0,34))
    # plot_supergridpaper_energymix("ssp2-26", 2003, false, false, false,
    #     indices=[2,3,5,6], runsuffix="_dec6",
    #     extrasuffix=", solar=mid, battery=mid", ylims=(0,34))
    # plot_supergridpaper_energymix("ssp1-45", 2017, false, false, false,
    #     indices=[2,3,5,6], runsuffix="_dec6",
    #     extrasuffix=", solar=mid, battery=mid", ylims=(0,34))
end

# Figure 3 in the supergrid paper.
# SG.plot_supergridpaper_energymix("ssp2-26", 2017, false, false, false, indices=[2,3,5,6], deletetechs=[1,2,6,7,12], hidetitle=true, runsuffix="_dec6")
function plot_supergridpaper_energymix(ssp="ssp2-26", datayear=2017, allowcsp=false,
            allownuclear=false, halftmcost=false;
            cap=0.001, indices=collect(1:6), deletetechs=[2,12],
            runsuffix="_dec6", extrasuffix="", hidetitle=false, figoptions...)
    path = "D:\\model runs\\"
    resultsfile = "$(path)supergridruns$runsuffix.jld2"
    oldresults = (runsuffix == "_apr14")
    # indices = (cap >= 0.1 || allownuclear) ? collect(1:6) : [2,3,5,6]
    if allownuclear
        indices = indices[indices .<= 3]
    end
    scen = ["R", "C", "S", "R-Hland", "C-Hland", "S-Hland"][indices]
    resultsnames = [
        runname(1, :none, cap, ssp, datayear, allowcsp, allownuclear, halftmcost; oldresults=oldresults),
        runname(1, :islands, cap, ssp, datayear, allowcsp, allownuclear, halftmcost; oldresults=oldresults),
        runname(1, :all, cap, ssp, datayear, allowcsp, allownuclear, halftmcost; oldresults=oldresults),
        runname(2, :none, cap, ssp, datayear, allowcsp, allownuclear, halftmcost; oldresults=oldresults),
        runname(2, :islands, cap, ssp, datayear, allowcsp, allownuclear, halftmcost; oldresults=oldresults),
        runname(2, :all, cap, ssp, datayear, allowcsp, allownuclear, halftmcost; oldresults=oldresults)
    ][indices] .* extrasuffix
    cap1000 = round(Int, 1000*cap)
    titletext = "CSP " * (allowcsp ? "" : "not ") * "allowed" *
        (halftmcost ? ", half transmission cost" : "") *
        (allownuclear ? ", nuclear allowed" : "") *
        ", $cap1000 g CO2/kWh" *
        (oldresults ? " [old results]" : " [$ssp, $datayear]")
    titletext = hidetitle ? "" : titletext
    plot_energymix(scen, resultsnames, resultsfile; size=(200+100*length(indices), 550),
        title=titletext, deletetechs=deletetechs, figoptions...)
    # plot_energymix(["default land", "high land"],
    #     ["regionset=Eurasia21, transmissionallowed=islands, nuclearallowed=false, carboncap=0.001, disabletechs=Symbol[:csp], islandindexes=UnitRange{Int64}[1:8, 9:15, 16:21], solar=low, battery=high",
    #      "regionset=Eurasia21, transmissionallowed=islands, nuclearallowed=false, carboncap=0.001, solarwindarea=2, disabletechs=Symbol[:csp], islandindexes=UnitRange{Int64}[1:8, 9:15, 16:21], solar=low, battery=high"],
    #     resultsfile; size=(400, 550), title="Solar low cost, battery high cost", deletetechs=[1,2,6,7,12])
end

function plot_energymix(scen, resultsnames, resultsfile; deletetechs=[], optionlist...)
    scenelec, demands, hoursperperiod, displayorder, techlabels = iew_getscenarioresults(scen, resultsnames, resultsfile)

    palette = [RGB([216,137,255]/255...), RGB([119,112,71]/255...), RGB([199,218,241]/255...), RGB([149,179,215]/255...),
        RGB([255,255,64]/255...), RGB([240,224,0]/255...), RGB([214,64,64]/255...), RGB([99,172,70]/255...), RGB([255,192,0]/255...),
        RGB([100,136,209]/255...), RGB([144,213,93]/255...), RGB([148,138,84]/255...), RGB([157,87,205]/255...)]

    deleteat!(palette, deletetechs)
    techlabels = permutedims(deleteat!(vec(techlabels), deletetechs))
    deleteat!(displayorder, deletetechs)

    display(demands'*hoursperperiod/1e6)
    display(techlabels)
    display(scen)
    display(collect(scenelec[displayorder,:]')/1e6)
    stackedbar(collect(scenelec[displayorder,:]')/1e6; label=string.(techlabels),
        size=(600,550), left_margin=20px, fg_legend = :transparent,
        xticks=(1:length(scen),scen), line=0, tickfont=12, legendfont=12, guidefont=12,
        color_palette=palette, ylabel="[PWh/year]", legend=:outertopright, optionlist...)
    xpos = (1:length(scen))'
    lab = fill("",(1,length(scen)))
    lab[1] = "demand"
    display(plot!([xpos; xpos], [zeros(length(scen))'; demands'*hoursperperiod/1e6], line=3, color=:black, label=lab))
    nothing
end

function regionalcosts(ssp="ssp2-34", allowcsp=false, allownuclear=false, halftmcost=false; cap=0.001, indices=collect(1:6), runsuffix="_apr14")
    path = "D:\\model runs\\"
    resultsfile = "$(path)supergridruns$runsuffix.jld2"
    # indices = (cap >= 0.1 || allownuclear) ? collect(1:6) : [2,3,5,6]
    if allownuclear
        indices = indices[indices .<= 3]
    end
    islandindexes=[1:8, 9:15, 16:21, 1:21]
    scen = ["R", "C", "S", "R-Hland", "C-Hland", "S-Hland"][indices]
    resultsnames = [runname(1, :none, cap, ssp, allowcsp, allownuclear, halftmcost),
                    runname(1, :islands, cap, ssp, allowcsp, allownuclear, halftmcost),
                    runname(1, :all, cap, ssp, allowcsp, allownuclear, halftmcost),
                    runname(2, :none, cap, ssp, allowcsp, allownuclear, halftmcost),
                    runname(2, :islands, cap, ssp, allowcsp, allownuclear, halftmcost),
                    runname(2, :all, cap, ssp, allowcsp, allownuclear, halftmcost)][indices]
    for resultname in resultsnames
        println(resultname, " ", resultsfile)
        r = loadresults(resultname, resultsfile=resultsfile)
        existinghydro = sumdimdrop(r.Electricity[:hydro,:x0], dims=1)
        elec = sum(r.Electricity[k,c] for k in r.sets.TECH for c in r.sets.CLASS[k])
        annualelec = sumdimdrop(elec, dims=1)
        demand = sumdimdrop(r.params[:demand], dims=2)
        cost_prod = r.Systemcost ./ (annualelec - existinghydro) * 1000
        cost_dem = r.Systemcost ./ (demand - existinghydro) * 1000
        regcost_prod = [sum(r.Systemcost[is]) ./ sum(annualelec[is] - existinghydro[is]) * 1000 for is in islandindexes]
        regcost_dem = [sum(r.Systemcost[is]) ./ sum(demand[is] - existinghydro[is]) * 1000 for is in islandindexes]
        # @show regcost_prod
        @show regcost_dem
        println()
    end
end

function iew_readscenariodata(resultname, resultsfile)
    println(resultname, " ", resultsfile)
    results = loadresults(resultname, resultsfile=resultsfile)
    @unpack TECH, REGION, CLASS, HOUR = results.sets
    hoursperperiod = results.hourinfo.hoursperperiod
    totaldemand = sum(results.params[:demand])
    totalelec = [sum(sum(results.Electricity[k,c]) for c in CLASS[k]) for k in TECH]

    techdisplayorder = [:nuclear, :coal, :wind, :offwind, :pv, :pvroof, :csp, :bioCCGT, :gasCCGT, :hydro, :bioGT, :gasGT, :battery]
    techlabels = [k for r=1:1, k in techdisplayorder]
    displayorder = [i for (i,k) in enumerate(TECH), d in techdisplayorder if d == k]

    return totalelec, totaldemand, hoursperperiod, displayorder, techlabels
end

function iew_getscenarioresults(scenarios, resultsnames, resultsfile)
    numscen = length(scenarios)
    scenelec = zeros(13,numscen)
    demands = zeros(numscen)
    hoursperperiod, displayorder, techlabels = nothing, nothing, nothing

    for (i,s) in enumerate(scenarios)
        println("Loading results: $s...")
        totalelec, totaldemand, hoursperperiod, displayorder, techlabels = iew_readscenariodata(resultsnames[i], resultsfile)
        scenelec[:,i] = totalelec
        demands[i] = totaldemand
    end
    return scenelec, demands, hoursperperiod, displayorder, techlabels
end

max_stars_repo_name: niclasmattsson/Supergrid

max_issues_repo_name: niclasmattsson/Supergrid

max_forks_repo_name: niclasmattsson/Supergrid


text:
@testset "Sample properties for arrays" begin
    @test Banyan.sample_axes(Base.ones(10,10)) == [1, 2]
    @test Banyan.sample_axes(Base.ones(10)) == [1]
    @test Banyan.sample_keys(Base.ones(10,10)) == [1, 2]
    @test Banyan.sample_keys(Base.ones(10)) == [1]

    increasing_array_1d = range(1,100,length=100)
    increasing_array_2d = reshape(range(1,100,length=100), (10,10))
    string_array_1d = Base.fill("abc", 10)
    string_array_2d = Base.fill("abc", (10,10))

    @test Banyan.sample_min(increasing_array_1d, 1) == 1
    @test Banyan.sample_max(increasing_array_1d, 1) == 100
    @test Banyan.sample_min(increasing_array_2d, 1) == 1
    @test Banyan.sample_max(increasing_array_2d, 1) == 10
    @test Banyan.sample_min(increasing_array_2d, 2) == 1
    @test Banyan.sample_max(increasing_array_2d, 2) == 91
    @test Banyan.sample_min(string_array_1d, 1) == orderinghash("abc")
    @test Banyan.sample_max(string_array_1d, 1) == orderinghash("abc")
    @test Banyan.sample_min(string_array_2d, 1) == orderinghash("abc")
    @test Banyan.sample_max(string_array_2d, 1) == orderinghash("abc")
    @test Banyan.sample_min(string_array_2d, 2) == orderinghash("abc")
    @test Banyan.sample_max(string_array_2d, 2) == orderinghash("abc")
end

@testset "Sample divisions for arrays" begin
    # Test outliers
    @test sample_divisions([5,6,7,8], 1) == [
        (5, 6)
        (6, 7)
        (7, 8)
        (8, 8)
    ]
    @test sample_divisions([5,6,7,8,8,8,8], 1) == [(5,8)]
    @test sample_divisions([5,6,7,8,8,8,20], 1) == [
        (5, 8)
        (8, 20)
    ]
    @test sample_divisions([5,6,7,8,8,8,8,8,20], 1) == [(5, 20)]

    # Test 2D arrays
    @test sample_divisions(reshape(repeat([5,6,7,8,8,8,20],4), (7,4)), 1) == [
        (5, 8)
        (8, 20)
    ]
    @test sample_divisions(reshape(repeat([5,6,7,8,8,8,20],4), (7,4)), 2) == [
        (5, 5)
    ]
    @test sample_divisions(reshape(repeat([5,6,7,8,8,8,8,8,20],4), (9,4)), 1) == [
        (5, 20)
    ]
end

@testset "Sample # of groups for arrays" begin
    @test Banyan.sample_max_ngroups(Base.ones(10), 1) == 1
    @test Banyan.sample_max_ngroups(Base.ones(10,10), 1) == 1
    @test Banyan.sample_max_ngroups(Base.ones(10,10), 2) == 1

    increasing_array_1d = range(1,100,length=100)
    increasing_array_2d = reshape(range(1,100,length=100), (10,10))
    @test Banyan.sample_max_ngroups(increasing_array_1d, 1) == 100
    @test Banyan.sample_max_ngroups(increasing_array_2d, 1) == 10
    @test Banyan.sample_max_ngroups(increasing_array_2d, 2) == 10
end

@testset "Sample percentile for arrays" begin
    increasing_array_1d = range(1,100,length=100)
    increasing_array_2d = reshape(range(1,100,length=100), (10,10))

    @test sample_percentile(increasing_array_1d, 1, 1, 20) == 0.2
    @test sample_percentile(increasing_array_1d, 1, 80, 100) == 0.21
    @test sample_percentile(increasing_array_1d, 1, 1, 100) == 1.0
    @test sample_percentile(increasing_array_1d, 1, 1, 1) == 0.01
    @test sample_percentile(increasing_array_1d, 1, -100, -10) == 0

    @test sample_percentile(increasing_array_2d, 1, 1, 20) == 1.0
    @test sample_percentile(increasing_array_2d, 1, 5, 10) == 0.6
    @test sample_percentile(increasing_array_2d, 1, 1, 1) == 0.1
    @test sample_percentile(increasing_array_2d, 1, -100, -10) == 0

    @test sample_percentile(["a", "b", "c", "d"], 1, orderinghash("a"), orderinghash("b")) == 0.5
    @test sample_percentile(["a", "b", "c", "d"], 1, orderinghash("b"), orderinghash("d")) == 0.75
    @test sample_percentile(["a", "b", "c", "d"], 1, orderinghash("e"), orderinghash("f")) == 0.0
    @test sample_percentile(["a", "b", "c", "d"], 1, orderinghash("d"), orderinghash("d")) == 0.25
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
@testset "Test the LS version of the OMP 0 relative error" begin

    ψ = randn(100, 20)
    u = randn(100)
    m, n = size(ψ)

    ctrue = ψ\u

    # Force a greedy approach over the entire set of indices
    # With LS
    idxls, cls, ϵls = lsomp(ψ, u; invert = true, verbose = true, ϵrel = eps())

    @test norm(ctrue[idxls]-cls)<1e-14

    # Residual error must decrease
    @test all(ϵls[2:end]-ϵls[1:end-1] .< 0) == true

    for k=1:n
        cverif = view(ψ,:, idxls[1:k])\u
        @test abs(ϵls[k+1] - norm(u - view(ψ,:, idxls[1:k])*cverif))<1e-14
    end
end

@testset "Test the LS version of the OMP given relative error" begin
    ψ = randn(100, 20)
    u = randn(100)
    m, n = size(ψ)

    ctrue = ψ\u

    ϵrel = 0.1*1.00 + 0.9*norm(u - ψ*ctrue)/norm(u)

    # Force a greedy approach over the entire set of indices
    # With LS
    idxls, cls, ϵls = lsomp(ψ, u; invert = true, verbose = true, ϵrel = ϵrel)

    # Residual error must decrease
    @test all(ϵls[2:end]-ϵls[1:end-1] .< 0) == true

    @test ϵls[end]<ϵrel*norm(u)
    @test ϵls[end-1]>ϵrel*norm(u)

    for k=1:size(idxls,1)
        cverif = view(ψ,:, idxls[1:k])\u
        @test abs(ϵls[k+1] - norm(u - view(ψ,:, idxls[1:k])*cverif))<1e-14
    end
end

@testset "Test the LS version of the OMP given a relative number of features" begin
    ψ = randn(100, 20)
    u = randn(100)
    m, n = size(ψ)
    maxterms = 10
    ctrue = ψ\u

    # Force a greedy approach over the entire set of indices
    # With LS
    idxls, cls, ϵls = lsomp(ψ, u; invert = true, verbose = true, maxterms = maxterms)

    # Residual error must decrease
    @test all(ϵls[2:end]-ϵls[1:end-1] .< 0) == true

    @test size(idxls, 1) == maxterms
    @test size(cls, 1) == maxterms
    @test size(ϵls, 1) == maxterms + 1

    for k=1:size(idxls,1)
        cverif = view(ψ,:, idxls[1:k])\u
        @test abs(ϵls[k+1] - norm(u - view(ψ,:, idxls[1:k])*cverif))<1e-14
    end
end

@testset "Test the QR version of the OMP" begin

    ψ = randn(100, 20)
    u = randn(100)
    m, n = size(ψ)
    ctrue = ψ\u

    # Force a greedy approach over the entire set of indices
    # With QR
    idxqr, cqr, ϵqr = qromp(ψ, u; invert = true, verbose = true, ϵrel = eps())

    @test norm(ctrue[idxqr]-cqr)<1e-14

    # Residual error must decrease
    @test all(ϵqr[2:end]-ϵqr[1:end-1] .< 0) == true

    for k=1:n
        cverif = view(ψ,:, idxqr[1:k])\u
        @test abs(ϵqr[k+1] - norm(u - view(ψ,:, idxqr[1:k])*cverif))<1e-14
    end
end

@testset "Test the QR version of the OMP given relative error" begin
    ψ = randn(100, 20)
    u = randn(100)
    m, n = size(ψ)

    ctrue = ψ\u

    ϵrel = 0.1*1.00 + 0.9*norm(u - ψ*ctrue)/norm(u)

    # Force a greedy approach over the entire set of indices
    # With QR
    idxqr, cqr, ϵqr = qromp(ψ, u; invert = true, verbose = true, ϵrel = ϵrel)

    # Residual error must decrease
    @test all(ϵqr[2:end]-ϵqr[1:end-1] .< 0) == true

    @test ϵqr[end]<ϵrel*norm(u)
    @test ϵqr[end-1]>ϵrel*norm(u)

    for k=1:size(idxqr,1)
        cverif = view(ψ,:, idxqr[1:k])\u
        @test abs(ϵqr[k+1] - norm(u - view(ψ,:, idxqr[1:k])*cverif))<1e-14
    end
end

@testset "Test the QR version of the OMP given a number of features" begin
    ψ = randn(100, 20)
    u = randn(100)
    m, n = size(ψ)

    maxterms = 10

    ctrue = ψ\u

    ϵrel = 0.1*1.00 + 0.9*norm(u - ψ*ctrue)/norm(u)

    # Force a greedy approach over the entire set of indices
    # With QR
    idxqr, cqr, ϵqr = qromp(ψ, u; invert = true, verbose = true, maxterms = maxterms)

    # Residual error must decrease
    @test all(ϵqr[2:end]-ϵqr[1:end-1] .< 0) == true

    @test size(idxqr, 1) == maxterms
    @test size(cqr, 1) == maxterms
    @test size(ϵqr, 1) == maxterms + 1

    for k=1:size(idxqr,1)
        cverif = view(ψ,:, idxqr[1:k])\u
        @test abs(ϵqr[k+1] - norm(u - view(ψ,:, idxqr[1:k])*cverif))<1e-14
    end
end

@testset "Check that LS and QR versions are consistent" begin

    ψ = randn(100, 20)
    u = randn(100)
    ctrue = ψ\u

    # Force a greedy approach over the entire set of indices
    # With LS
    idxls, cls, ϵls = lsomp(ψ, u; invert = true, verbose = true, ϵrel = eps())

    # With QR
    idxqr, cqr, ϵqr = qromp(ψ, u; invert = true, verbose = true, ϵrel = eps())

    @test norm(ctrue[idxls]-cls)<1e-14
    @test norm(ctrue[idxqr]-cqr)<1e-14
    # @test all(ϵls - ϵqr .>= 0) == true
end

max_forks_repo_name: mleprovost/QROMP