text:
using TriangularShapedCloudInterpolation, Test
using TriangularShapedCloudInterpolation: find_dim_bounds

function get_test_data()

    resx = 20
    resy = 20
    resz = 20

    pos_disk = zeros(100,3)
    pos_disk[:,1] = [48.27596607669735, 82.56551527383029, 59.674992478017664, 88.72867763530981, 83.74416420754775, 56.553478366219004, 79.40224633281969, 12.247702269226423, 24.76307597016234, 20.65396882675379, 39.739858815215, 53.62475250308947, 42.8799333224674, 99.30427626525548, 40.53984344173627, 98.39331484600464, 27.03099629126937, 17.450630381018904, 92.79680354167115, 1.3200570823053148, 12.117913180641903, 4.562368102053593, 87.02080768524914, 90.12759717061596, 58.105501738661246, 71.30583576368103, 59.18517757543411, 92.4420383395431, 77.6513847683014, 51.83257648918249, 40.932574768040574, 90.29451960422166, 78.40612915836837, 30.234070232039233, 21.904668208057522, 14.710524970038552, 24.965015822432825, 40.531297139835765, 72.08415683854568, 24.951503702876686, 7.952055758881382, 16.384912539215456, 80.23816701271183, 79.59108924359609, 47.50244919274409, 32.58640303654241, 66.67881322926301, 0.834989797604968, 0.09027604518616617, 54.00736832044688, 34.59985432603396, 60.22545028761639, 55.76372588019114, 98.74980384232335, 8.82358921426103, 70.09520172800605, 58.188177320186774, 23.293616266310014, 21.926567293568876, 70.19206131232734, 77.90738195484646, 33.45197083046363, 64.2135701482184, 13.908627644484106, 72.55638246300246, 32.8753721630465, 79.73303856613852, 98.13226337220782, 47.400892329979015, 84.61160817807271, 61.49369466119585, 56.51100321352167, 48.74279436208817, 93.7849625641569, 41.54486173607101, 91.50174126560717, 58.55910247426601, 11.277403796384554, 13.158869661372364, 76.35282256347016, 50.96104199434386, 8.369362954907489, 87.48237787005928, 94.36562189307767, 82.2596843584407, 98.90190443747824, 80.64028033152142, 82.27811763093183, 75.57697393340146, 49.29841329183786, 85.86136846154764, 15.365603830297726, 60.82742244892172, 2.508203108543694, 15.467284525598156, 36.68837053798115, 22.497454082363454, 65.82019419526065, 80.90694344251779, 72.382143332465]
    pos_disk[:,2] = [65.39480109245994, 70.10083828746288, 37.63810246458843, 25.55702502497532, 63.170200867430836, 26.222181031275493, 28.40076500939672, 93.16071978860619, 79.76899457558429, 18.37157757524963, 97.409994245956, 66.22547557732312, 3.6863624270635054, 42.32714426313824, 45.41216920871203, 60.662063203193426, 30.5873288310859, 48.17651887002696, 38.66593628945645, 44.28967590964359, 94.43491548979898, 71.33594937236028, 21.144350285384306, 2.739521704932346, 90.40108553908934, 80.66450870163537, 16.09113302185228, 19.793397425606617, 15.79947110891069, 43.30712856710876, 96.67147150204183, 43.327007552430196, 30.369994647083097, 3.169248891373466, 96.47167676442805, 19.22937093600179, 25.347315843710703, 78.5344212126635, 74.89827103456696, 27.750437567607023, 28.346649837442726, 47.52703983050135, 16.158439001847924, 88.32954929358989, 22.124037524686525, 68.89950571552734, 3.2528438884913635, 83.81769547396011, 73.51371688881345, 59.94019789672696, 81.63970724493946, 39.50110321488998, 13.167313760356691, 58.54206489611129, 92.92978138162051, 28.435264098330993, 47.97945193870354, 18.753540397736003, 4.54053096664464, 58.35089307393031, 64.10459221086693, 40.98522210242628, 43.04009709233625, 58.57504138145766, 86.3711262194778, 96.95621651132566, 40.868012086183136, 11.167131365588245, 12.883082941326052, 37.50868423333902, 33.44071881336519, 87.22231976163893, 92.87972467286691, 98.44586948532299, 23.01710259802792, 0.2786152416903054, 56.98600725065006, 80.50585591815833, 56.05859837567551, 37.004823789584826, 39.19021763001827, 26.493495860890803, 3.6690342929007302, 39.23192946424188, 35.65034152738238, 8.789465141105325, 72.65440204677299, 12.948159607879294, 3.906178425985596, 92.9524301408695, 53.22825001613241, 18.668048147191186, 96.47990059586134, 81.85824145723363, 95.71227928404382, 54.6000896405596, 11.601377078221532, 75.16641020695212, 4.5258015836853405, 49.917360882191545]
    pos_disk[:,3] = [48.564038903676106, 89.45478918925743, 21.750574371868048, 72.97546394350488, 92.40722189668385, 85.8526409023191, 97.64575845692376, 2.3616146277549266, 2.473875222815569, 21.7409350502628, 61.53377198562045, 48.71643708431306, 39.88321122267286, 60.74695452872412, 54.261336771570214, 13.569269939102725, 65.4715854674037, 59.4045857370967, 35.29067595647404, 30.990050711665027, 44.93465265641872, 64.2652654038083, 75.26209525677577, 52.58054850704568, 60.15092623262757, 56.56333940040141, 24.153665440506032, 46.070819063637146, 46.12586347763767, 37.26019208330016, 55.58778846491006, 68.93309729986655, 10.49844710270207, 39.885619439377116, 66.65103605318905, 77.07418628223249, 94.70412849426084, 9.630658643810076, 63.37844491906821, 47.75782956745871, 30.926343723448756, 4.698265680967562, 64.31720472014317, 98.18501866100497, 48.61407936708528, 19.731912607843462, 20.53968076824184, 40.175271529079005, 42.16525702461031, 15.573081317480785, 6.442716102471757, 21.94390759166074, 81.93648367461665, 80.91076751935535, 56.18631620596546, 75.27576968153708, 54.155301361972306, 43.934661356094054, 52.56565635272361, 66.29340009917834, 45.190758753661406, 48.49803399085015, 59.43112288352774, 2.5568744649344888, 69.67714967149436, 57.405196871226295, 18.566925726488904, 80.95842528647228, 65.9537989711509, 23.57553222620201, 4.394211962138184, 27.33168997679849, 51.31394335916977, 33.949265899040086, 74.43815123611886, 81.32197253371727, 7.762577162417839, 16.82475348550152, 93.06333496704026, 19.988963390752023, 94.36145800449829, 99.78290742237421, 5.889244713358632, 41.636983224154456, 46.84819297367977, 34.915457592428226, 47.16026814249889, 29.99799561223737, 91.96293960767035, 76.98025584903725, 32.59174420164748, 9.254245000034311, 54.80050340794167, 12.20906137713631, 51.516614996863154, 67.7763905236716, 44.975328537419415, 37.03769828711438, 93.66992368193591, 85.64487076244889]

    pos_disk_tsc = zeros(100,3)

    minx = minimum(pos_disk[:,1])
    maxx = maximum(pos_disk[:,1]) .* (1.0+1.e-6)
    dx   = -(minx - maxx) / resx
    pos_disk_tsc[:,1] = ( pos_disk[:,1] .- minx ) ./ dx

    minx = minimum(pos_disk[:,2])
    maxx = maximum(pos_disk[:,2]) .* (1.0+1.e-6)
    dx   = -(minx - maxx) / resy
    pos_disk_tsc[:,2] = ( pos_disk[:,2] .- minx ) ./ dx

    minx = minimum(pos_disk[:,3])
    maxx = maximum(pos_disk[:,3]) .* (1.0+1.e-6)
    dx   = -(minx - maxx) / resz
    pos_disk_tsc[:,3] = ( pos_disk[:,3] .- minx ) ./ dx

    rho = [0.1246130632669793, 0.25603749180063606, 0.4104355802162478, 0.5931598708600541, 0.7359315218105704, 0.9689369460686097, 0.38135503989297215, 0.6098051354275797, 0.8966399232205338, 0.675299997710952, 0.8198456567622867, 0.3242744798898247, 0.95224138615869, 0.6453328808947691, 0.6499731928261774, 0.4560801179385685, 0.2612971821651784, 0.8321480898033555, 0.7039184050232261, 0.5742836118900867, 0.16190054048366975, 0.8396157492301837, 0.00495882480188814, 0.6626729098202218, 0.44384137850656913, 0.37192372788756467, 0.01656256851134086, 0.16297499417611583, 0.9276381381659411, 0.8168864360712718, 0.30341799416020554, 0.24432549815597615, 0.2688388653618805, 0.6030398523895424, 0.5972765756216902, 0.24097468402579625, 0.7765414700319739, 0.8530800659297886, 0.5174513300434507, 0.0884473508293464, 0.8708937420043863, 0.8102675598242022, 0.9003113337536319, 0.13741967135155697, 0.09364327246670578, 0.17586275769297766, 0.566290677990132, 0.6358595632217443, 0.1782401573015988, 0.9421812936110101, 0.39430832614614153, 0.9449762810442555, 0.8081868343342873, 0.5158987086852937, 0.9825967999353682, 0.07458178257485804, 0.531105293927749, 0.7644847775769983, 0.3363977184144482, 0.9856472002247592, 0.5497509310270898, 0.8404188015705514, 0.36011013391154734, 0.19751213986592764, 0.874313109870857, 0.7510510536602204, 0.8151290522574428, 0.6669419081022179, 0.3365491392426945, 0.4795433719694546, 0.04286941282476198, 0.4712729226925354, 0.21316914619587846, 0.16241915654940464, 0.39568415093468223, 0.32291797712076664, 0.6234497162571426, 0.3464479715219908, 0.33823090500874864, 0.8940854179779834, 0.016924314962348808, 0.22225644374345754, 0.19775635489355392, 0.3498267055644031, 0.8371335263231685, 0.5777647141166327, 0.7717583996200446, 0.5490746259789394, 0.6046389447280436, 0.4675402645092084, 0.04331278006207917, 0.9302632894778857, 0.8171806161499484, 0.7211135755250679, 0.3483940202513707, 0.005278043113299269, 0.308268002323842, 0.8770268369379157, 0.5503219680779388, 0.23706305557960938]

    return rho, pos_disk_tsc, resx, resy, resz

end

@testset "TSC" begin

    @testset "Dimension bounds" begin

        dimx, dimy, dimz = find_dim_bounds(1)

        @test dimx == 3
        @test dimy == 1
        @test dimz == 1

        dimx, dimy, dimz = find_dim_bounds(2)

        @test dimx == 3
        @test dimy == 3
        @test dimz == 1

        dimx, dimy, dimz = find_dim_bounds(3)

        @test dimx == 3
        @test dimy == 3
        @test dimz == 3

    end

    @testset "TSC positions" begin

        pos = zeros(3,100)
        pos[1,:] = [48.27596607669735, 82.56551527383029, 59.674992478017664, 88.72867763530981, 83.74416420754775, 56.553478366219004, 79.40224633281969, 12.247702269226423, 24.76307597016234, 20.65396882675379, 39.739858815215, 53.62475250308947, 42.8799333224674, 99.30427626525548, 40.53984344173627, 98.39331484600464, 27.03099629126937, 17.450630381018904, 92.79680354167115, 1.3200570823053148, 12.117913180641903, 4.562368102053593, 87.02080768524914, 90.12759717061596, 58.105501738661246, 71.30583576368103, 59.18517757543411, 92.4420383395431, 77.6513847683014, 51.83257648918249, 40.932574768040574, 90.29451960422166, 78.40612915836837, 30.234070232039233, 21.904668208057522, 14.710524970038552, 24.965015822432825, 40.531297139835765, 72.08415683854568, 24.951503702876686, 7.952055758881382, 16.384912539215456, 80.23816701271183, 79.59108924359609, 47.50244919274409, 32.58640303654241, 66.67881322926301, 0.834989797604968, 0.09027604518616617, 54.00736832044688, 34.59985432603396, 60.22545028761639, 55.76372588019114, 98.74980384232335, 8.82358921426103, 70.09520172800605, 58.188177320186774, 23.293616266310014, 21.926567293568876, 70.19206131232734, 77.90738195484646, 33.45197083046363, 64.2135701482184, 13.908627644484106, 72.55638246300246, 32.8753721630465, 79.73303856613852, 98.13226337220782, 47.400892329979015, 84.61160817807271, 61.49369466119585, 56.51100321352167, 48.74279436208817, 93.7849625641569, 41.54486173607101, 91.50174126560717, 58.55910247426601, 11.277403796384554, 13.158869661372364, 76.35282256347016, 50.96104199434386, 8.369362954907489, 87.48237787005928, 94.36562189307767, 82.2596843584407, 98.90190443747824, 80.64028033152142, 82.27811763093183, 75.57697393340146, 49.29841329183786, 85.86136846154764, 15.365603830297726, 60.82742244892172, 2.508203108543694, 15.467284525598156, 36.68837053798115, 22.497454082363454, 65.82019419526065, 80.90694344251779, 72.382143332465]
        pos[2,:] = [65.39480109245994, 70.10083828746288, 37.63810246458843, 25.55702502497532, 63.170200867430836, 26.222181031275493, 28.40076500939672, 93.16071978860619, 79.76899457558429, 18.37157757524963, 97.409994245956, 66.22547557732312, 3.6863624270635054, 42.32714426313824, 45.41216920871203, 60.662063203193426, 30.5873288310859, 48.17651887002696, 38.66593628945645, 44.28967590964359, 94.43491548979898, 71.33594937236028, 21.144350285384306, 2.739521704932346, 90.40108553908934, 80.66450870163537, 16.09113302185228, 19.793397425606617, 15.79947110891069, 43.30712856710876, 96.67147150204183, 43.327007552430196, 30.369994647083097, 3.169248891373466, 96.47167676442805, 19.22937093600179, 25.347315843710703, 78.5344212126635, 74.89827103456696, 27.750437567607023, 28.346649837442726, 47.52703983050135, 16.158439001847924, 88.32954929358989, 22.124037524686525, 68.89950571552734, 3.2528438884913635, 83.81769547396011, 73.51371688881345, 59.94019789672696, 81.63970724493946, 39.50110321488998, 13.167313760356691, 58.54206489611129, 92.92978138162051, 28.435264098330993, 47.97945193870354, 18.753540397736003, 4.54053096664464, 58.35089307393031, 64.10459221086693, 40.98522210242628, 43.04009709233625, 58.57504138145766, 86.3711262194778, 96.95621651132566, 40.868012086183136, 11.167131365588245, 12.883082941326052, 37.50868423333902, 33.44071881336519, 87.22231976163893, 92.87972467286691, 98.44586948532299, 23.01710259802792, 0.2786152416903054, 56.98600725065006, 80.50585591815833, 56.05859837567551, 37.004823789584826, 39.19021763001827, 26.493495860890803, 3.6690342929007302, 39.23192946424188, 35.65034152738238, 8.789465141105325, 72.65440204677299, 12.948159607879294, 3.906178425985596, 92.9524301408695, 53.22825001613241, 18.668048147191186, 96.47990059586134, 81.85824145723363, 95.71227928404382, 54.6000896405596, 11.601377078221532, 75.16641020695212, 4.5258015836853405, 49.917360882191545]
        pos[3,:] = [48.564038903676106, 89.45478918925743, 21.750574371868048, 72.97546394350488, 92.40722189668385, 85.8526409023191, 97.64575845692376, 2.3616146277549266, 2.473875222815569, 21.7409350502628, 61.53377198562045, 48.71643708431306, 39.88321122267286, 60.74695452872412, 54.261336771570214, 13.569269939102725, 65.4715854674037, 59.4045857370967, 35.29067595647404, 30.990050711665027, 44.93465265641872, 64.2652654038083, 75.26209525677577, 52.58054850704568, 60.15092623262757, 56.56333940040141, 24.153665440506032, 46.070819063637146, 46.12586347763767, 37.26019208330016, 55.58778846491006, 68.93309729986655, 10.49844710270207, 39.885619439377116, 66.65103605318905, 77.07418628223249, 94.70412849426084, 9.630658643810076, 63.37844491906821, 47.75782956745871, 30.926343723448756, 4.698265680967562, 64.31720472014317, 98.18501866100497, 48.61407936708528, 19.731912607843462, 20.53968076824184, 40.175271529079005, 42.16525702461031, 15.573081317480785, 6.442716102471757, 21.94390759166074, 81.93648367461665, 80.91076751935535, 56.18631620596546, 75.27576968153708, 54.155301361972306, 43.934661356094054, 52.56565635272361, 66.29340009917834, 45.190758753661406, 48.49803399085015, 59.43112288352774, 2.5568744649344888, 69.67714967149436, 57.405196871226295, 18.566925726488904, 80.95842528647228, 65.9537989711509, 23.57553222620201, 4.394211962138184, 27.33168997679849, 51.31394335916977, 33.949265899040086, 74.43815123611886, 81.32197253371727, 7.762577162417839, 16.82475348550152, 93.06333496704026, 19.988963390752023, 94.36145800449829, 99.78290742237421, 5.889244713358632, 41.636983224154456, 46.84819297367977, 34.915457592428226, 47.16026814249889, 29.99799561223737, 91.96293960767035, 76.98025584903725, 32.59174420164748, 9.254245000034311, 54.80050340794167, 12.20906137713631, 51.516614996863154, 67.7763905236716, 44.975328537419415, 37.03769828711438, 93.66992368193591, 85.64487076244889]

        pos_tsc = get_tsc_positions(pos, 20)

        rho, pos_disk_tsc, resx, resy, resz = get_test_data()

        @test pos_tsc ≈ pos_disk_tsc
    end

    @testset "TSC interpolation 1D" begin

        isolated = false
        wraparound = false

        rho, pos_disk_tsc, resx, resy, resz = get_test_data()

        @test_nowarn TSCInterpolation(rho,
                                pos_disk_tsc[:,1], resx,
                                average=false, isolated=true)

    end

    @testset "TSC interpolation 2D" begin

        isolated = false
        wraparound = false

        rho, pos_disk_tsc, resx, resy, resz = get_test_data()

        @test_nowarn TSCInterpolation(rho,
                                pos_disk_tsc[:,1:2], resx,
                                average=false, isolated=true)

    end

    @testset "TSC interpolation 3D" begin

        isolated = false
        wraparound = false

        rho, pos_disk_tsc, resx, resy, resz = get_test_data()

        tsc = TSCInterpolation(rho,
                                pos_disk_tsc, [resx, resy, resz],
                                average=true)

        @test tsc[ 1,  1,  1] ≈ 0.0
        @test tsc[15, 12, 13] ≈ 0.9856472002247592
        @test tsc[20, 20, 20] ≈ 0.0

        @test_nowarn TSCInterpolation(rho,
                                pos_disk_tsc, [resx, resy, resz],
                                average=false, wraparound=true)

        @test_nowarn TSCInterpolation(rho,
                                pos_disk_tsc, [resx, resy, resz],
                                average=false, isolated=true)

        @test_nowarn TSCInterpolation(rho,
                                pos_disk_tsc, resx,
                                average=false, isolated=true)

        @test_nowarn TSCInterpolation(rho,
                                pos_disk_tsc[:,1], resx,
                                average=false, isolated=true)

    end

end


text:
function _find_parent_proj(currmod)
    parent = parentmodule(currmod)
    (currmod == parent) && return parent
    _is_proj(parent) && return parent
    _find_parent_proj(parent)
end

function gen_sub_proj(currmod::Module, parentmod = _find_parent_proj(currmod))

    # ---------------------------------------------------------------------
    # Check parenthood
    _is_proj(currmod) && error("Module ($(currmod)) is already a ProjAssitant project.")
    !_is_proj(parentmod) && error("Parent module ($(parentmod)) is not a ProjAssitant project.")

    # ---------------------------------------------------------------------
    # Check symbols
    _check_symbols(mod, [:_PROJ_ASSIST_CONF])

    # ---------------------------------------------------------------------
    # top funs
    @eval currmod begin
        const _PROJ_ASSIST_CONF = Dict{Symbol, Any}()
        _PROJ_ASSIST_CONF[:PARENT_PROJ] = $(parentmod)
        _PROJ_ASSIST_CONF[:TOP_PROJ] = $(topproj)($(parentmod))
        nothing
    end
end

# ---------------------------------------------------------------------
macro gen_sub_proj()
    quote $(gen_sub_proj)(@__MODULE__) end
end

macro gen_sub_proj(dirkw)
    # get dir
    k, parent = dirkw.args
    validarg = Meta.isexpr(dirkw, :(=))
    validarg &= (k == :parent) && (parent isa Symbol)
    !validarg && error("An expression `parent=mod::Module` is expected, got ", dirkw)

    quote $(gen_sub_proj)(@__MODULE__, $(esc(parent))) end
end

ext: jl


text:
import MCLightTransport.scatterHG
using Plots

## Plotting Henyey-Greenstein scattering function
rnd = 0.0:0.001:1
plot( rnd, scatterHG.(rnd, 0.99), label="g = 0.99")
plot!(rnd, scatterHG.(rnd, 0.90), label="g = 0.90")
plot!(rnd, scatterHG.(rnd, 0.80), label="g = 0.80")
plot!(rnd, scatterHG.(rnd, 0.50), label="g = 0.50")
plot!(rnd, scatterHG.(rnd, 0.10), label="g = 0.10")
plot!(rnd, scatterHG.(rnd, 0.00), label="g = 0.00")
plot!(legend=:bottomright)

##
plot( rnd[1:end-1], log10.(diff(scatterHG.(rnd, 0.99))), label="g = 0.99")
plot!(rnd[1:end-1], log10.(diff(scatterHG.(rnd, 0.90))), label="g = 0.90")
plot!(rnd[1:end-1], log10.(diff(scatterHG.(rnd, 0.80))), label="g = 0.80")
plot!(rnd[1:end-1], log10.(diff(scatterHG.(rnd, 0.50))), label="g = 0.50")
plot!(rnd[1:end-1], log10.(diff(scatterHG.(rnd, 0.10))), label="g = 0.10")
plot!(rnd[1:end-1], log10.(diff(scatterHG.(rnd, 0.00))), label="g = 0.00")
plot!(legend=:bottomright)

##
henyey_greenstein(cosθ, g) = 0.5*(1 - g^2)/(1 + g^2 - 2*g*cosθ)^(3/2)


text:
# * Seniority

"""
    Seniority(ν)

Seniority is an extra quantum number introduced by Giulio Racah (1943)
to disambiguate between terms belonging to a subshell with a specific
occupancy, that are assigned the same term symbols. For partially
filled f-shells (in ``LS`` coupling) or partially filled ``9/2``
shells (in ``jj`` coupling), seniority alone is not enough to
disambiguate all the arising terms.

The seniority number is defined as the minimum occupancy number `ν ∈
n:-2:0` for which the term first appears, e.g. the ²D term first
occurs in the d¹ configuration, then twice in the d³ configuration
(which will then have the terms ₁²D and ₃²D).
"""
struct Seniority
    ν::Int
end

Base.isless(a::Seniority, b::Seniority) = a.ν < b.ν
Base.iseven(s::Seniority) = iseven(s.ν)
Base.isodd(s::Seniority) = isodd(s.ν)

# This is the notation by Giulio Racah, p.377:
# - Racah, G. (1943). Theory of complex spectra. III. Physical Review,
#   63(9-10), 367–382. http://dx.doi.org/10.1103/physrev.63.367
Base.show(io::IO, s::Seniority) = write(io, to_subscript(s.ν))

istermvalid(term, s::Seniority) =
    iseven(multiplicity(term)) ⊻ iseven(s)

# This is due to the statement "For partially filled f-shells,
# seniority alone is not sufficient to distinguish all possible
# states." on page 24 of
#
# - Froese Fischer, C., Brage, T., & Jönsson, P. (1997). Computational
#   Atomic Structure : An Mchf Approach. Bristol, UK Philadelphia, Penn:
#   Institute of Physics Publ.
#
# This is too strict, i.e. there are partially filled (ℓ ≥ f)-shells
# for which seniority /is/ enough, but I don't know which, so better
# play it safe.
assert_unique_classification(orb::Orbital, occupation, term::Term, s::Seniority) =
    istermvalid(term, s) &&
    (orb.ℓ < 3 || occupation == degeneracy(orb))

function assert_unique_classification(orb::RelativisticOrbital, occupation, J::HalfInteger, s::Seniority)
    ν = s.ν
    # This is deduced by looking at Table A.10 of
    #
    # - Grant, I. P. (2007). Relativistic Quantum Theory of Atoms and
    #   Molecules : Theory and Computation. New York: Springer.
    istermvalid(J, s) &&
        !((orb.j == half(9) && (occupation == 4 || occupation == 6) &&
           (J == 4 || J == 6) && ν == 4) ||
          orb.j > half(9)) # Again, too strict.
end

# * Intermediate terms, seniority
"""
    struct IntermediateTerm{T,S}

Represents a term together with its extra disambiguating quantum number(s), labelled by `ν`.

The term symbol (`::T`) can either be a [`Term`](@ref) (for ``LS``-coupling) or a
`HalfInteger` (for ``jj``-coupling).

The disambiguating quantum number(s) (`::S`) can be anything as long as they are sortable
(i.e. implementing `isless`). It is up to the user to pick a scheme that is suitable for
their application. See "[Disambiguating quantum numbers](@ref)" in the manual for discussion
on how it is used in AtomicLevels.

See also: [`Term`](@ref), [`Seniority`](@ref)

# Constructors

    IntermediateTerm(term, ν)

Constructs an intermediate term with the term symbol `term` and disambiguating quantum
number(s) `ν`.

# Properties

To access the term symbol and the disambiguating quantum number(s), you can use the
`.term :: T` and `.ν :: S` (or `.nu :: S`) properties, respectively. E.g.:

```jldoctest
julia> it = IntermediateTerm(T"2D", 2)
₍₂₎²D

julia> it.term, it.ν
(²D, 2)

julia> it = IntermediateTerm(5//2, Seniority(2))
₂5/2

julia> it.term, it.nu
(5/2, ₂)
```
"""
struct IntermediateTerm{T,S}
    term::T
    ν::S
    IntermediateTerm(term::T, ν::S) where {T<:Union{Term,<:HalfInteger}, S} =
        new{T,S}(term, ν)
    IntermediateTerm(term::Real, ν) =
        IntermediateTerm(convert(HalfInt, term), ν)
end

Base.getproperty(it::IntermediateTerm, s::Symbol) = s == :nu ? getfield(it, :ν) : getfield(it, s)
Base.propertynames(::IntermediateTerm, private::Bool=false) = (:term, :ν, :nu)

function Base.show(io::IO, iterm::IntermediateTerm{<:Any,<:Integer})
    write(io, "₍")
    write(io, to_subscript(iterm.ν))
    write(io, "₎")
    show(io, iterm.term)
end

function Base.show(io::IO, iterm::IntermediateTerm)
    show(io, iterm.ν)
    show(io, iterm.term)
end

Base.isless(a::IntermediateTerm, b::IntermediateTerm) =
    a.ν < b.ν ||
    a.ν == b.ν && a.term < b.term

"""
    intermediate_terms(orb::Orbital, w::Int=one(Int))

Generates all [`IntermediateTerm`](@ref) for a given non-relativstic
orbital `orb` and occupation `w`.

# Examples

```jldoctest
julia> intermediate_terms(o"2p", 2)
3-element Array{IntermediateTerm{Term,Seniority},1}:
 ₀¹S
 ₂¹D
 ₂³P
```

The preceding subscript is the seniority number, which indicates at
which occupancy a certain term is first seen, cf.

```jldoctest
julia> intermediate_terms(o"3d", 1)
1-element Array{IntermediateTerm{Term,Seniority},1}:
 ₁²D

julia> intermediate_terms(o"3d", 3)
8-element Array{IntermediateTerm{Term,Seniority},1}:
 ₁²D
 ₃²P
 ₃²D
 ₃²F
 ₃²G
 ₃²H
 ₃⁴P
 ₃⁴F
```

In the second case, we see both `₁²D` and `₃²D`, since there are two
ways of coupling 3 `d` electrons to a `²D` symmetry.
"""
function intermediate_terms(orb::Union{<:Orbital,<:RelativisticOrbital}, w::Int=one(Int))
    g = degeneracy(orb)
    w > g÷2 && (w = g - w)
    ts = terms(orb, w)
    its = map(unique(ts)) do t
        its = IntermediateTerm{typeof(t),Seniority}[]
        previously_seen = 0

        # We have to loop in reverse, since odd occupation numbers
        # should go from 1 and even from 0.
        for ν ∈ reverse(w:-2:0)
            nn = count_terms(orb, ν, t) - previously_seen
            previously_seen += nn
            append!(its, repeat([IntermediateTerm(t, Seniority(ν))], nn))
        end
        its
    end
    sort(vcat(its...))
end

"""
    intermediate_terms(config)

Generate the intermediate terms for each subshell of `config`.

# Examples

```jldoctest
julia> intermediate_terms(c"1s 2p3")
2-element Array{Array{IntermediateTerm{Term,Seniority},1},1}:
 [₁²S]
 [₁²Pᵒ, ₃²Dᵒ, ₃⁴Sᵒ]

julia> intermediate_terms(rc"3d2 5g3")
2-element Array{Array{IntermediateTerm{HalfIntegers.Half{Int64},Seniority},1},1}:
 [₀0, ₂2, ₂4]
 [₁9/2, ₃3/2, ₃5/2, ₃7/2, ₃9/2, ₃11/2, ₃13/2, ₃15/2, ₃17/2, ₃21/2]
```
"""
function intermediate_terms(config::Configuration)
    map(config) do (orb,occupation,state)
        intermediate_terms(orb,occupation)
    end
end

assert_unique_classification(orb, occupation, it::IntermediateTerm) =
    assert_unique_classification(orb, occupation, it.term, it.ν)

export IntermediateTerm, intermediate_terms, Seniority


text:
# Script to generate Figures 16 and 17
using TestImages, Plots, SparseArrays, LinearAlgebra, Wavelets, MultiscaleGraphSignalTransforms
include("auxilaries.jl")

# Set up the resolution and display size
gr(dpi=200, size=(800,600))

# Load the cameraman image and subsample it to make it as 128x128 image
img = testimage("camera");
matrix = convert(Array{Float64,2}, img)[1:4:512,1:4:512]
m, n = size(matrix)

# Fig. 16a
display(heatmap(matrix, ratio=1, yaxis =:flip, showaxis = false, ticks = false,
    color = :grays, clim = (0,1), colorbar = false))
savefig("cameraman.pdf")
savefig("cameraman.png")

# Now, let's compute the weight matrix based on the Gaussian affinity
# of the small windows (or radius r) around pixels and the pixel locations.
# Set up the key parameters
r = 5; # specify radius of neighbors
σ = 0.007
# Do the weight matrix computation
W007 = image_Gaussian_affinity(matrix, r, σ)

# Preprocess to generate G (GraphSig struct) and GP (GraphPart struct)
# Note that GraphSig requires a matrix data even if it is just a one vector, i.e.,
# f = matrix[:] does not work!
G007 = GraphSig(W007, f = reshape(matrix, (length(matrix), 1)))
GP007 = partition_tree_fiedler(G007)
dmatrix007 = ghwt_analysis!(G007, GP=GP007)

# Construct or search the specific basis
# Haar
BS_haar007 = bs_haar(GP007)
dvec_haar007 = dmatrix2dvec(dmatrix007, GP007, BS_haar007)

# eGHWT
dvec_eghwt007, BS_eghwt007 = eghwt_bestbasis(dmatrix007, GP007)

# Generate Fig. 17a
top_vectors_plot2(dvec_eghwt007, m, n, BS_eghwt007, GP007)
savefig("cameraman_eghwt09_sigma007.pdf")
savefig("cameraman_eghwt09_sigma007.png")

# Now change the σ value
σ = 0.07
W07 = image_Gaussian_affinity(matrix, r, σ)
# Preprocess to generate G (GraphSig struct) and GP (GraphPart struct)
G07 = GraphSig(W07, f = reshape(matrix, (length(matrix), 1)))
GP07 = partition_tree_fiedler(G07)
dmatrix07 = ghwt_analysis!(G07, GP=GP07)

# Construct or search the specific basis
# Haar
BS_haar07 = bs_haar(GP07)
dvec_haar07 = dmatrix2dvec(dmatrix07, GP07, BS_haar07)
# eGHWT
dvec_eghwt07, BS_eghwt07 = eghwt_bestbasis(dmatrix07, GP07)
# Generate Fig. 17b
top_vectors_plot2(dvec_eghwt07, m, n, BS_eghwt07, GP07)
savefig("cameraman_eghwt09_sigma07.pdf")
savefig("cameraman_eghwt09_sigma07.png")
# Finally, do the classical Haar transform in Wavelets.jl
dvec_classichaar = dwt(matrix, wavelet(WT.haar))

# Fig. 16b (Approximation error plot)
DVEC = [ dvec_classichaar[:], dvec_haar07[:], dvec_eghwt07[:],
         dvec_haar007[:], dvec_eghwt007[:] ]
T = [ "Classical Haar", "Graph Haar (σ = 0.07)",
      "eGHWT (σ = 0.07)", "Graph Haar (σ = 0.007)",
      "eGHWT (σ = 0.007)"]
L = [ (:dashdot,:orange), (:dashdot, :red), (:dashdot, :black),
      (:solid, :red), (:solid, :black) ]
approx_error2(DVEC, T, L, 0.5)
savefig("cameraman_approx_error.pdf")
savefig("cameraman_approx_error.png")


text:
let
  # One step of dmrg (dmrg itself
  # leads to an error during precompilation)
  # Add `precompile(Tuple{typeof(dmrg),MPO,MPS,Sweeps})`
  # to precompile.jl by hand.
  N = 4
  sites = siteinds("S=1", N)
  ampo = AutoMPO()
  for j=1:N-1
    ampo .+= ("Sz",j,"Sz",j+1)
    ampo .+= (0.5,"S+",j,"S-",j+1)
    ampo .+= (0.5,"S-",j,"S+",j+1)
  end
  H = MPO(ampo,sites)
  psi0 = randomMPS(sites, 10)
  sweeps = Sweeps(1)
  maxdim!(sweeps, 10)
  cutoff!(sweeps, 1E-11)
  PH = ProjMPO(H)
  which_decomp = "automatic"
  obs = NoObserver()
  quiet = false
  eigsolve_tol = 1e-14
  eigsolve_krylovdim = 3
  eigsolve_maxiter = 1
  eigsolve_verbosity = 0
  ishermitian = true
  eigsolve_which_eigenvalue = :SR
  psi = copy(psi0)
  N = length(psi)
  ITensors.position!(PH, psi0, 1)
  energy = 0.0
  sw = 1
  for sw in 1:nsweep(sweeps)
    # This loop causes problems for precompilation
    #for (b, ha) in sweepnext(N)
    b, ha = 1, 1
      ITensors.position!(PH, psi, b)
      phi = psi[b] * psi[b+1]
      vals, vecs = ITensors.KrylovKit.eigsolve(PH, phi, 1, eigsolve_which_eigenvalue;
                                               ishermitian = ishermitian,
                                               tol = eigsolve_tol,
                                               krylovdim = eigsolve_krylovdim,
                                               maxiter = eigsolve_maxiter)
      energy, phi = vals[1], vecs[1]
      ortho = ha == 1 ? "left" : "right"
      drho = nothing
      spec = replacebond!(psi, b, phi; maxdim = maxdim(sweeps, sw),
                                       mindim = mindim(sweeps, sw),
                                       cutoff = cutoff(sweeps, sw),
                                       eigen_perturbation = drho,
                                       ortho = ortho,
                                       which_decomp = which_decomp)
      measure!(obs; energy = energy,
                    psi = psi,
                    bond = b,
                    sweep = sw,
                    half_sweep = ha,
                    spec = spec,
                    quiet = quiet)
    #end
  end
end


text:
# DESCRIPTION
# Runge-Kutta 4 integration scheme.
#
# INPUT
# - fdyn: name of system dynamics function
# - finp: name of control input function
# - t: time
# - dt: time step
# - dx_curr: current state derivative vector
# - x_curr: current state vector
# - u0: initial control vector
#
# OUTPUT
# - x_new: new state vector after integration
# - dx_new: new state derivative vector after integration
#
# ------------------------------------------------------------------------------

function rk4(fdyn,finp,t,dt,dx_curr,x_curr,u0)
    #
    u=finp(t)+u0
    xd=fdyn(dx_curr,x_curr,u,t)
    dx_new=xd
    k1=dt*xd
    #
    u=finp(t+0.5*dt)+u0
    xd=fdyn(dx_curr,x_curr+0.5*k1,u,t+0.5*dt)
    k2=dt*xd
    #
    xd=fdyn(dx_curr,x_curr+0.5*k2,u,t+0.5*dt)
    k3=dt*xd
    #
    u=finp(t+dt)+u0
    xd=fdyn(dx_curr,x_curr+k3,u,t+dt)
    k4=dt*xd
    # weighted average
    x_new=x_curr+k1/6+k2/3+k3/3+k4/6
    #
    return x_new, dx_new
end

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-GenHel

max_forks_repo_name: umbertosaetti/J-GenHel


text:
#####
##### `ArrayView`
#####

"""
    ArrayView{T} <: AbstractVector{T}

A Julia struct that wraps a MessagePack byte buffer to provide an immutable
view of the MessagePack Array stored within the wrapped byte buffer.

This type is intended to be utilized via [`unpack`](@ref). For example, a call
to `arr = unpack(bytes, ArrayView{Dict{String,Int32}})` will generally return
a value more quickly than `arr = unpack(bytes, Vector{Dict{String,Int32}})`;
the latter will perform full deserialization immediately while the former will
only scan over `bytes` to tag the positions of `arr`'s elements, deferring the
actual deserialization of these elements to the time of their access via
`arr[index]`.

Note that `ArrayView` does not implement any form of caching - repeat accesses
of the same element will re-deserialize the element upon every access.
"""
struct ArrayView{T,B<:AbstractVector{UInt8},S<:Tuple} <: AbstractVector{T}
    bytes::B
    positions::Vector{UInt64}
    strict::S
    ArrayView{T}(bytes::B, positions, strict::S=()) where {T,B,S} = new{T,B,S}(bytes, positions, strict)
end

Base.IndexStyle(::Type{<:ArrayView}) = Base.IndexLinear()

Base.size(arr::ArrayView) = (length(arr.positions),)

Base.@propagate_inbounds function Base.getindex(arr::ArrayView{T}, i::Int) where {T}
    @boundscheck checkbounds(arr, i)
    @inbounds start = arr.positions[i]
    @inbounds stop = i == length(arr) ? length(arr.bytes) : arr.positions[i + 1]
    @inbounds current_bytes = view(arr.bytes, start:stop)
    return unpack(current_bytes, T; strict=arr.strict)
end

#####
##### `MapView`
#####

"""
    MapView{K,V} <: AbstractDict{K,V}

Similar to [`ArrayView`](@ref), but provides an immutable view to a MessagePack
Map rather than a MessagePack Array.

This type is intended to be utilized via [`unpack`](@ref) in the same manner as
`ArrayView`, and is similarly implements a "delay-deserialization-until-access"
mechanism.
"""
struct MapView{K,V,B<:AbstractVector{UInt8},S<:Tuple} <: AbstractDict{K,V}
    bytes::B
    positions::Dict{K,UnitRange{UInt64}}
    strict::S
    MapView{K,V}(bytes::B, positions, strict::S=()) where {K,V,B,S} = new{K,V,B,S}(bytes, positions, strict)
end

Base.length(m::MapView) = length(m.positions)

function _get_by_position(m::MapView{K,V}, position) where {K,V}
    return unpack(view(m.bytes, position), V; strict=m.strict)
end

Base.get(m::MapView, key) = _get_by_position(m, m.positions[key])

Base.haskey(m::MapView, key) = haskey(m.positions, key)

function Base.get(m::MapView, key, default)
    haskey(m, key) && return get(m, key)
    return default
end

function Base.get(default::Base.Callable, m::MapView, key)
    haskey(m, key) && return get(m, key)
    return default()
end

function Base.iterate(m::MapView, prev...)
    x = iterate(m.positions, prev...)
    x isa Nothing && return nothing
    (key, position), state = x
    return key => _get_by_position(m, position), state
end

Base.keys(m::MapView) = keys(m.positions)

Base.values(m::MapView) = (_get_by_position(m, i) for i in values(m.positions))


text:
module FastExpm

using LinearAlgebra
using SparseArrays
using CUDA
using CUDA.CUSPARSE
export fastExpm

"""
    fastExpm(A)
    fastExpm(A; threshold=1e-6)
    fastExpm(A; nonzero_tol=1e-14)
    fastExpm(A; threshold=1e-6, nonzero_tol=1e-14)

 This function efficiently implements matrix exponential for sparse and full matrices.
 This code is based on scaling, taylor series and squaring.
 Currently works only on the CPU

 Two optional keyword arguments are used to speed up the computation and preserve sparsity.
 [1] threshold determines the threshold for the Taylor series (default 1e-6): e.g. fastExpm(A, threshold=1e-10)
 [2] nonzero_tol strips elements smaller than nonzero_tol at each computation step to preserve sparsity (default 1e-14) ,e.g. fastExpm(A, nonzero_tol=1e-10)
 The code automatically switches from sparse to full if sparsity is below 25% to maintain speed.

 This code was originally developed by Ilya Kuprov (http://spindynamics.org/) and has been adapted by F. Mentink-Vigier (fmentink@magnet.fsu.edu)
 and Murari Soundararajan (murari@magnet.fsu.edu)
 If you use this code, please cite
  - H. J. Hogben, M. Krzystyniak, G. T. P. Charnock, P. J. Hore and I. Kuprov, Spinach – A software library for simulation of spin dynamics in large spin systems, J. Magn. Reson., 2011, 208, 179–194.
  - I. Kuprov, Diagonalization-free implementation of spin relaxation theory for large spin systems., J. Magn. Reson., 2011, 209, 31–38.
"""
function fastExpm(A::AbstractMatrix;threshold=1e-6,nonzero_tol=1e-14)
    mat_norm=norm(A,Inf);
    scaling_factor = nextpow(2,mat_norm); # Native routine, faster
    A = A./scaling_factor;
    delta=1;
    rows = LinearAlgebra.checksquare(A); # Throws exception if not square

    # Run Taylor series procedure on the CPU
    if nnz_ext(A)/(rows^2)>0.25 || rows<64
        if issparse(A)
            A=_dense(A);
        end
        P= one(A);#=_dense((1.0+0*im)*I,(rows,rows));=# next_term=P; n=1;
    else
        A=sparse(A);
        P= _sparse(((1.0+0*im)*I,(rows,rows)),A); next_term=P; n=1;
    end

    #show(P)
    while delta>threshold
        # Compute the next term
        if issparse(next_term)
            next_term=(1/n)*A*next_term;
            #Eliminate small elements
            next_term=droptolerance!(next_term, nonzero_tol);
            if nnz_ext(next_term)/length(next_term)>0.25 && issparse(next_term)
                next_term=_dense(next_term);
            end
        else
            next_term=(1/n)*next_term*A;
        end
        delta=norm(next_term,Inf);
        #Add to the total and increment the counter
        P .+= next_term; n=n+1;
    end
    #show(P)
    #Squaring of P to generate correct P
    for n=1:log2(scaling_factor)
        P=P*P;
        if issparse(P)
            if nnz_ext(P)/length(P)<0.25
                P = droptolerance!(P, nonzero_tol);
            else
                P=_dense(P);
            end
        end
    end
    return P
end

function droptolerance!(A::Union{Matrix,CuArray}, tolerance)
    A .= tolerance*round.((1/tolerance).*A)
end
function droptolerance!(A::SparseMatrixCSC, tolerance)
    droptol!(A, tolerance) # Native routine, faster
end

function nnz_ext(A::Matrix)
    count(x->x>0, abs.(A))
end
function nnz_ext(A::SparseMatrixCSC)
    nnz(A) # Native routine, faster
end
function nnz_ext(A::CuArray)
    sum(abs.(A).>0)
end
function nnz_ext(A::AbstractCuSparseMatrix)
    nnz(A) # Native routine, faster
end

function _dense(A::SparseMatrixCSC)
    Matrix(A)
end
function _dense(A::AbstractCuSparseMatrix)
    cu(A)
end
function _sparse(A,B)
    if isa(B,AbstractCuSparseMatrix)
        cu(sparse(A...))
    else
        sparse(A...)
    end
end
function norm(A::CuSparseMatrixCSC,B) #Assume Inf norm
    A = CuSparseMatrixCSR(A)
    columns(n,A) = A[n,:];
    @CUDA.allowscalar rowsums = sum.(columns.(1:size(A,1),(A,)))
    A = CuSparseMatrixCSC(A)
    maximum(rowsums)
end
end # module


text:
function subspace_mul!(st::AbstractVector{T}, comspace, U, subspace, offset=0) where T
    if size(U, 1) == 4
        return subspace_mul4x4!(st, comspace, U, subspace, offset)
    elseif size(U, 1) == 8
        return subspace_mul8x8!(st, comspace, U, subspace, offset)
    else
        # force compiler to specialize on comspace size
        indices = StrideArray{Int}(undef, (StaticInt{length(comspace)}(), ))
        @simd ivdep for i in eachindex(indices)
            indices[i] = comspace[i] + 1
        end
        return subspace_mul_generic!(st, indices, U, subspace, offset)
    end
end

function subspace_mul!(st::AbstractMatrix, comspace, U::AbstractMatrix{T}, subspace, offset=0) where T
    # NOTE: avx version seems to be faster at 8 already
    if size(U, 1) == 4
        return subspace_mul4x4!(st, comspace, U, subspace, offset)
    else
        # force compiler to specialize on comspace size
        indices = StrideArray{Int}(undef, (StaticInt{length(comspace)}(), ))
        @simd ivdep for i in eachindex(indices)
            indices[i] = comspace[i] + 1
        end
        return subspace_mul_generic!(st, indices, U, subspace, offset)
    end
end

# NOTE: avx version won't work with array type that does not provide pointer
# we restrict the state for that to be Array only for now
# thus we need a function to dispatch different methods
function subspace_mul_generic!(st::AbstractVector{T}, indices, U, subspace, offset=0) where T
    # TODO: this should be size(U, 1) x nthreads
    y = similar(st, (size(U, 1), ))
    idx = similar(indices)

    @inbounds @threads for k in subspace
        for i in 1:size(U, 1)
            idx[i] = indices[i] + k + offset
        end

        for i in 1:size(U, 1)
            y[i] = zero(T)
            for j in 1:size(U, 2)
                y[i] += U[i, j] * st[idx[j]]
            end
        end

        for i in 1:size(U, 1)
            st[idx[i]] = y[i]
        end
    end
    return st
end

# NOTE: for other non-standard type, we should manually specialize a subspace routine
#       similar to the complex number case.
# TODO: specialize for HwReal
function subspace_mul_generic!(st::AbstractMatrix{T}, indices, U::AbstractMatrix, subspace, offset=0) where T
    # TODO: this should be size(U, 1) x nthreads
    y = similar(st, (size(U, 1), ))
    idx = similar(indices)

    @inbounds @threads for k in subspace
        for i in 1:size(U, 1)
            idx[i] = k + indices[i] + offset
        end

        for b in 1:size(st, 1)
            for i in 1:size(U, 1)
                y[i] = zero(T)
                for j in 1:size(U, 2)
                    y[i] += U[i, j] * st[b, idx[j]]
                end
            end

            for i in 1:size(U, 1)
                st[b, idx[i]] = y[i]
            end
        end
    end
    return st
end

function subspace_mul_generic!(st::Vector{Complex{T}}, indices, U::AbstractMatrix, subspace, offset=0) where {T <: Base.HWReal}
    D = StrideArrays.static_length(indices)
    y_re = StrideArray{T}(undef, (D, ))
    y_im = StrideArray{T}(undef, (D, ))
    U_re = StrideArray{T}(undef, (D, D))
    U_im = StrideArray{T}(undef, (D, D))
    idx = StrideArray{eltype(indices)}(undef, (D, ))
    Ur = reinterpret(reshape, T, U)

    @inbounds @simd ivdep for i ∈ 1:length(U_re)
        U_re[i] = Ur[2i-1]
        U_im[i] = Ur[2i]
    end

    str = reinterpret(reshape, T, st)

    @inbounds @threads for k in subspace
        for j in 1:size(U, 2)
            idx[j] = indices[j] + k + offset
        end

        @avx for i in axes(U, 1)
            y_re_i = zero(T)
            y_im_i = zero(T)
            for j in axes(U, 2)
                idx_j = idx[j]
                y_re_i += U_re[i, j] * str[1, idx_j] - U_im[i, j] * str[2, idx_j]
                y_im_i += U_re[i, j] * str[2, idx_j] + U_im[i, j] * str[1, idx_j]
            end
            y_re[i] = y_re_i
            y_im[i] = y_im_i
        end

        @avx for i in axes(U, 1)
            idx_i = indices[i] + k + offset
            str[1, idx_i] = y_re[i]
            str[2, idx_i] = y_im[i]
        end
    end
    return st
end

function subspace_mul4x4!(st::AbstractVector{T}, comspace, U, subspace, offset=0) where T
    Base.Cartesian.@nextract 4 indices i -> comspace[i] + 1

    Base.Cartesian.@nextract 4 U i->begin
        Base.Cartesian.@nextract 4 U_i j->U[i, j]
    end

    @inbounds @threads for k in subspace
        Base.Cartesian.@nextract 4 idx i-> k + indices_i + offset

        Base.Cartesian.@nexprs 4 i -> begin
            y_i = zero(T)
            Base.Cartesian.@nexprs 4 j -> begin
                y_i += U_i_j * st[idx_j]
            end
        end

        Base.Cartesian.@nexprs 4 i -> begin
            st[idx_i] = y_i
        end
    end
    return st
end

function subspace_mul8x8!(st::AbstractVector{T}, comspace, U, subspace, offset=0) where T
    Base.Cartesian.@nextract 8 indices i -> comspace[i] + 1

    Base.Cartesian.@nextract 8 U i->begin
        Base.Cartesian.@nextract 8 U_i j->U[i, j]
    end

    @inbounds @threads for k in subspace
        Base.Cartesian.@nextract 8 idx i-> k + indices_i + offset
        Base.Cartesian.@nexprs 8 i -> begin
            y_i = zero(T)
            Base.Cartesian.@nexprs 8 j -> begin
                y_i += U_i_j * st[idx_j]
            end
        end

        Base.Cartesian.@nexprs 8 i -> begin
            st[idx_i] = y_i
        end
    end
    return st
end

function subspace_mul4x4!(st::AbstractMatrix{T}, comspace, U::AbstractMatrix, subspace, offset=0) where T
    Base.Cartesian.@nextract 4 indices i -> comspace[i] + 1

    Base.Cartesian.@nextract 4 U i->begin
        Base.Cartesian.@nextract 4 U_i j->U[i, j]
    end

    @inbounds @threads for k in subspace
        Base.Cartesian.@nextract 4 idx i-> k + indices_i + offset

        for b in 1:size(st, 1)
            Base.Cartesian.@nexprs 4 i -> begin
                y_i = zero(T)
                Base.Cartesian.@nexprs 4 j -> begin
                    y_i += U_i_j * st[b, idx_j]
                end
            end

            Base.Cartesian.@nexprs 4 i -> begin
                st[b, idx_i] = y_i
            end
        end
    end
    return st
end


text:
using Mill, JSON, BSON, Flux, JsonGrinder, Test

using JsonGrinder: DictEntry, suggestextractor, schema
using Mill: reflectinmodel

@testset "basic behavior testing" begin
	j1 = JSON.parse("""{"a": 4, "b": {"a":[1,2,3],"b": 1},"c": { "a": {"a":[1,2,3],"b":[4,5,6]}}}""",inttype=Float64)
	j2 = JSON.parse("""{"a": 4, "c": { "a": {"a":[2,3],"b":[5,6]}}}""")
	j3 = JSON.parse("""{"a": 4, "b": {"a":[1,2,3],"b": 1}}""")
	j4 = JSON.parse("""{"a": 4, "b": {}}""")
	j5 = JSON.parse("""{"b": {}}""")
	j6 = JSON.parse("""{}""")

	sch = JsonGrinder.schema([j1,j2,j3,j4,j5,j6])

	@test sch["a"].counts == Dict(4 => 4)
	@test sch["a"].updated == 4
	@test sch["b"].updated == 4
	@test sch["b"]["a"].updated == 2
	@test sch["b"]["a"].l == Dict(3 => 2)
	@test sch["b"]["a"].items.counts == Dict(1 => 2, 2 => 2, 3 => 2)
	@test sch["b"]["a"].items.updated == 6
	@test sch["b"]["b"].counts == Dict(1 => 2)
	@test sch["b"]["b"].updated == 2
	@test sch["c"].updated == 2
	@test sch["c"]["a"].updated == 2
	@test sch["c"]["a"]["a"].updated == 2
	@test sch["c"]["a"]["a"].l == Dict(2 => 1, 3 => 1)
	@test sch["c"]["a"]["a"].items.counts == Dict(1 => 1, 2 => 2, 3 => 2)
	@test sch["c"]["a"]["a"].items.updated == 5
	@test sch["c"]["a"]["b"].updated == 2
	@test sch["c"]["a"]["b"].l == Dict(2 => 1, 3 => 1)
	@test sch["c"]["a"]["b"].items.counts == Dict(4 => 1, 5 => 2, 6 => 2)
	@test sch["c"]["a"]["b"].items.updated == 5
end

@testset "testing empty arrays" begin
	j1 = JSON.parse("""{"a": []}""")
	j2 = JSON.parse("""{"a": [{"a":1},{"b":2}]}""")
	j3 = JSON.parse("""{"a": [{"a":1,"b":3},{"b":2,"a" : 1}]}""")
	j4 = JSON.parse("""{"a": [{"a":2,"b":3}]}""")

	sch1 = JsonGrinder.schema([j1])
	sch2 = JsonGrinder.schema([j1,j2,j3])
	sch3 = JsonGrinder.schema([j2,j3,j1])

	@test sch1.updated == 1
	@test sch1["a"].updated == 1
	@test isnothing(sch1["a"].items)

	@test sch2.updated == sch3.updated
	@test sch2["a"].l == sch3["a"].l
	@test sch2["a"].updated == sch3["a"].updated
	@test sch2["a"].items["a"].updated == sch3["a"].items["a"].updated
	@test sch2["a"].items["a"].counts == sch3["a"].items["a"].counts
end

@testset "testing schema merging" begin
	j1 = JSON.parse("""{"a": [{"a":1},{"b":2}]}""")
	j2 = JSON.parse("""{"a": [{"a":1,"b":3},{"b":2,"a" : 1}]}""")
	j3 = JSON.parse("""{"a": [{"a":2,"b":3}]}""")
	j4 = JSON.parse("""{"a": []}""")
	j5 = JSON.parse("""{}""")
	j6 = JSON.parse("""{"a": [{"a":1,"b":3},{"b":2,"a" : 1}], "b": 1}""")

	# todo: otestovat jak funguje newentry s víceprvkovám polem
	sch1 = JsonGrinder.schema([j1,j2,j3])
	sch2 = JsonGrinder.schema([j4,j5,j6])

	sch = JsonGrinder.schema([j1,j2,j3,j4,j5,j6])
	sch_merged = merge(sch1, sch2)

	@test sch.updated == sch_merged.updated
	@test sch["a"].updated == sch_merged["a"].updated
	@test sch["a"].l == sch_merged["a"].l
	@test sch["a"].items["a"].updated == sch_merged["a"].items["a"].updated
	@test sch["a"].items["a"].counts == sch_merged["a"].items["a"].counts
	@test sch["a"].items["b"].updated == sch_merged["a"].items["b"].updated
	@test sch["a"].items["b"].counts == sch_merged["a"].items["b"].counts
	@test sch["b"].updated == sch_merged["b"].updated
	@test sch["b"].counts == sch_merged["b"].counts
end

@testset "bson and symbol keys testing" begin
	b1 = IOBuffer()
	j1 = BSON.bson(b1, Dict(:a=>4, :b=>Dict(:a=>[1,2,3], :b=>1), :c=>Dict(:a=>Dict(:a=>[1,2,3],:b=>[4,5,6]))))
	b2 = IOBuffer()
	j2 = BSON.bson(b2, Dict(:a=>4,:c=>Dict(:a=>Dict(:a=>[2,3],:b=>[5,6]))))
	b3 = IOBuffer()
	j3 = BSON.bson(b3, Dict(:a=>4,:b=>Dict(:a=>[1,2,3],:b=>1)))
	b4 = IOBuffer()
	j4 = BSON.bson(b4, Dict(:a=>4,:b=>Dict()))
	b5 = IOBuffer()
	j5 = BSON.bson(b5, Dict(:b=>Dict()))
	b6 = IOBuffer()
	j6 = BSON.bson(b6, Dict())
	bs = [(seek(b, 0); BSON.load(b)) for b in [b1,b2,b3,b4,b5,b6]]
	sch = JsonGrinder.schema(bs)

	@test sch["a"].counts == Dict(4 => 4)
	@test sch["a"].updated == 4
	@test sch["b"].updated == 4
	@test sch["b"]["a"].updated == 2
	@test sch["b"]["a"].l == Dict(3 => 2)
	@test sch["b"]["a"].items.counts == Dict(1 => 2, 2 => 2, 3 => 2)
	@test sch["b"]["a"].items.updated == 6
	@test sch["b"]["b"].counts == Dict(1 => 2)
	@test sch["b"]["b"].updated == 2
	@test sch["c"].updated == 2
	@test sch["c"]["a"].updated == 2
	@test sch["c"]["a"]["a"].updated == 2
	@test sch["c"]["a"]["a"].l == Dict(2 => 1, 3 => 1)
	@test sch["c"]["a"]["a"].items.counts == Dict(1 => 1, 2 => 2, 3 => 2)
	@test sch["c"]["a"]["a"].items.updated == 5
	@test sch["c"]["a"]["b"].updated == 2
	@test sch["c"]["a"]["b"].l == Dict(2 => 1, 3 => 1)
	@test sch["c"]["a"]["b"].items.counts == Dict(4 => 1, 5 => 2, 6 => 2)
	@test sch["c"]["a"]["b"].items.updated == 5
end

@testset "equals and hash test" begin
	j1 = JSON.parse("""{"a": []}""")
	j2 = JSON.parse("""{"a": [{"a":1},{"b":2}]}""")
	j3 = JSON.parse("""{"a": [{"a":1,"b":3},{"b":2,"a" : 1}]}""")
	j4 = JSON.parse("""{"a": [{"a":2,"b":3}]}""")

	sch0 = JsonGrinder.DictEntry()
	sch1 = JsonGrinder.schema([j1, j2, j3, j4])
	sch2 = JsonGrinder.schema([j1, j2, j3, j4])
	sch3 = JsonGrinder.schema([j1, j2, j3])

	@test hash(sch1) === hash(sch2)
	@test hash(sch1) !== hash(sch3)
	@test hash(sch1) !== hash(sch0)
	@test hash(sch3) !== hash(sch0)

	@test sch1 == sch2
	@test sch1 != sch3
	@test sch1 != sch0
	@test sch3 != sch0
end

@testset "show" begin
	j1 = JSON.parse("""{"a": []}""")
	j2 = JSON.parse("""{"a": [{"a":1},{"b":2}]}""")
	j3 = JSON.parse("""{"a": [{"a":1,"b":3},{"b":2,"a" : 1}]}""")
	j4 = JSON.parse("""{"a": [{"a":2,"b":3}], "b": 5}""")

	sch1 = JsonGrinder.schema([j1, j2, j3, j4])
	buf = IOBuffer()
	Base.show(buf, sch1)
	str_repr = String(take!(buf))
	@test str_repr ==
"""
: [Dict] (updated = 4)
  ├── a: [List] (updated = 4)
  │    └── : [Dict] (updated = 5)
  │          ├── a: [Scalar - Int64], 2 unique values, updated = 4
  │          └── b: [Scalar - Int64], 2 unique values, updated = 4
  └── b: [Scalar - Int64], 1 unique values, updated = 1
"""
end

@testset "testing schema merging with max keys" begin
	j1 = JSON.parse("""{"a": [{"a":1},{"b":2}]}""")
	j2 = JSON.parse("""{"a": [{"a":1,"b":3},{"b":2,"a":1}]}""")
	j3 = JSON.parse("""{"a": [{"a":2,"b":3}]}""")
	j4 = JSON.parse("""{"a": []}""")
	j5 = JSON.parse("""{}""")
	j6 = JSON.parse("""{"a": [{"a":1,"b":3},{"b":2,"a":1}], "b":1}""")
	j7 = JSON.parse("""{"a": [{"a":4,"b":5},{"b":6,"a":7}], "b":2}""")
	j8 = JSON.parse("""{"a": [{"a":9,"b":10},{"b":11,"a":12}], "b":2}""")
	j9 = JSON.parse("""{"a": [{"a":4,"b":3},{"b":2,"a":2}], "b":2}""")
	j10 = JSON.parse("""{"a": [{"a":11,"b":12},{"b":13,"a":14}], "b":2}""")
	j11 = JSON.parse("""{"a": [{"a":7,"b":5},{"b":6,"a":6}], "b":2}""")

	JsonGrinder.updatemaxkeys!(6)
	# todo: otestovat jak funguje newentry s víceprvkovám polem
	sch1 = JsonGrinder.schema([j1,j2,j3,j4,j5,j11])
	sch2 = JsonGrinder.schema([j6,j7,j8,j9,j10])

	sch = JsonGrinder.schema([j1,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11])
	sch_merged = merge(sch1, sch2)

	@test sch == sch_merged
end

@testset "Sample synthetic" begin
	j1 = JSON.parse("""{"a": []}""")
	j2 = JSON.parse("""{"a": [{"a":1},{"b":2}]}""")
	j3 = JSON.parse("""{"a": [{"a":1,"b":3},{"b":2,"a" : 1}]}""")
	j4 = JSON.parse("""{"a": [{"a":2,"b":3}]}""")

	sch1 = JsonGrinder.schema([j1, j2, j3, j4])
	@test JsonGrinder.sample_synthetic(sch1) == Dict(:a=>[Dict(:a=>2,:b=>2), Dict(:a=>2,:b=>2)])
end

@testset "Merge empty lists" begin
	j1 = JSON.parse("""{"a": [{"a":1},{"b":2}], "b": []}""")
	j2 = JSON.parse("""{"a": [{"a":3},{"b":4}], "b": []}""")
	j3 = JSON.parse("""{"a": [{"a":1},{"b":3}], "b": []}""")
	j4 = JSON.parse("""{"a": [{"a":2},{"b":4}], "b": [1]}""")

	sch = JsonGrinder.schema([j1, j2, j3, j4])
	sch123 = JsonGrinder.schema([j1, j2, j3])
	sch12 = JsonGrinder.schema([j1, j2])
	sch3 = JsonGrinder.schema([j3])
	sch4 = JsonGrinder.schema([j4])
	sch_merged123 = merge(sch12, sch3)
	sch_merged1234 = merge(sch12, sch3, sch4)
	@test sch == sch_merged1234
	@test sch123 == sch_merged123
end


text:
#Gravity model of commuting in Kenya

#Calculate interaction force for gravity model --- assuming Urban pop
function TupleDist(x1::Tuple{Float64,Float64},x2::Tuple{Float64,Float64})
    return sqrt((x1[1]-x2[1])^2 + (x1[2]-x2[2])^2 )
end

function get_transport_matrix(filename)
    KenyaTbl,n = get_Kenyadata(filename)
    transport_matrix = zeros(n,n)
    for i = 1:n
        for j = 1:n
            if i != j
            transport_matrix[i,j] = KenyaTbl[:Urban][i]*KenyaTbl[:Urban][j]/(TupleDist(KenyaTbl[:Loc][i],KenyaTbl[:Loc][j])^2)
        end
        end
    end
    #Normalise the columns
    for j = 1:n
        transport_matrix[:,j] = transport_matrix[:,j]/sum(transport_matrix[:,j])
    end
    return transport_matrix,n
end

function get_transport_matrix(KenyaTbl,n)
    transport_matrix = zeros(n,n)
    for i = 1:n
        for j = 1:n
            if i != j
            transport_matrix[i,j] = KenyaTbl[:Urban][i]*KenyaTbl[:Urban][j]/(TupleDist(KenyaTbl[:Loc][i],KenyaTbl[:Loc][j])^2)
        end
        end
    end
    #Normalise the columns
    for j = 1:n
        transport_matrix[:,j] = transport_matrix[:,j]/sum(transport_matrix[:,j])
    end
    return transport_matrix
end


text:
# This file is a part of PicoScopes.jl, licensed under the MIT License (MIT).

__precompile__(true)

module PicoScopes

end # module

max_stars_repo_name: mppmu/PicoScopes

max_issues_repo_name: mppmu/PicoScopes

max_forks_repo_name: mppmu/PicoScopes


text:
using QuantumStatistics
using Printf

rtol = [-4, -6, -8, -10, -12]
Λ = [100, 1000, 10000, 100000, 1000000]
# rtol = [-12]
# Λ = [1000]

for lambda in Λ
    for err in rtol
        # dlr = DLR.dlr(:fermi, lambda, 10.0^err)
        dlr = DLR.dlr(:corr, lambda, 10.0^err)
        filename = "basis/dlr$(lambda)_1e$err.dlr"
        open(filename, "w") do io
            for r in 1:length(dlr[:ω])
                @printf(io, "%5i  %32.17g  %32.17g  %8i\n", r, dlr[:ω][r], dlr[:τ][r], dlr[:ωn][r])
            end
        end
    end
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# Methods operating on different special matrix types

# Interconversion between special matrix types
convert{T}(::Type{Bidiagonal}, A::Diagonal{T}) =
    Bidiagonal(A.diag, zeros(T, size(A.diag,1)-1), true)
convert{T}(::Type{SymTridiagonal}, A::Diagonal{T}) =
    SymTridiagonal(A.diag, zeros(T, size(A.diag,1)-1))
convert{T}(::Type{Tridiagonal}, A::Diagonal{T}) =
    Tridiagonal(zeros(T, size(A.diag,1)-1), A.diag, zeros(T, size(A.diag,1)-1))

function convert(::Type{Diagonal}, A::Union{Bidiagonal, SymTridiagonal})
    if !iszero(A.ev)
        throw(ArgumentError("matrix cannot be represented as Diagonal"))
    end
    Diagonal(A.dv)
end

function convert(::Type{SymTridiagonal}, A::Bidiagonal)
    if !iszero(A.ev)
        throw(ArgumentError("matrix cannot be represented as SymTridiagonal"))
    end
    SymTridiagonal(A.dv, A.ev)
end

convert{T}(::Type{Tridiagonal}, A::Bidiagonal{T}) =
    Tridiagonal(A.isupper ? zeros(T, size(A.dv,1)-1) : A.ev, A.dv,
                A.isupper ? A.ev:zeros(T, size(A.dv,1)-1))

function convert(::Type{Bidiagonal}, A::SymTridiagonal)
    if !iszero(A.ev)
        throw(ArgumentError("matrix cannot be represented as Bidiagonal"))
    end
    Bidiagonal(A.dv, A.ev, true)
end

function convert(::Type{Diagonal}, A::Tridiagonal)
    if !(iszero(A.dl) && iszero(A.du))
        throw(ArgumentError("matrix cannot be represented as Diagonal"))
    end
    Diagonal(A.d)
end

function convert(::Type{Bidiagonal}, A::Tridiagonal)
    if iszero(A.dl)
        return Bidiagonal(A.d, A.du, true)
    elseif iszero(A.du)
        return Bidiagonal(A.d, A.dl, false)
    else
        throw(ArgumentError("matrix cannot be represented as Bidiagonal"))
    end
end

function convert(::Type{SymTridiagonal}, A::Tridiagonal)
    if A.dl != A.du
        throw(ArgumentError("matrix cannot be represented as SymTridiagonal"))
    end
    SymTridiagonal(A.d, A.dl)
end

function convert(::Type{Tridiagonal}, A::SymTridiagonal)
    Tridiagonal(copy(A.ev), A.dv, A.ev)
end

function convert(::Type{Diagonal}, A::AbstractTriangular)
    if full(A) != diagm(diag(A))
        throw(ArgumentError("matrix cannot be represented as Diagonal"))
    end
    Diagonal(diag(A))
end

function convert(::Type{Bidiagonal}, A::AbstractTriangular)
    fA = full(A)
    if fA == diagm(diag(A)) + diagm(diag(fA, 1), 1)
        return Bidiagonal(diag(A), diag(fA,1), true)
    elseif fA == diagm(diag(A)) + diagm(diag(fA, -1), -1)
        return Bidiagonal(diag(A), diag(fA,-1), false)
    else
        throw(ArgumentError("matrix cannot be represented as Bidiagonal"))
    end
end

convert(::Type{SymTridiagonal}, A::AbstractTriangular) =
    convert(SymTridiagonal, convert(Tridiagonal, A))

function convert(::Type{Tridiagonal}, A::AbstractTriangular)
    fA = full(A)
    if fA == diagm(diag(A)) + diagm(diag(fA, 1), 1) + diagm(diag(fA, -1), -1)
        return Tridiagonal(diag(fA, -1), diag(A), diag(fA,1))
    else
        throw(ArgumentError("matrix cannot be represented as Tridiagonal"))
    end
end

# Constructs two method definitions taking into account (assumed) commutativity
# e.g. @commutative f{S,T}(x::S, y::T) = x+y is the same is defining
#     f{S,T}(x::S, y::T) = x+y
#     f{S,T}(y::T, x::S) = f(x, y)
macro commutative(myexpr)
    @assert myexpr.head===:(=) || myexpr.head===:function # Make sure it is a function definition
    y = copy(myexpr.args[1].args[2:end])
    reverse!(y)
    reversed_call = Expr(:(=), Expr(:call,myexpr.args[1].args[1],y...), myexpr.args[1])
    esc(Expr(:block, myexpr, reversed_call))
end

for op in (:+, :-)
    SpecialMatrices = [:Diagonal, :Bidiagonal, :Tridiagonal, :Matrix]
    for (idx, matrixtype1) in enumerate(SpecialMatrices) # matrixtype1 is the sparser matrix type
        for matrixtype2 in SpecialMatrices[idx+1:end] # matrixtype2 is the denser matrix type
            @eval begin # TODO quite a few of these conversions are NOT defined
                ($op)(A::($matrixtype1), B::($matrixtype2)) = ($op)(convert(($matrixtype2), A), B)
                ($op)(A::($matrixtype2), B::($matrixtype1)) = ($op)(A, convert(($matrixtype2), B))
            end
        end
    end

    for  matrixtype1 in (:SymTridiagonal,)                      # matrixtype1 is the sparser matrix type
        for matrixtype2 in (:Tridiagonal, :Matrix) # matrixtype2 is the denser matrix type
            @eval begin
                ($op)(A::($matrixtype1), B::($matrixtype2)) = ($op)(convert(($matrixtype2), A), B)
                ($op)(A::($matrixtype2), B::($matrixtype1)) = ($op)(A, convert(($matrixtype2), B))
            end
        end
    end

    for matrixtype1 in (:Diagonal, :Bidiagonal) # matrixtype1 is the sparser matrix type
        for matrixtype2 in (:SymTridiagonal,)   # matrixtype2 is the denser matrix type
            @eval begin
                ($op)(A::($matrixtype1), B::($matrixtype2)) = ($op)(convert(($matrixtype2), A), B)
                ($op)(A::($matrixtype2), B::($matrixtype1)) = ($op)(A, convert(($matrixtype2), B))
            end
        end
    end

    for matrixtype1 in (:Diagonal,)
        for (matrixtype2,matrixtype3) in ((:UpperTriangular,:UpperTriangular),
                                          (:UnitUpperTriangular,:UpperTriangular),
                                          (:LowerTriangular,:LowerTriangular),
                                          (:UnitLowerTriangular,:LowerTriangular))
            @eval begin
                ($op)(A::($matrixtype1), B::($matrixtype2)) = ($op)(($matrixtype3)(A), B)
                ($op)(A::($matrixtype2), B::($matrixtype1)) = ($op)(A, ($matrixtype3)(B))
            end
        end
    end
    for matrixtype in (:SymTridiagonal,:Tridiagonal,:Bidiagonal,:Matrix)
        @eval begin
            ($op)(A::AbstractTriangular, B::($matrixtype)) = ($op)(full(A), B)
            ($op)(A::($matrixtype), B::AbstractTriangular) = ($op)(A, full(B))
        end
    end
end

A_mul_Bc!(A::AbstractTriangular, B::QRCompactWYQ) = A_mul_Bc!(full!(A),B)
A_mul_Bc!(A::AbstractTriangular, B::QRPackedQ) = A_mul_Bc!(full!(A),B)
A_mul_Bc(A::AbstractTriangular, B::Union{QRCompactWYQ,QRPackedQ}) = A_mul_Bc(full(A), B)


text:
# ==========================================================================
# Main interface for linear Support Vector Machines

"""
`svm(X, y⃗; nargs...) → SVM`

Description
============

Fits a Support Vector Machine (SVM) to the given training data `X` and `y⃗`.
Depending on the specified parameters in `nargs`, the SVM can be trained to
perform various forms of regression or classification tasks.
There are a number of parameters that can be configured to the individual needs,
such as the type of lossfunction, the type of regularization, or the algorithm
used for training the SVM in the primal/dual.

Usage
======

    svm(X, y; solver = DualCD(), kernel = ScalarProductKernel(), loss = HingeLoss(), regtype = L2Penalty, bias = 1, C = 1, dual = false, maxiter = 1000, ftol = 1.0e-6, xtol = 1.0e-6, grtol = 1.0e-9, callback = nothing, show_trace = false, nargs...)

    svm(X, Y; nargs...) do t, alpha, v, g
      # this is the callback function
    end

Arguments
==========

- **`X`** : The input matrix of the training set, with the observations as columns.
It can be either a dense or sparse matrix of floatingpoint numbers.
Note that it is generally advised to center and scale your data before training a SVM.

- **`y⃗`** : The target vector of the training set. It has to have the same length as `X` has columns.
In the case of classification, it can be either a vector with elements yᵢ ∈ {-1.0, 1.0}
for two class prediction, or a vector with elements yᵢ ∈ N for multiclass prediction.

- **`solver`** : The concrete algorithm to train the SVM. Depending on the algorithm
this will either solve the linear SVM in the primal or dual formulation.
In contrast to other frameworks, the algorithm does not determine if the solution
of the primal or dual problem is returned (see Details).

- **`kernel`** : The kernel that should be used. The ScalarProductKernel denotes the special
case of a linear SVM and allows for additional functionality and solvers that are specialized
for the linear case.

- **`loss`** : The utilized loss function. The typical loss functions for SVMs are
`L1HingeLoss()` (L1-SVM), `L2HingeLoss()` (L2-SVM), or `SmoothedL1HingeLoss(h)` for classification,
and `EpsilonInsLoss(e)` or `L2EpsilonInsLoss(e)` for support vector regression.
Note that in general each solvers is only able to deal with a subset of those
loss functions and some combinations might not be supported.
For example the implementation of `DualCD` does only support
`HingeLoss` and `L2HingeLoss`.

- **`regtype`** : The type of regularization that should be used. In general this can either
be `L1Penalty` or `L2Penalty`. Note that not all solver support `L1Penalty`.

- **`bias`** : The scaling factor of the bias. If set to 0, no intercept will be fitted.

- **`C`** : The penalty parameter. It plays a role similar to λ⁻¹ in Empirical Risk Minimization.
Conceptually, it denotes the trade off between model complexity and the training error.
Larger C will increase the penalization of training errors and thus lead to behavior
more similar to a hard margin classifier.

- **`dual`** : Boolean that specifies if the solution to the dual problem should be returned
instead of the solution to the primal (default = false).

- **`iterations`** : Specifies the maximum number of iterations before the solver exits with the
current (and probably suboptimal) solution.

- **`ftol`** : Specifies the tolerance for the change of the objective function value.

- **`xtol`** : Specifies the tolerance for the change of the solution.

- **`grtol`** : Specifies the tolerance for the change of the gradient norm.

- **`callback`** : The optional callback function with signature `f(i, w, v, G)`. If a callback is specified, then it will be called
every iteration with the following four parameters in order:
1) `i` … the current iteration number,
2) `w` … the current coefficients (i.e. α if dual=true, or θ if dual=false),
3) `v` … the current objective value
4) `G` … the current gradient.
Note: Some solver support early stopping if the callback function returns the symbol `:Exit`.

- **`show_trace`** : Instead of (or additional to) the callback function, the user can enable the output of the training information to STDOUT.

- **`nargs...`** : Additional named arguments that are passed unchecked to the specified solver.
This functionality can be used to pass around special arguments that the library does not natively implement.

Details
========

In contrast to other frameworks, the `dual` parameter does not influence whether the dual or
the primal problem is solved (this is instead specified implicitly by the `solver` parameter),
but the result the user is interested in.
It is generally recommended to leave `dual = false`, unless the user is either explicitly
interested in the support vectors, or the dataset is sparse and high-dimensional.
Note: Not all solver are able to natively provide both solutions. If the algorithm is unable
to provide the desired solution type then the solution will be converted accordingly, if feasible.

Returns
========

Either a subtype of `DualSVM` or `PrimalSVM`

Author(s)
==========

- Christof Stocker (email: c.stocker@hci-kdd.org)

Examples
=========

### Iris Dataset (adapted example from SVM.jl by John Myles White)

    using KSVM
    using RDatasets
    iris = dataset("datasets", "iris")
    X = convert(Array, iris[1:100, 1:2])'  # The observations have to be in the columns
    y = iris[1:100, :Species]
    train = bitrand(size(X,2))             # Split training and testset
    model = svm(X[:, train], y[train])     # Fit the linear SVM
    acc = accuracy(model, X[:, ~train], y[~train])

### Random Data (adapted example from Regression.jl by Dahua Lin)

    using KSVM
    d = 2             # Number of input dimensions
    n = 100           # Number of pbservations
    w = randn(d+1)    # generate the weight vector
    X = randn(d, n)   # generate input features
    y = sign(X'w[1:d] + w[d+1] + 0.2 * randn(n))  # generate (noisy) response
    model = svm(X, y) # Fit the linear SVM
    y_hat = classify(model, X)

References
===========

- Hsieh, Cho-Jui, et al. "A dual coordinate descent method for large-scale
linear SVM." Proceedings of the 25th international conference on Machine learning.
ACM, 2008. DOI=10.1145/1390156.1390208, http://doi.acm.org/10.1145/1390156.1390208
"""
@inline function svm{TReg<:Penalty}(
    X::AbstractMatrix, y⃗::AbstractVector;
    kernel::Kernel = ScalarProductKernel(),
    solver::SvmSolver = DualCD(),
    loss::Loss = L2HingeLoss(),
    regtype::Type{TReg} = L2Penalty,
    C::Real = 1.,
    nargs...)
  spec = CSVM(kernel = kernel,
              loss = loss,
              regtype = regtype,
              C = C)
  fit(spec, X, y⃗; solver = solver, nargs...)
end

# ==========================================================================
# Special method overload to support the "do" notation for the callback
# res = svm(X, Y) do t, alpha, v, g
#   println("Iteration $t: $v")
# end

@inline function svm(callback::Function,
             X::AbstractMatrix, y⃗::AbstractVector;
             nargs...)
  svm(X, y⃗; callback = callback, nargs...)
end


text:
using ..Frames: Rotation

# F1 -> F2
function transform(ep, rv, _, ::S, ::F1, ::C, ::S, ::F2, ::C) where {S, F1, F2, C}
    rot = Rotation(F1(), F2(), ep)
    State(ep, rot(rv)..., frame=F2(), body=C())
end

# S1 -> S2
function transform(ep, rv, _, ::S1, ::F, ::C, ::S2, ::F, ::C) where {S1, S2, F, C}
    State(Epoch{S2()}(ep), rv..., frame=F(), body=C())
end

# C1 -> C2
function transform(ep, rv, eph, ::S, ::F, ::C1, ::S, ::F, ::C2) where {S, F, C1, C2}
    rot = Rotation(F(), icrf, ep)
    rv_icrf = rot(rv)
    rv_body = state(eph, ep, C1(), C2())
    rv′ = inv(rot)(rv_icrf .- rv_body)
    State(ep, rv′..., frame=F(), body=C2())
end

# F1 -> F2, S1 -> S2
function transform(ep, rv, _, ::S1, ::F1, ::C, ::S2, ::F2, ::C) where {S1, S2, F1, F2, C}
    rot = Rotation(F1(), F2(), ep)
    State(Epoch{S2()}(ep), rot(rv)..., frame=F2(), body=C())
end

# F1 -> F2, C1 -> C2
function transform(ep, rv, eph, ::S, ::F1, ::C1, ::S, ::F2, ::C2) where {S, F1, F2, C1, C2}
    rot1 = Rotation(F1(), icrf, ep)
    rot2 = Rotation(icrf, F2(), ep)
    rv_icrf = rot1(rv)
    rv_body = state(eph, ep, C1(), C2())
    rv′ = rot2(rv_icrf .- rv_body)
    State(ep, rv′..., frame=F2(), body=C2())
end

# S1 -> S2, C1 -> C2
function transform(ep, rv, eph, ::S1, ::F, ::C1, ::S2, ::F, ::C2) where {S1, S2, F, C1, C2}
    rot = Rotation(F(), icrf, ep)
    rv_icrf = rot(rv)
    rv_body = state(eph, ep, C1(), C2())
    rv′ = inv(rot)(rv_icrf .- rv_body)
    State(Epoch{S2()}(ep), rv′..., frame=F(), body=C2())
end

# F1 -> F2, S1 -> S2, C1 -> C2
function transform(ep, rv, eph, ::S1, ::F1, ::C1, ::S2, ::F2, ::C2) where {S1, S2, F1, F2, C1, C2}
    rot1 = Rotation(F1(), icrf, ep)
    rot2 = Rotation(icrf, F2(), ep)
    rv_icrf = rot1(rv)
    rv_body = state(eph, ep, C1(), C2())
    rv′ = rot2(rv_icrf .- rv_body)
    State(Epoch{S2()}(ep), rv′..., frame=F2(), body=C2())
end

const _scale = scale
const _frame = frame
const _body = body

function State(s::AbstractState, eph::AbstractEphemeris;
               frame::AbstractFrame=frame(s),
               scale::TimeScale=scale(s),
               body::CelestialBody=body(s))
    inscale = _scale(s)
    inframe = _frame(s)
    inbody = _body(s)
    inframe == frame && inscale == scale && inbody == body && return s
    ep = epoch(s)
    rv = state(s)
    return transform(ep, rv, eph, inscale, inframe, inbody, scale, frame, body)
end

function KeplerianState(s::AbstractState, eph::AbstractEphemeris;
                        frame::AbstractFrame=frame(s),
                        scale::TimeScale=scale(s),
                        body::CelestialBody=body(s))
    KeplerianState(State(s, eph, frame=frame, scale=scale, body=body))
end


text:
A = sprand(100, 100, 0.2)
Am = Matrix(A)
l = QuadLoss()
rx = QuadReg()
ry = QuadReg()
obs = [(f[1], f[2]) for f in zip(findnz(A)[1], findnz(A)[2])]

gg = GGLRM(A, l, rx, ry, 10, obs=obs)
ggs = GGLRM(A, l, rx, ry, 10, obs=obs)

fit!(gg); fit_sparse!(ggs);
ggobj = whole_objective(gg, gg.X'gg.Y)
ggsobj = whole_objective(ggs, ggs.X'ggs.Y)
@test_approx_eq_eps(ggobj, ggsobj, 0.1*length(A))


text:
"""PART 7: ARITHMETIC WITH SCALED DIMENSIONS.
S99x109_P38x41
"""
module FixPoint

#= export =#
export UNITY, TWO, INF_BAD
export half, round_decimals, mult_and_add, x_over_n, xn_over_d, badness, badness, print_scaled

#= constant =#
"[#101] `2^{16}`, rep 1.00000"
const UNITY =   65536
"[#101] `2^{17}, rep 2.00000"
const TWO =     131072

"[#108]: infinitely bad value"
const INF_BAD = 10000

#= global var =#
"[#104] has arithmetic overflow occurred recently?"
arith_error = false
"[#104] amount subtracted to get an exact division"
tex_remainder = 0

#= functions =#
"""[#105] 455, 716, 1240.
Note: `07777777777L == 1073741823L`
"""
nx_plus_y(n, x, y) = mult_and_add(n, x, y, 0o77_7777_7777)

"""[#105] 1240.
Note: `017777777777L == 2147483647L`
"""
mult_integers(n, x) = mult_and_add(n, x, 0, 0o177_7777_7777)

"""[#100] x div 2

div(x, 2, RoundNearestTiesUp)
"""
function half(x::Integer)
    if isodd(x)
        (x + 1) / 2
    else
        x / 2
    end
end

"""[#102]: create a scaled integer from a given decimal fraction.
{ dig[i], | i = k-1, 0 ≤ k ≤ 17 }
"""
function round_decimals(k, dig) end

"[#103]: prints scaled real, rounded to five digits"
function print_scaled(s) end

"[#105]: return `nx + y`"
function mult_and_add(n, x, y, max_ans) end

"[#106]: `x / n`"
function x_over_n(x, n) end

"[#107]: `x * (n / d)`"
function xn_over_d(x, n, d) end

"""[#108]: compute the “badness” of glue.
given t >= 0
when a total t is supposed to be made from amounts that sum to s.

+ badness = 100 * (t/s)^3
+ badness(t + 1,s) ≥ badness(t,s) ≥ badness(t,s + 1)
"""
function badness(t, s) end

end # FixPoint


text:
push!(LOAD_PATH, "../src/")

using Documenter
using SemiLagrangian
using Plots

ENV["GKSwstype"] = "100"

makedocs(
    modules = [SemiLagrangian],
    sitename = "SemiLagrangian.jl",
    authors = "Yves Mocquard, Pierre Navaro and Nicolas Crouseilles",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        mathengine = MathJax(
            Dict(
                :TeX =>
                    Dict(:equationNumbers => Dict(:autoNumber => "AMS"), :Macros => Dict()),
            ),
        ),
        canonical = "https://juliavlasov.github.io/SemiLagrangian.jl",
        assets = String[],
    ),
    doctest = false,
    pages = [
        "Home" => "index.md",
        "Quickstart" => "quickstart.md",
        "Two dimensions" => "modele_2d.md",
        "API" => "api.md",
        "Types" => "types.md",
        "Functions" => "functions.md",
        "Contents" => "contents.md",
    ],
)

deploydocs(;
    branch = "gh-pages",
    devbranch = "master",
    repo = "github.com/JuliaVlasov/SemiLagrangian.jl",
)


text:
mutable struct SimulationProgressMessenger{T, A, D, Δ} <: Function
    wall_time₀ :: T  # Wall time at simulation start
    wall_time⁻ :: T  # Wall time at previous calback
       adv_cfl :: A
       dif_cfl :: D
            Δt :: Δ
end

SimulationProgressMessenger(Δt) =
    SimulationProgressMessenger(
                      1e-9 * time_ns(),
                      1e-9 * time_ns(),
                      AdvectiveCFL(Δt),
                      DiffusiveCFL(Δt),
                      Δt)

get_Δt(Δt) = Δt
get_Δt(wizard::TimeStepWizard) = wizard.Δt

function (pm::SimulationProgressMessenger)(simulation)
    model = simulation.model

    i, t = model.clock.iteration, model.clock.time

    progress = 100 * (t / simulation.stop_time)

    current_wall_time = 1e-9 * time_ns() - pm.wall_time₀
    time_since_last_callback = 1e-9 * time_ns() - pm.wall_time⁻
    wall_time_per_step = time_since_last_callback / simulation.iteration_interval
    pm.wall_time⁻ = 1e-9 * time_ns()

    u_max = maximum(abs, model.velocities.u)
    v_max = maximum(abs, model.velocities.v)
    w_max = maximum(abs, model.velocities.w)
    ν_max = maximum(abs, model.diffusivities.νₑ)

    @info @sprintf("[%06.2f%%] iteration: % 6d, time: % 10s, Δt: % 10s, wall time: % 8s (% 8s / time step)",
                    progress, i, prettytime(t), prettytime(get_Δt(pm.Δt)), prettytime(current_wall_time), prettytime(wall_time_per_step))

    @info @sprintf("          └── u⃗_max: (%.2e, %.2e, %.2e) m/s, CFL: %.2e, ν_max: %.2e m²/s, νCFL: %.2e",
                    u_max, v_max, w_max, pm.adv_cfl(model), ν_max, pm.dif_cfl(model))

    @info ""

    return nothing
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to conversion from quaternions to Euler angle and axis.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/conversions/quaternion_to_angleaxis.jl
# ==================================================

# Functions: quat_to_angleaxis
# ----------------------------

@testset "Quaternion => Euler angle and axis (Float64)" begin
    T = Float64

    q = Quaternion(cosd(T(75 / 2)), 0, sind(T(75 / 2)), 0)
    av = quat_to_angleaxis(q)
    @test eltype(av) === T
    @test av.a ≈ 75 * pi / 180
    @test av.v ≈ [0, 1, 0]

    q = Quaternion(cosd(T(225 / 2)), 0, sind(T(225 / 2)), 0)
    av = quat_to_angleaxis(q)
    @test eltype(av) === T
    @test av.a ≈ 135 * pi / 180
    @test av.v ≈ [0, -1, 0]
end

@testset "Quaternion => Euler angle and axis (Float32)" begin
    T = Float32

    q = Quaternion(cosd(T(75 / 2)), 0, sind(T(75 / 2)), 0)
    av = quat_to_angleaxis(q)
    @test eltype(av) === T
    @test av.a ≈ T(75 * pi / 180)
    @test av.v ≈ [0, 1, 0]

    q = Quaternion(cosd(T(225 / 2)), 0, sind(T(225 / 2)), 0)
    av = quat_to_angleaxis(q)
    @test eltype(av) === T
    @test av.a ≈ T(135 * pi / 180)
    @test av.v ≈ [0, -1, 0]
end

@testset "Quaternion => Euler angle and axis (Special cases)" begin
    q = Quaternion{Float64}(1, 0, 0, 0)
    av = quat_to_angleaxis(q)
    @test eltype(av) == Float64
    @test av.a == 0
    @test av.v == [0, 0, 0]

    q = Quaternion{Float32}(-1, 0, 0, 0)
    av = quat_to_angleaxis(q)
    @test eltype(av) == Float32
    @test av.a == 0
    @test av.v == [0, 0, 0]
end


text:
using SeisDetective
using SeisIO
using SeisIO.RandSeis
using Statistics
using Test

include("test_trigger.jl")


text:
# Tests for Profitability DEA Models
@testset "ProfitabilityDEAModel" begin

    ## Test Profitability DEA Model with Zofío and Prieto (2006) data
    X = [5 3; 2 4; 4 2; 4 8; 7 9]
    Y = [7 4; 10 8; 8 10; 5 4; 3 6]
    W = [2 1; 2 1; 2 1; 2 1; 2 1.0]
    P = [3 2; 3 2; 3 2; 3 2; 3 2.0]

    # alpha = 0.5 CRS equals Input Oriented CRS
    deaprofbl = deaprofitability(X, Y, W, P)

    @test typeof(deaprofbl) == ProfitabilityDEAModel

    @test nobs(deaprofbl) == 5
    @test ninputs(deaprofbl) == 2
    @test noutputs(deaprofbl) == 2

    @test efficiency(deaprofbl) ≈ [0.388;
                                   1.000;
                                   0.765;
                                   0.250;
                                   0.159] atol = 1e-3
    @test efficiency(deaprofbl, :CRS) ≈ [0.636;
                                1.000;
                                1.000;
                                0.250;
                                0.261] atol = 1e-3
    @test efficiency(deaprofbl, :VRS) ≈ [0.682;
                                1.000;
                                1.000;
                                0.250;
                                0.360] atol = 1e-3
    @test efficiency(deaprofbl, :Scale) ≈ [0.933;
                                1.000;
                                1.000;
                                1.000;
                                0.725] atol = 1e-3
    @test efficiency(deaprofbl, :Allocative) ≈ [0.610;
                                1.000;
                                0.765;
                                1.000;
                                0.609] atol = 1e-3

    @test efficiency(deaprofitability(targets(deaprofbl, :X), targets(deaprofbl, :Y), W, P)) ≈ ones(5) atol=1e-7

    # Check defaults
    @test efficiency(deaprofitability(X, Y, W, P, alpha = 0.5)) == efficiency(deaprofbl)
    @test efficiency(deaprofbl, :Economic) == efficiency(deaprofbl)

    # Print
    show(IOBuffer(), deaprofbl)

    # Test errors
    @test_throws DimensionMismatch deaprofitability([1; 2 ; 3], [4 ; 5], [1; 1; 1], [4; 5]) #  Different number of observations
    @test_throws DimensionMismatch deaprofitability([1; 2; 3], [4; 5; 6], [1; 2; 3; 4], [4; 5; 6]) # Different number of observation in input prices
    @test_throws DimensionMismatch deaprofitability([1; 2; 3], [4; 5; 6], [1; 2; 3], [4; 5; 6; 7]) # Different number of observation in output prices
    @test_throws DimensionMismatch deaprofitability([1 1; 2 2; 3 3], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3], [4; 5; 6]) # Different number of input prices and inputs
    @test_throws DimensionMismatch deaprofitability([1; 2; 3], [4 4; 5 5; 6 6], [1; 2; 3], [4 4 4; 5 5 5; 6 6 6]) # Different number of oputput prices and outputs
    @test_throws ArgumentError efficiency(deaprofbl, :Error)
    @test_throws ArgumentError normfactor(deaprofitability(X, Y, W, P)) # ERROR: ProfitabilityDEAModel has no normalization factor

    # ------------------
    # Test Vector and Matrix inputs and outputs
    # ------------------

    # Inputs is Matrix, Outputs is Vector
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6	8]
    Y = [1; 1; 1; 1; 1; 1; 1; 1]
    W = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1]
    P = [1; 1; 1; 1; 1; 1; 1; 1]

    @test efficiency(deaprofitability(X, Y, W, P)) ≈ ( sum(Y .* P, dims = 2) ./ sum(X .* W, dims = 2) )  / 0.25 atol = 1e-5

    # Inputs is Vector, Output is Matrix
    X = [1; 1; 1; 1; 1; 1; 1; 1]
    Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5]
    W = [1; 1; 1; 1; 1; 1; 1; 1]
    P = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1]

    @test efficiency(deaprofitability(X, Y, W, P)) ≈ ( sum(Y .* P, dims = 2) ./ sum(X .* W, dims = 2) )  / 14 atol = 1e-5

    # Inputs is Vector, Output is Vector
    X = [2; 4; 8; 12; 6; 14; 14; 9.412]
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353]
    W = [1; 1; 1; 1; 1; 1; 1; 1]
    P = [1; 1; 1; 1; 1; 1; 1; 1]

    @test efficiency(deaprofitability(X, Y, W, P)) ≈ ( sum(Y .* P, dims = 2) ./ sum(X .* W, dims = 2) )  / 1.25 atol = 1e-5

end


text:
using Test, MLToolkit.Scripting
using ArgParse: ArgParseSettings, @add_arg_table!

@testset "Scripting" begin
    @testset "Addition for NamedTuple" begin
        t1 = (x=1, y=2)
        t2 = (x=3, y=4)
        tsum = t1 + t2
        @test tsum == (x=4, y=6)
        t = reduce(+, (t1, t2))
        @test t == tsum
    end

    @testset "find_latestdir" begin
        @test DATETIME_FMT == "ddmmyyyy-H-M-S"
        @test find_latestdir(@__DIR__) == "20051994-12-34-56"
    end

    @testset "parse_toml" begin
        tomlpath = joinpath(@__DIR__, "Test.toml")

        argdict = parse_toml(tomlpath, (:level1 => "a1", :level2 => "b1",))
        @test argdict[:f1] == 1
        @test argdict[:f2] == 1
        @test argdict[:f3] == 1

        argdict = parse_toml(tomlpath, (:level1 => "a1", :level2 => "b2",))
        @test argdict[:f1] == 1
        @test argdict[:f2] == 1
        @test argdict[:f3] == 2

        @test_throws AssertionError parse_toml(tomlpath, (:level1 => "a2", :level2 => "b1",))

        argdict = parse_toml(tomlpath, (:level1 => "a2", :level2 => "b1", :level3 => "c1"))
        @test argdict[:f1] == 1
        @test argdict[:f2] == 2
        @test argdict[:f3] == 3
        @test argdict[:f4] == 1

        argdict = parse_toml(tomlpath, (:level1 => "a2", :level2 => "b2",))
        @test argdict[:f1] == 1
        @test argdict[:f2] == 2
        @test argdict[:f3] == 4
    end

    @testset "parse_argstr" begin
        s = ArgParseSettings()

        @add_arg_table! s begin
            "--a"
                arg_type = Int64
                required = true
            "--b"
                arg_type = Float64
                required = true
            "--c"
                arg_type = String
                default = nothing
        end

        argstr = "--a 1 --b 2.0"
        argdict = parse_argstr(argstr, s; as_symbols=true)
        @test argdict[:a] == 1
        @test argdict[:b] == 2.0
        @test argdict[:c] == nothing

        argdict = parse_argstr(argstr * " --c c", s; as_symbols=true)
        @test argdict[:c] == "c"
    end

    @testset "Scripting.check" begin
        x = 1

        @test @jupyter(:(x + 1), default=3) == 2
        @test @script(:(x + 1), default=3) == 3
        @test @tb(:(x + 1), default=3) == 2
        @test @wb(:(x + 1), default=3) == 2

        vc = [1 Inf 3; NaN 5 6]
        vm1 = [1 2 3; 3 4 5]
        vm2 = [11 22 33; 44 55 66]

        check_strs = Scripting.checknumerics_strs(vc, vm1, vm2; vcheckname="vc", vmonitornames=["vm1", "vm2"])
        @test check_strs[1] == "vc[2,1] = NaN\n  vm1[2,1] = 3\n  vm2[2,1] = 44"
        @test check_strs[2] == "vc[1,2] = Inf\n  vm1[1,2] = 2\n  vm2[1,2] = 22"

        @warn "`@checknumerics` is not tested."
    end
end


text:
# wrappers around apf::Field and apf::Numbering creation that deal with field
# naming and recording of which fields are attached to which Julia mesh
# objects (even though they are all attached to the same apf::Mesh object)

#------------------------------------------------------------------------------
# Internal reference counting functions
mutable struct ApfRef
  m::Ptr{Void}  # the mesh
  count::Int
end

global const FieldRefs = Dict{Ptr{Void}, ApfRef}()
global const NumberingRefs = Dict{Ptr{Void}, ApfRef}()

function pushRef(f_ptr::Ptr{Void}, m_ptr::Ptr{Void}, Refs::Dict)

  if haskey(Refs, f_ptr)
    @assert Refs[f_ptr].m == m_ptr
    Refs[f_ptr].count += 1
  else
    Refs[f_ptr] = ApfRef(m_ptr, 1)
  end

  apf.pushMeshRef(m_ptr)  # need to keep the mesh alive so the mesh doesn't try
                          # to free the field

  return nothing
end

function popRef(f_ptr::Ptr{Void}, Refs::Dict, destructor::Function)

  @assert haskey(Refs, f_ptr)

  obj = Refs[f_ptr]
  obj.count -= 1

  if obj.count == 0
    destructor(f_ptr)
    delete!(Refs, f_ptr)
  end

  apf.popMeshRef(obj.m)

  return nothing
end

#------------------------------------------------------------------------------
# Front ends for Fields and Numberings

function pushFieldRef(f_ptr::Ptr{Void}, m_ptr::Ptr{Void})

  pushRef(f_ptr, m_ptr, FieldRefs)

  return nothing
end

function popFieldRef(f_ptr::Ptr{Void})

  popRef(f_ptr, FieldRefs, apf.destroyField)

  return nothing
end

function countFieldRef(f_ptr::Ptr{Void})

  if haskey(FieldRefs, f_ptr)
    return FieldRefs[f_ptr].count
  else
    return 0
  end
end

function pushNumberingRef(f_ptr::Ptr{Void}, m_ptr::Ptr{Void})

  pushRef(f_ptr, m_ptr, NumberingRefs)

  return nothing
end

function popNumberingRef(f_ptr::Ptr{Void})

  popRef(f_ptr, NumberingRefs, apf.destroyNumbering)

  return nothing
end

function countNumberingRef(f_ptr::Ptr{Void})

  if haskey(NumberingRefs, f_ptr)
    return NumberingRefs[f_ptr].count
  else
    return 0
  end
end

#------------------------------------------------------------------------------
# User-visible API
# Extend the apf API replacting the mesh pointer argumetn with a PumiMesh
import apf.createPackedField

"""
  Wrapper around apf.createPackedFields that does the required bookkeeping
  (reference counting and mesh association)

  **Inputs**

   * mesh: PumiMesh object
   * fieldname: the requested name for the field.  If a field with this name
                already exists, a suffix of "_i`, where i is an integer
                will be added to generate a unique field name.
   * ncomponents: number of components on each node
   * fshape: an apf::FieldShape*

  **Outputs**

   * f_ptr: the newly created apf::Field*
"""
function createPackedField(mesh::PumiMesh, fieldname::AbstractString,
                           ncomponents::Integer, fshape::Ptr{Void}=C_NULL)

  # createPackedField does weird things if the name already exists, so find
  # a name that doesn't
  f_ptr = apf.findField(mesh.m_ptr, fieldname)
  if f_ptr != C_NULL
    idx = 1
    while true  # this scares me
      fname2 = string(fieldname, "_", idx)
      f_ptr = apf.findField(mesh.m_ptr, fname2)
      if f_ptr == C_NULL
        fieldname = fname2
        break
      end
      idx += 1
    end
  end

  f_ptr = apf.createPackedField(mesh.m_ptr, fieldname, ncomponents, fshape)

  # add to reference count and list of associated fields
  pushFieldRef(f_ptr, mesh.m_ptr)

  push!(mesh.fields.user, f_ptr)

  return f_ptr
end

"""
  Takes a field that exists on one Julia mesh object and attach it to
  another Julia mesh object as an original field).
  The field will then be attached to both mesh objects at the same time.

  **Inputs**

   * mesh: mesh to attach to
   * f_ptr: the apf::Field*
"""
function attachOrigField(mesh::PumiMesh, f_ptr::Ptr{Void})

  if f_ptr in mesh.fields.orig
    error("cannot attach a field that is already associated with the mesh")
  end

  @assert mesh.m_ptr == apf.getFieldMesh(f_ptr)

  pushFieldRef(f_ptr, mesh.m_ptr)
  push!(mesh.fields.orig, f_ptr)
end

"""
  Like `attachOrigField`, but attaches the field as a user field.
"""
function attachUserField(mesh::PumiMesh, f_ptr::Ptr{Void})

  if f_ptr in mesh.fields.user
    error("cannot attach a field that is already associated with the mesh")
  end

  @assert mesh.m_ptr == apf.getFieldMesh(f_ptr)

  pushFieldRef(f_ptr, mesh.m_ptr)
  push!(mesh.fields.user, f_ptr)
end

"""
  Attach a vector of fields as original fields to a mesh

  **Inputs**

   * mesh: the mesh to attach the fields to
   * fields: AbstractVector of fields.
"""
function attachOrigFields(mesh::PumiMesh, fields::AbstractVector{Ptr{Void}})

  for f in fields
    attachOrigField(mesh, f)
  end

  return nothing
end

import apf.destroyField

function destroyField(mesh::PumiMesh, f_ptr::Ptr{Void})

  if f_ptr == C_NULL
    return nothing
  end

  idx = findfirst(mesh.fields.user, f_ptr)
  arr = mesh.fields.user
  if idx == 0
    idx = findfirst(mesh.fields.orig, f_ptr)
    arr = mesh.fields.orig
  end
  @assert idx != 0

  arr[idx] = C_NULL

  popFieldRef(f_ptr)

  if length(arr) > 50
    shrinkVector(arr, C_NULL)
  end

  return nothing
end

#------------------------------------------------------------------------------
# Numberings

import apf: createNumberingJ, destroyNumbering

function createNumberingJ(mesh::PumiMesh, fieldname::AbstractString,
                           ncomponents::Integer, fshape::Ptr{Void})

  f_ptr = apf.findNumbering(mesh.m_ptr, fieldname)
  if f_ptr != C_NULL
    idx = 1
    while true  # this scares me
      fname2 = string(fieldname, "_", idx)
      f_ptr = apf.findNumbering(mesh.m_ptr, fname2)
      if f_ptr == C_NULL
        fieldname = fname2
        break
      end
      idx += 1
    end
  end

  f_ptr = apf.createNumberingJ(mesh.m_ptr, fieldname, fshape, ncomponents)

  # add to reference count and list of associated fields
  pushNumberingRef(f_ptr, mesh.m_ptr)

  push!(mesh.numberings.user, f_ptr)

  return f_ptr
end

"""
"""
function attachOrigNumbering(mesh::PumiMesh, f_ptr::Ptr{Void})

  if f_ptr in mesh.numberings.orig
    error("cannot attach a field that is already associated with the mesh")
  end

  @assert mesh.m_ptr == apf.getNumberingMesh(f_ptr)

  pushNumberingRef(f_ptr, mesh.m_ptr)
  push!(mesh.numberings.orig, f_ptr)
end

"""
"""
function attachUserNumbering(mesh::PumiMesh, f_ptr::Ptr{Void})

  if f_ptr in mesh.numberings.user
    error("cannot attach a field that is already associated with the mesh")
  end

  @assert mesh.m_ptr == apf.getNumberingMesh(f_ptr)

  pushNumberingRef(f_ptr, mesh.m_ptr)
  push!(mesh.numberings.user, f_ptr)
end

"""
"""
function attachOrigNumberings(mesh::PumiMesh, numberings::AbstractVector{Ptr{Void}})

  for f in numberings
    attachOrigNumbering(mesh, f)
  end

  return nothing
end

import apf.destroyNumbering

function destroyNumbering(mesh::PumiMesh, f_ptr::Ptr{Void})

  if f_ptr == C_NULL
    return nothing
  end

  idx = findfirst(mesh.numberings.user, f_ptr)
  arr = mesh.numberings.user
  if idx == 0
    idx = findfirst(mesh.numberings.orig, f_ptr)
    arr = mesh.numberings.orig
  end
  @assert idx != 0

  arr[idx] = C_NULL

  popNumberingRef(f_ptr)

  if length(arr) > 50
    shrinkVector(arr, C_NULL)
  end

  return nothing
end

#------------------------------------------------------------------------------
# Utilities

"""
  Shrinks a vector by removing all entries that have the given value

  **Inputs**

   * A: the vector
   * val: the value to remove
"""
function shrinkVector(A::AbstractVector, val)

  # indices to keep
  idxs = Vector{Int}(0)
  for i=1:length(A)
    if A[i] != val
      push!(idxs, i)
    end
  end

  for i=1:length(idxs)
    A[i] = A[idx]
  end

  resize!(A, length(idxs))

  return nothing
end

"""
  Adds reference counts for all fields currently attached to the mesh as
  `orig` fields.

  **Inputs**

   * mesh: a PumiMesh object
"""
function recordAllFields(mesh::PumiMesh)

  m_ptr = mesh.m_ptr
  nfields = apf.countFields(m_ptr)
  for i=0:(nfields-1)
    f_ptr = apf.getField(m_ptr, i)
    pushFieldRef(f_ptr, m_ptr)
    push!(mesh.fields.orig, f_ptr)
  end

  return nothing
end

function recordAllNumberings(mesh::PumiMesh)

  m_ptr = mesh.m_ptr
  nfields = apf.countNumberings(m_ptr)
  for i=0:(nfields-1)
    f_ptr = apf.getNumbering(m_ptr, i)
    attachOrigNumbering(mesh, f_ptr)
    #pushNumberingRef(f_ptr, m_ptr)
    #push!(mesh.numberings.orig, f_ptr)
  end

  return nothing
end

import Base.collect

function collect(obj::AttachedData)

  f_arr = Vector{Ptr{Void}}(0)
  for val in obj.orig
    if val != C_NULL
      push!(f_arr, val)
    end
  end

  for val in obj.user
    if val !=  C_NULL
      push!(f_arr, val)
    end
  end

  return f_arr
end


text:
module PiCamera
using PyCall

export picamera, init_pygame, demo

const picamera = PyNULL()
const pygame = PyNULL()
const np = PyNULL()
const PiRGBArray = PyNULL()
const Image = PyNULL()

function __init__()
    pushfirst!(PyVector(pyimport("sys")."path"), "")
    copy!(picamera, pyimport("picamera"))
    copy!(pygame, pyimport("pygame"))
    copy!(np, pyimport("numpy"))
    picamera_array=pyimport("picamera.array")
    copy!(PiRGBArray, picamera_array.PiRGBArray)
    copy!(Image, pyimport("PIL.Image"))
end

# use frame buffer console
ENV["SDL_VIDEODRIVER"] = "fbcon"

function init_pygame()
    pygame.init()
end

py"""
from PIL import Image
import numpy as np

def np2pil(x):
    return Image.fromarray(x.astype(np.uint8))
"""

"""
convert np array to pil data
"""
function np2pil(x)
    py"np2pil"(x)
end

function clip_image(x)
    py"clip_image"(x)
end

"""
display image taken by picamera
"""
function demo()
    pyutils = pyimport("utils")
    labels = pyutils.load_labels()
    classifier = pyimport("classifier")

    multiplier = 1.0
    input_size = 224
    mobilenetversion = classifier.VERSION
    modelname = "mobilenet_$(mobilenetversion)_$(multiplier)_$input_size"
    pyutils.download_checkpoint(multiplier, input_size)
    checkpoint = joinpath(pyutils.SAVEDIR, modelname, "$(modelname).tflite")
    model = classifier.TFModel(checkpoint)

    pygame.init()
    screen_size = (
        pygame.display.Info().current_w,
        pygame.display.Info().current_h
    )
    newsize = (min(screen_size...), min(screen_size...))
    screen = pygame.display.set_mode(screen_size)
    screen.fill((128,128,128)) # set background as gray
    pygame.display.update()
    # set destination where to display image
    dest = [
        (screen_size[1] - newsize[1])/2,
        (screen_size[2] - newsize[2])/2
    ]
    dest = map(x -> floor(Int, x), dest)
    dest[1] -= ceil(Int, dest[2]/3)
    @pywith picamera.PiCamera() as camera begin
        @pywith PiRGBArray(camera) as stream begin
            #camera.resolution = (640, 640)
            start = time()
            while true
                for event in pygame.event.get()
                    if event.type == pygame.KEYDOWN && event.key == pygame.K_ESCAPE
                        pygame.display.quit()
                        pygame.quit()
                    end
                end
                camera.capture(stream, "rgb", use_video_port=true)
                stream.seek(0)

                npImg = classifier.clip_image(stream.array)
                pilImg = np2pil(npImg).resize(newsize)
                top_k, prediction = model(npImg)
                elapsed = time() - start
                @show 1/elapsed
                canvas = classifier.create_canvas(labels, top_k, prediction, newsize[2], elapsed)
                start = time()
                dispImg = hcat([np.array(pilImg), canvas]...)
                dispImg = permutedims(dispImg, (2,1,3))
                img = pygame.surfarray.make_surface(dispImg)
                screen.blit(img, dest)
                stream.truncate(0)
                pygame.display.update()
            end
        end
    end
end
end # end of module

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using aoc2019.password: check, check_strict

first = 240298
last = 784956

p1 = check(first, last)
@assert p1 == 1150
p2 = check_strict(first, last)
@assert p2 == 748

print("-----------------------------------------------------------------------\n")
print("secure container -- part one\n    candidates : $p1\n")
print("secure container -- part two\n    candidates : $p2\n")
print("-----------------------------------------------------------------------\n")


text:
let
    @test C_LIGHT     == 299792458.0
    @test AU          == 1.49597870700e11
    @test MJD_ZERO    == 2400000.5
    @test MJD2000     == 51544.0
    @test GPS_TAI     == -19.0
    @test TAI_GPS     == -GPS_TAI
    @test TT_TAI      == 32.184
    @test TAI_TT      == -TT_TAI
    @test GPS_TT      == GPS_TAI + TAI_TT
    @test TT_GPS      == -GPS_TT
    @test GPS_ZERO    == 44244.0
    @test R_EARTH     == 6.378136300e6
    @test WGS84_a     == 6378137.0
    @test WGS84_f     == 1.0/298.257223563
    @test GM_EARTH    == 3.986004415e14
    @test e_EARTH     == 8.1819190842622e-2
    @test J2_EARTH    == 0.0010826358191967
    @test OMEGA_EARTH == 7.292115146706979e-5
    @test GM_SUN      == 132712440041.939400*1e9
    @test R_SUN       == 6.957*1e8
    @test P_SUN       == 4.560E-6
    @test GM_MOON     == 4902.800066*1e9
    @test GM_MERCURY  == 22031.780000*1e9
    @test GM_VENUS    == 324858.592000*1e9
    @test GM_MARS     == 42828.37521*1e9
    @test GM_JUPITER  == 126712764.8*1e9
    @test GM_SATURN   == 37940585.2*1e9
    @test GM_URANUS   == 5794548.6*1e9
    @test GM_NEPTUNE  == 6836527.100580*1e9
    @test GM_PLUTO    == 977.000000*1e9
end


text:
using Logging
include("yarntests.jl")
open("/tmp/ellytest.log", "a") do io
    global_logger(SimpleLogger(io, Logging.Info))
    run_managed_yarn_clustermanager()
end


text:
using PrinciplesCTImaging: Ellipse, isinside

@testset "points inside ellipse" begin
    e = Ellipse(0.0, 0.0, 1.0, 1.0, 0.0, 1.0)
    @test isinside(e, 0.0, 0.0) == true
    @test isinside(e, 1.0, 0.0) == true
    @test isinside(e, 1.0, 1.0) == false
end

using PrinciplesCTImaging: SheppLoganPhantom, construct

@testset "SheppLogan phantom construction" begin
    pixels = 128
    I = construct(SheppLoganPhantom, pixels)
    @test I[1, 1] ≈ 0.0
    @test I[25, 50] ≈ 1.02
    @test I[50, 50] ≈ 1.0
end