text:
# Squares building

W,(WW,EW,FW) = Lar.cuboid([1,1], true)
square = (W,EW)

squares = Lar.Struct([  square,  Lar.t(-1, -1),  Lar.s(3, 3),  square  ])

V,EV = Lar.struct2lar(squares)

model = CAGD.Model(V);
CAGD.addModelCells!(model, 1, EV);

# Lines adding
lines = CAGD.Model([
    -2.0  1.0  0.0  3.0
     0.0  3.0  3.0  0.0
]);
CAGD.addModelCells!(lines, 1, [[1, 2], [3, 4]]);

# Merge the models
CAGD.uniteModels!(model, lines)

# Planar Arrangment

A_model = CAGD.planar_arrangement(model, sparse(Int8[0,0,0,0,1,1,1,1,0,0]));

A_model.G
Matrix(A_model.T[1])
Matrix(A_model.T[2])

# Display Models

if todisplay
    VV = [[k] for k=1:size(A_model.G ,2)];
    triangulated_faces = CAGD.triangulate2D(A_model);
    FVs = convert(Array{Lar.Cells}, triangulated_faces);

    GL.VIEW(GL.GLExplode(A_model.G, FVs, 1.0, 1.0, 1.0, 99, 1));
    GL.VIEW(GL.GLExplode(A_model.G, FVs, 1.2, 1.2, 1.2, 99, 1));
end;


text:
module CardgenController
using Genie, Genie.Router, Genie.Requests, Genie.Renderer
using Cards, CardUtils, Verhoeff
using Random

const round_colors = ["maroon", "navy", "DarkMagenta"]

autoinc() = let i=-1; () -> i+=1 end
numgen = autoinc()

wrap_contenttype(data, type) = respond(WebRenderable(
    body=data,
    status=200,
    headers=Dict("Content-Type" => type)))

gencard_ps(id, round) = mk_card(gen_rand_array(MersenneTwister(stripcheck(id))),
  round_colors[round], num=id, round=round)

gencard_pdf = ps2pdf ∘ gencard_ps

gencard_rg(gc, type) = () -> verhoeff_verify(payload(:id)) ?
  wrap_contenttype(gc(payload(:id), payload(:round)), type) :
  Router.error("This card's check digit is bad; it can't possibly exist",
               "text/html", Val(404))

gencard_ps_r    = gencard_rg(gencard_ps,  "application/postscript")
gencard_pdf_r   = gencard_rg(gencard_pdf, "application/pdf")

genseq_rg(sym)   = () -> redirect(sym,
                                  id=verhoeff_gencheck(numgen()),
                                  round=payload(:round))

addcheck_rg(sym) = () -> redirect(sym,
                                  id=verhoeff_gencheck(payload(:id)),
                                  round=payload(:round))
end

ext: jl

max_stars_repo_name: tuckertwo/ECbingo

max_issues_repo_name: tuckertwo/ECbingo

max_forks_repo_name: tuckertwo/ECbingo


text:
function _count_α(α)
    return (size(α)[1]-1)*size(α)[2]
end

function _count_params(model)
    return sum(length.(params.(model)))
end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
using Documenter, Fishmetica

makedocs(modules=[Fishmetica],
        doctest=true)

deploydocs(deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/picestode/Fishmetica.git",
    julia  = "0.6.2",
    osname = "linux")

max_stars_repo_name: picestcode/Fishmetica

max_issues_repo_name: picestcode/Fishmetica

max_forks_repo_name: picestcode/Fishmetica


text:
# This file contains the functions for heat transfer between plane geometries

" Far-field heat transfer between two semi-infinite media,classical approximation "
function farfield_transfer(em1 , em2, T1, T2)

    return sigma*(T2^4-T1^4)/(1.0/em1 + 1.0/em2 - 1.0)

end

"""
    incoherent_total_heat_transfer_w(b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer, T1,T2,w;tolkx=1e-6)

compute the total monocromatic radiative heat flux between planar surfaces in the far field,
considereing total incoherent radiation. This should only be considered correct (less than 5 % error) when the separation
distance between both objects is larger than 100*λ_wien.

#Arguments:
- `b1 :: LayerOrMultiLayer`: Multilayer description of body 1
- `b2 :: LayerOrMultiLayer`: Multilayer description of body 2
- `gap :: Layer`: Material and thickness of gap separation
- `T1`: temperature of body 1
- `T2`: temperature of body 2
- `w` : frequency in radHz
- `tolkx = 1e-6`: relative error on integration over kx
"""
function incoherent_total_heat_transfer_w(b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer, T1,T2,w;tolkx=1e-6)
   T1==0.0 ? l1=0.0 : l1 = lambda_wien(T1)
   T2==0.0 ? l2=0.0 : l2 = lambda_wien(T2)

   @assert gap.thickness >= 100*max(l1,l2)  "the gap size should at least be 100 times Wien's wavelength"

    em1(kx) = emissivity_kx_w([Layer(gap.material) ; b1], kx, w)
    em2(kx) = emissivity_kx_w([Layer(gap.material) ; b2], kx, w)

    tn(kx)  = kx*em1(kx)*em2(kx)*0.25
    td(kx)  = 1 - ( 1- em1(kx)*0.5)*( 1- em2(kx)*0.5)
    t(kx)   = tn(kx)/td(kx)
    t2(u)   = t(u*w/c0)

    (val,err) = quadgk(t2, 0,1; rtol=tolkx)
    return 2/(w/c0)*val*(planck(w,T1) - planck(w,T2))
end

"""
    incoherent_total_heat_transfer(b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer, T1, T2, w1, w2; tolkx=1e-6)

compute the total radiative heat flux between planar surfaces in the far field,
considereing total incoherent radiation. This should only be considered correct (less than 5 % error) when the separation
distance between both objects is larger than 100*λ_wien.

#Arguments:
- `b1 :: LayerOrMultiLayer`: Multilayer description of body 1
- `b2 :: LayerOrMultiLayer`: Multilayer description of body 2
- `gap :: Layer`: Material and thickness of gap separation
- `T1`: temperature of body 1
- `T2`: temperature of body 2
- `w` : frequency in radHz
- `tolkx = 1e-6`: relative error on integration over kx
- `tolw = 1e-6`: relative error on integration over the frequency w
"""
function incoherent_total_heat_transfer(b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer, T1,T2,w1,w2;tolkx=1e-6, tolw = 1e-6)
    u1 = w1*ħ/kb
    u2 = w2*ħ/kb
    ht(w) = incoherent_total_heat_transfer_w(b1 ,b2, gap , T1,T2, w ;tolkx=tolkx)
    ht2(u) = ht(u*kb/ħ)
    (val,err) = quadgk(ht2, u1 , u2 ; rtol=tolw)
    return val*kb/ħ
end

"""
    incoherent_total_heat_transfer(b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer, T1, T2; tolkx=1e-6)

compute integration over ω from 0 to +∞
"""
function incoherent_total_heat_transfer(b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer, T1,T2;tolkx=1e-6, tolw = 1e-6)
    ht(w) = incoherent_total_heat_transfer_w(b1 ,b2, gap , T1,T2, w ;tolkx=tolkx)
    ht2(u) = ht(u*kb/ħ)
    ht3(t) = ht2(t/(1-t))/(1-t)^2
    (val,err) = quadgk(ht3, 0 , 1 ; rtol=tolw)
    return val*kb/ħ
end

" Evanescent contribution to radiative heat transfer "
function transmission_kx_w(:: Evanescent ,b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer ,pol :: Polarization , kx ,w)

    gap2 = Layer(gap.material,0.0)
    b1   = prepend!([b1;],[gap2])
    b2   = prepend!([b2;],[gap2])

    (r_21,t)  = rt(b1,pol,kx,w) :: Tuple{Complex{Float64},Complex{Float64}}
    (r_23,t)  = rt(b2,pol,kx,w) :: Tuple{Complex{Float64},Complex{Float64}}

    k2z = compute_kz(kx,permittivity(gap.material,w),w) :: Complex{Float64}
    exp_val1 = exp(-2.0*imag(k2z)*gap.thickness)         :: Float64
    exp_val2 = exp(2.0*im*k2z*gap.thickness)            :: Complex{Float64}

    return 4.0*exp_val1*imag(r_21)*imag(r_23)/abs(1.0-r_21*r_23*exp_val2)^2

end

" Propagative contribution to radiative heat transfer "
function transmission_kx_w(:: Propagative, b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer ,pol :: Polarization , kx ,w)

    gap2 = Layer(gap.material,0.0)
    b1 = prepend!([b1;],[gap2])
    b2 = prepend!([b2;],[gap2])

    T_21      = 0.0
    T_23      = 0.0
    (r_21,t_21)  = rt(b1,pol,kx,w) :: Tuple{Complex{Float64},Complex{Float64}}
    (r_23,t_23)  = rt(b2,pol,kx,w) :: Tuple{Complex{Float64},Complex{Float64}}

    k2z = compute_kz(kx,permittivity(gap.material,w),w)  :: Complex{Float64}
    exp_val2 = exp(2.0*im*k2z*gap.thickness)             :: Complex{Float64}

    if imag(permittivity(substrate(b1),w)) == 0.0
        eps1,k1z,eps2,k2z = outer_media(b1,kx,w)
        T_21 = power_t(pol,eps1,eps2, k1z,k2z,t_21)
    end
    if imag(permittivity(substrate(b2),w)) == 0.0
        eps1,k1z,eps2,k2z = outer_media(b2,kx,w)
        T_23 = power_t(pol,eps1,eps2, k1z,k2z,t_23)
    end

    return (1.0 - abs(r_21)^2 - T_21)*(1.0 - abs(r_23)^2 - T_23)/abs(1.0-r_21*r_23*exp_val2)^2

end

" Monocromatic transmission "
function transmission_w(field :: Evanescent ,b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer ,pol :: Polarization ,w,tol)

    t_kx_w(kx)  = kx*transmission_kx_w(field ,b1,b2,gap,pol,kx,w)
    t2_kx_w(u)  = t_kx_w(u*w/c0)
    t3_kx_w(v)  = t2_kx_w(1 + v/(1-v))/(1-v)^2

   (val,err) = quadgk(t3_kx_w, 0,1; rtol=tol)
    return val*(w/c0)
end

function transmission_w(field :: Propagative ,b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer ,pol :: Polarization ,w,tol)
    t_kx_w(kx)  = kx*transmission_kx_w(field ,b1,b2,gap,pol,kx,w)
    t2_kx_w(u)  = t_kx_w(u*w/c0)

    (val,err) = quadgk(t2_kx_w, 0,1; rtol=tol)
    return val*(w/c0)
end

" Monocromatic heat flux "
function heat_flux_w(field :: TotalField ,b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer ,pol :: Polarization ,w, T;tol=1e-2)
    t_w  = transmission_w(field,b1,b2,gap,pol,w,tol)
    be_w = bose_einstein(w,T)
    return be_w*t_w/4.0/pi^2
end

" Heat flux "
function heat_flux(field :: TotalField ,b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer ,pol :: Polarization , T; w1=1e13,w2=1e15, tol1=1e-8 , tol2=1e-8)
    q_w(w)  = heat_flux_w(field ,b1 , b2 , gap  ,pol , w ,T ; tol=tol1)
    q1_w(u) = q_w(u*kb*T/ħ)
    u1=w1*ħ/kb/T
    u2=w2*ħ/kb/T

    val :: Float64  = 0.0
    err :: Float64  = 0.0
    (val,err) = quadgk(q1_w, u1 , u2 ; rtol=tol2)

    return val*(kb*T)/ħ
end

" Monocromatic net heat transfer "
function heat_transfer_w(field :: TotalField ,b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer ,pol :: Polarization ,w, T1,T2;toler=1e-2)
    return (bose_einstein(w,T1)-bose_einstein(w,T2))*transmission_w(field,b1,b2,gap,pol,w,toler)/4.0/pi^2
end

" Net heat transfer "
function heat_transfer(field :: TotalField ,b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer ,pol :: Polarization, T1,T2 ; tol1=1e-8,tol2=1e-8)
    return heat_flux(field,b1,b2,gap,pol,T2 ; tol1=tol1 , tol2=tol2) - heat_flux(field,b1,b2,gap,pol,T1 ; tol1=tol1 , tol2=tol2)
end

function total_heat_transfer_w(b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer, T1,T2,w;tolkx=1e-6,tolw=1e-6)
    valt = 0.0
    cnt  = 0
    q    = zeros(Float64,5)
    for f in ( Evanescent(), Propagative())
        for p in (te(),tm())
            cnt += 1
            valt  += heat_transfer_w(f ,b1 ,b2, gap ,p ,w, T1,T2;toler=tolkx)
         end
    end
    return valt
end

"Total net heat transfer"
function total_heat_transfer(b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer, T1,T2;tolkx=1e-6,tolw=1e-6)
    valt = 0.0
    cnt  = 0
    q    = zeros(Float64,5)
    for f in ( Evanescent(), Propagative())
        for p in (te(),tm())
            cnt  += 1
            ht(w) = heat_transfer_w(f ,b1 ,b2, gap ,p ,w, T1,T2;toler=tolkx)
            ht2(u) = ht(u*kb/ħ)
            ht3(t) = ht2(t/(1-t))/(1-t)^2
            (val,err) = quadgk(ht3, 0 , 1 , rtol=tolw)
            valt  += val
            q[cnt] = val*kb/ħ
         end
    end
    q[5]=valt*kb/ħ

    return q
end

function total_heat_transfer(b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer, T1,T2,w1,w2;tolkx=1e-6,tolw=1e-6)
    valt = 0.0
    cnt  = 0
    q    = zeros(Float64,5)
    u1 = w1*ħ/kb
    u2 = w2*ħ/kb
    for f in (Evanescent(),Propagative())
        for p in (te(),tm())
            cnt  += 1
            ht(w) = heat_transfer_w(f ,b1 ,b2, gap ,p ,w, T1,T2;toler=tolkx)
            ht2(u) = ht(u*kb/ħ)
            (val,err) = quadgk(ht2, u1 , u2 ; rtol=tolw)
            valt  += val
            q[cnt] = val*kb/ħ
         end
    end
    q[5]=valt*kb/ħ

    return q
end

function total_heat_transfer_double(b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer, T1,T2,w1,w2;tol=1e-6)
    valt = 0.0
    cnt  = 0
    q    = zeros(Float64,5)
    u1 = w1*ħ/kb
    u2 = w2*ħ/kb
    for f in (Evanescent(), Propagative())
        for p in (te(),tm())
            cnt  += 1
            val = integrand_double(f, p, b1, b2, gap, T1,T2,w1,w2;tol=tol)
            valt  += val
            q[cnt] = val*(kb/ħ)^3/c0^2/4.0/pi^2
         end
    end
    q[5]=valt*(kb/ħ)^3/c0^2/4.0/pi^2

    return q
end

function integrand_double(field :: Propagative, pol :: Polarization, b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer, T1,T2,w1,w2;tol=1e-6)
    u1 = w1*ħ/kb
    u2 = w2*ħ/kb
    integr(x) = (bose_einstein(x[1]*kb/ħ,T1) - bose_einstein(x[1]*kb/ħ,T2))*x[2]*x[1]^2*transmission_kx_w(field ,b1,b2,gap,pol,x[2]*x[1]*kb/ħ/c0,x[1]*kb/ħ)
    (val,err) = hcubature(integr,[u1,0.0],[u2,1.0]; rtol=tol)
    return val
end

function integrand_double(field :: Evanescent, pol :: Polarization , b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer, T1,T2,w1,w2;tol=1e-6)
    u1 = w1*ħ/kb
    u2 = w2*ħ/kb
    integr(x)  = (bose_einstein(x[1]*kb/ħ,T1) - bose_einstein(x[1]*kb/ħ,T2))*x[2]*x[1]^2*transmission_kx_w(field ,b1,b2,gap,pol,x[2]*x[1]*kb/ħ/c0,x[1]*kb/ħ)
    integr2(y) = integr([y[1],1 + y[2]/(1-y[2])])/(1-y[2])^2
    (val,err)  = hcubature(integr2,[u1,0.0],[u2,1.0]; rtol=tol)
    return val
end

function total_transmission_kx_w(b1 :: LayerOrMultiLayer, b2 :: LayerOrMultiLayer, gap :: Layer, kx, w)
    valt = 0.0
    for p in (te(),tm())
        if kx <= w/c0
            valt += transmission_kx_w(Propagative(), b1, b2 , gap  ,p , kx ,w)
        else
            valt += transmission_kx_w(Evanescent(), b1, b2 , gap  ,p , kx ,w)
        end
    end
    return valt
end

function total_transmission_map(b1 :: LayerOrMultiLayer,
                                b2 :: LayerOrMultiLayer,
                                gap :: Layer,
                                kx :: AbstractArray, w :: AbstractArray)

    t = zeros(length(w),length(kx))
    for j = 1:length(kx)
        for i = 1:length(w)
            t[i,j] = total_transmission_kx_w(b1,b2,gap,kx[j],w[i])
        end
    end
    return t
end


text:
function fn()
    return Array{Float32}([1,2,3,4,5])
end

max_stars_repo_name: paulmelis/blender-julia-test

max_issues_repo_name: paulmelis/blender-julia-test

max_forks_repo_name: paulmelis/blender-julia-test


text:
# Email from Jean-Michel:
#
# I made a little tar file (-czf) of cs32 grid files (the same as in
# https://github.com/MITgcm/MITgcm/blob/master/verification/tutorial_held_suarez_cs/input
# with a litte matlab script that load them.
#
# Each file (big-endian, 64 bits) contains 18 reccords of 33 x 33 values, with:
# rec  | MITgcm | Description
#   #  |  name  |
#   1  : xC    <->   longitude (in -180:+180) at grid-cell center (cca)
#   2  : yC    <->   latitude  at grid-cell center (cca)
#   3  : dxF   <->   dx (in m) at center (cca)
#   4  : dyF   <->   dy (in m) at center (cca)
#   5  : rAc   <->   Area (in m^2) at center ( zA^cca )
#   6  : xG    <->   longitude (in -180:+180) at grid-cell corner (ffa)
#   7  : yG    <->   latitude  at grid-cell corner (ffa)
#   8  : dxV   <->   dx (in m) at grid-cell corner (ffa)
#   9  : dyU   <->   dy (in m) at grid-cell corner (ffa)
#  10  : rAz   <->   Area (in m^2) at grid-cell corner  ( zA^ffa )
#  11  : dxC   <->   dx (in m) at U pt (fca)
#  12  : dyC   <->   dy (in m) at V pt (cfa)
#  13  : rAw   <->   Area (in m^2) at U pt ( zA^fca )
#  14  : rAs   <->   Area (in m^2) at V pt ( zA^cfa )
#  15  : dxG   <->   dx (in m) at V pt (cfa)
#  16  : dyG   <->   dy (in m) at U pt (fca)
#  17  : cosine of the local angle between the grid orientation (e.g., U(1) velocity) and West-East direction at grid-cell center
#  18  : sine of the local angle between the grid orientation (e.g., U(1) velocity) and West-East direction at grid-cell center
#
# For all cca fields (1 to 5), the last row and last column (i.e., [33, :] and [:, 33]) are
# not used and are just zero.
# And for fields 11, 13, 16 : last column [:,33] is just zero.
# And for fields 12, 14, 15 : last row [33, :] is just zero.

using JLD2

const N_records = 18

function read_cubed_sphere_face(filepath)
    binary_data = read(filepath)

    # Calculate the number of Float64 elements in the binary file.
    N_numbers = Int(length(binary_data) / 8)

    # Calculate the size of the grid (number of faces)
    N = Nx = Ny = √(N_numbers / N_records) |> Int

    # Load binary data into a big array.
    cubed_sphere_face_data = Array{Float64}(undef, N_numbers)
    read!(filepath, cubed_sphere_face_data)

    # Reverse byte order since binary data is big-endian.
    cubed_sphere_face_data = bswap.(cubed_sphere_face_data)

    ## Move data into 2D arrays/fields.

    inds(n, N) = UnitRange((n-1) * N^2 + 1, n * N^2)
    contents(bin, n, N) = reshape(bin[inds(n, N)], N, N)

    cubed_sphere_variables = (
          λᶜᶜᵃ = contents(cubed_sphere_face_data,  1, N),
          φᶜᶜᵃ = contents(cubed_sphere_face_data,  2, N),
         Δxᶜᶜᵃ = contents(cubed_sphere_face_data,  3, N),
         Δyᶜᶜᵃ = contents(cubed_sphere_face_data,  4, N),
         Azᶜᶜᵃ = contents(cubed_sphere_face_data,  5, N),
          λᶠᶠᵃ = contents(cubed_sphere_face_data,  6, N),
          φᶠᶠᵃ = contents(cubed_sphere_face_data,  7, N),
         Δxᶠᶠᵃ = contents(cubed_sphere_face_data,  8, N),
         Δyᶠᶠᵃ = contents(cubed_sphere_face_data,  9, N),
         Azᶠᶠᵃ = contents(cubed_sphere_face_data, 10, N),
         Δxᶠᶜᵃ = contents(cubed_sphere_face_data, 11, N),
         Δyᶜᶠᵃ = contents(cubed_sphere_face_data, 12, N),
         Azᶠᶜᵃ = contents(cubed_sphere_face_data, 13, N),
         Azᶜᶠᵃ = contents(cubed_sphere_face_data, 14, N),
         Δxᶜᶠᵃ = contents(cubed_sphere_face_data, 15, N),
         Δyᶠᶜᵃ = contents(cubed_sphere_face_data, 16, N),
        cosᶜᶜᵃ = contents(cubed_sphere_face_data, 17, N),
        sinᶜᶜᵃ = contents(cubed_sphere_face_data, 18, N)
    )

    return cubed_sphere_variables
end

#####
##### Convert grid_cs32.face*.bin files to JLD2
#####

cubed_sphere_filepath = "cubed_sphere_32_grid.jld2"

jldopen(cubed_sphere_filepath, "w") do file
    for (face_number, face) in enumerate(["001", "002", "003", "004", "005", "006"])
        filepath = "grid_cs32.face$face.bin"
        @info "Converting $filepath..."

        cubed_sphere_vars = read_cubed_sphere_face(filepath)

        for (var_name, var) in pairs(cubed_sphere_vars)
            file["face$face_number/$var_name"] = var
        end
    end
end

@info "Generated $cubed_sphere_filepath"


text:
using StochasticArithmetic
using StochasticArithmetic.EFT
using Test, Jive
using Formatting
using LinearAlgebra
using Random

Random.seed!(42)

@onlyonce begin
    const x1 = [-6.340663515282668472e-01, -1.032478304663339008e+00, -9.398547076618038787e+00, 8.775666080119043144e+00,
                8.154021046137667206e-01, -2.872197835897075890e+00, -4.130201373187722957e+00, -3.800523733887637978e-01,
                3.767316625622183501e+00, -7.693130974005489620e+00, 3.789995392417770503e+00, 2.332943819408263633e-02,
                1.117290458364406547e+00, -1.029424506285847585e+01, -1.972399240364107298e+00, 2.409324890261351548e+00,
                1.514218404563302478e+01, -9.803395024842704863e-01, 3.624014860482689393e+00, 3.454279166512129517e+00,
                2.150407406181448877e+00, 8.380973701960218181e-01, -4.718497388831256245e+00, 1.092357091710888817e+01,
                6.201621778907238092e-01, -2.431850975835210704e+00, -9.587862761576261272e-01, -4.744452383420612485e+00,
                -1.928632137138599045e+00, -1.514994713586201414e+00, -2.183297984924969626e+00, 1.424576004184546996e+01,
                9.692913616904123231e+00, -1.547104695550167275e+00, 4.204554274223821331e+00, 9.347482098474469980e-01,
                7.399023959120087923e-02, -4.191315817777256925e+00, 8.475563783170365184e+00, -1.721899852461975300e+00,
                -7.486679611973100279e+00, -5.787527477378286850e-02, 5.560569808434617123e+00, 9.797972685704890716e-02,
                -1.002061593959254537e+00, -3.012224958934606622e+00, -8.084325856937959465e-01,-1.905132695346669347e-01,
                5.188198608642882181e-01, -7.601644981994752914e-01, -2.451000778662956403e+00, -9.693293689094003573e-01,
                2.882117571206391560e+00, 2.269843085532857518e+00, 2.958122951288597147e+00, 2.630063037097047030e+00,
                -6.039968388496026463e+00, -4.780642850877452332e+00, 2.900723065170244808e-01, -1.902925253860966448e+00,
                -1.317310182590558298e+01, 3.598067326176326297e-01, 9.089778308264061524e-01, -3.252711812214260867e+00,
                -2.310740125115115173e+00, -1.617256866537321480e+00, 2.935716535944118633e+00, -3.234285498797106584e+00,
                3.237522392256278359e+00, 4.234360820593143337e+00, -8.323039886529393083e-02, 3.392892532427555174e+00,
                4.675163615855799232e-02, -2.586710584045738504e+00, 7.263314423520003116e-01, -1.805428377059409950e+00,
                5.540332712174778074e-01, -1.680124146046180877e+00, 3.384520587773187650e-01, -9.285586956857446950e+00,
                1.226698449940088764e+01, 1.955258288381735765e+00, -1.751307418060096488e+00, 2.822752021460029237e+00,
                3.284769252355420832e+00, -1.435237503511320289e+00, -1.640731710755447281e+01, -5.073426138128573903e+00,
                2.769683815669464200e-01, -2.364893385542928161e+00, -2.852226279751059579e+00, 9.442997895825611110e-02,
                1.613660165004466629e+00, -7.535249859542684225e-01, 1.875042460291437951e+00, -1.204955334660282817e+00,
                -1.312339352451679320e+00, 1.029598837959142621e+01, 1.181600593671352062e+00, -1.298538943700622506e+00]

    const y1 = [2.904271758197305431e+00, -7.586923460587383872e+00, -2.433948194217723793e-01, 7.211094113495739588e-01,
                3.682896270200847866e+00, -1.300501673307355999e+00, -1.156385074689293901e+01, 3.527526896240007925e+00,
                6.138604794505512219e-01, 3.789599979884540115e+00, 4.964840840608681027e+00, 1.003350838878128926e+00,
                4.318670472060466281e+00, 1.014498387743745011e+00, -2.234343190283604841e+00, -6.082736396317145466e-01,
                2.676279050222233336e-01, -2.474834547943665974e+00, -1.502899068454430775e+00, -1.853840253722779163e+00,
                -2.690201016281298152e+00, -2.762974865876188346e+00, -1.691254623016295255e+00, 7.037535285788961217e-03,
                8.664535245217100945e-01, -1.152357491478533280e+00, 1.028831455069202150e+00, -1.999749826926913443e-01,
                -4.083433509619152524e+00, 8.094079296833567305e-01, 1.548542269669347382e+00, 1.980881016229119940e+00,
                8.859099623710096072e-01, 1.071154721221595274e+01, -9.323551904078861141e-02, 1.071342262133416101e+00,
                -4.657941295025913031e-01, -1.417360280450415067e+00, 2.973792467793110017e+00, -5.333088668455037151e-01,
                3.506402422142353981e+00, -2.626404402475698419e+01, -5.821956212486361082e+00, -2.081525484012217042e+00,
                1.166141951259942111e+01, -2.097572645415919368e-01, -4.170975773332724401e+00, -1.088414007670470252e+00,
                1.437306693315153927e+00, 2.525907646776998305e+00, -3.752587514423713078e+00, 1.547393157729613611e+00,
                5.679913072629369508e-02, -4.458331660254577500e+00, -1.152536727954037943e+00, 6.959027574408890082e+00,
                -4.529530640152803933e+00, 2.850439979915360689e-01, 1.551842054260149303e-02, -2.499408640456699970e+00,
                8.025705275120998294e+00, -2.339259319371728374e+00, 1.690685407197963785e-01, 1.773766221659853759e+00,
                8.418125128249426270e-02, -1.240164986459572516e+00, 1.244333273760872993e+00, 2.229039874696919110e+00,
                1.578366388587486036e+01, 4.010353714690342741e+00, 4.066297377969993221e+01, -3.590142776920160483e-01,
                -5.317967102854119776e-01, -5.773767334043824473e-01, -1.684750638812332113e+00, 3.036554271013544026e-01,
                -1.915078321632401881e+00, 1.230841242088420717e+00, -1.360733255768588945e+00, -2.194897100045668736e+00,
                -3.147503241618139214e-01, -6.491785562479393867e+00, 8.732579404206023410e-01, -1.197064788553268588e+00,
                -1.051560432555771030e+00, 3.857770155380043642e-01, 1.780083978520605070e+00, 6.985067850240034293e-01,
                5.260718423653707454e+00, 2.393421981789398967e-01, 1.226572487102020936e+00, -4.390511466037616550e-01,
                4.857518183853492544e-01, 9.529898633767956539e+00, 1.640668058073694624e+01, 1.437731392011400100e+00,
                2.081762488795357413e+00, -1.677143583187016507e+00, 3.714458425965566635e-01, -1.788480931953309927e+00]

const x2 = [-4.703899089072115913e+01, -2.302208194432241726e+03, -1.564726108041984262e+05, 1.526300899208613373e+04,
            -2.132307857407913616e+02, 2.307518097958435537e+05, -5.192572714570803782e+01, 2.325599833211251966e+00,
            2.817494622670464651e+03, 1.708710020563955254e+02, 6.269522630918198303e+01, 3.687443650750188340e+01,
            1.038604765677104957e+01, -1.915975573744401117e+03, 1.168276776869518835e+04, 2.454491144842770009e+05,
            -5.592615927834867051e+01, 1.850307202352605032e+02, -3.838926631532003521e+05, 4.581024635636574203e+03,
            -4.418458866400893044e+02, -1.291545011255889767e+04, -2.048725390018565307e+01, 1.204863349219630436e+03,
            9.760401372457109259e-01, 1.410078033770042748e+04, -4.813837462145154859e+04, -1.437299278084431897e+01,
            -5.155141433885022707e+03, -2.255002085683215992e+05, 5.356346283817102085e+04, 3.667837584637414693e+04,
            -7.216672700119137880e+04, -1.259108454694217698e+03, 4.419830333598044092e+01, -3.585915213658949360e+05,
            3.816928321706653549e+01, 2.135141050205000557e+00, 1.174394551478064532e+05, -2.710374703222661275e+00,
            1.611291396027762385e+05, -1.357925813584457501e+03, -1.242968716437326115e+01, 7.973254168167831267e+03,
            -1.866706364323428602e+02, 2.042441529597737571e+01, 2.158858184523954282e+00, 2.467055660032614028e+01,
            5.375520883445242362e+04, -1.144212376192913938e+00, 5.994613017286798140e+02, 3.193496564759223055e+00,
            -4.148337081370820911e-01, -4.005110071830330076e-01, 1.993341984143252921e+04, -1.147244379776030421e+01,
            6.119929339861500672e-01, -6.092377012266528027e+02, 2.411059939457267774e+03, 9.559736934852315926e+01,
            -4.805764149222165997e+02, 4.009221133304038085e+04, 8.398867979989321775e-02, 5.743855858417433069e+01,
            -4.980311032420013362e+02, 1.509343053835127830e+04, 1.454679269769394365e-02, 2.854821043634199995e+01,
            -4.072883091111050291e-01, 5.108044955045782842e+01, -4.071539984133301004e+03, 8.751867981322608614e+01,
            -1.285130410507232591e+05, -3.354125317288757287e+03, 8.135505746482945000e-01, 4.723687613018444154e+01,
            1.189127895476845197e+03, -2.595180439442086026e+03, -3.753511234127828402e+02, 6.828869753373305151e+00,
            -2.424235269921588838e+02, -1.527378814146540753e+00, -2.753704293718250701e+05, 6.474506937411601548e+00,
            -1.077862672807065770e+06, -4.139018884967313352e+04, -1.418115653149918671e+03,-1.042964882907427171e+00,
            1.508774032010480548e+00, -1.001265318103480295e+00, 2.624649797909716753e+00, -1.368302828037557447e+04,
            -6.551342612757771349e+02, -1.855860091170598389e+04, 3.231487400789144715e+01, 9.365283668012376438e+03,
            1.078202387444531450e+04, 1.383253991429271628e+04, 1.406923835053669469e+00, -5.288835748811105963e-01]

const y2 = [-1.252871778887357568e+00, 1.078962521009007787e+03, -3.174894745849149331e+04, -1.530373455916264093e-02,
            4.784478272815556821e+00, -1.192086818719253206e+05, 3.668476944821762942e+01, 1.952404524146996501e+00,
            4.765260911366479377e+03, -3.289741608770166970e-01, -4.226193799391714911e+00, -1.824747767003190901e+01,
            1.352880351874526710e+00, 2.528843202632653629e+03, -1.591117252770609957e+04, -7.948333831513273440e-01,
            1.763730594001834406e+00, -4.142264664747529190e-01, 1.066477328704419342e+05, 3.159332045246057532e+03,
            -1.553414963997087783e+00, 6.473297499671324395e-01, -5.969956700965419749e+00, -1.432466768506058770e+03,
            3.846365749072568452e+00, 2.790519287535665249e+00, 2.712376392123455826e+00, 3.495804741318556119e-01,
            -1.976446781200929399e+04, 1.520383644699746801e+05, 2.197669377729083135e+04, -1.069924380149832541e+04,
            -1.364616481795411164e+04, 2.396787353506200924e-01, 1.576007089234049907e+01, -1.640112041994733736e+06,
            5.423184829884636571e+00, -2.854528605476235548e+00, 1.147783795765330961e+00, -1.108619081995681377e+00,
            -1.290633352716010995e+05, 8.452162075538846020e-01, -1.339879936618119416e+02, -1.094006209459130559e+03,
            2.012798587947189688e+00, -9.383522798994409220e-01, -1.090080505163032232e+00, 2.617736898339717300e+01,
            1.163689330823294954e+00, -1.997489359474004189e+01, -1.497673101858773350e+02, 2.910347668348660832e+01,
            6.358777943232427621e-01, -3.223860439180939785e+00, -1.020470547974752958e-02, -2.573297990895500220e+01,
            1.219805387956967024e+00, 8.361151282404189260e+02, 3.744876783652914298e+01, -2.384432846910274861e+01,
            4.299066016316168692e+00, -3.581774420101832220e+04, 6.850586911579179095e+00, 1.991718463645072712e-01,
            -7.417107650666663403e+02, -7.336543124226103828e+03, -5.102552005287688264e-01, 4.026602472160422508e+01,
            3.045402731323658330e+01, 6.217979244611119327e+01, 4.537584779899785644e+03, -6.591459843464705592e+00,
            3.949078543721616708e+00, -1.752538562659680388e-01, 5.175050101609877107e-01, 9.134566870551368822e+01,
            2.127631922755461868e+00, -3.851289922958672785e-01, -5.430969234878427443e+01, 2.599543161451561613e+00,
            2.418636072495988856e+02, 2.920924458066446405e+00, -3.498182839331332798e+04, 1.412711102346127312e+01,
            4.450907201869644923e+05, -4.828728797199047129e+00, 4.874121684200859050e+01, 6.351064829043707505e+00,
            -1.052720771587801307e+00, 4.862959846339910541e+01, -5.724214959936722025e-01, 1.159505461894495326e+00,
            1.220304567286626707e+02, -8.614815388556389397e-01, -2.020054259416569076e-01, 7.288236435882137130e+02,
            -4.939913027080400365e-01, 2.953432423721218220e+04, -5.565878606086428704e+00, 6.137681992750880777e-01]
end

@testset "StochasticArithmetic" begin
    @testset "EFT" begin

        @testset "twoSum" begin
            let (x,y) = twoSum(0.1, 0.3)
                @test x ≈  0.4
                @test y ≈ -2.7755575615628914e-17
            end

            let (x,y) = twoSum(0.1f0, 0.2f0)
                @test x ≈  0.3f0
                @test y ≈ -7.450581f-9
            end
        end

        @testset "twoProd" begin
            let (x, y) = twoProd(0.1, 5.)
                @test x ≈ 0.5
                @test y ≈ 2.7755575615628914e-17
            end

            let (x, y) = twoProd(0.1f0, 5.f0)
                @test x ≈ 0.5f0
                @test y ≈ 7.450581f-9
            end
        end

        @testset "twoDiv" begin
            let (x, y) = twoDiv(1., 3.)
                @test x ≈ 0.3333333333333333
                @test y ≈ 1.850371707708594e-17
            end

            let (x, y) = twoDiv(1.f0, 3.f0)
                @test x ≈  0.333333f0
                @test y ≈ -9.934108f-9
            end
        end

        @testset "split" begin
            function checkSplit(a)
                binRepr(x) = format("{:b}", Int(significand(x)*2^(53)))
                (x, y) = EFT.split(a)

                xTrim = replace(binRepr(x), r"0+$" => "")
                @test startswith(binRepr(a), xTrim)

                n = length(xTrim)
                aEnd = (binRepr(a))[n+1:end]
                aEnd = replace(aEnd, r"^0+" => "")
                aEnd = replace(aEnd, r"0+$" => "")
                yTrim = replace(binRepr(y), r"0+$" => "")
                @test aEnd == yTrim
            end

            checkSplit(0.1)
        end
    end

    @testset "rounding modes" begin
        @testset "directed" begin
            for (x,y) in [(0.1, 0.3), (-0.31, 0.0125)]
                @test +(DWN, x, y) < +(UP, x, y)
                @test *(DWN, x, y) < *(UP, x, y)
                @test /(DWN, x, y) < /(UP, x, y)
            end
        end

        @testset "random" begin
            for (x,y) in [(0.1, 0.3), (-0.31, 0.0125)]
                @test all(+(RND, x, y) in (+(UP, x, y), +(DWN, x, y)) for _ in 1:100)
                @test all(*(RND, x, y) in (*(UP, x, y), *(DWN, x, y)) for _ in 1:100)
                @test all(/(RND, x, y) in (/(UP, x, y), /(DWN, x, y)) for _ in 1:100)
            end
        end
    end

    for T in (SFloat64, SFloat32)
        @testset "$T" begin
            @testset "4 ops" begin
                let x = T(0.1) + T(0.3)
                    @test x isa T
                    @test value(x) ≈ 0.4
                end

                let x = T(0.1) - T(0.3)
                    @test x isa T
                    @test value(x) ≈ -0.2
                end

                let x = T(0.1) * T(0.3)
                    @test x isa T
                    @test value(x) ≈ 0.03
                end

                let x = T(0.1) / T(0.3)
                    @test x isa T
                    @test value(x) ≈ 0.1/0.3
                end

                let x = -T(0.1)
                    @test x isa T
                    @test value(x) ≈ -0.1
                end
            end

            @testset "math" begin
                for x in [-0.1, 0.1]
                    ax = abs(T(x))
                    @test ax isa T
                    @test value(ax) ≈ abs(x)
                end

                let x = zero(T(42))
                    @test x isa T
                    @test value(x) == 0
                end

                let x = one(T(42))
                    @test x isa T
                    @test value(x) == 1
                end
            end

            @testset "isless" begin
                @test T(1) < 2.
                @test 1    < T(2)
                @test T(1) < T(2)

                @test T(1) <= 2.
                @test 1    <= T(2)
                @test T(1) <= T(2)

                @test 2.   > T(1)
                @test T(2) > 1
                @test T(2) > T(1)

                @test 2.   >= T(1)
                @test T(2) >= 1
                @test T(2) >= T(1)
            end

            @testset "linsolve" begin
                let
                    tol(SFloat64) = 1e-15
                    tol(SFloat32) = 1e-6

                    A = Float64[1 2 3;
                                4 5 6;
                                7 8 10]
                    B = Float64[1 ; 2 ; 3]
                    Xref = A \ B

                    function check(Xsto)
                        value(Xsto[1]) ≈ Xref[1]     || return false
                        value(Xsto[2]) ≈ Xref[2]     || return false
                        abs(value(Xsto[3])) < tol(T) || return false
                        true
                    end
                    @test all(check(T.(A) \ T.(B)) for _ in 1:10)
                end
            end
        end
    end

    @testset "mixed" begin
        using StochasticArithmetic: det_type, sto_type
        @test det_type(SFloat64) == Float64
        @test det_type(SFloat32) == Float32
        @test det_type(Float32)  == Float32
        @test det_type(Int32)    == Int32

        @test sto_type(Float64)  == SFloat64
        @test sto_type(Float32)  == SFloat32

        @test promote_type(SFloat64, Float32) == SFloat64
        @test promote_type(SFloat32, Float32) == SFloat32
        @test promote_type(SFloat32, Int)     == SFloat32
        @test promote_type(SFloat32, Float64) == SFloat64

        let
            x32 = SFloat32(32)
            x64 = SFloat64(64)

            y32 = 32.f0
            y64 = 64.0

            # Sto + Sto (mixed prec)

            ## S32 + S64 -> S64
            x = x32 + x64
            @test x isa SFloat64
            @test value(x) ≈ 96
            x = x64 + x32
            @test x isa SFloat64
            @test value(x) ≈ 96

            # Sto + Flo (same prec)

            ## S64 + F64 -> S64
            x = x64 + y64
            @test x isa SFloat64
            @test value(x) ≈ 128
            x = y64 + x64
            @test x isa SFloat64
            @test value(x) ≈ 128

            ## S32 + F32 -> S32
            x = x32 + y32
            @test x isa SFloat32
            @test value(x) ≈ 64
            x = y32 + x32
            @test x isa SFloat32
            @test value(x) ≈ 64

            # Sto + Flo (mixed prec)

            ## S64 + F32 -> S64
            x = x64 + y32
            @test x isa SFloat64
            @test value(x) ≈ 96
            x = y32 + x64
            @test x isa SFloat64
            @test value(x) ≈ 96

            ## F64 + S32 -> S64
            x = x32 + y64
            @test x isa SFloat64
            @test value(x) ≈ 96
            x = y64 + x32
            @test x isa SFloat64
            @test value(x) ≈ 96

            # Comparisons
            @test x32 < x64
            @test x32 < y64
            @test y32 < x64
        end
    end

    @testset "reliable_digits" begin
        let
            mu, s = (@reliable_digits dot(SFloat64.(x1), SFloat64.(y1)))
            @test mu ≈ dot(x1, y1)
            @test s  ≈ 13.37        atol = 0.01
        end

        let
            mu, s = (@reliable_digits dot(SFloat32.(x1), SFloat32.(y1)))
            @test mu ≈ dot(Float32.(x1), Float32.(y1))
            @test s  ≈ 4.87         atol = 0.01
        end

        let
            mu, s = (@reliable_digits dot(SFloat64.(x2), SFloat64.(y2)))
            @test mu ≈ dot(x2, y2)  atol = 0.01
            @test s  ≈ 3.53         atol = 0.01
        end

        let
            mu, s = (@reliable_digits dot(SFloat32.(x2), SFloat32.(y2)))
            @test s  ≈ -0.25        atol = 0.01
        end
    end
end


text:
# Bruker I/O functions

"""     read_bruker_binary(fname)
Bruker FID, and processed frequency-domain data files (fid, 1r, 1i) contain a flat
binary representation of their data points, each of which is an Int32."""
function read_bruker_binary(fname)
    reinterpret(Int32, open(read, fname))
end

function parse_float_list(m)
    lines = split(m)
    # The first line of a match is "(0..n)"; we don't need this.
    [parse(Float64, s) for s in split(m)[2:end]]
end

function parse_or(::Type{T}, s, default::T) where T
    v = tryparse(T, s)
    v === nothing ? default : v
end

filters = [ ( Set(["SW", "O1", "SFO1", "SF", "BF1"]),
              s -> parse(Float64, s) ),
            ( Set(["TD", "NS", "DS", "SI"]),
              s -> parse(Int, s) ),
            ( Set(["D", "P"]),
              parse_float_list),
            ( Set(["PULPROG"]),
              x -> strip(x)[2:end-1] ),
            ( Set(["AUTOPOS"]),
              x -> parse_or(Int, strip(x)[2:end-1], 0) )
]

function read_params(file)
    contents = read(file, String)
    matches = eachmatch(r"##\$?(.*?)=\s?([^#]+)"s, contents)
    res = Dict(m.captures[1] => parse_param(m.captures[1], m.captures[2]) for m in matches)

    # a few little tweaks
    if "O1" in keys(res) && !("O1P" in keys(res))
        res["O1P"] = res["O1"] / res["SFO1"]
    end
    res
end

function read_intrng(file)
    lines = try
        readlines(file)
    catch
        return []
    end
    if length(lines) > 2 && strip(lines[1])[1] == 'A'
        [tuple(map(s -> parse(Float64, s), split(line)[1:2])...) for line in lines[3:end]]
    elseif length(lines) > 1 && strip(lines[1])[1] == 'P'
        [tuple(map(s -> parse(Float64, s), split(line))...) for line in lines[2:end]]
    else
        []
    end
end

function parse_param(param, val)
    for (names, fun) in filters
        if param in names
            return fun(val)
        end
    end
    return strip(val)
end

function ProcessedSpectrum(path :: AbstractString, procno :: Int)
    re_ft = float(read_bruker_binary(joinpath(path, "1r")))
    im_ft = float(read_bruker_binary(joinpath(path, "1i")))
    params = read_params(joinpath(path, "proc"))
    title = read(joinpath(path, "title"), String)
    intrng = read_intrng(joinpath(path, "intrng"))
    ProcessedSpectrum(re_ft, im_ft, params, intrng, procno, title)
end

ProcessedSpectrum(path::AbstractString) = ProcessedSpectrum(path, parse(Int, basename(path)))

Spectrum(path :: AbstractString, procnos :: AbstractArray{Int}, default_proc :: Int) = begin
    fid = float(read_bruker_binary(joinpath(path, "fid")))
    acqu = read_params(joinpath(path, "acqu"))
    name = basename(dirname(path))
    expno = parse(Int, basename(path))
    procs = Dict()
    for procno in procnos
        proc_path = joinpath(path, "pdata", string(procno))
        procs[procno] = ProcessedSpectrum(proc_path, procno)
    end
    Spectrum(fid, acqu, procs, default_proc, name, expno)
end

Spectrum(path :: AbstractString, procno :: Int) = Spectrum(path, [procno], procno)

Spectrum(path :: AbstractString, procnos :: AbstractArray{Int}) = Spectrum(path, procnos, minimum(procnos))

Spectrum(path :: AbstractString) = begin
    procnos = [parse(Int,n) for n in readdir(joinpath(path, "pdata"))]
    Spectrum(path, procnos)
end


text:
function downloadGISSMOedata(molecule_name::String,
    simulation_label,
    entry_label,
    coupling_info_folder::String,
    edata_save_folder::String,
    eData_URL::String;
    unique_cs_tol::T = 1e-6,
    zero_tol_sigdigits = 6,
    keep_zip_file = false) where T <: Real

    # Download zip file.
    entry_save_folder = joinpath(edata_save_folder, "$(entry_label)_$(simulation_label)")
    mkpath(entry_save_folder)

    file_name = "nmredata.zip"
    dest = joinpath(entry_save_folder, file_name)
    t = @task begin; isfile(dest) || download(eData_URL, dest); end
    schedule(t); wait(t)

    #q = filter(x->endswith(x, ".zip"), readdir(entry_save_folder))

    # unzip.
    InfoZIP.unzip(dest, entry_save_folder)

    # remove zip file.
    if !keep_zip_file
        rm(dest)
    end

    # get SDF filename.
    tmp = filter(isdir, readdir(entry_save_folder; join=true))
    @assert length(tmp) == 1 # should be only one folder, which is the unzipped one.
    folder_path = tmp[1]

    tmp = filter(x->endswith(x, ".sdf"), readdir(folder_path))
    @assert length(tmp) == 1 # should be exactly 1 SDF file.
    SDF_name = tmp[1]
    SDF_path = joinpath(folder_path, SDF_name)

    file_strings = readlines(SDF_path)

    #### find the section on chemical shifts.
    cs_string = "> <NMREDATA_ASSIGNMENT>"

    H_IDs, H_css = parsecssfromNMReDATA(file_strings)

    ##### find the section for J-coupling. J_string = "> <NMREDATA_J>"
    J_IDs, J_vals = GISSMOReader.parseJsfromNMReDATA(file_strings)

    # store.
    mkpath(coupling_info_folder)
    save_path = joinpath(coupling_info_folder, "$(entry_label)_$(simulation_label).json")

    cs_dict = cstodict(H_IDs, H_css)
    J_dict = Jcouplingtodict(J_IDs, J_vals)
    # dict1 = Dict("H_IDs" => H_IDs, "H_css" => H_css,
    #             "J_IDs" => J_IDs, "J_vals" => J_vals)
    dict1 = Dict("chemical shift" => cs_dict, "J-coupling" => J_dict)
    stringdata = JSON.json(dict1)

    open(save_path, "w") do f
        JSON3.pretty(f, stringdata)
        println(f)
    end

    return save_path
end

function Jcouplingtodict(J_IDs::Vector{Tuple{Int,Int}}, J_vals::Vector{T}) where T

    return collect( Dict("ID1" => J_IDs[i][1], "ID2" => J_IDs[i][2], "value" => J_vals[i]) for i = 1:length(J_vals) )
end

function cstodict(H_IDs::Vector{Int}, H_css::Vector{T}) where T

    return collect( Dict("ID" => H_IDs[i], "value" => H_css[i]) for i = 1:length(H_css) )
end

# come back to this later.
function tempfunc()
    #J_string = "> <NMREDATA_J>"
    J_IDs, J_vals = parseJsfromNMReDATA(file_strings)

    ### attempt on auto-grouping.

    J_IDs_sys, J_vals_sys, H_IDs_sys = partitionJcouplings(J_IDs, J_vals)

    css_sys = getcssys(H_IDs_sys, H_IDs, H_css)

    𝐽_IDs_sys = getJsys(J_IDs_sys, J_vals_sys, H_IDs_sys)

    H_singlets, cs_singlets,
        cs_singlets_compact = createsingletsysems(H_IDs, H_IDs_sys,
                                    H_css;
                                    zero_tol = unique_cs_tol)
    #
    ### remove spin groups that only have one unique chem shift, but have J-coupling between the atoms.
    # Each of such groups are are really one singlet group.
    removeisolatedcs!(css_sys,
                    cs_singlets,
                    cs_singlets_compact,
                    H_singlets,
                    H_IDs_sys,
                    J_IDs_sys,
                    J_vals_sys,
                    𝐽_IDs_sys;
                    zero_tol_sigdigits = zero_tol_sigdigits)
    #
    println("Entry and simulation: ", "$(entry_label)_$(simulation_label)")
    println("cs_singlets: ", cs_singlets)
    println("cs_singlets_compact: ", cs_singlets_compact)
    println("H_singlets: ", H_singlets)

    println("css_sys: ", css_sys)
    println("J_vals_sys: ", J_vals_sys)
    println("𝐽_IDs_sys: ", 𝐽_IDs_sys)
    println()
    #### end code for moving spin group to singlet group.

    #
    cs_LUT, p_cs_sys = constructLUTcss(css_sys; zero_tol = unique_cs_tol)

    cs_len_sys = getcslengthfromLUT(cs_LUT)

    # ### cs and J-values upper and lower bounds.
    # perturblbfunc = xx->additiveperturbation(xx, δ_lb)
    # perturbubfunc = xx->additiveperturbation(xx, δ_ub)
    #
    # #
    # J_lb_sys = getJperturbed(J_vals_sys, perturblbfunc)
    # J_ub_sys = getJperturbed(J_vals_sys, perturbubfunc)
    #
    # cs_lb_sys = getJperturbed(css_sys, perturblbfunc)
    # cs_ub_sys = getJperturbed(css_sys, perturbubfunc)

    # store.
    save_path = joinpath(JLD_save_folder, "$(entry_label)_$(simulation_label).jld")

    JLD.save(save_path,
        "css_sys", css_sys,
        "J_IDs_sys", 𝐽_IDs_sys,
        "J_vals_sys", J_vals_sys,
        "cs_singlets", cs_singlets,
        "cs_LUT", cs_LUT)
end


text:
struct DictColumnTable <: AbstractColumns
    schema::Schema
    values::OrderedDict{Symbol, AbstractVector}
end

"""
    Tables.dictcolumntable(x) => Tables.DictColumnTable

Take any Tables.jl-compatible source `x` and return a `DictColumnTable`, which
can be thought of as a `OrderedDict` mapping column names as `Symbol`s to `AbstractVector`s.
The order of the input table columns is preserved via the `Tables.schema(::DictColumnTable)`.

For "schema-less" input tables, `dictcolumntable` employs a "column unioning" behavior,
as opposed to inferring the schema from the first row like `Tables.columns`. This
means that as rows are iterated, each value from the row is joined into an aggregate
final set of columns. This is especially useful when input table rows may not include
columns if the value is missing, instead of including an actual value `missing`, which
is common in json, for example. This results in a performance cost tracking all seen
values and inferring the final unioned schemas, so it's recommended to use only when
needed.
"""
function dictcolumntable(x)
    if columnaccess(x)
        cols = columns(x)
        names = columnnames(cols)
        sch = schema(cols)
        out = OrderedDict(nm => getcolumn(cols, nm) for nm in names)
    else
        r = rows(x)
        L = Base.IteratorSize(typeof(r))
        len = Base.haslength(r) ? length(r) : 0
        sch = schema(r)
        if sch !== nothing
            names, types = sch.names, sch.types
            out = OrderedDict{Int, AbstractVector}(i => allocatecolumn(types[i], len) for i = 1:length(types))
            for (i, row) in enumerate(r)
                eachcolumns(add!, sch, row, out, L, i)
            end
            out = OrderedDict(names[k] => v for (k, v) in out)
        else
            names = Symbol[]
            seen = Set{Symbol}()
            out = OrderedDict{Symbol, AbstractVector}()
            for (i, row) in enumerate(r)
                for nm in columnnames(row)
                    push!(seen, nm)
                    val = getcolumn(row, nm)
                    if haskey(out, nm)
                        col = out[nm]
                        if typeof(val) <: eltype(col)
                            add!(val, 0, nm, col, L, i)
                        else # widen column type
                            new = allocatecolumn(promote_type(eltype(col), typeof(val)), length(col))
                            i > 1 && copyto!(new, 1, col, 1, i - 1)
                            add!(new, val, L, i)
                            out[nm] = new
                        end
                    else
                        push!(names, nm)
                        if i == 1
                            new = allocatecolumn(typeof(val), len)
                            add!(new, val, L, i)
                            out[nm] = new
                        else
                            new = allocatecolumn(Union{Missing, typeof(val)}, len)
                            add!(new, val, L, i)
                            out[nm] = new
                        end
                    end
                end
                for nm in names
                    if !(nm in seen)
                        col = out[nm]
                        if !(eltype(col) >: Missing)
                            new = allocatecolumn(Union{Missing, eltype(col)}, len)
                            i > 1 && copyto!(new, 1, col, 1, i - 1)
                            out[nm] = new
                        end
                    end
                end
                empty!(seen)
            end
            sch = Schema(collect(keys(out)), eltype.(values(out)))
        end
    end
    return DictColumnTable(sch, out)
end

istable(::Type{DictColumnTable}) = true
columnaccess(::Type{DictColumnTable}) = true
columns(x::DictColumnTable) = x
schema(x::DictColumnTable) = getfield(x, :schema)
columnnames(x::DictColumnTable) = getfield(x, :schema).names
getcolumn(x::DictColumnTable, i::Int) = getfield(x, :values)[columnnames(x)[i]]
getcolumn(x::DictColumnTable, nm::Symbol) = getfield(x, :values)[nm]

struct DictRowTable
    names::Vector{Symbol}
    types::Dict{Symbol, Type}
    values::Vector{Dict{Symbol, Any}}
end

isrowtable(::Type{DictRowTable}) = true
schema(x::DictRowTable) = Schema(getfield(x, :names), [getfield(x, :types)[nm] for nm in getfield(x, :names)])

struct DictRow <: AbstractRow
    names::Vector{Symbol}
    row::Dict{Symbol, Any}
end

columnnames(x::DictRow) = getfield(x, :names)
getcolumn(x::DictRow, i::Int) = get(getfield(x, :row), columnnames(x)[i], missing)
getcolumn(x::DictRow, nm::Symbol) = get(getfield(x, :row), nm, missing)

Base.IteratorSize(::Type{DictRowTable}) = Base.HasLength()
Base.length(x::DictRowTable) = length(getfield(x, :values))
Base.IteratorEltype(::Type{DictRowTable}) = Base.HasEltype()
Base.eltype(x::DictRowTable) = DictRow

function Base.iterate(x::DictRowTable, st=1)
    st > length(x) && return nothing
    return DictRow(x.names, x.values[st]), st + 1
end

"""
    Tables.dictrowtable(x) => Tables.DictRowTable

Take any Tables.jl-compatible source `x` and return a `DictRowTable`, which
can be thought of as a `Vector` of `OrderedDict` rows mapping column names as `Symbol`s to values.
The order of the input table columns is preserved via the `Tables.schema(::DictRowTable)`.

For "schema-less" input tables, `dictrowtable` employs a "column unioning" behavior,
as opposed to inferring the schema from the first row like `Tables.columns`. This
means that as rows are iterated, each value from the row is joined into an aggregate
final set of columns. This is especially useful when input table rows may not include
columns if the value is missing, instead of including an actual value `missing`, which
is common in json, for example. This results in a performance cost tracking all seen
values and inferring the final unioned schemas, so it's recommended to use only when
the union behavior is needed.
"""
function dictrowtable(x)
    names = Symbol[]
    seen = Set{Symbol}()
    types = OrderedDict{Symbol, Type}()
    r = rows(x)
    L = Base.IteratorSize(typeof(r))
    out = Vector{OrderedDict{Symbol, Any}}(undef, Base.haslength(r) ? length(r) : 0)
    for (i, drow) in enumerate(r)
        row = OrderedDict{Symbol, Any}(nm => getcolumn(drow, nm) for nm in columnnames(drow))
        add!(row, 0, :_, out, L, i)
        if isempty(names)
            for (k, v) in row
                push!(names, k)
                types[k] = typeof(v)
            end
            seen = Set(names)
        else
            for nm in names
                if haskey(row, nm)
                    T = types[nm]
                    v = row[nm]
                    if !(typeof(v) <: T)
                        types[nm] = Union{T, typeof(v)}
                    end
                else
                    types[nm] = Union{Missing, types[nm]}
                end
            end
            for (k, v) in row
                if !(k in seen)
                    push!(seen, k)
                    push!(names, k)
                    # we mark the type as Union{T, Missing} here because
                    # we're at least on the 2nd row, and we didn't see
                    # this column in the 1st row, so its value will be
                    # `missing` for that row
                    types[k] = Union{typeof(v), Missing}
                end
            end
        end
    end
    return DictRowTable(names, types, out)
end


text:
function count_insertions(seed::Dict, L::Int64)

nseq = length(seed)
ins = zeros((nseq,L))
pointseq = 0

for name in keys(seed)
	s = seed[name]
	pointseq += 1
	pos = L
	for k in 0:length(s)-1
		if isuppercase(s[length(s)-k]) || s[length(s)-k] == '-'
			go = true
			j = 1
			delta = 0
			while go
				if length(s)-k-j <= 0
					go = false
				elseif  isuppercase(s[length(s)-k-j]) || s[length(s)-k-j] == '-'
					go = false
				elseif s[length(s)-k-j] == '.'
					j += 1
				elseif islowercase(s[length(s)-k-j])
					delta += 1
					j += 1
				end
			end
			ins[pointseq,pos] = delta
			pos -= 1
		end
	end
end

return ins

end

function maximum_likelihood(ins, L, nseq)

maxit = 1e8
tol = 1e-6
eta = 1e-3
lo0 = 0.0
le0 = 0.0
lambda_o = zeros(L)
lambda_e = zeros(L)

for i in 1:L
	allz = false
	f0 = sum(x -> x == 0, ins[:,i]) / nseq
	m = mean(ins[:,i])
	if f0 == 1.0 && le0 != 0.0 && lo0 != 0.0
		lambda_e[i] = le0
		lambda_o[i] = lo0
		continue
	end
	if f0 == 1.0
		allz = true
		f0 = 1.0 - 1e-3
	end
	lo = 1.0
	le = 1.0
	eps = 1.0
	it = 1
	while eps > tol && it < maxit
		dLdle = ( exp(-lo - le) *((1.0 - exp(-le))^-2) )/(1.0 + exp(-lo) * ( (1.0-exp(-le))^-1 )) - m + (1.0-f0) - (2.0/nseq)*le
		le = le + eta * dLdle
		dLdlo = ( exp(-lo) * ((1.0-exp(-le))^-1) ) / (1.0 + exp(-lo)* ( (1.0 -exp(-le))^-1 )) - (1.0 - f0) - (2.0/nseq) *lo
		lo = lo + eta * dLdlo
		eps = max(abs(dLdle), abs(dLdlo))
		it += 1
	end
	lambda_o[i] = lo
	lambda_e[i] = le
	if allz && le0 == 0.0 && lo0 == 0.0
		le0 = le
		lo0 = lo
	end
	@printf("site: %i it: %2.1e eps: %.2e f0: %.2f m: %.2f ins pen: ( %.3f, %.3f )\n", i, it, eps, f0, m,  lo, le)
end

return lambda_o, lambda_e

end

function infer_ins_pen(seed::Dict, L::Int64)

	ins = count_insertions(seed, L)
	lo, le = maximum_likelihood(ins, L, length(seed))

	return lo, le
end


text:
function test_convergence(bquad, sampler; D=3)
    rng = Random.MersenneTwister(42)
    prior = MvNormal(0.5 * ones(D))
    likelihood = MvNormal(7.0 * ones(D))
    logintegrand = x -> logpdf(likelihood, x)
    m = BayesModel(prior, logintegrand)
    pZ, x = bquad(m, sampler; nsamples=100)
    trueZ = exp(-0.5 * (logdet(cov(prior) + cov(likelihood)) + D * log(2π)))
    @test mean(pZ) ≈ trueZ rtol = 1e-2
end


text:
module TestSatisfiesConstraints

import Wordle
import Test: @test, @test_throws

constraints = Wordle.Constraints()
trial = Wordle.Trial(Wordle.Candidate("alone"), Wordle.Result("⬜⬜⬜⬜⬜"))
Wordle.update_constraints!(constraints, trial)

@test !Wordle.satisfies_constraints(Wordle.Candidate("alone"), constraints)
@test !Wordle.satisfies_constraints(Wordle.Candidate("alive"), constraints)
@test Wordle.satisfies_constraints(Wordle.Candidate("trust"), constraints)
@test Wordle.satisfies_constraints(Wordle.Candidate("birth"), constraints)

# TODO: Add lots more cases here

end # module


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

JLLWrappers.@generate_main_file_header("Libksba")
JLLWrappers.@generate_main_file("Libksba", UUID("52a58d30-3731-5a3f-8361-26ffb4f63669"))
end  # module Libksba_jll


text:
using Plots, StatsBase, Statistics, LaTeXStrings, JLD

function ColoringTheNode(i, j, Network, Num)
    check = false
    if Network[i, j] == 0
        Network[i, j] = Num
        check = true
        Num += 1
    end
    return Network, Num, check
end

function FindNeighbors(i, j, Network, dim)
    Numbereds = 0
    colorList = []
    for neighbor in ([0,1], [1,0], [0,-1], [-1,0])
        if j + neighbor[2] >= 1 && j + neighbor[2] <= dim && Network[([i, j] + neighbor)...] > 0
            push!(colorList,Network[([i, j] + neighbor)...])
            Numbereds += 1
        end
    end
    push!(colorList, Network[i, j])
    return colorList, Numbereds #colorList contain the node and its neighbors
end

function JoinColors(Network, dim, colorList)
    mincolor = min(colorList...)
    for i in 1:dim
        for j in 1:dim
            if Network[i, j] in colorList
                Network[i, j] = mincolor
            end
        end
    end
    return Network
end

function InitialNetwork(dim, P)
    Network = sample([-1000,0], Weights([1-P, P]),(dim,dim))
    return Network
end

function PercolationCheck(dim,P)
    Num = 2
    maxnum = 10000
    Network = InitialNetwork(dim, P)
    Network[1,:] = fill(1,dim)
    Network[end,:] = fill(maxnum,dim)
    TotalNets = []
    for i in 2:dim-1
        for j in 1:dim
            Network, Num, check = ColoringTheNode(i, j, Network, Num)
            if check == true
                colorList, Numbereds = FindNeighbors(i, j, Network, dim)
                if Numbereds == 0
                    continue
                elseif Numbereds == 1
                    Network[i, j] = colorList[1]
                else
                    Network = JoinColors(Network, dim, colorList)
                end
            end
            push!(TotalNets,copy(Network))
        end
    end
    return TotalNets
end

Plist = [0.5,0.55,0.6,0.65,0.7]
for P in Plist
    Data = PercolationCheck(100,P)
    save("Percolation$P.jld", "data", Data)
    anim = Animation()
    for i in 110:100:length(Data)
        plt = heatmap(Data[i], c = cgrad(:copper, 20), legend = false, border = :none, title = "Percolation for P = $P")
        frame(anim, plt)
    end
    gif(anim, "../../Figs/Q4/Percolation$P.gif")
end


text:
using GreenFunc
using Test, StaticArrays, LinearAlgebra, Printf, Statistics, Lehmann, CompositeGrids
using JLD2, FileIO
using CodecZlib

if isempty(ARGS)
    include("test_Green.jl")
    #include("interpolate.jl")
else
    include(ARGS[1])
end


text:
module SymbolicRegression

# Types
export Population,
    PopMember,
    HallOfFame,
    Options,

    #Functions:
    RunSR,
    SRCycle,
    calculateParetoFrontier,
    countNodes,
    copyNode,
    printTree,
    stringTree,
    evalTreeArray,
    node_to_symbolic,
    symbolic_to_node,
    custom_simplify,
    simplifyWithSymbolicUtils,
    combineOperators,

    #Operators:
    plus,
    sub,
    mult,
    square,
    cube,
    pow,
    div,
    logm,
    logm2,
    logm10,
    sqrtm,
    neg,
    greater,
    relu,
    logical_or,
    logical_and

using Printf: @printf
using Distributed

include("ProgramConstants.jl")
include("Operators.jl")
include("Options.jl")
include("Dataset.jl")
include("Equation.jl")
include("LossFunctions.jl")
include("Utils.jl")
include("EvaluateEquation.jl")
include("MutationFunctions.jl")
include("InterfaceSymbolicUtils.jl")
include("CustomSymbolicUtilsSimplification.jl")
include("SimplifyEquation.jl")
include("PopMember.jl")
include("HallOfFame.jl")
include("CheckConstraints.jl")
include("Mutate.jl")
include("Population.jl")
include("RegularizedEvolution.jl")
include("SingleIteration.jl")
include("ConstantOptimization.jl")

function RunSR(X::AbstractMatrix{T}, y::AbstractVector{T};
        niterations::Int=10,
        weights::Union{AbstractVector{T}, Nothing}=nothing,
        varMap::Union{Array{String, 1}, Nothing}=nothing,
        options::Options=Options()) where {T<:Real}

    testConfiguration(options)

    if size(X)[1] > 10000
        if !options.batching
            println("Note: you are running with more than 10,000 datapoints. You should consider turning on batching (`options.batching`), and also if you need that many datapoints. Unless you have a large amount of noise (in which case you should smooth your dataset first), generally < 10,000 datapoints is enough to find a functional form.")
        end
    end

    dataset = Dataset(X, y,
                     weights=weights,
                     varMap=varMap)

    if dataset.weighted
        avgy = sum(dataset.y .* dataset.weights)/sum(dataset.weights)
        baselineMSE = MSE(dataset.y, ones(T, dataset.n) .* avgy, dataset.weights)
    else
        avgy = sum(dataset.y)/dataset.n
        baselineMSE = MSE(dataset.y, ones(T, dataset.n) .* avgy)
    end

    #TODO - remove this as redundant
    nfeatures = dataset.nfeatures

    # 1. Start a population on every process
    allPops = Future[]
    # Set up a channel to send finished populations back to head node
    channels = [RemoteChannel(1) for j=1:options.npopulations]
    bestSubPops = [Population(dataset, baselineMSE, npop=1, options=options, nfeatures=nfeatures) for j=1:options.npopulations]
    hallOfFame = HallOfFame(options)
    actualMaxsize = options.maxsize + maxdegree
    frequencyComplexity = ones(T, actualMaxsize)
    curmaxsize = 3
    if options.warmupMaxsize == 0
        curmaxsize = options.maxsize
    end

    for i=1:options.npopulations
        future = @spawnat :any Population(dataset, baselineMSE, npop=options.npop, nlength=3, options=options, nfeatures=nfeatures)
        push!(allPops, future)
    end

    # # 2. Start the cycle on every process:
    @sync for i=1:options.npopulations
        @async allPops[i] = @spawnat :any SRCycle(dataset, baselineMSE, fetch(allPops[i]), options.ncyclesperiteration, curmaxsize, copy(frequencyComplexity)/sum(frequencyComplexity), verbosity=options.verbosity, options=options)
    end
    println("Started!")
    cycles_complete = options.npopulations * niterations
    if options.warmupMaxsize != 0
        curmaxsize += 1
        if curmaxsize > options.maxsize
            curmaxsize = options.maxsize
        end
    end

    last_print_time = time()
    num_equations = 0.0
    print_every_n_seconds = 5
    equation_speed = Float32[]

    for i=1:options.npopulations
        # Start listening for each population to finish:
        @async put!(channels[i], fetch(allPops[i]))
    end

    while cycles_complete > 0
        @inbounds for i=1:options.npopulations
            # Non-blocking check if a population is ready:
            if isready(channels[i])
                # Take the fetch operation from the channel since its ready
                cur_pop = take!(channels[i])
                bestSubPops[i] = bestSubPop(cur_pop, topn=options.topn)
                # bestSubPops[i] = bestSubPopParetoDominating(cur_pop, topn=options.topn)

                #Try normal copy...
                bestPops = Population([member for pop in bestSubPops for member in pop.members])

                for member in cur_pop.members
                    size = countNodes(member.tree)
                    frequencyComplexity[size] += 1
                    # println(member, hallOfFame.members[size])
                    actualMaxsize = options.maxsize + maxdegree
                    if size < actualMaxsize && member.score < hallOfFame.members[size].score
                        hallOfFame.members[size] = copyPopMember(member)
                        hallOfFame.exists[size] = true
                    end
                end

                # Dominating pareto curve - must be better than all simpler equations
                dominating = calculateParetoFrontier(dataset, hallOfFame, options)
                open(options.hofFile, "w") do io
                    println(io,"Complexity|MSE|Equation")
                    for member in dominating
                        println(io, "$(countNodes(member.tree))|$(member.score)|$(stringTree(member.tree, options, varMap=dataset.varMap))")
                    end
                end
                cp(options.hofFile, options.hofFile*".bkup", force=true)

                # Try normal copy otherwise.
                if options.migration
                    for k in rand(1:options.npop, round(Int, options.npop*options.fractionReplaced))
                        to_copy = rand(1:size(bestPops.members)[1])
                        cur_pop.members[k] = PopMember(
                            copyNode(bestPops.members[to_copy].tree),
                            bestPops.members[to_copy].score)
                    end
                end

                if options.hofMigration && size(dominating)[1] > 0
                    for k in rand(1:options.npop, round(Int, options.npop*options.fractionReplacedHof))
                        # Copy in case one gets used twice
                        to_copy = rand(1:size(dominating)[1])
                        cur_pop.members[k] = PopMember(
                           copyNode(dominating[to_copy].tree), dominating[to_copy].score
                        )
                    end
                end

                # TODO: Turn off this async when debugging - any errors in this code
                #         are silent.
                # begin
                @async begin
                    allPops[i] = @spawnat :any let
                        tmp_pop = SRCycle(dataset, baselineMSE, cur_pop, options.ncyclesperiteration, curmaxsize, copy(frequencyComplexity)/sum(frequencyComplexity), verbosity=options.verbosity, options=options)
                        @inbounds @simd for j=1:tmp_pop.n
                            tmp_pop.members[j].tree = simplifyTree(tmp_pop.members[j].tree, options)
                            tmp_pop.members[j].tree = combineOperators(tmp_pop.members[j].tree, options)
                            tmp_pop.members[j].tree = simplifyWithSymbolicUtils(tmp_pop.members[j].tree, options)
                            if rand() < 0.1 && options.shouldOptimizeConstants
                                tmp_pop.members[j] = optimizeConstants(dataset, baselineMSE, tmp_pop.members[j], options)
                            end
                        end
                        tmp_pop = finalizeScores(dataset, baselineMSE, tmp_pop, options)
                        tmp_pop
                    end
                    put!(channels[i], fetch(allPops[i]))
                end

                cycles_complete -= 1
                cycles_elapsed = options.npopulations * niterations - cycles_complete
                if options.warmupMaxsize != 0 && cycles_elapsed % options.warmupMaxsize == 0
                    curmaxsize += 1
                    if curmaxsize > options.maxsize
                        curmaxsize = options.maxsize
                    end
                end
                num_equations += options.ncyclesperiteration * options.npop / 10.0
            end
        end
        sleep(1e-3)
        elapsed = time() - last_print_time
        #Update if time has passed, and some new equations generated.
        if elapsed > print_every_n_seconds && num_equations > 0.0
            # Dominating pareto curve - must be better than all simpler equations
            current_speed = num_equations/elapsed
            average_over_m_measurements = 10 #for print_every...=5, this gives 50 second running average
            push!(equation_speed, current_speed)
            if length(equation_speed) > average_over_m_measurements
                deleteat!(equation_speed, 1)
            end
            average_speed = sum(equation_speed)/length(equation_speed)
            curMSE = baselineMSE
            lastMSE = curMSE
            lastComplexity = 0
            if options.verbosity > 0
                @printf("\n")
                @printf("Cycles per second: %.3e\n", round(average_speed, sigdigits=3))
                cycles_elapsed = options.npopulations * niterations - cycles_complete
                @printf("Progress: %d / %d total iterations (%.3f%%)\n", cycles_elapsed, options.npopulations * niterations, 100.0*cycles_elapsed/(options.npopulations*niterations))
                @printf("Hall of Fame:\n")
                @printf("-----------------------------------------\n")
                @printf("%-10s  %-8s   %-8s  %-8s\n", "Complexity", "MSE", "Score", "Equation")
                @printf("%-10d  %-8.3e  %-8.3e  %-.f\n", 0, curMSE, 0f0, avgy)
            end

            #TODO - call pareto function!
            actualMaxsize = options.maxsize + maxdegree
            for size=1:actualMaxsize
                if hallOfFame.exists[size]
                    member = hallOfFame.members[size]
                    curMSE = EvalLoss(member.tree, dataset, options)
                    numberSmallerAndBetter = 0
                    for i=1:(size-1)
                        hofMSE = EvalLoss(hallOfFame.members[i].tree, dataset, options)
                        if (hallOfFame.exists[size] && curMSE > hofMSE)
                            numberSmallerAndBetter += 1
                            break
                        end
                    end
                    betterThanAllSmaller = (numberSmallerAndBetter == 0)
                    if betterThanAllSmaller
                        delta_c = size - lastComplexity
                        delta_l_mse = log(curMSE/lastMSE)
                        score = convert(Float32, -delta_l_mse/delta_c)
                        if options.verbosity > 0
                            @printf("%-10d  %-8.3e  %-8.3e  %-s\n" , size, curMSE, score, stringTree(member.tree, options, varMap=dataset.varMap))
                        end
                        lastMSE = curMSE
                        lastComplexity = size
                    end
                end
            end
            debug(options.verbosity, "")
            last_print_time = time()
            num_equations = 0.0
        end
    end
    return hallOfFame
end

end #module SR


text:
# ------- R 2 Julia Conversion ---------------
#function saveModels(d::Dict,fname::String)
#    j = JSON.json([ d])
#    open(fname,"w") do f
#        write(f,j)
#    end
#end

function expandM(di::Dict)
    d2=deepcopy(di)
    if d2[:modelName]==:occ
        d2[:dist] = Poisson()
        d2[:lnk] = LogLink()
        d2[:Buyer_Pos_P1_is1] = true
    end
    if d2[:modelName]==:dolocc
        d2[:dist] = Gamma()
        d2[:lnk] = LogLink()
        d2[:Buyer_Pos_P1_is1] = true
    end
    if d2[:modelName]==:pen
        d2[:dist] = Bernoulli()
        d2[:lnk] = LogitLink()
        d2[:Buyer_Pos_P1_is1] = false
    end
    return d2
end

function readModels(fname::String)
    if isfile(fname)
        d=Dict()
        j = JSON.parsefile(fname)[1]
        #for (key, value) in dict
        for (k, v) in j
            d2=Dict()
            if Symbol(k)==:mocc
                d2[:dist] = Poisson()
                d2[:lnk] = LogLink()
                d2[:Buyer_Pos_P1_is1] = true
            end
            if Symbol(k)==:mdolocc
                d2[:dist] = Gamma()
                d2[:lnk] = LogLink()
                d2[:Buyer_Pos_P1_is1] = true
            end
            if Symbol(k)==:mpen
                d2[:dist] = Bernoulli()
                d2[:lnk] = LogitLink()
                d2[:Buyer_Pos_P1_is1] = false
            end

            if typeof(v) == Array{Any,1}
                d[Symbol(k)] = [Symbol(vx) for vx in v]
            else
                for (k1, v1) in v
                    if typeof(v1) in [ Array{Any,1} ]
                        v2=[Symbol(v1x) for v1x in v1]
                    else
                        v2=Symbol(v1)
                    end
                    d2[Symbol(k1)] = v2
                end
                d[Symbol(k)] = d2
            end
        end
        for (k,v) in d
            if k in [:idolocc, :ipen, :iocc]
                d[k] = expandM(v)
            end
        end
        return d
    else
        return "file not found!"
    end
end

read_dfx(root::String) = readtable(root*"/dfx.csv")
save_dfx(root::String,dfx::DataFrame) = writetable(root*"/dfx.csv", dfx)
read_modelsDict(root::String) = readModels(root*"/modelsDict.json")
function save_modelsDict(root::String, d::Dict)
    j = JSON.json([d])
    open(root*"/modelsDict.json","w") do f
        write(f,j)
    end
end
read_dfd(root::String) = readtable(root*"/dfd.csv",header=true);
save_dfd(root::String,dfd::DataFrame) = writetable(root*"/dfd.csv", dfd)

read_dfr(root::String) = readtable(root*"/scored.csv",header=true);
save_dfr(root::String,dfr::DataFrame) = writetable(root*"/scored.csv", dfr)

# ------- R 2 Julia Conversion - END ---------------

#include("trd.jl")
#function Base.lowercase(df::DataFrame)
function lowercase!(df::DataFrame)
    for (i,v) in enumerate(names(df))   if v != Symbol(lowercase(string(v)))  rename!(df,v,Symbol(lowercase(string(v)))) end end
end
function Base.lowercase(arr::Array{Symbol})    return  Symbol[Symbol(lowercase(string(s))) for s in arr] end
function Base.lowercase(arr::Array{String})    return  String[lowercase(string(s)) for s in arr] end
function Base.lowercase(arr::Array{String})    return  String[lowercase(string(s)) for s in arr] end
function Base.lowercase(s::Symbol) return Symbol(lowercase(string(s))) end
function Base.lowercase(d::OrderedDict)
    vout=OrderedDict()
    for (k, v) in d  vout[lowercase(k)] = typeof(v) in [Bool,Float64,Int64] ? v : lowercase(v) end
    return vout
end

function grep(v::String,vars::Array{Symbol})   #function grep(v::AbstractString) filter(x->ismatch(Regex("$(v)"), string(x)), vars) end
    filter(x->ismatch(Regex("$(v)"), string(x)), vars)
end
function grep1(v::String,vars::Array{Symbol})   #function grep(v::AbstractString) filter(x->ismatch(Regex("$(v)"), string(x)), vars) end
    lst = filter(x->ismatch(Regex("$(v)"), string(x)), vars)
    return length(lst)>0 ? lst[1] : nothing
end

function grep(varr::Array{String},vars::Array{Symbol})
    arr_out=Symbol[]
    for v in varr
        res = filter(x->ismatch(Regex("$(v)"), string(x)), vars)
        arr_out = vcat(arr_out, res)
    end
    return arr_out
end
#function grep(varr::Array{String},vars::Array{Symbol}) return grep(convert(Array{String},varr),vars) end  #grep(["M","pan"],names(df_in))

function df2dict(idf::DataFrame)
    dout=OrderedDict()
    for n in names(idf)
       dout[n] = idf[1,Symbol(n)]
    end
    return dout
end

function df2dict(row::DataFrameRow)
    dout=OrderedDict()
    for n in names(row)
       dout[n] = row[Symbol(n)]
    end
    return dout
end

"""
app.cfg
P2_Competitor = true
pvalue_lvl = 0.20
excludedKeys = key1,key2,key3
exposed_flag_var = new_exposed_flag
sigLevel = 0.2
random_demos =
random_campaigns =
dropvars = goop,snop
scoring_vars =
TotalModelsOnly=false
"""

function loadCFG(d::OrderedDict,fn::String)
    function cvt(v::String, t::DataType)
        if t in [Bool]
            return lowercase(v) =="true" ? true : lowercase(v) =="false" ? false : NA
        elseif t in [nothing]
            return NA
        elseif t in [Array{String,1}]
            #return AbstractString[string(x) split(v,",")]
            return String[string(x) for x in split(replace(v,"\"",""),",")]
        elseif t in [Array{Symbol,1}]
            return Symbol[Symbol(x) for x in split(replace(v,"\"",""),",")]
        elseif t in [Symbol]
            return Symbol(v)
        elseif t in [String]
            return string(v)
        elseif t in [Float64]
            return parse(Float64,v)
        else
            return "__ooops__ TYPE not found!"
        end
    end
    if isfile(fn)
        f = open(fn)
        for ln in eachline(f)
            ln=replace(ln,"\n","")
            sidx=search(ln,"=").start-1
            eidx=search(ln,"=").stop+1
            if sidx>0
                k=strip(ln[1:sidx])
                v=strip(ln[eidx:end])
                if (length(k)>0) & (length(v)>0)
                    println("LINE:: ~"*k*"~"*v*"~")
                    #println("DEF::       ______ ",d[Symbol(k)],"  ----> ",typeof(d[Symbol(k)]) )
                    newk=Symbol(k)
                    dv=get(d,newk,nothing)
                    if dv!=nothing
                        newv = cvt(v,typeof(dv))
                        #println("IN::        ______ ",Symbol(k),"  ----> ",ct, " ....... ", typeof(ct))
                        #println( d[Symbol(k)],"|",v," {",typeof(d[Symbol(k)]),"}  ----->  ",newk,"|",newv,"  {",typeof(newv),"}")
                        d[newk] = newv
                    else
                        println(newk," in app.cfg, but not defaults")
                        #d[newk] = v
                    end
                else
                    #println("CFG ERROR :: invalid length : "*"~"*k*"~"*v*"~",typeof(d[Symbol(k)]))
                end
            else
                if length(ln)>0   println("CFG ERROR :: invalid cfg parameter : ~",ln,"~") end
            end
        end
        close(f)
    else
        println("Using Default configuration")
    end
    return d
end

# Retype DataFrame Columns
# Ensure/Convert DataFrame Column Types
function retypeBools!(dfd::DataFrame,vars::Array{Symbol}=Symbol[])
    if length(vars)==0 vars=names(dfd) end
    for c in vars
        if typeof(dfd[c]) in [DataArray{Any,1}, DataArray{String,1}]  #DataArray{Bool,1}
            println(c," ~ ",typeof(dfd[c]))
            if  typeof(dfd[c]) == DataArray{Any,1}
                if unique(map(x-> typeof(x),dfd[c])) == [Bool]
                    println("Converting Col")
                    dfd[c] = convert(Array{Bool}, dfd[c])
                else
                    println("Type = Any; but not Bool")
                end
            else
                if true in unique(isna(dfd[c]))
                    println("String Col has NA's")
                else
                    ul=convert(Array{String},unique(map(x -> lowercase(x), dfd[c] )))
                    if (length(setdiff(ul,["true","false"]))==0) & (length(intersect(ul,["true","false"]))>0)
                        println("Mapping String Col")
                        dfd[c]=  map(x-> lowercase(x)=="true" ? true :  false , dfd[c])
                        dfd[c] = convert(Array{Bool}, dfd[c])
                    else
                        println("String col not Bool")
                    end
                end
            end
        end
    end
end

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
cmap_tum = ColorScheme(
    reverse(range(
        RGB{Float64}(1,       1,       1),
        RGB{Float64}(0, 101/255, 189/255),
        length=255,
    ))
)

cmap_tum2 = ColorScheme(
    [
        range(
            RGB{Float64}(  0/255,  51/255,  89/255),
            RGB{Float64}(100/255, 160/255, 200/255),
            length=32,
        )
        range(
            RGB{Float64}(100/255, 160/255, 200/255),
            RGB{Float64}(152/255, 198/255, 234/255),
            length=32,
        )
        range(
            RGB{Float64}(152/255, 198/255, 234/255),
            RGB{Float64}(218/255, 215/255, 203/255),
            length=64,
        )
        # range(
        #     RGB{Float64}(218/255, 215/255, 203/255),
        #     RGB{Float64}(162/255, 173/255,   0/255),
        #     length=64,
        # )
        # range(
        #     RGB{Float64}(162/255, 173/255,   0/255),
        #     RGB{Float64}(227/255, 114/255,  34/255),
        #     length=64,
        # )
        range(
            RGB{Float64}(218/255, 215/255, 203/255),
            RGB{Float64}(227/255, 114/255,  34/255),
            length=128,
        )
    ]
)


text:
#MiscTools: Collection of miscellaneous tools under permissive license
#-------------------------------------------------------------------------------
module MiscTools

include("base.jl")

#==Exported symbols
===============================================================================#
export poweroftwo #returns true if integer is a power of two
export writeword #Write data to a word boundary

#==Un-"exported" symbols
================================================================================
	DataWord
==#

#==Other interface tools (symbols not exported to avoid collisions):
================================================================================
	...
==#

end #MiscTools

#Last line


text:
T=PrepTest()

function VSMModel(A::Array{Float64,2},nq::Int)
    DQ=NewDataMatrix()
    DQ.A=A
    DQ.NumQueries=nq
    (ma,na)=size(DQ.A)
    D=DQ.A[:,DQ.NumQueries+1:na]
    Q=DQ.A[:,1:DQ.NumQueries]
    (md,nd)=size(D)
    (mq,nq)=size(Q)
    Costheta_VSM=zeros(nq,nd)
    #Set a tol level
    tol=linspace(0.01,0.99,20)
    lt=length(tol)
    VSM_TDM_Result=TMResults()
    VSM_TDM_Result=InitResults!(VSM_TDM_Result,nq,lt)
    qindex=zeros(nq)
    for z=1:1:nq
        qindex[z]=z
        for i=1:nd
            Costheta_VSM[z:z,i:i]=CosTheta(Q[:,z],D[:,i])
        end
        for k=1:1:lt
            (VSM_TDM_Result.Dr[z,k],
             VSM_TDM_Result.Dt[z,k],
             VSM_TDM_Result.Nr[z,k],
             VSM_TDM_Result.Rec[z,k],
             VSM_TDM_Result.Prec[z,k])=
             FindResults(k,z,Costheta_VSM[z,:],tol[k])
        end
    end
    (VSM_TDM_Result.RecFinal,VSM_TDM_Result.PrecFinal)=average_RecPrec(VSM_TDM_Result.Rec,VSM_TDM_Result.Prec,qindex)
    return VSM_TDM_Result.RecFinal,VSM_TDM_Result.PrecFinal
end

function VSMModel(Q_C::Corpus,D_C::Corpus)
    #Q_c and D_c are the query and Document corpuses.
    PreProcess!(Q_C)
    PreProcess!(D_C)
    #Compulsory update the lexicon, to identify the keywords
    update_lexicon!(Q_C)
    update_lexicon!(D_C)
    #Get the lexicon of the documents.
    DocLex=lexicon(D_C)
    #Cant use the above as query vectors since, the queries must also be wrt
    #Documents lexicon, and hence the dimensions will match
    D=full(tdm(D_C))
    (md,nd)=size(D)
    nq=length(Q_C.documents)
    Q=zeros(md,nq)
    for j=1:nq
        Q[:,j]=dtv(Q_C.documents[j],DocLex)'
    end
    QD=[Q D]
    r=VSMModel(QD,nq)
    return r
end

function VSMModel(QueryNum::Int,A::Array{Float64,2},NumQueries::Int)
    DQ=NewDataMatrix()
    DQ.A=A
    DQ.NumQueries=NumQueries
    (ma,na)=size(DQ.A)
    D=DQ.A[:,DQ.NumQueries+1:1063]
    Q=DQ.A[:,QueryNum]
    Qt=Q'
    (nq,mq)=size(Qt)
    (md,nd)=size(D)
    Costheta_VSM=zeros(nq,nd)
    #Set a tol level
    tol=linspace(0.01,0.99,20)
    lt=length(tol)
    VSM_1Query_Result=TMResults()
    VSM_1Query_Result=InitResults!(VSM_1Query_Result,nq,lt)
    qindex=zeros(nq)
    for z=1:1:nq
        qindex[z]=z
        for i=1:nd
            Costheta_VSM[z:z,i:i]=CosTheta(Q[:,z],D[:,i])
        end
        for k=1:1:lt
            (VSM_1Query_Result.Dr[z,k],
             VSM_1Query_Result.Dt[z,k],
             VSM_1Query_Result.Nr[z,k],
             VSM_1Query_Result.Rec[z,k],
             VSM_1Query_Result.Prec[z,k])=
             FindResults(k,z,Costheta_VSM[z,:],tol[k])
        end
    end
    (VSM_1Query_Result.RecFinal,VSM_1Query_Result.PrecFinal)=average_RecPrec(VSM_1Query_Result.Rec,VSM_1Query_Result.Prec,qindex)
    return VSM_1Query_Result
end


text:
abstract type Model end
struct Zerodominant <: Model end #trait that highlights the regulation mechanism of a type
# struct Onedominant <: Model end
## future work: regulation network for conductances. Onedominant is not used for now.

## trait functions
model(::Type{<:RegIonChannel}) = Zerodominant()
#model(::Type{<:Array{IonChannel}}) = Onedominant()

##trait dispatch
regul_dyn(xs::T, Ca) where {T} = regul_dyn(model(T), xs, Ca)
regul_dyn(::Zerodominant, xs, Ca) = regul_indep(xs, Ca)
#regul_dyn(::Onedominant, xs, Ca) = regul_network(xs, Ca)

function regul_indep(ch::RegIonChannel, Ca)
    g = get_g(ch.ionch)
    R = get_mRna_name(ch)
    τ_ch = ch.τion
    parameters = @parameters $τ_ch Ca_tgt τg
    states = @variables $g(t) $R(t)
    conductance, mRna = states[1], states[2]
    τch = parameters[1]
    eqs = [D(mRna) ~ (Ca_tgt - Ca) / τch,
        D(conductance) ~ (mRna - conductance) / τg]
    defaultmap = [conductance => getfield(ch.ionch, get_g(ch.ionch)), mRna => ch.Rion]
    return eqs, states, parameters, defaultmap
end

#function regul_network(conductances::RegNetwork, Ca)

function Regulated(i::IonChannel)
    if get_name(i) == :Leak
        return i
    else
        R = getfield(i, get_g(i)) #set initial value for R same as g
        return NeuronBuilder.RegIon(i, R, Symbol(:τ, get_name(i)))
    end
end

external_params(ch::RegIonChannel) = (external_params(ch.ionch)..., :Ca_tgt, :τg)
get_mRna_name(ch::RegIonChannel) = Symbol(:mRna, :_, get_name(ch.ionch))
ionic_current(ch::RegIonChannel, sys::ODESystem) = ionic_current(ch.ionch, sys)
calcium_current(ch::RegIonChannel, sys::ODESystem) = calcium_current(ch.ionch, sys)


text:
using JuMP, SCS
const DD = MOD.DefaultDistance()

@testset "Test projections distance on vector sets" begin
    for n in [1, 10] # vector sizes
        v = rand(n)
        for s in (MOI.Zeros(n), MOI.Nonnegatives(n), MOI.Reals(n))
            πv = MOD.projection_on_set(DD, v, s)
            @test MOD.distance_to_set(DD, πv, s) ≈ 0 atol=eps(Float64)
        end
    end
end

@testset "Test projection distance on scalar sets" begin
    v = 5 * randn()
    for s in (MOI.EqualTo(v), MOI.LessThan(v), MOI.GreaterThan(v))
        @test MOD.distance_to_set(DD, v, s) ≈ 0 atol=eps(Float64)
        @test MOD.projection_on_set(DD, v, s) ≈ v
        vm = -v
        πvm = MOD.projection_on_set(DD, vm, s)
        if vm < 0 && !isa(s, MOI.LessThan)
            @test MOD.distance_to_set(DD, vm, s) > 0
            @test abs(vm - πvm) ≈ MOD.distance_to_set(DD, vm, s)
        elseif vm > 0 && !isa(s, MOI.GreaterThan)
            @test MOD.distance_to_set(DD, vm, s) > 0
            @test abs(vm - πvm) ≈ MOD.distance_to_set(DD, vm, s)
        end
    end
end

@testset "Trivial projection on vector cones" begin
    #testing POS
    @test MOD.projection_on_set(DD, -ones(5), MOI.Nonnegatives(5)) ≈ zeros(5)
    @test MOD.projection_on_set(DD, ones(5), MOI.Nonnegatives(5)) ≈ ones(5)
end

@testset "Trivial projection gradient on vector cones" begin
    #testing POS
    @test MOD.projection_gradient_on_set(DD, -ones(5), MOI.Nonnegatives(5)) ≈ zeros(5,5)
    @test MOD.projection_gradient_on_set(DD, ones(5), MOI.Nonnegatives(5)) ≈  Matrix{Float64}(LinearAlgebra.I, 5, 5)

    # testing SOC
    @test MOD.projection_gradient_on_set(DD, [1.0; ones(1)], MOI.SecondOrderCone(1)) ≈ [1.0  0.0;
                                                                                0.0  1.0]
    @test MOD.projection_gradient_on_set(DD, [0.0; ones(4)], MOI.SecondOrderCone(4)) ≈ [0.5   0.25  0.25  0.25  0.25;
                                                                                0.25  0.5   0.0   0.0   0.0;
                                                                                0.25  0.0   0.5   0.0   0.0;
                                                                                0.25  0.0   0.0   0.5   0.0;
                                                                                0.25  0.0   0.0   0.0   0.5]

    # testing SDP trivial
    # eye4 = Matrix{Float64}(LinearAlgebra.I, 4, 4)
    eye4 = [1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]  # symmetrical PSD triangle format
    eye10 = Matrix{Float64}(LinearAlgebra.I, 10, 10)
    @test MOD.projection_gradient_on_set(DD, eye4, MOI.PositiveSemidefiniteConeTriangle(4)) ≈ eye10
end

@testset "Non-trivial joint projection" begin
    c2 = MOI.SecondOrderCone(5)
    v1 = rand(Float64, 10)
    v2 = rand(Float64, 5)
    c1 = MOI.PositiveSemidefiniteConeTriangle(5)
    output_1 = MOD.projection_on_set(DD, v1, c1)
    output_2 = MOD.projection_on_set(DD, v2, c2)
    output_joint = MOD.projection_on_set(DD, [v1, v2], [c1, c2])
    @test output_joint ≈ [output_1; output_2]
end

@testset "Non-trivial Block projection gradient" begin
    v1 = rand(Float64, 15)
    v2 = rand(Float64, 5)
    c1 = MOI.PositiveSemidefiniteConeTriangle(6)
    c2 = MOI.SecondOrderCone(5)

    output_1 = MOD.projection_gradient_on_set(DD, v1, c1)
    output_2 = MOD.projection_gradient_on_set(DD, v2, c2)
    output_joint = MOD.projection_gradient_on_set(DD, [v1, v2], [c1, c2])
    @test output_joint ≈ BlockDiagonal([output_1, output_2])
end

@testset "Exponential Cone Projections" begin
    function det_case_exp_cone(v; dual=false)
        v = dual ? -v : v
        if MOD.distance_to_set(DD, v, MOI.ExponentialCone()) < 1e-8
            return 1
        elseif MOD.distance_to_set(DD, -v, MOI.DualExponentialCone()) < 1e-8
            return 2
        elseif v[1] <= 0 && v[2] <= 0 #TODO: threshold here??
            return 3
        else
            return 4
        end
    end

    function _test_proj_exp_cone_help(x, tol; dual=false)
        cone = dual ? MOI.DualExponentialCone() : MOI.ExponentialCone()
        model = Model()
        set_optimizer(model, optimizer_with_attributes(
            SCS.Optimizer, "eps" => 1e-10, "max_iters" => 10000, "verbose" => 0))
        @variable(model, z[1:3])
        @variable(model, t)
        @objective(model, Min, t)
        @constraint(model, sum((x-z).^2) <= t)
        @constraint(model, z in cone)
        optimize!(model)
        z_star = value.(z)
        px = MOD.projection_on_set(DD, x, cone)
        if !isapprox(px, z_star, atol=tol)
            # error("Exp cone projection failed:\n x = $x\nMOD: $px\nJuMP: $z_star
            #        \nnorm: $(norm(px - z_star))")
            return false
       end
       return true
    end

    Random.seed!(0)
    n = 3
    atol = 1e-7
    case_p = zeros(4)
    case_d = zeros(4)
    for _ in 1:100
        x = randn(3)

        case_p[det_case_exp_cone(x; dual=false)] += 1
        @test _test_proj_exp_cone_help(x, atol; dual=false)

        case_d[det_case_exp_cone(x; dual=true)] += 1
        @test _test_proj_exp_cone_help(x, atol; dual=true)
    end
    @test all(case_p .> 0) && all(case_d .> 0)
end

@testset "Power Cone Projections" begin
    function det_case_pow_cone(x, α; dual=false)
        v = dual ? -x : x
        s = MOI.PowerCone(α)
        if MOD._in_pow_cone(v, s)
            return 1
        elseif MOD._in_pow_cone(v, MOI.dual_set(s))
            return 2
        elseif abs(v[3]) <= 1e-8
            return 3
        else
            return 4
        end
    end

    function _test_proj_pow_cone_help(x, α, tol; dual=false)
        cone = dual ? MOI.DualPowerCone(α) : MOI.PowerCone(α)
        model = Model()
        set_optimizer(model, optimizer_with_attributes(
            SCS.Optimizer, "eps" => 1e-10, "max_iters" => 10000, "verbose" => 0))
        @variable(model, z[1:3])
        @variable(model, t)
        @objective(model, Min, t)
        @constraint(model, sum((x-z).^2) <= t)
        @constraint(model, z in cone)
        optimize!(model)
        z_star = value.(z)
        px = MOD.projection_on_set(DD, x, cone)
        if !isapprox(px, z_star, atol=tol)
            error("x = $x\nα = $α\nnorm = $(norm(px - z_star))\npx=$px\ntrue=$z_star")
            return false
       end
       if !MOD._in_pow_cone(px, cone)
           error("x = $x\nα = $α\nnorm = $(norm(px - z_star))\npx=$px\ntrue=$z_star")
           return false
       end
       return true
    end

    Random.seed!(0)
    n = 3
    atol = 2e-7
    case_p = zeros(4)
    case_d = zeros(4)
    for _ in 1:100
        x = randn(3)
        for α in [rand(0.05:0.05:0.95); 0.5]

            # Need to get some into case 3
            if rand(1:10) == 1
                x[3] = 0
            end

            case_p[det_case_pow_cone(x, α; dual=false)] += 1
            @test _test_proj_pow_cone_help(x, α, atol; dual=false)

            case_d[det_case_pow_cone(x, α; dual=true)] += 1
            @test _test_proj_pow_cone_help(x, α, atol; dual=true)
        end
    end
    @test all(case_p .> 0) && all(case_d .> 0)
end


text:
# function: 　将带ＩＤ的分拆的文件合成为一个文件，并转换成ＰT格式的RAW?文件。

Base.MainInclude.include( "common/fileSplitMerge.jl" )
using Base.Filesystem
###$$$$$$$$$$$$$$$$$$$$
##自动根据输入的文件名，将对应文件名后面加上.001,.002,.003 .....,以此类推，最多到.999的所有文件合并，写入到srcfile指定的目标文件中。
#srcfile="/media/wang/705fc396-8c76-4812-9d0b-d17382d9dfc7/backup/t/app/fiberproc/sjz004/高邑元氏20191126/20191126T033756.lc"
srcfile=ENV["splitedFileRootName"]
try
		;MergeRestorefile(srcfile )
	catch e
		println(e)
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
function err(c)
# err(c)
#
# Interpret the return codes generated by the AMA routines.
#
# The return code c = 2 is used by AMA_schur.m but not by AMA_eig.m.

# Original author: Gary Anderson
# Original file downloaded from:
# http://www.federalreserve.gov/Pubs/oss/oss4/code.html
# Adapted for Dynare by Dynare Team.
#
# This code is in the public domain and may be used freely.
# However the authors would appreciate acknowledgement of the source by
# citation of any of the following papers:
#
# Anderson, G. and Moore, G.
# "A Linear Algebraic Procedure for Solving Linear Perfect Foresight
# Models."
# Economics Letters, 17, 1985.
#
# Anderson, G.
# "Solving Linear Rational Expectations Models: A Horse Race"
# Computational Economics, 2008, vol. 31, issue 2, pages 95-113
#
# Anderson, G.
# "A Reliable and Computationally Efficient Algorithm for Imposing the
# Saddle Point Property in Dynamic Models"
# Journal of Economic Dynamics and Control, 2010, vol. 34, issue 3,
# pages 472-489

    if(c==1)  e="AndersonMoore: unique solution."
elseif(c==2)  e="AndersonMoore: roots not correctly computed by real_schur."
elseif(c==3)  e="AndersonMoore: too many big roots."
elseif(c==35) e="AndersonMoore: too many big roots, and q(:,right) is singular."
elseif(c==4)  e="AndersonMoore: too few big roots."
elseif(c==45) e="AndersonMoore: too few big roots, and q(:,right) is singular."
elseif(c==5)  e="AndersonMoore: q(:,right) is singular."
elseif(c==61) e="AndersonMoore: too many exact shiftrights."
elseif(c==62) e="AndersonMoore: too many numeric shiftrights."
else          e="AndersonMoore: return code not properly specified"
end

return e

end


text:
#=
Given an array of time intervals (start, end) for classroom lectures (possibly overlapping), find the minimum number of rooms required.

For example, given [(30, 75), (0, 50), (60, 150)], you should return 2.
=#

function number_of_rooms_required(lecture_intervals::Array{Tuple{Int64,Int64},1})
    if length(lecture_intervals) == 0
        n_room = 0
    else
        room_reserved = [[0:0]]
        for lecture in lecture_intervals
            room_reserved = reserve_room(lecture, room_reserved)
        end
        n_room = length(room_reserved)
    end
    return n_room
end

function reserve_room(lecture::Tuple{Int64,Int64}, room_reserved::Array{Array{UnitRange{Int64},1},1})
    isbooked = false
    room = 1
    while !isbooked
        if room == length(room_reserved) && !isavailable(lecture, room_reserved[room])
            push!(room_reserved, [lecture[1]:lecture[end]])
            isbooked = true
        elseif isavailable(lecture, room_reserved[room])
            push!(room_reserved[room], lecture[1]:lecture[end])
            sort!(room_reserved[room])
            isbooked = true
        else
            room = room + 1
        end
    end
    return room_reserved
end

function isavailable(lecture::Tuple{Int64,Int64}, room_reserved::Array{UnitRange{Int64},1})
    room_available = false
    for lecture_in_room in 1:length(room_reserved)
        # before the first lecture?
        if lecture_in_room == 1 && minimum(room_reserved[lecture_in_room]) >= lecture[end]
            room_available = true
        elseif lecture_in_room == length(room_reserved) && maximum(room_reserved[lecture_in_room]) <= lecture[1] # after last lecture
            room_available = true
        elseif lecture_in_room != 1 && (minimum(room_reserved[lecture_in_room]) >= lecture[end] && maximum(room_reserved[lecture_in_room-1]) <= lecture[1])
            room_available = true
        end
    end
    return room_available
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem