text:
include("common.jl")

example = """
7,4,9,5,11,17,23,2,0,14,21,24,10,16,13,6,15,25,12,22,18,20,8,19,3,26,1

22 13 17 11  0
 8  2 23  4 24
21  9 14 16  7
 6 10  3 18  5
 1 12 20 15 19

 3 15  0  2 22
 9 18 13 17  5
19  8  7 25 23
20 11 10 24  4
14 21 16 12  6

14 21 17 24  4
10 16 15  9 19
18  8 23 26 20
22 11 13  6  5
 2  0 12  3  7
"""

struct Bingo
    order::Vector{Int64}
    size::Tuple{Int64,Int64}
    boards::Vector{Matrix{Int64}}
end

function _gobble!(x::AbstractVector{<:AbstractString})
    while true
        isempty(x) && break
        isempty(@inbounds(first(x))) || break
        popfirst!(x)
    end
    return x
end

function Base.parse(::Type{Bingo}, input::AbstractString)
    lines = split(input, "\n")
    line = popfirst!(lines)
    order = map(x -> parse(Int, x), split(line, ","))

    _gobble!(lines)
    boards = Matrix{Int64}[]
    while !isempty(lines)
        rows = Vector{Int64}[]

        while true
            isempty(lines) && break
            line = popfirst!(lines)
            isempty(line) && break  # End of this board
            push!(rows, parse.(Int, split(line)))
        end

        # Assemble board
        push!(boards, vcat(transpose.(rows)...))

        _gobble!(lines)
    end

    size_ = only(unique(map(size, boards)))
    return Bingo(order, size_, boards)
end

"""
Get the winning (losing) score if `winning` is true (false).
"""
function get_score(bingo::Bingo, winning::Bool)
    markers = [zeros(bingo.size) for _ in 1:length(bingo.boards)]
    row_markers = [zeros(bingo.size[1]) for _ in 1:length(bingo.boards)]
    col_markers = [zeros(bingo.size[2]) for _ in 1:length(bingo.boards)]

    completed_boards = Set{Int64}()

    for number in bingo.order
        @inbounds for i_board in 1:length(bingo.boards)
            # If this board is done, skip it :)
            i_board in completed_boards && continue

            board = bingo.boards[i_board]
            marker = markers[i_board]
            row_marker = row_markers[i_board]
            col_marker = col_markers[i_board]

            for index in findall(==(number), board)
                marker[index] = 1
                i_row, i_col = Tuple(index)
                row_marker[i_row] += 1
                col_marker[i_col] += 1
            end

            # Figure out if board is complete. If so, compute the score.
            complete = (
                any(==(bingo.size[1]), row_marker) ||
                any(==(bingo.size[2]), col_marker)
            )
            complete || continue

            # We have finished!
            push!(completed_boards, i_board)

            # If we're looking for the losing board, only return if all are complete.

            !winning && (length(completed_boards) < length(bingo.boards)) && continue

            # Compute and return the score.
            unmarked_sum = sum(board .* (marker .< 1))
            return unmarked_sum * number
        end
    end
end

get_winning_score(bingo::Bingo) = get_score(bingo, true)
get_losing_score(bingo::Bingo) = get_score(bingo, false)

bingo = parse(Bingo, example)
@assert get_winning_score(bingo) == 4512

get_winning_score(parse(Bingo, load_input(4)))

# Part 2
@assert get_losing_score(parse(Bingo, example)) == 1924
get_losing_score(parse(Bingo, load_input(4)))


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("DataScienceTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

# ## Preliminary steps## Let's start as with the previous tutorial:
using MLJ
using StableRNGs
import DataFrames

@load RidgeRegressor pkg=MultivariateStats

rng = StableRNG(6616) # for reproducibility
x1 = rand(rng, 300)
x2 = rand(rng, 300)
x3 = rand(rng, 300)
y = exp.(x1 - x2 -2x3 + 0.1*rand(rng, 300))
X = DataFrames.DataFrame(x1=x1, x2=x2, x3=x3)

test, train = partition(eachindex(y), 0.8);

# In this tutorial we will show how to generate a model from a network; there are two approaches:# * using the `@from_network` macro# * writing the model in full## the first approach should usually be the one considered as it's simpler.## Generating a model from a network allows subsequent composition of that network with other tasks and tuning of that network.## ## Using the `@from_network` macro## Let's define a simple network## _Input layer_
Xs = source(X)
ys = source(y)

# _First layer_
std_model = Standardizer()
stand = machine(std_model, Xs)
W = MLJ.transform(stand, Xs)

box_model = UnivariateBoxCoxTransformer()
box_mach = machine(box_model, ys)
z = MLJ.transform(box_mach, ys)

# _Second layer_
ridge_model = RidgeRegressor(lambda=0.1)
ridge = machine(ridge_model, W, z)
ẑ = predict(ridge, W)

# _Output_
ŷ = inverse_transform(box_mach, ẑ)

# No fitting has been done thus far, we have just defined a sequence of operations.## As we show next, a learning network needs to be exported to create a new stand-alone model type. Instances of that type can be bound with data in a machine, which can then be evaluated, for example. Somewhat paradoxically, one can wrap a learning network in a certain kind of machine, called a learning network machine, before exporting it, and in fact, the export process actually requires us to do so. Since a composite model type does not yet exist, one constructs the machine using a "surrogate" model, whose name indicates the ultimate model supertype (Deterministic, Probabilistic, Unsupervised or Static). This surrogate model has no fields.
surrogate = Deterministic()
mach = machine(surrogate, Xs, ys; predict=ŷ)

fit!(mach)
predict(mach, X[test[1:5], :])

# To form a model out of that network is easy using the `@from_network` macro.## Having defined a learning network machine, mach, as above, the following code defines a new model subtype WrappedRegressor <: Supervised with a single field regressor
@from_network mach begin
    mutable struct CompositeModel
        regressor=ridge_model
    end
end

# The macro defines a constructor CompositeModel and attributes a name to the# different models; the ordering / connection between the nodes is inferred# from `ŷ` via the `<= ŷ`.## **Note**: had the model been probabilistic (e.g. `RidgeClassifier`) you would have needed to add `is_probabilistic=true` at the end.
cm = machine(CompositeModel(), X, y)
res = evaluate!(cm, resampling=Holdout(fraction_train=0.8, rng=51),
                measure=rms)
round(res.measurement[1], sigdigits=3)

# ## Defining a model from scratch## An alternative to the `@from_network`, is to fully define a new model with its `fit` method:
mutable struct CompositeModel2 <: DeterministicNetwork
    std_model::Standardizer
    box_model::UnivariateBoxCoxTransformer
    ridge_model::RidgeRegressor
end

function MLJ.fit(m::CompositeModel2, verbosity::Int, X, y)
    Xs = source(X)
    ys = source(y)
    W = MLJ.transform(machine(m.std_model, Xs), Xs)
    box = machine(m.box_model, ys)
    z = MLJ.transform(box, ys)
    ẑ = predict(machine(m.ridge_model, W, z), W)
    ŷ = inverse_transform(box, ẑ)
    mach = machine(Deterministic(), Xs, ys; predict=ŷ)
    fit!(mach, verbosity=verbosity - 1)
    return mach()
end

mdl = CompositeModel2(Standardizer(), UnivariateBoxCoxTransformer(),
                      RidgeRegressor(lambda=0.1))
cm = machine(mdl, X, y)
res = evaluate!(cm, resampling=Holdout(fraction_train=0.8), measure=rms)
round(res.measurement[1], sigdigits=3)

# Either way you now have a constructor to a  model which can be used as a stand-alone object, tuned and composed as you would with any basic model.
# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

ext: jl


text:
"compute reaction source term `dC/dt`"
function wdot_func(reaction, T, C, S0, h_mole; get_qdot=false)

    @inbounds _kf = @. @view(reaction.Arrhenius_coeffs[:, 1]) * exp(
        @view(reaction.Arrhenius_coeffs[:, 2]) * log(T) -
        @view(reaction.Arrhenius_coeffs[:, 3]) * (4184.0 / R / T),
    )

    for i in reaction.index_three_body
        @inbounds _kf[i] *= dot(@view(reaction.efficiencies_coeffs[:, i]), C)
    end

    for (j, i) in enumerate(reaction.index_falloff)
        @inbounds A0, b0, Ea0 = reaction.Arrhenius_0[j, :]
        @inbounds k0 = A0 * exp(b0 * log(T) - Ea0 * 4184.0 / R / T)
        @inbounds Pr =
            k0 * dot(@view(reaction.efficiencies_coeffs[:, i]), C) / _kf[i]
        lPr = log10(Pr)
        _kf[i] *= (Pr / (1 + Pr))

        # reference:
        # http://web.mit.edu/2.62/cantera/doc/html/classCantera_1_1Troe4.html#a38aa787421d426dfd0a587fd6fc8108e
        if reaction.index_falloff_Troe[j] > 0
            k = reaction.index_falloff_Troe[j]
            @inbounds F_cent =
                (1 - reaction.Troe_[k, 1]) * exp(-T / reaction.Troe_[k, 4]) +
                reaction.Troe_[k, 1] * exp(-T / reaction.Troe_[k, 2]) +
                exp(-reaction.Troe_[k, 3] / T)

            lF_cent = log10(F_cent)
            _C = -0.4 - 0.67 * lF_cent
            N = 0.75 - 1.27 * lF_cent
            @inbounds f1 = (lPr + _C) / (N - 0.14 * (lPr + _C))
            @inbounds _kf[i] *= exp(log(10.0) * lF_cent / (1 + f1^2))
        end
    end

    @inbounds ΔS_R = reaction.vk' * S0 / R
    @inbounds ΔH_RT = reaction.vk' * h_mole / (R * T)
    @inbounds Keq =
        @. exp(ΔS_R - ΔH_RT + log(one_atm / R / T) * reaction.vk_sum)
    @inbounds _kr = @. _kf / Keq * reaction.is_reversible

    for i = 1:reaction.n_reactions
        @inbounds for j in reaction.i_reactant[i]
            @inbounds _kf[i] *= C[j]^reaction.reactant_orders[j, i]
        end
        if reaction.is_reversible[i]
            @inbounds for j in reaction.i_product[i]
                @inbounds _kr[i] *= C[j]^reaction.product_stoich_coeffs[j, i]
            end
        end
    end

    if get_qdot
        return _kf - _kr
    end
    return reaction.vk * (_kf - _kr)
end
export wdot_func


text:
# algorithm/z-test.jl

using Test: @testset

list = [
"general/pogm_restart.jl"
"general/ncg.jl"
"general/ogm_ls.jl"
"general/poweriter.jl"
]

for file in list
	@testset "$file" begin
		include(file)
	end
end


text:
#
# Basic definitions
#

function check_data_type(data)
  if typeof(data) <: Array && length(data) > 0
    if keytype(data[1]) == String && valtype(data[1]) <: Any
      return true
    end
  end
  return false
end

function stan(
  model::Stanmodel,
  data=Void,
  ProjDir=pwd();
  summary=true,
  diagnostics=false,
  CmdStanDir=CMDSTAN_HOME)

  old = pwd()

  println()
  if length(model.model) == 0
    println("\nNo proper model specified in \"$(model.name).stan\".")
    println("This file is typically created from a String passed to Stanmodel().\n")
    return
  end
  try
    cd(string(Pkg.dir("$(ProjDir)")))
    isfile("$(model.name)_build.log") && rm("$(model.name)_build.log")
    isfile("$(model.name)_run.log") && rm("$(model.name)_run.log")

    #println("Current working dir: $(pwd())")
    #println("Moving to dir: $(CmdStanDir)")

    cd(CmdStanDir)
    local tmpmodelname::String
    tmpmodelname = Pkg.dir(model.tmpdir, model.name)
    if @static is_windows() ? true : false
      tmpmodelname = replace(tmpmodelname*".exe", "\\", "/")
    end
    #println("Current working dir: $(pwd())")
    #println("Compile using make $(tmpmodelname)")
    run(pipeline(`make $(tmpmodelname)`, stderr="$(tmpmodelname)_build.log"))

    cd(model.tmpdir)
    if data != Void && check_data_type(data)
      if length(data) == model.nchains
        for i in 1:model.nchains
          if length(keys(data[i])) > 0
            update_R_file("$(model.name)_$(i).data.R", data[i])
          end
        end
      else
        for i in 1:model.nchains
          if length(keys(data[1])) > 0
            if i == 1
              println("\nLength of data array is not equal to nchains,")
              println("all chains will use the first data dictionary.")
            end
            update_R_file("$(model.name)_$(i).data.R", data[1])
          end
        end
      end
    end
    for i in 1:model.nchains
      model.id = i
      model.data_file ="$(model.name)_$(i).data.R"
      if isa(model.method, Sample)
        model.output.file = model.name*"_samples_$(i).csv"
        isfile("$(model.name)_samples_$(i).csv") && rm("$(model.name)_samples_$(i).csv")
        if diagnostics
          model.output.diagnostic_file = model.name*"_diagnostics_$(i).csv"
          isfile("$(model.name)_diagnostics_$(i).csv") && rm("$(model.name)_diagnostics_$(i).csv")
        end
      elseif isa(model.method, Optimize)
        isfile("$(model.name)_optimize_$(i).csv") && rm("$(model.name)_optimize_$(i).csv")
        model.output.file = model.name*"_optimize_$(i).csv"
      elseif isa(model.method, Variational)
        isfile("$(model.name)_variational_$(i).csv") && rm("$(model.name)_variational_$(i).csv")
        model.output.file = model.name*"_variational_$(i).csv"
      elseif isa(model.method, Diagnose)
        isfile("$(model.name)_diagnose_$(i).csv") && rm("$(model.name)_diagnose_$(i).csv")
        model.output.file = model.name*"_diagnose_$(i).csv"
      end
      model.command[i] = cmdline(model)
    end
    println()
    run(pipeline(par(model.command), stdout="$(model.name)_run.log"))
  catch e
    println(e)
    cd(old)
    return
  end

  local samplefiles = String[]
  local res = Dict[]
  local ftype

  if isa(model.method, Sample)
    ftype = diagnostics ? "diagnostics" : "samples"
    for i in 1:model.nchains
      push!(samplefiles, "$(model.name)_$(ftype)_$(i).csv")
    end
    if isa(model.method, Sample) && summary
      stan_summary(par(samplefiles), CmdStanDir=CmdStanDir)
    end
    #cd(pwd()*"1"); println(pwd())
    res = read_stanfit_samples(model, diagnostics)
  elseif isa(model.method, Variational)
    ftype = "variational"
    for i in 1:model.nchains
      push!(samplefiles, "$(model.name)_$(ftype)_$(i).csv")
    end
    if isa(model.method, Variational) && summary
      stan_summary(par(samplefiles), CmdStanDir=CmdStanDir)
    end
    #cd(pwd()*"1"); println(pwd())
    res = read_stanfit_variational_samples(model)
  elseif isa(model.method, Optimize)
    res = read_stanfit(model)
  elseif isa(model.method, Diagnose)
    res = read_stanfit(model)
  else
    println("Unknown method.")
  end

  cd(old)
  res
end

function update_R_file{T<:Any}(file::String, dct::Dict{String, T}; replaceNaNs::Bool=true)
  isfile(file) && rm(file)
  strmout = open(file, "w")

  str = ""
  for entry in dct
    str = "\""*entry[1]*"\""*" <- "
    val = entry[2]
    if length(val)==1 && length(size(val))==0
      # Scalar
      str = str*"$(val)\n"
    elseif length(val)==1 && length(size(val))==1
      # Single element vector
      str = str*"$(val[1])\n"
    elseif length(val)>1 && length(size(val))==1
      # Vector
      str = str*"structure(c("
      for i in 1:length(val)
        str = str*"$(val[i])"
        if i < length(val)
          str = str*", "
        end
      end
      str = str*"), .Dim=c($(length(val))))\n"
    elseif length(val)>1 && length(size(val))>1
      # Array
      str = str*"structure(c("
      for i in 1:length(val)
        str = str*"$(val[i])"
        if i < length(val)
          str = str*", "
        end
      end
      dimstr = "c"*string(size(val))
      str = str*"), .Dim=$(dimstr))\n"
    end
    write(strmout, str)
  end
  close(strmout)
end

function stan_summary(file::String; CmdStanDir=CMDSTAN_HOME)
  try
    pstring = Pkg.dir("$(CmdStanDir)", "bin", "stansummary")
    cmd = `$(pstring) $(file)`
    print(open(readstring, cmd, "r"))
  catch e
    println(e)
  end
end

function stan_summary(filecmd::Cmd; CmdStanDir=CMDSTAN_HOME)
  try
    pstring = Pkg.dir("$(CmdStanDir)", "bin", "stansummary")
    cmd = `$(pstring) $(filecmd)`
    print(open(readstring, cmd, "r"))
  catch e
    println(e)
  end
end

function read_stanfit(model::Stanmodel)

  ## Collect the results of a chain in an array ##

  chainarray = Dict[]

  ## Each chain dictionary can contain up to 5 types of results ##

  result_type_files = ["samples", "diagnostics", "optimize", "diagnose", "variational"]
  rtdict = Dict()

  ## tdict contains the arrays of values ##
  tdict = Dict()

  for i in 1:model.nchains
    for res_type in result_type_files
      if isfile("$(model.name)_$(res_type)_$(i).csv")
        ## A result type file for chain i is present ##
        ## Result type diagnose needs special treatment ##
        instream = open("$(model.name)_$(res_type)_$(i).csv")
        if res_type == "diagnose"
          tdict = Dict()
          str = readstring(instream)
          sstr = split(str)
          tdict = merge(tdict, Dict(:stan_major_version => [parse(Int, sstr[4])]))
          tdict = merge(tdict, Dict(:stan_minor_version => [parse(Int, sstr[8])]))
          tdict = merge(tdict, Dict(:stan_patch_version => [parse(Int, sstr[12])]))
          if tdict[:stan_minor_version][1] >= 10
            # Stan minor version >= 10
            sstr_lp = sstr[79]
            sstr_lp = parse(Float64, split(sstr_lp, '=')[2])
            tdict = merge(tdict, Dict(:lp => [sstr_lp]))
            tdict = merge(tdict, Dict(:var_id => [parse(Int, sstr[90])]))
            tdict = merge(tdict, Dict(:value => [parse(Float64, sstr[91])]))
            tdict = merge(tdict, Dict(:model => [parse(Float64, sstr[92])]))
            tdict = merge(tdict, Dict(:finite_dif => [parse(Float64, sstr[93])]))
            tdict = merge(tdict, Dict(:error => [parse(Float64, sstr[94])]))
          else
            instream = open("$(model.name)_$(res_type)_$(i).csv")
            skipchars(instream, isspace, linecomment='#')
            i = 0
            skipchars(instream, isspace, linecomment='#')
            while true
              i += 1
              line = normalize_string(readline(instream), newline2lf=true)
              if i == 1
                tdict = merge(tdict, Dict(:lp => [float(split(line[1:(length(line)-1)], "=")[2])]))
              elseif i == 3
                sa = split(line[1:(length(line)-1)])
                tdict = merge(tdict, Dict(:var_id => [parse(Int, sa[1])], :value => [float(sa[2])]))
                tdict = merge(tdict, Dict(:model => [float(sa[3])], :finite_dif => [float(sa[4])]))
                tdict = merge(tdict, Dict(:error => [float(sa[5])]))
              end
              if eof(instream)
                close(instream)
                break
              end
              skipchars(instream, isspace, linecomment='#')
            end
          end
        else
          #println("Type of result file is $(res_type)")
          tdict = Dict()
          skipchars(instream, isspace, linecomment='#')
          line = normalize_string(readline(instream), newline2lf=true)
          idx = split(line[1:length(line)-1], ",")
          index = [idx[k] for k in 1:length(idx)]
          #res_type == "optimize" && println(index)
          j = 0
          skipchars(instream, isspace, linecomment='#')
          while true
            j += 1
            #skipchars(instream, isspace, linecomment='#')
            line = normalize_string(readline(instream), newline2lf=true)
            flds = Float64[]
            #res_type == "optimize" && println(line)
            if eof(instream) && length(line) < 2
              #println("EOF detected")
              close(instream)
              #return(tdict)
              break
            else
              flds = float(split(line[1:length(line)-1], ","))
              #res_type == "optimize" && println(flds)
              for k in 1:length(index)
                if j ==1
                  tdict = merge(tdict, Dict(index[k] => [flds[k]]))
                else
                  tdict[index[k]] = push!(tdict[index[k]], flds[k])
                end
              end
            end
          end
        end

      end

      ## End of processing result type file ##
      ## If any keys were found, merge it in the rtdict ##

      if length(keys(tdict)) > 0
        #println("Merging $(res_type) with keys $(keys(tdict))")
        rtdict = merge(rtdict, Dict(res_type => tdict))
        tdict = Dict()
      end
    end

    ## If rtdict has keys, push it to the chain array ##

    if length(keys(rtdict)) > 0
      #println("Pushing the rtdict with keys $(keys(rtdict))")
      push!(chainarray, rtdict)
      rtdict = Dict()
    end
  end
  chainarray
end

#### Create a Mamba::Chains for results or warmup samples

function read_stanfit_warmup_samples(m::Stanmodel, diagnostics=false)
  if !m.method.save_warmup
    println("Saving of warmup samples not requested in StanModel. Returning samples.")
  end
  Stan.read_stanfit_samples(m::Stanmodel, diagnostics, true)
end

function read_stanfit_samples(m::Stanmodel, diagnostics=false, warmup_samples=false)

  global a3d, monitors, index, idx, indvec, ftype

  ftype = diagnostics ? "diagnostics" : "samples"

  for i in 1:m.nchains
    if isfile("$(m.name)_$(ftype)_$(i).csv")
      noofsamples = 0
      instream = open("$(m.name)_$(ftype)_$(i).csv")
      skipchars(instream, isspace, linecomment='#')
      line = normalize_string(readline(instream), newline2lf=true)
      idx = split(line[1:length(line)-1], ",")
      index = [idx[k] for k in 1:length(idx)]
      if length(m.monitors) == 0
        indvec = 1:length(index)
      else
        indvec = findin(index, m.monitors)
      end
      if m.method.save_warmup
        noofsamples = m.method.num_samples+m.method.num_warmup
      else
        noofsamples = m.method.num_samples
      end
      if i == 1
        a3d = fill(0.0, noofsamples, length(indvec), m.nchains)
      end
      skipchars(instream, isspace, linecomment='#')
      for j in 1:noofsamples
        skipchars(instream, isspace, linecomment='#')
        line = normalize_string(readline(instream), newline2lf=true)
        if eof(instream) && length(line) < 2
          close(instream)
          break
        else
          flds = float(split(line[1:length(line)-1], ","))
          flds = reshape(flds[indvec], 1, length(indvec))
          a3d[j,:,i] = flds
        end
      end
    end
  end

  if m.method.save_warmup
    if warmup_samples
      sr = getindex(a3d, [1:1:m.method.num_warmup;],
        [1:size(a3d, 2);], [1:size(a3d, 3);])
      Chains(sr, start=1, thin=1, names=idx[indvec],
        chains=[i for i in 1:m.nchains])
    else
      sr = getindex(a3d, [m.method.num_warmup+1:m.method.thin:size(a3d, 1);],
        [1:size(a3d, 2);], [1:size(a3d, 3);])
      Chains(sr, start=m.method.num_warmup+1, thin=m.method.thin, names=idx[indvec],
        chains=[i for i in 1:m.nchains])
    end
  else
    sr = getindex(a3d, [1:m.method.thin:size(a3d, 1);],
      [1:size(a3d, 2);], [1:size(a3d, 3);])
    Chains(sr, start=1, thin=m.thin, names=idx[indvec], chains=[i for i in 1:m.nchains])
  end
end

function read_stanfit_variational_samples(m::Stanmodel)

  local a3d, monitors, index, idx, indvec, ftype
  ftype = "variational"

  for i in 1:m.nchains
    if isfile("$(m.name)_$(ftype)_$(i).csv")
      instream = open("$(m.name)_$(ftype)_$(i).csv")
      skipchars(instream, isspace, linecomment='#')
      line = normalize_string(readline(instream), newline2lf=true)
      idx = split(line[1:length(line)-1], ",")
      index = [idx[k] for k in 1:length(idx)]
      if length(m.monitors) == 0
        indvec = 1:length(index)
      else
        indvec = findin(index, m.monitors)
      end
      if i == 1
        a3d = fill(0.0, m.method.output_samples, length(indvec), m.nchains)
      end
      skipchars(instream, isspace, linecomment='#')
      for j in 1:m.method.output_samples
        skipchars(instream, isspace, linecomment='#')
        line = normalize_string(readline(instream), newline2lf=true)
        if eof(instream) && length(line) < 2
          close(instream)
          break
        else
          flds = float(split(line[1:length(line)-1], ","))
          flds = reshape(flds[indvec], 1, length(indvec))
          a3d[j,:,i] = flds
        end
      end
    end
  end

  sr = getindex(a3d, [1:1:size(a3d, 1);], [1:size(a3d, 2);], [1:size(a3d, 3);])
  Chains(sr, start=1, thin=1, names=idx[indvec], chains=[i for i in 1:m.nchains])

end

"Recursively parse the model to construct command line"

function cmdline(m)
  cmd = ``
  if isa(m, Stanmodel)
    # Handle the model name field for unix and windows
    cmd = @static is_unix() ? `./$(getfield(m, :name))` : `cmd /c $(getfield(m, :name)).exe`

    # Method (sample, optimize, variational and diagnose) specific portion of the model
    cmd = `$cmd $(cmdline(getfield(m, :method)))`

    # Common to all models
    cmd = `$cmd $(cmdline(getfield(m, :random)))`
    cmd = `$cmd init=$(getfield(m, :init).init)`

    # Data file required?
    if length(m.data_file) > 0
      cmd = `$cmd id=$(m.id) data file=$(m.data_file)`
    end

    # Output options
    cmd = `$cmd output`
    if length(getfield(m, :output).file) > 0
      cmd = `$cmd file=$(string(getfield(m, :output).file))`
    end
    if length(getfield(m, :output).diagnostic_file) > 0
      cmd = `$cmd diagnostic_file=$(string(getfield(m, :output).diagnostic_file))`
    end
    cmd = `$cmd refresh=$(string(getfield(m, :output).refresh))`
  else
    # The 'recursive' part
    #println(lowercase(string(typeof(m))))
    if isa(m, Algorithm)
      cmd = `$cmd algorithm=$(split(lowercase(string(typeof(m))), '.')[end])`
    elseif isa(m, Engine)
      cmd = `$cmd engine=$(split(lowercase(string(typeof(m))), '.')[end])`
    elseif isa(m, Diagnostics)
      cmd = `$cmd test=$(split(lowercase(string(typeof(m))), '.')[end])`
    else
      cmd = `$cmd $(split(lowercase(string(typeof(m))), '.')[end])`
    end
    #println(cmd)
    for name in fieldnames(m)
      #println("$(name) = $(getfield(m, name)) ($(typeof(getfield(m, name))))")
      if  isa(getfield(m, name), String) || isa(getfield(m, name), Tuple)
        cmd = `$cmd $(name)=$(getfield(m, name))`
      elseif length(fieldnames(typeof(getfield(m, name)))) == 0
        if isa(getfield(m, name), Bool)
          cmd = `$cmd $(name)=$(getfield(m, name) ? 1 : 0)`
        else
          if name == :metric || isa(getfield(m, name), DataType)
            cmd = `$cmd $(name)=$(split(lowercase(string(typeof(getfield(m, name)))), '.')[end])`
          else
            cmd = `$cmd $(name)=$(getfield(m, name))`
          end
        end
      else
        cmd = `$cmd $(cmdline(getfield(m, name)))`
      end
    end
  end
  cmd
end


text:
export to_latex
function to_latex(self::Node)::String
    return to_latex(self, :(=), :left)
end

function to_latex(self::Node{<:Any}, prev_operator::Symbol, pos::Symbol)::String
    curr_operator = Symbol(data(self))
    result = "$(to_latex(args(self)[1], curr_operator, :left)) $(data(self)) $(to_latex(args(self)[2], curr_operator, :right))"

    if curr_operator == prev_operator
        acc = Base.operator_associativity(curr_operator)
        if acc == :left != pos
            return "(" * result * ")"
        elseif acc == :right != pos
            return "(" * result * ")"
        else  # :none
            return result
        end
    end

    curr_precedence = Base.operator_precedence(curr_operator)
    prev_precedence = Base.operator_precedence(prev_operator)

    if curr_precedence > prev_precedence  # if the current operation would be performed last
        return result
    elseif curr_precedence == prev_precedence
        return result
    elseif curr_precedence < prev_precedence
        return "(" * result * ")"
    end
end

function to_latex(self::Node{<:AAtom}, prev_operator::Symbol, pos::Symbol)::String
    return data(self)
end

function to_latex(self::Node{<:ACall}, prev_operator::Symbol, pos::Symbol)::String
    return "$(to_latex(args(self)[1]))($(to_latex(args(self)[2])))"
end


text:
nlevels = 2

nvar = 20

nfactors = [1, 5]

fassign = ones(Int, nvar, 2)
fassign[1:4, 2] = ones(Int, 4)
fassign[5:8, 2] = 2 * ones(Int, 4)
fassign[9:12, 2] = 3 * ones(Int, 4)
fassign[13:16, 2] = 4 * ones(Int, 4)
fassign[17:20, 2] = 5 * ones(Int, 4)

flags = [2, 2]

varlags = 2 * ones(Int, nvar)

varcoefs = zeros(nvar, 1 + nlevels)
varcoefs[:, 2] = 0.5 * ones(nvar)
varcoefs[1, 2] = 1.0
varcoefs[:, 3] = 0.1 * ones(nvar)
varcoefs[1, 3] = 1.0
varcoefs[5, 3] = 1.0
varcoefs[9, 3] = 1.0
varcoefs[13, 3] = 1.0
varcoefs[17, 3] = 1.0

varlagcoefs = ones(nvar, 2)
varlagcoefs[:, 1] = 0.5 * varlagcoefs[:, 1]
varlagcoefs[:, 2] = 0.25 * varlagcoefs[:, 2]

fcoefs = Any[]
fmat = [0.85 -0.3][:, :]
push!(fcoefs, fmat)
fmat = [0.5 0.05
    0.2 -0.1
    0.3 -0.1
    0.4 -0.1
    0.5 -0.1]
push!(fcoefs, fmat)

fvars = Any[]
fmat = [1.0]
push!(fvars, fmat)
fmat = [1.0, 1.0, 1.0, 1.0, 1.0]
push!(fvars, fmat)

varvars = 0.05 * ones(nvar);

hdfm = HDFM(nlevels = nlevels,
    nvar = nvar,
    nfactors = nfactors,
    fassign = fassign,
    flags = flags,
    varlags = varlags,
    varcoefs = varcoefs,
    varlagcoefs = varlagcoefs,
    fcoefs = fcoefs,
    fvars = fvars,
    varvars = varvars)

ssmodel = convertHDFMtoSS(hdfm)

num_obs = 100
data_y, data_z, data_β = simulateSSModel(num_obs, ssmodel::SSModel)

#flags = [3, 3]

#varlags = [3, 3, 3, 3, 3, 3, 3, 3, 3]

hdfmpriors = HDFMStruct(nlevels = nlevels,
    nfactors = nfactors,
    factorassign = fassign,
    factorlags = flags,
    errorlags = varlags,
    ndraws = 1000,
    burnin = 50)

results = OW2LevelEstimator(data_y, hdfmpriors)

medians = Any[]
quant33 = Any[]
quant66 = Any[]
stds = Any[]

j = 4
for i in 1:size(results.F)[1]
    push!(stds, std(results.F[i, j, :]))
    push!(quant33, quantile(results.F[i, j, :], 0.33))
    push!(quant66, quantile(results.F[i, j, :], 0.66))
    push!(medians, median(results.F[i, j, :]))
end

plot(data_β[:, 1+j])
plot!(results.means.F[:, j])
plot!(medians)
plot!(quant33)
plot!(quant66)
plot!(results.means.F[:, j] - stds)
plot!(results.means.F[:, j] + stds)

vardecomp = vardecomp2level(data_y, results.means.F, reshape(results.means.B, 3, nvar)', fassign)
vardecomp2 = vardecomp2level(data_y, data_β[:, 2:7], varcoefs, fassign)


text:
"""
    KSparseLMO{T}(K::Int, right_hand_side::T)

LMO for the K-sparse polytope:
```
C = B_1(τK) ∩ B_∞(τ)
```
with `τ` the `right_hand_side` parameter.
The LMO results in a vector with the K largest absolute values
of direction, taking values `-τ sign(x_i)`.
"""
struct KSparseLMO{T} <: LinearMinimizationOracle
    K::Int
    right_hand_side::T
end

function compute_extreme_point(lmo::KSparseLMO{T}, direction; kwargs...) where {T}
    K = min(lmo.K, length(direction))
    K_indices = sortperm(direction[1:K], by=abs, rev=true)
    K_values = direction[K_indices]
    for idx in K+1:length(direction)
        new_val = direction[idx]
        # new greater value: shift everything right
        if abs(new_val) > abs(K_values[1])
            K_values[2:end] .= K_values[1:end-1]
            K_indices[2:end] .= K_indices[1:end-1]
            K_indices[1] = idx
            K_values[1] = new_val
            # new value in the interior
        elseif abs(new_val) > abs(K_values[K])
            # NOTE: not out of bound since unreachable with K=1
            j = K - 1
            while abs(new_val) > abs(K_values[j])
                j -= 1
            end
            K_values[j+1:end] .= K_values[j:end-1]
            K_indices[j+1:end] .= K_indices[j:end-1]
            K_values[j+1] = new_val
            K_indices[j+1] = idx
        end
    end
    v = spzeros(T, length(direction))
    for (idx, val) in zip(K_indices, K_values)
        v[idx] = -lmo.right_hand_side * sign(val)
    end
    return v
end

function convert_mathopt(
    lmo::KSparseLMO{T},
    optimizer::OT;
    dimension::Integer,
    kwargs...,
) where {T,OT}
    τ = lmo.right_hand_side
    n = dimension
    K = min(lmo.K, n)
    MOI.empty!(optimizer)
    x = MOI.add_variables(optimizer, n)
    tinf = MOI.add_variable(optimizer)
    MOI.add_constraint(optimizer, MOI.VectorOfVariables([tinf; x]), MOI.NormInfinityCone(n + 1))
    MOI.add_constraint(optimizer, tinf, MOI.LessThan(τ))
    t1 = MOI.add_variable(optimizer)
    MOI.add_constraint(optimizer, MOI.VectorOfVariables([t1; x]), MOI.NormOneCone(n + 1))
    MOI.add_constraint(optimizer, t1, MOI.LessThan(τ * K))
    return MathOptLMO(optimizer)
end

"""
    BirkhoffPolytopeLMO

The Birkhoff polytope encodes doubly stochastic matrices.
Its extreme vertices are all permutation matrices of side-dimension `dimension`.
"""
struct BirkhoffPolytopeLMO <: LinearMinimizationOracle end

function compute_extreme_point(
    ::BirkhoffPolytopeLMO,
    direction::AbstractMatrix{T};
    kwargs...,
) where {T}
    n = size(direction, 1)
    n == size(direction, 2) ||
        DimensionMismatch("direction should be square and matching BirkhoffPolytopeLMO dimension")
    res_mat = Hungarian.munkres(direction)
    m = spzeros(Bool, n, n)
    (rows, cols, vals) = SparseArrays.findnz(res_mat)
    @inbounds for i in eachindex(cols)
        m[rows[i], cols[i]] = vals[i] == 2
    end
    m = convert(SparseArrays.SparseMatrixCSC{Float64,Int64}, m)
    return m
end

function compute_extreme_point(
    lmo::BirkhoffPolytopeLMO,
    direction::AbstractVector{T};
    kwargs...,
) where {T}
    nsq = length(direction)
    n = isqrt(nsq)
    return compute_extreme_point(lmo, reshape(direction, n, n); kwargs...)
end

function convert_mathopt(
    ::BirkhoffPolytopeLMO,
    optimizer::OT;
    dimension::Integer,
    kwargs...,
) where {OT}
    n = dimension
    MOI.empty!(optimizer)
    (x, _) = MOI.add_constrained_variables(optimizer, fill(MOI.Interval(0.0, 1.0), n * n))
    xmat = reshape(x, n, n)
    for idx in 1:n
        # column constraint
        MOI.add_constraint(
            optimizer,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(n), xmat[:, idx]), 0.0),
            MOI.EqualTo(1.0),
        )
        # row constraint
        MOI.add_constraint(
            optimizer,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(n), xmat[idx, :]), 0.0),
            MOI.EqualTo(1.0),
        )
    end
    return MathOptLMO(optimizer)
end

"""
    ScaledBoundLInfNormBall(lower_bounds, upper_bounds)

Polytope similar to a L-inf-ball with shifted bounds or general box constraints.
Lower- and upper-bounds are passed on as abstract vectors, possibly of different types.
For the standard L-inf ball, all lower- and upper-bounds would be -1 and 1.
"""
struct ScaledBoundLInfNormBall{T, VT1 <: AbstractVector{T}, VT2 <: AbstractVector{T}} <: LinearMinimizationOracle
    lower_bounds::VT1
    upper_bounds::VT2
end

function compute_extreme_point(lmo::ScaledBoundLInfNormBall, direction; kwargs...)
    v = copy(lmo.lower_bounds)
    for i in eachindex(direction)
        if direction[i] * lmo.upper_bounds[i] < direction[i] * lmo.lower_bounds[i]
            v[i] = lmo.upper_bounds[i]
        end
    end
    return v
end

"""
    ScaledBoundL1NormBall(lower_bounds, upper_bounds)

Polytope similar to a L1-ball with shifted bounds.
It is the convex hull of two scaled and shifted unit vectors for each axis (shifted to the center of the polytope, i.e., the elementwise midpoint of the bounds).
Lower and upper bounds are passed on as abstract vectors, possibly of different types.
For the standard L1-ball, all lower and upper bounds would be -1 and 1.
"""
struct ScaledBoundL1NormBall{T, VT1 <: AbstractVector{T}, VT2 <: AbstractVector{T}} <: LinearMinimizationOracle
    lower_bounds::VT1
    upper_bounds::VT2
end

function compute_extreme_point(lmo::ScaledBoundL1NormBall, direction; kwargs...)
    idx = 0
    lower = false
    val = zero(eltype(direction))
    if length(direction) != length(lmo.upper_bounds)
        throw(DimensionMismatch())
    end
    @inbounds for i in eachindex(direction)
        scale_factor = lmo.upper_bounds[i] - lmo.lower_bounds[i]
        scaled_dir = direction[i] * scale_factor
        if scaled_dir > val
            val = scaled_dir
            idx = i
            lower = true
        elseif -scaled_dir > val
            val = -scaled_dir
            idx = i
            lower = false
        end
    end
    # compute midpoint for all coordinates, replace with extreme coordinate on one
    # TODO use smarter array type if bounds are FillArrays
    v = (lmo.lower_bounds + lmo.upper_bounds) / 2
    # handle zero direction
    idx = max(idx, 1)
    v[idx] = ifelse(lower, lmo.lower_bounds[idx], lmo.upper_bounds[idx])
    return v
end


text:
# TODO add Kappa and other relevant distributions

"""
Temperature in eV
mass in Kg
"""
function Maxwellian(T, m; v_drift=zeros())
        Kb = 8.617333262145e-5
        v_th = sqrt(2*Kb*T/m)

        function P(x,v)
            if !(v isa Array)
                v = [v]
            end

            if v_drift == zeros()
                v_drift = zeros(length(v))
            end

            if length(v) != length(v_drift)
                error("v and v_drift should have the same length")
            end

            v_ = sqrt(sum(v .^2))
            v_drift_ = sqrt(sum(v_drift.^2))

            (π*v_th^2)^(-3/2) * exp(-(v_ - v_drift_)^2/v_th^2)
        end
end


text:
using ArrayLayouts, FillArrays, Test
import ArrayLayouts: DenseColumnMajor, AbstractStridedLayout, AbstractColumnMajor, DiagonalLayout

@testset "Matrix * Vector" begin
    @testset "eltype" begin
        v = MulAdd(1,zeros(Int,2,2), zeros(Float64,2),0,zeros(2))
        A = MulAdd(1,zeros(Int,2,2), zeros(Float64,2,2), 0,zeros(2,2))
        @test eltype(v) == eltype(A) == Float64
        @test @inferred(axes(v)) == (@inferred(axes(v,1)),) == (Base.OneTo(2),)
        @test @inferred(size(v)) == (@inferred(size(v,1)),) == (2,)
        @test @inferred(axes(A)) == (@inferred(axes(A,1)),@inferred(axes(A,2))) == (Base.OneTo(2),Base.OneTo(2))
        @test @inferred(size(A)) == (@inferred(size(A,1)),size(A,2)) == (2,2)

        Ã = materialize(A)
        @test Ã isa Matrix{Float64}
        @test Ã == zeros(2,2)

        c = similar(v)
        fill!(c,NaN)
        @test copyto!(c,v) == zeros(2)
        fill!(c,NaN)
        c .= v
        @test c == zeros(2)
    end

    @testset "gemv Float64" begin
        for A in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:),
                  view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)),
            b in (randn(5), view(randn(5),:), view(randn(5),1:5), view(randn(9),1:2:9))
            c = similar(b);

            muladd!(1.0,A,b,0.0,c)
            @test all(c .=== A*b .=== BLAS.gemv!('N', 1.0, A, b, 0.0, similar(c)))
            copyto!(c, MulAdd(1.0,A,b,0.0,c))
            @test all(c .=== A*b .=== BLAS.gemv!('N', 1.0, A, b, 0.0, similar(c)))
            c .= MulAdd(1.0,A,b,0.0,c)
            @test all(c .=== A*b .=== BLAS.gemv!('N', 1.0, A, b, 0.0, similar(c)))

            M = MulAdd(2.0, A,b, 0.0, c)
            @test M isa MulAdd{<:AbstractColumnMajor,<:AbstractStridedLayout,<:AbstractColumnMajor}
            c .= M
            @test all(c .=== BLAS.gemv!('N', 2.0, A, b, 0.0, similar(c)))
            copyto!(c, M)
            @test all(c .=== BLAS.gemv!('N', 2.0, A, b, 0.0, similar(c)))

            c = copy(b)
            M = MulAdd(1.0, A,b, 1.0, c)
            @test M isa MulAdd{<:AbstractColumnMajor,<:AbstractStridedLayout,<:AbstractColumnMajor}
            copyto!(c, M)
            @test all(c .=== BLAS.gemv!('N', 1.0, A, b, 1.0, copy(b)))
        end
    end

    @testset "gemv mixed array types" begin
        (A, b, c) = (randn(5,5), randn(5), 1.0:5.0)
        d = similar(b)
        d .= MulAdd(1.0,A,b,1.0,c)
        @test all(d .=== BLAS.gemv!('N', 1.0, A, b, 1.0, Vector{Float64}(c)))
    end

    @testset "gemv Complex" begin
        for T in (ComplexF64,),
                A in (randn(T,5,5), view(randn(T,5,5),:,:)),
                b in (randn(T,5), view(randn(T,5),:))
            c = similar(b);

            c .= MulAdd(one(T), A, b, zero(T), c)
            @test all(c .=== BLAS.gemv!('N', one(T), A, b, zero(T), similar(c)))

            c = copy(b)
            muladd!(2one(T),A,b,one(T),c)
            @test all(c .=== BLAS.gemv!('N', 2one(T), A, b, one(T), copy(b)))
        end
    end
end

@testset "Matrix * Matrix" begin
    @testset "gemm" begin
        for A in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:),
                  view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)),
            B in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:),
                      view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5))
            C = similar(B);

            C .= MulAdd(1.0,A,B,0.0,C)
            @test all(C .=== BLAS.gemm!('N', 'N', 1.0, A, B, 0.0, similar(C)))

            C .= MulAdd(2.0,A,B,0.0,C)
            @test all(C .=== BLAS.gemm!('N', 'N', 2.0, A, B, 0.0, similar(C)))

            C = copy(B)
            C .= MulAdd(2.0,A,B,1.0,C)
            @test all(C .=== BLAS.gemm!('N', 'N', 2.0, A, B, 1.0, copy(B)))
        end
    end

    @testset "gemm Complex" begin
        for T in (ComplexF64,),
                A in (randn(T,5,5), view(randn(T,5,5),:,:)),
                B in (randn(T,5,5), view(randn(T,5,5),:,:))
            C = similar(B);

            C .= MulAdd(one(T),A,B,zero(T),C)
            @test all(C .=== BLAS.gemm!('N', 'N', one(T), A, B, zero(T), similar(C)))

            C .= MulAdd(2one(T),A,B,zero(T),C)
            @test all(C .=== BLAS.gemm!('N', 'N', 2one(T), A, B, zero(T), similar(C)))

            C = copy(B)
            C .= MulAdd(one(T),A,B,one(T),C)
            @test all(C .=== BLAS.gemm!('N', 'N', one(T), A, B, one(T), copy(B)))
        end
    end

    @testset "gemm mixed array types" begin
        (A, B, C) = (randn(5,5), randn(5,5), reshape(1.0:25.0,5,5))
        D = similar(B)
        D .= MulAdd(1.0,A,B,1.0, C)
        @test all(D .=== BLAS.gemm!('N', 'N', 1.0, A, B, 1.0, Matrix{Float64}(C)))
    end
end

@testset "adjtrans" begin
    @testset "gemv adjtrans" begin
        for A in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:),
                  view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)),
                    b in (randn(5), view(randn(5),:), view(randn(5),1:5), view(randn(9),1:2:9)),
                    Ac in (transpose(A), A')
            c = similar(b);

            c .= MulAdd(1.0,Ac,b,0.0,c)
            @test all(c .=== BLAS.gemv!('T', 1.0, A, b, 0.0, similar(c)))

            c .= MulAdd(2.0,Ac,b,0.0,c)
            @test all(c .=== BLAS.gemv!('T', 2.0, A, b, 0.0, similar(c)))

            c = copy(b)
            c .= MulAdd(1.0,Ac,b,1.0,c)
            @test all(c .=== BLAS.gemv!('T', 1.0, A, b, 1.0, copy(b)))
        end
    end

    @testset "gemv adjtrans mixed types" begin
        (A, b, c) = (randn(5,5), randn(5), 1.0:5.0)
        for Ac in (transpose(A), A')
            d = similar(b)
            d .= MulAdd(1.0,Ac,b,1.0, c)
            @test all(d .=== BLAS.gemv!('T', 1.0, A, b, 1.0, Vector{Float64}(c)))

            d .= MulAdd(3.0,Ac,b,2.0, c)
            @test all(d .=== BLAS.gemv!('T', 3.0, A, b, 2.0, Vector{Float64}(c)))
        end
    end

    @testset "gemv adjtrans Complex" begin
        for T in (ComplexF64,),
                A in (randn(T,5,5), view(randn(T,5,5),:,:)),
                b in (randn(T,5), view(randn(T,5),:)),
                (Ac,trans) in ((transpose(A),'T'), (A','C'))
            c = similar(b);

            c .= MulAdd(one(T),Ac,b,zero(T),c)
            @test all(c .=== BLAS.gemv!(trans, one(T), A, b, zero(T), similar(c)))

            c = copy(b)
            c .= MulAdd(3one(T),Ac,b,2one(T),c)
            @test all(c .=== BLAS.gemv!(trans, 3one(T), A, b, 2one(T), copy(b)))
        end
    end

    @testset "gemm adjtrans" begin
        for A in (randn(5,5), view(randn(5,5),:,:)),
            B in (randn(5,5), view(randn(5,5),1:5,:))
            for Ac in (transpose(A), A')
                C = copy(B)
                C .= MulAdd(3.0, Ac, B, 2.0, C)
                @test all(C .=== BLAS.gemm!('T', 'N', 3.0, A, B, 2.0, copy(B)))
            end
            for Bc in (transpose(B), B')
                C = copy(B)
                C .= MulAdd(3.0, A, Bc, 2.0, C)
                @test all(C .=== BLAS.gemm!('N', 'T', 3.0, A, B, 2.0, copy(B)))
            end
            for Ac in (transpose(A), A'), Bc in (transpose(B), B')
                C = copy(B)
                d = similar(C)
                d .= MulAdd(3.0, Ac, Bc, 2.0, C)
                @test all(d .=== BLAS.gemm!('T', 'T', 3.0, A, B, 2.0, copy(B)))
            end
        end
    end

    @testset "symv adjtrans" begin
        A = randn(100,100)
        x = randn(100)
        for As in (Symmetric(A), Hermitian(A), Symmetric(A)', Symmetric(view(A,:,:)',:L), view(Symmetric(A),:,:))
            @test all( muladd!(1.0,As,x,0.0,similar(x)) .=== BLAS.symv!('U', 1.0, A, x, 0.0, similar(x)) )
        end

        for As in (Symmetric(A,:L), Hermitian(A,:L), Symmetric(A,:L)', Symmetric(view(A,:,:)',:U), view(Hermitian(A,:L),:,:))
            @test all( muladd!(1.0,As,x,0.0,similar(x)) .=== BLAS.symv!('L', 1.0, A, x, 0.0, similar(x)) )
        end

        T = ComplexF64
        A = randn(T,100,100)
        x = randn(T,100)
        for As in (Symmetric(A), transpose(Symmetric(A)), Symmetric(transpose(view(A,:,:)),:L), view(Symmetric(A),:,:))
            @test all( muladd!(one(T),As,x,zero(T),similar(x)) .=== BLAS.symv!('U', one(T), A, x, zero(T), similar(x)) )
        end

        for As in (Symmetric(A,:L), transpose(Symmetric(A,:L)), Symmetric(transpose(view(A,:,:)),:U), view(Symmetric(A,:L),:,:))
            @test all( muladd!(one(T),As,x,zero(T),similar(x)) .=== BLAS.symv!('L', one(T), A, x, zero(T), similar(x)) )
        end

        for As in(Hermitian(A), Hermitian(A)', view(Hermitian(A),:,:))
            @test all( muladd!(one(T),As,x,zero(T),similar(x)) .=== BLAS.hemv!('U', one(T), A, x, zero(T), similar(x)) )
        end

        for As in(Hermitian(A,:L), Hermitian(A,:L)', view(Hermitian(A,:L),:,:))
            @test all( muladd!(one(T),As,x,zero(T),similar(x)) .=== BLAS.hemv!('L', one(T), A, x, zero(T), similar(x)) )
        end
    end
end

@testset "Mul" begin
    @testset "Mixed types" begin
        A = randn(5,6)
        b = rand(Int,6)
        c = Array{Float64}(undef, 5)
        d = similar(c)
        copyto!(d, MulAdd(1, A, b, 0.0, d))
        @test d == copyto!(similar(d), MulAdd(1, A, b, 0.0, d)) ≈ A*b
        @test copyto!(similar(d), MulAdd(1, A, b, 1.0, d)) ≈ A*b + d

        @test all((similar(d) .= MulAdd(1, A, b, 1.0, d)) .=== copyto!(similar(d), MulAdd(1, A, b, 1.0, d)))

        B = rand(Int,6,4)
        C = Array{Float64}(undef, 5, 4)

        @test_throws DimensionMismatch materialize!(MulAdd(1.0,A,B,0.0,similar(C,4,4)))

        A = randn(Float64,20,22)
        B = randn(ComplexF64,22,24)
        C = similar(B,20,24)
        @test copyto!(similar(C), MulAdd(1.0, A, B, 0.0, C)) ≈ A*B
    end

    @testset "no allocation" begin
        A = randn(5,5); x = randn(5); y = randn(5); c = similar(y);
        muladd!(2.0, A, x, 3.0, y)
        @test @allocated(muladd!(2.0, A, x, 3.0, y)) == 0
        VA = view(A,:,1:2)
        vx = view(x,1:2)
        vy = view(y,:)
        muladd!(2.0, VA, vx, 3.0, vy)
        @test @allocated(muladd!(2.0, VA, vx, 3.0, vy)) == 0
    end

    @testset "BigFloat" begin
        A = BigFloat.(randn(5,5))
        x = BigFloat.(randn(5))
        @test A*x == muladd!(1.0,A,x,0.0,copy(x))
        @test_throws UndefRefError muladd!(1.0,A,x,0.0,similar(x))
    end

    @testset "Scalar * Vector" begin
        A, x =  [1 2; 3 4] , [[1,2],[3,4]]
        @test muladd!(1.0,A,x,0.0, 1.0*x) == A*x
    end
end

@testset "Lmul/Rmul" begin
    @testset "tri Lmul" begin
        @testset "Float * Float vector" begin
            A = randn(Float64, 100, 100)
            x = randn(Float64, 100)

            L = Lmul(UpperTriangular(A),x)
            @test size(L) == (size(L,1),) == (100,)
            @test axes(L) == (axes(L,1),) == (Base.OneTo(100),)
            @test eltype(L) == Float64
            @test length(L) == 100

            @test similar(L) isa Vector{Float64}
            @test similar(L,Int) isa Vector{Int}

            @test all(copy(Lmul(UpperTriangular(A),x)) .===
                        materialize!(Lmul(UpperTriangular(A),copy(x))) .===
                        copyto!(similar(x),Lmul(UpperTriangular(A),x)) .===
                        UpperTriangular(A)*x .===
                        BLAS.trmv!('U', 'N', 'N', A, copy(x)))
            @test all(copyto!(similar(x),Lmul(UnitUpperTriangular(A),x)) .===
                        UnitUpperTriangular(A)*x .===
                        BLAS.trmv!('U', 'N', 'U', A, copy(x)))
            @test all(copyto!(similar(x),Lmul(LowerTriangular(A),x)) .===
                        LowerTriangular(A)*x .===
                        BLAS.trmv!('L', 'N', 'N', A, copy(x)))
            @test all(materialize!(Lmul(UnitLowerTriangular(A),copy(x))) .===
                        UnitLowerTriangular(A)x .===
                        BLAS.trmv!('L', 'N', 'U', A, copy(x)))

            @test all(materialize!(Lmul(UpperTriangular(A)',copy(x))) .===
                        materialize!(Lmul(LowerTriangular(A'),copy(x))) .===
                        UpperTriangular(A)'x .===
                        BLAS.trmv!('U', 'T', 'N', A, copy(x)))
            @test all(materialize!(Lmul(UnitUpperTriangular(A)',copy(x))) .===
                        materialize!(Lmul(UnitLowerTriangular(A'),copy(x))) .===
                        UnitUpperTriangular(A)'*x .===
                        BLAS.trmv!('U', 'T', 'U', A, copy(x)))
            @test all(materialize!(Lmul(LowerTriangular(A)',copy(x))) .===
                        materialize!(Lmul(UpperTriangular(A'),copy(x))) .===
                        LowerTriangular(A)'*x .===
                        BLAS.trmv!('L', 'T', 'N', A, copy(x)))
            @test all(materialize!(Lmul(UnitLowerTriangular(A)',copy(x))) .===
                        materialize!(Lmul(UnitUpperTriangular(A'),copy(x))) .===
                        UnitLowerTriangular(A)'*x .===
                        BLAS.trmv!('L', 'T', 'U', A, copy(x)))
        end

        @testset "Float * Complex vector"  begin
            T = ComplexF64
            A = randn(T, 100, 100)
            x = randn(T, 100)

            @test all((y = copy(x); y .= Lmul(UpperTriangular(A),y) ) .===
                        UpperTriangular(A)x .===
                        BLAS.trmv!('U', 'N', 'N', A, copy(x)))
            @test all((y = copy(x); y .= Lmul(UnitUpperTriangular(A),y) ) .===
                        UnitUpperTriangular(A)x .===
                        BLAS.trmv!('U', 'N', 'U', A, copy(x)))
            @test all((y = copy(x); y .= Lmul(LowerTriangular(A),y) ) .===
                        LowerTriangular(A)x .===
                        BLAS.trmv!('L', 'N', 'N', A, copy(x)))
            @test all((y = copy(x); y .= Lmul(UnitLowerTriangular(A),y) ) .===
                        UnitLowerTriangular(A)x .===
                        BLAS.trmv!('L', 'N', 'U', A, copy(x)))
            @test LowerTriangular(A')  == UpperTriangular(A)'

            @test all((y = copy(x); y .= Lmul(transpose(UpperTriangular(A)),y) ) .===
                        (similar(x) .= Lmul(LowerTriangular(transpose(A)),x)) .===
                        transpose(UpperTriangular(A))x .===
                        BLAS.trmv!('U', 'T', 'N', A, copy(x)))
            @test all((y = copy(x); y .= Lmul(transpose(UnitUpperTriangular(A)),y) ) .===
                        (similar(x) .= Lmul(UnitLowerTriangular(transpose(A)),x)) .===
                        transpose(UnitUpperTriangular(A))x .===
                        BLAS.trmv!('U', 'T', 'U', A, copy(x)))
            @test all((y = copy(x); y .= Lmul(transpose(LowerTriangular(A)),y) ) .===
                        (similar(x) .= Lmul(UpperTriangular(transpose(A)),x)) .===
                        transpose(LowerTriangular(A))x .===
                        BLAS.trmv!('L', 'T', 'N', A, copy(x)))
            @test all((y = copy(x); y .= Lmul(transpose(UnitLowerTriangular(A)),y) ) .===
                        (similar(x) .= Lmul(UnitUpperTriangular(transpose(A)),x)) .===
                        transpose(UnitLowerTriangular(A))x .===
                        BLAS.trmv!('L', 'T', 'U', A, copy(x)))

            @test all((y = copy(x); y .= Lmul(UpperTriangular(A)',y) ) .===
                        (similar(x) .= Lmul(LowerTriangular(A'),x)) .===
                        UpperTriangular(A)'x .===
                        BLAS.trmv!('U', 'C', 'N', A, copy(x)))
            @test all((y = copy(x); y .= Lmul(UnitUpperTriangular(A)',y) ) .===
                        (similar(x) .= Lmul(UnitLowerTriangular(A'),x)) .===
                        UnitUpperTriangular(A)'x .===
                        BLAS.trmv!('U', 'C', 'U', A, copy(x)))
            @test all((y = copy(x); y .= Lmul(LowerTriangular(A)',y) ) .===
                        (similar(x) .= Lmul(UpperTriangular(A'),x)) .===
                        LowerTriangular(A)'x .===
                        BLAS.trmv!('L', 'C', 'N', A, copy(x)))
            @test all((y = copy(x); y .= Lmul(UnitLowerTriangular(A)',y) ) .===
                        (similar(x) .= Lmul(UnitUpperTriangular(A'),x)) .===
                        UnitLowerTriangular(A)'x .===
                        BLAS.trmv!('L', 'C', 'U', A, copy(x)))
        end

        @testset "Float * Float Matrix" begin
            A = randn(Float64, 100, 100)
            x = randn(Float64, 100, 100)

            @test UpperTriangular(A)*x ≈ (similar(x) .= Lmul(UpperTriangular(A), x))
        end

        @testset "adjtrans" begin
            for T in (Float64, ComplexF64)
                A = randn(T,100,100)
                b = randn(T,100)

                @test all(copy(Lmul(UpperTriangular(A)', b)) .=== UpperTriangular(A)'b)
                @test all(copy(Lmul(UnitUpperTriangular(A)', b)) .=== UnitUpperTriangular(A)'b)
                @test all(copy(Lmul(LowerTriangular(A)', b)) .=== LowerTriangular(A)'b)
                @test all(copy(Lmul(UnitLowerTriangular(A)', b)) .=== UnitLowerTriangular(A)'b)

                @test all(copy(Lmul(transpose(UpperTriangular(A)), b)) .=== transpose(UpperTriangular(A))b)
                @test all(copy(Lmul(transpose(UnitUpperTriangular(A)), b)) .=== transpose(UnitUpperTriangular(A))b)
                @test all(copy(Lmul(transpose(LowerTriangular(A)), b)) .=== transpose(LowerTriangular(A))b)
                @test all(copy(Lmul(transpose(UnitLowerTriangular(A)), b)) .=== transpose(UnitLowerTriangular(A))b)

                B = randn(T,100,100)

                @test all(copy(Lmul(UpperTriangular(A)', B)) .=== UpperTriangular(A)'B)
                @test all(copy(Lmul(UnitUpperTriangular(A)', B)) .=== UnitUpperTriangular(A)'B)
                @test all(copy(Lmul(LowerTriangular(A)', B)) .=== LowerTriangular(A)'B)
                @test all(copy(Lmul(UnitLowerTriangular(A)', B)) .=== UnitLowerTriangular(A)'B)

                @test all(copy(Lmul(transpose(UpperTriangular(A)), B)) .=== transpose(UpperTriangular(A))B)
                @test all(copy(Lmul(transpose(UnitUpperTriangular(A)), B)) .=== transpose(UnitUpperTriangular(A))B)
                @test all(copy(Lmul(transpose(LowerTriangular(A)), B)) .=== transpose(LowerTriangular(A))B)
                @test all(copy(Lmul(transpose(UnitLowerTriangular(A)), B)) .=== transpose(UnitLowerTriangular(A))B)
            end

            for T in (Float64, ComplexF64)
                A = big.(randn(T,100,100))
                b = big.(randn(T,100))

                @test all(copy(Lmul(UpperTriangular(A)', b)) ≈ UpperTriangular(A)'b)
                @test all(copy(Lmul(UnitUpperTriangular(A)', b)) ≈ UnitUpperTriangular(A)'b)
                @test all(copy(Lmul(LowerTriangular(A)', b)) ≈ LowerTriangular(A)'b)
                @test all(copy(Lmul(UnitLowerTriangular(A)', b)) ≈ UnitLowerTriangular(A)'b)

                @test all(copy(Lmul(transpose(UpperTriangular(A)), b)) ≈ transpose(UpperTriangular(A))b)
                @test all(copy(Lmul(transpose(UnitUpperTriangular(A)), b)) ≈ transpose(UnitUpperTriangular(A))b)
                @test all(copy(Lmul(transpose(LowerTriangular(A)), b)) ≈ transpose(LowerTriangular(A))b)
                @test all(copy(Lmul(transpose(UnitLowerTriangular(A)), b)) ≈ transpose(UnitLowerTriangular(A))b)

                B = big.(randn(T,100,100))

                @test all(copy(Lmul(UpperTriangular(A)', B)) ≈ UpperTriangular(A)'B)
                @test all(copy(Lmul(UnitUpperTriangular(A)', B)) ≈ UnitUpperTriangular(A)'B)
                @test all(copy(Lmul(LowerTriangular(A)', B)) ≈ LowerTriangular(A)'B)
                @test all(copy(Lmul(UnitLowerTriangular(A)', B)) ≈ UnitLowerTriangular(A)'B)

                @test all(copy(Lmul(transpose(UpperTriangular(A)), B)) ≈ transpose(UpperTriangular(A))B)
                @test all(copy(Lmul(transpose(UnitUpperTriangular(A)), B)) ≈ transpose(UnitUpperTriangular(A))B)
                @test all(copy(Lmul(transpose(LowerTriangular(A)), B)) ≈ transpose(LowerTriangular(A))B)
                @test all(copy(Lmul(transpose(UnitLowerTriangular(A)), B)) ≈ transpose(UnitLowerTriangular(A))B)
            end
        end
    end

    @testset "tri Rmul" begin
        for T in (Float64, ComplexF64)
            A = randn(T, 100,100)
            B = randn(T, 100,100)
            R = Rmul(copy(A), UpperTriangular(B))
            @test size(R) == (size(R,1),size(R,2)) == (100,100)
            @test axes(R) == (axes(R,1),axes(R,2)) == (Base.OneTo(100),Base.OneTo(100))
            @test eltype(R) == T
            @test length(R) == 100^2

            @test similar(R) isa Matrix{T}
            @test similar(R,Int) isa Matrix{Int}

            R2 = deepcopy(R)
            @test all(BLAS.trmm('R', 'U', 'N', 'N', one(T), B, A) .=== copyto!(similar(R2), R2) .=== materialize!(R))
            @test R.A ≠ A
            @test all(BLAS.trmm('R', 'U', 'T', 'N', one(T), B, A) .=== copy(Rmul(A, transpose(UpperTriangular(B)))) .=== A*transpose(UpperTriangular(B)))
            @test all(BLAS.trmm('R', 'U', 'N', 'U', one(T), B, A) .=== copy(Rmul(A, UnitUpperTriangular(B))).=== A*UnitUpperTriangular(B))
            @test all(BLAS.trmm('R', 'U', 'T', 'U', one(T), B, A) .=== copy(Rmul(A, transpose(UnitUpperTriangular(B)))) .=== A*transpose(UnitUpperTriangular(B)))
            @test all(BLAS.trmm('R', 'L', 'N', 'N', one(T), B, A) .=== copy(Rmul(A, LowerTriangular(B))).=== A*LowerTriangular(B))
            @test all(BLAS.trmm('R', 'L', 'T', 'N', one(T), B, A) .=== copy(Rmul(A, transpose(LowerTriangular(B)))) .=== A*transpose(LowerTriangular(B)))
            @test all(BLAS.trmm('R', 'L', 'N', 'U', one(T), B, A) .=== copy(Rmul(A, UnitLowerTriangular(B))).=== A*UnitLowerTriangular(B))
            @test all(BLAS.trmm('R', 'L', 'T', 'U', one(T), B, A) .=== copy(Rmul(A, transpose(UnitLowerTriangular(B)))) .=== A*transpose(UnitLowerTriangular(B)))

            if T <: Complex
                @test all(BLAS.trmm('R', 'U', 'C', 'N', one(T), B, A) .=== copy(Rmul(A, UpperTriangular(B)')) .=== A*UpperTriangular(B)')
                @test all(BLAS.trmm('R', 'U', 'C', 'U', one(T), B, A) .=== copy(Rmul(A, UnitUpperTriangular(B)')) .=== A*UnitUpperTriangular(B)')
                @test all(BLAS.trmm('R', 'L', 'C', 'N', one(T), B, A) .=== copy(Rmul(A, LowerTriangular(B)')) .=== A*LowerTriangular(B)')
                @test all(BLAS.trmm('R', 'L', 'C', 'U', one(T), B, A) .=== copy(Rmul(A, UnitLowerTriangular(B)')) .=== A*UnitLowerTriangular(B)')
            end
        end

        T = Float64
        A = big.(randn(100,100))
        B = big.(randn(100,100))
        @test materialize!(Rmul(copy(A),UpperTriangular(B))) == A*UpperTriangular(B)
    end

    @testset "Diagonal and SymTridiagonal" begin
        A = randn(5,5)
        B = Diagonal(randn(5))
        @test MemoryLayout(typeof(B)) == DiagonalLayout{DenseColumnMajor}()

        @test A*B == materialize!(Rmul(copy(A),B))
        @test B*A == materialize!(Lmul(B,copy(A)))
    end
end

@testset "MulAdd" begin
    A = randn(5,5)
    B = randn(5,4)
    C = randn(5,4)
    b = randn(5)
    c = randn(5)

    M = MulAdd(2.0,A,B,3.0,C)
    @test size(M) == size(C)
    @test size(M,1) == size(C,1)
    @test size(M,2) == size(C,2)
    @test_broken size(M,3) == size(C,3)
    @test length(M) == length(C)
    @test axes(M) == axes(C)
    @test eltype(M) == Float64
    @test materialize(M) ≈ 2.0A*B + 3.0C

    @test_throws DimensionMismatch materialize(MulAdd(2.0,A,randn(3),1.0,B))
    @test_throws DimensionMismatch materialize(MulAdd([1,2],A,B,[1,2],C))
    @test_throws DimensionMismatch materialize(MulAdd(2.0,A,B,3.0,randn(3,4)))
    @test_throws DimensionMismatch materialize(MulAdd(2.0,A,B,3.0,randn(5,5)))

    B = randn(5,5)
    C = randn(5,5)
    @test materialize(MulAdd(2.0,Diagonal(A),Diagonal(B),3.0,Diagonal(C))) ==
          materialize!(MulAdd(2.0,Diagonal(A),Diagonal(B),3.0,Diagonal(copy(C)))) == 2.0Diagonal(A)*Diagonal(B) + 3.0*Diagonal(C)
    @test_broken materialize(MulAdd(2.0,Diagonal(A),Diagonal(B),3.0,Diagonal(C))) isa Diagonal

    @test materialize(MulAdd(1.0, Eye(5), A, 3.0, C)) == materialize!(MulAdd(1.0, Eye(5), A, 3.0, copy(C))) == A + 3.0C
    @test materialize(MulAdd(1.0, A, Eye(5), 3.0, C)) == materialize!(MulAdd(1.0, A, Eye(5), 3.0, copy(C))) == A + 3.0C
end


text:
using ModelingToolkit: operation, istree, arguments, variable
import DomainSets

# Method of lines discretization scheme
struct MOLFiniteDifference{T,T2} <: DiffEqBase.AbstractDiscretization

end

MOLFiniteDifference(args...) = throw(ErrorException("MOLFiniteDifference has been moved to the MethodOfLines packgage, please use MethodOfLines.MOLFiniteDifference instead"))

# Constructors. If no order is specified, both upwind and centered differences will be 2nd order
MOLFiniteDifference(dxs, time; upwind_order = 1, centered_order = 2, grid_align=center_align) = throw(ErrorException("MOLFiniteDifference has been moved to the MethodOfLines packgage, please install with `Pkg.add(\"MethodOfLines\")` and use MethodOfLines.MOLFiniteDifference instead"))

function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discretization::DiffEqOperators.MOLFiniteDifference)
    throw(ErrorException("Symbolic Method of Lines discretization has been moved to the MethodOfLines packgage, please install this instead with `Pkg.add(\"MethodOfLines\")`"))
end

function SciMLBase.discretize(pdesys::ModelingToolkit.PDESystem,discretization::DiffEqOperators.MOLFiniteDifference)
    throw(ErrorException("Symbolic Method of Lines discretization has been moved to the MethodOfLines packgage, please install this instead with `Pkg.add(\"MethodOfLines\")`"))
end


text:
using Test
using Knet.Ops20: logp, logsumexp, softmax, ∇softmax, logsoftmax, ∇logsoftmax, nll, bce, logistic
using Knet.Ops20_gpu: _cudnnSoftmaxForward, _cudnnSoftmaxBackward
using Knet.KnetArrays: KnetArray
using CUDA: CUDA, functional
using CUDA.CUDNN: CUDNN_SOFTMAX_FAST, CUDNN_SOFTMAX_ACCURATE, CUDNN_SOFTMAX_LOG
using AutoGrad: grad, gradcheck, @gcheck, Param

@testset "loss" begin
    for f in (logp, logsumexp)
        a = rand(10,10)
        @test gradcheck(f,a)
        @test gradcheck(f,a,kw=(:dims=>1,))
        @test gradcheck(f,a,kw=(:dims=>2,))
        if CUDA.functional()
            k = KnetArray(a)
            @test gradcheck(f,k)
            @test gradcheck(f,k,kw=(:dims=>1,))
            @test gradcheck(f,k,kw=(:dims=>2,))
            @test isapprox(f(a),f(k))
            @test isapprox(f(a,dims=1),f(k,dims=1))
            @test isapprox(f(a,dims=2),f(k,dims=2))
        end

        a = rand(10,10,10)
        @test gradcheck(f,a)
        @test gradcheck(f,a,kw=(:dims=>1,))
        @test gradcheck(f,a,kw=(:dims=>2,))
        @test gradcheck(f,a,kw=(:dims=>3,))
        @test gradcheck(f,a,kw=(:dims=>(1,2),))
        @test gradcheck(f,a,kw=(:dims=>(3,2),))
        @test gradcheck(f,a,kw=(:dims=>(1,3),))

        if CUDA.functional()
            k = KnetArray(a)
            @test gradcheck(f,k)
            @test isapprox(f(a),f(k))
            for d in [1,2,3]
                @test gradcheck(f,k,kw=(:dims=>d,))
                @test isapprox(f(a,dims=d),f(k,dims=d))
            end
            for dims in [(1,2), (1,3), (3,1), [1,3,2]]
                @test isapprox(f(a,dims=dims),f(k,dims=dims))
            end
        end

        a = rand(10,10, 10)
        for d in [1, 2, 3, (1,2), (1,3), [2,3], [1,2,3]]
            @test softmax(a, dims=d) ≈ exp.(logsoftmax(a, dims=d))
            @test all(sum(softmax(a, dims=d), dims=d) .≈ 1)
            if CUDA.functional()
                k = KnetArray(a)
                @test softmax(k, dims=d) ≈ exp.(logsoftmax(k, dims=d))
                @test all(Array(sum(softmax(k, dims=d), dims=d)) .≈ 1)
            end
        end
    end

    a = randn(10,10)
    indices = rand(1:10,10)
    t01 = rand((0,1),100)
    t11 = rand((-1,1),100)
    @test gradcheck(nll, a, indices, kw=(:dims=>1,), args=1)
    @test gradcheck(nll, a, indices, kw=(:dims=>2,), args=1)
    @test gradcheck(logistic,vec(a),t11, args=1)
    @test gradcheck(bce,vec(a),t01, args=1)
    if CUDA.functional()
        k = KnetArray(a)
        @test gradcheck(nll, k, indices, kw=(:dims=>1,), args=1)
        @test gradcheck(nll, k, indices, kw=(:dims=>2,), args=1)
        @test gradcheck(logistic,vec(k),t11, args=1)
        @test gradcheck(bce,vec(k),t01, args=1)
        @test isapprox(nll(k, indices, dims=1), nll(a, indices, dims=1))
        @test isapprox(nll(k, indices, dims=2), nll(a, indices, dims=2))
        @test isapprox(logistic(vec(k),t11), logistic(vec(a),t11))
        @test isapprox(bce(vec(k),t01), bce(vec(a),t01))
    end

    # Issue 439: highorder derivatives
    x = randn(3,4); y1 = softmax(x,dims=1); y2 = logsoftmax(x,dims=1); dy = randn(3,4)
    @test @gcheck softmax(Param(x),dims=1)
    @test @gcheck ∇softmax(Param(x),Param(y1),Param(dy),dims=1)
    @test @gcheck logsoftmax(Param(x),dims=1)
    @test @gcheck ∇logsoftmax(Param(x),Param(y2),Param(dy),dims=1)
    if CUDA.functional()
        x = KnetArray(x); y1 = KnetArray(y1); y2 = KnetArray(y2); dy = KnetArray(dy)
        @test isapprox(softmax(x,dims=1), _cudnnSoftmaxForward(x,algo=CUDNN_SOFTMAX_FAST))
        @test isapprox(softmax(x,dims=1), _cudnnSoftmaxForward(x,algo=CUDNN_SOFTMAX_ACCURATE))
        @test isapprox(logsoftmax(x,dims=1), _cudnnSoftmaxForward(x,algo=CUDNN_SOFTMAX_LOG))
        @test isapprox(∇softmax(x,y1,dy,dims=1), _cudnnSoftmaxBackward(y1,dy,algo=CUDNN_SOFTMAX_FAST))
        @test isapprox(∇softmax(x,y1,dy,dims=1), _cudnnSoftmaxBackward(y1,dy,algo=CUDNN_SOFTMAX_ACCURATE))
        @test isapprox(∇logsoftmax(x,y2,dy,dims=1), _cudnnSoftmaxBackward(y2,dy,algo=CUDNN_SOFTMAX_LOG))
        @test @gcheck _cudnnSoftmaxForward(Param(x),algo=CUDNN_SOFTMAX_FAST)
        @test @gcheck _cudnnSoftmaxForward(Param(x),algo=CUDNN_SOFTMAX_ACCURATE)
        @test @gcheck _cudnnSoftmaxForward(Param(x),algo=CUDNN_SOFTMAX_LOG)
        @test @gcheck _cudnnSoftmaxBackward(Param(y1),Param(dy),algo=CUDNN_SOFTMAX_FAST)
        @test @gcheck _cudnnSoftmaxBackward(Param(y1),Param(dy),algo=CUDNN_SOFTMAX_ACCURATE)
        @test @gcheck _cudnnSoftmaxBackward(Param(y2),Param(dy),algo=CUDNN_SOFTMAX_LOG)

        # Broken example from Alkan's notebook:
        f(w,x,y) = nll(w*x,y)
        ∇f = grad(f)
        ∇fj(w,x,y,j) = ∇f(w,x,y)[j]
        ∇∇fj = grad(∇fj)
        a = rand(10,10); b = rand(10,10); c = rand(1:10,10)
        A = KnetArray(a); B = KnetArray(b); C = c
        @test isapprox(f(a,b,c), f(A,B,C))
        @test isapprox(∇f(a,b,c), ∇f(A,B,C))
        i = 10; j = 20; d = 1e-4
        @test isapprox(∇∇fj(a,b,c,i), ∇∇fj(A,B,C,i))
    end
end


text:
include("SLPDE.jl")
include("etdrk4.jl")
include("evaluate_lambda.jl")

srand(0)

println("Testing Local Allen--Cahn")

n = 25
U0 = zeros(4n, 8n-1);
U0[1:n, 1:2n-1] = sphrandn(Float64, n, n)/n

SLPDE = SemiLinearPDE(LaplaceBeltrami(0.1), NonlinearOperator(u->u-u^3), U0)

IMAX = 6

U = Vector{Matrix{Float64}}(IMAX)
for i = 1:IMAX
    U[i] = ETDRK4(SLPDE, 1.0, 10*2^(i-1))
    println("Done i = ", i)
end

err = zeros(IMAX-1)
for i = 1:IMAX-1
    err[i] = vecnorm(U[i] - U[end], Inf)
    println("i = ",i,", and the ∞-norm absolute error: ", err[i])
end

for i = 1:IMAX-2
    println(err[i]/err[i+1])
end

println("Testing Nonlocal Allen--Cahn")

SLPDE = SemiLinearPDE(NonlocalLaplaceBeltrami(0.1, -0.5, 0.2), NonlinearOperator(u->u-u^3), U0)

U = Vector{Matrix{Float64}}(IMAX)
for i = 1:IMAX
    U[i] = ETDRK4(SLPDE, 1.0, 10*2^(i-1))
    println("Done i = ", i)
end

err = zeros(IMAX-1)
for i = 1:IMAX-1
    err[i] = vecnorm(U[i] - U[end], Inf)
    println("i = ",i,", and the ∞-norm absolute error: ", err[i])
end

for i = 1:IMAX-2
    println(err[i]/err[i+1])
end

println("Testing Local Gray--Scott")

n = 25
U0 = zeros(4n, 8n-1), zeros(4n, 8n-1);
U0[1][1:n, 1:2n-1] = sphrandn(Float64, n, n)/n
U0[2][1:n, 1:2n-1] = sphrandn(Float64, n, n)/n

SLPDE = SemiLinearPDE((LaplaceBeltrami(sqrt(0.00002)), LaplaceBeltrami(sqrt(0.00001))), (NonlinearOperator((u,v)->0.04*(1-u)-u*v^2), NonlinearOperator((u,v)->-0.1*v+u*v^2)), U0)

U = Vector{Tuple{Matrix{Float64},Matrix{Float64}}}(IMAX)
for i = 1:IMAX
    U[i] = ETDRK4(SLPDE, 1.0, 10*2^(i-1))
    println("Done i = ", i)
end

err = zeros(IMAX-1)
for i = 1:IMAX-1
    err[i] = vecnorm(U[i][1] - U[end][1], Inf) + vecnorm(U[i][2] - U[end][2], Inf)
    println("i = ",i,", and the ∞-norm absolute error: ", err[i])
end

for i = 1:IMAX-2
    println(err[i]/err[i+1])
end

println("Testing Nonlocal Gray--Scott")

SLPDE = SemiLinearPDE((NonlocalLaplaceBeltrami(sqrt(0.00002), -0.5, 0.2), NonlocalLaplaceBeltrami(sqrt(0.00001), -0.5, 0.2)), (NonlinearOperator((u,v)->0.04*(1-u)-u*v^2), NonlinearOperator((u,v)->-0.1*v+u*v^2)), U0)

U = Vector{Tuple{Matrix{Float64},Matrix{Float64}}}(IMAX)
for i = 1:IMAX
    U[i] = ETDRK4(SLPDE, 1.0, 10*2^(i-1))
    println("Done i = ", i)
end

err = zeros(IMAX-1)
for i = 1:IMAX-1
    err[i] = vecnorm(U[i][1] - U[end][1], Inf) + vecnorm(U[i][2] - U[end][2], Inf)
    println("i = ",i,", and the ∞-norm absolute error: ", err[i])
end

for i = 1:IMAX-2
    println(err[i]/err[i+1])
end


text:
using LinearAlgebraicRepresentation
using Plasm
using DataStructures
import SparseArrays.spzeros
Lar = LinearAlgebraicRepresentation

v0 = [0. 1. 2. 3. 4. 5. 6. 7. 8.]

v1 = collect(0.:12)

c0 = Lar.larGridSkeleton([1,1,1])(0)

m0 = (v0, c0)

c1 = [[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8]]
v1 = hcat(collect(0.:12)...)

m0 = (v0,c0)

m1 = (v1,c1)

Plasm.view(Lar.larModelProduct( m1,m0 ))

m11 = Lar.larModelProduct( m0, m1 )

m12 = Lar.larModelProduct( m1, m0 )

Lar.Struct([ m11,m12 ])

m2 = Lar.struct2lar(Lar.Struct([ m11,m12 ]))

Plasm.view(m2)

m21 = Lar.larModelProduct( m2,m1 )

Plasm.view(m21)

m20 = Lar.larModelProduct( m1,m1 )

Plasm.view(m20)

m20 = Lar.larModelProduct( m20,m0 )

m20 = Lar.larModelProduct( m1,m1 )

m22 = Lar.larModelProduct( m20,m0 )

Plasm.view(m22)

Lar.struct2lar(Lar.Struct([ m22,m21 ]))

m32 = Lar.struct2lar(Lar.Struct([ m22,m21 ]))

Plasm.view(m32)

println(typeof(m32))

Vm32, FVm32 = m32

#VV
VVm32 = [[i] for i=1:size(Vm32)[2]]
Plasm.view( Plasm.numbering(.6)((Vm32,[VVm32, FVm32])) )

# Plasm.view( Plasm.numbering(.6)((V,[VV, EV])) )

# Plasm.view( Plasm.numbering(0.6)((m32[1], [m32[2]])))


text:
# This file was generated, do not modify it. # hide
post = m(; x, λ) | (; y)

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
include(joinpath(pwd(), "src/models/pendulum.jl"))
include(joinpath(pwd(), "src/constraints/free_time.jl"))

optimize = true

# Free-time model with additive noise
model_ft = Pendulum(2, 2, 2, 1.0, 0.1, 0.5, 0.25, 9.81)

function fd(model::Pendulum, x⁺, x, u, w, h, t)
    midpoint_implicit(model, x⁺, x, u, w, u[end]) - w
end

# Horizon
T = 51

# Bounds
tf = 2.0
h0 = tf / (T-1) # timestep

ul, uu = control_bounds(model_ft, T, [-3.0; 0.0], [3.0; h0])

# Initial and final states
x1 = [0.0; 0.0]
xT = [π; 0.0]

xl, xu = state_bounds(model_ft, T, x1 = x1, xT = xT)

# Objective (minimum time)
obj = PenaltyObjective(1.0, model_ft.m)

# Time step constraints
con_free_time = free_time_constraints(T)

# Problem
prob = trajectory_optimization_problem(model_ft,
			   obj,
			   T,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
			   con = con_free_time
               )

# Initialization
X0 = linear_interp(x1, xT, T) # linear interpolation for states
U0 = [ones(model_ft.m) for t = 1:T-1]

# Pack trajectories into vector
Z0 = pack(X0, U0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
if optimize
    @time Z̄ = solve(prob, copy(Z0))
    @save joinpath(@__DIR__, "sol.jld2") Z̄
else
    println("Loading solution...")
    @load joinpath(@__DIR__, "sol.jld2") Z̄
end

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
# p29.jl - solve Poisson equation on the unit disk
#         (compare p16.jl and p28.jl)

# Laplacian in polar coordinates:
N = 25; N2 = Int((N-1)/2);
(D,r) = cheb(N); D2 = D^2;
D1 = D2[2:N2+1,2:N2+1]; D2 = D2[2:N2+1,N:-1:N2+2];
E1 =  D[2:N2+1,2:N2+1]; E2 =  D[2:N2+1,N:-1:N2+2];
M = 20; dt = 2*pi/M; t = dt*(1:M); M2 = Int(M/2);
D2t = toeplitz([-pi^2/(3*dt^2)-1/6; @. .5*(-1)^(2:M)/sin(dt*(1:M-1)/2)^2]);
R = diagm(1./r[2:N2+1]); Z = zeros(M2,M2); I = eye(M2);
L = kron(D1+R*E1,eye(M)) + kron(D2+R*E2,[Z I;I Z]) + kron(R^2,D2t);

# Right-hand side and solution for u:
(rr,tt) = (r[2:N2+1]',t);
f = @. -rr^2*sin(tt/2)^4 + sin(6*tt)*cos(tt/2)^2; u = L\f[:];

# Reshape results onto 2D grid and plot them:
u = reshape(u,M,N2); u = [zeros(M+1) u[[M;1:M],:]];
(rr,tt) = (r[1:N2+1],t[[M;1:M]]);
(xx,yy) = @. (cos(tt)*rr',sin(tt)*rr');
clf();
surf(xx,yy,u);
gca()[:view_init](40,-70);
xlim(-1,1); ylim(-1,1); zlim(-.01,.05);
xlabel("x"); ylabel("y"); zlabel("z");


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

using ValueShapes
using Test

using Distributions
using Random

@testset "const_value_dist" begin
    @test @inferred(ConstValueDist(4.2)) isa Distribution{Univariate}
    @test @inferred(ConstValueDist([1, 2, 3])) isa Distribution{Multivariate}
    @test @inferred(ConstValueDist([1 2; 3 4])) isa Distribution{Matrixvariate}

    @test @inferred(length(ConstValueDist([4.2]))) == 1
    @test @inferred(length(ConstValueDist([1, 2, 3]))) == 3
    @test @inferred(length(ConstValueDist([1 2; 3 4]))) == 4

    @test @inferred(pdf(ConstValueDist(4.2), 4.2)) == Inf
    @test @inferred(logpdf(ConstValueDist(4.2), 4.2)) == Inf
    @test @inferred(pdf(ConstValueDist(4.2), 3.7)) == 0
    @test @inferred(logpdf(ConstValueDist(4.2), 3.7)) == -Inf
    @test @inferred(pdf(ConstValueDist(4.2), [4.2, 3.7])) == [Inf, 0.0]
    @test @inferred(logpdf(ConstValueDist(4.2), [4.2, 3.7])) == [Inf, -Inf]

    @test @inferred(pdf(ConstValueDist([1, 2, 3]), [1, 2, 3])) == Inf
    @test @inferred(logpdf(ConstValueDist([1, 2, 3]), [1, 2, 3])) == Inf
    @test @inferred(pdf(ConstValueDist([1, 2, 3]), [2, 3, 4])) == 0
    @test @inferred(logpdf(ConstValueDist([1, 2, 3]), [2, 3, 4])) == -Inf
    @test (pdf(ConstValueDist([1, 2, 3]), hcat([1, 2, 3], [2, 3, 4]))) == [Inf, 0.0]
    @test (logpdf(ConstValueDist([1, 2, 3]), hcat([1, 2, 3], [2, 3, 4]))) == [Inf, -Inf]

    @test @inferred(pdf(ConstValueDist([1 2; 3 4]), [1 2; 3 4])) == Inf
    @test @inferred(logpdf(ConstValueDist([1 2; 3 4]), [1 2; 3 4])) == Inf
    @test @inferred(pdf(ConstValueDist([1 2; 3 4]), [4 5; 6 7])) == 0
    @test @inferred(logpdf(ConstValueDist([1 2; 3 4]), [4 5; 6 7])) == -Inf

    @test @inferred(insupport(ConstValueDist(4.2), 4.2)) == true
    @test @inferred(insupport(ConstValueDist(4.2), 4.19)) == false

    @test @inferred(insupport(ConstValueDist([1, 2, 3]), [1, 2, 3])) == true
    @test @inferred(insupport(ConstValueDist([1, 2, 3]), [2, 3, 4])) == false

    @test @inferred(insupport(ConstValueDist([1 2; 3 4]), [1 2; 3 4])) == true
    @test @inferred(insupport(ConstValueDist([1 2; 3 4]), [4 5; 6 7])) == false

    @test @inferred(rand(ConstValueDist(4.2))) == 4.2
    @test @inferred(rand(ConstValueDist(4.2), 3)) == [4.2, 4.2, 4.2]

    @test @inferred(rand(ConstValueDist([1, 2, 3]))) == [1, 2, 3]
    @test @inferred(rand(ConstValueDist([1, 2, 3]), 2)) == hcat([1, 2, 3], [1, 2, 3])

    @test @inferred(rand(ConstValueDist([1 2; 3 4]))) == [1 2; 3 4]
    @test @inferred(rand(ConstValueDist([1 2; 3 4]), 2)) == [[1 2; 3 4], [1 2; 3 4]]

    univariate_cvd = @inferred(ConstValueDist(Int64(42)))

    shape = varshape(univariate_cvd)

    @test @inferred(totalndof(shape)) == 0
    @test @inferred(size(univariate_cvd)) == ()
    @test @inferred(length(univariate_cvd)) == 1

    @test @inferred(minimum(univariate_cvd)) == 42
    @test @inferred(maximum(univariate_cvd)) == 42

    @test @inferred(pdf(univariate_cvd, 42)) == Inf

    @test @inferred(cdf(univariate_cvd, 41.999)) == 0
    @test @inferred(cdf(univariate_cvd, 42)) == 1
    @test @inferred(cdf(univariate_cvd, Inf)) == 1

    @test @inferred(mean(univariate_cvd)) == 42

    @test @inferred(mode(univariate_cvd)) == 42

    @test @inferred(eltype(univariate_cvd)) == Int64

    μ1, μ2 = rand(2)
    cvd_from_named_tuple = @inferred(ConstValueDist((a=μ1, b=μ2)))

    @test @inferred(log(@inferred(pdf(cvd_from_named_tuple, (a=μ1, b=μ2))))) == @inferred(logpdf(cvd_from_named_tuple, (a=μ1, b=μ2))) == Inf

    @test @inferred(insupport(cvd_from_named_tuple, (a=μ1, b=μ2))) == true
    @test @inferred(insupport(cvd_from_named_tuple, (a=μ1+eps(μ1), b=μ2+eps(μ2)))) == false
    @test @inferred(insupport(cvd_from_named_tuple, (a=μ1+eps(μ1), b=μ2))) == false
    @test @inferred(insupport(cvd_from_named_tuple, (a=μ1, b=μ2+eps(μ2)))) == false

    n_samples = 100
    samples = @inferred(rand(cvd_from_named_tuple, n_samples))
    emptied_samples = @inferred(similar(samples))

    @test emptied_samples != samples
    @test samples == @inferred(fill((a=μ1, b=μ2), n_samples))
    @test @inferred(rand!(cvd_from_named_tuple, emptied_samples)) == samples
    @test samples == emptied_samples
end


text:
using FreeTypeAbstraction: height_insensitive_boundingbox

"""
Calculates the exact boundingbox of a Scene/Plot, without considering any transformation
"""
raw_boundingbox(x::Atomic) = data_limits(x)

rootparent(x) = rootparent(parent(x))
rootparent(x::Scene) = x

# function raw_boundingbox(x::Annotations)
#     bb = raw_boundingbox(x.plots)
#     inv(modelmatrix(rootparent(x))) * bb
# end

raw_boundingbox(x::Combined) = raw_boundingbox(x.plots)
boundingbox(x) = raw_boundingbox(x)

function combined_modelmatrix(x)
    m = Mat4f0(I)
    while true
        m = modelmatrix(x) * m
        if parent(x) !== nothing && parent(x) isa Combined
            x = parent(x)
        else
            break
        end
    end
    return m
end

function modelmatrix(x)
    t = transformation(x)
    transformationmatrix(t.translation[], t.scale[], t.rotation[])
end

function boundingbox(x::Atomic)
    bb = raw_boundingbox(x)
    return combined_modelmatrix(x) * bb
end

boundingbox(scene::Scene) = raw_boundingbox(scene)
function raw_boundingbox(scene::Scene)
    if scene[OldAxis] !== nothing
        return raw_boundingbox(scene[OldAxis])
    elseif scene.limits[] !== automatic
        return scene_limits(scene)
    elseif cameracontrols(scene) == EmptyCamera()
        # Empty camera means this is a parent scene that itself doesn't display anything
        return raw_boundingbox(scene.children)
    else
        plots = plots_from_camera(scene)
        children = filter(scene.children) do child
            child.camera == scene.camera
        end
        return raw_boundingbox([plots; children])
    end
end

function raw_boundingbox(plots::Vector)
    isempty(plots) && return FRect3D()
    plot_idx = iterate(plots)
    bb = FRect3D()
    while plot_idx !== nothing
        plot, idx = plot_idx
        plot_idx = iterate(plots, idx)
        # isvisible(plot) || continue
        bb2 = boundingbox(plot)
        isfinite_rect(bb) || (bb = bb2)
        isfinite_rect(bb2) || continue
        bb = union(bb, bb2)
    end
    return bb
end

function project_widths(matrix, vec)
    pr = project(matrix, vec)
    zero = project(matrix, zeros(typeof(vec)))
    return pr - zero
end

function rotate_bbox(bb::FRect3D, rot)
    points = decompose(Point3f0, bb)
    FRect3D(Ref(rot) .* points)
end

function gl_bboxes(gl::GlyphCollection)
    scales = gl.scales.sv isa Vec2f0 ? (gl.scales.sv for _ in gl.extents) : gl.scales.sv
    map(gl.extents, gl.fonts, scales) do ext, font, scale
        unscaled_hi_bb = height_insensitive_boundingbox(ext, font)
        hi_bb = FRect2D(
            Makie.origin(unscaled_hi_bb) * scale,
            widths(unscaled_hi_bb) * scale
        )
    end
end

function boundingbox(glyphcollection::GlyphCollection, position::Point3f0, rotation::Quaternion)

    if isempty(glyphcollection.glyphs)
        return FRect3D(position, Vec3f0(0, 0, 0))
    end

    chars = glyphcollection.glyphs
    glyphorigins = glyphcollection.origins
    glyphbbs = gl_bboxes(glyphcollection)

    bb = FRect3D()
    for (char, charo, glyphbb) in zip(chars, glyphorigins, glyphbbs)
        # ignore line breaks
        # char in ('\r', '\n') && continue

        charbb = rotate_bbox(FRect3D(glyphbb), rotation) + charo + position
        if !isfinite_rect(bb)
            bb = charbb
        else
            bb = union(bb, charbb)
        end
    end
    !isfinite_rect(bb) && error("Invalid text boundingbox")
    bb
end

function boundingbox(layouts::AbstractArray{<:GlyphCollection}, positions, rotations)

    if isempty(layouts)
        FRect3D((0, 0, 0), (0, 0, 0))
    else
        bb = FRect3D()
        broadcast_foreach(layouts, positions, rotations) do layout, pos, rot
            if !isfinite_rect(bb)
                bb = boundingbox(layout, pos, rot)
            else
                bb = union(bb, boundingbox(layout, pos, rot))
            end
        end
        !isfinite_rect(bb) && error("Invalid text boundingbox")
        bb
    end
end

function boundingbox(x::Text{<:Tuple{<:GlyphCollection}})
    boundingbox(
        to_value(x[1]),
        to_ndim(Point3f0, to_value(x[:position]), 0),
        convert_attribute(to_value(x[:rotation]), key"rotation"())
    )
end

function boundingbox(x::Text{<:Tuple{<:AbstractArray{<:GlyphCollection}}})
    boundingbox(
        to_value(x[1]),
        to_ndim.(Point3f0, to_value(x[:position]), 0),
        convert_attribute(to_value(x[:rotation]), key"rotation"())
    )
end

#     if x.space[] == :data
#         if x[1][] isa AbstractArray
#             bb = FRect3D()
#             for
#         else
#             boundingbox(x, to_value(x[1]))
#         end
#     elseif x.space[] == :screen
#         data_limits(x)
#     end
# end

# boundingbox(x::Text) = boundingbox(x, to_value(x[1]))

# function boundingbox(
#         text::String, position, textsize;
#         font = "default", align = (:left, :bottom), rotation = 0.0
#     )
#     return boundingbox(
#         text, position, textsize,
#         to_font(font), to_align(align), to_rotation(rotation)
#     )
# end

"""
Calculate an approximation of a tight rectangle around a 2D rectangle rotated by `angle` radians.
This is not perfect but works well enough. Check an A vs X to see the difference.
"""
function rotatedrect(rect::Rect{2}, angle)
    ox, oy = rect.origin
    wx, wy = rect.widths
    points = @SMatrix([
        ox oy;
        ox oy+wy;
        ox+wx oy;
        ox+wx oy+wy;
    ])
    mrot = @SMatrix([
        cos(angle) -sin(angle);
        sin(angle) cos(angle);
    ])
    rotated = mrot * points'

    rmins = minimum(rotated, dims = 2)
    rmaxs = maximum(rotated, dims = 2)

    return Rect2D(rmins..., (rmaxs .- rmins)...)
end

function quaternion_to_2d_angle(quat)
    # this assumes that the quaternion was calculated from a simple 2d rotation as well
    return 2acos(quat[4]) * (signbit(quat[1]) ? -1 : 1)
end

function boundingbox(
        text::String, position, textsize, font,
        align, rotation, model, justification, lineheight;
        # use the font's ascenders and descenders for the bounding box
        # this means that a string's boundingbox doesn't change in the vertical
        # dimension when characters change (for example numbers during an animation)
        # this is not wanted in most cases because of the jitter it creates when
        # the boundingbox slightly changes size in each frame (in MakieLayout mostly)
        use_vertical_dimensions_from_font = true
    )
    atlas = get_texture_atlas()
    N = length(text)
    ctext_state1 = iterate(text)
    ctext_state1 === nothing && return FRect3D()

    # call the layouting algorithm to find out where all the glyphs end up
    # this is kind of a doubling, maybe it could be avoided if at creation all
    # positions would be populated in the text object, but that seems convoluted
    if position isa VecTypes
        glyph_offsets = layout_text(text, textsize, font, align,
            rotation, model, justification, lineheight)
        position = [go .+ to_ndim(Point3f0, position, 0) for go in glyph_offsets]
    end

    bbox = Ref(FRect3D())            # use ref to avoid julia#15276
    ctext_state = Ref{Union{Nothing,typeof(ctext_state1)}}(ctext_state1)

    let position = position
        broadcast_foreach(1:N, rotation, font, textsize) do i, rotation, font, scale
            c, text_state = ctext_state[]::typeof(ctext_state1)
            ctext_state[] = iterate(text, text_state)

            if !(c in ('\r', '\n'))
                bb_unitspace = if use_vertical_dimensions_from_font
                    height_insensitive_boundingbox(
                        FreeTypeAbstraction.get_extent(font, c), font)
                else
                    inkboundingbox(FreeTypeAbstraction.get_extent(font, c))
                end

                scaled_bb = bb_unitspace * scale

                # TODO this only works in 2d
                rot_2d_radians = quaternion_to_2d_angle(rotation)
                rotated_bb = rotatedrect(scaled_bb, rot_2d_radians)

                # bb = rectdiv(bb, 1.5)
                shifted_bb = FRect3D(rotated_bb) + position[i]
                bboxc = bbox[]
                if !isfinite_rect(bboxc)
                    bbox[] = shifted_bb
                else
                    bbox[] = union(bboxc, shifted_bb)
                end
            end
        end
    end
    return bbox[]
end


text:
#!/usr/bin/env julia

using Gadfly
using RDatasets

tests = [
    "points",
    "colored_points",
    "vstack",
    "hstack",
    "stacked_discrete_histogram",
    "dodged_discrete_histogram",
    "array_aesthetics"
]

backends = {
    "svg" => name -> SVG("$(name).svg", 6inch, 3inch),
    "d3"  => name -> D3("$(name).js",   6inch, 3inch),
    "png" => name -> PNG("$(name).png", 6inch, 3inch),
    "ps"  => name -> PS("$(name).ps",   6inch, 3inch),
    "pdf" => name -> PDF("$(name).pdf", 6inch, 3inch)
}

function run_tests()
    for test in tests
        for (backend_name, backend) in backends
            println(STDERR, "Rendering $(test) on $(backend_name) backend.")
            p = evalfile("$(test).jl")
            draw(backend(test), p)
        end
    end
end

run_tests()


text:
using Documenter
using CatmullClark

makedocs(
    sitename="CatmullClark.jl",
    modules=[CatmullClark],
    pages=["index.md", "examples.md", "reference.md", "devnotes.md"],
)

deploydocs(
    repo="github.com/wherrera10/CatmullClark.jl.git",
)


text:
using Compose
using Colors
using DataFrames
using Random: shuffle!

const TernBack = RGB(253 / 255, 253 / 255, 241 / 255)

const TernPalette = distinguishable_colors(66, Color[TernBack, RGB(0, 0, 0), RGB(0 / 255, 168 / 255, 45 / 255)])[3:end]
const TernColorblind =
    distinguishable_colors(66, Color[TernBack, RGB(0, 0, 0), colorant"DodgerBlue4"], transform = deuteranopic)[3:end]

"""
    multiternary(data::DataFrame, cols::AbstractArray{Symbol}, clscol::Union{Symbol, Missing}=missing; palette=TernPalette, norm=1.0)

Uses Compose.jl to draw a simple multi-ternary diagram from the `data`, a `DataFrame`.
The columns are drawn in the order in `cols` and the vertices labeled with the names
of the columns.
"""
function multiternary(
    data::DataFrame,
    cols::AbstractArray,
    clscol::Any = missing;
    title = missing,
    palette = TernPalette,
    font = "Verdana",
    fontsize = 5.0,
    norm = 1.0,
    withcount = true,
)
    # Draw an individual ternary diagram for cola and colb in data
    function ternary(data, cola, colb, v1, L, θ)
        # Compute the intersection of the line (P1:P2) with (P3:P4)
        function xy(P1, P2, P3, P4)
            s =
                (P1[2] * (P2[1] - P3[1]) + P2[2] * P3[1] - P2[1] * P3[2] + P1[1] * (P3[2] - P2[2])) /
                (P2[2] * (P3[1] - P4[1]) + P1[2] * (P4[1] - P3[1]) + (P1[1] - P2[1]) * (P3[2] - P4[2]))
            return (P3[1] + s * (P4[1] - P3[1]), P3[2] + s * (P4[2] - P3[2]))
        end
        # From a, b on [0,1] compute P1, P2, P3 and P4
        function cc(a, b, va, vb, vc)
            return xy(vb .+ a .* (va .- vb), vc .+ a .* (va .- vc), va .+ b .* (vb .- va), vc .+ b .* (vb .- vc))
        end
        v3 = (v1[1] + L * cos(θ - π / 6), v1[2] - L * sin(θ - π / 6))
        v2 = (v1[1] + L * cos(θ + π / 6), v1[2] - L * sin(θ + π / 6))
        # 0.1 increment lines
        vx = [[v1 .+ f .* (v2 .- v1), v1 .+ f .* (v3 .- v1)] for f = 0.1:0.1:0.91]
        append!(vx, [[v2 .+ f .* (v3 .- v2), v2 .+ f .* (v1 .- v2)] for f = 0.1:0.1:0.91])
        append!(vx, [[v3 .+ f .* (v1 .- v3), v3 .+ f .* (v2 .- v3)] for f = 0.1:0.1:0.91])
        # Data points
        datapts, labels, range = [], [], ismissing(clscol) ? (1:1) : levels(data[:, clscol])
        ilvl = 1
        for lvl in range
            f(row) = ismissing(clscol) || (row[clscol] == lvl)
            pts = [cc(r[cola] / norm, r[colb] / norm, v3, v2, v1) for r in filter(f, eachrow(data))]
            if length(pts) > 0
                pa, pb, col = [pt[1] for pt in pts], [pt[2] for pt in pts], palette[ilvl]
                push!(
                    datapts,
                    (
                        context(),
                        Compose.circle(pa, pb, [L * 0.003]),
                        Compose.stroke(col),
                        Compose.fill("transparent"),
                        Compose.clip([v3, v2, v1]),
                    ),
                )
                if !ismissing(clscol)
                    push!(
                        labels,
                        (
                            context(),
                            Compose.text(0.02, 0.54 + ilvl * 0.03, withcount ? "$lvl - $(length(pts))" : "$lvl"),
                            Compose.fill(col),
                            Compose.font(font),
                            Compose.fontsize(0.8 * fontsize),
                        ),
                    )
                end
                ilvl = (ilvl % length(palette)) + 1
            end
        end
        shuffle!(datapts) # randomize datapoint order
        if length(labels) > 15
            last = labels[end]
            labels = labels[1:12]
            push!(
                labels,
                (last[1], Compose.text(0.02, 0.54 + 13 * 0.03, "…"), Compose.fill(colorant"black"), last[4], last[5]),
            )
            push!(
                labels,
                (last[1], Compose.text(0.02, 0.54 + 14 * 0.03, last[2].primitives[1].value), last[3], last[4], last[5]),
            )
        end
        return (
            context(),
            (context(), datapts...),
            (context(), Compose.polygon([v1, v2, v3]), Compose.stroke("black"), Compose.fill("transparent")),
            (context(), Compose.line(vx), Compose.stroke("darkgray")),
            (context(), Compose.polygon([v1, v2, v3]), Compose.stroke("transparent"), Compose.fill(TernBack)),
            (context(), labels...),
        )
    end
    # Label the ends of the axes
    function label(v0, labels, L)
        xl, yl, ll = Float64[], Float64[], String[]
        function ss(l)
            if startswith(l,"Element(")
                try
                    parse(Element,l[9:end-1]).symbol
                catch
                    titlecase(l)
                end
            else
                titlecase(l)
            end
        end
        for i in eachindex(labels)
            θ = -π + π * i / 3
            vv = v(v0, θ, L)
            push!(xl, vv[1])
            push!(yl, vv[2])
            push!(ll, ss(labels[i]))
        end
        return (
            context(),
            Compose.text(xl, yl, ll, [hcenter], [vcenter]),
            Compose.fill("black"),
            Compose.font(font),
            Compose.fontsize(fontsize),
        )
    end
    v(v0, θ, off = 0.02) = (v0[1] + off * cos(θ), v0[2] - off * sin(θ))
    center, ops, ncols = (0.5, 0.5), [], min(length(cols), 6)
    for i = 2:ncols
        sa, sb = cols[i-1], cols[i]
        θ = -π / 2 + π / 3 * (i - 2)
        append!(ops, ternary(data, sa, sb, v(center, θ), 0.42, θ))
    end
    if !ismissing(title)
        append!(
            ops,
            (
                context(),
                Compose.text(0.5, 0.03, title, hcenter, vtop),
                Compose.fill("black"),
                Compose.font(font),
                Compose.fontsize(1.3 * fontsize),
            ),
        )
    end
    Compose.compose(context(0.05, 0.0, 0.9, 0.9), (context(), ops...), (context(), label(center, string.(cols[1:ncols]), 0.48)))
end

"""
    multiternary(
        zep::Zeppelin;
        rows = missing,
        omit = [ n"C", n"O" ],
        palette = TernPalette,
        fontsize = 12pt,
        font = "Verdana",
    )

Plot the elemental data in `zep` to a multi-ternary diagram.
"""
function multiternary(
        zep::Zeppelin;
        rows = missing,
        omit = [ n"C", n"O" ],
        palette = TernPalette,
        fontsize = 12pt,
        font = "Verdana",
)
    # Determine which elements to plot...
    zd = ismissing(rows) ? zep.data : zep.data[rows, :]
    df=DataFrame(Elm=Symbol[], Mean=Float64[])
    foreach(elm->push!(df, [ convert(Symbol,elm), mean(zd[:,convert(Symbol, elm)])]),  filter(elm->!(elm in omit), zep.elms))
    sort!(df, :Mean, rev=true)
    elms = df[:,:Elm][1:min(size(df,1),6)]
    multiternary(zd, elms, :CLASS, title=zep.header["DESCRIPTION"], palette=palette, norm=100.0, fontsize=fontsize, font=font)
end

max_stars_repo_name: NicholasWMRitchie/NeXLParticle

max_issues_repo_name: NicholasWMRitchie/NeXLParticle

max_forks_repo_name: NicholasWMRitchie/NeXLParticle