text:
module ImageTools
# Copyright © 2014-2015 Vardan Andriasyan, Yauhen Yakimovich, Artur Yakimovich.
#
#  MIT license.
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including without
# limitation the rights to use, copy, modify, merge, publish, distribute,
# sublicense, and/or sell copies of the Software, and to permit persons to whom
# the Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

export Morph
export graythresh,
im2bw,
imshow,
imread,
bwlabel,
jet,
hsv,
label2rgb

include("imagetools/morph.jl")

# import Tk
import FileIO
import Images
import Images: Gray, properties, data, label_components
# import FixedPointNumbers
import ImageView # view images for users of the Julia command-line interface (non Juno or IJulia) https://juliaimages.org/latest/install/#Displaying-images-1
import Colors
import BinDeps

################################################################
function graythresh(img)
  # graythresh - Based on the original Paper: N. Otsu, "A Threshold Selection
  # Method from Gray-Level Histograms" 1979. Calculates a threshold of a grayscale
  # image, which can be used downstream to convert a grayscale image to binary
  # image. One input argument is required - image (expected image format obtained
  # by Images.imread, see Images package  for more information). E.g. to convert
  # an array to image use grayim(array) or colorim(array).

  #Check whether the input image is empty
  if isempty(img)
    error("Image is empty");
  end
  #Check whether the input image is gray

  if !(eltype(img) <: Gray)
    error("Input Image should be grayscale")
  end

  # raw --> rawview in Images # https://github.com/JuliaImages/ImageCore.jl/issues/83

  if isa(Images.rawview(img),Array{UInt16,2})
    # Convert image to 8bit and return it's raw values to compute the histogram
    # https://juliaimages.org/latest/conversions_views/#Views-for-%22converting%22-between-fixed-point-and-raw-representations-1
    image_array = map(Images.BitShift(Uint8,8),Images.rawview(img));
  elseif isa(Images.rawview(img),Array{Uint8,2})
    image_array = Images.rawview(img);
  else
    @warn("Input Image is neither Uint8 or Uint16");
    image_array = Images.rawview(img);
  end
  #image_array = rawview(img)

  expected_number_of_bins = 2^8;

  (range, counts) = hist(reshape(image_array,length(image_array)),expected_number_of_bins);

  #get real number of bins to be used below
  number_of_bins = size(counts)[1];
  #Next the algorithm implementation according to the original paper
  p = counts / sum(counts);
  ω = cumsum(p);
  μ = cumsum(p .* (1:number_of_bins));
  μ_t = μ[end];

  sigma_b_squared = (μ_t * ω - μ).^2 ./ (ω .* (1 - ω));

  # Find the location of the maximum value of sigma_b_squared.
  # The maximum may extend over several bins, so average together the
  # locations.  If maxval is NaN, meaning that sigma_b_squared is all NaN,
  # then return 0.
  maxvalue = maximum(sigma_b_squared);
  if isfinite(maxvalue)
    index = findfirst(sigma_b_squared, maxvalue);

    # Normalize the threshold to the range [0, 1].
    threshold = (index - 1) / (expected_number_of_bins - 1);
  else
    threshold = 0.0;
  end

  return threshold
end
################################################################
function im2bw(img, threshold)
  # im2bw - converts image into a binary image. Input: image (expected image
  # format obtained by Images.imread,  see Images package  for more information)
  # and threshold (floating point number). E.g. to convert an array to image use
  # grayim(array) or colorim(array).

  if !(eltype(img) <: Gray)
    error("Input Image should be grayscale")
  end

  if typeof(threshold) != Float64 && typeof(threshold) != Float32 && typeof(threshold) != Float16
    error("Threshold must be of types Float64, Float32 or Float16");
  end
  if threshold > 1 || threshold < 0
    error("Threshold must be between 0 and 1");
  end

  m,n = size(img);
  blacknWhite = falses(n,m)
  imageData = data(img);

  for i = 1:n
    for j = 1:m
      if imageData[j,i] > threshold
        blacknWhite[i,j] = true;
      end
    end
  end

  return Images.Image(blacknWhite, colorspace = "Binary")
end
################################################################
function imshow(image)
  # wrapper for ImageView.view, need to add convertion of image to array and back
  ImageView.view(image)
end
################################################################
function imread(path)
  # wrapper for FileIO.load with added image retrieval from url functionality
  if (occursin(r"http://.*", path) || occursin(r"https://.*", path) || occursin(r"ftp://.*", path)|| occursin(r"smb://.*", path))
    imageContainer = nothing
    # By default, do silent download.
    url = path
    path = tempname()
    cmd = BinDeps.download_cmd(url, path)
    run(pipe(cmd, devnull, devnull));
    imageContainer = path

  else
    imageContainer = path;
  end
  image = FileIO.load(imageContainer);
  return image
end
################################################################
function bwlabel(inputImage,connectivity)
  # wrapper for Images.label_components

  #define pixel connectivity matrices for appropriate inputs
  if connectivity == 4
    connectivityMatrix = bool([0 1 0;1 1 1; 0 1 0])
  elseif connectivity == 8
    connectivityMatrix = bool([1 1 1;1 1 1; 1 1 1]);
  else
    error("Invalid pixel connectivity. Pixel connectivity can be either 4 or 8");
  end

  labeledImage = label_components(inputImage,connectivityMatrix);

  return labeledImage;
end

################################################################
function jet(numberOfColors::Int64)
  #Generate matlab-like jet colormap with specified amount of colors
  return jetColorMap = Colors.RGB{Float64}[
  Colors.RGB(
  clamp(min(4*ellement - 1.5, -4*ellement + 4.5) ,0.0,1.0),
  clamp(min(4*ellement - 0.5, -4*ellement + 3.5) ,0.0,1.0),
  clamp(min(4*ellement + 0.5, -4*ellement + 2.5) ,0.0,1.0))
  for ellement in range(0.0, stop=1.0, length=numberOfColors)];

  end

  ################################################################
  function hsv(numberOfColors::Int64)
    #Generate matlab-like hsv colormap with specified amount of colors
    return hsvColorMap = range(Colors.HSV(0,1,1), stop=Colors.HSV(330,1,1), length=numberOfColors);

  end
  ################################################################
  function label2rgb(labeledMatrix, inputColorMap = "jet",backgroundColor = [1 1 1], isShuffled = "noshuffle")

    #check if the labeledMatrix contains only integers
    if !isinteger(labeledMatrix)
      error("labeledMatrix is not an Integer matrix");
    end

    #change indexing of the labeled Matrix to start from 1
    labeledMatrix = labeledMatrix +1;
    numberOfEllements =  maximum(labeledMatrix); # number of connected ellements in the Matrix

    if inputColorMap == "jet"
      currentColormap =  jet(numberOfEllements);
    elseif inputColorMap == "hsv"
      currentColormap =  hsv(numberOfEllements);
      #check if custom colormap is loaded
    elseif is(typeof(inputColorMap),Array{Colors.RGB,1})

      if length(inputColorMap) == maximum(labeledMatrix)
        currentColormap = inputColorMap;
      else
        error("Input Colormap size is inconsistent with the Labeled Image");
      end
    else
      error("Invalid Input Colormap");
    end

    #convert input backgroundColor to RGB type
    backgroundColor = Colors.RGB(backgroundColor[1],backgroundColor[2],backgroundColor[3]);

    if  isShuffled == "shuffle"
      #shuffle the colormap
      currentColormap  = shuffle(currentColormap);
    elseif isShuffled == "noshuffle"
      #do nothing
    else
      error("isShuffled can be either shuffle or noshuffle ");
    end

    if (currentColormap[1] != backgroundColor)
      currentColormap = [backgroundColor,currentColormap[1:end-1]];
    end

    indexedImage = Images.ImageCmap(labeledMatrix, currentColormap);
    outputRGB = convert(Images.Image, indexedImage);

    return outputRGB;

  end

  ################################################################
end #End of ImageTools


text:
# ---
# title: 653. Two Sum IV - Input is a BST
# id: problem653
# author: zhwang
# date: 2022-02-23
# difficulty: Easy
# categories: Tree
# link: <https://leetcode.com/problems/two-sum-iv-input-is-a-bst/description/>
# hidden: true
# ---
#
# Given the `root` of a Binary Search Tree and a target number `k`, return
# _`true` if there exist two elements in the BST such that their sum is equal to
# the given target_.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/21/sum_tree_1.jpg)
#
#
#
#     Input: root = [5,3,6,2,4,null,7], k = 9
#     Output: true
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/21/sum_tree_2.jpg)
#
#
#
#     Input: root = [5,3,6,2,4,null,7], k = 28
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: root = [2,1,3], k = 4
#     Output: true
#
#
# **Example 4:**
#
#
#
#     Input: root = [2,1,3], k = 1
#     Output: false
#
#
# **Example 5:**
#
#
#
#     Input: root = [2,1,3], k = 3
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[1, 104]`.
#   * `-104 <= Node.val <= 104`
#   * `root` is guaranteed to be a **valid** binary search tree.
#   * `-105 <= k <= 105`
#
#
## @lc code=start
using LeetCode

function find_target_inorder_traversal(root::TreeNode, k::Int)::Bool
    stack, nums = TreeNode{Int}[], Real[-Inf32] ## lowest bound
    while !isempty(stack) || !isnothing(root)
        while !isnothing(root)
            push!(stack, root)
            root = root.left
        end
        root = pop!(stack)
        while nums[end] + root.val > k
            pop!(nums)
        end
        root.val + nums[end] == k && return true
        push!(nums, root.val)
        root = root.right
    end
    false
end
## @lc code=end

## @lc test=start
@testset "653.two-sum-iv-input-is-a-bst.jl" begin
    @test find_target_inorder_traversal(TreeNode{Int}([5,3,6,2,4,nothing,7]), 9) == true
    @test find_target_inorder_traversal(TreeNode{Int}([5,3,6,2,4,nothing,7]), 28) == false
    @test find_target_inorder_traversal(TreeNode{Int}([2,1,3]), 4) == true
    @test find_target_inorder_traversal(TreeNode{Int}([2,1,3]), 1) == false
    @test find_target_inorder_traversal(TreeNode{Int}([2,1,3]), 3) == true
end
## @lc test=end

## @lc info=start
# link: [solution to 653](https://leetcode-cn.com/problems/two-sum-iv-input-is-a-bst/solution/653-liang-shu-zhi-he-iv-shu-ru-bst-by-re-2h4o/)
## @lc info=end

ext: jl


text:
# HS71
# Polynomial objective and constraints
# min x1 * x4 * (x1 + x2 + x3) + x3
# st  x1 * x2 * x3 * x4 >= 25
#     x1^2 + x2^2 + x3^2 + x4^2 = 40
#     1 <= x1, x2, x3, x4 <= 5
#
# Initial x = (1,5,5,1)
#  (1.000..., 4.743..., 3.821..., 1.379...)

using JuMP, SNOPT7

m = Model(with_optimizer(SNOPT7.Optimizer, print_level=0,system_information="yes"))

x0 = [1,5,5,1]

@variable(m, 1 <= x[i=1:4] <= 5)
@NLobjective(m, Min, x[1]*x[4]*(x[1]+x[2]+x[3]) + x[3])
@NLconstraint(m, x[1]*x[2]*x[3]*x[4] >= 25)
@NLconstraint(m, sum(x[i]^2 for i=1:4) == 40)

JuMP.optimize!(m)

println(JuMP.value.(x))

objval = JuMP.objective_value(m)
println("Final objective: $objval")

println(JuMP.termination_status(m))


text:
using CSV
using DataFrames
using StatsPlots
using Statistics
using Dates

const FULL_FILE_PROCESSED = joinpath("data","Melbourne_housing_FULL_processed_2.csv")

const SOLD_METHODS = ["S", "SP", "PI", "PN", "SN", "SA", "SS"]

"""
    currency_format(num)::String

Takes a number, returns it as a string in currency (\$AUD) form. Taken from [here](https://stackoverflow.com/questions/52213829/in-julia-insert-commas-into-integers-for-printing-like-python-3-6/52213830#52213830).

# Examples
```jldoctest
julia>a = 10000
10000
julia>currency_format(a)
\$10,000
```
"""
currency_format(num)::String = "\$" * replace(num |> string, r"(?<=[0-9])(?=(?:[0-9]{3})+(?![0-9]))" => ",")

function aggregate_date_to_months!(df::AbstractDataFrame, col=:Date)
    transform!(df, col => ByRow(x->DateTime(Dates.year(x), Dates.month(x))) => col)
end
"""
    plot_time_sales(df::AbstractDataFrame)

Returns a line graph of house sales over time, aggregated by month. Only considers houses with a Method in SOLD_METHODS.
"""
function plot_time_sales(df::AbstractDataFrame)
    df_time_sales = combine(groupby(filter!(:Method => x->x in SOLD_METHODS, select(df, :Date, :Method)) |> aggregate_date_to_months!, :Date), nrow => :count)

    xticks_range = DateTime(2015, 12):Dates.Month(4):DateTime(2017, 12)
    xticks=(xticks_range, [Dates.format(tick, "yyyy-mm") for tick in xticks_range])

    return @df df_time_sales plot(:Date, :count, seriestype=:line, xticks=xticks, xlims=(DateTime(2015, 10), DateTime(2018, 3)), linecolor=:red, linewidth=5, label=false, xlabel="Month of Sale", ylabel="Number of Properties Sold", plot_title="Properties Sold By Month")
end

"""
    plot_time_price(df::AbstractDataFrame)

Returns a line graph of mean sale price over time, aggregated by months. Two more lines appear on this chart; they are `mean + std` and `mean - std` respectively.
"""
function plot_time_price(df::AbstractDataFrame)
    df_time_price = combine(groupby(dropmissing!(filter!(:Method => x->x in SOLD_METHODS, select(df, :Date, :Method, :Price)) |> aggregate_date_to_months!, :Price), :Date), nrow => :count, :Price => mean => :Price_mean, :Price => std => :Price_std)

    transform!(df_time_price, [:Price_mean, :Price_std] => ByRow((x, y) -> (x - y, x + y)) => [:Price_mean_nstd, :Price_mean_pstd])

    plot_mean = @df df_time_price plot(:Date, :Price_mean, linewidth=5, seriestype=:line, linecolor=:green, label="Mean")

    yticks_range = 400000:400000:2000000
    yticks = (yticks_range, currency_format.(yticks_range))

    xticks_range = DateTime(2016, 4):Dates.Month(4):DateTime(2017, 12)
    xticks=(xticks_range, [Dates.format(tick, "yyyy-mm") for tick in xticks_range])

    return @df df_time_price plot(plot_mean, :Date, [:Price_mean_nstd, :Price_mean_pstd], linewidth=1, seriestype=:line, linestyle=:dash, linecolor=:blue, ylims=(0, 2000000), yticks=yticks, xticks=xticks, label="", xlabel="Month of Sale", ylabel="Mean Sale Price, w/ St. Dev. (\$AUD)", plot_title="Sale Price Over Time")
end

"""
    plot_time_type(df::AbstractDataFrame)

Returns a line graph of the proportion of houses sold in relation to other properties over time, aggregated by month.
"""
function plot_time_type(df::AbstractDataFrame)
    df_time_type = filter!(:Method => x->x in SOLD_METHODS, select(df, :Date, :Method, :Type) |> aggregate_date_to_months!)

    count_dates(df) = combine(groupby(df, :Date), :Type => x->count((a->a=="h"), x) / length(x))

    df_time_type_all = df_time_type |> count_dates

    xticks_range = DateTime(2016, 2):Dates.Month(4):DateTime(2018, 3)
    xticks=(xticks_range, [Dates.format(tick, "yyyy-mm") for tick in xticks_range])

    yticks_range = 0.5:0.05:0.8
    yticks = (yticks_range, 100 .* yticks_range)

    return @df df_time_type_all plot(:Date, :Type_function, seriestype=:line, linewidth=3, linecolor=:orange, label="", xticks=xticks, xlims = (DateTime(2016, 1), DateTime(2018, 3)), ylims=(0.45, 0.85), yticks=yticks, xlabel="Month of Sale", ylabel="Percentage of Properties Sold as Houses", plot_title = "Hice as a Proportion of Sold Properties", fillrange=0.45, fillcolor=:brown, tickdirection=:out)
end

generate_task3_df(filename=FULL_FILE_PROCESSED) = transform!(DataFrame(CSV.File(FULL_FILE_PROCESSED)), :Date => ByRow(x->Date(x, dateformat"d/m/y")) => :Date) |> sort!


text:
# 2020.04.01, Антонюк В.А.

include("HeuristicAlgorithm.jl")

function UpperSum(A::Array{Int,2})
    S,(H,W) = 0, size(A)
    @assert W == H
    for i=1:H, j=i+1:W  S += A[i,j] end
    return S
end

ReOrder(M::Matrix{Int}, V::Vector{Int}) = [M[i,j] for i ∈ V, j ∈ V]

function RecursiveIteration!(ScoreMatrix::Matrix{Int}, ToDrop=Vector{Int}())
    (H,W) = size(ScoreMatrix)
    @assert W == H
    TrueDropped = copy(ToDrop)
    TrueIndices!(TrueDropped)
    if length(TrueDropped) < H
        Summed = zeros(Int, H-length(TrueDropped))
        Dropped = 0
        for i = 1:H
            i ∈ TrueDropped && continue
            Dropped += 1
            for j = 1:W
                j ∈ TrueDropped && continue
                Summed[Dropped] += ScoreMatrix[i,j]
            end
        end
        (Val,Idx) = findmin(Summed)
        while Idx != nothing
            push!(ToDrop,Idx)
            if length(ToDrop) < H
                RecursiveIteration!(ScoreMatrix, ToDrop)
            else
                Dropped = copy(ToDrop)
                TrueIndices!(Dropped)
                LoopCount = 0
                while true
                    LoopCount += 1
                    OldDropped = copy(Dropped)
                    print("$Dropped : ")
                    print(UpperSum(ReOrder(ScoreMatrix, Dropped)))
                    OptimizeRanking!(Dropped,ScoreMatrix)
                    LoopCount < 100 || break
                    Dropped != OldDropped || break
                    print(" -> ")
                end
                println("\n")
            end
            pop!(ToDrop)
            Idx = findnext(x->x==Val, Summed, Idx+1)
        end
    end
end

max_stars_repo_name: vaa-msu/ranking-tests

max_issues_repo_name: vaa-msu/ranking-tests

max_forks_repo_name: vaa-msu/ranking-tests


text:
"""
Atomic solution for an SDE.

### Parameters

* `DT`: data type
* `TT`: time step type
* `AT`: array type
* `IT`: internal variable types

### Fields

* `t`: time of current time step
* `t̅`: time of previous time step
* `q`: current solution of q
* `q̅`: previous solution of q
* `q̃`: compensated summation error of q
* `p`: current solution of p
* `p̅`: previous solution of p
* `p̃`: compensated summation error of p
* `ΔW`: Wiener process driving the stochastic process q
* `ΔZ`: Wiener process driving the stochastic process q
* `K`:  integer parameter defining the truncation of the increments of the Wiener process (for strong solutions)
"""
mutable struct AtomicSolutionPSDE{DT <: Number, TT <: Real, AT <: AbstractArray{DT}, IT <: NamedTuple} <: AtomicSolution{DT,TT,AT}
    t::TT
    t̅::TT

    q::AT
    q̅::AT
    q̃::AT

    p::AT
    p̅::AT
    p̃::AT

    ΔW::AT
    ΔZ::AT

    K::Int

    internal::IT

    function AtomicSolutionPSDE{DT,TT,AT,IT}(nd, nm, internal::IT) where {DT,TT,AT,IT}
        new(zero(TT), zero(TT),
            AT(zeros(DT, nd)), AT(zeros(DT, nd)), AT(zeros(DT, nd)),
            AT(zeros(DT, nd)), AT(zeros(DT, nd)), AT(zeros(DT, nd)),
            AT(zeros(DT, nm)), AT(zeros(DT, nm)), 0, internal)
    end

    function AtomicSolutionPSDE{DT,TT,AT,IT}(t::TT, q::AT, p::AT, ΔW::AT, ΔZ::AT, internal::IT) where {DT,TT,AT,IT}
        new(zero(t), zero(t),
            zero(q), zero(q), zero(q),
            zero(p), zero(p), zero(p),
            zero(ΔW), zero(ΔZ), 0,
            internal)
    end
end

AtomicSolutionPSDE(::Type{DT}, ::Type{TT}, ::Type{AT}, nd, nm, internal::IT=NamedTuple()) where {DT,TT,AT,IT} = AtomicSolutionPSDE{DT,TT,AT,IT}(nd, nm, internal)
AtomicSolutionPSDE(t::TT, q::AT, p::AT, ΔW::AT, ΔZ::AT, internal::IT=NamedTuple()) where {DT, TT, AT <: AbstractArray{DT}, IT} = AtomicSolutionPSDE{DT,TT,AT,IT}(t, q, p, ΔW, ΔZ, internal)

function Solutions.set_solution!(asol::AtomicSolutionPSDE, sol)
    t, q, p = sol
    asol.t  = t
    asol.q .= q
    asol.p .= p
end

function set_increments!(asol::AtomicSolutionPSDE, incs)
    ΔW, ΔZ = incs
    asol.ΔW .= ΔW
    asol.ΔZ .= ΔZ
end

function Solutions.get_solution(asol::AtomicSolutionPSDE)
    (asol.t, asol.q, asol.p)
end

function get_increments(asol::AtomicSolutionPSDE)
    (asol.ΔW, asol.ΔZ)
end

function get_increments!(asol::AtomicSolutionPSDE, ΔW, ΔZ)
    ΔW .= asol.ΔW
    ΔZ .= asol.ΔZ
end

function Common.reset!(asol::AtomicSolutionPSDE, Δt)
    asol.t̅  = asol.t
    asol.q̅ .= asol.q
    asol.p̅ .= asol.p
    asol.t += Δt
end

function Solutions.update!(asol::AtomicSolutionPSDE{DT}, y::Vector{DT}, z::Vector{DT}) where {DT}
    for k in eachindex(y,z)
        update!(asol, y[k], z[k], k)
    end
end

function Solutions.update!(asol::AtomicSolutionPSDE{DT}, y::DT, z::DT, k::Union{Int,CartesianIndex}) where {DT}
    asol.q[k], asol.q̃[k] = compensated_summation(y, asol.q[k], asol.q̃[k])
    asol.p[k], asol.p̃[k] = compensated_summation(z, asol.p[k], asol.p̃[k])
end


text:
using Images, Colors, GeometryTypes
using Reactive, FileIO, GLVisualize
using GLAbstraction, GeometryTypes, GLWindow, ImageFiltering
window = glscreen()
slider, slider_s = widget(Signal(1f0), range = 0.1f0:0.1f0:1.0f0, window)

f(x,y) = (x,y) == (0,0) ? 0 : -x*y/(x^2 + y^2)

function graph_polar(rad)
    rv = linspace(0,rad,10)
    θv = linspace(0,2π,10)
    xv = Float32[r*cos(θ) for r=rv, θ=θv]
    yv = Float32[r*sin(θ) for r=rv, θ=θv]
    zv = Float32[f(r*cos(θ),r*sin(θ)) for r = rv, θ = θv];
    return (xv,yv,zv)
end

function graph_cartesian(a)
    zv = Float32[f(x,y) for x=linspace(-a,a,20),y=linspace(-a,a,20)];
    return zv
end

startvalue = graph_polar(1.0)
task, surfsig = async_map(graph_polar, startvalue, slider_s)

image_renderable = visualize(
    (map(first, surfsig), map(x-> getindex(x, 2), surfsig), map(last, surfsig)), :surface,
    color_map = GLVisualize.default(Vector{RGBA}, Style(:default)),
    wireframe = true, color = nothing, color_norm = Vec2f0(-0.4, 0.4)
)

w = widths(value(boundingbox(slider)))
h = round(Int, w[2])

slider_screen = Screen(window, area = map(window.area) do a
    SimpleRectangle(0, 0, a.w, h)
end)

image_screen = Screen(window, area = map(window.area) do a
    SimpleRectangle(0, h, a.w, a.h-h)
end)

_view(slider, slider_screen, camera = :fixed_pixel)
_view(image_renderable, image_screen, camera = :perspective)

renderloop(window)


text:
using Plots
include("generate_sim_data.jl")

#num_obs = repeat([30, 100, 1000], 3);  σ_ccf = vcat( fill(1e-7, 3), fill(1e-10, 3), fill(0, 3) )
num_obs = repeat([40, 80, 160, 320, 640, 1280], 1);  σ_ccf = vcat( fill(1e-9, 6) )
 results = Any[]
 label_list = String[]
 for i in 1:length(num_obs)
    (rvs, ccfs) = make_simulated_rvs_and_ccfs_to_test_scalpels(num_rvs=num_obs[i], σ_ccf=σ_ccf[i])
    push!(results, rms_clean_rvs_vs_num_basis_scalpels(rvs,ccfs, max_num_basis=8) )
    str = "N_obs=" * string(num_obs[i]) * " σ=" * string(σ_ccf[i])
    push!(label_list,  str )
 end

num_basis_vectors = collect(0:(length(first(results))-1))
 plt = plot()
 map(i->scatter!(num_basis_vectors,results[i], label=label_list[i], color=i), 1:length(results) )
 map(i->plot!(num_basis_vectors,results[i], label=:none, color=i), 1:length(results) )
 xlabel!("Number of basis vectors used")
 ylabel!("RMS RVs post Scalpels")
 display(plt)
 #yaxis!(:log10)


text:
#Some tools for Particle in Fourier
function eval_fourier_modes{T}(arg::T, N::Int)
    mode=Array(Complex{T},N)
    psin1=(cos(arg)+ im*sin(arg))
    psin=1.0+im*0.0
    for kdx=1:N
       psin=psin*psin1
       @inbounds mode[kdx]=psin
     end
     return mode
end

function eval_fourier_modes{T<:AbstractFloat}(arg::T, N::Int,
                           mode::Array{Complex{T},1})
    psin1=(cos(arg)+ im*sin(arg))
    psin=1.0+im*0.0
    @inbounds for kdx=1:N
       psin=psin*psin1
       @inbounds mode[kdx]=psin
     end
     nothing
end

function eval_fourier_modes0{T<:AbstractFloat}(arg::T, N::Int,
                           mode::Array{Complex{T},1})
    const psin1=(cos(arg)+ im*sin(arg))
    psin=1.0+im*0.0
    mode[1]=psin
    @inbounds for kdx=2:(N+1)
       psin=psin*psin1
       @inbounds mode[kdx]=psin
     end
     nothing
end

function eval_fourier_modes2{T<:AbstractFloat}(arg::T, N::Int)
    mode=Array(Complex{T},2*N+1)
    N0=N+1
    mode[N0]=1.0
    psin1=(cos(arg)+ im*sin(arg))
    psin=1.0+im*0.0
    @inbounds for kdx=1:N
       psin=psin*psin1
       @inbounds mode[N0+kdx]=psin
                 mode[N0-kdx]=conj(psin)
     end
    return mode
end
function eval_fourier_modes2{T<:AbstractFloat}(arg::T, N::Int,
                           mode::Array{Complex{T},1})
    N0::Int=N+1
    mode[N0]=1.0
    const psin1=(cos(arg)+ im*sin(arg))::Complex{T}
    psin=(1.0+im*0.0)::Complex{T}
    @inbounds for kdx=1:N
       psin=psin*psin1
       @inbounds mode[N0+kdx]=psin
                 mode[N0-kdx]=conj(psin)
     end
     nothing
end
# function eval_fourier_modes2{T<:AbstractFloat}(arg::T, N::Int,
#                            mode::Array{Complex{T},1})
#     N0=N+1
#     mode[N0]=1.0
#     mode[N0+1]=(cos(arg)- im*sin(arg))::Complex{T}
#     mode[N0-1]=conj(mode[N0+1])
#     @inbounds for kdx=2:N
#        mode[N0+kdx]=mode[N0+kdx-1]*mode[N0]
#        mode[N0-kdx]=conj(mode[N0-kdx])
#      end
#      nothing
# end

# General structure for a 3d PIF field solver
type pif3d
  Nx1::Int
  Nx2::Int
  Nx3::Int
  L1::Float64
  L2::Float64
  L3::Float64
  N1::Int
  N2::Int
  N3::Int
  function pif3d(Nx1,Nx2,Nx3,L1,L2,L3)
    new(Nx1,Nx2,Nx3,L1,L2,L3,Nx1+1,Nx2*2+1,Nx3*2+1)
  end
  function pif3d(Nx1::Int,Nx2::Int,Nx3::Int,L::Float64)
    pif3d(Nx1,Nx2,Nx3,L,L,L)
  end
end

#
# function rhs_particle0{T <:AbstractFloat}(self::pif3d,
#            x1n::Array{T,1},x2n::Array{T,1}, x3n::Array{T,1},  wn::Array{T,1})
#   Np=length(x1n)
#   @assert Np==length(x2n)
#   @assert Np==length(x3n)
#   @assert Np==length(wn)
#
#
#   rhs=zeros(Complex{T},self.N1,self.N2,self.N3)
#
#   k1=2*pi/self.L1
#   k2=2*pi/self.L2
#   k3=2*pi/self.L3
#
#   @inbounds begin
#   for pdx=1:Np
#     #Julia is column major
#
#     for kdx=1:self.N3
#        for jdx=1:self.N2
#         @simd  for idx=1:self.N1
#           arg::Complex{T}=k1*x1n[pdx]+k2*x2n[pdx]+k3*x3n[pdx]
#            rhs[idx,jdx,kdx]+=(cos(arg)+ im*sin(arg))*wn[pdx]
#         end
#       end
#     end
#
#   end
#   end
#   return conj(rhs)/self.L1/self.L2/self.L3
# end

function rhs_particle{T <:AbstractFloat}(self::pif3d,
           x1n::Array{T,1},x2n::Array{T,1}, x3n::Array{T,1},  wn::Array{T,1})
  Np=length(x1n)
  @assert Np==length(x2n)
  @assert Np==length(x3n)
  @assert Np==length(wn)

  rhs=zeros(Complex{T},self.N1,self.N2,self.N3)

  k1=2*pi/self.L1
  k2=2*pi/self.L2
  k3=2*pi/self.L3

  @inbounds begin
  for pdx=1:Np
    #Julia is column major
    eval_fourier_modes0(x1n[pdx]*k1, self.Nx1,psi1)
    eval_fourier_modes2(x2n[pdx]*k2, self.Nx2,psi2)
    eval_fourier_modes2(x3n[pdx]*k3, self.Nx3,psi3)

    #Introduce weight
    @simd for idx=1:self.N1
    psi1[idx]*=wn[pdx]
    end

    for kdx=1:self.N3
       for jdx=1:self.N2
        psi23::Complex{T}=psi2[jdx]*psi3[kdx]
        @simd  for idx=1:self.N1
           rhs[idx,jdx,kdx]+=psi1[idx]*psi23
        end
      end
    end
    #Accumulate
    # for idx=1:self.N1
    #   for jdx=1:self.N2
    #     @simd for kdx=1:self.N3
    #       rhs[idx,jdx,kdx]+=psi1[idx]*psi2[jdx]*psi3[kdx]
    #     end
    #   end
    # end
  end
  end
  return conj(rhs)/self.L1/self.L2/self.L3
end

function rhs_particle2{T <:AbstractFloat}(self::pif3d,
           x1n::Array{T,1},x2n::Array{T,1}, x3n::Array{T,1},  wn::Array{T,1})
  Np=length(x1n)
  @assert Np==length(x2n)
  @assert Np==length(x3n)
  @assert Np==length(wn)

  rhs=zeros(Complex{T},self.N1,self.N2,self.N3)
  psi1=Array(Complex{T},self.N1)
  psi2=Array(Complex{T},self.N2)
  psi3=Array(Complex{T},self.N3)

  k1=2*pi/self.L1
  k2=2*pi/self.L2
  k3=2*pi/self.L3

  @inbounds begin
  for pdx=1:Np
    #Julia is column major
    eval_fourier_modes0(x1n[pdx]*k1, self.Nx1,psi1)
    eval_fourier_modes2(x2n[pdx]*k2, self.Nx2,psi2)
    eval_fourier_modes2(x3n[pdx]*k3, self.Nx3,psi3)

    #Introduce weight
    psi1*=wn[pdx]
    for kdx=1:self.N3
       BLAS.ger!(psi3[kdx],psi1,psi2,rhs[:,:,kdx])
    end
    #for idx=1:self.N1
    #  BLAS.ger!(psi1[idx],psi2,psi3,rhs[idx,:,:])
    #end

  end
  end
  return conj(rhs)/self.L1/self.L2/self.L3
end

function eval_scalar{T}(self::pif3d,
           x1n::Array{T,1},x2n::Array{T,1}, x3n::Array{T,1},
           F::Array{Complex{T},3})
  Np=length(x1n)
  @assert Np==length(x2n)
  @assert Np==length(x3n)

  Fn=zeros(Complex{T},Np)
  psi1=Array(Complex{T},self.N1)
  psi2=Array(Complex{T},self.N2)
  psi3=Array(Complex{T},self.N3)
  k1=2*pi/self.L1
  k2=2*pi/self.L2
  k3=2*pi/self.L3

  for pdx=1:Np
    eval_fourier_modes0(x1n[pdx]*k1, self.Nx1,psi1)
    eval_fourier_modes2(x2n[pdx]*k2, self.Nx2,psi2)
    eval_fourier_modes2(x3n[pdx]*k3, self.Nx3,psi3)
    for kdx=1:self.N3
       for jdx=1:self.N2
        psi23::Complex{T}=psi2[jdx]*psi3[kdx]
        Fn[pdx]+=psi23*F[1,jdx,kdx]
        @simd  for idx=2:self.N1
          psin::Complex{T}=psi1[idx]*psi23
           Fn[pdx]+=(real(psin)*real(F[idx,jdx,kdx])
                    -imag(psin)*imag(F[idx,jdx,kdx]))/2.0
        end
      end
    end

  end
  return real(Fn)
end

function eval_scalar{T}(self::pif3d,
           x1n::Array{T,1},x2n::Array{T,1}, x3n::Array{T,1},
           F::Array{Array{Complex{T},3}})
  Np=length(x1n)
  @assert Np==length(x2n)
  @assert Np==length(x3n)

  D=length(F)
  Fn=zeros(Complex{T},Np,D)
  psi1=Array(Complex{T},self.N1)
  psi2=Array(Complex{T},self.N2)
  psi3=Array(Complex{T},self.N3)
  k1=2*pi/self.L1
  k2=2*pi/self.L2
  k3=2*pi/self.L3

  for pdx=1:Np
    eval_fourier_modes0(x1n[pdx]*k1, self.Nx1,psi1)
    eval_fourier_modes2(x2n[pdx]*k2, self.Nx2,psi2)
    eval_fourier_modes2(x3n[pdx]*k3, self.Nx3,psi3)
    for kdx=1:self.N3
       for jdx=1:self.N2
        psi23::Complex{T}=psi2[jdx]*psi3[kdx]
        for ddx=1:D
            Fn[pdx,ddx]+=psi23*F[ddx][1,jdx,kdx]
            end
            for ddx=1:D
          @simd  for idx=2:self.N1
            psin::Complex{T}=psi1[idx]*psi23
            Fn[pdx,ddx]+=(real(psin)*real(F[ddx][idx,jdx,kdx])
                          -imag(psin)*imag(F[ddx][idx,jdx,kdx]))/2.0
          end
        end
      end
    end

  end
  return real(Fn)
end

function getK(self::pif3d)
  K1=collect(0:self.Nx1)*self.L1*2*pi
  K2=collect(-self.Nx2:self.Nx2)*self.L2*2*pi
  K3=collect(-self.Nx3:self.Nx3)*self.L3*2*pi
  return K1,K2,K3
end

function getK2(self::pif3d)
  K1,K2,K3=getK(self)
  K1=reshape(K1,self.N1,1,1)
  K2=reshape(K2,1,self.N2,1)
  K3=reshape(K3,1,1,self.N3)
  return K1,K2,K3
end

function solve_poisson!{T}(self::pif3d,
                      rho::Array{Complex{T},3}, Phi::Array{Complex{T},3})
    K1,K2,K3=getK(self)
     for kdx=1:self.N3
        for jdx=1:self.N2
          for idx=1:self.N1
           if (idx!=1 && jdx!=self.Nx2+1  && kdx!=self.Nx3+1 )

            Phi[idx,jdx,kdx]=1/(K1[idx]^2 +K2[jdx]^2 +K3[kdx]^2)
          end
         end
       end
     end
end
function solve_poisson{T}(self::pif3d, rho::Array{Complex{T},3})
  Phi=zeros(rho)
  solve_poisson!(self,rho,Phi)
  return Phi
end

function solve_gauss!{T}(self::pif3d, rho::Array{Complex{T},3}, E1,E2,E3   )
  Phi=solve_poisson(self,rho)
  K1,K2,K3=getK2(self)
  E1=-Phi*im.*K1
  E2=-Phi*im.*K2
  E3=-Phi*im.*K3
end

# Gradient of a scalar function
function gradient{T}(self::pif3d, F::Array{Complex{T},3})
  K1,K2,K3=getK2(self)
  DF1=F*im.*K2
  DF2=F*im.*K2
  DF3=F*im.*K3
  return DF1,DF2,DF3
end

function div{T}(self::pif3d, F::Array{Complex{T},3})
  DF1,DF2,DF3=gradient(pif3d, F)
  return DF1.+DF2.+DF3
end

function curl{T}(self::pif3d, F1::Array{Complex{T},3},
              F2::Array{Complex{T},3},F3::Array{Complex{T},3} )
    K1,K2,K3=getK2(self)
    return im*(K2.*F3 - K3.*F2), im*(K3.*F1 - K1.*F3), im*(K1.*F2 - K2.*F1)
end

Np=Int(1e5)
x1n=rand(Float64,Np)
x2n=rand(Float64,Np)
x3n=rand(Float64,Np)

# psix1=Array(Complex128,Nx1*2+1)
# @time eval_fourier_modes2(x1n[1]*2*pi,Nx1,psix1)
# @time psix1=eval_fourier_modes2(x1n[1]*2*pi,Nx1)

# psix1*reshape(psix1,1,Nf1)
wn=x1n

self=pif3d(10,20,1, 1.0)

# wn=f0(x1n,x2n,x3n,vn)./g0(x1n,x2n,x3n,vn)
@time rho=rhs_particle(self,x1n,x2n,x3n,wn)/Np
@time rho=rhs_particle2(self,x1n,x2n,x3n,wn)/Np

sum(abs(rho-rho2)[:])
Phi=zeros(rho)
Phi=solve_poisson(self, rho)

E1,E2,E3=gradient(self,-Phi)

E=[E1,E2,E3]

E[1,2][1]
@time (E1n=eval_scalar(self,x1n,x2n,x3n,E2);E2n=eval_scalar(self,x1n,x2n,x3n,E2);E3n=eval_scalar(self,x1n,x2n,x3n,E3))

@time En=eval_scalar(self,x1n,x2n,x3n,[E1,E2,E3])

k0=0.5
dt=0.1
tmax=15

L=2*pi/k0

function f0(xi1,xi2,xi3,v)
  return exp(-0.5.*v.^2)/sqrt(2*pi).*(1+0.1.*cos(xi1))
end

function g0(xi1,xi2,xi3,v)
  return exp(-0.5.*v.^2)/sqrt(2*pi)
end

max_stars_repo_name: ameresj/VlasovPIF

max_issues_repo_name: ameresj/VlasovPIF

max_forks_repo_name: ameresj/VlasovPIF


text:
using Documenter, PolynomialTestSystems

makedocs(
    sitename = "PolynomialTestSystems.jl",
    pages = [
        "PolynomialTestSystems" => "index.md",
    ]
)

deploydocs(
    repo   = "github.com/JuliaHomotopyContinuation/PolynomialTestSystems.jl.git",
)


text:
immutable Protocol{N,C,U}
    insertread::Int
    insertpos::UnitRange{Int}
    cellpos::NTuple{N,UnitRange{Int}}
    umipos::NTuple{N,UnitRange{Int}}
    readnames::NTuple{N,String}
    readlengths::NTuple{N,Int}
end

# const config_path = joinpath(Pkg.dir("FASTQDemultiplexer"),"config","protocols")
const config_path = joinpath(@__DIR__,"..","config","protocols")
const protocol_config_paths = Dict(:marsseq => joinpath(config_path,"marsseq.yml"),
                                   :dropseq => joinpath(config_path,"dropseq.yml"),
                                   :x10 => joinpath(config_path,"10x.yml"))

function Protocol(prot::Symbol)

    allowed_protocols = keys(protocol_config_paths)

    if prot in allowed_protocols
        return Protocol(protocol_config_paths[prot])
    else
        error("Unknown protocol: $prot, use one of $allowed_protocol")
    end
end

function Protocol(yamlfile::String)
    config = YAML.load_file(yamlfile)

    nreads = length(config["reads"])

    cellpos = Array{UnitRange{Int}}(nreads)
    umipos  = Array{UnitRange{Int}}(nreads)
    readnames = Array{String}(nreads)
    readlengths = Array{Int}(nreads)
    insertpos = 1:0
    insertread = 0
    insertdefined = false

    # parse the reads
    for (i,read) in enumerate(config["reads"])

        if haskey(read,"cell")
            # TODO add safer parsing?
            cellpos[i] = eval(parse(read["cell"]))
        else
            cellpos[i] = 1:0
        end

        if haskey(read,"umi")
            umipos[i] = eval(parse(read["umi"]))
        else
            umipos[i] = 1:0
        end

        if haskey(read,"insert")
            if insertdefined
                error("The field insert can be defined in only one read")
            end
            insertpos = eval(parse(read["insert"]))
            insertread = i
            insertdefined = true
        end

        if haskey(read,"name")
            readnames[i]=read["name"]
        else
            error("Name of the read is required")
        end

        if haskey(read,"length")
            readlengths[i]=read["length"]
        else
            error("Length of the read is required")
        end

    end

    C = sum(map(length,cellpos))
    U = sum(map(length,umipos))

    return Protocol{nreads,C,U}(insertread,
                                insertpos,
                                (cellpos...),
                                (umipos...),
                                (readnames...),
                                (readlengths...))

end

cellidtoname(id::UInt,::Protocol{N,C,U}) where {N,C,U} = idtoname(id,C)

umiidtoname(id::UInt,::Protocol{N,C,U}) where {N,C,U} = idtoname(id,U)

max_stars_repo_name: pwl/FASTQDemultiplexer

max_issues_repo_name: pwl/FASTQDemultiplexer

max_forks_repo_name: pwl/FASTQDemultiplexer


text:
export StructLogicCircuit, PlainStructLogicCircuit,
    PlainStructLogicLeafNode, PlainStructLogicInnerNode,
    PlainStructLiteralNode, PlainStructConstantNode, PlainStructTrueNode, PlainStructFalseNode,
    PlainStruct⋀Node, PlainStruct⋁Node,
    vtree, vtree_safe, prime, sub

#####################
# Logic circuits that are structured,
# meaning that each conjunction is associated with a vtree node.
#####################

"Root of the structure logic circuit node hierarchy"
abstract type StructLogicCircuit <: LogicCircuit end

"Root of the plain structure logic circuit node hierarchy"
abstract type PlainStructLogicCircuit <: StructLogicCircuit end

"A plain structured logical leaf node"
abstract type PlainStructLogicLeafNode <: PlainStructLogicCircuit end

"A plain structured logical inner node"
abstract type PlainStructLogicInnerNode <: PlainStructLogicCircuit end

"A plain structured logical literal leaf node, representing the positive or negative literal of its variable"
mutable struct PlainStructLiteralNode <: PlainStructLogicLeafNode
    literal::Lit
    vtree::Vtree
    data
    counter::UInt32
    PlainStructLiteralNode(l,v) = begin
        @assert lit2var(l) ∈ v
        new(l, v, nothing, 0)
    end
end

"""
A plain structured logical constant leaf node, representing true or false.
These are the only structured nodes that don't have an associated vtree node (cf. SDD file format)
"""
abstract type PlainStructConstantNode <: PlainStructLogicInnerNode end

"A plain structured logical true constant. Never construct one, use `structtrue` to access its unique instance"
mutable struct PlainStructTrueNode <: PlainStructConstantNode
    data
    counter::UInt32
end

"A plain structured logical false constant.  Never construct one, use `structfalse` to access its unique instance"
mutable struct PlainStructFalseNode <: PlainStructConstantNode
    data
    counter::UInt32
end

"A plain structured logical conjunction node"
mutable struct PlainStruct⋀Node <: PlainStructLogicInnerNode
    prime::PlainStructLogicCircuit
    sub::PlainStructLogicCircuit
    vtree::Vtree
    data
    counter::UInt32
    PlainStruct⋀Node(p,s,v) = begin
        @assert isinner(v) "Structured conjunctions must respect inner vtree node"
        @assert isconstantgate(p) || varsubset_left(vtree(p),v) "$p does not go left in $v"
        @assert isconstantgate(s) || varsubset_right(vtree(s),v) "$s does not go right in $v"
        new(p,s, v, nothing, 0)
    end
end

"A plain structured logical disjunction node"
mutable struct PlainStruct⋁Node <: PlainStructLogicInnerNode
    children::Vector{PlainStructLogicCircuit}
    vtree::Vtree # could be leaf or inner
    data
    counter::UInt32
    PlainStruct⋁Node(c,v) = new(c, v, nothing, 0)
end

"The unique plain structured logical true constant"
const structtrue = PlainStructTrueNode(nothing, 0)

"The unique splain tructured logical false constant"
const structfalse = PlainStructFalseNode(nothing, 0)

#####################
# traits
#####################

@inline GateType(::Type{<:PlainStructLiteralNode}) = LiteralGate()
@inline GateType(::Type{<:PlainStructConstantNode}) = ConstantGate()
@inline GateType(::Type{<:PlainStruct⋀Node}) = ⋀Gate()
@inline GateType(::Type{<:PlainStruct⋁Node}) = ⋁Gate()

#####################
# methods
#####################

@inline constant(n::PlainStructTrueNode)::Bool = true
@inline constant(n::PlainStructFalseNode)::Bool = false
@inline children(n::PlainStruct⋁Node) = n.children
@inline children(n::PlainStruct⋀Node) = [n.prime,n.sub]

"Get the vtree corresponding to the argument"
@inline vtree(n::StructLogicCircuit) = n.vtree

"Get the vtree corresponding to the argument, or nothing if the node has no vtree"
@inline vtree_safe(n::PlainStructLogicInnerNode) = vtree(n)
@inline vtree_safe(n::PlainStructLiteralNode) = vtree(n)
@inline vtree_safe(::PlainStructConstantNode) = nothing

# alias some SDD terminology: primes and subs
"Get the prime, that is, the first conjunct"
@inline prime(n) = n.prime

"Get the sub, that is, the second and last conjunct"
@inline sub(n) = n.sub

# StructLogicCircuit has a default argument for respects: its root's vtree
respects_vtree(circuit::StructLogicCircuit) =
    respects_vtree(circuit, vtree(circuit))

#####################
# constructors and compilation
#####################

conjoin(arguments::Vector{<:PlainStructLogicCircuit};
        reuse=nothing, use_vtree=nothing) =
        conjoin(arguments...; reuse, use_vtree)

function conjoin(a1::PlainStructLogicCircuit,
                 a2::PlainStructLogicCircuit;
                 reuse=nothing, use_vtree=nothing)
    reuse isa PlainStruct⋀Node && reuse.prime == a1 && reuse.sub == a2 && return reuse
    !(use_vtree isa Vtree) && (reuse isa PlainStructLogicCircuit) &&  (use_vtree = reuse.vtree)
    if isconstantgate(a1) && isconstantgate(a2) && !(use_vtree isa Vtree)
        # constant nodes don't have a vtree: resolve to a constant
        return PlainStructLogicCircuit(istrue(a1) && istrue(a2))
    end
    !(use_vtree isa Vtree) && (use_vtree = find_inode(vtree_safe(a1), vtree_safe(a2)))
    return PlainStruct⋀Node(a1, a2, use_vtree)
end

function disjoin(arguments::Vector{<:PlainStructLogicCircuit};
                 reuse=nothing, use_vtree=nothing)
    @assert length(arguments) > 0
    reuse isa PlainStruct⋁Node && reuse.children == arguments && return reuse
    !(use_vtree isa Vtree) && (reuse isa PlainStructLogicCircuit) &&  (use_vtree = reuse.vtree)
    if all(isconstantgate, arguments) && !(use_vtree isa Vtree)
        # constant nodes don't have a vtree: resolve to a constant
        return PlainStructLogicCircuit(any(constant, arguments))
    end
    !(use_vtree isa Vtree) && (use_vtree = mapreduce(vtree_safe, lca, arguments))
    return PlainStruct⋁Node(arguments, use_vtree)
end

# Syntactic sugar for compile with a vtree
(t::Tuple{<:Type,<:Vtree})(arg) = compile(t[1], t[2], arg)
(t::Tuple{<:Vtree,<:Type})(arg) = compile(t[2], t[1], arg)

compile(n::StructLogicCircuit, args...) =
    compile(typeof(n), root(vtree(n)), args...)

# claim `PlainStructLogicCircuit` as the default `LogicCircuit` implementation that has a vtree
compile(::Type{StructLogicCircuit}, args...) =
    compile(PlainStructLogicCircuit, args...)

compile(vtree::Vtree, arg) =
    compile(StructLogicCircuit,vtree,arg)

compile(::Type{<:PlainStructLogicCircuit}, ::Vtree, b::Bool) =
    compile(PlainStructLogicCircuit, b)

compile(::Type{<:PlainStructLogicCircuit}, b::Bool) =
    b ? structtrue : structfalse

compile(::Type{<:PlainStructLogicCircuit}, vtree::Vtree, l::Lit) =
    PlainStructLiteralNode(l,find_leaf(lit2var(l),vtree))

function compile(::Type{<:PlainStructLogicCircuit}, vtree::Vtree, circuit::LogicCircuit)
    f_con(n) = compile(PlainStructLogicCircuit, constant(n))
    f_lit(n) = compile(PlainStructLogicCircuit, vtree, literal(n))
    f_a(n, cns) = conjoin(cns...) # note: this will use the LCA as vtree node
    f_o(n, cns) = disjoin(cns) # note: this will use the LCA as vtree node
    foldup_aggregate(circuit, f_con, f_lit, f_a, f_o, PlainStructLogicCircuit)
end

fully_factorized_circuit(::Type{StructLogicCircuit}, vtree::Vtree) =
    fully_factorized_circuit(PlainStructLogicCircuit, vtree)

function fully_factorized_circuit(::Type{<:PlainStructLogicCircuit}, vtree::Vtree)
    f_leaf(l) = begin
        v = variable(l)
        pos = compile(PlainStructLogicCircuit, vtree, var2lit(v))
        neg = compile(PlainStructLogicCircuit, vtree, -var2lit(v))
        pos | neg
    end
    f_inner(i,cs) = disjoin([conjoin(cs)])
    c = foldup_aggregate(vtree, f_leaf, f_inner, PlainStructLogicCircuit)
    disjoin([c]) # "bias term"
end

isstruct_decomposable(::PlainStructLogicCircuit) = true


text:
# highest_plate_idx should be greater-than-or-equal-to the largest
# value found in any field of any element of `cuts`.
@timeit TIMER function _reachable_plate_types(
	cuts :: Vector{NTuple{3, P}}, highest_plate_idx :: P
) where {P}
	children = [P[] for _ in 1:highest_plate_idx]
	for (pp, fc, sc) in cuts
		!iszero(fc) && push!(children[pp], fc)
		!iszero(sc) && push!(children[pp], sc)
	end
	reached = falses(highest_plate_idx)
	num_reached = zero(P)
	visit_list = P[one(P)]
	next_of_list = one(P)
	while next_of_list <= length(visit_list)
		pp = visit_list[next_of_list]
		if !reached[pp]
			num_reached += one(P)
			reached[pp] = true
			append!(visit_list, [c for c in children[pp] if !reached[c]])
		end
		next_of_list += one(P)
	end
	return num_reached, reached
end

# Works for both cuts and extractions, both have the parent plate as their
# first field.
@timeit TIMER function _reachable_carves(carves, reachable_plates)
	# We start from the premise the extractions are reachable because this is
	# more probable than the opposite.
	qt_rc = length(carves) # TODO: get the correct number type for this?
	rc = trues(qt_rc)
	for (carve_idx, carve) in pairs(carves)
		pp = first(carve) # gets the parent plate
		if !reachable_plates[pp]
			qt_rc -= one(qt_rc)
			rc[carve_idx] = false
		end
	end
	return qt_rc, rc
end

@timeit TIMER function _reachable(
	np :: Vector{Tuple{P, D}},
	cuts :: Vector{NTuple{3, P}},
	highest_plate_idx :: P
) where {D, P}
	qt_rp, rp = _reachable_plate_types(cuts, highest_plate_idx)
	qt_rc, rc = _reachable_carves(cuts, rp)
	qt_re, re = _reachable_carves(np, rp)
	return qt_re, re, qt_rc, rc, qt_rp, rp
end

@timeit TIMER function _print_unreachable_plates(
	num_reached, reachable, pli_lwb, debug
)
	total_num_plates = length(pli_lwb)
	@assert num_reached <= total_num_plates
	@assert length(reachable) == total_num_plates
	!debug && return
	qt_unreachable_plate_types = total_num_plates - num_reached
	@show qt_unreachable_plate_types
	# This loop increase the log files by orders of magnitude. Useful for
	# extensive debug, but not for experiment logs (in which hundreds
	# of files are created).
	#=
	iob = IOBuffer()
	write(iob, "START_UNREACHABLE_PLATES\n")
	if qt_unreachable_plate_types > 0
		for (plate_idx, isreachable) in zip(1:total_num_plates, reachable)
			if !isreachable
				l, w, b = pli_lwb[plate_idx]
				write(iob, "i $plate_idx l $l w $w b $b\n")
			end
		end
	end
	write(iob, "END_UNREACHABLE_PLATES\n")
	print(read(seekstart(iob), String))
	=#
	return
end

@timeit TIMER function _purge_unreachable!(
	bp :: ByproductPPG2KP{D, S, P}, model, qt_re, re, qt_rc, rc, qt_rp, rp
) where {D, S, P}
	@assert length(bp.np) == length(re)
	@assert length(bp.np) == length(model[:picuts])
	@assert length(bp.cuts) == length(rc)
	@assert length(bp.cuts) == length(model[:cuts_made])
	@assert length(bp.pli_lwb) == length(rp)
	@assert length(bp.pli_lwb) == length(model[:plate_cons])

	# The variables can be deleted first without problem, they make reference
	# to the plates, but the plates have no reference back to the variables.
	qt_re < length(bp.np) && (bp = _delete_vars!(bp, model, re, :picuts))
	qt_rc < length(bp.cuts) && (bp = _delete_vars!(bp, model, rc, :cuts_made))
	# All the code below is just to take care of the case in which some
	# plate types (i.e., constraints) are not reachable (this because all
	# bp.np and bp.cuts elements may need update in this case).
	qt_rp == length(bp.pli_lwb) && return bp

	# NOTE: the code below leaves the positions corresponding to the old index of
	# deleted plates undefined. The assumption is that they will not be accessed,
	# i.e., that no cut or extraction that is kept by this method will ever refer
	# to a removed plate (i.e., that the unreachable plates are, in fact, not
	# reachable from reachable cuts or extractions).
	# `new_plate_idxs[y] == x` where `y` is the old index and `x` the new one.
	new_plate_idxs = Vector{P}(undef, length(bp.pli_lwb))
	next_new_idx = one(P)
	for (old_idx, isreachable) in zip(keys(new_plate_idxs), rp)
		if isreachable
			new_plate_idxs[old_idx] = next_new_idx
			next_new_idx += one(P)
		end
	end

	# Now, for the kept cuts and extractions, we update them to refer to the new
	# id/position of the plates (not the old ones).
	for (cut_idx, (pp, fc, sc)) in pairs(bp.cuts)
		@assert rp[pp]
		@assert rp[fc]
		@assert iszero(sc) || rp[sc]
		new_pp, new_fc = new_plate_idxs[pp], new_plate_idxs[fc]
		new_sc = iszero(sc) ? sc : new_plate_idxs[sc]
		if pp != new_pp || fc != new_fc || sc != new_sc
			bp.cuts[cut_idx] = (new_pp, new_fc, new_sc)
		end
	end
	for (extraction_idx, (pp, ep)) in pairs(bp.np)
		@assert rp[pp]
		new_pp = new_plate_idxs[pp]
		pp != new_pp && (bp.np[extraction_idx] = (new_pp, ep))
	end

	# Then we remove the constraints and the new indexes of cuts and
	# extractions are now correct.
	kept_cons, blot_cons = _partition_by_bits(rp, model[:plate_cons])
	#@show @which JuMP.delete(model, blot_cons)
	#@show @which MOI.delete(backend(model), index.(blot_cons))
	#MOI.delete(backend(model), index.(blot_cons))
	JuMP.delete(model, blot_cons)
	model[:plate_cons] = kept_cons
	deleteat!(bp.pli_lwb, .!rp)

	@assert bp.first_vertical_cut_idx <= length(bp.cuts) + 1
	@assert issorted(
		@view bp.cuts[1:(bp.first_vertical_cut_idx-1)]; by = c -> c[PARENT]
	)
	# If there are no vertical cuts this creates an empty array, which is
	# sorted by definition.
	@assert	issorted(
		@view bp.cuts[bp.first_vertical_cut_idx:end]; by = c -> c[PARENT]
	)
	return bp
end

@timeit TIMER function _handle_unreachable!(bp, model, debug, purge)
	qt_re, re, qt_rc, rc, qt_rp, rp = _reachable(
		bp.np, bp.cuts, lastindex(bp.pli_lwb)
	)
	# The unreachable variables are not printed because they can be derived
	# from the set of unreachable plates. Every unreachable variable is just a
	# variable/cut in which the parent plate is an unreachable plate.
	_print_unreachable_plates(qt_rp, rp, bp.pli_lwb, debug)
	@assert qt_re <= length(bp.np) && qt_re >= zero(qt_re)
	@assert qt_rc <= length(bp.cuts) && qt_rc >= zero(qt_rc)
	@assert qt_rp <= length(bp.pli_lwb) && qt_rp >= zero(qt_rp)
	if debug
		println("qt_reachable_extractions = $qt_re")
		println("qt_reachable_cuts = $qt_rc")
		println("qt_reachable_plate_types = $qt_rp")
	end
	if purge && (qt_re + qt_rc + qt_rp) > 0
		bp = _purge_unreachable!(bp, model, qt_re, re, qt_rc, rc, qt_rp, rp)
	end
	return bp
end


text:
mutable struct World
    x_min
    x_max
    v_min
    v_max
    Pset::Vector{Polygonic}
    function World(x_min, x_max, v_min, v_max, Pset)
        new(x_min, x_max, v_min, v_max, Pset)
    end
end

@inline function isValid(this::World, s_q::Vec4f)
    # check if the sampled point is inside the world"
    # @inbounds !(this.x_min[1]<s_q[1]<this.x_max[1]) && return false
    # @inbounds !(this.x_min[2]<s_q[2]<this.x_max[2]) && return false
    if !(this.x_min[1]<s_q[1]<this.x_max[1])
        # println("x wrong ", this.x_min[1], " ", s_q[1], " ", this.x_max[1])
        return false
    end
    if !(this.x_min[2]<s_q[2]<this.x_max[2])
        # println("y wrong ", this.x_min[2], " ", s_q[2], " ", this.x_max[2])
        return false
    end
    # @inbounds !(this.v_min[1]<s_q[3]<this.v_max[1]) && return false
    # @inbounds !(this.v_min[2]<s_q[4]<this.v_max[2]) && return false
    if !(this.v_min[1]<s_q[3]<this.v_max[1])
        # println("v1 wrong ", s_q[3])
        return false
    end
    if !(this.v_min[2]<s_q[4]<this.v_max[2])
        # println("v2 wrong ", s_q[4])
        return false
    end
    for P in this.Pset
        if isInside(P, Vec2f(s_q[1:2]))
            # println("pos wrong ", s_q[1:2])
            return false
        end
    end
    return true
end

@inline function isValid(this::World, q_set::Vector{Vec4f})
    # check validity for multiple points.
    # will be used for piecewize path consited of multiple points
    for q in q_set
        !isValid(this, q) && return false
    end
    return true
end

@inline function isIntersect(this::World, q1::Vec4f, q2::Vec4f)
    for P in this.Pset
        isIntersect(P, q1, q2) && return true
    end
    return false
end

function show(this::World)
    p1 = [this.x_min[1], this.x_min[2]]
    p2 = [this.x_min[1], this.x_max[2]]
    p3 = [this.x_max[1], this.x_max[2]]
    p4 = [this.x_max[1], this.x_min[2]]
    # axis("scaled")
    plot([p1[1], p2[1]], [p1[2], p2[2]], "k-")
    plot([p2[1], p3[1]], [p2[2], p3[2]], "k-")
    plot([p3[1], p4[1]], [p3[2], p4[2]], "k-")
    plot([p4[1], p1[1]], [p4[2], p1[2]], "k-")
    for P in this.Pset
        show(P)
    end
end


text:
"""
    MinTime <: AbstractSequencingAlgorithm

Algorithm that decides the sequence based on minimizing the time from the current
stop to the next stop.
"""
struct MinTime <: AbstractSequencingAlgorithm end

function solve(alg::MinTime, route::AbstractRoute, times::TravelTimes; initial_stop=station(route))
    T = times.travel_times
    return solve(alg, route, T, initial_stop=initial_stop)
end

function solve(alg::MinTime, route::AbstractRoute, T::AbstractMatrix;
               initial_stop=station(route))

    S = stops(route)
    num_stops = length(S)

    # salida
    route_id = name(route)
    seq = Vector{BareStop}()
    out = Sequence(route_id, seq)

    # encontrar la base y guardarla en el vector de salida
    x = initial_stop |> name
    push!(out, BareStop(x))

    # para cada parada, encontrar la siguiente con el menor tiempo
    while length(out) < num_stops
        s = sort(view(T, x, :))

        i = 2 # ignoro la posicion actual (tiempo = 0)
        y = s.dicts[1].keys[i]
        while (y ∈ out) # iterar hasta encontrar uno que no haya visitado
            i += 1
            i > length(s.dicts[1]) && return out
            y = s.dicts[1].keys[i]
        end
        push!(out, BareStop(y))
        x = y
    end

    return out
end


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

Must provide either `loads_kw` or [`doe_reference_name` and `city`] or `doe_reference_name`.

When only `doe_reference_name` is provided the `Site.latitude` and `Site.longitude` are used to look up the ASHRAE climate zone, which determines the appropriate DoE Commercial Reference Building profile.

When using the [`doe_reference_name` and `city`] option, choose `city` from one of the
cities used to represent the ASHRAE climate zones:
- Albuquerque
- Atlanta
- Baltimore
- Boulder
- Chicago
- Duluth
- Fairbanks
- Helena
- Houston
- LosAngeles
- Miami
- Minneapolis
- Phoenix
- SanFrancisco
- Seattle
and `doe_reference_name` from:
- FastFoodRest
- FullServiceRest
- Hospital
- LargeHotel
- LargeOffice
- MediumOffice
- MidriseApartment
- Outpatient
- PrimarySchool
- RetailStore
- SecondarySchool
- SmallHotel
- SmallOffice
- StripMall
- Supermarket
- Warehouse
- FlatLoad

Each `city` and `doe_reference_name` combination has a default `annual_kwh`, or you can provide your
own `annual_kwh` or `monthly_totals_kwh` and the reference profile will be scaled appropriately.
"""
mutable struct ElectricLoad  # mutable to adjust (critical_)loads_kw based off of (critical_)loads_kw_is_net
    loads_kw::Array{Real,1}
    year::Int
    critical_loads_kw::Array{Real,1}
    loads_kw_is_net::Bool
    critical_loads_kw_is_net::Bool

    function ElectricLoad(;
        loads_kw::Union{Missing, Array{<:Real,1}} = missing,
        year::Int = 2020,
        doe_reference_name::Union{Missing, String} = missing,
        city::String = "",
        annual_kwh::Union{Real, Nothing} = nothing,
        monthly_totals_kwh::Array{<:Real,1} = Real[],
        critical_loads_kw::Union{Missing, Array{Real,1}} = missing,
        loads_kw_is_net::Bool = true,
        critical_loads_kw_is_net::Bool = false,
        critical_load_pct::Real = 0.5,
        latitude::Float64,
        longitude::Float64
        )

        if !ismissing(loads_kw)
            if ismissing(critical_loads_kw)
                critical_loads_kw = critical_load_pct * loads_kw
            end
            return new(
                loads_kw,
                year,
                critical_loads_kw,
                loads_kw_is_net,
                critical_loads_kw_is_net
            )

        elseif !ismissing(doe_reference_name)
            # NOTE: must use year that starts on Sunday with DOE reference doe_ref_profiles
            if year != 2017
                @warn "Changing ElectricLoad.year to 2017 because DOE reference profiles start on a Sunday."
            end
            year = 2017
            loads_kw = BuiltInElectricLoad(city, doe_reference_name, latitude, longitude, annual_kwh=annual_kwh)
            if ismissing(critical_loads_kw)
                critical_loads_kw = critical_load_pct * loads_kw
            end
            return new(
                loads_kw,
                year,
                critical_loads_kw,
                loads_kw_is_net,
                critical_loads_kw_is_net
            )

        else
            error("Cannot construct ElectricLoad. You must provide either loads_kw, [doe_reference_name, city],
                  or [doe_reference_name, latitude, longitude].")
        end
    end
end

function BuiltInElectricLoad(
    city::String,
    buildingtype::String,
    latitude::Float64,
    longitude::Float64;
    annual_kwh::Union{Float64,Nothing}=nothing
    )
    lib_path = joinpath(dirname(@__FILE__), "..", "..", "data")
    annual_loads = Dict(
        "Albuquerque" => Dict(
            "fastfoodrest" => 193235,
            "fullservicerest" => 367661,
            "hospital" => 8468546,
            "largehotel" => 2407649,
            "largeoffice" => 6303595,
            "mediumoffice" => 884408,
            "midriseapartment" => 269734,
            "outpatient" => 1678720,
            "primaryschool" => 1070908,
            "retailstore" => 505417,
            "secondaryschool" => 2588879,
            "smallhotel" => 755373,
            "smalloffice" => 87008,
            "stripmall" => 497132,
            "supermarket" => 1947654,
            "warehouse" => 228939,
            "flatload" => 500000
        ),
        "Atlanta" => Dict(
            "fastfoodrest" => 197467,
            "fullservicerest" => 353750,
            "hospital" => 9054747,
            "largehotel" => 2649819,
            "largeoffice" => 6995864,
            "mediumoffice" => 929349,
            "midriseapartment" => 285349,
            "outpatient" => 1672434,
            "primaryschool" => 1128702,
            "retailstore" => 543340,
            "secondaryschool" => 2849901,
            "smallhotel" => 795777,
            "smalloffice" => 90162,
            "stripmall" => 529719,
            "supermarket" => 2092966,
            "warehouse" => 223009,
            "flatload" => 500000
        ),
        "Baltimore" => Dict(
            "fastfoodrest" => 192831,
            "fullservicerest" => 341893,
            "hospital" => 8895223,
            "largehotel" => 2534272,
            "largeoffice" => 6836130,
            "mediumoffice" => 945425,
            "midriseapartment" => 273225,
            "outpatient" => 1623103,
            "primaryschool" => 1077312,
            "retailstore" => 510257,
            "secondaryschool" => 2698987,
            "smallhotel" => 767538,
            "smalloffice" => 86112,
            "stripmall" => 504715,
            "supermarket" => 2018760,
            "warehouse" => 229712,
            "flatload" => 500000
        ),
        "Boulder" => Dict(
            "fastfoodrest" => 189092,
            "fullservicerest" => 334005,
            "hospital" => 8281865,
            "largehotel" => 2313151,
            "largeoffice" => 6127030,
            "mediumoffice" => 884726,
            "midriseapartment" => 255428,
            "outpatient" => 1621950,
            "primaryschool" => 1018424,
            "retailstore" => 504256,
            "secondaryschool" => 2441588,
            "smallhotel" => 736174,
            "smalloffice" => 84900,
            "stripmall" => 495018,
            "supermarket" => 1956244,
            "warehouse" => 243615,
            "flatload" => 500000
        ),
        "Chicago" => Dict(
            "fastfoodrest" => 189558,
            "fullservicerest" => 333659,
            "hospital" => 8567087,
            "largehotel" => 2402021,
            "largeoffice" => 6369028,
            "mediumoffice" => 972772,
            "midriseapartment" => 265528,
            "outpatient" => 1587062,
            "primaryschool" => 1045477,
            "retailstore" => 513106,
            "secondaryschool" => 2568086,
            "smallhotel" => 759657,
            "smalloffice" => 86224,
            "stripmall" => 506886,
            "supermarket" => 2025507,
            "warehouse" => 245750,
            "flatload" => 500000
        ),
        "Duluth" => Dict(
            "fastfoodrest" => 183713,
            "fullservicerest" => 318867,
            "hospital" => 8134328,
            "largehotel" => 2231678,
            "largeoffice" => 6036003,
            "mediumoffice" => 1032533,
            "midriseapartment" => 256393,
            "outpatient" => 1534322,
            "primaryschool" => 982163,
            "retailstore" => 532503,
            "secondaryschool" => 2333466,
            "smallhotel" => 752284,
            "smalloffice" => 83759,
            "stripmall" => 500979,
            "supermarket" => 1980986,
            "warehouse" => 256575,
            "flatload" => 500000
        ),
        "Fairbanks" => Dict(
            "fastfoodrest" => 182495,
            "fullservicerest" => 314760,
            "hospital" => 7899166,
            "largehotel" => 2181664,
            "largeoffice" => 5956232,
            "mediumoffice" => 1267132,
            "midriseapartment" => 271840,
            "outpatient" => 1620270,
            "primaryschool" => 986128,
            "retailstore" => 573411,
            "secondaryschool" => 2344790,
            "smallhotel" => 831480,
            "smalloffice" => 86614,
            "stripmall" => 545421,
            "supermarket" => 2033295,
            "warehouse" => 285064,
            "flatload" => 500000
        ),
        "Helena" => Dict(
            "fastfoodrest" => 185877,
            "fullservicerest" => 325263,
            "hospital" => 8068698,
            "largehotel" => 2246239,
            "largeoffice" => 6003137,
            "mediumoffice" => 930630,
            "midriseapartment" => 252659,
            "outpatient" => 1568262,
            "primaryschool" => 994496,
            "retailstore" => 534933,
            "secondaryschool" => 2357548,
            "smallhotel" => 729797,
            "smalloffice" => 84219,
            "stripmall" => 503504,
            "supermarket" => 1969137,
            "warehouse" => 252245,
            "flatload" => 500000
        ),
        "Houston" => Dict(
            "fastfoodrest" => 210283,
            "fullservicerest" => 383987,
            "hospital" => 9634661,
            "largehotel" => 3050370,
            "largeoffice" => 7539308,
            "mediumoffice" => 972535,
            "midriseapartment" => 335063,
            "outpatient" => 1756541,
            "primaryschool" => 1258146,
            "retailstore" => 589419,
            "secondaryschool" => 3421024,
            "smallhotel" => 863952,
            "smalloffice" => 98508,
            "stripmall" => 577987,
            "supermarket" => 2225265,
            "warehouse" => 221593,
            "flatload" => 500000
        ),
        "LosAngeles" => Dict(
            "fastfoodrest" => 188857,
            "fullservicerest" => 352240,
            "hospital" => 8498389,
            "largehotel" => 2458786,
            "largeoffice" => 6642784,
            "mediumoffice" => 846742,
            "midriseapartment" => 248028,
            "outpatient" => 1565008,
            "primaryschool" => 1095263,
            "retailstore" => 486188,
            "secondaryschool" => 2584380,
            "smallhotel" => 751880,
            "smalloffice" => 86655,
            "stripmall" => 491972,
            "supermarket" => 1935886,
            "warehouse" => 182085,
            "flatload" => 500000
        ),
        "Miami" => Dict(
            "fastfoodrest" => 224494,
            "fullservicerest" => 448713,
            "hospital" => 10062043,
            "largehotel" => 3437188,
            "largeoffice" => 8002063,
            "mediumoffice" => 1021224,
            "midriseapartment" => 424956,
            "outpatient" => 1929148,
            "primaryschool" => 1426635,
            "retailstore" => 635086,
            "secondaryschool" => 4074081,
            "smallhotel" => 972090,
            "smalloffice" => 108279,
            "stripmall" => 675793,
            "supermarket" => 2260929,
            "warehouse" => 202082,
            "flatload" => 500000
        ),
        "Minneapolis" => Dict(
            "fastfoodrest" => 188368,
            "fullservicerest" => 330920,
            "hospital" => 8425063,
            "largehotel" => 2378872,
            "largeoffice" => 6306693,
            "mediumoffice" => 1005875,
            "midriseapartment" => 267383,
            "outpatient" => 1582701,
            "primaryschool" => 1022667,
            "retailstore" => 539203,
            "secondaryschool" => 2498647,
            "smallhotel" => 774571,
            "smalloffice" => 85921,
            "stripmall" => 511567,
            "supermarket" => 2034650,
            "warehouse" => 249332,
            "flatload" => 500000
        ),
        "Phoenix" => Dict(
            "fastfoodrest" => 216088,
            "fullservicerest" => 389739,
            "hospital" => 9265786,
            "largehotel" => 2990053,
            "largeoffice" => 7211666,
            "mediumoffice" => 1004988,
            "midriseapartment" => 378378,
            "outpatient" => 1849358,
            "primaryschool" => 1289084,
            "retailstore" => 593924,
            "secondaryschool" => 3503727,
            "smallhotel" => 881843,
            "smalloffice" => 104583,
            "stripmall" => 590954,
            "supermarket" => 2056195,
            "warehouse" => 241585,
            "flatload" => 500000
        ),
        "SanFrancisco" => Dict(
            "fastfoodrest" => 183953,
            "fullservicerest" => 317124,
            "hospital" => 7752817,
            "largehotel" => 2206880,
            "largeoffice" => 6085403,
            "mediumoffice" => 792199,
            "midriseapartment" => 229671,
            "outpatient" => 1394447,
            "primaryschool" => 1009369,
            "retailstore" => 449025,
            "secondaryschool" => 2327074,
            "smallhotel" => 698095,
            "smalloffice" => 78132,
            "stripmall" => 455802,
            "supermarket" => 1841655,
            "warehouse" => 185889,
            "flatload" => 500000
        ),
        "Seattle" => Dict(
            "fastfoodrest" => 184142,
            "fullservicerest" => 318741,
            "hospital" => 7912504,
            "largehotel" => 2212410,
            "largeoffice" => 6019271,
            "mediumoffice" => 878390,
            "midriseapartment" => 237242,
            "outpatient" => 1434195,
            "primaryschool" => 983498,
            "retailstore" => 455854,
            "secondaryschool" => 2282972,
            "smallhotel" => 693921,
            "smalloffice" => 79716,
            "stripmall" => 460449,
            "supermarket" => 1868973,
            "warehouse" => 210300,
            "flatload" => 500000
        ),
    )
    default_buildings = [
        "FastFoodRest",
        "FullServiceRest",
        "Hospital",
        "LargeHotel",
        "LargeOffice",
        "MediumOffice",
        "MidriseApartment",
        "Outpatient",
        "PrimarySchool",
        "RetailStore",
        "SecondarySchool",
        "SmallHotel",
        "SmallOffice",
        "StripMall",
        "Supermarket",
        "Warehouse",
        "FlatLoad",
    ]
    if !(buildingtype in default_buildings)
        error("buildingtype $(buildingtype) not in $(default_buildings).")
    end

    if isempty(city)
        city = find_ashrae_zone_city(latitude, longitude)
    end

    if isnothing(annual_kwh)
        annual_kwh = annual_loads[city][lowercase(buildingtype)]
    end
    # TODO implement BuiltInElectricLoad scaling based on monthly_totals_kwh

    profile_path = joinpath(lib_path, string("Load8760_norm_" * city * "_" * buildingtype * ".dat"))
    normalized_profile = vec(readdlm(profile_path, '\n', Float64, '\n'))

    load = [annual_kwh * ld for ld in normalized_profile]
end

function find_ashrae_zone_city(lat, lon)::String
    file_path = joinpath(dirname(@__FILE__), "..", "..", "data", "climate_cities.shp")
    table = Shapefile.Table(file_path)
    geoms = Shapefile.shapes(table)
    # TODO following for loop is relatively slow
    for (row, geo) in enumerate(geoms)
        g = length(geo.points)
        nodes = zeros(g, 2)
        edges = zeros(g, 2)
        for (i,p) in enumerate(geo.points)
            nodes[i,:] = [p.x, p.y]
            edges[i,:] = [i, i+1]
        end
        edges[g, :] = [g, 1]
        edges = convert(Array{Int64,2}, edges)
        # shapefiles have longitude as x, latitude as y
        if inpoly2([lon, lat], nodes, edges)[1]
            return table.city[row]
        end
        GC.gc()
    end
    @info "Could not find latitude/longitude in U.S. Using geometrically nearest city."
    cities = [
        (city="Miami", lat=25.761680, lon=-80.191790),
        (city="Houston", lat=29.760427, lon=-95.369803),
        (city="Phoenix", lat=33.448377, lon=-112.074037),
        (city="Atlanta", lat=33.748995, lon=-84.387982),
        (city="LasVegas", lat=36.1699, lon=-115.1398),
        (city="LosAngeles", lat=34.052234, lon=-118.243685),
        (city="SanFrancisco", lat=37.3382, lon=-121.8863),
        (city="Baltimore", lat=39.290385, lon=-76.612189),
        (city="Albuquerque", lat=35.085334, lon=-106.605553),
        (city="Seattle", lat=47.606209, lon=-122.332071),
        (city="Chicago", lat=41.878114, lon=-87.629798),
        (city="Boulder", lat=40.014986, lon=-105.270546),
        (city="Minneapolis", lat=44.977753, lon=-93.265011),
        (city="Helena", lat=46.588371, lon=-112.024505,),
        (city="Duluth", lat=46.786672, lon=-92.100485),
        (city="Fairbanks", lat=59.0397, lon=-158.4575),
    ]
    min_distance = 0.0
    nearest_city = ""
    for (i, c) in enumerate(cities)
        distance = sqrt((lat - c.lat)^2 + (lon - c.lon)^2)
        if i == 1
            min_distance = distance
            nearest_city = c.city
        elseif distance < min_distance
            min_distance = distance
            nearest_city = c.city
        end
    end
    return nearest_city
end

max_stars_repo_name: NLaws/REoptLite

max_issues_repo_name: NLaws/REoptLite

max_forks_repo_name: NLaws/REoptLite


text:
@testset "344.reverse-string.jl" begin
    s = ['h', 'e', 'l', 'l', 'o']
    reverse_string!(s)
    @test s == ['o', 'l', 'l', 'e', 'h']

    s = ['H', 'a', 'n', 'n', 'a', 'h']
    reverse_string!(s)
    @test s == ['h', 'a', 'n', 'n', 'a', 'H']
end


text:
@inline function get_first_bit_sign(x::VInt16, ::Val{N}) where N
    mysign(x << (16 - N - 1))
end
@inline function get_second_bit_sign(x::VInt16, ::Val{N}) where N
    mysign(x << (16 - N - 2))
end
@inline function get_first_bit_sign(x::VInt32, ::Val{N}) where N
    mysign(x << (32 - N - 1))
end
@inline function get_second_bit_sign(x::VInt32, ::Val{N}) where N
    mysign(x << (32 - N - 2))
end
@inline function get_quarter_angle(x, ::Val{N}) where N
    x & (one(x) << N - one(x)) - one(x) << (N - 1)
end
@inline mysign(x) = vifelse(x >= zero(x), one(x), -one(x))

"""
$(SIGNATURES)

Fixed point cos
"""
@inline function fpcos(phase::P, bits::Val{N}) where {N, P <: Union{VInt16, VInt32}}
    first_bit_sign = get_first_bit_sign(phase, bits)
    second_bit_sign = get_second_bit_sign(phase, bits)
    quarter_angle = get_quarter_angle(phase, bits)
    A = calc_A(quarter_angle, bits)
    B = calc_B(quarter_angle, bits)

    second_bit_sign * (first_bit_sign * A + B)
end

"""
$(SIGNATURES)

Fixed point sin
"""
@inline function fpsin(phase::P, bits::Val{N}) where {N, P <: Union{VInt16, VInt32}}
    first_bit_sign = get_first_bit_sign(phase, bits)
    second_bit_sign = get_second_bit_sign(phase, bits)
    quarter_angle = get_quarter_angle(phase, bits)
    A = calc_A(quarter_angle, bits)
    B = calc_B(quarter_angle, bits)

    second_bit_sign * (A - first_bit_sign * B)
end

"""
$(SIGNATURES)

Fixed point sin and cos
"""
@inline function fpsincos(x::P, bits::Val{N}) where {N, P <: Union{VInt16, VInt32}}
    (fpsin(x, bits), fpcos(x, bits))
end


text:
# using PyPlot
using MatrixColorings

#####
# Example evaluating a numerical a Jacobian.  Uses the Brusselator in IVPtestsuite.jl
####

const N = 10^4 # active gridpoints
dof = 2N # degrees of freedom (boundary points are not free)
dx = 1/(N+1)
x = dx:dx:1-dx
T = Float64

const alpha = 1/50
const gamma = alpha/dx^2

# BC
const one_ = one(T)
const three = 3*one_
ubc1(t) = one_
ubc2(t) = one_
vbc1(t) = three
vbc2(t) = three

# immutable Pair
#     i::Int
#     j::Int
# end
# inds(i) = Pair(2i-1, 2i)
inds(i) = (2i-1, 2i)
function fn!(t,y, dydt)
    # the ode function dydt = f(t,y)
    #
    # Note the u and v are staggered:
    # u = y[1:2:end], v=y[2:2:end]

    # setup for first step
    i = 1
    iu, iv = inds(i)

    ui0 = ubc1(t); vi0 = vbc1(t) # BC
    ui1 = y[iu];   vi1 = y[iv]
    ui2 = y[iu+2]; vi2 = y[iv+2]
    while i<=N
        dydt[iu] = 1 + ui1^2*vi1 - 4*ui1 + gamma*(ui0 - 2*ui1 + ui2)
        dydt[iv] = 3*ui1 - ui1^2*vi1     + gamma*(vi0 - 2*vi1 + vi2)
        # setup for next step:
        i += 1
        iu, iv = inds(i)
        ui0 = ui1; vi0 = vi1
        ui1 = ui2; vi1 = vi2
        if i<N
            ui2 = y[iu+2]
            vi2 = y[iv+2]
        else # on last step do BC
            ui2 = ubc2(t)
            vi2 = vbc2(t)
        end
    end
    return nothing
end
fn!(;T_::Type=T, dof_=dof) = zeros(T_,dof_) # returns array for storage for dydt

# Analytic Jacobian: it is a banded matrix with upper and lower
# bandwidth 2, c.f. W&H p.148, but Julia does not support banded
# matrices yet.  Thus use a sparse matrix.
function jac!(t,y,dfdy)
    # The Jacobian of fn
    ii = 1 # direct index into dfdy.nzval
    for i=1:N # iterate over double-columns
        iu, iv = inds(i)
        ui = y[iu]
        vi = y[iv]
        # iu: column belonging to d/du(i)
        if iu-2>=1
            # du'(i-1)/du(i)
            dfdy.nzval[ii] = gamma
            ii +=1
        end
        if iv-2>=1
            # dv'(i-1)/du(i)==0
            dfdy.nzval[ii] = 0
            ii +=1
        end
        # du'(i)/du(i)
        dfdy.nzval[ii]   = 2*ui*vi - 4 -2*gamma
        ii += 1
        # dv'(i)/du(i)
        dfdy.nzval[ii] = 3 - 2*ui*vi
        ii +=1
        if  iu+2<=2N
            # du'(i+1)/du(i)
            dfdy.nzval[ii] = gamma
            ii +=1
        end

        # iv: column belonging to d/dv(i)
        if iv-2>=1
            # dv'(i-1)/dv(i)
            dfdy.nzval[ii] = gamma
            ii +=1
        end
        # du'(i)/dv(i)
        dfdy.nzval[ii] = ui^2
        ii +=1
        # dv'(i)/dv(i)
        dfdy.nzval[ii]   = -ui^2 -2*gamma
        ii +=1
        if  iu+2<=2N
            # du'(i+1)/dv(i)==0
            dfdy.nzval[ii] = 0
            ii +=1
        end
        if iv+2<=2N
            # dv'(i+1)/dv(i)
            dfdy.nzval[ii] = gamma
            ii +=1
        end
    end
    return nothing
end
function jac!(;T_::Type=T, dof_=dof)
    # returns an initialized sparse matrix
    B = (ones(T_,dof_-2), ones(T_,dof_-1), ones(T_,dof_), ones(T_,dof_-1), ones(T_,dof_-2))
    spdiagm(B, (-2,-1,0,1,2))
end

####
# Numerical Jacobians
####

# Non-coloring method:
######################
# Copied & adapted from DASSL:
# generate a function that computes approximate Jacobian using forward
# finite differences
if VERSION < v"0.4.0-dev"
    # TODO: remove these two for 0.4
    rowvals(S::SparseMatrixCSC) = S.rowval
    nzrange(S::SparseMatrixCSC, col::Integer) = S.colptr[col]:(S.colptr[col+1]-1)
end
function numjac!(t, y, jac)
    # updates jac in-place

    edelta  = eps(4.0)*1e7 # trial and error
    dof = length(y)
    f1 = similar(y)
    tmpy = similar(y) # d+dy
    f0 = similar(y)
    fn!(t, y, f0)
    for i=1:length(y)
        copy!(tmpy, y)
        tmpy[i] += edelta
        fn!(t, tmpy, f1)
        for j=1:dof
            tmpy[j] = (f1[j]-f0[j])/edelta
        end
        for nz in nzrange(jac,i)
            nonzeros(jac)[nz] = tmpy[rowvals(jac)[nz]]
        end
        # which, for large matrices, is much faster than:
        # jac[:,i] = (f1-f0)/edelta
    end
    return nothing
end

# Coloring method:
function numjac_c!(t, y, jac, S)
    # updates jac in-place
    nc = size(S,2) # number of colors

    edelta  = eps(4.0)*1e7 # trial and error
    dof = length(y)
    f0 = similar(y)
    fn!(t, y, f0)
    ddy = deepcopy(y) # d+dy
    f1 = similar(y)   # f(d+dy)
    for c in 1:nc
        add_delta!(ddy, S, c, edelta)  # adds a delta to all indices of color c
        fn!(t, ddy, f1)
        for j=1:dof
            f1[j] = (f1[j]-f0[j])/edelta
        end
        recover_jac!(jac, f1, S, c) # recovers the columns of jac which correspond c
        remove_delta!(ddy, S, c, edelta) # remove delta again
    end
    return nothing
end

function numjac_c!(t, y, jac)
    # updates jac in-place, uses its sparsity structure for coloring
    S = color_jac(jac) # seed matrix
    numjac_c!(t, y, jac, S)
    return nothing
end

####
# Compare
####

t = 0.0 # not dependent on t
# initialize
dydt = fn!()
dfdy_ana = jac!()
dfdy_num = jac!()
dfdy_num_c = jac!()

y = fn!()
#IC
u0(x) = 1 + 0.5*sin(2π*x)
v0(x) = 3 * ones(length(x))
y[1:2:2N] = u0(x)
y[2:2:2N] = v0(x)
# warm ups:
@time 1
jac!(t, y, dfdy_ana)
numjac!(t, y, dfdy_num)

#ys = (y, rand(T, dof)*3+0.5)
for y in (y,)
    # calculate, time and compare
    fn!(t, y, dydt)
    println("@time for analytic:")
    @time jac!(t, y, dfdy_ana)
    println("@time for numeric with no coloring:")
    @time numjac!(t, y, dfdy_num)
    println("@time for numeric with coloring:")
    @time numjac_c!(t, y, dfdy_num_c)
    S = color_jac(dfdy_num_c) # seed matrix
    println("@time for numeric with coloring with provided seed matrix:")
    @time numjac_c!(t, y, dfdy_num_c,S)

    @show maximum(dfdy_ana-dfdy_num)
    @show maximum(dfdy_ana-dfdy_num_c)
end


text:
"""

# m_separation

$(SIGNATURES)

Computes the m_separation between 2 sets of nodes conditioned on a third set.

### Required arguments
```julia
m_separation(
* `d::DAG`                             : DAG
* `f::SymbolList`                      : First vertex or set
*  s::SymbolList`                      : Second vertex or set
)
```

### Keyword arguments
```julia
* `c::SymbolListOrNothing=nothing`     : Conditioning set
* `debug=false`                        : Trace execution
```

### Returns
```julia
* `res::Bool`                          : Boolean result of test
```

# Extended help

### Example

### m_separation between mechanics and statistics, conditioning on algebra
```julia
using StructuralCausalModels, CSV

df = DataFrame!(CSV.File(scm_path("..", "data", "marks.csv"));

d = OrderedDict(
  :mechanics => [:vectors, :algebra],
  :vectors => [:algebra],
  :analysis => [:algebra],
  :statistics => [:algebra, :analysis]
);

dag = DAG("marks", d, df);
m_separation(marks, [:statistics], [:mechanics]; c=[:algebra]))
```
### Acknowledgements

Original author:                       Giovanni M. Marchetti

Translated to Julia:                   Rob J Goedman

### License

The R package ggm is licensed under License: GPL-2.

The Julia translation is licenced under: MIT.

Part of the API, exported.
"""
function m_separation(d::DAG, f::SymbolList, s::SymbolList;
  c::SymbolListOrNothing=nothing, debug=false)

  if typeof(f) == Symbol
    f = [f]
  end
  if typeof(s) == Symbol
    s = [s]
  end
  if typeof(c) == Symbol
    c = [c]
  end

  if isnothing(c)
    m = setdiff(d.vars, vcat(f, s))
    debug && println("m = $m")
    ar = maximize(ribbon_graph(d, m=m))
    debug && println("ar = $ar")
  else
    m = setdiff(d.vars, vcat(f, s, c))
    debug && println("m = $m, c = $c")
    ar = maximize(ribbon_graph(d, m=m, c=c))
    debug && println("ar = $ar")
  end

  all(Array(ar[f, s]) + Array(ar[s, f]) .== 0)
end

export
  m_separation


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

mutable struct PolicyDefinitionReference <: SwaggerModel
    policyDefinitionId::Any # spec type: Union{ Nothing, String } # spec name: policyDefinitionId
    parameters::Any # spec type: Union{ Nothing, ParameterValues } # spec name: parameters
    policyDefinitionReferenceId::Any # spec type: Union{ Nothing, String } # spec name: policyDefinitionReferenceId
    groupNames::Any # spec type: Union{ Nothing, Vector{String} } # spec name: groupNames

    function PolicyDefinitionReference(;policyDefinitionId=nothing, parameters=nothing, policyDefinitionReferenceId=nothing, groupNames=nothing)
        o = new()
        validate_property(PolicyDefinitionReference, Symbol("policyDefinitionId"), policyDefinitionId)
        setfield!(o, Symbol("policyDefinitionId"), policyDefinitionId)
        validate_property(PolicyDefinitionReference, Symbol("parameters"), parameters)
        setfield!(o, Symbol("parameters"), parameters)
        validate_property(PolicyDefinitionReference, Symbol("policyDefinitionReferenceId"), policyDefinitionReferenceId)
        setfield!(o, Symbol("policyDefinitionReferenceId"), policyDefinitionReferenceId)
        validate_property(PolicyDefinitionReference, Symbol("groupNames"), groupNames)
        setfield!(o, Symbol("groupNames"), groupNames)
        o
    end
end # type PolicyDefinitionReference

const _property_map_PolicyDefinitionReference = Dict{Symbol,Symbol}(Symbol("policyDefinitionId")=>Symbol("policyDefinitionId"), Symbol("parameters")=>Symbol("parameters"), Symbol("policyDefinitionReferenceId")=>Symbol("policyDefinitionReferenceId"), Symbol("groupNames")=>Symbol("groupNames"))
const _property_types_PolicyDefinitionReference = Dict{Symbol,String}(Symbol("policyDefinitionId")=>"String", Symbol("parameters")=>"ParameterValues", Symbol("policyDefinitionReferenceId")=>"String", Symbol("groupNames")=>"Vector{String}")
Base.propertynames(::Type{ PolicyDefinitionReference }) = collect(keys(_property_map_PolicyDefinitionReference))
Swagger.property_type(::Type{ PolicyDefinitionReference }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_PolicyDefinitionReference[name]))}
Swagger.field_name(::Type{ PolicyDefinitionReference }, property_name::Symbol) =  _property_map_PolicyDefinitionReference[property_name]

function check_required(o::PolicyDefinitionReference)
    (getproperty(o, Symbol("policyDefinitionId")) === nothing) && (return false)
    true
end

function validate_property(::Type{ PolicyDefinitionReference }, name::Symbol, val)
end


text:
"""
    colourprint(str; foreground = :light_green, bold = false)

Auxiliary function for printing out text in colour in the terminal.
Available foreground colours are specified by the Crayons.jl package and are:
black, red, green, yellow, blue, magenta, cyan, light_gray, default, dark_gray, light_red, light_green, light_yellow, light_blue, light_magenta, light_cyan and white.
"""
function colourprint(str; foreground = :light_green, bold = false)

    print(Crayon(foreground = foreground, bold = bold), @sprintf("%s", str), Crayon(reset = true))

end


text:
# csg.jl

import Base.hypot

function volume(f, x_min,y_min,z_min,x_max,y_max,z_max, scale)
	x_rng = x_max - x_min
	y_rng = y_max - y_min
	z_rng = z_max - z_min

	nx, ny, nz = int(scale*x_rng), int(scale*y_rng), int(scale*z_rng)

	vol = zeros(Float64, (nx,ny,nz))

	for i = 1:nx, j = 1:ny, k = 1:nz
		x = x_min + x_rng*(i/nx)
		y = y_min + y_rng*(j/ny)
		z = z_min + z_rng*(k/nz)

		vol[i,j,k] = f(x,y,z)
	end

	vol
end
export volume

function hypot(x,y,z)
	sqrt(x^2 + y^2 + z^2)
end

function sphere(x,y,z, sx,sy,sz,sr)
	hypot(x-sx, y-sy, z-sz) - sr
end
export sphere

function cylinderX(x,y,z, cy,cz,cr,cxmin,cxmax)
	vr = hypot(y-cy,z-cz) - cr
	vx = max(cxmin-x, x-cxmax)
	max(vr,vx)
end
export cylinderX

function cylinderY(x,y,z, cx,cz,cr,cymin,cymax)
	vr = hypot(x-cx,z-cz) - cr
	vy = max(cymin-y, y-cymax)
	max(vr,vy)
end
export cylinderY

function cylinderZ(x,y,z, cx,cy,cr,czmin,czmax)
	vr = hypot(x-cx,y-cy) - cr
	vz = max(czmin-z, z-czmax)
	max(vr,vz)
end
export cylinderZ

function box(x,y,z, x_min,y_min,z_min,x_max,y_max,z_max)
	vx = max(x_min-x, x-x_max)
	vy = max(y_min-y, y-y_max)
	vz = max(z_min-z, z-z_max)
	max(vx, vy, vz)
end
export box


text:
"""
    SeisWindowPatch(in,out;<keyword arguments>)
Window a seis file using header words.

# Arguments
- `in::String`: filename of input
- `out::String`: filename of output
- `key=[]`: note that windowing along the time axis is achieved by using the key "t".
- `minval=[]`
- `maxval=[]`
- `it_nt=9e9`

*Credits: AS, FC, 2017*
"""
function SeisWindowPatch(in::String,out::String;key=[],minval=[],maxval=[],it_nt=9e9)
  println("processing patch ",out)
  DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
  extent = ReadTextHeader(in)
  tmin = extent.o1
  tmax = extent.o1 + extent.d1*(extent.n1-1)
  for ikey=1:length(key)
	   if key[ikey] == "t"
	      tmin = minval[ikey]
	      tmax = maxval[ikey]
	   end
  end

  #itmin is the sample at which tmin is saved
  itmin = convert(Int32,round((tmin - extent.o1)/extent.d1) + 1)

  if (itmin < 1)
	   itmin = 1
  end
  itmax = convert(Int32,round((tmax - extent.o1)/extent.d1) + 1)
  if (itmax) > extent.n1
	   itmax = extent.n1
  end

  SeisWindowHeadersPatch(in, out; key=key, minval=minval, maxval=maxval, tmin=tmin,
                      tmax=tmax,it_nt=it_nt)
  FetchTracesPatch(in,out;itmin=itmin,itmax=itmax ,key=key, minval=minval, maxval=maxval)
  tmp = join(["tmp_SeisWindow_",string(round(Int,rand()*100000))])

 #>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 #@compat SeisProcessHeaders(out, tmp, [UpdateHeader],[Dict(:itmin=>itmin,:itmax=>itmax)])

  itrace_in = 1
  itrace_out = 1
  ntrace=1000000
  NX = GetNumTraces(out)
  while itrace_in <= NX
    nx = NX - itrace_in + 1
	  ntrace = nx > ntrace ? ntrace : nx
	  h1 = SeisReadHeaders(out, group="some", itrace=itrace_in,
                                 ntrace=ntrace)
    ext=ReadTextHeader(out)

	  num_traces_in = length(h1)
	  h2=UpdateHeaderPatch(h1;itrace=itrace_in,itmin=itmin,itmax=itmax)
	  h1 = copy(h2)
	  num_traces_out = length(h1)
	  SeisWriteHeaders(tmp, h1, itrace=itrace_out)
    itrace_in += num_traces_in
	  itrace_out += num_traces_out
  end

  filename_h_tmp = join([DATAPATH tmp "@headers@"])
  filename_h_out = join([DATAPATH out "@headers@"])
  cp(filename_h_tmp,filename_h_out,force=true);
  rm(filename_h_tmp);
    #rm(tmp);
end

########################################################################################
#####################################################################################
function FetchTracesPatch(in::String, out::String; ntrace=500, itmin=round(Int,1),
    itmax=round(Int,9999999999),key=[], minval=[], maxval=[])

    NX = GetNumTraces(out)

    itrace = round(Int,1)
    filename_data_in = ParseDataName(in)

    DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
    filename_data_out = join([DATAPATH out "@data@"])

    stream_in  = open(filename_data_in)

    stream_out  = open(filename_data_out,"w")
    extent = ReadTextHeader(in)
    nt = extent.n1

    if itmin < 1
      itmin = round(Int,1)
    end
    if itmax > nt
      itmax = round(Int,nt)
    end

    nhead = length(fieldnames(Header))

    while itrace <= NX
      if (itrace > 1)
	       stream_out = open(filename_data_out,"a")
	    end
      nx = NX - itrace + 1

      ntrace = nx > ntrace ? ntrace : nx
      h  = SeisReadHeaders(out,group="some",itrace=itrace,ntrace=ntrace)
      d = zeros(Float32,itmax-itmin+1,ntrace)

      SeekTracesPatch!(d,stream_in,h,itmin,itmax,nt,round(Int,ntrace))

    write(stream_out,d)

    close(stream_out)
    itrace += ntrace
  end
  close(stream_in)
end

###############################################################################

function SeekTracesPatch!(d::AbstractArray{T,2}, stream_in::IOStream,
                        h::Array{Header,1},itmin,itmax,nt,ntrace) where T

    d1 = zeros(Float32,nt)
    for ix = 1 : ntrace
	     position = 4*nt*(h[ix].tracenum-1)
	     seek(stream_in,position)
         d1 = read!(stream_in,Array{Float32}(undef, nt))
         for it = itmin : itmax
	        d[it - itmin + 1,ix] = d1[it]
	     end
    end
    nothing
end

##############################################################################
#function SeekTracesPatch{T}(d::AbstractArray{T,2}, h::Array{Header,1},out::String, stream_in::IOStream,stream_out::IOStream, itmin,itmax,nt,ntrace;key=[],minval=[],maxval=[])
# println(itmin," ",itmax," ",nt)
#    d1 = zeros(Float32,nt)
#extent_out=ReadTextHeader(out)
#nx2=extent_out.n3
#nx3=extent_out.n4
#nx4=extent_out.n5
#    for ix = 1 : ntrace
#	position = 4*nt*(h[ix].tracenum-1)
#	seek(stream_in,Int(position))

#	d1 = read(stream_in,Float32,nt)
#	for it = itmin : itmax
#	    d[it - itmin + 1,ix] = d1[it]

#	end

#        itrace_out=(h[ix].imx-minval[2])*nx2*nx3*nx4+(h[ix].imy-minval[3])*nx3*nx4+(h[ix].ihx-minval[4])*nx4+h[ix].ihy-minval[5]+1

#        position_out=4*nt*(itrace_out-1)
#println(out," ",h[ix].imx," ",minval[2]," ",itrace_out)
#        seek(stream_out,Int(position_out))
#        write(stream_out,convert(Array{Float32,1},d[1:itmax-itmin+1,ix]))
#    end
#    nothing
#end

##############################################################################
##############################################################################
function UpdateHeaderPatch(h;itrace=1,itmin=1,itmax=9e9)
    ot = (itmin-1)*h[1].d1 + h[1].o1

    nt = itmax - itmin + 1
    for j = 1 : length(h)
	h[j].o1 = ot
	h[j].n1 = nt
	h[j].tracenum = j + itrace - 1
    end
    return h
end


text:
# This file was generated, do not modify it. # hide
rng = StableRNG(566)
train, test = partition(eachindex(y), 0.7, shuffle=true, rng=rng)
test[1:3]


text:
"""
(D)_{ij} = ∑_k (x_k^{(i)} - x_k^{(j)})^2
"""

function square_metric(x)
	N = size(x, 2)
	D = zeros(eltype(x), N, N)
	@inbounds for i in 1:N
		for j in i+1:N
			dx = norm(view(x, :, i) .- view(x, :, j))
			D[i,j] = D[j,i] = dx * dx
		end
	end
	D
end

function square_metric_back(x, sym_dD)
	dD = 0.5*(sym_dD + sym_dD')  # Cholesky decomposition only use lower triangle, reconstruct the true dD
	m, n = size(x)
	dx = zeros(eltype(x), m, n)
  @inbounds for i in 1:m
    for j in 1:n
      for k in 1:n
        dx[i,j] += 4*dD[j,k]*(x[i,j]-x[i,k])
      end
    end
  end
  dx
end

function square_metric(x, xo)
	N1 = size(x, 2)
	N2 = size(xo, 2)
	D = zeros(eltype(x), N1, N2)
	@inbounds for i in 1:N1
		for j in 1:N2
			dx = norm(view(x, :, i) .- view(xo, :, j))
			D[i,j] = dx * dx
		end
	end
	D
end

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
"""
    tripcolor(x, y, z, t; zmin=nothing, zmax=nothing)

Draw a pseudocolor plot of unstructured triangular data. `x`, `y`, and `z` are
one-dimensional input arrays of the same length (the number of points). `t` is an integer
array of size `(3,nt)`, where `nt` is the number of triangles. The coordinates of the `i`th
vertex of the `j`th triangle are given by `(x[t[i,j]], y[t[i,j]])`. Function values at
vertices are specified by the `z` array

GR's current colormap (set using `GR.setcolormap`) is used. The colormap bounds can be set
by specifying values for `zmin` and `zmax`. If these parameters are set to `nothing`, then
the minimum and maximum values of the `z` array will be used.
"""
function tripcolor(x, y, z, t; zmin=nothing, zmax=nothing)
    c = rasterize(x, y, z, t; zmin, zmax, yflip=true)
    GR.drawimage(0,1,0,1,size(c,1),size(c,2),c)
end

function rasterize(x, y, z, t; yflip=false, px=nothing, py=nothing, zmin=nothing, zmax=nothing)
    _,_,px_dsp,py_dsp = GR.inqdspsize()
    isnothing(px) && (px = px_dsp)
    isnothing(py) && (py = py_dsp)

    c = zeros(UInt32,px,py)

    # Precompute encoded RGBA values of colormap values
    cmap = zeros(UInt32,256)
    for i=0:255
        # Pixel is encoded as 0xAABBGGRR, GR.inqcolor returns 0x00BBGGRR, so set opaque
        cmap[begin+i] = 0xff000000 | GR.inqcolor(1000+i)
    end

    xmin,xmax = extrema(x)
    ymin,ymax = extrema(y)
    w = xmax-xmin
    h = ymax-ymin
    isnothing(zmin) && (zmin = minimum(z))
    isnothing(zmax) && (zmax = maximum(z))

    nt = size(t,2)
    for it=1:nt
        xt = @views x[t[:,it]]
        yt = @views y[t[:,it]]
        zt = @views z[t[:,it]]

        detT = (yt[2]-yt[3])*(xt[1]-xt[3]) + (xt[3]-xt[2])*(yt[1]-yt[3])

        xmint,xmaxt = extrema(xt)
        ymint,ymaxt = extrema(yt)

        x0 = floor(Int,(xmint-xmin)/w*px)
        x1 = ceil(Int,(xmaxt-xmin)/w*px)
        y0 = floor(Int,(ymint-ymin)/h*py)
        y1 = ceil(Int,(ymaxt-ymin)/h*py)

        # Loop over all pixels in triangle's bounding box
        for i=x0:x1-1, j=y0:y1-1
            # Cartesian coordinates of pixel (i,j) (zero-indexed)
            cx = xmin + i/px*w
            cy = ymin + j/py*h
            # Convert from Cartesian to barycentric coordinates
            λ1 = ((yt[2] -yt[3])*(cx-xt[3]) + (xt[3]-xt[2])*(cy-yt[3]))/detT
            λ2 = ((yt[3] -yt[1])*(cx-xt[3]) + (xt[1]-xt[3])*(cy-yt[3]))/detT
            λ3 = 1 - λ1 - λ2
            tol = 1e-12
            # If pixel is outside of triangle, move on to next pixel
            (λ1>-tol && λ2<1+tol && λ2>-tol && λ2<1+tol && λ3>-tol && λ3<1+tol) || continue
            zval = zt[1]*λ1 + zt[2]*λ2 + zt[3]*λ3
            zval = clamp(zval,zmin,zmax)
            c[i+1,yflip ? py-j : j+1] = cmap[begin+getcolorind(zval,zmin,zmax)-1000]
        end
    end
    c
end


text:
module VarianceReducedSG

using LinearAlgebra
using ProximalOperators
using SparseArrays
using Unrolled
using Printf
using StatsBase

include("gradientstorage.jl")
include("vralgorithms.jl")
include("loggers.jl")

################################################################################
# Main loop
################################################################################
export iterate!

function iterate!(
		alg::VRAlgorithm, niter, logger=NoLog(), nlogs=1; warmstart=false)

	iterperlog = Int(ceil(niter/nlogs))
	nstages = Int(ceil(max(1,nlogs-1)*iterperlog/expectedstagelen(alg)))

	initialize!(logger)
	!warmstart && initialize!(alg)

	log!(logger, alg, 0, 0)

	iter = 0
	for stage = 1:nstages
		stagelen = stageupdate!(alg, iter, stage)

		for _ = 1:stagelen
			iter += 1
			update!(alg, iter, stage)

			iter % iterperlog == 0 && log!(logger, alg, iter, stage)
		end
	end

	finalize!(logger)
end

end # module