text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Profiling time and memory for some of the functions in ACASXProblem and SA.
"""
module ACASX_SA_Timing

export acasx_sa_timing

using ExprSearch.SA
using DerivationTrees
using RLESUtils, SwapBuffers, Observers, Loggers, StringUtils

using GrammarExpts
using ACASXProblem

function acasx_sa_timing(problem=nothing, outdir::AbstractString="./ACASX_SA_Timing"; seed=1,
                  logfileroot::AbstractString="acasx_sa_log",

                  runtype::Symbol=:nmacs_vs_nonnmacs,
                  data::AbstractString="dasc",
                  manuals::AbstractString="dasc_manual",
                  clusterdataname::AbstractString="josh1",

                  maxsteps::Int64=20,
                  T1::Float64=12.184,
                  alpha::Float64=0.99976,
                  n_epochs::Int64=50,
                  n_starts::Int64=1,

                  loginterval::Int64=100,
                  vis::Bool=true,
                  observer::Observer=Observer())

  if problem == nothing
    problem = ACASXClustering(runtype, data, manuals, clusterdataname)
  end
  p = SAESParams(maxsteps, T1, alpha, n_epochs, n_starts, observer)

  grammar = get_grammar(problem)
  tree_params = DerivTreeParams(grammar, p.maxsteps)
  s_buffer = SwapBuffer(SA.SAState(DerivationTree(tree_params)),
                        SA.SAState(DerivationTree(tree_params)))
  result = SAESResult()

  j = 1
  @printeval @time s = active(s_buffer)
  @printeval @time SA.initialize!(s, problem, result) #initialize randomly

  T = p.T1
  i = 1
  @printeval @time s = active(s_buffer)
  @printeval @time sp = inactive(s_buffer)
  @printeval @time SA.neighbor!(sp, s, problem, result)

  @printeval @time iscomplete(sp.tree)
  @printeval @time SA.accept(s, sp, T)
  @printeval @time swap!(s_buffer)

  @printeval @time SA.initialize!(s, problem) #initialize randomly
  @printeval @time SA.update!(result, s)

  @printeval @time copy!(sp, s)
  @printeval @time SA.perturb!(sp.tree)

  @printeval @time rand!(s.tree, typemax(Int64))
  @printeval @time s.expr = get_expr(s.tree)
  @show s.expr
  @printeval @time s.fitness = get_fitness(problem, s.expr)
  @printeval @time f = to_function(problem, s.expr)
  @printeval @time predicts = map(f, problem.Dl.records)
  @printeval @time _, _, ent_post = ACASXProblem.entropy_metrics(predicts, problem.Dl.labels)

  @printeval @time predicts = Array(Bool, length(problem.Dl.records))
  @printeval @time begin
    records = problem.Dl.records
    for i = 1:length(records)
      predicts[i] = f(records[i])
    end
  end
  @printeval @time predicts[1] = f(records[1])

  v1 = rand(50)
  v2 = rand(50)
  @printeval @time ACASXProblem.abs_diff_lte(v1, v2, 0.1);
  @printeval @time ACASXProblem.sign_(v1,v2);

  members = rand(1:100, 30)
  @printeval @time Dl_sub = problem.Dl[members]

  problem
end

end #module


text:
#=
    Copyright (c) Yuki Ono.
    Licensed under the MIT License.
=#

using SendEML
using Test

function make_simple_mail_text()::String
    text = """From: a001 <a001@ah62.example.jp>
Subject: test
To: a002@ah62.example.jp
Message-ID: <b0e564a5-4f70-761a-e103-70119d1bcb32@ah62.example.jp>
Date: Sun, 26 Jul 2020 22:01:37 +0900
User-Agent: Mozilla/5.0 (Windows NT 10.0; Win64; x64; rv:78.0) Gecko/20100101
 Thunderbird/78.0.1
MIME-Version: 1.0
Content-Type: text/plain; charset=utf-8; format=flowed
Content-Transfer-Encoding: 7bit
Content-Language: en-US

test"""
    replace(text, "\n" => "\r\n")
end

function make_folded_mail()::Vector{UInt8}
    text = """From: a001 <a001@ah62.example.jp>
Subject: test
To: a002@ah62.example.jp
Message-ID:
 <b0e564a5-4f70-761a-e103-70119d1bcb32@ah62.example.jp>
Date:
 Sun, 26 Jul 2020
 22:01:37 +0900
User-Agent: Mozilla/5.0 (Windows NT 10.0; Win64; x64; rv:78.0) Gecko/20100101
 Thunderbird/78.0.1
MIME-Version: 1.0
Content-Type: text/plain; charset=utf-8; format=flowed
Content-Transfer-Encoding: 7bit
Content-Language: en-US

test"""
    Vector{UInt8}(replace(text, "\n" => "\r\n"))
end

function make_folded_end_date()::Vector{UInt8}
    text = """From: a001 <a001@ah62.example.jp>
Subject: test
To: a002@ah62.example.jp
Message-ID:
 <b0e564a5-4f70-761a-e103-70119d1bcb32@ah62.example.jp>
User-Agent: Mozilla/5.0 (Windows NT 10.0; Win64; x64; rv:78.0) Gecko/20100101
 Thunderbird/78.0.1
MIME-Version: 1.0
Content-Type: text/plain; charset=utf-8; format=flowed
Content-Transfer-Encoding: 7bit
Content-Language: en-US
Date:
 Sun, 26 Jul 2020
 22:01:37 +0900
"""
    Vector{UInt8}(replace(text, "\n" => "\r\n"))
end

function make_folded_end_message_id()::Vector{UInt8}
    text = """From: a001 <a001@ah62.example.jp>
Subject: test
To: a002@ah62.example.jp
User-Agent: Mozilla/5.0 (Windows NT 10.0; Win64; x64; rv:78.0) Gecko/20100101
 Thunderbird/78.0.1
MIME-Version: 1.0
Content-Type: text/plain; charset=utf-8; format=flowed
Content-Transfer-Encoding: 7bit
Content-Language: en-US
Date:
 Sun, 26 Jul 2020
 22:01:37 +0900
Message-ID:
 <b0e564a5-4f70-761a-e103-70119d1bcb32@ah62.example.jp>
"""
    Vector{UInt8}(replace(text, "\n" => "\r\n"))
end

function make_simple_mail()::Vector{UInt8}
    Vector{UInt8}(make_simple_mail_text())
end

function make_invalid_mail()::Vector{UInt8}
    Vector{UInt8}(replace(make_simple_mail_text(), "\r\n\r\n" => ""))
end

function make_test_send_cmd(expected::String)::Function
    cmd -> begin
        @test expected == cmd
        cmd
    end
end

# ! Avoid String(): -> String(copy())
# String constructor truncates data. #32528
# https://github.com/JuliaLang/julia/issues/32528

function get_header_line(header::Vector{UInt8}, name::String)::String
    match(name * r":[\s\S]+?\r\n(?=([^ \t]|$))", String(copy(header)), 1).match
end

function get_date_line(header::Vector{UInt8})::String
    get_header_line(header, "Date")
end

function get_message_id_line(header::Vector{UInt8})::String
    get_header_line(header, "Message-ID")
end

@testset "SendEML" begin
    @testset "make_simple_mail_text" begin
        text = make_simple_mail_text()
        lines = split(text, "\r\n")
        @test length(lines) == 13
        @test lines[1] == "From: a001 <a001@ah62.example.jp>"
        @test lines[13] == "test"
    end

    @testset "make_simple_mail" begin
        mail = make_simple_mail()
        @test typeof(mail) == Vector{UInt8}
    end

    @testset "get_header_line" begin
        mail = make_simple_mail()
        @test "Date: Sun, 26 Jul 2020 22:01:37 +0900\r\n" == get_date_line(mail)
        @test "Message-ID: <b0e564a5-4f70-761a-e103-70119d1bcb32@ah62.example.jp>\r\n" == get_message_id_line(mail)

        f_mail = make_folded_mail()
        @test "Date:\r\n Sun, 26 Jul 2020\r\n 22:01:37 +0900\r\n" == get_date_line(f_mail)
        @test "Message-ID:\r\n <b0e564a5-4f70-761a-e103-70119d1bcb32@ah62.example.jp>\r\n" == get_message_id_line(f_mail)

        e_date = make_folded_end_date()
        @test "Date:\r\n Sun, 26 Jul 2020\r\n 22:01:37 +0900\r\n" == get_date_line(e_date)

        e_message_id = make_folded_end_message_id()
        @test "Message-ID:\r\n <b0e564a5-4f70-761a-e103-70119d1bcb32@ah62.example.jp>\r\n" == get_message_id_line(e_message_id)
    end

     @testset "match_header" begin
        function match(s1::String, s2::String)::Bool
            SendEML.match_header(Vector{UInt8}(s1), Vector{UInt8}(s2))
        end

        @test match("Test:", "Test:") == true
        @test match("Test:   ", "Test:") == true
        @test match("Test: xxx", "Test:") == true

        @test match("", "Test:") == false
        @test match("T", "Test:") == false
        @test match("Test", "Test:") == false
        @test match("Xest:", "Test:") == false

        @test_throws ErrorException match("Test: xxx", "")
    end

    @testset "find_cr" begin
        mail = make_simple_mail()
        @test SendEML.find_cr(mail, 1) == 34
        @test SendEML.find_cr(mail, 35) == 49
        @test SendEML.find_cr(mail, 59) == 75
    end

    @testset "find_lf" begin
        mail = make_simple_mail()
        @test SendEML.find_lf(mail, 1) == 35
        @test SendEML.find_lf(mail, 36) == 50
        @test SendEML.find_lf(mail, 60) == 76
    end

    @testset "find_all_lf" begin
        mail = make_simple_mail()
        indices = SendEML.find_all_lf(mail)

        @test indices[1] == 35
        @test indices[2] == 50
        @test indices[3] == 76

        @test indices[end-2] == 391
        @test indices[end-1] == 416
        @test indices[end] == 418
    end

    @testset "get_lines" begin
        mail = make_simple_mail()
        lines = SendEML.get_lines(mail)

        @test length(lines) == 13

        @test lines[1] == b"From: a001 <a001@ah62.example.jp>\r\n"
        @test lines[2] == b"Subject: test\r\n"
        @test lines[3] == b"To: a002@ah62.example.jp\r\n"

        @test lines[end-2] == b"Content-Language: en-US\r\n"
        @test lines[end-1] == b"\r\n"
        @test lines[end] == b"test"
    end

    @testset "is_not_update" begin
        @test SendEML.is_not_update(true, true) == false
        @test SendEML.is_not_update(true, false) == false
        @test SendEML.is_not_update(false, true) == false
        @test SendEML.is_not_update(false, false) == true
    end

    @testset "make_now_date_line" begin
        line = SendEML.make_now_date_line()
        @test startswith(line, "Date: ")
        @test endswith(line, SendEML.CRLF)
        @test length(line) <= 80
    end

    @testset "make_random_message_id_line" begin
        line = SendEML.make_random_message_id_line()
        @test startswith(line, "Message-ID: ")
        @test endswith(line, SendEML.CRLF)
        @test length(line) <= 80
    end

    @testset "is_wsp" begin
        @test SendEML.is_wsp(UInt8(' ')) == true
        @test SendEML.is_wsp(UInt8('\t')) == true
        @test SendEML.is_wsp(UInt8('\0')) == false
        @test SendEML.is_wsp(UInt8('a')) == false
        @test SendEML.is_wsp(UInt8('b')) == false
    end

    @testset "first_byte" begin
        @test SendEML.first_byte(SendEML.DATE_BYTES, UInt8('0')) == UInt8('D')
        @test SendEML.first_byte(SendEML.MESSAGE_ID_BYTES, UInt8('0')) == UInt8('M')
        @test SendEML.first_byte(Vector{UInt8}(), UInt8('0')) == UInt8('0')
    end

    @testset "first_char" begin
        @test SendEML.first_char("Date:", '0') == 'D'
        @test SendEML.first_char("Message-ID:", '0') == 'M'
        @test SendEML.first_char("", '0') == '0'
    end

    @testset "is_folded_line" begin
        function match(chars::Vararg{Char})::Bool
            SendEML.is_folded_line(map(UInt8, collect(Iterators.flatten(chars))))
        end

        @test match(' ', 'a', 'b') == true
        @test match('\t', 'a', 'b') == true
        @test match('\0', 'a', 'b') == false
        @test match('a', 'a', ' ') == false
        @test match('b', 'a', '\t') == false
    end

    @testset "replace_date_line" begin
        f_mail = make_folded_mail()
        lines = SendEML.get_lines(f_mail)
        new_lines = SendEML.replace_date_line(lines)
        @test lines != new_lines

        new_mail = SendEML.concat_bytes(new_lines)
        @test f_mail != new_mail
        @test get_date_line(f_mail) != get_date_line(new_mail)
        @test get_message_id_line(f_mail) == get_message_id_line(new_mail)
    end

    @testset "replace_message_id_line" begin
        f_mail = make_folded_mail()
        lines = SendEML.get_lines(f_mail)
        new_lines = SendEML.replace_message_id_line(lines)
        @test lines != new_lines

        new_mail = SendEML.concat_bytes(new_lines)
        @test f_mail != new_mail
        @test get_message_id_line(f_mail) != get_message_id_line(new_mail)
        @test get_date_line(f_mail) == get_date_line(new_mail)
    end

    @testset "replace_header" begin
        mail = make_simple_mail()
        date_line = get_date_line(mail)
        mid_line = get_message_id_line(mail)

        repl_header_noupdate = SendEML.replace_header(mail, false, false)
        @test mail == repl_header_noupdate

        repl_header = SendEML.replace_header(mail, true, true)
        @test mail != repl_header

        function replace_header(header::Vector{UInt8}, update_date::Bool, update_message_id::Bool)::Tuple{String, String}
            r_header = SendEML.replace_header(header, update_date, update_message_id)
            @test header != r_header
            (get_date_line(r_header), get_message_id_line(r_header))
        end

        (r_date_line, r_mid_line) = replace_header(mail, true, true)
        @test date_line != r_date_line
        @test mid_line != r_mid_line

        (r_date_line, r_mid_line) = replace_header(mail, true, false)
        @test date_line != r_date_line
        @test mid_line == r_mid_line

        (r_date_line, r_mid_line) = replace_header(mail, false, true)
        @test date_line == r_date_line
        @test mid_line != r_mid_line

        f_mail = make_folded_mail()
        (f_date_line, f_mid_line) = replace_header(f_mail, true, true)
        @test count(c -> (c == '\n'), f_date_line) == 1
        @test count(c -> (c == '\n'), f_mid_line) == 1
    end

    @testset "concat_bytes" begin
        mail = make_simple_mail()
        lines = SendEML.get_lines(mail)
        new_mail = SendEML.concat_bytes(lines)
        @test mail == new_mail
    end

    @testset "combine_mail" begin
        mail = make_simple_mail()
        (header, body) = SendEML.split_mail(mail)
        new_mail = SendEML.combine_mail(header, body)
        @test mail == new_mail
    end

    import SendEML.CR
    import SendEML.LF

    @testset "has_next_lf_cr_lf" begin
        match = (bytes::Vector{UInt8}, idx::Int) -> SendEML.has_next_lf_cr_lf(bytes, idx)
        zero = UInt8(0)

        @test match([CR, LF, CR, LF], 1) == true
        @test match([zero, CR, LF, CR, LF], 2) == true

        @test match([CR, LF, CR, LF], 2) == false
        @test match([CR, LF, CR, zero], 1) == false
        @test match([CR, LF, CR, LF, zero], 2) == false
    end

    @testset "find_empty_line" begin
        mail = make_simple_mail()
        @test SendEML.find_empty_line(mail) == 415

        invalid_mail = make_invalid_mail()
        @test SendEML.find_empty_line(invalid_mail) === nothing
    end

    @testset "split_mail" begin
        mail = make_simple_mail()
        header_body = SendEML.split_mail(mail)
        @test header_body !== nothing

        (header, body) = header_body
        @test mail[1:414] == header
        @test mail[(415 + 4):end] == body

        invalid_mail = make_invalid_mail()
        @test SendEML.split_mail(invalid_mail) === nothing
    end

    @testset "replace_mail" begin
        mail = make_simple_mail()
        repl_mail_noupdate = SendEML.replace_mail(mail, false, false)
        @test mail == repl_mail_noupdate

        repl_mail = SendEML.replace_mail(mail, true, true)
        @test mail != repl_mail
        @test mail[end-100:end] == repl_mail[end-100:end]

        invalid_mail = make_invalid_mail()
        @test isnothing(SendEML.replace_mail(invalid_mail, true, true))
    end

    @testset "get_and_map_settings" begin
        text = SendEML.make_json_sample()
        settings = SendEML.map_settings(SendEML.get_settings_from_text(text))

        @test settings.smtp_host == "172.16.3.151"
        @test settings.smtp_port == 25
        @test settings.from_address == "a001@ah62.example.jp"
        @test settings.to_addresses == ["a001@ah62.example.jp", "a002@ah62.example.jp", "a003@ah62.example.jp"]
        @test settings.eml_files == ["test1.eml", "test2.eml", "test3.eml"]
        @test settings.update_date == true
        @test settings.update_message_id == true
        @test settings.use_parallel == false
    end

    @testset "is_last_reply" begin
        @test SendEML.is_last_reply("250-First line") == false
        @test SendEML.is_last_reply("250-Second line") == false
        @test SendEML.is_last_reply("250-234 Text beginning with numbers") == false
        @test SendEML.is_last_reply("250 The last line") == true
    end

    @testset "is_positive_reply" begin
        @test SendEML.is_positive_reply("200 xxx") == true
        @test SendEML.is_positive_reply("300 xxx") == true
        @test SendEML.is_positive_reply("400 xxx") == false
        @test SendEML.is_positive_reply("500 xxx") == false
        @test SendEML.is_positive_reply("xxx 200") == false
        @test SendEML.is_positive_reply("xxx 300") == false
    end

    @testset "replace_crlf_dot" begin
        @test SendEML.replace_crlf_dot("TEST") == "TEST"
        @test SendEML.replace_crlf_dot("CRLF") == "CRLF"
        @test SendEML.replace_crlf_dot(SendEML.CRLF) == SendEML.CRLF
        @test SendEML.replace_crlf_dot(".") == "."
        @test SendEML.replace_crlf_dot("$(SendEML.CRLF).") == "<CRLF>."
    end

    @testset "send_hello" begin
        SendEML.send_hello(make_test_send_cmd("EHLO localhost"))
    end

    @testset "send_from" begin
        SendEML.send_from(make_test_send_cmd("MAIL FROM: <a001@ah62.example.jp>"), "a001@ah62.example.jp")
    end

    @testset "send_rcpt_to" begin
        count = 1
        test_func = cmd -> begin
            @test "RCPT TO: <a00$(count)@ah62.example.jp>" == cmd
            count += 1
            cmd
        end

        SendEML.send_rcpt_to(test_func, ["a001@ah62.example.jp", "a002@ah62.example.jp", "a003@ah62.example.jp"])
    end

    @testset "send_data" begin
        SendEML.send_data(make_test_send_cmd("DATA"))
    end

    @testset "send_crlf_dot" begin
        SendEML.send_crlf_dot(make_test_send_cmd("$(SendEML.CRLF)."))
    end

    @testset "send_quit" begin
        SendEML.send_quit(make_test_send_cmd("QUIT"))
    end

    @testset "send_rset" begin
        SendEML.send_rset(make_test_send_cmd("RSET"))
    end

    @testset "check_settings" begin
        function check_no_key(key::String)
            json = SendEML.make_json_sample()
            no_key = replace(json, key => "X-$key")
            SendEML.check_settings(SendEML.get_settings_from_text(no_key))
        end

        @test_throws ErrorException check_no_key("smtpHost")
        @test_throws ErrorException check_no_key("smtpPort")
        @test_throws ErrorException check_no_key("fromAddress")
        @test_throws ErrorException check_no_key("toAddresses")
        @test_throws ErrorException check_no_key("emlFiles")

        @test_nowarn check_no_key("updateDate")
        @test_nowarn check_no_key("updateMessageId")
        @test_nowarn check_no_key("useParallel")
    end

    @testset "proc_json_file" begin
        @test_throws ErrorException SendEML.proc_json_file("__test__")
    end

    import JSON

    @testset "check_json_value" begin
        function check(json::String, type::Type)
            SendEML.check_json_value(JSON.parse(json), "test", type)
        end

        function check_error(json::String, type::Type, expected::String)
            try
                check(json, type)
                @test false
            catch e
                @test isa(e, ErrorException)
                @test e.msg == expected
            end
        end

        json = """{"test": "172.16.3.151"}"""
        @test_nowarn check(json, String)
        @test_throws ErrorException check(json, Int)
        check_error(json, Bool, "test: Invalid type: 172.16.3.151")

        json = """{"test": 172}"""
        @test_nowarn check(json, Int)
        @test_throws ErrorException check(json, String)
        check_error(json, Bool, "test: Invalid type: 172")

        json = """{"test": true}"""
        @test_nowarn check(json, Bool)
        @test_nowarn check(json, Int) # true <=> 1
        check_error(json, String, "test: Invalid type: true")

        json = """{"test": false}"""
        @test_nowarn check(json, Bool)
        @test_nowarn check(json, Int) # false <=> 0
        check_error(json, String, "test: Invalid type: false")

        json = """{"test": ["172.16.3.151", "172.16.3.152", "172.16.3.153"]}"""
        @test_nowarn check(json, Vector{String})
        check_error(json, String, "test: Invalid type: Any[\"172.16.3.151\", \"172.16.3.152\", \"172.16.3.153\"]")

        json = """{"test": ["172.16.3.151", "172.16.3.152", 172]}"""
        check_error(json, Vector{String}, "test: Invalid type: Any[\"172.16.3.151\", \"172.16.3.152\", 172]")
    end
end

ext: jl

max_stars_repo_name: ykon/sendeml-jl

max_issues_repo_name: ykon/sendeml-jl

max_forks_repo_name: ykon/sendeml-jl


text:
using Documenter
using DarkSky

makedocs(
    modules = [DarkSky],
    doctest = false
)

deploydocs(
    deps   = Deps.pip("pygments", "mkdocs", "mkdocs-material", "python-markdown-math"),
    repo   = "github.com/ellisvalentiner/DarkSky.jl.git",
    julia  = "0.7"
)


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

name = "MPSolve"
version = v"3.2.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/robol/MPSolve.git", "65be6767119a3df911943d7aaacdad4646fa6173")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/MPSolve/
export CPPFLAGS="-I${prefix}/include"
./autogen.sh
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
cd src/libmps/
make
make install
exit
"""

platforms = supported_platforms()
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libmps", :libmps)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="GMP_jll", uuid="781609d7-10c4-51f6-84f2-b8444358ff6d"), v"6.1.2"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module Chebyshev

using FFTW
using Base: copy!

import ..Utils: dz_partials, value, Dual, ComplexDual, Partials, valtype

import Base: *, \
struct Transform{T,I,D}
    dct!::FFTW.r2rFFTWPlan{T,(3,),true,D}
end

function plan_transform!(x::Vector{T}) where T
    Transform{T,true,1}(FFTW.plan_r2r!(x, FFTW.REDFT00))
end

function plan_transform(x::Vector{T}) where T
    Transform{T,false,1}(FFTW.plan_r2r!(x, FFTW.REDFT00))
end

function plan_transform!(x::Array{T,2}) where T
    Transform{T,true,2}(FFTW.plan_r2r!(x, FFTW.REDFT00,1))
end

function plan_transform(x::Array{T,2}) where T
    Transform{T,false,2}(FFTW.plan_r2r!(x, FFTW.REDFT00,1))
end

(C::Transform{T,true}  * x::Vector{T}) where T = transform!(x, C.dct!)
(C::Transform{T,false} * x::Vector{T}) where T = transform(x, C.dct!)
(C::Transform{T,true}  \ A::Vector{T}) where T = inv_transform!(A, C.dct!)
(C::Transform{T,false} \ A::Vector{T}) where T = inv_transform(A, C.dct!)

(C::Transform{Complex{T},true} * x::Vector{Complex{S}}) where {T,S <:Dual{R,V,N}} where {R,V,N}  = transform!(x, Array{Complex{V}}(undef,length(x),N+1), C.dct!)
(C::Transform{Complex{T},false} * x::Vector{Complex{S}}) where {T,S <:Dual{R,V,N}} where {R,V,N}  = transform(x, Array{Complex{V}}(undef,length(x),N+1),C.dct!)
(C::Transform{Complex{T},true} \ A::Vector{Complex{S}}) where {T,S <:Dual{R,V,N}} where {R,V,N}= inv_transform!(A, Array{Complex{V}}(undef,length(x),N+1),C.dct!)
(C::Transform{Complex{T},false} \ A::Vector{Complex{S}}) where {T,S <:Dual{R,V,N}} where {R,V,N} = inv_transform(A, Array{Complex{V}}(undef,length(x),N+1),C.dct!)

#=
(C::Transform{T,true} * x::Vector{S}) where {T<:Real,S <:Dual{R} where R} = transform!(x, C.dct!)
(C::Transform{T,false} * x::Vector{S}) where {T<:Real,S <:Dual{R} where R} = transform(x, C.dct!)
(C::Transform{T,true} \ A::Vector{S}) where {T<:Real,S <:Dual{R} where R} = inv_transform!(A, C.dct!)
(C::Transform{T,false} \ A::Vector{S}) where {T<:Real,S <:Dual{R} where R} = inv_transform(A, C.dct!)
=#

"""
    Chebyshev.nodes(N, T = Float64)

Gives `N` extrema Chebyshev points of type `T` in [-1 ,1].
"""
function nodes(N, T = Float64)
    T[cos(θ) for θ in range(π, 0, length=N)]
end

"""
    Chebyshev.firstkind(C, s[, offset = 0])

Evaluate a Chebyshev series of the first kind.

# Arguments

- `C`: a vector with the coefficients of the series
- `s`: a scalar or an array of evaluation points ∈ [-1, 1]
- `offset`: an optional starting index of the series summation
"""
firstkind(A, ss, offset = 0) = [firstkind(A, s, offset) for s in ss]
function firstkind(A::AbstractArray{C}, s::S, offset = 0)::C where {C,S <: Real}
    -1 ≤ s ≤ 1 || throw(DomainError("s ∉ [-1,1]"))
    T₋ = s
    T  = one(S)

    for i in 1:offset
        T₋, T = T, 2s*T - T₋
    end

    N = length(A)
    f = zero(C)
    for i in 1+offset:N
        f += A[i]*T
        T₋, T = T, 2s*T - T₋
    end
    f
end

"""
    Chebyshev.secondkind(C, s[, offset = 0])

Evaluate a Chebyshev series of the second kind.

# Arguments

- `C`: a vector with the coefficients of the series
- `s`: a scalar or an array of evaluation points ∈ [-1, 1]
- `offset`: an optional starting index of the series summation
"""
secondkind(C, ss, offset = 0) = [secondkind(C, s, offset) for s in ss]
function secondkind(A::AbstractArray{C}, s::S, offset = 0)::C where {C,S <: Real}
    -1 ≤ s ≤ 1 || throw(DomainError("s ∉ [-1,1]"))
    U₋ = zero(s)
    U  = one(S)

    for i in 1:offset
        U₋, U = U, 2s*U - U₋
    end

    N = length(A)
    f = zero(C)
    for i in 1+offset:N
        f += A[i]*U
        U₋, U = U, 2s*U - U₋
    end
    f
end

"""
    Chebyshev.transform(x[, plan!])

Compute the discrete Chebyshev transform of `x`.

# Arguments

- `x`: samples of a function distributed along [extrema Chebyshev nodes](@ref Chebyshev.nodes)
- `plan!`: an optional pre-planned in-place DCT-I used to compute the transform
"""
function transform(x, plan! = FFTW.plan_r2r!(x, FFTW.REDFT00))
    C = copy(x)
    transform!(C, plan!)
end

"""
    Chebyshev.transform!(A,[ x, plan!])

Compute the discrete Chebyshev transform of `x` in-place.

# Arguments

- `A`: the output vector (also the input vector if `x` is omitted)
- `x`: an optional input vector with samples of a function distributed
  along [extrema Chebyshev nodes](@ref Chebyshev.nodes)
- `plan!`: an optional pre-planned in-place DCT-I used to compute the transform
"""
function transform!(x, plan! = FFTW.plan_r2r!(x, FFTW.REDFT00))
    N = length(x)
    if N != length(plan!)
        throw(BoundsError("`x` must have the same size as the preplanned DCT"))
    end

    plan!*x

    s = 1/(N-1)

    @inbounds begin
        for n in 2:2:N
            x[n-1] *= s
            x[n]   *= -s
        end
        if isodd(N)
            x[N] *= s
        end

        x[1] /= 2
        x[N] /= 2
    end

    x
end

function transform!(A::T, x::T, plan! = FFTW.plan_r2r!(x, FFTW.REDFT00)) where T
    copy!(A, x)
    transform!(A, plan!)
end

function transform!(x::AbstractArray{V,2}, plan! = FFTW.plan_r2r!(x, FFTW.REDFT00,1)) where V
    N = size(x,1)
    if size(x) != size(plan!)
        throw(BoundsError("`x` must have the same size as the preplanned DCT"))
    end

    plan!*x

    s = 1/(N-1)

    @inbounds begin
        for n in 2:2:N
            x[n-1,:] *= s
            x[n,:]   *= -s
        end
        if isodd(N)
            x[N,:] *= s
        end

        x[1,:] /= 2
        x[N,:] /= 2
    end

    x
end

transform!(x::Vector{Complex{S}}, xhat = Array{Complex{V}}(undef,length(x),N+1),
                    plan! = FFTW.plan_r2r!(xhat, FFTW.REDFT00,1)) where S <:Dual{T,V,N} where {T,V,N} =
   _transform!(x,transform!,xhat,plan!)

"""
    Chebyshev.inv_transform(A[, plan!])

Perform the inverse Chebyshev transform.

This is the inverse of [`Chebyshev.transform`](@ref).

# Fields

- `A`: the coefficients of a Chebyshev series
- `plan!`: an optional pre-planned in-place DCT-I used to compute the transform
"""
function inv_transform(A, plan! = FFTW.plan_r2r!(A, FFTW.REDFT00))
    x = copy(A)
    inv_transform!(x, plan!)
end

"""
    Chebyshev.inv_transform!(x, [C, plan!)

Perform the inverse Chebyshev transform in place.

This is the inverse of [`Chebyshev.transform!`](@ref)

# Arguments

- `x`: the output vector to store the reconstructed function with
  points distributed along extrema Chebyshev nodes (also the input
  vector if `A` is omitted)
- `A`: an optional input vector with coefficients of a Chebyshev series
- `plan!`: an optional pre-planned in-place DCT-I used to compute the transform
"""
function inv_transform!(A, plan! = FFTW.plan_r2r!(A, FFTW.REDFT00,1))
    N = length(A)
    if N != length(plan!)
        throw(BoundsError("`A` must have the same size as the preplanned DCT"))
    end

    @inbounds begin
        for n in 2:2:N-2
            A[n]   *= -0.5
            A[n+1] *=  0.5
        end
        if iseven(N)
            A[N] *= -1
        else
            A[N-1] *= -0.5
        end
    end

    plan!*A
end

function inv_transform!(x::T, A::T, plan! = FFTW.plan_r2r!(A, FFTW.REDFT00)) where T
    copy!(x, A)
    Chebyshev.inv_transform!(x, plan!)
end

function inv_transform!(A::AbstractArray{V,2}, plan! = FFTW.plan_r2r!(A, FFTW.REDFT00,1)) where {V}
    N = size(A,1)
    if size(A) != size(plan!)
      throw(BoundsError("`x` must have the same size as the preplanned DCT"))
    end

    @inbounds begin
        for n in 2:2:N-2
            A[n,:]   *= -0.5
            A[n+1,:] *=  0.5
        end
        if iseven(N)
            A[N,:] *= -1
        else
            A[N-1,:] *= -0.5
        end
    end

    plan!*A
end

inv_transform!(x::Vector{Complex{S}}, xhat = Array{Complex{V}}(undef,length(x),N+1),
                    plan! = FFTW.plan_r2r!(xhat, FFTW.REDFT00,1)) where S <:Dual{T,V,N} where {T,V,N} =
   _transform!(x,inv_transform!,xhat,plan!)

function _transform!(x::Vector{Complex{S}}, f! = transform!, xhat = Array{Complex{V}}(undef,length(x),N+1),
                    plan! = FFTW.plan_r2r!(xhat, FFTW.REDFT00,1)) where S <:Dual{T,V,N} where {T,V,N}

    xz = first.(dz_partials.(x))
    xzstar = last.(dz_partials.(x))

    xhat[:,1] = value.(x)
    for j in eachindex(x)
        xhat[j,2:N÷2+1] = xz[j]
        xhat[j,N÷2+2:N+1] = xzstar[j]
    end
    f!(xhat,xhat,plan!)
    for j in eachindex(x)
        x[j] = ComplexDual{T}(xhat[j,1],Partials(tuple(xhat[j,2:N÷2+1]...)),
                                        Partials(tuple(xhat[j,N÷2+2:N+1]...)))
    end

    return x
end

"""
    clenshaw_curtis_weights(N)

Returns the `N` weights used in the Clenshaw Curtis quadrature
"""
function clenshaw_curtis_weights(N)
    n = N - 1
    w = zeros(N)

    # Unlike the documentation, we move the 1/(N-1)
    # factor inside the r vector.
    w[1] = 1.0/n
    for k in 2:2:n
        w[k+1] = (1.0 + (-1)^k)/(1 - k^2)/n
    end

    # Scale by S
    w[1] *= 2; w[N] *= 2
    # Multiply by D
    FFTW.r2r!(w, FFTW.REDFT00)
    # Scale by S⁻¹
    w[1] /= 2; w[N] /= 2

    return w
end

end


text:
module Bisect_Tests

using Compat
using Compat.Test
using IntervalArithmetic
using EAGOBranchBound

S = BnBSolver()
B = BnBModel([Interval(1.0,2.0),Interval(3.0,4.0),Interval(3.0,5.0),Interval(3.0,6.0)])
N = [Interval(1.0,2.0),Interval(3.0,4.0),Interval(4.0,5.0),Interval(4.0,6.0)]
nx = 2

X1,X2 = EAGOBranchBound.Bisect_Rel(S,B,N)
Y1,Y2 = EAGOBranchBound.Bisect_Abs_Imp(S,B,N,nx)
Z1,Z2 = EAGOBranchBound.Bisect_Rel_Imp(S,B,N,nx)
@test X1[1] == Interval(1.0,1.5)
@test X2[1] == Interval(1.5,2.0)
@test X2[2:4] == [Interval(3.0,4.0),Interval(4.0,5.0),Interval(4.0,6.0)]
@test Y1[1:3] == [Interval(1.0,2.0),Interval(3.0,4.0),Interval(4.0,5.0)]
@test Y1[4] == Interval(4.0,5.0)
@test Y2[4] == Interval(5.0,6.0)
@test Z1[1:3] == [Interval(1.0,2.0),Interval(3.0,4.0),Interval(4.0,5.0)]
@test Z1[4] == Interval(4.0,5.0)
@test Z2[4] == Interval(5.0,6.0)

end


text:
using DashBootstrapComponents, DashHtmlComponents

PLOTLY_LOGO = "https://images.plot.ly/logo/new-branding/plotly-logomark.png";

search_bar = dbc_row(
    [
        dbc_col(dbc_input(type = "search", placeholder = "Search")),
        dbc_col(
            dbc_button("Search", color = "primary", className = "ms-2", n_clicks = 0),
            width = "auto",
        ),
    ],
    className = "g-0 ms-auto flex-nowrap mt-3 mt-md-0",
    align = "center",
);

navbar = dbc_navbar(
    dbc_container([
        html_a(
            # Use row and col to control vertical alignment of logo / brand
            dbc_row(
                [
                    dbc_col(html_img(src = PLOTLY_LOGO, height = "30px")),
                    dbc_col(dbc_navbarbrand("Navbar", className = "ms-2")),
                ],
                align = "center",
                className = "g-0",
            ),
            href = "https://plotly.com",
            style = Dict("textDecoration" => "none"),
        ),
        dbc_navbartoggler(id = "navbar-toggler", n_clicks = 0),
        dbc_collapse(search_bar, id = "navbar-collapse", is_open = false, navbar = true),
    ]),
    color = "dark",
    dark = true,
)

# add callback for toggling the collapse on small screens
callback!(
    app,
    Output("navbar-collapse", "is_open"),
    Input("navbar-toggler", "n_clicks"),
    State("navbar-collapse", "is_open"),
) do n, is_open
    return n > 0 ? is_open == 0 : is_open
end;

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
using Documenter, ChiralWaveguide

makedocs(
  sitename = "ChiralWaveguide",
  modules = [ChiralWaveguide],
  format = Documenter.HTML(assets = ["assets/favicon.ico"]),
  pages = [
    "index.md",
    # "theory.md",
    "Examples" => [
      "examples/CoherentInput.md",
      "examples/SinglePhotonScattering.md",
      "examples/CustomModes.md",
      "examples/CustomSystems.md"
    ],
    "api.md"
  ]
)

deploydocs(
  repo = "github.com/KKleinbeck/ChiralWaveguide.jl.git",
)


text:
using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkitStandardLibrary.Blocks: _clamp, _dead_zone

@parameters t

@testset "Limiter" begin
    @testset "Constant" begin
        @named c = Constant(; k=1)
        @named int = Integrator(; k=1)
        @named sat = Limiter(; y_min=-0.6, y_max=0.8)
        @named model = ODESystem([
                connect(c.output, int.input),
                connect(int.output, sat.input),
            ],
            t,
            systems=[int, c, sat],
        )
        sys = structural_simplify(model)
        prob = ODEProblem(sys, [int.x=>1.0], (0.0, 1.0))

        sol = solve(prob, Rodas4())
        @test sol[int.output.u][end] ≈ 2
        @test sol[sat.output.u][end] ≈ 0.8
    end

    @testset "Sine" begin
        y_min, y_max = -0.3, 0.5
        @named source = Sine(; frequency=1/2)
        @named lim = Limiter(; y_max=y_max, y_min=y_min)
        @named int = Integrator(; k=1)
        @named iosys = ODESystem([
                connect(source.output, lim.input),
                connect(lim.output, int.input),
            ],
            t,
            systems=[source, lim, int],
        )
        sys = structural_simplify(iosys)

        prob = ODEProblem(sys, Pair[], (0.0, 10.0))

        sol = solve(prob, Rodas4())
        @test all(abs.(sol[lim.output.u]) .<= 0.5)
        @test all(isapprox.(sol[lim.output.u], _clamp.(sol[source.output.u], y_min, y_max), atol=1e-2))

        # Plots.plot(sol; vars=[source.output.u, lim.output.u])
        # Plots.scatter(sol[source.output.u], sol[lim.output.u])
        # Plots.scatter!(sol[source.output.u], _clamp.(sol[source.output.u], y_min, y_max))
    end
end

@testset "DeadZone" begin
    @testset "Constant" begin
        @named c = Constant(; k=1)
        @named int = Integrator(; k=1)
        @named dz = DeadZone(; u_min=-2, u_max=1)
        @named model = ODESystem([
                connect(c.output, int.input),
                connect(int.output, dz.input),
            ],
            t,
            systems=[int, c, dz],
        )
        sys = structural_simplify(model)
        prob = ODEProblem(sys, [int.x=>1.0], (0.0, 1.0))
        sol = solve(prob, Rodas4())

        @test all(sol[int.output.u][end] .≈ 2)
    end

    @testset "Sine" begin
        u_min, u_max = -2, 1
        @named source = Sine(; amplitude=3, frequency=1/2)
        @named dz = DeadZone(; u_min=u_min, u_max=u_max)
        @named int = Integrator(; k=1)
        @named model = ODESystem([
                connect(source.output, dz.input),
                connect(dz.output, int.input),
            ],
            t,
            systems=[int, source, dz],
        )
        sys = structural_simplify(model)
        prob = ODEProblem(sys, [int.x=>1.0], (0.0, 10.0))
        sol = solve(prob, Rodas4())

        @test all(sol[dz.output.u] .<= 2)
        @test all(sol[dz.output.u] .>= -1)
        @test all(isapprox.(sol[dz.output.u], _dead_zone.(sol[source.output.u], u_min, u_max), atol=1e-2))

        # Plots.plot(sol; vars=[source.output.u, dz.output.u])
        # Plots.scatter(sol[source.output.u], sol[dz.output.u])
        # Plots.scatter!(sol[source.output.u], _dead_zone.(sol[source.output.u], u_min, u_max))
    end
end

@testset "SlewRateLimiter" begin
    @named source = Sine(; frequency=1/2)
    @named rl = SlewRateLimiter(; rising=1, falling=-1, Td=0.001, y_start=-1/3)
    @named iosys = ODESystem([
        connect(source.output, rl.input),
        ],
        t,
        systems=[source, rl],
    )
    sys = structural_simplify(iosys)

    prob = ODEProblem(sys, Pair[], (0.0, 10.0))

    tS = 0.01
    sol = solve(prob, Rodas4(), saveat=tS, abstol=1e-10, reltol=1e-10)
    @test all(abs.(sol[rl.output.u]) .<= 0.51)
    @test all(-1 - 1e-5 .<= diff(sol[rl.output.u]) ./ tS .<= 1 + 1e-5) # just an approximation
end


text:
module TextSearch
    import Base: broadcastable
    import StatsBase: fit, predict
    using SparseArrays
    using LinearAlgebra
    include("textconfig.jl")
    include("dvec.jl")
    include("io.jl")
    include("basicmodels.jl")
    include("distmodel.jl")
    include("entmodel.jl")
    include("svecutils.jl")
    include("invindex.jl")
    include("rocchio.jl")
    include("neardup.jl")
end


text:
using Documenter, ErrorfreeArithmetic

makedocs(
    modules = [ErrorfreeArithmetic],
    format = :html,
    sitename = "ErrorfreeArithmetic",
    pages = [
        "Package" => "index.md",
        "Errorfree Arithmetic" => "intro.md",
        "Basic usage" => "usage.md"
    ]
)

deploydocs(
    repo   = "github.com/JeffreySarnoff/ErrorfreeArithmetic.jl.git",
    target = "build"
)


text:
# Autogenerated wrapper script for fts_jll for powerpc64le-linux-gnu
export libfts

JLLWrappers.@generate_wrapper_header("fts")
JLLWrappers.@declare_library_product(libfts, "libfts.so.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libfts,
        "lib/libfts.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Translation table for the types
#
# The table does not include types `array` and `object`, which are handled differently
_type_translator = Dict{String,String}(
    "integer" => "Int",
    "string" => "String",
    "number" => string(typeof(1.0)),
    "boolean" => "Bool",
    "null" => "Nothing",
)
"""
    set_translation(pairs::Pair{String, String}...)

The translation table from JSONSchema type to Julia type is defined initially by

    Dict{String,String}(
        "integer" => "Int",
        "string" => "String",
        "number" => string(typeof(1.0)),
        "boolean" => "Bool",
        "null" => "Nothing",
    )

The function `set_translation()` can be used to change one or more translations.

# Example

```jldoctest
julia> set_translation("integer" => "Int16")
Dict{String, String} with 5 entries:
  "string"  => "String"
  "number"  => "Float64"
  "integer" => "Int16"
  "null"    => "Nothing"
  "boolean" => "Bool"

julia> set_translation("integer" => "Int32", "number" => "Float32")
Dict{String, String} with 5 entries:
  "string"  => "String"
  "number"  => "Float32"
  "integer" => "Int32"
  "null"    => "Nothing"
  "boolean" => "Bool"

julia> set_translation("integer" => "Integer", "number" => "Number", "string" => "AbstractString")
  Dict{String, String} with 5 entries:
    "string"  => "AbstractString"
    "number"  => "Number"
    "integer" => "Integer"
    "null"    => "Nothing"
    "boolean" => "Bool"
```
"""
function set_translation(pairs::Pair{String, String}...)
    for (k,v) in pairs
        if haskey(_type_translator, k)
            _type_translator[k] = v
        else
            push!(_type_translator, k => v)
        end
    end
    return _type_translator
end

"""
    pascal_case(s::AbstractString)

Return a Pascal case version of the string `s`.

Pascal case changes, if necessary, every initial word character to uppercase, removes
every word separator ` `, `-`, `_`, and `.`, and does not change the case of the remaining
characters.

# Example

```jldoctest
julia> pascal_case("foo_bar")
"FooBar"

julia> pascal_case("foo_bar.baz")
"FooBarBaz"

julia> pascal_case("fOO Bar-bAz")
"FOOBarBAz"
```
"""
pascal_case(s::AbstractString) =
    replace(titlecase(s, strict=false), r" |-|_|\." => "")

function generate_constructor(json_schema::JSON3.Object, ::Val{:oneOf})
    nothing
end

function generate_constructor(json_schema::JSON3.Object, ::Val{:allOf})
    nothing
end

function generate_constructor(json_schema::JSON3.Object, ::Val{:anyOf})
    nothing
end

function generate_constructor(json_schema::JSON3.Object, ::Val{:enum})
    nothing
end

function generate_types(json_schema::JSON3.Object;
        gen_folder::String=".", module_name::String="", filename::String="")

    if Symbol("\$id") ∉ keys(json_schema)
        throw(ArgumentError("json schema missing `:\$id` keyword. No module created."))
    end

    if module_name == ""
        generated_string = generate_types("", json_schema)
    else
        generated_string = "module $module_name\n"
        generated_string *= generate_types("", json_schema)
        generated_string *= "\nend # end module"
    end

    if filename == ""
        if module_name != ""
            filename = "$module_name.jl"
        else
            filename = "JSONSchemaTypes.jl"
        end
    end

    open(joinpath(gen_folder, filename), "w") do io
        write(io, generated_string)
    end

    if module_name == ""
        @info "Types succesfully created in " *
        "$(joinpath(gen_folder, filename))"
    else
        @info "Module $module_name succesfully created in " *
            "$(joinpath(gen_folder, filename))"
    end
    return nothing
end

function generate_types(object_name::String, json_schema::JSON3.Object)

    if (:type, :properties) ⊈ keys(json_schema)
        throw(ArgumentError("Empty schema, with no `:properties` and `:type` keywords."))
    end

    if json_schema[:type] != "object"
        throw(ArgumentError("Expecting `object` type; got `$(json_schema[:type])`."))
    end

    inner_objects = Dict{String, JSON3.Object}()

    if object_name == ""
        object_name =
            :title ∈ keys(json_schema) ?
                json_schema[:title] :
                "Struct$(rand(1000:9999))"
    end

    include_docstring = false
    generated_docstring = "\n\"\"\"\n    mutable struct $object_name\n"
    if haskey(json_schema, :description)
        include_docstring = true
        generated_docstring *= "\n$(json_schema[:description])\n"
    end

    generated_struct = "\nmutable struct $object_name\n"

    include_docstring_fields = false
    generated_docstring_fields = ""
    for (k, v) in json_schema[:properties]
        if v[:type] in keys(_type_translator)
            generated_struct *= "    $(string(k))::$(_type_translator[v[:type]])\n"
        elseif v[:type] == "array"
            generated_struct *= "    $(string(k))::Array{$(_type_translator[v[:items][:type]])}\n"
        elseif v[:type] == "object"
            object_name = pascal_case(string(k))
            generated_struct *= "    $(string(k))::$(object_name)\n"
            push!(inner_objects, object_name => v)
        end
        generated_docstring_fields *= "   `$(string(k))`:"
        if haskey(v, :description)
            include_docstring_fields = true
            generated_docstring_fields *= " $(v[:description])\n"
        else
            generated_docstring_fields *= "\n"
        end
    end
    generated_struct *= "end\n" # end mutable struct

    if include_docstring_fields == true
        include_docstring == true
        generated_docstring *= "\nFields:\n" * generated_docstring_fields
    end
    if include_docstring == true
        generated_struct = generated_docstring * "\"\"\"" * generated_struct
    end

    if length(inner_objects) > 0
        for (k,v) in inner_objects
            generated_struct *= generate_types(k, v)
        end
    end

    return generated_struct
end


text:
export SymMatrix, MomentMatrix, getmat, moment_matrix, symmetric_setindex!

using SemialgebraicSets

"""
    mutable struct MomentMatrix{T, MT <: AbstractMonomial, MVT <: AbstractVector{MT}} <: AbstractMeasureLike{T}
        Q::SymMatrix{T}
        x::MVT
        support::Union{Nothing, AlgebraicSet}
    end

Measure ``\\nu`` represented by the moments of the monomial matrix ``x x^\\top`` in the symmetric matrix `Q`.
The set of points that are zeros of all the polynomials ``p`` such that ``\\mathbb{E}_{\\nu}[p] = 0`` is stored in the field `support` when it is computed.
"""
mutable struct MomentMatrix{T, MT <: AbstractMonomial, MVT <: AbstractVector{MT}} <: AbstractMeasureLike{T}
    Q::SymMatrix{T}
    x::MVT
    support::Union{Nothing, AlgebraicSet}
end
MomentMatrix{T, MT, MVT}(Q::SymMatrix{T}, x::MVT) where {T, MT, MVT} = MomentMatrix{T, MT, MVT}(Q, x, nothing)

MP.variables(μ::MomentMatrix) = variables(μ.x)
MP.nvariables(μ::MomentMatrix) = nvariables(μ.x)

function MomentMatrix{T}(f::Function, x::AbstractVector{MT}, σ) where {T, MT<:AbstractMonomial}
    MomentMatrix{T, MT, monovectype(x)}(trimat(T, f, length(x), σ), x)
end
MomentMatrix{T}(f::Function, x::AbstractVector, σ) where T = MomentMatrix{T}(f, monomial.(x), σ)
function MomentMatrix{T}(f::Function, x::AbstractVector) where T
    σ, X = sortmonovec(x)
    MomentMatrix{T}(f, X, σ)
end
MomentMatrix(f::Function, x) = MomentMatrix{Base.promote_op(f, Int, Int)}(f, x)
moment_matrix(f::Function, x) = MomentMatrix(f, x)

"""
    moment_matrix(μ::Measure, x)

Creates a matrix the moment matrix for the moment matrix  ``x x^\\top`` using the moments of `μ`.
"""
function moment_matrix(μ::Measure{T}, X) where T
    function getmom(i, j)
        x = X[i] * X[j]
        for m in moments(μ)
            if monomial(m) == x
                return moment_value(m)
            end
        end
        throw(ArgumentError("μ does not have the moment $(x)"))
    end
    MomentMatrix{T}(getmom, X)
end

function MomentMatrix(Q::AbstractMatrix{T}, x, σ) where T
    MomentMatrix{T}((i,j) -> Q[σ[i], σ[j]], x)
end
function MomentMatrix(Q::AbstractMatrix, x)
    σ, X = sortmonovec(x)
    MomentMatrix(Q, X, σ)
end
moment_matrix(Q::AbstractMatrix, x) = MomentMatrix(Q, x)

function getmat(μ::MomentMatrix)
    Matrix(μ.Q)
end

function measure(ν::MomentMatrix)
    n = length(ν.x)
    measure(ν.Q.Q, [ν.x[i] * ν.x[j] for i in 1:n for j in 1:i])
end


text:
function make_estim_decorator(breaks=[0.01, 0.05, 0.1], sym='*'; wrapper=identity)
  @assert issorted(breaks)

  function estim_decorator(s, pval)
    pval >= 0 || @error "p value = $pval, but it needs to be non-negative"

    i0 = findfirst(pval .<= breaks)
    i = isnothing(i0) ? length(breaks) + 1 : i0

    deco = sym^(length(breaks) - (i - 1))
    if deco != ""
      deco = wrapper(deco)
    end

    "$s"*deco
  end
end


text:
using QuadGK

sup = (-1,1)
f1(x) = 3/4*(1-x^2)
f2(x) = x < 0 ? x+1 : 1-x

expect(f,support) = quadgk((x) -> x*f(x),support...)[1]

println("Mean 1: ", expect(f1,sup))
println("Mean 2: ", expect(f2,sup))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
p_cutoff = 0.5
pp = p
W = pp[:, 2:end];
W[findall(abs.(W) .< p_cutoff)] .= 0.0
pp[:, 2:end] .= W;

for i_exp in 1:n_exp
    cbi(pp, i_exp)
end
cbp(pp, iter)


text:
"""Experimental package that evaluates pseudoinverses such as `inv(A'A)A'b` lazily, as if `A\b`.
Note that many normal matrix operations will fail, because only a few sample methods have been
defined.
"""
module LazyPseudoinverse

using LinearAlgebra, PDMats
import Base: *, \, promote_op, inv, collect, Matrix, size # overloads
import LinearAlgebra: matprod
import PDMats.PDMat

export MatrixInnerProd, adjointmul, InvMatrixInnerProd, LazyLeftPseudoinverse

abstract type AbstractMatrixInnerProd{T} <: AbstractMatrix{T} end
abstract type AbstractInvMatrixInnerProd{T} <: AbstractMatrix{T} end

"""
`MatrixInnerProd(A::AbstractMatrix)`

Construct a matrix inner product `A'A` from `A`. This stores `A` as a field and
can act like a positive definite matrix. For example, `inv` and `\` are defined,
although `inv` constructs a lazy InvMatrixInnerProd.
Use `collect` to materialize an actual matrix. It can also be converted to
`PDMat(MatrixInnerProd(A))` if `using PDMats`.
"""
struct MatrixInnerProd{T,S<:AbstractMatrix{T}} <: AbstractMatrixInnerProd{T}
    A::S
end

Base.collect(AtA::MatrixInnerProd) = adjointmul(AtA.A',AtA.A)
Base.Matrix(AtA::MatrixInnerProd) = collect(AtA)
Base.size(AtA::MatrixInnerProd) = (size(AtA.A,2),size(AtA.A,2))
# When user shows MatrixInnerProd, print the result (A'A) as if normal, even though it's really lazy
Base.print_array(io::IO, A::MatrixInnerProd) = Base.print_array(io, adjointmul(A.A',A.A))
\(AtA::MatrixInnerProd, B::AbstractVector) = cholesky(collect(AtA)) \ B
\(AtA::MatrixInnerProd, B::AbstractMatrix) = cholesky(collect(AtA)) \ B
*(AtA::MatrixInnerProd, B::AbstractVector) = collect(AtA)*B
*(AtA::MatrixInnerProd, B::AbstractMatrix) = collect(AtA)*B
# Only a few methods defined here, so many matrix operations will fail, e.g. A'A+B

"""
`adjointmul(A', B) performs inner product of `A` and `B`, `A'B`. This is originally `*`
from `matmul.jl` (which has been type-pirated for LazyPseudoinverse.jl). To actually
complete `A'*B`, use `adjointmul(A', B)`
"""
function adjointmul(A::Adjoint,B::AbstractMatrix) # A'*B
    TS = promote_op(matprod, eltype(A), eltype(B)) # this is what * normally does in matmul.jl:151
    mul!(similar(B, TS, (size(A,1), size(B,2))), A, B)
end

"""
`A'*B` is overloaded to perform a normal matrix multiplication, except when fed `A'A` which
returns a `MatrixInnerProd` which stores `A` and acts like positive definite matrix.
Warning: This is type piracy of `Base.*`, be careful!
"""
function *(A::Adjoint,B::AbstractMatrix) # A'*B
    if A.parent === B # might be a bit looser to do ==
        MatrixInnerProd(A.parent)
    else
        adjointmul(A, B) # do a regular matmul
    end
end
PDMats.PDMat(A::MatrixInnerProd) = PDMat(adjointmul(A.A',A.A)) # sample conversion to positive definite matrix

"""
`InvMatrixInnerProd(A::AbstractMatrix)`

Construct an inverse matrix inner product `inv(A'A`) from `A`. This stores `A` as a field and
can act like an inverted positive definite matrix.
Use `collect` to materialize an actual matrix. It can also be converted to
`PDMat(MatrixInnerProd(A))` when `using PDMats`.
"""
struct InvMatrixInnerProd{T,S<:AbstractMatrix{T}} <: AbstractInvMatrixInnerProd{T}
    A::S
end
Base.collect(invAtA::InvMatrixInnerProd) = inv(cholesky(adjointmul(invAtA.A',invAtA.A)))
Base.Matrix(invAtA::InvMatrixInnerProd) = collect(invAtA)
Base.size(invAtA::InvMatrixInnerProd) = (size(invAtA.A,2),size(invAtA.A,2))
# When user shows InvMatrixInnerProd, print the result inv(A'A) as computed regularly:
Base.print_array(io::IO, A::InvMatrixInnerProd) = Base.print_array(io, inv(adjointmul(A.A',A.A)))
PDMats.PDMat(A::InvMatrixInnerProd) = PDMat(inv(adjointmul(A.A',A.A))) # Cholesky decomposition

inv(AtA::MatrixInnerProd) = InvMatrixInnerProd(AtA.A) # lazy inverse

function *(invAtA::InvMatrixInnerProd,B::Adjoint) # inv(A'A)*B'
    if invAtA.A === B.parent  # B = A, so inv(A'A)*A'
        LazyLeftPseudoinverse(invAtA.A) # return the pseudo-inverse, ready to be right-multiplied
    else
        cholesky(adjointmul(invAtA.A',invAtA.A))\B # use linear solve instead of explicit inverse
    end
end

# Ensure that stuff like inv(A'A)B, or B*inv(A'A) still works as normal, although to get
# useful functionality many other methods need to be defined
*(invAtA::InvMatrixInnerProd,B::AbstractVector) = cholesky(adjointmul(invAtA.A',invAtA.A))\B
*(invAtA::InvMatrixInnerProd,B::AbstractMatrix) = cholesky(adjointmul(invAtA.A',invAtA.A))\B
*(B::AbstractMatrix,invAtA::InvMatrixInnerProd) = B*collect(invAtA)

"""
`LazyLeftPseudoinverse(A::AbstractMatrix)`

Construct equivalent of `inv(A'A)A'` from `A`. This stores `A` as a field and
can act (and displays) like a pseudo-inverse `pinv(A)`. But when actually multiplied
against `b`, it computes `A\\b` using QR.
Use `collect` to materialize an actual matrix.
"""
struct LazyLeftPseudoinverse{T,S<:AbstractMatrix{T}} <: AbstractMatrix{T}
    A::S
end

Base.size(pending::LazyLeftPseudoinverse) = reverse(size(pending.A))
Base.print_array(io::IO, pending::LazyLeftPseudoinverse) = Base.print_array(io, pinv(pending.A))
*(pending::LazyLeftPseudoinverse, B::AbstractVector) = pending.A \ B
*(pending::LazyLeftPseudoinverse, B::AbstractMatrix) = pending.A \ B
*(B::AbstractMatrix,pending::LazyLeftPseudoinverse) = B*pinv(pending.A)
Base.collect(lazyLeftpinv::LazyLeftPseudoinverse) = pinv(lazyLeftpinv.A)
Base.Matrix(lazyLeftpinv::LazyLeftPseudoinverse) = collect(lazyLeftpinv)

end # module


text:
###
### **acb_dirichlet.h** -- Dirichlet L-functions, Riemann zeta and related functions
###

### Roots of unity
#ni arbcall"void acb_dirichlet_roots_init(acb_dirichlet_roots_t roots, ulong n, slong num, slong prec)"
#ni arbcall"void acb_dirichlet_roots_clear(acb_dirichlet_roots_t roots)"
#ni arbcall"void acb_dirichlet_root(acb_t res, const acb_dirichlet_roots_t roots, ulong k, slong prec)"

### Truncated L-series and power sums
arbcall"void acb_dirichlet_powsum_term(acb_ptr res, arb_t log_prev, ulong * prev, const acb_t s, ulong k, int integer, int critical_line, slong len, slong prec)"
arbcall"void acb_dirichlet_powsum_sieved(acb_ptr res, const acb_t s, ulong n, slong len, slong prec)"
arbcall"void acb_dirichlet_powsum_smooth(acb_ptr res, const acb_t s, ulong n, slong len, slong prec)"

### Riemann zeta function
arbcall"void acb_dirichlet_zeta(acb_t res, const acb_t s, slong prec)"
arbcall"void acb_dirichlet_zeta_jet(acb_t res, const acb_t s, int deflate, slong len, slong prec)"
arbcall"void acb_dirichlet_zeta_bound(mag_t res, const acb_t s)"
arbcall"void acb_dirichlet_zeta_deriv_bound(mag_t der1, mag_t der2, const acb_t s)"
arbcall"void acb_dirichlet_eta(acb_t res, const acb_t s, slong prec)"
arbcall"void acb_dirichlet_xi(acb_t res, const acb_t s, slong prec)"

### Riemann-Siegel formula
arbcall"void acb_dirichlet_zeta_rs_f_coeffs(acb_ptr f, const arb_t p, slong n, slong prec)"
arbcall"void acb_dirichlet_zeta_rs_d_coeffs(arb_ptr d, const arb_t sigma, slong k, slong prec)"
arbcall"void acb_dirichlet_zeta_rs_bound(mag_t err, const acb_t s, slong K)"
arbcall"void acb_dirichlet_zeta_rs_r(acb_t res, const acb_t s, slong K, slong prec)"
arbcall"void acb_dirichlet_zeta_rs(acb_t res, const acb_t s, slong K, slong prec)"
arbcall"void acb_dirichlet_zeta_jet_rs(acb_t res, const acb_t s, slong len, slong prec)"

### Hurwitz zeta function
arbcall"void acb_dirichlet_hurwitz(acb_t res, const acb_t s, const acb_t a, slong prec)"

### Hurwitz zeta function precomputation
#ni arbcall"void acb_dirichlet_hurwitz_precomp_init(acb_dirichlet_hurwitz_precomp_t pre, const acb_t s, int deflate, ulong A, ulong K, ulong N, slong prec)"
#ni arbcall"void acb_dirichlet_hurwitz_precomp_init_num(acb_dirichlet_hurwitz_precomp_t pre, const acb_t s, int deflate, double num_eval, slong prec)"
#ni arbcall"void acb_dirichlet_hurwitz_precomp_clear(acb_dirichlet_hurwitz_precomp_t pre)"
arbcall"void acb_dirichler_hurwitz_precomp_choose_param(ulong * A, ulong * K, ulong * N, const acb_t s, double num_eval, slong prec)"
arbcall"void acb_dirichlet_hurwitz_precomp_bound(mag_t res, const acb_t s, ulong A, ulong K, ulong N)"
#ni arbcall"void acb_dirichlet_hurwitz_precomp_eval(acb_t res, const acb_dirichlet_hurwitz_precomp_t pre, ulong p, ulong q, slong prec)"

### Stieltjes constants
#ni arbcall"void acb_dirichlet_stieltjes(acb_t res, const fmpz_t n, const acb_t a, slong prec)"

### Dirichlet character evaluation
#ni arbcall"void acb_dirichlet_chi(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, ulong n, slong prec)"
#ni arbcall"void acb_dirichlet_chi_vec(acb_ptr v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv, slong prec)"
#ni arbcall"void acb_dirichlet_pairing(acb_t res, const dirichlet_group_t G, ulong m, ulong n, slong prec)"
#ni arbcall"void acb_dirichlet_pairing_char(acb_t res, const dirichlet_group_t G, const dirichlet_char_t a, const dirichlet_char_t b, slong prec)"

### Dirichlet character Gauss, Jacobi and theta sums
#ni arbcall"void acb_dirichlet_gauss_sum_naive(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)"
#ni arbcall"void acb_dirichlet_gauss_sum_factor(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)"
#ni arbcall"void acb_dirichlet_gauss_sum_order2(acb_t res, const dirichlet_char_t chi, slong prec)"
#ni arbcall"void acb_dirichlet_gauss_sum_theta(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)"
#ni arbcall"void acb_dirichlet_gauss_sum(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)"
#ni arbcall"void acb_dirichlet_gauss_sum_ui(acb_t res, const dirichlet_group_t G, ulong a, slong prec)"
#ni arbcall"void acb_dirichlet_jacobi_sum_naive(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi1, const dirichlet_char_t chi2, slong prec)"
#ni arbcall"void acb_dirichlet_jacobi_sum_factor(acb_t res,  const dirichlet_group_t G, const dirichlet_char_t chi1, const dirichlet_char_t chi2, slong prec)"
#ni arbcall"void acb_dirichlet_jacobi_sum_gauss(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi1, const dirichlet_char_t chi2, slong prec)"
#ni arbcall"void acb_dirichlet_jacobi_sum(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi1,  const dirichlet_char_t chi2, slong prec)"
#ni arbcall"void acb_dirichlet_jacobi_sum_ui(acb_t res, const dirichlet_group_t G, ulong a, ulong b, slong prec)"
#ni arbcall"void acb_dirichlet_chi_theta_arb(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, const arb_t t, slong prec)"
#ni arbcall"void acb_dirichlet_ui_theta_arb(acb_t res, const dirichlet_group_t G, ulong a, const arb_t t, slong prec)"
arbcall"ulong acb_dirichlet_theta_length(ulong q, const arb_t t, slong prec)"
#ni arbcall"void acb_dirichlet_qseries_powers_naive(acb_t res, const arb_t x, int p, const ulong * a, const acb_dirichlet_powers_t z, slong len, slong prec)"
#ni arbcall"void acb_dirichlet_qseries_powers_smallorder(acb_t res, const arb_t x, int p, const ulong * a, const acb_dirichlet_powers_t z, slong len, slong prec)"

### Discrete Fourier transforms
#ni arbcall"void acb_dirichlet_dft_conrey(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec)"
#ni arbcall"void acb_dirichlet_dft(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec)"

### Dirichlet L-functions
#ni arbcall"void acb_dirichlet_root_number_theta(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)"
#ni arbcall"void acb_dirichlet_root_number(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)"
#ni arbcall"void acb_dirichlet_l_hurwitz(acb_t res, const acb_t s, const acb_dirichlet_hurwitz_precomp_t precomp, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)"
#ni arbcall"void acb_dirichlet_l_euler_product(acb_t res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)"
#ni arbcall"void _acb_dirichlet_euler_product_real_ui(arb_t res, ulong s, const signed char * chi, int mod, int reciprocal, slong prec)"
#ni arbcall"void acb_dirichlet_l(acb_t res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)"
#ni arbcall"void acb_dirichlet_l_vec_hurwitz(acb_ptr res, const acb_t s, const acb_dirichlet_hurwitz_precomp_t precomp, const dirichlet_group_t G, slong prec)"
#ni arbcall"void acb_dirichlet_l_jet(acb_ptr res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, int deflate, slong len, slong prec)"
#ni arbcall"void _acb_dirichlet_l_series(acb_ptr res, acb_srcptr s, slong slen, const dirichlet_group_t G, const dirichlet_char_t chi, int deflate, slong len, slong prec)"
#ni arbcall"void acb_dirichlet_l_series(acb_poly_t res, const acb_poly_t s, const dirichlet_group_t G, const dirichlet_char_t chi, int deflate, slong len, slong prec)"

### Hardy Z-functions
#ni arbcall"void acb_dirichlet_hardy_theta(acb_ptr res, const acb_t t, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)"
#ni arbcall"void acb_dirichlet_hardy_z(acb_t res, const acb_t t, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)"
#ni arbcall"void _acb_dirichlet_hardy_theta_series(acb_ptr res, acb_srcptr t, slong tlen, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)"
#ni arbcall"void acb_dirichlet_hardy_theta_series(acb_poly_t res, const acb_poly_t t, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)"
#ni arbcall"void _acb_dirichlet_hardy_z_series(acb_ptr res, acb_srcptr t, slong tlen, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)"
#ni arbcall"void acb_dirichlet_hardy_z_series(acb_poly_t res, const acb_poly_t t, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)"

### Gram points
#ni arbcall"void acb_dirichlet_gram_point(arb_t res, const fmpz_t n, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)"

### Riemann zeta function zeros
#ni arbcall"ulong acb_dirichlet_turing_method_bound(const fmpz_t p)"
arbcall"int _acb_dirichlet_definite_hardy_z(arb_t res, const arf_t t, slong * pprec)"
#ni arbcall"void _acb_dirichlet_isolate_gram_hardy_z_zero(arf_t a, arf_t b, const fmpz_t n)"
#ni arbcall"void _acb_dirichlet_isolate_rosser_hardy_z_zero(arf_t a, arf_t b, const fmpz_t n)"
#ni arbcall"void _acb_dirichlet_isolate_turing_hardy_z_zero(arf_t a, arf_t b, const fmpz_t n)"
#ni arbcall"void acb_dirichlet_isolate_hardy_z_zero(arf_t a, arf_t b, const fmpz_t n)"
arbcall"void _acb_dirichlet_refine_hardy_z_zero(arb_t res, const arf_t a, const arf_t b, slong prec)"
#ni arbcall"void acb_dirichlet_hardy_z_zero(arb_t res, const fmpz_t n, slong prec)"
#ni arbcall"void acb_dirichlet_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, slong prec)"
#ni arbcall"void acb_dirichlet_zeta_zero(acb_t res, const fmpz_t n, slong prec)"
#ni arbcall"void acb_dirichlet_zeta_zeros(acb_ptr res, const fmpz_t n, slong len, slong prec)"
#ni arbcall"void _acb_dirichlet_exact_zeta_nzeros(fmpz_t res, const arf_t t)"
arbcall"void acb_dirichlet_zeta_nzeros(arb_t res, const arb_t t, slong prec)"
arbcall"void acb_dirichlet_backlund_s(arb_t res, const arb_t t, slong prec)"
arbcall"void acb_dirichlet_backlund_s_bound(mag_t res, const arb_t t)"
#ni arbcall"void acb_dirichlet_zeta_nzeros_gram(fmpz_t res, const fmpz_t n)"
#ni arbcall"slong acb_dirichlet_backlund_s_gram(const fmpz_t n)"

### Riemann zeta function zeros (Platt's method)
arbcall"void acb_dirichlet_platt_scaled_lambda(arb_t res, const arb_t t, slong prec)"
#ni arbcall"void acb_dirichlet_platt_scaled_lambda_vec(arb_ptr res, const fmpz_t T, slong A, slong B, slong prec)"
#ni arbcall"void acb_dirichlet_platt_multieval(arb_ptr res, const fmpz_t T, slong A, slong B, const arb_t h, slong J, slong K, slong sigma, slong prec)"
#ni arbcall"void acb_dirichlet_platt_multieval_threaded(arb_ptr res, const fmpz_t T, slong A, slong B, const arb_t h, slong J, slong K, slong sigma, slong prec)"
#ni arbcall"void acb_dirichlet_platt_ws_interpolation(arb_t res, arf_t deriv, const arb_t t0, arb_srcptr p, const fmpz_t T, slong A, slong B, slong Ns_max, const arb_t H, slong sigma, slong prec)"
#ni arbcall"slong _acb_dirichlet_platt_local_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, const fmpz_t T, slong A, slong B, const arb_t h, slong J, slong K, slong sigma_grid, slong Ns_max, const arb_t H, slong sigma_interp, slong prec)"
#ni arbcall"slong acb_dirichlet_platt_local_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, slong prec)"
#ni arbcall"slong acb_dirichlet_platt_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, slong prec)"
#ni arbcall"slong acb_dirichlet_platt_zeta_zeros(acb_ptr res, const fmpz_t n, slong len, slong prec)"


text:
# Implementations of CGD for linear games
using LinearAlgebra
using StaticArrays
using LinearMaps
using IterativeSolvers
include("AutoDiffUtils.jl")

function testinv!(x, A, b)
  x .= Matrix(A)\b
end

# Struct containing the information for cmd
struct CMD
    # First player
    # ############################################################
    # The present optimization variable for x
    x̄

    # An initial guess to use when solving for x
    x̄_guess

    # The dual potential's derivative that is used to map a dual variable to its primal
    Dxψ

    # The dual potential's Hessian
    Dxxψ

    # The inverse of the dual potential's Hessian, which is used precondition the conjugate gradient solver
    Dxxψ_inv

    # Gradient of x's loss function
    ∇_xf

    # Mixed Hessian of x's loss function
    D_xyf

    # Second player
    # ############################################################
    # The present optimization variable for y
    ȳ

    # An initial guess to use when solving for y
    ȳ_guess

    # The dual potential's derivative that is used to map a dual variable to its primal
    Dyψ

    # The dual potential's Hessian
    Dyyψ

    # The inverse of the dual potential's Hessian, which is used precondition the conjugate gradient solver
    Dyyψ_inv

    # Gradient of ys loss function
    ∇_yg

    # Mixed Hessian of ys loss function
    D_yxg

    inv_x::Vector{Bool}

    algorithm!
end

function CMD(x, xPotential, f, y, yPotential, g, algorithm=gmres!)
  return CMD(x, zero(x), xPotential[1], xPotential[3], xPotential[4], ∇_x(f), D_xy(f), y, zero(y), yPotential[1], yPotential[3], yPotential[4], ∇_y(g), D_yx(g), trues(1), algorithm)
end

# preconditioner class
struct Prec
  inverse
end

import LinearAlgebra.ldiv!
function ldiv!(y, P::Prec, x)
  y .= P.inverse * x
end

function ldiv!(P::Prec, x)
  return P.inverse * x
end

function step!(cmd::CMD)
  x = cmd.Dxψ(cmd.x̄)
  Dxxψ = cmd.Dxxψ(x)
  Dxxψ_inv = cmd.Dxxψ_inv(x)

  y = cmd.Dyψ(cmd.ȳ)
  Dyyψ = cmd.Dyyψ(y)
  Dyyψ_inv = cmd.Dyyψ_inv(y)

  function D_xyf!(u, v)
    u .= cmd.D_xyf(x, y, v)
  end

  function D_yxg!(u, v)
    u .= cmd.D_yxg(x, y, v)
  end

  ∇_xf = cmd.∇_xf(x, y)
  ∇_yg = cmd.∇_yg(x, y)

  # wrapping the mutation hvp into linear maps
  D_xyf = LinearMap{eltype(cmd.ȳ)}(D_xyf!, length(cmd.x̄), length(cmd.ȳ); ismutating=true, issymmetric=false, ishermitian=false)

  D_yxg = LinearMap{eltype(cmd.x̄)}(D_yxg!, length(cmd.ȳ), length(cmd.x̄); ismutating=true, issymmetric=false, ishermitian=false)

  # For the exact solve

  # cmd.x̄ .+= - Matrix(Dxxψ - Dxxψ * D_xyf * Dyyψ * D_yxg * Dxxψ) \ (Dxxψ * ∇_xf - Dxxψ * D_xyf * Dyyψ * ∇_yg)

  # cmd.ȳ .+= - Matrix(Dyyψ - Dyyψ * D_yxg * Dxxψ * D_xyf * Dyyψ) \ (Dyyψ * ∇_yg - Dyyψ * D_yxg * Dxxψ *∇_xf)

  # For timing comparison

  if cmd.inv_x[1] == true
    P = Prec(Dxxψ_inv)
    # P = cholesky(Matrix(Dxxψ))
    hist = cmd.algorithm!(cmd.x̄_guess, Dxxψ - Dxxψ * D_xyf * Dyyψ * D_yxg * Dxxψ, -(Dxxψ * ∇_xf - Dxxψ * D_xyf * Dyyψ * ∇_yg), log=true,
    tol=1e-9,
    Pl=P,
    )[2]
    cmd.ȳ_guess .= - ∇_yg - D_yxg * Dxxψ * cmd.x̄_guess

    cmd.inv_x[1] = false
  else
    P = Prec(Dyyψ_inv)
    # P = cholesky(Matrix(Dyyψ))
    hist = cmd.algorithm!(cmd.ȳ_guess, Dyyψ - Dyyψ * D_yxg * Dxxψ * D_xyf * Dyyψ, -(Dyyψ * ∇_yg - Dyyψ * D_yxg * Dxxψ *∇_xf), log=true,
    tol=1e-9,
    Pl=P,
    )[2]
    cmd.x̄_guess .= - ∇_xf - D_xyf * Dyyψ * cmd.ȳ_guess

    cmd.inv_x[1] = true
  end

  cmd.x̄ .+= cmd.x̄_guess
  cmd.ȳ .+= cmd.ȳ_guess

  return 4 + 2 * hist.iters
end

function iterate_primal(cmd::CMD)
  return cmd.Dxψ(cmd.x̄), cmd.Dyψ(cmd.ȳ)
end

function iterate_dual(cmd::CMD)
  return copy(cmd.x̄), copy(cmd.ȳ)
end

max_stars_repo_name: f-t-s/CMD

max_issues_repo_name: f-t-s/CMD

max_forks_repo_name: f-t-s/CMD


text:
abstract type AbstractZXDiagram{T, P} end

LightGraphs.nv(zxd::AbstractZXDiagram) = throw(MethodError(LightGraphs.nv, zxd))
LightGraphs.ne(zxd::AbstractZXDiagram) = throw(MethodError(LightGraphs.ne, zxd))
LightGraphs.degree(zxd::AbstractZXDiagram, v) = throw(MethodError(LightGraphs.degree, (zxd, v)))
LightGraphs.indegree(zxd::AbstractZXDiagram, v) = throw(MethodError(LightGraphs.indegree, (zxd, v)))
LightGraphs.outdegree(zxd::AbstractZXDiagram, v) = throw(MethodError(LightGraphs.outdegree, (zxd, v)))
LightGraphs.neighbors(zxd::AbstractZXDiagram, v) = throw(MethodError(LightGraphs.neighbors, (zxd, v)))
LightGraphs.outneighbors(zxd::AbstractZXDiagram, v) = throw(MethodError(LightGraphs.outneighbors, (zxd, v)))
LightGraphs.inneighbors(zxd::AbstractZXDiagram, v) = throw(MethodError(LightGraphs.inneighbors, (zxd, v)))
LightGraphs.rem_edge!(zxd::AbstractZXDiagram, args...) = throw(MethodError(LightGraphs.rem_edge!, (zxd, args...)))
LightGraphs.add_edge!(zxd::AbstractZXDiagram, args...) = throw(MethodError(LightGraphs.add_edge!, (zxd, args...)))
LightGraphs.has_edge(zxd::AbstractZXDiagram, args...) = throw(MethodError(LightGraphs.has_edge, (zxd, args...)))

Base.show(io::IO, zxd::AbstractZXDiagram) = throw(MethodError(Base.show, io, zxd))
Base.copy(zxd::AbstractZXDiagram) = throw(MethodError(Base.copy, zxd))

nqubits(zxd::AbstractZXDiagram) = throw(MethodError(ZXCalculus.nqubits, zxd))
spiders(zxd::AbstractZXDiagram) = throw(MethodError(ZXCalculus.spiders, zxd))
tcount(zxd::AbstractZXDiagram) = throw(MethodError(ZXCalculus.tcount, zxd))
get_inputs(zxd::AbstractZXDiagram) = throw(MethodError(ZXCalculus.get_inputs, zxd))
get_outputs(zxd::AbstractZXDiagram) = throw(MethodError(ZXCalculus.get_outputs, zxd))
scalar(zxd::AbstractZXDiagram) = throw(MethodError(ZXCalculus.scalar, zxd))
spider_sequence(zxd::AbstractZXDiagram) = throw(MethodError(ZXCalculus.spider_sequence, zxd))
round_phases!(zxd::AbstractZXDiagram) = throw(MethodError(ZXCalculus.round_phases!, zxd))

spider_type(zxd::AbstractZXDiagram, v) = throw(MethodError(ZXCalculus.spider_type, (zxd, v)))
phase(zxd::AbstractZXDiagram, v) = throw(MethodError(ZXCalculus.phase, (zxd, v)))
rem_spider!(zxd::AbstractZXDiagram, v) = throw(MethodError(ZXCalculus.rem_spider!, (zxd, v)))
rem_spiders!(zxd::AbstractZXDiagram, vs)= throw(MethodError(ZXCalculus.rem_spiders!, (zxd, vs)))
qubit_loc(zxd::AbstractZXDiagram, v) = throw(MethodError(ZXCalculus.qubit_loc, (zxd, v)))
column_loc(zxd::AbstractZXDiagram, v) = throw(MethodError(ZXCalculus.column_loc, (zxd, v)))
add_global_phase!(zxd::AbstractZXDiagram, p) = throw(MethodError(ZXCalculus.add_global_phase!, (zxd, p)))
add_power!(zxd::AbstractZXDiagram, n) = throw(MethodError(ZXCalculus.add_power!, (zxd, n)))
generate_layout!(zxd::AbstractZXDiagram, seq) = throw(MethodError(ZXCalculus.generate_layout!, (zxd, seq)))

set_phase!(zxd::AbstractZXDiagram, args...) = throw(MethodError(ZXCalculus.set_phase!, (zxd, args...)))
push_gate!(zxd::AbstractZXDiagram, args...) = throw(MethodError(ZXCalculus.push_gate!, (zxd, args...)))
pushfirst_gate!(zxd::AbstractZXDiagram, args...) = throw(MethodError(ZXCalculus.pushfirst_gate!, (zxd, args...)))
add_spider!(zxd::AbstractZXDiagram, args...) = throw(MethodError(ZXCalculus.add_spider!, (zxd, args...)))
insert_spider!(zxd::AbstractZXDiagram, args...) = throw(MethodError(ZXCalculus.insert_spider!, (zxd, args...)))


text:
using RMInsurance

using Distributions
using DataFrames
using Test

include("S2NonLife.jl")

println("Start S2NonLife test")

@test round(scr_prem_res, digits = 2) ≈ 210.95
@test round(scr_lapse, digits = 2) ≈ 0.0
@test round(scr_cat_fire, digits = 2) ≈ 60.00
@test round(scr_cat_liab, digits = 2) ≈ 239.65
@test round(scr_cat_man_made, digits = 2) ≈ 247.05
@test round(scr_cat_other, digits = 2) ≈ 0.0
@test round(scr_cat_nprop, digits = 2) ≈ 0.0
@test round(scr_cat, digits = 2) ≈ 247.05
@test round(scr_nl, digits = 2) ≈ 362.76

println("End S2NonLife test")


text:
###############################################################################
#
# Stopping can also be used for fixed point methods
# Example here concerns the AlternatingDirections Algorithm to find
# a feasible point in the intersection of 2 convex sets A and B.
# This algorithm relies on a fixed point argument, hence it stopped if it finds
# a fixed point.
#
# Example:
# A={ (x,y) | x=y} and B = {(x,y) | y=0}
# Clearly the unique intersection point is (0,0)
#
# Note that in this case the projection on A and the projection on B are trivial
#
# Takeaway: the 2nd scenario illustrates a situation where the algorithm stalls
# as it reached a personal success. (optimal_sub_pb is true)
#
###############################################################################
using LinearAlgebra, NLPModels, Stopping, Test

#Main algorithm
function AlternatingDirections(stp)

 xk = stp.current_state.x
 OK = update_and_start!(stp, cx = cons(stp.pb, x0))
 @show OK, xk

 while !OK

  #First projection
  xk1 = 0.5 * (xk[1] + xk[2]) * ones(2)
  #Second projection
  xk2 = [xk1[1],0.0]

  #check if we have a fixed point
  Fix = dot(xk-xk2,xk-xk2)
  if Fix <= min(eps(Float64),stp.meta.atol) stp.meta.suboptimal = true end
  #call the stopping
  OK = update_and_stop!(stp, x = xk2, cx = cons(stp.pb, xk2))

  xk = xk2
  @show OK, xk
 end

 return stp
end

# We model the problem using the NLPModels without objective function
#Formulate the problem with NLPModels
c(x) = [x[1] - x[2], x[2]]
lcon = [0.0, 0.0]
ucon = [0.0, 0.0]
nlp = ADNLPModel(x->0.0, zeros(2), c=c, lcon=lcon, ucon=ucon)

#1st scenario: we solve the problem
printstyled("1st scenario:\n")
#Prepare the Stopping
x0 = [0.0, 5.0]
state = NLPAtX(x0)
#Recall that for the optimality_check function x is the pb and y is the state
#Here we take the infinite norm of the residual.
stop = NLPStopping(nlp, state, optimality_check = (x,y) -> norm(y.cx,Inf))

AlternatingDirections(stop)
@show status(stop)
@test status(stop) == :Optimal

#2nd scenario: the user gives an irrealistic optimality condition
printstyled("2nd scenario:\n")
reinit!(stop, rstate = true, x = x0)
stop.meta.optimality_check = (x,y) -> norm(y.cx,Inf)+0.5

AlternatingDirections(stop)
#In this scenario, the algorithm stops because it attains a fixed point
#Hence, status is :SubOptimal.
@show status(stop)
@test status(stop) == :SubOptimal


text:
function demoPotLogArray()

    M1=rand(2,2)
    M2=rand(2,2)
    p1=PotArray([1 2],M1)
    p2=PotArray([2 3],M2)
    lp1=PotLogArray([1 2],log(M1))
    lp2=PotLogArray([2 3],log(M2))
    lA=FactorGraph([lp1 lp2])
    A=FactorGraph([p1 p2])

    marg,mess,normconst=sumprod([p1 p2],A)
    lmarg,lmess,lnormconst=sumprod([lp1 lp2],lA)

   println("Marginal table computed using Factor Graph running on standard PotArray, versus marginal table computed using Factor Graph running on LogPotArray:\n")

    println("log(PotArray) \t\t PotLogArray")
    for i=1:length(marg)
        println("\nvariable $i :")
        println("$([log(marg[i].content) lmarg[i].content])")
    end

end


text:
using LinearAlgebra
using Printf

struct RankInfo
    label::AbstractString
    count::Int
    score_values::Array{Int}
    reward::Int
    penalty::Int
end

mutable struct ScoreTransitionMatrix
    transition_matrix::Array{Float64, 2}
    rank_info::RankInfo
end

mutable struct ScoreSteadyState
    probabilities::Array{Float64}
    rank_info::RankInfo
end

function steady_state(transition_matrix::Array{Float64, 2})
    values = eigvals(transition_matrix')
    vectors = eigvecs(transition_matrix')
    for i=1:length(values)
        if abs(values[i] - 1.0) < 1e-9
            vec = real(vectors[:,i])
            return vec/sum(vec)
        end
    end
    return zeros(length(values))
end

function score_steady_state(score_transition_matrix::ScoreTransitionMatrix)
    steady_state_vector = steady_state(score_transition_matrix.transition_matrix)
    return ScoreSteadyState(steady_state_vector, score_transition_matrix.rank_info)
end

function rank_matrix_from_steady_states(steady_states::Array{ScoreSteadyState})
    count = sum([steady_states[i].rank_info.count for i=1:length(steady_states)])

    rank_matrix = Matrix{Float64}(I, count, count)

    i = 1
    for idx=1:length(steady_states)
        rank_info = steady_states[idx].rank_info
        probabilities = steady_states[idx].probabilities
        for offset=1:rank_info.count
            if i > 1
                for j=1:length(rank_info.score_values)
                    if rank_info.score_values[j] + rank_info.penalty < rank_info.score_values[1]
                        rank_matrix[i,i-1] += 0.5 * probabilities[j]
                    end
                end
                rank_matrix[i,i] -= rank_matrix[i,i-1]
            end
            if i < count
                for j=1:length(rank_info.score_values)
                    if rank_info.score_values[j] + rank_info.reward >= 100
                        rank_matrix[i,i+1] += 0.5 * probabilities[j]
                    end
                end
                rank_matrix[i,i] -= rank_matrix[i,i+1]
            end
            i = i + 1
        end
    end

    return rank_matrix
end

function report_score_steady_states(score_steady_states::Array{ScoreSteadyState})
    for i=1:length(score_steady_states)
        ss = score_steady_states[i]
        println(ss.rank_info.label, ": ")
        for j=1:length(ss.probabilities)
            println("  ", ss.rank_info.score_values[j], " => ", ss.probabilities[j])
        end
    end
end

function report_rank_steady_state(rank_steady_state)
    println("Rank steady state: ", rank_steady_state)

    @printf("C: %.2f%%\n", 100 * sum(rank_steady_state[1:3]))
    @printf("B: %.2f%%\n", 100 * sum(rank_steady_state[4:6]))
    @printf("A: %.2f%%\n", 100 * sum(rank_steady_state[7:9]))
    @printf("S: %.2f%%\n", 100 * rank_steady_state[10])
    @printf("S+: %.2f%%\n", 100 * rank_steady_state[11])
end

function default_score_transitions()
                            # 0  10  20  30  40  50  60  70  80  90
    c_minus_matrix = Float64[0.5  0  0.5  0   0   0   0   0   0   0; #  0
                             0.5  0   0  0.5  0   0   0   0   0   0; # 10
                              0  0.5  0   0  0.5  0   0   0   0   0; # 20
                              0   0  0.5  0   0  0.5  0   0   0   0; # 30
                              0   0   0  0.5  0   0  0.5  0   0   0; # 40
                              0   0   0   0  0.5  0   0  0.5  0   0; # 50
                              0   0   0   0   0  0.5  0   0  0.5  0; # 60
                              0   0   0   0   0   0  0.5  0   0  0.5;# 70
                              0   0   0   0   0   0   0  1.0  0   0; # 80
                              0   0   0   0   0   0   0  0.5 0.5  0 ]# 90

    c_minus = ScoreTransitionMatrix(c_minus_matrix, RankInfo("C-", 1, [10*i for i=0:9], 20, -10))

                      # 0   5  10   15  20  25  30  35  40  45  50  55  60  65  70  75  80  85  90  95
    c_matrix = Float64[ 0   0   0   0.5  0   0  0.5  0   0   0   0   0   0   0   0   0   0   0   0   0; #  0
                        0   0   0    0  0.5  0  0.5  0   0   0   0   0   0   0   0   0   0   0   0   0; #  5
                       0.5  0   0    0   0  0.5  0   0   0   0   0   0   0   0   0   0   0   0   0   0; # 10
                        0  0.5  0    0   0   0  0.5  0   0   0   0   0   0   0   0   0   0   0   0   0; # 15
                        0   0  0.5   0   0   0   0  0.5  0   0   0   0   0   0   0   0   0   0   0   0; # 20
                        0   0   0   0.5  0   0   0   0  0.5  0   0   0   0   0   0   0   0   0   0   0; # 25
                        0   0   0    0  0.5  0   0   0   0  0.5  0   0   0   0   0   0   0   0   0   0; # 30
                        0   0   0    0   0  0.5  0   0   0   0  0.5  0   0   0   0   0   0   0   0   0; # 35
                        0   0   0    0   0   0  0.5  0   0   0   0  0.5  0   0   0   0   0   0   0   0; # 40
                        0   0   0    0   0   0   0  0.5  0   0   0   0  0.5  0   0   0   0   0   0   0; # 45
                        0   0   0    0   0   0   0   0  0.5  0   0   0   0  0.5  0   0   0   0   0   0; # 50
                        0   0   0    0   0   0   0   0   0  0.5  0   0   0   0  0.5  0   0   0   0   0; # 55
                        0   0   0    0   0   0   0   0   0   0  0.5  0   0   0   0  0.5  0   0   0   0; # 60
                        0   0   0    0   0   0   0   0   0   0   0  0.5  0   0   0   0  0.5  0   0   0; # 65
                        0   0   0    0   0   0   0   0   0   0   0   0  0.5  0   0   0   0  0.5  0   0; # 70
                        0   0   0    0   0   0   0   0   0   0   0   0   0  0.5  0   0   0   0  0.5  0; # 75
                        0   0   0    0   0   0   0   0   0   0   0   0   0   0  0.5  0   0   0   0  0.5; # 80
                        0   0   0    0   0   0   0   0   0   0   0   0   0   0  0.5 0.5  0   0   0   0; # 85
                        0   0   0    0   0   0   0   0   0   0   0   0   0   0  0.5  0  0.5  0   0   0; # 90
                        0   0   0    0   0   0   0   0   0   0   0   0   0   0  0.5  0   0  0.5  0   0] # 95

    c = ScoreTransitionMatrix(c_matrix, RankInfo("C", 1, [5*i for i=0:19], 15, -10))

    c_plus_matrix = zeros(50, 50)

    for i=1:50
        if i+6 <= 50
            c_plus_matrix[i,i+6] = 0.5 # +12
        else
            c_plus_matrix[i,36] = 0.5 # 70
        end

        if i-5 >= 1
            c_plus_matrix[i,i-5] = 0.5 # -10
        else
            c_plus_matrix[i,16] = 0.5 # 30
        end
    end

    c_plus_and_b_minus = ScoreTransitionMatrix(c_plus_matrix, RankInfo("C+/B-", 2, [2*i for i=0:49], 12, -10))

                      # 0  10  20  30  40  50  60  70  80  90
    b_matrix = Float64[ 0  0.5  0  0.5  0   0   0   0   0   0; #  0
                       0.5  0  0.5  0   0   0   0   0   0   0; # 10
                        0  0.5  0  0.5  0   0   0   0   0   0; # 20
                        0   0  0.5  0  0.5  0   0   0   0   0; # 30
                        0   0   0  0.5  0  0.5  0   0   0   0; # 40
                        0   0   0   0  0.5  0  0.5  0   0   0; # 50
                        0   0   0   0   0  0.5  0  0.5  0   0; # 60
                        0   0   0   0   0   0  0.5  0  0.5  0; # 70
                        0   0   0   0   0   0   0  0.5  0  0.5;# 80
                        0   0   0   0   0   0   0  0.5 0.5  0];# 90

    b_thru_a_plus = ScoreTransitionMatrix(b_matrix, RankInfo("B/B+/A-/A/A+", 5, [10*i for i=0:9], 10, -10))

    s_matrix = zeros(100, 100)

    for i=1:100
        if i <= 40
            reward = 5
        else
            reward = 4
        end
        if i <= 80
            penalty = 5
        else
            penalty = 6
        end

        if i+reward <= 100
            s_matrix[i,i+reward] = 0.5
        else
            s_matrix[i,71] = 0.5
        end

        if i-penalty >= 1
            s_matrix[i,i-penalty] = 0.5
        else
            s_matrix[i,31] = 0.5
        end
    end

    s = ScoreTransitionMatrix(s_matrix, RankInfo("S", 1, [i for i=0:99], 4, -5))

    s_plus_matrix = zeros(100, 100)

    for i=1:100
        if i <= 40
            reward = 4
        elseif i <= 80
            reward = 3
        else
            reward = 2
        end
        if i <= 40
            penalty = 4
        else
            penalty = 5
        end

        if i+reward <= 100
            s_plus_matrix[i,i+reward] = 0.5
        else
            s_plus_matrix[i,100] = 0.5
        end

        if i-penalty >= 1
            s_plus_matrix[i,i-penalty] = 0.5
        else
            s_plus_matrix[i,31] = 0.5
        end
    end

    s_plus = ScoreTransitionMatrix(s_plus_matrix, RankInfo("S+", 1, [i for i=0:99], 2, -4))

#                          # 0  10  20  30  40  50  60  70  80  90 100
#   a_plus_matrix = Float64[ 0  0.5  0  0.5  0   0   0   0   0   0   0; #  0
#                           0.5  0  0.5  0   0   0   0   0   0   0   0; # 10
#                            0  0.5  0  0.5  0   0   0   0   0   0   0; # 20
#                            0   0  0.5  0  0.5  0   0   0   0   0   0; # 30
#                            0   0   0  0.5  0  0.5  0   0   0   0   0; # 40
#                            0   0   0   0  0.5  0  0.5  0   0   0   0; # 50
#                            0   0   0   0   0  0.5  0  0.5  0   0   0; # 60
#                            0   0   0   0   0   0  0.5  0  0.5  0   0; # 70
#                            0   0   0   0   0   0   0  0.5  0  0.5  0;# 80
#                            0   0   0   0   0   0   0   0  0.5  0  0.5;# 90
#                            0   0   0   0   0   0   0   0   0  0.5 0.5];# 100

#   a_plus = ScoreTransitionMatrix(a_plus_matrix, RankInfo("A+", 1, [10*i for i in 0:10], 10, -10))

    return ScoreTransitionMatrix[ c_minus, c, c_plus_and_b_minus, b_thru_a_plus, s, s_plus ]
end

function proposed_score_transitions()
                            # 0  10  20  30  40  50  60  70  80  90
    c_minus_matrix = Float64[0.5  0  0.5  0   0   0   0   0   0   0; #  0
                             0.5  0   0  0.5  0   0   0   0   0   0; # 10
                              0  0.5  0   0  0.5  0   0   0   0   0; # 20
                              0   0  0.5  0   0  0.5  0   0   0   0; # 30
                              0   0   0  0.5  0   0  0.5  0   0   0; # 40
                              0   0   0   0  0.5  0   0  0.5  0   0; # 50
                              0   0   0   0   0  0.5  0   0  0.5  0; # 60
                              0   0   0   0   0   0  0.5  0   0  0.5;# 70
                              0   0   0   0   0   0   0  1.0  0   0; # 80
                              0   0   0   0   0   0   0  0.5 0.5  0 ]# 90

    c_minus = ScoreTransitionMatrix(c_minus_matrix, RankInfo("C-", 1, [10*i for i=0:9], 20, -10))

                      #  5  10   15  20  25  30  35  40  45  50  55  60  65  70  75  80  85  90  95
    c_matrix0= Float64[  0   0    0  0.5  0  0.5  0   0   0   0   0   0   0   0   0   0   0   0   0; #  5
                         0   0    0   0  0.5 0.5  0   0   0   0   0   0   0   0   0   0   0   0   0; # 10
                        0.5  0    0   0   0  0.5  0   0   0   0   0   0   0   0   0   0   0   0   0; # 15
                         0  0.5   0   0   0   0  0.5  0   0   0   0   0   0   0   0   0   0   0   0; # 20
                         0   0   0.5  0   0   0   0  0.5  0   0   0   0   0   0   0   0   0   0   0; # 25
                         0   0    0  0.5  0   0   0   0  0.5  0   0   0   0   0   0   0   0   0   0; # 30
                         0   0    0   0  0.5  0   0   0   0  0.5  0   0   0   0   0   0   0   0   0; # 35
                         0   0    0   0   0  0.5  0   0   0   0  0.5  0   0   0   0   0   0   0   0; # 40
                         0   0    0   0   0   0  0.5  0   0   0   0  0.5  0   0   0   0   0   0   0; # 45
                         0   0    0   0   0   0   0  0.5  0   0   0   0  0.5  0   0   0   0   0   0; # 50
                         0   0    0   0   0   0   0   0  0.5  0   0   0   0  0.5  0   0   0   0   0; # 55
                         0   0    0   0   0   0   0   0   0  0.5  0   0   0   0  0.5  0   0   0   0; # 60
                         0   0    0   0   0   0   0   0   0   0  0.5  0   0   0   0  0.5  0   0   0; # 65
                         0   0    0   0   0   0   0   0   0   0   0  0.5  0   0   0   0  0.5  0   0; # 70
                         0   0    0   0   0   0   0   0   0   0   0   0  0.5  0   0   0   0  0.5  0; # 75
                         0   0    0   0   0   0   0   0   0   0   0   0   0  0.5  0   0   0   0  0.5; # 80
                         0   0    0   0   0   0   0   0   0   0   0   0   0  0.5 0.5  0   0   0   0; # 85
                         0   0    0   0   0   0   0   0   0   0   0   0   0  0.5  0  0.5  0   0   0; # 90
                         0   0    0   0   0   0   0   0   0   0   0   0   0  0.5  0   0  0.5  0   0] # 95

    c = ScoreTransitionMatrix(c_matrix0, RankInfo("C", 1, [5*i for i=1:19], 15, -10))

    c_plus_matrix0 = zeros(49, 49)

    for i=1:49
        if i+6 <= 49
            c_plus_matrix0[i,i+6] = 0.5 # +12
        else
            c_plus_matrix0[i,35] = 0.5 # 70
        end

        if i-5 >= 1
            c_plus_matrix0[i,i-5] = 0.5 # -10
        else
            c_plus_matrix0[i,15] = 0.5 # 30
        end
    end

    c_plus_and_b_minus = ScoreTransitionMatrix(c_plus_matrix0, RankInfo("C+/B-", 2, [2*i for i=1:49], 12, -10))

                      # 10  20  30  40  50  60  70  80  90
    b_matrix0= Float64[  0  0.5 0.5  0   0   0   0   0   0; # 10
                        0.5  0  0.5  0   0   0   0   0   0; # 20
                         0  0.5  0  0.5  0   0   0   0   0; # 30
                         0   0  0.5  0  0.5  0   0   0   0; # 40
                         0   0   0  0.5  0  0.5  0   0   0; # 50
                         0   0   0   0  0.5  0  0.5  0   0; # 60
                         0   0   0   0   0  0.5  0  0.5  0; # 70
                         0   0   0   0   0   0  0.5  0  0.5;# 80
                         0   0   0   0   0   0  0.5 0.5  0];# 90

    b_thru_a_plus = ScoreTransitionMatrix(b_matrix0, RankInfo("B/B+/A-/A/A+", 5, [10*i for i=1:9], 10, -10))

    s_matrix0 = zeros(99, 99)

    for i=1:99
        if i <= 39
            reward = 5
        else
            reward = 4
        end
        if i <= 79
            penalty = 5
        else
            penalty = 6
        end

        if i+reward <= 99
            s_matrix0[i,i+reward] = 0.5
        else
            s_matrix0[i,70] = 0.5
        end

        if i-penalty >= 1
            s_matrix0[i,i-penalty] = 0.5
        else
            s_matrix0[i,30] = 0.5
        end
    end

    s = ScoreTransitionMatrix(s_matrix0, RankInfo("S", 1, [i for i=1:99], 4, -5))

    s_plus_matrix0 = zeros(99, 99)

    for i=1:99
        if i <= 39
            reward = 4
        elseif i <= 79
            reward = 3
        else
            reward = 2
        end
        if i <= 39
            penalty = 4
        else
            penalty = 5
        end

        if i+reward <= 99
            s_plus_matrix0[i,i+reward] = 0.5
        else
            s_plus_matrix0[i,99] = 0.5
        end

        if i-penalty >= 1
            s_plus_matrix0[i,i-penalty] = 0.5
        else
            s_plus_matrix0[i,30] = 0.5
        end
    end

    s_plus = ScoreTransitionMatrix(s_plus_matrix0, RankInfo("S+", 1, [i for i=1:99], 2, -4))
#                          # 10  20  30  40  50  60  70  80  90 100
#   a_plus_matrix0= Float64[  0  0.5 0.5  0   0   0   0   0   0   0; # 10
#                            0.5  0  0.5  0   0   0   0   0   0   0; # 20
#                             0  0.5  0  0.5  0   0   0   0   0   0; # 30
#                             0   0  0.5  0  0.5  0   0   0   0   0; # 40
#                             0   0   0  0.5  0  0.5  0   0   0   0; # 50
#                             0   0   0   0  0.5  0  0.5  0   0   0; # 60
#                             0   0   0   0   0  0.5  0  0.5  0   0; # 70
#                             0   0   0   0   0   0  0.5  0  0.5  0;# 80
#                             0   0   0   0   0   0   0  0.5  0  0.5;# 90
#                             0   0   0   0   0   0   0   0  0.5 0.5];# 100

#   a_plus = ScoreTransitionMatrix(a_plus_matrix0, RankInfo("A+", 1, [10*i for i in 1:10], 10, -10))

    return ScoreTransitionMatrix[ c_minus, c, c_plus_and_b_minus, b_thru_a_plus, s, s_plus ]
end

function solve_and_report(score_transitions::Array{ScoreTransitionMatrix})
    score_steady_states = [score_steady_state(score_transitions[i]) for i in 1:length(score_transitions)]
    report_score_steady_states(score_steady_states)

    rank_matrix = rank_matrix_from_steady_states(score_steady_states)

    println("Rank transition matrix: ", rank_matrix)

    rank_steady_state = steady_state(rank_matrix)

    report_rank_steady_state(rank_steady_state)
end

println("\n *** Current Splatoon algorithm (demote when score < 0) *** ")
solve_and_report(default_score_transitions())

println("\n *** Proposed fix (demote when score = 0) *** ")
solve_and_report(proposed_score_transitions())

max_stars_repo_name: evanmiller/Splatoonalysis

max_issues_repo_name: evanmiller/Splatoonalysis

max_forks_repo_name: evanmiller/Splatoonalysis


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libsodium_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libsodium")
JLLWrappers.@generate_main_file("libsodium", UUID("a9144af2-ca23-56d9-984f-0d03f7b5ccf8"))
end  # module libsodium_jll


text:
using Test
using Yao
using QuAlgorithmZoo
using Random

function run_test(nbit::Int, depth_gen::Int, depth_disc::Int; g_lr=0.1, d_lr=0.2, niter=1000)
    qg = toy_qugan(rand_state(nbit), depth_gen, depth_disc)
    for info in QuGANGo!(qg, g_lr, d_lr, niter) end
    qg
end

# to fix
@testset "quantum circuit gan - opdiff" begin
    Random.seed!(2)
    N = 3
    target = rand_state(N)
    qcg = toy_qugan(target, 2, 2)
    grad = gradient(qcg)
    @test isapprox(grad, num_gradient(qcg), atol=1e-4)
    qg = run_test(3, 4, 4, g_lr=0.2, d_lr=0.5, niter=300)
    @test qg |> loss < 0.1
    qg = run_test(3, 4, 4, g_lr=Adam(lr=0.005), d_lr=Adam(lr=0.5), niter=1000)
    @test qg |> loss < 0.1
end


text:
function parse_paths(wires)
    paths = []

    for w in wires
        coord = (0, 0)
        path = []
        for p in w
            dir = if p[1] == 'U'
                (0, 1)
            elseif p[1] == 'R'
                (1, 0)
            elseif p[1] == 'D'
                (0, -1)
            elseif p[1] == 'L'
                (-1, 0)
            else
                throw(ErrorException("Invalid direction: " * p[1]))
            end
            for i = 1:(parse(Int, p[2:end]))
                coord = coord .+ dir
                push!(path, coord)
            end
        end
        push!(paths, path)
    end

    return paths
end

function shortest_distance(wires)::Int
    paths = parse_paths(wires)

    return minimum(sum.((x -> abs.(x)).(intersect(paths[1], paths[2]))))
end

function shortest_delay(wires)::Int
    paths = parse_paths(wires)
    intersections = intersect(paths[1], paths[2])

    delays = []
    for i in intersections
        s1 = findfirst(x -> x==i, paths[1])
        s2 = findfirst(x -> x==i, paths[2])
        push!(delays, s1+s2)
    end

    return minimum(delays)
end

wires = open("wires.csv") do f
    split.(readlines(f), ",")
end
println(shortest_distance(wires))
println(shortest_delay(wires))


text:
module TestInputNoise
using Test
using LinearAlgebra
using GaussianDistributions
using CovarianceFunctions

@testset "input noise" begin
    # set up GP
    l = 1/4
    k = CovarianceFunctions.Lengthscale(CovarianceFunctions.EQ(), l)
    G = Gaussian(k)

    # set up data
    n = 32 # number of observations
    x = randn(n) # uncorrupted input

    σx = 5e-2 # standard deviation of input noise
    σy = 1e-2 # standard deviation of ouput noise
    εx = σx .* randn(n) # input noise
    εy = σy .* randn(n) # output noise

    f(x) = sin(2π*x)
    y = @. f(x + εx) + εy # output

    C = G | (x, y, σx^2) # canonical conditional distribution
    Ci = G | (x, σx^2, y, σy^2) # NIGP conditional distribution

    xs = 2maximum(abs, x) # outside of data range, where regular posterior mean is constant,
    @test var(Ci)(xs) ≈ var(C)(xs) # NIGP variance is unchanged to regular variance
    xs = 2rand(n) .- 1
    nincrease = sum(var(C, xs) .< var(Ci, xs))
    @test nincrease > n * 3/4 # NIGP tends to increase uncertainty in data-rich regions

    # using Plots
    # plot(Ci, label = "NIGP")
    # plot!(C, label = "GP")
    # scatter!(x + εx, y + εy, label = "data")
    # xs = range(minimum(x) - std(x), maximum(x) + std(x), length = 128)
    # plot!(xs, f, label = "f(x)")
    # gui()
end

end

# plotting uncertainty
# xs = range(-3, 3, length = 256)
# plot(xs, var(Ci), label = "NIGP")
# plot!(xs, var(C), label = "GP")
# scatter!(x + εx, zero(x), label = "data")

# d = var(Ci, xs) .- var(C, xs)
# println([t for t in zip(x, d)])
# i = findfirst(<(0), d)
# println(xs[i])
# println(var(Ci)(xs[i]))
# println(var(C)(xs[i]))
# scatter!([xs[i]], [var(Ci)(xs[i])])
# gui()