text:
# functions required to start a solve

import PDESolver.solvePDE

"""
  This function creates and initializes the mesh, sbp, eqn, and opts objects

  Inputs:
    opts: options dictionary

  Outputs:
    mesh: an AbstractMesh.  The concrete type is determined by the options
          dictionary
    sbp: an AbstractOperator.  The concrete type is determined by the options
         dictionary
    eqn: an NSData object
    opts: the options dictionary
    pmesh: mesh used for preconditioning, can be same object as mesh
"""
function createObjects(opts::Dict)

  # make sure Euler options were set
  EulerEquationMod.checkOptions(opts)

  Tdim = opts["dimensions"]
  dofpernode = Tdim + 2

  sbp, mesh, pmesh, Tsol, Tres, Tmsh, mesh_time = createMeshAndOperator(opts, dofpernode)

  # create euler equation
  eqn = NSData_{Tsol, Tres, Tdim, Tmsh}(mesh, sbp, opts)

  # initialize physics module and populate any fields in mesh and eqn that
  # depend on the physics module
  init(mesh, sbp, eqn, opts, pmesh)

  return mesh, sbp, eqn, opts, pmesh
end

"""
  Given fully initialized mesh, sbp, eqn, opts, this function solves
  the Euler equations.  The 4 object should be obtained from createObjects().

  Specifically, it applies an initial condition and invokes a nonlinear
  solver according to the options dictionary.

  Inputs:
    mesh: an AbstractMesh
    sbp: an AbstractOperator
    eqn: an NSData
    opts: the options dictionary.  This must be the options dictionary returned
          by createObjects().  Changing values in the options dictionary after
          calling createObjects() results in undefined behavior.
    pmesh: mesh used for preconditioning, can be same object as mesh.
           default value of mesh

"""
function solvePDE(mesh::AbstractMesh, sbp::AbstractOperator, eqn::NSData, opts::Dict, pmesh::AbstractMesh=mesh)
  #delta_t = opts["delta_t"]   # delta_t: timestep for RK

  myrank = mesh.myrank

  # zero some arrays out in case the user forgot

  fill!(eqn.res, 0.0)
  fill!(eqn.res_vec, 0.0)
  res_vec = eqn.res_vec
  q_vec = eqn.q_vec       # solution at current timestep

  #TODO: move as many of these if statements into functions as possible

  # calculate residual of some other function for res_reltol0
  # TODO: add a boolean options here?
  Relfunc_name = opts["Relfunc_name"]
  if haskey(ICDict, Relfunc_name)
    @mpi_master println("\ncalculating residual for relative residual tolerance")
    Relfunc = ICDict[Relfunc_name]
    @mpi_master println("Relfunc = ", Relfunc)
    Relfunc(mesh, sbp, eqn, opts, q_vec)
    tmp = physicsRhs(mesh, sbp, eqn, opts, eqn.res_vec, (evalResidual,))
    res_real = real(eqn.res_vec)
    opts["res_reltol0"] = tmp
    @mpi_master println("res_reltol0 = ", tmp)

    saveSolutionToMesh(mesh, res_real)
    writeVisFiles(mesh, "solution_relfunc")
  end

  # populate u0 with initial condition
  @mpi_master println("\nEvaluating initial condition")
  ICfunc_name = opts["IC_name"]
  ICfunc = ICDict[ICfunc_name]
  @mpi_master println("ICfunc = ", ICfunc)
  ICfunc(mesh, sbp, eqn, opts, q_vec)

  if opts["calc_error"]
    @mpi_master println("\ncalculating error of file ",
                       opts["calc_error_infname"],
                      " compared to initial condition")

    # read in this processors portion of the solution
    vals = copy(eqn.q_vec)
    readSolutionFiles(mesh, sbp, eqn, opts, opts["calc_error_infname"], vals)

    err_vec = abs.(vals - eqn.q_vec)
    err = calcNorm(eqn, err_vec)

    # calculate avg mesh size
    h_avg = calcMeshH(mesh, sbp, eqn, opts)

    @mpi_master begin
      outname = opts["calc_error_outfname"]
      @mpi_master println("printint err = ", err, " to file ", outname)
      f = open(outname, "w")
      println(f, err, " ", h_avg)
      close(f)
    end

    # write visualization
    saveSolutionToMesh(mesh, vec(err_vec))
    writeVisFiles(mesh, "solution_error")
  end

  if opts["calc_trunc_error"]  # calculate truncation error
    @mpi_master println("\nCalculating residual for truncation error")
    tmp = physicsRhs(mesh, sbp, eqn, opts, eqn.res_vec, (evalResidual,))

    @mpi_master begin
      f = open("error_trunc.dat", "w")
      println(f, tmp)
      close(f)
    end
  end

  if opts["perturb_ic"]
    @mpi_master println("\nPerturbing initial condition")
    perturb_mag = opts["perturb_mag"]
    for i=1:mesh.numDof
      q_vec[i] += perturb_mag*rand()
    end
  end

  #writeSolutionFiles(mesh, sbp, eqn, opts, "IC")
  saveSolutionToMesh(mesh, q_vec)

  writeVisFiles(mesh, "solution_ic")
  if opts["calc_dt"]
    wave_speed = EulerEquationMod.calcMaxWaveSpeed(mesh, sbp, eqn, opts)
    @mpi_master println("max wave speed = ", wave_speed)
    @mpi_master println("min element size = ", mesh.min_el_size)
    delta_t = opts["CFL"]*mesh.min_el_size/wave_speed
    @mpi_master println("for a CFL of ", opts["CFL"], " delta_t = ", delta_t)
    opts["delta_t"] = delta_t
  end

  # this must be last, because all previous calculations need to use the
  # original initial condition
  if opts["is_restart"]
    loadRestartState(mesh, sbp, eqn, opts, pmesh)
  end

  call_nlsolver(mesh, sbp, eqn, opts, pmesh)
  postproc(mesh, sbp, eqn, opts)

  EulerEquationMod.cleanup(eqn.euler_eqn, opts)

  MPI.Barrier(mesh.comm)

  return mesh, sbp, eqn, opts
end  # end function

"""
  This function does post processing, if requested by the input options.
  Typical post processing includes calculation of errors, norms of important
  quantities, writing of files. etc.

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
"""
function postproc(mesh, sbp, eqn, opts)

  ##### Do postprocessing ######
  myrank = mesh.myrank
  @mpi_master println("\nDoing postprocessing")

  if opts["do_postproc"] && opts["solve"]
    exfname = opts["exact_soln_func"]
    if haskey(ICDict, exfname)
      exfunc = ICDict[exfname]
      q_exact = zeros(eqn.q_vec)
      exfunc(mesh, sbp, eqn, opts, q_exact)

      myrank = mesh.myrank
      q_diff = eqn.q_vec - q_exact
      saveSolutionToMesh(mesh, abs.(real(q_diff)))
      writeVisFiles(mesh, "solution_error")

      diff_norm = calcNorm(eqn, q_diff)
      @mpi_master println("solution error norm = ", diff_norm)
      # TODO: make this mesh.min_el_size?
      h_avg = calcMeshH(mesh, sbp, eqn, opts)

      # print to file
      @mpi_master begin
        outname = opts["calc_error_outfname"]
        f = open(outname, "w")
        println(f, diff_norm, " ", h_avg)
        close(f)
      end
    end  # end if haskey(ICname)
  end  # end if do_postproc

  return nothing
end


text:
# compare_classical_and_graph_discrepancies
#
# This script compares the values of the univariate classical and graph Stein
# discrepancies obtained when samples are drawn i.i.d. from their targets.
# We use non-uniform discrepancies with the best Stein factors drawn from the
# literature to enable comparison with a lower-bounding Wasserstein distance.

include("experiment_utils.jl")

using SteinDistributions: SteinGaussian, SteinUniform, gradlogdensity,
    supportlowerbound, supportupperbound, cdf
using SteinDiscrepancy: stein_discrepancy, wasserstein1d

# random number generator seed
# we allow this to be set on the commandline
@parseintcli seed "s" "seed" 7
# set the target distribution
@parsestringcli distname "t" "target" "uniform"
# Specify dimension of sampled vectors
d = 1
# Specify a maximum number of samples
n = 30000
if distname == "uniform"
    # Independent Uniform([0.0,1.0]) with best known Stein constants
    (c1,c2,c3) = (0.5,0.5,1.0)
    target = SteinUniform(d)
elseif distname == "gaussian"
    # Independent standard Gaussians with best known Stein constants
    (c1,c2,c3) = (1.0,4.0,2.0)
    target = SteinGaussian(d)
end
lowerbounds = [supportlowerbound(target, 1)]
upperbounds = [supportupperbound(target, 1)]
# define gradlogp
function gradlogp(x::Array{Float64,1})
    gradlogdensity(target, x)
end
# define targetcdf
function targetcdf(x::Float64)
    cdf(target, x)
end
# setup solvers for stein optimization problem
#graphsolver = classicalsolver = nothing
graphsolver = classicalsolver = "gurobi"
#graphsolver = "clp" # Only solves LPs
# Draw samples from target distribution
X = @setseed rand(target, n)
# Sample sizes at which optimization problem will be solved
ns = vcat(100:100:min(1000,n), 2000:1000:min(10000,n), 20000:10000:n)

# Solve optimization problem at each sample size
for i in ns
    @printf("[Beginning n=%d]\n", i)
    # Create program with data subset
    Xi = X[1:i,:]

    classicalobjective = graphobjective = nothing
    classicalsolvetime = graphsolvetime = nothing
    if classicalsolver != nothing
        @printf("[Beginning classical solver n=%d]\n", i)
        classicalresult = stein_discrepancy(points=Xi,
                                            gradlogdensity=gradlogp,
                                            solver=classicalsolver,
                                            method="classical",
                                            c1=c1,
                                            c2=c2,
                                            c3=c3,
                                            supportlowerbounds=lowerbounds,
                                            supportupperbounds=upperbounds)
        classicalobjective = classicalresult.objectivevalue
        classicalsolvetime = classicalresult.solvetime
    end
    if graphsolver != nothing
        @printf("[Beginning graph solver n=%d]\n", i)
        graphresult = stein_discrepancy(points=Xi,
                                        gradlogdensity=gradlogp,
                                        solver=graphsolver,
                                        c1=c1,
                                        c2=c2,
                                        c3=c3,
                                        supportlowerbounds=lowerbounds,
                                        supportupperbounds=upperbounds)
        graphobjective = graphresult.objectivevalue
        graphsolvetime = graphresult.solvetime
    end

    @printf("[Computing wasserstein distance n=%d]\n", i)
    (wassdist, error) = wasserstein1d(Xi, targetcdf=targetcdf)

    # save data
    instance_data = Dict{Any, Any}(
        "distname" => distname,
        "n" => i,
        "d" => d,
        "seed" => seed,
        "graphobjective" => graphobjective,
        "graphsolvetime" => graphsolvetime,
        "classicalobjective" => classicalobjective,
        "classicalsolvetime" => classicalsolvetime,
        "wassdist" => wassdist
    )
    save_json(
        instance_data;
        dir="compare_classical_and_graph_discrepancies",
        distname=distname,
        n=i,
        d=d,
        seed=seed
    )
end

@printf("COMPLETE!")

ext: jl


text:
using FastBEAST
using StaticArrays
using LinearAlgebra
using Printf

function OneoverRkernel(testpoint::SVector{3,T}, sourcepoint::SVector{3,T}) where T
    if isapprox(testpoint, sourcepoint, rtol=eps()*1e1)
        return 0.0
    else
        return 1.0 / (norm(testpoint - sourcepoint))
    end
end

function assembler(kernel, testpoints, sourcepoints)
    kernelmatrix = zeros(
        promote_type(eltype(testpoints[1]),eltype(sourcepoints[1])),
        length(testpoints),
        length(sourcepoints)
    )

    for i = 1:length(testpoints)
        for j = 1:length(sourcepoints)
            kernelmatrix[i,j] = kernel(testpoints[i], sourcepoints[j])
        end
    end
    return kernelmatrix
end

function assembler(kernel, matrix, testpoints, sourcepoints)
    for i = 1:length(testpoints)
        for j = 1:length(sourcepoints)
            matrix[i,j] = kernel(testpoints[i], sourcepoints[j])
        end
    end
end

##
N =  1000
NT = N

spoints = [@SVector rand(3) for i = 1:N]
tpoints = 0.1*[@SVector rand(3) for i = 1:NT] + [1.0*SVector(3.5, 3.5, 3.5) for i = 1:NT]

@views OneoverRkernelassembler(matrix, tdata, sdata) = assembler(
    OneoverRkernel,
    matrix,
    tpoints[tdata],
    spoints[sdata]
)

stree = create_tree(spoints, BoxTreeOptions(nmin=50))
ttree = create_tree(tpoints, BoxTreeOptions(nmin=50))
kmat = assembler(OneoverRkernel, tpoints, spoints)
hmat = HMatrix(OneoverRkernelassembler, ttree, stree, compressor=:naive, T=Float64)

@printf("Accuracy test: %.2e\n", estimate_reldifference(hmat,kmat))
@printf("Compression rate: %.2f %%\n", compressionrate(hmat)*100)

##
N = 3000
NT = N

spoints = [@SVector rand(3) for i = 1:N]

##
@views OneoverRkernelassembler(matrix, tdata, sdata) = assembler(
    OneoverRkernel,
    matrix,
    spoints[tdata],
    spoints[sdata]
)

stree = create_tree(spoints, BoxTreeOptions(nmin=100))
kmat = assembler(OneoverRkernel, spoints, spoints)
hmat = HMatrix(OneoverRkernelassembler, stree, stree, compressor=:aca, T=Float64)

@printf("Accuracy test: %.2e\n", estimate_reldifference(hmat,kmat))
@printf("Compression rate: %.2f %%\n", compressionrate(hmat)*100)

##
totalspoints = 0
for child in stree.children
    totalspoints += length(child.data)
end

@assert totalspoints == N

##
function hmatrix3d_benchmark(N)
    spoints = [@SVector rand(3) for i = 1:N]

    @views OneoverRkernelassembler(matrix, tdata, sdata) = assembler(
        OneoverRkernel,
        matrix,
        spoints[tdata],
        spoints[sdata]
    )

    stree = create_tree(spoints, nmin=400)
    @time hmat = HMatrix(
        OneoverRkernelassembler,
        stree,
        stree,
        compressor=:aca,
        tol=1e-4,
        T=Float64
    )

    @printf("Memory usage: %.2f GiB\n", nnz(hmat)*8/(1024*1024*1024.0))
    @printf("Compression rate: %.2f %%\n", compressionrate(hmat)*100)
    return nothing #stree, hmat
end

## Example1 a) Random distribution
N = 3000
spoints = [@SVector rand(3) for i = 1:N];

## Example2 a) Not evenly distributed sphere
spoints = [SVector((sin(i)*cos(j),sin(i)*sin(j),cos(i))) for j=0:0.1:2*pi for i = 0:0.1:pi];

## Example2 b) Not evenly distributed sphere like shape
spoints = [SVector((sin(i)*cos(j),sin(i)*sin(j),i^2*cos(i))) for j=0:0.1:2*pi for i = 0:0.1:pi];

##
@views OneoverRkernelassembler(matrix, tdata, sdata) = assembler(OneoverRkernel, matrix, spoints[tdata], spoints[sdata])
stree = create_tree(spoints, BoxTreeOptions(nmin=400))
kmat = assembler(OneoverRkernel, spoints, spoints)
@time hmat = HMatrix(OneoverRkernelassembler, stree, stree, compressor=:aca, T=Float64)

@printf("Accuracy test: %.2e\n", estimate_reldifference(hmat,kmat))
@printf("Compression rate: %.2f %%\n", compressionrate(hmat)*100)

stree = create_tree(spoints, KMeansTreeOptions(iterations=100,nchildren=2,nmin=10))
kmat = assembler(OneoverRkernel, spoints, spoints)
@time hmat = HMatrix(OneoverRkernelassembler, stree, stree, compressor=:aca, T=Float64)

@printf("Accuracy test: %.2e\n", estimate_reldifference(hmat,kmat))
@printf("Compression rate: %.2f %%\n", compressionrate(hmat)*100)

## Plot of first level KMeans clustering
scatter(
    [spoints[stree.children[1].data[i]][1] for i=1:length(stree.children[1].data)],
    [spoints[stree.children[1].data[i]][2] for i=1:length(stree.children[1].data)],
    [spoints[stree.children[1].data[i]][3] for i=1:length(stree.children[1].data)]
)
for j = 2:2
    scatter!(
        [spoints[stree.children[j].data[i]][1] for i=1:length(stree.children[j].data)],
        [spoints[stree.children[j].data[i]][2] for i=1:length(stree.children[j].data)],
        [spoints[stree.children[j].data[i]][3] for i=1:length(stree.children[j].data)]
    )
end
scatter!()

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: sbadrian/FastBEAST

max_forks_repo_name: sbadrian/FastBEAST


text:
struct TargetLoss
    policy::Float32
    reward::Float32
    value::Float32
end

gettotal(x::TargetLoss) = x.policy + x.reward + x.value
getpolicy(x::TargetLoss) = x.policy
getreward(x::TargetLoss) = x.reward
getvalue(x::TargetLoss) = x.value

function Base.:+(a::TargetLoss, b::TargetLoss)
    TargetLoss(
        a.policy + b.policy,
        a.reward + b.reward,
        a.value + b.value,
    )
end


text:
export ReplaySubject, ReplaySubjectFactory

import DataStructures: CircularBuffer, capacity
import Base: show, similar

"""
    ReplaySubject(::Type{D}, size::Int) where D
    ReplaySubject(::Type{D}, size::Int, factory::F) where { D, F <: AbstractSubjectFactory }
    ReplaySubject(::Type{D}, size::Int, subject::S) where { D, S }

A variant of Subject that "replays" or emits old values to new subscribers.
It buffers a set number of values and will emit those values immediately to any new subscribers
in addition to emitting new values to existing subscribers.

See also: [`ReplaySubjectFactory`](@ref), [`Subject`](@ref), [`SubjectFactory`](@ref)
"""
function ReplaySubject end

"""
    ReplaySubjectFactory(size::Int, factory::F) where { F <: AbstractSubjectFactory }
    ReplaySubjectFactory(size::Int; scheduler::H = AsapScheduler()) where { H <: AbstractScheduler }

A variant of SubjectFactory that creates an instance of ReplaySubject.

See also: [`SubjectFactory`](@ref), [`AbstractSubjectFactory`](@ref), [`ReplaySubject`](@ref), [`Subject`](@ref)
"""
function ReplaySubjectFactory end

struct ReplaySubjectInstance{D, S} <: AbstractSubject{D}
    subject :: S
    buffer  :: CircularBuffer{D}
end

Base.show(io::IO, ::ReplaySubjectInstance{D, S}) where { D, S } = print(io, "ReplaySubjectInstance($D, $S)")

Base.similar(subject::ReplaySubjectInstance{D, S}) where { D, S } = ReplaySubject(D, capacity(subject.buffer), similar(subject.subject))

function ReplaySubject(::Type{D}, size::Int) where D
    return ReplaySubject(D, size, SubjectFactory(AsapScheduler()))
end

function ReplaySubject(::Type{D}, size::Int, factory::F) where { D, F <: AbstractSubjectFactory }
    return ReplaySubject(D, size, create_subject(D, factory))
end

function ReplaySubject(::Type{D}, size::Int, subject::S) where { D, S }
    return as_replay_subject(D, as_subject(S), size, subject)
end

as_replay_subject(::Type{D},  ::InvalidSubjectTrait,    size::Int, subject)    where D          = throw(InvalidSubjectTraitUsageError(subject))
as_replay_subject(::Type{D1}, ::ValidSubjectTrait{D2},  size::Int, subject)    where { D1, D2 } = throw(InconsistentSubjectDataTypesError{D1, D2}(subject))
as_replay_subject(::Type{D},  ::ValidSubjectTrait{D},   size::Int, subject::S) where { D, S }   = ReplaySubjectInstance{D, S}(subject, CircularBuffer{D}(size))

##

function on_next!(subject::ReplaySubjectInstance{D}, data::D) where D
    push!(subject.buffer, data)
    next!(subject.subject, data)
end

function on_error!(subject::ReplaySubjectInstance, err)
    error!(subject.subject, err)
end

function on_complete!(subject::ReplaySubjectInstance)
    complete!(subject.subject)
end

##

function on_subscribe!(subject::ReplaySubjectInstance, actor)
    for value in subject.buffer
        next!(actor, value)
    end
    return subscribe!(subject.subject, actor)
end

##

struct ReplaySubjectFactoryInstance{ F <: AbstractSubjectFactory } <: AbstractSubjectFactory
    factory :: F
    size    :: Int
end

Base.show(io::IO, subject::ReplaySubjectFactoryInstance{F}) where F = print(io, "ReplaySubjectFactoryInstance($F, size = $(subject.size))")

create_subject(::Type{L}, factory::ReplaySubjectFactoryInstance) where L = ReplaySubject(L, factory.size, factory.factory)

function ReplaySubjectFactory(size::Int, factory::F) where { F <: AbstractSubjectFactory }
    return ReplaySubjectFactoryInstance(factory, size)
end

function ReplaySubjectFactory(size::Int; scheduler::H = AsapScheduler()) where { H <: AbstractScheduler }
    return ReplaySubjectFactoryInstance(SubjectFactory{H}(scheduler), size)
end

##


text:
mutable struct DEOptions{SType,TstopType,CType}
    saveat::SType
    tstops::TstopType
    save_everystep::Bool
    dense::Bool
    timeseries_errors::Bool
    dense_errors::Bool
    save_on::Bool
    save_end::Bool
    callback::CType
    verbose::Bool
    advance_to_tstop::Bool
    stop_at_next_tstop::Bool
end

abstract type AbstractSundialsIntegrator <: DiffEqBase.AbstractODEIntegrator end

mutable struct CVODEIntegrator{uType,pType,memType,solType,algType,fType,UFType,JType,oType,toutType,sizeType,tmpType,LStype,Atype} <: AbstractSundialsIntegrator
    u::uType
    p::pType
    t::Float64
    tprev::Float64
    mem::memType
    LS::LStype
    A::Atype
    sol::solType
    alg::algType
    f::fType
    userfun::UFType
    jac::JType
    opts::oType
    tout::toutType
    tdir::Float64
    sizeu::sizeType
    u_modified::Bool
    tmp::tmpType
    uprev::tmpType
    flag::Cint
    just_hit_tstop::Bool
    event_last_time::Int
end

function (integrator::CVODEIntegrator)(t::Number,deriv::Type{Val{T}}=Val{0}) where T
    out = similar(integrator.u)
    integrator.flag = @checkflag CVodeGetDky(integrator.mem, t, Cint(T), out)
    out
end

function (integrator::CVODEIntegrator)(out,t::Number,
                                          deriv::Type{Val{T}}=Val{0}) where T
    integrator.flag = @checkflag CVodeGetDky(integrator.mem, t, Cint(T), out)
end

mutable struct ARKODEIntegrator{uType,pType,memType,solType,algType,fType,UFType,JType,oType,toutType,sizeType,tmpType,LStype,Atype} <: AbstractSundialsIntegrator
    u::uType
    p::pType
    t::Float64
    tprev::Float64
    mem::memType
    LS::LStype
    A::Atype
    sol::solType
    alg::algType
    f::fType
    userfun::UFType
    jac::JType
    opts::oType
    tout::toutType
    tdir::Float64
    sizeu::sizeType
    u_modified::Bool
    tmp::tmpType
    uprev::tmpType
    flag::Cint
    just_hit_tstop::Bool
    event_last_time::Int
end

function (integrator::ARKODEIntegrator)(t::Number,deriv::Type{Val{T}}=Val{0}) where T
    out = similar(integrator.u)
    integrator.flag = @checkflag ARKodeGetDky(integrator.mem, t, Cint(T), out)
    out
end

function (integrator::ARKODEIntegrator)(out,t::Number,
                                          deriv::Type{Val{T}}=Val{0}) where T
    integrator.flag = @checkflag ARKodeGetDky(integrator.mem, t, Cint(T), out)
end

mutable struct IDAIntegrator{uType,duType,pType,memType,solType,algType,fType,UFType,JType,oType,toutType,sizeType,sizeDType,tmpType,LStype,Atype} <: AbstractSundialsIntegrator
    u::uType
    du::duType
    p::pType
    t::Float64
    tprev::Float64
    mem::memType
    LS::LStype
    A::Atype
    sol::solType
    alg::algType
    f::fType
    userfun::UFType
    jac::JType
    opts::oType
    tout::toutType
    tdir::Float64
    sizeu::sizeType
    sizedu::sizeDType
    u_modified::Bool
    tmp::tmpType
    uprev::tmpType
    flag::Cint
    just_hit_tstop::Bool
    event_last_time::Int
end

function (integrator::IDAIntegrator)(t::Number,deriv::Type{Val{T}}=Val{0}) where T
    out = similar(integrator.u)
    integrator.flag = @checkflag IDAGetDky(integrator.mem, t, Cint(T), out)
    out
end

function (integrator::IDAIntegrator)(out,t::Number,
                                          deriv::Type{Val{T}}=Val{0}) where T
    integrator.flag = @checkflag IDAGetDky(integrator.mem, t, Cint(T), out)
end

###  Error check (retcode)

DiffEqBase.check_error(integrator::AbstractSundialsIntegrator) =
    interpret_sundials_retcode(integrator.flag)

DiffEqBase.postamble!(integrator::AbstractSundialsIntegrator) = nothing
# No-op postamble! to make check_error! (and hence iterator interface
# implemented in DiffEqBase) work.

### Iterator interface

@inline function DiffEqBase.step!(integrator::AbstractSundialsIntegrator)
  if integrator.opts.advance_to_tstop
    while integrator.tdir*(integrator.t-top(integrator.opts.tstops)) < -1e6eps()
        tstop = top(integrator.opts.tstops)
        set_stop_time(integrator,tstop)
        integrator.tprev = integrator.t
        if !(typeof(integrator.opts.callback.continuous_callbacks)<:Tuple{})
            integrator.uprev .= integrator.u
        end
        solver_step(integrator,tstop)
        integrator.t = first(integrator.tout)
        check_error!(integrator) != :Success && return
        handle_callbacks!(integrator)
        check_error!(integrator) != :Success && return
    end
 else
      integrator.tprev = integrator.t
      if !(typeof(integrator.opts.callback.continuous_callbacks)<:Tuple{})
          integrator.uprev .= integrator.u
      end
      if !isempty(integrator.opts.tstops)
          tstop = top(integrator.opts.tstops)
          set_stop_time(integrator,tstop)
          solver_step(integrator,tstop)
      else
          solver_step(integrator,1.0) # fake tstop
      end
      integrator.t = first(integrator.tout)
      check_error!(integrator) != :Success && return
      handle_callbacks!(integrator)
      check_error!(integrator) != :Success && return
  end
  handle_tstop!(integrator)
  nothing
end


text:
#
# Parrot.jl/src/prospect.jl
#
# By Elijah Mathews (me@elijahmathews.com)
# openpgp4fpr:ac1d3fb1e8a5eb7d14bd587b2932c725055a90d8
#

using PyCall

"""
    buildobs(...)

Generate a observation dictionary for the purposes of generating
training and validation data from [Prospector](https://github.com/bd-j/prospector).
"""
function buildobs(; snr=10, ldist=10, extras...)

    # The obs dictionary, empty for now.
    obs = Dict()

    # We will not use any obs filters for now.
    obs["filters"] = nothing

    # Since no filters, set "phot_wave" to nothing too.
    obs["phot_wave"] = nothing

    # No observations here, so set these to nothing.
    obs["wavelength"] = nothing
    obs["spectrum"] = nothing
    obs["unc"] = nothing
    obs["mask"] = nothing

    return obs

end

"""
    buildmodel(...)

Generate an SedModel object from [Prospector](https://github.com/bd-j/prospector)
for the purposes of generating training and validation data.
"""
function buildmodel(; object_redshift = nothing, ldist=10, fixed_metallicity = nothing, add_duste = false, extras...)

    # Required Python packages.
    sedmodel = pyimport("prospect.models.sedmodel")

    # Set template to parametric star formation history.
    py"""

    from prospect.models.templates import TemplateLibrary

    def get_template_library(key):

        return TemplateLibrary[key]

    """
    modelparams = py"get_template_library"("parametric_sfh")

    # Change to fixed stellar metallicity.
    modelparams["logzsol"]["isfree"] = false
    modelparams["logzsol"]["init"] = 0

    # Change to fixed redshift.
    modelparams["zred"]["isfree"] = false
    modelparams["zred"]["init"] = 0

    # Change to fixed stellar mass formed.
    modelparams["mass"]["isfree"] = false
    modelparams["mass"]["init"] = 1

    # Instantiate model object using modelparams.
    model = sedmodel.SedModel(modelparams)

    return model, modelparams

end

"""
    buildsps(...)

Generate a CSPSpecBasis object for the purposes of generating training
and validation data from [Prospector](https://github.com/bd-j/prospector).
"""
function buildsps(; zcontinuous = 1, extras...)

    sources = pyimport("prospect.sources")

    sps = sources.CSPSpecBasis(; zcontinuous=zcontinuous, extras...)

    return sps

end

"""
    generatedata(amount::Integer)

Generate training and validation data from [Prospector](https://github.com/bd-j/prospector)
(using [FSPS](https://github.com/cconroy20/fsps)). Returns a tuple:

    (λarray, params, Iarray)

where `λarray` is an `amount × 1` array consisting of the wavelengths in the spectra
(where `M` is the number of wavelengths), `params` is an `amount × 3` array consisting
of stellar population parameters (first column is the diffuse dust coefficient
`dust2`, second column is the universe age at galaxy's lookback time `tage`, and
the third column is the star formation timescale `tau` for an exponentially
declining star formation history), and `Iarray` is an `amount × M` array consisting
of the `log10` spectra in units of maggies.
"""
function generatedata(amount::Integer)

    # Generate requisite model objects.
    obs = buildobs()
    model, modelparams = buildmodel()
    sps = buildsps()

    # Grab wavelengths.
    λarray = sps.wavelengths

    # Initialize output arrays.
    farray = zeros(amount, length(λarray))
    params = zeros(amount, 3)

    # Loop to generate data.
    for i in 1:amount

        # Sample the stellar population parameters and save them.
        params[i,1] = modelparams["dust2"]["prior"].sample()[1]
        params[i,2] = modelparams["tage"]["prior"].sample()[1]
        params[i,3] = modelparams["tau"]["prior"].sample()[1]

        # Construct the array θ of stellar population parameters.
        θ = params[i,:]

        # Generate the spectrum and other stuff.
        spec, _, _ = model.sed(θ; obs = obs, sps = sps)

        # Take log10(...) of the spetrum.
        farray[i,:] = log10.(spec)

    end

    return λarray, params, farray

end


text:
using Statistics
function BoundByQuantiles!(data, margin=0.005)
for j = 1:size(data,2)
    q = quantile(data[:,j],margin)
    data[:,j] = max.(q,data[:,j])
    q = quantile(data[:,j],1.0-margin)
    data[:,j] = min.(q,data[:,j])
end
end

max_stars_repo_name: Vkamalvanshi/Econometrics

max_issues_repo_name: Vkamalvanshi/Econometrics

max_forks_repo_name: Vkamalvanshi/Econometrics


text:
"""
Calculate the moving average on an array of b-values
"""
function  _moving_average(arrb::AbstractVector{Float64}, nwin::Int)

    bavrg = fill(NaN,length(arrb))

    for i in 1:(length(arrb) - (nwin - 1))
        bavrg[i] = mean( @view arrb[i:(i+nwin-1)] )
    end

    return bavrg
end

"""
    bMcstep(mag, bincenters, ΔM, minevents)

Calculate b-value for each possible cut-off magnitude

# Arguments
*  `mags`: array storing the magnitude values of an incomplete catalog.
* `bincenters`: array storing the central value of the magnitude bins.
* `ΔM`: bin width.
* `minevents`: minimum number of events allowed when calculating b-values.

# Return
* `Bvals`: 1D array with the b values calculated varying the cut-off magnitude
* `Stds`: 1D array with the uncertainties associated with the calculated b values
"""
function _bMcstep(mag::AbstractVector{T}, bincenters::AbstractVector{Float64}, ΔM::S, minevents::Int) where {T,S<:Real}

    N = length(bincenters)
    bvals = fill(NaN,N)
    stds  = fill(NaN,N)

    for i in 1:N
        Mc = bincenters[i]
        submag = @view mag[mag .>= Mc-ΔM/2]

        if length(submag) >= minevents
            fitgr = fit_mle(GRclassic, submag, Mc-ΔM/2)
            bvals[i] = bvalue(fitgr)
            stds[i] = fit_stderr(fitgr, submag)
        else
            break
        end
    end

    return bvals, stds
end

"""
    mbs(mags, ΔM; winsize=5, minevents=50)

Compute the magnitude of completeness with the b−value stability (MBS)
method.

# Arguments
* `mags`: array storing the magnitude values of an incomplete catalog.
* `ΔM`: bin width to be used to group the magnitude values.`

# Keyword arguments
* `winsize`: window length used for the moving average. Default is 5 magnitude bins
* `minevents`: minimum number of events allowed when calculating b-values. Defaul is
   50 events.

# References:
* Cao and Gao, 2002, Temporal variations of seismic b-values beneath northeastern
  japan island arc, Geophysical Research Letters, 29, https://doi.org/10.1029/2001GL013775
* Woessner and Wiemer (2005), Assessing the quality of earthquake catalogues:
  Estimating the magnitude of completeness and its uncertainty, Bulletin of the Seismological Society of America,
  95, https://doi.org/10.1785/0120040007
"""
function mbs(mag::AbstractVector{T}, ΔM::S; winsize::Int=5, minevents::Int=50) where {T,S<:Real}

    bincenters = histmag(mag,ΔM)[1]

    # Calculate b-with magnitude
    bvals, stds = _bMcstep(mag, bincenters, ΔM, minevents)

    baverages = _moving_average(bvals,winsize)

    # If fBave is in in between the error estimate of the b-value of the first
    # cut-off magnitude, take it as guess
    goodSteps = abs.(baverages .- bvals) .<= stds
    latest_good = findfirst(goodSteps)

    if !isnothing(latest_good)
        return bincenters[latest_good]
    else
        return NaN
    end

end


text:
#=
gaussian_smoothing:
- Julia version: 1.1.0
- Author: Jirik
- Date: 2019-04-08
=#
using LarSurf
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm, SparseArrays
using Pandas
using Seaborn

V = [[0 0]; [2 0]; [1 2]; [3 2]; [4 2]; [4 0] ]'
EVch = spzeros(Int8, 7, 6)
EVch[1,1] = 1
EVch[1,2] = 1
EVch[2,2] = 1
EVch[2,3] = 1
EVch[3,1] = 1
EVch[3,3] = 1
EVch[4,3] = 1
EVch[4,4] = 1
EVch[5,2] = 1
EVch[5,4] = 1
EVch[6,2] = 1
EVch[6,6] = 1
EVch[7,4] = 1
EVch[7,6] = 1
Matrix(EVch)

FVch = spzeros(Int8, 3, 6)
FVch[1, 1] = 1
FVch[1, 2] = 1
FVch[1, 3] = 1
FVch[2, 2] = 1
FVch[2, 3] = 1
FVch[2, 4] = 1
# FVch[3, 2] = 1
# FVch[3, 4] = 1
# FVch[3, 6] = 1

FEch = spzeros(Int8, 3, 7)
FEch[1, 1] = 1
FEch[1, 2] = 1
FEch[1, 3] = 1
FEch[2, 2] = 1
FEch[2, 4] = 1
FEch[2, 5] = 1

function eye(sz)
    arr = spzeros(Int8, sz, sz)
    for i=1:sz
        arr[i,i] = 1
    end
    return arr
end

function setIzero!(arr)
    for i=1:size(arr)[1]
        arr[i,i] = 0
    end
    return arr
end

function getDiag(arr)
    sz = size(arr)[1]
    arro = zeros(eltype(arr), sz)
    for i=1:sz
        arro[i] = arr[i,i]
    end
    return arro
end

function make_full(arr)
    arrout = Array{eltype(arr), 2}(undef, size(arr)...)
    for i=1:size(arr)[1]
        for j=1:size(arr)[2]
            arrout[i,j] = arr[i,j]
        end
    end
    return arrout
end
#neighboor vertexes

function smoothing(V, EVch, k=0.35)
#     Matrix(setIzero!(EVch' * EVch))

    neighboorNumber = getDiag(EVch' * EVch)
    neighboors = setIzero!(EVch' * EVch)

#     targetV =  neighboors * V' ./ neighboorNumber
#     diff = targetV - V'
#     newV = (V' + k * diff)'

    targetV = V * neighboors ./ neighboorNumber'
    println("targetV shape: ", size(targetV))
    diff = targetV - V
    println("diff shape: ", size(diff))
    newV = V + k * diff

    return make_full(newV)
end

smoothing(V, EVch)

## Read data from file
# xystep = 1
# zstep = 1
# threshold = 4000;
# pth = Pio3d.datasets_join_path("medical/orig/sample-data/nrn4.pklz")

xystep = 10
zstep = 5
threshold = 10
pth = Pio3d.datasets_join_path("medical/orig/3Dircadb1.1/MASKS_DICOM/liver")

datap = Pio3d.read3d(pth)
data3d_full = datap["data3d"]
data3d = data3d_full[1:zstep:end, 1:xystep:end, 1:xystep:end];
segmentation = data3d .> threshold
println("File read finished, data size: ", size(segmentation))

tim = time()
block_size = [5,5,5]
filtered_bigFV, Flin, (bigV, tmodel) = LarSurf.get_surface_grid_per_block(segmentation, block_size)
bigVV, bigEV, bigFV, bigCV = tmodel
println("Surface extracted")

tim_prev = tim
tim = time()
println("surface extracted in: ", tim - tim_prev)

# Plasm.View((bigV,[bigVV, bigEV, filtered_bigFV]))
# Plasm.View((bigV,[bigVV, filtered_bigFV]))

# Triangulation
filtered_triangulated_bigFV = Lar.quads2triangles(filtered_bigFV::Lar.Cells)::Lar.Cells

EV = []
FV = filtered_bigFV
FV = filtered_triangulated_bigFV

nVerticesPerFV = length(FV[1])
if nVerticesPerFV == 4
    for f in FV
        push!(EV, [[f[1],f[2]],[f[3],f[4]],  [f[1],f[3]],[f[2],f[4]]])
    end
elseif nVerticesPerFV == 3
    for f in FV
        push!(EV, [[f[1],f[2]], [f[2],f[3]], [f[1],f[3]]])
    end
end
doubleedges = sort(cat(EV))
doubleedges = convert(Lar.Cells, doubleedges)
EV = [doubleedges[k] for k=1:2:length(doubleedges)]
aEV = LarSurf.ll2array(EV)
kEV = LarSurf.characteristicMatrix(aEV, size(bigV)[2])

# nrn4

# notAllBigEV =

# kEV = Lar.characteristicMatrix(EV)
# newBigV = smoothing(bigV, kEV, 0.3)
# Plasm.View((newBigV * 100,[filtered_bigFV]))
# Plasm.View((newBigV * 100,[bigVV, bigEV, filtered_bigFV]))
tim = time()

newBigV = bigV
newBigV = smoothing(newBigV, kEV, 0.6)
newBigV = smoothing(newBigV, kEV, -0.7)

tim_prev = tim
tim = time()
println("smoothing done in: ", tim - tim_prev)

# Plasm.View((newBigV,[filtered_bigFV]))
Plasm.View((newBigV,[filtered_triangulated_bigFV]))


text:
@testset "94.binary-tree-inorder-traversal.jl" begin
    @test preorder_traversal(TreeNode{Int}([1])) == [1]
    @test preorder_traversal(nothing) == Int[]
    @test preorder_traversal(TreeNode{Int}([1, nothing, 2, 3])) == [1, 2, 3]
    @test preorder_traversal(TreeNode{Int}([1, 2, 3, 4, 5, 6])) == [1, 2, 4, 5, 3, 6]
    @test preorder_traversal(TreeNode{Int}([3, 9, 4, nothing, nothing, 5, 7])) ==
        [3, 9, 4, 5, 7]
end


text:
module MarketTechnicals

using Reexport
using Statistics
using StatsBase

@reexport using TimeSeries

# movingaverages.jl
export sma, ema, kama, env
# levels.jl
export floorpivots, woodiespivots
# momentum.jl
export rsi, macd, cci, roc, adx, stochasticoscillator, chaikinoscillator, aroon, vortex,
       trix, massindex, dpo, kst, ichimoku, moneyflowindex, tsi, ultimateoscillator, williamsr,
       awesomeoscillator
# candlesticks.jl
export doji
# utilities.jl
export typical
# volatility.jl
export bollingerbands, truerange, atr, keltnerbands, chaikinvolatility, donchianchannels
# volume.jl
export obv, vwap, adl, chaikinmoneyflow, forceindex, easeofmovement, volumepricetrend

include("utilities.jl")
include("candlesticks.jl")
include("levels.jl")
include("movingaverages.jl")
include("momentum.jl")
include("volatility.jl")
include("volume.jl")

# for user customization
RC_FILE = joinpath(@__DIR__, ".rc.jl")
if !ispath(RC_FILE)
  touch(RC_FILE)
end

include(RC_FILE)

end


text:
abstract type AbstractResnet end
struct Resnet50 <: AbstractResnet end
struct Resnet200 <: AbstractResnet end

#####
##### Resnet 50 from Metalhead
#####

struct ResidualBlock
  conv_layers
  norm_layers
  shortcut
end

Flux.@functor ResidualBlock

function ResidualBlock(
        filters,
        kernels::Array{Tuple{Int,Int}},
        pads::Array{NTuple{4,Int}},
        strides::Array{Tuple{Int,Int}},
        shortcut = identity
    )
    conv_layers = [
        Conv(kernels[i-1], filters[i-1] => filters[i], pad = pads[i-1], stride = strides[i-1], init = Flux.glorot_normal)
        for i in 2:length(filters)
    ]

    norm_layers = [BatchNorm(filters[i]) for i in 2:length(filters)]
    ResidualBlock(Tuple(conv_layers), Tuple(norm_layers), shortcut)
end

function ResidualBlock(
        filters,
        kernels::Array{Int},
        pads::Array{Int},
        strides::Array{Int},
        shortcut = identity
    )
    ResidualBlock(
        filters,
        [(i,i) for i in kernels],
        [(i,i,i,i) for i in pads],
        [(i,i) for i in strides],
        shortcut
    )
end

function (block::ResidualBlock)(input)
    value = input
    for i in 1:length(block.conv_layers)-1
        value = relu.((block.norm_layers[i])((block.conv_layers[i])(value)))
    end
    relu.(((block.norm_layers[end])((block.conv_layers[end])(value))) + block.shortcut(input))
end

function Bottleneck(filters::Int, downsample::Bool = false, res_top::Bool = false)
    if(!downsample && !res_top)
        return ResidualBlock(
            [4 * filters, filters, filters, 4 * filters],
            [1,3,1],
            [0,1,0],
            [1,1,1]
        )
    elseif(downsample && res_top)
        return ResidualBlock(
            [filters, filters, filters, 4 * filters],
            [1,3,1],
            [0,1,0],
            [1,1,1],
            Chain(
                Conv(
                    (1,1),
                    filters=>4 * filters,
                    pad = (0,0,0,0),
                    stride = (1,1),
                    init = Flux.glorot_normal,
                ),
                BatchNorm(4 * filters)
           )
        )
    else
        shortcut = Chain(
            Conv((1,1), 2 * filters=>4 * filters, pad = (0,0,0,0), stride = (2,2), init = Flux.glorot_normal),
            BatchNorm(4 * filters)
        )
        return ResidualBlock(
            [2 * filters, filters, filters, 4 * filters],
            [1,3,1],
            [0,1,0],
            [1,1,2],
            shortcut
        )
    end
end

_layers(::Resnet50) = [3, 4, 6, 3]
_layers(::Resnet200) = [3, 24, 36, 3]

function _resnet(version::AbstractResnet)
    layers = _layers(version)
    layer_arr = []

    push!(layer_arr, Conv((7,7), 3=>64, pad = (3,3,3,3), stride = (2,2), init = Flux.glorot_normal))
    push!(layer_arr, x -> maxpool(x, (3,3), pad = (1,1), stride = (2,2)))

    initial_filters = 64
    for i in 1:length(layers)
        push!(layer_arr, Bottleneck(initial_filters, true, i==1))
        for j in 2:layers[i]
            push!(layer_arr, Bottleneck(initial_filters))
        end
        initial_filters *= 2
    end

    push!(layer_arr, x -> meanpool(x, (7,7)))
    push!(layer_arr, x -> reshape(x, :, size(x,4)))
    push!(layer_arr, (Dense(2048, 1000)))
    push!(layer_arr, x -> log.(max.(x, Float32(1e-9)))),
    push!(layer_arr, softmax)

    return UnstableChain(layer_arr)
end

function resnet_training(version::T, batchsize = 16; kw...) where {T <: AbstractResnet}
    Random.seed!(123455)
    X = rand(Float32, 224, 224, 3, batchsize)
    Y = rand(Float32, 1000, batchsize)

    f = ForwardLoss(_resnet(version), Flux.crossentropy)
    kw = (optimizer = nGraph.SGD(Float32(0.001)),)

    return Actualizer(f, X, Y; optimizer = nGraph.SGD(Float32(0.001)))
end

function resnet_inference(version::T, batchsize = 16) where {T <: AbstractResnet}
    X = rand(Float32, 224, 224, 3, batchsize)
    return Actualizer(_resnet(version), X)
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
#----------------------------------------------------------------------------------------------------------------------------------------
# Implements solution to the diffusion equation for the fluence in the steady-state and time-domain in a semi-infinite medium [1][2][3].
#
# [1] Alwin Kienle and Michael S. Patterson, "Improved solutions of the steady-state and the time-resolved diffusion equations
#     for reflectance from a semi-infinite turbid medium," J. Opt. Soc. Am. A 14, 246-254 (1997)
# [2] Martelli, Fabrizio, et al. Light propagation through biological tissue and other diffusive media: theory, solutions and software.
#     Vol. 10. No. 3.824746. Bellingham: SPIE press, 2010.
#----------------------------------------------------------------------------------------------------------------------------------------

""" Structure containing inputs for simulating the fluence under the diffusion approximation in the semi-infinite space."""
struct DAsemiinf_params{T <: AbstractFloat} <: DiffusionParameters
    ρ::T                                 # distance away from isotropic point source (cm)
    μa::T                                # absorption coefficient (cm⁻¹)
    μsp::T                               # reduced scattering coefficient (cm⁻¹)
    n_med::T                             # medium's index of refraction
    n_ext::T                             # external medium's index of refraction
    ω::T                                 # modulation frequency (1/ns)
    z::T                                 # z depth of detector within medium (cm)
    t::Union{T, AbstractVector{T}}       # time vector (ns)

    # do not need to provide these arguments (calculated from previous inputs)
    D::T                                 # Diffusion coefficient
    ν::T                                 # Speed of light in medium (cm/ns)
    z0::T                                # isotroptic source depth
    zb::T                                # Extrapolation length
    function DAsemiinf_params{T}(ρ::T, μa::T, μsp::T, n_med::T, n_ext::T, ω::T, z::T, t::Union{T, AbstractVector{T}}) where {T <: AbstractFloat}
        @assert ρ >= zero(T) "ρ must greater than or equal to 0"
        @assert μa >= zero(T) "μa must greater than or equal to 0"
        @assert μsp > zero(T) "μsp must greater than 0"
        @assert n_med > zero(T) "n_med must be positive"
        @assert n_ext > zero(T) "n_ext must be positive"
        @assert all(t .> zero(T)) "t must be positive"
        D = D_coeff(μsp)
        return new{T}(ρ, μa, μsp, n_med, n_ext, ω, z, t, D, ν_coeff(n_med), z0_coeff(μsp), zb_coeff(A_coeff(n_med / n_ext), D))
    end
end

"""
    Generator function for DAsemiinf_params structure.

Provides default parameters to use in the semi-infinite space fluence calculation in either the CW, FD, or TD.

# Keyword Arguments
- `ρ`: source-detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)
- `n_med`: medium's index of refraction
- `n_ext`: external medium's index of refraction (air or detector)
- `ω`: modulation frequency (1/ns)
- `z`: the z-depth orthogonal from the boundary (cm)
- `t`: the time vector (ns).

# Examples
```
julia> data = DAsemiinf_params() # return default parameters
julia> data = DAsemiinf_params(ρ = 1.5) # return ρ = 1.5 with the rest of the parameters given by defaults
julia> fluence_DA_semiinf_CW(data) # can then be used with corresponding functions
```
"""
function DAsemiinf_params(;
    ρ::T = 1.0,
    μa::T = 0.1,
    μsp::T = 10.0,
    n_med::T = 1.0,
    n_ext::T = 1.0,
    ω::T = 0.0,
    z::T = 0.0,
    t::Union{T, AbstractVector{T}} = 1.0
) where {T<:AbstractFloat}
    return DAsemiinf_params{T}(ρ, μa, μsp, n_med, n_ext, ω, z, t)
end

#------------------------------
# Steady-State Fluence
#------------------------------

"""
    fluence_DA_semiinf_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0.0)

Compute the steady-state fluence in a semi-infinite geometry.

# Arguments
- `ρ`: the source detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Keyword arguments
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `z`: the z-depth orthogonal from the boundary (cm)

# Examples
```
julia> fluence_DA_semiinf_CW(1.0, 0.1, 10.0)
julia> fluence_DA_semiinf_CW(1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, z = 0.0)
```
"""
function fluence_DA_semiinf_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0.0)
    params = DiffusionKernelParams(μsp, n_med, n_ext)
    μeff = sqrt(3 * μa * μsp)
    return _kernel_fluence_DA_semiinf_CW(μeff, ρ, z, params.z0, params.zb, params.D)
end
"""
    fluence_DA_semiinf_CW(data::DiffusionParameters)

Wrapper to `fluence_DA_semiinf_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0.0...)`

# Examples
```
julia> data = DAsemiinf_params(ρ = 1.0) # use structure to generate inputs
julia> fluence_DA_semiinf_CW(data) # then call the function
```
"""
function fluence_DA_semiinf_CW(data::DiffusionParameters)
    μeff = sqrt(3 * data.μa * data.μsp)
    return _kernel_fluence_DA_semiinf_CW(μeff, data.ρ, data.z, data.z0, data.zb, data.D)
end

function _kernel_fluence_DA_semiinf_CW(μeff, ρ, z, z0, zb, D)
    ϕ = exp(-μeff * sqrt(ρ^2 + (z - z0)^2)) / (sqrt(ρ^2 + (z - z0)^2))
    ϕ -= exp(-μeff * sqrt(ρ^2 + (z + 2 * zb + z0)^2)) / (sqrt(ρ^2 + (z + 2 * zb + z0)^2))
    return ϕ / (4 * π * D)
end

#------------------------------
# Time-Domain Fluence
#------------------------------

"""
    fluence_DA_semiinf_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0.0)

Compute the time-domain fluence in a semi-infinite medium.

# Arguments
- `t`: the time vector (ns).
- `ρ`: the source detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Optional arguments
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `z`: the z-depth orthogonal from the boundary (cm)

# Examples
```
julia> fluence_DA_semiinf_TD(0.1:0.1:2.0, 1.0, 0.1, 10.0)
julia> fluence_DA_semiinf_TD(0.2:0.2:1.0, 1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, z = 0.0)
```
"""
function fluence_DA_semiinf_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0.0)
    params = DiffusionKernelParams(μsp, n_med, n_ext)
    return map(t -> _kernel_fluence_DA_semiinf_TD(params.D, params.ν, t, μa, z, params.z0, ρ, params.zb), t)
end

"""
    fluence_DA_semiinf_TD(t, data::DiffusionParameters)

Wrapper to `fluence_DA_semiinf_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0.0)`

# Examples
```
julia> data = DAsemiinf_params() # use structure to generate inputs
julia> fluence_DA_semiinf_TD(data) # then call the function
```
"""
function fluence_DA_semiinf_TD(data::DiffusionParameters)
    return map(t -> _kernel_fluence_DA_semiinf_TD(data.D, data.ν, t, data.μa, data.z, data.z0, data.ρ, data.zb), data.t)
end

@inline function _kernel_fluence_DA_semiinf_TD(D, ν, t, μa, z, z0, ρ, zb)
    tmp1 = 4 * D * ν * t
    ϕ = ν * exp(-μa * ν * t) / (π * tmp1)^(3/2)
    ϕ *= (exp(-((z - z0)^2 + ρ^2) / tmp1) - exp(-((z + z0 + 2*zb)^2 + ρ^2) / tmp1))
    return ϕ
end

#------------------------------
# Frequency-Domain Fluence
#------------------------------
"""
fluence_DA_semiinf_FD(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0.0, ω = 1.0)

Compute the frequency domain fluence in a semi-infinite geometry.

# Arguments
- `ρ`: the source detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Keyword arguments
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `z`: the z-depth orthogonal from the boundary (cm)
- `ω`: the modulation frequency (1/ns)

# Examples
```
julia> fluence_DA_semiinf_FD(1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0, z = 0.0, ω = 1.0)
```
"""
function fluence_DA_semiinf_FD(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0.0, ω = 1.0)
    params = DiffusionKernelParams(μsp, n_med, n_ext)
    μa_complex = μa + ω * im / params.ν
    μeff = sqrt(3 * μa_complex * data.μsp)
    return _kernel_fluence_DA_semiinf_CW(μeff, ρ, z, params.z0, params.zb, params.D)
end
"""
    fluence_DA_semiinf_FD(data::DiffusionParameters)

Wrapper to `fluence_DA_semiinf_FD(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0.0, ω = 1.0)`

# Examples
```
julia> data = DAsemiinf_params(ω = 1.0) # use structure to generate inputs
julia> fluence_DA_semiinf_FD(data) # then call the function
```
"""
function fluence_DA_semiinf_FD(data::DiffusionParameters)
    μa_complex = data.μa + data.ω * im / data.ν
    μeff = sqrt(3 * μa_complex * data.μsp)
    return _kernel_fluence_DA_semiinf_CW(μeff, data.ρ, data.z, data.z0, data.zb, data.D)
end

#------------------------------
# Time-Domain Flux
#------------------------------
"""
    flux_DA_semiinf_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0)

Compute the time-domain flux (D*∂ϕ(t)/∂z @ z = 0) from a semi-infinite medium from Eqn. 36 Contini 97 (m = 0).

# Arguments
- `t`: the time vector (ns).
- `ρ`: the source detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Keyword arguments
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction

# Examples
```
julia> flux_DA_semiinf_TD(0.1:0.1:2.0, 1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0)
```
"""
function flux_DA_semiinf_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0)
    params = DiffusionKernelParams(μsp, n_med, n_ext)
    z3m = - params.z0
    z4m = 2 * params.zb + params.z0

    return map(t -> _kernel_flux_DA_semiinf_TD(params.D, params.ν, t, μa, z3m, z4m, ρ), t)
 end

 """
    flux_DA_semiinf_TD(t, data::DiffusionParameters)

Wrapper to `flux_DA_semiinf_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0.0)`

# Examples
```
julia> data = DAsemiinf_params() # use structure to generate inputs
julia> flux_DA_semiinf_TD(data) # then call the function
```
"""
function flux_DA_semiinf_TD(data::DiffusionParameters)
    z3m = - data.z0
    z4m = 2 * data.zb + data.z0

    return map(t -> _kernel_flux_DA_semiinf_TD(data.D, data.ν, t, data.μa, z3m, z4m, data.ρ), data.t)
end

@inline function _kernel_flux_DA_semiinf_TD(D, ν, t, μa, z3m, z4m, ρ)
    tmp1 = 4 * D * ν * t
    Rt = -(ρ^2 / (tmp1))
    Rt -= μa * ν * t
    Rt = -exp(Rt) / (2 * (4 * π * D * ν)^(3/2) * sqrt(t * t * t * t * t))
    Rt *= (z3m * exp(-(z3m^2 / (tmp1))) - z4m * exp(-(z4m^2 / (tmp1))))

    return Rt
end

#------------------------------
# Steady-State Flux
#------------------------------
"""
    flux_DA_semiinf_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0)

Compute the steady-state flux (D*∂ϕ(ρ)/∂z @ z = 0) from a semi-infinite medium.

# Arguments
- `ρ`: the source detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)
- `ndet`: the boundary's index of refraction (air or detector)
- `nmed`: the sample medium's index of refraction

# Examples
```
julia> flux_DA_semiinf_CW(1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0)
```
"""
function flux_DA_semiinf_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0)
    params = DiffusionKernelParams(μsp, n_med, n_ext)
    return params.D * ForwardDiff.derivative(z -> fluence_DA_semiinf_CW(ρ, μa, μsp, n_med = n_med, n_ext = n_ext, z = z), 0.0)
end

"""
    flux_DA_semiinf_CW(data::DiffusionParameters)

Wrapper to `flux_DA_semiinf_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0)`

# Examples
```
julia> data = DAsemiinf_params(ρ = 1.0) # use structure to generate inputs
julia> flux_DA_semiinf_CW(data) # then call the function
```
"""
function flux_DA_semiinf_CW(data)
    return flux_DA_semiinf_CW(data.ρ, data.μa, data.μsp, n_ext = data.n_ext, n_med = data.n_med)
end

max_issues_repo_name: heltonmc/LightPropagation


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Boundary{P,Q}(topology)

The boundary relation from rank `P` to smaller rank `Q` for
a given `topology`.
"""
struct Boundary{P,Q,D,T<:Topology} <: TopologicalRelation
  topology::T
end

function Boundary{P,Q}(topology) where {P,Q}
  D = paramdim(topology)
  T = typeof(topology)

  @assert D ≥ P > Q "invalid boundary relation"

  Boundary{P,Q,D,T}(topology)
end

# --------------
# FULL TOPOLOGY
# --------------

function (∂::Boundary{D,0,D,T})(ind::Integer) where {D,T<:FullTopology}
  collect(connec4elem(∂.topology, ind))
end

# --------------
# GRID TOPOLOGY
# --------------

# vertices of hexahedron on 3D grid
function (∂::Boundary{3,0,3,T})(ind::Integer) where {T<:GridTopology}
  t = ∂.topology
  i, j, k = elem2cart(t, ind)
  i1 = cart2corner(t, i  , j  , k  )
  i2 = cart2corner(t, i+1, j  , k  )
  i3 = cart2corner(t, i+1, j+1, k  )
  i4 = cart2corner(t, i  , j+1, k  )
  i5 = cart2corner(t, i  , j  , k+1)
  i6 = cart2corner(t, i+1, j  , k+1)
  i7 = cart2corner(t, i+1, j+1, k+1)
  i8 = cart2corner(t, i  , j+1, k+1)
  [i1, i2, i3, i4, i5, i6, i7, i8]
end

# vertices of quadrangle on 3D grid
function (∂::Boundary{2,0,3,T})(ind::Integer) where {T<:GridTopology}
  @error "not implemented"
end

# vertices of quadrangle on 2D grid
function (∂::Boundary{2,0,2,T})(ind::Integer) where {T<:GridTopology}
  t = ∂.topology
  i, j = elem2cart(t, ind)
  i1 = cart2corner(t, i  , j  )
  i2 = cart2corner(t, i+1, j  )
  i3 = cart2corner(t, i+1, j+1)
  i4 = cart2corner(t, i  , j+1)
  [i1, i2, i3, i4]
end

# vertices of segment on 2D grid
function (∂::Boundary{1,0,2,T})(ind::Integer) where {T<:GridTopology}
  t = ∂.topology
  nx, ny = size(t)

  if ind ≤ nx*ny # vertical edges
    i1 = elem2corner(t, ind)
    i2 = i1 + 1
  elseif ind ≤ nx*(ny+1) # last vertical edges
    i, j = elem2cart(t, ind - nx)
    i1 = cart2corner(t, i, j+1)
    i2 = i1 + 1
  else # horizontal edges
    i1 = ind - nx*(ny+1)
    i, j = corner2cart(t, i1)
    i2 = cart2corner(t, i, j+1)
  end

  [i1, i2]
end

# vertices of segment on 1D grid
function (∂::Boundary{1,0,1,T})(ind::Integer) where {T<:GridTopology}
  i1 = ind
  i2 = ind + 1
  [i1, i2]
end

# -------------------
# HALF-EDGE TOPOLOGY
# -------------------

function (∂::Boundary{2,1,2,T})(elem::Integer) where {T<:HalfEdgeTopology}
  t = ∂.topology
  l = loop(half4elem(t, elem))
  v = CircularVector(l)
  [edge4pair(t, (v[i], v[i+1])) for i in 1:length(v)]
end

function (∂::Boundary{2,0,2,T})(elem::Integer) where {T<:HalfEdgeTopology}
  loop(half4elem(∂.topology, elem))
end

function (∂::Boundary{1,0,2,T})(edge::Integer) where {T<:HalfEdgeTopology}
  e = half4edge(∂.topology, edge)
  [e.head, e.half.head]
end


text:
module TestPackage

include("Integration.jl")

export left_riemann
export right_riemann
export trapezoid

end


text:
norm(x) = sqrt(mapreduce(x->x*x, +, x))
synchronize(device::KA.Device) = wait(Event(device))


text:
using Test

macro timeinclude(x)
    :(print($x); print("\t"); @time include($x))
end

@testset "Knet" begin

#                            commit v147 v146 555 686 ee1 c6a 0b5 e3a 9.2 8.3 8.3 6cb 6cb 8.6 6cb 6cb
#                           machine v100 v100 dy3 rzr  ci rzr  sc  sc tig ai5 ai4 tr5 tr4 aws osx os4
@timeinclude("kptr.jl")          #    15   13  30   2   ?   4   3   8  16   1   1   0   0  20   0   0
@timeinclude("gpu.jl")           #     9    7  31   2   3   1   2   2   6   1   1   0   0   2   0   0
@timeinclude("distributions.jl") #     1    1  26   1   1   1   1   1   2   1   1   2   1   3   3   2
@timeinclude("dropout.jl")       #    17    6  31   3   6   1   2   8   5                   2
@timeinclude("gcnode.jl")        #     5    5  39   1
@timeinclude("jld.jl")           #    12   12  45   2  11  11   9  26
@timeinclude("statistics.jl")    #    18   15  48   7
@timeinclude("bmm.jl")           #    41   36  65   9  10   3   9
@timeinclude("serialize.jl")     #     7    6  70  20  16   2  11   1
@timeinclude("loss.jl")          #    85   40  74  17  32  13  19  20  10                   4
@timeinclude("cuarray.jl")       #    20   21  74  18  25
@timeinclude("update.jl")        #    54   26  78  64 316  61  82  60  61  29  26 100  22  72  25  23
@timeinclude("linalg.jl")        #   101   64  85  17  33  17  22  22  62  24  14  22   7  28  33  19
@timeinclude("batchnorm.jl")     #    47   33  92  19  32  14  23  22  93
@timeinclude("ops20.jl")         #    82   41  94
@timeinclude("karray.jl")        #   179   84 114  27  49  20  27  21  55  19  12   -   -  21   -   0
@timeinclude("conv.jl")          #   112   59 114  47  71  33  48  51 107  22  12  62  47  26  44  16
@timeinclude("rnn.jl")           #    40   39 126  43  95  52  54  41  81                  12
@timeinclude("reduction.jl")     #   132   65 137  42  67  37  52  55 106  40  21  29  11  57  55  29
@timeinclude("unary.jl")         #   627  225 174  76 172  86 103 103 122  42   6  36   4  56  67  11
@timeinclude("binary.jl")        #   222  164 234  89 124  64  75  80  56  34  19 491 119  51  53  25
#TODO include("data.jl")
#TODO include("hyperopt.jl")
#TODO include("progress.jl")
#TODO include("random.jl")
#TODO include("train.jl")
#TODO include("uva.jl")

end


text:
export Clock

struct Clock
	interval::Real
	fixed::Bool
	start_running::Bool
	Clock(interval=1, fixed=false, start_running=false) = interval >= 0 ? new(interval, fixed, start_running) : error("interval must be non-negative")
end

function show(io::IO, ::MIME"text/html", clock::Clock)
    # We split the HTML string into multiple files, but you could also write all of this into a single (long) string 🎈
	cb = read(joinpath(PKG_ROOT_DIR, "assets", "clock_back.svg"), String)
	cf = read(joinpath(PKG_ROOT_DIR, "assets", "clock_front.svg"), String)
	cz = read(joinpath(PKG_ROOT_DIR, "assets", "clock_zoof.svg"), String)
	js = read(joinpath(PKG_ROOT_DIR, "assets", "clock.js"), String)
	css = read(joinpath(PKG_ROOT_DIR, "assets", "clock.css"), String)

	result = """
	<clock class='$(clock.fixed ? " fixed" : "")$(clock.start_running ? "" : " stopped")'>
		<analog>
			<back>$(cb)</back>
			<front>$(cf)</front>
			<zoof style="opacity: 0">$(cz)</zoof>
		</analog>
		<button></button>
		speed:
		<input type="number" value="$(clock.interval)"  min=0 step=any lang="en-001">

	</clock>

    """
    write(io, result)
end

get(clock::Clock) = 1


text:
# testing multiple agents coming together
# 2 agents considered are cubeprop systems
# they attach to a box and fly together
include("../src/plotSol.jl")
include("../src/simulateDiff.jl")
include("../src/OrientationConversion.jl")

clearconsole()
j1 = Nothing
j2 = Nothing
j3 = Nothing
## Agent 1
mb = 0.155; mp = 0.02;

Ib =  [ 0.0009    0.0000   -0.0000;
        0.0000    0.0053    0.0000;
       -0.0000    0.0000    0.0054]#*(10^-3)

Ip =    [ 0.0029    0.0000    0.0000;
          0.0000    0.5824    0.0000;
          0.0000    0.0000    0.5805]*(10^-4)

InFrame = InertialFrameAsRB()
QuadCube_1 = RigidBody(mb,Ib,2)
Props_1 = RigidBody(mp,Ip,3)

x0Cube_1 = [[-0.5;0.0;0.0];[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(QuadCube_1,x0Cube_1)
x0Props_1 = [[-0.5;0.0;0.5];[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(Props_1,x0Props_1)

axis = [0.0 0.0 1.0][:]
rjCube = [0.0 0.0 0.5][:]
rjProp = [0.0 0.0 0.0][:]

j1_1 = Joint(InFrame,QuadCube_1,zeros(3),[0.5;zeros(2)])
j1_2 = Joint(QuadCube_1,Props_1,rjCube,rjProp,type = "Revolute",axis = axis);#,jointTorque = [0.0;0.0;0.0001])

j1 = [j1_1,j1_2];

## Agent 2
QuadCube_2 = RigidBody(mb,Ib,2)
Props_2 = RigidBody(mp,Ip,3)

x0Cube_2 = [[0.5;0.0;0.0];[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(QuadCube_2,x0Cube_2)
x0Props_2 = [[0.5;0.0;0.5];[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(Props_2,x0Props_2)

axis = [0.0 0.0 1.0][:]
rjCube = [0.0 0.0 0.5][:]
rjProp = [0.0 0.0 0.0][:]

j2_1 = Joint(InFrame,QuadCube_2,[-0.5;zeros(2)],zeros(3))
j2_2 = Joint(QuadCube_2,Props_2,rjCube,rjProp,type = "Revolute",axis = axis);#,jointTorque = [0.0;0.0;0.0001])

j2 = [j2_1,j2_2];
## Box
Box = RigidBody(mb,Ib,2);
x0Box = [zeros(3);[1;zeros(3)];zeros(3);zeros(4)]; # Box cm is at the inertial frame
initialiseRigidBody!(Box,x0Box);

j3_1 = Joint(InFrame, Box, zeros(3), zeros(3));

j3 = [j3_1];
## Joining bodies together
# join cubeprop with a box
joinBodies!(j1,j3,2,2,[0.25;zeros(2)],[-0.25;zeros(2)]);
# println(length(j1))
println(j1[end].RB2.bodyID)
# join cubepropbox with another cubeprop
joinBodies!(j1,j2,4,2,[0.25;zeros(2)],[-0.25;zeros(2)]);
println(j1[end].RB2.bodyID)

# joining two cubeprops together
# joinBodies!(j1,j2,2,2,[0.5;zeros(2)],[-0.5;zeros(2)]);
# println(j1[3].RB2.bodyID)
# println(j1[5].RB2.bodyID)

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
###
### define wrapper with ArrayInterface.dimnames
###
@testset "dimension permutations" begin
    a = ones(2, 2, 2)
    perm = PermutedDimsArray(a, (3, 1, 2))
    mview = view(perm, :, 1, :)
    madj = mview'
    vview = view(madj, 1, :)
    vadj = vview'

    @test @inferred(ArrayInterface.to_parent_dims(typeof(a))) == (1, 2, 3)
    @test @inferred(ArrayInterface.to_parent_dims(typeof(perm))) == (3, 1, 2)
    @test @inferred(ArrayInterface.to_parent_dims(typeof(mview))) == (1, 3)
    @test @inferred(ArrayInterface.to_parent_dims(typeof(madj))) == (2, 1)
    @test @inferred(ArrayInterface.to_parent_dims(typeof(vview))) == (2,)
    @test @inferred(ArrayInterface.to_parent_dims(typeof(vadj))) == (2, 1)
    @test @inferred(ArrayInterface.to_parent_dims(typeof(vadj), static(1))) == 2
    @test @inferred(ArrayInterface.to_parent_dims(typeof(vadj), 1)) == 2
    @test @inferred(ArrayInterface.to_parent_dims(typeof(vadj), static(3))) == 2
    @test @inferred(ArrayInterface.to_parent_dims(typeof(vadj), 3)) == 2

    @test @inferred(ArrayInterface.from_parent_dims(a)) == (1, 2, 3)
    @test @inferred(ArrayInterface.from_parent_dims(typeof(perm))) == (2, 3, 1)
    @test @inferred(ArrayInterface.from_parent_dims(typeof(mview))) == (1, 0, 2)
    @test @inferred(ArrayInterface.from_parent_dims(typeof(madj))) == (2, 1)
    @test @inferred(ArrayInterface.from_parent_dims(typeof(vview))) == (0, 1)
    @test @inferred(ArrayInterface.from_parent_dims(typeof(vadj))) == (2,)
    @test @inferred(ArrayInterface.from_parent_dims(typeof(vadj), static(1))) == 2
    @test @inferred(ArrayInterface.from_parent_dims(typeof(vadj), 1)) == 2

    @test_throws DimensionMismatch ArrayInterface.to_parent_dims(typeof(vadj), 0)
    @test_throws DimensionMismatch ArrayInterface.to_parent_dims(typeof(vadj), static(0))

    @test_throws DimensionMismatch ArrayInterface.from_parent_dims(typeof(vadj), 0)
    @test_throws DimensionMismatch ArrayInterface.from_parent_dims(typeof(vadj), static(0))

    colormat = reinterpret(reshape, Float64, [(R=rand(), G=rand(), B=rand()) for i ∈ 1:100])
    @test @inferred(ArrayInterface.from_parent_dims(typeof(colormat))) === (static(2),)
    @test @inferred(ArrayInterface.to_parent_dims(typeof(colormat))) === (static(0), static(1),)

    Rr = reinterpret(reshape, Int32, ones(4))
    @test @inferred(ArrayInterface.from_parent_dims(typeof(Rr))) === (static(2),)
    @test @inferred(ArrayInterface.to_parent_dims(typeof(Rr))) === (static(0), static(1),)

    Rr = reinterpret(reshape, Int64, ones(4))
    @test @inferred(ArrayInterface.from_parent_dims(typeof(Rr))) === (static(1),)
    @test @inferred(ArrayInterface.to_parent_dims(typeof(Rr))) === (static(1),)

    Sr = reinterpret(reshape, Complex{Int64}, zeros(2, 3, 4))
    @test @inferred(ArrayInterface.from_parent_dims(typeof(Sr))) === (static(0), static(1), static(2))
    @test @inferred(ArrayInterface.to_parent_dims(typeof(Sr))) === (static(2), static(3))
end

@testset "order_named_inds" begin
    n1 = (static(:x),)
    n2 = (n1..., static(:y))
    n3 = (n2..., static(:z))
    @test @inferred(ArrayInterface.find_all_dimnames(n1, (), (), :)) == ()
    @test @inferred(ArrayInterface.find_all_dimnames(n1, (static(:x),), (2,), :)) == (2,)
    @test @inferred(ArrayInterface.find_all_dimnames(n2, (static(:x),), (2,), :)) == (2, :)
    @test @inferred(ArrayInterface.find_all_dimnames(n2, (static(:y),), (2,), :)) == (:, 2)
    @test @inferred(ArrayInterface.find_all_dimnames(n2, (static(:y), static(:x)), (20, 30), :)) == (30, 20)
    @test @inferred(ArrayInterface.find_all_dimnames(n2, (static(:x), static(:y)), (30, 20), :)) == (30, 20)
    @test @inferred(ArrayInterface.find_all_dimnames(n3, (static(:x), static(:y)), (30, 20), :)) == (30, 20, :)

    @test_throws ErrorException ArrayInterface.find_all_dimnames(n2, (static(:x), static(:y), static(:z)), (30, 20, 40), :)
end

@testset "ArrayInterface.dimnames" begin
    d = (static(:x), static(:y))
    x = NamedDimsWrapper(d, ones(Int64, 2, 2))
    y = NamedDimsWrapper((static(:x),), ones(2))
    z = NamedDimsWrapper((:x, static(:y)), ones(2))
    r1 = reinterpret(Int8, x)
    r2 = reinterpret(reshape, Int8, x)
    r3 = reinterpret(reshape, Complex{Int}, x)
    r4 = reinterpret(reshape, Float64, x)
    w = Wrapper(x)
    dnums = ntuple(+, length(d))
    @test @inferred(ArrayInterface.has_dimnames(x)) == true
    @test @inferred(ArrayInterface.has_dimnames(z)) == true
    @test @inferred(ArrayInterface.has_dimnames(ones(2, 2))) == false
    @test @inferred(ArrayInterface.has_dimnames(Array{Int,2})) == false
    @test @inferred(ArrayInterface.has_dimnames(typeof(x))) == true
    @test @inferred(ArrayInterface.has_dimnames(typeof(view(x, :, 1, :)))) == true
    @test @inferred(ArrayInterface.dimnames(x)) === d
    @test @inferred(ArrayInterface.dimnames(w)) === d
    @test @inferred(ArrayInterface.dimnames(r1)) === d
    @test @inferred(ArrayInterface.dimnames(r2)) === (static(:_), d...)
    @test @inferred(ArrayInterface.dimnames(r3)) === Base.tail(d)
    @test @inferred(ArrayInterface.dimnames(r4)) === d
    @test @inferred(ArrayInterface.ArrayInterface.dimnames(z)) === (:x, static(:y))
    @test @inferred(ArrayInterface.dimnames(parent(x))) === (static(:_), static(:_))
    @test @inferred(ArrayInterface.dimnames(reshape(x, (1, 4)))) === d
    @test @inferred(ArrayInterface.dimnames(reshape(x, :))) === (static(:_),)
    @test @inferred(ArrayInterface.dimnames(x')) === reverse(d)
    @test @inferred(ArrayInterface.dimnames(y')) === (static(:_), static(:x))
    @test @inferred(ArrayInterface.dimnames(PermutedDimsArray(x, (2, 1)))) === reverse(d)
    @test @inferred(ArrayInterface.dimnames(PermutedDimsArray(x', (2, 1)))) === d
    @test @inferred(ArrayInterface.dimnames(view(x, :, 1))) === (static(:x),)
    @test @inferred(ArrayInterface.dimnames(view(x, :, 1)')) === (static(:_), static(:x))
    @test @inferred(ArrayInterface.dimnames(view(x, :, :, :))) === (static(:x), static(:y), static(:_))
    @test @inferred(ArrayInterface.dimnames(view(x, :, 1, :))) === (static(:x), static(:_))
    @test @inferred(ArrayInterface.dimnames(x, ArrayInterface.One())) === static(:x)
    @test @inferred(ArrayInterface.dimnames(parent(x), ArrayInterface.One())) === static(:_)
    @test @inferred(ArrayInterface.known_dimnames(Iterators.flatten(1:10))) === (:_,)
    @test @inferred(ArrayInterface.known_dimnames(Iterators.flatten(1:10), static(1))) === :_
    @test @inferred(ArrayInterface.known_dimnames(z)) === (nothing, :y)
    @test @inferred(ArrayInterface.known_dimnames(reshape(x, (1, 4)))) === (:x, :y)
    @test @inferred(ArrayInterface.known_dimnames(r1)) === (:x, :y)
    @test @inferred(ArrayInterface.known_dimnames(r2)) === (:_, :x, :y)
    @test @inferred(ArrayInterface.known_dimnames(r3)) === (:y,)
    @test @inferred(ArrayInterface.known_dimnames(r4)) === (:x, :y)
    @test @inferred(ArrayInterface.known_dimnames(w)) === (:x, :y)
    @test @inferred(ArrayInterface.known_dimnames(reshape(x, :))) === (:_,)
    @test @inferred(ArrayInterface.known_dimnames(view(x, :, 1)')) === (:_, :x)
end

@testset "to_dims" begin
    x = NamedDimsWrapper(static((:x, :y)), ones(2, 2))
    y = NamedDimsWrapper(static((:x, :y, :a, :b, :c, :d)), ones(6))

    @test @inferred(ArrayInterface.to_dims(x, :)) == Colon()
    @test @inferred(ArrayInterface.to_dims(x, 1)) == 1
    @testset "small case" begin
        @test @inferred(ArrayInterface.to_dims(x, (:x, :y))) == (1, 2)
        @test @inferred(ArrayInterface.to_dims(x, (:y, :x))) == (2, 1)
        @test @inferred(ArrayInterface.to_dims(x, :x)) == 1
        @test @inferred(ArrayInterface.to_dims(x, :y)) == 2
        @test_throws DimensionMismatch ArrayInterface.to_dims(x, static(:z))  # not found
        @test_throws DimensionMismatch ArrayInterface.to_dims(x, :z)  # not found
    end

    @testset "large case" begin
        @test @inferred(ArrayInterface.to_dims(y, :x)) == 1
        @test @inferred(ArrayInterface.to_dims(y, :a)) == 3
        @test @inferred(ArrayInterface.to_dims(y, :d)) == 6
        @test_throws DimensionMismatch ArrayInterface.to_dims(y, :z) # not found
    end
end

@testset "methods accepting ArrayInterface.dimnames" begin
    d = (static(:x), static(:y))
    x = NamedDimsWrapper(d, ones(2, 2))
    y = NamedDimsWrapper((static(:x),), ones(2))
    @test @inferred(size(x, first(d))) == size(parent(x), 1)
    @test @inferred(ArrayInterface.size(y')) == (1, size(parent(x), 1))
    @test @inferred(axes(x, first(d))) == axes(parent(x), 1)
    @test strides(x, :x) == ArrayInterface.strides(parent(x))[1]
    @test @inferred(ArrayInterface.axes_types(x, static(:x))) <: Base.OneTo{Int}
    @test ArrayInterface.axes_types(x, :x) <: Base.OneTo{Int}
    @test @inferred(ArrayInterface.axes_types(LinearIndices{2,NTuple{2,Base.OneTo{Int}}})) <: NTuple{2,Base.OneTo{Int}}
    CI = CartesianIndices{2,Tuple{Base.OneTo{Int},UnitRange{Int}}}
    @test @inferred(ArrayInterface.axes_types(CI, static(1))) <: Base.OneTo{Int}

    x[x=1] = [2, 3]
    @test @inferred(getindex(x, x=1)) == [2, 3]
    y = NamedDimsWrapper((:x, static(:y)), ones(2, 2))
    # FIXME this doesn't correctly infer the output because it can't infer
    @test getindex(y, x=1) == [1, 1]
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   This file contains an intermediate representation of the state vector. RV is
#   a tuple with the position and velocity. This functions are used when
#   converting a state vector to Keplerian elements and vice-versa.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Schwarz, R (2014). Memorandum No. 2: Cartesian State Vectors to
#       Keplerian Orbit Elements. Available at www.rene-schwarz.com
#
#       https://downloads.rene-schwarz.com/dc/category/18
#       (Accessed on 2017-08-09).
#
#   [2] Vallado, D. A., McClain, W. D (2013). Fundamentals of Astrodynamics
#       and Applications. Microcosm Press.
#
#   [3] Kuga, H. K., Carrara, V., Rao, K. R (2005). Introdução à Mecânica
#       Orbital. 2ª ed. Instituto Nacional de Pesquisas Espaciais.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export kepler_to_rv, rv_to_kepler

"""
    kepler_to_rv(a::T1, e::T2, i::T3, Ω::T4, ω::T5, f::T6) where {T1, T2, T3, T4, T5, T6}
    kepler_to_rv(k::KeplerianElements)

Convert the Keplerian elements (`a`, `e`, `i`, `Ω`, `ω`, and `f`) to a Cartesian
representation (position vector `r` and velocity vector `v`). The Keplerian
elements can also be passed inside an instance of the `KeplerianElements`
structure.

# Args

- `a`: Semi-major axis [m].
- `e`: Eccentricity.
- `i`: Inclination [rad].
- `Ω`: Right ascension of the ascending node [rad].
- `ω`: Argument of perigee [rad].
- `f`: True anomaly [rad].

# Returns

- The position vector represented in the inertial reference frame [m].
- The velocity vector represented in the inertial reference frame [m].

# Remarks

This algorithm was adapted from **[1]** and **[2]**(p. 37-38).

# References

- **[1]**: Schwarz, R (2014). Memorandum No. 2: Cartesian State Vectors to
    Keplerian Orbit Elements. Available at www.rene-schwarz.com
- **[2]**: Kuga, H. K., Carrara, V., Rao, K. R (2005). Introdução à Mecânica
    Orbital. 2ª ed. Instituto Nacional de Pesquisas Espaciais.
"""
function kepler_to_rv(
    a::T1,
    e::T2,
    i::T3,
    Ω::T4,
    ω::T5,
    f::T6
) where {T1, T2, T3, T4, T5, T6}
    T = promote_type(T1, T2, T3, T4, T5, T6)

    # Check eccentricity.
    !(0 <= e < 1) && throw(ArgumentError("Eccentricity must be in the interval [0,1)."))

    # Auxiliary variables.
    sin_f, cos_f = sincos(f)

    # Compute the geocentric distance.
    r = a * (1 - e^2) / (1 + e * cos_f)

    # Compute the position vector in the orbit plane, defined as:
    #   - The X axis points towards the perigee;
    #   - The Z axis is perpendicular to the orbital plane (right-hand);
    #   - The Y axis completes a right-hand coordinate system.
    r_o = SVector{3, T}(r * cos_f, r * sin_f, 0)

    # Compute the velocity vector in the orbit plane without perturbations.
    n = angvel(a, e, i, :J0)
    v_o = n * a / sqrt(1 - e^2) * SVector{3, T}(-sin_f, e + cos_f, 0)

    # Compute the matrix that rotates the orbit reference frame into the
    # inertial reference frame.
    Dio = angle_to_dcm(-ω, -i, -Ω, :ZXZ)

    # Compute the position and velocity represented in the inertial frame.
    r_i = Dio * r_o
    v_i = Dio * v_o

    return r_i, v_i
end

kepler_to_rv(k::KeplerianElements) = kepler_to_rv(k.a, k.e, k.i, k.Ω, k.ω, k.f)

"""
    rv_to_kepler(r_i::AbstractVector, v_i::AbstractVector, t::Number = 0)
    rv_to_kepler(x::Number, y::Number, z::Number, vx::Number, vy::Number, vz::Number, t::Number = 0)

Convert a Cartesian representation (position vector `r_i` [m] and velocity
vector `v_i` [m/s²]) to the Keplerian elements. Optionally, the user can specify
the epoch of the returned elements using the parameter `t`. It it is omitted,
then it default to 0.

The input vectors can also be passed component by component:

    r_i = [x,   y,  z]
    v_i = [vx, vy, vz]

# Returns

An instance of the structure `KeplerianElements` [SI units].

# Remarks

The algorithm was adapted from **[1]**.

The special cases are treated as follows:

- **Circular and equatorial**: the right ascension of the ascending node and the
    argument of perigee are set to 0. Hence, the true anomaly is equal to the
    true longitude.
- **Elliptical and equatorial**: the right ascension of the ascending node is
    set to 0. Hence, the argument of perigee is equal to the longitude of
    periapsis.
- **Circular and inclined**: the argument of perigee is set to 0. Hence, the
    true anomaly is equal to the argument of latitude.

# References

- **[1]**: Schwarz, R (2014). Memorandum No. 2: Cartesian State Vectors to
    Keplerian Orbit Elements. Available at www.rene-schwarz.com
"""
function rv_to_kepler(
    r_i::AbstractVector{T1},
    v_i::AbstractVector{T2},
    t::Number = 0
) where {T1, T2}
    # Check inputs.
    length(r_i) != 3 && error("The vector r_i must have 3 dimensions.")
    length(v_i) != 3 && error("The vector v_i must have 3 dimensions.")

    # Obtain the type of the output elements.
    T = float(promote_type(T1, T2))

    @inbounds begin
        # Position and velocity vector norms and auxiliary dot products.
        r2 = T(dot(r_i, r_i))
        v2 = T(dot(v_i, v_i))

        r  = sqrt(r2)
        v  = sqrt(v2)

        rv = T(dot(r_i, v_i))

        Tm0 = T(m0)

        # Angular momentum vector.
        h_i = r_i × v_i
        h   = norm(h_i)

        # Vector that points to the right ascension of the ascending node (RAAN).
        n_i = SVector{3}(0, 0, 1) × h_i
        n   = norm(n_i)

        # Eccentricity vector.
        e_i = ((v2 - Tm0 / r) * r_i - rv * v_i ) / Tm0

        # Orbit energy.
        ξ = v2 / 2 - Tm0 / r

        # Eccentricity
        # ============

        ecc = norm(e_i)

        # Semi-major axis
        # ===============

        if abs(ecc) <= 1.0 - 1e-6
            a = -Tm0 / (2ξ)
        else
            error("Could not convert the provided Cartesian values to Kepler elements.\n" *
                  "The computed eccentricity was not between 0 and 1");
        end

        # Inclination
        # ===========

        cos_i = h_i[3] / h
        cos_i = abs(cos_i) > 1 ? sign(cos_i) : cos_i
        i     = acos(cos_i)

        # Check the type of the orbit to account for special cases
        # ======================================================================

        # Equatorial
        # ----------------------------------------------------------------------

        if abs(n) <= 1e-6

            # Right Ascension of the Ascending Node.
            # ======================================

            Ω = T(0)

            # Equatorial and elliptical
            # ------------------------------------------------------------------

            if abs(ecc) > 1e-6

                # Argument of Perigee
                # ===================

                cos_ω = e_i[1] / ecc
                cos_ω = abs(cos_ω) > 1 ? sign(cos_ω) : cos_ω
                ω     = acos(cos_ω)

                (e_i[3] < 0) && (ω = T(2π) - ω)

                # True anomaly
                # ============

                cos_v = dot(e_i, r_i) / (ecc * r)
                cos_v = abs(cos_v) > 1 ? sign(cos_v) : cos_v
                v     = acos(cos_v)

                (rv < 0) && (v = T(2π) - v)

            # Equatorial and circular
            # ------------------------------------------------------------------

            else
                # Argument of Perigee
                # ===================

                ω = T(0)

                # True anomaly
                # ============

                cos_v = r_i[1] / r
                cos_v = abs(cos_v) > 1 ? sign(cos_v) : cos_v
                v     = acos(cos_v)

                (r_i[2] < 0) && (v = T(2π) - v)
            end

        # Inclined
        # ----------------------------------------------------------------------
        else

            # Right Ascension of the Ascending Node.
            # ======================================

            cos_Ω = n_i[1] / n
            cos_Ω = abs(cos_Ω) > 1 ? sign(cos_Ω) : cos_Ω
            Ω     = acos(cos_Ω)

            (n_i[2] < 0) && (Ω = T(2π) - Ω)

            # Circular and inclined
            # ------------------------------------------------------------------

            if abs(ecc) < 1e-6

                # Argument of Perigee
                # ===================

                ω = T(0)

                # True anomaly
                # ============

                cos_v = dot(n_i, r_i) / (n*r)
                cos_v = abs(cos_v) > 1 ? sign(cos_v) : cos_v
                v     = acos(cos_v)

                (r_i[3] < 0) && (v = T(2π) - v)
            else

                # Argument of Perigee
                # ===================

                cos_ω = dot(n_i, e_i) / (n * ecc)
                cos_ω = abs(cos_ω) > 1 ? sign(cos_ω) : cos_ω
                ω     = acos(cos_ω)

                (e_i[3] < 0) && (ω = T(2π) - ω)

                # True anomaly
                # ============

                cos_v = dot(e_i, r_i) / (ecc * r)
                cos_v = abs(cos_v) > 1 ? sign(cos_v) : cos_v
                v     = acos(cos_v)

                (rv < 0) && (v = T(2π) - v)
            end
        end
    end

    # Return the Keplerian elements.
    # ==============================

    return KeplerianElements(t, a, ecc, i, Ω, ω, v)
end

function rv_to_kepler(x::Number,  y::Number,  z::Number,
                      vx::Number, vy::Number, vz::Number, t::Number = 0)
    # Create the position and velocity vectors.
    r_i = SVector{3}( x, y, z)
    v_i = SVector{3}(vx,vy,vz)

    # Compute the Keplerian orbit elements.
    return rv_to_kepler(r_i, v_i, t)
end


text:
#############################################################################
# Joulia
# A Large-Scale Spatial Power System Model for Julia
# See https://github.com/JuliaEnergy/Joulia.jl
#############################################################################
# This file contains functions and structs for lines

zbase(voltage::Number) = (voltage *1E3)^2 / (500 * 1E6)

"""
    calc_pmax(df::DataFrame)

Calculate thermal limit from line data including TRM.
"""
function calc_pmax(df::DataFrame)
	L = df[1]
	voltage = map(zbase, df[:Voltage])
	reactance = Dict(zip(L, df[:Reactance] ./ df[:Circuits] ./ voltage))
	resistance = Dict(zip(L,  df[:Resistance] ./ df[:Circuits] ./ voltage))

	pmax = Dict(zip(L,df[:ThermalLimit] .* df[:Circuits] .* (1-(20/100)) ))
	return reactance, resistance, pmax
end

mutable struct Lines
    id
    from
    to
	reactance
	resistance
	pmax

    function Lines(df::DataFrame)
		df[:id] = Symbol.(df[:id])
		reactance, resistance, pmax = calc_pmax(df)
        dict = df_to_dict_with_id(df)

        self =  new(df[:id],
                   dict[:From],
                   dict[:To],
				   reactance,
				   resistance,
                   pmax)

        return self
    end
end

"""
    create_incidence(lines::Lines, nodes::Nodes)

Create the incidence matrix from 'lines' and 'nodes'.
"""
function create_incidence(lines::Lines, nodes::Nodes)
	incidence = zeros(Int, length(lines.id), length(nodes.id))
	for (i,l) in enumerate(lines.id)
	    incidence[i, findfirst(nodes.id .== Symbol(lines.from[l]))] = 1
	    incidence[i, findfirst(nodes.id .== Symbol(lines.to[l]))] = -1
	end
	return incidence
end

"""
    calc_h_b(lines::Lines, nodes::Nodes)

Calculate the H and B matrices from 'lines' and 'nodes'.
"""
function calc_h_b(lines::Lines, nodes::Nodes)

	incidence = create_incidence(lines, nodes)
	bvector = [lines.reactance[l] /
		(lines.reactance[l].^2 .+ lines.resistance[l].^2) for l in lines.id]
	h_matrix = bvector .* incidence
	b_matrix = h_matrix' *incidence

	h = Dict((l,n) => h_matrix[i,j]
		for (i,l) in enumerate(lines.id), (j,n) in enumerate(nodes.id))

	b = Dict((n,nn) => b_matrix[i,j]
		for (i,n) in enumerate(nodes.id), (j,nn) in enumerate(nodes.id))

	return h,b
end


text:
using Labels
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
label = Label()

label.var[:race] = "Race/ethnicity"
label.val[:race] = Dict(1=>"White",2=>"Black",3=>"Hispanic",4=>"Other")
label.lblname[:race] = :race

@test varlab(label,:race) == "Race/ethnicity"
@test vallab(label,:race,1) == "White"


text:
using Test

include("one_euro_filter.jl")
no_errors(x) = true

@test no_errors(smoothing_factor(1.0, 1.0))
@test no_errors(exponential_smoothing(0.5, 1.0, 1.0))
@test no_errors(one_euro_filter(0.0, 0.0, 1.0, 1.0))
@test no_errors(one_euro_filter(0.0, 0.0, 1.0, 1.0, 1.0))

@test no_errors(smoothing_factor.(1.0, [1.0, 2.0]))
@test no_errors(exponential_smoothing.([0.5, 0.6], [1.0, 2.0], [1.0, 2.0]))
@test no_errors(one_euro_filter(0.0, [0.0], 1.0, [1.0]))
@test no_errors(one_euro_filter(0.0, [0.0], 1.0, [1.0], [1.0]))

max_stars_repo_name: jaantollander/OneEuroFilter

max_issues_repo_name: jaantollander/OneEuroFilter

max_forks_repo_name: jaantollander/OneEuroFilter


text:
using .util_convexopt

is_initiated = false

function __init__()
    is_initiated = true
    py"""
    # reference: https://github.com/pathak22/pyflow.git
    import numpy as np
    import pyflow

    #find flow from img2 to img1
    def py_flow(img1,img2,alpha=0.5, ratio=0.75, minWidth=20, nOuterFPIterations=7, nInnerFPIterations=1, nSORIterations=30, colType=1):
        img1 = np.array(img1)
        img2 = np.array(img2)
        height, width = img1.shape

        im1 = np.array(img1).reshape(height,width,1).copy(order='C')
        im2 = np.array(img2).reshape(height,width,1).copy(order='C')

        u, v, im2W = pyflow.coarse2fine_flow(
            im1, im2, alpha, ratio, minWidth, nOuterFPIterations, nInnerFPIterations,
            nSORIterations, colType)

        return u, v, im2W
    """
end

function get_flows(u::Array{T,3}, alpha=0.5, ratio=0.75, minWidth=20, nOuterFPIterations=7, nInnerFPIterations=1, nSORIterations=30, colType=1)::Array{T,4} where {T<:AbstractFloat}
    #ref: https://github.com/pathak22/pyflow.git
    @suppress begin
        #HACK!
        if !is_initiated
            __init__()
        end
        height, width, frames = size(u)
        result = zeros(height, width, 2, frames)
        for i in 1:frames-1
            flow_x, flow_y, im1Warped = py"py_flow"(u[:,:,i+1],u[:,:,i], alpha=alpha, ratio=ratio, minWidth=minWidth, nOuterFPIterations=nOuterFPIterations, nInnerFPIterations=nInnerFPIterations, nSORIterations=nSORIterations, colType=colType)
            result[:,:,1,i] = flow_x
            result[:,:,2,i] = flow_y
        end
        return result
    end
end

function sub2ind(array_shape, rows, cols)
    out = (cols.-1).*array_shape[1] .+ rows
    return convert.(Int64, out)
end

function compute_warping_operator(flow)
    # ref: https://github.com/HendrikMuenster/flexBox/blob/79abc7285703911cca2653434cca5bcefc79c722/operators/generateMatrix/warpingOperator.m
    """
    u: gray image
    v: flow [H, W, 2]
    """
    H, W, dims = size(flow)
    shape = (H,W)
    XX, YY = repeat(collect(1:W)', H, 1), repeat(collect(1:H), 1, W)

    targetPoint2Mat = XX + flow[:,:,1]
    targetPoint1Mat = YY + flow[:,:,2]

    x1 = floor.(targetPoint1Mat).-1
    x2 = x1 .+ 1
    x3 = x1 .+ 2
    x4 = x1 .+ 3

    y1 = floor.(targetPoint2Mat).-1
    y2 = y1 .+ 1
    y3 = y1 .+ 2
    y4 = y1 .+ 3

    v2 = targetPoint1Mat .- x2
    v1 = targetPoint2Mat .- y2

    indicator = ones(H,W)
    indicator = (x1 .> 0).*indicator
    indicator = (x4 .<= W).*indicator
    indicator = (y1 .> 0).*indicator
    indicator = (y4 .<= H).*indicator

    indicator = indicator .> 0
    v1 = v1[indicator]
    v2 = v2[indicator]

    idxList = sub2ind(shape, YY, XX)
    idxList_ = collect(Iterators.flatten(idxList))

    #find the indices
    list_J  = vcat(
        idxList_[sub2ind(shape, x2[indicator], y2[indicator])],
        idxList_[sub2ind(shape, x3[indicator], y2[indicator])],
        idxList_[sub2ind(shape, x2[indicator], y3[indicator])],
        idxList_[sub2ind(shape, x3[indicator], y3[indicator])])

    list_J = collect(Iterators.flatten(list_J))

    list_val  = ( (1 .-v1) .* (1 .-v2))
    list_val = vcat(list_val,(v2 .* (1 .-v1)))
    list_val = vcat(list_val,(v1 .* (1 .-v2)))
    list_val = vcat(list_val,(v1 .* v2))

    list_val = collect(Iterators.flatten(list_val))

    list_I = idxList[indicator]
    list_I =collect(Iterators.flatten(repeat(list_I,1,4)))

    W_mat = dropzeros!(sparse(list_I, list_J, list_val, H*W, H*W))

    return LinearOperator(W_mat)
end

function _recon2d_tv_primaldual_flow(As,A_norm,∇_norm, W_list,u0s,bs,w_tv,w_flow,c,niter,p1,p2,p3, mask)
    #A variational reconstruction method for undersampled dynamic x-ray tomography based on physical motion models (Burger)#
    height, width, frames = size(u0s)

    u = u0s
    u_prev = similar(u)
    ubar = deepcopy(u)

    p_adjoint  = similar(p3)
    p3_adjoint = similar(p3)#zeros(height*width, frames)
    p_adjoint_prev = similar(p3)
    W_norm = maximum(map(W-> compute_opnorm(W), W_list))
    opsnorm = [A_norm, ∇_norm, W_norm]
    tau = c /sum(opsnorm)

    data_term = similar(bs)
    u_grad = similar(p2)
    u_temp = similar(u)
    divp2 = similar(p2)

    flow_term = similar(u)

    sigmas = map(n -> 1.0/(n*c), opsnorm)

    @info sigmas
    it = 1
    primal_gap = 100
    tol = 0.001
    while it < niter && primal_gap > tol
        u_prev .= u
        p_adjoint_prev .= p_adjoint
        for t=frames:-1:1
            A = As[t]
            #dual
            @views data_term[:,t] .= mul!(data_term[:,t], A, vec(ubar[:,:,t])) .- bs[:,t]
            @views p1[:,t] .= (p1[:,t] .+ sigmas[1] * data_term[:,t]) ./ (sigmas[1] + 1.0)# l2 norm
            @views mul!( view(p_adjoint[:,:,t], :), A', p1[:,t]) #projection?

            @views util_convexopt.grad!(u_grad[:,:,:,t], ubar[:,:,t])
            @views p2[:,:,:,t] .+= sigmas[2] .* u_grad[:,:,:,t]
            @views util_convexopt.proj_dual_iso!(p2[:,:,:,t], u_temp[:,:,t], w_tv)
            @views p_adjoint[:,:,t] .-= util_convexopt.div!(divp2[:,:,:,t], p2[:,:,:,t])

            if t < frames
               @views mul!(view(flow_term[:,:,t], :), W_list[t], vec(ubar[:,:,t]))
               @views flow_term[:,:,t].-=ubar[:,:,t+1]
               @views p3[:,:,t] .= (p3[:,:,t] .+ sigmas[3] * flow_term[:,:,t])
               @views util_convexopt.proj_dual_iso!(p3[:,:,t], w_flow)

               @views mul!(view(p3_adjoint[:,:,t], :), W_list[t]', vec(p3[:,:,t]))
               @views p_adjoint[:,:,t] .+= p3_adjoint[:,:,t]
            end
        end
        u .-= tau * p_adjoint
        # primal update
        u .= max.(u.-tau .* p_adjoint, 0.0) # positivity constraint
        u .*=mask

        # acceleration
        ubar .= 2 .* u .- u_prev

        if it % 50 == 0
            #du = u_prev - u
            primal_gap = sum(abs.(-p_adjoint+p_adjoint_prev + (u_prev-u)/tau)) / length(p_adjoint)
            #@info "primal gap:" primal_gap
            @info it primal_gap
        end
        it = it+1
    end
    return u
end

"""
    recon2d_tv_primaldual_flow(As, bs, u0s, niter, w_tv, w_flow, c=1.0)

Reconstruct a 2d image by total variation model using Primal Dual optimization method and optical flow estimation

# Args
As : List of forward opeartors
bs : List of projection data
u0s: Initial guess of images
niter: number of iterations
w_tv: weight for TV term
W_flow: weight for flow term
c : See 61 page in 2016_Chambolle,Pock_An_introduction_to_continuous_optimization_for_imagingActa_Numerica

"""
function recon2d_tv_primaldual_flow(A_list, bs::Array{R, 2}, u0s::Array{R, 3}, niter1::Int, niter2::Int, w_tv::Float64, w_flow::Float64, c, mask, alpha) where {R <: AbstractFloat}

    height,width,frames = size(u0s)
    v = zeros(height, width, 2, frames)
    u = u0s
    shape = (height*width, frames)
    A_norm = maximum(map(A->compute_opnorm(A), A_list))#compute_opnorm_block(A_list, shape)
    N = height #only ok for height=width
    Dfd = spdiagm(0 => -ones(N-1), 1=>ones(N-1))
    ∇ = kron(Dfd, sparse(I,N,N)) + kron(Dfd, sparse(I.*im, N, N))
    ∇_norm = abs(compute_opnorm(∇))

    p1 = zeros(size(bs))
    p2 = zeros(height, width, 2,frames)
    p3 = zeros(height, width, frames)
    tol = 0.001
    i = 1
    rmain = 100
    while i < niter1 && rmain > tol
        u_prev = u
        v_prev = v
        W_list = mapslices(f -> compute_warping_operator(f), v,dims=[1,2,3])

        #W_norm = map(W)=#compute_opnorm_block(W_list, shape)
        u = _recon2d_tv_primaldual_flow(A_list,A_norm,∇_norm, W_list,u,bs,w_tv,w_flow,c,niter2,p1,p2,p3, mask)

        v = get_flows(u, alpha)
        rmain = norm(u-u_prev)+norm(v-v_prev)
        i = i+1
        @info i rmain
    end
    return u
end

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
module Verlet

export VerletIntegrator
include("./types.jl")
include("./main.jl")

end


text:
using NonLinearReactionAdvectionDiffusionWithFrontData;
using NonLinearReactionAdvectionDiffusionWithFrontData: shishkin_mesh;
using LinearAlgebra;

u_l(t) = -8;                    # ГУ
u_r(t) =  4;                    #
qf(x) = 4*sin(3 * π * x);       # Коэффициент линейного усиления, который в обратной
# задаче необходимо определить, но при генерации априорной
# информации мы задаем некоторый коэффициент, который,
# собственно, после имея априорную информацию и будем определять.
ε = 0.01;                       # Малый параметр при старшей производной
a, b = 0, 1;                    # Область по X
t₀, T = 0, 0.50;                # Область по T
x_tp = 0.12;                    # Положение переходного слоя
M = 500;                        # Кол-во разбиений по X, T
τ = (T-t₀)/M;                   # шаг по T
Tₘ = [t₀ + m*τ for m in 0:M];   # Сетка по Т
ulₘ=    u_l.(Tₘ);               # Сеточные значения левого  ГУ
urₘ=    u_r.(Tₘ);               # Сеточные значения правого ГУ
nothing #hide

#' Новое поведение будет реализоваться с помощью передачи функции создания сетки внутрь `solve` в качестве аргумента.
#' Эта функция должна принимать только один аргумент — положение переходного слоя и формировать соответствующую сетку.
#' В пакете есть формирование кусочно-равномерной сетки со сгущениями на границе и на переходном слое.
#' Создадим замыкание этой функции, которое будет иметь нужную сигнатуру.
mshfrm(x_tp) = shishkin_mesh(a, b, x_tp, ε, 40, 0.5, 1.0, 1.0, 0.25);
Xₙ = mshfrm(x_tp); ;                            # Сетка по Х
N = length(Xₙ) - 1                              # Примем за N длину сетки, что получилась в итоге.
qₙ =    qf.(Xₙ);                                # Сеточные значения коэффициента линейного усиления
y₀ = u_init.(Xₙ, ε=ε, x_tp = x_tp);             # Начальные условия

# Подготовим массивы, выбросив граничные точки
# ведь тестируемая функция — для внутреннего использования
qq = qₙ[2:N];
y = y₀[2:N];

dl, d, du = NonLinearReactionAdvectionDiffusionWithFrontData.DRP_y(y, 1, Xₙ, N, ε, ulₘ, urₘ, qq)

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

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

tdjac = NonLinearReactionAdvectionDiffusionWithFrontData.∂DRP_∂y(y, 1, Xₙ, N, ε, ulₘ, urₘ, qq)
@test typeof(tdjac) <: Tridiagonal
@test tdjac == Tridiagonal( dl, d, du )

# Сравним с якобианом автодифференцирования
jac = NonLinearReactionAdvectionDiffusionWithFrontData.∂directRP_∂y(y, 1, Xₙ, N, ε, ulₘ, urₘ, qq)
@test isapprox(tdjac, jac)


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

using Markdown
using InteractiveUtils

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

# ╔═╡ be4635d0-9bc1-11eb-2efd-d500550f4620
begin
	using PlutoUI 	# use interactive HTML binding of variables
	using Plots 	# use Plot package to ... plot
	using Printf
	plotlyjs()		# use interactive javascript backend for Plots.jl
	# gr()				# use the more conventional GR backend for plotting
end

# ╔═╡ 955fe07b-dd41-447a-98d5-1866904eefc3
function CEven(n::Int,ε) # define even coefficients
	if n==0
		return 1.0
	elseif n == 1
		return 0.0
	end
	return CEven(n-2,ε) * (2.0*(n-2) + 1 - 2ε)/(n*(n-1))
end

# ╔═╡ a43cd1de-02d0-45e2-8610-dae931ab8c10
function COdd(n::Int,ε) # define odd coefficients
	if n==1
		return 1.0
	elseif n == 0
		return 0.0
	end
	return COdd(n-2,ε) * (2.0*(n-2) + 1 - 2ε)/(n*(n-1))
end

# ╔═╡ 460a5f0a-d776-4bfd-8e67-4cab0315fd85
numberOfTerms = 200

# ╔═╡ 66acd7ad-45a9-4af7-b940-4b7d2ca9691c
ψEven(x, ε) = sum( η -> CEven(η, ε) * x^η, 0:1:numberOfTerms) * exp(-x^2/2)

# ╔═╡ 3f79e204-a2a6-46a2-b910-a3ce919cf340
ψOdd(x, ε) = sum( η -> COdd(η, ε) * x^η, 0:1:numberOfTerms) * exp(-x^2/2)

# ╔═╡ f4eca874-f58f-422e-a31b-114bde7373cd
xPlot = collect(-5:0.01:5)

# ╔═╡ 4ec87017-4f8e-4d1e-b1e0-34d3f6ca62aa
@bind ε Slider(0.3:0.1:7.6; default=0.3, show_value=true)

# ╔═╡ c39942fc-dfad-4afe-abed-5f523c265caa
begin
	plotlyjs()
	plot1 = plot(xPlot, xPlot .^2, label="V(x) = x²")
		ylims!(plot1, (-4,4))
		title!(plot1, "Even wavefunctions")
		plot!(plot1, xPlot, ψEven.(xPlot, ε), label="ψ even, ε=$(@sprintf("%.2f",ε))")

end

# ╔═╡ 335d3458-a511-42bb-af11-ca18b11432c0
begin
	plotlyjs()
	plot2 = plot(xPlot, xPlot .^2, label="V(x) = x²")
	plot!(plot2, xPlot, ψOdd.(xPlot, ε), label="ψ odd, ε=$(@sprintf("%.2f",ε))")
	ylims!(plot2, (-4,4))
	title!(plot2, "Odd wavefunctions")
end

# ╔═╡ 98216100-19a0-46ca-802a-7cfe46e4e138
begin
	gr()
	@gif for ε ∈ 0.1:0.05:10.0
		plot(xPlot, xPlot .^2, label="V(x) = x²")
		ylims!((-4,4))
		title!("Even wavefunctions")
		plot!(xPlot, ψEven.(xPlot, ε), label="ψ even, ε=$(@sprintf("%.2f",ε))")
	end
end

# ╔═╡ b5fcc994-c771-46dd-9558-248c58ce08df
begin
	gr()
	@gif for ε ∈ 0.1:0.05:10.0
		plot(xPlot, xPlot .^2, label="V(x) = x²")
		ylims!((-4,4))
		title!("Odd wavefunctions")
		plot!(xPlot, ψOdd.(xPlot, ε), label="ψ odd, ε=$(@sprintf("%.2f",ε))")
	end
end

# ╔═╡ Cell order:
# ╠═be4635d0-9bc1-11eb-2efd-d500550f4620
# ╠═955fe07b-dd41-447a-98d5-1866904eefc3
# ╠═a43cd1de-02d0-45e2-8610-dae931ab8c10
# ╠═460a5f0a-d776-4bfd-8e67-4cab0315fd85
# ╠═66acd7ad-45a9-4af7-b940-4b7d2ca9691c
# ╠═3f79e204-a2a6-46a2-b910-a3ce919cf340
# ╠═f4eca874-f58f-422e-a31b-114bde7373cd
# ╠═4ec87017-4f8e-4d1e-b1e0-34d3f6ca62aa
# ╠═c39942fc-dfad-4afe-abed-5f523c265caa
# ╠═335d3458-a511-42bb-af11-ca18b11432c0
# ╠═98216100-19a0-46ca-802a-7cfe46e4e138
# ╠═b5fcc994-c771-46dd-9558-248c58ce08df

max_stars_repo_name: jtschneider/quantum-harmonic-oscillator

max_issues_repo_name: jtschneider/quantum-harmonic-oscillator

max_forks_repo_name: jtschneider/quantum-harmonic-oscillator