text:
using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.DGmethods
using CLIMA.DGmethods.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using CLIMA.VariableTemplates: flattenednames, Vars
using CLIMA.HydrostaticBoussinesq
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates
using CLIMA.VTK
using CLIMA.PlanetParameters: grav
using CLIMA.GeneralizedMinimalResidualSolver
using CLIMA.ColumnwiseLUSolver: ManyColumnLU, SingleColumnLU
using CLIMA.HydrostaticBoussinesq: AbstractHydrostaticBoussinesqProblem
import CLIMA.HydrostaticBoussinesq: ocean_init_aux!, ocean_init_state!,
                                    ocean_boundary_state!,
                                    CoastlineFreeSlip, CoastlineNoSlip,
                                    OceanFloorFreeSlip, OceanFloorNoSlip,
                                    OceanSurfaceNoStressNoForcing,
                                    OceanSurfaceStressNoForcing,
                                    OceanSurfaceNoStressForcing,
                                    OceanSurfaceStressForcing
import CLIMA.DGmethods: update_aux!, vars_state, vars_aux, VerticalDirection
using Test
using GPUifyLoops

const ArrayType = CLIMA.array_type()

HBModel   = HydrostaticBoussinesqModel

struct HomogeneousSimpleBox{T} <: AbstractHydrostaticBoussinesqProblem
  Lˣ::T
  Lʸ::T
  H::T
  τₒ::T
  fₒ::T
  β::T
end

HSBox = HomogeneousSimpleBox

@inline function ocean_boundary_state!(m::HBModel, p::HSBox, bctype, x...)
  if bctype == 1
    ocean_boundary_state!(m, CoastlineNoSlip(), x...)
  elseif bctype == 2
    ocean_boundary_state!(m, OceanFloorFreeSlip(), x...)
  elseif bctype == 3
    ocean_boundary_state!(m, OceanSurfaceStressNoForcing(), x...)
  end
end

# aux is Filled afer the state
function ocean_init_aux!(m::HBModel, P::HSBox, A, geom)
  FT = eltype(A)
  @inbounds y = geom.coord[2]

  Lʸ = P.Lʸ
  τₒ = P.τₒ
  fₒ = P.fₒ
  β  = P.β

  A.τ  = -τₒ * cos(y * π / Lʸ)
  A.f  =  fₒ + β * y

  A.ν = @SVector [m.νʰ, m.νʰ, m.νᶻ]
  A.κ = @SVector [m.κʰ, m.κʰ, m.κᶻ]
end

function ocean_init_state!(p::HSBox, state, aux, coords, t)
  @inbounds z = coords[3]
  @inbounds H = p.H

  state.u = @SVector [rand(),rand()]
  state.η = 0
  state.θ = 20
end

function main()
  CLIMA.init()
  mpicomm = MPI.COMM_WORLD

  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
             ll == "WARN"  ? Logging.Warn  :
             ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))

  brickrange = (xrange, yrange, zrange)
  topl = StackedBrickTopology(mpicomm, brickrange;
                              periodicity = (false, false, false),
                              boundary = ((1, 1), (1, 1), (2, 3)))
  dt = 50
  nout = ceil(Int64, tout / dt)
  dt = tout / nout

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                         )

  prob = HSBox{FT}(Lˣ, Lʸ, H, τₒ, fₒ, β)

  model = HBModel{FT}(prob, cʰ = cʰ)

  linearmodel = LinearHBModel(model)

  dg = OceanDGModel(model,
                    grid,
                    Rusanov(),
                    CentralNumericalFluxDiffusive(),
                    CentralNumericalFluxGradient())

  lineardg = DGModel(linearmodel,
                     grid,
                     Rusanov(),
                     CentralNumericalFluxDiffusive(),
                     CentralNumericalFluxGradient();
                     direction=VerticalDirection(),
                     auxstate=dg.auxstate)

  Q = init_ode_state(dg, FT(0); forcecpu=true)
  update_aux!(dg, model, Q, FT(0))

  linearsolver = SingleColumnLU() # ManyColumnLU()

  odesolver = ARK2GiraldoKellyConstantinescu(dg, lineardg, linearsolver, Q;
                                             dt = dt, t0 = 0,
                                             split_nonlinear_linear=false)

  cbvector = make_callbacks(step, nout, mpicomm, odesolver, dg, model, Q)

  eng0 = norm(Q)
  @info @sprintf """Starting
  norm(Q₀) = %.16e
  ArrayType = %s""" eng0 ArrayType

  solve!(Q, odesolver, nothing; timeend=timeend, callbacks=cbvector, adjustfinalstep=false)

  maxQ =  Vars{vars_state(model, FT)}(maximum(Q, dims=(1,3)))
  minQ =  Vars{vars_state(model, FT)}(minimum(Q, dims=(1,3)))

  @test maxQ.θ ≈ minQ.θ

  return nothing
end

function make_callbacks(step, nout, mpicomm, odesolver, dg, model, Q)
  starttime = Ref(now())
  cbinfo = GenericCallbacks.EveryXWallTimeSeconds(60, mpicomm) do (s=false)
    if s
      starttime[] = now()
    else
      energy = norm(Q)
      @info @sprintf("""Update
                     simtime = %.16e
                     runtime = %s
                     norm(Q) = %.16e""", ODESolvers.gettime(odesolver),
                     Dates.format(convert(Dates.DateTime,
                                          Dates.now()-starttime[]),
                                  Dates.dateformat"HH:MM:SS"),
                     energy)
    end
  end

  return (cbinfo,)
end

#################
# RUN THE TESTS #
#################
FT = Float64

const timeend = 3600 # s
const tout    = 60 * 60 # s

const N  = 4
const Nˣ = 20
const Nʸ = 20
const Nᶻ = 50
const Lˣ = 4e6   # m
const Lʸ = 4e6   # m
const H  = 400   # m

xrange = range(FT(0);  length=Nˣ+1, stop=Lˣ)
yrange = range(FT(0);  length=Nʸ+1, stop=Lʸ)
zrange = range(FT(-H); length=Nᶻ+1, stop=0)

const cʰ = sqrt(grav * H)
const cᶻ = 0

const τₒ = 1e-1  # (m/s)^2
const fₒ = 1e-4  # Hz
const β  = 1e-11 # Hz / m

const αᵀ = 2e-4  # (m/s)^2 / K
const νʰ = 5e3   # m^2 / s
const νᶻ = 5e-3  # m^2 / s
const κʰ = 1e3   # m^2 / s
const κᶻ = 1e-10  # m^2 / s

main()

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
"""
    Outpoint

Because a single transaction can include multiple outputs,
the Outpoint structure includes both a TXID (or `txid`) and an output `index`
number to refer to specific output.
- `txid`, the TXID of the transaction holding the output to spend.
  The TXID is a hash provided here in internal byte order.
- `index`, the output index number of the specific output to spend from the transaction.
  The first output is 0x00000000.
"""
struct Outpoint
    txid    :: Vector{UInt8}
    index   :: UInt32
end

"""
    Outpoint(io::IOBuffer) -> Outpoint

Parse an `Outpoint` from an IOBuffer
"""
function Outpoint(io::IOBuffer)
    txid = reverse(read(io, 32))
    index = read(io, UInt32)
    Outpoint(txid, index)
end

function Base.show(io::IO, prevout::Outpoint)
    println(io, "Outpoint: " * bytes2hex(prevout.txid) * ":" * string(prevout.index, base=10))
end

"""
    serialize(prevout::Outpoint) -> Vector{UInt8}

Returns the byte serialization of the outpoint
"""
function serialize(prevout::Outpoint)
    result = copy(prevout.txid)
    reverse!(result)
    append!(result, bytes(prevout.index, len=4, little_endian=true))
    return result
end

ext: jl


text:
module GeoFormatTypes

# Use the README as the module docs
@doc let
    path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    read(path, String)
end GeoFormatTypes

export GeoFormat

export CoordinateReferenceSystemFormat, EPSG, ProjString, CoordSys

export GeometryFormat, GeoJSON, KML

export MixedFormat, GML

export AbstractWellKnownText, WellKnownText, WellKnownText2, ESRIWellKnownText, WellKnownBinary

const PROJ_PREFIX = "+proj="
const EPSG_PREFIX = "EPSG:"
# TODO more verification that types are wrapping the right format.

# Traits for mixed crs/geometry formats.
abstract type FormatMode end
struct Geom <: FormatMode end
struct CRS <: FormatMode end
struct Mixed <: FormatMode end

"""
    val(f::GeoFormat)

Get the contained value of a GeoFormat type.
"""
function val end

"""
Abstract supertype for geospatial data formats
"""
abstract type GeoFormat end

# Convert from the same type does nothing.
Base.convert(::Type{T1}, source::T2) where {T1<:GeoFormat,T2<:T1} = source
# Convert uses the `mode` trait to distinguish crs form geometry conversion
Base.convert(target::Type{T1}, source::T2; kwargs...) where {T1<:GeoFormat,T2<:GeoFormat} = begin
    sourcemode = mode(source)
    targetmode = mode(target)
    convertmode = if targetmode isa Geom
        if sourcemode isa Union{Mixed,Geom}
            Geom() # Geom is the default if both formats are mixed
        else
            throw(ArgumentError("cannot convert $(typeof(source)) to $target"))
        end
    elseif targetmode isa CRS
        if sourcemode isa Union{Mixed,CRS}
            CRS()
        else
            throw(ArgumentError("cannot convert $(typeof(source)) to $target"))
        end
    else # targetmode isa Mixed
        # Mixed to Mixed defaults to Geom
        if sourcemode isa Union{Mixed,Geom}
            Geom()
        else
            CRS()
        end
    end
    convert(target, convertmode, source; kwargs...)
end

"""
Formats representing coordinate reference systems
"""
abstract type CoordinateReferenceSystemFormat <: GeoFormat end

"""
Formats representing geometries. These wrappers simply mark string
formats that may optionally be converted to Geoetry objects at a later point.
"""
abstract type GeometryFormat <: GeoFormat end

"""
Formats that may hold either or both coordinate reference systems and geometries.
"""
abstract type MixedFormat{X} <: GeoFormat end

val(x::GeoFormat) = x.val

mode(format::GeoFormat) = mode(typeof(format))
mode(::Type{<:GeometryFormat}) = Geom()
mode(::Type{<:CoordinateReferenceSystemFormat}) = CRS()
mode(::Type{<:MixedFormat}) = Mixed()
mode(::Type{<:MixedFormat{M}}) where M = M()

# Most GeoFormat types wrap String or have a constructor for string inputs
Base.convert(::Type{String}, input::GeoFormat) = val(input)
Base.convert(::Type{T}, input::AbstractString) where T <: GeoFormat = T(input)

"""
Proj string
"""
struct ProjString <: CoordinateReferenceSystemFormat
    val::String
    ProjString(input::String) = begin
        startswith(input, PROJ_PREFIX) ||
            throw(ArgumentError("Not a Proj string: $input does not start with $PROJ_PREFIX"))
        new(input)
    end
end

"""
Mapinfo CoordSys string
"""
struct CoordSys <: CoordinateReferenceSystemFormat
    val::String
end

"""
Well known text has a number of versions and standards, and can
represent coordinate reference systems or geometric data.
"""
abstract type AbstractWellKnownText{X} <: MixedFormat{X} end

"""
Well known text v1 following the OGC standard
"""
struct WellKnownText{X,T<:String} <: AbstractWellKnownText{X}
    mode::X
    val::T
end
WellKnownText(val) = WellKnownText(Mixed(), val)

"""
Well known text v2 following the new OGC standard
"""
struct WellKnownText2{X,T<:String} <: AbstractWellKnownText{X}
    mode::X
    val::T
end
WellKnownText2(val) = WellKnownText(Mixed(), val)

"""
Well known text following the ESRI standard
"""
struct ESRIWellKnownText{X,T<:String} <: AbstractWellKnownText{X}
    mode::X
    val::T
end

"""
Well known binary
"""
struct WellKnownBinary{X,T} <: MixedFormat{X}
    mode::X
    val::T
end
WellKnownBinary(val) = WellKnownBinary(Mixed(), val)

Base.convert(::Type{String}, input::WellKnownBinary) =
    error("`convert` to `String` is not defined for `WellKnownBinary`")

"""
EPSG code representing a coordinate reference system from the
EPSG spatial reference system registry.
"""
struct EPSG <: CoordinateReferenceSystemFormat
    val::Int
end

"""
Constructor for "EPSG:1234" string input
"""
EPSG(input::AbstractString) = begin
    startswith(input, EPSG_PREFIX) || throw(ArgumentError("String $input does no start with $EPSG_PREFIX"))
    code = parse(Int, input[findlast(EPSG_PREFIX, input).stop+1:end])
    EPSG(code)
end

Base.convert(::Type{Int}, input::EPSG) = val(input)
Base.convert(::Type{String}, input::EPSG) = string(EPSG_PREFIX, val(input))
Base.convert(::Type{EPSG}, input::Int) = EPSG(input)

"""
Keyhole Markup Language
"""
struct KML <: GeometryFormat
    val::String
end

"""
Geography Markup Language
"""
struct GML{X,T<:String} <: MixedFormat{X}
    mode::X
    val::T
end
GML(val) = GML(Mixed(), val)

"""
GeoJSON String or Dict
"""
struct GeoJSON{T} <: GeometryFormat
    val::T
end

end # module


text:
"""
The base type used by all terminals within a Koza tree.
"""
abstract KozaTerminal <: KozaNode

function compose_terminal(inputs::Vector{KozaInput}, def::AbstractString)
  label = def[1:Base.search(def, ':')-1]
  ty = def[Base.last(Base.search(def, "::"))+1:end]

  # Compose the unique type for this non-terminal.
  src = "type <: KozaTerminal
    parent::KozaParent
    constructor() = new()
    constructor(p::KozaParent) = new(p)
  end"
  t = anonymous_type(koza, src)

  # Compose the label function.
  eval(koza, Base.parse("label(::$(t)) = $(label)"))

  # Compose the "fresh" function.
  eval(koza, Base.parse("fresh(::$(t), p::KozaParent) = $(t)(p)"))

  # Compose the execution function.
  src = join(vcat(["::$(t)"], ["$(i.label)::$(i.ty)" for i in inputs]), ",")
  src = "execute($(src)) = $(label)"
  eval(koza, Base.parse(src))
  return t
end

"""
Produces a clone of a provided terminal node from a Koza tree.
"""
clone{T <: KozaTerminal}(n::T, p::KozaParent) = T(p)
isterminal(::KozaTerminal) = true
num_nodes(n::KozaTerminal) = 1
height(n::KozaTerminal) = 1
describe(n::KozaTerminal) = label(n)


text:
using InterpretableValidation
using Test
using LinearAlgebra
using Statistics
using Random

val = mvrandn(Random.GLOBAL_RNG, [0.5,0.5,0.5], [1,1,1,], 0.2*Array{Float64}(I, 3, 3), 100)
@test all(val .> 0.5) && all(val .<1.0)

val2 = mvrandn_μ(Random.GLOBAL_RNG,[0.75, 0.75, 0.75],  [0.5,0.5,0.5], [1,1,1,], 0.2*Array{Float64}(I, 3, 3), 100)

@test mean(val) < mean(val2)
@test abs(mean(val2) - 0.75) < 0.1


text:
include("HapiTests.jl")

using Hapi
using .HapiTests

Hapi.runtests()
HapiTests.runtests()


text:
using AntennaExpts, YagiAntProblem
using ExprSearch, DerivTreeVis
using ExprSearch.GP
using RLESUtils
import RLESTypes.SymbolTable
using NECPP

srand(1)
function test(N)
    problem = YagiAnt()
    grammar = get_grammar(problem)
    mda = min_depth_actions(grammar)
    for i = 1:N
        ind = rand(grammar, mda, 10)
        #derivtreevis(ind.derivtree, "treevis")
        expr = get_expr(ind.derivtree)
        f = open("test.nec","w")
        (gain, Z) = nec_run(problem, expr, f) do nec
            gain = nec_gain(nec, 0, 0, 0)
            imp_re = nec_impedance_real(nec, 0)
            imp_im = nec_impedance_imag(nec, 0)
            Z = Complex(imp_re, imp_im)
            (gain, Z)
        end
        close(f)
        Z0 = Complex(50.0, 0.0)
        v = vswr(Z, Z0)
        @show gain
        @show Z
        @show v
        v < 0.0 && break
    end
end


text:
## How to insert d3 into a window

"Load d3.js code into current window"
function inject_d3(w::Blink.AtomShell.Window)
    Blink.loadjs!(w, "https://cdnjs.cloudflare.com/ajax/libs/d3/3.5.16/d3.min.js")
end
inject_d3() = inject_d3(figure(gcf()))

## create an empty page for working with
## can't inject d3 into  blank window.
function empty_page(w::Blink.AtomShell.Window=figure(gcf()); kwargs...)
    f = tempname() * ".html"
    io = open(f, "w")
    tpl = Mustache.template_from_file(Pkg.dir("DThree", "tpl", "d3.html"))
    open(f, "w") do io
        Mustache.render(io, tpl; kwargs...)
    end
    loadfile(w, f)
end

"Pass javascript d3 calls to Blink window"
Blink.js(w::Blink.AtomShell.Window, q::D3) = js(w, JSString(q.render()))
Blink.js(q::D3) = js(figure(gcf()), q)

Blink.loadurl(url::AbstractString) = Blink.loadurl(figure(gcf()), url)
Blink.loadfile(fname::AbstractString) = Blink.loadfile(figure(gcf()), fname)

## ## Example
## w = Blink.Window()
## d3 = D3()
## loadurl(w, "http://www.julialang.org") # must!! load something first??
## inject_d3(w) # inject d3 library into page so that it can be manipulated
## js(w, d3.select("body").transition(2000).style("background-color", "black"))

## Example
## d3 = D3()
## w = figure()
## empty_page(w)  # create a blank page to work with with d3 files loaded
## js(w, d3.select("body").transition(2000).style("background-color", "black"))


text:
using Test, FinancialMonteCarlo,ArrayFire;
@show "Black Scholes Model"
S0=100.0;
K=100.0;
r=0.02;
T=1.0;
d=0.01;
D=90.0;

Nsim=40000;
Nstep=30;
sigma=0.2;
mc=MonteCarloConfiguration(Nsim,Nstep,FinancialMonteCarlo.AFMode());
mc1=MonteCarloConfiguration(Nsim,Nstep,FinancialMonteCarlo.AntitheticMC(),FinancialMonteCarlo.AFMode());
toll=0.8

rfCurve=ZeroRate(r);
FwdData=Forward(T)
EUData=EuropeanOption(T,K)
AMData=AmericanOption(T,K)
BarrierData=BarrierOptionDownOut(T,K,D)
AsianFloatingStrikeData=AsianFloatingStrikeOption(T)
AsianFixedStrikeData=AsianFixedStrikeOption(T,K)
Model=BlackScholesProcess(sigma,Underlying(S0,d));

display(Model)

@show FwdPrice=pricer(Model,rfCurve,mc,FwdData);
@show EuPrice=pricer(Model,rfCurve,mc,EUData);
@show AmPrice=pricer(Model,rfCurve,mc,AMData);
@show BarrierPrice=pricer(Model,rfCurve,mc,BarrierData);
@show AsianPrice1=pricer(Model,rfCurve,mc,AsianFloatingStrikeData);
@show AsianPrice2=pricer(Model,rfCurve,mc,AsianFixedStrikeData);

@test abs(FwdPrice-99.1078451563562)<toll
@test abs(EuPrice-8.43005524824866)<toll
@test abs(AmPrice-8.450489415187354)<toll
@test abs(BarrierPrice-7.5008664470880735)<toll
@test abs(AsianPrice1-4.774451704549382)<toll

@show FwdPrice=pricer(Model,rfCurve,mc1,FwdData);
@show EuPrice=pricer(Model,rfCurve,mc1,EUData);
@show AmPrice=pricer(Model,rfCurve,mc1,AMData);
@show BarrierPrice=pricer(Model,rfCurve,mc1,BarrierData);
@show AsianPrice1=pricer(Model,rfCurve,mc1,AsianFloatingStrikeData);
@show AsianPrice2=pricer(Model,rfCurve,mc1,AsianFixedStrikeData);
tollanti=0.6;
@test abs(FwdPrice-99.1078451563562)<tollanti
@test abs(EuPrice-8.43005524824866)<tollanti
@test abs(AmPrice-8.450489415187354)<tollanti
@test abs(BarrierPrice-7.5008664470880735)<tollanti
@test abs(AsianPrice1-4.774451704549382)<tollanti

EUDataPut=EuropeanOption(T,K,false)
AMDataPut=AmericanOption(T,K,false)
BarrierDataPut=BarrierOptionDownOut(T,K,D,false)
AsianFloatingStrikeDataPut=AsianFloatingStrikeOption(T,false)
AsianFixedStrikeDataPut=AsianFixedStrikeOption(T,K,false)

@show EuPrice=pricer(Model,rfCurve,mc,EUDataPut);
@show AmPrice=pricer(Model,rfCurve,mc,AMDataPut);
@show BarrierPrice=pricer(Model,rfCurve,mc,BarrierDataPut);
@show AsianPrice1=pricer(Model,rfCurve,mc,AsianFloatingStrikeDataPut);
@show AsianPrice2=pricer(Model,rfCurve,mc,AsianFixedStrikeDataPut);
tollPut=0.6;
@test abs(FwdPrice-99.1078451563562)<tollPut
@test abs(EuPrice-7.342077422567968)<tollPut
@test abs(AmPrice-7.467065889603365)<tollPut
@test abs(BarrierPrice-0.29071929104261723)<tollPut
@test abs(AsianPrice1-4.230547012372306)<tollPut
@test abs(AsianPrice2-4.236220218194027)<tollPut

#Check Nsim even for antithetic
#mc2=MonteCarloConfiguration(Nsim+1,Nstep);
#@test_throws(ErrorException,pricer(Model,rfCurve,mc2,AsianFixedStrikeData));

@show var_std_=variance(Model,rfCurve,mc,EUDataPut)
@show var_ant_=variance(Model,rfCurve,mc1,EUDataPut);
@test var_std_>var_ant_
@show var_std_=variance(Model,rfCurve,mc,[EUDataPut,AMDataPut])
@show var_ant_=variance(Model,rfCurve,mc1,[EUDataPut,AMDataPut])
#@test prod(var_std_.>=var_ant_);

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

@show IC1=confinter(Model,rfCurve,mc,EUDataPut,0.99);
@show IC2=confinter(Model,rfCurve,mc1,EUDataPut,0.99);
L1=IC1[2]-IC1[1]
L2=IC2[2]-IC2[1]
@test L1>L2

@show IC1=confinter(Model,rfCurve,mc,[EUDataPut,AMDataPut],0.99);
@show IC2=confinter(Model,rfCurve,mc1,[EUDataPut,AMDataPut],0.99);
L1=IC1[2][2]-IC1[2][1]
L2=IC2[2][2]-IC2[2][1]
#@test L1>L2


text:
"""
    struct VectorPEPOTensor{S<:EuclideanSpace} <: AbstractTensorMap{S,3,2}

Represents a vectorised PEPO state. The `tensor` field gives the `TensorMap` representing
the state, and the `unitary` field gives the unitary who's inverse unvectorises the PEPO.
"""
struct VectorPEPOTensor{
    S<:EuclideanSpace,
    I<:Sector,
    A<:Union{<:DenseMatrix,TensorKit.SectorDict{I,<:DenseMatrix}},
    F₁,
    F₂,
} <: AbstractPhysicalTensor{S,3,2, I, A, F₁, F₂}
    tensor::TensorMap{S,3,2,I,A,F₁,F₂}
    unitary #::TensorMap{S,1,2,I,A}
    # function VectorPEPOTensor(
    #     tensor::TensorMap{S,3,2,I,A,F₁,F₂}, unitary::TensorMap{S,1,2,I,A}
    # ) where {
    #     S<:EuclideanSpace,
    #     I<:Sector,
    #     A<:Union{<:DenseMatrix,TensorKit.SectorDict{I,<:DenseMatrix}},
    #     F₁,
    #     F₂,
    # }
    #     return new{S,I,A,F₁,F₂}(tensor, unitary)
    # end
end

const VPEPOTensor = VectorPEPOTensor

max_stars_repo_name: jack-dunham/PEPOKit

max_issues_repo_name: jack-dunham/PEPOKit

max_forks_repo_name: jack-dunham/PEPOKit


text:
using wiki2
using Documenter

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

makedocs(;
    modules=[wiki2],
    authors="liuxu <liuxu89@mail.ustc.edu.cn> and contributors",
    repo="https://github.com/liuxu89/wiki2/blob/{commit}{path}#{line}",
    sitename="wiki2",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://liuxu89.github.io/wiki2",
        assets=String["src/custom.css",],
    ),
    pages=[
        "Home" => "index.md",
        "供应链金融" => Any[
            "gongjin/index.md",
            "gongjin/20200918-226.md",
        ],
    ],
)

deploydocs(;
    repo="github.com/liuxu89/wiki2",
    devbranch="main",
)


text:
import RowEchelon

#=
Homework 5
Section 2.2
Problem 1
=#

# Initialize the matrix
A = [8 6 ; 5 4]

# Create the inverse
invA = inv(A)

println("The Inverse of A is:")
println(invA)


text:
module Phonetics

include("sound.jl")
export Sound

include("VowelDensity.jl")
export Formants, VowelSpace, plot, plot!, area, vdi

include("normalize.jl")
export neareyE, neareyI, lobanov, formantWiseLogMean, nearey1, logmeanI, formantBlindLogMean, nearey2, logmeanE

include("vowelplot.jl")
export vowelPlot, ellipsePts

include("acdist.jl")
export acdist, avgseq, distinctiveness, sound2mfcc

include("util.jl")
export generateFormants

include("upt.jl")
include("vptree.jl")
include("pnd.jl")
include("phnprb.jl")

export
    pnd,
    lev,
    upt,
    phnprb,
    radiusSearch,
    nneighbors,
    TextVPTree

end # module


text:
using Measures

const assumed_ppmm = 3.78 # equivalent to 96 DPI
const px = mm/assumed_ppmm


text:
#! /usr/bin/julia

# Rosetta Code, Loops/Infinite

while true
    println("SPAM")
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
# stats function

# get mean quality scores
function Q_mean(qs,RC = size(qs))
    mean_qs = zeros(Float64,RC[1])
   # mean_qs = Vector{UInt8}(undef,RC[1])
    @turbo for col in 1:RC[2]
        #mean_qs[col] = 0
        for row in 1:RC[1]
            mean_qs[row] += qs[row,col]
        end
    end
    return round.(mean_qs ./ RC[2])
end

# calculate expected errors
function EE(q,n = length(q))
    ee = 0.0
    @turbo for i in 1:n
        ee+=10.0^(-Int(q[i])/10.0)
    end
    return ee
end

# readers
function get_reader(file)
    if endswith(file,"gz")
        if contains(file,"fastq")
            reader = FASTQ.Reader(GzipDecompressorStream(open(file)))
        else
            reader = FASTA.Reader(GzipDecompressorStream(open(file)))
        end
    else
        if contains(file,"fastq")
            reader = open(FASTQ.Reader, file)
        else
            reader = open(FASTA.Reader, file)
        end
    end
    return reader
end

function get_files(dir,identifier; just_names = false)
    file_list = []
    sample_names =[]
    for (root, dirs, files) in walkdir(dir)
        for (i,path) in enumerate(joinpath.(root, files))
            if contains(path,identifier)
            push!(file_list,path)
            push!(sample_names,files[i])
            end
        end
    end
    return just_names ? sample_names : file_list
end

# dereplicate

function dereplicate(file; sample_name = "sample")
   # println(sample_name)
    try
    seqs = Vector{String}(undef,0) #LongDNASeq

    qs = Vector{Vector{UInt8}}(undef,0)
    reader = get_reader(file)
    for record in reader
        push!(seqs,sequence(record))
        push!(qs,quality(record))
    end
    close(reader)
    sp = sortperm(seqs)
    srt = seqs[sp]#sort(seqs) # map from rank to sequence/quality/id
    unqs = unique(srt)
    qsrt=qs[sp]
    counts = countmap(srt) # map from rank to abundance (much faster than table)
    mnq = Vector{Vector{Int}}(undef,length(unqs))
    abunds= Vector{Int}(undef,length(unqs))
    sample = fill(sample_name,length(unqs))
    row = 1
    for (i,seq) in enumerate(unqs)
        abund = counts[seq]
        #mnq[i] = vec(mean(hcat(qsrt[row:row+abund-1]),dims =2 ))[1]
        mnq[i] = Q_mean(hcat(qsrt[row:row+abund-1]...))
        abunds[i] = abund
        row += abund
    end
    df =DataFrame([sample,unqs,abunds,mnq],[:sample,:sequence,:count,:quality])
    sort!(df,:count,rev = true)
    unq_map = Dict([df.sequence[i] => i for i in 1:length(unqs)])
    read_map = Vector{Int}(undef,length(seqs))
    for i in 1:length(seqs)
        read_map[i] = unq_map[seqs[i]]
    end

    return (sample = sample_name,df =df, read_map = read_map, sequence = df.sequence,count = df.count,quality = df.quality)
    catch e
        println("Error encountered with file $(sample_name)")
        println(e)
    end
end

function dereplicate(dir, identifier)
    file_list = []
    sample_names =[]
    for (root, dirs, files) in walkdir(dir)
        for (i,path) in enumerate(joinpath.(root, files))
            if contains(path,identifier)
            push!(file_list,path)
            push!(sample_names,files[i])
            end
        end
    end
    n = length(file_list)
    dfs = []#Vector{DataFrame}(undef,0)
    for i in 1:n #  not sure why multi-threading slows this down

        try
            push!(dfs,dereplicate(file_list[i], sample_name = sample_names[i]))
        catch e
            println("Error encountered with file $(file_list[i])")
            println(e)
        end
    end
return vcat(dfs...)
end

#trim and filter

function trim(file,out_path ="my-out.fastq";trunc_quality=0, trunc_length=Inf,max_length=Inf,min_length=0,trim_left=0,
    trim_right=0,max_N =0,min_quality = 0, max_expected_errors = Inf)

    trim_left+=1
    reader = get_reader(file)
    println(out_path)
    writer = open(FASTQ.Writer, out_path)
    for record in reader
        seq = sequence(record)
        q = quality(record)
        if length(seq) <=max_length                 #max length
            seq = seq[trim_left:(end-trim_right)]     # trim left/right
            q = q[trim_left:(end-trim_right)]
            if trunc_quality > 0                    # trunc quality
                tq = findfirst(x-> x <trunc_quality,q)
                if tq !== nothing
                    seq = seq[1:tq]
                    q = q[1:tq]
                end
            end
            if trunc_length <= length(seq)          # trunc length
                seq = seq[1:trunc_length]
                q = q[1:trunc_length]
                if min_length <= length(seq)        #min length
                    if count(isambiguous, seq) <= max_N # max n
                        if minimum(q) >= min_quality # min quality
                            if EE(q) <= max_expected_errors # max expected errors
                                write(writer,FASTQ.Record("id", seq, q))
                            end
                        end
                    end
                end
            end
        end
    end
    close(writer)
    close(reader)
end

function trim(dir,identifier,out_dir;trunc_quality=0, trunc_length=Inf,max_length=Inf,min_length=0,trim_left=0,
    trim_right=0,max_N =0,min_quality = 0, max_expected_errors = Inf)
    file_list = []
    sample_names =[]
    for (root, dirs, files) in walkdir(dir)
        for (i,path) in enumerate(joinpath.(root, files))
            if contains(path,identifier)
            push!(file_list,path)
            push!(sample_names,files[i])
            end
        end
    end
    n = length(file_list)
    for i in 1:n

        try
        file = file_list[i]
        out_path = joinpath(splitpath(out_dir)..., sample_names[i][1:end-3])
        trim(file,out_path,trunc_quality=trunc_quality, trunc_length=trunc_length,
        max_length=max_length,min_length=min_length,trim_left=trim_left,
        trim_right=trim_right,max_N =max_N,min_quality = min_quality, max_expected_errors = max_expected_errors)
        catch e
            println("Error encountered with file $(file_list[i])")
            println(e)
        end
    end
end


text:
baremodule Recyclers

###
### Generic API functions and types
###

function get! end
function maybeget! end

function recycle! end
function recycling! end

function unsafe_empty! end

function unsafe_init! end
function unsafe_destroy! end

abstract type AbstractCache{T} end
abstract type AbstractRecycler{T} end

function cacheof end
function factoryof end

function var"@global" end

"""
    InternalPrelude

Internal module. Define types and utilities just enough for defining public types.
"""
module InternalPrelude
using ..Recyclers: @global
const var"@_global" = var"@global"
include("prelude.jl")
end  # module InternalPrelude

###
### Concrete public types
###

struct ShardedCache{T} <: AbstractCache{T}
    limit::UInt
    eltype::InternalPrelude.Val{T}
    pools::InternalPrelude.Vector{InternalPrelude.Vector{T}}  # TODO: blocked stack?
end

struct ShardedRecycler{T,Factory} <: AbstractRecycler{T}
    factory::Factory
    cache::ShardedCache{T}
end

InternalPrelude.@if_has_atomic begin
    mutable struct CentralizedCache{T} <: AbstractCache{T}
        # pad?
        InternalPrelude.@atomic top::Union{Nothing,InternalPrelude.StackNode{T}}
        # cache::shardedCache{StackNode{T}}

        InternalPrelude._CentralizedCache(
            ::Type{T},
            top::Union{Nothing,InternalPrelude.StackNode{T}},
        ) where {T} = new{T}(top)
    end

    struct CentralizedRecycler{T,Factory} <: AbstractRecycler{T}
        factory::Factory
        cache::CentralizedCache{T}
    end
end

"""
    Internal

Internal module that contains main implementations.
"""
module Internal

using ..Recyclers: AbstractCache, AbstractRecycler, Recyclers, cacheof, factoryof
using ..InternalPrelude: HAS_ATOMIC, StackNode, _CentralizedCache

import ..InternalPrelude: @_global
import ..Recyclers: ShardedCache, ShardedRecycler
if HAS_ATOMIC
    import ..Recyclers: CentralizedCache, CentralizedRecycler
end

include("utils.jl")
include("generic.jl")
include("sharded.jl")
include("global.jl")

if HAS_ATOMIC
    include("centralized.jl")
end

end  # module Internal

# Something like this could be useful?
#
#     if InternalPrelude.HAS_ATOMIC
#         const Cache = CentralizedCache
#         const Recycler = CentralizedRecycler
#     else
#         const Cache = ShardedCache
#         const Recycler = ShardedRecycler
#     end
#
# But it will be problematic if people rely on `Cache(; limit = ...)` etc. OTOH,
# using a factory or abstract class for `Recycler` is not great since it's hard
# to use concrete recycler type on field or element type. Maybe a better
# midground is to define `default_recycler_type` or something to return a
# default recommended type. It's more apparent that the users have to call the
# common API.

Internal.define_docstrings()

end  # baremodule Recyclers


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using DataFrames, CSV, Statistics
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.10"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	data = CSV.read(datadir("purchaseData.csv"))
	rowsKeep = .!(ismissing.(data.Grade) .& ismissing.(data.Price))
	data = data[rowsKeep, :]
end

# ╔═╡ 02c4946e-2445-11eb-3c51-b550096ab843
begin
	replace!(x -> ismissing(x) ? "QQ" : x, data.Name)
	replace!(x -> ismissing(x) ? "31/06/2008" : x, data.Date)

	grPr = by(dropmissing(data), :Grade, :Price=>x ->
		AvgPrice = round(mean(x), digits=-3))

	d = Dict(grPr[:,1] .=> grPr[:,2])
	nearIndx(v, x) = findmin(abs.(v.-x))[2]
	for i in 1:nrow(data)
		if ismissing(data[i, :Price])
			data[i, :Price] = d[data[i, :Grade]]
		end
		if ismissing(data[i, :Grade])
			data[i, :Grade] = grPr[ nearIndx(grPr[:,2], data[i, :Price]), :Grade]
		end
	end
end

# ╔═╡ 02c50672-2445-11eb-03bc-d763278b5f14
first(data, 5)

# ╔═╡ 02d285fe-2445-11eb-1150-e9fa017e6170
describe(data)

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.10"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═02c4946e-2445-11eb-3c51-b550096ab843
# ╠═02c50672-2445-11eb-03bc-d763278b5f14
# ╠═02d285fe-2445-11eb-1150-e9fa017e6170
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
struct Vertex
    m_position::AbstractArray{Float32}
    m_uv_coord::AbstractArray{Float32}
    m_normal::AbstractArray{Float32}
    m_color::AbstractArray{Float32}

    Vertex() = new([0, 0, 0], [0, 0], [0, 0, 0], [0, 0, 0, 1])
end

struct Texture
    m_texture::UInt32
    m_width::Int32
    m_height::Int32
    m_nr_channels::Int32
	m_data::Ptr{UInt8}
    m_shape::Int32

    Texture() = new(0, 0, 0, 0, C_NULL, 0)
end

struct Material

end


text:
function DiffEqBase.__solve(prob::AbstractSDDEProblem, # TODO: DiffEqBase.AbstractSDDEProblem
    alg::AbstractSDEAlgorithm, args...; # TODO: Method of steps???
    kwargs...)
    integrator = DiffEqBase.__init(prob, alg, args...; kwargs...)
    DiffEqBase.solve!(integrator)
    integrator.sol
end

function DiffEqBase.__init(prob::AbstractSDDEProblem,# TODO DiffEqBasee.AbstractSDDEProblem
    alg::Union{AbstractRODEAlgorithm,AbstractSDEAlgorithm},
    timeseries_init = typeof(prob.u0)[],
    ts_init = eltype(prob.tspan)[],
    ks_init = nothing,
    recompile::Type{Val{recompile_flag}} = Val{true};
    saveat = eltype(prob.tspan)[],
    tstops = eltype(prob.tspan)[],
    d_discontinuities = Discontinuity{eltype(prob.tspan)}[],
    save_idxs = nothing,
    save_everystep = isempty(saveat),
    save_noise = save_everystep && (typeof(prob.f) <: Tuple ?
                 DiffEqBase.has_analytic(prob.f[1]) : DiffEqBase.has_analytic(prob.f)),
    save_on = true,
    save_start = save_everystep || isempty(saveat) || typeof(saveat) <: Number ? true : prob.tspan[1] in saveat,
    save_end = save_everystep || isempty(saveat) || typeof(saveat) <: Number ? true : prob.tspan[2] in saveat,
    callback = nothing,
    dense = save_everystep && isempty(saveat),
    calck = (!isempty(setdiff(saveat, tstops)) || dense),
    dt = eltype(prob.tspan)(0),
    adaptive = StochasticDiffEq.isadaptive(getalg(alg)),
    gamma = 9 // 10, # TODO gamma_default(alg.alg) ?
    abstol = nothing,
    reltol = nothing,
    qmax = StochasticDiffEq.qmax_default(getalg(alg)),
    qmin = StochasticDiffEq.qmin_default(getalg(alg)),
    qoldinit = 1 // 10^4, fullnormalize = true,
    failfactor = 2,
    beta2 = StochasticDiffEq.beta2_default(getalg(alg)),
    beta1 = StochasticDiffEq.beta1_default(getalg(alg), beta2),
    delta = StochasticDiffEq.delta_default(getalg(alg)),
    maxiters = adaptive ? 1000000 : typemax(Int),
    dtmax = eltype(prob.tspan)((prob.tspan[end] - prob.tspan[1])),
    dtmin = typeof(one(eltype(prob.tspan))) <: AbstractFloat ? eps(eltype(prob.tspan)) :
            typeof(one(eltype(prob.tspan))) <: Integer ? 0 :
            eltype(prob.tspan)(1 // 10^(10)),
    internalnorm = DiffEqBase.ODE_DEFAULT_NORM,
    isoutofdomain = DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN,
    unstable_check = DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK,
    verbose = true,force_dtmin = false,
    timeseries_errors = true, dense_errors = false,
    advance_to_tstop = false,stop_at_next_tstop = false,
    initialize_save = true,
    progress = false,progress_steps = 1000,progress_name = "SDDE",
    progress_message = DiffEqBase.ODE_DEFAULT_PROG_MESSAGE,
    userdata = nothing,
    initialize_integrator = true,
    seed = UInt64(0), alias_u0 = false,
    #  Keywords for Delay problems (from DDE)
    discontinuity_interp_points::Int = 10,
    discontinuity_abstol = eltype(prob.tspan)(1 // Int64(10)^12),
    discontinuity_reltol = 0, kwargs...) where recompile_flag

    # alg = getalg(alg0);
    if typeof(prob.f) <: Tuple
        if any(mm != I for mm in prob.f.mass_matrix)
            error("This solver is not able to use mass matrices.")
        end
    elseif prob.f.mass_matrix != I && !alg_mass_matrix_compatible(getalg(alg))
        error("This solver is not able to use mass matrices.")
    end

    if !isempty(saveat) && dense
        @warn("Dense output is incompatible with saveat. Please use the SavingCallback from the Callback Library to mix the two behaviors.")
    end

    if typeof(prob.noise) <: NoiseProcess && prob.noise.bridge === nothing && adaptive
        error("Bridge function must be given for adaptivity. Either declare this function in noise process or set adaptive=false")
    end

    if !StochasticDiffEq.alg_compatible(prob, getalg(alg))
        error("The algorithm is not compatible with the chosen noise type. Please see the documentation on the solver methods")
    end

    if haskey(kwargs, :initial_order)
        @warn "initial_order has been deprecated. Please specify order_discontinuity_t0 in the DDEProblem instead."
        order_discontinuity_t0 = kwargs[:initial_order]
    else
        order_discontinuity_t0 = prob.order_discontinuity_t0
    end

    if haskey(kwargs, :minimal_solution)
        @warn "minimal_solution is ignored"
    end
    progress && @logmsg(-1,progress_name,_id = _id = :StochasticDiffEq,progress = 0)

    tType = eltype(prob.tspan)
    noise = prob.noise
    tspan = prob.tspan
    tdir = sign(tspan[end] - tspan[1])

    t = tspan[1]

    if !adaptive && iszero(dt) && isempty(tstops)
        error("Fixed timestep methods require a choice of dt or choosing the tstops")
    end
    if has_dependent_lags(prob)
        @warn "Discontinuities are not handled for dependent lags. Please consider it when evaluating the soution!"
    end

    @unpack f, g, u0, h, tspan, p, constant_lags, dependent_lags, neutral = prob
    # g = prob isa AbstractSDDEProblem ? prob.g : nothing # TODO DiffEqBase.AbstractSDDEProblem

    # no fixed-point iterations for stochastic delay problems,
    # and thus `dtmax` should match minimal lag
    if  has_constant_lags(prob) # && isconstrained(alg)
        dtmax = tdir * min(abs(dtmax), minimum(abs, constant_lags))
    end

    u, uprev = u_uprev(prob.u0, alias_u0 = alias_u0)

    uType = typeof(u)
    uBottomEltype = recursive_bottom_eltype(u)
    uBottomEltypeNoUnits = recursive_unitless_bottom_eltype(u)

    ks = Vector{uType}(undef, 0)

    uEltypeNoUnits = recursive_unitless_eltype(u)
    tTypeNoUnits   = typeof(one(tType))

    if abstol === nothing
        if uBottomEltypeNoUnits == uBottomEltype
            abstol_internal = real(convert(uBottomEltype, oneunit(uBottomEltype) * 1 // 10^2))
        else
            abstol_internal = real.(oneunit.(u) .* 1 // 10^2)
        end
    else
        abstol_internal = real.(abstol)
    end

    if reltol === nothing
        if uBottomEltypeNoUnits == uBottomEltype
            reltol_internal = real(convert(uBottomEltype, oneunit(uBottomEltype) * 1 // 10^2))
        else
            reltol_internal = real.(oneunit.(u) .* 1 // 10^2)
        end
    else
        reltol_internal = real.(reltol)
    end

    if isinplace(prob) && typeof(u) <: AbstractArray && eltype(u) <: Number && uBottomEltypeNoUnits == uBottomEltype # Could this be more efficient for other arrays?
        if !(typeof(u) <: ArrayPartition)
            rate_prototype = recursivecopy(u)
        else
            rate_prototype = similar(u, typeof.(oneunit.(recursive_bottom_eltype.(u.x)) ./ oneunit(tType))...)
        end
    else
        if uBottomEltypeNoUnits == uBottomEltype
            rate_prototype = u
        else # has units!
            rate_prototype = u / oneunit(tType)
        end
    end
    rateType = typeof(rate_prototype) ## Can be different if united

    if StochasticDiffEq.is_diagonal_noise(prob)
        noise_rate_prototype = rate_prototype
    else
        if prob isa AbstractSDDEProblem # TODO DiffEqBase.AbstractSDDEProblem
            noise_rate_prototype = copy(prob.noise_rate_prototype)
        else
            noise_rate_prototype = copy(prob.rand_prototype)
        end
    end

    # tstops_internal, saveat_internal, d_discontinuities_internal =
    #   tstop_saveat_disc_handling(tstops, saveat, d_discontinuities, tspan) # TODO add delays to discontinuities
    # retrieve time stops, time points at which solutions is saved, and discontinuities
    maximum_order = StochasticDiffEq.alg_order(getalg(alg))
    tstops_internal, saveat_internal, d_discontinuities_internal =
      tstop_saveat_disc_handling(tstops, saveat, d_discontinuities, tspan,
                                                order_discontinuity_t0, maximum_order,
                                                constant_lags, neutral)

    tracked_discontinuities = Discontinuity{tType}[]
    if order_discontinuity_t0 ≤ maximum_order
        push!(tracked_discontinuities, Discontinuity(tdir * t, order_discontinuity_t0))
    end

    callbacks_internal = CallbackSet(callback)

    max_len_cb = DiffEqBase.max_vector_callback_length(callbacks_internal)
    if max_len_cb isa DiffEqBase.VectorContinuousCallback
        callback_cache = DiffEqBase.CallbackCache(max_len_cb.len, uBottomEltype, uBottomEltype)
    else
        callback_cache = nothing
    end

    QT = tTypeNoUnits <: Integer ? typeof(qmin) : tTypeNoUnits

    if !(uType <: AbstractArray)
        rand_prototype = zero(u / u) # Strip units and type info
        randType = typeof(rand_prototype)
    else
        randElType = uBottomEltypeNoUnits # Strip units and type info
        if StochasticDiffEq.is_diagonal_noise(prob)
            if typeof(u) <: SArray
                rand_prototype = zero(u) # TODO: Array{randElType} for units
            else
                rand_prototype = (u .- u) ./ sqrt(oneunit(t))
            end
        elseif prob isa AbstractSDDEProblem # TODO DiffEqBase.AbstractSDDEProblem
            rand_prototype = false .* noise_rate_prototype[1,:]
        else
            rand_prototype = copy(prob.rand_prototype)
        end
        randType = typeof(rand_prototype) # Strip units and type info
    end

    _seed = iszero(seed) ? (iszero(prob.seed) ? rand(UInt64) : prob.seed) : seed

    if prob.noise === nothing
        rswm = StochasticDiffEq.isadaptive(getalg(alg)) ? RSWM(adaptivealg = :RSwM3) : RSWM(adaptivealg = :RSwM1)
        if isinplace(prob)
            if StochasticDiffEq.alg_needs_extra_process(getalg(alg))
                W = WienerProcess!(t, rand_prototype, rand_prototype,
                             save_everystep = save_noise,
                             rswm = rswm,
                             rng = Xorshifts.Xoroshiro128Plus(_seed))
            else
                W = WienerProcess!(t, rand_prototype,
                             save_everystep = save_noise,
                             rswm = rswm,
                             rng = Xorshifts.Xoroshiro128Plus(_seed))
            end
        else
            if StochasticDiffEq.alg_needs_extra_process(getalg(alg))
                W = WienerProcess(t, rand_prototype, rand_prototype,
                             save_everystep = save_noise,
                             rswm = rswm,
                             rng = Xorshifts.Xoroshiro128Plus(_seed))
            else
                W = WienerProcess(t, rand_prototype,
                             save_everystep = save_noise,
                             rswm = rswm,
                             rng = Xorshifts.Xoroshiro128Plus(_seed))
            end
        end
    else
        W = prob.noise
        if W.reset
            if W.t[end] != t
                reinit!(W, t)
            end
        # Reseed
            if typeof(W) <: NoiseProcess && W.reseed
                Random.seed!(W.rng, _seed)
            end
        elseif W.t[end] != t
            error("Starting time in the noise process is not the starting time of the simulation. The noise process should be re-initialized for repeated use")
        end
    end

    ts, timeseries, saveiter = solution_arrays(u, tspan, rate_prototype,
                    timeseries_init = timeseries_init, ts_init = ts_init, save_idxs = save_idxs, save_start = save_start)

    if !adaptive && save_everystep && tspan[2] - tspan[1] != Inf
        iszero(dt) ? steps = length(tstops) : steps = ceil(Int, internalnorm((tspan[2] - tspan[1]) / dt, tspan[1]))
        sizehint!(timeseries, steps + 1)
        sizehint!(ts, steps + 1)
    elseif save_everystep
        sizehint!(timeseries, 50)
        sizehint!(ts, 50)
    elseif !isempty(saveat_internal)
        sizehint!(timeseries, length(saveat_internal) + 1)
        sizehint!(ts, length(saveat_internal) + 1)
    else
        sizehint!(timeseries, 2)
        sizehint!(ts, 2)
    end

    alg_choice = Int[]
    if save_start && typeof(getalg(alg)) <: StochasticDiffEq.StochasticDiffEqCompositeAlgorithm
        copyat_or_push!(alg_choice, 1, 1)
    end

    # create a history function
    history = build_history_function(prob, alg, reltol_internal,
                                    rate_prototype, noise_rate_prototype, W, _seed, dense;
                                    dt = dt, adaptive = adaptive,
                                    internalnorm = internalnorm)
    f_with_history, g_with_history = wrap_functions_and_history(f, g, history)

    sde_integrator = history.integrator;

    cache = StochasticDiffEq.alg_cache(getalg(alg), prob, u, W.dW, W.dZ, p, rate_prototype, noise_rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, uprev, f_with_history, t, dt, Val{isinplace(prob)})

    # id = StochasticDiffEq.LinearInterpolationData(timeseries,ts)
    id = StochasticDiffEq.LinearInterpolationData(sde_integrator.sol.u, sde_integrator.sol.t)

    opts = StochasticDiffEq.SDEOptions(maxiters, save_everystep,
                      adaptive, abstol_internal,
                      reltol_internal, QT(gamma),
                      QT(qmax), QT(qmin),
                      QT(failfactor),
                      tType(dtmax), tType(dtmin), internalnorm, save_idxs,
                      tstops_internal, saveat_internal,
                      d_discontinuities_internal,
                      tstops, saveat, d_discontinuities,
                      userdata,
                      progress, progress_steps,
                      progress_name, progress_message,
                      timeseries_errors, dense_errors,
                      QT(beta1), QT(beta2),
                      convert.(uBottomEltypeNoUnits, delta),
                      QT(qoldinit),
                      dense, save_on, save_start, save_end, save_noise,
                      callbacks_internal, isoutofdomain, unstable_check,
                      verbose, calck, force_dtmin,
                      advance_to_tstop, stop_at_next_tstop)

    if typeof(getalg(alg)) <: StochasticDiffEq.StochasticDiffEqCompositeAlgorithm
      # TODO: DISCONNECT!!!!
        sol =  DiffEqBase.build_solution(prob, alg, sde_integrator.sol.t, sde_integrator.sol.u, W = W,
                                      destats = DiffEqBase.DEStats(0),
                                      calculate_error = false, alg_choice = alg_choice,
                                      interp = id, dense = dense, seed = _seed)
    # separate statistics of the integrator and the history
    else
      # TODO: DISCONNECT!!!!
        sol = DiffEqBase.build_solution(prob, alg, sde_integrator.sol.t, sde_integrator.sol.u, W = W,
                                      destats = DiffEqBase.DEStats(0),
                                      calculate_error = false,
                                      interp = id, dense = dense, seed = _seed)
    # separate statistics of the integrator and the history
    end

    if recompile_flag == true
        FType = typeof(f_with_history)
        GType = typeof(g_with_history)
        SolType = typeof(sol)
        cacheType = typeof(cache)
    else
        FType = Function
        GType = Function
        SolType = DiffEqBase.AbstractRODESolution
        cacheType = StochasticDiffEq.StochasticDiffEqCache
    end

    tprev = t
    dtcache = tType(dt)
    iter = 0
    u_modified = false
    eigen_est = one(uBottomEltypeNoUnits) / oneunit(tType) # rate/state = (state/time)/state = 1/t units
    EEst = tTypeNoUnits(1)
    just_hit_tstop = false
    isout = false
    accept_step = false
    force_stepfail = false
    last_stepfail = false
    event_last_time = 0
    vector_event_last_time = 1
    last_event_error = zero(uBottomEltypeNoUnits)
    dtchangeable = true
    q11 = tTypeNoUnits(1)
    success_iter = 0
    q = tTypeNoUnits(1)

    integrator = SDDEIntegrator{typeof(getalg(alg)),isinplace(prob),uType,uBottomEltype,tType,typeof(p),typeof(eigen_est),QT,uEltypeNoUnits,typeof(W),rateType,typeof(sol),typeof(cache),FType,GType,typeof(opts),typeof(noise),typeof(last_event_error),typeof(callback_cache),typeof(history),typeof(sde_integrator)}(f_with_history, g_with_history, noise, uprev, tprev,
                      order_discontinuity_t0, tracked_discontinuities,
                      t, u, p, tType(dt), tType(dt), tType(dt), dtcache, tspan[2], tdir,
                      just_hit_tstop, isout, event_last_time, vector_event_last_time, last_event_error, accept_step,
                      last_stepfail, force_stepfail, dtchangeable, u_modified, saveiter, getalg(alg), sol,
                      cache, callback_cache, tType(dt), W,
                      opts, iter, success_iter, eigen_est, EEst, q, QT(qoldinit), q11, history, sde_integrator)

    if initialize_integrator
        StochasticDiffEq.initialize_callbacks!(integrator, initialize_save)
        initialize!(integrator, integrator.cache)

        save_start && typeof(alg) <: StochasticDiffEq.StochasticDiffEqCompositeAlgorithm && copyat_or_push!(alg_choice, 1, integrator.cache.current)
    end

    StochasticDiffEq.handle_dt!(integrator)

    ## Modify the first dt for tstops
    StochasticDiffEq.modify_dt_for_tstops!(integrator)
    ### Needs to be done before first rand
    integrator.sqdt = integrator.tdir * sqrt(abs(integrator.dt))

    integrator.W.dt = integrator.dt
    DiffEqNoiseProcess.setup_next_step!(integrator.W)

    integrator
end

function DiffEqBase.solve!(integrator::SDDEIntegrator)
    @inbounds while !isempty(integrator.opts.tstops)
        while integrator.tdir * integrator.t < top(integrator.opts.tstops)
            loopheader!(integrator)
            if DiffEqBase.check_error!(integrator) != :Success
                return integrator.sol
            end
            StochasticDiffEq.perform_step!(integrator, integrator.cache)
            loopfooter!(integrator)
            if isempty(integrator.opts.tstops)
                break
            end
        end
        @inbounds handle_tstop!(integrator)
    end
    postamble!(integrator)

    f = typeof(integrator.sol.prob.f) <: Tuple ? integrator.sol.prob.f[1] : integrator.sol.prob.f

    if DiffEqBase.has_analytic(f)
        DiffEqBase.calculate_solution_errors!(integrator.sol;timeseries_errors = integrator.opts.timeseries_errors,dense_errors = integrator.opts.dense_errors)
    end
    if integrator.sol.retcode != :Default
        return integrator.sol
    end
    integrator.sol = DiffEqBase.solution_new_retcode(integrator.sol, :Success)
end

# function StochasticDiffEq.tstop_saveat_disc_handling(tstops, saveat, d_discontinuities, tspan)
function tstop_saveat_disc_handling(tstops, saveat, d_discontinuities, tspan, order_discontinuity_t0, alg_maximum_order, constant_lags, neutral)

    tType = eltype(tspan)
    t0, tf = tspan
    tdir = sign(tf - t0)
    tdir_t0 = tdir * t0
    tdir_tf = tdir * tf

  # add discontinuities propagated from initial discontinuity
    add_propagated_constant_lags = order_discontinuity_t0 ≤ alg_maximum_order && constant_lags !== nothing && !isempty(constant_lags)

    if add_propagated_constant_lags
        maxlag = tdir_tf - tdir_t0
        next_order = neutral ? order_discontinuity_t0 : order_discontinuity_t0 + 1//2
    end

  # time stops
    tstops_internal = BinaryMinHeap{tType}()
    if isempty(d_discontinuities) && !add_propagated_constant_lags && isempty(tstops) # TODO: Specialize more
        push!(tstops_internal, tdir_tf)
    else
        for t in tstops
            tdir_t = tdir * t
            tdir_t0 < tdir_t ≤ tdir_tf && push!(tstops_internal, tdir_t)
        end

        for t in d_discontinuities
            tdir_t = tdir * t
            tdir_t0 < tdir_t ≤ tdir_tf && push!(tstops_internal, tdir_t)
        end

  # add propagated discontinuities
        if add_propagated_constant_lags
            for lag in constant_lags
                if tdir * lag < maxlag
                    push!(tstops_internal, tdir * (t0 + lag))
                end
            end
        end

        push!(tstops_internal, tdir_tf)
    end

  # saving time points
    saveat_internal = BinaryMinHeap{tType}()
    if typeof(saveat) <: Number
        if (t0:saveat:tf)[end] == tf
            for t in (t0 + saveat):saveat:tf
                push!(saveat_internal, tdir * t)
            end
        else
            for t in (t0 + saveat):saveat:(tf - saveat)
                push!(saveat_internal, tdir * t)
            end
        end
    elseif !isempty(saveat)
        for t in saveat
            tdir_t = tdir * t
            tdir_t0 < tdir_t < tdir_tf && push!(saveat_internal, tdir_t)
        end
    end

  # discontinuities
    d_discontinuities_internal = BinaryMinHeap{Discontinuity{tType}}()
    if add_propagated_constant_lags
        sizehint!(d_discontinuities_internal.valtree, length(d_discontinuities) + length(constant_lags))
    else
        sizehint!(d_discontinuities_internal.valtree, length(d_discontinuities))
    end

    for d in d_discontinuities
        tdir_t = tdir * d.t

        if tdir_t0 < tdir_t < tdir_tf && d.order ≤ alg_maximum_order + 1
            push!(d_discontinuities_internal, Discontinuity{tType}(tdir_t, d.order))
        end
    end

    if add_propagated_constant_lags
        for lag in constant_lags
            if tdir * lag < maxlag
                push!(d_discontinuities_internal, Discontinuity{tType}(tdir * (t0 + lag), next_order))
            end
        end
    end

    tstops_internal, saveat_internal, d_discontinuities_internal
end


text:
using Qsosed, Test

@testset "Compton test" begin
    parameters = [2.5, 100, 50, 0, 0]
    ear = 10 .^ range(-4, 2, length=50)
    @test compute_compton_photons_per_bin(ear, parameters) == donthcomp(ear, parameters)
end


text:
"""
固定書式ファイルからデータフレームへ入力するプログラム

このプログラムを fixedformat.jl という名前で保存しておき，
使用するときに

include("fixedformat.jl")

でインクルードする。

使用法-1

fixedformat(fn; settype=true)

データファイルの 1, 2 行目に読み取り方を示す情報を記述しておく場合

データファイルの1行目に，カンマで区切った変数名
2行目には，読み込む変数の開始桁と変数の型を下記のように指定する
i は整数値 欄内に符号があってもよい
f は実数値 欄内に符号，小数点があってもよい
s は文字列
x は読み飛ばす
文字のある桁位置から次の文字のある桁位置の前までを文字の種類に応じて読む
以下の例では
1〜5   の 5 文字を整数値として読む
6〜13  の 8 文字を実数値として読む
14〜21 の 8 文字を文字列として読む
32〜35 の 4 文字を読み飛ばす...　ことを表している。
データは デフォルトでは 3 行目から
i1, f1, s1, i2, i3, f2
i    f       s       x   i i  f
1234567890123456789012345678901234567890
67897978.7827897289078329738290878696899
7897807890.80780780780780676786767967867

使用例-1-1

include("fixedformat.jl") # include はセッションごとに 1 回でよい
fn = "sample1.dat"
df = fixedformat(fn)

  i1: (i)  1 - 5 [5]
  f1: (f)  6 - 13 [8]
  s1: (s)  14 - 21 [8]
SKIP: (x)  22 - 25 [4]
  i2: (i)  26 - 27 [2]
  i3: (i)  28 - 30 [3]
  f2: (f)  31 - 39 [9]

4 rows × 6 columns

   i1     f1         s1        i2     i3     f2
   Int64  Float64    String    Int64  Int64  Float64
1  12345  6.78901e7  45678901  67     890    1.23457e8
2  67897  978.783    89728907  73     829    8.78697e7
3  78978  7890.8     78078078  67     678    6.76797e8

使用例-1-2

include("fixedformat.jl") # include はセッションごとに 1 回でよい
df = fixedformat(fn, settype=false)

settype=false を指定すると，各列の変数型を Any のままにする。

4 rows × 6 columns

   i1     f1         s1        i2   i3     f2
   Any    Any        Any       Any  Any    Any
1  12345  6.78901e7  45678901  67   890    1.23457e8
2  67897  978.783    89728907  73   829    8.78697e7
3  78978  7890.8     78078078  67   678    6.76797e8

使用法-2

include("fixedformat.jl") # include はセッションごとに 1 回でよい
fixedformat(fn; datarow=1, settype=true, names=String[], begins=Int[], ends=Int[], types=Char[])

変数名 names，開始桁位置 begins，終了桁位置 ends，変数の型 types をベクトルで与える場合

使用例-1 と同じことをするための関数呼び出しを例示する。
なお，上と同じデータファイルを使用する場合は，データの始まる行を datarow で指定しなければならない。
names は "name1" のようなダブルクオートで囲んだ文字列，types は 'i' のような文字型（シングルクオートで囲んだ一文字）で与えることに注意。

df = fixedformat(fn, datarow=3,
    names=["Int1", "Float1", "String", "Int2", "Int3", "Float2"],
    begins=[1, 6, 14, 26, 28, 31],
    ends=[5, 13, 21, 27, 30, 39],
    types=['i', 'f', 's', 'i', 'i', 'f'])

"""

using DataFrames
function getformat(format)
    begins = Int[]
    types  = Char[]
    for (b, t) in enumerate(format)
        isnothing(indexin(t, ['i', 'f', 's', 'x'])[1]) && continue
        append!(begins, b)
        append!(types,  t)
    end
    # println(begins)
    # println(types)
    ends = begins .- 1
    popfirst!(ends)
    push!(ends, length(format))
    (begins, ends, types)
end

function getdata(txt, begins, ends, types)
    data = Any[]
    for (b, e, t) in zip(begins, ends, types)
        if t == 'i'
            append!(data, parse(Int, txt[b:e]))
        elseif t == 'f'
            append!(data, parse(Float64, txt[b:e]))
        elseif t == 's'
            append!(data, [txt[b:e]])
        end
    end
    DataFrame(reshape(data, 1, :), :auto)
end

function fixedformat(fn; datarow=1, settype=true, names=String[], begins=Int[], ends=Int[], types=Char[])
    f = open(fn, "r")
    n = length(names)
    if n == 0
        names = split(readline(f), r"[, ]+")
        format = readline(f)
        # println(names)
        # println(format)
        (begins, ends, types) = getformat(format)
    else
        length(names) == length(begins) == length(ends) == length(types) || return "names, begins, ends, types の要素数が違います"
        for _ in 1:datarow-1
            readline(f)
        end
    end
    n = length(names)
    w = max(maximum(length.(names)), 5)
    j = 0
    for i in 1:length(begins)
        if types[i] != 'x'
            j += 1
            name = ((" " ^ (w)) * names[j])[end-w+1:end]
        else
            name = ((" " ^ (w)) * "SKIP")[end-w+1:end]
        end
        println("$name: ($(types[i]))  $(begins[i]) - $(ends[i]) [$(ends[i] - begins[i] + 1)]")
    end
    df = DataFrame[]
    while (txt = readline(f)) != ""
        if df == []
            df = getdata(txt, begins, ends, types)
            first = false
        else
            append!(df, getdata(txt, begins, ends, types))
        end
        # println(data)
    end
    close(f)
    rename!(df, names)
    settype || return df
    for i in 1:length(names)
        type = typeof(df[1, i])
        df[!, i] = type.(df[:, i])
        # println("$i, $type")
    end
    df
end

max_stars_repo_name: r-de-r/stats

max_issues_repo_name: r-de-r/stats

max_forks_repo_name: r-de-r/stats


text:
function rod_image(xcen, ycen, length, pa, fwhm; nx=256)
    # returns image of a rod with some (x,y) position, length,
    #   position angle, and FWHM in an (nx,nx) image.
    rodimage = zeros(nx,nx)

    x=0
    y=0
    sig = fwhm/2.355
    dtor = π/180
    # -------- define a unit vector in direction of rod at position angle pa
    ux = sin(pa*dtor)   #  0 deg is up
    uy = cos(pa*dtor)   # 90 deg to the right

    for i=1:nx
        for j=1:nx
            x=i-nx/2+xcen
            y=j-nx/2+ycen

            # -------- distance parallel and perpendicular to
            dpara =  ux*x + uy*y
            dperp = -uy*x + ux*y

            if abs(dpara)-length <0
                dpara= 0
            end
            dpara = abs(dpara)
            dpara = min(abs(dpara-length),dpara)

            rodimage[i,j] = exp(-(dperp^2+dpara^2)/(2*sig^2))
        end
    end
    rodimage ./= sqrt(sum(rodimage.^2))
    return rodimage
end

@testset "std_use" begin
    a = rand(256,256)
    ref = fink_filter_hash(1,8)
    orig = eqws_compute(a,ref)
    @test size(orig)[1] == 2452

    rod_test = rod_image(1,1,30,35,8)
    filter_hash = fink_filter_hash(1,8)
    test = eqws_compute(rod_test,filter_hash)
    #check well sampled images sum layer by layer to 1
    @test abs(sum(test[3:end])-2) < 1e-6
end

@testset "p_power" begin
    filter_hash = fink_filter_hash(1,8)
    test_hash = fink_filter_hash_p(1,8)

    a = rand(256,256)
    orig = eqws_compute(a,filter_hash)
    new = eqws_compute_p(a,test_hash)
    @test orig ≈ new
end

@testset "conv_maps" begin
    a = rand(256,256)
    filter_hash = fink_filter_hash(1,8)
    d = eqws_compute_convmap(a,filter_hash)
    d2 = sum.(d)
    d1 = eqws_compute(a,filter_hash)
    @test maximum(abs.(d1.-d2)) < 1e-15
end

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
function all_corners_mark(r::Robot)
    coordinates = []
    cornerwalk(r, coordinates)
    paint!(r)
    go_back!(r, coordinates)
end

function paint!(r::Robot)
    for side in (Nord, Ost, Sud, West)
        moves!(r, side)
        putmarker!(r)
    end
end

function cornerwalk(r::Robot, coordinates::Vector{Any})
    while !(isborder(r, West) && isborder(r, Sud))
        if isborder(r, West) == false
            move!(r, West)
            push!(coordinates, West)
        end
        if isborder(r, Sud) == false
            move!(r, Sud)
            push!(coordinates, Sud)
        end
    end
end

all_corners_mark(r)


text:
"""
n = DMatrix!(n)
n = DMatrix!(n, (:mode=>:symbolic))
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 1/15/2017, Last Modified: 12/06/2019 \n
--------------------------------------------------------------------------------------\n
"""
function DMatrix!(n::NLOpt{T}, kwargs...) where { T <: Number } #TODO make IMatrix and option

    kw = Dict(kwargs)

    mode = get(kw, :mode, :defaut)

    if mode == :defaut

        check_ts = maximum.(n.ocp.ts)

        if maximum(check_ts) < 10 * eps()
            error("""
                ts is full of zeros: make sure that you call create_intervals() first to calculate ts!
                NOTE: This may have occured because  (:finalTimeDV => true) and the final time dv is not working properly!!
                """)
        end

        n.ocp.DMatrix = [zeros((n.ocp.Nck[int]),(n.ocp.Nck[int]+1)) for int in 1:n.ocp.Ni];
        DM = [zeros((n.ocp.Nck[int]),(n.ocp.Nck[int])) for int in 1:n.ocp.Ni];
        n.ocp.IMatrix = [zeros((n.ocp.Nck[int]),(n.ocp.Nck[int]+1)) for int in 1:n.ocp.Ni];

        for int in 1:n.ocp.Ni

            D = polyDiff(n.ocp.ts[int]) # +1 is already appended onto ts
            n.ocp.IMatrix[int] = zeros( (n.ocp.Nck[int]), (n.ocp.Nck[int]+1) )
            n.ocp.DMatrix[int] = D[1:end-1, :] # [Nck] X [Nck+1]

            if n.s.ocp.integrationScheme == :lgrImplicit

                DM = n.ocp.DMatrix[int][:,2:end]   # [Nck] X [Nck]

                n.ocp.IMatrix[int] = inv(DM)       # I = inv(D[:,2:N_k+1])

            end

        end

    elseif mode == :symbolic # for validation only, too slow otherwise

        error("""
            Cannot precompile with SymPy
            so this fucntion was turned off for typical use!!
            -> do a (using SymPy) in NLOptControl.jl then remove this error message and rerun
            """)
#=
        Dsym = [Array{Any}(undef, n.ocp.Nck[int],n.ocp.Nck[int]+1) for int in 1:n.ocp.Ni];
        n.ocp.DMatrix = [Array{Any}(undef, n.ocp.Nck[int],n.ocp.Nck[int]+1) for int in 1:n.ocp.Ni]
        test = [Array{Any}(undef, n.ocp.Nck[int]+1) for int in 1:n.ocp.Ni];
        val = 1; # since this is always = 1 this funtion is useful for testing, without scaling the problem from [-1,1] this was useful becuase tf was a design variable
        tf = Sym("tf")
        createIntervals!(n, tf); # gives symbolic expression
        for int in 1:n.ocp.Ni
            for i in 1:n.ocp.Nck[int]+1
                test[int][i] =  n.ocp.ts[int][i](tf=>val)
            end
        end
        for int in 1:n.ocp.Ni
            for idx in 1:n.ocp.Nck[int]+1
                for j in 1:n.ocp.Nck[int]
                    f = lagrange_basis_poly(tf, test[int], idx)
                    Dsym[int][j,idx] = diff(f,tf) # symbolic differentiation --> slow but useful # TODO include this in test functions
                    n.ocp.DMatrix[int][j,idx] = Dsym[int][j,idx](tf=>test[int][j])
                end
            end
        end
=#

    end

    return nothing

end

"""
n = create_intervals(n)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 12/23/2017, Last Modified: 12/06/2019 \n
--------------------------------------------------------------------------------------\n
"""
function createIntervals!(n::NLOpt{T}) where { T <: Number }
    tm = range(-1 , 1 ; length = n.ocp.Ni + 1)  # create mesh points
    di = 2 / n.ocp.Ni                           # interval size
    # go through each mesh interval creating time intervals; map [tm[i-1],tm[i]] --> [-1,1]
    n.ocp.ts = [[scale_tau(n.ocp.tau[int],tm[int],tm[int+1]);di*int-1] for int in 1:n.ocp.Ni]
    n.ocp.ws = [scale_w(n.ocp.w[int],tm[int],tm[int+1]) for int in 1:n.ocp.Ni]
    return nothing
end

max_stars_repo_name: JuliaMPC/NLOptControl

max_issues_repo_name: JuliaMPC/NLOptControl

max_forks_repo_name: JuliaMPC/NLOptControl


text:
#

using MultithreadingBenchmarks

function work(i, v, n)
    out = v
    for i in 1:n
        mul = @eval (a) -> a*$v
        out = Base.invokelatest(mul, out)
    end
    return out
end

#@time work(1,2,100)
#
#using PProf
#PProf.clear(); @pprof @sync begin
#    Threads.@spawn work(1, 2, 100)
#    Threads.@spawn work(1, 2, 100)
#    Threads.@spawn work(1, 2, 100)
#    Threads.@spawn work(1, 2, 100)
#end

MultithreadingBenchmarks.measure_work(work)


text:
include("polyagamma_sampler.jl")
#update the posterior, append sample to chain and set the value as a sample from the posterior
function _mcmc_iter_pg(
    parameters::AbstractParameters,
    latents::Vector{<:AbstractLatent},
    responses_val::Vector{Float64},
    W_val::Vector{Float64};
    sampling = true,
    )
        parameters.posterior = __posterior(parameters, latents, responses_val, W_val)
        _chain_append_and_set_val!(parameters; sampling = sampling)
    return parameters::AbstractParameters
end

function _mcmc_iter_pg(
    parameters::AbstractParameters,
    latents::Vector{<:AbstractLatent},
    responses_val::Vector{Float64},
    W_val::SharedArrays.SharedVector{Float64};
    sampling = true,
    )
        parameters.posterior = __posterior(parameters, latents, responses_val, W_val)
        _chain_append_and_set_val!(parameters; sampling = sampling)
    return parameters::AbstractParameters
end

function mcmc_iter_pg!(
    item::AbstractItem,
    examinees::Vector{<:AbstractExaminee},
    responses::Vector{<:AbstractResponse},
    W::Vector{PolyaGammaSample};
    sampling = true,
    already_sorted = false,
)
    if !already_sorted
        sort!(examinees, by = e -> e.idx)
        sort!(responses, by = r -> r.examinee_idx)
    end
    update_posterior!(item, examinees, responses, W; already_sorted = true)
    _chain_append_and_set_val!(item.parameters; sampling = sampling)
    return nothing
end

function _calibrate_item_pg(
    parameters::AbstractParameters,
    latents::Vector{<:AbstractLatent},
    response_vals::Vector{Float64};
    mcmc_iter::Int64 = 4_000,
    sampling::Bool = true,
    )
    for iter = 1:mcmc_iter
        W = map( l -> _generate_w(parameters, l), latents)
        parameters = _mcmc_iter_pg(parameters, latents, response_vals, W; sampling = sampling)
    end
    return parameters::AbstractParameters
end

function calibrate_item_pg!(
    item::AbstractItem,
    examinees::Vector{<:AbstractExaminee},
    responses::Vector{Response};
    mcmc_iter::Int64 = 4_000,
    sampling::Bool = true,
    already_sorted::Bool = false
)

    if !already_sorted
        sort!(examinees, by = e -> e.idx)
        sort!(responses, by = r -> r.examinee_idx)
        already_sorted = true
    end
    empty_chain!(item)
    # chain = SharedArrays.SharedArray{Float64}(2,mcmc_iter)
    # # extract `mcmc_iter` samples from the polyagamma and from theta conditional posterior
    # Distributed.@sync Distributed.@distributed for iter in 1:mcmc_iter
    #     W = generate_w(item, examinees)
    #     chain[:,iter] = rand(posterior(item, examinees, responses, W))
    # end
    # item.parameters.chain = [c[:] for c in eachcol(chain)]
    # same as, but without val update (faster)
    parameters = _calibrate_item_pg(
        item.parameters,
        get_latents(examinees),
        map( r -> r.val, responses),
        mcmc_iter = mcmc_iter,
        sampling = sampling
        )
    update_estimate!(item; sampling = sampling)
    return nothing
end


text:
# Autogenerated wrapper script for libtirpc_jll for x86_64-linux-musl
export libtirpc

JLLWrappers.@generate_wrapper_header("libtirpc")
JLLWrappers.@declare_library_product(libtirpc, "libtirpc.so.3")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libtirpc,
        "lib/libtirpc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module ToyPkg

greet() = print("Hello World 0.2.0!")

end # module


text:
module SemsFDTD

export SimData
export approx_stepsize
export new_sim
export timestep

"""
general notes

in 2d and 3d, want to operate as close to courant limit as possible

"""

"""
    SimData

    Struct containing the matrices required to run the FDTD simulation
"""
struct SimData{T<:AbstractFloat}
    e_x::AbstractArray{T,3}
    e_y::AbstractArray{T,3}
    e_z::AbstractArray{T,3}
    # H field definitions
    h_x::AbstractArray{T,3}
    h_y::AbstractArray{T,3}
    h_z::AbstractArray{T,3}
    # E field coefficients
    coeff_e_xe::AbstractArray{T,3}
    coeff_e_xh::AbstractArray{T,3}
    coeff_e_ye::AbstractArray{T,3}
    coeff_e_yh::AbstractArray{T,3}
    coeff_e_ze::AbstractArray{T,3}
    coeff_e_zh::AbstractArray{T,3}
    # H field coefficients
    coeff_h_xe::AbstractArray{T,3}
    coeff_h_xh::AbstractArray{T,3}
    coeff_h_ye::AbstractArray{T,3}
    coeff_h_yh::AbstractArray{T,3}
    coeff_h_ze::AbstractArray{T,3}
    coeff_h_zh::AbstractArray{T,3}
    # ABC X axis Planes
    abc_xp0_ty::AbstractArray{T,2}
    abc_xp0_tz::AbstractArray{T,2}
    abc_xp1_ty::AbstractArray{T,2}
    abc_xp1_tz::AbstractArray{T,2}
    # ABC Y axis Planes
    abc_yp0_tx::AbstractArray{T,2}
    abc_yp0_tz::AbstractArray{T,2}
    abc_yp1_tx::AbstractArray{T,2}
    abc_yp1_tz::AbstractArray{T,2}
    # ABC Z axis planes
    abc_zp0_tx::AbstractArray{T,2}
    abc_zp0_ty::AbstractArray{T,2}
    abc_zp1_tx::AbstractArray{T,2}
    abc_zp1_ty::AbstractArray{T,2}
end

# freespace impedance
vacuumZ = 376.730314 # ohms

# Courant limit in 3d space is 1/sqrt(3)
courant_limit = 1.0 / sqrt( 3.0 );

# first order ABC coefficient
abc_coeff = ( courant_limit - 1.0 ) / ( courant_limit + 1.0 )

# speed of light
c = 299792458

"""
    approx_stepsize - Give the approximate discrete spatial and time step sizes
    for a given minimum wavelength
"""
function approx_stepsize( min_wavelength::Float64 )
    delta_s = min_wavelength / 20
    delta_t = ( delta_s * courant_limit ) / c;
    return ( delta_s, delta_t )
end

"""
    new_sim - construct and initialize a new simulation assuming
"""
function new_sim( fprecision::DataType, dim_x::Int, dim_y::Int, dim_z::Int )::SimData
    newsim = SimData( zeros( fprecision, dim_x - 1, dim_y, dim_z ), # e_x
                      zeros( fprecision, dim_x, dim_y - 1, dim_z ), # e_y
                      zeros( fprecision, dim_x, dim_y, dim_z - 1 ), # e_z
                      zeros( fprecision, dim_x, dim_y - 1, dim_z - 1 ), # h_x
                      zeros( fprecision, dim_x - 1, dim_y, dim_z - 1 ), # h_y
                      zeros( fprecision, dim_x - 1, dim_y - 1, dim_z ), # h_z
                      zeros( fprecision, dim_x - 1, dim_y, dim_z ), # coeff_e_xe
                      zeros( fprecision, dim_x - 1, dim_y, dim_z ), # coeff_e_xh
                      zeros( fprecision, dim_x, dim_y - 1, dim_z ), # coeff_e_ye
                      zeros( fprecision, dim_x, dim_y - 1, dim_z ), # coeff_e_yh
                      zeros( fprecision, dim_x, dim_y, dim_z - 1 ), # coeff_e_ze
                      zeros( fprecision, dim_x, dim_y, dim_z - 1 ), # coeff_e_zh
                      zeros( fprecision, dim_x, dim_y - 1, dim_z - 1 ), # coeff_h_xe
                      zeros( fprecision, dim_x, dim_y - 1, dim_z - 1 ), # coeff_h_xh
                      zeros( fprecision, dim_x - 1, dim_y, dim_z - 1 ), # coeff_h_ye
                      zeros( fprecision, dim_x - 1, dim_y, dim_z - 1 ), # coeff_h_yh
                      zeros( fprecision, dim_x - 1, dim_y - 1, dim_z ), # coeff_h_ze
                      zeros( fprecision, dim_x - 1, dim_y - 1, dim_z ), # coeff_h_zh
                      zeros( fprecision, dim_y - 1, dim_z ), # abc_xp0_ty
                      zeros( fprecision, dim_y, dim_z - 1 ), # abc_xp0_tz
                      zeros( fprecision, dim_y - 1, dim_z ), # abc_xp1_ty
                      zeros( fprecision, dim_y, dim_z - 1 ), # abc_xp1_tz
                      zeros( fprecision, dim_x - 1, dim_z ), # abc_yp0_tx
                      zeros( fprecision, dim_x, dim_z - 1 ), # abc_yp0_tz
                      zeros( fprecision, dim_x - 1, dim_z ), # abc_yp1_tx
                      zeros( fprecision, dim_x, dim_z - 1 ), # abc_yp1_tz
                      zeros( fprecision, dim_x - 1, dim_y ), # abc_zp0_tx
                      zeros( fprecision, dim_x, dim_y - 1 ), # abc_zp0_ty
                      zeros( fprecision, dim_x - 1, dim_y ), # abc_zp1_tx
                      zeros( fprecision, dim_x, dim_y - 1 ) ) # abc_zp1_ty

    # initialize e-field coefficients to free space for now
    newsim.coeff_e_xe .= 1.0
    newsim.coeff_e_xh .= courant_limit * vacuumZ

    return newsim
end

"""
    update_e_field -
"""
function update_e_field( sim::SimData )
    # x component
    sim.e_x[ :, begin + 1 : end - 1, begin + 1 : end - 1 ] =
        ( sim.coeff_e_xe[ :, begin + 1 : end - 1, begin + 1 : end - 1 ] .*
          sim.e_x[ :, begin + 1 : end - 1, begin + 1 : end - 1 ] ) .+
        sim.coeff_e_xh[ :, begin + 1 : end - 1, begin + 1 : end - 1 ] .*
        ( ( sim.h_z[ :, begin + 1 : end, begin + 1 : end - 1 ]
         .- sim.h_z[ :, begin : end - 1, begin + 1 : end - 1 ] )
        .- ( sim.h_y[ :, begin + 1 : end - 1, begin + 1 : end ]
          .- sim.h_y[ :, begin + 1 : end - 1, begin : end - 1 ] ) )
    # y component
    sim.e_y[ begin + 1 : end - 1, :, begin + 1 : end - 1 ] =
        ( sim.coeff_e_ye[ begin + 1 : end - 1, :, begin + 1 : end - 1 ] .*
          sim.e_y[ begin + 1 : end - 1, :, begin + 1 : end - 1 ] ) .+
        sim.coeff_e_yh[ begin + 1 : end - 1, :, begin + 1 : end - 1 ] .*
        ( ( sim.h_x[ begin + 1 : end - 1, :, begin + 1 : end ]
         .- sim.h_x[ begin + 1 : end - 1, :, begin : end - 1 ] )
        .- ( sim.h_z[ begin + 1 : end, :, begin + 1 : end - 1 ]
          .- sim.h_z[ begin : end - 1, :, begin + 1 : end - 1 ] ) )
    # z component
    sim.e_z[ begin + 1 : end - 1, begin + 1 : end - 1, : ] =
        ( sim.coeff_e_ze[ begin + 1 : end - 1, begin + 1 : end - 1, : ] .*
          sim.e_z[ begin + 1 : end - 1, begin + 1 : end - 1, : ] ) .+
        sim.coeff_e_zh[  begin + 1 : end - 1, begin + 1 : end - 1, : ] .*
        ( ( sim.h_y[ begin + 1 : end, begin + 1 : end - 1, : ]
         .- sim.h_y[ begin : end - 1, begin + 1 : end - 1, : ] )
        .- ( sim.h_x[ begin + 1 : end - 1, begin + 1 : end, : ]
          .- sim.h_x[ begin + 1 : end - 1, begin : end - 1, : ] ) )

    return nothing
end

"""
    update_h_field -
"""
function update_h_field( sim::SimData )
    # x component
    sim.h_x[:] = ( sim.coeff_h_xh .* sim.h_x ) .+ sim.coeff_h_xe .*
        ( ( sim.e_y[ :, :, begin + 1 : end ] .- sim.e_y[ :, :, begin : end - 1 ] )
        .- ( sim.e_z[ :, begin + 1 : end, : ] .- sim.e_z[ :, begin : end - 1, : ] ) )
    # y component
    sim.h_y[:] = ( sim.coeff_h_yh .* sim.h_y ) .+ sim.coeff_h_ye .*
        ( ( sim.e_z[ begin + 1 : end, :, : ] .- sim.e_z[ begin : end - 1, :, : ] )
        .- ( sim.e_x[ :, :, begin + 1 : end ] .- sim.e_x[ :, :, begin : end - 1 ] ) )
    # z component
    sim.h_z[:] = ( sim.coeff_h_zh .* sim.h_z ) .+ sim.coeff_h_ze .*
        ( ( sim.e_x[ :, begin + 1 : end, : ] .- sim.e_x[ :, begin : end - 1, : ] )
        .- ( sim.e_y[ begin + 1 : end, :, : ] .- sim.e_y[ begin : end - 1, :, : ] ) )

    return nothing
end

"""
    enforce_abc - enforce first order ABC at edges of box
"""
function enforce_abc( sim::SimData )
    # x plane 0
    sim.e_y[ begin, :, : ] = sim.abc_xp0_ty .+ abc_coeff * ( sim.e_y[ begin + 1, :, : ] - sim.e_y[ begin, :, : ] )
    sim.abc_xp0_ty[:] = sim.e_y[ begin + 1, :, : ]
    sim.e_z[ begin, :, : ] = sim.abc_xp0_tz .+ abc_coeff * ( sim.e_z[ begin + 1, :, : ] - sim.e_z[ begin, :, : ] )
    sim.abc_xp0_tz[:] = sim.e_z[ begin + 1, :, : ]
    # x plane 1
    sim.e_y[ end, :, : ] = sim.abc_xp1_ty .+ abc_coeff * ( sim.e_y[ end - 1, :, : ] - sim.e_y[ end, :, : ] )
    sim.abc_xp1_ty[:] = sim.e_y[ end - 1, :, : ]
    sim.e_z[ end, :, : ] = sim.abc_xp1_tz .+ abc_coeff * ( sim.e_z[ end - 1, :, : ] - sim.e_z[ end, :, : ] )
    sim.abc_xp1_tz[:] = sim.e_z[ end - 1, :, : ]
    # y plane 0
    sim.e_x[ :, begin, : ] = sim.abc_yp0_tx .+ abc_coeff * ( sim.e_x[ :, begin + 1, : ] - sim.e_x[ :, begin, : ] )
    sim.abc_yp0_tx[:] = sim.e_x[ :, begin + 1, : ]
    sim.e_z[ :, begin, : ] = sim.abc_yp0_tz .+ abc_coeff * ( sim.e_z[ :, begin + 1, : ] - sim.e_z[ :, begin, : ] )
    sim.abc_yp0_tz[:] = sim.e_z[ :, begin + 1, : ]
    # y plane 1
    sim.e_x[ :, end, : ] = sim.abc_yp1_tx .+ abc_coeff * ( sim.e_x[ :, end - 1, : ] - sim.e_x[ :, end, : ] )
    sim.abc_yp1_tx[:] = sim.e_x[ :, end - 1, : ]
    sim.e_z[ :, end, : ] = sim.abc_yp1_tz .+ abc_coeff * ( sim.e_z[ :, end - 1, : ] - sim.e_z[ :, end, : ] )
    sim.abc_yp1_tz[:] = sim.e_z[ :, end - 1, : ]
    # z plane 0
    sim.e_x[ :, :, begin ] = sim.abc_zp0_tx .+ abc_coeff * ( sim.e_x[ :, :, begin + 1 ] - sim.e_x[ :, :, begin ] )
    sim.abc_zp0_tx[:] = sim.e_x[ :, :, begin + 1 ]
    sim.e_y[ :, :, begin ] = sim.abc_zp0_ty .+ abc_coeff * ( sim.e_y[ :, :, begin + 1 ] - sim.e_y[ :, :, begin ] )
    sim.abc_zp0_ty[:] = sim.e_y[ :, :, begin + 1 ]
    # z plane 1
    sim.e_x[ :, :, end ] = sim.abc_zp1_tx .+ abc_coeff * ( sim.e_x[ :, :, end - 1 ] - sim.e_x[ :, :, end ] )
    sim.abc_zp1_tx[:] = sim.e_x[ :, :, end - 1 ]
    sim.e_y[ :, :, end ] = sim.abc_zp1_ty .+ abc_coeff * ( sim.e_y[ :, :, end - 1 ] - sim.e_y[ :, :, end ] )
    sim.abc_zp1_ty[:] = sim.e_y[ :, :, end - 1 ]

    return nothing
end

"""
    timestep - advance the simulation one timestep
"""
function timestep( sim::SimData )
    # enforce_pec()

    update_h_field( sim )
    update_e_field( sim )

    # handle any field dumping here

    # handle any sources here - TODO: before/after abc?

    enforce_abc( sim )

    return nothing
end

"""
    enforce_pec - apply zero electric field (PEC) to specified regions
        region argument expected in form of [ [ 1:2, 2:3, 3:4 ], [ 5:6, 7:8, 9:10 ] ]
"""
function enforce_pec( sim::SimData, regions::Vector{Vector{UnitRange{Int64}}} )
    for region in regions
        sim.e_x( region ) .= 0.0
        sim.e_y( region ) .= 0.0
        sim.e_z( region ) .= 0.0
    end

    return nothing
end

"""
    apply_source_node -
        TODO: does this make sense as a function
"""
function apply_source_node()
    return nothing
end

"""
    apply_array_nodes -
        TODO: does this make sense as a function
"""
function apply_array_nodes()
    return nothing
end

"""
    gaussian_src -

"""
function gaussian_src( t::Float64, mean::Float64, σ::Float64 )
    return exp( -( t - mean )^2 / ( 2 * σ^2 ) ) / ( sqrt( 2 * π ) * σ )
end

"""
    gaussian_energy -
"""
function gaussian_energy( σ::Float64 )
    return 1 / ( 2 * σ * sqrt( π ) )
end

"""
    gaussian_derivative_src -
"""
function guassian_derivative( t::Float64, mean::Float64, σ::Float64 )
    return -( t - mean ) * exp( -( t - mean )^2 / ( 2 * σ^2 ) ) / ( σ^3 * sqrt( 2 * π ) )
end

"""
    gaussian_derivative_energy -
"""
function gaussian_derivative_energy( σ::Float64 )
    return 1 / ( 4 * σ^3 * sqrt( π ) )
end

"""
"""
function quartic_poly_src( t::Float64, duration::Float64 )
    τ = 1 - 2 * ( t / duration )
    if abs( τ ) > 1
        return 0
    else
        return ( 1 - τ^2 )^4
    end
end

function quartic_poly_derivative( t::Float64, duration::Float64 )
    τ = 1 - 2 * ( t / duration )
    if abs( τ ) > 1
        return 0
    else
        return -8 * τ * ( 1 - t^2 )^3
    end
end

function quartic_poly_sec_deriv( t::Float64, duration::Float64 )
    τ = 1 - 2 * ( t / duration )
    if abs( τ ) > 1
        return 0
    else
        return 48 * τ^2 * ( 1 - τ^2 )^2 - 8 * ( 1 - τ^2 )^3
    end
end

end # module

max_stars_repo_name: asandwhich/SemsFDTD

max_issues_repo_name: asandwhich/SemsFDTD

max_forks_repo_name: asandwhich/SemsFDTD


text:
import LevelJump

function test_basic()
    path = mktempdir()
    ldb = LevelJump.LevelDB(path)
    LevelJump.push!(ldb, "cheese", "yes")
end

test_basic()


text:
using OpenQuantumTools, Test, OrdinaryDiffEq, QuadGK

H = DenseHamiltonian([(s) -> 1.0], [σi], unit=:ħ)
u0 = PauliVec[1][1]
coupling = ConstantCouplings(["Z"], unit=:ħ)

# test for OhmicBath
bath = Ohmic(1e-4, 4, 16)
annealing = Annealing(H, u0; coupling=coupling, bath=bath)
tf = 10
U = solve_unitary(annealing, tf, alg=Tsit5(), abstol=1e-8, reltol=1e-8)
sol = solve_redfield(annealing, tf, InplaceUnitary(U), alg=Tsit5(),
    reltol=1e-6)
cfun(t) = correlation(t, bath)
f(t) = quadgk(cfun, 0, t)[1]
γ = real(quadgk(f, 0, tf)[1])
@test sol(10)[1, 2] ≈ exp(-4 * γ) * 0.5 atol = 1e-5 rtol = 1e-5

sol = solve_redfield(annealing, tf, InplaceUnitary(U), vectorize=true,
    alg=TRBDF2(), reltol=1e-6)
@test sol(10)[2] ≈ exp(-4 * γ) * 0.5 atol = 1e-5 rtol = 1e-5

# test for CustomBath
H = DenseHamiltonian([(s) -> 0.0], [σi], unit=:ħ)
u0 = PauliVec[1][1]
coupling = ConstantCouplings(["Z"], unit=:ħ)
tf = 20
cfun(x) = x <= 20 ? 1e-4 : 0
bath = CustomBath(correlation=cfun)
annealing = Annealing(H, u0; coupling=coupling, bath=bath)
U = solve_unitary(annealing, tf, alg=Tsit5(), abstol=1e-8, reltol=1e-8);
sol = solve_redfield(
    annealing,
    tf,
    InplaceUnitary(U);
    alg=Tsit5(),
    abstol=1e-8,
    reltol=1e-8,
);
@test sol(tf)[1, 2] ≈ exp(-4 * 0.02) * 0.5 atol = 1e-4 rtol = 1e-4

# test for InstPulseCallback
cbu = InstPulseCallback([0.5 * tf], (c, x) -> c .= σx * c)
U = solve_unitary(
    annealing,
    tf,
    alg=Tsit5(),
    abstol=1e-8,
    reltol=1e-8,
    callback=cbu,
);
cb = InstPulseCallback([0.5 * tf], (c, x) -> c .= σx * c * σx)
sol = solve_redfield(
    annealing,
    tf,
    InplaceUnitary(U),
    alg=Tsit5(),
    abstol=1e-8,
    reltol=1e-8,
    callback=cb,
)
@test sol(tf)[1, 2] ≈ 0.5 atol = 1e-5 rtol = 1e-5

# test for CustomCouplings
H = DenseHamiltonian([(s) -> 1.0], [σz], unit=:ħ)
Uz(t) = exp(-1.0im * t * σz)
u0 = PauliVec[3][1]
coupling = ConstantCouplings(["X"], unit=:ħ)
bath = Ohmic(1e-4, 4, 16)
annealing = Annealing(H, u0; coupling=coupling, bath=bath)
tf = 10
sol = solve_redfield(annealing, tf, Uz, alg=Tsit5(), reltol=1e-6)

iH = DenseHamiltonian([(s) -> 0.0], [σz], unit=:ħ)
icoupling = CustomCouplings([(s) -> Uz(10 * s)' * σx * Uz(10 * s)], unit=:ħ)
iannealing = Annealing(iH, u0; coupling=icoupling, bath=bath)
isol = solve_redfield(iannealing, tf, (t) -> σi, alg=Tsit5(), reltol=1e-6)

@test sol[end][1,1] ≈ isol[end][1,1] atol = 1e-6 rtol = 1e-6

# Test of non-positivity check
β = 4
T = β_2_temperature(β)
η = 0.1
fc= 10/(2π)
bath = Ohmic(η, fc, T)

Hp = 0.5*σz⊗σi - 0.7*σi⊗σz + 0.3*σz⊗σz
Hd = standard_driver(2)
H = DenseHamiltonian([(s)->1-s, (s)->s], [-Hd, Hp], unit=:ħ)

tf = 20
ρ0 = (σi+σx)⊗(σi+σx)/4
coupling = ConstantCouplings([σz⊗σi, σi⊗σz], unit=:ħ)
annealing = Annealing(H, ρ0, bath=bath, coupling=coupling)

U = solve_unitary(annealing, tf, alg = Tsit5(), abstol=1e-7, reltol=1e-7)

redfield_sol = solve_redfield(annealing, tf, U, alg = Tsit5(), abstol=1e-7, reltol=1e-7, callback=PositivityCheckCallback())

@test redfield_sol.t[end] ≈ 3.45278167 atol = 1e-5 rtol = 1e-5

max_stars_repo_name: huangjingzheng/OpenQuantumTools

max_issues_repo_name: huangjingzheng/OpenQuantumTools

max_forks_repo_name: huangjingzheng/OpenQuantumTools


text:
module Power

using PowerSystems
using JuMP, Gurobi
using Random
using Distributions

export pᵈᵗAllTime, qᵈᵗAllTime, solve_LMP, bus, branch

Random.seed!(1234)
##########################################
# case22 info
##########################################
DATA_DIR = raw"./data/case22.m"
system_data = System(DATA_DIR)
# sorted bus and branch
bus = get_components(Bus, system_data) |> collect;
sort!(bus, by = x->x.number); # ascending order
branch = get_components(Branch, system_data) |> collect;
sort!(branch, by = branch->branch.arc.from.number);
# for each node, find the idx of the subsequent edge
π = [findall(x -> x.arc.from.number == i, branch)  for i in 1:length(bus)]
# voltage limit
Umin = [bus_i.voltage_limits.min^2 for bus_i in bus];
Umax = [bus_i.voltage_limits.max^2 for bus_i in bus];

##########################################
# case22 additional configuration
##########################################
# time slots
T = 24
# standard load ( mean = 1.0 )
standard_load = [0.80092226, 0.75151295, 0.72181773, 0.70716663, 0.7085614,
                0.72807473, 0.81076389, 0.90212761, 1.02712352, 1.12649447,
                1.17364172, 1.1766259 , 1.08926146, 1.08708576, 1.09357851,
                1.10167761, 1.13332739, 1.15798613, 1.15985351, 1.21728262,
                1.21873949, 1.16554919, 1.03838863, 0.90243688];
# MVA
# In 2019, per-capita average energy consumption is 4.905 MWh in China
# 100,000 people city power load： 10^5 * 4.905 / 365 / 24 = 56 MW
# power factor = 0.95, Q : P ≈ 0.3
bus_average_P_load = 10 * [ 0.000,  0.168,  0.168,  0.338,  0.146,  0.105,  0.088,  0.144,
                 0.193,  0.144,  0.163,  0.163,  0.821,  0.347,  0.347,  0.803,
                 0.496,  0.496,  0.438,  0.373,  0.373,  0.310];
bus_average_Q_load = 3 * [ 0.000,  0.209,  0.209,  0.373,  0.125,  0.142,  0.117,  0.186,
                 0.259,  0.186,  0.195,  0.195,  0.717,  0.301,  0.301,  0.701,
                 0.478,  0.478,  0.389,  0.360,  0.360,  0.294 ];
# traditional loads, (bus, T)
pᵈᵗAllTime =  (bus_average_P_load * standard_load') .* rand(Normal(1, 0.2), (length(bus), T))
qᵈᵗAllTime =  (bus_average_Q_load * standard_load') .* rand(Normal(1, 0.2), (length(bus), T))
# Generator limit / MVA
pᵍmin = zeros(length(bus))
pᵍmax = ones(length(bus)) * 20
pᵍmax[1] = Inf # infinite input
qᵍmin = ones(length(bus)) * -5
qᵍmax = ones(length(bus)) * 5
qᵍmax[1] = Inf # infinite input
# Generator cost coefficient
a = 0.3 * ones(length(bus))
b = 150 * ones(length(bus))
# buying price of power market $ / MVA
ρₘ = 140
# branch capacity limit / MVA
Sˡ = ones(length(branch)) * 30;

"""
- @description: solve an LMP problem
- @param {pᵈᵗ} traditional active power demand for each bus
- @param {qᵈᵗ} traditional reactive power demand for each bus
- @param {pᵈʰ} hydrogen electrolyser active power demand for each bus
- @param {qᵈʰ} hydrogen electrolyser reactive power demand for each bus
- @return {λ}  LMP price for each bus
"""
function solve_LMP(pᵈᵗ, qᵈᵗ, pᵈʰ,  qᵈʰ)
    # power load
    pᵈ = pᵈᵗ + pᵈʰ
    qᵈ = qᵈᵗ + qᵈʰ
    # optimization problem
    model = Model(Gurobi.Optimizer)
    set_silent(model)
    JuMP.set_optimizer_attribute(model, "QCPDual", 1)
    # decision variables
    @variable(model, pᵍ[1:length(bus)]); #  generator active power for each bus
    @variable(model, qᵍ[1:length(bus)]); # generator reactive power for each bus
    @variable(model, U[1:length(bus)] >= 0); # squared voltage for each bus
    @variable(model, I[1:length(branch)] >= 0); # squared current for each branch
    @variable(model, Pˡ[1:length(branch)] >= 0); # active power of branch
    @variable(model, Qˡ[1:length(branch)] >= 0); # reactive power of branch
    # constraints
    @constraint(model, balance[k=1:length(branch)], Pˡ[k] + pᵍ[branch[k].arc.to.number] - branch[k].r * I[k] == sum(Pˡ[π[branch[k].arc.to.number]]) + pᵈ[branch[k].arc.to.number] )
    @constraint(model, [k=1:length(branch)], Qˡ[k] + qᵍ[branch[k].arc.to.number] - branch[k].x * I[k] == sum(Qˡ[π[branch[k].arc.to.number]]) + qᵈ[branch[k].arc.to.number] )
    @constraint(model, [k=1:length(branch)], U[branch[k].arc.to.number] == U[branch[k].arc.from.number] - 2 * (branch[k].r * Pˡ[k] + branch[k].x * Qˡ[k]) + (branch[k].r^2 + branch[k].x^2) * I[k] )
    @constraint(model, [k=1:length(branch)], [I[k] + U[branch[k].arc.from.number], 2Pˡ[k], 2Qˡ[k], I[k] - U[branch[k].arc.from.number]] ∈ SecondOrderCone())
    @constraint(model, pᵍmin .<= pᵍ .<= pᵍmax)
    @constraint(model, qᵍmin .<= qᵍ .<= qᵍmax)
    @constraint(model, Umin .<= U .<= Umax)
    @constraint(model, [k=1:length(branch)], [Sˡ[k], Pˡ[k], Qˡ[k]] ∈ SecondOrderCone())
    @constraint(model, [k=1:length(branch)], Pˡ[k] - branch[k].r * I[k] >= 0 )
    @constraint(model, [k=1:length(branch)], Qˡ[k] - branch[k].x * I[k] >= 0 )
    # objective
    @objective(model, Min, sum(a .* pᵍ .^ 2 + b .* pᵍ) + ρₘ * sum(Pˡ[π[1]]) )

    optimize!(model)
    println("【Power Result】:", raw_status(model))
    # if termination_status(model) == MOI.OPTIMAL
    #     println("Solution is optimal")
    # elseif termination_status(model) == MOI.TIME_LIMIT && has_values(model)
    #     println("Solution is suboptimal due to a time limit, but a primal solution is available")
    #     println(pᵈʰ)
    # else
    #     println("The model was not solved correctly.")
    #     println(pᵈʰ)
    # end

    return dual.(balance)
end

end


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

module Meshes

using IterTools
using StaticArrays
using CircularArrays
using SpecialFunctions
using LinearAlgebra
using RecipesBase

import Base: ==, +, -
import Random

# numerical tolerances
include("tolerances.jl")

# basic types
include("vectors.jl")
include("points.jl")
include("angles.jl")

# geometries and meshes
include("geometries.jl")
include("connectivities.jl")
include("conventions.jl")
include("mesh.jl")

# algorithms
include("sampling.jl")
include("discretization.jl")
include("boundingboxes.jl")

# utilities
include("utils.jl")

# plot recipes
include("plotrecipes/polytopes.jl")
include("plotrecipes/chains.jl")
include("plotrecipes/polyareas.jl")

export
  # points
  Point, Point1, Point2, Point3, Point1f, Point2f, Point3f,
  embeddim, coordtype, coordinates,

  # vectors
  Vec, Vec1, Vec2, Vec3, Vec1f, Vec2f, Vec3f,

  # angles
  ∠,

  # geometries
  Geometry,
  GeometrySet,
  embeddim, paramdim, coordtype,
  measure, area, volume, boundary,

  # primitives
  Primitive,
  Box, Ball, Sphere, Cylinder,
  center, radius, height, sides,

  # polytopes
  Polytope, Polygon, Polyhedron,
  Segment, Triangle, Quadrangle,
  Pyramid, Tetrahedron, Hexahedron,
  Chain, PolyArea,
  vertices, nvertices,
  faces, facets,
  intersecttype,
  windingnumber, chains,
  isclosed, issimple, hasholes,
  angles, innerangles, close!, open!,
  orientation, bridge,

  # connectivities
  Connectivity,
  polytopetype, connect, materialize,

  # conventions
  OrderingConvention,
  GMSH,
  connectivities,

  # meshes
  Mesh,
  CartesianGrid, UnstructuredMesh,
  faces, elements, spacing,

  # sampling
  SamplingMethod,
  RegularSampling,
  sample,

  # discretization
  DiscretizationMethod,
  FIST,
  discretize,

  # bounding boxes
  boundingbox,

  # utililities
  signarea,
  sideof

end # module


text:
using BayesianDataFusion
using Base.Test
using Compat

addprocs(2)
@everywhere using BayesianDataFusion

rows = Int32[ 1:200; 151:350 ]
cols = Int32[ 151:350; 1:2:399 ]

A = ParallelSBM(rows, cols, workers()[1:2])

########        copyto test         ########
addprocs(2)
@everywhere using BayesianDataFusion
A2 = BayesianDataFusion.copyto(BayesianDataFusion.nonshared(A), workers()[end-1:end])
@test A.m == A2.m
@test A.n == A2.n
@test A.numblocks == A2.numblocks
@test A2.pids     == workers()[end-1:end]
@test length(A.tmp)   == length(A2.tmp)
@test length(A.sems)  == length(A2.sems)
@test length(A.sbms)  == length(A2.sbms)
@test length(A.logic) == length(A2.logic)

########    solving cg after copy   #######
x    = rand(size(A2,2))
beta = BayesianDataFusion.cg_AtA(A2, x, 0.5)
Asp  = sparse(rows, cols, 1.0)
AA   = full(Asp' * Asp)
beta_e = (AA + eye(length(x))*0.5) \ x
@test_approx_eq beta beta_e

########     remote cg     #########
x    = rand(size(A,2))
Ans  = BayesianDataFusion.nonshared(A)
Fref = @spawnat 2 BayesianDataFusion.copyto(Ans, Int[3,4])
beta2 = fetch(@spawnat 2 BayesianDataFusion.cg_AtA_ref(Fref, x, 0.75, 1e-6, length(x)))
Asp   = sparse(rows, cols, 1.0)
AA    = full(Asp' * Asp)
beta2_e = (AA + eye(length(x))*0.75) \ x
@test_approx_eq beta2 beta2_e

########    solve_cg2 test     ########
rhs = rand(size(A,2), 3)
Y   = BayesianDataFusion.solve_cg2(RemoteRef[Fref], rhs, 0.5)
Ye  = (AA + eye(size(rhs,1))*0.5) \ rhs
@test_approx_eq Y Ye

########  Frefs_mul_B on SparseMatrixCSC  #########
Xm = rand(size(Asp,2), 5)
Asp_ref1 = @spawnat 2 fetch(Asp)
Asp_ref2 = @spawnat 3 fetch(Asp)
Ym  = BayesianDataFusion.Frefs_mul_B(RemoteRef[Asp_ref1, Asp_ref2], Xm)
Yme = Asp * Xm
@test_approx_eq Ym Yme

########  Frefs_mul_B on ParallelSBM  #########
Ym_psbm = BayesianDataFusion.Frefs_mul_B(RemoteRef[Fref], Xm)
@test_approx_eq Ym_psbm Yme

########  Frefs_t_mul_B on SparseMatrixCSC  #########
Xt  = rand(size(Asp,1), 5)
Yt  = BayesianDataFusion.Frefs_t_mul_B(RemoteRef[Fref], Xt)
Yte = At_mul_B(Asp, Xt)
@test_approx_eq Yt Yte

########  macau with several SparseMatrixCSC  #######
Y    = sprand(20, 10, 0.2)
feat = sprand(20, 5, 0.5)
rd   = RelationData(Y, class_cut = 0.5, feat1 = feat)
assignToTest!(rd.relations[1], 2)

res = macau(rd, burnin = 2, psamples = 2, verbose = false, compute_ff_size=0, cg_pids=Int[2,3,4], num_latent=5)

########  macau with ParallelSBM  ########
Y     = sprand(20, 10, 0.4)
rows, cols, vals = findnz(Y)
A     = ParallelSBM(rows, cols, workers()[1:2])
rd2   = RelationData(Y, class_cut = 0.5, feat1 = A)
assignToTest!(rd.relations[1], 2)

res2  = macau(rd2, burnin = 2, psamples = 2, verbose = false, compute_ff_size=0, cg_pids=Vector{Int}[ [2, 2,3], [4, 4,5] ], num_latent=5)

## using Dictionary for cg_pids
rd2   = RelationData(Y, class_cut = 0.5, feat1 = A)
assignToTest!(rd.relations[1], 2)
res3  = macau(rd2, burnin = 2, psamples = 2, verbose = false, compute_ff_size=0, cg_pids=@compat Dict(2 => [2,3], 4 => [4,5]) )