text:
module Visualize_GeometriesWithMaterial

using Modia3D

filename = joinpath(Modia3D.path, "objects", "fish", "SiameseTiger0.3ds")
#filename = joinpath(Modia3D.path, "objects","engine", "crank", "crank.obj")
#file = FileShape.convexFile(filename; scaleFactor=MVector{3,Float64}(4.0,4.0,4.0))

# Properties
mat1  = Modia3D.Material(color="LightBlue", transparency=0.5)
mat2  = Modia3D.Material(transparency=0.1)
font  = Modia3D.Font(fontFamily="TimesNewRoman", bold=false, charSize=0.2, color="Black")
rtext = [0.0, 0.0, 1.0]
ltext = [0.0, 0.0, -0.7]

# Objects3D
@assembly GeometriesWithMaterial begin
   world            = Modia3D.Object3D(visualizeFrame=false)
   box              = Modia3D.Object3D(world, Modia3D.Box(0.9,0.5,0.3,       material=mat1); r=[ 4.5,0.0, 0.0])
   sphere           = Modia3D.Object3D(world, Modia3D.Sphere(0.7,            material=mat1); r=[ 3.0,0.0, 0.0])
   cylinder         = Modia3D.Object3D(world, Modia3D.Cylinder(0.5,0.8,      material=mat1); r=[ 1.5,0.0, 0.0])
   cone             = Modia3D.Object3D(world, Modia3D.Cone(0.3,0.7,          material=mat1); r=[ 0.0,0.0, 0.0])
   capsule          = Modia3D.Object3D(world, Modia3D.Capsule(0.4,0.45,      material=mat1); r=[-1.5,0.0, 0.0])
   spring           = Modia3D.Object3D(world, Modia3D.Spring(0.3,0.7,        material=mat1); r=[ 4.5,0.0,-2.5])
   gearWheel        = Modia3D.Object3D(world, Modia3D.GearWheel(0.5,0.8,     material=mat1); r=[ 3.0,0.0,-2.5])
   pipe             = Modia3D.Object3D(world, Modia3D.Pipe(0.5,0.8,          material=mat1); r=[ 1.5,0.0,-2.5])
   beam             = Modia3D.Object3D(world, Modia3D.Beam(0.4,0.5,0.3,      material=mat1); r=[ 0.0,0.0,-2.5])
   fileMesh         = Modia3D.Object3D(world, Modia3D.FileMesh(filename,4.0, material=mat2); r=[-1.5,0.0,-2.5])

   # Place text above the shapes
   boxText              = Modia3D.Object3D(box             , Modia3D.TextShape("Box"             ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   sphereText           = Modia3D.Object3D(sphere          , Modia3D.TextShape("Sphere"          ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   cylinderText         = Modia3D.Object3D(cylinder        , Modia3D.TextShape("Cylinder"        ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   coneText             = Modia3D.Object3D(cone            , Modia3D.TextShape("Cone"            ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   capsuleText          = Modia3D.Object3D(capsule         , Modia3D.TextShape("Capsule"         ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   springText           = Modia3D.Object3D(spring          , Modia3D.TextShape("Spring"          ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   gearWheelText        = Modia3D.Object3D(gearWheel       , Modia3D.TextShape("GearWheel"       ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   pipeText             = Modia3D.Object3D(pipe            , Modia3D.TextShape("Pipe"            ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   beamText             = Modia3D.Object3D(beam            , Modia3D.TextShape("Beam"            ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   fileMeshText         = Modia3D.Object3D(fileMesh        , Modia3D.TextShape("FileMesh"        ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   fileMeshTextSb       = Modia3D.Object3D(fileMesh        , Modia3D.TextShape("(.obj, .3ds, .dxf, .stl)"        ; font=font); r=ltext, visualizeFrame=Modia3D.False)
end

Modia3D.visualizeAssembly!( GeometriesWithMaterial(sceneOptions = Modia3D.SceneOptions(visualizeFrames=true,
                                                                                defaultFrameLength=0.7,
                                                                                enableContactDetection=false)) )

println("... success of Visualize_GeometriesWithMaterial.jl!")

end


text:
function update!(mc::MarkovChain{T, M},
                 mcmc::MCMC{T, M},
                 Σ::Array{Float64, 2},
                 σ::Float64;
                 condition_on_network::Bool=false,
                 event_batches::Int64=1,
                 transmission_rate_cache::Union{Nothing, TransmissionRateCache}=nothing) where {
                 T <: DiseaseStateSequence,
                 M <: TNILM}
  # Initialize
  new_tr_cache = transmission_rate_cache
  new_events = mc.events[end]
  new_events_array = new_events[convert(DiseaseStates, T)[2:end]]
  new_rp = mc.risk_parameters[end]
  new_network = mc.transmission_network[end]
  new_lposterior = mc.log_posterior[end]
  # Randomize event time augmentation
  event_indices = findall(new_events_array[:] .> -Inf)
  aug_order = sample(event_indices, length(event_indices), replace=false)
  @debug "$(length(aug_order)) events to augment"
  if event_batches < 0
    @error "Cannot have negative amount of event batches"
  end
  if event_batches > length(aug_order)
    @warn "More event batches than there are events to augment ($(event_batches) > $(length(aug_order))), setting to maximum ($(length(aug_order)))"
    event_batches = length(aug_order)
  end
  batch_size = fld(length(aug_order), event_batches)
  @debug "Performing data augmentation in batches of $batch_size events at a time"
  for i = 1:(1 + event_batches)
    if i == 1
      # Propose new risk parameters
      proposed_events = new_events
      proposed_events_array = new_events_array
      proposed_rp = generate(RiskParameters{T}, new_rp, Σ)
      @debug "Generating new TransmissionRateCache for TN-ILM parameter proposal"
      proposed_tr_cache = TransmissionRateCache(size(new_events_array, 1))
    else
      proposed_events_array = copy(new_events_array)
      proposed_events = Events{T}(proposed_events_array)
      for j = (batch_size*(i-2) + 1):min(batch_size*(i-1),length(aug_order))
        id, state_index = Tuple(CartesianIndices((individuals(new_events),
                                            length(convert(DiseaseStates, T)[2:end])))[aug_order[j]])
        new_state = convert(DiseaseStates, T)[state_index+1]
        time = new_events[new_state][id]

        # Conditioning on network means that only event times valid under current network will be proposed.
        # This is useful for models which may have additional contributions to the posterior based on network,
        # and require more modest proposals (e.g. phylodynamic models).
        if condition_on_network
          proposed_event = generate(Event,
                                    Event{T}(time, id, new_state),
                                    σ,
                                    mcmc.event_extents,
                                    mcmc.event_observations,
                                    proposed_events,
                                    new_network)
        else
          proposed_event = generate(Event,
                                    Event{T}(time, id, new_state),
                                    σ,
                                    mcmc.event_extents,
                                    mcmc.event_observations,
                                    proposed_events)
        end
        proposed_events_array = reshape([proposed_events_array[1:(aug_order[j]-1)]
                                         _time(proposed_event)
                                         proposed_events_array[(aug_order[j]+1):end]],
                                        size(new_events_array))
        proposed_events = Events{T}(proposed_events_array)
      end
      proposed_rp = new_rp
      proposed_tr_cache = new_tr_cache
    end
    proposed_lprior = logprior(proposed_rp, mcmc.risk_priors)
    # Based on the logprior and competiting MCMC iteration, this loglikelihood is required for acceptance
    # Calculating this in advance allows us to cut loglikelihood calculation short if it goes below threshold
    ll_acceptance_threshold = log(rand()) + new_lposterior - proposed_lprior
    if ll_acceptance_threshold < Inf
      proposed_llikelihood = loglikelihood(proposed_rp,
                                           mcmc.risk_functions,
                                           proposed_events,
                                           mcmc.population,
                                           mcmc.starting_states,
                                           transmission_rates_output = false,
                                           transmissions_output = false,
                                           early_decision_value = ll_acceptance_threshold,
                                           transmission_rate_cache = proposed_tr_cache)[1]
      proposed_lposterior = proposed_lprior + proposed_llikelihood
    else
      proposed_llikelihood = -Inf
      proposed_lposterior = -Inf
    end
    if proposed_llikelihood >= ll_acceptance_threshold
      @debug "MCMC proposal accepted (acceptance probability = $(round(min(1.0, exp(proposed_lposterior - new_lposterior)), digits=3)))"
      new_rp = proposed_rp
      new_tr_cache = proposed_tr_cache
      new_events = proposed_events
      new_events_array = proposed_events_array
      new_lposterior = proposed_lposterior
    else
      @debug "next!: MCMC proposal rejected (acceptance probability = $(round(min(1.0, exp(proposed_lposterior - new_lposterior)), digits=3)))"
    end
  end
  mc.iterations += 1
  tnr, tx = loglikelihood(new_rp,
                          mcmc.risk_functions,
                          new_events,
                          mcmc.population,
                          mcmc.starting_states,
                          loglikelihood_output = false,
                          transmission_rate_cache = new_tr_cache)[[2; 3]]
  new_network = generate(TransmissionNetwork, tnr, mcmc, tx)
  push!(mc.events, new_events)
  push!(mc.transmission_network, new_network)
  push!(mc.risk_parameters, new_rp)
  push!(mc.log_posterior, new_lposterior)
  return mc
end

function update!(mc::MarkovChain{T, M},
                 mcmc::MCMC{T, M},
                 Σrp::Array{Float64, 2},
                 Σsm::Array{Float64, 2},
                 σ::Float64;
                 event_batches::Int64=1,
                 transmission_rate_cache::Union{Nothing, TransmissionRateCache}=nothing,
                 DA_couple::Bool=false) where {
                 T <: DiseaseStateSequence,
                 M <: PhyloILM}
  # Initialize
  local tnr
  new_tr_cache = transmission_rate_cache
  new_events = mc.events[end]
  new_events_array = convert(Array{Float64, 2}, new_events)
  new_rp = mc.risk_parameters[end]
  new_sm = mc.substitution_model[end]
  new_network = mc.transmission_network[end]
  new_lposterior = mc.log_posterior[end]
  # Randomize event time augmentation
  event_indices = findall(new_events_array[:] .> -Inf)
  aug_order = sample(event_indices, length(event_indices), replace=false)
  @debug "$(length(aug_order)) events to augment"
  if event_batches < 0
    @error "Cannot have negative amount of event batches"
  end
  if event_batches > length(aug_order)
    @warn "More event batches than there are events to augment ($(event_batches) > $(length(aug_order))), setting to maximum ($(length(aug_order)))"
    event_batches = length(aug_order)
  end
  batch_size = fld(length(aug_order), event_batches)
  @debug "Performing data augmentation in batches of $batch_size events at a time"
  steps = DA_couple ? 1 + event_batches : 1 + event_batches + 1
  for i = 1:steps
    if i == 1
      # Propose new risk parameters
      proposed_events = new_events
      proposed_events_array = new_events_array
      proposed_rp = generate(RiskParameters{T}, new_rp, Σrp)
      proposed_sm = generate(mcmc.substitution_model, new_sm, Σsm)
      @debug "Generating new TransmissionRateCache for TN-ILM parameter proposal"
      proposed_tr_cache = TransmissionRateCache(size(new_events_array, 1))
    elseif i < (event_batches + 2)
      proposed_events_array = copy(new_events_array)
      proposed_events = Events{T}(proposed_events_array)
      t_ids = Int64[]
      for j = (batch_size*(i-2) + 1):min(batch_size*(i-1),length(aug_order))
        id, state_index = Tuple(CartesianIndices((individuals(new_events),
                                            length(convert(DiseaseStates, T)[2:end])))[aug_order[j]])
        # Generate new transmission source for transmission events
        state_index == 1 && push!(t_ids, id)
        new_state = convert(DiseaseStates, T)[state_index+1]
        time = new_events[new_state][id]
        proposed_event = generate(Event,
                                  Event{T}(time, id, new_state),
                                  σ,
                                  mcmc.event_extents,
                                  mcmc.event_observations,
                                  proposed_events,
                                  new_network)
        proposed_events_array = reshape([proposed_events_array[1:(aug_order[j]-1)]
                                         _time(proposed_event)
                                         proposed_events_array[(aug_order[j]+1):end]],
                                        size(new_events_array))
        proposed_events = Events{T}(proposed_events_array)
      end
      proposed_rp = new_rp
      proposed_sm = new_sm
      proposed_tr_cache = new_tr_cache
    elseif i == (event_batches + 2)
      proposed_events = new_events
      proposed_events_array = new_events_array
      proposed_rp = new_rp
      proposed_sm = new_sm
      proposed_tr_cache = new_tr_cache
    end
    proposed_lprior = logprior(proposed_rp, mcmc.risk_priors)
    proposed_lprior += logprior(proposed_sm, mcmc.substitution_model_prior)
    # Based on the logprior and competiting MCMC iteration, this loglikelihood is required for acceptance
    # Calculating this in advance allows us to cut loglikelihood calculation short if it goes below threshold
    ll_acceptance_threshold = log(rand()) + new_lposterior - proposed_lprior
    if ll_acceptance_threshold < Inf
      proposed_llikelihood, tnr, tx = loglikelihood(proposed_rp,
                                                    mcmc.risk_functions,
                                                    proposed_events,
                                                    mcmc.population,
                                                    mcmc.starting_states,
                                                    early_decision_value = ll_acceptance_threshold,
                                                    transmissions_output = true,
                                                    transmission_rate_cache = proposed_tr_cache)
      if proposed_llikelihood > ll_acceptance_threshold
        if DA_couple
          if i == 1 || length(t_ids) == 0
            proposed_network = new_network
          elseif 1 < i < event_batches + 2
            proposed_network = generate(TransmissionNetwork, new_network, tnr, mcmc, t_ids)
            @debug "Generating a Transmission Network proposal" ProposedNetwork = proposed_network
          end
        else
          if i < event_batches + 2
            proposed_network = new_network
          elseif i == event_batches + 2
            proposed_network = generate(TransmissionNetwork, new_network, tnr, mcmc, sample(tx))
            @debug "Generating a transmission network proposal" ProposedNetwork = proposed_network
          end
        end
        @debug "Generating a Tree for proposed network and event times"
        tree, obs_nodes = generate(Tree, proposed_events, mcmc.event_observations, proposed_network)
        leaf_data = Dict{Int64, RNASeq}()
        for k = eachindex(obs_nodes)
          if !isnothing(obs_nodes[k])
            leaf_data[obs_nodes[k]] = mcmc.event_observations.seq[k]
          end
        end
        proposed_llikelihood += PhyloModels.loglikelihood(tree, proposed_sm, leaf_data)
        proposed_lposterior = proposed_lprior + proposed_llikelihood
      end
    else
      proposed_llikelihood = -Inf
      proposed_lposterior = -Inf
    end
    if proposed_llikelihood >= ll_acceptance_threshold
      @debug "MCMC proposal accepted (acceptance probability = $(round(min(1.0, exp(proposed_lposterior - new_lposterior)), digits=3)))"
      new_rp = proposed_rp
      new_sm = proposed_sm
      new_tr_cache = proposed_tr_cache
      new_events = proposed_events
      new_events_array = proposed_events_array
      new_network = proposed_network
      new_lposterior = proposed_lposterior
    else
      @debug "MCMC proposal rejected (acceptance probability = $(round(min(1.0, exp(proposed_lposterior - new_lposterior)), digits=3)))"
    end
  end
  mc.iterations += 1
  push!(mc.events, new_events)
  push!(mc.transmission_network, new_network)
  push!(mc.risk_parameters, new_rp)
  push!(mc.substitution_model, new_sm)
  push!(mc.log_posterior, new_lposterior)
  return mc
end

ext: jl


text:
println("------------------------------------")
println("mps_test.jl")
println("------------------------------------")

time_initial = time()

# using Revise
# using TensorLabXD
# using MPSLabXD
# using KrylovKit
# V_virtual = SU₂Space(0=>5,1/2=>5,1=>5);
# V_physical = SU₂Space(1/2=>1);
# V_boundary = SU₂Space(0=>1);
# A1 = TensorMap(randn, V_boundary*V_physical, V_virtual);
# A2 = TensorMap(randn,V_virtual*V_physical,V_virtual);
# A3 = TensorMap(randn,V_virtual*V_physical,V_boundary);
# L = 6;
# s = [[1.0] for i in 1:L];

@timedtestset "finite mps" begin
    # mps1 = FiniteMPS{SU₂Space}([A1,A2,A2,A2,A2,A3],s,[:non_canonical]);
    # @test field(mps1) == ℂ
    # @test spacetype(mps1) == SU₂Space == Rep[SU₂]
    # @test sectortype(mps1) == SU2Irrep == Irrep[SU₂]
    # make_left_canonical(mps1);
    # for i in 1:L
    #     @test mps1.Ts[i]'*mps1.Ts[i] ≈ isomorphism(domain(mps1.Ts[i]),domain(mps1.Ts[i]))
    # end
    # @test norm(mps1) ≈ 1.0
    # mps2 = FiniteMPS{SU₂Space}([A1,A2,A2,A2,A2,A2,A2,A3],s,[:non_canonical]);
    # make_right_canonical(mps2);
    # for i in 1:L
    #     T_tmp = TensorMap(undef,codomain(mps2.Ts[i],1),codomain(mps2.Ts[i],1))
    #     @tensor T_tmp[a;d] = mps2.Ts[i][a,b;c]*conj(mps2.Ts[i][d,b;c])
    #     @test T_tmp ≈ isomorphism(codomain(mps2.Ts[i],1),codomain(mps2.Ts[i],1))
    # end
    # @test norm(mps2) ≈ 1.0
end

@timedtestset "infinite mps" begin

end
time_final = time()
printstyled("Finished sector tests in ",
            string(round(time_final-time_initial; sigdigits=3)),
            " seconds."; bold = true, color = Base.info_color())
println()


text:
# This file includes portions from JuliaMath/Combinatorics.jl in modified form.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module CombinationsIterator

using IterTools

export ModuleCombinationsIterator
export Combinations

"""
* Combinations
"""
const ModuleCombinationsIterator = ""

#The combinations iterator
struct combinations
    n::Int
    t::Int
end

function Base.iterate(c::combinations, s=[min(c.t - 1, i) for i in 1:c.t])
    if c.t == 0 # special case to generate 1 result for t==0
        isempty(s) && return (s, [1])
        return
    end
    for i in c.t:-1:1
        s[i] += 1
        if s[i] > (c.n - (c.t - i))
            continue
        end
        for j in i + 1:c.t
            s[j] = s[j - 1] + 1
        end
        break
    end
    s[1] > c.n - c.t + 1 && return
    (s, s)
end

Base.length(c::combinations) = binomial(c.n, c.t)
Base.eltype(::Type{combinations}) = Vector{Int}

"""
Generate all Combinations of ``n`` elements from an indexable object ``a``. Because the number of Combinations can be very large, this function returns an iterator object.
 Use collect(Combinations(a, n)) to get an array of all Combinations.
"""
function Combinations(a, t::Integer)
    if t < 0
        # generate 0 Combinations for negative argument
        t = length(a) + 1
    end
    reorder(c) = [a[ci] for ci in c]
    (reorder(c) for c in combinations(length(a), t))
end

"""
Generate Combinations of the elements of ``a`` of all orders. Chaining of order iterators is eager, but the sequence at each order is lazy.
"""
Combinations(a) = Iterators.flatten([Combinations(a, k) for k = 1:length(a)])

#START-TEST-########################################################

using Test, SeqUtils

function test()
    C = collect(Combinations([2, 3, 5]))
    @testset "Combinations" begin
        @test C == [[2], [3], [5], [2, 3], [2, 5], [3, 5], [2, 3, 5]]
    end
end

function demo()
    Println([Combinations("abcd", 3)...])

    D = [2, 3, 4, 5]
    for c in Combinations(D)
        println(c)
    end

    println()
    C = collect(Combinations([2, 3, 5]))
    Println(C)
end

function perf() end

function main()
    test()
    demo()
    perf()
end

main()

end # module

#=

Array{Char,1}[['a', 'b', 'c'], ['a', 'b', 'd'], ['a', 'c', 'd'], ['b', 'c', 'd']]

[2][3]
[4]
[5]
[2, 3]
[2, 4]
[2, 5]
[3, 4]
[3, 5]
[4, 5]
[2, 3, 4]
[2, 3, 5]
[2, 4, 5]
[3, 4, 5]
[2, 3, 4, 5]

=#


text:
# an array that returns the indices from getindex

import Base: size, getindex, setindex!

"""
  Special array type such that `A[i] = i` in 1D and `A[i, j] = i` in 2D.
  This array can be stack-allocated.  It supportes the standard Array
  constructors.
"""
struct IdentityArray{T, N} <: DenseArray{T, N}
  dims::NTuple{N, Int}
  offset::Int  # compatability with ArrayViews

  function IdentityArray{T, N}(dims::NTuple{N, Int}, offset::Integer) where {T, N}
    return new(dims, offset)
  end

  function IdentityArray{T, N}(dims::NTuple{N, Int}) where {T, N}
    return new(dims, 0)
  end
end

size(A::IdentityArray) = A.dims

function setindex!(A::IdentityArray, v, i::Int)
  throw(ErrorException("setindex! not defind for IdentityArray"))
end

# construct from individual dimensions
# These constructors can be called as IdentityArray{T}(i, j) (which is not
# documented in Julia).
function IdentityArray{T}(i::Integer) where {T}
  return IdentityArray{T, 1}( (Int(i), ))
end

function IdentityArray{T}(i::Integer, j::Integer) where {T}
  return IdentityArray{T, 2}( (Int(i), Int(j)))
end

function IdentityArray{T}(i::Integer, j::Integer, k::Integer) where {T}
  return IdentityArray{T, 3}( (Int(i), Int(j), Int(k)) )
end

function IdentityArray{T}(i::Integer, j::Integer, k::Integer, l::Integer
                         ) where {T}
  return IdentityArray{T, 4}( (Int(i), Int(j), Int(k), Int(l)) )
end

function IdentityArray{T}(i::Integer, j::Integer, k::Integer, l::Integer,
                          m::Integer) where {T}
  return IdentityArray{T, 5}( (Int(i), Int(j), Int(k), Int(l), Int(m)) )
end

function IdentityArray{T}(i::Integer, j::Integer, k::Integer, l::Integer,
                          m::Integer, n::Integer) where {T}
  return IdentityArray{T, 6}( (Int(i), Int(j), Int(k), Int(l), Int(m), Int(n)) )
end

function IdentityArray{T}(i::Integer, j::Integer, k::Integer, l::Integer,
                          m::Integer, n::Integer, o::Integer) where {T}
  return IdentityArray{T, 7}( (Int(i), Int(j), Int(k), Int(l), Int(m), Int(n), Int(o)) )
end

# NTuple constructor
function IdentityArray{T}(d::NTuple{T, Int})
  return IdentityArray{T, N}(d)
end

#------------------------------------------------------------------------------
# Indexing

@inline function getindex(A::IdentityArray{T, 1}, i::Int) where {T}
  return T(i + A.offset)
end

@inline function getindex(A::IdentityArray{T, 2}, i::Integer, j::Integer) where {T}
  return T(i + A.offset) # for 2D arrays, make all columns the same
end

# it looks like Julia automatically implements linar indexing using the 2D method

# not doing the higher dimension yet

#------------------------------------------------------------------------------
# compatability with sview (limited to contiguous views, 2D parent array)

import ODLCommonTools.sview

function sview(A::IdentityArray{T, 2}, rng::UnitRange, i::Integer) where {T}
  dims = (Int(length(rng)), )
  offset = rng[1] - 1
  return IdentityArray{T, 1}(dims, offset)
end

function sview(A::IdentityArray{T, 2}, idx::Colon, i::Integer) where {T}

  dims = (Int(A.dims[1]),)
  return IdentityArray{T, 1}(dims)
end

#------------------------------------------------------------------------------
# FullFace
# This is the real reason for implementing IdentityArray

"""
  This special `DenseFace` contains all the nodes of the element in its stencil.
  This is useful for tricking `AssembleInterface` into assembling the entire
  Jacobian block instead of only the part in the stencil of the real `sbpface`.

  This type is stack-allocated.

"""
struct FullFace{T} <: SummationByParts.DenseFace{T}
  stencilsize::Int  # = numNodesPerElement
  perm::IdentityArray{Int, 2}
end

"""
  Outer constructor

  **Inputs**

   * numNodesPerElement
   * dim: dimension (2D or 3D)
"""
function FullFace(numNodesPerElement::Integer, dim::Integer)

  perm = IdentityArray{Int}(numNodesPerElement, dim + 1)
  return FullFace{Float64}(numNodesPerElement, perm)
end


text:
using Gen
using Gen: get_child
import LinearAlgebra
import CSV
import Random

using Statistics: mean

# TODO: figure out the correction for the probablity of picking the given node
# in the forward and backward proposal

#########################
# load airline data set #
#########################

function get_airline_dataset()
    df = CSV.read("$(@__DIR__)/airline.csv")
    xs = df[1]
    ys = df[2]
    xs .-= minimum(xs) # set x minimum to 0.
    xs /= maximum(xs) # scale x so that maximum is at 1.
    ys .-= mean(ys) # set y mean to 0.
    ys *= 4 / (maximum(ys) - minimum(ys)) # make it fit in the window [-2, 2]
    return (xs, ys)
end

################################
# abstract covariance function #
################################

"""
Node in a tree representing a covariance function
"""
abstract type Node end
abstract type LeafNode <: Node end
abstract type BinaryOpNode <: Node end

"""
Constant kernel
"""
struct Constant <: LeafNode
    param::Float64
end

eval_cov(node::Constant, x1, x2) = node.param

function eval_cov_mat(node::Constant, xs::Vector{Float64})
    n = length(xs)
    fill(node.param, (n, n))
end

"""
Linear kernel
"""
struct Linear <: LeafNode
    param::Float64
end

eval_cov(node::Linear, x1, x2) = (x1 - node.param) * (x2 - node.param)

function eval_cov_mat(node::Linear, xs::Vector{Float64})
    xs_minus_param = xs .- node.param
    xs_minus_param * xs_minus_param'
end

"""
Squared exponential kernel
"""
struct SquaredExponential <: LeafNode
    length_scale::Float64
end

function eval_cov(node::SquaredExponential, x1, x2)
    exp(-0.5 * (x1 - x2) * (x1 - x2) / node.length_scale)
end

function eval_cov_mat(node::SquaredExponential, xs::Vector{Float64})
    diff = xs .- xs'
    exp.(-0.5 .* diff .* diff ./ node.length_scale)
end

"""
Periodic kernel
"""
struct Periodic <: LeafNode
    scale::Float64
    period::Float64
end

function eval_cov(node::Periodic, x1, x2)
    freq = 2 * pi / node.period
    exp((-1/node.scale) * (sin(freq * abs(x1 - x2)))^2)
end

function eval_cov_mat(node::Periodic, xs::Vector{Float64})
    freq = 2 * pi / node.period
    abs_diff = abs.(xs .- xs')
    exp.((-1/node.scale) .* (sin.(freq .* abs_diff)).^2)
end

"""
Plus node
"""
struct Plus <: BinaryOpNode
    left::Node
    right::Node
end

Plus(left, right) = Plus(left, right)

function eval_cov(node::Plus, x1, x2)
    eval_cov(node.left, x1, x2) + eval_cov(node.right, x1, x2)
end

function eval_cov_mat(node::Plus, xs::Vector{Float64})
    eval_cov_mat(node.left, xs) .+ eval_cov_mat(node.right, xs)
end

"""
Times node
"""
struct Times <: BinaryOpNode
    left::Node
    right::Node
end

Times(left, right) = Times(left, right)

function eval_cov(node::Times, x1, x2)
    eval_cov(node.left, x1, x2) * eval_cov(node.right, x1, x2)
end

function eval_cov_mat(node::Times, xs::Vector{Float64})
    eval_cov_mat(node.left, xs) .* eval_cov_mat(node.right, xs)
end

"""
Compute a covariance matrix by evaluating the function on each pair of inputs
"""
function compute_cov_matrix(covariance_fn::Node, noise, xs)
    n = length(xs)
    cov_matrix = Matrix{Float64}(undef, n, n)
    for i=1:n
        for j=1:n
            cov_matrix[i, j] = eval_cov(covariance_fn, xs[i], xs[j])
        end
        cov_matrix[i, i] += noise
    end
    return cov_matrix
end

"""
Compute a covariance function by recursively computing covariance matrices
"""
function compute_cov_matrix_vectorized(covariance_fn, noise, xs)
    n = length(xs)
    eval_cov_mat(covariance_fn, xs) + Matrix(noise * LinearAlgebra.I, n, n)
end

function compute_log_likelihood(cov_matrix::Matrix{Float64}, ys::Vector{Float64})
    n = length(ys)
    logpdf(mvnormal, ys, zeros(n), cov_matrix)
end

"""
Return predictive log likelihood on new input values
"""
function predictive_ll(covariance_fn::Node, noise::Float64,
                       xs::Vector{Float64}, ys::Vector{Float64},
                       new_xs::Vector{Float64}, new_ys::Vector{Float64})
    (conditional_mu, conditional_cov_matrix) = compute_predictive(
        covariance_fn, noise, xs, ys, new_xs)
    logpdf(mvnormal, new_ys, conditional_mu, conditional_cov_matrix)
end

function compute_predictive(covariance_fn::Node, noise::Float64,
                            xs::Vector{Float64}, ys::Vector{Float64},
                            new_xs::Vector{Float64})
    n_prev = length(xs)
    n_new = length(new_xs)
    means = zeros(n_prev + n_new)
    cov_matrix = compute_cov_matrix(covariance_fn, noise, vcat(xs, new_xs))
    cov_matrix_11 = cov_matrix[1:n_prev, 1:n_prev]
    cov_matrix_22 = cov_matrix[n_prev+1:n_prev+n_new, n_prev+1:n_prev+n_new]
    cov_matrix_12 = cov_matrix[1:n_prev, n_prev+1:n_prev+n_new]
    cov_matrix_21 = cov_matrix[n_prev+1:n_prev+n_new, 1:n_prev]
    @assert cov_matrix_12 == cov_matrix_21'
    mu1 = means[1:n_prev]
    mu2 = means[n_prev+1:n_prev+n_new]
    conditional_mu = mu2 + cov_matrix_21 * (cov_matrix_11 \ (ys - mu1))
    conditional_cov_matrix = cov_matrix_22 - cov_matrix_21 * (cov_matrix_11 \ cov_matrix_12)
    conditional_cov_matrix = 0.5 * conditional_cov_matrix + 0.5 * conditional_cov_matrix'
    (conditional_mu, conditional_cov_matrix)
end

"""
Predict output values for some new input values
"""
function predict_ys(covariance_fn::Node, noise::Float64,
                    xs::Vector{Float64}, ys::Vector{Float64},
                    new_xs::Vector{Float64})
    (conditional_mu, conditional_cov_matrix) = compute_predictive(
        covariance_fn, noise, xs, ys, new_xs)
    mvnormal(conditional_mu, conditional_cov_matrix)
end

function compute_mse(covariance_fn, noise, xs_train, ys_train, xs_test, ys_test)
    (conditional_mu, _) = compute_predictive(covariance_fn, noise, xs_train, ys_train, xs_test)
    mean(sum((conditional_mu .- ys_test) .^ 2))
end

const CONSTANT = 1 # 0.2
const LINEAR = 2 # 0.2
const SQUARED_EXP = 3 # 0.2
const PERIODIC = 4 # 0.2
const PLUS = 5 # binary 0.1
const TIMES = 6 # binary 0.1
const BINARY_OPS = Set{Int}([PLUS, TIMES])

const node_dist = Float64[0.2, 0.2, 0.2, 0.2, 0.1, 0.1]

const node_type_to_num_children = Dict(
    CONSTANT => 0,
    LINEAR => 0,
    SQUARED_EXP => 0,
    PERIODIC => 0,
    PLUS => 2,
    TIMES => 2)

const max_branch = 2

@gen function pick_random_node(node::Node, cur::Int, depth::Int)
    if isa(node, LeafNode)
        (cur, depth)
    elseif @trace(bernoulli(0.5), :done => depth)
        (cur, depth)
    elseif @trace(bernoulli(0.5), :recurse_left => cur)
        @trace(pick_random_node(node.left, get_child(cur, 1, max_branch), depth+1))
    else
        @trace(pick_random_node(node.right, get_child(cur, 2, max_branch), depth+1))
    end
end

max_stars_repo_name: dendisuhubdy/Gen

max_issues_repo_name: kumsh/Gen

max_forks_repo_name: kumsh/Gen


text:
#using VoronoiCells
#using GeometricalPredicates

export read_cell_grid, read_n_cell_samples, read_cells, read_one_cellfile

# type of object written to data file
const CELLARRAY = 1
const CELLGRID = 2

type PoissonCellSamples
    _data::Array{Array{VoronoiCell,1},1}
end

Base.length(samples::PoissonCellSamples) = length(samples._data)
Base.getindex(samples::PoissonCellSamples, i) = samples._data[i]

function Base.write(fn::AbstractString, samples::PoissonCellSamples)
    strm = open(fn, "w")
    nsamp = length(samples)
#    write(strm,nsamp)
    for i in 1:nsamp
        cells = samples[1]
        write(strm,cells)
    end
    close(strm)
end

function Base.write(strm::IOStream,cells::Array{VoronoiCell,1})
    ncells = length(cells)
    write(strm,CELLARRAY)
    write(strm,ncells)
    for i in 1:ncells
        cell::VoronoiCell = cells[i]
        gx = getx(cell._generator)
        gy = gety(cell._generator)
        write(strm,gx)
        write(strm,gy)
        nverts = length(cell._verts)
        write(strm,nverts)
        @inbounds for j in 1:nverts
            vert = cell._verts[j]
            write(strm, getx(vert))
            write(strm, gety(vert))
        end
    end
end

function Base.write(strm::IOStream,gcells::VoronoiCellsA)
    write(strm,CELLGRID)
    write(strm,gcells._ngrid)
    write(strm,gcells._scale)
    write(strm,gcells._shift)
    write(strm,gcells._areascale)
    write(strm, gcells._cells)
end

function read_cell_grid(strm::IOStream)
    buf = zeros(Int,1)
    dbuf = zeros(Float64,3)
    read!(strm,buf)
    ngrid = buf[1]
    read!(strm,dbuf)
    scale = dbuf[1]
    shift = dbuf[2]
    areascale = dbuf[3]
    read!(strm,buf)
    data_type = buf[1]
    data_type == CELLARRAY || error("Unknown data type $data_type when reading cell grid")
    cells = read_cells(strm)
    gcells = cellstogrid(cells,ngrid, scale, shift, areascale)
    gcells
end

function Base.write(fn::AbstractString, cells::Array{VoronoiCell,1})
    samples = PoissonCellSamples(Array(Array{VoronoiCell,1},0))
    push!(samples._data,cells)
#    write_cells(fn,samples)
    write(fn,samples)
end

function Base.write(fn::AbstractString, gcells::VoronoiCellsA)
    strm = open(fn, "w")
    write(strm,gcells)
    close(strm)
end

function read_n_cell_samples(strm::IOStream)
    buf = zeros(Int,1)
    read!(strm,buf)
    nsamp = buf[1]
end

function read_cells(strm::IOStream)
    buf = zeros(Int,1)
    read!(strm,buf)
    ncells = buf[1]
    cells = Array(VoronoiCell,0)
    dbuf = zeros(Float64,1)
    for i in 1:ncells
        read!(strm,dbuf)
        gx = dbuf[1]
        read!(strm,dbuf)
        gy = dbuf[1]
        verts = Array(Point2D,0)
        read!(strm,buf)
        nverts = buf[1]
        for j in 1:nverts
            read!(strm,dbuf)
            vx = dbuf[1]
            read!(strm,dbuf)
            vy = dbuf[1]
            push!(verts, Point2D(vx,vy))
        end
        push!(cells,VoronoiCell(Point2D(gx,gy),verts))
    end
    cells
end

function read_cells_top(strm::IOStream)
    buf = zeros(Int,1)
    read!(strm,buf)
    data_type = buf[1]
    if data_type == CELLARRAY
        return read_cells(strm)
    elseif data_type == CELLGRID
        return read_cell_grid(strm)
    else
        error("read_cells_top: Unrecognized data type ", data_type)
    end
end

function read_one_cellfile(fn::AbstractString)
    strm = open(fn,"r")
#    nsamps = read_n_cell_samples(strm)
#    println("n samples ", nsamps)
#    read_cells(strm)
    read_cells_top(strm)
end

function ==(c1::VoronoiCell, c2::VoronoiCell)
    c1._generator == c2._generator || return false
    n1 = length(c1._verts)
    n2 = length(c1._verts)
    n1 == n2 || return false
    for i in 1:n1
        c1._verts[i] == c2._verts[i] || return false
    end
    true
end

# This can be used to check for equality of written and read cells with grid

function ==(cells1::VoronoiCellsA, cells2::VoronoiCellsA)
    length(cells1) != length(cells2) && return false
    for i in 1:length(cells1)
        cells1[i] != cells2[i] && return false
    end
    getshift(cells1) != getshift(cells2) && return false
    getscale(cells1) != getscale(cells2) && return false
    getareascale(cells1) != getareascale(cells2) && return false
    true
end


text:
# Command: New Session
# Refer to src/01_structs/05_Session.jl


text:
module Pickle

import Serialization
using Serialization: AbstractSerializer

using DataStructures
using SparseArrays

if VERSION < v"1.1"
    isnothing(::Nothing) = true
    isnothing(::Any) = false
end

islittle_endian() = Base.ENDIAN_BOM == 0x04030201

export load, store, loads, stores, AbstractPickle, Pickler

include("./readarg.jl")
include("./writearg.jl")
include("./opcode/opcode.jl")
using .OpCodes

include("./defaults.jl")
include("./mt_table.jl")

import Base: push!, pop!, first, isempty, show, setindex!, getindex, length

struct PickleStack
  meta::Stack{Stack{Any}}
end

function PickleStack()
  meta = Stack{Stack{Any}}()
  s = PickleStack(meta)
  mark!(s)
	return s
end

getstack(s::PickleStack) = first(s.meta)
@inline isempty(s::PickleStack) = isempty(getstack(s))
@inline push!(s::PickleStack, @nospecialize(x)) = push!(getstack(s), x)
@inline pop!(s::PickleStack) = pop!(getstack(s))
@inline pop2!(s::PickleStack) = (pop!(s), pop!(s))
@inline pop3!(s::PickleStack) = (pop!(s), pop!(s), pop!(s))
@inline first(s::PickleStack) = first(getstack(s))
@inline mark!(s::PickleStack) = push!(s.meta, Stack{Any}())
@inline unmark!(s::PickleStack) = reverse!(collect(pop!(s.meta)))
@inline updatefirst!(s::PickleStack, @nospecialize(x)) = (pop!(s); push!(s, x))
@inline maybeupdatefirst!(s::PickleStack, @nospecialize(x)) =
  objectid(first(s)) == objectid(x) ? nothing : updatefirst!(s, x)

struct Memo
  data::Dict{Int, Any}
  ref::IdDict{Any, Int}
end

function Memo(data)
  ref = IdDict{Any, Int}()
  for (k, v) in data
    ref[v] = k
  end
  return Memo(data, ref)
end

function setindex!(m::Memo, @nospecialize(value), key)
  setindex!(m.ref, key, value)
  setindex!(m.data, value, key)
end

length(m::Memo) = length(m.data)
getindex(m::Memo, key) = getindex(m.data, key)
hasref(m::Memo, @nospecialize(key)) = haskey(m.ref, key)
getref(m::Memo, @nospecialize(key)) = getindex(m.ref, key)

@inline function maybeupdate!(m::Memo, key, value)
  @nospecialize key value
  objectid(key) == objectid(value) ? nothing : update!(m, key, value)
end

@inline function update!(m::Memo, key, value)
  @nospecialize key value
  mid = m.ref[key]
  objectid(key) != objectid(value) && delete!(m.ref, key)
  setindex!(m, value, mid)
end

abstract type AbstractPickle <: AbstractSerializer end

"default pickle protocol version"
const DEFAULT_PROTO = 4

struct Pickler{PROTO} <: AbstractPickle
  memo::Memo
  stack::PickleStack
  mt::HierarchicalTable
end

function Pickler(proto=DEFAULT_PROTO, memo=Dict())
  mt = HierarchicalTable()
  Pickler{proto}(Memo(memo), PickleStack(), mt)
end

protocol(::Pickler{P}) where {P} = P
isbinary(pklr::Pickler) = protocol(pklr) >= 1

include("./deserializer.jl")
include("./serializer.jl")

# numpy hack
include("./np.jl")
include("sparse.jl")

include("./torch/torch.jl")
using .Torch

end # module


text:
using DiffEqBase: DEProblem, ODEProblem, DiscreteProblem, remake,
    parameterless_type
using .Stages: Stageable, AbstractSource, AbstractStage,
    StageIterator, StageState, StagedSolver, goto!
import .Stages: record!, current_result, is_finished

function time_type(is_discrete; times...)
    TT = promote_type((typeof(v) for (_, v) in times)...)

    if is_discrete && ! (TT <: Integer)
        msg = "Non integer time(s):"
        for (name, val) in times
            if ! (val isa Integer)
                msg = "$msg\n  $name::$(typeof(val)) = $val"
            end
        end
        error(msg)
    end

    return TT
end

"""
    LEProblem(phase_prob, t_attr; <keyword arguments>)

# Arguments
- `phase_prob`: Phase space dynamics represented in the form of
  `ODEProblem` or `DiscreteProblem` from DifferentialEquations.jl.
  `phase_prob.tspan` is ignored.
- `t_attr::Real`: Simulated time on the (presumed) attractor
  used for the solver.  It is used for computation of Lyapunov
  Exponents.  Roughly `t_attr / t_renorm` instantaneous exponents are
  sampled.
- `t_tran::Real`: Number of iterations to throw away to get rid of the
  effect from the transient dynamics.
- `t_renorm::Real`: Interval between orthonormalizations (renormalizations).
- `dim_lyap::Integer`: Number of Lyapunov exponents to be calculated.
  Default to the full system dimension.
- `Q0::Array`: The initial guess of the Gram-Schmidt "Lyapunov vectors".
- `tangent_dynamics::Function`: A vector field for solving phase space
   evolution *and* tangent space evolution together.  If this is not
   provided, `tangent_dynamics` is derived from `phase_prob.f`.  See
   also [`PhaseTangentDynamics`](@ref).
"""
struct LEProblem{PPr, TPr, TT} <: AbstractSource
    phase_prob::PPr
    tangent_prob::TPr
    t_tran::TT
    t_attr::TT
    t_renorm::TT

    function LEProblem{PPr, TPr}(
            phase_prob::PPr,
            tangent_prob::TPr,
            t_attr::Real;
            t_renorm::Real = 1,
            t_tran::Real = 1,
    ) where {PPr, TPr}

        TT = time_type((PPr <: DiscreteProblem);
                       t_tran = t_tran,
                       t_attr = t_attr,
                       t_renorm = t_renorm)

        new{PPr, TPr, TT}(
            phase_prob, tangent_prob,
            t_tran, t_attr, t_renorm)
    end
end

make_tangent_prob(phase_prob, tangent_dynamics, u0) =
    remake(phase_prob,
           f = tangent_dynamics,
           u0 = u0)

function validate_tangent_prob(phase_prob::DEProblem,
                               tangent_prob,
                               tangent_dynamics,
                               Q0)
    if tangent_prob !== nothing
        if tangent_dynamics !== nothing
            error("Both tangent_prob and tangent_dynamics are given. ",
                  "Only at most one of them can be specified.")
        end
    else
        x0 = phase_prob.u0
        u0 = similar(x0, (length(x0), size(Q0, 2) + 1))
        u0[:, 1] = x0  # ignored anyway
        u0[:, 2:end] = Q0
        if tangent_dynamics == nothing
            tangent_dynamics = fdiff_tangent_dynamics(phase_prob.f, u0)
        end
        tangent_prob = make_tangent_prob(phase_prob, tangent_dynamics, u0)
    end
    if ! is_semi_unitary(@view tangent_prob.u0[:, 2:end])
        error("Columns in Q0 are not orthonormal.\n",
              "Q0 = tangent_prob.u0[:, 2:end]")
    end
    return tangent_prob
end

function LEProblem(phase_prob::PPr;
                   t_attr::Real = error("Keyword argument",
                                        " `t_attr` is required."),
                   tangent_dynamics = nothing,
                   tangent_prob::Union{Void, DEProblem} = nothing,
                   dim_lyap = dimension(phase_prob),
                   Q0 = default_Q0(phase_prob, dimension(phase_prob), dim_lyap),
                   kwargs...) where {PPr <: DEProblem}
    tangent_prob = validate_tangent_prob(
        phase_prob, tangent_prob, tangent_dynamics, Q0)
    LEProblem{parameterless_type(PPr),
              parameterless_type(typeof(tangent_prob))}(
        phase_prob, tangent_prob, t_attr;
        kwargs...)
end

mutable struct PhaseRelaxer{Intr, T} <: AbstractStage
    integrator::Intr
    t_tran::T
    i::Int

    PhaseRelaxer(integrator::Intr, t_tran::T) where {Intr, T} =
        new{Intr, T}(integrator, t_tran, 0)
end

mutable struct LESolution{RecOS, RecFTLE, SS, MS, VV, CH}
    series::SS
    main_stat::MS
    ftle_history::VV
    convergence::CH
    num_orth::Int
    # TODO: Bundle ftle_history and num_orth in a struct.
    converged::Bool
end

const LESolRecOS = LESolution{true}
const LESolRecFTLE = LESolution{RecOS, true} where {RecOS}

function LESolution(dim_lyap::Int;
                    main_stat = VecMean,
                    add_stats = [],
                    history_type = Vector{Float64},
                    num_attr = nothing)
    RecOS = main_stat !== nothing || ! isempty(add_stats)
    RecFTLE = num_attr !== nothing
    if ! isa(main_stat, Union{OnlineStats.OnlineStat, Void})
        main_stat = main_stat(dim_lyap)
    end
    series = if RecOS
        OnlineStats.Series(main_stat, add_stats...)
    end
    ftle_history = if RecFTLE
        [history_type(dim_lyap) for _ in 1:num_attr]
        # TODO: generalize outer array type?
    end
    convergence = if RecOS && RecFTLE
        ConvergenceHistory(dim_lyap)
    end
    args = (series, main_stat, ftle_history, convergence)
    return LESolution{RecOS, RecFTLE, map(typeof, args)...}(
        args...,
        0,                      # num_orth
        false,                  # converged
    )
end

NullLESolution() = LESolution(0; main_stat = nothing)

# TODO: define record!(sol::LESolution, ftle)

"""
    get_dim_lyap(thing)::Int
"""
function get_dim_lyap end

abstract type Terminator end

abstract type AbstractRenormalizer{S} <: AbstractStage end

"""
    TangentRenormalizer(integrator; <keyword arguments>)

A type representing the main calculation of Lyapunov Exponents (LE).
This struct holds all temporary state required for LE calculation.
"""
mutable struct TangentRenormalizer{S <: LESolution,
                                   TMNR <: Terminator,
                                   Intr, T, V, M,
                                   } <: AbstractRenormalizer{S}
    integrator::Intr
    t_attr::T
    t_renorm::T
    inst_exponents::V
    i::Int
    phase_state::V
    tangent_state::M
    Q::M
    R::M
    sign_R::Vector{Bool}
    # TODO: Make sure that they result in concrete types

    sol::S
    tmnr::TMNR
    is_finished::Bool

    function TangentRenormalizer(
            integrator::Intr, t_attr, t_renorm,
            sol::S = NullLESolution(),
            tmnr::TMNR = NullTerminator();
            phase_state::V = init_phase_state(integrator),
            tangent_state::M = init_tangent_state(integrator),
            ) where {S, Intr,
                     TMNR <: Terminator,
                     V <: AbstractVector,
                     M <: AbstractMatrix}

        dim_lyap = get_dim_lyap(integrator)
        dim_phase = length(init_phase_state(integrator))
        @assert dim_lyap <= dim_phase
        @assert size(phase_state) == (dim_phase,)
        @assert size(tangent_state) == (dim_phase, dim_lyap)

        inst_exponents = zeros(eltype(phase_state), dim_lyap)
        Q = similar(tangent_state)
        R = similar(tangent_state, (0, 0))  # dummy
        sign_R = Array{Bool}(dim_lyap)

        T = promote_type(map(typeof, (t_attr, t_renorm))...)
        return new{S, TMNR, Intr, T, V, M}(
            integrator,
            t_attr,
            t_renorm,
            inst_exponents,
            0,  # i
            phase_state,
            tangent_state,
            Q,
            R,
            sign_R,
            sol,
            tmnr,
            false,  # is_finished
        )
    end
end

mutable struct MLERenormalizer{S <: LESolution,
                               TMNR <: Terminator,
                               Intr,
                               T <: Real,
                               E <: Real,
                               V <: AbstractVector,
                               M <: AbstractMatrix,
                               } <: AbstractRenormalizer{S}
    integrator::Intr
    t_attr::T
    t_renorm::T
    exponent::E
    inst_exponent::E
    i::Int
    phase_state::V
    tangent_state::M
    # TODO: Make sure that they result in concrete types

    sol::S
    tmnr::TMNR
    is_finished::Bool

    function MLERenormalizer(
            integrator::Intr, t_attr, t_renorm,
            sol::S = NullLESolution(),
            tmnr::TMNR = NullTerminator();
            phase_state::V = init_phase_state(integrator),
            tangent_state::M = init_tangent_state(integrator),
            ) where {S, Intr,
                     TMNR <: Terminator,
                     V <: AbstractVector,
                     M <: AbstractMatrix}

        if size(tangent_state) != (length(phase_state), 1)
            error("tangent_state must be an array of",
                  " size ($(length(phase_state)), 1);",
                  " given: $(size(tangent_state))")
        end

        E = eltype(phase_state)
        exponent = E(0)
        T = promote_type(map(typeof, (t_attr, t_renorm))...)
        new{S, TMNR, Intr, T, E, V, M}(
            integrator,
            t_attr,
            t_renorm,
            exponent,
            exponent,
            0,  # i
            phase_state,
            tangent_state,
            sol,
            tmnr,
            false,  # is_finished
        )
    end
end

"""
    phase_state(stage) :: Vector

Get current phase-space state stored in `stage`.
"""
phase_state(stage::TangentRenormalizer) = stage.phase_state
phase_state(stage::MLERenormalizer) = stage.phase_state

const LESolver = StagedSolver{<: LEProblem, <: LESolution}
const LESolverRecOS = StagedSolver{<: LEProblem, <: LESolRecOS}
const LESolverRecFTLE = StagedSolver{<: LEProblem, <: LESolRecFTLE}

"""
    LESolver(prob::LEProblem; record::Bool = false)

Create a solver object for a [`LEProblem`](@ref).  Record all
finite-time (instantaneous) Lyapunov exponents when `record = true` is
passed.
"""
function LESolver(prob::LEProblem;
                  phase_relaxer = PhaseRelaxer,
                  renormalizer = if get_dim_lyap(prob) == 1
                      MLERenormalizer
                  else
                      TangentRenormalizer
                  end,
                  kwargs...)
    stage_types = [
        phase_relaxer,
        renormalizer,
    ]
    return LESolver(prob, stage_types; kwargs...)
end

function LESolver(prob::LEProblem, stage_types::AbstractVector;
                  record::Bool = false,
                  integrator_options = [],
                  terminator = false,
                  terminator_options = [],
                  kwargs...)
    sol = LESolution(get_dim_lyap(prob);
                     num_attr = if record
                         ceil(Int, prob.t_attr / prob.t_renorm)
                     else
                         nothing
                     end,
                     kwargs...)
    tmnr = Terminator(terminator, prob, sol; terminator_options...)
    # Additional arguments to each `stage_types`:
    args = (prob, sol, tmnr, integrator_options)
    return StagedSolver(prob, sol, stage_types, args)
end

init(prob::LEProblem; kwargs...) = LESolver(prob; kwargs...)

num_orth(stage::AbstractRenormalizer) = stage.i
num_orth(stage::Stageable) = 0
num_orth(solver::LESolver) = num_orth(solver.state.stage)

Base.show(io::IO, solver::LESolver) =
    print(io,
          "#Orth.: ", num_orth(solver), ", ",
          "LEs: ", lyapunov_exponents(solver))

struct IntegratorError <: Exception
    retcode::Symbol
end
Base.showerror(io::IO, ex::IntegratorError) =
    print(io, "IntegratorError: Integration failed with retcode: $(ex.retcode)")


text:
texdir = @__DIR__()[1:end-2] * "tex/figs/"
jldir = @__DIR__() * "/figs/"

## FIGURE: illustration of the trapezoidal rule
using PyPlot
xkcd()
# http://nbviewer.jupyter.org/gist/cmundi/9041366
fig = figure(figsize=(5,5));
ax = axes()
x = linspace(0.5, 1.0, 100)
y = [0.55, 0.95]
p = plot(
         [y[1],y[1]], [0.0,1-y[1]^2], "k-",
         [y[2],y[2]], [0.0,1-y[2]^2], "k-",
         x,1-x.^2, "b",
         y, 1-y.^2, "r",
          )
xticks([])
yticks([])
ax[:set_xlim]([0.5, 1.0])
ax[:set_xlim]([0.25, 1.2])

annotate("error = area between curves",
         xy=[0.75, 1-0.75^2 - 0.03],
         arrowprops=Dict("arrowstyle"=>"->"),
         xytext=[0.1,0.1])

fig[:canvas][:draw]()
savefig(jldir * "/trapezoidal_rule.pdf")


text:
"""
    pet_to_raw(pet_file_path, sim_config, config_name)

AbstractString, PropDict, AbstractString -> Table

[WIP] Full simulation chain pet->stp->pss->raw
    based on simulated energy depositions contained in the HDF5 file
    found in <pet_file_path> and simulation settings given in <sim_config>.

The output name <config_name> is used to construct filenames for cached
    simulation files (currently the same as the simulation config basename)
"""
function pet_to_raw(sim_config::PropDict)
    # function pet_to_raw(pet_file_path::AbstractString, sim_config::PropDict, config_name::AbstractString)
    ## step 1: stepping information
    stp_table = pet_to_stp(sim_config)

    ## step 2: simulate pulses
    # det_meta = PropDicts.read(PropDict, sim_config.detector_metadata)
    # env = Environment(sim_config)
    # ps_simulator = PSSimulator(sim_config)
    # noise_model = NoiseModel(sim_config)

    # pss_table, pss_truth = stp_to_pss(stp_table, det_meta, env, ps_simulator, noise_model, config_name)
    # launch with sim config updated with stp table as new input file
    pss_table, pss_truth = stp_to_pss(load_config(stp_table, sim_config))

    ## step 3: simulate DAQ
    # elec_chain = ElecChain(sim_config)
    # trigger = Trigger(sim_config)
    # daq = DAQ(sim_config)

    # update sim config
    # sim_config = load_config(pss_table, sim_config)

    # pss_to_raw(pss_table, pss_truth, elec_chain, trigger, daq, noise_model)
    # launch with sim config updated with pss table as new input file
    pss_to_raw(load_config(pss_table, sim_config), pss_truth)
end

# for users
function simulate_raw(pet_input_fullpath::AbstractString, det_meta_fullpath::AbstractString, sim_config_filename::AbstractString)
    pet_to_raw(load_config(pet_input_fullpath, det_meta_fullpath, sim_config_filename))
end

function pet_to_pss(sim_config::PropDict)
    ## step 1: stepping information
    stp_table = pet_to_stp(sim_config)

    # pss_table, pss_truth = stp_to_pss(stp_table, det_meta, env, ps_simulator, noise_model, config_name)
    # launch with sim config updated with stp table as new input file
    pss_table, pss_truth = stp_to_pss(load_config(stp_table, sim_config))

    pss_table, pss_truth
end

function pet_to_pss(pet_input_fullpath::AbstractString, det_meta_fullpath::AbstractString, sim_config_filename::AbstractString)
    pet_to_pss(load_config(pet_input_fullpath, det_meta_fullpath, sim_config_filename))
end


text:
__precompile__()

module OrderedResponse

using StatsFuns
using DataFrames
using StatsModels
using Optim
using LinearAlgebra
using StatsBase
using NLSolversBase
using Distributions: _F1

function num_categories(y::AbstractVector{<:Integer})
    y0, L = extrema(y)
    length(unique(y)) == L - y0 + 1  || throw(error("y must be coded 1:n with all values present"))
    y0 == 1                          || throw(error("y must be coded 1:n with all values present"))
    return L
end

include("distributions.jl")
include("outer-wrapper.jl")
include("likelihood.jl")

# end module
end


text:
using Javis, Animations

function ground(args...)
    background("white")
    sethue("black")
end

function circle_with_color(pos, radius, color)
    sethue(color)
    circle(pos, radius, :fill)
end

# the safest option is to declare the Video first all the time
video = Video(600, 600)

# the points that describe our car
points = [
    Point(-200, 20),
    Point(-170, 20),
    Point(-160, 10),
    Point(-150, 10),
    Point(-140, 20),
    Point(140, 20),
    Point(150, 10),
    Point(160, 10),
    Point(170, 20),
    Point(200, 20),
    Point(180, -20),
    Point(100, -60),
    Point(80, -70),
    Point(60, -90),
    Point(0, -90),
    Point(-90, -40),
    Point(-190, -20),
]
npoints = length(points)

Background(1:220, ground)

# let the points appear one by one
objects = [
    Object(frame_start:200, (args...) -> circle(O, 10, :fill), points[i]) for
    (frame_start, i) in zip(1:2:(2 * npoints), 1:npoints)
]

# easiest to move canvas to draw at origin
act!(objects, Action(1:5, appear(:scale)))
act!(objects, Action(GFrames(100:120), anim_scale(0.5)))
act!(objects, Action(GFrames(190:200), disappear(:scale)))

# generate the bezier path
bezierpath = makebezierpath(points)
bezierpathpoly = bezierpathtopoly(bezierpath)

# let the bezier path appear and disappear in the end
bezier_object =
    Object((2 * npoints + 10):200, (args...) -> drawbezierpath(bezierpath, :stroke))

act!(bezier_object, Action(1:10, appear(:fade)))
act!(bezier_object, Action(GFrames(190:200), disappear(:fade)))

# let a red circle appear and follow the bezier path polygon
red_circle = Object(120:220, (args...) -> circle_with_color(first(points), 10, "red"))
act!(red_circle, Action(1:20, appear(:fade)))
act!(red_circle, Action(21:70, sineio(), follow_path(bezierpathpoly .- first(points))))
act!(red_circle, Action(71:80, disappear(:fade)))

render(video; pathname = joinpath(@__DIR__, "gifs/follow_bezier_path.gif"))


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
    OutageOutputs

Struct for storing `bau_critical_load_met` and `bau_critical_load_met_time_steps`.
This struct is attached to the BAUScenario s.t. the outputs can be passed back to the user.
It is mutable s.t. that it can be adjusted in BAUInputs.
"""
Base.@kwdef mutable struct OutageOutputs
    bau_critical_load_met::Bool = true
    bau_critical_load_met_time_steps::Int = 0
end

struct BAUScenario <: AbstractScenario
    settings::Settings
    site::Site
    pvs::Array{PV, 1}
    wind::Wind
    storage::Storage
    electric_tariff::ElectricTariff
    electric_load::ElectricLoad
    electric_utility::ElectricUtility
    financial::Financial
    generator::Generator
    dhw_load::DomesticHotWaterLoad
    space_heating_load::SpaceHeatingLoad
    existing_boiler::ExistingBoiler
    outage_outputs::OutageOutputs
end

function set_min_max_kw_to_existing(tech::AbstractTech)
    techdict = Dict(fn => getfield(tech, fn) for fn in fieldnames(typeof(tech)))
    techdict[:min_kw] = techdict[:existing_kw]
    techdict[:max_kw] = techdict[:existing_kw]
    eval(Meta.parse(string(typeof(tech)) * "(; $techdict...)"))
end

function bau_site(site::Site)
    Site(;
        latitude=site.latitude,
        longitude=site.longitude,
        land_acres=site.land_acres,
        roof_squarefeet=site.roof_squarefeet,
        min_resil_timesteps=0,
        mg_tech_sizes_equal_grid_sizes=site.mg_tech_sizes_equal_grid_sizes,
        node=site.node,
    )
end

"""
    BAUScenario(s::Scenario)

Constructor for BAUScenario (BAU = Business As Usual) struct.
- sets the PV and Generator max_kw values to the existing_kw values
- sets wind and storage max_kw values to zero
"""
function BAUScenario(s::Scenario)

    # set all PV.max_kw to existing_kw
    pvs = PV[]
    for pv in s.pvs
        if pv.existing_kw > 0
            push!(pvs, set_min_max_kw_to_existing(pv))
        end
    end

    # set Generator.max_kw to existing_kw
    generator = set_min_max_kw_to_existing(s.generator)

    # no existing wind
    wind = Wind(; max_kw=0)

    # no existing storage
    storage = Storage(Dict(:elec => Dict(:max_kw => 0)), s.financial)

    t0, tf = s.electric_utility.outage_start_time_step, s.electric_utility.outage_end_time_step
    #=
    When a deterministic grid outage is modeled we must adjust the BAU critical load profile to keep the problem
    feasible and to get the same ElectricTariff costs in both the optimal and BAU scenarios
    (because the BAU scenario may not have enough existing capacity to meet the critical load and because during an
    outage no grid costs are incurred).
    In the simplest case we set the BAU critical_loads_kw to zero during the outage.
    However, if the BAU scenario has existing Generator and/or PV we calculate how many time steps the critical load can
    be met and make the critical load non-zero for those time steps in order to show the most realistic dispatch results.
    This calculation requires the PV prod_factor_series and so it is done in BAUInputs.
    =#
    elec_load = deepcopy(s.electric_load)
    if tf > t0 && t0 > 0
        elec_load.critical_loads_kw[t0:tf] = zeros(tf-t0+1)  # set crit load to zero
    end
    outage_outputs = OutageOutputs()

    #=
    For random or uncertain outages there is no need to zero out the critical load but we do have to
    set the Site.min_resil_timesteps to zero s.t. the model is not forced to meet any critical load
    in the BAUScenario
    =#
    site = bau_site(s.site)

    return BAUScenario(
        s.settings,
        site,
        pvs,
        wind,
        storage,
        s.electric_tariff,
        elec_load,
        s.electric_utility,
        s.financial,
        generator,
        s.dhw_load,
        s.space_heating_load,
        s.existing_boiler,
        outage_outputs
    )
end


text:
export randSampling

"""
```julia
randSampling(A::Matrix, B=I::Matrix; k=0::Int)
```
Generate a random sampling random Matrix `S`, E(SS')=I.
ASS'B is typically used to approximate AB.

# Examples
```julia
S = randSampling(rand(2,3),rand(3,2),k=2)

3×2 Matrix{Float64}:
 0.0      0.0
 0.0      1.2439
 1.15342  0.0
```
"""
function randSampling(A::Matrix,B=I::Matrix; k=0::Int)
    if k <= 0
        k = floor(Int,size(A)[2]/2)
        print("k was <= $(0) or not provided, it is now set to $(k) by default")
    end

    m = size(A)[2]
    if size(B)[1] != m
        return "The matrices can not be multiplied "
    else
        Lens = [norm(A[:,i])*norm(B[i,:]) for i = 1:m ]
        Total = sum(Lens)
        S = spzeros(m,k)
        for (i,j) in enumerate([sample([1:m...], Weights(Lens) ) for i = 1 : k])
            S[j,i] = 1/sqrt(Lens[j]/Total*k) # normalizing so that E[SS']=I
        end
        return S
    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Libdl, p7zip_jll

@testset "p7zip_jll" begin
    @test isfile(p7zip_jll.p7zip_path)
end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
module Book

export xbook

include("xbook.jl")

end


text:
################################################
#
# Module to test Modia/src/symbolic/StateSelection.jl
#
# Author: Martin Otter, DLR-SR
#
################################################
module TestStateSelectionAlgorithm

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

include(joinpath("..", "..", "src", "symbolic", "StateSelection.jl"))
using .StateSelection

isAlgebraic(v::Int, A::Vector{Int}, Arev::Vector{Int}) = A[v] == 0 && Arev[v] == 0

"""
    checkSortedEquationGraph(eq,Vx,Vderx) - check result of sortedEquationGraph(..)

- eq   : Result of sortedEquationGraph(..)
- Vx   : Required value of eq.Vx
- Vderx: Required value of eq.Vderx
"""
function checkSortedEquationGraph(eq;nx=0, nc=0, Vx=Int[], Vderx=Int[])
    @test length(eq.Vx) + length(eq.Vmue) == nx
    @test length(eq.Er)                 == nx
    @test eq.nc    == nc
    @test eq.Vx    == Vx
    @test eq.Vderx == Vderx

   # Check correct indices
    for i in eachindex(eq.Vx)
        vx = eq.Vx[i]
        if vx != 0
            if eq.Arev[vx] > 0
                vxInt_index_isZero = eq.VxRev[ eq.Arev[vx] ] == 0
                @test vxInt_index_isZero == false
            end
        end
    end

    for i in eachindex(eq.Vderx)
        vderx = eq.Vderx[i]
        if vderx == 0
            if !isAlgebraic(eq.Vx[i], eq.A, eq.Arev)
                vx = eq.Vx[i]
                dervx = eq.A[vx]
                @test (dervx == 0) == false
                if dervx != 0
                    der_vx_index = eq.VxRev[ dervx ]
                    @test (der_vx_index == 0) == false
                end
            end
        end
    end
end

@testset "\nTest StateSelection algorithm" begin

    @testset "... Test two coupled inertias (all/subset/no unknowns can be solved for)" begin

    #= -----------------------------------------
    @model CoupledInertias1 begin
        J1_J = 1
        J1_tau = Float()
        J1_phi = Float(start = 1.0)
        J1_w = Float(start = 0.0)

        J2_J = 2
        J2_tau = Float()
        J2_phi = Float(start = 0.0)
        J2_w = Float(start = 0.0)
        @equations begin
            J1_w = der(J1_phi)
            J1_J * der(J1_w) = J1_tau

            J2_w = der(J2_phi)
            J2_J * der(J2_w) = J2_tau

            J2_phi = J1_phi
            J2_tau + J1_tau = 0
        end
    end

gives:

Assigned equations after index reduction:
   1:              der(J1_phi):  J1_w = der(J1_phi)
   2:                der(J1_w):  J1_J * der(J1_w) = J1_tau
   3:              der(J2_phi):  J2_w = der(J2_phi)
   4:                   J2_tau:  J2_J * der(J2_w) = J2_tau
   5:                         :  J2_phi = J1_phi
   6:                   J1_tau:  J2_tau + J1_tau = 0
   7:                         :  DER( J2_phi = J1_phi )
   8:             der2(J1_phi):  DER( J1_w = der(J1_phi) )
   9:                der(J2_w):  DER( J2_w = der(J2_phi) )
  10:             der2(J2_phi):  DER2( J2_phi = J1_phi )

Unassigned variables:
J1_phi, J1_w, J2_phi, J2_w

Unassigned equations:
J2_phi = J1_phi
DER(J2_phi = J1_phi)

NEW STATE SELECTION

vNames:
    1: J1_tau
    2: J1_phi                        A[2] = 7
    3: J1_w                          A[3] = 8
    4: J2_tau
    5: J2_phi                        A[5] = 9
    6: J2_w                          A[6] = 10
    7: der(J1_phi)                   A[7] = 11
    8: der(J1_w)
    9: der(J2_phi)                   A[9] = 12
   10: der(J2_w)
   11: der2(J1_phi)
   12: der2(J2_phi)

=#
# println("\n\n... Test two coupled inertias (all unknowns can be solved for)")

        Gorigin = [[3,7],[8,1],[6,9],[10,4],[5,2],[4,1],[9,7],[8,11],[10,12],[12,11]]
        BLT = Array{Int64,1}[[1],[3],[8,10,9,4,6,2],[5],[7]]
        assign = [6,0,0,4,0,0,1,2,3,9,8,10]
        Avar = [0,7,8,0,9,10,11,0,12,0,0,0]
        Bequ = [8,0,9,0,7,0,10,0,0,0]
        VNames = ["J1_tau",
          "J1_phi",
          "J1_w",
          "J2_tau",
          "J2_phi",
          "J2_w",
          "der(J1_phi)",
          "der(J1_w)",
          "der(J2_phi)",
          "der(J2_w)",
          "der2(J1_phi)",
          "der2(J2_phi)"]

        # First version where all unknowns can be solved for
        eqInit = StateSelection.SortedEquationGraph(Gorigin, BLT, assign, Avar, Bequ, VNames)
        Gsolvable = eqInit.Gunknowns

        eq = getSortedEquationGraph!(eqInit, Gsolvable)
        # printSortedEquationGraph(eq)
        checkSortedEquationGraph(eq; nx=2, nc=0, Vx=[2,7], Vderx=[0,11])

        # println("\n\n... Test two coupled inertias (only a subset of unknowns can be solved for)")
        eqInit = StateSelection.SortedEquationGraph(Gorigin, BLT, assign, Avar, Bequ, VNames)
        Gsolvable = Any[Any[3,7],Any[1],Any[6,9],Any[4],Any[2,5],Any[1,4],Any[7,9],[8,11],[10,12],[12,11]]
        eq = getSortedEquationGraph!(eqInit, Gsolvable)
        # printSortedEquationGraph(eq)
        checkSortedEquationGraph(eq; nx=2, nc=0, Vx=[5,9], Vderx=[0,12])

        # println("\n\n... Test two coupled inertias (no unknowns can be solved for)")
        Gsolvable = newRaggedIntMatrix(length(Gorigin))
        eq = getSortedEquationGraph(Gorigin, Gsolvable, BLT, assign, Avar, Bequ, VNames)
        # printSortedEquationGraph(eq)
        checkSortedEquationGraph(eq; nx=9, nc=4, Vx=[2, 5, 7, 9, 6, 3, 0, 0], Vderx=[0, 0, 11, 12, 10, 8, 4, 1])

    end

    @testset "... Test simple sliding mass model" begin

#=

  1: 0 = r - n*s
  2: 0 = v - der(r)
  3: 0 = m*der(v) - f - m*g - u
  4: 0 = n*f
  5: 0 = u + (c*s + d*der(s) + sf)*n
  9: 0 = der(sf) + sf - su(t)     # additional equation to test solving for der(sf)

  Variables
   1: s
   2: r
   3: v
   4: f
   5: u
  ------------
   6: der(s)
   7: der2(s)
   8: der(r)
   9: der2(r)
  10: der(v)
  11: sf
  12: der(sf)

   Equations
   1: (1)
   2: (2)
   3: (3)
   4: (4)
   5: (5)
   6: der(1)
   7: der2(1)
   8: der(2)
   9: (9)
=#

        assign = [0,0,0,3,5,0,4,0,7,8,0,9]
        A = [6,8,10,0,0,7,0,9,0,0,12,0]
        B = [6,8,0,0,0,7,0,0,0]
        VNames = ["s", "r", "v", "f", "u", "der(s)", "der2(s)",
          "der(r)", "der2(r)", "der(v)", "sf", "der(sf)"]

        BLT = [[5],
       [9],
       [7,8,3,4]]

        G = [[2,1],
     [3,8],
     [10,4,5],
     [4],
     [5,1,6],
     [8,6],
     [9,7],
     [10,9],
     [12,11]]

        Gsolvable = [[2],
             [3,8],
             [10,4,5],
             [],
             [5],
             [8],
             [9],
             [10,9],
             [12]]

        eqInit = SortedEquationGraph(G, BLT, assign, A, B, VNames)
        eq = getSortedEquationGraph!(eqInit, Gsolvable)
        #printSortedEquationGraph(eq)

        checkSortedEquationGraph(eq; nx=3, nc=0, Vx=[1,6,11], Vderx=[0,7,12])

    end

    @testset "... Test Multi-index DAE (with and without tearing)" begin

    #=
  1: 0 = u1(t) + x1 - x2
  2: 0 = u2(t) + x1 + x2 - x3 + der(x6)
  3: 0 = u3(t) + x1 + der(x3) - x4
  4: 0 = u4(t) + 2*der2(x1) + der2(x2) + der2(x3) + der(x4) + der3(x6)
  5: 0 = u5(t) + 3*der2(x1) + 2*der2(x2) + x5
  6: 0 = u6(t) + 2*x6 + x7
  7: 0 = u7(t) + 3*x6 + 4*x7
  8: 0 = u8(t) + x8 - sin(x8)

  Variables
   1: x1
   2: x2
   3: x3
   4: x4
   5: x5          (5)
   6: x6
   7: x7          (6)
   8: der(x1)
   9: der(x2)
  10: der(x3)
  11: der(x4)     (4)
  12: der(x6)
  13: der2(x1)    der2(1)
  14: der2(x2)    der2(2)
  15: der2(x3)    der(3)
  16: der2(x6)
  17: der3(x6)    der3(7)
  18: x8
  ------------
  19: der(x7)     der(6)
  20: der2(x7)    der2(6)
  21: der3(x7)    der3(6)

   Equations
   1: (1)
   2: (2)
   3: (3)
   4: (4)
   5: (5)
   6: (6)
   7: (7)
   8: (8)
   9: der(1)
  10: der2(1)
  11: der(2)
  12: der2(2)
  13: der(3)
  14: der(6)
  15: der(7)
  16: der2(6)
  17: der2(7)
  18: der3(6)
  19: der3(7)
=#

        assign = [0,0,0,0,5,0,6,0,0,0,4,0,10,12,13,0,19,8,14,16,18]
        A = [8,9,10,11,0,12,19,13,14,15,0,16,0,0,0,17,0,0,20,21,0]
        B = [9,11,13,0,0,14,15,0,10,0,12,0,0,16,17,18,19,0,0]
        VNames = ["x1", "x2", "x3", "x4", "x5", "x6", "x7",
          "der(x1)", "der(x2)", "der(x3)", "der(x4)", "der(x6)",
          "der2(x1)", "der2(x2)", "der2(x3)", "der2(x6)",
          "der3(x6)", "x8", "der(x7)", "der2(x7)", "der3(x7)"]

        BLT = [[8],
       [18,19],
       [10,12,13,4,5]]

        G = [[1,2],
     [1,2,3,12],
     [1,10,4],
     [13,14,15,11,17],
     [13,14,5,18],
     [6,7],
     [6,7],
     [18],
     [8,9],
     [13,14],
     [8,9,10,16],
     [13,14,15,17],
     [8,15,11],
     [12,19],
     [12,19],
     [16,20],
     [16,20],
     [17,21],
     [17,21]]

        # println("\n\n... Test Multi-Index DAE without tearing")
        Gsolvable = newRaggedIntMatrix(length(G))
        eqInit = SortedEquationGraph(G, BLT, assign, A, B, VNames)
        eq = getSortedEquationGraph!(eqInit, Gsolvable)
        # printSortedEquationGraph(eq)

        checkSortedEquationGraph(eq; nx=21, nc=12, Vx=[7, 6, 19, 12, 20, 16, 1, 2, 3, 8, 9, 10, 4, 18, 0],
                                           Vderx=[0, 0, 0, 0, 21, 17, 0, 0, 0, 13, 14, 15, 11, 0, 5])

        # println("\n\n... Test Multi-Index DAE WITH tearing")
        Gsolvable    = copy(G)
        Gsolvable[8] = fill(0, 0)  # "0 = u8(t) + x8 - sin(x8)" cannot be solved for x8
        eqInit = SortedEquationGraph(G, BLT, assign, A, B, VNames)
        eq = getSortedEquationGraph!(eqInit, Gsolvable)
        #printSortedEquationGraph(eq)

        checkSortedEquationGraph(eq; nx=8, nc=4, Vx=[7, 19, 20, 2, 9, 18],
                                         Vderx=[0, 0, 21, 0, 14, 0])

    end
end
end


text:
#! /usr/bin/julia

# Rosetta Code, Fibonacci n-step number sequences

type NFib{T<:Integer}
    n::T
    klim::T
    seeder::Function
end

type FState
    a::Array{BigInt,1}
    adex::Integer
    k::Integer
end

function Base.start{T<:Integer}(nf::NFib{T})
    a = nf.seeder(nf.n)
    adex = 1
    k = 1
    return FState(a, adex, k)
end

function Base.done{T<:Integer}(nf::NFib{T}, fs::FState)
    fs.k > nf.klim
end

function Base.next{T<:Integer}(nf::NFib{T}, fs::FState)
    f = sum(fs.a)
    fs.a[fs.adex] = f
    fs.adex = rem1(fs.adex+1, nf.n)
    fs.k += 1
    return (f, fs)
end

function fib_seeder{T<:Integer}(n::T)
    a = zeros(BigInt, n)
    a[1] = one(BigInt)
    return a
end

function fib{T<:Integer}(n::T, k::T)
    NFib(n, k, fib_seeder)
end

function luc_rc_seeder{T<:Integer}(n::T)
    a = zeros(BigInt, n)
    a[1] = 3
    a[2] = -1
    return a
end

function luc_rc{T<:Integer}(n::T, k::T)
    NFib(n, k, luc_rc_seeder)
end

function luc_seeder{T<:Integer}(n::T)
    a = -ones(BigInt, n)
    a[end] = big(n)
    return a
end

function luc{T<:Integer}(n::T, k::T)
    NFib(n, k, luc_seeder)
end

lo = 2
hi = 10
klim = 16

print("n-step Fibonacci for n = (", lo, ",", hi)
println(") up to k = ", klim, ":")
for i in 2:10
    print(@sprintf("%5d => ", i))
    for j in fib(i, klim)
        print(j, " ")
    end
    println()
end

println()
print("n-step Rosetta Code Lucas for n = (", lo, ",", hi)
println(") up to k = ", klim, ":")
for i in 2:10
    print(@sprintf("%5d => ", i))
    for j in luc_rc(i, klim)
        print(j, " ")
    end
    println()
end

println()
print("n-step MathWorld Lucas for n = (", lo, ",", hi)
println(") up to k = ", klim, ":")
for i in 2:10
    print(@sprintf("%5d => ", i))
    for j in luc(i, klim)
        print(j, " ")
    end
    println()
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
try
    files=__req_loaded_files
catch err
    # const prevents this from being overwritten by subsequent loads
    const global __req_loaded_files = Set{String}()
end

# load file if it hasn't been loaded previously in the session
function req(file::String)
    files = __req_loaded_files
    if !has(files, file)
        add(files, file)
        load(file)
    end
end

# force reload all req'd files next time they're req'd
clearreq() = del_all(__req_loaded_files)

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
using Sockets
include("header.jl")
include("strings.jl")

function play_roulette(p :: Player)
    result = rand(1:36)
    total_bet = 0
    winnings = 0
    print_dict(p, "roulette_0")
    printBalance(p)
    print_dict(p, "roulette_1")
    while true
        print_dict(p, "roulette_2")
        s = readline(p.socket)
        lines = split(s)
        if s == ""
            continue
        elseif s == "d"
            if total_bet > 0
                break
            else
                print_dict(p, "roulette_3")
                p.status = reception
                return
            end
        elseif size(lines,1) != 2
            print_dict(p, "repeat")
            continue
        else
            bet = tryparse(Int64, lines[1])

            if bet == nothing || bet < 0
                print_dict(p, "repeat")
                continue
            elseif bet > (p.balance - total_bet)
                print_dict(p, "roulette_4")
                continue
            elseif bet == 0
                print_dict(p, "roulette_3")
                p.status = reception
                return
            end

            number = tryparse(Int, lines[2])

            if number == nothing
                if lines[2] == "red"
                    if (result < 19 && result % 2 == 1) || (result > 19 && result % 2 == 0)
                        winnings += bet
                    else
                        winnings -= bet
                    end
                    total_bet += bet
                    continue
                elseif lines[2] == "black"
                    if (result < 19 && result % 2 == 0) || (result > 19 && result % 2 == 1)
                        winnings += bet
                    else
                        winnings -= bet
                    end
                    total_bet += bet
                    continue
                elseif lines[2] == "1-12"
                    if result <= 12
                        winnings += 3 * bet
                    else
                        winnings -= bet
                    end
                    total_bet += bet
                    continue
                elseif lines[2] == "13-24"
                    if result > 12 && result <= 24
                        winnings += 3 * bet
                    else
                        winnings -= bet
                    end
                    total_bet += bet
                    continue
                elseif lines[2] == "25-36"
                    if result > 24 && result <= 36
                        winnings += 3 * bet
                    else
                        winnings -= bet
                    end
                    total_bet += bet
                    continue
                else
                    print_dict(p, "repeat")
                    continue
                end
            elseif number >= 1 && number <= 36
                if number == result
                    winnings += 36 * bet
                else
                    winnings -= bet
                end

                total_bet += bet
                continue
            else
                print_dict(p, "repeat")
                continue
            end

        end
    end

    print_dict(p, "roulette_5")
    write(p.socket, "Congratulations the number is.. $result. Your total winnings are: $winnings\n")

    p.balance += winnings
end


text:
module MeshSteward

include("boxes.jl")
include("vselect.jl")
include("eselect.jl")
include("utilities.jl")
include("tetrahedra.jl")
include("triangles.jl")
include("quadrilaterals.jl")
include("lines.jl")
include("io.jl")
include("modification.jl")
include("mesh.jl")

# We can either use/import individual functions from MeshSteward like so:
# ```
# using MeshSteward: attach!
# ```
# or we can bring into our context all exported symbols as
# ```
# using MeshSteward.Exports
# ```
include("Exports.jl")

end # module


text:
padfac(alg::NUFT) = alg.padfac

function padimage(alg::NUFT, img)
    pf = padfac(alg)
    ny,nx = size(img)
    nnx = nextpow(2, pf*nx)
    nny = nextpow(2, pf*ny)
    nsx = nnx÷2-nx÷2
    nsy = nny÷2-ny÷2
    cimg = convert(Matrix{Complex{eltype(img)}}, img)
    return PaddedView(zero(eltype(cimg)), cimg,
                      (1:nnx, 1:nny),
                      (nsx+1:nsx+nx, nsy+1:nsy+ny)
                     )
end

function create_cache(alg::ObservedNUFT, img)
    pimg = padimage(alg.alg, img)
    # make nuft plan
    dx, dy = pixelsizes(img)
    plan = plan_nuft(alg, pimg, dx, dy)
    # get phases and pulse functions
    phases = make_phases(alg, img)

    return create_cache(alg, plan, phases, pimg)
end

function create_cache(alg::NUFT, img)
    pimg = padimage(alg, img)
    return NUFTCache(alg, nothing, nothing, pimg)
end

function update_cache(cache::NUFTCache, img)
    pimg = padimage(img, cache.alg)
    return @set cache.img = pimg
end

function nocachevis(m::ModelImage{M,I,<:NUFTCache}, u, v) where {M,I}
    alg = ObservedNUFT(m.cache.alg, vcat(u', v'))
    cache = create_cache(alg, m.image)
    m = @set m.cache = cache
    return visibilities(m, u, v)
end

function checkuv(uv, u, v)
    @assert u == @view(uv[1,:]) "Specified u don't match uv in cache. Did you pass the correct u,v?"
    @assert v == @view(uv[2,:]) "Specified v don't match uv in cache. Did you pass the correct u,v?"
end

function _visibilities(m::ModelImage{M,I,<:NUFTCache{A}},
                      u::AbstractArray,
                      v::AbstractArray) where {M,I,A<:ObservedNUFT}
    checkuv(m.cache.alg.uv, u, v)
    vis =  m.cache.plan*m.cache.img
    conj.(vis).*m.cache.phases
    #return vis
end

function _frule_vis(m::ModelImage{M,<:IntensityMap{<:ForwardDiff.Dual{T,V,P}},<:NUFTCache{A}}) where {M,T,V,P,A<:ObservedNUFT}
    S = typeof(ForwardDiff.value(first(m.cache.img)))
    p = m.cache.plan
    # Compute the fft
    buffer = similar(m.cache.img, S)
    buffer .= ForwardDiff.value.(m.cache.img)
    xtil = p*buffer
    out = similar(buffer, Complex{ForwardDiff.Dual{T,V,P}})
    # Now take the deriv of nuft
    ndxs = ForwardDiff.npartials(first(m.cache.img))
    dxtils = ntuple(ndxs) do n
        buffer .= ForwardDiff.partials.(m.cache.img, n)
        p * buffer
    end
    out = similar(xtil, Complex{ForwardDiff.Dual{T,V,P}})
    for i in eachindex(out)
        dual = getindex.(dxtils, i)
        prim = xtil[i]
        red = ForwardDiff.Dual{T,V,P}(real(prim), ForwardDiff.Partials(real.(dual)))
        imd = ForwardDiff.Dual{T,V,P}(imag(prim), ForwardDiff.Partials(imag.(dual)))
        out[i] = Complex(red, imd)
    end
    return out
end

function _visibilities(m::ModelImage{M,<:IntensityMap{<:ForwardDiff.Dual{T,V,P}},<:NUFTCache{A}},
    u::AbstractArray,
    v::AbstractArray) where {M,T,V,P,A<:ObservedNUFT}
    checkuv(m.cache.alg.uv, u, v)
    # Now reconstruct everything

    vis = _frule_vis(m)
    conj.(vis).*m.cache.phases
#return vis
end

function _visibilities(m::ModelImage{M,I,<:NUFTCache{A}},
                      u::AbstractArray,
                      v::AbstractArray) where {M,I,A<:NUFT}
    return nocachevis(m, u, v)
end

function amplitudes(m::ModelImage{M,I,C}, u::AbstractArray, v::AbstractArray) where {M,I,C<:NUFTCache}
    vis = visibilities(m, u, v)
    vis .= abs.(v)
    return vis
end

function bispectra(m::ModelImage{M,I,C},
                    u1::AbstractArray,
                    v1::AbstractArray,
                    u2::AbstractArray,
                    v2::AbstractArray,
                    u3::AbstractArray,
                    v3::AbstractArray) where {M,I,C<:NUFTCache}
    # TODO: Fix this so we can actually use a stored cache for closures
    n = length(u1)
    # In testing it is faster to concatenate everything into 1 vector and then
    # take an NFFT of that
    u = vcat(u1,u2,u3)
    v = vcat(v1,v2,v3)
    vis = nocachevis(m, u, v)
    @inbounds (@view(vis[1:n])).*(@view(vis[n+1:2n])).*(@view vis[2n+1:3n])
end

function closure_phases(m::ModelImage{M,I,C},
                        u1::AbstractArray,
                        v1::AbstractArray,
                        u2::AbstractArray,
                        v2::AbstractArray,
                        u3::AbstractArray,
                        v3::AbstractArray) where {M,I,C<:NUFTCache}
    # TODO: Fix this so we can actually use a stored cache for closures
    n = length(u1)
    # In testing it is faster to concatenate everything into 1 vector and then
    # take an NFFT of that
    u = vcat(u1,u2,u3)
    v = vcat(v1,v2,v3)
    vis = nocachevis(m, u, v)
    phase = @inbounds angle.(@view(vis[1:n])) .+ angle.(@view(vis[n+1:2n])) .+ angle.(@view vis[2n+1:3n])
    return phase
end

function logclosure_amplitudes(m::ModelImage{M,I,C},
                               u1::AbstractArray,
                               v1::AbstractArray,
                               u2::AbstractArray,
                               v2::AbstractArray,
                               u3::AbstractArray,
                               v3::AbstractArray,
                               u4::AbstractArray,
                               v4::AbstractArray,
                               ) where {M,I,C<:NUFTCache}
    # TODO: Fix this so we can actually use a stored cache for closures
    n = length(u1)
    u = vcat(u1,u2,u3,u4)
    v = vcat(v1,v2,v3,v4)
    vis = nocachevis(m, u, v)
    amp1 = @view vis[1:n]
    amp2 = @view vis[n+1:2n]
    amp3 = @view vis[2n+1:3n]
    amp4 = @view vis[3n+1:4n]

    lcamp = @. log(abs(amp1)*abs(amp2)/(abs(amp3)*abs(amp4)))
    return lcamp
end

"""
    NFFTAlg
Uses a non-uniform FFT to compute the visibilities.
You can optionally pass uv which are the uv positions you will
compute the NFFT at. This can allow for the NFFT plan to be cached improving
performance

# Fields
$(FIELDS)

"""
Base.@kwdef struct NFFTAlg <: NUFT
    """
    Amount to pad the image
    """
    padfac::Int = 1
    """
    Controls the accuracy of the NFFT usually don't need to change this
    """
    m::Int = 10
end
include(joinpath(@__DIR__, "nfft_alg.jl"))

"""
    DFTAlg
Uses a discrete fourier transform. This is not very efficient for larger images. In those cases
 NFFTAlg or FFTAlg are more reasonable. For small images this is a reasonable choice especially
since it's easy to define derivatives.
"""
struct DFTAlg <: NUFT end
include(joinpath(@__DIR__, "dft_alg.jl"))


text:
using StatsBase
using Distributions
using HypothesisTests

function gen_data!(x, y, d, te)
    rand!(d, x)
    rand!(d, y)
    n = length(x)
    @assert n == length(y)
    for i in 1:n
        y[i] += te
    end
end

threshold(x::Array, p::Real) = quantile(x, p)

# Winsorize above only
function winsorize(x::Array, τ::Real)
    n = length(x)

    x′ = copy(x)
    for i in 1:n
        x′[i] = min(x[i], τ)
    end

    x′
end

# Trim above only
function trim(x::Array, τ::Real)
    n = length(x)

    n_remove = 0
    for i in 1:n
        if x[i] > τ
            n_remove += 1
        end
    end

    n′ = n - n_remove
    x′ = similar(x, n′)

    j = 0
    for i in 1:n
        if x[i] <= τ
            j += 1
            x′[j] = x[i]
        end
    end

    x′
end

function unadjusted(x, y, τ_x, τ_y, τ_xy)
    x′ = copy(x)
    y′ = copy(y)
    x′, y′
end

function trim_separate_thresholds(x, y, τ_x, τ_y, τ_xy)
    x′ = trim(x, τ_x)
    y′ = trim(y, τ_y)
    x′, y′
end

function trim_threshold_from_A(x, y, τ_x, τ_y, τ_xy)
    x′ = trim(x, τ_x)
    y′ = trim(y, τ_x)
    x′, y′
end

function trim_shared_threshold(x, y, τ_x, τ_y, τ_xy)
    x′ = trim(x, τ_xy)
    y′ = trim(y, τ_xy)
    x′, y′
end

function winsorize_separate_thresholds(x, y, τ_x, τ_y, τ_xy)
    x′ = winsorize(x, τ_x)
    y′ = winsorize(y, τ_y)
    x′, y′
end

function winsorize_threshold_from_A(x, y, τ_x, τ_y, τ_xy)
    x′ = winsorize(x, τ_x)
    y′ = winsorize(y, τ_x)
    x′, y′
end

function winsorize_shared_threshold(x, y, τ_x, τ_y, τ_xy)
    x′ = winsorize(x, τ_xy)
    y′ = winsorize(y, τ_xy)
    x′, y′
end

function analyze(x′, y′)
    δ = mean(y′) - mean(x′)
    σ = sqrt(var(y′) / length(y′) + var(x′) / length(x′))
    p = pvalue(UnequalVarianceTTest(y′, x′))
    lower, upper = ci(UnequalVarianceTTest(y′, x′))
    δ, σ, p, lower, upper
end


text:
module ReturnBestOnMetric

using CatBoost

train_data = [[0, 3], [4, 1], [8, 1], [9, 1]]

train_labels = [0, 0, 1, 1]

eval_data = [[2, 1], [3, 1], [9, 0], [5, 3]]

eval_labels = [0, 1, 1, 0]

eval_dataset = Pool(; data=eval_data, label=eval_labels)

model = CatBoostClassifier(; learning_rate=0.03,
                           custom_metric=["Logloss", "AUC:hints=skip_train~false"])
fit!(model, train_data, train_labels; eval_set=eval_dataset, verbose=false)

for (k, v) in model.get_best_score()
    println("$k : $v")
end

end # module


text:
# Test file

using FractalTools
using Plots
pyplot()
# default(:size, (1900, 1000))

# Construct interpolation data
xi, dx, xf = 0, 0.1, 1
x = collect(xi : dx : xf)
y = sin.(2π * x)
z = cos.(2π * x)
u = 0.9
d = u * ones(length(x) - 1)
h = u * ones(length(x) - 1)
l = u * ones(length(x) - 1)
m = u * ones(length(x) - 1)

plt = plot(layout=(5,1))
for niter in 1 : 5
    interp = interpolate(x, y, d, niter=1)
    # interp = hiddenfif(x, y, z, d, h, l, m)
    xd = collect(xi : 0.1dx : xf)
    plot!(xd, getindex.(interp.(xd), 1), subplot=niter)
    scatter!(x, y, subplot=niter)
end
display(plt)


text:
#############################
## The functions in this file compute linear Volterra integral equations
## on the condition that the kernel can be supplied as a vector of polynomial coefficients.
## Its main use outside of constant kernel edge cases is testing.
## In most cases you will want to instead use the functions in src/triVolterraFullKernel.jl
#############################
##  triVolterraEQ1Solver : Solves Volterra integral equation of the first kind, i.e.
##                         g(x)= \int_0^(x) K(x,t) u(t) dt
##                                          OR
##                         g(x) =  \int_0^(1-x) K(x,t) u(t) dt (if flip = false)
##                         where g and the kernel K are given and u(t) is the unknown
##                         to be solved for. K must be monomial coefficient vector.
####
function triVolterraEQ1Solver(g,Kmon::Array{Float64,1},N::Integer,flip::Bool)
    V = triVolterraOpP01(Kmon,N,flip)
    x = Fun(identity,Jacobi(0,1,0..1))
    V = V-Multiplication(x, Jacobi(0,1,0..1))[1:N,1:N]*V
    if flip
        g = Fun(x->g(1-x),Jacobi(0,1, 0..1))
        u = Fun(Jacobi(0,1, 0..1), V\pad(g.coefficients,N))
    else
        u = Fun(Jacobi(0,1, 0..1), V\pad(g.coefficients,N))
    end
    return u
end
#############################
##  triVolterraEQ2Solver : Solves Volterra integral equation of the second kind, i.e.
##                         u(x) = g(x) + \int_0^(x) K(x,t) u(t) dt
##                                          OR
##                         u(x) = g(x) + \int_0^(1-x) K(x,t) u(t) dt (if flip = false)
##                         where g and the kernel K are given and u(x) is the unknown
##                         to be solved for. K must be monomial coefficient vector.
####
function triVolterraEQ2Solver(g,Kmon::Array{Float64,1},N::Integer,flip::Bool)
    V = triVolterraOpP01(Kmon,N,flip)
    x = Fun(identity,Jacobi(0,1,0..1))
    V = V-Multiplication(x, Jacobi(0,1,0..1))[1:N,1:N]*V
    if flip
        C01 = Conversion(Jacobi(0,1,0..1),Jacobi(1,1,0..1))[1:N,1:N]
        u = Fun(Jacobi(0,1, 0..1), (C01-reflectPabtoPba(N)*C01*V)\(C01*pad(g.coefficients,N)))
    else
        u = Fun(Jacobi(0,1, 0..1), (I-V)\pad(g.coefficients,N))
    end
    return u
end
#############################
##   triVolterraFun  :  Returns a Fun which evaluated at x gives the Volterra integral
##                      of the function f given as a Fun with the kernel K(x,y) given in monomial expansion.
##
##                      use "flip = false" for integral from 0 to 1-x
####
function triVolterraFunP01(f,Kmon::Array{Float64,1},N::Integer,flip::Bool)
    if length(Kmon)<2
        Kmon=pad(Kmon,2)
    end
    if flip
        x = Fun(identity,Jacobi(0,1,0..1))
        J = Multiplication(x, Jacobi(0,1,0..1))[1:N,1:N]
        VintFun = Fun(Jacobi(0,1, reverseorientation(0..1)), (I-J)*(triVolterraOpP01(Kmon,N,flip)*pad(f.coefficients,N)))
    else
        VintFun = Fun(JacobiWeight(0,1, Jacobi(0,1, 0..1)), 1/2*triVolterraOpP01(Kmon,N,flip)*pad(f.coefficients,N))
    end
    return VintFun
end
#############################
##   triVolterraOpP01:  Generates Volterra integral operator sum(J^mQEJ^n)
##                      when given monomial expansion of kernel K(x,y)
##                      on JacobiTriangle()
####
function triVolterraOpP01(Kmon::Array{Float64,1},N::Integer,flip::Bool)
    # generating algorithm does not work for k=1 so we pad to k=2
    if length(Kmon)<2
        Kmon=pad(Kmon,2)
    end
    d_K = length(Kmon)
    k = Int(ceil(sqrt(d_K)))
    basis = triJQEJbasislexP01(k,N,flip)
    operator = Kmon[1]*basis[1]
    for j=2:d_K
        operator=operator+Kmon[j]*basis[j]
    end
    return operator
end
#############################
##   triJQEJbasislexP01:  Generates full basis of J^lQEJ^k in a lexicographically sorted vector.
##                        basis[1] will be 00-element, basis[2] will be 01 element and so on.
##                        k>1 is the highest order of x and y in the kernel.
####
function triJQEJbasislexP01(k::Integer,N::Integer,flip::Bool)
    # generating algorithm does not work for k=1 so we pad to k=2
    if k<2
        k=k+1
    end
    S = Jacobi(0,1,0..1)
    x = Fun(identity,S)
    Jy = Multiplication(x, S)[1:N,1:N]
    if flip
        Jx = Multiplication(1-x, S)[1:N,1:N]
    else
        Jx = Jy
    end
    QEy=triQEygenP01(N)
    basis=[]
    for j=1:(k+1)^2
        push!(basis,zeros(N,N))
    end
    # some initialization for convenience
    basis[1]=QEy
    basis[2]=QEy*Jy
    basis[3]=Jx*QEy
    for i=2:k
        # compute y-only
        basis[i^2+1]=basis[(i-1)^2+1]*Jy
        # compute x-only
        basis[(i+1)^2-1]=Jx*basis[i^2-1]
        # compute y-dominant
        for j=1:Int(i-1)
            basis[i^2+j+1]=Jx*basis[i^2+j]
        end
        # compute x-y-balanced
        basis[i^2]=Jx*basis[i^2-i]
        # compute x-dominant
        for j=1:Int(i-1)
            basis[(i+1)^2-j-1]=basis[(i+1)^2-j]*Jy
        end
    end
    basis[(k+1)^2]=Jx*basis[k^2-k]
    return basis
end
#############################################################
## The stuff below is currently not in use outside of testing
#############################
##   triJQEJbasisgridP01: Generates full basis of J^lQEJ^k in a grid.
##                        k and l are the highest orders of x and y in the kernel.
##                        basis[m][n] will give the m+1,n+1 element.
##                        Primarily for testing purposes, not recommended outside of testing.
####
function triJQEJbasisgridgenP01(k::Integer,l::Integer,N::Integer)
    S = Jacobi(0,1,0..1)
    x = Fun(identity,S)
    J = Multiplication(x, S)[1:N,1:N]
    QEy=triQEygenP01(N)
    # careful: basis[1][1] actually relects the 0,0 element and so on.
    basis=[]
    ybasisV = [QEy]
    for m=1:k
        for n=2:l
            push!(ybasisV,ybasisV[n-1]*J)
        end
        push!(basis,ybasisV)
        ybasisV=[J*basis[m][1]]
    end
    return basis
end
#############################
##   triNoKernelDirect: Direct integration from 0 to 1-x with K(x,y)=1 using only Q and ApproxFun, skipping the use of
##                      E and commutation relations.
##                      Primarily for testing purposes, not recommended outside of testing.
####
function triNoKernelDirect(f,x::Float64)
    # Computes integral of f(y) from 0 to 1-x via Jacobi polynomials on triangle.
    T = JacobiTriangle()
    f̃ = Fun((x,y) -> f(y), T)
        N = length(f̃.coefficients)
        v = PseudoBlockArray(pad(f̃.coefficients,sum(1:N)), 1:N)
    Q = triQgen(N)
    Qf = Fun(JacobiWeight(0,1, Jacobi(0,1, 0..1)), (1/2*Q*v))
    return Qf(1-x)
end