text:
@testset "Water SPC/Fw test" begin
    qp = 1 # charge of a proton
    T = 298.16 # °K
    kb = 8.3144598e-3 # kJ/(K*mol)
    ϵOO = 0.1554253*4.184 # kJ
    σOO = 0.3165492 # nm
    ρ = 997/1.6747# Da/nm^3
    mO = 15.999 # Da
    mH = 1.00794 # Da
    mH2O = mO+2*mH
    N = 216
    L = (mH2O*N/ρ)^(1/3) # nm
    R = 0.9 # ≈3*σOO
    Rel = 0.49*L
    v_dev = sqrt(kb * T / mH2O) # nm/ps
    k_bond = 1059.162*4.184*1e2 # kJ/(mol*nm^2)
    k_angle = 75.90*4.184 # kJ/(mol*rad^2)
    rOH = 0.1012 # nm
    ∠HOH = 113.24 * pi / 180 # rad
    qH = 0.41 * qp
    qO = -0.84 * qp
    k_el = 138.935458 #
    τ = 0.5e-3 # ps

    jl_parameters = LennardJonesParameters(ϵOO, σOO, R)
    e_parameters = ElectrostaticParameters(k_el, Rel)
    spc_paramters = SPCFwParameters(rOH, ∠HOH, k_bond, k_angle)
    pbc = CubicPeriodicBoundaryConditions(L)

    @testset "Analyzing simulation result" begin
        t1 = 0τ
        t2 = 10τ

        r1 = SVector(L / 3, L / 3, 2 * L / 3)
        r2 = SVector(L / 3, 2 * L / 3, L / 3)
        r3 = SVector(2 * L / 3, L / 3, L / 3)
        v1 = SVector(0, 0, -v_dev)
        v2 = SVector(0, -v_dev, 0)
        v3 = SVector(-v_dev, 0, 0)
        p1 = MassBody(r1, v1, mH2O)
        p2 = MassBody(r2, v2, mH2O)
        p3 = MassBody(r3, v3, mH2O)

        bodies = [p1, p2, p3]
        water = WaterSPCFw(bodies, mH, mO, qH, qO,  jl_parameters, e_parameters, spc_paramters)
        simulation = NBodySimulation(water, (t1, t2), pbc, kb)

        result = run_simulation(simulation, VelocityVerlet(), dt=τ)

        cc = get_position(result, t2)
        ε = 0.01
        for i = 1:3
            indO = 3 * (i - 1) + 1
            @test (rOH - norm(cc[:,indO] - cc[:,indO + 1])) / rOH ≈ 0.0 atol = ε
            @test (rOH - norm(cc[:,indO] - cc[:,indO + 2])) / rOH ≈ 0.0 atol = ε
            ang = acos(dot(cc[:,indO] - cc[:,indO + 1], cc[:,indO] - cc[:,indO + 2]) / (norm(cc[:,indO] - cc[:,indO + 1]) * norm(cc[:,indO] - cc[:,indO + 2])))
            @test (ang - ∠HOH) / ∠HOH ≈ 0.0 atol = ε
        end

        e_tot_1 = total_energy(result, t1)
        e_tot_init = initial_energy(simulation)
        @test e_tot_1 == e_tot_init

        e_tot_2 = total_energy(result, t2)
        @test (e_tot_1 - e_tot_2)/e_tot_1 ≈ 0.0 atol = ε

        temperature_expected = (dot(v1, v1) + dot(v2, v2) + dot(v3, v3)) * (2 * mH + mO) / (kb * 21)
        temperature_1 = temperature(result, t1)
        @test (temperature_1 - temperature_expected)/temperature_expected ≈ 0.0 atol = ε

        count_plotting_data = 0
        for plotting_data ∈ result
            count_plotting_data += 1
            @test plotting_data[2] == result.solution.t[count_plotting_data]
        end
        @test count_plotting_data == length(result.solution.t)

        (ts, mean_square_displacement) = msd(result)
        @test mean_square_displacement[1] < mean_square_displacement[end]
    end

    @testset "Preparing data for thermostating" begin
        t1 = 0τ
        t2 = 10τ
        T0 = 275
        r1 = SVector(L / 3, L / 3, 2 * L / 3)
        r2 = SVector(L / 3, 2 * L / 3, L / 3)
        v1 = SVector(0, 0, -v_dev)
        v2 = SVector(0, -v_dev, 0)
        p1 = MassBody(r1, v1, mH2O)
        p2 = MassBody(r2, v2, mH2O)

        bodies = [p1, p2]
        water = WaterSPCFw(bodies, mH, mO, qH, qO,  jl_parameters, e_parameters, spc_paramters)

        thermostat = BerendsenThermostat(T0, 200τ)
        simulation = NBodySimulation(water, (t1, t2), pbc, kb)
        (ms_ber, kb_ber, n_ber, nc_ber, p_ber) = NBodySimulator.obtain_data_for_berendsen_thermostating(simulation)
        @test length(ms_ber) == n_ber
        @test kb_ber == kb
        @test n_ber == 6
        @test nc_ber == 4

        thermostat = NoseHooverThermostat(T0, 200τ)
        simulation = NBodySimulation(water, (t1, t2), pbc, kb)
        (ms_nh, kb_nh, n_nh, nc_nh, γind, p_nh) = NBodySimulator.obtain_data_for_nosehoover_thermostating(simulation)
        @test length(ms_nh) == n_nh
        @test kb_nh == kb
        @test n_nh == 6
        @test nc_nh == 4
        @test γind == 19
    end

    @testset "Water thermostating" begin
        ττ = 0.5e-4 # ps
        t1 = 0ττ
        t2 = 200ττ
        T0 = 275
        bodies = generate_bodies_in_cell_nodes(N, mH2O, v_dev, L)
        water = WaterSPCFw(bodies, mH, mO, qH, qO,  jl_parameters, e_parameters, spc_paramters);
        thermostat = LangevinThermostat(T0, 100)
        simulation = NBodySimulation(water, (t1, t2), pbc, thermostat, kb)
        result = run_simulation(simulation, EM(),  dt=τ)

        T2 = temperature(result, t2)
        ε = 1.0
        @test abs(T2 - T0) / T0 ≈ 0.0 atol = ε
    end

    @testset "Protein Data Bank file compilation" begin
        t1 = 0τ
        t2 = 2τ

        r1 = SVector(L / 3, L / 3, 2 * L / 3)
        r2 = SVector(L / 3, 2 * L / 3, L / 3)
        v1 = SVector(0, 0, -v_dev)
        v2 = SVector(0, -v_dev, 0)
        p1 = MassBody(r1, v1, mH2O)
        p2 = MassBody(r2, v2, mH2O)

        bodies = [p1, p2]
        water = WaterSPCFw(bodies, mH, mO, qH, qO,  jl_parameters, e_parameters, spc_paramters)
        simulation = NBodySimulation(water, (t1, t2), pbc, kb)
        result = run_simulation(simulation, VelocityVerlet(), dt=τ)

        io = IOBuffer()
        pdb_data = sprint(io -> NBodySimulator.write_pdb_data(io, result))
        splitted_data = split(pdb_data, '\n')

        hetatm_count = 0
        timestep_count = 0

        for s in splitted_data
            if length(s) >= 10 && s[1:10] == "REMARK 250"
                timestep_count += 1
            elseif length(s) >= 6 && s[1:6] == "HETATM"
                hetatm_count += 1
            end
        end

        molecule_number = hetatm_count / timestep_count / 3

        @test length(result.solution.t) == timestep_count
        @test molecule_number == length(bodies)
    end

    @testset "Load water molecules from a PDB file" begin

        pdb_data = IOBuffer("MODEL     1
REMARK 250 time=0.0000 picoseconds
HETATM    1  O   HOH     1      18.642  18.642   0.000  1.00  0.00
HETATM    2  H1  HOH     1      17.685  18.642   0.000  1.00  0.00
HETATM    3  H2  HOH     1      18.882  17.715   0.000  1.00  0.00
HETATM    4  O   HOH     2       0.000   0.000   3.107  1.00  0.00
HETATM    5  H1  HOH     2       0.957   0.000   3.107  1.00  0.00
HETATM    6  H2  HOH     2      -0.240   0.927   3.107  1.00  0.00
HETATM    7  O   HOH     3      18.642  18.642   6.214  1.00  0.00
HETATM    8  H1  HOH     3      17.685  18.642   6.214  1.00  0.00
HETATM    9  H2  HOH     3      18.882  17.715   6.214  1.00  0.00
ENDMDL
MODEL     2
REMARK 250 time=0.0005 picoseconds
HETATM    1  O   HOH     1      18.642  18.642  -0.000  1.00  0.00
HETATM    2  H1  HOH     1      17.685  18.642   0.000  1.00  0.00
HETATM    3  H2  HOH     1      18.882  17.715   0.000  1.00  0.00
HETATM    4  O   HOH     2      -0.000  -0.000   3.107  1.00  0.00
HETATM    5  H1  HOH     2       0.958  -0.000   3.107  1.00  0.00
HETATM    6  H2  HOH     2      -0.240   0.928   3.107  1.00  0.00
HETATM    7  O   HOH     3      18.642  18.642   6.214  1.00  0.00
HETATM    8  H1  HOH     3      17.685  18.642   6.214  1.00  0.00
HETATM    9  H2  HOH     3      18.881  17.715   6.214  1.00  0.00
ENDMDL")

        bodies = NBodySimulator.extract_from_pdb(pdb_data)
        water = WaterSPCFw(bodies, mH, mO, qH, qO,  jl_parameters, e_parameters, spc_paramters)
        t1 = 0τ
        t2 = 2τ
        simulation = NBodySimulation(water, (t1, t2), pbc, kb)
        result = run_simulation(simulation, VelocityVerlet(), dt=τ)

        ε = 0.00005

        @test bodies[1].O.r[1] - 1.8642 ≈ 0.0 atol = ε
        @test bodies[1].O.r[2] - 1.8642 ≈ 0.0 atol = ε
        @test bodies[1].O.r[3] - 0.0 ≈ 0.0 atol = ε

        @test bodies[2].H1.r[1] - 0.0957 ≈ 0.0 atol = ε
        @test bodies[2].H1.r[2] - 0.0 ≈ 0.0 atol = ε
        @test bodies[2].H1.r[3] - 0.3107 ≈ 0.0 atol = ε

        @test bodies[3].H2.r[1] - 1.8882 ≈ 0.0 atol = ε
        @test bodies[3].H2.r[2] - 1.7715 ≈ 0.0 atol = ε
        @test bodies[3].H2.r[3] - 0.6214 ≈ 0.0 atol = ε

        @test bodies[3].H2.v[1] - 0.0 ≈ 0.0 atol = ε
        @test bodies[3].H2.v[2] - 0.0 ≈ 0.0 atol = ε
        @test bodies[3].H2.v[3] - 0.0 ≈ 0.0 atol = ε
    end
end


text:
G = ssrand(2,3,4)

show_construction(G)

for ny = 1:3, nu = 1:3, nx = 1:3, proper=(true, false)
    G = ssrand(ny, nu, nx; proper)
    v = vec(G)
    G2 = vec2sys(v, ny, nu)
    @test G2 == G
end

ext: jl


text:
#
# Vadd
#
# Element wise addition of two vectors at a time in a chain (C=A+B; D=C+E; F=D+G)
# Asks the user to select a device at runtime
#
# History: Initial version based on vadd.c, written by Tim Mattson, June 2011
# Ported to C++ Wrapper API by Benedict Gaster, September 2011
# Updated to C++ Wrapper API v1.2 by Tom Deakin and Simon McIntosh-Smith, October 2012
# Ported to Python by Tom Deakin, July 2013
# Ported to Julia  by Jake Bolewski, Nov 2013

using OpenCL

# tolerance used in floating point comparisons
TOL = 1e-3

# length of vectors a, b, c
LENGTH = 1024

# Kernel: vadd
#
# To compute the elementwise sum c = a + b
#
# Input: a and b float vectors of length count
# Output c float vector of length count holding the sum a + b

kernelsource = "
__kernel void vadd(
    __global float* a,
    __global float* b,
    __global float* c,
    const unsigned int count)
{
    unsigned int i = get_global_id(0);
    if (i < count)
        c[i] = a[i] + b[i];
}
"

# create a compute context

# this selects the fastest opencl device available
# and creates a context and queue for using the
# the selected device
device, ctx, queue = cl.create_compute_context()

# create the compute program and build it
program = cl.Program(ctx, source=kernelsource) |> cl.build!

#create a, b, e, and g vectors and fill with random float values
#create empty vectors for c, d, and f
h_a = rand(Float32, LENGTH)
h_b = rand(Float32, LENGTH)
h_c = Vector{Float32}(undef, LENGTH)
h_d = Vector{Float32}(undef, LENGTH)
h_e = rand(Float32, LENGTH)
h_f = Vector{Float32}(undef, LENGTH)
h_g = rand(Float32, LENGTH)

# create the input (a,b,e,g) arrays in device memory and copy data from the host

# buffers can be passed memory flags:
# {:r = readonly, :w = writeonly, :rw = read_write (default)}

# buffers can also be passed flags for allocation:
# {:use (use host buffer), :alloc (alloc pinned memory), :copy (default)}

# Create the input (a, b, e, g) arrays in device memory and copy data from host
d_a = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=h_a)
d_b = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=h_b)
d_e = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=h_e)
d_g = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=h_g)
# Create the output (c, d, f) array in device memory
d_c = cl.Buffer(Float32, ctx, :w, LENGTH)
d_d = cl.Buffer(Float32, ctx, :w, LENGTH)
d_f = cl.Buffer(Float32, ctx, :w, LENGTH)

# create the kernel
vadd = cl.Kernel(program, "vadd")

# execute the kernel over the entire range of 1d input
# calling `queue` is blocking, it accepts the kernel, global / local work sizes,
# the the kernel's arguments.

# here we call the kernel with work size set to the number of elements and a local
# work size of nothing. This enables the opencl runtime to optimize the local size
# for simple kernels
queue(vadd, size(h_a), nothing, d_a, d_b, d_c, UInt32(LENGTH))

# an alternative syntax is to create an partial function to call
# by julia's getindex syntax for Kernel types.
# here the queue, global_size, and (optional) local_size are passed in which
# returns a partial function with these parameters set.
vadd[queue, size(h_e)](d_e, d_c, d_d, UInt32(LENGTH))
vadd[queue, size(h_g)](d_g, d_d, d_f, UInt32(LENGTH))

# copy back the results from the compute device
# copy!(queue, dst, src) follows same interface as julia's built in copy!
cl.copy!(queue, h_f, d_f)

# test the results
correct = 0
for i in 1:LENGTH
    tmp = h_a[i] + h_b[i] + h_e[i] + h_g[i]
    tmp -= h_f[i]
    if tmp^2 < TOL^2
        global correct += 1
    else
        println("tmp $tmp h_a $(h_a[i]) h_b $(h_b[i]) ",
                "h_e $(h_e[i]) h_g $(h_g[i]) h_f $(h_f[i])")
    end
end

# summarize results
println("3 vector adds to find F=A+B+E+G: $correct out of $LENGTH results were correct")


text:
"""
This function makes an elliptic solenoid and calculates the inductance.
    The function calculates the flux through the center of each loop's cross-section
        Lₛ = Φ/I (self inductance = flux per amp)
        Lₘ₂₁ = (Φ₂₁ ⋅ Φ₂₂/|Φ₂₂|) /I
            -mutual inductance from the 1ˢᵗ wire on 2ⁿᵈ is flux from 1 in the 2ⁿᵈ's cross section
            that is aligned with the flux from the second on itself. All per ampere
        Lₘ₂₁ = Lₘ₁₂
        Total L = ∑ᵢ∑ⱼ(Lᵢⱼ)
    The inputs are:
    N - Number of turns
    r₁ - one radius of the ellipse
    r₂ - Second radius of ellipse
    L - length of solenoid (axially)
    kwargs:
    MeasLayers - Number of concentric layers to calc field
                    The field is calculated over concentric ellipses (linearly spaced)
                    This is the number of concentric test point layers
    NLayers -  Number of test points per layer
    NPtsPath  - number of points that each wire loop will be discretized into,
    PlotOn deternimes if the coil will be plotted as it builds

"""
function eval_Solenoid_Induct(
    N, #Number of turns
    r₁, #Radius 1 of ellipse
    r₂,#Radius 2 of ellipse
    L; #Length of solenoid in meters
    DownSampleFac=1,PlotOn=false,NPtsPath=100,NLayers=20)

    if PlotOn
        pygui(true)
        gcf()
    end
    PointPath_SingleLoop =  MakeEllip(r₁, r₂;NPts=NPtsPath)

    LMat = zeros(N, N)
    if PlotOn
        scatter3D(PointPath_SingleLoop[:,1], PointPath_SingleLoop[:,2], PointPath_SingleLoop[:,3])
    end
    NewPP = MakeEllip(r₁,r₂;Center = [0, 0, 1/ N * L],NPts = NPtsPath)

    if PlotOn
        scatter3D(NewPP[:,1], NewPP[:,2], NewPP[:,3])
    end

    Mut, Self, SavedBSelfArr = Mutual_L_TwoLoops(PointPath_SingleLoop, NewPP;DownSampleFac=DownSampleFac,MeasLayers=NLayers,MinThreshold=1e-10,IncludeWireInduct=false,SaveΦ₁=true)
    LMat .+= Self .* OffDiagOnes(N, 1)
    LMat .+= Mut .* OffDiagOnes(N, 2)
    for i in 3:N
        NewPP = MakeEllip(r₁,r₂;Center = [0, 0, (i-1)/ N * L],NPts = NPtsPath)

        if PlotOn
            scatter3D(NewPP[:,1], NewPP[:,2], NewPP[:,3])
        end
        Mut = Mutual_L_TwoLoops(PointPath_SingleLoop, NewPP, SavedBSelfArr;DownSampleFac=DownSampleFac,MeasLayers=NLayers,MinThreshold=1e-10)
        println(Mut)

        LMat .+= Mut .* OffDiagOnes(N, i)

        println(i)
    end
    LMat[N,N] = LMat[1,1]
    return LMat, sum(LMat)

end

"""
This function makes an elliptic solenoid and calculates the field at the center of the solenoid
        Lₘ₂₁ = Lₘ₁₂
        Total L = ∑ᵢ∑ⱼ(Lᵢⱼ)
    The inputs are:
    N - Number of turns
    r₁ - one radius of the ellipse
    r₂ - Second radius of ellipse
    L - length of solenoid (axially)
    kwargs:

    NPts_Coil  - number of points that each wire loop will be discretized into,
"""
function EvalField_Centered(
    N, #Number of turns
    r₁, #Radius 1 of ellipse
    r₂,#Radius 2 of ellipse
    L; #Length of solenoid in meters
    NPts_Coil=100
    )

            FieldCentered = [
                BiotSav(
                    MakeEllip(
                                r₁, #Radius 1 of ellipse
                                r₂;#Radius 2 of ellipse
                                Center = [0, 0, (Wᵢ - 1) / N * L], #The solenoid's axis is in Z
                                NPts = NPts_Coil, #How discretized the windings are
                                ),
                        [0, 0, L/2];
                        MinThreshold = 0.001,
                        )[3]
                        for Wᵢ in 1:N
                        ]

    return sum(FieldCentered)
end


text:
function mygemm_pji!(C, A, B)
  ## Get the matrix sizes
  m, n = size(C)
  _, k = size(A)

  ## Check the matrix sizes
  @assert size(C) == (m, n)
  @assert size(A) == (m, k)
  @assert size(B) == (k, n)

  @inbounds begin
    for p = 1:k
      for j = 1:n
        for i = 1:m
          C[i, j] += A[i, p] * B[p, j]
        end
      end
    end
  end

  return C
end

function mygemm_jpi!(C, A, B)
  ## Get the matrix sizes
  m, n = size(C)
  _, k = size(A)

  ## Check the matrix sizes
  @assert size(C) == (m, n)
  @assert size(A) == (m, k)
  @assert size(B) == (k, n)

  @inbounds begin
    for j = 1:n
      for p = 1:k
        for i = 1:m
          C[i, j] += A[i, p] * B[p, j]
        end
      end
    end
  end

  return C
end

function mygemm_jip_pji!(C, A, B, mb, nb, kb)
  ## Get the matrix sizes
  m, n = size(C)
  _, k = size(A)

  ## Check the matrix sizes
  @assert size(C) == (m, n)
  @assert size(A) == (m, k)
  @assert size(B) == (k, n)

  @inbounds for J = 1:mb:m
    jrng = J:min(n, J+nb-1)
    for I = 1:nb:n
      irng = I:min(m, I+mb-1)
      for P = 1:kb:k
        prng = P:min(k, P+kb-1)
        @views mygemm_pji!(C[irng, jrng], A[irng, prng], B[prng, jrng])
      end
    end
  end
end

function mygemm_ji_pji!(C, A, B, mb, nb)
  ## Get the matrix sizes
  m, n = size(C)
  _, k = size(A)

  ## Check the matrix sizes
  @assert size(C) == (m, n)
  @assert size(A) == (m, k)
  @assert size(B) == (k, n)

  @inbounds for J = 1:mb:m
    J2 = min(n, J+nb-1)
    for I = 1:mb:m
      I2 = min(m, I+mb-1)
      @views mygemm_pji!(C[I:I2, J:J2], A[I:I2, :], B[:, J:J2])
    end
  end
end

# mr and nr need to be known at compile time!
# h/t to https://github.com/Sacha0/GemmDemo.jl for this realization
function mygemm_jpi_ji_pji!(C, A, B, mb, nb, kb, ::Val{mr}, ::Val{nr}) where {mr, nr}
  ## Get the matrix sizes
  m, n = size(C)
  _, k = size(A)

  ## Check the matrix sizes
  @assert size(C) == (m, n)
  @assert size(A) == (m, k)
  @assert size(B) == (k, n)

  @inbounds begin
    # Cache blocking
    for jb1 = 1:nb:n
      jb2 = min(n, jb1+nb-1)
      for pb1 = 1:kb:k
        pb2 = min(k, pb1+kb-1)
        for ib1 = 1:mb:m
          ib2 = min(n, ib1+mb-1)
          # Register Blocking
          for jr = jb1:nr:jb2
            for ir = ib1:mr:ib2
              #=
              @views mygemm_pji!(C[ir:ir+mr-1, jr:jr+nr-1],
                                 A[ir:ir+mr-1, pb1:pb2],
                                 B[pb1:pb2, jr:jr+nr-1])
              =#
              mygemm_pji!(C, A, B, ir, jr, pb1, pb2, Val(mr), Val(nr))
            end
          end
        end
      end
    end
  end
  return C
end

using SIMD
@inline function mygemm_pji!(C, A, B, ir, jr, p1, p2, ::Val{4}, ::Val{4})
  T = eltype(C)
  m = size(C, 1)
  ptrC = pointer(C, ir + (jr-1) * m)
  c_1 = vload(Vec{4, T}, ptrC + 0m*sizeof(T))
  c_2 = vload(Vec{4, T}, ptrC + 1m*sizeof(T))
  c_3 = vload(Vec{4, T}, ptrC + 2m*sizeof(T))
  c_4 = vload(Vec{4, T}, ptrC + 3m*sizeof(T))

  ptrA = pointer(A, ir)
  @inbounds for p = p1:p2
    a_p = vload(Vec{4, T}, ptrA + (p-1)*m*sizeof(T))

    b_pj = Vec{4,T}(B[p, jr+0])
    c_1 = fma(b_pj, a_p, c_1)

    b_pj = Vec{4,T}(B[p, jr+1])
    c_2 = fma(b_pj, a_p, c_2)

    b_pj = Vec{4,T}(B[p, jr+2])
    c_3 = fma(b_pj, a_p, c_3)

    b_pj = Vec{4,T}(B[p, jr+3])
    c_4 = fma(b_pj, a_p, c_4)
  end
  vstore(c_1, ptrC + 0m*sizeof(T))
  vstore(c_2, ptrC + 1m*sizeof(T))
  vstore(c_3, ptrC + 2m*sizeof(T))
  vstore(c_4, ptrC + 3m*sizeof(T))

  return C
end

@inline function mygemm_pji!(C, A, B, ir, jr, p1, p2, ::Val{mr}, ::Val{4}) where mr
  @assert mod(mr, 4) == 0 && mr <= 12

  T = eltype(C)
  m = size(C, 1)

  ptrC = pointer(C, ir + (jr-1) * m)
  mr > 0 && (c_1a = vload(Vec{4, T}, ptrC + (0m + 0)*sizeof(T)))
  mr > 4 && (c_1b = vload(Vec{4, T}, ptrC + (0m + 4)*sizeof(T)))
  mr > 8 && (c_1c = vload(Vec{4, T}, ptrC + (0m + 8)*sizeof(T)))
  mr > 0 && (c_2a = vload(Vec{4, T}, ptrC + (1m + 0)*sizeof(T)))
  mr > 4 && (c_2b = vload(Vec{4, T}, ptrC + (1m + 4)*sizeof(T)))
  mr > 8 && (c_2c = vload(Vec{4, T}, ptrC + (1m + 8)*sizeof(T)))
  mr > 0 && (c_3a = vload(Vec{4, T}, ptrC + (2m + 0)*sizeof(T)))
  mr > 4 && (c_3b = vload(Vec{4, T}, ptrC + (2m + 4)*sizeof(T)))
  mr > 8 && (c_3c = vload(Vec{4, T}, ptrC + (2m + 8)*sizeof(T)))
  mr > 0 && (c_4a = vload(Vec{4, T}, ptrC + (3m + 0)*sizeof(T)))
  mr > 4 && (c_4b = vload(Vec{4, T}, ptrC + (3m + 4)*sizeof(T)))
  mr > 8 && (c_4c = vload(Vec{4, T}, ptrC + (3m + 8)*sizeof(T)))

  ptrA = pointer(A, ir)
  @inbounds for p = p1:p2
    mr > 0 && (a_pa = vload(Vec{4, T}, ptrA + ((p-1)*m + 0)*sizeof(T)))
    mr > 4 && (a_pb = vload(Vec{4, T}, ptrA + ((p-1)*m + 4)*sizeof(T)))
    mr > 8 && (a_pc = vload(Vec{4, T}, ptrA + ((p-1)*m + 8)*sizeof(T)))

    b_pj = Vec{4,T}(B[p, jr+0])
    mr > 0 && (c_1a = fma(b_pj, a_pa, c_1a))
    mr > 4 && (c_1b = fma(b_pj, a_pb, c_1b))
    mr > 8 && (c_1c = fma(b_pj, a_pc, c_1c))

    b_pj = Vec{4,T}(B[p, jr+1])
    mr > 0 && (c_2a = fma(b_pj, a_pa, c_2a))
    mr > 4 && (c_2b = fma(b_pj, a_pb, c_2b))
    mr > 8 && (c_2c = fma(b_pj, a_pc, c_2c))

    b_pj = Vec{4,T}(B[p, jr+2])
    mr > 0 && (c_3a = fma(b_pj, a_pa, c_3a))
    mr > 4 && (c_3b = fma(b_pj, a_pb, c_3b))
    mr > 8 && (c_3c = fma(b_pj, a_pc, c_3c))

    b_pj = Vec{4,T}(B[p, jr+3])
    mr > 0 && (c_4a = fma(b_pj, a_pa, c_4a))
    mr > 4 && (c_4b = fma(b_pj, a_pb, c_4b))
    mr > 8 && (c_4c = fma(b_pj, a_pc, c_4c))
  end

  mr > 0 && (vstore(c_1a, ptrC + (0m + 0)*sizeof(T)))
  mr > 4 && (vstore(c_1b, ptrC + (0m + 4)*sizeof(T)))
  mr > 8 && (vstore(c_1c, ptrC + (0m + 8)*sizeof(T)))
  mr > 0 && (vstore(c_2a, ptrC + (1m + 0)*sizeof(T)))
  mr > 4 && (vstore(c_2b, ptrC + (1m + 4)*sizeof(T)))
  mr > 8 && (vstore(c_2c, ptrC + (1m + 8)*sizeof(T)))
  mr > 0 && (vstore(c_3a, ptrC + (2m + 0)*sizeof(T)))
  mr > 4 && (vstore(c_3b, ptrC + (2m + 4)*sizeof(T)))
  mr > 8 && (vstore(c_3c, ptrC + (2m + 8)*sizeof(T)))
  mr > 0 && (vstore(c_4a, ptrC + (3m + 0)*sizeof(T)))
  mr > 4 && (vstore(c_4b, ptrC + (3m + 4)*sizeof(T)))
  mr > 8 && (vstore(c_4c, ptrC + (3m + 8)*sizeof(T)))

  return C
end

function mygemm_jpi_ji_pji_packed!(C, A, B, mb, nb, kb, ::Val{mr}, ::Val{nr}) where {mr, nr}
  ## Get the matrix sizes
  m, n = size(C)
  _, k = size(A)

  ## Check the matrix sizes
  @assert size(C) == (m, n)
  @assert size(A) == (m, k)
  @assert size(B) == (k, n)

  Bt = similar(B, nr, kb, div(nb, nr))
  At = similar(A, mr, kb, div(mb, mr))

  @inbounds begin
    # Cache blocking
    for jb1 = 1:nb:n
      jb2 = min(n, jb1+nb-1)
      for pb1 = 1:kb:k
        pb2 = min(k, pb1+kb-1)
        for (jcol, j) in enumerate(jb1:nr:jb2)
          Bt[:, 1:pb2-pb1+1, jcol] = B[pb1:pb2, j:j+nr-1]'
        end
        for ib1 = 1:mb:m
          ib2 = min(n, ib1+mb-1)
          for (icol, i) in enumerate(ib1:mr:ib2)
            At[:, 1:pb2-pb1+1, icol] = A[i:i+mr-1, pb1:pb2]
          end
          # Register Blocking
          for (jcol, jr) in enumerate(jb1:nr:jb2)
            for (icol, ir) in enumerate(ib1:mr:ib2)
              mygemm_pji_packed!(C, At, Bt, ir, jr, icol, jcol, pb2-pb1+1, kb, Val(mr), Val(nr))
            end
          end
        end
      end
    end
  end
  return C
end

@inline function mygemm_pji_packed!(C, At, Bt, ir, jr, icol, jcol, pend, kb, ::Val{4}, ::Val{4})
  T = eltype(C)
  m = size(C, 1)
  ptrC = pointer(C, ir + (jr-1) * m)
  c_1 = vload(Vec{4, T}, ptrC + 0m*sizeof(T))
  c_2 = vload(Vec{4, T}, ptrC + 1m*sizeof(T))
  c_3 = vload(Vec{4, T}, ptrC + 2m*sizeof(T))
  c_4 = vload(Vec{4, T}, ptrC + 3m*sizeof(T))

  mr = 4
  ptrAt = pointer(At, 1 + (icol-1) * mr * kb)
  @inbounds for p = 1:pend
    a_p = vload(Vec{4, T}, ptrAt + (p-1)*mr*sizeof(T))

    b_pj = Vec{4,T}(Bt[1, p, jcol])
    c_1 = fma(b_pj, a_p, c_1)

    b_pj = Vec{4,T}(Bt[2, p, jcol])
    c_2 = fma(b_pj, a_p, c_2)

    b_pj = Vec{4,T}(Bt[3, p, jcol])
    c_3 = fma(b_pj, a_p, c_3)

    b_pj = Vec{4,T}(Bt[4, p, jcol])
    c_4 = fma(b_pj, a_p, c_4)
  end
  vstore(c_1, ptrC + 0m*sizeof(T))
  vstore(c_2, ptrC + 1m*sizeof(T))
  vstore(c_3, ptrC + 2m*sizeof(T))
  vstore(c_4, ptrC + 3m*sizeof(T))

  return C
end

@inline function mygemm_pji_packed!(C, At, Bt, ir, jr, icol, jcol, pend, kb, ::Val{mr}, ::Val{4}) where mr
  @assert mod(mr, 4) == 0 && mr <= 12

  T = eltype(C)
  m = size(C, 1)

  ptrC = pointer(C, ir + (jr-1) * m)
  mr > 0 && (c_1a = vload(Vec{4, T}, ptrC + (0m + 0)*sizeof(T)))
  mr > 4 && (c_1b = vload(Vec{4, T}, ptrC + (0m + 4)*sizeof(T)))
  mr > 8 && (c_1c = vload(Vec{4, T}, ptrC + (0m + 8)*sizeof(T)))
  mr > 0 && (c_2a = vload(Vec{4, T}, ptrC + (1m + 0)*sizeof(T)))
  mr > 4 && (c_2b = vload(Vec{4, T}, ptrC + (1m + 4)*sizeof(T)))
  mr > 8 && (c_2c = vload(Vec{4, T}, ptrC + (1m + 8)*sizeof(T)))
  mr > 0 && (c_3a = vload(Vec{4, T}, ptrC + (2m + 0)*sizeof(T)))
  mr > 4 && (c_3b = vload(Vec{4, T}, ptrC + (2m + 4)*sizeof(T)))
  mr > 8 && (c_3c = vload(Vec{4, T}, ptrC + (2m + 8)*sizeof(T)))
  mr > 0 && (c_4a = vload(Vec{4, T}, ptrC + (3m + 0)*sizeof(T)))
  mr > 4 && (c_4b = vload(Vec{4, T}, ptrC + (3m + 4)*sizeof(T)))
  mr > 8 && (c_4c = vload(Vec{4, T}, ptrC + (3m + 8)*sizeof(T)))

  ptrAt = pointer(At, 1 + (icol-1) * mr * kb)
  @inbounds for p = 1:pend
    mr > 0 && (a_pa = vload(Vec{4, T}, ptrAt + ((p-1)*mr + 0)*sizeof(T)))
    mr > 4 && (a_pb = vload(Vec{4, T}, ptrAt + ((p-1)*mr + 4)*sizeof(T)))
    mr > 8 && (a_pc = vload(Vec{4, T}, ptrAt + ((p-1)*mr + 8)*sizeof(T)))

    b_pj = Vec{4,T}(Bt[1, p, jcol])
    mr > 0 && (c_1a = fma(b_pj, a_pa, c_1a))
    mr > 4 && (c_1b = fma(b_pj, a_pb, c_1b))
    mr > 8 && (c_1c = fma(b_pj, a_pc, c_1c))

    b_pj = Vec{4,T}(Bt[2, p, jcol])
    mr > 0 && (c_2a = fma(b_pj, a_pa, c_2a))
    mr > 4 && (c_2b = fma(b_pj, a_pb, c_2b))
    mr > 8 && (c_2c = fma(b_pj, a_pc, c_2c))

    b_pj = Vec{4,T}(Bt[3, p, jcol])
    mr > 0 && (c_3a = fma(b_pj, a_pa, c_3a))
    mr > 4 && (c_3b = fma(b_pj, a_pb, c_3b))
    mr > 8 && (c_3c = fma(b_pj, a_pc, c_3c))

    b_pj = Vec{4,T}(Bt[4, p, jcol])
    mr > 0 && (c_4a = fma(b_pj, a_pa, c_4a))
    mr > 4 && (c_4b = fma(b_pj, a_pb, c_4b))
    mr > 8 && (c_4c = fma(b_pj, a_pc, c_4c))
  end

  mr > 0 && (vstore(c_1a, ptrC + (0m + 0)*sizeof(T)))
  mr > 4 && (vstore(c_1b, ptrC + (0m + 4)*sizeof(T)))
  mr > 8 && (vstore(c_1c, ptrC + (0m + 8)*sizeof(T)))
  mr > 0 && (vstore(c_2a, ptrC + (1m + 0)*sizeof(T)))
  mr > 4 && (vstore(c_2b, ptrC + (1m + 4)*sizeof(T)))
  mr > 8 && (vstore(c_2c, ptrC + (1m + 8)*sizeof(T)))
  mr > 0 && (vstore(c_3a, ptrC + (2m + 0)*sizeof(T)))
  mr > 4 && (vstore(c_3b, ptrC + (2m + 4)*sizeof(T)))
  mr > 8 && (vstore(c_3c, ptrC + (2m + 8)*sizeof(T)))
  mr > 0 && (vstore(c_4a, ptrC + (3m + 0)*sizeof(T)))
  mr > 4 && (vstore(c_4b, ptrC + (3m + 4)*sizeof(T)))
  mr > 8 && (vstore(c_4c, ptrC + (3m + 8)*sizeof(T)))

  return C
end


text:
using Docile
using Lexicon
using CompilerTools

index = Index()
ct_mods = Docile.Collector.submodules(CompilerTools)

for mod in ct_mods
    update!(index, save("../doc/api/$mod.md",mod))
end
save("../doc/api/api-index.md", index)


text:
using BinaryBuilder

name = "OpenLibm"
version = v"0.7.0"
sources = [
    "https://github.com/JuliaMath/openlibm/archive/v$(version).tar.gz" =>
    "1699f773198018b55b12631db9c1801fe3ed191e618a1ee1be743f4570ae06a3",
]

script = raw"""
# Enter the funzone
cd ${WORKSPACE}/srcdir/openlibm-*

# Install into output
flags=("prefix=${prefix}")

# Build ARCH from ${target}
flags+=("ARCH=${target%-*-*}")

# OpenLibm build system doesn't recognize our windows cross compilers properly
if [[ ${target} == *mingw* ]]; then
    flags+=("OS=WINNT")
fi

# Add `CC` override, since OpenLibm seems to think it knows best:
flags+=("CC=$CC")

# Build the library
make "${flags[@]}" -j${nproc}

# Install the library
make "${flags[@]}" install

install_license ${WORKSPACE}/srcdir/openlibm-*/LICENSE.md
"""

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

products = [
    LibraryProduct("libopenlibm", :libopenlibm)
]

dependencies = [
]

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


text:
__precompile__()
module ursa_functions

using Gallium
using Plots
using JuMP
using Clp
using Distributions

struct URSAState
    b::Real     #belief
    A::Int64    #allowed level of action
end

mutable struct URSA #<: MDP{URSAState, Int64}  #second parameter is the type of actions
    discount_factor::Float64
    Num_actions::Int64
    True_lambda::Int64
    Num_t_actions::Int64
    uncertain_factor::Float64
    N_steps::Int64
    N_iter::Int64
    URSA()=new()
    emit_prob_1::Array{Float64}#,Num_actions,Num_t_actions}
    emit_prob_0::Array{Float64}
    r_action_hostile::Array{Float64}
    r_action_civil::Array{Float64}
end

function URSA(df::Float64,Na::Int64,lam::Int64,Nta::Int64,u_f::Float64,n_step::Int64,n_iter::Int64)
    this=URSA()  #default discount factor, default three actions
    this.discount_factor=df
    this.Num_actions=Na
    this.True_lambda=lam
    this.Num_t_actions=Nta
    this.uncertain_factor=u_f
    this.N_steps=n_step
    this.N_iter=n_iter
    this.emit_prob_1=[0.4 0.6;0.45 0.55;0.7 0.3]
    this.emit_prob_0=[0.6 0.4;0.75 0.25;0.9 0.1]
    this.r_action_hostile=[0 0 0]
    this.r_action_civil=[-0.1 -0.3 -0.7]
    return this
end

mutable struct PBVI_solver
    N_B::Int64   #number of belief points
    ursa::URSA    #the ursa problem to be solved
    B::Array{Float64}   #set of belief point
    rho_alpha::Array{Array{Float64,1}}   #corresponding alpah vector for the belief dependent reward
    alpha::Array{Array{Float64,1}}
    Value::Array{Float64,1}
    opt_action::Array{Int64}
    p_opt::Array{Array{Float64,2}}
    PBVI_solver()=new()
end

function entropy_ursa(x::Float64)
    @assert (x>=0 && x<=1)
    if x==0 || x==1
        return 0
    else
        return x*log(x)+(1-x)*log(1-x)
    end
end

function PBVI_solver(N_B::Int64,ursa::URSA)
    this=PBVI_solver()
    this.N_B=N_B
    this.ursa=ursa
    db=0.998/(N_B-1)
    this.B=collect(0.001:db:0.999)
    slope=log.(this.B)-log.(1-this.B)
    entro=entropy_ursa.(this.B)
    #now solve for the rho_alpha vector
    this.rho_alpha=Array{Array{Float64,1}}(N_B)
    this.alpha=Array{Array{Float64,1}}(N_B)
    this.Value=Array{Float64,1}(N_B)
    this.p_opt=Array{Array{Float64,2}}(N_B,ursa.Num_actions)

    for k=1:N_B
        A=[-1 1;
        1-this.B[k] this.B[k]]
        b=[slope[k];entro[k]]
        this.rho_alpha[k]=A\b
        @assert abs(dot_alpha_b(this.rho_alpha[k],this.B[k])-entro[k])<10.0^(-6.0)
    end
    initialize_alpha!(this)
    return this
end
#return dot product between alpha and belief state b
dot_alpha_b(alpha::Array{Float64},b::Float64)=alpha[1]*(1-b)+alpha[2]*b

function initialize_alpha!(pbvi_solver::PBVI_solver)
    #pbvi_solver.alpha=Array{Array{Float64,1},1}(pbvi_solver.N_B)
    #pbvi_solver.Value=Array{Float64,1}(pbvi_solver.N_B)
    m_0=minimum(pbvi_solver.ursa.r_action_civil)+entropy_ursa(0.5)
    m_1=minimum(pbvi_solver.ursa.r_action_hostile)+entropy_ursa(0.5)
    for i=1:pbvi_solver.N_B
        pbvi_solver.alpha[i]=[m_0;m_1]/(1-pbvi_solver.ursa.discount_factor)
    end
    #=for i=1:pbvi_solver.N_B
        b=pbvi_solver.B[i]
        reward=-Inf
        max_action=0
        for j=1:pbvi_solver.ursa.Num_actions  #when we have escalation of action, here needs modification
            reward_j=get_reward(b,j,pbvi_solver)
            if reward_j>reward
                reward=reward_j
                max_action=j
            end
        end
        pbvi_solver.alpha[i]=pbvi_solver.rho_alpha[i]+[pbvi_solver.ursa.r_action_civil[max_action];pbvi_solver.ursa.r_action_hostile[max_action]]
        #initial alpha= rho_alpha+r_alpha
    end =#
end

function get_reward(b::Float64, a::Int64, pbvi_solver::PBVI_solver)
    temp=zeros(pbvi_solver.N_B)
    for i=1:pbvi_solver.N_B
        temp[i]=dot_alpha_b(pbvi_solver.rho_alpha[i],b)
    end
    reward_rho=maximum(temp)
    @assert reward_rho<=entropy_ursa(b)+10.0^(-6.0)
    reward_r=pbvi_solver.ursa.r_action_hostile[a]*b+pbvi_solver.ursa.r_action_civil[a]*(1-b)
    return reward_rho+reward_r
end

function solve!(pbvi_solver::PBVI_solver)
    N_action=pbvi_solver.ursa.Num_actions
    N_b=pbvi_solver.N_B
    N_ob=pbvi_solver.ursa.Num_t_actions
    gamma=pbvi_solver.ursa.discount_factor
    pbvi_solver.opt_action=Array{Int64}(N_b)
    #see PBVI paper section 3.1 pseudo code
    Gamma_a_star=Array{Array{Float64,1},2}(N_action,N_b)
    for i=1:N_action
        for j=1:N_b
            Gamma_a_star[i,j]=pbvi_solver.rho_alpha[j]+[pbvi_solver.ursa.r_action_civil[i];pbvi_solver.ursa.r_action_hostile[i]]
        end
    end

    for n_iter=1:pbvi_solver.ursa.N_iter
        Gamma_a_o=Array{Array{Float64,1},3}(N_action,N_ob,N_b)
        for i=1:N_action
            for j=1:N_ob
                for k=1:N_b
                    Gamma_a_o[i,j,k]=gamma*[pbvi_solver.ursa.emit_prob_0[i,j];pbvi_solver.ursa.emit_prob_1[i,j]].*pbvi_solver.alpha[k]
                end
            end
        end

        Gamma_a_b=Array{Array{Float64,1},2}(N_action,N_b)
        temp_alpha=Array{Array{Float64,1},3}(N_action,N_b,N_ob)
        for i=1:N_action
            for j=1:N_b
                Gamma_a_b[i,j]=Gamma_a_star[i,j]
                for k=1:N_ob
                    temp=get_argmax_alpha(Gamma_a_o[i,k,:],pbvi_solver.B[j])
                    @assert size(temp[1])==(2,)
                    Gamma_a_b[i,j]+=temp[1]
                    temp_alpha[i,j,k]=temp[1]
                end
            end
        end

        #pbvi_solver.Value=zeros(N_b)
        for i=1:N_b
            pbvi_solver.alpha[i],pbvi_solver.opt_action[i]=get_argmax_alpha(Gamma_a_b[:,i],pbvi_solver.B[i])
            pbvi_solver.Value[i]=dot_alpha_b(pbvi_solver.alpha[i],pbvi_solver.B[i])
        end
    end #end the iteration loop
    Plots.plot(pbvi_solver.B,pbvi_solver.Value,xlabel="Belief",ylabel="Value",label="Value function")
    Plots.plot!(pbvi_solver.B,pbvi_solver.opt_action, label="Optimal actions")
end

function get_argmax_alpha(array_alpha::Array{Array{Float64,1},1},b::Float64)
    max_dot=-Inf
    max_index=0
    for (i,alpha) in enumerate(array_alpha)
        dot=dot_alpha_b(alpha,b)
        if dot>max_dot
            max_dot=dot
            max_index=i
        end
    end
    return deepcopy(array_alpha[max_index]), max_index
end

function robust_solve!(pb::PBVI_solver)
    N_action=pb.ursa.Num_actions
    N_b=pb.N_B
    N_ob=pb.ursa.Num_t_actions
    gamma=pb.ursa.discount_factor

    for n_iter=1:pb.ursa.N_iter
        pb.opt_action=Array{Int64}(N_b)

        Lam_til=Array{Array{Float64,1},1}()
        for i=1:N_b
            Lam_b_til=Array{Array{Float64,1},1}()
            for j=1:pb.ursa.Num_actions
                pb.p_opt[i,j], alpha_z_opt= get_p_alpha(pb,pb.B[i],j)
                #line 12 of algorithm 3 robust DP backup
                alpha_opt=pb.rho_alpha[i]+[pb.ursa.r_action_civil[j];pb.ursa.r_action_hostile[j]]
                temp=zeros(2)
                for k=1:pb.ursa.Num_t_actions
                    temp+=pb.p_opt[i,j][k,:].*alpha_z_opt[k]
                end
                temp*=pb.ursa.discount_factor
                alpha_opt+=temp
                push!(Lam_b_til,deepcopy(alpha_opt))
            end
            @assert size(Lam_b_til)==(pb.ursa.Num_actions,)
            argmax_alpha, pb.opt_action[i]=get_argmax_alpha(Lam_b_til,pb.B[i])
            @assert size(argmax_alpha)==(2,)
            pb.Value[i]=dot_alpha_b(argmax_alpha,pb.B[i])
            push!(Lam_til, deepcopy(argmax_alpha))
        end
        pb.alpha=deepcopy(Lam_til)
    end
    Plots.plot(pb.B,pb.Value,xlabel="Belief",ylabel="Value",label="Robust Value function")
    Plots.plot!(pb.B,pb.opt_action, label="Optimal actions")
end

function get_p_alpha(pb::PBVI_solver,b::Float64,a::Int64)
    m=Model(solver=ClpSolver())
    @variable(m, U[1:pb.ursa.Num_t_actions])
    @variable(m, p[1:pb.ursa.Num_t_actions,1:2]) #last dim represents lambda, p's bound depends on a
    @variable(m, rhs_5[1:pb.N_B,1:pb.ursa.Num_t_actions])
    @objective(m, Min, sum(U))
    #first constraint in (5)
    for i=1:pb.N_B
        for j=1:pb.ursa.Num_t_actions
            @constraint(m, rhs_5[i,j]==(1-b)*p[j,1]*pb.alpha[i][1]+b*p[j,2]*pb.alpha[i][2])
            @constraint(m, U[j]>=rhs_5[i,j])
        end
    end
    #second constraint in (5)
    for i=1:pb.ursa.Num_t_actions
        @constraint(m, p[i,1]>=pb.ursa.emit_prob_0[a,i]*pb.ursa.uncertain_factor)
        @constraint(m, p[i,1]<=pb.ursa.emit_prob_0[a,i]/pb.ursa.uncertain_factor)
        @constraint(m, p[i,2]>=pb.ursa.emit_prob_1[a,i]*pb.ursa.uncertain_factor)
        @constraint(m, p[i,2]<=pb.ursa.emit_prob_1[a,i]/pb.ursa.uncertain_factor)
    end
    #normalization constraint
    @constraint(m, sum(p[:,1])==1.0)
    @constraint(m, sum(p[:,2])==1.0)
    status=solve(m)
    #now retrieve the needed variables
    U_opt=getvalue(U)
    alpha_opt=Array{Array{Float64,1},1}(pb.ursa.Num_t_actions)
    rhs_5_opt=getvalue(rhs_5)
    for j=1:pb.ursa.Num_t_actions
        for i=1:pb.N_B
            if rhs_5_opt[i,j]>=U_opt[j]-10.0^(-3.0)
                alpha_opt[j]=pb.alpha[i]  #retrieve the optimal alpha that make the constraint tight
            end
        end
    end
    return getvalue(p), alpha_opt
end

#pay attention!!! We are not actually solving for the finite horizon problem
#this solved value function is an approximate for the infiniy horizon problem
# we mix the use of N_iter and N_steps
function policy(b::Float64,pb::PBVI_solver,flag::String)
    diff=abs.(pb.B-b)
    nearest=indmin(diff) # index for the belief
    opt_action=pb.opt_action[nearest]
    if flag=="robust"
        return opt_action,pb.p_opt[nearest,opt_action]
    end
    if flag=="non robust"||flag=="chance"
        return pb.opt_action[nearest],[pb.ursa.emit_prob_0[opt_action,:] pb.ursa.emit_prob_1[opt_action,:]]
    end
end

function adversarial_p(b::Float64,pb_rob::PBVI_solver,a::Int64)
    #get the adversarial prob given belief, robust solved policy (to retrieve the p_opt) and an action
    diff=abs.(pb_rob.B-b)
    nearest=indmin(diff) # index for the belief
    #opt_action=pb.opt_action[nearest]
    return pb_rob.p_opt[nearest,a]
end

function belief_update(b::Float64, p::Array{Float64,2}, z::Int64)
    b_up=p[z,2]*b/(p[z,2]*b+p[z,1]*(1-b))
    return b_up
end

function interp_value(b::Float64, pb::PBVI_solver)
    #use linear interpolation to get the value at anywhere [0,1]
    @assert b<=1.0
    @assert b>=0
    if b<=pb.B[1]
        delta_b=pb.B[2]-pb.B[1]
        delta_value=pb.Value[2]-pb.Value[1]
        slope=delta_value/delta_b
        delta=b-pb.B[1]
        return pb.Value[1]+slope*delta
    end
    if b>=pb.B[end]
        delta_b=pb.B[end]-pb.B[end-1]
        delta_value=pb.Value[end]-pb.Value[end-1]
        slope=delta_value/delta_b
        delta=b-pb.B[end]
        return pb.Value[end]+slope*delta
    end
    for i=1:pb.N_B-1
        if b>=pb.B[i]&& b<=pb.B[i+1]  #b between i and i+1
            delta_b=pb.B[i+1]-pb.B[i]
            delta_value=pb.Value[i+1]-pb.Value[i]
            slope=delta_value/delta_b
            delta=b-pb.B[i]
            return pb.Value[i]+slope*delta
        end
    end
end

function chance_solve!(pb::PBVI_solver,epsilon::Float64,sigma2::Float64)
    #epsilon is the percentile
    #note that the get_chance_threshold only works for 2 observations
    N_action=pb.ursa.Num_actions
    N_b=pb.N_B
    N_ob=pb.ursa.Num_t_actions
    gamma=pb.ursa.discount_factor
    for n_iter=1:pb.ursa.N_iter
        pb.opt_action=Array{Int64}(N_b)
        for i=1:N_b
            value_a=zeros(pb.ursa.Num_actions)  #record the value of a at that belief point b
            for j=1:pb.ursa.Num_actions
                value_a[j]= get_chance_threshold(pb,pb.B[i],j,epsilon,sigma2)
            end
            pb.Value[i],pb.opt_action[i]=findmax(value_a)
        end
    end
    Plots.plot(pb.B,pb.Value,xlabel="Belief",ylabel="Value",label="Chance Value function")
    Plots.plot!(pb.B,pb.opt_action, label="Optimal actions")
end

function get_chance_threshold(pb::PBVI_solver,b::Float64,a::Int64,epsilon::Float64,sigma2::Float64)
    #solve a chance optimization
    #epsilon is the percentile, sigma2 is the variance
    reward=get_reward(b,a,pb)
    gamma=pb.ursa.discount_factor
    #check if this matrix is correct
    p=[pb.ursa.emit_prob_0[a,1] pb.ursa.emit_prob_1[a,1];
    pb.ursa.emit_prob_0[a,2] pb.ursa.emit_prob_1[a,2]]
    b1=belief_update(b,p,1)  #if receive obs1
    b2=belief_update(b,p,2)  #if receive obs2
    V1=interp_value(b1,pb)
    V2=interp_value(b2,pb)
    p11=pb.ursa.emit_prob_0[a,1]
    p21=pb.ursa.emit_prob_0[a,2]
    p12=pb.ursa.emit_prob_1[a,1]
    p22=pb.ursa.emit_prob_1[a,2]
    mean_lhs=(1-b)*(p11*V1+p21*V2)+b*(p12*V1+p22*V2)
    mean_lhs=gamma*mean_lhs
    var_lhs=gamma^2*sigma2*(b^2*(V1-V2)^2+(1-b)^2*(V1-V2)^2)
    y=mean_lhs-sqrt(var_lhs)*erfinv(1.0-epsilon)
    return reward+y
end

function ursa_target_action_sample(u::URSA,a::Int64)
    if u.True_lambda==1
        if !isprobvec(u.emit_prob_1[a,:])
            print(u.emit_prob_1[a,:])
        end
        @assert abs(sum(u.emit_prob_1[a,:])-1)<=10.0^(-5.0)
        dis=Categorical(u.emit_prob_1[a,:])
    elseif u.True_lambda==0
        if !isprobvec(u.emit_prob_0[a,:])
            print(u.emit_prob_0[a,:])
        end
        @assert abs(sum(u.emit_prob_0[a,:])-1)<=10.0^(-5.0)
        dis=Categorical(u.emit_prob_0[a,:])
    end
    sample_action=rand(dis)
    @assert sample_action>0
    return sample_action
end

function simulate_true_reward_ursa(i_b::Float64, pb::PBVI_solver, flag::String, flag_ursa::String, u_true::URSA)
    @assert flag=="robust"||flag=="non robust"||flag=="chance"
    @assert flag_ursa=="norminal"||flag_ursa=="random"||flag_ursa=="deter"
    N=pb.ursa.N_steps
    a=zeros(Int64, N)
    a_t=zeros(Int64, N)
    beliefs=zeros(N)
    beliefs[1]=i_b
    cumu_reward=0.0
    imme_reward=zeros(N)
    p_opt=Array{Array{Float64,2},1}(N)
    for i=1:N
        a[i], p_opt[i]=policy(beliefs[i],pb,flag)
        if flag_ursa=="random"
            u=sample_ursa(pb)
            a_t[i]=ursa_target_action_sample(u,a[i])
        elseif flag_ursa=="norminal"
            a_t[i]=ursa_target_action_sample(pb.ursa,a[i])
        elseif flag_ursa=="deter"
            a_t[i]=ursa_target_action_sample(u_true,a[i])
        end
        r=get_true_reward(beliefs[i], a[i], pb)
        imme_reward[i]=r
        cumu_reward+=r*pb.ursa.discount_factor^(i-1)
        if i<=N-1
            @assert a_t[i]>0
            beliefs[i+1]=belief_update(beliefs[i], p_opt[i], a_t[i])
        end
        #println("Step=",i,", belief=",round(beliefs[i],3),", action=",a[i],", t_action=",a_t[i],", i_reward=",round(r,3),", c_reward=",round(cumu_reward,3))
    end
    return cumu_reward
end

function simulate_true_reward_ursa_adversary(i_b::Float64, pb::PBVI_solver, flag::String, pb_rob::PBVI_solver)
    @assert flag=="robust"||flag=="non robust"||flag=="chance"
    N=pb.ursa.N_steps
    a=zeros(Int64, N)
    a_t=zeros(Int64, N)
    beliefs=zeros(N)
    beliefs[1]=i_b
    cumu_reward=0.0
    imme_reward=zeros(N)
    p_opt=Array{Array{Float64,2},1}(N)
    p_true=Array{Array{Float64,2},1}(N)
    for i=1:N
        a[i], p_opt[i]=policy(beliefs[i],pb,flag)
        p_true[i]=adversarial_p(beliefs[i],pb_rob,a[i])  #get an adversarial probability
        dis=Categorical(p_true[i][:,pb.ursa.True_lambda+1])
        a_t[i]=rand(dis)
        r=get_true_reward(beliefs[i], a[i], pb)
        imme_reward[i]=r
        cumu_reward+=r*pb.ursa.discount_factor^(i-1)
        if i<=N-1
            @assert a_t[i]>0
            beliefs[i+1]=belief_update(beliefs[i], p_opt[i], a_t[i])
        end
        #println("Step=",i,", belief=",round(beliefs[i],3),", action=",a[i],", t_action=",a_t[i],", i_reward=",round(r,3),", c_reward=",round(cumu_reward,3))
    end
    return cumu_reward
end

function get_true_reward(b::Float64, a::Int64, pbvi_solver::PBVI_solver)
    if pbvi_solver.ursa.True_lambda==1
        reward_rho=log(b)
        reward_r=pbvi_solver.ursa.r_action_hostile[a]
    elseif pbvi_solver.ursa.True_lambda==0
        reward_rho=log(1-b)
        reward_r=pbvi_solver.ursa.r_action_civil[a]
    end
    return reward_rho+reward_r
end

function sample_ursa(pb::PBVI_solver)
    u=URSA(pb.ursa.discount_factor,pb.ursa.Num_actions,pb.ursa.True_lambda,pb.ursa.Num_t_actions,pb.ursa.uncertain_factor,pb.ursa.N_steps,pb.ursa.N_iter)
    #now u already has emit and rewards
    dis=Categorical([0.5;0.5])
    for i=1:pb.ursa.Num_actions
        rand_ind=rand(dis)
        range=u.emit_prob_0[i,rand_ind]/u.uncertain_factor-u.emit_prob_0[i,rand_ind]*u.uncertain_factor
        u.emit_prob_0[i,rand_ind]*=u.uncertain_factor
        u.emit_prob_0[i,rand_ind]+=rand()*range
        u.emit_prob_0[i,rand_ind]=min(u.emit_prob_0[i,rand_ind],1.0)
        if rand_ind==1
            u.emit_prob_0[i,2]=1-u.emit_prob_0[i,1]
        elseif rand_ind==2
            u.emit_prob_0[i,1]=1-u.emit_prob_0[i,2]
        end
        @assert abs(u.emit_prob_0[i,1]+u.emit_prob_0[i,2]-1)<=10.0^(-5.0)
    end

    for i=1:pb.ursa.Num_actions
        rand_ind=rand(dis)
        range=u.emit_prob_1[i,rand_ind]/u.uncertain_factor-u.emit_prob_1[i,rand_ind]*u.uncertain_factor
        u.emit_prob_1[i,rand_ind]*=u.uncertain_factor
        u.emit_prob_1[i,rand_ind]+=rand()*range
        u.emit_prob_1[i,rand_ind]=min(u.emit_prob_1[i,rand_ind],1.0)
        if rand_ind==1
            u.emit_prob_1[i,2]=1-u.emit_prob_1[i,1]
        elseif rand_ind==2
            u.emit_prob_1[i,1]=1-u.emit_prob_1[i,2]
        end
        @assert abs(u.emit_prob_1[i,1]+u.emit_prob_1[i,2]-1)<=10.0^(-5.0)
    end
    return u
end

export URSA, URSAState, PBVI_solver, adversarial_p, belief_update, dot_alpha_b, entropy_ursa, get_argmax_alpha
export get_p_alpha, get_reward, get_true_reward, initialize_alpha!, policy, robust_solve!, sample_ursa
export simulate_true_reward_ursa, simulate_true_reward_ursa_adversary, solve!, ursa_target_action_sample
export get_chance_threshold, interp_value, chance_solve!
end


text:
using RegistratorServerless
using Test

@testset "RegistratorServerless.jl" begin
    @testset "parse-trigger-issue.jl" begin
        @test RegistratorServerless.parse_trigger_issue_title("register MyOrganization/MyPackage.jl") == Dict(:repo => "MyOrganization/MyPackage.jl")
        @test RegistratorServerless.parse_trigger_issue_title("register MyOrganization/MyPackage.jl branch=my/custom-branch_name") == Dict(:repo => "MyOrganization/MyPackage.jl", :branch => "my/custom-branch_name")
    end
end


text:
set_bigfloat_precision(128)
function c(x)
	if x == 2
		return BigFloat(0.0)
	else
		return sqrt( (BigFloat(1.0) + c(x-1)) / BigFloat(2.0 ) )
	end
end

function s(x)
	if x == 2
		return BigFloat(1.0)
	else
		return sqrt( (BigFloat(1.0) - c(x-1)) / BigFloat(2.0 ))
	end
end

P(x) = BigFloat( 2.0) ^ (x - BigFloat(1.0)) * s(x)

data = 2
for i = 2:2*128
	@printf("%i: %1e\n",i,P(data))
	data = data +1
end

max_stars_repo_name: Magikis/Uniwersity

max_issues_repo_name: Magikis/Uniwersity

max_forks_repo_name: Magikis/Uniwersity


text:
using CacheVariables, BSON, Test

## Add data directory, define data file path
dirpath = joinpath(@__DIR__, "data")
isdir(dirpath) && rm(dirpath; recursive=true)
mkdir(dirpath)
path = joinpath(dirpath, "test.bson")

## Test save functionality
@testset "Save" begin
    # 1. verify log message
    @test_logs (:info, "Saving to $path\nx\ny\nz") (@cache path begin
        x = collect(1:3)
        y = 4
        z = "test"
	"final output"
    end)

    rm(path)
    out = @cache path begin
        x = collect(1:3)
        y = 4
        z = "test"
	"final output"
    end

    # 3. did variables enter workspace correctly?
    @test x == [1, 2, 3]
    @test y == 4
    @test z == "test"
    @test out == "final output"
end

## Test load functionality
@testset "Load" begin
    # 1. set all variables to nothing
    x = nothing
    y = nothing
    z = nothing
    out = nothing

    # 2. file exists: load variables from it
    # verify log message
    @test_logs (:info, "Loading from $path\nx\ny\nz") (@cache path begin
        x = collect(1:3)
        y = 4
        z = "test"
	"final output"
    end)

    # load variables
    out = @cache path begin
        x = collect(1:3)
        y = 4
        z = "test"
	"final output"
    end

    # 3. did variables enter workspace correctly?
    @test x == [1, 2, 3]
    @test y == 4
    @test z == "test"
    @test out == "final output"
end

## Test overwrite behavior
@testset "Overwrite" begin
    # 1. change file contents
    bson(path; x=nothing, y=nothing, z=nothing, ans=nothing)

    # 2. add `true` to @cache call to overwrite
    # validate log message
    @test_logs (:info, "Overwriting $path\nx\ny\nz") (@cache path begin
        x = collect(1:3)
        y = 4
        z = "test"
	"final output"
    end true)

    # overwrite data file
    overwrite = true
    @cache path begin
        x = collect(1:3)
        y = 4
        z = "test"
	"final output"
    end overwrite

    # 3. check that correct data was overwritten
    x = nothing
    y = nothing
    z = nothing
    out = nothing
    out = @cache path begin
        x = collect(1:3)
        y = 4
        z = "test"
	"final output"
    end

    @test x == [1, 2, 3]
    @test y == 4
    @test z == "test"
    @test out == "final output"
end

## Clean up
rm(dirpath; recursive=true)


text:
# Internal helpers for defining the `add!` field of an `InplaceableThunk`

_update!(x, y) = x + y
_update!(x::Array{T,N}, y::AbstractArray{T,N}) where {T,N} = x .+= y

_update!(x, ::Zero) = x
_update!(::Zero, y) = y
_update!(::Zero, ::Zero) = Zero()

function _update!(x::NamedTuple, y, p::Symbol)
    y = extern(y)
    yp = getproperty(y, p)
    xp = getproperty(x, p)
    new_xp = _update!(xp, yp)
    new = NamedTuple{(p,)}((new_xp,))
    return merge(x, new)
end

"""
    _checked_rrule

like `rrule` but throws an error if the `rrule` is not defined.
Rather than returning `nothing`
"""
function _checked_rrule(f, args...; kwargs...)
    r = rrule(f, args...; kwargs...)
    r isa Nothing && _throw_checked_rrule_error(f, args...; kwargs...)
    return r
end

@noinline function _throw_checked_rrule_error(f, args...; kwargs...)
    io = IOBuffer()
    print(io, "can't differentiate `", f, '(')
    join(io, map(arg->string("::", typeof(arg)), args), ", ")
    if !isempty(kwargs)
        print(io, ";")
        join(io, map(((k, v),)->string(k, "=", v), kwargs), ", ")
    end
    print(io, ")`; no matching `rrule` is defined")
    throw(ArgumentError(String(take!(io))))
end


text:
macro size(z, s); esc(:(@assert (size($z) == $s) string(summary($z),!=,$s))); end

macro sizes(s); quote
    Ex = size(s.srcembed.w,1)
    Vx = size(s.srcembed.w,2)
    Ey = size(s.tgtembed.w,1)
    Vy = size(s.tgtembed.w,2)
    Dx = s.encoder.direction + 1
    Hx = Int(s.encoder.hiddenSize)
    Hy = Int(s.decoder.hiddenSize)
    Lx = Int(s.encoder.numLayers)
    Ly = Int(s.decoder.numLayers)
    @assert Hx == Hy && Ly == Lx*Dx
end |> esc; end

max_stars_repo_name: denizyuret/research

max_issues_repo_name: denizyuret/research

max_forks_repo_name: denizyuret/research


text:
using Test
using EFTfitter
using BAT
using Distributions
using IntervalSets

Test.@testset "Package EFTfitter" begin
    include("test_datatypes.jl")
    include("test_utils.jl")
    include("test_inputs/test_inputs.jl")
    include("test_ranking.jl")
end


text:
# ---
# title: 559. Maximum Depth of N-ary Tree
# id: problem559
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Tree, Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/maximum-depth-of-n-ary-tree/description/>
# hidden: true
# ---
#
# Given a n-ary tree, find its maximum depth.
#
# The maximum depth is the number of nodes along the longest path from the root
# node down to the farthest leaf node.
#
# _Nary-Tree input serialization  is represented in their level order traversal,
# each group of children is separated by the null value (See examples)._
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2018/10/12/narytreeexample.png)
#
#
#
#     Input: root = [1,null,3,2,4,null,5,6]
#     Output: 3
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2019/11/08/sample_4_964.png)
#
#
#
#     Input: root = [1,null,2,3,4,5,null,null,6,7,null,8,null,9,10,null,null,11,null,12,null,13,null,null,14]
#     Output: 5
#
#
#
#
# **Constraints:**
#
#   * The depth of the n-ary tree is less than or equal to `1000`.
#   * The total number of nodes is between `[0, 104]`.
#
#
## @lc code=start
using LeetCode

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


text:
## Token interface to a container.  A token is the address
## of an item in a container.  The token has two parts: the
## container and the item's address.  The address is of type
## AbstractSemiToken.

module Tokens

abstract AbstractSemiToken

immutable IntSemiToken <: AbstractSemiToken
    address::Int
end

end


text:
# ---
# title: 1425. Constrained Subsequence Sum
# id: problem1425
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/constrained-subsequence-sum/description/>
# hidden: true
# ---
#
# Given an integer array `nums` and an integer `k`, return the maximum sum of a
# **non-empty** subsequence of that array such that for every two
# **consecutive** integers in the subsequence, `nums[i]` and `nums[j]`, where `i
# < j`, the condition `j - i <= k` is satisfied.
#
# A  _subsequence_  of an array is obtained by deleting some number of elements
# (can be zero) from the array, leaving the remaining elements in their original
# order.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [10,2,-10,5,20], k = 2
#     Output: 37
#     Explanation: The subsequence is [10, 2, 5, 20].
#
#
# **Example 2:**
#
#
#
#     Input: nums = [-1,-2,-3], k = 1
#     Output: -1
#     Explanation: The subsequence must be non-empty, so we choose the largest number.
#
#
# **Example 3:**
#
#
#
#     Input: nums = [10,-2,-10,-5,20], k = 2
#     Output: 23
#     Explanation: The subsequence is [10, -2, -5, 20].
#
#
#
#
# **Constraints:**
#
#   * `1 <= k <= nums.length <= 10^5`
#   * `-10^4 <= nums[i] <= 10^4`
#
#
## @lc code=start
using LeetCode

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


text:
# module TypedExpressions

# export econvert
# export AbstractExpr, ECall, EArgs, EAssignment, ETuple, EKeyword

abstract type AbstractExpr end

isa_expr(ex::Expr, ::Type{<:AbstractExpr}) = false
isa_expr(ex::T, ::Type{K}) where {T<:AbstractExpr, K<:AbstractExpr} = T == K

args(ex::T) where T<:AbstractExpr = ex.args
args(ex::Expr) = ex.args

""" ECall """
struct ECall <: AbstractExpr
    fun::Any
    args::Vector{Any}
end
ECall(ex::Expr) = ECall(ex.args[1], ex.args[2:end])
ECall(ex::ECall) = ex
isa_expr(ex::Expr, ::Type{ECall}) = ex.head == :call

""" ETuple """
struct ETuple <: AbstractExpr
    args
end
ETuple(ex::Expr) = ETuple(ex.args)
ETuple(ex::ETuple) = ex
isa_expr(ex::Expr, ::Type{ETuple}) = ex.head == :tuple

""" EAssignment """
struct EAssignment <: AbstractExpr
    lhs
    rhs
end
EAssignment(ex::Expr) = EAssignment(ex.args[1], ex.args[2])
EAssignment(ex::EAssignment) = ex
isa_expr(ex::Expr, ::Type{EAssignment}) = ex.head == :(=)
args(ex::EAssignment) = [ex.lhs; ex.rhs]

""" EUsing """
struct EUsing <: AbstractExpr
    args
end
EUsing(ex::Expr) = EUsing(ex.args)
EUsing(ex::EUsing) = ex
isa_expr(ex::Expr, ::Type{EUsing}) = ex.head == :using

""" EKeyword """
struct EKeyword <: AbstractExpr
    lhs
    rhs
end
EKeyword(ex::Expr) = EKeyword(ex.args[1], ex.args[2])
EKeyword(ex::EKeyword) = ex
isa_expr(ex::Expr, ::Type{EKeyword}) = ex.head == :kw
args(ex::EKeyword) = [ex.lhs; ex.rhs]

""" ERef """
struct ERef <: AbstractExpr
    lhs
    rhs
end
ERef(ex::Expr) = ERef(ex.args[1], ex.args[2])
ERef(ex::ERef) = ex
isa_expr(ex::Expr, ::Type{ERef}) = (ex.head == :ref) && (length(ex.args) == 2)
args(ex::ERef) = [ex.lhs; ex.rhs]

""" EDot """
struct EDot <: AbstractExpr
    lhs
    rhs
end
EDot(ex::Expr) = EDot(ex.args[1], ex.args[2])
EDot(ex::ERef) = ex
isa_expr(ex::Expr, ::Type{EDot}) = (ex.head == :.) && (length(ex.args) == 2)
args(ex::EDot) = [ex.lhs; ex.rhs]

#

const EExpr = [ECall, ETuple, ERef, EDot, EAssignment, EUsing, EKeyword]#order matters (pobably a bug)
const EArgs = Union{Expr, ECall, ETuple, EUsing}
const EHanded = Union{EAssignment, EKeyword, ERef, EDot}

# I should maybe loop at each node instead of at the top
function econvert(ex)
    for E in EExpr
        ex = econvert(E, ex)
    end
    ex
end

econvert(::Type{<:AbstractExpr}, x) = x

function econvert(::Type{T}, ex::EArgs) where T <: AbstractExpr
    ex = econvert(T, ex, ex.args)
    isa_expr(ex,T) ? T(ex) : ex
end
function econvert(::Type{T}, ex, args::Vector) where T <: AbstractExpr
    for i in 1:length(ex.args)
        ex.args[i] = econvert(T, ex.args[i])
    end
    ex
end
econvert(::Type{T}, ex, args) where T <: AbstractExpr = ex

function econvert(::Type{T}, ex::EHanded) where T <: AbstractExpr
    ex = econvert(T, ex, ex.rhs)
    ex = econvert(T, ex, ex.lhs)
    isa_expr(ex,T) ? T(ex) : ex
end

#end


text:
struct Target{W,T}
    v::NTuple{W,Core.VecElement{T}}
end
@inline Target() = Target(LoopVectorization.SIMDPirates.extract_data(LoopVectorization.vbroadcast(LoopVectorization.pick_vector_width_val(Float64), 0.0)))
@inline Target(::Type{T}) where {T} = Target(LoopVectorization.SIMDPirates.extract_data(LoopVectorization.vbroadcast(LoopVectorization.pick_vector_width_val(T), zero(T))))
@inline Base.:(+)(t::Target, s) = Target(LoopVectorization.addscalar(t.v, s))
@inline Base.:(+)(s, t::Target) = Target(LoopVectorization.addscalar(t.v, s))
@inline Base.:(+)(t::Target{WT,T}, v::LoopVectorization.SVec{WV,T}) where {WT,WV,T} = Target(LoopVectorization.vadd(t.v, LoopVectorization.SIMDPirates.extract_data(v)))
@inline Base.:(+)(t::Target{WT,T}, v::NTuple{WV,Core.VecElement{T}}) where {WT,WV,T} = Target(LoopVectorization.vadd(t.v, v))
@inline Base.:(+)(v::LoopVectorization.SVec{WV,T}, t::Target{WT,T}) where {WT,WV,T} = Target(LoopVectorization.vadd(t.v, LoopVectorization.SIMDPirates.extract_data(v)))
@inline Base.:(+)(v::NTuple{WV,Core.VecElement{T}}, t::Target{WT,T}) where {WT,WV,T} = Target(LoopVectorization.vadd(t.v, v))
@inline Base.:(+)(t₁::Target, t₂::Target) = Target(LoopVectorization.vadd(t₁.v, t₂.v))
@inline Base.sum(t::Target) = LoopVectorization.vsum(t.v)
@inline VectorizationBase.extract_data(t::Target) = t.v


text:
using DataStructures

function getAllAdaptersInBag(input)
    allAdapters::SortedDict{Int, Vector{Int}} = SortedDict{Int, Vector{Int}}()
    for line in readlines(input)
        allAdapters[parse(Int,line)] = []
    end
    return allAdapters
end

function addSocketAndDevice(allAdaptersInBag)
    allAdapters = allAdaptersInBag
    allAdapters[0] = []
    highestRating = findmax(collect(keys(allAdapters)))[1]
    allAdapters[highestRating + 3] = []
    return allAdapters
end

function createTree(allAdapters)
    allkeys = sort(collect(keys(allAdapters)))
    for i in length(allkeys):-1:1
        key = allkeys[i]
        kids = allAdapters[key]
        j = i-1
        while (j >= 1) && (key - allkeys[j] <= 3)
            push!(kids, allkeys[j])
            j -= 1
        end
        #println(key, " ", kids)
    end
    return allAdapters, last(allkeys)
end

function permutationsForRange(range)
    diff = range[1] - range[2]
    diff == 1 && return 1
    diff == 2 && return 2
    diff == 3 && return 4
    diff == 4 && return 7
    diff == 5 && return 13
    diff > 5 && return error("Haven't figured it out yet...")
end

function findNextConsecutiveRangeAndCalcuatePermutations(startIndex, allKeys, allAdapters)
    firstAdapterWithKids = 0
    lastAdapterWithKids = 0
    nextStart = startIndex
    permutations = 1
    #println("startIndex =" , startIndex)
    for i in startIndex:-1:1
        rating = allKeys[i]
        kids = allAdapters[rating]
        println("i, adapter:", i, " ", rating, " ", kids)

        if length(kids) > 1
            if firstAdapterWithKids == 0 # found first
                firstAdapterWithKids = rating
                nextStart -= 1
                #println("found first")
            else # found possibly last
                nextStart -= 1
                #println("found possibly last")
            end
        else
            if (firstAdapterWithKids == 0) # thank you, next
                nextStart -= 1
                #println("Thank you, next")
                break
            else # previus was last
                lastAdapterWithKids = allKeys[i+1]
                println("lastawk=", lastAdapterWithKids)
                max = firstAdapterWithKids
                lastKids = allAdapters[lastAdapterWithKids]
                println("lastkids=", lastKids)
                min = findmin(lastKids)[1]
                range = (max,min)
                permutations = permutationsForRange(range)
                println(range, " ", permutations)
                nextStart -= 1
                break
            end
        end
    end
    #println("permutations=",permutations )
    return nextStart, permutations
end

function solveProblemB()
    allAdaptersInBag = getAllAdaptersInBag("data.txt")
    allAdapters = addSocketAndDevice(allAdaptersInBag)
    adapterTree, rootRating = createTree(allAdapters)

    allKeys = sort(collect(keys(allAdapters)))
    totalPermutations = 1
    index = findmax(allKeys)[2]

    while index > 1
        nextIndex, permutations = findNextConsecutiveRangeAndCalcuatePermutations(index, allKeys, allAdapters)
        totalPermutations *= permutations
        index = nextIndex
    end
    println("Arrangements=", totalPermutations)
end

solveProblemB()


text:
"""
    MyAwesomePackage

Here you are, you found an awesome package :-)

"""
module MyAwesomePackage

export plusTwo

"""
    plusTwo(x)

Sum the numeric "2" to whatever it receives as input

A more detailed explanation can go here, altought I guess it is not needed in this case

# Arguments
* `x`: The amount to which we want to add 2

# Notes
* Notes can go here

# Examples
```julia
julia> five = plusTwo(3)
5
```
"""
plusTwo(x) = return x+2

end # module

max_stars_repo_name: sylvaticus/MyAwesomePackage

max_issues_repo_name: sylvaticus/MyAwesomePackage

max_forks_repo_name: sylvaticus/MyAwesomePackage


text:
part1(input) = sum([1401, 1191, 1154, 1034, 950][c - UInt8('0')] for c in read(input)[1:2:end])

part2(input) = sum([6206821033, 5617089148, 5217223242, 4726100874, 4368232009][c - UInt8('0')] for c in read(input)[1:2:end])


text:
# DEPRECATED Please use https://github.com/JuliaPOMDP/ParticleFilters.jl instead
# maintained by @ebalaban
# Particle-based belief and state distribution types

"""
DEPRECATED: Please use https://github.com/JuliaPOMDP/ParticleFilters.jl instead

Belief particle type that contains a state and its probability
"""
mutable struct Particle{T} # should be immutable?
    state::T # particle state
    weight::Float64 # particle prob
end

"""
DEPRECATED: Please use https://github.com/JuliaPOMDP/ParticleFilters.jl instead

ParticleBelief type. Fields:

    -particles: a vector particles (state and weight) in the belief
    -probs_dict: dictionary maps states to weights (used in pdf only)
    -probs_arr: weights of each particle in particles (for convenience)
    -keep_dict: boolean flag for keeping track of dictionary
"""
mutable struct ParticleBelief{T}
    particles::Vector{Particle{T}} # array of particles
    probs_dict::Dict{T, Float64} # dict state => prob
    probs_arr::Vector{Float64} # particle weights
    keep_dict::Bool # for more efficient pdf
end
function ParticleBelief{S}(particles::Vector{Particle{S}})
    warn("POMDPToolbox.ParticleBelief is deprecated, please use the beliefs from https://github.com/JuliaPOMDP/ParticleFilters.jl instead.")
    return ParticleBelief(particles, Dict{S,Float64}(), Float64[p.weight for p in particles], false)
end
const ParticleDistribution{T} = ParticleBelief{T}

function pdf{S}(b::ParticleBelief{S}, s::S)
    if !b.keep_dict
        # search through the particle array (very slow)
        w = 0.0
        for (i, p) in enumerate(b.particles)
            if p.state == s
                w += p.weight
            end
        end
        return w
    else
        # tracking probs in dict (more efficient)
        haskey(b.probs_dict, s) ? (return b.probs_dict[s]) : (return 0.0)
    end
    return 0.0
end

function rand{S}(rng::AbstractRNG, b::ParticleBelief{S})
    cat = Weights(b.probs_arr)
    k = sample(rng, cat)
    s = b.particles[k].state
    return s
end
rand{S}(rng::AbstractRNG, b::ParticleBelief{S}, s::S) = rand(rng, b)

mean(b::ParticleBelief) = sum(p.weight*p.state for p in b.particles)/sum(p.weight for p in b.particles)

function mode{T}(b::ParticleBelief{T}) # don't know if this is efficient
    d = Dict{T, Float64}()
    best_weight = first(b.particles).weight
    most_likely = first(b.particles).state
    for p in b.particles
        if haskey(d, p.state)
            d[p.state] += p.weight
        else
            d[p.state] = p.weight
        end
        if d[p.state] > best_weight
            best_weight = d[p.state]
            most_likely = p.state
        end
    end
    return most_likely
end

"""
DEPRECATED: Please use https://github.com/JuliaPOMDP/ParticleFilters.jl instead

Updater for ParticleBelief that implements
the sampling importance resampling (SIR) algorithm. Fields:

    -pomdp: problem model
    -n: number of particles for this updater
    -rng: random number generator for generating state and observations from pomdp model
    -keep_dict: boolean for keeping a dictionary of particles (for pdf only)
"""
mutable struct SIRParticleUpdater <: Updater
    pomdp::POMDP # POMDP model
    n::Int64 # number of particles
    rng::AbstractRNG
    keep_dict::Bool
    function SIRParticleUpdater(pomdp::POMDP, n::Int64, rng::AbstractRNG, keep_dict::Bool)
        warn("POMDPToolbox.SIRParticleUpdater is deprecated. Please use https://github.com/JuliaPOMDP/ParticleFilters.jl instead.")
        new(pomdp, n, rng, keep_dict)
    end
end
function SIRParticleUpdater(pomdp::POMDP, n::Int64; rng::AbstractRNG=Base.GLOBAL_RNG, keep_dict::Bool=true)
    SIRParticleUpdater(pomdp, n, rng, keep_dict)
end

function create_belief(up::SIRParticleUpdater)
    st = state_type(typeof(up.pomdp))
    particles = Array{Particle{st}}(up.n)
    w = 1.0 / up.n
    probs = Dict{st, Float64}()
    probs_arr = zeros(up.n)
    return ParticleBelief(particles, probs, probs_arr, up.keep_dict)
end

function initialize_belief(up::SIRParticleUpdater, d::Any, b::ParticleBelief = create_belief(up))
    w = 1.0 / up.n # start with a uniform weight
    for i = 1:up.n
        # sample a state from initial dist
        s = rand(up.rng, d)
        # add particle
        p = Particle(s, w)
        b.particles[i] = p
        # update dict
        if haskey(b.probs_dict, s)
            b.probs_dict[s] += w
        else
            b.probs_dict[deepcopy(s)] = w
        end
    end
    return b
end

function update{A,O}(bu::SIRParticleUpdater, bold::ParticleBelief, a::A, o::O, bnew::ParticleBelief=create_belief(bu))

    rng = bu.rng
    particles = bold.particles
    pomdp = bu.pomdp

    for i = 1:bu.n
        # step particles forward
        s = particles[i].state
        sp = generate_s(pomdp, s, a, rng)
        # compute obs likelihoods
        od = observation(pomdp, s, a, sp)
        bnew.probs_arr[i] = pdf(od, o)
        particles[i].state = sp
    end

    particles = normalize!(particles)
    bnew.probs_arr /= sum(bnew.probs_arr)

    cat = Weights(bnew.probs_arr)

    w = 1.0 / bu.n
    bu.keep_dict ? (empty!(bnew.probs_dict)) : (nothing)
    # resample
    for i = 1:bu.n
        k = sample(rng, cat)
        sp = particles[k].state
        bnew.particles[i] = Particle(deepcopy(sp), w) # XXX slow?
        if bu.keep_dict
            haskey(bnew.probs_dict, sp) ? (bnew.probs_dict[sp] += w) : (bnew.probs_dict[deepcopy(sp)] = w)
        end
    end
    fill!(bnew.probs_arr, w)
    return bnew
end

function normalize!{s}(particles::Vector{Particle{s}})
    prob_sum = 0.0
    for p in particles
        prob_sum += p.weight
    end
    for p in particles
        p.weight /= prob_sum
    end
    return particles
end


text:
## This is our versions of Table with the loop variable having lexical scope.
## It is disabled because function names would clobber the same ones used in
## the dynamic scoping version in table_dyn.jl
##
## Table sets a symbol sym to a sequence of values, evaluating expr
## each time and building a list of the results.
##
## We replace sym everywhere it appears in expr by a gensym so that we don't
## clobber other bindings of sym.
##
## If there are few replacements to make, another approach appears to be faster (by about a factor of 2):
## Instead of substituting a unique symbol, we record the positions in expr at which sym appears.
## Then, at each iteration, we subsitute the value directly into (a copy of) expr and evaluate.
##
## But, if expr is large and sym appears many times, using the gensym approach is much faster.
##
## We should add some logic to choose between the strategies.
##
## To start, we use types representing strategies <: TableEval. We first implement strategies
## that use setsymval. Next we will implement those that substitute the value directly.
## Then we need a criterion for dispatch. Probably just the number of subsitutions is enough
## for zeroeth order. But, the setsymval strategy is best if we only have one strategy. It
## seems to be at worst about half as fast as direct substitution, wheras the latter can get
## arbitrarily bad if there are many substitutions.
##

abstract type TableEval end

struct TableSetEval{T,V} <: TableEval
    expr::T
    sym::V
end

struct TableNoSetEval{T} <: TableEval
    expr::T
end

@sjdoc NTable """
    NTable(expr,[imax])

return a `List` of `imax` copies of `expr`.

    NTable(expr,[i,imax])

returns a `List` of `expr` evaluated `imax` times with `i` set successively to `1` through `imax`.

    NTable(expr,iter)

`iter` can be any standard iterator.

    NTable(expr,iter1,iter2,...)

is equivalent to `NTable(NTable(expr,iter2),iter1)...`

This example calls an anonymous Julia function
```
f = :( g(x) = x^2 )
NTable( f(i), [i,10])
```
"""

@mkapprule NTable

@doap function NTable(inexpr, iters...)
    if length(iters) > 1
        riters = reverse(iters)
        mxout = mxpr(:NTable,inexpr,riters[1])
        for i in 2:length(riters)
            mxout = mxpr(:NTable,mxout,riters[i])
        end
        return mxout
    end
    expr = is_Mxpr(inexpr,:Jxpr) ? eval(inexpr) : deepcopy(inexpr)
    table(expr, iters)
end

function table(expr, iters)
    length(iters) == 1 && return mxpr(:List, intable(expr, make_sjiter(iters[1])))
end

function intable(expr, iter::SJIter1)
    tableloop1(expr, iter.imax)
end

function intable(expr, iter::SJIter2)
    (sym, nexpr) = localize_variable(iter.i,expr)
    tableloop2(sym, nexpr, iter.imax)
end

function intable(expr, iter::SJIterList)
    (sym, nexpr) = localize_variable(iter.i,expr)
    tableloop_list(sym, nexpr, iter.list)
end

function evaltable(model::TableSetEval, val)
    setandeval(model.expr, model.sym, val)
end

function evaltable(model::TableNoSetEval)
    nosetandeval(model.expr)
end

function setandeval(expr,sym,val)
    setsymval(sym,val)
    unsetfixed(expr)
    r = doeval(expr)
    setfixed(r) # trying to gain a little efficiency
    setcanon(r) # usually fails in some more or less edge case
    r
end

function nosetandeval(expr)
    unsetfixed(expr)
    r = doeval(expr)
    setfixed(r) # trying to gain a little efficiency
    setcanon(r) # usually fails in some more or less edge case
    r
end

function intable(expr0,iter::SJIter3{T,V}) where {T<:Real,V<:Real}
    (sym, expr) = localize_variable(iter.i,expr0)
    imax = doeval(iter.imax) # maybe this should be done earlier. When iter is created ?
    imin = doeval(iter.imin)
    args = newargs(imax-imin+1)
    if iter.i == expr
        @inbounds for i in imin:imax
            args[i-imin+1] = i
        end
    else
        model = TableSetEval(expr,sym)
        _table3(model,imin,imax,args)
    end
    return args
end

function _table3(model,imin,imax,args)
    for i in imin:imax
        args[i-imin+1] = evaltable(model,i)
    end
end

function intable(expr0, iter::SJIter4{T,V,W}) where {T<:Real, V<:Real, W<:Real}
    (sym, expr) = localize_variable(iter.i,expr0)
    imax = iter.imax  # why do we not evaluate here
    imin = iter.imin
    di = iter.di
    args = newargs(iter.num_iters)
    j::Int = 0
    if iter.i == expr
        @inbounds for i in imin:di:imax
            j += 1
            args[j] = i
        end
    else
        model = TableSetEval(expr,sym)
        @inbounds for i in imin:di:imax
            j += 1
            args[j] = evaltable(model,i)
        end
    end
    return args
end

# Symbolic value and increment
function intable(expr0, iter::SJIter4{T,V,W}) where {T, V, W}
    (sym, expr) = localize_variable(iter.i,expr0)
    imax = iter.imax
    imin = iter.imin
    di = iter.di
    num_iters = iter.num_iters
    args = newargs(iter.num_iters)
    model = TableSetEval(expr,sym)
    val = deepcopy(imin)
    @inbounds for i in 1:num_iters
        args[i] = evaltable(model,val)
        val = doeval(mxpr(:Plus, val, di))
    end
    args
end

function tableloop1(expr, imax)
    args = newargs(imax)
    model = TableNoSetEval(expr)
    @inbounds  for i=1:imax
        args[i] = evaltable(model)
    end
    args
end

function tableloop2(sym, expr, imax)
    args = newargs(imax)
    model = TableSetEval(expr,sym)
    @inbounds  for i=1:imax
        args[i] = evaltable(model,i)
    end
    args
end

function table3(expr,iter::SJIter3{T,V}) where {T<:Real,V<:Real}
    imax = doeval(iter.imax) # maybe this should be done earlier. When iter is created ?
    imin = doeval(iter.imin)
    args = newargs(imax-imin+1)
    if iter.i == expr
        @inbounds for i in imin:imax
            args[i-imin+1] = i
        end
    else
        model = TableSetEval(expr,sym)
        @inbounds for i in imin:imax
            args[i-imin+1] =  evaltable(model,i)
        end
    end
    return args
end

function tableloop_list(sym, expr, list)
    len = length(list)
    args = newargs(len)
    model = TableSetEval(expr,sym)
    for i in 1:len
        args[i] = evaltable(model,list[i])
    end
    args
end


text:
module MortalityTables
using HTTP
using Transducers
using OffsetArrays
using Parsers
using QuadGK
using Requires
using UnPack
using XMLDict
using Pkg.Artifacts

include("MetaData.jl")
include("death_distribution.jl")
include("MortalityTable.jl")
include("XTbML.jl")
include("get_SOA_table.jl")
include("parameterized_models.jl")

export MortalityTable,
    survival,
    decrement,
    omega,
    TableMetaData,
    SelectMortality,
    UltimateMortality,
    Balducci,
    Uniform,
    Constant,
    DeathDistribution,
    get_SOA_table,
    get_SOA_table!,
    Makeham, Gompertz, MakehamGompertz,
    hazard,cumhazard,
    mortality_vector

# lazy load part of the package
function __init__()
    @require CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" include("CSV.jl")
end

end # module


text:
import Gaugefields:Traceless_antihermitian_add!

struct Wilsonclover_data
end

abstract type Wilsontype_FermiAction{Dim,Dirac,fermion,gauge} <: FermiAction{Dim,Dirac,fermion,gauge}
end

struct WilsonFermiAction{Dim,Dirac,fermion,gauge,hascloverterm} <: Wilsontype_FermiAction{Dim,Dirac,fermion,gauge} #FermiAction{Dim,Dirac,fermion,gauge}
    hascovnet::Bool
    covneuralnet::Union{Nothing,CovNeuralnet{Dim}}
    diracoperator::Dirac
    hascloverterm::Bool
    clover_data::Union{Nothing,Wilsonclover_data}
    _temporary_fermionfields::Vector{fermion}
    _temporary_gaugefields::Vector{gauge}

    function WilsonFermiAction(D::Dirac_operator{Dim},hascovnet,covneuralnet,parameters_action) where Dim
        hascloverterm = check_parameters(parameters_action,"hascloverterm",false)
        if hascloverterm
            error("not implemented yet")
        else
            clover_data = nothing
        end

        num = 6

        x = D._temporary_fermi[1]
        xtype = typeof(x)
        _temporary_fermionfields = Array{xtype,1}(undef,num)
        for i=1:num
            _temporary_fermionfields[i] = similar(x)
        end

        Utemp = D.U[1]
        Utype = typeof(Utemp)
        numU = 2
        _temporary_gaugefields = Array{Utype,1}(undef,numU)
        for i=1:numU
            _temporary_gaugefields[i] = similar(Utemp)
        end

        return new{Dim,typeof(D),xtype,Utype,hascloverterm}(hascovnet,covneuralnet,D,
                            hascloverterm,
                            clover_data,
                        _temporary_fermionfields,
                        _temporary_gaugefields)

    end
end

function evaluate_FermiAction(fermi_action::Wilsontype_FermiAction{Dim,Dirac,fermion,gauge},U,ϕ::AbstractFermionfields) where {Dim,Dirac,fermion,gauge,hascloverterm}
    W = fermi_action.diracoperator(U)
    η = fermi_action._temporary_fermionfields[1]
    solve_DinvX!(η,W',ϕ)
    Sf = dot(η,η)
    return real(Sf)
end

function calc_UdSfdU!(UdSfdU::Vector{<: AbstractGaugefields},fermi_action::Wilsontype_FermiAction{Dim,Dirac,fermion,gauge} ,U::Vector{<: AbstractGaugefields},ϕ::AbstractFermionfields) where  {Dim,Dirac,fermion,gauge,hascloverterm}
    #println("------dd")
    W = fermi_action.diracoperator(U)
    WdagW = DdagD_Wilson_operator(W)
    X = fermi_action._temporary_fermionfields[end]
    Y = fermi_action._temporary_fermionfields[2]
    #X = (D^dag D)^(-1) ϕ
    #
    #println("Xd ",X[1,1,1,1,1,1])
    solve_DinvX!(X,WdagW,ϕ)
    #println("X ",X[1,1,1,1,1,1])
    clear_U!(UdSfdU)

    calc_UdSfdU_fromX!(UdSfdU,Y,fermi_action,U,X)
    #println("----aa--")
    set_wing_U!(UdSfdU)
end

function calc_UdSfdU_fromX!(UdSfdU::Vector{<: AbstractGaugefields},Y,fermi_action::WilsonFermiAction{Dim,Dirac,fermion,gauge,hascloverterm} ,U,X;coeff = 1) where {Dim,Dirac,fermion,gauge,hascloverterm}
    W = fermi_action.diracoperator(U)
    mul!(Y,W,X)
    #set_wing_fermion!(Y)

    temp0_f = fermi_action._temporary_fermionfields[3]
    temp1_f = fermi_action._temporary_fermionfields[4]
    temp0_g = fermi_action._temporary_gaugefields[1]

    for μ=1:Dim
        #!  Construct U(x,mu)*P1

        # U_{k,μ} X_{k+μ}
        Xplus = shift_fermion(X,μ)

        #@time mul!(temp0_f,U[μ],X)

        mul!(temp0_f,U[μ],Xplus)

        # (r-γ_μ) U_{k,μ} X_{k+μ}
        mul!(temp1_f,view(W.rminusγ,:,:,μ),temp0_f)

        # κ (r-γ_μ) U_{k,μ} X_{k+μ}
        mul!(temp0_f,W.hopp[μ],temp1_f)

        # κ ((r-γ_μ) U_{k,μ} X_{k+μ}) ⊗ Y_k
        mul!(temp0_g,temp0_f,Y')

        add_U!(UdSfdU[μ],-coeff,temp0_g)

        #!  Construct P2*U_adj(x,mu)
        # Y_{k+μ}^dag U_{k,μ}^dag
        Yplus = shift_fermion(Y,μ)
        mul!(temp0_f,Yplus',U[μ]')

        # Y_{k+μ}^dag U_{k,μ}^dag*(r+γ_μ)
        mul!(temp1_f,temp0_f,view(W.rplusγ,:,:,μ))

        # κ Y_{k+μ}^dag U_{k,μ}^dag*(r+γ_μ)
        mul!(temp0_f,W.hopm[μ],temp1_f)

        # X_k ⊗ κ Y_{k+μ}^dag U_{k,μ}^dag*(r+γ_μ)
        mul!(temp0_g,X,temp0_f)

        add_U!(UdSfdU[μ],coeff,temp0_g)

    end

end

function calc_p_UdSfdU!(p,fermi_action::WilsonFermiAction{Dim,Dirac,fermion,gauge,hascloverterm} ,U::Vector{<: AbstractGaugefields},ϕ::AbstractFermionfields,coeff = 1) where  {Dim,Dirac,fermion,gauge,hascloverterm}
    #println("------dd")
    W = fermi_action.diracoperator(U)
    WdagW = DdagD_Wilson_operator(W)
    X = fermi_action._temporary_fermionfields[end]
    Y = fermi_action._temporary_fermionfields[2]
    #X = (D^dag D)^(-1) ϕ
    #
    #println("Xd ",X[1,1,1,1,1,1])
    solve_DinvX!(X,WdagW,ϕ)
    #println("X ",X[1,1,1,1,1,1])
    #clear_U!(UdSfdU)

    calc_p_UdSfdU_fromX!(p,Y,fermi_action,U,X,coeff= coeff)
    #println("----aa--")
    #set_wing_U!(UdSfdU)
end

function calc_p_UdSfdU_fromX!(p,Y,fermi_action::WilsonFermiAction{Dim,Dirac,fermion,gauge,hascloverterm} ,U,X;coeff = 1) where {Dim,Dirac,fermion,gauge,hascloverterm}
    W = fermi_action.diracoperator(U)
    mul!(Y,W,X)
    #set_wing_fermion!(Y)

    temp0_f = fermi_action._temporary_fermionfields[3]
    temp1_f = fermi_action._temporary_fermionfields[4]
    temp0_g = fermi_action._temporary_gaugefields[1]

    for μ=1:Dim
        #!  Construct U(x,mu)*P1

        # U_{k,μ} X_{k+μ}
        Xplus = shift_fermion(X,μ)

        #@time mul!(temp0_f,U[μ],X)

        mul!(temp0_f,U[μ],Xplus)

        # (r-γ_μ) U_{k,μ} X_{k+μ}
        mul!(temp1_f,view(W.rminusγ,:,:,μ),temp0_f)

        # κ (r-γ_μ) U_{k,μ} X_{k+μ}
        mul!(temp0_f,W.hopp[μ],temp1_f)

        # κ ((r-γ_μ) U_{k,μ} X_{k+μ}) ⊗ Y_k
        mul!(temp0_g,temp0_f,Y')

        Traceless_antihermitian_add!(p[μ],-coeff,temp0_g)

        #add_U!(UdSfdU[μ],-coeff,temp0_g)

        #!  Construct P2*U_adj(x,mu)
        # Y_{k+μ}^dag U_{k,μ}^dag
        Yplus = shift_fermion(Y,μ)
        mul!(temp0_f,Yplus',U[μ]')

        # Y_{k+μ}^dag U_{k,μ}^dag*(r+γ_μ)
        mul!(temp1_f,temp0_f,view(W.rplusγ,:,:,μ))

        # κ Y_{k+μ}^dag U_{k,μ}^dag*(r+γ_μ)
        mul!(temp0_f,W.hopm[μ],temp1_f)

        # X_k ⊗ κ Y_{k+μ}^dag U_{k,μ}^dag*(r+γ_μ)
        mul!(temp0_g,X,temp0_f)

        Traceless_antihermitian_add!(p[μ],coeff,temp0_g)

        #add_U!(UdSfdU[μ],coeff,temp0_g)

    end

end

function gauss_sampling_in_action!(η::AbstractFermionfields,U,fermi_action::Wilsontype_FermiAction{Dim,Dirac,fermion,gauge} ) where {Dim,Dirac,fermion,gauge,hascloverterm}
    #gauss_distribution_fermion!(η)
    gauss_distribution_fermion!(η,rand)
end

using InteractiveUtils

function sample_pseudofermions!(ϕ::AbstractFermionfields,U,fermi_action::Wilsontype_FermiAction{Dim,Dirac,fermion,gauge} ,ξ::AbstractFermionfields) where {Dim,Dirac,fermion,gauge,hascloverterm}
    W = fermi_action.diracoperator(U)
    mul!(ϕ,W',ξ)
    set_wing_fermion!(ϕ)
end


text:
module HSL

using LinearAlgebra
using SparseArrays

if isfile(joinpath(dirname(@__FILE__), "..", "deps", "deps.jl"))
  include("../deps/deps.jl")
else
  error("HSL library not properly installed. Please run Pkg.build(\"HSL\")")
end

# definitions applicable to all packages
const data_map = Dict{Type, Type}(Float32 => Cfloat,
                                  Float64 => Cdouble,
                                  ComplexF32 => Cfloat,
                                  ComplexF64 => Cdouble)

# package-specific definitions
if @isdefined libhsl_ma97
  include("hsl_ma97.jl")
end

end


text:
mutable struct LevenbergMarquardt <: OptimizationMethod
  epsfcn::Float64
  xtol::Float64
  gtol::Float64
  useCostFunctionsJacobin::Bool
end

LevenbergMarquardt() = LevenbergMarquardt(EPS_VAL, 1.0e-8, 1.0e-8, false)

## Levenberg Marquardt Methods ##
function minimize!(lm::LevenbergMarquardt, p::Problem, endCriteria::EndCriteria)
  reset!(p)
  x = p.currentValue
  initCostValues = values!(p, x)

  m = length(initCostValues)
  n = length(x)
  if lm.useCostFunctionsJacobin
    initJacobin = zeros(m, n)
    get_jacobin!(p.costFunction, initJacobin, x)
  end
  fvec = zeros(m)
  diag_ = zeros(n)
  mode = 1
  factor_ = 1.0
  nprint = 0
  info_ = 0
  nfev = 0
  fjac = zeros(m, n)
  ldfjac = m
  ipvt = ones(Int, n)
  qtf = zeros(n)
  wa1 = zeros(n)
  wa2 = zeros(n)
  wa3 = zeros(n)
  wa4 = zeros(m)

  function fcn!(::Int, n::Int, _x::Vector{Float64}, _fvec::Vector{Float64})
    xt = _x[1:n]

    if test(p.constraint, xt)
      tmp = values!(p, xt)
      _fvec[1:length(tmp)] = tmp
    else
      _fvec = copy(initCostValues)
    end

    return _fvec
  end

  # function jacFcn!{I <: Integer}(m::I, n::I, x::Float64, fjac::Float64, ::I)
  #   xt = fill(x + n, n)
  # end

  # function fcn2(x::Vector{Float64})
  #   if test(p.constraint, x)
  #     fvec = values!(p, x)
  #   else
  #     fvec = copy(initCostValues)
  #   end
  #
  #   return fvec
  # end

  # jacFcn! = Calculus.jacobian(fcn2)

  # TODO check requirements, see levenbergmarquardt.cpp
  lmdif!(m, n, x, fvec, endCriteria.functionEpsilon, lm.xtol, lm.gtol, endCriteria.maxIterations, lm.epsfcn, diag_, mode, factor_, nprint, info_, nfev, fjac, ldfjac,
      ipvt, qtf, wa1, wa2, wa3, wa4, fcn!)

  # lmdif2!(n, m, x, mode, factor_, info_, lm.epsfcn, endCriteria.functionEpsilon, lm.xtol, lm.gtol, endCriteria.maxIterations, fcn!)
  # println(x)
  # res = lmfit(fcn2, xx, endCriteria.maxIterations)

  # current value is already set
  p.functionValue = QuantLib.value(p.costFunction, x)

  return p
end


text:
using Random
# using POMDPSimulators
using MCTS

function simple_mdps()
    pp = simple_pp()
    PMDPs.PMDPg(pp), PMDPs.PMDPe(pp)
end

@testset "policy_tools.jl" begin
    mg, me = simple_mdps()
    s₀ = rand(POMDPs.initialstate(mg))

    policy = PMDPs.get_VI_policy(me)
    @test isa(policy, Policy)
    @test isa(action(policy,s₀), PMDPs.Action)

    planner = PMDPs.get_MCTS_planner(mg)
    @test isa(planner, Policy)
    @test isa(action(planner,s₀), PMDPs.Action)

    params_mcts = Dict(pairs( (solver= DPWSolver, n_iterations=50, depth=30, exploration_constant=40.0, enable_state_pw = true, keep_tree=true, show_progress=false)))
    planner2 = PMDPs.get_MCTS_planner(mg; params_mcts)
    @test isa(planner2, Policy)
    @test isa(action(policy,s₀), PMDPs.Action)

end


text:
# Solves the Poisson equation using the Mixed finite element method
using Revise
using AdFem
using DelimitedFiles
using SparseArrays
using PyPlot

n = 50
mmesh = Mesh(n, n, 1/n, degree = BDM1)

testCase = [
    ((x, y)->begin
        x * (1-x) * y * (1-y)
    end,
    (x, y)->begin
        -2x*(1-x) -2y*(1-y)
    end),  # test case 1
    ((x, y)->begin
        x^2 * (1-x) * y * (1-y)^2
    end,
    (x, y)->begin
        2*x^2*y*(1 - x) + 2*x^2*(1 - x)*(2*y - 2) - 4*x*y*(1 - y)^2 + 2*y*(1 - x)*(1 - y)^2
    end) ,# test case 2
    (
        (x, y)->x*y,
        (x, y)->0.0
    ), # test case 3
    (
        (x, y)->x^2 * y^2 + 1/(1+x^2),
        (x, y)->2*x^2 + 8*x^2/(x^2 + 1)^3 + 2*y^2 - 2/(x^2 + 1)^2
    ), # test case 4
]

for k = 1:4
    @info "TestCase $k..."
    ufunc, ffunc = testCase[k]

    A = compute_fem_bdm_mass_matrix1(mmesh)
    B = compute_fem_bdm_div_matrix1(mmesh)
    C = [A -B'
        -B spzeros(mmesh.nelem, mmesh.nelem)]

    gD = bcedge(mmesh)
    t1 = eval_f_on_boundary_edge(ufunc, gD, mmesh)
    g = compute_fem_traction_term1(t1, gD, mmesh)
    t2 = eval_f_on_gauss_pts(ffunc, mmesh)
    f = compute_fvm_source_term(t2, mmesh)
    rhs = [-g; f]

    sol = C\rhs
    u = sol[mmesh.ndof+1:end]
    close("all")
    figure(figsize=(15, 5))
    subplot(131)
    title("Reference")
    xy = fvm_nodes(mmesh)
    x, y = xy[:,1], xy[:,2]
    uf = ufunc.(x, y)
    visualize_scalar_on_fvm_points(uf, mmesh)
    subplot(132)
    title("Numerical")
    visualize_scalar_on_fvm_points(u, mmesh)
    subplot(133)
    title("Absolute Error")
    visualize_scalar_on_fvm_points( abs.(u - uf) , mmesh)
    savefig("bdm$k.png")
end


text:
# Author: Mathias Louboutin, mlouboutin3@gatech.edu
# Date: December 2020

using LinearAlgebra

N = 10
A = Diagonal(1:N)

function proj(x)
    xp = deepcopy(x)
    xp[xp .< 0] .= 0
    return xp
end

x0 = 10 .+ 10 .* rand(N)
b = A*x0

function obj(x)
    fun = .5*norm(A*x - b)^2
    grad = A'*(A*x - b)
    return fun, grad
end

for ls in [BackTracking, HagerZhang, StrongWolfe]
    opt = spg_options(maxIter=1000, progTol=1f-6, optTol=1f-6, verbose=2)
    sol = spg(obj, randn(N), proj, opt, ls())

    @show norm(sol.x - x0)/(norm(x0) + norm(sol.x))
    @show sol.ϕ/sol.ϕ_trace[1]

    @test sol.ϕ/sol.ϕ_trace[1]  < 1e-9
    @test norm(sol.x - x0)/(norm(x0) + norm(sol.x)) < 1f-4
end


text:
# Range: 1 : 10,000
# Time:  13.548s

# Written by Alex Vear in 2019
# Public domain.  No rights reserved.

startval = 1
endval = 10000

for num in startval:endval
  sqr = num^2
  for i in startval:endval
    tri = (i*(i+1))/2
    if sqr == tri
      println("Square Number: $sqr Triangle Number: $tri Squared: $num Triangled: $i")
    end
  end
end

max_stars_repo_name: axvr/codedump

max_issues_repo_name: axvr/codedump

max_forks_repo_name: axvr/codedump


text:
module BouncingCapsulesSimulation

using  ModiaLang
import Modia3D
using  Modia3D.ModiaInterface
using  Unitful

vmatRed   = Modia3D.VisualMaterial(color="Red")
vmatGreen = Modia3D.VisualMaterial(color="Green")
vmatBlue  = Modia3D.VisualMaterial(color="Blue")
vmatGrey  = Modia3D.VisualMaterial(color="Grey", transparency=0.5)

BouncingCapsules = Model(
    boxHeigth = 0.1,
    gravField = UniformGravityField(g=9.81, n=[0, 0, -1]),
    world = Object3D(feature=Scene(gravityField=:gravField,
                                   visualizeFrames=false,
                                   defaultFrameLength=0.2,
                                   visualizeBoundingBox = true,
                                   enableContactDetection=true,
                                   visualizeContactPoints=false)),
    worldFrame = Object3D(parent=:world, feature=Visual(shape=CoordinateSystem(length=0.5))),
    ground = Object3D(parent=:world,
                      translation=:[0.0, 0.0, -boxHeigth/2],
                      feature=Solid(shape=Box(lengthX=5.0, lengthY=3.0, lengthZ=:boxHeigth),
                                    visualMaterial=vmatGrey,
                                    solidMaterial="Steel",
                                    collision=true)),
    frameX = Object3D(parent=:world,
                      translation=:[-1.0, 0.0, 1.0],
                      rotation=:[-90*u"°", 0.0, -90*u"°"],
                      feature=Visual(shape=CoordinateSystem(length=0.5))),
    capsuleX = Object3D(feature=Solid(shape=Capsule(axis=1, diameter=0.4, length=1.0),
                                      visualMaterial=vmatRed,
                                      solidMaterial="Steel",
                                      collision=true)),
    jointX = FreeMotion(obj1=:frameX, obj2=:capsuleX,
                        r=Var(init=[0.0, 0.0, 0.0]),
                        rot=Var(init=[0.0, -60*u"°", 0.0]),
                        v=Var(init=[0.0, 1.0, 0.0])),
    frameY = Object3D(parent=:world,
                      translation=:[0.0, 0.0, 1.0],
                      rotation=:[90*u"°", 90*u"°", 0.0],
                      feature=Visual(shape=CoordinateSystem(length=0.5))),
    capsuleY = Object3D(feature=Solid(shape=Capsule(axis=2, diameter=0.4, length=1.0),
                                      visualMaterial=vmatGreen,
                                      solidMaterial="Steel",
                                      collision=true)),
    jointY = FreeMotion(obj1=:frameY, obj2=:capsuleY,
                        r=Var(init=[0.0, 0.0, 0.0]),
                        rot=Var(init=[0.0, 0.0, -60*u"°"]),
                        v=Var(init=[0.0, 0.0, 1.0])),
    frameZ = Object3D(parent=:world,
                      translation=:[1.0, 0.0, 1.0],
                      feature=Visual(shape=CoordinateSystem(length=0.5))),
    capsuleZ = Object3D(feature=Solid(shape=Capsule(axis=3, diameter=0.4, length=1.0),
                                      visualMaterial=vmatBlue,
                                      solidMaterial="Steel",
                                      collision=true)),
    jointZ = FreeMotion(obj1=:frameZ, obj2=:capsuleZ,
                        r=Var(init=[0.0, 0.0, 0.0]),
                        rot=Var(init=[-60*u"°", 0.0, 0.0]),
                        v=Var(init=[1.0, 0.0, 0.0]))
)

bouncingCapsules = @instantiateModel(buildModia3D(BouncingCapsules), unitless=true, log=false, logStateSelection=false, logCode=false)

#@show bouncingCapsules.parameterExpressions
#@show bouncingCapsules.parameters

stopTime = 1.7
tolerance = 1e-8
requiredFinalStates = [-0.7964014792034407, 1.078911579441921, 0.23035036015771998, -0.10969534535624982, 0.4797449365169936, -0.08174820851332286, 4.48271590323064, 1.340697608488989, -1.5392203655926255, 2.4243895450782937, 0.15531689730526202, -0.03915705292205236, 0.22862079147109046, -0.7970862046473598, 1.0806110815754162, -0.06460242891125959, -0.11431031161481991, 0.4851429973131076, 3.385569081576411, 0.1916283718909546, -1.6236830622142888, -0.06474401604033102, 2.4344556240272657, 0.2553659880626085, 1.0818397597553406, 0.22945140558263866, -0.7966137267737096, 0.4853616788402741, -0.07085422346062514, -0.09972138522555073, -1.5640086933709134, 0.19867779732620508, 3.386466331485274, 0.21808527847860343, -0.05743721435841237, 2.4419014448796816]
simulate!(bouncingCapsules, stopTime=stopTime, tolerance=tolerance, log=true, logStates=true, logEvents=true, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plot(bouncingCapsules, [("jointX.r", "jointY.r", "jointZ.r") ("jointX.rot", "jointY.rot", "jointZ.rot")
                        ("jointX.v", "jointY.v", "jointZ.v") ("jointX.w", "jointY.w", "jointZ.w")], figure=1)

end