text:
# updates.jl
#
# Defines the functions that perform the diagonal update, and also
# that build the linked list and operator cluster update

lsize = 0
nullt = (0,0,0) #a null tuple

############################ FUNCTIONS ######################################

#Diagonal update
function DiagonalUpdate()

#define the Metropolis probability as a constant
#https://pitp.phas.ubc.ca/confs/sherbrooke2012/archives/Melko_SSEQMC.pdf
#equation 1.43
P_h = h_x*nSpin/(h_x*nSpin +2.0*J_*nBond) #J=1.0 tested only

    spin_prop = copy(spin_left)  #the propagated spin state
    for i = 1:2*M  #size of the operator list

       if operator_list[i,1] == -2
           spin_prop[operator_list[i,2]] = xor(spin_prop[operator_list[i,2]],1) #spinflip

       else
           flag = false
           while flag == false
               rr = rand()
               if P_h > rr #probability to choose a single-site operator
                   operator_list[i,1] = -1
                   site = rand(1:nSpin)
                   operator_list[i,2] = site
                   flag = true
                   #println(i," site ",site)
               else
                   bond = rand(1:nBond)
                   if spin_prop[bond_spin[bond,1]] == spin_prop[bond_spin[bond,2]] #spins must be the same
                       operator_list[i,1] = bond_spin[bond,1]
                       operator_list[i,2] = bond_spin[bond,2]
                       flag = true
                       #println(i," bond ",bond)
                   end#if
                   #println(P_h," ",rr," bond")
               end
           end #while

        end#if
    end #for

    #DEBUG
    if spin_prop != spin_right  #check the spin propagation for error
        println("Basis state propagation error: DiagonalUpdate")
    end

end #DiagonalUpdate

#############################################################################

#LinkedList
function LinkedList()

    #initialize linked list data structures
    global LinkList = zeros(Int,0)  #needed for cluster update
    global LegType = zeros(Int,0)

    #A diagonal bond operator has non trivial associates for cluster building
    #Associates = zeros((Int,Int,Int),0)
    global Associates = []

    First = zeros(Int,0)  #scope is this function only

    #The first N elements of the linked list are the spins of the LHS basis state
    for i = 1:nSpin
        push!(First,i)
        push!(LinkList,-99) #we don't yet know what these link to
        push!(LegType,spin_left[i])
        push!(Associates,nullt) #no nontrivial associates for train wf spins
    end #i

    spin_prop = copy(spin_left)  #the propagated spin state

    #Now, add the 2M operators to the linked list.  Each has either 2 or 4 legs
    for i = 1:2*M  #size of the operator list

       if operator_list[i,1] == -2  #off-diagonal site operator
           site = operator_list[i,2]
           #lower or left leg
           push!(LinkList,First[site])
           push!(LegType,spin_prop[site]) #the spin of the vertex leg
           spin_prop[site] = xor(spin_prop[site],1) #spinflip
           current_link = size(LinkList)
           LinkList[First[site]] = current_link[1] #completes backwards link
           First[site] = current_link[1] + 1
           push!(Associates,nullt)
           #upper or right leg
           push!(LinkList,-99) #we don't yet know what this links to
           push!(LegType,spin_prop[site]) #the spin of the vertex leg
           push!(Associates,nullt)

       elseif operator_list[i,1] == -1  #diagonal site operator
           site = operator_list[i,2]
           #lower or left leg
           push!(LinkList,First[site])
           push!(LegType,spin_prop[site]) #the spin of the vertex leg
           current_link = size(LinkList)
           LinkList[First[site]] = current_link[1] #completes backwards link
           First[site] = current_link[1] + 1
           push!(Associates,nullt)
           #upper or right leg
           push!(LinkList,-99) #we don't yet know what this links to
           push!(LegType,spin_prop[site]) #the spin of the vertex leg
           push!(Associates,nullt)

       else  #diagonal bond operator
           #lower left
           site1 = operator_list[i,1]
           push!(LinkList,First[site1])
           push!(LegType,spin_prop[site1]) #the spin of the vertex leg
           current_link = size(LinkList)
           LinkList[First[site1]] = current_link[1] #completes backwards link
           First[site1] = current_link[1] + 2
           vertex1 = current_link[1]
           push!(Associates,(vertex1+1,vertex1+2,vertex1+3))
           #lower right
           site2 = operator_list[i,2]
           push!(LinkList,First[site2])
           push!(LegType,spin_prop[site2]) #the spin of the vertex leg
           current_link = size(LinkList)
           LinkList[First[site2]] = current_link[1] #completes backwards link
           First[site2] = current_link[1] + 2
           push!(Associates,(vertex1,vertex1+2,vertex1+3))
           #upper left
           push!(LinkList,-99) #we don't yet know what this links to
           push!(LegType,spin_prop[site1]) #the spin of the vertex leg
           push!(Associates,(vertex1,vertex1+1,vertex1+3))
           #upper right
           push!(LinkList,-99) #we don't yet know what this links to
           push!(LegType,spin_prop[site2]) #the spin of the vertex leg
           push!(Associates,(vertex1,vertex1+1,vertex1+2))

       end #if

    end #i

    #The last N elements of the linked list are the final spin state
    for i = 1:nSpin
        push!(LinkList,First[i])
        push!(LegType,spin_prop[i])
        current_link = size(LinkList)
        LinkList[First[i]] = current_link[1]
        push!(Associates,nullt)
    end #i

    global lsize = size(LinkList)

    #DEBUG
    if spin_prop != spin_right
        println("Basis state propagation error: LINKED LIST")
    end

end #LinkedList

#############################################################################

#ClusterUpdate
function ClusterUpdate()

    #lsize is the size of the linked list
    in_cluster=zeros(Int,lsize[1])

    cstack = zeros(Int,0)  #This is the stack of vertices in a cluster

    ccount = 0 #cluster number counter
    for i = 1:lsize[1]

        #Add a new leg onto the cluster
        if (in_cluster[i] == 0 && Associates[i] == nullt)

            ccount+=1
            push!(cstack,i)
            in_cluster[cstack[end]] = ccount

            flip = rand(Bool) #flip a coin for the SW cluster flip
            if flip == true
                LegType[cstack[end]] =  xor(LegType[cstack[end]],1) #spinflip
            end

            while isempty(cstack) == false

                leg = LinkList[cstack[end]]
                pop!(cstack)
                #println("leg ",leg," ",cstack)

                if in_cluster[leg] == 0

                    in_cluster[leg] = ccount; #add the new leg and flip it
                    if flip == true
                        #println("gonna flip ",leg)
                        LegType[leg] =  xor(LegType[leg],1)
                    end
                    #now check all associates and add to cluster
                    assoc = Associates[leg] #a 3-element array
                    if assoc != nullt
                        push!(cstack,assoc[1])
                        in_cluster[assoc[1]] = ccount
                        push!(cstack,assoc[2])
                        in_cluster[assoc[2]] = ccount
                        push!(cstack,assoc[3])
                        in_cluster[assoc[3]] = ccount
                        if flip == true
                            LegType[assoc[1]] =  xor(LegType[assoc[1]],1)
                            LegType[assoc[2]] =  xor(LegType[assoc[2]],1)
                            LegType[assoc[3]] =  xor(LegType[assoc[3]],1)
                        end
                    end #if
                end #if in_cluster == 0
            end #while
        end #if
    end #for i

    #DEBUG
    #for i = 1:lsize[1]
    #    println(i," ",LegType[i]," ",in_cluster[i])
    #end

    #map back basis states and operator list
    ocount = 0
    for i = 1:nSpin
        spin_left[i] = LegType[i]  #left basis state
        ocount += 1
    end

    ocount += 1  #next on is leg nSpin + 1
    for i = 1:2*M
        if operator_list[i,1] != -2 && operator_list[i,1] != -1
            ocount += 4
        else
            if LegType[ocount] == LegType[ocount + 1]  #diagonal
                operator_list[i,1] = -1
                #println("DCHANGE ",i," ",ocount)
            else
                operator_list[i,1] = -2 #off-diagonal
                #println("OCHANGE ",i," ",ocount)
            end
            ocount += 2
        end
    end

    for i = 1:nSpin
        spin_right[i] = LegType[lsize[1] - nSpin + i]  #left basis state
    end

end #ClusterUpdate


text:
using AlmostBandedMatrices

n = 10
r = 2
V = randn(n,r)
A = A_old = SemiSeparableMatrix(LowRankMatrix(randn(n,r), randn(r,r), V),
                    LowRankMatrix(randn(n,r), randn(r,r), V),0,0)

A = deepcopy(A_old)
      A₁ = A - (A[:,1]*A[1,:]')/A[1,1]

A = deepcopy(A_old)
      v = A[:,1];
      v ./= v[1]
      v .*= 1
      v[1] = 0
      L = I - v*[1 zeros(1,n-1)]
      L*A

A = deepcopy(A_old)
      # A.L.Σ .= A.L.Σ*(I - A.L.V[1,:]*(A.U.Σ'*A.U.U[1,:])'/A[1,1])
      A.L.Σ .= A.L.Σ*(I - A.L.V[1,:]*A.U.U[1,:]'*A.U.Σ/A[1,1])
      A
A = deepcopy(A_old)
A.L.Σ*(I - A.L.V[1,:]*(A.U.Σ*A.U.U[1,:])'/A[1,1])

A.L.V[1,:]*A.U.U[1,:]'
A.L.V[1,:]
A.U.U[1,:]
A[2,1]

A.L[2,1]

A[2,1]

A.L[2,1]

A.L.U*A.L.Σ*A.L.V[1,:]

A - (A.L.U*A.L.Σ*A.L.V[1,:]*A[1,:]')/A[1,1]    - A₁
A - (A.L.U*A.L.Σ*A.L.V[1,:]*(A.U.U[1:1,:]*A.U.Σ*A.U.V'))/A[1,1] - A₁
(A.L.U*A.L.Σ*A.L.V') - (A.L.U*A.L.Σ*A.L.V[1,:]*(A.U.V*A.U.Σ'*A.U.U[1,:])')/A[1,1] - A₁
A.L.U*A.L.Σ*(I - A.L.V[1,:]*(A.U.Σ'*A.U.U[1,:])'/A[1,1])*A.U.V'- A₁
A.L.V[1,:]
(A.L.V')

(I - A.L.V[1,:]*(A.U.Σ*A.U.U[1,:])'/A[1,1])

(I - A.L.V[1,:]*(A.U.Σ*A.U.U[1,:])'/A[1,1])

A.L.V[1,:]

(A.U.U[1:1,:]*A.U.Σ*A.U.V')
(A.U.V*A.U.Σ*A.U.U[1,:])
(A.U.U[1:1,:]*A.U.Σ*A.U.V')

A[1,:]
A.L.U

A.L.U*A.L.Σ*A.L.V[1,:]

A = deepcopy(A_old)
      v = A[:,1];
      v ./= v[1]
      v .*= 1
      v[1] = 0
      v
      A.L.V[1,:] ./= A[1,1]
       # A.L.V[2:n,:] .-= (A.U.U[1:1,:] * A.U.Σ * A.U.V[2:n,:]') .* A.L.V[1,1]
       A.L.V[2:n,1] .-= A.U.U[1,1] * A.U.Σ[1] * A.L.V[2:n,1].* A.L.V[1,1]
       A.L.V[2:n,2] .-= A.U.U[1,2] * A.U.Σ[2] * A.L.V[2:n,2].* A.L.V[1,1]
      # A.L.V[2:n,1] .-= A.U.U[1,1] * A.U.Σ[1] * A.L.V[2:n,1].* A.L.V[1,1]
      # A.L.Σ[1] .*= (1 - A.U.U[1,1] * A.U.Σ[1] * A.L.V[1,1])
      A

A = randn(n,n)
      v = A[:,1];
      v ./= v[1]
      v .*= 1
      v[1] = 0
      L = I - v*[1 zeros(1,n-1)]
      L*A

A - (A[:,1]*A[1,:]')/A[1,1]

A - (A[:,1]*[A[1,1] zeros(1,n-1)])/A[1,1]

(A.U.U[1:1,:] * A.U.Σ * A.U.V[2:n,:]')

A.L.V[2:n,:]
A.L.V[1,1]

(A.U.U[1:1,:] * A.U.Σ * A.U.V[2:n,:]')

(A.U.U[1,:]' * A.U.Σ * A.U.V[2:n,:]')
A.U.U[1,:]

A.U.U[1,:]' * A.U.Σ * A.V

L = I - v*[1 zeros(1,n-1)]
      L*A_old

A

A = deepcopy(A_old)
  v = A[:,1];
  v ./= v[1]
  v .*= 1
  v[1] = 0
  A11 = A[1,1]
  A.L.Σ .-= A.U.Σ
  A.L.V[1,1] ./= A[1,1]

  A.L.V

L = I - v*[1 zeros(1,n-1)]
  A - L*A_old

A.U.V

A.L.V[2:n,:] .-= A.U.U[1,1] .* A.L.V[1,1] .* A.U.V[2:n,:]

A = deepcopy(A_old)
   v = A[:,1];
   v ./= v[1]
   v .*= 1
   v[1] = 0
   v
   A.L.V[1,:] ./= A[1,1]
   for j = 2:n
      A.L.V[j,:] .-= A[1,j]*A.L.V[1,:]
   end
   L = I - v*[1 zeros(1,n-1)]
   A - L*A_old

A = deepcopy(A_old)
      v = A[:,1];
      v ./= v[1]
      v .*= 1
      v[1] = 0
      v
      A.L.V[1,:] ./= A[1,1]
      # A.L.V[2:n,:] .-= A.U.U[1,1] * A.U.Σ[1] * A.U.V[2:n,1].* A.L.V[1,1]
      # A.L.V[2:n,1] .-= A.U.U[1,1] * A.U.Σ[1] * A.L.V[2:n,1].* A.L.V[1,1]
      A.L.Σ[1] .*= (1 - A.U.U[1,1] * A.U.Σ[1] * A.L.V[1,1])
      L = I - v*[1 zeros(1,n-1)]
      A - L*A_old

L*A_old

A

A.U.U[1,1] * A.U.Σ[1] .* A.L.V[1,1]

A.L.Σ[1]

A.L.V

A[1,2:n]

A = deepcopy(A_old)
   v = A[:,1];
   v ./= v[1]
   v .*= 1
   v[1] = 0
   v
   A.L.V[1,:] ./= A[1,1]
   A.L.V[2:n,:] .-= A.U.U[1,1] .* A.L.V[1,1] .* A.U.V[2:n,:]
   A
   L = I - v*[1 zeros(1,n-1)]
   A - L*A_old

A.L.V[2:n,:]

A.U.V[:,1]

A.U.U[1,:] .* A.L.V[1,:]
A.L.V[1,:]
A.U.U

v

v

A

j = 1

n = 100;
   u = randn(n); v = randn(n);
   w = randn(n); z = randn(n);

   A = triu(u*v') + tril(w*z')
   v = A[:,1];
   v ./= v[1]
   v[1] = 0
   L = (I - v*[1 zeros(1,length(v)-1)])

ext: jl


text:
# Part 1
parseitem(item) = first(rsplit(item, " "; limit=2)) |> x -> split(x, " "; limit=2)
parsecontents(c) = c[begin:end-1] |> x -> x == "no other bags" ? [] : parseitem.(split(x, ", "))
parserule(rule) = split(rule, " bags contain ") |> r -> r[begin] => parsecontents(r[end])
rules = Dict(parserule.(readlines("day07_input.txt")))
listcontains(list, colour) = (colour ∈ list) || any(map(c -> rulecontains(rules[c], colour), list))
rulecontains(ruleset, colour) = listcontains(map(x -> last(x), ruleset), colour)
println("Part 1: $(sum(rulecontains.(values(rules), "shiny gold"))) bag colours contain ≥1 SGB")

# Part 2
bagcount(colour) = reduce(+, [parse(Int32, first(c)) * bagcount(last(c)) for c in rules[colour]]; init=1)
println("Part 2: SGB contains $(bagcount("shiny gold") - 1) bags")

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
@testset "contact" begin
    @testset "HalfSpace3D" begin
        Random.seed!(4)
        frame = CartesianFrame3D()
        point = Point3D(frame, rand(), rand(), rand())
        normal = FreeVector3D(frame, 0., 0., 1.)
        halfspace = HalfSpace3D(point, normal)

        for i = 1 : 100
            x = Point3D(frame, randn(), randn(), randn())
            @test point_inside(halfspace, x) == (x.v[3] <= point.v[3])
            ϕ, normal = detect_contact(halfspace, x)
            @test ϕ == separation(halfspace, x)
            @test isapprox(normal.v, ForwardDiff.gradient(xyz -> separation(halfspace, Point3D(frame, xyz)), x.v))
        end
    end
end


text:
"""
    couple_models(aero::Steady, stru::TypicalSection)

Create an aerostructural model using a steady aerodynamics model and a
two-degree of freedom typical section model.  This model introduces the
freestream velocity ``U`` and air density ``\\rho`` as additional parameters.
"""
couple_models(aero::Steady, stru::TypicalSection) = (aero, stru)

"""
    couple_models(aero::QuasiSteady, stru::TypicalSection)

Create an aerostructural model using a quasi-steady aerodynamics model and a
two-degree of freedom typical section model.  This model introduces the
freestream velocity ``U`` and air density ``\\rho`` as additional parameters.
"""
couple_models(aero::QuasiSteady, stru::TypicalSection) = (aero, stru)

# --- Traits --- #

number_of_additional_parameters(::Type{QuasiSteady{0}}, ::Type{TypicalSection}) = 2
coupling_inplaceness(::Type{QuasiSteady{0}}, ::Type{TypicalSection}) = OutOfPlace()
coupling_rate_jacobian_type(::Type{QuasiSteady{0}}, ::Type{TypicalSection}) = Zeros()
coupling_state_jacobian_type(::Type{QuasiSteady{0}}, ::Type{TypicalSection}) = Constant()
coupling_parameter_jacobian_type(::Type{QuasiSteady{0}}, ::Type{TypicalSection}) = Nonlinear()
coupling_time_gradient_type(::Type{QuasiSteady{0}}, ::Type{TypicalSection}) = Zeros()

number_of_additional_parameters(::Type{QuasiSteady{1}}, ::Type{TypicalSection}) = 2
coupling_inplaceness(::Type{QuasiSteady{1}}, ::Type{TypicalSection}) = OutOfPlace()
coupling_rate_jacobian_type(::Type{QuasiSteady{1}}, ::Type{TypicalSection}) = Zeros()
coupling_state_jacobian_type(::Type{QuasiSteady{1}}, ::Type{TypicalSection}) = Constant()
coupling_parameter_jacobian_type(::Type{QuasiSteady{1}}, ::Type{TypicalSection}) = Nonlinear()
coupling_time_gradient_type(::Type{QuasiSteady{1}}, ::Type{TypicalSection}) = Zeros()

number_of_additional_parameters(::Type{QuasiSteady{2}}, ::Type{TypicalSection}) = 2
coupling_inplaceness(::Type{QuasiSteady{2}}, ::Type{TypicalSection}) = OutOfPlace()
coupling_rate_jacobian_type(::Type{QuasiSteady{2}}, ::Type{TypicalSection}) = Constant()
coupling_state_jacobian_type(::Type{QuasiSteady{2}}, ::Type{TypicalSection}) = Constant()
coupling_parameter_jacobian_type(::Type{QuasiSteady{2}}, ::Type{TypicalSection}) = Nonlinear()
coupling_time_gradient_type(::Type{QuasiSteady{2}}, ::Type{TypicalSection}) = Zeros()

# --- Methods --- #

function get_coupling_inputs(aero::QuasiSteady{0}, stru::TypicalSection, dx, x, p, t)
    # extract state variables
    h, θ, hdot, θdot = x
    # extract aerodynamic, structural, and aerostructural parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, U, ρ = p
    # local freestream velocity components
    u, v = section_steady_velocities(U, θ)
    # calculate aerodynamic loads
    N, A, M = quasisteady0_loads(a, b, ρ, a0, α0, cd0, cm0, u, v)
    # lift is approximately equal to the normal force
    L = N
    # return inputs
    return SVector(L, M)
end

function get_coupling_inputs(aero::QuasiSteady{1}, stru::TypicalSection, dx, x, p, t)
    # extract state variables
    h, θ, hdot, θdot = x
    # extract aerodynamic, structural, and aerostructural parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, U, ρ = p
    # local freestream velocity components
    u, v, ω = section_velocities(U, θ, hdot, θdot)
    # calculate aerodynamic loads
    N, A, M = quasisteady1_loads(a, b, ρ, a0, α0, cd0, cm0, u, v, ω)
    # lift is approximately equal to the normal force
    L = N
    # return inputs
    return SVector(L, M)
end

function get_coupling_inputs(aero::QuasiSteady{2}, stru::TypicalSection, dx, x, p, t)
    # extract state variables
    dh, dθ, dhdot, dθdot = dx
    # extract state variables
    h, θ, hdot, θdot = x
    # extract aerodynamic, structural, and aerostructural parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, U, ρ = p
    # local freestream velocity components
    u, v, ω = section_velocities(U, θ, hdot, θdot)
    # local freestream accelerations
    udot, vdot, ωdot = section_accelerations(dhdot, dθdot)
    # calculate aerodynamic loads
    N, A, M = quasisteady2_loads(a, b, ρ, a0, α0, cd0, cm0, u, v, ω, vdot, ωdot)
    # lift is approximately equal to the normal force
    L = N
    # return inputs
    return SVector(L, M)
end

# --- Performance Overloads --- #

function get_coupling_state_jacobian(aero::QuasiSteady{0}, stru::TypicalSection, p)
    # extract aerodynamic, structural, and aerostructural parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, U, ρ = p
    # return jacobian
    return quasisteady0_section_state_jacobian(a, b, ρ, a0, U)
end

function get_coupling_parameter_jacobian(aero::QuasiSteady{0}, stru::TypicalSection, dx, x, p, t)

    # extract state variables
    h, θ, hdot, θdot = x
    # extract aerodynamic, structural, and aerostructural parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, U, ρ = p
    # local freestream velocity components
    u, v = section_steady_velocities(U, θ)
    u_U, v_U = section_steady_velocities_U(θ)
    # calculate loads
    N_a, A_a, M_a = quasisteady0_loads_a(a, b, ρ, a0, α0, u, v)
    N_b, A_b, M_b = quasisteady0_loads_b(a, b, ρ, a0, α0, cd0, cm0, u, v)
    N_a0, A_a0, M_a0 = quasisteady0_loads_a0(a, b, ρ, α0, u, v)
    N_α0, A_α0, M_α0 = quasisteady0_loads_α0(a, b, ρ, a0, α0, u, v)
    N_cd0, A_cd0, M_cd0 = quasisteady0_loads_cd0(b, ρ, u)
    N_cm0, A_cm0, M_cm0 = quasisteady0_loads_cm0(b, ρ, u)

    N_kh, A_kh, M_kh = 0, 0, 0
    N_kθ, A_kθ, M_kθ = 0, 0, 0
    N_m, A_m, M_m = 0, 0, 0
    N_Sθ, A_Sθ, M_Sθ = 0, 0, 0
    N_Iθ, A_Iθ, M_Iθ = 0, 0, 0

    N_u, A_u, M_u = quasisteady0_loads_u(a, b, ρ, a0, α0, cd0, cm0, u, v)
    N_v, A_v, M_v = quasisteady0_loads_v(a, b, ρ, a0, α0, u, v)
    N_U = N_u * u_U + N_v * v_U
    M_U = M_u * u_U + M_v * v_U

    N_ρ, A_ρ, M_ρ = quasisteady0_loads_ρ(a, b, ρ, a0, α0, cd0, cm0, u, v)

    # return jacobian
    return @SMatrix [
        N_a N_b N_a0 N_α0 N_cd0 N_cm0 N_kh N_kθ N_m N_Sθ N_Iθ N_U N_ρ;
        M_a M_b M_a0 M_α0 M_cd0 M_cm0 M_kh M_kθ M_m M_Sθ M_Iθ M_U M_ρ
        ]
end

function get_coupling_state_jacobian(aero::QuasiSteady{1}, stru::TypicalSection, p)
    # extract aerodynamic, structural, and aerostructural parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, U, ρ = p
    # return jacobian
    return quasisteady1_section_state_jacobian(a, b, ρ, a0, U)
end

function get_coupling_parameter_jacobian(aero::QuasiSteady{1}, stru::TypicalSection, dx, x, p, t)

    # extract state variables
    h, θ, hdot, θdot = x
    # extract aerodynamic, structural, and aerostructural parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, U, ρ = p
    # local freestream velocity components
    u, v, ω = section_velocities(U, θ, hdot, θdot)
    u_U, v_U, ω_U = section_velocities_U(θ)

    # calculate loads
    N_a, A_a, M_a = quasisteady1_loads_a(a, b, ρ, a0, α0, u, v, ω)
    N_b, A_b, M_b = quasisteady1_loads_b(a, b, ρ, a0, α0, cd0, cm0, u, v, ω)
    N_a0, A_a0, M_a0 = quasisteady1_loads_a0(a, b, ρ, a0, α0, u, v, ω)
    N_α0, A_α0, M_α0 = quasisteady1_loads_α0(a, b, ρ, a0, α0, u, v, ω)
    N_cd0, A_cd0, M_cd0 = quasisteady1_loads_cd0(b, ρ, u)
    N_cm0, A_cm0, M_cm0 = quasisteady1_loads_cm0(b, ρ, u)

    N_kh, A_kh, M_kh = 0, 0, 0
    N_kθ, A_kθ, M_kθ = 0, 0, 0
    N_m, A_m, M_m = 0, 0, 0
    N_Sθ, A_Sθ, M_Sθ = 0, 0, 0
    N_Iθ, A_Iθ, M_Iθ = 0, 0, 0

    N_u, A_u, M_u = quasisteady1_loads_u(a, b, ρ, a0, α0, cd0, cm0, u, v, ω)
    N_v, A_v, M_v = quasisteady1_loads_v(a, b, ρ, a0, α0, u, v, ω)
    N_U = N_u * u_U + N_v * v_U# + L_ω * ω_U
    M_U = M_u * u_U + M_v * v_U# + M_ω * ω_U

    N_ρ, A_ρ, M_ρ = quasisteady1_loads_ρ(a, b, ρ, a0, α0, cd0, cm0, u, v, ω)

    # return jacobian
    return @SMatrix [
        N_a N_b N_a0 N_α0 N_cd0 N_cm0 N_kh N_kθ N_m N_Sθ N_Iθ N_U N_ρ;
        M_a M_b M_a0 M_α0 M_cd0 M_cm0 M_kh M_kθ M_m M_Sθ M_Iθ M_U M_ρ
        ]
end

function get_coupling_rate_jacobian(aero::QuasiSteady{2}, stru::TypicalSection, p)
    # extract aerodynamic, structural, and aerostructural parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, U, ρ = p
    # return jacobian
    return quasisteady2_section_rate_jacobian(a, b, ρ)
end

function get_coupling_state_jacobian(aero::QuasiSteady{2}, stru::TypicalSection, p)
    # extract aerodynamic, structural, and aerostructural parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, U, ρ = p
    # return jacobian
    return quasisteady2_section_state_jacobian(a, b, ρ, a0, U)
end

function get_coupling_parameter_jacobian(aero::QuasiSteady{2}, stru::TypicalSection, dx, x, p, t)
    # extract state variables
    dh, dθ, dhdot, dθdot = dx
    # extract state variables
    h, θ, hdot, θdot = x
    # extract aerodynamic, structural, and aerostructural parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, U, ρ = p
    # local freestream velocity components
    u, v, ω = section_velocities(U, θ, hdot, θdot)
    u_U, v_U, ω_U = section_velocities_U(θ)
    # local freestream accelerations
    udot, vdot, ωdot = section_accelerations(dhdot, dθdot)

    # calculate loads
    N_a, A_a, M_a = quasisteady2_loads_a(a, b, ρ, a0, α0, u, v, ω, vdot, ωdot)
    N_b, A_b, M_b = quasisteady2_loads_b(a, b, ρ, a0, α0, cd0, cm0, u, v, ω, vdot, ωdot)
    N_a0, A_a0, M_a0 = quasisteady2_loads_a0(a, b, ρ, a0, α0, u, v, ω)
    N_α0, A_α0, M_α0 = quasisteady2_loads_α0(a, b, ρ, a0, α0, u, v, ω)
    N_cd0, A_cd0, M_cd0 = quasisteady2_loads_cd0(b, ρ, u)
    N_cm0, A_cm0, M_cm0 = quasisteady2_loads_cm0(b, ρ, u)

    N_kh, A_kh, M_kh = 0, 0, 0
    N_kθ, A_kθ, M_kθ = 0, 0, 0
    N_m, A_m, M_m = 0, 0, 0
    N_Sθ, A_Sθ, M_Sθ = 0, 0, 0
    N_Iθ, A_Iθ, M_Iθ = 0, 0, 0

    N_u, A_u, M_u = quasisteady2_loads_u(a, b, ρ, a0, α0, cd0, cm0, u, v, ω)
    N_v, A_v, M_v = quasisteady2_loads_v(a, b, ρ, a0, α0, u, v, ω)
    N_U = N_u * u_U + N_v * v_U# + L_ω * ω_U
    M_U = M_u * u_U + M_v * v_U# + M_ω * ω_U

    N_ρ, A_ρ, M_ρ = quasisteady2_loads_ρ(a, b, ρ, a0, α0, cd0, cm0, u, v, ω, vdot, ωdot)

    # return jacobian
    return @SMatrix [
        N_a N_b N_a0 N_α0 N_cd0 N_cm0 N_kh N_kθ N_m N_Sθ N_Iθ N_U N_ρ;
        M_a M_b M_a0 M_α0 M_cd0 M_cm0 M_kh M_kθ M_m M_Sθ M_Iθ M_U M_ρ
        ]
end

# --- Convenience Methods --- #

function set_additional_parameters!(padd, aero::QuasiSteady, stru::TypicalSection; U, rho)

    padd[1] = U
    padd[2] = rho

    return padd
end

function separate_additional_parameters(aero::QuasiSteady, stru::TypicalSection, padd)

    return (U = padd[1], rho = padd[2])
end

# --- Plotting --- #

@recipe function f(aero::QuasiSteady, stru::TypicalSection, dx, x, y, p, t)

    framestyle --> :origin
    grid --> false
    xlims --> (-1.0, 1.0)
    ylims --> (-1.5, 1.5)
    label --> @sprintf("t = %6.3f", t)

    # extract state variables
    h, θ, hdot, θdot = x

    # extract parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, U, ρ = p

    xplot = [-(0.5 + a*b)*cos(θ),    (0.5 - a*b)*cos(θ)]
    yplot = [ (0.5 + a*b)*sin(θ)-h, -(0.5 - a*b)*sin(θ)-h]

    return xplot, yplot
end

# --- Internal Methods --- #

function section_steady_velocities(U, θ)

    u = U
    v = U*θ

    return SVector(u, v)
end

function section_steady_velocities_U(θ)

    u_U = 1
    v_U = θ

    return SVector(u_U, v_U)
end

function section_steady_velocities_θ(U)

    u_θ = 0
    v_θ = U

    return SVector(u_θ, v_θ)
end

function quasisteady0_section_state_jacobian(a, b, ρ, a0, U)
    L_θ = a0*ρ*U^2*b
    M_θ = (b/2 + a*b)*L_θ
    return @SMatrix [0 L_θ 0 0; 0 M_θ 0 0]
end

function quasisteady1_section_state_jacobian(a, b, ρ, a0, U)
    tmp1 = a0*ρ*U*b
    tmp2 = pi*ρ*b^3
    d1 = b/2 - a*b
    d2 = b/2 + a*b
    L_θ = tmp1*U
    L_hdot = tmp1
    L_θdot = tmp1*d1 + tmp2*U/b
    M_θ = d2*L_θ
    M_hdot = d2*L_hdot
    M_θdot = -tmp2*U + d2*L_θdot
    return @SMatrix [0 L_θ L_hdot L_θdot; 0 M_θ M_hdot M_θdot]
end

quasisteady2_section_state_jacobian(a, b, ρ, a0, U) = quasisteady1_section_state_jacobian(a, b, ρ, a0, U)

function quasisteady2_section_rate_jacobian(a, b, ρ)
    N_hddot, A_hddot, M_hddot = quasisteady2_loads_vdot(a, b, ρ)
    N_θddot, A_θddot, M_θddot = quasisteady2_loads_ωdot(a, b, ρ)
    return @SMatrix [0 0 N_hddot N_θddot; 0 0 M_hddot M_θddot]
end


text:
# module Lupin

function hreshape(V::Vector{T}; dims) where T
    A = Array{T}(undef, dims)
    (m, n) = dims
    for i=axes(A, 1), j=axes(A, 2)
        A[i,j] = V[(i-1)*n+j]
    end
    A
end

# similar Base.Iterators.partition but returns Matrix
function partition(v::Vector{T}; cols::Int, fill)::Matrix{T} where T
    len = length(v)
    n = div(len, cols)
    if iszero(len % cols)
        hreshape(v, dims=(n, cols))
    else
        A = hreshape(v[1:n*cols], dims=(n, cols))
        vcat(A, [v[n*cols+1:end]... Base.fill(fill, cols - len % cols)...])
    end
end

function partition(r::AbstractRange; kwargs...)
    partition(collect(r); kwargs...)
end

# partition by condition
function partition(f::Function, array::Array{T})::Tuple{Array{T}, Array{T}} where {T}
    (getindex(array, findall(f, array)),
     getindex(array, findall(!f, array)))
end

function partition(f::Function, r::AbstractRange{T})::Tuple{Vector{T}, Vector{T}} where {T}
    (getindex(r, findall(f, r)),
     getindex(r, findall(!f, r)))
end

# module Lupin


text:
@testset "Periodic extrapolation" begin
    xmax = 7
    f(x) = sin((x-3)*2pi/xmax - 1)
    A = Float64[f(x) for x in 1:xmax] # Does not include the periodic image

    itp0 = interpolate(A, BSpline(Constant(Periodic())))
    itp1 = interpolate(A, BSpline(Linear(Periodic())))
    itp2 = interpolate(A, BSpline(Quadratic(Periodic(OnCell()))))
    itp3 = interpolate(A, BSpline(Cubic(Periodic(OnCell()))))

    etp0 = extrapolate(itp0, Periodic())
    etp1 = extrapolate(itp1, Periodic())
    etp2 = extrapolate(itp2, Periodic())
    etp3 = extrapolate(itp3, Periodic())

    for x in -xmax:.4:2*xmax
        @test etp0(x) ≈ f(x) atol=0.5
        @test etp1(x) ≈ f(x) atol=0.1
        @test etp2(x) ≈ f(x) atol=0.01
        @test etp3(x) ≈ f(x) atol=0.003
    end
end


text:
module DeSolver

#Dependencies
using Plots; gr()

include("functions.jl")
include("odemethods.jl")

export odesolve
export RK4
end # module


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -13.110901270484293 <= q <= 11.380819522292354)

                     add_NL_constraint(m, :(gelu(-0.7773677229963925 + -0.31823111196481735*gelu(0.13898320039278467 + -0.4573031975765818*gelu(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + -0.027363793963717686*gelu(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + -0.7118578003946046*gelu(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))) + -0.05869662098445039*gelu(0.7177153834386609 + 0.046197114005355644*gelu(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + 0.8341427099559482*gelu(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + -0.5961814570241306*gelu(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))) + -0.31920840200501477*gelu(-0.05621054225465549 + -0.4685096496776473*gelu(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + -0.723697578017108*gelu(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + 0.643520408356284*gelu(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3])))) + gelu(-0.2800668423712529 + 0.5241807567635033*gelu(0.13898320039278467 + -0.4573031975765818*gelu(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + -0.027363793963717686*gelu(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + -0.7118578003946046*gelu(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))) + 0.9218858277065758*gelu(0.7177153834386609 + 0.046197114005355644*gelu(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + 0.8341427099559482*gelu(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + -0.5961814570241306*gelu(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))) + 0.0068491484591501894*gelu(-0.05621054225465549 + -0.4685096496776473*gelu(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + -0.723697578017108*gelu(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + 0.643520408356284*gelu(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3])))) + gelu(-0.0923150493931284 + -0.18688964829603716*gelu(0.13898320039278467 + -0.4573031975765818*gelu(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + -0.027363793963717686*gelu(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + -0.7118578003946046*gelu(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))) + 0.4699561510424055*gelu(0.7177153834386609 + 0.046197114005355644*gelu(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + 0.8341427099559482*gelu(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + -0.5961814570241306*gelu(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))) + -0.26284635568255643*gelu(-0.05621054225465549 + -0.4685096496776473*gelu(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + -0.723697578017108*gelu(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + 0.643520408356284*gelu(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Crux
using Test
using POMDPModels

## Spaces
mdp = SimpleGridWorld()

s1 = DiscreteSpace(5)
@test s1 isa AbstractSpace
@test s1 isa DiscreteSpace
@test s1.N == 5
@test type(s1) == Bool
@test Crux.dim(s1) == (5,)
@test s1.vals == [1,2,3,4,5]
@test tovec(4, s1) == Bool[0,0,0,1,0]

s2 = DiscreteSpace(5, [:a, :b, :c, :d, :e])
@test s2 isa AbstractSpace
@test s2 isa DiscreteSpace
@test s2.N == 5
@test type(s2) == Bool
@test Crux.dim(s2) == (5,)
@test s2.vals == [:a, :b, :c, :d, :e]
@test tovec(:d, s2) == Bool[0,0,0,1,0]

s3 = ContinuousSpace((3,4), Float64)
@test s3 isa AbstractSpace
@test s3 isa ContinuousSpace
@test s3.dims == (3,4)
@test type(s3) == Float64
@test Crux.dim(s3) == (3,4)
@test tovec(zeros(3,4), s3) == zeros(3, 4)

s3 = ContinuousSpace((3,4), Float64, 1f0, 2f0)
@test s3 isa AbstractSpace
@test s3 isa ContinuousSpace
@test s3.dims == (3,4)
@test type(s3) == Float64
@test Crux.dim(s3) == (3,4)
@test tovec(zeros(3,4),s3) == -0.5*ones(3, 4)

s4 = state_space(mdp)
@test s4 isa ContinuousSpace
@test Crux.dim(s4) == (2,)


text:
function build_astar_path(came_from::Dict, t::Integer, v::Integer, t0::Integer)
    path = Path()
    (τ, u) = (t, v)
    pushfirst!(path, (τ, u))
    while τ > t0
        (τ, u) = came_from[(τ, u)]
        pushfirst!(path, (τ, u))
    end
    return path
end

function temporal_astar(
    g::AbstractGraph{V},
    s::Integer,
    d::Integer,
    t0::Integer;
    edge_indices::Dict,
    edge_weights::Vector{W},
    heuristic=v -> 0.0,
    reservation::Reservation=Reservation(),
    conflict_price=Inf,
) where {V,W}
    T = Int
    safe_conflict_price = conflict_price * (conflict_price < Inf)

    # Initialize storage
    came_from = Dict{Tuple{T,V},Tuple{T,V}}()
    distance = Dict{Tuple{T,V},W}()
    conflicts = Dict{Tuple{T,V},Int}()
    queue = PriorityQueue{Tuple{T,V},Float64}()

    # Add first node to storage
    distance_s = zero(W)
    conflicts_s = Int((t0, s) in reservation)
    if conflicts_s == 0 || conflict_price < Inf
        priority_s = safe_conflict_price * conflicts_s + heuristic(s)
        distance[(t0, s)] = distance_s
        conflicts[(t0, s)] = conflicts_s
        enqueue!(queue, (t0, s), priority_s)
    end

    # Explore
    nodes_explored = 0
    while !isempty(queue)
        (t, v) = dequeue!(queue)
        nodes_explored += 1
        (v == d) && return build_astar_path(came_from, t, v, t0)

        for w in outneighbors(g, v)
            heur_w = heuristic(w)

            e_vw = edge_indices[(v, w)]
            weight_vw = edge_weights[e_vw]
            conflict_vw = Int((t + 1, w) in reservation)

            distance_w = distance[(t, v)] + weight_vw
            conflicts_w = conflicts[(t, v)] + conflict_vw

            if conflicts_w == 0 || conflict_price < Inf
                old_distance_w = get(distance, (t + 1, w), Inf)
                old_conflicts_w = get(conflicts, (t + 1, w), typemax(Int) ÷ 2)

                old_cost_w = (safe_conflict_price * old_conflicts_w + old_distance_w)
                cost_w = (safe_conflict_price * conflicts_w + distance_w)

                if cost_w < old_cost_w
                    came_from[(t + 1, w)] = (t, v)
                    distance[(t + 1, w)] = distance_w
                    conflicts[(t + 1, w)] = conflicts_w
                    queue[(t + 1, w)] = cost_w + heur_w
                end
            end
        end
    end
    return Path()
end


text:
using PyCall, Conda

JAXLIB="0.1.46"

if PyCall.conda
  Conda.add("pip")
  pip = joinpath(Conda.BINDIR, "pip")
  run(`$pip install jaxlib==$JAXLIB`)
else
  try
    pyimport("jaxlib")
  catch e
    e isa PyCall.PyError || rethrow(e)
    error("""
        Python Dependencies not installed
        Please either:
         - Rebuild PyCall to use Conda, by running in the julia REPL:
           - `using Pkg; ENV["PYTHON"]=""; Pkg.build("PyCall"); Pkg.build("XLA")`
         - Or install the depencencies, eg by running pip
        	 - `pip install jaxlib==$JAXLIB`
      	""")
  end
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Private functions and macros.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

const _html_alignment = Dict(
    :l => "left",
    :L => "left",
    :c => "center",
    :C => "center",
    :r => "right",
    :R => "right"
)

function _html_text_alignment_dict(alignment::Symbol)
    if (alignment == :n) || (alignment == :N)
        return Dict{String, String}()
    elseif haskey(_html_alignment, alignment)
        return Dict{String, String}("text-align" => _html_alignment[alignment])
    else
        return Dict{String, String}("text-align" => _html_alignment[:r])
    end
end

function _styled_html(
    tag::String,
    text::String,
    style::Dict{String,String} = Dict{String,String}();
    class::String = ""
)
    str_class = isempty(class) ? "" : " class = \"" * class * "\""

    # If there is no keys in the style dictionary, just return the tag.
    if isempty(style)
        return "<" * tag * str_class * ">" * text * "</" * tag * ">"
    else
        # Create the style string.
        style_str = ""

        # We must sort the keys so that we can provide stable outputs.
        v   = collect(values(style))
        k   = collect(keys(style))
        ind = sortperm(collect(keys(style)))
        vk  = @view k[ind]
        vv  = @view v[ind]

        num_styles = length(vk)

        @inbounds for i = 1:num_styles
            # If the value is empty, then just continue.
            value = vv[i]
            isempty(value) && continue

            style_str *= vk[i] * ": "
            style_str *= value * ";"

            i != num_styles && (style_str *= " ")
        end

        return "<" * tag * str_class * " style = \"" * style_str * "\">" * text * "</" * tag * ">"
    end
end


text:
function pxycat(pxy)
    # concat parsing results from multiple cpus
    # input [(p1,x1,y1), (p2,x2,y2), ...]
    # output (p, x, y)
    (mycat(map(z->z[1], pxy)),
     mycat(map(z->z[2], pxy)),
     mycat(map(z->z[3], pxy)))
end

function mycat(a::Array)
    # concat elements of a along their last dimension
    mx = size(a[1])[1:end-1]
    nx = 0
    for x in a
        @assert size(x)[1:end-1] == mx
        nx += size(x)[end]
    end
    b = similar(a[1], tuple(mx..., nx))
    nb = 0
    for x in a
        copy!(b, nb+1, x, 1, length(x))
        nb += length(x)
    end
    return b
end


text:
# generate benchmarks for documentation
function generate_docbenchmark()
    a = tensor_dict[(3, 1, Float64)]
    A = tensor_dict[(3, 2, Float64)]
    As = symtensor_dict[(3, 2, Float64)]
    AA = tensor_dict[(3, 4, Float64)]
    AAs = symtensor_dict[(3, 4, Float64)]

    # arrays and voigt forms
    b   = Array(a)
    B   = Array(A);   Bv = reshape(B, (9, ))
    Bs  = Array(As);  Bsv = reshape(Bs, (9, ))
    BB  = Array(AA);  BBv = reshape(BB, (9, 9))
    BBs = Array(AAs); BBsv = reshape(BBs, (9, 9))

    open(joinpath(dirname(@__FILE__), "../docs/src/benchmarks.md"), "w") do file
        function printheader(head)
            println(file, "| **$(head)** | | | |")
        end
        function printrow(op)
            pretty = (t) -> BenchmarkTools.prettytime(BenchmarkTools.time(minimum(t)))
            speedup = (ta, tt) -> round(10*BenchmarkTools.time(minimum(ta))/BenchmarkTools.time(minimum(tt)))/10
            println(file, "| `$(op)` | $(pretty(tt)) | $(pretty(ta)) | ×$(speedup(ta, tt)) |")
        end

        print(file, """
                    # Benchmarks

                    Here are some benchmark timings for tensors in 3 dimensions. For comparison
                    the timings for the same operations using standard Julia `Array`s are also
                    presented.

                    In the table below, `a` denotes a vector, `A`, `As` denotes second order
                    non-symmetric and symmetric tensors and `AA`, `AAs` denotes fourth order
                    non-symmetric and symmetric tensors respectively.

                    | Operation  | `Tensor` | `Array` | speed-up |
                    |:-----------|---------:|--------:|---------:|
                    """)

        printheader("Single contraction")
        tt = @benchmark $a ⋅ $a
        ta = @benchmark $b ⋅ $b
        printrow("a ⋅ a")

        tt = @benchmark $A ⋅ $a
        ta = @benchmark $B * $b
        printrow("A ⋅ a")

        tt = @benchmark $A ⋅ $A
        ta = @benchmark $B * $B
        printrow("A ⋅ A")

        tt = @benchmark $As ⋅ $As
        ta = @benchmark $Bs * $Bs
        printrow("As ⋅ As")

        printheader("Double contraction")
        tt = @benchmark $A ⊡ $A
        ta = @benchmark $B ⋅ $B
        printrow("A ⊡ A")

        tt = @benchmark $As ⊡ $As
        ta = @benchmark $Bsv ⋅ $Bsv
        printrow("As ⊡ As")

        tt = @benchmark $AA ⊡ $A
        ta = @benchmark $BBv * $Bv
        printrow("AA ⊡ A")

        tt = @benchmark $AA ⊡ $AA
        ta = @benchmark $BBv * $BBv
        printrow("AA ⊡ AA")

        tt = @benchmark $AAs ⊡ $AAs
        ta = @benchmark $BBsv * $BBsv
        printrow("AAs ⊡ AAs")

        tt = @benchmark $As ⊡ $AAs ⊡ $As
        ta = @benchmark $Bsv ⋅ ($BBsv * $Bsv)
        printrow("As ⊡ AAs ⊡ As")

        printheader("Outer product")
        tt = @benchmark $a ⊗ $a
        ta = @benchmark $b * $(b)'
        printrow("a ⊗ a")

        tt = @benchmark $A ⊗ $A
        ta = @benchmark $Bv * $(Bv)'
        printrow("A ⊗ A")

        tt = @benchmark $As ⊗ $As
        ta = @benchmark $Bsv * $(Bsv)'
        printrow("As ⊗ As")

        printheader("Other operations")
        tt = @benchmark det($A)
        ta = @benchmark det($B)
        printrow("det(A)")

        tt = @benchmark det($As)
        ta = @benchmark det($Bs)
        printrow("det(As)")

        tt = @benchmark inv($A)
        ta = @benchmark inv($B)
        printrow("inv(A)")

        tt = @benchmark inv($As)
        ta = @benchmark inv($Bs)
        printrow("inv(As)")

        tt = @benchmark norm($a)
        ta = @benchmark norm($b)
        printrow("norm(a)")

        tt = @benchmark norm($A)
        ta = @benchmark norm($(Bv))
        printrow("norm(A)")

        tt = @benchmark norm($As)
        ta = @benchmark norm($(Bsv))
        printrow("norm(As)")

        tt = @benchmark norm($AA)
        ta = @benchmark norm($(BBv[:]))
        printrow("norm(AA)")

        tt = @benchmark norm($AAs)
        ta = @benchmark norm($(BBsv[:]))
        printrow("norm(AAs)")

        tt = @benchmark $a × $a
        ta = @benchmark $b × $b
        printrow("a × a")

        # print versioninfo
        println(file, """

                      The benchmarks are generated by
                      [`benchmark_doc.jl`](https://github.com/Ferrite-FEM/Tensors.jl/blob/master/benchmark/benchmark_doc.jl)
                      on the following system:

                      ```
                      julia> versioninfo()
                      """)
        versioninfo(file)
        println(file, "```")
    end
end


text:
module logistic_regression
export sigmoid, h, J, init_theta, grad, accuracy

sigmoid(z) = 1 ./ (1 + exp(-z))
h(theta, X) = sigmoid(theta' * X)

function J(theta, X, y)
    hh = h(theta, X)
    return sum(-y .* log(hh) - (1-y) .* log(1-hh)) / size(X, 2)
end

init_theta(n) = randn(n, 1) * 0.01
grad(theta, X, y) = X * (h(theta, X) - y)' / size(X,2)
accuracy(theta, X, y) = size(y[:,(h(theta, X) .> 0.5) .== (y .> 0.5)], 2) / size(y, 2)
end

max_stars_repo_name: ilkerkesen/ufldl-tutorial

max_issues_repo_name: ilkerkesen/ufldl-tutorial

max_forks_repo_name: ilkerkesen/ufldl-tutorial


text:
# Import common packages and functions
include("_obj_func_definition.jl")

"""Example calibration function.

Illustrate model calibration using the BlackBoxOptim package.
"""
function calibrate(sn, v_id, climate, calib_data)

    # Fitness of model is dependent on next node.
    ins = inlets(sn, v_id)
    next_node_id = outlets(sn, v_id)[1]

    # Recurse through and calibrate all nodes upstream
    if !isempty(ins)
        for nid in ins
            calibrate(sn, nid, climate, calib_data)
        end
    end

    this_node = sn[v_id]

    # Create new optimization function (see definition inside `_obj_func_definition.jl`)
    opt_func = x -> obj_func(x, climate, sn, v_id, next_node_id, calib_data)

    # Get node parameters (default values and bounds)
    p_names, x0, param_bounds = param_info(this_node; with_level=false)
    opt = bbsetup(opt_func; SearchRange=param_bounds,
                  Method=:adaptive_de_rand_1_bin_radiuslimited,
                  MaxTime=2400.0,  # time in seconds to spend
                  TraceInterval=30.0,
                  PopulationSize=75,
                  )

    res = bboptimize(opt)

    bs = best_candidate(res)
    @info "Calibrated $(v_id) ($(this_node.name)), with score: $(best_fitness(res))"
    @info "Best Params:" collect(bs)

    # Update node with calibrated parameters
    update_params!(this_node, bs...)

    return res, opt
end

v_id, node = sn["406219"]
@info "Starting calibration..."
res, opt = calibrate(sn, v_id, climate, hist_data)

# Stream

best_params = best_candidate(res)

@info best_fitness(res)
@info best_params

using Plots

update_params!(node, best_params...)
dam_id, dam_node = sn["406000"]
Streamfall.run_node!(sn, dam_id, climate; extraction=hist_dam_releases)

h_data = hist_data["406000"]
n_data = dam_node.level

nnse_score = Streamfall.NNSE(h_data, n_data)
nse_score = Streamfall.NSE(h_data, n_data)
rmse_score = Streamfall.RMSE(h_data, n_data)

@info "Downstream Dam Level NNSE:" nnse_score
@info "Downstream Dam Level RMSE:" rmse_score

reset!(dam_node)

nse = round(nse_score, digits=4)
rmse = round(rmse_score, digits=4)

plot(h_data,
     legend=:bottomleft,
     title="Calibrated IHACRES\n(NSE: $(nse); RMSE: $(rmse))",
     label="Historic", xlabel="Day", ylabel="Dam Level [mAHD]")

plot!(n_data, label="IHACRES")

savefig("calibration_ts_comparison.png")

# 1:1 Plot
scatter(h_data, n_data, legend=false,
        markerstrokewidth=0, markerstrokealpha=0, alpha=0.2)
plot!(h_data, h_data, color=:red, markersize=.1, markerstrokewidth=0,
      xlabel="Historic [mAHD]", ylabel="IHACRES [mAHD]", title="Historic vs Modelled")

savefig("calibration_1to1.png")

# NNSE: 0.9643; RMSE: 1.43553
# d: 84.28015146853407
# d2: 2.4224106535469145
# e: 0.8129590022893607
# f: 2.579276454391652
# a: 5.923379062122229
# b: 0.0989925603647026
# storage_coef: 1.8613364808233752  # gw storage factor
# alpha: 0.7279050097363565


text:
"""
Module for visualizing water networks, and potentially other graphics
"""
module WaterModelsAnalytics
# imports
import WaterModels
const _WM = WaterModels
const _IM = WaterModels._IM

# importing only the viridis color scheme for now; might allow the color scheme to be
# selected by the user at some point, in which case will need to import all of ColorSchemes
import ColorSchemes.viridis
import ColorTypes.HSV
import Memento
import Statistics.mean
using Printf

import Interpolations
const _ITP = Interpolations

import Plots
import DataFrames
using LaTeXStrings

# Python imports.
import PyCall
import PyCall: PyObject, @pycall
const pgv = PyCall.PyNULL()
const wntr = PyCall.PyNULL()
const wntrctrls = PyCall.PyNULL()
const wntr_vis = PyCall.PyNULL()
const stack_cbar = PyCall.PyNULL()
const collate_viz = PyCall.PyNULL()

# Create our module-level logger (this will get precompiled).
const _LOGGER = Memento.getlogger(@__MODULE__)

function __init__()
    # alert about warnings generated during the build of the package (warnings are not shown
    # in the REPL)
    _build_log = joinpath(dirname(dirname(@__FILE__)), "deps", "build.log")

    if isfile(_build_log) && occursin("Warning:", read(_build_log, String))
        @warn("Warnings were generated during the last build of WaterModelsAnalytics.jl:  please check the build log at $_build_log")
    end

    # Register the module-level logger at runtime so users can access the logger via
    # `getlogger(WaterModelsAnalytics)` NOTE: If this line is not included, then the
    # precompiled `WaterModelsAnalytics._LOGGER` will not be registered at runtime.
    Memento.register(_LOGGER)

    try
        # Import WNTR-related components.
        PyCall.pyimport("wntr")
        copy!(wntr, PyCall.pyimport("wntr"))
        copy!(wntrctrls, PyCall.pyimport("wntr.network.controls"))
        wntr_vis_path = joinpath(dirname(pathof(@__MODULE__)), "python")

        PyCall.py"""
        import sys; sys.path.insert(0, $(wntr_vis_path))
        """

        copy!(wntr_vis, PyCall.pyimport("wntr_vis"))
        copy!(stack_cbar, wntr_vis.stack_cbar)
        copy!(collate_viz, wntr_vis.collate_viz)
    catch
        error("Python installation is missing the \"wntr\" module.")
    end

    try
        # Import pygraphviz-related components.
        PyCall.pyimport("pygraphviz")
        copy!(pgv, PyCall.pyimport("pygraphviz"))
    catch
        error("Python installation is missing the \"pygraphviz\" module.")
    end
end

"Suppresses information and warning messages output by WaterModelsAnalytics. For more
        fine-grained control, use the Memento package."
function silence()
    Memento.info(_LOGGER, "Suppressing information and warning messages for "
                 * "the rest of this session. Use the Memento package for more "
                 * "fine-grained control of logging.")
    Memento.setlevel!(Memento.getlogger(_IM), "error")
    Memento.setlevel!(Memento.getlogger(_WM), "error")
end

"Allows the user to set the logging level without the need to add Memento."
function logger_config!(level)
    Memento.config!(Memento.getlogger("WaterModelsAnalytics"), level)
end

include("analysis/utility.jl")
include("analysis/simulation.jl")
include("analysis/validation.jl")
include("analysis/visualization.jl")
include("analysis/pump_bep.jl")

include("graph/graph.jl")
include("graph/output.jl")

include("plots/pumps.jl")

export build_graph
export write_graph
export colorbar
export write_visualization
export write_multi_time_viz

export initialize_wntr_network

export initialize_wntr_network
export update_wntr_controls
export simulate_wntr

export get_node_dataframe
export get_tank_dataframe
export get_pipe_dataframe
export get_short_pipe_dataframe
export get_valve_dataframe
export get_pump_dataframe

# export compare_tank_head
export compare_tank_level

export calc_pump_bep!
export plot_pumps

end # module WaterModelsAnalytics


text:
using BayesianDataFusion
using Base.Test

using Compat

rows = [1:3; 2:4; 1:4]
cols = [1,1,1, 2,2,2, 3,3,3,3]
m = SparseBinMatrix(rows, cols)
@test size(m) == (4,3)

m2 = m[ [true, false, true, false], : ]
@test size(m2) == (2,3)
@test length(m2.rows) == 5
rows2 = [1:2; 2; 1:2]
cols2 = [1,1, 2, 3,3]
@test m2.rows == rows2
@test m2.cols == cols2

## test with random matrix
A = sprand(100, 50, 0.2)
A.nzval[:] = 1.0
I,J,V = findnz(A)
Asbm = SparseBinMatrix(I, J)
z = zeros(Bool, size(A,1))
z[1:20] = true
z[40]   = true
z[60:80]= true

Az = A[z,:]
Asbmz = Asbm[z,:]
Az2 = sparse(Asbmz.rows, Asbmz.cols, 1.0)
@test vecnorm(Az - Az2) < 1e-7


text:
struct NACA4{T<:AbstractFloat, SYM, STE} <: NACABase{T}
    m::T
    p::T
    t::T
end

function NACA4(::Type{T}, num::AbstractString, sharp_TE::Bool) where {T<:AbstractFloat}
    if length(num) != 4
        throw(ArgumentError("invalid 4-dig NACA airfoil: NACA$num"))
    end
    t = parse(T, num[3:4]) / 100
    m = parse(T, num[1:1]) / 100
    p = parse(T, num[2:2]) / 10
    NACA4{T, m==0, sharp_TE}(m,p,t)
end

function y_t(a::NACA4{T, SYM, false}, x::T)::T where {T<:AbstractFloat,SYM}
    5*a.t*(.2969*√x-.1260*x-.3516*x^2+.2843*x^3-.1015*x^4)
end
function y_t(a::NACA4{T, SYM, true}, x::T)::T where {T<:AbstractFloat,SYM}
    5*a.t*(.2969*√x-.1260*x-.3516*x^2+.2843*x^3-.1036*x^4)
end

function dy_t(a::NACA4{T, SYM, false}, x::T)::T where {T<:AbstractFloat,SYM}
    5*a.t*(.2969/(2*√x)-.1260-.3516*2*x+.2843*3*x^2-.1015*4*x^3)
end

function dy_t(a::NACA4{T, SYM, true}, x::T)::T where {T<:AbstractFloat,SYM}
    5*a.t*(.2969/(2*√x)-.1260-.3516*2*x+.2843*3*x^2-.1036*4*x^3)
end

function y_c(a::NACA4{T}, x::T)::T where {T<:AbstractFloat}
    if x < a.p
        a.m/a.p^2 * (2*a.p*x - x^2)
    else
        a.m/(1-a.p)^2 * (1 - 2*a.p + 2*a.p*x -x^2)
    end
end

function dy_c(a::NACA4{T}, x::T)::T where {T<:AbstractFloat}
    if x < a.p
        2 * a.m/a.p^2 * (a.p - x)
    else
        2 * a.m/(1-a.p)^2 * (a.p*x -x)
    end
end

function y_upper(a::NACA4{T,true}, x::T)::T where {T<:AbstractFloat}
    y_t(a, x)
end

function y_lower(a::NACA4{T,true}, x::T)::T where {T<:AbstractFloat}
    -y_t(a, x)
end

function dy_upper(a::NACA4{T,true}, x::T)::T where {T<:AbstractFloat}
    dy_t(a, x)
end

function dy_lower(a::NACA4{T,true}, x::T)::T where {T<:AbstractFloat}
    -dy_t(a, x)
end

function fy_upper(a::NACA4{T,true}, x::T)::T where {T<:AbstractFloat}
    y_t(a, x), dy_t(a, x)
end

function fy_lower(a::NACA4{T,true}, x::T)::T where {T<:AbstractFloat}
    -y_t(a, x), -dy_t(a, x)
end

function y_upper(a::NACA4{T,false}, x::T)::T where {T<:AbstractFloat}
    y_c(a,x) + y_t(a, x) * cos(atan(dy_c(a,x)))
end

function y_lower(a::NACA4{T,false}, x::T)::T where {T<:AbstractFloat}
    y_c(a,x) - y_t(a, x) * cos(atan(dy_c(a,x)))
end

function x_upper(a::NACA4{T,false}, x::T)::T where {T<:AbstractFloat}
    x - y_t(a, x) * sin(atan(dy_c(a,x)))
end

function x_lower(a::NACA4{T,false}, x::T)::T where {T<:AbstractFloat}
    x + y_t(a, x) * sin(atan(dy_c(a,x)))
end


text:
@testset "Opening models from BIGG" begin

    sbmlmodel = load_model(model_paths["iJO1366.xml"])
    @test sbmlmodel isa SBMLModel
    @test n_reactions(sbmlmodel) == 2583

    matlabmodel = load_model(model_paths["iJO1366.mat"])
    @test matlabmodel isa MATModel
    @test n_reactions(matlabmodel) == 2583

    jsonmodel = load_model(model_paths["iJO1366.json"])
    @test jsonmodel isa JSONModel
    @test n_reactions(jsonmodel) == 2583

    @test Set(lowercase.(reactions(sbmlmodel))) ==
          Set("r_" .* lowercase.(reactions(matlabmodel)))
    @test Set(lowercase.(reactions(sbmlmodel))) ==
          Set("r_" .* lowercase.(reactions(jsonmodel)))
end


text:
"""
    TypeG()

Lorem ipsum dolor sit amet.
"""
struct TypeG

end


text:
struct CouplingSampler{Tprob<:Real, PersonIdx<:Integer, GroupIdx<:Integer}
  group_ids::Vector{GroupIdx} # TODO: needed only for reverse grouping computation, forward can be computed on-the-fly
  grouping::PopulationGrouping{PersonIdx}
  matrix_sampler::MatrixAliasSampler{GroupIdx, Float64}
end

function CouplingSampler(group_ids::AbstractVector{GroupIdx}, coupling_weights::Matrix{T} where T<:Real, PersonIdx::DataType=UInt32) where GroupIdx<:Integer
  num_groups, M = size(coupling_weights)
  @assert num_groups == M

  @assert maximum(group_ids) <= num_groups <= typemax(GroupIdx)
  @assert minimum(group_ids) > 0

  @assert maximum(coupling_weights) < Inf
  @assert minimum(coupling_weights) >= 0

  grouping = PopulationGrouping(group_ids, num_groups, PersonIdx)
  matrix_sampler = MatrixAliasSampler(coupling_weights, GroupIdx)
  CouplingSampler{eltype(coupling_weights), PersonIdx, GroupIdx}(group_ids, grouping, matrix_sampler)
end

groupof(s::CouplingSampler, person_idx::Integer) = s.group_ids[person_idx]
getgroup(coupling::CouplingSampler, group_idx::Integer) = getgroup(coupling.grouping, group_idx)

function sample(rng::AbstractRNG, coupling::CouplingSampler, source_id::Integer)
  # retrieve the source group
  source_group_id = groupof(coupling, source_id)

  # alias sample the group id
  target_group_id = sample(rng, coupling.matrix_sampler, source_group_id)

  # retrieve the group
  target_group = getgroup(coupling.grouping, target_group_id)
  @assert length(target_group) > 0 "empty population group $target_group_id, grouping = $(coupling.grouping)"

  # sample from the group
  rand(rng, target_group)
end


text:
using Test

function check(n)
    d = reverse(digits(n))
    increasing = true
    double = false
    for i in 1:length(d)-1
        increasing = increasing && d[i] <= d[i+1]
        double = double || d[i] == d[i+1]
    end
    return increasing && double
end

function run(s)
    # println(s)
    n1, n2 = parse.(Int, (split(s, "-")))
    total = 0
    return sum(check.(n1:n2))
end

@test check(111111) == true
@test check(113456) == true
@test check(223450) == false
@test check(123789) == false

##############################

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time*1000)")
    println(res)
end

main()


text: Dict{String, Any}("table" => Dict{String, Any}())


text:
# (coeff, featurenum, lag, windowlen) for each factor and possibly an intercept.
"""
A moving average linear regression model for time series prediction. This is
similar to but more general than a lagged factor model since the factors can
be moving averages over a range of lags rather than simply a single lag.
In theory, this could "smoothen out" noise in the dependence from one feature
to another while still allowing non-moving average (lagged factors) by setting
the length of the moving window to 1. The `minlag` and `maxlag` limits
how far back in time we can look.
"""
mutable struct MovingAverageLITRM <: AbstractLinearITRM
    nfactors::Int
    nfeatures::Int
    minlag::Int
    maxlag::Int
    maxcoef::Float64
    mincoef::Float64
    hasintercept::Bool
    floorcoefs::Bool
    ndigits::Int # When printing
    params::Union{Nothing,Vector{Float64}}
    featurenames::Vector{String}
    MovingAverageLITRM(nfactors::Int, nfeatures::Int;
        minlag::Int = 3, maxlag::Int = 10,
        maxcoef = 1e3, mincoef = -maxcoef,
        hasintercept = true, floorcoefs = false, ndigits = 2,
        params = nothing,
        featurenames = String[]) = begin
        new(nfactors, nfeatures, minlag, maxlag, maxcoef, mincoef,
            hasintercept, floorcoefs, ndigits, params, featurenames)
    end
end

function parambounds(m::MovingAverageLITRM)
    bounds = Tuple{Float64, Float64}[]
    for _ in 1:nfactors(m)
        # Order of params: coef, feature, startlag, windowlen
        push!(bounds, (m.mincoef, m.maxcoef))
        # Since we will floor(Int, v) the following we add almost 1.0 to:
        push!(bounds, (1.0, float(m.nfeatures) + AlmostRoundUpTo1))
        push!(bounds, (float(m.minlag), float(m.maxlag) + AlmostRoundUpTo1))
        # Window length 1.0 means it is just a lagged factor with no
        # moving average window.
        push!(bounds, (1.0, float(m.maxlag) - float(m.minlag) + 1 + AlmostRoundUpTo1))
    end
    hasintercept(m) && push!(bounds, (m.mincoef, m.maxcoef))
    bounds
end

function isvalidparams(m::MovingAverageLITRM, ps::Vector{Float64})
    spans = windowspans(m, ps)
    (minlag(m) <= minimum(map(first, spans))) &&
    (maximum(map(last, spans)) <= maxlag(m)) &&
    all(t -> 1 <= t[2] <= (maxlag(m) - minlag(m) + 1),  spans)
end

windowspans(m::MovingAverageLITRM, ps::Vector{Float64} = params(m)) =
    Tuple{Int, Int, Int}[getfactor(m, fi, ps)[[3,4,5]] for fi in 1:nfactors(m)]

startofwindows(m::MovingAverageLITRM, ps::Vector{Float64} = params(m)) =
    map(first, windowspans(m, ps))

windowlengths(m::MovingAverageLITRM, ps::Vector{Float64} = params(m)) =
    map(t -> t[2], windowspans(m, ps))

endofwindows(m::MovingAverageLITRM, ps::Vector{Float64} = params(m)) =
    map(last, windowspans(m, ps))

# If we have params the smallest (since they are positive in the params vector)
# lag value is our horizon. If we don't have params we use the
actualminlag(m::MovingAverageLITRM, ps::Vector{Float64} = params(m)) =
    hasparams(m) ? minimum(startofwindows(m, ps)) : minlag(m)

actualmaxlag(m::MovingAverageLITRM, ps::Vector{Float64} = params(m)) =
        hasparams(m) ? maximum(endofwindows(m, ps)) : maxlag(m)

horizon(m::MovingAverageLITRM) = actualminlag(m)

function predictforrow(m::MovingAverageLITRM, X::AbstractMatrix{N}, row::Int,
    ps::Vector{Float64} = params(m)) where {N<:Number}
    prediction = hasintercept(m) ? ps[end] : 0.0
    for fi in 1:nfactors(m)
        coef, featureidx, winstart, winlen, wend = getfactor(m, fi, ps)
        #@show (row, winstart, winlen, wend, actualmaxlag(m), maxlag(m), ps)
        if winlen > 1
            windowavg = mean(X[(row-wend):(row-winstart), featureidx])
        else
            windowavg = X[row-winstart, featureidx]
        end
        prediction += (coef * windowavg)
    end
    return prediction
end

function windowfromparams(m::MovingAverageLITRM, ws::N1, wl::N2) where {N1<:Number,N2<:Number}
    wstart = floor(Int, ws)
    wend = min(maxlag(m), wstart + floor(Int, wl) - 1)
    wlen = wend - wstart + 1
    return wstart, wlen, wend
end

function getfactor(m::MovingAverageLITRM, fi::Int, ps::Vector{Float64} = params(m))
    i = 1 + (fi-1)*4
    c, f, ws, wl = ps[i:(i+3)]
    coef = m.floorcoefs ? floor(Int, c) : c
    return coef, floor(Int, f), windowfromparams(m, ws, wl)...
end

modelname(m::MovingAverageLITRM) =
    "MovingAverage factor LITRM (Linear Interpretable Timeseries Regression Model)"

function writefactor(io::IO, m::MovingAverageLITRM, fi::Int)
    c, f, ws, wlen, wend = getfactor(m, fi)
    coef = r(m, c)
    if coef != 0
        op = (coef > 0.0) ? "+" : "-"
        write(io, " $op $(abs(coef)) * ")
        f = "$(featurename(m, f))[T-$(ws)]"
        s = (wlen > 1) ? "MAvg($(f), $(wlen))" : f
        write(io, s)
    end
end


text:
const neghalflog2pi = -.5log(2pi)

#######################################
### loglik == Natural Parameter space
#######################################

"""
    loglik( g::GaussNP, x )

Exact log likelihood of a Gaussian expressed in the Natural Parameter space.
By default x is assumed to be of size D x N but if axis=1 then it will consider
it is given as N x D (where N=number of points, D=dimension).
"""
function loglik( g::GaussNP, x::AbstractArray{Float};
                 axis::Int=2 )::Union{Float,Vector{Float}}
    precmu   = g.theta1
    sqrtprec = chol(-g.theta2)
    tmp      = sqrtprec'\precmu
    fun(x)   = ( dot(x,g.theta2*x) + 2dot(x,precmu) )/2
    result   = sum(neghalflog2pi + log.(diag(sqrtprec))) - norm(tmp)^2/2
    result  += funarray(fun, x, length(precmu), axis)
end

"""
    loglik( g::DGaussNP, x )

Exact log likelihood of a Diagonal Gaussian expressed in the Natural Parameter
space. By default x is assumed to be of size D x N but if axis=1 then it will
consider it is given as N x D (where N=number of points, D=dimension).
"""
function loglik( g::DGaussNP, x::AbstractArray{Float};
                 axis::Int=2 )::Union{Float,Vector{Float}}
    precmu   = g.theta1
    sqrtprec = sqrt.(-g.theta2)
    tmp      = precmu./sqrtprec
    fun(x)   = ( dot(x, g.theta2 .* x) + 2dot(x, precmu) )/2
    result   = sum(neghalflog2pi + log.(sqrtprec)) - norm(tmp)^2/2
    result  += funarray(fun, x, length(precmu), axis)
end

#####################################
### loglik == Mean Parameter space
#####################################

"""
    loglik( g::GaussMP, x )

Exact log likelihood of a Gaussian expressed in the Mean Parameter space.
By default x is assumed to be of size D x N but if axis=1 then it will consider
it is given as N x D (where N=number of points, D=dimension).
"""
function loglik( g::GaussMP, x::AbstractArray{Float};
                 axis::Int=2 )::Union{Float,Vector{Float}}
    sqrtcov = chol(cov(g))
    fun(x)  = -norm(sqrtcov'\(x-g.mu1))^2/2
    result  = sum(neghalflog2pi - log.(diag(sqrtcov)))
    result += funarray(fun, x, length(g.mu1), axis)
end

"""
    loglik( g::DGaussMP, x )

Exact log likelihood of a Diagonal Gaussian expressed in the Mean Parameter
space. By default x is assumed to be of size D x N but if axis=1 then it will
consider it is given as N x D (where N=number of points, D=dimension).
"""
function loglik( g::DGaussMP, x::AbstractArray{Float};
                 axis::Int=2 )::Union{Float,Vector{Float}}
    stds    = std(g)
    fun(x)  = -norm((x-g.mu1)./stds)^2/2
    result  = sum(neghalflog2pi - log.(stds))
    result += funarray(fun, x, length(g.mu1), axis)
end

#######################################
### uloglik == Natural Parameter space
#######################################

"""
    uloglik( g::GaussNP, x )

Loglikelihood of a GaussNP without the constants (numerically more stable as
it does not require computing determinants and inversions). This is useful
in a sampler where the normalising constant can be neglected.
It corresponds to -0.5(<x, Px> -2<x,Pmu>).
By default x is assumed to be of size D x N but if axis=1 then it will consider
it is given as N x D (where N=number of points, D=dimension).
"""
function uloglik( g::GaussNP, x::AbstractArray{Float};
                  axis::Int=2 )::Union{Float,Vector{Float}}
    fun(x) = dot(x, g.theta1+g.theta2*x/2)
    funarray(fun, x, length(g), axis)
end

"""
    uloglik( g::DGaussNP, x )

Loglikelihood of a DGaussNP without the constants (numerically more stable as
it does not require computing determinants and inversions). This is useful
in a sampler where the normalising constant can be neglected.
It corresponds to -0.5(<x, Px> -2<x,Pmu>).
By default x is assumed to be of size D x N but if axis=1 then it will consider
it is given as N x D (where N=number of points, D=dimension).
"""
function uloglik( g::DGaussNP, x::AbstractArray{Float};
                  axis::Int=2 )::Union{Float,Vector{Float}}
    fun(x) = dot(x, g.theta1+g.theta2.*x/2)
    funarray(fun, x, length(g), axis)
end

#####################################
### uloglik == Mean Parameter space
#####################################

"""
    uloglik( g::GaussMP, x )

Loglikelihood of a GaussMP without the constants (numerically more stable as
it does not require computing determinants and inversions). This is useful
in a sampler where the normalising constant can be neglected.
It corresponds to -0.5(<x, Px> -2<x,Pmu>).
"""
function uloglik( g::GaussMP, x::AbstractArray{Float};
                  axis::Int=2 )::Union{Float,Vector{Float}}
    S = cov(g)
    fun(x) = dot(x, S\(g.mu1-x/2))
    funarray(fun, x, length(g), axis)
end

"""
    uloglik( g::DGaussMP, x )

Loglikelihood of a DGaussMP without the constants (numerically more stable as
it does not require computing determinants and inversions). This is useful
in a sampler where the normalising constant can be neglected.
It corresponds to -0.5(<x, Px> -2<x,Pmu>).
"""
function uloglik( g::DGaussMP, x::AbstractArray{Float};
                  axis::Int=2 )::Union{Float,Vector{Float}}
    S = cov(g)
    fun(x) = dot(x, (g.mu1-x/2)./S)
    funarray(fun, x, length(g), axis)
end

################
### gradloglik
################

# NOTE at the moment only handles point inputs

gradloglik(g::GaussNP,  x::Vector{Float}) = g.theta1 + g.theta2*x
gradloglik(g::DGaussNP, x::Vector{Float}) = g.theta1 + g.theta2 .* x
gradloglik(g::GaussMP,  x::Vector{Float}) = cov(g)\(g.mu1 - x)
gradloglik(g::DGaussMP, x::Vector{Float}) = (g.mu1 - x)./cov(g)


text:
module Shoco

using Shoco_jll

export compress, decompress

function compress(s::AbstractString)
    isempty(s) && return ""
    # The output should be no longer than the input
    compressed = Vector{UInt8}(undef, sizeof(s))
    # The function modifies `compressed` and returns the number of bytes written
    nbytes = ccall((:shoco_compress, libshoco), Csize_t,
                   (Ptr{Cchar}, Csize_t, Ptr{UInt8}, Csize_t),
                   s, 0, compressed, sizeof(s))
    nbytes > 0 || error("Compression failed for input $s")
    resize!(compressed, nbytes)
    return String(compressed)
end

function decompress(s::AbstractString)
    isempty(s) && return ""
    # The decompressed string will be at most twice as long as the input
    decompressed = Vector{UInt8}(undef, 2 * sizeof(s))
    nbytes = ccall((:shoco_decompress, libshoco), Csize_t,
                   (Ptr{Cchar}, Csize_t, Ptr{UInt8}, Csize_t),
                   s, sizeof(s), decompressed, 2 * sizeof(s))
    nbytes > 0 || error("Decompression failed for input $s")
    resize!(decompressed, nbytes)
    return String(decompressed)
end

end # module


text:
using StabTrussTopOpt
using Test

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


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

name = "Readline"
version = v"8.1.1"

# Collection of sources required to build Readline
sources = [
    ArchiveSource("https://ftp.gnu.org/gnu/readline/readline-$(version.major).$(version.minor).tar.gz",
                  "f8ceb4ee131e3232226a17f51b164afc46cd0b9e6cef344be87c65962cb82b02"),
    FileSource("https://ftp.gnu.org/gnu/readline/readline-$(version.major).$(version.minor)-patches/readline81-001",
               "682a465a68633650565c43d59f0b8cdf149c13a874682d3c20cb4af6709b9144"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/readline-*/

atomic_patch -p0 ${WORKSPACE}/srcdir/readline81-001
# Patch from https://aur.archlinux.org/cgit/aur.git/tree/fix_signal.diff?h=mingw-w64-readline
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-fix_signal.patch

export CPPFLAGS="-I${includedir}"
if [[ "${target}" == *-mingw* ]]; then
    # Only on Windows the library name embeds the ABI version
    export NCURSES_ABI_VER=6
    export LDFLAGS="-L${libdir}"
fi

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --with-curses
make -j${nproc} SHLIB_LIBS="-lncurses${NCURSES_ABI_VER}"
make install
install_license COPYING
"""

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

# The products that we will ensure are always built
products = Product[
    LibraryProduct(["libhistory", "libhistory8"], :libhistory),
    LibraryProduct(["libreadline", "libreadline8"], :libreadline),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Ncurses_jll"),
]

# Build the tarballs
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
using ArrayConstructors
using Compat.Test

@testset "constructors that take a scalar value and dims" begin
    @test Array(2.0, (2,2)) == fill(2.0, (2,2))
    @test Array{Float64}(2.0f0, (2,2)) == fill(2.0, (2,2))
    @test Array{Float64,2}(2.0f0, (2,2)) == fill(2.0, (2,2))
    @test Vector(2.0, (2,)) == fill(2.0, (2,))
    @test Matrix(2.0, (2,2)) == fill(2.0, (2,2))
end

@testset "constructors that take a non-scalar value and dims" begin
    A = rand(Float32, 2, 2)
    @test Array(A, (2,2)) == fill(A, (2,2))
    @test Array{Matrix{Float64}}(A, (2,2)) == fill(A, (2,2))
    @test Array{Matrix{Float64},2}(A, (2,2)) == fill(A, (2,2))
    @test Matrix(A, (2,2)) == fill(A, (2,2))

    # test for no aliasing
    B = Array(A, (2,2))
    @test all(x -> x !== A, B)
    @test all(x -> x  == A, B)
end

@testset "constructors that take an iterable with size" begin
    v = rand(Float32, 2)
    A = rand(Float32, 2, 2)
    @test Array(v) == v && Array(v) !== v
    @test Vector(v) == v && Vector(v) !== v
    @test Array(A) == A && Array(A) !== A
    @test Matrix(A) == A && Matrix(A) !== A
    @test Array{Float64}(v) == collect(Float64, v)
    @test Vector{Float64}(v) == collect(Float64, v)
    @test Array{Float64}(A) == collect(Float64, A)
    @test Matrix{Float64}(A) == collect(Float64, A)
    @test Array{Float64,1}(v) == collect(Float64, v)
    @test Array{Float64,2}(A) == collect(Float64, A)
end


text:
function load_result_data(file_name::String)
    file_name = get_res_full_path(file_name);
    if isfile(file_name)
        return FileIO.load(file_name)[data_key];
    end
    error("$(file_name) not found");
end

function load_result_description(file_name::String)
    file_name = get_desc_full_path(file_name);
    if isfile(file_name)
        return FileIO.load(file_name)[desc_key];
    end
    error("$(file_name) not found");
end

import FileIO;
import Rath2017PhD.Paths;
include("params.jl");
function rename!(old_name::String, new_name::String; overwrite = false)
    old_desc_name = get_desc_full_path(old_name);
    new_desc_name = get_desc_full_path(new_name);
    old_name = get_res_full_path(old_name)
    new_name = get_res_full_path(new_name);
    if isfile(old_name)
        !overwrite && isfile(new_name) && error("File $(new_name)"*
            " already exist!!! use overwrite keyword to obviously overwrite!!!")
        # moving res
        FileIO.mv(old_name, new_name, remove_destination = true);
        # moving desc
        if isfile(old_desc_name)
            FileIO.mv(old_desc_name, new_desc_name, remove_destination = true);
        else
            FileIO.save(new_desc_name, desc_key, get_default_description());
        end
    else
        error("$(old_name) not found");
    end
end

function save_result(file_res::String, data, description::String; force = false)
    res_path = get_res_full_path(file_res);
    desc_path = get_desc_full_path(file_res);
    !force && isfile(res_path) && error("File $(res_path)"*
        " already exist!!! use force keyword to overwrite!!!")
    FileIO.save(res_path, data_key, data);
    FileIO.save(desc_path, desc_key, description * get_default_description());
    println("Saved as $res_path");
end
save_result(file_name::String, data; description = get_default_description(), force = false) =
    save_result(file_name, data, description; force = force);

get_default_description() = " || Modified at $(now())";

function parse_basename(file_name, ext)
    !is_basename(file_name) && warn("$(file_name)"*
        " contain a parent dir!!!"*
        " We'll take only the basename $(basename(file_name))!!!")
    return add_ext(basename(file_name), ext);
end

exist_file(file_name) = isfile(get_res_full_path(file_name));

get_res_full_path(file_name) = joinpath(Paths.results_path, parse_basename(file_name, res_file_ext));
get_desc_full_path(file_name) = joinpath(Paths.descriptions_path, parse_basename(file_name, desc_file_ext));
list_results() = get_ext_files(Paths.results_path, res_file_ext);
list_desc_files() = get_ext_files(Paths.descriptions_path, res_file_ext);

function add_ext(file_name, ext_to_add)
    ext = extension(file_name);
    if ext != ""
        return replace(file_name, ext => ext_to_add);
    else
        return file_name *= ext_to_add;
    end
end

extension(file_name::String) = try match(r"\.[A-Za-z0-9]+$", file_name).match catch "" end

is_basename(basename) = dirname(basename) == "";

function status()
    res_names = basename.(list_results());
    for (i, name) in enumerate(res_names)
        desc = missing_desc_file_err_msg;
        try desc = load_result_description(name); catch end
        println("----------------------------------------------------");
        println(i, ": ", name)
        println("--------------------------------------");
        println();
        println(desc);
        println();
    end
end

get_ext_files(dir_path, ext) = filter((f) -> endswith(f, ext), readdir(dir_path));

function remove!(file_name)
    res_path = get_res_full_path(file_name);
    desc_path = get_desc_full_path(file_name);
    rm(res_path; force = true);
    rm(desc_path; force = true);
    println("$res_path was deleted!!!")
end

function edit_description!(file_name, description)
    res_path = get_res_full_path(file_name);
    desc_path = get_desc_full_path(file_name);
    if exist_file(file_name)
        FileIO.save(desc_path, desc_key, description * get_default_description());
    else
        error("$file_name not found!!!")
    end
end


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# Gaston test framework
#
# Note: once julia has its own test framework, we might want to migrate to it.
#
# There are three kinds of tests:
#   - Tests that should error -- things that Gaston should not allow, such as
#     building a figure with an invalid plotstyle. In these cases, Gaston
#     should issue an error().
#     For these tests, we use the test_error macro.
#   - Tests that should succeed -- things that should work without producing
#     an error.
#     For these tests, we use the test_success macro.
#   - Test that certain functions return certain values. These tests help to
#     verify the functions are working correctly, not just erroring or
#     succeeding.
#     For these tests, we use the test_success macro in combination with
#     @assert.
#
# In all cases, when a test fails, the expression under test is echoed to
# the screen to help identify it.
#
# Note: the test macros are specific to Gaston and are not meant to be used
# for testing other programs. Specifically, since the macros are not
# hygienic (as far as I understand), the expression passed to them could
# mess with the test statistics.

macro test_error(ex)
    quote
        global testnumber
        global testsrun
        global testspassed
        run(`sleep 0.3`)
        testnumber = testnumber + 1
        testsrun = testsrun + 1
        s = string("Test number ", testnumber, ". Error expected. Result: ")
        try
            $ex
            println(string(s, "Success (Test failed.)"))
            println($(string(ex)))
        catch
            println(string(s, "Error (Test passed.)"))
            testspassed = testspassed + 1
        end
    end
end

macro test_success(ex)
    quote
        global testnumber
        global testsrun
        global testspassed
        run(`sleep 0.4`)
        testnumber = testnumber + 1
        testsrun = testsrun + 1
        s = string("Test number ", testnumber, ". Success expected. Result: ")
        try
            $ex
            println(string(s, "Success (Test passed.)"))
            testspassed = testspassed + 1
        catch
            println(string(s, "Error (Test failed.)"))
            println($(string(ex)))
        end
    end
end

function run_tests_error(ini)
    testnumber = ini
    testspassed = 0
    testsrun = 0

    # high-level functions: figure
    @test_error figure(-1)
    @test_error figure("invalid")
    @test_error figure(1.0)
    @test_error figure(1:2)
    # high-level functions: closefigure
    @test_error closefigure(-1)
    @test_error closefigure("invalid")
    @test_error closefigure(1.0)
    @test_error closefigure(1:2)
    # high-level functions: plot and histogram
    for op = (:plot, :histogram)
        @test_error op("linewidth")
        @test_error op(0:10,0:11)
        @test_error op(0:10+im*0:10)
        @test_error op(0:10,"legend",0)
        @test_error op(0:10,"plotstyle","invalid")
        @test_error op(0:10,"marker","invalid")
        @test_error op(0:10,"marker",0)
        @test_error op(0:10,"linewidth","b")
        @test_error op(0:10,"linewidth",im)
        @test_error op(0:10,"pointsize","b")
        @test_error op(0:10,"pointsize",im)
        @test_error op(0:10,"title",0)
        @test_error op(0:10,"xlabel",0)
        @test_error op(0:10,"ylabel",0)
        @test_error op(0:10,"zlabel","z")
        @test_error op(0:10,"axis","invalid")
    end
    @test_error histogram(0:10,"plotstyle","lines")
    @test_error histogram(0:10,"bins")
    @test_error histogram(0:10,"bins","invalid")
    @test_error histogram(0:10,"bins",-1)
    @test_error histogram(0:10,"bins",3.1)
    @test_error histogram(0:10,"norm","invalid")
    @test_error histogram(0:10,"norm",-1)
    @test_error histogram(0:10,"marker","ecircle")
    @test_error histogram(0:10,"pointsize",1)
    @test_error histogram(0:10,"axis","loglog")
    # high-level functions: imagesc
    z = rand(3,3)
    @test_error imagesc(1:5,z)
    @test_error imagesc("title","none")
    @test_error imagesc(z,"clim",3)
    @test_error imagesc(1:3)
    @test_error imagesc(z,"none","none")
    # high-level functions: print
    @test_error begin
        closeall()
        printfigure(1)
    end
    @test_error begin
        closeall()
        plot(1,sin(-3:0.1:3))
        printfigure(1,"none")
    end
    @test_error begin
        closeall()
        plot(1,sin(-3:0.1:3))
        printfigure(2)
    end
    # high-level functions: surf
    @test_error surf("property")
    @test_error surf("title","none")
    @test_error surf(1:10)
    @test_error surf(1:10,1:10)
    @test_error begin
        z = rand(5,5)
        surf(1:5,1:4,z)
    end
    @test_error surf(1:10,(x,y)->x*y)
    # addcoords
    @test_error addcoords(im*(1:10))
    @test_error addcoords(["a" "b"])
    @test_error addcoords(1:10,1:10,1:10)
    @test_error addcoords(1:10,1:11)
    @test_error addcoords(1:2,1:3,[[1 2 3 4],[1 2 3 4]])
    # set_terminal
    @test_error set_terminal("none")
    # @osx_only @test_error set_terminal("aqua") # no support on recent OSX
    ## tests that should fail, but (still) don't
    ## commented out because gnuplot barfs all over the screen
    #@test_error plot(0:10,"color","nonexistant")
    #@test_error plot(0:10,"box","invalid")

    return testsrun, testspassed
end

function run_tests_success(ini)
    testnumber = ini
    testspassed = 0
    testsrun = 0

    closeall()
    # high-level functions: closefigure
    @test_success closefigure()
    @test_success @assert 0 == closefigure()
    @test_success @assert 0 == closefigure(10)
    # high-level functions: figure
    @test_success @assert 1 == figure()
    @test_success @assert 1 == figure(1)
    @test_success @assert 3 == figure(3)
    @test_success begin
        closeall()
        @assert 1 == figure(1)
        @assert 2 == figure()
        @assert 5 == figure(5)
        closefigure(2)
        @assert 2 == figure()
        closeall()
    end
    closeall()
    # high-level functions: plot
    @test_success begin
        plot(0:10)
        closeall()
    end
    @test_success begin
        plot(sin(-3:0.01:3),
            "legend", "sine",
            "plotstyle", "lines",
            "color","blue",
            "marker","ecircle",
            "linewidth",2,
            "pointsize",1.1,
            "title","test plot 1",
            "xlabel","x",
            "ylabel","y",
            "box","inside horizontal left top",
            "axis","loglog")
        closeall()
    end
    # high-level functions: histogram
    @test_success begin
        histogram(rand(1000))
        closeall()
    end
    @test_success begin
        histogram(randn(1000),
            "bins",100,
            "norm",1,
            "color","blue",
            "linewidth",2,
            "title","test histogram",
            "xlabel","x",
            "ylabel","y",
            "box","inside horizontal left top")
        closeall()
    end
    # high-level functions: imagesc
    z = rand(5,6)
    @test_success begin
        imagesc(z,"title","test imagesc 1","xlabel","xx","ylabel","yy")
        closeall()
    end
    @test_success begin
        imagesc(1:5,z,"title","test imagesc 2","xlabel","xx","ylabel","yy")
        closeall()
    end
    @test_success begin
        imagesc(1:5,1:6,z,"title","test imagesc 3","xlabel","xx","ylabel","yy")
        closeall()
    end
    # high-level functions: surf
    @test_success surf(rand(10,10))
    @test_success surf(2:11,rand(10,10))
    @test_success surf(0:9,2:11,rand(10,10))
    @test_success surf(0:9,2:11,(x,y)->x*y)
    @test_success surf(0:9,2:11,(x,y)->x*y,"title","test","plotstyle","pm3d")
    # high-level functions: print
    @test_success begin
        closeall()
        plot(1,sin(-3:0.1:3))
        sleep(0.1)
        set_filename("/dev/null")
        printfigure()
        closeall()
    end
    @test_success begin
        closeall()
        plot(2,sin(-3:0.1:3))
        sleep(0.1)
        set_filename("/dev/null")
        printfigure(2,"png")
        closeall()
    end
    @test_success begin
        closeall()
        plot(3,sin(-3:0.1:3))
        sleep(0.1)
        set_filename("/dev/null")
        printfigure(3,"eps")
        closeall()
    end
    @test_success begin
        closeall()
        plot(4,sin(-3:0.1:3))
        sleep(0.1)
        set_filename("/dev/null")
        printfigure(4,"pdf")
        closeall()
    end
    @test_success begin
        closeall()
        plot(5,sin(-3:0.1:3))
        sleep(0.1)
        set_filename("/dev/null")
        set_print_size("640,480")
        printfigure("svg")
        closeall()
    end
    @test_success begin
        closeall()
        plot(6,sin(-3:0.1:3))
        sleep(0.1)
        set_filename("/dev/null")
        set_print_size("640,480")
        printfigure("gif")
        closeall()
    end
    # type instantiation
    @test_success CurveConf()
    # setting terminal
    @test_success set_terminal("x11")
    @test_success set_terminal("wxt")
    #@osx_only @test_success set_terminal("aqua") # no support on recent OSX

    return testsrun, testspassed
end

function gaston_tests()
    println("Running tests...")
    (total,passed) = run_tests_error(0)
    (total1,passed1) = run_tests_success(total)
    s = println(string("Tests run: ", string(total+total1)))
    s = println(string("Tests passed: ", string(passed+passed1)))
end


text:
module Avatar

export generate

using HTTP
using Images
using JSON

""""""
function fetch(palette_id)
    response = HTTP.request("GET", "http://www.colourlovers.com/api/palette/$palette_id?format=json")
    json = JSON.parse(String(response.body))
    palette = json[1]["colors"]
    return map(c -> parse(Colorant, "#" * c), palette)
end

""""""
function generate(cells_x=5, cells_y=5, cell_size_x=80, cell_size_y=80,
                    palette_id=379413, out="avatar.png")
    palette = fetch(palette_id)
    img = rand(palette, cells_x, cells_y)
    img = repeat(img, inner=(cell_size_x, cell_size_y))
    save(out, img)
end

end # module Avatar


text:
module Caesar

# import RoME: initfg # collision on RoME.initfg() since no parameters are given in both RoME and Caesar
import Distributions: Normal
import RoME: getRangeKDEMax2D, getLastPose
import IncrementalInference: batchSolve!, getSample, initfg

using Reexport

@reexport using RoME
@reexport using IncrementalInference
@reexport using KernelDensityEstimate
@reexport using Distributions

using
  DelimitedFiles,
  Distributed,
  Statistics,
  LinearAlgebra,
  IncrementalInference,
  Graphs,
  TransformUtils,
  CoordinateTransformations,
  Rotations,
  JSON,
  JSON2,
  FileIO,
  DataStructures,
  ProgressMeter,
  ImageMagick,
  ImageCore,
  DocStringExtensions,
  CloudGraphs, # TODO: will be movedd to DFG
  Neo4j, # TODO: will be movedd to DFG
  Mongoc, # TODO: will be movedd to DFG
  Unmarshal,
  YAML,
  FFTW

export
  GenericInSituSystem,  # insitu components
  makeGenericInSituSys,
  InSituSystem,
  makeInSituSys,
  triggerPose,
  poseTrigAndAdd!,
  processTreeTrackersUpdates!,
  advOdoByRules,
  progressExamplePlot,
  plotPoseDict,
  plotTrckStep,
  SLAMWrapper,

  # servers
  tcpStringSLAMServer,
  tcpStringBRTrackingServer,

  # save and load data
  saveSlam, # TODO deprecate
  loadSlam, # TODO deprecate
  haselement,

  # user functions
  identitypose6fg,
  projectrbe,
  hasval,

  # repeats from RoME and IIF
  # initfg,
  # addNode!,
  # addFactor!,

  # CloudGraphs helper functions
  insertnodefromcv!,
  checkandinsertedges!,
  getbinarraymongo,
  gettopoint,
  getdotwothree,
  bin2arr,
  fetchsubgraph!,
  getVertNeoIDs!,
  insertrobotdatafirstpose!,
  tryunpackalltypes!,
  fetchrobotdatafirstpose,
  getExVertexNeoIDs,
  db2jld,

  # Robot Utils
  getRangeKDEMax2D,

  # would be CloudGraphs calls
  hasBigDataElement,
  getBigDataElement,
  removeNeo4jID,

  # solver service SLAMinDB
  getcredentials,
  startSlamInDb,
  runSlamInDbOnSession,
  slamindb,
  convertdb,
  resetconvertdb,
  getmaxfactorid,

  # webserver
  SolverStatus,
  CaesarConfig,
  IterationStatistics,
  VisualizationConfig,

  # multisession utils
  multisessionquery,
  parsemultisessionqueryresult!,
  getLandmOtherSessNeoIDs,
  getAllLandmarkNeoIDs,
  getLocalSubGraphMultisession,
  findExistingMSConstraints,
  getprpt2kde,
  rmInstMultisessionPriors!,
  removeMultisessions!,

  # sas-slam
  CBFFilterConfig,
  CZTFilter,
  prepCZTFilter,
  getCBFFilter2Dsize,
  constructCBFFilter2D!,
  CBF2D_DelaySum!,
  MatchedFilter,
  SASBearing2D,
  PackedSASBearing2D,
  compare,
  SASDebug,
  reset!,
  prepMF,
  loadConfigFile,
  prepareRangeModel,
  prepareSAS2DFactor,
  wrapRad,
  phaseShiftSingle!,
  liebf!,
  SASDebug

NothingUnion{T} = Union{Nothing, T}

include("BearingRangeTrackingServer.jl")

include("SlamServer.jl")
include("DataUtils.jl")
include("UserFunctions.jl")

# Configuration
include("config/CaesarConfig.jl")

# using CloudGraphs
include("cloudgraphs/SolverStatus.jl")
include("cloudgraphs/IterationStatistics.jl")
include("cloudgraphs/CloudGraphIntegration.jl") # Work in progress code
include("cloudgraphs/ConvertGeneralSlaminDB.jl")
include("cloudgraphs/slamindb.jl")
include("cloudgraphs/MultisessionUtils.jl")
include("cloudgraphs/FoveationUtils.jl")

# ZMQ server and endpoints
include("zmq/ZmqCaesar.jl")

# Multisession operation
include("multisession/Multisession.jl")

# SAS-SLAM
include("beamforming/czt.jl")
include("beamforming/CBF.jl")
include("beamforming/MatchedFilter.jl")
include("beamforming/SASBearing2D.jl")
include("beamforming/SASUtils.jl")

end


text:
module FakePackage

include("subfolder1/file1.jl")
include("subfolder2/file2.jl")

end


text:
existingfile(path) = isfile(path) ? path : nothing

_tomlpath(dir, candidates) =
    something(existingfile.(joinpath.(dir, candidates))...,
              joinpath(dir, candidates[2]))

projecttomlpath(dir) = _tomlpath(dir, ("JuliaProject.toml", "Project.toml"))
manifesttomlpath(dir) = _tomlpath(dir, ("JuliaManifest.toml", "Manifest.toml"))

function pkgat(path::AbstractString) :: PkgId
    prj = TOML.parsefile(projecttomlpath(path))
    return PkgId(UUID(prj["uuid"]), prj["name"])
end

function pkgspecof(path::AbstractString; prefer_https=false)
    prj = TOML.parsefile(projecttomlpath(path))
    if prefer_https
        url = vcslinktoroot(path=path)
    else
        url = strip(read(git_cmd(`config remote.origin.url`, path), String))
    end
    tree_sha = strip(read(git_cmd(`rev-parse "HEAD^{tree}"`, path), String))
    spec = Pkg.PackageSpec(
        name = prj["name"],
        uuid = prj["uuid"],
        url = url,
    )
    @set! tree_hash(spec) = Base.SHA1(tree_sha)
    return spec
end

tree_hash(spec) =
    isdefined(spec, :tree_hash) ? spec.tree_hash : spec.repo.tree_sha

Setfield.set(spec, ::typeof(@lens tree_hash(_)), value) =
    if isdefined(spec, :tree_hash)
        # Julia 1.2
        @set spec.tree_hash = value
    else
        # Julia 1.1
        @set spec.repo.tree_sha = value
    end

function tempbak(name::AbstractString)
    newname, io = mktemp(dirname(name))
    close(io)
    rename(name, newname)
    return newname
end

ensuredir(path) = isdir(path) ? path : dirname(path)

git_cmd(args::Cmd=``, dir::AbstractString=".") =
    setenv(`git --no-pager $args`; dir=ensuredir(dir))

git_is_clean(path::AbstractString) =
    isempty(read(git_cmd(`status --short --untracked-files=no`, path)))

vcslinktoroot(args...; kwargs...) =
    pyimport("vcslinks").root(args...; kwargs...)

vcslinktocommit(args...; kwargs...) =
    pyimport("vcslinks").commit(args...; kwargs...)

function commitmessage(fullrev, uppkgid::PkgId, from, title, comment)
    link = vcslinktocommit(fullrev, path=from)
    subject = strip(
        read(
            git_cmd(`show --format=format:%s --no-patch $fullrev`, from),
            String,
        ),
    )
    title = rstrip(something(title, "Update: $(uppkgid.name)"))
    comment = rstrip(something(comment, ""))
    footer = """
    Using commit:
    $subject
    $link
    """
    return join([title, comment, footer], "\n\n")
end

function _is_api(f)
    for m in methods(f)
        if m.file === Symbol(joinpath(@__DIR__, "api.jl"))
            return true
        end
    end
    return false
end


text:
#brexTools.jl contains functions for input/output (saving a reading results), and various other tools.

#Function gen_filename generates the name of folder where results are saved. This is based on time of script execution, and the name of the parameter file.
function genFolderName(parfile)
    #If parfile ends with ".jl", do not use that part in generating folder name
    #(this does not handle other file types, but parameter files should be .jl and
    #it is not a problem in Unix if the folder name contains a dot)
    if parfile[end-2:end] == ".jl"
        parfileloc = parfile[1:end-3]
    else
        parfileloc = parfile
    end

    #Get year, month, day, hour, minute, second as strings, prepend zeros as necessary
    #(for correct sorting when some numbers are less than 10)
    dat = Dates.now()
    y = prep0(year(dat),4)
    m = prep0(month(dat),2)
    d = prep0(day(dat),2)
    h = prep0(hour(dat),2)
    mi = prep0(minute(dat),2)
    s = prep0(second(dat),2)

    #Construct foldername and return it.
    #format is yymmdd_hhmmss_parfile
    foldername = y*m*d*"_"*h*mi*s*"_"*parfileloc
    return foldername
end

#Function prep0 converts number num to a string, and adds leading zeroes such that the length of the string is at least minlength
function prep0(num,minlength)
    num_str = string(num)
    #prepend zeroes
    for i=1:(minlength-length(num_str))
       num_str = "0"*num_str
    end
    return num_str
end

#Function saveAll saves stationary equilibria and transition paths in JLD (binary) format
function saveAll(foldername,SE,TP)
    #For saving complex data types use JLD format (need using JLD)
    #The file needs to have suffix .jld so the type is correctly inferred.
    save("$foldername/SE.jld","SE",SE)
    save("$foldername/TP.jld","TP",TP)

    #Below: saving data in text format (to be used for simulation if desired)
    #saveInFile(V,"test.out",foldername,text = false,replace = true)
    #loadFromFile!(V,"test.out",foldername,text = false)
end

#Functions saveInFile can be used to manually save data in binary or text format, not using JLD format. This can be useful if we want to visually inspect the output of the program (such as simulation results), or import it in some other program for further analysis.

#Function saveInFile saves variable x in folder foldername
#By default data is saved in binary form (this is usually faster).
#If text = true, the data is saved in text form so we can read it easily
function saveInFile(x,filename,foldername;text = false,replace = true)
    fullpath = "$foldername/$filename"

    #If the file exists and replace = true (default), delete it and continue
    #If the file exists and replace = false, write a warning message and exit the function
    if isfile(fullpath)
       if replace
            rm(fullpath)
        else
           println("Warning: File $fullpath already exists but replace = false, so is not replaced.")
            return nothing
        end
    end
    io = open(fullpath,"w")
    if text #text output
            show(io,x)
        else #binary output
            write(io,x)
    end
    close(io)
end

#Function loadFromFile loads a given variable from file
#text has to be the same as when the variable was saved
function loadFromFile!(x,filename,foldername;text = false)
    fullpath = "$foldername/$filename"

    #If text = true then the original data was saved in text form using function show,
    #Implementation of reading text data not done yet
    #A way to do this is to use package DelimitedFiles to save data and then load them
    #(low priority change)
    if text
       error("text = true not yet implemented in function loadFromFile.")
    end
    io = open(fullpath,"r")
            read!(io,x)
    close(io)
end


text:
"""
Flow.Format consists uniquely of optional methods for making it easier to
visualize the structure of the Flow model. It requires the user to install and
import packages not in the GuillotineModels dependencies in order to have
access to its methods. The `flow2digraph` method needs LightGraphs, and the
other three (`flow2file`, `flow2tikzpic`, `flow2pdf`) need LightGraphs,
TikzGraphs, and TikzGraphs.

The methods signatures are:

    flow2digraph(num_nodes, edges) :: LightGraphs.SimpleDiGraph
    flow2tikzpic(nodes, edges) :: TikzPictures.TikzPicture
    function flow2file(nodes, edges, out :: TikzPictures.SaveType) :: Nothing
    function flow2pdf(nodes, edges, fname) :: Nothing

The `Node` and `Edge` types are defined by Flow.Enumeration and the Vectors
of all nodes and all edges of the model are returned by
`GuillotineModels.build_model(::Val{:Flow}, ...)`.
"""
module Format

# CONSIDER MIGRATING TO GraphRecipes
# (https://github.com/JuliaPlots/GraphRecipes.jl/blob/master/README.md)
# pros: seems more mature, use the Plots library
# cons: seems even less documented

using ..Enumeration

import Requires.@require

function __init__()
	@require LightGraphs="093fc24a-ae57-5d10-9952-331d41423f4d" begin
		export flow2digraph
		function flow2digraph(
			num_nodes :: N,
			edges :: Vector{Edge{N, E}}
		) :: LightGraphs.SimpleDiGraph{E} where {N, E}
			g = LightGraphs.SimpleDiGraph{E}(num_nodes)

			for e in edges
				println(e.indx, ": ", e.head, " -> ", e.tail)
				LightGraphs.add_edge!(g, e.head, e.tail)
			end

			g
		end

		@require TikzPictures="37f6aa50-8035-52d0-81c2-5a1d08754b2d" begin
		@require TikzGraphs="b4f28e30-c73f-5eaf-a395-8a9db949a742" begin
			export flow2tikzpic, flow2file, flow2pdf
			function flow2tikzpic(
				nodes :: Vector{Node{S, N}},
				edges :: Vector{Edge{N, E}}
			) :: TikzPictures.TikzPicture where {S, N, E}
				@assert issorted(nodes, by = n -> n.idx)
				d = Dict{Tuple{E, E}, String}()
				for e in edges
					arc = convert.(E, (e.head, e.tail))
					if !haskey(d, arc)
						d[arc] = "$(e.indx)+$(e.back)"
					else
						d[arc] *= ", \\\\ $(e.indx)+$(e.back)"
					end
				end

				TikzGraphs.plot(
					flow2digraph(length(nodes), edges),
					labels = map(n -> "$(n.idx):$(n.par):$(n.per):$(n.ori)", nodes),
					TikzGraphs.Layouts.SimpleNecklace(),
					options="scale=3, every node/.style={align=center}",
					edge_labels = d
				)
			end

			function flow2file(
				nodes :: Vector{Node{S, N}},
				edges :: Vector{Edge{N, E}},
				out :: TikzPictures.SaveType
			) :: Nothing where {S, N, E}
				TikzPictures.save(out, flow2tikzpic(nodes, edges))
				nothing
			end

			function flow2pdf(
				nodes :: Vector{Node{S, N}},
				edges :: Vector{Edge{N, E}},
				fname :: String
			) :: Nothing where {S, N, E}
				flow2file(nodes, edges, TikzPictures.PDF(fname))
				nothing
			end
		end # TikzPictures
		end # TikzGraphs
	end # LightGraphs
end # __init__
end # module