text:
import GeometryTypes

@enum DIR north=1 west=2 south=3 east=4
@enum POS northWest=1 northEast=2 southWest=3 southEast=4

type QuadTreeElement
  "Index to parent element"
  parent::Nullable{Int}
  "Index to child element"
  northWest::Nullable{Int}
  "Index to child element"
  northEast::Nullable{Int}
  "Index to child element"
  southWest::Nullable{Int}
  "Index to child element"
  southEast::Nullable{Int}
  "Subdivision level (one-based)"
  level::Int
  "Own Index"
  index::Int

  bbLeftBottomIndex::vertex_index
  bbRightBottomIndex::vertex_index
  bbLeftTopIndex::vertex_index
  bbRightTopIndex::vertex_index
end

type QuadTree{T}

  elements::Array{QuadTreeElement, 1}
  values::Array{Nullable{T}, 1}
  vertices::Array{Point, 1}

  function QuadTree(bb::GeometryTypes.SimpleRectangle{Float64})
    @assert GeometryTypes.width(bb) == GeometryTypes.height(bb)
    @assert GeometryTypes.width(bb) > 0
    @assert GeometryTypes.height(bb) > 0
    bl = Point(GeometryTypes.origin(bb))
    br = bl + GeometryTypes.width(bb) * GeometryTypes.Vec(1.0,0.0)
    tr = br + GeometryTypes.height(bb) * GeometryTypes.Vec(0.0,1.0)
    tl = bl + GeometryTypes.height(bb) * GeometryTypes.Vec(0.0,1.0)
    initialVertices = Array{Point, 1}()
    push!(initialVertices, bl)
    push!(initialVertices, br)
    push!(initialVertices, tl)
    push!(initialVertices, tr)
    root = QuadTreeElement(
      Nullable{Int}(),
      Nullable{Int}(),
      Nullable{Int}(),
      Nullable{Int}(),
      Nullable{Int}(),
      1,
      1,
      1, 2, 3, 4
    )
    new([root], [Nullable{T}()], initialVertices)
  end
end

function has_child(qt::QuadTree, elIndex::ElIndex)
  return !isnull(qt.elements[elIndex].northWest)
end

function push_new_element!{T}(qt::QuadTree{T}, bbbl::vertex_index, bbbr::vertex_index, bbtl::vertex_index, bbtr::vertex_index, newLevel::Int, newIndex::Int, elIndex::Int)
  nwEl = QuadTreeElement(
    Nullable{Int}(elIndex),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Int}(),
    newLevel,
    newIndex,
    bbbl, bbbr, bbtl, bbtr
  )
  push!(qt.elements, nwEl)
  push!(qt.values, Nullable{T}())
end

function get_children(qt::QuadTree, elIndex::Int)
  qtEl = qt.elements[elIndex]
  if isnull(qtEl.northWest)
    children = Nullable{Array{Int, 1}}()
  else
    children = Nullable{Array{Int, 1}}(
                [get(qtEl.northWest),
                get(qtEl.northEast),
                get(qtEl.southWest),
                get(qtEl.southEast)])
  end
  return children
end

"""
    find_neighbour(qt::QuadTree, elIndex::ElIndex, dir::DIR)

Computes the index of the neigbour in direction `dir`.
"""
function find_neighbour(qt::QuadTree, eli::ElIndex, dir::DIR)
  if isnull(qt.elements[eli].parent)
    return Nullable{ElIndex}()
  end

  pos::POS = get_leave_dir(qt, eli)
  elp = qt.elements[eli].parent

  if dir == north
    if pos == southWest
      return (qt.elements[get(elp)].northWest)
    end
    if pos == southEast
      return (qt.elements[get(elp)].northEast)
    end

    neighbour = find_neighbour(qt, get(elp), north)
    if isnull(neighbour)
      return Nullable{ElIndex}()
    end

    if !has_child(qt, get(neighbour))
      return neighbour
    end
    if pos == northEast
      return qt.elements[get(neighbour)].southEast
    end
    if pos == northWest
      return qt.elements[get(neighbour)].southWest
    end
  end

  if dir == south
    if pos == northWest
      return (qt.elements[get(elp)].southWest)
    end
    if pos == northEast
      return (qt.elements[get(elp)].southEast)
    end

    neighbour = find_neighbour(qt, get(elp), south)
    if isnull(neighbour)
      return Nullable{ElIndex}()
    end

    if !has_child(qt, get(neighbour))
      return neighbour
    end
    if pos == southEast
      return qt.elements[get(neighbour)].northEast
    end
    if pos == southWest
      return qt.elements[get(neighbour)].northWest
    end
  end

  if dir == west
    if pos == northEast
      return (qt.elements[get(elp)].northWest)
    end
    if pos == southEast
      return (qt.elements[get(elp)].southWest)
    end

    neighbour = find_neighbour(qt, get(elp), west)
    if isnull(neighbour)
      return Nullable{ElIndex}()
    end

    if !has_child(qt, get(neighbour))
      return neighbour
    end
    if pos == northWest
      return qt.elements[get(neighbour)].northEast
    end
    if pos == southWest
      return qt.elements[get(neighbour)].southEast
    end
  end

  if dir == east
    if pos == northWest
      return (qt.elements[get(elp)].northEast)
    end
    if pos == southWest
      return (qt.elements[get(elp)].southEast)
    end

    neighbour = find_neighbour(qt, get(elp), east)
    if isnull(neighbour)
      return Nullable{ElIndex}()
    end

    if !has_child(qt, get(neighbour))
      return neighbour
    end
    if pos == northEast
      return qt.elements[get(neighbour)].northWest
    end
    if pos == southEast
      return qt.elements[get(neighbour)].southWest
    end
  end
end

"""
   get_half_vertex(qt::QuadTree, elIndex::Int, neighbour::ElIndex, dir::DIR)

Returns the vertex index of the half-grid point of element `elIndex` in
direction `dir`.
"""
function get_half_vertex(qt::QuadTree, elIndex::Int, neighbour::ElIndex, dir::DIR)
  neighbourEl = qt.elements[neighbour]
  if dir == north
    seChild = qt.elements[get(neighbourEl.southEast)]
    swChild = qt.elements[get(neighbourEl.southWest)]
    assert(seChild.bbLeftBottomIndex == swChild.bbRightBottomIndex)
    return seChild.bbLeftBottomIndex
  elseif dir == south
    neChild = qt.elements[get(neighbourEl.northEast)]
    nwChild = qt.elements[get(neighbourEl.northWest)]
    assert(nwChild.bbRightTopIndex == neChild.bbLeftTopIndex)
    return nwChild.bbRightTopIndex
  elseif dir == west
    neChild = qt.elements[get(neighbourEl.northEast)]
    seChild = qt.elements[get(neighbourEl.southEast)]
    assert(neChild.bbRightBottomIndex == seChild.bbRightTopIndex)
    return neChild.bbRightBottomIndex
  elseif dir == east
    nwChild = qt.elements[get(neighbourEl.northWest)]
    swChild = qt.elements[get(neighbourEl.southWest)]
    assert(nwChild.bbLeftBottomIndex == swChild.bbLeftTopIndex)
    return nwChild.bbLeftBottomIndex
  end
end

"""
    get_new_center(qt::QuadTree, elIndex::Int, dir::DIR)

If the boundary in direction `dir` is currently undivided, it adds a new
vertex at its center and returns its index. If the boundary is divided and
that vertex hence exists, it only returns its index.
"""
function get_new_center(qt::QuadTree, elIndex::Int, dir::DIR)

  # get vertex indices of requested boundary line
  if dir == north
    v1, v2 = qt.elements[elIndex].bbLeftTopIndex, qt.elements[elIndex].bbRightTopIndex
  elseif dir == south
    v1, v2 = qt.elements[elIndex].bbLeftBottomIndex, qt.elements[elIndex].bbRightBottomIndex
  elseif dir == west
    v1, v2 = qt.elements[elIndex].bbLeftBottomIndex, qt.elements[elIndex].bbLeftTopIndex
  elseif  dir == east
    v1, v2 = qt.elements[elIndex].bbRightBottomIndex, qt.elements[elIndex].bbRightTopIndex
  end

  # check if neighbour is there and if it is divided
  neighbour = find_neighbour(qt, elIndex, dir)
  if (isnull(neighbour) || !has_child(qt, get(neighbour)))
    newCenter = 0.5*(qt.vertices[v1] + qt.vertices[v2])
    push!(qt.vertices, newCenter)
    return length(qt.vertices)
  end

  # need to return center index
  return get_half_vertex(qt, elIndex, get(neighbour), dir)
end

"""
    subdivide!(qt::QuadTree, elIndex::Int, childrenCreated)

Subdivides quadtree element denoted by `elIndex`.

# Remarks
* If necessary, its neighbours are subdivided as well in order to balance the
  tree.
* The element needs to be a leave, i.e. it cannot have any children.
* `childrenCreated(qt::QuadTree, el::QuadTreeElement)` will be called for each cell that was subdivided
"""
function subdivide!(qt::QuadTree, elIndex::Int, childrenCreated)
  qtEl = qt.elements[elIndex]
  @assert(isnull(qtEl.northWest))
  @assert(isnull(qtEl.northEast))
  @assert(isnull(qtEl.southWest))
  @assert(isnull(qtEl.southEast))

  # we need a new center point
  bl = qt.vertices[qtEl.bbLeftBottomIndex]
  br = qt.vertices[qtEl.bbRightBottomIndex]
  tl = qt.vertices[qtEl.bbLeftTopIndex]
  tr = qt.vertices[qtEl.bbRightTopIndex]
  newCenter = 0.25*(bl+br+tl+tr)
  push!(qt.vertices, newCenter)
  newCenterIndex = length(qt.vertices)

  # get boundary center points
  northCenterIndex = get_new_center(qt, elIndex, north)
  southCenterIndex = get_new_center(qt, elIndex, south)
  westCenterIndex = get_new_center(qt, elIndex, west)
  eastCenterIndex = get_new_center(qt, elIndex, east)

  newLevel = qtEl.level + 1
  newIndex = length(qt.elements) + 1

  # northWest
  push_new_element!(qt, westCenterIndex, newCenterIndex, qtEl.bbLeftTopIndex, northCenterIndex, newLevel, newIndex, elIndex)
  qtEl.northWest = Nullable{Int}(newIndex)
  newIndex = newIndex + 1

  # northEast
  push_new_element!(qt, newCenterIndex, eastCenterIndex, northCenterIndex, qtEl.bbRightTopIndex, newLevel, newIndex, elIndex)
  qtEl.northEast = Nullable{Int}(newIndex)
  newIndex = newIndex + 1

  # southWest
  push_new_element!(qt, qtEl.bbLeftBottomIndex, southCenterIndex, westCenterIndex, newCenterIndex, newLevel, newIndex, elIndex)
  qtEl.southWest = Nullable{Int}(newIndex)
  newIndex = newIndex + 1

  # southEast
  push_new_element!(qt, southCenterIndex, qtEl.bbRightBottomIndex, newCenterIndex, eastCenterIndex, newLevel, newIndex, elIndex)
  qtEl.southEast = Nullable{Int}(newIndex)
  newIndex = newIndex + 1

  childrenCreated(qt, elIndex)

  # look at this parent's neighbors and check if they need subdiving
  # we only admit a difference in 1 to neighbour's subdivision level
  # that means we need to subdivide parent's neighbour if it has
  # no children
  par = qtEl.parent
  if (!isnull(par))
    parpar = (qt.elements[get(par)]).parent
    if (!isnull(parpar))
      children = get_children(qt, get(parpar))
      if !isnull(children)
        for child in get(children)
          if (child != elIndex) && (!has_child(qt, child))
            subdivide!(qt, child, childrenCreated)
          end
        end
      end
    end
  end
end

function get_element_bounding_box(qt::QuadTree, qtEl::QuadTreeElement)
  bl = qt.vertices[qtEl.bbLeftBottomIndex]
  tr = qt.vertices[qtEl.bbRightTopIndex]
  width = (tr - bl)[1]
  height = (tr - bl)[2]
  bb = GeometryTypes.SimpleRectangle(bl[1], bl[2], width, height)

  return bb
end

function get_element_bounding_box(qt::QuadTree, elIndex::ElIndex)
  qtEl = qt.elements[elIndex]

  return get_element_bounding_box(qt, qtEl)
end

function get_child_positions_from_coordinate(qt::QuadTree, qtEl::ElIndex, coord::Point)
  bb = get_element_bounding_box(qt, qtEl)
  cx, cy = coord
  xbl, ybl, w, h = bb.x, bb.y, bb.x + bb.w, bb.y + bb.h
  @assert(cx >= bb.x && cx <= bb.x + bb.w && cy >= bb.y && cy <= bb.y + bb.h)
  tol = 1e-5
  retPositions = Array{POS, 1}()
  #print_with_color(:blue, "Comparing ($cx, $cy) with ($(bb.x + 0.5 * bb.w), $(bb.y + 0.5 * bb.h))")
  if cy > bb.y + 0.5 * bb.h - tol
    if cx > bb.x + 0.5 * bb.w - tol
      push!(retPositions, northWest)
    end
    if cx < bb.x + 0.5 * bb.w + tol
      push!(retPositions, northEast)
    end
  end
  if cy < bb.y + 0.5 * bb.h + tol
    if cx < bb.x + 0.5 * bb.w + tol
      push!(retPositions, southWest)
    end
    if cx > bb.x + 0.5 * bb.w - tol
      push!(retPositions, southEast)
    end
  end

  return retPositions
end

function query(qt::QuadTree, testFunction, curElIndex::ElIndex, els::Array{ElIndex, 1})
  if !testFunction(get_element_bounding_box(qt, curElIndex))
    return els
  end

  children = get_children(qt, curElIndex)
  if (!isnull(children))
    ch = get(children)
    nw = query(qt, testFunction, ch[1], Array{ElIndex, 1}())
    ne = query(qt, testFunction, ch[2], Array{ElIndex, 1}())
    sw = query(qt, testFunction, ch[3], Array{ElIndex, 1}())
    se = query(qt, testFunction, ch[4], Array{ElIndex, 1}())

    return vcat(els, nw, ne, sw, se)
  end

  return vcat(els, [curElIndex])
end

function query_line(qt::QuadTree, ls::GeometryTypes.LineSegment{Point})
  # taken from
  # http://stackoverflow.com/questions/99353/how-to-test-if-a-line-segment-intersects-an-axis-aligned-rectange-in-2d
  function line_intersects_rectangle(r::GeometryTypes.SimpleRectangle{Float64})
    ((x1, y1), (x2, y2)) = ls
    xbl, ybl, xtr, ytr = r.x, r.y, r.x + r.w, r.y + r.h
    f(x, y) = (y2 - y1) * x + (x1 - x2) * y + (x2 * y1 - x1 * y2)
    s1 = f(xbl, ybl)
    s2 = f(xbl, ytr)
    s3 = f(xtr, ybl)
    s4 = f(xtr, ytr)
    if (s1>0 && s2>0 && s3>0 && s4>0) || (s1<0 && s2<0 && s3<0 && s4<0)
      return false
    end
    if (x1 > xtr && x2 > xtr)
      return false
    end
    if (x1 < xbl && x2 < xbl)
      return false
    end
    if (y1 > ytr && y2 > ytr)
      return false
    end
    if (y1 < ybl && y2 < ybl)
      return false
    end
    return true
  end

  return query(qt, line_intersects_rectangle, 1, Array{ElIndex, 1}())
end

function refine_line_to_level(
  qt::QuadTree,
  ls::GeometryTypes.LineSegment{Point},
  targetLevel::Int,
  childrenCreated)
  function filtFunc(eli::ElIndex)
    throwout = !has_child(qt, eli) && qt.elements[eli].level < targetLevel
    return !has_child(qt, eli) && qt.elements[eli].level < targetLevel
  end

  leaves = query_line(qt, ls)
  filter!(filtFunc, leaves)
  while(!isempty(leaves))
    for li in leaves
      if (!has_child(qt, li) && qt.elements[li].level < targetLevel)
        subdivide!(qt, li, childrenCreated)
      end
    end
    leaves = query_line(qt, ls)
    filter!(filtFunc, leaves)
  end
end

"""
    get_element_dir(qt::QuadTree, el::ElIndex)

Gets the relative position of the leave wrt parent.

#Remarks
# If the element is the root element, it returns `northWest`
"""
function get_leave_dir(qt::QuadTree, eli::ElIndex)
  el = qt.elements[eli]
  if (isnull(el.parent))
    return northWest
  end
  elp = qt.elements[get(el.parent)]

  if (!isnull(elp.northWest) && get(elp.northWest) == eli)
    return northWest
  end
  if (!isnull(elp.northEast) && get(elp.northEast) == eli)
    return northEast
  end
  if (!isnull(elp.southWest) && get(elp.southWest) == eli)
    return southWest
  end
  if (!isnull(elp.southEast) && get(elp.southEast) == eli)
    return southEast
  end
  error("Quadtree corruption: cannot locate element $eli in parent node.")
end

function get_leave_from_pos(qt::QuadTree, eli::ElIndex, pos::POS)
  @assert(has_child(qt, eli))
  el = qt.elements[eli]

  if pos == northWest
    return get(el.northWest)
  elseif pos == northEast
    return get(el.northEast)
  elseif pos == southWest
    return get(el.southWest)
  elseif pos == southEast
    return get(el.southEast)
  end
end


text:
function getvalue(A,line,col)
	A[sum(1:line-1)+col]
end

function getmax()
end

function getmax(A,N, line,col)
	if line==N
		return getvalue(A,line,col)
	else
		# left
		a = getmax(A,N, line+1, col)
		b = getmax(A,N, line+1, col+1)
		return getvalue(A,line,col) + max(a,b)
	end
end

function main()
	A = [75 95 64 17 47 82 18 35 87 10 20 04 82 47 65 19 01 23 75 03 34 88 02 77 73 07 63 67 99 65 04 28 06 16 70 92 41 41 26 56 83 40 80 70 33 41 48 72 33 47 32 37 16 94 29 53 71 44 65 25 43 91 52 97 51 14 70 11 33 28 77 73 17 78 39 68 17 57 91 71 52 38 17 14 91 43 58 50 27 29 48 63 66 04 68 89 53 67 30 73 16 69 87 40 31 04 62 98 27 23 09 70 98 73 93 38 53 60 04 23]
	N = 15 #lines
	assert( sum(1:N) == length(A))

	println(getmax(A,N,1,1))
end

main()

ext: jl

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
#==============================================================================#
# AWSCredentials.jl
#
# Load AWS Credentials from:
# - EC2 Instance Profile,
# - Environment variables, or
# - ~/.aws/credentials file.
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

export AWSCredentials,
       localhost_is_lambda,
       localhost_is_ec2,
       localhost_maybe_ec2,
       aws_user_arn,
       aws_account_number,
       check_credentials

"""
When you interact with AWS, you specify your [AWS Security Credentials](http://docs.aws.amazon.com/general/latest/gr/aws-security-credentials.html)
to verify who you are and whether you have permission to access the resources that you are requesting.
AWS uses the security credentials to authenticate and authorize your requests.

The fields `access_key_id` and `secret_key` hold the access keys used to authenticate API requests
(see [Creating, Modifying, and Viewing Access Keys](http://docs.aws.amazon.com/IAM/latest/UserGuide/id_credentials_access-keys.html#Using_CreateAccessKey)).

[Temporary Security Credentials](http://docs.aws.amazon.com/IAM/latest/UserGuide/id_credentials_temp.html) require the extra session `token` field.

The `user_arn` and `account_number` fields are used to cache the result of the [`aws_user_arn`](@ref) and [`aws_account_number`](@ref) functions.

AWSCore searches for credentials in a series of possible locations and stop as soon as it finds credentials.
The order of precedence for this search is as follows:

1. Passing credentials directly to the `AWSCredentials` constructor
2. [Environment variables](https://docs.aws.amazon.com/cli/latest/userguide/cli-configure-envvars.html)
3. Shared credential file [(~/.aws/credentials)](http://docs.aws.amazon.com/cli/latest/userguide/cli-config-files.html)
4. AWS config file [(~/.aws/config)](http://docs.aws.amazon.com/cli/latest/userguide/cli-config-files.html)
5. Assume Role provider via the aws config file
6. Instance metadata service on an Amazon EC2 instance that has an IAM role configured.

Once the credentials are found, the method by which they were accessed is stored in the `renew` field
and the DateTime at which they will expire is stored in the `expiry` field.
This allows the credentials to be refreshed as needed using [`check_credentials`](@ref).
If `renew` is set to `nothing`, no attempt will be made to refresh the credentials.

Any renewal function is expected to return `nothing` on failure or a populated `AWSCredentials` object on success.
The `renew` field of the returned `AWSCredentials` will be discarded and does not need to be set.

To specify the profile to use from `~/.aws/credentials`, do, for example, `AWSCredentials(profile="profile-name")`.
"""
mutable struct AWSCredentials
    access_key_id::String
    secret_key::String
    token::String
    user_arn::String
    account_number::String
    expiry::DateTime
    renew::Union{Function, Nothing}

    function AWSCredentials(access_key_id,secret_key,
                            token="", user_arn="", account_number="";
                            expiry=typemax(DateTime),
                            renew=nothing)
        new(access_key_id, secret_key, token, user_arn, account_number, expiry, renew)
    end
end

function AWSCredentials(;profile=nothing)
    creds = nothing
    renew = Nothing

    # Define our search options
    functions = [
        env_instance_credentials,
        () -> dot_aws_credentials(profile),
        () -> dot_aws_config(profile),
        instance_credentials,
    ]

    # Loop through our search locations until we get credentials back
    for f in functions
        renew = f
        creds = renew()
        creds === nothing || break
    end

    creds === nothing && error("Can't find AWS credentials!")
    creds.renew = renew

    if debug_level > 0
        display(creds)
        println()
    end

    return creds
end

will_expire(cr::AWSCredentials) = now(UTC) >= cr.expiry - Minute(5)

"""
    check_credentials(cr::AWSCredentials; force_refresh::Bool=false)

Checks current AWSCredentials, refreshing them if they are soon to expire.
If force_refresh is `true` the credentials will be renewed immediately.
"""
function check_credentials(cr::AWSCredentials; force_refresh::Bool=false)
    if force_refresh || will_expire(cr)
        if debug_level > 0
            println("Renewing credentials... ")
        end
        renew = cr.renew

        if renew !== nothing
            new_creds = renew()

            new_creds === nothing && error("Can't find AWS credentials!")
            copyto!(cr, new_creds)

            # Ensure renewal function is not overwritten by the new credentials
            cr.renew = renew
        else
            if debug_level > 0
                println("Credentials cannot be renewed...")
            end
        end
    end

    return cr
end

function Base.show(io::IO,c::AWSCredentials)
    println(io, string(c.user_arn,
                       c.user_arn == "" ? "" : " ",
                       "(",
                       c.account_number,
                       c.account_number == "" ? "" : ", ",
                       c.access_key_id,
                       c.secret_key == "" ? "" : ", $(c.secret_key[1:3])...",
                       c.token == "" ? "" : ", $(c.token[1:3])..."),
                       ")")
end

function Base.copyto!(dest::AWSCredentials, src::AWSCredentials)
    for f in fieldnames(typeof(dest))
        setfield!(dest, f, getfield(src, f))
    end
end
import Base: copy!
Base.@deprecate copy!(dest::AWSCredentials, src::AWSCredentials) copyto!(dest, src)

"""
Is Julia running in an AWS Lambda sandbox?
"""
localhost_is_lambda() = haskey(ENV, "LAMBDA_TASK_ROOT")

"""
Is Julia running on an EC2 virtual machine?
"""
function localhost_is_ec2()

    if localhost_is_lambda()
        return false
    end

    if isfile("/sys/hypervisor/uuid") &&
       String(read("/sys/hypervisor/uuid",3)) == "ec2"
        return true
    end

    if isfile("/sys/devices/virtual/dmi/id/product_uuid")
        try
            # product_uuid is not world readable!
            # https://patchwork.kernel.org/patch/6461521/
            # https://github.com/JuliaCloud/AWSCore.jl/issues/24
            if String(read("/sys/devices/virtual/dmi/id/product_uuid")) == "EC2"
                return true
            end
        catch
        end
    end

    return false
end

localhost_maybe_ec2() = localhost_is_ec2() ||
                        isfile("/sys/devices/virtual/dmi/id/product_uuid")

"""
    aws_user_arn(::AWSConfig)

Unique
[Amazon Resource Name]
(http://docs.aws.amazon.com/IAM/latest/UserGuide/id_users.html)
for configrued user.

e.g. `"arn:aws:iam::account-ID-without-hyphens:user/Bob"`
"""
function aws_user_arn(aws::AWSConfig)
    creds = aws[:creds]
    if creds.user_arn == ""
        r = Services.sts(aws, "GetCallerIdentity", [])
        creds.user_arn = r["Arn"]
        creds.account_number = r["Account"]
    end
    return creds.user_arn
end

"""
    aws_account_number(::AWSConfig)

12-digit [AWS Account Number](http://docs.aws.amazon.com/general/latest/gr/acct-identifiers.html).
"""
function aws_account_number(aws::AWSConfig)
    creds = aws[:creds]
    if creds.account_number == ""
        aws_user_arn(aws)
    end
    return creds.account_number
end

"""
    ec2_metadata(key)

Fetch [EC2 meta-data]
(http://docs.aws.amazon.com/AWSEC2/latest/UserGuide/ec2-instance-metadata.html)
for `key`.
"""
function ec2_metadata(key)

    @assert localhost_maybe_ec2()

    String(http_get("http://169.254.169.254/latest/meta-data/$key").body)
end

function instance_credentials()
    if haskey(ENV, "AWS_CONTAINER_CREDENTIALS_RELATIVE_URI")
        return ecs_instance_credentials()
    elseif localhost_maybe_ec2()
        return ec2_instance_credentials()
    else
        return nothing
    end
end

"""
Load [Instance Profile Credentials]
(http://docs.aws.amazon.com/AWSEC2/latest/UserGuide/iam-roles-for-amazon-ec2.html#instance-metadata-security-credentials)
for EC2 virtual machine.
"""
function ec2_instance_credentials()

    @assert localhost_maybe_ec2()

    info  = ec2_metadata("iam/info")
    info  = LazyJSON.value(info)

    name  = ec2_metadata("iam/security-credentials/")
    creds = ec2_metadata("iam/security-credentials/$name")
    new_creds = LazyJSON.value(creds)

    if debug_level > 0
        print("Loading AWSCredentials from EC2 metadata... ")
    end

    expiry = DateTime(strip(new_creds["Expiration"], 'Z'))

    AWSCredentials(new_creds["AccessKeyId"],
                   new_creds["SecretAccessKey"],
                   new_creds["Token"],
                   info["InstanceProfileArn"];
                   expiry = expiry)
end

"""
Load [ECS Task Credentials]
(http://docs.aws.amazon.com/AmazonECS/latest/developerguide/task-iam-roles.html)
"""
function ecs_instance_credentials()

    @assert haskey(ENV, "AWS_CONTAINER_CREDENTIALS_RELATIVE_URI")

    uri = ENV["AWS_CONTAINER_CREDENTIALS_RELATIVE_URI"]

    new_creds = String(http_get("http://169.254.170.2$uri").body)
    new_creds = LazyJSON.value(new_creds)

    if debug_level > 0
        print("Loading AWSCredentials from ECS metadata... ")
    end

    expiry = DateTime(strip(new_creds["Expiration"], 'Z'))

    AWSCredentials(new_creds["AccessKeyId"],
                   new_creds["SecretAccessKey"],
                   new_creds["Token"],
                   new_creds["RoleArn"];
                   expiry = expiry)
end

"""
Load Credentials from [environment variables]
(http://docs.aws.amazon.com/cli/latest/userguide/cli-environment.html)
`AWS_ACCESS_KEY_ID`, `AWS_SECRET_ACCESS_KEY` etc.
(e.g. in Lambda sandbox).
"""
function env_instance_credentials()

    if haskey(ENV, "AWS_ACCESS_KEY_ID")
        if debug_level > 0
            print("Loading AWSCredentials from ENV[\"AWS_ACCESS_KEY_ID\"]... ")
        end

        return AWSCredentials(
            ENV["AWS_ACCESS_KEY_ID"],
            ENV["AWS_SECRET_ACCESS_KEY"],
            get(ENV, "AWS_SESSION_TOKEN", ""),
            get(ENV, "AWS_USER_ARN", "");
            renew = env_instance_credentials
        )
    else
        return nothing
    end
end

using IniFile

function dot_aws_credentials_file()
    get(ENV, "AWS_SHARED_CREDENTIALS_FILE", joinpath(homedir(), ".aws", "credentials"))
end

"""
Try to load Credentials from [AWS CLI ~/.aws/credentials file]
(http://docs.aws.amazon.com/cli/latest/userguide/cli-config-files.html)
"""
function dot_aws_credentials(profile = nothing)
    creds = nothing
    credential_file = dot_aws_credentials_file()

    ini = nothing
    if isfile(credential_file)
        ini = read(Inifile(), credential_file)
        key, key_id, token = aws_get_credential_details(
            profile === nothing ? aws_get_profile() : profile,
            ini,
            false
        )

        if key !== :notfound
            creds = AWSCredentials(key_id, key, token)
        end
    end

    return creds
end

dot_aws_config_file() = get(ENV, "AWS_CONFIG_FILE", joinpath(homedir(), ".aws", "config"))

"""
Try to load Credentials or assume a role via the [AWS CLI ~/.aws/config file]
(http://docs.aws.amazon.com/cli/latest/userguide/cli-config-files.html)
"""
function dot_aws_config(profile = nothing)
    creds = nothing
    config_file = dot_aws_config_file()

    ini = nothing
    if isfile(config_file)
        ini = read(Inifile(), config_file)
        p = profile === nothing ? aws_get_profile() : profile
        key, key_id, token = aws_get_credential_details(p, ini, true)

        if key !== :notfound
            creds = AWSCredentials(key_id, key, token)
        else
            creds = aws_get_role(p, ini)
        end
    end

    return creds
end

function aws_get_role_details(profile::AbstractString, ini::Inifile)
    if debug_level > 0
        println("Loading \"$profile\" Profile from " *
                dot_aws_config_file() * "... ")
    end

    role_arn = get(ini, profile, "role_arn")
    source_profile = get(ini, profile, "source_profile")

    profile = "profile $profile"
    role_arn = get(ini, profile, "role_arn", role_arn)
    source_profile = get(ini, profile, "source_profile", source_profile)

    (source_profile, role_arn)
end

function aws_get_credential_details(profile::AbstractString, ini::Inifile, config::Bool)
    if debug_level > 0
        filename = config ? dot_aws_config_file() : dot_aws_credentials_file()
        println("Loading \"$profile\" AWSCredentials from " * filename
                * "... ")
    end

    key_id = get(ini, profile, "aws_access_key_id")
    key = get(ini, profile, "aws_secret_access_key")
    token = get(ini, profile, "aws_session_token", "")

    if config
        profile = "profile $profile"
        key_id = get(ini, profile, "aws_access_key_id", key_id)
        key = get(ini, profile, "aws_secret_access_key", key)
        token = get(ini, profile, "aws_session_token", token)
    end

    (key, key_id, token)
end

function aws_get_profile()
    get(ENV, "AWS_DEFAULT_PROFILE", get(ENV, "AWS_PROFILE", "default"))
end

function aws_get_region(profile::AbstractString, ini::Inifile)
    region = get(ENV, "AWS_DEFAULT_REGION", "us-east-1")

    region = get(ini, profile, "region", region)
    region = get(ini, "profile $profile", "region", region)
end

function aws_get_role(role::AbstractString, ini::Inifile)
    source_profile, role_arn = aws_get_role_details(role, ini)
    source_profile === :notfound && return nothing

    if debug_level > 0
        println("Assuming \"$source_profile\"... ")
    end
    credentials = nothing

    for f in [dot_aws_credentials, dot_aws_config]
        credentials = f(source_profile)
        credentials === nothing || break
    end

    credentials === nothing && return nothing

    config = AWSConfig(:creds=>credentials, :region=>aws_get_region(source_profile, ini))

    role = Services.sts(
        config,
        "AssumeRole",
        RoleArn=role_arn,
        RoleSessionName=replace(role, r"[^\w+=,.@-]" => s"-"),
    )
    role_creds = role["Credentials"]

    AWSCredentials(role_creds["AccessKeyId"],
        role_creds["SecretAccessKey"],
        role_creds["SessionToken"];
        expiry = unix2datetime(role_creds["Expiration"]))
end

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


text:
using Test

using SimpleDirectMediaLayer
SDL2 = SimpleDirectMediaLayer

include("../src/objects.jl")
include("../assets/configs.jl")

# Objects tests
kZeroVel = Vector2D(0,0)

@testset "isColliding(paddle, line)" begin
l = LineSegment(WorldPos(-10,1), WorldPos(-10,-1))
p = Paddle(WorldPos(0,0), kZeroVel, 20)
@test isColliding(p, l, 0)
@test isColliding(p, l, 10)
end

@testset "willCollide(ball, paddle, dt)" begin
@testset "dt" begin
b = Ball(WorldPos(-10,-1), Vector2D(0,2))
p = Paddle(WorldPos(0,0), kZeroVel, 30)
@test !willCollide(b, p, 0)
@test willCollide(b, p, 1)
@test willCollide(b, p, .5)
@test !willCollide(b, p, .4)
p = Paddle(WorldPos(0,0), kZeroVel, 5)  # too short
@test !willCollide(b, p, 1)

b = Ball(WorldPos(0,0), Vector2D(0,-2))
p = Paddle(WorldPos(0,-100), kZeroVel, 10)
@test !willCollide(b, p, 1)
end

@testset "collision detection" begin
b = Ball(WorldPos(-10,9), Vector2D(0,2))
p = Paddle(WorldPos(-8,10), kZeroVel, 5)
@test willCollide(b, p, 1)
b = Ball(WorldPos(-10,11), Vector2D(0,2))
@test !willCollide(b, p, 1)

b = Ball(WorldPos(0,-195), Vector2D(0,-200))
p = Paddle(WorldPos(0,-200), kZeroVel, 5)
@test willCollide(b, p, 1)
b = Ball(WorldPos(0,-200), Vector2D(0,-200000))
@test willCollide(b, p, 1)
@test willCollide(b, p, 0)
b = Ball(WorldPos(0,0), Vector2D(0,-200))
@test willCollide(b, p, 1)
b = Ball(WorldPos(0,0), Vector2D(0,-200))
@test willCollide(b, p, 1)
end

end

max_stars_repo_name: NHDaly/pongClone

max_issues_repo_name: NHDaly/pongClone

max_forks_repo_name: NHDaly/pongClone


text:
module WaveletMatrices

export WaveletMatrix, getindex, rank, select, freq

import Base: lastindex, length, size, sizeof

using IndexableBitVectors
import IndexableBitVectors: rank

include("build.jl")

"""
Indexable sequence.

`WaveletMatrix` is a data structure like the wavelet tree, which supports fast rank/select
queries for a sequence of alphabets. Elements (typed `T`) are encoded in `w` bits and therefore
a sequence of unsigned integers can be stored in a space-efficient manner if `w` is small.
The actual memory space is determined according to the underlying bit vector type.
The default bit vector type is `SucVector`, which requires 5/4 bits per bit, so the total
size will be about `w * 5/4 * length` bits.

See (Claude et al, 2012, doi:10.1007/978-3-642-34109-0_18) for more details.
"""
struct WaveletMatrix{w,T<:Unsigned,B<:AbstractBitVector} <: AbstractVector{T}
    bits::NTuple{w,B}
    nzeros::NTuple{w,Int}
    sps::Vector{Int}
    function WaveletMatrix{w,T,B}(bits) where {w,T<:Unsigned,B<:AbstractBitVector}
        @assert 1 ≤ w ≤ sizeof(T) * 8 ≤ 64
        @assert length(bits) == w
        nzeros = Int[]
        for bv in bits
            @assert length(bits[1]) == length(bv)
            push!(nzeros, rank0(bv, length(bv)))
        end
        if w ≤ 16
            # size of lookup table ≤ 512KiB (= sizeof(Int) * 2^16)
            alphabetsize = 2^w
            sps = Vector{Int}(undef, alphabetsize)
            for a in 0:alphabetsize-1
                sps[a+1] = locate_sp(T(a), bits, nzeros)
            end
        else
            sps = Int[]
        end
        return new(tuple(bits...), tuple(nzeros...), sps)
    end
    function WaveletMatrix{w}(data::AbstractVector{T}; destructive::Bool=false) where {w, T<:Unsigned}
       bits = build(default_bitvector, data, w, destructive)
       return WaveletMatrix{w,T,default_bitvector}(bits)
    end
    function WaveletMatrix(data::AbstractVector{T}, w::Integer=sizeof(T) * 8) where {T <: Unsigned}
       return WaveletMatrix{w}(data)
    end
end

const default_bitvector = SucVector

length(wm::WaveletMatrix) = length(wm.bits[1])
size(wm::WaveletMatrix) = (length(wm),)

function sizeof(wm::WaveletMatrix{w}) where w
    s = 0
    for d in 1:w
        s += sizeof(wm.bits[d])
    end
    s += sizeof(wm.nzeros)
    s += sizeof(wm.sps)
    return s
end

@inline function Base.getindex(wm::WaveletMatrix{w,T}, i::Int) where {w,T}
    if i < 0 || lastindex(wm) < i
        throw(BoundsError(i))
    end
    ret = T(0)
    @inbounds begin
        for d in 1:w-1
            bits = wm.bits[d]
            bit = bits[i]
            ret = ret << 1 | bit
            if bit
                i = wm.nzeros[d] + rank1(bits, i)
            else
                i = rank0(bits, i)
            end
        end
        bits = wm.bits[w]
        bit = bits[i]
        ret = ret << 1 | bit
    end
    return ret
end

@inline function Base.getindex(wm::WaveletMatrix{w,T,SucVector}, i::Int) where {w,T}
    if i < 0 || lastindex(wm) < i
        throw(BoundsError(i))
    end
    ret = T(0)
    @inbounds begin
        for d in 1:w-1
            bits = wm.bits[d]
            bit, rnk1 = IndexableBitVectors.accrank1(bits, i)
            ret = ret << 1 | bit
            if bit
                i = wm.nzeros[d] + rnk1
            else
                i = i - rnk1
            end
        end
        bits = wm.bits[w]
        bit = bits[i]
        ret = ret << 1 | bit
    end
    return ret
end

@inline Base.getindex(wm::WaveletMatrix{w,T}, i::Integer) where {w,T} = getindex(wm, convert(Int, i))

function rank(a::Unsigned, wm::WaveletMatrix{w}, i::Int) where w
    i = clamp(i, 0, lastindex(wm))
    if i == 0
        return 0
    end

    if !isempty(wm.sps)
        # use precomputed sp
        sp = wm.sps[a+1]
        ep = i
        @inbounds for d in 1:w
            bit = (a >> (w - d)) & 1 == 1
            ep = rank(bit, wm.bits[d], ep)
            if bit
                ep += wm.nzeros[d]
            end
        end
    else
        sp, ep = 0, i
        # scan from the most significant bit to the least significant bit
        @inbounds for d in 1:w
            bit = (a >> (w - d)) & 1 == 1
            sp = rank(bit, wm.bits[d], sp)
            ep = rank(bit, wm.bits[d], ep)
            if bit
                nz = wm.nzeros[d]
                sp += nz
                ep += nz
            end
        end
    end
    return ep - sp
end

rank(a::Unsigned, wm::WaveletMatrix, i::Integer) = rank(a, wm, convert(Int, i))

function select(a::Unsigned, wm::WaveletMatrix, j::Int)
    if j ≤ 0
        return 0
    end
    # binary search: j ∈ (rank(l), rank(u)]
    l = 0
    u = length(wm)
    rank_u = rank(a, wm, u)
    while j ≤ rank_u
        m = div(l + u, 2)
        if l == m
            return u
        end
        rank_m = rank(a, wm, m)
        if j ≤ rank_m
            u = m
            rank_u = rank_m
        else
            l = m
        end
    end
    return 0
end

select(a::Unsigned, wm::WaveletMatrix, j::Integer) = select(a, wm, convert(Int, j))

function freq(a::Unsigned, wm::WaveletMatrix, i::Integer, j::Integer)
    return j < i ? 0 : rank(a, wm, j) - rank(a, wm, i - 1)
end

end # module


text:
function plot_int_found(LN, lnglobal, z1, z2 ; marker = "b.", sleeptime = 0.01)

	figure(1) ; clf()
	LNPLOT=[]
	for i = 1:length(LN)
		push!(LNPLOT, LN.yn[i])
		i!=length(LN) && push!(LNPLOT, nadirs(LN)[i])
	end
	plot(map(x->x[1], LNPLOT), map(x->x[2], LNPLOT), "r--")
	plot(map(x->x[1], LN.yn), map(x->x[2], LN.yn), "r^")
	plot(map(x->x[1], nadirs(LN)), map(x->x[2], nadirs(LN)), "rs", markersize="4")
	plot([z1], [z2], marker)

	λ = LN.λ
	WS = sum((z1,z2).*λ)
	plot([LN.yn[1][1], (WS-λ[2]*LN.yn[end][2])/λ[1]], [(WS-λ[1]*LN.yn[1][1])/λ[2], LN.yn[end][2]], "k--")

	ax = gca()
    # ax[:set_xlim]([-1+min(LN.yn[1][1], z1), 1+max(LN.yn[end][1], z1+1)])
	# ax[:set_ylim]([-1+min(LN.yn[end][2]), 1+max(LN.yn[1][2], z2+1)])
	ax[:set_xlim]([LN.yn[1][1] - 1, LN.yn[end][1]+1])
	ax[:set_ylim]([LN.yn[end][2]-1, LN.yn[1][2] + 1])
	show()

	# figure(2)
	# LNPLOT=[]
	# for i = 1:length(lnglobal)
	# 	push!(LNPLOT, lnglobal.yn[i])
	# 	i!=length(lnglobal) && push!(LNPLOT, nadirs(lnglobal)[i])
	# end
	# plot(map(x->x[1], LNPLOT), map(x->x[2], LNPLOT), "r--")
	# plot(map(x->x[1], lnglobal.yn), map(x->x[2], lnglobal.yn), "r^", markersize="2")
	# plot(map(x->x[1], nadirs(lnglobal)), map(x->x[2], nadirs(lnglobal)), "rs", markersize="2")
	# # plot([0, z1+(lnglobal.λ[2]/lnglobal.λ[1])*z2 ], [z2+(lnglobal.λ[1]/lnglobal.λ[2])*z1, 0], "b--")
	# plot([z1], [z2], marker, markersize="2")
	# show()

	sleep(sleeptime)
end

function plot_pareto_branch(LN, z1, z2, bound1, bound2 ; sleeptime = 1)

	figure(1)
	clf()

	LNPLOT=[]
	for i = 1:length(LN)
		push!(LNPLOT, LN.yn[i])
		i!=length(LN) && push!(LNPLOT, nadirs(LN)[i])
	end
	plot(map(x->x[1], LNPLOT), map(x->x[2], LNPLOT), "r--")
	plot(map(x->x[1], LN.yn), map(x->x[2], LN.yn), "r^")
	plot(map(x->x[1], nadirs(LN)), map(x->x[2], nadirs(LN)), "rs", markersize="4")
	# plot([0, z1+(LN.λ[2]/LN.λ[1])*z2 ], [z2+(LN.λ[1]/LN.λ[2])*z1, 0], "b--")
	!isnan(bound1) && plot([bound1, bound1], [LN.yn[1][2], LN.yn[end][2]])
	!isnan(bound2) && plot([LN.yn[1][1], LN.yn[end][1]], [bound2, bound2])
	plot([z1], [z2], "bo")

	λ = LN.λ
	WS = sum((z1,z2).*λ)
	plot([LN.yn[1][1], (WS-λ[2]*LN.yn[end][2])/λ[1]], [(WS-λ[1]*LN.yn[1][1])/λ[2], LN.yn[end][2]], "k--")

	ax = gca()
    ax[:set_xlim]([-1+min(LN.yn[1][1], z1), 1+max(LN.yn[end][1], z1)])
	ax[:set_ylim]([-1+min(LN.yn[end][2]), 1+max(LN.yn[1][2])])
	show()

	sleep(sleeptime)
end

function plotdualbound(sl::Vector{Segment{T}}) where T
	for s in sl
		if T == Min
			plot([s.p1.x, s.p2.x, s.c.x, s.c.x, s.p1.x, s.p1.x], [s.p1.y, s.p2.y, s.p2.y, s.c.y, s.c.y, s.p1.y], "g-", linewidth=1)
		else
			plot([s.p1.x, s.p2.x, s.p2.x, s.c.x, s.c.x, s.p1.x], [s.p1.y, s.p2.y, s.c.y, s.c.y, s.p1.y, s.p1.y], "g-", linewidth=1)
		end
		plot([s.p1.x, s.p2.x], [s.p1.y, s.p2.y], "gs", markersize=3)
		plot([s.c.x, s.c.x], [s.c.y, s.c.y], "kx", markersize=3)
	end
end


text:
# Pointers with address space information

# Address spaces

export AS

module AS

const Generic  = 0
const Global   = 1
const Region   = 2
const Local    = 3
const Constant = 4
const Private  = 5

end # module


text:
abstract type IntegratorCache{DT,D} end

abstract type ODEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type DAEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type IODEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type IDAEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type PODEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type PDAEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end

create_integrator_cache(integrator::Integrator) = error("create_integrator_cache()! not implemented for ", typeof(integrator))

initialize!(::Integrator, ::AtomicSolution) = nothing

integrate_step!(integrator::Integrator, ::AtomicSolution) = error("integrate_step()! not implemented for ", typeof(integrator))


text:
"""
    curl(A::AbstractArray{T,N}, dx=ones(T,3)) where {T,N}

Calculate 2nd order cell-centered ∇×A where `A` is a 4D array of size (3, nx, ny, nz) and
`dx` is a vector of grid intervals in each dimension.
"""
function curl(A::AbstractArray{T,N}, dx=ones(T,3)) where {T,N}
   @assert N == 4 && length(dx) == 3 "Input vector shall be indexed in 3D!"

   @views Ax, Ay, Az = A[1,:,:,:], A[2,:,:,:], A[3,:,:,:]

   B = zeros(T, size(A))
   @views Bx, By, Bz = B[1,:,:,:], B[2,:,:,:], B[3,:,:,:]

   invdx = @. inv(2*dx)

   if any(==(1), size(A)) # 2D
      if size(A,4) == 1
         @inbounds for j in 2:size(A,3)-1, i in 2:size(A,2)-1
            ∂Ay∂x = (-Ay[i-1,j,1] + Ay[i+1,j,1]) * invdx[1]
            ∂Az∂x = (-Az[i-1,j,1] + Az[i+1,j,1]) * invdx[1]
            ∂Ax∂y = (-Ax[i-1,j,1] + Ax[i+1,j,1]) * invdx[2]
            ∂Az∂y = (-Az[i-1,j,1] + Az[i+1,j,1]) * invdx[2]
            ∂Ax∂z = (-Ax[i-1,j,1] + Ax[i+1,j,1]) * invdx[3]
            ∂Ay∂z = (-Ay[i-1,j,1] + Ay[i+1,j,1]) * invdx[3]

            Bx[i,j,1] = ∂Az∂y - ∂Ay∂z
            By[i,j,1] = ∂Ax∂z - ∂Az∂x
            Bz[i,j,1] = ∂Ay∂x - ∂Ax∂y
         end
      elseif size(A,3) == 1
         @inbounds for k in 2:size(A,4)-1, i in 2:size(A,2)-1
            ∂Ay∂x = (-Ay[i-1,1,k] + Ay[i+1,1,k]) * invdx[1]
            ∂Az∂x = (-Az[i-1,1,k] + Az[i+1,1,k]) * invdx[1]
            ∂Ax∂y = (-Ax[i-1,1,k] + Ax[i+1,1,k]) * invdx[2]
            ∂Az∂y = (-Az[i-1,1,k] + Az[i+1,1,k]) * invdx[2]
            ∂Ax∂z = (-Ax[i-1,1,k] + Ax[i+1,1,k]) * invdx[3]
            ∂Ay∂z = (-Ay[i-1,1,k] + Ay[i+1,1,k]) * invdx[3]

            Bx[i,1,k] = ∂Az∂y - ∂Ay∂z
            By[i,1,k] = ∂Ax∂z - ∂Az∂x
            Bz[i,1,k] = ∂Ay∂x - ∂Ax∂y
         end
      end
   else # 3D
      @inbounds for k in 2:size(A,4)-1, j in 2:size(A,3)-1, i in 2:size(A,2)-1
         ∂Ay∂x = (-Ay[i-1,j,  k  ] + Ay[i+1,j  ,k  ]) * invdx[1]
         ∂Az∂x = (-Az[i-1,j,  k  ] + Az[i+1,j  ,k  ]) * invdx[1]
         ∂Ax∂y = (-Ax[i  ,j-1,k  ] + Ax[i  ,j+1,k  ]) * invdx[2]
         ∂Az∂y = (-Az[i  ,j-1,k  ] + Az[i  ,j+1,k  ]) * invdx[2]
         ∂Ax∂z = (-Ax[i  ,j,  k-1] + Ax[i  ,j  ,k+1]) * invdx[3]
         ∂Ay∂z = (-Ay[i  ,j,  k-1] + Ay[i  ,j  ,k+1]) * invdx[3]

         Bx[i,j,k] = ∂Az∂y - ∂Ay∂z
         By[i,j,k] = ∂Ax∂z - ∂Az∂x
         Bz[i,j,k] = ∂Ay∂x - ∂Ax∂y
      end
   end
   B
end

"""
    gradient(A::AbstractArray{T,N}, dx=ones(T, 3)) where {T,N}

Calculate 2nd order cell-centered ∇A where `A` is a scalar array and `dx` is a vector of
grid intervals in each dimension.
!!! warning
    The current implementation has issues at the boundary if gradient is taken multiple
    times.
"""
function gradient(A::AbstractArray{T,N}, dx=ones(T, N)) where {T,N}
   @assert N < 4 "$N dimension array A detected!"
   @assert N == length(dx) "The array A shall have the same dimension as dx!"
   @assert all(!=(1), size(A)) "No singular dimension is allowed!"

   if N == 3 # 3D
      B = zeros(T, 3, size(A)[1], size(A)[2], size(A)[3])
      invdx = @. inv(2*dx)
      @views Bx, By, Bz = B[1,:,:,:], B[2,:,:,:], B[3,:,:,:]

      @inbounds for k in 2:size(A,3)-1, j in 2:size(A,2)-1, i in 2:size(A,1)-1
         ∂A∂x = (-A[i-1,j  ,k  ] + A[i+1,j,  k  ]) * invdx[1]
         ∂A∂y = (-A[i  ,j-1,k  ] + A[i,  j+1,k  ]) * invdx[2]
         ∂A∂z = (-A[i  ,j  ,k-1] + A[i,  j,  k+1]) * invdx[3]

         Bx[i,j,k] = ∂A∂x
         By[i,j,k] = ∂A∂y
         Bz[i,j,k] = ∂A∂z
      end
   elseif N == 2 # 2D
      B = zeros(T, 2, size(A)[1], size(A)[2])
      invdx = @. inv(2*dx)
      @views Bx, By = B[1,:,:], B[2,:,:]

      @inbounds for j in 2:size(A,2)-1, i in 2:size(A,1)-1
         ∂A∂x = (-A[i-1,j  ] + A[i+1,j  ]) * invdx[1]
         ∂A∂y = (-A[i  ,j-1] + A[i  ,j+1]) * invdx[2]

         Bx[i,j] = ∂A∂x
         By[i,j] = ∂A∂y
      end
   else # 1D
      B = zeros(T, size(A)[1])
      invdx = @. inv(2*dx)
      @views Bx = B[:]

      @inbounds for i in 2:size(A,1)-1
         ∂A∂x = (-A[i-1] + A[i+1]) * invdx[1]

         Bx[i] = ∂A∂x
      end
   end
   B
end

"""
    divergence(A::AbstractArray{T,N}, dx=ones(T, 3)) where {T,N}

Calculate 2nd order cell-centered ∇⋅A where `A` is a 4D array of size (3, nx, ny, nz) and
`dx` is a vector of grid intervals in each dimension.
"""
function divergence(A::AbstractArray{T,N}, dx=ones(T, 3)) where {T,N}
   @assert N == 4 && length(dx) == 3 "Input vector shall be indexed in 3D!"
   @assert all(!=(1), size(A)) "Input vector must be from 3D data!"

   @views Ax, Ay, Az = A[1,:,:,:], A[2,:,:,:], A[3,:,:,:]

   B = zeros(T, size(A)[2:end])

   invdx = @. inv(2*dx)

   @inbounds for k in 2:size(A,4)-1, j in 2:size(A,3)-1, i in 2:size(A,2)-1
      ∂Ax∂x = (-Ax[i-1,j  ,k  ] + Ax[i+1,j  ,k  ]) * invdx[1]
      ∂Ay∂y = (-Ay[i  ,j-1,k  ] + Ay[i  ,j+1,k  ]) * invdx[2]
      ∂Az∂z = (-Az[i  ,j  ,k-1] + Az[i  ,j  ,k+1]) * invdx[3]

      B[i,j,k] = ∂Ax∂x + ∂Ay∂y + ∂Az∂z
   end
   B
end


text:
vandercorput(num::Integer, base::Integer) = sum(d * Float64(base) ^ -ex for (ex,d) in enumerate(digits(num, base=base)))

function quasirandom(x, y;order=:none, squish=0.9)

	if order == :none
		grps = unique(x)
	else
		grps = order
	end

	nGroups = length(grps)
	width= squish/nGroups

	k = Array{UnivariateKDE}(undef, nGroups)
	maxDensity = zeros(nGroups)

	xTemp = zeros(length(x))
	for i in 1:nGroups
		xTemp[findall(elements -> elements == grps[i], x)] .= i
	end

	q = zeros(length(x))
	d = zeros(length(x))
	for i in 1:nGroups
		k[i] = kde(y[x .== grps[i]])
		d[x .== grps[i]] .= pdf(k[i], y[x .== grps[i]])
		maxDensity[i] = maximum(d[x .== grps[i]])
		q[x .== grps[i]] .= vandercorput.(1:length(q)
			, 2)[competerank(y[x .== grps[i]])]
	end

	x_jitter = xTemp .+ width./maxDensity[nGroups] .*  (q .- 0.5) .* 2 .* d

	return x_jitter .- 0.5
end

function quasirandomScatter(x, y; order=:none, squish=0.9, kwargs...)

	plt = scatter(x, y; markercolor=:white, markerstrokewidth=0,
		xlim=(0, length(unique(x))), margin=2mm)
	scatter!(quasirandom(x, y;order=order, squish), y; kwargs...)

	return plt
end


text:
using Turing, Random, CSV, DataFrames

function run_example()
    Random.seed!(12)

    # Reproducing the pymc3 implementation found here:
    # https://www.ritchievink.com/blog/2018/10/09/build-facebooks-prophet-in-pymc3-bayesian-time-series-analyis-with-generalized-additive-models/

    # downloaded from https://raw.githubusercontent.com/facebook/prophet/master/examples/example_wp_log_peyton_manning.csv
    url = "https://raw.githubusercontent.com/facebook/prophet/master/examples/example_wp_log_peyton_manning.csv"
    p = joinpath("data", basename(url))

    if !isfile(p)
        mkpath(dirname(p))
        download(url, p)
    end
    df = CSV.read(p)

    df[!, :yscale] = df[!, :y] ./ maximum(df[!, :y])

    min_ds, max_ds = extrema(df.ds)
    df[!, :t] = (df.ds - min_ds) / (max_ds - min_ds)

    # plot(df.ds, df.yscale)

    nchangepoints = 10

    # Proportion of history in which trend changepoints will be estimated.
    changepoints_range = .8
    # The standard deviation of the prior on the growth.
    growth_prior_scale = 5
    # The scale of the Laplace prior on the delta vector.
    changepoints_prior_scale = .05

    @model trendmodel(t, y, ::Type{TV}=Vector{Float64}) where TV = begin  # <= CHANGED
        # num observations
        nobs = length(t)
        s =  collect(range(1, stop=nobs*changepoints_range, length=nchangepoints))
        A = [(t[i] >= s[j])*1 for i=1:nobs, j=1:nchangepoints]

        k ~ Normal(0.0, growth_prior_scale)

        # rate of change
        delta = TV(undef, nchangepoints)  # <= CHANGED
        delta ~ [Laplace(0.0, changepoints_prior_scale)]

        # offset
        m ~ Normal(0.0, 5.0)

        gamma = -s.*delta

        g = (k .+ A*delta).*t + (m .+ A * gamma)

        # The number of observations.
        sd ~ Truncated(Cauchy(0.0, .5), 0.000001, Inf)
        for n in 1:nobs
            y[n] ~ Normal(g[n], sd)
        end
    end

    model = trendmodel(df.t, df.yscale)

    chain = sample(model, NUTS(1500, 200, 0.65), progress=true);
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM.Testing

#=
Two rings, RING1 = inner, RING2 = outer, RINGS combined mesh. Set T=1.0 for
inner ring and T=2.0 for outer ring, measure temperature from middle of ring.
Results are calculated using Code Aster for comparison.
=#
@testset "test 3d heat, two rings, and compare to CA solution" begin
    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
    mesh = aster_read_mesh(meshfile, "RINGS_UNION")

    rings = Problem(Heat, "RINGS", 1)
#   rings.elements = create_elements(mesh; element_type=:Tet4)
    rings.elements = create_elements(mesh, "RING1", "RING2")
    update!(rings.elements, "temperature thermal conductivity", 1.0)
    bc_inner = Problem(Dirichlet, "INNER SURFACE", 1, "temperature")
    bc_inner.elements = create_elements(mesh, "RING1_INNER")
    bc_outer = Problem(Dirichlet, "OUTER SURFACE", 1, "temperature")
    bc_outer.elements = create_elements(mesh, "RING2_OUTER")
    update!(bc_inner, "temperature 1", 1.0)
    update!(bc_outer, "temperature 1", 2.0)
    info("# of elements in RING1_INNER = ", length(bc_inner.elements))
    info("# of elements in RING2_OUTER = ", length(bc_outer.elements))
    solver = LinearSolver(rings, bc_inner, bc_outer)
    solver()

    temp_jf = rings("temperature", 0.0)

    fn = Pkg.dir("JuliaFEM") * "/test/testdata/rings.rmed"
    results = RMEDFile(fn)
    nodes = aster_read_nodes(results)
    temp_ca = aster_read_data(results, "TEMP")

    passed = true
    for j in sort(collect(keys(temp_jf)))
        X = nodes[j]
        T1 = temp_jf[j]
        T2 = temp_ca[j]
        rtol = norm(T1-T2) / max(T1,T2)
        @printf "% 5i : %8.5f %8.5f %8.5f | %8.5f %8.5f | %8.5e\n" j X... T1 T2 rtol
        passed &= rtol < 1.0e-12
    end
    @test passed
end


text:
using Dates
import Dates.value
import Dates.days

# Define a date-time with microsecond accuracy.
# Note that this doesn't include all features of a `DateTime`.
# Includes code derived from Julia base/dates/*.

struct Microsecond <: Dates.TimePeriod
    value::Int64
end

Dates.slotparse(slot::Dates.Slot{Microsecond}, x, locale) =
    !ismatch(r"[^0-9\s]",x) ? slot.parser(Base.parse(Float64,"."*x)*1e6) : throw(SLOTERROR)
# Dates.slotformat(slot::Slot{Microsecond},dt,locale) = rpad(string(microsecond(dt)/1000.0)[3:end], slot.width, "0")
Dates.SLOT_RULE['z'] = Microsecond
Dates.periodisless(x::Microsecond,y::Microsecond) = value(x) < value(y)
Dates.periodisless(::Microsecond,::Hour)   = true
Dates.periodisless(::Microsecond,::Minute) = true
Dates.periodisless(::Microsecond,::Second) = true
Dates.periodisless(::Microsecond,::Millisecond) = true
Dates.periodisless(::Period,::Microsecond) = false
Dates.coarserperiod(::Type{Microsecond}) = (Millisecond, 1000)
Dates.toms(c::Microsecond) = div(value(c), 1000)
Dates.days(c::Microsecond) = div(value(c), 86_400_000_000)

# DateTime is a millisecond precision UTInstant interpreted by ISOCalendar
immutable DateTimeMicro <: TimeType
    instant::Dates.UTInstant{Microsecond}
    DateTimeMicro(instant::Dates.UTInstant{Microsecond}) = new(instant)
end

microsecond(dt::DateTimeMicro) = mod(value(dt),1000000)

Dates.days(dt::DateTimeMicro) = fld(value(dt),86400000000)
Dates.hour(dt::DateTimeMicro)   = mod(fld(value(dt),3600000000),24)
Dates.minute(dt::DateTimeMicro) = mod(fld(value(dt),60000000),60)
Dates.second(dt::DateTimeMicro) = mod(fld(value(dt),1000000),60)

function Base.string(dt::DateTimeMicro)
    y,m,d = yearmonthday(days(dt))
    h,mi,s = hour(dt),minute(dt),second(dt)
    yy = y < 0 ? @sprintf("%05i",y) : lpad(y,4,"0")
    mm = lpad(m,2,"0")
    dd = lpad(d,2,"0")
    hh = lpad(h,2,"0")
    mii = lpad(mi,2,"0")
    ss = lpad(s,2,"0")
    μs = microsecond(dt) == 0 ? "" : string(microsecond(dt)/1e6)[2:end]
    return "$yy-$mm-$(dd)T$hh:$mii:$ss$(μs)"
end
Base.show(io::IO,x::DateTimeMicro) = print(io,string(x))

DateTimeMicro(x::Int64) = DateTimeMicro(Dates.UTInstant(Microsecond(x)))

_c(x) = convert(Int64,x)
DateTimeMicro(y,m=1,d=1,h=0,mi=0,s=0,μs=0) = DateTimeMicro(_c(y),_c(m),_c(d),_c(h),_c(mi),_c(s),_c(μs))

"""
    DateTimeMicro(y, [m, d, h, mi, s, μs]) -> DateTimeMicro
Construct a `DateTimeMicro` type by parts. Arguments must be convertible to `Int64`.
"""
function DateTimeMicro(y::Int64=1,m::Int64=1,d::Int64=1,
                       h::Int64=0,mi::Int64=0,s::Int64=0,μs::Int64=0)
    0 < m < 13 || throw(ArgumentError("Month: $m out of range (1:12)"))
    0 < d < daysinmonth(y,m)+1 || throw(ArgumentError("Day: $d out of range (1:$(daysinmonth(y,m)))"))
    -1 < h < 24 || throw(ArgumentError("Hour: $h out of range (0:23)"))
    -1 < mi < 60 || throw(ArgumentError("Minute: $mi out of range (0:59)"))
    -1 < s < 60 || throw(ArgumentError("Second: $s out of range (0:59)"))
    -1 < μs < 1000000 || throw(ArgumentError("Millisecond: $μs out of range (0:999999)"))
    rata = μs + 1000000*(s + 60mi + 3600h + 86400 * Dates.totaldays(y,m,d))
    return DateTimeMicro(Dates.UTInstant(Microsecond(rata)))
end
# Convenience constructors from Periods
function DateTimeMicro(y::Year,m::Month=Month(1),d::Day=Day(1),
                       h::Hour=Hour(0),mi::Minute=Minute(0),
                       s::Second=Second(0),μs::Microsecond=Microsecond(0))
    return DateTimeMicro(value(y),value(m),value(d),
                         value(h),value(mi),value(s),value(μs))
end

DateTimeMicro(dt::AbstractString,format::AbstractString;locale::AbstractString="english") = DateTimeMicro(dt,DateFormat(format,locale))
DateTimeMicro(dt::AbstractString,df::DateFormat=DateFormat("m/d/y,H:M:S.z")) = DateTimeMicro(Dates.parse(dt,df)...)

Base.convert(::Type{DateTimeMicro},dt::DateTime) =
    DateTimeMicro(Dates.UTInstant(Microsecond(value(dt) * 1000)))
Base.convert(::Type{DateTime},dt::DateTimeMicro) =
    DateTime(Dates.UTInstant(Millisecond(div(value(dt), 1000))))
Base.promote_rule(::Type{DateTime},x::Type{DateTimeMicro}) = DateTimeMicro
Base.isless(x::DateTimeMicro,y::DateTimeMicro) = isless(value(x),value(y))


text:
include("tree_func.jl")
include("solve_relax.jl")

function branch_and_price(data, graph, approx, verbose)
    """Branch-and-Price algorithm wich can output an approximated solution
    or an otpimal solution."""

    if verbose; println("--------------------------"); end
    if verbose; println("BRANCH-AND-PRICE"); end
    start = now()

    node_pile = []
    branch_pool = []
    tour_pool = []
    global_upper_bound = Nothing
    global_best_node = Nothing

    if verbose; println("-------"); end
    if verbose; println("Creating the root ..."); end

    push!(tour_pool, create_initial_tour(data, graph))
    root = Dict(
        "type" => "root",
        "branch" => [],
        "tours" => [1],
        "lower_bound" => 0,
        "alpha" => []
        )
    global_upper_bound = first(tour_pool)["cost"]
    global_best_node = root

    solve_relax(root, data, graph, global_upper_bound, tour_pool, branch_pool, global_best_node)

    if approx
        if verbose; println("-------"); end
        if verbose; println("Solving master without relaxation ..."); end
        solve_master(root, data, tour_pool, branch_pool, false)
        tours = [tour_pool[j] for (i,j) in enumerate(global_best_node["tours"]) if global_best_node["alpha"][i]>=1-10^-5]
        tours = remove_multiple_service(tours, graph)
        stop = now()
        solving_time = stop-start
        return(tours, solving_time)
    end

    pushfirst!(node_pile, root)

    # Explore tree
    while length(node_pile)!=0

        node = popfirst!(node_pile)
        if verbose; println("-------"); end
        if verbose; println("Exploring ", node["type"], " node"); end

        if !cut_branch(node, global_upper_bound, data, verbose)

            if verbose; println("-------"); end
            (branch_index, branch_pool) = create_branch(node, tour_pool, branch_pool, data)
            (even_node, odd_node) = create_child_node(node, branch_index)
            if verbose; println("Branch on depot ", branch_pool[branch_index][1],
                " and on ", branch_pool[branch_index][2], " number ",
                branch_pool[branch_index][3]); end

            if verbose; println("-------"); end
            if verbose; println("Solving odd node relaxation ..."); end
            (global_upper_bound, tour_pool, branch_pool, global_best_node) = solve_relax(odd_node, data, graph, global_upper_bound, tour_pool, branch_pool, global_best_node)
            pushfirst!(node_pile, odd_node)

            if verbose; println("-------"); end
            if verbose; println("Solving even node relaxation ..."); end
            (global_upper_bound, tour_pool, branch_pool, global_best_node) = solve_relax(even_node, data, graph, global_upper_bound, tour_pool, branch_pool, global_best_node)

            pushfirst!(node_pile, even_node)
        end
    end

    tours = [tour_pool[j] for (i,j) in enumerate(global_best_node["tours"]) if global_best_node["alpha"][i]>=1-10^-5]
    if (global_best_node==root) | (tour_pool[1] in tours)
        error("No solutions for the problem. Try to use more vehicle or
        increase the capacity.")
    end

    stop = now()
    solving_time = stop-start
    if verbose; println("--------------------------"); end
    return(tours, solving_time)
end

max_stars_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_issues_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_forks_repo_name: TheoGuyard/Groomer-Optimization-Problem


text:
using Plots

using simulator.materials.spesifics
using simulator.targetfigures

# Simulation parameters start

save = false

n1 = SiO2_core_Sellmeier
n2 = Au_Johnson
n3 = Pd052_unloaded
n4 = Air
n5 = Air

d1 = 10e-9
d2 = 30e-9
d3 = 5e-9
d4 = 10e-9
d5 = 10e-9
d2smax = 40e-9
d3smax = 40e-9
d2ss = 1e-9
d3ss = 1e-10
d3smrd2sm = 0.4

λ = 1200e-9
λs = [x for x in 600e-9:1e-10:1400e-9]

Hcs = [x for x in 1e-2:1e-2:1e-1]
Pd_c = 0.52

θd = 45

description = "52 unloaded Pd on Au"
fig_path = "../rapport/figures/targetplots/"

# Simulation parameters end

fig_token = replace(description, " " => "_")

if save
    gr()
else
    plotly()
end
θ1 = θd * π / 180

# Plot singe line:
println("\nSingle line timing:")
@time begin
system_single = make_layered_tm_system([n1, n2, n3, n4, n5], [d1, d2, d3, d4, d5], θ1)
xs, ys = scann_singleparameter(system_single, λs)
end
fig_single = plot(xs .* 1e9, ys, title="$description, Single scan spectra", xaxis="λ [nm]", yaxis="Power reflection coefficient", legend=false)
display(fig_single)
if save savefig(fig_single, fig_path * "single_" * fig_token * ".pdf") end

# Plot hydrogen concentration scan:
println("\nScan hydrogen timing:")
@time begin
system_Hc = make_hc_system(n1, n2, n3, n4, n5, d1, d2, d3, d4, d5, θ1, Pd_c)
xs, ysr, ysλ, ysw, ysΔ = scan_plasmon_singleparameter(system_Hc, Hcs, λs, 2)
end
fig_d3r = plot(xs .* 1e2, ysr, title="$description, hydrogen scan", xaxis="Hydrogen concentration [%]", yaxis="Power reflection coefficient minima", legend=false)
fig_d3λ = plot(xs .* 1e2, ysλ .* 1e9, title="$description, hydrogen scan", xaxis="Hydrogen concentration [%]", yaxis="λ [nm]", legend=false)
fig_d3w = plot(xs .* 1e2, ysw .* 1e9, title="$description, hydrogen scan", xaxis="Hydrogen concentration [%]", yaxis="λ plasmon dip FWHM [nm]", legend=false)
fig_d3Δ = plot(xs .* 1e2, ysΔ, title="$description, hydrogen scan", xaxis="Hydrogen concentration [%]", yaxis="Plasmon peak power reflection Δ", legend=false)
display(fig_d3r)
if save savefig(fig_d3r, fig_path * "Hcr_" * fig_token * ".pdf") end
display(fig_d3λ)
if save savefig(fig_d3λ, fig_path * "Hcl_" * fig_token * ".pdf") end
display(fig_d3w)
if save savefig(fig_d3w, fig_path * "Hcw_" * fig_token * ".pdf") end
display(fig_d3Δ)
if save savefig(fig_d3Δ, fig_path * "Hcd_" * fig_token * ".pdf") end

# Plot d3 scan:
println("\nScan d3 timing:")
@time begin
system_d3 = make_d3_system(n1, n2, n3, n4, n5, d1, d2, d4, d5, θ1)
xs, ysr, ysλ, ysw, ysΔ = scan_plasmon_singleparameter(system_d3, [x for x in 0:d3ss:(d2*d3smrd2sm)], λs, 2)
end
fig_d3r = plot(xs .* 1e9, ysr, title="$description, Active layer scan", xaxis="d [nm]", yaxis="Power reflection coefficient minima", legend=false)
fig_d3λ = plot(xs .* 1e9, ysλ .* 1e9, title="$description, Active layer scan", xaxis="d [nm]", yaxis="λ [nm]", legend=false)
fig_d3w = plot(xs .* 1e9, ysw .* 1e9, title="$description, Active layer scan", xaxis="d [nm]", yaxis="λ plasmon dip FWHM [nm]", legend=false)
fig_d3Δ = plot(xs .* 1e9, ysΔ, title="$description, Active layer scan", xaxis="d [nm]", yaxis="Plasmon peak power reflection Δ", legend=false)
display(fig_d3r)
if save savefig(fig_d3r, fig_path * "d3r_" * fig_token * ".pdf") end
display(fig_d3λ)
if save savefig(fig_d3λ, fig_path * "d3l_" * fig_token * ".pdf") end
display(fig_d3w)
if save savefig(fig_d3w, fig_path * "d3w_" * fig_token * ".pdf") end
display(fig_d3Δ)
if save savefig(fig_d3Δ, fig_path * "d3d_" * fig_token * ".pdf") end

# Plot d2 on d3 scan:
println("\nScan d2 and d3 timing:")
@time begin
system_d2 = make_d2_dualparameter_system(n1, n2, n3, n4, n5, d1, d4, d5, θ1)
xs_d2, ys_d3R, _, ys_d3ν, ys_d3Δ = scan_plasmon_dualthicknesses(system_d2, 0, d2smax, d2ss, 0, d3smax, d3ss, λs, 1, d3smrd2sm)
end
fig_d2R = plot(xs_d2 .* 1e9, ys_d3R .* 1e9, title="$description, minimum plasmon dip", xaxis="Au thickness [nm]", yaxis="Pd thickness [nm]", legend=false)
display(fig_d2R)
if save savefig(fig_d2R, fig_path * "d2-d3-R_" * fig_token * ".pdf") end
fig_d2ν = plot(xs_d2 .* 1e9, ys_d3ν .* 1e9, title="$description, minimum plasmon FWHM", xaxis="Au thickness [nm]", yaxis="Pd thickness [nm]", legend=false)
display(fig_d2ν)
if save savefig(fig_d2ν, fig_path * "d2-d3-FWHM_" * fig_token * ".pdf") end
fig_d2R = plot(xs_d2 .* 1e9, ys_d3Δ .* 1e9, title="$description, maximum plasmon peak to peak", xaxis="Au thickness [nm]", yaxis="Pd thickness [nm]", legend=false)
display(fig_d2R)
if save savefig(fig_d2R, fig_path * "d2-d3-P2P_" * fig_token * ".pdf") end


text:
# @everywhere module_dir = "/Users/brandonbehring/Desktop/Leap_Frog_2019"
# @everywhere push!(LOAD_PATH, "/Users/brandonbehring/Desktop/Leap_Frog_2019")
# # @everywhere module_dir ="/Users/brandonbehring/Desktop/Leap_Frog_2019"
# # @everywhere push!(LOAD_PATH, $module_dir)
# # @everywhere thisDir = dirname(@__FILE__())
# # @everywhere any(path -> path == thisDir, LOAD_PATH) || push!(LOAD_PATH, thisDir)
# include("escape.jl")
# @everywhere using .escape
using ProgressMeter
using Printf
using MATLAB
using LinearAlgebra
using BSON: @save, @load
include("escape_num_dimer.jl")
    # push!(LOAD_PATH, "/Users/brandonbehring/Desktop/Leap_Frog_2019")
push!(LOAD_PATH, "/home/brandon_behring/Desktop/Leap_Frog_2019")
using ProgressMeter
using Printf
using MATLAB
using .escape_num_dimer

Energy=range(.3,stop=.49,length=20);
t_end = exp(5)
    # width = sqrt(3); height = 1
width = 2.25; height = 3
n_iter_P = 375;n_iter_Q =500; N = n_iter_P * n_iter_Q;
ArrP = range(-width, stop = width, length = n_iter_P)
ArrQ = range(0, stop = height, length = n_iter_Q)
mesh_P = [P for Q in ArrQ, P in ArrP]
mesh_Q = [Q for Q in ArrQ, P in ArrP]
mesh = [(P,Q ) for Q in ArrQ, P in ArrP]

@showprogress for H in Energy
println(H)
location = "/mnt/bdd38f66-9ece-451a-b915-952523c139d2/Escape/"
h = replace(@sprintf("%.8f_",H), "." => "_")
nQ = @sprintf("%d",n_iter_Q)
file_name = location * "Dimer_time_before_Escape_" * h * nQ * ".fig"
h_BSON = location*"escape_data" * h * nQ * ".bson"
# num_until_exit = @showprogress pmap(escape_num_dimer.escape_exit_num_dimer, mesh, t_end, Energy)

# logz = log1p.(num_until_exit)

# SAVE_DATA = Dict("escape_times" => num_until_exit, "mesh" => mesh, "ArrP" => ArrP, "ArrQ" => ArrQ, "Energy" => H)
# @save h_BSON SAVE_DATA

SAVE_DATA=Dict()

# escape_time=zeros(n_iter_Q,n_iter_P)
escape_time=zeros(n_iter_Q,n_iter_P)

# figure()
# @showprogress for k=1:N
#     escape_time[k]=escape_num_dimer.escape_exit_num_dimer( mesh_list[k], t_end, H)
#     SAVE_DATA[mesh_list[k]]= escape_time[k]

# end
@showprogress for k=1:n_iter_Q
    for j=1:n_iter_P
        Q2=mesh[k,j][1]; P2=mesh[k,j][2];
        p1=[Q2,P2]
        p2=[0,1]
        p3=[0,-1]
        tol=.1
        if norm(p1)<tol|| norm(p1-p2)<tol || norm(p1-p3)<tol
            escape_time[k,j]=t_end
        end
        escape_time[k,j]=escape_num_dimer.escape_exit_num_dimer( mesh[k,j], t_end, H)
        SAVE_DATA[(ArrQ[k],ArrP[j])]=escape_time[k,j]
        # current_color=Brandons_Colors[mod(k,length(Brandons_Colors))+1]
    end
end
@save h_BSON SAVE_DATA

logz =log1p.(escape_time)
h_title=@sprintf("%.8f",H)

mat"figure();set(gcf, 'Position',  [0, 0, 1500, 1500]); hold on;"
mat"title($h_title)"
    # # # # mat"plot($P_0,$Q_0 ,'b.','MarkerSize',3)"
    # # # # mat"plot($P_1,$Q_1 ,'r.','MarkerSize',3)"
    # # # # mat"plot($P_2,$Q_2 ,'g.','MarkerSize',3)"
    # # # # mat"plot($P_4,$Q_4 ,'c.','MarkerSize',3)"
    # # # # mat"plot($P_5,$Q_5 ,'y.','MarkerSize',3)"
    # mat"imagesc([-$width,$width],[-$height,$height ],$logz)"
mat"imagesc([-$width,$width],[0,$height ],$logz)"
mat"imagesc([$width,-$width],[0,-$height ],$logz)"

mat"colorbar"

mat"axis([ -$width,$width,-$height,$height ])"
    # # # mat"axis([ -1,1,-1,1])"
mat"savefig($file_name)"

mat"figure();set(gcf, 'Position',  [0, 0, 1500, 1500]); hold on;"
mat"title($h_title)"
    # # # # mat"plot($P_0,$Q_0 ,'b.','MarkerSize',3)"
    # # # # mat"plot($P_1,$Q_1 ,'r.','MarkerSize',3)"
    # # # # mat"plot($P_2,$Q_2 ,'g.','MarkerSize',3)"
    # # # # mat"plot($P_4,$Q_4 ,'c.','MarkerSize',3)"
    # # # # mat"plot($P_5,$Q_5 ,'y.','MarkerSize',3)"
    # mat"imagesc([-$width,$width],[-$height,$height ],$logz)"
mat"contourf([-$width,$width],[0,$height ],$logz)"
mat"contourf([$width,-$width],[0,-$height ],$logz)"

mat"colorbar"

mat"axis([ -$width,$width,-$height,$height ])"
    # # # mat"axis([ -1,1,-1,1])"
end


text:
mf = (x -> 100, x -> -100)
met = FunctionalRegret(mf)
@test met(rand(1, 2), [1, 2], rand(1, 1)) == [0.0, 200.0]


text:
# Install RMS v0.3.0
import Pkg
Pkg.add(Pkg.PackageSpec(url="https://github.com/ReactionMechanismGenerator/ReactionMechanismSimulator.jl", rev="v0.3.0"))
# Pkg.build("ReactionMechanismSimulator")

max_stars_repo_name: APIOxy/APIOxy

max_issues_repo_name: APIOxy/APIOxy

max_forks_repo_name: APIOxy/APIOxy


text:
@testset "Hyperref" begin
    st = raw"""
       Some string
       $$ x = x \label{eq 1}$$
       then as per \citet{amari98b} also this \citep{bardenet17} and
       \cite{amari98b, bardenet17}
       Reference to equation: \eqref{eq 1}.

       Then maybe some text etc.

       * \biblabel{amari98b}{Amari and Douglas., 1998} **Amari** and **Douglas**: *Why Natural Gradient*, 1998.
       * \biblabel{bardenet17}{Bardenet et al., 2017} **Bardenet**, **Doucet** and **Holmes**: *On Markov Chain Monte Carlo Methods for Tall Data*, 2017.
    """;

    F.def_GLOBAL_VARS!()
    F.def_GLOBAL_LXDEFS!()

    m = F.convert_md(st)

    h1 = F.refstring("eq 1")
    h2 = F.refstring("amari98b")
    h3 = F.refstring("bardenet17")

    @test haskey(F.PAGE_EQREFS,  h1)
    @test haskey(F.PAGE_BIBREFS, h2)
    @test haskey(F.PAGE_BIBREFS, h3)

    @test F.PAGE_EQREFS[h1]  == 1 # first equation
    @test F.PAGE_BIBREFS[h2] == "Amari and Douglas., 1998"
    @test F.PAGE_BIBREFS[h3] == "Bardenet et al., 2017"

    h = F.convert_html(m)

    @test isapproxstr(h, """

          Some string
          \\[ x = x \\]
          then as per Amari and Douglas., 1998  also this (Bardenet et al., 2017)  and
          Amari and Douglas., 1998, Bardenet et al., 2017
          Reference to equation: (1) .

          Then maybe some text etc.

  Amari and Douglas: Why Natural Gradient, 1998.

  Bardenet, Doucet and Holmes: On Markov Chain Monte Carlo Methods for Tall Data, 2017.

        """)
end

@testset "Href-space" begin
    st = raw"""
       A
       $$ x = x \label{eq 1}$$
       B
       C \eqref{eq 1}.
       and *B $E$*.
    """
    h = st |> seval
    @test occursin(raw"""\[ x = x \]""", h)
    @test occursin(raw"""(1).""", h)
    @test occursin(raw"""B \(E\).""", h)
end

@testset "Eqref" begin
    st = raw"""
        \newcommand{\E}[1]{\mathbb E\left[#1\right]}
        \newcommand{\eqa}[1]{\begin{eqnarray}#1\end{eqnarray}}
        \newcommand{\R}{\mathbb R}
        Then something like
        \eqa{ \E{f(X)} \in \R &\text{if}& f:\R\maptso\R}
        and then
        \eqa{ 1+1 &=& 2 \label{eq:a trivial one}}
        but further
        \eqa{ 1 &=& 1 \label{beyond hope}}
        and finally a \eqref{eq:a trivial one} and maybe \eqref{beyond hope}.
        """
    m = F.convert_md(st, collect(values(F.GLOBAL_LXDEFS)))

    @test F.PAGE_EQREFS[F.PAGE_EQREFS_COUNTER] == 3
    @test F.PAGE_EQREFS[F.refstring("eq:a trivial one")] == 2
    @test F.PAGE_EQREFS[F.refstring("beyond hope")] == 3

    h1 = F.refstring("eq:a trivial one")
    h2 = F.refstring("beyond hope")

    m == "
Then something like  \$\$\\begin{array}{c}  \\mathbb E\\left[ f(X)\\right] \\in \\mathbb R &\\text{if}& f:\\mathbb R\\maptso\\mathbb R\\end{array}\$\$ and then  \$\$\\begin{array}{c}  1+1 &=&2 \\end{array}\$\$ but further  \$\$\\begin{array}{c}  1 &=& 1 \\end{array}\$\$ and finally a  ({{href EQR $h1}}) and maybe  ({{href EQR $h2}}).
\n"
end


text:
# ---
# title: 1525. Number of Good Ways to Split a String
# id: problem1525
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String, Bit Manipulation
# link: <https://leetcode.com/problems/number-of-good-ways-to-split-a-string/description/>
# hidden: true
# ---
#
# You are given a string `s`, a split is called _good_  if you can split `s`
# into 2 non-empty strings `p` and `q` where its concatenation is equal to `s`
# and the number of distinct letters in `p` and `q` are the same.
#
# Return the number of _good_ splits you can make in `s`.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "aacaba"
#     Output: 2
#     Explanation: There are 5 ways to split "aacaba" and 2 of them are good.
#     ("a", "acaba") Left string and right string contains 1 and 3 different letters respectively.
#     ("aa", "caba") Left string and right string contains 1 and 3 different letters respectively.
#     ("aac", "aba") Left string and right string contains 2 and 2 different letters respectively (good split).
#     ("aaca", "ba") Left string and right string contains 2 and 2 different letters respectively (good split).
#     ("aacab", "a") Left string and right string contains 3 and 1 different letters respectively.
#
#
# **Example 2:**
#
#
#
#     Input: s = "abcd"
#     Output: 1
#     Explanation: Split the string as follows ("ab", "cd").
#
#
# **Example 3:**
#
#
#
#     Input: s = "aaaaa"
#     Output: 4
#     Explanation: All possible splits are good.
#
# **Example 4:**
#
#
#
#     Input: s = "acbadbaada"
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `s` contains only lowercase English letters.
#   * `1 <= s.length <= 10^5`
#
#
## @lc code=start
using LeetCode

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


text:
"""
    Code for convenience functions for calculating seprate CCFs for each chunk in an AbstractChunkList (see RvSpectMLBase)
Author: Eric Ford
Created: August 2020
"""

"""  `calc_order_ccfs_chunklist ( chunklist_timeseries, list_of_ccf_plans )`
Convenience function to compute separate CCFs for each chunk (potentially an order or view around one or two lines) in a spectrum.
CCF is evaluated using line list and mask_shape provided by the ccf plan for each chunk.
# Inputs:
- `chunklist_timeseries`:
- `list_of_ccf_plans`: ccf plans (one for each chunk)
# Optional Arguments:
- `assume_sorted`:  if true, skips checking the line_list is sorted by wavelength
# Return:
A 2-d array containing the CCF at each (velocity, chunk)
"""
function calc_order_ccfs_chunklist!(chunk_ccfs_out::AbstractArray{T1,2}, chunk_list::AbstractChunkList,
    plan_for_chunk::AbstractVector{PlanT} = BasicCCFPlan(); assume_sorted::Bool = false  ) where {
                        T1<:Real, PlanT<:AbstractCCFPlan }
    num_chunks = length(chunk_list)
    @assert length(plan_for_chunk) == num_chunks
    @assert size(chunk_ccfs_out,2) == num_chunks
    @assert assume_sorted || issorted( first(plan_for_chunk).line_list.λ )
    num_vels = maximum(map(chid->calc_length_ccf_v_grid(plan_for_chunk[chid]), 1:num_chunks ))
    @assert size(chunk_ccfs_out,1) == num_vels
    #flush(stdout); println("num_vels = ",num_vels," num_chunks = ", num_chunks, " length(chunk_list.data)=", length(chunk_list.data));     flush(stdout);
    for chid in 1:num_chunks
        # chunk_ccfs_out[:,i] .=
        #flush(stdout);         println("chid = ", chid);         flush(stdout);
        calc_ccf_chunk!(view(chunk_ccfs_out, 1:num_vels, chid ), chunk_list.data[chid], plan_for_chunk[chid], assume_sorted=assume_sorted )
    end
    return chunk_ccfs_out
end

"""  `calc_order_ccfs_chunklist ( chunklist_timeseries, list_of_ccf_plans )`
Convenience function to compute separate CCFs for each chunk (potentially an order or view around one or two lines) in a spectrum.
CCF is evaluated using line list and mask_shape provided by the ccf plan for each chunk.
# Inputs:
- `chunklist_timeseries`:
- `list_of_ccf_plans`: ccf plans (one for each chunk)
# Optional Arguments:
- `assume_sorted`:  if true, skips checking the line_list is sorted by wavelength
# Return:
A 2-d array containing the CCF at each (velocity, chunk)
"""
function calc_order_ccfs_chunklist(chunk_list::AbstractChunkList,
                plan_for_chunk::AbstractVector{PlanT} = BasicCCFPlan();
                Δfwhm::Real = 0, assume_sorted::Bool = false ) where {
                        PlanT<:AbstractCCFPlan }
    num_chunks = length(chunk_list)
    @assert length(plan_for_chunk) == num_chunks
    @assert assume_sorted || issorted( first(plan_for_chunk).line_list.λ )
    num_vels = maximum(map(chid->calc_length_ccf_v_grid(plan_for_chunk[chid]), 1:num_chunks ))
    if Δfwhm > 0
       this_plan_for_chunk = copy(plan_for_chunk)
       #increase_mask_fwhm!(this_plan_for_chunk,Δfwhm)
       map(chid->increase_mask_fwhm!(this_plan_for_chunk[chid],Δfwhm), 1:num_chunks )
    else
       this_plan_for_chunk = plan_for_chunk
    end

    chunk_ccfs_out = zeros( num_vels, num_chunks )
    calc_order_ccfs_chunklist!(chunk_ccfs_out, chunk_list, this_plan_for_chunk, assume_sorted=assume_sorted )
    return chunk_ccfs_out
end

"""  `calc_order_ccfs_chunklist ( chunklist_timeseries, list_of_ccf_plans )`
Convenience function to compute separate CCFs for each chunk (potentially an order or view around one or two lines) in a spectrum.
CCF is evaluated using line list and mask_shape provided by the ccf plan for each chunk.
# Inputs:
- `chunklist_timeseries`:
- `list_of_ccf_plans`: ccf plans (one for each chunk)
# Optional Arguments:
- `assume_sorted`:  if true, skips checking the line_list is sorted by wavelength
# Return:
A 2-d array containing the CCF at each (velocity, chunk)
"""
function calc_order_ccf_and_vars_chunklist!(chunk_ccfs_out::AbstractArray{T1,2}, chunk_ccf_vars_out::AbstractArray{T1,2}, chunk_list::AbstractChunkList,
    plan_for_chunk::AbstractVector{PlanT} = BasicCCFPlan(); assume_sorted::Bool = false ) where {
                        T1<:Real, PlanT<:AbstractCCFPlan }
    num_chunks = length(chunk_list)
    @assert length(plan_for_chunk) == num_chunks
    @assert size(chunk_ccfs_out,2) == num_chunks
    @assert assume_sorted || issorted( first(plan_for_chunk).line_list.λ )
    num_vels = maximum(map(chid->calc_length_ccf_v_grid(plan_for_chunk[chid]), 1:num_chunks ))
    @assert size(chunk_ccfs_out,1) == num_vels

    for chid in 1:num_chunks
        calc_ccf_and_var_chunk!(view(chunk_ccfs_out, 1:num_vels, chid ), view(chunk_ccf_vars_out, 1:num_vels, chid ),
            chunk_list.data[chid], plan_for_chunk[chid], assume_sorted=assume_sorted )
    end
    return (ccfs=chunk_ccfs_out, ccf_vars=chunk_ccf_vars_out)
end

"""  `calc_order_ccfs_chunklist ( chunklist_timeseries, list_of_ccf_plans )`
Convenience function to compute separate CCFs for each chunk (potentially an order or view around one or two lines) in a spectrum.
CCF is evaluated using line list and mask_shape provided by the ccf plan for each chunk.
# Inputs:
- `chunklist_timeseries`:
- `list_of_ccf_plans`: ccf plans (one for each chunk)
# Optional Arguments:
- `assume_sorted`:  if true, skips checking the line_list is sorted by wavelength
# Return:
A 2-d array containing the CCF at each (velocity, chunk)
"""
function calc_order_ccf_and_vars_chunklist(chunk_list::AbstractChunkList,
            plan_for_chunk::AbstractVector{PlanT} = BasicCCFPlan();
            Δfwhm::Real = 0, assume_sorted::Bool = false ) where {
                        PlanT<:AbstractCCFPlan }
    num_chunks = length(chunk_list)
    @assert length(plan_for_chunk) == num_chunks
    @assert assume_sorted || issorted( first(plan_for_chunk).line_list.λ )

    if Δfwhm > 0
       this_plan_for_chunk = copy(plan_for_chunk)
       #increase_mask_fwhm!(this_plan_for_chunk,Δfwhm)
       map(chid->increase_mask_fwhm!(this_plan_for_chunk[chid],Δfwhm), 1:num_chunks )
    else
       this_plan_for_chunk = plan_for_chunk
    end

    num_vels = maximum(map(chid->calc_length_ccf_v_grid(plan_for_chunk[chid]), 1:num_chunks ))
    chunk_ccfs_out = zeros( num_vels, num_chunks )
    chunk_ccf_vars_out = zeros( num_vels, num_chunks )

    calc_order_ccf_and_vars_chunklist!(chunk_ccfs_out, chunk_ccf_vars_out, chunk_list, this_plan_for_chunk, assume_sorted=assume_sorted )
    return (ccfs=chunk_ccfs_out, ccf_vars=chunk_ccf_vars_out)
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
# # Knapsack Problem - Branch and Price
# This example demonstrates how the branch and price algorithm can be used to
# solve a JuDGE model.
using JuMP, JuDGE, GLPK, Random
# We will randomly generate data for this problem, but to ensure a consistent
# result we first set the random seed.
Random.seed!(50)

# We specify that there are 5 possible expansion investments, and 10 items to
# select from at each node.
numinvest = 5
numitems = 10

# The tree is defined to have a depth and 4 and a degree of 3.
mytree = narytree(4,3)

# Find the number of nodes in the tree
totalnodes = JuDGE.count(mytree)

# ## Data
# Set up the investment cost data
investcost = zeros(totalnodes,numinvest)
for i = 1:totalnodes
  investcost[i,:] = ([1,1.8,3.5,6.8,13.5])*(1-((i-1)/(totalnodes*1.2)))
end

# Specify the investment volumes using a binary expansion, with an initial capacity
# of 0
investvol = [1,2,4,8,16]
initialcap = 0

# Randomly specify the item volumes and costs
itemvolume = zeros(totalnodes,numitems)
for i = 1:totalnodes
  itemvolume[i,:] = ((rand(numitems))*2) + collect(range(4,22,length = numitems))
end

itemcost = zeros(totalnodes,numitems)
for i = 1:totalnodes
  itemcost[i,:] = ((rand(numitems) .- 0.5)*2)*2# + collect(range(0.5,1,length = numitems))
end

# Create a function that maps the node, and data name to a value
nodes = collect(mytree)
function data(node, input)
  input[findall(x -> x == node, nodes)[1], :]
end

# ## Subproblem definitions
# Define the subproblems as a function of the node
JuDGE_SP_Solver = optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => 0, "mip_gap" => 0.0)
function sub_problems(node)
  model = Model(JuDGE_SP_Solver)
  @expansion(model, bag[1:numinvest])
  @expansioncosts(model, sum(data(node,investcost)[i] * bag[i] for i in  1:numinvest))
  @variable(model, y[1:numitems], Bin)
  @constraint(model, BagExtension ,sum( y[i]*data(node,itemvolume)[i] for i in 1:numitems) <=
  				initialcap + sum(bag[i]*investvol[i] for i in 1:numinvest))
  @objective(model, Min, sum(-data(node,itemcost)[i] * y[i] for i in 1:numitems))
  model
end

# ## Defining and solving the JuDGE model
# Set up the JuDGE model.
JuDGE_MP_Solver = optimizer_with_attributes((method=GLPK.INTERIOR) -> GLPK.Optimizer(),
							"msg_lev" => 0, "mip_gap" => 0.0)
judy = JuDGEModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_MP_Solver)

# Solve the JuDGE model using the branch and price algorithm, using a
# constraint branch and branching on the lowest lower bound nodes first
best=JuDGE.branch_and_price(judy,rlx_abstol=10^-7,inttol=10^-6,
				branch_method=JuDGE.constraint_branch,search=:lowestLB)

# ## Displaying the output
# Print the objective function value
println("Objective: "*string(best.bounds.UB))

# Create a custom function that formats the output for print_expansions(...),
# and print the expansions using this function
function format_output(s::Symbol,values)
  if s==:bag
	 return sum(values[i]*investvol[i] for i in 1:numinvest)
  end
  return nothing
end

JuDGE.print_expansions(best,format=format_output)

# ## Deterministic equivalent
# Set up and solve the deterministic equivalent problem
JuDGE_DE_Solver = optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => 2, "mip_gap" => 3)
deteq = DetEqModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_DE_Solver)
JuDGE.solve(deteq)
println("Deterministic Equivalent Objective: " * string(objective_value(deteq.problem)))


text:
module Cerberus

import DataStructures,
    DisjunctiveConstraints,
    Gurobi,
    Logging,
    MathOptInterface,
    Printf,
    SparseArrays
const MOI = MathOptInterface
const MOIU = MOI.Utilities

const VI = MOI.VariableIndex
const SV = MOI.SingleVariable
const VOV = MOI.VectorOfVariables
const SAT = MOI.ScalarAffineTerm{Float64}
const SAF = MOI.ScalarAffineFunction{Float64}
const VAT = MOI.VectorAffineTerm{Float64}
const VAF = MOI.VectorAffineFunction{Float64}
const ET = MOI.EqualTo{Float64}
const GT = MOI.GreaterThan{Float64}
const LT = MOI.LessThan{Float64}
const IN = MOI.Interval{Float64}
const CI = MOI.ConstraintIndex
const ZO = MOI.ZeroOne
const GI = MOI.Integer

"Sets allowed in linear constraints."
const _C_SETS = Union{ET,GT,LT}

"Allowed objective types."
const _O_FUNCS = Union{SV,SAF}

"Sets allowed as variable bounds."
const _V_BOUND_SETS = Union{ET,GT,LT,IN}

"Integrality sets allowed."
const _V_INT_SETS = Union{Nothing,ZO,GI}

include("problem.jl")
include("tree.jl")

include("algorithm/algorithm_state.jl")
include("algorithm/config.jl")
include("algorithm/results.jl")
include("algorithm/lp_model.jl")
include("algorithm/util.jl")
include("algorithm/branch_and_bound.jl")
include("algorithm/branching.jl")
include("algorithm/logging.jl")

include("algorithm/formulaters/disjunctive_formulaters.jl")
include("algorithm/formulaters/naive_big_m.jl")

include("MOI/optimizer.jl")
include("MOI/variable.jl")
include("MOI/objective.jl")
include("MOI/linear_constraint.jl")
include("MOI/disjunctive_constraint.jl")
include("MOI/solution.jl")

end # module


text: include("staticarrays.jl")


text:
# Example of logistic regration with Mosek, MosekTools, and JuMP
# See Mosek manual 5.4.4 Logistic regression

using JuMP, Mosek, MosekTools
using RDatasets, Plots

# Iris dataset:  Our goal is to predict whether the iris species is versicolor
# using the sepal length and width and petal length and width.
iris = dataset("datasets", "iris");
Y = [species == "versicolor" ? 1.0 : -1.0 for species in iris.Species];
# Create feature matrix with one column for each feature
X = hcat(ones(size(iris, 1)), iris.SepalLength, iris.SepalWidth, iris.PetalLength, iris.PetalWidth);

# model size
n, p = size(X)
# regularizer
λ = 1e-2
# create model
model = Model(Mosek.Optimizer)
# variables
@variable(model, β[1:p])
@variable(model, β0)
# auxiliary variables for conic constraints
@variable(model, t[1:n])
@variable(model, r)
@variable(model, u̅[1:n])
@variable(model, u̲[1:n])
@variable(model, v̅[1:n])
@variable(model, v̲[1:n])

# minimize negative likelihood with a regularizer
@objective(model, Min, sum(t) + λ*r)
# conic constraints
for i in 1:n
    if Y[i] ==  1
        @constraint(model, u̅[i] + v̅[i] <= 1)
        @constraint(model, [-β0-X[i,:]'*β-t[i],1,u̅[i]] in MOI.ExponentialCone())
        @constraint(model, [-t[i],1,v̅[i]] in MOI.ExponentialCone())
    else
        @constraint(model, u̲[i] + v̲[i] <= 1)
        @constraint(model, [β0+X[i,:]'*β-t[i],1,u̲[i]] in MOI.ExponentialCone())
        @constraint(model, [-t[i],1,v̲[i]] in MOI.ExponentialCone())
    end
end
@constraint(model, [r; β] in SecondOrderCone())

# optimize
optimize!(model)
status=termination_status(model)
_β=JuMP.value.(β)
_β0=JuMP.value.(β0)
_t=sort!(JuMP.value.(t))

σ = 1 ./ (1 .+ exp.(-_β0.-X * _β))
if minimum(Y) == -1
  ŷ = (σ .> 0.5) .+ 2 .- 3 # bring the prediction in [-1, 1] range
else
  ŷ = (σ .> 0.5)
end
acc = count(Y .== ŷ) / n
print("accuracy: ", acc)

perm = sortperm(vec(X * _β));
plot(1:n, (Y[perm] .+ 1)/2, st=:scatter)
plot!(1:n, σ[perm])


text:
function sdpt1_test()
    # min X[1,1] + X[2,2]    max y
    #     X[2,1] = 1         [0   y/2     [ 1  0
    #                         y/2 0    <=   0  1]
    #     X >= 0              y free
    # Optimal solution:
    #     ⎛ 1   1 ⎞
    # X = ⎜       ⎟           y = 2
    #     ⎝ 1   1 ⎠

    model = TestModel{Float64}()

    X = MOI.add_variables(model, 3)

    vov = MOI.VectorOfVariables(X)

    cX = MOI.add_constraint(model, vov, MOI.PositiveSemidefiniteConeTriangle(2))

    c = MOI.add_constraint(model, MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, X[2])], 0.0), MOI.EqualTo(1.0))

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [X[1], X[end]]), 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return  model
end

function sdpt2_test()
    #=
        min TR(X)
    s.t.
        X[2,1] = 1
        X in PSD
    =#

    model = TestModel{Float64}()

    X = MOI.add_variables(model, 3)

    vov = MOI.VectorOfVariables(X)

    cX = MOI.add_constraint(model, MOI.VectorAffineFunction{Float64}(vov), MOI.PositiveSemidefiniteConeTriangle(2))

    c = MOI.add_constraint(model, MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, X[2])], 0.0), MOI.EqualTo(1.0))

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [X[1], X[end]]), 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return  model
end

function sdpt3_test()
    # Problem SDP1 - sdo1 from MOSEK docs
    # From Mosek.jl/test/mathprogtestextra.jl, under license:
    #   Copyright (c) 2013 Ulf Worsoe, Mosek ApS
    #   Permission is hereby granted, free of charge, to any person obtaining a copy of this
    #   software and associated documentation files (the "Software"), to deal in the Software
    #   without restriction, including without limitation the rights to use, copy, modify, merge,
    #   publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons
    #   to whom the Software is furnished to do so, subject to the following conditions:
    #   The above copyright notice and this permission notice shall be included in all copies or
    #   substantial portions of the Software.
    #
    #     | 2 1 0 |
    # min | 1 2 1 | . X + x1
    #     | 0 1 2 |
    #
    #
    # s.t. | 1 0 0 |
    #      | 0 1 0 | . X + x1 = 1
    #      | 0 0 1 |
    #
    #      | 1 1 1 |
    #      | 1 1 1 | . X + x2 + x3 = 1/2
    #      | 1 1 1 |
    #
    #      (x1,x2,x3) in C^3_q
    #      X in C_psd
    #
    # The dual is
    # max y1 + y2/2
    #
    # s.t. | y1+y2    y2    y2 |
    #      |    y2 y1+y2    y2 | in C_psd
    #      |    y2    y2 y1+y2 |
    #
    #      (1-y1, -y2, -y2) in C^3_q
    #
    # The dual of the SDP constraint is rank two of the form
    # [γ, 0, -γ] * [γ, 0, γ'] + [δ, ε, δ] * [δ, ε, δ]'
    # and the dual of the SOC constraint is of the form (√2*y2, -y2, -y2)
    #
    # The feasible set of the constraint dual contains only four points.
    # Eliminating, y1, y2 and γ from the dual constraints gives
    # -ε^2 + -εδ + 2δ^2 + 1
    # (√2-2)ε^2 + (-2√2+2)δ^2 + 1
    # Eliminating ε from this set of equation give
    # (-6√2+4)δ^4 + (3√2-2)δ^2 + (2√2-3)
    # from which we find the solution
    # δ = √(1 + (3*√2+2)*√(-116*√2+166) / 14) / 2
    # which is optimal
    # ε = √((1 - 2*(√2-1)*δ^2) / (2-√2))
    # y2 = 1 - ε*δ
    # y1 = 1 - √2*y2
    # obj = y1 + y2/2
    # The primal solution is rank one of the form
    # X = [α, β, α] * [α, β, α]'
    # and by complementary slackness, x is of the form (√2*x2, x2, x2)
    # The primal reduces to
    #      4α^2+4αβ+2β^2+√2*x2= obj
    #      2α^2    + β^2+√2*x2 = 1 (1)
    #      8α^2+8αβ+2β^2+ 4 x2 = 1
    # Eliminating β, we get
    # 4α^2 + 4x2 = 3 - 2obj (2)
    # By complementary slackness, we have β = kα where
    # k = -2*δ/ε
    # Replacing β by kα in (1) allows to eliminate α^2 in (2) to get
    # x2 = ((3-2obj)*(2+k^2)-4) / (4*(2+k^2)-4*√2)
    # With (2) we get
    # α = √(3-2obj-4x2)/2
    # β = k*α

   model = TestModel{Float64}()

    X = MOI.add_variables(model, 6)
    x = MOI.add_variables(model, 3)

    vov = MOI.VectorOfVariables(X)

    cX = MOI.add_constraint(model, vov, MOI.PositiveSemidefiniteConeTriangle(3))
    cx = MOI.add_constraint(model, MOI.VectorOfVariables(x), MOI.SecondOrderCone(3))

    c1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1., 1, 1, 1], [X[1], X[3], X[end], x[1]]), 0.), MOI.EqualTo(1.))
    c2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1., 2, 1, 2, 2, 1, 1, 1], [X; x[2]; x[3]]), 0.), MOI.EqualTo(1/2))

    objXidx = [1:3; 5:6]
    objXcoefs = 2*ones(5)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([objXcoefs; 1.0], [X[objXidx]; x[1]]), 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end

function sdpt4_test()
    # Problem SDP1 - sdo1 from MOSEK docs
    # From Mosek.jl/test/mathprogtestextra.jl, under license:
    #   Copyright (c) 2013 Ulf Worsoe, Mosek ApS
    #   Permission is hereby granted, free of charge, to any person obtaining a copy of this
    #   software and associated documentation files (the "Software"), to deal in the Software
    #   without restriction, including without limitation the rights to use, copy, modify, merge,
    #   publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons
    #   to whom the Software is furnished to do so, subject to the following conditions:
    #   The above copyright notice and this permission notice shall be included in all copies or
    #   substantial portions of the Software.
    #
    #     | 2 1 0 |
    # min | 1 2 1 | . X + x1
    #     | 0 1 2 |
    #
    #
    # s.t. | 1 0 0 |
    #      | 0 1 0 | . X + x1 = 1
    #      | 0 0 1 |
    #
    #      | 1 1 1 |
    #      | 1 1 1 | . X + x2 + x3 = 1/2
    #      | 1 1 1 |
    #
    #      (x1,x2,x3) in C^3_q
    #      X in C_psd
    #
    # The dual is
    # max y1 + y2/2
    #
    # s.t. | y1+y2    y2    y2 |
    #      |    y2 y1+y2    y2 | in C_psd
    #      |    y2    y2 y1+y2 |
    #
    #      (1-y1, -y2, -y2) in C^3_q
    #
    # The dual of the SDP constraint is rank two of the form
    # [γ, 0, -γ] * [γ, 0, γ'] + [δ, ε, δ] * [δ, ε, δ]'
    # and the dual of the SOC constraint is of the form (√2*y2, -y2, -y2)
    #
    # The feasible set of the constraint dual contains only four points.
    # Eliminating, y1, y2 and γ from the dual constraints gives
    # -ε^2 + -εδ + 2δ^2 + 1
    # (√2-2)ε^2 + (-2√2+2)δ^2 + 1
    # Eliminating ε from this set of equation give
    # (-6√2+4)δ^4 + (3√2-2)δ^2 + (2√2-3)
    # from which we find the solution
    # δ = √(1 + (3*√2+2)*√(-116*√2+166) / 14) / 2
    # which is optimal
    # ε = √((1 - 2*(√2-1)*δ^2) / (2-√2))
    # y2 = 1 - ε*δ
    # y1 = 1 - √2*y2
    # obj = y1 + y2/2
    # The primal solution is rank one of the form
    # X = [α, β, α] * [α, β, α]'
    # and by complementary slackness, x is of the form (√2*x2, x2, x2)
    # The primal reduces to
    #      4α^2+4αβ+2β^2+√2*x2= obj
    #      2α^2    + β^2+√2*x2 = 1 (1)
    #      8α^2+8αβ+2β^2+ 4 x2 = 1
    # Eliminating β, we get
    # 4α^2 + 4x2 = 3 - 2obj (2)
    # By complementary slackness, we have β = kα where
    # k = -2*δ/ε
    # Replacing β by kα in (1) allows to eliminate α^2 in (2) to get
    # x2 = ((3-2obj)*(2+k^2)-4) / (4*(2+k^2)-4*√2)
    # With (2) we get
    # α = √(3-2obj-4x2)/2
    # β = k*α

    model = TestModel{Float64}()

    X = MOI.add_variables(model, 6)
    x = MOI.add_variables(model, 3)

    vov = MOI.VectorOfVariables(X)
    cX = MOI.add_constraint(model, MOI.VectorAffineFunction{Float64}(vov), MOI.PositiveSemidefiniteConeTriangle(3))
    cx = MOI.add_constraint(model, MOI.VectorOfVariables(x), MOI.SecondOrderCone(3))

    c1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1., 1, 1, 1], [X[1], X[3], X[end], x[1]]), 0.), MOI.EqualTo(1.))
    c2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1., 2, 1, 2, 2, 1, 1, 1], [X; x[2]; x[3]]), 0.), MOI.EqualTo(1/2))

    objXidx = [1:3; 5:6]
    objXcoefs = 2*ones(5)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([objXcoefs; 1.0], [X[objXidx]; x[1]]), 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end


text:
# TODO: add to documentation page and do tests

#######################################################################################################################################################################################################
#
# Changes to this struct
# General
#     2021-Aug-04: refactor the structure with constants, variables, and temporary cache
#     2021-Aug-04: add concentrations and characteristic curves altogether
#     2021-Aug-10: add CBC and PRO supoort
#     2021-Agu-10: add constructors within the structure rather than initialize it externally
#     2021-Sep-30: rename LeafBio to LeafBiophysics to be more specific
#     2021-Oct-19: sort variable to prognostic and dignostic catergories
#     2021-Oct-21: rename f_sense and K_SENES to brown and K_BROWN
#     2021-Nov-24: tease apart the characteristic absorption curves to HyperspectralAbsorption
#
#######################################################################################################################################################################################################
"""
$(TYPEDEF)

Struct that contains leaf biophysical traits used to run leaf reflection and transmittance.

# Fields
$(TYPEDFIELDS)
"""
mutable struct HyperspectralAbsorption{FT<:AbstractFloat}
    # parameters that do not change with time
    "Specific absorption coefficients of anthocynanin `[-]`"
    K_ANT::Vector{FT}
    "Specific absorption coefficients of senescent material (brown pigments) `[-]`"
    K_BROWN::Vector{FT}
    "Specific absorption coefficients of chlorophyll a and b `[-]`"
    K_CAB::Vector{FT}
    "Specific absorption coefficients of violaxanthin carotenoid `[-]`"
    K_CAR_V::Vector{FT}
    "Specific absorption coefficients of zeaxanthin carotenoid `[-]`"
    K_CAR_Z::Vector{FT}
    "Specific absorption coefficients of carbon-based constituents `[-]`"
    K_CBC::Vector{FT}
    "Specific absorption coefficients of water `[-]`"
    K_H₂O::Vector{FT}
    "Specific absorption coefficients of dry matter `[-]`"
    K_LMA::Vector{FT}
    "Specific absorption coefficients of protein `[-]`"
    K_PRO::Vector{FT}
    "Specific absorption coefficients of PS I and II `[-]`"
    K_PS::Vector{FT}
    "Refractive index `[-]`"
    NR::Vector{FT}
end

"""
    HyperspectralAbsorption{FT}(wls::WaveLengthSet; opti::String=OPTI_2021) where {FT<:AbstractFloat}

Constructor for [`HyperspectralAbsorption`](@ref), given
- `wls` [`WaveLengthSet`](@ref) type structure
- `opti` Path to leaf optical properties

---
# Examples
```julia
ha = HyperspectralAbsorption{FT}();
ha = HyperspectralAbsorption{FT}(WaveLengthSet{FT}(collect(400:50:2400)));
ha = HyperspectralAbsorption{FT}(WaveLengthSet{FT}(collect(400:50:2400)); opti=ClimaCache.OPTI_2017);
```
"""
HyperspectralAbsorption{FT}(wls::WaveLengthSet = WaveLengthSet{FT}(); opti::String=OPTI_2021) where {FT<:AbstractFloat} = (
    @unpack NΛ, NΛ_SIF, NΛ_SIFE, SΛ = wls;

    # read data from the MAT file
    _opti    = matread(opti)["optipar"];
    __nr     = _opti["nr"  ];
    __Klma   = _opti["Kdm" ];
    __Kcab   = _opti["Kab" ];
    __Kant   = _opti["Kant"];
    __Kh2o   = _opti["Kw"  ];
    __Kbrown = _opti["Ks"  ];
    __Kps    = _opti["phi" ];
    __KcaV   = _opti["KcaV"];
    __KcaZ   = _opti["KcaZ"];
    __lambda = _opti["wl"  ];

    # read protein and carbon-based constituents if exist (not in file OPTI_2017)
    if opti == OPTI_2017
        __Kpro = __Klma;
        __Kcbc = __Klma;
    else
        __Kpro = _opti["Kp"  ];
        __Kcbc = _opti["Kcbc"];
    end;

    # create data to parse
    _nr     = zeros(FT, NΛ);
    _Klma   = zeros(FT, NΛ);
    _Kcab   = zeros(FT, NΛ);
    _Kant   = zeros(FT, NΛ);
    _Kh2o   = zeros(FT, NΛ);
    _Kbrown = zeros(FT, NΛ);
    _Kps    = zeros(FT, NΛ);
    _KcaV   = zeros(FT, NΛ);
    _KcaZ   = zeros(FT, NΛ);
    _Kpro   = zeros(FT, NΛ);
    _Kcbc   = zeros(FT, NΛ);

    # fill in the data arrays
    @inbounds for _i in 1:NΛ
        _wo = findall( SΛ[_i] .<= __lambda .< SΛ[_i+1] );
        if length(_wo) < 1 @warn "Warning, some wavelengths out of bounds $(string(SΛ[_i]))" end;

        _nr[_i]     = mean(__nr[_wo]    );
        _Klma[_i]   = mean(__Klma[_wo]  );
        _Kcab[_i]   = mean(__Kcab[_wo]  );
        _Kant[_i]   = mean(__Kant[_wo]  );
        _Kh2o[_i]   = mean(__Kh2o[_wo]  );
        _Kbrown[_i] = mean(__Kbrown[_wo]);
        _Kps[_i]    = mean(__Kps[_wo]   );
        _KcaV[_i]   = mean(__KcaV[_wo]  );
        _KcaZ[_i]   = mean(__KcaZ[_wo]  );
        _Kpro[_i]   = mean(__Kpro[_wo]  );
        _Kcbc[_i]   = mean(__Kcbc[_wo]  );
    end;

    return HyperspectralAbsorption{FT}(_Kant, _Kbrown, _Kcab, _KcaV, _KcaZ, _Kcbc, _Kh2o, _Klma, _Kpro, _Kps, _nr)
);