text:
# working with images ...
#
# (c) A. Dominik, 2020

const UNKNOWN_CLASS = "unknwon_class"

"""
    function mk_image_minibatch(dir, batchsize; split=false, fr=0.2,
                                balanced=false, shuffle=true, train=true,
                                pre_load=false,
                                aug_pipl=nothing, pre_proc=nothing)

Return one or two iterable image-loader-objects that provides
minibatches of images. For training each minibatch is a tupel
`(x,y)` with x: 4-d-array with the minibatch of data and y:
vector of class IDs as Int.

### Arguments:
+ `dir`: base-directory of the image dataset. The first level of
        sub-dirs are used as class names.
+ `batchsize`: size of minibatches

### Keyword arguments:
+ `split`: return two iterators for training and validation
+ `fr`: split fraction
+ `balanced`: return balanced data (i.e. same number of instances
        for all classes). Balancing is achieved via oversampling
+ `shuffle`: if true, shuffle the images everytime the iterator
        restarts
+ `train`: if true, minibatches with (x,y) Tuples are provided,
        if false only x (for prediction)
+ `pre_load`: if `true` all images are loaded in advance;
        otherwise images are loaded on demand durng training.
        (option is *not implemented yet!*)
+ `aug_pipl`: augmentation pipeline for Augmentor.jl. Augmentation
        is performed before the pre_proc-function is applied
+ `pre_proc`: function with preprocessing
        and augmentation algoritms of type x = f(x). In contrast
        to the augmentation that modifies images, is `pre_proc`
        working on Arrays{Float32}.
+ `pre_load=false`: read all images from disk once when populating the
        loader (requires loads of memory, but speeds up training).
"""
function mk_image_minibatch(dir, batchsize; split=false, fr=0.5,
                            balanced=false, shuffle=true, train=true,
                            pre_load=false,
                            aug_pipl=nothing, pre_proc=nothing)

    i_paths = get_files_list(dir)
    i_n = length(i_paths)

    if train
        i_class_names = get_class_names(dir, i_paths)
        classes = unique(i_class_names)
        i_classes = [findall(x->x==c, classes)[1] for c in i_class_names]
    else
        i_class_names = fill(UNKNOWN_CLASS, i_n)
        classes =[UNKNOWN_CLASS]
        i_classes = fill(0, i_n)
    end

    if split                    # return train_loader, valid_loader
        ((xvld,yvld),(xtrn,ytrn)) = do_split(i_paths, i_classes, at=fr)
        if balanced
            (xtrn,ytrn) = do_balance(xtrn, ytrn)
            (xvld,yvld) = do_balance(xvld, yvld)
        end
        trn_loader = ImageLoader(dir, xtrn, ytrn, classes,
                            batchsize, shuffle, train,
                            aug_pipl, pre_proc,
                            pre_load)

        vld_loader = ImageLoader(dir, xvld, yvld, classes,
                            batchsize, shuffle, train,
                            aug_pipl, pre_proc,
                            pre_load)
        return trn_loader, vld_loader
    else
        xtrn, ytrn = i_paths, i_classes
        if balanced
            (xtrn, ytrn) = do_balance(i_paths, i_classes)
        end
        trn_loader = ImageLoader(dir, xtrn, ytrn, classes,
                            batchsize, shuffle, train,
                            aug_pipl, pre_proc,
                            pre_load)
        return trn_loader
    end
end

"""
    function get_class_labels(d::DataLoader)

Extracts a list of class labels from a DataLoader.
"""
function get_class_labels(dl::DataLoader)
    return dl.classes
end

"""
    struct ImageLoader <: DataLoader
        dir
        i_paths
        i_classes
        classes
        batchsize
        shuffle
        train
        aug_pipl
        pre_proc
        pre_load
        i_images
    end

Iterable image loader to provide minibatches of images as
4-d-arrays (x,y,rgb,mb).
"""
mutable struct ImageLoader <: DataLoader
    dir                 # root dir
    i_paths             # list of all file image_paths
    i_classes           # list of classs IDs for each image
    classes             # unique list of class names
    batchsize           #
    shuffle             # if true: shuffle for each start
    i_sequence          # actual sequrence of images to take
    train               # if true: inlcude y in minibatches
    aug_pipl            # Augmentor.jl pipeline
    pre_proc            # function to process one 3d image tensor (RGB)
    pre_load            # load all images on init
    i_images            # list of all images; nothing if not predolad.

    function ImageLoader(dir, i_paths, i_classes, classes,
                         batchsize, shuffle, train,
                         aug_pipl, pre_proc, pre_load)

        i_seq = collect(1:length(i_paths))
        if pre_load
            i_images = pre_load_images(i_paths)
        else
            i_images = nothing
        end
        return new(dir, i_paths, i_classes, classes,
                    batchsize, shuffle, i_seq,
                    train,
                    aug_pipl, pre_proc,
                    pre_load, i_images)
    end
end

# two iterate funs for ImageLoader
#
function Base.iterate(il::ImageLoader)

    if il.shuffle
        Random.shuffle!(il.i_sequence)
    end
    state = 1
    return iterate(il, state)
end

# state is the index of the next image after the currect batch.
#
function Base.iterate(il::ImageLoader, state)

    # check if empty:
    #
    if state > length(il.i_paths)
        return nothing
    end

    # range for next minibatch:
    #
    n = length(il.i_paths)
    mb_start = state
    mb_size = mb_start + il.batchsize > n ? n-mb_start+1 : il.batchsize
    state = mb_start+il.batchsize

    return mk_image_mb(il, mb_start, mb_size), state
end

function mk_image_mb(il, mb_start, mb_size)

    if mb_size < 1
        return nothing
    end

    # nice way:
    image_ns = mb_start:mb_start+mb_size-1

    # make image index form shuffled list entry:
    #
    image_ns = il.i_sequence[image_ns]

    # avoid broadcasting of iterator il:
    #
    mb_i = cat(read_one_image.(image_ns, Ref(il))..., dims=4)
    mb_y = UInt8.(il.i_classes[image_ns])

    if CUDA.functional()
        mb_i = KnetArray(mb_i)
    end

    if il.train
        return mb_i, mb_y
    else
        return mb_i
    end
end

function read_one_image(i, il)

    if il.pre_load
        img = il.i_images[i]
    else
        img = Images.RGB.(Images.load(il.i_paths[i]))
    end

    if il.aug_pipl isa Augmentor.ImmutablePipeline ||
        il.aug_pipl isa Augmentor.ImageOperation
        img = Augmentor.augment(img, il.aug_pipl)
    end

    if isa(img[1,1], Colors.RGB)
        img = Float32.(permutedims(Images.channelview(img), (3,2,1)))
    else
        img = Float32.(Images.channelview(img))
    end

    if il.pre_proc != nothing && il.pre_proc isa Function
        img = il.pre_proc(img)
    end
    return(img)
end

function Base.length(il::ImageLoader)
    n = length(il.i_paths) / il.batchsize
    return ceil(Int, n)
end

# return a list of all files realtive to dir:
#
function get_files_list(dir)

    image_paths = String[]
    for (root, dirs, files) in walkdir(dir)
        if !isempty(files)
            append!(image_paths, joinpath.(root, files))
        end
    end
    # filter out non-images:
    #
    fi = r"(\.png$)|(\.jpg$)|(\.jpeg$)|(\.gif$)"i
    image_paths = filter(i->occursin(fi, i), image_paths)
    return image_paths
end

# extract class names from top-level dir of list of filenames:
#
function get_class_names(dir, image_paths)

    classes = String[]
    for image_path in image_paths
        image_path = replace(image_path, dir=>"")
        dirs = splitpath(image_path)
        while length(dirs) > 0 && !occursin(r"[A-Za-z1-9]", dirs[1])
            popfirst!(dirs)
        end
        if length(dirs) < 2
            class = UNKNOWN_CLASS
        else
            class = dirs[1]
        end
        push!(classes, class)
    end

    return classes
end

# read all images into mem:
#
function pre_load_images(i_paths)

    images = []
    for path in i_paths
        push!(images, Images.RGB.(Images.load(path)))
    end
    return images
end

"""
    function image2array(img)

Take an image and return a 3d-array for RGB and a 2d-array for grayscale
images with the colour channels as last dimension.
"""
function image2array(img)

    ch = Images.channelview(img)

    if length(size(ch)) == 2
     ch = reshape(ch, :,size(ch)...)
    end

    arr = Float32.(permutedims(ch, (3,2,1)))

    if CUDA.functional()
        arr = KnetArray(arr)
    end

    return arr
end

"""
    function array2image(arr)

Take a 3d-array with colour channels as last dimension or a 2d-array
and return an array of RGB or of Gray as Image.
"""
function array2image(arr)

    if length(size(arr)) == 2
        arr = permutedims(arr, (2,1))
        itype = Images.Gray

    elseif length(size(arr)) == 3
        arr = permutedims(arr, (3,2,1))
        if size(arr)[1] == 1
            arr = reshape(arr, size(arr)[2], size(arr)[3])
            itype = Images.Gray
        else
            itype = Images.RGB
        end
    end
    return Images.colorview(itype, arr)
end

"""
    function array2RGB(arr)

Take a 3d-array with colour channels as last dimension or a 2d-array
and return always an array of RGB as Image.
"""
function array2RGB(arr)

    img = array2image(arr)
    return Images.RGB.(img)
end


text:
module Astrodynamics

    using Datetime

    importall Base

    export State
    export julian, gregorian, JD2000, JD1950, MJD
    export ecctomea, meantoecc, ecctotrue, truetoecc
    export propagate
    export planets

    immutable State
        rv::Vector
        t::DateTime
        frame::String
        body::String
    end

    include("constants.jl")
    include("time.jl")
    include("frames.jl")
    include("math.jl")
    include("kepler.jl")

    iss = State([8.59072560e+02, -4.13720368e+03, 5.29556871e+03, 7.37289205e+00, 2.08223573e+00, 4.39999794e-01],
        datetime(2013,3,18,12,0,0), "eci", "earth")

    function propagate(s::State; method::String="kepler")
       mu = planets[s.body]["mu"]
       if method == "kepler"
           ele = elements(s.rv, mu)
           tend = period(ele[1], mu)
           dt = [0:tend]
           return cartesian(kepler(ele, dt, mu), mu)
       end
    end
end

ext: jl


text:
function run_agent_simulation(simulation::AgentSimulation, simulation_years::Int64)

    total_horizon = get_total_horizon(get_case(simulation))
    rolling_horizon = get_rolling_horizon(get_case(simulation))

    installed_capacity =  zeros(simulation_years)

    capacity_forward_years = get_capacity_forward_years(simulation)

    # Set initial capacity market profits considering forward capacity auctions
    if get_markets(simulation)[:Capacity]
        initial_existing_projects = vcat(get_existing.(get_investors(simulation))...)
        initial_capacity_prices = [70000.0, 80000.0]              # initial capacity prices - arbitrarily selected here - #TODO: need some meachanism to generate these
        for y in 1:capacity_forward_years - 1
            for project in initial_existing_projects
                if get_end_life_year(project) >= y
                    for product in get_products(project)
                        update_initial_capacity_revenues!(project, product, initial_capacity_prices, y)
                    end
                end
            end
        end
    end

    sys_UC = get_system_UC(simulation)
    sys_ED = get_system_ED(simulation)

    investors = get_investors(simulation)

    average_capital_cost_multiplier = Statistics.mean(get_cap_cost_multiplier.(investors))

    services = get_system_services(sys_UC)

    for iteration_year = 1:simulation_years

        yearly_horizon = min(total_horizon - iteration_year + 1, rolling_horizon)

        println("Year $(iteration_year)")
        set_iteration_year!(simulation, iteration_year)

        active_projects = deepcopy(get_activeprojects(simulation))

        installed_capacity = update_installed_cap!(installed_capacity,
                                                   active_projects,
                                                   iteration_year,
                                                   simulation_years)

        create_investor_predictions(investors,
                                          active_projects,
                                          iteration_year,
                                          yearly_horizon,
                                          get_data_dir(get_case(simulation)),
                                          get_results_dir(simulation),
                                          average_capital_cost_multiplier,
                                          get_zones(simulation),
                                          get_lines(simulation),
                                          get_peak_load(simulation),
                                          get_solver(get_case(simulation)),
                                          get_parallel_investors(get_case(simulation)),
                                          get_parallel_scenarios(get_case(simulation))
                                          )

        for investor in investors
            run_investor_iteration(investor,
                                    active_projects,
                                    iteration_year,
                                    yearly_horizon,
                                    simulation_years,
                                    get_start_year(get_case(simulation)),
                                    capacity_forward_years,
                                    get_data_dir(get_case(simulation)),
                                    sys_UC,
                                    services,
                                    get_solver(get_case(simulation))
                                    )
        end

        update_simulation_derating_data!(simulation)

        #Get all existing projects to calculate realized profits for energy and REC markets.
        all_existing_projects = vcat(get_existing.(get_investors(simulation))...)

        # Get all projects which are expected to be online for the forward capacity market auction.
        capacity_market_year = iteration_year + capacity_forward_years - 1
        capacity_market_projects = Project[]

        for project in get_activeprojects(simulation)
            end_life_year = get_end_life_year(project)
            construction_year = get_construction_year(project)
            if end_life_year >= capacity_market_year && construction_year <= capacity_market_year
                push!(capacity_market_projects, project)
            end
        end

        installed_capacity = update_installed_cap!(installed_capacity,
                                                   all_existing_projects,
                                                   iteration_year,
                                                   simulation_years)

        println("Current Installed Capacity = $(round(installed_capacity[iteration_year])) MW")

        #Find which markets to simulate.
        markets = union(hcat(get_markets.(get_investors(simulation))...))

        #Create realzed market prices for existing projects.
        realized_market_prices,
        realized_capacity_factors,
        realized_reserve_up_perc,
        realized_reserve_down_perc,
        capacity_accepted_bids,
        rec_accepted_bids = create_realized_marketdata(simulation,
                                                             sys_UC,
                                                             sys_ED,
                                                             markets,
                                                             all_existing_projects,
                                                             capacity_market_projects,
                                                             capacity_forward_years,
                                                             iteration_year,
                                                             simulation_years,
                                                             get_solver(get_case(simulation)),
                                                             get_results_dir(simulation))

        existing_project_types = unique(get_type.(get_tech.(all_existing_projects)))

        #Update forecasts and realized profits of all existing projects for each investor.

        for investor in get_investors(simulation)
            if iteration_year < simulation_years
                update_forecast!(get_forecast(investor), get_annual_growth(simulation)[:, iteration_year], iteration_year)
            end

            projects = get_projects(investor)
            for (i, project) in enumerate(projects)

                update_realized_profits!(project,
                                         realized_market_prices,
                                         realized_capacity_factors,
                                         realized_reserve_up_perc,
                                         realized_reserve_down_perc,
                                         capacity_accepted_bids,
                                         rec_accepted_bids,
                                         get_hour_weight(simulation),
                                         iteration_year,
                                         capacity_forward_years)

                update_annual_cashflow!(project, iteration_year)

                retire_old!(projects,
                            i,
                            project,
                            sys_UC,
                            get_data_dir(get_case(simulation)),
                            iteration_year)

            end
        end

        println("COMPLETED YEAR $(iteration_year)")
    end

    final_portfolio = vcat(get_existing.(get_investors(simulation))...)

    for project in final_portfolio
        extrapolate_profits!(project, simulation_years)
    end

    FileIO.save(joinpath(get_results_dir(simulation), "simulation_data.jld2"), "simulation_data", simulation)

    return
end


text:
using VertexEliminationOrder
using Test
using Graphs

G = smallgraph("house")
@info "Testing minwidth with house graph"
@test minwidth!(G) == ([1, 5, 4, 3, 2], 2)

G = smallgraph("house")
@info "Testing minwidth with house graph"
@test minfill!(G) == ([5, 1, 4, 3, 2], 2)


text:
__precompile__(true)
module cNODE

    # libraries
    using Distributed
    using Flux
    using Flux.Optimise: update!, apply!
    using DiffEqFlux
    using Zygote: gradient
    using MLDataPattern: eachbatch, splitobs, shuffleobs
    using DifferentialEquations: ODEProblem, solve, Tsit5
    using Distributions: Normal, Uniform
    using Statistics, LinearAlgebra
    using SparseArrays, SharedArrays
    using DelimitedFiles

    # functions
    include("./module/generator.jl")
    export getParameters, getIC, getGLV
    export getSteadyState, getRewiredNetwork
    export generate_data

    include("./module/trainer.jl")
    export FitnessLayer
    export getModel
    export predict
    export loss
    export train_reptile

    include("./module/loader.jl")
    export import_data
    export split_data

    # Parameters for module
    const synthetic_data = ["Strength","Connectivity","Universality","Noise","Rewiring"]
    const real_data = ["Drosophila_Gut","Soil_Vitro","Soil_Vivo","Human_Gut","Human_Oral","Ocean"]
    const S = [.1,.15,.2]
    const C = [0.5,0.75,1]
    const U = [0.1,1,2]
    const E = [.001,.01,.025]
    const R = [.01,.1,.25]
    export synthetic_data, real_data
    export S, C, U, E, R

end


text:
# # Bag of Words (BOW) model

using Knet: Param, @diff, grad, params, KnetArray
using Random: shuffle!

# Data has one instance per line, class and sentence separated by `|||`, for example:
#
#     4 ||| A warm , funny , engaging film .
#
# `readdata(file)` reads the file and returns an array of (x,y) pairs where x is a sequence
# of word ids for a sentence, and y is a tag id. Uses `w2i` and `t2i` functions to turn
# words and tags to integers.

function readdata(file)
    data = []
    for line in eachline(file)
        tag, sentence = split(strip(lowercase(line)), " ||| ")
        wordids = w2i.(split(sentence))
        tagid = t2i(tag)
        push!(data, (wordids, tagid))
    end
    return data
end

# Before reading the training data, we initalize the word->id and tag->id dictionaries and
# functions and the unknown word tag `"<unk>"`.

wdict = Dict()
tdict = Dict()
w2i(x) = get!(wdict, x, 1+length(wdict))
t2i(x) = get!(tdict, x, 1+length(tdict))
UNK = w2i("<unk>")

# Load the training data and peek at the first instance:

trn = readdata("../data/classes/train.txt")
first(trn)

# Before reading the dev/test data, we change the word->id function to return UNK for
# unknown words and tag->id function to error for unknown tags.

w2i(x) = get(wdict, x, UNK)     # unk if not found
t2i(x) = tdict[x]               # error if not found

# Load the dev/test data and print the number of instances:

dev = readdata("../data/classes/dev.txt")
tst = readdata("../data/classes/test.txt")
length.((trn, dev, tst))

# Use KnetArray to initialize parameters on GPU, use Array to initialize on CPU:

## param(dims...) = Param(Array(0.01f0 * randn(Float32, dims...)))
param(dims...) = Param(KnetArray(0.01f0 * randn(Float32, dims...)))

# Initialize the parameters of the BOW model as global variables W and b.

nwords = length(wdict)
ntags = length(tdict)
W = param(ntags, nwords)
b = param(ntags)

# Here is the prediction function for the BOW model:

pred(words) = b .+ sum(W[:,words], dims=2)

# Here is the loss function for the BOW model:

function loss(words, tag)
    scores = pred(words)
    logprobs = scores .- log(sum(exp.(scores)))
    -logprobs[tag]
end

# We use the following to report accuracy during training:

accuracy(data) = sum(argmax(pred(x))[1] == y for (x,y) in data) / length(data)

# Here is the SGD training loop:

function train(; nepochs = 10, lr = 0.01)
    for epoch in 1:nepochs
        shuffle!(trn)
        for (x,y) in trn
            ∇loss = @diff loss(x,y)
            for p in params(∇loss)
                p .= p - lr * grad(∇loss, p)
            end
        end
        println((epoch = epoch, trn = accuracy(trn), dev = accuracy(dev)))
    end
end


text:
module Solaris

export my_tst

"""
My test function

"""
function my_tst(a)
	println(a)
end

end


text:
module DIPTutorials

export DIP3e

include("DIP3e.jl")
using .DIP3e

"""
A helper package for DIP tutorials: https://johnnychen94.github.io/Image-Processing-in-Julia/dev/

# Examples

```julia
using DIPTutorials, FileIO
imgfile = DIP3e.get_filepath("0101") # figure 1.1
img = load(imgfile)
```
"""
DIPTutorials

end


text:
#### AppVeyor

export AppVeyor

"""
    AppVeyor

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct AppVeyor
    "Julia version"
    julia_version::String
    "Checks"
    checks::String="continuous-integration/appveyor/branch"
end

function yaml(ci::AppVeyor)
    return "environment:
  matrix:
    - julia_version: $(ci.julia_version)

platform:
  - x64 # 64-bit

notifications:
  - provider: Email
    on_build_success: false
    on_build_failure: false
    on_build_status_changed: false

install:
  - ps: iex ((new-object net.webclient).DownloadString(\"https://raw.githubusercontent.com/JuliaCI/Appveyor.jl/version-1/bin/install.ps1\"))
  - ps: \"[System.Net.ServicePointManager]::SecurityProtocol = [System.Net.SecurityProtocolType]::Tls12\"
  - ps: (new-object net.webclient).DownloadFile(
        \"https://download.microsoft.com/download/2/E/C/2EC96D7F-687B-4613-80F6-E10F670A2D97/msmpisetup.exe\",
        \"C:\\projects\\MSMpiSetup.exe\")
  - C:\\projects\\MSMpiSetup.exe -unattend -minimal
  - set PATH=C:\\Program Files\\Microsoft MPI\\Bin;%PATH%

build_script:
  - echo \"%JL_BUILD_SCRIPT%\"
  - C:\\julia\\bin\\julia -e \"%JL_BUILD_SCRIPT%\"

test_script:
  - echo \"%JL_TEST_SCRIPT%\"
  - C:\\julia\\bin\\julia -e \"%JL_TEST_SCRIPT%\"

"
end

max_stars_repo_name: charleskawczynski/PkgTemplates

max_issues_repo_name: charleskawczynski/PkgTemplates

max_forks_repo_name: charleskawczynski/PkgTemplates


text:
import Base

include("utils.jl")

"""
Supertype for productions.
[`Production`](@ref) and other types are subtypes of this.
"""
abstract type AbstractProduction end

"""
A grammar production, consisting of a left side and a right side
"""
struct Production <: AbstractProduction
    "The left-hand side of the production"
    left::Union{AbstractString, AbstractArray}
    "The right-hand side of the production"
    right::Union{AbstractString, AbstractArray}
    function Production(left, right)
        l = nothing
        r = nothing
        if typeof(left) <: AbstractString && left ≠ nothing && ~isempty(left)
            l = left
        elseif (typeof(left) <: AbstractArray || typeof(left) <: AbstractSet || typeof(left) <: Tuple) && ~isempty(left) && all(x-> typeof(x) <: AbstractString && x ≠ nothing && ~isempty(x), left)
            l = collect(left)
        else
            throw(ArgumentError("The left-hand side is not a nonempty String, nor an iterable of nonempty strings."))
        end

        if typeof(right) <: AbstractString && right ≠ nothing && ~isempty(left)
            r = right
        elseif (typeof(right) <: AbstractArray || typeof(right) <: AbstractSet || typeof(right) <: Tuple) && ~isempty(right) && all(x-> typeof(x) <: AbstractString && x ≠ nothing && ~isempty(x), right)
            r = collect(right)
        else
            throw(ArgumentError("The right-hand side is not an iterable of nonempty strings"))
        end

        if "ε" ∈ r && length(r) ≠ 1
            throw(ArgumentError("The right-hand side contains ε but has more than one symbol"))
        end
        new(l, r)
    end
end

### Functions ###

"""
    parseproduction(input::AbstractString, iscontextfree::Bool = true)::Array{Production}
Returns an Array of [`Production`](@ref) obtained from the given string.
"""
function parseproduction(input::AbstractString, iscontextfree::Bool = true)::Array{Production}
    P = []
    for p ∈ filter(f::AbstractString -> f ≠ "",split(input, "\n"))
        l, r = split(p, "->")
        left = split(l)
        if iscontextfree
            if length(left) ≠ 1
                throw(ArgumentError("Production " * p * " marked as context free while it's not"))
            end
            left = left[1]
        end
        for right ∈ split(r, "|")
            push!(P, Production(left, split(right)))
        end
    end
    return P
end

"""
    suchthat(;left::AbstractString = nothing, right::AbstractString = nothing, rightlen::Int = nothing, right_is_suffix_of::AbstractString = nothing)::Array
A predicate (that is a one-argument function that retuns `true` or `false`) that is `true` weather the production
given as argument satisfies all the predicates given by the named arguments.
"""
function suchthat(;left::Union{AbstractString,Nothing} = nothing, right::Union{AbstractString,Nothing} = nothing, rightlen::Union{Int,Nothing} = nothing, right_is_suffix_of::Union{Iterable,Nothing} = nothing)
    c = []
    if left ≠ nothing
        push!(c, p::AbstractProduction -> p.left == left)
    end
    if right ≠ nothing
        push!(c, p::AbstractProduction -> p.right == [right])
    end
    if rightlen ≠ nothing
        push!(c, p::AbstractProduction -> length(p.right) == rightlen)
    end
    if right_is_suffix_of ≠ nothing
        d = collect(right_is_suffix_of)
        push!(c, p::AbstractProduction -> d[(end - (length(p.right) -1)):end] == collect(p.right))
    end
    return p::AbstractProduction -> all(cond(p) for cond ∈ c)
end

"""
    astype0(p::Production)::Production
Returns a new [`Production`](@ref) that is type 0
"""
astype0(p::Production)::Production = typeof(p.left) <: AbstractArray ? p : Production([p.left], p.right)

### Operators ###

Base.:(==)(x::Production, y::Production) = (x.left, x.right) == (y.left, y.right)

Base.:<(x::Production, y::Production) = (x.left, x.right) < (y.left, y.right)

Base.hash(x::Production) = Base.hash((x.left, x.right))

Base.show(io::IO, x::Production) = Base.show(io, string(x.left, "->", x.right))

Base.iterate(x::Production, i...) = Base.iterate((x.left, x.right), i...)

Base.isless(x::Production, y::Production) = Base.isless((x.left, x.right), (y.left, y.right))


text:
function load_scene_data(scene_full_path::AbstractString)
    labels = Int[]
    X = Float64[]
    Y = Float64[]
    Z = Float64[]
    B = Float64[]
    G = Float64[]
    L = Int[]
    asset_names = String[]
    skymap_name = "nothing"
    assets_fn = "nothing"
    optic_name = "nothing"
    println( string("Opening scene file: ", scene_full_path, " ..."))
    open(scene_full_path) do f
      for line in eachline(f)
        if length(line) > 1
          if occursin("//", line)
            continue
          elseif occursin("#", line)
            label, asset_name = split(line, ",")
            labels = push!(labels, parse(Int, label[2:end]) )
            asset_names = push!(asset_names, asset_name)
          elseif occursin("@", line)
            assets_fn = split(line, ",")[2]
          elseif occursin("*", line)
            skymap_name = split(line, ",")[2]
          elseif occursin("&", line)
            optic_name = split(line, ",")[2]
          else
            x_, y_, z_, beta_, gamma_, label_ = split(line, ",")
            x, y, z = parse(Float64, x_), parse(Float64, y_), parse(Float64, z_)
            beta, gamma = parse(Float64, beta_), parse(Float64, gamma_)
            l = round(Int, parse(Float64, label_))
            X = push!(X, x)
            Y = push!(Y, y)
            Z = push!(Z, z)
            B = push!(B, beta)
            G = push!(G, gamma)
            L = push!(L, l)
          end
        end
      end
    end
    return skymap_name, optic_name, assets_fn, labels, asset_names, L, X, Y, Z, B, G
end

function read_scene(scene_name::AbstractString)
    scene_full_path = string(scenefolder, scene_name, ".jay")
    skymap_name, optic_name, assets_fn, asset_labels, asset_names, asset_idx, x, y, z, beta, gamma = load_scene_data(scene_full_path)

    # load assets...
    xml_path = string(assetfolder, assets_fn, ".xml")
    xdoc = parse_file(xml_path)
    xroot = LightXML.root(xdoc)
    unique_mesh_names, mesh_idx, unique_palette_names, palette_idx = get_assets(xroot, asset_names)
    asset2mesh_idx = mesh_idx[asset_idx]
    asset2palette_idx = palette_idx[asset_idx]
    n = length(x)
    assets = Asset{Float64, Int}[]
    for i=1:n
        push!(assets, Asset(x[i], y[i], z[i], beta[i], gamma[i], asset2mesh_idx[i], asset2palette_idx[i]))
    end

    # load mesh data...
    meshes::Array{Mesh{Float64, Int, Array{Float64,1}, Array{Int64,1}},1} = []
    mesh_bvhs::Array{Bvh{Float64, Int, Array{Float64,1}, Array{Int64,1}},1} = []
    for mesh_name in unique_mesh_names
        mesh = read_mesh(mesh_name)
        bvh = read_bvh(mesh_name)
        meshes = push!(meshes, mesh)
        mesh_bvhs = push!(mesh_bvhs, bvh)
    end

    # load scene bvh (or make it if necessary)...
    if !isfile( string(bvhfolder, scene_name, ".bvh") )
        make_bvh(assets, meshes, mesh_bvhs)
    end
    scene_bvh::Bvh{Float64, Int, Array{Float64,1}, Array{Int,1}} = read_bvh(scene_name)

    # load the camera...
    camera = read_camera(optic_name)

    # read the skymap (optional, only if specified in scene description file)...
    skymap_path = string(skyfolder, skymap_name, ".csv")
    skymap = readdlm(skymap_path, ',', Float64)

    # load materials/shaders...
    palettes = Vector{Function}[]
    for palette_name in unique_palette_names
        xmlfn = string(palettefolder, palette_name, ".xml")
        shaders = load_shaders(xmlfn ; assert_nBands=camera.nBands) # shaders, or palette, is of type: Vector{Function}
        palettes = push!(palettes, shaders)
    end

    return assets, meshes, palettes, mesh_bvhs, scene_bvh, skymap, camera
end

max_stars_repo_name: martinvanleeuwen/RenderJay

max_issues_repo_name: martinvanleeuwen/RenderJay

max_forks_repo_name: martinvanleeuwen/RenderJay


text:
using TomoForward
using XfromProjections

img = zeros(128, 128)
img[40:60, 40:60] .= 1.0
H, W = size(img)

nangles = 90
detcount = Int(floor(size(img,1)*1.4))
proj_geom = ProjGeom(1.0, detcount, LinRange(0,pi,nangles+1)[1:nangles])

A = fp_op_parallel2d_strip(proj_geom, size(img, 1), size(img, 2))
p = reshape(A * vec(img), nangles, detcount)

u = zeros(size(img))
niter=100
@time recon2d_sirt!(u, A, p, niter)

# using PyPlot
# imshow(u, cmap="gray")
# title("reconstructed image")
# show()

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
"""
An MeshFileGeometry represents a mesh which is stored as the raw contents
of a file, rather than as a collection of points and vertices. This is useful for
transparently passing mesh files which we can't load in Julia directly to meshcat.

Supported formats:
    * .stl (ASCII and binary)
    * .obj
    * .dae (Collada)

For .obj and .dae files, only a single mesh geometry will be loaded, and any
material or texture properties will be ignored. To load an entire collection of
objects (complete with materials and textures) from an .obj or .dae file, see
MeshFileObject instead.
"""
struct MeshFileGeometry{S <: Union{String, Vector{UInt8}}}
    contents::S
    format::String
end

function MeshFileGeometry(filename)
    ext = lowercase(splitext(filename)[2])
    if ext ∈ (".obj", ".dae")
        MeshFileGeometry(open(f -> read(f, String), filename), ext[2:end])
    elseif ext == ".stl"
        MeshFileGeometry(open(read, filename), ext[2:end])
    else
        throw(ArgumentError("Unsupported extension: $ext. Only .obj, .dae, and .stl meshes can be used to construct MeshFileGeometry"))
    end
end

"""
A MeshFileObject is similar to a MeshFileGeometry, but rather than representing
a single geometry, it supports loading an entire object (with geometries,
materials, and textures), or even a collection of objects (for supported file
types).

Supported formats:
    * .obj
    * .dae (Collada)

Since mesh files may include references to other files (such as texture images),
the MeshFileObject also includes a `resources` dictionary mapping relative paths
(as specified in the mesh file) to the contents of the referenced files. Those contents are specified using data URIs, e.g.:

    data:image/png;base64,<base-64 encoded data here>

If you construct a MeshFileObject from a `.dae` or `.obj` file, the resources
dictionary will automatically be populated.
"""
struct MeshFileObject <: AbstractObject
    contents::String
    format::String
    mtl_library::Union{String, Nothing}
    resources::Dict{String, String}
end

function MeshFileObject(filename)
    ext = lowercase(splitext(filename)[2])
    if ext ∉ (".obj", ".dae")
        throw(ArgumentError("Unsupported extension: $ext. Only .obj and .dae meshes can be used to construct MeshFileObject"))
    end
    contents = open(f -> read(f, String), filename)
    format = ext[2:end]
    if ext == ".obj"
        mtl_library = load_mtl_library(contents, dirname(filename))
        resources = load_mtl_textures(mtl_library, dirname(filename))
    else
        mtl_library = nothing
        resources = load_dae_textures(contents, dirname(filename))
    end
    MeshFileObject(contents, format, mtl_library, resources)
end

function load_mtl_library(obj_contents, directory=".")
    libraries = String[]
    for line in eachline(IOBuffer(obj_contents))
        m = match(r"^mtllib (.*)$", line)
        if m !== nothing
            push!(libraries, open(f -> read(f, String), joinpath(directory, m.captures[1])))
        end
    end
    join(libraries, '\n')
end

data_uri(contents::Vector{UInt8}) = string("data:image/png;base64,", base64encode(contents))

function load_mtl_textures(mtl_contents, directory=".")
    textures = Dict{String, String}()
    for line in eachline(IOBuffer(mtl_contents))
        m = match(r"^map_[^ ]* (.*)$", line)
        if m !== nothing
            name = m.captures[1]
            contents = open(read, joinpath(directory, name))
            textures[name] = data_uri(contents)
        end
    end
    textures
end

function load_dae_textures(dae_contents, directory=".")
    # TODO: this is probably not a very robust parsing strategy,
    # but I don't have enough examples to work from to decide on a
    # better one, so I'm keeping it simple.
    r = r"\<image [^\>]*\>\s*<init_from\>([^\<]*)\<\/init_from\>"m
    textures = Dict{String, String}()
    for match in eachmatch(r, dae_contents)
        name = match.captures[1]
        textures[name] = data_uri(open(read, joinpath(directory, name)))
    end
    textures
end


text:
# Copyright (c) 2022 Dongdong Kong. All rights reserved.
# This work is licensed under the terms of the MIT license.
# For a copy, see <https://opensource.org/licenses/MIT>.

module nctools

# using NetCDF
export NcDim

using NCDatasets
export nc_open, close

using CFTime
using Dates

# include("Ipaper/Ipaper.jl")
using Ipaper
export Ipaper

include("nc_info.jl")
include("nc_dim.jl")
include("nc_date.jl")
include("nc_read.jl")
include("nc_write.jl")
include("ncvar_def.jl")
include("ncdim_def.jl")
include("ncatt_put.jl")

# CMIP
include("CMIP/CMIP.jl")
export CMIP

export nc_open, nc_close, nc_bands, nc_date, nc_info, ncinfo
export nc_dim, nc_dims, ncvar_dim, nc_size, nc_cellsize
export nc_read, nc_write, nc_write!
export nc_date, nc_calendar

export ncdim_def, ncvar_def
export ncatt_put, ncatt_get, ncatt_del

export names

end # module


text:
# SPDX-License-Identifier: MIT

"""
    abstract type Semiring <: Number end

Abstract base type for all semirings.
"""
abstract type Semiring <: Number end

"""
    abstract type Semifield <: Semiring end

Abstract base type for all semifields.
"""
abstract type Semifield <: Semiring end

Base.zero(x::T) where T<:Semiring = zero(T)
Base.one(x::T) where T<:Semiring = one(T)
Base.convert(::Type{T}, x::T) where T<:Semiring = x
Base.convert(T::Type{<:Number}, x::Semiring) = T(x.val)
Base.promote_rule(x::Type{T}, y::Type{<:Real}) where T <: Semiring = T
Base.show(io::IO, x::Semiring) = print(io, x.val)


text:
using Word2Vec
using LinearAlgebra
using MultivariateStats
using Plots

function loadNonStopWords(fileName::String)::Tuple{Vector{String},Dict{String,Int}}
    file = open(fileName)
    nonStopWords = Vector{String}()
    nonStopWordIDs = Dict{String,Int}()
    for phrase in readlines(file)
        if phrase != ""
            push!(nonStopWords, phrase)
            nonStopWordIDs[phrase] = length(nonStopWords)
        end
    end
    close(file)
    (nonStopWords, nonStopWordIDs)
end

function normalizeVectors(wordVectors::WordVectors, nonStopWords::Vector{String})::Array{Float64,2}
    numDims, numWords = size(wordVectors)
    numVectors = length(nonStopWords)
    result = Array{Float64,2}(undef, (numDims, numVectors))
    for i in 1:numVectors
        vector = get_vector(wordVectors, nonStopWords[i])
        w = 1.0 / norm(vector)
        result[:,i] = vector * w
    end
    result
end

function loadPhraseIDGroups(fileName::String, nonStopWordIDs::Dict{String,Int})::Vector{Vector{Int}}
    file = open(fileName)
    phraseIDGroups = Vector{Vector{Int}}()
    for line in readlines(file)
        if length(line) == 0
            continue
        end
        phrases = split(line)
        idGroup = Vector{Int}()
        for phrase in phrases
            if haskey(nonStopWordIDs, phrase)
                push!(idGroup, nonStopWordIDs[phrase])
            end
        end
        push!(phraseIDGroups, idGroup)
    end
    close(file)
    phraseIDGroups
end

function getCardinalities(nonStopWords::Vector{String}, nonStopWordIDs::Dict{String,Int}, phraseIDGroups::Vector{Vector{Int}})::Vector{Int}
    numPhrases = length(nonStopWords)
    idMap = Dict{Int,Int}()
    for (newID, phrase) in enumerate(nonStopWords)
        idMap[nonStopWordIDs[phrase]] = newID
    end

    phraseFreqs = fill(0, numPhrases)
    for idGroup in phraseIDGroups
        for id in idGroup
            phraseFreqs[idMap[id]] += 1
        end
    end

    phraseFreqs
end

function loadAspects(fileName::String, nonStopWordIDs::Dict{String,Int})::Tuple{Vector{Set{Int}},Vector{Int}}
    file = open(fileName)
    aspects = Vector{Set{Int}}()
    numPhrases = length(nonStopWordIDs)
    aspectIDs = Vector{Int}(undef, numPhrases)
    for line in readlines(file)
        fields = split(line)
        if length(fields) != 3
            continue
        end
        aspectID = parse(Int, fields[1])
        card = parse(Int, fields[2])
        phrase = fields[3]
        phraseID = nonStopWordIDs[phrase]
        if aspectID > length(aspects)
            @assert aspectID == length(aspects) + 1
            push!(aspects, Set{Int}())
        end
        push!(aspects[aspectID], phraseID)
        aspectIDs[phraseID] = aspectID
    end
    close(file)
    (aspects, aspectIDs)
end

function plotAspects(aspects::Vector{Set{Int}}, vectors::Array{Float64,2}, pca::PCA{Float64}, dim1::Integer, dim2::Integer)
    numAspects = length(aspects)
    numDims = size(vectors, 1)
    aspectVectors = Vector{Array{Float64,2}}(undef, numAspects)
    aspectFreqs = Vector{Vector{Int}}(undef, numAspects)
    for (aspectID, aspect) in enumerate(aspects)
        aspectSize = length(aspect)
        aspectPhraseIDs = [id for id in aspect]
        aspectVectors[aspectID] = Array{Float64,2}(undef, (numDims, aspectSize))
        for (localID, globalID) in enumerate(aspectPhraseIDs)
            aspectVectors[aspectID][:,localID] = vectors[:,globalID]
        end
    end
    result1 = []
    result2 = []
    for aspectID in 1:numAspects
        aspectPCA = transform(pca, aspectVectors[aspectID])
        push!(result1, aspectPCA[dim1,:])
        push!(result2, aspectPCA[dim2,:])
        #p = scatter(aspectPCA[1,:],aspectPCA[2,:], marker=:circle, linewidth=0)
        #plot!(p,xlabel="PC1",ylabel="PC2")
    end
    (result1, result2)
end

wordVectors = wordvectors("DualAspectTopicRepresentation/workspace/NIPS.awe.model.txt")
nonStopWords, nonStopWordIDs = loadNonStopWords("DualAspectTopicRepresentation/workspace/NIPS.nonStopWords.txt")
vectors = normalizeVectors(wordVectors, nonStopWords)
phraseIDGroups = loadPhraseIDGroups("DualAspectTopicRepresentation/workspace/NIPS.w2v.training.txt", nonStopWordIDs)
cardinalities = getCardinalities(nonStopWords, nonStopWordIDs, phraseIDGroups)
aspects1, aspectIDs1 = loadAspects("DualAspectTopicRepresentation/workspace/NIPS.aspects.1.txt", nonStopWordIDs)
#aspects2, aspectIDs2 = loadAspects("DualAspectTopicRepresentation/workspace/NIPS.aspects.2.txt", nonStopWordIDs)
#aspects3, aspectIDs3 = loadAspects("DualAspectTopicRepresentation/workspace/NIPS.aspects.3.txt", nonStopWordIDs)
pca = fit(PCA, vectors)
X, Y = plotAspects(aspects1, vectors, pca, 1, 2)
scatter(X,Y)
#X, Y = plotAspects(aspects2, vectors, pca, 1, 2)
#scatter(X,Y)
#X, Y = plotAspects(aspects3, vectors, pca, 1, 2)
#scatter(X,Y)

X, Y = plotAspects(aspects1, vectors, pca, 1, 3)
scatter(X,Y)
#X, Y = plotAspects(aspects2, vectors, pca, 1, 3)
#scatter(X,Y)
#X, Y = plotAspects(aspects3, vectors, pca, 1, 3)
#scatter(X,Y)

X, Y = plotAspects(aspects1, vectors, pca, 2, 3)
scatter(X,Y)
#X, Y = plotAspects(aspects2, vectors, pca, 2, 3)
#scatter(X,Y)
#X, Y = plotAspects(aspects3, vectors, pca, 2, 3)
#scatter(X,Y)

max_stars_repo_path: src/08 plot_aspects.jl

max_issues_repo_path: src/08 plot_aspects.jl

max_forks_repo_path: src/08 plot_aspects.jl


text:
using PSSFSS
using PSSFSS.Modes: choose_layer_modes!
#using PSSFSS.Elements: rectstrip, polyring
using PSSFSS.Layers: Layer
using PSSFSS.Sheets: RWGSheet
#using PSSFSS.PSSFSSLen: mm, cm, inch, mil
using PSSFSS.GSMs: Gblock, choose_gblocks
using PSSFSS.Constants: c₀, twopi
using Test
using Logging: Error, ConsoleLogger, default_metafmt, global_logger

testlogger = ConsoleLogger(stderr, Error,
                       meta_formatter=default_metafmt, show_limited=true,
                       right_justify=0)
oldlogger = global_logger(testlogger)

@testset "choose_layer_modes!1" begin
    strata =
    [
        Layer(name="Vacuum")
        Layer(name="Duroid 6035", ϵᵣ=3.6, tanδ=0.0013, width=20mil)
        polyring(sides=4, orient=45, a=[0.1902], b=[0.2169], units=inch,
                    s1=[0.5117, 0], s2=[0, 0.5117], ntri=100)
        Layer(name="Duroid 6035", ϵᵣ=3.6, tanδ=0.0013, width=20mil)
        Layer(name="Vacuum", width=0.2inch)
        Layer(name="Duroid 6035", ϵᵣ=3.6, tanδ=0.0013, width=20mil)
        polyring(sides=4, orient=45, a=[0.1902], b=[0.2169], units=inch,
                    s1=[0.5117, 0], s2=[0, 0.5117], ntri=100)
        Layer(name="Duroid 6035", ϵᵣ=3.6, tanδ=0.0013, width=20mil)
        Layer(name="Vacuum")
    ]

    islayer = map(x -> x isa Layer, strata)
    issheet = map(x -> x isa RWGSheet, strata)
    layers = convert(Vector{Layer}, strata[islayer])
    sheets = convert(Vector{RWGSheet}, strata[issheet])
    nl = length(layers)
    nj = nl - 1
    ns = length(sheets)
    sint = cumsum(islayer)[issheet] # sint[k] contains dielectric interface number of k'th sheet
    junc = zeros(Int, nj)
    junc[sint] = 1:ns #  junc[i] is the sheet number present at interface i, or 0 if no sheet is there
    fmin = 2.0 * 1e9  # minimum frequency [Hz]
    k0min = twopi * fmin / c₀
    fmax = 8.0 * 1e9  # minimum frequency [Hz]
    k0max = twopi * fmax / c₀
    dbmin = 25.0

    gbl = choose_gblocks(layers, sheets, junc, k0min)
    choose_layer_modes!(layers, sheets, junc, gbl, k0max, dbmin)
    @test [length(l.P) for l in strata if l isa Layer] == [2,0,0,18,0,0,2]
end

@testset "choose_layer_modes!2" begin
    strata =
    [
        Layer(name="Vacuum")
        Layer(name="Duroid 6035", ϵᵣ=3.6, tanδ=0.0013, width=20mil)
        polyring(sides=4, orient=45, a=[0.1902], b=[0.2169], units=inch,
                    s1=[0.5117, 0], s2=[0, 0.5117], ntri=100)
        Layer(name="Duroid 6035", ϵᵣ=3.6, tanδ=0.0013, width=20mil)
        Layer(name="Vacuum", width=0.1inch)
        Layer(name="Vacuum", width=0.1inch)
        Layer(name="Duroid 6035", ϵᵣ=3.6, tanδ=0.0013, width=20mil)
        polyring(sides=4, orient=45, a=[0.1902], b=[0.2169], units=inch,
                    s1=[0.9005, 0], s2=[0, 0.9005], ntri=100)
        Layer(name="Duroid 6035", ϵᵣ=3.6, tanδ=0.0013, width=20mil)
        Layer(name="Vacuum")
    ]

    islayer = map(x -> x isa Layer, strata)
    issheet = map(x -> x isa RWGSheet, strata)
    layers = convert(Vector{Layer}, strata[islayer])
    sheets = convert(Vector{RWGSheet}, strata[issheet])
    nl = length(layers)
    nj = nl - 1
    ns = length(sheets)
    sint = cumsum(islayer)[issheet] # sint[k] contains dielectric interface number of k'th sheet
    junc = zeros(Int, nj)
    junc[sint] = 1:ns #  junc[i] is the sheet number present at interface i, or 0 if no sheet is there
    fmin = 2.0 * 1e9  # minimum frequency [Hz]
    k0min = twopi * fmin / c₀
    fmax = 8.0 * 1e9  # minimum frequency [Hz]
    k0max = twopi * fmax / c₀
    dbmin = 25.0

    gbl = choose_gblocks(layers, sheets, junc, k0min)
    choose_layer_modes!(layers, sheets, junc, gbl, k0max, dbmin)
    @test [length(l.P) for l in strata if l isa Layer] == [2,0,0,2,0,0,0,2]
end

@testset "choose_layer_modes!3" begin
    strata =
    [
        Layer(name="Vacuum")
        Layer(name="Duroid 6035", ϵᵣ=3.6, tanδ=0.0013, width=20mil)
        polyring(sides=4, orient=45, a=[0.1902], b=[0.2169], units=inch,
                    s1=[0.5117, 0], s2=[0, 0.5117], ntri=100)
        Layer(name="Duroid 6035", ϵᵣ=3.6, tanδ=0.0013, width=20mil)
        Layer(name="Vacuum", width=0.1inch)
        Layer(name="Vacuum", width=0.1inch)
        Layer(name="Duroid 6035", ϵᵣ=3.6, tanδ=0.0013, width=20mil)
        polyring(sides=4, orient=45, a=[0.1902], b=[0.2169], units=inch,
                    s1=[0.5117, 0], s2=[0, 0.5117], ntri=100)
        Layer(name="Duroid 6035", ϵᵣ=3.6, tanδ=0.0013, width=20mil)
        Layer(name="Vacuum")
    ]

    fmin = 2.0 * 1e9  # minimum frequency [Hz]
    k0min = twopi * fmin / c₀
    fmax = 8.0 * 1e9  # minimum frequency [Hz]
    k0max = twopi * fmax / c₀
    dbmin = 25.0
    islayer = map(x -> x isa Layer, strata)
    issheet = map(x -> x isa RWGSheet, strata)
    layers = convert(Vector{Layer}, strata[islayer])
    sheets = convert(Vector{RWGSheet}, strata[issheet])
    nl = length(layers)
    nj = nl - 1
    ns = length(sheets)
    sint = cumsum(islayer)[issheet] # sint[k] contains dielectric interface number of k'th sheet
    junc = zeros(Int, nj)
    junc[sint] = 1:ns #  junc[i] is the sheet number present at interface i, or 0 if no sheet is there

    gbl = choose_gblocks(layers, sheets, junc, k0min)
    choose_layer_modes!(layers, sheets, junc, gbl, k0max, dbmin)
    #@test_logs [(:info,),(:info,),(:info,)] choose_layer_modes!(strata, gbl, k0max, dbmin)
    @test [length(l.P) for l in strata if l isa Layer] == [2,0,0,42,0,0,0,2]

end

global_logger(oldlogger)
nothing


text:
## Test types ##
module Test_types

using MultipleTesting
using Test

@testset "Types" begin

    @testset "PValues type" begin

        n = 10
        vals = rand(n)

        pv = PValues(vals)

        # basic vector functionality
        @test values(pv) == vals
        @test sum(pv) == sum(vals)
        @test length(pv) == n
        @test zero(pv) == zeros(eltype(vals), n)
        @test minimum(pv) == minimum(vals)
        @test maximum(pv) == maximum(vals)
        @test extrema(pv) == extrema(vals)
        @test [x for x in pv] == vals
        @test convert(PValues, vals) == PValues(vals)

        # min/max are taken from the type fields,
        # not recomputed from the values
        pvt = PValues(vals)
        @test_throws ErrorException pvt.min = 0.0
        @test_throws ErrorException pvt.max = 1.0

        # parametric type that keeps input float type
        for T in (Float16, Float32, Float64)
            vals = rand(T, n)
            pv = PValues(vals)
            @test eltype(pv) == T
            @test eltype(values(pv)) == T
        end

        # min/max cannot be passed to constructor
        @test_throws MethodError PValues(vals, 0.0, 1.0)

        # invalid p-values are caught
        @test_throws DomainError PValues([0.1, 0.2, 2.0])
        @test_throws DomainError PValues([0.1, 0.2, -2.0])

        # TODO discuss expected behaviour in edge cases
        @test_throws MethodError PValues(0.5)
        @test_throws ArgumentError PValues([])
        @test_throws TypeError PValues([0, 1])

    end

    @testset "PValues immutability" begin

        n = 10
        vals = rand(n)
        ref = deepcopy(vals)

        pv = PValues(vals)

        @test pv.values == ref
        @test extrema(pv) == extrema(pv.values)

        vals[1] = 0.0
        vals[2] = 1.0

        @test pv.values == ref
        @test pv.values != vals
        @test extrema(pv) == extrema(pv.values)

        @test_throws ErrorException pv[1] = 0.0

    end

end

end


text:
using NiLang.AD
export bennett_step!, bennett_solve
export SeismicState

struct SeismicState{MT}
    upre::MT
    u::MT
    φ::MT
    ψ::MT
    step::Base.RefValue{Int}
end
## solving gradient
function SeismicState(::Type{T}, nx::Int, ny::Int) where T
    SeismicState([zeros(T, nx+2, ny+2) for i=1:4]..., Ref(0))
end
Base.zero(x::SeismicState) = SeismicState(zero(x.upre), zero(x.u), zero(x.φ), zero(x.ψ), Ref(0))
Base.copy(x::SeismicState) = SeismicState(copy(x.upre), copy(x.u), copy(x.φ), copy(x.ψ), Ref(x.step[]))

@i function :(+=)(identity)(x::SeismicState, y::SeismicState)
    x.upre += y.upre
    x.u += y.u
    x.φ += y.φ
    x.ψ += y.ψ
    x.step[] += y.step[]
end

@i function bennett_step!(dest, src, param::AcousticPropagatorParams, srcloc, srcv, c)
    @routine begin
        d2 ← zero(param.DELTAT)
        d2 += param.DELTAT^2
    end
    dest.upre += src.u
    dest.step[] += src.step[] + 1
    i_one_step!(param, dest.u, src.u, src.upre,
        dest.φ, src.φ, dest.ψ, src.ψ, c)
    dest.u[srcloc...] += srcv[dest.step[]] * d2
    ~@routine
end

function bennett_solve(s0, gn; param, src, srcv, c, k, N, logger=NiLang.BennettLog())
    # forward execution
    d = Dict(1=>s0)
    bennett!(bennett_step!, d, param, src, srcv, c; k=k, N=N, logger=logger)
    # backward execution
    y = d[N+1]
    y = SeismicState([GVar(getfield(y, field), getfield(gn, field)) for field in fieldnames(SeismicState)[1:end-1]]..., Ref(y.step[]))
    d = GVar(d)
    d[N+1] = y
    _, gs, _, _, gv, gc = (~bennett!)(bennett_step!, d, param, src, GVar(srcv), GVar(c); k=k, N=N, logger=logger)
    return grad(gs[1]), grad(gv), grad(gc)
end


text:
# Test writing component of SegyIO
global key = ["GroupX";"GroupY";"SourceX";"SourceY"]
global dir = "../data/"
global file_filter = "test_write"

@testset "write" begin

	@testset "segy_write" begin
		block = SeisBlock(rand(Float32,10,10));
		set_header!(block, "dt", 6000);
		set_header!(block, :SourceX, 4000);
		set_header!(block, :SourceY, 4500);

		set_header!(block, :GroupX, 8000);
		set_header!(block, :GroupY, 8500);

		segy_write("../data/test_write.segy", block)

		s = segy_scan(dir,file_filter,key)

		@test size(s[1].data) == (10,10)
		@test get_header(s, "SourceX")[1] == 4000
                @test get_header(s, "SourceY")[1] == 4500
                @test get_header(s, "GroupX")[1] == 8000
                @test get_header(s, "GroupY")[1] == 8500
		@test get_header(s, "SourceX")[end] == 4000
                @test get_header(s, "SourceY")[end] == 4500
                @test get_header(s, "GroupX")[end] == 8000
                @test get_header(s, "GroupY")[end] == 8500
	end

	@testset "segy_write_append" begin
		block_append = SeisBlock(rand(Float32,10,10));

		set_header!(block_append, "dt", 6000);
                set_header!(block_append, :SourceX, 4000);
                set_header!(block_append, :SourceY, 4500);

                set_header!(block_append, :GroupX, 8000);
                set_header!(block_append, :GroupY, 9000);

		segy_write_append("../data/test_write.segy", block_append)

		s = segy_scan(dir,file_filter,key)

                @test size(s[1].data) == (10,20)
		@test get_header(s, "SourceX")[1] == 4000
                @test get_header(s, "SourceY")[1] == 4500
                @test get_header(s, "GroupX")[1] == 8000
                @test get_header(s, "GroupY")[1] == 8500
                @test get_header(s, "SourceX")[end] == 4000
                @test get_header(s, "SourceY")[end] == 4500
                @test get_header(s, "GroupX")[end] == 8000
                @test get_header(s, "GroupY")[end] == 9000
	end
end


text:
# Word Problems Logical Reasoning Tests
# http://www.graduatewings.co.uk/word-problems-logical-tests

# Questions 11-17
using AbstractLogic

"""
A cruise line is scheduling seven week-long voyages for the ship Freedom. Each
voyage will occur in exactly one of the first seven weeks of the season:
weeks 1 through 7. Each voyage will be to exactly one of four destinations:
Guadeloupe, Jamaica, Martinique, or Trinidad. Each destination will be scheduled
for at least one of the weeks. The following conditions apply to Freedom’s schedule:
"""

# exactly one of four destinations: Guadeloupe, Jamaica, Martinique, or Trinidad
logicset = logicalparse(
  ["w1, w2, w3, w4, w5, w6, w7  ∈ Guadeloupe, Jamaica, Martinique, Trinidad"])

# Each destination will be scheduled for at least one of the weeks.
logicset = logicalparse(["{{i}} = 'Guadeloupe' {{1,}}"], logicset)
logicset = logicalparse(["{{i}} = 'Jamaica' {{1,}}"], logicset)
logicset = logicalparse(["{{i}} = 'Martinique' {{1,}}"], logicset)
logicset = logicalparse(["{{i}} = 'Trinidad' {{1,}}"], logicset)

# Jamaica will not be its destination in week 4.
logicset = logicalparse(["w4 != 'Jamaica'"], logicset)

# Trinidad will be its destination in week 7.
logicset = logicalparse(["w7 = 'Trinidad'"], logicset)

# Freedom will make exactly two voyages to Martinique,
logicset = logicalparse(["{{i}} = 'Martinique' {{2}}"], logicset)

# and at least one voyage to Guadeloupe will occur in some week between those two voyages.
# Unfortunately this can be a bit tedious.
logicset = logicalparse(["w1,w7 = 'Martinique' ==> w2,w3,w4,w5,w6 |= 'Guadeloupe'"], logicset)
logicset = logicalparse(["w2,w7 = 'Martinique' ==> w3,w4,w5,w6 |= 'Guadeloupe'"], logicset)
logicset = logicalparse(["w3,w7 = 'Martinique' ==> w4,w5,w6 |= 'Guadeloupe'"], logicset)
logicset = logicalparse(["w4,w7 = 'Martinique' ==> w5,w6 |= 'Guadeloupe'"], logicset)
logicset = logicalparse(["w5,w7 = 'Martinique' ==> w6 |= 'Guadeloupe'"], logicset)

logicset = logicalparse(["w1,w6 = 'Martinique' ==> w2,w3,w4,w5 |= 'Guadeloupe'"], logicset)
logicset = logicalparse(["w2,w6 = 'Martinique' ==> w3,w4,w5 |= 'Guadeloupe'"], logicset)
logicset = logicalparse(["w3,w6 = 'Martinique' ==> w4,w5 |= 'Guadeloupe'"], logicset)
logicset = logicalparse(["w4,w6 = 'Martinique' ==> w5 |= 'Guadeloupe'"], logicset)

logicset = logicalparse(["w1,w5 = 'Martinique' ==> w2,w3,w4 |= 'Guadeloupe'"], logicset)
logicset = logicalparse(["w2,w5 = 'Martinique' ==> w3,w4 |= 'Guadeloupe'"], logicset)
logicset = logicalparse(["w3,w5 = 'Martinique' ==> w4 |= 'Guadeloupe'"], logicset)

logicset = logicalparse(["w1,w4 = 'Martinique' ==> w2,w3 |= 'Guadeloupe'"], logicset)
logicset = logicalparse(["w2,w4 = 'Martinique' ==> w3 |= 'Guadeloupe'"], logicset)

logicset = logicalparse(["w1,w3 = 'Martinique' ==> w2 |= 'Guadeloupe'"], logicset)

# Guadeloupe will be its destination in the week preceding any voyage it makes to Jamaica.
logicset = logicalparse(["{{i}} = 'Jamaica' ==> {{i-1}} = 'Guadeloupe'"], logicset)
logicset = logicalparse(["w1 != 'Jamaica'"], logicset)

# No destination will be scheduled for consecutive weeks.
logicset = logicalparse(["{{i}} != {{i+1}}"], logicset)

######
# 11. Which one of the following is an acceptable schedule of destinations for
# Freedom, in order from week 1 through week 7?

# (A) Guadeloupe, Jamaica, Martinique, Trinidad, Guadeloupe, Martinique, Trinidad
# I'll just check the first 4 weeks at first
checkfeasible(
  "w1 = 'Guadeloupe'; w2 = 'Jamaica'; w3 = 'Martinique'; w4 = 'Trinidad';
   w5 = 'Guadeloupe'; w6 = 'Martinique'; w7 = 'Trinidad'", logicset)
# true, 1 out of 50 possible combinations 'true'.

# (B) Guadeloupe, Martinique, Trinidad, Martinique, Guadeloupe, Jamaica, Trinidad
checkfeasible("w1 = 'Guadeloupe'; w2 = 'Martinique'; w3 = 'Trinidad'; w4 = 'Martinique';
 w5 = 'Guadeloupe'; w6 = 'Jamaica'; w7 = 'Trinidad'", logicset)
# false, 0 out of 50 possible combinations 'true'.

# (C) Jamaica, Martinique, Guadeloupe, Martinique, Guadeloupe, Jamaica, Trinidad
checkfeasible("w1 = 'Jamaica'; w2 = 'Martinique'; w3 = 'Guadeloupe'; w4 = 'Martinique';
 w5 = 'Guadeloupe'; w6 = 'Jamaica'; w7 = 'Trinidad'", logicset)
# false, 0 out of 50 possible combinations 'true'.

# (D) Martinique, Trinidad, Guadeloupe, Jamaica, Martinique, Guadeloupe, Trinidad
checkfeasible("w1 = 'Martinique'; w2 = 'Trinidad'; w3 = 'Guadeloupe'; w4 = 'Jamaica';
 w5 = 'Martinique'; w6 = 'Guadeloupe'; w7 = 'Trinidad'", logicset)
# false, 0 out of 50 possible combinations 'true'.

# (E) Martinique, Trinidad, Guadeloupe, Trinidad, Guadeloupe, Jamaica, Martinique
checkfeasible("w1 = 'Martinique'; w2 = 'Trinidad'; w3 = 'Guadeloupe'; w4 = 'Trinidad';
 w5 = 'Guadeloupe'; w6 = 'Jamaica'; w7 = 'Martinique'", logicset)
# false, 0 out of 50 possible combinations 'true'.

######
# 12. Which one of the following CANNOT be true about Freedom’s schedule of voyages?
# (A) Freedom makes a voyage to Trinidad in week 6
checkfeasible("w6 = 'Trinidad'", logicset)
# false, 0 out of 50 possible combinations 'true'.  - Answer

# (B) Freedom makes a voyage to Martinique in week 5.
checkfeasible("w5 = 'Martinique'", logicset)
# possible,  11 out of 25 possible combinations 'true'.

# (C) Freedom makes a voyage to Jamaica in week 6.
checkfeasible("w6 = 'Jamaica'", logicset)
# possible,  3 out of 25 possible combinations 'true'.

# (D) Freedom makes a voyage to Jamaica in week 3.
checkfeasible("w3 = 'Jamaica'", logicset)
# possible,  8 out of 25 possible combinations 'true'.

# (E) Freedom makes a voyage to Guadeloupe in week 3
checkfeasible("w3 = 'Guadeloupe'", logicset)
# possible,  4 out of 25 possible combinations 'true'.

######
# 13. If Freedom makes a voyage to Trinidad in week 5, which one of the following
# could be true?
logicsetreduced = logicalparse(["w5 = 'Trinidad'"], logicset)

#  (A) Freedom makes a voyage to Trinidad in week 1.
checkfeasible("w1 = 'Trinidad'", logicsetreduced)[2]
#false, 0 out of 3 possible combinations 'true'.

# (B) Freedom makes a voyage to Martinique in week 2.
checkfeasible("w2 = 'Martinique'", logicsetreduced)
# false, 0 out of 3 possible combinations 'true'.

# (C) Freedom makes a voyage to Guadeloupe in week 3.
checkfeasible("w3 = 'Guadeloupe'", logicsetreduced)
# false, 0 out of 3 possible combinations 'true'.

# (D) Freedom makes a voyage to Martinique in week 4.
checkfeasible("w4 = 'Martinique'", logicsetreduced)
# possible,  1 out of 3 possible combinations 'true'.

# (E) Freedom makes a voyage to Jamaica in week 6.
checkfeasible("w6 = 'Jamaica'", logicsetreduced)
# false, 0 out of 3 possible combinations 'true'.

######
# 14. If Freedom makes a voyage to Guadeloupe in week 1 and a voyage to Jamaica in
# week 5, which one of the following must be true?
logicsetreduced = logicalparse(["w1 = 'Guadeloupe'; w5 = 'Jamaica'"], logicset)

# (A) Freedom makes a voyage to Jamaica in week 2.
checkfeasible("w2 = 'Jamaica'", logicsetreduced)
# possible,  1 out of 3 possible combinations 'true'.

# (B) Freedom makes a voyage to Trinidad in week 2.
checkfeasible("w2 = 'Trinidad'", logicsetreduced)
# possible,  1 out of 3 possible combinations 'true'.

# (C) Freedom makes a voyage to Martinique in week 3.
checkfeasible("w3 = 'Martinique'", logicsetreduced)
# possible,  2 out of 3 possible combinations 'true'.

# (D) Freedom makes a voyage to Guadeloupe in week 6.
checkfeasible("w6 = 'Guadeloupe'", logicsetreduced)
# false, 0 out of 3 possible combinations 'true'.

# (E) Freedom makes a voyage to Martinique in week 6.
checkfeasible("w6 = 'Martinique'", logicsetreduced)
#true, 3 out of 3 possible combinations 'true'.

######
# 15. If Freedom makes a voyage to Guadeloupe in week 1 and to Trinidad in
# week 2, which one of the following must be true?
logicsetreduced = logicalparse(["w1 = 'Guadeloupe'; w2 = 'Trinidad'"], logicset)

# (A) Freedom makes a voyage to Martinique in week 3.
checkfeasible("w3 = 'Martinique'", logicsetreduced, force=true)
# true, 1 out of 1 possible combinations 'true'.

# (B) Freedom makes a voyage to Martinique in week 4.
checkfeasible("w4 = 'Martinique'", logicsetreduced, force=true)
# false, 0 out of 1 possible combinations 'true'.

# (C) Freedom makes a voyage to Martinique in week 5.
checkfeasible("w5 = 'Martinique'", logicsetreduced, force=true)
# false, 0 out of 1 possible combinations 'true'.

# (D) Freedom makes a voyage to Guadeloupe in week 3.
checkfeasible("w3 = 'Guadeloupe'", logicsetreduced, force=true)
# false, 0 out of 1 possible combinations 'true'.

# (E) Freedom makes a voyage to Guadeloupe in week 5
checkfeasible("w5 = 'Guadeloupe'", logicsetreduced, force=true)
# false, 0 out of 1 possible combinations 'true'.

######
# 16. If Freedom makes a voyage to Martinique in week 3, which one of the
# following could be an accurate list of Freedom’s destinations in week 4
# and week 5,respectively?
logicsetreduced = logicalparse(["w3 = 'Martinique'"], logicset)

# (A) Guadeloupe, Trinidad
checkfeasible("w4 = 'Guadeloupe'; w5 = 'Trinidad'", logicsetreduced, any=true)
# true,  1 out of 8 possible combinations 'true'.

# (B) Jamaica, Guadeloupe
checkfeasible("w4 = 'Jamaica'; w5 = 'Guadeloupe'", logicsetreduced, any=true)
# false, 0 out of 8 possible combinations 'true'.

# (C) Martinique, Trinidad
checkfeasible("w4 = 'Martinique'; w5 = 'Trinidad'", logicsetreduced, any=true)
# false, 0 out of 8 possible combinations 'true'.

# (D) Trinidad, Jamaica
checkfeasible("w4 = 'Trinidad'; w5 = 'Jamaica'", logicsetreduced, any=true)
# false, 0 out of 8 possible combinations 'true'.

# (E) Trinidad, Martinique
checkfeasible("w4 = 'Trinidad'; w5 = 'Martinique'", logicsetreduced, any=true)
# false, 0 out of 8 possible combinations 'true'.

######
# 17. Which one of the following must be true about Freedom’s schedule of voyages?

# (A) Freedom makes a voyage to Guadeloupe either in week 1 or else in week 2.
checkfeasible("w1|w2 = 'Guadeloupe'", logicset, force=true)
# false, 22 out of 25 possible combinations 'true'.
logicalparse(["w1,w2 != 'Guadeloupe'"], logicset)

# (B) Freedom makes a voyage to Martinique either in week 2 or else in week 3.
checkfeasible("w2|w3 = 'Martinique'", logicset, force=true)
# false, 12 out of 25 possible combinations 'true'.

# (C) Freedom makes at most two voyages to Guadeloupe.
checkfeasible("{{i}} = 'Guadeloupe' {{,2}}", logicset, force=true)
# false, 20 out of 25 possible combinations 'true'.

# (D) Freedom makes at most two voyages to Jamaica.
checkfeasible("{{i}} = 'Jamaica' {{,2}}", logicset, force=true)
# true, 25 out of 25 possible combinations 'true'.

# (E) Freedom makes at most two voyages to Trinidad.
checkfeasible("{{i}} = 'Trinidad' {{,2}}", logicset, force=true)
# false, 24 out of 25 possible combinations 'true'.


text:
using Documenter
using MultiplesOfPi

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

makedocs(;
    modules=[MultiplesOfPi],
    authors="Jishnu Bhattacharya",
    repo="https://github.com/jishnub/MultiplesOfPi.jl/blob/{commit}{path}#L{line}",
    sitename="MultiplesOfPi.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jishnub.github.io/MultiplesOfPi.jl",
        assets=String[],
    ),
    pages=[
        "Reference" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/jishnub/MultiplesOfPi.jl",
)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Qt5QuickControls2"
version = v"5.15.3"

# Collection of sources required to build qt5
sources = [
    ArchiveSource("https://download.qt.io/official_releases/qt/$(version.major).$(version.minor)/$version/submodules/qtquickcontrols2-everywhere-opensource-src-$version.tar.xz",
                  "3b95e9be58c3d918f4c4a805116459a403ae6c61e33edf71b319672c46af562f"),
]

script = raw"""
# We need to make sure qmake is a copy and not a symlink for the Qt installation paths to be correct:
prefixqmake=$(which qmake)
realqmake=$(readlink $prefixqmake)
rm $prefixqmake
cp $bindir/$realqmake $prefixqmake

case "$target" in
    *apple-darwin*)
        ln -s $libdir/QtCore.framework/Headers $includedir/QtCore
        export QT_MAC_SDK_NO_VERSION_CHECK=1
        ;;
esac

cd $WORKSPACE/srcdir
mkdir build
cd build/

qtsrcdir=`ls -d ../qtquickcontrols2-*`

qmake $qtsrcdir
make -j${nproc}

# It appends the root path here with the qmake dir, which is /workspace/artifacts/blah...
make INSTALL_ROOT=$PWD/dummyroot install
rsync -ua dummyroot/workspace/*/destdir/* $prefix

install_license $qtsrcdir/LICENSE.LGPL3

case "$target" in
    *apple-darwin*)
        rm $includedir/QtCore
        ;;
esac
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms_linux = [
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("armv7l", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="glibc"),
    Platform("x86_64", "freebsd"),
    Platform("powerpc64le", "linux"; libc="glibc"),
]
platforms_linux = expand_cxxstring_abis(platforms_linux)
platforms_win = expand_cxxstring_abis([Platform("x86_64", "windows"), Platform("i686", "windows")])
platforms_macos = [ Platform("x86_64", "macos"), Platform("aarch64", "macos") ]

# The products that we will ensure are always built
products = [
    FileProduct(["qml/QtQuick/Controls.2"], :quickcontrols2),
    LibraryProduct(["Qt5QuickControls2", "libQt5QuickControls2", "QtQuickControls2"], :libqt5quickcontrols2),
    LibraryProduct(["Qt5QuickTemplates2", "libQt5QuickTemplates2", "QtQuickTemplates2"], :libqt5quicktemplates2),
]

products_macos = [
    FileProduct(["qml/QtQuick/Controls.2"], :quickcontrols2),
    FrameworkProduct("QtQuickControls2", :libqt5quickcontrols2),
    FrameworkProduct("QtQuickTemplates2", :libqt5quicktemplates2),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Qt5Declarative_jll"),
]

include("../../fancy_toys.jl")

julia_compat = "1.6"

if any(should_build_platform.(triplet.(platforms_linux)))
    build_tarballs(ARGS, name, version, sources, script, platforms_linux, products, dependencies; preferred_gcc_version = v"7", julia_compat)
end
if any(should_build_platform.(triplet.(platforms_win)))
    build_tarballs(ARGS, name, version, sources, script, platforms_win, products, dependencies; preferred_gcc_version = v"8", julia_compat)
end
if any(should_build_platform.(triplet.(platforms_macos)))
    build_tarballs(ARGS, name, version, sources, script, platforms_macos, products_macos, dependencies; preferred_gcc_version = v"7", julia_compat)
end

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
# This file was generated, do not modify it. # hide
x = rand(10)
println(sum(x))
y = 5

max_stars_repo_name: mcognetta/franklindocs

max_issues_repo_name: mcognetta/franklindocs

max_forks_repo_name: mcognetta/franklindocs


text:
# Needs CXX, Julia-0.5
ENV["PATH"]="/usr/local/cuda-7.5/bin:"*ENV["PATH"]
ENV["LD_LIBRARY_PATH"]="/usr/local/cuda-7.5/lib64:/usr/local/cuda/nvvm/lib64:"*get(ENV,"LD_LIBRARY_PATH","")
ENV["AFMODE"] = "CUDA"  # this example works with CUDA
using ArrayFire, CUDArt, GLAbstraction, Cxx, ModernGL

type CUDAGLBuffer{T} <: GPUArray{T, 1}
    buffer::GLBuffer{T}
    graphics_resource::Ref{CUDArt.rt.cudaGraphicsResource_t}
    ismapped::Bool
end

function CUDAGLBuffer(buffer::GLBuffer, flag = 0)
    cuda_resource = Ref{CUDArt.rt.cudaGraphicsResource_t}(C_NULL)
    CUDArt.rt.cudaGraphicsGLRegisterBuffer(cuda_resource, buffer.id, flag)
    CUDAGLBuffer(buffer, cuda_resource, false)
end
function map_resource(buffer::CUDAGLBuffer)
    if !buffer.ismapped
        CUDArt.rt.cudaGraphicsMapResources(1, buffer.graphics_resource, C_NULL)
        buffer.ismapped = true;
    end
    nothing
end

function unmap_resource(buffer::CUDAGLBuffer)
    if buffer.ismapped
        CUDArt.rt.cudaGraphicsUnmapResources(1, buffer.graphics_resource, C_NULL)
        buffer.ismapped = false
    end
    nothing
end

function copy_from_device_pointer{T}(
        cuda_mem_ptr::Ptr{T},
        cuda_gl_buffer::CUDAGLBuffer,
    )
    map_resource(cuda_gl_buffer)
    buffersize = length(cuda_gl_buffer.buffer)*sizeof(eltype(cuda_gl_buffer.buffer))
    if cuda_gl_buffer.buffer.buffertype == GL_RENDERBUFFER
        array_ptr = Ref{CUDArt.rt.cudaArray_t}(C_NULL)
        CUDArt.rt.cudaGraphicsSubResourceGetMappedArray(array_ptr, cuda_gl_buffer.graphics_resource[], 0, 0)
        CUDArt.rt.cudaMemcpyToArray(array_ptr[], 0, 0, cuda_mem_ptr, buffersize, CUDArt.rt.cudaMemcpyDeviceToDevice)
    else
        opengl_ptr = Ref{Ptr{Void}}(C_NULL); size_ref = Ref{Csize_t}(buffersize)
        CUDArt.rt.cudaGraphicsResourceGetMappedPointer(opengl_ptr, size_ref, cuda_gl_buffer.graphics_resource[])
        CUDArt.rt.cudaMemcpy(opengl_ptr[], cuda_mem_ptr, buffersize, CUDArt.rt.cudaMemcpyDeviceToDevice)
    end
    unmap_resource(cuda_gl_buffer)
end

"""
 Gets the device pointer from the mapped resource
 Sets is_mapped to true
"""
function copy_to_device_pointer{T}(
        cuda_mem_ptr::Ptr{T},
        cuda_gl_buffer::CUDAGLBuffer,
    )
    map_resource(cuda_gl_buffer)
    is_mapped = true
    buffersize = length(cuda_gl_buffer.buffer)*sizeof(eltype(cuda_gl_buffer.buffer))
    if cuda_gl_buffer.buffer.buffertype == GL_RENDERBUFFER
        array_ptr = Ref{CUDArt.rt.cudaArray_t}(C_NULL);
        CUDArt.rt.cudaGraphicsSubResourceGetMappedArray(array_ptr, cuda_gl_buffer.graphics_resource[], 0, 0)
        CUDArt.rt.cudaMemcpyFromArray(cuda_mem_ptr, array_ptr[], 0, 0, buffersize, CUDArt.rt.cudaMemcpyDeviceToDevice)
    else
        opengl_ptr = Ref{Ptr{Void}}(C_NULL); size_ref = Ref{Csize_t}(buffersize)
        CUDArt.rt.cudaGraphicsResourceGetMappedPointer(opengl_ptr, size_ref, cuda_gl_buffer.graphics_resource[])
        CUDArt.rt.cudaMemcpy(cuda_mem_ptr, opengl_ptr, buffersize, CUDArt.rt.cudaMemcpyDeviceToDevice)
    end
    unmap_resource(cuda_gl_buffer)
end

# ArrayFire.AFArray
function Base.copy(source::ArrayFire.AFArray, target::CUDAGLBuffer)
    d_ptr = ArrayFire.af_device(source)
    copy_from_device_pointer(d_ptr, target)
end
function Base.copy(source::CUDAGLBuffer, target::ArrayFire.AFArray)
    d_ptr = ArrayFire.af_device(target)
    copy_to_device_pointer(d_ptr, target)
end

using GLVisualize, GeometryTypes, Colors
w=glscreen();@async renderloop(w)
cat = loadasset("cat.obj")
colors = RGBA{Float32}[RGBA{Float32}(rand(), rand(), rand(), 1) for i=1:length(vertices(cat))]
catmesh = GLNormalVertexcolorMesh(
    vertices=vertices(cat), faces=faces(cat), normals=normals(cat),
    color=colors
)
_view(visualize(catmesh))

vertsvec = reinterpret(Float32, vertices(cat), (3,length(vertices(cat))));
colorsvec = reinterpret(Float32, colors, (4,length(colors)));
af_vertices = AFArray(vertsvec)
af_colors = AFArray(colorsvec)

gl_verts = w.renderlist[1][1][:vertices]
gl_colors = w.renderlist[1][1][:color]

cu_gl_verts = CUDAGLBuffer(gl_verts)
cu_gl_colors = CUDAGLBuffer(gl_colors)
b = af_vertices .* 2f0
copy(b, cu_gl_verts)


text:
#= TODO:
- How to deal with f_increased (as state and preconstate shares the same x and x_previous vectors)
- Check whether this makes sense for other preconditioners than GradientDescent and L-BFGS
* There might be some issue of dealing with x_current and x_previous in MomentumGradientDescent
* Trust region based methods may not work because we assume the preconditioner calls perform_linesearch!
=#

abstract type AbstractNGMRES <: FirstOrderOptimizer end

# TODO: Enforce TPrec <: Union{FirstOrderoptimizer,SecondOrderOptimizer}?
struct NGMRES{IL, Tp,TPrec <: AbstractOptimizer,L} <: AbstractNGMRES
    alphaguess!::IL       # Initial step length guess for linesearch along direction xP->xA
    linesearch!::L        # Preconditioner moving from xP to xA (precondition x to accelerated x)
    manifold::Manifold
    nlprecon::TPrec       # Nonlinear preconditioner
    nlpreconopts::Options # Preconditioner options
    ϵ0::Tp                # Ensure A-matrix is positive definite
    wmax::Int             # Maximum window size
end

struct OACCEL{IL, Tp,TPrec <: AbstractOptimizer,L} <: AbstractNGMRES
    alphaguess!::IL       # Initial step length guess for linesearch along direction xP->xA
    linesearch!::L        # Linesearch between xP and xA (precondition x to accelerated x)
    manifold::Manifold
    nlprecon::TPrec       # Nonlinear preconditioner
    nlpreconopts::Options # Preconditioner options
    ϵ0::Tp                # Ensure A-matrix is positive definite
    wmax::Int             # Maximum window size
end

Base.summary(s::NGMRES) = "Nonlinear GMRES preconditioned with $(summary(s.nlprecon))"
Base.summary(s::OACCEL) = "O-ACCEL preconditioned with $(summary(s.nlprecon))"

"""
# N-GMRES
## Constructor
```julia
NGMRES(;
        alphaguess = LineSearches.InitialStatic(),
        linesearch = LineSearches.HagerZhang(),
        manifold = Flat(),
        wmax::Int = 10,
        ϵ0 = 1e-12,
        nlprecon = GradientDescent(
            alphaguess = LineSearches.InitialStatic(alpha=1e-4,scaled=true),
            linesearch = LineSearches.Static(),
            manifold = manifold),
        nlpreconopts = Options(iterations = 1, allow_f_increases = true),
      )
```

## Description
This algorithm takes a step given by the nonlinear preconditioner `nlprecon`
and proposes an accelerated step by minimizing an approximation of
the (\ell_2) residual of the gradient on a subspace spanned by the previous
`wmax` iterates.

N-GMRES was originally developed for solving nonlinear systems [1], and reduces to
GMRES for linear problems.
Application of the algorithm to optimization is covered, for example, in [2].

## References
[1] De Sterck. Steepest descent preconditioning for nonlinear GMRES optimization. NLAA, 2013.
[2] Washio and Oosterlee. Krylov subspace acceleration for nonlinear multigrid schemes. ETNA, 1997.
"""
function NGMRES(;manifold::Manifold = Flat(),
                alphaguess = LineSearches.InitialStatic(),
                linesearch = LineSearches.HagerZhang(),
                nlprecon = GradientDescent(
                    alphaguess = LineSearches.InitialStatic(alpha=1e-4,scaled=true), # Step length arbitrary,
                    linesearch = LineSearches.Static(),
                    manifold = manifold),
                nlpreconopts = Options(iterations = 1, allow_f_increases = true),
                ϵ0 = 1e-12, # ϵ0 = 1e-12  -- number was an arbitrary choice#
                wmax::Int = 10) # wmax = 10  -- number was an arbitrary choice to match L-BFGS field `m`
    @assert manifold == nlprecon.manifold
    NGMRES(alphaguess, linesearch, manifold, nlprecon, nlpreconopts, ϵ0, wmax)
end

"""
# O-ACCEL
## Constructor
```julia
OACCEL(;manifold::Manifold = Flat(),
       alphaguess = LineSearches.InitialStatic(),
       linesearch = LineSearches.HagerZhang(),
       nlprecon = GradientDescent(
           alphaguess = LineSearches.InitialStatic(alpha=1e-4,scaled=true),
           linesearch = LineSearches.Static(),
           manifold = manifold),
       nlpreconopts = Options(iterations = 1, allow_f_increases = true),
       ϵ0 = 1e-12,
       wmax::Int = 10)
```

## Description
This algorithm takes a step given by the nonlinear preconditioner `nlprecon`
and proposes an accelerated step by minimizing an approximation of
the objective on a subspace spanned by the previous
`wmax` iterates.

O-ACCEL is a slight tweak of N-GMRES, first presented in [1].

## References
[1] Riseth. Objective acceleration for unconstrained optimization. 2018.
"""
function OACCEL(;manifold::Manifold = Flat(),
                alphaguess = LineSearches.InitialStatic(),
                linesearch = LineSearches.HagerZhang(),
                nlprecon = GradientDescent(
                    alphaguess = LineSearches.InitialStatic(alpha=1e-4,scaled=true), # Step length arbitrary
                    linesearch = LineSearches.Static(),
                    manifold = manifold),
                nlpreconopts = Options(iterations = 1, allow_f_increases = true),
                ϵ0 = 1e-12, # ϵ0 = 1e-12  -- number was an arbitrary choice
                wmax::Int = 10) # wmax = 10  -- number was an arbitrary choice to match L-BFGS field `m`
    @assert manifold == nlprecon.manifold
    OACCEL(alphaguess, linesearch, manifold, nlprecon, nlpreconopts, ϵ0, wmax)
end

mutable struct NGMRESState{P,Tx,Te,T,eTx} <: AbstractOptimizerState where P <: AbstractOptimizerState
    # eTx is the eltype of Tx
    x::Tx                    # Reference to nlpreconstate.x
    x_previous::Tx           # Reference to nlpreconstate.x_previous
    x_previous_0::Tx         # Used to deal with assess_convergence of NGMRES
    f_x_previous::T
    f_x_previous_0::T        # Used to deal with assess_convergence of NGMRES
    f_xP::T                  # For tracing purposes
    grnorm_xP::T             # For tracing purposes
    s::Tx                    # Search direction for linesearch between xP and xA
    nlpreconstate::P         # Nonlinear preconditioner state
    X::Array{eTx,2}          # Solution vectors in the window
    R::Array{eTx,2}          # Gradient vectors in the window
    Q::Array{T,2}            # Storage to create linear system (TODO: make Symmetric?)
    ξ::Te                    # Storage to create linear system
    curw::Int                # Counter for current window size
    A::Array{T,2}            # Container for Aα = b
    b::Vector{T}             # Container for Aα = b
    xA::Vector{eTx}          # Container for accelerated step
    k::Int                   # Used for indexing where to put values in the Storage containers
    restart::Bool            # Restart flag
    g_tol::T                 # Exit tolerance to be checked after nonlinear preconditioner apply
    subspacealpha::Vector{T} # Storage for coefficients in the subspace for the acceleration step
    @add_linesearch_fields()
end

"Update storage Q[i,j] and Q[j,i] for `NGMRES`"
@inline function _updateQ!(Q, i::Int, j::Int, X, R, ::NGMRES)
    Q[j,i] = real(dot(R[:, j], R[:,i]))
    if i != j
        Q[i,j] = Q[j, i]          # TODO: Use Symmetric?
    end
end

"Update storage A[i,j] for `NGMRES`"
@inline function _updateA!(A, i::Int, j::Int, Q, ξ, η, ::NGMRES)
    A[i,j] = Q[i,j]-ξ[i]-ξ[j]+η
end

"Update storage ξ[i,:] for `NGMRES`"
@inline function _updateξ!(ξ, i::Int, X, x, R, r, ::NGMRES)
    ξ[i] = real(dot(vec(r), R[:,i]))
end

"Update storage b[i] for `NGMRES`"
@inline function _updateb!(b, i::Int, ξ, η, ::NGMRES)
    b[i] = η - ξ[i]
end

"Update value η for `NGMRES`"
@inline function _updateη(x, r, ::NGMRES)
    real(dot(r, r))
end

"Update storage Q[i,j] and Q[j,i] for `OACCEL`"
@inline function _updateQ!(Q, i::Int, j::Int, X, R, ::OACCEL)
    Q[i,j] = real(dot(X[:,i], R[:,j]))
    if i != j
        Q[j,i] = real(dot(X[:,j], R[:,i]))
    end
end

"Update storage A[i,j] for `OACCEL`"
@inline function _updateA!(A, i::Int, j::Int, Q, ξ, η, ::OACCEL)
    A[i,j] = Q[i,j]-ξ[i,1]-ξ[j,2]+η
end

"Update storage ξ[i,:] for `OACCEL`"
@inline function _updateξ!(ξ, i::Int, X, x, R, r, ::OACCEL)
    ξ[i,1] = real(dot(X[:,i], r))
    ξ[i,2] = real(dot(x, R[:,i]))
end

"Update storage b[i] for `OACCEL`"
@inline function _updateb!(b, i::Int, ξ, η, ::OACCEL)
    b[i] = η - ξ[i,1]
end

"Update value η for `OACCEL`"
@inline function _updateη(x, r, ::OACCEL)
    real(dot(x, r))
end

const ngmres_oaccel_warned = Ref{Bool}(false)
function initial_state(method::AbstractNGMRES, options, d, initial_x::AbstractArray{eTx}) where eTx
    if !(typeof(method.nlprecon) <: Union{GradientDescent,LBFGS})
        if !ngmres_oaccel_warned[]
            @warn "Use caution. N-GMRES/O-ACCEL has only been tested with Gradient Descent and L-BFGS preconditioning."
            ngmres_oaccel_warned[] = true
        end
    end

    nlpreconstate = initial_state(method.nlprecon, method.nlpreconopts, d, initial_x)
    # Manifold comment:
    # We assume nlprecon calls retract! and project_tangent! on
    # nlpreconstate.x and gradient(d)
    T = real(eTx)

    n = length(nlpreconstate.x)
    wmax = method.wmax
    X = Array{eTx}(undef, n, wmax)
    R = Array{eTx}(undef, n, wmax)
    Q = Array{T}(undef, wmax, wmax)
    ξ = if typeof(method) <: OACCEL
        Array{T}(undef, wmax, 2)
    else
        Array{T}(undef, wmax)
    end

    copyto!(view(X,:,1), nlpreconstate.x)
    copyto!(view(R,:,1), gradient(d))

    _updateQ!(Q, 1, 1, X, R, method)

    NGMRESState(nlpreconstate.x,          # Maintain current state in state.x. Use same vector as preconditioner.
                nlpreconstate.x_previous, # Maintain  in state.x_previous. Use same vector as preconditioner.
                similar(nlpreconstate.x), # Maintain state at the beginning of an iteration in state.x_previous_0. Used for convergence asessment.
                T(NaN),                   # Store previous f in state.f_x_previous
                T(NaN),                   # Store f value from the beginning of an iteration in state.f_x_previous_0. Used for convergence asessment.
                T(NaN),                   # Store value f_xP of f(x^P) for tracing purposes
                T(NaN),                   # Store value grnorm_xP of |g(x^P)| for tracing purposes
                similar(initial_x),       # Maintain current search direction in state.s
                nlpreconstate,            # State storage for preconditioner
                X,
                R,
                Q,
                ξ,
                1,                        # curw
                Array{T}(undef, wmax, wmax),     # A
                Array{T}(undef, wmax),           # b
                vec(similar(initial_x)),  # xA
                0,                        # iteration counter
                false,                    # Restart flag
                options.g_tol,            # Exit tolerance check after nonlinear preconditioner apply
                Array{T}(undef, wmax),           # subspacealpha
                @initial_linesearch()...)
end

nlprecon_post_optimize!(d, state, method) = update_h!(d, state.nlpreconstate, method)

nlprecon_post_accelerate!(d, state, method) = update_h!(d, state.nlpreconstate, method)

function nlprecon_post_accelerate!(d, state::NGMRESState{X,T},
                                   method::LBFGS)  where X where T
    state.nlpreconstate.pseudo_iteration += 1
    update_h!(d, state.nlpreconstate, method)
end

function update_state!(d, state::NGMRESState{X,T}, method::AbstractNGMRES) where X where T
    # Maintain a record of previous position, for convergence assessment
    copyto!(state.x_previous_0, state.x)
    state.f_x_previous_0 = value(d)

    state.k += 1
    curw = state.curw

    # Step 1: Call preconditioner to get x^P
    res = optimize(d, state.x, method.nlprecon, method.nlpreconopts, state.nlpreconstate)
    # TODO: Is project_tangent! necessary, or is it called by nlprecon before exit?
    project_tangent!(method.manifold, gradient(d), state.x)

    if any(.!isfinite.(state.x)) || any(.!isfinite.(gradient(d))) || !isfinite(value(d))
        @warn("Non-finite values attained from preconditioner $(summary(method.nlprecon)).")
        return true
    end

    # Calling value_gradient! in normally done on state.x in optimize or update_g! above,
    # but there are corner cases where we need this.
    state.f_xP, _g = value_gradient!(d, state.x)
    # Manifold start
    project_tangent!(method.manifold, gradient(d), state.x)
    # Manifold stop
    gP = gradient(d)
    state.grnorm_xP = g_residual(gP)

    if g_residual(gP) ≤ state.g_tol
        return false # Exit on gradient norm convergence
    end

    # Deals with update_h! etc for preconditioner, if needed
    nlprecon_post_optimize!(d, state, method.nlprecon)

    # Step 2: Do acceleration calculation
    η = _updateη(state.x, gP, method)

    for i = 1:curw
        # Update storage vectors according to method {NGMRES, OACCEL}
        _updateξ!(state.ξ, i, state.X, state.x, state.R, gP, method)
        _updateb!(state.b, i, state.ξ, η, method)
    end

    for i = 1:curw
        for j = 1:curw
            # Update system matrix according to method {NGMRES, OACCEL}
            _updateA!(state.A, i, j, state.Q, state.ξ, η, method)
        end
    end

    α = view(state.subspacealpha, 1:curw)
    Aview = view(state.A, 1:curw, 1:curw)
    bview = view(state.b, 1:curw)
    # The outer max is to avoid δ=0, which may occur if A=0, e.g. at numerical convergence
    δ = method.ϵ0*max(maximum(diag(Aview)), method.ϵ0)
    try
        α .= (Aview + δ*I) \ bview
    catch e
        @warn("Calculating α failed in $(summary(method)).")
        @warn("Exception info:\n $e")
        α .= NaN
    end
    if any(isnan, α)
        @warn("Calculated α is NaN in $(summary(method)). Restarting ...")
        state.s .= zero(eltype(state.s))
        state.restart = true
    else
        # xA = xP + \sum_{j=1}^{curw} α[j] * (X[j] - xP)
        state.xA .= (1.0-sum(α)).*vec(state.x) .+
            sum(state.X[:,k]*α[k] for k = 1:curw)

        state.s .=  reshape(state.xA, size(state.x)) .- state.x
    end

    # 3: Perform condition checks
    if real(dot(state.s, gP)) ≥ 0 || !isfinite(real(dot(state.s, gP)))
        # Moving from xP to xA is *not* a descent direction
        # Discard xA
        state.restart = true # TODO: expand restart heuristics
        lssuccess = true
        state.alpha = 0.0
    else
        state.restart = false

        # Update f_x_previous and dphi_0_previous according to preconditioner step
        # This may be used in perform_linesearch!/alphaguess! when moving from x^P to x^A
        # TODO: make this a function?
        state.f_x_previous = state.nlpreconstate.f_x_previous
        if typeof(method.alphaguess!) <: LineSearches.InitialConstantChange
            nlprec = method.nlprecon
            if isdefined(nlprec, :alphaguess!) &&
                typeof(nlprec.alphaguess!) <: LineSearches.InitialConstantChange
                method.alphaguess!.dϕ_0_previous[] = nlprec.alphaguess!.dϕ_0_previous[]
            end
        end
        # state.x_previous and state.x are dealt with by reference

        lssuccess = perform_linesearch!(state, method, ManifoldObjective(method.manifold, d))
        @. state.x = state.x + state.alpha * state.s
        # Manifold start
        retract!(method.manifold, state.x)
        # Manifold stop

        # TODO: Move these into `nlprecon_post_accelerate!` ?
        state.nlpreconstate.f_x_previous = state.f_x_previous
        if typeof(method.alphaguess!) <: LineSearches.InitialConstantChange
            nlprec = method.nlprecon
            if isdefined(nlprec, :alphaguess!) &&
                typeof(nlprec.alphaguess!) <: LineSearches.InitialConstantChange
                nlprec.alphaguess!.dϕ_0_previous[] = method.alphaguess!.dϕ_0_previous[]
            end
        end
        # Deals with update_h! etc. for preconditioner, if needed
        nlprecon_post_accelerate!(d, state, method.nlprecon)
    end
    #=
    Update x_previous and f_x_previous to be the values at the beginning
    of the N-GMRES iteration. For convergence assessment purposes.
    =#
    copyto!(state.x_previous, state.x_previous_0)
    state.f_x_previous = state.f_x_previous_0

    lssuccess == false # Break on linesearch error
end

function update_g!(d, state, method::AbstractNGMRES)
    # Update the function value and gradient
    # TODO: do we need a retract! on state.x here?
    value_gradient!(d, state.x)
    project_tangent!(method.manifold, gradient(d), state.x)

    if state.restart == false
        state.curw = min(state.curw + 1, method.wmax)
    else
        state.k = 0
        state.curw = 1
    end
    j = mod(state.k, method.wmax) + 1

    copyto!(view(state.X,:,j), vec(state.x))
    copyto!(view(state.R,:,j), vec(gradient(d)))

    for i = 1:state.curw
        _updateQ!(state.Q, i, j, state.X, state.R, method)
    end
end

function trace!(tr, d, state, iteration, method::AbstractNGMRES, options)
    dt = Dict()
    if options.extended_trace
        dt["x"] = copy(state.x)
        dt["g(x)"] = copy(gradient(d))
        dt["subspace-α"] = state.subspacealpha[1:state.curw-1]
        if state.restart == true
            dt["Current step size"] = NaN
        else
            dt["Current step size"] = state.alpha
            # This is a wasteful hack to get the previous values for debugging purposes only.
            xP = state.x .- state.alpha .* state.s
            dt["x^P"] = copy(xP)
            # TODO: What's a good way to include g(x^P) here without messing up gradient counts?
        end
    end
    dt["Restart"] = state.restart
    if state.restart == false
        dt["f(x^P)"] = state.f_xP
        dt["|g(x^P)|"] = state.grnorm_xP
    end

    g_norm = g_residual(d)
    update!(tr,
            iteration,
            value(d),
            g_norm,
            dt,
            options.store_trace,
            options.show_trace,
            options.show_every,
            options.callback)
end

function assess_convergence(state::NGMRESState, d, options)
    default_convergence_assessment(state, d, options)
end

function default_options(method::AbstractNGMRES)
    Dict(:allow_f_increases => true)
end


text:
using FRGRealTime
using Test
using FiniteDifferences
@testset "flowpp" begin
    @test FRGRealTime.flowpp(0.1, 1.0, 5.0, -2.0, 10.0)==central_fdm(5, 1)(k->FRGRealTime.loopfunpp(0.1, 1.0, k, -2.0, 10.0),5.0)
    @test FRGRealTime.flowpp(0.5, 1.0, 5.0, -2.0, 10.0)==central_fdm(5, 1)(k->FRGRealTime.loopfunpp(0.5, 1.0, k, -2.0, 10.0),5.0)
    @test FRGRealTime.flowpp(1.1, 1.0, 5.0, -2.0, 10.0)==central_fdm(5, 1)(k->FRGRealTime.loopfunpp(1.1, 1.0, k, -2.0, 10.0),5.0)
    @test FRGRealTime.flowpp(10.0, 1.0, 5.0, -2.0, 10.0)≈central_fdm(5, 1)(k->FRGRealTime.loopfunpp(10.0, 1.0, k, -2.0, 10.0),5.0)
end


text:
module HH

import BP

#in-place builder of the q x q HH matrix in mom. sp.
momsphhmat!(A::Matrix{Complex{Float64}}, kx0::Float64,ky::Float64) =
   momsphhmat!(A, kx0,ky, 1)

function momsphhmat!(A::Matrix{Complex{Float64}},
                     kx0::Float64,ky::Float64,
                     p::Int)
    q::Int = size(A,1)
    A[1,q] = -exp(-im*q*ky)
    A[q,1] = -exp(im*q*ky)
    #upper subdiagonal
    for i in diagind(A,1)
        A[i] = -one(Complex{Float64})
    end
    #lower subdiagonal
    for i in diagind(A,-1)
        A[i] = -one(Complex{Float64})
    end
    #main diagonal
    α = p/q
    for (j,i) in enumerate(diagind(A))
        A[i] = -2*cos(kx0 + 2*π*α*j)
    end
    nothing
end

# Computes the q energy levels E(p_y)
hhladder(q::Int) = hhladder(1,q)
function hhladder(p::Int, q::Int)
    ky = linspace(-π, π, 100)
    kx₀= linspace(-π/q, π/q, 100)
    E = Array(Float64, (length(ky),length(kx₀),q))
    H = zeros(Complex{Float64}, (q,q))
    for col=1:length(kx₀), row=1:length(ky)
        momsphhmat!(H, kx₀[col], ky[row], p)
        E[row,col,:] = eigvals(H)
    end
    return E
end
hhladder!(E::Array{Float64,3}) = hhladder!(E, 1)
function hhladder!(E::Array{Float64, 3}, p::Int)
    ky = linspace(-π, π, size(E,1))
    q = size(E,3)
    kx₀= linspace(-π/q, π/q, size(E,2))
    H = zeros(Complex{Float64}, (q,q))
    for col=1:length(kx₀), row=1:length(ky)
        momsphhmat!(H, kx₀[col], ky[row], p)
        E[row,col,:] = eigvals(H)
    end
    nothing
end

# ground state energy of the HH hamiltonian
hhgrstate!(ve::Matrix{Float64}, q::Int) = hhgrstate!(ve, 1, q)
function hhgrstate!(ve::Matrix{Float64}, p::Int, q::Int)
    ky = linspace(-π, π, size(ve,1))
    kx₀= linspace(-π/q, π/q, size(ve,2))
    H = zeros(Complex{Float64}, (q,q))
    for col=1:length(kx₀), row=1:length(ky)
        momsphhmat!(H, kx₀[col], ky[row], p)
        ve[row,col] = eigmin(H)
    end
    nothing
end

# zero point energy error
ηzpe(q::Int, κ::Float64) = ηzpe(1, q, κ)
ηzpe(p::Int, q::Int, κ::Float64) = (N=15; A =
   spzeros(Complex{Float64}, N^2,N^2); ηzpe(A, p, q, κ))

function ηzpe(M::SparseMatrixCSC{Complex{Float64},Int}, p::Int,
   q::Int, κ::Float64)
    N::Int = sqrt(size(M,1))
    α::Float64 = p/q
    gs = Array(Float64, 25,25)
    hhgrstate!(gs, p, q)
    e1 = mean(gs)
    et = e1 + 1/2*κ/(2π*α)
    BP.buildham_exact!(M, N,α,κ)
    er = real(eigs(M, nev=1, which=:SR, ritzvec=false)[1][1])
    return 4π*α/κ * (er - et)
end

ηzpe(q::Int, κs::Vector{Float64}) = vec(ηzpe(collect(q), κs))
ηzpe(qs::UnitRange{Int}, κ::Float64) = ηzpe(collect(qs), κ)
ηzpe(qs::Vector{Int}, κ::Float64) = vec(ηzpe(qs, collect(κ)))
ηzpe(qs::UnitRange{Int}, κs::Vector{Float64}) = ηzpe(collect(qs), κs)

function ηzpe(qs::Vector{Int}, κs::Vector{Float64})
    N=15
    A = spzeros(Complex{Float64}, N^2,N^2)
    η = Array(Float64, length(qs), length(κs))
    for col=1:length(κs), row=1:length(qs)
        η[row,col] = ηzpe(A, 1, qs[row], κs[col])
    end
    return η
end

# level error
ηlev(q::Int, κ::Float64) = ηlev(1, q, κ)
ηlev(p::Int, q::Int, κ::Float64) = (N=15; A =
   spzeros(Complex{Float64}, N^2,N^2); ηlev(A, p, q, κ))

function ηlev(M::SparseMatrixCSC{Complex{Float64},Int}, p::Int,
   q::Int, κ::Float64)
    N::Int = sqrt(size(M,1))
    α::Float64 = p/q
    BP.buildham_exact!(M, N,α,κ)
    # get first 4 levels of spectrum of HH+trap
    bilevel = real(eigs(M, nev=4, which=:SR, ritzvec=false)[1])
    #calculate level spacing
    return 2π*α/κ*diff(bilevel[1:2])[1] - 1
end
function ηlev(qs::UnitRange{Int}, κ::Float64, lowβ::Int, highβ::Int)
    N=15 # system size is NxN
    # initialize (sparse) Hamiltonian matrix
    M = spzeros(Complex{Float64}, N^2,N^2)
    diffs = Array(Float64, length(qs))
    for (i,q) in enumerate(qs)
        # construct disipationless HH+trap Hamiltonian in Landau gauge
        # q dependent, need to loop over all q!
        BP.buildham_exact!(M, N,1/q,κ)
        # get first 5 levels of spectrum of M
        levels = real(eigs(M, nev=5, which=:SR, ritzvec=false)[1])
        #calculate level spacing
        diffs[i] = 2π/(q*κ) * (levels[highβ+1] - levels[lowβ+1]) - 1.0
    end
    return diffs
end

ηlev(q::Int, κs::Vector{Float64}) = vec(ηlev([q], κs))
ηlev(qs::UnitRange{Int}, κ::Float64) = ηlev([qs], κ)
ηlev(qs::Vector{Int}, κ::Float64) = vec(ηlev(qs, [κ]))
ηlev(qs::UnitRange{Int}, κs::Vector{Float64}) = ηlev([qs], κs)

function ηlev(qs::Vector{Int}, κs::Vector{Float64})
    N=15
    A = spzeros(Complex{Float64}, N^2,N^2)
    η = Array(Float64, length(qs), length(κs))
    for col=1:length(κs), row=1:length(qs)
        η[row,col] = ηlev(A, 1, qs[row], κs[col])
    end
    return η
end

bwidth(qs::UnitRange{Int}) = bwidth([qs])
bwidth(qs::Vector{Int}) = [bwidth(q) for q in qs]
bwidth(q::Int) = bwidth(1, q)

function bwidth(p::Int, q::Int)
    gstate = Array(Float64, 25, 25)
    hhgrstate!(gstate, p, q)
    a,b = extrema(gstate)
    return b-a
end

bgap(qs::UnitRange{Int}) = bgap([qs])
bgap(qs::Vector{Int}) = [bgap(q) for q in qs]
bgap(q::Int) = bgap(1, q)

function bgap(p::Int, q::Int)
    ladder = Array(Float64, (25,25,q))
    hhladder!(ladder, p)
    e1 = mean(ladder[:,:,1])
    e2 = mean(ladder[:,:,2])
    return e2-e1
end
end #module

max_stars_repo_name: berceanu/topo-photon

max_issues_repo_name: berceanu/topo-photon

max_forks_repo_name: berceanu/topo-photon