text:
# graphstructure.jl
# Authors : Heitor Guimaraes and Luiz Ciafrino
# @brief: Module to handle basic graph manipulation

type simpleGraph
"""
@brief: This composite type is an abstraction of a simple undirected
	graph. It is constructed through a text file contain details about it's
	edges and vertices. To a specific use in search algorithms as BFS and DFS
	you can use a method to transform this abstraction in a more concrete
	data structure as Adjacency List or Adjacency Matrix.
"""
	# Properties of the Graph
	num_vertex::Int64
	num_edge::Int64
	Graph::Union{Array{Array{Int64},1}, Array{UInt8,2}, Array{Array{Array{Float64}}}}
	make_report::Function
	format_graph::ASCIIString
	is_weighted::Bool

	function simpleGraph(filename::ASCIIString, format::ASCIIString = "adjlist")
	"""
	@brief: Constructor of the composite type simpleGraph. This method
		read the text file and construct the Graph as a List with integer
		numbers. Also extract basic information about the graph to make
		a report if the user ask for it.
	@param filename::ASCIIString Name or path of the file contain the graph
		information. This function expect an ASCIIString.
	@param format::ASCIIString Option for specify the type of representation
		for the graph. Options are: 'adjlist' or 'adjmatrix'.
	"""
		# Instantiate the Class
		this = new()

		# Read the text file
		infile = open(filename)
		this.num_vertex = parse(Int64, readline(infile))

		lines = readlines(infile)
		this.num_edge = size(lines)[1]

		this.is_weighted = false
		if length(split(lines[2])) == 3
			this.is_weighted = true
		end

		this.format_graph = format

		# Return the selected structure
		if format == "adjlist"
		"""
		@brief: With this flag you can take the graph abstraction and transform
			into an Adjacency List.
		"""
			if !this.is_weighted
				# Pre-allocate the structure
				this.Graph = Array{Array{Int64}}(this.num_vertex)

				# Initialize with Zeros
				@inbounds @simd for i in 1:this.num_vertex
	   	 			this.Graph[i] = Array{Int64}(0)
	    		end

	    		# Fill the Adjacency List
				@inbounds for i in lines
					vertex = parse(Int64,split(i)[1])
					edge = parse(Int64,split(i)[2])
					push!(this.Graph[vertex], edge)
					push!(this.Graph[edge], vertex)
				end

			else
				# Pre-allocate the structure
				this.Graph = Array{Array{Array{Float64}}}(this.num_vertex)

				# Initialize with Zeros
				@inbounds @simd for i in 1:this.num_vertex
	   	 		this.Graph[i] = Array[Int64[], Float64[]]
	    	end

	    	# Fill the Adjacency List
				@inbounds for i in lines
					vertex = parse(Int64, split(i)[1])
					edge = parse(Int64, split(i)[2])
					weight = parse(Float64, split(i)[3])

					push!(this.Graph[vertex][1], edge)
					push!(this.Graph[edge][1], vertex)

					push!(this.Graph[vertex][2], weight)
					push!(this.Graph[edge][2], weight)
				end
			end

		elseif format == "adjmatrix"
		"""
		@brief: With this flag you can take the graph abstraction and transform
			into an Adjacency Matrix.
		"""
			if !this.is_weighted
				# Pre-allocate the structure and initialize with zeros
				this.Graph = zeros(UInt8, this.num_vertex, this.num_vertex)

				@inbounds for i in lines
					rnum = parse(Int64,split(i)[1])
					cnum = parse(Int64,split(i)[2])
					this.Graph[rnum, cnum] = 1
					this.Graph[cnum, rnum] = 1
				end
			else
				# Pre-allocate the structure and initialize with zeros
				this.Graph = zeros(Float64, this.num_vertex, this.num_vertex)

				@inbounds for i in lines
					rnum = parse(Int64,split(i)[1])
					cnum = parse(Int64,split(i)[2])
					weight = parse(Float64,split(i)[3])
					this.Graph[rnum, cnum] = weight
					this.Graph[cnum, rnum] = weight
				end
			end
		else
			error(""" Invalid format selected, options are:\n
				* adjlist\n
				* adjmatrix\n
				Example: graph = simpleGraph(\"mygraph.txt\", \"adjmatrix\")""")
		end

		this.make_report = function(outfilename::ASCIIString)
		"""
		@brief: Method to make a small report with usefull informations
			about the graph.
		@param outfilename String with the name of the report filename
			in the ASCIIString format.
		"""
			# Get the mean degree and the empirical distance.
			mean_degree, empirical_dist = graph_properties(this)

			# Open the file and write the information on it.
			outfile = open(outfilename, "w")
			write(outfile, string("# n = ", this.num_vertex,"\n"))
			write(outfile, string("# m = ", this.num_edge,"\n"))
			write(outfile, string("# d_medio = ", mean_degree,"\n"))

			for i in 1:length(empirical_dist)
				if empirical_dist[i] != 0
					write(outfile, string(i, " ", empirical_dist[i], "\n"))
				end
			end

			# Close and save the file.
			close(outfile)
		end

		close(infile)
		return this
	end
end

function graph_properties(G::simpleGraph)
	#= Calculate the mean degree of the graph and
	store the distance for empirical distance =#
	empirical_dist = zeros(Float64, G.num_edge - 1)
	mean_degree = 0.0
	sum = 0

	for i in 1:G.num_vertex
		degree = (G.format_graph == "adjmatrix") ? +(G.Graph[:,i]...) : length(G.Graph[i])
		if degree != 0
			sum += degree
			empirical_dist[degree] += 1
		end
	end

	mean_degree = (sum/G.num_edge)
	empirical_dist /= G.num_vertex
	return mean_degree, empirical_dist
end


text:
module Bucket

using ...SimpleTasks.Types

export download, upload

"""
    download(bucket::BucketService, remote_file::AbstractString,
local_file::Union{AbstractString, IO})

Download a remote file either to a new location `AbstractString` or a stream `IO`.
If the local_file is not specified, a stream to the download is returned
"""
function download(bucket::BucketService, remote_file::AbstractString,
    local_file::Union{AbstractString, IO, Void}=nothing)
    error("download with $bucket is not implemented")
end

"""
    upload(bucket::BucketService, local_file::Union{AbstractString, IO},
remote_file::AbstractString)

Upload a file `AbstractString` or a stream `IO` to the bucket service.
"""
function upload(bucket::BucketService, local_file::Union{AbstractString, IO},
    remote_file::AbstractString)
    error("upload with $bucket is not implemented")
end

"""
    delete(bucket::BucketService, remote_file::AbstractString)

Delete the remote file from bucket
"""
function delete(bucket::BucketService, remote_file::AbstractString)
    error("delete with $bucket is not implemented")
end

end # module Bucket

ext: jl

max_forks_repo_name: seung-lab/SimpleTasks


text:
function distance(x,y)
    return sqrt((x[1] - y[1])^2  + (x[2] - y[2])^2)
end

function read_data()
    try
        count = 0
        g = Dict()
        for line in eachline(STDIN)
            cur = split(strip(line))
            if !in(cur[1], Set(["NAME", "COMMENT", "TYPE", "DIMENSION:",
                                "EDGE_WEIGHT_TYPE", "NODE_COORD_SECTION", "EOF"]))
                g[cur[1]] = [[parse(Int,cur[2]), parse(Int,cur[3])], Any[]]
                count += 1
            end
        end
        for i in keys(g)
            for j in keys(g)
                if i != j
                    dist = distance(g[i][1], g[j][1])
                    push!(g[i][2], (dist, j))
                end
            end
        end
        for i in keys(g)
            sort!(g[i][2], lt=(x,y)->isless(x[1],y[1]))
        end
        return (count, g)
    catch err
        showerror(STDOUT, err, backtrace())
        return (0, Dict())
    end
end

function tsp(g, n)
    cur = "1"
    seen = Set()
    path = []
    count = 0
    dist = 0
    show = false
    for i in 1:n-1
        push!(seen, cur)
        push!(path, cur)
        count += 1
        d, cur = get_nearest_unseen(g[cur][2], seen, count, n)
        dist += d
        show = false
    end
    println(dist)
    println(path)
end

function get_nearest_unseen(ls, seen, num_seen, n)
    for entry in ls
        if  num_seen == n - 1 && entry[2] == "1"
            return entry
        elseif !in(entry[2], seen)
            return  entry
        end
    end
    println("EXITED")
    exit(1)
end

function main()
    n, graph = read_data()
    tsp(graph, n)
end

main()


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module AffineTransforms

import Base: *, \, eltype, ndims, scale, show
using Interpolations, Requires, Base.Cartesian

if !isdefined(:AbstractVecOrMat)
    typealias AbstractVecOrMat Union(AbstractVector, AbstractMatrix)
end

export
    # types
    AffineTransform,
    TransformedArray,
    # functions
    rotation2,
    rotation3,
    rotationmatrix,
    rotationparameters,
    tformeye,
    tformfwd,
    tformfwd!,
    tforminv,
    tforminv!,
    tformfwd3rd,
    tformrigid,
    tformrotate,
    tformscale,
    tformtranslate,
    transform,
    transform!

include("transform.jl")
include("tformedarrays.jl")

end


text:
module Planets

include("radius.jl")
export core_radius

include("accretion.jl")
export accretion_rate

include("habitable_zone.jl")
export habitable_zone

include("stellar_evolution.jl")
export stellar_evolution

end # module


text:
"""
    abstract type SubGrid{T,N} <: AbstractGrid{T,N} end

A subgrid of an underlying grid.
"""
abstract type SubGrid{T,N} <: SimpleLazyGrid{T,N} end

supergrid(g::SubGrid) = g.supergrid

include("indexed.jl")
include("masked.jl")
include("product.jl")
include("boundary.jl")

subgrid(grid::AbstractGrid, domain::Domain) = MaskedGrid(grid, domain)

function subgrid(grid::AbstractEquispacedGrid, domain::AbstractInterval)
    a = infimum(domain)
    b = supremum(domain)
    h = step(grid)
    idx_a = convert(Int, ceil( (a-grid[1])/step(grid))+1 )
    idx_b = convert(Int, floor( (b-grid[1])/step(grid))+1 )
    idx_a = max(idx_a, 1)
    idx_b = min(idx_b, length(grid))
    IndexSubGrid(grid, idx_a:idx_b, domain)
end

function subgrid(grid::ScatteredGrid, domain::Domain)
    mask = in.(grid, Ref(domain))
    points = grid.points[mask]
    ScatteredGrid(points, domain)
end

function subgrid(grid::ProductGrid, domain::ProductDomain)
    if ncomponents(grid) == ncomponents(domain)
        productgrid(map(subgrid, components(grid), components(domain))...)
    else
        MaskedGrid(grid, domain)
    end
end


text:
########################################################################
# Facial projection
########################################################################

"""
    face_project(M, F::AbstractFace, b::Vector)

Obtain the projection of the vector `b` onto the image of the
face `F` under the linear map `M`, i.e., ``MF``. In particular,
`c` solves the least-squares problem

    minimize_c  ½‖(MF)c-b‖²  subj to  Fc ∈ cone(F)

Each face `F` is aware of the required constraint of `c`.
"""
function face_project!(M::AbstractLinearOp, F::AbstractFace, b::Vector{Float64}, c::Vector{Float64}, r::Vector{Float64})
    # create linear map
    MF = M*F
    # solve the projection by lsmr
    lsmr!(c, MF, b)
    # compute residual
    r .= b - MF*c
    return nothing
end


text:
@testset "Intersection methods: TemplateHullIntersection" begin

    # default constructor without types
    @test TemplateHullIntersection() == TemplateHullIntersection{Float64, Vector{Float64}, Missing, Val{false}}(missing, Val{false}())

    # default constructor passing normal vector types
    @test TemplateHullIntersection{Float64, Vector{Float64}}() == TemplateHullIntersection{Float64, Vector{Float64}, Missing, Val{false}}(missing, Val{false}())
    @test TemplateHullIntersection{Float64, SVector{2, Float64}}() == TemplateHullIntersection{Float64, SArray{Tuple{2}, Float64, 1, 2}, Missing, Val{false}}(missing, Val{false}())

    # constructor with a given template
    td = TemplateHullIntersection(BoxDirections(2))
    @test setrep(td) == Union{HPolytope{Float64,LazySets.Arrays.SingleEntryVector{Float64}}, HPolytope{Float64,Vector{Float64}}}
end


text:
export
ruleVBLogNormalOut,
ruleSPLogNormalOutNPP

ruleSPLogNormalOutNPP(  msg_out::Nothing,
                        msg_m::Message{PointMass, Univariate},
                        msg_s::Message{PointMass, Univariate}) =
    Message(Univariate, LogNormal, m=deepcopy(msg_m.dist.params[:m]), s=deepcopy(msg_s.dist.params[:m]))

ruleVBLogNormalOut( dist_out::Any,
                    dist_m::ProbabilityDistribution{Univariate},
                    dist_s::ProbabilityDistribution{Univariate}) =
    Message(Univariate, LogNormal, m=unsafeMean(dist_m), s=unsafeMean(dist_s))


text:
#include("utils.jl")
#include("setdata.jl")

using Test
#==
@testset "dont used compare literals" begin
 @test_skip gelit([],Meta.parse("-P(x)"),[],Meta.parse("-P(x)")) == true
 @test_skip gelit([:x],Meta.parse("-P(x)"),[:x],Meta.parse("-P(x)")) == true
 @test_skip gelit([:x],Meta.parse("-P(x)"),[:y],Meta.parse("-P(y)")) == true
 @test_skip gelit([:x],Meta.parse("-P(x)"),[],Meta.parse("-P(x)")) == true
 @test_skip gelit([],Meta.parse("-P(x)"),[:x],Meta.parse("-P(x)")) == false

 @test_skip gelit([:x,:y],Meta.parse("-P(x,y)"),[:z,:w],Meta.parse("-P(w,z)")) == true
 @test_skip gelit( [:z,:w],Meta.parse("-P(w,f(z))"), [:x,:y],Meta.parse("-P(x,y)")) == false
 @test_skip gelit([:x,:y],Meta.parse("-P(x,y)"),[:z,:w],Meta.parse("-P(w,f(z))")) == true
end

@testset "dont used sort clause" begin
 @test_skip sortcls([:y, :x], Meta.parse("[+P(a),-Q(x,y),-R(a)]").args) == ([:y,:x],Meta.parse("[+P(a),-Q(x,y),-R(a)]").args)
 @test_skip sortcls([:y, :x], Meta.parse("[-R(b),+P(a),-Q(x,y),-R(a)]").args) == ([:y,:x],Meta.parse("[+P(a),-Q(x,y),-R(a),-R(b)]").args)
 @test_skip sortcls([:y, :x], Meta.parse("[-R(b),+P(a),-Q(x,y),+R(a)]").args) == ([:y,:x],Meta.parse("[+P(a), +R(a), -Q(x,y), -R(b)]").args)
end
==#

@testset  "iscap" begin
 @test iscap("abc") == false
 @test iscap("aBC") == false
 @test iscap("Abc") == true
 @test iscap("ABC") == true
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
export latentboxplot
export annotatedheatmap
export plothistory
export plotreconstruction

@userplot LatentBoxplot
@recipe function f(b::LatentBoxplot)
    if length(b.args) != 2 ||
            !(typeof(b.args[1]) <: AbstractMatrix) ||
            !(typeof(b.args[2]) <: AbstractVector)
        error("latentboxplot needs a matrix Z and a vector of labels")
    end

    Z = b.args[1]
    labels = b.args[2]

    @series begin
        seriestype := :boxplot
        legend := false
        bar_width := 0.5
        whisker_width := 0.25

        repeat(labels, outer=size(Z,2)), reshape(Z,:)
    end
end

@userplot AnnotatedHeatmap
@recipe function f(h::AnnotatedHeatmap; annotationargs=(:white,))
    y = h.args[1]              #Get the input arguments, stored in h.args
                               # - in this case there's only one
    typeof(y) <: AbstractMatrix || error("Pass a Matrix as the arg to heatmap")

    grid := false                      # turn off the background grid
    xaxis := false
    yaxis := false

    @series begin                      # the main series, showing the heatmap
        seriestype := :heatmap
        y
    end

    rows, cols = size(y)

    #horizontal lines
    for i in 0:cols         # each line is added as its own series, for clearer code
        @series begin
            seriestype := :path
            primary := false          # to avoid showing the lines in a legend
            linewidth := 2
            linecolor --> :white
            [i, i] .+ 0.5, [0, rows] .+ 0.5  # x and y values of lines
        end
    end

    for i in 0:rows
        @series begin
            seriestype := :path
            primary := false
            linewidth := 2
            linecolor --> :white
            [0, cols] .+ 0.5, [i,i] .+ 0.5
        end
    end

    @series begin
        seriestype := :scatter
        # make the points transparent - setting marker (or seriestype)
        # to :none doesn't currently work right
        markerstrokecolor := RGBA(0,0,0,0.)
        seriescolor := RGBA(0,0,0,0.)
        series_annotations := text.(round.(reshape(y,:), digits=2), annotationargs...)
        primary := false
        repeat(1:cols, inner = rows), repeat(1:rows, outer = cols)
    end
end

@userplot PlotHistory
@recipe function f(p::PlotHistory; exclude_keys=[:xrec])
    if !(typeof(p.args[1]) <: MVHistory)
        error("Can only plot MVHistories")
    end

    history = p.args[1]
    plot_keys = collect(keys(history))
    plot_keys = filter!(x -> !(x in exclude_keys), plot_keys)

    N = length(plot_keys)
    layout := grid(N, 1, heights=ones(N)./N)
    size := (500,1000)

    for (nr, key) in enumerate(plot_keys)
        x, y = get(history, key)
        if ndims(y[1]) == 2
            y = [z[:,1] for z in y]
        end
        y = hcat(y...)'

        @series begin
            seriestype := :path
            subplot := nr
            ylabel := key
            legend := false
            linewidth := 2
            x, y
        end
    end
end

@userplot PlotReconstruction
@recipe function f(p::PlotReconstruction)
    if (length(p.args) != 2) ||
            !(typeof(p.args[1]) <: AbstractMatrix) ||
            !(typeof(p.args[2]) <: AbstractMatrix) ||
            (size(p.args[1]) != size(p.args[2]))
        error("Need two matrices of same size as arguments")
    end

    X = p.args[1]
    Xrec = p.args[2]
    @series begin
        linewidth := 2
        legend := false
        Xrec
    end

    @series begin
        color := "gray"
        linewidth := 2
        legend := false
        X
    end
end


text:
# Time series analysis of 2020 worldwide coronavirus outbreak
#
# data based on John Hopkins CSSE repository
# https://github.com/CSSEGISandData/COVID-19

module C

using  Dates, DSP, PyPlot, DifferentialEquations, Polynomials

include("read_data.jl")

function smooth_log(x)
# make log-transform and smooth data
    lx = log10.(x)
    b = (1/2)*ones(2)     # flat (rectangular) window -> moving average
    lx = filtfilt(b,lx)
    10 .^ lx, lx
end

function local_fit(x,y,di)
# local slope of graph y(x) from regression in running window
# window width: 2*di+1
  len = length(x)
  slope = zeros(len)
  for i=1:len
    i1 = max(1,i-di)
    i2 = min(len,i+di)
    p = polyfit(x[i1:i2],y[i1:i2],1)
    slope[i] = coeffs(p)[2]
  end
  slope, log(10.) .* slope
end

# doubling time when slope is taken from log10-transformed data
doubling_time(slope) = log(2.)/log(10.)/slope

function nlys()

   today = " March 19"

   t_de,I_de,R_de,D_de = covid_ts("Germany",       Date(2020,2,23))
   t_it,I_it,R_it,D_it = covid_ts("Italy",         Date(2020,2,20))
   t_es,I_es,R_es,D_es = covid_ts("Spain",         Date(2020,2,25))
   t_fr,I_fr,R_fr,D_fr = covid_ts("France",        Date(2020,2,25))
   t_uk,I_uk,R_uk,D_uk = covid_ts("United Kingdom",Date(2020,2,27))
   t_at,I_at,R_at,D_at = covid_ts("Austria",       Date(2020,2,29))
   t_nl,I_nl,R_nl,D_nl = covid_ts("Netherlands",   Date(2020,2,23))
   t_ch,I_ch,R_ch,D_ch = covid_ts("Switzerland",   Date(2020,2,25))
   t_cn,I_cn,R_cn,D_cn = covid_ts("China",         Date(2020,1,22))
   t_jp,I_jp,R_jp,D_jp = covid_ts("Japan",         Date(2020,1,24))
   t_kr,I_kr,R_kr,D_kr = covid_ts("Korea, South",  Date(2020,2,17))
   t_us,I_us,R_us,D_us = covid_ts("US",            Date(2020,2,22))
   t_ir,I_ir,R_ir,D_ir = covid_ts("Iran",          Date(2020,2,19))
   t_il,I_il,R_il,D_il = covid_ts("Israel",        Date(2020,2,27))
   t_no,I_no,R_no,D_no = covid_ts("Norway",        Date(2020,2,28))
   t_ca,I_ca,R_ca,D_ca = covid_ts("Canada",        Date(2020,2,27))

  I_de, lI_de = smooth_log(I_de)
  I_it, lI_it = smooth_log(I_it)
  I_fr, lI_fr = smooth_log(I_fr)
  I_es, lI_es = smooth_log(I_es)
  I_jp, lI_jp = smooth_log(I_jp)
  I_cn, lI_cn = smooth_log(I_cn)
  I_kr, lI_kr = smooth_log(I_kr)
  I_ir, lI_ir = smooth_log(I_ir)
  I_us, lI_us = smooth_log(I_us)
  I_ch, lI_ch = smooth_log(I_ch)
  I_il, lI_il = smooth_log(I_il)
  I_at, lI_at = smooth_log(I_at)
  I_uk, lI_uk = smooth_log(I_uk)
  I_nl, lI_nl = smooth_log(I_nl)
  I_no, lI_no = smooth_log(I_no)
  I_ca, lI_ca = smooth_log(I_ca)

  di = 2 # length of running window: 2*di+1
  slope_de,growth_de = local_fit(t_de,lI_de,di)
  slope_it,growth_it = local_fit(t_it,lI_it,di)
  slope_es,growth_es = local_fit(t_es,lI_es,di)
  slope_fr,growth_fr = local_fit(t_fr,lI_fr,di)
  slope_jp,growth_jp = local_fit(t_jp,lI_jp,di)
  slope_cn,growth_cn = local_fit(t_cn,lI_cn,di)
  slope_kr,growth_kr = local_fit(t_kr,lI_kr,di)
  slope_ir,growth_ir = local_fit(t_ir,lI_ir,di)
  slope_us,growth_us = local_fit(t_us,lI_us,di)
  slope_ch,growth_ch = local_fit(t_ch,lI_ch,di)
  slope_at,growth_at = local_fit(t_at,lI_at,di)
  slope_il,growth_il = local_fit(t_il,lI_il,di)
  slope_uk,growth_uk = local_fit(t_uk,lI_uk,di)
  slope_nl,growth_nl = local_fit(t_nl,lI_nl,di)
  slope_no,growth_no = local_fit(t_no,lI_no,di)
  slope_ca,growth_ca = local_fit(t_ca,lI_ca,di)

  p_tot_de1 = polyfit(t_de[4:12],lI_de[4:12],1)  # fit full slope Germany
  coef = coeffs(p_tot_de1)[2]
  println("slope Germany ",coef, " T12 ", doubling_time(coef))
  lI_fit_de1= polyval(p_tot_de1, t_de)

  p_tot_de2 = polyfit(t_de[12:end],lI_de[12:end],1)  # fit full slope Germany
  coef = coeffs(p_tot_de2)[2]
  println("slope Germany ",coef, " T12 ", doubling_time(coef))
  lI_fit_de2= polyval(p_tot_de2, t_de)

  p_tot_it1 = polyfit(t_it[4:9],lI_it[4:9],1)  # fit slope 1st phase Italy
  coef = coeffs(p_tot_it1)[2]
  println("slope Italy Phase 1 ",coef, " T12 ", doubling_time(coef))
  lI_fit_it1 = polyval(p_tot_it1, t_it)

  p_tot_it = polyfit(t_it[10:end],lI_it[10:end],1)  # fit slope 2nd phase Italy
  coef = coeffs(p_tot_it)[2]
  println("slope Italy Phase 2 ",coef, " T12 ", doubling_time(coef))
  lI_fit_it = polyval(p_tot_it, t_it)

  p_tot_es1 = polyfit(t_es[6:18],lI_es[6:18],1)  # fit full slope Spain
  lI_fit_es1= polyval(p_tot_es1, t_es)
  p_tot_es2 = polyfit(t_es[19:24],lI_es[19:24],1)  # fit full slope Spain
  lI_fit_es2= polyval(p_tot_es2, t_es)

  p_tot_fr1 = polyfit(t_fr[3:14],lI_fr[3:14],1)  # fit full slope France
  lI_fit_fr1= polyval(p_tot_fr1, t_fr)
  p_tot_fr2 = polyfit(t_fr[14:end],lI_fr[14:end],1)  # fit full slope France
  lI_fit_fr2= polyval(p_tot_fr2, t_fr)

  p_tot_cn = polyfit(t_cn,lI_cn,1)  # fit full slope China
  lI_fit_cn= polyval(p_tot_cn, t_cn)

  p_tot_kr = polyfit(t_kr,lI_kr,1)  # fit full slope South Korean
  lI_fit_kr= polyval(p_tot_kr, t_kr)

  p_tot_jp = polyfit(t_jp,lI_jp,1)  # fit full slope Japan
  lI_fit_jp= polyval(p_tot_jp, t_jp)

  p_tot_ir = polyfit(t_ir,lI_ir,1)  # fit full slope Iran
  lI_fit_ir= polyval(p_tot_ir, t_ir)

  p_tot_us = polyfit(t_us,lI_us,1)  # fit full slope USA
  lI_fit_us= polyval(p_tot_us, t_us)

  p_tot_ch = polyfit(t_ch[11:22],lI_ch[11:22],1)  # fit full slope Switzerland
  lI_fit_ch= polyval(p_tot_ch, t_ch)

  p_tot_at = polyfit(t_at,lI_at,1)  # fit full slope Austria
  lI_fit_at= polyval(p_tot_at, t_at)

  p_tot_il = polyfit(t_il,lI_il,1)  # fit full slope Israel
  lI_fit_il= polyval(p_tot_il, t_il)

  p_tot_uk = polyfit(t_uk,lI_uk,1)  # fit full slope UK
  lI_fit_uk= polyval(p_tot_uk, t_uk)

  p_tot_no = polyfit(t_no,lI_no,1)  # fit full slope Norway
  lI_fit_no= polyval(p_tot_no, t_no)

  p_tot_nl = polyfit(t_nl[15:23],lI_nl[15:23],1)  # fit full slope Netherlands
  lI_fit_nl= polyval(p_tot_nl, t_nl)

  p_tot_ca = polyfit(t_ca,lI_ca,1)  # fit full slope Norway
  lI_fit_ca= polyval(p_tot_ca, t_ca)

  #majorformatter = matplotlib.dates.DateFormatter("%d.%m.%Y")
  #majorlocator = matplotlib.dates.DayLocator(interval=1)

# ********************************************
  fig = figure("Outbreak Europe",figsize=(12,8))
  clf()
  subplot(3,4,1)         # Germany
      plot(t_de,I_de,"r*")
      ylabel("active cases",fontsize=12)

      title("Germany,"* today)
  ax=subplot(3,4,5)
     # p1 = plot_date(Date(2020,2,17),y,linestyle="-",marker="None",label="test")
      semilogy(t_de,I_de,"r*")
      semilogy(t_de,10 .^ lI_fit_de1,"k")
      semilogy(t_de,10 .^ lI_fit_de2,"k")
      ax.yaxis.grid()
      ylim(5,5e4)
      ylabel("active cases",labelpad=20,fontsize=12)
   ax = subplot(3,4,9)
      plot(t_de,growth_de,"r*-")
      ylim(0,0.6)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      #ax.yaxis.grid(which="major", color="r", linestyle="-", linewidth=2)
      ax.yaxis.grid()
     # annotate("slope 0.107",xy=[20,2],
     # xycoords="axes fraction"
     # )
      #ax.xaxis.set_major_formatter(majorformatter)
      #fig.autofmt_xdate(bottom=0.2,rotation=30,ha="right")
      xlabel("time (days)",fontsize=12)
      ylabel("growth rate",labelpad=20,fontsize=12)

  subplot(3,4,2)          # Italy
      plot(t_it,I_it,"r*")
      #ylabel("infected Italy")
      title("Italy,"*today)
  ax=subplot(3,4,6)
      semilogy(t_it,I_it,"r*")
      semilogy(t_it,10 .^ lI_fit_it1,"k")
      semilogy(t_it,10 .^ lI_fit_it,"k")
      ylim(5,5e4)
      ax.yaxis.grid()
  ax=subplot(3,4,10)
      plot(t_it,growth_it,"r*-")
      ylim(0,0.6)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("slope")

  subplot(3,4,3)              # Spain
      plot(t_es,I_es,"r*")
      title("Spain,"*today)
  ax=subplot(3,4,7)
      semilogy(t_es,I_es,"r*")
      semilogy(t_es,10 .^ lI_fit_es1,"k")
      semilogy(t_es,10 .^ lI_fit_es2,"k")
      ax.yaxis.grid()
      ylim(5,5e4)
   ax= subplot(3,4,11)
      plot(t_es,growth_es,"r*-")
      ylim(0,0.6)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)

  subplot(3,4,4)              # France
      plot(t_fr,I_fr,"r*")
      title("France,"*today)
  ax=subplot(3,4,8)
      semilogy(t_fr,I_fr,"r*")
      semilogy(t_fr,10 .^ lI_fit_fr1,"k")
      semilogy(t_fr,10 .^ lI_fit_fr2,"k")
      ax.yaxis.grid()
      ylim(5,5e4)
  ax = subplot(3,4,12)
      plot(t_fr,growth_fr,"r*-")
      ylim(0,0.6)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")

  PyPlot.tight_layout()
  savefig("../figures/outbreak_europe.png")

# ********************************************
  fig=figure("Outbreak Asia",figsize=(12,8))
  clf()
  subplot(3,4,1)         # China
      plot(t_cn,I_cn,"r*")
      ylabel("active cases",fontsize=12)
      title("China,"*today)
  ax=subplot(3,4,5)
      semilogy(t_cn,I_cn,"r*")
      #plot(t_cn,lI_fit_cn,"k")
      ylim(1e2,1e5)
      ax.yaxis.grid()
      ylabel("active cases", labelpad=20,fontsize=12)
  ax = subplot(3,4,9)
      plot(t_cn,growth_cn,"r*-")
      ylim(-0.1,0.5)
      #yticks([-0.05,0,0.05,0.1,0.15,0.2])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      ylabel("growth",fontsize=12,labelpad=10)

  subplot(3,4,2)          # Japan
      plot(t_jp,I_jp,"r*")
      #ylabel("infected Italy")
      title("Japan,"*today)
  ax=subplot(3,4,6)
      semilogy(t_jp,I_jp,"r*")
      semilogy(t_jp,10 .^ lI_fit_jp,"k")
      ax.yaxis.grid()
      ylim(5,2e4)
      #ylabel("log10 infected Italy")
  ax=subplot(3,4,10)
      plot(t_jp,growth_jp,"r*-")
      ylim(-0.1,0.5)
      #yticks([-0.05,0,0.05,0.1,0.15,0.2])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")

  subplot(3,4,3)              # South Korean
      plot(t_kr,I_kr,"r*")
      #ylabel("infected Spain")
      title("South Korea,"*today)
  ax=subplot(3,4,7)
      semilogy(t_kr,I_kr,"r*")
      ax.yaxis.grid()
      ylim(5,2e4)
      #plot(t_kr,lI_fit_kr,"k")
      #ylabel("log10 infected Spain")
  ax= subplot(3,4,11)
      plot(t_kr,growth_kr,"r*-")
      ylim(-0.1,0.5)
      #yticks([-0.05,0,0.05,0.1,0.15,0.2])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")

  subplot(3,4,4)              # Iran
      plot(t_ir,I_ir,"r*")
      #ylabel("infected")
      title("Iran,"*today)
  ax=subplot(3,4,8)
      semilogy(t_ir,I_ir,"r*")
      semilogy(t_ir,10 .^ lI_fit_ir,"k")
      ax.yaxis.grid()
      ylim(5,2e4)
      #ylabel("log10 infected")
  ax = subplot(3,4,12)
      plot(t_ir,growth_ir,"r*-")
      ylim(-0.1,0.5)
      #yticks([-0.05,0,0.05,0.1,0.15,0.2])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")
  PyPlot.tight_layout()
  savefig("../figures/outbreak_Asia.png")

# ********************************************
  fig=figure("Outbreak World",figsize=(12,8))
  clf()
  subplot(3,4,1)         # USA
      plot(t_us,I_us,"r*")
      ylabel("active cases",fontsize=12)
      title("USA,"*today)
  ax=subplot(3,4,5)
      semilogy(t_us,I_us,"r*")
      semilogy(t_us,10 .^ lI_fit_us,"k")
      ax.yaxis.grid()
      ylim(5,2e4)
      ylabel("active cases",fontsize=12,labelpad=12)
  ax = subplot(3,4,9)
      plot(t_us,growth_us,"r*-")
      ax.yaxis.grid()
      ylim(0,0.6)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      xlabel("time (days)",fontsize=12)
      ylabel("growth",fontsize=12,labelpad=12)

  subplot(3,4,2)              # Canada
      plot(t_ca,I_ca,"r*")
      #ylabel("infected")
      title("Canada,"* today)
  ax=subplot(3,4,6)
      semilogy(t_ca,I_ca,"r*")
      semilogy(t_ca,10 .^ lI_fit_ca,"k")
      ax.yaxis.grid()
      ylim(5,2e4)
      #ylabel("log10 infected")
  ax = subplot(3,4,10)
      plot(t_ca,growth_ca,"r*-")
      ylim(0,0.6)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")

  subplot(3,4,3)              # UK
      plot(t_uk,I_uk,"r*")
      #ylabel("infected Spain")
      title("UK,"*today)
  ax=subplot(3,4,7)
      semilogy(t_uk,I_uk,"r*")
      semilogy(t_uk,10 .^ lI_fit_uk,"k")
      ax.yaxis.grid()
      ylim(5,2e4)
      #ylabel("log10 infected Spain")
  ax= subplot(3,4,11)
      plot(t_uk,growth_uk,"r*-")
      ax.yaxis.grid()
      ylim(0,0.6)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")

  subplot(3,4,4)         # Israel
      plot(t_il,I_il,"r*")
      title("Israel,"*today)
  ax=subplot(3,4,8)
      semilogy(t_il,I_il,"r*")
      semilogy(t_il,10 .^ lI_fit_il,"k")
      ax.yaxis.grid()
      ylim(1,1e3)
  ax = subplot(3,4,12)
      plot(t_il,growth_il,"r*-")
      ax.yaxis.grid()
      ylim(0,0.6)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      xlabel("time (days)",fontsize=12)

  PyPlot.tight_layout()
  savefig("../figures/outbreak_World.png")

# ********************************************
  fig=figure("Outbreak Small Europe",figsize=(12,8))
  clf()
  subplot(3,4,1)         # Netherlands
      plot(t_nl,I_nl,"r*")
      ylabel("active cases",fontsize=12)

      title("Netherlands,"* today)
  ax=subplot(3,4,5)
      semilogy(t_nl,I_nl,"r*")
      semilogy(t_nl,10 .^ lI_fit_nl,"k")
      ax.yaxis.grid()
      ylim(5,1e4)
      ylabel("active cases",labelpad=16,fontsize=12)
   ax = subplot(3,4,9)
      plot(t_nl,growth_nl,"r*-")
      ylim(0,0.6)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      #ax.yaxis.grid(which="major", color="r", linestyle="-", linewidth=2)
      ax.yaxis.grid()
     # annotate("growth 0.107",xy=[20,2],
     # xycoords="axes fraction"
     # )
      #ax.xaxis.set_major_formatter(majorformatter)
      #fig.autofmt_xdate(bottom=0.2,rotation=30,ha="right")
      xlabel("time (days)",fontsize=12)
      ylabel("growth",fontsize=12,labelpad=20)

  subplot(3,4,2)          # Norway
      plot(t_no,I_no,"r*")
      #ylabel("infected Italy")
      title("Norway,"*today)
  ax=subplot(3,4,6)
      semilogy(t_no,I_no,"r*")
      semilogy(t_no,10 .^ lI_fit_no,"k")
      ax.yaxis.grid()
      ylim(5,1e4)
  ax=subplot(3,4,10)
      plot(t_no,growth_no,"r*-")
      ylim(0,0.6)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")

  subplot(3,4,3)              # Switzerland
      plot(t_ch,I_ch,"r*")
      #ylabel("infected Spain")
      title("Switzerland,"* today)
  ax=subplot(3,4,7)
      semilogy(t_ch,I_ch,"r*")
      semilogy(t_ch,10 .^ lI_fit_ch,"k")
      ax.yaxis.grid()
      ylim(5,1e4)
      #ylabel("log10 infected Spain")
  ax= subplot(3,4,11)
      plot(t_ch,growth_ch,"r*-")
      ylim(0,0.6)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")

  subplot(3,4,4)              # Austria
      plot(t_at,I_at,"r*")
      #ylabel("infected")
      title("Austria,"* today)
  ax=subplot(3,4,8)
      semilogy(t_at,I_at,"r*")
      semilogy(t_at,10 .^ lI_fit_at,"k")
      ax.yaxis.grid()
      ylim(5,1e4)
      #ylabel("log10 infected")
  ax = subplot(3,4,12)
      plot(t_at,growth_at,"r*-")
      ylim(0,0.6)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")
  PyPlot.tight_layout()
  savefig("../figures/outbreak_europe1.png")

# ********************************************

   lt_it = log10.(t_it)
   p_tot_it = polyfit(lt_it,lI_it,1)  # fit full slope Germany
   println("loglog It ", coeffs(p_tot_it)[2])
   I_fit_it = 10 .^ polyval(p_tot_it, lt_it)

   lt_de = log10.(t_de)
   p_tot_de = polyfit(lt_de,lI_de,1)  # fit full slope Germany
   println("loglog It ", coeffs(p_tot_de)[2])
   I_fit_de = 10 .^ polyval(p_tot_de, lt_de)

   figure("ts")
   clf()
   subplot(2,1,1)
   loglog(t_it,I_it,"*-")
   loglog(t_it,I_fit_it,"-")

   subplot(2,1,2)
   loglog(t_de,I_de,"*-")
   loglog(t_de,I_fit_de,"-")
   #semilogy(t_es,I_es,"r*-")

end

# **********************************************************
# **********************************************************
# **********************************************************

function nlys_ItDe()
  today = " March 18"
  t_it,I_it,R_it,D_it = covid_ts("Italy",Date(2020,2,21))
  I_it, lI_it = smooth_log(I_it)
  di = 2
  slope_it,growth_it = local_fit(t_it,lI_it,di)

  p_tot_it1 = polyfit(t_it[3:8],lI_it[3:8],1)  # fit slope 1st phase Italy
  coef = coeffs(p_tot_it1)[2]
  println("slope Italy Phase 1 ",coef, " T12 ", doubling_time(coef))
  lI_fit_it1 = polyval(p_tot_it1, t_it)

  p_tot_it2 = polyfit(t_it[9:22],lI_it[9:22],1)  # fit slope 2nd phase Italy
  coef = coeffs(p_tot_it2)[2]
  println("slope Italy Phase 2 ",coef, " T12 ", doubling_time(coef))
  lI_fit_it2 = polyval(p_tot_it2, t_it)

  p_tot_it3 = polyfit(t_it[23:end],lI_it[23:end],1)  # fit slope 2nd phase Italy
  coef = coeffs(p_tot_it3)[2]
  println("slope Italy Phase 3 ",coef, " T12 ", doubling_time(coef))
  lI_fit_it3 = polyval(p_tot_it3, t_it)

  # ********************
  t_de,I_de,R_de,D_de = covid_ts("Germany",  Date(2020,2,24))
  t_de = t_de .+ 3 # bring to the same time axis
  I_de, lI_de = smooth_log(I_de)
  di = 2
  slope_de, growth_de = local_fit(t_de,lI_de,di)

  p_tot_de1 = polyfit(t_de[4:11],lI_de[4:11],1)  # fit full slope Germany
  coef = coeffs(p_tot_de1)[2]
  println("slope Germany ",coef, " T12 ", doubling_time(coef))
  lI_fit_de1= polyval(p_tot_de1, t_de)

  p_tot_de2 = polyfit(t_de[12:end],lI_de[12:end],1)  # fit full slope Germany
  coef = coeffs(p_tot_de2)[2]
  println("slope Germany ",coef, " T12 ", doubling_time(coef))
  lI_fit_de2= polyval(p_tot_de2, t_de)

   figure("Germany_Italy", figsize=(12,6))
   clf()
   subplot(1,2,1)
      semilogy(t_it[1:12],10 .^ lI_fit_it1[1:12],"k")
      semilogy(t_it[5:25],10 .^ lI_fit_it2[5:25],"k")
      semilogy(t_it[19:28],10 .^ lI_fit_it3[19:28],"k")
      semilogy(t_it,I_it,"b*", label="Italy", markersize=8)
      semilogy(t_de[1:15],10 .^ lI_fit_de1[1:15],"k")
      semilogy(t_de[8:end],10 .^ lI_fit_de2[8:end],"k")
      semilogy(t_de,I_de,"r*",label="Germany",markersize=8)
      #plot([10,10],[8,5e4],"grey",linewidth=1)
      #ylim(0.5,5)
      ylim(8,5e4)
      ylabel("active cases",fontsize=14)
      xlabel("time (Feb 21 - March 19)", fontsize=14)
      legend(loc="upper left",fancybox="true")
   ax = subplot(1,2,2)
      plot(t_it,growth_it,"b*-",label="Italy",markersize=8)
      plot(t_de,growth_de,"r*-",label="Germany",markersize=8)
      #plot([10,10],[0,0.6],"grey",linewidth=1)
      #ylim(0,0.22)
      #yticks([0,0.05,0.1,0.15,0.2])
      ylim(0,0.6)
      #yticks([0,0.05,0.1,0.15,0.2])
      ax.yaxis.grid()
      #xlabel("time (days)",fontsize=14)
      xlabel("time (Feb 21 - March 18)", fontsize=14)
      ylabel("growth rate",fontsize=14)
      legend(loc="upper right",fancybox="true")
  PyPlot.tight_layout()
  savefig("../figures/ItalyGermany_March19.png")

end

#country = "Italy"; date = Date(2020,2,20)

#nlys_country(country, date)
nlys_ItDe()

nlys()

end   # module

max_stars_repo_name: AlexRyabov/covid

max_issues_repo_name: AlexRyabov/covid

max_forks_repo_name: AlexRyabov/covid


text:
function solve_sai(U0, t_grid, dt, lambda::Number)

    grid_len = length(t_grid)
    s = t_grid[1]

    t_outer = t_grid .- t_grid'

    K0 = make_volterra_kernel(U0)
    pert = make_volterra_kernel(t_outer)

    #c1 = airyai(s) * exp(-lambda * sqrt(s))
    #c2 = -sqrt(s) * c1

    c1 = airyai(s + lambda)
    c2 = airyaiprime(s + lambda)

    v1 = ones(grid_len)
    v2 = U0[:,1]
    v3 = t_outer[:,1]

    v = (c2 * v1) + (c1 * v2) + (c1 * lambda * v3)

    A = I - dt/2 * K0
    B = -dt/2 * pert

    sai_path = (A + lambda * B) \ v

    trapz_matrix = make_volterra_kernel(ones(size(U0))) * dt/2

    return (
        SAi = trapz_matrix * sai_path,
        DSAi = sai_path
    )
end

function solve_phi(U0, t_grid, dt, lambda_grid; verbose=false)

    grid_len = length(t_grid)

    s = t_grid[1]

    t_outer = t_grid .- t_grid'

    K0 = make_volterra_kernel(U0)
    pert = make_volterra_kernel(t_outer)

    u_deriv = [zeros(grid_len - 1); 1]
    u_full = dt/2 * ([0; ones(grid_len - 1)] + [ones(grid_len - 1); 0]);

    u = [u_deriv u_full]

    v1 = ones(grid_len)
    v2 = U0[:,1]
    v3 = t_outer[:,1]

    A = I - dt/2 * K0
    B = -dt/2 * pert

    BA_inv = B / A
    uA_inv = A' \ u

    N = BA_inv

    Nv_curr = [v1 v2 v3]

    # the 2*3 matrix u^t A^-1 v for u1-u2, v1-v3 for each power
    coefs_v = Array{Complex{Float64}, 3}(undef, (grid_len, 2, 3))

    for pow in 1:grid_len

        coefs_v[pow,:,:] = uA_inv' * Nv_curr * (-1)^(pow + 1)
        Nv_curr = N * Nv_curr

     end

    log_coefs_v = log.(coefs_v)

    phi_deriv_grid = undef_copy(lambda_grid) * 1.0
    phi_full_grid = undef_copy(lambda_grid) * 1.0

    for (i, lambda) in pairs(lambda_grid)

        c1 = airyai(s + lambda)
        c2 = airyaiprime(s + lambda)

        log_lambda = log(Complex(lambda))

        # Make the array [0, log(lambda), 2 log(lambda), ...]
        log_lambda_pows = accumulate(+, [0; log_lambda * ones(grid_len - 1)])

        # Make the array [C_0(lambda) lambda^0, C_1(lambda) lambda^1, ...]
        sums_deriv = sum(exp.(log_lambda_pows .+ log_coefs_v[:,1,:]), dims = 1)
        sums_full = sum(exp.(log_lambda_pows .+ log_coefs_v[:,2,:]), dims = 1)

        if any(isnan.(sums_deriv)) | any(isnan.(sums_full))
            ss = log_lambda_pows .+ log_coefs_v + maxval_log_cum
            display(ss)
            display(cumsum(ss, dims = 1))
            display(sums)
            throw("NAN sums")
        end

        # If lamba grid is purely real, make phi values real
        if eltype(lambda_grid) <: Real
            #@assert maximum(abs.(imag.(sums_deriv))/abs.(sums_deriv)) < 1e-6
            #@assert maximum(abs.(imag.(sums_full))/abs.(sums_deriv)) < 1e-6

            sums_deriv = real.(sums_deriv)
            sums_full = real.(sums_full)
        end

        phi_deriv_grid[i] = (c2 * sums_deriv[1,1] + c1 * sums_deriv[1,2] + c1 * lambda * sums_deriv[1,3]) + c1
        phi_full_grid[i] = (c2 * sums_full[1,1] + c1 * sums_full[1,2] + c1 * lambda * sums_full[1,3]) + c1

        if isnan(phi_deriv_grid[i]) | isnan(phi_full_grid[i])
            println(i)
            throw("NAN phi")
        end
    end

    if verbose
        return (
            SAi = phi_full_grid,
            DSAi = phi_deriv_grid,
            Vcoefs = coefs_v,
            LogVcoefs = log_coefs_v,
            MaxLogs = maxval_logs)
    else
        return (
            SAi = phi_full_grid,
            DSAi = phi_deriv_grid)
    end
end;


text:
import PipeLayout:  dist
using PipeLayout

@testset "build grid topology with 4 neighbors" begin
    vert, hori, degree, width = 3, 4, 4, 5.0
    grid = RegGrid{vert,hori,degree}(width)
    topo = gridtopology(grid)
    nodes = topo.nodes
    arcs = topo.arcs

    narcs = 2*((vert - 1)*hori + vert*(hori - 1))
    @test length(nodes) == vert*hori
    @test length(arcs) == narcs
    for arc in arcs
        tail, head = arc
        @test dist(nodes[tail], nodes[head]) ≈ width
    end

    # bottom left node should be at origin
    nodemat = reshape(nodes, (vert, hori))
    @test nodemat[vert, 1].x ≈ 0.0
    @test nodemat[vert, 1].y ≈ 0.0
end

@testset "build grid topology with 8 neighbors" begin
    vert, hori, degree, width = 3, 4, 8, 5.0
    grid = RegGrid{vert,hori,degree}(width)
    topo = gridtopology(grid)
    nodes = topo.nodes
    arcs = topo.arcs

    narcs = 2*((vert - 1)*hori + vert*(hori - 1) + 2*(vert - 1)*(hori - 1))
    @test length(nodes) == vert*hori
    @test length(arcs) == narcs

    diag = sqrt(2.0) * width
    for arc in arcs
        tail, head = arc
        d = dist(nodes[tail], nodes[head])
        @test  d ≈ width || d ≈ diag
    end

    # bottom left node should be at origin
    nodemat = reshape(nodes, (vert, hori))
    @test nodemat[vert, 1].x ≈ 0.0
    @test nodemat[vert, 1].y ≈ 0.0
end


text:
Base.length(v::ValueTypes) = Int(@with_scope begin
    nv = NapiValue(v)
    if is_typedarray(nv)
        get_array_info(nv, :typedarray).length
    elseif is_arraybuffer(nv)
        get_array_info(nv, :arraybuffer).length
    elseif is_buffer(nv)
        get_array_info(nv, :buffer).length
    elseif is_dataview(nv)
        get_array_info(nv, :dataview).length
    elseif is_array(nv)
        @napi_call napi_get_array_length(nv::NapiValue)::UInt32
    else
        nv.length
    end
end)

function get_array_info(arr::NapiValue, array_type=:typedarray)
    data = Ref{Ptr{Cvoid}}()
    if array_type == :arraybuffer
        byte_length = @napi_call napi_get_arraybuffer_info(arr::NapiValue, data::Ptr{Ptr{Cvoid}})::Csize_t
        return NapiArrayBufferInfo(data[], byte_length)
    elseif array_type == :buffer
        byte_length = @napi_call napi_get_buffer_info(arr::NapiValue, data::Ptr{Ptr{Cvoid}})::Csize_t
        return NapiArrayBufferInfo(data[], byte_length)
    end

    typ = Ref{NapiTypedArrayType}()
    len = Ref{Csize_t}()
    arraybuffer = Ref{NapiValue}()
    byte_offset = Ref{Csize_t}()
    if array_type == :dataview
        @napi_call napi_get_dataview_info(
            arr::NapiValue,
            len::Ptr{Csize_t},
            data::Ptr{Ptr{Cvoid}},
            arraybuffer::Ptr{NapiValue},
            byte_offset::Ptr{Csize_t}
        )
    elseif array_type == :typedarray
        @napi_call napi_get_typedarray_info(
            arr::NapiValue,
            typ::Ptr{NapiTypedArrayType},
            len::Ptr{Csize_t},
            data::Ptr{Ptr{Cvoid}},
            arraybuffer::Ptr{NapiValue},
            byte_offset::Ptr{Csize_t}
        )
    end
    NapiTypedArrayInfo(typ[], len[], data[], arraybuffer[], byte_offset[])
end

function _napi_value(v::AbstractArray)
    n = length(v)
    nv = @napi_call napi_create_array_with_length(n::Csize_t)::NapiValue
    for i in 1:n
        nv[i-1] = v[i]
    end
    nv
end

arraybuffer_finalizer(ptr) = dereference(ptr)
const COPY_ARRAY = Ref(false)
set_copy_array(copy::Bool) = COPY_ARRAY[] = copy

function _napi_value(v::TypedCompatibleArray, typedarray_type; copy_array=COPY_ARRAY[])
    isnothing(typedarray_type) && return _napi_value(v)
    T = eltype(v)
    n = length(v)
    byte_length = sizeof(T) * n
    arraybuffer = if copy_array
        data = Ref{Ptr{Cvoid}}()
        ab = @napi_call napi_create_arraybuffer(byte_length::Csize_t, data::Ptr{Ptr{Cvoid}})::NapiValue
        @ccall memcpy(data[]::Ptr{Cvoid}, v::Ptr{Cvoid}, byte_length::Csize_t)::Ptr{Cvoid}
        ab
    else
        ptr = pointer(reference(v))
        @napi_call create_external_arraybuffer(ptr::Ptr{Cvoid}, byte_length::Csize_t)::NapiValue
    end
    nv = @napi_call napi_create_typedarray(
        typedarray_type::NapiTypedArrayType, n::Csize_t,
        arraybuffer::NapiValue, 0::Csize_t
    )::NapiValue
    nv
end

function napi_value(
    v::TypedCompatibleArray{T};
    copy_array=COPY_ARRAY[], typedarray_type=nothing
) where T
    isnothing(typedarray_type) || return _napi_value(v, typedarray_type; copy_array=copy_array)
    _napi_value(v, if T == Int8
        NapiTypes.napi_int8_array
    elseif T == UInt8
        NapiTypes.napi_uint8_array
    elseif T == Int16
        NapiTypes.napi_int16_array
    elseif T == UInt16
        NapiTypes.napi_uint16_array
    elseif T == Int32
        NapiTypes.napi_int32_array
    elseif T == UInt32
        NapiTypes.napi_uint32_array
    elseif T == Float32
        NapiTypes.napi_float32_array
    elseif T == Float64
        NapiTypes.napi_float64_array
    elseif T == Int64
        NapiTypes.napi_bigint64_array
    elseif T == UInt64
        NapiTypes.napi_biguint64_array
    else
        nothing
    end; copy_array=copy_array)
end
napi_value(v::AbstractArray) = _napi_value(v, nothing)

const UnsafelyWrapped = Dict{Ptr{Cvoid}, NodeObject}()
array_finalizer(arr) = if initialized()
    ptr = Ptr{Cvoid}(pointer(arr))
    if ptr in keys(UnsafelyWrapped)
        ref = UnsafelyWrapped[ptr]
        if dec_ref(ref) == 0
            deleteat!(UnsafelyWrapped, ptr)
        end
    end
end
function unsafe_wrap_array(nv::NapiValue, ::Type{T}, data::Ptr, len) where T
    data = Ptr{Cvoid}(data)
    ref = if data in keys(UnsafelyWrapped)
        ref = UnsafelyWrapped[data]
        inc_ref(ref)
        ref
    else
        NodeObject(nv)
    end
    arr = unsafe_wrap(Array{T}, Ptr{T}(data), len)
    finalizer(array_finalizer, arr)
end
Base.Array(v::ValueTypes) = value(Array, v)
value(::Type{T}, v::NapiValue) where T <: AbstractArray = T(value(Array, v))
value(
    ::Type{Array}, v::NapiValue;
    array_type=nothing
) = if array_type == :typedarray || isnothing(array_type) && is_typedarray(v)
    info = get_array_info(v, :typedarray)
    T = if info.type == NapiTypes.napi_int8_array
        Int8
    elseif info.type == NapiTypes.napi_uint8_array || info.type == NapiTypes.napi_uint8_clamped_array
        UInt8
    elseif info.type == NapiTypes.napi_int16_array
        Int16
    elseif info.type == NapiTypes.napi_uint16_array
        UInt16
    elseif info.type == NapiTypes.napi_int32_array
        Int32
    elseif info.type == NapiTypes.napi_uint32_array
        UInt32
    elseif info.type == NapiTypes.napi_float32_array
        Float32
    elseif info.type == NapiTypes.napi_float64_array
        Float64
    elseif info.type == NapiTypes.napi_bigint64_array
        Int64
    elseif info.type == NapiTypes.napi_biguint64_array
        UInt64
    else
        UInt8 # should not be here
    end
    unsafe_wrap_array(v, T, info.data, info.length)
elseif array_type == :arraybuffer || isnothing(array_type) && is_arraybuffer(v)
    info = get_array_info(v, :arraybuffer)
    unsafe_wrap_array(v, UInt8, info.data, info.length)
elseif array_type == :buffer || isnothing(array_type) && is_buffer(v)
    info = get_array_info(v, :buffer)
    unsafe_wrap_array(v, UInt8, info.data, info.length)
elseif array_type == :dataview || isnothing(array_type) && is_dataview(v)
    info = get_array_info(v, :dataview)
    unsafe_wrap_array(v, UInt8, info.data, info.length)
else
    len = length(v)
    [v[i] for i in 0:len-1]
end

Base.Tuple(v::ValueTypes) = value(Tuple, v)
value(::Type{T}, v::NapiValue) where T <: Tuple = T(value(Tuple, v))
value(::Type{Tuple}, v::NapiValue) = @with_scope begin
    len = length(v)
    Tuple(v[i] for i in 0:len-1)
end


text:
using Documenter, LibSerialPort, LibSerialPort.Lib

makedocs(
    sitename="LibSerialPort.jl",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    )
)

deploydocs(
    repo = "github.com/JuliaIO/LibSerialPort.jl.git",
)


text:
Base.@deprecate contiguous_view ContiguousView
Base.@deprecate strided_view StridedView


text:
type Node{T, D}
  data::T
  lft::Integer
  rgt::Integer
  children::OrderedDict{Ptr, Node{D, D}}
  parent::Union{Void,Node}
  key::Ptr
end
type NodeDuplicationException <: Exception
  parent::Node
  child::Node
end

type NodeIsRootException <: Exception
  node::Node
end

type NodeIterator{T}
  node::Node
  eltype::Type{T}
  task::Task
end

NodeIterator{T, D}(node::Node{T, D}, skipself=false) =
  NodeIterator(node,
    ((skipself || T == D) ? Node{D, D} : Union{Node{T, D}, Node{D, D}}),
    @task iternode(produce, node, skipself)
  )

Base.start(it::NodeIterator) = start(it.task)
Base.done(it::NodeIterator, state) = done(it.task, state)
Base.next(it::NodeIterator, state) = next(it.task, state)


text:
# Copyright (c) 2020 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export posterior

"""
`posterior(Nt,a,A,y)`

Computes the state posterior probabilities (``\\gamma``) for an HMM model with initial probability `a`, transition probability matrix `A`  and observation likelihoods `y`.
"""
function posterior(Nt,a::AbstractArray{T},A,y) where {T<:AbstractFloat}
  alpha, c = forward(Nt,a,A,y)
  beta = backward(Nt,A,c,y)
  alpha .*= beta
  return alpha
end

"""
`posterior(Nt,t2tr,A,y)`

Computes the state posterior probabilities (``\\gamma``) for an HMM model with constrained path `t2tr` (see [`forward`](@ref)), transition probability matrix `A`  and observation likelihoods `y`.
"""
function posterior(Nt,t2tr,A,y)
  alpha, c = forward(Nt,t2tr,A,y)
  beta = backward(Nt,A,c,t2tr,y)
  alpha .*= beta
  return alpha
end

export logposterior

"""
`logposterior(Nt,a,A,y)`

Computes the state posterior log-probabilities (``\\hat{\\gamma}``) for an HMM model with initial probability `a`, transition probability matrix `A`  and observation likelihoods `y`.
"""
function logposterior(Nt,a::AbstractArray{T},A,y) where {T<:AbstractFloat}
  alpha, logML = logforward(Nt,a,A,y)
  beta = logbackward(Nt,A,y)
  alpha .+= beta .- logML .* Nt
  return alpha
end

"""
`logposterior(Nt,t2tr,A,y)`

Computes the state posterior log-probabilities (``\\hat{\\gamma}``) for an HMM model with constrained path `t2tr` (see [`forward`](@ref)), transition probability matrix `A`  and observation likelihoods `y`.
"""
function logposterior(Nt,t2tr,A,y)
  alpha, logML = logforward(Nt,t2tr,A,y)
  beta = logbackward(Nt,A,t2tr,y)
  alpha .+= beta .- logML .* Nt
  return alpha
end


text:
using GoogleCloud
using Test

include("api.jl")
#include("storage.jl")


text:
using apo
using Test

@testset "apo.jl" begin
    # Write your tests here.
    @test true == true
end


text:
#!/usr/bin/env julia

function char_to_digit(c::Char)
    if isdigit(c)
        return Int(c) - Int('0')
    else
        throw(ArgumentError(@sprintf("%c is not a digit", c)))
    end
end

function main()
    s = string(BigInt(2) ^ BigInt(1000))
    digit_sum::Int64 = 0
    for i in 1:length(s)
        digit_sum += char_to_digit(s[i])
    end
    @printf("The sum of teh digits of 2 ^ 1000 is %d\n", digit_sum)
end

main()


text:
# Leaky Integrate and Fire Demo -- see section 26.5.4 of the Bayesian Reasoning and Machine Learning Textbook
# David Barber UCL 2015
# Uses the dbAutoDiff Julia package

PlotResults=true
useproc("CPU") # CPU is faster than GPU for this demo

T=100 # number of timesteps
a=Array(ADnode,T)
p=Array(ADnode,T)
v=Array(ADnode,T)
loss=Array(ADnode,T)

ThetaRest=-5.0
ThetaFired=-15.0
alpha=0.95

StartCode()  # This is where you start the code for your problem
W=ADvariable() # parameter that we want to learn
a[1]=ADnode() # A variable that we will not take the derivative of
p[1]=sigmoid(a[1])
v[1]=ADnode()
for t=2:T
    v[t]=ADnode()
    a[t]=(alpha*a[t-1]+W*v[t]+ThetaRest*(1-alpha)).*(1-v[t-1])+v[t-1]*ThetaFired
    p[t]=sigmoid(a[t]) # probability that neurons fire
    loss[t]=BinaryKullbackLeiblerLoss(v[t],p[t])
end
totalloss=mean(loss[2:T])
net=EndCode() # defines the graph

# instantiate parameter nodes and inputs:
N=20 # number of neurons at each time step
net.value[W]=0.001*randn(N,N)

# training data:
net.value[a[1]]=zeros(N,1)
for t=1:T
    net.value[v[t]]=0.99*(rand(N,1).>0.85)+0.001 # bound slighty away from 0/1 so that no NAN occur
end

net=compile(net) # compile and preallocate memory

@gpu CUDArt.init([0])
@gpu net=convert(net,"GPU")
#gradcheck(net) # only use a small network to check the gradient, otherwise this will take a long time

# Training:
println("Training: using $(net.gpu==true? "GPU" : "CPU") ")
nupdates=1000
parstoupdate=Parameters(net)
error=Array(Float64,0)
    velo=NesterovInit(net) # Nesterov velocity
    LearningRate=5.5
    for i=1:nupdates
    ADforward!(net)
    ADbackward!(net)
    push!(error,extract(net.value[net.FunctionNode])[1])
    printover("iteration $i: training loss = $(error[i])")
    for par in parstoupdate
        #GradientDescentUpdate!(net.value[par],net.gradient[par],LearningRate)
        NesterovGradientDescentUpdate!(net.value[par],net.gradient[par],velo[par],LearningRate,i/500)
    end
end

trainvalue=zeros(N,T)
testprob=zeros(N,T)
actpot=zeros(N,T)

@gpu net=convert(net,"CPU")

for t=1:T
    trainvalue[:,t]=net.value[v[t]]
    testprob[:,t]=net.value[p[t]]
    actpot[:,t]=net.value[a[t]]
end

if PlotResults
    figure(); imshow(trainvalue, interpolation="none"); title("training firings v_i(t)")
    figure(); imshow(testprob, interpolation="none"); title("test firing probability p_i(t)")
    figure(); imshow(actpot, interpolation="none"); title("test membrane potential a_i(t)")
    figure(); plot(error); title("training loss against iteration")
end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
# constructors definitions
# export non_caching_grad

for (f,df,ddf) in f_nargs_1
    @eval begin
        $f(e::E) where E <: Expression = Expression1($f,e)
        @inline (e::Expression1{typeof($f),E})(x,p=nothing) where {E} = setrefval(e,$f(e.e1(x,p)))

        @inline function (d::Gradient1{typeof($f),F1})(y,x,p=nothing,d0=1.) where {F1}
            setrefval(d,$df(frefval1(d)))
            d.d1(y,x,p,d0*refval(d))
            return
        end

        @inline function (h::Hessian11{typeof($f),H1,H11})(z,x,p=nothing,h0=1) where {H1,H11}
            h.h1(z,x,p,h0*refval(h))
            h.h11(z,x,p,h0*$ddf(frefval(h)))
            return
        end

        @inline non_caching_eval(e::Expression1{typeof($f),E},x,p=nothing) where {E} = $f(non_caching_eval(e.e1,x,p))
        @inline function non_caching_eval(d::Gradient1{typeof($f),D1},y,x,p=nothing,d0=1.)  where D1
            non_caching_eval(d.d1,y,x,p,d0*$df(frefval1(d)))
            return
        end
        # @inline non_caching_grad(e::Expression1{typeof($f),F},y,x,p=nothing,d0=1.) where F =
        #     non_caching_grad(e.e1,y,x,p,d0*$df(e.e1(x)))
    end
end

for (f,df1,df2,ddf11,ddf12,ddf22) in f_nargs_2
    @eval begin
        $f(e1::E1,e2::E2) where {E1,E2} = Expression2($f,e1,e2)

        @inline (e::Expression2{typeof($f),F1,F2})(x,p=nothing) where {F1,F2} = setrefval(e,$f(e.e1(x,p),e.e2(x,p)))
        @inline (e::Expression2{typeof($f),F1,F2})(x,p=nothing) where {F1<:Real,F2} = setrefval(e,$f(e.e1,e.e2(x,p)))
        @inline (e::Expression2{typeof($f),F1,F2})(x,p=nothing) where {F1,F2<:Real} = setrefval(e,$f(e.e1(x,p),e.e2))

        @inline function (d::Gradient2F1{typeof($f),F1,R})(y,x,p=nothing,d0=1.) where {F1,R}
            setrefval(d,$df2(d.a,frefval1(d)))
            d.d1(y,x,p,d0*refval(d))
            return
        end
        @inline function (d::Gradient2F2{typeof($f),F1,R})(y,x,p=nothing,d0=1.) where {F1,R}
            setrefval(d,$df1(frefval1(d),d.a))
            d.d1(y,x,p,d0*refval(d))
            return
        end
        @inline function (d::Gradient2{typeof($f),D1,D2})(y,x,p=nothing,d0 = 1) where {D1,D2}
            setrefval1(d,$df1(frefval1(d),frefval2(d)))
            setrefval2(d,$df2(frefval1(d),frefval2(d)))
            d.d1(y,x,p,d0*refval1(d))
            d.d2(y,x,p,d0*refval2(d))
            return
        end
        @inline function (h::Hessian11F1{typeof($f),H1,H11,R})(z,x,p=nothing,h0=1) where {H1,H11,R}
            h.h1(z,x,p,h0*refval(h))
            h.h11(z,x,p,h0*$ddf22(h.a,frefval(h)))
            return
        end
        @inline function (h::Hessian11F2{typeof($f),H1,H11,R})(z,x,p=nothing,h0=1) where {H1,H11,R}
            h.h1(z,x,p,h0*refval(h))
            h.h11(z,x,p,h0*$ddf11(frefval(h),h.a))
            return
        end
        @inline function (h::Hessian22{typeof($f),H1,H2,H11,H12,H22})(z,x,p=nothing,h0=1) where {H1,H2,H11,H12,H22}
            ddf12 = $ddf12(frefval1(h),frefval2(h))
            h.h1(z,x,p,h0*refval1(h))
            h.h2(z,x,p,h0*refval2(h))
            h.h11(z,x,p,h0*$ddf11(frefval1(h),frefval2(h)))
            h.h12(z,x,p,h0*ddf12)
            h.h21(z,x,p,h0*ddf12)
            h.h22(z,x,p,h0*$ddf22(frefval1(h),frefval2(h)))
            return
        end

        @inline non_caching_eval(e::Expression2{typeof($f),F1,F2},x,p=nothing) where {F1,F2} = $f(non_caching_eval(e.e1,x,p),non_caching_eval(e.e2,x,p))
        @inline non_caching_eval(e::Expression2{typeof($f),F1,F2},x,p=nothing) where {F1<:Real,F2} = $f(e.e1,non_caching_eval(e.e2,x,p))
        @inline non_caching_eval(e::Expression2{typeof($f),F1,F2},x,p=nothing) where {F1,F2<:Real} = $f(non_caching_eval(e.e1,x,p),e.e2)

        @inline function non_caching_eval(d::Gradient2F1{typeof($f),D1,R},y,x,p=nothing,d0=1.) where {D1,R}
            non_caching_eval(d.d1,y,x,p,d0*$df2(d.a,frefval1(d)))
            return
        end
        @inline function non_caching_eval(d::Gradient2F2{typeof($f),D1,R},y,x,p=nothing,d0=1.) where {D1,R}
            non_caching_eval(d.d1,y,x,p,d0*$df1(frefval1(d),d.a))
            return
        end
        @inline function non_caching_eval(d::Gradient2{typeof($f),D1,D2},y,x,p=nothing,d0 = 1) where {D1,D2}
            non_caching_eval(d.d1,y,x,p,d0*$df1(frefval1(d),frefval2(d)))
            non_caching_eval(d.d2,y,x,p,d0*$df2(frefval1(d),frefval2(d)))
            return
        end

        # @inline function non_caching_grad(e::Expression2{typeof($f),F1,F2},y,x,p=nothing,d0=1.) where {F1,F2}
        #     f1 = e.e1(x,p)
        #     f2 = e.e2(x,p)
        #     non_caching_grad(e.e1,y,x,p,d0*$df1(f1,f2))
        #     non_caching_grad(e.e2,y,x,p,d0*$df2(f1,f2))
        #     return
        # end
        # @inline non_caching_grad(e::Expression2{typeof($f),F1,F2},y,x,p=nothing,d0=1.) where {F1<:Real,F2} = non_caching_grad(e.e2,y,x,p,d0*$df2(e.e1,e.e2(x)))
        # @inline non_caching_grad(e::Expression2{typeof($f),F1,F2},y,x,p=nothing,d0=1.) where {F1,F2<:Real} = non_caching_grad(e.e1,y,x,p,d0*$df1(e.e1(x),e.e2))
    end
end

add_sum(e1::E,e2) where E <: Expression = add_sum(ExpressionSum([e1]),e2)
add_sum(e1::ExpressionSum{E,I},e2) where {E,I} = _add_sum(e1,e2) ? e1 : ExpressionSum(e1,[e2])
function _add_sum(e1::ExpressionSum{E,I},e2) where {E,I}
    if e2 isa eltype(e1.es)
        push!(e1.es,e2)
        return true
    else
        return _add_sum(inner(e1),e2)
    end
end
function _add_sum(e1::ExpressionSum{E,Nothing},e2) where E
    if e2 isa eltype(e1.es)
        push!(e1.es,e2)
        return true
    else
        return false
    end
end


text:
using LazySets, LazySets.Approximations, Test, LinearAlgebra, SparseArrays

import IntervalArithmetic
const IA = IntervalArithmetic
using IntervalArithmetic: IntervalBox

# ========================
# Optional dependencies
# ========================
import Distributions, Expokit, IntervalMatrices, Optim, TaylorModels
using IntervalMatrices: ±, IntervalMatrix
using TaylorModels: set_variables, TaylorModelN

# ==============================
# Non-exported helper functions
# ==============================
using LazySets: ispermutation
using LazySets.Arrays: isinvertible, inner,
                       is_cyclic_permutation, SingleEntryVector

global test_suite_basic = true
global test_suite_doctests = true
global test_suite_polyhedra = true
global test_suite_plotting = true

if (length(ARGS) == 0) || (ARGS[1] == "--default")
    # default test suite including doctests
elseif ARGS[1] == "--basic"
    # basic test suite
    test_suite_doctests = false
    test_suite_polyhedra = false
    test_suite_plotting = false
elseif ARGS[1] == "--polyhedra"
    # Polyhedra.jl test suite
    test_suite_doctests = false
    test_suite_polyhedra = true
    test_suite_plotting = false
elseif ARGS[1] == "--plot"
    # plotting test suite
    test_suite_doctests = false
    test_suite_polyhedra = false
    test_suite_plotting = true
elseif ARGS[1] == "--all"
    # complete test suite
    test_suite_doctests = true
    test_suite_polyhedra = true
    test_suite_plotting = true
else
    error("unknown parameter 1: $(ARGS[1])")
end

if test_suite_polyhedra || test_suite_plotting
    import Polyhedra
    using CDDLib # for tests that require CDDLib specific backend=...

    # fix namespace conflicts with Polyhedra
    using LazySets: dim, HalfSpace, Interval, Line, translate
end

if test_suite_basic

    # =========================
    # Testing utility functions
    # =========================
    @time @testset "LazySets.Util" begin include("unit_util.jl") end
    @time @testset "LazySets.Comparisons" begin include("unit_comparisons.jl") end
    @time @testset "LazySets.Interior" begin include("unit_interior.jl") end

    # =======================================
    # Testing interfaces to external packages
    # =======================================
    @time @testset "LazySets.CDDLib" begin include("unit_CDDLib.jl") end

    # =======================================
    # Testing types that inherit from LazySet
    # =======================================
    @time @testset "LazySets.Singleton" begin include("unit_Singleton.jl") end
    @time @testset "LazySets.Ball1" begin include("unit_Ball1.jl") end
    @time @testset "LazySets.Ball2" begin include("unit_Ball2.jl") end
    @time @testset "LazySets.Ballp" begin include("unit_Ballp.jl") end
    @time @testset "LazySets.BallInf" begin include("unit_BallInf.jl") end
    @time @testset "LazySets.Hyperrectangle" begin include("unit_Hyperrectangle.jl") end
    @time @testset "LazySets.Polygon" begin include("unit_Polygon.jl") end
    @time @testset "LazySets.Polytope" begin include("unit_Polytope.jl") end
    @time @testset "LazySets.Polyhedron" begin include("unit_Polyhedron.jl") end
    @time @testset "LazySets.Zonotope" begin include("unit_Zonotope.jl") end
    @time @testset "LazySets.ZeroSet" begin include("unit_ZeroSet.jl") end
    @time @testset "LazySets.EmptySet" begin include("unit_EmptySet.jl") end
    @time @testset "LazySets.Ellipsoid" begin include("unit_Ellipsoid.jl") end
    @time @testset "LazySets.Hyperplane" begin include("unit_Hyperplane.jl") end
    @time @testset "LazySets.HalfSpace" begin include("unit_HalfSpace.jl") end
    @time @testset "LazySets.Interval" begin include("unit_Interval.jl") end
    @time @testset "LazySets.LineSegment" begin include("unit_LineSegment.jl") end
    @time @testset "LazySets.Line" begin include("unit_Line.jl") end
    @time @testset "LazySets.Universe" begin include("unit_Universe.jl") end

    # =========================================
    # Testing types representing set operations
    # =========================================
    @time @testset "LazySets.Bloating" begin include("unit_Bloating.jl") end
    @time @testset "LazySets.Intersection" begin include("unit_Intersection.jl") end
    @time @testset "LazySets.ConvexHull" begin include("unit_ConvexHull.jl") end
    @time @testset "LazySets.ExponentialMap" begin include("unit_ExponentialMap.jl") end
    @time @testset "LazySets.LinearMap" begin include("unit_LinearMap.jl") end
    @time @testset "LazySets.MinkowskiSum" begin include("unit_MinkowskiSum.jl") end
    @time @testset "LazySets.CartesianProduct" begin include("unit_CartesianProduct.jl") end
    @time @testset "LazySets.ResetMap" begin include("unit_ResetMap.jl") end
    @time @testset "LazySets.SymmetricIntervalHull" begin include("unit_SymmetricIntervalHull.jl") end
    @time @testset "LazySets.concrete_convex_hull" begin include("unit_convex_hull.jl") end
    @time @testset "LazySets.Translation" begin include("unit_Translation.jl") end
    @time @testset "LazySets.AffineMap" begin include("unit_AffineMap.jl") end

    # ======================
    # Testing set interfaces
    # ======================
    @time @testset "LazySets.CompactSet" begin include("unit_CompactSet.jl") end

    # =========================================================
    # Testing other set types that do not inherit from LazySet
    # =========================================================
    @time @testset "LazySets.Complement" begin include("unit_Complement.jl") end
    @time @testset "LazySets.PolynomialZonotope" begin include("unit_PolynomialZonotope.jl") end
    @time @testset "LazySets.Rectification" begin include("unit_Rectification.jl") end
    @time @testset "LazySets.UnionSet" begin include("unit_UnionSet.jl") end

    # ===================
    # Concrete operations
    # ===================
    @time @testset "LazySets.samples" begin include("unit_samples.jl") end
    @time @testset "LazySets.area" begin include("unit_area.jl") end

    # =================================================================
    # Algorithms for approximation of convex sets using support vectors
    # =================================================================
    @time @testset "LazySets.Approximations.overapproximation" begin include("unit_overapproximate.jl") end
    @time @testset "LazySets.Approximations.underapproximation" begin include("unit_underapproximate.jl") end
    @time @testset "LazySets.Approximations.template_directions" begin include("unit_template_directions.jl") end
    @time @testset "LazySets.Approximations.box_approximation" begin include("unit_box_approximation.jl") end
    @time @testset "LazySets.Approximations.ballinf_approximation" begin include("unit_ballinf_approximation.jl") end
    @time @testset "LazySets.Approximations.radiusdiameter" begin include("unit_radiusdiameter.jl") end
    @time @testset "LazySets.Approximations.decompose" begin include("unit_decompose.jl") end
    @time @testset "LazySets.Approximations.distance" begin include("unit_distance.jl") end
    @time @testset "LazySets.Approximations.hausdorff_distance" begin include("unit_hausdorff_distance.jl") end

    # ========================
    # Testing method ambiguity
    # ========================
    include("check_method_ambiguity_binary.jl")
    @time @testset "LazySets.binary_operations" begin include("unit_binary_operations.jl") end

    # ====================================
    # Testing common API of all interfaces
    # (must be the last test because it
    #  loads Polyhedra.jl)
    # ====================================
    include("check_method_implementation.jl")
    @time @testset "LazySets.interfaces" begin include("unit_interfaces.jl") end
end

if test_suite_plotting
    import Plots
    using Plots: plot

    # fix namespace conflicts with Plots
    using LazySets: center, translate

    @time @testset "LazySets.plotting" begin include("unit_plot.jl") end
end

if test_suite_doctests
    using Documenter
    include("../docs/init.jl")
    @time @testset "LazySets.doctests" begin doctest(LazySets) end
end


text:
# Copyright 2019 Tobias Frilling
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
    gini_coeff(data::Vector{<:Real})

Calculates the Gini coefficient for the given data.
"""
function gini_coeff(data::Vector{<:Real})
    @assert all(x->x >= 0, data)
    y = sort(data)
    n = length(y)
    # cargo-culted from Wikipedia
    2 * sum([i * y[i] for i in 1:n]) / (n * sum(y)) - (n + 1) / n
end

"""
    gini_coeff(curve::ABCcurve)

Calculates the Gini coefficient for the given ABC curve.
"""
function gini_coeff(curve::ABCcurve)
    yield = curve.yield
    step::Float64 = curve.effort.step
    # simplified trapezoidal rule. The '- 1/2 * 2' is for the gini coefficient
    # this only works for equidistant steps
    ((sum(yield) - yield[end] / 2) * step - 1/2) * 2
end


text:
"""
    TransformerLogger(f, logger)
Preprocesses log messages, using the function `f`, before passing them to the
`logger` that is wrapped.
This can be used, for example, to truncate a log message.
to conditionally change the log level of logs from a given module
(which depending on the wrappped `logger`, might cause the message to be dropped).

The transforming function `f` is given a named tuple with the fields:
`(level, message, _module, group, id, file, line, kwargs)`
and should return the same.
See [`LoggingExtras.handle_message_args`](@ref) for more information on what each is.
"""
struct TransformerLogger{T<:AbstractLogger, F} <: AbstractLogger
    transform::F
    logger::T
end

function handle_message(transformer::TransformerLogger, args...; kwargs...)
    log_args = handle_message_args(args...; kwargs...)
    new_log_args = transformer.transform(log_args)

    args = Tuple(new_log_args)[1:end-1]
    kwargs = new_log_args.kwargs

    if comp_handle_message_check(transformer.logger, args...; kwargs...)
        handle_message(transformer.logger, args...; kwargs...)
    end
end

shouldlog(transformer::TransformerLogger, args...) = true

min_enabled_level(transformer::TransformerLogger) = BelowMinLevel

catch_exceptions(transformer::TransformerLogger) = catch_exceptions(transformer.logger)


text:
module TSPSolver

using Bonobo
using Graphs
using MetaGraphs
using DataStructures

include("utils.jl")
include("lb.jl")
include("heuristic.jl")

struct LONGEST_EDGE <: Bonobo.AbstractBranchStrategy end

mutable struct Node <: AbstractNode
    std :: BnBNodeInfo
    tour :: Union{Nothing, Vector{Int}}
    mst :: Union{Nothing, Vector{Edge{Int}}}
    fixed_edges :: Vector{Tuple{Int,Int}}
    disallowed_edges :: Vector{Tuple{Int,Int}}
end

struct Root
    g::MetaGraph
    cost::Matrix{Float64}
    function Root(g::MetaGraph)
        cost = zeros((nv(g), nv(g)))
        for i in 1:nv(g)
            for j in i+1:nv(g)
                weight = get_prop(g, i, j, :weight)
                cost[i,j] = weight
                cost[j,i] = weight
            end
        end
        new(g,cost)
    end
end

function Bonobo.get_branching_indices(root::Root)
    return 1:1
end

function Bonobo.evaluate_node!(tree::BnBTree{Node, Root}, node::Node)
    root = deepcopy(tree.root)
    extra_cost = 0.0
    for fixed_edge in node.fixed_edges
        extra_cost += fix_edge!(root, fixed_edge)
    end
    for disallowed_edge in node.disallowed_edges
        disallow_edge!(root, disallowed_edge)
    end

    # @show node.fixed_edges
    # @show node.disallowed_edges

    mst, lb = get_optimized_1tree(root.cost; runs=50)
    # @time mst, lb = get_1tree(weights(g))
    tour, ub = greedy(root.g)
    if tour !== nothing
        tour, ub = two_opt(tour, ub, root.cost)
    end
    lb += extra_cost
    ub += extra_cost
    node.mst = mst
    node.tour = tour
    # no tour can exist
    if isinf(ub)
        return NaN, NaN
    end
   # @show lb
   # @show ub
   # @show tree.incumbent
   # @show tree.lb
   # @show length(tree.nodes)
   # println("======================================================")
    return lb, ub
end

function Bonobo.get_relaxed_values(tree::BnBTree{Node, Root}, node::Node)
    return node.tour
end

function Bonobo.get_branching_variable(tree::BnBTree{Node, Root}, ::LONGEST_EDGE, node::Node)
    longest_len = 0.0
    longest_edge = -1
    tour_edges = Set{Tuple{Int, Int}}()
    longest_is_tour_edge = true
    if node.tour !== nothing
        for i in 1:length(node.tour)
            push!(tour_edges, extrema((node.tour[i], node.tour[mod1(i+1, end)])))
        end
    end
    for edge in node.mst
        edge_tpl = extrema((edge.src, edge.dst))
        if !(edge_tpl in node.fixed_edges) && !(edge_tpl in node.disallowed_edges)
            len = get_prop(tree.root.g, edge_tpl..., :weight)
            is_tour_edge = edge_tpl in tour_edges
            if (len > longest_len && !is_tour_edge) || longest_is_tour_edge
                longest_edge = edge
                longest_len = len
                longest_is_tour_edge = is_tour_edge
            end
        end
    end
    return longest_edge
end

function Bonobo.get_branching_nodes_info(tree::BnBTree{Node, Root}, node::Node, branching_edge::Edge)
    nodes_info = NamedTuple[]
    new_fixed_edges = deepcopy(node.fixed_edges)
    push!(new_fixed_edges, extrema((branching_edge.src, branching_edge.dst)))
    push!(nodes_info, (
        tour = nothing,
        mst = nothing,
        fixed_edges = new_fixed_edges,
        disallowed_edges = deepcopy(node.disallowed_edges),
    ))

    new_disallowed_edges = deepcopy(node.disallowed_edges)
    push!(new_disallowed_edges, extrema((branching_edge.src, branching_edge.dst)))
    push!(nodes_info, (
        tour = nothing,
        mst = nothing,
        fixed_edges = deepcopy(node.fixed_edges),
        disallowed_edges = new_disallowed_edges,
    ))
    return nodes_info
end

function optimize(input_path)
    g = simple_parse_tsp(input_path)

    bnb_model = Bonobo.initialize(
        traverse_strategy = Bonobo.BFS(),
        branch_strategy = LONGEST_EDGE(),
        Node = Node,
        root = Root(g),
        sense = :Min,
        Value = Vector{Int},
    )

    Bonobo.set_root!(bnb_model, (
        tour = nothing,
        mst = nothing,
        fixed_edges = Vector{Tuple{Int,Int}}(),
        disallowed_edges = Vector{Tuple{Int,Int}}()
    ))

    Bonobo.optimize!(bnb_model)

    return bnb_model
end

end


text:
import GPLinearODEMaker.powers_of_negative_one

"""
    pp_kernel(hyperparameters, δ, dorder; shift_ind=0)

Created by kernel_coder(). Requires 1 hyperparameters.
Likely created using pp_kernel_base() as an input.
Use with include("src/kernels/pp_kernel.jl").

# Arguments
- `hyperparameters::Vector`: The hyperparameter values. For this kernel, they should be `["λ"]`
- `δ::Real`: The difference between the inputs (e.g. `t1 - t2`)
- `dorder::Vector{<:Integer}`: How many times to differentiate with respect to the inputs and the `hyperparameters` (e.g. `dorder=[0, 1, 0, 2]` would correspond to differentiating once w.r.t the second input and twice w.r.t `hyperparameters[2]`)
- `shift_ind::Integer=0`: If changed, the index of which hyperparameter is the `δ` shifting one
"""
function pp_kernel(
    hyperparameters::Vector{<:Real},
    δ::Real,
    dorder::Vector{<:Integer};
    shift_ind::Integer=0)

    include_shift = shift_ind!=0
    @assert length(hyperparameters)==1+Int(include_shift) "hyperparameters is the wrong length"
    dorder_len = 3 + Int(include_shift)
    @assert length(dorder)==dorder_len "dorder is the wrong length"
    @assert maximum(dorder) < 3 "No more than two derivatives for each time or hyperparameter can be calculated"

    dorder2 = dorder[2]
    dorder[2] += dorder[1]

    if include_shift
        dorder[2] += dorder[shift_ind+2]
        δ += hyperparameters[shift_ind]
        hyperparameters_view = view(hyperparameters, 1:2 .!= shift_ind)
        dorder_view = view(dorder, 2:dorder_len)
        dorder_view = view(dorder_view, 1:(dorder_len-1) .!= (shift_ind+1))
    else
        hyperparameters_view = hyperparameters
        dorder_view = view(dorder, 2:dorder_len)
    end

    λ = hyperparameters_view[1]

    if abs(δ) > λ
        dorder[2] = dorder2
        return 0
    end

    dabs_δ = powers_of_negative_one(δ < 0)  # store derivative of abs()
    abs_δ = abs(δ)

    if dorder_view==[6, 2]
        func = -403200*abs_δ*dabs_δ^6/λ^9 + 1209600*dabs_δ^6*(1 - abs_δ/λ)/λ^8 - 21600*dabs_δ^5*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^7 + 7200*dabs_δ^5*(-5*dabs_δ/λ^2 - 32*abs_δ*dabs_δ/λ^3)/λ^6 - 720*dabs_δ^5*(10*dabs_δ/λ^3 + 96*abs_δ*dabs_δ/λ^4)/λ^5
    end

    if dorder_view==[5, 2]
        func = -19200*abs_δ^2*dabs_δ^5/λ^9 - 288000*dabs_δ^5*(1 - abs_δ/λ)^2/λ^7 - 120*(48*abs_δ^2/λ^4 + 10*abs_δ/λ^3)*dabs_δ^5/λ^5 + 1200*(-16*abs_δ^2/λ^3 - 5*abs_δ/λ^2)*dabs_δ^5/λ^6 - 3600*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*dabs_δ^5/λ^7 + 230400*abs_δ*dabs_δ^5*(1 - abs_δ/λ)/λ^8 - 6000*abs_δ*dabs_δ^4*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^7 + 1200*abs_δ*dabs_δ^4*(-5*dabs_δ/λ^2 - 32*abs_δ*dabs_δ/λ^3)/λ^6 + 12000*dabs_δ^4*(1 - abs_δ/λ)*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^6 - 4800*dabs_δ^4*(1 - abs_δ/λ)*(-5*dabs_δ/λ^2 - 32*abs_δ*dabs_δ/λ^3)/λ^5 + 600*dabs_δ^4*(1 - abs_δ/λ)*(10*dabs_δ/λ^3 + 96*abs_δ*dabs_δ/λ^4)/λ^4
    end

    if dorder_view==[4, 2]
        func = 38400*dabs_δ^4*(1 - abs_δ/λ)^3/λ^6 + 11520*abs_δ^2*dabs_δ^4*(1 - abs_δ/λ)/λ^8 - 57600*abs_δ*dabs_δ^4*(1 - abs_δ/λ)^2/λ^7 - 480*abs_δ^2*dabs_δ^3*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^7 - 2880*dabs_δ^3*(1 - abs_δ/λ)^2*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^5 + 1440*dabs_δ^3*(1 - abs_δ/λ)^2*(-5*dabs_δ/λ^2 - 32*abs_δ*dabs_δ/λ^3)/λ^4 - 240*dabs_δ^3*(1 - abs_δ/λ)^2*(10*dabs_δ/λ^3 + 96*abs_δ*dabs_δ/λ^4)/λ^3 + 120*(48*abs_δ^2/λ^4 + 10*abs_δ/λ^3)*dabs_δ^4*(1 - abs_δ/λ)/λ^4 + 240*(-16*abs_δ^2/λ^3 - 5*abs_δ/λ^2)*abs_δ*dabs_δ^4/λ^6 - 960*(-16*abs_δ^2/λ^3 - 5*abs_δ/λ^2)*dabs_δ^4*(1 - abs_δ/λ)/λ^5 - 1200*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*abs_δ*dabs_δ^4/λ^7 + 2400*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*dabs_δ^4*(1 - abs_δ/λ)/λ^6 + 3840*abs_δ*dabs_δ^3*(1 - abs_δ/λ)*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^6 - 960*abs_δ*dabs_δ^3*(1 - abs_δ/λ)*(-5*dabs_δ/λ^2 - 32*abs_δ*dabs_δ/λ^3)/λ^5
    end

    if dorder_view==[3, 2]
        func = -2880*dabs_δ^3*(1 - abs_δ/λ)^4/λ^5 - 2880*abs_δ^2*dabs_δ^3*(1 - abs_δ/λ)^2/λ^7 + 7680*abs_δ*dabs_δ^3*(1 - abs_δ/λ)^3/λ^6 + 360*dabs_δ^2*(1 - abs_δ/λ)^3*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^4 - 240*dabs_δ^2*(1 - abs_δ/λ)^3*(-5*dabs_δ/λ^2 - 32*abs_δ*dabs_δ/λ^3)/λ^3 + 60*dabs_δ^2*(1 - abs_δ/λ)^3*(10*dabs_δ/λ^3 + 96*abs_δ*dabs_δ/λ^4)/λ^2 - 60*(48*abs_δ^2/λ^4 + 10*abs_δ/λ^3)*dabs_δ^3*(1 - abs_δ/λ)^2/λ^3 + 360*(-16*abs_δ^2/λ^3 - 5*abs_δ/λ^2)*dabs_δ^3*(1 - abs_δ/λ)^2/λ^4 - 120*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*abs_δ^2*dabs_δ^3/λ^7 - 720*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*dabs_δ^3*(1 - abs_δ/λ)^2/λ^5 + 360*abs_δ^2*dabs_δ^2*(1 - abs_δ/λ)*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^6 - 1080*abs_δ*dabs_δ^2*(1 - abs_δ/λ)^2*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^5 + 360*abs_δ*dabs_δ^2*(1 - abs_δ/λ)^2*(-5*dabs_δ/λ^2 - 32*abs_δ*dabs_δ/λ^3)/λ^4 - 240*(-16*abs_δ^2/λ^3 - 5*abs_δ/λ^2)*abs_δ*dabs_δ^3*(1 - abs_δ/λ)/λ^5 + 960*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*abs_δ*dabs_δ^3*(1 - abs_δ/λ)/λ^6
    end

    if dorder_view==[2, 2]
        func = 96*dabs_δ^2*(1 - abs_δ/λ)^5/λ^4 + 320*abs_δ^2*dabs_δ^2*(1 - abs_δ/λ)^3/λ^6 - 480*abs_δ*dabs_δ^2*(1 - abs_δ/λ)^4/λ^5 - 20*dabs_δ*(1 - abs_δ/λ)^4*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^3 + 20*dabs_δ*(1 - abs_δ/λ)^4*(-5*dabs_δ/λ^2 - 32*abs_δ*dabs_δ/λ^3)/λ^2 - 10*dabs_δ*(1 - abs_δ/λ)^4*(10*dabs_δ/λ^3 + 96*abs_δ*dabs_δ/λ^4)/λ + 20*(48*abs_δ^2/λ^4 + 10*abs_δ/λ^3)*dabs_δ^2*(1 - abs_δ/λ)^3/λ^2 - 80*(-16*abs_δ^2/λ^3 - 5*abs_δ/λ^2)*dabs_δ^2*(1 - abs_δ/λ)^3/λ^3 + 120*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*dabs_δ^2*(1 - abs_δ/λ)^3/λ^4 - 120*abs_δ^2*dabs_δ*(1 - abs_δ/λ)^2*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^5 + 160*abs_δ*dabs_δ*(1 - abs_δ/λ)^3*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^4 - 80*abs_δ*dabs_δ*(1 - abs_δ/λ)^3*(-5*dabs_δ/λ^2 - 32*abs_δ*dabs_δ/λ^3)/λ^3 + 120*(-16*abs_δ^2/λ^3 - 5*abs_δ/λ^2)*abs_δ*dabs_δ^2*(1 - abs_δ/λ)^2/λ^4 + 120*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*abs_δ^2*dabs_δ^2*(1 - abs_δ/λ)/λ^6 - 360*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*abs_δ*dabs_δ^2*(1 - abs_δ/λ)^2/λ^5
    end

    if dorder_view==[1, 2]
        func = (1 - abs_δ/λ)^5*(10*dabs_δ/λ^3 + 96*abs_δ*dabs_δ/λ^4) + 20*abs_δ^2*(1 - abs_δ/λ)^3*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^4 - 10*abs_δ*(1 - abs_δ/λ)^4*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^3 + 10*abs_δ*(1 - abs_δ/λ)^4*(-5*dabs_δ/λ^2 - 32*abs_δ*dabs_δ/λ^3)/λ^2 - 5*(48*abs_δ^2/λ^4 + 10*abs_δ/λ^3)*dabs_δ*(1 - abs_δ/λ)^4/λ + 10*(-16*abs_δ^2/λ^3 - 5*abs_δ/λ^2)*dabs_δ*(1 - abs_δ/λ)^4/λ^2 - 10*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*dabs_δ*(1 - abs_δ/λ)^4/λ^3 - 40*(-16*abs_δ^2/λ^3 - 5*abs_δ/λ^2)*abs_δ*dabs_δ*(1 - abs_δ/λ)^3/λ^3 - 60*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*abs_δ^2*dabs_δ*(1 - abs_δ/λ)^2/λ^5 + 80*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*abs_δ*dabs_δ*(1 - abs_δ/λ)^3/λ^4
    end

    if dorder_view==[0, 2]
        func = (48*abs_δ^2/λ^4 + 10*abs_δ/λ^3)*(1 - abs_δ/λ)^5 + 10*(-16*abs_δ^2/λ^3 - 5*abs_δ/λ^2)*abs_δ*(1 - abs_δ/λ)^4/λ^2 + 20*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*abs_δ^2*(1 - abs_δ/λ)^3/λ^4 - 10*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*abs_δ*(1 - abs_δ/λ)^4/λ^3
    end

    if dorder_view==[6, 1]
        func = 28800*abs_δ*dabs_δ^6/λ^8 - 172800*dabs_δ^6*(1 - abs_δ/λ)/λ^7 + 3600*dabs_δ^5*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^6 - 720*dabs_δ^5*(-5*dabs_δ/λ^2 - 32*abs_δ*dabs_δ/λ^3)/λ^5
    end

    if dorder_view==[5, 1]
        func = 48000*dabs_δ^5*(1 - abs_δ/λ)^2/λ^6 - 120*(-16*abs_δ^2/λ^3 - 5*abs_δ/λ^2)*dabs_δ^5/λ^5 + 600*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*dabs_δ^5/λ^6 - 19200*abs_δ*dabs_δ^5*(1 - abs_δ/λ)/λ^7 + 600*abs_δ*dabs_δ^4*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^6 - 2400*dabs_δ^4*(1 - abs_δ/λ)*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^5 + 600*dabs_δ^4*(1 - abs_δ/λ)*(-5*dabs_δ/λ^2 - 32*abs_δ*dabs_δ/λ^3)/λ^4
    end

    if dorder_view==[4, 1]
        func = -7680*dabs_δ^4*(1 - abs_δ/λ)^3/λ^5 + 5760*abs_δ*dabs_δ^4*(1 - abs_δ/λ)^2/λ^6 + 720*dabs_δ^3*(1 - abs_δ/λ)^2*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^4 - 240*dabs_δ^3*(1 - abs_δ/λ)^2*(-5*dabs_δ/λ^2 - 32*abs_δ*dabs_δ/λ^3)/λ^3 + 120*(-16*abs_δ^2/λ^3 - 5*abs_δ/λ^2)*dabs_δ^4*(1 - abs_δ/λ)/λ^4 + 120*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*abs_δ*dabs_δ^4/λ^6 - 480*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*dabs_δ^4*(1 - abs_δ/λ)/λ^5 - 480*abs_δ*dabs_δ^3*(1 - abs_δ/λ)*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^5
    end

    if dorder_view==[3, 1]
        func = 720*dabs_δ^3*(1 - abs_δ/λ)^4/λ^4 - 960*abs_δ*dabs_δ^3*(1 - abs_δ/λ)^3/λ^5 - 120*dabs_δ^2*(1 - abs_δ/λ)^3*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^3 + 60*dabs_δ^2*(1 - abs_δ/λ)^3*(-5*dabs_δ/λ^2 - 32*abs_δ*dabs_δ/λ^3)/λ^2 - 60*(-16*abs_δ^2/λ^3 - 5*abs_δ/λ^2)*dabs_δ^3*(1 - abs_δ/λ)^2/λ^3 + 180*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*dabs_δ^3*(1 - abs_δ/λ)^2/λ^4 + 180*abs_δ*dabs_δ^2*(1 - abs_δ/λ)^2*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^4 - 120*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*abs_δ*dabs_δ^3*(1 - abs_δ/λ)/λ^5
    end

    if dorder_view==[2, 1]
        func = -32*dabs_δ^2*(1 - abs_δ/λ)^5/λ^3 + 80*abs_δ*dabs_δ^2*(1 - abs_δ/λ)^4/λ^4 + 10*dabs_δ*(1 - abs_δ/λ)^4*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^2 - 10*dabs_δ*(1 - abs_δ/λ)^4*(-5*dabs_δ/λ^2 - 32*abs_δ*dabs_δ/λ^3)/λ + 20*(-16*abs_δ^2/λ^3 - 5*abs_δ/λ^2)*dabs_δ^2*(1 - abs_δ/λ)^3/λ^2 - 40*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*dabs_δ^2*(1 - abs_δ/λ)^3/λ^3 - 40*abs_δ*dabs_δ*(1 - abs_δ/λ)^3*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^3 + 60*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*abs_δ*dabs_δ^2*(1 - abs_δ/λ)^2/λ^4
    end

    if dorder_view==[1, 1]
        func = (1 - abs_δ/λ)^5*(-5*dabs_δ/λ^2 - 32*abs_δ*dabs_δ/λ^3) + 5*abs_δ*(1 - abs_δ/λ)^4*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^2 - 5*(-16*abs_δ^2/λ^3 - 5*abs_δ/λ^2)*dabs_δ*(1 - abs_δ/λ)^4/λ + 5*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*dabs_δ*(1 - abs_δ/λ)^4/λ^2 - 20*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*abs_δ*dabs_δ*(1 - abs_δ/λ)^3/λ^3
    end

    if dorder_view==[0, 1]
        func = (-16*abs_δ^2/λ^3 - 5*abs_δ/λ^2)*(1 - abs_δ/λ)^5 + 5*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*abs_δ*(1 - abs_δ/λ)^4/λ^2
    end

    if dorder_view==[6, 0]
        func = 28800*dabs_δ^6*(1 - abs_δ/λ)/λ^6 - 720*dabs_δ^5*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^5
    end

    if dorder_view==[5, 0]
        func = -9600*dabs_δ^5*(1 - abs_δ/λ)^2/λ^5 - 120*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*dabs_δ^5/λ^5 + 600*dabs_δ^4*(1 - abs_δ/λ)*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^4
    end

    if dorder_view==[4, 0]
        func = 1920*dabs_δ^4*(1 - abs_δ/λ)^3/λ^4 - 240*dabs_δ^3*(1 - abs_δ/λ)^2*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^3 + 120*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*dabs_δ^4*(1 - abs_δ/λ)/λ^4
    end

    if dorder_view==[3, 0]
        func = -240*dabs_δ^3*(1 - abs_δ/λ)^4/λ^3 + 60*dabs_δ^2*(1 - abs_δ/λ)^3*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ^2 - 60*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*dabs_δ^3*(1 - abs_δ/λ)^2/λ^3
    end

    if dorder_view==[2, 0]
        func = 16*dabs_δ^2*(1 - abs_δ/λ)^5/λ^2 - 10*dabs_δ*(1 - abs_δ/λ)^4*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2)/λ + 20*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*dabs_δ^2*(1 - abs_δ/λ)^3/λ^2
    end

    if dorder_view==[1, 0]
        func = (1 - abs_δ/λ)^5*(5*dabs_δ/λ + 16*abs_δ*dabs_δ/λ^2) - 5*(1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*dabs_δ*(1 - abs_δ/λ)^4/λ
    end

    if dorder_view==[0, 0]
        func = (1 + 8*abs_δ^2/λ^2 + 5*abs_δ/λ)*(1 - abs_δ/λ)^5
    end

    dorder[2] = dorder2  # resetting dorder[2]
    return powers_of_negative_one(dorder2) * float(func)  # correcting for amount of t2 derivatives

end

return pp_kernel, 1  # the function handle and the number of kernel hyperparameters

max_stars_repo_name: christiangil/GPLinearODEMaker

max_issues_repo_name: christiangil/GPLinearODEMaker

max_forks_repo_name: christiangil/GPLinearODEMaker


text:
module CellQuadraturesTests

using Test
using Gridap

p = Point(1.0,1.1)
l = 10

ref_quad = TensorProductQuadrature(orders=(5,4))
quad2 = ConstantCellQuadrature(ref_quad,l)

coo_cell = coordinates(quad2)
wei_cell = weights(quad2)

coo = coordinates(ref_quad)
wei = weights(ref_quad)

for i in 1:quad2.length
  @test coo_cell[1] == coo
  @test wei_cell[1] == wei
end

@test isa(quad2,CellQuadrature)

@test string(quad2) == "CellQuadrature object"
s = "CellQuadrature object:\n dim: 2\n ncells: 10"
@test sprint(show,"text/plain",quad2) == s

end # module CellQuadraturesTests


text:
module BenchHomology

using BenchmarkTools
using Ripserer

include(joinpath(@__DIR__, "utils.jl"))

benchmarks = [
    (file="klein200.pts", dim_max=1, threshold=1.7, filtration=Rips),
    (file="sphere100.pts", dim_max=2, threshold=nothing, filtration=Rips),
    (file="lena1024.dipha", dim_max=2, threshold=nothing, filtration=Cubical),
]

suite = BenchmarkGroup()

for b in benchmarks
    name = b.file
    data = load_data(joinpath(@__DIR__, "data", b.file))
    if isnothing(b.threshold)
        filtration = b.filtration(data)
    else
        filtration = b.filtration(data; threshold=b.threshold)
    end
    suite[name] = @benchmarkable(
        ripserer($filtration; dim_max=$(b.dim_max), alg=:involuted), samples = 2
    )
end

end
BenchHomology.suite


text:
@info "Starting SINDy tests"
@testset "SINDy" begin
    # Create a nonlinear pendulum
    function pendulum(u, p, t)
        x = u[2]
        y = -9.81sin(u[1]) - 0.1u[2]^3 -0.2*cos(u[1])
        return [x;y]
    end

    u0 = [0.99π; -1.0]
    tspan = (0.0, 20.0)
    dt = 0.3
    prob = ODEProblem(pendulum, u0, tspan)
    sol = solve(prob, Tsit5(), saveat = dt)

    # Create the differential data
    DX = similar(sol[:,:])
    for (i, xi) in enumerate(eachcol(sol[:,:]))
        DX[:,i] = pendulum(xi, [], 0.0)
    end

    # Create a basis
    @variables u[1:2]
    # Lots of polynomials
    polys = Operation[1]
    for i ∈ 1:5
        push!(polys, u.^i...)
        for j ∈ 1:i-1
            push!(polys, u[1]^i*u[2]^j)
        end
    end

    # And some other stuff
    h = [cos(u[1]); sin(u[1]); u[1]*u[2]; u[1]*sin(u[2]); u[2]*cos(u[2]); polys]

    basis = Basis(h, u)

    opt = STRRidge(1e-2)
    Ψ = SINDy(sol[:,:], DX, basis, opt, maxiter = 2000)
    @test_nowarn set_threshold!(opt, 1e-2)

    @test length(Ψ) == 2
    @test size(Ψ) == (2, )
    @test parameters(Ψ) ≈ [1.0; -0.2; -9.81; -0.1]
    W = get_coefficients(Ψ)
    @test size(W) == (length(basis), 2)
    @test_nowarn get_aicc(Ψ)
    @test sum(get_sparsity(Ψ)) == 4
    @test sum(get_error(Ψ)) < 1e-10
    @test Ψ(u0) ≈ pendulum(u0, nothing, nothing)

    # Simulate
    estimator = ODEProblem(dynamics(Ψ), u0, tspan, parameters(Ψ))
    sol_ = solve(estimator,Tsit5(), saveat = dt)
    @test sol[:,:] ≈ sol_[:,:]

    opt = ADMM(1e-2, 0.7)
    Ψ = SINDy(sol[:,:], DX, basis, opt, maxiter = 5000)
    @test_nowarn set_threshold!(opt, 1e-2)

    # Simulate
    estimator = ODEProblem(dynamics(Ψ), u0, tspan, parameters(Ψ))
    sol_2 = solve(estimator,Tsit5(), saveat = dt)
    @test norm(sol[:,:] - sol_2[:,:], 2) < 2e-1

    opt = SR3(1e-2, 1.0)
    Ψ = SINDy(sol[:,:], DX, basis, opt, maxiter = 5000)
    @test_nowarn set_threshold!(opt, 0.1)

    # Simulate
    estimator = ODEProblem(dynamics(Ψ), u0, tspan, parameters(Ψ))
    sol_3 = solve(estimator,Tsit5(), saveat = dt)
    @test norm(sol[:,:] - sol_3[:,:], 2) < 1e-1

    # Now use the threshold adaptation
    opt = STRRidge(1e-2)
    λs = exp10.(-7:0.1:-1)
    Ψ = SINDy(sol[:,:], DX[:, :], basis, λs, opt, maxiter = 100)
    estimator = ODEProblem(dynamics(Ψ), u0, tspan, parameters(Ψ))
    sol_4 = solve(estimator, Tsit5(), saveat = dt)
    @test norm(sol[:,:] - sol_4[:,:], 2) < 1e-1

    # Check for errors
    @test_nowarn SINDy(sol[:,:], DX[1,:], basis, λs, opt, maxiter = 1)
    @test_nowarn SINDy(sol[:, :], DX[1, :], basis, λs, opt, maxiter = 1, denoise = true, normalize = true)

    # Check with noise
    X = sol[:, :] + 1e-3*randn(size(sol[:,:])...)
    opt = SR3(1e-1, 10.0)
    Ψ = SINDy(X, DX, basis, λs, opt, maxiter = 10000, denoise = true, normalize = true)

    estimator = ODEProblem(dynamics(Ψ), u0, tspan, parameters(Ψ))
    sol_5 = solve(estimator,Tsit5(), saveat = dt)
    @test norm(sol[:,:] - sol_5[:,:], 2) < 5e-1

    # Check sparse_regression
    X .= Array(sol)
    opt = SR3(1e-1, 1.0)
    maxiter = 5000
    θ = basis(X)
    Ξ1, iters_1 = sparse_regression(X, DX, basis, parameters(basis), [], maxiter, opt, true, true, eps())
    Ξ2 = similar(Ξ1)
    Ξ3 = similar(Ξ1)
    iters_2 = sparse_regression!(Ξ2, X, DX, basis, parameters(basis), [], maxiter, opt, true, true, eps())
    iters_3 = sparse_regression!(Ξ3, θ, DX, maxiter, opt, true, true, eps())

    @test iters_1 == iters_2 == iters_3
    @test Ξ1 ≈ Ξ2 ≈ Ξ3
    @test isapprox(norm(Ξ1'*θ - DX,2), 10.18; atol = 1e-1)
end


text:
include("../../geometry/sphere.jl")
include("../../physics/fluid.jl")
include("../../physics/dipoleforce.jl")
include("../../utils/simulator.jl")
include("../../utils/plot.jl")

fluid = FluidParams( 1. )

rs = 0.1

dipoles = [
    DipoleForce( SphericalCoord( π/2 - 0.2, π ), 0, 1 ),
    DipoleForce( SphericalCoord( π/2, π), π/2, 1)
]

# Generating a frame is computationally expensive,
# but each step forward in differential equation is computationally cheap.
# Higher steps reduce errors in numerical integration.
numFrames = 10 # number of frames in animation,
numSteps = 100 # number of steps forward in each frame,
endTime = 1 # time till when we want to simulate
adaptivePosMove = rs/10

# do the simulation
simulatedObjects = [ dipoles ]

for i in 2:numFrames
    newState = stepForward( last(simulatedObjects), fluid, endTime/numFrames/numSteps, numSteps, adaptivePosMove)
    nvortices = newState
    append!(simulatedObjects, [ newState ] )
    print(i, "\n")
end

# do the plitting/animation
animate(
    simulatedObjects,
    plotLatLon = true,
    # fieldGenerator = (p, vortices, pointforces, dipoles) -> sum( map( s -> velocityField(p, s, fluid), vcat( vortices, pointforces, dipoles ) ) ),
    resolution = (1000,600),
    makeMovie = true,
    cameraPosition = Vec3f0(-2, 0, 0),
    cameraLookingAt = Vec3f0(-1, 0, 0),
)


text:
"""
In a graph provide access to neighbors, i.e. the *adjacent nodes*,
for a specified node.

!!! note

    [`BiAdjacentIndex`](@ref) always includes both directions, i.e.
    for an edge `(u, v)` it is `u ∈ neighbors[v]` and `v ∈ neighbors[u]`.
"""
struct BiAdjacentIndex
    idx::Vector{Int}
    pi::Vector{Int}
end

@inline Base.getindex(c::BiAdjacentIndex, j::Int) =
    view(c.pi, c.idx[j]:c.idx[j+1]-1)

"""
    num_edges(::BiAdjacentIndex)

Actual number of (undirected) edges (**not counting** both directions).
"""
num_edges(c::BiAdjacentIndex) = length(c.pi) ÷ 2

num_nodes(c::BiAdjacentIndex) = length(c.idx) - 1

Base.length(c::BiAdjacentIndex) = num_nodes(c)

"""
    BiAdjacentIndex(n, head, tail)

Construct an index for `n` nodes of the (undirected) neighbors given
by the edges `zip(head, tail)`.
If `head[i] <= 0` then edge `i` is excluded.
"""
function BiAdjacentIndex(n::Int, head::Vector{Int}, tail::Vector{Int})
    @assert length(head) == length(tail)
    BiAdjacentIndex(n, () -> zip(head, tail))
end

"""
    BiAdjacentIndex(n, edges)

Same but for edges like `[(1, 2), (2, 3)]`
"""
BiAdjacentIndex(n::Int, edges::Vector{Tuple{Int, Int}}) =
    BiAdjacentIndex(n, () -> edges)

BiAdjacentIndex(n::Int, edges::Vector{Pair{Int, Int}}) =
    BiAdjacentIndex(n, () -> edges)

"""
    BiAdjacentIndex(n, m, iter::Function)

Provide an iterator over the edges.
If the first node of an edge is zero, skip that edge.
"""
function BiAdjacentIndex(n::Int, iter::Function)
    local m::Int = 0
    idx = zeros(Int, n+1)
    for (h::Int, t::Int) in iter()
        h <= 0 && continue
        m += 1
        idx[h] += 1
        idx[t] += 1
    end
    acc = 1                        # accumulate degrees ==> positions
    deg_i = 0
    deg_ii = idx[1]
    for i = 1:n
        idx[i] = acc
        acc += deg_i
        deg_i, deg_ii = deg_ii, idx[i+1]
    end
    idx[n+1] = acc
    pi = Vector{Int}(undef, 2m)
    @assert(idx[end] + deg_i == 2m + 1,
            "idx[$(length(idx))]: $(idx[end] + deg_i) != $(2m + 1)")
    for (i, (u, v)) in enumerate(iter())
        u <= 0 && continue
        pi[idx[u+1]] = v
        idx[u+1] += 1
        pi[idx[v+1]] = u
        idx[v+1] += 1
    end
    @assert(idx[end] == 2m + 1, "$(idx[end]) vs $(2m + 1)")
    return BiAdjacentIndex(idx, pi)
end

import Base.==

==(a::BiAdjacentIndex, b::BiAdjacentIndex) =
    a.pi == b.pi && a.idx == b.idx

node_degrees(b::BiAdjacentIndex) =
    diff(b.idx)

max_degree(neighs::BiAdjacentIndex)::Int =
    maximum((neighs.idx[i+1] - neighs.idx[i] for i in 1:length(neighs.idx)-1))

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
include("euler/euler.jl")
using .Math: fast_min

function compute(grid::Array{Int,2})::Int
    edge = size(grid)[1]
    for i ∈ edge - 1:-1:1
        grid[edge, i] += grid[edge, i + 1]
        grid[i, edge] += grid[i + 1, edge]
    end
    for i ∈ edge - 1:-1:1 for j ∈ edge - 1:-1:1 grid[i, j] += fast_min(grid[i + 1, j], grid[i, j + 1]) end end
    return grid[1, 1]
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
# 3D rk4 run

arg_dict = Dict{Any, Any}(
"physics" => "Euler",
"run_type" => 1,
"jac_type" => 2,
"order" => 1,
"use_DG" => true,
"Flux_name" => "RoeFlux",
"IC_name" => "ICExp",
"SRCname" => "SRCExp",
"dimensions" => 3,
"numBC" => 1,
"BC1" => [ 0, 1, 2, 3, 4, 5],
"BC1_name" => "ExpBC",
"delta_t" => 0.00005,
"t_max" => 500.000,
"smb_name" => "SRCMESHES/cube_benchmarksmall.smb",
"dmg_name" => ".null",
"res_abstol" => 1e-8,
"res_reltol" => 1e-20,
"step_tol" => 1e-10,
"itermax" => 3000,
"output_freq" => 100,
"step_tol" => -1.0,
"write_timing" => true,
"do_postproc" => true,
"exact_soln_func" => "ICExp",
"solve" => true,
)


text:
"""
    restrict(refspace, element1, element2)

Computes the restriction of a set of local shape functions on `element1` as linear
combinations of the set of local shape functions on `element2`. More precisely `restrict`
returns an `NxM` matrix `P` such that the `i`-th local shape ``g_i`` function on element2
can be written as:

``g_i = sum_{j=1}^{M} P_{ij} f_j``
"""
function restrict end


text:
import Base: show, eltype, length, size, getindex, ==, (*), iszero
import Base: Vector, Matrix, adjoint, hash

import LinearAlgebra: dot
export HVector, HMatrix, dot, hash

"""
An `HVector` is a homogeneous vector; two `HVector`s are equal
iff they are nonzero multiples of each other.

`HVector(a,b,c,...)` or `HVector([a,b,c,...])` creates a new `HVector `
with coordinates `(a,b,c,...)`.
"""
struct HVector{T}
    data::Vector{T}
    function HVector(dat::Vector{T}) where {T}
        data = dat .// one(T)
        S = typeof(sum(data))
        new{S}(_canonical(data))
    end

end
HVector(x...) = HVector(collect(x))
HVector(x::HVector) = HVector(x.data)

"""
An `HMatrix` is a homogenous matrix; two of these are equal if they are
nonzero multiples of each other.

Given a two-dimensional array `A`, `HMatrix(A)` creates a new homogeneous
matrix.
"""
struct HMatrix{T}
    data::Matrix{T}
    function HMatrix(dat::AbstractMatrix{T}) where {T}
        data = dat .// one(T)
        S = typeof(sum(data))
        M = new{S}(_canonical(data))
        return M
    end
end
HMatrix(x::HMatrix) = HMatrix(x.data)

function _canonical(x)
    if length(x) == 0 || iszero(x)
        return x
    end
    idx = findlast(x .!= 0)
    a = x[idx]
    return x .// a
end

HObject = Union{HVector,HMatrix}

eltype(x::HObject) = eltype(x.data)
length(x::HObject) = length(x.data)
size(x::HObject) = size(x.data)
getindex(x::HObject, k...) = x.data[k...]
(==)(x::HObject, y::HObject) = x.data == y.data
iszero(x::HObject) = iszero(x.data)

"""
The `dot` produce of two `HVector`s is either `0`
(if they are orthogonal) or `1` (otherwise).
"""
function dot(x::HVector, y::HVector)::Int
    d = dot(x.data, y.data)
    if d == 0
        return 0
    end
    return 1
end

function show(io::IO, x::HVector)
    print(io, "HVector(")
    n = length(x.data)
    for j = 1:n-1
        print(io, "$(x.data[j]), ")
    end
    print(io, x.data[end], ")")
end

function show(io::IO, x::HMatrix)
    print(io, "HMatrix: ")
    print(io, x.data)
end

function (*)(A::AbstractArray{T,2}, v::HVector) where {T}
    Av = A * v.data
    return HVector(Av)
end

function (*)(A::AbstractArray{T,2}, B::HMatrix) where {T}
    AB = A * B.data
    return HMatrix(AB)
end

(*)(A::HMatrix, v::HVector) = HVector(A.data * v.data)
(*)(A::HMatrix, B::HMatrix) = HMatrix(A.data * B.data)

adjoint(A::HMatrix) = HMatrix(adjoint(A.data))

hash(x::HObject, h::UInt = UInt(0)) = hash(x.data, h)

Vector(x::HVector) = copy(x.data)
Matrix(x::HMatrix) = copy(x.data)