text:
@testset "BatchJob" begin
    job = BatchJob("00000000-0000-0000-0000-000000000000")

    @testset "status_reason" begin
        @testset "not provided" begin
            patch = @patch function AWSBatch.Batch.describe_jobs(args...; kwargs...)
                Dict(
                    "jobs" => [
                        Dict()
                    ]
                )
            end

            apply(patch) do
                @test status_reason(job) === nothing
            end
        end

        @testset "provided" begin
            reason = "Essential container in task exited"
            patch = @patch function AWSBatch.Batch.describe_jobs(args...; kwargs...)
                Dict(
                    "jobs" => [
                        Dict("statusReason" => reason)
                    ]
                )
            end

            apply(patch) do
                @test status_reason(job) == reason
            end
        end
    end

    @testset "log_events" begin
        @testset "Stream not yet created" begin
            # When a AWS Batch job is first submitted the description of the job will not
            # contain a reference to a log stream
            patches = log_events_patches(log_stream_name=nothing)
            apply(patches) do
                @test log_events(job) === nothing
            end
        end
    end

    @testset "wait" begin
        # Generate a patch which returns the next status each time it is requested
        function status_patch(states)
            index = 1

            return @patch function AWSBatch.Batch.describe_jobs(args...; kwargs...)
                json = Dict(
                    "jobs" => [
                        Dict("status" => states[index])
                    ]
                )

                if index < length(states)
                    index += 1
                end

                return json
            end
        end

        @testset "success" begin
            # Encounter all states possible for a successful job
            states = ["SUBMITTED", "PENDING", "RUNNABLE", "STARTING", "RUNNING", "SUCCEEDED"]
            apply(status_patch(states)) do
                @test_log logger "info" r"^[\d-]+ status \w+" begin
                    @test wait(state -> state < AWSBatch.SUCCEEDED, job; delay=0.1) == true
                end
                @test status(job) == AWSBatch.SUCCEEDED
            end
        end

        @testset "failed" begin
            # Encounter all states possible for a failed job
            states = ["SUBMITTED", "PENDING", "RUNNABLE", "STARTING", "RUNNING", "FAILED"]
            apply(status_patch(states)) do
                @test_log logger "info" r"^[\d-]+ status \w+" begin
                    @test wait(state -> state < AWSBatch.SUCCEEDED, job; delay=0.1) == true
                end
                @test status(job) == AWSBatch.FAILED
            end
        end

        @testset "timeout" begin
            apply(status_patch(["SUBMITTED", "RUNNING", "SUCCEEDED"])) do
                started = time()
                @test_nolog logger "info" r".*" begin
                    @test_throws BatchJobError wait(
                        state -> state < AWSBatch.SUCCEEDED,
                        job;
                        delay=0.1,
                        timeout=0
                    )
                end
                duration = time() - started
                @test status(job) != AWSBatch.SUCCEEDED  # Requires a minimum of 3 states
                @test duration < 1  # Less than 1 second
            end
        end
    end
end


text:
function swap(X, swap; dim = false, row_only = false, inv_perm = false, stop = false)
  @assert dim !== false
  perm = collect(1:length(dim))
  perm[swap[1]] = swap[2]
  perm[swap[2]] = swap[1]
  return permutesystems(X, perm, dim, row_only, inv_perm, stop)
end

ext: jl


text:
"""
    @pyconst ex

Equivalent to `Py(ex)` but always returns the exact same Julia object.

That is, if `foo() = @pyconst ex` then `foo() === foo()`.

The expression `ex` is evaluated the first time the code is run.

If `ex` is a string literal, the string is interned.

Do not use this macro at the top level of a module. Instead, use `pynew()` and `pycopy!()`.
"""
macro pyconst(ex)
    x = pynew()
    val = esc(ex)
    if ex isa String
        val = :($pystr_intern!($pystr($val)))
    else
        val = :($Py($val))
    end
    :(pyisnull($x) ? pycopy!($x, $val) : $x)
end
export @pyconst


text:
"""
    IndependentMOKernel(k::Kernel) <: Kernel

A Multi-Output kernel which assumes each output is independent of the other.
"""
struct IndependentMOKernel{Tkernel<:Kernel} <: MOKernel
    kernel::Tkernel
end

function (κ::IndependentMOKernel)((x, px)::Tuple{Any, Int}, (y, py)::Tuple{Any, Int})
    if px == py
        return κ.kernel(x, y)
    else
        return 0.0
    end
end

function kernelmatrix(k::IndependentMOKernel, x::MOInput, y::MOInput)
    @assert x.out_dim == y.out_dim
    temp = k.kernel.(x.x, permutedims(y.x))
    return cat((temp for _ in 1:y.out_dim)...; dims=(1,2))
end

function Base.show(io::IO, k::IndependentMOKernel)
    print(io, string("Independent Multi-Output Kernel\n\t", string(k.kernel)))
end


text:
module JLTest

_TESTCTX = :(testContext)

#Used to store test case state as testcase is executed
type TestContext
  name             #test case name
  curTest          #current test
  numRun           #num tests run
  numPassed        #num tests passed
  numFailed        #num assertions failed
  numErrors        #num assertions with errors
  numSkipped       #num tests skipped
  setUpCase        #before test case
  tearDownCase     #after test case
  setUp            #before test
  tearDown         #after test
  preAssert        #before assertion
  postAssert       #after assertion
  noprint          #disable outpt
  userData         #for user customization
  TestContext(desc) = new (desc,"",0,0,0,0,0,()->nothing, ()->nothing, ()->nothing, ()->nothing,(args...)->nothing, (args...)->nothing, false, nothing)
end

function printMsg(tc, args...;nl='\n')
  if tc.noprint;
  else
    print(args...,nl)
  end
end

function printTestReport(tc::JLTest.TestContext)
  local rep = "Run: $(tc.numRun) | Passed: $(tc.numPassed) | Failed: $(tc.numFailed) | Errors: $(tc.numErrors) | Skipped: $(tc.numSkipped)"
  local topLeft = repeat("=",int((length(rep) - length(tc.name) - 2)/2))
  local topRight = repeat("=", length(rep) - length(topLeft) - length(tc.name) - 2)
  local borderBot = repeat("=",length(rep))
  printMsg(tc,topLeft," ",tc.name," ", topRight)
  printMsg(tc,rep)
  printMsg(tc,borderBot)
end

function setUp(tc::JLTest.TestContext)
  try
    tc.numRun += 1
    tc.setUp()
  catch ex
    printMsg(tc,"Exception in setUp for ", tc.desc, " : ", ex)
  end
end

function doTest(tc::JLTest.TestContext, assertion::String, test::Function, args...)
  if !(test(args...))
    tc.numFailed += 1
    printMsg(tc,assertion,"(";nl="")
    for arg in args[1:end-1]
      printMsg(tc,arg,", ";nl="")
    end
    printMsg(tc,args[end], ") failed in ",tc.curTest)
    return false
  end
  return true
end

function handleException(tc::JLTest.TestContext, assertion, ex, bt)
  tc.numErrors += 1
  printMsg(tc,"Exception:", ex, " during ", assertion)
  if !tc.noprint
    Base.show_backtrace(STDOUT, bt)
  end
end

function tearDown(tc::JLTest.TestContext)
  try
    tc.tearDown()
  catch ex
    printMsg(tc,"Exception in tearDown for ", tc.desc, " : ", ex)
  end
end

export @quiet
macro quiet()
  quote
    $(esc(_TESTCTX)).noprint = true
  end
end

export @unquiet
macro unquiet()
  quote
    $(esc(_TESTCTX)).noprint = false
  end
end

macro assertion1(assertion,arg1,test,args...)
  local msg = ""
  if length(args) > 0
      msg = args[1]
  end
  quote
    local tc = $(esc(_TESTCTX))
    local a = $(esc(arg1))
    local success = false
    tc.preAssert(tc, $assertion,a)
    try
      success = doTest(tc, $assertion, $(esc(test)),a)
      if !success && $msg != ""
        printMsg(tc,$msg)
      end
    catch ex
      bt=catch_backtrace()
      handleException(tc, $assertion, ex,bt)
    end
    if success
      tc.postAssert(tc, $assertion,a)
    end
  end
end

macro assertion2(assertion,arg1,arg2,test,args...)
  local msg = ""
  if length(args) > 0
      msg = args[1]
  end
  quote
    local tc = $(esc(_TESTCTX))
    local a = $(esc(arg1))
    local b = $(esc(arg2))
    local success = false
    tc.preAssert(tc, $assertion, a, b)
    try
      success = doTest(tc, $assertion, $(esc(test)),a,b)
      if !success && $msg != ""
        printMsg(tc,$msg)
      end
    catch ex
      bt=catch_backtrace()
      handleException(tc, $assertion, ex,bt)
    end
    if success
      tc.postAssert(tc, $assertion, a, b)
    end
  end
end

export @casename
macro casename(str)
  quote
    $(esc(_TESTCTX)).name = $str
  end
end

export @testname
macro testname(str)
  quote
    $(esc(_TESTCTX)).curTest = $str
  end
end

export @setUp
#set a function to be called before each test starts
macro setUp(func)
  quote
    $(esc(_TESTCTX)).setUp =  $(esc(func))
    nothing
  end
end

export @tearDown
#set a function to be called after each test finishes
macro tearDown(func)
  quote
    $(esc(_TESTCTX)).tearDown =  $(esc(func))
    nothing
  end
end

export @setUpCase
#set a function to be called before testcase starts
macro setUpCase(func)
  quote
    $(esc(_TESTCTX)).setUpCase =  $(esc(func))
    nothing
  end
end

export @tearDownCase
#set a function to be called after testcase finishes
macro tearDownCase(func)
  quote
    $(esc(_TESTCTX)).tearDownCase =  $(esc(func))
    nothing
  end
end

export @setPreAssert
#set a callback function to be called before each assert
#function should have signiture f(testCtx, assertion, args...)
#where testCtx is the testcase's TestContext, assertion is a string with name of assertion
#(e.g., "assertEqual") and args are values that will be passed to the assertion
#One possible use is in conjunction with userData to capture fine grained stats
macro setPreAssert(func)
  quote
    $(esc(_TESTCTX)).preAssert =  $(esc(func))
    nothing
  end
end

export @setPostAssert
#set a callback function to be called after testcase finishes.
#See setPreAssert for details
macro setPostAssert(func)
  quote
    $(esc(_TESTCTX)).postAssert =  $(esc(func))
    nothing
  end
end

export @setUserData
#store some data in the test context. Possibly for collecting custom test stats.
macro setUserData(data)
  quote
    $(esc(_TESTCTX)).userData =  $(esc(data))
    nothing
  end
end

export @getUserData
#get user data
macro getUserData()
  quote
    $(esc(_TESTCTX)).userData
  end
end

export @assertEqual
macro assertEqual(val1,val2,args...)
  test = :((a,b)->(a == b))
  :( @assertion2("assertEqual", $(esc(val1)), $(esc(val2)), $(esc(test)), $(args...)) )
end

export @assertAproxEqual
macro assertAproxEqual(val1,val2,args...)
  local prec = 1E-4
  if length(args) > 0 && isa(args[1], Number)
    prec = args[1]
    args = args[2:end]
  end
  test = :((a,b)->(maximum(abs(a - b)) < $prec))
  :( @assertion2("assertAproxEqual", $(esc(val1)), $(esc(val2)), $(esc(test)), $(args...)) )
end

export @assertNotEqual
macro assertNotEqual(val1,val2,args...)
  test = :((a,b)->(a != b))
  :( @assertion2("assertNotEqual", $(esc(val1)), $(esc(val2)), $(esc(test)), $(args...)) )
end

export @assertLess
#val1 < val2
macro assertLess(val1,val2,args...)
  test = :((a,b)->(a < b))
  :( @assertion2("assertEqual", $(esc(val1)), $(esc(val2)), $(esc(test)), $(args...)) )
end

export @assertLessEqual
#val1 <= val2
macro assertLessEqual(val1,val2,args...)
  test = :((a,b)->(a <= b))
  :( @assertion2("assertLessEqual", $(esc(val1)), $(esc(val2)), $(esc(test)), $(args...)) )
end

export @assertGreater
#val1 > val2
macro assertGreater(val1,val2,args...)
  test = :((a,b)->(a > b))
  :( @assertion2("assertGreater", $(esc(val1)), $(esc(val2)), $(esc(test)), $(args...)) )
end

export @assertGreaterEqual
#val1 >= val2
macro assertGreaterEqual(val1,val2,args...)
  test = :((a,b)->(a >= b))
  :( @assertion2("assertGreater", $(esc(val1)), $(esc(val2)), $(esc(test)), $(args...)) )
end

export @assertIs
#an obj1 is the same as obj2
macro assertIs(obj1,obj2,args...)
  test = :((a,b)->(a === b))
  :( @assertion2("assertIs", $(esc(obj1)), $(esc(obj2)), $(esc(test)), $(args...)) )
end

export @assertIsNot
#an obj1 is not the same as obj2
macro assertIsNot(obj1,obj2,args...)
  test = :((a,b)->(a !== b))
  :( @assertion2("assertIsNot", $(esc(obj1)), $(esc(obj2)), $(esc(test)), $(args...)) )
end

export @assertIn
#an obj not in collection
macro assertIn(obj,collection,args...)
  test = :((a,b)->(a in b))
  :( @assertion2("assertIn", $(esc(obj)), $(esc(collection)), $(esc(test)), $(args...)) )
end

export @assertNotIn
#an obj not in a collection
macro assertNotIn(obj,collection,args...)
  test = :((a,b)->!(a in b))
  :( @assertion2("assertNotIn", $(esc(obj)), $(esc(collection)), $(esc(test)), $(args...)) )
end

export @assertItemsEqual
#two collections equal ignoring order difference
macro assertItemsEqual(col1,col2,args...)
  test = :((a,b)->(sort(a) == sort(b)))
  :( @assertion2("assertItemsEqual", $(esc(col1)), $(esc(col2)), $(esc(test)), $(args...)) )
end

export @assertIsA
#an obj is a type
macro assertIsA(obj,typ,args...)
  test = :((a,b)->isa(a,b))
  :( @assertion2("assertIsA", $(esc(obj)), $(esc(typ)), $(esc(test)), $(args...)) )
end

export @assertIsNotA
#an obj is not a type
macro assertIsNotA(obj,typ,args...)
  test = :((a,b)->!isa(a,b))
  :( @assertion2("assertIsNotA", $(esc(obj)),$(esc(typ)), $(esc(test)), $(args...)) )
end

export @assertMatches
#a regex matches a string
macro assertMatches(regex,str,args...)
  test = :((a,b)->ismatch(a,b))
  :( @assertion2("assertMatches", $(esc(regex)), $(esc(str)), $(esc(test)), $(args...)) )
end

export @assertNotMatches
macro assertNotMatches(regex,str,args...)
  test = :((a,b)->!ismatch(a,b))
  :( @assertion2("assertNotMatches", $(esc(regex)), $(esc(str)), $(esc(test)), $(args...)) )
end

export @assertSameType
#type of val1 is same as type of val2
macro assertSameType(val1,val2,args...)
  test = :((a,b)->(typeof(a) == typeof(b)))
  :( @assertion2("assertSameType", $(esc(val1)), $(esc(val2)), $(esc(test)), $(args...)) )
end

export @assertTypeOf
#type of val1 is typ
macro assertTypeOf(val,typ,args...)
  test = :((a,b)->(typeof(a) == b))
  :( @assertion2("assertTypeOf", $(esc(val)), $(esc(typ)), $(esc(test)), $(args...)) )
end

export @assertStrictlyEqual
#val1 == val2 and type of val1 is same as type of val2
macro assertStrictlyEqual(val1,val2,args...)
  test = :((a,b)->(a == b && typeof(a) == typeof(b)))
  :( @assertion2("assertStrictlyEqual", $(esc(val1)), $(esc(val2)), $(esc(test)), $(args...)) )
end

export @assertTrue
macro assertTrue(expr,args...)
  test = :((a)->(isa(a,Bool) && a))
  :( @assertion1("assertTrue", $(esc(expr)), $(esc(test)), $(args...)) )
end

export @assertFalse
macro assertFalse(expr,args...)
  test = :((a)->(isa(a,Bool) && !a))
  :( @assertion1("assertFalse", $(esc(expr)), $(esc(test)), $(args...)) )
end

export @testFailed(msg)
#Fail and print msg
macro testFailed(msg)
  test = :((a)->false)
  :( @assertion1("Failed", $(esc(msg)), $(esc(test))) )
end

export @expectFailures
#Use to assert that n failures are expected at the point where macro appears
#and deduct from failure count
macro expectFailures(n)
  quote
    if $(esc(_TESTCTX)).numFailed != $(esc(n))
      printMsg($(esc(_TESTCTX)), "Error: ", $(esc(n))," expected failures but ", $(esc(_TESTCTX)).numFailed, " actual")
    else
      $(esc(_TESTCTX)).numFailed -= $(esc(n))
    end
  end
end

export @expectErrors
#Use to assert that n errors are expected at the point where macro appears
#and deduct from error count
macro expectErrors(n)
  quote
    if $(esc(_TESTCTX)).numFailed != $(esc(n))
      printMsg($(esc(_TESTCTX)), "Error: ", $(esc(n))," expected errors but ", $(esc(_TESTCTX)).numErrors, " actual")
    else
      $(esc(_TESTCTX)).numErrors -= $(esc(n))
    end
  end
end

STARTFAIL = :startFail
STARTERR = :startErr

export @suppressAll
#usage @suppressAll blk - suppress all failures, errors and msgs in blk
#usage @suppressAll nFailure blk - expect nFailures and suppress all errors and msgs in blk
#usage @suppressAll nFailure nErrors blk - expect nFailures and nErrors and suppress all msgs in blk
macro suppressAll(args...)
  l = length(args)
  if  l == 1
    nFailures = nErrors = -1
    block = args[1]
  elseif l == 2
    (nFailures,block) = args
    nErrors = -1
  elseif l == 3
    (nFailures,nErrors, block)  = args
  else
    error("wrong number arguements to suppressAll: $l")
  end
  quote
    local tc = $(esc(_TESTCTX))
    $(esc(STARTFAIL)) = tc.numFailed
    $(esc(STARTERR)) = tc.numErrors
    @quiet()
    $(esc(block))
    @unquiet()
    @expectFailures($nFailures >= 0 ? $nFailures : (tc.numFailed - $(esc(STARTFAIL))))
    @expectErrors($nErrors >= 0 ? $nErrors : (tc.numErrors - $(esc(STARTFAIL))))
  end

end

export @assertThrows
macro assertThrows(args...)
  if length(args) == 1
    exceps = []
    f = args[1]
  else
    exceps = args[1:end-1]
    f = args[end]
  end
  quote
    sexceps = $exceps
    local tc = $(esc(_TESTCTX))
    try
      $(esc(f))
      tc.numFailed += 1
      printMsg(tc,"assertThrows: Expected exception", length($exceps) > 0 ? "in $sexceps " : " ","Test: ", tc.curTest)
    catch ex
      tex = typeof(ex).name.name
      if length($exceps) > 0 && !( tex in $exceps)
        tc.numFailed += 1
        printMsg(tc, "assertThrows: Expected exceptions in $sexceps got $tex. Test: ", tc.curTest)
        if !tc.noprint
          Base.show_backtrace(STDOUT, catch_backtrace())
        end
      end
    end
  end
end

export @testreport
macro testreport()
  quote
    printTestReport($(esc(_TESTCTX)))
  end
end

export @testcase
#Create test case around a block
macro testcase(args...)
  if length(args) == 1
    enabled = true
    block = args[1]
  else
    enabled = args[1]
    block = args[2]
  end
  if enabled
    quote
      let $(esc(_TESTCTX)) = JLTest.TestContext("No Name Test Case")
        local tc = $(esc(_TESTCTX))
        tc.setUpCase()
        $(esc(block))
        tc.tearDownCase()
        if tc.numErrors == 0 && tc.numFailed == 0
          printMsg(tc, tc.name, " Ok")
          if tc.numSkipped > 0
            printMsg(tc, tc.numSkipped, " Tests Skipped")
          end
        else
          printMsg(tc, tc.name, " FAILED! ", tc.numFailed, " failures ", tc.numErrors, " errors")
        end
      end
    end
  end
end

export @test
macro test(block)
  quote
    let tc = $(esc(_TESTCTX))
      numRun = tc.numRun += 1
      tc.curTest = "Test $numRun"
      local before = tc.numFailed + tc.numErrors
      tc.setUp()
      try
        $(esc(block))
      catch ex
        bt=catch_backtrace()
        handleException(tc,tc.curTest,ex,bt)
      end
      tc.tearDown()
      local after = tc.numFailed + tc.numErrors
      if before == after
        tc.numPassed +=1
      else
        printMsg(tc, tc.curTest, " Failed!")
      end
      tc.curTest = ""
    end
  end
end

export @testskip
#Place in front of a test to unconditionally skip it
macro testskip(args...)
  local skip, blk
  if length(args) == 1
    skip = true
    blk = args[1]
  elseif length(args) == 2
    (skip, blk) = args
  else
    error("@testskip accepts at most two arguments")
  end

  if skip
    quote
      $(esc(_TESTCTX)).numSkipped+=1
    end
  else
      :(@test($(blk)))
  end
end

end #JLTest

max_stars_repo_name: smangano/JLTest

max_issues_repo_name: smangano/JLTest

max_forks_repo_name: smangano/JLTest


text:
module TestQAQC

  using GBIF
  using Test

  set = occurrences()

  # Use single queries
  l1 = length(set)
  filter!(have_no_issues, set)
  @test l1 >= length(set)

  # Reset queries
  allrecords!(set)
  @test length(set) == l1

  # Using next! and complete! after filter! should work correctly
  filter!(have_no_issues, set)
  next!(set)
  @test length(set) == 2*l1

  # Using next! and complete! after filter! should work correctly
  filter!(have_no_issues, set)
  allrecords!(set)
  next!(set)
  @test length(set) == 3*l1

  # Multiple filters at once
  allrecords!(set)
  filter!.([have_ok_coordinates, have_neither_zero_coordinates, have_a_date], set)
  @test length(set) != 3*l1

  # Deprecation for qualitycontrol!
  # @test_warn "deprecated" qualitycontrol!(set; verbose=false)

end


text:
module HydrostaticBoussinesq

export HydrostaticBoussinesqModel, HydrostaticBoussinesqProblem, OceanDGModel

using StaticArrays
using LinearAlgebra: I, dot, Diagonal
using ..VariableTemplates
using ..MPIStateArrays
using ..DGmethods: init_ode_state
using ..PlanetParameters: grav
using ..Mesh.Filters: CutoffFilter, apply!, ExponentialFilter
using ..Mesh.Grids: polynomialorder, VerticalDirection

using ..DGmethods.NumericalFluxes: Rusanov, CentralGradPenalty,
                                   CentralNumericalFluxDiffusive,
                                   CentralNumericalFluxNonDiffusive

import ..DGmethods.NumericalFluxes: update_penalty!, numerical_flux_diffusive!,
                                    NumericalFluxNonDiffusive

import ..DGmethods: BalanceLaw, vars_aux, vars_state, vars_gradient,
                    vars_diffusive, vars_integrals, flux_nondiffusive!,
                    flux_diffusive!, source!, wavespeed,
                    boundary_state!, update_aux!, update_aux_diffusive!,
                    gradvariables!, init_aux!, init_state!,
                    LocalGeometry, indefinite_stack_integral!,
                    reverse_indefinite_stack_integral!, integrate_aux!,
                    DGModel, nodal_update_aux!, diffusive!,
                    copy_stack_field_down!, create_state

×(a::SVector, b::SVector) = StaticArrays.cross(a, b)
∘(a::SVector, b::SVector) = StaticArrays.dot(a, b)

abstract type OceanBoundaryCondition end
struct CoastlineFreeSlip             <: OceanBoundaryCondition end
struct CoastlineNoSlip               <: OceanBoundaryCondition end
struct OceanFloorFreeSlip            <: OceanBoundaryCondition end
struct OceanFloorNoSlip              <: OceanBoundaryCondition end
struct OceanSurfaceNoStressNoForcing <: OceanBoundaryCondition end
struct OceanSurfaceStressNoForcing   <: OceanBoundaryCondition end
struct OceanSurfaceNoStressForcing   <: OceanBoundaryCondition end
struct OceanSurfaceStressForcing     <: OceanBoundaryCondition end

abstract type HydrostaticBoussinesqProblem end

struct HydrostaticBoussinesqModel{P,T} <: BalanceLaw
  problem::P
  c₁::T
  c₂::T
  c₃::T
  αᵀ::T
  νʰ::T
  νᶻ::T
  κʰ::T
  κᶻ::T
end

HBModel   = HydrostaticBoussinesqModel
HBProblem = HydrostaticBoussinesqProblem

function OceanDGModel(bl::HBModel, grid, numfluxnondiff, numfluxdiff,
                      gradnumflux; kwargs...)
  vert_filter = CutoffFilter(grid, polynomialorder(grid)-1)
  exp_filter  = ExponentialFilter(grid, 1, 8)

  modeldata = (vert_filter = vert_filter, exp_filter=exp_filter)

  return DGModel(bl, grid, numfluxnondiff, numfluxdiff, gradnumflux;
                 kwargs..., modeldata=modeldata)
end

# If this order is changed check the filter usage!
function vars_state(m::HBModel, T)
  @vars begin
    u::SVector{2, T}
    η::T # real a 2-D variable TODO: should be 2D
    θ::T
  end
end

# If this order is changed check update_aux!
function vars_aux(m::HBModel, T)
  @vars begin
    w::T
    pkin_reverse::T # ∫(-αᵀ θ) # TODO: remove me after better integral interface
    w_reverse::T               # TODO: remove me after better integral interface
    pkin::T         # ∫(-αᵀ θ)
    wz0::T          # w at z=0
    θʳ::T           # SST given    # TODO: Should be 2D
    f::T            # coriolis
    τ::T            # wind stress  # TODO: Should be 2D
    # κ::SMatrix{3, 3, T, 9} # diffusivity tensor (for convective adjustment)
    κᶻ::T
  end
end

function vars_gradient(m::HBModel, T)
  @vars begin
    u::SVector{2, T}
    θ::T
  end
end

function vars_diffusive(m::HBModel, T)
  @vars begin
    ∇u::SMatrix{3, 2, T, 6}
    ∇θ::SVector{3, T}
  end
end

function vars_integrals(m::HBModel, T)
  @vars begin
    ∇hu::T
    αᵀθ::T
  end
end

@inline function flux_nondiffusive!(m::HBModel, F::Grad, Q::Vars,
                                    A::Vars, t::Real)
  @inbounds begin
    u = Q.u # Horizontal components of velocity
    η = Q.η
    θ = Q.θ
    w = A.w   # vertical velocity
    pkin = A.pkin
    v = @SVector [u[1], u[2], w]
    Ih = @SMatrix [ 1 -0;
                   -0  1;
                   -0 -0]

    # ∇ • (u θ)
    F.θ += v * θ

    # ∇h • (g η)
    F.u += grav * η * Ih

    # ∇h • (- ∫(αᵀ θ))
    F.u += grav * pkin * Ih

    # ∇h • (v ⊗ u)
    # F.u += v * u'

  end

  return nothing
end

@inline wavespeed(m::HBModel, n⁻, _...) = abs(SVector(m.c₁, m.c₂, m.c₃)' * n⁻)

# We want not have jump penalties on η (since not a flux variable)
function update_penalty!(::Rusanov, ::HBModel, n⁻, λ, ΔQ::Vars,
                         Q⁻, A⁻, Q⁺, A⁺, t)
  ΔQ.η = -0

  #=
  θ⁻ = Q⁻.θ
  u⁻ = Q⁻.u
  w⁻ = A⁻.w
  @inbounds v⁻ = @SVector [u⁻[1], u⁻[2], w⁻]
  n̂_v⁻ = n⁻∘v⁻

  θ⁺ = Q⁺.θ
  u⁺ = Q⁺.u
  w⁺ = A⁺.w
  @inbounds v⁺ = @SVector [u⁺[1], u⁺[2], w⁺]
  n̂_v⁺ = n⁻∘v⁺

  # max velocity
  # n̂∘v = (abs(n̂∘v⁺) > abs(n̂∘v⁻) ? n̂∘v⁺ : n̂∘v⁻

  # average velocity
  n̂_v = (n̂_v⁻ + n̂_v⁺) / 2

  ΔQ.θ = ((n̂_v > 0) ? 1 : -1) * (n̂_v⁻ * θ⁻ - n̂_v⁺ * θ⁺)
  # ΔQ.θ = abs(n̂_v⁻) * θ⁻ - abs(n̂_v⁺) * θ⁺
  =#

  return nothing
end

@inline function flux_diffusive!(m::HBModel, F::Grad, Q::Vars, D::Vars,
                                 A::Vars, t::Real)
  ν = Diagonal(@SVector [m.νʰ, m.νʰ, m.νᶻ])
  F.u -= ν * D.∇u

  κ = Diagonal(@SVector [m.κʰ, m.κʰ, A.κᶻ])
  F.θ -= κ * D.∇θ

  return nothing
end

@inline function gradvariables!(m::HBModel, G::Vars, Q::Vars, A, t)
  G.u = Q.u
  G.θ = Q.θ

  return nothing
end

@inline function diffusive!(m::HBModel, D::Vars, G::Grad, Q::Vars,
                            A::Vars, t)
  D.∇u = G.u
  D.∇θ = G.θ

  return nothing
end

@inline function source!(m::HBModel{P}, source::Vars, Q::Vars, A::Vars,
                         t::Real) where P
  @inbounds begin
    u = Q.u # Horizontal components of velocity
    f = A.f
    wz0 = A.wz0

    # f × u
    source.u -= @SVector [-f * u[2], f * u[1]]

    source.η += wz0
  end

  return nothing
end

@inline function integrate_aux!(m::HBModel, integrand::Vars, Q::Vars, A::Vars)
  αᵀ = m.αᵀ
  integrand.αᵀθ = -αᵀ * Q.θ
  integrand.∇hu = A.w # borrow the w value from A...

  return nothing
end

function update_aux!(dg::DGModel, m::HBModel, Q::MPIStateArray, t::Real)
  MD = dg.modeldata

  # required to ensure that after integration velocity field is divergence free
  vert_filter = MD.vert_filter
  # Q[1] = u[1] = u, Q[2] = u[2] = v
  apply!(Q, (1, 2), dg.grid, vert_filter, VerticalDirection())

  exp_filter = MD.exp_filter
  # Q[4] = θ
  apply!(Q, (4,), dg.grid, exp_filter, VerticalDirection())

  return true
end

function update_aux_diffusive!(dg::DGModel, m::HBModel, Q::MPIStateArray, t::Real)
  A  = dg.auxstate

  # store ∇ʰu as integrand for w
  # update vertical diffusivity for convective adjustment
  function f!(::HBModel, Q, A, D, t)
    @inbounds begin
      A.w = -(D.∇u[1,1] + D.∇u[2,2])

      # κʰ = m.κʰ
      # A.κ = @SMatrix [κʰ -0 -0; -0 κʰ -0; -0 -0 κᶻ]
      D.∇θ[3] < 0 ? A.κᶻ = 1000 * m.κᶻ : A.κᶻ = m.κᶻ
    end

    return nothing
  end
  nodal_update_aux!(f!, dg, m, Q, t; diffusive=true)

  # compute integrals for w and pkin
  indefinite_stack_integral!(dg, m, Q, A, t) # bottom -> top
  reverse_indefinite_stack_integral!(dg, m, A, t) # top -> bottom

  # project w(z=0) down the stack
  # Need to be consistent with vars_aux
  # A[1] = w, A[5] = wz0
  copy_stack_field_down!(dg, m, A, 1, 5)

  return true
end

function ocean_init_aux! end
function init_aux!(m::HBModel, A::Vars, geom::LocalGeometry)
  return ocean_init_aux!(m, m.problem, A, geom)
end

function ocean_init_state! end
function init_state!(m::HBModel, Q::Vars, A::Vars, coords, t)
  return ocean_init_state!(m.problem, Q, A, coords, t)
end

@inline function boundary_state!(nf, m::HBModel, Q⁺::Vars, A⁺::Vars, n⁻,
                                 Q⁻::Vars, A⁻::Vars, bctype, t, _...)
  return ocean_boundary_state!(m, bctype, nf, Q⁺, A⁺, n⁻, Q⁻, A⁻, t)
end
@inline function boundary_state!(nf, m::HBModel,
                                 Q⁺::Vars, D⁺::Vars, A⁺::Vars,
                                 n⁻,
                                 Q⁻::Vars, D⁻::Vars, A⁻::Vars,
                                 bctype, t, _...)
  return ocean_boundary_state!(m, bctype, nf, Q⁺, D⁺, A⁺, n⁻, Q⁻, D⁻, A⁻, t)
end

@inline function ocean_boundary_state!(::HBModel, ::CoastlineFreeSlip,
                                       ::Union{Rusanov, CentralGradPenalty},
                                       Q⁺, A⁺, n⁻, Q⁻, A⁻, t)
  return nothing
end

@inline function ocean_boundary_state!(::HBModel, ::CoastlineFreeSlip,
                                       ::CentralNumericalFluxDiffusive, Q⁺,
                                       D⁺, A⁺, n⁻, Q⁻, D⁻, A⁻, t)
  D⁺.∇u = -D⁻.∇u

  D⁺.∇θ = -D⁻.∇θ

  return nothing
end

@inline function ocean_boundary_state!(::HBModel, ::CoastlineNoSlip,
                                       ::Union{Rusanov, CentralGradPenalty},
                                       Q⁺, A⁺, n⁻, Q⁻, A⁻, t)
  Q⁺.u = -Q⁻.u

  return nothing
end

@inline function ocean_boundary_state!(::HBModel, ::CoastlineNoSlip,
                                       ::CentralNumericalFluxDiffusive, Q⁺,
                                       D⁺, A⁺, n⁻, Q⁻, D⁻, A⁻, t)
  Q⁺.u = -Q⁻.u

  D⁺.∇θ = -D⁻.∇θ

  return nothing
end

@inline function ocean_boundary_state!(m::HBModel, ::OceanFloorFreeSlip,
                                       ::Union{Rusanov, CentralGradPenalty},
                                       Q⁺, A⁺, n⁻, Q⁻, A⁻, t)
  A⁺.w = -A⁻.w

  return nothing
end

@inline function ocean_boundary_state!(m::HBModel, ::OceanFloorFreeSlip,
                                       ::CentralNumericalFluxDiffusive, Q⁺,
                                       D⁺, A⁺, n⁻, Q⁻, D⁻, A⁻, t)
  A⁺.w = -A⁻.w
  D⁺.∇u = -D⁻.∇u

  D⁺.∇θ = -D⁻.∇θ

  return nothing
end

@inline function ocean_boundary_state!(m::HBModel, ::OceanFloorNoSlip,
                                       ::Union{Rusanov, CentralGradPenalty},
                                       Q⁺, A⁺, n⁻, Q⁻, A⁻, t)
  Q⁺.u = -Q⁻.u
  A⁺.w = -A⁻.w

  return nothing
end

@inline function ocean_boundary_state!(m::HBModel, ::OceanFloorNoSlip,
                                       ::CentralNumericalFluxDiffusive, Q⁺,
                                       D⁺, A⁺, n⁻, Q⁻, D⁻, A⁻, t)

  Q⁺.u = -Q⁻.u
  A⁺.w = -A⁻.w

  D⁺.∇θ = -D⁻.∇θ

  return nothing
end

@inline function ocean_boundary_state!(m::HBModel, ::Union{
                                       OceanSurfaceNoStressNoForcing,
                                       OceanSurfaceStressNoForcing,
                                       OceanSurfaceNoStressForcing,
                                       OceanSurfaceStressForcing},
                                       ::Union{Rusanov, CentralGradPenalty},
                                       Q⁺, A⁺, n⁻, Q⁻, A⁻, t)
  return nothing
end

@inline function ocean_boundary_state!(m::HBModel,
                                       ::OceanSurfaceNoStressNoForcing,
                                       ::CentralNumericalFluxDiffusive,
                                       Q⁺, D⁺, A⁺, n⁻, Q⁻, D⁻, A⁻, t)
  D⁺.∇u = -D⁻.∇u

  D⁺.∇θ = -D⁻.∇θ

  return nothing
end

@inline function ocean_boundary_state!(m::HBModel,
                                       ::OceanSurfaceStressNoForcing,
                                       ::CentralNumericalFluxDiffusive,
                                       Q⁺, D⁺, A⁺, n⁻, Q⁻, D⁻, A⁻, t)
  τ = A⁻.τ
  D⁺.∇u = -D⁻.∇u + 2 * @SMatrix [ -0 -0;
                                  -0 -0;
                                  τ / 1000 -0]

  D⁺.∇θ = -D⁻.∇θ

  return nothing
end

@inline function ocean_boundary_state!(m::HBModel,
                                       ::OceanSurfaceNoStressForcing,
                                       ::CentralNumericalFluxDiffusive,
                                       Q⁺, D⁺, A⁺, n⁻, Q⁻, D⁻, A⁻, t)
  D⁺.∇u = -D⁻.∇u

  θ  = Q⁻.θ
  θʳ = A⁻.θʳ
  λʳ = m.problem.λʳ
  D⁺.∇θ = -D⁻.∇θ + 2 * @SVector [-0, -0, λʳ * (θʳ - θ)]

  return nothing
end

@inline function ocean_boundary_state!(m::HBModel,
                                       ::OceanSurfaceStressForcing,
                                       ::CentralNumericalFluxDiffusive,
                                       Q⁺, D⁺, A⁺, n⁻, Q⁻, D⁻, A⁻, t)
  τ = A⁻.τ
  D⁺.∇u = -D⁻.∇u + 2 * @SMatrix [ -0 -0;
                                  -0 -0;
                                  τ / 1000 -0]

  θ  = Q⁻.θ
  θʳ = A⁻.θʳ
  λʳ = m.problem.λʳ
  D⁺.∇θ = -D⁻.∇θ + 2 * @SVector [-0, -0, λʳ * (θʳ - θ)]

  return nothing
end

end

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
const PyJuliaRawValue_Type__ref = Ref(PyNULL)
PyJuliaRawValue_Type() = begin
    ptr = PyJuliaRawValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaBaseValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.RawValue",
                base = base,
                repr = pyjlraw_repr,
                str = pyjlraw_str,
                getattro = pyjlraw_getattro,
                setattro = pyjlraw_setattro,
                call = pyjlraw_call,
                as_mapping = (
                    length = pyjlraw_length,
                    subscript = pyjlraw_getitem,
                    ass_subscript = pyjlraw_setitem,
                ),
                methods = [
                    (name = "__dir__", flags = Py_METH_NOARGS, meth = pyjlraw_dir),
                    (name = "_jl_any", flags = Py_METH_NOARGS, meth = pyjlraw_toany),
                ],
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaRawValue_Type__ref[] = ptr
    end
    ptr
end

PyJuliaRawValue_New(x) = PyJuliaValue_New(PyJuliaRawValue_Type(), x)

pyjlraw_repr(xo::PyPtr) =
    try
        x = PyJuliaValue_GetValue(xo)
        s = "<jl $(repr(x))>"
        PyUnicode_From(s)
    catch err
        PyErr_SetJuliaError(err)
        return PyNULL
    end

pyjlraw_str(xo::PyPtr) =
    try
        x = PyJuliaValue_GetValue(xo)
        s = string(x)
        PyUnicode_From(s)
    catch err
        PyErr_SetJuliaError(err)
        return PyNULL
    end

pyjl_attr_py2jl(k::String) = replace(k, r"_[b]+$" => (x -> "!"^(length(x) - 1)))
pyjl_attr_jl2py(k::String) = replace(k, r"!+$" => (x -> "_" * "b"^length(x)))

pyjlraw_getattro(xo::PyPtr, ko::PyPtr) = begin
    # Try generic lookup first
    ro = PyObject_GenericGetAttr(xo, ko)
    if isnull(ro) && PyErr_IsSet(PyExc_AttributeError())
        PyErr_Clear()
    else
        return ro
    end
    # Now try to get the corresponding property
    x = PyJuliaValue_GetValue(xo)
    k = PyUnicode_AsString(ko)
    isempty(k) && PyErr_IsSet() && return PyNULL
    k = pyjl_attr_py2jl(k)
    try
        v = getproperty(x, Symbol(k))
        PyJuliaRawValue_New(v)
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjlraw_setattro(xo::PyPtr, ko::PyPtr, vo::PyPtr) = begin
    # Try generic lookup first
    ro = PyObject_GenericSetAttr(xo, ko, vo)
    if ism1(ro) && PyErr_IsSet(PyExc_AttributeError())
        PyErr_Clear()
    else
        return ro
    end
    # Now try to set the corresponding property
    x = PyJuliaValue_GetValue(xo)
    k = PyUnicode_AsString(ko)
    isempty(k) && PyErr_IsSet() && return Cint(-1)
    k = pyjl_attr_py2jl(k)
    ism1(PyObject_Convert(vo, Any)) && return Cint(-1)
    v = takeresult(Any)
    try
        setproperty!(x, Symbol(k), v)
        Cint(0)
    catch err
        PyErr_SetJuliaError(err)
        Cint(-1)
    end
end

pyjlraw_dir(xo::PyPtr, _::PyPtr) = begin
    fo = PyObject_GetAttrString(PyJuliaBaseValue_Type(), "__dir__")
    isnull(fo) && return PyNULL
    ro = PyObject_CallNice(fo, PyObjectRef(xo))
    Py_DecRef(fo)
    isnull(ro) && return PyNULL
    x = PyJuliaValue_GetValue(xo)
    ks = try
        collect(map(string, propertynames(x)))
    catch err
        Py_DecRef(ro)
        PyErr_SetJuliaError(err)
        return PyNULL
    end
    for k in ks
        ko = PyUnicode_From(pyjl_attr_jl2py(k))
        isnull(ko) && (Py_DecRef(ro); return PyNULL)
        err = PyList_Append(ro, ko)
        Py_DecRef(ko)
        ism1(err) && (Py_DecRef(ro); return PyNULL)
    end
    return ro
end

pyjlraw_call(fo::PyPtr, argso::PyPtr, kwargso::PyPtr) = begin
    f = PyJuliaValue_GetValue(fo)
    if isnull(argso)
        args = Vector{Any}()
    else
        ism1(PyObject_Convert(argso, Vector{Any})) && return PyNULL
        args = takeresult(Vector{Any})
    end
    if isnull(kwargso)
        kwargs = Dict{Symbol,Any}()
    else
        ism1(PyObject_Convert(kwargso, Dict{Symbol,Any})) && return PyNULL
        kwargs = takeresult(Dict{Symbol,Any})
    end
    try
        x = f(args...; kwargs...)
        PyJuliaRawValue_New(x)
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjlraw_length(xo::PyPtr) =
    try
        x = PyJuliaValue_GetValue(xo)
        Py_ssize_t(length(x))
    catch err
        PyErr_SetJuliaError(err)
        Py_ssize_t(-1)
    end

pyjlraw_getitem(xo::PyPtr, ko::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)
    if PyTuple_Check(ko)
        ism1(PyObject_Convert(ko, Tuple)) && return PyNULL
        k = takeresult(Tuple)
        try
            PyJuliaRawValue_New(x[k...])
        catch err
            PyErr_SetJuliaError(err)
            PyNULL
        end
    else
        ism1(PyObject_Convert(ko, Any)) && return PyNULL
        k = takeresult(Any)
        try
            PyJuliaRawValue_New(x[k])
        catch err
            PyErr_SetJuliaError(err)
            PyNULL
        end
    end
end

pyjlraw_setitem(xo::PyPtr, ko::PyPtr, vo::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)
    ism1(PyObject_Convert(vo, Any)) && return PyNULL
    v = takeresult(Any)
    if PyTuple_Check(ko)
        ism1(PyObject_Convert(ko, Tuple)) && return PyNULL
        k = takeresult(Tuple)
        try
            x[k...] = v
            Cint(0)
        catch err
            PyErr_SetJuliaError(err)
            Cint(-1)
        end
    else
        ism1(PyObject_Convert(ko, Any)) && return PyNULL
        k = takeresult(Any)
        try
            x[k] = v
            Cint(0)
        catch err
            PyErr_SetJuliaError(err)
            Cint(-1)
        end
    end
end

pyjlraw_toany(xo::PyPtr, ::PyPtr) = PyJuliaValue_From(PyJuliaValue_GetValue(xo))


text:
"""
    setReplacementOperator(method::Function, varArgs...;
                           genj::GeneticJuliaStructure = GenJ,
                           needsComparison::Bool = false)

Sets the replacement operator of the population for `genj`.
This function constructs a structure of type `ReplacementOperator` and adds it
to the main structure.

# Arguments
- `method::Function`: method that will tell which individuals of the population and
    the offspring will remain for the next generation.
- `varArgs`: list of arguments for the replacement method. See `"Examples"`.

# Keyword Arguments
- `genj::GeneticJuliaStructure = GenJ`: the main structure. For `code user`: don't modify
    the default value unless multiple experiments are going to be run.
- `needsComparison::Bool = false`: tells wether the user wants the system to compare
    between the individuals in order to select the ones that will replace the old
    population.
- `eliteSize::Integer`: number of fittest individuals to be maintained if there
    is elitism.

# Examples
```jldoctest
julia> function ExampleReplacement(population::Array{Individual}, offspring::Array{Individual}, arg1, arg2)
    #···Code···
    return selectedIndexesFromPopulation, selectedIndexesFromOffspring
end
ExampleReplacement (generic function with 1 method)

julia> setReplacementOperator(ExampleReplacement, arg1, arg2)
ReplacementOperator(false, 0, ExampleReplacement, Any[arg1, arg2])
```

```jldoctest
julia> function ExampleReplacement(population::Array{Individual}, offspring::Array{Individual}, compare::Function, arg1, arg2)
    #···Code···
    return selectedIndexesFromPopulation, selectedIndexesFromOffspring
end
ExampleReplacement (generic function with 1 method)

julia> setReplacementOperator(ExampleReplacement, arg1, arg2, needsComparison=true, eliteSize=1)
ReplacementOperator(true, 1, ExampleReplacement, Any[arg1, arg2])
```

See also: [`ReplacementOperator`](@ref)
"""
function setReplacementOperator(method::Function, varArgs...;
                                genj::GeneticJuliaStructure = GenJ,
                                needsComparison::Bool = false,
                                eliteSize::Integer = 0)

    varArgs = [i for i in varArgs]
    genj._replacementOp = ReplacementOperator(needsComparison, eliteSize, method, varArgs)
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
"""
get_affine_map(x, ξ)
Get (A,b) such that ξ = A * x + b is the affine
mapping from the simplex with vertices x to the simplex of vertices ξ.
"""
function get_affine_map(x, ξ)
    @assert length(x[1]) == length(ξ[1]) "dimension mismatch"
    n = length(x[1])
    T = eltype(x[1])

    B = zeros(T, n*(n+1), n*(n+1))
    L = zeros(T, n*(n+1))

    for i in 1:length(x)
        for j in 1:n
            row = (i-1) * n + j
            B[row, n*(j-1)+1:n*j] = x[i]
            L[row] = ξ[i][j]
            B[row, n * n + j] = one(T)
        end
    end
    X = B\L
    return reshape(X[1:n*n], (n, n)), X[n*n+1:end]
end

import Base.ht_keyindex2!

mutable struct ScalarWrapper{T}
    x::T
end

@inline Base.getindex(s::ScalarWrapper) = s.x
@inline Base.setindex!(s::ScalarWrapper, v) = s.x = v
Base.copy(s::ScalarWrapper{T}) where {T} = ScalarWrapper{T}(copy(s.x))


text:
module StructuringElements

export strel, strel_type, strel_size, strel_ndims
export SEMask, MorphologySE, SEOffset
export is_symmetric, require_symmetric_strel
export strel_split
export strel_box, SEBox, SEBoxArray
export strel_diamond, SEDiamond, SEDiamondArray
export strel_chain, strel_product, SEChain, SEChainArray

# TODO(johnnychen94): remove ImageCore dependency
using ImageCore: coords_spatial

using OffsetArrays
using OffsetArrays: centered

abstract type MorphologySE{N} end
abstract type MorphologySEArray{N} <: AbstractArray{Bool,N} end

OffsetArrays.centered(A::MorphologySEArray) = A

"""
    SEMask{N}()

A (holy) trait type for representing structuring element as connectivity mask. This
connectivity mask SE is a bool array where `true` indicates that pixel position is connected
to the center point.

```jldoctest; setup=:(using ImageMorphology)
julia> se = centered(Bool[0 1 0; 1 1 1; 0 1 0]) # commonly known as C4 connectivity
3×3 OffsetArray(::Matrix{Bool}, -1:1, -1:1) with eltype Bool with indices -1:1×-1:1:
 0  1  0
 1  1  1
 0  1  0

julia> strel_type(se)
ImageMorphology.StructuringElements.SEMask{2}()
```

See also [`SEOffset`](@ref StructuringElements.SEOffset) for the displacement offset
representation. More details can be found on he documentation page [Structuring
Element](@ref concept_se).
"""
struct SEMask{N} <: MorphologySE{N} end

"""
    SEOffset{N}()

A (holy) trait type for representing structuring element as displacement offsets. This
displacement offsets SE is an array of `CartesianIndex` where each element stores the
displacement offset from the center point.

```jldoctest; setup=:(using ImageMorphology)
julia> se = [CartesianIndex(-1, 0), CartesianIndex(0, -1), CartesianIndex(1, 0), CartesianIndex(0, 1)]
4-element Vector{CartesianIndex{2}}:
 CartesianIndex(-1, 0)
 CartesianIndex(0, -1)
 CartesianIndex(1, 0)
 CartesianIndex(0, 1)

julia> strel_type(se)
ImageMorphology.StructuringElements.SEOffset{2}()
```

See also [`SEMask`](@ref StructuringElements.SEMask) for the connectivity mask representation.
More details can be found on he documentation page [Structuring Element](@ref concept_se).
"""
struct SEOffset{N} <: MorphologySE{N} end

include("strel.jl")
# SE constructors
include("strel_chain.jl")
include("strel_box.jl")
include("strel_diamond.jl")

include("utils.jl")

end


text:
using VectorizedRNG: local_rng

@inline function Random.rand!(A::AbstractStrideArray, args::Vararg{Any,K}) where {K}
    rand!(local_rng(), A, args...)
end
@inline function Random.randn!(A::AbstractStrideArray, args::Vararg{Any,K}) where {K}
    randn!(local_rng(), A, args...)
end
## ignore type...
@inline Random.rand!(A::AbstractStrideArray, ::Type{T}) where {T} = rand!(local_rng(), A)

function rand_expr(expr, args...)
    # @show expr.args
    N = length(expr.args)
    n = 2
    randtypes = (:Float32,:Float64,:Int,:Int32,:Int64,:UInt,:UInt32,:UInt64)
    arg2 = expr.args[2]
    if arg2 ∈ randtypes
        T = arg2
        n += 1
    elseif arg2 isa Expr && arg2.head === :($)
        T = esc(first(arg2.args))
        n += 1
    else
        T = Float64
    end
    s = Expr(:tuple)
    for i in n:N
        aᵢ = expr.args[i]
        if aᵢ isa Integer
            push!(s.args, static_expr(Int(aᵢ)))
        else
            push!(s.args, esc(aᵢ))
        end
    end
    f! = Symbol(expr.args[1], :!)
    array = :(StrideArray{$T}(undef, $s))
    # we call PtrArray so that `rand!` doesn't force heap allocation.
    call = Expr(:call, f!, Expr(:call, :PtrArray, :array))
    for arg ∈ args
        push!(call.args, esc(arg))
    end
    quote
        array = $array
        b = array.data
        # we add GC.@preserve to save the array.
        GC.@preserve b $call
        array
    end
end

macro StrideArray(expr, args...)
    rand_expr(expr, args...)
end


text:
in_channel = 10
out_channel = 5
num_V = 6
num_E = 7
adj = [0. 1. 0. 0. 0. 0.;
       1. 0. 0. 1. 1. 1.;
       0. 0. 0. 0. 0. 1.;
       0. 1. 0. 0. 1. 0.;
       0. 1. 0. 1. 0. 1.;
       0. 1. 1. 0. 1. 0.]
ne = [[2], [1,4,5,6], [6], [2,5], [2,4,6], [2,3,5]]

struct NewLayer <: MessagePassing
    weight
end
NewLayer(m, n) = NewLayer(randn(m,n))

(l::NewLayer)(fg) = propagate(l, fg, :add)

X = Array(reshape(1:num_V*in_channel, in_channel, num_V))
fg = FeaturedGraph(adj, X)

l = NewLayer(out_channel, in_channel)

@testset "msgpass" begin
    @testset "no message or update" begin
        fg_ = l(fg)

        @test graph(fg_) === adj
        @test size(node_feature(fg_)) == (in_channel, num_V)
        @test size(edge_feature(fg_)) == (in_channel, 2*num_E)
        @test size(global_feature(fg_)) == (0,)
    end

    GeometricFlux.message(l::NewLayer, x_i, x_j, e_ij) = l.weight * x_j
    @testset "message function" begin
        fg_ = l(fg)

        @test graph(fg_) === adj
        @test size(node_feature(fg_)) == (out_channel, num_V)
        @test size(edge_feature(fg_)) == (out_channel, 2*num_E)
        @test size(global_feature(fg_)) == (0,)
    end

    GeometricFlux.update(l::NewLayer, m, x) = l.weight * x + m
    @testset "message and update" begin
        fg_ = l(fg)

        @test graph(fg_) === adj
        @test size(node_feature(fg_)) == (out_channel, num_V)
        @test size(edge_feature(fg_)) == (out_channel, 2*num_E)
        @test size(global_feature(fg_)) == (0,)
    end
end


text:
################################################################################
## addition

macro ejectsolutions()
  esc(quote
    (isnan(lower_result) || isnan(upper_result)) && return nan(U)
    (upper_result == lower_result) && return upper_result
  end)
end

@universal function add(a::Ubound, b::Unum)
  lower_result = resolve_lower(a.lower + b)
  upper_result = resolve_upper(a.upper + b)

  @ejectsolutions

  (is_ulp(lower_result) && is_ulp(upper_result)) && return resolve_as_utype!(lower_result, upper_result)
  return B(lower_result, upper_result)
end

@universal add(a::Unum, b::Ubound) = add(b, a)

@universal function add(a::Ubound, b::Ubound)
  lower_result = resolve_lower(a.lower + b.lower)
  upper_result = resolve_upper(a.upper + b.upper)

  @ejectsolutions

  (is_ulp(lower_result) && is_ulp(upper_result)) ? resolve_as_utype!(lower_result, upper_result) : B(lower_result, upper_result)
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
@testset "air_density" begin
  ad = @inferred air_density(25.0,100.0) # Tair, pressure
  # regression test
  @test ≈(ad, 1.168, atol =0.001)
end

@testset "pressure_from_elevation" begin
  pressure = @inferred pressure_from_elevation(500.0, 25.0) # elev, Tair
  # regression test
  @test ≈(pressure, 95.681, atol =0.001)
  #
  pressure = @inferred pressure_from_elevation(500.0, 25.0, 1.5) # elev, Tair, VPD
  # regression test
  @test ≈(pressure, 95.717, atol =0.001)
end

@testset "virtual_temp" begin
  Tair,pressure,VPD = 25.0,100.0,1.5
  vt = @inferred virtual_temp(Tair,pressure,VPD)
  # regression test
  @test ≈(vt, 26.9, atol =0.1)
end

@testset "kinematic_viscosity" begin
  Tair,pressure = 25.0,100.0
  vis = @inferred kinematic_viscosity(Tair,pressure)
  # regression test
  @test ≈(vis, 1.58e-5, atol =1e-7)
end

@testset "air_density" begin
  ad = @inferred air_density(25.0,100.0) # Tair, pressure
  # regression test
  @test ≈(ad, 1.168, atol =0.001)
end

@testset "dew_point" begin
  VPD = 1.5
  Tair = 25.0
  accuracy = 1e-2
  #@descend_code_warntype VPD_to_e(VPD,Tair)
  ea = @inferred VPD_to_e(VPD,Tair)
  Td = @inferred dew_point_from_e(ea,accuracy = accuracy)
  @test ≈(ea, Esat_from_Tair(Td), atol = accuracy)
  @test ≈(dew_point(Tair, VPD; accuracy = accuracy), Td, atol = accuracy)
  Td1 = @test_logs (:warn,r"set to 1")  @inferred dew_point(Tair, VPD; accuracy = 1.2)
  @test ≈(Td1, Td, atol = 1.0)
end

@testset "wetbulb temperature" begin
  Tair = 25.0
  VPD = 1.0
  pressure = 100.0
  accuracy = 1e-2
  #
  gamma  = @inferred psychrometric_constant(Tair,pressure)
  ea     = VPD_to_e(VPD,Tair)
  Tw = @inferred wetbulb_temp_from_e_Tair_gamma(ea, Tair, gamma; accuracy = accuracy)
  #@test ≈(ea, Esat_from_Tair(Tw) - gamma* (Tair - Tw), atol = accuracy)
  @test ≈(ea, Esat_from_Tair(Tw) - 0.93*gamma* (Tair - Tw), atol = accuracy)
  @test ≈(wetbulb_temp(Tair, pressure, VPD; accuracy = accuracy), Tw, atol = accuracy)
  Tw1 = @test_logs (:warn,r"set to 1") @inferred wetbulb_temp(Tair, pressure, VPD; accuracy = 1.2)
  @test ≈(Tw1, Tw, atol = 1.0)
end


text:
# Solve a simple 1D diffusion problem (without boundary conditions)
#
# ARGS[1] is the size of the system
#
# This script can be run with an arbitrary number of processes:
# julia ./exKSP.jl 40
# or
# mpirun -np 2 julia ./exKSP.jl 40
# or
# mpirun -np 4 julia ./exKSP.jl 40
#
# and so on

using MPI
using PETSc2
PetscInitialize()

npts = parse(Int, ARGS[1])

myrank = MPI.Comm_rank(MPI.COMM_WORLD)
commsize = MPI.Comm_size(MPI.COMM_WORLD)

opts = Dict{String, String}(
  "-ksp_monitor" => "",
  "-malloc_dump" => "",
  "-pc_type" => "bjacobi",
  "-sub_pc_type" => "ilu",
  "-sub_pc_factor_levels" => "0",
  "-ksp_gmres_modifiedgramschmidt" => "",
)

PetscSetOptions(opts)

stencil = PetscScalar[1.0, -2.0, 1.0]
n_off_diag = div(length(stencil), 2)  # number of points off the main diagonal

# create matrix
A = PetscMat(npts, npts, MATMPIAIJ, MPI.COMM_WORLD)

# create the vectors and populate the right hand side
x = PetscVec(npts, VECMPI, MPI.COMM_WORLD)
b = PetscVec(npts, VECMPI, MPI.COMM_WORLD)

# preallocate A (this is a slight overestimate)
lcl_indices = local_indices(b)
npts_local = length(lcl_indices)

d_nnz = zeros(PetscInt, npts_local)
o_nnz = zeros(PetscInt, npts_local)

if commsize == 1
  fill!(d_nnz, length(stencil))
else
  o_nnz[1] = n_off_diag + 1
  fill!(d_nnz, length(stencil))
  o_nnz[end] = n_off_diag + 1
end

MatMPIAIJSetPreallocation(A, 0, d_nnz, 0, o_nnz)

# populate the matrix
idxm = PetscInt[0]
idxn = zeros(PetscInt, length(stencil))
if myrank == 0
  # do the first row
  idxm[1] = 1
  fill!(idxn, -1)  # negative entries are ignored
  start_idx = n_off_diag + 1
  idx = 1
  for i = start_idx:length(stencil)
    idxn[i] = idx
    idx += 1
  end

  set_values1!(A, idxm, idxn, stencil, ADD_VALUES)
end

if myrank == (commsize-1)
  # do the last row
  idxm[1] = npts
  fill!(idxn, -1)

  idx = npts - (length(stencil) - n_off_diag) + 1
  for i=1:(length(stencil) - n_off_diag)
    idxn[i] = idx
    idx += 1
  end

  set_values1!(A, idxm, idxn, stencil, ADD_VALUES)
end

# do everything else
first_row = lcl_indices[1]
last_row = lcl_indices[end]
if myrank == 0
  first_row += n_off_diag
end
if myrank == (commsize - 1)
  last_row -= n_off_diag
end

for i = first_row:last_row
  idxm[1] = i
  idxn[1] = i-1
  idxn[2] = i
  idxn[3] = i+1
  set_values1!(A, idxm, idxn, stencil, ADD_VALUES)
end

# start assembly
assembly_begin(A, MAT_FINAL_ASSEMBLY)

# populate right hand side
b_tmp = VecGetArray(b)
for i=1:length(b_tmp)
  b_tmp[i] = lcl_indices[i]  # arbitrary
end
VecRestoreArray(b, b_tmp)

# create the ksp solver
ksp = KSP(MPI.COMM_WORLD)
SetFromOptions(ksp)
SetOperators(ksp, A, A)
SetUp(ksp)
SetTolerances(ksp, 1e-12, 1e-12, PETSc2.PETSC_DEFAULT, PETSc2.PETSC_DEFAULT)

# end matrix assembly
assembly_end(A, MAT_FINAL_ASSEMBLY)

# view the matrix
PetscView(A)

# solve the system
KSPSolve(ksp, b, x)

# view the solution
PetscView(x)

PetscDestroy(A)
PetscDestroy(x)
PetscDestroy(b)
PetscDestroy(ksp)
PetscFinalize()


text:
using SpikingNetworks
parseParameters"""
τ=20
σ=0.015"""
parseEquations"""
dv/dt=(1.1-v+σ*ξ(t))/τ"""
net=init"""
N=25
v=rand(N)"""
spike"""
v.>1
5"""
reset"v[i]=0"
net.dt=0.1
run!(net,500/net.dt)

using MatlabPlot
t,neuron=rasterPlot(net)
mplot(t,neuron,".k")


text:
## To Lower Case ##

rem{T<:BioUnit}(x::T, y::Type{UInt8}) = rem(UInt8(x), y)

function lowercase!{T<:BioUnit}(s::AbstractArray{T})
  for i = 1:length(s)
    if 'A' <= Char(s[i]) <= 'Z'
      s[i] += 32
    end
  end
  return s
end

lowercase{T<:BioUnit}(s::AbstractArray{T}) = lowercase!(copy(s))

## To Upper Case ##

function uppercase!{T<:BioUnit}(s::AbstractArray{T})
  for i = 1:length(s)
    if 'a' <= Char(s[i]) <= 'z'
      s[i] -= 32
    end
  end
  return s
end

uppercase{T<:BioUnit}(s::AbstractArray{T}) = uppercase!(copy(s))

## Upper and Lower for BioUnit ##

uppercase{T<:BioUnit}(x::T) = convert(T,uppercase(Char(x)))
lowercase{T<:BioUnit}(x::T) = convert(T,lowercase(Char(x))) # Faster but unsafe: convert( T , x | 32 )

## Matching ##

# _unsafe_ascii makes a faster conversion, but violates immutability
_unsafe_ascii{T<:BioUnit}(seq::Vector{T}) = ASCIIString( reinterpret(UInt8,seq) )

ismatch{T<:BioUnit}(r::Regex, s::Vector{T}) = ismatch(r, _unsafe_ascii(s))

match{T<:BioUnit}(re::Regex, s::Vector{T}, idx::Integer) = match(re,_unsafe_ascii(s),idx)
match{T<:BioUnit}(r::Regex, s::Vector{T}) = match(r, _unsafe_ascii(s), 1)

search{T<:BioUnit}(s::Vector{T}, r::Regex, idx::Integer) = search(_unsafe_ascii(s),r,idx)
search{T<:BioUnit}(s::Vector{T}, r::Regex) = search(_unsafe_ascii(s),r,1)

search{T<:BioUnit}(s::Vector{T}, c, i) = search(_unsafe_ascii(s),c,i)
search{T<:BioUnit}(s::Vector{T}, c) = search(_unsafe_ascii(s),c,1)

eachmatch{T<:BioUnit}(re::Regex, s::Vector{T}, ovr::Bool) = Base.RegexMatchIterator(re,_unsafe_ascii(s),ovr)
eachmatch{T<:BioUnit}(re::Regex, s::Vector{T}) = Base.RegexMatchIterator(re,_unsafe_ascii(s),false)

## Some String Function ##

repeat{T<:BioUnit}(s::Vector{T}, n::Int) = vec(repmat(s,n,1))

repeat{T<:BioUnit}(x::T, n::Int) = fill(x,n)

function swap!{T<:BioUnit,Tp<:Union{Integer, Char},Tc<:Union{Integer, Char}}(x::AbstractArray{T},p::Tp,c::Tc)
  pat = convert(T,p)
  cha = convert(T,c)
  for i in 1:length(x)
    if x[i] == pat
      x[i] = cha
    end
  end
  return x
end

swap{T<:BioUnit,Tp<:Union{Integer, Char},Tc<:Union{Integer, Char}}(x::AbstractArray{T},p::Tp,c::Tc) = swap!(copy(x),p,c)

function replace{T<:BioUnit}(s::Vector{T}, pattern, r, limit::Integer)
    cases = 0
    range = search(s,pattern)
    st = range.start
    en = range.stop
    if st==0
      return(s)
    else
      res = T[]
      pre = 1
    end
    while st != 0
	if cases == limit; break; end
        if (st-1) >= pre
	  append!(res,s[pre:(st-1)])
	end
	isa(r,Vector) ? append!(res,r) : push!(res,r)
	pre = en + 1
        range = search(s,pattern,pre)
	st = range.start
	en = range.stop
        cases += 1
    end
    if pre<=length(s); append!(res,s[pre:end]) end
    res
end

replace{T<:BioUnit}(s::Vector{T}, pat, r) = replace(s, pat, r, -1)


text:
using ArgParse

struct Commandline
    command_ok::Bool                 # Boolean the check if the line of command is valid
    cpu_time::Int                    # Allocated CPU time (defaults to 5min)
    seed::Int                        # Random seed (defaults to 0, in this case the current time value will be used as seed)
    instance_path::String            # Instance path
    output_path::String              # Output path
end

# Constructor of struct Commandline
function Commandline()
    """
    Read command line arguments.
    """

    s = ArgParseSettings()

    @add_arg_table s begin
        "instance_path"
            help = "Path to instance."
            arg_type = String
            required = true
        "--cpu_time","-t"
            help = "CPU time in seconds (defaults to 300s)."
            arg_type = Int
            default = 300
        "--sol"
            help = "File where to output the solution statistics (defaults to the instance file name prepended with 'sol-')."
            arg_type = String
        "--seed"
            help = "Random seed (defaults to 0, in this case the current time value will be used as seed)"
            arg_type = Int
            default = 0
    end

    args = parse_args(s)

    # If no output path is specified then save solution in "outputs" folder with the same name as the instance
    if args["sol"] == nothing
        instance_name = split(args["instance_path"],'/')[end]
        args["sol"] = "../outputs/sol-"*instance_name
    end

    c = Commandline(true, args["cpu_time"], args["seed"], args["instance_path"], args["sol"])
    return c
end

max_stars_repo_name: BSAraujo/QAP-Metaheuristics

max_issues_repo_name: BSAraujo/QAP-Metaheuristics

max_forks_repo_name: BSAraujo/QAP-Metaheuristics


text:
using SyntheticNetworks
using EmbeddedGraphs
using LightGraphs
using SimpleTraits
using Statistics

# defining parameters for graph initialization

n = 20
n0 = 10
p = 0.2
q = 0.3
r = 1. / 3.
s = 0.1
u = 0.0 # currently not supported

function synthetic_network(n, n0, p, q, r, s, u)
    RPG = RandomPowerGrid(n, n0, p, q, r, s, u)
    network = generate_graph(RPG) # network is an EmbeddedGraph object
    return network.graph
end

g = synthetic_network(n, n0, p, q, r, s, u)

## statistics

# average node degree
function average_node_degree(g::SimpleGraph)
    sum(degree(g))/nv(g)
end

average_node_degree(g)

# global clustering coefficient
global_clustering_coefficient(g)

# average shortest path length
mean(floyd_warshall_shortest_paths(g).dists)

## assortativity

function assortativity(g::AbstractGraph{T}) where T
    P = promote_type(Int64, T) # at least Int64 to reduce risk of overflow
    nue  = ne(g)
    sjk = sj = sk = sjs = sks = zero(P)
    for d in edges(g)
        j = P(outdegree(g, src(d)) - 1)
        k = P(indegree(g, dst(d)) - 1)
        sjk += j*k
        sj  += j
        sk  += k
        sjs += j^2
        sks += k^2
    end
    return assortativity_coefficient(g, sjk, sj, sk, sjs, sks, nue)
end

#=
assortativity coefficient for directed graphs:
see equation (21) in M. E. J. Newman: Mixing patterns in networks, Phys. Rev. E 67, 026126 (2003),
http://arxiv.org/abs/cond-mat/0209450
=#
@traitfn function assortativity_coefficient(g::::IsDirected, sjk, sj, sk, sjs, sks, nue)
    return (sjk - sj*sk/nue) / sqrt((sjs - sj^2/nue)*(sks - sk^2/nue))
end

#=
assortativity coefficient for undirected graphs:
see equation (4) in M. E. J. Newman: Assortative mixing in networks, Phys. Rev. Lett. 89, 208701 (2002),
http://arxiv.org/abs/cond-mat/0205405/
=#
@traitfn function assortativity_coefficient(g::::(!IsDirected), sjk, sj, sk, sjs, sks, nue)
    return (sjk/nue - ((sj + sk)/(2*nue))^2) / ((sjs + sks)/(2*nue) - ((sj + sk)/(2*nue))^2)
end

assortativity(g)

# end
##

# for plotting

open("julia_networks.txt", "a") do io
    for i in 1:1000
        graph = synthetic_network(n, n0, p, q, r, s, u)
        node_degree = average_node_degree(graph)
        assortativity_value = assortativity(graph)
        average_shortes_path_length = mean(floyd_warshall_shortest_paths(graph).dists)
        clust_coefficient = global_clustering_coefficient(graph)
        write(io, string(node_degree), "\t", string(assortativity_value), "\t", string(average_shortes_path_length), "\t", string(clust_coefficient), "\n")
        #write(io, string(i))
    end
end


text:
using MPIReco

# multi-gradient is a special case of multi-patch
@testset "multi-gradient in-memory reconstruction" begin
  # low gradient with patch no 3
  b = MultiMPIFile([joinpath(datadir, "measurements", "20211226_203921_MultiGradient", "1.mdf"),
                    joinpath(datadir, "measurements", "20211226_203921_MultiGradient", "4.mdf")])
  bSFs = MultiMPIFile([joinpath(datadir, "calibrations", "13.mdf"), joinpath(datadir, "calibrations", "14.mdf")])
  names = names = (:color, :x, :y, :z, :time)

  c1 = reconstruction(bSFs, b;
			    SNRThresh=2, frames=1, lambd=0.003, minFreq=80e3,
			    recChannels=1:2,iterations=3, roundPatches=false)
  @test axisnames(c1) == names
  @test axisvalues(c1) == (1:1, -26.0u"mm":1.0u"mm":26.0u"mm", -26.75u"mm":1.0u"mm":26.25u"mm", 0.0u"mm":1.0u"mm":0.0u"mm", 0.0u"ms":0.6528u"ms":0.0u"ms")
  im1 = reverse(c1[1,:,:,1,1]',dims=1)
  exportImage(joinpath(imgdir, "MultiGradient1.png"), im1)

  #Iabs = abs.(im1)
  #Iabs[Iabs .< 0.2*maximum(Iabs)] .= 0
  #Icolored = colorview(Gray, Iabs./maximum(Iabs))
  #FileIO.save("$imgdir/MultiGradient1.png", Icolored )

  # low gradient with all 3 high gradient patches
  dirs = ["1.mdf", "2.mdf", "3.mdf", "4.mdf"]
  b = MultiMPIFile(joinpath.(datadir, "measurements", "20211226_203921_MultiGradient",  dirs))

  c2 = reconstruction(bSFs, b;
			    SNRThresh=2, frames=1, lambd=0.003, minFreq=80e3,
			    recChannels=1:2,iterations=3, roundPatches=false)
  @test axisnames(c2) == names
  # TODO wo kommt dieser Floatingpoint Fehler rein? Lässt der sich vermeiden?
  @test axisvalues(c2) == (1:1, -26.0u"mm":1.0u"mm":26.0u"mm", -26.500000000000004u"mm":1.0u"mm":30.499999999999996u"mm", 0.0u"mm":1.0u"mm":0.0u"mm", 0.0u"ms":0.6528u"ms":0.0u"ms")
  im2 = reverse(c2[1,:,:,1,1]',dims=1)
  exportImage(joinpath(imgdir, "MultiGradient2.png"), im2)

  #Iabs = abs.(im2)
  #Iabs[Iabs .< 0.2*maximum(Iabs)] .= 0
  #Icolored = colorview(Gray, Iabs./maximum(Iabs))
  #FileIO.save("$imgdir/MultiGradient2.png", Icolored )

  # low gradient only
  b = MultiMPIFile([joinpath(datadir, "measurements", "20211226_203921_MultiGradient", "1.mdf")])
  bSFs = MultiMPIFile([joinpath(datadir, "calibrations", "13.mdf")])

  c3 = reconstruction(bSFs, b;
			    SNRThresh=2, frames=1, lambd=0.003, minFreq=80e3,
			    recChannels=1:2,iterations=3, roundPatches=false)
  im3 = reverse(c3[1,:,:,1,1]',dims=1)
  exportImage(joinpath(imgdir, "MultiGradient3.png"), im3)

  #Iabs = abs.(im3)
  #Iabs[Iabs .< 0.2*maximum(Iabs)] .= 0
  #Icolored = colorview(Gray, Iabs./maximum(Iabs))
  #FileIO.save("$imgdir/MultiGradient3.png", Icolored )
end


text:
#=
        Constants for the problem, just change nb_teams, the others depend on it
=#
const nb_teams = 10
const nb_weeks = nb_teams - 1
const nb_periods = div(nb_teams, 2)

# For better readability, we will use a tables constaining the variables
nb_var = (nb_teams * nb_weeks * nb_periods)
variables = reshape(1:nb_var, nb_teams, nb_weeks, nb_periods)

"""
    function read_results(file_name::String)

    Will displays the results contained in the file file_name, in a readable way
"""
function read_results(file_name::String)
    results = Vector{Int}(undef,nb_var)
	file = open(file_name, "r")

	lines = readlines(file)
	if lines[1] == "SAT" || lines[1] == "sat"
		lines = split(lines[2]," ")
		for i in 1:nb_var
			results[i]=parse(Int,lines[i])
		end
    end
    tables = reshape(results, nb_teams, nb_weeks, nb_periods)
    for w in 1:nb_weeks
        println("WEEK n°", w)
        for t in 1:nb_teams
            for p in 1:nb_periods
                if tables[t,w,p] == variables[t,w,p]
                    println("Team n°", t, " play during period : ", p)
                end
            end
        end
        println()
        println()
    end
    close(file)
end

read_results("res.res")

max_stars_repo_name: remi-garcia/ProjetCSP

max_issues_repo_name: remi-garcia/ProjetCSP

max_forks_repo_name: remi-garcia/ProjetCSP


text:
using Documenter, JuliaReachDevDocs

DocMeta.setdocmeta!(JuliaReachDevDocs, :DocTestSetup,
                    :(using JuliaReachDevDocs); recursive=true)

makedocs(;
    sitename = "JuliaReachDevDocs",
    authors = "Marcelo Forets, Christian Schilling",
    modules = [JuliaReachDevDocs],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets = ["assets/juliareach.css"]),
    pages = [
        "Home" => "index.md",
        "Adding a new package" => "newpkg.md",
        "Git workflow" => "workflow.md",
        "Releases and bots" => "release.md",
        "Coding guidelines" => "guidelines.md",
        "About" => "about.md"
    ],
    strict = true
)

deploydocs(
    repo = "github.com/JuliaReach/JuliaReachDevDocs",
    push_preview = true
)

max_stars_repo_name: JuliaReach/JuliaReachDevDocs

max_issues_repo_name: JuliaReach/JuliaReachDevDocs

max_forks_repo_name: JuliaReach/JuliaReachDevDocs


text:
module Manifolds

import ManifoldsBase:
    @trait_function,
    _access_nested,
    _get_basis,
    _injectivity_radius,
    _inverse_retract,
    _inverse_retract!,
    _read,
    _retract,
    _retract!,
    _write,
    active_traits,
    allocate,
    allocate_coordinates,
    allocate_result,
    allocate_result_type,
    allocation_promotion_function,
    array_value,
    base_manifold,
    check_point,
    check_size,
    check_vector,
    copy,
    copyto!,
    default_inverse_retraction_method,
    default_retraction_method,
    default_vector_transport_method,
    decorated_manifold,
    distance,
    dual_basis,
    embed,
    embed!,
    exp,
    exp!,
    get_basis,
    get_basis_default,
    get_basis_diagonalizing,
    get_basis_orthogonal,
    get_basis_orthonormal,
    get_basis_vee,
    get_component,
    get_coordinates,
    get_coordinates!,
    get_coordinates_diagonalizing,
    get_coordinates_diagonalizing!,
    get_coordinates_orthogonal,
    get_coordinates_orthonormal,
    get_coordinates_orthogonal!,
    get_coordinates_orthonormal!,
    get_coordinates_vee!,
    get_embedding,
    get_iterator,
    get_vector,
    get_vector!,
    get_vector_diagonalizing,
    get_vector_diagonalizing!,
    get_vector_orthogonal,
    get_vector_orthonormal,
    get_vector_orthogonal!,
    get_vector_orthonormal!,
    get_vectors,
    gram_schmidt,
    hat,
    hat!,
    injectivity_radius,
    _injectivity_radius,
    injectivity_radius_exp,
    inner,
    isapprox,
    is_point,
    is_vector,
    inverse_retract,
    inverse_retract!,
    _inverse_retract,
    _inverse_retract!,
    inverse_retract_caley!,
    inverse_retract_embedded!,
    inverse_retract_nlsolve!,
    inverse_retract_pade!,
    inverse_retract_polar!,
    inverse_retract_project!,
    inverse_retract_qr!,
    inverse_retract_softmax!,
    log,
    log!,
    manifold_dimension,
    mid_point,
    mid_point!,
    norm,
    number_eltype,
    number_of_coordinates,
    parallel_transport_along,
    parallel_transport_along!,
    parallel_transport_direction,
    parallel_transport_direction!,
    parallel_transport_to,
    parallel_transport_to!,
    parent_trait,
    power_dimensions,
    project,
    project!,
    representation_size,
    retract,
    retract!,
    retract_caley!,
    retract_exp_ode!,
    retract_pade!,
    retract_polar!,
    retract_project!,
    retract_qr!,
    retract_softmax!,
    set_component!,
    vector_space_dimension,
    vector_transport_along, # just specified in Euclidean - the next 5 as well
    vector_transport_along_diff,
    vector_transport_along_project,
    vector_transport_along!,
    vector_transport_along_diff!,
    vector_transport_along_project!,
    vector_transport_direction,
    vector_transport_direction_diff,
    vector_transport_direction!,
    vector_transport_direction_diff!,
    vector_transport_to,
    vector_transport_to_diff,
    vector_transport_to_project,
    vector_transport_to!,
    vector_transport_to_diff!,
    vector_transport_to_project!, # some overwrite layer 2
    _vector_transport_direction,
    _vector_transport_direction!,
    _vector_transport_to,
    _vector_transport_to!,
    vee,
    vee!,
    zero_vector,
    zero_vector!,
    CotangentSpace,
    TangentSpace
import Base:
    copyto!,
    convert,
    foreach,
    identity,
    in,
    inv,
    isempty,
    length,
    ndims,
    showerror,
    size,
    transpose

using Base.Iterators: repeated
using Distributions
using Einsum: @einsum
using HybridArrays
using Kronecker
using Graphs
using LinearAlgebra
using ManifoldsBase:
    ℝ,
    ℂ,
    ℍ,
    AbstractBasis,
    AbstractDecoratorManifold,
    AbstractInverseRetractionMethod,
    AbstractManifold,
    AbstractManifoldPoint,
    AbstractNumbers,
    AbstractOrthogonalBasis,
    AbstractOrthonormalBasis,
    AbstractPowerManifold,
    AbstractPowerRepresentation,
    AbstractRetractionMethod,
    AbstractTrait,
    AbstractVectorTransportMethod,
    AbstractLinearVectorTransportMethod,
    ApproximateInverseRetraction,
    ApproximateRetraction,
    CachedBasis,
    CayleyRetraction,
    CayleyInverseRetraction,
    ComplexNumbers,
    ComponentManifoldError,
    CompositeManifoldError,
    CotangentSpaceType,
    CoTFVector,
    DefaultBasis,
    DefaultOrthogonalBasis,
    DefaultOrthonormalBasis,
    DefaultOrDiagonalizingBasis,
    DiagonalizingBasisData,
    DiagonalizingOrthonormalBasis,
    DifferentiatedRetractionVectorTransport,
    EmbeddedManifold,
    EmptyTrait,
    ExponentialRetraction,
    FVector,
    IsIsometricEmbeddedManifold,
    IsEmbeddedManifold,
    IsEmbeddedSubmanifold,
    IsExplicitDecorator,
    LogarithmicInverseRetraction,
    ManifoldsBase,
    NestedPowerRepresentation,
    NestedReplacingPowerRepresentation,
    TraitList,
    NLSolveInverseRetraction,
    ODEExponentialRetraction,
    OutOfInjectivityRadiusError,
    PadeRetraction,
    PadeInverseRetraction,
    ParallelTransport,
    PolarInverseRetraction,
    PolarRetraction,
    PoleLadderTransport,
    PowerManifold,
    PowerManifoldNested,
    PowerManifoldNestedReplacing,
    ProjectedOrthonormalBasis,
    ProjectionInverseRetraction,
    ProjectionRetraction,
    ProjectionTransport,
    QuaternionNumbers,
    QRInverseRetraction,
    QRRetraction,
    RealNumbers,
    ScaledVectorTransport,
    SchildsLadderTransport,
    SoftmaxRetraction,
    SoftmaxInverseRetraction,
    TangentSpaceType,
    TCoTSpaceType,
    TFVector,
    TVector,
    ValidationManifold,
    ValidationMPoint,
    ValidationTVector,
    VectorSpaceType,
    VeeOrthogonalBasis,
    @invoke_maker,
    _euclidean_basis_vector,
    combine_allocation_promotion_functions,
    default_inverse_retraction_method,
    geodesic,
    merge_traits,
    next_trait,
    number_system,
    real_dimension,
    rep_size_to_colons,
    shortest_geodesic,
    size_to_tuple,
    trait
using Markdown: @doc_str
using Random
using RecipesBase
using RecipesBase: @recipe, @series
using Colors: RGBA
using Requires
using SimpleWeightedGraphs: AbstractSimpleWeightedGraph, get_weight
using StaticArrays
using Statistics
using StatsBase
using StatsBase: AbstractWeights
using RecursiveArrayTools: ArrayPartition

include("utils.jl")

include("product_representations.jl")
include("differentiation/differentiation.jl")
include("differentiation/riemannian_diff.jl")
include("differentiation/embedded_diff.jl")

# Main Meta Manifolds
include("manifolds/ConnectionManifold.jl")
include("manifolds/MetricManifold.jl")
include("manifolds/VectorBundle.jl")
include("groups/group.jl")

# Features I: Which are extended on Meta Manifolds
include("distributions.jl")
include("projected_distribution.jl")
include("statistics.jl")

# Meta Manifolds II: Products
include("manifolds/ProductManifold.jl")

METAMANIFOLDS = [
    AbstractManifold,
    AbstractDecoratorManifold,
    AbstractPowerManifold,
    PowerManifoldNested,
    PowerManifoldNestedReplacing,
    ProductManifold,
    TangentSpaceAtPoint,
    ValidationManifold,
    VectorBundle,
]

# Features II: That require metas
include("atlases.jl")
include("cotangent_space.jl")

# Meta Manifolds II: Power Manifolds
include("manifolds/PowerManifold.jl")
include("manifolds/GraphManifold.jl")

#
# Manifolds
#
include("manifolds/Euclidean.jl")
include("manifolds/Lorentz.jl")

include("manifolds/CenteredMatrices.jl")
include("manifolds/CholeskySpace.jl")
include("manifolds/Circle.jl")
include("manifolds/Elliptope.jl")
include("manifolds/FixedRankMatrices.jl")
include("manifolds/GeneralizedGrassmann.jl")
include("manifolds/GeneralizedStiefel.jl")
include("manifolds/Grassmann.jl")
include("manifolds/Hyperbolic.jl")
include("manifolds/MultinomialDoublyStochastic.jl")
include("manifolds/MultinomialSymmetric.jl")
include("manifolds/ProbabilitySimplex.jl")
include("manifolds/PositiveNumbers.jl")
include("manifolds/ProjectiveSpace.jl")
include("manifolds/Rotations.jl")
include("manifolds/SkewHermitian.jl")
include("manifolds/Spectrahedron.jl")
include("manifolds/Stiefel.jl")
include("manifolds/StiefelEuclideanMetric.jl")
include("manifolds/StiefelCanonicalMetric.jl")
include("manifolds/Sphere.jl")
include("manifolds/SphereSymmetricMatrices.jl")
include("manifolds/Symmetric.jl")
include("manifolds/SymmetricPositiveDefinite.jl")
include("manifolds/SymmetricPositiveDefiniteLinearAffine.jl")
include("manifolds/SymmetricPositiveDefiniteLogCholesky.jl")
include("manifolds/SymmetricPositiveDefiniteLogEuclidean.jl")
include("manifolds/SymmetricPositiveSemidefiniteFixedRank.jl")
include("manifolds/Tucker.jl")
include("manifolds/Symplectic.jl")
include("manifolds/SymplecticStiefel.jl")

# Product or power based manifolds
include("manifolds/Torus.jl")
include("manifolds/Multinomial.jl")
include("manifolds/Oblique.jl")
include("manifolds/EssentialManifold.jl")

#
# Group Manifolds
include("groups/GroupManifold.jl")

# a) generics
include("groups/addition_operation.jl")
include("groups/multiplication_operation.jl")
include("groups/connections.jl")
include("groups/metric.jl")
include("groups/group_action.jl")
include("groups/group_operation_action.jl")
include("groups/validation_group.jl")
include("groups/product_group.jl")
include("groups/semidirect_product_group.jl")

# Special Group Manifolds
include("groups/general_linear.jl")
include("groups/special_linear.jl")
include("groups/translation_group.jl")
include("groups/special_orthogonal.jl")
include("groups/circle_group.jl")
include("groups/heisenberg.jl")

include("groups/translation_action.jl")
include("groups/rotation_action.jl")

include("groups/special_euclidean.jl")

@doc raw"""
    Base.in(p, M::AbstractManifold; kwargs...)
    p ∈ M

Check, whether a point `p` is a valid point (i.e. in) a [`AbstractManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.AbstractManifold)  `M`.
This method employs [`is_point`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions.html#ManifoldsBase.is_point) deactivating the error throwing option.
"""
Base.in(p, M::AbstractManifold; kwargs...) = is_point(M, p, false; kwargs...)

@doc raw"""
    Base.in(p, TpM::TangentSpaceAtPoint; kwargs...)
    X ∈ TangentSpaceAtPoint(M,p)

Check whether `X` is a tangent vector from (in) the tangent space $T_p\mathcal M$, i.e.
the [`TangentSpaceAtPoint`](@ref) at `p` on the [`AbstractManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.AbstractManifold)  `M`.
This method uses [`is_vector`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions.html#ManifoldsBase.is_vector) deactivating the error throw option.
"""
function Base.in(X, TpM::TangentSpaceAtPoint; kwargs...)
    return is_vector(base_manifold(TpM), TpM.point, X, false; kwargs...)
end

function __init__()
    @require FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" begin
        using .FiniteDifferences
        include("differentiation/finite_differences.jl")
    end

    @require OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" begin
        using .OrdinaryDiffEq: ODEProblem, AutoVern9, Rodas5, solve
        include("differentiation/ode.jl")
    end

    @require NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" begin
        using .NLsolve: NLsolve
        include("nlsolve.jl")
    end

    @require Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" begin
        using .Test: Test
        include("tests/tests_general.jl")
        export test_manifold
        include("tests/tests_group.jl")
        export test_group, test_action
    end

    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
        using RecipesBase: @recipe, @series
        using Colors: RGBA
        include("recipes.jl")
    end

    @require RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" begin
        @require Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" begin
            using .RecipesBase: @recipe, @series
            using Colors: RGBA
            include("recipes.jl")
        end
    end

    return nothing
end

#
export CoTVector, AbstractManifold, AbstractManifoldPoint, TVector
export AbstractSphere, AbstractProjectiveSpace
export Euclidean,
    ArrayProjectiveSpace,
    ArraySphere,
    CenteredMatrices,
    CholeskySpace,
    Circle,
    Elliptope,
    EssentialManifold,
    FixedRankMatrices,
    GeneralizedGrassmann,
    GeneralizedStiefel,
    Grassmann,
    HeisenbergGroup,
    Hyperbolic,
    Lorentz,
    MultinomialDoubleStochastic,
    MultinomialMatrices,
    MultinomialSymmetric,
    Oblique,
    PositiveArrays,
    PositiveMatrices,
    PositiveNumbers,
    PositiveVectors,
    ProbabilitySimplex,
    ProjectiveSpace,
    Rotations,
    SkewHermitianMatrices,
    SkewSymmetricMatrices,
    Spectrahedron,
    Sphere,
    SphereSymmetricMatrices,
    Stiefel,
    SymmetricMatrices,
    SymmetricPositiveDefinite,
    SymmetricPositiveSemidefiniteFixedRank,
    Symplectic,
    SymplecticStiefel,
    SymplecticMatrix,
    Torus,
    Tucker
export HyperboloidPoint, PoincareBallPoint, PoincareHalfSpacePoint, SVDMPoint, TuckerPoint
export HyperboloidTVector,
    PoincareBallTVector, PoincareHalfSpaceTVector, UMVTVector, TuckerTVector
export AbstractNumbers, ℝ, ℂ, ℍ

# decorator manifolds
export AbstractDecoratorManifold
export IsIsometricEmbeddedManifold, IsEmbeddedManifold, IsEmbeddedSubmanifold
export IsDefaultMetric, IsDefaultConnection, IsMetricManifold, IsConnectionManifold
export ValidationManifold, ValidationMPoint, ValidationTVector, ValidationCoTVector
export CotangentBundle,
    CotangentSpaceAtPoint, CotangentBundleFibers, CotangentSpace, FVector
export AbstractPowerManifold,
    AbstractPowerRepresentation,
    ArrayPowerRepresentation,
    NestedPowerRepresentation,
    NestedReplacingPowerRepresentation,
    PowerManifold
export ProductManifold
export GraphManifold, GraphManifoldType, VertexManifold, EdgeManifold
export ProjectedPointDistribution, ProductRepr, TangentBundle, TangentBundleFibers
export TangentSpace, TangentSpaceAtPoint, VectorSpaceAtPoint, VectorSpaceType, VectorBundle
export VectorBundleFibers
export AbstractVectorTransportMethod,
    DifferentiatedRetractionVectorTransport, ParallelTransport, ProjectedPointDistribution
export PoleLadderTransport, SchildsLadderTransport
export ProductVectorTransport
export AbstractAffineConnection,
    AbstractConnectionManifold, ConnectionManifold, LeviCivitaConnection
export AbstractCartanSchoutenConnection,
    CartanSchoutenMinus, CartanSchoutenPlus, CartanSchoutenZero
export AbstractMetric,
    RiemannianMetric,
    LorentzMetric,
    EmbeddedManifold,
    EuclideanMetric,
    LinearAffineMetric,
    LogCholeskyMetric,
    LogEuclideanMetric,
    MinkowskiMetric,
    PowerMetric,
    ProductMetric,
    RealSymplecticMetric,
    ExtendedSymplecticMetric,
    CanonicalMetric,
    MetricManifold
export AbstractAtlas, RetractionAtlas
export AbstractVectorTransportMethod, ParallelTransport, ProjectionTransport
export AbstractRetractionMethod,
    CayleyRetraction,
    ExponentialRetraction,
    QRRetraction,
    PolarRetraction,
    ProjectionRetraction,
    SoftmaxRetraction,
    ODEExponentialRetraction,
    PadeRetraction,
    ProductRetraction,
    PowerRetraction
export AbstractInverseRetractionMethod,
    ApproximateInverseRetraction,
    ApproximateLogarithmicMap,
    CayleyInverseRetraction,
    LogarithmicInverseRetraction,
    QRInverseRetraction,
    PolarInverseRetraction,
    ProjectionInverseRetraction,
    SoftmaxInverseRetraction
export AbstractEstimationMethod,
    GradientDescentEstimation,
    CyclicProximalPointEstimation,
    GeodesicInterpolation,
    GeodesicInterpolationWithinRadius,
    ExtrinsicEstimation
export CachedBasis,
    DefaultBasis,
    DefaultOrthogonalBasis,
    DefaultOrthonormalBasis,
    DiagonalizingOrthonormalBasis,
    InducedBasis,
    ProjectedOrthonormalBasis
export ComponentManifoldError, CompositeManifoldError
export ×,
    allocate,
    allocate_result,
    base_manifold,
    bundle_projection,
    change_metric,
    change_metric!,
    change_representer,
    change_representer!,
    check_point,
    check_vector,
    christoffel_symbols_first,
    christoffel_symbols_second,
    christoffel_symbols_second_jacobian,
    convert,
    complex_dot,
    decorated_manifold,
    det_local_metric,
    distance,
    dual_basis,
    einstein_tensor,
    embed,
    embed!,
    exp,
    exp!,
    flat,
    flat!,
    gaussian_curvature,
    geodesic,
    get_default_atlas,
    get_component,
    get_embedding,
    grad_euclidean_to_manifold,
    grad_euclidean_to_manifold!,
    hat,
    hat!,
    identity_element,
    identity_element!,
    induced_basis,
    incident_log,
    injectivity_radius,
    inner,
    inverse_local_metric,
    inverse_retract,
    inverse_retract!,
    isapprox,
    is_default_connection,
    is_default_metric,
    is_group_manifold,
    is_identity,
    is_point,
    is_vector,
    kurtosis,
    local_metric,
    local_metric_jacobian,
    log,
    log!,
    log_local_metric_density,
    manifold_dimension,
    metric,
    mean,
    mean!,
    mean_and_var,
    mean_and_std,
    median,
    median!,
    mid_point,
    mid_point!,
    minkowski_metric,
    moment,
    norm,
    normal_tvector_distribution,
    number_eltype,
    one,
    power_dimensions,
    project,
    project!,
    projected_distribution,
    real_dimension,
    ricci_curvature,
    ricci_tensor,
    representation_size,
    retract,
    retract!,
    riemann_tensor,
    set_component!,
    sharp,
    sharp!,
    shortest_geodesic,
    skewness,
    std,
    sym_rem,
    symplectic_inverse_times,
    symplectic_inverse_times!,
    submanifold,
    submanifold_component,
    submanifold_components,
    uniform_distribution,
    var,
    vector_space_dimension,
    vector_transport_along,
    vector_transport_along!,
    vector_transport_direction,
    vector_transport_direction!,
    vector_transport_to,
    vector_transport_to!,
    vee,
    vee!,
    zero_vector,
    zero_vector!
# Lie group types & functions
export AbstractGroupAction,
    AbstractGroupOperation,
    ActionDirection,
    AdditionOperation,
    CircleGroup,
    GeneralLinear,
    GroupManifold,
    GroupOperationAction,
    Identity,
    InvariantMetric,
    LeftAction,
    LeftInvariantMetric,
    MultiplicationOperation,
    ProductGroup,
    ProductOperation,
    RealCircleGroup,
    RightAction,
    RightInvariantMetric,
    RotationAction,
    SemidirectProductGroup,
    SpecialEuclidean,
    SpecialLinear,
    SpecialOrthogonal,
    TranslationGroup,
    TranslationAction
export AbstractInvarianceTrait
export IsMetricManifold, IsConnectionManifold
export IsGroupManifold,
    HasLeftInvariantMetric, HasRightInvariantMetric, HasBiinvariantMetric
export adjoint_action,
    adjoint_action!,
    affine_matrix,
    apply,
    apply!,
    apply_diff,
    apply_diff!,
    base_group,
    center_of_orbit,
    has_approx_invariant_metric,
    compose,
    compose!,
    direction,
    exp_lie,
    exp_lie!,
    group_manifold,
    geodesic,
    get_coordinates_lie,
    get_coordinates_lie!,
    get_coordinates_orthogonal,
    get_coordinates_orthonormal,
    get_coordinates_orthogonal!,
    get_coordinates_orthonormal!,
    get_vector_diagonalizing!,
    get_vector_lie,
    get_vector_lie!,
    get_vector_orthogonal,
    get_vector_orthonormal,
    get_coordinates_vee!,
    has_biinvariant_metric,
    has_invariant_metric,
    identity_element,
    identity_element!,
    inv,
    inv!,
    inverse_apply,
    inverse_apply!,
    inverse_apply_diff,
    inverse_apply_diff!,
    inverse_translate,
    inverse_translate!,
    inverse_translate_diff,
    inverse_translate_diff!,
    lie_bracket,
    lie_bracket!,
    log_lie,
    log_lie!,
    optimal_alignment,
    optimal_alignment!,
    screw_matrix,
    switch_direction,
    translate,
    translate!,
    translate_diff,
    translate_diff!
# Orthonormal bases
export AbstractBasis,
    AbstractOrthonormalBasis,
    DefaultOrthonormalBasis,
    DiagonalizingOrthonormalBasis,
    ProjectedOrthonormalBasis,
    CachedBasis,
    DiagonalizingBasisData,
    ProductBasisData,
    PowerBasisData
export OutOfInjectivityRadiusError
export get_basis,
    get_coordinates, get_coordinates!, get_vector, get_vector!, get_vectors, number_system

# atlases and charts
export get_point, get_point!, get_parameters, get_parameters!

end # module

max_stars_repo_name: JuliaNLSolvers/ManifoldMuseum

max_issues_repo_name: JuliaNLSolvers/ManifoldMuseum

max_forks_repo_name: JuliaNLSolvers/ManifoldMuseum


text:
using PyPlot
using Images

## Load an image from a file, rescale its dynamic to [0,1], turn it into a grayscale image and resize it to size n x n.

function load_image(name, n = -1, flatten = 1, resc = 1, grayscale = 1)

  if ndims(PyPlot.imread(name)) >=3
    f = PyPlot.imread(name)[:, :, 1:3]
  end

  f = PyPlot.imread(name)

  if grayscale == 1
    if (flatten == 1) & (ndims(f) > 2)
      f = sum(f, dims = 3)   # We must be careful, this returns an array of dimension 512*512*1 while
      f = f[:, :]     # in the python file np.sum(f, axis=2) returns an array of dimension 512*512

    end
  end
  if resc == 1
    f = rescale(f)  # rescale must be defined in the file general.jl of the package NtToolBox
                    # and we must find a way to import it to this file signal.jl
  end
  if n > 0
    if ndims(f) == 3
        f = Images.imresize(f,(n,n,3))
    elseif ndims(f) == 2
      f = Images.imresize(f, (n, n))
    end
  end
  return f
end

##

## Use nearest neighbor interpolation for the display.

function imageplot(f, str = "", sbpt = [])    # Function to check.

  if sbpt != []
    subplot(sbpt[1], sbpt[2], sbpt[3])
  end
  imshow(f, interpolation = "nearest")
  set_cmap("gray")
  axis("off")
  if str != ""
    title(str)
  end
end

function snr(x, y)
  """
  snr - signal to noise ratio

     v = snr(x,y);

   v = 20*log10( norm(x(:)) / norm(x(:)-y(:)) )

     x is the original clean signal (reference).
     y is the denoised signal.

  Copyright (c) 2014 Gabriel Peyre
  """
  return 20 * log10(norm(x[:]) / norm(x[:] - y[:]))
end

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

function resize_image(f, N)
	P = size(f,1);
	# add 1 pixel boundary
	g = f;
	g = cat(2, g, g[:,1,:]);
	g = cat(1, g, g[1,:,:]);
	# interpolate
	t = range(1, stop=P, length=N);
	ti = floor(t); tj = ceil(t);
	fi = t-floor(t); fj = 1-fi;
	h = zeros(N,N,size(f,3));
	for s=1:size(f,3)
	    h[:,:,s] = g[ti,ti,s] .* (fj*fj') + g[tj,tj,s] .* (fi*fi') + g[ti,tj,s] .* (fj*fi') + g[tj,ti,s] .* (fi*fj');
	end
	return h;
end

###################################################
function trim_dim(x)
	# remove empty trailing dimensions
	if size(x,3)==1
		x = x[:,:,1];
		if size(x,2)==1
			x = x[:,1];
		end
	end
	return x;
end

###################################################
function subsampling(x,d=1,p=2)

	# downsampling - subsampling along dimension d
	#
	#   y = subsampling(x,d,p);
	#
	#   default is p==2, d==1
	#
	#   Copyright (c) 2009 Gabriel Peyre

	if d==1
	        y = x[1:p:end,:,:];
	elseif d==2
	        y = x[:,1:p:end,:];
	elseif d==3
	        y = x[:,:,1:p:end];
	else
	        error("Not implemented");
	end
	y = trim_dim(y);

	return y;

end # subsampling

###################################################
function upsampling(x,d=1,p=2)

	# upsampling - add p zeros between samples along dimension d
	#
	#   y = upsampling(x,d,p);
	#
	#   default is p==2, d==1
	#
	#   Copyright (c) 2009 Gabriel Peyre

	if d==1
	        y = zeros(p*size(x,1),size(x,2),size(x,3));
	        y[1:p:end,:,:] = x;
	elseif d==2
	        y = zeros(size(x,1),p*size(x,2),size(x,3));
	        y[:,1:p:end,:] = x;
	elseif d==3
	        y = zeros(size(x,1),size(x,2),p*size(x,3));
	        y[:,:,1:p:end] = x;
	else
	        error("Not implemented");
	end
	y = trim_dim(y);

	return y;
end # upsampling

###################################################
function cconvol(x,h,d=1)

	# cconvol - spatial domain circular convolution
	#
	#   y = cconv(x,h,d);
	#
	#   Circular convolution along dimension |d|.
	#
	#   If p=length(h) is odd then h((p+1)/2) corresponds to the zero position.
	#   If p is even, h(p/2) corresponds to the zero position.
	#
	#   Copyright (c) 2009 Gabriel Peyre

	if d==2
	    return ( cconvol( x',h, 1 )' );
	elseif d==3
		# permute do not exists in python
	    return permute( cconvol( permute(x,[3, 2, 1]),h), [3 2 1]);
	end

	p = length(h);
	if mod(p,2)==0
	    pc = p/2;
	else
	    pc = (p+1)/2;
	end

	y = zeros(size(x));
	for i=1:length(h)
	    y = y + h[i]*circshift(x,i-pc);
	end
	return y;

end # cconvol

#####################################################
# function plot_wavelet(fW, Jmin=0)
#     #    plot_wavelet - plot wavelets coefficients.
#     #    U = plot_wavelet(fW, Jmin):
#     #    Copyright (c) 2014 Gabriel Peyre
#
#     function rescaleWav(A)
#         v = maximum( abs(A[:]) );
#         B = copy(A)
#         if v > 0
#             B = .5 + .5 * A / v;
#         end
#         return B;
#     end
#     function rescale(x,a=0,b=1)
# 		m = minimum(x[:]);
# 		M = maximum(x[:]);
# 		if M-m<1e-10
# 		    y = x;
# 		else
# 		    y = (b-a) * (x-m)/(M-m) + a;
# 		end
# 		return y;
# 	  end ## end of rescale
#     ##
#     n = size(fW,1);
#     Jmax = Int(log2(n) - 1);
#     U = copy(fW);
#     for j = Jmax:-1:Jmin
#     	L = Array{Int64,1}(1:2^j); H = Array{Int64,1}((2^j)+1:2^(j+1));
#         U[L,H] = rescaleWav( U[L,H] );
#         U[H,L] = rescaleWav( U[H,L] );
#         U[H,H] = rescaleWav( U[H,H] );
#     end
#     # coarse scale
#     L = Array{Int64,1}(1:2^Jmin);
#     U[L,L] = rescale(U[L,L]);
#     # plot underlying image
#     imageplot(U);
#     # display crosses
#     for j=Jmax:-1:Jmin
#         plot([0, 2^(j+1)], [2^j, 2^j], "r");
#         plot([2^j, 2^j], [0, 2^(j + 1)], "r");
#     end
#     # display box
#     plot([0, n], [0, 0], "r");
#     plot([0, n], [n, n], "r");
#     plot([0, 0], [0, n], "r");
#     plot([n, n], [0, n], "r");
#     return U;
#
# end # plot_wavelet

#########
function plot_wavelet(fW, Jmin = 0)
    """
        plot_wavelet - plot wavelets coefficients.

        U = plot_wavelet(fW, Jmin):

        Copyright (c) 2014 Gabriel Peyre
    """
    function rescaleWav(A)
        v = maximum(abs(A[:]))
        B = copy(A)
        if v > 0
            B = .5 + .5 .* A ./ v
        end
        return B
    end
    ##
    n = size(fW)[2]
    Jmax = Int(log2(n) - 1)
    U = copy(fW)
    for j in Jmax:-1:Jmin
        U[1:2^j,    2^j+1:2 ^
            (j + 1)] = rescaleWav(U[1:2^j, 2^j+1:2^(j + 1)])
        U[2^j+1:2^(j+1), 1:2 ^
          j] = rescaleWav(U[2^j+1:2^(j+1), 1:2^j])
        U[2^j+1:2^(j+1), 2^j+1:2^(j+1)] = (
            rescaleWav(U[2^j+1:2^(j+1), 2^j+1:2^(j+1)]))
    end

    # coarse scale
    U[1:2^Jmin, 1:2^Jmin] = rescale(U[1:2^Jmin, 1:2^Jmin])

    # plot underlying image
    imageplot(U)
    # display crosses
    for j in Jmax:-1:Jmin
        plot([0, 2^(j+1)], [2^j, 2^j], "r")
        plot([2^j, 2^j], [0, 2^(j+1)], "r")
    end
    # display box
    plot([0, n], [0, 0], "r")
    plot([0, n], [n, n], "r")
    plot([0, 0], [0, n], "r")
    plot([n, n], [0, n], "r")
    return U
end

########

#########################################
function subselectdim(f,sel,d)
	g = [];
	if d==1
	    g = f[sel,:,:];
	elseif d==2
	    g = f[:,sel,:];
	elseif d==3
	    g = f[:,:,sel];
	end
	return trim_dim(g);
end

#########################################
function subassign(f,sel,g)
	d = ndims(f);
	if d==1
	    f[sel] = g;
	elseif d==2
	    f[sel,sel] = g;
	elseif d==3
	    f[sel,sel,sel] = g;
	end
	return f;
end # subassign

#########################################
function subselect(f,sel)
	d = ndims(f)
    if d==1
        return f[sel];
    elseif d==2
        return f[sel,sel];
    elseif d==3
        return f[sel,sel,sel];
	end
	return [];
end

##############################################################
function perform_wavortho_transf(f0,Jmin,dir,h)

	# perform_wavortho_transf - compute orthogonal wavelet transform
	#
	#   fw = perform_wavortho_transf[f,Jmin,dir,options);
	#
	#   You can give the filter in options.h.
	#
	#   Works in arbitrary dimension.
	#
	#   Copyright (c) 2009 Gabriel Peyre

	n = size(f0,1);
	Jmax = Int(log2(n)-1);

	g = [0; h[length(h):-1:2]] .* (-1).^(1:length(h)); # Je pense qu'il faut changer [0, h[length(h):-1:2]] en [0; h[length(h):-1:2]]

	f = copy(f0)
	if dir==1
	    ### FORWARD ###
	    for j=Jmax:-1:Jmin
	        sel = 1:2^(j+1);
	        a = subselect(f,sel);
	        for d=1:ndims(f)
	            a = cat(d, subsampling(cconvol(a,h,d),d), subsampling(cconvol(a,g,d),d) );
	        end
	        f = subassign(f,sel,a);
	    end
	else
	    ### FORWARD ###
	    for j=Jmin:Jmax
	        sel = 1:2^(j+1);
	        a = subselect(f,sel);
	        for d=1:ndims(f)
	            w = subselectdim(a,2^j+1:2^(j+1),d);
	            a = subselectdim(a,1:2^j,d);
	            a = cconvol(upsampling(a,d),reverse(h),d) + cconvol(upsampling(w,d),reverse(g),d);
	        end
	        f = subassign(f,sel,a);
	    end
	end
	return f;
end # perform_wavortho_transf

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

function div(g)
  """
      Compute a finite difference approximation of the gradient of a 2D vector field, assuming periodic BC.
  """

  S = size(g)
  s0 = [[S[1]]; collect(1:S[1] - 1)]
  s1 = [[S[2]]; collect(1:S[2] - 1)]
  f = (g[:, :, 1] - g[s0, :, 1]) + (g[:, : , 2] - g[:, s1, 2])
  return f
end

function grad(f)
    """
        Compute a finite difference approximation of the gradient of a 2D image, assuming periodic BC.
    """
    S = size(f)
    s0 = [collect(2:S[1]); [1]]
    s1 = [collect(2:S[2]); [1]]
    g = cat(3, f[s0, :] - f, f[:, s1] - f)
    return g
end

function bilinear_interpolate(im, x, y)

    # x0 = Array{Int64,1}(floor(x))
    x0 = [Int(floor(x))]
    x1 = x0 + 1
    # y0 = Array{Int64,1}(floor(y))
    y0 = [Int(floor(y))]
    y1 = y0 + 1

    x0 = clamp(x0, 1, size(im)[2])
    x1 = clamp(x1, 1, size(im)[2])
    y0 = clamp(y0, 1, size(im)[1])
    y1 = clamp(y1, 1, size(im)[1])

    Ia = im[ y0, x0 ]
    Ib = im[ y1, x0 ]
    Ic = im[ y0, x1 ]
    Id = im[ y1, x1 ]

    wa = (x1 - x) .* (y1 - y)
    wb = (x1 - x) .* (y - y0)
    wc = (x - x0) .* (y1 - y)
    wd = (x - x0) .* (y - y0)

    return wa.*Ia + wb.*Ib + wc.*Ic + wd.*Id
end


text:
using Monkey
using Test

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


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

name = "SCIP"
version = v"0.1.2"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://scip.zib.de/download/release/scipoptsuite-7.0.2.tgz", "f81b5a2c1c0eb949cf06bd50f42826e55284fa1269a6f28a92ac1a06d9c93a03")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd scipoptsuite*
mkdir build
cd build/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DPAPILO=0 -DZIMPL=OFF -DGCG=0 -DIPOPT=0 -DREADLINE=OFF -DBOOST=off ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("i686", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("armv7l", "linux"; libc="glibc"),
    Platform("powerpc64le", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="musl"),
    Platform("x86_64", "linux"; libc="musl"),
    Platform("aarch64", "linux"; libc="musl"),
    Platform("armv7l", "linux"; libc="musl"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "freebsd"),
    Platform("i686", "windows"),
    Platform("x86_64", "windows"),
]

platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libscip", :libscip),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="bliss_jll", uuid="508c9074-7a14-5c94-9582-3d4bc1871065")),
    Dependency(PackageSpec(name="GMP_jll", uuid="781609d7-10c4-51f6-84f2-b8444358ff6d"), v"6.1.2"),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
]

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

max_stars_repo_name: Affie/Yggdrasil

max_issues_repo_name: Affie/Yggdrasil

max_forks_repo_name: Affie/Yggdrasil