text:
function matrix_mul(y, x, iflag::Ti; A::Matrix)

    # forward
    if iflag == 1
       mul!(y, A, x)

    # adjoint
    elseif iflag == 2
       mul!(y, transpose(A), x)

    else
       error("undefined operation")
    end

    return y
end

# set up the problem
for T in (Float32, Float64, ComplexF32, ComplexF64)
end

m = 100; n = 100;
A  = rand(T, m, n)
x0 = rand(T, n)
d  = A * x0

(x, his) = cgls(matrix_mul!, y; op_params=params, maxIter=m)

"""
   modified for memory saving
"""
function cgls(A::Tf, b::Vector{Tv}; op_params::NamedTuple=NamedTuple(),
              x0=[], shift=0.0, tol=1e-6, maxIter=50, print_flag=true) where {Tf<:Union{Function, Matrix}, Tv<:Number}

    # determine the type of A(matrix or operator)
    if isa(A, Matrix)
       explicitA = true

    elseif isa(A, Function)
       explicitA = false

    else
       error("A can be either a matrix or a function")
    end

    # initial guess of the unknowns
    if explicitA

       (m, n) = size(A)

       # provided initial guess
       if x0 != 0.0
          x = copy(x0)

          r = b - A * x
          s = A' * r - shift * x

       # zero vector as initial guess
       else
          x = zeros(eltype(b), n)
          r = copy(b)
          s = A' * r
       end

    else

       m = length(b)

       # provided initial guess
       if x0 != 0.0
          x = copy(x0)
          r = b - A(x, 1; op_params...)
          s = A(r, 2; op_params...) - shift * x
          n = length(s)
       # zero vector as initial guess
       else
          r = copy(b)
          s = A(b, 2; params...)
          n = length(s)
          x = zeros(eltype(b), n)
       end

    end

    data_fitting = dot(r,r)
    constraint   = 0.0
    cost0 = data_fitting + constraint

    # initialize some intermidiate vectors
    p = copy(s)
    norms0= vecnorm(s)
    gamma = norms0^2
    normx = vecnorm(x)
    xmax  = normx     # keep the initial one
    resNE = 1.0

    gamma0= copy(gamma)
    delta = 0.0

    # iteration counter and stop condition
    k = 0
    run_flag = true

    if print_flag
       header = "  k         data_fitting           constraint             normx                resNE"
       println(""); println(header);
       @printf("%3.0f %20.10e %20.10e %20.10e %20.10e\n", k, data_fitting, constraint, normx, resNE);
    end

    while (k < maxIter) && run_flag

          k = k + 1

          if explicitA
             q = A * p
          else
             q = A(p, 1; params...)
          end

          delta = (vecnorm(q))^2 + shift * (vecnorm(p))^2
          indefinite = delta <= 0 ? true : false
          delta      = delta == 0.? eps(): delta

          alpha = gamma / delta

          x = x + alpha * p
          r = r - alpha * q

          data_fitting = dot(r,r)
          constraint   = shift * dot(x, x)
          cost = data_fitting + constraint

          if explicitA
             s = A' * r - shift * x
          else
             s = A(r, 2; params...) - shift * x
          end

          norms  = norm(s)
          gamma0 = copy(gamma)
          gamma  = norms^2
          beta   = gamma / gamma0

          p = s + beta * p

          # check the stopping crietia
          normx = norm(x)
          xmax  = normx > xmax ? normx : xmax
          if norms <= norms0 * tol || normx * tol >= 1.0
             run_flag = false
          end

          # print information
          # resNE = norms / norms0
          resNE = cost  / cost0
          if print_flag
             @printf("%3.0f %20.10e %20.10e %20.10e %20.10e\n", k, data_fitting, constraint, normx, resNE);
          end

    end

    return x, resNE

end

"""
   CGLS algorithm for linear operators whose input and output can be save in memory
"""
function cgls(A::Tf, b::Vector{Tv}; op_params::NamedTuple=NamedTuple(),
              x0=[], shift=0.0, tol=1e-6, maxIter=50, print_flag=true) where {Tf<:Union{Function, Matrix}, Tv<:Number}

    # determine the type of A(matrix or operator)
    if isa(A, Matrix)
       explicitA = true

    elseif isa(A, Function)
       explicitA = false

    else
       error("A can be either a matrix or a function")
    end

    # initial guess of the unknowns
    if explicitA

       (m, n) = size(A)

       # provided initial guess
       if x0 != 0.0
          x = copy(x0)
          r = b - A * x
          s = A' * r - shift * x
       # zero vector as initial guess
       else
          x = zeros(eltype(b), n)
          r = copy(b)
          s = A' * r
       end

    else

       m = length(b)

       # provided initial guess
       if x0 != 0.0
          x = copy(x0)
          r = b - A(x, 1; op_params...)
          s = A(r, 2; op_params...) - shift * x
          n = length(s)
       # zero vector as initial guess
       else
          r = copy(b)
          s = A(b, 2; params...)
          n = length(s)
          x = zeros(eltype(b), n)
       end

    end

    data_fitting = dot(r,r)
    constraint   = 0.0
    cost0 = data_fitting + constraint

    # initialize some intermidiate vectors
    p = copy(s)
    norms0= vecnorm(s)
    gamma = norms0^2
    normx = vecnorm(x)
    xmax  = normx     # keep the initial one
    resNE = 1.0

    gamma0= copy(gamma)
    delta = 0.0

    # iteration counter and stop condition
    k = 0
    run_flag = true

    if print_flag
       header = "  k         data_fitting           constraint             normx                resNE"
       println(""); println(header);
       @printf("%3.0f %20.10e %20.10e %20.10e %20.10e\n", k, data_fitting, constraint, normx, resNE);
    end

    while (k < maxIter) && run_flag

          k = k + 1

          if explicitA
             q = A * p
          else
             q = A(p, 1; params...)
          end

          delta = (vecnorm(q))^2 + shift * (vecnorm(p))^2
          indefinite = delta <= 0 ? true : false
          delta      = delta == 0.? eps(): delta

          alpha = gamma / delta

          x = x + alpha * p
          r = r - alpha * q

          data_fitting = dot(r,r)
          constraint   = shift * dot(x, x)
          cost = data_fitting + constraint

          if explicitA
             s = A' * r - shift * x
          else
             s = A(r, 2; params...) - shift * x
          end

          norms  = norm(s)
          gamma0 = copy(gamma)
          gamma  = norms^2
          beta   = gamma / gamma0

          p = s + beta * p

          # check the stopping crietia
          normx = norm(x)
          xmax  = normx > xmax ? normx : xmax
          if norms <= norms0 * tol || normx * tol >= 1.0
             run_flag = false
          end

          # print information
          # resNE = norms / norms0
          resNE = cost  / cost0
          if print_flag
             @printf("%3.0f %20.10e %20.10e %20.10e %20.10e\n", k, data_fitting, constraint, normx, resNE);
          end

    end

    return x, resNE

end


text:
using PromiseExtractor
using Base.Test

# write your own tests here
@test 1 == 1

@test promise_for("""

    <link rel="stylesheet" type="text/css" href="example.css">

""", "http://localhost") == ["example.css","http://localhost/example.jpg","example.js"]

ext: jl


text:
include("board.jl")
using BM
gameType = "standard"
board = startGame(gameType)
global states = Array(Dict{Tuple{Int64,Int64},Tuple{Char,Char}},0)
push!(states,board.state)
seed = time()
calculationTime = 100 #Max calculation time permitted in seconds
max_moves = 500
#global wins = Dict('w'=>Dict(),'b'=>Dict())
#global plays = Dict('w'=>Dict(),'b'=>Dict())
global wins = Dict()
global plays = Dict()
global C = 1.4

function update(state)
  push!(states,state)
end

function get_play()
  global max_depth,states,wins,plays
  global max_depth = 0
  state = states[size(states)[1]]
  player = currentPlayer(states)
  legal = legalMovesPlayer(state,player)

  if size(legal)[1] == 0
    return
  elseif size(legal)[1] == 1
    return legal[1]
  end

  games = 0
  bTime = time()
  while time()-bTime < calculationTime
  #  println(time()-bTime)
    run_simulation()
    games += 1

    print(STDOUT, "\u1b[1G")   # go to first column
    print(STDOUT,time()-bTime)
    print(STDOUT, "\u1b[K")    # clear the rest of the line
  end

  moves_states = collect( (p,next_state(state,p)) for p in legal)

  println("$games, $(time()-bTime)")

  (percent_wins,moveIndex) = findmax(
  get(wins,(player,S),0)/get(plays,(player,S),1)
  for (p,S) in moves_states )
  move = moves_states[moveIndex][1]

  #=Print stats?=#
  println("Maximum depth searched: $(max_depth)")
  println("with %wins $percent_wins : move:")
  return move
end

function run_simulation() #2,3,5
  global max_depth,states,plays,wins
  visited_states = Set()
  states_copy = deepcopy(states)
  state = states_copy[size(states_copy)[1]]
  player = currentPlayer(states_copy)
  #println(player)
  expandP = true
  for i in 2:(max_moves+1)
    state = states_copy[size(states_copy)[1]]
    player = currentPlayer(states_copy)
    legalC  = legalMovesPlayer(state,player)
    moves_states = collect( (p,next_state(state,p)) for p in legalC)

    if all(collect(get(plays,(p,S),0)!=0 for (p,S) in moves_states))#Have stats
      log_total = log(
        sum(plays[(player,S)] for (p,S) in moves_states)
      )
    #  println("have all stats")
      (value,moveIndex)=findmax(
      collect(wins[(player,S)]/plays[(player,S)] + C*sqrt(log_total / plays[(player,S)]) for (p,S) in moves_states))

      (move,state) = moves_states[moveIndex]
    else #Random choice
      move = legalC[mod(Int(ceil(seed))*i,size(legalC)[1])+1] #Choose a "random" move for CURRENT PLAYER from legalC seeded from 'seed'
      state = next_state(state, move)
    end

    push!(states_copy,state)

    if expandP && !( (player,state) in keys(plays) )
      expandP = false
      plays[(player,state)] = 0
      wins[(player,state)] = 0
      if i > max_depth
        max_depth = i
      end
   end

    push!(visited_states,(player,state))
    player = currentPlayer(states)
    winner = BM.winner(states[size(states)[1]])
    if winner!='?' #The current state has a winner
      #println(winner)
      break #Stop searching
    end

  end

  for (player,state) in visited_states
    if ! ( (player,state) in keys(plays))
      continue
    end
    winner = BM.winner(state)
    plays[(player,state)] += 1
    if player == winner
      wins[(player,state)] += 1
    end
  end
end

max_stars_repo_name: kinoroy/Shogi

max_issues_repo_name: kinoroy/Shulia

max_forks_repo_name: kinoroy/Shulia


text:
using TreeverseAlgorithm, Test, ForwardDiff
using NiLang
using NiLang.AD: GVar
@testset "treeverse" begin
    struct P3{T}
        x::T
        y::T
        z::T
    end

    Base.zero(::Type{P3{T}}) where T = P3(zero(T), zero(T), zero(T))
    Base.zero(::P3{T}) where T = P3(zero(T), zero(T), zero(T))

    @inline function Base.:(+)(a::P3, b::P3)
        P3(a.x + b.x, a.y + b.y, a.z + b.z)
    end

    @inline function Base.:(/)(a::P3, b::Real)
        P3(a.x/b, a.y/b, a.z/b)
    end

    @inline function Base.:(*)(a::Real, b::P3)
        P3(a*b.x, a*b.y, a*b.z)
    end

    function lorentz(t, y, θ)
        P3(10*(y.y-y.x), y.x*(27-y.z)-y.y, y.x*y.y-8/3*y.z)
    end

    function rk4_step(f, t, y, θ; Δt)
        k1 = Δt * f(t, y, θ)
        k2 = Δt * f(t+Δt/2, y + k1 / 2, θ)
        k3 = Δt * f(t+Δt/2, y + k2 / 2, θ)
        k4 = Δt * f(t+Δt, y + k3, θ)
        return y + k1/6 + k2/3 + k3/3 + k4/6
    end

    function rk4(f, y0::T, θ; t0, Δt, Nt) where T
        history = zeros(T, Nt+1)
        history[1] = y0
        y = y0
        for i=1:Nt
            y = rk4_step(f, t0+(i-1)*Δt, y, θ; Δt=Δt)
            @inbounds history[i+1] = y
        end
        return history
    end

    @i @inline function :(+=)(identity)(Y::P3, X::P3)
        Y.x += X.x
        Y.y += X.y
        Y.z += X.z
    end

    @i @inline function :(+=)(*)(Y::P3, a::Real, X::P3)
        Y.x += a * X.x
        Y.y += a * X.y
        Y.z += a * X.z
    end

    @i @inline function :(+=)(/)(Y::P3, X::P3, b::Real)
        Y.x += X.x/b
        Y.y += X.y/b
        Y.z += X.z/b
    end

    @i function lorentz!(y!::P3{T}, t, y::P3{T}, θ) where T
        @routine @invcheckoff begin
            @zeros T a b c b_a ab αc ac
            a += y.x
            b += y.y
            c += y.z
            b_a += b-a
            ab += a * b
            αc += (8/3) * c
            c -= 27
            ac += a * c
        end
        y!.x += 10 * b_a
        y!.y -= ac + b
        y!.z += ab - αc
        ~@routine
    end

    @i function rk4_step!(f, y!::T, y::T, θ; Δt, t) where T
        @routine @invcheckoff begin
            @zeros T k1 k2 k3 k4 o1 o2 o3 o4 yk1 yk2 yk3
            f(o1, t, y, θ)
            k1 += Δt * o1
            yk1 += y
            yk1 += k1 / 2
            t += Δt/2
            f(o2, t, yk1, θ)
            k2 += Δt * o2
            yk2 += y
            yk2 += k2 / 2
            f(o3, t, yk2, θ)
            k3 += Δt * o3
            yk3 += y
            yk3 += k3
            t += Δt/2
            f(o4, t, yk3, θ)
            k4 += Δt * o4
        end
        y! += y
        y! += k1 / 6
        y! += k2 / 3
        y! += k3 / 3
        y! += k4 / 6
        ~@routine
    end

    @i function rk4!(f, history, y0::T, θ; t0, Δt, Nt) where T
        history[1] += y0
        @invcheckoff @inbounds for i=1:Nt
            rk4_step!(f, history[i+1], history[i], θ; Δt=Δt, t=t0+(i-1)*Δt)
        end
    end

    @i function iloss!(out, f, history, y0, θ; t0, Δt, Nt)
        rk4!((@const f), history, y0, θ; t0=t0, Δt=Δt, Nt=Nt)
        out += history[end].x
    end
    @i function i_step_fun(state2, state)
        rk4_step!((@const lorentz!), state2.:2, state.:2, (); Δt=3e-3, t=state[1])
        state2.:1 += state.:1 + 3e-3
    end

    function step_fun(x)
        i_step_fun((0.0, zero(x[2])), x)[1]
    end

    function backward(x, g)
        if g===nothing
            g = (0.0, P3(1.0, 0.0, 0.0))
        end
        y = step_fun(x)
        _, gs = (~i_step_fun)(
            (GVar(y[1], g[1]), P3(GVar(y[2].x, g[2].x), GVar(y[2].y, g[2].y), GVar(y[2].z,g[2].z))),
            (GVar(x[1]), GVar(x[2])))
        NiLang.AD.grad(gs[1]), NiLang.AD.grad(gs[2])
    end

    @testset "treeverse gradient" begin
        x0 = P3(1.0, 0.0, 0.0)
        for N in [20, 120, 126]
            g_fd = ForwardDiff.gradient(x->rk4(lorentz, P3(x...), nothing; t0=0.0, Δt=3e-3, Nt=N)[end].x, [x0.x, x0.y, x0.z])
            g_tv = treeverse(step_fun, backward, (0.0, x0); δ=4, N=N)
            @test g_fd ≈ [g_tv[2].x, g_tv[2].y, g_tv[2].z]
        end
    end
end


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

JLLWrappers.@generate_main_file_header("basiclu")
JLLWrappers.@generate_main_file("basiclu", UUID("5629d00f-ba42-508e-84a8-8193befe9d4f"))
end  # module basiclu_jll


text:
## softmax

const CUDNN_SOFTMAX_FAST = (UInt32)(0)
const CUDNN_SOFTMAX_ACCURATE = (UInt32)(1)
const CUDNN_SOFTMAX_LOG = (UInt32)(2)

const CUDNN_SOFTMAX_MODE_INSTANCE = (UInt32)(0)
const CUDNN_SOFTMAX_MODE_CHANNEL = (UInt32)(1)

## batchnorm

const CUDNN_BATCHNORM_PER_ACTIVATION = (UInt32)(0)
const CUDNN_BATCHNORM_SPATIAL = (UInt32)(1)
const CUDNN_BN_MIN_EPSILON = 1.0e-5


text:
using IPIL8
using Test

@testset "IPIL8.jl" begin
    include("direct.jl");
end


text:
struct Opps

    #Coordinates
    xcoord      :: Array{Float64,1}
    ycoord      :: Array{Float64,1}

    aux_vex     :: Array{Float64,1}

    #Operators
    Dx          :: SparseMatrixCSC{Float64,Int64}
    Dxx         :: SparseMatrixCSC{Float64,Int64}
    Dy          :: SparseMatrixCSC{Float64,Int64}
    Dyy         :: SparseMatrixCSC{Float64,Int64}

    #RHS Differential Operator
    DiffO_1     :: SparseMatrixCSC{Float64,Int64}

    #Boundary Conditions
    BC_Mat      ::SparseMatrixCSC{Float64,Int64}

    #Absorption Matrix
    Abs_Mat      ::SparseMatrixCSC{Float64,Int64}

    #Angular momentum
    Lz      ::SparseMatrixCSC{Float64,Int64}

end


text:
# tag transforms from image

using YAML
using AprilTags, Images, ImageMagick
using Rotations, CoordinateTransformations
using StaticArrays
using MeshCat
using GeometryTypes

using PyCall
@pyimport numpy as np
@pyimport cv2

include(joinpath(Pkg.dir("Caesar"),"examples","wheeled","racecar","cameraUtils.jl"))
include(joinpath(Pkg.dir("Caesar"),"examples","wheeled","racecar","racecarUtils.jl"))
include(joinpath(Pkg.dir("Caesar"),"examples","wheeled","racecar","visualizationUtils.jl"))

cfg = loadConfig()

# fx = fy = 300.0
# cx,cy = (320,240)
k1,k2 = (0.0,0.0)
cx, cy = cfg[:intrinsics][:cx], cfg[:intrinsics][:cy]
fx = fy = cfg[:intrinsics][:fx]
# k1,k2 = cfg[:intrinsics][:k1], cfg[:intrinsics][:k2]

camK = [fx 0 cx; 0 fy cy; 0 0 1]
tagsize = 0.172

# tag extrinsic rotation
Rx = RotX(-pi/2)
Rz = RotZ(-pi/2)
bTc= LinearMap(Rz) ∘ LinearMap(Rx)

detector = AprilTagDetector()

image = load(joinpath(Pkg.dir("Caesar"),"examples","wheeled","racecar","example_tag_image_4.jpeg"))
tags = detector(image)

freeDetector!(detector)

imgd = deepcopy(image)
foreach(tag->drawTagBox!(imgd,tag, width = 5, drawReticle = false), tags)
imgd

# visualization

vis = Visualizer()
open(vis)

setobject!(vis["home"], Triad(1.0))

tidx = 1
currtag = ""
bTt = LinearMap(Translation(0.0,0,0) ∘ Quat(1.0,0,0,0))
for tidx in 1:length(tags)

currtag = "tag$(tags[tidx].id)"

# get relative tag transform
Q,T, = getAprilTagTransform(tags[tidx],
                          camK,
                          k1,
                          k2,
                          tagsize )
cTt = (T ∘ LinearMap(Q))

# tTl = LinearMap(RotY(pi/2)) ∘ LinearMap(RotX(-pi/2))
bTt = bTc ∘ cTt # ∘ tTl

# get Pose2Pose2 tag orientation transform
bP2t = getTagPP2(bTt)
# dt = convert(RotXYZ, bP2t.linear).theta3

drawTagDetection(vis, currtag, Q, T, bTc, bP2t, posename=:x5)

end

#


text:
module TestLmul

include("preamble.jl")
using LinearAlgebra

@testset begin
    @test lmul!!(1, 2) === 2

    A = LowerTriangular(reshape(1:4, 2, 2))
    B = ones(2, 2)
    AB = A * B
    @test lmul!!(A, B) === B
    @test B == AB
    @test_inferred lmul!!(A, B)

    A = LowerTriangular(collect(Float64, reshape(1:4, 2, 2)))
    B = ones(Int, 2, 2)
    AB = A * B
    @test lmul!!(A, B) :: Matrix{Float64} == AB
    @test_inferred lmul!!(A, B)
end

end  # module


text:
using Shapefile, GeoInterface
using Test

test_tuples = [
    (
        path="shapelib_testcases/test.shp",
        geomtype=Shapefile.Polygon,
        coordinates=Array{Array{Array{Array{Float64,1},1},1},1}[Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[20.0,20.0],[20.0,30.0],[30.0,30.0],[20.0,20.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[100.0,0.0],[100.0,100.0],[0.0,100.0],[0.0,0.0]]]],Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[150.0,150.0],[160.0,150.0],[180.0,170.0],[150.0,150.0]]]],Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[150.0,150.0],[160.0,150.0],[180.0,170.0],[150.0,150.0]]]]],
        bbox=Shapefile.Rect(0.0, 0.0, 180.0, 170.0),
        ),(
        path="shapelib_testcases/test0.shp",
        geomtype=Missing,
        coordinates=nothing,
        bbox=Shapefile.Rect(0.0, 0.0, 10.0, 20.0),
    ),(
        path="shapelib_testcases/test1.shp",
        geomtype=Shapefile.Point,
        coordinates=Array{Float64,1}[[1.0,2.0],[10.0,20.0]],
        bbox=Shapefile.Rect(1.0, 2.0, 10.0, 20.0),
    ),(
        path="shapelib_testcases/test2.shp",
        geomtype=Shapefile.PointZ,
        coordinates=Array{Float64,1}[[1.0,2.0,3.0],[10.0,20.0,30.0]],
        bbox=Shapefile.Rect(1.0, 2.0, 10.0, 20.0),
    ),(
        path="shapelib_testcases/test3.shp",
        geomtype=Shapefile.PointM,
        coordinates=Array{Float64,1}[[1.0,2.0],[10.0,20.0]],
        bbox=Shapefile.Rect(1.0, 2.0, 10.0, 20.0),
    ),(
        path="shapelib_testcases/test4.shp",
        geomtype=Shapefile.MultiPoint,
        coordinates=Array{Array{Float64,1},1}[Array{Float64,1}[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],Array{Float64,1}[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],Array{Float64,1}[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]],
        bbox=Shapefile.Rect(1.15, 2.25, 24.15, 25.25),
    ),(
        path="shapelib_testcases/test5.shp",
        geomtype=Shapefile.MultiPointZ,
        coordinates=Array{Array{Float64,1},1}[Array{Float64,1}[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],Array{Float64,1}[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],Array{Float64,1}[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]],
        bbox=Shapefile.Rect(1.15, 2.25, 24.15, 25.25),
    ),(
        path="shapelib_testcases/test6.shp",
        geomtype=Shapefile.MultiPointM,
        coordinates=Array{Array{Float64,1},1}[Array{Float64,1}[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],Array{Float64,1}[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],Array{Float64,1}[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]],
        bbox=Shapefile.Rect(1.15, 2.25, 24.15, 25.25),
    ),(
        path="shapelib_testcases/test7.shp",
        geomtype=Shapefile.Polyline,
        coordinates=Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]],
        bbox=Shapefile.Rect(0.0, 0.0, 100.0, 100.0),
    ),(
        path="shapelib_testcases/test8.shp",
        geomtype=Shapefile.PolylineZ,
        coordinates=Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]],
        bbox=Shapefile.Rect(0.0, 0.0, 100.0, 100.0),
    ),(
        path="shapelib_testcases/test9.shp",
        geomtype=Shapefile.PolylineM,
        coordinates=Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]],
        bbox=Shapefile.Rect(0.0, 0.0, 100.0, 100.0),
    ),(
        path="shapelib_testcases/test10.shp",
        geomtype=Shapefile.Polygon,
        coordinates=Array{Array{Array{Array{Float64,1},1},1},1}[Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]]],Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]]],Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]]],Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]]],
        bbox=Shapefile.Rect(0.0, 0.0, 100.0, 100.0),
    ),(
        path="shapelib_testcases/test11.shp",
        geomtype=Shapefile.PolygonZ,
        coordinates=Array{Array{Array{Array{Float64,1},1},1},1}[Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]]],Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]]],Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]]],Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]]],
        bbox=Shapefile.Rect(0.0, 0.0, 100.0, 100.0),
    ),(
        path="shapelib_testcases/test12.shp",
        geomtype=Shapefile.PolygonM,
        coordinates=Array{Array{Array{Array{Float64,1},1},1},1}[Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]]],Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]]],Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]]],Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]]],
        bbox=Shapefile.Rect(0.0, 0.0, 100.0, 100.0),
    ),(
        path="shapelib_testcases/test13.shp",
        geomtype=Shapefile.MultiPatch,
        coordinates=nothing,
        bbox=Shapefile.Rect(0.0, 0.0, 100.0, 100.0),
    )
]

@testset "Shapefile" begin

for test in test_tuples
    shp = open(joinpath(@__DIR__, test.path)) do fd
        read(fd, Shapefile.Handle)
    end
    shapes = unique(map(typeof, shp.shapes))
    @test length(shapes) == 1
    @test shapes[1] == test.geomtype
    @test eltype(shp.shapes) == Union{test.geomtype,Missing}
    # missing and MultiPatch are not covered by the GeoInterface
    if !(test.geomtype <: Union{Missing,Shapefile.MultiPatch})
        @test GeoInterface.coordinates.(shp.shapes) == test.coordinates
    end
    @test shp.MBR == test.bbox
end

# Test all .shx files; the values in .shx must match the .shp offsets
for test in test_tuples

    offsets = Int32[]
    contentlens = Int32[]

    # Get the shapefile's record offsets and contentlens
    shp = open(joinpath(@__DIR__, test.path)) do fd
        seek(fd, 32)
        shapeType = read(fd, Int32)
        seek(fd, 100)
        jltype = Shapefile.SHAPETYPE[shapeType]

        push!(offsets, position(fd))
        while (!eof(fd))

            num = bswap(read(fd, Int32))
            rlength = bswap(read(fd, Int32))
            shapeType = read(fd, Int32)
            if shapeType !== Int32(0)
                read(fd, jltype)
            end

            # records the offest after this geometry record
            push!(offsets, position(fd))
        end
    end
    contentlens = diff(offsets)
    offsets = offsets[1:end-1]

    # Match the Index values to Shapefile offsets
    shx =
        open(joinpath(@__DIR__, replace(test.path, r".shp$" => ".shx"))) do fd
            shx = read(fd, Shapefile.IndexHandle)
            for sIdx = 1:lastindex(shx.indices)
                @test shx.indices[sIdx].offset * 2 == offsets[sIdx]
                @test shx.indices[sIdx].contentLen * 2 + 8 == contentlens[sIdx]
            end
        end

end

include("table.jl")

end  # @testset "Shapefile"


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -16.957145096415093 <= q <= 17.94995800016443)

                     add_NL_constraint(m, :(sigmoid(0.510608367415756 + 0.13576703913191768*sigmoid(-0.023366318069798275 + -0.6972645804796413*sigmoid(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + 0.623708383159082*sigmoid(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.5526763214522168*sigmoid(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + 0.3585964327428499*sigmoid(-0.17574028037161238 + -0.933468976339102*sigmoid(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.824606977579692*sigmoid(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.3761458682864913*sigmoid(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + 0.2511872166134532*sigmoid(0.9643857894710246 + 0.9348568947078797*sigmoid(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.4027229019076284*sigmoid(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.6455757662896615*sigmoid(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3])))) + sigmoid(-0.7493006135688143 + 0.2483293387413994*sigmoid(-0.023366318069798275 + -0.6972645804796413*sigmoid(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + 0.623708383159082*sigmoid(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.5526763214522168*sigmoid(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + -0.7456041150082044*sigmoid(-0.17574028037161238 + -0.933468976339102*sigmoid(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.824606977579692*sigmoid(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.3761458682864913*sigmoid(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + 0.49415154189847854*sigmoid(0.9643857894710246 + 0.9348568947078797*sigmoid(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.4027229019076284*sigmoid(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.6455757662896615*sigmoid(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3])))) + sigmoid(0.3105901915386422 + 0.5384893864274902*sigmoid(-0.023366318069798275 + -0.6972645804796413*sigmoid(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + 0.623708383159082*sigmoid(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.5526763214522168*sigmoid(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + 0.5405413918049242*sigmoid(-0.17574028037161238 + -0.933468976339102*sigmoid(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.824606977579692*sigmoid(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.3761458682864913*sigmoid(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + -0.18067211353899282*sigmoid(0.9643857894710246 + 0.9348568947078797*sigmoid(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.4027229019076284*sigmoid(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.6455757662896615*sigmoid(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(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:
module TestSIMD

using FLoops
using Test

function test_at_simd()
    xs = 1:10
    @floop begin
        s = 0
        @simd for x in xs
            s += x
        end
    end
    @test s == sum(1:10)
end

function test_at_simd_ivdep()
    xs = 1:10
    ys = similar(xs)
    @floop begin
        @simd ivdep for (i, x) in enumerate(xs)
            @inbounds ys[i] = 2x
        end
    end
    @test ys == 2 .* xs
end

function test_Base_at_simd()
    xs = 1:10
    @floop begin
        s = 0
        Base.@simd for x in xs
            s += x
        end
    end
    @test s == sum(1:10)
end

function test_Base_at_simd_ivdep()
    xs = 1:10
    ys = similar(xs)
    @floop begin
        Base.@simd ivdep for (i, x) in enumerate(xs)
            @inbounds ys[i] = 2x
        end
    end
    @test ys == 2 .* xs
end

module TestFakeSIMD

using FLoops
using Test

macro simd(args...)
    esc(args[end])
end

function test_fake_at_simd()
    xs = 1:10
    @floop begin
        s = 0
        @simd for x in xs
            s += x
        end
    end
    @test s == sum(1:10)
end

function test_fake_at_simd_ivdep()
    xs = 1:10
    ys = similar(xs)
    @floop begin
        @simd ivdep for (i, x) in enumerate(xs)
            @inbounds ys[i] = 2x
        end
    end
    @test ys == 2 .* xs
end

end  # module TestFakeSIMD

function test_TestFakeSIMD_at_simd()
    xs = 1:10
    @floop begin
        s = 0
        TestFakeSIMD.@simd for x in xs
            s += x
        end
    end
    @test s == sum(1:10)
end

function test_TestFakeSIMD_at_simd_ivdep()
    xs = 1:10
    ys = similar(xs)
    @floop begin
        Base.@simd ivdep for (i, x) in enumerate(xs)
            @inbounds ys[i] = 2x
        end
    end
    @test ys == 2 .* xs
end

end  # module


text:
module Hiker
export answer

answer() = 6 * 9

end

max_stars_repo_name: cyber-dojo-start-points/julia-test

max_issues_repo_name: cyber-dojo-start-points/julia-test

max_forks_repo_name: cyber-dojo-start-points/julia-test


text:
function *(x::Sonum8,y::Sonum8)
    if signbit(x)
        if signbit(y)
            @inbounds return TableMul8S[Int(-x)+1,Int(-y)+1]
        else
            @inbounds return -TableMul8S[Int(-x)+1,Int(y)+1]
        end
    else
        if signbit(y)
            @inbounds return -TableMul8S[Int(x)+1,Int(-y)+1]
        else
            @inbounds return TableMul8S[Int(x)+1,Int(y)+1]
        end
    end
end

function +(x::Sonum8,y::Sonum8)
    if signbit(x)
        if signbit(y)   # -a-b = -(a+b)
            @inbounds return -TableAdd8S[Int(-x)+1,Int(-y)+1]
        else            # -a+b = b-a

            # anti-symmetric: check for size
            a = Int(y)+1
            b = Int(-x)+1

            if a > b
                @inbounds return -TableSub8S[b,a]
            else
                @inbounds return TableSub8S[a,b]
            end

        end
    else
        if signbit(y)   # a-b
            #return TableSub8[Int(x)+1,Int(-y)+1]

            # anti-symmetric: check for size
            a = Int(x)+1
            b = Int(-y)+1

            if a > b
                @inbounds return -TableSub8S[b,a]
            else
                @inbounds return TableSub8S[a,b]
            end
        else
            @inbounds return TableAdd8S[Int(x)+1,Int(y)+1]
        end
    end
end

function -(x::Sonum8,y::Sonum8)
    if signbit(x)
        if signbit(y)   # -a--b = b-a
            # return TableSub8[Int(-y)+1,Int(-x)+1]

            # anti-symmetric: check for size
            a = Int(-y)+1
            b = Int(-x)+1

            if a > b
                @inbounds return -TableSub8S[b,a]
            else
                @inbounds return TableSub8S[a,b]
            end

        else            # -a-b = -(a+b)
            @inbounds return -TableAdd8S[Int(-x)+1,Int(y)+1]
        end
    else
        if signbit(y)   # a--b = a+b
            @inbounds return TableAdd8S[Int(x)+1,Int(-y)+1]
        else            # a-b
            # return TableSub8[Int(x)+1,Int(y)+1]

            # anti-symmetric: check for size
            a = Int(x)+1
            b = Int(y)+1

            if a > b
                @inbounds return -TableSub8S[b,a]
            else
                @inbounds return TableSub8S[a,b]
            end
        end
    end
end

function /(x::Sonum8,y::Sonum8)
    if signbit(x)
        if signbit(y)
            @inbounds return TableDiv8[Int(-x)+1,Int(-y)+1]
        else
            @inbounds return -TableDiv8[Int(-x)+1,Int(y)+1]
        end
    else
        if signbit(y)
            @inbounds return -TableDiv8[Int(x)+1,Int(-y)+1]
        else
            @inbounds return TableDiv8[Int(x)+1,Int(y)+1]
        end
    end
end

# This increases the error as two roundings are applied
# function /(x::Sonum8,y::Sonum8)
#     return x*inv(y)
# end

function sqrt(x::Sonum8)
    if signbit(x)
        return notareal(Sonum8)
    else
        @inbounds return ListSqrt8[Int(x)+1]
    end
end

function inv(x::Sonum8)
    if signbit(x)
        @inbounds return -ListInv8[Int(-x)+1]
    else
        @inbounds return ListInv8[Int(x)+1]
    end
end


text:
using NonLinearReactionAdvectionDiffusionWithFrontData;
using NonLinearReactionAdvectionDiffusionWithFrontData: f1, f2;
using NonLinearReactionAdvectionDiffusionWithFrontData: ARP_y, ∂ARP_∂y;
using LinearAlgebra;

#' Сначала, сгенирируем априорные данные, на увеличенном количестве узлов.
u_l(t) = -8 #+ cos(2*π * t);
u_r(t) =  4 #+ (1 + sin(2*π * t));
qf(x) = 4*sin(3 * π * x);        # Коэффициент линейного усиления, который в обратной
                                # задаче необходимо определить, но при генерации априорной
                                # информации мы задаем некоторый коэффициент, который,
                                # собственно, после имея априорную информацию и будем определять.
ε = 0.2;                        # Малый параметр при старшей производной
a, b = 0, 1;                    # Область по X
t₀, T = 0.0, 0.28;              # Область по T
N, M = 50, 80;                  # Увеличенное Кол-во разбиений по X, T
h = (b-a)/N;                    # шаг по X
τ = (T-t₀)/M;                   # шаг по T
Xₙ = [a  + n*h for n in 0:N];   # Сетка по Х
Tₘ = [t₀ + m*τ for m in 0:M];   # Сетка по Т
qₙ =      qf.(Xₙ);               # Сеточные значения коэффициента линейного усиления
ulₘ=    u_l.(Tₘ);               # Сеточные значения левого  ГУ
urₘ=    u_r.(Tₘ);               # Сеточные значения правого ГУ
u₀ = u_init.(Xₙ);               # Начальные условия
#
u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);
nothing #hide

#' ## Генерация априорной информации
ϕl = phidetermination(qₙ, ulₘ, Xₙ, N, Tₘ, M);                       # Левый вырожденный корень
ϕr = phidetermination(qₙ, urₘ, Xₙ, N, Tₘ, M, reverseX = true);      # Правый вырожденный корень
ϕ = Φ(ϕl, ϕr, N, M);                                                # Серединный корень
f1_data = f1(ϕ, u, Xₙ, N, M);                                       # Положение переходного слоя
f2_data = f2(f1_data, u, Xₙ, N, M);                                 # Значение функции на переходном слое

# Подготовим массивы, выбросив граничные точки
# ведь тестируемая функция — для внутреннего использования
# И развернем массивы по времени
qq = qₙ[2:N];
y = u₀[2:N];
U = u[2:N, :];
U = reverse(U, dims=2);
f1_data = reverse(f1_data);
f2_data = reverse(f2_data);
Tₘ = reverse(Tₘ);

dl, d, du = NonLinearReactionAdvectionDiffusionWithFrontData.ARP_y(y, 1, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qq, U, f1_data, f2_data)

@test length(dl) == N - 2
@test length(d)  == N - 1
@test length(du) == N - 2

@test count(x -> x == 0, dl) == 0
@test count(x -> x == 0, d)  == 0
@test count(x -> x == 0, du) == 0

tdjac = NonLinearReactionAdvectionDiffusionWithFrontData.∂ARP_∂y(y, 1, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qq, U, f1_data, f2_data)
@test typeof(tdjac) <: Tridiagonal
@test tdjac == Tridiagonal( dl, d, du )

# Сравним с якобианом автодифференцирования
# Передадим **любой** априорный параметр для аппроксимации дельта-функции
# он никак не участвует в вычислении якобиана,
# но при использовании якобиана через автодифференцирование, требуется, чтобы он вызывался ТОЧНО так же
# как и функция, якобиан от которой вычисляется.
jac = NonLinearReactionAdvectionDiffusionWithFrontData.∂adjointRP_∂y(y, 1, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qq, U, f1_data, f2_data, 0.0001)
@test isapprox(tdjac, jac)


text:
using BinaryBuilder

name = "mongo-c-driver-builder"
version = v"1.14.0"

# Collection of sources required to build mongo-c-driver-builder
sources = [
    "https://zlib.net/zlib-1.2.11.tar.gz" =>
    "c3e5e9fdd5004dcb542feda5ee4f0ff0744628baf8ed2dd5d66f8ca1197cb1a1",

    "https://www.openssl.org/source/openssl-1.1.1b.tar.gz" =>
    "5c557b023230413dfb0756f3137a13e6d726838ccd1430888ad15bfb2b43ea4b",

    "https://github.com/cyrusimap/cyrus-sasl/releases/download/cyrus-sasl-2.1.27/cyrus-sasl-2.1.27.tar.gz" =>
    "26866b1549b00ffd020f188a43c258017fa1c382b3ddadd8201536f72efb05d5",

    "https://github.com/mongodb/mongo-c-driver/releases/download/1.14.0/mongo-c-driver-1.14.0.tar.gz" =>
    "ebe9694f7fa6477e594f19507877bbaa0b72747682541cf0cf9a6c29187e97e8",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd zlib-1.2.11/
./configure --prefix=$prefix
make -j${nproc}
make install
cd ..
cd openssl-1.1.1b/
./config --prefix=$prefix
make -j${nproc}
make install_sw
cd ..
cd cyrus-sasl-2.1.27/
./configure --prefix=$prefix --host=$target
make -j${nproc}
make install
cd ..
cd mongo-c-driver-1.14.0/
mkdir cmake-build
cd cmake-build/
cmake -DENABLE_AUTOMATIC_INIT_AND_CLEANUP=OFF -DCMAKE_BUILD_TYPE=Release -DENABLE_SSL=OPENSSL -DENABLE_STATIC=ON -DENABLE_MONGOC=ON -DENABLE_BSON=ON -DENABLE_ZLIB=BUNDLED -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain ..
make -j${nproc}
make install
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "padlock", :padlock),
    LibraryProduct(prefix, "libssl", :libssl),
    LibraryProduct(prefix, "afalg", :afalg),
    LibraryProduct(prefix, "capi", :capi),
    LibraryProduct(prefix, "libz", :libz),
    LibraryProduct(prefix, "libsasl2", :libsasl2),
    LibraryProduct(prefix, "libcrypto", :libcrypto),
    LibraryProduct(prefix, "libbson", :libbson),
    LibraryProduct(prefix, "libmongoc", :libmongoc),
    LibraryProduct(prefix, "libcrammd5", :libcrammd5),
    LibraryProduct(prefix, "libdigestmd5", :libdigestmd5),
    LibraryProduct(prefix, "libplain", :libplain),
    LibraryProduct(prefix, "libanonymous", :libanonymous)
]

# Dependencies that must be installed before this package can be built
dependencies = []

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

max_stars_repo_name: felipenoris/mongo-c-driver-builder

max_issues_repo_name: felipenoris/mongo-c-driver-builder

max_forks_repo_name: felipenoris/mongo-c-driver-builder


text:
module test2

using Test
using Example

@test Example.f() == 42

end


text:
#==============================================================================#
# DLM.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/dlm-2018-01-12.normal.json
#==============================================================================#

__precompile__()

module DLM

using AWSCore

"""
    using AWSSDK.DLM.create_lifecycle_policy
    create_lifecycle_policy([::AWSConfig], arguments::Dict)
    create_lifecycle_policy([::AWSConfig]; ExecutionRoleArn=, Description=, State=, PolicyDetails=)

    using AWSCore.Services.dlm
    dlm([::AWSConfig], "POST", "/policies", arguments::Dict)
    dlm([::AWSConfig], "POST", "/policies", ExecutionRoleArn=, Description=, State=, PolicyDetails=)

# CreateLifecyclePolicy Operation

Creates a policy to manage the lifecycle of the specified AWS resources. You can create up to 100 lifecycle policies.

# Arguments

## `ExecutionRoleArn = ::String` -- *Required*
The Amazon Resource Name (ARN) of the IAM role used to run the operations specified by the lifecycle policy.

## `Description = ::String` -- *Required*
A description of the lifecycle policy. The characters ^[0-9A-Za-z _-]+\$ are supported.

## `State = "ENABLED" or "DISABLED"` -- *Required*
The desired activation state of the lifecycle policy after creation.

## `PolicyDetails = [ ... ]` -- *Required*
The configuration of the lifecycle policy.

Target tags cannot be re-used across lifecycle policies.
```
 PolicyDetails = [
        "ResourceTypes" =>  ["VOLUME", ...],
        "TargetTags" =>  [[
            "Key" => <required> ::String,
            "Value" => <required> ::String
        ], ...],
        "Schedules" =>  [[
            "Name" =>  ::String,
            "TagsToAdd" =>  [[
                "Key" => <required> ::String,
                "Value" => <required> ::String
            ], ...],
            "CreateRule" =>  [
                "Interval" => <required> ::Int,
                "IntervalUnit" => <required> "HOURS",
                "Times" =>  [::String, ...]
            ],
            "RetainRule" =>  ["Count" => <required> ::Int]
        ], ...]
    ]
```

# Returns

`CreateLifecyclePolicyResponse`

# Exceptions

`InvalidRequestException`, `LimitExceededException` or `InternalServerException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/dlm-2018-01-12/CreateLifecyclePolicy)
"""
@inline create_lifecycle_policy(aws::AWSConfig=default_aws_config(); args...) = create_lifecycle_policy(aws, args)

@inline create_lifecycle_policy(aws::AWSConfig, args) = AWSCore.Services.dlm(aws, "POST", "/policies", args)

@inline create_lifecycle_policy(args) = create_lifecycle_policy(default_aws_config(), args)

"""
    using AWSSDK.DLM.delete_lifecycle_policy
    delete_lifecycle_policy([::AWSConfig], arguments::Dict)
    delete_lifecycle_policy([::AWSConfig]; policyId=)

    using AWSCore.Services.dlm
    dlm([::AWSConfig], "DELETE", "/policies/{policyId}/", arguments::Dict)
    dlm([::AWSConfig], "DELETE", "/policies/{policyId}/", policyId=)

# DeleteLifecyclePolicy Operation

Deletes the specified lifecycle policy and halts the automated operations that the policy specified.

# Arguments

## `policyId = ::String` -- *Required*
The identifier of the lifecycle policy.

# Returns

`DeleteLifecyclePolicyResponse`

# Exceptions

`ResourceNotFoundException`, `InternalServerException` or `LimitExceededException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/dlm-2018-01-12/DeleteLifecyclePolicy)
"""
@inline delete_lifecycle_policy(aws::AWSConfig=default_aws_config(); args...) = delete_lifecycle_policy(aws, args)

@inline delete_lifecycle_policy(aws::AWSConfig, args) = AWSCore.Services.dlm(aws, "DELETE", "/policies/{policyId}/", args)

@inline delete_lifecycle_policy(args) = delete_lifecycle_policy(default_aws_config(), args)

"""
    using AWSSDK.DLM.get_lifecycle_policies
    get_lifecycle_policies([::AWSConfig], arguments::Dict)
    get_lifecycle_policies([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.dlm
    dlm([::AWSConfig], "GET", "/policies", arguments::Dict)
    dlm([::AWSConfig], "GET", "/policies", <keyword arguments>)

# GetLifecyclePolicies Operation

Gets summary information about all or the specified data lifecycle policies.

To get complete information about a policy, use [GetLifecyclePolicy](@ref).

# Arguments

## `policyIds = [::String, ...]`
The identifiers of the data lifecycle policies.

## `state = "ENABLED", "DISABLED" or "ERROR"`
The activation state.

## `resourceTypes = ["VOLUME", ...]`
The resource type.

## `targetTags = [::String, ...]`
The target tags.

Tags are strings in the format `key:value`.

## `tagsToAdd = [::String, ...]`
The tags to add to the resources.

Tags are strings in the format `key:value`.

These tags are added in addition to the AWS-added lifecycle tags.

# Returns

`GetLifecyclePoliciesResponse`

# Exceptions

`ResourceNotFoundException`, `InvalidRequestException`, `InternalServerException` or `LimitExceededException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/dlm-2018-01-12/GetLifecyclePolicies)
"""
@inline get_lifecycle_policies(aws::AWSConfig=default_aws_config(); args...) = get_lifecycle_policies(aws, args)

@inline get_lifecycle_policies(aws::AWSConfig, args) = AWSCore.Services.dlm(aws, "GET", "/policies", args)

@inline get_lifecycle_policies(args) = get_lifecycle_policies(default_aws_config(), args)

"""
    using AWSSDK.DLM.get_lifecycle_policy
    get_lifecycle_policy([::AWSConfig], arguments::Dict)
    get_lifecycle_policy([::AWSConfig]; policyId=)

    using AWSCore.Services.dlm
    dlm([::AWSConfig], "GET", "/policies/{policyId}/", arguments::Dict)
    dlm([::AWSConfig], "GET", "/policies/{policyId}/", policyId=)

# GetLifecyclePolicy Operation

Gets detailed information about the specified lifecycle policy.

# Arguments

## `policyId = ::String` -- *Required*
The identifier of the lifecycle policy.

# Returns

`GetLifecyclePolicyResponse`

# Exceptions

`ResourceNotFoundException`, `InternalServerException` or `LimitExceededException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/dlm-2018-01-12/GetLifecyclePolicy)
"""
@inline get_lifecycle_policy(aws::AWSConfig=default_aws_config(); args...) = get_lifecycle_policy(aws, args)

@inline get_lifecycle_policy(aws::AWSConfig, args) = AWSCore.Services.dlm(aws, "GET", "/policies/{policyId}/", args)

@inline get_lifecycle_policy(args) = get_lifecycle_policy(default_aws_config(), args)

"""
    using AWSSDK.DLM.update_lifecycle_policy
    update_lifecycle_policy([::AWSConfig], arguments::Dict)
    update_lifecycle_policy([::AWSConfig]; policyId=, <keyword arguments>)

    using AWSCore.Services.dlm
    dlm([::AWSConfig], "PATCH", "/policies/{policyId}", arguments::Dict)
    dlm([::AWSConfig], "PATCH", "/policies/{policyId}", policyId=, <keyword arguments>)

# UpdateLifecyclePolicy Operation

Updates the specified lifecycle policy.

# Arguments

## `policyId = ::String` -- *Required*
The identifier of the lifecycle policy.

## `ExecutionRoleArn = ::String`
The Amazon Resource Name (ARN) of the IAM role used to run the operations specified by the lifecycle policy.

## `State = "ENABLED" or "DISABLED"`
The desired activation state of the lifecycle policy after creation.

## `Description = ::String`
A description of the lifecycle policy.

## `PolicyDetails = [ ... ]`
The configuration of the lifecycle policy.

Target tags cannot be re-used across policies.
```
 PolicyDetails = [
        "ResourceTypes" =>  ["VOLUME", ...],
        "TargetTags" =>  [[
            "Key" => <required> ::String,
            "Value" => <required> ::String
        ], ...],
        "Schedules" =>  [[
            "Name" =>  ::String,
            "TagsToAdd" =>  [[
                "Key" => <required> ::String,
                "Value" => <required> ::String
            ], ...],
            "CreateRule" =>  [
                "Interval" => <required> ::Int,
                "IntervalUnit" => <required> "HOURS",
                "Times" =>  [::String, ...]
            ],
            "RetainRule" =>  ["Count" => <required> ::Int]
        ], ...]
    ]
```

# Returns

`UpdateLifecyclePolicyResponse`

# Exceptions

`ResourceNotFoundException`, `InvalidRequestException`, `InternalServerException` or `LimitExceededException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/dlm-2018-01-12/UpdateLifecyclePolicy)
"""
@inline update_lifecycle_policy(aws::AWSConfig=default_aws_config(); args...) = update_lifecycle_policy(aws, args)

@inline update_lifecycle_policy(aws::AWSConfig, args) = AWSCore.Services.dlm(aws, "PATCH", "/policies/{policyId}", args)

@inline update_lifecycle_policy(args) = update_lifecycle_policy(default_aws_config(), args)

end # module DLM

#==============================================================================#
# End of file
#==============================================================================#


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -7.9284408379108005 <= q <= 14.377266175441108)

                     add_NL_constraint(m, :(gelu(0.4343485145630286 + -0.5735104656090533*gelu(-0.47832195115882303 + 0.6358617391786257*$(x[1]) + -0.4037373919488174*$(x[2])) + -0.40833038450807946*gelu(0.4113265195541662 + -0.7994001067331897*$(x[1]) + -0.3401124060181533*$(x[2])) + 0.09714914815406805*gelu(-0.3249145933093094 + 0.9313960721440679*$(x[1]) + 0.2643917347767659*$(x[2])) + -0.19623354791767778*gelu(0.34085895078435247 + -0.33650259644885194*$(x[1]) + -0.9436140496760688*$(x[2])) + -0.3237276355524821*gelu(0.8342164127508305 + -0.6721074785776127*$(x[1]) + 0.009028597341863787*$(x[2]))) + gelu(-0.11905941659917563 + 0.17224590981590104*gelu(-0.47832195115882303 + 0.6358617391786257*$(x[1]) + -0.4037373919488174*$(x[2])) + 0.8406566605626429*gelu(0.4113265195541662 + -0.7994001067331897*$(x[1]) + -0.3401124060181533*$(x[2])) + 0.24719777641241691*gelu(-0.3249145933093094 + 0.9313960721440679*$(x[1]) + 0.2643917347767659*$(x[2])) + -0.24308200319021767*gelu(0.34085895078435247 + -0.33650259644885194*$(x[1]) + -0.9436140496760688*$(x[2])) + -0.39557038452022475*gelu(0.8342164127508305 + -0.6721074785776127*$(x[1]) + 0.009028597341863787*$(x[2]))) + gelu(0.967203829156619 + 0.3785230261408685*gelu(-0.47832195115882303 + 0.6358617391786257*$(x[1]) + -0.4037373919488174*$(x[2])) + -0.23118999327456224*gelu(0.4113265195541662 + -0.7994001067331897*$(x[1]) + -0.3401124060181533*$(x[2])) + -0.5197073892915105*gelu(-0.3249145933093094 + 0.9313960721440679*$(x[1]) + 0.2643917347767659*$(x[2])) + 0.6005003791404913*gelu(0.34085895078435247 + -0.33650259644885194*$(x[1]) + -0.9436140496760688*$(x[2])) + 0.15463985115916623*gelu(0.8342164127508305 + -0.6721074785776127*$(x[1]) + 0.009028597341863787*$(x[2]))) + gelu(-0.18340968019959902 + 0.01980109053282808*gelu(-0.47832195115882303 + 0.6358617391786257*$(x[1]) + -0.4037373919488174*$(x[2])) + 0.8767310690037822*gelu(0.4113265195541662 + -0.7994001067331897*$(x[1]) + -0.3401124060181533*$(x[2])) + 0.42243485881005904*gelu(-0.3249145933093094 + 0.9313960721440679*$(x[1]) + 0.2643917347767659*$(x[2])) + -0.1694649191201152*gelu(0.34085895078435247 + -0.33650259644885194*$(x[1]) + -0.9436140496760688*$(x[2])) + 0.5278959282338329*gelu(0.8342164127508305 + -0.6721074785776127*$(x[1]) + 0.009028597341863787*$(x[2]))) + gelu(0.7043455967116397 + -0.48930630508951234*gelu(-0.47832195115882303 + 0.6358617391786257*$(x[1]) + -0.4037373919488174*$(x[2])) + 0.7375930151357304*gelu(0.4113265195541662 + -0.7994001067331897*$(x[1]) + -0.3401124060181533*$(x[2])) + -0.9188673755139103*gelu(-0.3249145933093094 + 0.9313960721440679*$(x[1]) + 0.2643917347767659*$(x[2])) + 0.5801272291202126*gelu(0.34085895078435247 + -0.33650259644885194*$(x[1]) + -0.9436140496760688*$(x[2])) + 0.06743755305487698*gelu(0.8342164127508305 + -0.6721074785776127*$(x[1]) + 0.009028597341863787*$(x[2]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m


text:
using LibDeflate
using Test
using CodecZlib

@testset "DEFLATE" begin
    include("deflate.jl")
end

@testset "gzip" begin
    include("gzip.jl")
end

@testset "zlib" begin
    include("zlib.jl")
end


text:
module ElasticAnisotropy

export getRmatrix, OrthoModel, TransIsomodel, IsoModel, rotateTensor, Tensor2VoigtMatrix, VoigtMatrix2Tensor, generateVoigtMatrix

using LinearAlgebra

struct OrthoModel
    moduli::NamedTuple{(:E_1, :E_2, :E_3, :ν_12, :ν_23, :ν_13, :G_12, :G_23, :G_13), Tuple{Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64}}
end

struct TransIsoModel
    moduli::NamedTuple{(:E_i,:E_a,:ν_i,:ν_ia,:G_a), Tuple{Float64,Float64,Float64,Float64,Float64}}
end

struct IsoModel
    moduli::NamedTuple{(:E,:ν), Tuple{Float64,Float64}}
end

function generateVoigtMatrix(p::OrthoModel)
    E_1, E_2, E_3, ν_12, ν_23, ν_13, G_12, G_23, G_13 = p.moduli
    ν_21 = ν_12 * E_2 / E_1
    ν_32 = ν_23 * E_3 / E_2
    ν_31 = ν_13 * E_3 / E_1
    D = (1 - ν_12*ν_21-ν_23*ν_32-ν_31*ν_13-2*ν_12*ν_23*ν_31)/(E_1 * E_2 * E_3)
    C = [(1-ν_23*ν_32)/(E_2*E_3*D) (ν_21+ν_31*ν_23)/(E_2*E_3*D) (ν_31+ν_21*ν_32)/(E_2*E_3*D) 0 0 0; (ν_12+ν_13*ν_32)/(E_1*E_3*D) (1-ν_31*ν_13)/(E_1*E_3*D) (ν_32+ν_31*ν_12)/(E_1*E_3*D) 0 0 0; (ν_13+ν_12*ν_23)/(E_1*E_2*D) (ν_23+ν_13*ν_21)/(E_1*E_2*D) (1-ν_12*ν_21)/(E_1*E_2*D) 0 0 0; 0 0 0 G_23 0 0; 0 0 0 0 G_13 0; 0 0 0 0 0 G_12]
    return C
end

function generateVoigtMatrix(p::TransIsoModel)
    E_i, E_a, ν_i, ν_ia, G_a = p.moduli
    G_i = E_i / (2.0*(1. + ν_i))
    ν_ai = ν_ia * E_a / E_i
    D = ((1.0+ν_i)*(1.0-ν_i-2*ν_ia * ν_ai))/(E_i^2*E_a)
    a_ii = (1-ν_ia*ν_ai)/(E_i*E_a*D)
    a_ai = (1-ν_i^2)/(E_i^2*D)
    b_ii = (ν_i+ν_ia*ν_ai)/(E_i*E_a*D)
    b_ai = (ν_ia*(1+ν_i))/(E_i^2*D)
    c_ii = E_i/(2*(1+ν_i))
    c_ai = G_a
    C = [a_ii b_ii b_ai 0 0 0; b_ii a_ii b_ai 0 0 0; b_ai b_ai a_ai 0 0 0; 0 0 0 c_ai 0 0; 0 0 0 0 c_ai 0; 0 0 0 0 0 c_ii]
    return C
end

function generateVoigtMatrix(p::IsoModel)
    E, ν = p.moduli
    prefac = E/((1+ν)*(1-2*ν))
    a_ii = 1-ν
    a_ij = ν
    b_ii = (1-2*ν)/2
    C = prefac .* [a_ii a_ij a_ij 0 0 0; a_ij a_ii a_ij 0 0 0; a_ij a_ij a_ii 0 0 0; 0 0 0 b_ii 0 0; 0 0 0 0 b_ii 0; 0 0 0 0 0 b_ii]
    return C
end

function getCijkl(C,i,j,k,l)
    function getvoigtindex(i,j)
        if i == j
            return i
        elseif (i == 2 && j==3) || (i == 3 && j==2)
            return 4
        elseif (i == 1 && j==3) || (i == 3 && j==1)
            return 5
        elseif (i == 1 && j==2) || (i == 2 && j==1)
            return 6
        end
    end
    m = getvoigtindex(i,j)
    n = getvoigtindex(k,l)
    return C[m,n]
end

function getCij(c,i,j)
    function gettensorindex(i)
        if (i >0) && (i < 4)
            return i, i
        elseif i == 4
            return 2, 3
        elseif i == 5
            return 1, 3
        elseif i == 6
            return 1, 2
        end
    end
    k, l = gettensorindex(i)
    m, n = gettensorindex(j)
    return c[k,l,m,n]
end

function VoigtMatrix2Tensor(C)
    cijkl=zeros((3,3,3,3));
    for i in 1:3
        for j in 1:3
            for k in 1:3
                for l in 1:3
                    cijkl[i,j,k,l] = getCijkl(C,i, j,k,l)
                end
            end
        end
    end
    return cijkl
end

function transform_tensor(c,R,i,j,k,l)
    c_new = 0.0
    for r in 1:3
        for s in 1:3
            for t in 1:3
                for u in 1:3
                    c_new +=R[i,r]*R[j,s]*R[k,t]*R[l,u]*c[r,s,t,u]
                end
            end
        end
    end
    return c_new
end

function getRmatrix(n)
    e3 = [n[1], n[2], n[3]]
    e1 = [n[3]/(n[1]*sqrt(1+(n[3]/n[1])^2)), 0, -1.0/sqrt(1+(n[3]/n[1])^2)]
    e2 = cross(e3,e1)
    R = [e1[1] e2[1] e3[1]; e1[2] e2[2] e3[2]; e1[3] e2[3] e3[3]]
    return R
end

function rotateTensor(cijkl, n::Array{Float64,1})
    function getRmatrix(n)
        e3 = [n[1], n[2], n[3]]
        e1 = [n[3]/(n[1]*sqrt(1+(n[3]/n[1])^2)), 0, -1.0/sqrt(1+(n[3]/n[1])^2)]
        e2 = cross(e3,e1)
        R = [e1[1] e2[1] e3[1]; e1[2] e2[2] e3[2]; e1[3] e2[3] e3[3]]
        return R
    end
    R=getRmatrix(n)
    return rotateTensor(cijkl,R)
end

function rotateTensor(cijkl, R::Array{Float64,2})
    cijkl_new=zeros((3,3,3,3));
    for i in 1:3
        for j in 1:3
            for k in 1:3
                for l in 1:3
                    cijkl_new[i,j,k,l] = transform_tensor(cijkl,R,i,j,k,l)
                end
            end
        end
    end
    return cijkl_new
end

function Tensor2VoigtMatrix(cijkl)
    C_new = zeros(6,6)
    for j in 1:6
        for i in 1:6
            C_new[i,j] = getCij(cijkl,i,j)
        end
    end
    return C_new
end

end


text:
__precompile__()

module MathProgBase

include(joinpath(dirname(@__FILE__),"SolverInterface","SolverInterface.jl"))
using .SolverInterface

include(joinpath(dirname(@__FILE__),"HighLevelInterface","HighLevelInterface.jl"))
using .HighLevelInterface
export linprog, mixintprog, quadprog, buildlp, solvelp

end


text:
using Plots
plotly()
using brainvisionloader

filename = "../data/example"
EEG = brainvisionloader.loadFile(filename)

toi     = range(EEG.markers[2,3]-EEG.srate,stop=EEG.markers[2,3]+EEG.srate);
y       = EEG.data[25,Int64.(toi)]
x       = collect(1:length(y)) .-EEG.srate
xtick   = -EEG.srate:(EEG.srate/4):EEG.srate

plot(x, y,
    linewidth=2,
    alpha=0.6,
    title="EEG Trace",
    legend=:none,
    xlabel = "Time [ms]",
    ylabel = "Amplitude [μV]",
    size = (1024,480),
    xticks = xtick,
    );
gui()
savefig("./example.png")


text:
using GPUCompiler

if !@isdefined(TestRuntime)
    include("../util.jl")
end

# create a native test compiler, and generate reflection methods for it

function native_job(@nospecialize(func), @nospecialize(types); kernel::Bool=false, kwargs...)
    source = FunctionSpec(func, Base.to_tuple_type(types), kernel)
    target = NativeCompilerTarget()
    params = TestCompilerParams()
    CompilerJob(target, source, params), kwargs
end

function native_code_typed(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kwargs...)
    GPUCompiler.code_typed(job; kwargs...)
end

function native_code_warntype(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kwargs...)
    GPUCompiler.code_warntype(io, job; kwargs...)
end

function native_code_llvm(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kwargs...)
    GPUCompiler.code_llvm(io, job; kwargs...)
end

function native_code_native(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kwargs...)
    GPUCompiler.code_native(io, job; kwargs...)
end

# aliases without ::IO argument
for method in (:code_warntype, :code_llvm, :code_native)
    native_method = Symbol("native_$(method)")
    @eval begin
        $native_method(@nospecialize(func), @nospecialize(types); kwargs...) =
            $native_method(stdout, func, types; kwargs...)
    end
end

# simulates codegen for a kernel function: validates by default
function native_code_execution(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kernel=true, kwargs...)
    GPUCompiler.compile(:asm, job; kwargs...)
end

module LazyCodegen
    using LLVM
    using GPUCompiler

    import ..native_job

    # We have one global JIT and TM
    const jit = Ref{LLVM.OrcJIT}()
    const tm  = Ref{LLVM.TargetMachine}()

    function __init__()
        optlevel = LLVM.API.LLVMCodeGenLevelDefault

        tm[] = GPUCompiler.JITTargetMachine(optlevel=optlevel)
        LLVM.asm_verbosity!(tm[], true)

        jit[] = LLVM.OrcJIT(tm[]) # takes ownership of tm
        atexit() do
            dispose(jit[])
        end
    end

    import GPUCompiler: deferred_codegen_jobs, CompilerJob
    mutable struct CallbackContext
        job::CompilerJob
        stub::String
        compiled::Bool
    end

    const outstanding = IdDict{CallbackContext, Nothing}()

    # Setup the lazy callback for creating a module
    function callback(orc_ref::LLVM.API.LLVMOrcJITStackRef, callback_ctx::Ptr{Cvoid})
        orc = LLVM.OrcJIT(orc_ref)
        cc = Base.unsafe_pointer_to_objref(callback_ctx)::CallbackContext

        @assert !cc.compiled
        job = cc.job

        ir, entry = GPUCompiler.codegen(:llvm, job; validate=false)
        entry_name = name(entry)

        jitted_mod = compile!(orc, ir)

        addr = addressin(orc, jitted_mod, entry_name)
        ptr  = pointer(addr)

        cc.compiled = true
        delete!(outstanding, cc)

        # 4. Update the stub pointer to point to the recently compiled module
        set_stub!(orc, cc.stub, ptr)

        # 5. Return the address of the implementation, since we are going to call it now
        return UInt64(reinterpret(UInt, ptr))
    end

    @generated function deferred_codegen(::Val{f}, ::Val{tt}) where {f,tt}
        job, _ = native_job(f, tt)

        cc = CallbackContext(job, String(gensym(:trampoline)), false)
        outstanding[cc] = nothing

        c_callback = @cfunction(callback, UInt64, (LLVM.API.LLVMOrcJITStackRef, Ptr{Cvoid}))

        orc = jit[]
        initial_addr = callback!(orc, c_callback, pointer_from_objref(cc))
        create_stub!(orc, cc.stub, initial_addr)
        addr = address(orc, cc.stub)
        trampoline = pointer(addr)
        id = Base.reinterpret(Int, trampoline)

        deferred_codegen_jobs[id] = job

        quote
            ccall("extern deferred_codegen", llvmcall, Ptr{Cvoid}, (Ptr{Cvoid},), $trampoline)
        end
    end
end


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

"""
    BAT.MultimodalCauchy <: Distribution{Multivariate,Continuous}

*Experimental feature, not part of stable public API.*

The Multimodal Cauchy Distribution (see
[Caldwell et al.](https://arxiv.org/abs/1808.08051) for definition).

Assumes two bimodal peaks, each in its own dimension.

Constructors:

* ```BAT.MultimodalCauchy(; μ::Real=1, σ::Float64=0.2, n::Integer=4)```

Arguments:

* `μ::Real`: The location parameter used for the two bimodal peaks.
* `σ::Float64`: The scale parameter shared among all components.
* `n::Int`: The number of dimensions.

Fields:

$(TYPEDFIELDS)

!!! note

    All fields of `MultimodalCauchy` are considered internal and subject to
    change without deprecation.
"""
struct MultimodalCauchy{M<:MixtureModel, P<:Product} <: Distribution{Multivariate,Continuous}
    bimodals::M
    σ::Float64
    n::Int
    dist::P
end

function MultimodalCauchy(;μ::Real=1, σ::Float64=0.2, n::Integer=4)
    @argcheck n > 1 "Minimum number of dimensions for MultimodalCauchy is 2"
    mixture_model = MixtureModel([Cauchy(-μ, σ), Cauchy(μ, σ)])
    dist = _construct_dist(mixture_model, σ, n)
    MultimodalCauchy(mixture_model, σ, n, dist)
end

Base.size(d::MultimodalCauchy) = size(d.dist)
Base.length(d::MultimodalCauchy) = length(d.dist)
Base.eltype(d::MultimodalCauchy) = eltype(d.dist)

Statistics.mean(d::MultimodalCauchy) = Distributions.mean(d.dist)
Statistics.var(d::MultimodalCauchy) = Distributions.var(d.dist)
Statistics.cov(d::MultimodalCauchy) = Distributions.cov(d.dist)

function StatsBase.params(d::MultimodalCauchy)
    (
        vcat(d.bimodals.components[1].μ, d.bimodals.components[2].μ, zeros(d.n-2)),
        [d.σ for i in 1:d.n],
        d.n
    )
end

function Distributions._logpdf(d::MultimodalCauchy, x::AbstractArray)
    Distributions._logpdf(d.dist, x)
end

function Distributions._rand!(rng::AbstractRNG, d::MultimodalCauchy, x::AbstractVector)
    Distributions._rand!(rng, d.dist, x)
end

function _construct_dist(mixture_model::MixtureModel, σ::Real, n::Integer)
    vector_of_dists = vcat(mixture_model, mixture_model, [Cauchy(0, σ) for i in 3:n])
    dist = product_distribution(vector_of_dists)
    return dist
end


text:
function calculate_u1u2_spectrum!(Ef,u,cplane)
    ux = u.c.x
    uy = u.c.y
    # Initialize the shells to zeros
    fill!(Ef,0)
    maxdk2d = max(KX[2],KY[2])
    nshells = min(NX,NY÷2)
    @inbounds for j in YRANGE
        conjFactX=1.0
        for i in XRANGE
            k = fsqrt(KX[i]^2 + KY[j]^2)
            n = round(Int, k/maxdk2d) + 1
            if n <= nshells
                magsq = abs2(ux[i,j,cplane]) + abs2(uy[i,j,cplane])
                ee = 0.5*conjFactX * magsq;
                Ef[n]+=ee;
            end
            conjFactX=2.0;
        end
    end
end

function calculate_u1u2_spectrum(u,cplane,s::AbstractSimulation)
    Ef = zeros(min(NX, NY÷2))
    calculate_u1u2_spectrum!(Ef,u,1)
    return Ef
end

function calculate_abs2!(Ef,u,cplane)
    ux = u.c.x
    uy = u.c.y
    # Initialize the shells to zeros
    fill!(Ef,0)
    @inbounds for j in YRANGE
        conjFactX=1.0
        for i in XRANGE
            magsq = abs2(ux[i,j,cplane]) + abs2(uy[i,j,cplane])
            ee = 0.5*conjFactX * magsq;
            Ef[i,j]=ee;
            conjFactX=2.0;
        end
    end
end

function calculate_u1u2_mag_spectrum!(Ef,u,cplane)
    ux = u.c.x
    uy = u.c.y
    # Initialize the shells to zeros
    fill!(Ef,0)
    maxdk2d = max(KX[2],KY[2])
    nshells = min(NX,NY÷2)
    @inbounds for j in YRANGE
        conjFactX=1.0
        for i in XRANGE
            k = fsqrt(KX[i]^2 + KY[j]^2)
            n = round(Int, k/maxdk2d) + 1
            if n <= nshells
                mag = fsqrt(abs2(ux[i,j,cplane]) + abs2(uy[i,j,cplane]))
                ee = conjFactX * mag;
                Ef[n]+=ee;
            end
            conjFactX=2.0;
        end
    end
end


text:
using GeometryBasics
using TriangleIntersect
using Test

a = Point(0, 0, 0)
b = Point(1, 0, 0)
c = Point(0, 1, 0)
t = Triangle(a, b, c)

#=
         45° .
           ⟋ |
s = d√2  ⟋   |
       ⟋     | d
     ⟋       |
45°⟋_________|
        d

Fig. 1: Diagram of a 45-45-90 triangle.

This is a top view.
The top corner is (x0, y0, 0), the vertical projection of the initial point (x0, y0, z0) onto the triange t
The hypotenuse is the path of the ray as viewed from above.
The bottom left corner is the intersection point of the ray and triangle. =#

#=
    . 30°
    |\
    | \
    |  \  2s
s√3 |   \
    |    \
    |_____\ 60°
       s

Fig. 2: Diagram of a 30-60-90 triangle.

This is a side view.
The top corner is the initial point (x0, y0, z0).
The bottom side is the same as the hypotenuse of Fig. 1
The intersection distance is the length of the hypotenuse. =#

d = 0.1
s = √2d
(x0, y0, z0) = (d, d, √3s)
# Scale the Vec just to make sure normalization is working.
r = Ray(Point(x0, y0, z0), .12345 * Vec(d, d, -√3s))
@test intersect(r, t) == Intersection(Point(x0 + d, y0 + d, 0), 2s, true)

# yeah, more tests are needed :)


text:
using VARmodels
using Base.Test

# write your own tests here
i = rand(1:2000)
run(`Rscript varest.R $i`)
C = readcsv("coefs.csv")
S = readcsv("sigma.csv")
resids = readcsv("residuals.csv")

data = readcsv("data.csv")
data = log(data)
est = varEstimate(data[(1+i):(500+i),:], 2, "Const")
residuals = (est.Y-est.X*est.C)
@test all(abs(C[2:size(C)[1],:]-est.C).<1e-10)
@test all(abs(residuals-resids[2:size(resids)[1],:]).<1e-10)
@test all(abs(S[2:size(S)[1],:]-est.Σ).<1e-10)
run(`rm coefs.csv`)
run(`rm residuals.csv`)
run(`rm sigma.csv`)


text:
using PkgLicenses
using Test

# Make sure that we actually have all the licenses we claim to
for lic in keys(LICENSES)
    @test isa(readlicense(lic), String)
end


text:
export AstrobeeSE3
export init_traj_straightline

mutable struct AstrobeeSE3 <: DynamicsModel
  # state: r v p ω
  x_dim
  u_dim
  clearance

  # Parameters that can be updated
  f::Vector
  A::Vector
  B
end

function AstrobeeSE3()
  x_dim = 12
  u_dim = 6
  clearance = 0.03
  AstrobeeSE3(x_dim, u_dim, clearance, [], [], [])
end

function SCPParam(model::AstrobeeSE3, fixed_final_time::Bool)
  convergence_threshold = 0.01
  SCPParam(fixed_final_time, convergence_threshold)
end

function SCPParam_GuSTO(model::AstrobeeSE3)
  Δ0 = 10.
  ω0 = 1.
  ω_max = 1.0e10
  ε = 1.0e-6
  ρ0 = 0.01
  ρ1 = 0.05
  β_succ = 2.
  β_fail = 0.5
  γ_fail = 5.

  SCPParam_GuSTO(Δ0, ω0, ω_max, ε, ρ0, ρ1, β_succ, β_fail, γ_fail)
end

function SCPParam_Mao(model::AstrobeeSE3)
  ρ = [0.; 0.25; 0.9]
  Δ_u0 = 1000.
  λ = 1.
  α = 2.
  SCPParam_Mao(ρ, Δ_u0, λ, α)
end

function SCPParam_TrajOpt(model::AstrobeeSE3)
  μ0 = 1.
  s0 = 10.
  c = 10.
  τ_plus = 2.
  τ_minus = 0.5
  k = 5.
  ftol = 0.01
  xtol = 0.01
  ctol = 0.01
  max_penalty_iteration = 5
  max_convex_iteration = 5
  max_trust_iteration = 5

  SCPParam_TrajOpt(μ0, s0, c, τ_plus, τ_minus, k, ftol, xtol, ctol, max_penalty_iteration, max_convex_iteration, max_trust_iteration)
end

######
# CVX
######
function cost_true(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3}) where T
  U,N = traj.U, SCPP.N
  Jm = 0
  # Trapezoidal
  for k in 2:N
    Jm += sum(1/2*dtp*(U[j,k-1]^2 + U[j,k]^2) for j = 1:u_dim)
  end
  return Jm
end

function cost_true_convexified(traj, traj_prev::Trajectory, OAP::A) where A <: OptAlgorithmProblem{Astrobee3D{T}, AstrobeeSE3, E} where {T,E}
  cost_true(traj, traj_prev, SCPP)
end

#############################
# Trajectory Initializations
#############################
function init_traj_nothing(TOP::TrajectoryOptimizationProblem{Astrobee3D{T}, AstrobeeSE3, E}) where {T,E}
  model, x_init, goal_set = TOP.PD.model, TOP.PD.x_init, TOP.PD.goal_set
  x_dim, u_dim, N, tf_guess = model.x_dim, model.u_dim, TOP.N, TOP.tf_guess

  goal_final = get_first_goal_at_time(goal_set, tf_guess)
  x_goal = goal_final.params.point

  X = repmat(0.5(x_init + x_goal),1,N)
  U = zeros(u_dim,N)
  Trajectory(X, U, tf_guess)
end

function init_traj_straightline(TOP::TrajectoryOptimizationProblem{Astrobee3D{T}, AstrobeeSE3, E}) where {T,E}
  model, x_init, goal_set = TOP.PD.model, TOP.PD.x_init, TOP.PD.goal_set
  x_dim, u_dim, N, tf_guess = model.x_dim, model.u_dim, TOP.N, TOP.tf_guess

  # Set last state to the center of the goal defined latest in time
  t_goal_final = last(goal_set.goals)[1]
  x_goal = zeros(x_dim)
  for goal in values(inclusive(goal_set.goals, searchequalrange(goal_set.goals, t_goal_final)))
    x_goal[goal.ind_coordinates] = center(goal.params)
  end

  X = hcat(range(x_init, stop=x_goal, length=N)...)
  U = zeros(u_dim, N)
  Trajectory(X, U, tf_guess)
end

####################
# Constraint-adding
####################
function initialize_model_params!(SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}, traj_prev::Trajectory) where {T,E}
  N, robot, model = SCPP.N, SCPP.PD.robot, SCPP.PD.model
  x_dim, u_dim = model.x_dim, model.u_dim
  Xp, Up = traj_prev.X, traj_prev.U

  model.f, model.A, model.B = [], [], B_dyn(Xp[:,1],robot,model)
  for k = 1:N
    push!(model.f, f_dyn(Xp[:,k],Up[:,k],robot,model))
    push!(model.A, A_dyn(Xp[:,k],robot,model))
  end
end

function update_model_params!(SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}, traj_prev::Trajectory) where {T,E}
  N, robot, model = SCPP.N, SCPP.PD.robot, SCPP.PD.model
  Xp, Up, f, A = traj_prev.X, traj_prev.U, model.f, model.A

  for k = 1:N
    update_f!(f[k], Xp[:,k], Up[:,k], robot, model)
    update_A!(A[k], Xp[:,k], robot, model)
  end
end

macro scp_shortcut_AstrobeeSE3(traj, traj_prev, SCPP)
  quote
    X, U, Tf = $(esc(traj)).X, $(esc(traj)).U, $(esc(traj)).Tf
    Xp, Up, Tfp, dtp = $(esc(traj_prev)).X, $(esc(traj_prev)).U, $(esc(traj_prev)).Tf, $(esc(traj_prev)).dt
    robot, model, WS, x_init, goal_set = $(esc(SCPP)).PD.robot, $(esc(SCPP)).PD.model, $(esc(SCPP)).WS, $(esc(SCPP)).PD.x_init, $(esc(SCPP)).PD.goal_set
    x_dim, u_dim, N, dh = model.x_dim, model.u_dim, $(esc(SCPP)).N, $(esc(SCPP)).dh
    X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh
  end
end

## Dynamics constraints
function dynamics_constraints(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}, k::Int, i::Int) where {T,E}
  # Where i is the state index, and k is the timestep index
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_AstrobeeSE3(traj, traj_prev, SCPP)
  fp, Ap, Bp = get_f(k, model), get_A(k, model), get_B(k, model)

  # p previous
  Xkp, Ukp, Xk = X[:,k-1], U[:,k-1], X[:,k]
  fpkp, Apkp, Bpkp, Xpkp, Upkp = get_f(k-1, model), get_A(k-1, model), get_B(k-1, model), Xp[:,k-1], Up[:,k-1]

  # Just Trapezoidal rule
  Uk = U[:,k]
  fpk, Apk, Bpk, Xpk, Upk = get_f(k, model), get_A(k, model), get_B(k, model), Xp[:,k], Up[:,k]
  return (Xkp-Xk) + 1/2*dtp.*(fpkp + Apkp*(Xkp-Xpkp) + Bpkp*(Ukp-Upkp) +
                              fpk  + Apk *(Xk-Xpk)   + Bpk *(Uk-Upk))
end

# Get current dynamics structures for a time step
get_f(k::Int, model::AstrobeeSE3) = model.f[k]
get_A(k::Int, model::AstrobeeSE3) = model.A[k]
get_B(k::Int, model::AstrobeeSE3) = model.B

# Generate full dynamics structures for a time step
function f_dyn(x::Vector, u::Vector, robot::Robot, model::AstrobeeSE3)
  x_dim = model.x_dim
  f = zeros(x_dim)
  update_f!(f, x, u, robot, model)
  return f
end

function update_f!(f, x::Vector, u::Vector, robot::Robot, model::AstrobeeSE3)
  r, v, p, w = x[1:3], x[4:6], x[7:9], x[10:12]
  F, M = u[1:3], u[4:6]

  f[1:3] = v
  f[4:6] = 1/robot.mass*F

  # SO(3)
  f[7:9] = mrp_derivative(p,w)
  f[10:12] = robot.Jinv*(M - cross(w,robot.J*w))
end

function A_dyn(x::Vector, robot::Robot, model::AstrobeeSE3)
  x_dim = model.x_dim
  A = zeros(x_dim, x_dim)
  A[1:6,1:6] = kron([0 1; 0 0], Eye(3))
  update_A!(A, x, robot, model)
  return A
end

function update_A!(A, x::Vector, robot::Robot, model::AstrobeeSE3)
  Jxx, Jyy, Jzz = diag(robot.J)
  px, py, pz = x[7:9]
  wx, wy, wz = x[10:12]

  # See scripts/symbolic_math.m for derivation of these equations
  A[7,7] = (px*wx)/2+(py*wy)/2+(pz*wz)/2
  A[7,8] = wz/2+(px*wy)/2-(py*wx)/2
  A[7,9] = (px*wz)/2-wy/2-(pz*wx)/2
  A[7,10] = px^2/4-py^2/4-pz^2/4+1/4
  A[7,11] = (px*py)/2-pz/2
  A[7,12] = py/2+(px*pz)/2

  A[8,7] = (py*wx)/2-(px*wy)/2-wz/2
  A[8,8] = (px*wx)/2+(py*wy)/2+(pz*wz)/2
  A[8,9] = wx/2+(py*wz)/2-(pz*wy)/2
  A[8,10] = pz/2+(px*py)/2
  A[8,11] = -px^2/4+py^2/4-pz^2/4+1/4
  A[8,12] = (py*pz)/2-px/2

  A[9,7] = wy/2-(px*wz)/2+(pz*wx)/2
  A[9,8] = (pz*wy)/2-(py*wz)/2-wx/2
  A[9,9] = (px*wx)/2+(py*wy)/2+(pz*wz)/2
  A[9,10] = (px*pz)/2-py/2
  A[9,11] = px/2+(py*pz)/2
  A[9,12] = -px^2/4-py^2/4+pz^2/4+1/4

  A[10,11] =  (Jyy-Jzz)*wz/Jxx
  A[10,12] =  (Jyy-Jzz)*wy/Jxx
  A[11,10] = -(Jxx-Jzz)*wz/Jyy
  A[11,12] = -(Jxx-Jzz)*wx/Jyy
  A[12,10] =  (Jxx-Jyy)*wy/Jzz
  A[12,11] =  (Jxx-Jyy)*wx/Jzz
end

function B_dyn(x::Vector, robot::Robot, model::AstrobeeSE3)
  x_dim, u_dim = model.x_dim, model.u_dim
  B = zeros(x_dim, u_dim)
  B[1:6,1:3] = 1/robot.mass*kron([0;1], Eye(3))
  B[10:12,4:6] = robot.Jinv   # SO(3)
  return B
end

## Convex state inequality constraints
function csi_translational_velocity_bound(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_AstrobeeSE3(traj, traj_prev, SCPP)
  return sum(X[3+j,k]^2 for j = 1:3)  - robot.hard_limit_vel^2
end

function csi_angular_velocity_bound(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_AstrobeeSE3(traj, traj_prev, SCPP)
  return sum(X[9+j,k]^2 for j = 1:3)  - robot.hard_limit_ω^2
end

## Convex control inequality constraints
function cci_translational_accel_bound(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,x_init,x_goal,x_dim,u_dim,N,dh = @scp_shortcut_AstrobeeSE3(traj, traj_prev, SCPP)
  return 1/robot.mass^2*sum(U[j,k]^2 for j = 1:3)  - robot.hard_limit_accel^2
end

function cci_angular_accel_bound(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_AstrobeeSE3(traj, traj_prev, SCPP)
  return sum(sum(robot.Jinv[j,i]*U[3+i,k] for i=1:3)^2 for j = 1:3) - robot.hard_limit_α^2
end

## Nonconvex state inequality constraints
function ncsi_obstacle_avoidance_constraints(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_AstrobeeSE3(traj, traj_prev, SCPP)

  rb_idx, env_idx = 1, i
  env_ = WS.btenvironment_keepout

  clearance = model.clearance

  r = X[1:3,k]
  dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r,env_idx)

  # See Eq. 12b in "Convex optimization for proximity maneuvering of a spacecraft with a robotic manipulator"
  return clearance - dist
end

## Nonconvex state inequality constraints (convexified)
function ncsi_obstacle_avoidance_constraints_convexified(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_AstrobeeSE3(traj, traj_prev, SCPP)

  rb_idx, env_idx = 1, i
  env_ = WS.btenvironment_keepout

  clearance = model.clearance

  r0 = get_workspace_location(traj_prev, SCPP, k)
  dist, xbody, xobs = BulletCollision.distance(env_, rb_idx, r0, env_idx)

  if dist < SCPP.param.obstacle_toggle_distance
    r = get_workspace_location(traj, SCPP, k)

    # See Eq. 12b in "Convex optimization for proximity maneuvering of a spacecraft with a robotic manipulator"
    nhat = dist > 0 ?
      (xbody-xobs)./norm(xbody-xobs) :
      (xobs-xbody)./norm(xobs-xbody)

    return (clearance - (dist + nhat'*(r-r0)))
  else
    return 0.
  end
end

## State trust region inequality constraints
function stri_state_trust_region(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_AstrobeeSE3(traj, traj_prev, SCPP)
  return sum((X[j,k]-Xp[j,k])^2 for j = 1:x_dim)
end

## Trust region inequality constraints
function ctri_control_trust_region(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_AstrobeeSE3(traj, traj_prev, SCPP)
  return sum((U[j,k]-Up[j,k])^2 for j = 1:u_dim)
end

function get_workspace_location(traj, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}, k::Int, i::Int=0) where {T,E}
  return traj.X[1:3,k]
end

# TODO(ambyld): Possibly make a custom version of this for each algorithm
function SCPConstraints(SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}) where {T,E}
  model = SCPP.PD.model
  x_dim, u_dim, N = model.x_dim, model.u_dim, SCPP.N
  WS = SCPP.WS
  goal_set, tf_guess = SCPP.PD.goal_set, SCPP.tf_guess

  SCPC = SCPConstraints()

  ## Dynamics constraints
  add_constraint_category!(SCPC.dynamics, dynamics_constraints, :array, 2:N)

  ## State init constraints
  add_constraint_category!(SCPC.state_init_eq, sie_init_constraints, :scalar, 1, 1:x_dim)

  ## State boundary condition constraints
  # add_constraint_category!(SCPC.convex_state_boundary_condition_eq, csbce_goal_constraints, get_first_goal_at_time(goal_set, tf_guess), :array)
  for goal in values(inclusive(goal_set.goals, searchsortedfirst(goal_set.goals, tf_guess), searchsortedlast(goal_set.goals, tf_guess)))
    if typeof(goal.params) == PointGoal
      add_constraint_category!(SCPC.convex_state_boundary_condition_eq, csbce_goal_constraints, goal, :array)
    else
      add_constraint_category!(SCPC.convex_state_boundary_condition_ineq, csbci_goal_constraints, goal, :array)
    end
  end

  ## Convex state inequality constraints
  add_constraint_category!(SCPC.convex_state_ineq, csi_translational_velocity_bound, :scalar, 1:N)
  add_constraint_category!(SCPC.convex_state_ineq, csi_angular_velocity_bound, :scalar, 1:N)

  ## Nonconvex state equality constraints
  nothing

  ## Nonconvex state inequality constraints
  env_ = WS.btenvironment_keepout
  N_obs = length(WS.btenvironment_keepout.convex_env_components)
  add_constraint_category!(SCPC.nonconvex_state_ineq, ncsi_obstacle_avoidance_signed_distance, :scalar, 1:N, 1:N_obs)

  ## Nonconvex state equality constraints (convexified)
  nothing

  ## Nonconvex state inequality constraints (convexified)
  add_constraint_category!(SCPC.nonconvex_state_convexified_ineq, ncsi_obstacle_avoidance_signed_distance_convexified, :scalar, 1:N, 1:N_obs)

  ## Convex control equality constraints
  nothing

  ## Convex control inequality constraints
  add_constraint_category!(SCPC.convex_control_ineq, cci_translational_accel_bound, :scalar, 1:N-1)
  add_constraint_category!(SCPC.convex_control_ineq, cci_angular_accel_bound, :scalar, 1:N-1)

  ## State trust region ineqality constraints
  add_constraint_category!(SCPC.state_trust_region_ineq, stri_state_trust_region, :scalar, 1:N)

  ## Constrol trust region inequality constraints
  add_constraint_category!(SCPC.control_trust_region_ineq, ctri_control_trust_region, :scalar, 1:N-1)

  return SCPC
end

# TODO(ambyld): Generalize this
function trust_region_ratio_gusto(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}) where {T,E}
  # Where i is the state index, and k is the timestep index
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_AstrobeeSE3(traj, traj_prev, SCPP)
  fp, Ap = model.f, model.A
  num,den = 0, 0
  env_ = WS.btenvironment_keepout

  for k in 1:N-1
    linearized = fp[k] + Ap[k]*(X[:,k]-Xp[:,k])
    num += norm(f_dyn(X[:,k],U[:,k],robot,model) - linearized)
    den += norm(linearized)
  end

  clearance = model.clearance

  for k in 1:N
    r0 = get_workspace_location(traj_prev, SCPP, k)
    r = get_workspace_location(traj, SCPP, k)
    for (rb_idx,body_point) in enumerate(env_.convex_robot_components)
      for (env_idx,convex_env_component) in enumerate(env_.convex_env_components)
        dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r0,env_idx)
        nhat = dist > 0 ?
          (xbody-xobs)./norm(xbody-xobs) :
          (xobs-xbody)./norm(xobs-xbody)
        linearized = clearance - (dist + nhat'*(r-r0))

        dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r,env_idx)

        num += abs((clearance-dist) - linearized)
        den += abs(linearized)
      end
    end
  end
  return num/den
end

function trust_region_ratio_trajopt(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}) where {T,E}
  # Where i is the state index, and k is the timestep index
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_AstrobeeSE3(traj, traj_prev, SCPP)
  fp = model.f
  num, den = 0, 0
  env_ = WS.btenvironment_keepout
  dt = traj.dt

  for k in 1:N-1
    ϕ_old = norm(fp[k] - (Xp[:,k]-Xp[:,k])/dtp, 1)
    ϕ_new = norm(f_dyn(X[:,k],U[:,k],robot,model) - (X[:,k]-X[:,k])/dt, 1)
    ϕ_hat_new = norm(dynamics_constraints(traj,traj_prev,SCPP,k,0), 1)
    num += (ϕ_old-ϕ_new)
    den += (ϕ_old-ϕ_hat_new)
  end

  clearance = model.clearance

  for k in 1:N
    r0 = get_workspace_location(traj_prev, SCPP, k)
    r = get_workspace_location(traj, SCPP, k)
    for (rb_idx,body_point) in enumerate(env_.convex_robot_components)
      for (env_idx,convex_env_component) in enumerate(env_.convex_env_components)
        dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r0,env_idx)
        ϕ_old = clearance-dist

        dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r,env_idx)
        nhat = dist > 0 ?
          (xbody-xobs)./norm(xbody-xobs) :
          (xobs-xbody)./norm(xobs-xbody)
        ϕ_new = clearance-dist
        ϕ_hat_new = clearance - (dist + nhat'*(r-r0))

        num += (ϕ_old-ϕ_new)
        den += (ϕ_old-ϕ_hat_new)
      end
    end
  end

  return num/den
end

function trust_region_ratio_mao(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}) where {T,E}
  # Where i is the state index, and k is the timestep index
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_AstrobeeSE3(traj, traj_prev, SCPP)
  num,den = 0, 0

  cost_true_prev = cost_true(traj_prev, traj_prev, SCPP)
  cost_true_new = cost_true(traj, traj, SCPP)
  for k in 1:N-1
    cost_true_prev += norm(dynamics_constraints(traj_prev, traj_prev, SCPP, k, 0),Inf)
    cost_true_new += norm(dynamics_constraints(traj, traj, SCPP, k, 0),Inf)
  end
  cost_linearized_new = cost_true(traj, traj, SCPP)
  return (cost_true_prev-cost_true_new)/(cost_true_prev-cost_linearized_new)
end

##################
# Shooting Method
##################
function get_dual_cvx(prob::Convex.Problem, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}, solver) where {T,E}
	if solver == "Mosek"
		return -MathProgBase.getdual(prob.model)[1:SCPP.PD.model.x_dim]
	else
		return []
	end
end

function get_dual_jump(SCPS::SCPSolution, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}) where {T,E}
	@show -MathProgBase.getconstrduals(SCPS.solver_model.internalModel)[1:SCPP.PD.model.x_dim]
end

###########
# Dynamics
###########
function interpolate_traj(traj::Trajectory, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}, dt_min=0.1) where {T,E}
  N, robot, model = SCPP.N, SCPP.PD.robot, SCPP.PD.model
  x_dim, u_dim = model.x_dim, model.u_dim
  X, U, dt = traj.X, traj.U, traj.dt

  Nstep = ceil(Int, dt/dt_min)
  Nfull = Nstep*(N-1) + 1
  dtfull = dt/Nstep

  Ufull = Matrix(u_dim, Nfull-1)
  Xfull = Matrix(x_dim, Nfull)
  for k in 1:N-1
    istart = Nstep*(k-1)+1
    Ufull[:, istart:Nstep*k] = repmat(U[:,k], 1, Nstep)

    Xfull[:,istart] = X[:,k]
    Ustep = Ufull[:,istart]
    for i = istart:istart+(Nstep-1)
      k1 = f_dyn(Xfull[:,i],Ustep,robot,model)
      x2 = Xfull[:,i] + 0.5*dtfull*k1
      k2 = f_dyn(x2,Ustep,robot,model)
      x3 = Xfull[:,i] + 0.5*dtfull*k2
      k3 = f_dyn(x3,Ustep,robot,model)
      x4 = Xfull[:,i] + dtfull*k3
      k4 = f_dyn(x4,Ustep,robot,model)
      Xfull[:,i+1] = Xfull[:,i] + 1/6*dtfull*(k1 + 2*k2 + 2*k3 + k4)
    end
  end
  Xfull[:,end] = X[:,end]

  return Trajectory(Xfull, Ufull, Tf)
end

function dynamics_constraint_satisfaction(traj::Trajectory, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}) where {T,E}
  N, robot, model = SCPP.N, SCPP.PD.robot, SCPP.PD.model
  WS = SCPP.WS
  X, U = traj.X, traj.U
  Tf, dt = traj.Tf, traj.dt

  J = 0
  for k in 1:N-1
    J += norm((X[:,k+1]-X[:,k])/dt - f_dyn(X[:,k],U[:,k],robot,model),1)
  end
  return J
end

function verify_collision_free(traj::Trajectory, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3, E}) where {T,E}
  model, WS = SCPP.PD.model, SCPP.WS
  N = SCPP.N

  rb_idx = 1
  env_ = WS.btenvironment_keepout

  clearance = model.clearance

  for env_idx = 1:length(env_.convex_env_components)
    for k = 1:N
      r = get_workspace_location(traj, SCPP, k)
      dist, xbody, xobs = BulletCollision.distance(env_, rb_idx, r, env_idx)
      if dist < 0
        return false, k, dist
      end
    end
  end

  return true, 0, 0.
end

max_stars_repo_name: StanfordASL/GuSTO

max_issues_repo_name: StanfordASL/GuSTO

max_forks_repo_name: StanfordASL/GuSTO


text:
### A Pluto.jl notebook ###
# v0.12.4

using Markdown
using InteractiveUtils

# ╔═╡ ad01571d-3bc3-4d11-851d-5a00e477d6df
md"""
# Condicionales

En Julia, la sintaxis

```julia
if *condición 1*
    *opción 1*
elseif *condición 2*
    *opción 2*
else
    *opción 3*
end
```

Nos permite eventualmente evaluar una de nuestras opciones.

Por ejemplo, tal vez queremos implementar la prueba de FizzBuzz: Dado un número N, imprime "Fizz" si N es divisible entre 3, "Buzz" si N es divisible entre 5, y "FizzBuzz" si N es divisible entre ambos 3 y 5. En cualquier otro caso, imprimo el número mismo.
"""

# ╔═╡ 7dc5bb25-f3f3-4fdf-91f7-c6c5ce8c3934
N = 3

# ╔═╡ ec5720b1-de5c-4eba-aa37-9c2248eb9b01
if (N % 3 == 0) & (N % 5 == 0)
    println("FizzBuzz")
elseif N % 3 == 0
    println("Fizz")
elseif N % 5 == 0
    println("Buzz")
else
    println(N)
end

# ╔═╡ c0e38468-a793-4946-ab22-28b302c50c42
md"""
Ahora digamos que queremos regresar el mayor número de ambos. Escoge tus propios x y y
"""

# ╔═╡ 7f187a5d-561c-4010-9c8e-82fd27de0fed
x = 2

# ╔═╡ d4b84b42-193e-11eb-3a56-b5c792a3bf2e
y = 3

# ╔═╡ f3b473cf-fdbb-43ad-b5b6-54aa039f1ea0
if x > y
    x
else
    y
end

# ╔═╡ a1e393f1-fd92-43f0-b8f5-4f5a760fd690
md"""
Para el último bloque, podemos usar el operador ternario, con la sintaxis

```julia
a ? b : c
```

que equivale a

```julia
if a
    b
else
    c
end
```
"""

# ╔═╡ 9f672275-efe8-4474-b8a2-34502e3301e7
(x > y) ? x : y

# ╔═╡ 903a6dea-57ef-4141-9a93-830cb0f4cc2f
md"""
Un truco relacionado es la evaluación de corto-circuito

```julia
a && b
```
"""

# ╔═╡ 4d2656be-7402-413e-b178-145eb7d73209
(x > y) && println(x)

# ╔═╡ 768514c7-a170-4a9a-8c09-e614d86825d6
(x < y) && println(y)

# ╔═╡ ab52e37a-d650-45a9-a6a7-857d047533df
md"""
Cuando escribimos `a && b`, `b` se ejecuta sólo si `a` se evalúa a `true`.

Si `a` se evalúa a `false`, la expresión `a && b` regresa `false`
"""

# ╔═╡ e4a7bf0b-8ce3-4848-a657-382c71a196f4
md"""
### Ejercicios

5.1 Reescribe FizzBuzz sin usar `elseif`.
"""

# ╔═╡ 9555b539-88b9-434d-b60b-9593dfb76571

# ╔═╡ 27737899-7776-438b-b616-ee7bf2ed901d
md"""
5.2 Reescribe FizzBuzz usando el operador ternario.
"""

# ╔═╡ 0699b9c0-88f3-4bcd-a16c-cc1313a68a96

# ╔═╡ Cell order:
# ╟─ad01571d-3bc3-4d11-851d-5a00e477d6df
# ╠═7dc5bb25-f3f3-4fdf-91f7-c6c5ce8c3934
# ╠═ec5720b1-de5c-4eba-aa37-9c2248eb9b01
# ╟─c0e38468-a793-4946-ab22-28b302c50c42
# ╠═7f187a5d-561c-4010-9c8e-82fd27de0fed
# ╠═d4b84b42-193e-11eb-3a56-b5c792a3bf2e
# ╠═f3b473cf-fdbb-43ad-b5b6-54aa039f1ea0
# ╟─a1e393f1-fd92-43f0-b8f5-4f5a760fd690
# ╠═9f672275-efe8-4474-b8a2-34502e3301e7
# ╟─903a6dea-57ef-4141-9a93-830cb0f4cc2f
# ╠═4d2656be-7402-413e-b178-145eb7d73209
# ╠═768514c7-a170-4a9a-8c09-e614d86825d6
# ╟─ab52e37a-d650-45a9-a6a7-857d047533df
# ╟─e4a7bf0b-8ce3-4848-a657-382c71a196f4
# ╠═9555b539-88b9-434d-b60b-9593dfb76571
# ╟─27737899-7776-438b-b616-ee7bf2ed901d
# ╠═0699b9c0-88f3-4bcd-a16c-cc1313a68a96

max_stars_repo_name: miguelraz/JuliaTutorials

max_issues_repo_name: miguelraz/JuliaTutorials

max_forks_repo_name: miguelraz/JuliaTutorials