text:
module Judo

import Base: start, next, done, display, writemime
import JSON
import Mustache

include("walkdir.jl")
include("harvest.jl")
include("collate.jl")

# An iterator for the parse function: parsit(source) will iterate over the
# expressiosn in a string.
type ParseIt
    value::String
end

function parseit(value::String)
    ParseIt(value)
end

function start(it::ParseIt)
    1
end

function next(it::ParseIt, pos)
    (ex,newpos) = Base.parse(it.value, pos)
    ((it.value[pos:(newpos-1)], ex), newpos)
end

function done(it::ParseIt, pos)
    pos > length(it.value)
end

# A special dummy module in which a documents code is executed.
module WeaveSandbox
end

# A super-simple pandoc interface.
#
# Args:
#   input: Input string.
#   infmt: Input format.
#   outfmt: Output format.
#   args: Additional arguments appended to the pandoc command.
#
# Returns:
#   A string containing the output from pandoc.
#
function pandoc(input::String, infmt::Symbol, outfmt::Symbol, args::String...)
    cmd = ByteString["pandoc",
                     "--from=$(string(infmt))",
                     "--to=$(string(outfmt))"]
    for arg in args
        push!(cmd, arg)
    end
    pandoc_out, pandoc_in, proc = readandwrite(Cmd(cmd))
    write(pandoc_in, input)
    close(pandoc_in)
    readall(pandoc_out)
end

# Convert a string to a structure interpretable by pandoc json parser.
function string_to_pandoc_json(s::String)
    out = {}
    for mat in eachmatch(r"\S+", s)
        if !isempty(out)
            push!(out, {"t" => "Space", "c" => {}})
        end
        push!(out, {"t" => "Str", "c" => mat.match})
    end
    out
end

# A Display implementation that renders multimedia into the document currently
# being processed.
type WeaveDoc <: Display
    # JSON respresentation of the document being built up.
    blocks::Array

    # Queued output from display calls
    display_blocks::Array

    # Document name (used for naming any external files generated)
    name::String

    # Output formal
    outfmt::Symbol

    # Current figure number
    fignum::Int

    # Redirected stdout
    stdout_read
    stdout_write

    # Output directory
    outdir::String

    function WeaveDoc(name::String, outfmt::Symbol, stdout_read, stdout_write,
                      outdir::String)
        new({}, {}, name, outfmt, 1, stdout_read, stdout_write, outdir)
    end
end

# Display functions that render output and insert them into the document.

const supported_mime_types =
    { MIME"text/html",
      MIME"image/svg+xml",
      MIME"image/png",
      MIME"text/latex",
      MIME"text/vnd.graphviz",
      MIME"text/plain" }

function display(doc::WeaveDoc, data)
    for m in supported_mime_types
        if mimewritable(m(), data)
            display(doc, m(), data)
            break
        end
    end
end

function display(doc::WeaveDoc, m::MIME"text/plain", data)
    block =
        {"t" => "CodeBlock",
         "c" => {{"", {"output"}, {}}, stringmime(m, data)}}
    push!(doc.display_blocks, block)
end

function display(doc::WeaveDoc, m::MIME"text/latex", data)
    # latex to dvi
    input_path, input = mktemp()
    writemime(input, m, data)
    flush(input)
    seek(input, 0)
    latexout_dir = mktempdir()
    run(`latex -output-format=dvi -output-directory=$(latexout_dir) $(input_path)` |> SpawnNullStream())
    rm(input_path)

    # dvi to svg
    latexout_path = "$(latexout_dir)/$(basename(input_path)).dvi"
    output = readall(`dvisvgm --stdout --no-fonts $(latexout_path)` .> SpawnNullStream())
    run(`rm -rf $(latexout_dir)`)

    display(doc, MIME("image/svg+xml"), output)
end

function display(doc::WeaveDoc, m::MIME"image/svg+xml", data)
    filename = @sprintf("%s_figure_%d.svg", doc.name, doc.fignum)
    out = open(joinpath(doc.outdir, filename), "w")
    writemime(out, m, data)
    close(out)

    alttext = @sprintf("Figure %d", doc.fignum)
    figurl = filename
    caption = ""

    if doc.outfmt == :html || doc.outfmt == :html5
        block =
            {"t" => "RawBlock",
             "c" => {"html",
               """

                 <object alt="$(alttext)" data="$(figurl)" type="image/svg+xml">
                 </object>

                   $(caption)

               """}}
    else
        block =
            {"t" => "Para",
             "c" => {
                {"t" => "Image",
                 "c" => {{{"Str" => alttext}},
                 {figurl, ""}}}}}
    end

    doc.fignum += 1
    push!(doc.display_blocks, block)
end

function display(doc::WeaveDoc, m::MIME"image/png", data)
    filename = @sprintf("%s_figure_%d.png", doc.name, doc.fignum)
    out = open(joinpath(doc.outdir, filename), "w")
    writemime(out, m, data)
    close(out)

    alttext = @sprintf("Figure %d", doc.fignum)
    figurl = filename
    caption = ""

    block =
        {"t" => "Para",
         "c" => {
            {"t" => "Image",
             "c" => {{{"t" => "Str", "c" => alttext}},
                     {figurl, ""}}}}}

    doc.fignum += 1
    push!(doc.display_blocks, block)
end

function display(doc::WeaveDoc, m::MIME"text/vnd.graphviz", data)
    output, input, proc = readandwrite(`dot -Tsvg`)
    writemime(input, m, data)
    close(input)
    display(doc, MIME("image/svg+xml"), readall(output))
end

function display(doc::WeaveDoc, m::MIME"text/html", data)
    block = {"t" => "RawBlock",
             "c" =>  {"html",
                string("
\n", stringmime(m, data), "\n
")}}
    push!(doc.display_blocks, block)
end

# This is maybe an abuse. TODO: This is going to be a problem.
writemime(io, m::MIME"text/vnd.graphviz", data::String) = write(io, data)
writemime(io, m::MIME"image/vnd.graphviz", data::String) = write(io, data)
writemime(io, m::MIME"image/svg+xml", data::String) = write(io, data)
writemime(io, m::MIME"text/latex", data::String) = write(io, data)

# Turn YAML frontmatter parsed by pandoc into simplified julia structure.
#
# Args:
#   pandoc_metadata: YAML fontmatter parsed by pandoc then by JSON
#
# Returns:
#   A (String => String) dictionary with key value pairs.
#
function flatten_pandoc_metadata(pandoc_metadata::Dict)
    metadata = Dict{String, String}()
    for (key, val) in pandoc_metadata["unMeta"]
        if val["t"] == "MetaInlines"
            metadata[key] = flatten_pandoc_metainlines(val["c"])
        elseif val["t"] == "MetaString"
            metadata[key] = val["c"]
        end
    end
    return metadata
end

function flatten_pandoc_metainlines(pandoc_metainlines::Vector)
    buf = IOBuffer()
    for val in pandoc_metainlines
        if val["t"] == "Str"
            write(buf, val["c"])
        elseif val["t"] == "Space"
            write(buf, " ")
        end
    end
    return takebuf_string(buf)
end

# Strip parsed pandoc JSON to a plain text representation.
function plaintext(blocks::Array)
    buf = IOBuffer()
    atstart = true
    for block in blocks
        if block["t"] == "Space"
            write(buf, " ")
        elseif block["t"] == "Str"
            write(buf, block["c"])
        elseif block["t"] == "Cite"
            write(buf, plaintext(block["c"][2]))
        end
        atstart = false
    end
    return takebuf_string(buf)
end

# Transform a annotated markdown file into a variety of formats.
#
# This reads markdown input, optionally prefixed with a YAML metadata document,
# and transforms in into a document in the desired output format, while
# executing code blocks and inserting results into the document.
#
# Args:
#   input: An input source.
#   output: Where the resulting document should be written.
#   outfmt: One of the output formats supported by pandoc:
#       Eg. markdown, rst, html, json, latex.
#   name: A short name for the document used for naming output files.
#
# Returns:
#   A pair (metadata, sections) where
#   `metadata` is metadata parsed from the YAML front matter if any, and
#   `sections` is an array of section/subsection names in the order in which
#   they occur. Each section name is a pair (level, name) where `level` an
#   integer giving the section level (1 is a section, 2 a sub-section, etc).
#
function weave(input::IO, output::IO;
               outfmt=:html5, name="judo", template=nothing,
               toc::Bool=false, outdir::String=".", dryrun::Bool=false,
               keyvals::Dict=Dict(),
               pandocargs=nothing)
    input_text = readall(input)

    # first pandoc pass
    pandoc_metadata, document =
        JSON.parse(pandoc(input_text, :markdown, :json))
    metadata = flatten_pandoc_metadata(pandoc_metadata)
    pandoc_metadata = pandoc_metadata["unMeta"]
    prev_stdout = STDOUT
    stdout_read, stdout_write = redirect_stdout()
    doc = WeaveDoc(name, outfmt, stdout_read, stdout_write, outdir)
    pushdisplay(doc)

    sections = {}

    for block in document
        if block["t"] == "Header"
            level = block["c"][1]
            headername = plaintext(block["c"][3])
            push!(sections, (level, headername))
            if !dryrun
                push!(doc.blocks, process_block(block))
            end
        elseif dryrun
            continue
        elseif block["t"] == "CodeBlock"
            try
                process_code_block(doc, block)
            catch err
                error("Error processing codeblock in document $(name):\n$(string(err))")
            end
        else
            push!(doc.blocks, process_block(block))
        end
    end

    if dryrun
        return metadata, sections
    end

    popdisplay(doc)
    redirect_stdout(prev_stdout)
    #Base.reinit_stdio()

    # second pandoc pass
    buf = IOBuffer()

    # splice the document's creation date into pandoc's metadata
    pandoc_metadata["today"] =
        {"t" => "MetaInlines",
         "c" => {{"t" => "Str",
                  "c" => strip(readall(`date`))}}}

    JSON.print(buf, {{"unMeta" => pandoc_metadata}, doc.blocks})

    args = {}
    for (k, v) in keyvals
        push!(args, "--variable=$(k):$(v)")
    end

    if template != nothing
        push!(args, "--template=$(template)")
    end

    if toc
        push!(args, "--toc")
    end

    if pandocargs != nothing
        push!(args, pandocargs)
    end

    write(output, pandoc(takebuf_string(buf), :json, outfmt, args...))

    metadata, sections
end

# Evaluate the truthiness of a code block attribute.
function codeblock_keyval_bool(keyvals::Dict, key, default::Bool)
    if haskey(keyvals, key)
        val = lowercase(strip(keyvals[key]))
        @assert val == "true" || val == "false"
        val == "true"
    else
        default
    end
end

# Process a text paragraph.
#
# Args:
#   block: Block parsed from the json document representation.
#
# Modifies:
#    block, prior to inserting it into doc
#
# Returns:
#   block
#
#function process_block(block::Dict)
    #if block["t"] == "Code"
        #(id, classes, keyvals), text = block["c"]
        #if "julia" in classes
            #out = eval(WeaveSandbox, parse(text))
            #block["c"] = string(out)
        #end
    #else
        #block["c"] = process_block(block["c"])
    #end
    #return block
#end

#function process_block(block::Array)
    #for i in 1:length(block)
        #block[i] = process_block(block[i])
    #end
    #block
#end

function process_block(block::Any)
    block
end

# Code block classes supported.
const code_block_classes = Set(
    ["julia",
     "graphviz",
     "latex",
     "svg"])

# Process code blocks.
#
# This consists of (optioclasses =ally) executing code and inserting the results into
# the document.
#
# Args:
#   doc: Document being generated.
#   block:: Code block.
#
function process_code_block(doc::WeaveDoc, block::Dict)
    (id, classes, keyvals_array), text = block["c"]
    keyvals = [k => v for (k, v) in keyvals_array]

    # Options are:
    #
    #  hide: Don't show the code block. (defaut: false)
    #  execute: Do excute the code block. (default: true)
    #  display: Do display output. (default: true)
    #  results:
    #    none (default)
    #    block
    #    expression

    keyvals["hide"]    = codeblock_keyval_bool(keyvals, "hide",    false)
    keyvals["execute"] = codeblock_keyval_bool(keyvals, "execute", true)
    keyvals["display"] = codeblock_keyval_bool(keyvals, "display", true)
    keyvals["results"] = get(keyvals, "results", "block")

    if isempty(classes) || classes[1] == "julia" ||
        !in(classes[1], code_block_classes)

        if keyvals["results"] == "none" || keyvals["results"] == "block"
            result = nothing
            if keyvals["execute"]
                for (cmd, ex) in parseit(strip(text))
                    result = safeeval(ex)
                end
            end

            if keyvals["results"] == "block" && result != nothing
                display(result)
            end
            output_text = text
        elseif keyvals["results"] == "expression"
            if keyvals["execute"]
                output_text = IOBuffer()
                for (cmd, ex) in parseit(strip(text))
                    println(output_text, cmd)
                    println(output_text, "## ", string(safeeval(ex)), "\n")
                end
            else
                output_text = text
            end
        end

        if !keyvals["hide"]
            push!(doc.blocks,
            {"t" => "CodeBlock", "c" => {{id, {"julia", classes...}, keyvals_array}, output_text}})
        end

        if keyvals["display"]
            flush(doc.stdout_read)
            flush(doc.stdout_write)
            if nb_available(doc.stdout_read) > 0
                stdout_output = readavailable(doc.stdout_read)
                if !isempty(stdout_output)
                    display("text/plain", stdout_output)
                end
            end
            append!(doc.blocks, doc.display_blocks)
        end
        empty!(doc.display_blocks)
    else
        if !keyvals["hide"]
            push!(doc.blocks, block)
        end

        class = classes[1]
        if keyvals["display"]
            if class == "graphviz"
                display("text/vnd.graphviz", text)
            elseif class == "latex"
                display("text/latex", text)
            elseif class == "svg"
                display("image/svg+xml", text)
            end
            append!(doc.blocks, doc.display_blocks)
        end
    end
end

# Evaluate an expression and return its result and a string.
function safeeval(ex::Expr)
    eval(WeaveSandbox, ex)
end

end


text:
#=
Various routines to discretize AR(1) processes

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date : 2014-04-10 23:55:05

References
----------

http://quant-econ.net/jl/finite_markov.html
=#

std_norm_cdf{T <: Real}(x::T) = 0.5 * erfc(-x/sqrt(2))
std_norm_cdf{T <: Real}(x::Array{T}) = 0.5 .* erfc(-x./sqrt(2))

"""
Tauchen's (1996) method for approximating AR(1) process with finite markov chain

The process follows

    y_t = μ + ρ y_{t-1} + ε_t,

where ε_t ~ N (0, σ^2)

##### Arguments

- `N::Integer`: Number of points in markov process
- `ρ::Real` : Persistence parameter in AR(1) process
- `σ::Real` : Standard deviation of random component of AR(1) process
- `μ::Real(0.0)` : Mean of AR(1) process
- `n_std::Integer(3)` : The number of standard deviations to each side the process
should span

##### Returns

- `y::Vector{Real}` : Nodes in the state space
- `Π::Matrix{Real}` Matrix transition probabilities for Markov Process

"""
function tauchen(N::Integer, ρ::Real, σ::Real, μ::Real=0.0, n_std::Integer=3)
    # Get discretized space
    a_bar = n_std * sqrt(σ^2 / (1 - ρ^2))
    y = linspace(-a_bar, a_bar, N)
    d = y[2] - y[1]

    # Get transition probabilities
    Π = zeros(N, N)
    for row = 1:N
        # Do end points first
        Π[row, 1] = std_norm_cdf((y[1] - ρ*y[row] + d/2) / σ)
        Π[row, N] = 1 - std_norm_cdf((y[N] - ρ*y[row] - d/2) / σ)

        # fill in the middle columns
        for col = 2:N-1
            Π[row, col] = (std_norm_cdf((y[col] - ρ*y[row] + d/2) / σ) -
                           std_norm_cdf((y[col] - ρ*y[row] - d/2) / σ))
        end
    end

    # NOTE: I need to shift this vector after finding probabilities
    #       because when finding the probabilities I use a function
    #       std_norm_cdf that assumes its input argument is distributed
    #       N(0, 1). After adding the mean E[y] is no longer 0, so
    #       I would be passing elements with the wrong distribution.
    #
    #       It is ok to do after the fact because adding this constant to each
    #       term effectively shifts the entire distribution. Because the
    #       normal distribution is symmetric and we just care about relative
    #       distances between points, the probabilities will be the same.
    #
    #       I could have shifted it before, but then I would need to evaluate
    #       the cdf with a function that allows the distribution of input
    #       arguments to be [μ/(1 - ρ), 1] instead of [0, 1]

    y .+= μ / (1 - ρ) # center process around its mean (wbar / (1 - rho))

    return y, Π
end

"""
Rouwenhorst's method to approximate AR(1) processes.

The process follows

    y_t = μ + ρ y_{t-1} + ε_t,

where ε_t ~ N (0, σ^2)

##### Arguments
- `N::Integer` : Number of points in markov process
- `ρ::Real` : Persistence parameter in AR(1) process
- `σ::Real` : Standard deviation of random component of AR(1) process
- `μ::Real(0.0)` :  Mean of AR(1) process

##### Returns

- `y::Vector{Real}` : Nodes in the state space
- `Θ::Matrix{Real}` Matrix transition probabilities for Markov Process

"""
function rouwenhorst(N::Integer, ρ::Real, σ::Real, μ::Real=0.0)
    σ_y = σ / sqrt(1-ρ^2)
    p  = (1+ρ)/2
    Θ = [p 1-p; 1-p p]
    ψ = sqrt(N-1) * σ_y
    m = μ / (1 - ρ)

    return rouwenhorst(p, p, m, ψ, N)
end

function rouwenhorst(p::Real, q::Real, m::Real, Δ::Real, n::Integer)
    if n == 2
        return Real[m-Δ, m+Δ], [p 1-p; 1-q q]
    else
        _, θ_nm1 = rouwenhorst(p, q, m, Δ, n-1)
        θN = p    *[θ_nm1 zeros(n-1, 1); zeros(1, n)] +
             (1-p)*[zeros(n-1, 1) θ_nm1; zeros(1, n)] +
             q    *[zeros(1, n); zeros(n-1, 1) θ_nm1] +
             (1-q)*[zeros(1, n); θ_nm1 zeros(n-1, 1)]

        θN[2:end-1, :] ./= 2

        return linspace(m-Δ, m+Δ, n), θN
    end
end

ext: jl


text:
function testNHInertiaGravityWave()
close all

% Grid
nx=60;
ny=2;
nz=10;
Param.Lx=300*1.e3;
Param.Ly=2*1.e3;
x0=0;
y0=0;
Param.H=10*1.e3;

Boundary.WE='Period';
Boundary.BT='Period';
Param.hS='';
Param.Grid=CartGrid(nx,ny,Param.Lx,Param.Ly,x0,y0,OrientFaceCart,Boundary,Param);

Param.Grid.nz=nz;
Param.Grid.zP=zeros(nz,1);
Param.Grid.z=zeros(nz+1,1);
Param.Grid.dz=Param.H/nz;
Param.Grid.zP(1)=Param.Grid.dz/2;
for i=2:nz
  Param.Grid.zP(i)=Param.Grid.zP(i-1)+Param.Grid.dz;
end
for i=2:nz+1
  Param.Grid.z(i)=Param.Grid.z(i-1)+Param.Grid.dz;
end

% Model
Param.ProfRho='InertiaGravityWave';
Param.ProfTheta='InertiaGravityWave';
Param.ProfVel='Const';
Param.Damping=false;
Param.Coriolis=false;
Param.Buoyancy=true;
Param.HyperD=0; %1.e15*(nx/30)^3.2;
Param.xC=Param.Lx/3;
Param.a=5*1.e3;
Param.Th0=300;
Param.DeltaTh=1.e-2;
Param.uMax=0;
Param.vMax=0;
Param.NBr=1.e-2;
Param.Equation='Compressible';
Param.TopoS='';
Param.NumV=5;
Param.RhoPos=1;
Param.uPos=2;
Param.vPos=3;
Param.wPos=4;
Param.ThPos=5;
Param.Thermo='Energy';

% Physical parameters
Param.Grav=9.81;
Param.Coriolis=false;
Param.cS=360;
Param.Grav=9.81d0;
Param.Cpd=1004.0d0;
Param.Cvd=717.0d0;
Param.Rd=Param.Cpd-Param.Cvd;
Param.p0=1.0d5;
Param.Cpv=1885.0d0;
Param.Gamma=Param.Cpd/Param.Cvd;
Param.kappa=Param.Rd/Param.Cpd;

% Output
Param.Flat=false;
Param.level=1;
Param.fig=1;
Param.SliceXZ.Type='XZ';
Param.SliceXZ.y=y0+Param.Ly/ny/2;
Param.SliceXZ.iy=1;

% Discretization

OrdPoly=3;
[CG,Param]=Discretization(OrdPoly,Param);

% Initial conditions

U=zeros(CG.NumG,nz,Param.NumV);
U(:,:,Param.RhoPos)=Project(@fRho,CG,Param);
[U(:,:,Param.uPos),U(:,:,Param.vPos)]=ProjectVec(@fVel,CG,Param);
U(:,:,Param.ThPos)=Project(@fTheta,CG,Param).*U(:,:,Param.RhoPos);
if strcmp(Param.Thermo,'Energy')
  U(:,:,Param.ThPos)=PotToEnergy(U,CG,Param);
end
ThB=Project(@fThetaBGrd,CG,Param);
Param.fig=PlotCG(U(:,:,Param.ThPos)./U(:,:,Param.RhoPos)-ThB...
  ,CG,@TransCart,@Topo,Param,Param.fig,Param.SliceXZ);

% Integration
CFL=0.125;
dtau=.40;
time=0;

%IntMethod='Rosenbrock'; %'RungeKutta';
IntMethod='RungeKutta';
nIter=20000;
%v = VideoWriter ('Galewsky.avi');
%open (v);
Param.RK=RungeKuttaMethod('RK4');
Param.ROS=RosenbrockMethod('ROSRK3');
switch IntMethod
  case 'Rosenbrock'
    for i=1:nIter
      U=Rosenbrock(U,dtau,@FcnNHCurlVec,@Jac,CG,Param);
      time=time+dtau;
      if mod(i,2000)==0
        Param.fig=PlotCG(U(:,:,Param.ThPos)./U(:,:,Param.RhoPos)-ThB...
          ,CG,@TransCart,@Topo,Param,Param.fig,Param.SliceXZ);
%         Param.fig=PlotCG(U(:,:,Param.uPos),CG,@TransCart,@Topo,Param,Param.fig,Param.SliceXZ);
%         Param.fig=PlotCG(U(:,:,Param.vPos),CG,@TransCart,@Topo,Param,Param.fig,Param.SliceXZ);
%         W=zeros(size(U,1),nz+1);
%         W(:,2:nz+1)=U(:,:,Param.wPos);
%         Param.fig=PlotCG(0.5*(W(:,1:nz)+W(:,2:nz+1)),CG,@TransCart...
%           ,@Topo,Param,Param.fig,Param.SliceXZ);
      end
    end
  case 'RungeKutta'
    for i=1:nIter
      U=RungeKuttaExplicit(U,dtau,@FcnNHCurlVec,CG,Param);

      time=time+dtau;
      if mod(i,1000)==0
        if strcmp(Param.Thermo,'Energy')
          Th=EnergyToPot(U,CG,Param);
        else
          Th=U(:,:,Param.ThPos);
        end
        Param.fig=PlotCG(Th./U(:,:,Param.RhoPos)-ThB...
          ,CG,@TransCart,@Topo,Param,Param.fig,Param.SliceXZ);
        Param.fig=PlotCG(U(:,:,Param.uPos),CG,@TransCart,@Topo,Param,Param.fig,Param.SliceXZ);
        W=zeros(size(U,1),nz+1);
        W(:,2:nz+1)=U(:,:,Param.wPos);
        Param.fig=PlotCG(0.5*(W(:,1:nz)+W(:,2:nz+1)),CG,@TransCart...
          ,@Topo,Param,Param.fig,Param.SliceXZ);
%         frame=getframe(gcf);
%         writeVideo (v, frame);
      end
    end
end
%fig=PlotCG(U(:,1),CG,@JacobiSphere2,Param,fig);
%close(v)
end


text:
################################################################################
# Toy Example
################################################################################
using Plots
T = Float64

# Define the dynamics of the system
p = 2 # Number of players
model = DoubleIntegratorGame(p=p) # game with 3 players with unicycle dynamics
# model = UnicycleGame(p=p) # game with 3 players with unicycle dynamics
# model = BicycleGame(p=p) # game with 3 players with unicycle dynamics
n = model.n
m = model.m

# Define the horizon of the problem
N = 20 # N time steps
dt = 0.1 # each step lasts 0.1 second
probsize = ProblemSize(N,model) # Structure holding the relevant sizes of the problem

# Define the objective of each player
# We use a LQR cost
Q = [Diagonal(10*ones(SVector{model.ni[i],T})) for i=1:p] # Quadratic state cost
# Q = [Diagonal(10*ones(SVector{model.ni[1],T})), Diagonal(1*ones(SVector{model.ni[2],T}))] # Quadratic state cost
# Q = [Diagonal(1*ones(SVector{model.ni[1],T})), Diagonal(10*ones(SVector{model.ni[2],T}))] # Quadratic state cost
R = [Diagonal(0.1*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
# R = [Diagonal(0.1*ones(SVector{model.mi[1],T})), Diagonal(0.01*ones(SVector{model.mi[2],T}))] # Quadratic control cost
# R = [Diagonal(0.01*ones(SVector{model.mi[1],T})), Diagonal(0.1*ones(SVector{model.mi[2],T}))] # Quadratic control cost
# Desrired state
xf = [SVector{model.ni[1],T}([2, 0.0,0,0]),
      SVector{model.ni[2],T}([2, 0.0,0,0]),
      # SVector{model.ni[3],T}([3,-0.4,0,0]),
      # SVector{model.ni[4],T}([3,+0.8,0,0]),
      ]
# Desired control
uf = [zeros(SVector{model.mi[i],T}) for i=1:p]

# # Objectives of the game
game_obj = GameObjective(Q,R,xf,uf,N,model)
radius = 0.5*ones(p)
# μ = 4.0*ones(p)
# μ = 4.0*[1.0, 0.1]
# μ = 4.0*[0.1, 1.0]
# add_collision_cost!(game_obj, radius, μ)

# Define the constraints that each player must respect
game_con = GameConstraintValues(probsize)
# Add collision avoidance
radius = 0.19
add_collision_avoidance!(game_con, radius)

# # Add control bounds
# u_max =  5*ones(SVector{m,T})
# u_min = -5*ones(SVector{m,T})
# add_control_bound!(game_con, u_max, u_min)

# # Add wall constraint
# walls = [Wall([-1.0,-1.75], [1.0,0.25], [1.,-1.]/sqrt(2))]
# add_wall_constraint!(game_con, walls)

# # Add circle constraint
# xc = [1., 2., 3.]
# yc = [1., 2., 3.]
# radius = [0.1, 0.2, 0.3]
# add_circle_constraint!(game_con, xc, yc, radius)

# Define the initial state of the system
x0 = SVector{model.n,T}([
    0.0, 0.0,# 0.5, #0.0,
   -0.2, 0.2,# 0.7, #0.6,
    0.0, 0.0,# 0.0, #0.0,
    0.0, 0.0,# 0.0, #0.0,
    ])

# Define the Options of the solver
opts = Options()
opts.ls_iter = 15
opts.outer_iter = 20
opts.inner_iter = 20
opts.ρ_0 = 1e0
opts.reg_0 = 1e-7
opts.α_dual = 1.0
opts.λ_max = 1.0*1e7
opts.ϵ_dyn = 1e-6
opts.ϵ_sta = 1e-6
opts.ϵ_con = 1e-6
opts.ϵ_opt = 1e-6
opts.regularize = true
# opts.αx_dual = ones(p)
opts.αx_dual = [sqrt(10), 1/sqrt(10)]
s = 0.1
opts.αx_dual = [sqrt(s), 1/sqrt(s)]
# opts.αx_dual = [1/sqrt(10), sqrt(10)]
# Define the game problem
prob = GameProblem(N,dt,x0,model,opts,game_obj,game_con)

prob.game_con.state_conval[1][1]
# Solve the problem
@time newton_solve!(prob)
# @profiler newton_solve!(prob)

plot!(prob.model, prob.pdtraj.pr)
plot!(prob.stats)

residual!(prob, prob.pdtraj)
# regularize_residual!(prob.core, prob.opts, prob.pdtraj, prob.pdtraj)
residual_jacobian!(prob, prob.pdtraj)
# regularize_residual_jacobian!(prob)
prob.core.res

Δtraj = - \(lu(prob.core.jac), prob.core.res)
set_traj!(prob.core, prob.Δpdtraj, Δtraj)
prob.Δpdtraj
Δ_step(prob.Δpdtraj, 1.0)

vecs = eigvecs(Matrix(prob.core.jac))
vals = eigvals(Matrix(prob.core.jac))
vals_norm = norm.(vals)
min_val = minimum(vals_norm)
i_min = argmin(vals_norm)
min_vec = vecs[:,i_min]

mean(norm.((prob.core.jac*min_vec - min_val*min_vec)))

function get_min_eig(prob::GameProblem)
    residual!(prob, prob.pdtraj)
    res = eigs(prob.core.jac, which=:SM, nev=1)
    min_val = res[1]
    min_vec = res[2]

    return res[1:2]
end

function evolve(prob::GameProblem, dt::T) where {T}
    eig_pdtraj = PrimalDualTraj(prob.probsize, dt)
    for j = 1:10
        min_val, min_vec = get_min_eig(prob)
        # @show min_val
        residual!(prob)
        @show mean(abs.(prob.core.res))
        set_traj!(prob.core, eig_pdtraj, real.(min_vec)[:,1])
        scale = mean(norm.(real.(min_vec)))
        α = 1e-2/scale
        update_traj!(prob.pdtraj, prob.pdtraj, α, eig_pdtraj)
        prob.opts.shift = 0
        newton_solve!(prob)
        plot!(prob.model, prob.pdtraj.pr)
    end
    return nothing
end

prob_copy = deepcopy(prob)
evolve(prob_copy, dt)

prob.game_con.state_conval[1]
prob.game_con.state_conval[1][1].λ
prob.game_con.state_conval[1][1].vals
prob.game_con.state_conval[1][1].jac[1]
prob.game_con.state_conval[1][1].con
prob.game_con.state_conval[2][1].λ
prob.game_con.state_conval[2][1].vals
prob.game_con.state_conval[2][1].con

eigs(prob.core.jac, which=:LM)
res = eigs(prob.core.jac, which=:SM, nev=1)
res[1]
res[2]
min_vec + res[2]

size(prob.core.jac)
mat = Matrix(log.(10, abs.(prob.core.jac)))
heatmap(mat, color = :greys) #which will also give you a color legend
heatmap(mat[190:290,1:100], color = :greys) #which will also give you a color legend

# kick out of inner loop if not progress made
# add probsize to game_con struct
# add constructor for multiple radius of collision

function traj_dist(prob::GameProblem ;M::Int=10)
    traj = []
    for j = 1:M
        prob.opts.seed = j
        newton_solve!(prob)
        deepcopy(prob.pdtraj.pr)
        push!(traj, deepcopy(prob.pdtraj.pr))
    end
    plotset(prob.model, traj)
    return traj
end

function plotset(model::AbstractGameModel, traj)
    plt = plot(aspect_ratio=:equal)
    M = length(traj)
    N = length(traj[1])
    col = [:blue, :orange, :green]
    for j = 1:M
        for i = 1:model.p
            xi = [state(traj[j][k])[model.pz[i][1]] for k=1:N]
            yi = [state(traj[j][k])[model.pz[i][2]] for k=1:N]
            plot!(xi, yi, color=col[i], label=false)
            scatter!(xi, yi, color=col[i], label=false)
        end
    end
    display(plt)
    return nothing
end

function solver_progress(prob::GameProblem; M::Int=10, II::Int=20, OI::Int=20)
    anim = @animate for ii = 1:Int(ceil(II/2)):II
        prob.opts.outer_iter=OI
        prob.opts.inner_iter=ii
        traj_dist(prob, M=M)
    end
    gif(anim, "Documents/solver_progress_con.gif", fps = 2)
    return nothing
end

prob.opts.amplitude_init = 1e1
prob.opts.f_init = x -> (rand(x) .- 0.5)
prob.opts.outer_iter = 0
prob.opts.inner_iter = 0
traj_dist(prob, M=100)

prob_copy = deepcopy(prob)
traj_dist(prob_copy, M=100)
solver_progress(prob_copy, M=50, II=20, OI=20)
newton_solve!(prob_copy)
prob_copy.stats

# function balance_dual_update!(game_con::GameConstraintValues, game_obj::GameObjective,
#     pdtraj::PrimalDualTraj)
#     probsize = game_obj.probsize
#     N = probsize.N
#     p = probsize.p
#     pu = probsize.pu
#
#     cost_gradient!(game_obj, pdtraj)
#     @show 111
#     for i = 1:p
#         for k = 2:N
#             n_obj = length(game_obj.E[i])
#             for j = 1:n_obj
#                 ∇x_Jik = game_obj.E[i][j].cost[k].q
#                 n_con = length(game_con.state_conlist[i].constraints)
#                 for l = 1:n_con
#                     conval = game_con.state_conval[i][l]
#                     TrajectoryOptimization.evaluate!(conval, pdtraj.pr)
#                     TrajectoryOptimization.jacobian!(conval, pdtraj.pr)
#                     TrajectoryOptimization.cost_expansion!(conval)
#                     ∇x_Cik = conval.jac[k-1]# need to use enumarate to coincide k+1 and k
#                     @show i
#                     @show -(∇x_Cik*∇x_Cik')*∇x_Cik*∇x_Jik
#                 end
#             end
#         end
#     end
#
#     return nothing
# end

colmult = CollisionMultiplier(probsize)
get_collision_multiplier!(colmult, prob.game_con)
get_balance!(colmult, prob.game_con, prob.pdtraj)
balance_dual!(colmult, prob.game_con)
unbalance_dual!(colmult, prob.game_con)

prob.game_con.state_conval[1]

balance_dual_update!(prob.game_con, prob.game_obj, prob.pdtraj)

prob.game_obj

prob.game_con.state_conlist[1]
prob.game_con.state_conval[1][1]

function residual2!(prob::GameProblem{KN,n,m,T,SVd,SVx}, pdtraj::PrimalDualTraj{KN,n,m,T,SVd}) where {KN,n,m,T,SVd,SVx}
	N = prob.probsize.N
	p = prob.probsize.p
    pu = prob.probsize.pu
	core = prob.core
	model = prob.model
	game_obj = prob.game_obj

	# Initialization
	prob.core.res .= 0.0
    stamp = VStamp()
	∇dyn = zeros(MMatrix{n,(n+m),T,n*(n+m)})

	# Cost
	cost_gradient!(game_obj, pdtraj)
	for i = 1:p
		# State cost
		for k = 1:N
			stampify!(stamp, :opt, i, :x, 1, k)
			n_obj = length(game_obj.E[i])
			for j = 1:n_obj
				valid(stamp, N, p) ? add2sub(core.res_sub[stamp], game_obj.E[i][j].cost[k].q) : nothing
			end
		end
		# Control Cost
		for k = 1:N-1
			stampify!(stamp, :opt, i, :u, i, k)
			n_obj = length(game_obj.E[i])
			for j = 1:n_obj
				valid(stamp, N, p) ? add2sub(core.res_sub[stamp], game_obj.E[i][j].cost[k].r[pu[i]]) : nothing
			end
		end
	end
	# # Dynamics penalty
	# for k = 1:N-1
	# 	∇dynamics!(∇dyn, model, pdtraj, k)
	# 	for i = 1:p
	# 		λik = pdtraj.du[i][k]
	# 		stampify!(stamp, :opt, i, :x, 1, k)
	# 		valid(stamp, N, p) ? add2sub(core.res_sub[stamp], ∇dyn[:,core.dyn[:x][1]]'*λik) : nothing
	# 		stampify!(stamp, :opt, i, :u, i, k)
	# 		valid(stamp, N, p) ? add2sub(core.res_sub[stamp], ∇dyn[:,core.dyn[:u][i]]'*λik) : nothing
	# 		stampify!(stamp, :opt, i, :x, 1, k+1)
	# 		valid(stamp, N, p) ? add2sub(core.res_sub[stamp], -λik) : nothing
	# 	end
	# end

	# # Constraints
	# constraint_residual!(prob, pdtraj)

	# # Dynamics
    # for k = 1:N-1
    #     stampify!(stamp, :dyn, 1, :x, 1, k)
    #     valid(stamp, N, p) ? add2sub(core.res_sub[stamp], dynamics_residual(model, pdtraj, k)) : nothing # shouldn't allocate
    # end
    return nothing
end

prob.opts.γ


text:
m = Model(GLPK.Optimizer)

@variable(m, x ≥ 0)

lista = []

push!(lista,@expression(m, x-1))
push!(lista,@expression(m, x-2))
push!(lista,@expression(m, x-3))


text:
include("../expmap.jl")
include("../rotqua.jl")

R = eye(3);
q = rotqua(R)
println("-------------------------------------------------------------")
println("for R = ", R)
println("q = ", q)

R = expmap([1.0 1.0 1.0])
q = rotqua(R)
println("-------------------------------------------------------------")
println("for R = ", R)
println("q = ", q)

max_stars_repo_name: tkoziara/edyna

max_issues_repo_name: tkoziara/edyna

max_forks_repo_name: tkoziara/edyna


text:
using SciMLBenchmarks
SciMLBenchmarks.weave_file(".","test.jmd")


text:
import Base: rand
export RPT∇, STL∇

# ReParametrisation Trick (RPT) gradient estimator. Assumes that `rand(q(ϕ))` is
# reparametrisable.
function RPT∇(rng::AbstractRNG, q, ϕ, log_π̃_over_q)
    elbo, back = Zygote.pullback(ϕ->log_π̃_over_q(ϕ, rand(rng, q(ϕ))), ϕ)
    ∂ϕ = first(back(1.0))
    return elbo, ∂ϕ
end

# Assumes single-sample. Multi-sample objectives should be viewed via AuxVI. No global
# parameters to learn for now. Include them later when this actually works.
function STL∇(rng::AbstractRNG, q, ϕ, log_π̃_over_q)

    # Compute approximation to elbo.
    z, back_z = Zygote.pullback(ϕ->rand(rng, q(ϕ)), ϕ)
    elbo, back_elbo = Zygote.pullback(log_π̃_over_q, ϕ, z)

    # Compute reverse-pass, dropping gradients appropriately.
    ∂ϕ_score, ∂z = back_elbo(1.0)
    ∂ϕ = first(back_z(∂z))

    return elbo, ∂ϕ, ∂ϕ_score
end


text:
__precompile__()

module ErgodicControl

import Base.normalize!
using SpecialFunctions

# export functions I've made
export ErgodicManager, decompose!, reconstruct, decompose
export ergodic_score, control_score, total_score
export control_effort
export centroid, covariance
export TrajectoryManager, dynamics!, sample_trajectory
export kmeans_trajectory, random_trajectory
export assign_step, mean_step

export
    Initializer,
    initialize,
    CornerConstantInitializer,
    GreedyInitializer,
    PointInitializer,
    DirectionInitializer

export
    x_min, y_min, z_min,
    x_max, y_max, z_max,
    x_size, y_size, z_size,
    x_cells, y_cells, z_cells

# to make some things easier
const T2F = NTuple{2, Float64}    # x, y
const MF = Matrix{Float64}
const VMF = Vector{MF}
const VF = Vector{Float64}
const VVF = Vector{VF}
const VVVF = Vector{VVF}

# math-type stuff I might need
include("math.jl")
include("tsp.jl")
include("lq.jl")
include("lqr.jl")

# ergodic manager
include("ergodic_manager/domain.jl")
include("ergodic_manager/gaussian.jl")
include("ergodic_manager/circle.jl")
include("ergodic_manager/ergodic_manager.jl")
include("ergodic_manager/r2.jl")
include("ergodic_manager/r3.jl")
include("ergodic_manager/r2t.jl")
include("ergodic_manager/se2.jl")
include("ergodic_manager/examples.jl")

# trajectory manager stuff
include("trajectory_manager/trajectory_manager.jl")
include("trajectory_manager/dynamics.jl")
include("trajectory_manager/initializers/initializer.jl")
include("trajectory_manager/descender.jl")

# trajectory generation
include("trajectory_generation/trajectory_generation.jl")

end # module


text:
struct Cglobal{T}
	handle::Ptr{T}

	Cglobal{T}(sym::Symbol, libs::Clibrary...) where {T} = new{T}(reinterpret(Ptr{T}, _dlsym(sym, libs...)))
end
Cglobal{T}(lib::Clibrary, sym::Symbol) where {T} = Cglobal{T}(sym, lib)

Base.getindex(g::Cglobal{T}) where {T} = unsafe_load(g.handle)
Base.setindex!(g::Cglobal{T}, val::T) where {T} = unsafe_store!(g.handle, val)


text:
n = 10 # rand(1:10)

@test matrixdepot("poisson", Float64, n) == matrixdepot("poisson", n)

A = Matrix(matrixdepot("poisson", n))

n_mesh, n = n, n^2

B = zeros(Float64, n, n)

i = 0
for i_mesh in 1 : n_mesh
    for j_mesh in 1 : n_mesh
        global i = i + 1

        if i_mesh > 1
            j = i - n_mesh
            B[i,j] = -1
        end

        if j_mesh > 1
            j = i - 1
            B[i,j] = -1
        end

        j = i
        B[i,j] = 4

        if j_mesh < n_mesh
            j = i + 1
            B[i,j] = -1
        end

        if i_mesh < n_mesh
            j = i + n_mesh
            B[i,j] = -1
        end
    end
end

@test A ≈ B

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


text:
"""
    struct AbstractGegenbauerBasis{P} <: AbstractMultipleOrthogonalBasis{P}
        polynomials::Vector{P}
    end

Orthogonal polynomial with respect to the univariate weight function ``w(x) = (1 - x^2)^{\\alpha - 1/2}`` over the interval ``[-1, 1]``.
"""
abstract type AbstractGegenbauerBasis{P} <: AbstractMultipleOrthogonalBasis{P} end

even_odd_separated(::Type{<:AbstractGegenbauerBasis}) = true
reccurence_second_coef(::Type{<:AbstractGegenbauerBasis}, degree) = 0

"""
    struct LegendreBasis{P} <: AbstractGegenbauerBasis{P}
        polynomials::Vector{P}
    end

Orthogonal polynomial with respect to the univariate weight function ``w(x) = 1`` over the interval ``[-1, 1]``.
"""
struct LegendreBasis{P} <: AbstractGegenbauerBasis{P}
    polynomials::Vector{P}
end

polynomial_type(::Type{<:LegendreBasis}, V::Type) = MP.polynomialtype(V, Float64)

reccurence_first_coef(::Type{<:LegendreBasis}, degree) = (2degree - 1)
reccurence_third_coef(::Type{<:LegendreBasis}, degree) = -(degree - 1)
reccurence_deno_coef(::Type{<:LegendreBasis}, degree) = degree

degree_one_univariate_polynomial(::Type{<:LegendreBasis}, variable::MP.AbstractVariable) = MA.@rewrite(variable + 0)

max_stars_repo_name: JuliaAlgebra/PolynomialBases

max_issues_repo_name: JuliaAlgebra/PolynomialBases

max_forks_repo_name: JuliaAlgebra/PolynomialBases


text:
struct Figure11
    θ_est::Vector{Vector{Float64}}
    CI::Vector{Vector{Float64}}
    pardates::Vector{Float64}
    parnames::Vector{String}
    colors::Vector{Symbol}
    plotlayout::Vector{Int}

    Hs::Vector{Float64}
    Hsdates::Vector{Float64}

    power::Matrix{Float64}
    Ω::Vector{Float64}
    specdates::Vector{Float64}

    winddirection::Vector{Float64}
    winddates::Vector{Float64}

    function Figure11()
        @load "dat/anonomous_data.jld"
        @load "dat/wavefits.jld"
        wavetime_days = wavetime ./ 24 ./ 3600
        windtime_days = windtime ./ 24 ./ 3600
        spectrogram, Ω, specdates = PaperAnalysis.compute_spectrogram(waveseries,wavetime_days)

        θ_est = [[f.fits.θ_est[i] for f ∈ fits] for i ∈ 1:9]
        CI = [[1.96 * sqrt(f.fits.varθ[i,i]) for f ∈ fits] for i ∈ 1:9]

        all(f.fits.converged for f in fits) || error("Some fits failed to converge")

        plotlayout = [2,3,5,7,4,6,6,8,8]
        colors = parametercolors()
        parnames = ["\\alpha", "\\omega_p", "\\gamma", "r", "\\phi_m", "\\beta", "\\nu", "\\sigma_l", "\\sigma_r"]
        new(θ_est, CI, wavetime_days[fitstart:fitend], parnames, colors, plotlayout,
            4std.(waveseries), wavetime_days,
            spectrogram, Ω, specdates,
            winddirection,windtime_days)
    end
end

@recipe function f(f::Figure11)
    layout := (4,2)
    colorbar := false
    @series begin
        subplot := 1
        label := L"H_s"
        ylims := (1.5,5.5)
        xticks := (2:4,map(i->"day $i",2:4))
        xlims := (2-10/24, 5)
        f.Hsdates, f.Hs
    end
    @series begin
        seriestype := :heatmap
        subplot := 3
        f.specdates, f.Ω, 10 .* log10.(f.power)
    end
    @series begin
        seriescolor := :black
        subplot := 4
        label := "wind direction"
        f.winddates, mod2pi_nan.(deg2rad.(f.winddirection))
    end
    for i ∈ 1:9
        @series begin
            subplot := f.plotlayout[i]
            if f.plotlayout[i] == 3
                ylims := (0,2)
            end
            yerror := f.CI[i]
            seriestype := :scatter
            markersize := 3
            seriescolor := f.colors[i]
            markerstrokecolor := f.colors[i]
            markerstrokewidth := 1
            label := L"%$(f.parnames[i])"
            xticks := (2:4,map(i->"day $i",2:4))
            xlims := (2-10/24, 5)
            if i == 5
                return f.pardates, mod2pi.(f.θ_est[i])
            else
                return f.pardates, f.θ_est[i]
            end
        end
    end
end


text:
#=
  Test of cumulative constraint in Julia ConstraintSolver.jl

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP, Printf
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

#
# The cumulative constraint is in constraints_utils.jl
#
function cumulative_test(problem, print_solutions=true,all_solutions=false,time_limit=600)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS, # <-
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf, # <-

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>time_limit,

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    n             = problem[:n]
    problem_name  = problem[:name]
    durations     = problem[:durations]
    resources     = problem[:resources]
    start_time_lb = problem[:start_time_domain][1]
    start_time_ub = problem[:start_time_domain][2]
    end_time_lb   = 0
    end_time_ub   = start_time_ub+maximum(durations)
    limit_lb  = problem[:limit_domain][1]
    limit_ub  = problem[:limit_domain][2]

    println("$problem_name (n:$n)")
    println("durations:$durations")
    println("resources:$resources")
    println("start_time_lb:$start_time_lb")
    println("start_time_ub:$start_time_ub")
    println("end_time_lb:$end_time_lb")
    println("end_time_ub:$end_time_ub")
    println("limit_lb:$limit_lb")
    println("limit_ub:$limit_ub")

    @variable(model, start_time_lb <= start_times[1:n] <= start_time_ub, Int)
    @variable(model, end_time_lb   <= end_times[1:n]   <= end_time_ub,   Int)
    @variable(model, limit_lb      <= limit            <= limit_ub,      Int)
    @variable(model, start_time_lb <= max_time         <= end_time_ub,   Int)

    for i in 1:n
        @constraint(model,end_times[i] == start_times[i] + durations[i])
    end
    @constraint(model, end_times .<= max_time)
    cumulative(model, start_times, durations, resources, limit)

    # This don't work: MethodError: no method matching iterate(::Nothing)
    # @objective(model, Min, max_time)
    # @objective(model,Min,limit)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                start_timesx = convert.(Integer,JuMP.value.(start_times; result=sol))
                end_timesx = convert.(Integer,JuMP.value.(end_times; result=sol))
                max_timex = convert.(Integer,JuMP.value.(max_time; result=sol))
                limitx = convert.(Integer,JuMP.value.(limit; result=sol))
                println("start_times:$start_timesx")
                println("durations  :$durations")
                println("end_times  :$end_timesx")
                println("resources  :$resources")
                # println("max_time   :$max_timex")
                println("limit      :$limitx")

                for i in 1:n
                    println("Task $i: $(start_timesx[i])..$(end_timesx[i]) (duration:$(durations[i]) resource:$(resources[i]))")
                end

                time_line = minimum(start_timesx):maximum(end_timesx)
                println("\nTime line: (tasks active (resource))")
                for t in time_line
                    @printf "time %3d: " t
                    total_resources = 0
                    for i in 1:length(start_timesx)
                        if start_timesx[i] <= t && t < start_timesx[i]+durations[i]
                            total_resources += resources[i]
                            @printf " %3d (%2d) " i resources[i]
                        else
                            @printf " %3s  %3s " "" ""
                        end
                    end
                    println("  total_resources: $total_resources")
                end
            end
        end
    else
        println("status:$status")
    end

    return status
end

furniture_moving = Dict(
                        :name               => "Furniture moving",
                        :tasks              => ["piano", "chair", "bed", "table"],
                        :n                  => 4,
                        :durations          => [30,10,15,15],
                        :resources          => [3,1,3,2],
                        :start_time_domain  => [0,60],
                        :limit_domain       => [1,3]
)

# From my Picat model cumulative_decomp.pi
test1 = Dict(
             :name => "cumulative_decomp.pi",
             :tasks => [1,2,3,4,5],
             :n     => 5,
             :durations => [3,9,10,6,2],
             :resources => [1,2,1,1,3],
             :start_time_domain => [1,10],
             :limit_domain => [1,3] # [1,8]
)

# From cumulative_test2.pi
# Originally from (simplified)
# http://stackoverflow.com/questions/25170868/search-heuristics-for-cumulatives
# Note: This is currently (ConstraintSolver.jl 0.5.1) too hard for a 10 min time limit
test2 = Dict(
            :name => "cumulative_test2.pi",
            :n     => 20,
            :tasks => collect(1:20),
            :durations => [318, 246, 797, 238, 251, 279, 987, 847, 426, 787,6, 681,465,46, 3, 427,956,657,113,251],
            :resources => [1 for _ in 1:20],
            :start_time_domain => [1,1000],
            :limit_domain => [1,1]

)

# instances = [furniture_moving,test1,test2]
instances = [furniture_moving,test1]
for instance in instances
    println("Problem $(instance[:name])")
    @time cumulative_test(instance)
    println("\n")
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# Run demonstration
include("conf2022/cuboid_example.jl")

# Create plots and save them in folder results/plots
include("conf2022/plot_results.jl")


text:
function (eos::EnergyEquation{<:Murnaghan1st})(v)
    @unpack v0, b0, b′0, e0 = getparam(eos)
    x, y = b′0 - 1, (v0 / v)^b′0
    return e0 + b0 / b′0 * v * (y / x + 1) - v0 * b0 / x
end
function (eos::EnergyEquation{<:BirchMurnaghan2nd})(v)
    @unpack v0, b0, e0 = getparam(eos)
    f = ToEulerianStrain(v0)(v)
    return e0 + 9b0 * v0 * f^2 / 2
end
function (eos::EnergyEquation{<:BirchMurnaghan3rd})(v)
    @unpack v0, b0, b′0, e0 = getparam(eos)
    f = ToEulerianStrain(v0)(v)
    return e0 + 9b0 * v0 * f^2 / 2 * (1 + f * (b′0 - 4))
end
function (eos::EnergyEquation{<:BirchMurnaghan4th})(v)
    @unpack v0, b0, b′0, b″0, e0 = getparam(eos)
    f = ToEulerianStrain(v0)(v)
    h = b″0 * b0 + b′0^2
    return e0 + 3b0 * v0 / 8 * f^2 * ((9h - 63b′0 + 143) * f^2 + 12f * (b′0 - 4) + 12)
end
function (eos::EnergyEquation{<:PoirierTarantola2nd})(v)
    @unpack v0, b0, e0 = getparam(eos)
    f = ToNaturalStrain(v0)(v)
    return e0 + 9b0 * v0 * f^2 / 2
end
function (eos::EnergyEquation{<:PoirierTarantola3rd})(v)
    @unpack v0, b0, b′0, e0 = getparam(eos)
    f = ToNaturalStrain(v0)(v)
    return e0 + 9b0 * v0 * f^2 / 2 * ((2 - b′0) * f + 1)
end
function (eos::EnergyEquation{<:PoirierTarantola4th})(v)
    @unpack v0, b0, b′0, b″0, e0 = getparam(eos)
    f = ToNaturalStrain(v0)(v)
    h = b″0 * b0 + b′0^2
    return e0 + 9b0 * v0 * f^2 * (3f^2 * (h + 3b′0 + 3) + 4f * (b′0 + 2) + 4)
end
function (eos::EnergyEquation{<:Vinet})(v)
    @unpack v0, b0, b′0, e0 = getparam(eos)
    x, y = 1 - (v / v0)^_⅓, 3 / 2 * (b′0 - 1)
    return e0 + 9b0 * v0 / y^2 * (1 + (x * y - 1) * exp(x * y))
end
function (eos::EnergyEquation{<:AntonSchmidt})(v)
    @unpack v0, b0, b′0, e∞ = getparam(eos)
    n₊₁ = -b′0 / 2 + 1
    x = v / v0
    return e∞ + b0 * v0 / n₊₁ * x^n₊₁ * (log(x) - 1 / n₊₁)
end


text:
# Tests based on the VRT format tutorial found at http://www.gdal.org/gdal_vrttut.html

# To create a VRT dataset that is a clone of an existing dataset use the CreateCopy() method.
driver = GDAL.getdriverbyname("VRT")
src_dataset = GDAL.openshared("data/utmsmall.tif", GDAL.GA_ReadOnly)
vrt_dataset = GDAL.createcopy(driver, "tmp/utmsmall.vrt", src_dataset, 0,
                              C_NULL, progressfunc, C_NULL)

GDAL.setmetadataitem(vrt_dataset, "SourceAgency", "JuliaGeo", "")

GDAL.close(src_dataset)
GDAL.close(vrt_dataset)

# checking the contents of the VRT XML
vrt_xml = readstring("tmp/utmsmall.vrt")
@test contains(vrt_xml, "<VRTDataset rasterXSize=\"100\" rasterYSize=\"100\">")
@test contains(vrt_xml, "<MDI key=\"SourceAgency\">JuliaGeo</MDI>")
@test contains(vrt_xml, "<SourceFilename relativeToVRT=\"0\">data/utmsmall.tif</SourceFilename>")
rm("tmp/utmsmall.vrt")


text:
using ColloidFluctuatingField
using Test

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

max_stars_repo_name: phyjonas/ColloidFluctuatingField


text:
function ethmult(a::Integer, b::Integer)
    r = 0
    while a > 0
        r += b * !even(a)
        a = halve(a)
        b = double(b)
    end
    return r
end

@show ethmult(17, 34)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@testset "ch10_puddle_world2" begin
    dt = 0.1
    # environment
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    # id of landmark must start from 0 with 1 step
    landmarks =
        [Landmark([-4.0, 2.0], 0), Landmark([2.0, -3.0], 1), Landmark([3.0, 3.0], 2)]
    envmap = Map()
    push!(envmap, landmarks)
    world = World(xlim, ylim)
    push!(world, envmap)
    # robot side
    initial_pose = [0.0, 0.0, 0.0]
    estimator = KalmanFilter(envmap, initial_pose)
    agent = EstimatorAgent(0.2, 10.0 * pi / 180, dt, estimator)
    robot = RealRobot(
        initial_pose,
        agent,
        RealCamera(landmarks, distance_bias_rate_stddev = 0.0, direction_bias_stddev = 0.0);
        color = "red",
    )
    push!(world, robot)
    # goal
    goal = Goal(-3.0, -3.0)
    push!(world, goal)
    # puddles
    push!(world, Puddle([-2.0, 0.0], [0.0, 2.0], 0.1))
    push!(world, Puddle([-0.5, -2.0], [2.5, 1.0], 0.1))

    anim = @animate for i = 1:50
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        z = observations(robot.sensor_, robot.pose_; noise = false, bias = false)
        p = draw(world, annota)
        v, ω = decision(agent, z, envmap)
        state_transition(robot, v, ω, dt; move_noise = true, vel_bias_noise = true)
    end
    if GUI
        gif(anim, "ch10_puddle_world2.gif", fps = 20)
    end
end

@testset "ch10_puddle_world4" begin
    dt = 0.1
    # environment
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    # id of landmark must start from 0 with 1 step
    landmarks = [
        Landmark([-4.0, 2.0], 0),
        Landmark([2.0, -3.0], 1),
        Landmark([4.0, 4.0], 2),
        Landmark([-4.0, -4.0], 3),
    ]
    envmap = Map()
    push!(envmap, landmarks)
    world = PuddleWorld(xlim, ylim)
    push!(world, envmap)
    # goal
    goal = Goal(-3.0, -3.0)
    push!(world, goal)
    # robot side
    initial_pose = [2.0, 2.0, 0.0]
    estimator = KalmanFilter(envmap, initial_pose)
    agent = PuddleIgnoreAgent(0.2, 10.0 * pi / 180, dt, estimator, goal)
    robot = RealRobot(
        initial_pose,
        agent,
        RealCamera(landmarks, distance_bias_rate_stddev = 0.0, direction_bias_stddev = 0.0);
        color = "red",
    )
    push!(world, robot)
    # puddles
    push!(world, Puddle([-2.0, 0.0], [0.0, 2.0], 0.1))
    push!(world, Puddle([-0.5, -2.0], [2.5, 1.0], 0.1))

    anim = @animate for i = 1:200
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        # t
        update_status(world)
        z = observations(robot.sensor_, robot.pose_; noise = false, bias = false)
        p = draw(world, annota)

        # t+1
        v, ω = decision(agent, z)
        state_transition(robot, v, ω, dt; move_noise = true, vel_bias_noise = true)
    end
    if GUI
        gif(anim, "ch10_puddle_world4.gif", fps = 10)
    end
end

@testset "ch10_policy_evaluation4" begin
    pe = PolicyEvaluator([0.1, 0.1, pi / 20], Goal(-3.0, -3.0))
    init_value(pe)
    init_policy(pe)
    init_state_transition_probs(pe, 0.1, 10)
end

@testset "ch10_policy_evaluation7" begin
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    world = PuddleWorld(xlim, ylim)
    push!(world, Puddle([-2.0, 0.0], [0.0, 2.0], 0.1))
    push!(world, Puddle([-0.5, -2.0], [2.5, 1.0], 0.1))

    sampling_num = 10
    pe = PolicyEvaluator([0.1, 0.1, pi / 20], Goal(-3.0, -3.0), dt = 0.1)
    init_value(pe)
    init_policy(pe)
    init_state_transition_probs(pe, 0.1, sampling_num)
    init_depth(pe, world, sampling_num)

    Δ = 1e100
    sweep_num = 0
    while Δ > 1.0 && sweep_num < 10
        Δ = policy_evaluation_sweep(pe)
        sweep_num += 1
    end
    println("$(sweep_num): Δ = $(Δ)")
end

@testset "ch10_dp1_value_iteration" begin
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    world = PuddleWorld(xlim, ylim)
    push!(world, Puddle([-2.0, 0.0], [0.0, 2.0], 0.1))
    push!(world, Puddle([-0.5, -2.0], [2.5, 1.0], 0.1))

    sampling_num = 10
    pe = PolicyEvaluator([0.1, 0.1, pi / 20], Goal(-3.0, -3.0), dt = 0.1)
    init_value(pe)
    init_policy(pe)
    init_state_transition_probs(pe, 0.1, sampling_num)
    init_depth(pe, world, sampling_num)

    Δ = 1e100
    sweep_num = 0
    while Δ > 1.0 && sweep_num < 10
        Δ = value_iteration_sweep(pe)
        sweep_num += 1
    end
    println("$(sweep_num): Δ = $(Δ)")
end

@testset "ch10_dp2_policy_iteration" begin
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    world = PuddleWorld(xlim, ylim)
    push!(world, Puddle([-2.0, 0.0], [0.0, 2.0], 0.1))
    push!(world, Puddle([-0.5, -2.0], [2.5, 1.0], 0.1))

    sampling_num = 10
    dp = PolicyEvaluator([0.1, 0.1, pi / 20], Goal(-3.0, -3.0), dt = 0.1)
    init_value(dp)
    init_policy(dp)
    init_state_transition_probs(dp, 0.1, sampling_num)
    init_depth(dp, world, sampling_num)

    action_switch_rate = 1.0
    sweep_num = 0
    while action_switch_rate > 0.05 && sweep_num < 10
        action_switch_rate = policy_iteration_sweep(dp)
        sweep_num += 1
    end
    println("$(sweep_num)th: $(action_switch_rate)")
end

max_stars_repo_name: soblin/JuliaProbo

max_issues_repo_name: soblin/JuliaProbo

max_forks_repo_name: soblin/JuliaProbo


text:
function cvt_to_binary_string(i::Int64)
	return bitstring(i)
end

function cvt_to_binary_string(ui::UInt64)
	return bitstring(ui)
end

function cvt_to_int(binary_string::String)
    return parse(UInt64, binary_string; base=2)
end

function cvt_to_int(ui::UInt64)
    return cvt_to_int(cvt_to_binary_string(ui))
end

function cvt_to_binary_string(bit_array::BitArray{1})
    binary_string = ""
    for i in 1:length(bit_array)
        if bit_array[i]
            binary_string *= "1"
        else
            binary_string *= "0"
        end
    end
    return binary_string
end

function cvt_to_uint(pieces_array::Array{UInt64,1})
    pieces_uint = UInt64(0)
    for piece in pieces_array
        pieces_uint |= piece
    end
    return pieces_uint
end

function cvt_to_uint(binary_string::String)
    return UInt(cvt_to_int(binary_string))
end

function cvt_to_uint(i::Int64)
    return UInt(i)
end

function cvt_to_uint(bit_array::BitArray{1})
    return cvt_to_uint(cvt_to_binary_string(bit_array))
end

function cvt_to_bitarray(pieces_array::Array{UInt64,1})
    return cvt_to_bitarray(cvt_to_uint(pieces_array))
end

function cvt_to_bitarray(i::UInt64)
    return cvt_to_bitarray(cvt_to_binary_string(i))
end

function cvt_to_bitarray(s::String)
    bit_array = falses(64)
    for i = 1:64
        if s[i] == '1'
            bit_array[i] = true
        end
    end
    return bit_array
end

# julia> Bobby.generate_pgn_square_to_uint()
# Dict{String,UInt64} with 64 entries:
#   "d7" => 0x0010000000000000
#   "e4" => 0x0000000008000000
#   ⋮    => ⋮
function gen_pgn_square_to_uint_dict()
    pgn_squares = Dict{String,UInt64}()

    i = 0
    for rank in 8:-1:1
        for file in 1:8
           square = "0"^i * "1" * "0"^(63 - i)
           square_int = cvt_to_int(square)
           pgn_coordinate = string(Char(Int('a')+(file - 1))) * string(rank)
           push!(pgn_squares, pgn_coordinate=>square_int)
           i += 1
        end
    end

    return pgn_squares
end
const PGN2UINT = gen_pgn_square_to_uint_dict()

# julia> squares = Bobby.generate_pgn_square_to_int()
# Dict{String,Int64} with 64 entries:
#   "d7" => 11
#   "a5" => 24
#   ⋮    => ⋮
function gen_pgn_square_to_int_dict()
    pgn_squares = Dict{String,Int64}()
    int_squares = Dict{Int64,String}()

    i = 1
    for rank in 8:-1:1
        for file in 1:8
           pgn_coordinate = string(Char(Int('a')+(file - 1))) * string(rank)
           push!(pgn_squares, pgn_coordinate=>i)
           push!(int_squares, i=>pgn_coordinate)
           i += 1
        end
    end

    return pgn_squares, int_squares
end
const PGN2INT, INT2PGN = gen_pgn_square_to_int_dict()
const A1 = PGN2UINT["a1"]
const B1 = PGN2UINT["b1"]
const C1 = PGN2UINT["c1"]
const D1 = PGN2UINT["d1"]
const E1 = PGN2UINT["e1"]
const F1 = PGN2UINT["f1"]
const G1 = PGN2UINT["g1"]
const H1 = PGN2UINT["h1"]
const A8 = PGN2UINT["a8"]
const B8 = PGN2UINT["b8"]
const C8 = PGN2UINT["c8"]
const D8 = PGN2UINT["d8"]
const E8 = PGN2UINT["e8"]
const F8 = PGN2UINT["f8"]
const G8 = PGN2UINT["g8"]
const H8 = PGN2UINT["h8"]

# julia> Bobby.generate_int_to_uint()
# Dict{Int64,UInt64} with 64 entries:
#   2  => 0x4000000000000000
#   11 => 0x0020000000000000
#   ⋮  => ⋮
function gen_int_to_uint_dict()
    squares_int_uint = Dict{Int64,UInt64}()

    for k in keys(PGN2INT)
        push!(squares_int_uint, PGN2INT[k]=>PGN2UINT[k])
    end

    return squares_int_uint
end
const INT2UINT = gen_int_to_uint_dict()

function gen_uint_to_pgn_dict()
    uint_pgn = Dict{UInt64,String}()

    for k in keys(INT2PGN)
        push!(uint_pgn, INT2UINT[k] => INT2PGN[k])
    end
    return uint_pgn
end
const UINT2PGN = gen_uint_to_pgn_dict()

function change_color(color::String)
    if color == "white"
        return "black"
    else
        return "white"
    end
end


text:
function _get_bounds(T, default_lb, default_ub, vars)
  cell_lb = T[]
  cell_ub = T[]
  for v in vars
    contains(string(v), "InPin") && continue
    contains(string(v), "OutPin") && continue
    hasmetadata(v, ModelingToolkit.VariableOutput) && continue
    lower = hasmetadata(v, VariableLowerBound) ? getmetadata(v, VariableLowerBound) : default_lb
    upper = hasmetadata(v, VariableUpperBound) ? getmetadata(v, VariableUpperBound) : default_ub
    push!(cell_lb, lower)
    push!(cell_ub, upper)
  end
  return cell_lb, cell_ub
end

function print_cell_info(cell, wiring, train_u0)
  println("--------------- MTK Cell ---------------")
  println("in:                         $(wiring.n_in)")
  println("out:                        $(wiring.n_out)")
  println("# neurons:                  $(wiring.n_total)")
  println("# input-neuron synapses:    $(Int(sum(wiring.sens_mask)))")
  println("# neuron-neuron synapses:   $(Int(sum(wiring.syn_mask)))")
  println("# params:                   $(length(cell.p))")
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    eltypes = (N0f8, N0f16, Float32, Float64)        # eltypes of parametric colors
    pctypes = (Gray, RGB, AGray, GrayA, ARGB, RGBA)  # parametric colors
    cctypes = (Gray24, AGray32, RGB24, ARGB32)       # non-parametric colors
    realtypes = (Float16, Float32, Float64, Int)     # types for mixed Normed/Real operations
    # Constructors
    for R in realtypes
        for T in eltypes
            precompile(Tuple{Type{Gray{T}},R})
            precompile(Tuple{Type{AGray{T}},R})
            precompile(Tuple{Type{GrayA{T}},R})
            precompile(Tuple{Type{AGray{T}},R,R})
            precompile(Tuple{Type{GrayA{T}},R,R})
            precompile(Tuple{Type{RGB{T}},R,R,R})
            precompile(Tuple{Type{RGBA{T}},R,R,R})
            precompile(Tuple{Type{RGBA{T}},R,R,R,R})
            precompile(Tuple{Type{ARGB{T}},R,R,R})
            precompile(Tuple{Type{ARGB{T}},R,R,R,R})
        end
        precompile(Tuple{Type{Gray},R})
        precompile(Tuple{Type{AGray},R})
        precompile(Tuple{Type{GrayA},R})
        precompile(Tuple{Type{AGray},R,R})
        precompile(Tuple{Type{GrayA},R,R})
        precompile(Tuple{Type{RGB},R,R,R})
        precompile(Tuple{Type{RGBA},R,R,R})
        precompile(Tuple{Type{RGBA},R,R,R,R})
        precompile(Tuple{Type{ARGB},R,R,R})
        precompile(Tuple{Type{ARGB},R,R,R,R})
        precompile(Tuple{Type{Gray24},R})
        precompile(Tuple{Type{AGray32},R})
        precompile(Tuple{Type{AGray32},R,R})
        precompile(Tuple{Type{RGB24},R,R,R})
        precompile(Tuple{Type{ARGB32},R,R,R})
        precompile(Tuple{Type{ARGB32},R,R,R,R})
    end
    precompile(Tuple{Type{Gray24},UInt32,Type{Val{true}}})
    precompile(Tuple{Type{AGray32},UInt32,Type{Val{true}}})
    precompile(Tuple{Type{RGB24},UInt32,Type{Val{true}}})
    precompile(Tuple{Type{ARGB32},UInt32,Type{Val{true}}})
    # LCHab is used by distinguishable_colors
    for R in (Float32, Float64)
        precompile(Tuple{Type{LCHab},R,R,R})
    end
    # convert
    for T1 in eltypes, T2 in eltypes
        precompile(Tuple{typeof(convert),Type{Gray{T1}},Gray{T2}})
        precompile(Tuple{typeof(convert),Type{RGB{T1}},RGB{T2}})
        precompile(Tuple{typeof(convert),Type{RGB{T1}},Gray{T2}})
        precompile(Tuple{typeof(convert),Type{RGBA{T1}},RGB{T2}})
        precompile(Tuple{typeof(convert),Type{RGBA{T1}},RGBA{T2}})
        precompile(Tuple{typeof(convert),Type{RGBA{T1}},ARGB{T2}})
        precompile(Tuple{typeof(convert),Type{RGBA{T1}},Gray{T2}})
        precompile(Tuple{typeof(convert),Type{ARGB{T1}},RGB{T2}})
        precompile(Tuple{typeof(convert),Type{ARGB{T1}},RGBA{T2}})
        precompile(Tuple{typeof(convert),Type{ARGB{T1}},ARGB{T2}})
        precompile(Tuple{typeof(convert),Type{ARGB{T1}},Gray{T2}})
    end
    # traits
    for T in eltypes, C in pctypes
        precompile(Tuple{typeof(eltype),Type{C}})
        precompile(Tuple{typeof(eltype),Type{C{T}}})
        precompile(Tuple{typeof(alpha),C{T}})
        C <: AbstractGray && precompile(Tuple{typeof(gray),C{T}})
        C <: AbstractRGB || continue
        precompile(Tuple{typeof(red),C{T}})
        precompile(Tuple{typeof(green),C{T}})
        precompile(Tuple{typeof(blue),C{T}})
    end
    for C in cctypes
        precompile(Tuple{typeof(eltype),Type{C}})
        precompile(Tuple{typeof(alpha),C})
        C <: AbstractGray && precompile(Tuple{typeof(gray),C})
        C <: AbstractRGB || continue
        precompile(Tuple{typeof(red),C})
        precompile(Tuple{typeof(green),C})
        precompile(Tuple{typeof(blue),C})
    end
    precompile(Tuple{typeof(eltype),Type{AbstractGray{T} where T<:Fractional}})
    precompile(Tuple{typeof(eltype),Type{AbstractRGB{T} where T<:Fractional}})
    # ccolor typically gets compiled as part of other things
    # hash
    for T in eltypes, C in pctypes
        precompile(Tuple{typeof(hash),Type{C{T}},UInt})
    end
    for C in cctypes
        precompile(Tuple{typeof(hash),Type{C},UInt})
    end
    # rand
    for T in eltypes, C in pctypes
        precompile(Tuple{typeof(rand),Type{C{T}},Tuple{Int}})
        precompile(Tuple{typeof(rand),Type{C{T}},Tuple{Int,Int}})
        precompile(Tuple{typeof(rand),Type{C{T}},Tuple{Int,Int,Int}})
        precompile(Tuple{typeof(rand),Type{C{T}},Int})
        precompile(Tuple{typeof(rand),Type{C{T}},Int,Int})
        precompile(Tuple{typeof(rand),Type{C{T}},Int,Int,Int})
    end
    for C in cctypes
        precompile(Tuple{typeof(rand),Type{C},Tuple{Int}})
        precompile(Tuple{typeof(rand),Type{C},Tuple{Int,Int}})
    end
    precompile(Tuple{typeof(rand),Type{Gray{Bool}},Tuple{Int}})
    precompile(Tuple{typeof(rand),Type{Gray{Bool}},Tuple{Int,Int}})
    precompile(Tuple{typeof(rand),Type{Gray{Bool}},Tuple{Int,Int,Int}})
    # show
    for IO in (IOBuffer, IOContext{IOBuffer}, IOContext{Base.TTY})
        for T in eltypes, C in pctypes
            precompile(Tuple{typeof(show),IO,C{T}})
        end
        for C in cctypes
            precompile(Tuple{typeof(show),IO,C})
        end
        precompile(Tuple{typeof(show),IO,Gray{Bool}})
    end
    # FIXME the following do not yet "work", meaning you can issue these precompile directives but no value comes of it.
    # Possibly this is https://github.com/JuliaLang/julia/pull/32705, but the actual cause is unknown.
    # # permutedims
    # for T in eltypes, C in pctypes, n = 2:4
    #     precompile(permutedims, (Array{C{T},n}, Tuple{Vararg{Int,n}}))
    # end
    # # broadcast
    # for T1 in eltypes, T2 in eltypes, n = 2:4
    #     precompile(broadcast, (Type{T1}, Array{Gray{T2},n}))
    #     precompile(broadcast, (Type{Gray{T1}}, Array{T2,n}))
    #     precompile(broadcast, (Type{Gray{T1}}, Array{Gray{T2},n}))
    #     precompile(broadcast, (Type{RGB{T1}}, Array{RGB{T2},n}))
    #     precompile(Broadcast.materialize, (Broadcast.Broadcasted{Broadcast.DefaultArrayStyle{n},Nothing,Type{T1},Tuple{Array{Gray{T2},n}}},))
    #     precompile(Broadcast.materialize, (Broadcast.Broadcasted{Broadcast.DefaultArrayStyle{n},Nothing,Type{Gray{T1}},Tuple{Array{T2,n}}},))
    #     precompile(Broadcast.materialize, (Broadcast.Broadcasted{Broadcast.DefaultArrayStyle{n},Nothing,Type{Gray{T1}},Tuple{Array{Gray{T2},n}}},))
    #     precompile(Broadcast.materialize, (Broadcast.Broadcasted{Broadcast.DefaultArrayStyle{n},Nothing,Type{RGB{T1}},Tuple{Array{RGB{T2},n}}},))
    # end

    # mapc, reducec, and mapreducec are not really precompilable,
    # since they will be specialized for f
end


text:
const Counter = let
    init() = 0
    op(n, _) = n + 1
    value(n) = n
    (; init, op, value)
end

struct FrequencyAnalysis
    options::Dict{Symbol, Any}
end

to_nothings(v) = fill(nothing, axes(v))

function (f::FrequencyAnalysis)(entry::Entry)
    positional = copy(entry.positional)
    push!(positional, map(to_nothings, first(positional)))
    labels = copy(entry.labels)
    labels[length(positional)] = "count"
    augmented_entry = Entry(entry; positional, labels)
    return groupreduce(Counter, augmented_entry)
end

"""
    frequency()

Compute a frequency table of the arguments.
"""
frequency() = transformation(FrequencyAnalysis(Dict{Symbol, Any}()))


text:
# module Bukdu

type Conn
    status::Int
    resp_header::Dict{String,String}
    resp_body::Any
    params::Assoc
    query_params::Assoc
    private::Dict{Symbol,String}
    assigns::Dict{Symbol,String}

    function Conn(status::Int, resp_header::Dict{String,String}, resp_body::Any)
        new(status, resp_header, resp_body, Assoc(), Assoc(), Dict{Symbol,String}(), Dict{Symbol,String}())
    end

    function Conn(status::Int, resp_header::Dict{String,String}, resp_body::Any, params::Assoc, query_params::Assoc, private::Dict{Symbol,String}, assigns::Dict{Symbol,String})
        new(status, resp_header, resp_body, params, query_params, private, assigns)
    end
end

const CONN_NOT_FOUND = Conn(404, Dict("Content-Type"=>"text/html"), "not found")


text:
# Ref:https://machinelearningmastery.com/least-squares-generative-adversarial-network/

using Base.Iterators: partition
using Flux
using Flux.Optimise: update!
using Flux: mse
using Images
using MLDatasets
using Statistics
using Parameters: @with_kw
using Random
using Printf
using CUDAapi
using Zygote
if has_cuda()		# Check if CUDA is available
    @info "CUDA is on"
    import CuArrays		# If CUDA is available, import CuArrays
    CuArrays.allowscalar(false)
end

@with_kw struct HyperParams
    batch_size::Int = 128
    latent_dim::Int = 100
    epochs::Int = 20
    verbose_freq::Int = 1000
    output_x::Int = 6        # No. of sample images to concatenate along x-axis
    output_y::Int = 6        # No. of sample images to concatenate along y-axis
    lr_dscr::Float64 = 0.0002
    lr_gen::Float64 = 0.0002
end

function generator(args)
    n_nodes = 256*7*7
    return Chain(Dense(args.latent_dim, n_nodes),BatchNorm(n_nodes),x -> leakyrelu.(x,0.2f0), x-> reshape(x, 7,7,256, size(x,2)),
                 ConvTranspose((4,4), 256=>128, stride=(2,2), pad=(1,1)), BatchNorm(128), x -> leakyrelu.(x,0.2f0),
                 ConvTranspose((4,4), 128=>64, stride=(2,2), pad = (1,1)), BatchNorm(64), x -> leakyrelu.(x,0.2f0),
                 Conv((7,7), 64=>1, stride=1, pad=3, tanh)) |> gpu
end

function discriminator(args)
    return  Chain(Conv((4,4), 1=>64, stride=2, pad=1), BatchNorm(64), x -> leakyrelu.(x,0.2f0),
                  Conv((4,4), 64=>128, stride=2, pad=1), BatchNorm(64), x -> leakyrelu.(x,0.2f0),
                  x->reshape(x, :, size(x,4)), Dense(6272, 1, sigmoid))|> gpu
end

function load_data(hparams)
    # Load MNIST dataset
    images, labels = MLDatasets.MNIST.traindata(Float32)
    # Normalize to [-1, 1] and convert it to WHCN
    image_tensor = permutedims(reshape(@.(2f0 * images - 1f0), 28, 28, 1, :), (2, 1, 3, 4))
    # Partition into batches
    data = [image_tensor[:,:,:,r] |> gpu for r in partition(1:60000, hparams.batch_size)]
    return data
end

# Loss functions
function discr_loss(real_output, fake_output)
    real_loss = mean(mse.(real_output, 1f0))
    fake_loss = mean(mse.(fake_output, 0f0))
    return (real_loss + fake_loss)
end

generator_loss(fake_output) = mean(mse.(fake_output, 1f0))

function train_discr(discr, original_data, fake_data, opt_discr)
    ps = Flux.params(discr)
    loss, back = Zygote.pullback(ps) do
                      discr_loss(discr(original_data), discr(fake_data))
    end
    grads = back(1f0)
    update!(opt_discr, ps, grads)
    return loss
end

Zygote.@nograd train_discr

function train_gan(gen, discr, original_data, opt_gen, opt_discr, hparams)
    noise = randn!(similar(original_data, (hparams.latent_dim, hparams.batch_size))) |> gpu
    loss = Dict()
    ps = Flux.params(gen)
    loss["gen"], back = Zygote.pullback(ps) do
                          fake_ = gen(noise)
                          loss["discr"] = train_discr(discr, original_data, fake_, opt_discr)
                          generator_loss(discr(fake_))
    end
    grads = back(1f0)
    update!(opt_gen, ps, grads)
    return loss
end

function create_output_image(gen, fixed_noise, hparams)
    @eval Flux.istraining() = false
    fake_images = @. cpu(gen(fixed_noise))
    @eval Flux.istraining() = true
    image_array = dropdims(reduce(vcat, reduce.(hcat, partition(fake_images, hparams.output_y))); dims=(3, 4))
    image_array = @. Gray(image_array + 1f0) / 2f0
    return image_array
end

function train()
    hparams = HyperParams()

    data = load_data(hparams)

    fixed_noise = [randn(hparams.latent_dim, 1) |> gpu for _=1:hparams.output_x*hparams.output_y]

    # Discriminator
    dscr = discriminator(hparams)

    # Generator
    gen =  generator(hparams) |> gpu

    # Optimizers
    opt_dscr = ADAM(hparams.lr_dscr)
    opt_gen = ADAM(hparams.lr_gen)

    isdir("output")||mkdir("output")

    # Training
    train_steps = 0
    for ep in 1:hparams.epochs
        @info "Epoch $ep"
        for x in data
            # Update discriminator and generator
            loss = train_gan(gen, dscr, x, opt_gen, opt_dscr, hparams)

            if train_steps % hparams.verbose_freq == 0
                @info("Train step $(train_steps), Discriminator loss = $(loss["discr"]), Generator loss = $(loss["gen"])")
                # Save generated fake image
                output_image = create_output_image(gen, fixed_noise, hparams)
                save(@sprintf("output/gan_steps_%06d.png", train_steps), output_image)
            end
            train_steps += 1
        end
    end

    output_image = create_output_image(gen, fixed_noise, hparams)
    save(@sprintf("output/gan_steps_%06d.png", train_steps), output_image)
end

cd(@__DIR__)
train()


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    compute2dcoords

"""
    compute2dcoords(mol::MolGraph) -> InternalCoordinates

Compute 2D coordinates of the molecular graph.
"""
function compute2dcoords(mol::MolGraph)
    graph = mol.graph
    fragments = InternalCoordinates[]

    # Extract scaffolds
    scaffoldnodes = biconnected_component(graph) # not 2-edge connected
    scaffolds = nodesubgraph(scaffoldnodes)
    for scaffold in connected_component(scaffolds)
        if is_outerplanar(nodesubgraph(scaffold))
            push!(fragments, outerplaner_embed2d(scaffold))
        else
            push!(fragments, cartesian_embed2d(scaffold))
        end
        # TODO: macrocycle templates
    end

    # Chains
    chainnodes = setdiff(nodeset(graph), scaffoldnodes)
    chains = nodesubgraph(chainnodes)
    for chain in connected_component(chains)
        push!(fragments, chain_embed2d(chain))
    end

    # TODO: Merge fragments

    # TODO: Boundary check and avoid overlap

    return coords
end

"""
    chain_embed2d(graph::UndirectedGraph; kwargs...) -> InternalCoords

Return a 2D embedding of chain(=tree graph).
"""
function chain_embed2d(graph)
    nodes = Set(chains)
    while !isempty(nodes)
        longest = longestpath(chain)
        push!(path, chain(longest))
        setdiff!(nodes, longest)
    end
end

"""
    is_outerplanar(mol::GraphMol) -> Bool

Return whether the graph is outerplanar or not
"""
function is_outerplanar(mol::GraphMol)
    # TODO: generalize and move to graph.planarity
    # The given graph should be a biconnected graph

    # A biconnected braph G is K4 minor if the graphs have 2n-3 and more edges.
    edgecount(mol) > nodecount(mol) * 2 - 3 && return false

    # A biconnected graph G is K2,3 minor if and only if any pair of basic
    # cycles has two or more edges in common.
    seen = Set{Set{Int}}()
    for i in findall(mol[bond_ringmem] .> 1)
        m = mol[bond_ringmem][i]
        m in seen && return false
        push!(seen, m)
    end
    return true
end

"""
    outerplaner_embed2d(graph::UndirectedGraph; kwargs...) -> InternalCoords

Return a 2D embedding of the outerplanar graph.

A 2D embedding of an outerplanar graph can be easily determined by Hamiltonian
path traversal.
"""
function outerplaner_embed2d(mol::VecterMol)
    coords = internalcoords()
    root = pop!(nodekeys(mol))
    stack = [root]
    pred = Dict(root => root)
    while !isempty(stack)
        c = pop!(stack)

        p1 = pred[c]
        p2 = pred[p1]
        p3 = pred[p2]
        r = pop!(copy(mol[:RingMem][c]))
        ringsize = pop!(copy(mol[:RingSize][c]))

        for (nbr, bond) in neighbors(mol, c)
            if length(mol[bond_ringmem][bond]) == 1 && !(nbr in keys(pred))
                pred[nbr] = c
            end
        end

        if !isempty(isec3)
            # elongate ring
            angle = 1 - 2 / length(rings[isec3[1]])
            setnodes!(coords, c, [p1, p2, p3], [1,0, angle, 0.0])
        elseif !isempty(isec2)
            # branch ring
            angle = 1 - 2 / length(rings[isec2[1]])
            setnodes!(coords, c, [p1, p2, p3], [1,0, angle, 1.0])
        end
        # TODO: spiro
        nextlevels = []
        pop!(nbrs, pred[c])
        backtracked = length(nextlevels) == 0
        append!(stack, nextlevels)
        push!(done, c)
    end
    coords3d = cartesian(zmatrix)
    coords3d[:, 1:2]
end

"""
    cartesian_embed2d(graph::UndirectedGraph; kwargs...) -> Cartesian2D

Cartesian 2D embedding
"""
function cartesian_embed2d(graph)

end

function coords2d(mol::GraphMol, root)
    # Requires :Topology
    zmatrix = [
        -2 nothing nothing nothing nothing nothing nothing;
        -1 -2 1.0 nothing nothing nothing nothing;
        0 -1 1.0 -2 2/3 nothing 1
    ]
    stack = [root]
    done = []
    pred = Dict(root => 0, 0 => -1, -1 => -2)
    ringmap = Dict(0 => Set(), -1 => Set(), -2 => Set())
    cyclemap = Dict(i => c for (i, c) in enumerate(mol[:Cycle]))
    merge!(ringmap, cyclemap)
    rings = mol.annotation[:Topology].rings
    backtracked = false
    while length(stack) > 0
        c = pop!(stack)
        if c in done
            continue
        end
        p1 = pred[c]
        p2 = pred[p1]
        p3 = pred[p2]

        isec3 = collect(intersect(ringmap[c], ringmap[p3]))
        isec2 = collect(intersect(ringmap[c], ringmap[p2]))
        if !isempty(isec3)
            # elongate ring
            angle = 1 - 2 / length(rings[isec3[1]])
            zmatrix = vcat(zmatrix, [c p1 1.0 p2 angle p3 0.0])
        elseif !isempty(isec2)
            # branch ring
            angle = 1 - 2 / length(rings[isec2[1]])
            zmatrix = vcat(zmatrix, [c p1 1.0 p2 angle p3 1.0])
        elseif backtracked
            # branch chain
            zmatrix = vcat(zmatrix, [c p1 1.0 p2 2/3 p3 0.0])
        else
            # elongate chain
            zmatrix = vcat(zmatrix, [c p1 1.0 p2 2/3 p3 1.0])
        end
        # TODO: cis-trans specified
        # TODO: triple bond
        # TODO: spiro
        nextlevels = []
        for n in keys(neighbors(mol, c))
            if n == pred[c]
                continue
            end
            pred[n] = c
            if !(n in done)
                push!(nextlevels, n)
            end
        end
        backtracked = length(nextlevels) == 0
        append!(stack, nextlevels)
        push!(done, c)
    end
    coords3d = cartesian(zmatrix)
    coords3d[:, 1:2]
end


text:
using Combinatorics: combinations

test = """
    [[[0,[5,8]],[[1,7],[9,6]]],[[4,[1,2]],[[1,4],2]]]
    [[[5,[2,8]],4],[5,[[9,9],0]]]
    [6,[[[6,2],[5,6]],[[7,6],[4,7]]]]
    [[[6,[0,7]],[0,9]],[4,[9,[9,0]]]]
    [[[7,[6,4]],[3,[1,3]]],[[[5,5],1],9]]
    [[6,[[7,3],[3,2]]],[[[3,8],[5,7]],4]]
    [[[[5,4],[7,7]],8],[[8,3],8]]
    [[9,3],[[9,9],[6,[4,9]]]]
    [[2,[[7,7],7]],[[5,8],[[9,3],[0,2]]]]
    [[[[5,2],5],[8,[3,7]]],[[5,[7,5]],[4,4]]]
    """

function parse_input(input::AbstractString)
    lines = split(input, '\n'; keepempty = false)
    return eval(Expr(:vect, Meta.parse.(lines)...))
end

struct ArrayRef
    v::Array
    i
end

Base.lastindex(ref::ArrayRef) = lastindex(ref[])
Base.getindex(ref::ArrayRef) = ref.v[ref.i]
Base.getindex(ref::ArrayRef, i) = ArrayRef(ref[], i)
Base.setindex!(ref::ArrayRef, x) = ref.v[ref.i] = x

reffirst(v) = v[] isa Vector ? reffirst(v[1]) : v
reflast(v) = v[] isa Vector ? reflast(v[end]) : v

magnitude(v::Integer) = v
magnitude(v::AbstractVector) = 3magnitude(v[1]) + 2magnitude(v[2])

explode!(v) = explode!(v, 0, Ref(0), Ref(0))
explode!(v::Integer, kwargs...) = (v, false)
function explode!(v::AbstractVector, depth, left, right)
    if depth == 4
        left[] += v[1]
        right[] += v[2]
        return (0, true)
    end
    (v1r, ok) = explode!(v[1], depth + 1, left, reffirst(ArrayRef(v, 2)))
    ok && return ([v1r, v[2]], true)
    (v2r, ok) = explode!(v[2], depth + 1, reflast(ArrayRef(v, 1)), right)
    ok && return ([v[1], v2r], true)
    return (v, false)
end

split!(i::Integer) = i >= 10 ? ([i ÷ 2, i - i ÷ 2], true) : (i, false)
function split!(v::AbstractVector)
    (v1s, ok) = split!(v[1])
    ok && return ([v1s, v[2]], true)
    (v2s, ok) = split!(v[2])
    ok && return ([v[1], v2s], true)
    return (v, false)
end

function reduce!(v::AbstractVector)
    while true
        (v, ok) = explode!(v)
        ok && continue
        (v, ok) = split!(v)
        !ok && break
    end
    return v
end

add(v1, v2) = reduce!(deepcopy([v1, v2]))

q1(v) = magnitude(foldl(add, v))

function q2(v)
    pairs = combinations(v, 2)
    return map(pairs) do (a, b)
        s1 = add(a, b)
        s2 = add(b, a)
        return maximum(magnitude.([s1, s2]))
    end |> findmax |> first
end

@assert magnitude([[9,1], [1,9]]) == 129

@assert explode!([[[[[9,8],1],2],3],4]) == ([[[[0,9],2],3],4], true)
@assert explode!([7,[6,[5,[4,[3,2]]]]]) == ([7,[6,[5,[7,0]]]], true)
@assert explode!([[6,[5,[4,[3,2]]]],1]) == ([[6,[5,[7,0]]],3], true)
@assert explode!([[3,[2,[1,[7,3]]]],[6,[5,[4,[3,2]]]]]) ==
    ([[3,[2,[8,0]]],[9,[5,[4,[3,2]]]]], true)
@assert explode!([[3,[2,[8,0]]],[9,[5,[4,[3,2]]]]]) ==
    ([[3,[2,[8,0]]],[9,[5,[7,0]]]], true)

@assert split!(11) == ([5, 6], true)

@assert add([[[[4,3],4],4],[7,[[8,4],9]]], [1,1]) ==
    [[[[0,7],4],[[7,8],[6,0]]],[8,1]]
@assert add([[[0,[4,5]],[0,0]],[[[4,5],[2,6]],[9,5]]],
            [7,[[[3,7],[4,3]],[[6,3],[8,8]]]]) ==
                [[[[4,0],[5,4]],[[7,7],[6,0]]],[[8,[7,7]],[[7,9],[5,0]]]]

@assert foldl(add, [[i,i] for i in 1:4]) == [[[[1,1],[2,2]],[3,3]],[4,4]]
@assert foldl(add, [[i,i] for i in 1:5]) == [[[[3,0],[5,3]],[4,4]],[5,5]]
@assert foldl(add, [[i,i] for i in 1:6]) == [[[[5,0],[7,4]],[5,5]],[6,6]]
@assert foldl(add, [[[[0,[4,5]],[0,0]],[[[4,5],[2,6]],[9,5]]],
                    [7,[[[3,7],[4,3]],[[6,3],[8,8]]]],
                    [[2,[[0,8],[3,4]]],[[[6,7],1],[7,[1,6]]]],
                    [[[[2,4],7],[6,[0,5]]],[[[6,8],[2,8]],[[2,1],[4,5]]]],
                    [7,[5,[[3,8],[1,4]]]],
                    [[2,[2,2]],[8,[8,1]]],
                    [2,9],
                    [1,[[[9,3],9],[[9,0],[0,7]]]],
                    [[[5,[7,4]],7],1],
                    [[[[4,2],2],6],[8,7]]]) ==
                        [[[[8,7],[7,7]],[[8,6],[7,7]]],[[[0,7],[6,6]],[8,7]]]

let v = parse_input(test)
    @assert foldl(add, v) ==
        [[[[6,6],[7,6]],[[7,7],[7,0]]],[[[7,7],[7,7]],[[7,8],[9,9]]]]
    @assert q1(v) == 4140
    @assert q2(v) == 3993
end

v = parse_input(read("day18.in", String))
println("Q1: ", q1(v))
println("Q2: ", q2(v))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
include("utils.jl")

@testset "Product manifold" begin
    M1 = Sphere(2)
    M2 = Euclidean(2)
    Mse = ProductManifold(M1, M2)
    @test Mse == M1 × M2

    types = [Vector{Float64},
             SizedVector{5, Float64},
             MVector{5, Float64},
             Vector{Float32},
             SizedVector{5, Float32},
             MVector{5, Float32}]

    retraction_methods = [Manifolds.ProductRetraction(ManifoldsBase.ExponentialRetraction(), ManifoldsBase.ExponentialRetraction())]
    inverse_retraction_methods = [Manifolds.InverseProductRetraction(ManifoldsBase.LogarithmicInverseRetraction(), ManifoldsBase.LogarithmicInverseRetraction())]

    reshapers = (Manifolds.StaticReshaper(), Manifolds.ArrayReshaper())

    for T in types, reshaper in reshapers
        if reshaper == reshapers[2] && T != Vector{Float64}
            continue
        end
        shape_se = Manifolds.ShapeSpecification(reshaper, M1, M2)
        @testset "Type $T" begin
            pts_base = [convert(T, [1.0, 0.0, 0.0, 0.0, 0.0]),
                        convert(T, [0.0, 1.0, 0.0, 1.0, 0.0]),
                        convert(T, [0.0, 0.0, 1.0, 0.0, 0.1])]
            pts = map(p -> Manifolds.ProductArray(shape_se, p), pts_base)
            distr_M1 = Manifolds.uniform_distribution(M1, pts_base[1][1:3])
            distr_M2 = Manifolds.projected_distribution(M2, Distributions.MvNormal(zero(pts_base[1][4:5]), 1.0))
            distr_tv_M1 = Manifolds.normal_tvector_distribution(M1, pts_base[1][1:3], 1.0)
            distr_tv_M2 = Manifolds.normal_tvector_distribution(M2, pts_base[1][4:5], 1.0)
            test_manifold(Mse,
                          pts;
                          test_reverse_diff = isa(T, Vector),
                          test_musical_isomorphisms = true,
                          retraction_methods = retraction_methods,
                          inverse_retraction_methods = inverse_retraction_methods,
                          test_mutating_rand = isa(T, Vector),
                          point_distributions = [Manifolds.ProductPointDistribution(distr_M1, distr_M2)],
                          tvector_distributions = [Manifolds.ProductFVectorDistribution(distr_tv_M1, distr_tv_M2)])
        end
    end

    M3 = Manifolds.Rotations(2)
    Mser = ProductManifold(M1, M2, M3)
    shape_ser = Manifolds.ShapeSpecification(reshapers[1], M1, M2, M3)

    @testset "ShapeSpecification detailed reshapers" begin
        shape_ser_3_test = Manifolds.ShapeSpecification((reshapers[1], reshapers[2], reshapers[1]), M1, M2, M3)
        @test shape_ser_3_test.reshapers[1] == reshapers[1]
        @test shape_ser_3_test.reshapers[2] == reshapers[2]
        @test shape_ser_3_test.reshapers[3] == reshapers[1]
    end

    @testset "high-dimensional product manifold" begin
        Mhigh = Euclidean(100000)
        shape_high = Manifolds.ShapeSpecification(reshapers[2], M1, Mhigh)
        a = Manifolds.ProductArray(shape_high, collect(1.0:100003.0))
        b = a.parts[2].^2
        @test b[1] == 16
    end

    @test submanifold(Mser, 2) == M2
    @test submanifold(Mser, Val((1, 3))) == M1 × M3
    @test submanifold(Mser, 2:3) == M2 × M3
    @test submanifold(Mser, [1, 3]) == M1 × M3
    @inferred submanifold(Mser, Val((1, 3)))

    # testing the slower generic constructor
    Mprod4 = ProductManifold(M2, M2, M2, M2)
    shape4 = Manifolds.ShapeSpecification(reshapers[1], Mprod4.manifolds...)
    data_x4 = rand(8)
    @test Manifolds.ProductArray(shape4, data_x4).data === data_x4

    pts_sphere = [[1.0, 0.0, 0.0],
                  [0.0, 1.0, 0.0],
                  [0.0, 0.0, 1.0]]
    pts_r2 = [[0.0, 0.0],
              [1.0, 0.0],
              [0.0, 0.1]]
    angles = (0.0, π/2, 2π/3)
    pts_rot = [[cos(ϕ) sin(ϕ); -sin(ϕ) cos(ϕ)] for ϕ in angles]
    pts = [Manifolds.prod_point(shape_ser, p[1], p[2], p[3]) for p in zip(pts_sphere, pts_r2, pts_rot)]
    test_manifold(Mser,
                  pts,
                  test_forward_diff = false,
                  test_reverse_diff = false)

    @testset "prod_point" begin
        shape_se = Manifolds.ShapeSpecification(reshapers[1], M1, M2)

        Ts = SizedVector{3, Float64}
        Tr2 = SizedVector{2, Float64}
        T = SizedVector{5, Float64}
        pts_base = [convert(T, [1.0, 0.0, 0.0, 0.0, 0.0]),
                    convert(T, [0.0, 1.0, 0.0, 1.0, 0.0]),
                    convert(T, [0.0, 0.0, 1.0, 0.0, 0.1])]
        pts = map(p -> Manifolds.ProductArray(shape_se, p), pts_base)
        pts_sphere = [convert(Ts, [1.0, 0.0, 0.0]),
                      convert(Ts, [0.0, 1.0, 0.0]),
                      convert(Ts, [0.0, 0.0, 1.0])]
        pts_r2 = [convert(Tr2, [0.0, 0.0]),
                  convert(Tr2, [1.0, 0.0]),
                  convert(Tr2, [0.0, 0.1])]
        pts_prod = [Manifolds.prod_point(shape_se, p[1], p[2]) for p in zip(pts_sphere, pts_r2)]
        for p in zip(pts, pts_prod)
            @test isapprox(Mse, p[1], p[2])
        end
        for p in zip(pts_sphere, pts_r2, pts_prod)
            @test isapprox(M1, p[1], submanifold_component(p[3], 1))
            @test isapprox(M2, p[2], submanifold_component(p[3], 2))
        end
    end

    @testset "ProductRepr" begin
        Ts = SizedVector{3, Float64}
        Tr2 = SizedVector{2, Float64}
        pts_sphere = [convert(Ts, [1.0, 0.0, 0.0]),
                      convert(Ts, [0.0, 1.0, 0.0]),
                      convert(Ts, [0.0, 0.0, 1.0])]
        pts_r2 = [convert(Tr2, [0.0, 0.0]),
                  convert(Tr2, [1.0, 0.0]),
                  convert(Tr2, [0.0, 0.1])]

        pts = [ProductRepr(p[1], p[2]) for p in zip(pts_sphere, pts_r2)]
        test_manifold(Mse,
                      pts,
                      test_musical_isomorphisms = true,
                      test_tangent_vector_broadcasting = false,
                      test_forward_diff = false,
                      test_reverse_diff = false)
    end
end


text:
function ob(r,side,y)
    x = 0
    Flg = true
    while isborder(r,HorizonSide(side)) && !isborder(r,HorizonSide((side + 1) % 4))
        move!(r,HorizonSide((side + 1) % 4))
        x += 1
    end
    if isborder(r,HorizonSide(side)) && isborder(r,HorizonSide((side + 1) % 4))
        Flg = false
    end
    if !isborder(r,HorizonSide(side))
        move!(r,HorizonSide(side))
        y += 1
    end
    while !isborder(r,HorizonSide(side)) && isborder(r,HorizonSide((side + 3) % 4))
        move!(r,HorizonSide(side))
        y += 1
    end
    for i = 1:x
        move!(r,HorizonSide((side+3) % 4))
    end
    return Flg,y
end

function Vozvrat(r,side,y)
    while y!=0
        if !isborder(r,HorizonSide((side + 2) % 4))
            move!(r,HorizonSide(((side + 2) % 4)))
        else
            x = 0
            while isborder(r,HorizonSide((side + 2) % 4)) && !isborder(r,HorizonSide((side + 1) % 4))
                move!(r,HorizonSide((side + 1) % 4))
                x += 1
            end
            move!(r,HorizonSide((side + 2) % 4))
            while isborder(r,HorizonSide((side + 3) % 4))
                move!(r,HorizonSide((side + 2) % 4))
                y-=1
            end
            for i = 1:x
                move!(r,HorizonSide((side + 3) % 4))
            end
        end
        y -= 1
    end
end

function go(r,side)
    y=0
    Flg=true
    while !isborder(r,HorizonSide(side)) || Flg
        if !isborder(r,HorizonSide(side))
            move!(r,HorizonSide(side))
            y+=1
        end
        if isborder(r,HorizonSide(side))
            Flg,y=ob(r,side,y)
        end
    end
    putmarker!(r)
    Vozvrat(r,side,y)
end

function task_11(r)
    for side = 0:3
        go(r,side)
    end
end

task_11(r)


text:
include("utils.jl")

@testset "Cholesky Space" begin
    M = Manifolds.CholeskySpace(3)

    types = [
        Matrix{Float64},
        MMatrix{3, 3, Float64},
        Matrix{Float32},
    ]
    for T in types
        A(α) = [1. 0. 0.; 0. cos(α) sin(α); 0. -sin(α) cos(α)]
        ptsF = [#
            cholesky([1. 0. 0.; 0. 1. 0.; 0. 0. 1]).L,
            cholesky([2. 0. 0.; 0. 2. 0.; 0. 0. 1]).L,
            cholesky(A(π/6) * [1. 0. 0.; 0. 2. 0.; 0. 0. 1] * transpose(A(π/6))).L,
            ]
        pts = [convert(T, a) for a in ptsF]
        test_manifold(M, pts;
            test_injectivity_radius = false,
            test_vector_transport = true,
            test_forward_diff = false,
            test_reverse_diff = false,
            exp_log_atol_multiplier = 8.0,
        )
    end
    @testset "Test Error cases in is_manifold_point and is_tangent_vector" begin
        pt1f = zeros(2,3); # wrong size
        pt2f = [1. 0. 0.; 0. -1. 0.; 0. 0. 1.]; # nonpos diag
        pt3f = [2. 0. 1.; 0. 1. 0.; 0. 0. 4.]; # no lower and nonsym
        pt4 = [2. 0. 0.; 1. 2. 0.; 0. 0. 4.]
        @test !is_manifold_point(M, pt1f)
        @test_throws DomainError is_manifold_point(M, pt1f, true)
        @test !is_manifold_point(M, pt2f)
        @test_throws DomainError is_manifold_point(M, pt2f, true)
        @test !is_manifold_point(M,pt3f)
        @test_throws DomainError is_manifold_point(M,pt3f, true)
        @test is_manifold_point(M, pt4)
        @test !is_tangent_vector(M,pt3f, pt1f)
        @test_throws DomainError is_tangent_vector(M,pt3f, pt1f, true)
        @test !is_tangent_vector(M,pt4, pt1f)
        @test_throws DomainError is_tangent_vector(M,pt4, pt1f, true)
        @test !is_tangent_vector(M,pt4, pt3f)
        @test_throws DomainError is_tangent_vector(M,pt4, pt3f, true)
        @test is_tangent_vector(M,pt4, pt2f)
    end
end


text:
# This file was generated, do not modify it. # hide
using MLJ, MLJBase, PrettyPrinting, DataFrames, Statistics

X, y = @load_reduced_ames
@show size(X)
first(X, 3) |> pretty

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
using Test
using NQCBase
using CubeLDFAModel
using PyCall

py"""
import sys
sys.path.insert(0, ".")
"""

cube = pyimport("cube")

filename = "test.cube"

c = CubeLDFAModel.Cube(filename)
cube_object = cube.cube()
cube_object.read(filename)

@testset "Compare to cube.py" begin
    @test c.density == cube_object."density".reshape(cube_object.x_len, cube_object.y_len, cube_object.z_len)

    r = c.origin .+ 1e-8
    @test c(r) ≈ cube_object(au_to_ang.(r)...)

    for _=1:10
        r = c.origin
        r += rand() * c.cell.vectors[:,1]
        r += rand() * c.cell.vectors[:,2]
        r += rand() * c.cell.vectors[:,3]

        @test c(r) ≈ cube_object(au_to_ang.(r)...)
    end

    r = c.origin + c.cell.vectors[:,1] + c.cell.vectors[:,2] + c.cell.vectors[:,3] .- 1e-8
    @test c(r) ≈ cube_object(au_to_ang.(r)...)
end


text:
visualize_search(init_hist, "red")
visualize_search(x_hist)

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using Test
using BaremetalPi

################################################################################
#                                  C Macros
################################################################################

# Test if the C Macros for the SPI protocol are being computed correctly.
@testset "C Macros - SPI" begin
    @test BaremetalPi.SPI_IOC_RD_MODE          == 0x80016b01
    @test BaremetalPi.SPI_IOC_WR_MODE          == 0x40016b01
    @test BaremetalPi.SPI_IOC_RD_LSB_FIRST     == 0x80016b02
    @test BaremetalPi.SPI_IOC_WR_LSB_FIRST     == 0x40016b02
    @test BaremetalPi.SPI_IOC_RD_BITS_PER_WORD == 0x80016b03
    @test BaremetalPi.SPI_IOC_WR_BITS_PER_WORD == 0x40016b03
    @test BaremetalPi.SPI_IOC_RD_MAX_SPEED_HZ  == 0x80046b04
    @test BaremetalPi.SPI_IOC_WR_MAX_SPEED_HZ  == 0x40046b04
    @test BaremetalPi.SPI_IOC_MESSAGE(1)       == 0x40206b00
    @test BaremetalPi.SPI_IOC_MESSAGE(2)       == 0x40406b00
    @test BaremetalPi.SPI_IOC_MESSAGE(3)       == 0x40606b00
    @test BaremetalPi.SPI_IOC_MESSAGE(4)       == 0x40806b00
    @test BaremetalPi.SPI_IOC_MESSAGE(5)       == 0x40a06b00
    @test BaremetalPi.SPI_IOC_MESSAGE(6)       == 0x40c06b00
    @test BaremetalPi.SPI_IOC_MESSAGE(7)       == 0x40e06b00
    @test BaremetalPi.SPI_IOC_MESSAGE(8)       == 0x41006b00
    @test BaremetalPi.SPI_IOC_MESSAGE(9)       == 0x41206b00
    @test BaremetalPi.SPI_IOC_MESSAGE(10)      == 0x41406b00
end


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 = "XPA"
version = v"2.1.19"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/ericmandel/xpa.git", "c0452e139134d6d1677b5e6fda2ad5283c8ba4c7"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/xpa
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j$(nproc)
if [[ ${target} == *mingw* ]]; then
    # Target Windows specifically until https://github.com/ericmandel/xpa/pull/11 is merged and released
    # absolutely filthy hack edits generated Makefile manually
    sed -i 's/$(INSTALL_PROGRAM) $$i$(EXE)/$(INSTALL_PROGRAM) $$i/' Makefile
fi
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("xpamb", :xpamb),
    ExecutableProduct("xpaget", :xpaget),
    ExecutableProduct("xpainfo", :xpainfo),
    ExecutableProduct("xpans", :xpans),
    ExecutableProduct("xpaset", :xpaset),
    ExecutableProduct("xpaaccess", :xpaaccess)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

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


text:
"""
    AnovaEffect

Supertype for all ANOVA data items: `AnovaValue`, `AnovaFactor`, `AnovaResult`
"""
abstract type AnovaEffect end

Broadcast.broadcastable(a::T) where {T <: AnovaEffect} = (a,) # workaround for current behavior