text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

using JuLIP, ACE, LinearAlgebra
using JuLIP.MLIPs: combine

#--- radial basis test

# Pr = basis.pibasis.basis1p.J
# ACE.Export.export_ace(@__DIR__() * "/testpot_rbasis.ace", Pr, ntests=5)

function export_dimer_test(V, fname, species=:Al)
   # export tests with the dimer structure
   # n_atoms = 2
   # # type x y z
   # 0 0.0 0.0 -1.0
   # 0 1.0 2.0 3.0
   at = Atoms(species, [ JVecF(0.0, 0.0, -1.0), JVecF(1.0, 2.0, 3.0) ])
   set_pbc!(at, false)
   fptr = open(@__DIR__() * fname * "_dimer_test.dat", "w")
   println(fptr, "E=$(energy(V, at))")
   println(fptr, "natoms = 2")
   println(fptr, "# type x y z")
   for n = 1:2
      r = at.X[n]
      println(fptr, "0 $(r[1]) $(r[2]) $(r[3])")
   end
   close(fptr)
end

function export_trimer_test(V, fname, species=:Al)
   # export tests with the trimer structuress
   at = Atoms(species, [ JVecF(0.0, 0.0, -1.0),
                         JVecF(1.0, 2.0, 2.0),
                         JVecF(0.0, 1.0, 1.0)  ])
   set_pbc!(at, false)
   fptr = open(@__DIR__() * fname * "_trimer_test.dat", "w")
   println(fptr, "E=$(energy(V, at))")
   println(fptr, "natoms = 3")
   println(fptr, "# type x y z")
   for n = 1:3
      r = at.X[n]
      println(fptr, "0 $(r[1]) $(r[2]) $(r[3])")
   end
   close(fptr)
end

#--- first test
V = nothing
N = 1
for maxdeg = 8:8
   basis = ACE.Utils.rpi_basis(; species = :Al, N = N, maxdeg=maxdeg)
   global V = ACE.Random.randcombine(basis)
   # V.coeffs[1][3] = 0.0
   fname = "/testpot_ord=$(N)_maxn=$(maxdeg)"
   ACE.Export.export_ace(@__DIR__() * fname * ".acejl", V)
   JuLIP.save_dict(@__DIR__() * fname * ".json", write_dict(V))
   # ACE.Export.export_ace_tests(@__DIR__() * "/testpot_$(N)_test", V, 3)
   export_dimer_test(V, fname)
end

#--- copy the first test over to the ace evaluator

for maxdeg = 8:8
   fname = "testpot_ord=$(N)_maxn=$(maxdeg)"
   filelist = [ fname * ".acejl",
                fname * ".json",
                fname * "_dimer_test.dat", ]
                # fname * _test_1.dat",
                # "testpot_$(N)_test_2.dat",
                # "testpot_$(N)_test_3.dat",

                # "testpot_$(N)_rbasis.ace" ]
   for f in filelist
      try
         run(`mv ./scripts/$f /Users/ortner/gits/ace-evaluator/test/julia/`)
      catch
      end
   end
end

#--- export a 3-body potential with a single 3-b basis function
#    and

ord = 2
rcut = 6.5
maxdeg = 4
basis = ACE.Utils.rpi_basis(; species = :Al, N = ord, maxdeg=maxdeg,
                               D = SparsePSHDegree(wL = 1.0), rcut=rcut)
V = ACE.Random.randcombine(basis)
inner = V.pibasis.inner[1]
for b in keys(inner.b2iAA)
   # if ACE.order(b) == 2
   if all([b1.l == 0 for b1 in b.oneps])
      V.coeffs[1][inner.b2iAA[b]] = 0.0
   end
   # end
end

fname = "/testpot_ord=$(ord)_mini"
ACE.Export.export_ace(@__DIR__() * fname * ".acejl", V)
JuLIP.save_dict(@__DIR__() * fname * ".json", write_dict(V))
ACE.Export.export_ace_tests(@__DIR__() * fname * "_test", V, 1, nrepeat=2)
export_dimer_test(V, fname, :Al)
export_trimer_test(V, fname, :Al)

filelist = [ fname * ".acejl",
             fname * ".json",
             fname * "_dimer_test.dat",
             fname * "_trimer_test.dat",
             fname * "_test_1.dat", ]
for f in filelist
   try
      run(`mv ./scripts/$f /Users/ortner/gits/ace-evaluator/test/julia/`)
   catch
   end
end

#--- export a random many-body potential

for N in [2, 3, 5], maxdeg in [10, ]

   rcut = 6.5

   basis = ACE.Utils.rpi_basis(; species = :Al, N = N, maxdeg=maxdeg, rcut=rcut)
   V = ACE.Random.randcombine(basis)
   fname = "/testpot_ord=$(N)_maxn=$(maxdeg)"
   ACE.Export.export_ace(@__DIR__() * fname * ".acejl", V)
   JuLIP.save_dict(@__DIR__() * fname * ".json", write_dict(V))
   ACE.Export.export_ace_tests(@__DIR__() * fname * "_test", V, 1)
   export_dimer_test(V, fname, :Al)
   export_trimer_test(V, fname, :Al)

   filelist = [ fname * ".acejl",
                fname * ".json",
                fname * "_dimer_test.dat",
                fname * "_trimer_test.dat",
                fname * "_test_1.dat"]
   for f in filelist
      try
         run(`mv ./scripts/$f /Users/ortner/gits/ace-evaluator/test/julia/`)
      catch
      end
   end
end

#---
# a potential with pair potential and repulsive core

# TODO:
#  - add on also E0
#  - move pairpot rcut to 7.0

basis = ACE.Utils.rpi_basis(; species = :Al, N = 4, maxdeg=8, pin = 2)
VN = randcombine(basis; diff=2)
pairbasis = ACE.Utils.pair_basis(; species = :Al, maxdeg = 8, rcut = 5.0)
V2 = combine(pairbasis, rand(8) .* (1:8).^(-2))
V2rep = ACE.PairPotentials.RepulsiveCore(V2, 2.0, -0.1234)
V1 = OneBody(:Al => -1.5 + rand())
V = JuLIP.MLIPs.SumIP(V1, V2rep, VN)

fname = "/ace_reppair"
JuLIP.save_dict(@__DIR__() * fname * ".json", write_dict(V))
# Vl = read_dict(load_dict(@__DIR__() * fname * ".json"))
# V1, V2rep, VN = tuple(Vl.components...)
ACE.Export.export_ace(@__DIR__() * fname * ".acejl", VN, V2rep, V1)
ACE.Export.export_ace_tests(@__DIR__() * fname * "_test", V, 1, s=:Al)
export_dimer_test(V, fname, :Al)
export_trimer_test(V, fname, :Al)

filelist = [ fname * ".acejl",
             fname * ".json",
             fname * "_dimer_test.dat",
             fname * "_trimer_test.dat",
             fname * "_test_1.dat"]
for f in filelist
   try
      run(`mv ./scripts/$f /Users/ortner/gits/ace-evaluator/test/julia/`)
   catch
   end
end

#---
# Cas Si potential

fname = "/Si_B6_N7_18_lap_1.1_rep"
D = load_dict(@__DIR__() * fname * ".json")
Vsi = read_dict(D["IP"])
ACE.Export.export_ace(@__DIR__() * fname * ".acejl",
                      Vsi.components[3], Vsi.components[2], Vsi.components[1])
ACE.Export.export_ace_tests(@__DIR__() * fname * "_test", Vsi, 1, s=:Si)
export_dimer_test(Vsi, fname, :Si)
export_trimer_test(Vsi, fname, :Si)

filelist = [ fname * ".acejl",
             fname * ".json",
             fname * "_dimer_test.dat",
             fname * "_trimer_test.dat",
             fname * "_test_1.dat"]
for f in filelist
   try
      run(`mv ./scripts/$f /Users/ortner/gits/ace-evaluator/test/julia/`)
   catch
   end
end

#---

fname = "/Si_B6_N7_18_lap_1"
VN = Vsi.components[3]
V2 = Vsi.components[2]
V1 = Vsi.components[1]
ACE.Export.export_ace(@__DIR__() * fname * ".acejl", VN)
ACE.Export.export_ace_tests(@__DIR__() * fname * "_test", VN, 1, s=:Si)
export_dimer_test(VN, fname, :Si)
export_trimer_test(VN, fname, :Si)
filelist = [ fname * ".acejl",
             fname * ".json",
             fname * "_dimer_test.dat",
             fname * "_trimer_test.dat",
             fname * "_test_1.dat"]
for f in filelist
   try
      run(`mv ./scripts$f /Users/ortner/gits/ace-evaluator/test/julia/`)
   catch
   end
end

#---
#
# at = Atoms(:Si, [ JVecF(0.0, 0.0, -1.0),
#                       JVecF(1.0, 2.0, 2.0),
#                       JVecF(0.0, 1.0, 1.0)  ])
# set_pbc!(at, false)
#
# energy(VN, at)
#
# site_energies(VN, at)
#
# X = at.X
# [ norm(X[i] - X[j]) for i = 1:2 for j = i+1:3 ]
#
# VN
#
#
# r = 1.732051
# z = AtomicNumber(:Si)
#
# J = VN.pibasis.basis1p.J
# ACE.evaluate(J, r)
# J


text:
abstract type AbstractSourceParameters end
abstract type AbstractSinkParameters end
abstract type AbstractStreamRequest end
abstract type AbstractStreamResponse end

function getdatatypes()
    return ["Boolean","Int32","Int64","Decimal","Double","DateTime","Guid","String"]
end

ext: jl

max_stars_repo_name: SnowdenNeuroverse/NeuroJulia

max_issues_repo_name: SnowdenNeuroverse/NeuroJulia

max_forks_repo_name: SnowdenNeuroverse/NeuroJulia


text:
__precompile__() #First Created on Fri Aug 14 20:12:01 2015

"""

# GeoEfficiency Package
Introduce a fast and flexible tool to calculate in batch or individually the `geometrical efficiency`
for a set of common radiation detectors shapes (cylindrical, Bore-hole, Well-type) as seen form
a source. The source can be a point, a disc or even a cylinder.

# Quick Usage
*  geoEff()	: Calculate the geometrical efficiency (GE) for one geometrical setup; returns a numerical value represent the calculated GE.

*  calc() 	: Calculate the geometrical efficiency (GE) for one geometrical setup and display full information on the console.

*  calcN()	: Calculate the geometrical efficiency (GE) for geometrical setup(s) recusevily and display full information on the console until the user quit.

*  batch()	: Calculate  in ``batch mode`` the geometrical efficiency using data in the **`$(join(split(dataDir,"/travis")))`** folder.
   For more information see `batch`, `batchInfo`.

# Documentation and Updates

\t Repository:    [`GitHub.com`](https://github.com/DrKrar/GeoEfficiency.jl/)
\t Documentation: https://GeoEfficiency.GitHub.io/dev/index.html
\t                https://juliahub.com/docs/GeoEfficiency/
\t PDF_Manual:    https://GeoEfficiency.GitHub.io/dev/GeoEfficiency.jl.pdf

To use Julia pakage manger to check for and obtaining the latest stable vesrion

```julia
julia> import Pkg

julia> Pkg.update("GeoEfficiency")
````

"""
module GeoEfficiency
const GeoEfficiency_Version = v"0.9.4-dev"

export

 # Config

 # Input_Console

 # Physics_Model
	Point,
	source,
	Detector,
	CylDetector,
	BoreDetector,
	WellDetector,

# Input_Batch
	getDetectors,
	setSrcToPoint,
	typeofSrc,

 # Calculations
	geoEff,

 # Output_Interface
 	max_batch,
 	calc,
	calcN,
	batch,
	batchInfo

include("Config.jl") # configurable parameters of the program.
include("Error.jl")	# define error system for the package.
include("Input_Console.jl")
include("Physics_Model.jl")
include("Input_Batch.jl")
include("Calculations.jl")
include("Output_Interface.jl")

#------------------------ splash ---------------------------

printstyled("""
\n
\t **************************************************
\t **            -=) GeoEfficiency (=-             **
\t **  Accurate Geometrical Efficiency Calculator  **
\t **************************************************

\t Author:     Mohamed E. Krar [ DrKrar@gmail.com ]
\t Profile:    https://www.researchgate.net/profile/Mohamed_Krar3

\t Welcome to GeoEfficiency $GeoEfficiency_Version

\t Start a calculation by using:
\t julia> calc()

\t Type ? just after 'julia>' to open help system:
\t julia>? GeoEfficiency
\n
""", color = :blue, bold = true)

end #module

max_stars_repo_name: DrKrar/GeoEfficiency

max_issues_repo_name: DrKrar/GeoEfficiency

max_forks_repo_name: DrKrar/GeoEfficiency


text:
using QuantLib

import Base.getindex, Base.endof

type TimeGrid
  times::Vector{Float64}
  dt::Vector{Float64}
  mandatoryTimes::Vector{Float64}
end

function TimeGrid(times::Vector{Float64}, steps::Int)
  sortedUniqueTimes = sort(unique(times))

  lastTime = sortedUniqueTimes[end]

  # TODO check if steps is 0
  dtMax = lastTime / steps
  periodBegin = 0.0
  times = zeros(1)

  for t in sortedUniqueTimes
    periodEnd = t
    if periodEnd != 0.0
      nSteps = Int(floor((periodEnd - periodBegin) / dtMax + 0.5))
      nSteps = nSteps != 0 ? nSteps : 1
      dt = (periodEnd - periodBegin) / nSteps

      tempTimes = zeros(nSteps)
      for n=1:nSteps
        tempTimes[n] = periodBegin + n * dt
      end
    end
    periodBegin = periodEnd
    times = vcat(times, tempTimes)
  end

  dt = diff(times)

  return TimeGrid(times, dt, sortedUniqueTimes)
end

function TimeGrid(endTime::Float64, steps::Int)
  endTime > 0.0 || error("negative times not allowed")
  dt = endTime / steps
  times = zeros(steps + 1)
  for i in eachindex(times)
    times[i] = dt * (i - 1)
  end

  mandatoryTimes = [endTime]

  dtVec = fill(dt, steps)

  return TimeGrid(times, dtVec, mandatoryTimes)
end

getindex(tg::TimeGrid, i::Int) = tg.times[i]
endof(tg::TimeGrid) = endof(tg.times)

is_empty(tg::TimeGrid) = length(tg.times) == 0

function closest_index(tg::TimeGrid, t::Float64)
  # stuff
  res = searchsortedfirst(tg.times, t)
  if res == 1
    return 1
  elseif res == length(tg.times) + 1
    return length(tg.times)
  else
    dt1 = tg.times[res] - t
    dt2 = t - tg.times[res - 1]
    if dt1 < dt2
      return res
    else
      return res - 1
    end
  end
end

closest_time(tg::TimeGrid, t::Float64) = tg.times[closest_index(tg, t)]

function return_index(tg::TimeGrid, t::Float64)
  i = closest_index(tg, t)
  if QuantLib.Math.close_enough(t, tg.times[i])
    return i
  else
    error("this time grid is wrong, $i $t $(tg.times[i]) $(tg.times[i-1]) $(tg.times[i+1])")
  end
end


text:
module MarketDatas
using DataFrames
using TimeSeries
using SQLite

export MarketData, storeMarketData, loadMarketData, loadHistoryData
export storeHistoryData
export DataSource, Logs, CREST, EveCentral, Unknown
@enum DataSource Logs CREST EveCentral Unknown
type MarketData
  itemID :: Int64
  source :: DataSource
  timeFetched :: DateTime
  sells :: AbstractDataFrame
  buys :: AbstractDataFrame
end

# MARKET DATA SAVING
const schema = ["""
CREATE TABLE metadata (
id INTEGER PRIMARY KEY ASC,
itemID INTEGER,
timeFetched TEXT
source INTEGER
);""",
"""CREATE TABLE sells (
id INTEGER PRIMARY KEY ASC,
mID INTEGER,
OrderID INTEGER,
Price INTEGER,
Location INTEGER,
Quantity INTEGER,
MinQuantity INTEGER,
Duration INTEGER,
Issued TEXT
);""",
"""CREATE TABLE buys (
id INTEGER PRIMARY KEY ASC,
mID INTEGER,
OrderID INTEGER,
Price INTEGER,
Location INTEGER,
Quantity INTEGER,
MinQuantity INTEGER,
Duration INTEGER,
Issued TEXT
);
""",
"""
CREATE TABLE history_data (
  id INTEGER PRIMARY KEY ASC,
  ItemID INTEGER,
  Date TEXT,
  LowPrice Real,
  HighPrice Real,
  AvgPrice Real,
  Volume INTEGER,
  OrderCount INTEGER,
  CONSTRAINT one_record UNIQUE (date, itemID)
);
"""]

# Code to save and load market data #
function initializeDB(filename :: ASCIIString)
  db = SQLite.DB(filename)
  for e in schema
    query(db, e)
  end
  db
end

function storeHistoryData(db, itemID :: Int, data :: TimeArray)
  q = "INSERT OR REPLACE INTO history_data " *
      "(ItemID, Date, LowPrice, HighPrice, AvgPrice, Volume, OrderCount) " *
      "VALUES (?, ?, ?, ?, ?, ?, ?)"
  s = SQLite.Stmt(db, q)
  for row = 1:length(data)
    SQLite.bind!(s, 1, itemID)
    SQLite.bind!(s, 2, string(data[row].timestamp[1]))
    SQLite.bind!(s, 3, data[row]["LowPrice"].values[1])
    SQLite.bind!(s, 4, data[row]["HighPrice"].values[1])
    SQLite.bind!(s, 5, data[row]["AvgPrice"].values[1])
    SQLite.bind!(s, 6, round(data[row]["Volume"].values[1]))
    SQLite.bind!(s, 7, round(data[row]["OrderCount"].values[1]))
    SQLite.execute!(s)
  end
end

function storeHistoryData(db :: SQLite.DB, data :: Dict{Int64, TimeArray})
  query(db, "BEGIN TRANSACTION")
  for (idx, i) in data
    storeHistoryData(db, idx, i)
  end
  query(db, "END TRANSACTION")
end

function storePriceData(db, table, mID, data :: DataFrame)
  q = "INSERT INTO $table "*
      "(mID, OrderID, Price, Location, Quantity, MinQuantity, Duration, Issued)"*
      "VALUES (?, ?, ?, ?, ?, ?, ?, ?)"
  s = SQLite.Stmt(db, q)
  for row = 1:length(data[1])
    SQLite.bind!(s, 1, mID)
    SQLite.bind!(s, 2, data[:OrderID][row])
    SQLite.bind!(s, 3, data[:Price][row])
    SQLite.bind!(s, 4, data[:Location][row])
    SQLite.bind!(s, 5, data[:Quantity][row])
    SQLite.bind!(s, 6, data[:MinQuantity][row])
    SQLite.bind!(s, 7, Int(data[:Duration][row]))
    SQLite.bind!(s, 8, string(data[:Issued][row]))
    SQLite.execute!(s)
  end
end

function storeMarketData(db :: SQLite.DB, data :: MarketData)
  query(db, "INSERT INTO metadata (itemID, timeFetched, source) VALUES (?, ?, ?)",
        [data.itemID, string(data.timeFetched), Int(data.source)])
  r = query(db, "SELECT last_insert_rowid()")
  mID :: Int = get(r.data[1][1])
  storePriceData(db, "sells", mID, data.sells)
  storePriceData(db,"buys", mID, data.buys)
end
function storeMarketData(db :: SQLite.DB, data)
  SQLite.query(db, "BEGIN TRANSACTION")
  for (idx,i :: MarketData) in data
    storeMarketData(db, i)
  end
  SQLite.query(db, "END TRANSACTION")
end
function storeMarketData(dbname :: ASCIIString, data)
  local db
  if isfile(dbname)
    db = SQLite.DB(dbname)
  else
    db = initializeDB(dbname)
  end
  storeMarketData(db, data)
end

function loadHistoryData(db, itemID :: Int)
  df = DataFrame(query(db, "SELECT * FROM history_data WHERE ItemID = ?", [itemID]))
  #convertMarketHistory(df)
end

function loadHistoryData(db, itemID :: Int, since :: DateTime)
  local timestamps :: Vector{DateTime}
  df = DataFrame(query(db, "SELECT * FROM history_data WHERE ItemID = ? AND Date > datetime(?)", [itemID, string(since)]))
  #convertMarketHistory(df)
end

function loadMarketData(db :: SQLite.DB, itemID :: Int, fresh :: DateTime)
  metadata = DataFrame(query(db, "SELECT * FROM metadata WHERE ItemID = ? AND timeFetched > datetime(?) ORDER BY timeFetched DESC", [itemID, string(since)]))
  (rows, cols) = size(metadata)
  if rows == 0
    return MarketData(0, DateTime(2003, 5, 5), Unknown,
                      DataFrame(OrderID = [], Price = [], Location = [], Quantity = [], MinQuantity = [], Duration = [], Issued = []),
                      DataFrame(OrderID = [], Price = [], Location = [], Quantity = [], MinQuantity = [], Duration = [], Issued = []))
  else
    df = Dates.DateFormat("yyyy-mm-ddTHH:MM:SS")
    sells = DataFrame(query(db, "SELECT * FROM sells WHERE mID = ?", [Int(metadata[:id][1])]))
    buys = DataFrame(query(db, "SELECT * FROM buys WHERE mID = ?", [Int(metadata[:id][1])]))
    return MarketData(Int(metadata[:itemID][1]),
                      DataSource(metadata[:source][1]),
                      DateTime(metadata[:timeFetched][1]).
                      sells,
                      buys)
  end
end
end


text:
using FactorGraph
using OMEinsum
using Zygote
using Test

function bpcheck(f, args...; η = 1e-5, verbose = false)
    g = gradient(f, args...)
    all(x->x===(nothing), g) && error()
    dy_ref = 0
    for x in g
        x === nothing && continue
        x isa Tuple && (dy_ref += η * mapreduce(y -> y === nothing ? 0 : sum(abs2,y), +, x))
        x isa AbstractArray && (dy_ref += η * sum(abs2,x))
    end
    dy = f(args...) - f([gi === nothing ? arg : arg .- η .* gi for (arg, gi) in zip(args,g)]...)

    verbose && @show dy
    verbose && @show dy_ref

    isapprox(dy, dy_ref, rtol=1e-2, atol=1e-8)
end

#Base.conj(arr::AbstractArray{<:Any, 0}) = conj!(copy(arr))

@testset "einsum bp" begin
    a = randn(ComplexF64, 3,3)
    f2(a) = ein"ij,ik,il->jkl"(a, conj(a), a) |> sum |> real

    @test gradient(f2, a)[1] !== nothing
    @test bpcheck(f2, a)
    @test bpcheck(a->treecontract(((1,2),3), ((1,2), (2,3), (3,1)), (a, a, a), ())[] |> real, a)
end

@testset "optcontract" begin
    a = randn(ComplexF64, 3,3)
    @test optcontract(((1,2), (2,3)), (a, a), (1,3)) ≈ a*a
    @test bpcheck(a->ein"ij,ji->"(a, a) |> sum |> real, a)
    @test bpcheck(a->optcontract(((1,2), (2,3)), (a, a), (1,3)) |> sum |> real, a)
    @test bpcheck(a->optcontract(((1,2), (2,3), (3,4)), (a,a,a), (4,1)) |> sum |> real, a)
end


text:
module NumericalDataManipulation

include("assertion/assertion.jl")
include("strings/strings.jl")
include("types/types.jl")
include("commonmath/commonmath.jl")
include("common/common.jl")
include("interpolation/interpolation.jl")
include("dataio/dataio.jl")
include("data/data.jl")

end # module


text:
using Documenter

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

deploydocs(
     repo = "github.com/PeaceFounder/SynchronicBallot.jl.git",
 )


text:
using Pkg
Pkg.activate(joinpath("test_env"))
Pkg.add("GitHub")
Pkg.instantiate()

using GitHub

ORG, REPO, PR = ENV["org"], ENV["repo"], ENV["pullrequest"]
TEST_RESULTS_FILE = "$(ORG)_$(REPO)_$(PR).txt"

# Need to add GITHUB_AUTH to your .bashrc
myauth = GitHub.authenticate(ENV["GITHUB_AUTH"])

function create_gist(authentication)

    file_content = read(TEST_RESULTS_FILE, String)

    file_dict = Dict(TEST_RESULTS_FILE => Dict("content" => file_content))
    gist = Dict{String,Any}(
        "description" => "Test results",
        "public" => true,
        "files" => file_dict,
    )

    posted_gist = GitHub.create_gist(params = gist, auth = authentication)

    return posted_gist
end

function post_gist_url_to_pr(comment::String; kwargs...)
    api = GitHub.DEFAULT_API
    repo = get_repo(api, ORG, REPO; kwargs...)
    pull_request = get_pull_request(api, ORG, repo, parse(Int, PR); kwargs...)
    GitHub.create_comment(api, repo, pull_request, comment; kwargs...)
end

function get_repo(api::GitHub.GitHubWebAPI, org::String, repo_name::String; kwargs...)
    my_params = Dict(:visibility => "all")
    # return GitHub.repo(api, repo; params = my_params, kwargs...)
    return Repo(
        GitHub.gh_get_json(
            api,
            "/repos/$(org)/$(repo_name)";
            params = my_params,
            kwargs...,
        ),
    )
end

function get_pull_request(
    api::GitHub.GitHubWebAPI,
    org::String,
    repo::Repo,
    pullrequest_id;
    kwargs...,
)
    my_params = Dict(:sort => "popularity", :direction => "desc")
    pull_request = PullRequest(
        GitHub.gh_get_json(
            api,
            "/repos/$(org)/$(repo.name)/pulls/$(pullrequest_id)";
            params = my_params,
            kwargs...,
        ),
    )

    return pull_request
end

function get_comment_from_test_results()
    open(TEST_RESULTS_FILE, "r") do file
        text_to_match = r"tests passed"
        for line in readlines(file)
            if occursin(text_to_match, line)
                return "$(strip(line)): "
            end
        end
        return "Tests failed: "
    end
end

comment = get_comment_from_test_results()
post_gist_url_to_pr("$comment $(create_gist(myauth).html_url)"; auth = myauth)

max_stars_repo_name: geoffroyleconte/QPModelCPLEX

max_issues_repo_name: geoffroyleconte/QPModelCPLEX

max_forks_repo_name: geoffroyleconte/QPModelCPLEX


text:
using Test, CarlemanLinearization
using DynamicPolynomials

@testset "Kronecker power (symbolic)" begin
    @polyvar x[1:2]
    y = kron_pow(x, 2)
    @test findall(y, x[1]*x[2]) == [2, 3]
    @test findall(y, x[1]^2) == [1]
    @test findall(y, x[2]^2) == [4]
    @test findall(y, x[2]^3) == Int[]
end


text:
"""

	F = G ∪ H

"""

struct DisjointUnion <: Constraint

	F::Variable
	G::Variable
	H::Variable

end

import Base.show
function Base.show(io::IO, unionConst::DisjointUnion)
	print(io, unionConst.F.name * " = " * unionConst.G.name * " ∪ " * unionConst.H.name)
end

function DisjointUnion(model::Model, F, G, H)
	constraint = DisjointUnion(F, G, H)
	push!(F.linkedConstraint, constraint)
	push!(G.linkedConstraint, constraint)
	push!(H.linkedConstraint, constraint)

	push!(model.constraints, constraint)

	constraint
end

function filtrage!(unionConst::DisjointUnion)

	stop =  !(isempty(intersect(unionConst.G.lowerBound, unionConst.H.lowerBound)))
	stop = stop || !(isempty(setdiff(union(unionConst.G.lowerBound, unionConst.H.lowerBound), unionConst.F.upperBound)))
	stop = stop || !(isempty(setdiff(unionConst.F.lowerBound, union(unionConst.G.upperBound, unionConst.H.upperBound))))

	changeVariable = Vector{Variable}(undef, 3)
	nbChange = 0

	changeF = Change(unionConst.F)
	changeG = Change(unionConst.G)
	changeH = Change(unionConst.H)

	if stop
		return [], (changeF, changeG, changeH), true
	else

		#Set
		if !unionConst.F.isFixed

			rem = setdiff(unionConst.F.upperBound, union(unionConst.G.upperBound, unionConst.H.upperBound))
			add = setdiff(union(unionConst.G.lowerBound, unionConst.H.lowerBound), unionConst.F.lowerBound)

			setdiff!(unionConst.F.upperBound, rem)
			union!(unionConst.F.lowerBound, add)

			cardInf = max(unionConst.F.cardinalInf, length(unionConst.F.lowerBound), unionConst.G.cardinalInf, unionConst.H.cardinalInf)
			cardSup = min(unionConst.F.cardinalSup, length(unionConst.F.upperBound), unionConst.G.cardinalSup + unionConst.H.cardinalSup)

			changeF.cardAdded = cardInf - unionConst.F.cardinalInf
			changeF.cardRemoved = cardSup - unionConst.F.cardinalSup

			unionConst.F.cardinalInf += changeF.cardAdded
			unionConst.F.cardinalSup += changeF.cardRemoved

			if !(isempty(rem) && isempty(add))
				nbChange += 1
				changeVariable[nbChange] = unionConst.F

				if isFixed(unionConst.F)
					fix!(unionConst.F, changeF)
				end
			end

			for elt in add
				push!(changeF.added, elt)
			end
			for elt in rem
				push!(changeF.removed, elt)
			end
		end

		if !unionConst.G.isFixed

			rem = intersect(unionConst.G.upperBound, unionConst.F.lowerBound, unionConst.H.lowerBound)
			add = setdiff(setdiff(unionConst.F.lowerBound, unionConst.H.upperBound), unionConst.G.lowerBound)

			setdiff!(unionConst.G.upperBound, rem)
			union!(unionConst.G.lowerBound, add)

			cardInf = max(unionConst.G.cardinalInf, length(unionConst.G.lowerBound))
			cardSup = min(unionConst.G.cardinalSup, length(unionConst.G.upperBound), unionConst.F.cardinalSup - unionConst.H.cardinalInf)

			changeG.cardAdded = cardInf - unionConst.G.cardinalInf
			changeG.cardRemoved = cardSup - unionConst.G.cardinalSup

			unionConst.G.cardinalInf += changeG.cardAdded
			unionConst.G.cardinalSup += changeG.cardRemoved

			if !(isempty(rem) && isempty(add))
				nbChange += 1
				changeVariable[nbChange] = unionConst.G

				if isFixed(unionConst.G)
					fix!(unionConst.G, changeG)
				end
			end

			for elt in add
				push!(changeG.added, elt)
			end
			for elt in rem
				push!(changeG.removed, elt)
			end
		end

		if !unionConst.H.isFixed
			rem = intersect(unionConst.H.upperBound, unionConst.F.lowerBound, unionConst.G.lowerBound)
			add = setdiff(setdiff(unionConst.F.lowerBound, unionConst.G.upperBound), unionConst.H.lowerBound)

			setdiff!(unionConst.H.upperBound, rem)
			union!(unionConst.H.lowerBound, add)

			cardInf = max(unionConst.H.cardinalInf, length(unionConst.H.lowerBound))
			cardSup = min(unionConst.H.cardinalSup, length(unionConst.H.upperBound), unionConst.F.cardinalSup - unionConst.G.cardinalInf)

			changeH.cardAdded = cardInf - unionConst.H.cardinalInf
			changeH.cardRemoved = cardSup - unionConst.H.cardinalSup

			unionConst.H.cardinalInf += changeH.cardAdded
			unionConst.H.cardinalSup += changeH.cardRemoved

			if !(isempty(rem) && isempty(add))
				nbChange += 1
				changeVariable[nbChange] = unionConst.H

				if isFixed(unionConst.H)
					fix!(unionConst.H, changeH)
				end
			end

			for elt in add
				push!(changeH.added, elt)
			end
			for elt in rem
				push!(changeH.removed, elt)
			end
		end

		#Error
		stop = stop || !(unionConst.F.cardinalInf <= unionConst.F.cardinalSup)
		stop = stop || !(unionConst.G.cardinalInf <= unionConst.G.cardinalSup)
		stop = stop || !(unionConst.H.cardinalInf <= unionConst.H.cardinalSup)

		changes = (changeF, changeG, changeH)
		for c in changes
			if !(isempty(c.added))
				#println("$unionConst ajoute $(c.added) à $(c.var)")
			end
			if !(isempty(c.removed))
				#println("$unionConst enlève $(c.removed) à $(c.var)")
			end
		end

		return changeVariable[1:nbChange], (changeF, changeG, changeH), stop
	end

end


text:
"""
    Cylinder Geoemtry (with homogeneous distribution)

The Cylinder with heigh `h`, that goes from [-h/2, h/2], and Radius `R`
"""
function Atom(
    geometry::Cylinder,
    N::Int64,
    R::Union{Real,Integer},
    h::Union{Real,Integer};
    createFunction = ftn_AtomsOnCylinder::Function,
)
    @debug "start: Shape - Cylinder"

    dimensions = 3
    ρ = N/( h*π*R^2 )
    rₘᵢₙ = get_rₘᵢₙ(ρ)
    if rₘᵢₙ ≥ R / 10
        rₘᵢₙ = R / 100
    end

    r = get_atoms(dimensions, N, rₘᵢₙ; createFunction, R, h)
    sizes = (R=Float64(R), h=Float64(h))
    @debug "end  : Shape - Cylinder"
    return Atom(Cylinder(), r, N, sizes)
end

function Atom(geometry::Cylinder, r::Matrix, R::Union{Real,Integer}, h::Union{Real,Integer})
    N = size(r, 2) # remember to use each collum as a atom position
    return Atom(Cylinder(), r, N, (R=Float64(R), h=Float64(h)) )
end

function ftn_AtomsOnCylinder(; kwargs...)
    R, h = kwargs[:R], kwargs[:h]

    r,θ = R*rand(), 2π*rand()

    x = sqrt(r)*cos(θ)
    y = sqrt(r)*sin(θ)
    z = -h*rand() + h/2

    return [x, y, z]
end


text:
module RobustDMD

using LinearAlgebra
using Random
using Printf
using Optim
using Munkres
using StatsBase

include("DMD_BFGS.jl")
include("DMD_Func.jl")
include("DMD_Loss.jl")
include("DMD_SVRG.jl")
include("DMD_Type.jl")
include("DMD_Util.jl")

end


text:
# Steam Buffer
# ============
#
# Stream buffer type.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"""
A buffer in which pushing to the back, popping from the front (shift),
and indexing is efficient.

It's constructed as a circular buffer that resizes when it runs out of space.
"""
type StreamBuffer{T}
    data::Vector{T}
    first_idx::Int
    last_idx::Int
    len::Int

    function StreamBuffer()
        return new(Array(T, 16), 1, 0, 0)
    end
end

function Base.length{T}(buf::StreamBuffer{T})
    return buf.len
end

function Base.isempty{T}(buf::StreamBuffer{T})
    return buf.last_idx == 0
end

function Base.endof{T}(buf::StreamBuffer{T})
    return length(buf)
end

function Base.push!{T}(buf::StreamBuffer{T}, x::T)
    if length(buf) == length(buf.data)
        newdata = Array(T, 2 * length(buf.data))
        len = length(buf)
        for i in 1:len
            newdata[i] = buf[i]
        end
        buf.data = newdata
        buf.first_idx = 1
        buf.last_idx = len
    end

    buf.last_idx += 1
    if buf.last_idx > length(buf.data)
        buf.last_idx -= length(buf.data)
    end
    @inbounds buf.data[buf.last_idx] = x
    buf.len += 1
    return x
end

function Base.shift!{T}(buf::StreamBuffer{T})
    if buf.last_idx == 0
        throw(BoundsError())
    end

    @inbounds x = buf.data[buf.first_idx]

    if buf.first_idx == buf.last_idx
        buf.last_idx = 0
        buf.first_idx = 1
    elseif buf.first_idx == length(buf.data)
        buf.first_idx = 1
    else
        buf.first_idx += 1
    end

    buf.len -= 1
    return x
end

function Base.getindex{T}(buf::StreamBuffer{T}, idx::Int)
    if idx > length(buf) || idx < 1
        throw(BoundsError())
    end

    idx = buf.first_idx + idx - 1
    if idx > length(buf.data)
        idx -= length(buf.data)
    end
    @inbounds x = buf.data[idx]
    return x
end


text:
## ExoplanetsSysSim/src/window_function.jl
## (c) 2018 Darin Ragozzine

# Gather and prepare the window function data
module WindowFunction

export setup_window_function, get_window_function_data, get_window_function_id, eval_window_function

#using DataArrays
using DataFrames
#using CSV
using FileIO
using ExoplanetsSysSim
import ExoplanetsSysSim.SimulationParameters

# Object to hold window function data
struct window_function_data
  window_func_array::Array{Float64,3}      # Value of window function (window_function_id, duration_id, period_id).  Maybe rename to wf_value or data?
  wf_durations_in_hrs::Array{Float64,1}    # Boundaries for duration bins in window_func_array.  Maybe rename to durations?
  wf_periods_in_days::Array{Float64,1}     # Boundaries for periods bins in window_func_array.  Maybe rename to periods?
  sorted_quarter_strings::Array{Int64,1}   # TODO OPT: Is there a reason to keep this?  Maybe rename to quarter_strings?
  allsortedkepids::Array{Int64,1}          # value is Kepler ID.  Index is same as index to window_function_id_arr
  window_function_id_arr::Array{Int64,1}   # value is index for window_func_array.  Index is same as index to allsortedkepids
  default_wf_id::Int64                     # Index corresponding to the default window function
  already_warned::Array{Bool,1}            # Wether we've already thrown a warning about this kepid
end

function window_function_data()
  window_function_data( Array{Float64,3}(undef,0,0,0), Array{Float64,1}(undef,0),Array{Float64,1}(undef,0), Array{Int64,1}(undef,0),Array{Int64,1}(undef,0),Array{Int64,1}(undef,0), 0, falses(0) )
end

win_func_data = window_function_data()

function setup(sim_param::SimParam; force_reread::Bool = false)
  global win_func_data
  if haskey(sim_param,"read_window_function") && !force_reread
     return win_func_data
  end
  window_function_filename = convert(String,joinpath(dirname(pathof(ExoplanetsSysSim)), "..","data", convert(String,get(sim_param,"window_function","DR25topwinfuncs.jld2")) ) )
  setup(window_function_filename)
  add_param_fixed(sim_param,"read_window_function",true)
  @assert( size(win_func_data.window_func_array,2) == length(win_func_data.wf_durations_in_hrs) )
  @assert( size(win_func_data.window_func_array,3) == length(win_func_data.wf_periods_in_days) )
  @assert( size(win_func_data.window_func_array,1) >= maximum(win_func_data.window_function_id_arr) )
  @assert( size(win_func_data.window_func_array,1) >= win_func_data.default_wf_id )
  return win_func_data
end

function setup(filename::String)
# Reads in the window function data collected from the Kepler Completeness Products
# see Darin Ragozzine's get/cleanDR25winfuncs.jl

  if occursin(r".jld2$",filename)
    try
      wfdata = load(filename)
      window_func_array = wfdata["window_func_array"]
      wf_durations_in_hrs = wfdata["wf_durations_in_hrs"]  # TODO OPT DETAIL: Should we convert units to days here?
      wf_periods_in_days = wfdata["wf_periods_in_days"]
      sorted_quarter_strings = wfdata["sorted_quarter_strings"]
      allsortedkepids = wfdata["allsortedkepids"]
      window_function_id_arr = wfdata["window_function_id_arr"]
      already_warned = falses(length(allsortedkepids))
      global win_func_data = window_function_data(window_func_array, wf_durations_in_hrs, wf_periods_in_days, sorted_quarter_strings,
                                                allsortedkepids, window_function_id_arr, maximum(window_function_id_arr), already_warned )

    catch
      error(string("# Failed to read window function data > ", filename," < in jld2 format."))
    end
  end

  return win_func_data
end

setup_window_function(sim_param::SimParam; force_reread::Bool = false) = setup(sim_param, force_reread=force_reread)
setup_window_function(filename::String; force_reread::Bool = false) = setup(filename, force_reread=force_reread)

function get_window_function_data()::window_function_data
   #global win_func_data
   return win_func_data
end

function get_window_function_id(kepid::Int64; use_default_for_unknown::Bool = true)::Int64
  # takes the quarter string from the stellar catalog and determines the window function id
  # from DR25topwinfuncs.jld2 made by Darin Ragozzine's cleanDR25winfuncs.jl script.
  no_win_func_available::Int64 = -1        # hardcoding this in, should match convention in window function input file

  idx = searchsortedfirst(win_func_data.allsortedkepids,kepid) # all Kepler kepids are in allsortedkepids
  wf_id = win_func_data.window_function_id_arr[idx]

  if wf_id == no_win_func_available && use_default_for_unknown
    # if a target is observed for less than 4 quarters, then it won't have a corresponding
    # window function in this list, so throw a warning and use the last window_function_id
    # which corresponds to an "averaged" window function
    if !win_func_data.already_warned[idx]
       win_func_data.already_warned[idx] = true
       if sum(win_func_data.already_warned) < 20
          @warn "Window function data is not avaialble for kepid $kepid, using default."
       end
    end
    wf_id = win_func_data.default_wf_id
  end
  # TODO SCI DETAIL IMPORTANT? This does not include TPS timeouts or MESthresholds (see DR25 Completeness Products)

  return wf_id
end

function calc_period_idx(P::Float64)::Int64
  @assert(P>zero(P))
  idx = searchsortedlast(win_func_data.wf_periods_in_days,P)
  if idx == 0
     return 1
  elseif idx<length(win_func_data.wf_periods_in_days)
     if P-win_func_data.wf_periods_in_days[idx]>win_func_data.wf_periods_in_days[idx+1]-P
        idx += 1
     end
  end
  return idx   # TODO IMPORTANT: IMPLEMENT / TEST
end

function calc_duration_idx(D::Float64)::Int64
  # NOTE IMPORTANT: Currently assumes we left wf data in hours, so deal with that conversion here
  @assert(D>=zero(D)) ##### Make sure this function is still doing the right thing if D = 0!
  hours_in_day = 24
  idx = searchsortedlast(win_func_data.wf_durations_in_hrs,D*hours_in_day)
  if idx == 0
     return 1
  elseif idx<length(win_func_data.wf_durations_in_hrs)
     if D*hours_in_day-win_func_data.wf_durations_in_hrs[idx]>win_func_data.wf_durations_in_hrs[idx+1]-D*hours_in_day
        idx += 1
     end
  end
  return idx   # TODO IMPORTANT: IMPLEMENT / TEST
end

function eval_window_function(wf_idx::Int64=-1; Duration::Float64=0., Period::Float64=0.)::Float64
  D_idx = calc_duration_idx(Duration)
  P_idx = calc_period_idx(Period)
  wf = eval_window_function(wf_idx,D_idx,P_idx)
  # TODO IMPORTANT: Improve way deal with missing wf values for some durations. Interpolate?
  while wf<=zero(wf) && D_idx<length(win_func_data.wf_durations_in_hrs)
     D_idx += 1
     wf = eval_window_function(wf_idx,D_idx,P_idx)
  end
  return wf
end

function eval_window_function(wf_idx::Int64, D_idx::Int64, P_idx::Int64)::Float64
   global win_func_data
   #@assert(1<=wf_idx<maximum(win_func_data.window_function_id_arr))
   #@assert(1<=P_idx<=length(win_func_data.wf_periods_in_days))
   #@assert(1<=D_idx<=length(win_func_data.wf_durations_in_hrs))
   return win_func_data.window_func_array[wf_idx,D_idx,P_idx]
end

end  # module WindowFunction


text:
using JuMP

m = Model()

type Recipe
	name::AbstractString
	time::Float64
	outs::Dict{AbstractString, Float64}
	ins::Dict{AbstractString, Float64}
	Recipe(name, time, speed) = new(name, time/speed, Dict{AbstractString, Float64}(), Dict{AbstractString, Float64}())
end

function production(r::Recipe, nfacs)
	@printf "%s: %0.1f\n" r.name nfacs
	@printf "\tProduces\n"
	for (k,v) in r.outs
		@printf "\t\t%0.2f %s per min / %0.2f per sec\n" 60v * nfacs k v*nfacs
	end
	@printf "\tConsumes\n"
	for (k,v) in r.ins
		@printf "\t\t%0.2f %s per min / %0.2f per sec\n" 60v * nfacs k v*nfacs
	end
end

@enum Rnum RED GREEN WIRE

recips = Dict{Rnum, Recipe}()

macro IN(rnum, iname, iqty)
	:(recips[$rnum].ins[$iname] = $iqty / recips[$rnum].time)
end

macro OUT(rnum, iname, iqty)
	:(recips[$rnum].outs[$iname] = $iqty / recips[$rnum].time)
end

recips[GREEN] = Recipe("Green", 0.5, 1.25)
@IN GREEN "Iron" 1
@IN GREEN "Wire" 3
@OUT GREEN "Green" 1

recips[RED] = Recipe("Red", 6.0, 1.25)
@IN RED "Green" 2
@IN RED "Plastic" 2
@IN RED "Wire" 4
@OUT RED "Red" 1

recips[WIRE] = Recipe("Wire", 0.5, 1.25)
@IN WIRE "Copper" 1
@OUT WIRE "Wire" 2

@variable(m, REDs >= 0) #, Int)
@variable(m, GREENs >= 0) #, Int)
@variable(m, WIREs >= 0) #, Int)

macro FIN(Fn, F, R) :($Fn * recips[$F].ins[$R]) end
macro FOUT(Fn, F, R) :($Fn * recips[$F].outs[$R] )end

@constraint(m, REDs == 6)

@constraint(m, @FIN(REDs, RED, "Green") == @FOUT(GREENs, GREEN, "Green"))
@constraint(m, @FIN(REDs, RED, "Wire") + @FIN(GREENs, GREEN, "Wire") == @FOUT(WIREs, WIRE, "Wire"))

@objective(m, Max, sum([REDs, GREENs, WIREs]))
solve(m)

production(recips[RED], getvalue(REDs))
production(recips[GREEN], getvalue(GREENs))
production(recips[WIRE], getvalue(WIREs))

max_stars_repo_name: lawless-m/factorio-utils

max_issues_repo_name: lawless-m/factorio-utils

max_forks_repo_name: lawless-m/factorio-utils


text:
import ClimaCore:
    Fields,
    Domains,
    Topologies,
    Meshes,
    DataLayouts,
    Operators,
    Geometry,
    Spaces

using OrdinaryDiffEq: ODEProblem, solve, SSPRK33

import Logging
import TerminalLoggers
Logging.global_logger(TerminalLoggers.TerminalLogger())
const CI = !isnothing(get(ENV, "CI", nothing))

const FT = Float64

a = FT(0.0)
b = FT(1.0)
n = 10
α = FT(0.1)

domain = Domains.IntervalDomain(
    Geometry.ZPoint{FT}(a),
    Geometry.ZPoint{FT}(b),
    boundary_names = (:bottom, :top),
)
mesh = Meshes.IntervalMesh(domain, nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(mesh)
T = Fields.zeros(FT, cs)

# Solve Heat Equation: ∂_t T = α ∇²T
function ∑tendencies!(dT, T, _, t)

    bcs_bottom = Operators.SetValue(FT(0.0))
    bcs_top = Operators.SetGradient(Geometry.WVector(FT(1.0)))

    gradc2f = Operators.GradientC2F(bottom = bcs_bottom, top = bcs_top)
    divf2c = Operators.DivergenceF2C()

    return @. dT = α * divf2c(gradc2f(T))
end

@show ∑tendencies!(similar(T), T, nothing, 0.0);

# Solve the ODE operator
Δt = 0.02

prob = ODEProblem(∑tendencies!, T, (0.0, 10.0))
sol = solve(
    prob,
    SSPRK33(),
    dt = Δt,
    saveat = 10 * Δt,
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);

ENV["GKSwstype"] = "nul"
using ClimaCorePlots, Plots
Plots.GRBackend()

dirname = "heat"
path = joinpath(@__DIR__, "output", dirname)
mkpath(path)

anim = Plots.@animate for u in sol.u
    Plots.plot(u, xlim = (0, 1))
end
Plots.mp4(anim, joinpath(path, "heat.mp4"), fps = 10)
Plots.png(Plots.plot(sol.u[end], xlim = (0, 1)), joinpath(path, "heat_end.png"))

function linkfig(figpath, alt = "")
    # buildkite-agent upload figpath
    # link figure in logs if we are running on CI
    if get(ENV, "BUILDKITE", "") == "true"
        artifact_url = "artifact://$figpath"
        print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
    end
end

linkfig(
    relpath(joinpath(path, "heat_end.png"), joinpath(@__DIR__, "../..")),
    "Heat End Simulation",
)


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you 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.

"""
    Arrow.FixedSizeList

An `ArrowVector` where each element is a "fixed size" list of some kind, like a `NTuple{N, T}`.
"""
struct FixedSizeList{T, A <: AbstractVector} <: ArrowVector{T}
    arrow::Vector{UInt8} # need to hold a reference to arrow memory blob
    validity::ValidityBitmap
    data::A
    ℓ::Int
    metadata::Union{Nothing, Base.ImmutableDict{String,String}}
end

Base.size(l::FixedSizeList) = (l.ℓ,)

@propagate_inbounds function Base.getindex(l::FixedSizeList{T}, i::Integer) where {T}
    @boundscheck checkbounds(l, i)
    S = Base.nonmissingtype(T)
    X = ArrowTypes.ArrowKind(ArrowTypes.ArrowType(S))
    N = ArrowTypes.getsize(X)
    Y = ArrowTypes.gettype(X)
    if X !== T && !(l.validity[i])
        return missing
    else
        off = (i - 1) * N
        if X === T && isbitstype(Y)
            tup = _unsafe_load_tuple(NTuple{N, Y}, l.data, off + 1)
        else
            tup = ntuple(j->l.data[off + j], N)
        end
        return ArrowTypes.fromarrow(T, tup)
    end
end

function _unsafe_load_tuple(::Type{NTuple{N,T}}, bytes::Vector{UInt8}, i::Integer) where {N,T}
    x = Ref(bytes, i)
    y = Ref{NTuple{N,T}}()
    ArrowTypes._unsafe_cast!(y, x, N)
    return y[]
end

@propagate_inbounds function Base.setindex!(l::FixedSizeList{T}, v::T, i::Integer) where {T}
    @boundscheck checkbounds(l, i)
    if v === missing
        @inbounds l.validity[i] = false
    else
        N = ArrowTypes.getsize(ArrowTypes.ArrowKind(ArrowTypes.ArrowType(Base.nonmissingtype(T))))
        off = (i - 1) * N
        foreach(1:N) do j
            @inbounds l.data[off + j] = v[j]
        end
    end
    return v
end

# lazy equal-spaced flattener
struct ToFixedSizeList{T, N, A} <: AbstractVector{T}
    data::A # A is AbstractVector of (AbstractVector or AbstractString)
end

origtype(::ToFixedSizeList{T, N, A}) where {T, N, A} = eltype(A)

function ToFixedSizeList(input)
    NT = ArrowTypes.ArrowKind(Base.nonmissingtype(eltype(input))) # typically NTuple{N, T}
    return ToFixedSizeList{ArrowTypes.gettype(NT), ArrowTypes.getsize(NT), typeof(input)}(input)
end

Base.IndexStyle(::Type{<:ToFixedSizeList}) = Base.IndexLinear()
Base.size(x::ToFixedSizeList{T, N}) where {T, N} = (N * length(x.data),)

Base.@propagate_inbounds function Base.getindex(A::ToFixedSizeList{T, N}, i::Integer) where {T, N}
    @boundscheck checkbounds(A, i)
    a, b = fldmod1(i, N)
    @inbounds x = A.data[a]
    return @inbounds x === missing ? ArrowTypes.default(T) : x[b]
end

# efficient iteration
@inline function Base.iterate(A::ToFixedSizeList{T, N}, (i, chunk, chunk_i, len)=(1, 1, 1, length(A))) where {T, N}
    i > len && return nothing
    @inbounds y = A.data[chunk]
    @inbounds x = y === missing ? ArrowTypes.default(T) : y[chunk_i]
    if chunk_i == N
        chunk += 1
        chunk_i = 1
    else
        chunk_i += 1
    end
    return x, (i + 1, chunk, chunk_i, len)
end

arrowvector(::FixedSizeListKind, x::FixedSizeList, i, nl, fi, de, ded, meta; kw...) = x

function arrowvector(::FixedSizeListKind{N, T}, x, i, nl, fi, de, ded, meta; kw...) where {N, T}
    len = length(x)
    validity = ValidityBitmap(x)
    flat = ToFixedSizeList(x)
    if eltype(flat) == UInt8
        data = flat
        S = origtype(flat)
    else
        data = arrowvector(flat, i, nl + 1, fi, de, ded, nothing; kw...)
        S = withmissing(eltype(x), NTuple{N, eltype(data)})
    end
    return FixedSizeList{S, typeof(data)}(UInt8[], validity, data, len, meta)
end

function compress(Z::Meta.CompressionType, comp, x::FixedSizeList{T, A}) where {T, A}
    len = length(x)
    nc = nullcount(x)
    validity = compress(Z, comp, x.validity)
    buffers = [validity]
    children = Compressed[]
    if eltype(A) == UInt8
        push!(buffers, compress(Z, comp, x.data))
    else
        push!(children, compress(Z, comp, x.data))
    end
    return Compressed{Z, typeof(x)}(x, buffers, len, nc, children)
end

function makenodesbuffers!(col::FixedSizeList{T, A}, fieldnodes, fieldbuffers, bufferoffset, alignment) where {T, A}
    len = length(col)
    nc = nullcount(col)
    push!(fieldnodes, FieldNode(len, nc))
    @debug 1 "made field node: nodeidx = $(length(fieldnodes)), col = $(typeof(col)), len = $(fieldnodes[end].length), nc = $(fieldnodes[end].null_count)"
    # validity bitmap
    blen = nc == 0 ? 0 : bitpackedbytes(len, alignment)
    push!(fieldbuffers, Buffer(bufferoffset, blen))
    @debug 1 "made field buffer: bufferidx = $(length(fieldbuffers)), offset = $(fieldbuffers[end].offset), len = $(fieldbuffers[end].length), padded = $(padding(fieldbuffers[end].length, alignment))"
    bufferoffset += blen
    if eltype(A) === UInt8
        blen = ArrowTypes.getsize(ArrowTypes.ArrowKind(Base.nonmissingtype(T))) * len
        push!(fieldbuffers, Buffer(bufferoffset, blen))
        @debug 1 "made field buffer: bufferidx = $(length(fieldbuffers)), offset = $(fieldbuffers[end].offset), len = $(fieldbuffers[end].length), padded = $(padding(fieldbuffers[end].length, alignment))"
        bufferoffset += padding(blen, alignment)
    else
        bufferoffset = makenodesbuffers!(col.data, fieldnodes, fieldbuffers, bufferoffset, alignment)
    end
    return bufferoffset
end

function writebuffer(io, col::FixedSizeList{T, A}, alignment) where {T, A}
    @debug 1 "writebuffer: col = $(typeof(col))"
    @debug 2 col
    writebitmap(io, col, alignment)
    # write values array
    if eltype(A) === UInt8
        n = writearray(io, UInt8, col.data)
        @debug 1 "writing array: col = $(typeof(col.data)), n = $n, padded = $(padding(n, alignment))"
        writezeros(io, paddinglength(n, alignment))
    else
        writebuffer(io, col.data, alignment)
    end
    return
end

max_stars_repo_name: jonalm/arrow-julia

max_issues_repo_name: jonalm/arrow-julia

max_forks_repo_name: jonalm/arrow-julia


text:
"""
$(TYPEDEF)

A system of difference equations.

# Fields
$(FIELDS)

# Example

```
using ModelingToolkit

@parameters σ ρ β
@variables t x(t) y(t) z(t) next_x(t) next_y(t) next_z(t)

eqs = [next_x ~ σ*(y-x),
       next_y ~ x*(ρ-z)-y,
       next_z ~ x*y - β*z]

de = DiscreteSystem(eqs,t,[x,y,z],[σ,ρ,β])
```
"""
struct DiscreteSystem <: AbstractTimeDependentSystem
    """The differential equations defining the discrete system."""
    eqs::Vector{Equation}
    """Independent variable."""
    iv::Sym
    """Dependent (state) variables. Must not contain the independent variable."""
    states::Vector
    """Parameter variables. Must not contain the independent variable."""
    ps::Vector
    """Array variables."""
    var_to_name
    """Control parameters (some subset of `ps`)."""
    ctrls::Vector
    """Observed states."""
    observed::Vector{Equation}
    """
    Name: the name of the system
    """
    name::Symbol
    """
    systems: The internal systems. These are required to have unique names.
    """
    systems::Vector{DiscreteSystem}
    """
    default_u0: The default initial conditions to use when initial conditions
    are not supplied in `DiscreteSystem`.
    """
    default_u0::Dict
    """
    default_p: The default parameters to use when parameters are not supplied
    in `DiscreteSystem`.
    """
    default_p::Dict
    function DiscreteSystem(discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name, systems, default_u0, default_p)
        check_variables(dvs,iv)
        check_parameters(ps,iv)
        new(discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name, systems, default_u0, default_p)
    end
end

"""
    $(TYPEDSIGNATURES)

Constructs a DiscreteSystem.
"""
function DiscreteSystem(
                   eqs::AbstractVector{<:Equation}, iv, dvs, ps;
                   controls = Num[],
                   observed = Num[],
                   systems = DiscreteSystem[],
                   name=nothing,
                   default_u0=Dict(),
                   default_p=Dict(),
                   defaults=_merge(Dict(default_u0), Dict(default_p)),
                  )
    name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
    eqs = collect(eqs)
    iv′ = value(iv)
    dvs′ = value.(dvs)
    ps′ = value.(ps)
    ctrl′ = value.(controls)

    if !(isempty(default_u0) && isempty(default_p))
        Base.depwarn("`default_u0` and `default_p` are deprecated. Use `defaults` instead.", :ODESystem, force=true)
    end
    defaults = todict(defaults)
    defaults = Dict(value(k) => value(v) for (k, v) in pairs(defaults))

    var_to_name = Dict()
    process_variables!(var_to_name, defaults, dvs′)
    process_variables!(var_to_name, defaults, ps′)

    sysnames = nameof.(systems)
    if length(unique(sysnames)) != length(sysnames)
        throw(ArgumentError("System names must be unique."))
    end
    DiscreteSystem(eqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, name, systems, default_u0, default_p)
end

"""
    $(TYPEDSIGNATURES)

Generates an DiscreteProblem from an DiscreteSystem.
"""
function DiffEqBase.DiscreteProblem(sys::DiscreteSystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    eval_module = @__MODULE__,
                                    eval_expression = true,
                                    kwargs...)
    dvs = states(sys)
    ps = parameters(sys)
    eqs = equations(sys)
    eqs = linearize_eqs(sys, eqs)
    # defs = defaults(sys)
    t = get_iv(sys)
    u0 = varmap_to_vars(u0map,dvs)
    rhss = [eq.rhs for eq in eqs]
    u = dvs
    p = varmap_to_vars(parammap,ps)

    f_gen = generate_function(sys; expression=Val{eval_expression}, expression_module=eval_module)
    f_oop, _ = (@RuntimeGeneratedFunction(eval_module, ex) for ex in f_gen)
    f(u,p,t) = f_oop(u,p,t)
    DiscreteProblem(f,u0,tspan,p;kwargs...)
end

function linearize_eqs(sys, eqs=get_eqs(sys); return_max_delay=false)
    unique_states = unique(operation.(states(sys)))
    max_delay = Dict(v=>0.0 for v in unique_states)

    r = @rule ~t::(t -> istree(t) && any(isequal(operation(t)), operation.(states(sys))) && is_delay_var(get_iv(sys), t)) => begin
        delay = get_delay_val(get_iv(sys), first(arguments(~t)))
        if delay > max_delay[operation(~t)]
            max_delay[operation(~t)] = delay
        end
        nothing
    end
    SymbolicUtils.Postwalk(r).(rhss(eqs))

    if any(values(max_delay) .> 0)

        dts = Dict(v=>Any[] for v in unique_states)
        state_ops = Dict(v=>Any[] for v in unique_states)
        for v in unique_states
            for eq in eqs
                if isdifferenceeq(eq) && istree(arguments(eq.lhs)[1]) && isequal(v, operation(arguments(eq.lhs)[1]))
                    append!(dts[v], [operation(eq.lhs).dt])
                    append!(state_ops[v], [operation(eq.lhs)])
                end
            end
        end

        all(length.(unique.(values(state_ops))) .<= 1) || error("Each state should be used with single difference operator.")

        dts_gcd = Dict()
        for v in keys(dts)
            dts_gcd[v] = (length(dts[v]) > 0) ? first(dts[v]) : nothing
        end

        lin_eqs = [
            v(get_iv(sys) - (t)) ~ v(get_iv(sys) - (t-dts_gcd[v]))
            for v in unique_states if max_delay[v] > 0 && dts_gcd[v]!==nothing for t in collect(max_delay[v]:(-dts_gcd[v]):0)[1:end-1]
        ]
        eqs = vcat(eqs, lin_eqs)
    end
    if return_max_delay return eqs, max_delay end
    eqs
end

function get_delay_val(iv, x)
    delay = x - iv
    isequal(delay > 0, true) && error("Forward delay not permitted")
    return -delay
end

check_difference_variables(eq) = check_operator_variables(eq, Difference)

function generate_function(
        sys::DiscreteSystem, dvs = states(sys), ps = parameters(sys);
        kwargs...
    )
    eqs = equations(sys)
    foreach(check_difference_variables, eqs)
    # substitute x(t) by just x
    rhss = [eq.rhs for eq in eqs]

    u = map(x->time_varying_as_func(value(x), sys), dvs)
    p = map(x->time_varying_as_func(value(x), sys), ps)
    t = get_iv(sys)

    build_function(rhss, u, p, t; kwargs...)
end


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/iup/include/cdiup.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@static Sys.iswindows() ? (const libiupcd_ = "iupcd") : (const libiupcd_ = "libiupcd")  # Name of IUP shared lib.

function cdContextIup()
    ccall((:cdContextIup, libiupcd_), Ptr{cdContext}, (), )
end

# Was (alone) in cdiup_h.jl
const CD_IUP = cdContextIup


text:
include("kepler_drift_solver.jl")

# Takes a single kepler step, with reverse drift, calling Wisdom & Hernandez solver
#
function kepler_drift_step!(gm::T,h::T,state0::Array{T,1},state::Array{T,1},drift_first::Bool) where {T <: Real}
# compute beta, r0,  get x/v from state vector & call correct subroutine
x0 = zeros(eltype(state0),3)
v0 = zeros(eltype(state0),3)
for k=1:3
  x0[k]=state0[k+1]
  v0[k]=state0[k+4]
end
s0=state0[11]
iter = kep_drift_ell_hyp!(x0,v0,gm,h,s0,state,drift_first)
return
end

function kepler_drift_step!(gm::T,h::T,state0::Array{T,1},state::Array{T,1},jacobian::Array{T,2},drift_first::Bool,d::Derivatives{T}) where {T <: Real}
# compute beta, r0,  get x/v from state vector & call correct subroutine
x0 = zeros(eltype(state0),3)
v0 = zeros(eltype(state0),3)
for k=1:3
  x0[k]=state0[k+1]
  v0[k]=state0[k+4]
end
s0=state0[11]
iter = kep_drift_ell_hyp!(x0,v0,gm,h,s0,state,jacobian,drift_first,d)
return
end


text:
using Catlab, Catlab.Theories
using Catlab.CategoricalAlgebra
using Catlab.Graphs
using Catlab.Graphics
using Catlab.Graphics.Graphviz
#import Catlab.Graphics.Graphviz: to_graphviz, to_graphviz_property_graph

using Colors
draw(g) = to_graphviz(g, node_labels=true, edge_labels=true)

GraphvizGraphs.to_graphviz(f::ACSetTransformation; kw...) =
  to_graphviz(GraphvizGraphs.to_graphviz_property_graph(f; kw...))

function GraphvizGraphs.to_graphviz_property_graph(f::ACSetTransformation; kw...)
  pg = GraphvizGraphs.to_graphviz_property_graph(dom(f); kw...)
  vcolors = hex.(range(colorant"#0021A5", stop=colorant"#FA4616", length=nparts(codom(f), :V)))
  ecolors = hex.(range(colorant"#6C9AC3", stop=colorant"#E28F41", length=nparts(codom(f), :E)))
  hex.(colormap("Oranges", nparts(codom(f), :V)))
  for v in vertices(dom(f))
    fv = f[:V](v)
    set_vprops!(pg, v, Dict(:color => "#$(vcolors[fv])"))
  end
  for e in edges(dom(f))
    fe = f[:E](e)
    set_eprops!(pg, e, Dict(:color => "#$(ecolors[fe])"))
  end
  pg
end

to_graphviz(Catlab.Graphs.BasicGraphs.TheoryGraph)

# # The Category of Graphs
#
# The Theory of Graphs is given by the following Schema:
# ```julia
# @present TheoryGraph(FreeSchema) begin
#   V::Ob
#   E::Ob
#   src::Hom(E,V)
#   tgt::Hom(E,V)
# end

# """ Abstract type for graphs, aka directed multigraphs.
# """
# @abstract_acset_type AbstractGraph <: HasGraph

# """ A graph, also known as a directed multigraph.
# """
# @acset_type Graph(TheoryGraph, index=[:src,:tgt]) <: AbstractGraph
# ```
#
# That is all we need to do to generate the functor category [TheoryGraph, FinSet].
# Catlab knows how to take a finitely presented category and generate all the data structures
# that you need to represent functors into FinSet and natural transformations between those functors.
# Note: the index=[:src, :tgt] keyword argument tells Catlab that you want to have an efficient index
# the preimages of those morphisms. in this example, we want to be able to find the incoming and
# outgoing edges of a vertex in O(1) time.

# ## Creating some Graphs

# Once you have fixed the schema (aka indexing category or theory), you can make some instances.
# Catlab has a DSL for specifying instances of any schema. It is called `@acset`.
# In order to specify a Functor F=(F₀, F₁) into FinSet, you need to provide some data.
#     1. For every A:Ob(C), you need F₀(A):FinSet
#     2. For every f:A→B, you need to specify a FinFunction F₁(f):F₀(A)→F₀(B)
# If the theory C has some equations, the data you provide would have to also satisfy those equations.
# The theory of graphs has no equations, so there are no constraints on the data you provide,
# except for those that come from functoriality.

e = @acset Graphs.Graph begin
    V = 2
    E = 1
    src = [1]
    tgt = [2]
end

draw(e)

# a wedge is two edges that share a target
w = @acset Graphs.Graph begin
    V = 3
    E = 2
    src=[1,3]
    tgt=[2,2]
end

draw(w)

# The CSet API generalizes the traditional Graph API

parts(w, :V)  # vertex set

parts(w,:E) # edge set

w[:src] # source map

w[:tgt] # target map

incident(w, 1, :src) # edges out of vertex 1

incident(w, 2, :tgt) # edges into vertex 2

w[incident(w, 2, :tgt), :src] # vertices that are the source of edges whose target is vertex 2

w[incident(w, 1, :src), :tgt] # vertices that are the target of edges whose src is vertex 1

# ### Exercise:
# 1. Use the @acset macro to make a graph with at least 5 vertices
# 2. Draw the graph
# 3. Compute in neighbors and out neighbors and make sure they match your expectations.
# 4. Write a function that computes the 2-hop out-neighbors of a vertex.

# # Graph Homomorphisms
# We can construct some graph homomorphisms between our two graphs.
# What data do we need to specify?

ϕ = ACSetTransformation(e,w,E=[1], V=[1,2])

is_natural(ϕ)

# The ACSetTransformation constructor does not automatically validate that the naturality squares commute!

ϕᵦ = ACSetTransformation(e,w,E=[1], V=[3,2])
is_natural(ϕᵦ)

# Our ϕᵦ in't natural because the edge map e₁ ↦ e₁ is not consistent with our vertex map, which sends v₁ ↦ v₃ and v₂ ↦ v₂. We can fix this by sending e₁ to e₂

ϕᵦ′ = ACSetTransformation(e,w,E=[2], V=[3,2])
is_natural(ϕᵦ′)

# So how does Catlab store the data of the natural transformation?
# the domain
ϕ.dom
# the codomain
ϕ.codom
# the components
ϕ.components
# you can see the components using standard indexing notation with the object names. Notice that while CSets are indexed by morphisms, natural transformations are indexed by objects.
ϕ[:V]
ϕ[:E]
# We can check the  naturality squares ourselves
# The sources are preserved: `src ⋅ ϕᵥ == ϕₑ ⋅ src`
ϕ[:V](dom(ϕ)[:,:src]) == codom(ϕ)[collect(ϕ[:E]), :src]
# The targets are preserved: `tgt ⋅ ϕᵥ == ϕₑ ⋅ tgt`
ϕ[:V](dom(ϕ)[:,:tgt]) == codom(ϕ)[collect(ϕ[:E]), :tgt]

# This approach generalizes to the following:
#
# ```julia
# function is_natural(α::ACSetTransformation{S}) where {S}
#    X, Y = dom(α), codom(α)
#    for (f, c, d) in zip(hom(S), dom(S), codom(S))
#      Xf, Yf, α_c, α_d = subpart(X,f), subpart(Y,f), α[c], α[d]
#      all(Yf[α_c(i)] == α_d(Xf[i]) for i in eachindex(Xf)) || return false
#    end
#    return true
# end
# ```
#
# Notice how we iterate over the homs in the schema category S `(f, c, d) in zip(hom(S), dom(S), codom(S))` We get one universally quantified equation `all(Yf[α_c(i)] == α_d(Xf[i]) for i in eachindex(Xf))` for each morphism in the indexing category
#
# ### Exercise:
# 1. Take your graph from the previous exercise and construct a graph homomorphism from the wedge (w) into it.
# 2. Check that the naturality equations are satisfied.
# 3. Explain why we don't need to specify any data for the source and target morphisms in TheoryGraph when definining a graph homomorphism

# ## Finding Homomorphisms Automatically
# As you saw in the previous exercise, constructing a natural transformation can be quite tedious. We want computers to automate tedious things for us. So we use an algorithm to enumerate all the homomorphisms between two CSets.

# CSet homomorphisms f:A→B are ways of finding a part of B that is shaped like A. You can view this as pattern matching. The graph A is the pattern and the graph B is the data. A morphism f:A→B is a collection of vertices and edges in B that is shaped like A. Note that you can ask Catlab to enforce constraints on the homomorphisms it will find including computing monic (injective) morphisms by passing the keyword `monic=true`. A monic morphism into B is a subobject of B.  You can pass `iso=true` to get isomorphisms.

t = @acset Graphs.Graph begin
    V = 3
    E = 3
    src = [1,2,1]
    tgt = [2,3,3]
end

draw(t)

T = @acset Graphs.Graph begin
    V = 6
    E = 9
    src = [1,2,1, 3, 1,5,2,2,4]
    tgt = [2,3,3, 4, 4,6,5,6,6]
end

draw(T)

# The simplest pattern in a graph is just a single edge and each homomorphism ϕ:e→G is a single edge in G.
length(homomorphisms(e, T, monic=true)) == nparts(T,:E) # number of edges

length(homomorphisms(t, T))

# We can define this helper function to print out all the homomorphisms between graphs. Because all our graphs are simple, we only need to look at the vertex components.
graphhoms(g,h) = begin
    map(homomorphisms(g,h)) do ϕ
        collect(ϕ[:V])
    end
end

graphhoms(t, T)

# Homs ϕ:t→T are little triangles in T, but homs ϕ:T→t are colorings of T with 3 colors. The fact that there are edges in t that are missing, means that it provides constraints on what graphs have morphisms into it. For example, there are no morphisms T→t.

graphhoms(T, t)

# The reason we don't have a morphism into t is vertices 2,3,4,5 aren't arranged into a triangle. We can relax those constraints by adding loops to the codomain. Loops in the codomain allow you to merge adjacent vertices when you construct the homomorphism.

add_loops!(g) = add_parts!(g, :E, nparts(g,:V), src=parts(g,:V), tgt=parts(g,:V))
add_loops(g) = begin
    h = copy(g)
    add_loops!(h)
    return h
end

draw(add_loops(t))

# Once we add loops, then we have morphisms.

length(homomorphisms(T,add_loops(t)))

# ## Bipartite Graphs
# Many computer science problems involve graphs that have two types of vertices. For example, when matching students to classes, you might represent the students as one type of vertex and the classes as another type of vertex. Then the edges (s,c) represent "student s is enrolled in class c". In this scenario there can be no edges from a class vertex to another class vertex, or from a student vertex to a student vertex. Graphs for which there exists such a classification are called bipartite graphs. In Category Theory, we shift from thinking about graphs with properties to graph homomorphisms that witness that property and think of bipartitioned graphs.

# First we construct a bipartite graph:
sq = apex(product(add_loops(e), add_loops(e)))
rem_parts!(sq, :E, [1,5,6,8,9])
draw(sq)

# We will use the symmetric edge graph to identify the bipartitions of this graph.
esym = @acset Graphs.Graph begin
    V = 2
    E = 2
    src = [1,2]
    tgt = [2,1]
end

draw(id(esym))
# There are two ways to bipartition sq.
graphhoms(sq, esym)

# This comes from the fact that esym has 2 automorphisms!
graphhoms(esym, esym)

# the first coloring
draw(homomorphisms(sq, esym)[1])

# but we can also swap the roles of the colors
draw(homomorphisms(sq, esym)[2])

# ### Exercise:
# 1. Construct a graph representation of a checkerboard
# 2. Draw the two bipartitions of the checkerboard

# We can generalize the notion of Bipartite graph to any number of parts. I like to call Kₖ the k-coloring classifier because homomorphims into α:G → Kₖ are k-colorings of G.

clique(k::Int) = begin
    Kₖ = Graphs.Graph(k)
    for i in 1:k
        for j in 1:k
            if j ≠ i
                add_parts!(Kₖ, :E, 1, src=i, tgt=j)
            end
        end
    end
    return Kₖ
end

K₃ = clique(3)
draw(id(K₃))

# Our graph T is not 2-colorable,
length(homomorphisms(T, esym))

# but we can use 3 colors to color T.
draw(homomorphism(T, K₃))

# ### Exercise:
# 1. Find a graph that is not 3-colorable
# 2. Find a graph that is not 4-colorable

# ## Homomorphisms in [C, Set] are like Types
# Any graph can play the role of the codomain. If you pick a graph that is incomplete, you get a more constrained notion of coloring where there are color combinations that are forbidden.
triloop = @acset Graphs.Graph begin
    V = 3
    E = 3
    src = [1,2,3]
    tgt = [2,3,1]
end

draw(id(triloop))

# With this graph, we can pick out directed 3-cycles in a graph like T2,
T2 = @acset Graphs.Graph begin
    V = 6
    E = 6
    src = [1,2,3,4,5,6]
    tgt = [2,3,1,5,6,4]
end
graphhoms(T2, triloop)

# and we can draw those cyclic roles with colors
draw(homomorphisms(T2, triloop)[1])

T3 = @acset Graphs.Graph begin
    V = 6
    E = 7
    src = [1,2,3,4,5,6, 2]
    tgt = [2,3,1,5,6,4, 4]
end
graphhoms(T3, triloop)

# Using the colors as shown:
draw(id(triloop))

# We can see our coloring of T3:
draw(homomorphisms(T3, triloop)[1])


text:
#############################################################################
# Joulia
# A Large-Scale Spatial Power System Model for Julia
# See https://github.com/JuliaEnergy/Joulia.jl
#############################################################################
# Example: ELMOD-DE

using Joulia

using DataFrames, CSV
using Gurobi

# data load for 2015 sample data
# see http://doi.org/10.5281/zenodo.1044463
pp_df = CSV.read("data_2015/power_plants.csv")
avail_con_df = CSV.read("data_2015/avail_con.csv")
prices_df = CSV.read("data_2015/prices.csv")

storages_df = CSV.read("data_2015/storages.csv")

lines_df = CSV.read("data_2015/lines.csv")

load_df = CSV.read("data_2015/load.csv")
nodes_df = CSV.read("data_2015/nodes.csv")
exchange_df = CSV.read("data_2015/exchange.csv")

res_df = CSV.read("data_2015/res.csv")
avail_pv = CSV.read("data_2015/avail_pv.csv")
avail_windon = CSV.read("data_2015/avail_windon.csv")
avail_windoff = CSV.read("data_2015/avail_windoff.csv")

avail = Dict(:PV => avail_pv,
	:WindOnshore => avail_windon,
	:WindOffshore => avail_windoff,
	:global => avail_con_df)

# generation of Joulia data types
pp = PowerPlants(pp_df, avail=avail_con_df, prices=prices_df)
storages = Storages(storages_df)
lines = Lines(lines_df)
nodes = Nodes(nodes_df, load_df, exchange_df)
res = RenewableEnergySource(res_df, avail)

# generation of the Joulia model
elmod = JouliaModel(pp, res, storages, nodes, lines)

# sclicing the data in weeks with the first full week starting at hour 49
slices = week_slices(49)

# running the Joulia model for week 30 using the Gurobi solver
results = run_model(elmod, slices[30], solver=GurobiSolver())

max_stars_repo_name: JuliaEnergy/JouliaExamples

max_issues_repo_name: JuliaEnergy/JouliaExamples

max_forks_repo_name: JuliaEnergy/JouliaExamples


text:
n = rand(1:5)
A = matrixdepot("wathen", n)

@test typeof(A) <: SparseMatrixCSC
@test issym(full(A))
@test isposdef(full(A))

println("'wathen' passed test...")


text:
using FluxGAN
using MLDatasets
using Flux
using BSON: @save

const n = parse(Int, ARGS[1])
const m = parse(Int, ARGS[2])
const skip = parse(Int, ARGS[3])

images = MNIST.traintensor(Float32)

const img_size = size(images)[1:end-1]
const img_num  = size(images)[end]
const img_dim  = *(img_size...)

# centering pixels, mapping to [-1,1]; and flattening image tensors
images = reshape((@. 2f0 * images - 1f0), (img_dim, img_num))

const d = 100

generator = Chain(
    Dense(d, 1200, relu),
    Dense(1200, 1200, relu),
    Dense(1200, img_dim, tanh_fast)
)

discriminator = Chain(
    Maxout(() -> Dense(img_dim, 240), 5),
    Maxout(() -> Dense(240, 240), 5),
    Dense(240, 1)
)

# learning rate
const η = 0.0002f0

# minibatch size

model = GAN(generator, discriminator;
    minibatch=m,
    η_gen=η,
    η_dscr=η,
    latent_dim=d,
    img_size=img_size
)

train!(model, images, iterations=n, skip=skip)

# uncomment below to save model when done training
# @save "models/mnist_goodfellow_MLP_test.bson" model

println("\nsaving generated images...\n")

output_dir = "images/MNIST"

info = ["MLP", "n", n, "m", m]

layout = (x=5, y=5)

image_grid(model, output_dir,
    layout=layout,
    file_info=info,
    img_res=150
)

println("finished!\n")

max_stars_repo_name: aarontrowbridge/GAN

max_issues_repo_name: aarontrowbridge/GAN

max_forks_repo_name: aarontrowbridge/GAN


text:
using Distributed
using Dates
@everywhere using Serialization
@everywhere using DelimitedFiles
@everywhere using JLD2

const jobs = RemoteChannel(()->Channel{Tuple}(32))
const res = RemoteChannel(()->Channel{Tuple}(32))

# count the number of jobs
num = 0
function make_jobs(ns::Array{Int, 1}, σs::Array{Int, 1})
    for n in ns
        for σ in σs
            put!(jobs, (n, σ))
        end
    end
    global num = length(ns) * length(σs)
end

@everywhere function do_work(jobs, res, folder, resfolder, nrep; adaptive = false)
    while true
        n, σ = take!(jobs)
        exec_time = @elapsed begin
            # no sampling rate
            @load folder * "/$(nrep)_XM_$(n)_$(σ).jld2" X M
            acc = cum_run_KTHSE(X, M, σ = 1.0*σ)
        end
        writedlm("$(resfolder)/$(nrep)_exectime_$(n)_$(σ).txt", exec_time)
        writedlm("$(resfolder)/$(nrep)_acc_$(n)_$(σ).txt", acc)
        put!(res, ("$(n)_$(σ)", exec_time, myid()))
    end
end

# main function
nrep, parent_folder, folder = ARGS
@async make_jobs(collect(15:5:50),
                 collect(1:4))
# timestamp = Dates.now()
# folder = "data_$(nf)_$(nr)_$(zoom)_$(timestamp)"
# mkdir(folder)

for p in workers()
    remote_do(do_work, p, jobs, res, parent_folder, folder, nrep; adaptive = true)
end

while num > 0
    job_id, exec_time, where = take!(res)
    println("$job_id finished in $(round(exec_time; digits=2)) secs on worker $where")
    global num -= 1
end

max_stars_repo_name: szcf-weiya/BaxterAlgorithms

max_issues_repo_name: szcf-weiya/BaxterAlgorithms

max_forks_repo_name: szcf-weiya/BaxterAlgorithms


text:
include("Base.jl")
include("Dict.jl")
include("JSON.jl")
include("Series.jl")
include("Show.jl")


text:
#
#  TIMESTEP
#

"""
	gettime(ts::FixedTimestep)

Return the time (year) represented by Timestep `ts`
"""
function gettime(ts::FixedTimestep{FIRST, STEP, LAST}) where {FIRST, STEP, LAST}
	return FIRST + (ts.t - 1) * STEP
end

"""
	gettime(ts::VariableTimestep)

Return the time (year) represented by Timestep `ts`
"""
function gettime(ts::VariableTimestep)
	return ts.current
end

"""
	is_time(ts::AbstractTimestep, t::Int)

Return true or false, true if the current time (year) for `ts` is `t`
"""
function is_time(ts::AbstractTimestep, t::Int)
	return gettime(ts) == t
end

"""
	is_first(ts::AbstractTimestep)

Return true or false, true if `ts` is the first timestep to be run.
"""
function is_first(ts::AbstractTimestep)
	return ts.t == 1
end

"""
	is_timestep(ts::AbstractTimestep, t::Int)

Return true or false, true if `ts` timestep is step `t`.
"""
function is_timestep(ts::AbstractTimestep, t::Int)
	return ts.t == t
end

"""
	is_last(ts::FixedTimestep)

Return true or false, true if `ts` is the last timestep to be run.
"""
function is_last(ts::FixedTimestep{FIRST, STEP, LAST}) where {FIRST, STEP, LAST}
	return gettime(ts) == LAST
end

"""
	is_last(ts::VariableTimestep)

Return true or false, true if `ts` is the last timestep to be run.  Note that you may
run `next_timestep` on `ts`, as ths final timestep has not been run through yet.
"""
function is_last(ts::VariableTimestep{TIMES}) where {TIMES}
	return gettime(ts) == TIMES[end]
end

function finished(ts::FixedTimestep{FIRST, STEP, LAST}) where {FIRST, STEP, LAST}
	return gettime(ts) > LAST
end

function finished(ts::VariableTimestep{TIMES}) where {TIMES}
	return gettime(ts) > TIMES[end]
end

function next_timestep(ts::FixedTimestep{FIRST, STEP, LAST}) where {FIRST, STEP, LAST}
	if finished(ts)
		error("Cannot get next timestep, this is last timestep.")
	end
	return FixedTimestep{FIRST, STEP, LAST}(ts.t + 1)
end

function next_timestep(ts::VariableTimestep{TIMES}) where {TIMES}
	if finished(ts)
		error("Cannot get next timestep, this is last timestep.")
	end
	return VariableTimestep{TIMES}(ts.t + 1)
end

function Base.:-(ts::FixedTimestep{FIRST, STEP, LAST}, val::Int) where {FIRST, STEP, LAST}
	if val != 0 && is_first(ts)
		error("Cannot get previous timestep, this is first timestep.")
	elseif ts.t - val <= 0
		error("Cannot get requested timestep, precedes first timestep.")
	end
	return FixedTimestep{FIRST, STEP, LAST}(ts.t - val)
end

function Base.:-(ts::VariableTimestep{TIMES}, val::Int) where {TIMES}
	if val != 0 && is_first(ts)
		error("Cannot get previous timestep, this is first timestep.")
	elseif ts.t - val <= 0
		error("Cannot get requested timestep, precedes first timestep.")
	end
	return VariableTimestep{TIMES}(ts.t - val)
end

function Base.:-(ts::TimestepValue, val::Int)
	return TimestepValue(ts.value; offset = ts.offset - val)
end

function Base.:-(ts::TimestepIndex, val::Int)
	return TimestepIndex(ts.index - val)
end

function Base.:+(ts::FixedTimestep{FIRST, STEP, LAST}, val::Int) where {FIRST, STEP, LAST}
	if finished(ts)
		error("Cannot get next timestep, this is last timestep.")
	elseif gettime(ts) + val > LAST + 1
		error("Cannot get requested timestep, exceeds last timestep.")
	end
	return FixedTimestep{FIRST, STEP, LAST}(ts.t + val)
end

function Base.:+(ts::VariableTimestep{TIMES}, val::Int) where {TIMES}
	if finished(ts)
		error("Cannot get next timestep, this is last timestep.")
	elseif gettime(ts) + val > TIMES[end] + 1
		error("Cannot get requested timestep, exceeds last timestep.")
	end
	new_ts = VariableTimestep{TIMES}(ts.t + val)
end

function Base.:+(ts::TimestepValue, val::Int)
	return TimestepValue(ts.value; offset = ts.offset + val)
end

function Base.:+(ts::TimestepIndex, val::Int)
	return TimestepIndex(ts.index + val)
end

# Colon support
function Base.:(:)(start::T, step::Int, stop::T) where {T<:TimestepIndex}
	indices = [start.index:step:stop.index...]
	return TimestepIndex.(indices)
end

function Base.:(:)(start::T, stop::T) where {T<:TimestepIndex}
	return Base.:(:)(start, 1, stop)
end

#
#  CLOCK
#

function Clock(time_keys::Vector{Int})
    last = time_keys[end]

    if isuniform(time_keys)
        first, stepsize = first_and_step(time_keys)
        return Clock{FixedTimestep}(first, stepsize, last)
    else
        last_index = findfirst(isequal(last), time_keys)
        times = (time_keys[1:last_index]...,)
        return Clock{VariableTimestep}(times)
    end
end

function timestep(c::Clock)
	return c.ts
end

function time_index(c::Clock)
	return c.ts.t
end

"""
	gettime(c::Clock)

Return the current time of the timestep held by the `c` clock.
"""
function gettime(c::Clock)
	return gettime(c.ts)
end

function advance(c::Clock)
	c.ts = next_timestep(c.ts)
	nothing
end

function finished(c::Clock)
	return finished(c.ts)
end

function Base.reset(c::Clock)
	c.ts = c.ts - (c.ts.t - 1)
	nothing
end

function timesteps(c::Clock)
	c = deepcopy(c)
	timesteps = AbstractTimestep[]
	while !finished(c)
		push!(timesteps, timestep(c))
		advance(c)
	end
	return timesteps
end


text:
using Documenter
using BioStructures

makedocs(
    sitename = "BioStructures.jl",
    pages = [
        "Home"         => "index.md",
        "Documentation"=> "documentation.md",
        "Examples"     => "examples.md",
        "API"          => "api.md",
    ]
)

deploydocs(repo="github.com/BioJulia/BioStructures.jl.git")


text:
#=
P36 (**) Determine the prime factors of a given positive integer (2).
Construct a list containing the prime factors and their multiplicity.
Example:
* (prime-factors-mult 315)
((3 2) (5 1) (7 1))
Hint: The problem is similar to problem P13.
=#

SMALLEST_PRIME = 2

# From P35
function next_prime(n)
    # Brute-force incremental check
    isprime(n) ? n : next_prime(n + 1)
end

# From P35
function _prime_factors(n, x, pf)
    if isprime(n)
        vcat(pf, n)
    else
        if mod(n, x) == 0
            _prime_factors(div(n, x), SMALLEST_PRIME, vcat(pf, x))
        else
            _prime_factors(n, next_prime(x + 1), pf)
        end
    end
end

# From P35
function prime_factors(n)
    # Starting with the smallest prime, 2
    isprime(n) ? n : _prime_factors(n, SMALLEST_PRIME, [])
end

function _prime_factors_multi(pf, n, count, counters)
    if isempty(pf)
        vcat(counters, Vector{}[[n, count]])
    else
        if pf[1] == n
            _prime_factors_multi(pf[2:end], n, count + 1, counters)
        else
            new_counters = vcat(counters, Vector{}[[n, count]])
            _prime_factors_multi(pf[2:end], pf[1], 1, new_counters)
        end
    end
end

function prime_factors_multi(n)
    pf = prime_factors(n)
    _prime_factors_multi(pf[2:end], pf[1], 1, Vector{}[])
end

@assert prime_factors_multi(315) == Vector{}[[3, 2], [5, 1], [7, 1]]

println("Tests passed: JL-36.jl")


text:
# Index of entries in D matrix (sign included)
struct D_Index
	sign::Int64
	μ::Int64
	m::Int64
end

# Equation (1.1) - forms the covariant matrix D(Q)(indices only)
function rotation_D_matrix(L::Integer)
	if L<0
		error("Orbital type shall be represented as a positive integer!")
	end
    D = Array{D_Index}(undef, 2 * L + 1, 2 * L + 1)
    for i = 1 : 2 * L + 1
        for j = 1 : 2 * L + 1
            D[j,i] = D_Index(1, i - 1 - L, j - 1 - L);
        end
    end
	return D
end

# Equation (1.1) - forms the covariant matrix D(Q)(indices only)
function rotation_D_matrix_ast(L::Integer)
	if L<0
		error("Orbital type shall be represented as a positive integer!")
	end
    D = Array{D_Index}(undef, 2 * L + 1, 2 * L + 1)
    for i = 1 : 2 * L + 1
        for j = 1 : 2 * L + 1
            D[i,j] = D_Index((-1)^(i+j), -(i - 1 - L), -(j - 1 - L));
        end
    end
	return D
end

# Equation (1.2) - vector value coupling coefficients
# ∫_{SO3} D^{ll}_{μμmm} D^*(Q) e^t dQ -> 2L+1 column vector
function vec_cou_coe(rotc::Rot3DCoeffs{T},
					   ll::StaticVector{N},
	                   mm::StaticVector{N},
					   μμ::StaticVector{N},
					   L::Integer, t::Integer) where {T,N}
	if t > 2L + 1 || t <= 0
		error("Rotation D matrix has no such column!")
	end
	Z = zeros(2L + 1)
	D = rotation_D_matrix_ast(L)
	Dt = D[:,t]   # D^* ⋅ e^t
	μt = [Dt[i].μ for i in 1:2L+1]
	mt = [Dt[i].m for i in 1:2L+1]
	LL = [ll; L]
	for i = 1:(2L + 1)
		MM = [mm; mt[i]]
		KK = [μμ; μt[i]]
		Z[i] = Dt[i].sign * rotc(LL, MM, KK)
	end
	return Z
end

# Equation (1.5) - possible set of mm w.r.t. index ll & vector k
function collect_m(ll::StaticVector{N}, k::T) where {N,T}
	A = CartesianIndices(ntuple(i -> -ll[i]:ll[i], length(ll)));
	B = Array{typeof(A[1].I)}(undef, 1, prod(size(A)))
	t = 0
	for i in A
		if prod(sum(i.I) .+ k) == 0
			t = t + 1
			B[t] = i.I
		end
	end
	B = [SVector(i) for i in B[1:t]]
	return B
end

function gramian(A::Rot3DCoeffs{T}, ll::StaticVector{N},
	                 L::Integer) where {T,N}
	LenM = 0
	D = rotation_D_matrix_ast(L)
	D1 = D[:,1]
	μt = [D1[i].μ for i in 1:2L+1]
	μ_list = collect_m(ll,μt)
	Z = fill(zeros(2L + 1), (length(μ_list), length(μ_list)))
	for t = 1:2L+1
		Dt = D[:,t]
		mt = [Dt[i].m for i in 1:2L+1]
		m_list = collect_m(ll,mt)
		for (im, mm) in enumerate(m_list), (iμ, μμ) in enumerate(μ_list)
			Z[iμ, im + LenM] = vec_cou_coe(A, ll, mm, μμ, L, t)
		end
		LenM += length(m_list)
	end
	return Z' * Z, Z, μ_list
end

function rc_basis_all(A::Rot3DCoeffs{T}, ll::StaticVector{N},
	                  L::Integer) where {N,T}
	G, C, μ_list = gramian(A, ll, L)
	S = svd(G)
	rk = rank(G; rtol =  1e-8)
	Urcpi = fill(zeros(2L + 1), (rk, length(μ_list)))
	U = S.U[:, 1:rk]
	Sigma = S.S[1:rk]
	Urcpi = C * U * Diagonal(sqrt.(Sigma))^(-1)
	return Urcpi', μ_list
end

# Equation (1.12) - Gramian over nn
function gramian_nn(A::Rot3DCoeffs,
			     nn::StaticVector{N},
			     ll::StaticVector{N},
				 L::Integer) where {N}
	Uri, μ_list = rc_basis_all(A, ll, L)
	G = zeros(size(Uri)[1], size(Uri)[1])
	for σ in permutations(1:N)
       if (nn[σ] != nn) || (ll[σ] != ll); continue; end
       for (iU1, mm1) in enumerate(μ_list), (iU2, mm2) in enumerate(μ_list)
          if mm1[σ] == mm2
             for i1 = 1:size(Uri)[1]
				 for i2 = 1:size(Uri)[1]
                 	G[i1, i2] += Uri[i1, iU1] * Uri[i2, iU2]'
				end
             end
          end
       end
    end
    return G, Uri, μ_list
end

## Equation (1.13) - LI coefficients(& corresponding μ) over nn, ll
function yvec_symm_basis(A::Rot3DCoeffs,
				         nn::StaticVector{N},
					     ll::StaticVector{N},
					      L::Integer) where {N}
	G, C, μ_list= gramian_nn(A, nn, ll, L)
	S = svd(G)
	rk = rank(G; rtol =  1e-8)
	Urcpi = fill(zeros(2L + 1), (rk, length(μ_list)))
	U = S.U[:, 1:rk]
	Sigma = S.S[1:rk]
	Urcpi = C' * U * Diagonal(sqrt.(Sigma))
	return Urcpi', μ_list
end


text:
using Yao, YaoExtensions

@const_gate sqrtX = ComplexF64[0.5+0.5im 0.5-0.5im; 0.5-0.5im 0.5+0.5im]
@const_gate sqrtY = ComplexF64[0.5+0.5im -0.5-0.5im; 0.5+0.5im 0.5+0.5im]
@const_gate sqrtW = ComplexF64[0.5+0.5im -sqrt(1/2)*im; sqrt(1/2) 0.5+0.5im]

function to_yao_block(layout::RQCLayout{VT, PT, SCGate}, s = collect(vertices(layout));
        haar = false) where {VT, PT}
    vs = sort!(intersect(collect(vertices(layout)), s))
    qubit_map = gen_qubit_map(layout, s)
    N = length(vs)
    D = depth(layout)
    blk = chain(N)
    single_prev = [I2 for i = 1:N]
    for d in 1:D
        for v in vs
            q = qubit_map[v]
            if gates(layout, v)[d] === Single
                if haar
                    M = rand(Haar(2), 2)
                    g = matblock(M)
                else
                    single_prev[q] == I2 && (g2 = rand((sqrtX, sqrtY, sqrtW)))
                    single_prev[q] == sqrtX && (g2 = rand((sqrtY, sqrtW)))
                    single_prev[q] == sqrtY && (g2 = rand((sqrtX, sqrtW)))
                    single_prev[q] == sqrtW && (g2 = rand((sqrtX, sqrtY)))
                    g = chain(Rz(rand()*2π), g2, Rz(rand()*2π))
                end
                push!(blk, put(q => g))
            elseif gates(layout, v)[d] === Noise
                g = rand((I2, X, Y, Z))
                push!(blk, put(q => g))
            elseif gates(layout, v)[d] === FSim
                u = partner(layout, v, depth_to_cycle(d))
                u in vs || continue
                q1, q2 = q, qubit_map[u]
                push!(blk, put((q1, q2) => FSimGate(π/2, π/6)))
            end
        end
    end
    return blk
end


text:
module HP

using DocStringExtensions

import ..Delayed

export Prob

"""
$(TYPEDEF)
$(TYPEDFIELDS)

State the weighted probability of an option, for use with PChoice.
"""
struct Prob
    probability::Float64
    option
end

struct PChoice <: Delayed.AbstractSwitch
    choice::Delayed.Param
    options::Vector

    @doc """
    $(TYPEDSIGNATURES)

    Choose from a list of options with weighted probabilities

    # Arguments
    - `label`               : Label
    - `probability_options` : Array of Prob objects

    Example:

    ```julia
    HP.PChoice(:a, [
        Prob(0.1, 0),
        Prob(0.2, 1),
        Prob(0.7, 2),
    ])
    ```

    Note that the Prob probability weights must sum to 1.
    """
    function PChoice(label::Symbol, probability_options::Vector{Prob})
        if isempty(probability_options)
            throw(ArgumentError("$(label): you must give at least one choice"))
        end
        probabilities = [o.probability for o in probability_options]
        if sum(probabilities) != 1
            throw(ArgumentError("Prob probability weights must sum to 1"))
        end
        options = [o.option for o in probability_options]

        return new(Delayed.Param(label, Delayed.CategoricalIndex(probabilities)), options)
    end
end

struct Choice <: Delayed.AbstractSwitch
    choice::Delayed.Param
    options::Vector

    @doc """
    $(TYPEDSIGNATURES)

    Randomly choose which option will be extracted, and also supply the list of options.
    The elements of options can themselves be [nested] stochastic expressions. In this case, the
    stochastic choices that only appear in some of the options become conditional parameters.

    Example:

    ```julia
    Dict(
        :example => HP.Choice(:example, [1.0, 0.9, 0.8, 0.7]),
    )
    ```

    Example of conditional paramaters:

    ```julia
    Dict(
        :example => HP.Choice(:example, [
            (:case1, HP.Uniform(:param1, 0.0, 1.0)),
            (:case2, HP.Uniform(:param2, -10.0, 10.0)),
        ]),
    )
    ```

    ':param1' and ':param2' are examples of conditional parameters. Each of ':param1' and
    ':param2' only features in the returned sample for a particular value of ':example'. If
    ':example' is 0, then ':param1' is used but not ':param2'. If ':example' is 1, then
    ':param2' is used but not ':param1'.

    Example with nested arrays of different lengths:

    ```julia
    Dict(
        :example => HP.Choice(:example, [
            [HP.Normal(:d0_c0, 0.0, 5.0)],
            [HP.Normal(:d1_c0, 0.0, 5.0), HP.Normal(:d1_c1, 0.0, 5.0)],
            [
                HP.Normal(:d2_c0, 0.0, 5.0), HP.Normal(:d2_c1, 0.0, 5.0),
                HP.Normal(:d2_c2, 0.0, 5.0),
            ],
        ]),
    )
    ```

    Note that all labels (the symbol given as the first parameter to all the HP.* functions)
    must be unique. These labels identify the parts of the space that the optimiser learns from
    over iterations.
    """
    function Choice(label::Symbol, options::Vector)
        if isempty(options)
            throw(ArgumentError("$(label): you must give at least one choice"))
        end

        return new(Delayed.Param(label, Delayed.RandIndex(length(options))), options)
    end
end

struct Uniform <: Delayed.AbstractParam
    label::Symbol
    obj::Delayed.Uniform

    @doc """
    $(TYPEDSIGNATURES)

    Returns a value with uniform probability from between `low` and `high`.
    When optimising, this variable is constrained to a two-sided interval.

    ```julia
    Dict(
        :example => HP.Uniform(:example, 0.0, 1.0),
    )
    ```

    where `label` is the parameter and the returned value is uniformly distributed between
    `low` at 0.0 and `high` at 1.0
    """
    function Uniform(label::Symbol, low::Delayed.NestedFloat, high::Delayed.NestedFloat)
        return new(label, Delayed.Uniform(low, high))
    end
end

struct QuantUniform <: Delayed.AbstractParam
    label::Symbol
    obj::Delayed.QuantUniform

    @doc """
    $(TYPEDSIGNATURES)

    Returns a value uniformly between low and high, with a quantisation.
    When optimising, this variable is constrained to a two-sided interval.

    ```julia
    Dict(
        :example => HP.QuantUniform(:example, 0.0, 10.0, 2.0),
    )
    ```

    where `label` is the parameter and the returned value is uniformly distributed between
    `low` at 0.0 and `high` at 10.0, with the `q`uantisation set at 2.0.
    Valid sampled values would be 0.0, 2.0, 4.0, 6.0, 8.0 and 10.0.
    """
    function QuantUniform(
        label::Symbol, low::Delayed.NestedFloat, high::Delayed.NestedFloat,
        q::Delayed.NestedFloat
    )
        return new(label, Delayed.QuantUniform(low, high, q))
    end
end

struct Normal <: Delayed.AbstractParam
    label::Symbol
    obj::Delayed.Normal

    @doc """
    $(TYPEDSIGNATURES)

    Returns a real value that's normally-distributed with mean mu and standard deviation sigma.
    When optimizing, this is an unconstrained variable.

    ```julia
    Dict(
        :example => HP.Normal(:example, 4.0, 5.0),
    )
    ```
    """
    function Normal(label::Symbol, mu::Delayed.NestedFloat, sigma::Delayed.NestedFloat)
        return new(label, Delayed.Normal(mu, sigma))
    end
end

struct QuantNormal <: Delayed.AbstractParam
    label::Symbol
    obj::Delayed.QuantNormal

    @doc """
    $(TYPEDSIGNATURES)

    Returns a real value that's normally-distributed with mean mu and standard deviation sigma, with a quantisation.
    When optimizing, this is an unconstrained variable.

    ```julia
    Dict(
        :example => HP.QuantNormal(:example, 2., 0.5, 1.0),
    )
    ```

    In this example, the values are sampled normally first, and then `q`uantised in rounds of 1.0, so one
    only observes 1.0, 2.0, 3.0, etc, centered around 2.0.

    N.B. that due to rounding, the observed values will not follow exactly normal distribution, particularly
    if sigma is much smaller than quantisation.
    """
    function QuantNormal(
        label::Symbol, mu::Delayed.NestedFloat, sigma::Delayed.NestedFloat,
        q::Delayed.NestedFloat
    )
        return new(label, Delayed.QuantNormal(mu, sigma, q))
    end
end

struct LogNormal <: Delayed.AbstractParam
    label::Symbol
    obj::Delayed.LogNormal

    @doc """
    $(TYPEDSIGNATURES)

    Returns a value drawn according to exp(normal(mu, sigma)) so that the logarithm of the
    sampled value is normally distributed. When optimising, this variable is constrained to be
    positive.

    ```julia
    Dict(
        :example => HP.LogNormal(:example, log(3.0), 0.5),
    )
    ```

    In this example, the log normal distribution will be centred around 3. The distribution is
    not truncated.
    """
    function LogNormal(label::Symbol, mu::Delayed.NestedFloat, sigma::Delayed.NestedFloat)
        return new(label, Delayed.LogNormal(mu, sigma))
    end
end

struct LogQuantNormal <: Delayed.AbstractParam
    label::Symbol
    obj::Delayed.LogQuantNormal

    @doc """
    $(TYPEDSIGNATURES)

    Returns a value drawn according to exp(normal(mu, sigma)), with a quantisation, so that the
    logarithm of the sampled value is normally distributed. When optimising, this variable is
    constrained to be positive.

    ```julia
    Dict(
        :example => HP.LogQuantNormal(:example, log(1e-3), 0.5*log(10), log(sqrt(10))),
    )
    ```

    In this example, the log normal distribution will be centred around 1e-3, with stddev of sqrt(10) (in exp).
    The distribution is not truncated. The distinct values would therefore be in every power of `sqrt(10)`.
    """
    function LogQuantNormal(
        label::Symbol, mu::Delayed.NestedFloat, sigma::Delayed.NestedFloat,
        q::Delayed.NestedFloat
    )
        return new(label, Delayed.LogQuantNormal(mu, sigma, q))
    end
end

struct QuantLogNormal <: Delayed.AbstractParam
    label::Symbol
    obj::Delayed.QuantLogNormal

    @doc """
    $(TYPEDSIGNATURES)

    Returns a value drawn according to exp(normal(mu, sigma)), with a quantisation, so that the
    logarithm of the sampled value is normally distributed. When optimising, this variable is
    constrained to be positive.

    ```julia
    Dict(
        :example => HP.QuantLogNormal(:example, log(3.0), 0.5, 2.0),
    )
    ```

    In this example, the log normal distribution will be centred around 3. The distribution is
    not truncated. The values with be quantised to multiples of 2, i.e. 2.0, 4.0, 6.0, etc.
    """
    function QuantLogNormal(
        label::Symbol, mu::Delayed.NestedFloat, sigma::Delayed.NestedFloat,
        q::Delayed.NestedFloat
    )
        return new(label, Delayed.QuantLogNormal(mu, sigma, q))
    end
end

struct LogUniform <: Delayed.AbstractParam
    label::Symbol
    obj::Delayed.LogUniform

    @doc """
    $(TYPEDSIGNATURES)

    Returns a value drawn according to exp(uniform(low, high)) such that the logarithm
    of the return value is uniformly distributed.
    When optimizing, this variable is constrained to the interval [exp(low), exp(high)].

    ```julia
    Dict(
        :example => HP.LogUniform(:example, log(3.0), 1.0),
    )
    ```
    """
    function LogUniform(label::Symbol, low::Delayed.NestedFloat, high::Delayed.NestedFloat)
        return new(label, Delayed.LogUniform(low, high))
    end
end

struct LogQuantUniform <: Delayed.AbstractParam
    label::Symbol
    obj::Delayed.LogQuantUniform

    @doc """
    $(TYPEDSIGNATURES)

    Returns a value drawn according to exp(uniform(low, high)), quantised in log-space, such that
    the logarithm of the return value is uniformly distributed. The value is constrained to be positive.

    Suitable for searching logarithmically through a space while keeping the number of candidates
    bounded, e.g. searching learning rate through 1e-6 to 1e-1

    ```julia
    Dict(
        :example => HP.LogQuantUniform(:example, log(1e-6), log(1e-1), log(10)),
    )
    ```
    """
    function LogQuantUniform(
        label::Symbol, low::Delayed.NestedFloat, high::Delayed.NestedFloat,
        q::Delayed.NestedFloat
    )
        return new(label, Delayed.LogQuantUniform(low, high, q))
    end
end

struct QuantLogUniform <: Delayed.AbstractParam
    label::Symbol
    obj::Delayed.QuantLogUniform

    @doc """
    $(TYPEDSIGNATURES)

    Returns a value drawn according to exp(uniform(low, high)), with a quantisation q, such that
    the logarithm of the return value is uniformly distributed.

    Suitable for a discrete variable with respect to which the objective is "smooth" and gets
    smoother with the size of the value, but which should be bounded both above and below.

    ```julia
    Dict(
        :example => HP.QuantLogUniform(:example, log(3.0), log(10.0), 2.0),
    )
    ```
    """
    function QuantLogUniform(
        label::Symbol, low::Delayed.NestedFloat, high::Delayed.NestedFloat,
        q::Delayed.NestedFloat
    )
        return new(label, Delayed.QuantLogUniform(low, high, q))
    end
end

end # module HP


text:
include(joinpath(@__DIR__, "../src/DB.jl"))
include(joinpath(@__DIR__, "FLS.jl"))

db = init_db(reset=true)
@test add_fls(db, fls) == 1
@test get_fls(db, 1) == fls

@test add_fls(db, fls) == 1 # does not insert again

@test add_premodel(db, fls, xrel) == 1
@test get_premodel(db, 1) == (fls => xrel)


text:
function transformPolyline(value, index)
end

function decodePolyline(expr; precision=5)
end


text:
using ZWarpedArrays, Images, CoordinateTransformations
using Test

img0 = zeros(Normed{UInt16, 16}, 5,5,3,3)
imgpost = zeros(5,5,3,3)
normed1 = reinterpret(Normed{UInt16,16}, UInt16(2^8))
imgpost[3,3,:,:] .= normed1 #We hope the warped array looks like this
img0[3,3,1,:] .= normed1 #centered
img0[4,3,2,:] .= normed1 #shifted +1 in x
img0[4,4,3,:] .= normed1 #shifted +1 in x and y
tfms = [IdentityTransformation(); Translation(1,0); Translation(1,1)]
za  = ZWarpedArray(img0, tfms, Float32; correct_bias=false, sqrt_tfm=false);
@test eltype(za[:,:,2,2]) == Float32
za  = ZWarpedArray(img0, tfms, Float64; correct_bias=false, sqrt_tfm=false);
@test eltype(za[:,:,2,2]) == Float64

for t = 1:size(imgpost,4)
    for z = 1:size(imgpost,3)
        slw = za[:,:,z,t]
        slp = imgpost[:,:,z,t]
        @test !isnan(slw[3,3])
        for i in eachindex(slw)
            if !isnan(slw[i])
                @test isapprox(slw[i], slp[i])
            end
        end
    end
end

#test 3D (single stack)
img0t = img0[:,:,:,1]
imgpostt = imgpost[:,:,:,1]
za3  = ZWarpedArray(img0t, tfms, Float64; correct_bias=false, sqrt_tfm=false);

for z = 1:size(imgpostt,3)
    slw = za3[:,:,z]
    slp = imgpostt[:,:,z]
    @test !isnan(slw[3,3])
    for i in eachindex(slw)
        if !isnan(slw[i])
            @test isapprox(slw[i], slp[i])
        end
    end
end

##With bias correction
img0t .+= reinterpret(Normed{UInt16, 16}, UInt16(100))
za  = ZWarpedArray(img0t, tfms, Float64; correct_bias=true, sqrt_tfm=false);

for z = 1:size(imgpostt,3)
    slw = za[:,:,z]
    slp = imgpostt[:,:,z]
    @test !isnan(slw[3,3])
    for i in eachindex(slw)
        if !isnan(slw[i])
            @test isapprox(slw[i], slp[i])
        end
    end
end
##With bias correction and square root
za  = ZWarpedArray(img0t, tfms, Float64; correct_bias=true, sqrt_tfm=true);

for z = 1:size(imgpostt,3)
    slw = za[:,:,z]
    slp = imgpostt[:,:,z]
    @test !isnan(slw[3,3])
    for i in eachindex(slw)
        if !isnan(slw[i])
            @test isapprox(slw[i], sqrt(slp[i]))
        end
    end
end

meta = ImageMeta(img0t)
zm  = ZWarpedArray(meta, tfms, Float64; correct_bias=false, sqrt_tfm=false);
@test isa(zm, ImageMeta)
@test all(data(zm) .== za3)

aa = AxisArray(img0t)
zaa  = ZWarpedArray(aa, tfms, Float64; correct_bias=false, sqrt_tfm=false);
@test isa(zaa, AxisArray)
@test all(data(zaa) .== za3)


text:
ENV["GKSwstype"]="100"

using Test
using WaterWaves1D

param = ( ϵ  = 1/2, μ = 1)
paramX= ( N  = 2^8, L  = 10)
paramT= ( T  = 5, dt = 0.1)

init     = Init(x->exp.(-x.^2),x->0 .*x )
model    = WaterWaves(merge(param,paramX);verbose=false)
problem1 = Problem( model, init, merge(paramT,paramX); solver = RK4(paramX) )

# @testset "LoadSave" begin
#
#     dump  = convert( ProblemSave, problem1 )
#     pload = convert( Problem, dump )
#
#     save(problem1, "testsave")
#     pload = loadpb("testsave")
#
#     @test pload.model.kwargs == problem1.model.kwargs
#     @test pload.solver.Uhat == problem1.solver.Uhat
#
# end

@testset "Parameters" begin

    @test param.ϵ  == 0.5
    @test param.μ  == 1
    @test paramX.N  == 256
    @test paramX.L  == 10
    @test paramT.T  == 5
    @test paramT.dt == 0.1

end

solve!(problem1)

@testset "Test problem with water waves model" begin
    @test !any(isnan,problem1.data.U[end][1])
    @test !any(isnan,problem1.data.U[end][2])
end

WW2  = WWn(merge(param,paramX);δ=0.01,n=2,dealias=1,verbose=false)
problem2 = Problem(WW2, init, merge(paramX,paramT) )

solve!( problem2 )

@testset "Test problem with Pseudo-spectral model" begin
    @test !any(isnan,problem2.data.U[end][1])
    @test !any(isnan,problem2.data.U[end][2])
end


text:
include("odehelper.jl")

function basic!(u,p,t)
    u= 1.01*u
end

f(u,p,t) = 1.01*u
solve_plot_print(basic!, 1/2, (0.0,1.0))
f(u,p,t) = 10-2*u
solve_plot_print(f, 1, (0.0,1.0))

max_stars_repo_name: bkgsur/DifferentialEquationsWithJulia

max_issues_repo_name: bkgsur/DifferentialEquationsWithJulia

max_forks_repo_name: bkgsur/DifferentialEquationsWithJulia