text:
using Random: shuffle, seed!
using SymbolicUtils: getdepth, Rewriters

@testset "Chain, Postwalk and Fixpoint" begin
    @syms w z α::Real β::Real

    r1 = @rule ~x + ~x => 2 * (~x)
    r2 = @rule ~x * +(~~ys) => sum(map(y-> ~x * y, ~~ys));

    rset = Rewriters.Postwalk(Rewriters.Chain([r1, r2]))
    @test getdepth(rset) == typemax(Int)

    ex = 2 * (w+w+α+β)

    @eqtest rset(ex) == (((2 * w) + (2 * w)) + (2 * α)) + (2 * β)
    @eqtest Rewriters.Fixpoint(rset)(ex) == ((2 * (2 * w)) + (2 * α)) + (2 * β)
end

@testset "Numeric" begin
    @syms a::Integer b c d x::Real y::Number
    @eqtest simplify(x - y) == x + -1*y
    @eqtest simplify(x - sin(y)) == x + -1*sin(y)
    @eqtest simplify(-sin(x)) == -1*sin(x)
    @eqtest simplify(1 * x * 2) == 2 * x
    @eqtest simplify(1 + x + 2) == 3 + x
    @eqtest simplify(b*b) == b^2 # tests merge_repeats
    @eqtest simplify((a*b)^2) == a^2 * b^2
    @eqtest simplify((a*b)^c) == (a*b)^c

    @eqtest simplify(1x + 2x) == 3x
    @eqtest simplify(3x + 2x) == 5x

    @eqtest simplify(a + b + (x * y) + c + 2 * (x * y) + d)     == (3 * x * y) + a + b + c + d
    @eqtest simplify(a + b + 2 * (x * y) + c + 2 * (x * y) + d) == (4 * x * y) + a + b + c + d

    @eqtest simplify(a * x^y * b * x^d) == (a * b * (x ^ (d + y)))

    @eqtest simplify(a + b + 0*c + d) == a + b + d
    @eqtest simplify(a * b * c^0 * d) == a * b * d
    @eqtest simplify(a * b * 1*c * d) == a * b * c * d

    @test simplify(Term(one, [a])) == 1
    @test simplify(Term(one, [b+1])) == 1
    @test simplify(Term(one, [x+2])) == 1

    @test simplify(Term(zero, [a])) == 0
    @test simplify(Term(zero, [b+1])) == 0
    @test simplify(Term(zero, [x+2])) == 0
end

@testset "boolean" begin
    @syms a::Real b c

    @eqtest simplify(a < 0) == (a < 0)
    @eqtest simplify(0 < a) == (0 < a)
    @eqtest simplify((0 < a) | true) == true
    @eqtest simplify(true | (0 < a)) == true
    @eqtest simplify((0 < a) & true) == (0 < a)
    @eqtest simplify(true & (0 < a)) == (0 < a)
    @eqtest simplify(false & (0 < a)) == false
    @eqtest simplify((0 < a) & false) == false
    @eqtest simplify(Term{Bool}(!, [true])) == false
    @eqtest simplify(Term{Bool}(|, [false, true])) == true
    @eqtest simplify(cond(true, a,b)) == a
    @eqtest simplify(cond(false, a,b)) == b

    # abs
    @test simplify(substitute(cond(!(a < 0), a,-a), Dict(a=>-1))) == 1
    @test simplify(substitute(cond(!(a < 0), a,-a), Dict(a=>1))) == 1
    @test simplify(substitute(cond(a < 0, -a, a), Dict(a=>-1))) == 1
    @test simplify(substitute(cond(a < 0, -a, a), Dict(a=>1))) == 1
end

@testset "Pythagorean Identities" begin
    @syms a::Integer x::Real y::Number

    @test simplify(cos(x)^2 + 1 + sin(x)^2) == 2
    @test simplify(cos(y)^2 + 1 + sin(y)^2) == 2
    @test simplify(sin(y)^2 + cos(y)^2 + 1) == 2

    @eqtest simplify(1 + y + tan(x)^2) == sec(x)^2 + y
    @eqtest simplify(1 + y + cot(x)^2) == csc(x)^2 + y
end

@testset "Depth" begin
    @syms x
    R = Rewriters.Postwalk(Rewriters.Chain([@rule(sin(~x) => cos(~x)),
                                            @rule( ~x + 1 => ~x - 1)]))
    @eqtest R(sin(sin(sin(x + 1)))) == cos(cos(cos(x - 1)))
    #@eqtest R(sin(sin(sin(x + 1))), depth=2) == cos(cos(sin(x + 1)))
end

@testset "RuleRewriteError" begin
    pred(x) = error("Fail")

    @syms a b

    rs = Rewriters.Postwalk(Rewriters.Chain(([@rule ~x + ~y::pred => ~x])))
    @test_throws SymbolicUtils.RuleRewriteError rs(a+b)
    err = try rs(a+b) catch err; err; end
    @test sprint(io->Base.showerror(io, err)) == "Failed to apply rule ~x + ~(y::pred) => ~x on expression a + b"
end

@testset "Threading" begin
    @syms a b c d
    ex = (((0.6666666666666666 / (c / 1)) + ((1 * a) / (c / 1))) +
          (1.0 / (((1 * d) / (1 + b)) * (1 / b)))) +
          ((((1 * a) + (1 * a)) / ((2.0 * (d + 1)) / 1.0)) +
           ((((d * 1) / (1 + c)) * 2.0) / ((1 / d) + (1 / c))))
    @eqtest simplify(ex) == simplify(ex, threaded=true, thread_subtree_cutoff=3)
    @test SymbolicUtils.node_count(a + b * c / d) == 7
end

@testset "timerwrite" begin
    @syms a b c d
    expr1 = foldr((x,y)->rand([*, /])(x,y), rand([a,b,c,d], 100))
    SymbolicUtils.@timerewrite simplify(expr1)
end


text:
const QNBlock = Pair{QN,Int64}

const QNBlocks = Vector{QNBlock}

qn(qnblock::QNBlock) = qnblock.first

NDTensors.blockdim(qnblock::QNBlock) = qnblock.second

NDTensors.blockdim(qnblocks::QNBlocks,
                   b::Int) = blockdim(qnblocks[b])

qn(qnblocks::QNBlocks,b::Int) = qn(qnblocks[b])

NDTensors.nblocks(qnblocks::QNBlocks) = length(qnblocks)

function NDTensors.dim(qnblocks::QNBlocks)
  dimtot = 0
  for (_,blockdim) in qnblocks
    dimtot += blockdim
  end
  return dimtot
end

function Base.:-(qnb::QNBlock)
  return QNBlock(-qn(qnb),blockdim(qnb))
end

function Base.:+(qn1::QNBlock,qn2::QNBlock)
  qn(qn1) != qn(qn2) && error("Cannot add qn blocks with different qns")
  return QNBlock(qn(qn1),blockdim(qn1)+blockdim(qn2))
end

function Base.:-(qns::QNBlocks)
  qns_new = copy(qns)
  for i in 1:length(qns_new)
    qns_new[i] = -qns_new[i]
  end
  return qns_new
end

"""
A QN Index is an Index with QN block storage instead of
just an integer dimension. The QN block storage is a
vector of pairs of QNs and block dimensions.
The total dimension of a QN Index is the sum of the
dimensions of the blocks of the Index.
"""
const QNIndex = Index{QNBlocks}

const QNIndexVal = IndexVal{QNIndex}

hasqns(::QNIndex) = true

function have_same_qns(qnblocks::QNBlocks)
  qn1 = qn(qnblocks,1)
  for n in 2:nblocks(qnblocks)
    !have_same_qns(qn1,qn(qnblocks,n)) && return false
  end
  return true
end

function have_same_mods(qnblocks::QNBlocks)
  qn1 = qn(qnblocks,1)
  for n in 2:nblocks(qnblocks)
    !have_same_mods(qn1,qn(qnblocks,n)) && return false
  end
  return true
end

"""
    Index(qnblocks::Vector{Pair{QN, Int64}}; dir::Arrow = Out,
                                             tags = "",
                                             plev::Int = 0)

Construct a QN Index from a Vector of pairs of QN and block
dimensions.

Note: in the future, this may enforce that all blocks have the
same QNs (which would allow for some optimizations, for example
when constructing random QN ITensors).

# Example
```
Index([QN("Sz", -1) => 1, QN("Sz", 1) => 1]; tags = "i")
```
"""
function Index(qnblocks::QNBlocks; dir::Arrow = Out, tags = "", plev = 0)
  # TODO: make this a debug check?
  #have_same_qns(qnblocks) || error("When creating a QN Index, the QN blocks must have the same QNs")
  #have_same_mods(qnblocks) || error("When creating a QN Index, the QN blocks must have the same mods")
  return Index(rand(IDType), qnblocks, dir, tags, plev)
end

"""
    Index(qnblocks::Vector{Pair{QN, Int64}}, tags; dir::Arrow = Out,
                                                   plev::Int = 0)

Construct a QN Index from a Vector of pairs of QN and block
dimensions.

# Example
```
Index([QN("Sz", -1) => 1, QN("Sz", 1) => 1], "i"; dir = In)
```
"""
Index(qnblocks::QNBlocks,
      tags;
      dir::Arrow = Out,
      plev::Int = 0) = Index(qnblocks; dir = dir,
                                       tags = tags,
                                       plev = plev)

"""
    Index(qnblocks::Pair{QN, Int64}...; dir::Arrow = Out,
                                        tags = "",
                                        plev::Int = 0)

Construct a QN Index from a list of pairs of QN and block
dimensions.

# Example
```
Index(QN("Sz", -1) => 1, QN("Sz", 1) => 1; tags = "i")
```
"""
function Index(qnblocks::QNBlock...; dir::Arrow=Out,
                                     tags="",
                                     plev=0)
  return Index([qnblocks...]; dir = dir,
                              tags = tags,
                              plev = plev)
end

NDTensors.dim(i::QNIndex) = dim(space(i))

NDTensors.nblocks(i::QNIndex) = nblocks(space(i))

qn(ind::QNIndex,b::Int) = dir(ind)*qn(space(ind),b)

qnblocks(ind::QNIndex) = space(ind)

NDTensors.blockdim(ind::QNIndex,
                   b::Int) = blockdim(space(ind),b)

function qn(iv::QNIndexVal)
  i = ind(iv)
  v = val(iv)
  tdim = 0
  for b=1:nblocks(i)
    tdim += blockdim(i,b)
    (v <= tdim) && return qn(i,b)
  end
  error("qn: QNIndexVal out of range")
  return QN()
end

"""
    qnblocknum(ind::QNIndex,q::QN)

Given a QNIndex `ind` and QN `q`, return the
number of the block (from 1,...,nblocks(ind))
of the QNIndex having QN equal to `q`. Assumes
all blocks of `ind` have a unique QN.
"""
function qnblocknum(ind::QNIndex,q::QN)
  for b=1:nblocks(ind)
    if qn(ind,b) == q
      return b
    end
  end
  error("No block found with QN equal to $q")
  return 0
end

"""
    qnblockdim(ind::QNIndex, q::QN)

Given a QNIndex `ind` and QN `q`, return the
dimension of the block of the QNIndex having
QN equal to `q`. Assumes all blocks of `ind`
have a unique QN.
"""
qnblockdim(ind::QNIndex, q::QN) = blockdim(ind, qnblocknum(ind,q))

function Base.:*(dir::Arrow, qnb::QNBlock)
  return QNBlock(dir*qn(qnb), blockdim(qnb))
end

function Base.:*(dir::Arrow, qn::QNBlocks)
  qnR = copy(qn)
  for i in 1:nblocks(qnR)
    qnR[i] = dir*qnR[i]
  end
  return qnR
end

function Base.:*(qn1::QNBlock, qn2::QNBlock)
  return QNBlock(qn(qn1)+qn(qn2), blockdim(qn1)*blockdim(qn2))
end

function NDTensors.outer(qn1::QNBlocks, qn2::QNBlocks)
  qnR = ITensors.QNBlocks(undef,nblocks(qn1)*nblocks(qn2))
  for (i,t) in enumerate(Iterators.product(qn1,qn2))
    qnR[i] = prod(t)
  end
  return qnR
end

function NDTensors.outer(i1::QNIndex, i2::QNIndex;
                         dir = nothing,
                         tags = "",
                         plev::Int = 0)
  if isnothing(dir)
    if ITensors.dir(i1) == ITensors.dir(i2)
      dir = ITensors.dir(i1)
    else
      dir = Out
    end
  end
  newspace = dir * ((ITensors.dir(i1) * space(i1)) ⊗
                    (ITensors.dir(i2) * space(i2)))
  return Index(newspace;
               dir = dir,
               tags = tags,
               plev = plev)
end

function NDTensors.outer(i::QNIndex;
                         dir = nothing,
                         tags = "",
                         plev::Int = 0)
  if isnothing(dir)
    dir = ITensors.dir(i)
  end
  newspace = dir * (ITensors.dir(i) * space(i))
  return Index(newspace;
               dir = dir,
               tags = tags,
               plev = plev)
end

function Base.isless(qnb1::QNBlock, qnb2::QNBlock)
  return isless(qn(qnb1), qn(qnb2))
end

function NDTensors.permuteblocks(i::QNIndex, perm)
  qnblocks_perm = space(i)[perm]
  return replaceqns(i, qnblocks_perm)
end

function combineblocks(qns::QNBlocks)
  perm = sortperm(qns)
  qnsP = qns[perm]
  qnsC = [qnsP[1]]
  comb = Vector{Int}(undef,nblocks(qns))

  # Which block this is, after combining
  block_count = 1
  comb[1] = block_count
  for i in 2:nblocks(qnsP)
    if qn(qnsP[i]) == qn(qnsP[i-1])
      qnsC[block_count] += qnsP[i]
    else
      push!(qnsC,qnsP[i])
      block_count += 1
    end
    comb[i] = block_count
  end
  return qnsC,perm,comb
end

# Make a new Index with the specified qn blocks
function replaceqns(i::QNIndex,qns::QNBlocks)
  return Index(id(i),qns,dir(i),tags(i),plev(i))
end

function NDTensors.setblockdim!(i::QNIndex,newdim::Int,n::Int)
  qns = space(i)
  qns[n] = qn(qns[n]) => newdim
  return i
end

function setblockqn!(i::QNIndex,newqn::QN,n::Int)
  qns = space(i)
  qns[n] = newqn => blockdim(qns[n])
  return i
end

function Base.deleteat!(i::QNIndex,pos)
  deleteat!(space(i),pos)
  return i
end

function Base.resize!(i::QNIndex,n::Integer)
  resize!(space(i),n)
  return i
end

function combineblocks(i::QNIndex)
  qnsR,perm,comb = combineblocks(space(i))
  iR = replaceqns(i,qnsR)
  return iR,perm,comb
end

removeqns(i::QNIndex) = Index(id(i),
                              dim(i),
                              Neither,
                              tags(i),
                              plev(i))

function Base.show(io::IO,
                   i::QNIndex)
  idstr = "$(id(i) % 1000)"
  if length(tags(i)) > 0
    print(io,"(dim=$(dim(i))|id=$(idstr)|\"$(tagstring(tags(i)))\")$(primestring(plev(i)))")
  else
    print(io,"(dim=$(dim(i))|id=$(idstr))$(primestring(plev(i)))")
  end
  println(io," <$(dir(i))>")
  for (n,qnblock) in enumerate(space(i))
    print(io," $n: $qnblock")
    n < length(space(i)) && println(io)
  end
end

function HDF5.write(parent::Union{HDF5File, HDF5Group},
                    name::AbstractString,
                    B::QNBlocks)
  g = g_create(parent, name)
  attrs(g)["type"] = "QNBlocks"
  attrs(g)["version"] = 1
  write(g,"length",length(B))
  dims = [block[2] for block in B]
  write(g,"dims",dims)
  for n=1:length(B)
    write(g,"QN[$n]",B[n][1])
  end
end

function HDF5.read(parent::Union{HDF5File,HDF5Group},
                   name::AbstractString,
                   ::Type{QNBlocks})
  g = g_open(parent,name)
  if read(attrs(g)["type"]) != "QNBlocks"
    error("HDF5 group or file does not contain QNBlocks data")
  end
  N = read(g,"length")
  dims = read(g,"dims")
  B = QNBlocks(undef,N)
  for n=1:length(B)
    B[n] = QNBlock(read(g,"QN[$n]",QN),dims[n])
  end
  return B
end

ext: jl


text:
using Optim, LineSearches
stress_scale = 1.0e5
strain_scale = 1.0
T = 200.0
NT = 200

include("nnutil.jl")

# H0 = constant(H1/stress_scale)
testtype = "NeuralNetwork2D"

force_scales = [5.0]

force_scale = 5.0

H0 = [1335174.0968380707  326448.3267263398 0.0
       326448.3267263398 1326879.2022994285 0.0
            0.0                0.0     526955.763626241]/stress_scale

H0inv = inv(H0)

n_data = [100, 101, 102, 103, 104,  200, 201, 202, 203, 204]

porder = 2
# density 4.5*(1 - 0.25) + 3.2*0.25
fiber_fraction = 0.25
#todo
#fiber_fraction = 1.0
prop = Dict("name"=> testtype, "rho"=> 4.5*(1 - fiber_fraction) + 3.2*fiber_fraction, "nn"=>nn)

# DNS computaional domain
fiber_size = 5
# nx_f, ny_f = 40*fiber_size, 20*fiber_size
nx_f, ny_f = 80*fiber_size, 40*fiber_size

# nx_f, ny_f = 12, 4

# homogenized computaional domain
# number of elements in each directions
nx, ny = 10, 5
# number of subelements in one element in each directions
sx_f, sy_f = div(nx_f,nx), div(ny_f,ny)

ndofs = 2
fine_to_coarse = zeros(Int64, ndofs*(nx*porder+1)*(ny*porder+1))
for idof = 1:ndofs
    for iy = 1:ny*porder+1
        for ix = 1:nx*porder+1
            fine_to_coarse[ix + (iy - 1)*(nx*porder+1) + (idof-1)*(nx*porder+1)*(ny*porder+1)] =
            1 + (ix - 1) * sx_f + (iy - 1) * (nx_f*porder + 1) * sy_f + (nx_f*porder + 1)*(ny_f*porder + 1)*(idof - 1)
        end
    end
end

Δt = T/NT

nodes, _, _, _, _, _, _ = BoundaryCondition(n_data[1], nx, ny, porder)
elements = []
for j = 1:ny
    for i = 1:nx
        n = (nx*porder+1)*(j-1)*porder + (i-1)porder+1
        #element (i,j)
        if porder == 1
            #   4 ---- 3
            #
            #   1 ---- 2

            elnodes = [n, n + 1, n + 1 + (nx + 1), n + (nx + 1)]
        elseif porder == 2
            #   4 --7-- 3
            #   8   9   6
            #   1 --5-- 2
            elnodes = [n, n + 2, n + 2 + 2*(2*nx+1),  n + 2*(2*nx+1), n+1, n + 2 + (2*nx+1), n + 1 + 2*(2*nx+1), n + (2*nx+1), n+1+(2*nx+1)]
        else
            error("polynomial order error, porder= ", porder)
        end
        coords = nodes[elnodes,:]
        push!(elements,SmallStrainContinuum(coords,elnodes, prop, 3))
    end
end

#####################Step 1 preprocess data
function PreprocessData(n_data, NT)
    globdat_arr = []
    domain_arr = []
    obs_state_arr = []

    for (id, tid) in enumerate(n_data)

        nodes, EBC, g, gt, FBC, fext, ft = BoundaryCondition(tid, nx, ny, porder)
        domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext)
        neqs = domain.neqs
        globdat = GlobalData(zeros(domain.neqs),zeros(domain.neqs), zeros(domain.neqs),zeros(domain.neqs), neqs, gt, ft)
        assembleMassMatrix!(globdat, domain)

        obs_state = zeros(Float64, NT+1, domain.neqs)
        full_state_history, _ = read_data("Data/order$porder/$(tid)_$(force_scale)_$(fiber_size).dat")
        @assert length(full_state_history) == NT+1
        for i = 1:NT+1
            obs_state[i,:] = (full_state_history[i][fine_to_coarse])[domain.dof_to_eq]
        end

        push!(domain_arr, domain)
        push!(globdat_arr, globdat)
        push!(obs_state_arr, obs_state)
    end

    return globdat_arr, domain_arr, obs_state_arr
end

globdat_arr, domain_arr, obs_state_arr = PreprocessData(n_data, NT)

mutable struct Buffer
    J::Array{Float64}
    dJ::Array{Array{Float64}}
    dts::Array{Any,1}
    state::Array{Any,1}
    strain::Array{Any,1}
    stress::Array{Any,1}

    function Buffer(n::Int64, ntheta::Int64, NT::Int64, neqs::Array{Int64}, ngps::Int64, nstrain::Int64)
        J = zeros(n)
        dJ = [zeros(ntheta) for i = 1:n]
        dts = [[] for i = 1:n]
        state = [[] for i = 1:n]
        strain = [[] for i = 1:n]
        stress = [[] for i = 1:n]
        new(J, dJ, dts, state, strain, stress)
    end
end

function calculate_common!(theta, last_theta, buffer)
    @show " theta norm ", norm(theta), " last_theta norm ", norm(last_theta)
    if theta != last_theta

        copy!(last_theta, theta)

        Threads.@threads for i = 1:length(n_data)
        #for i = 1:length(n_data)
            #@show Threads.threadid()
            # todo: inplace
            @show n_data[i]
            buffer.J[i], buffer.dts[i], buffer.state[i], buffer.strain[i], buffer.stress[i] =
            ForwardNewmarkSolver(globdat_arr[i], domain_arr[i], theta, T, NT, strain_scale, stress_scale, obs_state_arr[i])

            #if buffer.J[i]  == Inf
            #    break
            #end
        end
    end
end

function f(theta, buffer, last_theta)
    calculate_common!(theta, last_theta, buffer)

    @show buffer.J
    J = sum(buffer.J)

    @show "function evaluation: |J|=", norm(J)

    return J
end

function g!(theta, storage, buffer, last_theta)
    calculate_common!(theta, last_theta, buffer)

    Threads.@threads for i = 1:length(n_data)
        # todo: inplace
        buffer.dJ[i] = BackwardNewmarkSolver(globdat_arr[i], domain_arr[i], theta, T, NT, buffer.dts[i], buffer.state[i], buffer.strain[i], buffer.stress[i], strain_scale, stress_scale, obs_state_arr[i])
    end

    storage[:] = sum(buffer.dJ)

    @show "gradient evaluation: |dJ/dtheta|=", norm( storage )

    storage[:] ./= max(1.0, norm(storage))
end

nstrain = 3
ngps_per_elem = length(domain_arr[1].elements[1].weights)
neles = domain_arr[1].neles

config = [9, 20, 20, 20, 20, 4]
start_id=5
initial_theta = convert_mat("nn2array", config,  "Data/piecewise/NNPreLSfit_$(idx)_$(H_function)_$(start_id).mat")

#initial_theta = 1.0e-3*rand((9+1)*20 + (20+1)*20 + (20+1)*20 + (20+1)*4)

#@load "Data/theta_ite_0_5.jld2" last_theta
#initial_theta = copy(last_theta)

neqs_arr = [domain_arr[i].neqs for i = 1:length(n_data)]
buffer = Buffer(length(n_data), length(initial_theta), NT, neqs_arr, neles*ngps_per_elem, nstrain) # Preallocate an appropriate buffer
last_theta = similar(initial_theta)
#algo = LBFGS(alphaguess = InitialStatic(0.01), linesearch=LineSearches.BackTracking(order=3))
algo = LBFGS(alphaguess = InitialPrevious(), linesearch=LineSearches.BackTracking(order=3))
#algo = GradientDescent(alphaguess = InitialStatic(), linesearch=LineSearches.BackTracking(order=2))

for i = 1:50
    println("************************** Outer Iteration = $i ************************** ")

    optimize(x -> f(x, buffer, last_theta),
        (stor, x) -> g!(x, stor, buffer, last_theta),
        initial_theta,
        algo,
        Optim.Options(iterations=100))

        initial_theta[:] = last_theta
        @save "Data/piecewise/theta_ite_idx$(idx)_from$(start_id)_$(i).jld2" last_theta
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Profile

function spawnmany(n)
    if n > 2
        m = n ÷ 2
        t = Threads.@spawn spawnmany(m)
        spawnmany(m)
        wait(t)
    end
end

@profile spawnmany(parse(Int, get(ENV, "NTASKS", "2000000")))

max_stars_repo_name: JuliaLabs/julia

max_issues_repo_name: JuliaLabs/julia

max_forks_repo_name: JuliaLabs/julia


text:
using QuasinormalModes
using SymEngine

# ------------------------------------------------------------------
# 1. Analytic Schwarzschild Black Hole
# ------------------------------------------------------------------

struct SchwarzschildData{N,T} <: QuadraticEigenvalueProblem{N,T}
    nIter::N
    x0::T

    vars::Tuple{Basic, Basic}
    exprs::Tuple{Basic, Basic}
end

function SchwarzschildData(nIter::N, x0::T, l::N, s::N) where {N,T}
    vars = @vars x ω

    λ0 = (-1 + (2*im)*ω + x*(4 - 3*x + (4*im)*(-2 + x)*ω))/((-1 + x)^2*x)
    S0 = (l + l^2 + (-1 + s^2)*(-1 + x) + (4*im)*(-1 + x)*ω + 4*(-2 + x)*ω^2)/((-1 + x)^2*x)

    return SchwarzschildData{N,T}(nIter, x0, vars, (λ0, S0))
end

QuasinormalModes.λ0(d::SchwarzschildData{N,T}) where {N,T} = d.exprs[1]
QuasinormalModes.S0(d::SchwarzschildData{N,T}) where {N,T}  = d.exprs[2]

QuasinormalModes.get_niter(d::SchwarzschildData{N,T}) where {N,T} = d.nIter
QuasinormalModes.get_x0(d::SchwarzschildData{N,T}) where {N,T} = d.x0

QuasinormalModes.get_ODEvar(d::SchwarzschildData{N,T}) where {N,T} = d.vars[1]
QuasinormalModes.get_ODEeigen(d::SchwarzschildData{N,T}) where {N,T} = d.vars[2]

# ------------------------------------------------------------------
# 2. Numeric Schwarzschild Black Hole
# ------------------------------------------------------------------

struct NSchwarzschildData{N,T} <: NumericAIMProblem{N,T}
    nIter::N
    x0::T
    l::N
    s::N
end

function NSchwarzschildData(nIter::N, x0::T, l::N, s::N) where {N,T}
    return NSchwarzschildData{N,T}(nIter, x0, l, s)
end

QuasinormalModes.λ0(::NSchwarzschildData{N,T}) where {N,T} = (x,ω) -> (-1 + (2*im)*ω + x*(4 - 3*x + (4*im)*(-2 + x)*ω))/((-1 + x)^2*x)
QuasinormalModes.S0(d::NSchwarzschildData{N,T}) where {N,T} = (x,ω) -> (d.l + d.l^2 + (-1 + d.s^2)*(-1 + x) + (4*im)*(-1 + x)*ω + 4*(-2 + x)*ω^2)/((-1 + x)^2*x)

QuasinormalModes.get_niter(d::NSchwarzschildData{N,T}) where {N,T} = d.nIter
QuasinormalModes.get_x0(d::NSchwarzschildData{N,T}) where {N,T} = d.x0

# ------------------------------------------------------------------
# 3. Constructing problems and caches
# ------------------------------------------------------------------

p_ana = SchwarzschildData(0x00030, Complex(BigFloat("0.43"), BigFloat("0.0")), 0x00000, 0x00000);
p_num = NSchwarzschildData(0x00030, Complex(0.43, 0.0), 0x00000, 0x00000);

c_ana = AIMCache(p_ana)
c_num = AIMCache(p_num)

# ------------------------------------------------------------------
# 4. Computing quasinormal modes
# ------------------------------------------------------------------

m_ana = computeEigenvalues(Serial(), p_ana, c_ana)

function printQNMs(qnms, cutoff, instab)
    println("-"^165)
    println("|", " "^36, "Re(omega)", " "^36, " ", " "^36, "Im(omega)", " "^36, "|")
    println("-"^165)

    for qnm in qnms
        if real(qnm) > cutoff && ( instab ? true : imag(qnm) < big"0.0" )
        println(real(qnm), "    ", imag(qnm))
        end
    end

    println("-"^165)

    return nothing
end

sort!(m_ana, by = x -> imag(x))
printQNMs(m_ana, 1.0e-10, false)

ev = computeEigenvalues(Serial(), p_num, c_num, Complex(0.22, -0.20), nls_xtol = 1.0e-10, nls_ftol = 1.0e-10)


text:
using Koala
using KoalaEnsembles
import UnicodePlots
using Base.Test

const X, y = load_ames();
const train, test = split(1:length(y), 0.8);
const forest = EnsembleRegressor()
showall(forest)
tree = forest.atom
tree.max_features = 4
forest.n = 50
forestM = Machine(forest, X, y, train)
showall(forestM)
fit!(forestM, train)
forest.weight_regularization = 0.5
fit_weights!(forestM)
err(forestM, test)
fit!(forestM, train, add=true)
err(forestM, test)
u,v = weight_regularization_curve(forestM, test, raw=false, range=linspace(0,1,21))
UnicodePlots.lineplot(u,v)
forest.weight_regularization = u[indmin(v)]
fit_weights!(forestM)
err(forestM, test)
forest.n = 200
fit!(forestM, train)
err(forestM, test)


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

name = "LinearElasticity"
version = v"5.0"

# Collection of sources required to build MMG
sources = [
    GitSource("https://github.com/ISCDtoolbox/LinearElasticity.git",
              "58ace131fd1af9293ded9f84f2a6447147a8ecdc"), # master @ May 6, 2021
    GitSource("https://github.com/ISCDtoolbox/Commons.git",
              "57f0b0ed46d45a12fa60c42bce2c6fa73c84bc30"), # master @ Jul 23, 2021
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/Commons
if [[ "${target}" == *mingw* ]]; then
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/Commons.mingw.patch"
elif [[ "${target}" == *freebsd* ]]; then
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/Commons.freebsd.patch"
else
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/Commons.patch"
fi
mkdir build
cd build
cmake .. \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DBUILD_SHARED_LIBS=ON
make -j${nproc}
make install

cd ${WORKSPACE}/srcdir/LinearElasticity
if [[ "${target}" == *mingw* ]]; then
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/LinearElasticity.mingw.patch"
else
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/LinearElasticity.patch"
fi
mkdir build
cd build
cmake .. \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DCMAKE_C_STANDARD_LIBRARIES="${LDFLAGS}" \
    -DBUILD_SHARED_LIBS=ON
make -j${nproc}
make install
install_license ../LICENSE.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libCommons", :libCommons),
    LibraryProduct("libElas", :libElas),
    ExecutableProduct("elastic", :elastic)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[]

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


text:
module Tutorial4
using FEniCS

T = 2.0            # final time
num_steps = 50     # number of time steps
dt = T / num_steps # time step size

# Create mesh and define function space
nx = ny = 30
mesh = RectangleMesh(Point([-2.0, -2.]), Point([2., 2.]), nx, ny)
V = FunctionSpace(mesh, "P", 1)

bc = DirichletBC(V, Constant(0), "on_boundary")

# Define initial value
u_0 = Expression("exp(-a*pow(x[0], 2) - a*pow(x[1], 2))",
                 degree=2, a=5)

u_n = interpolate(u_0, V)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)

F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)

# Create VTK file for saving solution
#vtkfile = File("heat_gaussian/solution.pvd")

# Time-stepping
u=FeFunction(V)
global t = 0
#fig = PyPlot.figure()
for n = 0:(num_steps-1)
    # Update current time
    global t += dt

    # Compute solution
    lvsolve(a,L,u,bc)

    # Save to file and plot solution
    # vtkfile << (u.pyobject, t)

    # Update previous solution
    assign(u_n,u)
end
end#module


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

#test-warlpiri-division.jl

###############################################################################
##COMPLETE DIVISION TABLE

#column key
#none_ pfew_ pone_ psome ptwo_ pmany pall_ junk_ nnone n_few n_one nsome n_two nmany n_all
wdiv = Utype[
#none_ / line
 junk_ none_ none_ none_ none_ none_ none_ junk_ junk_ none_ none_ none_ none_ none_ none_;
#pfew_ / line
 junk_ pf_pm pfew_ pfew_ pfew_ pfew_ none_ junk_ junk_ nm_nf n_few n_few n_few n_few none_;
#pone_ / line
 junk_ ps_pm pone_ pfew_ pfew_ pfew_ none_ junk_ junk_ nm_ns n_one n_few n_few n_few none_;
#psome / line
 junk_ ps_pm psome pf_ps pfew_ pfew_ none_ junk_ junk_ nm_ns nsome ns_nf n_few n_few none_;
#ptwo_ / line
 junk_ pmany ptwo_ psome pone_ pfew_ none_ junk_ junk_ nmany n_two nsome n_one n_few none_;
#pmany / line
 junk_ pmany pmany ps_pm ps_pm pf_pm none_ junk_ junk_ nmany nmany nm_ns nm_ns nm_nf none_;
#pall_ / line
 junk_ pall_ pall_ pall_ pall_ pall_ junk_ junk_ junk_ n_all n_all n_all n_all n_all junk_;
#junk_ / line
 junk_ junk_ junk_ junk_ junk_ junk_ junk_ junk_ junk_ junk_ junk_ junk_ junk_ junk_ junk_;
#nnone / line
 junk_ none_ none_ none_ none_ none_ none_ junk_ junk_ none_ none_ none_ none_ none_ none_;
#n_few / line
 junk_ nm_nf n_few n_few n_few n_few none_ junk_ junk_ pf_pm pfew_ pfew_ pfew_ pfew_ none_;
#n_one / line
 junk_ nm_ns n_one n_few n_few n_few none_ junk_ junk_ ps_pm pone_ pfew_ pfew_ pfew_ none_;
#nsome / line
 junk_ nm_ns nsome ns_nf n_few n_few none_ junk_ junk_ ps_pm psome pf_ps pfew_ pfew_ none_;
#n_two / line
 junk_ nmany n_two nsome n_one n_few none_ junk_ junk_ pmany ptwo_ psome pone_ pfew_ none_;
#nmany / line
 junk_ nmany nmany nm_ns nm_ns nm_nf none_ junk_ junk_ pmany pmany ps_pm ps_pm pf_pm none_;
#n_all / line
 junk_ n_all n_all n_all n_all n_all junk_ junk_ junk_ pall_ pall_ pall_ pall_ pall_ junk_
]
testop(/, wdiv)

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
@testset "931.minimum-falling-path-sum.jl" begin
    @test min_falling_path_sum([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) == 12
    @test min_falling_path_sum([
        [100, -42, -46, -41], [31, 97, 10, -10], [-58, -51, 82, 89], [51, 81, 69, -51]
    ]) == -36
end


text:
#=##############################################################################
# DESCRIPTION
    Methods for generation and evaluation of grids.
# AUTHORSHIP
  * Author    : Eduardo J Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Mar 2018
  * License   : MIT License
=###############################################################################

################################################################################
# GRID
################################################################################
"""
  `Grid(P_min, P_max, NDIVS)`

Generates an n-dimensional grid.

  **Arguments**
  * `P_min::Array{Float64,1}`   : Minimum point of the domain.
  * `P_max::Array{Float64,1}`   : Maximum point of the domain.
  * `NDIVS::Array{Int64,1}`     : Number of divisions in each coordinate.

  **Properties**
  * `dims::Int64`               : Number of dimensions.
  * `nnodes::Int64`             : Number of nodes in the grid.
  * `nodes::Array{Float64,2}`   : Matrix size (`nnodes`, `dims`) of node position.
  * `field` : Contains calculated fields formated as field[field_name] = Dict(
                            "field_name" => field_name::String,
                            "field_type" => "scalar" or "vector",
                            "entry_type" => "node" or "cell",
                            "field_data" => data
                            )
            where `data` is an array data[i] = [val1, val2, ...] containing
            this field values (scalar or vector) at each node in the grid.

NOTE: All indexing is done linearly, meaning that `nodes` is indexed from 1 to
      `nnodes`, and all data fields follow the same indexing.

NOTE2: `NDIVS` can either be an array of integers with NDIVS[i] indicating the
      number of divisions in the i-th coordinate, or it can be an array of
      sections (see `multidiscretize()` doc) with NDIVS[i] = [sec1, sec2, ...]
      indicating the discretization into sections in the i-th coordinate.
"""
mutable struct Grid <: AbstractGrid

  # User inputs
  P_min::Array{T,1} where {T<:Real}   # Minimum point of the domain
  P_max::Array{T,1} where {T<:Real}   # Maximum point of the domain
  NDIVS::Array{T,1} where {T<:Any}    # Number of divisions in each coordinate
  # Optional inputs
  loop_dim::Int64                     # Index of dimension to close in loop

  # Properties
  dims::Int64                         # Number of dimensions
  nnodes::Int64                       # Number of nodes
  ncells::Int64                       # Number of cells
  nodes::Array{T,2} where{T<:Real}    # Position of each node
  # bbox::Array{Int64, 1}               # Bounding box (cells in each dimension)
  field::Dict{String, Dict{String, Any}}  # Calculated fields

  # Internal data
  _ndivsnodes::Tuple                  # Number of nodes in each coordinate
  _ndivscells::Tuple                  # Number of cells in each coordinate
  _override_vtkcelltype::Int64        # Option for overriding vtk outputs

  Grid(P_min, P_max, NDIVS, loop_dim=0,
              dims=_calc_dims(P_min),
                nnodes=_calc_nnodes(NDIVS, loop_dim),
                ncells=_calc_ncells(NDIVS),
                nodes=_generate_grid(P_min, P_max, NDIVS, loop_dim),
                # bbox=_calc_ndivs(NDIVS),
                field=Dict{String, Dict{String, Any}}(),
              _ndivsnodes=Tuple(_calc_ndivsnodes(NDIVS, loop_dim)),
                _ndivscells=Tuple(_calc_ndivs(NDIVS)),
                _override_vtkcelltype=-1
      ) = _check(P_min, P_max, NDIVS, loop_dim) ? new(P_min, P_max, NDIVS,
              loop_dim,
              dims,
                nnodes,
                ncells,
                nodes,
                # bbox,
                field,
              _ndivsnodes,
                _ndivscells,
                _override_vtkcelltype
      ) : error("Logic error!")
end

"""
  `get_node(grid, i)`

Returns the position of the i-th node (1-indexed) in the grid
"""
function get_node(self::Grid, i::Int64)
  if i>self.nnodes
    error("Requested invalid node index $i; max is $(self.nnodes).")
  elseif i<1
    error("Invalid index $i (it must be greater than 0).")
  end

  return self.nodes[:, i]
end

"""
  `get_node(grid, coor)`

Returns the position of the node of subscript coordinates `coor` (1-indexed)
"""
function get_node(self::Grid, coor::Array{Int64,1})
  return get_node(self, LinearIndices(self._ndivsnodes)[coor...])
end

"""
  `get_cell(grid, i)`

Returns the nodes indices of i-th cell in the grid (1-indexed)
"""
function get_cell(self::Grid, i::Integer)
  if i>self.ncells
    error("Requested invalid cell index $i; max is $(self.ncells).")
  elseif i<1
    error("Invalid index $i (it must be greater than 0).")
  end
  dims = Tuple(d != 0 ? d : 1 for d in self._ndivscells)
  # return get_cell(self, collect(ind2sub(self._ndivscells, i)))
  return get_cell(self, collect(Tuple(CartesianIndices(dims)[i])))
end

"""
  `get_cell(grid, coor)`

Returns the node indices of the cell with subscript coordinates `coor`
(1-indexed). The format corresponds to VTK_HEXAHEDRON (=12) in 3D, VTK_QUAD (=9)
in 2D, or VTK_LINE (=3) in 1D---except that points are 1-indexed instead of
0-indexed.
"""
function get_cell(self::Grid, coor_in::AbstractVector{<:Integer})

  # Correct 0 coordinate in quasi-dimensions
  coor = [ self._ndivscells[i]==0 && coor_in[i]==0 ? 1 : coor_in[i]
                                                    for i in 1:length(coor_in) ]
  # ERROR CASES
  if length(coor)!=self.dims
    error("$(self.dims)-dimensional grid requires $(self.dims) coordinates,"*
            " got $(length(coor)).")
  end
  for (dim, i) in enumerate(coor)
    if i>self._ndivscells[dim]
      if i==1 && self._ndivscells[dim]==0
        nothing
      else
        error("Requested cell $coor but max cell in"*
              " $dim-dimension is $(self._ndivscells[dim])")
      end
    end
  end

  # Displacements according to VTK convention, or closed loop
  disps = ones(Int64, self.dims)
  lpdm = self.loop_dim
  if lpdm!=0 && coor[lpdm]==self._ndivscells[lpdm]
    disps[lpdm] = -self._ndivscells[lpdm] + 1
  end
  d1, d2, d3 = vcat(disps, ones(Int64, 3-self.dims))

  # Checks for quasi-dimensions (NDIVS[i]==0)
  qdims = length( [1 for ndiv in self._ndivscells if ndiv==0] )
  cdims = length(coor)
  dims = self.dims - qdims
  aux1 = zeros(Int64, qdims)

  if dims==1
    return [LinearIndices(self._ndivsnodes)[(coor+vcat([0], aux1))...],
            LinearIndices(self._ndivsnodes)[(coor+vcat([d1], aux1))...]]

  elseif dims==2
    return [LinearIndices(self._ndivsnodes)[(coor+vcat([0,0], aux1))...],
            LinearIndices(self._ndivsnodes)[(coor+vcat([d1,0], aux1))...],
            LinearIndices(self._ndivsnodes)[(coor+vcat([d1,d2], aux1))...],
            LinearIndices(self._ndivsnodes)[(coor+vcat([0,d2], aux1))...]]

  elseif dims==3
    return [LinearIndices(self._ndivsnodes)[(coor+vcat([0,0,0], aux1))...],
            LinearIndices(self._ndivsnodes)[(coor+vcat([d1,0,0], aux1))...],
            LinearIndices(self._ndivsnodes)[(coor+vcat([d1,d2,0], aux1))...],
            LinearIndices(self._ndivsnodes)[(coor+vcat([0,d2,0], aux1))...],
            LinearIndices(self._ndivsnodes)[(coor+vcat([0,0,d3], aux1))...],
            LinearIndices(self._ndivsnodes)[(coor+vcat([d1,0,d3], aux1))...],
            LinearIndices(self._ndivsnodes)[(coor+vcat([d1,d2,d3], aux1))...],
            LinearIndices(self._ndivsnodes)[(coor+vcat([0,d2,d3], aux1))...]]

  else
    error("Definition of $(self.ndims)-dimensional cells not implemented yet!")
  end
end

"Plots the grid on PyPlot"
function plot(grid::Grid; fig_name="gridplot", fontsize=15,
                          xlims=nothing, ylims=nothing, zlims=nothing,
                          labelcells=true, labelnodes=false, labelndivs=true,
                          title_str=nothing,
                          alpha=1.0)

  if grid.dims>3
    error("There is no plotting method for $(grid.dims)-dimensional grids")
  end

  fig = PyPlot.figure(fig_name)
  ax = fig.gca(projection="3d")
  vectors_to_plot = []

  nc = grid.ncells

  # Iterates over every child plotting them
  for i in 1:nc
    nodes = [vcat(get_node(grid, node), zeros(Float64, 3-grid.dims))
                                                  for node in get_cell(grid, i)]

    if labelcells
      center = vcat(get_cellcenter(grid, i), zeros(Float64, 3-grid.dims))
      PyPlot.text3D(center[1], center[2], center[3], "$i", fontsize=fontsize, color="g")
    end

    # Connects the nodes for visualization
    dims = grid.dims - length( [1 for ndiv in grid._ndivscells if ndiv==0] )
    if dims==3
      for j in 1:3
        push!(vectors_to_plot, [nodes[j], -nodes[j]+nodes[j+1]])
        push!(vectors_to_plot, [nodes[4+j], -nodes[4+j]+nodes[4+j+1]])
      end
      push!(vectors_to_plot, [nodes[4], -nodes[4]+nodes[1]])
      push!(vectors_to_plot, [nodes[8], -nodes[8]+nodes[5]])
      for j in 1:4
        push!(vectors_to_plot, [nodes[0+j], -nodes[0+j]+nodes[4+j]])
      end
    elseif dims==2
      for j in 1:3
        push!(vectors_to_plot, [nodes[j], -nodes[j]+nodes[j+1]])
      end
      push!(vectors_to_plot, [nodes[4], -nodes[4]+nodes[1]])
    else
      push!(vectors_to_plot, [nodes[1], -nodes[1]+nodes[2]])
    end
  end

  # Iterates over every dimension labeling the subdivision
  if labelndivs
    for i in 1:grid.dims
      for j in 1:grid._ndivscells[i]
        coor = ones(Int64, grid.dims)
        coor[i] = j
        p1 = get_node(grid, coor)
        if grid.loop_dim==i && j==grid._ndivscells[i]
          coor[i] = 1
        else
          coor[i] += 1
        end
        p2 = get_node(grid, coor)
        center = vcat((p1+p2)/2, zeros(Float64, 3-grid.dims))
        PyPlot.text3D(center[1], center[2], center[3], "$j", fontsize=fontsize,
                                                                      color="r")
      end
    end
  end

  # Plots
  xyzuvw = [[] for _ in 1:6]
  for vector in vectors_to_plot
      _vector = []
      append!(_vector, vector[1])
      append!(_vector, vector[2])
      for (i,val) in enumerate(_vector)
          push!(xyzuvw[i], val)
      end
  end
  x,y,z,u,v,w = [coor for coor in xyzuvw]
  ax.quiver(x,y,z, u,v,w, arrow_length_ratio=0.0, alpha=alpha);

  # Labels nodes
  if labelnodes
    for i in 1:grid.nnodes
      pos = vcat(get_node(grid, i), zeros(Float64, 3-grid.dims))
      PyPlot.text3D(pos[1], pos[2], pos[3], "$i", fontsize=fontsize, color="k")
    end
  end

  # Format axes
  for (lim, lims, label, lbl) in [(PyPlot.xlim, xlims, PyPlot.xlabel, "x"),
                                  (PyPlot.ylim, ylims, PyPlot.ylabel, "y"),
                                  (PyPlot.zlim, zlims, PyPlot.zlabel, "z")]
    if lims!=nothing; lim(lims); end;
    label(lbl)

  end

  # Legend
  handles = []
  for (flag, clr, lbl) in [(labelnodes, "k", "Node"),
                            (labelcells, "g", "Cell"),
                            (labelndivs, "r", "Coordinate")]
    if flag; push!(handles, patch.Patch(color=clr, label=lbl)); end;
  end
  if size(handles,1)>0; PyPlot.legend(handles=handles); end;

  if title_str!=nothing; PyPlot.title(title_str); end;

end

"""
  `lintransform!(grid::Grid, M::Array{Float64,2}, T::Array{Float64,1})`

Rotates and translates the grid by the rotation matrix `M` and translation
vector `T` (linear transformation).
"""
function lintransform!(grid::Grid, M::Array{P,2}, T::Array{P,1};
                        reset_fields::Bool=true) where{P<:Real}
  # Error cases
  if size(M, 1)!=grid.dims
    error("Invalid rotation matrix `M`."*
            "Expected $(grid.dims) dimensions, got $(size(M,1))")
  elseif length(T)!=grid.dims
    error("Invalid translation vector `T`."*
            "Expected $(grid.dims) dimensions, got $(length(T))")
  end

  if reset_fields; grid.field = Dict{String, Dict{String, Any}}(); end;

  invM = inv(M)
  for i in 1:grid.nnodes
    grid.nodes[:,i] = countertransform(grid.nodes[:,i], invM, T)
  end
end

"""
  `transform!(grid::Grid, f)`

Applies the space transformation given by function `f` to the grid, where the
position of every node is given to the function `f`.
"""
function transform!(grid::Grid, f; reset_fields::Bool=true)

    if reset_fields; grid.field = Dict{String, Dict{String, Any}}(); end;

    for i in 1:grid.nnodes
      grid.nodes[:,i] .= f(view(grid.nodes, :, i))
    end
end

"""
  `transform2!(grid::Grid, f)`

Applies the space transformation given by function `f` to the grid, where the
indices of every node is given to the function `f`.
"""
function transform2!(grid::Grid, f; reset_fields::Bool=true)

    if reset_fields; grid.field = Dict{String, Dict{String, Any}}(); end;

    for i in 1:grid.nnodes
      grid.nodes[:,i] = f(CartesianIndices(grid._ndivsnodes)[i])
    end
end

"""
  `transform3!(grid::Grid, f)`

Applies the space transformation given by function `f` to the grid, where the
both the position and indices of every node is given to the function `f`.
"""
function transform3!(grid::Grid, f; reset_fields::Bool=true)

    if reset_fields; grid.field = Dict{String, Dict{String, Any}}(); end;

    for i in 1:grid.nnodes
      grid.nodes[:,i] = f(grid.nodes[:,i], CartesianIndices(grid._ndivsnodes)[i])
    end
end

##### INTERNAL FUNCTIONS  ######################################################
function _check(P_min::Array{T,1} where {T<:Real},
                P_max::Array{T,1} where {T<:Real} ,
                NDIVS::Array{T,1} where {T<:Any},
                loop_dim::Int64)

  # Error cases
  if size(P_min)!=size(P_max)
    error("`P_min` and `P_max` must have the same dimensions "*
                                          "($(size(P_min))!=$(size(P_max)))")
  elseif loop_dim>size(P_min,1) || loop_dim<0
    error("Invalid loop dimension $loop_dim"*
                                      " in $(size(P_min,1))-dimensional grid.")

  elseif typeof(NDIVS)==Array{Int64,1}
    if _calc_dims(P_min)!=length(NDIVS)
      error("Division for each dimension must be given "*
                                          "$(length(P_min))!=$(length(NDIVS))")
    end
    for (i, this_div) in enumerate(NDIVS)
      if this_div<0
        error("Invalid division $this_div in $i-dimension.")
      end
    end
  end
  _ndivs = _calc_ndivs(NDIVS)
  for i in 1:_calc_dims(P_min)
    if _ndivs[i]==0
      if P_min[i]!=P_max[i]
        error("Value in `P_min` must be equal to `P_max` in quasi-dimension $i"*
                                            "($(P_min[i])!=$(P_max[i]))")
      end
    else
      if P_min[i]>=P_max[i]
        error("Values in `P_min` must be less than values in `P_max` "*"
                                            ($(P_min[i])>=$(P_max[i]))")
      end
    end
  end

  return true
end

function _calc_ndivs(NDIVS::Array{T,1} where {T<:Any})
  if typeof(NDIVS)==Array{Int64,1}
    return deepcopy(NDIVS)

  elseif typeof(NDIVS)==Array{Array{Tuple{Float64,Int64,Float64,Bool},1},1}
    return [ sum([sec[2] for sec in secs]) for secs in NDIVS ]

  else
      error("`NDIVS` expected to be type `Array{Int64,1}` or"*
              " `Array{Array{Tuple{Float64,Int64,Float64,Bool},1},1}`; got "*
              "$(typeof(NDIVS))")
  end

end

"Calculates total number of nodes given NDIVS"
function _calc_nnodes(NDIVS::Array{T,1} where {T<:Any}, loop_dim::Int64)
  ndivs = _calc_ndivsnodes(NDIVS, loop_dim)
  return prod([ elem==0 ? 1 : elem for elem in ndivs])
end

"Calculates number of nodes in each dimension"
function _calc_ndivsnodes(NDIVS::Array{T,1} where {T<:Any}, loop_dim::Int64)
  ndivs = _calc_ndivs(NDIVS)
  if loop_dim!=0
    ndivs[loop_dim] -= 1
  end
  return ndivs .+ 1
end

"Calculates the number of cells given NDIVS"
function _calc_ncells(NDIVS::Array{T,1} where {T<:Any})
  return prod([ elem==0 ? 1 : elem for elem in _calc_ndivs(NDIVS)])
end

function _calc_dims(P::Array{T,1} where {T<:Real})
  return length(P)
end

function _generate_grid(P_min::Array{T,1} where {T<:Real},
                        P_max::Array{T,1} where {T<:Real},
                        NDIVS::Array{T,1} where {T<:Any},
                        loop_dim::Int64)

  dims = _calc_dims(P_min)
  nnodes = _calc_nnodes(NDIVS, loop_dim)
  nodes = Array{Float64}(undef, dims, nnodes)
  ndivs = Tuple(_calc_ndivs(NDIVS))

  # Discretizes each coordinate according to NDIVS
  if typeof(NDIVS)==Array{Int64,1}
    spacing = [range(P_min[i], stop=P_max[i], length=ndivs[i]+1) for i in 1:dims]

  elseif typeof(NDIVS)==Array{Array{Tuple{Float64,Int64,Float64,Bool},1},1}
    spacing = [multidiscretize(x->x, P_min[i], P_max[i], NDIVS[i]) for i in 1:dims]

  else
    error("Not implemented yet")
  end

  # Creates the grid
  ind2 = 1
  ndivsnodes = Tuple(_calc_ndivsnodes(NDIVS, 0))
  for ind in 1:_calc_nnodes(NDIVS, 0)
    sub = CartesianIndices(ndivsnodes)[ind]

    if loop_dim==0 || sub[loop_dim]!=ndivsnodes[loop_dim]
      p = [spacing[dim][sub[dim]] for dim in 1:dims]
      nodes[:, ind2] = p
      ind2 += 1
    end
  end

  return nodes
end

##### END OF GRID ##############################################################


text:
#Function to measue the maximum dimensions (in each dimensions) of the protein
  export measure_prot
  function measure_prot(nome, side)

  ## Variables - Measurement of the protein(X, Y and Z axis)
    pdbfile = nome
    segment = "all"
    restype = "all"
    atomtype = "all"
    firstatom = "all"
    lastatom = "all"

  ## Open the file -
    file = open("$pdbfile","r+")

    ## variables to compute the size of the protein
    natoms = 0
    var = 0

    for line in eachline(file)

      data = split(line)
      consider = false
      if data[1] == "ATOM" || data[1] == "HEATOM"
        natoms = natoms + 1
      #  ss = data[11]
      #  rt = data[4]
      #  at = data[3]
        consider = true
      else
        continue
      end

      if firstatom != "all"
        if natoms < firstatom
          consider = false
        end
      end

      if lastatom != "all"
        if natoms > lastatom
          consider = false
        end
      end

      if segment != "all"
        if ss != segment
          consider = false
        end
      end

      if restype != "all"
        if rt != restype
          consider = false
        end
      end

      if atomtype != "all"
        if at != atomtype
          consider = false
        end
      end

      if consider == true
        global x = parse(Float64,data[7])
        global y = parse(Float64,data[8])
        global z = parse(Float64,data[9])
      end

      if var == 1
        if x < xmin
            xmin = x
        end
        if y < ymin
            ymin = y
        end
        if z < zmin
            zmin = z
        end
        if x > xmax
            xmax = x
        end
        if y > ymax
            ymax = y
        end
        if z > zmax
            zmax = z
        end
      else

        global   xmin = x
        global   ymin = y
        global   zmin = z
        global   xmax = x
        global   ymax = y
        global   zmax = z
        var = 1

      end

    end

    ## Box size
    bx = round(Int,((xmax - xmin) + side + side))
    by = round(Int,((ymax - ymin) + side + side))
    bz = round(Int,((zmax - zmin) + side + side))

    close(file)

    cfile = open("$pdbfile","r+")

    # charge calculation
    charge = 0
    nhis   = 0
    narg   = 0
    nlys   = 0
    nglu   = 0
    nasp   = 0

    #mass=0

    marker = 0

    for line in eachline(cfile)
      data = split(line)
      if data[1] == "ATOM" || data[1] == "HEATOM"
        if data[4]=="HIS" || data[4]=="HSD" && marker != data[6]
          nhis = nhis + 1
          marker = data[6]
        elseif data[4]=="ARG" && marker != data[6]
          narg = narg + 1
          marker = data[6]
        elseif data[4]=="LYS" && marker != data[6]
          nlys = nlys + 1
          marker = data[6]
        elseif data[4]=="GLU" && marker != data[6]
          nglu = nglu + 1
          marker = data[6]
        elseif data[4]=="ASP" && marker != data[6]
          nasp = nasp + 1
          marker = data[6]
        end
      end
    end
    close(cfile)

    charge = (nhis * 0)  + (narg * (+1)) + (nlys * (+1)) + (nglu * (-1)) + (nasp * (-1))

    return bx , by, bz, charge

  end


text:
export WaveletOp

"""
  WaveletOp(shape, wt=wavelet(WT.db2))

returns a `ẀaveletOp <: AbstractLinearOperator`, which performs a Wavelet transform on
a given input array.

# Arguments

* `shape`                 - size of the Array to transform
* (`wt=wavelet(WT.db2)`)  - Wavelet to apply
"""
function WaveletOp(T::Type, shape, wt=wavelet(WT.db2))
  return LinearOperator{T}(prod(shape), prod(shape), false, false
            , x->vec( dwt(reshape(x,shape),wt) )
            , nothing
            , y->vec( idwt(reshape(y,shape),wt) ) )
end


text:
using SExpressions
using Test
using Base.Iterators

include("parser-util.jl")
include("parser.jl")
include("lists.jl")
include("keywords.jl")


text:
using QuadGK

"""
    AMFMcomp(in::Vector{Float64}, out::Vector{Complex{Float64}})
Creates an AMFM component with parameters `a`, `ω`, and `φ`.
    ψ =  a(t) exp[j(∫ω(τ)dt + φ)]

The input `t` must be a vector of Float64 and the out `ψ` will be a vector of
Complex{Float64} the same length as `t`.

```julia
julia> a₀(t)=1; ω₀(t)=1; φ₀ = 0; ψ₀ = AMFMcomp(a₀,ω₀,φ₀)
AMFMcomp(a₀, ω₀, 0)

julia> t = Array(0.0:0.5:1.0); ψ₀(t)
3-element Array{Complex{Float64},1}:
  1.0 + 0.0im
 -1.0 + 1.2246467991473532e-16im
  1.0 - 2.4492935982947064e-16im
```
"""

#Define type AMFMcomp as a struct of IA, IF, and phase reference
struct AMFMcomp
#References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.
  a::Function
  ω::Function
  φ::Real
end

#Unspecifed phase references are assumed to be zero
AMFMcomp(a, ω) = AMFMcomp(a, ω, 0.0)

#Method for type AMFMcomp when called with a time index
function (ψ::AMFMcomp)(t::Vector{Float64})
#References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.
  a, ω, φ = ψ.a, ψ.ω, ψ.φ
  A = quadgk.(ω, 0.0, t)
  θ = [ A[i][1] for i in 1:length(t) ]
  a.(t) .* exp.( 1im .* ( θ .+ φ ) )
end

function (ψ::AMFMcomp)(t::Float64)
  return (ψ::AMFMcomp)([t::Float64])[1]
end

function (ψ::AMFMcomp)(t::StepRangeLen)
  return (ψ::AMFMcomp)(collect(t))
end

#Construction of type AMFMcomp from a cannonical triplet
function AMFMcomp(triplet::Tuple{Function, Function, Real})
#References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.
      return AMFMcomp(triplet[1],triplet[2],triplet[3])
end

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA


text:
@testset "Transforms" begin
  c = Composition(1,2,3)
  @test alrinv(alr(c)) == c
  @test clrinv(clr(c)) == c
end


text:
include("utils.jl")
include("earthquake_main.jl")

disp, state, v_bd = simulate()

# mode = "data"
mode = "inv"

if mode == "data"

    init(sess)
    disp0, state0, v_bd0 = run(sess, [disp, state, v_bd]);
    matwrite("data-eq.mat", Dict("disp"=>disp0, "state"=>state0, "v_bd"=>v_bd0))

    close("all")
    figure(figsize=[8, 6])
    subplot(311)
    plot(t, v_bd0[1:end-1,1], "k", label="V")
    ylabel("Velosity (m/s)")
    legend()
    xlim([t[1], t[end]])
    # gca().set_ylim(bottom=1e-5)
    gca().get_xaxis().set_ticklabels([])
    subplot(312)
    plot(t, state0[1:end-1,1], "k", label=L"\Psi")
    ylabel("State Variable")
    legend()
    xlim([t[1], t[end]])
    gca().get_xaxis().set_ticklabels([])
    subplot(313)
    plot(t, disp0[1:end-1,1], "k", label="U")
    ylabel("Slip (earthquick_main)")
    legend()
    xlim([t[1], t[end]])
    xlabel("Time (s)")
    savefig("earthquake-simulation.png", bbox_inches="tight")

end

if mode == "inv"

    disp0 = matread("data-eq.mat")["disp"]
    state0 = matread("data-eq.mat")["state"]
    v_bd0 = matread("data-eq.mat")["v_bd"]
    loss = mean((disp[:,1:m+1]-disp0[:,1:m+1])^2) +
            mean((state-state0)^2) +
            mean((v_bd-v_bd0)^2)

    loss = 1e10*loss
    sess = Session(); init(sess)

    @show run(sess, loss)
    # BFGS!(sess, loss)

## DEBUG
# gradients(loss, a)
# gradients(loss, b)
# gradients(loss, σn)
# gradients(loss, v0)
# gradients(loss, Dc)
# gradients(loss, η)
# gradients(loss, f0)
# gradients(loss, bd_left0)

# lineview(sess, a, loss, ones(n+1) * 0.01, ones(n+1) * 0.02)
# lineview(sess, b, loss,ones(n+1) * 0.02, ones(n+1) * 0.03)
# lineview(sess, σn, loss, ones(n+1) * 50.0, ones(n+1) * 30.0)
# lineview(sess, Dc, loss, 0.03, 1.0)
# lineview(sess, f0, loss, 0.6, 0.8)
# lineview(sess, v0, loss, 1e-6, 1e-4)
# lineview(sess, η, loss, 4.7434, 10.0)
# lineview(sess, bd_left0, loss, ones(n+1) * 0.3, ones(n+1) * 0.4)

gradview(sess, a, loss, ones(n+1) * 0.02)
# gradview(sess, b, loss, ones(n+1) * 0.03)
# gradview(sess, σn, loss, ones(n+1) * 30.0)
# gradview(sess, Dc, loss, 1.0)
# gradview(sess, f0, loss, 0.8)
# gradview(sess, v0, loss, 1e-4)
# gradview(sess, η, loss, 10.0)
# gradview(sess, bd_left0, loss, ones(n+1) * 0.3)

# gradview(sess, pl, loss, [1.2])

# # run(sess, gradients(loss, σn), b=>0.018*ones(n+1)

end

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
using Plots
using DelimitedFiles
using LinearAlgebra
using ProgressMeter
using LsqFit
using KernelDensity
using Statistics
using LaTeXStrings
using StatsBase
using Measures
using JLD2
include("../Analysis/AlphaShapes/AlphaShapes.jl")
include("../Analysis/analysis-utils.jl")

# using CSV rather than Julia specific JLD format
X50 = readdlm("../Data/50-X.csv",',')
X50 = zeros(size(X50)...,3)
X50[:,:,1] = readdlm("../Data/50-X.csv",',')
X50[:,:,2] = readdlm("../Data/50-Y.csv",',')
X50[:,:,3] = readdlm("../Data/50-theta.csv",',')

X100 = readdlm("../Data/100-X.csv",',')
X100 = zeros(size(X100)...,3)
X100[:,:,1] = readdlm("../Data/100-X.csv",',')
X100[:,:,2] = readdlm("../Data/100-Y.csv",',')
X100[:,:,3] = readdlm("../Data/100-theta.csv",',')

X200 = readdlm("../Data/200-X.csv",',')
X200 = zeros(size(X200)...,3)
X200[:,:,1] = readdlm("../Data/200-X.csv",',')
X200[:,:,2] = readdlm("../Data/200-Y.csv",',')
X200[:,:,3] = readdlm("../Data/200-theta.csv",',');

JLD2.@load "../Data/alpha-tau-pd-15.jld"
# data is accesible from the global variable "Data" like
# Data["100--0.14285714285714285-0.0"]
# 3×100 Array{Float64,2}:
#  3.47599  3.47702  3.47817  3.47946  3.4809   …  0.0  0.0  0.0  0.0  0.0  0.0
#  3.91475  3.91589  3.91717  3.9186   3.92019     0.0  0.0  0.0  0.0  0.0  0.0
#  3.94935  3.9507   3.95221  3.9539   3.95578     0.0  0.0  0.0  0.0  0.0  0.0
# where each row is the density distribution of one run with N=100, alpha = -0.14285714285714285, and tau = 0.0

runs = 3
Ps = [50,100,200]
n = 15
# alpha and tau values used for the data
Alpha = collect(linspace(-2.0,0.0,n))
tau = collect(linspace(0.0,30.0,n))

bins_density = Float64.(logspace(1e-4,6.0,100))
bins_time = 100

BX = bins_density
B_pdfs = [TimeBinDensity(BX,100,X50,length_scale=beetle_length),
          TimeBinDensity(BX,100,X100[500:end,:,:],length_scale=beetle_length),
          TimeBinDensity(BX,100,X200,length_scale=beetle_length)]

@info "Computing Phase Diagrams"
@info "Looping over Alpha-Tau"
BX = bins_density
prog = Progress((n+1)*(n+1)*length(Ps))
for p in 1:length(Ps)
    plot(layout=(length(Alpha)+1)*(length(tau)+1),size=(10000,10000))
    plot!(xaxis=nothing,yaxis=nothing,bordercolor="white",subplot=1)
    # tau axis label
    for j in 1:n
        annotate!(0.5,0.5,"$(round(tau[j],digits=3))",xaxis=nothing,yaxis=nothing,bordercolor="white",100,subplot=j+1)
    end
    # alpha axis label
    for j in 1:n
        annotate!(0.5,0.5,"$(round(-Alpha[j],digits=3))",xaxis=nothing,yaxis=nothing,bordercolor="white",100,subplot=j*(n+1)+1)
    end
    i = n+2
    for a in Alpha
        for b in tau
            pdfs = Data["$(Ps[p])-$a-$b"]
            e = std(pdfs,dims=1)[1,:]
            pdfs = mean(pdfs,dims=1)[1,:]
            plot!(BX,pdfs,ribbon=e,label="",c="red",subplot=i+1,title="")
            plot!(BX,mean(B_pdfs[p],dims=1)[1,:],ribbon=std(B_pdfs[p],dims=1)[1,:],label="",c="blue",subplot=i+1,fillalpha=0.1)
            i += 1
            next!(prog)
            ylims!(0.0,Inf) # min 0.0, but max found as needed
            if p == 1
                xlims!(0,2.0)
            else
                xlims!(0,3.5)
            end
            #plot!(xaxis=nothing,yaxis=nothing,bordercolor="white",subplot=i)
        end
        i += 1
    end
    savefig("alpha-tau-pd-$(Ps[p]).png")
end

JLD2.@load "N-pd.jld"

bins_density = Float64.(logspace(1e-4,6.0,100))
BX = bins_density

plot(layout=(1,3),size=(1200,300))
for (i,p) in enumerate(Ps)
    pdfs = Data["N=$p"]["pdfs"]
    plot!(BX,mean(pdfs,dims=1)[1,:],ribbon=std(pdfs,dims=1)[1,:],label="",c="red",subplot=i)
    if p == 200
        B_pdfs = TimeBinDensity(BX,100,X200,length_scale=beetle_length)
        plot!(BX,mean(B_pdfs,dims=1)[1,:],ribbon=std(B_pdfs,dims=1)[1,:],label="",c="blue",subplot=i)
        ylims!(0.0,2.,subplot=i)
    elseif p == 100
        B_pdfs = TimeBinDensity(BX,100,X100,length_scale=beetle_length)
        plot!(BX,mean(B_pdfs,dims=1)[1,:],ribbon=std(B_pdfs,dims=1)[1,:],label="",c="blue",subplot=i)
        ylims!(0.0,2.,subplot=i)
    elseif p == 50
        B_pdfs = TimeBinDensity(BX,100,X50,length_scale=beetle_length)
        plot!(BX,mean(B_pdfs,dims=1)[1,:],ribbon=std(B_pdfs,dims=1)[1,:],label="",c="blue",subplot=i)
        ylims!(0.0,5.,subplot=i)
    end
    xlims!(0,3.)
end

savefig("N-pd.svg")
savefig("N-pd.png")

max_stars_repo_path: Model Analysis/model-analysis.jl

max_stars_repo_name: harveydevereux/CUDA-Whirligigs

max_issues_repo_path: Model Analysis/model-analysis.jl

max_issues_repo_name: harveydevereux/CUDA-Whirligigs

max_forks_repo_path: Model Analysis/model-analysis.jl

max_forks_repo_name: harveydevereux/CUDA-Whirligigs


text:
#!/usr/bin/env julia
module SCUTIWEN
import LightXML: parse_string, root,
    child_elements, find_element,
    content
import TextWrap: wrap

function main()
    uri = "https://api.forismatic.com/api/1.0/?" *
        "method=getQuote&format=xml&lang=en"
    response = read(`curl $uri -s`, String)
    let
        xml = (response
               |> parse_string
               |> root
               |> child_elements
               |> collect)[1]

        quote_text = (xml["quoteText"][1]
                      |> content
                      |> strip
                      |> xs ->
                          wrap(xs; width = 60,
                               subsequent_indent = 1))
        quote_author = content(xml["quoteAuthor"][1])

        println("\"\x1B[94m\x1B[1m$quote_text\x1B[0m\"")
        !isempty(quote_author) &&
            println("\x1B[93m$quote_author\x1B[0m")
    end # let
end # function

function __init__()
    main()
end # function
end # module


text:
function create_car_wcs(::Type{WCSTransform}, cdelt, crpix, crval)
    wcs = WCSTransform(2;
        ctype = ["RA---CAR", "DEC--CAR"],
        cdelt = collect(cdelt),
        crpix = collect(crpix),
        crval = collect(crval))

    return wcs
end

# try to follow WCS.WCSTransform with degree units
function create_car_wcs(::Type{CarClenshawCurtis{T}}, cdelt, crpix, crval) where T
    return CarClenshawCurtis{T}(
        (cdelt[1], cdelt[2]),
        (crpix[1], crpix[2]),
        (crval[1], crval[2]),
        π/180)  # degree conversion
end
create_car_wcs(::Type{CarClenshawCurtis}, cdelt, crpix, crval) =
    create_car_wcs(CarClenshawCurtis{Float64}, cdelt, crpix, crval)

"""
    fullsky_geometry([W=CarClenshawCurtis], res; shape = nothing, dims = ())

Generates a full-sky geometry.

# Arguments:
- `proj=CarClenshawCurtis()`: [optional] projection
- `res`: resolution in radians. Passing a Number produces a square pixel.
Passing a tuple with (ΔRA, ΔDEC) produces a rectangular pixel.

# Keywords
- `shape::NTuple=nothing`: shape of the map. If not specified, will be computed.
- `dims::NTuple=()`: additional dimensions to append to the shape, such as (3,) for IQU
to generate a map with `(nx, ny, 3)`.

# Returns:
- `shape::Tuple, wcs::W`: a tuple containing the shape of the map and the WCS

# Examples
```julia-repl
julia> shape, wcs = fullsky_geometry(deg2rad(30/60))  # 30 arcmin pixel
((720, 361), WCSTransform(naxis=2,cdelt=[-0.5, 0.5],crval=[0.25, 0.0],crpix=[360.5, 181.0]))
```
"""
function fullsky_geometry(W::Type{<:AbstractWCSTransform}, res; shape = nothing, dims = ())
    if isnothing(shape)
        shape = (round.(Int, (2π, π) ./ res .+ (0, 1)))  # CAR has pixels on poles
    end
    nx, ny = shape
    res_in_radians = ustrip.(uconvert.(radian, res))
    resx, resy = res_in_radians

    @assert abs(resx * nx - 2π) < 1e-8 "Horizontal resolution does not evenly divide the sky; this is required for SHTs."
    @assert abs(resy * (ny - 1) - π) < 1e-8 "Vertical resolution does not evenly divide the sky; this is required for SHTs."

    # Note the reference point is shifted by half a pixel to keep
    # the grid in bounds, from ra=180+cdelt/2 to ra=-180+cdelt/2.
    wcs = create_car_wcs(W,
        (-360.0 / nx, 180.0 / (ny - 1)),      # cdelt
        (floor(nx / 2) + 0.5, (ny + 1) / 2),  # crpix
        (resy * 90 / π, 0.0)                  # crval
    )

    return (nx, ny, dims...), wcs
end

fullsky_geometry(W::Type{<:AbstractWCSTransform}, res::Number; shape = nothing, dims = ()) =
    fullsky_geometry(W, (res, res); shape = shape, dims = dims)

fullsky_geometry(res; shape = nothing, dims = ()) =
    fullsky_geometry(CarClenshawCurtis{Float64}, res; shape = shape, dims = dims)

# ONLY DOES CAR FOR NOW
function geometry(W::Type{<:AbstractWCSTransform}, bbox_coords, res)

    # get everything into radians
    res_in_radians = ustrip.(uconvert.(radian, res))
    resx, resy = res_in_radians
    @assert abs(2π / resx - round(2π / resx)) < 1e-8 "Horizontal resolution" *
        " does not evenly divide the sky; this is required for SHTs."
    @assert abs(2π / resy - round(2π / resy)) < 1e-8 "Vertical resolution" *
        " does not evenly divide the sky; this is required for SHTs."

    # unpack the bounding box coordinates and set up the shape
    pos1 = ustrip.(uconvert.(radian, bbox_coords[:,1]))
    pos2 = ustrip.(uconvert.(radian, bbox_coords[:,2]))
    Δra, Δdec = abs.(pos1 .- pos2)
    shape = round.(Int, (Δra, Δdec) ./ res_in_radians)

    # set up coordinates such that the [1, 1] coordinate is the first pos of bbox
    mid = (pos1 .+ pos2) ./ 2
    crval = (mid[1], 0.0)
    cdelt = abs.(res_in_radians) .* sign.(pos2 .- pos1)
    crpix = 1 .- (pos1 .- crval) ./ cdelt

    wcs = create_car_wcs(W,
        rad2deg.(cdelt),  # cdelt
        crpix,            # crpix
        rad2deg.(crval))  # crval

    return shape, wcs
end

geometry(W::Type{<:AbstractWCSTransform}, bbox_coords, res::Number) =
    geometry(W, bbox_coords, (res, res))


text:
module TransmuteDims

export TransmutedDimsArray, transmute, transmutedims, transmutedims!, TupleVector

#========== alla PermutedDimsArrays, mostly ==========#

struct TransmutedDimsArray{T,N,perm,iperm,AA<:AbstractArray} <: AbstractArray{T,N}
    parent::AA
end

"""
    TransmutedDimsArray(A, perm⁺)

This is just like `PermutedDimsArray`, except that `perm⁺` need not be a permutation:

* Where it contains `0`, this inserts a trivial dimension into the output, size 1.
  Anything outside `1:ndims(A)` is treated like `0` -- fitting with `size(A,99) == 1`,
  but also allowing `nothing`.

* Any number omitted from `perm⁺` is a dimension to be dropped, which must be of size 1.

* When a number appears twice in `perm⁺`, the result works like `Diagonal` in those dimensions.

By calling the type directly you get the simplest constructor.
Calling instead [`transmute(A, perm⁺)`](@ref) gives one which un-wraps nested objects,
picks simpler alternatives, etc. And [`transmutedims(A, perm⁺)`](@ref) is the eager version.

# Examples

```jldoctest; setup=:(using TransmuteDims)
julia> TransmutedDimsArray('a':'e', (2,1))  # like transpose
1×5 transmute(::StepRange{Char, Int64}, (0, 1)) with eltype Char:
 'a'  'b'  'c'  'd'  'e'

julia> TransmutedDimsArray('a':'e', (nothing, 1, nothing))  # two trivial dimensions
1×5×1 transmute(::StepRange{Char, Int64}, (0, 1, 0)) with eltype Char:
[:, :, 1] =
 'a'  'b'  'c'  'd'  'e'

julia> A = rand(10, 20, 30);

julia> B = TransmutedDimsArray(A, (3,0,1,2));  # like reshape + permute

julia> size(B)
(30, 1, 10, 20)

julia> B[3,1,1,2] == A[1,2,3] == A[1,2,3,1]  # implicit trailing dimensions
true

julia> B == TransmutedDimsArray(A, (3,4,1,2)) == permutedims(A[:,:,:,:], (3,4,1,2))
true

julia> TransmutedDimsArray(reshape(1:6,3,2), (1,1,2))  # generalised Diagonal
3×3×2 transmute(reshape(::UnitRange{Int64}, 3, 2), (1, 1, 2)) with eltype Int64:
[:, :, 1] =
 1  ⋅  ⋅
 ⋅  2  ⋅
 ⋅  ⋅  3

[:, :, 2] =
 4  ⋅  ⋅
 ⋅  5  ⋅
 ⋅  ⋅  6
```
"""
function TransmutedDimsArray(data::AT, perm) where {AT <: AbstractArray{T,M}} where {T,M}
    P = sanitise_zero(perm, Val(M))
    Q = invperm_zero(P, size(data))
    TransmutedDimsArray{T,length(perm),P,Q,AT}(data)
end

getperm(A::TransmutedDimsArray{T,N,P,Q}) where {T,N,P,Q} = P
getinvperm(A::TransmutedDimsArray{T,N,P,Q}) where {T,N,P,Q} = Q

Base.parent(A::TransmutedDimsArray) = A.parent

Base.size(A::TransmutedDimsArray{T,N,P}) where {T,N,P} =
    genperm_zero(size(parent(A)), P)

Base.axes(A::TransmutedDimsArray{T,N,P}) where {T,N,P} =
    genperm_zero(axes(parent(A)), P, Base.OneTo(1))

Base.similar(A::TransmutedDimsArray, T::Type, dims::Base.Dims) = similar(parent(A), T, dims)

Base.unsafe_convert(::Type{Ptr{T}}, A::TransmutedDimsArray{T}) where {T} =
    Base.unsafe_convert(Ptr{T}, parent(A))

# It's OK to return a pointer to the first element, and indeed quite
# useful for wrapping C routines that require a different storage
# order than used by Julia. But for an array with unconventional
# storage order, a linear offset is ambiguous---is it a memory offset
# or a linear index?
Base.pointer(A::TransmutedDimsArray, i::Integer) = throw(ArgumentError(
    "pointer(A, i) is deliberately unsupported for TransmutedDimsArray"))

Base.strides(A::TransmutedDimsArray{T,N,P}) where {T,N,P} =
    genperm_zero(strides(parent(A)), P, 0)

Base.dataids(A::TransmutedDimsArray) = Base.dataids(A.parent)

Base.unaliascopy(A::TransmutedDimsArray) = typeof(A)(Base.unaliascopy(A.parent))

@inline function Base.getindex(A::TransmutedDimsArray{T,N,perm,iperm}, I::Vararg{Int,N}) where {T,N,perm,iperm}
    @boundscheck checkbounds(A, I...)
    val = @inbounds getindex(A.parent, genperm_zero(I, iperm, FirstIndex())...)
    if unique_or_zero(Val(perm))
        val
    else
        ifelse(is_off_diag(Val(perm), I), zero(T), val)
    end
end

@inline function Base.setindex!(A::TransmutedDimsArray{T,N,perm,iperm}, val, I::Vararg{Int,N}) where {T,N,perm,iperm}
    @boundscheck checkbounds(A, I...)
    if !unique_or_zero(Val(perm)) && is_off_diag(Val(perm), I)
        iszero(val) || throw(ArgumentError(
            "cannot set off-diagonal entry $I to a nonzero value ($val)"))
    else
        @inbounds setindex!(A.parent, val, genperm_zero(I, iperm, FirstIndex())...)
    end
    val
end

# # Not entirely sure this is a good idea, but passing KW along...
# Base.@propagate_inbounds Base.getindex(A::TransmutedDimsArray; kw...) =
#     getindex(A.parent; kw...)
# Base.@propagate_inbounds Base.setindex!(A::TransmutedDimsArray, val; kw...) =
#     setindex!(A.parent, val; kw...)

struct FirstIndex end

function Base.to_indices(A, tup, I::Tuple{FirstIndex, Vararg{Any}})
    if isempty(tup)
        (1, Base.to_indices(A, (), Base.tail(I))...)
    else
        (first(first(tup)), Base.to_indices(A, Base.tail(tup), Base.tail(I))...)
    end
end

#========== Utils ==========#

@inline function sanitise_zero(perm, ::Val{M}) where {M}
    map(n -> n isa Integer && 0<n<=M ? Int(n) : 0, Tuple(perm))
end

@inline function genperm_zero(I::Tuple, perm::Dims{N}, gap=1) where {N}
    ntuple(d -> perm[d]==0 ? gap : I[perm[d]], Val(N))
end

@inline function invperm_zero(P::NTuple{N,Int}, S::NTuple{M,Integer}) where {N,M}
    Q = ntuple(M) do d
        w = ntuple(n -> P[n]==d ? n : 0, N)
        x = max_zero(w...)
        x >= 1 || S[d]==1 || throw(ArgumentError(
            "dimension $d is missing from trasmutation $P, which is not allowed when size(A, $d) = $(S[d]) != 1"))
        x
    end
end

@inline max_zero() = 0
@inline max_zero(xs...) = max(xs...)

@inline function is_off_diag(P::Tuple, I::Tuple)
    for a in 1:length(P)
        for b in 1:length(P)
            P[a] == P[b] || continue
            I[a] == I[b] || return true
        end
    end
    return false
end
@inline @generated function is_off_diag(::Val{P}, I::Tuple) where {P}  # for getindex
    out = []
    for a in 1:length(P), b in 1:length(P)
        P[a] == P[b] && push!(out,  :( I[$a] != I[$b] ))
    end
    Expr(:call, :|, out...)
end

@inline function increasing_or_zero(tup::Tuple, prev=0)  # strictly increasing! hence a reshape
    d = first(tup)
    (d != 0) & (d <= prev) && return false
    return increasing_or_zero(Base.tail(tup), max(d, prev))
end
@inline increasing_or_zero(::Tuple{}, prev=0) = true

@inline function unique_or_zero(P::Tuple)
    first(P) != 0 && first(P) in Base.tail(P) && return false
    return unique_or_zero(Base.tail(P))
end
@inline unique_or_zero(P::Tuple{}) = true
@generated unique_or_zero(::Val{P}) where {P} = unique_or_zero(P)

#========== Printing ==========#

function Base.showarg(io::IO, A::TransmutedDimsArray{T,N,perm}, toplevel) where {T,N,perm}
    print(io, "transmute(")
    Base.showarg(io, parent(A), false)
    print(io, ", ", perm, ')')
    toplevel && print(io, " with eltype ", eltype(A))
end

# This is from  julia/stdlib/v1.3/LinearAlgebra/src/diagonal.jl:109
# function Base.replace_in_print_matrix(A::Diagonal,i::Integer,j::Integer,s::AbstractString)
function Base.replace_in_print_matrix(A::TransmutedDimsArray,
        i::Integer,j::Integer,s::AbstractString)
    is_off_diag(getperm(A), (i,j)) ? Base.replace_with_centered_mark(s) : s
end
function Base.replace_in_print_matrix(A::SubArray{<:Any,2,<:TransmutedDimsArray},
        i::Integer,j::Integer,s::AbstractString) where {S}
    ijk = (i, j, A.indices[3:end]...)
    is_off_diag(getperm(A.parent), ijk) ? Base.replace_with_centered_mark(s) : s
end

#========== Types ==========#

using LinearAlgebra

LazyTranspose{AT} = Union{
    Transpose{T,AT} where T<:Number,
    Adjoint{T,AT} where T<:Real,
}

getperm(A::LazyTranspose{<:AbstractVector}) = (0,1)
getperm(A::LazyTranspose{<:AbstractMatrix}) = (2,1)
getinvperm(A::LazyTranspose{<:AbstractVector}) = (2,)
getinvperm(A::LazyTranspose{<:AbstractMatrix}) = (2,1)

LazyPermute{P,AT} = Union{
    PermutedDimsArray{T,N,P,Q,AT} where {T,N,Q},
    TransmutedDimsArray{T,N,P,Q,AT} where {T,N,Q},
}

getperm(A::PermutedDimsArray{T,N,P}) where {T,N,P} = P
getinvperm(A::PermutedDimsArray{T,N,P,Q}) where {T,N,P,Q} = Q

may_reshape(::Type) = false
may_reshape(::Type{<:DenseArray}) = true
may_reshape(::Type{<:Base.ReshapedArray}) = true

#========== New constructor transmute() ==========#

"""
    transmute(A, perm⁺)
    transmute(A, Val(perm⁺))

Gives a result `== TransmutedDimsArray(A, perm⁺)`, but:

* When `A` is a `PermutedDimsArray`, `Transpose{<:Number}`, `Adjoint{<:Real}`, or `Diagonal`,
  it will un-wrap this, and adjust `perm⁺` to work directly on `parent(A)`.

* When the permutation simply inserts or removes a trivial dimension, it prefers to `reshape`
  instead of wrapping. This is controlled by `may_reshape(::Type)`, by default true for
  `DenseArray`s and `ReshapedArray`s.

* When possible, it prefers to create a `Transpose` or `Diagonal`, instead of a `TransmutedDimsArray`.

* If the permutation is wrapped in `Val`, then computing its inverse (and performing sanity checks)
  is always done at compile-time.

# Examples

```jldoctest; setup=:(using TransmuteDims, Random; Random.seed!(42);)
julia> A = transpose(rand(Int8, 4, 2))
2×4 transpose(::Matrix{Int8}) with eltype Int8:
 115    99  0  57
  88  -105  3  76

julia> B = transmute(A, (3,2,1))  # unwraps transpose, and then reshapes
1×4×2 Array{Int8, 3}:
[:, :, 1] =
 115  99  0  57

[:, :, 2] =
 88  -105  3  76

julia> B == TransmutedDimsArray(A, (0,2,1))  # same values, different type
true

julia> transmute(A, (2,2,0,1))
4×4×1×2 transmute(::Matrix{Int8}, (1, 1, 0, 2)) with eltype Int8:
[:, :, 1, 1] =
 115   ⋅  ⋅   ⋅
   ⋅  99  ⋅   ⋅
   ⋅   ⋅  0   ⋅
   ⋅   ⋅  ⋅  57

[:, :, 1, 2] =
 88     ⋅  ⋅   ⋅
  ⋅  -105  ⋅   ⋅
  ⋅     ⋅  3   ⋅
  ⋅     ⋅  ⋅  76

julia> ans == transmute(B, (2,2,1,3))
true

julia> TransmuteDims.may_reshape(typeof(A))  # avoids reshaping wrappers
false

julia> TransmuteDims.may_reshape(typeof(B))
true

julia> transmute(('α', "β", :γ), (2,1))  # accepts tuples, too
1×3 transmute(TupleVector(::Tuple{Char, String, Symbol}), (0, 1)) with eltype Any:
 'α'  "β"  :γ
```
"""
function transmute end

# @static if VERSION > v"1.7.0-DEV.400" # not exactly right!
#     @inline Base.@aggressive_constprop transmute(data::AbstractArray, perm) = _transmute(data, sanitise_zero(perm, Val(ndims(data))))
# else
    @inline transmute(data::AbstractArray, perm) = _transmute(data, sanitise_zero(perm, Val(ndims(data))))
# end

# First dispatch is on perm, second is _transmute(::SomeArray, ::Any),
# un-wrapping all before running this:

@inline function _transmute(data::AT, P) where {AT <: AbstractArray{T,M}} where {T,M}
    N = length(P)
    Q = invperm_zero(P, size(data))
    if M == N && P == ntuple(identity, N)  # trivial case
        data
    elseif may_reshape(AT) && increasing_or_zero(P)
        S = map(d -> d==0 ? Base.OneTo(1) : axes(data,d), P)
        reshape(data, S)
    elseif M == 2 && P == (2,1) && T<:Number
        transpose(data)
    elseif M == 1 && P == (1,1)
        Diagonal(data)
    elseif may_reshape(AT) && N == 2 && P[1] == P[2]
        Diagonal(vec(data))
    else
        TransmutedDimsArray{T,length(P),P,Q,AT}(data)
    end
end

for _fun in [:_transmute, :_transmutedims]
    @eval begin

        @inline function $_fun(data::LazyTranspose{<:AbstractVector}, perm)
            new_perm = map(d -> d==1 ? 0 : d==2 ? 1 : d, perm)
            $_fun(parent(data), new_perm)
        end

        @inline function $_fun(data::LazyTranspose{<:AbstractMatrix}, perm)
            new_perm = map(d -> d==1 ? 2 : d==2 ? 1 : d, perm)
            $_fun(parent(data), new_perm)
        end

        @inline function $_fun(data::LazyPermute{inner}, perm) where {inner}
            new_perm = map(d -> d==0 ? 0 : inner[d], perm)
            $_fun(parent(data), new_perm)
        end

        @inline function $_fun(data::Diagonal, perm)
            new_perm = map(d -> d==1 ? 1 : d==2 ? 1 : 0, perm)
            $_fun(parent(data), new_perm)
        end

    end
end

#========== Version with Val(perm) ==========#

@generated function transmute(data::AbstractArray, ::Val{perm}) where {perm}
    _trex(:data, data, sanitise_zero(perm, Val(ndims(data))))
end

# Identical list of cases:

function _trex(ex, AT, P)
    M = ndims(AT)
    N = length(P)
    if M == N && P == 1:N  # trivial case
        return ex
    end
    @gensym sym
    T = eltype(AT)
    Q, checks = invperm_zero(P, M, sym)

    noreshape = if M == 2 && P == (2,1) && T<:Number
        :(transpose($sym))
    elseif M == 1 && P == (1,1) # "may_reshape(AT) && N == 2 && P[1] == P[2]" will be messy.
        :(Diagonal($sym))
    else
        :(TransmutedDimsArray{$T,$(length(P)),$P,$Q,$AT}($sym))
    end

    if increasing_or_zero(P)
        S = map(d -> d==0 ? :(Base.OneTo(1)) : :(axes($sym,$d)), P)
        withreshape = :(reshape($sym, ($(S...),)))
        return quote
            $sym = $ex
            $(checks...)
            if may_reshape($AT) # check this trait at run-time
                $withreshape
            else
                $noreshape
            end
        end
    else
        return quote
            $sym = $ex
            $(checks...)
            $noreshape
        end
    end
end

function invperm_zero(P::NTuple{N,Int}, M::Int, sym::Symbol) where {N}
    checks = []
    Q = ntuple(M) do d
        w = ntuple(n -> P[n]==d ? n : 0, N)
        str = "dimension $d is missing from transmutation $P, which is not allowed when size(A, $d) = != 1"
        x = max_zero(w...)
        x >= 1 || push!(checks, quote
            size($sym,$d)==1 || throw(ArgumentError($str))
        end)
        x
    end
    Q, checks
end

function _trex(ex, ::Type{AT}, perm) where {AT<:LazyTranspose{PT}} where {PT<:AbstractVector}
    new_perm = map(d -> d==1 ? 0 : d==2 ? 1 : d, perm)
    _trex(:(parent($ex)), PT, new_perm)
end

function _trex(ex, ::Type{AT}, perm) where {AT<:LazyTranspose{PT}} where {PT<:AbstractMatrix}
    new_perm = map(d -> d==1 ? 2 : d==2 ? 1 : d, perm)
    _trex(:(parent($ex)), PT, new_perm)
end

function _trex(ex, ::Type{AT}, perm) where {AT<:LazyPermute{inner,PT}} where {inner,PT}
    new_perm = map(d -> d==0 ? 0 : inner[d], perm)
    _trex(:(parent($ex)), PT, new_perm)
end

function _trex(ex, ::Type{AT}, perm) where {AT<:Diagonal{T,PT}} where {T,PT}
    new_perm = map(d -> d==1 ? 1 : d==2 ? 1 : 0, perm)
    _trex(:(parent($ex)), PT, new_perm)
end

#========== More eager transmutedims() ==========#

"""
    transmutedims(A, perm⁺)

This is an eager version of [`transmute`](@ref), which always returns a `DenseArray`,
but is not guaranteed to copy the data.

Like `transmute`, it knows to un-wrap `PermutedDimsArray`, `Transpose{<:Number}`, etc.

Like both `transmute` and [`TransmutedDimsArray`](@ref), it accepts in addition to perturbations,
values outside `1:ndims(A)` (which insert trivial dimensions), omitted values
(which like `dropdims` must be dimensions of size 1), and repeated values (which generalise `diagm`).

# Examples

```jldoctest; setup=:(using TransmuteDims)
julia> A = transmutedims(reshape(1:15,3,5), (3,2,1))  # A is a new Array
1×5×3 Array{Int64, 3}:
[:, :, 1] =
 1  4  7  10  13

[:, :, 2] =
 2  5  8  11  14

[:, :, 3] =
 3  6  9  12  15

julia> B = transmutedims(A, (2,3))  # drop A's first, trivial, dimension
5×3 Matrix{Int64}:
  1   2   3
  4   5   6
  7   8   9
 10  11  12
 13  14  15

julia> C = transmutedims(adjoint(B), (2,1,0)); summary(C)  # un-wraps adjoint
"5×3×1 Array{Int64, 3}"

julia> D = transmutedims(adjoint(B), (1,0,2)); summary(D)  # unwraps, then permutes
"3×1×5 Array{Int64, 3}"

julia> B[5,3] = 5030;  # B and C are reshaped views of A

julia> A[1,5,3]
5030

julia> C[5,3,1]
5030

julia> D[3,1,5]  # but D is not, it required a copy
15

julia> transmutedims((1, 2, 3.0, 4+5im))  # default perm = (0,1) for vectors (and tuples)
1×4 Matrix{Number}:
 1  2  3.0  4+5im

julia> transmutedims(B)  # default perm = (2,1) for matrices
3×5 Matrix{Int64}:
 1  4  7  10    13
 2  5  8  11    14
 3  6  9  12  5030
```
"""
@inline function transmutedims(data::AbstractArray, perm)
    P = sanitise_zero(perm, Val(ndims(data)))
    _transmutedims(data, P)
end

@inline transmutedims(data::AbstractMatrix) = _transmutedims(data, (2,1))

@inline transmutedims(data::AbstractVector) = _transmutedims(data, (0,1))

@inline function _transmutedims(data::AT, P) where {AT <: AbstractArray{T,M}} where {T,M}
    N = length(P)
    Q = invperm_zero(P, size(data))
    S = map(d -> d==0 ? Base.OneTo(1) : axes(data,d), P)
    if may_reshape(AT) && M == N && P == ntuple(identity, N)  # trivial case demands may_reshape
        data
    elseif may_reshape(AT) && increasing_or_zero(P)
        reshape(data, S)
    else
        copy(TransmutedDimsArray{T,N,P,Q,AT}(data))
    end
end

# This exists for testing, with integer data it calls @generated _densecopy_permuted!
function _float_transmutedims(data::AT, perm) where {AT <: AbstractArray{T,M}} where {T,M}
    P = sanitise_zero(perm, Val(ndims(data)))
    N = length(P)
    Q = invperm_zero(P, size(data))
    S = map(d -> d==0 ? Base.OneTo(1) : axes(data,d), P)
    if may_reshape(AT) && M == N && P == ntuple(identity, N)
        data
    elseif may_reshape(AT) && increasing_or_zero(P)
        reshape(data, S)
    else
        out = similar(data, float(T), S)
        copy!(out, TransmutedDimsArray{T,N,P,Q,AT}(data))
    end
end

"""
    transmutedims!(dst, src, perm⁺)

This is just `copy!(dst, transmute(src, perm⁺))`.
"""
@inline function transmutedims!(dst, src::AbstractArray, perm)
    P = sanitise_zero(perm, Val(ndims(src)))
    copy!(dst, _transmute(src, P))  # using _transmute might do one more reshape than ideal
end

#========== The rest ==========#

include("base.jl")

include("strided.jl") # this costs about 0.1s

include("chainrules.jl") # this costs about 0.2s

using Requires # without loading anything, this costs about 0.1s

@init @require GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7" begin
    # GPUArrays is still in Project.toml, to constrain versions
    include("gpu.jl") # this costs just over a second
end

@init @require OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" begin
    using .OffsetArrays
    may_reshape(::Type{<:OffsetArray{T,N,AT}}) where {T,N,AT} = may_reshape(AT)
end

@init @require Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
    include("tracker.jl")
end

end


text:
using ConstrainedDynamics
using LinearAlgebra
using StaticArrays

# Parameters
h = .1
r = 1.
b1 = Cylinder(r, h, h*r, color = RGBA(1., 0., 0.))

# Links
origin = Origin{Float64}()
link1 = Body(b1)

# Constraints
joint1 = EqualityConstraint(Floating(origin, link1))

links = [link1]
constraints = [joint1]
shapes = [b1]

mech = Mechanism(origin, links, constraints, g = 0., shapes = shapes)

axis = [0;0;1.]
speed = 2pi #*0
setVelocity!(link1, ω = speed*axis)

function nutation_control!(mechanism, k)
    if k==1
        setForce!(mechanism.bodies[1], F = SA[0;0;0.2] * 0, p = SA[0;1.;0])
    else
        setForce!(mechanism.bodies[1], F = szeros(3), p = szeros(3))
    end
    return
end


text:
using ProblemParser
using Utils

using DataStructures
using Dictionaries
using SplitApplyCombine

import IterTools

const file = find_input(@__FILE__)
const input = parse(FirstRest(Blocks(),
        nothing,
        LineMappings(Split(" -> "), Apply(Tuple), Apply(c -> c[begin])),
    ), slurp(file))

function solve(input, rounds)
    line = input[1]
    pairs = DefaultDict{Tuple{Char,Char},Int}(0)
    for tup in IterTools.partition(line, 2, 1)
        pairs[tup] = pairs[tup] = pairs[tup] + 1
    end
    for _ in 1:rounds
        next = DefaultDict{Tuple{Char,Char},Int}(0)
        for (pair, count) in pairs
            middle = input[2][pair]
            tup = pair[1], middle
            next[tup] = next[tup] + count
            tup = middle, pair[2]
            next[tup] = next[tup] + count
        end
        pairs = next
    end
    letters = groupsum(tup -> tup[1][1], tup -> tup[2], pairs)
    set!(letters, line[end], get(letters, line[end], 0) + 1)
    mi, ma = extrema(letters)
    ma - mi
end

assertequal(solve(input, 10), 3048)

assertequal(solve(input, 40), 3288891573057)

max_stars_repo_name: CmdQ/AoC


text:
using Pkg
using SafeTestsets
using Test
using SciMLBase

@test isempty(detect_ambiguities(SciMLBase))

const GROUP = get(ENV, "GROUP", "All")
const is_APPVEYOR = ( Sys.iswindows() && haskey(ENV,"APPVEYOR") )

@time begin
@time @safetestset "Existence functions" begin include("existence_functions.jl") end
@time @safetestset "Integrator interface" begin include("integrator_tests.jl") end
end


text:
"""
	MCSMethod

Abstract type for nesting the various methods that can be used to perform a model confidence set. Subtypes are: \n
	MCSBoot
	MCSBootLowRAM
The subtypes have entries in the help (?) menu.
"""
abstract type MCSMethod end

"""
	MCSBoot(data ; alpha, basecaseindex, kwargs...)

Method type for doing a MCS test using a dependent bootstrap. Relevant keyword arguments are: \n
    alpha=0.05 <- The confidence level to use in the test \n
	basecaseindex=1 <- The MCS does not have a natural basecase. However, the input data is a
		matrix of losses, not loss differentials, and the method itself considers every possible
		loss differential combination. This is far too many to check for determining
		the optimal block length for bootstrapping, so users who wish to auto-estimate
		the appropriate block length can specify a notional base case variable index
		for the purposes of computing a optimal block lengths of all loss differences
		relative to the base case.
    blocklength=0.0 <- Block length to use with the bootstrap. Default of 0.0 implies
        the block length will be optimally estimated from the data use the method deemed
        most appropriate by the DependentBootstrap package (typically the selection procedure
        of Patton, Politis, and White (2009)). \n
    numresample=1000 <- Number of resamples to use when bootstrapping. \n
    bootmethod=:stationary <- Bootstrap methodology to use, where the default is the
        stationary bootstrap of Politis and Romano (1993) \n
Note, see the DependentBootstrap docs for more information on valid keyword arguments
for the BootInput constructor.
"""
struct MCSBoot <: MCSMethod
	alpha::Float64
	bootinput::BootInput
	basecaseindex::Int
	function MCSBoot(alpha::Float64, bootinput::BootInput, basecaseindex::Int)
		!(0.0 < alpha < 0.5) && error("alpha = $(alpha). Invalid confidence level.")
		new(alpha, bootinput, basecaseindex)
	end
end
function MCSBoot(data ; alpha::Float64=0.05, basecaseindex::Int=1, kwargs...)::MCSBoot
    return MCSBoot(alpha, BootInput(loss_diff_base_case(data, basecaseindex), kwargs...), basecaseindex)
end

"""
	MCSBootLowRAM(data ; alpha::Float64=0.05, kwargs...)

Method type for doing an MCS test using a dependent bootstrap. This method type has an identical
constructor arguments to MCSBoot, so use ?MCSBoot for more detail.

WARNING: This method corresponds to a MCS algorithm that has double the runtime of MCSBoot, but uses about half as much RAM. The
vast majority of users will want to use MCSBoot.

WARNING: Results from MCSBootLowRAM are not guaranteed to be identical to those from MCSBoot in finite sample, although they
are very likely to both recommend the same set of models in the model confidence set.
"""
struct MCSBootLowRAM <: MCSMethod
	alpha::Float64
	bootinput::BootInput
	basecaseindex::Int
	function MCSBootLowRAM(alpha::Float64, bootinput::BootInput, basecaseindex::Int)
		!(0.0 < alpha < 0.5) && error("alpha = $(alpha). Invalid confidence level.")
		new(alpha, bootinput, basecaseindex)
	end
end
function MCSBootLowRAM(data ; alpha::Float64=0.05, basecaseindex::Int=1, kwargs...)::MCSBootLowRAM
    return MCSBootLowRAM(alpha, BootInput(loss_diff_base_case(data, basecaseindex), kwargs...), basecaseindex)
end

"loss_diff_base_case <- Local function for getting loss differences relative to notional base case for purposes of estimating the block length"
function loss_diff_base_case(l::Matrix{<:Number}, basecaseindex::Int)
	!(1 <= basecaseindex <= size(l, 2)) && error("basecaseindex = $(basecaseindex), which does not lie between 1 and $(size(l,2)) (inclusive)")
	lD = l[:, 1:size(l, 2) .!= basecaseindex] .- l[:, basecaseindex]
	return lD
end

"""
	MCSTest(inQF::Vector{Int}, outQF::Vector{Int}, pvalueQF::Vector{Float64}, inMT::Vector{Int}, outMT::Vector{Int}, pvalueMT::Vector{Float64})

Output type from performing the MCS test proposed in Hansen, Lunde, Nason (2011) "The Model Confidence Set", Econometrica, 79 (2), pp. 453-497. \n
The field names of this type have trailing "QF" or "MT", where QF corresponds to the quadratic form test (see section 3.1.1 of original paper),
while MT corresponds to the maximum t-stat test (see section 3.1.2 of original paper). \n
Note, in the fields of this type, the forecast models input to the MCS method are indicated via an integer counting up from 1 for
each forecast model. These integers correspond to the columns of the input loss matrix; see ?mcs for more info. \n
The fields of this type follow: \n
	inQF <- Models that are in the MCS via the quadratic form method
	outQF <- Models that are not in the MCS via the quadratic form method
	pvalueQF <- The cumulative p-values from the quadratic form method
	inMT <- Models that are in the MCS via the max t-stat method
	outMT <- Models that are not in the MCS via the max t-stat method
	pvalueMT <- The cumulative p-values from the max t-stat method
"""
struct MCSTest #Output from MCS
	inQF::Vector{Int}
	outQF::Vector{Int}
	pvalueQF::Vector{Float64}
	inMT::Vector{Int}
	outMT::Vector{Int}
	pvalueMT::Vector{Float64}
	MCSTest(inQF::Vector{Int}, outQF::Vector{Int}, pvalueQF::Vector{Float64}, inMT::Vector{Int}, outMT::Vector{Int}, pvalueMT::Vector{Float64}) = new(inQF, outQF, pvalueQF, inMT, outMT, pvalueMT)
end
function Base.show(io::IO, a::MCSTest)
	println("Model confidence test results:")
	println("    Quadratic form test models in MCS: $(a.inQF)")
	println("    Maximum t-stat test models in MCS: $(a.inMT)")
end

"""
	mcs(l::Matrix{<:Number}, method ; kwargs...)
	mcs(l::Matrix{<:Number} ; kwargs...)

This function implements the MCS test proposed in Hansen, Lunde, Nason (2011) "The Model Confidence Set", Econometrica, 79 (2), pp. 453-497. \n
Let x_k, k = 1, ..., K, denote K forecasts (from K different forecasting models), and y denote the forecast target.
Let L(., .) denote a loss function.
The first argument of mcs is a matrix where the kth column of the matrix is created by the operation: \n
L(x_k, y) \n
Note that unlike the Reality Check and SPA test, there is no base case for the MCS. \n
The second method argument determines which methodology to use. Currently, only MCSBoot is available and if
this input type is provided, the keyword arguments are not needed.
Alternatively, the user can omit the second argument, and then any keyword arguments will be passed to the
MCSBoot constructor. See ?MCSBoot for more detail. The most relevant keyword arguments are:
    alpha=0.05 <- The confidence level to use in the test \n
	basecaseindex=1 <- The MCS does not have a natural basecase. However, the input data is a
		matrix of losses, not loss differentials, and the method itself considers every possible
		loss differential combination. This is far too many to check for determining
		the optimal block length for bootstrapping, so users who wish to auto-estimate
		the appropriate block length can specify a notional base case variable index
		for the purposes of computing a optimal block lengths of all loss differences
		relative to the base case.
    blocklength=0.0 <- Block length to use with the bootstrap. Default of 0.0 implies
        the block length will be optimally estimated from the data use the method deemed
        most appropriate by the DependentBootstrap package (typically the selection procedure
        of Patton, Politis, and White (2009)). \n
    numresample=1000 <- Number of resamples to use when bootstrapping. \n
    bootmethod=:stationary <- Bootstrap methodology to use, where the default is the
        stationary bootstrap of Politis and Romano (1993) \n
For more detail on the bootstrap options, please see the docs for the DependentBootstrap package. \n
The output of a MCS test is of type MCSTest. Use ?MCSTest for more information. \n
Note, if you are hitting RAM limits, type ?ForecastEval.MCSBootLowRAM at the REPL for more detail on an alternative algorithm that
is also available. \n
Note, for any developers, the main algorithm (associated with MCSBoot) still has the following potential performance issues: \n
	ISSUE 1: Some of the temporary arrays in the loops could probably be eliminated \n
	ISSUE 2: For MCS method A, I think the loop over K could be terminated as soon as cumulative p-values are greater than method.alpha. Need to double check this. \n
	ISSUE 3: Need to add option to do just max(abs) method or just sum(sq) method (or both) \n
Comments or pull requests on these issues would be most welcome on the package github page.
"""
function mcs(l::Matrix{T}, method::MCSBoot)::MCSTest where {T<:Number}
	(N, K) = size(l)
	N < 2 && error("Input must have at least two observations")
	K < 2 && error("Input must have at least two models")
	T != Float64 && (l = Float64[ Float64(l[j, k]) for j = 1:size(l, 1), k = 1:size(l, 2) ]) #Convert input to Float64
	numResample = method.bootinput.numresample
	#Get bootstrap indices
	inds = dbootinds(Float64[], method.bootinput)
	#Get matrix of loss differential sample means
	lMuVec = mean(l, dims=1)
	lDMu = Float64[ lMuVec[k] - lMuVec[j] for j = 1:K, k = 1:K  ]
	#Get array of  bootstrapped loss differential sample means
	lDMuStar = Array{Float64}(undef, K, K, numResample) #This array is affected by ISSUE 1 above
	for m = 1:numResample
		lMuVecStar = mean(l[inds[m], :], dims=1)
		lDMuStar[:, :, m] = Float64[ lMuVecStar[k] - lMuVecStar[j] for j = 1:K, k = 1:K  ]
	end
	#Get variance estimates from bootstrapped loss differential sample means (note, we centre on lDMu since these are the population means for the resampled data)
	#Note, for efficiency, we only use varm to fill out the lower triangular matrix
	lDMuVar = ones(Float64, K, K)
	for j = 2:K
		for k = 1:j-1
			lDMuVar[j, k] = varm(vec(lDMuStar[j, k, :]), lDMu[j, k], corrected=false)
		end
	end
	ltri_to_utri!(lDMuVar)
	#Get original and re-sampled t-statistics
	tStatStar = Float64[ (lDMuStar[j, k, m] - lDMu[j, k]) / sqrt(lDMuVar[j, k]) for j = 1:K, k = 1:K, m = 1:numResample ]
	tStat = Float64[ lDMu[j, k] / sqrt(lDMuVar[j, k]) for j = 1:K, k = 1:K ]
	#Perform model confidence set method A
	inA = collect(1:K) #Models in MCS (start off with all models in MCS)
	outA = Array{Int}(undef, K) #Models not in MCS (start off with no models in MCS)
	pValA = ones(Float64, K) #p-values constructed in loop
	for k = 1:K-1
		bootSum = 0.5 * vec(sum(abs2, tStatStar[inA, inA, :], dims=(1,2)))
		origSum = 0.5 * sum(abs2, tStat[inA, inA])
		pValA[k] = mean(bootSum .> origSum)
		scalingTerm = length(inA) / (length(inA) - 1)
		lDAvgMu = scalingTerm * vec(mean(lDMu[inA, inA], dims=1))
		lDAvgMuStar = scalingTerm * dropdims(mean(lDMuStar[inA, inA, :], dims=1), dims=1)
		lDAvgMuVar = Float64[ varm(vec(lDAvgMuStar[k, :]), lDAvgMu[k], corrected=false) for k = 1:length(lDAvgMu) ]
		tStatInc = lDAvgMu ./ sqrt.(lDAvgMuVar)
		iRemove = argmax(tStatInc) #Find index in inA of model to be removed
		outA[k] = inA[iRemove] #Add model to be removed to excluded list
		deleteat!(inA, iRemove) #Remove model to be removed
	end
	pValA = accumulate(max, pValA)
	outA[end] = inA[1] #Finish constructing excluded models
	iCutOff = findfirst(pValA .>= method.alpha) #method.alpha < 1.0, hence there will always be at least one p-value > method.alpha
	inA = outA[iCutOff:end]
	outA = outA[1:iCutOff-1]
	#Perform model confidence method B
	inB = collect(1:K) #Models in MCS (start off with all models included)
	outB = Array{Int}(undef, K) #Models not in MCS (start off with no models in MCS)
	pValB = ones(Float64, K) #p-values constructed in loop
	for k = 1:K-1
		bootMax = Float64[ maximum(abs, tStatStar[inB, inB, m]) for m = 1:numResample ]
		#bootMax = vec(maximum(abs(tStatStar[inB, inB, :]), [1, 2])) #ALTERNATIVE METHOD
		origMax = maximum(tStat[inB, inB])
		pValB[k] = mean(bootMax .> origMax)
		scalingTerm = length(inB) / (length(inB) - 1)
		lDAvgMu = scalingTerm * vec(mean(lDMu[inB, inB], dims=1))
		lDAvgMuStar = scalingTerm * dropdims(mean(lDMuStar[inB, inB, :], dims=1), dims=1)
		lDAvgMuVar = Float64[ varm(vec(lDAvgMuStar[k, :]), lDAvgMu[k], corrected=false) for k = 1:length(lDAvgMu) ]
		tStatInc = lDAvgMu ./ sqrt.(lDAvgMuVar)
		iRemove = argmax(tStatInc) #Find index in inB of model to be removed
		outB[k] = inB[iRemove] #Add model to be removed to excluded list
		deleteat!(inB, iRemove) #Remove model to be removed
	end
	pValB = accumulate(max, pValB)
	outB[end] = inB[1] #Finish constructing excluded models
	iCutOff = findfirst(pValB .>= method.alpha) #method.alpha < 1.0, hence there will always be at least one p-value > method.alpha
	inB = outB[iCutOff:end]
	outB = outB[1:iCutOff-1]
	#Prepare the output
	mcsOut = MCSTest(inA, outA, pValA, inB, outB, pValB)
end
#Local function to shift lower triangular elements to upper triangular portion of matrix
function ltri_to_utri!(x::AbstractMatrix{T}) where {T<:Number}
	size(x, 1) != size(x, 2) && error("Input matrix  must be symmetric")
	for j = 2:size(x, 1)
		for k = 1:j-1
			x[k, j] = x[j, k]
		end
	end
	return(true)
end
#Keyword wrapper
mcs(l::Matrix{<:Number} ; kwargs...)::MCSTest = mcs(l, MCSBoot(l, kwargs...))


text:
using JuliennedArrays
using Documenter: deploydocs, makedocs

makedocs(sitename = "JuliennedArrays.jl", modules = [JuliennedArrays], doctest = false)
deploydocs(repo = "github.com/bramtayl/JuliennedArrays.jl.git")