text:
# Add packages

using Vann
using DataAssim
using ExcelReaders
using DataFrames
using PyPlot

# Settings

file_exper  = "C:/Work/Studies/vann/experiments.xlsx"
path_inputs = "C:/Work/Studies/vann/data/norway"
path_save   = "C:/Work/Studies/vann/filter_res"

df_exper = readxlsheet(DataFrame, file_exper, "Filter")

# Run one experiment

function run_single_experiment(path_inputs, path_save, df_exper, i_exper)

  # Folder for saving results

  path_save = joinpath(path_save, "experiment_" * string(i_exper))

  mkpath(path_save * "/tables")
  mkpath(path_save * "/figures")

  # Get settings

  opt = Dict()

  opt["epot_choice"]   = parse(df_exper[:epot_choice][i_exper])
  opt["snow_choice"]   = parse(df_exper[:snow_choice][i_exper])
  opt["hydro_choice"]  = parse(df_exper[:hydro_choice][i_exper])
  opt["filter_choice"] = parse(df_exper[:filter_choice][i_exper])
  opt["nens"]          = round(Int64, df_exper[:nens][i_exper])
  opt["warmup"]        = round(Int64, df_exper[:warmup][i_exper])
  opt["date_start"]    = df_exper[:date_start][i_exper]
  opt["date_stop"]     = df_exper[:date_stop][i_exper]

  opt["path_inputs"]   = path_inputs
  opt["path_save"]     = path_save

  # Run for all stations

  run_all_stations(opt)

end

# Run over all stations

function run_all_stations(opt)

  # Empty dataframes for summary statistics

  df_summary = DataFrame(Station = String[], NSE = Float64[], KGE = Float64[])

  # Loop over all watersheds

  dir_all = readdir(opt["path_inputs"])

  for dir_cur in dir_all

    try

      # Load data

      date, tair, prec, q_obs, frac = load_data("$(opt["path_inputs"])/$dir_cur")

      # Crop data

      date, tair, prec, q_obs = crop_data(date, tair, prec, q_obs, opt["date_start"], opt["date_stop"])

      # Compute potential evapotranspiration

      epot = eval(Expr(:call, opt["epot_choice"], date))

      # Initilize model

      tstep = 24.0

      st_snow = eval(Expr(:call, opt["snow_choice"], tstep, date[1], frac))
      st_hydro = eval(Expr(:call, opt["hydro_choice"], tstep, date[1]))

      # Run calibration

      param_snow, param_hydro = run_model_calib(st_snow, st_hydro, tair, prec, epot, q_obs;
                                                warmup = opt["warmup"])

      # Run model and filter

      st_snow = eval(Expr(:call, opt["snow_choice"], tstep, date[1], param_snow, frac))
      st_hydro = eval(Expr(:call, opt["hydro_choice"], tstep,  date[1], param_hydro))

      q_res = eval(Expr(:call, opt["filter_choice"], st_snow, st_hydro, prec, tair, epot, q_obs, opt["nens"]))

      # Add results to dataframe

      q_obs = round(q_obs, 2)
      q_sim = round(q_res[:, 1], 2)
      q_min = round(q_res[:, 2], 2)
      q_max = round(q_res[:, 3], 2)

      df_res = DataFrame(date = date, q_obs = q_obs, q_sim = q_sim, q_min = q_min, q_max = q_max)

      # Compute summary statistics

      station = dir_cur[1:end-5]
      nse_res = nse(q_sim[opt["warmup"]:end], q_obs[opt["warmup"]:end])
      kge_res = kge(q_sim[opt["warmup"]:end], q_obs[opt["warmup"]:end])

      push!(df_summary, [station nse_res kge_res])

      # Save results to txt file

      file_name = joinpath(opt["path_save"], "tables", dir_cur[1:end-5] * "_station.txt")

      writetable(file_name, df_res, quotemark = '"', separator = '\t')

      # Plot results

      ioff()

      fig = plt[:figure](figsize = (12,7))

      plt[:style][:use]("ggplot")

      plt[:plot](date, q_obs, linewidth = 1.2, color = "k", label = "Observed", zorder = 1)
      plt[:fill_between](date, q_max, q_min, facecolor = "r", edgecolor = "r", label = "Simulated", alpha = 0.55, zorder = 2)
      plt[:ylabel]("Runoff (mm/day)")

      plt[:legend]()

      file_name = joinpath(opt["path_save"], "figures", dir_cur[1:end-5] * "_station.png")

      savefig(file_name, dpi = 600)
      close(fig)

    catch

      info("Unable to run files in directory $(dir_cur)\n")

    end

  end

  writetable(string(opt["path_save"], "/summary_table.txt"), df_summary, quotemark = '"', separator = '\t')

end

# Run all experiments

for i_exper = 1:size(df_exper, 1)

  run_single_experiment(path_inputs, path_save, df_exper, i_exper)

end

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text:
function inferExact(nodepot, edgepot, es::EdgeStruct)

  nnodes = size(nodepot,1)
  nedges = size(edgepot,3)

  nodebel = zeros(size(nodepot))
  edgebel = zeros(size(edgepot))
  y = ones(nnodes)
  Z = 0
  i = 1
  while true

      pot = configPotential(y,nodepot,edgepot,es)

      for n = 1:nnodes
          nodebel[n,y[n]] = nodebel[n,y[n]]+pot;
      end

      for e = 1:nedges
          n1 = es.edgeends[e,1]
          n2 = es.edgeends[e,2]
          edgebel[y[n1],y[n2],e] = edgebel[y[n1],y[n2],e]+pot
      end

      Z += pot

      for i = 1:nnodes
          y[i] += 1
          if y[i] <= es.nstates[i]
              break
          else
              y[i] = 1;
          end
      end

      if  i == nnodes && y[end] == 1
          break;
      end
  end

  nodebel = nodebel./Z;
  edgebel = edgebel./Z;
  logZ = log(Z);

  return nodebel, edgebel, logZ, Z
end

ext: jl


text:
using UpdateJulia
using Documenter

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

makedocs(;
    modules=[UpdateJulia],
    authors="Lilith Orion Hafner <60898866+LilithHafner@users.noreply.github.com> and contributors",
    repo="https://github.com/LilithHafner/UpdateJulia.jl/blob/{commit}{path}#{line}",
    sitename="UpdateJulia.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://LilithHafner.github.io/UpdateJulia.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/LilithHafner/UpdateJulia.jl",
    devbranch="main",
)


text:
using ThinkJulia

# Exercise 4-1
🐢 = Turtle() # \:turtle:
T = 🐢
# @svg begin
#     forward(T, 100)
#     turn(T, -90)
#     forward(T, 100)
#     turn(T, -90)
#     forward(T, 100)
#     turn(T, -90)
#     forward(T, 100)
# end

# Exercise 4-2, 4-4
function square(t, len = 100)
    for i in 1:4
        forward(t, len)
        turn(t, -90)
    end
end

# Exercise 4-5
function polygon(t, n, len = 100)
    degree = 360/n;
    for i in 1:n
        forward(t, len)
        turn(t, -degree)
    end
end

# Exercise 4-6
function circle(t, r)
    # Approximate by polygon of n-sides
    # circumference = n * len = pi * 2r, better one with tan?
    # len = pi*2r/n
    # len
    #
    # Could be better.

    n = 100 # sides of polygon
    circumference = pi * 2r
    len = circumference / n
    degree = 360/n

    polygon(t, n, len)
end

# Exercise 4-7
function arc(t, r, angle)
    n = 100 # sides of polygon
    circumference = pi * 2r * angle/360
    len = circumference / n
    degree = angle/n

    #polygon(t, n, len)
    for i in 1:n
        forward(t, len)
        turn(t, -degree)
    end
end

@svg begin
    # Exercise 4-1
    # for i in 1:4
    #     forward(T, 100)
    #     turn(T, -90)
    # end

    # Exercise 4-3
    square(T)

    # Exercise 4-4
    square(T, 10)

    # Exercise 4-5
    polygon(T, 7, 80)

    # Exercise 4-6
    circle(T, 30)

    # Exercise 4-7
    arc(T, 120, 270)
end

# TODO
# Exercise 4-8
# Exercise 4-12

# Encapsulation

println("DONE")


text:
# Copies figure 18 to figures directory:

run(`cp ../static/figure_water_period_TTV_trim.png  ..`)
#run(`cp ../static/figure_water_content_versus_period_agol_et_al.png  ../figure_water_period_TTV_trim.png`)


text:
module FrostHelperFlagIfVariantTriggerModule

using ..Ahorn, Maple

@mapdef Trigger "FrostHelper/FlagIfVariantTrigger" FlagIfVariantPlacement(x::Integer, y::Integer, width::Integer=16, height::Integer=16, variant::String="Invincible", flag::String="", variantValue::String="true", inverted::Bool=false)

const placements = Ahorn.PlacementDict(
	"Flag If Variant (Frost Helper)" => Ahorn.EntityPlacement(
		FlagIfVariantPlacement,
		"rectangle"
	)
)

const variants = String[
	"DashAssist",
    "GameSpeed",
    "Hiccups",
    "InfiniteStamina",
    "Invincible",
    "InvisibleMotion",
    "LowFriction",
    "MirrorMode",
    "NoGrabbing",
    "PlayAsBadeline",
    "SuperDashing",
    "ThreeSixtyDashing",
]

const variantValuesBool = String[
	"true", "false"
]

const variantValuesGameSpeed = String[
	"160", "150", "140", "130", "120", "110", "100", "90", "80", "70", "60", "50"
]

function getEditingOptions(entity::FlagIfVariantPlacement)
	variant = get(entity.data, "variant", "Invincible")
	if variant == "GameSpeed"
		return Dict{String, Any}(
    		"variant" => variants,
			"variantValue" => variantValuesGameSpeed
		)
	else
		return Dict{String, Any}(
    		"variant" => variants,
			"variantValue" => variantValuesBool
		)
	end
end

#=
Ahorn.editingOptions(entity::FlagIfVariantPlacement) = Dict{String, Any}(
	"variant" => variants,
	"variantValue" => variantValuesBool
)
=#

Ahorn.editingOptions(entity::FlagIfVariantPlacement) = getEditingOptions(entity)
#=

function Ahorn.editingOptions(entity::FlagIfVariantPlacement)
	variant = Int(get(entity.data, "variant", "Invincible"))
	if variant == "GameSpeed"
		return Dict{String, Any}(
    		"variant" => variants,
			"variantValue" => variantValuesGameSpeed
		)
	else
		return Dict{String, Any}(
    		"variant" => variants,
			"variantValue" => variantValuesBool
		)
	end
end =#

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
@testset "nl_linearity" begin
    @testset "check simplification of formulae" begin
        # First term is true, we should choose `then`
        expr = :(ifelse(1 > 0, x[1] ^ 2, x[2] + 1))
        lin_expr = ReSHOP.LinearityExpr(expr)
        @test ReSHOP.convert_formula(lin_expr.c) == :(x[1] ^ 2)
        @test lin_expr.linearity == :nonlinear

        # First term is false, we should choose `else`
        expr = :(ifelse(1 < 0, x[1] ^ 2, x[2] + 1))
        lin_expr = ReSHOP.LinearityExpr(expr)
        @test ReSHOP.convert_formula(lin_expr.c) == :(x[2] + 1)
        @test lin_expr.linearity == :linear

        # First term isn't constant, we can't simplify
        expr = :(ifelse(1 > x[1], x[1] ^ 2, x[2] + 1))
        lin_expr = ReSHOP.LinearityExpr(expr)
        @test ReSHOP.convert_formula(lin_expr.c) == expr
        @test lin_expr.linearity == :nonlinear
    end
end


text:
export getProfileType, getColMap
export resumeFill, fillComplete
export insertLocalIndices, insertGlobalIndices

#### RowGraph methods ####

getRowMap(graph::CSRGraph) = graph.rowMap
getColMap(graph::CSRGraph) = graph.colMap
getDomainMap(graph::CSRGraph) = graph.domainMap
getRangeMap(graph::CSRGraph) = graph.rangeMap
getImporter(graph::CSRGraph) = graph.importer
getExporter(graph::CSRGraph) = graph.exporter

getGlobalNumRows(graph::CSRGraph) = numGlobalElements(getRowMap(graph))
getGlobalNumCols(graph::CSRGraph) = numGlobalElements(getColMap(graph))
getLocalNumRows(graph::CSRGraph) = numMyElements(getRowMap(graph))
getLocalNumCols(graph::CSRGraph) = numMyElements(getColMap(graph))

getGlobalNumEntries(graph::CSRGraph) = graph.globalNumEntries
getLocalNumEntries(graph::CSRGraph) = graph.nodeNumEntries

Base.@propagate_inbounds function getNumEntriesInGlobalRow(graph::CSRGraph{GID}, globalRow::Integer)::Integer where {GID <: Integer}
    localRow = lid(graph.rowMap, GID(globalRow))
    getNumEntriesInLocalRow(graph, localRow)
end

function getNumEntriesInLocalRow(graph::CSRGraph{GID, PID, LID}, localRow::Integer)::Integer where {GID, PID, LID <: Integer}
    @boundscheck if !hasRowInfo(graph)
        error("Row info is gone, cannot determine number of local entries in the local row $localRow")
    end
    @boundscheck if !myLID(graph.rowMap, LID(localRow))
        error("Row $localRow is not present in this graph")
    end
    info = getRowInfo(graph, LID(localRow))
    retVal = info.numEntries
    recycleRowInfo(info)
    retVal
end

getGlobalNumDiags(graph::CSRGraph) = graph.globalNumDiags
getLocalNumDiags(graph::CSRGraph) = graph.nodeNumDiags

getGlobalMaxNumRowEntries(graph::CSRGraph) = graph.globalMaxNumRowEntries
getLocalMaxNumRowEntries(graph::CSRGraph) = graph.nodeMaxNumRowEntries

hasColMap(graph::CSRGraph) = graph.colMap != nothing

isLowerTriangular(graph::CSRGraph) = graph.lowerTriangle
isUpperTriangular(graph::CSRGraph) = graph.upperTriangle

isGloballyIndexed(graph::CSRGraph) = graph.indicesType == GLOBAL_INDICES
isLocallyIndexed(graph::CSRGraph) = graph.indicesType == LOCAL_INDICES

isFillComplete(g::CSRGraph) = g.fillComplete

function getGlobalRowCopy!(copy::AbstractVector{<:Integer}, graph::CSRGraph{GID}, globalRow::GID) where GID

    rowInfo = getRowInfoFromGlobalRow(myGraph, GID(globalRow))
    viewRange = 1:rowInfo.numEntries

    @boundscheck if rowInfo.localRow == 0
        recycleRowInfo(rowInfo)
        return LID(0)
    end

    numElts = rowInfo.numEntries

    @inbounds if isGloballyIndexed(graph)
        indsView = getGlobalView(myGraph, rowInfo[viewRange])
        for i in LID(1):numElts
            copy[i] = indsView[i]
        end
    else
        colMap = getColMap(myGraph)
        indsView = getLocalView(myGraph, rowInfo)[viewRange]
        for i in LID(1):numElts
            copy[i] = gid(colMap, indsView[i])
        end
    end

    recycleRowInfo(rowInfo)
    numElts
end

function getLocalRowCopy!(copy::AbstractVector{<:Integer},
        graph::CSRGraph{GID, PID, LID}, localRow::LID)::LID where {GID, PID, LID}

    rowInfo = getRowInfo(graph, LID(localRow))
    viewRange = 1:rowInfo.numEntries

    @boundscheck if rowInfo.localRow == 0
        recycleRowInfo(rowInfo)
        return LID(0)
    end

    numElts = rowInfo.numEntries

    @inbounds if isLocallyIndexed(graph)
        indsView = getLocalView(graph, rowInfo[viewRange])
        for i in LID(1):numElts
            copy[i] = indsView[i]
        end
    else
        colMap = getColMap(graph)
        indsView = getGlobalView(graph, rowInfo)[viewRange]

        for i in LID(1):numElts
            copy[i] = lid(colMap, indsView[i])
        end
    end

    recycleRowInfo(rowInfo)
    numElts
end

#### DistObject methods ####
function checkSizes(source::RowGraph{GID, PID, LID},
        target::CSRGraph{GID, PID, LID}) where {GID, PID, LID}
    #T and E petra's don't do any checks
    true
end

function copyAndPermute(source::RowGraph{GID, PID, LID},
        target::CSRGraph{GID, PID, LID}, numSameIDs::LID,
        permuteToLIDs::AbstractArray{LID, 1}, permuteFromLIDs::AbstractArray{LID, 1}) where {
        GID, PID, LID}
    copyAndPermuteNoViewMode(source, target,
        numSameIDs, permuteToLIDs, permuteFromLIDs)
end

function copyAndPermute(source::CSRGraph{GID, PID, LID},
        target::CSRGraph{GID, PID, LID}, numSameIDs::LID,
        permuteToLIDs::AbstractArray{LID, 1}, permuteFromLIDs::AbstractArray{LID, 1}) where {
        GID, PID, LID}
    if isFillComplete(target)
        throw(InvalidStateError("Target cannot be fill complete"))
    end
    if isFillComplete(source)
            copyAndPermuteNoViewMode(source, target,
        numSameIDs, permuteToLIDs, permuteFromLIDs)
    else
        if length(permuteToLIDs) != length(permuteFromLIDs)
            throw(InvalidArgumentError("permuteToLIDs and "
                    * "permuteFromLIDs must have the same size"))
        end

        srcRowMap = getRowMap(source)
        tgtRowMap = getRowMap(target)

        #copy part
        for myid = 1:numSameIDs
            myGID = gid(srcRowMap, myid)
            row = getGlobalRowView(source, myGID)
            insertGlobalIndices(target, myGID, row)
        end

        #permute part
        for i = 1:length(permuteToLIDs)
            srcGID = gid(srcRowMap, permuteFromLIDs[i])
            tgtGID = gid(tgtRowMap, permuteToLIDs[i])
            row = getGlobalRowView(source, srcGID)
            insertGlobalIndices(target, tgtGID, row)
        end
    end
end

function copyAndPermuteNoViewMode(source::RowGraph{GID, PID, LID},
        target::CSRGraph{GID, PID, LID}, numSameIDs::LID,
        permuteToLIDs::AbstractArray{LID, 1}, permuteFromLIDs::AbstractArray{LID, 1}) where {
        GID, PID, LID}
    if length(permuteToLIDs) != length(premuteFromLIDs)
        throw(InvalidArgumentError("permuteToLIDs and "
                * "permuteFromLIDs must have the same size"))
    end

    srcRowMap = getRowMap(source)
    tgtRowMap = getRowMap(target)

    #copy part
    for myid = 1:numSameIDs
        myGID = gid(srcRowMap, myid)
        rowCopy = getGlobalRowCopy(sourceRowGraph, myGID)
        insertGlobalIndices(target, myGID, rowCopy)
    end

    #permute part
    for i = 1:length(permuteToLIDs)
        tgtGID = gid(tgtRowMap, permuteToLIDs[i])
        srcGID = gid(srcRowMap, permuteFromLIDs[i])
        rowCopy = getGlobalRowCopy(source, srcGID)
        insertGlobalIndices(target, tgtGID, rowCopy)
    end
end

function packAndPrepare(source::RowGraph{GID, PID, LID},
        target::CSRGraph{GID, PID, LID}, exportLIDs::AbstractArray{LID, 1},
        distor::Distributor{GID, PID, LID})::Array{GID, 1} where {
        GID, PID, LID}

    pack(source, exportLIDs, distor)
end

function unpackAndCombine(target::CSRGraph{GID, PID, LID},
        importLIDs::AbstractArray{LID, 1}, imports::AbstractArray,
        distor::Distributor{GID, PID, LID}, cm::CombineMode) where {
        GID, PID, LID}
    #should be caught else where
    @assert(isFillActive(target),
        "Import and Export operations require a fill active graph")

    tgtMap = getMap(target)

    for i = 1:length(importLIDs)
        row = imports[i]
        insertGlobalIndicesFiltered(target, gid(tgtMap, importLIDs[i]), row)
    end
end

#### CSRGraph methods ####
#TODO document the CSRGraph methods

"""
    insertLocalIndices(::CSRGraph{GID, PID, LID}, localRow::Integer, [numEntries::Integer,] inds::AbstractArray{<: Integer, 1})

Inserts the given local indices into the graph.  If `numEntries` is given,
only the first `numEntries` elements are inserted
"""
function insertLocalIndices(graph::CSRGraph{GID, PID, LID}, localRow::Integer,
        numEntries::Integer, inds::AbstractArray{<:Integer, 1}) where {GID, PID, LID <: Integer}
    insertLocalIndices(graph, LID(localRow), LID(numEntries), Array{LID, 1}(inds))
end
function insertLocalIndices(graph::CSRGraph{GID, PID, LID}, localRow::LID,
        numEntries::LID, inds::AbstractArray{LID, 1}) where {
        GID, PID, LID <: Integer}
    indicesView = view(inds, 1:numEntries)
    insertLocalIndices(graph, localRow, indsT)
end

function insertLocalIndices(graph::CSRGraph{GID, PID, LID}, localRow::Integer,
        inds::AbstractArray{<:Integer, 1}) where {GID, PID, LID <: Integer}
    insertLocalIndices(graph, LID(localRow), Array{LID, 1}(inds))
end
function insertLocalIndices(graph::CSRGraph{GID, PID, LID},
        localRow::LID, indices::AbstractArray{LID, 1}) where{
        GID, PID, LID <: Integer}
    if !isFillActive(graph)
        throw(InvalidStateError("insertLocalIndices requires that fill is active"))
    end
    if isGloballyIndexed(graph)
        throw(InvalidStateError("graph indices are global, use insertGlobalIndices(...) instead"))
    end
    if !hasColMap(graph)
        throw(InvalidStateError("Cannot insert local indices without a column map"))
    end
    if !myLID(getMap(graph), localRow)
        throw(InvalidArgumentError("Row does not belong to this process"))
    end
    if !hasRowInfo(graph)
        throw(InvalidStateError("Row information was deleted"))
    end

#    colMap = getColMap(graph)
#    badColIndices = [ind for ind in indices if myLID(colMap, ind)]
#    if length(badColIndices) != 0
#        throw(InvalidArgumentError(
#            "Attempting to insert entries in owned row $localRow, "
#            * "at the following column indices: $indices.\n"
#
#            * "Of those, the following indices are not in "
#            * "the column map on this process: $badColIndices.\n"
#
#            * "Since the graph has a column map already, it is "
#            * "invalid to insert entries at those locations"))
#    end

    insertLocalIndicesImpl(graph, localRow, indices)
end

"""
    insertGlobalIndices(::CSRGraph{GID, PID, LID}, localRow::Integer, [numEntries::Integer,] inds::AbstractArray{<: Integer, 1})

Inserts the given global indices into the graph.  If `numEntries` is given,
only the first `numEntries` elements are inserted
"""
function insertGlobalIndices(graph::CSRGraph{GID, PID, LID}, globalRow::Integer,
        numEntries::Integer, inds::AbstractArray{<: Integer, 1}) where {
        GID <: Integer, PID, LID <: Integer}
    insertGlobalIndices(graph, GID(globalRow), LID(numEntries), Array{GID, 1}(inds))
end
function insertGlobalIndices(graph::CSRGraph{GID, PID, LID}, globalRow::GID,
        numEntries::LID, inds::AbstractArray{GID, 1}) where {
        GID <: Integer, PID, LID <: Integer}
    indicesView = view(inds, 1:numEntries)
    insertGlobalIndices(graph, globalRow, indsT)
end
function insertGlobalIndices(graph::CSRGraph{GID, PID, LID}, globalRow::Integer,
        inds::AbstractArray{<: Integer, 1}) where {GID <: Integer, PID, LID <: Integer}
    insertGlobalIndices(graph, GID(globalRow), Array{GID, 1}(inds))
end
function insertGlobalIndices(graph::CSRGraph{GID, PID, LID}, globalRow::GID,
        indices::AbstractArray{GID, 1}) where {GID <: Integer, PID, LID <: Integer}
    if isLocallyIndexed(graph)
        throw(InvalidStateError("Graph indices are local, use insertLocalIndices()"))
    end
    if !hasRowInfo(graph)
        throw(InvalidStateError("Graph row information was deleted"))
    end
    if isFillComplete(graph)
        throw(InvalidStateError("Cannot call this method if the fill is not active"))
    end

    myRow = lid(graph.rowMap, globalRow)
    if myRow != 0
        if false #for debuggin
            if hasColMap(graph)
                colMap = graph.colMap

                badColInds = [index for index in indices
                                    if myGid(colMap, index)==0]
                if length(badColInds) != 0
                    throw(InvalidArgumentError("$(myPid(getComm(graph))): "
                        * "Attempted to insert entries in owned row $globalRow, "
                        * "at the following column indices: $indices.\n"

                        * "Of those, the following indices are not in the "
                        * "column Map on this process: $badColInds.\n"

                        * "Since the matrix has a column map already, it "
                        * "is invalid to insert entries at those locations"))
               end
           end
        end
        insertGlobalIndicesImpl(graph, myRow, indices)
    else
        append!(graph.nonlocalRow, indices)
    end
end

"""
    insertGlobalIndicesFiltered(::CSRGraph{GID, PID, LID}, localRow::Integer, inds::AbstractArray{<: Integer, 1})

As `insertGlobalIndices(...)` but filters by those present in
the column map, if present
"""
function insertGlobalIndicesFiltered(graph::CSRGraph{GID, PID, LID}, globalRow::GID,
        indices::AbstractArray{GID, 1}) where{GID, PID, LID}
    if isLocallyIndexed(graph)
        throw(InvalidStateError(
                "graph indices are local, use insertLocalIndices(...)"))
    end
    if !hasRowInfo(graph)
        throw(InvalidStateError("Graph row information was deleted"))
    end
    if isFillComplete(graph)
        throw(InvalidStateError("Cannot insert into fill complete graph"))
    end

    myRow = lid(graph.rowMap, globalRow)
    if myRow != 0
        #if column map present, use it to filter the entries
        if hasColMap(graph)
            colMap = getColMap(graph)
            indices = [index for index in indices if myLid(colMap, index)]
        end
        insertGlobalIndicesImpl(myRow, indices)
    else
        #nonlocal row
        append!(graph.nonlocals[globalRow], indices)
    end
end

function getGlobalView(graph::CSRGraph{GID, PID, LID}, rowInfo::RowInfo{LID}) where {GID <: Integer, PID, LID <: Integer}
    if rowInfo.allocSize > 0
        if length(graph.globalIndices1D) != 0
            range = rowInfo.offset1D : rowInfo.offset1D + rowInfo.allocSize
            view(graph.globalIndices1D, range)
        elseif length(graph.globalIndices2D[rowInfo.localRow]) == 0
            globalIndices2D[rowInfo.localRow]
        else
            GID[]
        end
    else
        GID[]
    end
end

@inline function getGlobalViewPtr(graph::CSRGraph{GID, PID, LID}, rowInfo::RowInfo{LID})::Tuple{Ptr{GID}, LID} where {GID <: Integer, PID, LID <: Integer}
    if rowInfo.allocSize > 0
        if length(graph.globalIndices1D) != 0
            return (pointer(graph.globalIndices1D, rowInfo.offset1D), rowInfo.allocSize)
        elseif length(graph.globalIndices2D[rowInfo.localRow]) != 0
            baseArray = graph.globalIndices2D[rowInfo.localRow]::Vector{GID}
            return (pointer(baseArray), GID(length(baseArray)))
        end
    end
    return (Ptr{GID}(C_NULL), 0)
end

function getLocalView(graph::CSRGraph{GID, PID, LID}, rowInfo::RowInfo{LID}) where {GID <: Integer, PID, LID <: Integer}
    if rowInfo.allocSize > 0
        if length(graph.localIndices1D) != 0
            range = rowInfo.offset1D : rowInfo.offset1D + rowInfo.allocSize-LID(1)
            return view(graph.localIndices1D, range)
        elseif length(graph.localIndices2D[rowInfo.localRow]) != 0
            baseArray = graph.localIndices2D[rowInfo.localRow]
            return view(baseArray, LID(1):LID(length(baseArray)))
        end
    end
    return LID
end

Base.@propagate_inbounds @inline function getLocalViewPtr(graph::CSRGraph{GID, PID, LID}, rowInfo::RowInfo{LID})::Tuple{Ptr{LID}, LID} where {GID <: Integer, PID, LID <: Integer}
    if rowInfo.allocSize > 0
        if length(graph.localIndices1D) != 0
            return (pointer(graph.localIndices1D, rowInfo.offset1D), rowInfo.allocSize)
        elseif length(graph.localIndices2D[rowInfo.localRow]) != 0
            baseArray::Array{LID, 1} = graph.localIndices2D[rowInfo.localRow]
            return (pointer(baseArray), LID(length(baseArray)))
        end
    end
    return (C_NULL, 0)
end

function getGlobalRowView(graph::CSRGraph{GID}, globalRow::GID)::AbstractArray{GID, 1} where {GID <: Integer}
    if isLocallyIndexed(graph)
        throw(InvalidArgumentError("The graph's indices are currently stored as local indices, so a view with global column indices cannot be returned.  Use getGlobalRowCopy(::CSRGraph) instead"))
    end

    rowInfo = getRowInfoFromGlobalRow(graph, globalRow)

    if rowInfo.localRow != 0 && rowInfo.numEntries > 0
        retVal = view(getGlobalView(graph, rowInfo), 1:rowInfo.numEntries)
    else
        retVal = GID[]
    end
    recycleRowInfo(rowInfo)
    retVal
end

Base.@propagate_inbounds function getGlobalRowViewPtr(graph::CSRGraph{GID, PID, LID}, globalRow::GID)::Tuple{Ptr{GID}, LID} where {GID <: Integer, PID <: Integer, LID <: Integer}
    if isLocallyIndexed(graph)
        throw(InvalidArgumentError("The graph's indices are currently stored as local indices, so a view with global column indices cannot be returned.  Use getGlobalRowCopy(::CSRGraph) instead"))
    end

    rowInfo = getRowInfoFromGlobalRow(graph, globalRow)

    isLocalRow = true
    @boundscheck isLocalRow = rowInfo.localRow != 0

    if isLocalRow && rowInfo.numEntries > 0
        retVal = (getGlobalViewPtr(graph, rowInfo)[1], rowInfo.numEntries)
    else
        retVal = (Ptr{GID}(C_NULL), LID(0))
    end
    recycleRowInfo(rowInfo)
    retVal
end

function getLocalRowView(graph::CSRGraph{GID}, localRow::GID)::AbstractArray{GID, 1} where {GID}
    rowInfo = getRowInfoFromLocalRowIndex(graph, localRow)

    retVal = getLocalRowView(graph, rowInfo)
    recycleRowInfo(rowInfo)
    retVal
end

function getLocalRowView(graph::CSRGraph{GID, PID, LID}, rowInfo::RowInfo{LID}
        )::AbstractArray{GID, 1} where {GID, PID, LID}

    if isGloballyIndexed(graph)
        throw(InvalidArgumentError("The graph's indices are currently stored as global indices, so a view with local column indices cannot be returned.  Use getLocalRowCopy(::CSRGraph) instead"))
    end

    if rowInfo.localRow != 0 && rowInfo.numEntries > 0
        view(getLocalView(graph, rowInfo), 1:rowInfo.numEntries)
    else
        LID[]
    end
end

resumeFill(graph::CSRGraph; plist...) = resumeFill(graph, Dict(plist))

function resumeFill(graph::CSRGraph, plist::Dict)
    if !hasRowInfo(graph)
        throw(InvalidStateError("Cannot resume fill of the CSRGraph, "
                * "since the graph's row information was deleted."))
    end

    clearGlobalConstants(graph)
    graph.plist = plist
    graph.lowerTriangle = false
    graph.upperTriangle = false
    graph.indicesAreSorted = true
    graph.noRedundancies = true
    graph.fillComplete = false
end

fillComplete(graph::CSRGraph; plist...) = fillComplete(graph, Dict(plist))

function fillComplete(graph::CSRGraph, plist::Dict)
    if graph.domainMap == nothing
        domMap = graph.rowMap
    else
        domMap = graph.domainMap
    end

    if graph.rangeMap == nothing
        ranMap = graph.colMap
    else
        ranMap = graph.rangeMap
    end

    fillComplete(graph, ranMap, domMap, plist)
end

function fillComplete(graph::CSRGraph{GID, PID, LID},
        domainMap::BlockMap{GID, PID, LID}, rangeMap::BlockMap{GID, PID, LID};
        plist...) where {GID, PID, LID}
    fillComplete(graph, Dict(plist))
end

function fillComplete(graph::CSRGraph{GID, PID, LID}, domainMap::BlockMap{GID, PID, LID}, rangeMap::BlockMap{GID, PID, LID}, plist::Dict) where {GID, PID, LID}
    if !isFillActive(graph) || isFillComplete(graph)
        throw(InvalidStateError("Graph fill state must be active to call fillComplete(...)"))
    end

    numProcs = numProc(getComm(graph))

    assertNoNonlocalInserts = get(plist, :noNonlocalChanges, false)

    mayNeedGlobalAssemble = !assertNoNonlocalInserts && numProcs > 1
    if mayNeedGlobalAssemble
        globalAssemble(graph)
    else
        if numProcs == 1 && length(graph.nonLocals) > 0
            throw(InvalidStateError("Only one process, but nonlocal entries are present"))
        end
    end

    setDomainRangeMaps(graph, domainMap, rangeMap)

    if !hasColMap(graph)
        makeColMap(graph)
    end

    makeIndicesLocal(graph)

    sortAndMergeAllIndices(graph, isSorted(graph), isMerged(graph))

    makeImportExport(graph)
    computeGlobalConstants(graph)
    fillLocalGraph(graph, plist)
    graph.fillComplete(true)

    checkInternalState(graph)
end

function makeColMap(graph::CSRGraph{GID, PID, LID}) where {GID, PID, LID}
    localNumRows = getLocalNumEntries(graph)

    #TODO look at FIXME on line 4898

    error, colMap = __makeColMap(graph, graph.domainMap)

    comm = getComm(graph)
    globalError = maxAll(comm, error)

    if globalError
        Base.error("makeColMap reports an error on at least one process")
    end

    graph.colMap = colMap

    checkInternalState(graph)
end

"""
    isSorted(::CSRGraph)

Whether the indices are sorted
"""
isSorted(graph::CSRGraph) = graph.indicesAreSorted

"""
    isMerged(::CSRGraph)

Whether duplicate column indices in each row have been merged
"""
isMerged(graph::CSRGraph) = graph.noRedundancies

"""
    setAllIndices(graph::CSRGraph{GID, PID, LID}, rowPointers::Array{LID, 1}, columnIndices::AbstractArray{LID, 1})

Sets the graph's data directly, using 1D storage
"""
function setAllIndices(graph::CSRGraph{GID, PID, LID},
        rowPointers::AbstractArray{LID, 1},columnIndices::Array{LID, 1}) where {
        GID, PID, LID <: Integer}

    localNumRows = getLocalNumRows(graph)

    if graph.colMap == nothing
        throw(InvalidStateError("The graph must have a "
                * "column map before calling setAllIndices"))
    end
    if length(rowPointers) != localNumRows + 1
        throw(InvalidArgumentError("length(rowPointers) = $(length(rowPointers)) "
                * "!= localNumRows+1 = $(localNumRows+1)"))
    end

    localNumEntries = rowPointers[localNumRows+1]

    graph.indicesType    = LOCAL_INDICES
    graph.pftype         = STATIC_PROFILE
    graph.localIndices1D = columnIndices
    graph.rowOffsets     = rowPointers
    graph.nodeNumEntries = localNumEntries
    graph.storageStatus  = STORAGE_1D_UNPACKED

    graph.localGraph     = LocalCSRGraph(columnIndices, rowPointers)

    checkInternalState(graph)
end

"""
    isStorageOptimized(::CSRGraph)

Whether the graph's storage is optimized
"""
function isStorageOptimized(graph::CSRGraph)
    length(graph.numRowEntries) == 0 && getLocalNumRows(graph) > 0
end

"""
    getProfileType(::CSRGraph)

Gets the profile type of the graph
"""
getProfileType(graph::CSRGraph) = graph.pftype


text:
export SimpleVectorObservable, stddev

type SimpleVectorObservable <: VectorObservable
  bins :: Vector{Vector{Float64}}
  num :: Int64
  sum :: Vector{Float64}
  sum2 :: Vector{Float64}
end

SimpleVectorObservable() = SimpleVectorObservable(Vector{Float64}[], 0, Float64[], Float64[])

function reset!(obs :: SimpleVectorObservable)
  obs.bins = Vector{Vector{Float64}}[]
  obs.num = 0
  obs.sum = Float64[]
  obs.sum2 = Float64[]
  return obs
end

count(obs::SimpleVectorObservable) = obs.num

function push!(obs :: SimpleVectorObservable, value::Vector)
  if obs.num == 0
    obs.num = 1
    push!(obs.bins, value)
    obs.sum = deepcopy(value)
    obs.sum2 = squared(value)
  else
    push!(obs.bins, value)
    obs.num += 1
    for i in 1:length(obs.bins[1])
      obs.sum[i] += value[i]
      obs.sum2[i] += squared(value[i])
    end
  end
  return obs
end

function mean(obs::SimpleVectorObservable)
  if obs.num > 0
    return obs.sum .* (1.0/obs.num)
  else
    return [NaN]
  end
end

function var(obs::SimpleVectorObservable)
  if obs.num  > 1
    v = (obs.sum2 - squared(obs.sum)/obs.num)/(obs.num-1)
    return map(maxzero, v)
  elseif obs.num == 1
    return fill(Inf, length(obs.sum))
  else
    return [NaN]
  end
end
stddev(obs::SimpleVectorObservable) = sqrt(var(obs))
stderror(obs::SimpleVectorObservable) = sqrt(var(obs)/count(obs))
function confidence_interval(obs::SimpleVectorObservable, confidence_rate :: Real)
  if count(obs) == 0
    return [Inf]
  elseif count(obs) == 1
    return fill(Inf, length(obs.sum))
  end
  q = 0.5 + 0.5confidence_rate
  correction = quantile( TDist(obs.num - 1), q)
  serr = stderror(obs)
  return correction * serr
end

function confidence_interval(obs::SimpleVectorObservable, confidence_rate_symbol::Symbol = :sigma1)
  n = parsesigma(confidence_rate_symbol)
  return confidence_interval(obs, erf(0.5n*sqrt(2.0)))
end

function merge!(obs::SimpleVectorObservable, other::SimpleVectorObservable)
  append!(obs.bins, other.bins)
  obs.num += other.num
  for i in 1:length(obs.bins[i])
    obs.sum[i] += other.sum[i]
    obs.sum2[i] += other.sum2[i]
  end
  return obs
end
merge(lhs::SimpleVectorObservable, rhs::SimpleVectorObservable) = merge!(deepcopy(lhs), rhs)

export SimpleVectorObservableSet
typealias SimpleVectorObservableSet MCObservableSet{SimpleVectorObservable}

function merge!(obs::SimpleVectorObservableSet, other::SimpleVectorObservableSet)
  obs_names = Set(keys(obs))
  union!(obs_names, Set(keys(other)))
  for name in obs_names
    if !haskey(obs, name)
      # in other only
      obs[name] = deepcopy(other[name])
    elseif haskey(other, name)
      # in both
      merge!(obs[name], other[name])

    # else
    # in obs only
    # NOTHING to do
    end
  end
  return obs
end
merge(lhs::SimpleVectorObservableSet, rhs::SimpleVectorObservableSet) = merge!(deepcopy(lhs), rhs)


text:
"""
```
oneletter(residue::String)
```

Function to return a one-letter residue code from the three letter code or residue name. The function is case-insensitive.

### Examples

```julia-repl
julia> oneletter("ALA")
"A"

julia> oneletter("Glutamic acid")
"E"

```

"""
function oneletter(residue::String)
  l = length(residue)
  code = uppercase(residue)
  if l > 1
    if l == 3
      ires = findfirst(r->r.three_letter_code == code, natural_aminoacids)
    else
      ires = findfirst(r->uppercase(r.name) == code, natural_aminoacids)
    end
    code = (ires == nothing ? nothing : natural_aminoacids[ires].one_letter_code)
  end
  return code
end


text:
###
### **arf.h** -- arbitrary-precision floating-point numbers
###

### Types, macros and constants

### Memory management
arbcall"void arf_init(arf_t x)"
arbcall"void arf_clear(arf_t x)"
arbcall"slong arf_allocated_bytes(const arf_t x)"

### Special values
arbcall"void arf_zero(arf_t res)"
arbcall"void arf_one(arf_t res)"
arbcall"void arf_pos_inf(arf_t res)"
arbcall"void arf_neg_inf(arf_t res)"
arbcall"void arf_nan(arf_t res)"
arbcall"int arf_is_zero(const arf_t x)"
arbcall"int arf_is_one(const arf_t x)"
arbcall"int arf_is_pos_inf(const arf_t x)"
arbcall"int arf_is_neg_inf(const arf_t x)"
arbcall"int arf_is_nan(const arf_t x)"
arbcall"int arf_is_inf(const arf_t x)"
arbcall"int arf_is_normal(const arf_t x)"
arbcall"int arf_is_special(const arf_t x)"
arbcall"int arf_is_finite(const arf_t x)"

### Assignment, rounding and conversions
arbcall"void arf_set(arf_t res, const arf_t x)"
arbcall"void arf_set_mpz(arf_t res, const mpz_t x)"
#ni arbcall"void arf_set_fmpz(arf_t res, const fmpz_t x)"
arbcall"void arf_set_ui(arf_t res, ulong x)"
arbcall"void arf_set_si(arf_t res, slong x)"
arbcall"void arf_set_mpfr(arf_t res, const mpfr_t x)"
#ns arbcall"void arf_set_fmpr(arf_t res, const fmpr_t x)"
arbcall"void arf_set_d(arf_t res, double x)"
arbcall"void arf_swap(arf_t x, arf_t y)"
arbcall"void arf_init_set_ui(arf_t res, ulong x)"
arbcall"void arf_init_set_si(arf_t res, slong x)"
arbcall"int arf_set_round(arf_t res, const arf_t x, slong prec, arf_rnd_t rnd)"
arbcall"int arf_set_round_si(arf_t res, slong x, slong prec, arf_rnd_t rnd)"
arbcall"int arf_set_round_ui(arf_t res, ulong x, slong prec, arf_rnd_t rnd)"
arbcall"int arf_set_round_mpz(arf_t res, const mpz_t x, slong prec, arf_rnd_t rnd)"
#ni arbcall"int arf_set_round_fmpz(arf_t res, const fmpz_t x, slong prec, arf_rnd_t rnd)"
arbcall"void arf_set_si_2exp_si(arf_t res, slong m, slong e)"
arbcall"void arf_set_ui_2exp_si(arf_t res, ulong m, slong e)"
#ni arbcall"void arf_set_fmpz_2exp(arf_t res, const fmpz_t m, const fmpz_t e)"
#ni arbcall"int arf_set_round_fmpz_2exp(arf_t res, const fmpz_t x, const fmpz_t e, slong prec, arf_rnd_t rnd)"
#ni arbcall"void arf_get_fmpz_2exp(fmpz_t m, fmpz_t e, const arf_t x)"
#ni arbcall"void arf_frexp(arf_t m, fmpz_t e, const arf_t x)"
#mo arbcall"double arf_get_d(const arf_t x, arf_rnd_t rnd)" # clashes with arf_get_si
#ns arbcall"void arf_get_fmpr(fmpr_t res, const arf_t x)"
arbcall"int arf_get_mpfr(mpfr_t res, const arf_t x, mpfr_rnd_t rnd)"
#ni arbcall"int arf_get_fmpz(fmpz_t res, const arf_t x, arf_rnd_t rnd)"
#mo arbcall"slong arf_get_si(const arf_t x, arf_rnd_t rnd)" # clashes with arf_get_d
#ni arbcall"int arf_get_fmpz_fixed_fmpz(fmpz_t res, const arf_t x, const fmpz_t e)"
#ni arbcall"int arf_get_fmpz_fixed_si(fmpz_t res, const arf_t x, slong e)"
arbcall"void arf_floor(arf_t res, const arf_t x)"
arbcall"void arf_ceil(arf_t res, const arf_t x)"
#ni arbcall"void arf_get_fmpq(fmpq_t res, const arf_t x)"

### Comparisons and bounds
arbcall"int arf_equal(const arf_t x, const arf_t y)"
arbcall"int arf_equal_si(const arf_t x, slong y)"
arbcall"int arf_cmp(const arf_t x, const arf_t y)"
arbcall"int arf_cmp_si(const arf_t x, slong y)"
arbcall"int arf_cmp_ui(const arf_t x, ulong y)"
arbcall"int arf_cmp_d(const arf_t x, double y)"
arbcall"int arf_cmpabs(const arf_t x, const arf_t y)"
arbcall"int arf_cmpabs_ui(const arf_t x, ulong y)"
arbcall"int arf_cmpabs_d(const arf_t x, double y)"
arbcall"int arf_cmpabs_mag(const arf_t x, const mag_t y)"
arbcall"int arf_cmp_2exp_si(const arf_t x, slong e)"
arbcall"int arf_cmpabs_2exp_si(const arf_t x, slong e)"
arbcall"int arf_sgn(const arf_t x)"
arbcall"void arf_min(arf_t res, const arf_t a, const arf_t b)"
arbcall"void arf_max(arf_t res, const arf_t a, const arf_t b)"
arbcall"slong arf_bits(const arf_t x)"
arbcall"int arf_is_int(const arf_t x)"
arbcall"int arf_is_int_2exp_si(const arf_t x, slong e)"
#ni arbcall"void arf_abs_bound_lt_2exp_fmpz(fmpz_t res, const arf_t x)"
#ni arbcall"void arf_abs_bound_le_2exp_fmpz(fmpz_t res, const arf_t x)"
arbcall"slong arf_abs_bound_lt_2exp_si(const arf_t x)"

### Magnitude functions
arbcall"void arf_get_mag(mag_t res, const arf_t x)"
arbcall"void arf_get_mag_lower(mag_t res, const arf_t x)"
arbcall"void arf_set_mag(arf_t res, const mag_t x)"
arbcall"void mag_init_set_arf(mag_t res, const arf_t x)"
arbcall"void mag_fast_init_set_arf(mag_t res, const arf_t x)"
arbcall"void arf_mag_set_ulp(mag_t res, const arf_t x, slong prec)"
arbcall"void arf_mag_add_ulp(mag_t res, const mag_t x, const arf_t y, slong prec)"
arbcall"void arf_mag_fast_add_ulp(mag_t res, const mag_t x, const arf_t y, slong prec)"

### Shallow assignment
arbcall"void arf_init_set_shallow(arf_t z, const arf_t x)"
arbcall"void arf_init_set_mag_shallow(arf_t z, const mag_t x)"
arbcall"void arf_init_neg_shallow(arf_t z, const arf_t x)"
arbcall"void arf_init_neg_mag_shallow(arf_t z, const mag_t x)"

### Random number generation
#ns arbcall"void arf_randtest(arf_t res, flint_rand_t state, slong bits, slong mag_bits)"
#ns arbcall"void arf_randtest_not_zero(arf_t res, flint_rand_t state, slong bits, slong mag_bits)"
#ns arbcall"void arf_randtest_special(arf_t res, flint_rand_t state, slong bits, slong mag_bits)"

### Input and output
arbcall"void arf_debug(const arf_t x)"
#mo arbcall"void arf_print(const arf_t x)" # clashes with Base.print
arbcall"void arf_printd(const arf_t x, slong d)"
#ns arbcall"void arf_fprint(FILE * file, const arf_t x)"
#ns arbcall"void arf_fprintd(FILE * file, const arf_t y, slong d)"
arbcall"char * arf_dump_str(const arf_t x)"
arbcall"int arf_load_str(arf_t x, const char * str)"
#ns arbcall"int arf_dump_file(FILE * stream, const arf_t x)"
#ns arbcall"int arf_load_file(arf_t x, FILE * stream)"

### Addition and multiplication
arbcall"void arf_abs(arf_t res, const arf_t x)"
arbcall"void arf_neg(arf_t res, const arf_t x)"
arbcall"int arf_neg_round(arf_t res, const arf_t x, slong prec, arf_rnd_t rnd)"
arbcall"int arf_add(arf_t res, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_add_si(arf_t res, const arf_t x, slong y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_add_ui(arf_t res, const arf_t x, ulong y, slong prec, arf_rnd_t rnd)"
#ni arbcall"int arf_add_fmpz(arf_t res, const arf_t x, const fmpz_t y, slong prec, arf_rnd_t rnd)"
#ni arbcall"int arf_add_fmpz_2exp(arf_t res, const arf_t x, const fmpz_t y, const fmpz_t e, slong prec, arf_rnd_t rnd)"
arbcall"int arf_sub(arf_t res, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_sub_si(arf_t res, const arf_t x, slong y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_sub_ui(arf_t res, const arf_t x, ulong y, slong prec, arf_rnd_t rnd)"
#ni arbcall"int arf_sub_fmpz(arf_t res, const arf_t x, const fmpz_t y, slong prec, arf_rnd_t rnd)"
arbcall"void arf_mul_2exp_si(arf_t res, const arf_t x, slong e)"
#ni arbcall"void arf_mul_2exp_fmpz(arf_t res, const arf_t x, const fmpz_t e)"
#mo arbcall"int arf_mul(arf_t res, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd)" # defined using #DEFINE in C which doesn't work in Julia
arbcall"int arf_mul_ui(arf_t res, const arf_t x, ulong y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_mul_si(arf_t res, const arf_t x, slong y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_mul_mpz(arf_t res, const arf_t x, const mpz_t y, slong prec, arf_rnd_t rnd)"
#ni arbcall"int arf_mul_fmpz(arf_t res, const arf_t x, const fmpz_t y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_addmul(arf_t z, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_addmul_ui(arf_t z, const arf_t x, ulong y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_addmul_si(arf_t z, const arf_t x, slong y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_addmul_mpz(arf_t z, const arf_t x, const mpz_t y, slong prec, arf_rnd_t rnd)"
#ni arbcall"int arf_addmul_fmpz(arf_t z, const arf_t x, const fmpz_t y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_submul(arf_t z, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_submul_ui(arf_t z, const arf_t x, ulong y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_submul_si(arf_t z, const arf_t x, slong y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_submul_mpz(arf_t z, const arf_t x, const mpz_t y, slong prec, arf_rnd_t rnd)"
#ni arbcall"int arf_submul_fmpz(arf_t z, const arf_t x, const fmpz_t y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_sosq(arf_t res, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd)"

### Summation
#ni arbcall"int arf_sum(arf_t res, arf_srcptr terms, slong len, slong prec, arf_rnd_t rnd)"

### Division
arbcall"int arf_div(arf_t res, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_div_ui(arf_t res, const arf_t x, ulong y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_ui_div(arf_t res, ulong x, const arf_t y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_div_si(arf_t res, const arf_t x, slong y, slong prec, arf_rnd_t rnd)"
arbcall"int arf_si_div(arf_t res, slong x, const arf_t y, slong prec, arf_rnd_t rnd)"
#ni arbcall"int arf_div_fmpz(arf_t res, const arf_t x, const fmpz_t y, slong prec, arf_rnd_t rnd)"
#ni arbcall"int arf_fmpz_div(arf_t res, const fmpz_t x, const arf_t y, slong prec, arf_rnd_t rnd)"
#ni arbcall"int arf_fmpz_div_fmpz(arf_t res, const fmpz_t x, const fmpz_t y, slong prec, arf_rnd_t rnd)"

### Square roots
arbcall"int arf_sqrt(arf_t res, const arf_t x, slong prec, arf_rnd_t rnd)"
arbcall"int arf_sqrt_ui(arf_t res, ulong x, slong prec, arf_rnd_t rnd)"
#ni arbcall"int arf_sqrt_fmpz(arf_t res, const fmpz_t x, slong prec, arf_rnd_t rnd)"
arbcall"int arf_rsqrt(arf_t res, const arf_t x, slong prec, arf_rnd_t rnd)"
arbcall"int arf_root(arf_t res, const arf_t x, ulong k, slong prec, arf_rnd_t rnd)"

### Complex arithmetic
arbcall"int arf_complex_mul(arf_t e, arf_t f, const arf_t a, const arf_t b, const arf_t c, const arf_t d, slong prec, arf_rnd_t rnd)"
arbcall"int arf_complex_mul_fallback(arf_t e, arf_t f, const arf_t a, const arf_t b, const arf_t c, const arf_t d, slong prec, arf_rnd_t rnd)"
arbcall"int arf_complex_sqr(arf_t e, arf_t f, const arf_t a, const arf_t b, slong prec, arf_rnd_t rnd)"

### Low-level methods
#ni arbcall"int _arf_get_integer_mpn(mp_ptr y, mp_srcptr xp, mp_size_t xn, slong exp)"
#ni arbcall"int _arf_set_mpn_fixed(arf_t z, mp_srcptr xp, mp_size_t xn, mp_size_t fixn, int negative, slong prec, arf_rnd_t rnd)"
arbcall"int _arf_set_round_ui(arf_t z, ulong x, int sgnbit, slong prec, arf_rnd_t rnd)"
#ni arbcall"int _arf_set_round_uiui(arf_t z, slong * fix, mp_limb_t hi, mp_limb_t lo, int sgnbit, slong prec, arf_rnd_t rnd)"
#ni arbcall"int _arf_set_round_mpn(arf_t z, slong * exp_shift, mp_srcptr x, mp_size_t xn, int sgnbit, slong prec, arf_rnd_t rnd)"


text:
"""
```
    kernelmatrix!(K::Matrix, κ::Kernel, X::Matrix; obsdim::Integer=2)
    kernelmatrix!(K::Matrix, κ::Kernel, X::Matrix, Y::Matrix; obsdim::Integer=2)
```
In-place version of `kernelmatrix` where pre-allocated matrix `K` will be overwritten with the kernel matrix.
"""
kernelmatrix!

function kernelmatrix!(
        K::Matrix,
        κ::Kernel,
        X::AbstractMatrix;
        obsdim::Int = defaultobs
        )
        @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))"
        if !check_dims(K,X,X,feature_dim(obsdim),obsdim)
            throw(DimensionMismatch("Dimensions of the target array K $(size(K)) are not consistent with X $(size(X))"))
        end
        map!(x->kappa(κ,x),K,pairwise(metric(κ),transform(κ,X,obsdim),dims=obsdim))
end

function kernelmatrix!(
        K::AbstractMatrix,
        κ::Kernel,
        X::AbstractMatrix,
        Y::AbstractMatrix;
        obsdim::Int = defaultobs
        )
        @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))"
        if !check_dims(K,X,Y,feature_dim(obsdim),obsdim)
            throw(DimensionMismatch("Dimensions $(size(K)) of the target array K are not consistent with X ($(size(X))) and Y ($(size(Y)))"))
        end
        map!(x->kappa(κ,x),K,pairwise(metric(κ),transform(κ,X,obsdim),transform(κ,Y,obsdim),dims=obsdim))
end

## Apply kernel on two reals ##
function _kernel(κ::Kernel, x::Real, y::Real)
    _kernel(κ, [x], [y])
end

## Apply kernel on two vectors ##
function _kernel(
        κ::Kernel,
        x::AbstractVector,
        y::AbstractVector;
        obsdim::Int = defaultobs
    )
    @assert length(x) == length(y) "x and y don't have the same dimension!"
    kappa(κ, evaluate(metric(κ),transform(κ,x),transform(κ,y)))
end

"""
```
    kernelmatrix(κ::Kernel, X::Matrix ; obsdim::Int=2)
    kernelmatrix(κ::Kernel, X::Matrix, Y::Matrix; obsdim::Int=2)
```
Calculate the kernel matrix of `X` (and `Y`) with respect to kernel `κ`.
`obsdim=1` means the matrix `X` (and `Y`) has size #samples x #dimension
`obsdim=2` means the matrix `X` (and `Y`) has size #dimension x #samples
"""
kernelmatrix

function kernelmatrix(
        κ::Kernel,
        X::AbstractVector{<:Real};
        obsdim::Int=defaultobs
        )
        kernelmatrix(κ,reshape(X,1,:),obsdim=2)
end

function kernelmatrix(
        κ::Kernel,
        X::AbstractMatrix;
        obsdim::Int = defaultobs
    )
    K = map(x->kappa(κ,x),pairwise(metric(κ),transform(κ,X,obsdim),dims=obsdim))
end

function kernelmatrix(
        κ::Kernel,
        X::AbstractMatrix,
        Y::AbstractMatrix;
        obsdim=defaultobs
    )
    @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))"
    if !check_dims(X,Y,feature_dim(obsdim),obsdim)
        throw(DimensionMismatch("X $(size(X)) and Y $(size(Y)) do not have the same number of features on the dimension : $(feature_dim(obsdim))"))
    end
    _kernelmatrix(κ,X,Y,obsdim)
end

@inline _kernelmatrix(κ,X,Y,obsdim) = map(x->kappa(κ,x),pairwise(metric(κ),transform(κ,X,obsdim),transform(κ,Y,obsdim),dims=obsdim))

"""
```
    kerneldiagmatrix(κ::Kernel, X::Matrix; obsdim::Int=2)
```
Calculate the diagonal matrix of `X` with respect to kernel `κ`
`obsdim=1` means the matrix `X` has size #samples x #dimension
`obsdim=2` means the matrix `X` has size #dimension x #samples
"""
function kerneldiagmatrix(
        κ::Kernel,
        X::AbstractMatrix;
        obsdim::Int = defaultobs
        )
        @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))"
        if obsdim == 1
            @compat eachrow(X) .|> x->_kernel(κ,x,x) #[@views _kernel(κ,X[i,:],X[i,:]) for i in 1:size(X,obsdim)]
        elseif obsdim == 2
            @compat eachcol(X) .|> x->_kernel(κ,x,x) #[@views _kernel(κ,X[:,i],X[:,i]) for i in 1:size(X,obsdim)]
        end
end

"""
```
    kerneldiagmatrix!(K::AbstractVector,κ::Kernel, X::Matrix; obsdim::Int=2)
```
In place version of `kerneldiagmatrix`
"""
function kerneldiagmatrix!(
        K::AbstractVector,
        κ::Kernel,
        X::AbstractMatrix;
        obsdim::Int = defaultobs
        )
        @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))"
        if length(K) != size(X,obsdim)
            throw(DimensionMismatch("Dimensions of the target array K $(size(K)) are not consistent with X $(size(X))"))
        end
        if obsdim == 1
            for i in eachindex(K)
                @inbounds @views K[i] = _kernel(κ, X[i,:],X[i,:])
            end
        else
            for i in eachindex(K)
                @inbounds @views K[i] = _kernel(κ,X[:,i],X[:,i])
            end
        end
        return K
end


text: Dict{String,Any}("million" => Dict{String,Any}("value" => "1000000","type" => "integer"))

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
"Holds the tableau of a diagonally implicit Runge-Kutta method."
struct TableauDIRK{T} <: AbstractTableauIRK{T}
    @HeaderTableau

    q::CoefficientsRK{T}

    function TableauDIRK{T}(q) where {T}
        @assert istril(q.a)
        @assert !(q.s==1 && q.a[1,1] ≠ 0)

        if q.s > 1 && istrilstrict(q.a)
            @warn "Initializing TableauDIRK with explicit tableau $(q.name).\nYou might want to use TableauERK instead."
        end

        new(q.name, q.o, q.s, q)
    end
end

function TableauDIRK(q::CoefficientsRK{T}) where {T}
    TableauDIRK{T}(q)
end

function TableauDIRK(name::Symbol, order::Int, a::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
    TableauDIRK{T}(CoefficientsRK(name, order, a, b, c))
end

# TODO function readTableauDIRKFromFile(dir::AbstractString, name::AbstractString)

"Parameters for right-hand side function of diagonally implicit Runge-Kutta methods."
mutable struct ParametersDIRK{DT, TT, ET <: ODE{DT,TT}, D, S} <: Parameters{DT,TT}
    equ::ET
    tab::TableauDIRK{TT}
    Δt::TT

    t::TT
    q::Vector{DT}
    V::Vector{Vector{DT}}
end

function ParametersDIRK(equ::ET, tab::TableauDIRK{TT}, Δt::TT) where {DT, TT, ET <: ODE{DT,TT}}
    ParametersDIRK{DT, TT, ET, equ.d, tab.s}(equ, tab, Δt, 0, zeros(DT, equ.d), create_internal_stage_vector(DT, equ.d, tab.s))
end

"""
Diagonally implicit Runge-Kutta integrator cache.
"""
struct IntegratorCacheDIRK{DT,D,S} <: ODEIntegratorCache{DT,D}
    q̃::Vector{DT}
    ṽ::Vector{DT}
    s̃::Vector{DT}
    y::Vector{DT}

    Q::Vector{Vector{DT}}
    V::Vector{Vector{DT}}
    Y::Vector{Vector{DT}}

    function IntegratorCacheDIRK{DT,D,S}() where {DT,D,S}
        Q = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        Y = create_internal_stage_vector(DT, D, S)
        new(zeros(DT,D), zeros(DT,D), zeros(DT,D), zeros(DT,D), Q, V, Y)
    end
end

"Diagonally implicit Runge-Kutta integrator."
struct IntegratorDIRK{DT, TT, PT <: ParametersDIRK{DT,TT},
                              ST,# <: NonlinearSolver{DT},
                              IT <: InitialGuessODE{DT,TT}, N, D, S} <: IntegratorRK{DT,TT}
    params::PT
    solver::ST
    iguess::IT
    cache::IntegratorCacheDIRK{DT,D,S}
end

function IntegratorDIRK(equation::ODE{DT,TT,FT,N}, tableau::TableauDIRK{TT}, Δt::TT) where {DT,TT,FT,N}
    D = equation.d
    M = equation.n
    S = tableau.q.s

    # create params
    params = ParametersDIRK(equation, tableau, Δt)

    # create solvers
    function create_nonlinear_solver(DT, N, params, i)
        # create solution vector for nonlinear solver
        x = zeros(DT, N)

        # create wrapper function f(x,b) that calls `function_stages!(x, b, params)`
        # with the appropriate `params`
        f = (x,b) -> function_stages!(x, b, params, i)

        # create nonlinear solver with solver type obtained from config dictionary
        s = get_config(:nls_solver)(x, f)
    end

    solvers = [create_nonlinear_solver(DT, D, params, i) for i in 1:S]

    # create initial guess
    iguess = InitialGuessODE(get_config(:ig_interpolation), equation, Δt)

    # create cache
    cache = IntegratorCacheDIRK{DT,D,S}()

    # create integrator
    IntegratorDIRK{DT, TT, typeof(params), typeof(solvers), typeof(iguess), N, D, S}(
                params, solvers, iguess, cache)

end

has_initial_guess(int::IntegratorDIRK) = true

"Initialise initial guess"
function initialize!(int::IntegratorDIRK, cache::IntegratorCacheDIRK)
    # initialise initial guess
    cache.t̅ = cache.t - timestep(int)

    int.params.equ.v(cache.t, cache.q, cache.v)

    initialize!(int.iguess, cache.t, cache.q, cache.v, cache.t̅, cache.q̅, cache.v̅)
end

"Compute initial guess for internal stages."
function initial_guess!(int::IntegratorDIRK, sol::AtomicSolutionODE)
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q, sol.v, sol.q̅, sol.v̅, int.cache.q̃, int.cache.ṽ, int.params.tab.q.c[i])
        for k in eachindex(int.cache.V[i], int.cache.ṽ)
            int.cache.V[i][k] = int.cache.ṽ[k]
        end
    end
    for i in eachstage(int)
        for k in eachdim(int)
            int.solver[i].x[k] = 0
            for j in eachstage(int)
                int.solver[i].x[k] += int.params.tab.q.a[i,j] * int.cache.V[j][k]
            end
        end
    end
end

function compute_stages!(x::Vector{ST}, Q::Vector{ST}, V::Vector{ST}, Y::Vector{ST},
                                    params::ParametersDIRK{DT,TT,ET,D,S}, i::Int) where {ST,DT,TT,ET,D,S}

    local tᵢ::TT

    @assert D == length(Q) == length(V) == length(Y)

    # copy x to Y and compute Q = q + Δt Y
    for k in 1:D
        Y[k] = x[k]
        Q[k] = params.q[k] + params.Δt * Y[k]
    end

    # compute V = v(Q)
    tᵢ = params.t + params.Δt * params.tab.q.c[i]
    params.equ.v(tᵢ, Q, V)
end

"Compute stages of fully implicit Runge-Kutta methods."
function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersDIRK{DT,TT,ET,D,S}, i::Int) where {ST,DT,TT,ET,D,S}
    local Q = zeros(ST,D)
    local V = zeros(ST,D)
    local Y = zeros(ST,D)

    local y1::ST
    local y2::ST

    compute_stages!(x, Q, V, Y, params, i)

    # compute b = - (Y-AV)
    for k in 1:D
        y1 = params.tab.q.a[i,i] * V[k]
        y2 = params.tab.q.â[i,i] * V[k]
        for j in 1:i-1
            y1 += params.tab.q.a[i,j] * params.V[j][k]
            y2 += params.tab.q.â[i,j] * params.V[j][k]
        end
        b[k] = - Y[k] + (y1 + y2)
    end
end

"Integrate ODE with diagonally implicit Runge-Kutta integrator."
function integrate_step!(int::IntegratorDIRK{DT,TT}, sol::AtomicSolutionODE{DT,TT}) where {DT,TT,N}
    int.params.t  = sol.t
    int.params.q .= sol.q

    # compute initial guess
    initial_guess!(int, sol)

    # reset cache
    reset!(sol, timestep(int))

    # consecutively solve for all stages
    for i in eachstage(int)
        # call nonlinear solver
        solve!(int.solver[i])

        # print solver status
        print_solver_status(int.solver[i].status, int.solver[i].params)

        # check if solution contains NaNs or error bounds are violated
        check_solver_status(int.solver[i].status, int.solver[i].params)

        # compute vector field at internal stages
        compute_stages!(int.solver[i].x, int.cache.Q[i], int.cache.V[i], int.cache.Y[i], int.params, i)

        # copy velocity field to parameters
        int.params.V[i] .= int.cache.V[i]
    end

    # compute final update
    update_solution!(sol.q, sol.q̃, int.cache.V, int.params.tab.q.b, int.params.tab.q.b̂, int.params.Δt)

    # update vector field for initial guess
    update!(int.iguess, sol.t, sol.q, sol.v)
end


text:
"""
    Wishart(ν, S)
```julia
ν::Real           degrees of freedom (whole number or a real number greater than p - 1)
S::AbstractPDMat  p x p scale matrix
```
The [Wishart distribution](http://en.wikipedia.org/wiki/Wishart_distribution)
generalizes the gamma distribution to ``p\\times p`` real, positive semidefinite
matrices ``\\mathbf{H}``.

If ``\\nu>p-1``, then ``\\mathbf{H}\\sim \\textrm{W}_p(\\nu, \\mathbf{S})``
has rank ``p`` and its probability density function is

```math
f(\\mathbf{H};\\nu,\\mathbf{S}) = \\frac{1}{2^{\\nu p/2} \\left|\\mathbf{S}\\right|^{\\nu/2} \\Gamma_p\\left(\\frac {\\nu}{2}\\right ) }{\\left|\\mathbf{H}\\right|}^{(\\nu-p-1)/2} e^{-(1/2)\\operatorname{tr}(\\mathbf{S}^{-1}\\mathbf{H})}.
```

If ``\\nu\\leq p-1``, then ``\\mathbf{H}`` is rank ``\\nu`` and it has
a density with respect to a suitably chosen volume element on the space of
positive semidefinite matrices. See [here](https://doi.org/10.1214/aos/1176325375).

For integer ``\\nu``, a random matrix given by

```math
\\mathbf{H} = \\mathbf{X}\\mathbf{X}^{\\rm{T}},
\\quad\\mathbf{X} \\sim \\textrm{MN}_{p,\\nu}(\\mathbf{0}, \\mathbf{S}, \\mathbf{I}_{\\nu})
```

has ``\\mathbf{H}\\sim \\textrm{W}_p(\\nu, \\mathbf{S})``.
For non-integer ``\\nu``, Wishart matrices can be generated via the
[Bartlett decomposition](https://en.wikipedia.org/wiki/Wishart_distribution#Bartlett_decomposition).
"""
struct Wishart{T<:Real, ST<:AbstractPDMat, R<:Integer} <: ContinuousMatrixDistribution
    df::T          # degree of freedom
    S::ST          # the scale matrix
    logc0::T       # the logarithm of normalizing constant in pdf
    rank::R        # rank of a sample
    singular::Bool # singular of nonsingular wishart?
end

#  -----------------------------------------------------------------------------
#  Constructors
#  -----------------------------------------------------------------------------

function Wishart(df::T, S::AbstractPDMat{T}) where T<:Real
    df > 0 || throw(ArgumentError("df must be positive. got $(df)."))
    p = dim(S)
    singular = df <= p - 1
    if singular
        isinteger(df) || throw(
            ArgumentError("df of a singular Wishart distribution must be an integer (got $df)")
        )
    end
    rnk::Integer = ifelse(singular, df, p)
    logc0 = wishart_logc0(df, S, rnk)
    _df, _logc0 = promote(df, logc0)
    Wishart{typeof(_df), typeof(S), typeof(rnk)}(_df, S, _logc0, rnk, singular)
end

function Wishart(df::Real, S::AbstractPDMat)
    T = Base.promote_eltype(df, S)
    Wishart(T(df), convert(AbstractArray{T}, S))
end

Wishart(df::Real, S::Matrix) = Wishart(df, PDMat(S))
Wishart(df::Real, S::Cholesky) = Wishart(df, PDMat(S))

#  -----------------------------------------------------------------------------
#  REPL display
#  -----------------------------------------------------------------------------

show(io::IO, d::Wishart) = show_multline(io, d, [(:df, d.df), (:S, d.S)])

#  -----------------------------------------------------------------------------
#  Conversion
#  -----------------------------------------------------------------------------

function convert(::Type{Wishart{T}}, d::Wishart) where T<:Real
    P = convert(AbstractArray{T}, d.S)
    Wishart{T, typeof(P), typeof(d.rank)}(T(d.df), P, T(d.logc0), d.rank, d.singular)
end
function convert(::Type{Wishart{T}}, df, S::AbstractPDMat, logc0, rnk, singular) where T<:Real
    P = convert(AbstractArray{T}, S)
    Wishart{T, typeof(P), typeof(rnk)}(T(df), P, T(logc0), rnk, singular)
end

#  -----------------------------------------------------------------------------
#  Properties
#  -----------------------------------------------------------------------------

insupport(::Type{Wishart}, X::AbstractMatrix) = ispossemdef(X)
function insupport(d::Wishart, X::AbstractMatrix)
    size(X) == size(d) || return false
    if d.singular
        return ispossemdef(X, rank(d))
    else
        return isposdef(X)
    end
end

dim(d::Wishart) = dim(d.S)
size(d::Wishart) = (p = dim(d); (p, p))
rank(d::Wishart) = d.rank
params(d::Wishart) = (d.df, d.S)
@inline partype(d::Wishart{T}) where {T<:Real} = T

mean(d::Wishart) = d.df * Matrix(d.S)

function mode(d::Wishart)
    r = d.df - dim(d) - 1
    r > 0 || throw(ArgumentError("mode is only defined when df > p + 1"))
    return Matrix(d.S) * r
end

function meanlogdet(d::Wishart)
    logdet_S = logdet(d.S)
    p = dim(d)
    v = logdet_S + p * oftype(logdet_S, logtwo)
    df = oftype(logdet_S, d.df)
    for i in 0:(p - 1)
        v += digamma((df - i) / 2)
    end
    return d.singular ? oftype(v, -Inf) : v
end

function entropy(d::Wishart)
    d.singular && throw(ArgumentError("entropy not defined for singular Wishart."))
    p = dim(d)
    df = d.df
    return -d.logc0 - ((df - p - 1) * meanlogdet(d) - df * p) / 2
end

#  Gupta/Nagar (1999) Theorem 3.3.15.i
function cov(d::Wishart, i::Integer, j::Integer, k::Integer, l::Integer)
    S = d.S
    return d.df * (S[i, k] * S[j, l] + S[i, l] * S[j, k])
end

function var(d::Wishart, i::Integer, j::Integer)
    S = d.S
    return d.df * (S[i, i] * S[j, j] + S[i, j] ^ 2)
end

#  -----------------------------------------------------------------------------
#  Evaluation
#  -----------------------------------------------------------------------------

function wishart_logc0(df::T, S::AbstractPDMat{T}, rnk::Integer) where {T<:Real}
    p = dim(S)
    if df <= p - 1
        return singular_wishart_logc0(p, df, S, rnk)
    else
        return nonsingular_wishart_logc0(p, df, S)
    end
end

function logkernel(d::Wishart, X::AbstractMatrix)
    if d.singular
        return singular_wishart_logkernel(d, X)
    else
        return nonsingular_wishart_logkernel(d, X)
    end
end

#  Singular Wishart pdf: Theorem 6 in Uhlig (1994 AoS)
function singular_wishart_logc0(p::Integer, df::T, S::AbstractPDMat{T}, rnk::Integer) where {T<:Real}
    logdet_S = logdet(S)
    h_df = oftype(logdet_S, df) / 2
    return -h_df * (logdet_S + p * oftype(logdet_S, logtwo)) - logmvgamma(rnk, h_df) + ((rnk * (rnk - p)) // 2) * oftype(logdet_S, logπ)
end

function singular_wishart_logkernel(d::Wishart, X::AbstractMatrix)
    p = dim(d)
    r = rank(d)
    L = eigvals(Hermitian(X), (p - r + 1):p)
    return ((d.df - (p + 1)) * sum(log, L) - tr(d.S \ X)) / 2
end

#  Nonsingular Wishart pdf
function nonsingular_wishart_logc0(p::Integer, df::T, S::AbstractPDMat{T}) where {T<:Real}
    logdet_S = logdet(S)
    h_df = oftype(logdet_S, df) / 2
    return -h_df * (logdet_S + p * oftype(logdet_S, logtwo)) - logmvgamma(p, h_df)
end

function nonsingular_wishart_logkernel(d::Wishart, X::AbstractMatrix)
    return ((d.df - (dim(d) + 1)) * logdet(cholesky(X)) - tr(d.S \ X)) / 2
end

#  -----------------------------------------------------------------------------
#  Sampling
#  -----------------------------------------------------------------------------

function _rand!(rng::AbstractRNG, d::Wishart, A::AbstractMatrix)
    if d.singular
        axes2 = axes(A, 2)
        r = rank(d)
        randn!(rng, view(A, :, axes2[1:r]))
        fill!(view(A, :, axes2[(r + 1):end]), zero(eltype(A)))
    else
        _wishart_genA!(rng, A, d.df)
    end
    unwhiten!(d.S, A)
    A .= A * A'
end

function _wishart_genA!(rng::AbstractRNG, A::AbstractMatrix, df::Real)
    # Generate the matrix A in the Bartlett decomposition
    #
    #   A is a lower triangular matrix, with
    #
    #       A(i, j) ~ sqrt of Chisq(df - i + 1) when i == j
    #               ~ Normal()                  when i > j
    #
    T = eltype(A)
    z = zero(T)
    axes1 = axes(A, 1)
    @inbounds for (j, jdx) in enumerate(axes(A, 2)), (i, idx) in enumerate(axes1)
        A[idx, jdx] = if i < j
            z
        elseif i > j
            randn(rng, T)
        else
            rand(rng, Chi(df - i + 1))
        end
    end
    return A
end

#  -----------------------------------------------------------------------------
#  Test utils
#  -----------------------------------------------------------------------------

function _univariate(d::Wishart)
    check_univariate(d)
    df, S = params(d)
    α = df / 2
    β = 2 * first(S)
    return Gamma(α, β)
end

function _rand_params(::Type{Wishart}, elty, n::Int, p::Int)
    n == p || throw(ArgumentError("dims must be equal for Wishart"))
    ν = elty(n - 1 + abs(10 * randn()))
    X = rand(elty, n, n)
    X .= 2 .* X .- 1
    S = X * X'
    return ν, S
end


text:
"""
`LinearKernel([ρ=1.0,[c=0.0]])`
The linear kernel is a Mercer kernel given by
```
    κ(x,y) = ρ²xᵀy + c
```
Where `c` is a real number
"""
struct LinearKernel{T,Tr,Tc<:Real} <: Kernel{T,Tr}
    transform::Tr
    metric::DotProduct
    c::Tc
    function LinearKernel{T,Tr,Tc}(transform::Tr,c::Tc) where {T,Tr<:Transform,Tc<:Real}
        return new{T,Tr,Tc}(transform,DotProduct(),c)
    end
end

function LinearKernel(ρ::T₁=1.0,c::T₂=zero(T₁)) where {T₁<:Real,T₂<:Real}
    LinearKernel{T₁,ScaleTransform{T₁},T₂}(ScaleTransform(ρ),c)
end

function LinearKernel(ρ::AbstractVector{T₁},c::T₂=zero(T₁)) where {T₁<:Real,T₂<:Real}
    LinearKernel{T₁,ARDTransform{T₁,length(ρ)},T₂}(ARDTransform(ρ),c)
end

function LinearKernel(t::Tr,c::T=zero(Float64)) where {Tr<:Transform,T<:Real}
    LinearKernel{eltype(t),Tr,T}(t,c)
end

params(k::LinearKernel) = (params(transform(k)),k.c)
opt_params(k::LinearKernel) = (opt_params(transform(k)),k.c)

@inline kappa(κ::LinearKernel, xᵀy::T) where {T<:Real} = xᵀy + κ.c

"""
`PolynomialKernel([ρ=1.0[,d=2.0[,c=0.0]]])`
The polynomial kernel is a Mercer kernel given by
```
    κ(x,y) = (ρ²xᵀy + c)^d
```
Where `c` is a real number, and `d` is a shape parameter bigger than 1
"""
struct PolynomialKernel{T,Tr,Tc<:Real,Td<:Real} <: Kernel{T,Tr}
    transform::Tr
    metric::DotProduct
    c::Tc
    d::Td
    function PolynomialKernel{T,Tr,Tc,Td}(transform::Tr,c::Tc,d::Td) where {T,Tr<:Transform,Tc<:Real,Td<:Real}
        return new{T,Tr,Tc,Td}(transform,DotProduct(),c,d)
    end
end

function PolynomialKernel(ρ::T₁=1.0,d::T₂=2.0,c::T₃=zero(T₁)) where {T₁<:Real,T₂<:Real,T₃<:Real}
    @check_args(PolynomialKernel, d, d >= one(T₁), "d >= 1")
    PolynomialKernel{T₁,ScaleTransform{T₁},T₂,T₃}(ScaleTransform(ρ),c,d)
end

function PolynomialKernel(ρ::AbstractVector{T₁},d::T₂=2.0,c::T₃=zero(T₁)) where {T₁<:Real,T₂<:Real,T₃<:Real}
    @check_args(PolynomialKernel, d, d >= one(T₂), "d >= 1")
    PolynomialKernel{T₁,ARDTransform{T₁,length(ρ)},T₂,T₃}(ARDTransform(ρ),c,d)
end

function PolynomialKernel(t::Tr,d::T₁=2.0,c::T₂=zero(eltype(T₁))) where {Tr<:Transform,T₁<:Real,T₂<:Real}
    @check_args(PolynomialKernel, d, d >= one(T₁), "d >= 1")
    PolynomialKernel{eltype(Tr),Tr,T₁,T₂}(t,c,d)
end

params(k::PolynomialKernel) = (params(transform(k)),k.d,k.c)
opt_params(k::PolynomialKernel) = (opt_params(transform(k)),k.d,k.c)

@inline kappa(κ::PolynomialKernel, xᵀy::T) where {T<:Real} = (xᵀy + κ.c)^(κ.d)


text:
function dft(x::Array{ArbComplex{P},1}) where {P}
    # to guarantee safe bounds
    length = size(x)[1]

    # we transfer the data into a C acb vector

    #(ALLOC x_c) Allocate the memory in C space
    x_c = ccall(@libarb(_acb_vec_init), Ptr{ArbComplex{P}}, (Int32,), length)
    #

    # copy x_c
    @inbounds for i in 1:length
                unsafe_store!(x_c, x[i], i)
              end
    #

    #(ALLOC transf_C)
    transf_c = ccall(@libarb(_acb_vec_init), Ptr{ArbComplex{P}}, (Int32,), length)
    #

    # we call the acb_dft
    # Signature: void acb_dft(acb_ptr w, acb_srcptr v, slong n, slong prec)
    ccall(@libarb(acb_dft), Cvoid, (Ref{ArbComplex{P}}, Ref{ArbComplex{P}}, Cint, Cint ), transf_c, x_c, length, P)

    #(DEALLOC x_c)
    ccall(@libarb(_acb_vec_clear), Cvoid, (Ref{ArbComplex{P}}, Cint ), x_c, length)
    #

    transf = zeros(ArbComplex{P}, length)
    @inbounds for i in 1:length
                transf[i] = unsafe_load(transf_c, i)
              end

    #(DEALLOC transf_c)
    ccall(@libarb(_acb_vec_clear), Cvoid, (Ref{ArbComplex{P}}, Cint ), transf_c, length)
    #

    return transf
end

dft(x::Array{ArbComplex,1}) = dft(Array{ArbComplex{workingprecision(ArbComplex)}, 1}(x))
dft(x::Array{ArbReal{P},1}) where {P} = dft(Array{ArbComplex{P}, 1}(x))
dft(x::Array{ArbFloat{P},1}) where {P} = dft(Array{ArbComplex{P}, 1}(x))
dft(x::Array{ArbReal,1}) = dft(Array{ArbComplex{workingprecision(ArbReal)}, 1}(x))
dft(x::Array{ArbFloat,1}) = dft(Array{ArbComplex{workingprecision(ArbFloat)}, 1}(x))

function inverse_dft(x::Array{ArbComplex{P},1}) where {P}
    # to guarantee safe bounds
    length = size(x)[1]

    # we transfer the data into a C acb vector

    #(ALLOC x_c) Allocate the memory in C space
    x_c = ccall(@libarb(_acb_vec_init), Ptr{ArbComplex{P}}, (Int32,), length)
    #

    # copy x_c
    @inbounds for i in 1:length
                unsafe_store!(x_c, x[i], i)
              end
    #

    #(ALLOC transf_C)
    transf_c = ccall(@libarb(_acb_vec_init), Ptr{ArbComplex{P}}, (Int32,), length)
    #

    # we call the acb_dft
    # Signature: void acb_dft(acb_ptr w, acb_srcptr v, slong n, slong prec)
    ccall(@libarb(acb_dft_inverse), Cvoid, (Ref{ArbComplex{P}}, Ref{ArbComplex{P}}, Cint, Cint ), transf_c, x_c, length, P)

    #(DEALLOC x_c)
    ccall(@libarb(_acb_vec_clear), Cvoid, (Ref{ArbComplex{P}}, Cint ), x_c, length)
    #

    transf = zeros(ArbComplex{P}, length)
    @inbounds for i in 1:length
                transf[i] = unsafe_load(transf_c, i)
              end

    #(DEALLOC transf_c)
    ccall(@libarb(_acb_vec_clear), Cvoid, (Ref{ArbComplex{P}}, Cint ), transf_c, length)
    #

    return transf
end

inverse_dft(x::Array{ArbComplex,1}) = inverse_dft(Array{ArbComplex{workingprecision(ArbComplex)}, 1}(x))
inverse_dft(x::Array{ArbReal{P},1}) where {P} = dft(Array{ArbComplex{P}, 1}(x))
inverse_dft(x::Array{ArbFloat{P},1}) where {P} = dft(Array{ArbComplex{P}, 1}(x))
inverse_dft(x::Array{ArbReal,1}) = dft(Array{ArbComplex{workingprecision(ArbReal)}, 1}(x))
inverse_dft(x::Array{ArbFloat,1}) = dft(Array{ArbComplex{workingprecision(ArbFloat)}, 1}(x))


text:
module SandboxTests
import ..Pkg # ensure we are using the correct Pkg

# Order-dependence in the tests, so we delay this until we need it
if Base.find_package("Preferences") === nothing
    @info "Installing Preferences for Pkg tests"
    Pkg.add("Preferences") # Needed for sandbox and artifacts tests
end

using Test
using UUIDs
using Pkg
using Preferences

using ..Utils
test_test(fn, name; kwargs...) = Pkg.test(name; test_fn=fn, kwargs...)
test_test(fn; kwargs...)       = Pkg.test(;test_fn=fn, kwargs...)

@testset "Basic `test` sandboxing" begin
    # also indirectly checks that test `compat` is obeyed
    temp_pkg_dir() do project_path; mktempdir() do tmp
        copy_test_package(tmp, "SandboxFallback2")
        proj = joinpath(tmp, "SandboxFallback2")
        Pkg.activate(proj)
        withenv("JULIA_PROJECT" => proj) do; test_test("Unregistered") do
            json = get(Pkg.Types.Context().env.manifest, UUID("682c06a0-de6a-54ab-a142-c8b1cf79cde6"), nothing)
            @test json !== nothing
            @test json.version == v"0.20.0"
            # test that the active project is the tmp one even though
            # JULIA_PROJECT might be set
            @test !haskey(ENV, "JULIA_PROJECT")
            @test Base.active_project() != proj
            @test Base.LOAD_PATH[1] == "@"
            @test startswith(Base.active_project(), Base.LOAD_PATH[2])
        end end
    end end
    # test dependencies should be preserved, when possible
    temp_pkg_dir() do project_path; mktempdir() do tmp
        copy_test_package(tmp, "Sandbox_PreserveTestDeps")
        Pkg.activate(joinpath(tmp, "Sandbox_PreserveTestDeps"))
        test_test("Foo") do
            x = get(Pkg.Types.Context().env.manifest, UUID("7876af07-990d-54b4-ab0e-23690620f79a"), nothing)
            @test x !== nothing
            @test x.version == v"0.4.0"
        end
    end end
end

@testset "Preferences sandboxing without test/Project.toml" begin
    # Preferences should be copied over into sandbox
    temp_pkg_dir() do project_path; mktempdir() do tmp
        copy_test_package(tmp, "Sandbox_PreservePreferences")
        Pkg.activate(joinpath(tmp, "Sandbox_PreservePreferences"))
        test_test() do
            uuid =  UUID("3872bf94-3adb-11e9-01dc-bf80c7641364")
            @test !Preferences.has_preference(uuid, "does_not_exist")
            @test Preferences.load_preference(uuid, "tree") == "birch"
            @test Preferences.load_preference(uuid, "default") === nothing
        end
    end end
end

@testset "Preferences sandboxing with test/Project.toml" begin
    # Preferences should be copied over into sandbox
    temp_pkg_dir() do project_path; mktempdir() do tmp
        copy_test_package(tmp, "Sandbox_PreservePreferences")
        Pkg.activate(joinpath(tmp, "Sandbox_PreservePreferences"))

        # Create fake test/Project.toml and test/LocalPreferences.toml
        open(joinpath(tmp, "Sandbox_PreservePreferences", "test", "Project.toml"), write=true) do io
            print(io, """
            [deps]
            Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
            """)
        end
        Preferences.set_preferences!(
            joinpath(tmp, "Sandbox_PreservePreferences", "test", "LocalPreferences.toml"),
            "Sandbox_PreservePreferences",
            "scent" => "juniper",
        )

        # This test should have completely different preferences
        test_test() do
            uuid =  UUID("3872bf94-3adb-11e9-01dc-bf80c7641364")
            @test !Preferences.has_preference(uuid, "does_not_exist")
            # `tree` has no mapping because we do not know anything about the project
            # in the root directory; test environments are not stacked!
            @test Preferences.load_preference(uuid, "tree") === nothing
            @test Preferences.load_preference(uuid, "scent") == "juniper"
            @test Preferences.load_preference(uuid, "default") === nothing
        end
    end end
end

@testset "Nested Preferences sandboxing" begin
    # Preferences should be copied over into sandbox
    temp_pkg_dir() do project_path; mktempdir() do tmp
        copy_test_package(tmp, "Sandbox_PreservePreferences")
        Pkg.activate(joinpath(tmp, "Sandbox_PreservePreferences"))
        test_test("Foo") do
            uuid =  UUID("48898bec-3adb-11e9-02a6-a164ba74aeae")
            @test !Preferences.has_preference(uuid, "does_not_exist")
            @test Preferences.load_preference(uuid, "toy") == "car"
            @test Preferences.load_preference(uuid, "tree") == "birch"
            @test Preferences.load_preference(uuid, "default") === nothing
        end
    end end
end

@testset "Basic `build` sandbox" begin
    temp_pkg_dir() do project_path; mktempdir() do tmp
        copy_test_package(tmp, "BasicSandbox")
        Pkg.activate(joinpath(tmp, "BasicSandbox"))
        Pkg.build()
    end end
end

end # module


text:
import IterativeSolvers: svdl
import Base: svd, dot, copy

export copy

# fix output of svds to make it like svd
function mysvds(args...; kwargs...)
  svdobj,_ = Base.svds(args...; kwargs...)
  return svdobj.U, svdobj.S, svdobj.Vt
end

function svd(A::AbstractArray, k::Int)
	if k==0 || k>=min(size(A)...)
		return svd(A)
	else
		u,s,v = mysvds(A, nsv = k)
		return u, s, v
		# inaccurate when k=1
		# usv, pf = svdl(A, k, vecs=:both)
		# return usv.U, usv.S, usv.Vt'
	end
end

# function svd(A::AbstractArray, k::Int)
# 	u,s,v = svd(A)
# 	return u[:,1:k], s[1:k], v[:,1:k]
# end

dot(x::Array{Float64,2}, y::Array{Float64,2}) = sum(x.*y)

# find a zero of f over [a, b]
function find_zero(f, a, b; tol=1e-9, maxiters=1000)
    f(a) < 0 || return a
    f(b) > 0 || return b
    for i=1:maxiters
        mid = a + (b-a)/2
        fmid = f(mid)
        if abs(fmid) < tol
            return mid
        end
        if f(mid) < 0
            a = mid
        else
            b = mid
        end
    end
    warn("hit maximum iterations in bisection search")
    return (b-a)/2
end

function copy(p::OptParams)
  newp = typeof(p)()
  for field in fieldnames(p)
    setfield!(newp, field, copy(getfield(p, field)))
  end
  newp
end

function copy(p::StepSizeRule)
  newp = typeof(p)()
  for field in fieldnames(p)
    setfield!(newp, field, copy(getfield(p, field)))
  end
  newp
end

max_stars_repo_name: madeleineudell/FirstOrderOptimization


text:
using JuMP
using SCIP
using Clp;
using Cbc
model = Model(SCIP.Optimizer)
# set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ON)
# set_optimizer_attribute(model, "maxNodes", 5)
month_num = 8
product_num = 1
month = 1:month_num
product = 1:product_num
# paramerters
q = [3000] # product fix cost
p = [60] # cost per product
d = [200, 200, 200, 200, 200, 200, 500, 500]' # demand of product by month
s_init = [100] # initial stock
h = [3] # stock cost per product month

@variable(model, x[i in product, j in month] >= 0)
@variable(model, s[i in product, j in month] >= 0)
@variable(model, y[i in product, j in month], Bin)

# 库存平衡
@constraint(model ,dem_sat1[i in product], s_init[i] + x[i, 1] == d[i, 1] + s[i, 1])
@constraint(model, dem_sat[i in product, j in month[2:end]], s[i, j - 1] + x[i, j] == d[i, j] + s[i, j] )
# 最大产量约束
@constraint(model, vub[i in product, j in month], x[i, j] <= sum(d[i, k]  for k in j:month_num) * y[i, j])
# 目标函数
@expression(model ,cost, sum(p[i] * x[i, j] + q[i] * y[i, j] for i in product, j in month) + sum(h[i] * s[i, j] for i in product, j in month[1:end - 1]) + sum(h[i] * s[i, month[end]] / 2 for i in product))
@objective(model, Min, cost)
optimize!(model)
# 显示结果
cost_val = value.(cost)
x_val = value.(x)
s_val = value.(s)
@show cost_val;
@show x_val;
@show s_val;

model2 = Model(SCIP.Optimizer);
x_val[1,1] = 600
x_val[1,4] = 500
order = 1:8;
orderquantity = [100, 200, 200, 200, 200, 200, 500, 500];
ordertime = [1,2,3,4,5,6,7,8];
@variable(model2, wx[t in month, p in order] >= 0);
@variable(model2, 0 <= wy[t in month, p in order] <= 1, Bin);
@constraint(model2, dem_sato[p in order], sum(wx[t, p] for t in 1:ordertime[p])  == orderquantity[p]);
# @constraint(model2, xlimit[t in month, p in order], wx[t, p] <= x_val[1,t])
@constraint(model2, xlimit[t in month, p in order], wx[t, p] <= x_val[1,t] * wy[t, p]);
@constraint(model2, xeq[t in month], sum(wx[t, p] for p in order) == x_val[1, t]);
@objective(model2, Min, sum(wy[t, p] * ordertime[p] for t in month, p in order));
optimize!(model2);
@show value.(wx);
@show value.(wy);

max_stars_repo_path: example/julia/exercise1.1-1 copy.jl

max_stars_repo_name: Behemoth-s/ProductPlanningByMIP

max_issues_repo_path: example/julia/exercise1.1-1 copy.jl

max_issues_repo_name: Behemoth-s/ProductPlanningByMIP

max_forks_repo_path: example/julia/exercise1.1-1 copy.jl

max_forks_repo_name: Behemoth-s/ProductPlanningByMIP


text:
# Test the creation of SACtr's
using Test
using SAC
@static if VERSION >= v"0.7-"
    import Statistics: mean
end

@testset "Construction" begin
    # Construction from array
    let npts = 1000, delta = 1, b = -5
        local t = rand(SAC.SACFloat, npts)
        local s = SACtr(t, delta)
        @test s.t == t
        @test s.delta == delta
        @test s.npts == npts

        s = SACtr(t, delta, b)
        @test s.b == b
        @test s.e ≈ b + (npts - 1)*delta
        @test s.depmax == maximum(t)
        @test s.depmin == minimum(t)
        @test s.depmen == mean(t)
    end

    # Construction without array
    let npts = rand(1:1000), b = 2, delta = 0.1
        @test SACtr(delta, npts).npts == npts
        @test SACtr(delta, npts).b ≈ 0.0
        @test SACtr(delta, npts, b).b ≈ b
    end

    ## Field modification and access
    # Single traces
    let a = rand(SAC.SACFloat, 100), delta = rand(), s = SACtr(a, delta)
        s[:kevnm] = "ABCD"
        @test s.kevnm == "ABCD"
        @test s[:kevnm] == "ABCD"
        s[:user0] = 1
        @test s.user0 ≈ 1
        @test s[:user0] ≈ 1
        s[:t] = 1:100
        @test s.t ≈ 1:100
        @test s.t ≈ s[:t]
        s[:t] = 2:101
        @test s.t ≈ 2:101
        @test s[:depmin] ≈ 2
        @test s[:depmax] ≈ 101
        @test s[:depmen] ≈ mean(2:101)
        local k = 1
        s[:t] = k
        @test s[:depmin] ≈ k
        @test s[:depmax] ≈ k
        @test all(s[:t] .== SAC.SACFloat(k))
        k = 1:100
        s[:t] = k
        @test all(s[:t] .== SAC.SACFloat.(k))
    end

    # Arrays of traces
    let N = 5, b = -10, delta = 2, npts = 100
        local a = [SACtr(delta, npts, b) for _ in 1:N]
        @test all(a[:b] .≈ b)
        @test all(a[:delta] .≈ delta)
        @test all(a[:npts] .== npts)
        a[:a] = 1:N
        @test all(a[:a] .== SAC.SACFloat.(1:N))
        a[:t] = [1:npts for _ in 1:N]
        @test all(a[:t] .== [SAC.SACFloat.(1:npts) for _ in 1:N])
        local k = -1
        a[:t] = k
        @test all(a[:t] .== [[SAC.SACFloat(k) for j in 1:npts] for i in 1:N])
        a[:kcmpnm] = "ABC"
        @test all(a[:kcmpnm] .== "ABC")
        a[:kcmpnm] = ["a", "b", "c", "d", "e"]
        @test all(a[:kcmpnm] .== ["a", "b", "c", "d", "e"])
    end
end


text:
@generated Base.convert(::Type{T}, x::InfExtendedReal{S}) where {T<:Real,S<:Real} = :(convert($(typeof(convert(T,zero(S)))), x.val))
Base.convert(::Type{Infinite}, x::InfExtendedReal{T}) where {T<:Real} = isinf(x) ? Infinite(signbit(x)) : throw(InexactError(:convert,Infinite,x))

(::Type{T})(x::InfExtendedReal) where {T<:AbstractFloat} = convert(T, x)

Base.widen(::Type{InfExtendedReal{T}}) where {T<:Real} = InfExtendedReal(widen(T))
Base.big(::Type{InfExtendedReal{T}}) where {T<:Real} = InfExtendedReal(big(T))
Base.big(x::InfExtendedReal) = convert(big(typeof(x)), x)

Base.float(::Type{InfExtendedReal{T}}) where {T} = float(T)
Base.float(::Type{InfExtendedReal}) = float(Infinite)


text:
using Base.Test
using QuantumOptics

mutable struct SpectralanalysisTestOperator <: Operator
end

@testset "spectralanalysis" begin

srand(0)

sprandop(b) = sparse(DenseOperator(b, rand(Complex128, length(b), length(b))))

# Test diagonalization
@test_throws ArgumentError eigenstates(SpectralanalysisTestOperator())
@test_throws bases.IncompatibleBases eigenstates(DenseOperator(GenericBasis(3), GenericBasis(4)))
@test_throws ArgumentError eigenenergies(SpectralanalysisTestOperator())
@test_throws bases.IncompatibleBases eigenenergies(DenseOperator(GenericBasis(3), GenericBasis(4)))

# Test hermitian diagonalization
b = GenericBasis(5)
a = randoperator(b)
H = (a+dagger(a))/2
U = expm(1im*H)
d = [-3, -2.6, -0.1, 0.0, 0.6]
D = DenseOperator(b, diagm(d))
Dsp = sparse(D)
@test eigenenergies(D) ≈ d
@test eigenenergies(Dsp, 3) ≈ d[1:3]

R = U*D*dagger(U)
Rsp = sparse(R)
@test eigenenergies((R+dagger(R))/2) ≈ d
@test eigenenergies((Rsp+dagger(Rsp))/2, 3) ≈ d[1:3]

E, states = eigenstates((R+dagger(R))/2, 3)
Esp, states_sp = eigenstates((Rsp+dagger(Rsp))/2, 3)
for i=1:3
    @test E[i] ≈ d[i]
    @test Esp[i] ≈ d[i]
    v = U.data[1,i]/states[i].data[1]
    @test states[i].data*v ≈ U.data[:,i]
    v = U.data[1,i]/states_sp[i].data[1]
    @test states_sp[i].data*v ≈ U.data[:,i]
end

# Test nonhermitian diagonalization
b = GenericBasis(5)
a = randoperator(b)
H = (a+dagger(a))/2
U = expm(1im*H)
d = [-3+0.2im, -2.6-0.1im, -0.1+0.5im, 0.0, 0.6+0.3im]
D = DenseOperator(b, diagm(d))
Dsp = sparse(D)
@test eigenenergies(D; warning=false) ≈ d
@test eigenenergies(Dsp, 3; warning=false) ≈ d[1:3]

R = U*D*dagger(U)
Rsp = sparse(R)
@test eigenenergies(R; warning=false) ≈ d
@test eigenenergies(Rsp, 3; warning=false) ≈ d[1:3]

E, states = eigenstates(R, 3; warning=false)
Esp, states_sp = eigenstates(Rsp, 3; warning=false)
for i=1:3
    @test E[i] ≈ d[i]
    @test Esp[i] ≈ d[i]
    v = U.data[1,i]/states[i].data[1]
    @test states[i].data*v ≈ U.data[:,i]
    v = U.data[1,i]/states_sp[i].data[1]
    @test states_sp[i].data*v ≈ U.data[:,i]
end

# Test simdiag
spinbasis = SpinBasis(1//2)
sx = sigmax(spinbasis)
sy = sigmay(spinbasis)
sz = sigmaz(spinbasis)
twospinbasis = spinbasis ⊗ spinbasis
Sx = full(sum([embed(twospinbasis, i, sx) for i=1:2]))/2.
Sy = full(sum([embed(twospinbasis, i, sy) for i=1:2]))/2.
Sz = full(sum([embed(twospinbasis, i, sz) for i=1:2]))/2.
Ssq = Sx^2 + Sy^2 + Sz^2
d, v = simdiag([Sz, Ssq])
@test d[1] == [-1.0, 0, 0, 1.0]
@test d[2] ≈ [2, 0.0, 2, 2]
@test_throws ErrorException simdiag([Sx, Sy])

threespinbasis = spinbasis ⊗ spinbasis ⊗ spinbasis
Sx3 = full(sum([embed(threespinbasis, i, sx) for i=1:3])/2.)
Sy3 = full(sum([embed(threespinbasis, i, sy) for i=1:3])/2.)
Sz3 = full(sum([embed(threespinbasis, i, sz) for i=1:3])/2.)
Ssq3 = Sx3^2 + Sy3^2 + Sz3^2
d3, v3 = simdiag([Ssq3, Sz3])
dsq3_std = eigenenergies(full(Ssq3))
@test diagm(dsq3_std) ≈ v3'*Ssq3.data*v3

fockbasis = FockBasis(4)
@test_throws ErrorException simdiag([Sy3, Sz3])
@test_throws ErrorException simdiag([full(destroy(fockbasis)), full(create(fockbasis))])

end # testset


text: println("Hello, testing MySDPSolver")

max_stars_repo_name: TomaszWasTaken/MySDPSolver

max_issues_repo_name: TomaszWasTaken/MySDPSolver

max_forks_repo_name: TomaszWasTaken/MySDPSolver


text:
struct Scale{T, N} <: Bijector{N}
    a::T
end

Base.:(==)(b1::Scale{<:Any, N}, b2::Scale{<:Any, N}) where {N} = b1.a == b2.a

function Scale(a::Union{Real,AbstractArray}; dim::Val{D} = Val(ndims(a))) where D
    return Scale{typeof(a), D}(a)
end

# field is a numerical parameter
function Functors.functor(::Type{<:Scale{<:Any,N}}, x) where N
    function reconstruct_scale(xs)
        return Scale{typeof(xs.a),N}(xs.a)
    end
    return (a = x.a,), reconstruct_scale
end

up1(b::Scale{T, N}) where {N, T} = Scale{T, N + 1}(a)

(b::Scale)(x) = b.a .* x
(b::Scale{<:AbstractMatrix, 1})(x::AbstractVecOrMat) = b.a * x
(b::Scale{<:AbstractMatrix, 2})(x::AbstractMatrix) = b.a * x
(ib::Inverse{<:Scale})(y) = Scale(inv(ib.orig.a))(y)
(ib::Inverse{<:Scale{<:AbstractVector}})(y) = Scale(inv.(ib.orig.a))(y)
(ib::Inverse{<:Scale{<:AbstractMatrix, 1}})(y::AbstractVecOrMat) = ib.orig.a \ y
(ib::Inverse{<:Scale{<:AbstractMatrix, 2}})(y::AbstractMatrix) = ib.orig.a \ y

# We're going to implement custom adjoint for this
logabsdetjac(b::Scale{<:Any, N}, x) where {N} = _logabsdetjac_scale(b.a, x, Val(N))

_logabsdetjac_scale(a::Real, x::Real, ::Val{0}) = log(abs(a))
_logabsdetjac_scale(a::Real, x::AbstractVector, ::Val{0}) = fill(log(abs(a)), length(x))
_logabsdetjac_scale(a::Real, x::AbstractVector, ::Val{1}) = log(abs(a)) * length(x)
_logabsdetjac_scale(a::Real, x::AbstractMatrix, ::Val{1}) = fill(log(abs(a)) * size(x, 1), size(x, 2))
_logabsdetjac_scale(a::Real, x::AbstractMatrix, ::Val{2}) = log(abs(a)) * length(x)
_logabsdetjac_scale(a::Real, x::AbstractArray{<:AbstractMatrix}, ::Val{2}) = map(x) do x
    _logabsdetjac_scale(a, x, Val(2))
end
_logabsdetjac_scale(a::AbstractVector, x::AbstractVector, ::Val{1}) = sum(x -> log(abs(x)), a)
_logabsdetjac_scale(a::AbstractVector, x::AbstractMatrix, ::Val{1}) = fill(sum(x -> log(abs(x)), a), size(x, 2))
_logabsdetjac_scale(a::AbstractVector, x::AbstractMatrix, ::Val{2}) = sum(x -> log(abs(x)), a)
_logabsdetjac_scale(a::AbstractVector, x::AbstractArray{<:AbstractMatrix}, ::Val{2}) = map(x) do x
    _logabsdetjac_scale(a, x, Val(2))
end
_logabsdetjac_scale(a::AbstractMatrix, x::AbstractVector, ::Val{1}) = logabsdet(a)[1]
_logabsdetjac_scale(a::AbstractMatrix, x::AbstractMatrix{T}, ::Val{1}) where {T} = logabsdet(a)[1] * ones(T, size(x, 2))
_logabsdetjac_scale(a::AbstractMatrix, x::AbstractMatrix, ::Val{2}) = logabsdet(a)[1]
function _logabsdetjac_scale(
    a::AbstractMatrix,
    x::AbstractArray{<:AbstractMatrix},
    ::Val{2},
)
    map(x) do x
        _logabsdetjac_scale(a, x, Val(2))
    end
end


text:
__precompile__(true)

module DAWG

using LegacyStrings
using SparseArrays

include("dawg.jl")

export dictionaryMatch, SparseMatrixDAWG, DoubleArrayDAWG,
    makeSparseMatrixDAWG, makeDoubleArrayDAWG

end # module