text:
# Test gausslaguerre.jl
import FastGaussQuadrature.getUQ
import FastGaussQuadrature.asyAiry
import FastGaussQuadrature.asyAirygen
import FastGaussQuadrature.asyBessel
import FastGaussQuadrature.asyBesselgen
import FastGaussQuadrature.asyBulk
import FastGaussQuadrature.asyBulkgen
import FastGaussQuadrature.polyAsyRH
import FastGaussQuadrature.polyAsyRHgen

# Test integration
tol = 3.e-13
tolGen = 2.e-9
# Possibly set srand(...) to get the same results for integration of a random polynomial.
ns = [42; 251; 5000];
for wei = 1:3
    if wei == 1; alpha = 0.0;qm = 1.0; m = 1;
    elseif wei == 2; alpha = pi+0.1; qm = 1.0; m = 1;
    elseif wei == 3; alpha = 0.67; qm = 0.34; m = 4;
    end

    for ni = 1:length(ns)
        n = ns[ni]
	# Make sure the integrand does not overflow at x[n]
        d = min(2*n-1, ceil(Int64, log(realmax(Float64)-4*n)/n ) )
        ra = rand(d+1)
        function fct(x)
            if length(x) == 1
                return sum(ra.*(x.^(0:d)))
            end
            result = NaN*x
            for i = 1:length(x)
                result[i] = sum(ra.*(x[i].^(0:d)))
            end
            return result
        end
        fcti(x) = fct(x).*exp(-qm*x.^m).*x^(alpha)
        # Test integration of a random polynomial which should be integrated exactly:
        exa = quadgk(fcti, 0, Inf, reltol = tol)[1]
        for methi = 1:ifelse(ni == 1, 1, 5)
            if (ni == 1)
                # Golub-Welsch is the accurate algorithm for low n but does not allow m or qm != 1.
                if wei == 3; continue; end
                println("    α = $alpha, n = $ni, Golub-Welsch")
                x, w = gausslaguerre( n, alpha )
            elseif methi == 1
                println("    α = $alpha, n = $ni, default")
                x, w = gausslaguerre( n, alpha, "default", qm, m)
            elseif methi == 2
                if wei == 3; continue; end # General weight
                println("    α = $alpha, n = $ni, RH")
                x, w = gausslaguerre( n, alpha, "RH", qm, m)
            elseif methi == 3
                if wei == 3; continue; end # General weight
                println("    α = $alpha, n = $ni, RHW")
                x, w = gausslaguerre( n, alpha, "RHW", qm, m)
            elseif methi == 4
                println("    α = $alpha, n = $ni, gen")
                x, w = gausslaguerre( n, alpha, "gen", qm, m)
            elseif methi == 5
                println("    α = $alpha, n = $ni, genW")
                x, w = gausslaguerre( n, alpha, "genW", qm, m)
            end
            # Multiply the tolerance by sqrt(n) because the error on each node/meaningful weight should be about 1e-14
            @test abs(dot(x,w) - gamma((alpha+2)/m)*qm^(-(alpha+2)/m)/m )/gamma((alpha+2)/m)/qm^(-(alpha+2)/m)*m <= sqrt(n)*ifelse(wei == 3, tolGen, tol)
            @test abs(dot(x.^2,w) - gamma((alpha+3)/m)*qm^(-(alpha+3)/m)/m )/gamma((alpha+3)/m)/qm^(-(alpha+3)/m)*m <= sqrt(n)*ifelse(wei == 3, tolGen, tol)
            if (n == 42) && (wei == 1)
                @test abs(x[37] - 98.388267163326702) < tol
                @test abs(w[7] - 0.055372813167092) < tol
            elseif (n == 251) && (wei == 1)
                @test abs(x[37] - 13.309000189442097) < tol
                @test abs(w[3] - 0.050091759039996) < tol
            end
            @test abs(dot(fct(x), w) - exa)/abs(exa) < sqrt(n)*ifelse(wei == 3, tolGen, tol)
	    end
    end
end

# Test whether higher-order terms correspond between "RH(W)", "gen(W)" and exact results computed in high precision.
n = NaN
tolHot = eps(Float64)*10
np = 200
tolEx = np*4e-15
T = ceil(Int64, 34/log(np) )
@test T == 7
for alpha = [0.0; 4.15]
    println("    α = $alpha")
    UQ = getUQ(alpha, 1.0, 1, T+1)

    z = 0.95 # > 3.7/4
    mnxi = 2*np*( sqrt(z)*sqrt(1 - z) - acos(sqrt(z) ) )
    fn = (mnxi*3im/2)^(2/3)
    aRHa = asyAiry(np, 4*np*z, alpha, T)
    @test abs(asyAirygen(np, z, alpha, T, 1.0, 1, UQ, fn, true) -aRHa)./abs(aRHa) < tolHot
    # Not using Q has an explicit substraction of ceil(3*T/2)-order poles, so loosen the tolerance.
    @test abs(asyAirygen(np, z, alpha, T, 1.0, 1, UQ, fn, false, mnxi*1im/np) -aRHa)./abs(aRHa) < 300*tolHot

    z = 0.01 # < sqrt(np)/4/np
    mnxi = 2*np*( sqrt(z)*sqrt(1 - z) - acos(sqrt(z) ) )
    aRHe = asyBessel(np, 4*np*z, alpha, T)
    @test abs(asyBesselgen(np, z, alpha, T, 1.0, 1, UQ, mnxi + pi*np, true) -aRHe)./abs(aRHe) < 20*tolHot
    @test abs(asyBesselgen(np, z, alpha, T, 1.0, 1, UQ, mnxi + pi*np, false) -aRHe)./abs(aRHe) < 20*tolHot

    z = 0.35
    mnxi = 2*np*( sqrt(z)*sqrt(1 - z) - acos(sqrt(z) ) )
    aRHu = asyBulk(np, 4*np*z, alpha, T)
    @test abs(asyBulkgen(np, z, alpha, T, 1.0, 1, UQ, mnxi) -aRHu)./abs(aRHu) < tolHot

    a = alpha
    factorp = (1/3840*a^10 - 5/2304*a^9 + 11/2304*a^8 + 7/1920*a^7 - 229/11520*a^6 + 107/34560*a^5 + 2653/103680*a^4 - 989/155520*a^3 -         3481/311040*a^2 + 139/103680*a + 9871/6531840)/np^5
    factorp = factorp + (1/384*a^8 - 1/96*a^7 + 1/576*a^6 + 43/1440*a^5 - 5/384*a^4 - 23/864*a^3 + 163/25920*a^2 + 31/6480*a - 139/155520)/np^4
    factorp = factorp + (1/48*a^6 - 1/48*a^5 - 1/24*a^4 + 5/144*a^3 + 1/36*a^2 - 1/144*a - 31/6480)/np^3
    factorp = factorp + (1/8*a^4 + 1/12*a^3 - 1/24*a^2 + 1/72)/np^2
    factorp = factorp + (1/2*a^2 + 1/2*a + 1/6)/np^1 +1
    # Add an asymptotic approximation of gamma(alpha+n+1)/gamma(n+1)
    factorp = (4*np)^(-alpha/2-1/2)*sqrt(1/2/pi)/sqrt(factorp)*sqrt(np^alpha*(1+alpha*(alpha+1)/2/np + alpha*(3*alpha^3+2*alpha^2-3*alpha-2)/24/np^2 + alpha^2*(alpha+1)^2*(alpha^2-3*alpha+2)/48/np^3 + alpha*(15*alpha^7-60*alpha^6-10*alpha^5+192*alpha^4 -25*alpha^3-180*alpha^2+20*alpha+48)/5760/np^4 + alpha^2*(alpha+1)^2*(3*alpha^6-31*alpha^5+109*alpha^4-125*alpha^3-88*alpha^2+276*alpha-144)/11520/np^5) )
    if alpha == 0.0
        @test norm([+6.975493821675553e-05  ,  -8.294589753829749e-05/1im]-vec(UQ[1,:,6,1,2]))/norm(UQ[1,:,6,1,2]) < tolHot
        @test norm([-6.975493821675553e-05  ,  +1.401453805748335e-05/1im]-vec(UQ[1,:,6,1,1]))/norm(UQ[1,:,6,1,1]) < tolHot*10
        # We choose to test the constants most prone to errors.
        @test norm([-0.01196102063075393    ,  +0.1904949571852236/1im]-vec(UQ[1,:,6,3,4]))/norm(UQ[1,:,6,1,4]) < tolHot*1e3
        @test norm([+1.309628097160176e-05  ,  -0.001326131973043531/1im]-vec(UQ[1,:,6,3,3]))/norm(UQ[1,:,6,1,3]) < tolHot*10

        # Compare with results obtained from a high precision calculation of the Laguerre polynomials with standard normalisation.
        @test abs(-4.9327289759964085492326419516911524182405e163/exp(0.95*4*np/2)/factorp -aRHa)/abs(aRHa) < tolEx
        @test abs(-1.6224971613979232593218524827791333158448e59/exp(0.35*4*np/2)/factorp -aRHu)/abs(aRHu) < tolEx
        @test abs(-3.9144031220161646274486036731469804751328/exp(0.01*4*np/2)/factorp -aRHe)/abs(aRHe) < tolEx
    else
        @test abs(5/49152*alpha^8 - 35/49152*alpha^7 + 67/49152*alpha^6 + 11/36864*alpha^5 - 1529/589824*alpha^4 + 1891/2359296*alpha^3 + 26827/26542080*alpha^2 - 109/524288*alpha - 190867/1698693120  -UQ[1,1,4,1,2])/abs(UQ[1,1,4,1,2]) < 6*tolHot
        @test abs(-1/24576*alpha^8 - 17/24576*alpha^7 - 107/147456*alpha^6 + 1751/983040*alpha^5 + 1535/1179648*alpha^4 - 853/884736*alpha^3 - 7583/26542080*alpha^2 + 14131/141557760*alpha - 12829/1274019840 -1im*4^alpha*UQ[1,2,4,2,1])/abs(4^alpha*UQ[1,2,4,2,1]) < tolHot
        @test abs(-23/259187003712000*alpha^21 + 353/76796149248000*alpha^20 - 333791/4561691265331200*alpha^19 + 87631/600222534912000*alpha^18 + 24928777/4001483566080000*alpha^17 - 205079/4803701760000*alpha^16 - 121008773/882680198400000*alpha^15 + 221720011/117690693120000*alpha^14 - 11808705347/16991593819200000*alpha^13 - 119431297679/3485455142400000*alpha^12 + 100119618541/1742727571200000*alpha^11 + 3526062948193/12357522777600000*alpha^10 - 83939498901133/123575227776000000*alpha^9 - 48402313811119/41191742592000000*alpha^8 + 4645283107358809/1050389436096000000*alpha^7 - 91277362389443/40399593696000000*alpha^6 - 49328830551721603/7622617782780000000*alpha^5 + 23985609069512791/1960101715572000000*alpha^4 + 2660362470248081/3587114250720000000*alpha^3 - 779839198817230403/60980942262240000000*alpha^2 + 361055235303455299/3347022171893400000000*alpha + 91136762727842675269/36817243890827400000000 - 1im*4^alpha*UQ[1,2,4,5,3])/abs(4^alpha*UQ[1,2,4,5,3]) < 3e-11 # First term = 8.9e-14 times alpha^21=9.5e12, and some integers in the formula cannot be represented in Int64, so lower accuracy can be expected
        @test abs(1/362880*alpha^13 - 1/453600*alpha^12 - 71/435456*alpha^11 + 1361/7257600*alpha^10 + 625/217728*alpha^9 - 21569/7257600*alpha^8 - 1337/64800*alpha^7 + 291101/21772800*alpha^6 + 9329/170100*alpha^5 + 3224237/59875200*alpha^4 - 272927/3628800*alpha^3 - 6823639/20528640*alpha^2 + 138793/2721600*alpha + 40018577/179625600 - UQ[1,1,3,4,4])/abs(UQ[1,1,3,4,4]) < 40*tolHot

        @test abs(-3.7536142288312445691754139973055317263278e162/exp(0.95*4*np/2)/factorp -aRHa)/abs(aRHa) < tolEx
        @test abs(1.3002334410655323780033551365563985116651e59/exp(0.35*4*np/2)/factorp -aRHu)/abs(aRHu) < tolEx
        @test abs(-2395.51952258921326919097391744358588135/exp(0.01*4*np/2)/factorp -aRHe)/abs(aRHe) < tolEx
    end

    for n = ceil(Int64, linspace(111, 5000, 5) )
        for z = linspace(1e-4, 0.99, 10)
            aRH = polyAsyRH(n, 4*n*z, alpha, T)
            @test abs(polyAsyRHgen(n, 4*n*z, alpha, T, 1.0, 1, UQ) -aRH)/abs(aRH) < 2*tolEx
        end
    end
end


text:
module XCB

using DocStringExtensions

@template (FUNCTIONS, METHODS, MACROS) =
    """
    $(DOCSTRING)
    $(TYPEDSIGNATURES)
    """

@template TYPES =
    """
    $(DOCSTRING)
    $(TYPEDEF)
    $(TYPEDSIGNATURES)
    $(TYPEDFIELDS)
    """

using Parameters
using WindowAbstractions
import WindowAbstractions: set_title,
                           set_icon_title,
                           set_extent,
                           extent,
                           terminate_window!,
                           map_window,
                           unmap_window,
                           wait_for_event,
                           poll_for_event,
                           attach_graphics_context!,
                           get_window,
                           get_window_symbol,

                           MouseEvent,
                           MouseState,
                           KeyEvent,
                           KeyModifierState,
                           KeyContext,
                           KeyCombination,
                           EventDetails,
                           KeySymbol

include(joinpath(@__DIR__, "..", "gen", "Libxcb.jl"))
include(joinpath(@__DIR__, "..", "gen", "Libxkb.jl"))

const xcb = Libxcb
using .Libxcb
using .Libxkb

include("exceptions.jl")
include("connection.jl")
include("window.jl")
include("inputs.jl")
include("xkb.jl")
include("window_handler.jl")
include("graphics.jl")
include("testing.jl")
include("events.jl")

export xcb,
       key,
       Connection,
       Setup,
       check,
       check_flush,
       flush,
       XCBWindow,
       GraphicsContext,
       XWindowHandler,
       set_callbacks!,
       set_attributes,
       get_window,
       get_window_symbol,
       event_details_xkb,
       @check,
       @flush,
       keystroke_info,
       keycode,
       key_name,
       keysym_string,
       key_symbol_from_keysym_string,
       key_event_from_name,
       unsafe_load_event,
       send_event,
       event_xcb,
       action_xcb,
       state_xcb,
       event_type_xcb,
       button_xcb

end # module

ext: jl


text:
using Documenter
using MieScattering

makedocs(
    sitename = "MieScattering",
    authors="Adam Fekete <adam.fekete@unamur.be> and contributors",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        canonical="https://fekad.github.io/MieScattering.jl"
    ),
    modules = [MieScattering],
    pages=[
        "Home" => "index.md",
        "Mie" => "mie.md",
        "Rayleigh" => "rayleigh.md"
    ]
)

deploydocs(
    repo = "github.com/fekad/MieScattering.jl",
    devbranch = "main"
)

max_stars_repo_name: fekad/jlmie

max_issues_repo_name: fekad/jlmie

max_forks_repo_name: fekad/jlmie


text:
export bcd

"""
function bcd(param,w;wtol,maxIter,maxIterQP,storeInterm,out,pwr)

Block coordinate descent method for solving

min_{M,w} || M*w - d ||_{Gamma,pwr} + lambda ||M||_0 s.t. w in Omega_W, M in Omega_M

Input:

    1. param       - StrainRecon param
    2. w           - initial starting point w for BCD

Keyword arguments:

    3. wtol        - tolerance for breaking out of BCD when |w_old - w_current| < wtol
    4. maxIter     - maximum number of BCD trials
    5. maxIterQP   - maximum number of iterations for solving quadratic problem in findw (passed to asqp)
    6. storeInterm - return last (M,w) from each BCD trial (storeInterm = 0 or 2)
                     return all (M,w) from each BCD trial (storeInterm = 1)
    7. out         - print intermediate results of findw (passed to asqp) if out > 1, else print nothing
    8. pwr           - exponent of the likelihood l_p norm

Output:

    1. Mopt - optimal M (storeInterm = 0 or 2) or all intermediate history of M (1)
    2. wopt - optimal w (storeInterm = 0 or 2) or all intermediate history of w (1)
    3. his  - history of misfit, |M_old-M_current|, |w_old-w_current| for optimal (M,w) (storeInterm = 0 or 2)
              history of misfit, |M_old-M_current|, |w_old-w_current| for all intermediate history of (M,w) (1)

"""
function bcd(param,w;wtol::Float64=1e-3, maxIter::Int=10, maxIterQP::Int=5, storeInterm = 0, out::Int=1, pwr::Float64=param.pwr)

	n = length(w)
    if storeInterm > 1
        storeInterm = 0
    end

	if storeInterm == 1
		ws      = zeros(n,maxIter+1)
		ws[:,1] = w
		Ms      = falses(size(param.d,1),n,maxIter)
	end

	his = zeros(maxIter,3)

	if out>1
        @printf "=== Sequential Quadratic Program (n=%d) ===\n" n
        @printf "%4s\t%3s\t\t%6s\t\t%6s\n" "iter" "Jc" "|dM|" "|dW|"
	end

	iter =1; flag = -1
	wOld = copy(w); M = falses(size(param.d,1),n); MOld = copy(M)
	while iter < maxIter

		his[iter,1] = misfit(param,M,w,pwr=pwr)[1]

		M     = findM(param,w,pwr=pwr)

		his[iter,2] = sum(abs.(M-MOld))
		if his[iter,2]==0
			flag=0
			# break
		end

		w =  findw(param,M,maxIterQP=maxIterQP,pwr=pwr)
		# evaluate objective function
		his[iter,3] = norm(w-wOld);

		(out>1) &&	@printf "%d\t%1.2e\t%d\t%1.2e\n" iter his[iter,1] his[iter,2] his[iter,3]

		if his[iter,3] < wtol
			flag = 0
			break
		end

		if storeInterm == 1
			ws[:,iter+1] = w;
			Ms[:,:,iter] = M
		end
		wOld = copy(w)
        MOld = copy(M)
		iter +=1
	end

    # need to fix misleading print statement
	if out>=0
        if flag==-1
            @printf "bcd iterated maxIter(=%d) times but reached only |w - wOld|=%1.2e (wtol=%1.2e)\n" maxIter his[iter,3] wtol
        elseif flag==0 && out>=1
            @printf "bcd achieved |w - wOld| = %1.2e < %1.2e = wtol at iteration %d.\n" his[iter,3] wtol iter
        end
    end

	if storeInterm == 1
		return Ms[:,:,1:iter],ws[:,1:iter+1], his[1:iter,:]
	else
		return M,w, his[1:iter,:]
	end
end


text:
function variable_converter_filter_voltage(pm::_PM.AbstractWModels; kwargs...)
    variable_converter_filter_voltage_magnitude_sqr(pm; kwargs...)
    variable_converter_filter_voltage_cross_products(pm; kwargs...)
end

function variable_converter_internal_voltage(pm::_PM.AbstractWModels; kwargs...)
    variable_converter_internal_voltage_magnitude_sqr(pm; kwargs...)
    variable_converter_internal_voltage_cross_products(pm; kwargs...)
end
"""
Creates lossy converter model between AC and DC side
```
pconv_ac[i] + pconv_dc[i] == a + b*I + c*Isq
```
"""
function constraint_converter_losses(pm::_PM.AbstractWModels, n::Int, i::Int, a, b, c, plmax)
    pconv_ac = _PM.var(pm, n,  :pconv_ac, i)
    pconv_dc = _PM.var(pm, n,  :pconv_dc, i)
    iconv = _PM.var(pm, n,  :iconv_ac, i)
    iconv_sq = _PM.var(pm, n,  :iconv_ac_sq, i)

    JuMP.@constraint(pm.model, pconv_ac + pconv_dc == a + b*iconv + c*iconv_sq)
end
"""
Converter transformer constraints
```
p_tf_fr ==  g/(tm^2)*w_fr + -g/(tm)*wr + -b/(tm)*wi)
q_tf_fr == -b/(tm^2)*w_fr +  b/(tm)*wr + -g/(tm)*wi)
p_tf_to ==  g*w_to + -g/(tm)*wr     + -b/(tm)*(-wi))
q_tf_to == -b*w_to +  b/(tm)*wr     + -g/(tm)*(-wi))
```
"""
function constraint_conv_transformer(pm::_PM.AbstractWRModels, n::Int, i::Int, rtf, xtf, acbus, tm, transformer)
    ptf_fr = _PM.var(pm, n,  :pconv_tf_fr, i)
    qtf_fr = _PM.var(pm, n,  :qconv_tf_fr, i)
    ptf_to = _PM.var(pm, n,  :pconv_tf_to, i)
    qtf_to = _PM.var(pm, n,  :qconv_tf_to, i)

    w = _PM.var(pm, n,  :w, acbus)  # vm^2
    wf = _PM.var(pm, n,  :wf_ac, i)   # vmf * vmf
    wrf = _PM.var(pm, n,  :wrf_ac, i) # vm*vmf*cos(va-vaf) =  vmf*vm*cos(vaf-va)
    wif = _PM.var(pm, n,  :wif_ac, i) # vm*vmf*sin(va-vaf) = -vmf*vm*sin(vaf-va)

    ztf = rtf + im*xtf
    if transformer
        ytf = 1/(rtf + im*xtf)
        gtf = real(ytf)
        btf = imag(ytf)
        c1, c2, c3, c4 = ac_power_flow_constraints_w(pm, gtf, btf, w, wf, wrf, wif, ptf_fr, ptf_to, qtf_fr, qtf_to, tm)
        constraint_voltage_product_converter(pm, wrf, wif, wf, w)
    else
        pcon, qcon = constraint_lossless_section(pm, w, wf, wrf, wif, ptf_fr, ptf_to, qtf_fr, qtf_to)
    end
end
"constraints for a voltage magnitude transformer + series impedance"
function ac_power_flow_constraints_w(pm::_PM.AbstractWRModels, g, b, w_fr, w_to, wr, wi, p_fr, p_to, q_fr, q_to, tm)
    c1 = JuMP.@constraint(pm.model, p_fr ==  g/(tm^2)*w_fr + -g/(tm)*wr + -b/(tm)*wi)
    c2 = JuMP.@constraint(pm.model, q_fr == -b/(tm^2)*w_fr +  b/(tm)*wr + -g/(tm)*wi)
    c3 = JuMP.@constraint(pm.model, p_to ==  g*w_to + -g/(tm)*wr     + -b/(tm)*(-wi))
    c4 = JuMP.@constraint(pm.model, q_to == -b*w_to +  b/(tm)*wr     + -g/(tm)*(-wi))
    return c1, c2, c3, c4
end
"""
Converter reactor constraints
```
p_pr_fr ==  g/(tm^2)*w_fr + -g/(tm)*wr + -b/(tm)*wi)
q_pr_fr == -b/(tm^2)*w_fr +  b/(tm)*wr + -g/(tm)*wi)
p_pr_to ==  g*w_to + -g/(tm)*wr     + -b/(tm)*(-wi))
q_pr_to == -b*w_to +  b/(tm)*wr     + -g/(tm)*(-wi))
```
"""
function constraint_conv_reactor(pm::_PM.AbstractWRModels, n::Int, i::Int, rc, xc, reactor)
    pconv_ac = _PM.var(pm, n,  :pconv_ac, i)
    qconv_ac = _PM.var(pm, n,  :qconv_ac, i)
    ppr_to = - pconv_ac
    qpr_to = - qconv_ac
    ppr_fr = _PM.var(pm, n,  :pconv_pr_fr, i)
    qpr_fr = _PM.var(pm, n,  :qconv_pr_fr, i)

    wf = _PM.var(pm, n,  :wf_ac, i)
    wc = _PM.var(pm, n,  :wc_ac, i)
    wrc = _PM.var(pm, n,  :wrc_ac, i)
    wic = _PM.var(pm, n,  :wic_ac, i)

    zc  = rc  + im*xc

    if reactor
        yc = 1/(zc)
        gc = real(yc)
        bc = imag(yc)
        constraint_voltage_product_converter(pm, wrc, wic, wf, wc)
        c1, c2, c3, c4 = ac_power_flow_constraints_w(pm, gc, bc, wf, wc, wrc, wic, ppr_fr, ppr_to, qpr_fr, qpr_to, 1)
    else
        pcon, qcon = constraint_lossless_section(pm, wf, wc, wrc, wic, ppr_fr, ppr_to, qpr_fr, qpr_to)
    end
end
"""
Converter filter constraints

```
p_pr_fr + p_tf_to == 0
q_pr_fr + q_tf_to + -bv*filter*wf == 0
```
"""
function constraint_conv_filter(pm::_PM.AbstractWModels, n::Int, i::Int, bv, filter)
    ppr_fr = _PM.var(pm, n,  :pconv_pr_fr, i)
    qpr_fr = _PM.var(pm, n,  :qconv_pr_fr, i)
    ptf_to = _PM.var(pm, n,  :pconv_tf_to, i)
    qtf_to = _PM.var(pm, n,  :qconv_tf_to, i)

    wf = _PM.var(pm, n,  :wf_ac, i)

    JuMP.@constraint(pm.model, ppr_fr + ptf_to == 0 )
    JuMP.@constraint(pm.model, qpr_fr + qtf_to + -bv*filter*wf == 0)
end

function constraint_lossless_section(pm::_PM.AbstractWModels, w_fr, w_to, wr, wi, p_fr, p_to, q_fr, q_to)
    JuMP.@constraint(pm.model, w_fr ==  w_to)
    JuMP.@constraint(pm.model, wr   ==  w_fr)
    JuMP.@constraint(pm.model, wi   ==  0)

    pcon = JuMP.@constraint(pm.model, p_fr + p_to == 0)
    qcon = JuMP.@constraint(pm.model, q_fr + q_to == 0)
    return pcon, qcon
end

function add_converter_voltage_setpoint(sol, pm::_PM.AbstractWModels)
    _PM.add_setpoint!(sol, pm, "convdc", "vmconv", :wc_ac; scale = (x,item) -> sqrt(x))
    _PM.add_setpoint!(sol, pm, "convdc", "vmfilt", :wf_ac; scale = (x,item) -> sqrt(x))
end
"""
LCC firing angle constraints
```
qconv_ac >= Q1 + (pconv_ac-P1) * (Q2-Q1)/(P2-P1)

P1 = cos(0) * Srated
Q1 = sin(0) * Srated
P2 = cos(pi) * Srated
Q2 = sin(pi) * Srated
```
"""
function constraint_conv_firing_angle(pm::_PM.AbstractWModels, n::Int, i::Int, S, P1, Q1, P2, Q2)
    pc = _PM.var(pm, n,  :pconv_ac, i)
    qc = _PM.var(pm, n,  :qconv_ac, i)
    coeff = (Q2-Q1)/(P2-P1)

    JuMP.@constraint(pm.model, qc >= Q1 + (pc-P1) * coeff )
end

################ TNEP constraints #####################
function variable_converter_filter_voltage_ne(pm::_PM.AbstractWRModels; kwargs...)
    variable_converter_filter_voltage_magnitude_sqr_ne(pm; kwargs...)
    variable_converter_filter_voltage_cross_products_ne(pm; kwargs...)
end

function variable_converter_internal_voltage_ne(pm::_PM.AbstractWRModels; kwargs...)
    variable_converter_internal_voltage_magnitude_sqr_ne(pm; kwargs...)
    variable_converter_internal_voltage_cross_products_ne(pm; kwargs...)
end

"""
Creates lossy converter model between AC and DC side

```
pconv_ac[i] + pconv_dc[i] == a*z + b*I + c*Isq
```
"""
function constraint_converter_losses_ne(pm::_PM.AbstractWModels, n::Int, i::Int, a, b, c, plmax)
    pconv_ac = _PM.var(pm, n, :pconv_ac_ne, i)
    pconv_dc = _PM.var(pm, n, :pconv_dc_ne, i)
    iconv = _PM.var(pm, n, :iconv_ac_ne, i)
    iconv_sq = _PM.var(pm, n, :iconv_ac_sq_ne, i)
    z = _PM.var(pm, n, :conv_ne, i)

    JuMP.@constraint(pm.model, pconv_ac + pconv_dc == a*z + b*iconv + c*iconv_sq)
end

"""
Converter transformer constraints

```
p_tf_fr ==  g/(tm^2)*w_fr + -g/(tm)*wr + -b/(tm)*wi)
q_tf_fr == -b/(tm^2)*w_fr +  b/(tm)*wr + -g/(tm)*wi)
p_tf_to ==  g*w_to + -g/(tm)*wr     + -b/(tm)*(-wi))
q_tf_to == -b*w_to +  b/(tm)*wr     + -g/(tm)*(-wi))
```
"""
function constraint_conv_transformer_ne(pm::_PM.AbstractWRModels, n::Int, i::Int, rtf, xtf, acbus, tm, transformer)
    ptf_fr = _PM.var(pm, n, :pconv_tf_fr_ne, i)
    qtf_fr = _PM.var(pm, n, :qconv_tf_fr_ne, i)
    ptf_to = _PM.var(pm, n, :pconv_tf_to_ne, i)
    qtf_to = _PM.var(pm, n, :qconv_tf_to_ne, i)

    w = _PM.var(pm, n, :w, acbus)  # vm^2
    w_du = _PM.var(pm, n, :w_du, i)
    wf = _PM.var(pm, n, :wf_ac_ne, i)   # vmf * vmf
    wrf = _PM.var(pm, n, :wrf_ac_ne, i) # vm*vmf*cos(va-vaf) =  vmf*vm*cos(vaf-va)
    wif = _PM.var(pm, n, :wif_ac_ne, i) # vm*vmf*sin(va-vaf) = -vmf*vm*sin(vaf-va)
    z = _PM.var(pm, n, :conv_ne)[i]
    ztf = rtf + im*xtf
    if transformer
        ytf = 1/(rtf + im*xtf)
        gtf = real(ytf)
        btf = imag(ytf)
        c1, c2, c3, c4 = ac_power_flow_constraints_w(pm, gtf, btf, w_du, wf, wrf, wif, ptf_fr, ptf_to, qtf_fr, qtf_to, tm)

        _IM.relaxation_equality_on_off(pm.model, w, w_du, z)
        JuMP.@constraint(pm.model, w_du >= z*JuMP.lower_bound(w))
        JuMP.@constraint(pm.model, w_du <= z*JuMP.upper_bound(w))
        constraint_voltage_product_converter_ne(pm, wrf, wif, wf, w, z)
    else
        pcon, qcon = constraint_lossless_section_ne(pm, w_du, wf, wrf, wif, ptf_fr, ptf_to, qtf_fr, qtf_to)
        _IM.relaxation_equality_on_off(pm.model, w, w_du, z)
        JuMP.@constraint(pm.model, w_du >= z*JuMP.lower_bound(w))
        JuMP.@constraint(pm.model, w_du <= z*JuMP.upper_bound(w))
    end
end

"""
Converter reactor constraints

```
p_pr_fr ==  g/(tm^2)*w_fr + -g/(tm)*wr + -b/(tm)*wi)
q_pr_fr == -b/(tm^2)*w_fr +  b/(tm)*wr + -g/(tm)*wi)
p_pr_to ==  g*w_to + -g/(tm)*wr     + -b/(tm)*(-wi))
q_pr_to == -b*w_to +  b/(tm)*wr     + -g/(tm)*(-wi))
```
"""
function constraint_conv_reactor_ne(pm::_PM.AbstractWRModels, n::Int, i::Int, rc, xc, reactor)
    pconv_ac = _PM.var(pm, n, :pconv_ac_ne, i)
    qconv_ac = _PM.var(pm, n, :qconv_ac_ne, i)
    ppr_to = - pconv_ac
    qpr_to = - qconv_ac
    ppr_fr = _PM.var(pm, n, :pconv_pr_fr_ne, i)
    qpr_fr = _PM.var(pm, n, :qconv_pr_fr_ne, i)

    wf = _PM.var(pm, n, :wf_ac_ne, i)
    wc = _PM.var(pm, n, :wc_ac_ne, i)
    wrc = _PM.var(pm, n, :wrc_ac_ne, i)
    wic = _PM.var(pm, n, :wic_ac_ne, i)
    z = _PM.var(pm, n, :conv_ne)[i]
    zc  = rc  + im*xc

    if reactor
        yc = 1/(zc)
        gc = real(yc)
        bc = imag(yc)
        c1, c2, c3, c4 = ac_power_flow_constraints_w(pm, gc, bc, wf, wc, wrc, wic, ppr_fr, ppr_to, qpr_fr, qpr_to, 1)
        constraint_voltage_product_converter_ne(pm, wrc, wic, wf, wc, z)
    else
        pcon, qcon = constraint_lossless_section_ne(pm, wf, wc, wrc, wic, ppr_fr, ppr_to, qpr_fr, qpr_to)
    end
end

"""
Converter filter constraints

```
p_pr_fr + p_tf_to == 0
q_pr_fr + q_tf_to + -bv*filter*wf == 0
```
"""
function constraint_conv_filter_ne(pm::_PM.AbstractWModels, n::Int, i::Int, bv, filter)
    ppr_fr = _PM.var(pm, n, :pconv_pr_fr_ne, i)
    qpr_fr = _PM.var(pm, n, :qconv_pr_fr_ne, i)
    ptf_to = _PM.var(pm, n, :pconv_tf_to_ne, i)
    qtf_to = _PM.var(pm, n, :qconv_tf_to_ne, i)
    wf = _PM.var(pm, n, :wf_ac_ne, i)

    JuMP.@constraint(pm.model, ppr_fr + ptf_to == 0 )
    JuMP.@constraint(pm.model, qpr_fr + qtf_to + (-bv) * filter * (wf) == 0)
end

function constraint_lossless_section_ne(pm::_PM.AbstractWModels, w_fr, w_to, wr, wi, p_fr, p_to, q_fr, q_to)
    JuMP.@constraint(pm.model, w_fr ==  w_to)
    JuMP.@constraint(pm.model, wr   ==  w_fr)
    JuMP.@constraint(pm.model, wi   ==  0)

    pcon = JuMP.@constraint(pm.model, p_fr + p_to == 0)
    qcon = JuMP.@constraint(pm.model, q_fr + q_to == 0)
    return pcon, qcon
end
function add_converter_voltage_setpoint_ne(sol, pm::_PM.AbstractWModels)
    _PM.add_setpoint!(sol, pm, "convdc", "vmconv", :wc_ac_ne; scale = (x,item) -> sqrt(x))
    _PM.add_setpoint!(sol, pm, "convdc", "vmfilt", :wf_ac_ne; scale = (x,item) -> sqrt(x))
end

"""
LCC firing angle constraints

```
qconv_ac >= Q1 + (pconv_ac-P1) * (Q2-Q1)/(P2-P1)

P1 = cos(0) * Srated
Q1 = sin(0) * Srated
P2 = cos(pi) * Srated
Q2 = sin(pi) * Srated
```
"""
function constraint_conv_firing_angle_ne(pm::_PM.AbstractWModels, n::Int, i::Int, S, P1, Q1, P2, Q2)
    pc = _PM.var(pm, n, :pconv_ac_ne, i)
    qc = _PM.var(pm, n, :qconv_ac_ne, i)
    coeff = (Q2-Q1)/(P2-P1)

    JuMP.@constraint(pm.model, qc >= Q1 + (pc-P1) * coeff )
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -23.072345451695163 <= q <= 18.38163265505858)

                     add_NL_constraint(m, :(sigmoid(-0.6895881478464707 + 0.42824131807966515*sigmoid(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4])) + -0.8758217593243995*sigmoid(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4])) + -0.3819126296826152*sigmoid(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4])) + -0.7523924735752416*sigmoid(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))) + sigmoid(0.5903696565764314 + 0.6045010784464115*sigmoid(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4])) + 0.19992590084675577*sigmoid(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4])) + -0.13493555263327828*sigmoid(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4])) + 0.0007730145672497635*sigmoid(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))) + sigmoid(-0.987583954230804 + 0.6297203029880079*sigmoid(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4])) + -0.939181196020844*sigmoid(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4])) + 0.23747324715066886*sigmoid(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4])) + 0.2556344288953589*sigmoid(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))) + sigmoid(-0.44105426128556546 + -0.39452594471110514*sigmoid(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4])) + 0.16620534874676274*sigmoid(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4])) + 0.6245292013722641*sigmoid(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4])) + 0.9983459126930239*sigmoid(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module Forces

  initial_point(side,N) =
     [ [-side/2 + rand()*side, -side/2+rand()*side] for i in 1:N ]

  # Energies
  upair(x,y) = (x[1]-y[1])^2 + (x[2]-y[2])^2
  function utotal(p)
    u = 0.
    for i in 1:length(p)-1
      for j in i+1:length(p)
        u += upair(p[i],p[j])
      end
    end
    u
  end

  # Forces
  function forcepair(x,y)
    dx = x[1]-y[1]
    dy = x[2]-y[2]

    upair = dx^2 + dy^2

    dudx1 = -2*dx
    dudx2 = -2*dy
    return upair, (dudx1, dudx2)
  end

  # Compute all forces, modifies vector f
  function forces!(p,f)
    u = 0.
    for i in 1:length(f)
      f[i] .= 0
    end
    for i in 1:length(p)-1
      for j in i+1:length(p)
        up, fp = forcepair(p[i],p[j])
        u += up
        f[i] .= f[i] .+ fp
        f[j] .= f[j] .- fp
      end
    end
    return u, f
  end

  export upair, utotal, forcepair, forces!, initial_point
end


text:
using Documenter
using MyPlayground

makedocs(
    sitename = "MyPlayground",
    modules = [MyPlayground],
    strict = true,
)

deploydocs(
    repo = "github.com/tkf/MyPlayground.jl",
    push_preview = true,
)


text:
# This file is part of RLinearAlgebra.jl

# using LinearAlgebra

"""
    LSLogFull <: LinSysSolverLog

A mutable structure that stores information about a randomized linear solver's behavior.
    The log assumes that the full linear system is available for processing. The goal of
    this log is usually for research, development or testing as it is unlikely that the
    entire residual vector is readily available.

# Fields
- `collection_rate::Int64`, the frequency with which to record information to append to the
    remaining fields, starting with the initialization (i.e., iteration 0).
- `resid_hist::Vector{Float64}`, retains a vector of numbers corresponding to the residual
    (uses the whole system to compute the residual)
- `resid_norm::Function`, a function that accepts a single vector argument and returns a
    scalar. Used to compute the residual size.
- `iterations::Int64`, the number of iterations of the solver.
- `converged::Bool`, a flag to indicate whether the system has converged by some measure

# Constructors
- Calling `LSLogFull()` sets `collection_rate = 1`, `resid_hist = Float64[]`,
    `resid_norm = norm` (Euclidean norm), `iterations = -1`, and `converged = false`.
- Calling `LSLogFull(cr::Int64)` is the same as calling `LSLogFull()` except
    `collection_rate = cr`.
"""
mutable struct LSLogFull <: LinSysSolverLog
    collection_rate::Int64
    resid_hist::Vector{Float64}
    resid_norm::Function
    iterations::Int64
    converged::Bool
end
LSLogFull() = LSLogFull(1, Float64[], norm, -1, false)
LSLogFull(cr::Int64) = LSLogFull(cr, Float64[], norm, -1, false)

# Common interface for update
function log_update!(
    log::LSLogFull,
    sampler::LinSysSampler,
    x::AbstractVector,
    samp::Tuple,
    iter::Int64,
    A::AbstractArray,
    b::AbstractVector
)
    # Update iteration counter
    log.iterations = iter

    # Check whether to record information
    if mod(iter, 1:log.collection_rate) == log.collection_rate

        # Push residual norm of current iterate
        push!(log.resid_hist, log.resid_norm(A * x - b))
    end

    return nothing
end


text:
using ClimateMachine.DGMethods.NumericalFluxes: NumericalFluxFirstOrder
import ClimateMachine.DGMethods.NumericalFluxes:
    numerical_volume_conservative_flux_first_order!,
    numerical_volume_fluctuation_flux_first_order!,
    ave,
    numerical_flux_first_order!,
    numerical_flux_second_order!,
    numerical_boundary_flux_second_order!
import ClimateMachine.BalanceLaws: wavespeed

struct CentralVolumeFlux <: NumericalFluxFirstOrder end
struct KGVolumeFlux <: NumericalFluxFirstOrder end
struct LinearKGVolumeFlux <: NumericalFluxFirstOrder end
struct VeryLinearKGVolumeFlux <: NumericalFluxFirstOrder end

function numerical_volume_fluctuation_flux_first_order!(
    ::NumericalFluxFirstOrder,
    model::DryAtmosModel,
    D::Grad,
    state_1::Vars,
    aux_1::Vars,
    state_2::Vars,
    aux_2::Vars,
)

    sources = model.physics.sources

    ntuple(Val(length(sources))) do s
        Base.@_inline_meta
        calc_fluctuation_component!(D, sources[s], state_1, state_2, aux_1, aux_2)
    end
    #=
    if fluctuation_gravity
        ρ_1, ρ_2 = state_1.ρ, state_2.ρ
        Φ_1, Φ_2 = aux_1.Φ, aux_2.Φ

        α = ave(ρ_1, ρ_2) * 0.5

        D.ρu -= α * (Φ_1 - Φ_2) * I
    end
    =#
end

function numerical_volume_conservative_flux_first_order!(
    ::CentralVolumeFlux,
    m::DryAtmosModel,
    F::Grad,
    state_1::Vars,
    aux_1::Vars,
    state_2::Vars,
    aux_2::Vars,
)
    FT = eltype(F)
    F_1 = similar(F)
    flux_first_order!(m, F_1, state_1, aux_1, FT(0), EveryDirection())

    F_2 = similar(F)
    flux_first_order!(m, F_2, state_2, aux_2, FT(0), EveryDirection())

    parent(F) .= (parent(F_1) .+ parent(F_2)) ./ 2
end

function numerical_volume_conservative_flux_first_order!(
    ::KGVolumeFlux,
    model::DryAtmosModel,
    F::Grad,
    state_1::Vars,
    aux_1::Vars,
    state_2::Vars,
    aux_2::Vars,
)
    eos = model.physics.eos
    parameters = model.physics.parameters

    ρ_1 = state_1.ρ
    ρu_1 = state_1.ρu
    ρe_1 = state_1.ρe
    ρq_1 = state_1.ρq
    u_1 = ρu_1 / ρ_1
    e_1 = ρe_1 / ρ_1
    q_1 = ρq_1 / ρ_1
    p_1 = calc_pressure(eos, state_1, aux_1, parameters)

    ρ_2 = state_2.ρ
    ρu_2 = state_2.ρu
    ρe_2 = state_2.ρe
    ρq_2 = state_2.ρq
    u_2 = ρu_2 / ρ_2
    e_2 = ρe_2 / ρ_2
    q_2 = ρq_2 / ρ_2
    p_2 = calc_pressure(eos, state_2, aux_2, parameters)

    ρ_avg = ave(ρ_1, ρ_2)
    u_avg = ave(u_1, u_2)
    e_avg = ave(e_1, e_2)
    q_avg = ave(q_1, q_2)
    p_avg = ave(p_1, p_2)

    F.ρ = ρ_avg * u_avg
    F.ρu = p_avg * I + ρ_avg * u_avg .* u_avg'
    F.ρe = ρ_avg * u_avg * e_avg + p_avg * u_avg
    F.ρq = ρ_avg * u_avg * q_avg
end

function numerical_volume_conservative_flux_first_order!(
    ::LinearKGVolumeFlux,
    model::DryAtmosModel,
    F::Grad,
    state_1::Vars,
    aux_1::Vars,
    state_2::Vars,
    aux_2::Vars,
)
    eos = model.physics.eos
    parameters = model.physics.parameters

    ρu_1 = state_1.ρu
    ρuᵣ = ρu_1 * 0
    p_1 = calc_linear_pressure(eos, state_1, aux_1, parameters)

    # grab reference state
    ρᵣ_1 = aux_1.ref_state.ρ
    pᵣ_1 = aux_1.ref_state.p
    ρeᵣ_1 = aux_1.ref_state.ρe

    # only ρu fluctuates in the non-pressure terms
    u_1 = ρu_1 / ρᵣ_1
    eᵣ_1 = ρeᵣ_1 / ρᵣ_1
    ρu_2 = state_2.ρu

    ρuᵣ = ρu_2 * 0
    p_2 = calc_linear_pressure(eos, state_2, aux_2, parameters)

    # grab reference state
    ρᵣ_2 = aux_2.ref_state.ρ
    pᵣ_2 = aux_2.ref_state.p
    ρeᵣ_2 = aux_2.ref_state.ρe

    # only ρu fluctuates in the non-pressure terms
    u_2 = ρu_2 / ρᵣ_2
    eᵣ_2 = ρeᵣ_2 / ρᵣ_2

    # construct averages
    ρᵣ_avg = ave(ρᵣ_1, ρᵣ_2)
    eᵣ_avg = ave(eᵣ_1, eᵣ_2)
    pᵣ_avg = ave(pᵣ_1, pᵣ_2)

    u_avg = ave(u_1, u_2)
    p_avg = ave(p_1, p_2)

    F.ρ = ρᵣ_avg * u_avg
    F.ρu = p_avg * I + ρuᵣ .* ρuᵣ' # the latter term is needed to determine size of I
    F.ρe = (ρᵣ_avg * eᵣ_avg + pᵣ_avg) * u_avg
end

# TODO: UPDATE EOS
#=
function calc_linear_pressure(eos::DryIdealGas{(:ρ, :ρu, :ρe)}, state, aux, params)
    ρ  = state.ρ
    ρe = state.ρe
    Φ  = aux.Φ
    γ  = calc_γ(eos, state, params)
    # full : (ρe - dot(ρu, ρu) / 2ρ - ρ * Φ)

    return (γ - 1) * (ρe - ρ * Φ - dot(ρu, ρuᵣ) / ρᵣ + ρ * dot(ρuᵣ, ρuᵣ) / (2 * ρᵣ^2) )
end
=#
function numerical_volume_conservative_flux_first_order!(
    ::VeryLinearKGVolumeFlux,
    model::DryAtmosModel,
    F::Grad,
    state_1::Vars,
    aux_1::Vars,
    state_2::Vars,
    aux_2::Vars,
)
    eos = model.physics.eos
    parameters = model.physics.parameters

    ## State 1 Stuff
    # unpack the perturbation state
    ρ_1 = state_1.ρ
    ρu_1 = state_1.ρu
    ρe_1 = state_1.ρe
    ρq_1 = state_1.ρq

    # grab reference state
    ρᵣ_1 = aux_1.ref_state.ρ
    ρuᵣ_1 = aux_1.ref_state.ρu
    ρeᵣ_1 = aux_1.ref_state.ρe
    ρqᵣ_1 = aux_1.ref_state.ρq
    pᵣ_1 = aux_1.ref_state.p

    # calculate pressure perturbation
    p_1 = calc_very_linear_pressure(eos, state_1, aux_1, parameters)

    # calculate u_1, q_1, e_1, and reference states
    u_1 = ρu_1 / ρᵣ_1 - ρ_1 * ρuᵣ_1 / (ρᵣ_1^2)
    q_1 = ρq_1 / ρᵣ_1 - ρ_1 * ρqᵣ_1 / (ρᵣ_1^2)
    e_1 = ρe_1 / ρᵣ_1 - ρ_1 * ρeᵣ_1 / (ρᵣ_1^2)

    uᵣ_1 = ρuᵣ_1 / ρᵣ_1
    qᵣ_1 = ρqᵣ_1 / ρᵣ_1
    eᵣ_1 = ρeᵣ_1 / ρᵣ_1

    ## State 2 Stuff
    # unpack the state perubation
    ρ_2 = state_2.ρ
    ρu_2 = state_2.ρu
    ρe_2 = state_2.ρe
    ρq_2 = state_2.ρq

    # grab reference state
    ρᵣ_2 = aux_2.ref_state.ρ
    ρuᵣ_2 = aux_2.ref_state.ρu
    ρeᵣ_2 = aux_2.ref_state.ρe
    ρqᵣ_2 = aux_2.ref_state.ρq
    pᵣ_2 = aux_2.ref_state.p

    # calculate pressure perturbation
    p_2 = calc_very_linear_pressure(eos, state_2, aux_2, parameters)

    # calculate u_2, q_2, e_2, and reference states
    u_2 = ρu_2 / ρᵣ_2 - ρ_2 * ρuᵣ_2 / (ρᵣ_2^2)
    q_2 = ρq_2 / ρᵣ_2 - ρ_2 * ρqᵣ_2 / (ρᵣ_2^2)
    e_2 = ρe_2 / ρᵣ_2 - ρ_2 * ρeᵣ_2 / (ρᵣ_2^2)

    uᵣ_2 = ρuᵣ_2 / ρᵣ_2
    qᵣ_2 = ρqᵣ_2 / ρᵣ_2
    eᵣ_2 = ρeᵣ_2 / ρᵣ_2

    # construct averages for perturbation variables
    ρ_avg = ave(ρ_1, ρ_2)
    u_avg = ave(u_1, u_2)
    e_avg = ave(e_1, e_2)
    q_avg = ave(q_1, q_2)
    p_avg = ave(p_1, p_2)

    # construct averages for reference variables
    ρᵣ_avg = ave(ρᵣ_1, ρᵣ_2)
    uᵣ_avg = ave(uᵣ_1, uᵣ_2)
    eᵣ_avg = ave(eᵣ_1, eᵣ_2)
    qᵣ_avg = ave(qᵣ_1, qᵣ_2)
    pᵣ_avg = ave(pᵣ_1, pᵣ_2)

    F.ρ = ρᵣ_avg * u_avg + ρ_avg * uᵣ_avg
    F.ρu = p_avg * I + ρᵣ_avg .* (uᵣ_avg .* u_avg' + u_avg .* uᵣ_avg')
    F.ρu += (ρ_avg .* uᵣ_avg) .* uᵣ_avg'
    F.ρe = (ρᵣ_avg * eᵣ_avg + pᵣ_avg) * u_avg
    F.ρe += (ρᵣ_avg * e_avg + ρ_avg * eᵣ_avg + p_avg) * uᵣ_avg
    F.ρq = ρᵣ_avg * qᵣ_avg * u_avg + (ρᵣ_avg * q_avg + ρ_avg * qᵣ_avg) * uᵣ_avg
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
#=
Example run from command line reading from the parameter file  examples/ia_example1.jl:
$  julia run.jl examples/ia_example1

This run produces the output CSV file:  examples/ia_example1.csv
Front end for src/inf_alleles.jl
This file contains functions for the infinite alleles model in support of the Nearly Neutral paper with title ???
  by Cheyenne Laue, ???, and Alden Wright
Author:  Alden H. Wright, Department of Computer Science, Univeristy of Montana, Missoula, MT 59812 USA
    alden.wright@umontana.edu
=#
using Dates
using Statistics

@doc """ function run_trials()
  Run multiple trials of the infinite alleles model.
  popsize_multiplier_list  is a list of ratios of N to sample size n
          For the nearly neutral paper, sampling is not used, so it is always the case that n==N.
  There are 2 possibilities for iterating over mutation rate.
   (1) mu_list_flag==true:    Then  mu_list  is a list of mutation rates.
   (2) mu_list_flag==false:   Then  N_mu_list  is a list of N_mu values.  N_mu is the product of N and mu.  This corresponds to constant theta.
"""
function ia_run_trials( mu_list_flag::Bool=false)
  #println("stream: ",stream)
  popsize_multiplier_list = [1]
  trial = 1
  N = N_list[1]
  n = Int(floor(N*(1//popsize_multiplier_list[1])))
  if !mu_list_flag
    tr = InfAlleles.trial_result( nn_simtype, n, N, N_mu_list[1], ngens, burn_in, dfe, dfe_str, use_poplist=use_poplist )
    ia_writeheader(Base.stdout, popsize_multiplier_list, N_list, N_mu_list, tr )
    ia_writeheader(stream, popsize_multiplier_list, N_list, N_mu_list, tr )
    for N_mu in N_mu_list
      for N in N_list
        for psize_m in popsize_multiplier_list
          n = Int(floor(N*(1//psize_m)))
          tr = InfAlleles.trial_result( nn_simtype, n, N, N_mu, ngens, burn_in, dfe, dfe_str, use_poplist=use_poplist )
          ia_run_trial( tr )
          ia_writerow(Base.stdout, trial, tr )
          ia_writerow(stream, trial, tr )
          trial += 1
          Base.flush(Base.stdout)
          Base.flush(stream)
        end
      end
    end
  else  # if mu_list_flag
    tr = InfAlleles.trial_result( nn_simtype, n, N, N*mu_list[1], ngens, burn_in, dfe, dfe_str, use_poplist=use_poplist )
    ia_writeheader(stream, popsize_multiplier_list, N_list, mu_list, tr, mu_list_flag=mu_list_flag )
    for mu in mu_list
      for N in N_list
        for psize_m in popsize_multiplier_list
          n = Int(floor(N*(1//psize_m)))
          N_mu = N*mu
          tr = InfAlleles.trial_result( nn_simtype, n, N, N_mu, ngens, burn_in, dfe, dfe_str, use_poplist=use_poplist )
          ia_run_trial( tr )
          ia_writerow(Base.stdout, trial, tr, mu_list_flag=mu_list_flag )
          ia_writerow(stream, trial, tr, mu_list_flag=mu_list_flag )
          trial += 1
          Base.flush(Base.stdout)
          Base.flush(stream)
        end
      end
    end
  end  # if !mu_list_flag
end

function ia_run_trial( tr::InfAlleles.trial_result )
  if tr.nn_simtype == 1
    if tr.use_poplist
      poplist = inf_alleles(tr,combine=false)
      if tr.n < tr.N
        poplist = map(x->sample_population(x,tr.n),poplist)
      end
      add_stats_to_trial_result!( tr, poplist )
    else
      inf_alleles(tr,combine=false)
    end
    #print_trial_result( tr )
    Base.flush(Base.stdout)
    return tr
  else
    println("nn_simtype ",nn_simtype," not implemented")
  end
end

function stddev(lst)
  N = length(lst)
  sqrt(mapreduce(x->x^2,+,lst)/(N-1) - sum(lst)^2/N/(N-1))
end

function add_stats_to_trial_result!( tr::InfAlleles.trial_result, poplist::Vector{Population} )
  pcounts = map(pop_counts64,poplist)
  richness_list = map(x->length(x),pcounts)
  #@assert isapprox( tr.average_richness, Statistics.mean(richness_list) )
  if !isapprox( tr.average_richness, Statistics.mean(richness_list) )
    println("warning: tr.average_richness not approximately Statistics.mean(richness_list)" )
  end
  tr.average_richness = Statistics.mean(richness_list)
  #@assert isapprox( tr.stderr_richness, stddev(richness_list)/sqrt(length(richness_list) ) )
  if !isapprox( tr.stderr_richness, stddev(richness_list)/sqrt(length(richness_list) ) )
    println("warning: tr.stderr_richness not approximately stddev(richness_list)/sqrt(length(richness_list) )" )
  end
  tr.stderr_richness = stddev(richness_list)/sqrt(length(richness_list))
  w_homoz_list = map(watterson_homozygosity,pcounts)
  #@assert isapprox( tr.w_homoz, Statistics.mean(w_homoz_list ) )
  if !isapprox( tr.w_homoz, Statistics.mean(w_homoz_list) )
    println("warning: tr.w_homoz not approximately Statistics.mean(richness_list)" )
  end
  tr.w_homoz = Statistics.mean(w_homoz_list)
  tr.stderr_w_homoz = stddev(w_homoz_list)/sqrt(length(w_homoz_list) )
  IQV_list = map(IQV,pcounts)
  #@assert isapprox( tr.IQV, Statistics.mean(IQV_list) )
  if !isapprox( tr.IQV, Statistics.mean(IQV_list) )
    println("warning: tr.IQV not approximately Statistics.mean(IQV_list)" )
  end
  tr.IQV = Statistics.mean(IQV_list)
  #@assert isapprox( tr.stderr_IQV, stddev(IQV_list)/sqrt(length(IQV_list)) )
  if !isapprox( tr.stderr_IQV, stddev(IQV_list)/sqrt(length(IQV_list)) )
    println("warning: tr.w_homoz not approximately stddev(IQV_list)/sqrt(length(IQV_list))" )
  end
  tr.stderr_IQV = stddev(IQV_list)/sqrt(length(IQV_list))
  #tr.entropy = map(entropy,pcounts)
  entropy_list = map(InfAlleles.entropy,pcounts)
  tr.entropy = Statistics.mean(entropy_list)
  tr.stderr_entropy = stddev(entropy_list)/sqrt(length(entropy_list))
end

@doc """ function ia_writeheader()
Infinite alleles version.
Write header line for output CSV file that corresponds to stream.
See the comments for run_simulation for a description of the parameters.
"""
function ia_writeheader(stream::IO, popsize_multiplier_list::Vector{Int64}, N_list::Vector{Int64},
    mu_list::Vector{Float64}, tr::InfAlleles.trial_result; mu_list_flag::Bool=false )
  dfe_params = [:dfe_adv_prob, :dfe_adv_alpha, :dfe_adv_beta, :dfe_disadv_prob, :dfe_disadv_alpha, :dfe_disadv_beta]
  param_strings = [
    "# $(string(Dates.today()))",
    (@isdefined(seed) ? "# random number seed: $(seed)" : "# no random number seed defined"),
    "# $((nn_simtype==1) ? "infinite alleles model" : "infinite_sites_model")",
    "# popsize_multiplier_list=$(popsize_multiplier_list)",
    "# N_list=$(N_list)",
    mu_list_flag ? "# mu_list=$(mu_list)" : # N_mu_list=$(N_mu_list)",
    "# ngens=$(tr.ngens)",
    "# burn_in=$(tr.burn_in)",
    "# dfe=$(tr.dfe)",
    "# dfe_str=$(tr.dfe_str)"
    ]
  for dp in dfe_params
    if isdefined(Main,dp)
      Base.push!(param_strings,"# $(string(dp))=$(eval(dp))")
    end
  end
  write(stream, join(param_strings, "\n"), "\n")
  first_heads = [
    "type",
    #  "n",
    "N", (mu_list_flag ? "mu" : "N_mu")]
  last_heads =
  [
    "mean_fitness",
    "expected_richness",
    "average_richness",
    "stderr_richness",
    "expected_w_heteroz",
    "w_heteroz",
    "stderr_heteroz",
    "IQV",
    "stderr_IQV",
    "entropy",
    "stderr_entropy"
  ]
  line = join(vcat( first_heads, last_heads), ",")
  write(stream, line, "\n")
end

@doc """ function ia_writerow
Infinite alleles version.
Writes one row of the output CSV file.
"""
function ia_writerow(stream::IO, trial::Int64, tr::InfAlleles.trial_result; mu_list_flag::Bool=false   )
  first = Any[
    (isdefined(Main,:type_str) ? type_str : trial),
    #tr.n,           # sample size
    tr.N,           # popsize
    mu_list_flag ? tr.N_mu/tr.N : tr.N_mu,
  ]
  if tr.nn_simtype == 0
    mid = Any[]
  elseif tr.nn_simtype == 1
    mid = Any[
      tr.mean_fitness,
      tr.expected_richness,
      tr.average_richness,
      tr.stderr_richness,
      1.0-tr.expected_w_homoz,
      1.0-tr.w_homoz,
      tr.stderr_w_homoz,
      tr.IQV,
      tr.stderr_IQV,
      tr.entropy,
      tr.stderr_entropy
    ]
  end
  line = join( vcat( first, mid ), "," )
  write(stream, line, "\n")
end

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
#
#--------------------------------------#
export diver
#--------------------------------------#
"""
 for q ∈ H¹₀(Ω)  (test  pressure)
 for u⃗ ∈ H¹₀(Ω)² (trial velocity)

 (q,∇⋅u⃗) = (J*q)' * Bᵥ * [Dx Dy][ux]
                                [uy]

 Bᵥ is mass matrix on velocity grid
 J  is the interpolation matrix
 from pressure grid to velocity grid.

"""
function diver(ux  ::Array
              ,uy  ::Array
              ,mshV::Mesh
              ,Jr,Js)

    uxdx,_    = grad(ux,msh)
    _   ,uydy = grad(uy,msh)

    div = uxdx + uydy

    Bdiv  = mass(div,msh)
    JBdiv = ABu(Js',Jr',Bdiv)

    return JBdiv # pressure grid
end
#--------------------------------------#
export diverᵀ
#--------------------------------------#
"""
 for p ∈ H¹₀(Ω)  (trial pressure)
 for v⃗ ∈ H¹₀(Ω)² (test  velocity)

 -(v⃗,∇p) = -([vx] .* [px])
            ([vy]    [py])

 integrating by parts,

 -(∇v⃗,p) = ([Dx Dy] * [vx])' Bᵥ * J * p
           (          [vy])

         + surface terms (vanish)

         = [vx vy] *  [Dx'] Bᵥ * J * p
                      [Dy']

"""
function diverᵀ(pr  ::Array
               ,mshV::Mesh
               ,Jr,Js)

    Jp  = ABu(Js,Jr,pr) # pr -> vel gird
    BJp = mass(Jp,mshV)

    qx,qy = gradᵀ(BJp,mshV)

    return qx,qy
end
#--------------------------------------#
export stokesOp
#--------------------------------------#
"""
 Schur complement + operator splitting
 operation for stokes solve

 EE = -DD AA^-1 * DD'
"""
function stokesOp(q   ::Array
                 ,mshV::Mesh
                 ,Mvx,Mvy,Jr,Js)

    # DD'
    qx,qy = diverᵀ(q,mshV,Jr,Js)

    # HHinv
    qx = approxHlmzInv(qx,mshV)
    qy = approxHlmzInv(qy,mshV)

    # DD
    Eq = diver(qx,qy,mshV,Jr,Js)

  return -Eq
end
#--------------------------------------#
export approxHlmzInv
#--------------------------------------#
function approxHlmzInv(u::Array,b0::Number
                      ,mshV::Mesh)

    v = gatherScatter(u,mshV)
    v = mask(v,Mvx)

    v = v .* mshV.Bi ./ b0

    v = gatherScatter(v,mshV)
    v = mask(v,Mvx)

    return v
end
#--------------------------------------#


text:
using Revise
using GeoPhyInv
using Statistics
#md using Plots; gr();
using Test

# This tutorial demonstrates how an instance of `SeisForwExpt` can be used iteratively
# to perform forward modeling using `update!` for
# various bundles of medium parameters.

# ### Grid
zgrid=range(-1000.0,stop=1000.0,length=201)
xgrid=range(-1000.0,stop=1000.0,length=201)
mgrid = [zgrid, xgrid];

# ### Medium #1
vpb = [1500., 3500.] # bounds for vp
rhob = [1500., 3500.] # density bounds
medium=Medium(mgrid);
update!(medium, [:vp,:rho], [vpb,rhob]);
fill!(medium);
update!(medium, [:vp,:rho], randn_perc=5); # add some random noise

# ### Medium #2
medium_new=similar(medium)
update!(medium_new, [:vp,:rho], randn_perc=5.) # add some noise
update!(medium_new, [:vp], constant_pert=500.) # perturb velocity

# ### AGeom
ageom=AGeom(medium.mgrid,:surf, SSrcs(3), Recs(30)); # surface seismic

# ### Plotting #1
#md p1=heatmap(medium, :vp)
#md scatter!(ageom, SSrcs())
#md scatter!(ageom, Recs())
#md p2=heatmap(medium_new, :vp)
#md scatter!(ageom, SSrcs())
#md scatter!(ageom, Recs())
#md plot(p1, p2, size=(800,300))

# ### SrcWav
tgrid = range(0.0,stop=1.0,length=1000) # generate a temporal grid
wav = ricker(10.0, tgrid, tpeak=0.25,); # Choose a source wavelet
srcwav = SrcWav(tgrid, ageom, [:p]) # initialize
update!(srcwav, [:p], wav) # distribute to all supersources

# ### SeisForwExpt
pa=SeisForwExpt(FdtdAcoustic(),medium=medium, ageom=ageom, srcwav=srcwav, tgrid=tgrid, verbose=true);

# ### Modeling #1
@time update!(pa);
d1=copy(pa[:data][:p])
#md p1=heatmap(pa[:data], grid=true);

# ### Change medium in `pa` without memory allocation
update!(pa, medium_new)

# ### Modeling #2
@time update!(pa);
d2=copy(pa[:data][:p])
#md p2=heatmap(pa[:data], grid=true);

# Test
@test d1≠d2

# ### Plotting #2
#md plot(p1,p2, size=(500, 300))


text:
# This file was generated, do not modify it. # hide
fig = figure(figsize=(15, 15))

subplot(221)
scatter(data.Solar_gen, data.Radiation_dir)
xlabel("Solar power (kW)", fontsize=14)
ylabel("Solar radiation - directional", fontsize=14)

subplot(222)
scatter(data.Solar_gen, data.Radiation_dif)
xlabel("Solar power (kW)", fontsize=14)
ylabel("Solar radiation - diffuse", fontsize=14)

subplot(223)
scatter(data.Solar_gen, data.Windspeed)
xlabel("Solar power (kW)", fontsize=14)
ylabel("Wind speed (m/s)", fontsize=14)

subplot(224)
scatter(data.Solar_gen, data.Temperature)
xlabel("Solar power (kW)", fontsize=14)
ylabel("Temperature (C)", fontsize=14)

savefig(joinpath(@OUTPUT, "solar_scatter.png"), bbox_inches="tight") # hide


text:
__precompile__(true)

module RefractiveIndex

using YAML, DataDeps

# package code goes here

struct Material
    λ::Array{Float64,1}
    n::Array{Float64,1}
    k::Array{Float64,1}
    nfun::Function
    kfun::Function
    reference::AbstractString
    comment::AbstractString
    specs::Dict{String,Any}
end

include("dispersion_formulae.jl")

export Material, getdata, getshelves, getbooks, getpages

function __init__()
    DATABASEURL = "https://github.com/polyanskiy/refractiveindex.info-database/archive/master.zip"

    RegisterDataDep("RefractiveIndexInfoDB",
        """
        RefractiveIndex.INFO website: © 2008-2018 Mikhail Polyanskiy
        refractiveindex.info database: public domain via CC0 1.0
        NO GUARANTEE OF ACCURACY - Use on your own risk
        """,
        [DATABASEURL],
        post_fetch_method=unpack)

    const global DATAROOT = joinpath(datadep"RefractiveIndexInfoDB",
                              "refractiveindex.info-database-master/database/data/")
end

getshelves() = shelves = readdir(DATAROOT)
getbooks(shelf) = readdir(joinpath(DATAROOT, shelf))

function getpages(shelf, book)
    contents = readdir(joinpath(DATAROOT, shelf, book))
    # Since some pages are located in a subdirectory we need to extract those
    for (i, content) in enumerate(contents)
        if isdir(joinpath(DATAROOT, shelf, book, content))
            dircontents = readdir(joinpath(DATAROOT, shelf, book, content))
            splice!(contents, i, joinpath.(content, dircontents))
        end
    end
    contents
end

parsearray(datastr) = "[" * datastr * "]" |> parse |> eval

function getdata(shelf, book, page)
    isempty(splitext(page)[2]) ? page *= ".yml" : nothing
    path = joinpath(DATAROOT, shelf, book, page)
    yaml = YAML.load(open(path))
    haskey(yaml, "REFERENCES") ? reference = yaml["REFERENCES"] : reference = ""
    haskey(yaml, "COMMENTS") ? comment = yaml["COMMENTS"] : comment = ""
    haskey(yaml, "SPECS") ? specs = yaml["SPECS"] : specs = Dict{String,Any}()

    # Default field contents for material
    λ = Float64[]
    n = Float64[]
    k = Float64[]
    nfun() = nothing
    kfun() = nothing

    for i in 1:length(yaml["DATA"])
        # Possible datatypes:
        #   -tablulated nk
        #   -tablulated k
        #   -formula 1-9
        datatype = yaml["DATA"][i]["type"]

        if any(datatype .== ["tabulated nk", "tabulated n", "tabulated k"])
            datastr = yaml["DATA"][i]["data"]
            data = parsearray(datastr)
        elseif any(datatype .== ["formula 1",
                                 "formula 2",
                                 "formula 3",
                                 "formula 4",
                                 "formula 5",
                                 "formula 6",
                                 "formula 7",
                                 "formula 8",
                                 "formula 9"])
            rangestr = yaml["DATA"][i]["range"]
            coeffstr = yaml["DATA"][i]["coefficients"]
            range = parsearray(rangestr)
            coeffs = parsearray(coeffstr)
        else
            error("Datatype $datatype unknown")
        end

        if datatype == "tabulated nk"
            λ = data[:,1]
            n = data[:,2]
            k = data[:,3]
            nfun = splinefun(λ, n, :n)
            kfun = splinefun(λ, k, :k)
            range = [maximum(λ), minimum(λ)]
        elseif datatype == "tabulated n"
            λ = data[:,1]
            n = data[:,2]
            nfun = splinefun(λ, n, :n)
            range = [minimum(λ), maximum(λ)]
        elseif datatype == "tabulated k"
            λ = data[:,1]
            k = data[:,2]
            kfun = splinefun(λ, k, :k)
            range = [minimum(λ), maximum(λ)]
        elseif datatype == "formula 1"
            nfun = eval(f1(range, coeffs))
        elseif datatype == "formula 2"
            nfun = eval(f2(range, coeffs))
        elseif datatype == "formula 3"
            nfun = eval(f3(range, coeffs))
        elseif datatype == "formula 4"
            nfun = eval(f4(range, coeffs))
        elseif datatype == "formula 5"
            nfun = eval(f5(range, coeffs))
        elseif datatype == "formula 6"
            nfun = eval(f6(range, coeffs))
        elseif datatype == "formula 7"
            nfun = eval(f7(range, coeffs))
        elseif datatype == "formula 8"
            nfun = eval(f8(range, coeffs))
        elseif datatype == "formula 9"
            nfun = eval(f9(range, coeffs))
        else
            error("Datatype \"$(datatype)\" not defined.")
        end
    end
    return Material(λ, n, k, nfun, kfun, reference, comment, specs)
end

end # module


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

JLLWrappers.@generate_main_file_header("assimp")
JLLWrappers.@generate_main_file("assimp", UUID("54ae6823-98c6-5a7c-8365-5a43b909f91f"))
end  # module assimp_jll


text:
## runtests.jl : run all test files

# using PyPlot   # *****
using GaussianRandomFields
using Test
using Suppressor
using DelimitedFiles
using Pkg

# test indexsets
include("test_covariance_functions.jl")
include("test_cholesky.jl")
include("test_spectral.jl")
include("test_karhunen_loeve.jl")
include("test_circulant_embedding.jl")
include("test_gaussian_random_fields.jl")
include("test_separable.jl")
include("test_fem.jl")
#include("test_plot.jl")    # *****
include("test_show.jl")


text:
using NumericalRationalFunctions
using LinearAlgebra
using Test

@testset "AAA" begin
    Z = range(-1, stop=1, length=20)
    F = exp.(Z)
    r, pol, res, zer, z, f, w, errvec = aaa(F, Z)
    @test length(z) == 7
    @test length(pol) == length(res) == length(zer) == 6
    @test r(Z) ≈ F atol=1e-14
    x = 0.12345
    @test r(x) ≈ exp(x) atol=1e-14
    @test norm(r(zer), Inf) < 1e-9
    @test errvec[end] < 1e-13
end

@testset "Laplace" begin
    u, err, f, Z, Zplot, A, inpoly = laplace(:iso, tol=1e-10, boundary=z->log(abs(z)))
    @test err < 1e-10
    u, err, f, Z, Zplot, A, inpoly = laplace(:circleL)
    @test err < 1e-6
    u, err, f, Z, Zplot, A, inpoly = laplace(complex.([[1, 0.5], 1+1im, -1+1im, -1]), boundary=[1, 0, 0, 0])
    @test err < 1e-6

    #=
    sx = sy = range(-1.5, stop=1.5, length=200)
    p = contour(sx, sy, (x, y)->begin
                    z = x + im*y
                    inpoly(z) == 1 ? u(z) : NaN
                end, levels=30, c=:viridis, aspect_ratio=1)
    nw = length(Zplot)
    for k in 1:nw-1
        kn = k+1
        plot!(p, [Zplot[k], Zplot[kn]], color=:black, lab=false, linewidth=2)
    end
    p
    =#
end


text:
function CellRef(n::AbstractString)
    @assert is_valid_cellname(n) "$n is not a valid CellRef."
    column_name, row_number = split_cellname(n)
    return CellRef(n, row_number, decode_column_number(column_name))
end

@inline CellRef(row::Int, col::Int) = CellRef(encode_column_number(col) * string(row))
@inline CellPosition(ref::CellRef) = CellPosition(row_number(ref), column_number(ref))
@inline row_number(p::CellPosition) = p.row
@inline column_number(p::CellPosition) = p.column
@inline CellRef(p::CellPosition) = CellRef(row_number(p), column_number(p))

"""
    decode_column_number(column_name::AbstractString) :: Int

Converts column name to a column number.

```julia
julia> XLSX.decode_column_number("D")
4
```

See also: `encode_column_number`.
"""
function decode_column_number(column_name::AbstractString) :: Int
    local result::Int = 0

    num_characters = length(column_name)

    iteration = 1
    for i in num_characters:-1:1
        column_char_as_int = Int(column_name[i])
        result += (26^(iteration-1)) * (column_char_as_int - 64) # From A to Z we have 26 values. 'A' Char is ASCII 65.
        iteration += 1
    end

    return result
end

"""
    encode_column_number(column_number::Int) :: String

Converts column number to a column name.

# Example

```julia
julia> XLSX.encode_column_number(4)
"D"
```

See also: `decode_column_number`.
"""
function encode_column_number(column_number::Int) :: String
    @assert column_number > 0 && column_number <= 16384 "Column number should be in the range from 1 to 16384."

    third_letter_sequence = div(column_number - 26 - 1, 26*26)
    column_number = column_number - third_letter_sequence*(26*26)

    second_letter_sequence = div(column_number - 1, 26) # 26^1
    column_number = column_number - second_letter_sequence*(26)

    first_letter_sequence = column_number # 26^0

    if third_letter_sequence > 0
        # result will have 3 letters
        return String([ Char(third_letter_sequence+64), Char(second_letter_sequence+64), Char(first_letter_sequence+64) ])

    elseif second_letter_sequence > 0
        # result will have 2 letters
        return String([ Char(second_letter_sequence+64), Char(first_letter_sequence+64) ])

    else
        # result will have 1 letter
        return String([ Char(first_letter_sequence+64) ])
    end
end

Base.string(c::CellRef) = c.name
Base.show(io::IO, c::CellRef) = print(io, string(c))
Base.:(==)(c1::CellRef, c2::CellRef) = c1.name == c2.name
Base.hash(c::CellRef) = hash(c.name)

const RGX_COLUMN_NAME = r"^[A-Z]?[A-Z]?[A-Z]$"
const RGX_CELLNAME = r"^[A-Z]+[0-9]+$"
const RGX_CELLRANGE = r"^[A-Z]+[0-9]+:[A-Z]+[0-9]+$"

function is_valid_column_name(n::AbstractString) :: Bool
    if !occursin(RGX_COLUMN_NAME, n)
        return false
    end

    column_number = decode_column_number(n)
    if column_number < 1 || column_number > 16384
        return false
    end

    return true
end

const RGX_CELLNAME_LEFT = r"^[A-Z]+"
const RGX_CELLNAME_RIGHT = r"[0-9]+$"

"""
    split_cellname(n::AbstractString) -> column_name, row_number

Splits a string representing a cell name to its column name and row number.

# Example

```julia
julia> XLSX.split_cellname("AB:12")
("AB:", 12)
```
"""
@inline function split_cellname(n::AbstractString)
    for (i, c) in enumerate(n)
        if isdigit(c)
            column_name = SubString(n, 1, i-1)
            row = parse(Int, SubString(n, i, length(n)))

            return column_name, row
        end
    end

    error("Couldn't split (column_name, row) for cellname $n.")
end

"""
    is_valid_cellname(n::AbstractString) :: Bool

Checks wether `n` is a valid name for a cell.

Cell names are bounded by `A1 : XFD1048576`.
"""
function is_valid_cellname(n::AbstractString) :: Bool

    if !occursin(RGX_CELLNAME, n)
        return false
    end

    column_name, row = split_cellname(n)

    if row < 1 || row > 1048576
        return false
    end

    if !is_valid_column_name(column_name)
        return false
    end

    return true
end

const RGX_CELLRANGE_START = r"^[A-Z]+[0-9]+"
const RGX_CELLRANGE_STOP = r"[A-Z]+[0-9]+$"

"""
    split_cellrange(n::AbstractString) -> start_name, stop_name

Splits a string representing a cell range into its cell names.

# Example

```julia
julia> XLSX.split_cellrange("AB12:CD24")
("AB12", "CD24")
```
"""
@inline function split_cellrange(n::AbstractString)
    s = split(n, ":")
    @assert length(s) == 2 "$n is not a valid cell range."
    return s[1], s[2]
end

function is_valid_cellrange(n::AbstractString) :: Bool

    if !occursin(RGX_CELLRANGE, n)
        return false
    end

    start_name, stop_name = split_cellrange(n)

    if !is_valid_cellname(start_name)
        return false
    end

    if !is_valid_cellname(stop_name)
        return false
    end

    return true
end

macro ref_str(ref)
    CellRef(ref)
end

function CellRange(r::AbstractString)
    @assert occursin(RGX_CELLRANGE, r) "Invalid cell range: $r."
    start_name, stop_name = split_cellrange(r)
    return CellRange(CellRef(start_name), CellRef(stop_name))
end

CellRange(start_row::Integer, start_column::Integer, stop_row::Integer, stop_column::Integer) = CellRange(CellRef(start_row, start_column), CellRef(stop_row, stop_column))

Base.string(cr::CellRange) = "$(string(cr.start)):$(string(cr.stop))"
Base.show(io::IO, cr::CellRange) = print(io, string(cr))
Base.:(==)(cr1::CellRange, cr2::CellRange) = cr1.start == cr2.start && cr2.stop == cr2.stop
Base.hash(cr::CellRange) = hash(cr.start) + hash(cr.stop)

macro range_str(cellrange)
    CellRange(cellrange)
end

"""
    Base.in(ref::CellRef, rng::CellRange) :: Bool

Checks wether `ref` is a cell reference inside a range given by `rng`.
"""
function Base.in(ref::CellRef, rng::CellRange) :: Bool
    top = row_number(rng.start)
    bottom = row_number(rng.stop)
    r = row_number(ref)

    if top <= r && r <= bottom
        left = column_number(rng.start)
        right = column_number(rng.stop)
        c = column_number(ref)

        if left <= c && c <= right
            return true
        end
    end

    return false
end

"""
    Base.issubset(subrng::CellRange, rng::CellRange)

Checks wether `subrng` is a cell range contained in `rng`.
"""
Base.issubset(subrng::CellRange, rng::CellRange) :: Bool = in(subrng.start, rng) && in(subrng.stop, rng)

function Base.size(rng::CellRange)
    top = row_number(rng.start)
    bottom = row_number(rng.stop)
    left = column_number(rng.start)
    right = column_number(rng.stop)

    return ( bottom - top + 1, right - left + 1 )
end

"""
    row_number(c::CellRef) :: Int

Returns the row number of a given cell reference.
"""
row_number(c::CellRef) :: Int = c.row_number

"""
    column_number(c::CellRef) :: Int

Returns the column number of a given cell reference.
"""
column_number(c::CellRef) :: Int = c.column_number

column_name(c::CellRef) :: String = encode_column_number(column_number(c))

"""
Returns (row, column) representing a `ref` position relative to `rng`.

For example, for a range "B2:D4", we have:

* "C3" relative position is (2, 2)

* "B2" relative position is (1, 1)

* "C4" relative position is (3, 2)

* "D4" relative position is (3, 3)

"""
function relative_cell_position(ref::CellRef, rng::CellRange)
    @assert ref ∈ rng "$ref is outside range $rng."

    top = row_number(rng.start)
    left = column_number(rng.start)

    r, c = row_number(ref), column_number(ref)

    return ( r - top + 1 , c - left + 1 )
end

#
# ColumnRange
#

Base.string(cr::ColumnRange) = "$(encode_column_number(cr.start)):$(encode_column_number(cr.stop))"
Base.show(io::IO, cr::ColumnRange) = print(io, string(cr))
Base.:(==)(cr1::ColumnRange, cr2::ColumnRange) = cr1.start == cr2.start && cr2.stop == cr2.stop
Base.hash(cr::ColumnRange) = hash(cr.start) + hash(cr.stop)
Base.in(column_number::Integer, rng::ColumnRange) = rng.start <= column_number && column_number <= rng.stop

function relative_column_position(column_number::Integer, rng::ColumnRange)
    @assert column_number ∈ rng "Column $column_number is outside range $rng."
    return column_number - rng.start + 1
end

@inline relative_column_position(ref::CellRef, rng::ColumnRange) = relative_column_position(column_number(ref), rng)

const RGX_COLUMN_RANGE = r"^[A-Z]?[A-Z]?[A-Z]:[A-Z]?[A-Z]?[A-Z]$"
const RGX_COLUMN_RANGE_START = r"^[A-Z]+"
const RGX_COLUMN_RANGE_STOP = r"[A-Z]+$"

"""
Returns tuple (column_name_start, column_name_stop).
"""
@inline function split_column_range(n::AbstractString)
    s = split(n, ":")
    return s[1], s[2]
end

function is_valid_column_range(r::AbstractString) :: Bool
    if !occursin(RGX_COLUMN_RANGE, r)
        return false
    end

    start_name, stop_name = split_column_range(r)

    if !is_valid_column_name(start_name) || !is_valid_column_name(stop_name)
        return false
    end

    return true
end

function ColumnRange(r::AbstractString)
    @assert is_valid_column_range(r) "Invalid column range: $r."
    start_name, stop_name = split_column_range(r)
    return ColumnRange(decode_column_number(start_name), decode_column_number(stop_name))
end

convert(::Type{ColumnRange}, str::AbstractString) = ColumnRange(str)
convert(::Type{ColumnRange}, column_range::ColumnRange) = column_range

column_bounds(r::ColumnRange) = (r.start, r.stop)
Base.length(r::ColumnRange) = r.stop - r.start + 1

<<<<<<< HEAD
# ColumnRange iterator
start(itr::ColumnRange) = itr.start
done(itr::ColumnRange, column_index::Int) = column_index > itr.stop
next(itr::ColumnRange, column_index::Int) = (encode_column_number(column_index), column_index + 1)

# CellRange iterator
start(rng::CellRange) = CellRefIteratorState(row_number(rng.start), column_number(rng.start))
done(rng::CellRange, state::CellRefIteratorState) = state.row > row_number(rng.stop)
=======
# ColumnRange iterator: element is a String with the column name, the state is the column number.
function Base.iterate(itr::ColumnRange, state::Int=itr.start)
    if state > itr.stop
        return nothing
    end
>>>>>>> upstream/master

    return encode_column_number(state), state + 1
end

<<<<<<< HEAD
# (i, state) = next(I, state)
function next(rng::CellRange, state::CellRefIteratorState)
    local next_state::CellRefIteratorState
    if state.col == column_number(rng.stop)
=======
# CellRange iterator: element is a CellRef, the state is a CellPosition.
function Base.iterate(rng::CellRange, state::CellPosition=CellPosition(rng.start))

    if row_number(state) > row_number(rng.stop)
        return nothing
    elseif column_number(state) == column_number(rng.stop)
>>>>>>> upstream/master
        # reached last column. Go to the next row.
        next_state = CellPosition(row_number(state) + 1, column_number(rng.start))
    else
        # go to the next column
        next_state = CellPosition(row_number(state), column_number(state) + 1)
    end

    return CellRef(state), next_state
end

function Base.length(rng::CellRange)
    (r, c) = size(rng)
    return r * c
end

#
# SheetCellRef, SheetCellRange, SheetColumnRange
#

Base.string(cr::SheetCellRef) = string(cr.sheet, "!", cr.cellref)
Base.show(io::IO, cr::SheetCellRef) = print(io, string(cr))
Base.:(==)(cr1::SheetCellRef, cr2::SheetCellRef) = cr1.sheet == cr2.sheet && cr2.cellref == cr2.cellref
Base.hash(cr::SheetCellRef) = hash(cr.sheet) + hash(cr.cellref)

Base.string(cr::SheetCellRange) = string(cr.sheet, "!", cr.rng)
Base.show(io::IO, cr::SheetCellRange) = print(io, string(cr))
Base.:(==)(cr1::SheetCellRange, cr2::SheetCellRange) = cr1.sheet == cr2.sheet && cr2.rng == cr2.rng
Base.hash(cr::SheetCellRange) = hash(cr.sheet) + hash(cr.rng)

Base.string(cr::SheetColumnRange) = string(cr.sheet, "!", cr.colrng)
Base.show(io::IO, cr::SheetColumnRange) = print(io, string(cr))
Base.:(==)(cr1::SheetColumnRange, cr2::SheetColumnRange) = cr1.sheet == cr2.sheet && cr2.colrng == cr2.colrng
Base.hash(cr::SheetColumnRange) = hash(cr.sheet) + hash(cr.colrng)

const RGX_SHEET_CELLNAME = r"^.+![A-Z]+[0-9]+$"
const RGX_SHEET_CELLRANGE = r"^.+![A-Z]+[0-9]+:[A-Z]+[0-9]+$"
const RGX_SHEET_COLUMN_RANGE = r"^.+![A-Z]?[A-Z]?[A-Z]:[A-Z]?[A-Z]?[A-Z]$"

const RGX_SHEET_CELLNAME_RIGHT = r"[A-Z]+[0-9]+$"
const RGX_SHEET_CELLRANGE_RIGHT = r"[A-Z]+[0-9]+:[A-Z]+[0-9]+$"
const RGX_SHEET_COLUMN_RANGE_RIGHT = r"[A-Z]?[A-Z]?[A-Z]:[A-Z]?[A-Z]?[A-Z]$"

function is_valid_sheet_cellname(n::AbstractString) :: Bool
    if !occursin(RGX_SHEET_CELLNAME, n)
        return false
    end

    cellname = match(RGX_SHEET_CELLNAME_RIGHT, n).match
    if !is_valid_cellname(cellname)
        return false
    end

    return true
end

function is_valid_sheet_cellrange(n::AbstractString) :: Bool
    if !occursin(RGX_SHEET_CELLRANGE, n)
        return false
    end

    cellrange = match(RGX_SHEET_CELLRANGE_RIGHT, n).match
    if !is_valid_cellrange(cellrange)
        return false
    end

    return true
end

function is_valid_sheet_column_range(n::AbstractString) :: Bool
    if !occursin(RGX_SHEET_COLUMN_RANGE, n)
        return false
    end

    column_range = match(RGX_SHEET_COLUMN_RANGE_RIGHT, n).match
    if !is_valid_column_range(column_range)
        return false
    end

    return true
end

const RGX_CELLNAME_RIGHT_FIXED = r"\$[A-Z]+\$[0-9]+$"
const RGX_SHEET_CELNAME_RIGHT_FIXED = r"\$[A-Z]+\$[0-9]+:\$[A-Z]+\$[0-9]+$"

function SheetCellRef(n::AbstractString)
    if is_valid_fixed_sheet_cellname(n)
        fixed_cellname = match(RGX_CELLNAME_RIGHT_FIXED, n).match
        cellname = replace(fixed_cellname, "\$" => "")
        sheetname = SubString(n, 1, length(n) - length(fixed_cellname) - 1)
        return SheetCellRef(sheetname, CellRef(cellname))
    else
        @assert is_valid_sheet_cellname(n) "$n is not a valid SheetCellRef."
        cellname = match(RGX_SHEET_CELLNAME_RIGHT, n).match
        sheetname = SubString(n, 1, length(n) - length(cellname) - 1)
        return SheetCellRef(sheetname, CellRef(cellname))
    end
end

function SheetCellRange(n::AbstractString)
    if is_valid_fixed_sheet_cellrange(n)
        fixed_cellrange = match(RGX_SHEET_CELNAME_RIGHT_FIXED, n).match
        cellrange = replace(fixed_cellrange, "\$" => "")
        sheetname = SubString(n, 1, length(n) - length(fixed_cellrange) - 1)
        return SheetCellRange(sheetname, CellRange(cellrange))
    else
        @assert is_valid_sheet_cellrange(n) "$n is not a valid SheetCellRange."
        cellrange = match(RGX_SHEET_CELLRANGE_RIGHT, n).match
        sheetname = SubString(n, 1, length(n) - length(cellrange) - 1)
        return SheetCellRange(sheetname, CellRange(cellrange))
    end
end

function SheetColumnRange(n::AbstractString)
    @assert is_valid_sheet_column_range(n) "$n is not a valid SheetColumnRange."
    column_range = match(RGX_SHEET_COLUMN_RANGE_RIGHT, n).match
    sheetname = SubString(n, 1, length(n) - length(column_range) - 1)
    return SheetColumnRange(sheetname, ColumnRange(column_range))
end

# Named ranges
const RGX_FIXED_SHEET_CELLNAME = r"^.+!\$[A-Z]+\$[0-9]+$"
const RGX_FIXED_SHEET_CELLRANGE = r"^.+!\$[A-Z]+\$[0-9]+:\$[A-Z]+\$[0-9]+$"

is_valid_fixed_sheet_cellname(s::AbstractString) = occursin(RGX_FIXED_SHEET_CELLNAME, s)
is_valid_fixed_sheet_cellrange(s::AbstractString) = occursin(RGX_FIXED_SHEET_CELLRANGE, s)


text:
# Pkg.clone("https://github.com/wookay/GlyphSheets.jl.git")
using GlyphSheets
import GlyphSheets.DSL: Move, Line, Cubic, Conic

library = create_library()
typeface = create_typeface(library, "Butler_Regular.otf")

glyph = create_glyph(typeface, 'B', 5pt)
outlines = outline_decompose(glyph)

close_typeface(typeface)
close_library(library)

using Gtk, Gtk.ShortNames, Graphics, Cairo
using Colors

win = Window("Gtk", 300, 300)
hbox = Box(:h)
setproperty!(hbox, :homogeneous, true)
push!(win, hbox)
canv = Canvas()
push!(hbox, canv)

function operate(cr, op)
    if op isa Move
        move_to(cr, op.to...)
    elseif op isa Line
        line_to(cr, op.to...)
    elseif op isa Conic
        curve_to(cr, op.control..., op.to..., op.to...)
    elseif op isa Cubic
        curve_to(cr, op.control1..., op.control2..., op.to...)
    end
end

function fill_paths(cr, outlines)
    save(cr)
    operate.(cr, outlines)
    set_source_rgba(cr, 0.9, 1, 0.9, 1)
    fill(cr)
    restore(cr)
end

function stroke_paths(cr, outlines)
    save(cr)
    set_line_width(cr, 3)
    colors = distinguishable_colors(length(outlines))
    last_pos = (0, 0)
    for (color, op) in zip(colors, outlines)
        move_to(cr, last_pos...)
        operate(cr, op)
        set_source_rgba(cr, red(color), green(color), blue(color), alpha(color))
        stroke(cr)
        last_pos = op.to
    end
    restore(cr)
end

@guarded draw(canv) do widget
    cr = getgc(canv)

    translate(cr, 0, glyph.slot.metrics.height+20)
    scale(cr, 1, -1)

    fill_paths(cr, outlines)
    stroke_paths(cr, outlines)
end # guarded

showall(win)

cond = Condition()
signal_connect(win, :destroy) do widget
    notify(cond)
end
wait(cond)


text:
function analyzeleadlag(;measure::DataFrame, measureinfo::NamedTuple,
  llname = PARAM[:llname])
  @unpack outputpath, resultid, returntype, datefrequency = measureinfo

  measure, llinfo = prepleadlag(; measure, measureinfo)
  llresults = analyzeleadlag(measure, measureinfo, llinfo)

  sort!(llresults, :betaxr)

  llresults |> CSV.write("$outputpath\\$(resultid)_$(llname)_results.csv")
  measure |> CSV.write("$outputpath\\$(resultid)_$(llname)_series.csv")

  return nothing
end

#this lags and leads the focal measure and return columns
function prepleadlag(;
  measure::DataFrame,
  measureinfo::NamedTuple,
  measureleadlags::Vector{Int} = PARAM[:llmeasureleadlags][measureinfo.datefrequency],
  returnwindows::Vector{Int} = PARAM[:llreturnwindows][measureinfo.datefrequency],
  additionalreturncolsidx = PARAM[:lladditionalreturncolsidx],
  additionalreturncols = additionalreturncolsidx |> keys |> collect,
  measurecolpatterns = PARAM[:llmeasurecolpatterns],
  detrendll = PARAM[:lldetrendll],
  additionallagcols = PARAM[:lladditionallagcols],
  controls= PARAM[:llcontrols])

  #NOTE: the below are a couple fast sanity checks
  #testshift(10)
  #testforwardreturnwindow(24)

  initialreturncols = [findcols(measure, :return); additionalreturncols]
  initialmeasurecols = [(p->findcols(measure, p)).(measurecolpatterns)...;]

  @assert !isempty(initialmeasurecols)
  @assert !isempty(initialreturncols)

  #the below will include shifted cols and variations on the return window
  measurecols = Symbol[]
  returncols = Symbol[]

  additionallagcols = [(s->findcols(measure, s)).(additionallagcols)...;]
  #fore each lead and lag and each measure col, apply the appropriate lagging
  for l ∈ measureleadlags
    prefix = Symbol(l>0 ? :N : :L, abs(l))
    for mcol ∈ initialmeasurecols
      shiftedmcol = Symbol(prefix, mcol)
      measure[!, shiftedmcol] = shift(measure[!, mcol], l)
      push!(measurecols, shiftedmcol) #record the newly created column
    end

    for lcol ∈ additionallagcols
      shiftedlcol = Symbol(prefix, lcol)
      measure[!, shiftedlcol] = shift(measure[!, lcol], l)
    end
  end

  #create the return windows
  for windowsize ∈ returnwindows
    prefix = Symbol(:FW, windowsize)
    for rcol ∈ initialreturncols
      returnwindowcol = Symbol(
        prefix, rcol ∈ additionalreturncols ? "_" : "", rcol)
      measure[!, returnwindowcol] = forwardreturnwindow(
        measure[!, rcol]; windowsize, logged=false)
      push!(returncols, returnwindowcol)
    end
  end

  #this will hold the lead lag pairs
  estimations = NamedTuple{(:rcol, :mcol, :ccols, :lag)}[]

  function extractlag(colname::String)
    prefix =  match(r"^L[0-9]+[A-Za-z]*_",colname)
    lagstring = match(r"[0-9]+", prefix)
    return parse(Int, lagstring)
  end
  extractlag(colname::Symbol) = extractlag(colname |> string)

  for rcol ∈ returncols, mcol ∈ measurecols
    local makepair
    strippedrcol = stripcolprefix(rcol)|>string
    strippedmcol = stripcolprefix(mcol)|>string

    #if the stripped versions match assume a pairing is desired
    if occursin(strippedrcol, strippedmcol)
      makepair=true
    #see if we explicitly want a pairing
    elseif (haskey(additionalreturncolsidx, strippedrcol |> Symbol) &&
        any((m->occursin(m|>string, strippedmcol)).(
          additionalreturncolsidx[strippedrcol |> Symbol])))
        makepair=true
    else
      makepair = false
    end
  #println("mcol: $mcol, rcol: $rcol,",
    #  " strippedmcol: $strippedmcol, strippedrcol: $strippedrcol",
    #  " makepair: $makepair")

    if makepair
      for control ∈ controls
        @unpack ccols, controlmatch = control
          if (controlmatch≡nothing) || (match(Regex(controlmatch), mcol |> string) !== nothing)
            if mcol ∈ ccols
              (length(ccols)==1) && continue #this means none of the controls are valid
              setdiff!(ccols, [mcol])
            end
            push!(estimations, (; rcol, mcol, ccols))
          end
        end
    end
  end

  #linearly detrends the results
  if detrendll
    @assert issorted(measure, [:date])

    originalestimations = estimations |> deepcopy
    for orig ∈ originalestimations
      @unpack rcol,mcol, ccols = orig

      #for now do not detrend the controls
      dtrcol, dtmcol = Symbol(:dtX,rcol), Symbol(:dtX,mcol)

      #do the actual detrending if necessary, being careful of missing values
      if dtrcol ∉ propertynames(measure)
        measure[!, dtrcol] = missings(Float64, nrow(measure))
        cleanmeasure = view(measure, completecases(measure, [rcol;]), :)
        cleanmeasure[:, dtrcol] .= detrend(cleanmeasure[:, rcol])
      end

      if dtmcol ∉ propertynames(measure)
        measure[!, dtmcol] = missings(Float64, nrow(measure))
        cleanmeasure = view(measure, completecases(measure, [mcol;]), :)
        cleanmeasure[:, dtmcol] .= detrend(cleanmeasure[:, mcol])
      end

      push!(estimations, (;rcol = dtrcol, mcol = dtmcol, ccols))
    end
  end

  llinfo = (; returncols, measurecols, measureleadlags, returnwindows, estimations)

  measure |> CSV.write(
    "$(PARAM[:testpath])\\$(PARAM[:analysismeasurefilename])_prepll.csv")
  return measure, llinfo
end

#function for simple linear detrending
function detrend(col; idx = 1:length(col) |> collect .|> Float64, idx1 = [idx ones(length(col))])
  @assert length(col) == length(idx) == size(idx1,1)

  m = cov(col, idx)/var(idx)
  b = mean(col) - m*mean(idx)

  chk = FMLM(idx1, col |> Vector{Float64})
  @assert m+10.0 ≈ chk.β[1]+10.0 "m=$m, chk.β[1]=$(chk.β[1])"
  @assert b+10.0 ≈ chk.β[2]+10.0

  return col .- (idx .* m .+ b)
end

function shift(v::AbstractVector{MFloat64}, i::Int)
  vout = circshift(v,-i)
  vout[ifelse(i>0, (end-i+1):end, 1:abs(i))] .= missing
  return vout
end

shift(v::AbstractVector, i::Int) = (i > 0 ?
  [v[(i+1):end]; missings(Float64,i);] :
  [missings(Float64, abs(i)); v[1:(end+i)];])

#averages returns over a forward window
function forwardreturnwindow(vin::AbstractVector{T}; windowsize, logged::Bool) where T
  v = deepcopy(vin)
  (windowsize==1) && return v #this is an identity scenario

  if !logged #compute the log returns
    @assert (vᵢ->vᵢ≡missing || vᵢ > -1.).(v) |> all
    v .= log.(v .+ 1.0)
  end

  #compute the window of lgoged returns
  cumlogret = v |> deepcopy |> Vector{Union{T,Missing}}

  for t ∈ 2:windowsize
    shiftedv = shift(v, t-1)
    @assert (shiftedv[(end  - t + 2) : end] .≡ missing) |> all
    cumlogret .+= shiftedv
  end

  logged && return cumlogret #unlog the returns
  return exp.(cumlogret) .- 1
end

function testshift(N=10)
  #throw("something is off, maybe with circshift")
  v = rand(N)
  @assert (shift(v,3) .≡ [v[4:end]; missings(3)]) |> all
  @assert (shift(v,-3) .≡ [missings(3); v[1:(end-3)]]) |> all
  @assert (shift(v,0) .≡ v) |> all

  vm = v |> Vector{MFloat64}
  @assert (shift(vm,3) .≡ [v[4:end]; missings(3)]) |> all "vm: $(shift(vm,3)); ver: $([v[4:end]; missings(3)])"
  @assert (shift(vm,-3) .≡ [missings(3); v[1:(end-3)]]) |> all
  @assert (shift(vm,0) .≡ v) |> all
end

function testforwardreturnwindow(N=24)
  v = rand(N) .- 0.5
  N1v = [v[2:end]; missing]
  N2v = [v[3:end]; missing; missing]

  approxormissing(x,y) = ((x≡y) | coalesce((x≈y),false))

  @assert (v .≈
    forwardreturnwindow(v,windowsize=1, logged=false)) |> all
  @assert approxormissing.((1 .+ v) .* (1 .+ N1v) .- 1,
    forwardreturnwindow(v,windowsize=2, logged=false)) |> all
  @assert approxormissing.((1 .+ v) .* (1 .+ N1v) .* (1 .+ N2v) .- 1,
    forwardreturnwindow(v,windowsize=3, logged=false)) |> all

  @assert (v .≈
    forwardreturnwindow(v,windowsize=1, logged=true)) |> all
  @assert approxormissing.((v .+ N1v),
    forwardreturnwindow(v,windowsize=2, logged=true)) |> all
  @assert approxormissing.((v .+ N1v .+ N2v),
    forwardreturnwindow(v,windowsize=3, logged=true)) |> all
end

#currently only analyzes the lag
function analyzeleadlag(measure, measureinfo, llinfo)

  #=function bsstatfunc(sim)
    r,m = sim[1], sim[2]

    ρ = cor(r,m)
    return Dict(
      :betaxr => ρ * std(r)/std(m),
      #:betamr => ρ * std(m)/std(r)
      )
  end

  nullar1statfunc(r,m) = bsstatfunc((r,m))=#
  function nullar1statfunc(Xy)
    X = @view Xy[:, 1:(end-1)]
    y = @view Xy[:, end]
    #println("X[1:3,:]", X[1:3,:])
    #println("r[1:3]", r[1:3])
    Dict(:betaxr => cholesky!(Symmetric(X'*X))\(X'*y),)
  end

  @unpack (neweylags, nwΣ,
    bsaggfunc, bsmeasurecolsregex, bsbenchmarkcolsregex, bssamples, bswidth,
      nullar1aggfunc, nullar1simulations) = measureinfo
  @unpack returncols, measurecols, measureleadlags, returnwindows, estimations = llinfo

  #create a container for the resutls
  results =   DataFrame(
    label=Symbol[], #label in case we concatenate this df with others
    controls=String[],
    measure=Symbol[], #measure field that we are benchmarking
    returnwindow=Symbol[], #benchmark field

    betaxr=MFloat64[],
      betaxrmw=MFloat64[],
      betaxrhac=MFloat64[],
      betaxrmbbp=MFloat64[],
      betaxrmbbz=MFloat64[],
      betaxrar1h0z=MFloat64[],
      betaxrar1h0p=MFloat64[],
      ibetaxrar1h0z=MFloat64[],
      ibetaxrar1h0p=MFloat64[],
      N=MInt[])

  #iterate through the cartesian product of column pairins
  for estimation ∈ estimations

    @unpack rcol, mcol, ccols = estimation
    teststr ="testing $rcol ~ $mcol + $ccols"

    occursin(r"^N", string(mcol)) && continue #an N would imply a lead column, not a lag column

    #make sure we have no missing entries
    jointnomissing = view(measure, completecases(measure, [rcol;mcol; ccols;]),
      [:date; rcol;mcol;ccols;])

    N=nrow(jointnomissing)
    r::Vector{Float64} = jointnomissing[:,rcol]
    #println([mcol, ccols])
    Xm::Matrix{Float64} = hcat(jointnomissing[:,[mcol; ccols]], ones(N)) |> Matrix{Float64}

    #pear = cor(r,m) #pearson correlation

    lmrm = FMLM(Xm, r)
    betaxr=lmrm.β[1]
    @assert ((lmrm.β[1] ≈ (qr(Xm)\r)[1]) || (lmrm.β[1] ≡ (qr(Xm)\r)[1]))
    betaxrmw = betaxr/modifiedwhiteΣ!(lmrm)[1,1]^0.5
    betaxrhac = betaxr/nwΣ(lmrm)[1,1]^0.5

    #to save time only use the bootstrap and simulation some of the time
    if occursin(bsmeasurecolsregex, string(mcol))
      bs = regmbb(Xm, r; bssamples, bswidth)

      #check dimensions of results line up
      @assert length(bs.beta.se) == length(bs.beta.p) == length(ccols) + 2 "
        bs.beta: $(bs.beta) ccols: $ccols"

      betaxrmbbz =  betaxr/bs.beta.se[1]
      betaxrmbbp =  bs.beta.p[1]

      println(teststr)
      ar1h0teststats = Dict(:betaxr => betaxr#=, :betamr => betamr=#)
      ar1h0 = comparear1undernull(Xm,r, nullar1statfunc,
        (M,s)->nullar1aggfunc(M,s,ar1h0teststats),
        Nsimulations=nullar1simulations)
      #println("$teststr: ar1 complete")

      @assert length(ar1h0.betaxr.se) == length(ar1h0.betaxr.p) == length(ccols) + 2

      betaxrar1h0z =  betaxr/ar1h0.betaxr.se[1]
      betaxrar1h0p =  ar1h0.betaxr.p[1]

      iar1h0 = comparearima11undernull(Xm, r, nullar1statfunc,
        (M,s)->nullar1aggfunc(M,s,ar1h0teststats),
        Nsimulations=nullar1simulations)
      #println("$teststr: iar1 complete")
      @assert length(iar1h0.betaxr.se) == length(iar1h0.betaxr.p) == length(ccols) + 2

      ibetaxrar1h0z =  betaxr/iar1h0.betaxr.se[1]
      ibetaxrar1h0p =  iar1h0.betaxr.p[1]

    else
      betaxrmbbp=missing;
      betaxrmbbz=missing;

      betaxrar1h0p=missing;
      betaxrar1h0z=missing;

      ibetaxrar1h0p=missing;
      ibetaxrar1h0z=missing;
    end

    push!(results,(;
      label=:predict,
      measure=mcol,
      returnwindow=rcol,
      controls = join(string.(ccols), ", "),
      betaxr, betaxrmw, betaxrhac, betaxrmbbz,  betaxrmbbp,
      betaxrar1h0z,  betaxrar1h0p, ibetaxrar1h0z,  ibetaxrar1h0p,
      N))
  end

  return  results

  #form the prediction pairs
  #first match on LP

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using GLFW

# Device callbacks
GLFW.SetJoystickCallback((joystick, event) -> println("$joystick $event"))
GLFW.SetMonitorCallback((monitor, event) -> println("$monitor $event"))

window = GLFW.CreateWindow(640, 480, "GLFW Callback Test")
GLFW.MakeContextCurrent(window)

# Window callbacks
GLFW.SetWindowCloseCallback(window, (_) -> println("window closed"))
GLFW.SetWindowPosCallback(window, (_, x, y) -> println("window position: $x, $y"))
GLFW.SetWindowSizeCallback(window, (_, w, h) -> println("window size: $h x $w"))
GLFW.SetWindowFocusCallback(window, (_, focused) -> println("window focus: $focused"))
GLFW.SetWindowIconifyCallback(window, (_, iconified) -> println("window iconify: $iconified"))
GLFW.SetWindowMaximizeCallback(window, (_, maximized) -> println("window maximized: $maximized"))

# Input callbacks
GLFW.SetKeyCallback(window, (_, key, scancode, action, mods) -> begin
	name = GLFW.GetKeyName(key, scancode)
	if name == nothing
		println("scancode $scancode ", action)
	else
		println("key $name ", action)
	end
end)

GLFW.SetCharModsCallback(window, (_, c, mods) -> println("char: $c, mods: $mods"))
GLFW.SetMouseButtonCallback(window, (_, button, action, mods) -> println("$button $action"))
GLFW.SetCursorPosCallback(window, (_, x, y) -> println("cursor: $x, $y"))
GLFW.SetScrollCallback(window, (_, xoff, yoff) -> println("scroll: $xoff, $yoff"))

GLFW.SetDropCallback(window, (_, paths) -> println(paths))

glClear() = ccall(@eval(GLFW.GetProcAddress("glClear")), Cvoid, (Cuint,), 0x00004000)

GC.gc() # Force garbage collection so that improper reference management is more apparent via crashes

try
	while !GLFW.WindowShouldClose(window)
		glClear()
		GLFW.SwapBuffers(window)
		GLFW.WaitEvents()
	end
finally
	GLFW.DestroyWindow(window)
end

max_stars_repo_name: leroy-walton/playgrounds

max_issues_repo_name: leroy-walton/playgrounds

max_forks_repo_name: leroy-walton/playgrounds


text:
using BinaryBuilder

name = "minimap2"
version = v"2.22"

sources = [
    GitSource("https://github.com/lh3/minimap2", "7358a1ead1adfa89a2d3d0e72ffddd05732f9850"),
    DirectorySource("./static")
]

script = raw"""
cd $WORKSPACE/srcdir/minimap2
sed 's/SO/${dlext}/g' $WORKSPACE/srcdir/patches/Makefile.diff | patch Makefile
INC=${includedir} make libminimap2.${dlext}
cp libminimap2.${dlext} ${prefix}/lib
"""

platforms = filter(supported_platforms()) do platform
    (Sys.islinux(platform) | Sys.isapple(platform)) && platform.tags["arch"] == "x86_64"
end

products = [
    LibraryProduct("libminimap2", :minimap2),
]

dependencies = [
    Dependency("Zlib_jll"),
]

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

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Description
#
#   Pre-defined formatters.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export ft_printf, ft_round

"""
    function ft_printf(ftv_str, [columns])

Apply the formats `ftv_str` (see the function `sprintf1` of the package
**Formatting.jl**) to the elements in the columns `columns`.

If `ftv_str` is a `Vector`, then `columns` must be also be a `Vector` with the
same number of elements. If `ftv_str` is a `String`, and `columns` is not
specified (or is empty), then the format will be applied to the entire table.
Otherwise, if `ftv_str` is a `String` and `columns` is a `Vector`, then the
format will be applied only to the columns in `columns`.

# Remarks

This formatter will be applied only to the cells that are of type `Number`.

"""
ft_printf(ftv_str::String) = ft_printf([ftv_str])
ft_printf(ftv_str::String, column::Int) = ft_printf(ftv_str, [column])
ft_printf(ftv_str::String, columns::AbstractVector{Int}) =
    ft_printf([ftv_str for i = 1:length(columns)], columns)

function ft_printf(ftv_str::Vector{String}, columns::AbstractVector{Int} = Int[])
    lc = length(columns)

    lc == 0 && (length(ftv_str) != 1) &&
    error("If columns is empty, then ftv_str must have only one element.")

    lc > 0 && (length(ftv_str) != lc) &&
    error("The vector columns must have the same number of elements of the vector ftv_str.")

    # By using the function `sprintf1` from the package Formatting.jl, it was
    # possible to reduce the time to print a 100x5 table of `Float64`s from 25s
    # to 0.04s. Thanks to @RalphAS for the tip!

    if lc == 0
        return Dict{Int,Function}(0 => (v,i) -> begin
            return typeof(v) <: Number ? sprintf1(ftv_str[1], v) : v
        end)
    else
        ft = Dict{Int,Function}()

        for i = 1:length(columns)
            push!(ft, columns[i] => (v,j) -> begin
                return typeof(v) <: Number ? sprintf1(ftv_str[i], v) : v
            end)
        end

        return ft
    end
end

"""
    function ft_round(digits, [columns])

Round the elements in the columns `columns` to the number of digits in `digits`.

If `digits` is a `Vector`, then `columns` must be also be a `Vector` with the
same number of elements. If `digits` is a `Number`, and `columns` is not
specified (or is empty), then the rounding will be applied to the entire table.
Otherwise, if `digits` is a `Number` and `columns` is a `Vector`, then the
elements in the columns `columns` will be rounded to the number of digits
`digits`.

"""
ft_round(digits::Int) = ft_round([digits])
ft_round(digits::Int, columns::AbstractVector{Int}) =
    ft_round([digits for i = 1:length(columns)], columns)

function ft_round(digits::AbstractVector{Int}, columns::AbstractVector{Int} = Int[])
    lc = length(columns)

    lc == 0 && (length(digits) != 1) &&
    error("If columns is empty, then digits must have only one element.")

    lc > 0 && (length(digits) != lc) &&
    error("The vector columns must have the same number of elements of the vector digits.")

    if lc == 0
        return Dict{Int,Function}(0 => (v,i) -> begin
                                      if applicable(round,v)
                                          return round(v; digits = digits[1])
                                      else
                                          return v
                                      end
                                  end)
    else
        ft = Dict{Int,Function}()

        for i = 1:length(columns)
            push!(ft, columns[i] => (v,j) -> begin
                      if applicable(round,v)
                          return round(v; digits = digits[i])
                      else
                          return v
                      end
                  end)
        end

        return ft
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "mdoodz"
version = v"0.5.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/tduretz/MDOODZ7.0.git", "2e18a7685a8e7ef79e910d3461b317d52b8ff079")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd MDOODZ*/
rm makefile
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DOMP=ON -DOPT=ON -DJULIA=ON
make
make install
"""

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

# due to hdf5 dependency
# also due to Find scripts are not working on windows
platforms = [
    Platform("x86_64", "linux"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("x86_64", "macos"),
    Platform("aarch64", "macos"),
    Platform("x86_64", "windows"),
    Platform("i686", "windows"),
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libmdoodz", :libmdoodz)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CXSparse_jll", uuid="c77e7b6a-7cf9-58ed-a396-e1da12b05d87")),
    Dependency(PackageSpec(name="HDF5_jll", uuid="0234f1f7-429e-5d53-9886-15a909be8d59")),
    Dependency(PackageSpec(name="SuiteSparse_jll", uuid="bea87d4a-7f5b-5778-9afe-8cc45184846c")),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"); platforms=filter(!Sys.isbsd, platforms)),
    Dependency(PackageSpec(name="LLVMOpenMP_jll", uuid="1d63c593-3942-5779-bab2-d838dc0a180e"); platforms=filter(Sys.isbsd, platforms)),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
module BindingsGenerator

using Clang
using Clang_jll

const LIBCLANG_INCLUDE = joinpath(dirname(Clang_jll.libclang_path), "..", "include", "clang-c") |> normpath
const MLIR_C_INCLUDE = joinpath(@__DIR__, "..", "julia/usr/include/mlir-c") |> normpath
const MLIR_C_HEADERS = [joinpath(MLIR_C_INCLUDE, header) for header in readdir(MLIR_C_INCLUDE) if endswith(header, ".h")]

println("Found headers:")
display(MLIR_C_HEADERS)
try
    mkdir("lib")
catch e
    println("/lib already exists. Writing into /lib.")
end

wc = init(; headers = MLIR_C_HEADERS,
          output_file = joinpath(@__DIR__, "lib/libmlir_api.jl"),
          common_file = joinpath(@__DIR__, "lib/libmlir_common.jl"),
          clang_includes = vcat(LIBCLANG_INCLUDE, CLANG_INCLUDE),
          clang_args = ["-I", joinpath(LIBCLANG_INCLUDE, ".."),
                        "-I", joinpath(MLIR_C_INCLUDE, "..")],
          header_wrapped = (root, current) -> root == current,
          header_library = x -> "libmlir",
          clang_diagnostics = false,
         )

run(wc)

end # module


text:
"""
    Origin{T} <: Node{T}

    Global reference frame

    id: always 0
    name: always :origin
    state: State; defaults to zero values
    shape: Shape; defaults to EmptyShape
"""
mutable struct Origin{T} <: Node{T}
    id::Int64
    name::Symbol
    state::State{T}
    shape::Shape{T}

    function Origin{T}(;
        name::Symbol=:origin,
        state=State{T}(),
        shape::Shape=EmptyShape()) where T
        new{T}(getGlobalID(), name, state, shape)
    end

    function Origin(body::Body{T}) where T
        new{T}(body.id, body.name, body.state, body.shape)
    end
end


text:
using PumasTutorials
PumasTutorials.weave_file(".","test.jmd")