text:
using Test
using Printf

tests = Test.DefaultTestSet("all")
Test.push_testset(tests)

dir = dirname(@__FILE__)
for f in readdir(dir)
    m = match(r"^test_(.*)\.jl$", f)
    fp = joinpath(dir, f)
    if m === nothing || !isfile(fp)
        continue
    end

    name = m.captures[1]
    @printf("%-12s ... ", name)

    ts = Test.DefaultTestSet(name)
    Test.push_testset(ts)
    try
        evalfile(fp)
    catch err
        res = Test.Error(:nontest_error, :(), err, catch_backtrace())
        Test.record(ts, res)
    end
    Test.pop_testset()
    Test.finish(ts)

    # lazy: use get_test_counts to set ts.anynonpass
    Test.get_test_counts(ts)
    ts.anynonpass || println("OK")
end

Test.pop_testset()
Test.finish(tests)


text:
#
# Just for testing Julia and C++/cocos2dx
#
# @author DeadbraiN
#
module RemoteWorldJson
  import Event
  import Connection
  import Creature
  import Client
  import RpcApi
  import FastApi
  import Config
  import Helper
  import Dots
  import JSON

  export create
  export start
  export stop

  export RemoteDataRT
  #
  # Amount of records (objects) in one JSON file
  #
  const AMOUNT_OF_RECORDS = 10000
  #
  # Name of the folder for JSON files
  #
  const JSON_FOLDER = "jsons"
  #
  # Contains data of for remote host, from where we displaying
  # world's region and shows it on a canvas.
  #
  type RemoteDataRT
    #
    # Required fields
    #
    cmdCon::Client.ClientConnection
    poolingCon::Client.ClientConnection
    width::UInt16
    height::UInt16
    fileIndex::Int
    diffs::Array{Any,1}
    #
    # Lazy loading fields
    #
    poolingBeforeCb::Function
    cmdAfterCb::Function
    ts::Float64
    poolingRequests::Int
    oldRequests::Int
    ips::Float64

    RemoteDataRT(
      cmdCon::Client.ClientConnection,
      poolingCon::Client.ClientConnection,
      width::UInt16,
      height::UInt16,
      fileIndex::Int,
      diffs::Array{Any, 1}
    ) = new(cmdCon, poolingCon, width, height, fileIndex, diffs)
  end
  #
  # Creates connection with remote host for display pixels from remote world.
  # @param cfg Global config type
  # @param host Remote host we are connecting to
  # @param cmdPort Remote port for commands
  # @param poolingPort Remote port for pooling requests
  # @return {RemoteDataRT}
  #
  function create(cfg::Config.ConfigData, host::Base.IPAddr, cmdPort::Int, poolingPort::Int)
    if !isdir(JSON_FOLDER) mkdir(JSON_FOLDER) end

    RemoteDataRT(
      Client.create(host, cmdPort),
      Client.create(host, poolingPort, true),
      UInt16(cfg.worldWidth),
      UInt16(cfg.worldHeight),
      0,
      Any[]
    )
  end
  #
  # Start displaying remote plane. It handles requests from remote
  # server with dots diffs. Response of first request will be all dots
  # of remote plane. All other requests will be runned from remove
  # server back to this client.
  # @param rd Remote world data object. See create()
  #
  function start(rd::RemoteDataRT)
    local color::Dots.Color

    rd.poolingBeforeCb = (sock::Base.TCPSocket, data::Array{Any, 1}, ans::Connection.Answer) -> _onDot(rd, data)
    rd.cmdAfterCb      = (ans::Connection.Answer) -> _onRegion(rd, ans)
    rd.ts              = time()
    rd.poolingRequests = 0
    rd.oldRequests     = 0

    Event.on(rd.poolingCon.observer, Connection.EVENT_BEFORE_RESPONSE, rd.poolingBeforeCb)
    Event.on(rd.cmdCon.observer, Connection.EVENT_AFTER_REQUEST, rd.cmdAfterCb)
    Client.request(rd.cmdCon, RpcApi.RPC_SET_WORLD_STREAMING)
  end
  #
  # Stops displaying organism's world. Closes the connection.
  # @param rd Remote world data object
  #
  function stop(rd::RemoteDataRT)
    Event.off(rd.poolingCon.observer, Connection.EVENT_BEFORE_RESPONSE, rd.poolingBeforeCb)
    Event.off(rd.cmdCon.observer, Connection.EVENT_AFTER_REQUEST, rd.cmdAfterCb)
    Client.stop(rd.cmdCon)
    Client.stop(rd.poolingCon)
  end
  #
  # Handler of remote server pooling request. It may contain two types
  # of data:
  #    x::Uint16, y::UInt16, color::UInt16, ips::UInt16 or
  #    ips::UInt16
  # @param rd Remote Data object
  # @param data Command related data
  #
  function _onDot(rd::RemoteDataRT, data::Array{Any, 1})
    local paramAmount::Int = length(data)
    if time() - rd.ts > 1.0
      rd.ts = time()
      rd.oldRequests = rd.poolingRequests
      rd.poolingRequests = 0
    end
    rd.poolingRequests += 1
    # TODO: it's better to check type of request by FastApi.API_XXX constants
    if paramAmount === 5 _onOrganismDot(rd, data)
    elseif paramAmount === 3 _onEnergyDot(rd, data)
    end
  end
  #
  # Request for drawing organism dot
  # @param rd Remote Data object
  # @param data Command related data
  #
  function _onOrganismDot(rd::RemoteDataRT, data::Array{Any, 1})
    local x::UInt16       = UInt16(data[1])
    local y::UInt16       = UInt16(data[2])
    local nibbles::UInt16 = UInt16(data[3])
    local infoBits::UInt8 = UInt8(data[5])
    local color::UInt16   = nibbles & 0x0fff # last 3 nibbles are color
    local dir::Int        = nibbles >> 12    # first nibble is direction
    local outOfBorder::Bool = infoBits & 0b1 > 0
    local pixelDiff::Dict

    local orgId::UInt = UInt(data[4]) # is an organism id

    local sourceX::UInt16 = x;
    local sourceY::UInt16 = y;
    local updateNumber::UInt32 = rd.fileIndex * AMOUNT_OF_RECORDS + length(rd.diffs);

    if color === Dots.INDEX_EMPTY
      pixelDiff = Dict("n" => updateNumber, "id" => orgId, "sx" => sourceX, "sy" => sourceY, "dx" => x, "dy" => y, "c" => color, "a" => "remove")
      push!(rd.diffs, pixelDiff)
    elseif dir === Dots.DIRECTION_NO
      pixelDiff = Dict("n" => updateNumber, "id" => orgId, "sx" => sourceX, "sy" => sourceY, "dx" => x, "dy" => y, "c" => color, "a" => "add")
      push!(rd.diffs, pixelDiff)
    # TODO: this scenario is not valid at the moment
    #elseif dir === Dots.DIRECTION_NO
    #  pixelDiff = Dict("id" => orgId, "sx" => sourceX, "sy" => sourceY, "dx" => x, "dy" => y, "c" => color, "a" => "color")
    #  push!(rd.diffs, pixelDiff)
    else
      if     dir === Dots.DIRECTION_UP
        if outOfBorder sourceY = rd.height
        else sourceY = y + 0x0001
        end
      elseif dir === Dots.DIRECTION_RIGHT
        if outOfBorder sourceX = 0x0001
        else sourceX = x - 0x0001
        end
      elseif dir === Dots.DIRECTION_DOWN
        if outOfBorder sourceY = 0x0001
        else sourceY = y - 0x0001
        end
      elseif dir === Dots.DIRECTION_LEFT
        if outOfBorder sourceX = rd.width
        else sourceX = x + 0x0001
        end
      end

      pixelDiff = Dict("n" => updateNumber, "id" => orgId, "sx" => sourceX, "sy" => sourceY, "dx" => x, "dy" => y, "a" => "move", "c" => color)
      push!(rd.diffs, pixelDiff)
    end

    if length(rd.diffs) > AMOUNT_OF_RECORDS
      local file::String = string(JSON_FOLDER, "/", lpad(rd.fileIndex, 4, "0"), ".json")
      Helper.save(JSON.json(rd.diffs), file, true)
      rd.diffs = []
      rd.fileIndex += 1
    end
  end

  #
  # Request for drawing energy dot. If we are here, then new energy block was created
  # @param rd Remote Data object
  # @param data Command related data
  #
  function _onEnergyDot(rd::RemoteDataRT, data::Array{Any, 1})

    local x::UInt16       = UInt16(data[1])
    local y::UInt16       = UInt16(data[2])
    local color::UInt16   = UInt16(data[3])
    local updateNumber::UInt32 = rd.fileIndex * AMOUNT_OF_RECORDS + length(rd.diffs);
    local pixelDiff::Dict

    if color === Dots.INDEX_EMPTY
      pixelDiff = Dict("n" => updateNumber, "id" => 0, "sx" => x, "sy" => y, "dx" => x, "dy" => y, "c" => color, "a" => "remove")
    else
      pixelDiff = Dict("n" => updateNumber, "id" => 0, "sx" => x, "sy" => y, "dx" => x, "dy" => y, "c" => color, "a" => "add")
    end
    push!(rd.diffs, pixelDiff)

    if length(rd.diffs) > AMOUNT_OF_RECORDS
      local file::String = string(JSON_FOLDER, "/", lpad(rd.fileIndex, 4, "0"), ".json")
      Helper.save(JSON.json(rd.diffs), file, true)
      rd.diffs = []
      rd.fileIndex += 1
    end

  end

  #
  # Handler of RpcApi.RPC_SET_WORLD_STREAMING request
  # @param data Answer object with region data
  #
  function _onRegion(rd::RemoteDataRT, ans::Connection.Answer)
    local region::RpcApi.Region = ans.data
    local blockIdx::Int
    local block::RpcApi.Block
    local offs::UInt8
    local org::RpcApi.Org
    local xOffs::UInt16
    local yOffs::UInt16

    keyFrameRegion = []

    for blockIdx = 0:length(region.blocks) - 1
      block = region.blocks[blockIdx + 1]
      yOffs = div(blockIdx, region.xBlocks) * RpcApi.BLOCK_SIZE
      xOffs = blockIdx % region.xBlocks * RpcApi.BLOCK_SIZE
      for offs in block.energy

        local x::UInt16       = xOffs + UInt16(offs % RpcApi.BLOCK_SIZE + 1)
        local y::UInt16       = yOffs + UInt16(div(offs, RpcApi.BLOCK_SIZE) + 1)

        pixel = Dict("id" => 0, "x" => x, "y" => y, "c" =>  UInt16(Dots.INDEX_ENERGY))
        push!(keyFrameRegion, pixel)
      end
      for org in block.orgs

        local x::UInt16       = xOffs + UInt16(org.offs % RpcApi.BLOCK_SIZE + 1)
        local y::UInt16       = yOffs + UInt16(div(org.offs, RpcApi.BLOCK_SIZE) + 1)

        pixel = Dict("id" => org.id, "x" => x, "y" => y, "c" =>  org.color)
        push!(keyFrameRegion, pixel)
      end
    end

    keyFrame = Dict("region" => keyFrameRegion,
                    "width" => region.width,
                    "height" => region.height)

    local fileName::String = string(JSON_FOLDER, "/keyframe.json")
    Helper.save(JSON.json(keyFrame), fileName, true)

    #
    # This command means "turn on server pooling"
    #
    Client.request(rd.poolingCon, UInt8(FastApi.API_UINT8), UInt8(0))
  end

end

ext: jl

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
include("../../../src/plotSol.jl")
include("../../../src/simulate.jl")
include("../../../src/linearize.jl")
include("../../../src/OrientationConversion.jl")
# include("../../../src/trim_kronLazy.jl");
include("DoublePendulumMinReal.jl")
using Revise
using BenchmarkTools
## clearing variables
# clearconsole()
j1 = Nothing
j2 = Nothing
##
m = 1.0; l = 1.0; # Mass and length of bar
# Assuming bar revolves about Y axis
I1 = [1 0 0; 0 m*l^2/12 0; 0 0 1]

# Testing Dynamics with Revolute Joint
R1 = RigidBody(m,I1,2)
R2 = RigidBody(m,I1,3)
RbI = InertialFrameAsRB()

# Perturbed from equilibrium
theta1 = deg2rad(30); theta2 = deg2rad(45);
x1_temp = [l/2*sin(theta1);0.0;-l/2*cos(theta1)];
x2_temp = 2*x1_temp + [l/2*sin(theta2);0.0;-l/2*cos(theta2)];
x0R1 = [x1_temp;[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(R1,x0R1)
x0R2 = [x2_temp;[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(R2,x0R2)

axisY = [0.0 1.0 0.0][:] # Axis about which bar is revolving

rj1 = [0.0 0.0 0.0][:] # Joint Location in body frame of first body
rj2 = -R1.x[1:3]

rj3 = x1_temp; # Joint location in body frame of pend1
rj4 = -x2_temp + 2*x1_temp; # Joint location in body frame of pend2

j1 = Joint(RbI,R1,rj1,rj2,type="Revolute",axis=axisY)
j2 = Joint(R1,R2,rj3,rj4,type="Revolute",axis=axisY)

# External Forces Definition
g = [0.0,0.0,-9.806]; # Gravity Vector.

j = [j1;j2];

## Linearization
# x0, u0 = getXU_0(j);
# A,B = linearize(x0, u0, j, g);

## Simulation
tEnd = 1.0
tInt = 0.01;

tSim, solFinal = simulate(tEnd,tSpan,j,g=g)
println("Simulation complete; plotting results:")

solPend1  = solFinal[1];
solPend2  = solFinal[2];
##
ωPend1, jointLoc1 = checkRevJointIn(solPend1, rj2)
ωPend1, ωPend2, ωRel2, jointLoc2 = checkRevJoint(solPend1, solPend2, rj3, rj4);
##
close("all");
# errors
plotErrNorm(tSim,solPend1.β)
plotErrNorm(tSim,solPend2.β)
plotPos(tSim,jointLoc1);
plotPos(tSim,jointLoc2);

# trajectory
plotPos(tSim, solPend1.r)
plotPos(tSim, solPend2.r)
## Compare with MinReal
r1Min, r2Min, v1Min, v2Min, ω1Min, ω2Min = DoublePendulumMinReal(m,l,theta1,theta2, tInt, tEnd);

# Plotting errors vs MinReal
close("all");
plotPos(tSim,r1Min-solPend1.r)
plotPos(tSim,r2Min-solPend2.r)
plotPos(tSim,v1Min-solPend1.v)
plotPos(tSim,v2Min-solPend2.v)
plotAngVel(tSim,ω1Min-ωPend1)
plotAngVel(tSim,ω2Min-ωPend2)

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
module GridCreation

export CFFT, CInterSample, mkCartVecs, cartVecsToPolarGrid

const CFFT = -1;
const CInterSample = -2;

# this function is lifted from pkg Fourier.
# A little copying is better than a little dependency.
function _fftrange(n)
    if iseven(n)
        return -(n ÷ 2):(n÷2-1)
    end
    return -(n ÷ 2):(n÷2)
end
"""
    mkCartVecs(dx, size[, center])

Create a 2D "grid" by returning a row vector x and column vector y.

The grid will have sample spacing dx, size[0] == size(x), and size[1] == size(y).

If size is an int, it is broadcast to [size,size].

The grid will be centered according to such that the center element contains zero.  Special
values for center include -1 (default, FFT-like center) and -2, for "interpixel"
sampling, where 0 is at size/2.  Otherwise, zero is at the centerth element.

The centering rule may not be different between x and y.

# Examples

## FFT-aligned grid

Note that the ceil rounded center element contains zero.
```julia-repl
julia> x,y = mkCartVecs(.1, 4); # implicit center=CFFT
julia> x
1×4 LinearAlgebra.Adjoint{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}:
 -0.2  -0.1  0.0  0.1
julia> collect(y)
4-element Array{Float64,1}:
 -0.2
 -0.1
  0.0
  0.1
```

## Intersample aligned grid
```julia-repl
x,y=mkCartVecs(.1, 4, center=CInterSample);

julia> x
1×4 LinearAlgebra.Adjoint{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}:
 -0.15  -0.05  0.05  0.15
```

## index-centered grid
```
julia> x,y=mkCartVecs(.1, 4, center=1);
julia> x
1×4 LinearAlgebra.Adjoint{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}:
 1.11022e-17  0.1  0.2  0.3
```
"""
function mkCartVecs(dx, size::Tuple{Integer,Integer}; center::Integer=CFFT)
    # fftrange produces a range object.  dot add shifts the range
    # and multiplication scales it.
    if center == CFFT
        X = (_fftrange(size[2])*dx)';
        Y = _fftrange(size[1])*dx;
    elseif center == CInterSample
        X = (_fftrange(size[2])*dx.+(dx/2))';
        Y = (_fftrange(size[1])*dx.+(dx/2));
    else
        X = _fftrange(size[2])*dx;
        Y = _fftrange(size[1])*dx;
        X = (X .- X[center])';
        Y = (Y .- Y[center]);
    end
    return X, Y;
end

function mkCartVecs(dx, size::Integer; center::Integer=CFFT)
    return mkCartVecs(dx, (size,size), center=center);
end

"""
    cartVecsToPolarGrid(X, Y)

Construct a polar grid (ρ,θ) from the row vector X and column vector Y.  Returns
a pair of 2D arrays containing the radial and azimuthal coordinates.

# Examples
```julia-repl
julia> X=(-2:1)'; Y=-2:1;
julia> ρ,θ=cartVecsToPolarGrid(X,Y);
julia> ρ
4×4 Array{Float64,2}:
 2.82843  2.23607  2.0  2.23607
 2.23607  1.41421  1.0  1.41421
 2.0      1.0      0.0  1.0
 2.23607  1.41421  1.0  1.41421
julia> θ
4×4 Array{Float64,2}:
 -2.35619  -2.03444  -1.5708  -1.10715
 -2.67795  -2.35619  -1.5708  -0.785398
  3.14159   3.14159   0.0      0.0
  2.67795   2.35619   1.5708   0.785398
```
See also: [`mkCartVecs`](@ref)
"""
function cartVecsToPolarGrid(X, Y)
    # I don't think allocating a meshgrid is avoidable here.
    # the rho matrix can be made with a dot on x and y,
    # but we would need a 2D view into both X and Y for atan
    # grid_x = [x for x in X, y in Y;
    # grid_y = [y for y in Y, x in X;
    ρ = @. sqrt(X^2 + Y^2);
    θ = [atan(y,x) for y in Y, x in X'];
    return ρ, θ
end

end

max_stars_repo_name: JuliaOptics/Gridding

max_issues_repo_name: JuliaOptics/Gridding

max_forks_repo_name: JuliaOptics/Gridding


text:
using WAVEFORM, SparseArrays

n = 101 .* [1;1]
d = 10.0 .*[1;1]
npml = [1 1; 1 1]
freq = 5.0
v = ones((n...)) .* 2000
f0 = 10.0
unit = "m/s"

H, dH, ddH = helm2d_7pt(n,d,npml,freq,v,f0,unit)    # 10201
H, dH, ddH = WAVEFORM.helm2d_chen2013(n,d,npml,freq,v,f0,unit)


text:
module NASAEarthObservations

using Glob

using Pkg.Artifacts

const rooturl = "neo.sci.gsfc.nasa.gov"

function _download_files(root, intermediate, target; dir = mktempdir())

    imgpath = "$root/$intermediate/$target"

    cd(dir) do
        run(```
            wget
                --no-parent
                --no-directories
                --no-verbose
                --reject="index.html*"
                --execute robots=off
                --include $intermediate/$target
                -r
                https://$imgpath
            ```
            )
    end
    return dir
end

download_files(dirname; root = rooturl) = _download_files(root, "archive", dirname)

"""
    observations(path)

Downloads (if necessary) and stores the data found at the path given, relative
to [`$rooturl`]($rooturl).

This method downloads the observations and then stores them as an artifact,
so that redownloading becomes unneccessary.
"""
function observations(path)

    if startswith(path, "/")
        @error "Path cannot start with `/`!" path
    elseif startswith(path, "archives")
        @error "Please give the path after archive!" path
    end

    mutable_artifacts_toml = joinpath(dirname(@__DIR__), "MutableArtifacts.toml")

    dir_name = endswith("/", path) ? path[1:end-1] : path

    dir_hash = artifact_hash(dir_name, mutable_artifacts_toml)

    if dir_hash === nothing

        dir_hash = create_artifact() do dir
            cd(dir) do
                _download_files(rooturl, "archive", path; dir = dir)
            end
        end

        bind_artifact!(
            mutable_artifacts_toml,
            dir_name,
            dir_hash;
            force=true
            )

    end

    return artifact_path(dir_hash)

end

export observations

function __init__()
    touch(joinpath(dirname(@__DIR__), "MutableArtifacts.toml"))
end

end # module


text:
function V=RungeKuttaHEVI(V,fSlow,dt,FcnSlow,FcnFast,L,U,p,q,Param)
Vn=V;
RK=Param.RK;
NN=Param.NN;
nV=Param.nV;
fVE=zeros([size(V) RK.nStage]);
fVI=zeros([size(V) RK.nStage]);
for iStage=1:RK.nStage
  V=Vn;
  if iStage>1
    for jStage=1:iStage-1
      V=V+dt*RK.ARKE(iStage,jStage)*fVE(:,:,:,:,jStage)+...
          dt*RK.ARKI(iStage,jStage)*fVI(:,:,:,:,jStage);
    end
    v=reshape(V,NN*nV,1);
    v(q)=U\(L\v(p));
    V=reshape(v,size(Vn));
  end
  fVE(:,:,:,:,iStage)=FcnSlow(V,Param)+fSlow;
  fVI(:,:,:,:,iStage)=FcnFast(V,Param);
end
V=Vn;
for iStage=1:RK.nStage
  V=V+dt*RK.bRKE(iStage)*fVE(:,:,:,:,iStage)+...
      dt*RK.bRKI(iStage)*fVI(:,:,:,:,iStage);
end
end


text:
using QuantumCumulants
using ModelingToolkit, OrdinaryDiffEq
using Test

@testset "two-level-laser" begin

N = 10

Δ = cnumbers((Symbol(:Δ_, i) for i=1:N)...)
g = cnumbers((Symbol(:g_, i) for i=1:N)...)
γ = cnumbers((Symbol(:γ_, i) for i=1:N)...)
ν = cnumbers((Symbol(:ν_, i) for i=1:N)...)
@cnumbers κ

h_cavity = FockSpace(:cavity)
h_atoms = [NLevelSpace(Symbol(:atom, i),(:g,:e)) for i=1:N]
h = tensor(h_cavity, h_atoms...)

a = Destroy(h, :a)
σ(i,j,k) = Transition(h,Symbol(:σ_,k),i,j,k+1)

ops = begin
    ops_ = [a'*a;[σ(:e,:e,k) for k=1:N];[a'*σ(:g,:e,k) for k=1:N]]
    for i=1:N
        for j=i+1:N
            push!(ops_, σ(:e,:g,i)*σ(:g,:e,j))
        end
    end
    ops_
end

n_eqs = div(N*(N-1),2) + 2N + 1

H = sum(Δ[i]*σ(:e,:e,i) for i=1:N) + sum(g[i]*(a'*σ(:g,:e,i) + a*σ(:e,:g,i)) for i=1:N)

J = [a;[σ(:g,:e,k) for k=1:N];[σ(:e,:g,k) for k=1:N]]
Jdagger = adjoint.(J)
rates = [κ,γ...,ν...]
he = meanfield(ops, H, J; Jdagger=Jdagger, rates=rates, simplify=true, order=2)

missed = find_missing(he)

subs = Dict(missed .=> 0)

he_nophase = substitute(he, subs)

@test isempty(find_missing(he_nophase))

eqs_mtk = equations(he_nophase)

sys = ODESystem(he_nophase)

u0 = zeros(ComplexF64, length(ops))
p0 = [κ => 1,
    (γ .=> 0.25 .* ones(N))...,
    (ν .=> 4 .* ones(N))...,
    (g .=> 1.5 .* ones(N))...,
    (Δ .=> ones(N))...,]

prob = ODEProblem(sys,u0,(0.0,10.0),p0,jac=true,sparse=true)

sol = solve(prob,RK4())

n = real.(sol[average(a'*a)])
pe = sol[average(σ(:e,:e,1))]

# Test with complete and custom filter
ϕ(::Destroy) = -1
ϕ(::Create) = 1
function ϕ(t::Transition)
    if t.i != t.j
        t.i == :e && return 1
        return -1
    else
        return 0
    end
end
function ϕ(q::QuantumCumulants.QMul)
    p = 0
    for arg ∈ q.args_nc
        p += ϕ(arg)
    end
    return p
end
ϕ(avg::Average) = ϕ(avg.arguments[1])
phase_invariant(x) = iszero(ϕ(x))

he_n = meanfield(a'*a, H, J; rates=rates)
complete!(he_n;filter_func=phase_invariant)

@test length(he_n.equations) == length(ops)
@test isempty(find_missing(he_n))

sys_comp = ODESystem(he_n)
prob_comp = ODEProblem(sys_comp,u0,(0.0,10.0),p0)

sol_comp = solve(prob_comp,RK4())

@test getindex.(sol.u, 1) ≈ getindex.(sol_comp.u, 1)

end # testset


text:
import Base:showerror
export ip_error
abstract type InplaceException <: Exception end

struct InplaceException1 <: InplaceException
    msg::String
end
showerror(io::IO, e::InplaceException1) = print(io, "InplaceException: ", e.msg)

struct InplaceException2 <: InplaceException
    blas_error::Exception
end

function showerror(io::IO, e::InplaceException2)
    println(io, "InplaceException: inplace assignment for these arguments has failed.")
    println(io, "Try some other combination, or try plain Julia.")
    println(io, "The original exception was:")
    showerror(io::IO, e.blas_error)
    stacktrace(catch_backtrace())
end

ip_error(msg::String) = throw(InplaceException1(msg))
ip_error(err::Exception) = throw(InplaceException2(err))


text:
using UnitJuMP
using GLPK

m = Model(GLPK.Optimizer)

const speedunit = u"m/s"
const timeunit = u"s"

@variable(m, x >= 0, speedunit)
@variable(m, y >= 0, timeunit)
@variable(m, z >= 0, u"km/s")

@variable(m, v[1:4], speedunit)

@variable(m, w, Bin)

max_speed = 2u"m/s"
a = 9.8u"ft/s^2"

con = @constraint(m, x + a * y <= max_speed)
@constraint(m, z + x <= max_speed)
@constraint(m, x + 0.5 * a * y  <= max_speed)

M = 0.5u"m/s"
@constraint(m, x <= M * w)

no_speed = 0.0u"m/s"
@constraint(m, -3x  + 2z <= no_speed)

@constraint(m, x == 0.15 * max_speed)

@constraint(m, con7[i=1:4], v[i] <= max_speed)

@constraint(m, con8, sum(v[i] for i=1:3) <= max_speed, u"ft/s")

@objective(m, Max, x + 0.2*z)

data = 4
@constraint(m, -3*data*x  + 2z <= no_speed)

optimize!(m)

value(m[:x])
value(m[:y])
value(m[:z])


text:
using Turing, Distributions, DifferentialEquations

# Import MCMCChain, Plots, and StatsPlots for visualizations and diagnostics.
using MCMCChains, Plots, StatsPlots

# Set a seed for reproducibility.
using Random
Random.seed!(14);

function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, γ, δ  = p
  du[1] = (α - β*y)x # dx =
  du[2] = (δ*x - γ)y # dy =
end
p = [1.5, 1.0, 3.0, 1.0]
u0 = [1.0,1.0]
prob1 = ODEProblem(lotka_volterra,u0,(0.0,10.0),p)
sol = solve(prob1,Tsit5())
plot(sol)

sol1 = solve(prob1,Tsit5(),saveat=0.1)
odedata = Array(sol1) + 0.8 * randn(size(Array(sol1)))
plot(sol1, alpha = 0.3, legend = false); scatter!(sol1.t, odedata')

Turing.setadbackend(:forwarddiff)

@model function fitlv(data, prob1)
    σ ~ InverseGamma(2, 3) # ~ is the tilde character
    α ~ truncated(Normal(1.5,0.5),0.5,2.5)
    β ~ truncated(Normal(1.2,0.5),0,2)
    γ ~ truncated(Normal(3.0,0.5),1,4)
    δ ~ truncated(Normal(1.0,0.5),0,2)

    p = [α,β,γ,δ]
    prob = remake(prob1, p=p)
    predicted = solve(prob,Tsit5(),saveat=0.1)

    for i = 1:length(predicted)
        data[:,i] ~ MvNormal(predicted[i], σ)
    end
end

model = fitlv(odedata, prob1)

# This next command runs 3 independent chains without using multithreading.
chain = mapreduce(c -> sample(model, NUTS(.65),1000), chainscat, 1:3)

plot(chain)

pl = scatter(sol1.t, odedata');

chain_array = Array(chain)
for k in 1:300
    resol = solve(remake(prob1,p=chain_array[rand(1:1500), 1:4]),Tsit5(),saveat=0.1)
    plot!(resol, alpha=0.1, color = "#BBBBBB", legend = false)
end
# display(pl)
plot!(sol1, w=1, legend = false)

@model function fitlv2(data, prob1) # data should be a Vector
    σ ~ InverseGamma(2, 3) # ~ is the tilde character
    α ~ truncated(Normal(1.5,0.5),0.5,2.5)
    β ~ truncated(Normal(1.2,0.5),0,2)
    γ ~ truncated(Normal(3.0,0.5),1,4)
    δ ~ truncated(Normal(1.0,0.5),0,2)

    p = [α,β,γ,δ]
    prob = remake(prob1, p=p)
    predicted = solve(prob,Tsit5(),saveat=0.1)

    for i = 1:length(predicted)
        data[i] ~ Normal(predicted[i][2], σ) # predicted[i][2] is the data for y - a scalar, so we use Normal instead of MvNormal
    end
end

model2 = fitlv2(odedata[2,:], prob1)

Threads.nthreads()

# This next command runs 3 independent chains with multithreading.
chain2 = sample(model2, NUTS(.45), MCMCThreads(), 5000, 3, progress=false)

pl = scatter(sol1.t, odedata');
chain_array2 = Array(chain2)
for k in 1:300
    resol = solve(remake(prob1,p=chain_array2[rand(1:12000), 1:4]),Tsit5(),saveat=0.1)
    # Note that due to a bug in AxisArray, the variables from the chain will be returned always in
    # the order it is stored in the array, not by the specified order in the call - :α, :β, :γ, :δ
    plot!(resol, alpha=0.1, color = "#BBBBBB", legend = false)
end
#display(pl)
plot!(sol1, w=1, legend = false)

function delay_lotka_volterra(du, u, h, p, t)
   x, y = u
   α, β, γ, δ = p
   du[1] = α * h(p, t-1; idxs=1) - β * x * y
   du[2] = -γ * y + δ * x * y
   return
end

p = (1.5,1.0,3.0,1.0)
u0 = [1.0; 1.0]
tspan = (0.0,10.0)
h(p, t; idxs::Int) = 1.0
prob1 = DDEProblem(delay_lotka_volterra,u0,h,tspan,p)

sol = solve(prob1,saveat=0.1)
ddedata = Array(sol)
ddedata = ddedata + 0.5 * randn(size(ddedata))

scatter(sol.t, ddedata'); plot!(sol)

Turing.setadbackend(:forwarddiff)
@model function fitlv(data, prob1)

    σ ~ InverseGamma(2, 3)
    α ~ Truncated(Normal(1.5,0.5),0.5,2.5)
    β ~ Truncated(Normal(1.2,0.5),0,2)
    γ ~ Truncated(Normal(3.0,0.5),1,4)
    δ ~ Truncated(Normal(1.0,0.5),0,2)

    p = [α,β,γ,δ]

    #prob = DDEProblem(delay_lotka_volterra,u0,_h,tspan,p)
    prob = remake(prob1, p=p)
    predicted = solve(prob,saveat=0.1)
    for i = 1:length(predicted)
        data[:,i] ~ MvNormal(predicted[i], σ)
    end
end;
model = fitlv(ddedata, prob1)

chain = sample(model, NUTS(.65), MCMCThreads(), 300, 3, progress=true)
plot(chain)

chain

pl = scatter(sol.t, ddedata')
chain_array = Array(chain)
for k in 1:100
    resol = solve(remake(prob1,p=chain_array[rand(1:450),1:4]),Tsit5(),saveat=0.1)
    # Note that due to a bug in AxisArray, the variables from the chain will be returned always in
    # the order it is stored in the array, not by the specified order in the call - :α, :β, :γ, :δ

    plot!(resol, alpha=0.1, color = "#BBBBBB", legend = false)
end
#display(pl)
plot!(sol)

using Zygote, DiffEqSensitivity
Turing.setadbackend(:zygote)
prob1 = ODEProblem(lotka_volterra,u0,(0.0,10.0),p)

@model function fitlv(data, prob)
    σ ~ InverseGamma(2, 3)
    α ~ truncated(Normal(1.5,0.5),0.5,2.5)
    β ~ truncated(Normal(1.2,0.5),0,2)
    γ ~ truncated(Normal(3.0,0.5),1,4)
    δ ~ truncated(Normal(1.0,0.5),0,2)
    p = [α,β,γ,δ]
    prob = remake(prob, p=p)

    predicted = solve(prob,saveat=0.1)
    for i = 1:length(predicted)
        data[:,i] ~ MvNormal(predicted[i], σ)
    end
end;
model = fitlv(odedata, prob1)
chain = sample(model, NUTS(.65),1000)

@model function fitlv(data, prob)
    σ ~ InverseGamma(2, 3)
    α ~ truncated(Normal(1.5,0.5),0.5,2.5)
    β ~ truncated(Normal(1.2,0.5),0,2)
    γ ~ truncated(Normal(3.0,0.5),1,4)
    δ ~ truncated(Normal(1.0,0.5),0,2)
    p = [α,β,γ,δ]
    prob = remake(prob, p=p)
    predicted = solve(prob,saveat=0.1,sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))
    for i = 1:length(predicted)
        data[:,i] ~ MvNormal(predicted[i], σ)
    end
end;
model = fitlv(odedata, prob1)
@time chain = sample(model, NUTS(.65),1000)

u0 = [1.0,1.0]
tspan = (0.0,10.0)
function multiplicative_noise!(du,u,p,t)
  x,y = u
  du[1] = p[5]*x
  du[2] = p[6]*y
end
p = [1.5,1.0,3.0,1.0,0.1,0.1]

function lotka_volterra!(du,u,p,t)
  x,y = u
  α,β,γ,δ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = δ*x*y - γ*y
end

prob_sde = SDEProblem(lotka_volterra!,multiplicative_noise!,u0,tspan,p)

ensembleprob = EnsembleProblem(prob_sde)
@time data = solve(ensembleprob,SOSRI(),saveat=0.1,trajectories=1000)
plot(EnsembleSummary(data))

Turing.setadbackend(:forwarddiff)
@model function fitlv(data, prob)
    σ ~ InverseGamma(2,3)
    α ~ truncated(Normal(1.3,0.5),0.5,2.5)
    β ~ truncated(Normal(1.2,0.25),0.5,2)
    γ ~ truncated(Normal(3.2,0.25),2.2,4.0)
    δ ~ truncated(Normal(1.2,0.25),0.5,2.0)
    ϕ1 ~ truncated(Normal(0.12,0.3),0.05,0.25)
    ϕ2 ~ truncated(Normal(0.12,0.3),0.05,0.25)
    p = [α,β,γ,δ,ϕ1,ϕ2]
    prob = remake(prob, p=p)
    predicted = solve(prob,SOSRI(),saveat=0.1)

    if predicted.retcode != :Success
        Turing.acclogp!(_varinfo, -Inf)
    end
    for j in 1:length(data)
        for i = 1:length(predicted)
            data[j][i] ~ MvNormal(predicted[i],σ)
        end
    end
end;

model = fitlv(data, prob_sde)
chain = sample(model, NUTS(0.25), 5000, init_theta = [1.5,1.3,1.2,2.7,1.2,0.12,0.12])
plot(chain)

if isdefined(Main, :TuringTutorials)
    Main.TuringTutorials.tutorial_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
end

max_stars_repo_name: leachim/TuringTutorials

max_issues_repo_name: leachim/TuringTutorials

max_forks_repo_name: leachim/TuringTutorials


text:
using Test
using Smooth
using Plots

# open("log.log","w") do io
#     redirect_stdout(io)
# end

@testset "CanvasModule.jl" begin
X_DIM = 1800
Y_DIM = 1800

canvas = Canvas(X_DIM, Y_DIM)
initialize(canvas,20,2000,10)

Points = collect(Iterators.flatten([canvas.clusters[i].points for i in 1:length(canvas.clusters)]))
X = [ Points[i][1] for i in 1:length(Points) ]
Y = [ Points[i][2] for i in 1:length(Points) ]

@test length(filter(x -> x >X_DIM,X)) == 0
@test length(filter(x -> x >Y_DIM,Y)) == 0

# gr()
plot( X, Y, seriestype = :scatter,xlims = (0,X_DIM), ylims = (0,Y_DIM))
savefig("canvas.png")

#Make a move
m = MobileUser(0,(100,100),canvas,1,0.3)
move(m,30)

end

@testset "LATP.jl" begin
    g = SimpleWeightedGraph(0)
    for i in 1:5
        add_vertex!(g)
        for j in 1:nv(g) - 1
            add_edge!(g,i,j,rand(1:15))
        end
    end
    latp_algorithm(g,1,3)
end


text:
# This script loads the raw data from all experiments, runs the protocols and
# stores some values that will be useful to run the tests and the extraction of
# the results from a simulated run in fields of type Base.RefValue.
# The results are saved and are loaded by FoodCachingExperiments at init.
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
using TOML, CSV, FoodCachingExperiments, DataFramesMeta, Distributions, BSON, CodecZstd
import FoodCachingExperiments: EXPERIMENTS, Experiment, TestSummary, Food,
       default_keycols, bsave
for foodtype in instances(Food)
    eval(:(import FoodCachingExperiments: $(Symbol(foodtype))))
end

include(joinpath(@__DIR__, "..", "examples", "dummy_model.jl"))

const DATADIR = joinpath(@__DIR__, "..", "data", "raw")

function load_experiment(name; datadir = DATADIR, run = true)
    meta = TOML.parsefile(joinpath(datadir, "$name.toml"))
    data = CSV.File(joinpath(datadir, "$name.csv"), pool = false) |> DataFrame
    tests = CSV.File(joinpath(datadir, "$(name)_tests.csv")) |> DataFrame |>
            processtests
    ex = Experiment{name}(Base.RefValue{Any}(),
                          Base.RefValue{Any}(),
                          haskey(meta, "major_finding") ?
                          meta["major_finding"] : (@warn("$name misses major_finding"); ""),
                          haskey(meta, "comments") ?
                          meta["comments"] : (@warn("$name misses comments"); ""),
                          sort!(data, default_keycols(data)),
                          tests, meta["nbirds"],
                          eval.(Meta.parse.(meta["foodtypes"])))
    if run
        data = run!(ex, fill(DummyModel(), ex.nbirds))
        res = results(ex, data)
        ex = dropref(ex)
        length(res) == length(target(ex)) || @warn "$name: result and target have different lengths"
        ex
    else
        ex
    end
end

function processtests(tests)
    result = Dict{String,TestSummary}()
    for g in groupby(tests, :id)
        result[g.id[1]] = TestSummary(
                                @with(g, pvalue.(:statistic, :df, :value)),
                                string.(g.key),
                                string.(g.df),
                                string.(g.label),
                                string.(g.statistic),
                                string.(g.value),
                                Base.RefValue{Any}(),
                                Base.RefValue{Any}()
                               )
    end
    result
end
pvalue(statistic, df::String, value) = pvalue(statistic, eval(Meta.parse(df)), value)
function pvalue(statistic, df, value)
    if typeof(value) == String
        if value[1] == '<'
            value = eval(Meta.parse(value[2:end])) - .01
        elseif value[1] == '>'
            value = eval(Meta.parse(value[2:end])) + .01
        else
            value = eval(Meta.parse(value))
        end
    end
    if statistic == "F"
        ccdf(FDist(df[1], df[2]), value)
    elseif statistic == "t"
        ccdf(FDist(1, df[1]), value^2)
    elseif statistic == "U" || statistic == "FisherExactTest" || statistic == "Binomial"
        value
    elseif statistic == "Chisq"
        ccdf(Chisq(df), value)
    else
        error("Don't know how to treat statistic $statistic.")
    end
end

function dropref(x::TestSummary)
    TestSummary(x.pvalues, x.keys, x.df, x.labels, x.statistic, x.value,
                isassigned(x.locals) ? x.locals[] : nothing,
                isassigned(x.accessor) ? x.accessor[] : nothing)
end
function dropref(x::Dict{String, TestSummary})
    Dict(id => dropref(test) for (id, test) in x)
end
function dropref(e::Experiment{name}) where name
    Experiment{name}(e.data_accessor[],
                     e.target[],
                     e.major_finding,
                     e.comments,
                     e.data,
                     dropref(e.tests),
                     e.nbirds,
                     e.foodtypes)
end

println("Processing experiments. This will take some time!")
for filename in readdir(DATADIR)
    name, ending = splitext(filename)
    name = Symbol(name)
    run = name ∉ FoodCachingExperiments.CLAYTON0103_EXPERIMENTS
    if ending == ".toml"
        println("$name ...")
        EXPERIMENTS[name] = load_experiment(name; run)
    end
end

# This is a dummy experiment to aggregate all the sequentially performed experiments
println("Clayton0103 ...")
EXPERIMENTS[:Clayton0103] = Experiment{:Clayton0103}(Base.RefValue{Any}(),
                                                     Base.RefValue{Any}(),
                                                     "",
                                                     "Running 01/03 experiments in one go.",
                                                     DataFrame(),
                                                     Dict{String, TestSummary}(),
                                                     16,
                                                     [Peanut, Mealworm, Cricket])
data = run!(EXPERIMENTS[:Clayton0103], fill(DummyModel(), 16))
res = results(EXPERIMENTS[:Clayton0103], data)
EXPERIMENTS[:Clayton0103] = dropref(EXPERIMENTS[:Clayton0103])
length(res) == length(target(EXPERIMENTS[:Clayton0103])) || @warn "Clayton0103: result and target have different lengths"
for e in FoodCachingExperiments.CLAYTON0103_EXPERIMENTS
    EXPERIMENTS[e] = dropref(EXPERIMENTS[e])
end

bsave(joinpath(DATADIR, "..", "processed", "experiments"), EXPERIMENTS)


text:
"
    createLog(params)

Initializes the log dictionary for iteration progress data storage.
"
function create_log(params)
    logVar = Dict()

    # iterate through the params type and create apropriate fields
    for name in fieldnames(params)
        strname = string(name) # otherwise we would have symbols and would have
        # to call dict[:name] instead of dict["name"]

        # if some of the fields in params are not array, create an array from
        # them so that future new values can be concatenated
        if supertype(typeof(getfield(params, name))) == AbstractFloat ||
            supertype(typeof(getfield(params, name))) == Signed
            logVar[strname] = [getfield(params, name)]
        else
            logVar[strname] = copy(getfield(params, name))
        end
    end

    return logVar
end

"
    updateLog!(logVar::Dict{Any,Any}, params)

Concatenates all arrays and floats from params to an existing log dictionary.
"
function update_log!(logVar::Dict{Any,Any}, params)
    # iterate through the params type and create apropriate fields
    for name in fieldnames(params)
        strname = string(name)

        # if some of the fields in params are not array, arrayize it
        if supertype(typeof(getfield(params, name))) == AbstractFloat ||
            supertype(typeof(getfield(params, name))) == Signed
            logVar[strname] = cat(1, logVar[strname], [getfield(params, name)])
        else
            logVar[strname] = cat(ndims(getfield(params, name))+1, logVar[strname], getfield(params, name))
        end
    end
end

"
    saveLog(logVar::Dict{Any, Any}, Y, params, priors, logdir::String; desc::String = "")

Saves the log dictionary and inputs to a .jld file specified in logdir.
Name of the log is date_time if unspecified, or desc.
"
function save_log(logVar::Dict{Any, Any}, Y, priors, logdir::String; desc::String = "")
    # if no description given, create it from current datetime
    if desc == ""
        desc = Dates.format(now(), "yyyymmdd_HHMMSS")
    end

    # create the dir
    dir = string(logdir, "/", desc)
    mkpath(dir)

    # save the data
    save(string(dir, "/log.jld"), logVar)
    save(string(dir, "/inputs.jld"), "Y", Y, "priors", priors)
end

"
    loadLog(path::String)

Saves the log dictionary to a .jld file specified in logdir.
Name of the log is date_time if unspecified, or desc.
"
function load_log(path::String)
    try
        logVar = load(string(path, "/log.jld"))
        inputs = load(string(path, "/inputs.jld"))
        return logVar, inputs["Y"], inputs["priors"]
    catch
        ls = readdir(path)
        msg = "The specified folder does not contain any log files but it contains the following:\n"
        for item in ls
            msg = string(msg, item, "\n")
        end
        error(msg)

        return
    end
end

"
    extractParams!(logVar::Dict{Any, Any}, t::Integer, params)\n

Extracts a single slice with index t from the log dictionary variable given an initialized params type variable.
"
function extract_params!(logVar::Dict, t::Integer, params)
    for name in fieldnames(params)
        strname = string(name)

        # we need just one time slice - time is always the last array dimension
        # but since there are arrays with different number of dimensions
        # we must slice the array dynamically - we cant use [:,t] or [:,:,t]
        dims = size(logVar[strname])
        sliceSize = 1 # how large is one slice
        for d in dims[1:end-1]
            sliceSize*=d
        end

        # this uses linear indices to get the data, then reshapes the final array
        if length(dims) == 1
            setfield!(params, name, logVar[strname][t])
        else
            setfield!(params, name, reshape(logVar[strname][(t-1)*sliceSize+1:t*sliceSize], dims[1:end-1]))
        end
    end

    return params
end


text:
using Plots, Test

@testset "Subplot sclicing" begin
    pl = @test_nowarn plot(
        rand(4, 8),
        layout = 4,
        yscale = [:identity :identity :log10 :log10],
    )
    @test pl[1][:yaxis][:scale] == :identity
    @test pl[2][:yaxis][:scale] == :identity
    @test pl[3][:yaxis][:scale] == :log10
    @test pl[4][:yaxis][:scale] == :log10
end


text:
function spharm( l, m, Gvec )

  Ng = size(Gvec)[2]

  r = zeros( Float64, Ng )
  for ig = 1:Ng
    r[ig] = sqrt( Gvec[1,ig]^2 + Gvec[2,ig]^2 + Gvec[3,ig]^2 )
  end
  # This will fail for very small r[ig]

  ylm = zeros( Float64, Ng )

  if l == 0
    for ig = 1:Ng
      ylm[ig] = 0.5*sqrt(1.0/pi)
    end
    return ylm

  elseif l == 1
    #
    if m == -1
      for ig = 1:Ng
        ylm[ig] = 0.5*sqrt(3/pi)*Gvec[2,ig]/r[ig]
      end
      return ylm
    #
    elseif m == 0
      for ig = 1:Ng
        ylm[ig] = 0.5*sqrt(3/pi)*Gvec[3,ig]/r[ig]
      end
      return ylm
    #
    elseif m == 1
      for ig = 1:Ng
        ylm[ig] = 0.5*sqrt(3/pi)*Gvec[1,ig]/r[ig]
      end
      return ylm
    end

  elseif l == 2
    #
    if m == -2
      for ig = 1:Ng
        ylm[ig] = 0.5*sqrt(15/pi)*Gvec[1,ig]*Gvec[2,ig]/r[ig]^2
      end
      return ylm
    #
    elseif m == -1
      for ig = 1:Ng
        ylm[ig] = 0.5*sqrt(15/pi)*Gvec[2,ig]*Gvec[3,ig]/r[ig]^2
      end
      return ylm
    #
    elseif m == 0
      for ig = 1:Ng
        ylm[ig] = 0.25*sqrt(5/pi)*( -Gvec[1,ig]^2 - Gvec[2,ig]^2 + 2*Gvec[3,ig]^2 )/r[ig]^2
      end
      return ylm
    #
    elseif m == 1
      for ig = 1:Ng
        ylm[ig] = 0.5*sqrt(15/pi)*Gvec[3,ig]*Gvec[1,ig]/r[ig]^2
      end
      return ylm
    #
    elseif m == 2
      for ig = 1:Ng
        ylm[ig] = 0.25*sqrt(15/pi)*(Gvec[1,ig]^2 - Gvec[2,ig]^2)/r[ig]^2
      end
      return ylm

    end

  end

end


text:
import Base: copyto!, copy, deepcopy_internal
using CUDA: CuArray
# include("karray.jl") ## KnetArray

# Avoid Array->Array to prevent base conflict.
copyto!(dest::KnetArray{T}, doffs::Integer, src::KnetArray{T}, soffs::Integer, n::Integer) where {T} = (copyto!(CuArray(dest), doffs, CuArray(src), soffs, n); dest)
copyto!(dest::KnetArray{T}, doffs::Integer, src::Array{T}, soffs::Integer, n::Integer) where {T} = (copyto!(CuArray(dest), doffs, src, soffs, n); dest)
copyto!(dest::Array{T}, doffs::Integer, src::KnetArray{T}, soffs::Integer, n::Integer) where {T} = (copyto!(dest, doffs, CuArray(src), soffs, n); dest)
copyto!(dest::KnetArray{T}, src::KnetArray{T}) where {T} = copyto!(dest, 1, src, 1, length(dest))
copyto!(dest::KnetArray{T}, src::Array{T}) where {T} = copyto!(dest, 1, src, 1, length(dest))
copyto!(dest::Array{T}, src::KnetArray{T}) where {T} = copyto!(dest, 1, src, 1, length(dest))

# Changing types
copyto!(dest::KnetArray{T}, doffs::Integer, src::KnetArray{S}, soffs::Integer, n::Integer) where {T,S} = copyto!(dest, doffs, convert(KnetArray{T},src), soffs, n)
copyto!(dest::KnetArray{T}, doffs::Integer, src::Array{S}, soffs::Integer, n::Integer) where {T,S} = copyto!(dest, doffs, convert(Array{T},src), soffs, n)
copyto!(dest::Array{T}, doffs::Integer, src::KnetArray{S}, soffs::Integer, n::Integer) where {T,S} = copyto!(dest, doffs, convert(KnetArray{T},src), soffs, n)
copyto!(dest::KnetArray{T}, src::KnetArray{S}) where {T,S} = copyto!(dest, convert(KnetArray{T},src))
copyto!(dest::KnetArray{T}, src::Array{S}) where {T,S} = copyto!(dest, convert(Array{T},src))
copyto!(dest::Array{T}, src::KnetArray{S}) where {T,S} = copyto!(dest, convert(KnetArray{T},src))

# This will make deepcopy work properly
copy(a::KnetArray) = copyto!(similar(a), a)
deepcopy_internal(x::KnetArray, s::IdDict) = (haskey(s,x) ? s[x] : copy(x))


text:
#ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля, на котором могут находиться также внутренние прямоугольные перегородки (все перегородки изолированы друг от друга, прямоугольники могут вырождаться в отрезки)

#РЕЗУЛЬТАТ: Робот - в исходном положении и в углах поля стоят маркеры

function main_glavn!(r::Robot)
    path = leftDown(r) #Робот в левом нижнем углу и массив path создан

    putmarker!(r)
    moveSide!(r, Nord)
    putmarker!(r)
    moveSide!(r, Ost)
    putmarker!(r)
    moveSide!(r, Sud)
    putmarker!(r)
    moveSide!(r, West)
    #Каждый угол помечен

    leftDown!(r, path) #Робот в исходном положении
end

function leftDown(r::Robot)
    path = []
    while isborder(r, West) == false || isborder(r, Sud) == false
          if isborder(r, West) == false
                move!(r, West)
                push!(path, true)
          else
                move!(r, Sud)
                push!(path, false)
          end
    end

    return reverse!(path)
end

function leftDown!(r::Robot, path)
    for i in path
          if i move!(r, Ost)
          else move!(r, Nord)
          end
    end
end

function moveSide!(r::Robot, side::HorizonSide)
    while isborder(r, side) == false
          move!(r, side)
    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

# simple 2d convolution with padding
function convolution(input, filter; padding=:same)
    input_r, input_c = size(input)
    filter_r, filter_c = size(filter)

    if padding == :same
        pad_r = (filter_r - 1) ÷ 2
        pad_c = (filter_c - 1) ÷ 2

        input_padded = zeros(input_r+(2*pad_r), input_c+(2*pad_c))
        for i in 1:input_r, j in 1:input_c
            input_padded[i+pad_r, j+pad_c] = input[i, j]
        end
        input = input_padded
        input_r, input_c = size(input)
    end

    result = zeros(input_r-filter_r+1, input_c-filter_c+1)
    result_r, result_c = size(result)

    for i in 1:result_r
        for j in 1:result_c
            for k in 1:filter_r
                for l in 1:filter_c
                    result[i,j] += input[i+k-1,j+l-1]*filter[k,l]
                end
            end
        end
    end

    return result
end

# gaussian kernel with same σ in both dimensions
function gaussian_kernel(σ::Real; l::Int = 4*ceil(Int,σ)+1)
    isodd(l) || throw(ArgumentError("length must be odd"))
    w = l>>1
    g = σ == 0 ? [exp(0/(2*oftype(σ, 1)^2))] : [exp(-x^2/(2*σ^2)) for x=-w:w]
    k = g/sum(g)
    return (k * k')
end

# gaussian kernel with different σs in both dimensions
function gaussian_kernel(
    σs::Tuple{Real, Real};
    l::Tuple{Int, Int} = (4*ceil(Int,σs[1])+1, 4*ceil(Int,σs[2])+1)
)
    all(isodd.(l)) || throw(ArgumentError("length must be odd"))
    w1 = l[1]>>1
    g1 = σs[1] == 0 ? [exp(0/(2*oftype(σs[1], 1)^2))] : [exp(-x^2/(2*σs[1]^2)) for x=-w1:w1]
    k1 = g1/sum(g1)

    w2 = l[2]>>1
    g2 = σs[2] == 0 ? [exp(0/(2*oftype(σs[2], 1)^2))] : [exp(-x^2/(2*σs[2]^2)) for x=-w2:w2]
    k2 = g2/sum(g2)

    return (k1 * k2')
end


text:
include("../common.jl")

build_libcurl(ARGS, "LibCURL")

# Build trigger: 0


text:
using ExcelFiles, DataFrames
using MetabolicFluxAnalysis

# conditions = excel sheet names
conditions = ["glucose","malate","malate-glucose","glycerol","pyruvate","gluconate","aspartate","glutamate succinate","fructose"]

map_fragments = Dict{String,Tuple{String,Vector{Integer}}}()

map_fragments["ALA-15"]  = ("Ala",[])
map_fragments["ALA-57"]  = ("Ala",[1,2,3])
map_fragments["ALA-85"]  = ("Ala",[1,3])

map_fragments["ASP-15"]  = ("Asp",[1,2,3,4])
map_fragments["ASP-57"]  = ("Asp",[1,2,3,4])
map_fragments["ASP-85"]  = ("Asp",[])
map_fragments["ASP159"] = ("Asp",[1,2,4])
map_fragments["ASP302"] = ("Asp",[1,2])

map_fragments["GLU-15"]  = ("Glu",[])
map_fragments["GLU-57"]  = ("Glu",[1,2,3,4,5])
map_fragments["GLU-85"]  = ("Glu",[1,2,4,5])
map_fragments["GLU159"]  = ("Glu",[1,2,4,5])
map_fragments["GLU302"]  = ("Glu",[1,2])

map_fragments["GLY-15"]  = ("Gly",[])
map_fragments["GLY-57"]  = ("Gly",[1,2])
map_fragments["GLY-85"]  = ("Gly",[2])

map_fragments["HIS-57"]  = ("His",[])
map_fragments["HIS159"] = ("His",[])

map_fragments["ILE-15"] = ("Ile",[])
map_fragments["ILE-85"] = ("Ile",[1,2,4,5,6])
map_fragments["ILE159"] = ("Ile",[1,2,4,5,6])

map_fragments["LEU-15"] = ("Leu",[])
map_fragments["LEU-85"] = ("Leu",[1,2,4,5,6])
map_fragments["LEU159"] = ("Leu",[1,2,4,5,6])

map_fragments["PHE-57"]  = ("Phe",[1,2,3,4,5,6,7,8,9])
map_fragments["PHE-85"]  = ("Phe",[1,2,3,4,5,6,7,9])
map_fragments["PHE159"]  = ("Phe",[1,2,3,4,5,6,7,9])
map_fragments["PHE302"]  = ("Phe",[])

map_fragments["PRO-85"]  = ("Pro",[])
map_fragments["PRO159"]  = ("Pro",[1,3,4,5])

map_fragments["SER-15"]  = ("Ser",[])
map_fragments["SER-57"]  = ("Ser",[1,2,3])
map_fragments["SER-85"]  = ("Ser",[1,3])
map_fragments["SER159"]  = ("Ser",[1,3])
map_fragments["SER302"]  = ("Ser",[])

map_fragments["THR-15"]  = ("Thr",[])
map_fragments["THR-57"]  = ("Thr",[1,2,3,4])
map_fragments["THR-85"]  = ("Thr",[])

map_fragments["TYR-15"]  = ("Tyr",[])
map_fragments["TYR-57"]  = ("Tyr",[])
map_fragments["TYR-85"]  = ("Tyr",[])
map_fragments["TYR159"]  = ("Tyr",[])

map_fragments["VAL-15"]  = ("Val",[])
map_fragments["VAL-57"]  = ("Val",[1,2,3,4,5])
map_fragments["VAL-85"]  = ("Val",[1,2,3,4])
map_fragments["VAL159"]  = ("Val",[1,2,3,4])
map_fragments["VAL302"]  = ("Val",[])

exp_m_glucose = Dict{Int64,String}()
exp_m_glucose[1] = "e1m1"
exp_m_glucose[2] = "e1m2"
exp_m_glucose[3] = "e2m1"
exp_m_glucose[4] = "e2m2"
exp_m_glucose[5] = "e2m3"
exp_m_mal = Dict{Int64,String}()
exp_m_mal[1] = "e1m1"
exp_m_mal[2] = "e1m2"
exp_m_mal[3] = "e2m1"
exp_m_mal[4] = "e2m2"
exp_m_mal[5] = "e2m3"
exp_m_malgluc = Dict{Int64,String}()
exp_m_malgluc[1] = "e1m1"
exp_m_malgluc[2] = "e1m2"
exp_m_malgluc[3] = "e2m1"
exp_m_malgluc[4] = "e2m2"
exp_m_malgluc[5] = "e2m3"
exp_m_glycerol = Dict{Int64,String}()
exp_m_glycerol[1] = "e1m1"
exp_m_glycerol[2] = "e1m2"
exp_m_glycerol[3] = "e2m1"
exp_m_glycerol[4] = "e2m2"
exp_m_glycerol[5] = "e2m3"
exp_m_pyruvate = Dict{Int64,String}()
exp_m_pyruvate[1] = "e1m1"
exp_m_pyruvate[2] = "e1m2"
exp_m_pyruvate[3] = "e2m1"
exp_m_pyruvate[4] = "e2m2"
exp_m_pyruvate[5] = "e2m3"
exp_m_pyruvate[6] = "e3m1"
exp_m_pyruvate[7] = "e3m2"
exp_m_pyruvate[8] = "e3m3"
exp_m_pyruvate[9] = "e4m1"
exp_m_pyruvate[10] = "e4m2"
exp_m_pyruvate[11] = "e4m3"
exp_m_gluconate = Dict{Int64,String}()
exp_m_gluconate[1] = "e1m1"
exp_m_gluconate[2] = "e1m2"
exp_m_gluconate[3] = "e2m1"
exp_m_gluconate[4] = "e2m2"
exp_m_gluconate[5] = "e2m3"
exp_m_aspartate = Dict{Int64,String}()
exp_m_aspartate[1] = "e1m1"
exp_m_aspartate[2] = "e1m2"
exp_m_aspartate[3] = "e1m3"
exp_m_aspartate[4] = "e2m1"
exp_m_aspartate[5] = "e2m2"
exp_m_aspartate[6] = "e2m3"
exp_m_glutamsuccin = Dict{Int64,String}()
exp_m_glutamsuccin[1] = "e1m1"
exp_m_glutamsuccin[2] = "e1m2"
exp_m_glutamsuccin[3] = "e1m3"
exp_m_glutamsuccin[4] = "e2m1"
exp_m_glutamsuccin[5] = "e2m2"
exp_m_glutamsuccin[6] = "e2m3"
exp_m_fructose = Dict{Int64,String}()
exp_m_fructose[1] = "e1m1"
exp_m_fructose[2] = "e1m2"
exp_m_fructose[3] = "e2m1"
exp_m_fructose[4] = "e2m2"
exp_m_fructose[5] = "e2m3"

conditions = ["glucose","malate","malate-glucose","glycerol","pyruvate","gluconate","aspartate","glutamate succinate","fructose"]

map_measurements = Dict{String,Dict{Int64,String}}()
map_measurements["glucose"] = exp_m_glucose
map_measurements["malate"] = exp_m_mal
map_measurements["malate-glucose"] = exp_m_malgluc
map_measurements["glycerol"] = exp_m_glycerol
map_measurements["pyruvate"] = exp_m_pyruvate
map_measurements["gluconate"] = exp_m_gluconate
map_measurements["aspartate"] = exp_m_aspartate
map_measurements["glutamate succinate"] = exp_m_glutamsuccin
map_measurements["fructose"] = exp_m_fructose

# parse everything:

measurements = Dict{ String, Vector{Tuple{SpeciesEMU,Tuple{Int64,Array{Float64}}}} }()

for zi in 1:length(conditions)

    condition_i = conditions[zi]
    measurements_i = map_measurements[ condition_i ]

    df = DataFrame(load("C:\\Users\\Thomas\\Desktop\\chubukov_data\\chubukov MDVs.xls", condition_i))

    # process the fragment mass description:
    list_of_fragments = Vector{Tuple{SpeciesEMU,Tuple{Int64,Array{Float64}}}}()

    for zj in 6:256
        str_frag = df[1][zj]
        str_frag_split = split(str_frag," ")
        frag = map_fragments[str_frag_split[1]]
        if(isempty(frag[2]))
            continue # fragment not clear..
        end

        str_mass = replace(str_frag_split[2],"m","")
        i_mass = parse(Int64,str_mass)

        m_values = Vector{Float64}()
        for ki in keys(measurements_i)
            # pos is (ki+1)
            push!( m_values , df[ki+1][zj] )
        end
        push!( list_of_fragments , ( constructSpeciesEMU(frag[1],frag[2]) , (i_mass,m_values) ) )
    end

    measurements[condition_i] = list_of_fragments
end

v1 = Vector([1,2,3])


text:
using Test
using QuantumOptics
using LinearAlgebra
import StochasticDiffEq

@testset "stochastic_schroedinger" begin

b_spin = SpinBasis(1//2)
sz = sigmaz(b_spin)
sm = sigmam(b_spin)
sp = sigmap(b_spin)
zero_op = 0*sz
γ = 0.5
noise_op = 0.5γ*sz

H = γ*(sp + sm)
Hs = [noise_op]

ψ0 = spindown(b_spin)
ρ0 = dm(ψ0)

dt = 1/30.0
T = [0:1.0:100;]
T_short = [0:dt:dt;]

# Test equivalence of stochastic schroedinger phase noise and master dephasing
Ntraj = 100
ρ_avg = [0*ρ0 for i=1:length(T)]
for i=1:Ntraj
    t, ψt = stochastic.schroedinger(T, ψ0, H, Hs; dt=dt,
        alg=StochasticDiffEq.EulerHeun())
    ρ_avg += dm.(ψt)./Ntraj
end
tout, ρt = timeevolution.master(T, ρ0, H, [sz]; rates=[0.25γ^2])

for i=1:length(tout)
    @test tracedistance(ρ_avg[i], ρt[i]) < 5dt
end

# Function definitions for schroedinger_dynamic
function fdeterm(t, psi)
    H
end
function fstoch_1(t, psi)
    [zero_op]
end
function fstoch_2(t, psi)
    [zero_op, zero_op, zero_op]
end

# Non-dynamic Schrödinger
tout, ψt4 = stochastic.schroedinger(T, ψ0, H, [zero_op, zero_op]; dt=0.1dt)
# Dynamic Schrödinger
tout, ψt1 = stochastic.schroedinger_dynamic(T, ψ0, fdeterm, fstoch_1; dt=0.1dt)
tout, ψt2 = stochastic.schroedinger_dynamic(T, ψ0, fdeterm, fstoch_2; noise_processes=3, dt=0.1dt)

# Test equivalence to Schrödinger equation with zero noise
# Test sharp equality for same algorithms
@test ψt1 == ψt2 == ψt4

tout, ψt_determ = timeevolution.schroedinger_dynamic(T, ψ0, fdeterm)
# Test approximate equality for different algorithms
for i=1:length(tout)
    @test norm(ψt1[i] - ψt2[i]) < 3dt
    @test norm(ψt1[i] - ψt_determ[i]) < 10dt
end

# Test remaining function calls for short times to test whether they work in principle
tout, ψt1 = stochastic.schroedinger(T, ψ0, H, noise_op; dt=dt,
        normalize_state=true,
        alg=StochasticDiffEq.EulerHeun())
tout, ψt2 = stochastic.schroedinger_dynamic(T, ψ0, fdeterm, fstoch_1; dt=dt,
        normalize_state=true)

for i=1:length(T)
    @test norm(ψt1[i]) ≈ 1.0
    @test norm(ψt2[i]) ≈ 1.0
end

end # testset


text:
abstract type AbstractBounds end
abstract type AbstractConstBounds <: AbstractBounds end
abstract type AbstractVarBounds <: AbstractBounds end

# TODO: isa for type checking
# TODO: use promote
# Point(x::Real, y::Real) = Point(promote(x,y)...);
@with_kw struct ConstSymBounds{T<:Real} <: AbstractConstBounds
    hi::T
    lo::T
    function ConstSymBounds{T}(hi::T, lo::T) where T<:Real
        hi > zero(hi) || error("upper bound needs to be positive")
        lo < zero(lo) || error("lower bound needs to be negative")
        @assert isequal(hi, -lo)
        new(hi, lo)
    end
end

ConstSymBounds(b::Real) = ConstSymBounds(b, -b)
ConstSymBounds() = ConstSymBounds(1, -1)

@with_kw struct ConstAsymBounds{T<:Real} <: AbstractConstBounds
    hi::T
    lo::T
    function ConstAsymBounds{T}(hi::T, lo::T) where T<:Real
        hi > zero(hi) || error("upper bound needs to be positive")
        lo < zero(lo) || error("lower bound needs to be negative")
        new(hi, lo)
    end
end

@with_kw struct VarSymBounds <: AbstractVarBounds
    hi::Array{Float64,1}
    lo::Array{Float64,1}
    function VarSymBounds(hi::Array{Float64,1}, lo::Array{Float64,1})
        all(hi .> 0) || error("upper bound needs to be positive")
        all(lo .< 0) || error("lower bound needs to be negative")
        @assert isequal(hi, -lo)
        new(hi, lo)
    end
end

@with_kw struct VarAsymBounds <: AbstractVarBounds
    hi::Array{Float64,1}
    lo::Array{Float64,1}
    function VarAsymBounds(hi::Array{Float64,1}, lo::Array{Float64,1})
        all(hi .> 0) || error("upper bound needs to be positive")
        all(lo .< 0) || error("lower bound needs to be negative")
        new(hi, lo)
    end
end

function get(b::AbstractConstBounds)
    hi = [b.hi]
    lo = [b.lo]
    return (hi = hi, lo = lo)
end

function get(b::AbstractVarBounds)
    hi = b.hi
    lo = b.lo
    return (hi = hi, lo = lo)
end

# function decayingbound(type::String, b0::Real, κ::Real, Δt::Real, tmax::Real)
#       type ∈ ["upper", "lower"] || error("Bound type needs to be upper or lower.")
#       type = type == "upper" ? 1 : -1
#       x = range(0, tmax, step = Δt)
#       type .* (exp.(b0 .- κ.*x))
# end

@with_kw struct DecayingBound
    bound::Array{Float64,1}

    function DecayingBound(type::String, b0::Real, κ::Real, Δt::Real, tmax::Real)
        type ∈ ["upper", "lower"] || error("Bound type needs to be upper or lower.")
        type = type == "upper" ? 1 : -1
        x = range(0, tmax, step = Δt)
        return new(type .* (exp.(b0 .- κ.*x)))
    end
end

@with_kw struct DecayingBounds <: AbstractVarBounds
    hi::Array{Float64,1}
    lo::Array{Float64,1}
    function DecayingBounds(upper::DecayingBound, lower::DecayingBound)
        return new(upper.bound, lower.bound)
    end
end


text:
function test_backup()
    p = LightDark1D()
    policy = MCVIPolicy(p)
    MCVI.initialize_updater!(policy)
    b0 = MCVI.initial_belief(p, 500, MersenneTwister(42))
    # s0 = initialstate(p)
    # sim = MCVISimulator(MersenneTwister(420), s0, 1)
    scratch = MCVI.Scratch(Vector{Float64}(undef, 50), zeros(50), zeros(50), zeros(50, 2))
    sim = MCVISimulator()

    n, _ = MCVI.backup(b0, policy, sim, p, 500, 1000, 500, scratch)
    MCVI.addnode!(policy.updater, n)
    policy.updater.root = n

    v1 = MCVI.evaluate(policy, sim, p)
    n, _ = MCVI.backup(b0, policy, sim, p, 500, 1000, 500, scratch)
    MCVI.addnode!(policy.updater, n)
    policy.updater.root = n

    v2 = MCVI.evaluate(policy, sim, p)
    return (v1 - 0.1 < v2)
end


text:
#Work in progress

# Examples of constraint functions which can be used, however,
# \/\/\/these won't pack into ProtoBuf yet -- need to write special converters
# for Neo4j DataBase storage of complicated types. See constraints hereafter
# for more standard types.
type UniPriorPose2D <: Singleton
  Z::Distributions.MvNormal
end
type GMMPriorPose2D <: Singleton
  Z::Array{Distributions.MvNormal,1}
  W::Array{Float64,1}
end
type KDEPriorPoint2D <: Singleton
  Z::BallTreeDensity
end
type KDERangePoint2D <: Pairwise
  Z::BallTreeDensity
end
# ^^^will work on these later


text:
module demo
using CLI

@command """
this is a demo CLI
"""

module math
using CLI

@command """
math calculations
"""

"""
add up
"""
@command function add(x::Int, y::Int; extra::Int=1)
    x + y + extra
end

"""
multiply two integers
"""
@command mul(x::Int, y::Int) = x * y

end # math

"""
plus two integers
"""
@command function plus(x::Int, y::Int)
    x + y
end

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    return main()
end

end

demo.julia_main(ARGS)


text:
#=
Generalized versions of the JuMPeR interface
Delegates to correct cutting plane algorithm via suppFcn
=#
typealias DDUSOracles Union{UMOracle, UIOracle, FBOracle, LCXOracle, UCSOracle}

# JuMPeR alerting us that we're handling this contraint
registerConstraint(w::DDUSOracles, rm::Model, ind::Int, prefs) =
    !get(prefs, :prefer_cuts, true) && error("Only cutting plane supported")

# JuMPeR wants us to generate a constraint for every uncertain constraint in inds
function generateCut(w::DDUSOracles, m::Model, rm::Model, inds::Vector{Int}, active=false)
    new_cons = Any[]
    rd = JuMPeR.getRobust(rm)
    for ind in inds
        con = JuMPeR.get_uncertain_constraint(rm, ind)
        cut_sense, xs, lhs_const = JuMPeR.build_cut_objective(rm, con, m.colVal)
        zstar, ustar = suppFcn(xs, w, cut_sense)
        lhs_of_cut = zstar + lhs_const

        # SUBJECT TO CHANGE: active cut detection
        if active
            push!(rd.activecuts[ind],
                JuMPeR.cut_to_scen(ustar,
                    JuMPeR.check_cut_status(con, lhs_of_cut, w.cut_tol) == :Active))
            continue
        end

        # Check violation
        if JuMPeR.check_cut_status(con, lhs_of_cut, w.cut_tol) != :Violate
            w.debug_printcut && JuMPeR.debug_printcut(rm ,m,w,lhs_of_cut,con,nothing)
            continue  # No violation, no new cut
        end

        # Create and add the new constraint
        new_con = JuMPeR.build_certain_constraint(m, con, ustar)
        w.debug_printcut && JuMPeR.debug_printcut(rm, m, w, lhs_of_cut, con, new_con)
        push!(new_cons, new_con)
    end
    return new_cons
end

# JuMPeR asking us for any reformulations we might want to make - we make none
generateReform(w::DDUSOracles, m::Model, rm::Model, inds::Vector{Int}) = 0


text:
# This file is part of GenericSchur.jl, released under the MIT "Expat" license

# The methods in this file are derived from LAPACK's ztrsyl.
# LAPACK is released under a BSD license, and is
# Copyright:
# Univ. of Tennessee
# Univ. of California Berkeley
# Univ. of Colorado Denver
# NAG Ltd.

function trsylvester!(A::UpperTriangular{T,S}, B::UpperTriangular{T,S},
                      C::StridedVecOrMat{T};
                      possign=false) where {T,S<:StridedMatrix{T}}
    trsylvester!(A.data,B.data,C,possign=possign)
end

"""
    trsylvester!(A,B,C) => X, σ

solve the Sylvester equation ``A X - X B = σ C``
for upper triangular and square `A` and `B`, overwriting `C`
and setting `σ` to avoid overflow.
"""
function trsylvester!(A::StridedMatrix{T},B::StridedMatrix{T},C::StridedVecOrMat{T}; possign=false) where {T}
    m = checksquare(A)
    n = checksquare(B)
    ((size(C,1) == m) && (size(C,2) == n)) || throw(DimensionMismatch(
        "dimensions of C $(size(C)) must match A, ($m,$m), and B, ($n,$n)"))

    scale = one(real(T))

    tiny = eps(real(T))
    small = safemin(real(T)) * m * n / tiny
    bignum = one(real(T)) / small
    smin = max(small, tiny * norm(A,Inf) * norm(B,Inf))

    isgn = possign ? one(real(T)) : -one(real(T))
    ierr = 0
    for l=1:n
        for k=m:-1:1
            # FIXME: these should be dotu(), but I can't find anything usable in stdlib::LA
            # WARNING: this relies on isempty() pass in checkindex() for view()
            suml = sum(view(A,k,k+1:m) .* view(C,k+1:m,l))
            sumr = sum(view(C,k,1:l-1) .* view(B,1:l-1,l))
            v = C[k,l] - (suml + isgn * sumr)
            scaloc = one(real(T))
            a11 = A[k,k] + isgn * B[l,l]
            da11 = abs1(a11)
            if da11 <= smin
                a11 = smin
                da11 = smin
                ier = 1
            end
            db = abs1(v)
            if (da11 < 1) && (db > 1)
                if (db > bignum * da11)
                    # println("scaling by $db")
                    scaloc = 1 / db
                end
            end
            x11 = (v * scaloc) / a11
            if scaloc != 1
                lmul!(scaloc, C)
                scale *= scaloc
            end
            C[k,l] = x11
        end
    end
    return C, scale
end

function adjtrsylvester!(A::UpperTriangular{T,S}, B::UpperTriangular{T,S},
                         C::StridedVecOrMat{T};
                         possign=false) where {T,S<:StridedMatrix{T}}
    adjtrsylvester!(A.data,B.data,C,possign=possign)
end

"""
    adjtrsylvester!(A,B,C) => X, σ

solve the Sylvester equation ``Aᴴ X - X Bᴴ = σ C``,
for upper triangular `A` and `B`, overwriting `C`
and setting `σ` to avoid overflow.
"""
function adjtrsylvester!(A::StridedMatrix{T},B::StridedMatrix{T},C::StridedVecOrMat{T}; possign=false) where {T}
    m = checksquare(A)
    n = checksquare(B)
    ((size(C,1) == m) && (size(C,2) == n)) || throw(DimensionMismatch(
        "dimensions of C $(size(C)) must match A, ($m,$m), and B, ($n,$n)"))

    scale = one(real(T))

    tiny = eps(real(T))
    small = safemin(real(T)) * m * n / tiny
    bignum = one(real(T)) / small
    smin = max(small, tiny * norm(A,Inf) * norm(B,Inf))

    isgn = possign ? one(real(T)) : -one(real(T))
    for l=n:-1:1
        for k=1:m
            suml = dot(view(A,1:k-1,k),view(C,1:k-1,l))
            # WARNING: this relies on isempty() pass in checkindex() for view()
            sumr = dot(view(C,k,l+1:n),view(B,l,l+1:n))
            v = C[k,l] - (suml + isgn * conj(sumr))
            scaloc = one(real(T))
            a11 = conj(A[k,k] + isgn * B[l,l])
            da11 = abs1(a11)
            if da11 < smin
                a11 = smin
                da11 = smin
                ier = 1
            end
            db = abs1(v)
            if (da11 < 1) && (db > 1)
                if (db > bignum * da11)
                    scaloc = 1 / db
                end
            end
            x11 = (v * scaloc) / a11
            if scaloc != 1
                lmul!(scaloc, C)
                scale *= scaloc
            end
            C[k,l] = x11
        end
    end
    return C, scale
end


text:
"""
    compactify(rational_to_compactify, radius_of_indifference)

    compactify(low=one_inteval_bound, high=other_interval_bound)

From all of the rational values that exist within ±`radius_of_indifference`
of the `rational_to_compactify`, the interval that obtains is uniquely determines
that rational for which (a) the denominator is that of least magnitude and
(b) the numerator is either uniquely given or, of those available, of least magnitude.

The first form suffers from inaccuracies due to rounding errors.

We are indifferent to the two rational values, source and result, as magnitudes.
We prefer to use the compactified value in calculations, as with it, overflow
is less likely, probably, with the next arithmetic operation.
"""
function compactify(;low::Q, high::Q) where {Q<:Rationals}
    if low < high <= 0 || 0 <= high < low
        low, high = high, low
    end
    lonum, loden, hinum, hiden = low.num, low.den, high.num, high.den
    if lonum == 0 || hinum == 0 || lonum < 0 == hinum > 0
        return zero(Q)
    end
    lonumsbit = signbit(lonum)
    if lonumsbit
        lonum, hinum = -lonum, -hinum
    end
    num, den = compact_integer(lonum, loden, hinum, hiden)
    return Q(lonumsbit ? -num : num, den)
end

# compute Q(midpoint - rad) and Q(midpoint + rad), approximately
# preserving the underlying integer type

function compact_reduce(midpoint::Q, radius::Q) where {Q<:Rationals}
    midnum, midden = midpoint.num, midpoint.den
    radnum, radden = radius.num, radius.den
    u = typemax(basetype(Q)) ÷ max(midden, midnum)
    # either `abs(au)` or `abs(bu)` is close to typemax of the underlying integer type
    au = midnum * u
    bu = midden * u
    v = (bu ÷ radden) * radnum
    # (au - v) // bu  <=  midpoint - rad  <  (au - v + 1) // bu
    lo3 = compact_integer(au - v, bu, au - v + 1, bu)
    # (au + v - 1) // bu  <=  midpoint + rad  <  (au + v) // bu
    hi3 = compact_integer(au + v - 1, bu, au + v, bu)
    tuple(lo3..., hi3...)
end

for (Q1, Q2) in ((:FastRational, :FastRational), (:FastRational, :Rational), (:Rational, :Rational))
  @eval begin
    function compactify(midpoint::$Q1{T}, radius::$Q2{T}) where {T<:Integer}
        mid, rad = float(midpoint), float(abs(radius))
        lo, hi = mid-rad, mid+rad
        lo, hi = compact_rational_constraints(mid, rad, lo, hi)
        num, den = T.(compact_rational(lo, hi))
        return $Q1{T}(num, den)
    end
  end
end

for Q in (:FastRational, :Rational)
  @eval begin
    function compactify(midpoint::$Q{T}, radius::F) where {T<:Integer, F<:Union{Float32, Float64}}
        mid = F(midpoint)
        lo, hi = mid-radius, mid+radius
        lo, hi = compact_rational_constraints(mid, radius, lo, hi)
        num, den = T.(compact_rational(lo, hi))
        return $Q{T}(num, den)
    end
  end
end

#=
    The algorithm for `compact_rational` is from Hiroshi Murakami's paper
    Calculation of rational numbers in an interval whose denominator is the smallest
    ACM Communications in Computer Algebra, Vol 48, No. 3, Issue 189, September 2014
    Hiroshi Murakami, Department of Mathematics and Information Sciences,
    Tokyo Metropolitan University, Tokyo, 192-0397, Japan
    mrkmhrsh@tmu.ac.jp, https://www.rs.tus.ac.jp/hide-murakami/index.html
=#
function compact_rational(lo::T, hi::T) where {T<:Real}
    if ceil(lo) <= floor(hi)                           # [lo,hi] contains some integer
        num, den = ceil(lo), one(T)                    # the CF expansion terminates here.
    else                                               # [lo,hi] contains no integer
        m = floor(lo)                                  #
        lo, hi = inv(hi-m), inv(lo-m)                  # the CF expansion continues.
        num, den = compact_rational(lo, hi)  # Recursive call is made here.
        num, den = num*m + den, num
    end
    return num, den
end

# implementation of previous algorithm using integer arithmetic
function compact_integer(a::T, b::T, c::T, d::T) where T<:Integer
    m, r = divrem(a, b)
    if r == 0
        m, one(T)
    else
        s = c - d * m
        if s >= d
            m+one(T), one(T)
        else
            num, den = compact_integer(d, s, b, r)
            num * m + den, num
        end
    end
end

function compact_rational(x::Q, y::Q) where {Q<:Rationals}
   compact_integer(x.num, x.den, y.num, y.den)
end

# convert input arguments to float and check for argument errors
function compact_rational_constraints(mid::T, rad::T, lo::T, hi::T) where {T<:Real}
    rad < eps(mid) && throw(ArgumentError("radius is less than eps(float(midpoint)))"))
    lo, hi = abs(lo), abs(hi)
    lo, hi = lo < hi ? (lo, hi) : (hi, lo)
    !iszero(lo) || throw(ArgumentError("lo == 0"))
    flo, fhi = float(lo), float(hi)
    return flo, fhi
end

#=
from the paper (section 1)
For a given positive real interval I = [α,β] (0 < α < β),
    we ﬁnd a rational number P/Q in I
        whose positive denominator Q is the smallest.
We may change the problem to:
    “For the interval I, ﬁnd (all) rational numbers in I whose positive denominator is the smallest”.
We may also change the problem to:
    “For the interval I, from all rational numbers in I whose positive denominator is the smallest,
        ﬁnd the one whose numerator is the smallest”,
        because when we ﬁnd Q the smallest positive denominator,
        the numerator P of the original problem is any integer in
        ceil(Qα),···,ﬂoor(Qβ) contained in [Qα, Qβ].
The rational number P/Q is uniquely determined from the real positive closed interval I [α,β] (0<α<β)
   can be found by the CF expansion method[1]. This algorithm is an analogy of the
   regular continued fraction expansion of a real number, and gives a solution very quickly
   even when the denominator of the solution is a huge integer. The count of arithmetic operations
   is linear to the bit length of the denominator of the solution.
=#


text:
mutable struct GFUngenerateState
    trace::GFTrace
    constraints::Any
    weight::Float64
    visitor::AddressVisitor
    params::Dict{Symbol,Any}
end

function GFUngenerateState(trace::GFTrace, constraints, params::Dict{Symbol,Any})
    GFUngenerateState(trace, constraints, 0., AddressVisitor(), params)
end

function addr(state::GFUngenerateState, dist::Distribution{T}, args, addr) where {T}
    visit!(state, addr)
    call::CallRecord = get_primitive_call(state.trace, addr)
    @assert call.args == args
    retval::T = call.retval
    if has_leaf_node(state.trace, addr)
        if retval != get_leaf_node(state.trace, addr)
            error("constraints and trace do not match at $addr")
        end
        state.weight += retval.score
    elseif has_internal_node(state.trace, addr)
        lightweight_got_internal_node_err(addr)
    end
    retval
end

function addr(state::GFUngenerateState, gen::Generator{T}, args, addr) where {T}
    visit!(state, addr)
    subtrace = get_subtrace(state.trace, addr)
    if has_internal_node(state.constraints, addr)
        constraints = get_internal_node(state.constraints, addr)
    elseif has_leaf_node(state.constraints, addr)
        lightweight_got_leaf_node_err(addr)
    else
        constraints = EmptyAssignment()
    end
    weight = ungenerate(gen, subtrace, constraints)
    call::CallRecord = get_call_record(subtrace)
    call.retval::T
end

function splice(state::GFUngenerateState, gf::GenFunction, args::Tuple)
    exec(gf, state, args)
end

function ungenerate(gf::GenFunction, trace::GFTrace, constraints)
    state = GFUngenerateState(trace, constraints, gf.params)
    exec(gf, state, args)
    state.weight
end

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
using Statistics
using Oceananigans.Architectures: arch_array
using Oceananigans.Fields: ReducedField, CenterField, ZFaceField, compute_at!, @compute
using Oceananigans.BoundaryConditions: fill_halo_regions!
using Oceananigans.Grids: halo_size

include("utils_for_runtests.jl")

archs = test_architectures()

"""
    correct_reduced_field_size(loc, arch, grid, dims, Tx, Ty, Tz)

Test that the ReducedField at `loc`ation on `arch`itecture and `grid`
and reduced along `dims` has size `(Tx, Ty, Tz)`.
"""
correct_reduced_field_size(loc, arch, grid, dims, Tx, Ty, Tz) =
    size(parent(ReducedField(loc, arch, grid; dims=dims))) == (Tx, Ty, Tz)

function correct_reduced_field_value_was_set(arch, grid, loc, dims, val::Number)
    f = ReducedField(loc, arch, grid; dims=dims)
    set!(f, val)
    return all(interior(f) .≈ val * arch_array(arch, ones(size(f))))
end

@testset "AbstractReducedFields" begin
    @info "Testing AbstractReducedFields..."

    N = (4, 6, 8)
    L = (2π, 3π, 5π)
    H = (1, 1, 1)

    @testset "ReducedField initialization" begin
        @info "  Testing ReducedField initialization..."
        for arch in archs, FT in float_types

            grid = RectilinearGrid(arch, FT, size=N, extent=L, halo=H, topology=(Bounded, Bounded, Bounded))

            @test correct_reduced_field_size((Center, Center, Center), arch, grid, 1,         1,               N[2] + 2 * H[2],     N[3] + 2 * H[3])
            @test correct_reduced_field_size((Face,   Center, Center), arch, grid, 1,         1,               N[2] + 2 * H[2],     N[3] + 2 * H[3])
            @test correct_reduced_field_size((Center, Face,   Center), arch, grid, 1,         1,               N[2] + 2 * H[2] + 1, N[3] + 2 * H[3])
            @test correct_reduced_field_size((Center, Face,   Face),   arch, grid, 1,         1,               N[2] + 2 * H[2] + 1, N[3] + 2 * H[3] + 1)
            @test correct_reduced_field_size((Center, Center, Center), arch, grid, 2,         N[1] + 2 * H[1], 1,                   N[3] + 2 * H[3])
            @test correct_reduced_field_size((Center, Center, Center), arch, grid, 2,         N[1] + 2 * H[1], 1,                   N[3] + 2 * H[3])
            @test correct_reduced_field_size((Center, Center, Center), arch, grid, 3,         N[1] + 2 * H[1], N[2] + 2 * H[2],     1)
            @test correct_reduced_field_size((Center, Center, Center), arch, grid, (1, 2),    1,               1,                   N[3] + 2 * H[3])
            @test correct_reduced_field_size((Center, Center, Center), arch, grid, (2, 3),    N[1] + 2 * H[1], 1,                   1)
            @test correct_reduced_field_size((Center, Center, Center), arch, grid, (1, 2, 3), 1,               1,                   1)
        end
    end

    reduced_dims = (1, 2, 3, (1, 2), (2, 3), (1, 3), (1, 2, 3))

    int_vals = Any[0, Int8(-1), Int16(2), Int32(-3), Int64(4)]
    uint_vals = Any[6, UInt8(7), UInt16(8), UInt32(9), UInt64(10)]
    float_vals = Any[0.0, -0.0, 6e-34, 1.0f10]
    rational_vals = Any[1//11, -23//7]
    other_vals = Any[π]
    vals = vcat(int_vals, uint_vals, float_vals, rational_vals, other_vals)

    @testset "Setting ReducedFields" begin
        @info "  Testing ReducedField setting..."
        for arch in archs, FT in float_types

            grid = RectilinearGrid(arch, FT, size=N, extent=L, halo=H, topology=(Periodic, Periodic, Bounded))

            for dims in reduced_dims, val in vals
                @test correct_reduced_field_value_was_set(arch, grid, (Center, Center, Center), dims, val)
            end

            for dims in reduced_dims
                reduced_field = ReducedField((Center, Center, Center), arch, grid, dims=dims)
                sz = size(reduced_field)
                A = rand(FT, sz...)
                set!(reduced_field, A)

                @test reduced_field[1, 1, 1] == A[1, 1, 1]

                fill_halo_regions!(reduced_field, arch)

                # No-flux boundary conditions at top and bottom
                @test reduced_field[1, 1, 0] == A[1, 1, 1]
                @test reduced_field[1, 1, grid.Nz+1] == A[1, 1, end]

                # Periodic boundary conditions in the horizontal directions
                @test reduced_field[1, 0, 1] == A[1, end, 1]
                @test reduced_field[1, grid.Ny+1, 1] == A[1, 1, 1]

                @test reduced_field[0, 1, 1] == A[end, 1, 1]
                @test reduced_field[grid.Nx+1, 1, 1] == A[1, 1, 1]
            end
        end
    end

    for arch in archs
        arch_str = string(typeof(arch))

        grid = RectilinearGrid(arch, topology = (Periodic, Periodic, Bounded),
                                      size = (2, 2, 2),
                                      x = (0, 2),
                                      y = (0, 2),
                                      z = (0, 2))

        w = ZFaceField(arch, grid)
        T = CenterField(arch, grid)

        trilinear(x, y, z) = x + y + z

        set!(T, trilinear)
        set!(w, trilinear)

        @compute Txyz = AveragedField(T, dims=(1, 2, 3))

        # Note: halo regions must be *filled* prior to computing an average
        # if the average within halo regions is to be correct.
        fill_halo_regions!(T, arch)
        @compute Txy = AveragedField(T, dims=(1, 2))

        fill_halo_regions!(T, arch)
        @compute Tx = AveragedField(T, dims=1)

        @compute wxyz = AveragedField(w, dims=(1, 2, 3))
        @compute wxy = AveragedField(w, dims=(1, 2))
        @compute wx = AveragedField(w, dims=1)

        Nx, Ny, Nz = grid.Nx, grid.Ny, grid.Nz

        @testset "Averaged fields [$arch_str]" begin
            @info "  Testing AveragedFields [$arch_str]"
            for FT in float_types

                @test Txyz[1, 1, 1] ≈ 3

                @test Array(interior(Txy))[1, 1, :] ≈ [2.5, 3.5]
                @test Array(interior(Tx))[1, :, :] ≈ [[2, 3] [3, 4]]

                @test wxyz[1, 1, 1] ≈ 3

                @test Array(interior(wxy))[1, 1, :] ≈ [2, 3, 4]
                @test Array(interior(wx))[1, :, :] ≈ [[1.5, 2.5] [2.5, 3.5] [3.5, 4.5]]

                # Test whether a race condition gets hit for averages over large fields
                big_grid = RectilinearGrid(arch, topology = (Periodic, Periodic, Bounded),
                                                  size = (256, 256, 128),
                                                  x = (0, 2),
                                                  y = (0, 2),
                                                  z = (0, 2))

                c = CenterField(arch, big_grid)
                c .= 1

                C = AveragedField(c, dims=(1, 2))

                # Test that the mean consistently returns 1 at every z for many evaluations
                results = [all(interior(mean!(C, C.operand)) .== 1) for i = 1:10] # warm up...
                results = [all(interior(mean!(C, C.operand)) .== 1) for i = 1:10] # the real deal
                @test mean(results) == 1.0
            end
        end

        @testset "Allocating reductions [$arch_str]" begin
            @info "  Testing allocating reductions"

            # Mean
            @test Txyz[1, 1, 1] == mean(T)
            @test interior(Txy) == interior(mean(T, dims=(1, 2)))
            @test interior(Tx) == interior(mean(T, dims=1))

            @test wxyz[1, 1, 1] == mean(w)
            @test interior(wxy) == interior(mean(w, dims=(1, 2)))
            @test interior(wx) == interior(mean(w, dims=1))

            # Maximum and minimum
            @test maximum(T) == maximum(interior(T))
            @test minimum(T) == minimum(interior(T))

            for dims in (1, (1, 2))
                @test interior(minimum(T; dims)) == minimum(interior(T); dims)
                @test interior(minimum(T; dims)) == minimum(interior(T); dims)
            end
        end

        @testset "Conditional computation of AveragedFields [$(typeof(arch))]" begin
            @info "  Testing conditional computation of AveragedFields [$(typeof(arch))]"
            for FT in float_types
                grid = RectilinearGrid(arch, FT, size=(2, 2, 2), extent=(1, 1, 1))
                c = CenterField(arch, grid)

                for dims in (1, 2, 3, (1, 2), (2, 3), (1, 3), (1, 2, 3))
                    C = AveragedField(c, dims=dims)

                    @test !isnothing(C.status)

                    # Test conditional computation
                    set!(c, 1)
                    compute_at!(C, FT(1)) # will compute
                    @test all(interior(C) .== 1)
                    @test C.status.time == FT(1)

                    set!(c, 2)
                    compute_at!(C, FT(1)) # will not compute because status == 1
                    @test C.status.time == FT(1)
                    @test all(interior(C) .== 1)

                    compute_at!(C, FT(2))
                    @test C.status.time == FT(2)
                    @test all(interior(C) .== 2)
                end
            end
        end
    end
end


text:
using OAuth

consumer_key = ENV["OAUTH_CONSUMER_KEY"]
consumer_secret = ENV["OAUTH_CONSUMER_SECRET"]
access_token = ENV["OAUTH_TOKEN"]
access_token_secret = ENV["OAUTH_TOKEN_SECRET"]

get_oauth(endpoint::String, options::Dict) = oauth_request_resource(endpoint, "GET", options, consumer_key,
consumer_secret, access_token, access_token_secret)

tweet_oauth(endpoint::String, options::Dict) = oauth_request_resource(endpoint,"POST", options, consumer_key, consumer_secret,
access_token, access_token_secret)

max_stars_repo_name: PseudoCodeNerd/julia-memebot-twitter

max_issues_repo_name: PseudoCodeNerd/julia-memebot-twitter

max_forks_repo_name: PseudoCodeNerd/julia-memebot-twitter


text:
function build_robot!(vis, model::Box;
    i=1,
    r=0.1,
    tl=1.0,
    box_color=Colors.RGBA(0.0, 0.0, 0.0, tl))

    setobject!(vis["box_$i"], GeometryBasics.Rect(
        Vec(-1.0 * r, -0.005, -1.0 * r),
        Vec(2.0 * r, 0.01, 2.0 * r)),
        MeshPhongMaterial(color=box_color))
end

function set_robot!(vis, model::Box, q; i=1)
    settransform!(vis["box_$i"],
        compose(Translation(q[1], 0.001 * i, q[2]), LinearMap(RotY(-1.0 * q[3]))))
end

function visualize!(vis, model::Box, q;
    i=1,
    r=0.1,
    tl=1.0,
    box_color=Colors.RGBA(0.0, 0.0, 0.0, tl),
    Δt=0.1,
    fixed_camera=true)

    default_background!(vis)

    build_robot!(vis, model,
        i=i,
        r=r,
        tl=tl,
        box_color= box_color)

    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

    T = length(q)
    for t = 1:T-1
        MeshCat.atframe(anim, t) do
            set_robot!(vis, model, q[t], i=i)
        end
    end

    if fixed_camera
        MeshCat.settransform!(vis["/Cameras/default"],
        MeshCat.compose(MeshCat.Translation(0.0, -50.0, -1.0), MeshCat.LinearMap(Rotations.RotZ(-pi / 2.0))))
        setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 25)
    end

    MeshCat.setanimation!(vis, anim)
end


text:
module CSL

end


text:
using DimensionalData, Test, Unitful
using DimensionalData: slicedims, basetypeof, formatdims, modetype, name

@dim TestDim "Testname"

@testset "dims creation macro" begin
    @test TestDim(1:10, Sampled(), NoMetadata()) == TestDim(1:10, Sampled(), NoMetadata())
    @test TestDim(1:10; mode=Categorical()) == TestDim(1:10, Categorical(), NoMetadata())
    @test DimensionalData.name(TestDim) == :Testname
    @test label(TestDim) == "Testname"
    @test val(TestDim(:testval)) == :testval
    @test metadata(TestDim(1, AutoMode(), Metadata(a=1))) == Metadata(a=1)
    @test units(TestDim) == nothing
    @test label(TestDim) == "Testname"
    @test eltype(TestDim(1)) == Int
    @test eltype(TestDim([1, 2, 3])) <: Int
    @test length(TestDim(1)) == 1
    @test length(TestDim([1, 2, 3])) == 3
    @test eachindex(TestDim(10:10:30)) == Base.OneTo(3)
    @test eachindex(TestDim(10)) == Base.OneTo(1)
    @test step(TestDim(1:2:3)) == 2
    @test step(TestDim(Base.OneTo(10), mode=NoIndex())) == 1
    @test step(TestDim(1:2:3; mode=Sampled(span=Regular(2)))) == 2
    @test first(TestDim(5)) == 5
    @test last(TestDim(5)) == 5
    @test first(TestDim(10:20)) == 10
    @test last(TestDim(10:20)) == 20
    @test firstindex(TestDim(5)) == 1
    @test lastindex(TestDim(5)) == 1
    @test firstindex(TestDim(10:20)) == 1
    @test lastindex(TestDim(10:20)) == 11
    @test eachindex(TestDim(10:20)) == Base.OneTo(11)
    @test size(TestDim(10:20)) == (11,)
    @test ndims(TestDim(10:20)) == 1
    @test ndims(TestDim(1)) == 0
    @test Array(TestDim(10:15)) == [10, 11, 12, 13, 14, 15]
    @test iterate(TestDim(10:20)) == iterate(10:20)
    @test TestDim(5.0:7.0)[2] == 6.0
    @test TestDim(10)[] == 10
    @test TestDim(5.0:7.0, Sampled(Ordered(), Regular(1.0), Points()))[At(6.0)] == 6.0
end

@testset "formatdims" begin
    A = [1 2 3; 4 5 6]
    @test formatdims(A, (X, Y)) == (X(Base.OneTo(2), NoIndex(), NoMetadata()),
                                    Y(Base.OneTo(3), NoIndex(), NoMetadata()))
    @test formatdims(zeros(3), Ti) == (Ti(Base.OneTo(3), NoIndex(), NoMetadata()),)
    @test formatdims(A, (:a, :b)) == (Dim{:a}(Base.OneTo(2), NoIndex(), NoMetadata()),
                                      Dim{:b}(Base.OneTo(3), NoIndex(), NoMetadata()))
    @test formatdims(51:100, :c) == (Dim{:c}(Base.OneTo(50), NoIndex(), NoMetadata()),)
    @test formatdims(A, (a=[:A, :B], b=(10.0:10.0:30.0))) ==
    (Dim{:a}([:A, :B], Categorical(Ordered()), NoMetadata()),
     Dim{:b}(10.0:10:30.0, Sampled(Ordered(), Regular(10.0), Points()), NoMetadata()))
    @test formatdims(A, (X([:A, :B]; metadata=Metadata(a=5)),
           Y(10.0:10.0:30.0, Categorical(Ordered()), Metadata("metadata"=>1)))) ==
          (X([:A, :B], Categorical(Ordered()), Metadata(a=5)),
           Y(10.0:10:30.0, Categorical(Ordered()), Metadata("metadata"=>1)))
    @test formatdims(zeros(3, 4),
        (Dim{:row}(Val((:A, :B, :C))),
         Dim{:column}(Val((-20, -10, 0, 10)), Sampled(), NoMetadata()))) ==
        (Dim{:row}(Val((:A, :B, :C)), Categorical(), NoMetadata()),
         Dim{:column}(Val((-20, -10, 0, 10)), Sampled(Ordered(),Irregular(),Points()), NoMetadata()))
    @test formatdims(A, (X(:, Sampled(Ordered(), Regular(), Points())), Y)) ==
        (X(Base.OneTo(2), Sampled(Ordered(), Regular(1), Points()), NoMetadata()), Y(Base.OneTo(3), NoIndex(), NoMetadata()))
end

@testset "Val" begin
    @test dims(Val{X}()) == Val{X}()
    @test mode(Val{X}()) == NoIndex()
end

@testset "AnonDim" begin
    @test val(AnonDim()) == Colon()
    @test mode(AnonDim()) == NoIndex()
    @test metadata(AnonDim()) == NoMetadata()
    @test name(AnonDim()) == :Anon
end

@testset "Basic dim and array initialisation and methods" begin
    a = ones(5, 4)

    @test_throws DimensionMismatch DimArray(a, X)
    @test_throws DimensionMismatch DimArray(a, (X, Y, Z))

    da = DimArray(a, (X(LinRange(140, 148, 5)), Y(LinRange(2, 11, 4))))
    dimz = dims(da)

    for args in ((dimz,), (dimz, (X(), Y())), (dimz, X(), Y()),
                 (dimz, (X, Y)), (dimz, X, Y),
                 (dimz, (1, 2)), (dimz, 1, 2))
        @test val(args...) == index(args...) == (LinRange(140, 148, 5), LinRange(2, 11, 4))
        @test name(args...) == (:X, :Y)
        @test units(args...) == (nothing, nothing)
        @test label(args...) == ("X", "Y")
        @test sampling(args...) == (Points(), Points())
        @test span(args...) == (Regular(2.0), Regular(3.0))
        @test locus(args...) == (Center(), Center())
        @test order(args...) == (Ordered(), Ordered())
        @test arrayorder(args...) == order(ArrayOrder, args...) == (ForwardArray(), ForwardArray())
        @test indexorder(args...) == order(IndexOrder, args...) == (ForwardIndex(), ForwardIndex())
        @test relation(args...) == order(Relation, args...) == (ForwardRelation(), ForwardRelation())
        @test bounds(args...) == ((140, 148), (2, 11))
        @test mode(args...) == (Sampled(Ordered(), Regular(2.0), Points()),
                             Sampled(Ordered(), Regular(3.0), Points()))
    end

    @test val(dimz, ()) == index(dimz, ()) == ()
    @test val(dimz, 1) == val(dimz, X) == val(dimz, X()) == val(dimz[1])
    @test order(Relation, dimz, 1) == order(Relation, dimz, X) == order(Relation, dimz, X()) == ForwardRelation()
    @test order(Relation, dimz, ()) == ()

    @test dims(dimz, Y) === dimz[2]
    @test slicedims(dimz, (2:4, 3)) ==
        ((X(LinRange(142,146,3); mode=Sampled(Ordered(), Regular(2.0), Points())),),
             (Y(8.0, mode=Sampled(Ordered(), Regular(3.0), Points())),))

    @test mode(X) == NoIndex()
    @test modetype(dims(dimz, X)) == Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points}
end

@testset "arbitrary dim names and Val index" begin
    dimz = formatdims(zeros(3, 4),
           (Dim{:row}(Val((:A, :B, :C))),
            Dim{:column}(Val((-20, -10, 0, 10)), Sampled(Ordered(),Regular(10),Points()), NoMetadata()))
    )
    @test name(dimz) == (:row, :column)
    @test label(dimz) == ("row", "column")
    @test basetypeof(dimz[1]) == Dim{:row}
    @test length(dims(dimz, :row)) == 3
    @test size(dims(dimz, :row)) == (3,)
    @test eltype(dims(dimz, :row)) == Symbol
    @test firstindex(dimz[1]) == 1
    @test lastindex(dimz[1]) == 3
    @test ndims(dimz[1]) == 1
    @test Array(dimz[1]) == [:A, :B, :C]

    @testset "specify dim with Symbol" begin
        # @test_throws ArgumentError arrayorder(dimz, :x)
        # TODO Does this make sense?
        @test arrayorder(dimz, :row) == ForwardArray()
        @test arrayorder(dimz, :column) == ForwardArray()
        @test bounds(dimz, :row) == (nothing, nothing)
        @test bounds(dimz, :column) == (-20, 10)
        @test index(dimz, :row) == (:A, :B, :C)
        @test index(dimz, :column) == (-20, -10, 0, 10)
    end

end

@testset "repeating dims of the same type is allowed" begin
    dimz = X((143, 145)), X((-38, -35))
    @test dims(dimz, X) === dimz[1]
    @test dims(dimz, (X, X)) === dimz
    @test dimnum(dimz, (X, X)) === (1, 2)
    @test hasdim(dimz, (X, X)) === (true, true)
    @test sortdims(dimz, (X, X)) === dimz
end

@testset "applying function on a dimension" begin
    d = X(0:0.01:2π)
    a = DimArray(cos, d)
    @test length(dims(a)) == 1
    @test typeof(dims(a)[1]) <: X
    @test a.data == cos.(d.val)
end


text:
"""
Central token type: a flyweight token plus its buffer
"""
struct Token{F<:FlyToken} <: AbstractToken
    fly :: F # category, size, offset, possibly content
    buffer :: String # memory with token text data or EMPTYSTRING
    Base.@propagate_inbounds function Token{F}(f::F, buffer::String) where F<:FlyToken
        @boundscheck isdirect(f) || checkulimit(f.ofs+f.len,usize(buffer))
        new(f,buffer)
    end
end

## GenericToken constructors

# is NOT redundant because here is Token not type-qualified
Base.@propagate_inbounds function Token(f::F, buffer::String) where F<:FlyToken
    Token{F}(f,buffer)
end

"""
Recommended general-purpose Token implementation.

Immutable token, contents either directly stored or buffer based.
Use this type if you expect many token sizes below 8, and want to reduce
memory consumptionand better data locality. It comes at the price of an
additional comparison for most token operations (casting to BToken or DirectFly).

# implementation hint:

Constructors using Strings as content source will encode contents with
up to 7 bytes directly. Only construction from BToken will use buffered content,
even if byte count is below 8. To guarantee content is always buffer based
(which is not the intention behind the Token type), construct using a
BToken.
"""
const HToken = Token{HybridFly}

function Base.show(io::IO,t::Token{T}) where T
    print(io, T <: HybridFly ? 'H' : 'B',Char(t.cat))
    Base.print_quoted(io, t)
end

"""
Immutable token, always buffer based.

Use this type if you expect mostly a token size above 7
"""
const BToken = Token{BufferFly}

HToken(t::BufferFly,s::String) = HToken(hf(t),s)
HToken(t::DirectFly) = @inbounds HToken(hf(t),EMPTYSTRING)
HToken(t::BToken) = @inbounds HToken(t.fly,t.buffer)
HToken(t::HToken) = t

BToken(t::BToken) = t
BToken(t::HToken) = isdirect(t) ? BToken(df(t.fly)) : BToken(bf(t.fly),t.buffer)

Token{F}(cat::Nibble, s::String) where F<:FlyToken = @inbounds Token{F}(cat,0%UInt32,usize(s),s)
Token{F}(cat::Nibble) where F<:FlyToken = @inbounds Token{F}(cat,EMPTYSTRING)
Token{F}(t::AbstractToken) where F<:FlyToken = Token{F}(t.cat,0%UInt32,t.len,substring(t))

Base.@propagate_inbounds function BToken(cat::Nibble,offset::UInt32, size::UInt64, s::String)
    size == 0 && return BToken(cat)
    BToken(BufferFly(cat,size)+offset,s)
end

Base.@propagate_inbounds function HToken(cat::Nibble,offset::UInt32, size::UInt64, s::String)
    size <= MAX_DIRECT_SIZE && return HToken(DirectFly(cat,offset,size,s))
    HToken(BufferFly(cat,size)+offset,s)
end

Token{T}(cat::Nibble,s::SubString{String}) where T =
  Token{T}(cat,s.offset,usize(s),s.string)

Token{T}(cat::Nibble,s::AbstractString)  where T =
  Token{T}(cat,string(s))

Base.@propagate_inbounds function BToken(offset::UInt32, size::UInt64,t::Union{Token,DirectFly})
    @boundscheck checkulimit(offset+size,t.len)
    @inbounds isdirect(t) ? BToken(t.cat,offset,size,string(t)) : BToken(t.cat,offset+offset(t),size,t.buffer)
end

Base.@propagate_inbounds HToken(offset::UInt32, size::UInt64,t::Union{Token,DirectFly}) =
    isdirect(t) ? HToken(DirectFly(offset,size,t)) : HToken(BToken(offset,size,t))

# conversions of standard types
Token{F}(v::Integer) where F<:FlyToken = Token{F}(T_INT,string(v))
Token{F}(v::Real) where F<:FlyToken =    Token{F}(T_REAL,string(v))
Token{F}(v::Char) where F<:FlyToken =    Token{F}(T_CHAR,v)
Token{F}(v::Bool) where F<:FlyToken=    Token{F}(T_KEY,string(v))
Token{F}(v::AbstractString) where F<:FlyToken= Token{F}(T_TEXT,v)
# was soll das? Token{F}(offset::UInt32, size::UInt64,s::String) = where F<:FlyToken Token{F}(T_TEXT,s)

DirectFly(t::Token{T}) where T = isdirect(t) ? df(t.fly) : DirectFly(t.cat,offset(t),t.len,t.buffer)

## special values for missing and nothing

HToken(::Nothing) = HToken(DIRECT_NOTHING)
BToken(::Nothing) = BToken(BufferFly(nothing))
HToken(::Missing) = HToken(DIRECT_MISSING)
BToken(::Missing) = BToken(BufferFly(missing))

## Token API methods ##

offset(t::Token) = offset(t.fly)

category(t::Token) = category(t.fly)

isdirect(t::Token) = isdirect(t.fly)

usize(t::Token) = usize(t.fly)

## ByteString API

"""
ByteString is a family of types which support the following API:

 * pointer(t) -> Ptr{UInt8}

 * usize(t) -> number of bytes as UInt64

 * byte(t,ofs::UInt32) -> UInt8 value at offset ofs of internal buffer

"""
const ByteString = Union{String,SubString{String},BToken,Vector{UInt8},Vector{Int8}}

usize(s::Union{Vector{UInt8},Vector{Int8}}) = length(s)%UInt64

"read codeunit with offset (aka 0-based index)"
Base.@propagate_inbounds function byte(s::ByteString, ofs::UInt32)
    @boundscheck checkbyteofs(ofs,s)
    b = GC.@preserve s unsafe_load(pointer(s)+ofs)
    return b
end

Base.@propagate_inbounds function byte(s::HToken, ofs::UInt32)
    isdirect(s) && return byte(s.fly,ofs)
    @boundscheck checkbyteofs(ofs,s)
    @inbounds byte(s.buffer,ofs+s.ofs)
end

## Base methods for tokens ##

Base.pointer(t::BToken) = pointer(t.buffer)+offset(t)

Base.cmp(a::HToken, b::DirectFly) = isdirect(a.fly) ? cmp(a.fly,b) : cmp_codeunits(a,b)

Base.cmp(a::DirectFly, b::HToken) = isdirect(b.fly) ? cmp(a,b.fly) : cmp_codeunits(a,b)

Base.cmp(a::HToken, b::HToken) = reinterpret(Int64,u64(a.fly)|u64(b.fly))>=0 ? cmp(a.fly,b.fly) : cmp_codeunits(a,b)

## substring implementations

Base.@propagate_inbounds function substring(offset::UInt32, size::UInt64, t::Token{F}) where F
    isdirect(t) && return substring(offset,size,df(t.fly))
    @boundscheck checkulimit(offset+size,t.len) # necessary to ensure substring is in bounds of s, not only of s.string
    @inbounds return substring(offset+t.offset,size,t.buffer)
end

## converts

Base.convert(::Type{Token{F}}, s::AbstractString) where F <:FlyToken = Token{F}(T_TEXT,s)

Base.convert(::Type{Token{F}}, t::AbstractToken) where F <:FlyToken = Token{F}(t)

Base.convert(::Type{SubString{String}}, t::BToken) = @inbounds SubString(offset(t),t.len,t.buffer)

Base.convert(::Type{SubString{String}}, t::HToken) = convert(SubString{String},BToken(t))

Base.convert(::Type{SubString{String}}, t::DirectFly) = convert(SubString{String},string(t))

Base.convert(::Type{String}, t::AbstractToken) = string(t)

function Base.codeunit(t::Token, i::Int)
    if isdirect(t.fly) # optimized away for BToken because isdirect is constantly false
        codeunit(df(t.fly),i)
    else
        @boundscheck checkbounds(t, i)
        @inbounds return codeunit(t.buffer, offset(bf(t.fly)) + i)
    end
end

function Base.hash(t::BToken, h::UInt)
    h += base.memhash_seed
    # note: use pointer(s) here (see #6058).
    s = t.buffer
    GC.@preserve s ccall(memhash, UInt, (Ptr{UInt8}, Csize_t, UInt32), pointer(s)+offset(t.fly), sizeof(t.fly), h % UInt32) + h
end

function Base.hash(t::HToken, h::UInt)
    isdirect(t) && return hash(df(t.fly),h)
    hash(BToken(t),h)
end

#=
function Base.iterate(t::Token, i::Integer=firstindex(t))
    i > ncodeunits(t) && return nothing
    @boundscheck checkulimiti <= 0 && @boundscheck checkbounds(t, i) # ?? wozu??
    o = offset(t)
    if isdirect(t)
        y = iterate(df(t.fly),i)
    else
        y = iterate(t.buffer, o + i)
    end
    c, i = y
    return c, i - o
end
=#

function Base.getindex(t::Token, i::Integer)
    @boundscheck checkbounds(t, i)
    @inbounds if isdirect(t)
        getindex(df(t.fly),i)
    else
        getindex(t.buffer, offset(bf(t.fly)) + i)
    end
end

function Base.read(io::IO, ::Type{HToken})
    cat_size = read(io,Packed31)
    size = bits4_30(cat_size)
    if size <= MAX_DIRECT_SIZE
        t = read(io,cat_size,DirectFly)
        HToken(t)
    else
        @inbounds HToken(Nibble(bits0_3(cat_size)), read(io,size,String))
    end
end

function Base.read(io::IO, ::Type{BToken})
    cat_size = read(io,Packed31)
    size = bits4_30(cat_size)
    @inbounds BToken(BufferFly(cat_size), read(io,size,String))
end

# serialize Token
function Base.write(io::IO, t::Token)
    tt = t.fly
    if isdirect(tt) # optimized away for BufferFly
        write(io,df(tt))
    else
        size = usize(bf(tt))
        p = Packed31(category(tt)%UInt8, size%UInt32)
        write(io,p)
        twrite(io,offset(bf(tt)),size,t.buffer)
    end
end


text:
using RNSpectra
using Test
using DataFrames,DelimitedFiles,Glob

@testset "RNSpectra.jl" begin
    f=Glob.glob("*.txt")
    dx=readdlm(f[1],',',skipstart=2)
    df=DataFrame(dx,:auto)
    Abn=df[:,2]
    λ = df[:,1]
    En = 1240*λ.^-1
    t=798e-7 # nm
    α = log(10)*Abn/t
    RNSpectra.toyozawa(α,En)

end