text:
function get_cm(C::ClimGrid, cm, surface, center_cs, ncolors)

    if isempty(cm)
      if C[10] == "pr" || C[10]=="huss" || C[10]=="clwvi"
          cm = cmocean.cm.deep

      elseif C[10]=="tasmax" || C[10]=="tasmin" || C[10]=="tas" || C[10]=="tmax" || C[10]=="tmin" || C[10] == "wbgtmean" || C[10] == "wbgtmax" || C[10]=="t2m" || C[10]=="tmean" || C[10]=="dc" || C[10]=="hfls"
          cm = "RdYlBu_r"

      elseif C[10]=="psl" || C[10]=="vp"
          cm = cmocean.cm.deep_r

      elseif C[10]=="ua" || C[10]=="uv"
          cm = cmocean.cm.balance

      elseif C[10]=="clt" # cloud fraction
          cm = "Blues_r"

      elseif C[10]=="rlut" # longwave radiation
          cm = cmocean.cm.balance

      elseif C[10]=="hurs" || C[10]=="evspsbl"# longwave radiation
          cm = cmocean.cm.haline_r
      else
        cm = "viridis"
      end
    end

    # overide colorscale if the user specify to center the colorscale
    if center_cs
        cm = "RdBu_r"
    end

    if surface == :contourf || surface == :contour
        N = 256
    elseif surface == :pcolormesh
        N = ncolors
    end
    cmap = mpl.cm.get_cmap(cm)
    colorlist = cmap(range(0, stop=1, length=N))
    # #
    cm = mpl.colors.LinearSegmentedColormap.from_list("cm_custom", colorlist, N)

    return cm

end

"""
    getcslimits(caxis, data, C)

Returns minimum and maximum values of the colorscale axis. Used internally by [`mapclimgrid`](@ref).
"""
function getcslimits(caxis, data, center_cs)

    if !isempty(caxis)
        vmin = caxis[1]
        vmax = caxis[2]
    else
        vmin=minimum(data[.!isnan.(data)])
        vmax=maximum(data[.!isnan.(data)])
    end

    if center_cs # Hack, we want to center the divergent colorscale
        vmin = -vmax
    end

    return vmin, vmax

end

"""
    titledef(C::ClimGrid)

Returns the title. Used internally by [`mapclimgrid`](@ref).
"""
function titledef(C::ClimGrid)
    if ndims(C[1]) > 2

        if typeof((C[1][Axis{:time}][1])) <: AbstractCFDateTime || typeof((C[1][Axis{:time}][1])) == DateTime
            begYear = string(Dates.year(C[1][Axis{:time}][1]))
            endYear = string(Dates.year(C[1][Axis{:time}][end]))
        elseif typeof((C[1][Axis{:time}][1])) == Dates.Year
            begYear = string(C[1][Axis{:time}][1])[1:4]
            endYear = string(C[1][Axis{:time}][end])[1:4]
        elseif typeof((C[1][Axis{:time}][1])) == Int
            begYear = string(C[1][Axis{:time}][1])[1:4]
            endYear = string(C[1][Axis{:time}][end])[1:4]
        end
        titlestr = string(C[3], " - ", C[5], " - ", C.variable, " - ", begYear, " - ", endYear)

    else
        titlestr = string(C[3], " - ", C[5], " - ", C.variable)
    end

    return titlestr

end

"""
    getunitslabel(C::ClimGrid)

Return verbose label for colorbar. Used internally by [`mapclimgrid`](@ref).
"""
function getunitslabel(C::ClimGrid)

    try
        standardname = C.varattribs["standard_name"]
    catch
        standardname = C.varattribs["long_name"]
    end

    units = Dict(["air_temperature" => "Air temperature",
    "specific_humidity" => "Specific humidity (%)",
    "precipitation" => "Precipitation",
    "precipitation_flux" => "Precipitation flux",
    "air_pressure_at_sea_level" => "Air pressure at sea level",
    "eastward_wind" => "Eastward wind",
    "daily maximum temperature" => "Air temperature"
    ])

    label = ""

    try
        label = string(units[standardname], " (", C[2], ")")
    catch
        label = C[2]
    end

    return label

end

function roundup(i,v)
    return round(Int, i/v)*v
end


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#
using DiscreteEvents

a = 1
b = 1
f(x, y) = x+y
g(x, y) = x == y
h(clk, x, y) = clk.time + x + y
incra(clk) = global a += 1

clk = Clock()

# 1st @event macro: timed events
@test repr(@macroexpand @event h(clk, a, b) at a) == ":(event!(clk, fun(h, clk, a, b), at, a))"
VERSION > v"1.5" && @test repr(@macroexpand @event h(clk, a, b) every a 10) == ":(event!(clk, fun(h, clk, a, b), every, a, n = 10))"

# 2nd @event macro: conditional events
@test repr(@macroexpand @event h(clk, a, b) g(x, y)) == ":(event!(clk, fun(h, clk, a, b), fun(g, x, y)))"
@test repr(@macroexpand @event h(clk, a, b) a==1) == ":(event!(clk, fun(h, clk, a, b), fun(==, a, 1)))"
@test repr(@macroexpand @event h(clk, a, b) a) == ":(event!(clk, fun(h, clk, a, b), a))"

@event incra(clk) every 1
@event stop!(clk) ()->a ≥ 2 && tau(clk) ≥ 5  # @event with an anonymous function condition
@run! clk 10
@test clk.time ≈ 5.01
@test clk.evcount == 5

@test repr(@macroexpand @periodic h(clk, a, b)) == ":(periodic!(clk, fun(h, clk, a, b)))"
@test repr(@macroexpand @periodic h(clk, a, b) 1) == ":(periodic!(clk, fun(h, clk, a, b), 1))"

@test repr(@macroexpand @process h(clk, a, b)) == ":(process!(clk, Prc(h, a, b)))"
@test repr(@macroexpand @process h(clk, a, b) 10) == ":(process!(clk, Prc(h, a, b), 10))"

@test repr(@macroexpand @delay clk a) == ":(delay!(clk, a))"
@test repr(@macroexpand @delay clk until 10) == ":(delay!(clk, until, 10))"

@test repr(@macroexpand @wait clk g(a,b)) == ":(wait!(clk, fun(g, a, b)))"

@test repr(@macroexpand @run! clk a) == ":(run!(clk, a))"

ext: jl


text:
# This file holds the definition of the functions pertaining to the Personnel
#   type.

# The functions of the Personnel type require the HistoryEvent and HistoryEntry
#   types.
requiredTypes = [ "personnel" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

# Load in the type aliases.
include( joinpath( typePath, "typeAliases.jl" ) )

# This function tests if a personnel member has the requested attribute.
export hasAttribute
function hasAttribute( person::Personnel, attr::AttributeType )
    return haskey( person.persData, Symbol( attr ) )
end # hasAttribute( person, attr )

# This funciton tests if a personnel member has a history of the requested
#   attribute.
# This test consists of test to see if the personnel member has a history
#   attribute at all, if it's the correct type (Dict{Symbol, History}), and the
#   history has an entry for the requested attribute.
export hasHistory
function hasHistory( person::Personnel, attr::AttributeType )
    return hasAttribute( person, :history ) &&
        isa( person[ :history ], Dict{Symbol, History} ) &&
        haskey( person[ :history ], Symbol( attr ) )
end  # hasHistory( person, attr )

# This function removes the respested attribute from the personnel member.
export removeAttribute!
function removeAttribute!( person::Personnel, attr::AttributeType )
    tmpAttr = Symbol( attr )

    if hasAttribute( person, tmpAttr )
        delete!( person.persData, tmpAttr )
    end  # if hasAttribute( person, tmpAttr )
end  # removeAttribute!( person, attr )

# This function adds a value to a vector attribute. If the attribute's value is
#   nothing, a vector will be created. If the attribute's value is not a vector,
#   nothing will happen.
export addValue
function addValue( person::Personnel, attr::AttributeType, value::String )
    tmpAttr = Symbol( attr )
    # If the attribute doesn't exist, if the attribute's value is nothing, or if
    #   the attribute is an empty Vector, turn it into an empty Vector{String}.
    if !hasAttribute( person, tmpAttr ) || ( person[ tmpAttr ] === nothing ) ||
        isempty( person[ tmpAttr ] )
        person[ tmpAttr ] = Vector{String}()
    end  # if !hasAttribute( person, tmpAttr ) ...

    # If the attribute isn't of Vector{String} type, do nothing.
    if !isa( person[ tmpAttr ], Vector{String} )
        return
    end  # if !isa( person[ tmpAttr ], ...

    if !in( value, person[ tmpAttr ] )
        push!( person[ tmpAttr ], value )
    end
end  # addValue( person, attr, value )

# This function displays a personnel record onto the entered IO stream.
export displayPersonnel
function displayPersonnel( io::IO, person::Personnel )
    linePrinted = false

    for attr in keys( person.persData )
        val = person[ attr ]

        if nothing !== val
            if linePrinted
                print( io, "\n" )
            else
                linePrinted = true
            end  # if linePrinted

            print( io, string( attr ), ": $(val)" )
        end  # if nothing !== val
    end  # for attr ...
end  # displayPersonnel( io, person )

# This function displays a personnel record onto the entered IO stream, with the
#   index key displayed first.
function displayPersonnel( io::IO, person::Personnel, idKey::AttributeType )
    tmpKey = Symbol( idKey )

    if hasAttribute( person, tmpKey )
        keyVal = person[ tmpKey ]
        println( io, string( tmpKey ), " (key): ",
            nothing === keyVal ? "Undefined" : keyVal )
    end  # if hasAttribute( person, tmpKey )

    for attr in keys( person.persData )
        val = person[ attr ]

        if ( tmpKey !== attr ) && ( nothing !== val )
            println( io, string( attr ), ": $val" )
        end  # if ( tmpKey !== attr ) ...
    end  # for attr ...
end  # displayPersonnel( io, person, idKey )

# This function displays a personnel record onto the entered IO stream, showing
#   the requested attributes.
function displayPersonnel( io::IO, person::Personnel, attrs::Array{Symbol} )
    for attr in attrs
        if hasAttribute( person, attr ) && ( nothing !== person[ attr ] )
            println( io, string( attr ), ": $(person[ attr ])" )
        end  # if hasField( person, attr )
    end  # for attr in attrs
end  # displayPersnnel( io, person, attrs )

# This function retrieves the requested attributes for a personnel record,
#   without needing to access the dictionary directly (syntactic sugar).
# If the personnel record does not contain the attributes, this function returns
#   nothing.
function Base.getindex( person::Personnel, attr::AttributeType )
    tmpAttr = Symbol( attr )

    if !hasAttribute( person, tmpAttr )
        return nothing
    end  # if !hasField( person, tmpAttr )

    return person.persData[ Symbol( tmpAttr ) ]
end  # Base.getindex( person, attr )

function Base.getindex( person::Personnel, attrs::Array{Symbol} )
    out = similar( attrs, Any )

    for ii in eachindex( attrs )
        out[ ii ] = person[ attrs[ ii ] ]
    end  # for ii in eachindex( attrs )

    return out
end  # Base.getindex( person, attrs )

# This function returns the value of the specified attribute at the given time.
# If the personnel record does not have a history for this attribute, the
#   attribute itself is returned if possible.
function Base.getindex( person::Personnel, attr::AttributeType,
    timestamp::T ) where T <: Real
    tmpAttr = Symbol( attr )
    # If the personnel record has no history for the attribute, return the
    #   regular person[ attribute ].
    if !hasHistory( person, tmpAttr )
        return person[ tmpAttr ]
    end  # if !hasHistory( person, tmpAttr )

    return person[ :history ][ tmpAttr ][ timestamp ]
end  # Base.getindex( person, attr, timestamp )

# This function sets the attribute for a personnel record, without needing to
#   access the dictionary directly (syntactic sugar). If the attribute does not
#   exist yet, it will be created.
function Base.setindex!( person::Personnel, data, attr::AttributeType )
    person.persData[ Symbol( attr ) ] = data
end  # Base.setindex!( person, data, attr )

# This function sets a historical value for the attribute for a personnel
#   record.
# If no history exists for this attribute, it will be created.
# If the history attribute already exists, but it is not of the correct type
#   (Dict{Symbol, History}), it will be overwritten.
function Base.setindex!( person::Personnel, data, attr::AttributeType,
    timestamp::T ) where T <: Real
    tmpAttr = Symbol( attr )

    # If the history attribute doesn't exist, or if it's the wrong type, create
    #   and (re)initialize it.
    if !hasAttribute( person, :history ) ||
        !isa( person[ :history ], Dict{Symbol, History} )
        person[ :history ] = Dict{Symbol, History}()
    end  # if !hasAttribute( person, :history ) || ...

    # If the history doesn't have an entry for the atrtibute, create it.
    if !haskey( person[ :history ], tmpAttr )
        person[ :history ][ tmpAttr ] = History( tmpAttr )
    end  # if !haskey( person[ :history ], tmpAttr )

    person[ :history ][ Symbol( attr ) ][ timestamp ] = data
end  # Base.setindex!( person, data, attr, timestamp )

function Base.show( io::IO, person::Personnel )
    displayPersonnel( io, person )
end  # Base.show( io, person )


text:
"""
$(DocStringExtensions.README)
"""
module Catalyst

using DocStringExtensions
using Reexport, ModelingToolkit
using ModelingToolkit: Symbolic, value, istree
@reexport using ModelingToolkit
import MacroTools
import Base: (==), merge!, merge
using Latexify, Requires

# as used in Catlab
const USE_GV_JLL = Ref(false)
function __init__()
    @require Graphviz_jll="3c863552-8265-54e4-a6dc-903eb78fde85" begin
      USE_GV_JLL[] = true
      let cfg = joinpath(Graphviz_jll.artifact_dir, "lib", "graphviz", "config6")
        if !isfile(cfg)
          Graphviz_jll.dot(path -> run(`$path -c`))
        end
      end
    end
  end

const ExprValues = Union{Expr,Symbol,Float64,Int}

include("expression_utils.jl")
include("reaction_network.jl")

# reaction network macro
export @reaction_network, @add_reactions

include("spatial_network.jl")
export @spatial_reaction_network
export graph_iterator

# functions to query network properties
include("networkapi.jl")
export species, params, reactions, speciesmap, paramsmap, numspecies, numreactions, numparams
export make_empty_network, addspecies!, addparam!, addreaction!
export dependants, dependents, substoichmat, prodstoichmat

# for Latex printing of ReactionSystems
include("latexify_recipes.jl")

# for making and saving graphs
import Base.Iterators: flatten
import DataStructures: OrderedDict
import Parameters: @with_kw_noshow
include("graphs.jl")
export Graph, savegraph

end # module


text:
# Copyright 2017 Stephen L. Smith and Frank Imeson
#
# 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.

#####################################################
#########  GTSP Utilities ###########################

""" tour type that stores the order array and the length of the tour
"""
mutable struct Tour
	tour::Array{Int64,1}
	cost::Int64
end

""" return the vertex before tour[i] on tour """
@inline function prev_tour(tour, i)
	i != 1 && return tour[i - 1]
	return tour[length(tour)]
end

######################################################
#############  Randomizing tour before insertions ####

""" some insertions break tie by taking first minimizer -- this
randomization helps avoid getting stuck choosing same minimizer """
function pivot_tour!(tour::Array{Int64,1})
	pivot = rand(1:length(tour))
	tour = [tour[pivot:end]; tour[1:pivot-1]]
end

function randomize_sets!(sets::Array{Any, 1}, sets_to_insert::Array{Int64, 1})
	for i in sets_to_insert
		shuffle!(sets[i])
	end
end

function findmember(num_vertices::Int64, sets::Array{Any, 1})
    """  create an array containing the set number for each vertex """
	member = zeros(Int64, num_vertices)
    num_verts = 0
    for i = 1:length(sets)
        set = sets[i]
        num_verts += length(set)
        for vertex in set
			if member[vertex] != 0
				error("vertex ", vertex, " belongs to more than one set")
			else
				member[vertex] = i
			end
        end
    end
    return member
end

struct Distsv
	set_vert::Array{Int64, 2}
	vert_set::Array{Int64,2}
	min_sv::Array{Int64, 2}
end

function set_vertex_dist(dist::Array{Int64, 2}, num_sets::Int, member::Array{Int64,1})
    """
	Computes the minimum distance between each set and each vertex
	Also compute the minimum distance from a set to a vertex, ignoring direction
	This is used in insertion to choose the next set.
	"""
    numv = size(dist, 1)
    dist_set_vert = typemax(Int64) * ones(Int64, num_sets, numv)
	mindist = typemax(Int64) * ones(Int64, num_sets, numv)
	dist_vert_set = typemax(Int64) * ones(Int64, numv, num_sets)

	for i = 1:numv
        for j = 1:numv
			set = member[j]
			if dist[j,i] < dist_set_vert[set, i]
				dist_set_vert[set, i] = dist[j,i]
			end
			if dist[j,i] < mindist[set, i]  # dist from set containing j to vertex i
				mindist[set, i] = dist[j, i]
			end
			set = member[i]
			if dist[j,i] < dist_vert_set[j, set]  # dist from j to set containing i
				dist_vert_set[j, set] = dist[j,i]
			end
			if dist[j,i] < mindist[set,j] # record as distance from set containing i to j
				mindist[set,j] = dist[j,i]
			end
		end
	end
    return Distsv(dist_set_vert, dist_vert_set, mindist)
end

function set_vertex_distance(dist::Array{Int64, 2}, sets::Array{Any, 1})
    """
	Computes the minimum distance between each set and each vertex
	"""
    numv = size(dist, 1)
    nums = length(sets)
    dist_set_vert = typemax(Int64) * ones(Int64, nums, numv)
	# dist_vert_set = typemax(Int64) * ones(Int64, numv, nums)
    for i = 1:nums
        for j = 1:numv
			for k in sets[i]
				newdist = min(dist[k, j], dist[j, k])
				dist_set_vert[i,j] > newdist && (dist_set_vert[i,j] = newdist)
			end
        end
    end
    return dist_set_vert
end

""" Find the set with the smallest number of vertices """
function min_set(sets::Array{Any, 1})
    min_size = length(sets[1])
	min_index = 1
    for i = 2:length(sets)
		set_size = length(sets[i])
		if set_size < min_size
			min_size = set_size
			min_index = i
		end
	end
	return min_index
end

############################################################
############ Trial acceptance criteria #####################

"""
decide whether or not to accept a trial based on simulated annealing criteria
"""
function accepttrial(trial_cost::Int64, current_cost::Int64, temperature::Float64)
    if trial_cost <= current_cost
        accept_prob = 2.0
	else
        accept_prob = exp((current_cost - trial_cost)/temperature)
    end
    return (rand() < accept_prob ? true : false)
end

"""
decide whether or not to accept a trial based on simple probability
"""
function accepttrial_noparam(trial_cost::Int64, current_cost::Int64, prob_accept::Float64)
    if trial_cost <= current_cost
        return true
	end
	return (rand() < prob_accept ? true : false)
end

###################################################
################ Tour checks ######################

"""  Compute the length of a tour  """
@inline function tour_cost(tour::Array{Int64,1}, dist::Array{Int64,2})
    tour_length = dist[tour[end], tour[1]]
    @inbounds for i in 1:length(tour)-1
    	tour_length += dist[tour[i], tour[i+1]]
    end
    return tour_length
end

"""
Checks if a tour is feasible in that it visits each set exactly once.
"""
function tour_feasibility(tour::Array{Int64,1}, membership::Array{Int64,1},
					      num_sets::Int64)
    length(tour) != num_sets && return false

    set_test = falses(num_sets)
    for v in tour
        set_v = membership[v]
		if set_test[set_v]
			return false  # a set is visited twice in the tour
		end
		set_test[set_v] = true
    end
    for visited_set in set_test
        !visited_set && return false
    end
    return true
end

#####################################################
#############  Incremental Shuffle ##################

@inline function incremental_shuffle!(a::AbstractVector, i::Int)
    j = i + floor(Int, rand() * (length(a) + 1 - i))
   	a[j], a[i] = a[i], a[j]
	return a[i]
end

""" rand_select for randomize over all minimizers """
@inline function rand_select(a::Array{Int64, 1}, val::Int)
	inds = Int[]
	@inbounds for i=1:length(a)
		a[i] == val && (push!(inds, i))
	end
	return rand(inds)
end


text:
using KhepriBlender
using Test

@testset "KhepriBlender.jl" begin
    # Write your tests here.
end


text:
"""
Generates data for plotting f on [a, b]
"""
function plt_data(f, a, b, numpts=1000)
    h = (b-a)/(numpts-1)
    y = zeros(1,numpts)
    x = collect(a:h:b)
    y = map(f, x)
    return x, y
end

"""
Displays graph of f on interval [a, b]
"""
function plt(f, a, b, numpts=100, label = "y")
    h = (b-a)/(numpts-1)
    y = zeros(1,numpts)
    x = collect(a:h:b)
    y = map(f, x)
    display(plot(x, y, label = label))
end

function plt!(f, a, b, numpts=100, label = "y")
    h = (b-a)/(numpts-1)
    y = zeros(1,numpts)
    x = collect(a:h:b)
    y = map(f, x)
    display(plot!(x, y, label = label))
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text: abstract type AbstractUncertainScalarKDE{T} <: AbstractEmpiricalValue end


text:
using Test
using LogicCircuits
import .LogicCircuits.IO:
   parse_sdd_file, CircuitFormatLines

@testset "SDD file parser tests" begin

  @test parse_sdd_file(zoo_sdd_file("random.sdd")) isa CircuitFormatLines

  circuit = load_logical_circuit(zoo_sdd_file("random.sdd"))

  @test circuit isa UnstLogicalΔ
  @test num_nodes(circuit) == 1676
  @test is_decomposable(circuit)
  @test !is_smooth(circuit)
  @test any(n -> n isa FalseNode, circuit)
  @test any(n -> n isa TrueNode, circuit)
  @test num_variables(circuit) == 30
  @test has_unique_literal_nodes(circuit)
  @test has_unique_constant_nodes(circuit)

  prop_circuit = propagate_constants(circuit)

  @test prop_circuit isa UnstLogicalΔ
  @test num_nodes(prop_circuit) == 1232
  @test is_decomposable(prop_circuit)
  @test !is_smooth(prop_circuit)
  @test !any(n -> n isa FalseNode, prop_circuit)
  @test !any(n -> n isa TrueNode, prop_circuit)
  @test num_variables(prop_circuit) == 30
  @test has_unique_literal_nodes(prop_circuit)
  @test has_unique_constant_nodes(prop_circuit)

  @test prop_circuit[end] === propagate_constants(prop_circuit)[end] # no new circuit created if no changes

  smooth_circuit = smooth(prop_circuit)

  @test smooth_circuit isa UnstLogicalΔ
  @test num_nodes(smooth_circuit) == 2056
  @test is_decomposable(smooth_circuit)
  @test is_smooth(smooth_circuit)
  @test !any(n -> n isa FalseNode, smooth_circuit)
  @test !any(n -> n isa TrueNode, smooth_circuit)
  @test num_variables(smooth_circuit) == 30
  @test has_unique_literal_nodes(smooth_circuit)
  @test has_unique_constant_nodes(smooth_circuit)

  @test smooth_circuit[end] === smooth(smooth_circuit)[end] # no new circuit created if no changes

  forgotten_circuit = forget(v -> (v > 16), circuit)

  @test forgotten_circuit isa UnstLogicalΔ
  @test num_nodes(forgotten_circuit) == 1648
  @test num_variables(forgotten_circuit) == 16
  @test is_decomposable(forgotten_circuit)
  @test !is_smooth(forgotten_circuit)
  @test has_unique_literal_nodes(forgotten_circuit)
  @test has_unique_constant_nodes(forgotten_circuit)

  @test forgotten_circuit[end] === forget(v -> (v > 16), forgotten_circuit)[end] # no new circuit created if no changes

  random_circuit = smooth(propagate_constants(forgotten_circuit))

  @test random_circuit isa UnstLogicalΔ
  @test num_nodes(random_circuit) == 1956
  @test is_decomposable(random_circuit)
  @test is_smooth(random_circuit)
  @test !any(n -> n isa FalseNode, random_circuit)
  @test !any(n -> n isa TrueNode, random_circuit)
  @test num_variables(random_circuit) == 16
  @test has_unique_literal_nodes(random_circuit)
  @test has_unique_constant_nodes(random_circuit)

end


text:
# DP formulation for independent set

# State is IntSet

"""Independent set transition function"""
function create_indepset_transition_function(g::Graph)
	return function indepset_transition_function(state::IntSet, variable::Int, value::Int)
		@assert(value == 0 || value == 1)
		new_state = IntSet(state)
		if value == 0
			delete!(new_state, variable)
		else # value == 1
			if !in(variable, state)
				return nothing # indicates false state
			end
			delete!(new_state, variable)
			setdiff!(new_state, out_neighbors(g, variable))
		end
		return new_state
	end
end

"""Independent set initial state"""
function create_indepset_initial_state(g::Graph)
	return IntSet(1:nv(g))
end

max_stars_repo_name: ctjandra/ddopt-ecp

max_issues_repo_name: ctjandra/ddopt-ecp

max_forks_repo_name: ctjandra/ddopt-ecp


text:
#################### Heidelberger and Welch Diagnostic ####################
"""
    heideldiag(x::Vector{T}; alpha::Real=0.05, eps::Real=0.1,
                etype=:imse, start::Integer=1, args...) where {T<:Real}
"""
function heideldiag(
    x::Vector{T};
    alpha::Real = 0.05,
    eps::Real = 0.1,
    etype = :imse,
    start::Integer = 1,
    args...,
) where {T<:Real}
    n = length(x)
    delta = trunc(Int, 0.10 * n)
    y = x[trunc(Int, n / 2):end]
    S0 = length(y) * mcse(y, etype; args...)^2
    i, pvalue, converged, ybar = 1, 1.0, false, NaN
    while i < n / 2
        y = x[i:end]
        m = length(y)
        ybar = mean(y)
        B = cumsum(y) - ybar * collect(1:m)
        Bsq = (B .* B) ./ (m * S0)
        I = sum(Bsq) / m
        pvalue = 1.0 - pcramer(I)
        converged = pvalue > alpha
        if converged
            break
        end
        i += delta
    end
    halfwidth = sqrt(2.0) * erfinv(1.0 - alpha) * mcse(y, etype; args...)
    passed = halfwidth / abs(ybar) <= eps
    [i + start - 2, converged, round(pvalue, digits = 4), ybar, halfwidth, passed]
end
"""
    heideldiag(c::AbstractChains; alpha::Real=0.05, eps::Real=0.1,
                etype=:imse, args...)

Compute the convergence diagnostic of Heidelberger and Welch for MCMC sampler output. The diagnostic is designed to assess convergence of posterior means estimated with autocorrelated samples and to determine whether a target degree of accuracy is achieved. A stationarity test is performed for convergence assessment by iteratively discarding 10% of the initial samples until the test p-value is non-significant and stationarity is concluded or until 50% have been discarded and stationarity is rejected, whichever occurs first. Then, a halfwidth test is performed by calculating the relative halfwidth of a posterior mean estimation interval as ``z_{1 - \\alpha / 2} \\hat{s} / |\\bar{\\theta}|``; where ``z`` is a standard normal quantile, ``\\hat{s}`` is the Monte Carlo standard error, and ``\\bar{\\theta}`` is the estimated posterior mean. If the relative halfwidth is greater than a target ratio, the test is rejected. Rejection of the stationarity or halfwidth test suggests that additional samples are needed.

Returns a `ChainSummary` type object with parameters contained in the rows of the `value` field, and numbers of burn-in sequences to discard, whether the stationarity tests are passed (1 = yes, 0 = no), their p-values (``p > \\alpha`` implies stationarity), posterior means, halfwidths of their ``(1 - \\alpha) 100\\%`` estimation intervals, and whether the halfwidth tests are passed (1 = yes, 0 = no) in the columns. Results are chain-specific.

* `c` : sampler output on which to perform calculations.

* `alpha` : significance level for evaluations of stationarity tests and calculations of relative estimation interval halfwidths.

* `eps` : target ratio for the relative halfwidths.

* `etype` : method for computing Monte Carlo standard errors. See `mcse()` for options.

* `args...` : additional arguments to be passed to the `etype` method.

* `x` : vector on which to perform calculations.

* `start` : ??

"""
function heideldiag(
    c::AbstractChains;
    alpha::Real = 0.05,
    eps::Real = 0.1,
    etype = :imse,
    args...,
)
    _, p, m = size(c.value)
    vals = Array{Float64}(undef, p, 6, m)
    for j = 1:p, k = 1:m
        vals[j, :, k] = heideldiag(
            c.value[:, j, k],
            alpha = alpha,
            eps = eps,
            etype = etype,
            start = c.range.start;
            args...,
        )
    end
    hdr =
        header(c) *
        "\nHeidelberger and Welch Diagnostic:\n" *
        "Target Halfwidth Ratio = $eps\nAlpha = $alpha\n"
    ChainSummary(
        vals,
        c.names,
        ["Burn-in", "Stationarity", "p-value", "Mean", "Halfwidth", "Test"],
        hdr,
    )
end

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
"""
Find all dipole in an activity map.

Determines the local maxima in a 3 dimensional array

#### Input

* s: Activity in 3d matrix
* window: Windowing to use in each dimension for min max filter
* x,y,z: Coordinates associated with s matrix

#### Output

* dips: An array of dipoles

"""
function find_dipoles(s::Array{T, 3}; window::Array{Int}=[6,6,6], x::AbstractVector{T}=1:size(s,1),
                      y::AbstractVector{T}=1:size(s,2), z::AbstractVector{T}=1:size(s,3)) where T <: Number

    @debug("Finding dipoles for 3d array")

    localmaxima_locations = findlocalmaxima(s)
    peak_values = s[localmaxima_locations]
    smallest_peak_value = (0.1 * maximum(peak_values))
    localmaxima_locations_subset = findlocalmaxima(s .>= smallest_peak_value)
    localmaxima_locations_subset = localmaxima_locations

    # Store dipoles in an array
    dips = Dipole[]

    for location in localmaxima_locations_subset
        push!(dips, Dipole("Unknown", location[1] * u"m", location[2] * u"m", location[3] * u"m", 0, 0, 0, 0, 0, s[location]))
    end

    left_side = s[x .< 0, :, :]
    a = CartesianIndices(size(left_side))[findall(left_side .== maximum(left_side))]
    a = a[1]
    xidx = a[1]
    yidx = a[2]
    zidx = a[3]
    push!(dips, Dipole("Unknown", x[xidx[1]] * u"m", y[yidx[1]] * u"m", z[zidx[1]] * u"m", 0, 0, 0, 0, 0, maximum(left_side)))

    right_side = s[x .> 0, :, :]
    a = CartesianIndices(size(right_side))[findall(right_side .== maximum(right_side))]
    a = a[1]
    xidx = a[1]
    yidx = a[2]
    zidx = a[3]
    x_tmp = x[x .> 0]
    push!(dips, Dipole("Unknown", x_tmp[xidx[1]] * u"m", y[yidx[1]] * u"m", z[zidx[1]] * u"m", 0, 0, 0, 0, 0, maximum(right_side)))

    # Sort dipoles by size
    vec(dips[sortperm([dip.size for dip in dips], rev=true)])
end

#######################################
#
# Find the best dipoles from selection
#
#######################################

"""
Find best dipole relative to reference location.

Finds the largest dipole within a specified distance of a reference location

#### Input

* ref: Reference coordinate or dipole
* dips: Dipoles to find the best dipole from
* maxdist: Maximum distance a dipole can be from the reference

#### Output

* dip: The best dipole

"""
function best_dipole(ref::Union{Coordinate, Dipole}, dips::Array{Dipole}; maxdist::Number=0.30, min_dipole_size::Real=-Inf, kwargs...)

    @info("Calculating best dipole for $(length(dips)) dipoles")

    dips = dips[findall([d.size > min_dipole_size for d in dips])]

    if length(dips) > 0

      # Find all dipoles within distance
      dists = [euclidean(ref, dip) for dip=dips]
      valid_dist = dists .< maxdist

      if sum(valid_dist) >= 2
          # Valid dipoles exist find the largest one
          sizes = [dip.size for dip =dips]
          bestdip = maximum(sizes[valid_dist])
          dip = dips[findall(sizes .== bestdip)]
          @debug("$(sum(valid_dist)) dipoles within $(maxdist) m. ")

      elseif sum(valid_dist) == 1
          # Return the one valid dipole
          dip = dips[findall(valid_dist)]
          @debug("Only one dipole within $(maxdist) m. ")

      else
          # No dipoles within distance
          # Take the closest
          bestdip = minimum(dists)
          dip = dips[findall(dists .== bestdip)]
          @debug("No dipole within $(maxdist) m. ")

      end
      @debug("Best = $(euclidean(ref, dip[1]))")

      return dip[1]

    else
      return NaN
    end
end


text:
# let statements allocate new variable bindings each time they run
# this is also an example of how to do closures
let state = 0
   global counter() = (state += 1)
end

println(counter())
println(counter())

# struct (composite type) example
struct Foo
     bar
     baz
 end

foo = Foo(1, 2)
println("bar: $(foo.bar), baz: $(foo.baz)")

# example of defining a new outer constructor method that only take 1 arg
Foo(x) = Foo(x, x)
foo = Foo(3)
println("bar: $(foo.bar), baz: $(foo.baz)")

# example of inner constructor method
struct OrderedPair
     x::Real
     y::Real
     OrderedPair(x,y) = x > y ? error("out of order") : new(x,y)
end

op1 = OrderedPair(1, 2)
println("(x, y) = ($(op1.x), $(op1.y))")

try
   op2 = OrderedPair(2, 1)
catch e
   println(e)
end

# conversion example
x = 2
println("typeof(x) = $(typeof(x))")
x = convert(Float64, x)
println("typeof(x) = $(typeof(x))")

# promotion - converting multiple different types to a common type
x = promote(1, 2.5, 3, 3//4)
println("promotion: $(x)")
x = promote(1 + 2im, 3//4)
println("promotion: $(x)")

# iteration - defined by iterate method
struct Squares
    count::Int
end
Base.iterate(S::Squares, state=1) = state > S.count ? nothing : (state*state, state+1)

for i in Squares(7)
    println(i)
end

println("is 25 a square number? $(25 in Squares(10))")

using Statistics
println("mean of the squares up to 100^2: $(mean(Squares(100)))")

max_stars_repo_name: suicidejack/julia-practice

max_issues_repo_name: suicidejack/julia-practice

max_forks_repo_name: suicidejack/julia-practice


text:
using Pkg
Pkg.activate(".")

#using LSDescent

using NLPModels, JuMP,  NLPModelsJuMP

using SolverCore, Logging

using LinearOperators

using LinearAlgebra

function test_algo(algo::Function, nlp, B, maxiter )
    reset!(nlp)

    logger = Logging.ConsoleLogger(stderr,Logging.Warn)
    Logging.with_logger(logger) do
        iter, f, normg, B, x = algo(nlp, B₀ = B, maxiter = maxiter)
        @show iter, f, normg
    end

end

function test_algo_stp(algo::Function, nlp, stp, B )
    reset!(nlp)

    logger = Logging.ConsoleLogger(stderr,Logging.Warn)
    Logging.with_logger(logger) do
        stp, B = algo(nlp, stp=stp,  B₀ = B)
        iter, f, normg =  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score
        @show iter, f, normg
    end

end

function test_wrapper_stp(algo::Function, nlp, stp; kwargs... )
    reset!(nlp)

    logger = Logging.ConsoleLogger(stderr,Logging.Warn)
    Logging.with_logger(logger) do
        stp, B = algo(nlp, stp=stp; kwargs...)
        iter, f, normg =  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score
        @show iter, f, normg
    end

end

# test all solvers with the well known Woods test function
using OptimizationProblems

#include("woods.jl")
#include("genrose.jl")
#nlp = MathOptNLPModel(PureJuMP.dixmaank(40), name="dixmaank")
#nlp = MathOptNLPModel(PureJuMP.dixmaang(100), name="dixmaang")
#nlp = MathOptNLPModel(PureJuMP.srosenbr(10), name="srosenbr")
nlp = MathOptNLPModel(PureJuMP.woods(8), name="woods")
#nlp = MathOptNLPModel(PureJuMP.genrose(100), name="genrose")
n = nlp.meta.nvar

include("bfgsSolver.jl")
include("FormuleN2.jl")

maxiter = 1600
scaling = false
#@info "Version avec matrice"
#bfgs(nlp, maxiter = maxiter)

include("Type.jl")

@info "Version encapsulée opérateur, formule O(n^3)"
include("FormuleNOp.jl")
B = InverseBFGSOperator(Float64, n; scaling = scaling);

test_algo(bfgs, nlp, B, maxiter )
#bfgs(nlp, B₀ = B, maxiter = maxiter)
#Bfull =  Matrix(B)

@info "Version encapsulée opérateur, formule O(n^2)"
include("FormuleN2Op.jl")
B = InverseBFGSOperator(Float64, n; scaling = scaling)

test_algo(bfgs, nlp, B, maxiter )
#bfgs(nlp, B₀ = B, maxiter = maxiter)
#Bfull2 =  Matrix(B)

@info "Version encapsulée opérateur, formule O(n^2), stp"
using Stopping
include("bfgsStop.jl")
stp = NLPStopping(nlp,
                  NLPAtX(nlp.meta.x0)  )
stp.meta.max_iter = maxiter

B = InverseBFGSOperator(Float64, n; scaling = scaling)

test_algo_stp(bfgs_Stop, nlp, stp, B )

@info "Version encapsulée opérateur, formule O(n^2), stp LS"
using OneDmin
include("bfgsStopLS.jl")
stp = NLPStopping(nlp,
                  NLPAtX(nlp.meta.x0)  )
stp.meta.max_iter = maxiter
reinit!(stp)

test_wrapper_stp(M_bfgs_StopLS, nlp, stp; scaling = scaling )

B = InverseBFGSOperator(Float64, n; scaling = scaling)
reinit!(stp)

test_algo_stp(bfgs_StopLS, nlp, stp, B )

@info "Version L-BFGS"
reinit!(stp)

test_wrapper_stp(L_bfgs_Stop, nlp, stp; scaling = scaling )

B = InverseLBFGSOperator(Float64, n, mem=160; scaling = scaling)
reinit!(stp)

test_algo_stp(bfgs_Stop, nlp, stp, B)

;
#bfgs(nlp, B₀ = B, maxiter = maxiter)
#BL = Matrix(B)

#norm(Bfull-BL)

max_stars_repo_name: vepiteski/LSDescent

max_issues_repo_name: vepiteski/LSDescent

max_forks_repo_name: vepiteski/LSDescent


text:
using TerminalLoggers: StickyMessages

@testset "Sticky messages without ANSI codes" begin
    buf = IOBuffer()
    # Without TTY, messages are just piped through
    stickies = StickyMessages(buf, ansi_codes=false)
    push!(stickies, :a=>"Msg\n")
    @test String(take!(buf)) == "Msg\n"
    push!(stickies, :a=>"Msg\n")
    @test String(take!(buf)) == "Msg\n"
    pop!(stickies, :a)
    @test String(take!(buf)) == ""
end

@testset "Sticky messages with ANSI codes" begin
    buf = IOBuffer()
    dsize = (20, 80) # Intentionally different from default of 25 rows
    # In TTY mode, we generate various escape codes.
    stickies = StickyMessages(IOContext(buf, :displaysize=>dsize), ansi_codes=true)
    push!(stickies, :a=>"Msg\n")
    @test String(take!(buf)) ==  #scroll    #csr    #pos   #msg #pos
                                "\e[20;1H\n\e[1;19r\e[20;1HMsg\e[19;1H"
    push!(stickies, :a=>"MsgMsg\n")
    @test String(take!(buf)) == #clear      #msgpos #msg #repos
                               "\e[20;1H\e[J\e[20;1HMsgMsg\e[19;1H"
    push!(stickies, :b=>"BBB\n")
    @test String(take!(buf)) ==
        #clear       #scroll   #csr    #pos    #msgs      #pos
        "\e[20;1H\e[J\e[19;1H\n\e[1;18r\e[19;1HMsgMsg\nBBB\e[18;1H"
    pop!(stickies, :a)
    @test String(take!(buf)) == #clear       #csr    #pos   #msg #pos
                                "\e[19;1H\e[J\e[1;19r\e[20;1HBBB\e[18;1H"
    pop!(stickies, :b)
    @test String(take!(buf)) == #clear       #csr    #pos
                                "\e[20;1H\e[J\e[1;20r\e[19;1H"
    pop!(stickies, :b)
    @test String(take!(buf)) == ""

    push!(stickies, :a=>"αβγ\n")
    @test String(take!(buf)) ==  #scroll    #csr    #pos   #msg #pos
                                "\e[20;1H\n\e[1;19r\e[20;1Hαβγ\e[19;1H"

    push!(stickies, :b=>"msg\n")
    take!(buf)

    # Remove all two messages
    empty!(stickies)
    @test String(take!(buf)) == #clear       #csr    #pos
                                "\e[19;1H\e[J\e[1;20r\e[18;1H"
end

@noinline func_barrier(f) = f()

@testset "Sticky messages with ANSI codes" begin
    buf = IOBuffer()
    dsize = (20, 80) # Intentionally different from default of 25 rows
    func_barrier() do
        # Hide stickies variable behind a function barrier to make sure(er)
        # that it can be GC'd.
        stickies = StickyMessages(IOContext(buf, :displaysize=>dsize), ansi_codes=true)
        push!(stickies, :a=>"a-msg\n")
        push!(stickies, :b=>"b-msg\n")
        take!(buf)
        nothing
    end
    # Hack to force StickyMessages finalizer
    GC.gc(true) # trigger finalizer
    for i=1:1000
        yield() # allow async cleanup
    end
    @test String(take!(buf)) == #clear       #csr    #pos
                                "\e[19;1H\e[J\e[1;20r\e[18;1H"
end


text:
println("こんにちは")
println("AtCoder")

max_stars_repo_name: sugitanishi/competitive-programming

max_issues_repo_name: sugitanishi/competitive-programming

max_forks_repo_name: sugitanishi/competitive-programming


text:
module SSHT

using ssht_jll

# File ssht_adjoint.h

# File ssht_core.h

# L is number of modes; 0 ≤ l < L
# verbosity: 0 ≤ v ≤ 5

"""
    SSHT.core_dh_inverse_sov!(f::AbstractArray{Complex{Float64},2},
                              flm::AbstractVector{Complex{Float64}},
                              L::Integer,
                              spin::Integer,
                              verbosity::Integer=0)

Evaluate spin-weighted spherical harmonic coefficients `flm` with
spin-weight `spin`. `L = lmax+1` is the number of modes in `flm`.

The arrays `flm` and `f` must have the following sizes:

    nphi = SSHT.sampling_dh_nphi(L)
    ntheta = SSHT.sampling_dh_ntheta(L)
    length(flm) == L^2
    size(f) == (nphi, ntheta)

See the [SSHT
reference](https://astro-informatics.github.io/ssht/c/html/ssht__core_8h.html)
for details.

See also: [`SSHT.core_dh_inverse_sov`](@ref),
[`SSHT.core_dh_inverse_sov_real!`](@ref),
[`SSHT.core_dh_forward_sov!`](@ref).
"""
function core_dh_inverse_sov!(f::AbstractArray{Complex{Float64},2}, flm::AbstractVector{Complex{Float64}}, L::Integer,
                              spin::Integer, verbosity::Integer=0)
    nphi = sampling_dh_nphi(L)
    ntheta = sampling_dh_ntheta(L)
    any(size(f) .< (nphi, ntheta)) && throw(DimensionMismatch())
    length(flm) < L^2 && throw(DimensionMismatch())
    ccall((:ssht_core_dh_inverse_sov, libssht), Cvoid, (Ptr{Complex{Float64}}, Ptr{Complex{Float64}}, Cint, Cint, Cint), f, flm, L,
          spin, verbosity)
    return f
end

"""
    f = SSHT.core_dh_inverse_sov(flm::AbstractVector{Complex{Float64}},
                                 L::Integer,
                                 spin::Integer,
                                 verbosity::Integer=0)
    f::Array{Complex{Float64},2}

Evaluate spin-weighted spherical harmonic coefficients `flm` with
spin-weight `spin`. `L = lmax+1` is the number of modes in `flm`.

The array `flm` must have the length `L^2`. The result `f` will have
the following size:

    nphi = SSHT.sampling_dh_nphi(L)
    ntheta = SSHT.sampling_dh_ntheta(L)
    size(f) == (nphi, ntheta)

See the [SSHT
reference](https://astro-informatics.github.io/ssht/c/html/ssht__core_8h.html)
for details.

See also: [`SSHT.core_dh_inverse_sov!`](@ref),
[`SSHT.core_dh_inverse_sov_real`](@ref),
[`SSHT.core_dh_forward_sov`](@ref).
"""
function core_dh_inverse_sov(flm::AbstractVector{Complex{Float64}}, L::Integer, spin::Integer, verbosity::Integer=0)
    nphi = sampling_dh_nphi(L)
    ntheta = sampling_dh_ntheta(L)
    f = Array{Complex{Float64}}(undef, nphi, ntheta)
    core_dh_inverse_sov!(f, flm, L, spin, verbosity)
    return f
end

"""
    SSHT.core_dh_inverse_sov_real!(f::AbstractArray{Float64,2},
                                   flm::AbstractVector{Complex{Float64}},
                                   L::Integer,
                                   spin::Integer,
                                   verbosity::Integer=0)

Evaluate spin-weighted spherical harmonic coefficients `flm` with
spin-weight `spin`. `L = lmax+1` is the number of modes in `flm`. The
result `f` is real-valued.

The arrays `flm` and `f` must have the following sizes:

    nphi = SSHT.sampling_dh_nphi(L)
    ntheta = SSHT.sampling_dh_ntheta(L)
    length(flm) == L^2
    size(f) == (nphi, ntheta)

See the [SSHT
reference](https://astro-informatics.github.io/ssht/c/html/ssht__core_8h.html)
for details.

See also: [`SSHT.core_dh_inverse_sov_real`](@ref),
[`SSHT.core_dh_inverse_sov!`](@ref),
[`SSHT.core_dh_forward_sov_real!`](@ref).
"""
function core_dh_inverse_sov_real!(f::AbstractArray{Float64,2}, flm::AbstractVector{Complex{Float64}}, L::Integer,
                                   verbosity::Integer=0)
    nphi = sampling_dh_nphi(L)
    ntheta = sampling_dh_ntheta(L)
    any(size(f) .< (nphi, ntheta)) && throw(DimensionMismatch())
    length(flm) < L^2 && throw(DimensionMismatch())
    ccall((:ssht_core_dh_inverse_sov_real, libssht), Cvoid, (Ptr{Complex{Float64}}, Ptr{Complex{Float64}}, Cint, Cint), f, flm, L,
          verbosity)
    return f
end

"""
    f = SSHT.core_dh_inverse_sov_real(flm::AbstractVector{Complex{Float64}},
                                      L::Integer,
                                      spin::Integer,
                                      verbosity::Integer=0)
    f::Array{Float64,2}

Evaluate spin-weighted spherical harmonic coefficients `flm` with
spin-weight `spin`. `L = lmax+1` is the number of modes in `flm`. The
result `f` is real-valued.

The array `flm` must have the length `L^2`. The result `f` will have
the following size:

    nphi = SSHT.sampling_dh_nphi(L)
    ntheta = SSHT.sampling_dh_ntheta(L)
    size(f) == (nphi, ntheta)

See the [SSHT
reference](https://astro-informatics.github.io/ssht/c/html/ssht__core_8h.html)
for details.

See also: [`SSHT.core_dh_inverse_sov_real!`](@ref),
[`SSHT.core_dh_inverse_sov`](@ref),
[`SSHT.core_dh_forward_sov_real`](@ref).
"""
function core_dh_inverse_sov_real(flm::AbstractVector{Complex{Float64}}, L::Integer, verbosity::Integer=0)
    nphi = sampling_dh_nphi(L)
    ntheta = sampling_dh_ntheta(L)
    f = Array{Float64}(undef, nphi, ntheta)
    core_dh_inverse_sov_real!(f, flm, L, verbosity)
    return f
end

"""
    SSHT.core_dh_forward_sov!(flm::AbstractVector{Complex{Float64}},
                              f::AbstractArray{Complex{Float64},2},
                              L::Integer,
                              spin::Integer,
                              verbosity::Integer=0)

Calculate spin-weighted spherical harmonic coefficients with
spin-weight `spin` `flm` from grid point values `f`. `L = lmax+1` is
the number of modes in `flm`.

The arrays `flm` and `f` must have the following sizes:

    nphi = SSHT.sampling_dh_nphi(L)
    ntheta = SSHT.sampling_dh_ntheta(L)
    length(flm) == L^2
    size(f) == (nphi, ntheta)

Use [`SSTH.sampling_elm2ind`](@ref) to access individual modes in the
`flm` array.

See the [SSHT
reference](https://astro-informatics.github.io/ssht/c/html/ssht__core_8h.html)
for details.

See also: [`SSHT.core_dh_forward_sov`](@ref),
[`SSHT.core_dh_forward_sov_real!`](@ref),
[`SSHT.core_dh_inverse_sov!`](@ref).
"""
function core_dh_forward_sov!(flm::AbstractVector{Complex{Float64}}, f::AbstractArray{Complex{Float64},2}, L::Integer,
                              spin::Integer, verbosity::Integer=0)
    nphi = sampling_dh_nphi(L)
    ntheta = sampling_dh_ntheta(L)
    length(flm) < L^2 && throw(DimensionMismatch())
    any(size(f) .< (nphi, ntheta)) && throw(DimensionMismatch())
    ccall((:ssht_core_dh_forward_sov, libssht), Cvoid, (Ptr{Complex{Float64}}, Ptr{Complex{Float64}}, Cint, Cint, Cint), flm, f, L,
          spin, verbosity)
    return flm
end

"""
    flm = SSHT.core_dh_forward_sov(fl:AbstractVector{Complex{Float64}},
                                   L::Integer,
                                   spin::Integer,
                                   verbosity::Integer=0)
    flm::Vector{Complex{Float64}}

Calculate spin-weighted spherical harmonic coefficients `flm` with
spin-weight `spin` from grid point values `f`. `L = lmax+1` is the
number of modes in `flm`.

The array `flm` will have length `L^2`. The input `f` must have the
following size:

    nphi = SSHT.sampling_dh_nphi(L)
    ntheta = SSHT.sampling_dh_ntheta(L)
    size(f) == (nphi, ntheta)

Use [`SSTH.sampling_elm2ind`](@ref) to access individual modes in the
`flm` array.

See the [SSHT
reference](https://astro-informatics.github.io/ssht/c/html/ssht__core_8h.html)
for details.

See also: [`SSHT.core_dh_forward_sov!`](@ref),
[`SSHT.core_dh_forward_sov_real`](@ref),
[`SSHT.core_dh_inverse_sov`](@ref).
"""
function core_dh_forward_sov(f::AbstractArray{Complex{Float64},2}, L::Integer, spin::Integer, verbosity::Integer=0)
    flm = Array{Complex{Float64}}(undef, L^2)
    core_dh_forward_sov!(flm, f, L, spin, verbosity)
    return flm
end

"""
    SSHT.core_dh_forward_sov_real!(flm::AbstractVector{Complex{Float64}},
                                   f::AbstractArray{Float64,2},
                                   L::Integer,
                                   spin::Integer,
                                   verbosity::Integer=0)

Calculate spin-weighted spherical harmonic coefficients `flm` with
spin-weight `spin` from real-valued grid point values `f`. `L =
lmax+1` is the number of modes in `flm`.

The arrays `flm` and `f` must have the following sizes:

    nphi = SSHT.sampling_dh_nphi(L)
    ntheta = SSHT.sampling_dh_ntheta(L)
    length(flm) == L^2
    size(f) == (nphi, ntheta)

Use [`SSTH.sampling_elm2ind`](@ref) to access individual modes in the
`flm` array.

See the [SSHT
reference](https://astro-informatics.github.io/ssht/c/html/ssht__core_8h.html)
for details.

See also: [`SSHT.core_dh_forward_sov_real`](@ref),
[`SSHT.core_dh_forward_sov!`](@ref),
[`SSHT.core_dh_inverse_sov_real!`](@ref).
"""
function core_dh_forward_sov_real!(flm::AbstractVector{Complex{Float64}}, f::AbstractArray{Float64,2}, L::Integer,
                                   verbosity::Integer=0)
    nphi = sampling_dh_nphi(L)
    ntheta = sampling_dh_ntheta(L)
    length(flm) < L^2 && throw(DimensionMismatch())
    any(size(f) .< (nphi, ntheta)) && throw(DimensionMismatch())
    ccall((:ssht_core_dh_forward_sov_real, libssht), Cvoid, (Ptr{Complex{Float64}}, Ptr{Complex{Float64}}, Cint, Cint), flm, f, L,
          verbosity)
    return flm
end

"""
    flm = SSHT.core_dh_inverse_sov_real(f::AbstractVector{Float64},
                                        L::Integer,
                                        spin::Integer,
                                        verbosity::Integer=0)
    f::Array{Complex{Float64},2}

Calculate spin-weighted spherical harmonic coefficients `flm` with
spin-weight `spin` from real-valued grid point values `f`. `L =
lmax+1` is the number of modes in `flm`.

The array `flm` must have the length `L^2`. The result `f` will have
the following size:

    nphi = SSHT.sampling_dh_nphi(L)
    ntheta = SSHT.sampling_dh_ntheta(L)
    size(f) == (nphi, ntheta)

Use [`SSTH.sampling_elm2ind`](@ref) to access individual modes in the
`flm` array.

See the [SSHT
reference](https://astro-informatics.github.io/ssht/c/html/ssht__core_8h.html)
for details.

See also: [`SSHT.core_dh_forward_sov_real!`](@ref),
[`SSHT.core_dh_forward_sov`](@ref),
[`SSHT.core_dh_inverse_sov_real`](@ref).
"""
function core_dh_forward_sov_real(f::AbstractArray{Float64,2}, L::Integer, verbosity::Integer=0)
    flm = Array{Complex{Float64}}(undef, L^2)
    core_dh_forward_sov_real!(flm, f, L, verbosity)
    return flm
end

# File ssht_dl.h

@enum dl_size_t begin
    DL_QUARTER = 0
    SSHT_DL_QUARTER_EXTENDED
    SSHT_DL_HALF
    SSHT_DL_FULL
end

@enum dl_method_t begin
    SSHT_DL_RISBO = 0
    SSHT_DL_TRAPANI
end

# File ssht_sampling.h

"""
    n = sampling_dh_n(L::Integer)
    n::Int

Total number of unique collocation points for `L = lmax+1` modes that.
(If the collocation scheme places points onto the poles, then not all
collocation points are unique.)

See also: [`SSHT.sampling_dh_nphi`](@ref),
[`SSHT.sampling_dh_ntheta`](@ref).
"""
sampling_dh_n(L::Integer) = Int(ccall((:ssht_sampling_dh_n, libssht), Cint, (Cint,), L))

"""
    nphi = sampling_dh_nphi(L::Integer)
    nphi::Int

Number of collocation points in the `phi` direction for `L = lmax+1`
modes.

See also: [`SSHT.sampling_dh_n`](@ref),
[`SSHT.sampling_dh_ntheta`](@ref).
"""
sampling_dh_nphi(L::Integer) = Int(ccall((:ssht_sampling_dh_nphi, libssht), Cint, (Cint,), L))

"""
    ntheta = sampling_dh_ntheta(L::Integer)
    ntheta::Int

Number of collocation points in the `theta` direction for `L = lmax+1`
modes.

See also: [`SSHT.sampling_dh_n`](@ref),
[`SSHT.sampling_dh_nphi`](@ref).
"""
sampling_dh_ntheta(L::Integer) = Int(ccall((:ssht_sampling_dh_ntheta, libssht), Cint, (Cint,), L))

"""
    phi = sampling_dh_p2phi(p::Integer, L::Integer)
    phi::Float64

Calculate the `phi` coordinate for point `p` in the `phi` direction
(`1 ≤ p ≤ nphi`).

See also: [`SSHT.sampling_dh_t2theta`](@ref),
[`SSHT.sampling_dh_nphi`](@ref).
"""
sampling_dh_p2phi(p::Integer, L::Integer) = Float64(ccall((:ssht_sampling_dh_p2phi, libssht), Cdouble, (Cint, Cint), p - 1, L))

"""
    theta = sampling_dh_p2theta(t::Integer, L::Integer)
    theta::Float64

Calculate the `theta` coordinate for point `t` in the `theta` direction
(`1 ≤ t ≤ ntheta`).

See also: [`SSHT.sampling_dh_p2phi`](@ref),
[`SSHT.sampling_dh_ntheta`](@ref).
"""
sampling_dh_t2theta(t::Integer, L::Integer) = Float64(ccall((:ssht_sampling_dh_t2theta, libssht), Cdouble, (Cint, Cint), t - 1, L))

"""
    w = sampling_weight_dh(theta_t::Real, L::Integer)
    w::Float64

Calculate the sampling weight of a point with the given `theta`
coordinate. This value is essentially `sin(theta) * dtheta`, possibly
modified to achive higher accuracy depending on how the collocation
points are spaced out in the `theta` direction.

Integrating a function over the sphere should be done as follows:

    # Input: Choose `L` and array `f`
    s = 0.0
    nphi = SSHT.sampling_dh_nphi(L)
    ntheta = SSHT.sampling_dh_ntheta(L)
    for p in 1:nphi, t in 1:ntheta
        phi = SSHT.sampling_dh_p2phi(p, L)
        theta = SSHT.sampling_dh_t2theta(t, L)
        # dtheta = π / ntheta
        dtheta = SSHT.sampling_weight_dh(theta, L)
        dphi = 2π / nphi
        s += f[p, t] * dtheta * dphi
    end
    # Output: s

See also: [`SSHT.sampling_dh_nphi`](@ref),
[`SSHT.sampling_dh_ntheta`](@ref), [`SSHT.sampling_dh_p2phi`](@ref),
[`SSHT.sampling_dh_t2theta`](@ref).
"""
function sampling_weight_dh(theta_t::Real, L::Integer)
    return Float64(ccall((:ssht_sampling_weight_dh, libssht), Cdouble, (Cdouble, Cint), theta_t, L))
end

sampling_mw_n(L::Integer) = Int(ccall((:ssht_sampling_mw_n, libssht), Cint, (Cint,), L))
sampling_mw_nphi(L::Integer) = Int(ccall((:ssht_sampling_mw_nphi, libssht), Cint, (Cint,), L))
sampling_mw_ntheta(L::Integer) = Int(ccall((:ssht_sampling_mw_ntheta, libssht), Cint, (Cint,), L))
sampling_mw_p2phi(p::Integer, L::Integer) = Float64(ccall((:ssht_sampling_mw_p2phi, libssht), Cdouble, (Cint, Cint), p - 1, L))
sampling_mw_t2theta(t::Integer, L::Integer) = Float64(ccall((:ssht_sampling_mw_t2theta, libssht), Cdouble, (Cint, Cint), t - 1, L))
function sampling_weight_mw(theta_t::Real, L::Integer)
    return Float64(ccall((:ssht_sampling_weight_mw, libssht), Cdouble, (Cdouble, Cint), theta_t, L))
end

"""
    ind = sampling_elm2ind(el::Integer, m::Integer)
    ind::Int

Calculate the mode array index `ind` for a given mode `l`, `m`. For `L
= lmax+1` modes, there are `L^2` modes in total.

See also: [`sampling_ind2elm`](@ref).
"""
function sampling_elm2ind(el::Int, m::Int)
    0 ≤ el || throw(DomainError(el, "Need 0 ≤ el"))
    -el ≤ m ≤ el || throw(DomainError(m, "Need -el ≤ m ≤ el"))
    return el^2 + el + m + 1
end
sampling_elm2ind(el::Integer, m::Integer) = sampling_elm2ind(Int(el), Int(m))

"""
    l, m = sampling_ind2elm(ind::Integer)
    l::Int
    m::Int

Calculate the mode numbers `l` and `m` from a given mode array index
`ind`. For `L = lmax+1` modes, there are `L^2` modes in total: `1 ≤
ind ≤ L^2`.

This function needs to evaluate a square root internally. If possible,
using [`sampling_elm2ind`](@ref) instead is slightly preferred.

See also: [`sampling_elm2ind`](@ref).
"""
function sampling_ind2elm(ind::Int)
    1 ≤ ind || throw(DomainError(ind, "Need 1 ≤ ind"))
    el = isqrt(ind - 1)
    m = ind - el^2 - el - 1
    return el, m
end
sampling_ind2elm(ind::Integer) = sampling_ind2elm(Int(ind))

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

"""
    eth!(ðflm::AbstractVector, flm::AbstractVector, L::Integer, spin::Integer)

Calculate the `ð` (eth) derivative of the spin-weighted spherical
harmonic coefficients `flm` with spin weight `s`. `L = lmax+1` is the
number of modes in `flm`. The result has spin weight `s+1`.

The arrays `ðflm` and `flm` must have length `L^2`.

See
[Wikipedia](https://en.wikipedia.org/wiki/Spin-weighted_spherical_harmonics)
for a definition of the `ð` (eth) operator.

See also: [`SSHT.eth`](@ref), [`SSHT.ethbar!`](@ref),
[`SSHT.core_dh_transform_sov!`](@ref),
[`SSHT.core_dh_inverse_sov!`](@ref).
"""
function eth!(ðflm::AbstractVector, flm::AbstractVector, L::Int, spin::Int)
    0 < L || throw(DomainError(L, "Need 0 < L"))
    T = typeof(abs(zero(eltype(flm))))
    # @assert T <: Real
    lmin = max(abs(spin), abs(spin + 1))
    lmax = L - 1
    for l in 0:min(lmin - 1, lmax), m in (-l):l
        # @assert (l - spin) * (l + spin + 1) ≤ 0
        ind = sampling_elm2ind(l, m)
        ðflm[ind] = 0
    end
    for l in lmin:lmax, m in (-l):l
        # @assert (l - spin) * (l + spin + 1) > 0
        ind = sampling_elm2ind(l, m)
        # <https://en.wikipedia.org/wiki/Spin-weighted_spherical_harmonics>
        # ð sYlm = + sqrt((l-s)*(l+s+1)) (s+1)Ylm
        ðflm[ind] = sqrt(T((l - spin) * (l + spin + 1))) * flm[ind]
    end
    return ðflm
end
eth!(ðflm::AbstractVector, flm::AbstractVector, L::Integer, spin::Integer) = eth!(ðflm, flm, Int(L), Int(spin))

"""
    ðflm = eth(flm::AbstractVector, L::Integer, spin::Integer)
    ðflm::AbstractVector

Calculate the `ð` (eth) derivative of the spin-weighted spherical
harmonic coefficients `flm` with spin weight `s`. `L = lmax+1` is the
number of modes in `flm`. The result has spin weight `s+1`.

The array `flm` must have length `L^2`. The result `ðflm` will also
have length `L^2`.

See
[Wikipedia](https://en.wikipedia.org/wiki/Spin-weighted_spherical_harmonics)
for a definition of the `ð` (eth) operator.

See also: [`SSHT.eth`](@ref), [`SSHT.ethbar!`](@ref),
[`SSHT.core_dh_transform_sov!`](@ref),
[`SSHT.core_dh_inverse_sov!`](@ref).
"""
eth(flm::AbstractVector, L::Integer, spin::Integer) = eth!(similar(flm), flm, L, spin)

"""
    ethbar!(ð̄flm::AbstractVector, flm::AbstractVector, L::Integer, spin::Integer)

Calculate the `ð̄` (eth-bar) derivative of the spin-weighted spherical
harmonic coefficients `flm` with spin weight `s`. `L = lmax+1` is the
number of modes in `flm`. The result has spin weight `s-1`.

The arrays `ð̄flm` and `flm` must have length `L^2`.

See
[Wikipedia](https://en.wikipedia.org/wiki/Spin-weighted_spherical_harmonics)
for a definition of the `ð̄` (eth-bar) operator.

See also: [`SSHT.ethbar`](@ref), [`SSHT.eth!`](@ref),
[`SSHT.core_dh_transform_sov!`](@ref),
[`SSHT.core_dh_inverse_sov!`](@ref).
"""
function ethbar!(ð̄flm::AbstractVector, flm::AbstractVector, L::Int, spin::Int)
    0 < L || throw(DomainError(L, "Need 0 < L"))
    T = typeof(abs(zero(eltype(flm))))
    # @assert T <: Real
    lmin = max(abs(spin), abs(spin - 1))
    lmax = L - 1
    for l in 0:min(lmin - 1, lmax), m in (-l):l
        # @assert (l + spin) * (l - spin + 1) ≤ 0
        ind = sampling_elm2ind(l, m)
        ð̄flm[ind] = 0
    end
    for l in lmin:lmax, m in (-l):l
        # @assert (l + spin) * (l - spin + 1) > 0
        ind = sampling_elm2ind(l, m)
        # <https://en.wikipedia.org/wiki/Spin-weighted_spherical_harmonics>
        # ð̄ sYlm = - sqrt((l+s)*(l-s+1)) (s-1)Ylm
        ð̄flm[ind] = -sqrt(T((l + spin) * (l - spin + 1))) * flm[ind]
    end
    return ð̄flm
end
ethbar!(ð̄flm::AbstractVector, flm::AbstractVector, L::Integer, spin::Integer) = ethbar!(ð̄flm, flm, Int(L), Int(spin))

"""
    ð̄flm = ethbar(flm::AbstractVector, L::Integer, spin::Integer)
    ð̄flm::AbstractVector

Calculate the `ð̄` (eth) derivative of the spin-weighted spherical
harmonic coefficients `flm` with spin weight `s`. `L = lmax+1` is the
number of modes in `flm`. The result has spin weight `s-1`.

The array `flm` must have length `L^2`. The result `ð̄flm` will also
have length `L^2`.

See
[Wikipedia](https://en.wikipedia.org/wiki/Spin-weighted_spherical_harmonics)
for a definition of the `ð` (eth) operator.

See also: [`SSHT.ethbar`](@ref), [`SSHT.eth!`](@ref),
[`SSHT.core_dh_transform_sov!`](@ref),
[`SSHT.core_dh_inverse_sov!`](@ref).
"""
ethbar(flm::AbstractVector, L::Integer, spin::Integer) = ethbar!(similar(flm), flm, L, spin)

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

export ash_grid_size, ash_nmodes
function ash_grid_size(lmax::Integer)
    0 ≤ lmax || throw(DomainError(lmax, "Need 0 ≤ lmax"))
    L = Int(lmax) + 1
    nphi = sampling_dh_nphi(L)
    ntheta = sampling_dh_ntheta(L)
    return (nphi, ntheta)::NTuple{2,Int}
end
function ash_nmodes(lmax::Integer)
    0 ≤ lmax || throw(DomainError(lmax, "Need 0 ≤ lmax"))
    L = Int(lmax) + 1
    return (L^2,)::NTuple{1,Int}
end

export ash_ntheta, ash_nphi, ash_thetas, ash_phis, ash_point_coord, ash_point_delta, ash_grid_as_phi_theta
ash_ntheta(lmax) = ash_grid_size(lmax)[2]::Int
ash_nphi(lmax) = ash_grid_size(lmax)[1]::Int
function ash_thetas(lmax::Integer)
    0 ≤ lmax || throw(DomainError(lmax, "Need 0 ≤ lmax"))
    L = Int(lmax) + 1
    ntheta = sampling_dh_ntheta(L)
    return [sampling_dh_t2theta(t, L) for t in 1:ntheta]
end
function ash_phis(lmax::Integer)
    0 ≤ lmax || throw(DomainError(lmax, "Need 0 ≤ lmax"))
    L = Int(lmax) + 1
    nphi = sampling_dh_nphi(L)
    return [sampling_dh_p2phi(p, L) for p in 1:nphi]
end
function ash_point_coord(ij::Union{CartesianIndex{2},NTuple{2,Int}}, lmax::Integer)
    0 ≤ lmax || throw(DomainError(lmax, "Need 0 ≤ lmax"))
    L = Int(lmax) + 1
    p, t = Tuple(ij)
    return sampling_dh_t2theta(t, L), sampling_dh_p2phi(p, L)
end
function ash_point_delta(ij::Union{CartesianIndex{2},NTuple{2,Int}}, lmax::Integer)
    0 ≤ lmax || throw(DomainError(lmax, "Need 0 ≤ lmax"))
    L = Int(lmax) + 1
    p, t = Tuple(ij)
    nphi = sampling_dh_nphi(L)
    theta = sampling_dh_t2theta(t, L)
    dtheta = sampling_weight_dh(theta, L) / sin(theta)
    dphi = 2π / nphi
    return dtheta, dphi
end
ash_grid_as_phi_theta(grid::AbstractMatrix) = grid

export ash_mode_index
function ash_mode_index(s::Integer, l::Integer, m::Integer, lmax::Integer)
    0 ≤ lmax || throw(DomainError(lmax, "Need 0 ≤ lmax"))
    abs(s) ≤ l ≤ lmax || throw(DomainError(l, "Need abs(s) ≤ l ≤ lmax"))
    -l ≤ m ≤ l || throw(DomainError(m, "Need -l ≤ m ≤ l"))
    return CartesianIndex(sampling_elm2ind(l, m))::CartesianIndex{1}
end
export ash_mode_numbers
function ash_mode_numbers(s::Int, ind::Union{CartesianIndex{1},NTuple{1,<:Int}}, lmax::Int)
    0 ≤ lmax || throw(DomainError(lmax, "Need 0 ≤ lmax"))
    l, m = sampling_ind2elm(ind[1])
    return (l, m)::NTuple{2,Int}
end
ash_mode_numbers(s::Integer, ind::NTuple{1,<:Integer}, lmax::Integer) = ash_mode_numbers(Int(s), NTuple{1,Int}(ind), Int(lmax))

export ash_transform!, ash_transform, ash_evaluate!, ash_evaluate
function ash_transform!(flm::AbstractArray{<:Complex}, f::AbstractMatrix{<:Complex}, s::Integer, lmax::Integer)
    L = lmax + 1
    return core_dh_forward_sov!(flm, f, L, s)
end
function ash_transform(f::AbstractMatrix{<:Complex}, s::Integer, lmax::Integer)
    L = lmax + 1
    return core_dh_forward_sov(f, L, s)
end

function ash_evaluate!(f::AbstractMatrix{<:Complex}, flm::AbstractArray{<:Complex}, s::Integer, lmax::Integer)
    L = lmax + 1
    return core_dh_inverse_sov!(f, flm, L, s)
end
function ash_evaluate(flm::AbstractArray{<:Complex}, s::Integer, lmax::Integer)
    L = lmax + 1
    return core_dh_inverse_sov(flm, L, s)
end

export ash_eth!, ash_eth, ash_ethbar!, ash_ethbar
function ash_eth!(ðflm::AbstractArray{<:Complex}, flm::AbstractArray{<:Complex}, s::Integer, lmax::Integer)
    L = lmax + 1
    return eth!(ðflm, flm, L, s)
end
function ash_eth(flm::AbstractArray{<:Complex}, s::Integer, lmax::Integer)
    L = lmax + 1
    return eth(flm, L, s)
end

function ash_ethbar!(ðflm::AbstractArray{<:Complex}, flm::AbstractArray{<:Complex}, s::Integer, lmax::Integer)
    L = lmax + 1
    return ethbar!(ðflm, flm, L, s)
end
function ash_ethbar(flm::AbstractArray{<:Complex}, s::Integer, lmax::Integer)
    L = lmax + 1
    return ethbar(flm, L, s)
end

end


text:
using MakieGallery, AbstractPlotting, GLMakie, Test

using MakieGallery: @block, @cell

database = MakieGallery.load_test_database()
tested_diff_path = joinpath(@__DIR__, "tested_different")
test_record_path = joinpath(@__DIR__, "test_recordings")
isdir(tested_diff_path) && rm(tested_diff_path, force = true, recursive = true)
mkpath(tested_diff_path)
isdir(test_record_path) && rm(test_record_path, force = true, recursive = true)
mkpath(test_record_path)
examples = MakieGallery.record_examples(test_record_path)
@test length(examples) == length(database)
MakieGallery.run_comparison(test_record_path, tested_diff_path)

empty!(database) # remove other examples
include("glmakie_tests.jl") # include GLMakie specific tests
# THese examples download additional data - don't want to deal with that!
for path in (tested_diff_path, test_record_path)
    rm(path, force = true, recursive = true)
    mkpath(path)
end

examples = MakieGallery.record_examples(test_record_path)
MakieGallery.run_comparison(test_record_path, tested_diff_path, maxdiff=0.00001)


text:
include("../common.jl")
build_julia(ARGS, v"1.7.0-beta3")

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
@testset "Euclidean graphs" begin
    N = 10
    d = 2
    g, weights, points = @inferred(euclidean_graph(N, d))
    @test nv(g) == N
    @test ne(g) == N * (N - 1) ÷ 2
    @test (d, N) == size(points)
    @test maximum(x -> x[2], weights) <= sqrt(d)
    @test minimum(x -> x[2], weights) >= 0
    @test maximum(points) <= 1
    @test minimum(points) >= 0.

    g, weights, points = @inferred(euclidean_graph(N, d, bc=:periodic))
    @test maximum(x -> x[2], weights) <= sqrt(d / 2)
    @test minimum(x -> x[2], weights) >= 0.
    @test maximum(points) <= 1
    @test minimum(points) >= 0.

    @test_throws DomainError euclidean_graph(points, L=0.01,  bc=:periodic)
    @test_throws ArgumentError euclidean_graph(points, bc=:badbc)

    # In our algorithm we ensure, that the resulting graph has the same
    # number of vertices as we have points. For this we have a special case,
    # where we have to insert vertices if the vertices with the highest indices are
    # isolated. This test ensures that this case is covered.
    @testset "Euclidean graph with vertex of highest index isolated" begin
        point1 = [1.0, 0.0, 0.0]
        point2 = [0.0, 1.0, 0.0]
        point3 = [0.0, 0.0, 100.0]
        matrix = hcat(point1, point2, point3)
        g, _ = euclidean_graph(matrix, cutoff= 5.0)
        @test has_edge(g, 1, 2) && ne(g) == 1
        @test nv(g) == 3
    end
end


text:
#= Mol Operations
Functions to process raw OxygenMols and extract tangential low level information from such.

=#

"""
    getbonds(mol)

Retrieves list of bonds in a molecule.
Bond list formatted as `(starting atom, ending atom, bond type)`.
All bonds are unique.

# Examples
```
julia> getbonds(smilestomol("c1ccccc1"))
6-element Array{Tuple,1}:
 (1, 2, 1.5)
 (1, 6, 1.5)
 (2, 3, 1.5)
 (3, 4, 1.5)
 (4, 5, 1.5)
 (5, 6, 1.5)
```
"""
function getbonds(mol::OxygenMol)::Array{Tuple,1}
    bondlist = []
    for (start, bonds) in enumerate(mol.adj)
        for (last, type) in bonds
            if last > start
                push!(bondlist, (start, last, type))
            end
        end
    end
    bondlist
end

"""
    getadjmatrix(mol)

Generates adjacency matrix of molecule, weighted by bond type

# Examples
```
julia> getadjmatrix(smilestomol("c1ccccc1"))
6×6 Array{Float32,2}:
 0.0  1.5  0.0  0.0  0.0  1.5
 1.5  0.0  1.5  0.0  0.0  0.0
 0.0  1.5  0.0  1.5  0.0  0.0
 0.0  0.0  1.5  0.0  1.5  0.0
 0.0  0.0  0.0  1.5  0.0  1.5
 1.5  0.0  0.0  0.0  1.5  0.0
```
"""
function getadjmatrix(::Type{T}, mol::OxygenMol)::Array{T,2} where {T<:AbstractFloat}
    adj = zeros(T, (length(mol.atoms), length(mol.atoms)))
    for (start, bonds) in enumerate(mol.adj)
        for (last, type) in bonds
            adj[start, last] = type
        end
    end
    adj
end

getadjmatrix(mol::OxygenMol) = getadjmatrix(Float32, mol::OxygenMol)

function _onehot(states::Array{T,1}, state::T)::Array{Float32,1} where {T}
    [Float32(state == s) for s in states]
end

"""
    atomtypes(mol)

Creates one-hot mapping of atoms in a molecule to common atoms.
Current valid atoms array is `["C", "N", "O", "S", "F", "P", "Cl", "Mg", "Na", "Br", "Ca", "Cu", "Pd", "I", "Al"]`.

# Examples
```
julia> atomtypes(smilestomol("c1ccccc1"))
6×15 Array{Float32,2}:
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
function atomtypes(::Type{T}, mol::OxygenMol)::Array{T,2} where {T<:AbstractFloat}
    validatoms = String["C", "N", "O", "S", "F", "P", "Cl", "Mg", "Na", "Br", "Ca", "Cu", "Pd", "I", "Al"]
    atoms = map(a -> _onehot(validatoms, a.symbol), mol.atoms)
    transpose(hcat(atoms...))
end

atomtypes(mol::OxygenMol) = atomtypes(Float32, mol::OxygenMol)

"""
    valencies(mol)

Creates one-hot mapping of atom valencies in a molecule.
Current maximum valid valencies (number of connected non-hydrogen atoms) is 6.

# Examples
```
julia> valencies(smilestomol("c1ccccc1"))
6×6 Array{Int64,2}:
 0  1  0  0  0  0
 0  1  0  0  0  0
 0  1  0  0  0  0
 0  1  0  0  0  0
 0  1  0  0  0  0
 0  1  0  0  0  0
```
"""
function valencies(mol::OxygenMol)::Array{Int,2}
    valcount = [length(bonds) for bonds in mol.adj]
    valences = map(a -> _onehot(collect(1:6), a), valcount)
    transpose(hcat(valences...))
end

"""
    hybridization(mol)

Creates one-hot mapping of atom hybridization states in a molecule.
Currently creates onehot mapping of bonding electron count. Hybridization states will be added soon.

# Examples
```
julia> hybridization(smilestomol("c1ccccc1"))
6×7 Array{Int64,2}:
 0  0  0  0  1  0  0
 0  0  0  0  1  0  0
 0  0  0  0  1  0  0
 0  0  0  0  1  0  0
 0  0  0  0  1  0  0
 0  0  0  0  1  0  0
```
"""
function hybridization(mol::OxygenMol)::Array{Int,2}
    #TODO: Retrieve hybridization states instead of bonding electrons.
    electrons(bonds) = sum([e for (_, e) in bonds])
    hybcount = [electrons(bonds) for bonds in mol.adj]
    hybrid = map(a -> _onehot(collect(1:.5:4), a), hybcount)
    transpose(hcat(hybrid...))
end

"""
    molfeatures(mol)

Creates atom-wise features matrix for a molecule based on atomtypes, valence, and hybridization.
Note: `molfeatures` and `getadjmatrix` do **not** return padded arrays.

# Examples
```
julia> molfeatures(smilestomol("c1ccccc1"))
6×28 Array{Float32,2}:
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
```
"""
function molfeatures(::Type{T}, mol::OxygenMol)::Array{T, 2} where {T <: AbstractFloat}
    atoms = atomtypes(T, mol)
    valence = convert(Array{T}, valencies(mol))
    hybrid = convert(Array{T}, hybridization(mol))
    hcat(atoms, valence, hybrid)
end

molfeatures(mol::OxygenMol) = molfeatures(Float32, mol::OxygenMol)


text:
using Test, MadNLP, MadNLPPardiso, MadNLPTests

testset = [
    [
        "Pardiso",
        ()->MadNLP.Optimizer(
            linear_solver=MadNLPPardiso,
            print_level=MadNLP.ERROR),
        [],
        @isdefined(MadNLPPardiso)
    ],
    [
        "PardisoMKL",
        ()->MadNLP.Optimizer(
            linear_solver=MadNLPPardisoMKL,
            print_level=MadNLP.ERROR),
        ["eigmina"]
    ]
]

@testset "MadNLPPardiso test" begin
    for (name,optimizer_constructor,exclude) in testset
        test_madnlp(name,optimizer_constructor,exclude)
    end
end


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

# Enum for escaped (multi-byte) keys such as the arrows or the home/end keys
@enum(Key,
    ARROW_LEFT = 1000,
    ARROW_RIGHT,
    ARROW_UP,
    ARROW_DOWN,
    DEL_KEY,
    HOME_KEY,
    END_KEY,
    PAGE_UP,
    PAGE_DOWN)

readbyte(stream::IO=stdin) = read(stream, Char)

# Read the next key from stdin. It is also able to read several bytes for
#   escaped keys such as the arrow keys, home/end keys, etc.
# Escaped keys are returned using the `Key` enum.
readkey(stream::Base.LibuvStream=stdin) = UInt32(_readkey(stream))
function _readkey(stream::Base.LibuvStream=stdin)
    c = readbyte(stream)

    # Escape characters
    if c == '\x1b'
        stream.buffer.size < 2 && return '\x1b'
        esc_a = readbyte(stream)
        esc_a == 'v' && return PAGE_UP  # M-v
        esc_a == '<' && return HOME_KEY # M-<
        esc_a == '>' && return END_KEY  # M->

        stream.buffer.size < 3 && return '\x1b'
        esc_b = readbyte(stream)

        if esc_a == '[' || esc_a == 'O'
            if esc_b >= '0' && esc_b <= '9'
                stream.buffer.size < 4 && return '\x1b'
                esc_c = readbyte(stream)
                if esc_c == '~'
                    esc_b == '1' && return HOME_KEY
                    esc_b == '4' && return END_KEY
                    esc_b == '3' && return DEL_KEY
                    esc_b == '5' && return PAGE_UP
                    esc_b == '6' && return PAGE_DOWN
                    esc_b == '7' && return HOME_KEY
                    esc_b == '8' && return END_KEY
                    return '\x1b'
                end
            else
                # Arrow keys
                esc_b == 'A' && return ARROW_UP
                esc_b == 'B' && return ARROW_DOWN
                esc_b == 'C' && return ARROW_RIGHT
                esc_b == 'D' && return ARROW_LEFT
                esc_b == 'H' && return HOME_KEY
                esc_b == 'F' && return END_KEY
                return '\x1b'
            end
        end
        esc_a == 'H' && return HOME_KEY
        esc_a == 'F' && return END_KEY
        return '\x1b'
    end
    c == '\x16' && return PAGE_DOWN  # C-v
    c == '\x10' && return ARROW_UP   # C-p
    c == '\x0e' && return ARROW_DOWN # C-n
    return c
end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using TimeWarp
 #using TimeWarp.WarpPlots
using TimeWarp.Datasets
using PyPlot # use pyplot
plt = PyPlot

# UCI data repository must be downloaded first, run TimeWarp.Datasets.download_ucr()
data, labels = ucr_traindata("Gun_Point");

# c = [ cl==1 ? :blue : :red for cl in class ]'
# plot(y, linecolor=c, legend=false)

avg, result = dba(data, FastDTW(1))

figure()
plt.plot(data,color="0.75")
plt.plot(avg,color="red")

if is_linux()
  plt.show()
end


text:
"""
    PriorBasis

An abstract basis used for sampling a prior Gaussian process, with weights.
"""
abstract type PriorBasis end

include("prior_basis/euclidean.jl")


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# simpleusercut.jl
#
# Solve a simple problem using user cuts. See documentation for more
# details about callbacks.
#############################################################################

using JuMP
using Gurobi

# We will use Gurobi, which requires that we manually set the attribute
# PreCrush to 1 if we have user cuts. We will also disable PreSolve, Cuts,
# and Heuristics so only our cut will be used
m = Model(solver=GurobiSolver(PreCrush=1, Cuts=0, Presolve=0, Heuristics=0.0))

# Define our variables to be inside a box, and integer
@variable(m, 0 <= x <= 2, Int)
@variable(m, 0 <= y <= 2, Int)

# Optimal solution is trying to go towards top-right corner (2.0, 2.0)
@objective(m, Max, x + 2y)

# We have one constraint that cuts off the top right corner
@constraint(m, y + x <= 3.5)

# Optimal solution of relaxed problem will be (1.5, 2.0)
# We can add a user cut that will cut of this fractional solution.

# We now define our callback function that takes one argument,
# the callback handle. Note that we can access m, x, and y because
# this function is defined inside the same scope
function mycutgenerator(cb)
    x_val = getvalue(x)
    y_val = getvalue(y)
    println("In callback function, x=$x_val, y=$y_val")

    # Allow for some impreciseness in the solution
    TOL = 1e-6

    # Check top right
    if y_val + x_val > 3 + TOL
        # Cut off this solution
        println("Fractional solution was in top right, cut it off")
        # Use the original variables
        @usercut(cb, y + x <= 3)
    end
end  # End of callback function

# Tell JuMP/Gurobi to use our callback function
addcutcallback(m, mycutgenerator)

# Solve the problem
solve(m)

# Print our final solution
println("Final solution: [ $(getvalue(x)), $(getvalue(y)) ]")


text:
function produce_epout!(epfields, epmat, alpha, scalefact, returnstatus, iter; drop_epfields)

    #= Scale back μ, s, av, va of epfields and lub, llb and Y =#
    scaleepfield!(scalefact, epfields)
    μ, σ = epfields.μ, epfields.s
    av, va = epfields.av, epfields.va
    sperm = isinf(alpha) ?
        sortperm(epmat.usperm) :
        eachindex(μ)
    μ, σ, av, va = μ[sperm], σ[sperm], av[sperm], va[sperm]

    sol = drop_epfields ? nothing : epfields
    return  EPOut(μ, σ, av, va, sol, returnstatus, iter)
end

produce_epout!(epmodel, returnstatus = UNSET_STATUS, iter = -1; drop_epfields = false) =
    produce_epout!(epmodel.epfields, epmodel.epmat, epmodel.alpha, epmodel.scalefact, returnstatus, iter; drop_epfields)


text:
function interactive_coupler_design(f₀ = 1.0)
    ## SETTINGS
    # Coupler
    n = 3
    K = "0.1"
    BW_REL = "0.3"
    POLYNOMIAL = :Butterworth

    # Plotting window
    _title = "Branch-Line Coupler Design"
    _res = (1920, 1080);
    _font = "Seaford";
    _fontsize = 20;
    _linewidth = 2;
    _color = :Hiroshige# TODO

    # Plots
    PTS = 1000
    fn_titles = ["Return losses", "Directivity", "Insertion Losses", "Coupling"]
    fn_labels = ["Γ", "D", "IL", "C"]
    f̄ = range(0.85, 1.15, length = PTS)
    f = f̄ .* f₀
    t = im .* tan.(π/4 .* f̄)

    ## SYNTHESISE FIRST COUPLER
    c = synthesise_coupler(n, K, BW_REL, POLYNOMIAL)
    funs = analysis_functions(c)

    ## MAKIE
    # Create figure
    set_theme!(font = _font,
        fontsize = _fontsize,
        linewidth = _linewidth,
        palette = ( colormap = _color,)); # color = _color
    fig = Figure(resolution = _res)

    # Title
    supertitle = Label(fig[1,1:10], _title, textsize = 40)

    # Grid of plots
    plot_grid = GridLayout()
    fig.layout[2,1:10] = plot_grid
    axs = [Axis(plot_grid[reverse(fldmod1(i,2))...],
                title = fn_titles[i],
                xlabel = "GHz",
                ylabel = "dB") for i in 1:4]
    # for (ax_, fn_, lb_) in zip(axs, funs, fn_labels)
    #     lines!(ax_, f, db.(fn_.(t)), label = lb_)
    #     axislegend(ax_)
    #     xlims!(ax_, f[1], f[end])
    # end

    # ## ADMITTANCE DISPLAY
    ab_display = GridLayout()
    fig.layout[3,1:10] = ab_display
    a_labels = [ Label(ab_display[1,i],
                    "",
                    textsize = 13,
                    tellwidth = false) for i in 1:10]
    b_labels = [ Label(ab_display[2,i],
                    "",
                    textsize = 13,
                    tellwidth = false) for i in 1:10]

    ## GENERATE FIRST PLOT
    _update_plot!(axs, a_labels, b_labels, c, fn_labels, t, f)

    ## PARAMETER INPUTS
    # n
    label_n = Label(fig[4,1], "n = ", halign = :right)
    input_n = Textbox(fig[4, 2],
                placeholder = "Enter number of branches",
                stored_string = "$(n)",
                validator = Int,
                tellwidth = false,
                halign = :left)
    # K
    label_k = Label(fig[5,1], "K = ", halign = :right)
    input_k = Textbox(fig[5,2],
                placeholder = "Enter a coupling factor",
                stored_string = K,
                validator = BigFloat,
                tellwidth = false,
                halign = :left)
    # BW_REL
    label_bwrel = Label(fig[6,1], "BW_REL = ", halign = :right)
    input_bwrel = Textbox(fig[6,2],
                    placeholder = "Enter a bandwidth factor",
                    stored_string = BW_REL,
                    validator = BigFloat,
                    tellwidth = false,
                    halign = :left)
    # POLYNOMIAL
    label_polynomial = Label(fig[7,1], "Polynomial = ", halign = :right)
    menu_polynomial = Menu(fig[7,2],
                            options = zip(["Butterworth", "Chebyshev"],
                                [:Butterworth,:Chebyshev]),
                            i_selected = 1,
                            selection = POLYNOMIAL,
                            halign = :left)

    ## EVENT-DRIVEN PLOTTING
    # n event
    on(input_n.stored_string) do new_n
        new_c = synthesise_coupler(parse(Int,to_value(new_n)),
                        to_value(input_k.stored_string),
                        to_value(input_bwrel.stored_string),
                        to_value(menu_polynomial.selection))
        _update_plot!(axs, a_labels, b_labels, new_c, fn_labels, t, f)
    end

    # K event
    on(input_k.stored_string) do new_k
        new_c = synthesise_coupler(parse(Int,to_value(input_n.stored_string)),
                        to_value(new_k),
                        to_value(input_bwrel.stored_string),
                        to_value(menu_polynomial.selection))
        _update_plot!(axs, a_labels, b_labels, new_c, fn_labels, t, f)
    end

    # BW event
    on(input_bwrel.stored_string) do new_bw
        new_c = synthesise_coupler(parse(Int,to_value(input_n.stored_string)),
                        to_value(input_k.stored_string),
                        to_value(new_bw),
                        to_value(menu_polynomial.selection))
        _update_plot!(axs, a_labels, b_labels, new_c, fn_labels, t, f)
    end

    # Polynomial event
    on(menu_polynomial.selection) do new_pol
        new_c = synthesise_coupler(parse(Int,to_value(input_n.stored_string)),
                        to_value(input_k.stored_string),
                        to_value(input_bwrel.stored_string),
                        to_value(new_pol))
        _update_plot!(axs, a_labels, b_labels, new_c, fn_labels, t, f)
    end

    return fig
end


text:
#=== write Simulation Parameters ===#
function wSPT(fid::IOStream)
    write(fid,"Output files:\n\n")
    write(fid,"SimuPara.bin: ") # SimuPara.bin
    write(fid,"simulation parameters file, read using function: ReadSimuPara  \n\n")
end


text:
inf(::Type{Float64}) = Inf64
inf(::Type{Float32}) = Inf32
inf(::Type{Float16}) = Inf16
# inf(::Type{T}) where {T <: ForwardDiff.Dual} = T(Inf)

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using LabelledArrays
using OrdinaryDiffEq
import OrdinaryDiffEq: solve
using ModelingToolkit
using Petri
import Petri: stripnullterms, fluxes, odefunc, quotesplat
using Test
using Plots

N(x) = sum(x)

@variables S,E,I,R, β,γ,μ
sir = Petri.Model([S,I,R],[(S+I, 2I), (I,R)])
seir = Petri.Model([S,E,I,R],[(S+I, E+I), (E,I), (I,R)])
seirs = Petri.Model([S,E,I,R],[(S+I, E+I), (E,I), (I,R), (R,S)])

@info "SIR Solving"
u0 = @LArray [100.0, 1, 0] (:S, :I, :R)
p = @LArray [0.35, 0.05] (:μ, :β)
fex = odefunc(sir, :sir)
f = eval(fex)
prob = ODEProblem(f,u0,(0.0,365.0),p)
sol = solve(prob,Tsit5())
plt = plot(sol, labels=LabelledArrays.symnames(typeof(sol[end]))|> collect)
savefig(plt, "img/sir_sol.png")
@test sol[end].S < 1
@test sol[end].I < 1
@test sol[end].R > 99

@info "SEIR Solving"
u0 = @LArray [100.0, 1, 0, 0] (:S, :E, :I, :R)
p = @LArray [0.35, 0.05, 0.05] (:μ, :β, :γ)
fex = odefunc(seir, :seir)
f = eval(fex)
prob = ODEProblem(f,u0,(0.0,365.0),p)
sol = solve(prob,Tsit5())
plt = plot(sol, labels=LabelledArrays.symnames(typeof(sol[end]))|> collect)
savefig(plt, "img/seir_sol.png")
@test sol[end].S < 1
@test sol[end].E < 1
@test sol[end].I < 1
@test sol[end].R > 99

@info "SEIRS Solving"
u0 = @LArray [100.0, 1, 0, 0] (:S, :E, :I, :R)
p = @LArray [0.35, 0.05, 0.07, 0.3] (:μ, :β, :γ, :η)
fex = odefunc(seirs, :seirs)
f = eval(fex)
prob = ODEProblem(f,u0,(0.0,365.0),p)
sol = solve(prob,Tsit5())
plt = plot(sol, labels=LabelledArrays.symnames(typeof(sol[end])) |> collect)
savefig(plt, "img/seirs_sol.png")
@test sol[end].S > 5
@test sol[end].E > 5
@test sol[end].I > 5
@test sol[end].R > 5
@test sol[end].S < 95
@test sol[end].E < 95
@test sol[end].I < 95
@test sol[end].R < 95


text:
using SimpleLockFiles
const SLF = SimpleLockFiles
using Test

@testset "SLF.jl" begin

    include("core_tests.jl")
    include("race_test.jl")

end


text:
origin(prim::HyperCube) = prim.origin
widths(prim::HyperCube{N,T}) where {N,T} = Vec{N,T}(prim.width)
maximum(prim::HyperCube{N,T}) where {N,T} = origin(prim)+widths(prim)
minimum(prim::HyperCube{N,T}) where {N,T} = origin(prim)


text:
include("./util.jl")

instructions = util.load_input(1)

floor = 0
basement = nothing

for (pos, i) in enumerate(instructions)
    if i == '('
        global floor += 1
    elseif i == ')'
        global floor -= 1
    else
        throw(DomainError("unknown instruction $i"))
    end

    if floor < 0 && isnothing(basement)
        global basement = pos
    end
end

println("Part 1: $floor")
println("Part 2: $basement")


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using RoboDojo

# ## Initial conditions
q1 = nominal_configuration(point_foot_quadruped)
height = 0.5
q1[3] += height
q1[6 + 3] += height
q1[9 + 3] += height
q1[12 + 3] += height
q1[15 + 3] += height

v1 = zeros(point_foot_quadruped.nq)
# v1= randn(point_foot_quadruped.nq)
# ## Time
h = 0.01
T = 100

# ## Simulator
s = Simulator(point_foot_quadruped, T,
    h=h)

# ## Simulate
simulate!(s, q1, v1)

using BenchmarkTools
@benchmark simulate!(s, q1, v1) setup=(s=s, q1=q1, v1=v1)

# ## Visualizer
vis = Visualizer()
open(vis)

# ## Visualize
visualize!(vis, s,
    fixed_camera=false)


text:
using Test, ColorSchemes, FileIO, Colors

function run_all_tests()
    # create a colorscheme from image file

    hokusai_test = extract(dirname(@__FILE__) * "/hokusai.jpg")

    @test length(hokusai_test) == 10

    # extract colors and weights from image file located here in the test directory
    c, w =  extract_weighted_colors(dirname(@__FILE__) * "/hokusai.jpg", 10, 10, 0.01; shrink = 4)

    @test length(c) == 10
    @test length(w) == 10

    # load scheme
    hok = ColorSchemes.hokusai

    @test length(hok) == 32

    # test that sampling schemes yield different values
    @test get(hokusai_test, 0.0) != get(hokusai_test, 0.5)

    # test sort
    @test sortcolorscheme(hokusai_test, rev=true) != sortcolorscheme(hokusai_test)

    # create weighted palette; there is some unpredictability here...
    csw = colorscheme_weighted(c, w, 37)
    @test 36 <= length(csw) <= 38

    # default is 50
    csw = colorscheme_weighted(c, w)
    @test length(csw) == 50

    # save as text
    colorscheme_to_text(ColorSchemes.hokusai, "hokusai_test_version", "hokusai_as_text.jl", comment="a test")

    @test filesize("hokusai_as_text.jl") > 2000 # TODO read it back in, for goodness sake...! :)

    open("hokusai_as_text.jl") do f
        lines = readlines(f)
        @test startswith(lines[4], "RGB{Float64}(0.085")
    end

    @testset "getinverse tests" begin
        cs = range(RGB(0, 0, 0), stop=RGB(1, 1, 1), length=5)
        @test getinverse(cs, cs[3]) == 0.5

        # Note that getinverse() takes the first closest match.
        cs = [RGB(0,0,0), RGB(1,1,1),
              RGB(0,0,0),
              RGB(0,0,0), RGB(1,1,1)];
        @test getinverse(cs, cs[3]) == 0

        # Note that getinverse() takes the left index when two are identical.
        cs = [RGB(1,1,1), RGB(0,0,0), RGB(0,0,0), RGB(0,1,1), RGB(1,1,0)];
        @test getinverse(cs, cs[2]) == 0.25
        cs = [RGB(0,0,0), RGB(0,0,0)];
        @test getinverse(cs, cs[1]) == 0

        cs = [RGB(0,0,0)];
        @test_throws MethodError getinverse(cs, cs[1])
        # (The above line throws for the same reason the below line does.
        #  If this behavior ever changes, so should `getinverse`.)
        @test_throws InexactError get(cs, 1.0, (1, 1))
    end

    @testset "convert_to_scheme tests" begin
        # Add color to a grayscale image.
        red_cs = range(RGB(0,0,0), stop=RGB(1,0,0), length=11)
        gray_img = range(RGB(0,0,0), stop=RGB(1,1,0), length=11)
        vs = [getinverse(gray_img, p) for p in red_cs]
        cs = [RGB(v,v,v) for v in vs]
        rcs = [get(red_cs, p) for p in vs]
        hcat(promote(red_cs,gray_img,cs,rcs)...)'
        new_img = convert_to_scheme(red_cs, gray_img)
        @test all(.≈(new_img, red_cs, atol=0.5))  # This is broken.. It should be way more specific. See next test.

        # Should be able to uniquely match each increasing color with the next
        # increasing color in the new scale.
        red_cs = range(RGB(0,0,0), stop=RGB(1,1,1))
        blue_scale_img = range(RGB(0,0,0), stop=RGB(0,0,1))
        new_img = convert_to_scheme(red_cs, blue_scale_img)
        @test_broken unique(new_img) == new_img
    end
end

function run_minimum_tests()

    # load scheme
    hok = ColorSchemes.hokusai

    @test length(hok) == 32

    # test sort
    @test sortcolorscheme(hok, rev=true) != sortcolorscheme(hok)

    # save as text
    colorscheme_to_text(hok, "hokusai_test_version", "hokusai_as_text.jl", comment="a test")

    @test filesize("hokusai_as_text.jl") > 2000

    open("hokusai_as_text.jl") do f
        lines = readlines(f)
        @test startswith(lines[4], "RGB{Float64}(0.085")
    end

    # convert an Array{T,2} to an RGB image
    tmp = get(ColorSchemes.leonardo, rand(10, 10))
    @test typeof(tmp) == Array{ColorTypes.RGB{Float64}, 2}

    # test conversion with default clamp
    x = [0.0 1.0 ; -1.0 2.0]
    y=get(ColorSchemes.leonardo, x)
    @test y[1,1] == y[2,1]
    @test y[1,2] == y[2,2]

    # test conversion with symbol clamp
    y2=get(ColorSchemes.leonardo, x, :clamp)
    @test y2 == y

    # test conversion with symbol extrema
    y2=get(ColorSchemes.leonardo, x, :extrema)
    @test y2[2,1] == y[1,1]   # Minimum now becomes one edge of ColorScheme
    @test y2[2,2] == y[1,2]   # Maximum now becomes other edge of ColorScheme
    @test y2[1,1] !== y2[2,1] # Inbetween values or now different

    # test conversion with manually supplied range
    y3=get(ColorSchemes.leonardo, x, (-1.0, 2.0))
    @test y3 == y2

    # test with steplen (#17)
    r  = range(0, stop=5, length=10)
    y  = get(ColorSchemes.leonardo, r)
    y2 = get(ColorSchemes.leonardo, collect(r))
    @test y == y2

    # test for specific value
    val = 0.2
    y   = get(ColorSchemes.leonardo, [val])
    y2  = get(ColorSchemes.leonardo, val)
    @test y2 == y[1]
end

if get(ENV, "COLORSCHEMES_KEEP_TEST_RESULTS", false) == "true"
        cd(mktempdir())
        @info("running tests in: $(pwd())")
        @info("...Keeping the results")
        run_minimum_tests()
        @info("Test images saved in: $(pwd())")
else
    mktempdir() do tmpdir
        cd(tmpdir) do
            @info("running tests in: $(pwd())")
            @info("but not keeping the results")
            @info("because you didn't do: ENV[\"COLORSCHEMES_KEEP_TEST_RESULTS\"] = \"true\"")
            run_minimum_tests()
            @info("Test images weren't saved. To see the test images, next time do this before running:")
            @info(" ENV[\"COLORSCHEMES_KEEP_TEST_RESULTS\"] = \"true\"")
            @info("These are the very brief tests.")
            @info("For more extensive testing, call `run_all_tests()`")
        end
    end
end


text:
using Zygote, LinearAlgebra
using ArnoldiMethod: partialschur, SR

function test_eigmin(A::AbstractMatrix)
    decomp, history = partialschur(A; nev = 1, which = SR())
    ev = minimum(real.(decomp.eigenvalues))
    return ev
end

A = rand(4, 4); A = A + A';

@show eigmin(A)
@show test_eigmin(A)
@show Zygote.gradient(test_eigmin, A)

max_stars_repo_name: ericphanson/Attempt-diff-eigmin

max_issues_repo_name: ericphanson/Attempt-diff-eigmin

max_forks_repo_name: ericphanson/Attempt-diff-eigmin


text:
# descriptor

mutable struct TensorDesc
    ptr::cudnnTensorDescriptor_t
end

unsafe_free!(td::TensorDesc) = cudnnDestroyTensorDescriptor(td.ptr)

Base.unsafe_convert(::Type{cudnnTensorDescriptor_t}, td::TensorDesc) = td.ptr

function TensorDesc(T::Type, size::NTuple{N,Integer}, strides::NTuple{N,Integer} = tuple_strides(size)) where N
    sz = Cint.(size) |> reverse |> collect
    st = Cint.(strides) |> reverse |> collect
    d = Ref{cudnnTensorDescriptor_t}()
    cudnnCreateTensorDescriptor(d)
    cudnnSetTensorNdDescriptor(d[], cudnnDataType(T), length(sz), sz, st)
    this = TensorDesc(d[])
    finalizer(unsafe_free!, this)
    return this
end

TensorDesc(a::CuArray) = TensorDesc(eltype(a), size(a), strides(a))

# wrappers

function cudnnAddTensor(A::CuArray{T,N}, C::CuArray{T,N}; alpha=1,
                        beta=1) where {T,N}
    aDesc = TensorDesc(A)
    cDesc = TensorDesc(C)
    cudnnAddTensor(handle(), Ref(T(alpha)), aDesc, A, Ref(T(beta)), cDesc, C)
    return C
end


text:
@testset "Deserialize" begin
    @testset "JSON" begin
        @testset "Nothing" begin
            @test deserialize("null", Genser.GenserJSON) === nothing
            @test deserialize(Nothing, "null", Genser.GenserJSON) === nothing
            @test deserialize(Union{Nothing,Int64}, "null", Genser.GenserJSON) === nothing
            @test_throws MethodError deserialize(Int64, "null", Genser.GenserJSON)
        end

        @testset "Simple values" begin
            @test deserialize("1", Genser.GenserJSON) === Int64(1)
            @test deserialize(Int32, "1", Genser.GenserJSON) === Int32(1)
            @test deserialize(Float16, "1", Genser.GenserJSON) === Float16(1)
            @test deserialize("1.5", Genser.GenserJSON) === Float64(1.5)
            @test deserialize(Float32, "1.5", Genser.GenserJSON) === Float32(1.5)
            @test deserialize("true", Genser.GenserJSON) === true
            @test deserialize("\"OKAY!\"", Genser.GenserJSON) === "OKAY!"
            @test deserialize(Char, "\"A\"", Genser.GenserJSON) === 'A'
        end

        @testset "Binary" begin
            @test deserialize(Vector{UInt8}, "[0,1,2]", Genser.GenserJSON) isa Vector{UInt8}
            @test deserialize(Vector{UInt8}, "[2,5,7]", Genser.GenserJSON) == [2,5,7]
        end

        @testset "UUID" begin
            id = Base.UUID("266945fa-8815-402a-ab1a-7701c5f53cc3")
            @test deserialize(Base.UUID, "\"266945fa-8815-402a-ab1a-7701c5f53cc3\"", Genser.GenserJSON) === id
        end

        @testset "Symbol" begin
            @test deserialize(Symbol, "\"tester\"", Genser.GenserJSON) === :tester
        end

        @testset "Arrays" begin
            @test deserialize("[1,2,3]", Genser.GenserJSON) == [1,2,3]
            @test deserialize(Vector{Int64}, "[1,2,3]", Genser.GenserJSON) == [1,2,3]
            @test eltype(deserialize(Vector{Int32}, "[1,2,3]", Genser.GenserJSON)) === Int32
            @test deserialize(Vector{Int32}, "[1,2,3]", Genser.GenserJSON) == [1,2,3]
        end

        @testset "Sets" begin
            @test deserialize(Set{Int64}, "[1,2,3]", Genser.GenserJSON) == Set([1,2,3])
            @test eltype(deserialize(Set{Int32}, "[1,2,3]", Genser.GenserJSON)) === Int32
            @test deserialize(Set{Int32}, "[1,2,3]", Genser.GenserJSON) == Set([1,2,3])
        end

        @testset "Dicts" begin
            v1 = Dict("a" => 1, "b" => 2, "c" => 3)
            v2 = Dict("a" => 1, "b" => true, "c" => "Apple")

            @test deserialize("{\"a\":1,\"b\":2,\"c\":3}", Genser.GenserJSON) == v1
            @test deserialize(typeof(v1), "{\"a\":1,\"b\":2,\"c\":3}", Genser.GenserJSON) == v1
            @test deserialize(Dict{Symbol,UInt8}, "{\"a\":1,\"b\":2,\"c\":3}", Genser.GenserJSON) == Dict(:a => 1, :b => 2, :c => 3)
            @test keytype(deserialize(Dict{Symbol,UInt8}, "{\"a\":1,\"b\":2,\"c\":3}", Genser.GenserJSON)) === Symbol
            @test valtype(deserialize(Dict{Symbol,UInt8}, "{\"a\":1,\"b\":2,\"c\":3}", Genser.GenserJSON)) === UInt8
            @test deserialize("{\"a\":1,\"b\":true,\"c\":\"Apple\"}", Genser.GenserJSON) == v2
            @test deserialize(typeof(v2), "{\"a\":1,\"b\":true,\"c\":\"Apple\"}", Genser.GenserJSON) == v2
        end

        @testset "Records" begin
            NamedJSONTuple = @NamedTuple{a::Int64,b::Bool,c::String}
            @test deserialize(NamedJSONTuple, "{\"a\":1,\"b\":true,\"c\":\"Apple\"}", Genser.GenserJSON) == (a=1,b=true,c="Apple")
            struct JSONDeserializeStruct
                a::Int64
                b::Bool
                c::String
            end
            @test deserialize(JSONDeserializeStruct, "{\"a\":1,\"b\":true,\"c\":\"Apple\"}", Genser.GenserJSON) == JSONDeserializeStruct(1,true,"Apple")

            @testset "Binary encoding" begin
                struct Base64DecodeStruct
                    data::Vector{UInt8}
                end

                Genser.gensertypefor(::Type{Base64DecodeStruct}) = GenserRecord{@NamedTuple{data::GenserBinaryValue{Encoding{:base64}}}}
                @test deserialize(Base64DecodeStruct, "{\"data\":\"YmFuYW5hcw==\"}", Genser.GenserJSON).data == Base64DecodeStruct(Vector{UInt8}("bananas")).data
            end
        end
    end
end