text:
using DataFrames
using CSV
using PyPlot
using Random
using StatsBase
using Statistics

using Distributions
using LabelledArrays
using LinearAlgebra

include("../get_home_dir.jl")

pygui(true)

df_raw = CSV.read(string(PROJECT_HOME,"/data/moflux_land_data_skipyear_hourly2.csv"), DataFrame);
df_raw[!,"RAIN"] *= 2;

avg_et_daily = mean(df_raw.LE/44200)*18.02/1000 * 60*60*24;
avg_rain_daily = mean(df_raw.RAIN)*24;

include(string(PROJECT_HOME,"/simulation_code/sim_vary_new_stomata_med_varyLayer.jl"));
include("time_averaging.jl")
include("mironov.jl")
include("tau_omega_funs.jl")

N = 24*365*1
istart = 24*365*2 + 1;
soil0 = 0.375;

deltaT = FT(60*60);
#alpha = FT(50);
#alpha = FT(20);

#alpha = FT(5);
#nsoil = FT(1.2);

#alpha = FT(50)#FT(50);
#nsoil = FT(1.6);

alpha = FT(10)#FT(50);
nsoil = FT(1.5);

parnames = ["Kplant","Zsoil","Capacitance","Ksoil",
"Terrain_slope","root_dist_factor",
"Vcmax","P63","g1"];
par_indices = [3,4,6,8,9,10,
				1,5,7 ];

function run_sim_3layer(vcmax_par::FT, k_frac::FT, k_plant::FT,
      z_soil::FT, weibB::FT, vol_factor::FT, g1::FT, k_soil::FT,
      smc_runoff::FT, exp_root_dist::FT,
      canopy_pvslope::FT, trunk_pvslope::FT)
	return run_sim_varyB(vcmax_par, k_frac, weibB, FT(4), k_plant,
    k_soil, z_soil, istart, N, soil0, vol_factor, 1e-5, 3, df_raw, g1, deltaT, alpha, nsoil,0,0,
    smc_runoff, exp_root_dist,canopy_pvslope, trunk_pvslope);
end

c1 = 1

pars0 = convert(Array{FT}, [40, 0.5, 12, 2000, 8, 1.0, 250,0.4e-6,
         0.35,2,1/20,1/20])
sim_res1 = run_sim_3layer(pars0...);
cm1 = mean(sim_res1[6][:,1:3],dims=2);
pdLWP = Float64.(replace(sim_res1[1].LWP_predawn, missing => NaN));

tsoil = sim_res1[1].T_SOIL .+ 273.15;
tair = sim_res1[1].T_AIR .+ 273.15;
tcan = (sim_res1[1].LW_OUT / 0.97 / K_STEFAN ) .^ 0.25;
laiM = sim_res1[1].LAI_modis;

vodA = 0.067;
vodB = 0.82;
vodC = 0.051;

tb1 = get_TB_2(sim_res1[2][:,1], tsoil, tcan, cm1, laiM, vodA, vodB, vodC);

mynormfac = 1/36;
synthetic_ll = zeros(1000)
for ri in 1:1000
	noisyH = tb1[1] + 1.3*randn(length(tb1[1]));
	noisyV = tb1[2] + 1.3*randn(length(tb1[2]));
	obsLL_H = logpdf(MvNormal(tb1[1],1.3),noisyH);
	obsLL_V = logpdf(MvNormal(tb1[2],1.3),noisyV);
	obsLL_total = mynormfac * (obsLL_H + obsLL_V);
	synthetic_ll[ri] = obsLL_total;
end

prior_min = [10, 0.01, 0.1, 500, 0.75,  0.01,100,1e-6,0.01,0.01];
prior_max = [120,0.9,  50, 3000, 10, 10,1000,    5e-4,50.0,5];
prior_dist = Product(Uniform.(log.(prior_min), log.(prior_max)));
constLL = logpdf(prior_dist,log.(pars0[1:10]));

#=
chaindata = CSV.read("assim_code/sample_likelihood/post_par.csv",DataFrame);
chaindata2 = CSV.read("assim_code/sample_likelihood/allobs_c2/post_par.csv",DataFrame);
chaindata3 = CSV.read("assim_code/sample_likelihood/allobs_c3/post_par.csv",DataFrame);

chaindata = CSV.read("assim_code/sample_likelihood/only6/post_par.csv",DataFrame);
chaindata2 = CSV.read("assim_code/sample_likelihood/only6_c2/post_par.csv",DataFrame);
chaindata3 = CSV.read("assim_code/sample_likelihood/only6_c3/post_par.csv",DataFrame);
=#

chaindata1 = CSV.read("assim_code/sample_like_2/all1/post_par.csv",DataFrame);
chaindata2 = CSV.read("assim_code/sample_like_2/all2/post_par.csv",DataFrame);
chaindata3 = CSV.read("assim_code/sample_like_2/all3/post_par.csv",DataFrame);

post_LL = vcat([x[2000:25:end,21] for x in [chaindata1,chaindata2,chaindata3]]...);
post_errvar = vcat([x[2000:25:end,22] for x in [chaindata1,chaindata2,chaindata3]]...);

#=
figure()
hist(synthetic_ll .+ constLL,density=1,label="Theoretical")
hist(post_LL,density=1,alpha=0.5,label="From MCMC")
xlabel("Log likelihood")
ylabel("Density")
legend()
=#

sse_mcmc = (-(post_LL .- constLL) - 243*log.(sqrt.(post_errvar)) .- 243/2*log(2*pi)) .* (2*post_errvar);
sse_theory = (-synthetic_ll .- 243*log.(1.3) .- 243/2*log(2*pi)) .* (2*1.3^2);
sse_mcmc2 = (-(post_LL .- constLL) .- 243*log.(1.3) .- 243/2*log(2*pi)) .* (2*1.3^2);

post_tab = vcat([x[2000:50:end,:] for x in [chaindata1,chaindata2,chaindata3]]...);
#=
figure()
parI = 7
hist(post_tab[:,parI])
xlim(log(prior_min[parI]),log(prior_max[parI]))
=#
prior_mean = (Array(chaindata1[1,1:10]) + Array(chaindata2[1,1:10]) + Array(chaindata3[1,1:10]))/3;
post_mean = mean(Array(post_tab[:,1:10]),dims=1);

post_samp = Array(post_tab[end,1:10]);

sim_res_prior = run_sim_3layer(convert(Array{FT},exp.(prior_mean))...,FT(1/20),FT(1/20));
sim_res_post = run_sim_3layer(convert(Array{FT},exp.(post_mean))...,FT(1/20),FT(1/20));
sim_res_samp = run_sim_3layer(convert(Array{FT},exp.(post_samp))...,FT(1/20),FT(1/20));

cm_prior = mean(sim_res_prior[6][:,1:3],dims=2);
cm_post = mean(sim_res_post[6][:,1:3],dims=2);
cm_samp = mean(sim_res_samp[6][:,1:3],dims=2);

N = 24*365*12
istart = 24*365*0 + 1;

#=
println("long 1")
sim_res1long = run_sim_3layer(pars0...);
cm1_long = mean(sim_res1long[6][:,1:3],dims=2);
#println("long 2")
#sim_res_post_long = run_sim_3layer(convert(Array{FT},exp.(post_mean))...,FT(1/20),FT(1/20));
println("long 3")
par_samp = convert(Array{FT},exp.(post_samp));
sim_res_samp_long = run_sim_3layer(convert(Array{FT},exp.(post_samp))...,FT(1/20),FT(1/20));

cm2_long = mean(sim_res_samp_long[ 6][:,1:3],dims=2);
=#
#parametrize as stomatal p63 and quantity abs(xylem p63)-abs(stomatal p63), which must be more than 0
psi_list = collect(-8:0.01:0);
wS = exp.(-(psi_list/2) .^ 2);
wX = exp.(-(psi_list/2.5) .^ 2);
#amount of water leaving via ET in this time step is proportional to wS*(psoil - psi_list)
#dV = k dP
#amount of water incoming is proportional to wX
#change in potential is inversely proportional to wX
#runaway will occur when there is a positive feedback between ET and dP
#i.e. if the change in
function plot_tree(x,y)
	leaf_levels = [12,16,18]
	trunk_level = 9
	root_levels = collect(-8:1:-1)
	figure()
	for i in 1:3
		plot([x[i],x[i+3],x[7]],[leaf_levels[i]+1,leaf_levels[i],trunk_level],
			 "ko-")
	end
	plot([1,1]*x[7],[0,trunk_level],"ko-")
	for i in 1:8
		plot([y[i],x[7+i],x[7]],[root_levels[i],root_levels[i],0],
			 "ko-")
	end
end

#=
figure()
plot(sim_res1[6][:,1:3],color="green",alpha=0.5)
plot(pdLWP,"ko")
plot(sim_res1[6][:,4:6],color="orange",alpha=0.5)
plot(sim_res1[6][:,7],color="blue",alpha=0.5)
plot(sim_res1[6][:,8:15],color="brown",alpha=0.5)

daylist = collect(1:length(sim_res1long[1].VPD))/24;
figure()
plot(daylist, sim_res1long[6][:,1:3],color="green",alpha=0.5)
#plot(pdLWP2,"ko")
plot(daylist, sim_res1long[6][:,4:6],color="orange",alpha=0.5)
plot(daylist, sim_res1long[6][:,7],color="blue",alpha=0.5)
plot(daylist,sim_res1long[6][:,8:15],color="brown",alpha=0.5)

figure()
plot(daylist[1:(24*5):end],get_daily(cm1_long,24*5))
xlabel("Time (days since 1 Jan 2005)")
ylabel("5-day average LWP (MPa)")

ylabel("LWP (MPa)")

=#
#=
pars1 = [3.1927648046462975,-2.6467226953772807,2.099458400464035,7.6636431524168485,2.200646613730597,-1.1787635223907147,5.74270506922163,-14.480768091978758,-1.1286712435746784,0.673368997239951,3.3204521717948343,-3.764499807981834,1.9399104710601602,7.329693071448318,2.3019307796709465,1.2814998731299516,5.5803917979572315,-14.278778878664784,-1.0954914081021205,0.7479192727655555,-826.554645047097,1.6900000000000002,0.05087957867351999,0.802651566774747,0.04858293737916987];
pars1b = convert(Array{FT}, exp.(pars1[1:10]))

post1 = convert(Array,CSV.read("assim_code/sample_output_mar9/postLeaf.csv",DataFrame));
plot(daylist,post1,"r",alpha=0.33)
plot(daylist,cm1_long,"k")
xlim(365,365*2)
ylabel("LWP (MPa)")
xlabel("Time (days since 1 Jan 2005)")
=#
figure()
subplot(2,2,1)
plot(sim_res1[1].ETmod[(244*24):(246*24)])
plot(sim_res1[1].ETmod[(253*24):(255*24)])
ylim(-0.001,0.008)
#xlabel("Time (hours)")
title("Model ET (mol/m2/s)")

subplot(2,2,2)
plot(sim_res1[1].LE[(244*24):(246*24)]/42400)
plot(sim_res1[1].LE[(253*24):(255*24)]/42400)
ylim(-0.001,0.008)
#xlabel("Time (hours)")
title("Flux tower ET (mol/m2/s)")

subplot(2,2,3)
plot(mean(sim_res1[10][(244*24):(246*24),:],dims=2))
plot(mean(sim_res1[10][(253*24):(255*24),:],dims=2))
#xlabel("Time (hours)")
title("Model stomatal\nconductance (mol/m2/s)")

subplot(2,2,4)
plot(cm1[(244*24):(246*24)])
plot(cm1[(253*24):(255*24)])
#xlabel("Time (hours)")
title("Model LWP (MPa)")

figure()
subplot(2,2,1)
plot(mean(sim_res1[8][(244*24):(246*24),:],dims=2))
plot(mean(sim_res1[8][(253*24):(255*24),:],dims=2))
#xlabel("Time (hours)")
title("APAR (umol/m2/s)")

subplot(2,2,2)
plot(sim_res1[1].RelHum[(244*24):(246*24)])
plot(sim_res1[1].RelHum[(253*24):(255*24)])
title("Relative humidity")

subplot(2,2,3)
plot(sim_res1[1].T_AIR[(244*24):(246*24)])
plot(sim_res1[1].T_AIR[(253*24):(255*24)])
title("Temperature (deg C)")

subplot(2,2,4)
plot(sim_res1[1].VPD[(244*24):(246*24)])
plot(sim_res1[1].VPD[(253*24):(255*24)])
#xlabel("Time (hours)")
title("VPD (Pa)")

sval = collect(FT,0.1:0.001:0.55);
scurve = VanGenuchten{FT}(stype = "Ozark",
                                           α = FT(10),
                                           n = FT(1.5),
                                          Θs = 0.55,
                                          Θr = 0.067);
spotval = [soil_p_25_swc(scurve,x) for x in sval];

figure()
plot(sim_res1[1].SMC,pdLWP,"o")
plot(sval,spotval)


text:
export BCQDLearner

"""
    BCQDLearner(;kwargs)

See paper: [Benchmarking Batch Deep Reinforcement Learning Algorithms](https://arxiv.org/abs/1910.01708).

# Keyword arguments

- `approximator`::[`ActorCritic`](@ref): used to get Q-values (Critic) and logits (Actor) of a state.
- `target_approximator`::[`ActorCritic`](@ref): similar to `approximator`, but used to estimate the target.
- `γ::Float32 = 0.99f0`, reward discount rate.
- `τ::Float32 = 0.005f0`, the speed at which the target network is updated.
- `θ::Float32 = 0.99f0`, regularization coefficient.
- `threshold::Float32 = 0.3f0`, determine whether the action can be used to calculate the Q value.
- `batch_size::Int=32`
- `update_freq::Int`: the frequency of updating the `approximator`.
- `update_step::Int = 0`
- `rng = Random.GLOBAL_RNG`
"""
mutable struct BCQDLearner{
    Aq<:ActorCritic,
    At<:ActorCritic,
    R<:AbstractRNG,
} <: AbstractLearner
    approximator::Aq
    target_approximator::At
    γ::Float32
    τ::Float32
    θ::Float32
    threshold::Float32
    batch_size::Int
    update_freq::Int
    update_step::Int
    rng::R
    # for logging
    actor_loss::Float32
    critic_loss::Float32
end

function BCQDLearner(;
    approximator::Aq,
    target_approximator::At,
    γ::Float32 = 0.99f0,
    τ::Float32 = 0.005f0,
    θ::Float32 = 1f-2,
    threshold::Float32 = 0.3f0,
    batch_size::Int = 32,
    update_freq::Int = 10,
    update_step::Int = 0,
    rng = Random.GLOBAL_RNG,
) where {Aq<:ActorCritic, At<:ActorCritic}
    copyto!(approximator, target_approximator)
    BCQDLearner(
        approximator,
        target_approximator,
        γ,
        τ,
        θ,
        threshold,
        batch_size,
        update_freq,
        update_step,
        rng,
        0.0f0,
        0.0f0,
    )
end

Flux.functor(x::BCQDLearner) = (Q = x.approximator, Qₜ = x.target_approximator),
y -> begin
    x = @set x.approximator = y.Q
    x = @set x.target_approximator = y.Qₜ
    x
end

function (learner::BCQDLearner)(env)
    s = state(env)
    s = Flux.unsqueeze(s, ndims(s) + 1)
    s = send_to_device(device(learner), s)
    q = learner.approximator.critic(s)
    prob = softmax(learner.approximator.actor(s), dims=1)
    mask = Float32.((prob ./ maximum(prob, dims=1)) .> learner.threshold)
    new_q = q .* mask .+ (1.0f0 .- mask) .* -1f8
    new_q |> vec |> send_to_host
end

function RLBase.update!(learner::BCQDLearner, batch::NamedTuple)
    AC = learner.approximator
    target_AC = learner.target_approximator
    γ, τ, θ = learner.γ, learner.τ, learner.θ
    batch_size = learner.batch_size
    D = device(AC)

    s, a, r, t, s′ = (send_to_device(D, batch[x]) for x in SARTS)
    a = CartesianIndex.(a, 1:batch_size)

    prob = softmax(AC.actor(s′))
    mask = Float32.((prob ./ maximum(prob, dims=1)) .> learner.threshold)
    q′ = AC.critic(s′)
    a′ = argmax(q′ .* mask .+ (1.0f0 .- mask) .* -1f8, dims=1)
    target_q = target_AC.critic(s′)

    target = r .+ γ .* (1 .- t) .* vec(target_q[a′])

    ps = Flux.params(AC)
    gs = gradient(ps) do
        # Critic loss
        q_t = AC.critic(s)
        critic_loss = Flux.Losses.huber_loss(q_t[a], target)

        # Actor loss
        logit = AC.actor(s)
        log_prob = -log.(softmax(logit, dims=1))
        actor_loss = mean(log_prob[a])

        ignore() do
            learner.actor_loss = actor_loss
            learner.critic_loss = critic_loss
        end

        actor_loss + critic_loss + θ * mean(logit .^ 2)
    end

    update!(AC, gs)

    # polyak averaging
    for (dest, src) in zip(
        Flux.params([learner.target_approximator]),
        Flux.params([learner.approximator]),
    )
        dest .= (1 - τ) .* dest .+ τ .* src
    end
end

ext: jl


text:
using Random
using JSON
using CPLEX
using JuMP
include("../src/JOH.jl")
include("../experiments/SSIT.jl")
include("../MDMKP/MDMKP.jl")

all_problems = MDMKP.load_folder()

for problem in all_problems
	method = SE.make_SSIT_methods(2)[1]
	i = 1

	m = MDMKP.create_MIPS_model(problem, with_optimizer(CPLEX.Optimizer,
		CPX_PARAM_SCRIND=0))
	JOH.Matheur.set_threads!(m, method.num_threads)
	JOH.Matheur.set_tolerance!(m, method.tolerances[i])
	time_limit = method.times[i]
	JOH.Matheur.set_time!(m, time_limit)

	JOH.Matheur.silent_optimize_slow!(m)
	optimize!(m)

	row = JOH.Matheur.get_DF_row(m, elapsed_time=0, index=i, getdettime=CPLEX.CPXgetdettime)
	# println(row)

	out = read(`top -bn1 -p $(getpid())`, String)
	println(split(split(out,  "\n")[end-1])[6])
	println()
end

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
matrix = rand(10,5)
fmv = FastBEAST.FullMatrixView(matrix, Vector(1:5), Vector(1:10), 10, 5)

@test matrix == fmv.matrix

v = rand(5)

@test matrix*v  == fmv*v

adjointfmv = adjoint(fmv)
adjoint_adjointfmv = adjoint(adjointfmv)

@test adjoint_adjointfmv == fmv

v2 = rand(10)
@test adjoint(matrix)*v2 == adjoint(fmv)*v2

rightmatrix = rand(5,10)
leftmatrix = rand(3,5)

lmv = FastBEAST.LowRankMatrixView(rightmatrix, leftmatrix, Vector(1:10), Vector(1:3), 3, 10)

@test lmv.rightmatrix == rightmatrix
@test lmv.leftmatrix == leftmatrix

@test lmv*v2  ==  leftmatrix*(rightmatrix*v2)

v3 = rand(3)

@test adjoint(lmv)*v3  ==  adjoint(rightmatrix)*(adjoint(leftmatrix)*v3)

A = rand(100,100)
B = rand(100,100)
using LinearAlgebra
@test opnorm(A) ≈ estimate_norm(A) atol=1e-3
@test opnorm(A-B)/opnorm(B) ≈ estimate_reldifference(A,B, tol=1e-6) atol=1e-3

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: JoshuaTetzner/FastBEAST

max_forks_repo_name: JoshuaTetzner/FastBEAST


text:
#EDAData: .sNp (Touchstone) file reader
#-------------------------------------------------------------------------------
#=Notes:
   -Since .sNp files don't store # of lines (read until EOF), the algorithm uses
    push!() to read in data files.  TODO: Change algorithm if push!(::Vector)
    implementation is too slow.
=#

#==Constants
===============================================================================#
const SNP_FSCALE_MAP = Dict("GHZ"=>1e9, "MHZ"=>1e6, "KHZ"=>1e3, "HZ"=>1)

#==Main data structures
===============================================================================#

mutable struct SNPReader
	r::DSVReader
	numports::Int
end

#==Helper functions
===============================================================================#
iscomment(line::AbstractString) = length(line) < 1 || '!' == line[1]

#Read mangitude/angle data
snp_read_ma(r::SNPReader) = read(r.r, DataFloat) * e^(im*read(r.r, DataFloat))

#Read dB/angle data
snp_read_db(r::SNPReader) = 10^(read(r.r, DataFloat)/20) * e^(im*read(r.r, DataFloat))

#Read real/imaginary data
snp_read_ri(r::SNPReader) = complex(read(r.r, DataFloat), read(r.r, DataFloat))

const SNP_READ_MAP = Dict("MA"=>snp_read_ma, "DB"=>snp_read_db, "RI"=>snp_read_ri)

#==Open/read/close functions
===============================================================================#
#TODO: Should "open" read all header info???
function Base.open(::Type{SNPReader}, path::String; numports=0)
	#TODO: autodetect numports from path
	if numports < 1
		error("Must specify number of ports > 0")
	end
	#Open with specific reader to ensure consistent behaviour:
	r = open(DSVReader, path)
	return SNPReader(r, numports)
end

Base.close(r::SNPReader) = close(r.r)

function Base.read(::Type{SNPReader}, path::String; numports=0)
	#Default open ... do definitions do not support kwargs.
	reader = open(SNPReader, path, numports=numports)
	try
		return readall(reader)
	finally
		close(reader)
	end
end

#==Main reader algorithm
===============================================================================#
#Returns DataF1 instead of DataFreq... in case step size is not constant.
function readall(r::SNPReader)
	str = String #WANTCONST
	x = DataFloat[]
	y = Array{Vector{Complex{DataFloat}}}(undef, r.numports,r.numports)

	for row = 1:r.numports, col = 1:r.numports
		y[row, col] = Complex{DataFloat}[]
	end

	line = ""
	while iscomment(line)
		line = lstrip(readline(r.r))
	end

	if !(length(line) < 1 || '#' == line[1])
		error("Options line not found where expected.")
	end

	line = strip(line[2:end])

	local xunit, nptype, datafmt, refres
	optreader = DSVReader(IOBuffer(line))
	try
		xunit = uppercase(read(optreader, str))
		nptype = uppercase(read(optreader, str))
		datafmt = uppercase(read(optreader, str))
		refres = uppercase(read(optreader, str))
		refres = read(optreader, DataFloat)
	catch e
		rethrow(e)
	finally
		close(optreader)
	end

	snp_read = SNP_READ_MAP[datafmt]

	#Move to data start:
	line = ""
	datastart = 0
	while iscomment(line)
		datastart = position(r.r.s)
		line = lstrip(readline(r.r))
	end
	seek(r.r.s, datastart)

	#Read all data:
	while !eof(r.r)
		#TODO: push! algorithm potentially slow...
		push!(x, read(r.r, DataFloat))

		for row = 1:r.numports, col = 1:r.numports
			push!(y[row, col], snp_read(r))
		end
	end

	x = x .* SNP_FSCALE_MAP[xunit]
	m = Array{DataF1}(undef, r.numports, r.numports)
	for row = 1:r.numports, col = 1:r.numports
		m[row, col] = DataF1(x, y[row, col])
	end

	istwoport = (2==r.numports)
	if istwoport
		#Special case: columns represent x11, x21, x12, x22... gross...
		x21 = m[1,2]
		m[1,2] = m[2,1]
		m[2,1] = x21
	end

	nwkwargs = Dict()
	nptype = Symbol(nptype)
	if :S == nptype
		push!(nwkwargs, :z0 => refres)
	end

	return Network(nptype, m; nwkwargs...)
end

#==High-level interface
===============================================================================#
read_snp(path::String, args...; kwargs...) = read(SNPReader, path, args...; kwargs...)

#Last line


text:
mutable struct DioOption
    f0_floor::Cdouble
    f0_ceil::Cdouble
    channels_in_octave::Cdouble
    frame_period::Cdouble  # msec
    speed::Cint  # (1, 2, ..., 12)
    allowed_range::Cdouble  # Threshold used for fixing the F0 contour.

    # allow uninitialized creation
    DioOption() = new()
end

function dio(x::Vector{Float64}, fs::Int; kwargs...)

    diooption = Ref(DioOption())
    @ccall libworld.InitializeDioOption(diooption::Ptr{DioOption})::Cvoid

    for (symbol, value) in kwargs
        setproperty!(diooption[], symbol, value)
    end
    frame_period = diooption[].frame_period

    n_dio_samples = @ccall libworld.GetSamplesForDIO(fs::Cint, length(x)::Cint, frame_period::Cdouble)::Cint

    # allocate f0 and temporal positions vectors
    f0 = Vector{Cdouble}(undef, n_dio_samples)
    temporal_positions = similar(f0)

    @ccall libworld.Dio(x::Ptr{Cdouble}, length(x)::Cint, fs::Cint, diooption::Ptr{DioOption},
        temporal_positions::Ptr{Cdouble}, f0::Ptr{Cdouble})::Cvoid

    f0, temporal_positions
end


text:
export eval_mcc_inf, df_fields, eval_runopt
export remove_branches!

#####################
##################### Evaluating inference of splits
#####################

#####################
##################### Evaluating MCC inference
#####################

function get_r(ρ::Float64, n::Int64, N::Int64, simtype::Symbol)
    if simtype == :kingman
        return ρ * n / N
    elseif simtype == :yule
        return ρ / N
    else
        @error "Unrecognized `simtype`."
    end
end

"""
    eval_naive_inf(N::Int64, n::Int64, ρ::Float64, simtype::Symbol;
        cutoff = 0.,
        Nrep = 1,
        sfields::Tuple = (:ρ,:cutoff),
        out = ""
    )
"""
function eval_naive_inf(N::Int64, n::Int64, ρ::Float64, simtype::Symbol;
    cutoff = 0.,
    Nrep = 1,
    sfields::Tuple = (:ρ,:cutoff),
    out = ""
)
    #
    args = Dict(:N=>N, :n=>n, :ρ=>ρ, :simtype=>simtype, :cutoff=>cutoff)
    #
    dat = DataFrame(df_fields())
    for f in sfields
        insertcols!(dat, f=>Any[])
    end
    #
    r = get_r(ρ, n, N, simtype)
    #
    function f(arg::ARGTools.ARG)
        t1, t2 = ARGTools.trees_from_ARG(arg)
        if cutoff != 0
            remove_branches!(t1, Distributions.Exponential(cutoff*N))
            remove_branches!(t2, Distributions.Exponential(cutoff*N))
        end
        init_splits = Dict(1=>SplitList(t1), 2=>SplitList(t2))
        trees = Dict(1=>t1, 2=>t2)
        MCCs, resolved_splits = computeMCCs!(trees; naive=true, output=:all)
        return MCCs, resolved_splits, init_splits
    end
    #
    for rep in 1:Nrep
        td = @timed _eval_mcc_inf(f, N, n, r, simtype=simtype)
        d = td[1]
        args[:time] = td[2]
        args[:bytes] = td[3]
        for f in sfields
            d[f] = args[f]
        end
        push!(dat, d)
    end
    #
    return dat
end

"""
    eval_real(N::Int, n::Int, ρ::Float64, simtype::Symbol;
        Nrep = 1,
        out = "",
    )
"""
function eval_real(N::Int, n::Int, ρ::Float64, simtype::Symbol;
    Nrep = 1,
    out = "",
)
    dat = DataFrame(df_fields())
    insertcols!(dat, :ρ=>Float64[])
    r = get_r(ρ, n, N, simtype)
    #
    function f(arg::ARGTools.ARG)
        t1, t2 = ARGTools.trees_from_ARG(arg)
        init_splits = Dict(1=>SplitList(t1), 2=>SplitList(t2))
        resolved_splits = Dict(1=>[], 2=>[])
        MCCs = Dict((1,2) => ARGTools.MCCs_from_arg(arg))
        return MCCs, resolved_splits, init_splits
    end
    #
    for rep in 1:Nrep
        td = @timed _eval_mcc_inf(f, N, n, r, simtype=simtype)
        d = td[1]
        d[:ρ] = ρ
        push!(dat, d)
    end
    #
    return dat
end

"""
    eval_runopt(γ::Real, N::Int64, n::Int64, ρ::Float64, simtype::Symbol;
    Tmin=0.001, dT=0.01, Tmax=1., Md=10, lk_sort=true,
    cutoff = 0.,
    Nrep = 1,
    preresolve = true,
    simulate_seqs = false,
    sfields::Tuple = (:ρ,:cutoff, :preresolve),
    out = "",
    verbose=false
)

Eval the performance of `SplitGraph.runopt` at inferring MCCs.
"""
function eval_runopt(γ::Real, N::Int64, n::Int64, ρ::Float64, simtype::Symbol;
    Tmin=0.001, dT=0.01, Tmax=1., Md=10, lk_sort=true,
    cutoff = 0.,
    Nrep = 1,
    preresolve = true,
    simulate_seqs = false,
    sfields::Tuple = (:ρ,:cutoff, :preresolve),
    out = "",
    verbose=false
)
    #
    args = Dict(:γ=>γ, :N=>N, :n=>n, :ρ=>ρ, :simtype=>simtype,
        :Md=>Md, :Tmin=>Tmin, :dT=>dT, :Tmax=>Tmax, :lk_sort=>lk_sort,
        :cutoff=>cutoff, :preresolve=>preresolve,
    )
    #
    dat = DataFrame(df_fields())
    for f in sfields
        insertcols!(dat, f=>Any[])
    end
    #
    r = get_r(ρ, n, N, simtype)
    Trange = reverse(Tmin:dT:Tmax)
    #
    function f(arg::ARGTools.ARG)
        let cutoff=cutoff, N=N
            t1, t2 = ARGTools.trees_from_ARG(arg)
            trees = Dict(1=>t1, 2=>t2)
            if cutoff != 0
                remove_branches!(t1, Distributions.Exponential(cutoff*N))
                remove_branches!(t2, Distributions.Exponential(cutoff*N))
            end
            init_splits = Dict(1=>SplitList(t1), 2=>SplitList(t2))
            oa = OptArgs(
                γ=γ,
                Tmin=Tmin, dT=dT, Tmax=Tmax,
                Md=Md,
                likelihood_sort=lk_sort,
                verbose=verbose,
            )
            try
                MCCs, resolved_splits = computeMCCs!(
                    trees, oa;
                    preresolve=preresolve, output=:all
                )
                return MCCs, resolved_splits, init_splits
            catch err
                mkpath("tmp")
                write_newick("tmp/tree1.nwk", trees[1])
                write_newick("tmp/tree2.nwk", trees[2])
                error(err)
            end
        end
    end
    #
    for rep in 1:Nrep
        td = @timed _eval_mcc_inf(f, N, n, r, simtype=simtype)
        d = td[1]
        args[:time] = td[2]
        args[:bytes] = td[3]
        for f in sfields
            d[f] = args[f]
        end
        push!(dat, d)
    end
    #
    return dat
end

function simulate_sequences!(trees::Dict, N, c = 0.5, L = 1000; polytomies=true)
    μ = c == 0 ? 1 / (N*L*0.2) * 1.25 : 1 / (N*L*c) * 1.25
    # Simulate sequences && prune branches with no mutations
    for (i,t) in trees
        TreeAlgs.Evolve.evolve!(t, L, μ;  seqkey = :selfseq)
        TreeTools.compute_mutations!(n -> n.data.dat[:selfseq], t, :realmuts)
        c != 0 && polytomies && TreeTools.delete_branches!(n->isempty(n.data.dat[:realmuts]), t)
    end
    #
    return nothing
end

function _eval_mcc_inf(rMCC, iMCC, t1::Tree)
    # Mean MCC sizes
    rmMCC = mean([length(x) for x in values(rMCC)])
    imMCC = mean([length(x) for x in values(iMCC)])
    # Common branches
    τc_p = 0.; τnc_p = 0.; τc_n = 0.; τnc_n = 0.
    νc_p = 0.; νnc_p = 0.; νc_n = 0.; νnc_n = 0.
    for n in Iterators.filter(x->!x.isroot, values(t1.lnodes))
        if TreeKnit.is_branch_in_mccs(n,iMCC) # Branch predicted to be shared with the other tree
            if TreeKnit.is_branch_in_mccs(n, rMCC) # Correct prediction
                τc_p += n.tau
                νc_p += 1.
            else
                τc_n += n.tau
                νc_n += 1.
            end
        else # Branch predicted to not be shared with the other tree
            if !TreeKnit.is_branch_in_mccs(n, rMCC) # Correct prediction
                τnc_p += n.tau
                νnc_p += 1.
            else
                τnc_n += n.tau
                νnc_n += 1.
            end
        end
    end
    νc_p /= (length(t1.lnodes)-1); νnc_p /= (length(t1.lnodes)-1); νc_n /= (length(t1.lnodes)-1); νnc_n /= (length(t1.lnodes)-1)
    T = sum(skipmissing(x.tau for x in values(t1.lnodes))) # Total branch length
    τc_p /= T; τnc_p /= T; τc_n /= T; τnc_n /= T

    return Dict(:τc_p => τc_p, :τnc_p => τnc_p, :τc_n => τc_n, :τnc_n => τnc_n,
                :νc_p => νc_p, :νnc_p => νnc_p, :νc_n => νc_n, :νnc_n => νnc_n,
                :mMCCsize=>imMCC, :nMCC=>length(iMCC))
end

function _eval_split_inf(true_splits, resolved_splits, init_splits)
    rs_p = 0
    rs_n = 0
    rs_init = length(init_splits) / length(true_splits)
    for rs in resolved_splits
        if in(rs, true_splits, usemask=false)
            rs_p += 1
        else
            rs_n += 1
        end
    end
    return Dict(:rs_p => rs_p / length(true_splits), :rs_n=> rs_n / length(true_splits),
                :rs_init => rs_init)

end

"""
    _eval_mcc_splits_inf(f, N::Int64, n0::Int64, r::Float64; v=true, simtype=:kingman)

Evaluate the inference of MCCs by function `f`. ARG simulation made with `(N,n0,r)`.
"""
function _eval_mcc_inf(f, N::Int64, n0::Int64, r::Float64; v=true, simtype=:kingman)

    arg = ARGTools.SimulateARG.simulate(N, r, n0, simtype=simtype)

    # Real MCCs
    rMCC = ARGTools.MCCs_from_arg(arg);
    trees = ARGTools.trees_from_ARG(arg)
    true_splits = Dict(i=>SplitList(t) for (i,t) in enumerate(trees))

    # Inferred MCCs
    iMCC, resolved_splits, init_splits = f(arg)
    t1 = trees[1]

    out = _eval_mcc_inf(rMCC, iMCC[1,2], t1)
    out_ = _eval_split_inf(true_splits[1], resolved_splits[1], init_splits[1])
    for (k,v) in out_
        out[k] = v
    end
    out[:N] = N
    out[:n0] = n0
    out[:r] = r
    return out
end

function eval_mcc_inf(f, N::Int64, n0::Int64, rrange;
        Nrep=25*ones(Int64, length(rrange)), v=true, simtype=:kingman)
    dat = DataFrame(df_fields())
    for (i,r) in enumerate(rrange)
        v && print("r=$(round(ρ,digits=5)), $i/$(length(rrange))                                                     \r")
        Nr = Nrep[i]
        for rep in 1:Nr
            push!(dat, _eval_mcc_inf(f, N, n0, r, v=v, simtype=simtype))
        end
    end
    return dat
end
function eval_mcc_inf(f, N::Int64, n0::Int64, r::Float64;
        Nrep=25, v=true, simtype=:kingman)
    dat = DataFrame(df_fields())
    for rep in 1:Nrep
        push!(dat, _eval_mcc_inf(f, N, n0, r, v=v, simtype=simtype))
    end
    return dat
end

"""
    df_fields()

Fields returned by `_eval_mcc_inf`:
- `N`, `n0`, `r`
- `ν(n)c_p(n)`: Fraction of inferred (non-)common branches, true (p) or false (n)
- `τ(n)c_p(n)`: Length of inferred (non-)common branches, true (p) or false (n) (normalized)
- `mMCCsize`: mean inferred MCC size
- `nMCC`: number of inferred MCCs
"""
function df_fields()
    return (N = Int64[], n0 = Int64[], r = Float64[], # Simulation parameters
        νc_p = Float64[], νnc_p = Float64[], # Fraction of correctly inferred common and non-common branches
        νc_n = Float64[], νnc_n = Float64[], # Fraction of incorrectly inferred common and non-common branches
        τc_p = Float64[], τnc_p = Float64[], # Total length of correctly inferred common and non-common branches
        τc_n = Float64[], τnc_n = Float64[], # Total length of incorrectly inferred common and non-common branches
        rs_p = Float64[], rs_n = Float64[], # Number of (in)-correctly resolved splits, scaled by total number of true splits
        rs_init = Float64[], # Initial number of splits scaled by total number of true splits
        mMCCsize = Float64[], # Mean MCC size (inferred)
        nMCC = Int64[], # Number of MCCs (inferred)
        )
end

#####################
##################### REMOVING BRANCHES
#####################
"""
    remove_branches!(n::TreeNode, p::Distribution)
"""
function remove_branches!(n::TreeNode, p::Distribution)
    if !ismissing(n.tau) && n.tau < rand(p)
        if !n.isleaf
            nr = delete_node!(n)
            for c in nr.child
                remove_branches!(c, p)
            end
        else
            n.tau = 0.
        end
    else
        for c in n.child
            remove_branches!(c, p)
        end
    end
end

"""
    remove_branches!(t::Tree, p)

Stochastically remove branches from `t` using probability distribution `p`
"""
function remove_branches!(t::Tree, p)
    remove_branches!(t.root, p)
    node2tree!(t, t.root)
    nothing
end
remove_branches!(t::Tree, c::Real, N::Int) = remove_branches!(t, Exponential(c*N))

max_stars_repo_name: PierreBarrat/OldTreeKnit

max_issues_repo_name: PierreBarrat/OldTreeKnit

max_forks_repo_name: PierreBarrat/OldTreeKnit


text:
module Day5
export get_inputs, get_solution1, get_solution2

using ..Advent2020.Utils: zero_based

## Input getting
function get_inputs()
    test_input1 = [
        "FBFBBFFRLR"
        "BFFFBBFRRR"
        "FFFBBBFRRR"
        "BBFFBBFRLL"
    ]
    test_input2 = nothing
    test_output1 = maximum([357, 70*8+7, 14*8+7, 102*8+4])
    test_output2 = nothing
    data = readlines(joinpath(@__DIR__, "input.txt"))
    return (; test_input1, test_input2, test_output1, test_output2, data)
end

## Solution functions
replace_n(str, pairs...) = String(replace(collect(str), pairs...))

binary_string(str; zero='0', one='1') = parse(Int, replace_n(str, zero=>'0', one=>'1'), base=2)

function parse_boarding_pass(pass)
    row = binary_string(pass[1:7], zero='F', one='B')
    col = binary_string(pass[8:10], zero='L', one='R')
    seat = row*8 + col
    return row, col, seat
end

# Part 1
function get_solution1(data)
    seats = last.(parse_boarding_pass.(data))
    return maximum(seats)
end

# Part 2
get_solution2(::Nothing) = nothing
function get_solution2(data)
    seats = sort(last.(parse_boarding_pass.(data)))
    seat_2 = false
    seat_1 = false
    for seat in minimum(seats)+1:maximum(seats)
        is_seat = seat in seats
        if seat_2 && !seat_1 && is_seat
            return seat - 1
        end
        seat_2 = seat_1
        seat_1 = is_seat
    end
    return nothing
end

end


text:
# perf or animation ?
const INFO = false
const record_animation = false

using Printf
using StaticArrays
#using Makie
using Statistics
using DelimitedFiles

# the type for the domain
struct Domain
    x_start::Float64
    x_end::Float64
end

# Initial solution
function InitSolution!(V::Array{Float64, 2}, domain::Domain)
    nx = size(V, 2)

    hx = (domain.x_end - domain.x_start) / nx
    xc = 0.5  * (domain.x_start + domain.x_end) + nx/8 * hx

    for i=1:nx
        x = domain.x_start + (i-0.5) * hx
        if abs(x - xc) < 0.2
            V[1, i] = 1
        else
            V[1, i] = 0
        end

        V[2, i] = 0
    end
end

# Model
function F!(f, V::Array{Float64, 2}, i::Int64)
    f[1] = V[2, i]
    if isapprox(V[1, i], 0., atol=eps(Float64), rtol=0)
        f[2] = 0
    else
        f[2] = V[2, i]^2 / V[1, i] + 0.5 * 9.81 * V[1, i]^2
    end
end

# Scheme
function LaxFriedrich!(V::Array{Float64, 2}, Vold::Array{Float64, 2}, lambdas::Array{Float64, 1}, dt, dx)
    nx = size(V, 2)
    Cx = dt / dx

    flux1 = @MVector [0., 0.]
    flux2 = @MVector [0., 0.]

    F!(flux1, Vold, 1)

    # left boundary
    V[2, 1] += Cx * (flux1[2] - lambdas[1] * Vold[2, 1])

    # interior edges
    for i=1:nx-1
        F!(flux2, Vold, i+1)
        lambda = max(lambdas[i], lambdas[i+1])

        for k=1:2
            flux = 0.5 * Cx * ((flux2[k] + flux1[k]) - lambda * (Vold[k, i+1] - Vold[k, i]))
            V[k, i] -= flux
            V[k, i+1] += flux
        end
        flux1, flux2 = flux2, flux1
    end

    # right boundary
    V[2, nx] -= Cx * (flux1[2] + lambdas[nx] * Vold[2, nx])
end

# Computes time step satisfying the CFL condition
function UpdateCFL!(lambdas::Array{Float64, 1}, V::Array{Float64, 2}, dx::Float64)
    nx = size(V, 2)

    M = 0
    for i=1:nx
        if isapprox(V[1, i], 0., atol=eps(Float64), rtol=0)
            lambdas[i] = 0.
        else
            lambdas[i] = abs(V[2, i] / V[1, i]) + sqrt(9.81 * V[1, i])
            M = max(M, lambdas[i])
        end
    end

    return 0.5 * dx / M
end

# initialization
const domain = Domain(0., 1.)
const Nx = 16384

# loop in time
const T = 2.0
const dx = (domain.x_end - domain.x_start) / Nx

let t=0
    V = Array{Float64, 2}(undef, 2, Nx)
    InitSolution!(V, domain)

    Vold = Array{Float64, 2}(undef, 2, Nx)
    lambdas = Array{Float64, 1}(undef, Nx)

    function update_to_time(tframe::Float64)
        while(t < tframe)
            dt = min(UpdateCFL!(lambdas, V, dx), tframe - t)
            @. Vold = V
            LaxFriedrich!(V, Vold, lambdas, dt, dx)
            t += dt
            INFO && @printf("\t%.5f / %.5f\n", t, tframe)
        end
    end

    if record_animation
        nframe_per_second = 1000
        nframe = convert(Int64, ceil(nframe_per_second * T)) + 1

        scene = Scene(resolution = (1920, 1080))
        x = range(domain.x_start, stop=domain.x_end, length=Nx)

        plt = lines!(scene, x, view(V, 1, :), color=:blue, linewidth=2)[end]
        record(scene, "./height1d.mp4", range(t, stop=T, length=nframe), framerate=30) do tframe
            println(tframe, " / ", T)
            update_to_time(tframe)
            plt[2] = view(V, 1, :)
        end
    else
        print("Warming up... ")
        update_to_time(eps(Float64))
        println("Done")

        t = 0
        InitSolution!(V, domain)

        println("Initial time t = ", t)
        @time update_to_time(T)
        println("End of simulation, t = ", t)

        println("mean(h) = ", mean(V[1, :]))
        writedlm("sol-cpu", V[1, :], "\n")
    end
end

max_stars_repo_name: aitzkora/BenchmarksPythonJuliaAndCo

max_issues_repo_name: aitzkora/BenchmarksPythonJuliaAndCo

max_forks_repo_name: aitzkora/BenchmarksPythonJuliaAndCo


text:
module _GaussLegendre

# http://www.johndcook.com/blog/2012/06/03/calculating-pi-with-agm-and-mpmath/
# http://numbers.computation.free.fr/Constants/Pi/iterativePi.html
# https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm

include("pi_iterative_gauss-legendre_int.jl")
include("pi_iterative_gauss-legendre_float.jl")

end


text:
module test_stdlib_random

using Test
using Random: Random, Sampler, MersenneTwister

const rand_nums = let
    if VERSION >= v"1.8.0-DEV.575"
        [5, 1]
    elseif VERSION >= v"1.6"
        [5, 8]
    else
        [8, 2]
    end
end

Random.seed!(0)
@test rand(Sampler(MersenneTwister, 1:10)) == rand_nums[1]
@test rand(Sampler(MersenneTwister, 1:10)) == rand_nums[2]

Random.seed!(0)
@test rand(Sampler(MersenneTwister, 1:10)) == rand_nums[1]
@test rand(Sampler(MersenneTwister, 1:10)) == rand_nums[2]

end # module test_stdlib_random


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function paczynski(u::AbstractFloat)
    absu = abs(u)
    if absu >= 1e5
        # You can verify that:
        #   u=1e5; (u*u + 2)/(u*sqrt(u*u + 4)) === 1.0
        return copysign(1.0, u)
    elseif absu <= 1e-8
        # You can verify that:
        #   u = 1e-8; (u*u + 2)/(u*sqrt(u*u + 4)) === inv(u)
        inv(u)
    else
        u2 = u*u
        # Using `hypot' in place of the square root would be an overkill since
        #   1e-8 < |u| < 1e5
        return (u2 + 2)/(u*sqrt(u2 + 4))
    end
end

"""
    paczynski(u) -> amplification

### Purpose ###

Calculate gravitational microlensing amplification of a point-like source by a
single point-like lens.

### Explanation ###

Return the [gravitational
microlensing](https://en.wikipedia.org/wiki/Gravitational_microlensing)
amplification of a point-like source by a single point-like lens, using
Paczyński formula

```math
A(u) = \\frac{u^2 + 2}{u\\sqrt{u^2 + 4}}
```

where \$u\$ is the projected distance between the lens and the source in units
of [Einstein radii](https://en.wikipedia.org/wiki/Einstein_radius).

In order to speed up calculations for extreme values of \$u\$, the following
asyntotic expressions for \$A(u)\$ are used:

```math
A(u) =
\\begin{cases}
 1/u & |u| \\ll 1 \\\\
 \\text{sgn}(u) & |u| \\gg 1
\\end{cases}
```

### Arguments ###

* `u`: projected distance between the lens and the source, in units of Einstein
  radii

### Output ###

The microlensing amplification for the given distance.

### Example ###

Calculate the microlensing amplification for \$u = 10^{-10}\$, \$10^{-1}\$,
\$1\$, \$10\$, \$10^{10}\$:

```jldoctest
julia> using AstroLib

julia> paczynski.([1e-10, 1e-1, 1, 10, 1e10])
5-element Array{Float64,1}:
  1.0e10
 10.037461005722337
  1.3416407864998738
  1.0001922892047386
  1.0
```

### Notes ###

The expression of \$A(u)\$ of microlensing amplification has been given by
Bohdan Paczyński in

* Paczynski, B. 1986, ApJ, 304, 1.
  DOI:[10.1086/164140](http://dx.doi.org/10.1086/164140),
  Bibcode:[1986ApJ...304....1P](http://adsabs.harvard.edu/abs/1986ApJ...304....1P)

The same expression was actually found by Albert Einstein half a century
earlier:

* Einstein, A. 1936, Science, 84, 506.
  DOI:[10.1126/science.84.2188.506](http://dx.doi.org/10.1126/science.84.2188.506),
  Bibcode:[1936Sci....84..506E](http://adsabs.harvard.edu/abs/1936Sci....84..506E)
"""
paczynski(u::Real) = paczynski(float(u))


text:
"""
    EGPpower(σ, ξ, κ)

*EGPpower* corresponds to the first extended GP model of Naveau et al. (2016), with the power law distribution G(v) = v^κ.

It is a three parameters family: κ controls the shape of the lower tail, σ is a scale parameter, and ξ controls the rate of upper tail decay.

```julia
EGPpower(σ, ξ, κ)   # EGP of Naveau et al. (2016) (type 1) with scale parameter σ, rate of upper tail decay ξ and shape of the lower tail κ.

params(d)           # Get the parameters, i.e. (σ, ξ, κ)
```

Reference :

* Naveau, P., Huser, R., Ribereau, P., and Hannart, A. (2016), Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection, Water Resour. Res., 52, 2753– 2769, doi:10.1002/2015WR018552.
"""
struct EGPpower{T<:Real} <: ContinuousUnivariateDistribution
    σ::T    # scale parameter
    ξ::T    # rate of upper tail decay
    κ::T    # shape of the lower tail

    function EGPpower{T}(σ::T, ξ::T, κ::T) where {T <: Real}
        new{T}(σ, ξ, κ)
    end
end

function EGPpower(σ::T, ξ::T, κ::T; check_args=true) where {T <: Real}
    check_args && @check_args(EGPpower, σ > zero(σ) && κ > zero(κ))
    return EGPpower{T}(σ, ξ, κ)
end

EGPpower(σ::Real, ξ::Real, κ::Real) = EGPpower(promote(σ, ξ, κ)...)
EGPpower(σ::Integer, ξ::Integer, κ::Integer) = EGPpower(float(σ), float(ξ), float(κ))

minimum(d::EGPpower) = 0.0
maximum(d::EGPpower) = Inf * (d.ξ >= 0) - ( d.σ / d.ξ ) * (d.ξ < 0)
insupport(d::EGPpower, x::Real) = minimum(d) <= x <= maximum(d)

#### Parameters

scale(d::EGPpower) = d.σ
decay(d::EGPpower) = d.ξ    # noms à vérifier...
shape(d::EGPpower) = d.κ

params(d::EGPpower) = (d.σ, d.ξ, d.κ)
partype(::EGPpower{T}) where {T} = T

#### Evaluation

function logpdf(d::EGPpower{T}, x::Real) where T<:Real
    μ = 0
    (σ, ξ, κ) = params(d)

    pd = GeneralizedPareto(μ, 1, ξ)

    lg(v::Real) = log(κ) + (κ-1)*log(v)

    p = -log(σ) + Distributions.logpdf(pd, x/σ) + lg(Distributions.cdf(pd, x/σ))

    return p
end

pdf(d::EGPpower, x::Real) = exp(logpdf(d, x))

function logcdf(d::EGPpower{T}, x::Real) where T<:Real
    # À vérifier
    μ = 0
    (σ, ξ, κ) = params(d)

    pd = GeneralizedPareto(μ, 1, ξ)

    #G(v::Real, κ::Real) = v^κ
    #lG(v::Real, κ::Real) = κ*log(v)

    p = κ*Distributions.logcdf(pd, x/σ)

    return p
end

cdf(d::EGPpower, x::Real) = exp(logcdf(d, x))

function quantile(d::EGPpower{T}, p::Real) where T<:Real
    @assert zero(p)<p<one(p) "the quantile level should be between 0 and 1."

    (σ, ξ, κ) = params(d)

    invG(p::Real, κ::Real) = p^(1/κ)

    x = (σ/ξ)*(((1 - invG(p,κ))^(-ξ)) - 1)

    return x
end

#### Sampling

function rand(rng::AbstractRNG, d::EGPpower)
    # Generate a Float64 random number uniformly in (0,1].
    u = 1 - rand(rng)

    invG(p::Real, κ::Real) = p^(1/κ)

    return (d.σ/d.ξ)*(((1 - invG(u,d.κ))^(-d.ξ)) - 1)
end

sampler(d::EGPpower) = d    # à vérifier


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.7.0-DEV
#     language: julia
#     name: julia-1.7
# ---

# %%
@show VERSION
@time using DifferentialEquations
@time using Plots
@time using ProgressMeter

# %% [markdown]
# Wave equaiton: $\displaystyle
# \frac{\partial^2 u}{\partial t^2} =
# \frac{\partial^2 u}{\partial x^2} +
# \frac{\partial^2 u}{\partial y^2}
# $

# %%
# the wave equation on the square with free boundary condition

#f!(du, v, u, p, t) = @. du = v

function g!(dv, v, u, p, t)
    h⁻² = p[1]
    @. @views dv[2:end-1,2:end-1] = h⁻²*(
        u[1:end-2,2:end-1] + u[3:end,2:end-1] + u[2:end-1,1:end-2] + u[2:end-1,3:end] - 4u[2:end-1,2:end-1])
    @. @views dv[1  ,2:end-1] = h⁻²*(u[2,    2:end-1] + u[1,  1:end-2] + u[1,  3:end] - 3u[1,  2:end-1])
    @. @views dv[end,2:end-1] = h⁻²*(u[end-1,2:end-1] + u[end,1:end-2] + u[end,3:end] - 3u[end,2:end-1])
    @. @views dv[2:end-1,1  ] = h⁻²*(u[2:end-1,2    ] + u[1:end-2,1  ] + u[3:end,1  ] - 3u[2:end-1,1  ])
    @. @views dv[2:end-1,end] = h⁻²*(u[2:end-1,end-1] + u[1:end-2,end] + u[3:end,end] - 3u[2:end-1,end])
    dv[1  ,1  ] = h⁻²*(u[1  ,2    ] + u[2    ,1  ] - 2u[1  ,1  ])
    dv[end,1  ] = h⁻²*(u[end,2    ] + u[end-1,1  ] - 2u[end,1  ])
    dv[1  ,end] = h⁻²*(u[1  ,end-1] + u[2    ,end] - 2u[1  ,end])
    dv[end,end] = h⁻²*(u[end,end-1] + u[end-1,end] - 2u[end,end])
end

x = y = range(-1, 1; length=201)
h = step(x)
p = (1/h^2,)

f(x, y) = exp(-(x^2+y^2)/2)
u0 = f.(x', y)
v0 = zero(u0)
tspan = (0.0, 10.0)

#@time prob = DynamicalODEProblem(g!, f!, v0, u0, tspan, p)
@time prob = SecondOrderODEProblem(g!, v0, u0, tspan, p)
sol = nothing; GC.gc()
@time sol = solve(prob; saveat=0.02)
sol = nothing; GC.gc()
@time sol = solve(prob; sageat=0.02)
sol = nothing; GC.gc()
@time sol = solve(prob; saveat=0.02);

# %%
x = y = range(-1, 1; length=201)
h = step(x)
p = (1/h^2,)

f(x, y) = exp(-(x^2+y^2)/2)
u0 = f.(x', y)
v0 = zero(u0)
tspan = (0.0, 20.0)

#@time prob = DynamicalODEProblem(g!, f!, v0, u0, tspan, p)
@time prob = SecondOrderODEProblem(g!, v0, u0, tspan, p)
sol = nothing; GC.gc()
@time sol = solve(prob; saveat=0.02);

# %%
umin = minimum(minimum.(u.x[2] for u in sol.u))
umax = maximum(maximum.(u.x[2] for u in sol.u))
umin, umax

# %%
L = 201
ts = range(tspan...; length=L) |> r -> [fill(r[1], 20); r; fill(r[end], 20)]
prog = Progress(length(ts), 0)
@gif for t in ts
    heatmap(x, y, sol(t).x[2]'; clim=(umin, umax), colorbar=false, size=(300, 316), c=:CMRmap)
    title!("t = $(round(t, digits=1))"; titlefontsize=12, aspectratio=1)
    next!(prog)
end

# %%
L = 201
ts = range(tspan...; length=L) |> r -> [fill(r[1], 20); r; fill(r[end], 20)]
prog = Progress(length(ts), 0)
@gif for t in ts
    surface(x, y, sol(t).x[2]'; zlim=(umin, umax), colorbar=false, c=:CMRmap, camera=(30, 60))
    title!("t = $(round(t, digits=1))"; titlefontsize=12)
    next!(prog)
end

# %%
L = 201
ts = range(tspan...; length=L) |> r -> [fill(r[1], 20); r; fill(r[end], 20)]
prog = Progress(length(ts), 0)
@gif for t in ts
    wireframe(x[begin:4:end], y[begin:4:end],
        sol(t).x[2][begin:4:end, begin:4:end]';
        zlim=(umin, umax), colorbar=false, c=:CMRmap, camera=(30, 60))
    title!("t = $(round(t, digits=1))"; titlefontsize=12)
    next!(prog)
end

# %%
g(x, y) = exp(-x^2/2-y^2/4)
u0 = g.(x', y)
v0 = zero(u0)
tspan = (0.0, 20.0)

@time prob = SecondOrderODEProblem(g!, v0, u0, tspan, p)
sol = nothing; GC.gc()
@time sol = solve(prob; saveat=0.02)

umax = maximum(maximum.(u.x[2] for u in sol.u))
umin = minimum(minimum.(u.x[2] for u in sol.u))

L = 201
ts = range(tspan...; length=L) |> r -> [fill(r[1], 20); r; fill(r[end], 20)]
prog = Progress(length(ts), 0)
@gif for t in ts
    surface(x, y, sol(t).x[2]'; zlim=(umin, umax), colorbar=false, c=:CMRmap, camera=(30, 60))
    title!("t = $(round(t, digits=1))"; titlefontsize=12)
    next!(prog)
end

# %%
L = 201
ts = range(tspan...; length=L) |> r -> [fill(r[1], 20); r; fill(r[end], 20)]
prog = Progress(length(ts), 0)
@gif for t in ts
    wireframe(x[begin:4:end], y[begin:4:end],
        sol(t).x[2][begin:4:end, begin:4:end]';
        zlim=(umin, umax), colorbar=false, c=:CMRmap, camera=(30, 60))
    title!("t = $(round(t, digits=1))"; titlefontsize=12)
    next!(prog)
end

# %%

max_stars_repo_path: 0003/Wave equation on the square 2021-06-09 Part 1.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

relaxed_infeas_tols = (tol_infeas = 1e-9,)
insts = OrderedDict()
insts["minimal"] = [
    ((0.7, 4, 1e-3, true, true),),
    ((0.7, 4, 1e-3, false, true), :SOCExpPSD),
    ((2.0, 4, 1e-3, true, false), nothing, relaxed_infeas_tols),
    ((2.0, 4, 1e-3, false, false), :SOCExpPSD, relaxed_infeas_tols),
    ]
insts["fast"] = [
    ((1.0, 2, 1e-3, true, false),),
    ((1.0, 2, 1e-3, false, false), :SOCExpPSD),
    ((4.0, 2, 1e-3, true, false),),
    ((4.0, 2, 1e-3, false, false), :SOCExpPSD),
    ]
insts["various"] = vcat(insts["minimal"], insts["fast"])
return (ContractionJuMP, insts)


text:
using PackageCompiler
function PackageCompiler.create_sysimage(
    packages;
    sysimage_path::String,
    project::String = dirname(Base.active_project()),
    precompile_execution_file::Union{String,Vector{String}} = String[],
    precompile_statements_file::Union{String,Vector{String}} = String[],
    incremental::Bool = true,
    filter_stdlibs::Bool = false,
    cpu_target::String = NATIVE_CPU_TARGET,
    script::Union{Nothing,String} = nothing,
    sysimage_build_args::Cmd = ``,
    include_transitive_dependencies::Bool = true,
    # Internal args
    base_sysimage::Union{Nothing,String} = nothing,
    julia_init_c_file = nothing,
    version = nothing,
    soname = nothing,
    compat_level::String = "major",
)
    write(sysimage_path, "test")
    return
end


text:
module JLFzf

import fzf_jll
using Pipe
import REPL
import REPL.LineEdit

export inter_fzf, read_repl_hist, insert_history_to_repl

"""
    read_repl_hist()

Find corresponding history file, read it, split by history records
and return resulted `Vector{String}` in reverse order.
"""
function read_repl_hist()
    rx = r"(^|\n)# time:[^\n]*\n"m
    # markers for history modes, used by `insert_history_to_repl()`
    rxj = r"^# mode: julia\n"m => ""
    rxh = r"^# mode: help\n"m => "?\n"
    rxp = r"^# mode: pkg\n"m => "]\n"
    rxs = r"^# mode: shell\n"m => ";\n"
    @pipe open(REPL.find_hist_file()) |>
          read |>
          String |>
          replace(_, r"\n$" => "") |>           # remove last new line
          replace(_, r"^\t"m => "") |>          # remove leading tabs
          replace(_, rxj) |>                    # replace mode with corresponding marker
          replace(_, rxh) |>                    #
          replace(_, rxp) |>                    #
          replace(_, rxs) |>                    #
          split(_, rx, keepempty = false) |>    # split by history records
          reverse |>
          unique                     # keep unique entries
end

"""
    inter_fzf(in_str::String, args...)

Run interactive fzf and return selected result, by default `in_str` should use `\n`s as the delimiter,
return selected string.
Additional arguments `args` for `fzf` are allowed.
"""
function inter_fzf(in_str::String, args...)
    fzf_jll.fzf() do exe
        if length(args) == 0
            return read(
                pipeline(Cmd(`$exe`, ignorestatus = true), stdin = IOBuffer(in_str)),
                String,
            ) |> chomp
        else
            return read(
                pipeline(
                    Cmd(`$exe $(args)`, ignorestatus = true),
                    stdin = IOBuffer(in_str),
                ),
                String,
            ) |> chomp
        end
    end
end

"""
    inter_fzf(ary::AbstractArray, args...)

Run interactive fzf with an array of inputs, return selected string.
Additional arguments `args` for `fzf` are allowed.
"""
function inter_fzf(ary::AbstractArray, args...)
    inter_fzf(join(ary, '\0'), args...)
end

function edit_insert_and_state_transition(mistate, line, mode)
    LineEdit.edit_insert(mistate, line)
    iobuffer = LineEdit.buffer(mistate)
    LineEdit.transition(mistate, Base.active_repl.interface.modes[mode]) do
        prompt_state = LineEdit.state(mistate, Base.active_repl.interface.modes[mode])
        prompt_state.input_buffer = copy(iobuffer)
    end
    nothing
end

"""
    insert_history_to_repl(mistate, line)

Paste found by `fzf` history `line` into REPL `mistate`.
"""
function insert_history_to_repl(mistate, line)
    rxh = r"^\?\n"
    rxp = r"^\]\n"
    rxs = r"^;\n"
    LineEdit.edit_clear(mistate)
    if occursin(rxh, line)
        edit_insert_and_state_transition(mistate, replace(line, rxh => ""), 3)
    elseif occursin(rxp, line)
        edit_insert_and_state_transition(mistate, replace(line, rxp => ""), 6)
    elseif occursin(rxs, line)
        edit_insert_and_state_transition(mistate, replace(line, rxs => ""), 2)
    else
        LineEdit.edit_insert(mistate, line)
    end
    LineEdit.refresh_line(mistate)
end

end


text:
module test
using ForwardDiff
#println("Loading the Module")

export TestFunc, Func1, Func3, Func2
export @ab, @u

export foo!, mySquare
foo!(x) = (x*=2) # multiply entries of x by 2 inplace

function TestFunc()
    println("In the function")
end

function Func1(x)
     println("This அம்ம function returns 4")
     return 4
end

function Func2(x)

    f(x::Vector) =  prod(log, x) * sum(sqrt, x);
    x = rand(1,10)
    g = x -> ForwardDiff.gradient(f, x); # g = ∇f

    x=100 .*rand(100)
    #return string(sin.(x))*"#"*string(g(x)) #"அம்மா"
    x = collect(10:.01:30)
    y = log.(x .+ sin.(x))
    z = (x .+ cos.(x))

    return string(y)*"#"*string(z) #"அம்மா"
end

function Func3(x)
     return (rand(1,10))
end

function mySquare(a)
       return sum(a)
end

macro ab(n)
    quote
    println($(esc(n)))
    "அம்மா"
    end
end

macro u(n, body)
  quote
    for i = 1:$(esc(n))
      println($(esc(body)))
    end
      "அம்மா"
  end
end

end

max_stars_repo_name: ResearchersHelp/Ayudha

max_issues_repo_name: ResearchersHelp/Ayudha

max_forks_repo_name: PiResearchTech/Ayudha


text:
using Documenter
using ExploringOceanData

makedocs(
    sitename = "ExploringOceanData",
    authors = "Andre",
    format = Documenter.HTML(collapselevel = 1, mathengine = MathJax3()),
    modules = [ExploringOceanData],
    pages = [
    "Home" => "index.md",
    ],
)

deploydocs(repo = "github.com/sandreza/ExploringOceanData.jl.git", devbranch = "main")


text:
using HyperbolicPlane, Plots

"""
`star_david(rad=1)` draw a Star of David in the hyperbolic plane.
"""
function star_david(r::Real = 1.0)
    tlist = [pi/2   pi/2+2pi/3   pi/2+4pi/3]
    points = [ HPoint(r,t) for t in tlist]
    T = HTriangle(points...)
    draw([T,T'])
    finish()
end


text:
module DirectStepIntegration

using Dierckx
using DifferentialEquations
using Plots

export tkris4, tkris5, razavi,rk_bi

function cvs(omgd,t)
	c=cos(omgd*t);
	s=sin(omgd*t);
	return c,s
end

function ltd_dqc(tz,p)
# Cho ra bieu dien luc tac dung dung quy cach, phu hop voi khoang thoi gian ma ket cau duoc phan tich. Vi du, luc tac dung luc dau chi duoc cho den 5 giay, nhung lai muon xet he o thoi diem 6 giay, khi do tu giay thu 5 tro di thi phai dat gia tri cua luc bang 0.
	h=tz[1];
	nstep=Int(tz[2]);
	hns=h*nstep;
	siz=size(p[2],1);
	luc=p[2];
	if size(p[1],1)==1
		ti=linspace(0,p[1]*(siz-1),siz);
	else
		ti=p[1];
	end
	tp=maximum(ti);

	if hns>tp
		noas=Int(ceil((hns-tp)/h));
		luc=[luc;zeros(noas,size(p[2],2))];
		ti=[ti;[hns-(i-1)*h for i=noas:-1:1]]
	else
		ind=searchsorted(ti,hns).stop;
		if hns>ti[ind]
			luc[ind+1,:]=luc[ind,:]+(luc[ind+1,:]-luc[ind,:])*(hns-ti[ind])/(ti[ind+1]-ti[ind]);
			ti=[ti[1:ind];hns];
			luc=luc[1:ind+1,:];
			ti=ti[1:ind+1]
		else
			luc=luc[1:ind,:];
			ti=ti[1:ind];
		end
	end
	return luc,ti
end

function tkris4(dth=([2 0;0 1],[0 0; 0 0],[96 -32;-32 32]),tz=(0.1,20),ic=([0;0],[0;0]),p=(0.1,ones(15,1)*[0 100]));
# Giai phuong trinh dao dong HE TUYEN TINH theo phuong phap Tkris
# dth: Dac trung cua he dth=(m,c,k) trong do m la ma tran khoi luong, c la ma tran he so can nhot va k la ma tran do cung
# tz: vecto co cac thong tin ve thoi gian, nhu sau tz=(h,nstep) trong do h la do lon buoc chia theo thoi gian va nstep la so buoc chia theo thoi gian.
# ic: Dieu kien ban dau ic=(cvbd,vtbd) trong do cvbd va vtbd tuong ung la vecto chuyen vi ban dau va vecto van toc ban dau
# p: Vecto tai trong p=([t_1; t_2; ...; t_n],[p(t_1);p(t_2);p(t_3);...;p(t_n)]) trong do p(t_i) la gia tri cua tai trong o thoi diem t_i. Gia tri t_1 bao gio cung bang 0. Khi cac thoi diem nay cach deu nhau, chi can vao 1 gia tri la delta_t o vecto dau tien (vi du bang 0.02).

	function xdhs_rt(ced,pse);
		hsc,hse,hsd=ced;
		ps,pe=pse;
		rn=-((8h1m2+2(cm+3mc)*h2+2(km+6mk+4c2)*h3/5+(2ck+kc)*h4/3+h5k2/7)*hsc+(4h1mc+(3mk+c2)*h2+(4ck+kc)*h3/5+h4k2/6)*hsd+(4h1mk+h2ck+h3k2/5)*hse-(h1m1+h2c1/5+h3k1/30)*ps-(3h1m1+4h2c1/5+h3k1/6)*pe);
		sn=-((6m2+2(cm+2mc)*h+(km+3mk+3c2)*h2/2+(3ck+2kc)*h3/5+h4k2/6)*hsc+(3mc+(2mk+c2)*h+(3ck+kc)*h2/4+h3k2/5)*hsd+(3mk+h1ck+h2k2/4)*hse-(m+h1c1/4+h2k1/20)*ps-(2m+3h1c1/4+h2k1/5)*pe);
		b=[rn;sn]
		aob=a\b;
		return aob[1:siz2], aob[siz2+1:2*siz2]
	end

	m,c,k=dth;
	h=tz[1];
	nstep=Int(tz[2]);
	hns=h*nstep;

	siz2=size(p[2],2);
	pt,t=ltd_dqc(tz,p);
	xxtt=[Spline1D(t,pt[:,i]) for i=1:siz2];

	u=zeros(nstep+1,siz2);
	ud=zeros(nstep+1,siz2);
	udd=zeros(nstep+1,siz2);

	u[1,:],ud[1,:]=ic;

	c2=c^2;
	c3=c^3;
	ck=c*k;
	cm=c*m;

	k2=k^2;
	kc=k*c;
	km=k*m;

	m2=m^2;
	mc=m*c;
	mk=m*k;

	h2=h^2;
	h3=h^3;
	h4=h^4;
	h5=h^5;
	h6=h^6;
	h7=h^7;

	h1c1=h*c;
	h1ck=h*ck;
	h1m1=h*m;
	h1mc=h*mc;
	h1mk=h*mk;
	h1m2=h*m2;
	h2c1=h2*c;
	h2ck=h2*ck;
	h2k1=h2*k;
	h2k2=h2*k2;
	h2m2=h2*m2;
	h3k1=h3*k;
	h3k2=h3*k2;
	h3m2=h3*m2;
	h4k2=h4*k2;
	h5k2=h5*k2;
	h6k2=h6*k2;
	h7k2=h7*k2;

	q11=144h3m2/5+8(cm+mc)*h4+4(3km+3mk+4c2)*h5/7+(ck+kc)*h6/2+h7k2/9;

	q12=18h2m2+12(2cm+3mc)*h3/5+(km+2mk+2c2)*h4+(4ck+3kc)*h5/7+h6k2/8;
	q21=q12';

	q22=12h1m2+9(cm+mc)*h2/2+3(2km+2mk+3c2)*h3/5+(ck+kc)*h4/2+h5k2/7;

	a=[q11 q12; q21 q22];

    for i=1:nstep
    	ps=[xxtt[j](h*(i-1)) for j=1:siz2];
    	pe=[xxtt[j](h*i) for j=1:siz2];

        hse=u[i,:];
        hsd=ud[i,:];
        hsc=(m\(ps-c*hsd-k*hse))/2;
        ced=(hsc,hse,hsd);
        pse=(ps,pe);
        hsa,hsb=xdhs_rt(ced,pse);
        u[i+1,:]=hsa*h4+hsb*h3+hsc*h2+hsd*h+hse;
        ud[i+1,:]=4hsa*h3+3hsb*h2+2hsc*h+hsd;
        udd[i+1,:]=m\(pe-c*ud[i+1,:]-k*u[i+1,:]);
    end
    return u,ud,udd
end

function tkris5(dth=([2 0;0 1],[0 0; 0 0],[96 -32;-32 32]),tz=(0.1,20),ic=([0;0],[0;0]),p=(0.1,ones(15,1)*[0 100]))
# Giai phuong trinh dao dong HE TUYEN TINH theo phuong phap Tkris
# dth: Dac trung cua he dth=(m,c,k) trong do m la ma tran khoi luong, c la ma tran he so can nhot va k la ma tran do cung
# tz: vecto co cac thong tin ve thoi gian, nhu sau tz=(h,nstep) trong do h la do lon buoc chia theo thoi gian va nstep la so buoc chia theo thoi gian.
# ic: Dieu kien ban dau ic=(cvbd,vtbd) trong do cvbd va vtbd tuong ung la vecto chuyen vi ban dau va vecto van toc ban dau
# p: Vecto tai trong p=([t_1; t_2; ...; t_n],[p(t_1);p(t_2);p(t_3);...;p(t_n)]) trong do p(t_i) la gia tri cua tai trong o thoi diem t_i. Gia tri t_1 bao gio cung bang 0. Khi cac thoi diem nay cach deu nhau, chi can vao 1 gia tri la delta_t o vecto dau tien (vi du bang 0.02).

	function mt4(fed,pse);
		hsf,hse,hsd=fed;
		pe=pse[2];
		return mt1\(pe-k*(hsf+hse*h+hsd*h2)-c*(2hsd*h+hse)-m*(2hsd))/h;
	end;

	function mt11(fed,pse);
		hsf,hse,hsd=fed;
		return k*hsd+3c*mt4(fed,pse);
	end;

	function mt12(fed,pse);
		hsf,hse,hsd=fed;
		ps,pe=pse;
		return (k*hse+2c*hsd+6m*mt4(fed,pse)+(ps-pe)/h);
	end;

	function mt13(fed,pse);
		hsf,hse,hsd=fed;
		ps=pse[1];
		return k*hsf+c*hse+2m*hsd-ps;
	end;

	function xdhs_tk(fed,pse);

		sh1=(k2*mt4(fed,pse))*h6/9;
		sh2=(k'*mt11(fed,pse)+5c'*k*mt4(fed,pse))*h5/8;
		sh3=(k'*mt12(fed,pse)+5c'*mt11(fed,pse)+mt5'*k*mt4(fed,pse))*h4/7;
		sh4=(k'*mt13(fed,pse)+5c'*mt12(fed,pse)+mt5'*mt11(fed,pse)+mt6'*k*mt4(fed,pse))*h3/6;
		sh5=(5c'*mt13(fed,pse)+mt5'*mt12(fed,pse)+mt6'*mt11(fed,pse)+mt7'*k*mt4(fed,pse))*h2/5;
		sh6=(mt5'*mt13(fed,pse)+mt6'*mt12(fed,pse)+mt7'*mt11(fed,pse))*h/4;
		sh7=(mt6'*mt13(fed,pse)+mt7'*mt12(fed,pse))/3;
		sh8=(mt7'*mt13(fed,pse))/(2h);

		rn=-(sh1+sh2+sh3+sh4+sh5+sh6+sh7+sh8);

		sh1=(k2*mt4(fed,pse))*h5/8;
		sh2=(k'*mt11(fed,pse)+mt8'*k*mt4(fed,pse))*h4/7;
		sh3=(k'*mt12(fed,pse)+mt8'*mt11(fed,pse)+mt9'*k*mt4(fed,pse))*h3/6;
		sh4=(k'*mt13(fed,pse)+mt8'*mt12(fed,pse)+mt9'*mt11(fed,pse)+mt10'*k*mt4(fed,pse))*h2/5;
		sh5=(mt8'*mt13(fed,pse)+mt9'*mt12(fed,pse)+mt10'*mt11(fed,pse))*h/4;
		sh6=(mt9'*mt13(fed,pse)+mt10'*mt12(fed,pse))/3;
		sh7=(mt10'*mt13(fed,pse))/(2h);

		sn=-(sh1+sh2+sh3+sh4+sh5+sh6+sh7);

		b=[rn ;sn]

		aob=a\b;
		return aob[1:siz2], aob[siz2+1:2*siz2]

	end

	m,c,k=dth;
	h=tz[1];
	nstep=Int(tz[2]);
	hns=h*nstep;

	siz2=size(p[2],2);
	pt,t=ltd_dqc(tz,p);
	xxtt=[Spline1D(t,pt[:,i]) for i=1:siz2];

	u=zeros(nstep+1,siz2);
	ud=zeros(nstep+1,siz2);
	udd=zeros(nstep+1,siz2);

	u[1,:],ud[1,:]=ic;

	c2=c^2;
	c3=c^3;
	c4=c^4;
	c5=c^5;

	h2=h^2;
	h3=h^3;
	h4=h^4;
	h5=h^5;
	h6=h^6;
	h7=h^7;
	h8=h^8;
	h9=h^9;

	k2=k^2;
	k3=k^3;
	k4=k^4;
	k5=k^5;

	m2=m^2;
	m3=m^3;
	m4=m^4;

	c1k1=c*k;
	c1k2=c*k2;
	c1k3=c*k3;
	c2k1=c2*k;
	c2k2=c2*k2;
	c3k1=c3*k;

	kcck=k'*c+c'*k;

	mt1=6m+3h*c+h2*k;
	mt2=mt1\((-20m-5h*c-h2*k)*h2);
	mt3=mt1\((-12m-4h*c-h2*k)*h);
	mt5=(20m+k*mt2);
	mt6=(3c*mt2);
	mt7=(6m*mt2);
	mt8=(4c+k*mt3);
	mt9=(12m+3c*mt3);
	mt10=(6m*mt3);

	sh1=k2*h8/11;
	sh2=kcck*h7/2;
	sh3=(k'*mt5+mt5'*k+25c2)*h6/9;
	sh4=(k'*mt6+mt6'*k+5c'*mt5+5mt5'*c)*h5/8;
	sh5=(k'*mt7+mt7'*k+5c'*mt6+5mt6'*c+mt5'*mt5)*h4/7;
	sh6=(5c'*mt7+5mt7'*c+mt5'*mt6+mt6'*mt5)*h3/6;
	sh7=(mt6'*mt6+mt5'*mt7+mt7'*mt5)*h2/5;
	sh8=(mt6'*mt7+mt7'*mt6)*h/4;
	sh9=(mt7'*mt7)/3;

	q11=sh1+sh2+sh3+sh4+sh5+sh6+sh7+sh8+sh9;

	sh1=k2*h6/9;
	sh2=(k'*mt8+mt8'*k)*h5/8;
	sh3=(k'*mt9+mt9'*k+mt8'*mt8)*h4/7;
	sh4=(k'*mt10+mt10'*k+mt8'*mt9+mt9'*mt8)*h3/6;
	sh5=(mt8'*mt10+mt10'*mt8+mt9'*mt9)*h2/5;
	sh6=(mt9'*mt10+mt10'*mt9)*h/4;
	sh7=(mt10'*mt10)/3;

	q22=sh1+sh2+sh3+sh4+sh5+sh6+sh7;

	sh1=k2*h7/10;
	sh2=(k'*mt8+5c'*k)*h6/9;
	sh3=(k'*mt9+5c'*mt8+mt5'*k)*h5/8;
	sh4=(k'*mt10+5c'*mt9+mt5'*mt8+mt6'*k)*h4/7;
	sh5=(5c'*mt10+mt5'*mt9+mt6'*mt8+mt7'*k)*h3/6;
	sh6=(mt5'*mt10+mt6'*mt9+mt7'*mt8)*h2/5;
	sh7=(mt6'*mt10+mt7'*mt9)*h/4;
	sh8=(mt7'*mt10)/3;

	q12=sh1+sh2+sh3+sh4+sh5+sh6+sh7+sh8;
	q21=q12';

	a=[q11 q12; q21 q22];

    for i=1:nstep
    	ps=[xxtt[j](h*(i-1)) for j=1:siz2];
    	pe=[xxtt[j](h*i) for j=1:siz2];

        hsf=u[i,:];
        hse=ud[i,:];
        hsd=m\(ps-c*hse-k*hsf)/2;
        fed=(hsf,hse,hsd);
        pse=(ps,pe);
        hsa,hsb=xdhs_tk(fed,pse);
        hsc=mt1\((pe-m*(20hsa*h3+12hsb*h2+2hsd)-c*(5hsa*h4+4hsb*h3+2hsd*h+hse)-k*(hsa*h5+hsb*h4+hsd*h2+hse*h+hsf))/h);
        u[i+1,:]=hsa*h5+hsb*h4+hsc*h3+hsd*h2+hse*h+hsf;
        ud[i+1,:]=5hsa*h4+4hsb*h3+3hsc*h2+2hsd*h+hse;
        udd[i+1,:]=m\(pe-c*ud[i+1,:]-k*u[i+1,:]);
    end
    return u,ud,udd
end

function razavi(dth=([2 0;0 1],[0 0; 0 0],[96 -32;-32 32]),tz=(0.1,20),ic=([0;0],[0;0]),p=(0.1,ones(15,1)*[0 100]))
# Giai phuong trinh dao dong HE TUYEN TINH theo phuong phuong phap Razavi
# dth: Dac trung cua he dth=(m,c,k) trong do m la khoi luong, c la he so can nhot va k la do cung
# tz: vecto co cac thong tin ve thoi gian, nhu sau tz=[h;nstep] trong do h la do lon buoc chia theo thoi gian va nstep la so buoc chia theo thoi gian.
# ic: Dieu kien ban dau ic=(cvbd,vtbd) trong do cvbd va vtbd tuong ung la vecto chuyen vi ban dau va vecto van toc ban dau
# p: Vecto tai trong p=([t_1; t_2; ...; t_n],[p(t_1);p(t_2);p(t_3);...;p(t_n)]) trong do p(t_i) la gia tri cua tai trong o thoi diem t_i. Gia tri t_1 bao gio cung bang 0. Khi cac thoi diem nay cach deu nhau, chi can vao 1 gia tri la delta_t o vecto dau tien (vi du bang 0.02).
	function mths_ra(ced,pse);
		ps,pe=pse;
		hsc,hse,hsd=ced;
	    a=[12mh2+4ch3+kh4 6mh+3ch2+kh3; 4mh3+ch4+kh5/5  3mh2+ch3+kh4/4];
	    b=[pe-(2m+2ch+kh2)*hsc-(c+kh)*hsd-k*hse;(ps+pe)*h/2-(2mh+ch2+kh3/3)*hsc-(ch+kh2/2)*hsd-kh*hse];
	    aob=a\b;
	    return aob[1:siz2], aob[siz2+1:2*siz2]
	end
	m,c,k=dth;
	h=tz[1];
	nstep=Int(tz[2]);

	siz2=size(p[2],2);
	pt,t=ltd_dqc(tz,p);
	xxtt=[Spline1D(t,pt[:,i]) for i=1:siz2];

	u=zeros(nstep+1,siz2);
	ud=zeros(nstep+1,siz2);
	udd=zeros(nstep+1,siz2);

	u[1,:],ud[1,:]=ic;

	h2=h^2;
	h3=h^3;
	h4=h^4;
	h5=h^5;
	mh=m*h;
	mh2=m*h2;
	mh3=m*h3;
	ch=c*h;
	ch2=c*h2;
	ch3=c*h3;
	ch4=c*h4;
	kh=k*h;
	kh2=k*h2;
	kh3=k*h3;
	kh4=k*h4;
	kh5=k*h5;

    for i=1:nstep
       	ps=[xxtt[j](h*(i-1)) for j=1:siz2];
    	pe=[xxtt[j](h*i) for j=1:siz2];

    	hsd=ud[i,:];
        hse=u[i,:];

	    hsc=(m\(ps-c*hsd-k*hse))/2;
        ced=(hsc,hse,hsd);

        pse=(ps,pe);

        hsa,hsb=mths_ra(ced,pse);
        u[i+1,:]=hsa*h4+hsb*h3+hsc*h2+hsd*h+hse;
        ud[i+1,:]=4hsa*h3+3hsb*h2+2hsc*h+hsd;
        udd[i+1,:]=m\(pe-c*ud[i+1,:]-k*u[i+1,:]);
    end
    return u,ud,udd
end

function newmark(dth=([2 0;0 1],[0 0; 0 0],[96 -32;-32 32]),tz=(0.1,20),ic=([0;0],[0;0]),p=(0.1,ones(15,1)*[0 100]),param=(1/6,1/2))

# Giai phuong trinh dao dong HE TUYEN TINH theo phuong phuong phap Newmark
# dth: Dac trung cua he dth=(m,c,k) trong do m la khoi luong, c la he so can nhot va k la do cung
# tz: vecto co cac thong tin ve thoi gian, nhu sau tz=[h;nstep] trong do h la do lon buoc chia theo thoi gian va nstep la so buoc chia theo thoi gian.
# ic: Dieu kien ban dau ic=(cvbd,vtbd) trong do cvbd va vtbd tuong ung la vecto chuyen vi ban dau va vecto van toc ban dau
# p: Vecto tai trong p=([t_1; t_2; ...; t_n],[p(t_1);p(t_2);p(t_3);...;p(t_n)]) trong do p(t_i) la gia tri cua tai trong o thoi diem t_i. Gia tri t_1 bao gio cung bang 0. Khi cac thoi diem nay cach deu nhau, chi can vao 1 gia tri la delta_t o vecto dau tien (vi du bang 0.02).

	m,c,k=dth;
	ps=p[1];
	h=tz[1];
	nstep=Int(tz[2]);
	h2=h^2;
	beta,gamma=param;
	bh=beta*h;
	bh2=beta*h2;
	hg2b=h*(gamma/(2beta)-1);
	gobh=gamma/bh;
	gob=gamma/beta;

	siz2=size(p[2],2);
	pt,t=ltd_dqc(tz,p);
	xxtt=[Spline1D(t,pt[:,i]) for i=1:siz2];

	u=zeros(nstep+1,siz2);
	ud=zeros(nstep+1,siz2);
	udd=zeros(nstep+1,siz2);

	u[1,:],ud[1,:]=ic;

    udd[1,:]=m\(p[2][1,:]-c*ud[1,:]-k*u[1,:]);
    k_hat=k+gobh*c+m/bh2;
    a=m/bh+gob*c;
    b=m/(2beta)+hg2b*c;
    for i=1:nstep
       	ps=[xxtt[j](h*(i-1)) for j=1:siz2];
    	pe=[xxtt[j](h*i) for j=1:siz2];

        delta_p=pe-ps;
        delta_p_hat=delta_p+a*ud[i,:]+b*udd[i,:];
        delta_u=k_hat\delta_p_hat;
        delta_udot=gobh*delta_u-gob*ud[i,:]-hg2b*udd[i,:];
        delta_uddot=delta_u/bh2-ud[i,:]/bh-udd[i,:]/(2beta);
        u[i+1,:]=u[i,:]+delta_u;
        ud[i+1,:]=ud[i,:]+delta_udot;
        udd[i+1,:]=udd[i,:]+delta_uddot;
    end
    return u,ud,udd
end

function wilson(dth=([2 0;0 1],[0 0; 0 0],[96 -32;-32 32]),tz=(0.1,20),ic=([0;0],[0;0]),p=(0.1,ones(15,1)*[0 100]),param=1.4)
# Giai phuong trinh dao dong HE TUYEN TINH theo phuong phuong phap Newmark
# dth: Dac trung cua he dth=[m;c;k] trong do m la khoi luong, c la he so can nhot va k la do cung
# tz: vecto co cac thong tin ve thoi gian, nhu sau tz=[h;nstep] trong do h la do lon buoc chia theo thoi gian va nstep la so buoc chia theo thoi gian.
# ic: Dieu kien ban dau ic=[cvbd;vtbd] trong do cvbd va vtbd tuong ung la chuyen vi ban dau va van toc ban dau
# p: Vecto tai trong p=[p(0);p(t1);p(t2);...;p(t_n)] trong do p(t_i) la gia tri cua tai trong o thoi diem t_i
# param: khai bao he so theta cua phuong phap
	m,c,k=dth;
	ps=p[1];
	h=tz[1];
	nstep=Int(tz[2]);
	h2=h^2;
	theta=param;
	thh=theta*h;
	thh2=thh^2;
    hs1=6m/thh2+3c/thh+k;

	siz2=size(p[2],2);
	pt,t=ltd_dqc(tz,p);
	xxtt=[Spline1D(t,pt[:,i]) for i=1:siz2];

	u=zeros(nstep+1,siz2);
	ud=zeros(nstep+1,siz2);
	udd=zeros(nstep+1,siz2);

	u[1,:],ud[1,:]=ic;

    udd[1,:]=m\(p[2][1,:]-c*ud[1,:]-k*u[1,:]);

    for i=1:nstep
       	ps=[xxtt[j](h*(i-1)) for j=1:siz2];
    	pe=[xxtt[j](h*i) for j=1:siz2];
    	pith=ps*(1-theta)+pe*theta;
        shtd1=m*(6u[i,:]/thh2+6ud[i,:]/thh+2udd[i,:])+c*(3u[i,:]/thh+2ud[i,:]+thh*udd[i,:]/2)+pith;
        uith=hs1\shtd1;
        uddith=6*(uith-u[i,:]-thh*ud[i,:]-thh2*udd[i,:]/3)/thh2;
        udd[i+1,:]=udd[i,:]+(uddith-udd[i,:])/theta;
        ud[i+1,:]=ud[i,:]+h*(udd[i,:]+udd[i+1,:])/2;
        u[i+1,:]=u[i,:]+h*ud[i,:]+h2*udd[i,:]/3+h2*udd[i+1,:]/6;
    end
    return u,ud,udd
end

function rk_bi(dth,ccd,tz,ic,ztn_raw);
	function ftcd(t,u,du)
		temp1=Int((length(u)-1)/2)-1;
		temp2=length(u)-2;
		ptt=[];
		for i=1:siz2
			a=xxtt[i](t);
			ptt=[ptt;a]
		end
		ztn=xxztn(t);
		x1=u[1:temp1];
		x2=u[temp1+1:temp1+1];
		v1=u[temp1+2:temp2];
		v2=u[temp2+1:temp2+1];
		w=u[temp2+2:temp2+2];
		du[1:temp1]=v1;
		du[temp1+1:temp1+1]=v2;
		du[temp1+2:temp2]=m\(-c*(v1-iot*v2)-k*(x1-iot*x2))-iot*ztn;
		du[temp2+1:temp2+1]=-(alp*ke*x2+(1-alp)*ke*uy*w-c_s*(delt1*v1-v2)-k_s*(delt1*x1-x2))/m_b-ztn;
		du[temp2+2:temp2+2]=(v2-gamm*abs(v2).*abs(w).^(n-1).*w-beta1*v2.*abs(w).^n)/uy;
	end

	m,c,k=dth;
	alp,beta1,gamm,n,uy,ke=ccd;
	h=tz[1];
	nstep=Int(tz[2])

	siz2=size(m,2);
	iot=ones(siz2,1);
	delt1=[1 zeros(1,siz2-1)];
	sil=(ztn_raw[1],-ztn_raw[2]*iot'*m);
	sil
	silt,t=ltd_dqc(tz,sil);
	ztnt,t1=ltd_dqc(tz,ztn_raw);
	xxtt=[Spline1D(t,silt[:,i]) for i=1:siz2];
	xxztn=Spline1D(t1,ztnt[:]);

	hns=h*nstep;

	u0=ic;

	prob=ODEProblem(ftcd,u0,(0.0,hns));

	sol=solve(prob,Tsit5());
	return sol;
end

end


text:
############################################################################################
### Types and Constructor
############################################################################################
"""
Connects to `Optim.jl` as the optimization backend.

# Constructor

    SemOptimizerOptim(;
        algorithm = LBFGS(),
        options = Optim.Options(;f_tol = 1e-10, x_tol = 1.5e-8),
        kwargs...)

# Arguments
- `algorithm`: optimization algorithm.
- `options::Optim.Options`: options for the optimization algorithm

# Usage
All algorithms and options from the Optim.jl library are available, for more information see
the Optim.jl online documentation.

# Examples
```julia
my_optimizer = SemOptimizerOptim()

# hessian based optimization with backtracking linesearch and modified initial step size
using Optim, LineSearches

my_newton_optimizer = SemOptimizerOptim(
    algorithm = Newton(
        ;linesearch = BackTracking(order=3),
        alphaguess = InitialHagerZhang()
    )
)
```

# Extended help

## Interfaces
- `algorithm(::SemOptimizerOptim)`
- `options(::SemOptimizerOptim)`

## Implementation

Subtype of `SemOptimizer`.
"""
mutable struct SemOptimizerOptim{A, B} <: SemOptimizer
    algorithm::A
    options::B
end

SemOptimizerOptim(;
    algorithm = LBFGS(),
    options = Optim.Options(;f_tol = 1e-10, x_tol = 1.5e-8),
    kwargs...) =
    SemOptimizerOptim(algorithm, options)

############################################################################################
### Recommended methods
############################################################################################

update_observed(optimizer::SemOptimizerOptim, observed::SemObserved; kwargs...) = optimizer

############################################################################################
### additional methods
############################################################################################

algorithm(optimizer::SemOptimizerOptim) = optimizer.algorithm
options(optimizer::SemOptimizerOptim) = optimizer.options


text:
function functional_test()
    Logging.info( "Manager Functional Tests passing. Good job.")
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
module FiniteHorizonValueIteration

using POMDPs
using FiniteHorizonPOMDPs
using POMDPModelTools
import POMDPs: Solver, solve, Policy, action, value
import POMDPLinter: @POMDP_require, @req, @subreq
import ProgressMeter: @showprogress

export
    FiniteHorizonSolver,
    FiniteHorizonPolicy,
    solve,
    action,
    value

include("valueiteration.jl")
include("solver.jl")

end


text:
"""
`SVPMRunTime(pfield, t, dt)`

[description]

Psuedocode:
SVPMRunTime
    Add particles on body
    Calc: U_slip = U(pfield) + U_inf
    Calc VS Strength
    Diffuse vorticity
end

# ARGUMENTS
* `pfield::ParticleField`   : SimpleVPM ParticleField
* `t::Real`                 : Current time
* `dt::Real`                : Time step
* `body::Array{Array{T,1}}` : Matrix of points representing physical boundary

# OUTPUTS
* `breakflag::Bool`       : Flag for stopping SVPM simulation
"""
function SVPMRunTime(pfield::SVPM.ParticleField,
                     t::Real,
                     dt::Real)
                     #body::Array{Array{T,1},1}) where{T<:Real}

    # === Add Particles to field around Boundary ===

    # === Calculate U_slip ===

    # === Calculate Vortex Sheet strength ===

    # === Diffuse vorticity ===

    # === Returns ===

    return true
end


text:
immutable WindowsPath <: AbstractPath
    parts::Tuple
    drive::ASCIIString
end

function WindowsPath(str::AbstractString)
    if isempty(str)
        return WindowsPath(tuple("."), "")
    end

    drive, path = splitdir(str)
    tokenized = split(path, WIN_PATH_SEPARATOR)

    if isempty(tokenized[1])
        tokenized[1] = WIN_PATH_SEPARATOR
    end

    return WindowsPath(tuple(tokenized...), drive)
end

# The following should be implemented in the concrete types
==(a::WindowsPath, b::WindowsPath) = parts(a) == parts(b) && drive(a) == drive(b)
Base.string(path::WindowsPath) = joinpath(parts(path)...)
parts(path::WindowsPath) = path.parts
drive(path::WindowsPath) = path.drive

function isabs(path::WindowsPath)
    if parts(path[1]) == WIN_PATH_SEPARATOR && !isempty(drive(path))
        return true
    else
        return false
    end
end

function root(path::WindowsPath)
    if parts(path)[1] == WIN_PATH_SEPARATOR
        return WIN_PATH_SEPARATOR
    else
        return ""
    end
end

expanduser(path::WindowsPath) = path


text:
@testset "recv_special.jl" begin
    @testset begin receiver = TCPReceiver(4000)
        @test window_size(receiver) == 4000
        @test ackno(receiver) === nothing
        @test unassembled_bytes(receiver) == 0
        segment_arrives(receiver; with_syn=true, result=SegmentArrives_Result_OK)
        @test ackno(receiver).val == 1
        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 0
    end
    @testset begin receiver = TCPReceiver(2358)
        isn = rand(UInt32(0):typemax(UInt32))

        segment_arrives(receiver; data="hello", seqno=isn + 1, result=SegmentArrives_Result_NOT_SYN)

        @test read!(stream_out(receiver)) == ""
        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 0

        segment_arrives(receiver; with_syn=true, seqno=isn, result=SegmentArrives_Result_OK)
        @test ackno(receiver) == WrappingInt32(isn + 1)
    end

    @testset begin receiver = TCPReceiver(2358)
        isn = rand(UInt32(0):typemax(UInt32))

        segment_arrives(receiver; with_syn=true, data="Hello, CS144!", seqno=isn, result=SegmentArrives_Result_OK)

        @test read!(stream_out(receiver)) == "Hello, CS144!"
        @test unassembled_bytes(receiver) == 0
        @test ackno(receiver) == WrappingInt32(isn + 14)
    end

    @testset begin receiver = TCPReceiver(2358)
        isn = rand(UInt32(0):typemax(UInt32))

        segment_arrives(receiver; with_syn=true, seqno=isn, result=SegmentArrives_Result_OK)

        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 0
        @test ackno(receiver) == WrappingInt32(isn + 1)

        segment_arrives(receiver; with_syn=true, seqno=isn + 1, result=SegmentArrives_Result_OK)

        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 0
        @test !stream_out(receiver).has_eof

        segment_arrives(receiver; with_syn=true, seqno=isn + 5, result=SegmentArrives_Result_OK)

        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 0

        @test !stream_out(receiver).has_eof
    end

    @testset begin receiver = TCPReceiver(4000)
        isn = rand(UInt32(0):typemax(UInt32))
        text = "Here's a null byte:" * '\0' * "and it's gone.";
        segment_arrives(receiver; with_syn=true, seqno=isn, result=SegmentArrives_Result_OK)
        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 0

        segment_arrives(receiver; data=text, seqno=isn + 1, result=SegmentArrives_Result_OK)

        @test ackno(receiver) == WrappingInt32(isn + 35)
        @test read!(stream_out(receiver)) == text
        @test !stream_out(receiver).has_eof
    end

    @testset begin receiver = TCPReceiver(2358)
        isn = rand(UInt32(0):typemax(UInt32))

        segment_arrives(receiver; with_syn=true, data="Goodbye, CS144!", with_fin=true, seqno=isn, result=SegmentArrives_Result_OK)

        @test read!(stream_out(receiver)) == "Goodbye, CS144!"
        @test unassembled_bytes(receiver) == 0
        @test ackno(receiver) == WrappingInt32(isn + 17)
        @test eof(stream_out(receiver))
    end

    @testset begin receiver = TCPReceiver(2358)
        isn = rand(UInt32(0):typemax(UInt32))

        segment_arrives(receiver; with_syn=true, seqno=isn, result=SegmentArrives_Result_OK)
        segment_arrives(receiver; data="oodbye, CS144!", with_fin=true, seqno=isn + 2, result=SegmentArrives_Result_OK)

        @test ackno(receiver) == WrappingInt32(isn + 1)
        @test !stream_out(receiver).has_eof
        @test read!(stream_out(receiver)) == ""

        segment_arrives(receiver; data="G", seqno=isn + 1, result=SegmentArrives_Result_OK)

        @test read!(stream_out(receiver)) == "Goodbye, CS144!"
        @test ackno(receiver) == WrappingInt32(isn + 17)
        @test eof(stream_out(receiver))
    end

end

main()

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
for method = methods_list
  str = "Using " * String(method)
  @testset "$str" begin
      @testset "bad local_alg" begin
          @test_logs (:warn, r"Using local_alg .* may result in wrong output\.") get_endpoint(
              [3., 2., 2.1],
              1,
              f_3p_1im_dep,
              method;
              loss_crit = 9.,
              local_alg = :LN_BOBYQA
              )

          @test_throws ArgumentError get_endpoint(
              [3., 2., 2.1],
              1,
              f_3p_1im_dep,
              method;
              loss_crit = 9.,
              local_alg = :xxx
              )
      end

      @testset "warn for LN_NELDERMEAD" begin
          @test_logs (:warn, r"Close-to-zero parameters found when using :LN_NELDERMEAD\.") get_endpoint(
              [3., 0., 2.1],
              1,
              f_3p_1im_dep,
              method;
              loss_crit = 90.,
              local_alg = :LN_NELDERMEAD
              )
      end
  end
end


text:
module test_ctriangulate_catch_exit
using Triangulate

function test()
    nodes=Matrix{Cdouble}([1.0 0.0 ; 0.0 1.0 ;])
    triin=Triangulate.CTriangulateIO()
    triout=Triangulate.CTriangulateIO()
    vorout=Triangulate.CTriangulateIO()

    triin.numberofpoints=Cint(size(nodes,2))
    triin.pointlist=pointer(nodes)
    rc=triangulate("Q",triin,triout,vorout)
    rc==1
end
end


text:
type ScalarContext <: DirectContext  # de-vectorized scalar for-loop
end

##########################################################################
#
#   array access helpers
#
##########################################################################

# const

get_value(r::Number, ::Int) = r
get_value(r::Number, ::Int, ::Int) = r

# vector reader

get_value{T<:Number}(r::Array{T}, i::Int) = r[i]
get_value{T<:Number}(r::Array{T}, i::Int, j::Int) = r[i, j]
get_value{T<:Number}(r::AbstractArray{T}, i::Int) = r[i]
get_value{T<:Number}(r::AbstractArray{T}, i::Int, j::Int) = r[i, j]

##########################################################################
#
#   Assignment composition
#
##########################################################################

function compose_kernel(ctx::ScalarContext, ex::TAssign, lhs_info, rhs_info, i::Symbol)

    lhs_kernel = compose_lhs_kernel(ctx, ex.lhs, lhs_info, i)
    rhs_kernel = compose_rhs_kernel(ctx, ex.rhs, rhs_info, i)
    assignment(lhs_kernel, rhs_kernel)
end

function compose_kernel(ctx::ScalarContext, ex::TAssign, lhs_info, rhs_info, i::Symbol, j::Symbol)

    lhs_kernel = compose_lhs_kernel(ctx, ex.lhs, lhs_info, i, j)
    rhs_kernel = compose_rhs_kernel(ctx, ex.rhs, rhs_info, i, j)
    assignment(lhs_kernel, rhs_kernel)
end

##########################################################################
#
#   LHS composition
#
##########################################################################

# general setup routine for symbol replacement

function with_sym_replacement(ex::TGeneralVar)
    if isa(ex, TVar)
        (nothing, ex.name)
    else
        t = gensym("t")
        pre = assignment(t, ju_expr(ex))
        (pre, t)
    end
end

# TVar

function setup_lhs(ctx::ScalarContext, ex::TVar)
    (nothing, nothing)
end

@compat function length_getter(ctx::ScalarContext, lhs::TVar, info::Void)
    fun_call(:length, lhs.name)
end

@compat function size2d_getter(ctx::ScalarContext, lhs::TVar, info::Void)
    fun_call(qname(:to_size2d), fun_call(:size, lhs.name))
end

@compat function compose_lhs_kernel(ctx::ScalarContext, ex::TVar, info::Void, i::Symbol)
    :( $(ex.name)[$i] )
end

@compat function compose_lhs_kernel(ctx::ScalarContext, ex::TVar, info::Void, i::Symbol, j::Symbol)
    :( $(ex.name)[$i, $j] )
end

# TRef1D

function setup_lhs(ctx::ScalarContext, ex::TRef1D)
    with_sym_replacement(ex.host)
end

function length_getter(ctx::ScalarContext, lhs::TRef1D, h::Symbol)
    length_inference(h, lhs.rgn)
end

function compose_lhs_kernel(ctx::ScalarContext, ex::TRef1D, h::Symbol, i::Symbol)
    :( $(h)[$(indexer(ex.rgn, i))] )
end

# TRefCol

function setup_lhs(ctx::ScalarContext, ex::TRefCol)
    (hpre, h) = with_sym_replacement(ex.host)
    if isa(ex.icol, Int)
        icol = ex.icol
        pre = hpre
    else
        icol = gensym("icol")
        pre = code_block(hpre,
            assignment(icol, fun_call(:convert, :Int, ex.icol)) )
    end
    (pre, (h, icol))
end

function length_getter(ctx::ScalarContext, lhs::TRefCol, info)
    baselen = fun_call(:size, info[1], 1)
    length_inference_(baselen, lhs.rrgn)
end

function compose_lhs_kernel(ctx::ScalarContext, ex::TRefCol, info, i::Symbol)
    h, icol = info
    :( $(h)[$(indexer(ex.rrgn, i)), $icol] )
end

# TRefRow

function setup_lhs(ctx::ScalarContext, ex::TRefRow)
    (hpre, h) = with_sym_replacement(ex.host)
    if isa(ex.irow, Int)
        irow = ex.irow
        pre = hpre
    else
        irow = gensym("irow")
        pre = code_block(hpre,
            assignment(irow, fun_call(:convert, :Int, ex.irow)) )
    end
    (pre, (h, irow))
end

function length_getter(ctx::ScalarContext, lhs::TRefRow, info)
    baselen = fun_call(:size, info[1], 2)
    length_inference_(baselen, lhs.crgn)
end

function compose_lhs_kernel(ctx::ScalarContext, ex::TRefRow, info, i::Symbol)
    h, irow = info
    :( $(h)[$irow, $(indexer(ex.crgn, i))] )
end

# TRef2D

function setup_lhs(ctx::ScalarContext, ex::TRef2D)
    with_sym_replacement(ex.host)
end

function size2d_getter(ctx::ScalarContext, lhs::TRef2D, h::Symbol)
    mat_size_inference(h, lhs.rrgn, lhs.crgn)
end

function compose_lhs_kernel(ctx::ScalarContext, ex::TRef2D, h::Symbol, i::Symbol, j::Symbol)
    :( $(h)[$(indexer(ex.rrgn, i)), $(indexer(ex.crgn, j))] )
end

# TGeneralRef1

function setup_lhs(ctx::ScalarContext, ex::TGeneralRef1)
    hpre, h = with_sym_replacement(ex.host)
    # insert code to convert to integer index
    I = gensym("I")
    pre = code_block(hpre,
        assignment(I, fun_call(:(Base.to_index), ex.i)) )
    (pre, (h, I))
end

function length_getter(ctx::ScalarContext, lhs::TGeneralRef1, info)
    h, I = info
    length_inference(I)
end

function compose_lhs_kernel(ctx::ScalarContext, ex::TGeneralRef1, info, i::Symbol)
    h, I = info
    :( $(h)[$(I)[$(i)]] )
end

##########################################################################
#
#   RHS composition
#
##########################################################################

# TNum

function setup_rhs(ctx::ScalarContext, ex::TNum)
    init = nothing
    siz = :()
    ty = fun_call(:typeof, ex.val)
    final = nothing
    info = nothing
    (init, siz, ty, final, info)
end

@compat function compose_rhs_kernel(ctx::ScalarContext, ex::TNum, info::Void, i::Symbol)
    ex.val
end

@compat function compose_rhs_kernel(ctx::ScalarContext, ex::TNum, info::Void, i::Symbol, j::Symbol)
    ex.val
end

# TVar

function setup_rhs(ctx::ScalarContext, ex::TVar)
    init = nothing
    siz = fun_call(:size, ex.name)
    ty = fun_call(:eltype, ex.name)
    final = nothing
    info = nothing
    (init, siz, ty, final, info)
end

@compat function compose_rhs_kernel(ctx::ScalarContext, ex::TVar, info::Void, i::Symbol)
    fun_call(qname(:get_value), ex.name, i)
end

@compat function compose_rhs_kernel(ctx::ScalarContext, ex::TVar, info::Void, i::Symbol, j::Symbol)
    fun_call(qname(:get_value), ex.name, i, j)
end

# TScalarVar

function setup_rhs(ctx::ScalarContext, ex::TScalarVar)
    init = nothing
    siz = :()
    ty = fun_call(:typeof, ex.name)
    final = nothing
    info = nothing
    (init, siz, ty, final, info)
end

@compat function compose_rhs_kernel(ctx::ScalarContext, ex::TScalarVar, info::Void, i::Symbol)
    ex.name
end

@compat function compose_rhs_kernel(ctx::ScalarContext, ex::TScalarVar, info::Void, i::Symbol, j::Symbol)
    ex.name
end

# other TGeneralVar (e.g. TScalarRef1, TGeneralRef2)

function setup_rhs(ctx::ScalarContext, ex::TGeneralVar)
    t = gensym("t")
    init = assignment(t, ju_expr(ex))
    siz = fun_call(:size, t)
    ty = fun_call(:eltype, t)
    final = nothing
    (init, siz, ty, final, t)
end

function compose_rhs_kernel(ctx::ScalarContext, ex::TGeneralVar, t::Symbol, i::Symbol)
    fun_call(qname(:get_value), t, i)
end

function compose_rhs_kernel(ctx::ScalarContext, ex::TGeneralVar, t::Symbol, i::Symbol, j::Symbol)
    fun_call(qname(:get_value), t, i, j)
end

# TScalarRef1

compose_rhs_kernel(ctx::ScalarContext, ex::TScalarRef1, t::Symbol, i::Symbol) = t

# TGeneralRef1

function setup_rhs(ctx::ScalarContext, ex::TGeneralRef1)
    (init, info) = setup_lhs(ctx, ex)
    (h, I) = info
    siz = fun_call(:size, I)
    ty = fun_call(:eltype, h)
    final = nothing
    (init, siz, ty, final, info)
end

function compose_rhs_kernel(ctx::ScalarContext, ex::TGeneralRef1, info, i::Symbol)
    (h, I) = info
    :( $(h)[$(fun_call(qname(:get_value), I, i))] )
end

# TRef1D, TRefCol, TRefRow

function setup_rhs(ctx::ScalarContext, ex::TRef1D)
    (init, h) = setup_lhs(ctx, ex)
    siz = vec_size_inference(h, ex.rgn)
    ty = fun_call(:eltype, h)
    final = nothing
    (init, siz, ty, final, h)
end

function setup_rhs(ctx::ScalarContext, ex::TRefCol)
    (init, info) = setup_lhs(ctx, ex)
    h = info[1]
    siz = col_size_inference(h, ex.rrgn)
    ty = fun_call(:eltype, h)
    final = nothing
    (init, siz, ty, final, info)
end

function setup_rhs(ctx::ScalarContext, ex::TRefRow)
    (init, info) = setup_lhs(ctx, ex)
    h = info[1]
    siz = row_size_inference(h, ex.crgn)
    ty = fun_call(:eltype, h)
    final = nothing
    (init, siz, ty, final, info)
end

@compat compose_rhs_kernel(ctx::ScalarContext, ex::Union{TRef1D, TRefCol, TRefRow}, info,
    i::Symbol) = compose_lhs_kernel(ctx, ex, info, i)

# TRef2D

function setup_rhs(ctx::ScalarContext, ex::TRef2D)
    (init, h) = setup_lhs(ctx, ex)
    siz = mat_size_inference(h, ex.rrgn, ex.crgn)
    ty = fun_call(:eltype, h)
    final = nothing
    (init, siz, ty, final, h)
end

compose_rhs_kernel(ctx::ScalarContext, ex::TRef2D, h::Symbol,
    i::Symbol, j::Symbol) = compose_lhs_kernel(ctx, ex, h, i, j)

# TMap

function setup_args(ctx::ScalarContext, f::Symbol, args)
    # setup arguments

    arg_setups = [setup_rhs(ctx, a) for a in args]

    arg_inits = [s[1] for s in arg_setups]
    arg_sizes = [s[2] for s in arg_setups]
    arg_types = [s[3] for s in arg_setups]
    arg_finals = [s[4] for s in arg_setups]
    arg_infos = [s[5] for s in arg_setups]

    # integrate

    init = code_block(arg_inits...)
    siz = ewise_size_inference(arg_sizes...)
    ty = result_type_inference(f, arg_types...)
    final = code_block(arg_finals...)

    (init, siz, ty, final, arg_infos)
end

function setup_rhs(ctx::ScalarContext, ex::TMap)
    setup_args(ctx, ex.fun, ex.args)
end

function compose_rhs_kernel(ctx::ScalarContext, ex::TMap, arg_infos, i::Symbol)
    na = length(ex.args)
    arg_kernels = [compose_rhs_kernel(ctx, ex.args[ia], arg_infos[ia], i) for ia in 1 : na]
    fun_call(ex.fun, arg_kernels...)
end

function compose_rhs_kernel(ctx::ScalarContext, ex::TMap, arg_infos, i::Symbol, j::Symbol)
    na = length(ex.args)
    arg_kernels = [compose_rhs_kernel(ctx, ex.args[ia], arg_infos[ia], i, j) for ia in 1 : na]
    fun_call(ex.fun, arg_kernels...)
end

##########################################################################
#
#   main body compilation
#
##########################################################################

# Scalar mode (the result is sure to be a scalar)

function compile(ctx::ScalarContext, mode::ScalarMode, ex::TAssign)
    lhs = ex.lhs
    rhs = ex.rhs
    @assert isa(lhs, TVar) || isa(lhs, TScalarVar)
    assignment(ju_expr(lhs), ju_expr(rhs))
end

# Element-wise (1D or Scalar)

function compile(ctx::ScalarContext, mode::EWiseMode{0}, ex::TAssign)

    lhs = ex.lhs
    rhs = ex.rhs

    @assert isa(lhs, TVar) || isa(lhs, TRef)

    # setup rhs

    (rhs_init, rhs_siz, rhs_ty, rhs_final, rhs_info) = setup_rhs(ctx, rhs)

    # setup lhs

    @gensym siz ty len i
    if isa(lhs, TVar)
        init_lhs = code_block(
            assignment(siz, rhs_siz),
            assignment(ty, rhs_ty),
            assignment(lhs.name, fun_call(:Array, ty, siz))
        )

        lhs_pre = init_lhs
        lhs_info = nothing
    else
        (lhs_pre, lhs_info) = setup_lhs(ctx, lhs)
    end

    get_len = assignment(len, length_getter(ctx, lhs, lhs_info))

    kernel = compose_kernel(ctx, ex, lhs_info, rhs_info, i)

    main_loop = for_statement(i, 1, len, kernel)

    # integrate

    flatten_code_block(
        rhs_init,
        assignment(siz, rhs_siz),
        if_statement( :($siz == ()),
            code_block( # scalar
                assignment(ju_expr(lhs), ju_expr(rhs))
            ),
            flatten_code_block( # ewise 1D
                lhs_pre,
                get_len,
                main_loop,
                rhs_final
            )
        )
    )
end

# Element-wise 1D

function compile(ctx::ScalarContext, mode::EWiseMode{1}, ex::TAssign)

    lhs = ex.lhs
    rhs = ex.rhs

    @assert(isa(lhs, TVar) || isa(lhs, TRef))

    # setup rhs

    (rhs_init, rhs_siz, rhs_ty, rhs_final, rhs_info) = setup_rhs(ctx, rhs)

    # setup lhs

    if isa(lhs, TVar)
        siz = gensym("siz")
        ty = gensym("ty")

        init_lhs = code_block(
            assignment(siz, rhs_siz),
            assignment(ty, rhs_ty),
            assignment(lhs.name, fun_call(:Array, ty, siz))
        )

        lhs_pre = init_lhs
        lhs_info = nothing
    else
        (lhs_pre, lhs_info) = setup_lhs(ctx, lhs)
    end

    len = gensym("len")
    get_len = assignment(len, length_getter(ctx, lhs, lhs_info))

    # kernel & main loop

    i = gensym("i")
    kernel = compose_kernel(ctx, ex, lhs_info, rhs_info, i)

    main_loop = for_statement(i, 1, len, kernel)

    # integrate

    flatten_code_block(
        rhs_init,
        lhs_pre,
        get_len,
        main_loop,
        rhs_final)
end

# Element-wise 2D

function compile(ctx::ScalarContext, mode::EWiseMode{2}, ex::TAssign)

    lhs = ex.lhs
    rhs = ex.rhs

    @assert(isa(lhs, TVar) || isa(lhs, TRef))

    # setup rhs

    (rhs_init, rhs_siz, rhs_ty, rhs_final, rhs_info) = setup_rhs(ctx, rhs)

    # setup lhs

    if isa(lhs, TVar)
        ty = gensym("ty")
        siz = gensym("siz")
        init_lhs = code_block(
            assignment(siz, rhs_siz),
            assignment(ty, rhs_ty),
            assignment(lhs.name, fun_call(:Array, ty, siz))
        )

        lhs_pre = init_lhs
        lhs_info = nothing
    else
        (lhs_pre, lhs_info) = setup_lhs(ctx, lhs)
    end

    m = gensym("m")
    n = gensym("n")
    get_mn = Expr(:(=), :( ($m, $n) ), size2d_getter(ctx, lhs, lhs_info))

    # kernel & main loop

    i = gensym("i")
    j = gensym("j")
    kernel = compose_kernel(ctx, ex, lhs_info, rhs_info, i, j)

    main_loop = for_statement(j, 1, n, for_statement(i, 1, m, kernel))

    # integrate

    flatten_code_block(
        rhs_init,
        lhs_pre,
        get_mn,
        main_loop,
        rhs_final)
end

##########################################################################
#
#   reduction composition & compilation
#
##########################################################################

@compat function setup_reduc(ctx::ScalarContext, r::Union{TReduc, TColwiseReduc, TRowwiseReduc},
    ty::Symbol, siz::Symbol, s::Symbol, rlen::Symbol)

    init_args, infer_siz, infer_ty, final_args, args_info = setup_args(ctx, r.fun, r.args)

    tf = TFun{r.fun}()

    pre = flatten_code_block(
        init_args,
        assignment(siz, infer_siz),
        assignment(ty, infer_ty)
    )

    rinit = reduc_initializer(tf, ty)
    rempty = reduc_emptyval(tf, ty)
    rget = reduc_result_getter(tf, s, rlen)
    (pre, final_args, rinit, rempty, rget, args_info)
end

function compose_fold_kernel(ctx::ScalarContext, f::Symbol, s::Symbol, args, args_info, idxspec...)
    na = length(args)
    @assert na == 1 || na == 2
    tf = TFun{f}()

    if na == 1
        x1 = gensym("x1")
        rupdate = reduc_updater(tf, s, x1)
        a1_kernel = compose_rhs_kernel(ctx, args[1], args_info[1], idxspec...)

        flatten_code_block(
            assignment(x1, a1_kernel),
            rupdate
        )
    else
        x1 = gensym("x1")
        x2 = gensym("x2")
        rupdate = reduc_updater(tf, s, x1, x2)

        a1_kernel = compose_rhs_kernel(ctx, args[1], args_info[1], idxspec...)
        a2_kernel = compose_rhs_kernel(ctx, args[2], args_info[2], idxspec...)

        flatten_code_block(
            assignment(x1, a1_kernel),
            assignment(x2, a2_kernel),
            rupdate
        )
    end
end

function compose_fold_loop(ctx::ScalarContext, rfun::Symbol, s::Symbol, len::Symbol, args, args_info)

    @gensym i
    kernel = compose_fold_kernel(ctx, rfun, s, args, args_info, i)
    for_statement(i, 1, len, kernel)
end

function compose_fold_loop2(ctx::ScalarContext, rfun::Symbol, s::Symbol, m::Symbol, n::Symbol, args, args_info)

    @gensym i j
    kernel = compose_fold_kernel(ctx, rfun, s, args, args_info, i, j)
    for_statement(j, 1, n, for_statement(i, 1, m, kernel))
end

function compile(ctx::ScalarContext, mode::ReducMode, ex::TAssign)

    lhs = ex.lhs
    rhs = ex.rhs
    @assert isa(lhs, TVar) || isa(lhs, TScalarVar)
    @assert isa(rhs, TReduc)

    # generate symbols
    @gensym rlen siz ty
    s = lhs.name # compiler has ensured no alias between lhs and rhs

    # setup
    (pre, final_args, rinit, rempty, rget, args_info) = setup_reduc(ctx, rhs, ty, siz, s, rlen)

    if rget == s
        rpost = nothing
    else
        rpost = assignment(s, rget)
    end

    # main loop
    amode = rhs.arg_mode
    if isa(amode, EWiseMode{0}) || isa(amode, EWiseMode{1})
        main = compose_fold_loop(ctx, rhs.fun, s, rlen, rhs.args, args_info)
    else
        @assert isa(amode, EWiseMode{2})
        m = gensym("m")
        n = gensym("n")
        main_loop = compose_fold_loop2(ctx, rhs.fun, s, m, n, rhs.args, args_info)
        main = code_block(
            :( ($m, $n) = ($siz) ),
            main_loop
        )
    end

    # integrate

    flatten_code_block(
        pre,
        assignment(rlen, fun_call(qname(:to_length), siz)),
        if_statement( :( $rlen > 0 ),
            flatten_code_block(  # not empty
                assignment(s, rinit),
                main,
                rpost
            ),
            flatten_code_block(  # empty
                assignment(s, rempty)
            )
        ),
        final_args
    )
end

# partial reduction

function compile(ctx::ScalarContext, mode::ColwiseReducMode, ex::TAssign)

    lhs = ex.lhs
    rhs = ex.rhs
    @assert isa(lhs, TVar) || isa(lhs, TRef)
    @assert isa(rhs, TColwiseReduc)

    # generate symbols
    @gensym siz ty m n i j s

    # setup rhs
    (pre, final_args, rinit, rempty, rget, args_info) = setup_reduc(ctx, rhs, ty, siz, s, m)

    if rget == s
        rpost = nothing
    else
        rpost = assignment(s, rget)
    end

    # setup lhs
    if isa(lhs, TVar)
        lhs_pre = assignment(lhs.name, fun_call(:Array, ty, 1, n))
        lhs_info = nothing
    else
        (lhs_pre, lhs_info) = setup_lhs(ctx, lhs)
    end

    lhs_kernel = compose_lhs_kernel(ctx, lhs, lhs_info, j)

    # main loop
    kernel = compose_fold_kernel(ctx, rhs.fun, s, rhs.args, args_info, i, j)
    main_loop = for_statement(j, 1, n, flatten_code_block(
        assignment(s, rinit),
        for_statement(i, 1, m, kernel),
        rpost,
        assignment(lhs_kernel, s)
    ))

    # integrate

    flatten_code_block(
        pre,
        :( ($m, $n) = Devectorize.to_size2d($siz) ),
        lhs_pre,
        if_statement( :( $m > 0 ),
            flatten_code_block(  # not empty
                main_loop
            ),
            flatten_code_block(  # empty
                assignment(s, rempty),
                for_statement(j, 1, n, assignment(lhs_kernel, s))
            )
        ),
        final_args
    )
end

function compile(ctx::ScalarContext, mode::RowwiseReducMode, ex::TAssign)

    lhs = ex.lhs
    rhs = ex.rhs
    @assert isa(lhs, TVar) || isa(lhs, TRef)
    @assert isa(rhs, TRowwiseReduc)

    # generate symbols
    @gensym siz ty m n i j s

    # setup rhs
    (pre, final_args, rinit, rempty, rget, args_info) = setup_reduc(ctx, rhs, ty, siz, s, n)

    # setup lhs

    if isa(lhs, TVar)
        lhs_pre = assignment(lhs.name, fun_call(:Array, ty, m, 1))
        lhs_info = nothing
    else
        (lhs_pre, lhs_info) = setup_lhs(ctx, lhs)
    end
    lhs_kernel = compose_lhs_kernel(ctx, lhs, lhs_info, i)

    # loops

    kernel = compose_fold_kernel(ctx, rhs.fun, s, rhs.args, args_info, i, j)

    init_loop = for_statement(i, 1, m, flatten_code_block(
        assignment(s, rinit),
        kernel,
        assignment(lhs_kernel, s)
    ))

    main_loop = for_statement(j, 2, n,
        for_statement(i, 1, m, flatten_code_block(
            assignment(s, lhs_kernel),
            kernel,
            assignment(lhs_kernel, s)
        ))
    )

    if rget == s
        post_loop = nothing
    else
        post_loop = for_statement(i, 1, m, code_block(
            assignment(s, lhs_kernel),
            assignment(s, rget),
            assignment(lhs_kernel, s)
        ))
    end

    # integrate

    flatten_code_block(
        pre,
        :( ($m, $n) = Devectorize.to_size2d($siz) ),
        lhs_pre,
        if_statement( :( $m > 0 ),
            flatten_code_block(  # not empty
                assignment(j, 1),
                init_loop,
                main_loop,
                post_loop
            ),
            flatten_code_block(  # empty
                assignment(s, rempty),
                for_statement(i, 1, m, assignment(lhs_kernel, s))
            )
        ),
        final_args
    )
end

##########################################################################
#
#   code-generating macros
#
##########################################################################

function dump_devec(ex::Expr)
    println(compile(ScalarContext(), ex))
end

macro devec(assign_ex)
    esc(compile(ScalarContext(), assign_ex))
end

macro inspect_devec(assign_ex)
    let code__ = compile(ScalarContext(), assign_ex)
        println("$assign_ex ==>")
        println(code__)
        esc(code__)
    end
end


text:
struct ExponentialSampler <: Sampleable{Univariate,Continuous}
    scale::AbstractFloat
end

rand(rng::AbstractRNG, s::ExponentialSampler) = s.scale * randexp(rng)

struct ExponentialLogUSampler <: Sampleable{Univariate,Continuous}
    scale::AbstractFloat
end

rand(rng::AbstractRNG, s::ExponentialLogUSampler) = -s.scale * log(rand(rng))