text:
# ===============================================
# [1] Concatenate matrices of diferent countries.
# ===============================================

function GSgroups(name,range,nv,h,nrep,cut)
	nc = length(range)
	aux1 = Array{Float32,3}(undef,nv,h,nrep*nc)
	aux2 = Array{Float32,3}(undef,nv,h,nrep*nc)
	aux3 = Array{Float32,3}(undef,nv,h,nrep*nc)
	aux4 = Array{Float32,3}(undef,nv,h,nrep*nc)
	for i in 1:nc
		nam = Symbol.(name.* string(i));
		ex =:(n1 = $nam[3].irfgs);	eval(ex);
		aux1[:,:,nrep*(i-1)+1:nrep*i] = n1;
		ex =:(n1 = $nam[3].fevgs); 	eval(ex);
		aux2[:,:,nrep*(i-1)+1:nrep*i] = n1;
		ex =:(n1 = $nam[3].irfnf); 	eval(ex);
		aux3[:,:,nrep*(i-1)+1:nrep*i] = n1;
		ex =:(n1 = $nam[3].fevnf);  eval(ex);
		aux4[:,:,nrep*(i-1)+1:nrep*i] = n1;
	end
	reddim(x) = dropdims(mean(x, dims = 3), dims = 3)
	Qntls  = GShock.Qntls
	crep = nrep*cut
	crep2 = nrep*(nc-cut)
	qq = [0.16 0.50 0.84]
	g1 = (aux1[:,:,1:crep],aux2[:,:,1:crep],aux3[:,:,1:crep],aux4[:,:,1:crep])
	g2 = (aux1[:,:,crep+1:end],aux2[:,:,crep+1:end],aux3[:,:,crep+1:end],aux4[:,:,crep+1:end])

	ecx = (
		IrfGS = (Mean = reddim(g1[1]), Qntls = Qntls(g1[1], crep, qq, nv, h) ),
		FevGS = (Mean = reddim(g1[2]), Qntls = Qntls(g1[2], crep, qq, nv, h) ),
		IrfNF = (Mean = reddim(g1[3]), Qntls = Qntls(g1[3], crep, qq, nv, h) ),
		FevNF = (Mean = reddim(g1[4]), Qntls = Qntls(g1[4], crep, qq, nv, h) ),
	)
	dcx = (
		IrfGS = (Mean = reddim(g2[1]), Qntls = Qntls(g2[1], crep2, qq, nv, h) ),
		FevGS = (Mean = reddim(g2[2]), Qntls = Qntls(g2[2], crep2, qq, nv, h) ),
		IrfNF = (Mean = reddim(g2[3]), Qntls = Qntls(g2[3], crep2, qq, nv, h) ),
		FevNF = (Mean = reddim(g2[4]), Qntls = Qntls(g2[4], crep2, qq, nv, h) ),
	)
	return ecx, dcx
end
# ===============================================
# [2] Function for graphs.
# ===============================================
function GSGraph(model,	name, labels; colg = :sienna, subdir = "Countries", varI=1, varF = 0)
    nf = nfields(model.FevGS.Qntls)
	if varF==0
		varF = size(model.FevGS.Mean)[1];
	end
    if ~isdir(".//Figures")
        mkdir("Figures");
    end
    if ~isdir(".//Figures//$subdir")
        mkdir(".//Figures//$subdir");
    end
    # ----------------------------------
    # FEVD GLOBAL shocks
	name1 = ".//Figures//$subdir//FEV_" * name;
	data  = model.FevGS;
	ModelGraph(data,nf,varI,varF,name1,labels,colg)
	# ----------------------------------
    # IRF GLOBAL SHOCKS
	name1 = ".//Figures//$subdir//IRF_" * name;
	data  = model.IrfGS;
	ModelGraph(data,nf,varI,varF,name1,labels,colg)
	# ----------------------------------
    # FEVD Non Fundamental shocks
	name1 = ".//Figures//$subdir//FEVNF_" * name;
	data  = model.FevNF;
	ModelGraph(data,nf,varI,varF,name1,labels,colg)
	# ----------------------------------
    # IRF Non Fundamental SHOCKS
	name1 = ".//Figures//$subdir//IRFNF_" * name;
	data  = model.IrfNF;
	ModelGraph(data,nf,varI,varF,name1,labels,colg)
end

function ModelGraph(data,nf,varI,varF,name,labels,colg)
	quint = cat(dims = 3, [data.Qntls[i] for i = 1:nf]...)
	quint = permutedims(quint, [2, 3, 1])
	h, n, nvar = size(quint)
	nrows = Int(ceil((varF-varI+1) / 3))
	if nrows==1
		plot(layout = (nrows, 3),size=(1200,400),title =["Global Output" "Commodity Price" "BAA spread"])
	else
		plot(layout = (nrows, 3))
	end
	j = 1;
	for i = varI:varF
		meanG = data.Mean[i, :]
		line  = quint[:, 2, i]
		LB 	  = line .- quint[:, 1, i]
		UB    = quint[:, 3, i] .- line
		plot!(1:h,
			line;
			ribbon = (LB, UB),
			tickfontsize = 6,
			subplot = j,
			c = [colg],
			xlabel=labels[i],
			w = 1.5,
			fillalpha = 0.2,
			legend = false,
		)
		scatter!(1:h,
			meanG,
			markersize = 1.25,
			c = [colg],
			markershape = :x,
			tickfontsize = 6,
			subplot = j,
		)
		j+=1;
	end
	savefig(name);
end

function GraphAux(raw,name)
	nvar = size(raw)[1];
	h    = size(raw)[2];
	p = plot(legendfontsize=10,layout=(2,3),size=(1200,800),title =["GDP" "Consumption" "Investment" "Trade" "REER" "Monetary Policy"])
	for i in 1:nvar
	        data = raw[i,:,:];
	        if i!=2
	                plot!(1:h,data, label="",
	                        color= [:red :slateblue4], markersize = 3.5,
	                        markershape = [:none :x], w = [2 0.1],
	                        style = [:solid :dot],markerstrokewidth= 0.1,
	                        subplot = i, grid = :false,
	                );
	        else
	                plot!(1:h,data, label=["Global Shock" "Non-fundamental"],
	                        color= [:red :slateblue4], markersize = 3.5,
	                        markershape = [:none :x ], w = [2 0.1],
	                        style = [:solid :dot],markerstrokewidth= 0.1,
	                        subplot=i, fg_legend = :transparent,
							bg_legend = :transparent, grid = :false,
	                );
	        end
	end
	savefig(name);
end

function myplot(data,h,mylabel)
	p1 = plot(1:h,data,
        	label=mylabel, color= [:red :blue :red :blue],
        	markersize = 2.5, markershape = [:none :none :x :x],
        	w = [2 2 0.1 0.1], style = [:solid :solid :dot :dot],
        	markerstrokewidth= 0.1, fg_legend = :transparent,
			bg_legend = :transparent,
        );
	return p1;
end
# ===============================================
# [3] HTML tables
# ===============================================

function ToHtml(file, matt, colnames)
	io = open(file, "w");
	println(io,"
");
	println(io," ", "h", " 	")
	for x in colnames
		println(io," ",x, String(" 	"))
	end
	T, cols = size(matt);
	for i in 1:T
		println(io,"
")
		println(io," ", string(i-1), " 	")
		for j in 1:cols
			println(io," ", string(matt[i,j]), " 	")
		end
		println(io,"\n
")
	end
	println(io,"
")
	close(io);
end

max_stars_repo_name: ChrisVeJa/GlobalShock

max_issues_repo_name: ChrisVeJa/GlobalShock

max_forks_repo_name: ChrisVeJa/GlobalShock


text:
# function visualize!(vis, model::PlanarPush, q, u;
#         r = 1.0,
#         r_pusher = 0.1,
#         Δt = 0.1,
#         contact_points = false,
#         force_vis = false)
#
# 	default_background!(vis)
#
#     r_box = r - r_pusher
#
#     setobject!(vis["box"], GeometryBasics.Rect(Vec(-1.0 * r_box,
# 		-1.0 * r_box,
# 		-1.0 * r_box),
# 		Vec(2.0 * r_box, 2.0 * r_box, 2.0 * r_box)),
# 		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))
#
#     setobject!(vis["pusher"],
#         Cylinder(Point(0.0, 0.0, 0.0), Point(0.0, 0.0, r_box), r_pusher),
#         MeshPhongMaterial(color = RGBA(1.0, 153.0 / 255.0, 51.0 / 255.0, 1.0)))
#
#     if contact_points
#         for i = 1:4
#             setobject!(vis["contact$i"], GeometryBasics.Sphere(Point3f0(0),
#                 convert(Float32, 0.02)),
#                 MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0.0, 1.0)))
#         end
#     end
#
#     if force_vis
# 	       force_vis = ArrowVisualizer(vis[:force])
#     	setobject!(force_vis, MeshPhongMaterial(color=RGBA(1.0, 0.0, 0.0, 1.0)))
#
#         unorm = maximum(abs.(hcat(u...)))
#
#     	us = u[1] / unorm
#
#     	settransform!(force_vis,
#     				Point(q[1][4] - us[1], q[1][5] - us[1], 0.0),
#     				Vec(us[1], us[1], 0.0),
#     				shaft_radius=0.025,
#     				max_head_radius=0.05)
#     end
#
#     anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))
#
# 	T = length(q)
#     for t = 1:T-1
#         MeshCat.atframe(anim, t) do
#
#             if force_vis
#     			if t < T-1
#     				if norm(u[t]) < 1.0e-6
#     					setvisible!(vis[:force], false)
#     				else
#     					setvisible!(vis[:force], true)
#
#     					us = u[t] / unorm
#     					settransform!(force_vis,
#     								Point(q[t+1][4] - us[1], q[t+1][5] - us[2], 0.0),
#     								Vec(us[1], us[2], 0.0),
#                                     shaft_radius=0.025,
#                     				max_head_radius=0.05)
#     				end
#     			end
#             end
#
#             settransform!(vis["box"],
# 				compose(Translation(q[t+1][1], q[t+1][2], r), LinearMap(RotZ(q[t+1][3]))))
#
#             settransform!(vis["pusher"], Translation(q[t+1][4], q[t+1][5], r))
#
#             if contact_points
#                 for i = 1:4
#                     settransform!(vis["contact$i"],
#                         Translation(([q[t+1][1:2]; 0.0] + RotZ(q[t+1][3]) * [contact_corner_offset[i]; 0.0])...))
#                 end
#             end
#         end
#     end
#
# 	settransform!(vis["/Cameras/default"],
# 		compose(Translation(0.0, 0.0, 50.0), LinearMap(RotZ(0.5 * pi) * RotY(-pi/2.5))))
# 	setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 50)
#
#
#     MeshCat.setanimation!(vis, anim)
# end

function _create_planar_push!(vis, model::PlanarPush;
        i = 1,
        r = 1.0,
        r_pusher = 0.1,
        tl = 1.0,
        box_color = RGBA(0.0, 0.0, 0.0, tl),
        pusher_color = RGBA(0.5, 0.5, 0.5, tl))

    r_box = r - r_pusher

    setobject!(vis["box_$i"], GeometryBasics.Rect(Vec(-1.0 * r_box,
		-1.0 * r_box,
		-1.0 * r_box),
		Vec(2.0 * r_box, 2.0 * r_box, 2.0 * r_box)),
		MeshPhongMaterial(color = box_color))

    setobject!(vis["pusher_$i"],
        Cylinder(Point(0.0, 0.0, 0.0), Point(0.0, 0.0, r_box), r_pusher),
        MeshPhongMaterial(color = pusher_color))
end

function _set_planar_push!(vis, model::PlanarPush, q;
    i = 1)
    settransform!(vis["box_$i"],
		compose(Translation(q[1], q[2], 0.01 * i), LinearMap(RotZ(q[3]))))
    settransform!(vis["pusher_$i"], Translation(q[4], q[5], 0.01 * i))
end

function visualize!(vis, model::PlanarPush, q;
        i = 1,
        r = 1.0,
        r_pusher = 0.1,
        tl = 1.0,
        box_color = RGBA(0.0, 0.0, 0.0, tl),
        pusher_color = RGBA(0.5, 0.5, 0.5, tl),
        Δt = 0.1)

	default_background!(vis)

    _create_planar_push!(vis, model,
        i = i,
        r = r,
        r_pusher = r_pusher,
        tl = tl,
        box_color = box_color,
        pusher_color = pusher_color)

    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

	T = length(q)
    for t = 1:T-1
        MeshCat.atframe(anim, t) do
            _set_planar_push!(vis, model, q[t])
        end
    end

	settransform!(vis["/Cameras/default"],
		compose(Translation(0.0, 0.0, 50.0), LinearMap(RotZ(0.5 * pi) * RotY(-pi/2.5))))
	setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 50)

    MeshCat.setanimation!(vis, anim)
end

ext: jl

max_issues_repo_name: thowell/DirectMotionPlanning


text:
mutable struct RMSProp <: Optimizer
    xparams::Vector{XVariable}
    w::Vector
    lr::AbstractFloat
    ϵ::AbstractFloat
    inertia::AbstractFloat
    name::String
    function RMSProp(xparams::Vector{XVariable}; lr=1e-2, inertia=0.99, epsilon=1e-8)
        n = length(xparams)
        w = Vector(undef,n)
        for i = 1:n
            c , θ = xparams[i]
            w[i] = Zeros(typeof(θ.value), θ.shape)
        end
        new(xparams, w, lr, epsilon, inertia, "RMSProp")
    end
end

function Base.show(io::IO, O::RMSProp)
    print("RMSProp(lr=$(O.lr), ϵ=$(O.ϵ), inertia=$(O.inertia))")
end

function update!(O::RMSProp; clipfn::Function=LPInfNormClip, clipvalue=10.0)
    w = O.w
    ϵ = O.ϵ
    μ = - O.lr
    ρ = O.inertia

    for i = 1:length(O.xparams)
        c , θ = O.xparams[i]
        ∇ = clipfn(setNanInfZero(θ.delta), clipvalue)
        @. w[i] = ρ * w[i] + (1-ρ) * ∇ * ∇
        @. θ.value += μ / (sqrt(w[i])+ϵ) * ∇
    end
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
using UnsteadyFlowSolvers
using Test

include("ldvmLinTest.jl")


text:
using JuMP
using JuMP._Derivatives
using LinearAlgebra
using Test

const ForwardDiff = JuMP.ForwardDiff

struct ΦEvaluator <: MOI.AbstractNLPEvaluator end

@testset "Derivatives" begin
    NO_USER_OPS = JuMP._Derivatives.UserOperatorRegistry()

    ex = :(sin(x[1]^2) + cos(x[2] * 4) / 5 - 2.0)

    nd, const_values = expr_to_nodedata(ex)
    adj = adjmat(nd)

    # @show nd

    storage = zeros(length(nd))
    partials_storage = zeros(length(nd))
    reverse_storage = zeros(length(nd))

    x = [2.0, 3.0]
    # @show x
    fval = forward_eval(
        storage,
        partials_storage,
        nd,
        adj,
        const_values,
        [],
        x,
        [],
        [],
        [],
        NO_USER_OPS,
    )
    true_val = sin(x[1]^2) + cos(x[2] * 4) / 5 - 2.0
    @test fval ≈ true_val

    grad = zeros(2)
    reverse_eval(reverse_storage, partials_storage, nd, adj)
    reverse_extract(grad, reverse_storage, nd, adj, [], 1.0)

    true_grad = [2 * x[1] * cos(x[1]^2), -4 * sin(x[2] * 4) / 5]
    @test grad ≈ true_grad

    # Testing view
    xx = [1.0, 2.0, 3.0, 4.0, 5.0]
    xv = @view xx[2:3]
    fval_view = forward_eval(
        storage,
        partials_storage,
        nd,
        adj,
        const_values,
        [],
        xv,
        [],
        [],
        [],
        NO_USER_OPS,
    )
    @test fval_view == fval

    grad_view = zeros(2)
    reverse_eval(reverse_storage, partials_storage, nd, adj)
    reverse_extract(grad_view, reverse_storage, nd, adj, [], 1.0)
    @test grad_view == grad

    grad_view = zeros(5)
    gv = @view grad_view[2:2:end]
    reverse_eval(reverse_storage, partials_storage, nd, adj)
    reverse_extract(gv, reverse_storage, nd, adj, [], 1.0)
    @test gv == grad

    # subexpressions

    nd_outer = [NodeData(SUBEXPRESSION, 1, -1)]
    li = list_subexpressions(nd_outer)
    @test li == [1]
    li, li_individual =
        order_subexpressions(Vector{NodeData}[nd_outer], Vector{NodeData}[nd])
    @test li == [1]
    @test li_individual[1] == [1]

    nd_outer2 = [
        NodeData(CALL, operator_to_id[:+], -1),
        NodeData(SUBEXPRESSION, 2, 1),
        NodeData(SUBEXPRESSION, 1, 1),
    ]
    li = list_subexpressions(nd_outer2)
    @test li == [1, 2]
    li, li_individual = order_subexpressions(
        Vector{NodeData}[nd_outer2],
        Vector{NodeData}[nd, nd_outer],
    )
    @test li == [1, 2]
    @test li_individual[1] == [1, 2]

    @testset "order_subexpressions with nested subexpressions" begin
        expr_order, expr_order_individual = order_subexpressions(
            [
                [NodeData(SUBEXPRESSION, 1, -1)],
                [NodeData(SUBEXPRESSION, 1, -1)],
            ],
            [[NodeData(SUBEXPRESSION, 2, -1)], NodeData[]],
        )
        @test expr_order == [2, 1]
        @test expr_order_individual[1] == [2, 1]
        @test expr_order_individual[2] == [2, 1]
    end

    adj_outer = adjmat(nd_outer)
    outer_storage = zeros(1)
    outer_storage_partials = zeros(1)
    fval = forward_eval(
        outer_storage,
        outer_storage_partials,
        nd_outer,
        adj_outer,
        [],
        [],
        x,
        [fval],
        [],
        [],
        NO_USER_OPS,
    )
    @test fval ≈ true_val

    outer_reverse_storage = zeros(1)
    fill!(grad, 0.0)
    subexpr_output = zeros(1)
    reverse_eval(
        outer_reverse_storage,
        outer_storage_partials,
        nd_outer,
        adj_outer,
    )
    reverse_extract(
        grad,
        outer_reverse_storage,
        nd_outer,
        adj_outer,
        subexpr_output,
        1.0,
    )
    @assert subexpr_output[1] == 1.0
    reverse_eval(reverse_storage, partials_storage, nd, adj)
    reverse_extract(grad, reverse_storage, nd, adj, [], subexpr_output[1])
    @test grad ≈ true_grad

    ex = :((1 / x[1])^x[2] - x[3])

    nd, const_values = expr_to_nodedata(ex)
    adj = adjmat(nd)

    storage = zeros(length(nd))
    partials_storage = zeros(length(nd))
    reverse_storage = zeros(length(nd))

    x = [2.5, 3.5, 1.0]
    # @show x
    fval = forward_eval(
        storage,
        partials_storage,
        nd,
        adj,
        const_values,
        [],
        x,
        [],
        [],
        [],
        NO_USER_OPS,
    )
    true_val = (1 / x[1])^x[2] - x[3]
    @test fval ≈ true_val

    grad = zeros(3)
    reverse_eval(reverse_storage, partials_storage, nd, adj)
    reverse_extract(grad, reverse_storage, nd, adj, [], 1.0)

    true_grad = [-x[2] * x[1]^(-x[2] - 1), -((1 / x[1])^x[2]) * log(x[1]), -1]
    @test grad ≈ true_grad

    # logical expressions
    ex = :(x[1] > 0.5 && x[1] < 0.9)
    nd, const_values = expr_to_nodedata(ex)
    adj = adjmat(nd)
    storage = zeros(length(nd))
    partials_storage = zeros(length(nd))
    fval = forward_eval(
        storage,
        partials_storage,
        nd,
        adj,
        const_values,
        [],
        [1.5],
        [],
        [],
        [],
        NO_USER_OPS,
    )
    @test fval == 0
    fval = forward_eval(
        storage,
        partials_storage,
        nd,
        adj,
        const_values,
        [],
        [0.6],
        [],
        [],
        [],
        NO_USER_OPS,
    )
    @test fval == 1

    ex = :(ifelse(x[1] >= 0.5 || x[1] <= 0.1, x[1], 5))
    nd, const_values = expr_to_nodedata(ex)
    adj = adjmat(nd)
    storage = zeros(length(nd))
    partials_storage = zeros(length(nd))
    reverse_storage = zeros(length(nd))
    grad = zeros(1)
    fval = forward_eval(
        storage,
        partials_storage,
        nd,
        adj,
        const_values,
        [],
        [1.5],
        [],
        [],
        [],
        NO_USER_OPS,
    )
    @test fval == 1.5
    reverse_eval(reverse_storage, partials_storage, nd, adj)
    reverse_extract(grad, reverse_storage, nd, adj, [], 1.0)
    @test grad[1] == 1

    fval = forward_eval(
        storage,
        partials_storage,
        nd,
        adj,
        const_values,
        [],
        [-0.1],
        [],
        [],
        [],
        NO_USER_OPS,
    )
    @test fval == -0.1
    fill!(grad, 0)
    reverse_eval(reverse_storage, partials_storage, nd, adj)
    reverse_extract(grad, reverse_storage, nd, adj, [], 1.0)
    @test grad[1] == 1

    fval = forward_eval(
        storage,
        partials_storage,
        nd,
        adj,
        const_values,
        [],
        [0.2],
        [],
        [],
        [],
        NO_USER_OPS,
    )
    @test fval == 5
    fill!(grad, 0)
    reverse_eval(reverse_storage, partials_storage, nd, adj)
    reverse_extract(grad, reverse_storage, nd, adj, [], 1.0)
    @test grad[1] == 0

    ex = :(sqrt(x[1]))
    nd, const_values = expr_to_nodedata(ex)
    adj = adjmat(nd)
    storage = zeros(length(nd))
    partials_storage = zeros(length(nd))
    reverse_storage = zeros(length(nd))
    grad = zeros(1)
    fval = forward_eval(
        storage,
        partials_storage,
        nd,
        adj,
        const_values,
        [],
        [-1.5],
        [],
        [],
        [],
        NO_USER_OPS,
    )
    @test isnan(fval)
    reverse_eval(reverse_storage, partials_storage, nd, adj)
    reverse_extract(grad, reverse_storage, nd, adj, [], 1.0)
    @test isnan(grad[1])

    # parameters
    nd = [NodeData(PARAMETER, 1, -1)]
    adj = adjmat(nd)
    storage = zeros(length(nd))
    partials_storage = zeros(length(nd))
    fval = forward_eval(
        storage,
        partials_storage,
        nd,
        adj,
        [],
        [105.2],
        [-0.1],
        [],
        [],
        [],
        NO_USER_OPS,
    )
    @test fval == 105.2

    # abs
    ex = :(abs(x[1]))

    nd, const_values = expr_to_nodedata(ex)
    adj = adjmat(nd)

    storage = zeros(length(nd))
    partials_storage = zeros(length(nd))
    reverse_storage = zeros(length(nd))
    fval = forward_eval(
        storage,
        partials_storage,
        nd,
        adj,
        const_values,
        [],
        [2.0],
        [],
        [],
        [],
        NO_USER_OPS,
    )
    @test fval == 2.0
    grad = zeros(1)
    reverse_eval(reverse_storage, partials_storage, nd, adj)
    reverse_extract(grad, reverse_storage, nd, adj, [], 1.0)
    @test grad[1] == 1.0

    fval = forward_eval(
        storage,
        partials_storage,
        nd,
        adj,
        const_values,
        [],
        [-2.0],
        [],
        [],
        [],
        NO_USER_OPS,
    )
    @test fval == 2.0
    grad = zeros(1)
    reverse_eval(reverse_storage, partials_storage, nd, adj)
    reverse_extract(grad, reverse_storage, nd, adj, [], 1.0)
    @test grad[1] == -1.0

    # https://github.com/jump-dev/JuMP.jl/issues/855
    ex = :(ifelse(x[1] <= 3.0, (x[1] - 2.0)^2, 2 * log(x[1] - 2.0) + 1.0))

    nd, const_values = expr_to_nodedata(ex)
    adj = adjmat(nd)

    storage = zeros(length(nd))
    partials_storage = zeros(length(nd))
    reverse_storage = zeros(length(nd))
    fval = forward_eval(
        storage,
        partials_storage,
        nd,
        adj,
        const_values,
        [],
        [-1.0],
        [],
        [],
        [],
        NO_USER_OPS,
    )
    @test fval == 9.0
    grad = zeros(1)
    reverse_eval(reverse_storage, partials_storage, nd, adj)
    reverse_extract(grad, reverse_storage, nd, adj, [], 1.0)
    @test grad[1] == -6.0

    fval = forward_eval(
        storage,
        partials_storage,
        nd,
        adj,
        const_values,
        [],
        [2.0],
        [],
        [],
        [],
        NO_USER_OPS,
    )
    @test fval == 0.0
    grad = zeros(1)
    reverse_eval(reverse_storage, partials_storage, nd, adj)
    reverse_extract(grad, reverse_storage, nd, adj, [], 1.0)
    @test grad[1] == 0.0

    function test_linearity(ex, testval, IJ = [], indices = [])
        nd, const_values = expr_to_nodedata(ex)
        adj = adjmat(nd)
        linearity = classify_linearity(nd, adj, [])
        @test linearity[1] == testval
        idxset = Coloring.IndexedSet(100)
        edgelist = compute_hessian_sparsity(
            nd,
            adj,
            linearity,
            idxset,
            Array{Set{Tuple{Int,Int}}}(undef, 0),
            Array{Vector{Int}}(undef, 0),
        )
        if linearity[1] != NONLINEAR
            @test length(edgelist) == 0
        elseif length(IJ) > 0
            @test IJ == edgelist
        end
        if length(indices) > 0
            ix = sort(collect(compute_gradient_sparsity(nd)))
            compute_gradient_sparsity!(idxset, nd)
            ix2 = sort(collect(idxset))
            empty!(idxset)
            @test ix == indices
            @test ix == ix2
        end
    end

    test_linearity(
        :(sin(x[1]^2) + cos(x[2] * 4) - 2.0),
        NONLINEAR,
        Set([(2, 2), (1, 1)]),
        [1, 2],
    )
    test_linearity(:(3 * 4 * (x[1] + x[2])), LINEAR)
    test_linearity(
        :(x[3] * x[2]),
        NONLINEAR,
        Set([(3, 2), (3, 3), (2, 2)]),
        [2, 3],
    )
    test_linearity(:(3 + 4), CONSTANT)
    test_linearity(:(sin(3) + x[1]), LINEAR)
    test_linearity(
        :(cos(x[3]) * sin(3) + x[1]),
        NONLINEAR,
        Set([(3, 3)]),
        [1, 3],
    )
    test_linearity(:(x[1] - 3x[2]), LINEAR)
    test_linearity(:(-x[1]), LINEAR)
    test_linearity(:(+x[1]), LINEAR)
    test_linearity(:(x[1]^x[2]), NONLINEAR, Set([(2, 2), (1, 1), (2, 1)]))
    test_linearity(:(x[1] / 3 + x[2]), LINEAR)
    test_linearity(
        :(3 / (x[1] * x[2])),
        NONLINEAR,
        Set([(2, 2), (1, 1), (2, 1)]),
    )
    test_linearity(:(1 / (x[1] + 3)), NONLINEAR)
    test_linearity(:(ifelse(x[1] <= 1, x[1], x[2])), PIECEWISE_LINEAR, Set([]))
    test_linearity(:(ifelse(x[1] <= 1, x[1]^2, x[2])), NONLINEAR, Set([(1, 1)]))
    test_linearity(:(ifelse(1 <= 1, 2, 3)), CONSTANT)
    test_linearity(:(1 / ifelse(x[1] < 1, x[1], 0)), NONLINEAR, Set([(1, 1)]))

    function test_registered_function_linearity()
        # See issue https://github.com/jump-dev/JuMP.jl/issues/2340
        foo(x, y) = x + y
        r = JuMP._Derivatives.UserOperatorRegistry()
        register_multivariate_operator!(
            r,
            :foo,
            JuMP._UserFunctionEvaluator(2, foo),
        )
        nd, _ = expr_to_nodedata(:(foo(x[1], x[1]^2)), r)
        linearity = classify_linearity(nd, adjmat(nd), Any[])
        @test linearity == [NONLINEAR, LINEAR, NONLINEAR, LINEAR, CONSTANT]
    end
    test_registered_function_linearity()

    function MOI.eval_objective(::ΦEvaluator, x)
        @assert length(x) == 2
        return (1 / 3) * x[2]^3 - 2x[1]^2
    end
    function MOI.eval_objective_gradient(::ΦEvaluator, grad, x)
        grad[1] = -4x[1]
        return grad[2] = x[2]^2
    end
    r = _Derivatives.UserOperatorRegistry()
    register_multivariate_operator!(r, :Φ, ΦEvaluator())
    register_univariate_operator!(r, :c, cos, x -> -sin(x), x -> -cos(x))
    Φ(x, y) = MOI.eval_objective(ΦEvaluator(), [x, y])
    ex = :(Φ(x[2], x[1] - 1) * c(x[3]))
    nd, const_values = expr_to_nodedata(ex, r)
    @test _Derivatives.has_user_multivariate_operators(nd)
    adj = adjmat(nd)
    storage = zeros(length(nd))
    partials_storage = zeros(length(nd))
    reverse_storage = zeros(length(nd))
    x = [2.0, 3.0, 4.0]
    fval = forward_eval(
        storage,
        partials_storage,
        nd,
        adj,
        const_values,
        [],
        x,
        [],
        zeros(2),
        zeros(2),
        r,
    )
    true_val = Φ(x[2], x[1] - 1) * cos(x[3])
    @test fval ≈ true_val
    grad = zeros(3)
    reverse_eval(reverse_storage, partials_storage, nd, adj)
    reverse_extract(grad, reverse_storage, nd, adj, [], 1.0)

    true_grad = [
        cos(x[3]) * (x[1] - 1)^2,
        -4cos(x[3]) * x[2],
        -sin(x[3]) * Φ(x[2], x[1] - 1),
    ]
    @test grad ≈ true_grad

    # dual forward test
    function dualforward(ex, x; ignore_nan = false)
        nd, const_values = expr_to_nodedata(ex)
        adj = adjmat(nd)
        forward_storage = zeros(length(nd))
        partials_storage = zeros(length(nd))
        reverse_storage = zeros(length(nd))

        fval = forward_eval(
            forward_storage,
            partials_storage,
            nd,
            adj,
            const_values,
            [],
            x,
            [],
            [],
            [],
            NO_USER_OPS,
        )
        grad = zeros(length(x))
        reverse_eval(reverse_storage, partials_storage, nd, adj)
        reverse_extract(grad, reverse_storage, nd, adj, [], 1.0)

        zero_ϵ = zero(ForwardDiff.Partials{1,Float64})
        forward_storage_ϵ = fill(zero_ϵ, length(nd))
        partials_storage_ϵ = fill(zero_ϵ, length(nd))
        x_values_ϵ = fill(ForwardDiff.Partials((1.0,)), length(x))
        reverse_storage_ϵ = fill(zero_ϵ, length(nd))
        output_ϵ = fill(zero_ϵ, length(x))
        fval_ϵ = forward_eval_ϵ(
            forward_storage,
            forward_storage_ϵ,
            partials_storage,
            partials_storage_ϵ,
            nd,
            adj,
            x_values_ϵ,
            [],
            NO_USER_OPS,
        )
        reverse_eval_ϵ(
            output_ϵ,
            reverse_storage,
            reverse_storage_ϵ,
            partials_storage,
            partials_storage_ϵ,
            nd,
            adj,
            [],
            [],
            2.0,
            zero_ϵ,
        )
        @test fval_ϵ[1] ≈ dot(grad, ones(length(x)))

        # compare with running dual numbers
        _epsilon(x::ForwardDiff.Dual{Nothing,Float64,1}) = x.partials[1]

        forward_dual_storage =
            zeros(ForwardDiff.Dual{Nothing,Float64,1}, length(nd))
        partials_dual_storage =
            zeros(ForwardDiff.Dual{Nothing,Float64,1}, length(nd))
        output_dual_storage =
            zeros(ForwardDiff.Dual{Nothing,Float64,1}, length(x))
        reverse_dual_storage =
            zeros(ForwardDiff.Dual{Nothing,Float64,1}, length(nd))

        x_dual = [ForwardDiff.Dual(x[i], 1.0) for i in 1:length(x)]
        fval = forward_eval(
            forward_dual_storage,
            partials_dual_storage,
            nd,
            adj,
            const_values,
            [],
            x_dual,
            [],
            [],
            [],
            NO_USER_OPS,
        )
        reverse_eval(reverse_dual_storage, partials_dual_storage, nd, adj)
        reverse_extract(
            output_dual_storage,
            reverse_dual_storage,
            nd,
            adj,
            [],
            ForwardDiff.Dual(2.0, 0.0),
        )
        for k in 1:length(nd)
            @test _epsilon(forward_dual_storage[k]) ≈ forward_storage_ϵ[k][1]

            if !(isnan(_epsilon(partials_dual_storage[k])) && ignore_nan)
                @test _epsilon(partials_dual_storage[k]) ≈
                      partials_storage_ϵ[k][1]
            else
                @test !isnan(forward_storage_ϵ[k][1])
            end
            if !(isnan(_epsilon(reverse_dual_storage[k])) && ignore_nan)
                @test _epsilon(reverse_dual_storage[k]) ≈
                      reverse_storage_ϵ[k][1] / 2
            else
                @test !isnan(reverse_storage_ϵ[k][1])
            end
        end
        for k in 1:length(x)
            if !(isnan(_epsilon(output_dual_storage[k])) && ignore_nan)
                @test _epsilon(output_dual_storage[k]) ≈ output_ϵ[k][1]
            else
                @test !isnan(output_ϵ[k][1])
            end
        end
    end

    dualforward(:(sin(x[1]^2) + cos(x[2] * 4) / 5 - 2.0), [1.0, 2.0])
    dualforward(:(sin(x[1]^x[2]) + cos(x[2] * 4) / 5 - 2.0), [1.0, 2.0])
    dualforward(:(sin(x[1]^3) + cos(x[1] * x[2] * 4) / 5 - 2.0), [1.0, 0.0])
    dualforward(
        :(x[1] * x[2]),
        [3.427139283036299e-206, 1.0],
        ignore_nan = true,
    )
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
## 调包
using CurveFit
using CSV
using Interpolations
using PyCall
using PyPlot
using DataFrames

## 调 python 包
@pyimport numpy as np
pygui(true)

## 原始数据
"""
原始数据输入
(进料组成, 塔顶产品组成, 塔釜组成) 的质量分数 (ω_f, ω_d, ω_w)
进料温度 t_f
回流流量 L
产品流量 D
"""

# (进料组成, 塔顶产品组成, 塔釜组成) 的质量分数
ω_f = 16.2419e-2
ω_d = 89.8466e-2
ω_w = 6.0228e-2

# 进料温度 / C
t_f = 33.9

# 回流流量 / (L * h ^ -1)
L = 3.6
# 产品流量 / (L * h ^ -1)
D = 1.5

# 质量分数 -> 摩尔分数 换算函数
f(ω_i) = (ω_i / 46.07) / (ω_i / 46.07 + (1 - ω_i) / 18.02)

# (进料组成, 塔顶产品组成, 塔釜组成) 摩尔分数
x_f = f(ω_f)
x_d = f(ω_d)
x_w = f(ω_w)

## 对角线
x = y = np.linspace(0, 1)

## 平衡线
# EtOH-H2O 相平衡数据
EtOH = CSV.read("EtOH-H2O.csv", DataFrame)
# EtOH-H2O 相平衡组成插值
interp_y = Interpolations.LinearInterpolation(EtOH.x, EtOH.y)
# 相平衡线
yy = interp_y(x)

## 精馏段操作线
R = L / D
f_r(x) = (R * x + x_d) / (R + 1)

## q 线
# 泡点温度插值
interp_t = Interpolations.LinearInterpolation(EtOH.x, EtOH.t)
# 泡点温度
t_s = interp_t(x_f)

# 乙醇摩尔热容
c_pm_EtOH(t) = (1.56 * 10e-2 * t + 2.012) * 46.07
# 水摩尔热容
c_pm_H2O(t) = (2.143 * 10e-4 * t + 4.198) * 18.02
# 溶液热容 / (kJ * kmol ^ -1 * C ^ -1)
t_m = (t_f + t_s) / 2
c_pm = x_f * c_pm_EtOH(t_m) + (1 - x_f) * c_pm_H2O(t_m)

# 乙醇摩尔气化潜热
r_m_EtOH(t) = 113 * (243 - t) ^ 0.4218 * 46.07
# 水的摩尔气化潜热
r_m_H2O(t) = 445.6 * (374 - t) ^ 0.3003 * 18.02
# 溶液气化潜热 / (kJ * kmol ^ -1)
r_m = x_f * r_m_EtOH(t_f) + (1 - x_f) * r_m_EtOH(t_f)

# q 值
q = 1 + c_pm * (t_s - t_f) / r_m

# q 线方程
f_q(x) = (q * x - x_f) / (q - 1)

## 交点
# 精馏段操作线和对角线交点
point_a = [x_d, x_d]
# 提馏段操作线和对角线的交点
point_b = [x_w, x_w]
# 精馏段操作线和 y 轴的交点
point_c = [0, f_r(0)]
# q 线和精馏段操作线的交点
point_d = [((R + 1) * x_f + (q - 1) * x_d) / (R + q), (R * x_f + q * x_d) / (R + q)]
# q 线和对角线的交点
point_f = [x_f, x_f]

## 提馏段操作线方程
p_d = CurveFit.linear_fit([point_b[1], point_d[1]], [point_b[2], point_d[2]])
f_d(x) = p_d[1] + p_d[2] * x

## 绘图
# 对角线
PyPlot.plot(x, y, "-")
# 相平衡线
PyPlot.plot(x, yy, "-")
# 精馏段操作线 a-c
PyPlot.plot([point_a[1], point_c[1]], [point_a[2], point_c[2]], "-")
# q 线 d-f
PyPlot.plot([point_d[1], point_f[1]], [point_d[2], point_f[2]], "-")
# 提馏段操作线 b-d
PyPlot.plot([point_b[1], point_d[1]], [point_b[2], point_d[2]], "-")
# x_w 垂直线
PyPlot.plot([x_w, x_w], [0, x_w], "b-")
# x_d 垂直线
PyPlot.plot([x_d, x_d], [0, x_d], "b-")
# x_f 垂直线
PyPlot.plot([x_f, x_f], [0, x_f], "b-")

# 绘图设置
PyPlot.xlabel("\$ x \$")
PyPlot.ylabel("\$ y \$")
PyPlot.xlim([0, 1])
PyPlot.ylim([0, 1])

## 塔板数计算
# 迭代点
point_1 = [x_d, x_d]
point_2 = zeros(2)
# 塔板数
num = 0
# 迭代记数
n = 1
# 迭代计算 & 绘图
println("开始迭代")
while true
    global point_1, point_2, num, n
    if point_1[1] < x_w
        break
    end
    println("迭代第 $(n) 次")
    n += 1
    for x_i in np.linspace(0, 1, 10000000)
        if abs(interp_y(x_i) - point_1[2]) < 0.00001 * point_1[2]
            point_temp = [x_i, point_1[2]]
            point_2 = zeros(2)
            if point_d[1] < x_i
                # 在精馏段操作线上
                point_2 = [x_i, f_r(x_i)]
            elseif point_b[1] < x_i < point_d[1]
                # 在提馏段操作线上
                point_2 = [x_i, f_d(x_i)]
            elseif x_i < point_b[1]
                # 在对角线上
                point_2 = [x_i, x_i]
            end
            num += 1
            PyPlot.plot([point_1[1], point_temp[1]], [point_1[2], point_temp[2]], "b-")
            PyPlot.plot([point_2[1], point_temp[1]], [point_2[2], point_temp[2]], "b-")
            # 迭代后新点
            point_1 = point_2
            point_2 = zeros(2)
            break
        end
    end
end

# 理论塔板数结果
println("理论塔板数 = $(num)")

max_stars_repo_path: Experiment 07/partical_reflux.jl

max_issues_repo_path: Experiment 07/partical_reflux.jl

max_forks_repo_path: Experiment 07/partical_reflux.jl


text:
# Functions for large C4 correlation job
export correlate_pair, autocorrelate, diag_chunks, offdiag_chunks, preprocess2, get_blocks,
        stack_auto, stack_corr, correlate_big, stack_all

# correlation helper functions
function correlate_pair(src::FFTData, rec::FFTData, pref::String="CORR", params::Dict=params)
    try
        C = correlate(src, rec, params["maxlag"])
        cc_medianmute!(C, 10.) # remove correlation windows with high noise
        stack!(C)
        name = join(split(name_corr(C), ".")[1:2],".")
        save_named_corr(C,"$(params["rootdir"])/$pref/$name/$(C.comp)")
    catch e
        println(e)
    end
end

function autocorrelate(station::String, list::Array{String,1}=fft_list_100, params::Dict=params)
    """ Filter channels by station, load and correlate """
    channels = filter(x -> occursin(station, x), list)
    ffts = map(x -> load_fft(x, string(x[end-7:end-5])), channels)
    pairs = vec([collect(x) for x in Iterators.product(1:length(ffts),1:length(ffts))])
    map(pair -> correlate_pair(ffts[pair[1]], ffts[pair[2]], "AUTOCORR", params), pairs)
end

function diag_chunks(chunk::Array{String,1}, prefix::String = "CORR", filt_dist::Bool = true, params::Dict=params)
    ffts  = map(x -> load_fft(x, string(x[end-7:end-5])), chunk)
    pairs = vec([collect(x) for x in Iterators.product(1:length(ffts),1:length(ffts))])
    # maybe remove this filter so we get all components
    #filter!(x -> (x[1] < x[2]) || (x[1]>x[2]), pairs) # filter upper diagonal
    if filt_dist; filter!(x -> get_dist(ffts[x[1]], ffts[x[2]]) <= 300., pairs); end # filter distances
    map(pair -> correlate_pair(ffts[pair[1]], ffts[pair[2]], prefix, params), pairs)
end

function offdiag_chunks(chunkpair::Array{Array{String, 1}}, prefix::String = "CORR", filt_dist::Bool = true,
                        params::Dict=params)
    sources   = map(x -> load_fft(x, string(x[end-7:end-5])), chunkpair[1]) # load
    receivers = map(x -> load_fft(x, string(x[end-7:end-5])), chunkpair[2])
    pairs = vec([collect(x) for x in Iterators.product(1:length(sources),1:length(receivers))]) # all pairs src -> rec
    if filt_dist; filter!(x -> get_dist(sources[x[1]].loc, receivers[x[2]].loc) <= 300., pairs); end # filt distances
    map(pair -> correlate_pair(sources[pair[1]], receivers[pair[2]], prefix, params), pairs)
end

function get_blocks(paths::Array{String,1}, params::Dict=params)
    # Chunk raw ffts into blocks. 30 is best size for IO, but distribute across cores if more cores available
    chunks = collect(Iterators.partition(paths, 30))#convert(Int64, minimum([30, ceil(length(paths)/params["num_procs"])]))))
    off_chunks = vec([collect(x) for x in Iterators.product(1:length(chunks),1:length(chunks))])
    filter!(chunk_pair -> chunk_pair[1] < chunk_pair[2], off_chunks)
    # map chunk indices to filenames
    off_chunk_names = map(chunk -> [collect(chunks[chunk[1]]), collect(chunks[chunk[2]])], off_chunks)
    return chunks, off_chunk_names
end

# preprocessing
function preprocess2(file::String, accelerometer::Bool=false, params::Dict=params, path::String=path)
    """
    Load raw seisdata file and process, saving to fft
    """
    try
        data = SeisData()
        if occursin(".ms", file)
            read_data!(data, "mseed", file)
        else # data is SeisIO native from NCEDC/IRIS
            read_data!(data, "seisio", file)
        end
        gaps = size(data.t[1])[1] # N-2 gaps (eg gaps = 12 tests for 10 gaps)
        pts = size(data.x[1])[1]
        fs_temp = data.fs[1]
        windows = u2d.(SeisIO.t_win(data.t[1], data.fs[1]) .* 1e-6)
        startend = windows[:,2] .- windows[:,1] .> Second(params["cc_len"])
        if gaps < 25 && any(startend) # If few gaps and sufficient (2/3) data present, data is clean
            try
                add_location(data, params["all_stations"])
                for samp_rate in params["samp_rates"]
                    S = SeisData()
                    try
                        S = process_raw(data, samp_rate)
                    catch e# trying to sample non 100Hz data at 100 Hz - skip resampling
                        println(e)
                        S = process_raw(data, data.fs[1])
                    end
                    R = RawData(S,params["cc_len"],params["cc_step"])
                    SeisNoise.detrend!(R)
                    bandpass!(R,params["freqmin"],params["freqmax"],zerophase=true)
                    # SMOOTHING LINE GOES HERE
                    SeisNoise.taper!(R)
                    FFT = nothing
                    if accelerometer # accelerometer - so we integrate
                        FFT = rfft_raw(R,1)
                    else # then its a seismometer
                        FFT = compute_fft(R)
                    end
                    coherence!(FFT,params["half_win"], params["water_level"])
                    try # save fft
                        root_fft = "ffts/$path/$(Int(samp_rate))/"
                        save_fft(FFT, joinpath(params["rootdir"], root_fft))
                    catch e
                        println(e)
                    end
                    #println("Successfully processed $(data.id[1])")
                end
                return data.id
            catch e
                println(e)
            end
        end
    catch e
        println(e)
    end
end

# stacking codes
function stack_auto(name::String, startdate::Date=startdate)
    month = Dates.monthname(startdate) # get month for filename
    yr = Dates.year(startdate)
    # stack autocorrelations

    CORROUT = expanduser(joinpath(params["rootdir"],"autocorrelations/$name/"))
    if !isdir(CORROUT)
        mkpath(CORROUT)
    end
    filename = joinpath(CORROUT,"$(yr)_$(month)_$name.h5") # get output filename

    # Get list of files to save
    autocorr_list = glob("AUTOCORR/$name*/*/*.jld2", params["rootdir"])
    components = ["EE", "EN", "EZ", "NE", "NN", "NZ", "ZE", "ZN", "ZZ"]

    C = load_corr(autocorr_list[1], convert(String, split(autocorr_list[1],"/")[end-1])) # sample autocorr for meta
    h5open(filename, "cw") do file
        if !haskey(read(file), "meta")
            write(file, "meta/corr_type", C.corr_type)
            write(file, "meta/cc_len", C.cc_len)
            write(file, "meta/cc_step", C.cc_step)
            write(file, "meta/whitened", C.whitened)
            write(file, "meta/time_norm", C.time_norm)
            write(file, "meta/notes", C.notes)
            write(file, "meta/maxlag", C.maxlag)
            write(file, "meta/starttime", Dates.format(u2d(C.t[1]), "yyyy-mm-dd HH:MM:SS"))
            write(file, "meta/samp_freq", C.fs)
            write(file, "meta/lat", C.loc.lat)
            write(file, "meta/lon", C.loc.lon)
            write(file, "meta/el", C.loc.el)
            write(file, "meta/dep", C.loc.dep)
            write(file, "meta/az", C.loc.az)
            write(file, "meta/inc", C.loc.inc)
        end
        # stack per month
        for comp in components
            #comp_files = filter(f -> convert(String, split(f, "/")[end-1]) == comp, autocorr_list)
            comp_files = glob("AUTOCORR/$name*/$comp/*.jld2", params["rootdir"]) # uni-component files
            autocorrs = [load_corr(f, comp) for f in comp_files]

            autocorr_mean = SeisNoise.stack(sum(autocorrs), allstack=true, stacktype=mean)
            autocorr_pws = SeisNoise.stack(sum(autocorrs), allstack=true, stacktype=pws)
            autocorr_robust = SeisNoise.stack(sum(autocorrs), allstack=true, stacktype=robuststack)
            # save into file
            write(file, "$rec/$comp/linear", autocorr_mean.corr[:])
            write(file, "$rec/$comp/pws", autocorr_pws.corr[:])
            write(file, "$rec/$comp/robust", autocorr_robust.corr[:])
        end

        # Following code saves each correlation individually
        # for (ind, name) in enumerate(autocorr_list)
        #     # load file
        #     comp = string(split(name, "/")[end-1]) # get component
        #     Cl = load_corr(name, comp)
        #     # write
        #     write(file, "$comp/$(Cl.id)", Cl.corr[:])
        # end
    end
end

function stack_corr(name::String, startdate::Date=startdate, prefix::String = "CORR_20HZ")
    month = Dates.monthname(startdate) # get month for filename
    yr = Dates.year(startdate)
    # stack autocorrelations

    CORROUT = expanduser(joinpath(params["rootdir"], "correlations/$name/"))
    if !isdir(CORROUT)
        mkpath(CORROUT)
    end
    filename = joinpath(CORROUT,"$(yr)_$(month)_$name.h5") # get output filename
    # r = r"(?<=$(name)).*"
    components = ["EE","EN","EZ", "NE", "NN","NZ", "ZE", "ZN", "ZZ"]

    receivers = glob("$prefix/$name/*/*", params["rootdir"])
    receivers = Set([join(split(x, ".")[end-2:end-1], ".") for x in receivers]) # just get reciever names
    println(receivers)
    corr_list = glob("$prefix/$name*/*/*", params["rootdir"])
    C = load_corr(corr_list[1], convert(String, split(corr_list[1],"/")[end-1]))
    source_loc = GeoLoc(lat = C.loc.lat, lon = C.loc.lon, el = C.loc.el)
    h5open(filename, "cw") do file
        if !haskey(read(file), "meta")
            write(file, "meta/corr_type", C.corr_type)
            write(file, "meta/cc_len", C.cc_len)
            write(file, "meta/cc_step", C.cc_step)
            write(file, "meta/whitened", C.whitened)
            write(file, "meta/time_norm", C.time_norm)
            write(file, "meta/notes", C.notes)
            write(file, "meta/maxlag", C.maxlag)
            write(file, "meta/starttime", Dates.format(u2d(C.t[1]), "yyyy-mm-dd HH:MM:SS"))
            write(file, "meta/samp_freq", C.fs)
            write(file, "meta/lat", C.loc.lat)
            write(file, "meta/lon", C.loc.lon)
            write(file, "meta/el", C.loc.el)
            write(file, "meta/dep", C.loc.dep)
            write(file, "meta/az", C.loc.az)
            write(file, "meta/inc", C.loc.inc)
        end
        for rec in receivers
            try
                # sample_r = glob("CORR/$tf*$rec*/ZZ/*")[1]
                # Cr = load_corr(sample_r, "ZZ")
                rec_station = split(rec, ".")[2]
                #rec_loc = LLE_geo(rec_station, all_stations)
                rec_loc = load_fft(glob("ffts/*/*/$rec*BHZ*", params["rootdir"])[1], "BHZ").loc # get location from fft - perhaps a faster way to do this.
                if !haskey(read(file), "$rec/meta")
                    write(file, "$rec/meta/lon", rec_loc.lon)
                    write(file, "$rec/meta/lat", rec_loc.lat)
                    write(file, "$rec/meta/el", rec_loc.el)
                    write(file, "$rec/meta/dist", get_dist(source_loc, rec_loc))
                    write(file, "$rec/meta/azi", get_azi(source_loc, rec_loc))
                    write(file, "$rec/meta/baz", get_baz(source_loc, rec_loc))
                end
                for comp in components
                    # load correlations for this receiver by component
                    files = glob("$prefix/$name/$comp/*$rec*", params["rootdir"])
                    corrs = [load_corr(f, comp) for f in files]

                    # implement various stacktypes
                    corr_mean = SeisNoise.stack(sum(corrs), allstack=true, stacktype=mean)
                    corr_pws = SeisNoise.stack(sum(corrs), allstack=true, stacktype=pws)
                    corr_robust = SeisNoise.stack(sum(corrs), allstack=true, stacktype=robuststack)
                    # save into file
                    write(file, "$rec/$comp/linear", corr_mean.corr[:])
                    write(file, "$rec/$comp/pws", corr_pws.corr[:])
                    write(file, "$rec/$comp/robust", corr_robust.corr[:])
                end
            catch e
                println(e)
            end
        end
    end
end

function stack_all(params::Dict=params)
    autocorr_names = [join(split(split(auto, "/")[end], ".")[1:2],".") for auto in glob("AUTOCORR/*", params["rootdir"])]
    @eval @everywhere autocorr_names = $autocorr_names
    Sauto = @elapsed robust_pmap(x -> stack_auto(x, startdate), autocorr_names)

    corr_names = [join(split(split(auto, "/")[end], ".")[1:2],".") for auto in glob("CORR_20HZ/*", params["rootdir"])]
    @eval @everywhere corr_names = $corr_names
    S20 = @elapsed robust_pmap(x -> stack_corr(x, startdate), corr_names)

    corr_names_lf = [join(split(split(auto, "/")[end], ".")[1:2],".") for auto in glob("CORR_1HZ/*", params["rootdir"])]
    @eval @everywhere corr_names_lf = $corr_names_lf
    S1 = @elapsed robust_pmap(x -> stack_corr(x, startdate), corr_names_lf)
    println("Stacking Completed for autocorrs and interstation cross correlations in $(Sauto+S20+S1) seconds")
end

# wrapper function for correlations
function correlate_big(dd::Date, startdate::Date = startdate, params::Dict = params)
    """ Wrapper Function: Computes autocorrelations for a specific day"""
    aws = params["aws"]
    yr = Dates.year(dd)
    path = join([Dates.year(dd),lpad(Dates.dayofyear(dd),3,"0")],"_") # Yeilds "YEAR_JDY"
    @eval @everywhere path = $path

    ###### WRAP THIS MESS #######
    # get filepaths for source stations - would be faster if they're available
    scedc_files = nothing
    if !isdir("root/scedc_path/$yr/"); mkpath("root/scedc_path/$yr/"); end
    try # most filelists are stored on seisbasin
        s3_get_file(aws, "seisbasin", "scedc_path/$yr/$path.csv", "root/scedc_path/$yr/$path.csv")
        scedc_files = DataFrame(CSV.File("root/scedc_path/$yr/$path.csv")).Path
    catch e # in case file not found/present we use SCEDC lookup
        println(e)
        scedc_files = get_scedc_files(dd, aws)
    end
    println(length(scedc_files), " is the length of scedc_files")
    #########################

    # filepaths for iris, ncedc, and node data
    iris_path  = S3Path("s3://seisbasin/iris_waveforms/$yr/$path/", config=aws)
    ncedc_path = S3Path("s3://seisbasin/ncedc_waveforms/$yr/$path/", config=aws)
    node_path  = S3Path("s3://seisbasin/continuous_waveforms/$yr/$path/", config=aws)
    iris_query, ncedc_query, node_query = convert.(String, readdir(iris_path)), convert.(String, readdir(ncedc_path)),
                                        convert.(String, readdir(node_path))
    filelist_basin = vcat(joinpath.("iris_waveforms/$yr/$path/", iris_query), joinpath.("ncedc_waveforms/$yr/$path/", ncedc_query),
                                        joinpath.("continuous_waveforms/$yr/$path/", node_query))
    println("There are $(length(filelist_basin)) node files and $(length(scedc_files)) SCEDC files available for $path.")

    # download scedc and seisbasin data
    ec2download(aws, "scedc-pds", scedc_files, joinpath(params["rootdir"],"data"))
    ec2download(aws, "seisbasin", filelist_basin, joinpath(params["rootdir"],"data"))
    println("Download complete!")

    # preprocess - broadbands and seismometers are processed separately
    allf = glob("data/*/$yr/$path/*", params["rootdir"])
    println(allf)
    println(readdir())
    println(readdir("$(params["rootdir"])/data"))
    accelerometers = filter(x -> isfile(x), joinpath.("/$(params["rootdir"])/data/continuous_waveforms/$yr/$path/", node_query))
    broadbands = setdiff(allf, accelerometers) # get the rest of the files

    @eval @everywhere accelerometers, broadbands = $accelerometers, $broadbands
    T_b = @elapsed robust_pmap(f -> preprocess2(f, false, params, path), broadbands)
    if length(accelerometers) != 0 # save time on pmap allocation overhead
        T_a = @elapsed robust_pmap(f -> preprocess2(f, true, params, path), accelerometers) # integrate accelerometers
        println("Preprocessing Completed in $(T_a+T_b) seconds.")
    else
        println("Preprocessing Completed in $T_b seconds.")
    end

    # autocorrelations
    fft_list_100 = glob("ffts/$path/100/*.jld2", params["rootdir"])
    println("FFT list 100 is $(length(fft_list_100)) stations long")
    fft_100_stations = unique([join(split(elt, ".")[1:2],".") for elt in fft_list_100]) # get stations to iterate over
    @eval @everywhere fft_100_stations = $fft_100_stations
    pmap(x -> autocorrelate(x, fft_list_100, params), fft_100_stations)

    fft_paths_20 = sort(glob("ffts/$path/20/*", params["rootdir"])) # alphabetic sort for all stations so that diagonal files are made correctly
    fft_paths_1 = sort(glob("ffts/$path/1/*", params["rootdir"])) # alphabetic sort for all stations so that diagonal files are made correctly
    println("There are $(length(fft_paths_20)) fft datas for the 20 HZ correlations")
    # run 20 HZ correlations
    chunks_20HZ, off_chunk_names_20HZ = get_blocks(fft_paths_20, params);

    # correlate diagonal chunks
    @eval @everywhere chunks_20HZ = $chunks_20HZ
    T20D = @elapsed robust_pmap(chunk -> diag_chunks(convert(Array,chunk), "CORR_20HZ", true, params), chunks_20HZ)

    # correlate off-diagonal chunks
    @eval @everywhere off_chunk_names_20HZ = $off_chunk_names_20HZ
    T20O = @elapsed robust_pmap(chunk -> offdiag_chunks(chunk, "CORR_20HZ", true, params), off_chunk_names_20HZ) # run mega correlations

    # run low freq correlations
    chunks_1HZ, off_chunk_names_1HZ = get_blocks(fft_paths_1, params);

    # correlate diagonal chunks
    @eval @everywhere chunks_1HZ = $chunks_1HZ
    T1D = @elapsed robust_pmap(chunk -> diag_chunks(convert(Array,chunk), "CORR_1HZ", false, params), chunks_1HZ)

    # correlate off-diagonal chunks
    @eval @everywhere off_chunk_names_1HZ = $off_chunk_names_1HZ
    T1O = @elapsed robust_pmap(chunk -> offdiag_chunks(chunk, "CORR_1HZ", false, params), off_chunk_names_1HZ) # run mega correlations

    # Correlation summary
    println("All $(length(glob("CORR_*/*/*/$path*", params["rootdir"]))) Inter-station Correlations computed in $(T20D + T20O + T1D + T1O) seconds")
    try; rm("$(params["rootdir"])/data/continuous_waveforms/", recursive=true); catch e; println(e); end
end


text:
@testset "Derivatives" begin
    @testset "Normal finite-differences" begin
        B = FiniteDifferences(5,1.0)
        D = Derivative(axes(B,1))

        ∇ = applied(*, B', D, B)
        ∇² = applied(*, B', D', D, B)

        ∂ = materialize(∇)
        T = materialize(∇²)

        @test ∂ isa Tridiagonal
        @test T isa SymTridiagonal

        @test all(diag(∂) .== 0)
        @test all(diag(∂,1) .== 0.5)
        @test all(diag(∂,-1) .== -0.5)

        @test all(diag(T) .== -2)
        @test all(diag(T,1) .== 1)
        @test all(diag(T,-1) .== 1)
    end

    @testset "Restricted bases" begin
        N = 10
        ρ = 1.0

        R = StaggeredFiniteDifferences(N, ρ)
        r = axes(R, 1)
        D = Derivative(r)

        ∇ = R'D*R
        ∇² = R'D'D*R

        for (sela,selb) in Iterators.product([1:10, 3:6, 8:10, 5:10, 4:5], [1:10, 3:6, 8:10, 5:10, 4:5])
            # In case we have partially overlapping restrictions, the
            # bandwidth will not be the full three-point stencil, and if we
            # have no overlap at all, the bandwidth will be zero.
            overlap = sela ∩ selb
            expected_bandwidth = isempty(overlap) ? 0 : min(3, length(overlap)+1)

            Ra = R[:,sela]
            Rb = R[:,selb]

            ∂ = Ra'D*Rb
            @test ∂ == Matrix(∇)[sela,selb]
            if sela == selb
                @test bandwidths(∂) == (1,1)
            end
            @test ∂ isa BandedMatrix
            @test length(bandrange(∂)) == expected_bandwidth

            ∂² = Ra'D'*D*Rb
            @test ∂² == Matrix(∇²)[sela,selb]
            if sela == selb
                @test bandwidths(∂) == (1,1)
            end
            @test ∂² isa BandedMatrix
            @test length(bandrange(∂²)) == expected_bandwidth
        end
    end

    @testset "Implicit finite differences" begin
        N = 10
        ρ = 0.3
        Z = 1.0
        R = ImplicitFiniteDifferences(N, ρ)
        r = axes(R, 1)
        @testset "Singular origin=$(singular_origin)" for singular_origin = [true, false]
            D̃ = Derivative(r)
            D = singular_origin ? CoulombDerivative(D̃, Z, 0) : D̃

            λ,δβ₁ = if singular_origin
                √3 - 2, -Z*ρ*inv(12 - 10*Z*ρ)
            else
                0,0
            end

            ∂ = R'D*R
            @test ∂ isa CompactBases.ImplicitDerivative

            @test ∂.Δ isa Tridiagonal
            @test ∂.Δ[1,1] ≈ λ/2ρ
            @test all(iszero, ∂.Δ.d[2:end])
            @test all(∂.Δ.dl .≈ -1/2ρ)
            @test all(∂.Δ.du .≈ 1/2ρ)

            @test ∂.M isa SymTridiagonal
            @test ∂.M[1,1] ≈ (4+λ)/6
            @test all(∂.M.dv[2:end] .≈ 4/6)
            @test all(∂.M.ev .≈ 1/6)

            ∂² = R'D'D*R
            @test ∂ isa CompactBases.ImplicitDerivative

            @test ∂².Δ isa SymTridiagonal
            @test ∂².Δ[1,1] ≈ -2*(1+δβ₁)/ρ^2
            @test all(∂².Δ.dv[2:end] .≈ -2/ρ^2)
            @test all(∂².Δ.ev .≈ 1/ρ^2)

            @test ∂².M isa SymTridiagonal
            @test ∂².M[1,1] ≈ (-1/2)*(-(10-2δβ₁)/6)
            @test all(∂².M.dv[2:end] .≈ (-1/2)*(-10/6))
            @test all(∂².M.ev .≈ (-1/2)*(-1/6))
        end

        @testset "Restricted bases" begin
            R̃ = R[:, 1:5]
            D = Derivative(axes(R̃, 1))

            @test_throws DimensionMismatch R̃'D*R
            @test_throws DimensionMismatch R'D*R̃
            @test_throws DimensionMismatch R̃'D'D*R
            @test_throws DimensionMismatch R'D'D*R̃

            ∂ = R̃'D*R̃
            @test ∂ isa CompactBases.ImplicitDerivative
            @test size(∂) == (5,5)

            ∂² = R̃'D'D*R̃
            @test ∂² isa CompactBases.ImplicitDerivative
            @test size(∂²) == (5,5)

            @warn "Need to test correctness of ImplicitDerivative in restricted bases"
        end
    end
end

@testset "Derivative accuracy" begin
    Ns = 2 .^ (5:20)

    d = 1.0
    f,g,h,a,b = derivative_test_functions(d)

    @testset "kind = $B" for (order,B) in [(2,FiniteDifferences),
                                           (4,ImplicitFiniteDifferences)]
        @info "$B derivative accuracy"
        hs,ϵg,ϵh,ϵh′,pg,ph,ph′ = compute_derivative_errors(Ns, f, g, h, 1) do N
            B(range(a, stop=b, length=N+2)[2:end-1])
        end

        @test isapprox(pg, order, atol=0.03) || pg > order
        @test isapprox(ph, order, atol=0.03) || ph > order
    end
    @testset "kind = StaggeredFiniteDifferences" begin
        dd = b-a
        @testset "$(label)" for (order,label,fun) = [(2.0, "Uniform", N -> (N,dd/(N+1))),
                                                     (2.0, "Non-uniform", N -> (dd*(1 .- cos.(π*range(0,stop=1,length=N+2)[2:end-1]))/2,))]
            @info "$(label) StaggeredFiniteDifferences derivative accuracy"
            hs,ϵg,ϵh,ϵh′,pg,ph,ph′ = compute_derivative_errors(Ns, x -> f(x-dd/2), x -> g(x-dd/2), x -> h(x-dd/2), 1) do N
                StaggeredFiniteDifferences(fun(N)...)
            end

            @test isapprox(pg, order, atol=0.03) || pg > order
            @test isapprox(ph, order, atol=0.03) || ph > order
        end
    end
end

@testset "Particle-in-a-box" begin
    @testset "Eigenvalues convergence rate" begin
        Ns = 2 .^ (7:14)
        nev = 3
        L = 1.0
        @testset "kind = $B" for (order,B) in [(2,FiniteDifferences),
                                               (4,ImplicitFiniteDifferences)]
            @info "$B particle-in-a-box eigenvalues convergence rate"
            ϵλ,slopes,elapsed =
                compute_diagonalization_errors(N -> B(N, L/(N+1)), test_particle_in_a_box,
                                               Ns, L, nev,
                                               verbosity=1)
            for (p,o) in zip(slopes, [order,order-1,order-1])
                @test isapprox(p, o, atol=0.1) || p > o
            end
        end
    end
end

# This tests the discretization of the Laplacian, especially near
# the origin where the potential term is singular.
@testset "Hydrogen bound states" begin
    @testset "Eigenvalues and eigenvectors" begin
        rₘₐₓ = 300
        ρ = 0.25
        N = ceil(Int, rₘₐₓ/ρ + 1/2)

        R = StaggeredFiniteDifferences(N, ρ)
        r = axes(R, 1)

        D = CoulombDerivative(Derivative(Base.axes(R,1)), 1, 0)
        Tm = (R'D'D*R) / -2
        @test Tm isa SymTridiagonal

        V = R'*QuasiDiagonal((r -> -inv(r)).(r))*R

        H = Tm + V

        ee = eigen(H)

        n = 1:10
        λₐ = -inv.(2n.^2)

        abs_error = ee.values[n] - λₐ
        rel_error = abs_error ./ abs.(1e-10 .+ abs.(λₐ))

        @test abs_error[1] < 3e-5
        @test all(rel_error .< 1e-3)

        r̃ = locs(R)
        # Table 2.2 Foot (2005)
        Rₐ₀ = [2exp.(-r̃),
               -(1/2)^1.5 * 2 * (1 .- r̃/2).*exp.(-r̃/2), # Why the minus?
               (1/3)^1.5 * 2 * (1 .- 2r̃/3 .+ 2/3*(r̃/3).^2).*exp.(-r̃/3)]
        expected_errors = [1e-3,1e-3,2e-3]

        for i = 1:3
            v = ee.vectors[:,i]
            N = norm(v)*√ρ
            # The sign from the diagonalization is arbitrary; make max lobe positive
            N *= sign(v[argmax(abs.(v))])
            abs_error = v/N .- r̃.*Rₐ₀[i]
            @test norm(abs_error)/abs(1e-10+norm(r̃.*Rₐ₀[i])) < expected_errors[i]
        end
    end

    @testset "Eigenvalues convergence rate" begin
        Ns = 2 .^ (7:14)
        nev = 3
        rₘₐₓ = 100.0
        Z = 1.0
        ℓ = 0

        @testset "kind = $B" for (order,B) in [(2,StaggeredFiniteDifferences),
                                               (4,ImplicitFiniteDifferences)]
            @info "$B hydrogen eigenvalues convergence rate"
            ϵλ,slopes,elapsed = compute_diagonalization_errors(test_singular_scheme, Ns, Z, ℓ, nev, verbosity=1) do N
                ρ = B == StaggeredFiniteDifferences ? rₘₐₓ/(N+0.5) : rₘₐₓ/(N+1)
                B(N, ρ)
            end
            for p in slopes
                @test isapprox(p, order, atol=0.04) || p > order
            end
        end
    end
end


text:
## binning algorithms

struct EqualSize
    "Number of bins."
    n::Int

    @doc """
        EqualSize(; n::Int=10)

    Create binning algorithm with `n` bins of the same size.
    """
    function EqualSize(; n::Int=10)
        n > 0 || error("number of bins must be positive")
        return new(n)
    end
end

struct EqualMass
    "Number of bins."
    n::Int

    @doc """
        EqualMass(; n::Int=10)

    Create binning algorithm with `n` bins of (approximately) equal mass.
    """
    function EqualMass(; n::Int=10)
        n > 0 || error("number of bins must be positive")
        return new(n)
    end
end

function bins(::Any, ::Any, algorithm::EqualSize)
    return StatsBase.Histogram(range(0, 1; length=algorithm.n + 1))
end
function bins(probabilities, ::Any, algorithm::EqualMass)
    return StatsBase.Histogram(StatsBase.nquantile(probabilities, algorithm.n))
end

function binindex(bins::StatsBase.Histogram{<:Real,1}, x::Real)
    index = StatsBase.binindex(bins, x)
    nbins = length(bins.weights)
    edges = first(bins.edges)
    return if index == 0 && x == first(edges)
        1
    elseif index == nbins + 1 && x == last(edges)
        nbins
    else
        index
    end
end

function max_nonempty_bins(bins::StatsBase.Histogram{<:Real,1}, n::Int)
    return min(length(bins.weights), n)
end


text:
################################################################################
#
#   module LatticePhysics_LatticeConstruction
#   -> LatticePhysics_Base
#   -> LinearAlgebra
#
#   --> NAIVE CONSTRUCTION OF LATTICES
#           - periodic pattern of unitcells
#           - semi-periodic pattern of unitcells
#           - open pattern of unitcells
#           - by bond distance to an origin site
#           - in a shape around an origin site
#
#   --> ADDING NNN TO LATTICE
#
#   --> LATTICE CONSTRUCTION ALGORITHM FROM SPIN SPIRALS PAPER
#
################################################################################

# module start
module LatPhysLatticeConstruction

# uses the Base package
using LatPhysBase

# FUNCTIONS TO DEFINE INDICES IN ARRAYS
include("index_functions.jl")

# HELPER FUNCTIONS
include("helper_functions.jl")

# CONSTRUCTION OF LATTICES FROM UNITCELLS
# BY PUTTING TOGETHER UNITCELLS
# (currently: only periodic and open patterns)
include("construction_uc_patterns_periodic.jl")
include("construction_uc_patterns_open.jl")

# CONSTRUCTION OF LATTICES FROM UNITCELLS
# BY SPREADING OUTWARD OVER BONDS
include("construction_by_bond_distance.jl")

# CONSTRUCTION OF LATTICES FROM UNITCELLS
# BY SPREADING OUTWARD INSIDE A GIVEN SHAPE
include("construction_in_shape.jl")

# FUNCTIONS TO CONSTRUCT NEIGHBORS / BONDS

# Construction of NNN neighbors by bond distance
include("neighbors_NNN_by_bonds.jl")

# module end
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

@stochastic_model integer begin
    @stage 1 begin
        @decision(integer, x₁, Bin)
        @decision(integer, x₂, Bin)
        @objective(integer, Max, 6*x₁ + 4*x₂)
    end
    @stage 2 begin
        @uncertain T₁ T₂ h
        @recourse(integer, y₁, Bin)
        @recourse(integer, y₂, Bin)
        @recourse(integer, y₃, Bin)
        @recourse(integer, y₄, Bin)
        @objective(integer, Max, 3*y₁ + 7*y₂ + 9*y₃ + 6*y₄)
        @constraint(integer, 2*y₁ + 4*y₂ + 5*y₃ + 3*y₄ <= h - T₁*x₁ - T₂*x₂)
    end
end

ξ₁ = @scenario T₁ = 2. T₂ = 4. h = 10. probability = 0.25
ξ₂ = @scenario T₁ = 4. T₂ = 3. h = 15. probability = 0.75


text:
using SuiteSparseGraphBLAS
function triangle_centrality1(A)
    T = mul(A, A', (+, *), mask=A)
    y = reduce(+, T, dims=2)
    k = reduce(+, y)
    return (3 * mul(A, y) - 2 * mul(T, y) .+ y) ./ k
end

function triangle_centrality2(A)
    T = mul(A, A', (+, pair), mask=A, Descriptor(structural_mask=true))
    y = reduce(+, T, dims=2)
    k = reduce(+, y)
    return map(*, eadd(eadd(3.0 * mul(A, y), 2.0 * mul(T, y, (+, second)), -), y), 1/k)
end


text:
# onehot array
struct OneHotArray{K, N, var"N+1", A<:AbstractArray{OneHot{K}, N}} <: AbstractOneHotArray{var"N+1"}
    onehots::A
    function OneHotArray{K, N, var"N+1", A}(onehots::A) where {K, N, var"N+1", A<:AbstractArray{OneHot{K}, N}}
        @assert K isa UInt32
        @assert N+1 == var"N+1"
        return new{K, N, var"N+1", A}(onehots)
    end
end

OneHotArray(onehots::A) where {K, A<:AbstractArray{OneHot{K}}} = OneHotArray{K, ndims(onehots), ndims(onehots)+1, A}(onehots)
OneHotArray{K}(indices::A) where {K, A<:AbstractArray{<:Union{Int32, UInt32}}} = OneHotArray(reinterpret(OneHot(K), indices))
OneHotArray{K}(xs) where K = OneHotArray(OneHot(K).(xs))
OneHotArray(k, xs) where {A<:AbstractArray{<:Integer}} = OneHotArray{UInt32(k)}(xs)

OneHotArray(k, o::OneHot) = OneHotArray([OneHot(k, o)])
OneHotArray(o::OneHot) = OneHotArray([o])

OneHotArray{K, N}(xs::AbstractArray{T, N}) where {K, N, T} = OneHotArray{K, N, N+1}(xs)
function OneHotArray{K, N, var"N+1"}(xs::AbstractArray{T, N}) where {K, N, var"N+1", T}
    @assert N+1 == var"N+1"
    return OneHotArray(K, xs)
end

const OneHotVector{K} = OneHot{K}
const OneHotMatrix{K} = OneHotArray{K, 1}

onehotsize(::OneHotArray{K}) where K = Int(K)

Base.reinterpret(::Type{T}, oa::OneHotArray) where T = reinterpret(T, parent(oa))

# array interface
Base.size(oa::OneHotArray{K}) where K = (onehotsize(oa), size(oa.onehots)...)

Base.@propagate_inbounds function Base.getindex(oa::OneHotArray{K, N}, i, is::Vararg{Int, N}) where {K, N}
    @boundscheck checkbounds(oa, i, is...)
    return oa.onehots[is...][i]
end

Base.@propagate_inbounds function Base.getindex(oa::OneHotArray{K, N}, i::Colon, is::Vararg{Int, N}) where {K, N}
    @boundscheck checkbounds(oa, i, is...)
    return oa.onehots[is...]
end

Base.@propagate_inbounds function Base.getindex(oa::OneHotArray{K}, i::Colon, is...) where {K}
    @boundscheck checkbounds(oa, i, is...)
    return OneHotArray(oa.onehots[is...])
end

Base.@propagate_inbounds function Base.getindex(oa::OneHotArray{K}, i::Integer, is::Vararg{Colon, N}) where {K, N}
    @boundscheck checkbounds(oa, i, is...)
    return map(x->x[i], oa.onehots)
end

Base.similar(o::OneHotArray, ::Type{T}, dims::Dims{N}) where {T, N} = similar(o.onehots, T, dims)

Base.parent(o::OneHotArray) = o.onehots

# printing

function Base.summary(io::IO, oa::OneHotArray)
    join(io, size(oa), 'x')
    print(io, " OneHotArray{", onehotsize(oa), ", ", ndims(oa), ", ")
    Base.showarg(io, oa.onehots, true)
    print(io, "}")
end

# cat

is1(::Val{V}) where V = is1(V)
is1(v) = isone(v)
predecessor(::Val{V}) where V = Val(V-1)
predecessor(v) = Val(v - 1)

Base.vcat(xss::OneHot{K}...) where K = cat(xss...; dims=Val(1))
Base.hcat(xss::OneHot{K}...) where K = cat(xss...; dims=Val(2))

function Base.cat(xss::OneHot{K}...; dims) where K
    if is1(dims)
        @warn "concat OneHot{$K} along dimension 1."
        return Base._cat(Val(1), xss...)
    else
        yss = reshape(collect(xss), reverse(Base.rdims(predecessor(dims), axes(xss))))
        return OneHotArray(yss)
    end
end

Base.vcat(xss::OneHotArray{K}...) where K = cat(xss...; dims=Val(1))
Base.hcat(xss::OneHotArray{K}...) where K = cat(xss...; dims=Val(2))

function Base.cat(xss::OneHotArray{K}...; dims) where K
    if is1(dims)
        @warn "concat OneHotArray{$K} along dimension 1."
        return Base._cat(Val(1), xss...)
    else
        sdims = predecessor(dims)
        xidss = map(parent, xss)
        ret = cat(xidss...; dims=sdims)
        return OneHotArray(ret)
    end
end

# reshape

"reshape the underlying onehot array"
function ohreshape(parent::OneHotArray{K}, dims) where K
    onehots = parent.onehots
    return OneHotArray(reshape(onehots, dims))
end

function Base.reshape(parent::OneHotArray{K}, dims::Dims) where K
    isequal(prod(dims), length(parent)) || throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size $(length(parent))"))
    return isequal(K, first(dims)) ?
        ohreshape(parent, Base.tail(dims)) :
        Base._reshape(parent, dims)
end

Base.reshape(parent::OneHotArray{K}, dims::Tuple{Vararg{Union{Colon, Int}}}) where K = reshape(parent, Base._reshape_uncolon(parent, dims))

# adapt

import Adapt: adapt, adapt_structure
adapt_structure(T, oa::OneHotArray) = OneHotArray(adapt(T, oa.onehots))


text:
module DiffusionMoments

using ArgCheck
using DataFrames
using LsqFit
using Query
using StatsBase
using Statistics
using DataFramesMeta

# import Query: @where

# exportしたい関数一覧
export
    # msd.jl
    ensemble_time_average_msd,
    time_average_msd,
    ensemble_msd,
    ensemble_tamsd,
    # mme.jl
    mme,
    moment,
    # fit.jl
    fit,
    estimate_α

include("msd.jl")
include("_msd.jl")
include("mme.jl")
include("fit.jl")
include("utilities.jl")

end

max_stars_repo_name: SosUts/DiffusionMoments

max_issues_repo_name: SosUts/DiffusionMoments

max_forks_repo_name: SosUts/DiffusionMoments


text:
# Functions for performing discretization

using Discretizers
# using PyCall

# @pyimport astroML.density_estimation as de

function getfrequencies(values::Array{Float64,2}, mode)
	binids, numberofbins = getbinids(values, mode)
	frequencies = zeros(Int, (numberofbins, 1))
	for binid in binids
		frequencies[binid] += 1
	end
	return frequencies
end

function getjointfrequencies(valuesX::Array{Float64,2}, valuesY::Array{Float64,2}, mode)
	n = size(valuesX)[1]
	binidsX, numberofbinsX = getbinids(valuesX, mode)
	binidsY, numberofbinsY = getbinids(valuesY, mode)
	frequencies = zeros(Int, (numberofbinsY, numberofbinsX))
	for i in 1:n
		frequencies[binidsY[i], binidsX[i]] += 1
	end
	return frequencies
end
function getjointfrequencies(valuesX::Array{Float64,2}, valuesY::Array{Float64,2}, valuesZ::Array{Float64,2}, mode)
	n = size(valuesX)[1]
	binidsX, numberofbins = getbinids(valuesX, mode)
	binidsY = getbinids(valuesY, mode)[1]
	binidsZ = getbinids(valuesZ, mode)[1]
	frequencies = zeros(Int, (numberofbins, numberofbins, numberofbins))
	for i in 1:n
		frequencies[binidsZ[i], binidsY[i], binidsX[i]] += 1
	end
	return frequencies
end

function getbinids(values::Array{Float64,2}, mode)
	# This is a separate function because might want to offer different methods in the future
	function getnumberofbins(values)
		# size(values)[1] is n
		return round(Int, sqrt(size(values)[1]))
	end

	numberofbins = 0 # So can be accessed outside the loops
	n, numberofdimensions = size(values)
	binids = zeros(Int, size(values))
	for i in 1:numberofdimensions
		# If values are all the same, assign them all to bin 1
		min, max = extrema(values)
		if min == max
			numberofbins += 1
			binids[:, i:i] += convert(Array{Int}, values) + 1
		elseif mode == "uniformwidth"
			numberofbins += getnumberofbins(values)
			binids[:, i:i] += encode(LinearDiscretizer(binedges(DiscretizeUniformWidth(numberofbins), values)), values)
		elseif mode == "uniformwidth2"
			numberofbins += getnumberofbins(values) * 2
			binids[:, i:i] += encode(LinearDiscretizer(binedges(DiscretizeUniformWidth(numberofbins), values)), values)
		elseif mode == "uniformcount"
			numberofbins = getnumberofbins(values)
			try
				binids[:, i:i] += encode(LinearDiscretizer(binedges(DiscretizeUniformCount(numberofbins), reshape(values, length(values)))), values)
			catch
				binids[:, i:i] += encode(LinearDiscretizer(binedges(DiscretizeUniformWidth(numberofbins), values)), values)
			end
		elseif mode == "bayesianblocks"
			# edges = de.bayesian_blocks(reshape(values, length(values)))
			# binsizes = de.histogram(values, edges)[1]
			# binsizes = length(binsizes) == 0 ? [n] : binsizes
			# # The highest index in the current bin
			# highestindex = 0
			# # j is the index of the current bin
			# for j in 1:length(binsizes)
			# 	highestindex += binsizes[j]
			# 	binids[i:i, 1:highestindex] += 1
			# end
			# numberofbins = length(edges) - 1 <= 1 ? 1 : length(edges) - 1
		# TODO: Handle "mode doesn't exist" error
		end
	end
	return binids, numberofbins
end


text:
#=
  https://www.youtube.com/watch?v=e1Ykk_CqKTY&t=458s

  Probabilistic Programming: What It Is and How It Works - Noel Welsh

  We can see either 1, 2, or 3 cats.
  There are 3 different enticements:

  - Milkshake
  - Fish
  - Nothing

  And there are different probabilities how many cats there are given
  an enticement, see below.

  Now: We see 3 cats, what is the probability that it's a milkshake?

  The video got the following (for 3 cats):
   - milkshake: 0.42
   - fish: 0.04
   - nothing: 0.03

  Normalized to percentage (from the video):

  0.42/(0.42 + 0.04 + 0.03) milkshake
      0.85714285714285714286
  0.04/(0.42 + 0.04 + 0.03) fish
      0.081632653061224489796
  0.03/(0.42 + 0.04 + 0.03)  nothing
      0.061224489795918367347

  Here are two models.
  Model2 seems to be the correct one (i.e. corresponds to the video),
  i.e. where the enticements are mutual exclusive.
  In model1 we accept that fish and milkshake can coexist.

  See ~/cplint/cat.pl
      ~/blog/cat.blog
      ~/psi/cat.psi
=#

using Turing # , StatsPlots, DataFrames
include("jl_utils.jl")
# using ReverseDiff, Zygote, Tracker
# Turing.setadbackend(:reversediff)
# Turing.setadbackend(:zygote)
# Turing.setadbackend(:tracker)

# This is not correct
@model function cat1()
    NoCat = 0
    OneCat = 1
    TwoCats = 2
    ThreeCats = 3

    # probability of an enticement in the garden
    milkshake      ~ Bernoulli(0.6)
    fish           ~ Bernoulli(0.1)
    no_enticement  ~ Bernoulli(0.3)
    println("milkshake:$milkshake fish:$fish no_enticement:$no_enticement")

    # Number of cats per enticement
    vs = [NoCat,OneCat,TwoCats,ThreeCats]
    vs_cat ~ DiscreteUniform(NoCat,ThreeCats)
    if !no_enticement && milkshake
        cat ~ Categorical(simplex([0.0,0.1,0.2,0.7]))
    elseif !no_enticement && fish
        cat ~ Categorical(simplex([0.0,0.2,0.3,0.4]))
    elseif no_enticement && !fish && !milkshake
        cat ~ Categorical(simplex([0.0,0.6,0.3,0.1]))
    else
        cat ~ Categorical(simplex([1.0,0.0,0.0,0.0]))
    end
    println("cat before:$cat  (vs[cat]: $(vs[cat]))")
    true ~ Dirac(vs_cat == vs[cat]) # || begin Turing.@addlogprob! -Inf;  end
    true ~ Dirac(vs_cat == ThreeCats) #  || begin Turing.@addlogprob! -Inf;  end
    println("cat after:$cat")

end

#=
model = cat1()
num_chains = 4
# chains = sample(model, Prior(), MCMCThreads(), 1000, num_chains)
chains = sample(model, MH(), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, MH(), 1000)

display(chains)
=#

# This is correct (according to the video)
#=
Summary Statistics
            parameters      mean       std   naive_se      mcse         ess      rhat
                Symbol   Float64   Float64    Float64   Float64     Float64   Float64

                  cats    3.0000    0.0000     0.0000    0.0000         NaN       NaN
            enticement    1.2009    0.5322     0.0013    0.0090   2218.9994    1.0026
       enticement_fish    0.0783    0.2686     0.0007    0.0043   2703.4937    1.0026
  enticement_milkshake    0.8600    0.3470     0.0009    0.0058   2345.7527    1.0017
    enticement_nothing    0.0617    0.2407     0.0006    0.0041   2176.3555    1.0048
=#
enticements = []
enticements_dict = Dict(1=>"milkshake",2=>"fish",3=>"no_enticement")
@model cat2() = begin

    milkshake = 1
    fish = 2
    no_enticement = 3

    # Priors
    # probability of an enticement in the garden
    enticement ~ Categorical(simplex([0.6,0.1,0.3]))
    # Note: This seems a little too complex...
    enticement_milkshake ~ Bernoulli(0.5)
    enticement_fish      ~ Bernoulli(0.5)
    enticement_nothing   ~ Bernoulli(0.5)
    cats ~ DiscreteUniform(1,3)
    if enticement == milkshake
        cats ~ Categorical(simplex([0.1,0.2,0.7]))
        enticement_milkshake ~ Bernoulli(1.0)
        enticement_fish      ~ Bernoulli(0.0)
        enticement_nothing   ~ Bernoulli(0.0)
    elseif enticement == fish
        cats ~ Categorical(simplex([0.2,0.4,0.4]))
        enticement_milkshake ~ Bernoulli(0.0)
        enticement_fish      ~ Bernoulli(1.0)
        enticement_nothing   ~ Bernoulli(0.0)
    else
        cats ~ Categorical(simplex([0.6,0.3,0.1]))
        enticement_milkshake ~ Bernoulli(0.0)
        enticement_fish      ~ Bernoulli(0.0)
        enticement_nothing   ~ Bernoulli(1.0)
    end

    # println("enticement:$(enticements_dict[$enticement]) cats:$cats")

    # We observe 3 cats
    true ~ Dirac(cats == 3) # || begin Turing.@addlogprob! -Inf; return; end
    # push!(enticements,[enticements_dict[enticement],cats])

    return enticement
end

model = cat2()
num_chains = 4
# chains = sample(model, Prior(), MCMCThreads(), 1000, num_chains)
chains = sample(model, MH(), MCMCThreads(), 40_000, num_chains)
# chains = sample(model, MH(), 40_000)

# chains = sample(model, PG(20), MCMCThreads(), 1000, num_chains)

display(chains)

println("\nEnticement: milkshake:1 fish:2 no_enticement:3")
gen = generated_quantities(model, chains)
show_var_dist_pct(gen, 120)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using Pkg
tempdir = mktempdir()
Pkg.activate(tempdir)
Pkg.develop(PackageSpec(path = joinpath(@__DIR__, "..")))
Pkg.add(["BenchmarkTools", "PkgBenchmark", "MathOptInterface"])
Pkg.resolve()

using Convex: Convex, ProblemDepot
using BenchmarkTools
using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

const SUITE = BenchmarkGroup()

problems = [
    "constant_fix!_with_complex_numbers",
    "affine_dot_multiply_atom",
    "affine_hcat_atom",
    "affine_trace_atom",
    "exp_entropy_atom",
    "exp_log_perspective_atom",
    "socp_norm_2_atom",
    "socp_quad_form_atom",
    "socp_sum_squares_atom",
    "lp_norm_inf_atom",
    "lp_maximum_atom",
    "sdp_and_exp_log_det_atom",
    "sdp_norm2_atom",
    "sdp_lambda_min_atom",
    "sdp_sum_largest_eigs",
    "mip_integer_variables",
]

SUITE["formulation"] = ProblemDepot.benchmark_suite(problems) do problem
    model = MOIU.MockOptimizer(MOIU.Model{Float64}())
    return Convex.load_MOI_model!(model, problem)
end


text:
module TestPointSpreadFunctions

using Test
using PointSpreadFunctions
using PointSpreadFunctions:
    AbstractPSF,
    finite_and_nonnegative,
    finite_and_positive,
    throw_bad_argument,
    to_float,
    to_int
import PointSpreadFunctions:
    check_structure,
    getfwhm,
    immutable,
    isvalid,
    parameters
using PointSpreadFunctions.Fitting:
    objfun

struct MyPSF <: AbstractPSF{3}
    prm::NTuple{3,Float64}
    MyPSF(a::Real, b::Real, c::Real) =
        new(map(x -> to_float(Float64, x), (a,b,c)))
end

struct MyPSF2 <: AbstractPSF{0}; end

parameters(P::MyPSF) = P.prm
parameters(P::MyPSF2) = ()
getfwhm(P::MyPSF) = P[1]
getfwhm(P::MyPSF2) = 1.0
# Purposely not extend isvalid and check_structure for MyPSF2 to check
# default implementations.
isvalid(P::MyPSF) = (_check_fwhm(P) && _check_b(P) && _check_c(P))
function check_structure(P::MyPSF)
    _check_fwhm(P) || throw_bad_argument("invalid FWHM (",getfwhm(P),")")
    _check_b(P) || throw_bad_argument("invalid parameter B (",_get_b(P),")")
    _check_c(P) || throw_bad_argument("invalid parameter C (",_get_c(P),")")
end

_get_b(P::MyPSF) = P[2]
_get_c(P::MyPSF) = P[3]

_check_fwhm(P::MyPSF) = finite_and_nonnegative(getfwhm(P))
_check_b(P::MyPSF) = finite_and_positive(_get_b(P))
_check_c(P::MyPSF) = (c = _get_c(P); isfinite(c) && 0 < c < 1)

@testset "PSF models and fitting" begin

    for T in (Int8, Int16, Int32, Int64)
        @test to_int(T(1)) === one(Int)
    end
    for T in (Float32, Float64)
        @test to_float(T, 1) === one(T)
    end

    @test finite_and_positive(NaN) == false
    @test finite_and_positive(+Inf) == false
    @test finite_and_positive(-1) == false
    @test finite_and_positive(0) == false
    @test finite_and_positive(π) == true

    @test finite_and_nonnegative(NaN) == false
    @test finite_and_nonnegative(+Inf) == false
    @test finite_and_nonnegative(-1) == false
    @test finite_and_nonnegative(0) == true
    @test finite_and_nonnegative(π) == true

    @test isvalid(MyPSF2()) == true
    @test check_structure(MyPSF2()) == true

    @test_throws ArgumentError throw_bad_argument("oops!")
    @test_throws ArgumentError throw_bad_argument("oops! (", π, ")")

    @test_throws ArgumentError check_structure(MyPSF(-1.0,1,0.3))
    @test_throws ArgumentError check_structure(AiryPSF(-1.0))
    @test_throws ArgumentError check_structure(AiryPSF(1.0, 1.1))
    @test_throws ArgumentError check_structure(GaussianPSF(-1.0))
    @test_throws ArgumentError check_structure(CauchyPSF(Inf))
    @test_throws ArgumentError check_structure(MoffatPSF(0.0,0.5))
    @test_throws ArgumentError check_structure(MoffatPSF(1.0,-0.2))
    @test isvalid(AiryPSF(1.0, 1.1)) == false
    @test isvalid(AiryPSF(1.0, 0.3)) == true
    @test isvalid(GaussianPSF(-1.0)) == false
    @test isvalid(GaussianPSF(1.0)) == true
    @test isvalid(CauchyPSF(NaN)) == false
    @test isvalid(CauchyPSF(1.0)) == true
    @test isvalid(MoffatPSF(0.0,0.5)) == false
    @test isvalid(MoffatPSF(1.0,-0.2)) == false
    @test isvalid(MoffatPSF(1.0,0.4)) == true

    # Tests for any PSF models.
    for psf in (AiryPSF(1.1), AiryPSF(1.2, 0.3),
                CauchyPSF(1.8), GaussianPSF(2.7),
                MoffatPSF(1.4, 0.5))
        # Check peak PSF is 1 at (0,0).
        @test psf(0) == 1
        @test psf(0, 0) == 1

        # Check FWHM.
        let fwhm = getfwhm(psf)
            @test psf(fwhm/2) ≈ 0.5
            @test psf(fwhm/2, 0) ≈ 0.5
            @test psf(0, fwhm/2) ≈ 0.5
        end

        # Check show extension.
        io = IOBuffer()
        show(io, psf)
        str = String(take!(io))
        tmp = repr(typeof(psf))
        # Note: `findlast(char,str)` is not available until Julia 1.3 and
        # is not in `Compat` so we have to do it ourselves.
        i = findlast(==('.'), tmp)
        pfx = (i === nothing ? tmp : tmp[i+1:end])
        @test startswith(str, pfx*"(")
        @test endswith(str, ")")

        for (x,y) in ((1.23, -0.12), (-0.7, sqrt(2)))
            @test psf(x, y) ≈ psf(hypot(x,y))
            for T in (Float32, Float64)
                @test psf(T, x, y) ≈ T(psf(x,y))
                @test psf(T, hypot(x,y)) ≈ T(psf(hypot(x,y)))
                @test psf(T, x, y) ≈ psf(T, hypot(x,y))
            end
        end
    end

    # Specific tests for AiryPSF.
    let psf = AiryPSF(1.0)
        @test psf[:] === (1.0, 0.0)
        @test getfwhm(psf) ≈ PointSpreadFunctions.AIRY_FWHM
        @test findzeros(psf, 1:5) ≈ [PointSpreadFunctions.AIRY_FIRST_ZERO,
                                     PointSpreadFunctions.AIRY_SECOND_ZERO,
                                     PointSpreadFunctions.AIRY_THIRD_ZERO,
                                     PointSpreadFunctions.AIRY_FOURTH_ZERO,
                                     PointSpreadFunctions.AIRY_FIFTH_ZERO] atol=0.01
    end

    # Test fitting.
    x = 1:15
    y = 1:17
    x0 = 6.2
    y0 = 7.3
    α = 3.1
    σ = 0.03
    xinit = 6
    yinit = 7
    roi = (2:last(x)-1, :)
    for (psf, initialpsf) in (
        (GaussianPSF(5.4), GaussianPSF(5)),
        (CauchyPSF(5.4), CauchyPSF(4)),
        (AiryPSF(5.4, 0.33), AiryPSF(5, 0.3)),
        (MoffatPSF(5.2, 1.5), MoffatPSF(5, 1)),
    )

        mdl = psf((x .- x0), (y .- y0)')
        dat = α*mdl + σ*randn(size(mdl))
        wgt = ones(eltype(dat), size(dat))

        if psf isa Union{AiryPSF,GaussianPSF,CauchyPSF,MoffatPSF}
            # Check that PSF has radial symmetry.
            @test psf(hypot.((x .- x0), (y .- y0)')) ≈ mdl
        end

        # Check conversion of arguments in objective function with given
        # alpha.
        @test objfun(     dat, 1, psf, (0, 0)) == objfun(     dat, 1.0, psf, (0.0, 0.0))
        @test objfun(wgt, dat, 1, psf, (0, 0)) == objfun(wgt, dat, 1.0, psf, (0.0, 0.0))

        # Check conversion of arguments in objective function with
        # automatic alpha.
        @test objfun(     dat, psf, (0, 0)) == objfun(     dat, psf, (0.0, 0.0))
        @test objfun(wgt, dat, psf, (0, 0)) == objfun(wgt, dat, psf, (0.0, 0.0))

        # Check that the nonnegative constraints works.
        @test objfun(     -dat, psf, (0, 0), true) == 0
        @test objfun(wgt, -dat, psf, (0, 0), true) == 0

        # Check the versions of the objective function that behave as if
        # the model is zero everywhere.
        @test objfun(     dat) ≈ objfun(     dat, 0, psf, (0, 0))
        @test objfun(wgt, dat) ≈ objfun(wgt, dat, 0, psf, (0, 0))
        @test objfun(     dat, false) ≈ objfun(     dat, true)
        @test objfun(wgt, dat, false) ≈ objfun(wgt, dat, true)

        # Fit all PSF parameters.
        for (psf1, pos1) in (
            PointSpreadFunctions.fit(initialpsf, (xinit,yinit), dat;
                                     nonnegative=true),
            PointSpreadFunctions.fit(initialpsf, (xinit,yinit), wgt, dat;
                                     nonnegative=true),
            PointSpreadFunctions.fit(initialpsf, (xinit,yinit), dat, roi;
                                     nonnegative=true),
            PointSpreadFunctions.fit(initialpsf, (xinit,yinit), wgt, dat, roi;
                                     nonnegative=true),
        )
            @test pos1[1] ≈ x0     rtol=0.05 atol=0.0
            @test pos1[2] ≈ y0     rtol=0.05 atol=0.0
            for i in 1:length(psf[:])
                @test psf1[i] ≈ psf[i] rtol=0.1  atol=0.0
            end
        end

        # Fit only the PSF position.
        for (psf1, pos1) in (
            PointSpreadFunctions.fit(immutable(psf), (xinit,yinit), dat;
                                     nonnegative=true),
            PointSpreadFunctions.fit(immutable(psf), (xinit,yinit), wgt, dat;
                                     nonnegative=true),
            PointSpreadFunctions.fit(immutable(psf), (xinit,yinit), dat, roi;
                                     nonnegative=true),
            PointSpreadFunctions.fit(immutable(psf), (xinit,yinit), wgt, dat, roi;
                                     nonnegative=true),
        )
            @test pos1[1] ≈ x0     rtol=0.05 atol=0.0
            @test pos1[2] ≈ y0     rtol=0.05 atol=0.0
            @test psf1 === psf
        end
    end
end

end # module


text:
let
	logfile = "tmp_out.log"
	errfile = "tmp_err.log"
	redirect_to_files(logfile,errfile) do
		limit = 2
		msg = "This message should be printed"
		@log_info(limit,3,msg)
		@log_info(limit,1+2,msg)
		@log_info(limit,limit-1,msg) # nothing should be printed
		@log_info(limit,limit,msg) # nothing should be printed
		@log_info(limit,3,split(msg," ")...)
		@log_info(limit,3)
	end
	@test countlines(logfile) == 4
	rm(logfile)
	rm(errfile)
end


text:
struct TimeHomogeneousForwardCorrelation end

function evolved_matrices(::Type{TimeHomogeneousForwardCorrelation}, fwdCorrelation::Matrix{Float64})
  numberOfRates = size(fwdCorrelation)[1]
  correlations = Matrix{Float64}[zeros(numberOfRates, numberOfRates) for i = 1:numberOfRates]
  # correlations = fill(zeros(numberOfRates, numberOfRates), numberOfRates)

  @inbounds @simd for k in eachindex(correlations)
    # proper diagonal values
    for i = k:numberOfRates
      correlations[k][i, i] = 1.0
    end
    # copy only time homogeneous values
    for i = k:numberOfRates, j = k:i
      correlations[k][i,j] = correlations[k][j,i] = fwdCorrelation[i-(k-1), j-(k-1)]
    end
  end
  return correlations
end


text:
#  ------------------------------------------------------------------------------------------
#  Copyright (c) Microsoft Corporation. All rights reserved.
#  Licensed under the MIT License (MIT). See LICENSE in the repo root for license information.
#  ------------------------------------------------------------------------------------------

using Serialization
export BlobRef, BatchFuture, BlobFuture, fetch!, remote_reduction#, fetchreduce_batch
import Base.fetch

#######################################################################################################################
# Futures

mutable struct BlobRef
    name::Union{String, Tuple}
end

mutable struct BatchFuture
    container::String
    blob
end

function fetch(arg::BatchFuture)
    if haskey(ENV, "AZ_BATCH_TASK_SHARED_DIR")
        iostream = open(join([ENV["AZ_BATCH_TASK_SHARED_DIR"], "/", arg.blob.name]), "r")
    else
        iostream = open(arg.blob.name, "r")
    end
    data = deserialize(iostream)
    close(iostream)
    return data
end

function fetch!(arg::BatchFuture)
    if haskey(ENV, "AZ_BATCH_TASK_SHARED_DIR")
        iostream = open(join([ENV["AZ_BATCH_TASK_SHARED_DIR"], "/", arg.blob.name]), "r")
    else
        iostream = open(arg.blob.name, "r")
    end
    arg.blob = deserialize(iostream)
    close(iostream)
    return arg.blob
end

mutable struct BlobFuture
    container::String
    blob
    client_index::Integer
end

function fetch(arg::BlobFuture)

    num_files = length(arg.blob.name)
    out_files = []
    for blob in arg.blob.name

        # Fetch blob and add to collection
        if haskey(ENV, "AZ_BATCH_TASK_SHARED_DIR")
            iostream = open(join([ENV["AZ_BATCH_TASK_SHARED_DIR"], "/", blob]), "r")
        else
            iostream = open(blob, "r")
        end
        push!(out_files, deserialize(iostream))
        close(iostream)
    end

    if num_files > 1
        data = tuple(out_files...)
    else
        data = out_files[1]
    end
    return data
end

function fetch!(arg::BlobFuture)

    num_files = length(arg.blob.name)
    out_files = []
    for blob in arg.blob.name

        # Fetch blob and add to collection
        if haskey(ENV, "AZ_BATCH_TASK_SHARED_DIR")
            iostream = open(join([ENV["AZ_BATCH_TASK_SHARED_DIR"], "/", blob]), "r")
        else
            iostream = open(blob, "r")
        end
        push!(out_files, deserialize(iostream))
        close(iostream)
    end

    if num_files > 1
        arg.blob = tuple(out_files...)
    else
        arg.blob = out_files[1]
    end
    return arg.blob
end


text:
#################################################################################################
#
# demo2: Convert a image from .png to .jpg format (with compression)
#
#################################################################################################

# Julia filename with full path
inputfile = joinpath(Pkg.dir("OpenCV"), "./test/images/lena.png")
outfile = joinpath(Pkg.dir("OpenCV"), "./test/images/lena.jpeg")

cxx"""
void imageConversion(const char *inputfile, const char *outfile)
{
    const std::string ifname = inputfile;
    const std::string ofname = outfile;

    // Open video file
     cv::Mat img = cv::imread(ifname);

     if (img.empty())
     {
          std::cout<< "Cannot read the file!";
     }

     //cv::cvtColor(img, img, cv::COLOR_BGR2GRAY);

     cv::namedWindow("Original", cv::WINDOW_AUTOSIZE);
     cv::imshow("Original", img);  // jpeg, actually same image displayed in both windows,

     std::vector<int> compression_params; // vector that stores compression params
     compression_params.push_back(cv::IMWRITE_JPEG_QUALITY);
     compression_params.push_back(98); // specify % compression quality

     //compression_params.push_back(cv::IMWRITE_PNG_COMPRESSION);
     //compression_params.push_back(9);

     // not working!
     bool Success = cv::imwrite(ofname, img, compression_params);
     //

     if(!Success)
     {
         std::cout<< "Failed to save the image!";
     }

     cv::namedWindow("Modified", cv::WINDOW_AUTOSIZE);
     cv::imshow("Modified", img); //  jpeg

     cv::waitKey(0);

     cv::destroyAllWindows();
}
"""

@cxx imageConversion(pointer(inputfile), pointer(outfile))