text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl,md
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.0
#   kernelspec:
#     display_name: Julia 1.4.0
#     language: julia
#     name: julia-1.4
# ---

using PyPlot, NNlib,  ImageFiltering, Images, MEngProject, MEngProject.LamKernels, MEngProject.Laminart, MEngProject.Utils

σ_2 = 0.5
θ = pi/3
K = 12;

# + pycharm={"name": "#%%\n"}
for k=1:1:2K
    θ = π*(k-1)/K
    Utils.plot_rb(LamKernels.kern_d_pv(σ_2, θ))
end

# + pycharm={"name": "#%%\n"}
for k=1:1:2K
    θ = π*(k-1)/K
    Utils.plot_rb(LamKernels.kern_d_ph(σ_2, θ))
end

# + pycharm={"name": "#%%\n"}
for k=1:1:2K
    θ = π*(k-1)/K
    Utils.plot_rb(LamKernels.kern_d_mv(σ_2, θ))
end

# + pycharm={"name": "#%%\n"}
for k=1:1:2K
    θ = π*(k-1)/K
    Utils.plot_rb(LamKernels.kern_d_mh(σ_2, θ))
end

# + pycharm={"name": "#%%\n"}
for k=1:1:2K
    θ = π*(k-1)/K
    Utils.plot_rb(LamKernels.kern_d_p(σ_2, θ))
end

# + pycharm={"name": "#%%\n"}

# + pycharm={"name": "#%%\n"}
for k=1:1:2K
    θ = π*(k-1)/K
    Utils.plot_rb(LamKernels.kern_d_m(σ_2, θ))
end

# + pycharm={"name": "#%%\n"}
for k=1:1:2K
    θ = π*(k-1)/K
    Utils.plot_rb(LamKernels.kern_d(σ_2, θ))
end

# + pycharm={"name": "#%%\n"}
for k=1:1:2K
    θ = π*(k-1)/K
    Utils.plot_rb(LamKernels.kern_d(σ_2, θ))
end

# + pycharm={"name": "#%%\n"}
for k=1:1:2K
    θ = π*(k-1)/K
    Utils.plot_rb(LamKernels.kern_b(σ_2, θ))
end

# + pycharm={"name": "#%%\n"}

# + pycharm={"name": "#%%\n"}

# + pycharm={"name": "#%%\n"}

# +

Gray.(-(LamKernels.kern_d(σ_2, θ)-LamKernels.kern_d_p(σ_2, θ)))

# + pycharm={"name": "#%%\n"}
@manipulate for k=1:1:2K
    θ = π*(k-1)/K
    Gray.(3 * LamKernels.kern_d(σ_2, θ))
end

# + pycharm={"name": "#%%\n"}
@manipulate for k=1:1:2K
    θ = π*(k-1)/K
    Gray.(3 * LamKernels.kern_b(σ_2, θ))
end

# + pycharm={"name": "#%%\n"}

# + pycharm={"name": "#%%\n"}

# + pycharm={"name": "#%%\n"}

# + pycharm={"name": "#%%\n"}

# + pycharm={"name": "#%%\n"}

# + pycharm={"name": "#%%\n"}

# + pycharm={"name": "#%%\n"}

# + pycharm={"name": "#%%\n"}

# + pycharm={"name": "#%%\n"}

# -


text:
# TODO: Replace by macro
# Be a function: y = f(x, args...) that depends uppon
# a parameter vector θ, for x dependent on θ and ϕ.
# Provides a routine for computing y, dy/dθ and dy/dϕ,
# given x, dx/dθ and dx/dϕ.
# We expect as argument a function:
# Jf(y, df/dθ, df/dx, x, args...)
# that is able to evaluate the function output y and to
# compute its derivatives df/dθ, df/dx.
function sensitivity_equation(Jf::Function, jacobian_buffer)
    function f(y_extended, x_extended, args...; dydx=jacobian_buffer)
        # Expand
        x, dxdθ, dxdϕ = x_extended
        y, dydθ, dydϕ = y_extended
        # Evaluate and compute derivatives
        Jf(y, dydx, dydθ, x, args...)
        # Using the chain rule:
        # dy/dθ += dy/dx*dx/dθ
        Base.LinAlg.BLAS.gemm!('N', 'N', 1.0, dydx, dxdθ, 1.0, dydθ)
        # dy/dϕ = dy/dx*dx/dxϕ
        A_mul_B!(dydϕ, dydx, dxdϕ)
    end
    return f
end

ext: jl


text:
using AWSCore, AWSSDK, Dates, Sockets

function spot_wrangle(no_instances::Integer, spot_price::AbstractFloat, sgrp_name::String, sgrp_desc::String, skeys::String, zone::String, ami::String, instance_type::String)
    #ssh permissions
    WANip=chomp(read(pipeline(`dig +short myip.opendns.com @resolver1.opendns.com`), String))
    permissions = [
        [
            "FromPort" => 22,
            "IpProtocol" => "tcp",
            "IpRanges" => [
                [
                    "CidrIp" => WANip*"/32",
                    "Description" => "SSH access"
                ]
            ],
            "ToPort" => 22
        ]
    ]

    @info "Configuring AWS credentials..."
    aws=aws_config()

    @info "Getting security group info..."
    sgrp=Dict()

    try
        sgrp=AWSSDK.EC2.create_security_group(aws,GroupDescription=sgrp_desc,GroupName=sgrp_name)
        inret=AWSSDK.EC2.authorize_security_group_ingress(aws, GroupId=sgrp["groupId"], IpPermissions=permissions)
        outret=AWSSDK.EC2.authorize_security_group_egress(aws, GroupId=sgrp["groupId"], IpPermissions=permissions)
    catch error
        if error.code == "InvalidGroup.Duplicate"
            sgrp=AWSSDK.EC2.describe_security_groups(aws,GroupName=sgrp_name)["securityGroupInfo"]["item"]
        else
            throw(error)
        end
    end

    #Launch specs
    LaunchSpec = [
        "ImageId" => ami,
        "InstanceType" => instance_type,
        "SecurityGroupId" => [
            sgrp["groupId"]
        ],
        "Placement" => [
            "AvailabilityZone" => zone
        ],
        "KeyName" => skeys
    ]

    @info "Requesting spot instances..."

    spot_reqs=AWSSDK.EC2.request_spot_instances(aws, LaunchSpecification=LaunchSpec, InstanceCount=no_instances, SpotPrice=spot_price, Type="one-time")

    instance_public_ips=Vector()

    @info "Waiting for instances to become active..."

    if no_instances==1
        spot_req_Id=spot_reqs["spotInstanceRequestSet"]["item"]["spotInstanceRequestId"]
        active=false
        while !active
            AWSSDK.EC2.describe_spot_instance_requests(aws, SpotInstanceRequestId=spot_req_Id)["spotInstanceRequestSet"]["item"]["state"] == "active" ? active=true : sleep(10)
        end
        instance_Id=AWSSDK.EC2.describe_spot_instance_requests(aws, SpotInstanceRequestId=spot_req_Id)["spotInstanceRequestSet"]["item"]["instanceId"]
        address=AWSSDK.EC2.describe_instances(aws, InstanceId=instance_Id)["reservationSet"]["item"]["instancesSet"]["item"]["ipAddress"]
        push!(instance_public_ips,address)
    else
        for spot_req in spot_reqs["spotInstanceRequestSet"]["item"]
            spot_req_Id = spot_req["spotInstanceRequestId"]
            active=false
            while !active
                AWSSDK.EC2.describe_spot_instance_requests(aws, SpotInstanceRequestId=spot_req_Id)["spotInstanceRequestSet"]["item"]["state"] == "active" ? active=true : sleep(10)
            end
            instance_Id = AWSSDK.EC2.describe_spot_instance_requests(aws, SpotInstanceRequestId=spot_req_Id)["spotInstanceRequestSet"]["item"]["instanceId"]
            address=AWSSDK.EC2.describe_instances(aws, InstanceId=instance_Id)["reservationSet"]["item"]["instancesSet"]["item"]["ipAddress"]
            push!(instance_public_ips,address)
        end
    end

    return instance_public_ips
end

function get_cheapest_zone(instance_type::String)
    arguments=[
        "ProductDescription"=>["Linux/UNIX"]
        "InstanceType"=>[instance_type]
        "StartTime"=>Dates.now(Dates.UTC)
    ]
    history=AWSSDK.EC2.describe_spot_price_history(aws_config(), arguments)

    price_dict=Dict{String,Float64}()
    for record in history["spotPriceHistorySet"]["item"]
        price_dict[record["availabilityZone"]]=parse(Float64, record["spotPrice"])
    end

    return sort(collect(price_dict), by=x->x[2])[1]
end


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

# usage: dump_to_hdf5.jl [mps file source] [hdf5 file dest]

import HDF5
import MathProgBase
import GLPK  # Used for reading MPS files.
import GLPKMathProgInterface
using SparseArrays

"""
```
minimize objective_vector' * x + objective_constant

s.t. constraint_matrix[1:num_equalities, :] * x =
     right_hand_side[1:num_equalities]

     constraint_matrix[(num_equalities + 1):end, :] * x >=
     right_hand_side[(num_equalities + 1):end]

     variable_lower_bound <= x <= variable_upper_bound
```
"""
mutable struct LinearProgrammingProblem
  variable_lower_bound::Vector{Float64}
  variable_upper_bound::Vector{Float64}
  objective_vector::Vector{Float64}
  objective_constant::Float64
  constraint_matrix::SparseMatrixCSC{Float64,Int64}
  right_hand_side::Vector{Float64}
  num_equalities::Int64
end

mutable struct TwoSidedLpProblem
  variable_lower_bound::Vector{Float64}
  variable_upper_bound::Vector{Float64}
  constraint_lower_bound::Vector{Float64}
  constraint_upper_bound::Vector{Float64}
  constraint_matrix::SparseArrays.SparseMatrixCSC{Float64,Int64}
  objective_vector::Vector{Float64}
end

function transform_to_standard_form(lp::TwoSidedLpProblem)
  is_equality_row = lp.constraint_lower_bound .== lp.constraint_upper_bound
  is_geq_row = .!is_equality_row .& isfinite.(lp.constraint_lower_bound)
  is_leq_row = .!is_equality_row .& isfinite.(lp.constraint_upper_bound)

  @assert !any(is_geq_row .& is_leq_row)

  num_equalities = sum(is_equality_row)

  if num_equalities +
     sum(is_geq_row) +
     sum(is_leq_row) != length(lp.constraint_lower_bound)
    error("Not all constraints have finite bounds on at least one side.")
  end

  # Flip the signs of the leq rows in place.
  for idx in 1:SparseArrays.nnz(lp.constraint_matrix)
    if is_leq_row[lp.constraint_matrix.rowval[idx]]
      lp.constraint_matrix.nzval[idx] *= -1
    end
  end

  new_row_to_old = [findall(is_equality_row); findall(.!is_equality_row)]
  if new_row_to_old != 1:size(lp.constraint_matrix, 1)
    row_permute_in_place(lp.constraint_matrix, invperm(new_row_to_old))
  end

  right_hand_side = copy(lp.constraint_lower_bound)
  right_hand_side[is_leq_row] .= .-lp.constraint_upper_bound[is_leq_row]
  permute!(right_hand_side, new_row_to_old)

  return LinearProgrammingProblem(
    lp.variable_lower_bound,
    lp.variable_upper_bound,
    lp.objective_vector,
    0.0,
    lp.constraint_matrix,
    right_hand_side,
    num_equalities,
  )
end

function read_mps_to_standard_form(filename::String)
  if !endswith(filename, ".mps.gz") && !endswith(filename, ".mps")
    error("Invalid name: `filename` needs to end with `.mps.gz` or `.mps`.")
  end

  model = MathProgBase.LinearQuadraticModel(GLPKMathProgInterface.GLPKSolverMIP(
    msg_lev = GLPK.MSG_OFF,
  ))
  MathProgBase.loadproblem!(model, filename)

  return transform_to_standard_form(TwoSidedLpProblem(
    MathProgBase.getvarLB(model),
    MathProgBase.getvarUB(model),
    MathProgBase.getconstrLB(model),
    MathProgBase.getconstrUB(model),
    MathProgBase.getconstrmatrix(model),
    MathProgBase.getobj(model),
  ))
end

lp = read_mps_to_standard_form(ARGS[1])

HDF5.h5open(ARGS[2], "w") do file
  write(file, "variable_lower_bound", lp.variable_lower_bound)
  write(file, "variable_upper_bound", lp.variable_upper_bound)
  write(file, "right_hand_side", lp.right_hand_side)
  write(file, "objective_vector", lp.objective_vector)
  write(file, "objective_constant", lp.objective_constant)
  write(file, "num_equalities", lp.num_equalities)
  # The matrix in zero-based CSC format.
  write(file, "constraint_matrix_data", lp.constraint_matrix.nzval)
  write(file, "constraint_matrix_indices", lp.constraint_matrix.rowval .- 1)
  write(file, "constraint_matrix_indptr", lp.constraint_matrix.colptr .- 1)
end


text:
using GOEnrichment
using Test

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


text:
#Visualization

#Colormap
using Color
cm = Uint32[Color.convert(Color.RGB24,c) for c in flipud(Color.colormap("RdBu"))]
#cm = Uint32[Color.convert(Color.RGB24,c) for c in Color.colormap("RdBu")]

#Visualize common matrix
function visualize_M(data::AbstractMatrix,
                     x, y;
                     g=1,
                     title="",
                     dmin=0.0,
                     dmax=:find)

    ny, nx = size(data)
    xs = g +1
    xe = nx - g + 1
    ye = ny - g + 1

    xmin = x[xs]
    xmax = x[xe]
    ymin =  y[xs]
    ymax =  y[ye]

    hdata = data[xs:ye, xs:xe]
    p = FramedPlot()

    #look or set min and max values
    if dmin == :find
        hmin = minimum(hdata)
    elseif isa(dmin, Real)
        hmin = dmin
    end

    if dmax == :find
        hmax = maximum(hdata)
    elseif isa(dmax, Real)
        hmax = dmax
    end

    clims = (minimum(hdata), maximum(hdata))

    #make image
    img = Winston.data2rgb(hdata, clims, cm)
    add(p, Image((xmin, xmax), (ymin, ymax), img;))
    setattr(p, xrange=(xmin, xmax))
    setattr(p, yrange=(ymin, ymax))
    setattr(p, title=title)

    return p
end

function visualize(hyd::data2d)

    #density
    p1 = visualize_M(hyd.rho,
                     hyd.x, hyd.y,
                     g=3,
                     title="density \\rho",
                     dmin = :find,
                     dmax = :find)

    #pressure
    p2 = visualize_M(hyd.press,
                     hyd.x, hyd.y,
                     g=3,
                     title="pressure p",
                     dmin = 0.0,
                     dmax = :find)

    #velocity
    hdata = sqrt(hyd.velx.^2.0 .+ hyd.vely.^2.0)
    #hdata = hyd.velx[xs:ye, xs:xe]

    p3 = visualize_M(hdata,
                     hyd.x, hyd.y,
                     g=3,
                     title="velocity (u^2 + v^2)^{1/2}",
                     dmin = :find,
                     dmax = :find)

    #internal energy
    p4 = visualize_M(hyd.press,
                     hyd.x, hyd.y,
                     g=3,
                     title="internal energy ε",
                     dmin = 0.0,
                     dmax = :find)

    ######
    #additional window for 1-dim figs

    #diagonal slice
    xy = Float64[sqrt(hyd.x[i]^2+hyd.y[i]^2.0) for i = 1:hyd.nx]
    rhoxy = Float64[hyd.rho[hyd.nx-i+1,i] for i = hyd.nx:-1:1]
    velxy = Float64[sqrt(hyd.velx[hyd.nx-i+1,i]^2.0 + hyd.velx[hyd.nx-i+1,i]^2.0) for i = hyd.nx:-1:1]

    velxxy = Float64[hyd.velx[hyd.nx-i+1,i] for i = hyd.nx:-1:1]
    velyxy = Float64[hyd.vely[hyd.nx-i+1,i] for i = hyd.nx:-1:1]

    pressxy = Float64[hyd.press[hyd.nx-i+1,i] for i = hyd.nx:-1:1]
    epsxy = Float64[hyd.eps[hyd.nx-i+1,i] for i = hyd.nx:-1:1]

    p11 =  plot(xy, rhoxy, "r",title="diagonal slice")
    p11 = oplot(xy, abs(velxxy), "b")
    p11 = oplot(xy, abs(velyxy), "c")
    p11 = oplot(xy, pressxy, "g")
    #p11 = oplot(hyd.y, epsxy, "k")

    #middle slice
    mid = int(hyd.nx/2)
    p12 = plot(hyd.y, hyd.rho[:,mid], "r",title="middle slice")
    p12 = oplot(hyd.y, hyd.vely[:,mid], "b")
    p12 = oplot(hyd.y, hyd.press[:,mid], "g")
    #p12 = oplot(hyd.y, hyd.eps[:,mid], "k")

    t = Table(2,3)
    t[1,1] = p1
    t[1,2] = p2
    t[2,1] = p3
    t[2,2] = p4

    t[1,3] = p11
    t[2,3] = p12

    display(t)

end

max_stars_repo_name: natj/hydro

max_issues_repo_name: natj/hydro

max_forks_repo_name: natj/hydro


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5]);F = [5184*x1^4-12528*x1^3*x3+100*x1^2*x2^2+1640*x1^2*x2*x3+5348*x1^2*x2*x5+14293*x1^2*x3^2-11644*x1^2*x3*x5-12960*x1^2*x4*x5+11233*x1^2*x5^2-320*x1*x2^3-2624*x1*x2^2*x3+2272*x1*x2^2*x5-1660*x1*x2*x3*x4-8178*x1*x2*x3*x5-180*x1*x2*x4*x5-13612*x1*x3^2*x4+25970*x1*x3*x4*x5-7482*x1*x3*x5^2+1278*x1*x4*x5^2+256*x2^4+2656*x2^2*x3*x4+288*x2^2*x4*x5+2209*x2^2*x5^2-8460*x2*x4*x5^2+4042*x2*x5^3+6889*x3^2*x4^2+1494*x3*x4^2*x5+8181*x4^2*x5^2-7740*x4*x5^3+1849*x5^4-13824*x1^3+16704*x1^2*x3-9024*x1*x2*x5+17280*x1*x4*x5-8256*x1*x5^2+9216*x1^2, 200*x1^2*x2+1640*x1^2*x3+5348*x1^2*x5-960*x1*x2^2-5248*x1*x2*x3+4544*x1*x2*x5-1660*x1*x3*x4-8178*x1*x3*x5-180*x1*x4*x5+1024*x2^3+5312*x2*x3*x4+576*x2*x4*x5+4418*x2*x5^2-8460*x4*x5^2+4042*x5^3-9024*x1*x5, -12528*x1^3+1640*x1^2*x2+28586*x1^2*x3-11644*x1^2*x5-2624*x1*x2^2-1660*x1*x2*x4-8178*x1*x2*x5-27224*x1*x3*x4+25970*x1*x4*x5-7482*x1*x5^2+2656*x2^2*x4+13778*x3*x4^2+1494*x4^2*x5+16704*x1^2, -12960*x1^2*x5-1660*x1*x2*x3-180*x1*x2*x5-13612*x1*x3^2+25970*x1*x3*x5+1278*x1*x5^2+2656*x2^2*x3+288*x2^2*x5-8460*x2*x5^2+13778*x3^2*x4+2988*x3*x4*x5+16362*x4*x5^2-7740*x5^3+17280*x1*x5, 5348*x1^2*x2-11644*x1^2*x3-12960*x1^2*x4+22466*x1^2*x5+2272*x1*x2^2-8178*x1*x2*x3-180*x1*x2*x4+25970*x1*x3*x4-14964*x1*x3*x5+2556*x1*x4*x5+288*x2^2*x4+4418*x2^2*x5-16920*x2*x4*x5+12126*x2*x5^2+1494*x3*x4^2+16362*x4^2*x5-23220*x4*x5^2+7396*x5^3-9024*x1*x2+17280*x1*x4-16512*x1*x5];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
module test_http_multipart_basic

using Test
using HTTP: Multipart

part1 = Multipart("filename", IOBuffer())
@test eof(part1)

end # module test_http_multipart_basic

module test_http_multipart_scanner

using Test
using HTTP: Multipart

mutable struct FormScanner
    data::Vector{UInt8}
    boundary::String
    FormScanner(data::Vector{UInt8}, boundary::String) = new(data, boundary)
end

const CR = 0x0d
const LF = 0x0a
const ContentDisposition = "Content-Disposition: form-data;"

function rstripcr(vec::Vector{UInt8})::Vector{UInt8}
    if !isempty(vec) && vec[end] == CR
        vec[1:end-1]
    else
        vec
    end
end

function rstripcrlf(vec::Vector{UInt8})::Vector{UInt8}
    if !isempty(vec) && vec[end] == LF
        rstripcr(vec[1:end-1])
    else
        vec
    end
end

function content_disposition_fields(io::IOBuffer)::Dict{String,String}
    fields = Dict{String,String}()
    str = String(readuntil(io, LF))
    h = rstripcr(readuntil(io, LF))
    if isempty(h)
    else
        (k, v) = split(String(h), ": ")
        fields[k] = v
        readuntil(io, LF)
    end
    pat = r" (?P<key>[^\"]*)=\"(?P<val>[^\"]*)\""
    for m in eachmatch(pat, str)
        fields[m[:key]] = m[:val]
    end
    fields
end

function scan(s::FormScanner)::Vector{Pair{String,Any}}
    ps = Vector{Pair{String, Any}}()
    io = IOBuffer(s.data, read=true, truncate=false)
    while !eof(io)
        chunk = readuntil(io, s.boundary)
        if !isempty(chunk)
            chunkio = IOBuffer(chunk)
            skipchars(==(Char(CR)), chunkio)
            skipchars(==(Char(LF)), chunkio)
            readuntil(chunkio, ContentDisposition)
            contentdisposition = rstripcrlf(read(chunkio))
            if !isempty(contentdisposition)
                contentio = IOBuffer(contentdisposition)
                fields = content_disposition_fields(contentio)
                name = fields["name"]
                if haskey(fields, "filename")
                    filename = fields["filename"]
                    contenttype = get(fields, "Content-Type", "")
                    # contenttransferencoding
                    multipart = Multipart(filename, contentio, contenttype)
                    push!(ps, Pair(name, multipart))
                else
                    value = read(contentio)
                    push!(ps, Pair(name, String(value)))
                end
            end
        end
    end
    ps
end

# https://discourse.julialang.org/t/http-multipart-form-data-processing-by-server/24076/3
reqbody = """
-----------------------------182023285717490760841965583652
Content-Disposition: form-data; name="image"; filename="file1.jpg"
Content-Type: image/jpeg

......JFIF.............C..........
-----------------------------182023285717490760841965583652
Content-Disposition: form-data; name="num"

2
-----------------------------182023285717490760841965583652--
"""
boundary = "---------------------------182023285717490760841965583652"
scanner = FormScanner(Vector{UInt8}(reqbody), string("--", boundary))
body_params = scan(scanner)

@test body_params[1][1] == "image"
multipart = body_params[1][2]
@test multipart isa Multipart
@test String(read(multipart.data)) == "......JFIF.............C.........."

@test body_params[2] == ("num" => "2")

end # module test_http_multipart_scanner


text:
type BitFlipMutationDefinition <: MutationDefinition
  stage::AbstractString
  rate::Float

  BitFlipMutationDefinition() = new("", 0.01)
  BitFlipMutationDefinition(rate::Float) = new("", rate)
  BitFlipMutationDefinition(stage::AbstractString) = new(stage, 0.01)
  BitFlipMutationDefinition(stage::AbstractString, rate::Float) =
    new(stage, rate)
end

type BitFlipMutation <: Mutation
  stage::AbstractString
  rate::Float
  representation::representation.BitVectorRepresentation

  BitFlipMutation(stage::AbstractString, rate::Float, rep::representation.BitVectorRepresentation) =
    new(stage, rate, rep)
end

function compose!(d::BitFlipMutationDefinition, sp::Species)
  d.stage = d.stage == "" ? genotype(sp).label : d.stage
  BitFlipMutation(d.stage, d.rate, sp.stages[d.stage].representation)
end

"""
Performs bit-flip mutation on a fixed or variable length chromosome of binary
digits, by flipping 1s to 0s and 0s to 1s at each point within the chromosome
with a given probability, equal to the mutation rate.

**Parameters:**

* `stage::AbstractString`, the name of the developmental stage that this
  operator should be applied to. Defaults to the genotype if no stage is
  specified.
* `rate::Float`, the probability of a bit flip at any given index.
  Defaults to 0.01 if no rate is provided.
"""
bit_flip() = bit_flip(0.01)
bit_flip(rate::Float) = BitFlipMutationDefinition(rate)
bit_flip(stage::AbstractString) = BitFlipMutationDefinition(stage)
bit_flip(stage::AbstractString, rate::Float) =
  BitFlipMutationDefinition(stage, rate)

function mutate!(o::BitFlipMutation, input::Vector{Int})
  for i in 1:length(input)
    if rand() <= o.rate
      input[i] = input[i] == 1 ? 0 : 1
    end
  end
  return input
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

name = "SPIRV_Headers"
version = v"1.5.3"

# Collection of sources required to build this package
sources = [
    GitSource("https://github.com/KhronosGroup/SPIRV-Headers.git",
              "c0df742ec0b8178ad58c68cff3437ad4b6a06e26"),
]

# Bash recipe for building across all platforms
script = raw"""
cd SPIRV-Headers
install_license LICENSE

CMAKE_FLAGS=()

# Install things into $prefix
CMAKE_FLAGS+=(-DCMAKE_INSTALL_PREFIX=${prefix})

cmake -B build -S . -GNinja ${CMAKE_FLAGS[@]}
ninja -C build -j ${nproc} install
"""

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

# The products that we will ensure are always built
products = [
    FileProduct("include/spirv/spir-v.xml", :spirv_xml)
]

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

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
global srand_core=zero(UInt32)
function srand(x::Integer)
	global srand_core=UInt32(x)
end
function crand()
	global srand_core=(srand_core*1103515245+12345)
	return (srand_core>>16)&32767
end
"[l,r]"
function crandom(l::Int,r::Int)
	return l+crand()%(r-l+1)
end
"[l,r]中与l,r同奇偶"
function crandom_even(l::Int,r::Int)
	return l+(crand()%((r-l)>>1+1))<<1
end


text:
using Printf
using LinearAlgebra

include("../common/PWGrid_v03.jl")
include("../common/gen_lattice_pwscf.jl")

function test_main()
    LatVecs = gen_lattice_sc(8.0)

    pw = PWGrid( 15.0, LatVecs )
    Ns = pw.Ns
    @printf("Sampling points: (%d,%d,%d)\n", Ns[1], Ns[2], Ns[3])

    actual = prod(pw.Ns)/pw.gvecw.Ngwx
    theor = 1/(4*pi*0.25^3/3)
    @printf("Actual, theor: %10.5f %10.5f\n", actual, theor)

end

@time test_main()


text:
module Util
export  dyadicdetailindex,
        dyadicdetailrange,
        dyadicscalingrange,
        dyadicdetailn,
        ndyadicscales,
        maxdyadiclevel,
        tl2dyadiclevel,
        dyadiclevel2tl,
        # non-dyadic
        detailindex,
        detailrange,
        detailn,
        maxtransformlevels,
        #
        mirror,
        upsample,
        downsample,
        iscube,
        isdyadic,
        sufficientpoweroftwo,
        wcount,
        stridedcopy!,
        isvalidtree,
        maketree,
        makewavelet,
        testfunction

# WAVELET INDEXING AND SIZES

# Non-dyadic
# detail coef at level l location i (i starting at 1) -> vector index
detailindex(arraysize::Integer, l::Integer, i::Integer) = round(Int, arraysize/2^l+i)
detailindex(x::AbstractArray, l::Integer, i::Integer) = detailindex(size(x,1), l ,i)
# the range of detail coefs at level l
detailrange(arraysize::Integer, l::Integer) = round(Int, (arraysize/2^l+1)) : round(Int, arraysize/2^(l-1))
detailrange(x::AbstractArray, l::Integer) = detailrange(size(x,1), l)
# number of detail coefs at level l
detailn(arraysize::Integer, l::Integer) = round(Int, arraysize/2^l)
detailn(x::AbstractArray, l::Integer) = detailn(size(x,1), l)
# max levels to transform
maxtransformlevels(x::AbstractArray) = maxtransformlevels(minimum(size(x)))
function maxtransformlevels(arraysize::Integer)
    arraysize > 1 || return 0
    tl = 0
    while sufficientpoweroftwo(arraysize, tl)
        tl += 1
    end
    return tl - 1
end

# Dyadic
# detail coef at level j location i (i starting at 1) -> vector index
dyadicdetailindex(j::Integer,i::Integer) = 2^j+i
# the range of detail coefs at level j
dyadicdetailrange(j::Integer) = (2^j+1):(2^(j+1))
# the range of scaling coefs at level j
dyadicscalingrange(j::Integer) = 1:2^j
# number of detail coefs at level j
dyadicdetailn(j::Integer) = 2^j
# number of scales of dyadic length signal (n=2^J)
ndyadicscales(n::Integer) = round(Int, log2(n))
ndyadicscales(x::AbstractArray) = ndyadicscales(size(x,1))
# the largest detail level
maxdyadiclevel(n::Integer) = ndyadicscales(n)-1
maxdyadiclevel(x::AbstractArray) = maxdyadiclevel(size(x,1))
# convert number of transformed levels L to minimum level j
tl2dyadiclevel(n::Integer, L::Integer) = ndyadicscales(n)-L
tl2dyadiclevel(x::AbstractVector, L::Integer) = tl2dyadiclevel(length(x), L)
# convert maximum level j to number of transformed levels L
dyadiclevel2tl(arg...) = tl2dyadiclevel(arg...)

# UTILITY FUNCTIONS

# are all dimensions equal length?
function iscube(x::AbstractArray)
    for i = 1:ndims(x)
        size(x,1) != size(x,i) && return false
    end
    return true
end
# are all dimensions dyadic length?
function isdyadic(x::AbstractArray)
    for i = 1:ndims(x)
        isdyadic(size(x,i)) || return false
    end
    return true
end
isdyadic(n::Integer) = (n == 2^(ndyadicscales(n)))

# To perform a level L transform, the size of the signal in each dimension
# must have 2^L as a factor.
function sufficientpoweroftwo(x::AbstractArray, L::Integer)
    for i = 1:ndims(x)
        sufficientpoweroftwo(size(x,i), L) || return false
    end
    return true
end
sufficientpoweroftwo(n::Integer, L::Integer) = (n%(2^L) == 0)

# mirror of filter
mirror(f::AbstractVector{<:Number}) = f .* (-1).^(0:length(f)-1)
# upsample
function upsample(x::AbstractVector, sw::Int=0)
    @assert sw==0 || sw==1
    n = length(x)
    y = zeros(eltype(x), n<<1)
    sw -= 1

    for i = 1:n
        @inbounds y[i<<1 + sw] = x[i]
    end
    return y
end
# downsample
function downsample(x::AbstractVector, sw::Int=0)
    @assert sw==0 || sw==1
    n = length(x)
    @assert n%2 == 0
    y = zeros(eltype(x), n>>1)
    sw -= 1

    for i in eachindex(y)
        @inbounds y[i] = x[i<<1 + sw]
    end
    return y
end

# count coefficients above threshold t (>=), excluding coefficients in levels < level
# where level -1 is the x[1] coefficient
function wcount(x::AbstractVector, t::Real=0; level::Int=-1)
    @assert level >= -1
    c = 0
    si = 1
    level >= 0 && (si = 1 + 2^level)
    @inbounds for i = si:length(x)
        if abs(x[i]) >= t
            c += 1
        end
    end
    return c
end
# count coefficients above threshold t (>=)
function wcount(x::AbstractArray, t::Real=0)
    c = 0
    @inbounds for i in eachindex(x)
        if abs(x[i]) >= t
            c += 1
        end
    end
    return c
end

# inplace circular shift of vector by shift, such that anew[i]=aold[i-shift] (mod length(a))
function circshift!(a::AbstractVector, shift::Integer)
    atype = typeof(a)
    s = length(a)
    shift = mod(shift,s)
    shift == 0 && return a
    shift = ifelse(s>>1 < shift, shift-s, shift)  # shift a smaller distance if possible
    if shift < 0
        tmp = a[1:-shift]
        for i = 1:s+shift
            @inbounds a[i] = a[i-shift]
        end
        a[s+1+shift:s] = tmp
    else
        tmp = a[s+1-shift:s]
        for i = s:-1:1+shift
            @inbounds a[i] = a[i-shift]
        end
        a[1:shift] = tmp
    end
    return a::atype
end
# out of place circular shift of vector by shift, such that b[i]=a[i-shift] (mod length(a))
function circshift!(b::AbstractVector, a::AbstractVector, shift::Integer)
    @assert length(a) == length(b)
    atype = typeof(a)
    s = length(a)
    shift = mod(shift,s)
    shift == 0 && return copy!(b,a)
    shift = ifelse(s>>1 < shift, shift-s, shift)  # shift a smaller distance if possible
    if shift < 0
        for i = 1:s+shift
            @inbounds b[i] = a[i-shift]
        end
        sh = s+shift
        for i = s+shift+1:s
            @inbounds b[i] = a[i-sh]
        end
    else
        for i = s:-1:1+shift
            @inbounds b[i] = a[i-shift]
        end
        sh = -s+shift
        for i = shift:-1:1
            @inbounds b[i] = a[i-sh]
        end
    end
    return b::atype
end

# put odd elements into first half, even into second half
function split!(a::AbstractVector{T}) where T<:Number
    n = length(a)
    nt = n>>2 + (n>>1)%2
    tmp = Vector{T}(nt)
    split!(a, n, tmp)
    return a
end

# split only the range 1:n
function split!(a::AbstractVector{T}, n::Integer, tmp::Vector{T}) where T<:Number
    @assert n <= length(a)
    @assert n%2 == 0
    n == 2 && return a
    nt = n>>2 + (n>>1)%2
    @assert nt <= length(tmp)

    for i=1:nt # store evens
        @inbounds tmp[i] = a[i<<1]
    end
    for i=1:n>>1  # odds to first part
        @inbounds a[i] = a[(i-1)<<1 + 1]
    end
    for i=0:nt - 1  # evens to end
        @inbounds a[n-i] = a[n - i<<1]
    end
    copy!(a,n>>1+1,tmp,1,nt)
    return a
end

# out of place split from a to b, only the range 1:n
function split!(b::AbstractVector{T}, a::AbstractVector{T}, n::Integer) where T<:Number
    @assert n <= length(a) && n <= length(b)
    @assert n%2 == 0
    if n == 2
        b[1] = a[1]
        b[2] = a[2]
        return b
    end

    h = n>>1
    for i=1:h  # odds to b
        @inbounds b[i] = a[(i-1)<<1 + 1]
    end
    for i=h+1:n  # evens to b
        @inbounds b[i] = a[(i - h)<<1]
    end

    return b
end
# out of place split from a to b, only the range a[ia:inca:ia+(n-1)*inca] to b[1:n]
function split!(b::AbstractVector{T}, a::AbstractArray{T}, ia::Integer, inca::Integer, n::Integer) where T<:Number
    @assert ia+(n-1)*inca <= length(a) && n <= length(b)
    @assert n%2 == 0
    if n == 2
        b[1] = a[ia]
        b[2] = a[ia + inca]
        return b
    end

    h = n>>1
    inca2 = inca<<1
    for i=1:h  # odds to b
        @inbounds b[i] = a[ia + (i-1)*inca2]
    end
    iainca = ia + inca
    hp1 = h + 1
    for i=h+1:n  # evens to b
        @inbounds b[i] = a[iainca + (i - hp1)*inca2]
    end

    return b
end

# inverse the operation of split!
function merge!(a::AbstractVector{T}) where T<:Number
    n = length(a)
    nt = n>>2 + (n>>1)%2
    tmp = Vector{T}(nt)
    merge!(a, n, tmp)
    return a
end

# merge only the range 1:n
function merge!(a::AbstractVector{T}, n::Integer, tmp::Vector{T}) where T<:Number
    @assert n <= length(a)
    @assert n%2 == 0
    n == 2 && return a
    nt = n>>2 + (n>>1)%2
    @assert nt <= length(tmp)

    copy!(tmp,1,a,n>>1+1,nt)
    for i=nt-1:-1:0  # evens from end
        @inbounds a[n - 2*i] = a[n-i]
    end
    for i=n>>1:-1:1  # odds from first part
        @inbounds a[(i-1)<<1 + 1] = a[i]
    end
    for i=nt:-1:1 # retrieve evens
        @inbounds a[i<<1] = tmp[i]
    end
    return a
end

# out of place merge from a to b, only the range 1:n
function merge!(b::AbstractVector{T}, a::AbstractVector{T}, n::Integer) where T<:Number
    @assert n <= length(a) && n <= length(b)
    @assert n%2 == 0
    if n == 2
        b[1] = a[1]
        b[2] = a[2]
        return b
    end

    h = n>>1
    for i=1:h  # odds to b
        @inbounds b[(i-1)<<1 + 1] = a[i]
    end
    for i=h+1:n  # evens to b
        @inbounds b[(i - h)<<1] = a[i]
    end

    return b
end
# out of place merge from a to b, only the range a[1:n] to b[ib:incb:ib+(n-1)*incb]
function merge!(b::AbstractArray{T}, ib::Integer, incb::Integer, a::AbstractVector{T}, n::Integer) where T<:Number
    @assert n <= length(a) && ib+(n-1)*incb <= length(b)
    @assert n%2 == 0
    if n == 2
        b[ib] = a[1]
        b[ib + incb] = a[2]
        return b
    end

    h = n>>1
    incb2 = incb<<1
    for i=1:h  # odds to b
        @inbounds b[ib + (i-1)*incb2] = a[i]
    end
    ibincb = ib + incb
    hp1 = h + 1
    for i=h+1:n  # evens to b
        @inbounds b[ibincb + (i - hp1)*incb2] = a[i]
    end

    return b
end

function stridedcopy!(b::AbstractVector{T}, a::AbstractArray{T}, ia::Integer, inca::Integer, n::Integer) where T<:Number
    @assert ia+(n-1)*inca <= length(a) && n <= length(b)

    @inbounds for i = 1:n
        b[i] = a[ia + (i-1)*inca]
    end
    return b
end
function stridedcopy!(b::AbstractArray{T}, ib::Integer, incb::Integer, a::AbstractVector{T}, n::Integer) where T<:Number
    @assert ib+(n-1)*incb <= length(b) && n <= length(a)

    @inbounds for i = 1:n
        b[ib + (i-1)*incb] = a[i]
    end
    return b
end

# wavelet packet transforms WPT
# valid if 0 nodes have 0 children and length+1 is dyadic
# and has a nodes corresponding every transform level
function isvalidtree(x::AbstractVector, b::BitVector)
    ns = maxtransformlevels(x)
    nb = length(b)
    nb == 2^(ns)-1 || return false
    @assert (2^(ns-1)-1)<<1+1 <= nb

    for i in 1:2^(ns-1)-1
        @inbounds if !b[i] && (b[i<<1] || b[i<<1+1])
            return false
        end
    end
    return true
end
# return a tree (BitVector)
# s=:full, all nodes for first L levels equal 1, others 0
# s=:dwt, nodes corresponding to a dwt for first L levels equal 1, others 0
maketree(x::Vector, s::Symbol=:full) = maketree(length(x), maxtransformlevels(x), s)
function maketree(n::Int, L::Int, s::Symbol=:full)
    ns = maxtransformlevels(n)
    nb = 2^(ns)-1
    @assert 0 <= L <= ns
    @assert (2^(ns-1)-1)<<1+1 <= nb

    b = BitArray(nb)
    fill!(b, false)

    t = true
    if s == :full
        for i in 1:2^(L)-1
            @inbounds b[i] = t
        end
    elseif s == :dwt
        for i in 1:L
            @inbounds b[2^(i-1)] = t
        end
    else
        throw(ArgumentError("uknown symbol"))
    end
    return b
end

# return scaling and wavelet functions and location vector, made from filter h
# iterated with a cascade algorithm with N steps
makewavelet(h, arg...) = makewavelet(h.qmf, arg...)
function makewavelet(h::AbstractVector, N::Integer=8)
    @assert N>=0
    sc = norm(h)
    h = h*sqrt(2)/sc
    phi = copy(h)
    psi = mirror(reverse(h))

    for i=1:N
        phi = conv(upsample(phi), h)
        psi = conv(upsample(psi), h)
    end
    phi = phi[1:end-2^(N)+1]
    psi = psi[1:end-2^(N)+1]
    return scale!(phi,sc/sqrt(2)), scale!(psi,sc/sqrt(2)), linspace(0,length(h)-1,length(psi))
end

# return a vector of test function values on [0,1), see
#Donoho, D.L.; I.M. Johnstone (1994), "Ideal spatial adaptation by wavelet shrinkage," Biometrika, vol. 81, pp. 425–455.
function testfunction(n::Int, ft::AbstractString)
    @assert n >= 1
    f = Vector{Float64}(n)
    range = 0:1/n:1-eps()
    i = 1
    if ft=="Blocks"
        tj = [0.1,0.13,0.15,0.23,0.25,0.4,0.44,0.65,0.76,0.78,0.81]
        hj = [4,-5,3,-4,5,-4.2,2.1,4.3,-3.1,2.1,-4.2]
        for t in range
            f[i] = 0
            for j = 1:11
                f[i] += hj[j]*(1+sign(t-tj[j]))/2
            end
            i += 1
        end
    elseif ft=="Bumps"
        tj = [0.1,0.13,0.15,0.23,0.25,0.4,0.44,0.65,0.76,0.78,0.81]
        hj = [4,5,3,4,5,4.2,2.1,4.3,3.1,5.1,4.2]
        wj = [0.005,0.005,0.006,0.01,0.01,0.03,0.01,0.01,0.005,0.008,0.005]
        for t in range
            f[i] = 0
            for j = 1:11
                f[i] += hj[j]/(1+abs((t-tj[j])/wj[j]))^4
            end
            i += 1
        end
    elseif ft=="HeaviSine"
        for t in range
            f[i] = 4*sin(4*pi*t) - sign(t-0.3) - sign(0.72-t)
            i += 1
        end
    elseif ft=="Doppler"
        for t in range
            f[i] = sqrt(t*(1-t))*sin(2*pi*1.05/(t+0.05))
            i += 1
        end
    else
        throw(ArgumentError("unknown test function"))
    end
    return f
end

end # module


text:
using QuickCheck2
using Base.Test

@testset "EXAMPLES:" begin

## Examples for @test_forall macro. Returns 4 Passes; 1 Fail; 1 Error
@testset "@test_forall macro" begin
	@test_forall x in -1:1, x*(x-1)*(x+1) == 0 # test passed
	@test_forall x in -1:1, x*(x-2)*(x+1) == 0 # test failed
	@test_forall i in 1:2, x in -1:3, y in x+1:4, (x+i)^2<(y+i)^2 # test passed
	@test_forall x in [0,1,2], y in x:4, (y+4>2*x)==true # test passed
	@test_forall y in x:4, x in 0:2, (y+4>2*x)==true # bad design: throws UndefVarError: x not defined, because iteration over y comes earlier and refers to value of x
	@test_forall x in ["a","b"], y in ["z","w"], x*y in Set(["az","aw","bz","bw"]) # test passed
end

## Examples for @test_formany macro. Returns 3 Passes; 2 Fails; 1 Error
@testset "@test_formany macro" begin
	@test_formany 100>x>10, x+10<y<<1000, y-5<z<Inf,z>x+5 # test passed
	@test_formany 100>x>10, x+10<y<<1000, y-5<z<Inf, z>x+6 # test failed
	@test_formany ntests = 1000 Inf>x>0,Inf>y>0, x<y-->log(x)<log(y) # test passed
    @test_formany ntests = 100 maxtests = 100 0<x<10, 0<y<100, x^2>y^2-->x>y #will usually fail
    @test_formany ntests = 100 maxtests = 100000 0<x<10, 0<y<100, x^2>y^2-->x>y #will usually pass
	@test_formany -Inf<x<Inf, -Inf<y<Inf x<y-->log(x)<log(y) 1000 # test throws DomainError
	@test_formany 0<a<1,-100<x<100,-100<y<100,(a*x+(1-a)*y)^3<a*x^3+(1-a)*y^3 # test failed: x^3 is not convex on [-100,100]
	@test_formany 0<a<1,0<x<100,0<y<100, (a*x+(1-a)*y)^3<=a*x^3+(1-a)*y^3 # test passed
end

## Examples for @test_exists macro. Returns 5 Passes; 1 Fail
@testset "@test_exists macro" begin
	@test_exists ntests = 1000 -10<x<10, x^2>99  # test passed
	@test_exists ntests = 100000 0<x<10000, x>9999 # test passed
	@test_exists ntests = 100 0<x<10000, x>9999 # fails almost surely
	@test_exists 0<x<2*pi,x<y<2*pi,y<z<2*pi, sin(x)<sin(y)<sin(z) # test passes
	@test_exists ntests = 10000 0<x<2*pi,0<y<2*pi,0<z<2*pi, x<y<z --> sin(x)<sin(y)<sin(z) # test passes almost surely
	@test_exists 0<x1<2,x1<x2<2, x1<y<x2, (x1-1)^2*(x1+1)^2<(y-1)^2*(y+1)^2<(x2-1)^2*(x2+1)^2 # test passed
end

end # EXAMPLES


text:
"""
    selectCrossPoint(genotype::STGPGenotype, rng::Random.AbstractRNG;
                     fixedDepth::Integer = -1, minDepth::Integer = -1,
                     maxDepth::Integer = -1,
                     allowedType::Union{DataType, Union} = Any,
                     forbidden::Array{Integer} = Array{Integer}(undef, 0))

Selects a point in the genotype where the crossover is going to be performed
by taking some considerations into account, such as the maximum depth
and the type allowed.

# Arguments
- `genotype::CGPGenotype`: genotype from which the cross point is going to be
    selected.
- `rng::Random.AbstractRNG`: random number generator for random number consistency
    along an experiment.

# Keyword Arguments
- `fixedDepth::Integer = -1`: specifies the exact depth at where the crosspoint
    must be. If this argument is set to -1, it will be ignored.
- `minDepth::Integer = -1`: specifies the minimum depth at where the crosspoint
    can be. If this argument is set to -1, it will be ignored.
- `maxDepth::Integer = -1`: specifies the maximum depth at where the crosspoint
    can be. If this argument is set to -1, it will be ignored.
- `allowedType::Union{DataType, Union} = Any`: specifies the allowed type that the
    crosspoint must have.
- `forbidden::Array{Integer} = Array{Integer}(undef, 0)`: array of the indexes
    of the nodes that are forbidden to be crosspoints because they have already
    been proved not to have a compatible crosspoint in the other parent.
"""
function selectCrossPoint(genotype::STGPGenotype, rng::Random.AbstractRNG;
                          fixedDepth::Integer = -1, minDepth::Integer = -1,
                          maxDepth::Integer = -1,
                          allowedType::Union{DataType, Union} = Any,
                          forbidden::Array{Integer} = Array{Integer}(undef, 0))

    if fixedDepth != -1 && minDepth != -1
        minDepth = -1
    end

    representation = genotype._representation
    lenRep = length(representation)

    # Sets depths for every node
    depths = Array{Integer}(undef, lenRep)
    visitedFromRoot = zeros(Integer, 1, 2)
    nodeIndex = 1
    depths[nodeIndex] = 0
    parent = 0
    depth = 0

    while nodeIndex != 0

        if nodeIndex > size(visitedFromRoot)[1]
            visitedFromRoot = vcat(visitedFromRoot, zeros(Integer, 1, 2))
            depths[nodeIndex] = depth
            visitedFromRoot[end, 2] = parent
        end

        if visitedFromRoot[nodeIndex, 1] < getArity(representation[nodeIndex])
            parent = nodeIndex
            depth = depths[nodeIndex] + 1
            nodeIndex = size(visitedFromRoot)[1] + 1
        else
            nodeIndex = visitedFromRoot[nodeIndex, 2]
            if nodeIndex != 0
                visitedFromRoot[nodeIndex, 1] += 1
            end
        end
    end

    # Sets heights for every node
    heights = Array{Integer}(undef, lenRep)
    heights[end] = 0
    stack = [lenRep]
    currIndex = lenRep - 1

    while currIndex > 0

        if typeof(representation[currIndex]) <: TerminalNode
            heights[currIndex] = 0
        elseif typeof(representation[currIndex]) <: FunctionNode
            childrenIndexes = Array{Integer}(undef, 0)
            for i=1:getArity(representation[currIndex])
                push!(childrenIndexes, pop!(stack))
            end
            heights[currIndex] = maximum(heights[childrenIndexes]) + 1
        end

        push!(stack, currIndex)
        currIndex -= 1
    end

    # Chooses cross point
    indexes = Array{Integer}(undef, 0)

    if allowedType == Any
        indexes = collect(1:lenRep)
    else
        for i=1:lenRep
            if getType(representation[i]) <: allowedType
                push!(indexes, i)
            end
        end
    end

    indexes = Random.shuffle(rng, filter(x->!in(x, forbidden), indexes))
    pointIndex = 0

    if maxDepth != -1
        if fixedDepth != -1
            indexes = union(filter(x->depths[x]==fixedDepth, indexes),
                            filter(x->heights[x]+fixedDepth<=maxDepth, indexes))
        elseif minDepth != -1
            indexes = union(filter(x->depths[x]>=minDepth, indexes),
                            filter(x->heights[x]+minDepth<=maxDepth, indexes))
        end
    end

    if !isempty(indexes)
        for i in indexes
            random = rand(rng)
            prob = typeof(representation[i]) <: FunctionNode ? 0.9 : 0.1

            if random < prob
                pointIndex = i
                break
            end
        end

        if pointIndex == 0
            pointIndex = rand(rng, indexes)
        end

        # Subtree obtention
        subtree = Array{Node}(undef, 0)
        visitedFromRoot = zeros(Integer, 1, 2)
        push!(subtree, representation[pointIndex])
        nodeIndex = pointIndex
        parent = 0

        while nodeIndex != 0

            if (nodeIndex-pointIndex+1) > size(visitedFromRoot)[1]
                visitedFromRoot = vcat(visitedFromRoot, zeros(Integer, 1, 2))
                visitedFromRoot[end, 2] = parent
                push!(subtree, representation[nodeIndex])
            end

            if visitedFromRoot[nodeIndex-pointIndex+1, 1] < getArity(representation[nodeIndex])
                parent = nodeIndex
                nodeIndex = size(visitedFromRoot)[1] + pointIndex
            else
                nodeIndex = visitedFromRoot[nodeIndex-pointIndex+1, 2]
                if nodeIndex != 0
                    visitedFromRoot[nodeIndex-pointIndex+1, 1] += 1
                end
            end
        end

        return subtree, pointIndex, depths[pointIndex], getType(representation[pointIndex])
    else
        return nothing
    end
end

"""
    onePointCross(parent1::STGPGenotype, parent2::STGPGenotype,
                  gpExperimentInfo::STGPInfo, rng::Random.AbstractRNG)

Performs a crossover between two parent trees in which the crosspoints are at the
same depth. First, a random crosspoint is chosen in the first parent, and in the
second parent, a crosspoint of the same type and depth is searched. If there is
no such crosspoint in the second parent, the process repeats and a different
crosspoint is chosen in the first parent. If crossover is not possible between the
two parents (due to type and depth constrictions) the parents are returned as children.

`Self-provided Arguments` are provided by the library, so only `User Arguments` must be provided.

# Self-provided Arguments
- `parent1::CGPGenotype`: first parent for the crossover operation.
- `parent2::CGPGenotype`: first parent for the crossover operation.
- `gpExperimentInfo::CGPInfo`: information about the GP experiment.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

# User Arguments
None

# Returns
An array that contains the two resulting children.
"""
function onePointCross(parent1::STGPGenotype, parent2::STGPGenotype,
                       gpExperimentInfo::STGPInfo, rng::Random.AbstractRNG)

    maxDepth = gpExperimentInfo._maxTreeDepth
    forbidden = Array{Integer}(undef, 0)
    cantCross = false
    lenRep = length(getRepresentation(parent1))

    subtree1, crossPoint1, pointDepth1, allowedType = selectCrossPoint(parent1, rng)
    returned = selectCrossPoint(parent2, rng, fixedDepth=pointDepth1, maxDepth=maxDepth,
        allowedType=allowedType)
    while returned == nothing
        push!(forbidden, crossPoint1)
        if length(forbidden) == lenRep
            cantCross = true
            break
        end
        subtree1, crossPoint1, pointDepth1, allowedType = selectCrossPoint(parent1,
            rng, forbidden=forbidden)
        returned = selectCrossPoint(parent2, rng, fixedDepth=pointDepth1,
            maxDepth=maxDepth, allowedType=allowedType)
    end

    if !cantCross
        child1Tree = mergeSubtrees(parent1, returned[1], crossPoint1, length(subtree1))
        child2Tree = mergeSubtrees(parent2, subtree1, returned[2], length(returned[1]))
    else
        child1Tree = deepcopy(getRepresentation(parent1))
        child2Tree = deepcopy(getRepresentation(parent2))
    end

    return [STGPGenotype(child1Tree), STGPGenotype(child2Tree)]
end

"""
    oneChildSubtreeCross(parent1::STGPGenotype, parent2::STGPGenotype,
                         gpExperimentInfo::STGPInfo, rng::Random.AbstractRNG)

Performs the classical subtree cross between two parents but resulting in only
one child. Selects a random crosspoint in the first parent and then selects
another crosspoint in the second parent, knowing that it should not generate a
child that violates the maximum depth established and that it must be of the same
type as the crosspoint selected in the first parent. If crossover is not possible between the
two parents (due to type and depth constrictions) the parents are returned as children.

`Self-provided Arguments` are provided by the library, so only `User Arguments` must be provided.

# Self-provided Arguments
- `parent1::CGPGenotype`: first parent for the crossover operation.
- `parent2::CGPGenotype`: first parent for the crossover operation.
- `gpExperimentInfo::CGPInfo`: information about the GP experiment.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

# User Arguments
None

# Returns
An array that contains the resulting child.
"""
function oneChildSubtreeCross(parent1::STGPGenotype, parent2::STGPGenotype,
                              gpExperimentInfo::STGPInfo, rng::Random.AbstractRNG)

    maxDepth = gpExperimentInfo._maxTreeDepth
    forbidden = Array{Integer}(undef, 0)
    cantCross = false
    lenRep = length(getRepresentation(parent1))

    subtree1, crossPoint1, pointDepth1, allowedType = selectCrossPoint(parent1, rng)
    returned = selectCrossPoint(parent2, rng, minDepth=pointDepth1, maxDepth=maxDepth,
               allowedType=allowedType)
    while returned == nothing
        push!(forbidden, crossPoint1)
        if length(forbidden) == lenRep
            cantCross = true
            break
        end
        subtree1, crossPoint1, pointDepth1, allowedType = selectCrossPoint(parent1,
            rng, forbidden=forbidden)
        returned = selectCrossPoint(parent2, rng, minDepth=pointDepth1,
            maxDepth=maxDepth, allowedType=allowedType)
    end
    if !cantCross
        childTree = mergeSubtrees(parent1, returned[1], crossPoint1, length(subtree1))
    else
        childTree = deepcopy(getRepresentation(parent1))
    end

    return [STGPGenotype(childTree)]
end

"""
    subtreeCross(parent1::STGPGenotype, parent2::STGPGenotype,
                 gpExperimentInfo::STGPInfo, rng::Random.AbstractRNG)

Performs the classical subtree cross between two parents. Selects a random
crosspoint in the first parent and then selects another crosspoint in the second
parent, knowing that it should not generate a child that violates the maximum
depth established and that it must be of the same type as the crosspoint selected
in the first parent. Then, repeats the same operation but taking the second parent
as it was the first one. If crossover is not possible between the
two parents (due to type and depth constrictions) the parents are returned as children.

`Self-provided Arguments` are provided by the library, so only `User Arguments` must be provided.

# Self-provided Arguments
- `parent1::CGPGenotype`: first parent for the crossover operation.
- `parent2::CGPGenotype`: first parent for the crossover operation.
- `gpExperimentInfo::CGPInfo`: information about the GP experiment.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

# User Arguments
None

# Returns
An array that contains the two resulting children.
"""
function subtreeCross(parent1::STGPGenotype, parent2::STGPGenotype,
                      gpExperimentInfo::STGPInfo, rng::Random.AbstractRNG)

    child1 = oneChildSubtreeCross(parent1, parent2, gpExperimentInfo, rng)[1]
    child2 = oneChildSubtreeCross(parent2, parent1, gpExperimentInfo, rng)[1]

    return [child1, child2]
end

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
# simple Julia APIs over Kubernetes Swagger interface

function sel(label::String, op::Symbol)
    @assert op === :exists
    label
end
sel(label::String, op::Symbol, items::String...) = label * " " * string(op) * " (" * join(items, ",") * ")"
sel(cnd::String...) = join(cnd, ", ")

_delopts(; kwargs...) = Typedefs.MetaV1.DeleteOptions(; preconditions=Typedefs.MetaV1.Preconditions(; kwargs...), kwargs...)

function get(ctx::KuberContext, O::Symbol, name::String; kwargs...)
    isempty(ctx.apis) && set_api_versions!(ctx)

    kapi = ctx.modelapi[O]
    apictx = kapi.api(ctx.client)

    try
        apicall = eval(Symbol("read$O"))
        return apicall(apictx, name; kwargs...)
    catch ex
        isa(ex, UndefVarError) || rethrow()
        apicall = eval(Symbol("readNamespaced$O"))
        return apicall(apictx, name, ctx.namespace; kwargs...)
    end
end

function get(ctx::KuberContext, O::Symbol; label_selector=nothing, namespace::Union{String,Nothing}=ctx.namespace)
    isempty(ctx.apis) && set_api_versions!(ctx)

    kapi = ctx.modelapi[O]
    apictx = kapi.api(ctx.client)

    try
        apiname = "list$O"
        (namespace === nothing) && (apiname *= "ForAllNamespaces")
        apicall = eval(Symbol(apiname))
        return apicall(apictx; labelSelector=label_selector)
    catch ex
        isa(ex, UndefVarError) || rethrow()
        (namespace === nothing) && rethrow()
        apicall = eval(Symbol("listNamespaced$O"))
        return apicall(apictx, namespace; labelSelector=label_selector)
    end
end

function put!(ctx::KuberContext, v::T) where {T<:SwaggerModel}
    vjson = convert(Dict{String,Any}, v)
    put!(ctx, Symbol(vjson["kind"]), vjson)
end

function put!(ctx::KuberContext, O::Symbol, d::Dict{String,Any})
    isempty(ctx.apis) && set_api_versions!(ctx)

    kapi = ctx.modelapi[O]
    apictx = kapi.api(ctx.client)

    try
        apicall = eval(Symbol("create$O"))
        return apicall(apictx, d)
    catch ex
        isa(ex, UndefVarError) || rethrow()
        apicall = eval(Symbol("createNamespaced$O"))
        return apicall(apictx, ctx.namespace, d)
    end
end

function delete!(ctx::KuberContext, v::T; kwargs...) where {T<:SwaggerModel}
    vjson = convert(Dict{String,Any}, v)
    kind = vjson["kind"]
    name = vjson["metadata"]["name"]
    delete!(ctx, Symbol(kind), name; kwargs...)
end

function delete!(ctx::KuberContext, O::Symbol, name::String; kwargs...)
    isempty(ctx.apis) && set_api_versions!(ctx)

    kapi = ctx.modelapi[O]
    apictx = kapi.api(ctx.client)
    params = [apictx, name]

    try
        apicall = eval(Symbol("delete$O"))
        return apicall(params...; kwargs...)
    catch ex
        isa(ex, UndefVarError) || rethrow()
        apicall = eval(Symbol("deleteNamespaced$O"))
        push!(params, ctx.namespace)
        return apicall(params...; kwargs...)
    end
end

function update!(ctx::KuberContext, v::T, patch, patch_type) where {T<:SwaggerModel}
    vjson = convert(Dict{String,Any}, v)
    kind = vjson["kind"]
    name = vjson["metadata"]["name"]
    update!(ctx, Symbol(kind), name, patch, patch_type)
end

function update!(ctx::KuberContext, O::Symbol, name::String, patch, patch_type)
    isempty(ctx.apis) && set_api_versions!(ctx)

    kapi = ctx.modelapi[O]
    apictx = kapi.api(ctx.client)

    try
        apicall = eval(Symbol("patch$O"))
        return apicall(apictx, name, patch; _mediaType=patch_type)
    catch ex
        isa(ex, UndefVarError) || rethrow()
        apicall = eval(Symbol("patchNamespaced$O"))
        return apicall(apictx, name, ctx.namespace, patch; _mediaType=patch_type)
    end
end

"""
Fetch logs of the specified Pod

Parameters:
- name::String (required)

Keyword Args:
- container::String : The container for which to stream logs. Defaults to only container if there is one container in the pod.
- follow::Bool : Follow the log stream of the pod. Defaults to false.
- limitBytes::Int32 : If set, the number of bytes to read from the server before terminating the log output. This may not display a complete final line of logging, and may return slightly more or slightly less than the specified limit.
- pretty::String : If 'true', then the output is pretty printed.
- previous::Bool : Return previous terminated container logs. Defaults to false.
- sinceSeconds::Int32 : A relative time in seconds before the current time from which to show logs. If this value precedes the time a pod was started, only logs since the pod start will be returned. If this value is in the future, no logs will be returned. Only one of sinceSeconds or sinceTime may be specified.
- sinceTime::String : An RFC3339 timestamp from which to show logs. If this value precedes the time a pod was started, only logs since the pod start will be returned. If this value is in the future, no logs will be returned. Only one of sinceSeconds or sinceTime may be specified.
- tailLines::Int32 : If set, the number of lines from the end of the logs to show. If not specified, logs are shown from the creation of the container or sinceSeconds or sinceTime
- timestamps::Bool : If true, add an RFC3339 or RFC3339Nano timestamp at the beginning of every line of log output. Defaults to false.

Returns: String of all log entries, one per line
"""
get_logs(ctx::KuberContext, pod_name::String; kwargs...) = get(ctx, :PodLog, pod_name; kwargs...)


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

JLLWrappers.@generate_main_file_header("PDAL")
JLLWrappers.@generate_main_file("PDAL", UUID("a8197b14-d70b-5660-b10f-8b1ebb62825c"))
end  # module PDAL_jll


text:
function TriaGrid(nrows, ncols)::Grid{Tria,NonPeriodic}
    return Grid(Tria, NonPeriodic, nrows, ncols)
end

function PeriodicTriaGrid(nrows, ncols)::Grid{Tria,Periodic}
    @assert iseven(ncols)
    return Grid(Tria, Periodic, nrows, ncols)
end


text:
struct ConstrainedAquisition{A<:EGOAquisition,N} <: EGOAquisition
    base::A
    cons::NTuple{N,Function}
end

update_parameters!(a::ConstrainedAquisition, args...) = update_parameters!(a.base, args...)

function objective(a::ConstrainedAquisition, model::Kriging, x)
    for g in a.cons
        if g(x) > eps()
            return 0.0
        end
    end
    return objective(a.base, model, x)
end


text:
function iteration_time(filename, iter)
    file = jldopen(filename)
    time = file["timeseries/t/$iter"]
    close(file)
    return time
end

function field_data(filename, field_name, iter)
    file = jldopen(filename)
    data = file["timeseries/$field_name/$iter"]
    close(file)
    return data
end

function RegularCartesianGrid(filename)
    file = jldopen(filename)
    Nx = file["grid/Nx"]
    Ny = file["grid/Ny"]
    Nz = file["grid/Nz"]
    Lx = file["grid/Lx"]
    Ly = file["grid/Ly"]
    Lz = file["grid/Lz"]
    close(file)

    grid = RegularCartesianGrid(size=(Nx, Ny, Nz), x=(0, Lx), y=(0, Ly), z=(0, Lz),
                                topology=(Periodic, Periodic, Bounded))

    return grid
end

function iterations(filename)
    file = jldopen(filename)
    iters = parse.(Int, keys(file["timeseries/t"]))
    close(file)

    return iters
end


text:
if !@isdefined USE_GPU
    const USE_GPU = true
end

using DataFrames
using CSV
using ImplicitGlobalGrid
import MPI
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
if !ParallelStencil.is_initialized()
    @static if USE_GPU
        @init_parallel_stencil(CUDA, Float64, 2)
    else
        @init_parallel_stencil(Threads, Float64, 2)
    end
end
include("part2.jl")

const FILENAME = joinpath("benchmark-results", "part2_semi_implicit_vs_explicit_experiment_results.csv")

const NX, NY = begin
    local nx, ny
    nx = 2^11
    ny = (nx ÷ 4)
    (nx, ny) .+ 1
end
@assert NX == (2^11 + 1)
@assert NY == (2^9 + 1)

# Pr_{liquid metals} ≈ 0.004 to 0.03
# Pr_{air} ≈ 0.7
# Pr_{water} ≈ 1.7 to 12
# Pr_{rock} ≈ practically infinite, i.e. 10^24
for Pr in 10 .^(-3.0:1.0:1.0)
    for beta in [0.0, 0.5, 1.0]
        println("Running for Pr=$(Pr) and beta=$(beta).")
        opt = SimIn_t()
        opt.nx = NX
        opt.ny = NY
        opt.Pr = Pr
        opt.tol = 1.0e-7
        opt.beta = beta
        opt.ttot = 0.005

        res = navier_stokes_2D(opt=opt)

        bench_row = (nx=NX, ny=NY, Pr=Pr, beta=beta, t_elapsed=res.t_elapsed, timed_iters=res.timed_iters)

        bench_df =
            if isfile(FILENAME)
                DataFrame(CSV.File(FILENAME))
            else
                DataFrame()
            end
        push!(bench_df, bench_row)

        open(FILENAME, "w") do io
            CSV.write(io, bench_df)
        end
    end
end


text:
using Printf
using NPZ
include("spectre-defense/util.jl")
include("spectre-defense/quantum_filters.jl")

input_folder = ARGS[1]
output_path = ARGS[2]
target_label = ARGS[3]
eps = parse(Int16,ARGS[4])

reps = npzread(input_folder * "$(target_label).npy")'
n = size(reps)[2]
removed = round(Int, 1.5*eps)

println("Running quantum filter")
quantum_poison_ind = .! rcov_auto_quantum_filter(reps, eps)
poison_removed = sum(quantum_poison_ind[end-eps+1:end])
clean_removed = removed - poison_removed
@printf("%d poisons removed, %d cleans removed\n", poison_removed, clean_removed)
npzwrite(output_path, quantum_poison_ind)

max_stars_repo_name: SewoongLab/backdoor-suite

max_issues_repo_name: SewoongLab/backdoor-suite

max_forks_repo_name: SewoongLab/backdoor-suite


text:
using POMDPs

POMDPs.add("POMDPToolbox")

Pkg.clone("https://github.com/zsunberg/LightDarkPOMDPs.jl")


text:
#
# Plot Circular Restricted Three-body Problem trajectories
#

"""
Plot the orbital positions of a CR3BP orbit.
"""
function plotpositions(traj::Trajectory{<:CircularRestrictedThreeBodyOrbit}; lengthunit = lengthunit(first(traj).system), exclude_z = false, kwargs...)

    all_true_or_false(vec) = all(vec) || all(map(v->!v, vec))
    @assert all_true_or_false(map(orbit -> coordinateframe(orbit.state), traj) isa Union{Synodic, Inertial}) "All coordinate frames within a trajectory must be identical for plotting!"
    @assert all_true_or_false(map(orbit -> orbit isa NormalizedCartesianState, traj)) "If one orbit `state` is normalized, all must be!"
    @assert all([isequal(traj[i].system, traj[i+1].system) for i ∈ 1:length(traj)-1]) "All orbits with the trajectory must have the same CR3BP system."

    isnormalized = first(traj).state isa NormalizedCartesianState
    frame        = coordinateframe(first(traj).state)

    pos = transpose(hcat(position_vector.(traj)...))

    if !isnormalized
        pos = ustrip.(lengthunit, pos)
        LU  = string(lengthunit)
    else
        LU  = string(normalized_length_unit(first(traj).system))
    end

    defaults = (; title        = "$(string(frame)) Positions" * (first(traj).system.name == "" ? "" : " w.r.t. $(first(traj).system.name) Barycenter"),
                  xlabel       = "X ($LU)",
                  ylabel       = "Y ($LU)",
                  zlabel       = "Z ($LU)",
                  label        = "Orbit Position",
                  formatter    = :plain,
                  grid         = :on,
                  linewidth    = 2,
                  dpi          = 150)

    options = merge(defaults, kwargs)

    if !exclude_z
        return plot(pos[:,1], pos[:,2], pos[:,3]; options...)
    else
        return plot(pos[:,1], pos[:,2]; options...)
    end
end

"""
Plot the orbital positions of a CR3BP orbit to the last plot.
"""
function plotpositions!(fig, traj::Trajectory{<:CircularRestrictedThreeBodyOrbit}; lengthunit = lengthunit(first(traj).system), exclude_z = false, kwargs...)

    all_true_or_false(vec) = all(vec) || all(map(v->!v, vec))
    @assert all_true_or_false(map(orbit -> coordinateframe(orbit.state), traj) isa Union{Synodic, Inertial}) "All coordinate frames within a trajectory must be identical for plotting!"
    @assert all_true_or_false(map(orbit -> orbit isa NormalizedCartesianState, traj)) "If one orbit `state` is normalized, all must be!"
    @assert all([isequal(traj[i].system, traj[i+1].system) for i ∈ 1:length(traj)-1]) "All orbits with the trajectory must have the same CR3BP system."

    isnormalized = first(traj).state isa NormalizedCartesianState
    frame        = coordinateframe(first(traj).state)

    pos = transpose(hcat(position_vector.(traj)...))

    if !isnormalized
        pos = ustrip.(lengthunit, pos)
        LU  = string(lengthunit)
    else
        LU  = string(normalized_length_unit(first(traj).system))
    end

    defaults = (; linewidth = 2, label = :none)
    options = merge(defaults, kwargs)

    if !exclude_z
        return plot!(fig, pos[:,1], pos[:,2], pos[:,3]; options...)
    else
        return plot!(fig, pos[:,1], pos[:,2]; options...)
    end
end

"""
Plot the orbital velocities of a CR3BP orbit.
"""
function plotvelocities(traj::Trajectory{<:CircularRestrictedThreeBodyOrbit}; velocityunit = velocityunit(first(traj).system), exclude_z = false, kwargs...)

    all_true_or_false(vec) = all(vec) || all(map(v->!v, vec))
    @assert all_true_or_false(map(orbit -> coordinateframe(orbit.state), traj) isa Union{Synodic, Inertial}) "All coordinate frames within a trajectory must be identical for plotting!"
    @assert all_true_or_false(map(orbit -> orbit isa NormalizedCartesianState, traj)) "If one orbit `state` is normalized, all must be!"
    @assert all([isequal(traj[i].system, traj[i+1].system) for i ∈ 1:length(traj)-1]) "All orbits with the trajectory must have the same CR3BP system."

    isnormalized = first(traj).state isa NormalizedCartesianState
    frame        = coordinateframe(first(traj).state)

    vel = transpose(hcat(velocity_vector.(traj)...))

    if !isnormalized
        vel = ustrip.(velocityunit, vel)
        VU  = string(velocityunit)
    else
        VU  = string(normalized_length_unit(first(traj).system) / normalized_time_unit(first(traj).system))
    end

    defaults = (; title        = "$(string(frame)) Velocities" * (first(traj).system.name == "" ? "" : " w.r.t. $(first(traj).system.name) Barycenter"),
                  xlabel       = "X ($VU)",
                  ylabel       = "Y ($VU)",
                  zlabel       = "Z ($VU)",
                  label        = "Orbit Velocity",
                  formatter    = :plain,
                  grid         = :on,
                  linewidth    = 2,
                  dpi          = 150)

    options = merge(defaults, kwargs)

    if !exclude_z
        return plot(vel[:,1], vel[:,2], vel[:,3]; options...)
    else
        return plot(vel[:,1], vel[:,2]; vel...)
    end
end

"""
Plot the orbital velocities of a CR3BP orbit to the last plot.
"""
function plotvelocities!(fig, traj::Trajectory{<:CircularRestrictedThreeBodyOrbit}; velocityunit = velocityunit(first(traj).system), exclude_z = false, kwargs...)

    all_true_or_false(vec) = all(vec) || all(map(v->!v, vec))
    @assert all_true_or_false(map(orbit -> coordinateframe(orbit.state), traj) isa Union{Synodic, Inertial}) "All coordinate frames within a trajectory must be identical for plotting!"
    @assert all_true_or_false(map(orbit -> orbit isa NormalizedCartesianState, traj)) "If one orbit `state` is normalized, all must be!"
    @assert all([isequal(traj[i].system, traj[i+1].system) for i ∈ 1:length(traj)-1]) "All orbits with the trajectory must have the same CR3BP system."

    isnormalized = first(traj).state isa NormalizedCartesianState
    frame        = coordinateframe(first(traj).state)

    vel = transpose(hcat(velocity_vector.(traj)...))

    if !isnormalized
        vel = ustrip.(velocityunit, vel)
        VU  = string(velocityunit)
    else
        VU  = string(normalized_length_unit(first(traj).system) / normalized_time_unit(first(traj).system))
    end

    defaults = (; linewidth = 2, label = :none)
    options = merge(defaults, kwargs)

    if !exclude_z
        return plot!(fig, vel[:,1], vel[:,2], vel[:,3]; options...)
    else
        return plot!(fig, vel[:,1], vel[:,2]; vel...)
    end
end

"""
Plot the nondimensional positions of
CR3BP celestial bodies.
"""
function plotbodies(system::CircularRestrictedThreeBodySystem; bodies = ("", ""), normalize = true, exclude_z = true, kwargs...)

    μ = normalized_mass_parameter(system)

    if all(isempty, bodies)
        if occursin("-", system.name)
            bodies = (string.(split(system.name, "-"))...,)
        else
            bodies = ("Body 1", "Body 2")
        end
    end

    if normalize
        r₁ = [[-μ],  [0], [0]]
        r₂ = [[1-μ], [0], [0]]
        LU = string(normalized_length_unit(system))
    else
        r₁ = [[ustrip(length_unit(system), -μ * normalized_length_unit(system))],  [0], [0]]
        r₂ = [[ustrip(length_unit(system), 1-μ * normalized_length_unit(system))], [0], [0]]
        LU = string(lengthunit(system))
    end

    defaults = (; title = "CR3BP Bodies", xlabel = "X ($LU)", ylabel = "Y ($LU)", zlabel = "Z ($LU)")
    options  = merge(defaults, kwargs)

    if exclude_z
        fig = scatter(r₁[1:2]...; label=bodies[1], options...)
        scatter!(fig, r₂[1:2]...; label=bodies[2], options...)
        plot!(fig; options...)
        return fig
    else
        fig = scatter(r₁...; markersize=10, markercolor=:lightblue, label="Body 1")
        scatter!(fig, r₂...; markersize=6, markercolor=:gray, label="Body 2")
        plot!(fig; options...)
        return fig
    end

end

"""
Plot the nondimensional positions of
CR3BP celestial bodies to the last figure.
"""
function plotbodies!(fig, system::CircularRestrictedThreeBodySystem;
                     bodies = ("", ""),
                     normalize = true,
                     exclude_z = true,
                     kwargs...)

    μ = normalized_mass_parameter(system)

    if all(isempty, bodies)
        if occursin("-", system.name)
            bodies = (string.(split(system.name, "-"))...,)
        else
            bodies = ("Body 1", "Body 2")
        end
    end

    if normalize
        r₁ = [[-μ],  [0], [0]]
        r₂ = [[1-μ], [0], [0]]
        LU = string(normalized_length_unit(system))
    else
        r₁ = [[ustrip(length_unit(system), -μ * normalized_length_unit(system))],  [0], [0]]
        r₂ = [[ustrip(length_unit(system), 1-μ * normalized_length_unit(system))], [0], [0]]
        LU = string(lengthunit(system))
    end

    if exclude_z
        scatter!(fig, r₁[1:2]...; label = bodies[1], kwargs...)
        scatter!(fig, r₂[1:2]...; label = bodies[2], kwargs...)
        return fig
    else
        scatter!(fig, r₁...; label = bodies[1], kwargs...)
        scatter!(fig, r₂...; label = bodies[2], kwargs...)
        return fig
    end

end


text:
module HallOfFameModule

import ..CoreModule: CONST_TYPE, MAX_DEGREE, Node, Options, Dataset, string_tree
import ..EquationUtilsModule: compute_complexity
import ..PopMemberModule: PopMember, copy_pop_member
import ..LossFunctionsModule: eval_loss
using Printf: @sprintf

""" List of the best members seen all time in `.members` """
mutable struct HallOfFame
    members::Array{PopMember,1}
    exists::Array{Bool,1} #Whether it has been set

    # Arranged by complexity - store one at each.
end

"""
    HallOfFame(options::Options)

Create empty HallOfFame. The HallOfFame stores a list
of `PopMember` objects in `.members`, which is enumerated
by size (i.e., `.members[1]` is the constant solution).
`.exists` is used to determine whether the particular member
has been instantiated or not.
"""
function HallOfFame(options::Options)
    actualMaxsize = options.maxsize + MAX_DEGREE
    return HallOfFame(
        [
            PopMember(
                Node(convert(CONST_TYPE, 1)),
                1.0f9,
                1.0f9;
                parent=-1,
                deterministic=options.deterministic,
            ) for i in 1:actualMaxsize
        ],
        [false for i in 1:actualMaxsize],
    )
end

"""
    calculate_pareto_frontier(dataset::Dataset{T}, hallOfFame::HallOfFame,
                            options::Options) where {T<:Real}
"""
function calculate_pareto_frontier(
    dataset::Dataset{T}, hallOfFame::HallOfFame, options::Options
)::Array{PopMember,1} where {T<:Real}
    # Dominating pareto curve - must be better than all simpler equations
    dominating = PopMember[]
    actualMaxsize = options.maxsize + MAX_DEGREE
    for size in 1:actualMaxsize
        if !hallOfFame.exists[size]
            continue
        end
        member = hallOfFame.members[size]
        # We check if this member is better than all members which are smaller than it and
        # also exist.
        betterThanAllSmaller = true
        for i in 1:(size - 1)
            if !hallOfFame.exists[i]
                continue
            end
            simpler_member = hallOfFame.members[i]
            if member.loss >= simpler_member.loss
                betterThanAllSmaller = false
                break
            end
        end
        if betterThanAllSmaller
            push!(dominating, copy_pop_member(member))
        end
    end
    return dominating
end

"""
    calculate_pareto_frontier(X::AbstractMatrix{T}, y::AbstractVector{T},
                            hallOfFame::HallOfFame, options::Options;
                            weights=nothing, varMap=nothing) where {T<:Real}

Compute the dominating Pareto frontier for a given hallOfFame. This
is the list of equations where each equation has a better loss than all
simpler equations.
"""
function calculate_pareto_frontier(
    X::AbstractMatrix{T},
    y::AbstractVector{T},
    hallOfFame::HallOfFame,
    options::Options;
    weights=nothing,
    varMap=nothing,
) where {T<:Real}
    return calculate_pareto_frontier(
        Dataset(X, y; weights=weights, varMap=varMap), hallOfFame, options
    )
end

function string_dominating_pareto_curve(hallOfFame, baselineMSE, dataset, options, avgy)
    output = ""
    curMSE = Float64(baselineMSE)
    lastMSE = curMSE
    lastComplexity = 0
    output *= "Hall of Fame:\n"
    output *= "-----------------------------------------\n"
    output *= @sprintf(
        "%-10s  %-8s   %-8s  %-8s\n", "Complexity", "loss", "Score", "Equation"
    )

    dominating = calculate_pareto_frontier(dataset, hallOfFame, options)
    for member in dominating
        complexity = compute_complexity(member.tree, options)
        if member.loss < 0.0
            throw(
                DomainError(
                    member.loss,
                    "Your loss function must be non-negative. To do this, consider wrapping your loss inside an exponential, which will not affect the search (unless you are using annealing).",
                ),
            )
        end
        # User higher precision when finding the original loss:
        relu(x) = x < 0 ? 0 : x
        curMSE = member.loss

        delta_c = complexity - lastComplexity
        ZERO_POINT = 1e-10
        delta_l_mse = log(abs(curMSE / lastMSE) + ZERO_POINT)
        score = convert(Float32, -delta_l_mse / delta_c)
        output *= @sprintf(
            "%-10d  %-8.3e  %-8.3e  %-s\n",
            complexity,
            curMSE,
            score,
            string_tree(member.tree, options, varMap=dataset.varMap)
        )
        lastMSE = curMSE
        lastComplexity = complexity
    end
    output *= "\n"
    return output
end

end


text:
using ForwardDiff
using LinearAlgebra
include("poisson_linear_system.jl")
include("gradient_descent.jl")
include("learned_jacobi.jl")

function optimize_parameters(
    numelements,
    wavenumber,
    guessparams;
    lr = 1.0,
    gditers = 20,
)
    numiter = numelements^2

    xL, xR = -1.0, 1.0
    bcleft, bcright = 0.0, 0.0

    stepsize = (xR - xL) / numelements
    numpts = numelements + 1
    xrange = range(xL, stop = xR, length = numpts)

    rhs = cos.(wavenumber * xrange)

    op = -poisson_linear_system(numpts)
    exactsol = solve_direct(op, rhs, bcleft, bcright, stepsize)

    func(x) = LearnedJacobi.error_for_parameters(
        bcleft,
        bcright,
        rhs,
        stepsize,
        numiter,
        exactsol,
        x,
    )
    grad(x) = ForwardDiff.gradient(func, x)

    optimumparams = run_gradient_descent_iterations(func, grad, p0, lr, gditers)

    itersol = LearnedJacobi.run_parametric_iterations(
        bcleft,
        bcright,
        rhs,
        stepsize,
        numiter,
        optimumparams,
    )

    return optimumparams, itersol, exactsol
end

numelements = 32
# jacobiparams = [0.5, 0.0, 0.5, 0.5]
p0 = [0.0, 0.0, 0.0, 0.0]
wavenumber = pi / 2

optimumparams, itersol, exactsol = optimize_parameters(numelements,wavenumber,p0)

max_stars_repo_name: ArjunNarayanan/ExploreJacobi

max_issues_repo_name: ArjunNarayanan/ExploreJacobi

max_forks_repo_name: ArjunNarayanan/ExploreJacobi


text:
# 1D CA using Agents.jl
module CA1D
using Agents

mutable struct Agent{T<:Integer, Y<:AbstractString} <: AbstractAgent
  id::T
  pos::Tuple{T, T}
  status::Y
end

mutable struct Model{T<:AbstractVector, Y<:AbstractDict, X<:AbstractSpace} <: AbstractModel
  space::X
  agents::T  #Array{AbstractAgent}
  scheduler::Function
  rules::Y
end

mutable struct Space{T<:Integer, Y<:AbstractVector} <: AbstractSpace
  dimensions::Tuple{T, T}
  space::SimpleGraph
  agent_positions::Y  # an array of arrays for each grid node
end

"""
    build_model(;rules::Dict, ncols=101)

Builds a 1D cellular automaton. `rules` is a dictionary with this format: `Dict("000" => "0")`. `ncols` is the number of cells in a 1D CA.
"""
function build_model(;rules::Dict, ncols::Integer=101)
  gridsize=(ncols,1)
  agents = [Agent(i, (i,1), "0") for i in 1:gridsize[1]]
  agent_positions = [Array{Integer}(undef, 0) for i in 1:gridsize[1]]
  mygrid = Space(gridsize, grid(gridsize, true, false), agent_positions)
  # mygrid = Space(gridsize, grid(gridsize, false, false), agent_positions)  # this is for when there space is not toroidal
  model = Model(mygrid, agents, as_added, rules)
  return model
end

function ca_step!(model)
  agentnum = nagents(model)
  new_status = Array{String}(undef, agentnum)
  for agid in 1:agentnum
    agent = model.agents[agid]
    center = agent.status
    if agent.id == agentnum
      # right = "0"  # this is for when there space is not toroidal
      right =  model.agents[1].status
      left = model.agents[agent.id-1].status
    elseif agent.id == 1
      # left = "0"  # this is for when there space is not toroidal
      left = model.agents[agentnum].status
      right = model.agents[agent.id+1].status
    else
      left = model.agents[agent.id-1].status
      right = model.agents[agent.id+1].status
    end
    rule = left*center*right
    new_status[agid] = model.rules[rule]
  end
  for ss in 1:agentnum
    model.agents[ss].status = new_status[ss]
  end
end

"""
    ca_run(model::AbstractModel, runs::Integer, filename::String="CA_1D")

Runs a 1D cellular automaton.
"""
function ca_run(model::AbstractModel, runs::Integer, filename::String="CA_1D")
  data = step!(dummystep, CA1D.ca_step!, model, runs, [:pos, :status], collect(1:runs))
  visualize_1DCA(data, model, :pos, :status, runs, savename=filename)
end

end  # module


text:
"""
stock.jl
Functions for structured stock-level model components
"""

function harvest!(effort::Float64, stockdb::StockDB, stockassumptions::StockAssumptions)
  """
  This function will generate fishing mortality based on a specified integer level of effort, and assumed age specific catchabilities.
  """
  harvest_size = fill(0, size(stockassumptions.catchability))
  if effort > 0
    for i = 1:length(stockassumptions.catchability)
      harvest_size[i] = rand(Poisson(stockdb.population[end,i]*stockassumptions.catchability[i]*effort))
      stockdb.population[end, i] -= harvest_size[i]
    end
  end
  push!(stockdb.harvest, DataArray(harvest_size))
end

function ageadults!(stockdb::StockDB, stockassumptions::StockAssumptions, carryingcapacity::Float64, bump=0::Int64)
  """
  This function will apply transition probabilities to the current adult population.
  """
  stock_size = fill(0, size(stockassumptions.naturalmortality))
  if isnan(stockassumptions.mortalitycompensation)
    compensation_factor = 1
  else
    compensation_factor = 2*(cdf(Normal(carryingcapacity, carryingcapacity/stockassumptions.mortalitycompensation), sum(stockdb.population[end,:][1,])))
  end
  @assert(0.01 < compensation_factor < 1.99, "Population regulation has failed, respecify simulation parameters")
  for i = 1:(length(stockassumptions.naturalmortality)-1)
    stock_size[i+1] = rand(Binomial(stockdb.population[end,i], 1-(stockassumptions.naturalmortality[i]*compensation_factor)))
  end
  stock_size[end] += rand(Binomial(stockdb.population[end,end], 1-(stockassumptions.naturalmortality[end]*compensation_factor)))
  stock_size[1] += bump
  push!(stockdb.population, DataArray(stock_size))
end


text:
module DecisionMakingEnvironments

import Base:length, push!
import StatsBase: sample, sample!
using Distributions
using RecipesBase

abstract type AbstractDecisionMakingProblem end

"""
    sample

For a Bandit problem this samples the reward given an action and posibily context.
For a squential decision making problem this samples the next state, observation, reward, and terminal check given an action and current state.
"""
function sample end
function sample! end

"""
    sample_trajectory!

records a trajectory for a sequential decision making problem
"""
function sample_trajectory! end

"""
    sample_objective

    returns a sample of the objective function given a policy
"""
function sample_objective end

"""
    meta_information

returns problem meta-information useful for initializing a solver
"""
function meta_information end

"""
    solve(prob::AbstractDecisionMakingProblem, alg)

solves the decision making problem using the supplied algorithm. Returns a policy and any other algorithm specific properties
"""
function solve end

export AbstractDecisionMakingProblem
export sample, sample!, sample_trajectory!, sample_objective, meta_information, solve

export AbstractBanditProblem, BanditProblem, ContextualBanditProblem
export AbstractMDP, MDP, POMDP, SequentialProblem
export BanditExperience, ContextualBanditExperience, Trajectory
export length, push!, finish!, finished

# export create_simple_discrete_bandit, create_simple_discrete_contextualbandit
# export create_simple_chain, create_simple_chain_finitetime
# export create_finitetime_cartpole

include("bandits/bandit.jl")
include("mdps/mdp.jl")
include("history.jl")

end

# What is a decision making problem
# Optimization problem that looks to find a policy that maximizes return.
# Solution space: policy that maps states (possibly a history) to a distribtuion over the action space. No restrictions here.
# Objective function - expectation so must be able to sample.
# need a way to sample next states
# need a method to sample rewards
# need a way to identify if a squence has terminated and should be reset
# needs a method of getting the initial state or sampling one
# needs a discount function. Default to 1 for finite horizon.
# state space
# action space

# meta information:
# problem type: finite horizon, POMDP, infinite horizon, Bandit, Contextual Bandit
#

# decision making problem: {State space, action space, transition, reward function, initial state distribution, discount function}

# sample(prob, inputs...) get next state reward, termination
# sample_trajectory(prob, pi) record sequence of state actions rewards, logprobs
# sample_objective
# meta_info
# solve(prob, alg) return policy and other algorithms outputs

# problems I want to solve:
# MDP discrete / continuous
# POMDP discrete / continuous
#

# state space, action type, obs type
# meta-info

# space description
# tuple each element contains upper and lower bounds


text:
# ---
# title: 1250. Check If It Is a Good Array
# id: problem1250
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math
# link: <https://leetcode.com/problems/check-if-it-is-a-good-array/description/>
# hidden: true
# ---
#
# Given an array `nums` of positive integers. Your task is to select some subset
# of `nums`, multiply each element by an integer and add all these numbers. The
# array is said to be  **good  **if you can obtain a sum of `1` from the array
# by any possible subset and multiplicand.
#
# Return `True` if the array is **good  **otherwise return `False`.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [12,5,7,23]
#     Output: true
#     Explanation: Pick numbers 5 and 7.
#     5*3 + 7*(-2) = 1
#
#
# **Example 2:**
#
#
#
#     Input: nums = [29,6,10]
#     Output: true
#     Explanation: Pick numbers 29, 6 and 10.
#     29*1 + 6*(-3) + 10*(-1) = 1
#
#
# **Example 3:**
#
#
#
#     Input: nums = [3,6]
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 10^5`
#   * `1 <= nums[i] <= 10^9`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
const Dat = Float32
using Printf, Statistics, LinearAlgebra, DelimitedFiles, Plots
# pyplot()

run_id = "run_ds_40"

outdir = "../output/$run_id"
function load_array(A, A_name, extension, isave, outd)
    fid=open("$outd/$(isave)_$(A_name).$(extension)", "r"); read!(fid, A); close(fid)
end

Info = readdlm("$outdir/0_infos.inf")

Lx    = Info[1]
Ly    = Info[2]
nx    = Int(Info[3])
ny    = Int(Info[4])
isave = Int(Info[5])
err   = Info[6]

Mask = zeros(Dat, nx, ny)
xcp  = zeros(Dat, nx, 1 )
ycp  = zeros(Dat, ny, 1 )
F    = zeros(Dat, nx, ny)

load_array(Mask, "M"   , "res", 0, outdir)
load_array(xcp , "xcp" , "res", 0, outdir)
load_array(ycp , "ycp" , "res", 0, outdir)

load_array(F   , "F"   , "res", isave, outdir)

# visu
F[Mask.==0] .= NaN
xcp = xcp[:,1]
ycp = ycp[:,1]

plt = heatmap(xcp, ycp, log10.(F'), aspect_ratio=1, xlims=(xcp[1], xcp[end]), ylims=(ycp[1], ycp[end]), c=:viridis, framestyle=:box, title="log10 ||Flux|| (ϵ=$err)")

savefig(plot(plt, dpi=300), joinpath(@__DIR__, "$outdir/o$(length(xcp))_$isave.png"))

display(plt)


text:
struct HamiltonianProblem{iip} <: DiffEqBase.AbstractDynamicalODEProblem end

function HamiltonianProblem(H,q0,p0,tspan,p=nothing;kwargs...)
  iip = (typeof(q0) <: AbstractArray) && !(typeof(q0) <: SArray)
  HamiltonianProblem{iip}(H,q0,p0,tspan,p;kwargs...)
end

struct PhysicsTag end

function HamiltonianProblem{T}(H,p0,q0,tspan,p=nothing;kwargs...) where T
    if T == false
        if typeof(q0) <: Number
            dp = function (v,x,p,t)
                ForwardDiff.derivative(x->-H(v, x, p), x)
            end
            dq = function (v,x,p,t)
                ForwardDiff.derivative(v->H(v, x, p), v)
            end

            return ODEProblem(DynamicalODEFunction{T}(dp,dq), ArrayPartition(p0,q0), tspan, p; kwargs...)
        else
            dp = function (v,x,p,t)
                ForwardDiff.gradient(x->-H(v, x, p), x)
            end
            dq = function (v,x,p,t)
                ForwardDiff.gradient(v->H(v, x, p), v)
            end

            return ODEProblem(DynamicalODEFunction{T}(dp,dq), ArrayPartition(p0,q0), tspan, p; kwargs...)
        end
    else
        let cfg = ForwardDiff.GradientConfig(PhysicsTag(), p0), cfg2 = ForwardDiff.GradientConfig(PhysicsTag(), q0)
            dp = function (dv,v,x,p,t)
                fun2 = x->-H(v, x, p)
                ForwardDiff.gradient!(dv, fun2, x, cfg2, Val{false}())
            end
            dq = function (dx,v,x,p,t)
                fun1 = v-> H(v, x, p)
                ForwardDiff.gradient!(dx, fun1, v, cfg, Val{false}())
            end
            return ODEProblem(DynamicalODEFunction{T}(dp,dq), ArrayPartition(p0,q0), tspan, p; kwargs...)
        end
    end
end


text:
using LinearAlgebra

A = [2 3 0;
     1 2 4;
     0 2 1]

println("A+A=", A+A)
println("A-A=", A-A)
println("A*A=", A*A)
println("5*A=", 5*A)


text:
using Plots, LaTeXStrings
gr(colorbar=false,size=(600,150),legend=false,grid=false,xticks=false,yticks=false,axis=false)

using FourierGPE

## Initialize simulation
L = (60.0,)
N = (512,)
sim = Sim(L,N)
@unpack_Sim sim

## set simulation parameters
μ = 25.0
g = -0.01
γ = 0.0
Ns = 200
ξs = 2/abs(g)/Ns
us = 20
tf = π |> Float64
Nt = 150
ti = 0.0
t = LinRange(ti,tf,Nt)

## soliton wavefunction
x = X[1]
ψs(x) = sqrt(Ns/2ξs)*sech(x/ξs)*exp(im*us*x)
ψi = ψs.(x)
ϕi = kspace(ψi,sim);

## Set all fields
@pack_Sim! sim

## Evolve in k space
sol = runsim(sim)

## plot

y = abs2.(xspace(sol[Nt-76],sim))
plot(x,y,fill=(0, 0.2))

anim = @animate for i=1:Nt-8
    ψ = xspace(sol[i],sim)
    y = abs2.(ψ)
    plot(x,y,c=c3,fill=(0, 0.4,c3))
end

gif(anim, "./examples/brightsoliton.gif", fps = 25)


text:
module Descartes

import Base: *, union, diff, intersect
import GeometryBasics: AbstractMesh,
                       HyperRectangle,
                       Vec,
                       Mesh

using StaticArrays,
      Meshing,
      LinearAlgebra,
      CoordinateTransformations,
      Rotations

include("types.jl")
include("constructors.jl")
include("operations.jl")
include("frep.jl")
include("hyperrectangles.jl")
include("meshes.jl")

# 3d primitives
export Cuboid, Cylinder, Sphere, Piping

# 2d primitives
export Square, Circle

# transforms
export translate, rotate

# operations
export Shell, LinearExtrude

# CSG
export CSGUnion, RadiusedCSGUnion, CSGDiff, CSGIntersect

# frep
export FRep

end # module