text:
#=**************************************************************************

This script imports and runs an ozone NOx formaldehyde mechanism, producing 3
text files to be used in training physics constrained neural networks.

Written in Julia Version 1.2.0, tested up to 1.4.2

Written by Obin Sturm (posturm@ucdavis.edu)

The inverse problem tries to solve for S according to

  A*S = ∆C

 where S is the integrated reaction rates over a given timestep, and ∆C is the
 change in concentration.

The output files are "C.txt", "S.txt" and "delC.txt" for use in machine
learning research and inverse problems.  This was first used in a paper
introducing a framework for conserving mass in machine learning applications
for air quality and climate modeling: https://doi.org/10.5194/gmd-13-4435-2020

 C and delC correspond to the following species:
         1) O3
         2) NO
         3) NO2
         4) HCHO
         5) HO2
         6) HO2H
         7) HO.
         8) O
         9) HNO3
        10) CO
        11) H2

 S corresponds to the following reactions:
         1) NO2 + HV = NO + O
         2) O + O2 = O3
         3) O3 + NO = NO2 + O2
         4) HCHO + HV = 2 HO2. + CO
         5) HCHO + HV = H2 + CO
         6) HCHO + HO. = HO2. + CO + H2O
         7) HO2. + NO = HO. + NO2
         8) HO. + NO2 = HNO3
         9) HO2H + HV = 2 HO.
        10) HO2H + HO. = H2O + HO2
The stoichiometry matrix A, which can also be interpreted as the directed,
weighted incidence matrix of the species-reaction network, is

A =           [0  1 -1  0  0  0  0  0  0  0;
               1  0 -1  0  0  0 -1  0  0  0;
              -1  0  1  0  0  0  1 -1  0  0;
               0  0  0 -1 -1 -1  0  0  0  0;
               0  0  0  2  0  1 -1  0  0  1;
               0  0  0  0  0  0  0  0 -1 -1;
               0  0  0  0  0 -1  1 -1  2 -1;
               1 -1  0  0  0  0  0  0  0  0;
               0  0  0  0  0  0  0  1  0  0;
               0  0  0  1  1  1  0  0  0  0;
               0  0  0  0  1  0  0  0  0  0]
****************************************************************************=#

mechanism = include("assign3_driver.jl") # Include photochemical mechanism module
using DelimitedFiles

mechanism = include("assign3_driver.jl")
C,S,J = mechanism.gas(0.0,24.0*60.0,6.0,5000) # Run 5000 full days, sampling every 6 minutes
# delC = diff([C ; zeros(1,11)],dims = 1)

  open("S.txt", "w") do io
                         writedlm(io, S,',')
                   end
  open("C.txt", "w") do io
                          writedlm(io, C,',')
                   end
  open("J.txt", "w") do io
                         writedlm(io, J,',')
                   end


text:
include("eigenspectrum.jl")

function harmonicpotential{T<:Real}(x::T)
    kx = 1.0
    kx * (x^2) / 2
end

function harmonicpotential{T<:Real}(x::T, y::T)
    kx = 1.0
    ky = 1.0
    ( kx * (x^2) / 2 ) + (ky * (y^2) / 2 )
end

function harmonicpotential{T<:Real}(x::T, y::T, z::T)
    kx = 1.0
    ky = 1.0
    kz = 1.0
    ( kx * (x^2) / 2 ) + (ky * (y^2) / 2 ) + (kz * (z^2) / 2 )
end

function kineticenergy{T<:Real}(px::T)
    m = 1.0
    px^2 / (2*m)
end

function kineticenergy{T<:Real}(px::T, py::T)
    m = 1.0
    (px^2 + py^2) / (2*m)
end

function kineticenergy{T<:Real}(px::T, py::T, pz::T)
    m = 1.0
    (px^2 + py^2 + pz^2) / (2*m)
end

function solveschrodingereqn()
    pdesc = ProblemDescription(0.5, [4.0], [Int16(128)], 0.4, Int8(20))
    sysdesc = generatesystemdescriptor(pdesc)

    ψ = StateSet(Complex128, 20, 128)

    T = discretizeOperators(kineticenergy, sysdesc, pdesc, false)
    V = discretizeOperators(harmonicpotential, sysdesc, pdesc, true)

    # convergence constants
    ϵ = pdesc.timestep
    χ = 0.5
    γ = pdesc.tolerence

    # second order factorization
    eV = exp(-0.5*ϵ*V)
    eVˈ = exp(-0.5*ϵ*χ*V)
    eT = exp(-1.0*ϵ*T)
    eTˈ = exp(-1.0*ϵ*χ*T)
    τ = eV*eT*eV
    τˈ = eVˈ*eTˈ*eVˈ

    multiply!(τ, ψ, sysdesc)
    sysdesc.normenergies = log(orthonormalize!(ψ))/(-2.0*ϵ)

    H = expectationvalue(ψ, V)
    momentumspace!(ψ, sysdesc)
    H += expectationvalue(ψ, T)
    positionspace!(ψ, sysdesc)
    δH = 0.0
    δE = 0.0

    print_with_color(:green, "Iteration:")
    println(1)
    for h in H
        println(abs(h))
    end
    println("")

    iter = 1

    for i=1:100
        ψˈ = deepcopy(ψ)

        multiply!(τ, ψ, sysdesc)
        E = log(orthonormalize!(ψ))/(-2.0*ϵ)
        rmsE = norm(E)/sqrt(pdesc.nstates)
        δE = norm(E - sysdesc.normenergies)/sqrt(pdesc.nstates)

        multiply!(τˈ,ψˈ, sysdesc)
        Eˈ = log(orthonormalize!(ψˈ))/(-2.0*ϵ*χ)
        rmsEˈ = norm(Eˈ)/sqrt(pdesc.nstates)
        δEˈ = norm(Eˈ - sysdesc.normenergies)/sqrt(pdesc.nstates)

        H = expectationvalue(ψ, V)
        momentumspace!(ψ, sysdesc)
        H += expectationvalue(ψ, T)
        positionspace!(ψ, sysdesc)
        δH = norm(E - H)/sqrt(pdesc.nstates)

        Hˈ = expectationvalue(ψˈ, V)
        momentumspace!(ψˈ, sysdesc)
        Hˈ += expectationvalue(ψˈ, T)
        positionspace!(ψˈ, sysdesc)
        δHˈ = norm(Eˈ - Hˈ)/sqrt(pdesc.nstates)

        if rmsE<rmsEˈ || δE<δEˈ
            χ = sqrt(χ)
            sysdesc.normenergies = E
            sysdesc.energies = H
        else
            ϵ *= χ
            χ = χ^2
            ψ = ψˈ
            sysdesc.normenergies = Eˈ
            sysdesc.energies = Hˈ
            δE = δEˈ
            δH = δHˈ
        end

        if δE <= γ*δH
            iter = i+1
            break
        end

        eV = exp(-0.5*ϵ*V)
        eVˈ = exp(-0.5*ϵ*χ*V)
        eT = exp(-1.0*ϵ*T)
        eTˈ = exp(-1.0*ϵ*χ*T)
        τ = eV*eT*eV
        τˈ = eVˈ*eTˈ*eVˈ

        print_with_color(:green, "Iteration:")
        println(i+1)
        for h in H
            println(abs(h))
        end
        println("")

    end

    fs = open("output.txt", "w")
    serialize(fs, pdesc)
    serialize(fs, V)
    serialize(fs, ψ)
    serialize(fs, abs(H))
    serialize(fs, δE/δH)
    serialize(fs, iter)
    close(fs)
end

function discretegradient{T<:FloatingPoint}(operator::Function, syd::SystemDescription{T}, pdesc::ProblemDescription{T})
    grid = syd.pcoordinategrid
    cres = syd.cres
    t = [finitedifferencegradient(operator, p, cres) for p in grid]
    QuantumOperator(convert(Array{Complex{T}},reshape(t, convert(Array{Int64,1},pdesc.discretesteps)...)), true)
end

function finitedifferencegradient{T<:FloatingPoint}(operator::Function, point::Array{T,1}, step::Array{T,1})
    dim = length(point)
    grad = zeros(Complex{T}, dim)
    for i=1:dim
        k = zeros(step)
        k[i] = step[i]/2.0
        grad[i] = (operator((point+k)...) - operator((point-k)...))/(step[i])
    end
    norm(grad)^2
end

function solveschrodingereqn4()
    pdesc = ProblemDescription(0.5, [4.0], [Int16(64)], 0.05, Int8(4))
    sysdesc = generatesystemdescriptor(pdesc)

    ψ = StateSet(Complex128, 4, 64)

    T = discretizeOperators(kineticenergy, sysdesc, pdesc, false)
    V = discretizeOperators(harmonicpotential, sysdesc, pdesc, true)
    ΔV = discretegradient(harmonicpotential, sysdesc, pdesc)

    # convergence constants
    ϵ = pdesc.timestep
    χ = 0.5
    γ = pdesc.tolerence

    # second order factorization
    eV = exp(-(1.0/6.0)*ϵ*V)
    eVˈ = exp(-(1.0/6.0)*ϵ*χ*V)
    eT = exp(-0.5*ϵ*T)
    eTˈ = exp(-0.5*ϵ*χ*T)
    eΔV = exp(-(2.0/3.0)*ϵ*(V + ((1.0/48.0)*(ϵ^2)*ΔV)))
    eΔVˈ = exp(-(2.0/3.0)*ϵ*χ*(V + ((1.0/48.0)*((ϵ*χ)^2)*ΔV)))
    τ = eV*eT*eΔV*eT*eV
    τˈ = eVˈ*eTˈ*eΔVˈ*eTˈ*eVˈ

    multiply!(τ, ψ, sysdesc)
    sysdesc.normenergies = log(orthonormalize!(ψ))/(-2.0*ϵ)

    H = expectationvalue(ψ, V)
    momentumspace!(ψ, sysdesc)
    H += expectationvalue(ψ, T)
    positionspace!(ψ, sysdesc)
    δH = 0.0
    δE = 0.0

    print_with_color(:green, "Iteration:")
    println(1)
    for h in H
        println(abs(h))
    end
    println("")

    iter = 1
    for i=1:100
        ψˈ = deepcopy(ψ)

        multiply!(τ, ψ, sysdesc)
        E = log(orthonormalize!(ψ))/(-2.0*ϵ)
        rmsE = norm(E)/sqrt(pdesc.nstates)
        δE = norm(E - sysdesc.normenergies)/sqrt(pdesc.nstates)

        multiply!(τˈ,ψˈ, sysdesc)
        Eˈ = log(orthonormalize!(ψˈ))/(-2.0*ϵ*χ)
        rmsEˈ = norm(Eˈ)/sqrt(pdesc.nstates)
        δEˈ = norm(Eˈ - sysdesc.normenergies)/sqrt(pdesc.nstates)

        H = expectationvalue(ψ, V)
        momentumspace!(ψ, sysdesc)
        H += expectationvalue(ψ, T)
        positionspace!(ψ, sysdesc)
        δH = norm(E - H)/sqrt(pdesc.nstates)

        Hˈ = expectationvalue(ψˈ, V)
        momentumspace!(ψˈ, sysdesc)
        Hˈ += expectationvalue(ψˈ, T)
        positionspace!(ψˈ, sysdesc)
        δHˈ = norm(Eˈ - Hˈ)/sqrt(pdesc.nstates)

        if rmsE<rmsEˈ || δE<δEˈ
            χ = sqrt(χ)
            sysdesc.normenergies = E
            sysdesc.energies = H
        else
            ϵ *= χ
            χ = χ^2
            ψ = ψˈ
            sysdesc.normenergies = Eˈ
            sysdesc.energies = Hˈ
            δE = δEˈ
            δH = δHˈ
        end

        if δE <= γ*δH
            iter = i
            break
        end

        eV = exp(-(1.0/6.0)*ϵ*V)
        eVˈ = exp(-(1.0/6.0)*ϵ*χ*V)
        eT = exp(-0.5*ϵ*T)
        eTˈ = exp(-0.5*ϵ*χ*T)
        eΔV = exp(-(2.0/3.0)*ϵ*(V + ((1.0/48.0)*(ϵ^2)*ΔV)))
        eΔVˈ = exp(-(2.0/3.0)*ϵ*χ*(V + ((1.0/48.0)*((ϵ*χ)^2)*ΔV)))
        τ = eV*eT*eΔV*eT*eV
        τˈ = eVˈ*eTˈ*eΔVˈ*eTˈ*eVˈ

        print_with_color(:green, "Iteration:")
        println(i+1)
        for h in H
            println(abs(h))
        end
        println("")

    end

    fs = open("output.txt", "w")
    serialize(fs, pdesc)
    serialize(fs, V)
    serialize(fs, ψ)
    serialize(fs, abs(H))
    serialize(fs, δE/δH)
    serialize(fs, iter)
    close(fs)
end

@time solveschrodingereqn()

ext: jl


text:
module Groups

using GroupsCore
using ThreadsX
import KnuthBendix
import KnuthBendix: AbstractWord, Alphabet, Word
import KnuthBendix: alphabet
import Random

import OrderedCollections: OrderedSet

export Alphabet, AutomorphismGroup, FreeGroup, FreeGroup, FPGroup, FPGroupElement, SpecialAutomorphismGroup
export alphabet, evaluate, word, gens

include("types.jl")
include("hashing.jl")
include("normalform.jl")
include("autgroups.jl")

include("groups/sautFn.jl")
include("groups/mcg.jl")

include("wl_ball.jl")
end # of module Groups


text:
"""
$(TYPEDEF)

Represents a difference operator.

# Fields
$(FIELDS)

# Examples

```jldoctest
julia> using Symbolics

julia> @variables t;

julia> D = Difference(t; dt=0.01)
Difference(t; dt=0.01)
```
"""
struct Difference <: Function
    """Fixed Difference"""
    t
    dt
    Difference(t; dt) = new(value(t), dt)
end
(D::Difference)(t) = Term{symtype(t)}(D, [t])
(D::Difference)(t::Num) = Num(D(value(t)))
SymbolicUtils.promote_symtype(::Difference, t) = t

Base.show(io::IO, D::Difference) = print(io, "Difference(", D.t, "; dt=", D.dt, ")")

Base.:(==)(D1::Difference, D2::Difference) = isequal(D1.t, D2.t) && isequal(D1.dt, D2.dt)
Base.isequal(D1::Difference, D2::Difference) = isequal(D1.t, D2.t) && isequal(D1.dt, D2.dt)
Base.hash(D::Difference, u::UInt) = hash(D.dt, hash(D.t, xor(u, 0x055640d6d952f101)))


text:
using CSV
using DataFrames
using DataStructures
using JSON
const Path = Vector{String}

function parserange(str::AbstractString)
  parts = split(str, ":")
  if 2 == length(parts)
    b = parse(Float64, parts[1])
    e = parse(Float64, parts[2])
    b:e
  elseif 3 == length(parts)
    b = parse(Float64, parts[1])
    s = parse(Float64, parts[2])
    e = parse(Float64, parts[3])
    b:s:e
  else
    nothing
  end
end

function findranges(dict::AbstractDict{String,T} where T)
  paths = Vector{Path}()
  for key in keys(dict)
    val = dict[key]
    if isa(val, AbstractString)
      if nothing != parserange(val)
        push!(paths, Path([key]))
      end
    elseif isa(val, AbstractDict)
      append!(paths, map(path->insert!(copy(path), 1, key), findranges(dict[key])))
    end
  end
  paths
end

getbypath(dict::AbstractDict{T} where T<:AbstractString, path::Path) = getbypath(dict[path[1]], path[2:end])
getbypath(val, ::Path) = val

function setbypath!(dict::AbstractDict{T} where T<:AbstractString, path::Path, value::Real)
  if length(path) > 1
    setbypath!(dict[path[1]], path[2:end], value)
  else
    dict[only(path)] = value
  end
end

function main()
  @assert length(ARGS) > 0 "JSON file needed"

  json = JSON.parsefile(ARGS[1], dicttype=OrderedDict)
  workdir = length(ARGS)>1 ? ARGS[2] : splitext(ARGS[1])[1]

  rangepaths = findranges(json) |> sort
  ranges = map(x->getbypath(json, x) |> parserange, rangepaths)

  @info "ranges are" rangepaths ranges

  @assert length(ranges) == 2 "we support grids now only"
  xpath, ypath = rangepaths
  xrange, yrange = ranges

  df = DataFrame([String[], Float64[], Float64[]], [:path, Symbol(join(xpath, "_")), Symbol(join(ypath, "_"))])

  for (x,xval) in enumerate(xrange)
    for (y,yval) in enumerate(yrange)
      subdir = joinpath(workdir, "grid_$(x-1)_$(y-1)")
      mkpath(subdir)
      specific_json = deepcopy(json)
      setbypath!(specific_json, xpath, xval)
      setbypath!(specific_json, ypath, yval)

      open(joinpath(subdir, "params_experiment.json"),"w") do f
        JSON.print(f, specific_json)
      end
      push!(df, ("grid_$(x-1)_$(y-1)", xval, yval))
    end
  end
  CSV.write(joinpath(workdir, "parameters_map.csv"), df)

  open(joinpath(workdir, "template.json"),"w") do f
    JSON.print(f, json, 2)
  end

end

main()


text:
#!/usr/bin/env julia
#
# @license Apache-2.0
#
# Copyright (c) 2022 The Stdlib Authors.
#
# 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.

import BenchmarkTools

# Benchmark variables:
repeats = 3;
samples = 1e6;
count = 0;

"""
	print_version()

Prints the TAP version.

# Examples

``` julia
julia> print_version()
```
"""
function print_version()
	@printf( "TAP version 13\n" );
end

"""
	print_summary( total, passing )

Print the benchmark summary.

# Arguments

* `total`: total number of tests
* `passing`: number of passing tests

# Examples

``` julia
julia> print_summary( 3, 3 )
```
"""
function print_summary( total, passing )
	@printf( "#\n" );
	@printf( "1..%d\n", total ); # TAP plan
	@printf( "# total %d\n", total );
	@printf( "# pass  %d\n", passing );
	@printf( "#\n" );
	@printf( "# ok\n" );
end

"""
	print_results( iterations, elapsed )

Print benchmark results.

# Arguments

* `iterations`: number of iterations
* `elapsed`: elapsed time (in seconds)

# Examples

``` julia
julia> print_results( 1000000, 0.131009101868 )
```
"""
function print_results( iterations, elapsed )
	rate = iterations / elapsed

	@printf( "  ---\n" );
	@printf( "  iterations: %d\n", iterations );
	@printf( "  elapsed: %0.9f\n", elapsed );
	@printf( "  rate: %0.9f\n", rate );
	@printf( "  ...\n" );
end

"""
	benchmark( setup, expr )

Run a benchmark.

# Notes

* Benchmark results are returned as a two-element array: [ iterations, elapsed ].
* The number of iterations is not the true number of iterations. Instead, an 'iteration' is defined as a 'sample', which is a computed estimate for a single evaluation.
* The elapsed time is in seconds.

# Arguments

* `setup`: setup expression
* `expr`: expression to benchmark

# Examples

``` julia
julia> out = benchmark( :(), :( sin( 3.14 ) ) );
```
"""
function benchmark( setup, expr )
	t = eval( :( BenchmarkTools.@benchmark $expr samples=$samples setup=($setup) ) )

	# Compute the total "elapsed" time and convert from nanoseconds to seconds:
	s = sum( t.times ) / 1.0e9;

	# Determine the number of "iterations":
	iter = length( t.times );

	# Return the results:
	[ iter, s ];
end

"""
	bench( name, setup, expr )

Run a named benchmark.

# Arguments

* `name`: benchmark name (suffix)
* `setup`: setup expression
* `expr`: expression to benchmark

# Examples

``` julia
julia> bench( "sin", :(), :( sin( 3.14 ) ) );
```
"""
function bench( name, setup, expr )
	for i in 1:repeats
		@printf( "# julia::%s\n", name );
		global count += 1;
		results = benchmark( setup, expr );
		print_results( results[ 1 ], results[ 2 ] );
		@printf( "ok %d benchmark finished\n", count );
	end
end

"""
	main()

Run benchmarks.

# Examples

``` julia
julia> main();
```
"""
function main()
	print_version();

	name = "TODO";
	setup = :();
	stmt = :( TODO );
	bench( name, setup, stmt );

	print_summary( count, count );
end

main();

max_stars_repo_name: mhmdaminraeisi/stdlib

max_issues_repo_name: mhmdaminraeisi/stdlib

max_forks_repo_name: mhmdaminraeisi/stdlib


text:
using Test, Mods

sample = joinpath(@__DIR__, "sample")
example = readlines(sample)

input = joinpath(@__DIR__, "input")
lines = readlines(input)

function prepare(input)
    now = parse(Int, input[1])
    schedule = map(id -> id == "x" ? missing : parse(Int, id), split(input[2], ","))
    return (now, schedule)
end

function busses(schedule)
    return [(id, idx - 1) for (idx, id) in enumerate(schedule) if !ismissing(id)]
end

function ids(busses)
    return first.(busses)
end

function offsets(busses)
    return last.(busses)
end

function depart(now, schedule)
    wait, id = minimum(map(id -> (mod(-1 * now, id), id), filter(!ismissing, schedule)))
    return id * wait
end

function solver(schedule)
    targets = busses(schedule)
    (period, t) = (1, 0)
    for idx in 1:length(targets)
        (candidate, offset) = targets[idx]
        while ((t + offset) % candidate != 0)
            t += period
        end
        period *= candidate
    end
    return t
end

(now, schedule) = prepare(example)

@test depart(now, schedule) == 295
@test solver(schedule) == 1068781

(now, schedule) = prepare(lines)

p1 = depart(now, schedule)
p2 = solver(schedule)

@assert(p1 == 104)
@assert(p2 == 842186186521918)

println("-----------------------------------------------------------------------")
println("shuttle search -- part one :: $p1")
println("shuttle search -- part two :: $p2")
println("-----------------------------------------------------------------------")


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

using Plots;

"""
    democoons()

Construction of a bilinearly blended Coons surface.

# Examples:
```julia
julia> democoons()
```
"""
function democoons()
# Boundary curve 1
pnts = [0.0  3.0  4.5  6.5 8.0 10.0;
        0.0  0.0  0.0  0.0 0.0  0.0;
        2.0  2.0  7.0  4.0 7.0  9.0];
crv1 = nrbmak(pnts, vec([0 0 0 1/3 0.5 2/3 1 1 1]));

# Boundary curve 2
pnts= [ 0.0  3.0  5.0  8.0 10.0;
       10.0 10.0 10.0 10.0 10.0;
        3.0  5.0  8.0  6.0 10.0];
crv2 = nrbmak(pnts, vec([0 0 0 1/3 2/3 1 1 1]));

# Boundary curve 3
pnts= [0.0 0.0 0.0 0.0;
       0.0 3.0 8.0 10.0;
       2.0 0.0 5.0 3.0];
crv3 = nrbmak(pnts, vec([0 0 0 0.5 1 1 1]));

# Boundary curve 4
pnts= [10.0 10.0 10.0 10.0 10.0;
       0.0   3.0  5.0  8.0 10.0;
       9.0   7.0  7.0 10.0 10.0];
crv4 = nrbmak(pnts, vec([0 0 0 0.25 0.75 1 1 1]));

srf = nrbcoons(crv1, crv2, crv3, crv4);

# Draw the surface
nrbplot(srf,[20;20],
        title = "Construction of a bilinearly blended Coons surface.",
        c = :roma,
        legend = nothing,
        linewidth = 0.5,
        linecolor = :black,
        camera = (-60,30));
end # democoons


text:
using PtFEM, Test

data = Dict(
  # Plane(ndim, nst, nxe, nye, nip, direction, finite_element(nod, nodof), axisymmetric)
  :struc_el => Plane(2, 3, 3, 2, 4, :y, Quadrilateral(4, 2), false),
  :properties => [1.0e6 0.3;],
  :x_coords => [0.0, 10.0, 20.0, 30.0],
  :y_coords => [0.0, -5.0, -10.0],
  :support => [
    (1, [0 1]),
    (2, [0 1]),
    (3, [0 0]),
    (6, [0 0]),
    (9, [0 0]),
    (10, [0 1]),
    (11, [0 1]),
    (12, [0 0])
    ],
  :fixed_freedoms => [
    (1, 2, -1.0e-5),
    (4, 2, -1.0e-5)
    ]
)

@time m, dis_df, fm_df = p51(data)

@test m.loads[0:12] ≈ [0.0,-1.0000000000000035e-5,-5.152429570719046e-6,
  8.101475481120316e-8,-1.0000000000000072e-5,1.5820948250551495e-6,
  -4.593608962323964e-6,1.2407010196358915e-7,1.257758473173503e-6,
  1.4721403994132775e-6,1.953453349357416e-7,2.8154845809451437e-7,
  3.4748952771616183e-7]


text:
# v0.6

using Primes

function mersennefactor(p::Int)::Int
    q = 2p + 1
    while true
        if log2(q) > p / 2
            return -1
        elseif q % 8 in (1, 7) && Primes.isprime(q) && powermod(2, p, q) == 1
            return q
        end
    q += 2p
    end
end

for i in filter(Primes.isprime, push!(collect(1:20), 929))
    mf = mersennefactor(i)
    if mf != -1 println("M$i = ", mf, " × ", (big(2) ^ i - 1) ÷ mf)
    else println("M$i is prime") end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Underscores
using Test

@testset "Underscores Examples" begin
    @test [1,"a",2.0] == @_ map(_, [1,"a",2.0])

    strs = ["ab", "ca", "ad"]
    @test [true, false, true] == @_ map(startswith(_, "a"), strs)

    data = [(x="a", y=1),
            (x="b", y=2),
            (x="c", y=3)]

    @test data[1:1] == @_ filter(startswith(_.x, "a"), data)
    @test data[2:3] == @_ filter(_.y >= 2, data)

    # Multiple args
    @test [0,0] == @_ map(_-_, [1,2])

    # Use with piping and __
    @test [1] == @_ data |>
                    filter(startswith(_.x, "a"), __) |>
                    map(_.y, __)

    @test [3,2,1] == @_ data |>
                        sort(__, by=_.x, rev=true) |>
                        map(_.y, __)

    @test [1] == (@_ map(_.y, __) ∘
                     filter(startswith(_.x, "a"), __))(data)

    @test [0,0,0] == @_ data |> map(_1.y + _2, __, [-1,-2,-3])

    # Use with piping and lazy versions of map and filter
    Filter(f) = x->filter(f,x)
    Map(f) = x->map(f,x)

    @test [1] == @_ data |>
                    Filter(startswith(_.x, "a")) |>
                    Map(_.y)

    @test [1] == @_(Map(_.y) ∘ Filter(startswith(_.x, "a")))(data)
end

@testset "Underscores lowering" begin
    cleanup! = Base.remove_linenums!
    lower(ex) = cleanup!(Underscores.lower_underscores(ex))

    # Simple cases
    # _
    @test lower(:(f(_))) == cleanup!(:(f((_1,)->_1)))
    @test lower(:(f(g(h(_))))) == cleanup!(:(f(((_1,)->g(h(_1))))))
    # __
    @test lower(:(f(__))) == cleanup!(:((__1,)->f(__1)))
    @test lower(:(f(g(h(__))))) == cleanup!(:((__1,)->f(g(h(__1)))))
    # Tricky case when function name is itself an underscore expression.
    # We want the function name to be also expanded here.
    @test lower(:(g(_)(x))) == cleanup!(:(((_1,)->g(_1))(x)))

    # Multiple arguments
    # _
    @test lower(:(f(_,a))) == cleanup!(:(f(((_1,)->_1), a)))
    @test lower(:(f(a,_))) == cleanup!(:(f(a, ((_1,)->_1))))
    @test lower(:(f(_,_))) == cleanup!(:(f(((_1,)->_1), ((_1,)->_1))))
    @test lower(:(f(g(_,_)))) == cleanup!(:(f(((_1,)->g(_1,_1)))))
    # __
    @test lower(:(f(__,__))) == cleanup!(:((__1,)->f(__1, __1)))

    # Numbered arguments
    # _
    @test lower(:(f(_1))) == cleanup!(:(f((_1,)->_1)))
    @test lower(:(f(_2))) == cleanup!(:(f((_1,_2)->_2)))
    @test lower(:(f(_2+_1))) == cleanup!(:(f((_1,_2)->_2+_1)))
    # __
    @test lower(:(f(__1))) == cleanup!(:((__1,)->f(__1)))
    @test lower(:(f(__2))) == cleanup!(:((__1,__2)->f(__2)))
    @test lower(:(f(__2+__1))) == cleanup!(:((__1,__2)->f(__2+__1)))

    # Cute subscript-numbered arguments
    @test lower(:(f(_₁))) == cleanup!(:(f((_1,)->_1)))
    @test lower(:(f(_₂))) == cleanup!(:(f((_1,_2)->_2)))
    @test lower(:(f(__₁))) == cleanup!(:((__1,)->f(__1)))

    # Can't mix numbered and non-numbered placeholders
    @test_throws ArgumentError lower(:(f(_+_1)))
    @test_throws ArgumentError lower(:(f(__+__1)))

    # piping and composition
    # _
    @test lower(:(f(_) |> g(_) |> h(_))) ==
          cleanup!(:(f((_1,)->_1) |>
                     g((_1,)->_1) |>
                     h((_1,)->_1)))
    @test lower(:(f(_) ∘ g(_) ∘ h(_))) ==
          cleanup!(:(f((_1,)->_1) ∘
                     g((_1,)->_1) ∘
                     h((_1,)->_1)))
    @test lower(:(f(_) <| g(_) <| h(_))) ==
          cleanup!(:(f((_1,)->_1) <|
                     g((_1,)->_1) <|
                     h((_1,)->_1)))
    # __
    @test lower(:(f(_, __) |> g(_, __))) ==
        cleanup!(:(((__1,)->f((_1,)->_1, __1)) |>
                   ((__1,)->g((_1,)->_1, __1))))

    # Keyword arguments
    @test lower(:(f(x, k=_+1))) == cleanup!(:(f(x, k=((_1,)->_1+1))))
    @test lower(:(f(x; k=_+1))) == cleanup!(:(f(x; k=((_1,)->_1+1))))

    # Broadcast notation
    @test lower(:(f.(_)))   == cleanup!(:(f.((_1,)->_1)))
    @test lower(:(f.(_,_))) == cleanup!(:(f.((_1,)->_1, (_1,)->_1)))
    @test lower(:(f.(__)))  == cleanup!(:((__1,)->f.(__1)))
    @test lower(:((_).(x)))  == cleanup!(:(((_1,)->_1).(x)))

    # Random sample of other syntax
    @test lower(:([_]))  == cleanup!(:([(_1,)->_1]))
    @test lower(:((f(_), g(_))))  == cleanup!(:(((_1,)->f(_1)), ((_1,)->g(_1))))
    @test lower(:([__])) == cleanup!(:((__1,)->[__1]))
    @test lower(:(__.x)) == cleanup!(:((__1,)->__1.x))

    # do syntax is disabled for now because desired behaviour is not entirely
    # clear. See #4
    @test_throws ErrorException lower(:(f() do
                                            body
                                        end))
end


text:
using CondaPkg
using Test

# Only run the gc tests on CI (because it's annoying to do it locally)
testgc = get(ENV, "CI", "") == "true"

status() = sprint(io -> CondaPkg.status(io=io))

@testset "CondaPkg" begin

    @testset "add/remove channel" begin
        @test !occursin("conda-forge", status())
        CondaPkg.add_channel("conda-forge")
        @test occursin("conda-forge", status())
        CondaPkg.rm_channel("conda-forge")
        @test !occursin("conda-forge", status())
    end

    @testset "install python" begin
        @test !occursin("python", status())
        CondaPkg.add("python", version="==3.10.2")
        CondaPkg.resolve()
        CondaPkg.withenv() do
            pythonpath = joinpath(CondaPkg.envdir(), Sys.iswindows() ? "python.exe" : "bin/python")
            @test isfile(pythonpath)
        end
        @test occursin("python", status())
    end

    @testset "install/remove python package" begin
        # verify package isn't already installed
        @test !occursin("six", status())
        CondaPkg.withenv() do
            @test_throws Exception run(`python -c "import six"`)
        end

        # install package
        CondaPkg.add("six", version="==1.16.0")
        @test occursin("six", status())
        @test occursin("(==1.16.0)", status())
        CondaPkg.withenv() do
            run(`python -c "import six"`)
        end
        @test occursin("v1.16.0", status())

        # remove package
        CondaPkg.rm("six")
        @test !occursin("six", status())
        CondaPkg.withenv() do
            @test_throws Exception run(`python -c "import six"`)
        end
    end

    @testset "pip install/remove python package" begin
        # verify package isn't already installed
        @test !occursin("six", status())
        CondaPkg.withenv() do
            @test_throws Exception run(`python -c "import six"`)
        end

        # install package
        CondaPkg.add_pip("six", version="==1.16.0")
        @test occursin("six", status())
        @test occursin("(==1.16.0)", status())
        CondaPkg.withenv() do
            run(`python -c "import six"`)
        end
        @test occursin("v1.16.0", status())

        # remove package
        CondaPkg.rm_pip("six")
        @test !occursin("six", status())
        CondaPkg.withenv() do
            @test_throws Exception run(`python -c "import six"`)
        end
    end

    @testset "install/remove executable package" begin
        # install package, verify that executable exists
        CondaPkg.add("curl")
        curl_path = CondaPkg.which("curl")
        @test isfile(curl_path)

        # test uninstall
        CondaPkg.rm("curl")
        CondaPkg.resolve(force=true)
        @test !isfile(curl_path)
    end

    if testgc
        @testset "gc()" begin
            # verify that micromamba clean runs without errors
            CondaPkg.gc()
        end
    end

    @testset "validation" begin
        @test_throws Exception CondaPkg.add("!invalid!package!")
        @test_throws Exception CondaPkg.add_pip("!invalid!package!")
        @test_throws Exception CondaPkg.add("valid-package", version="*invalid*version*")
        @test_throws Exception CondaPkg.add_pip("valid-package", version="*invalid*version*")
        @test_throws Exception CondaPkg.add_channel("")
        @test !occursin("valid", status())
    end

    @testset "PkgREPL" begin
        # The functions in CondaPkg.PkgREPL are not API
        # but calling them is the closest thing to calling
        # the Pkg REPL.

        # add
        CondaPkg.PkgREPL.add([" foo >=1.2 "])
        @test occursin("foo", status())
        @test occursin("(>=1.2)", status())

        # rm
        CondaPkg.PkgREPL.rm([" foo "])
        @test !occursin("foo", status())

        # add channel
        CondaPkg.PkgREPL.add([" foo-channel "], mode=:channel)
        @test occursin("foo-channel", status())

        # rm channel
        CondaPkg.PkgREPL.rm([" foo-channel "], mode=:channel)
        @test !occursin("foo-channel", status())

        # add pip
        CondaPkg.PkgREPL.add([" foo ~=1.3 "], mode=:pip)
        @test occursin("foo", status())
        @test occursin("(~=1.3)", status())

        # rm pip
        CondaPkg.PkgREPL.rm([" foo "], mode=:pip)
        @test !occursin("foo", status())

        # status
        # TODO: capture the output and check it equals status()
        CondaPkg.PkgREPL.status()

        # resolve
        CondaPkg.PkgREPL.resolve()

        # gc
        testgc && CondaPkg.PkgREPL.gc()

        # run
        # TODO: capture the output and check it contains "Python 3.10.2"
        CondaPkg.PkgREPL.run(["python", "--version"])
    end

end


text:
using CompScienceMeshes
using StaticArrays
using Test

for T in [Float32, Float64]
    local p1 = point(T, 1.0, 0.0, 0.0)
    local p2 = point(T, 10.0, 0.0, 1.0)
    local p3 = point(T, 0.0, 2.0, 5.0)

    ch = simplex(p1,p2,p3)

    dom = domain(ch)
    q, u = quadpoints(dom, 7)
    @which neighborhood(ch, q[1])

    pw = CompScienceMeshes.quadpoints(ch, 7);

    #t1 = sum(w)
    t1 = sum(q[2] for q in pw)
    t2 = T(0.5)*norm(cross(p1-p3, p2-p3))
    @test 2 * (t1-t2) / (t1+t2) + 1  ≈ 1

    p1 = point(T, 1.0, 0.0, 0.0)
    p2 = point(T, 2.0, 0.0, 0.0)
    seg = simplex(p1,p2)

    wpvs = quadpoints(x->T(x[1]^2), [seg], (20,))
    global r = T(0.0)
    for wpv in wpvs[1,1]
        w = wpv.weight
        v = wpv.value
        r += w * v
    end

    @test r ≈ T(1/3*(2^3-1^3))
end

max_stars_repo_name: krcools/CompScienceMeshes

max_issues_repo_name: krcools/CompScienceMeshes

max_forks_repo_name: krcools/CompScienceMeshes


text:
using Test
using Sqlite3Stats
using DataFrames
using Sqlite3Stats.SQLite

@testset "Functions" begin
    db = SQLite.DB()

    SQLite.execute(db, "create table Numbers (val float, otherval float)")
    x = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
    for v in x
        y = 11.0 - v
        SQLite.execute(db, "insert into Numbers (val, otherval) values ($(v), $(y))")
    end

    result = DBInterface.execute(db, "select val, otherval from Numbers") |> DataFrame

    Sqlite3Stats.register_functions(db)

    result = DBInterface.execute(db, "select val, otherval from Numbers") |> DataFrame
    @test size(result) == (10, 2)

    @testset "Q1" begin
        result =
            DBInterface.execute(db, "select count(val), Q1(val) as MYQ1 from Numbers") |>
            DataFrame
        @test result[!, "MYQ1"] == [3.25]
    end

    @testset "Q2" begin
        result =
            DBInterface.execute(db, "select count(val), Q2(val) as MYQ2 from Numbers") |>
            DataFrame
        @test result[!, "MYQ2"] == [5.5]
    end

    @testset "MEDIAN" begin
        result =
            DBInterface.execute(
                db,
                "select count(val), MEDIAN(val) as MYMEDIAN from Numbers",
            ) |> DataFrame
        @test result[!, "MYMEDIAN"] == [5.5]
    end

    @testset "Q3" begin
        result =
            DBInterface.execute(db, "select count(val), Q3(val) as MYQ3 from Numbers") |>
            DataFrame
        @test result[!, "MYQ3"] == [7.75]
    end

    @testset "QUANTILE" begin
        result =
            DBInterface.execute(
                db,
                "select QUANTILE(val, 0.25) as MYRESULT from Numbers",
            ) |> DataFrame
        @test result[!, "MYRESULT"] == [3.25]

        result =
            DBInterface.execute(
                db,
                "select QUANTILE(val, 0.50) as MYRESULT from Numbers",
            ) |> DataFrame
        @test result[!, "MYRESULT"] == [5.5]

        result =
            DBInterface.execute(
                db,
                "select QUANTILE(val, 0.75) as MYRESULT from Numbers",
            ) |> DataFrame
        @test result[!, "MYRESULT"] == [7.75]
    end

    @testset "COV" begin
        result =
            DBInterface.execute(db, "select COV(val, val) as MYCOV from Numbers") |>
            DataFrame
        @test result[!, "MYCOV"] == [9.166666666666666]
    end

    @testset "COR" begin
        result =
            DBInterface.execute(db, "select COR(val, val) as MYCOR from Numbers") |>
            DataFrame
        @test result[!, "MYCOR"] == [1.0]

        result =
            DBInterface.execute(db, "select COR(val, otherval) as MYCOR from Numbers") |>
            DataFrame
        @test result[!, "MYCOR"] == [-1.0]
    end

    @testset "CORSPEARMAN" begin
        result =
            DBInterface.execute(db, "select CORSPEARMAN(val, val) as MYCOR from Numbers") |>
            DataFrame
        @test result[!, "MYCOR"] == [1.0]

        result =
            DBInterface.execute(
                db,
                "select CORSPEARMAN(val, otherval) as MYCOR from Numbers",
            ) |> DataFrame
        @test result[!, "MYCOR"] == [-1.0]
    end

    @testset "CORKENDALL" begin
        result =
            DBInterface.execute(db, "select CORKENDALL(val, val) as MYCOR from Numbers") |>
            DataFrame
        @test result[!, "MYCOR"] == [1.0]

        result =
            DBInterface.execute(
                db,
                "select CORKENDALL(val, otherval) as MYCOR from Numbers",
            ) |> DataFrame
        @test result[!, "MYCOR"] == [-1.0]
    end

    @testset "MAD" begin
        result =
            DBInterface.execute(db, "select MAD(val) as MYMAD from Numbers") |> DataFrame
        @test result[!, "MYMAD"] == [3.7065055462640046]
    end

    @testset "IQR" begin
        result =
            DBInterface.execute(db, "select IQR(val) as MYIQR from Numbers") |> DataFrame
        @test result[!, "MYIQR"] == [4.5]
    end

    @testset "SKEWNESS" begin
        result =
            DBInterface.execute(db, "select SKEWNESS(val) as MYSKEW from Numbers") |>
            DataFrame
        @test result[!, "MYSKEW"] == [0.0]
    end

    @testset "KURTOSIS" begin
        result =
            DBInterface.execute(db, "select KURTOSIS(val) as MYKURT from Numbers") |>
            DataFrame
        @test result[!, "MYKURT"] == [-1.2242424242424244]
    end

    @testset "GEOMEAN" begin
        result =
            DBInterface.execute(db, "select GEOMEAN(val) as MYGEOMEAN from Numbers") |>
            DataFrame
        @test result[!, "MYGEOMEAN"] == [4.528728688116766]
    end

    @testset "HARMMEAN" begin
        result =
            DBInterface.execute(db, "select HARMMEAN(val) as MYHARMMEAN from Numbers") |>
            DataFrame
        @test result[!, "MYHARMMEAN"] == [3.414171521474055]
    end

    @testset "MODE" begin
        result = DBInterface.execute(db, "select MODE(val) as m from Numbers") |> DataFrame
        @test result[!, "m"] == [1.0]

        result =
            DBInterface.execute(db, "select MODE(otherval) as m from Numbers") |> DataFrame
        @test result[!, "m"] == [10.0]
    end

    @testset "MAXAD" begin
        result =
            DBInterface.execute(db, "select MAXAD(val, val) as m from Numbers") |> DataFrame
        @test result[!, "m"] == [0.0]

        result =
            DBInterface.execute(db, "select MAXAD(val, otherval) as m from Numbers") |>
            DataFrame
        @test result[!, "m"] == [9.0]
    end

    @testset "MEANAD" begin
        result =
            DBInterface.execute(db, "select MEANAD(val, val) as m from Numbers") |>
            DataFrame
        @test result[!, "m"] == [0.0]

        result =
            DBInterface.execute(db, "select MEANAD(val, otherval) as m from Numbers") |>
            DataFrame
        @test result[!, "m"] == [5.0]
    end

    @testset "MSD" begin
        result =
            DBInterface.execute(db, "select MSD(val, val) as m from Numbers") |> DataFrame
        @test result[!, "m"] == [0.0]

        result =
            DBInterface.execute(db, "select MSD(val, otherval) as m from Numbers") |>
            DataFrame
        @test result[!, "m"] == [33.0]
    end

    SQLite.close(db)
end

@testset "Weighted Functions" begin
    db = SQLite.DB()

    SQLite.execute(db, "create table Numbers (val float, w float)")
    x = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
    w = [
        0.01818181818181818,
        0.03636363636363636,
        0.05454545454545454,
        0.07272727272727272,
        0.09090909090909091,
        0.10909090909090909,
        0.12727272727272726,
        0.14545454545454545,
        0.16363636363636364,
        0.18181818181818182,
    ]
    for i = 1:10
        v = x[i]
        y = w[i]
        SQLite.execute(db, "insert into Numbers (val, w) values ($(v), $(y))")
    end

    result = DBInterface.execute(db, "select val, w from Numbers") |> DataFrame

    Sqlite3Stats.register_functions(db)

    result = DBInterface.execute(db, "select val, w from Numbers") |> DataFrame
    @test size(result) == (10, 2)

    @testset "WMEAN" begin
        result =
            DBInterface.execute(db, "select WMEAN(val, w) as MYRESULT from Numbers") |>
            DataFrame
        @test result[!, "MYRESULT"] == [7.0]
    end

    @testset "WMEDIAN" begin
        result =
            DBInterface.execute(db, "select WMEDIAN(val, w) as MYRESULT from Numbers") |>
            DataFrame
        @test result[!, "MYRESULT"] == [7.0]
    end

    @testset "ENTROPY" begin
        result =
            DBInterface.execute(db, "select ENTROPY(w) as MYRESULT from Numbers") |>
            DataFrame
        @test result[!, "MYRESULT"] == [2.151281720651836]
    end

    SQLite.close(db)

end

@testset "Linear Regression" begin

    db = SQLite.DB()

    SQLite.execute(db, "create table Numbers (x float, y float)")
    x = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
    y = 2.0 .* x
    for i = 1:10
        u = x[i]
        v = y[i]
        SQLite.execute(db, "insert into Numbers (x, y) values ($(u), $(v))")
    end

    Sqlite3Stats.register_functions(db)

    @testset "Raw function" begin
        x = [1.0, 2.0, 3.0, 4.0, 5.0]
        y = [2.0, 4.0, 6.0, 8.0, 10.0]
        reg = Sqlite3Stats.linear_regression(x, y)
        @test reg[1] == 0.0
        @test reg[2] == 2.0
    end

    @testset "LINSLOPE" begin
        result =
            DBInterface.execute(db, "select LINSLOPE(x, y) as MYRESULT from Numbers") |>
            DataFrame
        @test result[!, "MYRESULT"] == [2.0]
    end

    @testset "LININTERCEPT" begin
        result =
            DBInterface.execute(db, "select LININTERCEPT(x, y) as MYRESULT from Numbers") |>
            DataFrame
        @test result[!, "MYRESULT"] == [0.0]
    end

    SQLite.close(db)
end

@testset "Normal Distribution" begin

    tol = 0.001

    db = SQLite.DB()

    Sqlite3Stats.register_functions(db)

    SQLite.execute(db, "create table numbers (NUM1 float, NUM2 float)")
    for i = 1:10
        a = rand()
        b = rand() * 10
        SQLite.execute(db, "insert into numbers(num1, num2) values ($a, $b)")
    end

    @testset "QNORM" begin
        result =
            DBInterface.execute(
                db,
                "select QNORM(0.025, 0, 1) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"], [-1.95996], atol = tol)
    end

    @testset "PNORM" begin
        result =
            DBInterface.execute(
                db,
                "select PNORM(1.9599639845400576, 0, 1) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"], [0.0584451], atol = tol)
    end

    @testset "RNORM" begin
        result =
            DBInterface.execute(
                db,
                "select RNORM(0, 1) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test result[!, "MYRESULT"][1] < 10.0
        @test result[!, "MYRESULT"][1] > -10.0
    end

    @testset "RNORM" begin
        result =
            DBInterface.execute(
                db,
                "select RNORM(10, 1) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test result[!, "MYRESULT"][1] < 30.0
        @test result[!, "MYRESULT"][1] > -10.0
    end

    SQLite.close(db)

end

@testset "T Distribution" begin

    tol = 0.001

    db = SQLite.DB()

    Sqlite3Stats.register_functions(db)

    SQLite.execute(db, "create table numbers (NUM1 float, NUM2 float)")
    for i = 1:10
        a = rand()
        b = rand() * 10
        SQLite.execute(db, "insert into numbers(num1, num2) values ($a, $b)")
    end

    @testset "PT" begin
        result =
            DBInterface.execute(
                db,
                "select PT(1.9599639845400576,30) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], 0.970326, atol = tol)
    end

    @testset "QT" begin
        result =
            DBInterface.execute(
                db,
                "select QT(0.025, 30) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], -2.04227, atol = tol)
    end

    @testset "RT" begin
        result =
            DBInterface.execute(
                db,
                "select RT(30) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test result[!, "MYRESULT"][1] < 30.0
        @test result[!, "MYRESULT"][1] > -10.0
    end

    SQLite.close(db)
end

@testset "ChiSquare Distribution" begin

    tol = 0.001

    db = SQLite.DB()

    Sqlite3Stats.register_functions(db)

    SQLite.execute(db, "create table numbers (NUM1 float, NUM2 float)")
    for i = 1:10
        a = rand()
        b = rand() * 10
        SQLite.execute(db, "insert into numbers(num1, num2) values ($a, $b)")
    end

    @testset "PCHISQ" begin
        result =
            DBInterface.execute(
                db,
                "select PCHISQ(30, 30) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], 0.534346, atol = tol)
    end

    @testset "QCHISQ" begin
        result =
            DBInterface.execute(
                db,
                "select QCHISQ(0.05, 30) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], 18.4927, atol = tol)
    end

    @testset "RCHISQ" begin
        result =
            DBInterface.execute(
                db,
                "select RCHISQ(30) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test result[!, "MYRESULT"][1] < 100.0
        @test result[!, "MYRESULT"][1] >= 0.0
    end

    SQLite.close(db)
end

@testset "F Distribution" begin

    tol = 0.001

    db = SQLite.DB()

    Sqlite3Stats.register_functions(db)

    SQLite.execute(db, "create table numbers (NUM1 float, NUM2 float)")
    for i = 1:10
        a = rand()
        b = rand() * 10
        SQLite.execute(db, "insert into numbers(num1, num2) values ($a, $b)")
    end

    @testset "PF" begin
        result =
            DBInterface.execute(
                db,
                "select PF(0.1109452, 3, 5) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], 0.0499999, atol = tol)
    end

    @testset "QF" begin
        result =
            DBInterface.execute(
                db,
                "select QF(0.05, 3, 5) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], 0.1109452, atol = tol)
    end

    @testset "RF" begin
        result =
            DBInterface.execute(
                db,
                "select RF(3, 5) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test result[!, "MYRESULT"][1] < 100.0
        @test result[!, "MYRESULT"][1] >= 0.0
    end

    SQLite.close(db)
end

@testset "Poisson Distribution" begin

    tol = 0.001

    db = SQLite.DB()

    Sqlite3Stats.register_functions(db)

    SQLite.execute(db, "create table numbers (NUM1 float, NUM2 float)")
    for i = 1:10
        a = rand()
        b = rand() * 10
        SQLite.execute(db, "insert into numbers(num1, num2) values ($a, $b)")
    end

    @testset "PPOIS" begin
        result =
            DBInterface.execute(
                db,
                "select PPOIS(1000000, 1000000) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], 0.5, atol = tol)
    end

    @testset "QPOIS" begin
        result =
            DBInterface.execute(
                db,
                "select QPOIS(0.50, 5) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], 5.0, atol = tol)
    end

    @testset "RPOIS" begin
        result =
            DBInterface.execute(
                db,
                "select RPOIS(5) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test result[!, "MYRESULT"][1] < 100.0
        @test result[!, "MYRESULT"][1] >= 0.0
    end

    SQLite.close(db)
end

@testset "Binomial Distribution" begin

    tol = 0.001

    db = SQLite.DB()

    Sqlite3Stats.register_functions(db)

    SQLite.execute(db, "create table numbers (NUM1 float, NUM2 float)")
    for i = 1:10
        a = rand()
        b = rand() * 10
        SQLite.execute(db, "insert into numbers(num1, num2) values ($a, $b)")
    end

    @testset "PBINOM" begin
        result =
            DBInterface.execute(
                db,
                "select PBINOM(5, 10, 0.5) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], 0.62304687, atol = tol)
    end

    @testset "QBINOM" begin
        result =
            DBInterface.execute(
                db,
                "select QBINOM(0.50, 10, 0.5) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], 5.0, atol = tol)
    end

    @testset "RBINOM" begin
        result =
            DBInterface.execute(
                db,
                "select RBINOM(10, 0.5) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test result[!, "MYRESULT"][1] <= 10.0
        @test result[!, "MYRESULT"][1] >= 0.0
    end

    SQLite.close(db)
end

@testset "Uniform Distribution" begin

    tol = 0.001

    db = SQLite.DB()

    Sqlite3Stats.register_functions(db)

    SQLite.execute(db, "create table numbers (NUM1 float, NUM2 float)")
    for i = 1:10
        a = rand()
        b = rand() * 10
        SQLite.execute(db, "insert into numbers(num1, num2) values ($a, $b)")
    end

    @testset "PUNIF" begin
        result =
            DBInterface.execute(
                db,
                "select PUNIF(5, 0, 10) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], 0.5, atol = tol)
    end

    @testset "QUNIF" begin
        result =
            DBInterface.execute(
                db,
                "select QUNIF(0.5, 5, 10) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], 7.5, atol = tol)
    end

    @testset "RUNIF" begin
        result =
            DBInterface.execute(
                db,
                "select RUNIF(0, 10) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test result[!, "MYRESULT"][1] <= 10.0
        @test result[!, "MYRESULT"][1] >= 0.0
    end

    SQLite.close(db)
end

@testset "Exponential Distribution" begin

    tol = 0.001

    db = SQLite.DB()

    Sqlite3Stats.register_functions(db)

    SQLite.execute(db, "create table numbers (NUM1 float, NUM2 float)")
    for i = 1:10
        a = rand()
        b = rand() * 10
        SQLite.execute(db, "insert into numbers(num1, num2) values ($a, $b)")
    end

    @testset "PEXP" begin
        result =
            DBInterface.execute(
                db,
                "select PEXP(5, 10) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], 0.39346934, atol = tol)
    end

    @testset "QEXP" begin
        result =
            DBInterface.execute(
                db,
                "select QEXP(0.5, 10) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], 6.93147180, atol = tol)
    end

    @testset "REXP" begin
        result =
            DBInterface.execute(
                db,
                "select REXP(10.0) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test result[!, "MYRESULT"][1] <= 100.0
        @test result[!, "MYRESULT"][1] >= 0.0
    end

    SQLite.close(db)
end

@testset "Beta Distribution" begin

    tol = 0.001

    db = SQLite.DB()

    Sqlite3Stats.register_functions(db)

    SQLite.execute(db, "create table numbers (NUM1 float, NUM2 float)")
    for i = 1:10
        a = rand()
        b = rand() * 10
        SQLite.execute(db, "insert into numbers(num1, num2) values ($a, $b)")
    end

    @testset "PBETA" begin
        result =
            DBInterface.execute(
                db,
                "select PBETA(0.5, 0.12, 1) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], 0.92018, atol = tol)
    end

    @testset "QBETA" begin
        result =
            DBInterface.execute(
                db,
                "select QBETA(0.5, 0.12, 1) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], 0.00310, atol = tol)
    end

    @testset "RBETA" begin
        result =
            DBInterface.execute(
                db,
                "select RBETA(0.12, 1.0) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test result[!, "MYRESULT"][1] <= 1.0
        @test result[!, "MYRESULT"][1] >= 0.0
    end

    SQLite.close(db)
end

@testset "Cauchy Distribution" begin

    tol = 0.001

    db = SQLite.DB()

    Sqlite3Stats.register_functions(db)

    SQLite.execute(db, "create table numbers (NUM1 float, NUM2 float)")
    for i = 1:10
        a = rand()
        b = rand() * 10
        SQLite.execute(db, "insert into numbers(num1, num2) values ($a, $b)")
    end

    @testset "PCAUCHY" begin
        result =
            DBInterface.execute(
                db,
                "select PCAUCHY(0, 0, 1) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], 0.5, atol = tol)
    end

    @testset "QCAUCHY" begin
        result =
            DBInterface.execute(
                db,
                "select QCAUCHY(0.05, 0, 1) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test isapprox(result[!, "MYRESULT"][1], -6.313751, atol = tol)
    end

    @testset "RCAUCHY" begin
        result =
            DBInterface.execute(
                db,
                "select RCAUCHY(0, 1) as MYRESULT from numbers limit 1",
            ) |> DataFrame
        @test result[!, "MYRESULT"][1] isa Number
        @test result[!, "MYRESULT"][1] isa Number
    end

    SQLite.close(db)
end


text:
#=%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example : Find the distribution of a tracer given:              %
%       (a) the pathways described by A,                          %
%       (b) interior sources and sinks given by dC,               %
%           that best fits observations, Cobs,                    %
%   and (c) inequality constraints on the tracer concentration.   %
%                                                                  %
% Mathematically, minimize J = (C-Cobs)^T W (C-Cobs) subject to    %
%                         AC = d + Gamma u                         %
%  where u is the estimated change in surface concentration.    %
%
% See Supplementary Section 2, Gebbie & Huybers 2011.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% =#

using Revise, TMI, GoogleDrive
using PyPlot, PyCall, Test
#, Distributions, LinearAlgebra,  Zygote, ForwardDiff, Optim

TMIversion = "modern_90x45x33_GH10_GH12"
A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion)

# first guess of change to surface boundary conditions
# ocean values are 0
#u = zerosurfaceboundary(γ)
u = (;surface = zerosurfaceboundary(γ))
uvec = vec(u)

# take synthetic, noisy observations
y, W⁻, ctrue = synthetic_observations(TMIversion,"θ",γ)

# a first guess: observed surface boundary conditions are perfect.
# set surface boundary condition to the observations.
# below surface = 0 % no internal sinks or sources.
#b = getsurfaceboundary(y)
b = (;surface = getsurfaceboundary(y))

# get sample J value
F,G = costfunction_gridded_obs(uvec,Alu,b,u,y,W⁻,γ)
fg!(F,G,x) = costfunction_gridded_obs!(F,G,x,Alu,b,u,y,W⁻,γ)

fg(x) = costfunction_gridded_obs(x,Alu,b,u,y,W⁻,γ)
f(x) = fg(x)[1]
J₀,gJ₀ = fg(uvec)

#### gradient check ###################
# check with forward differences
ϵ = 1e-3
ii = rand(1:sum(γ.wet[:,:,1]))
println("Location for test =",ii)
δu = copy(uvec); δu[ii] += ϵ
∇f_finite = (f(δu) - f(uvec))/ϵ
println(∇f_finite)

fg!(J₀,gJ₀,(uvec+δu)./2) # J̃₀ is not overwritten
∇f = gJ₀[ii]
println(∇f)

# error less than 10 percent?
println("Percent error ",100*abs(∇f - ∇f_finite)/abs(∇f + ∇f_finite))
#### end gradient check #################

# filter the data with an Optim.jl method
iterations = 15
out = steadyclimatology(uvec,fg!,iterations)

# reconstruct by hand to double-check.
ũ = unvec(u,out.minimizer)

# apply to the boundary conditions
b̃ = adjustboundarycondition(b,ũ)

# reconstruct tracer map
c₀ = steadyinversion(Alu,b,γ)
c̃  = steadyinversion(Alu,b̃,γ)

Δc̃ = c̃ - ctrue
Δc₀ = c₀ - ctrue

# plot the difference
level = 15 # your choice 1-33
depth = γ.depth[level]
cntrs = -1:0.05:1
label = "Optimized misfit: Δc̃"
# Help: needs work with continents and labels
planviewplot(Δc̃, depth, cntrs, titlelabel=label)

label = "First guess misfit: Δc₀"
planviewplot(Δc₀, depth, cntrs, titlelabel=label)

# future step: box minimization to eliminate unreasonable temperatures.
# should be made more generic for other tracers
# From MATLAB:
#lbT = -2.*ones(Nsfc,1); % temperature lower bound: can not freeze.
#ubT = 40.*ones(Nsfc,1); % ad-hoc temperature upper bound: 40 C.


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

# Check to make sure types are imported properly
@test RadioMenu <: TerminalMenus.AbstractMenu

# Invalid Menu Params
@test_throws ErrorException RadioMenu(["one"])
@test_throws ErrorException RadioMenu(["one", "two", "three"], pagesize=1)

# Constructor
@test RadioMenu(["one", "two", "three"]).pagesize == 3
@test RadioMenu(string.(1:30), pagesize=-1).pagesize == 30
@test RadioMenu(string.(1:4), pagesize=10).pagesize == 4
@test RadioMenu(string.(1:100)).pagesize == 10

radio_menu = RadioMenu(string.(1:20))
@test TerminalMenus.options(radio_menu) == string.(1:20)
radio_menu.selected = 2
TerminalMenus.cancel(radio_menu)
@test radio_menu.selected == -1
@test TerminalMenus.header(radio_menu) == ""

# Output
TerminalMenus.config() # Use default chars
CONFIG = TerminalMenus.CONFIG

radio_menu = RadioMenu(string.(1:10))
buf = IOBuffer()
TerminalMenus.writeLine(buf, radio_menu, 1, true)
@test String(take!(buf)) == string(CONFIG[:cursor], " 1")
TerminalMenus.config(cursor='+')
TerminalMenus.writeLine(buf, radio_menu, 1, true)
@test String(take!(buf)) == "+ 1"
TerminalMenus.config(charset=:unicode)
TerminalMenus.writeLine(buf, radio_menu, 1, true)
@test String(take!(buf)) == string(CONFIG[:cursor], " 1")

# Test using stdin
radio_menu = RadioMenu(string.(1:10))
@test simulateInput(3, radio_menu, :down, :down, :enter)

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
let
    mdp = SimpleGridWorld()

    pomdp = FullyObservablePOMDP(mdp)

    @test observations(pomdp) == states(pomdp)
    @test statetype(pomdp) == obstype(pomdp)

    @test observations(pomdp) == states(pomdp)
    @test statetype(pomdp) == obstype(pomdp)

    s_po = initialstate(pomdp, MersenneTwister(1))
    s_mdp = initialstate(mdp, MersenneTwister(1))
    @test s_po == s_mdp

    solver = ValueIterationSolver(max_iterations = 100)
    mdp_policy = solve(solver, mdp)
    pomdp_policy = solve(solver, UnderlyingMDP(pomdp))
    @test mdp_policy.util == pomdp_policy.util

    is = initialstate(mdp, MersenneTwister(3))
    for (sp, o, r) in stepthrough(pomdp,
                               FunctionPolicy(o->:left),
                               PreviousObservationUpdater(),
                               is, is, "sp,o,r",
                               rng=MersenneTwister(2),
                               max_steps=10)
        @test sp == o
    end
end


text:
export num_rho_adaptions
# set initial values of rhoVec
function set_rho_vec!(ws::COSMO.Workspace{T}) where {T <: AbstractFloat}
	m = ws.p.model_size[1]

	ws.ρ = ws.settings.rho
	ws.ρvec = ws.ρ * ones(T, m)

	apply_constraint_rho_scaling!(ws.ρvec, ws.row_ranges, ws.p.C.sets, ws.settings.RHO_MIN, ws.settings.RHO_EQ_OVER_RHO_INEQ)

	push!(ws.rho_updates, ws.ρ)
	return nothing
end

"Scale active constraints by `RHO_EQ_OVER_RHO_INEQ` and set loose constraints to `RHO_MIN`."
function apply_constraint_rho_scaling!(ρvec::AbstractVector{T}, row_ranges::Array{UnitRange{Int}, 1}, sets::Vector{AbstractConvexSet}, RHO_MIN::T, RHO_EQ_OVER_RHO_INEQ::T) where {T <: AbstractFloat}
	for (k, C) in enumerate(sets)
		if C isa ZeroSet
			for row in row_ranges[k]
				ρvec[row] *= RHO_EQ_OVER_RHO_INEQ
			end
		elseif C isa Nonnegatives
			row_offset = row_ranges[k].start - 1
			for (j, val) in enumerate(C.constr_type)
				if val == true
					ρvec[row_offset + j] = RHO_MIN
				end
			end

		elseif C isa Box
			row_offset = row_ranges[k].start - 1
			_rho_scale_box_bounds!(ρvec, C.constr_type, row_offset, RHO_MIN::T, RHO_EQ_OVER_RHO_INEQ::T)
		end
	end
	return nothing
end

function _rho_scale_box_bounds!(ρvec::AbstractVector{T}, constr_type::Vector{Int}, row_offset::Int, RHO_MIN::T, RHO_EQ_OVER_RHO_INEQ::T) where {T <: AbstractFloat}
	@inbounds for j = eachindex(constr_type)
		# loose inequality constraint
		if constr_type[j] == -1
			ρvec[j + row_offset] = RHO_MIN
		# equality constraint
		elseif constr_type[j] == 1
			ρvec[j + row_offset] *= RHO_EQ_OVER_RHO_INEQ
		end
	end
end

# adapt rhoVec based on residual ratio
function adapt_rho_vec!(ws::COSMO.Workspace{T}) where {T <: AbstractFloat}
	settings = ws.settings
	# compute normalized residuals based on the working variables (dont unscale)
	ignore_scaling = true
	r_prim::T, r_dual::T = calculate_residuals!(ws, ignore_scaling)
	max_norm_prim::T, max_norm_dual::T  = max_res_component_norm(ws, ignore_scaling)

	r_prim = r_prim / (max_norm_prim + T(1e-10))
	r_dual = r_dual / (max_norm_dual + T(1e-10))

	new_rho = ws.ρ * sqrt(r_prim / (r_dual + T(1e-10)))
	new_rho = min(max(new_rho, settings.RHO_MIN), settings.RHO_MAX)
	# only update rho if significantly different than current rho
	if (new_rho > settings.adaptive_rho_tolerance * ws.ρ) || (new_rho < (one(T) ./ settings.adaptive_rho_tolerance) * ws.ρ)
		update_rho_vec!(new_rho, ws)
		return true
	else
		return false
	end
end

function update_rho_vec!(new_rho::T, ws::COSMO.Workspace{T}) where {T <: AbstractFloat}

	ws.ρ     = new_rho
	ws.ρvec .= new_rho

	apply_constraint_rho_scaling!(ws.ρvec, ws.row_ranges, ws.p.C.sets, ws.settings.RHO_MIN, ws.settings.RHO_EQ_OVER_RHO_INEQ)

	# log rho updates
	push!(ws.rho_updates, new_rho)

	if ws.settings.verbose_timing
		ws.times.factor_update_time += @elapsed update_rho!(ws.kkt_solver, ws.ρvec)
	else
		update_rho!(ws.kkt_solver, ws.ρvec)
	end

	return nothing
end

"Return the number of rho adaptions so far."
num_rho_adaptions(rho_updates::Vector{T}) where {T <: AbstractFloat} = length(rho_updates) - 1
num_rho_adaptions(model::COSMO.Model{T}) where {T <: AbstractFloat} = num_rho_adaptions(model.rho_updates)


text:
#!/usr/bin/env julia

include("src/MLToolkit.jl")

# only run anything if it's being run from the command line
if !isinteractive()
	MLToolkit.run()
end

max_stars_repo_name: coreywoodfield/mltoolkit-julia

max_issues_repo_name: coreywoodfield/mltoolkit-julia

max_forks_repo_name: coreywoodfield/mltoolkit-julia


text:
module AutoViz

using Reexport
using Parameters
using AutomotiveDrivingModels

@reexport using Colors
@reexport using Cairo

import Reel
Reel.set_output_type("gif")

# using Reactive
# import Gtk
# using PGFPlots

export
        DEFAULT_CANVAS_WIDTH,
        DEFAULT_CANVAS_HEIGHT,
        render!

const DEFAULT_CANVAS_WIDTH = 1000
const DEFAULT_CANVAS_HEIGHT = 600

include(Pkg.dir("AutoViz", "src", "colorscheme.jl"))
include(Pkg.dir("AutoViz", "src", "rendermodels.jl"))
include(Pkg.dir("AutoViz", "src", "camera.jl"))

include(Pkg.dir("AutoViz", "src", "render_roadways.jl"))
include(Pkg.dir("AutoViz", "src", "render_vehicles.jl"))
include(Pkg.dir("AutoViz", "src", "render_scenes.jl"))

include(Pkg.dir("AutoViz", "src", "overlays.jl"))

include(Pkg.dir("AutoViz", "src", "reel_drive.jl"))
# include(Pkg.dir("AutoViz", "src", "gtk.jl"))

end # module


text:
# automatically generated by the FlatBuffers compiler, do not modify

# module: Example

import FlatBuffers

FlatBuffers.@UNION(AnyAmbiguousAliases, (
    Nothing,
    Monster,
    Monster,
    Monster,
))


text:
module Floodfill
include("./environments.jl")
#===
Flood-fill (node):
    1. Set Q to the empty queue or stack.
    2. Add node to the end of Q.
    3. While Q is not empty:
    4.   Set n equal to the first element of Q.
    5.   Remove first element from Q.
    6.   If n is Inside:
            Set the n
            Add the node to the west of n to the end of Q.
            Add the node to the east of n to the end of Q.
            Add the node to the north of n to the end of Q.
            Add the node to the south of n to the end of Q.
    7. Continue looping until Q is exhausted.
    8. Return.
===#
function floodfill(_phi::Array, N, L, filled, beginx, beginy, filled_index, indexI = nothing, multiprocess = true)
    # 始点は平面全体の縁を一周すべき
    # -> 閉曲線が2箇所で境界に接していたりすると、その領域のみで色塗り(符合つけ)が終わってしまうから。
    point_que = [(beginx, beginy)]
    closed_zero = L * 2 * 1.42 / (N + 1)
    # println("The lattice size: ",size(_phi), " the beginning point: ", point_que)
    if indexI !== nothing
        bounse_x = size(_phi[:, 1])[1] + 1
        bounse_min_x = 0
    else
        bounse_x = N + 1
        bounse_min_x = 0
    end
    bounse_y = size(_phi)[2] + 1
    bounse_min_y = 0

    STEP = 2^1
    # STEP=1
    while length(point_que) > 0
        _x = copy(point_que[1][1])
        _y = copy(point_que[1][2])
        # 現在の格子点の符合を反転させる
        if _phi[_x, _y] <= 0
            _phi[_x, _y] *= (-1)
        end
        filled[filled_index] = point_que[1]
        popfirst!(point_que)
        filled_index += 1
        # 下
        if bounse_min_x < _x - STEP < bounse_x && !((_x - STEP, _y) in filled) && abs(_phi[_x-STEP, _y]) > closed_zero
            if !((_x - STEP, _y) in point_que) && _phi[_x-STEP, _y] <= 0
                append!(point_que, [(_x - STEP, _y)]) # 下側の格子点をqueueに積む
            end
        end
        # 左
        if bounse_min_y < _y - STEP < bounse_y && !((_x, _y - STEP) in filled) && abs(_phi[_x, _y-STEP]) > closed_zero
            if !((_x, _y - STEP) in point_que) && _phi[_x, _y-STEP] <= 0
                append!(point_que, [(_x, _y - STEP)]) # 左側の格子点をqueueに積む
            end
        end
        # 上
        if bounse_min_x < _x + STEP < bounse_x && !((_x + STEP, _y) in filled) && abs(_phi[_x+STEP, _y]) > closed_zero
            if !((_x + STEP, _y) in point_que) && _phi[_x+STEP, _y] <= 0
                append!(point_que, [(_x + STEP, _y)]) # 上側の格子点をqueueに積む
            end
        end
        # 右
        if _y + STEP < bounse_y && !((_x, _y + STEP) in filled) && abs(_phi[_x, _y+STEP]) > closed_zero
            if !((_x, _y + STEP) in point_que) && _phi[_x, _y+STEP] <= 0
                append!(point_que, [(_x, _y + STEP)]) # 右側の格子点をqueueに積む
            end
        end
    end
    # end
    if STEP == 2
        _phi = assign_signs(_phi, STEP, N, L, multiprocess)
    end
    return _phi, filled_index#, filled
end
precompile(floodfill, (Array, Int, Float64, Array, Int, Int, Int, Int, Int))

## 今のところSTEP=2の場合のみ対応
function assign_signs(_phi, STEP, N, L, multiprocess = true)
    # 閉曲線内部が「-」である。デフォが
    # Int((N-1)/100) # 再帰回数
    loops = Int(log2(STEP))
    steps_signed_grid = STEP
    steps_unsigned_grid = Int(STEP / 2)
    # if JULIA_MULTI_PROCESS
    if multiprocess
        while loops > 0
            if STEP != 1
                Threads.@threads for i = 1:steps_signed_grid:length(_phi[1, :]) # 各行を一つ飛ばしで。
                    for j = 1:steps_signed_grid:length(_phi[:, 1])
                        # if j != length(_phi[:, 1]) # ケツでない
                        if j + steps_signed_grid <= length(_phi[:, 1]) && j + steps_unsigned_grid <= length(_phi[:, 1])
                            # 掛けたらマイナス->境界を示す. かつ 両方の距離を足したらgrid点の間の2倍の距離になる(数値誤差を考慮?できてる?)
                            if _phi[i, j] * _phi[i, j+steps_signed_grid] <= 0
                                if abs(_phi[i, j]) <= abs(_phi[i, j+steps_signed_grid]) # jの方が近い
                                    if _phi[i, j] < 0 # jが内
                                        _phi[i, j+steps_unsigned_grid] *= (-1) # 外側にある
                                    end
                                else # j+2の方が近い
                                    if _phi[i, j] >= 0 # jが外
                                        _phi[i, j+steps_unsigned_grid] *= (-1)
                                    end
                                end
                            elseif _phi[i, j] > 0 && _phi[i, j+steps_signed_grid] > 0# 外側
                                _phi[i, j+steps_unsigned_grid] *= (-1)
                            end
                        end
                    end
                end
                Threads.@threads for j = 1:Int(steps_signed_grid / 2):length(_phi[:, 1])
                    for i = 1:steps_signed_grid:length(_phi[:, 1])
                        # if i != length(_phi[:, 1]) # ケツでない
                        if i + steps_signed_grid <= length(_phi[:, 1]) && i + steps_unsigned_grid <= length(_phi[:, 1])
                            if _phi[i, j] * _phi[i+steps_signed_grid, j] <= 0
                                if abs(_phi[i, j]) <= abs(_phi[i+steps_signed_grid, j]) # iの方が近い
                                    if _phi[i, j] < 0 # iが内
                                        _phi[i+steps_unsigned_grid, j] *= (-1) # 外
                                    end
                                else # i+2の方が近い
                                    # if _phi[i+steps_signed_grid, j] < 0 # i+steps_signed_gridが内
                                    if _phi[i, j] >= 0 # iが外
                                        _phi[i+steps_unsigned_grid, j] *= (-1) # 外
                                    end
                                end
                            elseif _phi[i, j] > 0 && _phi[i+steps_signed_grid, j] > 0# 外側
                                _phi[i+steps_unsigned_grid, j] *= (-1)
                            end
                        end
                    end
                end
            end
            steps_signed_grid = steps_signed_grid >= 4 ? deepcopy(Int(steps_signed_grid / 2)) : undef
            steps_unsigned_grid = steps_unsigned_grid >= 2 ? deepcopy(Int(steps_unsigned_grid / 2)) : undef
            loops -= 1
        end
    else
        while loops > 0
            if STEP != 1
                for i = 1:steps_signed_grid:length(_phi[1, :]) # 各行を一つ飛ばしで。
                    for j = 1:steps_signed_grid:length(_phi[:, 1])
                        # if j != length(_phi[:, 1]) # ケツでない
                        if j + steps_signed_grid <= length(_phi[:, 1]) && j + steps_unsigned_grid <= length(_phi[:, 1])
                            # 掛けたらマイナス->境界を示す. かつ 両方の距離を足したらgrid点の間の2倍の距離になる(数値誤差を考慮?できてる?)
                            if _phi[i, j] * _phi[i, j+steps_signed_grid] <= 0
                                if abs(_phi[i, j]) < abs(_phi[i, j+steps_signed_grid]) # jの方が近い
                                    if _phi[i, j] < 0 # jが内
                                        _phi[i, j+steps_unsigned_grid] *= (-1) # 外側にある
                                    end
                                else # j+2の方が近い
                                    if _phi[i, j] >= 0 # jが外
                                        _phi[i, j+steps_unsigned_grid] *= (-1)
                                    end
                                end
                            elseif _phi[i, j] > 0 && _phi[i, j+steps_signed_grid] > 0# 外側
                                _phi[i, j+steps_unsigned_grid] *= (-1)
                            end
                        end
                    end
                end
                for j = 1:Int(steps_signed_grid / 2):length(_phi[:, 1])
                    for i = 1:steps_signed_grid:length(_phi[:, 1])
                        # if i != length(_phi[:, 1]) # ケツでない
                        if i + steps_signed_grid <= length(_phi[:, 1]) && i + steps_unsigned_grid <= length(_phi[:, 1])
                            if _phi[i, j] * _phi[i+steps_signed_grid, j] <= 0
                                if abs(_phi[i, j]) < abs(_phi[i+steps_signed_grid, j]) # jの方が近い
                                    if _phi[i, j] < 0 # jが内
                                        _phi[i+steps_unsigned_grid, j] *= (-1)
                                    end
                                else # j+2の方が近い
                                    if _phi[i, j] >= 0 # jが外
                                        _phi[i+steps_unsigned_grid, j] *= (-1)
                                    end
                                end
                            elseif _phi[i, j] > 0 && _phi[i+steps_signed_grid, j] > 0# 外側
                                _phi[i+steps_unsigned_grid, j] *= (-1)
                            end
                        end
                    end
                end
            end
            steps_signed_grid = steps_signed_grid >= 4 ? deepcopy(Int(steps_signed_grid / 2)) : undef
            steps_unsigned_grid = steps_unsigned_grid >= 2 ? deepcopy(Int(steps_unsigned_grid / 2)) : undef
            loops -= 1
        end
    end

    return _phi
end

function signining_field(_phi::Array, N, L, multiprocess = true)
    _phi .*= (-1)
    filled = Array{Tuple{Int64,Int64}}(undef, N * N) # N=100だと12倍速! N=200だと60倍速!
    filled_index = 1
    beginx = 1
    beginy = N
    indexI = 1
    _phi = floodfill(_phi, N, L, filled, beginx, beginy, filled_index, indexI, multiprocess)
    return _phi
end
precompile(signining_field, (Array, Int, Float64))
export signining_field
end

max_stars_repo_name: jabelic/Signed-distance-function

max_issues_repo_name: jabelic/Signed-distance-function

max_forks_repo_name: jabelic/Signed-distance-function


text:
module Losses

#==============================================================================#
#                          MODULE OPTIONS & FLAGS                              #
const DEBUG_MODE   = false
#==============================================================================#
const VERBOSE_MODE = true
#------------------------------------------------------------------------------#

using Printf
import Random, AutoGrad, Distributions
using ..Concepts
using ..Utilities.BatchUtils
# export BatchFactory

# struct SGD end

####################### HELPER METHODS ###################
# vectorized sigmoid function
σ(z) = 1.0 ./ (1.0 .+ exp.(-z))
##########################################################

# function Poisson()
#     L(x,y,c,ρ) = sum(exp.(x) .- y .* x) + sum(ρ .* (x .- c).^2)
#     return L;
# end

# function Logistic()
#     L(x,y,c,ρ) = -sum(y .* log.(σ.(x)) .+ (1 .- y) .* log.(1 .- σ.(x))) .+  sum(ρ .* (x .- c).^2);

#     return L;
# end

export Loss,
  train,
  sgd_train,
  subgrad_train,
  negative_binomial_train

struct Loss{T} <: AbstractLoss where T<:Any
  function Loss{T}() where T<:Any
    @abstract_instance
    return new{T}()
  end

  #    function Loss{AbstractPoisson}(of::Union{AbstractPoisson,Type{Val{:Poisson}}})
  function Loss{T}(of::T) where T<:ExponentialFamily
    return new{T}()
  end
end

# kind of hackish.. the second one just defined in the first.
const Loss(of::Union{T,Symbol}) where T<:ExponentialFamily =
  typeof(of) <: ExponentialFamily ? Loss{T}(of) : Loss(type_conversion(ExponentialFamily, of))

loss_logistic(x,y,c,ρ) = -sum(y .* log.(σ.(x)) .+ (1 .- y) .* log.(1 .- σ.(x))) .+  sum(ρ .* (x .- c).^2);

 #==============================================================================#
#                             Gaussian Loss                                    #
#==============================================================================#
function Concepts.provide(loss::Loss{AbstractGaussian})
  #TODO
end

function Concepts.evaluate(loss::Loss{AbstractGaussian},
                  x,y,c,ρ)
  #TODO
  return sum(0.5 .* (y .- x).^2) + sum(ρ .* (x .- c).^2);
end

function grad(loss::Loss{AbstractGaussian} ,
              x,y,c,ρ)
  #TODO
  return -(y .- x) .+ (2*ρ) .* (x .- c)
end
#==============================================================================#
#                         Bernoulli (Logistic) Loss                            #
#==============================================================================#
function Concepts.provide(loss::Loss{AbstractBernoulli})
  L(x,y,c,ρ) = -sum(y .* log.(σ.(x)) .+ (1 .- y) .* log.(1 .- σ.(x))) .+  sum(ρ .* (x .- c).^2);
  return L
end

function Concepts.evaluate(loss::Loss{AbstractBernoulli},
                           x,y,c,ρ)
  return  -sum(y .* log.(σ.(x)) .+ (1 .- y) .* log.(1 .- σ.(x)))
  .+ sum(ρ .* (x .- c).^2);
end

function grad(loss::Loss{AbstractBernoulli},
              x,y,c,ρ)
  ex = exp.(x)
  inv_ex1 = 1 ./(ex .+ 1)
  return inv_ex1 .* (-y + (1 .-y) .* ex) .+ (2*ρ) .* (x .- c)
end

#==============================================================================#
#                             Poisson Loss                                     #
#==============================================================================#
@overload
function Concepts.provide(loss::Loss{AbstractPoisson})
  L(x,y,c,ρ) = sum(exp.(x) .- y .* x) + sum(ρ .* (x .- c).^2)
  return L;
end

function Concepts.evaluate(loss::Loss{AbstractPoisson},
                           x,y,c,ρ)
  return sum(exp.(x) .- y .* x) + sum(ρ .* (x .- c).^2)
end

function grad(loss::Loss{AbstractPoisson},
              x,y,c,ρ)
  return exp.(x) .- y .+ (2*ρ) .* (x .- c)
end

#==============================================================================#
#                              Gamma Loss                                      #
#==============================================================================#
function Concepts.provide(loss::Loss{AbstractGamma})
  L(x,y,c,ρ) = sum(-x .* y .- log.(-x))+ sum(ρ .* (x .- c).^2)
  return L
end

# function Concepts.evaluate(loss::Loss{AbstractGamma},x,y,c,ρ)
#     return sum(x .* y .- log.(x))+ sum(ρ .* (x .- c).^2)
# end

function Concepts.evaluate(loss::Loss{AbstractGamma},x,y,c,ρ)
  return sum(y .* exp.(x) .- x) + sum(ρ .* (x .- c).^2)
end

# function _evaluate(loss::Loss{AbstractGamma},x,y,c,ρ)
#   x₋ = x[findall(a -> a < 0, x)]
#   y₋ = y[findall(a -> a < 0, x)]
#   x₊ = x[findall(a -> a > 0, x)]
#   y₊ = x[findall(a -> a > 0, x)]
#   return (-1) * sum(x₋ .* y₋ .+ log.(-x₋)) + (-1) * sum(-x₊ .* y₊ .+ log.(x₊)) + sum(ρ .* (x .- c).^2)
# end

## Use the reciprocal link instead of the negative reciprocal link
# function grad(loss::Loss{AbstractGamma}, x, y, c, ρ)
#     return y .- (1 ./ x) .+ (2*ρ) .* (x .- c)
# end

function grad(loss::Loss{AbstractGamma}, x, y, c, ρ)
  return y .* exp.(x) .- 1 + (2*ρ) .* (x .- c)
end

function subgrad(loss::Loss{AbstractGamma}, x, y, c, ρ)
  ∇ = zeros(length(x))
  pos_id = findall(a -> a > 0, x)
  neg_id = findall(a -> a < 0, x)
  ∇[neg_id] = (-1) .* (y[neg_id] - (1 ./ x[neg_id])) .+ (2*ρ) .* (x[neg_id] .- c[neg_id])
  ∇[pos_id] = (-1) .* (-y[pos_id] + (1 ./ x[pos_id])) .+ (2*ρ) .* (x[pos_id] .- c[pos_id])
  return ∇
end

function grad_logistic(x,y,c,ρ)
  ex = exp.(x)
  inv_ex1 = 1 ./(ex .+ 1);
  return inv_ex1 .* (-y + (1 .-y) .* ex) .+ (2*ρ) .* (x .- c);
  # return (-y .* inv_ex1 + (1 .- y) .* (ex .* inv_ex1)) .+ (2*ρ) .* (x.-c);
end

#==============================================================================#
#                         Negative Binomial Loss                               #
#==============================================================================#
function Concepts.provide(loss::Loss{AbstractNegativeBinomial})
  ## to implement
  return nothing
end

function Concepts.evaluate(loss::Loss{AbstractNegativeBinomial}, x, y, c, ρ; r_estimate)
  return sum(y .* exp.(x) - r_estimate .* log.(1 .- exp.(-exp.(x)))) .+ sum(ρ .* (x .- c).^2)
  # originally it is. However, this is a constrained optimization problem.
  # return sum(-y .* x .- r_estimate .* log.(1 .- exp.(x))) + sum(ρ .* (x .- c).^2)
end

function grad(loss::Loss{AbstractNegativeBinomial}, x, y, c, ρ; r_estimate = nothing)
  return y .* exp.(x) .- r_estimate .* exp.(x) ./ (exp.(exp.(x)) .- 1)  + (2*ρ) .* (x .- c)
end

function train(loss;fx,y,c,ρ,γ=0.02,iter=20,verbose=false)
  DEBUG_MODE && @info "Gradient Descent with Autograd"
  ∇ = AutoGrad.grad(loss);
  curFx = fx;
  for i = 1:iter
    curFx = curFx .- γ .* ∇(curFx,y,c,ρ);
  end
  return curFx;
end

function train(native_loss::Loss{T};
               fx, y, c, ρ, γ=0.02, iter=20, verbose=false, subgrad = false) where T<:ExponentialFamily
  DEBUG_MODE && @info "Gradient Descent with native differentitaion"
  curFx = fx;
  for i = 1:iter
    curFx .-= γ * grad(native_loss, curFx, y, c, ρ);
    # curFx .-= γ * grad(native_loss, curFx, y, c, ρ);
    # if project == true
    #   curFx = abs.(curFx)
    # end
    if verbose == true
      @printf("loss:%f\n", Concepts.evaluate(native_loss,curFx,y,c,ρ ))
    end
  end
  return curFx;
end

# specialization for negative binomial loss
function negative_binomial_train(;fx, y, c, ρ, γ=0.02, iter=20, verbose=false, r_estimate = nothing)
  DEBUG_MODE && @info "Gradient Descent with native differentitaion"
  curFx = fx;
  for i = 1:iter
    # @show("here")
    curFx .-= γ * grad(Loss{AbstractNegativeBinomial}(), curFx, y, c, ρ;r_estimate = r_estimate);
    if verbose == true
      @printf("loss:%f\n", Concepts.evaluate(Loss{AbstractNegativeBinomial}(),
                                             curFx,y,c,ρ; r_estimate = r_estimate))
    end
  end
  return curFx;
end

function subgrad_train(native_loss::Loss{T};
                       fx, y, c, ρ, γ=0.02, iter=20, verbose=false) where T<:ExponentialFamily
  DEBUG_MODE && @info "Gradient Descent with native differentitaion"
  curFx = fx;
  for i = 1:iter
    curFx .-= γ * subgrad(native_loss, curFx, y, c, ρ);
    if verbose == true
      @printf("loss:%f\n",_evaluate(native_loss,curFx,y,c,ρ ))
    end
  end
  return curFx;
end

function sgd_train(native_loss::Loss{T};
                   fx, y, c, ρ, α, ρ₁, ρ₂, batch_size, epoch) where T<:ExponentialFamily
  n = length(fx)
  curFx = fx
  batch = BatchFactory{SequentialScan}(size = batch_size)
  initialize(batch, fx)
  s = zeros(batch_size)
  r = zeros(batch_size)
  ŝ = zeros(batch_size)
  r̂ = zeros(batch_size)
  for i in 1:epoch
    while has_next(batch)
      cur_batch = next(batch)
      ∇ₛₐₘₚₗₑ = grad(native_loss, curFx[cur_batch], y[cur_batch], c[cur_batch], ρ)
      s .= (ρ₁ .* s) .+ (1 .- ρ₁) .* ∇ₛₐₘₚₗₑ
      r .= (ρ₂ .* r) .+ (1 .- ρ₂) .* (∇ₛₐₘₚₗₑ.^2)
      ŝ .= s ./ (1 - ρ₁^i)
      r̂ .= r ./ (1 - ρ₂^i)
      # @show(r̂)
      curFx[cur_batch] = curFx[cur_batch] - α ./ sqrt.(r̂) .* ŝ
    end
    # @show(Concepts.evaluate(native_loss,curFx,y,c,ρ))
    reset(batch)
  end
  return curFx
end

function sgd_subgrad_train(native_loss::Loss{T};
                           fx, y, c, ρ, α, ρ₁, ρ₂, batch_size, epoch) where T<:ExponentialFamily
  n = length(fx)
  curFx = fx
  batch = BatchFactory{SequentialScan}(size = batch_size)
  initialize(batch, fx)
  s = zeros(batch_size)
  r = zeros(batch_size)
  ŝ = zeros(batch_size)
  r̂ = zeros(batch_size)
  for i in 1:epoch
    while has_next(batch)
      cur_batch = next(batch)
      ∇ₛₐₘₚₗₑ = subgrad(native_loss, curFx[cur_batch], y[cur_batch], c[cur_batch], ρ)
      s .= (ρ₁ .* s) .+ (1 .- ρ₁) .* ∇ₛₐₘₚₗₑ
      r .= (ρ₂ .* r) .+ (1 .- ρ₂) .* (∇ₛₐₘₚₗₑ.^2)
      ŝ .= s ./ (1 - ρ₁^i)
      r̂ .= r ./ (1 - ρ₂^i)
      # @show(r̂)
      curFx[cur_batch] = curFx[cur_batch] - α ./ sqrt.(r̂) .* ŝ
    end
    # @show(Concepts.evaluate(native_loss,curFx,y,c,ρ))
    reset(batch)
  end
  return curFx
end

end


text:
function InnerProduct(in1,in2)

	return convert(Float32,real(sum(conj(in1[:]).*in2[:])))

end


text:
@testset "NVTX" begin

# markers

NVTX.mark("test")

# ranges

NVTX.@range "test" begin
end

NVTX.@range function test()
end
test()

@eval test2() = nothing

NVTX.@range function Main.test2(::Int)
end

NVTX.@range function Main.test2(::T) where T
end

NVTX.@range test3() = nothing

end


text:
export IrisEnv, fitness, new_gen!

import RDatasets

"Iris BERL environment"
mutable struct IrisEnv <: BERLenv
    name::String
    memory::Array # Array of storable data
    X::Array{Float64, 2}
    y::Array{Float64, 2}
end

environments["iris"]=IrisEnv

function IrisEnv(cfg::Dict)
    cfg["env"]="Iris"
    cfg["n_in"]=4
    cfg["n_out"]=3
    iris = RDatasets.dataset("datasets", "iris")
    X = convert(Matrix, iris[:, 1:4])'
    X = X ./ maximum(X; dims=2) # normalize
    r = iris[:, 5].refs # labels
    # 1-hot encoding:
    y = zeros(maximum(r), size(X, 2))
    for i in 1:length(r)
        y[r[i], i] = 1.0
    end
    IrisEnv("iris", [], X, y)
end

"Returns the fitness of one Cambrian Individual in this environment"
function fitness(indiv::Cambrian.Individual, env::IrisEnv)
    accuracy = 0.0
    for i in 1:size(env.X, 2)
        out = process(indiv, env.X[:, i])
        if argmax(out) == argmax(env.y[:, i])
            accuracy += 1
        end
    end
    f = [accuracy / size(env.X, 2)]
    try
        Integer.(round.(f))
    catch
        println(f, accuracy, size(env.X, 2))
    end
    f
end

"Sets up the env for the next generation"
function new_gen!(env::IrisEnv, cfg::Dict=nothing)

end


text:
module DCMQ
using AMQPClient
using DICOM

export consumer_loop, publish_nifti

function publish(channel, routing_key, dcmout; uri="")
    if tag"Pixel Data" in keys(dcmout)
        delete!(dcmout, tag"Pixel Data")
    end
    io = IOBuffer()
    dcm_write(io, dcmout)
    headers = Dict{String,Any}()
    if uri != ""
        headers["uri"] = uri
    end
    msg_out = Message(io.data, headers=headers)
    basic_publish(channel, msg_out; exchange="dicom", routing_key=routing_key)
end

function consumer_loop(server, queue, methods, dcmhandler; port = AMQPClient.AMQP_DEFAULT_PORT)
    login = "guest"
    password = "guest"
    auth_params = Dict{String,Any}("MECHANISM"=>"AMQPLAIN", "LOGIN"=>login, "PASSWORD"=>password)

    conn = connection(;virtualhost="/", host=server, port=port, auth_params=auth_params)
    chan = channel(conn, AMQPClient.UNUSED_CHANNEL, true)
    queue_declare(chan, queue)
    for method = methods
        queue_bind(chan, queue, "amq.topic", method)
    end

    println("entering polling loop")
    while true
        maybe_msg = basic_get(chan, queue, false)
        # check if we got a message
        if !isnothing(maybe_msg)
            println("message received")
            msg = maybe_msg
            try
                io = IOBuffer(msg.data)
                d = dcm_parse(io)
                uri = extractHeaders(msg)["uri"]
                dcmhandler(chan, d, uri)
                basic_ack(chan, msg.delivery_tag)
            catch e
                println(e)
                basic_reject(chan, msg.delivery_tag; requeue=true)
                throw(e)
            end
        else
            sleep(1)
        end
    end

    if isopen(conn)
        close(conn)
        # close is an asynchronous operation. To wait for the negotiation to complete:
        AMQPClient.wait_for_state(conn, AMQPClient.CONN_STATE_CLOSED)
    end
end

"Return AMQP messsage headers as dictionary."
function extractHeaders(msg::AMQPClient.Message)
    d = Dict()
    for x ∈ msg.properties[:headers].data
        key = String(copy(x.name.data))
        value = String(copy(x.val.fld.data))
        d[key] = value
    end
    return d
end

end # module


text:
# using Revise

using Test

using IncrementalInference, Distributions
using DistributedFactorGraphs
using Statistics
# going to introduce two new constraint types
import Base: convert
import IncrementalInference: getSample

##

mutable struct DevelopPrior{T <: SamplableBelief} <: AbstractPrior
  x::T
end
getSample(cf::CalcFactor{<:DevelopPrior}, N::Int=1) = (reshape(rand(cf.factor.x, N),1,N), )

mutable struct DevelopLikelihood{T <: SamplableBelief} <: AbstractRelativeRoots
  x::T
end

getSample(cf::CalcFactor{<:DevelopLikelihood}, N::Int=1) = (reshape(rand(cf.factor.x, N),1,N), )
function (cf::CalcFactor{<:DevelopLikelihood})(meas, wXi, wXj)
  #
  return meas - (wXj - wXi)
end

N  = 100
fg = initfg()

##

@testset "test populate factor graph with a multi-hypothesis factor..." begin

##

v1 = addVariable!(fg, :x1, ContinuousScalar, N=N)

pr = DevelopPrior(Normal(10.0,1.0))
f1 = addFactor!(fg,[:x1],pr)

ensureAllInitialized!(fg)

pts = evalFactor(fg, f1, v1.label, N=N)

@test sum(abs.(pts .- 1.0) .< 5) < 30
@test sum(abs.(pts .- 10.0) .< 5) > 30

v2 = addVariable!(fg, :x2, ContinuousScalar, N=N)
pp = DevelopLikelihood(Normal(100.0,1.0))
f2 = addFactor!(fg, [:x1;:x2], pp)

ensureAllInitialized!(fg)

@test abs(Statistics.mean(getVal(fg, :x2))-110.0) < 10.0

v3 = addVariable!(fg, :x3, ContinuousScalar, N=N)
v4 = addVariable!(fg, :x4, ContinuousScalar, N=N)

ppMH = DevelopLikelihood(Normal(90.0,1.0))
f3 = addFactor!(fg, [:x2;:x3;:x4], ppMH, multihypo=[1.0;0.5;0.5])

# @test IIIF._getCCW(f3).hypoverts == [:x3, :x4]
@test sum(abs.(IIF._getCCW(f3).hypotheses.p[1] .- 0.0)) < 0.1  # 1.0 becomes 0.0 for computational convenience
@test sum(abs.(IIF._getCCW(f3).hypotheses.p[2:3] .- 0.5)) < 0.1

initManual!(fg, :x2, 1*ones(1,100))
initManual!(fg, :x3, 2*ones(1,100))
initManual!(fg, :x4, 3*ones(1,100))

##

end

@testset "Test multi-hypothesis factor convolution exploration" begin

##

pts = approxConv(fg, :x2x3x4f1, :x2, N=N)

@test 99 < sum(pts .<= -70.0)

pts = approxConv(fg, :x2x3x4f1, :x3, N=N)

@test 15 < sum(pts .== 3.0) < 75

pts = approxConv(fg, :x2x3x4f1, :x4, N=N)

@test 15 < sum(pts .== 2.0) < 75

##

end

##

println("Packing converters")

mutable struct PackedDevelopPrior <: PackedInferenceType
  x::String
  PackedDevelopPrior() = new()
  PackedDevelopPrior(x) = new(x)
end
function convert(::Type{PackedDevelopPrior}, d::DevelopPrior)
  PackedDevelopPrior(convert(PackedSamplableBelief, d.x))
end
function convert(::Type{DevelopPrior}, d::PackedDevelopPrior)
  DevelopPrior(convert(SamplableBelief, d.x))
end

mutable struct PackedDevelopLikelihood <: PackedInferenceType
  x::String
  PackedDevelopLikelihood() = new()
  PackedDevelopLikelihood(x) = new(x)
end
function convert(::Type{PackedDevelopLikelihood}, d::DevelopLikelihood)
  PackedDevelopLikelihood(convert(PackedSamplableBelief, d.x))
end
function convert(::Type{DevelopLikelihood}, d::PackedDevelopLikelihood)
  DevelopLikelihood(convert(SamplableBelief, d.x))
end

##

@testset "test packing and unpacking the data structure" begin

##

topack = getSolverData(getFactor(fg,:x1f1))
dd = convert(PackedFunctionNodeData{PackedDevelopPrior},topack)
unpacked = convert(FunctionNodeData{CommonConvWrapper{DevelopPrior}},dd)

@test abs(IIF._getCCW(unpacked).usrfnc!.x.μ - 10.0) < 1e-10
@test abs(IIF._getCCW(unpacked).usrfnc!.x.σ - 1.0) < 1e-10

fct = getFactor(fg, :x2x3x4f1)
@show typeof(fct)
topack = getSolverData(fct) # f3
dd = convert(PackedFunctionNodeData{PackedDevelopLikelihood},topack)
unpacked = convert(FunctionNodeData{CommonConvWrapper{DevelopLikelihood}},dd)

# @test IIF._getCCW(unpacked).hypoverts == Symbol[:x3; :x4]
@test sum(abs.(IIF._getCCW(unpacked).hypotheses.p[1] .- 0.0)) < 0.1
@test sum(abs.(IIF._getCCW(unpacked).hypotheses.p[2:3] .- 0.5)) < 0.1

##

end

##

# start a new factor graph
N = 200
fg = initfg()

##

@testset "test tri-modal factor..." begin

##

v1 = addVariable!(fg, :x1, ContinuousScalar, N=N)

pr = DevelopPrior(Normal(10.0,1.0))
f1 = addFactor!(fg,[:x1],pr)

ensureAllInitialized!(fg)

# Juno.breakpoint("/home/dehann/.julia/v0.5/IncrementalInference/src/ApproxConv.jl",121)

pts = approxConv(fg, Symbol(f1.label), :x1, N=N)

@test sum(abs.(pts .- 1.0) .< 5) < 30
@test sum(abs.(pts .- 10.0) .< 5) > 30

v2 = addVariable!(fg, :x2, ContinuousScalar, N=N)
pp = DevelopLikelihood(Normal(100.0,1.0))
f2 = addFactor!(fg, [:x1;:x2], pp)

ensureAllInitialized!(fg)

@test abs(Statistics.mean(getVal(fg, :x2))-110.0) < 10.0

v3 = addVariable!(fg, :x3, ContinuousScalar, N=N)
v4 = addVariable!(fg, :x4, ContinuousScalar, N=N)
v5 = addVariable!(fg, :x5, ContinuousScalar, N=N)

ppMH = DevelopLikelihood(Normal(90.0,1.0))
f3 = addFactor!(fg, [:x2;:x3;:x4;:x5], ppMH, multihypo=[1.0,0.333,0.333,0.334])

# @test IIF._getCCW(f3).hypoverts == [:x3, :x4]
@test sum(abs.(IIF._getCCW(f3).hypotheses.p[1] .- 0.0)) < 0.1  # 1.0 becomes 0.0 for computational convenience
@test sum(abs.(IIF._getCCW(f3).hypotheses.p[2] .- 0.333)) < 0.001
@test sum(abs.(IIF._getCCW(f3).hypotheses.p[3] .- 0.333)) < 0.001
@test sum(abs.(IIF._getCCW(f3).hypotheses.p[4] .- 0.334)) < 0.001

initManual!(fg, :x2 ,1*ones(1,100))
initManual!(fg, :x3 ,2*ones(1,100))
initManual!(fg, :x4 ,3*ones(1,100))
initManual!(fg, :x5 ,4*ones(1,100))

# solve for certain idx
pts = approxConv(fg, :x2x3x4x5f1, :x2, N=N)

@test 0.95*N < sum(pts .<= -70.0)

# solve for one of uncertain variables
pts = approxConv(fg, :x2x3x4x5f1, :x3, N=N)

@test 0.1*N < sum(80 .< pts .< 100.0) < 0.5*N
@test 0.1*N < sum(pts .== 3.0) < 0.5*N
@test 0.1*N < sum(pts .== 4.0) < 0.5*N

# 0.7 to accomodate bad-init null hypo
@test 0.5*N <= sum(70 .< pts .< 110.0) + sum(pts .== 3.0) + sum(pts .== 4.0)

# solve for one of uncertain variables
pts = approxConv(fg, :x2x3x4x5f1, :x4, N=N)

@test 0.1*N < sum(80 .< pts .< 100.0) < 0.5*N
@test 0.1*N < sum(pts .== 2.0) < 0.5*N
@test 0.1*N < sum(pts .== 4.0) < 0.5*N

@test 0.5*N <= sum(80 .< pts .< 100.0) + sum(pts .== 2.0) + sum(pts .== 4.0)

# solve for one of uncertain variables
pts = approxConv(fg, :x2x3x4x5f1, :x5, N=N)

@test 0.1*N < sum(80 .< pts .< 100.0) < 0.5*N
@test 0.1*N < sum(pts .== 2.0) < 0.5*N
@test 0.1*N < sum(pts .== 3.0) < 0.5*N

@test 0.5*N <= sum(80 .< pts .< 100.0) + sum(pts .== 2.0) + sum(pts .== 3.0)

##

end

@testset "test multihypo api numerical tolerance, #1086" begin

##

fg = initfg()

addVariable!(fg, :x0, ContinuousEuclid{1})
addVariable!(fg, :x1a, ContinuousEuclid{1})
addVariable!(fg, :x1b, ContinuousEuclid{1})
addFactor!(fg, [:x0;:x1a;:x1b], LinearRelative(Normal()), multihypo=[1; 0.5;0.4999999999999])
addFactor!(fg, [:x0;:x1a;:x1b], LinearRelative(Normal()), multihypo=[1; 0.5;0.5000000000001])

##

end

#