text:
# Downward initialization through a cascading model
# these functions need to be consolidated refactored and deprecated accordingly
# see #910

## ============================================================================================
## Down init sibling priority order
## ============================================================================================

"""
$SIGNATURES

Test waiting order between siblings for cascading downward tree initialization.

Notes
- State machine function 8j.

DevNotes
- FIXME, this guy is never called during WIP of 459 dwnMsg consolidation
- FIXME, something wrong with CSM sequencing, https://github.com/JuliaRobotics/IncrementalInference.jl/issues/602#issuecomment-682114232
- This might be replaced with 4-stroke tree-init.
"""
function dwnInitSiblingWaitOrder_StateMachine(csmc::CliqStateMachineContainer)

  prnt_ = getParent(csmc.tree, csmc.cliq)
  if 0 == length(prnt_)
    # go to 7
    return determineCliqNeedDownMsg_StateMachine
  end
  prnt = prnt_[1]

  opt = getSolverParams(csmc.cliqSubFg) # csmc.dfg

  # now get the newly computed message from the appropriate container
  # make sure this is a pull model #674 (pull msg from source/prnt)
  # FIXME must be consolidated as part of #459
  # dwinmsgs = getfetchCliqueInitMsgDown(getCliqueData(prnt), from=:dwnInitSiblingWaitOrder_StateMachine)
  dwinmsgs = fetchDwnMsgConsolidated(prnt)
  dwnkeys_ = collect(keys(dwinmsgs.belief))
  infocsm(csmc, "8j, dwnInitSiblingWaitOrder_StateMachine, dwinmsgs keys=$(dwnkeys_)")

  # add downward belief prop msgs
  msgfcts = addMsgFactors!(csmc.cliqSubFg, dwinmsgs, DownwardPass)
  ## update solveable dims with newly avaiable init information
  # determine if more info is needed for partial
  sdims = getCliqVariableMoreInitDims(csmc.cliqSubFg, csmc.cliq)
  infocsm(csmc, "8j, dwnInitSiblingWaitOrder_StateMachine, sdims=$(sdims)")
  updateCliqSolvableDims!(csmc.cliq, sdims, csmc.logger)

  _dbgCSMSaveSubFG(csmc, "fg_DWNCMN_8j")

  ## FIXME, if this new, and all sibling clique's solvableDim are 0, then go back to waitChangeOnParentCondition_StateMachine

  # go to 8o.i
  return testDirectDwnInit_StateMachine
end

"""
    $SIGNATURES

Can this clique initialize directly from available down message info?

Notes
- State machine function nr. 8o.i
- Assume must have parent since this only occurs after 8j.
"""
function testDirectDwnInit_StateMachine(csmc::CliqStateMachineContainer)

  prnt = getParent(csmc.tree, csmc.cliq)[1]
  dwinmsgs = fetchDwnMsgConsolidated(prnt)
  dwnkeys_ = collect(keys(dwinmsgs.belief))
  # NOTE, only use separators, not all parent variables (DF ???)
  # dwnkeys_ = lsf(csmc.cliqSubFg, tags=[:__DOWNWARD_COMMON__;]) .|> x->ls(csmc.cliqSubFg, x)[1]
  # @assert length(intersect(dwnkeys, dwnkeys_)) == length(dwnkeys) "split dwnkeys_ is not the same, $dwnkeys, and $dwnkeys_, separators: $(getCliqSeparatorVarIds(csmc.cliq))"

  # priorize solve order for mustinitdown with lowest dependency first
  # follow example from issue #344
  # mustwait = false
  if length(intersect(dwnkeys_, getCliqSeparatorVarIds(csmc.cliq))) == 0
    # can directly use DOWNWARD_COMMON
    infocsm(csmc, "8j, dwnInitSiblingWaitOrder_StateMachine, no can do, must wait for siblings to update parent first.")
    # mustwait = true

    # go to 8o.iv
    return sibsDwnPriorityInit_StateMachine
  end

  # go to 8o.ii
  return testDelayOrderDwnInit_StateMachine
end

"""
    $SIGNATURES

Return true if this clique's down init should be delayed on account of prioritization among sibling separators.

Notes
- State machine function nr. 8o.ii
- process described in issue #344

Dev Notes
- not priorizing order yet (TODO), just avoiding unsolvables at this time.
- Very closely related to 8o.iii -- refactor likely (NOTE).
- should precompute `allinters`.
- # FIXME ON FIRE, this is doing `getCliqueStatus` of siblings without following channels.
  - NOTE This is only used in old tree-init, not part of active regular solve code.
  - See #954
"""
function testDelayOrderDwnInit_StateMachine(csmc::CliqStateMachineContainer)

  prnt = getParent(csmc.tree, csmc.cliq)[1]
  dwinmsgs = fetchDwnMsgConsolidated(prnt)
  dwnkeys = collect(keys(dwinmsgs.belief))

  tree = csmc.tree
  cliq = csmc.cliq
  logger = csmc.logger

  # when is a cliq upsolved
  solvedstats = Symbol[:upsolved; :marginalized; :uprecycled]

  # safety net double check
  cliqst = getCliqueStatus(cliq)
  if cliqst in solvedstats
    infocsm(csmc, "getSiblingsDelayOrder -- clique status should not be here with a solved cliqst=$cliqst")
    # go to 8o.iii
    return testPartialNeedsDwnInit_StateMachine
  end

  # get siblings separators
  sibs = getCliqSiblings(tree, cliq, true)
  ids = map(s->s.index, sibs)
  len = length(sibs)
  sibidx = collect(1:len)[ids .== cliq.index][1]
  seps = getCliqSeparatorVarIds.(sibs)
  lielbls = setdiff(ids, cliq.index)
  # get intersect matrix of siblings (should be exactly the same across siblings' csm)
  allinters = Array{Int,2}(undef, len, len)
  dwninters = Vector{Int}(undef, len)
  infocsm(csmc, "getSiblingsDelayOrder -- number siblings=$(len), sibidx=$sibidx")

  # sum matrix with all "up solved" rows and columns eliminated
  fill!(allinters, 0)
  for i in 1:len
    for j in i:len
      if i != j
        allinters[i,j] = length(intersect(seps[i],seps[j]))
      end
    end
    dwninters[i] = length(intersect(seps[i], dwnkeys))
  end

  # sum "across/over" rows, then columns (i.e. visa versa "along" columns then rows)
  rows = sum(allinters, dims=1)
  cols = sum(allinters, dims=2)

  infocsm(csmc, "getSiblingsDelayOrder -- allinters=$(allinters), getSiblingsDelayOrder -- rows=$(rows), getSiblingsDelayOrder -- rows=$(cols)")

  # is this clique a non-zero row -- i.e. sum across columns? if not, no further special care needed
  if cols[sibidx] == 0
    infocsm(csmc, "getSiblingsDelayOrder -- cols[sibidx=$(sibidx))] == 0, no special care needed")
    # go to 8o.iii
    return testPartialNeedsDwnInit_StateMachine
  end

  # now determine if initializing from below or needdownmsg
  if cliqst in Symbol[:needdownmsg;]
    # be super careful about delay (true) vs pass (false) at this point -- might be partial too TODO
    # return true if delay beneficial to initialization accuracy

    # find which siblings this cliq epends on
    symm = allinters + allinters'
    maskcol = 0 .< symm[:,sibidx]
    # lenm = length(maskcol)
    stat = Vector{Symbol}(undef, len)
    stillbusymask = fill(false, len)

    # get each sibling status (entering atomic computation segment -- until wait command)
    stat .= getCliqueStatus.(sibs) #[maskcol]

    ## (long down chain case)
    # need different behaviour when all remaining siblings are blocking with :needdownmsg
    remainingmask = stat .== :needdownmsg
    if sum(remainingmask) == length(stat)
      infocsm(csmc, "getSiblingsDelayOrder -- all blocking: sum(remainingmask) == length(stat), stat=$stat")
      # pick sibling with most overlap in down msgs from parent
      # list of similar length siblings
      candidates = dwninters .== maximum(dwninters)
      if candidates[sibidx]
        # must also pick minimized intersect with other remaing siblings
        maxcan = collect(1:len)[candidates]
        infocsm(csmc, "getSiblingsDelayOrder -- candidates=$candidates, maxcan=$maxcan, rows=$rows")
        if rows[sibidx] == minimum(rows[maxcan])
          infocsm(csmc, "getSiblingsDelayOrder -- FORCE DOWN INIT SOLVE ON THIS CLIQUE: $(cliq.index), $(getLabel(cliq))")
          # go to 8o.iii
          return testPartialNeedsDwnInit_StateMachine
        end
      end
      infocsm(csmc, "getSiblingsDelayOrder -- not a max and should block")
      # go to 8o.iv
      return sibsDwnPriorityInit_StateMachine
    end

    # still busy solving on branches, so potential to delay
    for i in 1:len
      stillbusymask[i] = maskcol[i] && !(stat[i] in solvedstats)
    end
    infocsm(csmc, "getSiblingsDelayOrder -- busy solving: maskcol=$maskcol, stillbusy=$stillbusymask")

    # Too blunt -- should already have returned false by this point perhaps
    if 0 < sum(stillbusymask)
      # yes something to delay about
      infocsm(csmc, "getSiblingsDelayOrder -- yes delay, stat=$stat symm=$symm")

      # go to 8o.iv
      return sibsDwnPriorityInit_StateMachine
    end
  end

  infocsm(csmc, "getSiblingsDelayOrder -- default will not delay")
  # carry over default from partial init process

  # go to 8o.iii
  return testPartialNeedsDwnInit_StateMachine
end

"""
    $SIGNATURES

Return true if both, i.) this clique requires more downward information, ii.) more
downward message information could potentially become available.

Notes
- State machine function nr. 8o.iii
- Delay initialization to the last possible moment.

Dev Notes:
- # FIXME ON FIRE, this is doing `getCliqueStatus` of siblings without following channels.
  - NOTE This is only used in old tree-init, not part of active regular solve code.
  - See #954
- Determine clique truely isn't able to proceed any further:
  - should be as self reliant as possible (using clique's status as indicator)
  - OBSOLETE ??
    - change status to :mustinitdown if have only partial beliefs so far:
    - combination of status, while partials belief siblings are not :mustinitdown
"""
function testPartialNeedsDwnInit_StateMachine(csmc::CliqStateMachineContainer)
  #
  prnt = getParent(csmc.tree, csmc.cliq)[1]
  dwinmsgs = fetchDwnMsgConsolidated(prnt)

  tree = csmc.tree
  cliq = csmc.cliq
  logger = csmc.logger

  # which incoming messages are partials
  hasPartials = Dict{Symbol, Int}()
  for (sym, tmsg) in dwinmsgs.belief
    # assuming any down message per label that is not partial voids further partial consideration
    if sum(tmsg.inferdim) > 0
      if !haskey(hasPartials, sym)
        hasPartials[sym] = 0
      end
      hasPartials[sym] += 1
    end
  end
  partialKeys = collect(keys(hasPartials))

  ## determine who might be able to help init this cliq
  # check sibling separator sets against this clique's separator
  sibs = getCliqSiblings(tree, cliq)

  infocsm(csmc, "getCliqSiblingsPartialNeeds -- CHECK PARTIAL")
  # identify which cliques might have useful information
  localsep = getCliqSeparatorVarIds(cliq)
  seps = Dict{Int, Vector{Symbol}}()
  for si in sibs
    # @show getLabel(si)
    mighthave = intersect(getCliqSeparatorVarIds(si), localsep)
    if length(mighthave) > 0
      seps[si.index] = mighthave
      if getCliqueStatus(si) in [:initialized; :null; :needdownmsg]
        # partials treated special -- this is slightly hacky
        if length(intersect(localsep, partialKeys)) > 0 && length(mighthave) > 0
          # this sibling might have info to delay about
          setCliqueDrawColor!(cliq,"magenta")
          # go to 8o.iv
          return sibsDwnPriorityInit_StateMachine
        end
      end
    end
  end
  # determine if those cliques will / or will not be able to provide more info
  # when does clique change to :mustinitdown
  # default

  # go to 8e.ii.
  return tryDwnInitCliq_StateMachine
end

"""
$SIGNATURES

Return true there is no other sibling that will make progress.

Notes
- State machine function nr. 8o.iv
- Relies on sibling priority order with only one "currently best" option that will force progress in global upward inference.
- Return false if one of the siblings is still busy

DevNotes
- # FIXME ON FIRE, calling getCliqueStatus of siblings directly
  - Only used for older tree init, not part of regular up-down solving code
  - see #954
- Best is likely if parent actually made determination on which child to solve first in this case.
"""
function sibsDwnPriorityInit_StateMachine(csmc::CliqStateMachineContainer)

  tree = csmc.tree
  cliq = csmc.cliq
  prnt_ = getParent(tree, cliq)
  prnt = prnt_[1]

  dwinmsgs = fetchDwnMsgConsolidated(prnt)
  dwnkeys_ = collect(keys(dwinmsgs.belief))

  solvord = getCliqSiblingsPriorityInitOrder( csmc.tree, prnt, dwinmsgs, csmc.logger )

  # noOneElse = areSiblingsRemaingNeedDownOnly(csmc.tree, csmc.cliq)
  noOneElse = true
  stillbusylist = [:null; :initialized;]
  if 0 < length(prnt_)
    for si in getChildren(tree, prnt)
      # are any of the other siblings still busy?
      if si.index != cliq.index && getCliqueStatus(si) in stillbusylist
        noOneElse = false
      end
    end
  end
  # nope, everybody is waiting for something to change -- proceed with forcing a cliq solve

  infocsm(csmc, "8j, dwnInitSiblingWaitOrder_StateMachine, $(prnt.index), $noOneElse, solvord = $solvord")

  if csmc.cliq.index != solvord[1]
    # TODO, is this needed? fails hex if included with correction :needdowninit --> :needdownmsg
    # if dwinmsgs.status == :initialized && getCliqueStatus(csmc.cliq) == :needdownmsg  # used to be :needdowninit
    #   # go to 7e
    #   return slowWhileInit_StateMachine
    # end

    infocsm(csmc, "8j, dwnInitSiblingWaitOrder_StateMachine, must wait on change.")
    # remove all message factors
    # remove msg factors previously added
    fctstorm = deleteMsgFactors!(csmc.cliqSubFg, [:__DOWNWARD_COMMON__])
    infocsm(csmc, "8j, dwnInitSiblingWaitOrder_StateMachine, removing factors $fctstorm")

    # go to 8c
    return waitChangeOnParentCondition_StateMachine
  end

  # go to 8e.ii.
  return tryDwnInitCliq_StateMachine
end


text:
# using Lar
using SparseArrays

#=== GEOMETRY ===#
cubes = Lar.Model([
	-0.0797872 -0.1833684 0.0219132 -0.0816681 0.4465752 0.3429939 0.5482755 0.4446942 0.7187785 0.5249728 0.5799289 0.3861232 0.8613671 0.6675614 0.7225174 0.5287117 -0.2241867 0.1406597 -0.0449854 0.319861 -0.12314 0.2417064 0.0560613 0.4209077
	0.2943071 0.3960075 0.8243972 0.9260975 0.2118999 0.3136002 0.7419899 0.8436902 0.5392904 0.4004408 0.7762076 0.6373579 0.5812714 0.4424218 0.8181886 0.6793389 0.9338056 1.1130069 0.7580241 0.9372254 0.5985085 0.7777098 0.422727 0.6019283
	-0.1823453 0.344017 -0.2647526 0.2616098 -0.0628419 0.4635205 -0.1452491 0.3811132 0.1230172 0.2656057 0.1649982 0.3075868 0.3577031 0.5002916 0.3996841 0.5422726 0.4075841 0.5086308 0.0722869 0.1733336 0.6373704 0.7384171 0.3020733 0.40312
]);

#=== TOPOLOGY : 1-CELLS ===#
I = Array{Int64}([1, 5, 9, 1, 6, 10, 2, 5, 11, 2, 6, 12, 3, 7, 9, 3, 8, 10, 4, 7, 11, 4, 8, 12, 13, 17, 21, 13, 18, 22, 14, 17, 23, 14, 18, 24, 15, 19, 21, 15, 20, 22, 16, 19, 23, 16, 20, 24, 25, 29, 33, 25, 30, 34, 26, 29, 35, 26, 30, 36, 27, 31, 33, 27, 32, 34, 28, 31, 35, 28, 32, 36, ]);
J = Array{Int64}([1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 14, 14, 14, 15, 15, 15, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, ]);
K = Array{Int8}([-1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, -1, -1, -1, 1, 1, -1, 1, -1, 1, 1, 1, 1, 1, -1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, -1, -1, -1, 1, 1, -1, 1, -1, 1, 1, 1, 1, 1, -1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, -1, -1, -1, 1, 1, -1, 1, -1, 1, 1, 1, 1, 1, ]);

Lar.addModelCells!(cubes, 1, SparseArrays.sparse(I, J, K, 36, 24))

#=== TOPOLOGY : 2-CELLS ===#
I = Array{Int64}([1, 3, 1, 4, 2, 3, 2, 4, 1, 5, 1, 6, 2, 5, 2, 6, 3, 5, 3, 6, 4, 5, 4, 6, 7, 9, 7, 10, 8, 9, 8, 10, 7, 11, 7, 12, 8, 11, 8, 12, 9, 11, 9, 12, 10, 11, 10, 12, 13, 15, 13, 16, 14, 15, 14, 16, 13, 17, 13, 18, 14, 17, 14, 18, 15, 17, 15, 18, 16, 17, 16, 18, ]);
J = Array{Int64}([1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17, 18, 18, 19, 19, 20, 20, 21, 21, 22, 22, 23, 23, 24, 24, 25, 25, 26, 26, 27, 27, 28, 28, 29, 29, 30, 30, 31, 31, 32, 32, 33, 33, 34, 34, 35, 35, 36, 36, ]);
K = Array{Int8}([1, 1, -1, 1, 1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, 1, 1, -1, 1, 1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, 1, 1, -1, 1, 1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, ]);

Lar.addModelCells!(cubes, 2, SparseArrays.sparse(I, J, K, 18, 36))

#=== TOPOLOGY : 3-CELLS ===#
I = Array{Int64}([]);
J = Array{Int64}([]);
K = Array{Int8}([]);

Lar.addModelCells!(cubes, 3, SparseArrays.sparse(I, J, K, 0, 18))

ext: jl


text:
using Documenter, Dubins

makedocs(
    modules = [Dubins],
    format = :html,
    sitename = "Dubins",
    authors = "Kaarthik Sundar",
    pages = [
        "Home" => "index.md",
        "API Documentation" => "api.md",
        "Library" => "library.md"
    ]
)

deploydocs(
    deps = nothing,
    make = nothing,
    target = "build",
    repo = "github.com/kaarthiksundar/Dubins.jl.git",
    julia = "1.1"
)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct ProvidersApi <: SwaggerApi
    client::Swagger.Client
end

"""

Gets the specified resource provider.
Param: resourceProviderNamespace::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: expand::String
Return: Provider
"""
function _swaggerinternal_providersGet(_api::ProvidersApi, resourceProviderNamespace::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Provider, "/subscriptions/{subscriptionId}/providers/{resourceProviderNamespace}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceProviderNamespace", resourceProviderNamespace)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function providersGet(_api::ProvidersApi, resourceProviderNamespace::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_providersGet(_api, resourceProviderNamespace, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function providersGet(_api::ProvidersApi, response_stream::Channel, resourceProviderNamespace::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_providersGet(_api, resourceProviderNamespace, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets the specified resource provider at the tenant level.
Param: resourceProviderNamespace::String (required)
Param: api_version::String (required)
Param: expand::String
Return: Provider
"""
function _swaggerinternal_providersGetAtTenantScope(_api::ProvidersApi, resourceProviderNamespace::String, api_version::String; expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Provider, "/providers/{resourceProviderNamespace}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceProviderNamespace", resourceProviderNamespace)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function providersGetAtTenantScope(_api::ProvidersApi, resourceProviderNamespace::String, api_version::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_providersGetAtTenantScope(_api, resourceProviderNamespace, api_version; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function providersGetAtTenantScope(_api::ProvidersApi, response_stream::Channel, resourceProviderNamespace::String, api_version::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_providersGetAtTenantScope(_api, resourceProviderNamespace, api_version; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets all resource providers for a subscription.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: top::Int32
Param: expand::String
Return: ProviderListResult
"""
function _swaggerinternal_providersList(_api::ProvidersApi, api_version::String, subscriptionId::String; top=nothing, expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ProviderListResult, "/subscriptions/{subscriptionId}/providers", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$top", top)  # type Int32
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function providersList(_api::ProvidersApi, api_version::String, subscriptionId::String; top=nothing, expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_providersList(_api, api_version, subscriptionId; top=top, expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function providersList(_api::ProvidersApi, response_stream::Channel, api_version::String, subscriptionId::String; top=nothing, expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_providersList(_api, api_version, subscriptionId; top=top, expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets all resource providers for the tenant.
Param: api_version::String (required)
Param: top::Int32
Param: expand::String
Return: ProviderListResult
"""
function _swaggerinternal_providersListAtTenantScope(_api::ProvidersApi, api_version::String; top=nothing, expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ProviderListResult, "/providers", ["azure_auth"])
    Swagger.set_param(_ctx.query, "\$top", top)  # type Int32
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function providersListAtTenantScope(_api::ProvidersApi, api_version::String; top=nothing, expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_providersListAtTenantScope(_api, api_version; top=top, expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function providersListAtTenantScope(_api::ProvidersApi, response_stream::Channel, api_version::String; top=nothing, expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_providersListAtTenantScope(_api, api_version; top=top, expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Registers a subscription with a resource provider.
Param: resourceProviderNamespace::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Provider
"""
function _swaggerinternal_providersRegister(_api::ProvidersApi, resourceProviderNamespace::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Provider, "/subscriptions/{subscriptionId}/providers/{resourceProviderNamespace}/register", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceProviderNamespace", resourceProviderNamespace)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function providersRegister(_api::ProvidersApi, resourceProviderNamespace::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_providersRegister(_api, resourceProviderNamespace, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function providersRegister(_api::ProvidersApi, response_stream::Channel, resourceProviderNamespace::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_providersRegister(_api, resourceProviderNamespace, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Unregisters a subscription from a resource provider.
Param: resourceProviderNamespace::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Provider
"""
function _swaggerinternal_providersUnregister(_api::ProvidersApi, resourceProviderNamespace::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Provider, "/subscriptions/{subscriptionId}/providers/{resourceProviderNamespace}/unregister", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceProviderNamespace", resourceProviderNamespace)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function providersUnregister(_api::ProvidersApi, resourceProviderNamespace::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_providersUnregister(_api, resourceProviderNamespace, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function providersUnregister(_api::ProvidersApi, response_stream::Channel, resourceProviderNamespace::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_providersUnregister(_api, resourceProviderNamespace, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export providersGet, providersGetAtTenantScope, providersList, providersListAtTenantScope, providersRegister, providersUnregister


text:
using LinearAlgebra:dot

"""
    struct LinearV <: AbstractVApproximator{Int}
        features::Array{Float64, 2}
        weights::Vector{Float64}
    end

Using a matrix `features` to represent each state along with a vector of `weights`.

See more details at Section (9.4) on Page 205 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
struct LinearV <: AbstractVApproximator{Int}
    features::Array{Float64, 2}
    weights::Vector{Float64}
end

(learner::LinearV)(s::Int) = @views dot(learner.features[s, :], learner.weights)

function update!(learner::LinearV, s::Int, e::Float64)
    for i in 1:length(learner.weights)
        learner.weights[i] += learner.features[s, i] * e
    end
end


text:
## lattice.jl : implementation of rank-1 lattice rules
#
# Basic implementation of rank-1 lattice rules to generate points in [0,1]
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

abstract type AbstractLatticeRule{s} <: AbstractRNG end

## LatticeRule32 ##
mutable struct LatticeRule32{s} <: AbstractLatticeRule{s}
    z::Vector{UInt32}
    n::UInt32
end

"""
    LatticeRule32(z, s, n)
    LatticeRule32(z, s)
    LatticeRule32(z)

Returns a rank-1 lattice rule in `s` dimensions with generating vector `z` and at most `n` points.

When no maximum number of points `n` is provided, we assume `n=2^32`. When no number of dimensions `s` is provided, we assume `s=length(z)`.

!!! info

    Technically, we return an extensible lattice sequence where the `k`-th point is transformed using the radical inverse function. This has the advantage that we can add points to the set without changing the already computed points.

More generating vectors can be found online [here](https://web.maths.unsw.edu.au/~fkuo/lattice/index.html) or [here](https://people.cs.kuleuven.be/~dirk.nuyens/qmc-generators/).

# Examples
```jldoctest; setup = :(using MultilevelEstimators; import Random; Random.seed!(1))
julia> lattice_rule = LatticeRule32([0x00000001,0x00000022], 2, 0x00000037) # Fibonacci lattice rule
LatticeRule32{2}

julia> get_point(lattice_rule, 2)
2-element Array{Float64,1}:
 0.75
 0.5

julia> get_point(lattice_rule, 56) # returns random points beyond n
2-element Array{Float64,1}:
 0.23603334566204692
 0.34651701419196046
```
See also: [`get_point`](@ref), [`ShiftedLatticeRule`](@ref)
"""
LatticeRule32(z::Vector{UInt32}) = LatticeRule32(z, length(z))
LatticeRule32(z::Vector{UInt32}, s::Integer) = LatticeRule32(z, s, typemax(UInt32))
LatticeRule32(z::Vector{UInt32}, s::Integer, n::Integer) = check_args(z, s, n) && LatticeRule32{s}(view(z, 1:s), convert(UInt32, n))

"""
    LatticeRule32(str, s, n)
    LatticeRule32(str, s)
    LatticeRule32(str)

Returns a rank-1 lattice rule in `s` dimensions with generating vector `z` read from the file `str` and at most `n` points.

# Examples
```jldoctest; setup = :(using MultilevelEstimators)
julia> z_file = joinpath(dirname(pathof(MultilevelEstimators)), "generating_vectors", "K_3600_32.txt");

julia> lattice_rule = LatticeRule32(z_file, 16)
LatticeRule32{16}

julia> get_point(lattice_rule, 123)
16-element Array{Float64,1}:
 0.3828125
 0.2890625
 0.1484375
 0.3515625
 0.6015625
 0.1171875
 0.4296875
 0.4609375
 0.2265625
 0.6484375
 0.9609375
 0.6796875
 0.0546875
 0.2265625
 0.1328125
 0.1640625

```
See also: [`get_point`](@ref), [`ShiftedLatticeRule`](@ref)
"""
LatticeRule32(s::Int) = check_arg(s) && LatticeRule32(joinpath(@__DIR__(), "generating_vectors", "K_3600_32.txt"), s)
LatticeRule32(str::String, s::Integer) = LatticeRule32(vec(readdlm(str, UInt32)), s)
LatticeRule32(str::String, s::Integer, n::Integer) = LatticeRule32(vec(readdlm(str, UInt32)), s, n)

show(io::IO, lattice_rule::LatticeRule32{s}) where s = print(io, string("LatticeRule32{", s, "}"))

## ShiftedLatticeRule ##
"""
    ShiftedLatticeRule(lat::AbstractLatticeRule)

Returns a shifted rank-1 lattice rule based on the lattice rule `lat`.

# Examples
```jldoctest; setup = :(using MultilevelEstimators; import Random; Random.seed!(1))
julia> z_file = joinpath(dirname(pathof(MultilevelEstimators)), "generating_vectors", "K_3600_32.txt");

julia> lattice_rule = LatticeRule32(z_file, 16)
LatticeRule32{16}

julia> shifted_lattice_rule = ShiftedLatticeRule(lattice_rule)
ShiftedLatticeRule{LatticeRule32{16}}

julia> get_point(shifted_lattice_rule, 0)
16-element Array{Float64,1}:
 0.23603334566204692
 0.34651701419196046
 0.3127069683360675
 0.00790928339056074
 0.4886128300795012
 0.21096820215853596
 0.951916339835734
 0.9999046588986136
 0.25166218303197185
 0.9866663668987996
 0.5557510873245723
 0.43710797460962514
 0.42471785049513144
 0.773223048457377
 0.2811902322857298
 0.20947237319807077

```
See also: [`LatticeRule32`](@ref), [`get_point`](@ref)
"""
struct ShiftedLatticeRule{L<:AbstractLatticeRule, V<:AbstractVector{<:Real}}
    lattice_rule::L
    shift::V
end

ShiftedLatticeRule(lattice_rule::AbstractLatticeRule{s}) where s = ShiftedLatticeRule(lattice_rule, rand(s))

show(io::IO, s::ShiftedLatticeRule) = print(io, string("ShiftedLatticeRule{", s.lattice_rule, "}"))

## get_point ##
@inline function get_point!(x::Vector{<:AbstractFloat}, lattice_rule::LatticeRule32, k::UInt32)
    ϕ = convert(UInt32, reversebits(xor(k, k >> 1)))
	copyto!(x, ϕ * 2.0^(-32) * lattice_rule.z .% 1)
end

"""
    get_point(lat::AbstractLatticeRule, k::Integer)

Get the `k`-th point of the lattice rule `lat`.

```jldoctest; setup = :(using MultilevelEstimators)
julia> lattice_rule = LatticeRule32([0x00000001,0x00000022], 2, 0x00000037) # Fibonacci lattice rule
LatticeRule32{2}

julia> get_point(lattice_rule, 2)
2-element Array{Float64,1}:
 0.75
 0.5

```
See also: [`LatticeRule32`](@ref), [`ShiftedLatticeRule`](@ref)
"""
@inline get_point(lattice_rule::LatticeRule32{s}, k::UInt32) where s = k > lattice_rule.n ? rand(Float32, s) : get_point!(Vector{Float32}(undef, s), lattice_rule, k)

@inline get_point(lattice_rule::LatticeRule32{s}, k::Int) where s = k > lattice_rule.n ? rand(s) : get_point!(Vector{Float64}(undef, s), lattice_rule, convert(UInt32, k))

@inline function get_point!(x::Vector{<:AbstractFloat}, shifted_lattice_rule::ShiftedLatticeRule, k::UInt32)
    ϕ = convert(UInt32, reversebits(xor(k, k >> 1)))
	copyto!(x, (ϕ * 2.0^(-32) * shifted_lattice_rule.lattice_rule.z .+ shifted_lattice_rule.shift) .% 1)
end

@inline get_point(shifted_lattice_rule::ShiftedLatticeRule{<:LatticeRule32{s}}, k::UInt32) where s = k > shifted_lattice_rule.lattice_rule.n ? rand(Float32, s) : get_point!(Vector{Float32}(undef, s), shifted_lattice_rule, k)

@inline get_point(shifted_lattice_rule::ShiftedLatticeRule{<:LatticeRule32{s}}, k::Int) where s = k > shifted_lattice_rule.lattice_rule.n ? rand(s) : get_point!(Vector{Float64}(undef, s), shifted_lattice_rule, convert(UInt32, k))

## reverse_bits ##
reversebits(n::U) where U<:Unsigned = parse(U, reverse(bitstring(n)), base=2)

## input checking ##
function check_args(z::Vector{UInt32}, s::Integer, n::Integer)
	check_larger_than(LatticeRule32, s, "number of dimensions s", 0)
	check_ordered(LatticeRule32, s, length(z) + 1, "number of dimensions s", "or equal to the length of the generating vector z")
	check_ordered(LatticeRule32, n, typemax(UInt32) + 1, "maximum number of points n", "or equal to 2^32, consider implementing a LatticeRule64 type")
end

check_arg(s::Integer) = check_ordered(LatticeRule32, s, 3601, "number of dimensions s", "or equal to 3600, please supply your own generating vector z to proceed")


text:
using SafeTestsets
@safetestset "My form tests" begin include("my_form_test.jl") end


text:
using OpenSSLCrypto.Crypto
using Base.Test
using Codecs

s = Crypto.hmacsha256_digest(
    "GET\nelasticmapreduce.amazonaws.com\n/\nAWSAccessKeyId=AKIAIOSFODNN7EXAMPLE&Action=DescribeJobFlows&SignatureMethod=HmacSHA256&SignatureVersion=2&Timestamp=2011-10-03T15%3A19%3A30&Version=2009-03-31",
    "wJalrXUtnFEMI/K7MDENG/bPxRfiCYEXAMPLEKEY")

sb64 = bytestring(Codecs.encode(Base64, s))

#println("Output bytes : $sl, b64 : $sb64")

@test "i91nKc4PWAt0JJIdXwz9HxZCJDdiy6cf/Mj6vPxyYIs=" == sb64
println("PASSED : HMAC-SHA256")

s = Crypto.hmacsha1_digest(
    "GET\n\n\nTue, 27 Mar 2007 19:36:42 +0000\n/johnsmith/photos/puppy.jpg",
    "wJalrXUtnFEMI/K7MDENG/bPxRfiCYEXAMPLEKEY")

sb64 = bytestring(Codecs.encode(Base64, s))

#println("Output bytes : $sl, b64 : $sb64")

@test "bWq2s1WEIj+Ydj0vQ697zp+IXMU=" == sb64
println("PASSED : HMAC-SHA1")

d="Just a really long string......Just a really long string......Just a really long string......Just a really long string......Just a really long string......Just a really long string......"
sb = bytes2hex(Crypto.md5(d))

@test "5c0614d60fb48f8f318c3578ca959120" == sb
println("PASSED : MD5")

repeat="hello"
xtimes=100000
real_long_string = ^(repeat, xtimes)

iob = IOBuffer(length(repeat) * xtimes)
write(iob, real_long_string)

seekstart(iob)

iomd5 = bytes2hex(Crypto.md5(iob))

smd5 = bytes2hex(Crypto.md5(real_long_string))

@test iomd5 == smd5
println("PASSED : IO MD5")


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the associated Legendre functions.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Holmes, S. A. and W. E. Featherstone, 2002. A unified approach to the
#       Clenshaw summation and the recursive computation of very high degree and
#       order normalised associated Legendre functions. Journal of Geodesy,
#       76(5), pp. 279-299.
#
#       For more info.: http://mitgcm.org/~mlosch/geoidcookbook/node11.html
#
#   [2] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
#   [3] Schmidt, A (1917). Erdmagnetismus, Enzykl. Math. Wiss., 6, pp. 265–396.
#
#   [4] Winch, D. E., Ivers, D. J., Turner, J. P. R., Stening R. J (2005).
#       Geomagnetism and Schmidt quasi-normalization. Geophysical Journal
#       International, 160(2), pp. 487-504.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export denormalize_legendre, legendre!, legendre

export legendre_fully_normalized!, legendre_fully_normalized
export legendre_schmidt_quasi_normalized!, legendre_schmidt_quasi_normalized
export legendre_conventional!, legendre_conventional

"""
    legendre!([N,] P::AbstractMatrix, ϕ::Number, ph_term::Bool = false, n_max::Integer = -1, m_max::Integer = -1)

Compute the associated Legendre function `P_n,m[cos(ϕ)]`. The maximum degree and
order that will be computed are given by the parameters `n_max` and `m_max`. If
they are negative, then the dimensions of matrix `P` will be used.

The result will be stored at matrix `P`.

The optional parameter `N` can be used to select the normalization. The
following values are valid:

- `Val(:full)`: Compute the fully normalized associated Legendre function (see
    [`legendre_fully_normalized!`](@ref)).
- `Val(:schmidt)`: Compute the Schmidt quasi-normalized associated Legendre
    function (see [`legendre_schmidt_quasi_normalized!`](@ref)).
- `Val(:conv)`: Compute the conventional associated Legendre function (see
    [`legendre_conventional!`](@ref)).

If `N` is omitted, then the full normalization will be used.

If `ph_term` is set to `true`, then the Condon-Shortley phase term `(-1)ᵐ` will
be included. If `ph_term` is not present, then it defaults to `false`.
"""
function legendre!(
    P::AbstractMatrix,
    ϕ::Number,
    ph_term::Bool = false,
    n_max::Integer = -1,
    m_max::Integer = -1
)
    return legendre_fully_normalized!(P, float(ϕ), ph_term, n_max, m_max)
end

function legendre!(
    ::Val{:full},
    P::AbstractMatrix,
    ϕ::Number,
    ph_term::Bool = false,
    n_max::Integer = -1,
    m_max::Integer = -1
)
    return legendre_fully_normalized!(P, float(ϕ), ph_term, n_max, m_max)
end

function legendre!(
    ::Val{:schmidt},
    P::AbstractMatrix,
    ϕ::Number,
    ph_term::Bool = false,
    n_max::Integer = -1,
    m_max::Integer = -1
)
    return legendre_schmidt_quasi_normalized!(P, float(ϕ), ph_term, n_max, m_max)
end

function legendre!(
    ::Val{:conv},
    P::AbstractMatrix,
    ϕ::Number,
    ph_term::Bool = false,
    n_max::Integer = -1,
    m_max::Integer = -1
)
    return legendre_conventional!(P, float(ϕ), ph_term, n_max, m_max)
end

"""
    legendre([N,] ϕ::Number, n_max::Integer, m_max::Integer = -1, ph_term::Bool = false)

Compute the associated Legendre function `P_n,m[cos(ϕ)]`. The maximum degree
that will be computed is `n_max` and the maximum order is `m_max`. Notice that
if `m_max` is higher than `n_max` or negative, than it is set to `n_max`.

The optional parameter `N` can be used to select the normalization. The
following values are valid:

- `Val(:full)`: Compute the fully normalized associated Legendre function (see
    [`legendre_fully_normalized`](@ref)).
- `Val(:schmidt)`: Compute the Schmidt quasi-normalized associated Legendre
    function (see [`legendre_schmidt_quasi_normalized`](@ref)).
- `Val(:conv)`: Compute the conventional associated Legendre function (see
    [`legendre_conventional`](@ref)).

If `N` is omitted, then the full normalization will be used (`Val(:full)`).

If `ph_term` is set to `true`, then the Condon-Shortley phase term `(-1)ᵐ` will
be included. If `ph_term` is not present, then it defaults to `false`.

# Returns

A matrix with the Legendre associated functions `P_n,m[cos(ϕ)]`.
"""
function legendre(
    ϕ::Number,
    n_max::Integer,
    m_max::Integer = -1,
    ph_term::Bool = false
)
    return legendre_fully_normalized(float(ϕ), n_max, m_max, ph_term)
end

function legendre(
    ::Val{:full},
    ϕ::Number,
    n_max::Integer,
    m_max::Integer = -1,
    ph_term::Bool = false
)
    return legendre_fully_normalized(float(ϕ), n_max, m_max, ph_term)
end

function legendre(
    ::Val{:schmidt},
    ϕ::Number,
    n_max::Integer,
    m_max::Integer = -1,
    ph_term::Bool = false
)
    return legendre_schmidt_quasi_normalized(float(ϕ), n_max, m_max, ph_term)
end

function legendre(
    ::Val{:conv},
    ϕ::Number,
    n_max::Integer,
    m_max::Integer = -1,
    ph_term::Bool = false
)
    return legendre_conventional(float(ϕ), n_max, m_max, ph_term)
end

################################################################################
#                Fully Normalized Associated Legendre Functions
################################################################################

"""
    legendre_fully_normalized!(P::AbstractMatrix, ϕ::Number, ph_term::Bool = false, n_max::Integer = -1, m_max::Integer = -1)

Compute the fully normalized associated Legendre function `P_n,m[cos(ϕ)]`. The
maximum degree and order that will be computed are given by the parameters
`n_max` and `m_max`. If they are negative, then the dimensions of matrix `P`
will be used:

    maximum degree -> number of rows
    maximum order  -> number of columns

The result will be stored at matrix `P`.

If `ph_term` is set to `true`, then the Condon-Shortley phase term `(-1)ᵐ` will
be included. If `ph_term` is not present, then it defaults to `false`.

# Remarks

This algorithm was based on **[1]**. Our definition of fully normalized associated
Legendre function can be seen in **[2, p. 546]**. The conversion is obtained by:

                 _                     -
                |  (n-m)! . k . (2n+1)  |      k = 1 if m  = 0
    K_n,m = sqrt| --------------------- |,     k = 2 if m != 0
                |         (n+m)!        |
                 -                     -
    _
    P_n,m = P_n,m * K_n,m,

          _
    where P_n,m is the fully normalized Legendre associated function.

# References

- **[1]** Holmes, S. A. and W. E. Featherstone, 2002. A unified approach to the
    Clenshaw summation and the recursive computation of very high degree and
    order normalised associated Legendre functions. Journal of Geodesy,
    76(5), pp. 279-299. For more info.:
    http://mitgcm.org/~mlosch/geoidcookbook/node11.html

- **[2]** Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
function legendre_fully_normalized!(
    P::AbstractMatrix{T},
    ϕ::Number,
    ph_term::Bool = false,
    n_max::Integer = -1,
    m_max::Integer = -1
) where T<:AbstractFloat

    # Obtain the maximum degree and order that must be computed.
    n_max, m_max = _get_degree_and_order(P, n_max, m_max)

    # Auxiliary variables to improve code performance.
    s, c = sincos(T(ϕ))

    # The sine must be always positive. In fact, `s` was previously computed
    # using `sqrt(1 - c^2)`. However, we had numerical problems for very small
    # angles that lead to `cos(ϕ) = 1`.
    if (s < 0)
        s = -s
    end

    s_fact = !ph_term ? +s : -s

    @inbounds for n in 0:n_max
        # Starting values.
        if n == 0
            P[0+1,0+1] = 1
            continue
        elseif n == 1
            P[1+1,0+1] = +sqrt(T(3)) * c

            if m_max > 0
                P[1+1,1+1] = +sqrt(T(3)) * s_fact
            end

            continue
        end

        aux_n = T((2n - 1) * (2n + 1))

        for m in 0:n

            if n == m
                P[n+1,n+1] = s_fact * sqrt(T(2n + 1) / T(2n)) * P[n-1+1,n-1+1]
            else
                aux_nm = T(n - m) * T(n + m)
                a_nm   = sqrt(aux_n / aux_nm)
                b_nm   = sqrt(
                    T((2n + 1) * (n + m - 1) * (n - m - 1)) /
                    T(aux_nm * (2n - 3))
                )

                # We assume that the matrix is not initialized. Hence, we must
                # not access elements on the upper triangle.
                if m != n - 1
                    P[n+1,m+1] = a_nm * c * P[n-1+1,m+1] - b_nm * P[n-2+1,m+1]
                else
                    P[n+1,m+1] = a_nm * c * P[n-1+1,m+1]
                end
            end

            # Check if the maximum desired order has been reached.
            m == m_max && break
        end
    end

    return nothing
end

"""
    legendre_fully_normalized(ϕ::T, n_max::Integer, m_max::Integer = -1, ph_term::Bool = false) where T<:AbstractFloat

Compute the fully normalized associated Legendre function `P_n,m[cos(ϕ)]`. The
maximum degree that will be computed is `n_max` and the maximum order is
`m_max`. Notice that if `m_max` is higher than `n_max` or negative, than it is
set to `n_max`.

If `ph_term` is set to `true`, then the Condon-Shortley phase term `(-1)ᵐ` will
be included. If `ph_term` is not present, then it defaults to `false`.

# Returns

A matrix with the Legendre associated functions `P_n,m[cos(ϕ)]`.

# Remarks

This algorithm was based on [1]. Our definition of fully normalized associated
Legendre function can be seen in [2, p. 546]. The conversion is obtained by:

                 _                     -
                |  (n-m)! . k . (2n+1)  |      k = 1 if m  = 0
    K_n,m = sqrt| --------------------- |,     k = 2 if m != 0
                |         (n+m)!        |
                 -                     -
    _
    P_n,m = P_n,m * K_n,m,

          _
    where P_n,m is the fully normalized Legendre associated function.

# References

- **[1]** Holmes, S. A. and W. E. Featherstone, 2002. A unified approach to the
    Clenshaw summation and the recursive computation of very high degree and
    order normalised associated Legendre functions. Journal of Geodesy,
    76(5), pp. 279-299. For more info.:
    http://mitgcm.org/~mlosch/geoidcookbook/node11.html

- **[2]** Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
function legendre_fully_normalized(
    ϕ::T,
    n_max::Integer,
    m_max::Integer = -1,
    ph_term::Bool = false
) where T<:AbstractFloat
    (n_max < 0) && throw(ArgumentError("n_max must be positive."))

    if ((m_max < 0) || (m_max > n_max))
        m_max = n_max
    end

    P = zeros(T, n_max + 1, m_max + 1)
    legendre_fully_normalized!(P, ϕ, ph_term)

    return P
end

################################################################################
#            Schmidt Quasi-Normalized Associated Legendre Functions
################################################################################

"""
    legendre_schmidt_quasi_normalized!(P::AbstractMatrix, ϕ::Number, ph_term::Bool = false, n_max::Integer = -1, m_max::Integer = -1)

Compute the Schmidt quasi-normalized associated Legendre function
`P_n,m[cos(ϕ)]` [3,4]. The maximum degree and order that will be computed are
given by the parameters `n_max` and `m_max`. If they are negative, then the
dimensions of matrix `P` will be used:

    maximum degree -> number of rows
    maximum order  -> number of columns

The result will be stored at matrix `P`.

If `ph_term` is set to `true`, then the Condon-Shortley phase term `(-1)ᵐ` will
be included. If `ph_term` is not present, then it defaults to `false`.

# Remarks

This algorithm was based on **[1, 2]**. The conversion is obtained by:

                 _           -
                |     (n-m)!  |    k = 1 if m  = 0
    K_n,m = sqrt| k. -------- |,   k = 2 if m != 0
                |     (n+m)!  |
                 -           -

    =
    P_n,m = P_n,m * K_n,m,

          =
    where P_n,m is the quasi-normalized normalized Legendre associated function.

# References

- **[1]** Schmidt, A (1917). Erdmagnetismus, Enzykl. Math. Wiss., 6, pp.
    265–396.

- **[2]** Winch, D. E., Ivers, D. J., Turner, J. P. R., Stening R. J (2005).
    Geomagnetism and Schmidt quasi-normalization. Geophysical Journal
    International, 160(2), pp. 487-504.
"""
function legendre_schmidt_quasi_normalized!(
    P::AbstractMatrix{T},
    ϕ::Number,
    ph_term::Bool = false,
    n_max::Integer = -1,
    m_max::Integer = -1
) where T<:AbstractFloat

    # Obtain the maximum degree and order that must be computed.
    n_max, m_max = _get_degree_and_order(P, n_max, m_max)

    # Auxiliary variables to improve code performance.
    s, c = sincos(T(ϕ))

    # The sine must be always positive. In fact, `s` was previously computed
    # using `sqrt(1 - c^2)`. However, we had numerical problems for very small
    # angles that lead to `cos(ϕ) = 1`.
    if (s < 0)
        s = -s
    end

    s_fact = !ph_term ? +s : -s

    @inbounds for n in 0:n_max
        # Starting values.
        if n == 0
            P[0+1,0+1] = 1
            continue

        elseif n == 1
            P[1+1,0+1] = +c

            if m_max > 0
                P[1+1,1+1] = +s_fact
            end

            continue
        end

        aux_n = T(2n - 1) # -> sqrt( (2n-1)*(2n-1) )

        for m in 0:n

            if m == n
                P[n+1,n+1] = s_fact * sqrt(aux_n / T(2n)) * P[n-1+1,n-1+1]
            else
                aux_nm = sqrt(T((n - m) * (n + m)))
                a_nm   = aux_n / aux_nm
                b_nm   = sqrt(T((n + m - 1) * (n - m - 1))) / aux_nm

                # We assume that the matrix is not initialized. Hence, we must not
                # access elements on the upper triangle.
                if m != n - 1
                    P[n+1,m+1] = a_nm * c * P[n-1+1,m+1] - b_nm * P[n-2+1,m+1]
                else
                    P[n+1,m+1] = a_nm * c * P[n-1+1,m+1]
                end
            end

            # Check if the maximum desired order has been reached.
            m == m_max && break
        end
    end

    return nothing
end

"""
    legendre_schmidt_quasi_normalized(ϕ::T, n_max::Integer, m_max::Integer = -1, ph_term::Bool = false) where T<:AbstractFloat

Compute the Schmidt quasi-normalized associated Legendre function
`P_n,m[cos(ϕ)]`. The maximum degree that will be computed is `n_max` and the
maximum order is `m_max`. Notice that if `m_max` is higher than `n_max` or
negative, than it is set to `n_max`.

If `ph_term` is set to `true`, then the Condon-Shortley phase term `(-1)ᵐ` will
be included. If `ph_term` is not present, then it defaults to `false`.

# Returns

A matrix with the Legendre associated functions `P_n,m[cos(ϕ)]`.

# Remarks

This algorithm was based on **[1, 2]**. The conversion is obtained by:

                 _           -
                |     (n-m)!  |    k = 1 if m  = 0
    K_n,m = sqrt| k. -------- |,   k = 2 if m != 0
                |     (n+m)!  |
                 -           -

    =
    P_n,m = P_n,m * K_n,m,

          =
    where P_n,m is the quasi-normalized normalized Legendre associated function.

# References

- **[1]** Schmidt, A (1917). Erdmagnetismus, Enzykl. Math. Wiss., 6, pp.
    265–396.

- **[2]** Winch, D. E., Ivers, D. J., Turner, J. P. R., Stening R. J (2005).
    Geomagnetism and Schmidt quasi-normalization. Geophysical Journal
    International, 160(2), pp. 487-504.
"""
function legendre_schmidt_quasi_normalized(
    ϕ::T,
    n_max::Integer,
    m_max::Integer = -1,
    ph_term::Bool = false
) where T<:AbstractFloat
    (n_max < 0) && throw(ArgumentError("n_max must be positive."))

    if ((m_max < 0) || (m_max > n_max))
        m_max = n_max
    end

    P = zeros(T, n_max + 1, m_max + 1)
    legendre_schmidt_quasi_normalized!(P, ϕ, ph_term)

    return P
end

################################################################################
#                  Conventional Associated Legendre Function
################################################################################

"""
    legendre_conventional!(P::AbstractMatrix, ϕ::Number, ph_term::Bool = false, n_max::Integer = -1, m_max::Integer = -1)

Compute the conventional associated Legendre function `P_n,m[cos(ϕ)]`. The
maximum degree and order that will be computed are given by the parameters
`n_max` and `m_max`. If they are negative, then the dimensions of matrix `P`
will be used:

    maximum degree -> number of rows
    maximum order  -> number of columns

The result will be stored at matrix `P`.

If `ph_term` is set to `true`, then the Condon-Shortley phase term `(-1)ᵐ` will
be included. If `ph_term` is not present, then it defaults to `false`.
"""
function legendre_conventional!(
    P::AbstractMatrix{T},
    ϕ::Number,
    ph_term::Bool = false,
    n_max::Integer = -1,
    m_max::Integer = -1
) where T<:AbstractFloat
    # Obtain the maximum degree and order that must be computed.
    n_max, m_max = _get_degree_and_order(P, n_max, m_max)

    # Auxiliary variables to improve code performance.
    s, c = sincos(T(ϕ))

    # The sine must be always positive. In fact, `s` was previously computed
    # using `sqrt(1 - c^2)`. However, we had numerical problems for very small
    # angles that lead to `cos(ϕ) = 1`.
    if (s < 0)
        s = -s
    end

    s_fact = !ph_term ? +s : -s

    @inbounds for n in 0:n_max
        # Starting values.
        if n == 0
            P[0+1,0+1] = 1
            continue
        elseif n == 1
            P[1+1,0+1] = +c

            if m_max > 0
                P[1+1,1+1] = +s_fact
            end

            continue
        end

        aux_n = T(2n - 1) # -> sqrt( (2n-1)*(2n-1) )

        for m in 0:n

            if n == m
                P[n+1,n+1] = s_fact * aux_n * P[n-1+1,n-1+1]
            else
                aux_nm = T(n - m) # -> sqrt( (n-m)*(n-m) )
                a_nm   = aux_n / aux_nm
                b_nm   = T(n + m - 1) / aux_nm  # -> sqrt( (n+m-1)*(n+m-1) ) / aux_nm

                # We assume that the matrix is not initialized. Hence, we must
                # not access elements on the upper triangle.
                if m != n-1
                    P[n+1,m+1] = a_nm * c * P[n-1+1,m+1] - b_nm * P[n-2+1,m+1]
                else
                    P[n+1,m+1] = a_nm * c *P[n-1+1,m+1]
                end
            end

            # Check if the maximum desired order has been reached.
            m == m_max && break
        end
    end

    return nothing
end

"""
    legendre_conventional(ϕ::T, n_max::Integer, m_max::Integer = -1, ph_term::Bool = false) where T<:AbstractFloat

Compute the conventional associated Legendre function `P_n,m[cos(ϕ)]`. The
maximum degree that will be computed is `n_max` and the maximum order is
`m_max`. Notice that if `m_max` is higher than `n_max` or negative, than it is
set to `n_max`.

If `ph_term` is set to `true`, then the Condon-Shortley phase term `(-1)ᵐ` will
be included. If `ph_term` is not present, then it defaults to `false`.

# Returns

A matrix with the Legendre associated functions `P_n,m[cos(ϕ)]`.
"""
function legendre_conventional(
    ϕ::T,
    n_max::Integer,
    m_max::Integer = -1,
    ph_term::Bool = false
) where T<:AbstractFloat
    (n_max < 0) && throw(ArgumentError("n_max must be positive."))

    if ((m_max < 0) || (m_max > n_max))
        m_max = n_max
    end

    P = zeros(T, n_max + 1, m_max + 1)
    legendre_conventional!(P, ϕ, ph_term)

    return P
end

################################################################################
#                                   Private
################################################################################

"""
    _get_degree_and_order(P, n_max, m_max)

Return the maximum degree and order to compute the Legendre associated functions
given the matrix `P` and the configuration values `n_max` and `m_max`.
"""
@inline function _get_degree_and_order(P, n_max, m_max)
    # Get the size of the matrix.
    rows, cols = size(P)

    # If the order or degree is less than 0, then the user wants to use all the
    # available memory.
    if n_max < 0
        n_max = rows - 1
    end

    if m_max < 0
        m_max = (cols <= rows) ? cols - 1 : n_max
    end

    # Make sure that the degree and order fits the matrix.
    if n_max > rows - 1
        n_max = rows - 1
    end

    if ((m_max > cols - 1) || (m_max > n_max))
        m_max = min(cols - 1, n_max)
    end

    return n_max, m_max
end


text:
include("../src/geom_module.jl")

module ptest

using Random
using BenchmarkTools
using ..geom
using PyPlot
PyPlot.pygui(true)

# poly_tri = geom.create_polygon([0., 10., 5.], [0., 0., 5])
poly_L_shape = geom.create_polygon(
    geom.MutablePoint,
    [0., 10, 10, 5, 5, 0],
    [0., 0, 5, 5, 15, 15])

poly_3 = geom.create_polygon(
    geom.MutablePoint,
    [0., 1, 0],
    [0., 1, 2])

poly_4_convex = geom.create_polygon(
    geom.MutablePoint,
    [0., 1, 0, -1],
    [0., 1, 2, 1])

poly_4_concave = geom.create_polygon(
    geom.MutablePoint,
    [-1, 0., 1, 0],
    [1, 0., 1, 0.7])

poly_5_convex = geom.create_polygon(
    geom.MutablePoint,
    [0., 1, 1.71, 0.5, -0.2],
    [0., 0, 0.5, 1.21, 0.5])

poly_5_concave = geom.create_polygon(
    geom.MutablePoint,
    [0., 1, 0.21, 0.5, -0.2],
    [0., 0, 0.5, 1.21, 0.5])

# PyPlot.figure()
# PyPlot.grid()
# PyPlot.plot(
#     geom.x(poly_5_concave.point_array),
#     geom.y(poly_5_concave.point_array),
#     linewidth = 2,
#     marker = "o")
#
PyPlot.figure()
PyPlot.grid()
PyPlot.plot(
    geom.x(poly_5_convex.point_array),
    geom.y(poly_5_convex.point_array),
    linewidth = 2,
    marker = "o")

tri_array = geom.triangulate(poly_5_convex)

# L = length(poly_L_shape.point_array)
# b = falses(L, L)
#
# for ii = 1:L
#     for jj = 1:L
#         b[ii, jj] = geom.is_diagonal_helper(ptest.poly_L_shape, ii, jj)
#     end
# end

# geom.is_diagonal_helper(ptest.poly_L_shape, 2, 5)
# geom.is_diagonal_helper(ptest.poly_L_shape, 2, 6)
# geom.is_diagonal_helper(ptest.poly_L_shape, 3, 6)

#=
dmesh = geom.build_delaunay_mesh!(poly_L_shape.point_array)

# (ptri_x, ptri_y) = geom.discretize_for_plot(poly_tri)
(pL_x, pL_y) = geom.discretize_for_plot(poly_L_shape)

tri_vec = geom.triangle_vector_for_plot(dmesh)
(init_tri_x, init_tri_y) = geom.discretize_for_plot(
  [dmesh.vertex_array[1], dmesh.vertex_array[2], dmesh.vertex_array[3]], true)

# r = geom.rotation_matrix( deg2rad(30.0) )
#
# v_xy = r * [ptest.pL_x ptest.pL_y]'
# pL_x = v_xy[1, :]
# pL_y = v_xy[2, :]

PyPlot.figure()
PyPlot.grid()
PyPlot.plot(pL_x, pL_y, color = "k", alpha = 0.7, marker = ".", markersize = 25, linewidth = 3)
PyPlot.plot(init_tri_x, init_tri_y, marker = "o", linewidth = 2)

for kk = 1:length(tri_vec)
    (p_x, p_y) = geom.discretize_for_plot(tri_vec[kk])
    PyPlot.plot(p_x, p_y, marker = "o", linewidth = 1, color = "r", alpha = 0.7)
end

=#

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
module EnsembleAnalysis

using SciMLBase, Statistics, RecursiveArrayTools

# Getters
get_timestep(sim,i) = (getindex(sol,i) for sol in sim)
get_timepoint(sim,t) = (sol(t) for sol in sim)
function componentwise_vectors_timestep(sim,i)
  arr = [get_timestep(sim,i)...]
  if typeof(arr[1]) <: AbstractArray
   return vecarr_to_vectors(VectorOfArray(arr))
  else
   return arr
  end
 end
 function componentwise_vectors_timepoint(sim,t)
   arr = [get_timepoint(sim,t)...]
   if typeof(arr[1]) <: AbstractArray
    return vecarr_to_vectors(VectorOfArray(arr))
   else
    return arr
   end
  end

# Timestep statistics
timestep_mean(sim,i) = componentwise_mean(get_timestep(sim,i))
timestep_mean(sim,::Colon) = timeseries_steps_mean(sim)
function timestep_median(sim,i)
  arr = componentwise_vectors_timestep(sim,i)
  if typeof(first(arr)) <: AbstractArray
    return reshape([median(x) for x in arr],size(sim[1][i])...)
  else
    return median(arr)
  end
end
timestep_median(sim,::Colon) = timeseries_steps_median(sim)
function timestep_quantile(sim,q,i)
  arr = componentwise_vectors_timestep(sim,i)
  if typeof(first(arr)) <: AbstractArray
    return reshape([quantile(x,q) for x in arr],size(sim[1][i])...)
  else
    return quantile(arr,q)
  end
end
timestep_quantile(sim,q,::Colon) = timeseries_steps_quantile(sim,q)
timestep_meanvar(sim,i) = componentwise_meanvar(get_timestep(sim,i))
timestep_meanvar(sim,::Colon) = timeseries_steps_meanvar(sim)
timestep_meancov(sim,i,j) = componentwise_meancov(get_timestep(sim,i),get_timestep(sim,j))
timestep_meancov(sim,::Colon,::Colon) = timeseries_steps_meancov(sim)
timestep_meancor(sim,i,j) = componentwise_meancor(get_timestep(sim,i),get_timestep(sim,j))
timestep_meancor(sim,::Colon,::Colon) = timeseries_steps_meancor(sim)
timestep_weighted_meancov(sim,W,i,j) = componentwise_weighted_meancov(get_timestep(sim,i),get_timestep(sim,j),W)
timestep_weighted_meancov(sim,W,::Colon,::Colon) = timeseries_steps_weighted_meancov(sim,W)

function timeseries_steps_mean(sim)
  DiffEqArray([timestep_mean(sim,i) for i in 1:length(sim[1])],sim[1].t)
end
function timeseries_steps_median(sim)
  DiffEqArray([timestep_median(sim,i) for i in 1:length(sim[1])],sim[1].t)
end
function timeseries_steps_quantile(sim,q)
  DiffEqArray([timestep_quantile(sim,q,i) for i in 1:length(sim[1])],sim[1].t)
end
function timeseries_steps_meanvar(sim)
  m,v = timestep_meanvar(sim,1)
  means = [m]
  vars = [v]
  for i in 2:length(sim[1])
    m,v = timestep_meanvar(sim,i)
    push!(means,m)
    push!(vars,v)
  end
  DiffEqArray(means,sim[1].t),DiffEqArray(vars,sim[1].t)
end
function timeseries_steps_meancov(sim)
  reshape([timestep_meancov(sim,i,j) for i in 1:length(sim[1]) for j in 1:length(sim[1])],length(sim[1]),length(sim[1]))
end
function timeseries_steps_meancor(sim)
  reshape([timestep_meancor(sim,i,j) for i in 1:length(sim[1]) for j in 1:length(sim[1])],length(sim[1]),length(sim[1]))
end
function timeseries_steps_weighted_meancov(sim,W)
  reshape([timestep_meancov(sim,W,i,j) for i in 1:length(sim[1]) for j in 1:length(sim[1])],length(sim[1]),length(sim[1]))
end

timepoint_mean(sim,t) = componentwise_mean(get_timepoint(sim,t))
function timepoint_median(sim,t)
  arr = componentwise_vectors_timepoint(sim,t)
  if typeof(first(arr)) <: AbstractArray
    return reshape([median(x) for x in arr],size(sim[1][1])...)
  else
    return median(arr)
  end
end
function timepoint_quantile(sim,q,t)
  arr = componentwise_vectors_timepoint(sim,t)
  if typeof(first(arr)) <: AbstractArray
    return reshape([quantile(x,q) for x in arr],size(sim[1][1])...)
  else
    return quantile(arr,q)
  end
end
timepoint_meanvar(sim,t) = componentwise_meanvar(get_timepoint(sim,t))
timepoint_meancov(sim,t1,t2) = componentwise_meancov(get_timepoint(sim,t1),get_timepoint(sim,t2))
timepoint_meancor(sim,t1,t2) = componentwise_meancor(get_timepoint(sim,t1),get_timepoint(sim,t2))
timepoint_weighted_meancov(sim,W,t1,t2) = componentwise_weighted_meancov(get_timepoint(sim,t1),get_timepoint(sim,t2),W)

function EnsembleSummary(sim::SciMLBase.AbstractEnsembleSolution{T,N},
                         t=sim[1].t;quantiles=[0.05,0.95]) where {T,N}
  if typeof(sim[1]) <: DESolution
    m,v = timeseries_point_meanvar(sim,t)
    qlow = timeseries_point_quantile(sim,quantiles[1],t)
    qhigh = timeseries_point_quantile(sim,quantiles[2],t)
  else
    m,v = timeseries_steps_meanvar(sim)
    qlow = timeseries_steps_quantile(sim,quantiles[1])
    qhigh = timeseries_steps_quantile(sim,quantiles[2])
  end

  trajectories = length(sim)
  EnsembleSummary{T,N,typeof(t),typeof(m),typeof(v),typeof(qlow),typeof(qhigh)}(t,m,v,qlow,qhigh,trajectories,sim.elapsedTime,sim.converged)
end

function timeseries_point_mean(sim,ts)
  DiffEqArray([timepoint_mean(sim,t) for t in ts],ts)
end
function timeseries_point_median(sim,ts)
  DiffEqArray([timepoint_median(sim,t) for t in ts],ts)
end
function timeseries_point_quantile(sim,q,ts)
  DiffEqArray([timepoint_quantile(sim,q,t) for t in ts],ts)
end
function timeseries_point_meanvar(sim,ts)
  m,v = timepoint_meanvar(sim,first(ts))
  means = [m]
  vars = [v]
  for t in Iterators.drop(ts,1)
    m,v = timepoint_meanvar(sim,t)
    push!(means,m)
    push!(vars,v)
  end
  DiffEqArray(means,ts),DiffEqArray(vars,ts)
end
function timeseries_point_meancov(sim,ts1,ts2)
  reshape([timepoint_meancov(sim,t1,t2) for t1 in ts1 for t2 in ts2],length(ts1),length(ts2))
end
function timeseries_point_meancor(sim,ts1,ts2)
  reshape([timepoint_meancor(sim,t1,t2) for t1 in ts1 for t2 in ts2],length(ts1),length(ts2))
end
function timeseries_point_weighted_meancov(sim,W,ts1,ts2)
  reshape([timepoint_meancov(sim,W,t1,t2) for t1 in ts1 for t2 in ts2],length(ts1),length(ts2))
end

function componentwise_mean(A)
  x0 = first(A)
  n = 0
  mean = zero(x0)
  for x in A
    n += 1
    if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
      mean .+= x
    else
      mean += x
    end
  end
  if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
    mean ./= n
  else
    mean /= n
  end
  mean
end

# Welford algorithm
# https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
function componentwise_meanvar(A;bessel=true)
  x0 = first(A)
  n = 0
  mean = zero(x0)
  M2 = zero(x0)
  delta = zero(x0)
  delta2 = zero(x0)
  for x in A
    n += 1
    if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
      delta .= x .- mean
      mean .+= delta./n
      delta2 .= x .- mean
      M2 .+= delta.*delta2
    else
      delta = x .- mean
      mean += delta./n
      delta2 = x .- mean
      M2 += delta.*delta2
    end
  end
  if n < 2
    return NaN
  else
    if bessel
      if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
        M2 .= M2 ./ (n .- 1)
      else
        M2 = M2 ./ (n .- 1)
      end
    else
      if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
        M2 .= M2 ./ n
      else
        M2 = M2 ./ n
      end
    end
    return mean,M2
  end
end

function componentwise_meancov(A,B;bessel=true)
  x0 = first(A)
  y0 = first(B)
  n = 0
  meanx = zero(x0)
  meany = zero(y0)
  C = zero(x0)
  dx = zero(x0)
  for (x,y) in zip(A,B)
    n += 1
    if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
      dx .= x .- meanx
      meanx .+= dx./n
      meany .+= (y.-meany)./n
      C .+= dx .* (y .- meany)
    else
      dx = x .- meanx
      meanx += dx./n
      meany += (y.-meany)./n
      C += dx .* (y .- meany)
    end
  end
  if n < 2
    return NaN
  else
    if bessel
      if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
        C .= C ./ (n .- 1)
      else
        C = C ./ (n .- 1)
      end
    else
      if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
        C .= C ./ n
      else
        C = C ./ n
      end
    end
    return meanx,meany,C
  end
end

function componentwise_meancor(A,B;bessel=true)
  mx,my,cov = componentwise_meancov(A,B;bessel=bessel)
  mx,vx = componentwise_meanvar(A;bessel=bessel)
  my,vy = componentwise_meanvar(B;bessel=bessel)
  if typeof(vx) <: AbstractArray
    vx .= sqrt.(vx)
    vy .= sqrt.(vy)
  else
    vx = sqrt.(vx)
    vy = sqrt.(vy)
  end
  mx,my,cov./(vx.*vy)
end

function componentwise_weighted_meancov(A,B,W;weight_type=:reliability)
  x0 = first(A)
  y0 = first(B)
  w0 = first(W)
  n = 0
  meanx = zero(x0)
  meany = zero(y0)
  wsum = zero(w0)
  wsum2 = zero(w0)
  C = zero(x0)
  dx = zero(x0)
  for (x,y,w) in zip(A,B,W)
    n += 1
    if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
      wsum .+= w
      wsum2 .+= w.*w
      dx .= x .- meanx
      meanx .+= (w ./ wsum) .* dx
      meany .+= (w ./ wsum) .* (y .- meany)
      C .+= w .* dx .* (y .- meany)
    else
      wsum += w
      wsum2 += w.*w
      dx = x .- meanx
      meanx += (w ./ wsum) .* dx
      meany += (w ./ wsum) .* (y .- meany)
      C += w .* dx .* (y .- meany)
    end
  end
  if n < 2
    return NaN
  else
    if weight_type == :population
      if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
        C .= C ./ wsum
      else
        C = C ./ wsum
      end
    elseif weight_type == :reliability
      if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
        C .= C ./ (wsum .- wsum2 ./ wsum)
      else
        C = C ./ (wsum .- wsum2 ./ wsum)
      end
    elseif weight_type == :frequency
      if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
        C .= C ./ (wsum .- 1)
      else
        C = C ./ (wsum .- 1)
      end
    else
      error("The weight_type which was chosen is not allowed.")
    end
    return meanx,meany,C
  end
end

export get_timestep, get_timepoint, apply_timestep, apply_timepoint,
       componentwise_vectors_timestep, componentwise_vectors_timepoint

export componentwise_mean, componentwise_meanvar

export timestep_mean, timestep_median, timestep_quantile, timestep_meanvar,
       timestep_meancov, timestep_meancor, timestep_weighted_meancov

export timeseries_steps_mean, timeseries_steps_median, timeseries_steps_quantile,
       timeseries_steps_meanvar, timeseries_steps_meancov,
       timeseries_steps_meancor, timeseries_steps_weighted_meancov

export timepoint_mean, timepoint_median, timepoint_quantile,
       timepoint_meanvar, timepoint_meancov,
       timepoint_meancor, timepoint_weighted_meancov

export timeseries_point_mean, timeseries_point_median, timeseries_point_quantile,
       timeseries_point_meanvar, timeseries_point_meancov,
       timeseries_point_meancor, timeseries_point_weighted_meancov

end


text:
using DifferentialEquations
using MATLAB
using PyCall
pygui(:qt)
using PyPlot
pygui(true)
include("leap_frog_definitions.jl")
include("PSS_Definitions_Dimer_X.jl")
include("orbit_plots_leapfrog_dimer.jl")

H=.075
h=1/(2*H)
max_time=20
a=(2+h-2*sqrt(h+1))/h
Y=Yfind_Aref(1e-3,0,H);
Yfind(h)=sqrt(h/(2h+1));
P=Yfind(H)

q10=0
# q2=0
# p2=P

# q20=.3587
q20=.16
# q20=0.0
p20=0

p10=P1_poly(q20,p20,H)

Q10=(q10+q20)/sqrt(2)
Q20=(q10-q20)/sqrt(2)
P10=(p10[1]+p20)/sqrt(2)
P20=(p10[1]-p20)/sqrt(2);
println(p10)

u0=[Q10;P20;Q20;P10;0;0] # X P Q Y
P0=[Q10,Q20]; Q0=[P10,P20];
tspan = (0.0,max_time);

println(Aref_Hamiltonian(Q0,P0,0))
prob = ODEProblem(Eq_of_M_LAB_FRAME,u0,tspan)
t,A=solve(prob,Tsit5(),reltol=1e-14,abstol=1e-14,maxiters=1e15);

N=length(t.t)
Z=A[1,:]+im *A[2,:];
W=A[3,:]+im *A[4,:];

Q1=A[1,:]
Q2=A[3,:]
P1=A[4,:]
P2=A[2,:];

# figure()
# plot(Q1,P1,"b",linewidth=.1)
# plot(Q2,P2,"r",linewidth=.1)
# title("Aref coordinates");

q1=(Q1+Q2)/sqrt(2)
q2=(Q1-Q2)/sqrt(2)
p1=(P1+P2)/sqrt(2)
p2=(P1-P2)/sqrt(2)

figure()
plot(q1,p1,"b",linewidth=.1)
plot(q2,p2,"r",linewidth=.1)
title("Dimer coordinates");

Centeroid=A[5,:]+im *A[6,:];

lin_impulse=im*(1+a )*ones(N);
z_1pos=.5*(Centeroid+Z+lin_impulse+W);
z_1neg=.5*(Centeroid+Z-lin_impulse-W);
z_2pos=.5*(Centeroid-Z+lin_impulse-W);
z_2neg=.5*(Centeroid-Z-lin_impulse+W);

z_1posX=real(z_1pos)
z_1posY=imag(z_1pos)

z_2posX=real(z_2pos)
z_2posY=imag(z_2pos)

z_1negX=real(z_1neg)
z_1negY=imag(z_1neg)

z_2negX=real(z_2neg)
z_2negY=imag(z_2neg)

figure()
plot(z_1posX,z_1posY,c="b",linewidth=1)
plot(z_1negX,z_1negY,c="r",linewidth=1)
plot(z_2posX,z_2posY,c="b",linewidth=.5,linestyle="--")
plot(z_2negX,z_2negY,c="r",linewidth=.5,linestyle="--")
ylim(-3.2,3.2)
axis("equal")
title("Lab coordinates");

Q_PSS,P_PSS=PSS_function(q20,p20, H,  1e5)
mat"figure();set(gcf, 'Position',  [0, 0, 1500, 1500]); hold on;"
mat"xlabel('Q')"
mat"ylabel('P')"
Q_fix=sqrt(6)/3

mat"axis([-3, 3, -3, 3])"
mat"plot($Q_fix,0,'b.','MarkerSize',30)"
mat"plot(-$Q_fix,0,'b.','MarkerSize',30)"
mat"plot(0,$P,'r.','MarkerSize',30)"
mat"plot(0,-$P,'r.','MarkerSize',30)"
mat"plot($Q_PSS,$P_PSS,'k.','MarkerSize',5); hold on;"


text:
"""
    deploy(
        source, [dir="."], [targets...=html];
        versioned=true,
        named=false,
        force=false,
        label=nothing,
    )

Build the `source` using the given `targets` list in the `dir` directory.
`source` can be either a `Module` or a `String` representing a `Project.toml`
path.

Keyword arguments can be used to control resulting directory structure.

{#keywords}
  - `versioned` and `named`.

    These keywords will place the built files in either a versioned subdirectory,
    or a named subdirectory of `dir`, or both (with name superceding version).

    The values for `name` and `version` are taken from those provided in the
    project's `Project.toml` file. If these values are not specified then the
    "deployment" will fail.

  - `force` will remove the calculated build path prior to building
    if it already exists.

  - `label` specifies a temporal folder name to copy the finished build to.
    This can be used to build a "tracking" version of documentation such as a
    "dev" or "stable" that changes over time while still retaining the exact
    versioned builds.

## Examples

In the following examples our project will be the `Publish` package. This can
be switched out for any other project source, such as a Julia package or a
simple `Project.toml` file.

```julia
deploy(Publish, "build")
```

writes the output to the `"build"` subdirectory of the current directory.
There will be a `build/<version>` folder containing HTML content.

```julia
deploy(Publish, "build", pdf)
```

does the same as above, but build the [`pdf`](#) output instead.

```julia
deploy(Publish, "all-docs", pdf, html)
```

or build everything at once.

The keyword arguments control other aspects of the build, as shown
[above](#keywords). For example,

```julia
deploy(Publish, "ecosystem"; named=true)
```

would build `Publish` documentation to an `ecosystem/Publish/<version>`
subdirectory.
"""
function deploy(
    source,
    dir::AbstractString=".",
    targets...=html;
    versioned::Bool=true,
    named::Bool=false,
    force::Bool=false,
    label::AbstractString="",
    root::AbstractString="/",
    kws...
)
    startswith(root, '/') || error("'root' keyword must be an absolute path.")
    p = Project(source)
    name = named ? p.env["name"] : ""
    version = versioned ? p.env["version"] : ""
    parts = filter(!isempty, [dir, name, version])
    path = joinpath(parts...)
    force && rm(path; recursive=true, force=true)
    if isdir(path)
        @warn "'$path' already exists. Use force to overwrite it."
    else
        for target in targets
            target(p, path)
        end
    end
    if !isempty(label)
        ## Build the project for the given `label` as well.
        to = joinpath(filter(!isempty, [dir, name, label])...)
        rm(to; recursive=true, force=true)
        for target in targets
            target(p, to)
        end
    end
    if versioned
        ## Find the current versions.
        dir = joinpath(filter(!isempty, [dir, name])...)
        versions = sort!(filter(s->!startswith(s, '.'), readdir(dir)); rev=true)
        ## Write a version.js file containing the list of all versions.
        io = IOBuffer()
        println(io, "var PUBLISH_VERSIONS = [")
        for v in versions
            println(io, " "^4, "[", repr(v), ",", repr(abspath(joinpath(root, dir, v, "index.html"))), "],")
        end
        println(io, "];")
        ## Create/update versions.js file for every built version.
        for v in versions
            file = joinpath(dir, v, "versions.js")
            open(file, "w") do f
                println(f, "var PUBLISH_ROOT = '$(abspath(joinpath(root, dirname(file))))';")
                println(f, "var PUBLISH_VERSION = $(repr(v));")
                write(f, seekstart(io))
            end
        end
    end
    return source
end


text:
using MPI
using MPIStragglers

MPI.Init()
const comm = MPI.COMM_WORLD
const root = 0
const rank = MPI.Comm_rank(comm)
# const nworkers = MPI.Comm_size(comm) - 1
isroot() = MPI.Comm_rank(comm) == root
const data_tag = 0
const control_tag = 1
const nelems = 100 # number of 64-bit floats sent in each direction per iteration and worker

function shutdown(pool::StragglerPool)
    for i in pool.ranks
        MPI.Isend(zeros(1), i, control_tag, comm)
    end
end

function root_main()

    # parse command-line-arguments
    nworkers = MPI.Comm_size(comm) - 1
    println((nworkers, rank, root))
    print("starting with arguments $ARGS\n")
    nsamples = parse(Int, ARGS[1])
    0 < nsamples || throw(ArgumentError("The number of samples must be positive"))
    filename = ARGS[2]

    # setup
    pool = StragglerPool(nworkers)
    sendbuf = Vector{Float64}(undef, nelems)
    isendbuf = Vector{Float64}(undef, nworkers*length(sendbuf))
    recvbuf = Vector{Float64}(undef, nworkers*nelems)
    irecvbuf = copy(recvbuf)

    # warmup (wait for all workrs)
    print("warming up\n")
    epoch = 1
    kmap!(sendbuf, recvbuf, isendbuf, irecvbuf, nworkers, epoch, pool, comm; tag=data_tag)

    # output file
    f = open(filename, "a")
    write(f, "nworkers, nwait, nMB_m2w, nMB_w2m, rtt\n")
    nMB_m2w = nelems * 64/8 / 1e6 # MB sent from master to workers
    nMB_w2m = nelems * 64/8 / 1e6 # MB sent from workers to master

    # collect rtt samples
    print("starting test\n")
    for k in 1:nworkers
        for _ in 1:nsamples
            epoch += 1
            rtt = @elapsed kmap!(sendbuf, recvbuf, isendbuf, irecvbuf, k, epoch, pool, comm; tag=data_tag)
            write(f, "$nworkers, $k, $nMB_m2w, $nMB_w2m, $rtt\n")
        end
    end

    print("tests done\n")
    close(f)
    shutdown(pool)
end

function worker_main()
    recvbuf = Vector{Float64}(undef, nelems)
    sendbuf = Vector{Float64}(undef, nelems)
    crreq = MPI.Irecv!(zeros(1), root, control_tag, comm)
    while true
        rreq = MPI.Irecv!(recvbuf, root, data_tag, comm)
        index, _ = MPI.Waitany!([crreq, rreq])
        if index == 1 # exit message on control channel
            break
        end
        sendbuf .= recvbuf
        sreq = MPI.Isend(sendbuf, root, data_tag, comm)
    end
end

if isroot()
    root_main()
else
    worker_main()
end
MPI.Barrier(comm)


text:
using MacroTools, SymEngine
import MacroTools.postwalk

"""
    abs_to_relu!(ex::Expr)
    abs_to_relu(ex::Expr)
function to convert an `abs(x)` expression to `relu(x) + relu(-x)`. Mutating and nonmutating versions.
"""
function abs_to_relu(ex::Expr)
    check_for_abs(ex)
    x = ex.args[2] # the thing getting abs'ed
    return :(relu($x) + relu(-$x))
end
function abs_to_relu!(ex::Expr)
    check_for_abs(ex)
    x = ex.args[2] # the thing getting abs'ed
    ex.args = [:+, :(relu($x)), :(relu(-$x))]
end

function check_for_abs(ex::Expr)
    ex.args[1] == :abs || throw(ArgumentError("Not an absolute value expression"))
    @assert length(ex.args) == 2 "Malformed expression. `abs` can take only one argument. Got $ex"
    return nothing
end

function check_for_maxmin(ex::Expr)
    ex.args[1] ∈ (:max, :min) || throw(ArgumentError("function is neither max nor min. Got $(ex.args[1])"))
    @assert length(ex.args) == 3 "max/min_to_abs can't handle more than two inputs at the moment."
    return nothing
end

function maxmin_to_abs(ex::Expr)
    check_for_maxmin(ex)
    a, b = ex.args[2:3]
    if ex.args[1] == :max
        return :(0.5*($a + $b + abs($a-$b)))
    elseif ex.args[1] == :min
        return :(0.5*($a + $b - abs($a-$b)))
    end
end

function max0torelu(ex::Expr)
    @assert(length(ex.args)==3)
    a, b = ex.args[2:3]
    if (a == 0)
        return :(relu($b))
    elseif (b == 0)
        return :(relu($a))
    end
end

to_relu_expression(not_ex) = not_ex
function to_relu_expression(ex::Expr)
    ex.head != :call && return ex
    if ex.args[1] ∈ (:max, :min)
        if (ex.args[1] == :max) & ((ex.args[2] == 0) | (ex.args[3] == 0))
            ex = max0torelu(ex)
        else
            ex = maxmin_to_abs(ex)
        end
    elseif ex.args[1] == :abs
        ex = abs_to_relu(ex) end
    for i in 1:length(ex.args)
        ex.args[i] = to_relu_expression(ex.args[i])
    end
    return ex
end

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

"""
Construct a symbolic expression for a line between the points (x₀, 1.0) on the left and (x₁, 0) on the right.
`pos_unit` has positive slope while `neg_unit` has negative slope. Note that due to the left-to-right assumption

    neg_unit(x₀, x₁) == pos_unit(x₁, x₀)
"""
neg_unit(x0, x1) = :($(1/(x0-x1)) * (x - $x1))
"""
Construct a symbolic expression for a line between the points (x₀, 0.0) on the left and (x₁, 1.0) on the right.
`pos_unit` has positive slope while `neg_unit` has negative slope. Note that due to the left-to-right assumption

    neg_unit(x₀, x₁) == pos_unit(x₁, x₀)
"""
pos_unit(x0, x1) = :($(1/(x1-x0)) * (x - $x0))

"""
    closed_form_piecewise_linear(pts)::Expr

Constructs a closed-form piecewise linear expression from an ordered (left to right) sequence of points.
The method is inspired by the paper by Lum and Chua (cite) that considers a piecewise linear function of the form:
`f(x) = Σᵢ(yᵢ⋅gᵢ(xᵢ))` where `gᵢ(xⱼ) = δᵢⱼ`
Here in the 1D case, `gᵢ = max(0, yᵢ*min(L1, L2))`, where `L1` and `L2` are the lines "ramping up" towards xᵢ
and "ramping down" away from xᵢ. The function returns an `Expr` based on a variable `x`.
This can be turned into a callable function `f` by running something like `eval(:(f(x) = \$expression_of_x))`.

# Example
    julia> pts = [(0,0), (1,1), (2, 0)]
    3-element Array{Tuple{Int64,Int64},1}:
     (0, 0)
     (1, 1)
     (2, 0)

    julia> closed_form_piecewise_linear(pts)
    :(max(0, 0 * (-1.0 * (x - 1))) + max(0, 1 * min(1.0 * (x - 0), -1.0 * (x - 2))) + max(0, 0 * (1.0 * (x - 1))))
"""
function closed_form_piecewise_linear(pts)
    n = length(pts)
    x, y = first.(pts), last.(pts) # split the x and y coordinates
    G = []
    for i in 2:n-1
        x0, x1, x2 = x[i-1:i+1] # consider the "triangulation" of points x0,x1,x2
        L1 = pos_unit(x0, x1) # x0-x1 is an increasing linear unit
        L2 = neg_unit(x1, x2) # x1-x2 is a decreasing linear unit
        gᵢ = :($(y[i]) * max(0, min($L1, $L2)))
        push!(G, gᵢ)
    end
    # first and last points are special cases that ignore the min
    g₀ = :($(y[1]) * max(0, $(neg_unit(x[1], x[2]))))
    gᵣ = :($(y[end]) * max(0, $(pos_unit(x[end-1], x[end]))))
    # Order doesn't matter now but for our debugging purposes earlier we enforce sequential ordering.
    pushfirst!(G, g₀)
    push!(G, gᵣ)
    return :(+$(G...))
end

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

function make_expr_dict(ex)
    D = Dict()
    ex = postwalk(ex) do e
        if is_relu_expr(e)
            return get!(D, e, Symbol("z$(length(D)+1)"))
        end
        return e
    end
    D[ex] = Symbol("z$(length(D)+1)")
    return D
end

is_relu_expr(ex) = ex isa Expr && ex.head == :call && ex.args[1] == :relu

simplify(ex::Expr) = postwalk(e -> _simplify(e), ex)
_simplify(s) = s
_simplify(e::Expr) = Meta.parse(string(expand(Basic(e))))

# thinking/scripting:
Base.occursin(x, y) = x == y
function Base.occursin(needle::Union{Symbol, Expr}, haystack::Expr)
    needle == haystack && return true
    for arg in haystack.args
        occursin(needle, arg) && return true
    end
    return false
end

# Type piracy again:
Base.Expr(B::Basic) = Meta.parse(string(B))

function get_symbols(ex::Union{Expr, Symbol})
    syms = Symbol[]
    ops = (:*, :+, :-, :relu)
    postwalk(e -> e isa Symbol && e ∉ ops ? push!(syms, e) : nothing, ex)
    unique(syms)
end

function count_min_max(expr)
    """ count number of min and max in expr"""
    c = [0, 0]
    if expr == :min
        c[1] = 1
    elseif expr == :max
        c[2] = 1
    elseif expr isa Expr
        for arg in expr.args
            c += count_min_max(arg)
        end
    end
    return c
end


text:
module Utility

export PrintMatrix
include("Utility/PrintMatrix.jl")

end


text:
using FactCheck2TestSets
using Base.Test

include("test_comments.jl")
include("test_parsercombinator.jl")
include("test_jlparse.jl")
include("test_migration.jl")
include("test_sugar.jl")


text:
info("Compare JuLIP EAM Implementation against ASE EAM Implementation")

info("Loading libraries...")
using PyCall, JuLIP, ASE
@pyimport ase.calculators.eam as eam

pot_file = Pkg.dir("JuLIP") * "/data/w_eam4.fs"

info("Generate the ASE potential")
eam4_ase = eam.EAM(potential=pot_file) |> ASECalculator

info("Generate low-accuracy JuLIP potential")
eam4_jl1 = EAM(pot_file)
info("Generate medium-accuracy JuLIP potential")
eam4_jl2 = EAM(pot_file; s = 1e-4)
info("Generate high-accuracy JuLIP potential")
eam4_jl3 = EAM(pot_file; s = 1e-6)

info("Generate test configurations")
at1 = rattle!(bulk(:W, cubic=true) * 3, 0.1)
at2 = deleteat!(bulk(:W, cubic=true) * 3, 1)
at1_ase = ASEAtoms(at1)
at2_ase = ASEAtoms(at2)

info("Test 1")
println("   Low Accuracy energy error:",
      (energy(eam4_ase, at1_ase) - energy(eam4_jl1, at1)) / length(at1))
println("Medium Accuracy energy error:",
      (energy(eam4_ase, at1_ase) - energy(eam4_jl2, at1)) / length(at1))
println("  High Accuracy energy error:",
      (energy(eam4_ase, at1_ase) - energy(eam4_jl3, at1)) / length(at1))

info("Test 2")
println("   Low Accuracy energy error:",
      (energy(eam4_ase, at2_ase) - energy(eam4_jl1, at2)) / length(at2))
println("Medium Accuracy energy error:",
      (energy(eam4_ase, at2_ase) - energy(eam4_jl2, at2)) / length(at2))
println("  High Accuracy energy error:",
      (energy(eam4_ase, at2_ase) - energy(eam4_jl3, at2)) / length(at2))


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ f66e452c-97db-11eb-23fc-bd95da10e2d6
using Revise, FourierTools, Plots, FFTW, TestImages, PlutoUI, ImageShow

# ╔═╡ 8ea180a2-97dc-11eb-2276-43d5c21fc168
md"## Minimal Example"

# ╔═╡ 3b40d714-97dc-11eb-1d06-d3077fede734
begin
	arr = [1.0 2.0 3.0; 4.0 5.0 6.0]
	shift!(arr, (1, 2))
end

# ╔═╡ 48e01ee8-97dc-11eb-2330-93089d4b14b4
begin
	arr2 = [0, 1.0, 0.0, 1.0]
	shift!(arr2, 0.5)
end

# ╔═╡ 882699d8-97dc-11eb-1d71-131242e03e5f
md"# Subpixel shift"

# ╔═╡ 94433620-97dc-11eb-34df-6b93b05824cf
begin
	f(x) = cos(4π * x / 30)
	x1 = 1:30
	x2 = x1 .+ 3
end

# ╔═╡ f682339c-97dc-11eb-14c6-59570cc41c0f
@bind offset Slider(-5:0.01:5)

# ╔═╡ 3610c3aa-97dd-11eb-0053-733c17d8e7c3
md"offset=$offset"

# ╔═╡ 942ab232-97dc-11eb-18d2-d33374a89efa
begin
	y1 = f.(x1)
	y2 = f.(x2)
	y3 = shift(y2, tuple(offset))
end

# ╔═╡ d7d7fe18-97dc-11eb-1800-d1dc9ba78574
begin
	plot(y1, label="Original signal")
	plot!(y2, label="Shifted signal")
	plot!(y3, label="Fourier shifted with $offset")
end

# ╔═╡ 880650ec-97dc-11eb-37ef-451f1431c32f
md"## Image Shifting"

# ╔═╡ 5167d98e-97dc-11eb-38ef-ebe0005635d5
img = testimage("fabio_512_gray")

# ╔═╡ 5a1a8d06-97dc-11eb-3dec-452bd1fb2d65
Gray.(shift(Float64.(img), (100.12, 0.1)))

# ╔═╡ 7af79a8c-97dc-11eb-107e-2f25d880eb16
Gray.(shift(Float64.(img), (100.12, -100)))

# ╔═╡ Cell order:
# ╠═f66e452c-97db-11eb-23fc-bd95da10e2d6
# ╟─8ea180a2-97dc-11eb-2276-43d5c21fc168
# ╠═3b40d714-97dc-11eb-1d06-d3077fede734
# ╠═48e01ee8-97dc-11eb-2330-93089d4b14b4
# ╠═882699d8-97dc-11eb-1d71-131242e03e5f
# ╠═94433620-97dc-11eb-34df-6b93b05824cf
# ╠═3610c3aa-97dd-11eb-0053-733c17d8e7c3
# ╠═f682339c-97dc-11eb-14c6-59570cc41c0f
# ╠═942ab232-97dc-11eb-18d2-d33374a89efa
# ╠═d7d7fe18-97dc-11eb-1800-d1dc9ba78574
# ╟─880650ec-97dc-11eb-37ef-451f1431c32f
# ╠═5167d98e-97dc-11eb-38ef-ebe0005635d5
# ╠═5a1a8d06-97dc-11eb-3dec-452bd1fb2d65
# ╠═7af79a8c-97dc-11eb-107e-2f25d880eb16


text:
using GMRF
using SparseArrays, Random
using Test

m₁ = 15
m₂ = 15
m = m₁*m₂

V = collect(1:m)
A = sort(rand(1:m, 2))
B = setdiff(V, A)

order = 1
F = iGMRF(m₁,m₂,order, 10)
X = rand(F)
l = logpdf(F,X)
fc = fullconditionals(F,X)
l = fullcondlogpdf(F,X)
fc = getconditional(F, B, X[B])

order = 2
F = iGMRF(m₁,m₂,order, 10)
X = rand(F)
l = logpdf(F,X)
fc = fullconditionals(F,X)
l = fullcondlogpdf(F,X)
fc = getconditional(F, B, X[B])

Random.seed!(12)

@testset "GMRF.jl" begin
    include("graphstructure_test.jl")
    include("gridstructure_test.jl")
    include("igmrf_test.jl")
end


text:
function print_header(method::IPOptimizer)
    @printf "Iter     Lagrangian value Function value   Gradient norm    |==constr.|      μ\n"
end

function Base.show(io::IO, t::OptimizationState{<:Real, <:IPOptimizer})
    md = t.metadata
    @printf io "%6d   %-14e   %-14e   %-14e   %-14e   %-6.2e\n" t.iteration md["Lagrangian"] t.value t.g_norm md["ev"] md["μ"]
    if !isempty(t.metadata)
        for (key, value) in md
            key ∈ ("Lagrangian", "μ", "ev") && continue
            @printf io " * %s: %s\n" key value
        end
    end
    return
end

function Base.show(io::IO, tr::OptimizationTrace{<:Real, <:IPOptimizer})
    @printf io "Iter     Lagrangian value Function value   Gradient norm    |==constr.|      μ\n"
    @printf io "------   ---------------- --------------   --------------   --------------   --------\n"
    for state in tr
        show(io, state)
    end
    return
end

function trace!(tr, d, state, iteration, method::IPOptimizer, options, curr_time=time())
    dt = Dict()
    dt["Lagrangian"] = state.L
    dt["μ"] = state.μ
    dt["ev"] = abs(state.ev)
    dt["time"] = curr_time
    if options.extended_trace
        dt["α"] = state.alpha
        dt["x"] = copy(state.x)
        dt["g(x)"] = copy(state.g)
        dt["h(x)"] = copy(state.H)
        if !isempty(state.bstate)
            dt["gtilde(x)"] = copy(state.gtilde)
            dt["bstate"] = copy(state.bstate)
            dt["bgrad"] = copy(state.bgrad)
            dt["c"] = copy(state.constr_c)
        end
    end
    g_norm = norm(state.g, Inf) + norm(state.bgrad, Inf)
    Optim.update!(tr,
            iteration,
            value(d),
            g_norm,
            dt,
            options.store_trace,
            options.show_trace,
            options.show_every,
            options.callback)
end


text:
using IterableTables
using IteratorInterfaceExtensions
using DataFrames
using DataValues
using Test

@testset "DataFrames" begin

source_df = DataFrame(a=Union{Int,Missing}[1,2,3], b=Union{Float64,Missing}[1.,2.,3.], c=Union{String,Missing}["A","B","C"])

@test IteratorInterfaceExtensions.isiterable(source_df) == true

as_array = collect(IteratorInterfaceExtensions.getiterator(source_df))

@test length(as_array) == 3
@test as_array == [(a=DataValue(1),b=DataValue(1.),c=DataValue("A")),
    (a=DataValue(2),b=DataValue(2.),c=DataValue("B")),
    (a=DataValue(3),b=DataValue(3.),c=DataValue("C"))]

data = []
push!(data, [1,2,3])
push!(data, [1.,2.,3.])
push!(data, ["A","B","C"])
source_df_non_nullable = DataFrame(data, [:a,:b,:c])

non_nullable_as_array = collect(IteratorInterfaceExtensions.getiterator(source_df_non_nullable))

@test length(non_nullable_as_array) == 3
@test as_array == [(a=1,b=1.,c="A"),(a=2,b=2.,c="B"),(a=3,b=3.,c="C")]

source_array_non_nullable = [(a=1,b=1.,c="A"), (a=2,b=2.,c="B"), (a=3,b=3.,c="C")]
df = DataFrame(source_array_non_nullable)

@test size(df) == (3,3)
@test isa(df[:a], Array)
@test isa(df[:b], Array)
@test isa(df[:c], Array)
@test df[:a] == [1,2,3]
@test df[:b] == [1.,2.,3.]
@test df[:c] == ["A","B","C"]

source_array = [(a=DataValue(1),b=DataValue(1.),c=DataValue("A")), (a=DataValue(2),b=DataValue(2.),c=DataValue("B")), (a=DataValue(3),b=DataValue(3.),c=DataValue("C"))]
df = DataFrame(source_array)

@test size(df) == (3,3)
@test isa(df[:a], Vector{Union{Int,Missing}})
@test isa(df[:b], Vector{Union{Float64,Missing}})
@test isa(df[:c], Vector{Union{String,Missing}})
@test df[:a] == [1,2,3]
@test df[:b] == [1.,2.,3.]
@test df[:c] == ["A","B","C"]

end


text:
export ElementWiseFunctorType, ElementWiseFunctors
export get_num_args
@compat abstract type ElementWiseFunctorType{NArg} end

get_num_args{NArg}(::ElementWiseFunctorType{NArg}) = NArg

module ElementWiseFunctors
using ..Mocha
type Add <: ElementWiseFunctorType{2}
end
type Subtract <: ElementWiseFunctorType{2}
end
type Multiply <: ElementWiseFunctorType{2}
end
type Divide <: ElementWiseFunctorType{2}
end
end # module ElementWiseFunctors

############################################################
# Element-wise operation layer
############################################################
@defstruct ElementWiseLayer Layer (
  name :: AbstractString = "element-wise",
  operation :: ElementWiseFunctorType = ElementWiseFunctors.Add(),
  (tops :: Vector{Symbol} = Symbol[], length(tops) == 1),
  (bottoms :: Vector{Symbol} = Symbol[], length(bottoms) == get_num_args(operation)),
)
@characterize_layer(ElementWiseLayer,
  can_do_bp => true
)

type ElementWiseLayerState{Op<:ElementWiseFunctorType} <: LayerState
  layer      :: ElementWiseLayer
  blobs      :: Vector{Blob}
  blobs_diff :: Vector{Blob}
end

function setup(backend::Backend, layer::ElementWiseLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  @assert all(map(i -> size(inputs[i]) == size(inputs[1]), 2:length(inputs)))
  blobs = Blob[make_blob(backend, eltype(inputs[1]), size(inputs[1]))]
  if all(map(b -> isa(b, NullBlob), diffs))
    blobs_diff = Blob[NullBlob()]
  else
    blobs_diff = Blob[make_blob(backend, eltype(inputs[1]), size(inputs[1]))]
  end

  return ElementWiseLayerState{typeof(layer.operation)}(layer, blobs, blobs_diff)
end
function shutdown(backend::Backend, state::ElementWiseLayerState)
  map(destroy, state.blobs)
  map(destroy, state.blobs_diff)
end

for (functor, op) in ((ElementWiseFunctors.Add, (+)),
                      (ElementWiseFunctors.Subtract, (-)),
                      (ElementWiseFunctors.Multiply, (*)),
                      (ElementWiseFunctors.Divide, (/)))
  @eval begin
    # I'm getting the following warning unless I extract the for loop
    # as a separate function with clear type annotations.
    # Warning: could not attach metadata for @simd loop.
    function functor_impl{T}(::$functor, input1::Array{T}, input2::Array{T}, output::Array{T})
      len = length(input1)
      @simd for i = 1:len
        @inbounds output[i] = $op(input1[i], input2[i])
      end
    end
    function forward(backend::CPUBackend, state::ElementWiseLayerState{$functor},
        inputs::Vector{Blob})

      input1 = inputs[1].data
      input2 = inputs[2].data
      output = state.blobs[1].data
      functor_impl(state.layer.operation, input1, input2, output)
    end
  end
end

function backward(backend::Backend, state::ElementWiseLayerState{ElementWiseFunctors.Add},
    inputs::Vector{Blob}, diffs::Vector{Blob})

  for i = 1:length(diffs)
    if !isa(diffs[i], NullBlob)
      copy!(diffs[i], state.blobs_diff[1])
    end
  end
end
function backward(backend::CPUBackend, state::ElementWiseLayerState{ElementWiseFunctors.Subtract},
    inputs::Vector{Blob}, diffs::Vector{Blob})

  if !isa(diffs[1], NullBlob)
    copy!(diffs[1], state.blobs_diff[1])
  end
  if !isa(diffs[2], NullBlob)
    copy!(diffs[2], state.blobs_diff[1])
    BLAS.scal!(length(diffs[2]), convert(eltype(diffs[2]),-1), diffs[2].data, 1)
  end
end
function backward(backend::CPUBackend, state::ElementWiseLayerState{ElementWiseFunctors.Multiply},
    inputs::Vector{Blob}, diffs::Vector{Blob})

  for i = 1:length(diffs)
    if !isa(diffs[i], NullBlob)
      copy!(diffs[i], state.blobs_diff[1])
      Vec.mul!(diffs[i].data, inputs[i%2 + 1].data)
    end
  end
end
function backward(backend::CPUBackend, state::ElementWiseLayerState{ElementWiseFunctors.Divide},
    inputs::Vector{Blob}, diffs::Vector{Blob})

  if !isa(diffs[1], NullBlob)
    copy!(diffs[1], state.blobs_diff[1])
    Vec.div!(diffs[1].data, inputs[2].data)
  end
  if !isa(diffs[2], NullBlob)
    copy!(diffs[2], state.blobs_diff[1])
    Vec.mul!(diffs[2].data, state.blobs[1].data)
    Vec.div!(diffs[2].data, inputs[2].data)
    BLAS.scal!(length(diffs[2]), convert(eltype(diffs[2]),-1), diffs[2].data, 1)
  end
end


text:
module mlmetrics
using  Statistics

"""
    struct mse{T<:AbstractFloat}
        y::Vector
        ŷ::Vector
    end
    The mean squared error is the sum of the squared differences between the
    predicted and observed values.

    struct mae{T<:AbstractFloat}
        y::Vector
        ŷ::Vector
    end
    The mean absolute error is the sum of the absolute differences between the
    predicted and observed values.

    struct rmse{T<:AbstractFloat}
        y::Vector
        ŷ::Vector
    end
    The root mean squared error is the square root of the mean squared error.

    struct rmae{T<:AbstractFloat}
        y::Vector
        ŷ::Vector
    end
    The root mean absolute error is the square root of the mean absolute error.

    struct rmsle{T<:AbstractFloat}
        y::Vector
        ŷ::Vector
    end
    The root mean squared log error is the square root of the mean squared log error.

    struct r2score{T<:AbstractFloat}
        y::Vector
        ŷ::Vector
    end
    The R2 is a statistic that will give some information about the goodness of fit of a model.
"""

struct mse{T<:AbstractFloat}
    y::Vector
    ŷ::Vector

    function mse(y, ŷ)
        return mean((y - ŷ).^2)
    end
end

struct mae{T<:AbstractFloat}
    y::Vector
    ŷ::Vector

    function mae(y, ŷ)
        return mean(abs.(y - ŷ))
    end
end

struct rmse{T<:AbstractFloat}
    y::Vector
    ŷ::Vector

    function rmse(y, ŷ)
        return sqrt(mse(y, ŷ))
    end
end

struct rmae{T<:AbstractFloat}
    y::Vector
    ŷ::Vector

    function rmae(y, ŷ)
        return sqrt(mae(y, ŷ))
    end
end

struct rmsle{T<:AbstractFloat}
    y::Vector
    ŷ::Vector

    function rmsle(y, ŷ)
        return sqrt(mean((log.(y .+ 1) - log.(ŷ .+ 1)).^2))
    end
end

struct r2score{T<:AbstractFloat}
    y::Vector
    ŷ::Vector

    function r2score(y, ŷ)
        squares = sum(y.^2)
        residuals = sum(y - ŷ).^2
        return 1 - residuals/squares
    end
end

struct binarycrossentropy{T<:AbstractFloat}
    y::Vector
    ŷ::Vector

    function binarycrossentropy(y::Vector, ŷ::Vector)
        ε = 1e-10
        return -mean(y.*log2.(ŷ .+ ε) + (1 .- y).*log2.(1 .- ŷ .+ ε))
    end
end

struct accuracy{T<:AbstractFloat}
    y_true::AbstractArray
    y_predict::AbstractArray

    function accuracy(y_true::AbstractArray, y_predict::AbstractArray)
        return sum(y_true.==y_predict)/length(y_true)
    end

end

struct precision{T<:AbstractFloat}
    y_true::AbstractArray
    y_predict::AbstractArray

    function precision(y_true::AbstractArray, y_predict::AbstractArray)
        true_positives = sum(y_true.== 1 .== y_predict)
        true_negatives = sum(y_true.== 0 .== y_predict)
        false_positives = sum(y_true.!= 1 .== y_predict)
        false_negatives = sum(y_true.!= 0 .== y_predict)

        return true_positives/(true_positives + false_positives)
    end

end

struct recall{T<:AbstractFloat}
    y_true::AbstractArray
    y_predict::AbstractArray

    function recall(y_true::AbstractArray, y_predict::AbstractArray)
        true_positives = sum(y_true.== 1 .== y_predict)
        true_negatives = sum(y_true.== 0 .== y_predict)
        false_positives = sum(y_true.!= 1 .== y_predict)
        false_negatives = sum(y_true.!= 0 .== y_predict)

        return true_positives/(true_positives + false_negatives)
    end

end

struct f1_score{T<:AbstractFloat}
    y_true::AbstractArray
    y_predict::AbstractArray

    function f1_score(y_true::AbstractArray, y_predict::AbstractArray)
        num =  2 * precision(y_true, y_predict) * recall(y_true, y_predict)
        den = precision(y_true, y_predict) + recall(y_true, y_predict)
        return num/den
    end

end

end