text:
#-----------------------------------------------------------------------
#       DataType TEST
#-----------------------------------------------------------------------

tpfl    = Float64
par     = (1.4365276211950395e9,1.4365277e9,6e9,1e-14,1e-10)
N       = 10
np      = 5
X       = [rand(tpfl,4,np) for _ in 1:N]
Xtar    = [rand(tpfl,4) for _ in 1:N]

tcR =  squirrel.seval.TestCases(  par , N , np , X , Xtar )

@test tcR.par   == par
@test tcR.N     == N
@test tcR.ne    == np
@test tcR.X     == X
@test tcR.Xtar  == Xtar

Xc      = [rand(tpfl,4) for _ in 1:N]
Xsc     = [rand(tpfl,4) for _ in 1:N]
erh     = rand(tpfl,N)
erv     = rand(tpfl,N)
err     = rand(tpfl,N)
erhC    = rand(tpfl,N)
ervC    = rand(tpfl,N)
errC    = rand(tpfl,N)
Xc2     = [rand(tpfl,4) for _ in 1:N]
Xsc2    = [rand(tpfl,4) for _ in 1:N]
erh2    = rand(tpfl,N)
erv2    = rand(tpfl,N)
err2    = rand(tpfl,N)
erhC2   = rand(tpfl,N)
ervC2   = rand(tpfl,N)
errC2   = rand(tpfl,N)

tdR =  squirrel.seval.TestData( par , N , X , Xtar , Xc , Xsc , erh ,
                erv , err , erhC , ervC , errC , Xc2 , Xsc2 , erh2 ,
                erv2 , err2 , erhC2 , ervC2 , errC2 )

@test tdR.par   == par
@test tdR.N     == N
@test tdR.X     == X
@test tdR.Xtar  == Xtar
@test tdR.Xc    == Xc
@test tdR.Xsc   == Xsc
@test tdR.erh   == erh
@test tdR.erv   == erv
@test tdR.err   == err
@test tdR.erhC  == erhC
@test tdR.ervC  == ervC
@test tdR.errC  == errC
@test tdR.Xc2   == Xc2
@test tdR.Xsc2  == Xsc2
@test tdR.erh2  == erh2
@test tdR.erv2  == erv2
@test tdR.err2  == err2
@test tdR.erhC2 == erhC2
@test tdR.ervC2 == ervC2
@test tdR.errC2 == errC2

#-----------------------------------------------------------------------
#       DataType converter TEST
#-----------------------------------------------------------------------

tcRD = squirrel.seval.tcfl( tcR , Double64 )

@test typeof(tcRD.X[1][1,1])  == Double64
@test typeof(tcRD.Xtar[1][1]) == Double64

tcT = squirrel.seval.tc2tup( tcR )

@test tcT == (par,N,np,X,Xtar)
@test tcT == (tcR.par,tcR.N,tcR.ne,tcR.X,tcR.Xtar)

tdT = squirrel.seval.td2tup( tdR )

@test tdT == (par,N,X,Xtar,Xc,Xsc,erh,erv,err,erhC,ervC,errC,Xc2,Xsc2,
              erh2,erv2,err2,erhC2,ervC2,errC2)

@test tdT == (tdR.par,tdR.N,tdR.X,tdR.Xtar,tdR.Xc,tdR.Xsc,tdR.erh,
              tdR.erv,tdR.err,tdR.erhC,tdR.ervC,tdR.errC,tdR.Xc2,
              tdR.Xsc2,tdR.erh2,tdR.erv2,tdR.err2,tdR.erhC2,tdR.ervC2,
              tdR.errC2)

tcRb = squirrel.seval.tup2tc( tcT )

tdRb = squirrel.seval.tup2td( tdT )

@test tcRb == tcR

@test tdRb == tdR

#-----------------------------------------------------------------------
#       roteuler & angvec TEST
#-----------------------------------------------------------------------

Rot = squirrel.seval.roteuler()

v3 = rand(3)

@test (Rot[2]*(Rot[1]*v3)) ≈ v3 atol=1e-14

@test squirrel.seval.angvec( Rot[1]*[1.;0.;0.]  ,
                             Rot[1]*[0.;1.;0.]  ) ≈ π/2
@test squirrel.seval.angvec( Rot[1]*[1.;0.;0.]  ,
                             Rot[1]*[0.;0.;-1.] ) ≈ π/2
@test squirrel.seval.angvec( Rot[1]*[0.;-1.;0.] ,
                             Rot[1]*[0.;0.;1.]  ) ≈ π/2
@test squirrel.seval.angvec( Rot[1]*[1.;0.;0.]  ,
                             Rot[1]*[1.;1.;0.]  ) ≈ π/4

#-----------------------------------------------------------------------
#       angmar TEST
#-----------------------------------------------------------------------
Δψ0 = (π/180)*10

@test squirrel.seval.angmar( 0 , Δψ0 )         == false
@test squirrel.seval.angmar( π/2 , Δψ0 )       == false
@test squirrel.seval.angmar( π/2-Δψ0/2 , Δψ0 ) == false

@test squirrel.seval.angmar( Δψ0/2 , Δψ0 )
@test squirrel.seval.angmar( Δψ0 , Δψ0 )
@test squirrel.seval.angmar( 2*Δψ0 , Δψ0 )

@test squirrel.seval.angmar( π/2-Δψ0 , Δψ0 )
@test squirrel.seval.angmar( π/2-2*Δψ0 , Δψ0 )

#-----------------------------------------------------------------------
#       slchk TEST
#-----------------------------------------------------------------------

YE = idg(4)[1]

@test squirrel.seval.slchk(YE)

YE[:,4] = YE[:,4] + [-10^4;0.;0.;0.]

@test squirrel.seval.slchk(YE) == false

#-----------------------------------------------------------------------
#       vrgen TEST
#-----------------------------------------------------------------------

u = rand(3)
sc = 100.
v = squirrel.seval.vrgen( sc , Δψ0 , u )

@test sqrt(dot(v,v)) ≈ sc atol=5e-14
@test squirrel.seval.angmar(squirrel.seval.angvec(u,v),Δψ0)

#-----------------------------------------------------------------------
#       λiRscalc TEST
#-----------------------------------------------------------------------

v = u/sqrt(dot(u,u))

u = rand(3)
x = u/sqrt(dot(u,u))

@test squirrel.seval.λiRscalc( x , v , 0  , 100 ) ≈ 100 atol=5e-14
@test squirrel.seval.λiRscalc( x , x , 10 , 100 ) ≈ 90  atol=5e-14
@test squirrel.seval.λiRscalc( x , x , 50 , 100 ) ≈ 50  atol=5e-14

#-----------------------------------------------------------------------
#       tidc TEST
#-----------------------------------------------------------------------

λ = squirrel.seval.λiRscalc( x , v , 10 , 100 )
v = λ*v
Zi = squirrel.seval.tidc( x , v , λ , η )

@test Zi[2:4] ≈ x
@test Zi[6:8] ≈ v

#-----------------------------------------------------------------------
#       erv & erh TEST
#-----------------------------------------------------------------------

U =  [0.24981918733261255, 0.6256573311421523, 0.8453503802138926,
      0.3093630156499423]

ΔU = [0.000981491116663281, 0.0008969460327706518,
      0.00034983811644343675, 0.00011308545853025453]

@test   squirrel.seval.erV(U,ΔU) ≈
        dot(U[2:4],ΔU[2:4])/norm(U[2:4]) atol=tol

@test   √(squirrel.seval.erV(U,ΔU)^2 + squirrel.seval.erH(U,ΔU)^2) ≈
        norm(ΔU[2:4]) atol=tol

#-----------------------------------------------------------------------
#       pgen TEST
#-----------------------------------------------------------------------

ne    = 6
ntest = 5

P = squirrel.seval.pgen( 6e9 , η , 1e-14 , ne , Δψ0 )

@test size(P[1]) == (4,ne)
@test size(P[2]) == (4,)

#-----------------------------------------------------------------------
#       gen & main TEST
#-----------------------------------------------------------------------

nev = 2

nb	= 24
tol	= 1e-9
ξ	= 1e1

for k=1:10

tc = squirrel.seval.gen(3,η,ne)
typeof(tc) == squirrel.seval.TestCases

td = squirrel.seval.main(tc,squirrel.locator,η,nev,Float64,tol,ξ,nb,ne)

@test length(td.erh) ≈ nev
@test length(td.erv) ≈ nev
@test length(td.err) ≈ nev

Xer = [zeros(4) for _ in 1:nev]

for i=1:nev
    Xer[i][2:4] = (td.Xsc[i][2:4]-td.Xtar[i][2:4])./td.Xtar[i][2:4]
    @test Xer[i] ≈ zeros(4) atol=tol
end

end


text:
module DifferentiablePhysics2D

const DP2D = DifferentiablePhysics2D
export DP2D

import PhysicsPrimitives2D
import PhysicsPrimitives2D: PP2D
import StaticArrays
const SA = StaticArrays
import LinearAlgebra
const LA = LinearAlgebra
import ForwardDiff
const FD = ForwardDiff
import Requires

function __init__()
    Requires.@require Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" include("render.jl")
end

end

ext: jl


text:
# Copyright (C) 2021 Ilya Baran.  This program is distributed under the terms of the MIT license.

import BlackBoxOptim
import NearestNeighbors
import SpecialFunctions
import Statistics
import Distributions
using LinearAlgebra

# Generate plausible population
function generatePPop(blockWithOutputs::AbstractBlock,
    parametersAndBounds::AbstractArray{Tuple{String, Float64, Float64}},
    outputsAndStds::AbstractDict{String, Tuple{S, T}}, # Name to value and std
    count::Int;
    kwargs...) where {S, T}
    generatePPop([blockWithOutputs], parametersAndBounds, outputsAndStds, count; kwargs...)
end

# Generate plausible population with multiple block runs
function generatePPop(blocksWithOutputs::Vector{<:AbstractBlock},
                      parametersAndBounds::AbstractArray{Tuple{String, Float64, Float64}},
                      outputsAndStds::AbstractDict{String, Tuple{S, T}}, # Name to value and std.
                      count::Int;
                      MaxTime = 1000) where {S, T}
    boundBlocks = [BlockWithBindings(block, [name for (name, lower, upper) in parametersAndBounds]) for block in blocksWithOutputs];
    lower = [lower for (_, lower, _) in parametersAndBounds];
    upper = [upper for (_, _, upper) in parametersAndBounds];

    outputs = [computeOutputs(block) for block in blocksWithOutputs];

    (expectedValuesArray, stdsArray) = makeExpectedValuesAndStdsArrays(outputs, outputsAndStds);

    objective = x -> begin
        obj = 0.;
        for i in 1:length(boundBlocks)
            blockCopy = deepcopy(boundBlocks[i]);
            setParameters!(blockCopy, x);
            outputs = computeOutputs(blockCopy);
            outputs = unfold(outputs.values - expectedValuesArray[i]) ./ unfold(stdsArray[i]);
            obj += sum(max.(outputs .^ 2, 1) .- 1)
        end
        return obj
    end

    ppop = Matrix{Float64}(undef, length(lower), count);

    Threads.@threads for i in 1:count
        result = BlackBoxOptim.bboptimize(objective; SearchRange = collect(zip(lower, upper)),
                                          MaxTime = Threads.nthreads() * MaxTime / count, TargetFitness = 0.);
        ppop[:, i] = BlackBoxOptim.best_candidate(result);
    end

    # Print results
    outputResults = Vector(undef, count);
    Threads.@threads for i in 1:count
        allOutputs = "";
        for j in 1:length(boundBlocks)
            blockCopy = deepcopy(boundBlocks[j]);
            setParameters!(blockCopy, ppop[:, i]);
            local outputs = computeOutputs(blockCopy);
            blockOutputs = "";
            for k = 1:length(expectedValuesArray[j])
                v = outputs.values[k];
                min = expectedValuesArray[j][k] - stdsArray[j][k];
                max = expectedValuesArray[j][k] + stdsArray[j][k];
                if min < v < max
                    blockOutputs *= "   $(outputs.variables[k].name): $v in ($min, $max) \n";
                else
                    blockOutputs *= " * $(outputs.variables[k].name): $v OUT ($min, $max) \n";
                end
            end
            allOutputs = allOutputs * blockOutputs;
        end
        outputResults[i] = allOutputs;
    end

    for i in 1:count
        println("Plausible patient $i");
        println(outputResults[i]);
    end

    ppop;
end

# Generate plausible population with multiple block runs
function generatePPopFarthest(blocksWithOutputs::Vector{<:AbstractBlock},
                              parametersAndBounds::AbstractArray{Tuple{String, Float64, Float64}},
                              outputsAndStds::AbstractDict{String, Tuple{S, T}}, # Name to value and std.
                              count::Int;
                              MaxTime = 1000, DistanceFactor = 1.) where {S, T}
    boundBlocks = [BlockWithBindings(block, [name for (name, lower, upper) in parametersAndBounds]) for block in blocksWithOutputs];
    lower = [lower for (_, lower, _) in parametersAndBounds];
    upper = [upper for (_, _, upper) in parametersAndBounds];

    outputs = [computeOutputs(block) for block in blocksWithOutputs];

    (expectedValuesArray, stdsArray) = makeExpectedValuesAndStdsArrays(outputs, outputsAndStds);

    metricPPopSoFar::Vector{Vector{Float64}} = [];
    ppop = Matrix{Float64}(undef, length(lower), count);

    toMetric = p -> begin
        result = deepcopy(p);
        for i in 1:length(p)
            if lower[i] <= 0
                result[i] = (p[i] - lower[i]) / (upper[i] - lower[i]);
            else
                result[i] = log(p[i]);
            end
        end
        return result;
    end

    threads = Threads.nthreads();
    rounds::Integer = ceil(count / threads);

    objective = (x, idx) -> begin
        obj = 0.;
        for i in 1:length(boundBlocks)
            blockCopy = deepcopy(boundBlocks[i]);
            setParameters!(blockCopy, x);
            outputs = computeOutputs(blockCopy);
            outputs = unfold(outputs.values - expectedValuesArray[i]) ./ unfold(stdsArray[i]);
            obj += sum(max.(outputs .^ 2, 1) .- 1);
        end
        if length(metricPPopSoFar) > 0
            m = toMetric(x);
            closest = Inf64;
            for i in 1:length(metricPPopSoFar)
                dist = sum((m - metricPPopSoFar[i]) .^ 2);
                if mod(i + idx, threads) == 1
                    dist = dist * 2;
                end
                closest = min(closest, dist);
            end
            obj -= closest * DistanceFactor;
        end
        return obj;
    end

    for i in 1:rounds
        println("Starting round $i of $rounds");
        batch::Integer = min(threads, count - (i - 1) * threads);
        batchPpop = Matrix{Float64}(undef, length(lower), batch);
        Threads.@threads for j in 1:batch
            result = BlackBoxOptim.bboptimize((x) -> objective(x, j); SearchRange = collect(zip(lower, upper)),
                                              MaxTime = MaxTime / rounds);
            batchPpop[:, j] = BlackBoxOptim.best_candidate(result);
        end
        for j in 1:batch
            push!(metricPPopSoFar, toMetric(batchPpop[:, j]));
            ppop[:, (i - 1) * threads + j] = batchPpop[:, j];
        end
    end

    return ppop;
end

# Take a population and a set of blocks and output ranges.  Optionally filter the population to just those
# whose outputs are in the ranges.  Generate random new patients by sampling in a normal distribution around
# random patients and keeping those that satisfy the output constraints.  The covariance matrix is the
# covariance matrix of the input ppop, scaled by ratio.  We generate and try count new patients (and a subset of
# those are kept).  The return value includes the original (possibly filtered) ppop.
function expandPPop(blocksWithOutputs::Vector{<:AbstractBlock},
                    parametersAndBounds::AbstractArray{Tuple{String, Float64, Float64}},
                    outputsAndStds::AbstractDict{String, Tuple{S, T}}, # Name to value and std.
                    ppop::Matrix{Float64},
                    ratio::Number, # What fraction of the covariance to use
                    count::Int;
                    filter = false) where {S, T}

    boundBlocks = [BlockWithBindings(block, [name for (name, lower, upper) in parametersAndBounds]) for block in blocksWithOutputs];
    lower = [lower for (_, lower, _) in parametersAndBounds];
    upper = [upper for (_, _, upper) in parametersAndBounds];

    outputs = [computeOutputs(block) for block in blocksWithOutputs];

    (expectedValuesArray, stdsArray) = makeExpectedValuesAndStdsArrays(outputs, outputsAndStds);

    check = function(parameters)
        for j in 1:length(boundBlocks)
            blockCopy = deepcopy(boundBlocks[j]);
            setParameters!(blockCopy, parameters);
            outputs = computeOutputs(blockCopy);
            outputs = unfold(outputs.values - expectedValuesArray[j]) ./ unfold(stdsArray[j]);
            if any(abs.(outputs) .> 1)
                return false;
            end
        end
        return true;
    end

    isOk = [true for i = 1:size(ppop, 2)];
    if filter
        Threads.@threads for i in 1:size(ppop, 2)
            isOk[i] = check(ppop[:, i]);
        end
        ppop = ppop[:, isOk];
        println("Kept $(size(ppop, 2)) patients");
    end

    covariance = Statistics.cov(ppop, dims = 2);
    distribution = Distributions.MvNormal(covariance * ratio);

    candidates = Matrix{Float64}(undef, length(lower), count);
    for i in 1:count
        candidates[:, i] = ppop[:, rand(1:size(ppop, 2))] + rand(distribution);
        candidates[:, i] = max.(lower, min.(upper, candidates[:, i]));
    end

    isOk = [true for i = 1:count];

    Threads.@threads for i in 1:count
        isOk[i] = check(candidates[:, i]);
    end

    for i in 1:count
        if isOk[i]
            ppop = hcat(ppop, candidates[:, i]);
        end
    end

    return ppop;
end

function makeExpectedValuesAndStdsArrays(outputs::Vector{Variables}, outputsAndStds::AbstractDict{String, Tuple{S, T}}) where {S, T}
    expectedValuesArray = [];
    stdsArray = [];
    for output in outputs
        expectedValues = [];
        stds = [];
        for variable in output.variables
            push!(expectedValues, outputsAndStds[variable.name][1]);
            push!(stds, outputsAndStds[variable.name][2]);
        end
        push!(expectedValuesArray, expectedValues);
        push!(stdsArray, stds);
    end
    return (expectedValuesArray, stdsArray);
end

function subsamplePPop(ppop::Matrix, block::AbstractBlock, parameterNames::AbstractArray{String}, mean::Vector, covariance::Matrix, count::Integer)
    if length(parameterNames) != 0
        block = BlockWithBindings(block, parameterNames);
    end
    outputs = Array{Variables}(undef, size(ppop, 2));
    for i in 1:length(outputs)
        setParameters!(block, ppop[:, i]);
        outputs[i] = computeOutputs(block);
    end
    outputsMatrix = [outputs[j].values[i] for i=1:length(outputs[1].values), j=1:length(outputs)];
    return subsamplePPop(outputsMatrix, mean, covariance, count);
end

const NEIGHBORS = 5;

# outputs are columns of outputs.  Returns a vector of indices.
function subsamplePPop(outputs::Matrix, mean::Vector, covariance::Matrix, count::Integer)
    tree = NearestNeighbors.KDTree(outputs; leafsize = 10);
    idxs, _ = NearestNeighbors.knn(tree, outputs, NEIGHBORS + 1, true);
    radii = [norm(outputs[idx[end]] - outputs[idx[1]]) for idx in idxs];
    densities = NEIGHBORS ./ nBallVolume(radii, size(outputs, 1));
    scaledProbabilities = normalPDF(outputs, mean, covariance) ./ densities;
    scaledProbabilities = scaledProbabilities / sum(scaledProbabilities);
    return weightedSample(scaledProbabilities, count);
end

function weightedSample(weights::Vector, count::Integer)
    fAndIdx = [(weights[i] < 1e-14 ? -1e100 : log(rand()) / weights[i], i) for i in 1:length(weights)];
    fAndIdx = sort(fAndIdx, rev = true);
    count = min(count, length(weights));
    return sort([s[2] for s in fAndIdx[1:count]]);
end

function normalPDF(points::Matrix, mean::Vector, covariance::Matrix)
    points = points .- mean;
    result = exp.(-0.5 * sum(points .* (inv(covariance) * points), dims=1)[1,:]) / sqrt((2 * pi) ^ length(mean) * abs(det(covariance)));
    return result;
end

function nBallVolume(radius, d::Number)
    (pi ^ (d / 2) / SpecialFunctions.gamma(d / 2 + 1)) * radius .^ d
end

max_stars_repo_name: ilya-baran-personal/modelBlocks

max_issues_repo_name: ilya-baran-personal/modelBlocks

max_forks_repo_name: ilya-baran-personal/modelBlocks


text:
using Pkg
pkg"activate ."

include("Database.jl")
include("Wikipedia.jl")
include("Gameplay.jl")

using .Wikipedia, .Gameplay

articles = newgame(Gameplay.DIFFICULTY_EASY)

for article in articles
  println(article.title)
end

max_stars_repo_name: PacktPublishing/Julia-Projects

max_issues_repo_name: PacktPublishing/Julia-Projects

max_forks_repo_name: PacktPublishing/Julia-Projects


text:
# Run a chain with:
#  1). T_1, T_2, f with the constraint T_1^4(1-f)+f T_2^4 = T_eff^4 within the uncertainties.
#  2). Fit these to the transit transmission data (assuming an achromatic atmosphere).
#  3). See how much this affects the Spitzer IRAC band inferred depths.

# Load in requisite Julia packages:
using DelimitedFiles
using Statistics
using PyPlot
using Printf

# Read in the Flux vs. Teff:
star_grid = readdlm("phot_teff_v02.csv",',',skipstart=1)
# Indices of the "surface" fluxes:
jsurf = collect(7:2:37)  # These bands are: 1 Kepler,2 IRAC Ch1,3 IRAC Ch2,4 SPECULOOS
      # 5 g, 6 r, 7 i, 8 z, 9 y, 10 J, 11 H, 12 Ks, 13 W1, 14 W2, 15 W3, 16 W4.
lam_mod = [ 0.807, 3.6, 4.5, 0.9102, 0.4811, 0.6156, 0.7504, 0.8668, 0.9613,
            1.235, 1.662, 2.159, 3.368, 4.618, 12.082, 22.194]
isort = sortperm(lam_mod)
# Create an effective temperature grid:
tgrid = collect(2000.0:20.0:2980.0)

# Now, read in the transit transmission spectrum for each
# planet from Ducrot et al. (2020):
data_b = readdlm("transmission_spectrum_b.csv",',',skipstart=1)
data_c = readdlm("transmission_spectrum_c.csv",',',skipstart=1)
data_d = readdlm("transmission_spectrum_d.csv",',',skipstart=1)
data_e = readdlm("transmission_spectrum_e.csv",',',skipstart=1)
data_f = readdlm("transmission_spectrum_f.csv",',',skipstart=1)
data_g = readdlm("transmission_spectrum_g.csv",',',skipstart=1)
data_h = readdlm("transmission_spectrum_h.csv",',',skipstart=1)
# Define the spectral indices from Adam's grid for each planet.
# These are the bands for each planet in the transit transmission
# spectra from Elsa's paper:
jb = [1,8,4,10,12,2,3] # K2, LT (z'), SPECULOOS, J, NB2090 (~Ks), Ch1, Ch2
jc = [1,8,4,10,12,2,3]
jd = [1,8,4,10,2,3]
je = [1,8,4,10,2,3]
jf = [1,4,10,2,3]
jg = [1,4,10,2,3]
jh = [1,8,4,2,3]

function compute_tlam(param)
  # Computes the transit wavelength dependence:
  T1,T2,f,db,dc,dd,de,df,dg,dh = param
  # Interpolate the spectra for T1 and T2:
  it1 = ceil(Int64,(T1 - 2000)/980*49)
  it2 = ceil(Int64,(T2 - 2000)/980*49)
  # Now, find interpolated fluxes:
  f1 = star_grid[it1,jsurf]*(tgrid[it1+1]-T1)/20.0+star_grid[it1+1,jsurf]*(T1-tgrid[it1])/20.0
  f2 = star_grid[it2,jsurf]*(tgrid[it2+1]-T2)/20.0+star_grid[it2+1,jsurf]*(T2-tgrid[it2])/20.0
  # Compute the impact on transit depth (equation 1 in Rackham et al.).
  # We are assuming that *all* planets transit the same temperature,
  # so we require that f > 50%; otherwise the spots and surface flip roles:
  if f < 0.5
    tlam = 1 ./(1 .-  f .*(1 .- f2 ./ f1))
  else
    tlam = 1 ./(1 .- (1 - f) .*(1 .- f1 ./ f2))
  end
  return tlam
end

# Define a log likelihood function:
function log_likelihood(param)
  # Ten parameters: two temperatures, spot covering fraction, and seven depths.
  T1,T2,f,db,dc,dd,de,df,dg,dh = param
  # Okay, now interpolate effective temperatures. First check to see if in bounds:
  if (T1 < 2000.0) || (T1 > 2980.0) || (T2 < 2000.0) || (T2 > 2980.0) || (minimum(param) <= 0.0) || (f >= 1.0)
    return -Inf
  end
  # Compute the transit wavelength depedence:
  tlam = compute_tlam(param)
  # Compute the chi-squares:
  chisquare  = sum((data_b[:,4] .- tlam[jb] .* db).^2 ./ data_b[:,5].^2)
  chisquare += sum((data_c[:,4] .- tlam[jc] .* dc).^2 ./ data_c[:,5].^2)
  chisquare += sum((data_d[:,4] .- tlam[jd] .* dd).^2 ./ data_d[:,5].^2)
  chisquare += sum((data_e[:,4] .- tlam[je] .* de).^2 ./ data_e[:,5].^2)
  chisquare += sum((data_f[:,4] .- tlam[jf] .* df).^2 ./ data_f[:,5].^2)
  chisquare += sum((data_g[:,4] .- tlam[jg] .* dg).^2 ./ data_g[:,5].^2)
  chisquare += sum((data_h[:,4] .- tlam[jh] .* dh).^2 ./ data_h[:,5].^2)
  # That's it!
  # Add to it the constraint on Teff:
  chisquare += ((T1^4*(1-f)+T2^4*f)^0.25 - 2566)^2/25^2
  # Return the log likelihood:
  return -0.5*chisquare
end

# Compute the correction factors in the two IRAC bands:
function irac_correction(param)
  T1,T2,f,_ = param
# Computes the correction in the IRAC bands:
  if (T1 < 2000.0) || (T1 > 2980.0) || (T2 < 2000.0) || (T2 > 2980.0) || (minimum(param) <= 0.0) || (f >= 1.0)
    return -Inf
  end
  # Interpolate the spectra for T1 and T2:
  it1 = ceil(Int64,(T1 - 2000)/980*49)
  it2 = ceil(Int64,(T2 - 2000)/980*49)
  # Now, find interpolated fluxes:
  f1 = star_grid[it1,jsurf[2:3]]*(tgrid[it1+1]-T1)/20.0+star_grid[it1+1,jsurf[2:3]]*(T1-tgrid[it1])/20.0
  f2 = star_grid[it2,jsurf[2:3]]*(tgrid[it2+1]-T2)/20.0+star_grid[it2+1,jsurf[2:3]]*(T2-tgrid[it2])/20.0
  if f < 0.5
    return 1 ./(1 .- f .*(1 .- f2 ./ f1))
  else
    return 1 ./(1 .- (1-f) .*(1 .- f1 ./ f2))
  end
end

# Okay, now run a Markov chain:
using AffineInvariantMCMC

# Set up initial conditions for the chain:
nparam  = 10
nwalker = 30
astep   = 2.0
xstart = zeros(nparam,nwalker)
for i=1:nwalker
  param = zeros(nparam)
  param[1] = 2000.0 + 980.0 * rand()
  param[2] = 2000.0 + 980.0 * rand()
  param[3:10] = rand(8)
  xstart[:,i] = param
end

# Run the chain:
#nsample = 100000
#nthin   = 100
nsample = 10000
nthin   = 10
# First run burn-in for 2000 steps:
chain,llhoodvals = AffineInvariantMCMC.sample(log_likelihood,nwalker,xstart,2000,nthin,astep)
# Okay, now run the chain:
chain,llhoodvals = AffineInvariantMCMC.sample(log_likelihood,nwalker,chain[:,:,end],nsample,nthin,astep)
# Now flatten the chains:
flatchain, flatllhoodvals = AffineInvariantMCMC.flattenmcmcarray(chain, llhoodvals)
nflat = nwalker*div(nsample,nthin)
# Now, compute the IRAC band corrections:
irac_grid = zeros(2,nflat)
for i=1:nflat
  irac_grid[:,i] = irac_correction(flatchain[:,i])
end

# Print the results:
println("IRAC Ch1: ",@sprintf("%6.4f",mean(sqrt.(irac_grid[1,:]))),"+-",@sprintf("%6.4f",std(sqrt.(irac_grid[1,:]))))
println("IRAC Ch2: ",@sprintf("%6.4f",mean(sqrt.(irac_grid[2,:]))),"+-",@sprintf("%6.4f",std(sqrt.(irac_grid[2,:]))))

# Now, plot some results:
fig,axes = subplots(2,4)
ax = axes[1]
ax.errorbar(data_b[:,3],data_b[:,4],yerr=data_b[:,5],fmt="o",label="b")
ax.legend()
ax = axes[2]
ax.errorbar(data_c[:,3],data_c[:,4],yerr=data_c[:,5],fmt="o",label="c")
ax.legend()
ax = axes[3]
ax.errorbar(data_d[:,3],data_d[:,4],yerr=data_d[:,5],fmt="o",label="d")
ax.legend()
ax = axes[4]
ax.errorbar(data_e[:,3],data_e[:,4],yerr=data_e[:,5],fmt="o",label="e")
ax.legend()
ax = axes[5]
ax.errorbar(data_f[:,3],data_f[:,4],yerr=data_f[:,5],fmt="o",label="f")
ax.legend()
ax = axes[6]
ax.errorbar(data_g[:,3],data_g[:,4],yerr=data_g[:,5],fmt="o",label="g")
ax.legend()
ax = axes[7]
ax.errorbar(data_h[:,3],data_h[:,4],yerr=data_h[:,5],fmt="o",label="h")
ax.legend()

dsamp = zeros(nflat,7,16)
for i=1:nflat
  # Randomly choose posterior values:
  param = flatchain[:,i]
  tlam_samp = compute_tlam(param)
  for ip=1:7
     dsamp[i,ip,:] =  tlam_samp .* param[3+ip]
  end
end
for ip=1:7
  ax = axes[ip]
  depth_mean = zeros(16)
  depth_sig = zeros(16)
  for j=1:16
    depth_mean[j] = mean(dsamp[:,ip,j])
    depth_sig[j] = std(dsamp[:,ip,j])
  end
  ax.semilogx(lam_mod[isort],depth_mean[isort],color="g")
  ax.semilogx(lam_mod[isort],depth_mean[isort]-depth_sig[isort],color="g",linestyle=":")
  ax.semilogx(lam_mod[isort],depth_mean[isort]+depth_sig[isort],color="g",linestyle=":")
  ax.semilogx(lam_mod[isort],zeros(16) .+ mean(flatchain[3+ip,:]),color="r")
  ax.semilogx(lam_mod[isort],zeros(16) .+ (mean(flatchain[3+ip,:]) - std(flatchain[3+ip,:])),color="r",linestyle=":")
  ax.semilogx(lam_mod[isort],zeros(16) .+ (mean(flatchain[3+ip,:]) + std(flatchain[3+ip,:])),color="r",linestyle=":")
  ax.set_xlabel(L"$\lambda$ [micron]")
  ax.set_ylabel("Depth [%]")
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/domain_reduction.jl
# Contains subroutines used for domain reduction.
#############################################################################
"""
    trivial_filtering!

Excludes OBBT on variable indices that are tight for the solution of the relaxation.
"""
function trivial_filtering!(m::GlobalOptimizer)

    m._preprocess_termination_status = MOI.get(m.relaxed_optimizer, MOI.TerminationStatus())
    m._preprocess_result_status = MOI.get(m.relaxed_optimizer, MOI.PrimalStatus())
    valid_flag, feasible_flag = is_globally_optimal(m._preprocess_termination_status,
                                                    m._preprocess_result_status)

    if valid_flag && feasible_flag
        @inbounds for j = 1:length(m._obbt_working_lower_index)
            if m._obbt_working_lower_index[j]
                vi = m._relaxed_variable_index[j]
                diff = MOI.get(m.relaxed_optimizer, MOI.VariablePrimal(), vi)
                diff -= _lower_bound(FullVar, m, j)
                if abs(diff) <= _obbt_tolerance(m)
                    m._obbt_working_lower_index[j] = false
                end
            end
        end
        @inbounds for j = 1:length(m._obbt_working_upper_index)
            if m._obbt_working_upper_index[j]
                vi = m._relaxed_variable_index[j]
                diff = -MOI.get(m.relaxed_optimizer, MOI.VariablePrimal(), vi)
                diff += _upper_bound(FullVar, m, j)
                if abs(diff) <= _obbt_tolerance(m)
                    m._obbt_working_upper_index[j] = false
                end
            end
        end
    end

    return
end

"""
    bool_indx_diff

Utility function used to set vector of booleans z to x & ~y. Avoids the
generation of conversion of the BitArray created by broadcasting logical operators.
"""
function bool_indx_diff(z::Vector{Bool}, x::Vector{Bool}, y::Vector{Bool})
    for i = 1:length(z)
        @inbounds z[i] = (x[i] & ~y[i])
    end
    return
end

"""
    aggressive_filtering!

Excludes OBBT on variable indices after a search in a filtering direction.
"""
function aggressive_filtering!(m::GlobalOptimizer)

    # Initial filtering vector (negative one direction per remark in Gleixner2017)
    variable_number = m._working_problem._variable_num
    v = -ones(variable_number)

    # Copy prior index set (ignores linear and binary terms)
    obbt_variable_num = m._obbt_variable_num
    copyto!(m._old_low_index, m._obbt_working_lower_index)
    copyto!(m._old_upp_index, m._obbt_working_upper_index)
    copyto!(m._new_low_index, m._obbt_working_lower_index)
    copyto!(m._new_upp_index, m._obbt_working_upper_index)

    # Exclude unbounded directions
    for i = 1:obbt_variable_num
        m._new_low_index[i] &= !(_lower_bound(FullVar, m, i) == -Inf)
        m._new_upp_index[i] &= !(_upper_bound(FullVar, m, i) == Inf)
    end

    # Begin the main algorithm
    for k = 1:m.obbt_aggressive_max_iteration

        # Set index differences and vector for filtering direction
        bool_indx_diff(m._lower_indx_diff, m._old_low_index, m._new_low_index)
        bool_indx_diff(m._upper_indx_diff, m._old_upp_index, m._new_upp_index)

        for i = 1:obbt_variable_num
            if m._lower_indx_diff[i] && (v[i] < 0.0)
                v[i] = 0.0
            elseif m._upper_indx_diff[i] && (v[i] > 0.0)
                v[i] = 0.0
            end
        end

        # Termination Condition
        ((~any(m._new_low_index) & ~any(m._new_upp_index)) || (iszero(v))) && break
        if k >= 2
            if (count(m._lower_indx_diff) + count(m._upper_indx_diff)) < m.obbt_aggressive_min_dimension
                break
            end
        end

        # Set objective in OBBT problem to filtering vector
        MOI.set(m.relaxed_optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)
        saf = SAF(SAT.(v, m._relaxed_variable_index), 0.0)
        MOI.set(m.relaxed_optimizer, MOI.ObjectiveFunction{SAF}(), saf)

        # Optimizes the problem and if successful filter additional bounds
        MOI.optimize!(m.relaxed_optimizer)

        m._preprocess_termination_status = MOI.get(m.relaxed_optimizer, MOI.TerminationStatus())
        m._preprocess_result_status = MOI.get(m.relaxed_optimizer, MOI.PrimalStatus())
        valid_flag, feasible_flag = is_globally_optimal(m._preprocess_termination_status,
                                                        m._preprocess_result_status)

        if valid_flag && feasible_flag
            variable_primal = MOI.get(m.relaxed_optimizer, MOI.VariablePrimal(), m._relaxed_variable_index)
            copyto!(m._new_low_index, m._old_low_index)
            copyto!(m._new_upp_index, m._old_upp_index)
            for i = 1:obbt_variable_num
                vp_value = variable_primal[i]
                if m._old_low_index[i] && vp_value == m.lower_variable_bound[i]
                    m._new_low_index[i] = false
                end
                if m._old_upp_index[i] && vp_value == m.upper_variable_bound[i]
                    m._new_upp_index[i] = false
                end
            end
        else
            return false
        end
    end
    copyto!(m._obbt_working_lower_index, m._new_low_index)
    copyto!(m._obbt_working_upper_index, m._new_upp_index)
    return true
end

"""
    set_node_flag!
"""
function set_node_flag!(m::GlobalOptimizer)
    foreach(set_node_flag!, m._working_problem._nonlinear_constr)
    return nothing
end

"""
    set_reference_point!
"""
function set_reference_point!(m::GlobalOptimizer)

    evaluator = m._working_problem._relaxed_evaluator
    evaluator_x = evaluator.x
    current_xref = m._current_xref

    new_reference_point = false
    for node_i = 1:m._branch_variable_num
        solution_i = m._branch_to_sol_map[node_i]

        node_x = current_xref[node_i]
        solution_x = evaluator_x[solution_i]

        if node_x !== solution_x
            evaluator_x[solution_i] = node_x
            new_reference_point = true
        end
    end

    if new_reference_point
        foreach(x -> (x.has_value = false;), m._working_problem._nonlinear_constr)
    end
    fill!(evaluator.subexpressions_eval, false)

    return nothing
end

"""
    obbt!

Performs OBBT with filtering and greedy ordering as detailed in:
Gleixner, A.M., Berthold, T., Müller, B. et al. J Glob Optim (2017) 67: 731.
https://doi.org/10.1007/s10898-016-0450-4
"""
function obbt!(m::GlobalOptimizer)

    feasibility = true

    # set node and reference point if necessary then solve initial problem to
    # feasibility. This is repeated `obbt_repetitions` number of times in the
    # following fashion. Relax the problem, populate affine constraints, run
    # obbt which contracts variable bounds, delete affine constraints...
    # update variable bounds and repeat. TODO: Keep track of which variables
    # participate in which functions and only delete a constraint if a
    # variable participating in a nonlinear term changes it bounds.
    if m._obbt_performed_flag
        reset_relaxation!(m)
        set_first_relax_point!(m)
    end

    update_relaxed_problem_box!(m)
    if m._nonlinear_evaluator_created
        set_node!(m._working_problem._relaxed_evaluator, n)
        set_node_flag!(m)
        set_reference_point!(m)
    end
    relax_constraints!(m, 1)
    MOI.set(relaxed_optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(relaxed_optimizer)

    # Sets indices to attempt OBBT on
    obbt_variable_num = m._obbt_variable_num
    fill!(m._obbt_working_lower_index, true)
    fill!(m._obbt_working_upper_index, true)

    # Filters out any indicies with active bounds on variables
    # determined by solving the feasibility problem
    trivial_filtering!(m, n)

    # Applies an aggressive filter to eliminate indices that
    # cannot be tightened by obbt
    if m._obbt_aggressive_on
        feasibility = aggressive_filtering!(m, n)
    end

    # extracts info from relaxed problem (redundant if aggressive obbt is called)
    m._preprocess_termination_status = MOI.get(relaxed_optimizer, MOI.TerminationStatus())
    m._preprocess_result_status = MOI.get(relaxed_optimizer, MOI.PrimalStatus())
    valid_flag, feasible_flag = is_globally_optimal(m._preprocess_termination_status,
                                                    m._preprocess_result_status)

    if valid_flag && feasible_flag
        xLP = MOI.get(relaxed_optimizer, MOI.VariablePrimal(), m._relaxed_variable_index)
    else
        return false
    end

    # continue tightening bounds by optimization until all indices have been checked
    # or the node is empty and the problem is thus proven infeasible
    while (any(m._obbt_working_lower_index) || any(m._obbt_working_upper_index)) && !isempty(n)

        # Get lower value
        lower_indx = -1;     upper_indx = -1
        lower_value = Inf;   upper_value = Inf

        # min of xLP - yL on active
        if any(m._obbt_working_lower_index)
            for i = 1:obbt_variable_num
                if m._obbt_working_lower_index[i]
                    temp_value = _lower_solution(FullVar, m, i) - _lower_bound(FullVar, m, i)
                    # Need less than or equal to handle unbounded cases
                    if temp_value <= lower_value
                        lower_value = temp_value
                        lower_indx = i
                    end
                end
            end
        end

        # min of yU - xLP on active
        if any(m._obbt_working_upper_index)
            for i = 1:obbt_variable_num
                if m._obbt_working_upper_index[i]
                    temp_value = _upper_bound(FullVar, m, i) - _lower_solution(FullVar, m, i)
                    if temp_value <= upper_value
                        upper_value = temp_value
                        upper_indx = i
                    end
                end
            end
        end

        # default to upper bound if no lower bound is found, use maximum distance otherwise
        if lower_value <= upper_value && lower_indx > 0

            m._obbt_working_lower_index[lower_indx] = false
            var = SV(m._relaxed_variable_index[lower_indx])

            MOI.set(relaxed_optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
            MOI.set(relaxed_optimizer, MOI.ObjectiveFunction{SV}(), var)

            MOI.optimize!(m.relaxed_optimizer)
            m._preprocess_termination_status = MOI.get(relaxed_optimizer, MOI.TerminationStatus())
            m._preprocess_result_status = MOI.get(relaxed_optimizer, MOI.PrimalStatus())
            valid_flag, feasible_flag = is_globally_optimal(m._preprocess_termination_status,
                                                            m._preprocess_result_status)

            if valid_flag && feasible_flag
                xLP .= MOI.get(relaxed_optimizer, MOI.VariablePrimal(), m._relaxed_variable_index)

                node_index = branch_to_sol_map[lower_indx]
                updated_value = xLP[node_index]
                previous_value = n.lower_variable_bound[lower_indx]

                # if bound is improved update node and corresponding constraint update
                # the node bounds and the single variable bound in the relaxation
                # we assume branching does not occur on fixed variables and interval
                # constraints are internally bridged by EAGO. So the only L <= x
                # constraint in the model is a GreaterThan.
                if updated_value > previous_value && (updated_value - previous_value) > 1E-6
                    ci_list = MOI.get(relaxed_optimizer, MOI.ListOfConstraintIndices{SAF,LT}())
                    sv_geq_ci = m._node_to_sv_geq_ci[lower_indx]
                    MOI.set(relaxed_optimizer, MOI.ConstraintSet(), sv_geq_ci, GT(updated_value))
                    n.lower_variable_bound[lower_indx] = updated_value
                end

                if isempty(n)
                    feasibility = false
                    break
                end

            elseif valid_flag && !feasible_flag
                feasibility = false
                break

            else
                break
            end

        elseif upper_indx > 0

            m._obbt_working_upper_index[upper_indx] = false
            var = SV(m._relaxed_variable_index[upper_indx])
            MOI.set(relaxed_optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)
            MOI.set(relaxed_optimizer, MOI.ObjectiveFunction{SV}(), var)
            MOI.optimize!(relaxed_optimizer)
            m._preprocess_termination_status = MOI.get(m.relaxed_optimizer, MOI.TerminationStatus())
            m._preprocess_result_status = MOI.get(m.relaxed_optimizer, MOI.PrimalStatus())
            valid_flag, feasible_flag = is_globally_optimal(m._preprocess_termination_status,
                                                            m._preprocess_result_status)

            if valid_flag && feasible_flag
                xLP .= MOI.get(relaxed_optimizer, MOI.VariablePrimal(), m._relaxed_variable_index)
                node_index = branch_to_sol_map[upper_indx]
                updated_value = xLP[node_index]
                previous_value = n.upper_variable_bound[upper_indx]

                # if bound is improved update node and corresponding constraint update
                # the node bounds and the single variable bound in the relaxation
                # we assume branching does not occur on fixed variables and interval
                # constraints are internally bridged by EAGO. So the only U => x
                # constraint in the model is a LessThan.
                if updated_value < previous_value && (previous_value - updated_value) > 1E-6
                    sv_leq_ci = m._node_to_sv_leq_ci[upper_indx]
                    MOI.set(relaxed_optimizer, MOI.ConstraintSet(), sv_leq_ci, LT(updated_value))
                    n.upper_variable_bound[upper_indx] = updated_value
                end

                if isempty(n)
                    feasibility = false
                    break
                end

            elseif valid_flag && !feasible_flag
                feasibility = false
                break

            else
                break
            end

        else
            break
        end
        trivial_filtering!(m)
    end

    return feasibility
end

"""
    fbbt!

Performs feasibility-based bound tightening on a back-end constraint and returns `true` if it is feasible or
`false` if it is infeasible.
"""
function fbbt! end

function fbbt!(::Type{Val{:leq}}, m::Optimizer, f::AffineFunction)

    # compute full sum
    temp_sum = -f.constant
    @inbounds for k = 1:f.len
        aik, indx_k = f.terms[k]
        if aik !== 0.0
            aik_xL = aik*_lower_bound(FullVar, m, indx_k)
            aik_xU = aik*_upper_bound(FullVar, m, indx_k)
            temp_sum -= min(aik_xL, aik_xU)
        end
    end

    # subtract extra term, check to see if implied bound is better, if so update the node and
    # the working sum if the node is now empty then break
    @inbounds for k = 1:f.len
        aik, indx_k = f.terms[k]
        if aik !== 0.0
            xL = _lower_bound(FullVar, m, indx_k)
            xU = _upper_bound(FullVar, m, indx_k)
            aik_xL = aik*xL
            aik_xU = aik*xU
            temp_sum += min(aik_xL, aik_xU)
            xh = temp_sum/aik

            if aik > 0.0
                (xh < xL) && return false
                if xh > xL
                    _set_upper_bound(FullVar, m, indx_k, xh)
                end
            elseif aik < 0.0
                (xh > xU) && return false
                if xh < xU
                    _set_lower_bound(FullVar, m, indx_k, xh)
                end
            else
                temp_sum -= min(aik_xL, aik_xU)
                continue
            end

            aik_xL = aik*_lower_bound(FullVar, m, indx_k)
            aik_xU = aik*_upper_bound(FullVar, m, indx_k)
            temp_sum -= min(aik_xL, aik_xU)
        end
    end

    return true
end

function fbbt!(::Type{Val{:eq}}, m::Optimizer, f::AffineFunction)

    # compute full sum
    terms = f.terms
    temp_sum_leq = -f.constant
    temp_sum_geq = -f.constant

    @inbounds for k = 1:f.len
        aik, indx_k = terms[k]
        if aik !== 0.0
            aik_xL = aik*_lower_bound(FullVar, m, indx_k)
            aik_xU = aik*_upper_bound(FullVar, m, indx_k)
            temp_sum_leq -= min(aik_xL, aik_xU)
            temp_sum_geq -= max(aik_xL, aik_xU)
        end
    end

    # subtract extra term, check to see if implied bound is better, if so update the node and
    # the working sum if the node is now empty then break
    @inbounds for k = 1:f.len

        aik, indx_k = terms[k]
        if aik !== 0.0

            xL = _lower_bound(FullVar, m, indx_k)
            xU = _upper_bound(FullVar, m, indx_k)

            aik_xL = aik*xL
            aik_xU = aik*xU

            temp_sum_leq += min(aik_xL, aik_xU)
            temp_sum_geq += max(aik_xL, aik_xU)

            xh_leq = temp_sum_leq/aik
            xh_geq = temp_sum_geq/aik

            if aik > 0.0
                (xh_leq < xL) && return false
                if xh_leq > xL
                    upper_bounds[indx_k] = xh_leq
                end
                (xh_geq > xU) && return false
                if (xh_geq < xU)
                    lower_bounds[indx_k] = xh_geq
                end

            elseif aik < 0.0
                (xh_leq > xU) && return false
                if xh_leq < xU
                    lower_bounds[indx_k] = xh_leq
                end
                (xh_geq < xL) && return false
                if (xh_geq > xL)
                    upper_bounds[indx_k] = xh_geq
                end

            else
                temp_sum_leq -= min(aik_xL, aik_xU)
                temp_sum_geq -= max(aik_xL, aik_xU)
                continue

            end
            aik_xL = aik*_lower_bound(FullVar, m, indx_k)
            aik_xU = aik*_upper_bound(FullVar, m, indx_k)

            temp_sum_leq -= min(aik_xL, aik_xU)
            temp_sum_geq -= max(aik_xL, aik_xU)
        end
    end

    return true
end

cp_condition(m::Optimizer) = false

"""
    set_constraint_propagation_fbbt!

Performs bound tightening based on forward/reverse interval and/or McCormick passes. This routine
resets the current node with new interval bounds.
"""
function set_constraint_propagation_fbbt!(m::Optimizer)
    feasible_flag = true

    if m._nonlinear_evaluator_created
        evaluator = m._working_problem._relaxed_evaluator
        set_node!(m._working_problem._relaxed_evaluator, m._current_node)
        set_node_flag!(m)
        set_reference_point!(m)

        m._working_problem._relaxed_evaluator.is_first_eval = m._new_eval_constraint
        for constr in m._working_problem._nonlinear_constr
            if feasible_flag
                forward_pass!(evaluator, constr)
                feasible_flag &= reverse_pass!(evaluator, constr)
                evaluator.interval_intersect = true
            end
        end

        m._working_problem._relaxed_evaluator.is_first_eval = m._new_eval_constraint
        for constr in m._working_problem._nonlinear_constr
            if feasible_flag
                forward_pass!(evaluator, constr)
            end
        end

        evaluator.is_post = m.subgrad_tighten

        m._new_eval_constraint = false
        m._new_eval_objective = false

        retrieve_x!(m._current_xref, evaluator)
        m._current_node = retrieve_node(evaluator)
    end

    return feasible_flag
end


text:
import IntervalArithmetic: ∩, ∪, entireinterval

export ∩, ∪, entireinterval

for f in (:∩, :∪)
    @eval $f(a::NumberInterval, b::NumberInterval) =
        NumberInterval($f(Interval(a), Interval(b)))
end

for f in (:entireinterval, )
    @eval $f(a::NumberInterval) = NumberInterval($f(Interval(a)))
end


text:
using JuMP
using SeaPearl

struct Item
    id      :: Int
    value   :: Int
    weight  :: Int
end

mutable struct Solution
    content     :: AbstractArray{Bool}
    value       :: Int
    weight      :: Int
    optimality  :: Bool
end

struct InputData
    items               :: AbstractArray{Union{Item, Nothing}}
    sortedItems         :: AbstractArray{Union{Item, Nothing}}
    numberOfItems       :: Int
    capacity            :: Int
end

include("IOmanager.jl")

struct KnapsackVariableSelectionJuMP <: SeaPearl.AbstractVariableSelection{true} end
function (::KnapsackVariableSelectionJuMP)(model::SeaPearl.CPModel)
    i = 1
    while SeaPearl.isbound(model.variables[string(i)])
        i += 1
    end
    return model.variables[string(i)]
end

function solve_knapsack_JuMP(filename::String; benchmark=false)
    input = parseFile!(filename)

    permutation = sortperm(input.items; by=(x) -> x.value/x.weight, rev=true)

    n = input.numberOfItems

    model = Model(SeaPearl.Optimizer)

    ### Variable declaration ###
    @variable(model, 0 <= x[1:n] <= 1)

    ### Constraints ###
    @expression(model, weight_sum, input.items[permutation[1]].weight * x[1])
    for i in 2:n
        add_to_expression!(weight_sum, input.items[permutation[i]].weight, x[i])
    end
    @constraint(model, 0 <= weight_sum <= input.capacity)

    ### Objective ### minimize: -sum(v[i]*x_a[i])
    @expression(model, val_sum, input.items[permutation[1]].value * x[1])
    for i in 2:n
        add_to_expression!(val_sum, input.items[permutation[i]].value, x[i])
    end
    @objective(model, Min, -val_sum)

    # define the heuristic used for variable selection
    variableheuristic = KnapsackVariableSelectionJuMP()
    MOI.set(model, SeaPearl.MOIVariableSelectionAttribute(), variableheuristic)

    optimize!(model)
    status = MOI.get(model, MOI.TerminationStatus())

    println(model)
    println(status)
    println(has_values(model))
    println(value.(x))
    println(solutionFromJuMPWithoutDP(value.(x), input, permutation))
end

function solutionFromSeaPearl(SeaPearlSol::SeaPearl.Solution, input::InputData, permutation::Array{Int})
    taken = falses(input.numberOfItems)
    value = 0
    weight = 0
    for i in 1:input.numberOfItems
        if haskey(SeaPearlSol, "x_a[" * string(i) * "]")
            taken[permutation[i]] = convert(Bool, SeaPearlSol["x_a[" * string(i) * "]"])
            if taken[permutation[i]]
                value += input.items[permutation[i]].value
                weight += input.items[permutation[i]].weight
            end
        end
    end
    return Solution(taken, value, weight, false)
end

function solutionFromJuMPWithoutDP(jump_sol, input::InputData, permutation::Array{Int})
    taken = falses(input.numberOfItems)
    value = 0
    weight = 0
    for i in 1:input.numberOfItems
        taken[permutation[i]] = convert(Bool, jump_sol[i])
        if taken[permutation[i]]
            value += input.items[permutation[i]].value
            weight += input.items[permutation[i]].weight
        end
    end
    return Solution(taken, value, weight, false)
end

max_stars_repo_name: corail-research/SeaPearlZoo

max_issues_repo_name: corail-research/SeaPearlZoo

max_forks_repo_name: corail-research/SeaPearlZoo


text:
using SeisReconstruction
using LinearAlgebra
using Test

# test of Conjugate Gradients
nd = 1000
nm = 300
L = randn(nd,nm)
m = randn(nm)
d = L*m
d = d + 0.2*randn(nd)

# closed form solution
mu = 0.2
m1 = (L'*L + mu*Matrix(I,nm,nm))\(L'*d)

m0 = zeros(nm)
m2,cost = ConjugateGradients(d,[MatrixMultiplyOp],[Dict(:matrix=>L)],Niter=nm,mu=mu)

# test that quality factor between CG and closed form solution
# is greater than 50 Decibels
quality_factor = 10*log10(norm(m1[:],2)/norm(m2[:]-m1[:],2))
println("Quality factor = ",quality_factor)
@test quality_factor > 10.


text:
mutable struct AdvancedMatrixOperators

    bmo       :: BasicMatrixOperators

    T_DIVx_U    :: AbstractArray{Float64, 2}
    T_DIVy_V    :: AbstractArray{Float64, 2}
    T_DIVz_W    :: AbstractArray{Float64, 2}

    U_∂x_T      :: AbstractArray{Float64, 2}
    V_∂y_T      :: AbstractArray{Float64, 2}
    W_∂z_T      :: AbstractArray{Float64, 2}

    T_∂x_U      :: AbstractArray{Float64, 2}
    T_∂y_V      :: AbstractArray{Float64, 2}
    T_∂z_W      :: AbstractArray{Float64, 2}

    U_∂y_UV      :: AbstractArray{Float64, 2}
    V_∂x_UV      :: AbstractArray{Float64, 2}

    U_interp_T :: AbstractArray{Float64, 2}  # interpolation of T grid onto U grid
    V_interp_T :: AbstractArray{Float64, 2}  # interpolation of T grid onto V grid
    W_interp_T :: AbstractArray{Float64, 2}  # interpolation of T grid onto W grid

    T_interp_UV :: AbstractArray{Float64, 2}
    U_interp_UV :: AbstractArray{Float64, 2}
    V_interp_UV :: AbstractArray{Float64, 2}

    UV_interp_U :: AbstractArray{Float64, 2}
    UV_interp_V :: AbstractArray{Float64, 2}

    T_mask_T       :: AbstractArray{Float64, 2}
    U_mask_U       :: AbstractArray{Float64, 2}
    V_mask_V       :: AbstractArray{Float64, 2}
    W_mask_W       :: AbstractArray{Float64, 2}
    T_bordermask_T :: AbstractArray{Float64, 2}

    T_Δx_T :: AbstractArray{Float64, 2}
    T_Δy_T :: AbstractArray{Float64, 2}
    T_Δz_T :: AbstractArray{Float64, 2}

    U_Δx_U :: AbstractArray{Float64, 2}
    U_Δy_U :: AbstractArray{Float64, 2}
    U_Δz_U :: AbstractArray{Float64, 2}

    V_Δx_V :: AbstractArray{Float64, 2}
    V_Δy_V :: AbstractArray{Float64, 2}
    V_Δz_V :: AbstractArray{Float64, 2}

    W_Δx_W :: AbstractArray{Float64, 2}
    W_Δy_W :: AbstractArray{Float64, 2}
    W_Δz_W :: AbstractArray{Float64, 2}

    UV_Δx_UV :: AbstractArray{Float64, 2}
    UV_Δy_UV :: AbstractArray{Float64, 2}
    UV_Δz_UV :: AbstractArray{Float64, 2}

    T_invΔx_T :: AbstractArray{Float64, 2}
    T_invΔy_T :: AbstractArray{Float64, 2}
    T_invΔz_T :: AbstractArray{Float64, 2}

    U_invΔx_U :: AbstractArray{Float64, 2}
    U_invΔy_U :: AbstractArray{Float64, 2}
    U_invΔz_U :: AbstractArray{Float64, 2}

    V_invΔx_V :: AbstractArray{Float64, 2}
    V_invΔy_V :: AbstractArray{Float64, 2}
    V_invΔz_V :: AbstractArray{Float64, 2}

    W_invΔx_W :: AbstractArray{Float64, 2}
    W_invΔy_W :: AbstractArray{Float64, 2}
    W_invΔz_W :: AbstractArray{Float64, 2}

    UV_invΔx_UV :: AbstractArray{Float64, 2}
    UV_invΔy_UV :: AbstractArray{Float64, 2}
    UV_invΔz_UV :: AbstractArray{Float64, 2}

    T_Δvol_T    :: AbstractArray{Float64, 2}
    T_invΔvol_T :: AbstractArray{Float64, 2}

    function AdvancedMatrixOperators(;
        gd             :: PolelikeCoordinate.Grid,
        mask_T         :: AbstractArray{Float64, 3},
        bmo :: Union{Nothing, BasicMatrixOperators} = nothing,
    )

        # define a converter to make 2D variable repeat in z direction for Nz times
        cvt3_diagm = (x,) -> spdiagm( 0 => view(x, :) )

        Nx = gd.Nx
        Ny = gd.Ny
        Nz = gd.Nz
        if bmo == nothing
            println("Construct BMO")
            @time bmo = BasicMatrixOperators(Nx=Nx, Ny=Ny, Nz=Nz)
        end

        if length(mask_T) != bmo.T_pts
            throw(ErrorException("Length of mask_T does not conform"))
        end

        #println("Making masks")
        mask3_flat = view(mask_T,  :)

        onV_if_unblocked_north_onT = bmo.V_S_T  * mask3_flat
        onV_if_unblocked_south_onT = bmo.V_N_T  * mask3_flat
        onU_if_unblocked_east_onT  = bmo.U_W_T  * mask3_flat
        onU_if_unblocked_west_onT  = bmo.U_E_T  * mask3_flat
        onW_if_unblocked_up_onT    = bmo.W_DN_T * mask3_flat
        onW_if_unblocked_dn_onT    = bmo.W_UP_T * mask3_flat

        V_mask = onV_if_unblocked_north_onT .* onV_if_unblocked_south_onT
        U_mask = onU_if_unblocked_east_onT  .* onU_if_unblocked_west_onT
        W_mask = onW_if_unblocked_up_onT    .* onW_if_unblocked_dn_onT

        T_mask_T = spdiagm(0 => mask3_flat)
        V_mask_V = spdiagm(0 => V_mask)
        U_mask_U = spdiagm(0 => U_mask)
        W_mask_W = spdiagm(0 => W_mask)

        T_bordermask_T = spdiagm( 0 => (
               (bmo.T_N_T  * mask3_flat)
            .* (bmo.T_S_T  * mask3_flat)
            .* (bmo.T_E_T  * mask3_flat)
            .* (bmo.T_W_T  * mask3_flat)
            .* (bmo.T_UP_T * mask3_flat)
            .* (bmo.T_DN_T * mask3_flat)
        ))

        #=
        nomotionmask3_flat = view(nomotionmask3,  :)

        onV_if_unblocked_north_onT = bmo.V_S_T  * nomotionmask3_flat
        onV_if_unblocked_south_onT = bmo.V_N_T  * nomotionmask3_flat
        onU_if_unblocked_east_onT  = bmo.U_W_T  * nomotionmask3_flat
        onU_if_unblocked_west_onT  = bmo.U_E_T  * nomotionmask3_flat
        onW_if_unblocked_up_onT    = bmo.W_DN_T * nomotionmask3_flat
        onW_if_unblocked_dn_onT    = bmo.W_UP_T * nomotionmask3_flat

        V_fluxmask = onV_if_unblocked_north_onT .* onV_if_unblocked_south_onT
        U_fluxmask = onU_if_unblocked_east_onT  .* onU_if_unblocked_west_onT
        W_fluxmask = onW_if_unblocked_up_onT    .* onW_if_unblocked_dn_onT

        V_fluxmask_V = V_mask_V * spdiagm(0 => V_fluxmask)
        U_fluxmask_U = U_mask_U * spdiagm(0 => U_fluxmask)
        W_fluxmask_W = W_mask_W * spdiagm(0 => W_fluxmask)
        =#
        #println("Making sides")
        # ===== [ BEGIN face area and lengths on U V ] =====

        T_Δx_T = (gd.Δx_T |>  cvt3_diagm)
        T_Δy_T = (gd.Δy_T |>  cvt3_diagm)
        T_Δz_T = (gd.Δz_T |>  cvt3_diagm)

        U_Δx_U = (gd.Δx_U |>  cvt3_diagm)
        U_Δy_U = (gd.Δy_U |>  cvt3_diagm)
        U_Δz_U = (gd.Δz_U |>  cvt3_diagm)

        V_Δx_V = (gd.Δx_V |>  cvt3_diagm)
        V_Δy_V = (gd.Δy_V |>  cvt3_diagm)
        V_Δz_V = (gd.Δz_V |>  cvt3_diagm)

        W_Δx_W = (gd.Δx_W |>  cvt3_diagm)
        W_Δy_W = (gd.Δy_W |>  cvt3_diagm)
        W_Δz_W = (gd.Δz_W |>  cvt3_diagm)

        UV_Δx_UV = (gd.Δx_UV |>  cvt3_diagm)
        UV_Δy_UV = (gd.Δy_UV |>  cvt3_diagm)
        UV_Δz_UV = (gd.Δz_UV |>  cvt3_diagm)

        T_invΔx_T = (gd.Δx_T.^(-1) |>  cvt3_diagm)
        T_invΔy_T = (gd.Δy_T.^(-1) |>  cvt3_diagm)
        T_invΔz_T = (gd.Δz_T.^(-1) |>  cvt3_diagm)

        U_invΔx_U = (gd.Δx_U.^(-1) |>  cvt3_diagm)
        U_invΔy_U = (gd.Δy_U.^(-1) |>  cvt3_diagm)
        U_invΔz_U = (gd.Δz_U.^(-1) |>  cvt3_diagm)

        V_invΔx_V = (gd.Δx_V.^(-1) |>  cvt3_diagm)
        V_invΔy_V = (gd.Δy_V.^(-1) |>  cvt3_diagm)
        V_invΔz_V = (gd.Δz_V.^(-1) |>  cvt3_diagm)

        W_invΔx_W = (gd.Δx_W.^(-1) |>  cvt3_diagm)
        W_invΔy_W = (gd.Δy_W.^(-1) |>  cvt3_diagm)
        W_invΔz_W = (gd.Δz_W.^(-1) |>  cvt3_diagm)

        UV_invΔx_UV = (gd.Δx_UV.^(-1) |>  cvt3_diagm)
        UV_invΔy_UV = (gd.Δy_UV.^(-1) |>  cvt3_diagm)
        UV_invΔz_UV = (gd.Δz_UV.^(-1) |>  cvt3_diagm)

        T_Δvol_T = T_Δx_T * T_Δy_T * T_Δz_T
        T_invΔvol_T = T_invΔx_T * T_invΔy_T * T_invΔz_T

        # Δz is special. Need to clear NaN
        function clearNaN!(m)
            for i = 1:length(m.nzval)
                if isnan(m.nzval[i])
                    m.nzval[i] = 0
                end
            end
            dropzeros!(m)
        end

        # ===== [ END face area and lengths on T ] =====

        #println("Making derivatives")

        # ===== [ BEG making matrix ] =====
        # MAGIC!!

        T_DIVx_U = T_mask_T * T_invΔvol_T * (bmo.T_W_U  - bmo.T_E_U  ) * U_Δz_U * U_Δy_U  ; dropzeros!(T_DIVx_U);
        T_DIVy_V = T_mask_T * T_invΔvol_T * (bmo.T_S_V  - bmo.T_N_V  ) * V_Δz_V * V_Δx_V  ; dropzeros!(T_DIVy_V);
        T_DIVz_W = T_mask_T * T_invΔvol_T * (bmo.T_DN_W - bmo.T_UP_W ) * W_Δx_W * W_Δy_W  ; dropzeros!(T_DIVz_W);

        U_∂x_T = U_mask_U * U_invΔx_U * (bmo.U_W_T  - bmo.U_E_T)                 ; dropzeros!(U_∂x_T);
        V_∂y_T = V_mask_V * V_invΔy_V * (bmo.V_S_T  - bmo.V_N_T)                 ; dropzeros!(V_∂y_T);
        W_∂z_T = W_mask_W * W_invΔz_W * (bmo.W_DN_T - bmo.W_UP_T)                ; dropzeros!(W_∂z_T);

        T_∂x_U  = T_mask_T * T_invΔx_T * ( bmo.T_W_U - bmo.T_E_U )               ; dropzeros!(T_∂x_U);
        T_∂y_V  = T_mask_T * T_invΔy_T * ( bmo.T_S_V - bmo.T_N_V )               ; dropzeros!(T_∂y_V);
        T_∂z_W  = T_mask_T * T_invΔz_T * ( bmo.T_DN_W - bmo.T_UP_W )             ; dropzeros!(T_∂z_W);

        U_∂y_UV = U_mask_U * U_invΔy_U * (bmo.U_S_UV  - bmo.U_N_UV)              ; dropzeros!(U_∂y_UV);
        V_∂x_UV = V_mask_V * V_invΔx_V * (bmo.V_W_UV  - bmo.V_E_UV)              ; dropzeros!(V_∂x_UV);

        function selfDivision(m, ones_vec)
            local wgts = m * ones_vec
            m_t = transpose(m) |> sparse
            for (i, wgt) in enumerate(wgts)
                if wgt != 0
                    _beg = m_t.colptr[i]
                    _end = m_t.colptr[i+1]-1
                    m_t.nzval[_beg:_end] ./= wgt
                end
            end

            return dropzeros(transpose(m_t) |> sparse)
        end

        #println("Making interpolations part 1")
        ones_T  = ones(Float64, bmo.T_pts)
        ones_U  = ones(Float64, bmo.U_pts)
        ones_V  = ones(Float64, bmo.V_pts)
        ones_UV = ones(Float64, bmo.UV_pts)

        #println("Making interpolations part 2")

        U_interp_T = (bmo.U_W_T + bmo.U_E_T) * T_mask_T
        U_interp_T = selfDivision(U_interp_T, ones_T)

        V_interp_T = (bmo.V_S_T + bmo.V_N_T) * T_mask_T
        V_interp_T = selfDivision(V_interp_T, ones_T)

        W_interp_T = (bmo.W_DN_T + bmo.W_UP_T) * T_mask_T
        W_interp_T = selfDivision(W_interp_T, ones_T)

        # Notice the configuration might make T grid value
        # be produced on a masked grid if it is surrounded
        # by four valid UV pts. So, we need to put a T mask
        # at the end
        T_interp_UV = T_mask_T * (bmo.T_E_U + bmo.T_W_U) * (bmo.U_S_UV + bmo.U_N_UV)
        T_interp_UV = selfDivision(T_interp_UV, ones_UV)

        U_interp_UV = U_mask_U * (bmo.U_S_UV + bmo.U_N_UV)
        U_interp_UV = selfDivision(U_interp_UV, ones_UV)

        V_interp_UV = V_mask_V * (bmo.V_W_UV + bmo.V_E_UV)
        V_interp_UV = selfDivision(V_interp_UV, ones_UV)

        UV_interp_U = (bmo.UV_N_U + bmo.UV_S_U)
        UV_interp_U = selfDivision(UV_interp_U, ones_U)

        UV_interp_V = (bmo.UV_E_V + bmo.UV_W_V)
        UV_interp_V = selfDivision(UV_interp_V, ones_V)

        #U_fluxmask_interp_T = U_fluxmask_U * U_interp_T
        #V_fluxmask_interp_T = V_fluxmask_V * V_interp_T
        #W_fluxmask_interp_T = W_fluxmask_W * W_interp_T

        #println("Sum")

        sumΔvol_T = reshape(ones(Float64, bmo.T_pts), 1, :) * T_mask_T * T_Δvol_T

        #println("finished.")

        return new(

            bmo,

            T_DIVx_U,
            T_DIVy_V,
            T_DIVz_W,

            U_∂x_T,
            V_∂y_T,
            W_∂z_T,

            T_∂x_U,
            T_∂y_V,
            T_∂z_W,

            U_∂y_UV,
            V_∂x_UV,

            U_interp_T,
            V_interp_T,
            W_interp_T,

            T_interp_UV,
            U_interp_UV,
            V_interp_UV,

            UV_interp_U,
            UV_interp_V,

            T_mask_T,
            U_mask_U,
            V_mask_V,
            W_mask_W,
            T_bordermask_T,

            T_Δx_T,
            T_Δy_T,
            T_Δz_T,

            U_Δx_U,
            U_Δy_U,
            U_Δz_U,

            V_Δx_V,
            V_Δy_V,
            V_Δz_V,

            W_Δx_W,
            W_Δy_W,
            W_Δz_W,

            UV_Δx_UV,
            UV_Δy_UV,
            UV_Δz_UV,

            T_invΔx_T,
            T_invΔy_T,
            T_invΔz_T,

            U_invΔx_U,
            U_invΔy_U,
            U_invΔz_U,

            V_invΔx_V,
            V_invΔy_V,
            V_invΔz_V,

            W_invΔx_W,
            W_invΔy_W,
            W_invΔz_W,

            UV_invΔx_UV,
            UV_invΔy_UV,
            UV_invΔz_UV,

            T_Δvol_T,
            T_invΔvol_T,

        )
    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
function train!(;
    model, ds_train, ds_test, ds_validation,
    loss, optimizer, measures, recorded_measures,
    max_epochs, min_epochs, early_stopping_n, early_stopping_percentage,
)
    for i in 1:max_epochs
        @info "Epoch $i"
        train, test, validation = step!(
            model, ds_train, ds_test, ds_validation;
            loss_function=loss, optimizer, measures,
        )
        @show train
        @show test
        @show validation
        push!(
            recorded_measures,
            merge(prefix_labels.((train, test, validation), (:train_, :test_, :validation_))...),
        )

        let
            global fig, ax = subplots()
            losses = ["train_loss", "test_loss", "validation_loss"]
            plts = foreach(1:3, losses) do i, loss
                ax.plot(
                    getproperty(recorded_measures, loss),
                    label=replace(loss, '_' => ' '),
                )
            end
            ax.legend()
            display(fig)
        end

        if i >= min_epochs
            if argmin(recorded_measures[!, :test_loss]) <= i - early_stopping_n
                @info "Loss has not decreased in the last 10 epochs, stopping training"
                break
            elseif test.loss / train.loss > 1 + early_stopping_percentage / 100
                @info "Test loss more than $early_stopping_percentage% greater than training loss, stopping training"
                break
            end
        end
    end
end


text:
using FilteredGroupbyMacro
using DataFrames
using Test

df = DataFrame(
    a = repeat(1:10, inner = 10),
    b = repeat('a':'e', inner = 20),
    c = 1:100,
)

@testset "Base functionality" begin
    dt_1 = @by df[!, :b, d = :a .* :c]
    by_1 = by(df, :b, d = (:a, :c) => x -> x.a .* x.c)

    # do both versions do the same thing?
    @test dt_1 == by_1

    # test chaining of multiple groupings
    dt_2 = @by dt_1[!, :b, e = sum(:d .^ 2)]
    @test dt_2 == @by df[!, :b, d = :a .* :c][
        !, :b, e = sum(:d .^ 2)]

    # test filtering
    dt_3 = @by df[:a .< 8, :b, d = sum(:c)]
    by_3 = by(df[df.a .< 8, :], :b, d = (:c,) => x -> sum(x.c))
    @test dt_3 == by_3
end

@testset "In Function" begin
    function test(df)
        @by df[!, :a, csum = sum(:c)]
    end

    dt = test(df)

    function test2(df)
        localsumfunc = x -> sum(x)
        @by df[!, :a, csum = localsumfunc(:c)]
    end

    dt2 = test2(df)
end

@testset "Sorted keyword" begin
    df2 = DataFrame(a = 2:-1:1, b = 1:2)
    dt = @by df2[!, :a, c = 1 * :b]
    @test dt.a == [2, 1]
    @test dt.c == [1, 2]

    dt2 = @by df2[!, :a, c = 1 * :b; sort = true]
    @test dt2.a == [1, 2]
    @test dt2.c == [2, 1]
end

@testset "Assign syntax" begin
    df2 = DataFrame(a = [1, 1, 2, 2], b = 1:4)
    dt = @by df2[!, :a, c := sum(:b)]
    @test dt == DataFrame(a = [1, 1, 2, 2], b = 1:4, c = [3, 3, 7, 7])
end

# @testset "Single symbol" begin
#     df2 = DataFrame(a = 2:-1:1, b = 1:2)
#     dt = @by df2[!, :a, c = :b] # just :b
#     @test dt.a == [2, 1]
#     @test dt.c == [1, 2]
# end

max_stars_repo_name: jkrumbiegel/GroupbyIndexingMacro


text:
using Genie, Genie.Router

route("/hello") do
    "Hello World."
end

up()


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/>.

"""
    nrbextrude(curve::NURBS1D{I,F},vector::Vector{F}
              )::NURBS2D{I,F} where {I<:Integer,F<:AbstractFloat}

Construct a NURBS surface by extruding a NURBS curve.

Constructs a NURBS surface by extruding a NURBS curve along a defined vector.
The NURBS curve forms the U direction of the surface edge, and extruded along
the vector in the V direction. Note NURBS surfaces cannot be extruded.

# Arguments:
- `crv`: NURBS curve to extrude, see nrbmak.
- `vec`: Vector along which the curve is extruded.

# Output:
- `srf`: NURBS surface constructed.

# Examples:
```julia
julia> srf = nrbextrude(crv,vec)
```

Form a hollow cylinder by extruding a circle along the z-axis.

```julia
srf = nrbextrude(nrbcirc(),[0.0;0.0;1.0]);
```
"""
function nrbextrude(curve::NURBS1D{I,F},vector::Vector{F}
                   )::NURBS2D{I,F} where {I<:Integer,F<:AbstractFloat}
if length(vector) != 3
  throw(ArgumentError("Extrusion vector has to be a 3D vector!"));
end # if
coefs = cat(curve.coefs,vectrans(vector)*curve.coefs,dims=3);
return nrbmak(coefs,[curve.knots,vec([0.0 0.0 1.0 1.0])]);
end # nrbextrude


text:
##
## The flaming sheep, contributed by Geordan Rosario (geordan@csua.berkeley.edu)
##
"""
    function flaming_sheep()

The flaming sheep, contributed by Geordan Rosario (geordan@csua.berkeley.edu)

# Example

```jldoctest
julia> cowsay("Baaahhh!", cow=Cowsay.flaming_sheep)
 __________
< Baaahhh! >
 ----------
  \\            .    .     .
   \\      .  . .     `  ,
    \\    .; .  : .' :  :  : .
     \\   i..`: i` i.i.,i  i .
      \\   `,--.|i |i|ii|ii|i:
           UooU\\.'@@@@@@`.||'
           \\__/(@@@@@@@@@@)'
                (@@@@@@@@)
                `YY~~~~YY'
                 ||    ||
```
"""
function flaming_sheep(;eyes="oo", tongue="  ", thoughts="\\")
the_cow =
"""
  $thoughts            .    .     .
   $thoughts      .  . .     `  ,
    $thoughts    .; .  : .' :  :  : .
     $thoughts   i..`: i` i.i.,i  i .
      $thoughts   `,--.|i |i|ii|ii|i:
           U$(eyes)U\\.'@@@@@@`.||'
           \\__/(@@@@@@@@@@)'
                (@@@@@@@@)
                `YY~~~~YY'
                 ||    ||
"""
return the_cow
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    for T in (Bool, Int, Float32, Float64)
        for A in (Vector, Matrix)
            precompile(stridedpointer, (A{T},))
            precompile(stridedpointer, (LinearAlgebra.Adjoint{T,A{T}},))
        end
    end

    precompile(offset_ptr, (Symbol, Symbol, Char, Int, Int, Int, Int, Int, Bool, Int))
    precompile(vload_quote_llvmcall_core, (Symbol, Symbol, Symbol, Int, Int, Int, Int, Bool, Bool, Int))
    precompile(vstore_quote, (Symbol, Symbol, Symbol, Int, Int, Int, Int, Bool, Bool, Bool, Bool, Int))

    precompile(reset_features!, ())
    precompile(safe_topology_load!, ())
    precompile(redefine_attr_count, ())
    precompile(redefine_cache, (Int,))

    # precompile(_pick_vector_width, (Type, Vararg{Type,100}))
    # the `"NATIVE_PRECOMPILE_VECTORIZATIONBASE" ∈ keys(ENV)` isn't respected, seems
    # like it gets precompiled anyway given that the first condition is `true`.
    # if VERSION ≥ v"1.7.0-DEV.346" && "NATIVE_PRECOMPILE_VECTORIZATIONBASE" ∈ keys(ENV)
    #     set_features!()
    #     for T ∈ (Float32, Float64)
    #         W = pick_vector_width(T)
    #         precompile(>=, (Int, MM{W, 1, Int}))
    #         for op ∈ (-, Base.FastMath.sub_fast)
    #             precompile(op, (Vec{W, T}, ))
    #         end
    #         for op ∈ (+, -, *, Base.FastMath.add_fast, Base.FastMath.sub_fast, Base.FastMath.mul_fast)
    #             precompile(op, (Vec{W, T}, Vec{W, T}))
    #         end
    #         for op ∈ (VectorizationBase.vfmadd, VectorizationBase.vfmadd_fast)
    #             precompile(op, (Vec{W, T}, Vec{W, T}, Vec{W, T}))
    #         end
    #     end
    # end
end


text:
function show(io::IO, oc::Core.OpaqueClosure)
    A, R = typeof(oc).parameters
    show_tuple_as_call(io, Symbol(""), A; hasfirst=false)
    print(io, "::", R)
    print(io, "->◌")
end

function show(io::IO, ::MIME"text/plain", oc::Core.OpaqueClosure{A, R}) where {A, R}
    show(io, oc)
end

"""
    @opaque (args...) -> body

Marks a given closure as "opaque". Opaque closures capture the
world age of their creation (as opposed to their invocation).
This allows for more aggressive optimization of the capture
list, but trades off against the ability to inline opaque
closures at the call site, if their creation is not statically
visible.

!!! warning
    This interface is experimental and subject to change or removal without notice.
"""
macro opaque(ex)
    esc(Expr(:opaque_closure, ex))
end

max_stars_repo_name: pablosanjose/julia

max_issues_repo_name: pablosanjose/julia

max_forks_repo_name: pablosanjose/julia


text:
"""
    from_kwargs(convention!, ::Type{T}; kw...) where T

Convert keyword arguments to given option type `T` using `convention!`.
See also [`from_dict`](@ref).

# Convention

- `from_underscore_kwargs!`: use `_` to disambiguate subfields of the same name, this is the default behaviour.
- `from_field_kwargs!`: do not disambiguate subfields, errors if there are disambiguity
"""
function from_kwargs(convention!, ::Type{T}; kw...) where T
    d = OrderedDict{String, Any}()
    convention!(d, T; kw...)
    return from_dict(T, d)
end

"""
    from_kwargs(::Type{T}; kw...) where T

Convert keyword arguments to given option type `T` using the underscore convention.
"""
from_kwargs(::Type{T}; kw...) where T = from_underscore_kwargs(T; kw...)

"""
    from_underscore_kwargs(::Type{T}; kw...) where T

Convert keyword arguments to given option type `T` using the underscore convention.
"""
from_underscore_kwargs(::Type{T}; kw...) where T = from_kwargs(from_underscore_kwargs!, T; kw...)

"""
    from_field_kwargs(::Type{T}; kw...) where T

Convert keyword arguments to given option type `T` using the field keyword convention.
"""
from_field_kwargs(::Type{T}; kw...) where T = from_kwargs(from_field_kwargs!, T; kw...)

# NOTE: this is for compatibilty
"""
    from_kwargs!(d::AbstractDict{String}, ::Type{T}, prefix::Maybe{Symbol} = nothing; kw...) where T

Internal method for inserting keyword arguments to given dictionary object `d`. It will overwrite
existing keys in `d` if it is specified by keyword argument.
"""
function from_kwargs!(d::AbstractDict{String}, ::Type{T}, prefix::Maybe{Symbol} = nothing; kw...) where T
    return from_underscore_kwargs!(d, T, prefix; kw...)
end

function from_underscore_kwargs!(d::AbstractDict{String}, ::Type{T}, prefix::Maybe{Symbol} = nothing; kw...) where T
    validate_keywords(T, underscore_keywords(T); kw...)
    unsafe_from_underscore_kwargs!(d, T, prefix; kw...)
end

function unsafe_from_underscore_kwargs!(d::AbstractDict{String}, ::Type{T}, prefix::Maybe{Symbol} = nothing; kw...) where T
    return foreach_keywords!(d, T) do name, type
        key = underscore(prefix, name)

        from_kwargs_option_key!(d, type, name, key, kw) do field_d, field_type
            unsafe_from_underscore_kwargs!(field_d, field_type, key; kw...)
        end
    end
end

function from_field_kwargs!(d::AbstractDict{String}, ::Type{T}; kw...) where T
    validate_keywords(T, field_keywords(T); kw...)
    unsafe_from_field_kwargs!(d, T; kw...)
end

function unsafe_from_field_kwargs!(d::AbstractDict{String}, ::Type{T}; kw...) where T
    return foreach_keywords!(d, T) do name, type
        from_kwargs_option_key!(d, type, name, name, kw) do field_d, field_type
            unsafe_from_field_kwargs!(field_d, field_type; kw...)
        end
    end
end

function from_kwargs_option_key!(f, d::AbstractDict, ::Type{T}, name::Symbol, key::Symbol, kw) where T
    key_str = string(key)
    name_str = string(name)
    # shortcut
    if haskey(kw, key)
        d[name_str] = kw[key]
        return d
    end

    if is_option(T)
        field_d = OrderedDict{String, Any}()
        if haskey(d, name_str) && (d_value = d[name_str]) isa AbstractDict
            field_d = merge!(field_d, d_value)
        end

        f(field_d, T) # recurse into subfields
        if !isempty(field_d)
            d[name_str] = field_d
        end
    elseif T isa Union
        from_kwargs_option_key!(f, d, T.a, name, key, kw)
        from_kwargs_option_key!(f, d, T.b, name, key, kw)
    end
    return d
end

function validate_keywords(::Type{T}, keys = underscore_keywords(T); kw...) where T
    for (k, v) in kw
        k in keys || throw(InvalidKeyError(k, keys))
    end
    return
end

"""
    field_keywords(::Type{T}) where T

Return all the option type field names given `T`, error if there are duplicated sub-fields.
"""
function field_keywords(::Type{T}) where T
    return collect_field_keywords!(Symbol[], T, T)
end

"""
    underscore_keywords(::Type{T}) where T

Return keywords given `T` using the underscore convention.
"""
function underscore_keywords(::Type{T}) where T
    return collect_underscore_keywords!(Symbol[], T)
end

function foreach_keywords!(f, list, ::Type{T}) where T
    is_option(T) || return list

    for name in fieldnames(T)
        type = fieldtype(T, name)
        f(name, type)
    end
    return list
end

function underscore(prefix::Maybe{Symbol}, name)
    if prefix === nothing
        return name
    else
        return Symbol(prefix, :_, name)
    end
end

function collect_field_keywords!(list::Vector{Symbol}, ::Type{Top}, ::Type{T}) where {Top, T}
    return foreach_keywords!(list, T) do name, type
        if is_option(type)
            collect_field_keywords!(list, Top, type)
        elseif type isa Union
            name in list && error(msg)
            push!(list, name)
            # recurse into Union
            collect_field_keywords!(list, Top, type.a)
            collect_field_keywords!(list, Top, type.b)
        else
            name in list && throw(DuplicatedFieldError(name, Top))
            push!(list, name)
        end
    end
end

function collect_underscore_keywords!(list::Vector{Symbol}, ::Type{T}, prefix::Maybe{Symbol} = nothing) where T
    return foreach_keywords!(list, T) do name, type
        key = underscore(prefix, name)

        if is_option(type)
            collect_underscore_keywords!(list, type, key)
        elseif type isa Union
            push!(list, key)
            collect_underscore_keywords!(list, type.a, key)
            collect_underscore_keywords!(list, type.b, key)
        else
            push!(list, key)
        end
    end
end


text:
using Bolt
using ForwardDiff
using PyPlot
using BenchmarkTools

par = CosmoParams()
bg = Background(par)
ih = IonizationHistory(Peebles(), par, bg)
k_grid = quadratic_k(0.1bg.H₀, 1000bg.H₀, 100)
sf = source_grid(par, bg, ih, k_grid, BasicNewtonian())

ells = 100:50:1200
cl = cltt(ells, par, bg, ih, sf)

clf()
plt.plot(ells, cl .* ells.^2, "-")
ylabel(raw"$\ell^2 C_{\ell}^{TT}$")
xlabel(raw"$\ell$")
gcf()


text:
using Test

include("util.jl")

@testset "hypergraph_visualization" begin
    using HRGVisualization

    he1 = HyperEdge("S", ["1"], ["2"])
    he_dot = to_dot(he1)
    expected = """
    digraph HyperEdge {
        rankdir=LR
        labelloc=b
        color=white
        edge [arrowsize=0.5]
        s1[shape=point;label="";xlabel="1";fontsize=10]
        he[shape=box;label="S"]
        t1[shape=point;label="";xlabel="2";fontsize=10]
        {s1} -> he -> {t1}
    }
    """
    _test_normalized_strings_are_equal(he_dot,expected)
    _print_dot(he_dot)

    he2 = HyperEdge("X", ["1", "3"], ["2", "4"])
    he_dot = to_dot(he2)
    expected = """
    digraph HyperEdge {
        rankdir=LR
        labelloc=b
        color=white
        edge [arrowsize=0.5]
        s1[shape=point;label="";xlabel="1";fontsize=10]
        s2[shape=point;label="";xlabel="3";fontsize=10]
        he[shape=box;label="X"]
        t1[shape=point;label="";xlabel="2";fontsize=10]
        t2[shape=point;label="";xlabel="4";fontsize=10]
        {s1,s2} -> he -> {t1,t2}
    }
    """
    _test_normalized_strings_are_equal(he_dot,expected)
    _print_dot(he_dot)

    he3 = HyperEdge("Y", ["2", "4"], ["5"])
    he4 = HyperEdge("Z", ["5"], ["6", "7", "8"])
    hg = HyperGraph{String}(HyperEdge[he1, he2, he3, he4])

    hg_dot = to_dot(hg)
    expected = """
    digraph HyperGraph {
        rankdir=LR
        labelloc=b
        color=white
        edge [arrowsize=0.5]
        n1[shape=point;label="";xlabel="1";fontsize=10]
        n2[shape=point;label="";xlabel="2";fontsize=10]
        n3[shape=point;label="";xlabel="3";fontsize=10]
        n4[shape=point;label="";xlabel="4";fontsize=10]
        n5[shape=point;label="";xlabel="5";fontsize=10]
        n6[shape=point;label="";xlabel="6";fontsize=10]
        n7[shape=point;label="";xlabel="7";fontsize=10]
        n8[shape=point;label="";xlabel="8";fontsize=10]
        he1[shape=box;label="S"]
        he2[shape=box;label="X"]
        he3[shape=box;label="Y"]
        he4[shape=box;label="Z"]
        {n1} -> he1 -> {n2}
        {n1,n3} -> he2 -> {n2,n4}
        {n2,n4} -> he3 -> {n5}
        {n5} -> he4 -> {n6,n7,n8}
    }
    """
    _test_normalized_strings_are_equal(hg_dot,expected)
    _print_dot(hg_dot)
end


text:
#
# Object headers
#

const HM_NIL = 0x00
const HM_DATASPACE = 0x01
const HM_LINK_INFO = 0x02
const HM_DATATYPE = 0x03
const HM_FILL_VALUE_OLD = 0x04
const HM_FILL_VALUE = 0x05
const HM_LINK_MESSAGE = 0x06
const HM_EXTERNAL_FILE_LIST = 0x07
const HM_DATA_LAYOUT = 0x08
const HM_BOGUS = 0x09
const HM_GROUP_INFO = 0x0a
const HM_FILTER_PIPELINE = 0x0b
const HM_ATTRIBUTE = 0x0c
const HM_OBJECT_COMMENT = 0x0d
const HM_SHARED_MESSAGE_TABLE = 0x0e
const HM_OBJECT_HEADER_CONTINUATION = 0x0f
const HM_SYMBOL_TABLE = 0x10
const HM_MODIFICATION_TIME = 0x11
const HM_BTREE_K_VALUES = 0x12
const HM_DRIVER_INFO = 0x13
const HM_ATTRIBUTE_INFO = 0x14
const HM_REFERENCE_COUNT = 0x15

immutable ObjectStart
    signature::UInt32
    version::UInt8
    flags::UInt8
end
ObjectStart(flags::UInt8) = ObjectStart(OBJECT_HEADER_SIGNATURE, 2, flags)
define_packed(ObjectStart)

# Reads the start of an object including the signature, version, flags,
# and (payload) size. Returns the size.
function read_obj_start(io::IO)
    os = read(io, ObjectStart)
    os.signature == OBJECT_HEADER_SIGNATURE || throw(InvalidDataException())
    os.version == 2 || throw(UnsupportedVersionException())

    if (os.flags & OH_TIMES_STORED) != 0
        # Skip access, modification, change and birth times
        skip(io, 128)
    end
    if (os.flags & OH_ATTRIBUTE_PHASE_CHANGE_VALUES_STORED) != 0
        # Skip maximum # of attributes fields
        skip(io, 32)
    end

    read_size(io, os.flags)
end

immutable HeaderMessage
    msg_type::UInt8
    size::UInt16
    flags::UInt8
end
define_packed(HeaderMessage)


text:
# module for optimizing policy
# by value iteration

using FreqTables, NamedArrays

include(joinpath(split(@__FILE__, "src")[1], "src/model/goal/goal.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/model/agent/puddle_ignore_agent.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/common/state_transition/state_transition.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/model/puddle/puddle.jl"))

mutable struct DynamicProgramming
  pose_min
  pose_max
  widths
  goal
  index_nums
  indexes
  value_function
  is_final_state
  policy
  actions
  state_transition_probs
  depths
  delta_time
  puddle_coef

  # init
  function DynamicProgramming(widths, goal,
                              puddles,
                              delta_time,
                              sampling_num;
                              puddle_coef=100,
                              lower_left=[-4.0, -4.0],
                              upper_right=[4.0, 4.0])
    self = new()
    self.pose_min = [lower_left[1], lower_left[2], 0.0]
    self.pose_max = [upper_right[1], upper_right[2], 2*pi]
    self.widths = widths
    self.goal = goal

    self.index_nums = round.(Int64, (self.pose_max - self.pose_min)./self.widths)
    nx, ny, nt = self.index_nums[1], self.index_nums[2], self.index_nums[3]
    self.indexes = vec(collect(Base.product(0:nx-1, 0:ny-1, 0:nt-1)))
    self.value_function, self.is_final_state = init_value_function(self)
    self.policy = init_policy(self)

    actions_set = Set([Tuple(self.policy[i[1]+1, i[2]+1, i[3]+1, :]) for i in self.indexes])
    self.actions = [a for a in actions_set]
    self.state_transition_probs = init_state_transition_probs(self, delta_time, sampling_num)

    self.depths = depth_means(self, puddles, sampling_num)

    self.delta_time = delta_time
    self.puddle_coef = puddle_coef

    return self
  end
end

function depth_means(self::DynamicProgramming, puddles, sampling_num)
  # sampling
  dx = range(0, self.widths[1], length=sampling_num)
  dy = range(0, self.widths[2], length=sampling_num)
  samples = vec(collect(Base.product(dx, dy)))

  # mean of depth at each cell
  tmp = zeros(Tuple([self.index_nums[1], self.index_nums[2]]))
  for xy in collect(Base.product(0:self.index_nums[1]-1, 0:self.index_nums[2]-1))
    for s in samples
      center = self.pose_min + self.widths.*[xy[1], xy[2], 0] + [s[1], s[2], 0]
      for p in puddles
        tmp[xy[1]+1, xy[2]+1] += p.depth * inside(p, center)
      end
    end
    tmp[xy[1]+1, xy[2]+1] /= sampling_num^2
  end

  return tmp
end

function init_state_transition_probs(self::DynamicProgramming, delta_time, sampling_num)
  # sampling
  dx = range(0.001, self.widths[1]*0.999, length=sampling_num)
  dy = range(0.001, self.widths[2]*0.999, length=sampling_num)
  dt = range(0.001, self.widths[3]*0.999, length=sampling_num)
  samples = vec(collect(Base.product(dx, dy, dt)))

  # move sampled points and record delta index
  tmp = Dict()
  for a in self.actions
    for i_t in 0:self.index_nums[3]-1
      transition = []
      for s in samples
        # before transition
        before_pose = [s[1], s[2], s[3] + i_t*self.widths[3]] + self.pose_min
        before_index = [1, 1, i_t+1]

        # after transition
        after_pose = state_transition(a[1], a[2], delta_time, before_pose)
        after_index = get_index(self, after_pose)

        # index difference
        push!(transition, after_index-before_index)
      end
      # frequency of each index difference
      freq_tbl = freqtable(transition)
      freqs = [value for (name, value) in enumerate(freq_tbl)]
      diffs = [diff for diff in names(freq_tbl, 1)]

      # probability
      probs = freqs ./ sampling_num^3
      tmp[a, i_t+1] = [[d, p] for (d, p) in zip(diffs, probs)]
    end
  end

  return tmp
end

function get_index(self::DynamicProgramming, pose)
  return Int64.(floor.((pose - self.pose_min)./self.widths))
end

function init_policy(self::DynamicProgramming)
  tmp = zeros(Tuple([self.index_nums[1], self.index_nums[2], self.index_nums[3], 2]))

  for index in self.indexes
    center = self.pose_min + self.widths.*(index .+ 0.5)
    tmp[index[1]+1, index[2]+1, index[3]+1, :] .= policy(center, self.goal)
  end

  return tmp
end

function init_value_function(self::DynamicProgramming)
  v = Array{Float64}(undef, Tuple(self.index_nums))
  f = zeros(Tuple(self.index_nums))

  for index in self.indexes
    f[index[1]+1, index[2]+1, index[3]+1] = final_state(self, index)
    if f[index[1]+1, index[2]+1, index[3]+1] == true
      v[index[1]+1, index[2]+1, index[3]+1] = self.goal.value
    else
      v[index[1]+1, index[2]+1, index[3]+1] = -100.0
    end
  end

  return v, f
end

function final_state(self::DynamicProgramming, index)
  lower_left = self.pose_min + self.widths.*index
  upper_right = self.pose_min + self.widths.*(index.+1)

  corners = [[lower_left[1], lower_left[2]],
             [lower_left[1], upper_right[2]],
             [upper_right[1], lower_left[2]],
             [upper_right[1], upper_right[2]]]

  # check arrived at goal
  # all of corner point should be inside of goal area
  result = true
  for c in corners
    if inside(self.goal, c) == false
      result = false
    end
  end

  return result
end

# give penalty when robot moved out of space
# and then, return state inside space
function out_correction(self::DynamicProgramming, index)
  out_reward = 0.0

  # direction
  index[3] = (index[3] + self.index_nums[3])%self.index_nums[3]

  for i in 1:2
    if index[i] < 0
      index[i] = 0
      out_reward = -1e100
    elseif index[i] >= self.index_nums[i]
      index[i] = self.index_nums[i] - 1
      out_reward = -1e100
    end
  end

  return index, out_reward
end

function action_value(self::DynamicProgramming, action, index; out_penalty=true)
  value = 0.0

  for delta_prob in self.state_transition_probs[(action, index[3]+1)]
    after, out_reward = out_correction(self, [index[1]+1, index[2]+1, index[3]+1]+delta_prob[1])
    reward = -self.delta_time*self.depths[after[1]+1, after[2]+1]*self.puddle_coef - self.delta_time + out_reward*out_penalty
    value += (self.value_function[after[1]+1, after[2]+1, after[3]+1] + reward) * delta_prob[2]
  end

  return value
end

function policy_evaluation_sweep(self::DynamicProgramming)
  max_delta = 0.0

  for i in self.indexes
    if self.is_final_state[i[1]+1, i[2]+1, i[3]+1] != true
      q = action_value(self, Tuple(self.policy[i[1]+1, i[2]+1, i[3]+1, :]), i)

      delta = abs(self.value_function[i[1]+1, i[2]+1, i[3]+1] - q)

      if delta > max_delta
        max_delta = delta
      else
        max_delta = max_delta
      end

      self.value_function[i[1]+1, i[2]+1, i[3]+1] = q
    end
  end

  return max_delta
end

function value_iteration_sweep(self::DynamicProgramming)
  max_delta = 0.0

  for i in self.indexes
    if self.is_final_state[i[1]+1, i[2]+1, i[3]+1] != true
      max_q = -1e100
      max_a = nothing

      qs = [action_value(self, a, i) for a in self.actions]
      max_q = findmax(qs)[1]
      max_a = self.actions[argmax(qs)]

      delta = abs(self.value_function[i[1]+1, i[2]+1, i[3]+1] - max_q)
      if delta > max_delta
        max_delta = delta
      else
        max_delta = max_delta
      end

      self.value_function[i[1]+1, i[2]+1, i[3]+1] = max_q
      self.policy[i[1]+1, i[2]+1, i[3]+1, :] .= max_a
    end
  end

  return max_delta
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
module TestMessages

using ..Fixtures
using Log4jl
using Log4jl.Messages
import Log4jl: format, formatted, parameters
using FactCheck

msgInput = ["test", 10, 11.]

facts("Massages") do
    context("should have required methods implemented") do
        msg = Fixtures.IncompleteMessage()
        @fact_throws AssertionError format(msg)
        @fact_throws AssertionError parameters(msg)
        @fact_throws AssertionError formatted(msg)
    end
    context("can be created from string input") do
        msg = SimpleMessage(msgInput[1])
        @fact format(msg) --> msgInput[1]
        @fact parameters(msg) --> nothing
        @fact formatted(msg) --> msgInput[1]
        @fact parameters(SimpleMessage(msg)) --> Any[msg]
    end
    context("can be created from object input") do
        msg = ObjectMessage(msgInput)
        @fact format(msg) --> string(msgInput)
        @fact parameters(msg) --> Any[msgInput]
        @fact formatted(msg) --> string(msgInput)
    end
    context("can be created from parameterized string input") do
        msgPattern = "Test: {}, {}, {}."
        msg = ParameterizedMessage(msgPattern, msgInput...)
        @fact format(msg) --> msgPattern
        @fact parameters(msg) --> msgInput
        @fact formatted(msg) --> "Test: test, 10, 11.0."
        @fact parameters(ParameterizedMessage(msg)) --> Any[msg]
    end
    context("can be created from 'printf' formatted input") do
        msgPattern = "Test: %s, %d, %.1f"
        msg = PrintfFormattedMessage(msgPattern, msgInput...)
        @fact format(msg) --> msgPattern
        @fact parameters(msg) --> msgInput
        @fact formatted(msg) --> "Test: test, 10, 11.0"
        @fact parameters(PrintfFormattedMessage(msg)) --> Any[msg]
    end
end

end


text:
using MelGeneralizedCepstrums

import MelGeneralizedCepstrums: _mcep, _mgcep
import SPTK

function perf(f1, f2, params...; name="foo", n=100, N=1024)
    println("benchmark: $name")

    srand(98765)
    x = rand(N)

    @show params
    f1(x, params...)
    f2(x, params...)

    @time begin
        elapsed = @elapsed for i=1:n
            f1(x, params...)
        end
    end

    @time begin
        elapsed_sptk = @elapsed for i=1:n
            f2(x, params...)
        end
    end

    r = elapsed/elapsed_sptk
    println("$r x slower than SPTK implementation")
end

gc_enable(false)

gc()
perf(_mgcep, SPTK.mgcep, 25, 0.41, -0.1; name="mgcep", n=200, N=1024)

gc()
perf(_mcep, SPTK.mcep, 25, 0.41; name="mcep", n=2000, N=1024)

gc()
perf(mc2b, SPTK.mc2b, 0.41; name="mb2b", n=20000, N=25)

gc()
perf(b2mc, SPTK.b2mc, 0.41; name="b2mc", n=20000, N=25)

gc()
perf(freqt, SPTK.freqt, 30, 0.41; name="freqt", n=20000, N=25)


text:
using PyCall
# using ITensorPySCF
pyscf = pyimport("pyscf")
fci = pyimport("pyscf.fci")
using ITensors

#
# Create Simple Molecule
#

# mol = pyscf.gto.M(atom = "N 0 0 0; N 0 0 1;", basis = "sto3g", verbose = 3)
# mol = pyscf.gto.M(atom = "o 0 0 0; o 0 0 1;", basis = "sto3g", verbose = 3)
mol = pyscf.gto.M(atom = """C      0.00000    0.00000    0.00000
  H      0.00000    0.00000    1.08900
  H      1.02672    0.00000   -0.36300
  H     -0.51336   -0.88916   -0.36300
  H     -0.51336    0.88916   -0.36300""", basis = "sto3g", verbose = 3)

# Run HF
mf = pyscf.scf.RHF(mol).run()
println("RHF Energy (Ha): ", mf.e_tot)

# Create shorthands for 1- and 2-body integrals in MO basis
mo = mf.mo_coeff
n = size(mo, 1)
one_body = mo' * mf.get_hcore() * mo
two_body = reshape(mol.ao2mo(mf.mo_coeff, aosym = 1), n, n, n, n)

# FCI (i.e. exact diagonalization)
cisolver = fci.FCI(mf)
cisolver.kernel()
println("FCI Energy (Ha): ", cisolver.e_tot)

#
# Setup for MPS Calculation
#
t = one_body
V = 0.5 * permutedims(two_body, (3, 2, 1, 4))
n_occ = mf.mo_occ
e_nuclear = mf.energy_nuc()

# Create operators and MPO
hamiltonian = OpSum()
hamiltonian += e_nuclear, "Id", 1
for i in 1:n, j in 1:n
  hamiltonian .+= t[i, j], "Cdagup", i, "Cup", j
  hamiltonian .+= t[i, j], "Cdagdn", i, "Cdn", j
end
for i in 1:n, j in 1:n, k in 1:n, l in 1:n
  if norm(V[i, j, k, l]) > 1e-15
    if (i ≠ j) && (k ≠ l) # Otherwise the terms are exactly zero
      hamiltonian .+= V[i, j, k, l], "Cdagup", i, "Cdagup", j, "Cup", k, "Cup", l
      hamiltonian .+= V[i, j, k, l], "Cdagdn", i, "Cdagdn", j, "Cdn", k, "Cdn", l
    end
    hamiltonian .+= V[i, j, k, l], "Cdagup", i, "Cdagdn", j, "Cdn", k, "Cup", l
    hamiltonian .+= V[i, j, k, l], "Cdagdn", i, "Cdagup", j, "Cup", k, "Cdn", l
  end
end

s = siteinds("Electron", n; conserve_qns = true)
H = MPO(hamiltonian, s)

#
# DMRG
#

# Check that we can make an MPS properly
occupation_to_state(n) = n == 0 ? 1 : (n == 2 ? 4 : error("Occupation is $n"))
ψmf = MPS(s, occupation_to_state.(n_occ))
e_mf_mps = inner(ψmf', H, ψmf)
println("Energy Error from MF MPS (Ha) ", abs(e_mf_mps - mf.e_tot))

# Initialize our MPS
ψ0 = randomMPS(s, occupation_to_state.(n_occ); linkdims = 40)
@show inner(ψ0', H, ψ0)

# DMRG Settings
sweeps = Sweeps(10)
setmaxdim!(sweeps, 100, 300)
setcutoff!(sweeps, 1e-7)
setnoise!(sweeps, 1e-6, 1e-7, 1e-8, 0.0)

# Run DMRG
e, ψ = dmrg(H, ψ0, sweeps)
println("DMRG Error ", abs(e - cisolver.e_tot))


text:
function Base.sleep(sec::Real, token::CancellationToken)
    # Create a cancel source with a timeout
    timer_src = CancellationTokenSource(sec)

    timer_token = get_token(timer_src)

    # Create a cancel source that cancels either if the timeout source cancels,
    # or when the passed token cancels
    combined = CancellationTokenSource(timer_token, token)

    # Wait for the combined source to cancel
    wait(get_token(combined))

    if is_cancellation_requested(timer_src)
        return
    else
        throw(OperationCanceledException(token))
    end
end


text:
module Backends # Octo

struct UnsupportedError <: Exception
    msg
end

function backend(adapter::Module) # UnsupportedError
    sym = nameof(adapter)
    L = Symbol(sym, :Loader)
    mod = Main
    if isdefined(mod, L)
        getfield(mod, L)
    else
        try
            path = normpath(@__DIR__, "Backends", string(sym, ".jl"))
            Base.include(mod, path)
        catch err
            throw(UnsupportedError(string("error on ", adapter)))
        end
    end
end

end # module Octo.Backends


text:
@testset "SimilaritySearch for ManifoldLearning api" begin
    X, L = ManifoldLearning.scurve(segments=5)
    for t in [ExactEuclidean, ApproxEuclidean, ExactManhattan, ApproxChebyshev]  # some examples
        M = fit(Isomap, X, nntype=t, maxoutdim=2)
        Y = predict(M)
        n = size(X, 2)
        @test size(Y) == (2, n)
    end
end