text:
const one_t64      = (1.0, 0.0, 0.0)

const pi_4o1_t64    = (12.56637061435917200,  4.8985871965894130e-16, -1.1979079238873359e-32 )
const pi_2o1_t64    = ( 6.28318530717958600,  2.4492935982947064e-16, -5.9895396194366790e-33 )
const pi_1o1_t64    = ( 3.14159265358979300,  1.2246467991473532e-16, -2.9947698097183397e-33 )
const pi_1o2_t64    = ( 1.57079632679489660,  6.1232339957367660e-17, -1.4973849048591698e-33 )
const pi_15o32_t64  = ( 1.47262155637021560, -5.3616983752483473e-17, -6.3342637055057730e-34 )
const pi_7o16_t64   = ( 1.37444678594553450,  5.3578297462696700e-17,  2.3053216375801512e-34 )
const pi_1o3_t64    = ( 1.04719755119659790, -1.0720817664510910e-16, -9.9825660323944640e-34 )
const pi_1o4_t64    = ( 0.78539816339744830,  3.0616169978683830e-17, -7.4869245242958490e-34 )
const pi_1o6_t64    = ( 0.52359877559829890, -5.3604088322554550e-17, -4.9912830161972320e-34 )
const pi_1o8_t64    = ( 0.39269908169872414,  1.5308084989341915e-17, -3.7434622621479246e-34 )
const pi_1o12_t64   = ( 0.26179938779914946, -2.6802044161277275e-17, -2.4956415080986160e-34 )
const pi_1o16_t64   = ( 0.19634954084936207,  7.6540424946709580e-18, -1.8717311310739623e-34 )
const pi_1o24_t64   = ( 0.13089969389957473, -1.3401022080638637e-17, -1.2478207540493080e-34 )
const pi_1o32_t64   = ( 0.09817477042468103,  3.8270212473354790e-18, -9.3586556553698110e-35 )
const pi_1o64_t64   = ( 0.04908738521234052,  1.9135106236677394e-18, -4.6793278276849057e-35 )
const pi_1o128_t64  = ( 0.02454369260617026,  9.5675531183386970e-19, -2.3396639138424529e-35 )
const pi_1o256_t64  = ( 0.01227184630308513,  4.7837765591693480e-19, -1.1698319569212264e-35 )
const pi_1o512_t64  = ( 0.006135923151542565, 2.3918882795846740e-19, -5.8491597846061320e-36 )
const pi_1o1024_t64 = ( 0.0030679615757712823, 1.195944139792337e-19, -2.9245798923030660e-36 )
const pi_1o2048_t64 = ( 0.0015339807878856412, 5.979720698961686e-20, -1.4622899461515330e-36 )

const inv_pi_4o1_t64    = (   0.07957747154594767, -4.9196691687956215e-18, -2.6803590707232510e-34 )
const inv_pi_2o1_t64    = (   0.15915494309189535, -9.8393383375912430e-18, -5.3607181414465020e-34 )
const inv_pi_1o1_t64    = (   0.31830988618379070, -1.9678676675182486e-17, -1.0721436282893004e-33 )
const inv_pi_1o2_t64    = (   0.63661977236758140, -3.9357353350364970e-17, -2.1442872565786008e-33 )
const inv_pi_15o32_t64  = (   0.67906109052542010,  2.4632204570453420e-17,  9.9968069807037510e-34 )
const inv_pi_7o16_t64   = (   0.72756545413437870,  2.6011543692324530e-18,  1.3563477494445884e-34 )
const inv_pi_1o3_t64    = (   0.95492965855137200, -3.5248787942896340e-18, -1.3494297384832360e-34 )
const inv_pi_1o4_t64    = (   1.27323954473516280, -7.8714706700729940e-17, -4.2885745131572016e-33 )
const inv_pi_1o6_t64    = (   1.90985931710274400, -7.0497575885792670e-18, -2.6988594769664720e-34 )
const inv_pi_1o8_t64    = (   2.54647908947032550, -1.5742941340145989e-16, -8.5771490263144030e-33 )
const inv_pi_1o12_t64   = (   3.81971863420548800, -1.4099515177158535e-17, -5.3977189539329440e-34 )
const inv_pi_1o16_t64   = (   5.09295817894065100, -3.1485882680291977e-16, -1.7154298052628806e-32 )
const inv_pi_1o24_t64   = (   7.63943726841097600, -2.8199030354317070e-17, -1.0795437907865888e-33 )
const inv_pi_1o32_t64   = (  10.18591635788130200, -6.2971765360583950e-16, -3.4308596105257613e-32 )
const inv_pi_1o64_t64   = (  20.37183271576260400, -1.2594353072116790e-15, -6.8617192210515230e-32 )
const inv_pi_1o128_t64  = (  40.74366543152521000, -2.5188706144233580e-15, -1.3723438442103045e-31 )
const inv_pi_1o256_t64  = (  81.48733086305042000, -5.0377412288467160e-15, -2.7446876884206090e-31 )
const inv_pi_1o512_t64  = ( 162.97466172610083000, -1.0075482457693433e-14, -5.4893753768412180e-31 )
const inc_pi_1o1024_t64 = ( 325.94932345220167000, -2.0150964915386866e-14, -1.0978750753682436e-30 )
const inv_pi_1o2048_t64 = ( 651.89864690440330000, -4.0301929830773730e-14, -2.1957501507364872e-30 )

const sin_pi1o3_t64 = (0.8660254037844386, 5.0175421109034514e-17, -7.479771237866948e-34)
const cos_pi1o3_t64 = (0.5, 0.0, 0.0)
const sin_pi1o3_t64 = (1.7320508075688772, 1.0035084221806903e-16, -1.4959542475733896e-33)
const csc_pi1o3_t64 = (1.1547005383792515, 6.690056147871269e-17, -5.105953379741696e-33)
const sec_pi1o3_t64 = (2.0, 0.0, 0.0)
const cot_pi1o3_t64 = (0.5773502691896257, 3.3450280739356345e-17, -2.552976689870848e-33)

const sin_pi1o4_t64 = (0.7071067811865476, -4.833646656726457e-17, 2.0693376543497068e-33)
const cos_pi1o4_t64 = (0.7071067811865476, -4.833646656726457e-17, 2.0693376543497068e-33)
const sin_pi1o4_t64 = (1.0, 0.0, 0.0)
const csc_pi1o4_t64 = (1.4142135623730951, -9.667293313452913e-17, 4.1386753086994136e-33)
const sec_pi1o4_t64 = (1.4142135623730951, -9.667293313452913e-17, 4.1386753086994136e-33)
const cot_pi1o4_t64 = (1.0, 0.0, 0.0)

const sin_pi1o5_t64 = (0.5877852522924731, -7.93475083819002e-18, 3.2552371219323954e-34)
const cos_pi1o5_t64 = (0.8090169943749475, -2.716057601841253e-17, 1.3271626041907827e-33)
const sin_pi1o5_t64 = (0.7265425280053609, -1.3283798985961028e-17, 1.1900616500893826e-34)
const csc_pi1o5_t64 = (1.7013016167040798, 1.0886984134897833e-16, -1.0200333713443962e-33)
const sec_pi1o5_t64 = (1.2360679774997898, -1.0864230407365012e-16, 5.308650416763131e-33)
const cot_pi1o5_t64 = (1.3763819204711736, -4.287034122518378e-17, 9.00390429616504e-34)

const sin_pi1o6_t64 = (0.5, 0.0, 0.0)
const cos_pi1o6_t64 = (0.8660254037844386, 5.0175421109034514e-17, -7.479771237866948e-34)
const sin_pi1o6_t64 = (0.5773502691896257, 3.3450280739356345e-17, -2.552976689870848e-33)
const csc_pi1o6_t64 = (2.0, 0.0, 0.0)
const sec_pi1o6_t64 = (1.1547005383792515, 6.690056147871269e-17, -5.105953379741696e-33)
const cot_pi1o6_t64 = (1.7320508075688772, 1.0035084221806903e-16, -1.4959542475733896e-33)

const sin_pi1o7_t64 = (0.4338837391175581, 7.407189078946677e-20, -1.2028582333265511e-36)
const cos_pi1o7_t64 = (0.9009688679024191, -1.9762646853069492e-17, 9.897202619545834e-34)
const sin_pi1o7_t64 = (0.48157461880752866, -2.056669142326466e-17, -3.854197091963584e-34)
const csc_pi1o7_t64 = (2.3047648709624866, -1.062197674057547e-16, 4.368038423415315e-33)
const sec_pi1o7_t64 = (1.1099162641747424, -4.564997930888251e-17, -5.941617871131375e-34)
const cot_pi1o7_t64 = (2.0765213965723364, 2.0289117907376653e-16, -9.25927028463036e-33)

const sin_pi1o8_t64 = (0.3826834323650898, -1.0050772696461588e-17, -2.0605316302806695e-34)
const cos_pi1o8_t64 = (0.9238795325112867, 1.7645047084336677e-17, -5.044253732158682e-34)
const sin_pi1o8_t64 = (0.41421356237309503, 1.4349369327986523e-17, 1.0571873976798362e-33)
const csc_pi1o8_t64 = (2.613125929752753, -9.583375368676548e-17, -1.4209570724878703e-33)
const sec_pi1o8_t64 = (1.082392200292394, -5.563066290091912e-17, -5.9674442037560255e-34)
const cot_pi1o8_t64 = (2.414213562373095, 1.2537167179050217e-16, 4.1386753086994136e-33)

const sin_pi1o12_t64 = (0.25881904510252074, 2.287249500495561e-17, 1.0253436680177266e-33)
const cos_pi1o12_t64 = (0.9659258262890683, -2.5463971562308955e-17, 1.3193411347856062e-35)
const sin_pi1o12_t64 = (0.2679491924311227, 1.0671460244446628e-17, -4.4789707936399173e-35)
const csc_pi1o12_t64 = (3.8637033051562732, -1.0185588624923582e-16, 5.277364539142425e-35)
const sec_pi1o12_t64 = (1.035276180410083, 9.148998001982243e-17, 4.1013746720709066e-33)
const cot_pi1o12_t64 = (3.732050807568877, 1.0035084221806903e-16, -1.4959542475733896e-33)

const sin_pi1o15_t64 = (0.20791169081775934, -5.47375691962595e-18, -2.1547502853552625e-35)
const cos_pi1o15_t64 = (0.9781476007338057, -5.0904377976839195e-17, -1.154640558932597e-33)
const sin_pi1o15_t64 = (0.21255656167002213, -8.27484570190349e-19, 4.655808820295626e-36)
const csc_pi1o15_t64 = (4.809734344744131, -5.569426258316774e-17, 4.547378632737873e-33)
const sec_pi1o15_t64 = (1.0223405948650293, 1.1456987869876802e-17, 5.239963357056466e-34)
const cot_pi1o15_t64 = (4.704630109478455, -3.8014994127817193e-16, 9.374806675690387e-33)

const sin_pi1o16_t64 = (0.19509032201612828, -7.991079068461731e-18, 6.184627002422071e-34)
const cos_pi1o16_t64 = (0.9807852804032304, 1.8546939997825006e-17, -1.0696564445530757e-33)
const sin_pi1o16_t64 = (0.198912367379658, 8.391794477636538e-19, 2.5564457535700006e-35)
const csc_pi1o16_t64 = (5.125830895483013, -3.8148564901441474e-16, -1.4181052667649092e-32)
const sec_pi1o16_t64 = (1.0195911582083184, -2.2083968146912082e-17, -1.4352691750559562e-33)
const cot_pi1o16_t64 = (5.027339492125848, 2.95379181037367e-17, 2.7177182362115432e-33)

const sin_pi1o24_t64 = (0.1305261922200516, -9.607666023901524e-18, -6.716953849066469e-34)
const cos_pi1o24_t64 = (0.9914448613738104, 2.898989874454537e-17, -5.34967350433339e-34)
const sin_pi1o24_t64 = (0.13165249758739586, -7.917699157902454e-18, 7.983937455025082e-36)
const csc_pi1o24_t64 = (7.661297575540389, -2.5754942653453446e-16, 2.049433239355979e-32)
const sec_pi1o24_t64 = (1.0086289605801526, 1.092498757722754e-16, 1.8625688875143196e-33)
const cot_pi1o24_t64 = (7.59575411272515, -1.505044031166795e-18, -9.502964111090016e-35)

const sin_pi1o36_t64 = (0.08715574274765818, -6.189574214131301e-18, -2.9880893964585044e-34)
const cos_pi1o36_t64 = (0.9961946980917455, -1.2903694855897886e-17, -2.0850330180951652e-34)
const sin_pi1o36_t64 = (0.08748866352592401, -1.7475065678744648e-18, -6.573564498151167e-35)
const csc_pi1o36_t64 = (11.473713245669854, 8.409824359888741e-16, -4.2775072900091694e-32)
const sec_pi1o36_t64 = (1.0038198375433474, 2.988608889329276e-17, -1.592359231061579e-33)
const cot_pi1o36_t64 = (11.430052302761343, 8.11514658899593e-17, -7.42810733044527e-34)

const sin_pi1o48_t64 = (0.06540312923014306, 5.136633854723179e-18, 1.8634981210190512e-34)
const cos_pi1o48_t64 = (0.9978589232386035, 2.7660245814577704e-17, 1.4578690036224136e-33)
const sin_pi1o48_t64 = (0.06554346281523823, -6.244201962707419e-18, 3.670976186047138e-34)
const csc_pi1o48_t64 = (15.289788298678511, 7.267038125576664e-16, 2.6974131728868103e-32)
const sec_pi1o48_t64 = (1.0021456708072995, 6.035414166361829e-17, 2.9457074038261743e-33)
const cot_pi1o48_t64 = (15.25705168826554, -2.5905447056570125e-16, 1.9051151791377825e-32)

const sin_pi1o64_t64 = (0.049067674327418015, -6.79610372051828e-19, -4.4318868124718325e-35)
const cos_pi1o64_t64 = (0.9987954562051724, -1.2291693337075465e-17, 2.446844678649127e-34)
const sin_pi1o64_t64 = (0.049126849769467254, 9.097765655528944e-20, -4.7470025699432424e-36)
const csc_pi1o64_t64 = (20.380016247096112, 1.6969944228826224e-15, -2.1307108889014292e-32)
const sec_pi1o64_t64 = (1.0012059964703925, 5.748536721404889e-17, 2.6420799527691316e-33)
const cot_pi1o64_t64 = (20.355467624987188, -2.3792881581892444e-17, 5.288182696219774e-34)

const sin_pi1o128_t64 = (0.024541228522912288, -9.186849012577878e-20, 4.870614870446706e-36)
const cos_pi1o128_t64 = (0.9996988186962042, -2.985148640379975e-17, -1.9084787370733737e-33)
const sin_pi1o128_t64 = (0.024548622108925444, -5.838370447784443e-20, -6.964141459827823e-37)
const csc_pi1o128_t64 = (40.74775633446287, -2.5217154537342434e-15, -1.4049982810542148e-31)
const sec_pi1o128_t64 = (1.000301272041302, -6.001229236340269e-17, -3.315945536296742e-33)
const cot_pi1o128_t64 = (40.7354838720833, 1.6732015413007298e-15, 6.550337088915585e-32)

const sin_pi1o512_t64 = (0.006135884649154475, 9.054525748247493e-20, 1.626011313374532e-37)
const cos_pi1o512_t64 = (0.9999811752826011, 3.3568103522895585e-17, -1.4740132559368063e-35)
const sin_pi1o512_t64 = (0.006136000157623402, -1.746345263565922e-19, 4.277381971459752e-36)
const csc_pi1o512_t64 = (162.9756843844514, -1.222532369546018e-14, 1.3070853382863645e-31)
const sec_pi1o512_t64 = (1.0000188250717754, 9.449694210745273e-17, -7.377442746575774e-34)
const cot_pi1o512_t64 = (162.97261641324997, -3.0164596628737455e-15, -1.55236796149015e-31)

const sin_pi1o1024_t64 = (0.003067956762965976, 1.2690279085455925e-19, 5.287946424532839e-36)
const cos_pi1o1024_t64 = (0.9999952938095762, -1.966806428532219e-17, -6.305395509588348e-34)
const sin_pi1o1024_t64 = (0.003067971201422665, -8.446036355314405e-20, 4.918546909029024e-36)
const csc_pi1o1024_t64 = (325.9498347796924, 1.4061463896494086e-14, 4.991692857301745e-31)
const sec_pi1o1024_t64 = (1.000004706212572, 6.900021265655887e-17, -4.061965915398279e-33)
const cot_pi1o1024_t64 = (325.94830079770134, 1.3179926072070083e-14, -2.4528262320378532e-32)

const sin_pi1o2048_t64 = (0.0015339801862847657, -1.0467712971596958e-19, 2.4745604308673836e-36)
const cos_pi1o2048_t64 = (0.9999988234517019, 3.067881716126108e-17, -2.6247039963835992e-33)
const sin_pi1o2048_t64 = (0.0015339819910886664, 8.145819328928074e-20, 1.0849852064036664e-36)
const csc_pi1o2048_t64 = (651.8989025679381, 4.950041690164551e-14, 1.4246444414492298e-30)
const sec_pi1o2048_t64 = (1.0000011765496823, 1.0780476438354186e-17, -2.862837848637511e-35)
const cot_pi1o2048_t64 = (651.8981355773938, 2.724138996856417e-14, 4.74641023409796e-31)


text:
# TODO: Test arbitrary profile waveguide, for now have only tested simple waveguide
module Waveguides

export add_halfspace_waveguide,
add_halfspace_waveguides,
add_planar_waveguide,
add_planar_waveguides,
add_pc_waveguide,
add_pc_waveguides

using ...Bravais,
...BoundaryConditions,
...Shapes,
...DefineSimulation,
..ConstructionToolsBase,
..PhotonicCrystal

include("Waveguides/waveguides_planar.jl")
include("Waveguides/waveguides_halfspace.jl")
include("Waveguides/waveguides_pc.jl")

include("DefineSimulation/iss.jl")

end # module

ext: jl


text:
using ISA
using Test

@testset "ISA.jl" begin
    # Write your own tests here.
end

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA


text:
# This example simulates cluster synchronization using pinning control

import Cairo, Fontconfig
using LinearAlgebra
using DifferentialEquations
using PyPlot
using BlockArrays
using Graphs, GraphPlot
pygui(true)
close("all")

# Helper functions
getwave(x, i, d) = getindex.(x, (i - 1) * d + 1)

# Define node dynamics
h(x, a = -1 / 7, b = 2 / 7) = b * x + 1 / 2 * (a - b) * (abs(x + 1) - abs(x - 1))
function chua(dx, x, α = 9, β = 100 / 7, h = h)
    dx[1] = α * (x[2] - h(x[1]))
    dx[2] = x[1] - x[2] + x[3]
    dx[3] = -β * x[2]
end
d = 3

function lorenz(dx, x, σ = 10, β = 8 / 3, ρ = 35)
    dx[1] = σ * (x[2] - x[1])
    dx[2] = x[1] * (ρ - x[3]) - x[2]
    dx[3] = x[1] * x[2] - β * x[3]
end

# Define clusters
cls = [1:3, 4:6]
ks = length.(cls)
n = sum(ks)
l = length(ks)

# Define passivity index.
f = lorenz
Δ = diagm(ones(d)) * 10
P = diagm(ones(d))

# Define topology of the network
μ(A) = 1 / 2 * maximum(size(A)) * maximum(abs.(A))
Φ = BlockArray([
        -1.0 1 0 0 0 0
        1 -2 1 0 -1 1
        0 1 -1 0 1 -1
        0 0 0 -1 1 0
        0 -1 1 1 -2 1
        0 1 -1 0 1 -1
    ], [3, 3], [3, 3])

μs = zeros(l * l - l)
k = 1
for i = 1:l, j = 1:l
    global k
    if j != i
        μs[k] = μ(getblock(Φ, i, j))
        k += 1
    end
end
μm = maximum(μs)
δm = maximum(diag(Δ))

# Define initial constants
γs = [1, 1] * 200
Γiis = [diagm([zeros(k - 1); γ]) for (k, γ) in zip(ks, γs)]
Φiis = [getblock(Φ, i, i) for i = 1:l]
λs = [-maximum(eigvals(Φii - Γii)) for (Φii, Γii) in zip(Φiis, Γiis)]
βis = (δm + 2 * (l - 1) * μm) ./ λs
β = maximum(βis) * 1  # The condition is not statisfied.
αs = β * γs
E = copy(Φ)
for i = 1:l
    setblock!(E, getblock(E, i, i) * β, i, i)
end
E = [zeros(l, l) zeros(l, n); zeros(n, l) E]
H1 = zeros(n, l)
for (k, (j, α)) in zip(cumsum(ks), enumerate(αs))
    H1[k, j] = -α
end
H2 = zeros(n, n)
for (k, α) in zip(cumsum(ks), αs)
    H2[k, k] = α
end
H = [zeros(l, l + n); [H1 H2]]

# Define coupled network dynamics
function statefunc(dx, x, u, t, f = f, E = E, H = H, P = P)
    n = size(E, 1)
    d = size(P, 1)
    # Individual evolution
    for i = 1:n
        f(view(dx, (i-1)*d+1:i*d), view(x, (i-1)*d+1:i*d))
    end
    # Add coupling and control inputs
    dx .+= kron(E - H, P) * x
end

# Simulation network
tspan = (0.0, 100.0)
x0 = rand((n + l) * d) * 1e-3
sol = solve(ODEProblem(statefunc, x0, tspan), saveat = 0.001)
t = sol.t
x = sol.u

# Plot the results
ni = 1
nf = 50
graph = SimpleDiGraph(E - H)
gplot(graph, nodelabel = 1:nv(graph), layout = circular_layout) |> display
subplot(5, 2, 1)
plot(t, getindex.(x, 1))
subplot(5, 2, 2)
plot(getindex.(x, 1), getindex.(x, 2))
subplot(5, 2, 3)
plot(t, abs.(getwave(x, 1, d) - getwave(x, 3, d)))
subplot(5, 2, 4)
plot(t, abs.(getwave(x, 1, d) - getwave(x, 4, d)))
subplot(5, 2, 5)
plot(t, abs.(getwave(x, 1, d) - getwave(x, 5, d)))
subplot(5, 2, 6)
plot(t, abs.(getwave(x, 2, d) - getwave(x, 7, d)))
subplot(5, 2, 7)
plot(t, abs.(getwave(x, 2, d) - getwave(x, 7, d)))
subplot(5, 2, 8)
plot(t, abs.(getwave(x, 2, d) - getwave(x, 8, d)))
subplot(5, 2, 9)
plot(t, abs.(getwave(x, 1, d) - getwave(x, 2, d)))
subplot(5, 2, 10)
plot(t, abs.(getwave(x, 4, d) - getwave(x, 7, d)))
tight_layout()


text:
using Test

import DivPopGa.ClusteredGa as CLGA

@testset "Chromosome constructor with genes" begin
    genes = zeros(Float64, 10)
    c = CLGA.Chromosome(genes)

    @test length(c.genes) == 10
    @test c.clusterid == -1
    @test isinf(c.cost)
end

@testset "Create Chromosome" begin
    lower = [0.0, 1.0, 2.0, 3.0, 5.0]
    upper = [1.0, 2.0, 3.0, 4.0, 6.0]
    c = CLGA.Chromosome(lower, upper)
    for i in 1:length(upper)
        @test c.genes[i] <= upper[i]
        @test c.genes[i] >= lower[i]
    end
end

@testset "Distance of chromosomes" begin
    ch1 = CLGA.Chromosome([0.0, 0.0, 0.0], Inf64, -1)
    ch2 = CLGA.Chromosome([1.0, 1.0, 1.0], Inf64, -1)
    d = CLGA.distance(ch1, ch2)

    @test d isa Float64
    @test d == 3.0
end

@testset "Uniform crossover" begin
    ch1 = CLGA.Chromosome([1.0, 2.0, 3.0])
    ch2 = CLGA.Chromosome([10.0, 20.0, 30.0])

    crossfn = CLGA.makeuniformcrossover()

    child1 = crossfn(ch1, ch2)

    @test child1.genes[1] in [1.0, 10.0]
    @test child1.genes[2] in [2.0, 20.0]
    @test child1.genes[3] in [3.0, 30.0]
    @test isinf(child1.cost)
    @test child1.clusterid == -1
end

@testset "Weighted Crossover" begin
    ch1 = CLGA.Chromosome(
        [1.0, 2.0],           # genes
        Inf64,                # cost
        -1                    # cluster id
    )

    ch2 = CLGA.Chromosome(
        [10.0, 20.0],
        Inf64,
        -1
    )

    crossfn = CLGA.makeweightedcrossover()

    child1 = crossfn(ch1, ch2)
    child2 = crossfn(ch2, ch1)

    @test child1.genes[1] <= 10.0
    @test child1.genes[2] <= 20.0
    @test child1.genes[1] >= 1.0
    @test child1.genes[2] >= 2.0

    @test child2.genes[1] <= 10.0
    @test child2.genes[2] <= 20.0
    @test child2.genes[1] >= 1.0
    @test child2.genes[2] >= 2.0
end

@testset "Linear Crossover" begin
    ch1 = CLGA.Chromosome(
        [1.0, 2.0],           # genes
        Inf64,                # cost
        -1                    # cluster id
    )

    ch2 = CLGA.Chromosome(
        [1.0, 2.0],
        Inf64,
        -1
    )

    function costfn(x::Array{Float64, 1})::Float64
        return abs(x[1] - 3.141592) + abs(x[2] - 2.71828)
    end

    xfn = CLGA.makelinearcrossover(costfn)

    child = xfn(ch1, ch2)
    g = child.genes

    alternative1 = 0.5 * ch1.genes .+ 0.5 * ch2.genes
    alternative2 = 1.5 * ch1.genes .- 0.5 * ch2.genes
    alternative3 = 1.5 * ch2.genes .- 0.5 * ch1.genes

    @test g == alternative1 || g == alternative2 || g == alternative3
end

@testset "Normal mutation" begin
    mutfn = CLGA.makenormalmutation(
        1.0,          # std deviation
        1.0           # mutation probability
    )

    ch = CLGA.Chromosome(
        [1.0, 1.0], # genes
        Inf64,      # cost of chromosome
        -1         # cluster id
    )

    mutated = mutfn(ch)

    @test mutated isa CLGA.Chromosome
    @test mutated.genes[1] < 1 + 5
    @test mutated.genes[1] > 1 - 5
    @test mutated.genes[2] < 1 + 5
    @test mutated.genes[2] > 1 - 5
end

@testset "Random mutation (within range)" begin
    lower = [0.0, 0.0, 0.0, 0.0]
    upper = [9.0, 9.0, 9.0, 9.0]

    ch = CLGA.Chromosome(lower, upper)

    mutatefn = CLGA.makerandommutation(lower, upper, 0.50)

    newch = mutatefn(ch)

    @test newch isa CLGA.Chromosome
    @test all(x -> x == 1,  newch.genes .<= upper)
    @test all(x -> x == 1,  newch.genes .>= lower)
end

@testset "Tournament selection" begin
    pop = [
        CLGA.Chromosome([0.0, 0.0], 1, -1),
        CLGA.Chromosome([0.0, 0.0], 2, -1),
        CLGA.Chromosome([0.0, 0.0], 3, -1),
        CLGA.Chromosome([0.0, 0.0], 4, -1),
        CLGA.Chromosome([0.0, 0.0], 5, -1),
    ]
    chs = CLGA.tournamentselection(pop, 2)

    @test length(chs) == 2
    @test chs[1].cost < 5.0
    @test chs[2].cost < 5.0
end

@testset "KMeans Tournament selection" begin
    pop = [
        CLGA.Chromosome([0.0, 0.0], 1, 1),
        CLGA.Chromosome([0.0, 0.0], 2, 1),
        CLGA.Chromosome([0.0, 0.0], 3, 1),
        CLGA.Chromosome([0.0, 0.0], 4, 1),
        CLGA.Chromosome([0.0, 0.0], 5, 2),
    ]
    chs = CLGA.kmeanstournamentselection(pop, 2)

    @test length(chs) == 2
    @test chs[1].clusterid != chs[2].clusterid
end

@testset "Simulated KMeans Tournament selection" begin
    pop = [
        CLGA.Chromosome([0.0, 0.0], 1, -1),
        CLGA.Chromosome([0.0, 0.0], 2, -1),
        CLGA.Chromosome([1.0, 1.0], 3, -1),
    ]
    chs = CLGA.simulatedkmeanstournamentselection(pop, 3)

    @test length(chs) == 2
    @test CLGA.distance(chs[1], chs[2]) > 0.0
end

@testset "Pi and E with classical FPGA" begin
    function costfn(vals)
        return abs(vals[1] - 3.141592) + abs(vals[2] - 2.71828)
    end
    crossfn = CLGA.makelinearcrossover(costfn)
    lower = [-10.0, -10.0]
    upper = [10.0, 10.0]
    mutatefn = CLGA.makerandommutation(lower, upper, 0.1)
    result = CLGA.ga(
        100, # popsize
        10000, # generations
        lower, #lower
        upper,   #upper
         costfn, #costfunction
         crossfn, #crossoverfunction
         mutatefn, #mutation function
         CLGA.GA_TYPE_CLASSIC, # classical selection
         verbose = true
    )
    @test result isa CLGA.GAResult

    best = result.best
    @test best isa CLGA.Chromosome
    @test best.genes[1] > 2.0
    @test best.genes[1] < 5.0
    @test best.genes[2] > 1.0
    @test best.genes[2] < 4.0

    # Fine tuning
    finetuning_mutatefn = CLGA.makenormalmutation(1.0,  #stddev
                                                 0.05   #mutation prob
    )
    pop = result.population
    for i in 1:500
        pop = CLGA.generation(pop, costfn, crossfn, finetuning_mutatefn, CLGA.GA_TYPE_CLASSIC)
    end
    CLGA.calculatefitness(pop, costfn)
    sort!(pop, by = ch -> ch.cost)
    best = pop[1]
    @test best isa CLGA.Chromosome
    @test best.genes[1] > 3.0  # 3.141592
    @test best.genes[1] < 3.3
    @test best.genes[2] > 2.6  # 2.71828
    @test best.genes[2] < 2.9
end

@testset "Pi and E with simulated kmeans" begin
    function costfn(vals)
        return abs(vals[1] - 3.141592) + abs(vals[2] - 2.71828)
    end
    crossfn = CLGA.makelinearcrossover(costfn)
    lower = [-10.0, -10.0]
    upper = [10.0, 10.0]
    mutatefn = CLGA.makerandommutation(lower, upper, 0.1)
    result = CLGA.ga(
        100, # popsize
        10000, # generations
        lower, #lower
        upper,   #upper
         costfn, #costfunction
         crossfn, #crossoverfunction
         mutatefn, #mutation function
         CLGA.GA_TYPE_CLUSTER_SIM # simulated kmeans selection
    )

    best = result.best
    @test best isa CLGA.Chromosome
    @test best.genes[1] > 2.0
    @test best.genes[1] < 5.0
    @test best.genes[2] > 1.0
    @test best.genes[2] < 4.0

    # Fine tuning
    finetuning_mutatefn = CLGA.makenormalmutation(1.0,  #stddev
                                                 0.05   #mutation prob
    )
    pop = result.population
    for i in 1:500
        pop = CLGA.generation(pop, costfn, crossfn, finetuning_mutatefn, CLGA.GA_TYPE_CLUSTER_SIM)
    end
    CLGA.calculatefitness(pop, costfn)
    sort!(pop, by = ch -> ch.cost)
    best = pop[1]
    @test best isa CLGA.Chromosome
    @test best.genes[1] > 3.0  # 3.141592
    @test best.genes[1] < 3.3
    @test best.genes[2] > 2.6  # 2.71828
    @test best.genes[2] < 2.9
end

@testset "Pi and E with Kmeans clustering GA" begin
    function costfn(vals)
        return abs(vals[1] - 3.141592) + abs(vals[2] - 2.71828)
    end
    crossfn = CLGA.makelinearcrossover(costfn)
    mutatefn = CLGA.makenormalmutation(10.0,  #stddev
                                      0.10   #mutation prob
    )
    result = CLGA.ga(
        100, # popsize
        500, # generations
        [-10.0, -10.0], #lower
        [10.0, 10.0],   #upper
         costfn, #costfunction
         crossfn, #crossoverfunction
         mutatefn, #mutation function
         CLGA.GA_TYPE_CLUSTER # cluster based selection
    )

    best = result.best
    @test best isa CLGA.Chromosome
    @test best.genes[1] > 2.0
    @test best.genes[1] < 5.0
    @test best.genes[2] > 1.0
    @test best.genes[2] < 4.0

    # Fine tuning
    finetuning_mutatefn = CLGA.makenormalmutation(1.0,  #stddev
                                                 0.05   #mutation prob
    )
    pop = result.population
    for i in 1:500
        pop = CLGA.generation(pop, costfn, crossfn, finetuning_mutatefn, CLGA.GA_TYPE_CLUSTER)
    end
    CLGA.calculatefitness(pop, costfn)
    sort!(pop, by = ch -> ch.cost)
    best = pop[1]
    @test best isa CLGA.Chromosome
    @test best.genes[1] > 3.0  # 3.141592
    @test best.genes[1] < 3.3
    @test best.genes[2] > 2.6  # 2.71828
    @test best.genes[2] < 2.9
end

@testset "Pi and E with Kmeans clustering GA - with elitism" begin
    function costfn(vals)
        return abs(vals[1] - 3.141592) + abs(vals[2] - 2.71828)
    end
    crossfn = CLGA.makelinearcrossover(costfn)
    mutatefn = CLGA.makenormalmutation(10.0,  #stddev
                                      0.10   #mutation prob
    )
    result = CLGA.ga(
        100, # popsize
        500, # generations
        [-10.0, -10.0], #lower
        [10.0, 10.0],   #upper
         costfn, #costfunction
         crossfn, #crossoverfunction
         mutatefn, #mutation function
         CLGA.GA_TYPE_CLUSTER, # cluster based selection
         elitism = 1
    )

    best = result.best

    @test best isa CLGA.Chromosome
    @test best.genes[1] > 2.0
    @test best.genes[1] < 5.0
    @test best.genes[2] > 1.0
    @test best.genes[2] < 4.0
end

@testset "(Hybrid) Pi and E with Kmeans clustering GA - with elitism - then classical ga with initial pop" begin
    function costfn(vals)
        return abs(vals[1] - 3.141592) + abs(vals[2] - 2.71828)
    end
    crossfn = CLGA.makelinearcrossover(costfn)
    mutatefn = CLGA.makenormalmutation(10.0,  #stddev
                                      0.10   #mutation prob
    )
    result = CLGA.ga(
        100, # popsize
        10000, # generations
        [-10.0, -10.0], #lower
        [10.0, 10.0],   #upper
         costfn, #costfunction
         crossfn, #crossoverfunction
         mutatefn, #mutation function
         CLGA.GA_TYPE_CLUSTER, # cluster based selection
         elitism = 1
    )

    result_classical = CLGA.ga(
        100, # popsize
        10000, # generations
        [-10.0, -10.0], #lower
        [10.0, 10.0],   #upper
         costfn, #costfunction
         crossfn, #crossoverfunction
         mutatefn, #mutation function
         CLGA.GA_TYPE_CLASSIC, # classical FPGA
         elitism = 1,
         initialpopulation = result.population
    )
    best = result_classical.best

    @test best isa CLGA.Chromosome
    @test best.genes[1] > 2.0
    @test best.genes[1] < 5.0
    @test best.genes[2] > 1.0
    @test best.genes[2] < 4.0
end

@testset "Pi and E with classical FPGA and BLX-alpha crossover" begin
    function costfn(vals)
        return abs(vals[1] - 3.141592) + abs(vals[2] - 2.71828)
    end
    crossfn = CLGA.makeblxalphacrossover()
    mutatefn = CLGA.makenormalmutation(10.0,  #stddev
                                      0.10   #mutation prob
    )
    result = CLGA.ga(
        100, # popsize
        10000, # generations
        [-10.0, -10.0], #lower
        [10.0, 10.0],   #upper
         costfn, #costfunction
         crossfn, #crossoverfunction
         mutatefn, #mutation function
         CLGA.GA_TYPE_CLASSIC # classical selection
    )

    best = result.best
    @test best isa CLGA.Chromosome
    @test best.genes[1] > 2.0
    @test best.genes[1] < 5.0
    @test best.genes[2] > 1.0
    @test best.genes[2] < 4.0

end

@testset "hybridga() - pi and e" begin
    function costfn(vals)
        return abs(vals[1] - 3.141592) + abs(vals[2] - 2.71828)
    end
    crossfn = CLGA.makelinearcrossover(costfn)
    mutatefn = CLGA.makenormalmutation(10.0,  #stddev
                                      0.10   #mutation prob
    )
    result = CLGA.hybridga(
        100, # popsize
        [500, 500], # generations (500 for kmeans, 500 for classical)
        [-10.0, -10.0], #lower
        [10.0, 10.0],   #upper
         costfn, #costfunction
         crossfn, #crossoverfunction
         mutatefn, #mutation function
         elitism = 1
    )

    best = result[1]
    @test best isa CLGA.Chromosome
    @test best.genes[1] > 2.0
    @test best.genes[1] < 5.0
    @test best.genes[2] > 1.0
    @test best.genes[2] < 4.0
    @info best
end

max_stars_repo_name: jbytecode/divpopga

max_issues_repo_name: jbytecode/divpopga

max_forks_repo_name: jbytecode/divpopga


text:
# ---
# title: 486. Predict the Winner
# id: problem486
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming, Minimax
# link: <https://leetcode.com/problems/predict-the-winner/description/>
# hidden: true
# ---
#
# Given an array of scores that are non-negative integers. Player 1 picks one of
# the numbers from either end of the array followed by the player 2 and then
# player 1 and so on. Each time a player picks a number, that number will not be
# available for the next player. This continues until all the scores have been
# chosen. The player with the maximum score wins.
#
# Given an array of scores, predict whether player 1 is the winner. You can
# assume each player plays to maximize his score.
#
# **Example 1:**
#
#
#
#     Input: [1, 5, 2]
#     Output: False
#     Explanation: Initially, player 1 can choose between 1 and 2.
#     If he chooses 2 (or 1), then player 2 can choose from 1 (or 2) and 5. If player 2 chooses 5, then player 1 will be left with 1 (or 2).
#     So, final score of player 1 is 1 + 2 = 3, and player 2 is 5.
#     Hence, player 1 will never be the winner and you need to return False.
#
#
#
#
# **Example 2:**
#
#
#
#     Input: [1, 5, 233, 7]
#     Output: True
#     Explanation: Player 1 first chooses 1. Then player 2 have to choose between 5 and 7. No matter which number player 2 choose, player 1 can choose 233.
#     Finally, player 1 has more score (234) than player 2 (12), so you need to return True representing player1 can win.
#
#
#
#
# **Constraints:**
#
#   * 1 <= length of the array <= 20.
#   * Any scores in the given array are non-negative integers and will not exceed 10,000,000.
#   * If the scores of both players are equal, then player 1 is still the winner.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
## Copyright (c) 2015, Júlio Hoffimann Mendes <juliohm@stanford.edu>
##
## Permission to use, copy, modify, and/or distribute this software for any
## purpose with or without fee is hereby granted, provided that the above
## copyright notice and this permission notice appear in all copies.
##
## THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
## WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
## MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
## ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
## WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
## ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
## OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

# Julia doesn't support optional dependencies yet. We have
# to disable precompilation in order to avoid issues with
# missing components at runtime.
__precompile__(false)

module ImageQuilting

using ImageFiltering
using ImageMorphology
using LightGraphs
using LightGraphsFlows
using RecipesBase
using Base: @nexprs, @nloops, @nref
using Combinatorics: nthperm!
using StatsBase: sample, weights
using Primes: factor
using ProgressMeter: Progress, next!
using Hwloc: num_physical_cores

# optional dependencies
try
  using OpenCL
catch
  global cl = nothing
end

try
  using CLFFT
  global clfft = CLFFT
catch
  global clfft = nothing
end

# GeoStats.jl interface
importall GeoStatsBase
using GeoStatsDevTools

include("utils.jl")
include("utils_gpu.jl")
include("datatypes.jl")
include("plot_recipes.jl")
include("imfilter_cpu.jl")
include("imfilter_gpu.jl")
include("relaxation.jl")
include("tau_model.jl")
include("dijkstra_cut.jl")
include("boykov_kolmogorov_cut.jl")
include("simplex_transform.jl")
include("voxel_reuse.jl")
include("iqsim.jl")

include("geostats_api.jl")

export
  # functions
  iqsim,
  voxelreuse,

  # data types
  HardData,

  # geostats solver
  ImgQuilt,

  # deprecated
  SoftData

end


text:
#!/usr/bin/env julia

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "-M"
    help = "Print out the list of files this program will create."
    action = :store_true
    "config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskTracer.model
using DiskTracer.constants
using DiskTracer.geometry
using SpecialFunctions

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]
pars = convert_dict(config["parameters"], config["model"])

mol = mols[species*transition]

nr=128
rmax = 500 * AU
nz=128
zmax = 200 * AU
# Get the interpolator object for this particular setup.
interp = DiskTracer.model.get_grids(pars, mol, mol.nu_0, nr, rmax, nz, zmax)

using LaTeXStrings
using Plots
pyplot()

# pick a pixel (xprime, yprime), and trace a full ray through the disk

function plot_ray(xprime, yprime)

    v0 = 1.0 # km/s

    # Get a decent estimate of DeltaVmax
    # Calculate DeltaVmax for this ray
    rcyl_min = abs(xprime)
    DeltaVmax = sqrt(2 * kB * temperature(rcyl_min, pars)/mol.mol_weight + (pars.ksi * 1e5)^2) * 1e-5 # km/s

    # DeltaVmax = 0.5 # km/s
    rmax = 600 * AU
    nzp = 512
    zps = linspace(-600 * AU, 600 * AU, nzp)

    rcyls = Array{Float64}(nzp)
    zs = Array{Float64}(nzp)
    Ss = Array{Float64}(nzp)
    Us = Array{Float64}(nzp)
    Deltas = Array{Float64}(nzp)
    alphas = Array{Float64}(nzp)

    z0s = DiskTracer.geometry.get_zps(xprime, yprime, pars, v0)
    println("z0s ", z0s)

    for i=1:nzp
        zprime = zps[i]
        rcyl = DiskTracer.geometry.get_r_cyl(xprime, yprime, zprime, pars.incl)
        z = DiskTracer.geometry.get_z(xprime, yprime, zprime, pars.incl)

        rcyls[i] = rcyl
        zs[i] = z

        Deltav = v0 * 1e5 - DiskTracer.geometry.get_vlos(xprime, yprime, zprime, pars)
        DeltaV = sqrt(2 * kB * DiskTracer.model.temperature(rcyl, pars)/mol.mol_weight + (pars.ksi * 1e5)^2)

        Deltas[i] = exp(-Deltav^2/DeltaV^2)

        Us[i] = DiskTracer.model.Upsilon_nu(rcyl, z, mol.nu_0, pars, mol)
        alphas[i] = Us[i] * exp(-Deltav^2/DeltaV^2)

    end

    # get bounding zps
    zbounds = DiskTracer.geometry.get_bounding_zps(xprime, yprime, pars, v0, DeltaVmax, rmax)
    if length(zbounds) == 2
        zstart, zend = zbounds
        println("H bound verified merged: ", DiskTracer.geometry.verify_h_bound(xprime, yprime, zstart, zend, pars, 0.5))
    else
        z1start, z1end, z2start, z2end = zbounds
        println("H bound verified front-side: ", DiskTracer.geometry.verify_h_bound(xprime, yprime, z1start, z1end, pars, 0.5))
        println("H bound verified rear-side: ", DiskTracer.geometry.verify_h_bound(xprime, yprime, z2start, z2end, pars, 0.5))
    end

    # Call the routine that actually traces the ray, and returns zp, tau, and I quantities
    zpts, sol = DiskTracer.trace.trace_pixel(xprime, yprime, v0, mol, pars, interp)
    println("end tau, I: ", sol.u[end])
    println("fancy steps: ", length(zpts))

    # Also call our custom integrator to get tau
    zpts_c, tau_c, tot_i = DiskTracer.trace_simple.trace_pixel(xprime, yprime, v0, mol, pars, interp)
    println("num iterations: ", length(zpts_c))
    println("tot intensity: ", tot_i)
    # println("zpts_c: ", zpts_c)
    # println("tau_c: ", tau_c)

    p1 = plot(zps./AU, Us, xlabel=L"$z^\prime$", ylabel=L"$\Upsilon$") #, yaxis=(scale=:log10))

    p2 = plot(zps./AU, alphas, xlabel=L"$z^\prime$", ylabel=L"$\alpha$")
    # vline!((zbounds./AU)', line=([:red :green :blue :cyan]))

    p3 = plot(zps./AU, alphas, xlabel=L"$z^\prime$", ylabel=L"$\alpha$", yaxis=(scale=:log10))
    vline!((zbounds./AU)', line=([:red :green :blue :cyan]))

    p4 = plot(zpts_c./AU, tau_c, xlim=(200, 400), ylabel=L"$\tau$")
    plot!(zpts./AU, sol[1,:], xlim=(200, 400), ylabel=L"$\tau$")

    # p6 = plot(zpts./AU, sol[2,:], xlim=(-600, 600), ylabel=L"$I$")
    # p6 = plot(zpts./AU, zpts_inv./AU, xlim=(-600, 600))

    # println("final ", sol[2,end])

    p = plot(p1,p2,p3,p4, layout=(4,1), legend=false, size=(800,1200))

    savefig(p, "ray_trace_simple.png")

end

plot_ray(100.*AU, -100.0*AU)


text:
#+ term=true

y= [2, 5, 12]

#+ term=true

x = 1:10
d = Dict("Weave" => "testing")
y = [2, 4 ,8]

#+

x = [12, 10]
println(y)
println(x)

#+ echo=false

println("Results without code")
println(x)

#+ eval=false; tangle=false

y = randn(5)
println("Don't eval, but show code")

#+

y = 1:5
println(y)

#+ results="hidden"

a = "Don't print me"
println(a)

#+ results="as_is"

println("No markup for results.")

#'
#'
#' Test wrapping:
#'
#+

println(collect(0:10:1000))

#+ wrap=false

println(collect(0:10:1000))

#+ line_width=60

println(collect(0:10:1000))


text:
"""
Takes a time series and produces a matrix containing all the possible lags.

Those points which do not have lags of that magnitude are set to Null.
"""
function lagmatrix{T}(x::AbstractVector{T})
    n = length(x)
    X = fill(Nullable{T}(),n,n)
    for j in 1:n, i in 1:n
        X[i,j] = j-i>=0 ? x[j-i+1] : Nullable()
        end
    X
end

"""
Create a lagmatrix for any given Quantity
"""
function makelagmatrix(q::Quantity)
    lagmatrix(quantity(q))
end

function preparedata(p::Stage,d::Discharge,range,M)
    H = makelagmatrix(p)
    Q = quantity(d)
    Hm = validate(H,range,M)
    Qm = validate(Q,range,M)
    Hm,Qm
end

"""
Given a matrix formed by `lagmatrix`, a range of valid values desired and a number of lags,
return those data points which are valid data
"""
function validate{T}(X::Matrix{Nullable{T}},range,n::Int)
    X = X[:,range]
    i = 1
    while any(map(isnull,X[1:n,i]))
        i+=1
    end
    _,m = size(X)
    Y = zeros(T,n,m-i+1)
    for k in 1:n, j in i:m
        Y[k,j-i+1] = get(X[k,j])
    end
    Y
end

"""
The same idea, but give all the valid response values from a vector of discharges

Try not to pass in a range which is invalid for the vector
"""
function validate(x::AbstractVector,range,n::Int)
    if range.start < n
        return x[n:range.stop]
    else
        return x[range]
    end
end

include("statsutils.jl")


text:
###=============================================================================
#
#          FILE: Infra.jl
#
#         USAGE: include("../Infra.jl")
#
#   DESCRIPTION: Julia interface to launch containers through Azure VM
#
#       OPTIONS: ---
#  DEPENDENCIES: sshpass
#          BUGS: ---
#         NOTES: ---
#        AUTHOR: Raphael P. Ribeiro <raphaelpr01@gmail.com>
#  ORGANIZATION: GSD-UFAL
#       CREATED: 23-09-2015 17:16
#
###=============================================================================

module Infra

using Requests: get, parse
using Docker

###== Top-level variables ======================================================

const nprocs=16 # number of cpu cores available from server
const ssh_key=homedir()*"/.ssh/azkey"
const ssh_pubkey=homedir()*"/.ssh/azkey.pub"
const carray_dir= length(LOAD_PATH) == 3 ? LOAD_PATH[3]*"Infra/src" : Pkg.dir("Infra")*"/src" # Infra package directory

###=============================================================================

type GlobalSettings
    host::ASCIIString
    url::ASCIIString
    passwd::ASCIIString
    cpuset::Array{Int64,1} # number of cpu cores available from server
end

type Container # Abstraction for Docker container
          cid::AbstractString
          pid::Integer
          cpuset::Array{Int64,1}
          mem_size::Integer

          function Container(cid,pid,cpuset=[],mem_size=512)
              return new(cid,pid,cpuset,mem_size)
          end
end

const map_containers = Dict{Integer, Container}()

GlobalSettings() = GlobalSettings("","","",collect(0:nprocs-1))
const SETTINGS = GlobalSettings()

set_h(h) = SETTINGS.host=h
set_url(u) = SETTINGS.url=u
set_passwd(p) = SETTINGS.passwd=p
set_cpuset(c) = SETTINGS.cpuset=c

let next_key = 1    # cid -> container id
    global get_next_key
    function get_next_key() retval = next_key
        next_key += 1
        retval
    end
end

function parse_cpuset(n_of_cpus)
    if (n_of_cpus == 0 || length(SETTINGS.cpuset) < n_of_cpus)
        "",[]
    else
        s = strip(string(SETTINGS.cpuset[1:n_of_cpus]),['[',']'])
        a = SETTINGS.cpuset[1:n_of_cpus]
        deleteat!(SETTINGS.cpuset,1:n_of_cpus)
        s,a
    end
end

@doc """
### set_host(h::AbstractString,p::AbstractString)

Configures passwordless SSH connections at host `h` whose password is `p`.

This function calls the `cloud_setup.sh` script which requires `sshpass`.

```Example
set_host("cloudarray.cloudapp.net","password")
```
""" ->
function set_host(h::AbstractString,p::AbstractString)
    reply = success(`$(carray_dir)/cloud_setup.sh $h $p`) # set up ssh. if errors occurs, return false
    if (reply)
        set_h(h)
        set_url(h*":4243")
        set_passwd(p)
        true
    else
        error("There is an error during SSH configuration. Please see the log for more details: cloud_setup.log")
        false
    end
end

function get_port()
    response = get("http://$(SETTINGS.host):8000")
    parse(Int,join(map(Char,response.data)))
end

@doc """
### create_containers(n_of_containers::Integer, n_of_cpus::Integer, mem_size::Integer; tunnel::bool)

Launches Docker containers and adds them as Julia workers configured with passwordless SSH.

This function requires `sshpass` to be installed:

* Debian-based Linux distros as Ubuntu:

```
sudo apt-get install sshpass
```

* OS X through [macports](http://macports.org):

```
sudo port install sshpass
```

```Example
create_containers(2,3,1024) # 2 containers with 3 CPU Cores and 1gb RAM
create_containers(1,2,512)  # 1 container with 2 CPU Cores and 512mb RAM
```
""" ->
function create_containers(n_of_containers::Integer, n_of_cpus=0, mem_size=512;tunnel=false)
        reserved_mem=200 # reserved memory for initializing a worker into a container
        mem_size=mem_size+reserved_mem
        for i in 1:n_of_containers
            ssh_config = false
            key = get_next_key()
            port = 3000+get_port()
            n_of_cpus=parse_cpuset(n_of_cpus)
            # Creating a docker container at VM
            info("Creating container ($key)...")
            container = Docker.create_container("$(SETTINGS.url)","cloudarray:latest",memory=mem_size*(10^6),cpuSets=n_of_cpus[1],portBindings=[22,"$port"])
            Docker.start_container("$(SETTINGS.url)",container["Id"])
            info("Creating container ($key)... OK")
            # Configuring ssh without password (transfer public key to container)
            info("SSH configuration ($key)... ")
            while !ssh_config
                ssh_config = success(pipeline(`cat $ssh_pubkey`,`sshpass -p $(SETTINGS.passwd) ssh -o StrictHostKeyChecking=no -p $port root@$(SETTINGS.host) 'umask 077; mkdir -p ~/.ssh; cat >> ~/.ssh/authorized_keys'`)) # if ssh configuration is successful: return true or false
                if !ssh_config
                    info("SSH configuration ($key) failed! Trying again...")
                end
            end
            info("SSH configuration ($key)... OK")
            info("Adding worker ($key)...")
            pid = addprocs(["root@$(SETTINGS.host)"];tunnel=tunnel,sshflags=`-i $ssh_key -p $port`,dir="/opt/julia/bin",exename="/opt/julia/bin/julia")
            info("Adding worker ($key)... OK")
            map_containers[key] = Container(chomp(container["Id"]),pid[1],n_of_cpus[2],mem_size) # Adding Container to Dict
        end
end

@doc """
### delete_containers(args...)

Removes the specified container(s)/worker(s).

```Example
delete_containers(3)    # delete container 3
create_containers(1:5)  # delete from 1st to 5th container
create_containers(all)  # delete all containers
```
""" ->
function delete_containers(args...) #  (splat) variable number of arguments. Ex.: delete_containers(1,2,3) or delete_containers(1:3)
    containers_rmlist = Dict()
    if vcat(args...)[1] == all
        for i in collect(keys(map_containers))
            if haskey(map_containers, i) # container exist?
                container = map_containers[i]
                rmprocs(container.pid)
                delete!(map_containers,i)
                set_cpuset(collect(0:nprocs-1))
                Docker.remove_container("$(SETTINGS.url)","$(container.cid)")
            end
        end
    else
        for i in vcat(args...) # vcat -> concatenate to a array 1 dimension
            if haskey(map_containers, i) # container exist?
                container = map_containers[i]
                rmprocs(container.pid)
                delete!(map_containers,i)
                set_cpuset([container.cpuset;SETTINGS.cpuset])
                Docker.remove_container("$(SETTINGS.url)","$(container.cid)")
            end
        end
    end
end

@doc """
### containers()

Returns the list of all containers' processes identifiers (IDs).

```Example
containers()
```
""" ->
function containers()
    sort(collect(keys(map_containers)))
end

@doc """
### ncontainers()

Gets the number of available container processes.

```Example
ncontainers()
```
""" ->
function ncontainers()
    length(map_containers)
end

@doc """
### list_containers()

List container(s) as a sorted list.

```Example
list_containers()
```
""" ->
function list_containers()
    for key in sort(collect(keys(map_containers)))
           println("$key => $(map_containers[key])")
    end
end

@doc """
### mem_usage(key::Integer)

Returns the container memory usage.

```Example
mem_usage(number_of_container)
```
""" ->
function mem_usage(key::Integer)
    Docker.stats_container("$(SETTINGS.url)","$(map_containers[key].cid)")["memory_stats"]["usage"]/10^6
end

@doc """
### cpu_usage(key::Integer)

Returns the container CPU usage (%).

```Example
cpu_usage(number_of_container)
```
""" ->
function cpu_usage(key)
    stats = Docker.stats_container("$(SETTINGS.url)","$(map_containers[key].cid)")

    percpu_usage = stats["cpu_stats"]["cpu_usage"]["percpu_usage"]
    previousSystem = stats["precpu_stats"]["system_cpu_usage"]
    previousCPU = stats["precpu_stats"]["cpu_usage"]["total_usage"]
    totalUsage = stats["cpu_stats"]["cpu_usage"]["total_usage"]
    systemUsage = stats["cpu_stats"]["system_cpu_usage"]

    cpuPercent = 0.0
    cpuDelta = totalUsage - previousCPU
    systemDelta = systemUsage - previousSystem
    cpuPercent = (cpuDelta / systemDelta) * length(percpu_usage) * 100.0
    cpuPercent
end

@doc """
### io_usage(key::Integer)

Returns the number of kilobytes read and written by the cgroup.

```Example
io_usage(number_of_container)
```
""" ->
function io_usage(key::Integer)
    stats = Docker.stats_container("$(SETTINGS.url)","$(map_containers[key].cid)")["blkio_stats"]
    w = stats["io_service_bytes_recursive"][1]["value"]/10^3   # write
    r = stats["io_service_bytes_recursive"][2]["value"]/10^3   # read
    [w,r]
end

@doc """
### net_usage(key::Integer)

Returns networking TX/RX usage.

tx = number of bytes transmitted
rx = number of bytes reiceved

```Example
net_usage(number_of_container)
```
""" ->
function net_usage(key::Integer)
    stats = Docker.stats_container("$(SETTINGS.url)","$(map_containers[key].cid)")["networks"]["eth0"]
    tx = stats["tx_bytes"]
    rx = stats["rx_bytes"]
    [tx,rx]
end

end


text:
using Iris
using Test

@testset "PumpFunctions" begin
        F = @test_nowarn Iris.Common.PumpFunctions.PiecewiseConstant()
        @test F(Point(0,0)) ≈ 1
        F = @test_nowarn Iris.Common.PumpFunctions.PiecewiseConstant(.1)
        @test F(Point(0,0)) ≈ .1
end


text:
using LightGraphs, SimpleWeightedGraphs, DataStructures

function limited_dijkstra_shortest_paths(g::AbstractGraph,
    srcs::Vector{U},
    distmx::AbstractMatrix{T}=weights(g);
    allpaths=false,
    trackvertices=false,
    limit::T=zero(T)
    ) where T <: Real where U <: Integer

    nvg = nv(g)
    dists = fill(typemax(T), nvg)
    parents = zeros(U, nvg)
    visited = zeros(Bool, nvg)

    pathcounts = zeros(UInt64, nvg)
    preds = fill(Vector{U}(), nvg)
    H = PriorityQueue{U,T}()
    # fill creates only one array.

    dolimit = (limit > zero(T))

    for src in srcs
        dists[src] = zero(T)
        visited[src] = true
        pathcounts[src] = 1
        H[src] = zero(T)
    end

    closest_vertices = Vector{U}()  # Maintains vertices in order of distances from source
    sizehint!(closest_vertices, nvg)

    while !isempty(H)
        u = dequeue!(H)

        if trackvertices
            push!(closest_vertices, u)
        end

        d = dists[u] # Cannot be typemax if `u` is in the queue

        if dolimit
            if d > limit
                break
            end
        end

        for v in outneighbors(g, u)
            alt = d + distmx[u, v]

            if !visited[v]
                visited[v] = true
                dists[v] = alt
                parents[v] = u

                pathcounts[v] += pathcounts[u]
                if allpaths
                    preds[v] = [u;]
                end
                H[v] = alt
            elseif alt < dists[v]
                dists[v] = alt
                parents[v] = u
                #615
                pathcounts[v] = pathcounts[u]
                if allpaths
                    resize!(preds[v], 1)
                    preds[v][1] = u
                end
                H[v] = alt
            elseif alt == dists[v]
                pathcounts[v] += pathcounts[u]
                if allpaths
                    push!(preds[v], u)
                end
            end
        end
    end

    if trackvertices
        for s in vertices(g)
            if !visited[s]
                push!(closest_vertices, s)
            end
        end
    end

    for src in srcs
        pathcounts[src] = 1
        parents[src] = 0
        empty!(preds[src])
    end

    dijkstrastate = LightGraphs.DijkstraState{T,U}(parents, dists, preds, pathcounts, closest_vertices)
    trim_dijkstra_state!(dijkstrastate, limit, skipinf=true)
    return dijkstrastate
end

# single vertice
limited_dijkstra_shortest_paths(g::AbstractGraph, src::Integer, distmx::AbstractMatrix=weights(g); allpaths=false, trackvertices=false, limit=0.) =
limited_dijkstra_shortest_paths(g, [src;], distmx; allpaths=allpaths, trackvertices=trackvertices, limit=limit)

function trim_dijkstra_state(dijkstra_state::LightGraphs.DijkstraState{T,U}, limit::T=zero(T); skipinf=true) where T <: Real where U <: Integer
    s = LightGraphs.DijkstraState(
        deepcopy(dijkstra_state.parents),
        deepcopy(dijkstra_state.dists),
        deepcopy(dijkstra_state.predecessors), # no need for deep copy because we are not modifying content of predecessors itself
        deepcopy(dijkstra_state.pathcounts),
        deepcopy(dijkstra_state.closest_vertices)
    )
    return trim_dijkstra_state!(s, limit; skipinf=true)
end

function trim_dijkstra_state!(dijkstra_state::LightGraphs.DijkstraState{T,U}, limit::T=zero(T); skipinf=true) where T <: Real where U <: Integer
    s = dijkstra_state
    dolimit = (limit > zero(T))

    if (dolimit)
        for i in 1:length(s.dists)
            d = s.dists[i]
            if d > limit && (d < typemax(T) || !skipinf)
                s.parents[i] = zero(U)
                s.dists[i] = typemax(T)
                s.predecessors[i] = Vector{U}()
                s.pathcounts[i] = zero(U)
            end
        end
    end
    return s
end

max_stars_repo_name: kymok/local-betweenness

max_issues_repo_name: kymok/local-betweenness

max_forks_repo_name: kymok/local-betweenness


text:
using .ArchGDAL

const GDAL_X_INDEX = ForwardIndex()
const GDAL_Y_INDEX = ReverseIndex()
const GDAL_BAND_INDEX = ForwardIndex()
const GDAL_X_ARRAY = ForwardArray()
const GDAL_Y_ARRAY = ReverseArray()
const GDAL_BAND_ARRAY = ForwardArray()
const GDAL_RELATION = ForwardRelation()
const GDAL_X_LOCUS = Start()
const GDAL_Y_LOCUS = Start()

export GDALarray, GDALstack

# Array ########################################################################

"""
    GDALarray(filename; kw...)

Load a file lazily using gdal. `GDALarray` will be converted to [`GeoArray`](@ref)
after indexing or other manipulations. `GeoArray(GDALarray(filename))` will do this
immediately.

`GDALarray`s are always 3 dimensional, and have `X`, `Y` and [`Band`](@ref) dimensions.

# Arguments

- `filename`: `String` pointing to a tif or other file that GDAL can load.

# Keywords

- `crs`: crs to use instead of the detected crs
- `mappedcrs`: CRS format like `EPSG(4326)` used in `Selectors` like `Between` and `At`, and
    for plotting. Can be any CRS `GeoFormat` from GeoFormatTypes.jl, like `WellKnownText`.
- `name`: `Symbol` name for the array.
- `dims`: `Tuple` of `Dimension`s for the array. Detected automatically, but can be passed in.
- `refdims`: `Tuple of` position `Dimension`s the array was sliced from.
- `missingval`: Value reprsenting missing values. Detected automatically when possible, but
    can be passed it.
- `metadata`: `Metadata` object for the array. Detected automatically as
    `Metadata{:GDAL}`, but can be passed in.

# Example

```julia
A = GDALarray("folder/file.tif"; mappedcrs=EPSG(4326))
# Select Australia using lat/lon coords, whatever the crs is underneath.
A[Y(Between(-10, -43), X(Between(113, 153)))
```
"""
struct GDALarray{T,N,F,D<:Tuple,R<:Tuple,Na<:Symbol,Me,Mi,S
                } <: DiskGeoArray{T,N,D,LazyArray{T,N}}
    filename::F
    dims::D
    refdims::R
    name::Na
    metadata::Me
    missingval::Mi
    size::S
end
function GDALarray(filename::AbstractString; kw...)
    isfile(filename) || error("file not found: $filename")
    _gdalread(filename) do raster
        GDALarray(raster, filename; kw...)
    end
end
function GDALarray(raster::AG.RasterDataset, filename, key=nothing;
    crs=nothing,
    mappedcrs=nothing,
    dims=dims(raster, crs, mappedcrs),
    refdims=(),
    name=Symbol(""),
    metadata=metadata(raster),
    missingval=missingval(raster)
)
    sze = size(raster)
    T = eltype(raster)
    N = length(sze)
    name = Symbol(name)
    GDALarray{T,N,typeof.((filename,dims,refdims,name,metadata,missingval,sze))...
             }(filename, dims, refdims, name, metadata, missingval, sze)
end

# AbstractGeoArray methods

"""
    Base.write(filename::AbstractString, ::Type{GDALarray}, A::AbstractGeoArray; kw...)

Write a [`GDALarray`](@ref) to file, `.tif` by default, but other GDAL drivers also work.

# Keywords

- `driver::String`: a GDAL driver name. Guessed from the filename extension by default.
- `compress::String`: GeoTIFF compression flag. "DEFLATE" by default.
- `tiled::Bool`: GeoTiff tiling. Defaults to `true`.

Returns `filename`.
"""
function Base.write(
    filename::AbstractString, ::Type{<:GDALarray}, A::AbstractGeoArray{T,2}; kw...
) where T
    all(hasdim(A, (XDim, Y))) || error("Array must have Y and X dims")

    correctedA = permutedims(A, (X(), Y())) |>
        a -> reorder(a, (X(GDAL_X_INDEX), Y(GDAL_Y_INDEX))) |>
        a -> reorder(a, GDAL_RELATION)
    checkarrayorder(correctedA, (GDAL_X_ARRAY, GDAL_Y_ARRAY))
    checkindexorder(correctedA, (GDAL_X_INDEX, GDAL_Y_INDEX))

    nbands = 1
    indices = 1
    _gdalwrite(filename, correctedA, nbands, indices; kw...)
end
function Base.write(
    filename::AbstractString, ::Type{<:GDALarray}, A::AbstractGeoArray{T,3}, kw...
) where T
    all(hasdim(A, (X, Y))) || error("Array must have Y and X dims")
    hasdim(A, Band()) || error("Must have a `Band` dimension to write a 3-dimensional array")

    correctedA = permutedims(A, (X(), Y(), Band())) |>
        a -> reorder(a, (X(GDAL_X_INDEX), Y(GDAL_Y_INDEX), Band(GDAL_BAND_INDEX))) |>
        a -> reorder(a, GDAL_RELATION)
    checkarrayorder(correctedA, (GDAL_X_ARRAY, GDAL_Y_ARRAY, GDAL_BAND_ARRAY))
    checkindexorder(correctedA, (GDAL_X_INDEX, GDAL_Y_INDEX, GDAL_BAND_INDEX))

    nbands = size(correctedA, Band())
    indices = Cint[1:nbands...]
    _gdalwrite(filename, correctedA, nbands, indices; kw...)
end

# AbstractGeoStack methods

"""
    GDALstack(filenames; keys, kw...)
    GDALstack(filenames...; keys, kw...)
    GDALstack(filenames::NamedTuple; kw...)

Convenience method to create a DiskStack  of [`GDALarray`](@ref) from `filenames`.

Load a stack of files lazily from disk.

# Arguments

- `filenames`: A NamedTuple of stack keys and `String` filenames, or a `Tuple`,
    `Vector` or splatted arguments of `String` filenames.

# Keyword arguments

- `keys`: Used as stack keys when a `Tuple`, `Vector` or splat of filenames are passed in.
- `window`: A `Tuple` of `Dimension`/`Selector`/indices that will be applied to the
    contained arrays when they are accessed.
- `metadata`: a `DimensionalData.Metadata` object.
- `childkwargs`: A `NamedTuple` of keyword arguments to pass to the `childtype` constructor.
- `refdims`: `Tuple` of  position `Dimension` the array was sliced from.

# Example

Create a `GDALstack` from four files, that sets the child arrays `mappedcrs` value
when they are loaded.

```julia
files = (:temp="temp.tif", :pressure="pressure.tif", :relhum="relhum.tif")
stack = GDALstack(files; childkwargs=(mappedcrs=EPSG(4326),))
stack[:relhum][Y(Contains(-37), X(Contains(144))
```
"""
GDALstack(args...; kw...) = DiskStack(args...; childtype=GDALarray, kw...)

withsource(f, ::Type{<:GDALarray}, filename::AbstractString, key...) = _gdalread(f, filename)

# DimensionalData methods for ArchGDAL types ###############################

function DD.dims(raster::AG.RasterDataset, crs=nothing, mappedcrs=nothing)
    gt = try
        AG.getgeotransform(raster) catch GDAL_EMPTY_TRANSFORM end
    xsize, ysize = size(raster)

    nbands = AG.nraster(raster)
    band = Band(1:nbands, mode=Categorical(Ordered()))
    crs = crs isa Nothing ? GeoData.crs(raster) : crs
    xy_metadata = Metadata{:GDAL}()

    # Output Sampled index dims when the transformation is lat/lon alligned,
    # otherwise use Transformed index, with an affine map.
    if _isalligned(gt)
        xstep = gt[GDAL_WE_RES]
        xmin = gt[GDAL_TOPLEFT_X]
        xmax = gt[GDAL_TOPLEFT_X] + xstep * (xsize - 1)
        xindex = LinRange(xmin, xmax, xsize)

        ystep = gt[GDAL_NS_RES] # A negative number
        ymax = gt[GDAL_TOPLEFT_Y] + ystep
        ymin = gt[GDAL_TOPLEFT_Y] + ystep * ysize
        yindex = LinRange(ymax, ymin, ysize)

        # Spatial data defaults to area/inteval
        xsampling, ysampling = if _gdalmetadata(raster.ds, "AREA_OR_POINT") == "Point"
            Points(), Points()
        else
            # GeoTiff uses the "pixelCorner" convention
            Intervals(GDAL_X_LOCUS), Intervals(GDAL_Y_LOCUS)
        end

        xmode = Projected(
            order=Ordered(GDAL_X_INDEX, GDAL_X_ARRAY, GDAL_RELATION),
            span=Regular(step(xindex)),
            sampling=xsampling,
            crs=crs,
            mappedcrs=mappedcrs,
        )
        ymode = Projected(
            order=Ordered(GDAL_Y_INDEX, GDAL_Y_ARRAY, GDAL_RELATION),
            sampling=ysampling,
            # Use the range step as is will be different to ystep due to float error
            span=Regular(step(yindex)),
            crs=crs,
            mappedcrs=mappedcrs,
        )
        x = X(xindex; mode=xmode, metadata=xy_metadata)
        y = Y(yindex; mode=ymode, metadata=xy_metadata)

        DimensionalData._formatdims(map(Base.OneTo, (xsize, ysize, nbands)), (x, y, band))
    else
        error("Rotated/transformed dimensions are not handled yet. Open a github issue for GeoData.jl if you need this.")
        # affinemap = geotransform2affine(geotransform)
        # x = X(affinemap; mode=TransformedIndex(dims=X()))
        # y = Y(affinemap; mode=TransformedIndex(dims=Y()))

        # formatdims((xsize, ysize, nbands), (x, y, band))
    end
end

function DD.metadata(raster::AG.RasterDataset, args...)
    band = AG.getband(raster.ds, 1)
    # color = AG.getname(AG.getcolorinterp(band))
    scale = AG.getscale(band)
    offset = AG.getoffset(band)
    # norvw = AG.noverview(band)
    path = first(AG.filelist(raster))
    units = AG.getunittype(band)
    upair = units == "" ? () : (:units=>units,)
    Metadata{:GDAL}(Dict(:filepath=>path, :scale=>scale, :offset=>offset, upair...))
end

function missingval(raster::AG.RasterDataset, args...)
    # We can only handle data where all bands have the same missingval
    band = AG.getband(raster.ds, 1)
    nodata = AG.getnodatavalue(band)
    return if nodata isa Nothing
        nothing
    else
        # Convert in case getnodatavalue is the wrong type.
        # This conversion should always be safe.
        if eltype(band) <: AbstractFloat && nodata isa Real
            convert(eltype(band), nodata)
        else
            nodata
        end
    end
end

# metadata(raster::RasterDataset, key) = begin
#     regex = Regex("$key=(.*)")
#     i = findfirst(f -> occursin(regex, f), meta)
#     if i isa Nothing
#         nothing
#     else
#         match(regex, meta[i])[1]
#     end
# end

crs(raster::AG.RasterDataset, args...) =
    WellKnownText(GeoFormatTypes.CRS(), string(AG.getproj(raster.ds)))

# Utils ########################################################################

function _gdalmetadata(dataset::AG.Dataset, key)
    meta = AG.metadata(dataset)
    regex = Regex("$key=(.*)")
    i = findfirst(f -> occursin(regex, f), meta)
    if i isa Nothing
        nothing
    else
        match(regex, meta[i])[1]
    end
end

function _gdalread(f, filename::AbstractString)
    AG.readraster(filename) do raster
        f(raster)
    end
end

function _gdalwrite(filename, A, nbands, indices;
    driver=AG.extensiondriver(filename), compress="DEFLATE", tiled=true
)
    tiledstring = tiled isa Bool ? (tiled ? "YES" : "NO") : tiled
    kw = (width=size(A, X()), height=size(A, Y()), nbands=nbands, dtype=eltype(A))
    gdaldriver = AG.getdriver(driver)
    if driver == "GTiff"
        options = ["COMPRESS=$compress", "TILED=$tiledstring"]
        AG.create(filename; driver=gdaldriver, options=options, kw...) do dataset
            _gdalsetproperties!(dataset, A)
            AG.write!(dataset, data(A), indices)
        end
    else
        # Create a  memory object and copy it to disk, as ArchGDAL.create
        # does not support direct creation of ASCII etc. rasters
        ArchGDAL.create(""; driver=AG.getdriver("MEM"), kw...) do dataset
            _gdalsetproperties!(dataset, A)
            AG.write!(dataset, data(A), indices)
            AG.copy(dataset; filename=filename, driver=gdaldriver) |> AG.destroy
        end
    end
    return filename
end

function _gdalsetproperties!(dataset, A)
    # Convert the dimensions to `Projected` if they are `Converted`
    # This allows saving NetCDF to Tiff
    # Set the index loci to the start of the cell for the lat and lon dimensions.
    # NetCDF or other formats use the center of the interval, so they need conversion.
    x = convertmode(Projected, DD.maybeshiftlocus(GDAL_X_LOCUS, dims(A, X)))
    y = convertmode(Projected, DD.maybeshiftlocus(GDAL_Y_LOCUS, dims(A, Y)))
    # Convert crs to WKT if it exists
    if !(crs(x) isa Nothing)
        AG.setproj!(dataset, convert(String, convert(WellKnownText, crs(x))))
    end
    # Get the geotransform from the updated lat/lon dims and write
    AG.setgeotransform!(dataset, _dims2geotransform(x, y))

    # Set the nodata value. GDAL can't handle missing. We could choose a default,
    # but we would need to do this for all possible types. `nothing` means
    # there is not missing value.
    # TODO define default nodata values for missing
    if (missingval(A) !== missing) && (missingval(A) !== nothing)
        bands = hasdim(A, Band) ? index(A, Band) : 1
        for i in bands
            AG.setnodatavalue!(AG.getband(dataset, i), missingval(A))
        end
    end

    dataset
end

# Create a GeoArray from a memory-backed dataset
function GeoArray(dataset::AG.Dataset;
    crs=nothing,
    mappedcrs=nothing,
    dims=dims(AG.RasterDataset(dataset), crs, mappedcrs),
    refdims=(),
    name=Symbol(""),
    metadata=metadata(AG.RasterDataset(dataset)),
    missingval=missingval(AG.RasterDataset(dataset))
)
    GeoArray(AG.read(dataset), dims, refdims, name, metadata, missingval)
end

# Create a memory-backed GDAL dataset from any AbstractGeoArray
function unsafe_ArchGDALdataset(A::AbstractGeoArray)
    width = size(A, X)
    height = size(A, Y)
    nbands = size(A, Band)

    dataset = AG.unsafe_create("tmp";
        driver=AG.getdriver("MEM"),
        width=width,
        height=height,
        nbands=nbands,
        dtype=eltype(A)
    )
    # write bands to dataset
    AG.write!(dataset, data(permutedims(A, (X, Y, Band))), Cint[1:nbands...])
    _gdalsetproperties!(dataset, A)
    dataset
end

function AG.Dataset(f::Function, A::AbstractGeoArray)
    dataset = unsafe_ArchGDALdataset(A)
    try
        f(dataset)
    finally
        AG.destroy(dataset)
    end
end

#= Geotranforms ########################################################################

See https://lists.osgeo.org/pipermail/gdal-dev/2011-July/029449.html

"In the particular, but common, case of a “north up” image without any rotation or
shearing, the georeferencing transform takes the following form" :
adfGeoTransform[0] /* top left x */
adfGeoTransform[1] /* w-e pixel resolution */
adfGeoTransform[2] /* 0 */
adfGeoTransform[3] /* top left y */
adfGeoTransform[4] /* 0 */
adfGeoTransform[5] /* n-s pixel resolution (negative value) */
=#

const GDAL_EMPTY_TRANSFORM = [0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
const GDAL_TOPLEFT_X = 1
const GDAL_WE_RES = 2
const GDAL_ROT1 = 3
const GDAL_TOPLEFT_Y = 4
const GDAL_ROT2 = 5
const GDAL_NS_RES = 6

_isalligned(geotransform) = geotransform[GDAL_ROT1] == 0 && geotransform[GDAL_ROT2] == 0

_geotransform2affine(gt) =
    AffineMap([gt[GDAL_WE_RES] gt[GDAL_ROT1]; gt[GDAL_ROT2] gt[GDAL_NS_RES]],
              [gt[GDAL_TOPLEFT_X], gt[GDAL_TOPLEFT_Y]])

function _dims2geotransform(x::X, y::Y)
    gt = zeros(6)
    gt[GDAL_TOPLEFT_X] = first(x)
    gt[GDAL_WE_RES] = step(x)
    gt[GDAL_ROT1] = 0.0
    gt[GDAL_TOPLEFT_Y] = first(y) - step(y)
    gt[GDAL_ROT2] = 0.0
    gt[GDAL_NS_RES] = step(y)
    return gt
end


text:
"""
ДАНО: Робот находится в произвольной клетке ограниченного прямоугольного поля без внутренних перегородок и маркеров.
РЕЗУЛЬТАТ: Робот — в исходном положении в центре прямого креста из маркеров, расставленных вплоть до внешней рамки.
"""

function cross!(r::Robot) #Основная функция
    for side in (Nord,West,Sud,Ost) #Перебор всех возможных направлений
        putmarkers!(r,side)
        move_by_markers!(r,inverse(side))
    end
    putmarker!(r)
end

function putmarkers!(r::Robot,side::HorizonSide)  #Идя в заданном направлении, Робот ставит маркеры, пока не упрётся в стену (то есть пока в заданном направлении нет стены)
    while isborder(r,side)==false
        move!(r,side)
        putmarker!(r)
    end
end

function move_by_markers!(r::Robot,side::HorizonSide)  #Робот Двигается по маркерам
    while ismarker(r)==true
        move!(r,side)
    end
end

inverse(side::HorizonSide) = HorizonSide(mod(Int(side)+2, 4))  #Функция задаёт направление, противоположное данному

"""
В строках 10, 18 и 25 начинаются описания вспомогательных функций.
"""

max_stars_repo_name: alexgruzberg/mirea-progs

max_issues_repo_name: alexgruzberg/mirea-progs

max_forks_repo_name: alexgruzberg/mirea-progs


text:
@inline function abs_dd_dd(x::Tuple{T,T}) where {T<:AbstractFloat}
    hi, lo = x
    if signbit(hi)
        hi = -hi
        lo = -lo
    end
    return hi, lo
end

@inline function neg_dd_dd(x::Tuple{T,T}) where {T<:AbstractFloat}
    hi, lo = x
    hi = -hi
    lo = -lo
    return hi, lo
end

@inline function negabs_dd_dd(x::Tuple{T,T}) where {T<:AbstractFloat}
    hi, lo = x
    if !signbit(hi)
        hi = -hi
        lo = -lo
    end
    return hi, lo
end

function inv_dd_dd(x::Tuple{T,T}) where {T<:AbstractFloat}
    return DWInvDW3(HI(x), LO(x))
end
function inv_dd_dd_fast(x::Tuple{T,T}) where {T<:AbstractFloat}
    return DWInvDW2(HI(x), LO(x))
end

#=
   slightly faster than inv_dd_dd_fast, without analytic relerr though
@inline function inv_dd_dd_fast(y::Tuple{T,T}) where {T<:AbstractFloat}
    xhi, xlo = one(T), zero(T)
    yhi, ylo = y
    hi = xhi / yhi
    uh, ul = mul_2(hi, yhi)
    lo = ((((xhi - uh) - ul) + xlo) - hi*ylo)/yhi
    hi, lo = add_2(hi, lo)
    return hi, lo
end
=#
@inline function square_dd_dd(a::Tuple{T,T}) where {T<:AbstractFloat}
    ahi, alo = HILO(a)
    zhi = ahi * ahi
    zlo = fma(ahi, ahi, -zhi)
    zlo += (ahi * alo) * 2
    zlo += alo * alo
    return zhi, zlo
end

@inline function cube_dd_dd(a::Tuple{T,T}) where {T<:AbstractFloat}
    zhi, zlo = square_dd_dd(a)
    zhi, zlo = mul_dddd_dd((zhi, zlo), a)
    return zhi, zlo
end

function sqrt_dd_dd(x::Tuple{T,T}) where {T<:AbstractFloat}
    iszero(HI(x)) && return x
    signbit(HI(x)) && throw(DomainError("sqrt(x) expects x >= 0"))

    half = T(0.5)
    dhalf = (half, zero(T))

    r = inv(sqrt(HI(x)))
    h = (HI(x) * half, LO(x) * half)

    r2 = mul_fpfp_dd(r, r)
    hr2 = mul_dddd_dd(h, r2)
    radj = sub_fpdd_dd(half, hr2)
    radj = mul_ddfp_dd(radj, r)
    r = add_fpdd_dd(r, radj)

    r2 = mul_dddd_dd(r, r)
    hr2 = mul_dddd_dd(h, r2)
    radj = sub_fpdd_dd(half, hr2)
    radj = mul_dddd_dd(radj, r)
    r = add_dddd_dd(r, radj)

    r = mul_dddd_dd(r, x)

    return r
end

#=
    from mpfun90.f
cuberoot(A) = A * inv(cuberootsquared(A))
invcuberootsquared(A) is found iteratively using Newton's method with a final approximation from Karp
   adjustx[k] = ((1- x[k]^3 * A^2)*x[k])/3
   x[k+1] = x[k] + adjustx[k]
    ...
   x[n-1] = x[n-2] + adjustx[n-2]
   x[n] = x[n-1] + adjustx[n-1]
   cuberoot(A) = (A * x[n]) + ((A - (A * x[n])^3)) * x[n] / 3)
# x := x - x * (x^3 - a) / (2*x^3 + y)
# x -= ( x - (z/(x*x)))*(1/3)
=#

function cbrt_dd_dd(a::Tuple{T,T}) where {T<:AbstractFloat}
    hi, lo = HILO(a)
    a2 = mul_dddd_dd(a,a)
    one1 = one(T)
    onethird = (0.3333333333333333, 1.850371707708594e-17)

    a_inv = inv_dd_dd(a)
    tmp = cbrt(HI(a_inv))
    # initial approximation to a^(-2/3)
    x = mul_fpfp_dd(tmp, tmp)

    x3 = mul_dddd_dd(x,x)
    x3 = mul_dddd_dd(x3, x)
    x3 = mul_dddd_dd(x3, a2)
    x3 = sub_fpdd_dd(one1, x3)
    x3 = mul_dddd_dd(x3, x)
    x3 = mul_dddd_dd(x3, onethird)

    x = add_dddd_dd(x, x3)

    x3 = mul_dddd_dd(x,x)
    x3 = mul_dddd_dd(x3, x)
    x3 = mul_dddd_dd(x3, a2)
    x3 = sub_fpdd_dd(one1, x3)
    x3 = mul_dddd_dd(x3, x)
    x3 = mul_dddd_dd(x3, onethird)

    x = add_dddd_dd(x, x3)

    x3 = mul_dddd_dd(x,x)
    x3 = mul_dddd_dd(x3, x)
    x3 = mul_dddd_dd(x3, a2)
    x3 = sub_fpdd_dd(one1, x3)
    x3 = mul_dddd_dd(x3, x)
    x3 = mul_dddd_dd(x3, onethird)

    x = add_dddd_dd(x, x3)

    ax = mul_dddd_dd(a, x)

    # return ax
    #                                      amax3
    #                    ax      +     (a  -  ax^3) * x /3
    #   cuberoot(A) = (A * x[n]) + ((A - (A * x[n])^3)) * x[n] / 3)

    ax3 = mul_dddd_dd(ax, ax)
    ax3 = mul_dddd_dd(ax3, ax)
    amax3 = sub_dddd_dd(a, ax3)
    amax3 = mul_dddd_dd(amax3, x)
    amax3 = mul_dddd_dd(amax3, onethird)

    ax = add_dddd_dd(ax, amax3)
    return ax
end


text:
function create_component_macros{N}(name::ASCIIString, sizes::NTuple{N, Int})
    unpack_name = symbol("unpack_" * name)
    unpack_name! = symbol("unpack_" * name * "!")
    pack_name = symbol("pack_" * name)
    @eval begin
        macro $unpack_name(ex)
            length($sizes) == length(ex.args[1].args) || throw(ArgumentError("Invalid number of output vars"))
            sizes_loc = $sizes
            exvec = quote end
            x_val = gensym("x")
            push!(exvec.args, :($x_val = $(ex.args[end])))
            length_cum = 1
            for (i, sym) in enumerate(ex.args[1].args)
                assert(typeof(sym) == Symbol)
                if $(sizes)[i] == 1
                    push!(exvec.args, :($(sym) = $(x_val)[$(length_cum)]))
                else
                    push!(exvec.args, :($(sym) = $(x_val)[$length_cum : $(length_cum + sizes_loc[i]-1)]))
                end
                length_cum += $(sizes)[i]
            end
            return esc(exvec)
        end

         macro $unpack_name!(ex)
            length($sizes) == length(ex.args[1].args) || throw(ArgumentError("Invalid number of output vars"))
            sizes_loc = $sizes
            exvec = quote end
            x_val = gensym("x")
            push!(exvec.args, :($x_val = $(ex.args[end])))
            length_cum = 1
            for (i, sym) in enumerate(ex.args[1].args)
                assert(typeof(sym) == Symbol)
                if $(sizes)[i] == 1
                    push!(exvec.args, :($(sym) = $(x_val)[$(length_cum)]))
                else
                    push!(exvec.args, :(@devec $(sym)[:] = $(x_val)[$length_cum : $(length_cum + sizes_loc[i]-1)]))
                end
                length_cum += $(sizes)[i]
            end
            return esc(exvec)
        end

        macro $pack_name(ex)
            length($sizes) == length(ex.args[2].args) || throw(ArgumentError("Invalid number of output vars"))
            sizes_loc = $sizes
            exvec = quote end
            length_cum = 1
            x = ex.args[1]
            for (i, sym) in enumerate(ex.args[2].args)
                assert(typeof(sym) == Symbol)
                if $(sizes)[i] == 1
                    push!(exvec.args, :($(x)[$(length_cum)] = $(sym)))
                else
                    push!(exvec.args, :($(x)[$length_cum : $(length_cum + sizes_loc[i]-1)] = $(sym)))
                end
                length_cum += $(sizes)[i]
            end
            return esc(exvec)
        end

    end
end


text:
@testset "Testing General submodule" begin

  fnMeasBruker = joinpath(datadir,"BrukerStore","20150915_102110_Wuerfelphantom_1_1","18")
  fnSMBruker = joinpath(datadir,"BrukerStore","20141121_130749_CalibrationScans_1_1","76")
  fnSM1DBruker = joinpath(datadir,"BrukerStore","20170807_142514_Service_1_1","89")

  fnMeasV1 = joinpath(tmpdir,"mdf","measurement_V1.mdf")
  fnMeasV2 = joinpath(tmpdir,"mdf","measurement_V2.mdf")
  fnSMV1 = joinpath(tmpdir,"mdf","systemMatrix_V1.mdf")
  fnSMV2 = joinpath(tmpdir,"mdf","systemMatrix_V2.mdf")
  fnSMV3 = joinpath(tmpdir,"mdf","systemMatrix_V3.mdf")
  fnSMV4 = joinpath(tmpdir,"mdf","systemMatrix_V4.mdf")
  fnSMV5 = joinpath(tmpdir,"mdf","systemMatrix_V5.mdf")
  fnSMV6 = joinpath(tmpdir,"mdf","systemMatrix_V6.mdf")
  fnSM1DV1 = joinpath(tmpdir,"mdf","systemMatrix1D_V1.mdf")
  fnSM1DV2 = joinpath(tmpdir,"mdf","systemMatrix1D_V2.mdf")

  @testset "Fields" begin
    measBruker = MPIFile(fnMeasBruker)
    @test typeof(measBruker) == BrukerFileMeas
    saveasMDF(fnMeasV2, measBruker)#, frames=1:100) <- TODO test this

    mdfv2 = MPIFile(fnMeasV2)
    @test typeof(mdfv2) <: MDFFileV2

    mdfv2InMemory = MDFv2InMemory(mdfv2)
    @test typeof(mdfv2InMemory) <: MDFv2InMemory

    for mdf in (measBruker, mdfv2, mdfv2InMemory)
      @info "Test $mdf"

      @test length(mdf) == 1

      @test studyName(mdf) == "Wuerfelphantom"
      @test studyNumber(mdf) == 1
      @test studyDescription(mdf) == "n.a."
      @test studyTime(mdf) == DateTime( "2015-09-15T10:21:10.992" )
      @test studyUuid(mdf) == studyUuid(measBruker) #UUID("fe9635e0-7cfc-44eb-9f5a-585013a2cb51") see issue #32 on github
      @test uuid_version(studyUuid(mdf)) == 4

      @test experimentName(mdf) == "fuenf (E18)"
      @test experimentNumber(mdf) == 18
      @test experimentDescription(mdf) == "fuenf (E18)"
      @test experimentSubject(mdf) == "Wuerfelphantom"
      @test experimentIsSimulation(mdf) == false
      @test experimentIsCalibration(mdf) == false
      @test experimentUuid(mdf) == experimentUuid(measBruker) # UUID("ef2110b5-cce8-4ca6-b041-38ea01254c47")
      @test uuid_version(experimentUuid(mdf)) == 4

      @test scannerFacility(mdf) == "Universitätsklinikum Hamburg Eppendorf"
      @test scannerOperator(mdf) == "nmrsu"
      @test scannerManufacturer(mdf) == "Bruker/Philips"
      @test scannerName(mdf) == "Preclinical MPI System"
      @test scannerTopology(mdf) == "FFP"

      @test tracerName(mdf) == ["Resovist"]
      @test tracerBatch(mdf) == ["0"]
      @test tracerVendor(mdf) == ["n.a."]
      @test tracerVolume(mdf) == [0.0]
      @test tracerConcentration(mdf) == [0.5]
      @test tracerInjectionTime(mdf) == [DateTime("2015-09-15T11:17:23.011")]

      @test acqStartTime(mdf) == DateTime("2015-09-15T11:17:23.011")
      @test acqGradient(mdf)[:,:,1,1] == [-1.25 0 0; 0 -1.25 0;0 0 2.5]
      @test acqFramePeriod(mdf) == 6.528E-4
      @test acqNumPeriodsPerFrame(mdf) == 1
      @test acqOffsetFieldShift(mdf)[:,1,1] == [0.0; 0.0; -0.0]

      @test dfNumChannels(mdf) == 3
      @test dfWaveform(mdf) == fill("sine", (1, 1))
      @test dfStrength(mdf)[:,:,1] == [0.014 0.014 0.0]
      @test dfPhase(mdf)[:,:,1] == [1.5707963267948966 1.5707963267948966 1.5707963267948966]
      @test dfBaseFrequency(mdf) == 2500000.0
      @test dfDivider(mdf)[:,1] == [102; 96; 99]
      @test dfCycle(mdf) == 6.528E-4

      @test rxNumChannels(mdf) == 3
      @test rxBandwidth(mdf) == 1250000.0
      @test rxNumSamplingPoints(mdf) == 1632
      @test rxDataConversionFactor(mdf) == repeat([1.0, 0.0], outer=(1,rxNumChannels(mdf)))
      @test size(rxTransferFunction(mdf)) == (817, 3)
      @test occursin("../../../transferFunction/tf.h5", rxTransferFunctionFileName(mdf))

      @test acqNumAverages(mdf) == 1
      @test acqNumFrames(mdf) == 500
      @test acqNumPeriodsPerFrame(mdf) == 1
      @test acqNumPeriods(mdf) == 500
      @test acqNumPatches(mdf) == 1
      @test acqNumPeriodsPerPatch(mdf) == 1

      @test size( measData(mdf) ) == (1632,3,1,500)
      @test size( measDataTDPeriods(mdf) ) == (1632,3,500)
      @test size( measDataTDPeriods(mdf, 101:200) ) == (1632,3,100)

      N = acqNumFrames(mdf)

      @test size(getMeasurements(mdf, numAverages=1,
                  spectralLeakageCorrection=false)) == (1632,3,1,500)

      @test size(getMeasurements(mdf, numAverages=10,
                  spectralLeakageCorrection=false)) == (1632,3,1,50)

      @test size(getMeasurements(mdf, numAverages=10, frames=1:100,
                  spectralLeakageCorrection=true)) == (1632,3,1,10)

      @test size(getMeasurementsFD(mdf, numAverages=10, frames=1:100)) == (817,3,1,10)

      @test size(getMeasurementsFD(mdf, numAverages=10, frames=1:100, loadasreal=true)) == (1634,3,1,10)

      @test size(getMeasurementsFD(mdf,frequencies=1:10, numAverages=10)) == (10,1,50)

    end
  end

  # Calibration File

  @testset "Calibration" begin
    smBruker = MPIFile(fnSMBruker)
    @test typeof(smBruker) == BrukerFileCalib
    saveasMDF(fnSMV2, smBruker)

    smv2 = MPIFile(fnSMV2)
    @test typeof(smv2) <: MDFFileV2

    smInMemory = MDFv2InMemory(smv2)
    @test typeof(smInMemory) <: MDFv2InMemory

    smBrukerPretendToBeMeas = MPIFile(fnSMBruker, isCalib=false)
    saveasMDF(fnSMV3, smBrukerPretendToBeMeas, applyCalibPostprocessing=true)

    smv3 = MPIFile(fnSMV3)
    @test typeof(smv3) <: MDFFileV2

    # Bruker specific test
    @test rawDataLengthConsistent(smBruker)

    for sm in (smBruker,smv2,smv3,smInMemory)
      @info "Test $sm"

      @test size( systemMatrixWithBG(sm) ) == (1959,817,3,1)
      @test size( systemMatrix(sm,1:10) ) == (1936,10)

      @test size(rxTransferFunction(sm)) == (817, 3)
      @test rxHasTransferFunction(sm) == true
      @test occursin("../../../transferFunction/tf.h5", rxTransferFunctionFileName(sm))

      @test measIsFourierTransformed(sm) == true
      @test measIsTFCorrected(sm) == false
      @test measIsFastFrameAxis(sm) == true
      @test measIsBGCorrected(sm) == false
      @test measFramePermutation(sm) == [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 89, 88, 87, 86, 85, 84, 83, 82, 81, 80, 79, 78, 77, 76, 75, 74, 73, 72, 71, 70, 69, 68, 67, 66, 65, 64, 63, 62, 61, 60, 59, 58, 57, 56, 55, 54, 53, 52, 51, 50, 49, 48, 47, 46, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 178, 177, 176, 175, 174, 173, 172, 171, 170, 169, 168, 167, 166, 165, 164, 163, 162, 161, 160, 159, 158, 157, 156, 155, 154, 153, 152, 151, 150, 149, 148, 147, 146, 145, 144, 143, 142, 141, 140, 139, 138, 137, 136, 135, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 267, 266, 265, 264, 263, 262, 261, 260, 259, 258, 257, 256, 255, 254, 253, 252, 251, 250, 249, 248, 247, 246, 245, 244, 243, 242, 241, 240, 239, 238, 237, 236, 235, 234, 233, 232, 231, 230, 229, 228, 227, 226, 225, 224, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 356, 355, 354, 353, 352, 351, 350, 349, 348, 347, 346, 345, 344, 343, 342, 341, 340, 339, 338, 337, 336, 335, 334, 333, 332, 331, 330, 329, 328, 327, 326, 325, 324, 323, 322, 321, 320, 319, 318, 317, 316, 315, 314, 313, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401, 445, 444, 443, 442, 441, 440, 439, 438, 437, 436, 435, 434, 433, 432, 431, 430, 429, 428, 427, 426, 425, 424, 423, 422, 421, 420, 419, 418, 417, 416, 415, 414, 413, 412, 411, 410, 409, 408, 407, 406, 405, 404, 403, 402, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 534, 533, 532, 531, 530, 529, 528, 527, 526, 525, 524, 523, 522, 521, 520, 519, 518, 517, 516, 515, 514, 513, 512, 511, 510, 509, 508, 507, 506, 505, 504, 503, 502, 501, 500, 499, 498, 497, 496, 495, 494, 493, 492, 491, 536, 537, 538, 539, 540, 541, 542, 543, 544, 545, 546, 547, 548, 549, 550, 551, 552, 553, 554, 555, 556, 557, 558, 559, 560, 561, 562, 563, 564, 565, 566, 567, 568, 569, 570, 571, 572, 573, 574, 575, 576, 577, 578, 579, 623, 622, 621, 620, 619, 618, 617, 616, 615, 614, 613, 612, 611, 610, 609, 608, 607, 606, 605, 604, 603, 602, 601, 600, 599, 598, 597, 596, 595, 594, 593, 592, 591, 590, 589, 588, 587, 586, 585, 584, 583, 582, 581, 580, 625, 626, 627, 628, 629, 630, 631, 632, 633, 634, 635, 636, 637, 638, 639, 640, 641, 642, 643, 644, 645, 646, 647, 648, 649, 650, 651, 652, 653, 654, 655, 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666, 667, 668, 712, 711, 710, 709, 708, 707, 706, 705, 704, 703, 702, 701, 700, 699, 698, 697, 696, 695, 694, 693, 692, 691, 690, 689, 688, 687, 686, 685, 684, 683, 682, 681, 680, 679, 678, 677, 676, 675, 674, 673, 672, 671, 670, 669, 714, 715, 716, 717, 718, 719, 720, 721, 722, 723, 724, 725, 726, 727, 728, 729, 730, 731, 732, 733, 734, 735, 736, 737, 738, 739, 740, 741, 742, 743, 744, 745, 746, 747, 748, 749, 750, 751, 752, 753, 754, 755, 756, 757, 801, 800, 799, 798, 797, 796, 795, 794, 793, 792, 791, 790, 789, 788, 787, 786, 785, 784, 783, 782, 781, 780, 779, 778, 777, 776, 775, 774, 773, 772, 771, 770, 769, 768, 767, 766, 765, 764, 763, 762, 761, 760, 759, 758, 803, 804, 805, 806, 807, 808, 809, 810, 811, 812, 813, 814, 815, 816, 817, 818, 819, 820, 821, 822, 823, 824, 825, 826, 827, 828, 829, 830, 831, 832, 833, 834, 835, 836, 837, 838, 839, 840, 841, 842, 843, 844, 845, 846, 890, 889, 888, 887, 886, 885, 884, 883, 882, 881, 880, 879, 878, 877, 876, 875, 874, 873, 872, 871, 870, 869, 868, 867, 866, 865, 864, 863, 862, 861, 860, 859, 858, 857, 856, 855, 854, 853, 852, 851, 850, 849, 848, 847, 892, 893, 894, 895, 896, 897, 898, 899, 900, 901, 902, 903, 904, 905, 906, 907, 908, 909, 910, 911, 912, 913, 914, 915, 916, 917, 918, 919, 920, 921, 922, 923, 924, 925, 926, 927, 928, 929, 930, 931, 932, 933, 934, 935, 979, 978, 977, 976, 975, 974, 973, 972, 971, 970, 969, 968, 967, 966, 965, 964, 963, 962, 961, 960, 959, 958, 957, 956, 955, 954, 953, 952, 951, 950, 949, 948, 947, 946, 945, 944, 943, 942, 941, 940, 939, 938, 937, 936, 981, 982, 983, 984, 985, 986, 987, 988, 989, 990, 991, 992, 993, 994, 995, 996, 997, 998, 999, 1000, 1001, 1002, 1003, 1004, 1005, 1006, 1007, 1008, 1009, 1010, 1011, 1012, 1013, 1014, 1015, 1016, 1017, 1018, 1019, 1020, 1021, 1022, 1023, 1024, 1068, 1067, 1066, 1065, 1064, 1063, 1062, 1061, 1060, 1059, 1058, 1057, 1056, 1055, 1054, 1053, 1052, 1051, 1050, 1049, 1048, 1047, 1046, 1045, 1044, 1043, 1042, 1041, 1040, 1039, 1038, 1037, 1036, 1035, 1034, 1033, 1032, 1031, 1030, 1029, 1028, 1027, 1026, 1025, 1070, 1071, 1072, 1073, 1074, 1075, 1076, 1077, 1078, 1079, 1080, 1081, 1082, 1083, 1084, 1085, 1086, 1087, 1088, 1089, 1090, 1091, 1092, 1093, 1094, 1095, 1096, 1097, 1098, 1099, 1100, 1101, 1102, 1103, 1104, 1105, 1106, 1107, 1108, 1109, 1110, 1111, 1112, 1113, 1157, 1156, 1155, 1154, 1153, 1152, 1151, 1150, 1149, 1148, 1147, 1146, 1145, 1144, 1143, 1142, 1141, 1140, 1139, 1138, 1137, 1136, 1135, 1134, 1133, 1132, 1131, 1130, 1129, 1128, 1127, 1126, 1125, 1124, 1123, 1122, 1121, 1120, 1119, 1118, 1117, 1116, 1115, 1114, 1159, 1160, 1161, 1162, 1163, 1164, 1165, 1166, 1167, 1168, 1169, 1170, 1171, 1172, 1173, 1174, 1175, 1176, 1177, 1178, 1179, 1180, 1181, 1182, 1183, 1184, 1185, 1186, 1187, 1188, 1189, 1190, 1191, 1192, 1193, 1194, 1195, 1196, 1197, 1198, 1199, 1200, 1201, 1202, 1246, 1245, 1244, 1243, 1242, 1241, 1240, 1239, 1238, 1237, 1236, 1235, 1234, 1233, 1232, 1231, 1230, 1229, 1228, 1227, 1226, 1225, 1224, 1223, 1222, 1221, 1220, 1219, 1218, 1217, 1216, 1215, 1214, 1213, 1212, 1211, 1210, 1209, 1208, 1207, 1206, 1205, 1204, 1203, 1248, 1249, 1250, 1251, 1252, 1253, 1254, 1255, 1256, 1257, 1258, 1259, 1260, 1261, 1262, 1263, 1264, 1265, 1266, 1267, 1268, 1269, 1270, 1271, 1272, 1273, 1274, 1275, 1276, 1277, 1278, 1279, 1280, 1281, 1282, 1283, 1284, 1285, 1286, 1287, 1288, 1289, 1290, 1291, 1335, 1334, 1333, 1332, 1331, 1330, 1329, 1328, 1327, 1326, 1325, 1324, 1323, 1322, 1321, 1320, 1319, 1318, 1317, 1316, 1315, 1314, 1313, 1312, 1311, 1310, 1309, 1308, 1307, 1306, 1305, 1304, 1303, 1302, 1301, 1300, 1299, 1298, 1297, 1296, 1295, 1294, 1293, 1292, 1337, 1338, 1339, 1340, 1341, 1342, 1343, 1344, 1345, 1346, 1347, 1348, 1349, 1350, 1351, 1352, 1353, 1354, 1355, 1356, 1357, 1358, 1359, 1360, 1361, 1362, 1363, 1364, 1365, 1366, 1367, 1368, 1369, 1370, 1371, 1372, 1373, 1374, 1375, 1376, 1377, 1378, 1379, 1380, 1424, 1423, 1422, 1421, 1420, 1419, 1418, 1417, 1416, 1415, 1414, 1413, 1412, 1411, 1410, 1409, 1408, 1407, 1406, 1405, 1404, 1403, 1402, 1401, 1400, 1399, 1398, 1397, 1396, 1395, 1394, 1393, 1392, 1391, 1390, 1389, 1388, 1387, 1386, 1385, 1384, 1383, 1382, 1381, 1426, 1427, 1428, 1429, 1430, 1431, 1432, 1433, 1434, 1435, 1436, 1437, 1438, 1439, 1440, 1441, 1442, 1443, 1444, 1445, 1446, 1447, 1448, 1449, 1450, 1451, 1452, 1453, 1454, 1455, 1456, 1457, 1458, 1459, 1460, 1461, 1462, 1463, 1464, 1465, 1466, 1467, 1468, 1469, 1513, 1512, 1511, 1510, 1509, 1508, 1507, 1506, 1505, 1504, 1503, 1502, 1501, 1500, 1499, 1498, 1497, 1496, 1495, 1494, 1493, 1492, 1491, 1490, 1489, 1488, 1487, 1486, 1485, 1484, 1483, 1482, 1481, 1480, 1479, 1478, 1477, 1476, 1475, 1474, 1473, 1472, 1471, 1470, 1515, 1516, 1517, 1518, 1519, 1520, 1521, 1522, 1523, 1524, 1525, 1526, 1527, 1528, 1529, 1530, 1531, 1532, 1533, 1534, 1535, 1536, 1537, 1538, 1539, 1540, 1541, 1542, 1543, 1544, 1545, 1546, 1547, 1548, 1549, 1550, 1551, 1552, 1553, 1554, 1555, 1556, 1557, 1558, 1602, 1601, 1600, 1599, 1598, 1597, 1596, 1595, 1594, 1593, 1592, 1591, 1590, 1589, 1588, 1587, 1586, 1585, 1584, 1583, 1582, 1581, 1580, 1579, 1578, 1577, 1576, 1575, 1574, 1573, 1572, 1571, 1570, 1569, 1568, 1567, 1566, 1565, 1564, 1563, 1562, 1561, 1560, 1559, 1604, 1605, 1606, 1607, 1608, 1609, 1610, 1611, 1612, 1613, 1614, 1615, 1616, 1617, 1618, 1619, 1620, 1621, 1622, 1623, 1624, 1625, 1626, 1627, 1628, 1629, 1630, 1631, 1632, 1633, 1634, 1635, 1636, 1637, 1638, 1639, 1640, 1641, 1642, 1643, 1644, 1645, 1646, 1647, 1691, 1690, 1689, 1688, 1687, 1686, 1685, 1684, 1683, 1682, 1681, 1680, 1679, 1678, 1677, 1676, 1675, 1674, 1673, 1672, 1671, 1670, 1669, 1668, 1667, 1666, 1665, 1664, 1663, 1662, 1661, 1660, 1659, 1658, 1657, 1656, 1655, 1654, 1653, 1652, 1651, 1650, 1649, 1648, 1693, 1694, 1695, 1696, 1697, 1698, 1699, 1700, 1701, 1702, 1703, 1704, 1705, 1706, 1707, 1708, 1709, 1710, 1711, 1712, 1713, 1714, 1715, 1716, 1717, 1718, 1719, 1720, 1721, 1722, 1723, 1724, 1725, 1726, 1727, 1728, 1729, 1730, 1731, 1732, 1733, 1734, 1735, 1736, 1780, 1779, 1778, 1777, 1776, 1775, 1774, 1773, 1772, 1771, 1770, 1769, 1768, 1767, 1766, 1765, 1764, 1763, 1762, 1761, 1760, 1759, 1758, 1757, 1756, 1755, 1754, 1753, 1752, 1751, 1750, 1749, 1748, 1747, 1746, 1745, 1744, 1743, 1742, 1741, 1740, 1739, 1738, 1737, 1782, 1783, 1784, 1785, 1786, 1787, 1788, 1789, 1790, 1791, 1792, 1793, 1794, 1795, 1796, 1797, 1798, 1799, 1800, 1801, 1802, 1803, 1804, 1805, 1806, 1807, 1808, 1809, 1810, 1811, 1812, 1813, 1814, 1815, 1816, 1817, 1818, 1819, 1820, 1821, 1822, 1823, 1824, 1825, 1869, 1868, 1867, 1866, 1865, 1864, 1863, 1862, 1861, 1860, 1859, 1858, 1857, 1856, 1855, 1854, 1853, 1852, 1851, 1850, 1849, 1848, 1847, 1846, 1845, 1844, 1843, 1842, 1841, 1840, 1839, 1838, 1837, 1836, 1835, 1834, 1833, 1832, 1831, 1830, 1829, 1828, 1827, 1826, 1871, 1872, 1873, 1874, 1875, 1876, 1877, 1878, 1879, 1880, 1881, 1882, 1883, 1884, 1885, 1886, 1887, 1888, 1889, 1890, 1891, 1892, 1893, 1894, 1895, 1896, 1897, 1898, 1899, 1900, 1901, 1902, 1903, 1904, 1905, 1906, 1907, 1908, 1909, 1910, 1911, 1912, 1913, 1914, 1958, 1957, 1956, 1955, 1954, 1953, 1952, 1951, 1950, 1949, 1948, 1947, 1946, 1945, 1944, 1943, 1942, 1941, 1940, 1939, 1938, 1937, 1936, 1935, 1934, 1933, 1932, 1931, 1930, 1929, 1928, 1927, 1926, 1925, 1924, 1923, 1922, 1921, 1920, 1919, 1918, 1917, 1916, 1915, 1, 90, 179, 268, 357, 446, 535, 624, 713, 802, 891, 980, 1069, 1158, 1247, 1336, 1425, 1514, 1603, 1692, 1781, 1870, 1959]

      @test size( calibSNR(sm) ) == (817,3,1)
      @test calibFov(sm) == [0.044; 0.044; 0.001]
      @test calibFovCenter(sm) == [0.0; -0.0; 0.0]
      @test calibSize(sm) == [44; 44; 1]
      @test calibOrder(sm) == "xyz"
      @test calibDeltaSampleSize(sm) == [0.0, 0.0, 0.0] #[0.001; 0.001; 0.001]
      @test calibMethod(sm) == "robot"

      @test size(filterFrequencies(sm, SNRThresh = 5)) == (147,)
      #@test size(filterFrequencies(sm, numUsedFreqs = 100)) == (100,) # not working

      @test size(getSystemMatrix(sm,1:10)) == (1936,10)
      @test size(getSystemMatrix(sm,1:10,loadasreal=true)) == (1936,20)
      @test size(getSystemMatrix(sm,1:10,bgCorrection=true)) == (1936,10)
      # test on the data level if the conversion was successful
      SNRThresh = 2
      freq = filterFrequencies(smBruker,SNRThresh=SNRThresh)
      SBruker = getSystemMatrix(smBruker,frequencies=freq)
      S = getSystemMatrix(sm,frequencies=freq)

      relativeDeviation = zeros(Float32,length(freq))
      for f in 1:length(freq)
        relativeDeviation[f] = norm(SBruker[:,f]-S[:,f])/norm(SBruker[:,f])
      end
      # test if relative deviation for most of the frequency components is below 0.003
      @test quantile(relativeDeviation,0.95)<0.003
    end

    # Next test checks if the cached system matrix is the same as the one loaded
    # from the raw data
    S_loadedfromraw = getMeasurementsFD(smBrukerPretendToBeMeas,
          frames=1:acqNumFGFrames(smBrukerPretendToBeMeas),sortFrames=true,
          spectralLeakageCorrection=false,transposed=true,tfCorrection=false)

    S_loadedfromproc = systemMatrix(smBruker)

    @test norm(vec(S_loadedfromraw-S_loadedfromproc)) / norm(vec(S_loadedfromproc)) < 1e-6

    # Compression of calib files
    saveasMDF(fnSMV4, smv2, SNRThresh=0.0)
    smv4 = MPIFile(fnSMV4)
    @test size(calibSNR(smv4), 1) == 817
    @test length(filterFrequencies(smv4)) == 817*3

    saveasMDF(fnSMV5, smv2, SNRThresh=4.0)
    smv5 = MPIFile(fnSMV5)
    @test size(calibSNR(smv5), 1) == 110
    @test length(filterFrequencies(smv5)) == 110*3
    @test length(filterFrequencies(smv5, SNRThresh=4.0)) == 194

    @test size(measData(smv5)) == (1959, 110, 3, 1)

    freq = filterFrequencies(smv5, SNRThresh=4.0)
    S1 = getSystemMatrix(smv5, freq)
    @test size(S1)  == (1936,194)

    S2 = getSystemMatrix(smv2, freq)
    @test S1  == S2

    saveasMDF(fnSMV6, smv2, SNRThresh=4.0, sparsityTrafoRedFactor=0.2)
    smv6 = MPIFile(fnSMV6)

    S1 = getSystemMatrix(smv6, freq)

    relativeDeviation = zeros(Float32,length(freq))
    for f in 1:length(freq)
      relativeDeviation[f] = norm(S1[:,f]-S2[:,f])/norm(S2[:,f])
    end
    # test if relative deviation for most of the frequency components is below 0.003
    @test quantile(relativeDeviation,0.95)<0.14

    # Calibration file 1D

    # Calibration File

    sm1DBruker = MPIFile(fnSM1DBruker)
    @test typeof(sm1DBruker) == BrukerFileCalib
    saveasMDF(fnSM1DV1, sm1DBruker)

    sm1D = MPIFile(fnSM1DV1)
    @test typeof(sm1D) <: MDFFileV2

    smInMemory1D = MDFv2InMemory(sm1D)
    @test typeof(smInMemory1D) <: MDFv2InMemory

    for sm in (sm1DBruker,sm1D, smInMemory1D)
      @info "Test $sm"

      @test size( systemMatrixWithBG(sm) ) == (67,52,3,1)
      @test size( systemMatrix(sm,1:10) ) == (60,10)
      @test size( systemMatrix(sm) ) == (60,52,3,1)

      @test size(calibSNR(sm)) == (52, 3, 1)

      @test rxNumSamplingPoints(sm) == 102
      @test rxNumFrequencies(sm) == 52
    end

    sm1DBrukerMeas = MPIFile(fnSM1DBruker, isCalib=false)
    saveasMDF(fnSM1DV2, sm1DBrukerMeas)

    sm1DMeas = MPIFile(fnSM1DV2)
    @test typeof(sm1DMeas) <: MDFFileV2

    sm1DMeasInMemory = MDFv2InMemory(sm1DMeas)
    @test typeof(sm1DMeasInMemory) <: MDFv2InMemory

    for sm in (sm1DBrukerMeas,sm1DMeas,sm1DMeasInMemory)
      @info "Test $sm"

      @test size(measData(sm)) == (102, 3, 1, 67)

      @test rxNumSamplingPoints(sm) == 102
      @test rxNumFrequencies(sm) == 52
    end
  end
end

@testset "Testing Block Averaging" begin
  measPath = joinpath(tmpdir,"mdf","measurement_V2.mdf")
  avgMeasPath = joinpath(tmpdir,"mdf","avg.mdf")
  mdf = MPIFile(measPath)

	numF = acqNumFrames(mdf)
	numA = acqNumAverages(mdf)
	numAverages = 500

	# test exceptions, where no averaging should be done
	@test_throws DomainError blockAverage(avgMeasPath, mdf, 41)

	# average data
	blockAverage(avgMeasPath, mdf, numAverages)
	mdfavg = MPIFile(avgMeasPath)

	# test averaged data
	@test mean(measData(mdf),dims=4) ≈ measData(mdfavg)

	# test new acquisition parameters
	@test acqNumFrames(mdfavg) == div(numF,numAverages)
	@test acqNumAverages(mdfavg) == numA*numAverages

	# test new measurement parameters
	@test [any(measIsBGFrame(mdf))] == measIsBGFrame(mdfavg)
end


text:
using MPIMeasurements
using Plots
using Unitful
using UnitfulRecipes
using Dates
using UUIDs

#pyplot()
#plotly()
#inspectdr()
default(show = true)

ENV["JULIA_DEBUG"] = "all"

# Add test configurations to path
testConfigDir = normpath(string(@__DIR__), "../test/TestConfigs")
addConfigurationPath(testConfigDir)

scannerName_ = "TestSimpleSimulatedScanner"
protocolName_ = "SimulatedMagneticFieldMeasurement"
operator = "Jonas"
filename_ = joinpath("./field.h5")

protocol = Protocol(protocolName_, scannerName_)
filename(protocol, filename_)
runProtocol(protocol)


text:
# This file was generated, do not modify it. # hide
Random.seed!(1515)

SPCA = @pipeline(Standardizer(),
                 PCA(pratio=1-1e-4))

spca = machine(SPCA, X)
fit!(spca)
W = transform(spca, X)
names(W)


text:
using AWSCore
using AWSS3
using AWSLambda

JL_VERSION_BASE="0.6"
JL_VERSION_PATCH="2"
JL_VERSION="$JL_VERSION_BASE.$JL_VERSION_PATCH"

image_name = "octech/$(replace(basename(pwd()), "_", "")):$JL_VERSION"

lambda_name = basename(pwd())

source_bucket = "awslambda.jl.deploy.$(aws_account_number(aws_config()))"

base_zip = "$(lambda_name)_$(VERSION)_$(AWSLambda.aws_lamabda_jl_version).zip"

if length(ARGS) == 0 || ARGS[1] == "build"

    cp("../../src/AWSLambda.jl", "AWSLambda.jl"; remove_destination=true)

    for f in ["AWSLambdaWrapper.jl",
              "module_jl_lambda_eval.jl",
              "Dockerfile",
              "lambda_config.py",
              "lambda_function.py"]

        cp("../jl_lambda_eval/$f", f; remove_destination=true)
    end

    run(`docker build
            --build-arg JL_VERSION_BASE=$JL_VERSION_BASE
            --build-arg JL_VERSION_PATCH=$JL_VERSION_PATCH
             -t $image_name .`)
end

if length(ARGS) > 0 && ARGS[1] == "shell"
    run(`docker run --rm -it -v $(pwd()):/var/host $image_name bash`)
end

if length(ARGS) > 0 && ARGS[1] == "zip"
    rm(base_zip; force=true)
    cmd = `zip --symlinks -r -9 /var/host/$base_zip .`
    run(`docker run --rm -it -v $(pwd()):/var/host $image_name $cmd`)
end

if length(ARGS) > 0 && ARGS[1] == "deploy"
    s3_create_bucket(source_bucket)
    AWSCore.Services.s3("PUT", "/$source_bucket/$base_zip",
                       headers=Dict("x-amz-acl" => "public-read"),
                       Body=read(base_zip))
    AWSLambda.deploy_jl_lambda_eval(bucket = source_bucket, base_zip = base_zip)
end