text:
include("neighbor.jl")
include("sortedOffsets.jl")
include("result/synthResult.jl")

#=
A patch that can be neither rectangular nor contiguous.
That is, it can be an irregular shape and sparse.

Shape
=======

Irregular i.e. non-centered

 Here continuguous but skewed with distinguished point in corner
 XXX
 XXC

Here, non-contiguous

X  X  X
X
X  C
X XXXXX

Here, non-rectangular but nearest to distinguished point.

 X
XCX
 X

Length
======

All ScatterPatch's need not have the same length.
Depending on the given images and the stage of the process,
the length of a ScatterPatch might not be as large as we desire.
The algorithm does not require the given tensors (usually images)
to be larger than a minimum size needed to create lengthy patches.
More discussion later, under "testing."

The usual size is defined by parameters.maxPatchSize

Disambiguation
==============

Originally called neighbors (plural) or neighborhood, but that denotes contiguous.
We still use the word neighbor for an element of a scatterPatch,
since neighbor (singular) does not denote adjacent.

Patches vary across and within passes of the algorithm
======================================================

During the Resynthesizer process, when repeating a point,
the point's scatterPatch might have more, different, closer neighbors
than on a prior (especially the first) pass.
IOW a scatterPatch is like a shotgun pattern in the first pass,
but a contiguous patch in later passes.

Also, since we are replacing the value of a point,
that affects the patches of its nearby points.

Distinguished point
===================

A scatterPatch is around a distinguished point.
Often, but not necessarily, the distinguished point is near the center.
An example when the distinguished point is not near the center:
when the synth region is at the edge of the target.

Defined in target, overlaid on corpus
=====================================

A ScatterPatch comprises target points.
The same pattern of a ScatterPatch is searched for in the corpus.
But the same pattern overlain on the corpus might be out of bounds,
i.e. is clipped.

Testing
=======

Case: vector of offsets shorter in length than  ScatterPatch.
Expect ScatterPatch no longer in length than vector of offsets.

Case: tensor not longer in length than ScatterPatch.
(e.g. small image)
Expect ScatterPatch no longer in length than tensor.

Case: tensor small and synth area large.
Because the algorithm starts with the synth area uninitialized,
in this case the context is shorter in length than ScatterPatch.
Expect during the startup, length of ScatterPatch is shorter than desired,
and only later reaches the desired ScatterPatch length.

Case: texture synthesis without "use border"
In this case, then entire target is being synthesized.
So there is never a context, the synth region covers the entire target.
In this case, the first probe for the first point of the first pass
is the best match.
IOW more or less choose a best match at random.
Expect that during startup, length of ScatterPatch grows from 1
to desired length of ScatterPatch
(as each synthesized point is added to the ScatterPatch of the next search.)

TODO review and create test cases for this discussion.
=#

#=
OLD alias for type, and a function that is not a constructor

# Define type
# Original was Vector{Neighbor}
ScatterPatch = Vector{Neighbor{DimensionCount, ValueType} where {DimensionCount, ValueType}}
=#

struct ScatterPatch{ValueType, DimensionCount}
    neighbors::Vector{Neighbor{ValueType, DimensionCount}}
end

#=
Outer constructor.
Constructs fixed length but uninitialized.
From an image having desired type.
=#
function ScatterPatch(image::MaskedImage{ValueType, DimensionCount}) where {ValueType, DimensionCount}
    #=
    Fixed length vector of neighbors, each neighbor is undef, not the reference to the Neighbor.
    =#
    # Fail: undef reference later:
    # neighbors = Vector{Neighbor{ValueType, DimensionCount}}(undef, parameters.maxPatchSize)

    # in vector comprehension, call inner constructor of Neighbor
    # Fail: Cannot `convert` an object of type UndefInitializer to an object of type CartesianIndex{2}
    # neighbors = [Neighbor{ValueType, DimensionCount}(undef, undef, undef) for _ = 1:parameters.maxPatchSize]
    neighbors = [Neighbor(image) for _ = 1:parameters.maxPatchSize]
    # call inner constructor
    return ScatterPatch(neighbors)
end

#=
Return a ScatterPatch
with elements both inside the synth, and outside i.e. in the context (if use_border).

Ensure each neighbor is in the target (the synth or its context)
i.e. not clipped
(and TODO and if the target has another mask, for an irregular shape,
is in the irregular shape. Usually transparency is used instead of a second mask.)

Of a given maximum length.  May be shorter than the given maximum.
The absolute minimum is one (the distinguished point itself), but thats not useful as a patch,
and its not reasonable to have a test case for it.

Neighbors in the ScatterPatch are the nearest distance to the given distinguishedSynthPoint.
(Since SortedOffsets is sorted by distance.)

The distinguishedSynthPoint is strictly in the synth region.
Other neighbors can be, but need not be, in the synth region.
On the first pass, the proportion of neighbors that are also in the synth region builds up.
But for example, the first point of the first pass has no neighbors in the synth regions.
Indeeed, if the context is empty (it can be, the synth region can cover the entire target)
Then there are no neighbors at all.

This is a heuristic: best matches are usually nearer.
If the algorithm did not require neighbors to be near,
it would still find matches,
but the algorithm would not do what the user expect: reconstruct plausibly.

ScatterPatch is ordered by distance but the algorithm does not require it.
This is distinct: not only not required to be near, also not required to be ordered.

When choosing neighbors,
filter so that any neighbor in the synth has a value, i.e. been synthesized already.
(Early on, they have been nulled or blackened.)
Use the current synthResult to determine what points have values.

The context region of the target is not mutated
so the value of (color of) a neighbor in the context is constant.

When parameters.withReplacement,
the synth region of the target is mutated.
But we compute a new ScatterPatch after each point is synthesized
(using the new value of any synth points just mutated.)

The value of neighbors is cached in Neighbor.

In original, called prepare_neighbors()
=#

#=
!!! Weirdness of the original algorithm.
The distinguishedSynthPoint is in the ScatterPatch (offset 0,0)
even though during the first pass it doesn't have a value.
And the comparePointsTargetWithCorpus() tested for and omitted the origin point.
Below, in the for loop, we won't push it again, it is filtered.

push!(result, Neighbor(CartesianIndex(0,0), distinguishedSynthPoint) )
count = 1
=#

#=
outer constructor from a MaskedImage.

This performs poorly, many allocations at push.
=#
#=
NOT USED
function ScatterPatch(
        targetImage::MaskedImage{ValueType, DimensionCount},
        distinguishedSynthPoint::CartesianIndex{DimensionCount},
        synthResult::SynthResult{DimensionCount},
        sortedOffsets::SortedOffsets{DimensionCount}
        ) where {ValueType, DimensionCount}
    # sortedOffsets are offsets to span the tensor of the distinguishedSynthPoint

    # empty vector of undef neighbors
    neighbors = Vector{Neighbor{ValueType, DimensionCount}}(undef, 0)

    count = 0
    for offset in sortedOffsets.offsets
        # wild: may be out of bounds
        wildPoint = pointAtOffset(distinguishedSynthPoint, offset)

        # TODO filtering side effect on wildPoint i.e. toroid
        if ! isPatchPointInTargetFiltered(
            targetImage, wildPoint,     # framed point
            synthResult)
            # assert wildPoint is not wild: is in bounds

            #=
            TODO this allocates.  Because it is a push?
            A better implementation would be to use a fixed length vector
            and use indexing but then return the length of the valid portion of the vector?
            =#
            push!(neighbors,
                Neighbor(
                    offset,
                    targetImage, wildPoint  # framed point
                ) )
            count += 1
            if count >= parameters.maxPatchSize
                break
            end
        end
    end

    # TODO fixed length Vector??
    #result[0] = Neighbor(CartesianIndex(0,0))

    # Except in edge use cases (e.g. synth large and target small), the patch is full size.
    if length(neighbors) != 10
        @debug "Patch small, length %d" length(neighbors)
    end

    @debug "Scatter patch" neighbors
    # call default constructor
    return ScatterPatch(neighbors)
end
=#

#=
Setter method for ScatterPatch.
Given an existing ScatterPatch, mutate its elements for a new set of neighbors.

This avoids allocating an reallocating the ScatterPatch vector
and the Neighbor elements of the vector.
=#
function prepareScatterPatch(
        scatterPatch::ScatterPatch{ValueType, DimensionCount},  # out, mutated
        targetImage::MaskedImage{ValueType, DimensionCount},
        distinguishedSynthPoint::CartesianIndex{DimensionCount},
        synthResult::SynthResult{DimensionCount},
        sortedOffsets::SortedOffsets{DimensionCount}
        )::UInt where {ValueType, DimensionCount}

    #=
    sortedOffsets are offsets to span the tensor of the distinguishedSynthPoint.
    [0,0] i.e. the origin offset is at sortedOffsets[1] !!!
    =#

    neighborIndex = 1

    # OLD do not exclude the distinguished point
    # for offset in sortedOffsets.offsets

    #=
    sortedOffsets might include the zero offset.
    IOW the distinguished point is its own neighbor.
    See sortedOffsets.jl
    =#
    for offset in sortedOffsets.offsets

        wildPoint = offsetPatchPoint(targetImage.image, distinguishedSynthPoint, offset)
        #OLD wildPoint = pointAtOffset(distinguishedSynthPoint, offset)
        # wild: may be out of bounds

        # TODO filtering side effect on wildPoint i.e. toroid
        if ! isPatchPointInTargetFiltered(
                targetImage,
                wildPoint,     # framed point
                synthResult)
            #=
            assert wildPoint is not wild: is in bounds
            assert wildPoint has a value
            =#

            scatterPatch.neighbors[neighborIndex].offset = offset
            scatterPatch.neighbors[neighborIndex].targetPoint = wildPoint  # framed point
            scatterPatch.neighbors[neighborIndex].value = targetImage.image[wildPoint]

            # @debug "Neighbor" neighborIndex scatterPatch.neighbors[neighborIndex]

            neighborIndex += 1
            if neighborIndex >= parameters.maxPatchSize
                # neighborIndex is at last element of vector
                break
            end
        end
    end

    # Except in edge use cases (e.g. synth large and target small), the patch is full size.
    if neighborIndex != parameters.maxPatchSize
        @debug "Patch small, length %d" neighborIndex
    end

    # @debug "Prepared scatter patch" scatterPatch
    # return count of valid points in synthPatch
    return neighborIndex
end

#=
For the Resynthesizer algorithm,
ensure a point in a ScatterPatch:
1) is in bounds of the target (TODO or can be toroidally wrapped to the target)
2) when is in the synth, has a value (synthesized already)
=#
function isPatchPointInTargetFiltered(targetImage, patchPoint, synthResult)::Bool
    if isPointInBoundsOfImage(targetImage, patchPoint)
        if ! doesPointHaveValue(synthResult, patchPoint)
            # in bounds, but is in synth region and has no color yet
            return  true
        else
            # in bounds and has a value
            return false  # not filtered, i.e. valid
        end
    else
        # out of bounds of target

        # TODO try wrap to toroid
        # if we succeed in wrapping, again ensure it has a value
        # return a new point

        return true
    end
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


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

# global vars

CELL_SIZE = 0
ROBOT_CURRENT_POS = [0, 0]

function mark_chess(r::Robot, n::Int = 0)
    global CELL_SIZE

    if n != 0
        CELL_SIZE = n
    end

    side = Right
    mark_row(r, side)
    while isborder(r, Up) == false
        move_decart!(r, Up)
        side = inverse(side)
        mark_row(r,side)
    end

end

function get_cell_size()
    global CELL_SIZE
    println("Введите размер клетки: ")
    CELL_SIZE = parse(Int, readline())
end

function mark_row(r::Robot, side::HorizonSide)
    putmarker_chess!(r)
    while isborder(r, side) == false
        move_decart!(r, side)
        putmarker_chess!(r)
    end
end

function putmarker_chess!(r::Robot)
    global ROBOT_CURRENT_POS
    if (mod(ROBOT_CURRENT_POS[1], 2 * CELL_SIZE)) < CELL_SIZE && (mod(ROBOT_CURRENT_POS[2], 2 * CELL_SIZE)) < CELL_SIZE || (mod(ROBOT_CURRENT_POS[1] + CELL_SIZE, 2 * CELL_SIZE)) < CELL_SIZE && (mod(ROBOT_CURRENT_POS[2], 2 * CELL_SIZE)) >= CELL_SIZE
        putmarker!(r)
    end
end

function move_decart!(r::Robot, side::HorizonSide)
    global ROBOT_CURRENT_POS
    if side == Up
        ROBOT_CURRENT_POS[2] += 1
    elseif side == Down
        ROBOT_CURRENT_POS[2] -= 1
    elseif side == Right
        ROBOT_CURRENT_POS[1] += 1
    else
        ROBOT_CURRENT_POS[1] -= 1
    end
    move!(r, side)
end

ext: jl


text:
function member_attr(node, attr)
    attr == "name" && haskey(node, "alias") && return node["name"] # special handling for aliases
    val = findfirst(".//$attr", node)
    isnothing(val) && error("Attribute $attr not found in node\n$node")
    val.content
end

function resolve_aliases!(collection::Dict, nodes)
    for node ∈ nodes
        if haskey(node, "alias")
            collection[node["name"]] = collection[node["alias"]]
        end
    end
end

function extract_type(param; include_pointer=true)
    type = split(replace(param.content, "const " => ""))[1]
    type_nostar = rstrip(type, ['*', ' '])
    endswith(type, "*") && include_pointer ? "Ptr{$(type_nostar)}" : type_nostar
end

extract_identifier(param) = findfirst("./name", param).content
getattr(node, attr; default = nothing) = haskey(node, attr) ? node[attr] : default

function parent_name(node)
    parel = node.parentelement
    parel.name == "command" ? command_name(parel) : parel.name == "type" && parel["category"] ∈ ["struct", "union"] ? struct_name(parel) : error("Unknown parent element:\n    $parel")
end

function command_name(node)
    isnothing(findfirst("proto", node)) && return command_name(node.parentelement)
    findfirst("proto/name", node).content
end

function struct_name(node)
    (!haskey(node, "category") || node["category"] ∉ ["struct", "union"]) && return struct_name(node.parentelement)
    node["name"]
end

is_constant(node) = any(split(node.content) .== "const")

equals(column, value) = x -> getproperty(x, column) == value

function dfmatch(f, df)
    index = findfirst(f, eachrow(df))
    isnothing(index) && error("No match")
    df[index, :]
end

dfmatch(df, column, value) = dfmatch(equals(column, value), df)
function dfmatches(f, df)
    indices = findall(f, eachrow(df))
    isnothing(indices) && error("No match")
    getindex.(Ref(df), indices, :)
end

dfmatches(df, column, value) = dfmatches(equals(column, value), df)


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCell
include("useful_routines.jl")

function reinitialization_error(distancefunction, nelmts, polyorder)
    L, W = 1.0, 1.0
    basis = TensorProductBasis(2, polyorder)
    numqp = required_quadrature_order(polyorder) + 4
    quad = tensor_product_quadrature(2, numqp)
    levelset = InterpolatingPolynomial(1, basis)

    mesh = CutCell.Mesh([0.0, 0.0], [L, W], [nelmts, nelmts], basis)
    levelsetcoeffs = CutCell.levelset_coefficients(distancefunction, mesh)
    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)

    refseedpoints, spatialseedpoints, seedcellids =
        CutCell.seed_zero_levelset(2, levelset, levelsetcoeffs, cutmesh)

    nodalcoordinates = CutCell.nodal_coordinates(cutmesh)
    signeddistance = CutCell.distance_to_zero_levelset(
        nodalcoordinates,
        refseedpoints,
        spatialseedpoints,
        seedcellids,
        levelset,
        levelsetcoeffs,
        cutmesh,
        1e-8,
        4.5,
    )

    err = uniform_mesh_L2_error(
        signeddistance',
        x -> distancefunction(x)[1],
        basis,
        quad,
        cutmesh.mesh,
    )
    cellmaps = CutCell.cell_maps(cutmesh.mesh)
    den = integral_norm_on_uniform_mesh(distancefunction, quad, cellmaps, 1)

    return err[1] / den[1]
end

function convergence_rate(v, dx)
    return diff(log.(v)) ./ diff(log.(dx))
end

xc = [0.5, 0.5]
rad = 0.25
polyorder = 2
powers = [2, 3, 4, 5]
nelmts = [2^i + 1 for i in powers]

dx = 1.0 ./ nelmts
err = [
    reinitialization_error(x -> circle_distance_function(x, xc, rad), ne, polyorder)
    for ne in nelmts
]
rate = convergence_rate(err, dx)
@test allapprox(rate, 2 * ones(length(rate)), 0.05)

corner = [0.5, 0.5]
err = [
    reinitialization_error(x -> corner_distance_function(x, corner), ne, polyorder)
    for ne in nelmts
]
rate = convergence_rate(err,dx)
@test all(rate .> 1.0)


text:
@testset "LU solver" begin
    n = 20
    m = 10
    A = rand(n, n)
    X = rand(n, m)
    B = rand(n, m)
    x = rand(n)
    b = rand(n)

    solver = RoboDojo.empty_solver(A)
    @test solver isa RoboDojo.LinearSolver

    solver = RoboDojo.lu_solver(A)
    RoboDojo.linear_solve!(solver, X, A, B)
    @test norm(A * X - B, Inf) < 1.0e-10

    solver = RoboDojo.lu_solver(A)
    RoboDojo.linear_solve!(solver, x, A, b)
    @test norm(A * x - b, Inf) < 1.0e-10

end


text:
using InelasticScattering
using Plots
λ = 440*1.e-7 #cm
ν̃ = 1/λ
T = 300 #k

n2 = InelasticScattering.getMolecularConstants(InelasticScattering.N₂(), (0.8));
compute_effective_coefficents!(ν̃, T, n2)
compute_energy_levels!(n2)
compute_σ_Rayl_coeff!(n2)
compute_σ_Rayl_VibRaman_coeff_hires!(T, n2)
compute_σ_VibRaman_coeff!(T, n2)
compute_σ_RoVibRaman_coeff!(T, n2)

o2 = InelasticScattering.getMolecularConstants(InelasticScattering.O₂(), (0.2));
compute_effective_coefficents!(ν̃, T, o2)
compute_energy_levels!(o2)
compute_σ_Rayl_coeff!(o2)
compute_σ_Rayl_VibRaman_coeff_hires!(T, o2)
compute_σ_VibRaman_coeff!(T, o2)
compute_σ_RoVibRaman_coeff!(T, o2)

plotLegend = false
plot_coefficients = 1
plot_cross_sections = 0
if plot_coefficients==1
    @show "YEAH"
    # N2
    # Rayleigh scattering, pure vibrational Raman scattering
    x = [n2.effCoeff.Δν̃_VibRaman_coeff_0to1 n2.effCoeff.Δν̃_VibRaman_coeff_1to0 0];
    y = [n2.effCoeff.σ_VibRaman_coeff_0to1 n2.effCoeff.σ_VibRaman_coeff_1to0 n2.effCoeff.σ_Rayl_coeff];
    @show(ν̃)
    @show(x)
    @show(y)
    #plot(x, y*1e50, color=:blue, seriestype=:sticks, legend=plotLegend)

    # Rayleigh hires scattering
    #x = [n2.effCoeff.Δν̃_Rayl_coeff_hires.-n2.effCoeff.Δν̃_Rayl_coeff_hires[0]];
    #y = [n2.effCoeff.σ_Rayl_coeff_hires];
    #@show(ν̃)
    #@show(x)
    #@show(y)
    #plot(x, y*1e50, color=:blue, seriestype=:sticks, legend=plotLegend)

    # Rotational Raman scattering
    x = [n2.effCoeff.Δν̃_RoRaman_coeff_JtoJp2 n2.effCoeff.Δν̃_RoRaman_coeff_JtoJm2];
    y = [n2.effCoeff.σ_RoRaman_coeff_JtoJp2 n2.effCoeff.σ_RoRaman_coeff_JtoJm2];
    @show(x)
    @show(y)
    #plot!(x, y * 1e50, color=:blue, seriestype=:sticks, legend=plotLegend)

    #Hires vibrational scattering
    x = [n2.effCoeff.Δν̃_VibRaman_coeff_0to1_hires]# n2.effCoeff.Δν̃_VibRaman_coeff_1to0_hires];
    y = [n2.effCoeff.σ_VibRaman_coeff_0to1_hires]# n2.effCoeff.σ_VibRaman_coeff_1to0_hires];
    @show(x)
    @show(y)
    plot(x, y*1e50, color=:blue, seriestype=:sticks, legend=plotLegend)

    # Rovibrational scattering (P)
    x = [n2.effCoeff.Δν̃_RoVibRaman_coeff_1to0_JtoJp2 n2.effCoeff.Δν̃_RoVibRaman_coeff_1to0_JtoJm2];
    y = [n2.effCoeff.σ_RoVibRaman_coeff_1to0_JtoJp2 n2.effCoeff.σ_RoVibRaman_coeff_1to0_JtoJm2];
    @show(x)
    @show(y)
    #plot!(x, y*1e50, color=:blue, seriestype=:sticks, legend=plotLegend)

    # Rovibrational scattering (R)
    x = [n2.effCoeff.Δν̃_RoVibRaman_coeff_0to1_JtoJp2 n2.effCoeff.Δν̃_RoVibRaman_coeff_0to1_JtoJm2];
    y = [n2.effCoeff.σ_RoVibRaman_coeff_0to1_JtoJp2 n2.effCoeff.σ_RoVibRaman_coeff_0to1_JtoJm2];
    @show(x)
    @show(y)
    plot!(x, y*1e50, color=:blue, seriestype=:sticks, legend=plotLegend)

    # O2
    # Rayleigh scattering, pure vibrational Raman scattering
    x = [o2.effCoeff.Δν̃_VibRaman_coeff_0to1 o2.effCoeff.Δν̃_VibRaman_coeff_1to0 0];
    y = [o2.effCoeff.σ_VibRaman_coeff_0to1 o2.effCoeff.σ_VibRaman_coeff_1to0 o2.effCoeff.σ_Rayl_coeff];
    @show(ν̃)
    @show(x)
    @show(y)
    #plot!(x, y*1e50, color=:orange, seriestype=:sticks, legend=plotLegend)

    # Rotational Raman scattering
    x = [o2.effCoeff.Δν̃_RoRaman_coeff_JtoJp2 o2.effCoeff.Δν̃_RoRaman_coeff_JtoJm2];
    y = [o2.effCoeff.σ_RoRaman_coeff_JtoJp2 o2.effCoeff.σ_RoRaman_coeff_JtoJm2];
    @show(x)
    @show(y)
    #plot!(x, y*1e50, color=:orange, seriestype=:sticks, legend=plotLegend)

    #Hires vibrational scattering
    x = [o2.effCoeff.Δν̃_VibRaman_coeff_0to1_hires]# o2.effCoeff.Δν̃_VibRaman_coeff_1to0_hires];
    y = [o2.effCoeff.σ_VibRaman_coeff_0to1_hires]# o2.effCoeff.σ_VibRaman_coeff_1to0_hires];
    @show(x)
    @show(y)
    plot!(x, y*1e50, color=:orange, seriestype=:sticks, legend=plotLegend)

    # Rovibrational scattering (P)
    x = [o2.effCoeff.Δν̃_RoVibRaman_coeff_1to0_JtoJp2 o2.effCoeff.Δν̃_RoVibRaman_coeff_1to0_JtoJm2];
    y = [o2.effCoeff.σ_RoVibRaman_coeff_1to0_JtoJp2 o2.effCoeff.σ_RoVibRaman_coeff_1to0_JtoJm2];
    @show(x)
    @show(y)
    #plot!(x, y*1e50, color=:orange, seriestype=:sticks, legend=plotLegend)

    # Rovibrational scattering (R)
    x = [o2.effCoeff.Δν̃_RoVibRaman_coeff_0to1_JtoJp2 o2.effCoeff.Δν̃_RoVibRaman_coeff_0to1_JtoJm2];
    y = [o2.effCoeff.σ_RoVibRaman_coeff_0to1_JtoJp2 o2.effCoeff.σ_RoVibRaman_coeff_0to1_JtoJm2];
    @show(x)
    @show(y)
    plot!(x, y*1e50, color=:orange, seriestype=:sticks, legend=plotLegend)

elseif plot_cross_sections==1
    # N2
    # Rayleigh scattering, pure vibrational Raman scattering
    x = [ν̃+n2.effCoeff.Δν̃_VibRaman_coeff_0to1 ν̃+n2.effCoeff.Δν̃_VibRaman_coeff_1to0 ν̃];
    y = [n2.effCoeff.σ_VibRaman_coeff_0to1*(ν̃+n2.effCoeff.Δν̃_VibRaman_coeff_0to1)^4 n2.effCoeff.σ_VibRaman_coeff_1to0*(ν̃+n2.effCoeff.Δν̃_VibRaman_coeff_1to0)^4 n2.effCoeff.σ_Rayl_coeff*(ν̃)^4];
    @show(ν̃)
    @show(x)
    @show(y)
    #plot((1e7) ./x, y*1e20, color=:blue, seriestype=:sticks, legend=plotLegend)

    # Rayleigh hires scattering
    x = [n2.effCoeff.Δν̃_Rayl_coeff_hires];
    y = [n2.effCoeff.σ_Rayl_coeff_hires];
    @show(ν̃)
    @show(x)
    @show(y)
    plot((1e7) ./x, y*1e50, color=:blue, seriestype=:sticks, legend=plotLegend)

    # Rotational Raman scattering
    x = [ν̃.+n2.effCoeff.Δν̃_RoRaman_coeff_JtoJp2' ν̃.+n2.effCoeff.Δν̃_RoRaman_coeff_JtoJm2'];
    y = [n2.effCoeff.σ_RoRaman_coeff_JtoJp2'.*(ν̃.+n2.effCoeff.Δν̃_RoRaman_coeff_JtoJp2').^4 n2.effCoeff.σ_RoRaman_coeff_JtoJm2'.*(ν̃.+n2.effCoeff.Δν̃_RoRaman_coeff_JtoJm2').^4];
    @show(x)
    @show(y)
    plot!((1e7)*1 ./x, y*1e20, color=:blue, seriestype=:sticks, legend=plotLegend)

    # Rovibrational scattering (P)
    x = [ν̃.+n2.effCoeff.Δν̃_RoVibRaman_coeff_1to0_JtoJp2' ν̃.+n2.effCoeff.Δν̃_RoVibRaman_coeff_1to0_JtoJm2'];
    y = [n2.effCoeff.σ_RoVibRaman_coeff_1to0_JtoJp2'.*(ν̃.+n2.effCoeff.Δν̃_RoVibRaman_coeff_1to0_JtoJp2').^4 n2.effCoeff.σ_RoVibRaman_coeff_1to0_JtoJm2'.*(ν̃.+n2.effCoeff.Δν̃_RoVibRaman_coeff_1to0_JtoJm2').^4];
    @show(x)
    @show(y)
    #plot!((1e7)*1 ./x, y*1e20, color=:blue, seriestype=:sticks, legend=plotLegend)

    # Rovibrational scattering (R)
    x = [ν̃.+n2.effCoeff.Δν̃_RoVibRaman_coeff_0to1_JtoJp2' ν̃.+n2.effCoeff.Δν̃_RoVibRaman_coeff_0to1_JtoJm2'];
    y = [n2.effCoeff.σ_RoVibRaman_coeff_0to1_JtoJp2'.*(ν̃.+n2.effCoeff.Δν̃_RoVibRaman_coeff_0to1_JtoJp2').^4 n2.effCoeff.σ_RoVibRaman_coeff_0to1_JtoJm2'.*(ν̃.+n2.effCoeff.Δν̃_RoVibRaman_coeff_0to1_JtoJm2').^4];
    @show(x)
    @show(y)
    #plot!((1e7)*1 ./x, y*1e20, color=:blue, seriestype=:sticks, legend=plotLegend)

    # O2
    # Rayleigh scattering, pure vibrational Raman scattering
    x = [ν̃+o2.effCoeff.Δν̃_VibRaman_coeff_0to1 ν̃+o2.effCoeff.Δν̃_VibRaman_coeff_1to0 ν̃];
    y = [o2.effCoeff.σ_VibRaman_coeff_0to1*(ν̃+o2.effCoeff.Δν̃_VibRaman_coeff_0to1)^4 o2.effCoeff.σ_VibRaman_coeff_1to0*(ν̃+o2.effCoeff.Δν̃_VibRaman_coeff_1to0)^4 o2.effCoeff.σ_Rayl_coeff*(ν̃)^4];
    #@show(ν̃)
    #@show(x)
    #@show(y)
    #plot!((1e7)*1 ./x, y*1e20, color=:orange, seriestype=:sticks, legend=plotLegend)

    # Rayleigh hires scattering
    x = [o2.effCoeff.Δν̃_Rayl_coeff_hires];
    y = [o2.effCoeff.σ_Rayl_coeff_hires];
    #@show(ν̃)
    #@show(x)
    #@show(y)
    #plot((1e7) ./x, y*1e50, color=:blue, seriestype=:sticks, legend=plotLegend)

    # Rotational Raman scattering
    x = [ν̃.+o2.effCoeff.Δν̃_RoRaman_coeff_JtoJp2' ν̃.+o2.effCoeff.Δν̃_RoRaman_coeff_JtoJm2'];
    y = [o2.effCoeff.σ_RoRaman_coeff_JtoJp2'.*(ν̃.+o2.effCoeff.Δν̃_RoRaman_coeff_JtoJp2').^4 o2.effCoeff.σ_RoRaman_coeff_JtoJm2'.*(ν̃.+o2.effCoeff.Δν̃_RoRaman_coeff_JtoJm2').^4];
    #@show(x)
    #@show(y)
    #plot!((1e7)*1 ./x, y*1e20, color=:orange, seriestype=:sticks, legend=plotLegend)

    #Hires vibrational scattering
    x = ν̃.+o2.effCoeff.Δν̃_VibRaman_coeff_0to1_hires# ν̃+o2.effCoeff.Δν̃_VibRaman_coeff_1to0_hires];
    y = o2.effCoeff.σ_VibRaman_coeff_0to1_hires.*(ν̃.+o2.effCoeff.Δν̃_VibRaman_coeff_0to1_hires).^4# o2.effCoeff.σ_VibRaman_coeff_1to0_hires];
    @show(x)
    @show(y)
    plot(1e7 ./x, y*1e20, color=:orange, seriestype=:sticks, legend=plotLegend)

    # Rovibrational scattering (P)
    x = [ν̃.+o2.effCoeff.Δν̃_RoVibRaman_coeff_1to0_JtoJp2' ν̃.+o2.effCoeff.Δν̃_RoVibRaman_coeff_1to0_JtoJm2'];
    y = [o2.effCoeff.σ_RoVibRaman_coeff_1to0_JtoJp2'.*(ν̃.+o2.effCoeff.Δν̃_RoVibRaman_coeff_1to0_JtoJp2').^4 o2.effCoeff.σ_RoVibRaman_coeff_1to0_JtoJm2'.*(ν̃.+o2.effCoeff.Δν̃_RoVibRaman_coeff_1to0_JtoJm2').^4];
    @show(x)
    @show(y)
    #plot!((1e7)*1 ./x, y*1e20, color=:orange, seriestype=:sticks, legend=plotLegend)

    # Rovibrational scattering (R)
    x = [ν̃.+o2.effCoeff.Δν̃_RoVibRaman_coeff_0to1_JtoJp2' ν̃.+o2.effCoeff.Δν̃_RoVibRaman_coeff_0to1_JtoJm2'];
    y = [o2.effCoeff.σ_RoVibRaman_coeff_0to1_JtoJp2'.*(ν̃.+o2.effCoeff.Δν̃_RoVibRaman_coeff_0to1_JtoJp2').^4 o2.effCoeff.σ_RoVibRaman_coeff_0to1_JtoJm2'.*(ν̃.+o2.effCoeff.Δν̃_RoVibRaman_coeff_0to1_JtoJm2').^4];
    @show(x)
    @show(y)
    plot!((1e7) ./x, y*1e20, color=:orange, seriestype=:sticks, legend=plotLegend)
end
#foreground_color_legend = nothing
#background_color_legend = nothing

max_stars_repo_name: RadiativeTransfer/InelasticScattering

max_issues_repo_name: RadiativeTransfer/InelasticScattering

max_forks_repo_name: RadiativeTransfer/InelasticScattering


text:
"""
    CompatDevTools

Featured functions:
 - [`compat_kick_start`](@ref) kick-start compat entries from a directory
"""
module CompatDevTools

import OrderedCollections
import REPL.TerminalMenus

if VERSION >= v"1.6.0"
    radio_menu(options; kwargs...) = TerminalMenus.RadioMenu(options; charset=:ascii, kwargs...)
else
    radio_menu(options; kwargs...) = TerminalMenus.RadioMenu(options; kwargs...)
end

is_deps(x) = startswith(x, "[[deps.") || startswith(x, "[deps]")
dep_match(x, dep) = x == "[[$dep]]" || x ==  "[[deps.$dep]]"

"""
    synchronize_compats(code_dir::AbstractString)

This function
 - Recursively finds all Project.toml files in `code_dir`
 - Collects the compat entries
 - Finds any inconsistent compat entries
 - Asks the user (via `REPL.TerminalMenus`) which
   version (if any) to update to, and modifies the
   Project.toml files accordingly.
"""
function synchronize_compats(code_dir::AbstractString)

    project_toml_files = [joinpath(root, f) for (root, dirs, files) in Base.Filesystem.walkdir(code_dir) for f in files if endswith(f, "Project.toml")]

    compat_entries = OrderedCollections.OrderedDict()
    for project_toml in project_toml_files
        contents = join(readlines(project_toml; keep=true))
        if !occursin("[compat]", contents)
            @warn "No compat entry found in $project_toml"
            continue
        end
        if !occursin("[compat]\n", contents)
            @warn "Compat header found, but no compat entry found in $project_toml"
            continue
        end
        compat_section = last(split(contents, "[compat]\n"))
        if occursin("\n\n", compat_section) # split by end of compat section
            compat_section = first(split(compat_section, "\n\n"))
        end
        compats = split(compat_section, "\n")
        filter!(x->!isempty(x), compats)
        compats = map(compats) do compat
            s = split(compat, "=")
            (strip(first(s)), replace(strip(last(s)), "\"" => ""), compat)
        end
        dir = dirname(project_toml)
        for (pkg, ver, line) in compats
            if haskey(compat_entries, pkg)
                ver in compat_entries[pkg] && continue
                push!(compat_entries[pkg], (ver, dir, line))
            else
                compat_entries[pkg] = [(ver, dir, line)]
            end
        end
    end

    inconsistent_compat_entries = OrderedCollections.OrderedDict()
    ice = inconsistent_compat_entries
    for k in keys(compat_entries)
        length(compat_entries[k]) == 1 && continue # only one env with this dep
        if length(unique(first.(compat_entries[k]))) ≠ 1
            ice[k] = compat_entries[k]
        end
    end
    @debug begin
        for k in keys(ice)
            println("$k = $(ice[k])")
        end
    end
    if length(keys(ice)) == 0
        @info "All compat entries are consistent! 🎉"
        return nothing
    end
    answers = OrderedCollections.OrderedDict()
    @info "$(length(keys(ice))) inconsistent compat entries found across $(length(project_toml_files)) toml files."
    for pkg in keys(ice)
        versions = getindex.(ice[pkg], 1)
        env_dirs = getindex.(ice[pkg], 2)
        lines = getindex.(ice[pkg], 3)
        entries = map(zip(versions, env_dirs)) do (ver, env_dir)
            "$ver ($(env_dir))"
        end
        options = [entries..., "Leave alone."]
        menu = radio_menu(options)
        msg = "Inconsistent compat entries found for $pkg. Select a compat entry (or leave alone):"
        choice = TerminalMenus.request(msg, menu)
        if !(choice == length(options))
            answers[pkg] = (env_dirs[choice], lines[choice], lines)
            @info "Changing compat entry for $pkg to $(lines[choice]) for all environments ($(env_dirs))"
        end
    end

    # Edit Project.toml files according to answers

    for project_toml in project_toml_files
        contents = join(readlines(project_toml; keep=true))
        new_contents = contents
        for pkg in keys(answers)
            all_compat_lines = last(answers[pkg])
            chosen_compat_line = getindex(answers[pkg], 2)
            for compat_line in all_compat_lines
                new_contents = replace(new_contents, compat_line => chosen_compat_line)
            end
        end
        open(project_toml, "w") do io
            print(io, new_contents)
        end
    end

end

"""
    compat_kick_start(code_dir::String, julia_version = "1.5.4")

Given a directory containing both a Project.toml and Manifest.toml,
print a string of a suggested compat entry for all Project.toml
dependencies.
"""
function compat_kick_start(code_dir::AbstractString, julia_version = "1")
    project = joinpath(code_dir, "Project.toml");
    manifest = joinpath(code_dir, "Manifest.toml");
    @assert isfile(project)
    @assert isfile(manifest)

    project_contents = readlines(project)
    deps = String[]
    collect_deps = false
    end_collection = false
    for line in project_contents
        end_collection && break
        if is_deps(line)
            collect_deps = true
        end
        if collect_deps
            if line == ""
                end_collection = true
            end
        end
        if collect_deps && !is_deps(line) && line ≠ ""
            push!(deps, first(split(line, " = ")))
        end
    end
    compat_entries = String[]
    manifest_contents = readlines(manifest)
    for dep in deps
        i = findfirst(x -> dep_match(x, dep), manifest_contents)
        i == nothing && continue
        j = findfirst(x -> x == "", manifest_contents[i:end])

        contents_block = manifest_contents[i:i+j-1]
        i_version_entry = findfirst(x->startswith(x, "version = "), contents_block)
        if i_version_entry == nothing
            continue # must be a stdlib
        end
        version_str = contents_block[i_version_entry]
        version_num = last(split(version_str, "version = "))

        push!(compat_entries, "$dep = $version_num")
    end
    push!(compat_entries, "julia = \"$julia_version\"")

    str = "[compat]\n"*join(compat_entries, "\n")
    print(str)
    return str
end

end # module


text:
@doc doc"""
    latexraw(arg)

Generate LaTeX equations from `arg`.

Parses expressions, ParameterizedFunctions, SymEngine.Base and arrays thereof.
Returns a string formatted for LaTeX.

# Examples

## using expressions
```jldoctest
expr = :(x/(y+x))
latexraw(expr)

# output

"\\frac{x}{y + x}"
```

```jldoctest
expr = Meta.parse("x/(y+x)")
latexraw(expr)

# output

"\\frac{x}{y + x}"
```

## using ParameterizedFunctions
```julia
using DifferentialEquations;
f = @ode_def feedback begin
         dx = y/c_1 - x
         dy = x^c_2 - y
       end c_1=>1.0 c_2=>1.0
latexraw(f)

# output

2-element Array{String,1}:
 "dx/dt = \\frac{y}{c_{1}} - x"
 "dy/dt = x^{c_{2}} - y"
```

## using SymEngine
```jldoctest
using SymEngine
@vars x y
symExpr = x + x + x*y*y
latexraw(symExpr)

# output

"2 \\cdot x + x \\cdot y^{2}"
```
"""
latexraw(args...; kwargs...) = process_latexify(args...; kwargs..., env=:raw)

function _latexraw(inputex::Expr; convert_unicode=true, kwargs...)
    ## Pass all arrays or matrices in the expr to latexarray
    inputex = postwalk(x -> x isa Expr && x.head in [:hcat, :vcat, :vect, :typed_vcat, :typed_hcat] ?
                       latexarray(expr_to_array(x); kwargs...)
                       : x,
                       inputex)

    recurseexp!(lstr::LaTeXString) = lstr.s
    function recurseexp!(ex)
        prevOp = Vector{Symbol}(undef, length(ex.args))
        fill!(prevOp, :none)
        for i in 1:length(ex.args)
            if isa(ex.args[i], Expr)
                length(ex.args[i].args) > 1 && ex.args[i].args[1] isa Symbol && (prevOp[i] = ex.args[i].args[1])
                ex.args[i] = recurseexp!(ex.args[i])
            elseif ex.args[i] isa AbstractArray
                ex.args[i] = latexraw(ex.args[i]; kwargs...)
            end
        end
        return latexoperation(ex, prevOp; convert_unicode=convert_unicode, kwargs...)
    end
    ex = deepcopy(inputex)
    str = recurseexp!(ex)
    convert_unicode && (str = unicode2latex(str))
    return LaTeXString(str)
end

function _latexraw(args...; kwargs...)
    @assert length(args) > 1 "latexify does not support objects of type $(typeof(args[1]))."
    _latexraw(args; kwargs...)
end
_latexraw(arr::Union{AbstractArray, Tuple}; kwargs...) = _latexarray(arr; kwargs...)
_latexraw(i::Nothing; kwargs...) = ""
_latexraw(i::SubString; kwargs...) = latexraw(Meta.parse(i); kwargs...)
_latexraw(i::SubString{LaTeXStrings.LaTeXString}; kwargs...) = i
_latexraw(i::Rational; kwargs...) = i.den == 1 ? latexraw(i.num; kwargs...) : latexraw(:($(i.num)/$(i.den)); kwargs...)
_latexraw(z::Complex; kwargs...) = LaTeXString("$(latexraw(z.re;kwargs...))$(z.im < 0 ? "-" : "+" )$(latexraw(abs(z.im);kwargs...))\\textit{i}")
#latexraw(i::DataFrames.DataArrays.NAtype) = "\\textrm{NA}"
_latexraw(str::LaTeXStrings.LaTeXString; kwargs...) = str

function _latexraw(i::Number; fmt=PlainNumberFormatter(), kwargs...)
    try i == Inf && return LaTeXString("\\infty") catch; end
    try i == -Inf && return LaTeXString("-\\infty") catch; end
    fmt isa String && (fmt = PrintfNumberFormatter(fmt))
    return fmt(i)
end

function _latexraw(i::Char; convert_unicode=true, kwargs...)
    LaTeXString(convert_unicode ? unicode2latex(string(i)) : string(i))
end

function _latexraw(i::Symbol; convert_unicode=true, kwargs...)
    str = string(i == :Inf ? :∞ : i)
    str = convertSubscript(str)
    convert_unicode && (str = unicode2latex(str))
    return LaTeXString(str)
end

function _latexraw(i::String; kwargs...)
    try
        ex = Meta.parse(i)
        return latexraw(ex; kwargs...)
    catch ParseError
        error("""
in Latexify.jl:
You are trying to create latex-maths from a `String` that cannot be parsed as
an expression.

`latexify` will, by default, try to parse any string inputs into expressions
and this parsing has just failed.

If you are passing strings that you want returned verbatim as part of your input,
try making them `LaTeXString`s first.

If you are trying to make a table with plain text, try passing the keyword
argument `latex=false`. You should also ensure that you have chosen an output
environment that is capable of displaying not-maths objects. Try for example
`env=:table` for a latex table or `env=:mdtable` for a markdown table.
""")
    end
end

_latexraw(i::Missing; kwargs...) = "\\textrm{NA}"


text:
macro prop(x,y)
	quote
		if is_CoTangent($x)
			($x).delta+=$y
		end
	end
end

function backprop(x::Variable) end
function backprop(x::Constant) end
function backprop(c::Call)
	backprop_dispatch(c,c.args...)
end
function backprop_dispatch{F}(c::Call{F},args...)
	#println("Warning: no backprop definition for $F on $(map(typeof,args))")
end
function backprop(x::Lift)
	for (a,b) in zip(x.args,x.delta)
		@prop(a,b)
	end
end

function backprop_dispatch(c::Call{:+},args...)
	map(x->@prop(x,c.delta),args)
end

function backprop_dispatch(c::Call{:*},A,B)
	@prop(A,c.delta*transpose(val(B)))
	@prop(B,transpose(val(A))*c.delta)
end

function backprop_dispatch(c::Call{:-},x)
	@prop(x,-c.delta)
end
function backprop_dispatch(c::Call{:-},x,y)
	@prop(x,c.delta)
	@prop(y,-c.delta)
end

function backprop_dispatch(c::Call{:/},A,B)
	@prop(A,c.delta/val(B))
	@prop(B,-c.delta.*val(A)/val(B)^2)
end

function backprop_dispatch(c::Call{:.^},A,n::Real)
	@prop(A,c.delta.*n.*val(A).^(n-1))
end

function backprop_dispatch(c::Call{:sigmoid},A)
	@prop(A, c.delta .* val(c) .* (1.0-val(c)))
end

function backprop_dispatch(c::Call{:rect_lin},A)
	@prop(A,!map(signbit,val(A)).*c.delta)
end

function backprop_dispatch(c::Call{:exp},A)
	@prop(A,c.delta .* val(c))
end

function backprop_dispatch(c::Call{:log},A)
	@prop(A,c.delta./val(A))
end

macro downcast_prop(x,y)
	quote
		if is_CoTangent($x)
			downcast_prop($x,$y)
		end
	end
end

function downcast_prop(c::CoTangent,x)
	#reduction_dims=filter(k->size(c.delta,k)==1,1:length(size(x)))
	reduction_dims=[]
	for i in 1:length(size(x))
		if size(c.delta,i)==1
			push!(reduction_dims,i)
		end
	end
	c.delta+=reshape(sum(x,reduction_dims),size(c.delta))
end
function downcast_prop{T<:Real}(c::CoTangent{T},x)
	c.delta+=sum(x)
end

function backprop_dispatch(c::Call{:.*},A,B)
	@downcast_prop(A,c.delta .* val(B))
	@downcast_prop(B,c.delta .* val(A))
end
function backprop_dispatch(c::Call{:.-},A,B)
	@downcast_prop(A,c.delta)
	@downcast_prop(B,-c.delta)
end
function backprop_dispatch(c::Call{:.-},A)
	@downcast_prop(A,-c.delta)
end
function backprop_dispatch(c::Call{:.+},args...)
	map(x->@downcast_prop(x,c.delta),args)
end

function backprop_dispatch(c::Call{:./},A,B)
	@downcast_prop(A,c.delta./val(B))
	@downcast_prop(B,-c.delta.*val(A)./val(B).^2)
end

macro upcast_prop(x,y)
	quote
		if is_CoTangent($x)
			($x).delta .+= $y
		end
	end
end
function backprop_dispatch(c::Call{:sum},A)
	@upcast_prop(A,Any[c.delta])
end
function backprop_dispatch(c::Call{:sum},A,dims)
	@upcast_prop(A,c.delta)
end

function backprop_dispatch(c::Call{:getindex},A,indices...)
	A.delta[indices...]+=c.delta
end
function backprop_dispatch(c::Call{:permutedims},A,perm)
	@prop(A,ipermutedims(c.delta,perm))
end


text:
using Documenter, Symbolics, SymbolicUtils

makedocs(
    sitename="Symbolics.jl",
    authors="Chris Rackauckas",
    modules=[Symbolics,SymbolicUtils],
    clean=true,doctest=false,
    format = Documenter.HTML(#analytics = "UA-90474609-3",
                             assets = ["assets/favicon.ico"]),
                             #canonical="https://mtk.sciml.ai/stable/"),
    pages=[
        "Home" => "index.md",
        "Tutorials" => Any[
            "tutorials/symbolic_functions.md",
            "tutorials/auto_parallel.md",
            "tutorials/converting_to_C.md"
        ],
        "Manual" => Any[
            "manual/variables.md",
            "manual/expression_manipulation.md",
            "manual/derivatives.md",
            "manual/arrays.md",
            "manual/build_function.md",
            "manual/functions.md",
            "manual/io.md",
            "manual/sparsity_detection.md",
            "manual/faq.md"
        ],
        "Comparison Against SymPy" => "comparison.md",
    ]
)

deploydocs(
   repo = "github.com/JuliaSymbolics/Symbolics.jl.git";
   push_preview = true
)


text:
    """
        mid_point(M, p, q, x)

    Compute the mid point between p and q. If there is more than one mid point
    of (not neccessarily miniizing) geodesics (i.e. on the sphere), the one nearest
    to x is returned.
    """
    mid_point(M::Manifold, p, q, x) = mid_point(M, p, q)
    mid_point!(M::Manifold, y, p, q, x) = mid_point!(M, y, p, q)

    function mid_point(M::Sphere, p, q, x)
        if isapprox(M,p,-q)
            X = log(M,p,x)/distance(M,p,x)*π
        else
            X = log(M,p,q)
        end
        return exp(M,p,0.5*X)
    end
    function mid_point!(M::Sphere, y, p, q, x)
        if isapprox(M,p,-q)
            X = log(M,p,x)/distance(M,p,x)*π
        else
            X = log(M,p,q)
        end
        y .= exp(M,p,0.5*X)
        return y
    end

    """
        mid_point(M, p, q)

    Compute the (geodesic) mid point of the two points `p` and `q` on the
    manfold `M`. If the geodesic is not unique, either a deterministic choice is taken or
    an error is raised depending on the manifold. For the deteministic choixe, see
    [`mid_point(M, p, q, x)`](@ref), the mid point closest to a third point
    `x`.
    """
    mid_point(M::Manifold, p, q) = exp(M, p, log(M, p, q), 0.5)
    mid_point!(M::Manifold, y, p, q) = exp!(M, y, p, log(M, p, q), 0.5)

    @doc raw"""
        reflect(M, f, x)

    reflect the point `x` from the manifold `M` at the point `f(x)` of the
    function $f\colon \mathcal M \to \mathcal M$, i.e.,

    ````math
        \operatorname{refl}_f(x) = \operatorname{refl}_{f(x)}(x),
    ````
    see also [`reflect`](@ref reflect(M::Manifold, p, x))`(M,p,x)`.
    """

    reflect(M::Manifold, pr::Function, x) = reflect(M::Manifold, pr(x), x)

    @doc raw"""
        reflect(M, p, x)

    reflect the point `x` from the manifold `M` at point `x`, i.e.

    ````math
        \operatorname{refl}_p(x) = \exp_p(-\log_p x).
    ````
    where exp and log denote the exponential and logarithmic map on `M`.
    """
    reflect(M::Manifold, p, x) = exp(M, p, -log(M, p, x))

    @doc raw"""
        sym_rem(x,[T=π])

    Compute symmetric remainder of `x` with respect to the interall 2*`T`, i.e.
    `(x+T)%2T`, where the default for `T` is $π$
    """
    function sym_rem(x::N, T = π) where {N<:Number}
        return (x ≈ T ? convert(N, -T) : rem(x, convert(N, 2 * T), RoundNearest))
    end


text:
module CachedSubVectors

using Base: @propagate_inbounds

export CachedSubVector
export locate!
import Base: size
import Base: getindex
import Base: setindex!
import Base: IndexStyle

mutable struct CachedSubVector{T,V<:AbstractArray{T,1}} <: AbstractArray{T,1}
  vector::V
  pini::Int
  pend::Int
end

function CachedSubVector(v::AbstractVector)
  CachedSubVector(v,1,length(v))
end

size(self::CachedSubVector) = (1+self.pend-self.pini,)

@propagate_inbounds function getindex(self::CachedSubVector, i::Int)
  @inbounds self.vector[self.pini+i-1]
end

@propagate_inbounds function setindex!(self::CachedSubVector,v, i::Int)
  @inbounds self.vector[self.pini+i-1] = v
end

IndexStyle(::Type{CachedSubVector{T,V}}) where {T,V} = IndexLinear()

function locate!(self::CachedSubVector,pini,pend)
  self.pini = pini
  self.pend = pend
end

end # module


text:
module MambaModels

using Reexport

@reexport using Mamba
@reexport using CSV, DataFrames, Distributions
@reexport using MCMCChains, StatsFuns, StatsPlots

end # module


text:
using Test
using ThreadGantt

@testset "stuff" begin
    @test 1 == 1
end


text:
const SHIFT = 1
const REDUCEL = 2
const REDUCER = 3

# calculates action id for labeled reduce
# even numbers >= 2
reducel(label::Int) = label << 1
reducel() = REDUCEL

 # odd numbers >= 3
reducer(label::Int) = (label << 1) | 1
reducer() = REDUCER

# retrieve action type 1, 2 or 3
acttype(action::Int) = action == 1 ? SHIFT : 2 + (action & 1)

# retrieve label id
tolabel(action::Int) = action >> 1

type State{T <: ParserType}
    step::Int
    score::Union{Float64,Var}
    top::Int
    right::Int
    left::State{T}
    lchild::State{T}
    rchild::State{T}
    lsibl::State{T}
    rsibl::State{T}
    tokens::Vector{Token}
    parser::DepParser
    prev::State{T}
    prevact::Int
    head::State{T}
    feat::Vector{Int}

    function State(step, score, top, right)
        new(step, score, top, right)
    end

    function State(step, score, top, right, left,
        lchild, rchild, lsibl, rsibl, tokens, parser, prev, prevact, head)
        new(step, score, top, right, left,
            lchild, rchild, lsibl, rsibl, tokens, parser, prev, prevact, head)
    end
end

# nullstate is a State object whose all State type fields
# reference itself. used instead of Nullable{State} or nothing::Void
function nullstate_(T::Type)
    s = State{T}(0, 0.0, 0, 0)
    s.left, s.lchild, s.rchild = s, s, s
    s.lsibl, s.rsibl, s.prevact = s, s, reducel(1)
    s.head = s
    return s
end
const labeled_ = nullstate_(Labeled)
const unlabeled_ = nullstate_(Unlabeled)

nullstate(::Type{Labeled}) = labeled_
nullstate(::Type{Unlabeled}) = unlabeled_

# Base.isnull(s::State) = s.step == 0
Base.isnull{T}(s::State{T}) = s === nullstate(T)

function State{T}(tokens::Vector{Token}, parser::DepParser{T})
    score = parser.model == Perceptron ? 0.0 : Var([0f0])
    ns = nullstate(T)
    State{T}(1, score, 0, 1, ns, ns, ns,
        ns, ns, tokens, parser, ns, reducel(1), ns)
end

# to retrieve result
function heads(s::State)
    @assert isfinal(s)
    res = fill(-1, length(s.tokens))
    st = s
    while !isnull(st.prev)
        !isnull(st.lchild) && ( res[st.lchild.top] = st.top )
        !isnull(st.rchild) && ( res[st.rchild.top] = st.top )
        st = st.prev
    end
    @assert all(h -> h >= 0, res)
    res
end

function labels(s::State)
    @assert isfinal(s)
    res = fill(-1, length(s.tokens))
    st = s
    while !isnull(st.prev)
        atype = acttype(st.prevact)
        atype == REDUCEL && ( res[st.lchild.top] = tolabel(st.prevact) )
        atype == REDUCER && ( res[st.rchild.top] = tolabel(st.prevact) )
        st = st.prev
    end
    @assert all(l -> l >= 0, res)
    res
end

tokenat(s::State, i::Int) = get(s.tokens, i, roottoken)
tokenat(s::State, t::State) = get(s.tokens, t.top, roottoken)

function labelat(s::State, child::State)
    st = s
    while st.step != child.step + 1
        isnull(st.prev) && break
        st = st.prev
    end
    tolabel(st.prevact)
end

function expand{T}(s::State{T}, act::Int)
    atype = acttype(act)
    top, right, left, lchild, rchild, lsibl, rsibl = begin
        if atype == SHIFT
            ns = nullstate(T)
            s.right, s.right+1, s, ns, ns, ns, ns
        elseif atype == REDUCEL
            s.top, s.right, s.left.left, s.left, s.rchild, s, s.rsibl
        elseif atype == REDUCER
            s.left.top, s.right, s.left.left, s.left.lchild, s, s.left.lsibl, s.left
        else
            throw("Invalid action: $(act).")
        end
    end
    score = s.parser.model(s, act)
    State{T}(s.step+1, score, top, right, left, lchild,
             rchild, lsibl, rsibl, s.tokens, s.parser, s, act, nullstate(T))
end

# check if buffer is empty
bufferisempty(s::State) = s.right > length(s.tokens)

# check if can perform Reduce, that is, length(stack) >= 2
reducible(s::State) = !isnull(s.left)

# check if the state is final state
isfinal(s::State) = bufferisempty(s) && !reducible(s)

function isvalid(s::State, act::Int)
    if act == SHIFT
        return !bufferisempty(s)
    elseif act == REDUCEL
        return reducible(s) && s.left.top != 0
    elseif act == REDUCER
        return reducible(s)
    else
        throw("invalid action: $act")
    end
end

function expand_reduce!(res::Vector{State{Labeled}}, s::State{Labeled}, act::Function)
    for label in 1:length(s.parser.labels)
        push!(res, expand(s, act(label)))
    end
end

function expand_reduce!(res::Vector{State{Unlabeled}}, s::State{Unlabeled}, act::Function)
    push!(res, expand(s, act()))
end

function expandpred{T}(s::State{T})
    isfinal(s) && return []
    res = State{T}[]
    if reducible(s)
        expand_reduce!(res, s, reducer)
        if s.left.top != 0
            expand_reduce!(res, s, reducel)
        end
    end
    bufferisempty(s) || push!(res, expand(s, SHIFT))
    res
end

function expand_reduce(s::State{Unlabeled}, act::Function)
    expand(s, act())
end

function expand_reduce(s::State{Labeled}, act::Function)
    child = act == reducel ? s.left.top :
            act == reducer ? s.top :
            throw("Action must be reduce")
    expand(s, act(tokenat(s, child).label))
end

function expandgold{T}(s::State{T})
    isfinal(s) && return []
    if !reducible(s)
        return [expand(s, SHIFT)]
    else
        s0 = tokenat(s, s)
        s1 = tokenat(s, s.left)
        if s1.head == s.top
            return [expand_reduce(s, reducel)]
        elseif s0.head == s.left.top
            if all(i -> tokenat(s, i).head != s.top, s.right:length(s.tokens))
                return [expand_reduce(s, reducer)]
            end
        end
    end
    return [expand(s, SHIFT)]
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
"""
    GreensMatrix{eltype, mattype}

`GreensMatrix` is a thin wrapper for `mattype` eturned by `greens(...)`,
`greens!(...)` and related iterators. The elements `G[i, j]` represent
`⟨cᵢ(k) cⱼ^†(l)⟩` where `k` and `l` are time slice indices contained in the
type.

Related to this is the `P::Permuted` type, returned by `swapop(::GreensMatrix)`.
It represents `P[i, j] = ⟨cᵢ^†(l) cⱼ(k)⟩ = δᵢⱼ δₖₗ - ⟨cⱼ(k) cᵢ^†(l)⟩`.
"""
struct GreensMatrix{T, M <: AbstractMatrix{T}} <: AbstractMatrix{T}
    k::Int
    l::Int
    val::M
end

"""
see GreensMatrix
"""
struct Permuted{T, GT <: GreensMatrix{T}} <: AbstractMatrix{T}
    x::GT
end

# function Base.show(io::IO, x::GreensMatrix)
#     println(io, "G[i, j] = cᵢ(", x.k, ") cⱼ^†(", x.l, ") =")
#     show(io, x.val)
# end
# function Base.show(io::IO, x::Permuted{<: GreensMatrix})
#     println(io, "G'[i, j] = cᵢ^†(", x.l, ") cⱼ(", x.k, ") =")
#     show(io, I[x.x.k, x.x.l] * I - transpose(x.x.val))
# end

Base.size(x::GreensMatrix) = size(x.val)
Base.size(x::Permuted) = size(x.x.val)
Base.getindex(x::GreensMatrix, i, j) = x.val[i, j]
Base.getindex(x::Permuted, i, j) = I[x.x.k, x.x.l] * I[i, j] - x.x.val[j, i]

"""
    swapop(G::GreensMatrix)

Permute the operator order of a `GreensMatrix` from ⟨cᵢ(k) cⱼ^†(l)⟩ to
⟨cᵢ^†(l) cⱼ(k)⟩.
"""
swapop(x::GreensMatrix) = Permuted(x)
swapop(x::Permuted) = x.x
Base.copy(x::GreensMatrix) = GreensMatrix(x.k, x.l, copy(x.val))
function Base.:(==)(a::GreensMatrix, b::GreensMatrix)
    a.k == b.k && a.l == b.l && a.val == b.val
end
function Base.isapprox(a::GreensMatrix, b::GreensMatrix; kwargs...)
    a.k == b.k && a.l == b.l && isapprox(a.val, b.val; kwargs...)
end

"""
    greens(mc::DQMC)

Obtain the current equal-time Green's function, i.e. the fermionic expectation
value of `Gᵢⱼ = ⟨cᵢcⱼ^†⟩`. The indices relate to sites and flavors, but the
exact meanign depends on the model. For the attractive Hubbard model
`G[i, j] = ⟨c_{i, ↑} c_{j, ↑}^†⟩ = ⟨c_{i, ↓} c_{j, ↓}^†⟩` due to symmetry.

Internally, `mc.stack.greens` is an effective Green's function. This method
transforms it to the actual Green's function by multiplying hopping matrix
exponentials from left and right.
"""
@bm greens(mc::DQMC) = GreensMatrix(0, 0, copy(_greens!(mc)))

"""
    greens!(mc::DQMC[; output=mc.stack.greens_temp, input=mc.stack.greens, temp=mc.stack.Ur])

Inplace version of `greens`.
"""
@bm function greens!(mc::DQMC; output=mc.stack.greens_temp, input=mc.stack.greens, temp=mc.stack.Ur)
    # TODO rework measurements to work well with StructArrays and remove this
    if isdefined(mc.stack, :complex_greens_temp)
        _greens!(mc, output, input, temp)
        copyto!(mc.stack.complex_greens_temp, output)
        GreensMatrix(0, 0, mc.stack.complex_greens_temp)
    else
        GreensMatrix(0, 0, _greens!(mc, output, input, temp))
    end
end

# Implementations with and without checkerboards.
# Note that we use the Trotter decomposition
# `exp(-Δτ(V + T)) = exp(-0.5 Δτ T) exp(-Δτ V) exp(-0.5 Δτ T)`
# alongside the cyclic property of determinants
# `1 / det(I + eT eV1 eT ... eT eVN eT) = 1 / det(I * eT eT eV1 ... eT eT eVN)`
# to increase analytic accuracy. When calculating the greens function we need to
# reverse the cyclic permutation, which happens here.

function _greens!(
        mc::DQMC_CBFalse, target::AbstractMatrix = mc.stack.greens_temp,
        source::AbstractMatrix = mc.stack.greens, temp::AbstractMatrix = mc.stack.Ur
    )
    eThalfminus = mc.stack.hopping_matrix_exp
    eThalfplus = mc.stack.hopping_matrix_exp_inv
    vmul!(temp, source, eThalfminus)
    vmul!(target, eThalfplus, temp)
    return target
end

function _greens!(
        mc::DQMC_CBTrue, target::AbstractMatrix = mc.stack.greens_temp,
        source::AbstractMatrix = mc.stack.greens, temp::AbstractMatrix = mc.stack.Ur
    )
    chkr_hop_half_minus = mc.stack.chkr_hop_half
    chkr_hop_half_plus = mc.stack.chkr_hop_half_inv
    copyto!(target, source)

    @inbounds @views begin
        for i in reverse(1:mc.stack.n_groups)
            vmul!(temp, target, chkr_hop_half_minus[i])
            copyto!(target, temp)
        end
        for i in reverse(1:mc.stack.n_groups)
            vmul!(temp, chkr_hop_half_plus[i], target)
            copyto!(target, temp)
        end
    end
    return target
end

# Same stuff with a specified time slice.

"""
    greens(mc::DQMC, l::Integer)

Calculates the equal-time greens function at a given slice index `l`, i.e.
`G_{ij}(l, l) = G_{ij}(l⋅Δτ, l⋅Δτ) = ⟨cᵢ(l⋅Δτ)cⱼ(l⋅Δτ)^†⟩`.

Note: This internally overwrites the stack variables `Ul`, `Dl`, `Tl`, `Ur`,
`Dr`, `Tr`, `curr_U`, `tmp1` and `tmp2`. All of those can be used as temporary
or output variables here, however keep in mind that other results may be
invalidated. (I.e. `G = greens!(mc)` would get overwritten.)
"""
@bm greens(mc::DQMC, slice::Integer) = GreensMatrix(slice, slice, copy(_greens!(mc, slice)))

"""
    greens!(mc::DQMC, l::Integer[; output=mc.stack.greens_temp, temp1=mc.stack.tmp1, temp2=mc.stack.tmp2])

Inplace version of `greens!(mc, l)`
"""
@bm function greens!(
        mc::DQMC, slice::Integer;
        output=mc.stack.greens_temp, temp1 = mc.stack.tmp1, temp2 = mc.stack.tmp2
    )
    GreensMatrix(slice, slice, _greens!(mc, slice, output, temp1, temp2))
end

function _greens!(
        mc::DQMC, slice::Integer, output::AbstractMatrix = mc.stack.greens_temp,
        temp1::AbstractMatrix = mc.stack.tmp1, temp2::AbstractMatrix = mc.stack.tmp2
    )
    calculate_greens(mc, slice, temp1)
    _greens!(mc, output, temp1, temp2)
end


text:
module TestContexts

using Test

using QXContexts

include("utils.jl")

@testset "Test contexts module" begin

    @testset "Test SliceIterator" begin
        dims = [2,3,4]
        si = SliceIterator(dims)
        @test length(si) == prod(dims)
        @test si[1] == [1,1,1]
        @test si[dims...] == [2,3,4]

        si = SliceIterator(dims, 1, 10)
        @test length(si) == 10

        si = SliceIterator(dims)
        @test length(si) == prod(dims)
        si2 = SliceIterator(si, 1, 5)
        @test length(si2) == 5

        si3 = SliceIterator(si2, 3, 5)
        @test length(si3) == 3
    end

    @testset "Test QXContext" begin
        cg = ComputeGraph(build_tree(sample_cmds), deepcopy(sample_tensors))
        ctx = QXContext(cg)

        # contract without setting view parameters
        set_open_bonds!(ctx, "0")
        output_0 = ctx()
        @test size(output_0) == (2,)

        set_open_bonds!(ctx, "1")
        output_1 = ctx()
        @test size(output_1) == (2,)

        # contract while summing over view parameters
        @test compute_amplitude!(ctx, "0") ≈ output_0
        @test compute_amplitude!(ctx, "1") ≈ output_1

        @test compute_amplitude!(ctx, "1", max_slices=0) ≈ output_1
    end
end

end


text:
function r2(x)
    nrows=size(x,1)
    out = zeros(nrows,1)
    for i = 1:nrows
        out[i,1] = sum((x[i,:]).^2.)^(0.5)
    end
    out
end

function f(w,alpha)
    if abs(w) <= 1/alpha
        return 2*alpha*sign(w)-2*alpha^2*w
    end
    return 0
end

function dJ(x,alpha)
    (nrows,ncols) = size(x)
    out = zeros(nrows,ncols)

    for i = 1:nrows
        l2xi  = sum(x[i,:].^2)^0.5 #euclidean distance of the i'th row of X
        if l2xi !=0  #if the distance is not 0
            for j = 1:ncols #set the columns of out to non-zero
                out[i,j] = f(l2xi,alpha)/l2xi*x[i,j]
            end
        end
    end
    out
end

function ZAP(Y::Array{Float64,2},A::Array{Float64,2};maxIter=500,kappa=.1,neu=.1,threshold=1e-6,verbose=false,Q=11,alpha=1)
    #initialize guess
    A_inv = pinv(A)
    guess = A_inv*Y

    if verbose
        print("\nMax Iterations: ",maxIter,"\nCurrent Iteration:\n")
    end

    iteration=0
    for iteration = 1:maxIter

        if mod(iteration ,int(Q))==0
            previous_guess = deepcopy(guess)

            if verbose
                print("\r"*string("",iteration ,"\tKappa = ",kappa,"\t",
                round(maximum([iteration /maxIter*100.,log(threshold,10)/log(kappa,10)*100.]),1),"%"))
            end
        end

        guess = guess - kappa * dJ(guess,alpha)
        guess = guess + A_inv*(Y-A*guess)
        if mod(iteration ,int(Q))==0 && sum(r2(guess))>=sum(r2(previous_guess))
            kappa=neu*kappa

            if kappa<= threshold
                break
            end
        end
    end
    if verbose
        print("\n\nDone.\n")
    end
    return guess

end

max_stars_repo_name: dahlend/CompressedSensing


text:
const defaultTransitions = ["attrition", "retirement"]
const specialNodes = ["", "external", "out"]

function validateTransition( mpSim::MPsim, transition::String )
    if lowercase( transition ) ∈ defaultTransitions
        return true
    end  # if lowercase( transition ) ∈ defaultTransitions

    # Transition name can be a recruitment, or another transition.
    return haskey( mpSim.recruitmentByName, transition ) ||
        haskey( mpSim.transitionsByName, transition )
end  # validateTransition( mpSim, transition )

function validateTransition( mpSim::MPsim, transition::NTuple{2,String} )
    sourceNode, targetNode = transition

    # If the source is external, the transition is a recruitment.
    if lowercase( sourceNode ) ∈ specialNodes
        return haskey( mpSim.recruitmentByTarget, targetNode )
    end  # if lowercase( sourceNode ) ∈ specialNodes

    # If the target is external, the transition is an OUT transition.
    if lowercase( targetNode ) ∈ specialNodes
        return true
    end  # if lowercase( targetNode ) ∈ specialNodes

    if !haskey( mpSim.transitionsBySource, sourceNode )
        return false
    end  # if !haskey( mpSim.transitionsBySource, sourceNode )

    return any( transition -> transition.targetNode == targetNode,
        mpSim.transitionsBySource[sourceNode] )
end  # validateTransition( mpSim, transition )

function validateTransition( mpSim::MPsim, transition::NTuple{3,String} )
    name, sourceNode, targetNode = transition

    # If the source is external, the transition is a recruitment.
    if lowercase( sourceNode ) ∈ specialNodes
        return haskey( mpSim.recruitmentByName, name ) &&
            any( recruitment -> recruitment.targetNode == targetNode,
                mpSim.recruitmentByName[name] )
    end  # if lowercase( sourceNode ) ∈ specialNodes

    # Default transitions (attrition and retirement).
    if lowercase( name ) ∈ defaultTransitions
        return lowercase( targetNode ) ∈ specialNodes
    end  # if lowercase( name ) ∈ defaultTransitions

    # If the target is external, the transition is OUT.
    if lowercase( targetNode ) ∈ specialNodes
        return any( transition -> ( transition.name == name ) &&
            ( transition.sourceNode == sourceNode ),
            mpSim.transitionsByTarget["OUT"] )
    end  # if lowercase( targetNode ) ∈ specialNodes

    return haskey( mpSim.transitionsByTarget, targetNode ) &&
        any( transition -> ( transition.name == name ) &&
            ( transition.sourceNode == sourceNode ),
            mpSim.transitionsByTarget[targetNode] )
end  # validateTransition( mpSim, transition )

function createTransitionFluxReport( simDB::SQLite.DB, idKey::String,
    transDBname::String, timeGrid::Vector{Float64}, transition::TransitionType )
    queryPartCmd = string( "SELECT count( `", idKey,
        "` ) counts FROM `", transDBname, "` WHERE\n    ",
        generateTransitionQuery( transition ), " AND\n    " )

    result = zeros( Int, length( timeGrid ) )

    for ii in eachindex( timeGrid )
        queryCmd = string( queryPartCmd, generateTimeFork( ii == 1 ?
            timeGrid[1] : timeGrid[ii - 1], timeGrid[ii] ) )
        result[ii] = DataFrame( DBInterface.execute( simDB, queryCmd ) )[1, 1]
    end  # for ii in eachindex( timeGrid )

    return result
end  # createTransitionFluxReport( simDB, idKey, transDBname, timeGrid,
     #    transition )

createTransitionFluxReport( mpSim::MPsim, timeGrid::Vector{Float64},
    transition::TransitionType ) = createTransitionFluxReport( mpSim.simDB,
    mpSim.idKey, mpSim.transDBname, timeGrid, transition )

#=
function createTransitionFluxReport( mpSim::MPsim, timeGrid::Vector{Float64},
    transition::TransitionType )

    queryPartCmd = string( "SELECT count( `", mpSim.idKey,
        "` ) counts FROM `", mpSim.transDBname, "` WHERE\n    ",
        generateTransitionQuery( transition, mpSim ), " AND\n    " )

    result = zeros( Int, length( timeGrid ) )

    for ii in eachindex( timeGrid )
        queryCmd = string( queryPartCmd, generateTimeFork( ii == 1 ?
            timeGrid[1] : timeGrid[ii - 1], timeGrid[ii] ) )
        result[ii] =
            DataFrame( DBInterface.execute( mpSim.simDB, queryCmd ) )[1, 1]
    end  # for ii in eachindex( timeGrid )

    return result

end  # createTransitionFluxReport( mpSim, timeGrid, transition )
=#

generateTransitionQuery( transition::String ) =
    string( "transition IS '", transition, "'" )
    # string( "transition IS '", transition, "' AND",
    # "\n    sourceNode IS NOT 'active' AND",
    # "\n    targetNode IS NOT 'active'" )

function generateTransitionQuery( transition::NTuple{2,String} )
    sourceNode = lowercase( transition[1] ) ∈ specialNodes ? "NULL" :
        string( "'", transition[1], "'" )
    targetNode = lowercase( transition[2] ) ∈ specialNodes ? "NULL" :
        string( "'", transition[2], "'" )
    return string( "sourceNode IS ", sourceNode, " AND",
        "\n    targetNode IS ", targetNode )
end  # generateTransitionQuery( transition )

generateTransitionQuery( transition::NTuple{3,String} ) =
    string( "transition IS '", transition[1], "' AND",
    "\n    ", generateTransitionQuery( transition[2:3] ) )

#=
generateTransitionQuery( transition::String, mpSim::MPsim ) =
    string( "transition IS '", transition, "' AND",
    "\n    ", mpSim.sNode, " IS NOT 'active' AND",
    "\n    ", mpSim.tNode, " IS NOT 'active'" )

function generateTransitionQuery( transition::NTuple{2,String}, mpSim::MPsim )

    sourceNode = lowercase( transition[1] ) ∈ specialNodes ? "NULL" :
        string( "'", transition[1], "'" )
    targetNode = lowercase( transition[2] ) ∈ specialNodes ? "NULL" :
        string( "'", transition[2], "'" )
    return string( mpSim.sNode, " IS ", sourceNode, " AND",
        "\n    ", mpSim.tNode, " IS ", targetNode )

end  # generateTransitionQuery( transition, mpSim )

generateTransitionQuery( transition::NTuple{3,String}, mpSim::MPsim ) =
    string( "transition IS '", transition[1], "' AND",
    "\n    ", generateTransitionQuery( transition[2:3], mpSim ) )
=#

generateTransitionName( transition::String ) =
    lowercase( transition ) ∈ defaultTransitions ? lowercase( transition ) :
    transition

function generateTransitionName( transition::NTuple{2,String} )
    sourceNode, targetNode = transition
    sourceNode = lowercase( sourceNode ) ∈ specialNodes ? "external" :
        sourceNode
    targetNode = lowercase( targetNode ) ∈ specialNodes ? "external" :
        targetNode
    return string( sourceNode, " => ", targetNode )
end  # generateTransitionName( transition )

generateTransitionName( transition::NTuple{3,String} ) =
    string( generateTransitionName( transition[1] ), ": ",
    generateTransitionName( transition[2:3] ) )


text:
using Test
using LinearAlgebra
using QCMaterial

import PyCall: pyimport
import Random

#@testset "ucc.num_theta" begin
#    n_qubit = 1
#    theta = 1.0
#    c = UCCQuantumCircuit(n_qubit)
#    add_parametric_multi_Pauli_rotation_gate!(
#            c.circuit, [1], [pauli_Y], theta)
#    println(c.thetas)
#    @test num_theta(c) == 1
#end

@testset "ucc.uccgsd" begin
    #Random.seed!(1)
    scipy_opt = pyimport("scipy.optimize")

    nsite = 2
    n_qubit = 2*nsite
    U = 1.0
    t = -0.01
    μ = 0

    ham = FermionOperator()
    for i in 1:nsite
        up = up_index(i)
        down = down_index(i)
        ham += FermionOperator("$(up)^ $(down)^ $(up) $(down)", -U)
    end

    for i in 1:nsite-1
        ham += FermionOperator("$(up_index(i+1))^ $(up_index(i))", t)
        ham += FermionOperator("$(up_index(i))^ $(up_index(i+1))", t)
        ham += FermionOperator("$(down_index(i+1))^ $(down_index(i))", t)
        ham += FermionOperator("$(down_index(i))^ $(down_index(i+1))", t)
    end

    for i in 1:nsite
        up = up_index(i)
        down = down_index(i)
        ham += FermionOperator("$(up)^  $(up) ", -μ)
        ham += FermionOperator("$(down)^ $(down)", -μ)
    end

    n_electron = 2
    @assert mod(n_electron, 2) == 0
    sparse_mat = get_number_preserving_sparse_operator(ham, n_qubit, n_electron);

    enes_ed = eigvals(sparse_mat.toarray());

    ham_jw = jordan_wigner(ham)

    circuit = uccgsd(n_qubit, orbital_rot=true)

    function cost(theta_list)
        state = create_hf_state(n_qubit, n_electron)
        update_circuit_param!(circuit, theta_list)
        update_quantum_state!(circuit, state)
        get_expectation_value(ham_jw, state)
    end

    theta_init = rand(num_theta(circuit))
    cost_history = Float64[]
    init_theta_list = theta_init
    push!(cost_history, cost(init_theta_list))

    method = "BFGS"
    options = Dict("disp" => true, "maxiter" => 200, "gtol" => 1e-5)
    callback(x) = push!(cost_history, cost(x))
    opt = scipy_opt.minimize(cost, init_theta_list, method=method, callback=callback)
    println("Eigval_vqe=", cost_history[end])
    EigVal_min = minimum(enes_ed)
    println("EigVal_min=", EigVal_min)
    @test abs(EigVal_min-cost_history[end]) < 1e-3
end

@testset "ucc.UCCQuantumCircuit" begin
    n_qubit = 4
    c = uccgsd(n_qubit, orbital_rot=true)
    c.thetas .= 1.0
    c_copy = copy(c)
    @test all(c_copy.thetas == c.thetas)
end

@testset "ucc.one_rotation_gate" begin
    using LinearAlgebra
    n_qubit = 1
    theta = 1e-5

    c = UCCQuantumCircuit(n_qubit)
    add_parametric_multi_Pauli_rotation_gate!(
            c.circuit, [1], [pauli_Y], theta)

    state = QulacsQuantumState(n_qubit, 0b1)
    update_quantum_state!(c, state)
    vec = get_vector(state)
    vec ≈ [0.5*theta, 1.0]
end

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
abstract type AbstractSampleMethod end
struct SampleUniform        <: AbstractSampleMethod end
struct SampleBinCenter      <: AbstractSampleMethod end

const SAMPLE_UNIFORM          = SampleUniform()
const SAMPLE_BIN_CENTER       = SampleBinCenter()


text:
using GaussianRandomFields;
using Plots;
using Plots.PlotMeasures;
using Colors;
using Netpbm;
using FileIO;

width = 1024
height = 1024
type = "CirculantEmbedding"

function generaterandomfield(type)
    if type == "Cholesky"
        cov = CovarianceFunction(2, Matern(1 / 2, 5 / 4))
        pts = range(0, stop = 1, length = 101)
        grf = GaussianRandomField(cov, Cholesky(), pts, pts)
    elseif type == "Spectral"
        cov = CovarianceFunction(2, Linear(1))
        pts = range(0, stop = 1, length = 51)
        grf = GaussianRandomField(cov, Spectral(), pts, pts)
    elseif type == "KarhunenLoeve"
        cov = CovarianceFunction(2, Whittle(.1))
        pts = range(0, stop = 1, length = 201)
        grf = GaussianRandomField(cov, KarhunenLoeve(500), pts, pts)
    elseif type == "CirculantEmbedding"
        cov = CovarianceFunction(2, Exponential(.5))
        pts = range(0, stop = 1, length = 1001)
        grf = GaussianRandomField(cov, CirculantEmbedding(), pts, pts, minpadding = 2001)
    end
end

function normalize(a)
    (minval, maxval) = extrema(a)
    valrange = maxval - minval

    function normalize0(v)
        (v - minval) / valrange
    end

    normalize0.(a)
end

println("Generating random field...")
@time grf = generaterandomfield(type)

println("Saving random field...")
data = sample(grf)

plot_options = (legend = false, c = :grays, framestyle = :none, size = (width,height), margin = -2mm, bottom_margin = -10mm, left_margin = -13mm)
heatmap_plot = heatmap(data; plot_options...)
contour_plot = contourf(data; plot_options...)
contour8_plot = contourf(data; levels = 8, plot_options...)
png(heatmap_plot, "heatmap.png")
png(contour_plot, "contour.png")
png(contour8_plot, "contour8.png")

normalized_data = reverse(normalize(data), dims = 1)
pgm_image_data = convert(Array{Gray}, normalized_data)
pgm_file = File(format"PGMBinary", "heatmap.pgm")
Netpbm.save(pgm_file, pgm_image_data)

println("Done")

max_stars_repo_name: ngyewch/fjage-sim-sentuator

max_issues_repo_name: ngyewch/fjage-sim-sentuator

max_forks_repo_name: ngyewch/fjage-sim-sentuator


text:
# Stochastic gradient optimizers

abstract type Optimizer end

function training_step!(opt,θ,G)
    optimizer_step!(opt,G)
    θ + opt.Δ
end

"""
Maximize likelihood with stochastic gradient

- dℓ(θ) returns loglikelihood and gradient
"""
function train(dℓ,θ0,opt;max_iter=1000, verbose=false, progress=false)
    θ = zeros(length(θ0))
    θ .= θ0[:]
    ℓ = zeros(max_iter)
    if progress
        p = Progress(max_iter,1)
    end
    for i in 1:max_iter
        ℓ[i],G = dℓ(θ)
        θ .= training_step!(opt,θ,G)
        if progress
            next!(p)#; showvalues = [(:iteration, i),(:likelihood,ℓ[i])])
        end
        if verbose
            println("Step: ",rpad(i,5,' ')," | Likelihood: ", ℓ[i])
        end
    end
    θ,ℓ
end

mutable struct SGD <: Optimizer
    η
    δ
    Δ::Array{Float64,1}
end

SGD(n::Int,η=0.001,δ=1.0) = SGD(η,δ,zeros(n))

function optimizer_step!(opt::SGD,G)
    opt.Δ .= opt.η*G
    opt.η = opt.η*opt.δ
    opt
end

mutable struct Momentum <: Optimizer
    η
    ρ
    δ
    Δ::Array{Float64,1}
end

Momentum(n::Int,η=0.001,ρ=0.9,δ=1.0) = Momentum(η,ρ,δ,zeros(n))

function optimizer_step!(opt::Momentum,G)
    opt.Δ = opt.ρ*opt.Δ + opt.η*G
    opt.η = opt.η*opt.δ
    opt
end

mutable struct Adam <: Optimizer
    η
    δ
    β1
    β2
    βp
    m
    v
    Δ::Array{Float64,1}
end

Adam(n::Int,η=0.001,δ=1.0,β1=0.9,β2=0.999) = Adam(η,δ,β1,β2,[β1;β2],zeros(n),zeros(n),zeros(n))

function optimizer_step!(opt::Adam,G)
    opt.m .= opt.β1*opt.m + (1-opt.β1)*G
    opt.v .= opt.β2*opt.v + (1-opt.β2)*(G.^2)
    opt.Δ .= opt.η*(opt.m/(1-opt.βp[1]))./(sqrt.(opt.v/(1-opt.βp[2])).+1e-8)
    opt.βp .*= [opt.β1,opt.β2]
    opt.η = opt.δ*opt.η
    opt
end


text:
using GLVisualize, GeometryTypes, Reactive

if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0f0,1f0,360))
end
N = 7
# generate some rotations
function rotation_func(t)
    t = (t == 0f0 ? 0.01f0 : t)
    Vec3f0[(sin(x/t), cos(y/(t/2f0)), sqrt(t+z^2)) for x=1:N, y=1:N, z=1:N]
end

# us Reactive.map to transform the timesignal signal into the arrow flow
flow = map(rotation_func, timesignal)

# create a visualisation
s=step(linspace(0.,1.,N))
vis = visualize(flow, scale=Vec3f0(s,s,0.001))
_view(vis, window)

if !isdefined(:runtests)
    renderloop(window)
end


text:
# Utilities for working with HSA queues

"""
An HSA queue.
Each queue is uniquely associated with an agent.
"""
mutable struct HSAQueue
    agent::HSAAgent
    queue::Ref{Ptr{HSA.Queue}}
end

const DEFAULT_QUEUES = IdDict{HSAAgent,HSAQueue}()

function HSAQueue(agent::HSAAgent)
    queue_size = Ref{UInt32}(0)
    getinfo(agent.agent, HSA.AGENT_INFO_QUEUE_MAX_SIZE, queue_size) |> check
    @assert queue_size[] > 0
    queue = HSAQueue(agent, Ref{Ptr{HSA.Queue}}())
    HSA.queue_create(agent.agent, queue_size[], HSA.QUEUE_TYPE_SINGLE,
                     C_NULL, C_NULL, typemax(UInt32), typemax(UInt32),
                     queue.queue) |> check

    hsaref!()
    finalizer(queue) do queue
        HSA.queue_destroy(queue.queue[]) |> check
        hsaunref!()
    end
    return queue
end

get_default_queue() = get_default_queue(get_default_agent())

"Gets the default queue for an agent, creating it if necessary"
function get_default_queue(agent::HSAAgent)
    if VERSION >= v"1.1"
        return get!(DEFAULT_QUEUES, agent) do
            HSAQueue(agent)
        end
    else
        if haskey(DEFAULT_QUEUES, agent)
            return DEFAULT_QUEUES[agent]
        else
            return HSAQueue(agent)
        end
    end
end


text:
export mipReg

function mipReg(m::Vector,mref,M::AbstractMesh;Iact=I, q::Int64=1)
	# Rc = e'* m.^3.*(1-m).^3
	Rc = sum(m.^q .* (1.0 .- m).^q)
	dR = q * m.^(q-1) .* (1.0 .- m).^q - q .* m.^q .*(m .- 1.0).^(q-1)
	# d2R = SparseMatrixCSC(Diagonal(max.(0.0 ,- 6.0*m.*(m .- 1.).^3 - 3.0 * m.^3 .* (2.0*m .- 2.) - 18.0*m.^2.0 .* (m .- 1.0).^2)))
	d2R = SparseMatrixCSC(Diagonal( (q*(q-1))*m.^(max(0,q-2)) .*(1.0 .- m).^q -
	                                (q*(q-1)) * m.^q .* (m .- 1.0).^(max(0,q-2))  -
								     2*(q*q)*m.^(q-1) .* (m .- 1.0).^(q-1)))
	d2R = max.(d2R,0.0)
	return Rc,dR,d2R
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
module MeanTests
using Base.Test;
using FactCheck;
using ValidatedNumerics;
using IntervalStatistic;
using Distributions;

srand(10)

in_3_sigma_interval(mu, sigma) = (i) -> (mu in i) && (i in mu + 3*@interval(-sigma, sigma))

facts("estimate average of standard (0, 1) distribution") do
    d = Normal()
    length = 100
    values = rand(d, length)
    mu, sigma = params(d)
    average = mean(values)

    context("by known variance") do
        result = IntervalStatistic.mean(values, IntervalStatistic.Mean.byKnownVariance(0.95, sigma*sigma))
        println(result)
        @fact result --> in_3_sigma_interval(mu, sigma)
    end

    context("by interquartile range") do
        result = IntervalStatistic.mean(values, IntervalStatistic.Mean.byInterQuartileRange(0.95))
        println(result)
        @fact result --> in_3_sigma_interval(mu, sigma)
    end

    context("by unknown variance") do
        result = IntervalStatistic.mean(values, IntervalStatistic.Mean.byUnknownVariance(0.95))
        println(result)
        @fact result --> in_3_sigma_interval(mu, sigma)
    end

    context("by standard error") do
        result = IntervalStatistic.mean(values, IntervalStatistic.Mean.byStandardError())
        println(result)
        @fact result --> in_3_sigma_interval(mu, sigma)
    end

    context("by mean abs deviation") do
        result = IntervalStatistic.mean(values, IntervalStatistic.Mean.byMeanAbsDeviation(0.95))
        println(result)
        @fact result --> in_3_sigma_interval(mu, sigma)
    end
end

facts("estimate average of normal mu=3, sigma=0.1 distribution") do
    d = Normal(3, 0.1)
    length = 100
    values = rand(d, length)
    mu, sigma = params(d)
    average = mean(values)

    context("by known variance") do
        result = IntervalStatistic.mean(values, IntervalStatistic.Mean.byKnownVariance(0.95, sigma*sigma))
        println(result)
        @fact result --> in_3_sigma_interval(mu, sigma)
    end

    context("by interquartile range") do
        result = IntervalStatistic.mean(values, IntervalStatistic.Mean.byInterQuartileRange(0.95))
        println(result)
        @fact result --> in_3_sigma_interval(mu, sigma)
    end

    context("by standard error") do
        result = IntervalStatistic.mean(values, IntervalStatistic.Mean.byStandardError())
        println(result)
        @fact result --> in_3_sigma_interval(mu, sigma)
    end

    context("by unknown variance") do
        result = IntervalStatistic.mean(values, IntervalStatistic.Mean.byUnknownVariance(0.95))
        println(result)
        @fact result --> in_3_sigma_interval(mu, sigma)
    end

    context("by mean abs deviation") do
        result = IntervalStatistic.mean(values, IntervalStatistic.Mean.byMeanAbsDeviation(0.95))
        println(result)
        @fact result --> in_3_sigma_interval(mu, sigma)
    end
end

end


text:
function vecadd1(a,b,c,N)
  for i = 1:N
    c = a + b
  end
end

function vecadd2(a,b,c,N)
  for i = 1:N, j=1:length(c)
    c[j] = a[j] + b[j]
  end
end

A = rand(2)
B = rand(2)
C = zeros(2)

time = @elapsed vecadd1(A,B,C,100000000)
println(time)
C = zeros(2)
time = @elapsed vecadd2(A,B,C,100000000)
println(time)


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

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

function sumtest_nobounds(a::Array{Float64, 1})
  sum = 0.0
  for idx=1:length(a)
    sum += a[idx]
  end
  sum
end

function sumtest_inbounds(a::Array{Float64, 1})
  sum = 0.0
  for idx=1:length(a)
    @inbounds sum += a[idx]
  end
  sum
end

function sumtest_prevar(a::Array{Float64, 1})
  l::Int64 = length(a)
  sum = 0.0
  for idx = 1:l
    @inbounds sum += a[idx]
  end
  sum
end

count = 10000

alist = [rand(1000) for idx = 1:count]

function testwith(f, aofa)
  for idx=1:count
    f(aofa[idx])
  end
end

function run()
@time testwith(sumtest_nobounds, alist)
@time testwith(sumtest_inbounds, alist)
@time testwith(sumtest_prevar, alist)
end

for idx = 1:10
  run()
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
@inline cos(x::Float32) = ccall("extern @air.cos.f32", llvmcall, Cfloat, (Cfloat,), x)
@inline cos_fast(x::Float32) = ccall("extern @air.fast_cos.f32", llvmcall, Cfloat, (Cfloat,), x)

@inline sin(x::Float32) = ccall("extern @air.sin.f32", llvmcall, Cfloat, (Cfloat,), x)
@inline sin_fast(x::Float32) = ccall("extern @air.fast_sin.f32", llvmcall, Cfloat, (Cfloat,), x)

@inline exp(x::Float32) = ccall("extern @air.exp.f32", llvmcall, Cfloat, (Cfloat,), x)
@inline exp_fast(x::Float32) = ccall("extern @air.fast_exp.f32", llvmcall, Cfloat, (Cfloat,), x)


text:
using Genie.Router

route("/") do
  serve_static_file("welcome.html")
end

route("/error500") do
  error_500("Something went wrong")
end

route("/error404") do
  error_404("the page you want")
end


text:
using Plots
pyplot()

x = 1:10
y = rand(10)
plot(x,y)

n = 10
x = [rand()+1 * randn(n) + 2i for i in 1:5]
y = [rand()+1 * randn(n) + i for i in 1:5]

f(v) = 1.96std(v) / sqrt(n)
xerr = map(f,x)
yerr = map(f,y)
x = map(mean, x)
y = map(mean, y)

plot(x,y, xerr = xerr, yerr = yerr, marker=stroke(2, :orange))


text:
module Example #hide
# # A first demo with fancy names
# This example does more or less the same things as the first example, but using non regular API function names.
#
# This example serves as a test since this project doesn't have a "testing" procedure yet. In this example,
# the equation ``u'(x) = 2`` with ``u(0) = 0`` is solved on the domain ``[0,1]`` using a backward finite
# difference scheme.
#
# Note that the way we achieve things in the document can be highly improved and the purpose of this example
# is only demonstrate some method calls to give an overview.
#
# To run this example, execute : `mpirun -n your_favorite_positive_integer julia example2.jl`

# Import package
using PetscWrap

# Initialize PETSc. Command line arguments passed to Julia are parsed by PETSc. Alternatively, you can
# also provide "command line arguments by defining them in a string, for instance
# `PetscInitialize("-ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always")` or by providing each argument in
# separate strings : `PetscInitialize(["-ksp_monitor_short", "-ksp_gmres_cgs_refinement_type", "refine_always")`
PetscInitialize()

# Number of mesh points and mesh step
n = 11
Δx = 1. / (n - 1)

# Create a matrix of size `(n,n)` and a vector of size `(n)`
A = create_matrix(n, n)
b = create_vector(n)

# We can then use command line options to set our matrix/vectors.
set_from_options!(A)
set_from_options!(b)

# Finish the set up
set_up!(A)
set_up!(b)

# Let's build the right hand side vector. We first get the range of rows of `b` handled by the local processor.
# The `rstart, rend = get_range(*)` methods differ a little bit from PETSc API since the two integers it
# returns are the effective Julia range of rows handled by the local processor. If `n` is the total
# number of rows, then `rstart ∈ [1,n]` and `rend` is the last row handled by the local processor.
b_start, b_end = get_range(b)

# Now let's build the right hand side vector. Their are various ways to do this, this is just one.
n_loc = length(b) ## Note that n_loc = b_end - b_start + 1...
b[b_start:b_end] = 2 * ones(n_loc)

# And here is the differentiation matrix. Rembember that PETSc.MatSetValues simply ignores negatives rows indices.
A_start, A_end = get_range(A)
for i in A_start:A_end
    A[i, i-1:i] = [-1. 1.] / Δx
end

# Set boundary condition (only the proc handling index `1` is acting)
(b_start == 1) && (b[1] = 0.)

# Assemble matrice and vector
assemble!(A)
assemble!(b)

# Set up the linear solver
ksp = create_ksp(A)
set_from_options!(ksp)
set_up!(ksp)

# Solve the system
x = solve(ksp, b)

# Print the solution (here x is still a `PetscVec`)
@show x

# Convert `PetscVec` to Julia `Array` (and play with it!)
array = vec2array(x)

# Free memory
destroy!(A)
destroy!(b)
destroy!(x)

# Note that it's also possible to build a matrix using the COO format as in `SparseArrays`:
M = create_matrix(3,3; auto_setup = true)
M_start, M_end = get_range(M)
I = [1, 1, 1, 2, 3]
J = [1, 3, 1, 3, 2]
V = [1, 2, 3, 4, 5]
k = findall(x -> M_start <= x <= M_end, I) # just a trick to allow this example to run in parallel
set_values!(M, I[k], J[k], V[k], ADD_VALUES)
assemble!(M)
@show M
destroy!(M)
# This is very convenient in sequential since you can fill the three vectors I, J, V in your code and decide only
# at the last moment if you'd like to use `SparseArrays` or `PetscMat`.

# Finalize Petsc
PetscFinalize()

end #hide


text:
@testset "visualization" begin

    @testset "Peristent plots" begin
        # All tests use the same dir to also catch if one plot overwrites another by default
        testdir = "test_persistent_plots"

        struct DummyGraph
            nv::Int
            np::Int
        end
        NaiveNASlib.nvertices(g::DummyGraph) = g.nv
        NaiveGAflux.nparams(g::DummyGraph) = g.np
        struct PlotTestCand <: AbstractCandidate
            fitness::Real
            graph::DummyGraph
            PlotTestCand(fitness, nv, np) = new(fitness, DummyGraph(nv, np))
        end
        NaiveGAflux.model(c::PlotTestCand, f=identity) = f(c.graph)
        NaiveGAflux.fitness(c::PlotTestCand) = c.fitness

        try
            import MemPool
            @testset "PlotFitness" begin

                p = PlotFitness((args...;kwargs...) -> [], testdir)
                @test !isdir(p.basedir)

                p(PlotTestCand.(1:3, [10, 20, 30], [100, 200, 300]))

                @test p.plt == [1, [3, 2]]
                @test p.best == [3]
                @test p.avg == [2]

                p(PlotTestCand.(2:4, [10, 20, 30], [100, 200, 300]))
                @test p.plt == [1, [3, 2], 2, [4, 3]]
                @test p.best == [3, 4]
                @test p.avg == [2, 3]

                p2 = PlotFitness((args...;kwargs...) -> [args...], testdir)
                @test p2.plt == [[1,2], [3 2; 4 3]]
                @test p2.best == p.best
                @test p2.avg == p.avg
            end

            @testset "ScatterPop" begin

                p = ScatterPop((args...;kwargs...) -> true, testdir)
                @test !isdir(p.basedir)

                @test p(PlotTestCand.(1:3, [10, 20, 30], [100, 200, 300]))
                @test p.data ==  [[10 1 100; 20 2 200; 30 3 300]]

                @test p(PlotTestCand.(2:4, [20, 30, 40], [200, 300, 400]))
                @test p.data ==  [[10 1 100; 20 2 200; 30 3 300], [20 2 200; 30 3 300; 40 4 400]]

                p2 = ScatterPop((args...;kwargs...) -> true, testdir)
                @test p2.data == p.data

                @test p2(PlotTestCand.(3:5, [30, 40, 50], [300, 400, 500]))
                @test p2.data ==  [[10 1 100; 20 2 200; 30 3 300], [20 2 200; 30 3 300; 40 4 400], [30 3 300; 40 4 400; 50 5 500]]
            end

            @testset "ScatterOpt" begin
                NaiveGAflux.opt(c::PlotTestCand) = fitness(c) > 2 ? ADAM(nparams(c) - fitness(c)) : Flux.Optimiser([ShieldedOpt(Descent(nparams(c) - fitness(c)))])

                p = ScatterOpt((args...;kwargs...) -> true, testdir)
                @test !isdir(p.basedir)

                @test p(PlotTestCand.(1:3, [10, 20, 30], [100, 200, 300]))
                @test p.data ==  [[1 99.0 Descent; 2 198.0 Descent; 3 297.0 ADAM]]

                @test p(PlotTestCand.(2:4, [20, 30, 40], [200, 300, 400]))
                @test p.data ==  [[1 99 Descent; 2 198 Descent; 3 297 ADAM], [2 198 Descent; 3 297 ADAM; 4 396 ADAM]]

                p2 = ScatterOpt((args...;kwargs...) -> true, testdir)
                @test p2.data == p.data

                @test p2(PlotTestCand.(3:5, [30, 40, 50], [300, 400, 500]))
                @test p2.data ==  [[1 99 Descent; 2 198 Descent; 3 297 ADAM], [2 198 Descent; 3 297 ADAM; 4 396 ADAM],[3 297 ADAM; 4 396 ADAM; 5 495 ADAM]]

                p3 = ScatterOpt((args...;kwargs...) -> true, testdir)
                @test p3(PlotTestCand.(3:5, [30, 40, 50], [300, 400, 500]))
                # Test that something was added
                @test length(p3.data) == 1 + length(p2.data)
                # What was added is just a copy paste of the last thing inserted to p2 so this lazy check works
                @test p3.data[end] == p2.data[end]

            end

        finally
            rm(testdir, force=true, recursive=true)
            MemPool.cleanup()
        end

        @testset "MultiPlot" begin
            struct MultiPlotTestMock
                nr::Int
                seen::Vector{Int}
                MultiPlotTestMock(nr) = new(nr, Int[])
            end
            NaiveGAflux.plotgen(p::MultiPlotTestMock, gen=length(p.seen)) = gen==0 ? p.nr : p.nr * p.seen[gen]
            function (p::MultiPlotTestMock)(pop)
                push!(p.seen, pop)
                return NaiveGAflux.plotgen(p)
            end

            pdata = []
            spy(x...) = push!(pdata, x...)
            p = MultiPlot(spy, MultiPlotTestMock.(2:4)...)

            @test pdata == 2:4

            @test p(3) == pdata ==[2,3,4,6,9,12]
            @test mapfoldl(pp -> pp.seen, vcat, p.plts) == [3,3,3]
        end

        @testset "CbAll" begin
            cb1 = false
            cb2 = false

            CbAll(pop -> cb1 = pop, pop -> cb2 = pop)(123)

            @test cb1 == 123
            @test cb2 == 123
        end
    end
end


text:
module ThreadsTiming

using ParallelUtilities
using Distributed

function initialize_serial(sleeptime)
    s = zeros(Int, 2_000)
    for ind in eachindex(s)
        sleep(sleeptime)
        s[ind] = ind
    end
    return s
end

function initializenode_threads(sleeptime)
    s = zeros(Int, 2_000)
    Threads.@threads for ind in eachindex(s)
        sleep(sleeptime)
        s[ind] = ind
    end
    return s
end

function main_threads(sleeptime)
    workers_node_pool = ParallelUtilities.workerpool_nodes()
    nw_nodes = nworkers(workers_node_pool)
    pmapreduce(x -> initializenode_threads(sleeptime), hcat, workers_node_pool, 1:nw_nodes)
end

function main_serial(sleeptime)
    workers_node_pool = ParallelUtilities.workerpool_nodes()
    nw_nodes = nworkers(workers_node_pool)
    mapreduce(x -> initialize_serial(sleeptime), hcat, 1:nw_nodes)
end

function compare_with_serial()
    # precompile
    main_serial(0)
    main_threads(0)
    # time
    println("Testing serial")
    A = @time main_serial(5e-3);
    println("Testing threads")
    B = @time main_threads(5e-3);

    println("Results match : ", A == B)
end

end


text:
import UUIDs: UUID
import .ExploreExpression: SymbolsState
import JSON: lower

"The building block of `Notebook`s. Contains code, output and reactivity data."
mutable struct Cell
    "because Cells can be reordered, they get a UUID. The JavaScript frontend indexes cells using the UUID."
    uuid::UUID
    code::String
    parsedcode::Any
    output_repr::Union{String, Nothing}
    error_repr::Union{String, Nothing}
    repr_mime::MIME
    runtime::Union{Missing,UInt64}
    symstate::SymbolsState
    module_usings::Set{Expr}
end

Cell(uuid, code) = Cell(uuid, code, nothing, nothing, nothing, MIME("text/plain"), missing, SymbolsState(), Set{Expr}())
Cell(code) = Cell(uuid1(), code)


text:
export solve_qp, solve_qp_batch

using Convex
using Zygote: @adjoint
using SCS
using LinearAlgebra

"""
Solving QP

Given `Q`, `p`, `G`, `h`, `A` and `b`, it solves following convex optimization problem:

```math
min 0.5*x'Qx + p'x
s.t Gx ≤ h
    Ax = b
```

# Parameters

- `Q`: array of size (`dim_x`, `dim_x`)
- `p`: array of size (`dim_x`,)
- `G`: array of size (`num_ineq`, `dim_x`), can be empty
- `h`: array of size (`num_ineq`,), can be empty
- `A`: array of size (`num_eq`, `dim_x`), can be empty
- `b`: array of size (`num_eq`,), can be empty

# Return

- `x`: array of size (`dim_x`,) which is the solution to the given QP
- `λs_`: adjoint values of inequality constraints
- `νs_`: adjoint values of equality constraints
"""
function solve_qp(Q, p, G, h, A, b, optimizer = SCS.Optimizer; kwargs...)
    !(isempty(G) ⊻ isempty(h)) ||
        throw(ArgumentError("G and h should be both empty or not empty"))
    !(isempty(A) ⊻ isempty(b)) ||
        throw(ArgumentError("A and b should be both empty or not empty"))
    # solve qp and return x
    dim_x = size(Q, 1)
    x = Variable(dim_x)
    constraints = Dict{String,Constraint}()
    if !isempty(A)
        constraints["Eq"] = (A * x == b)
    end
    if !isempty(G)
        constraints["InEq"] = (G * x ≤ h)
    end
    obj = 0.5 * quadform(x, Q) + p'x
    problem = minimize(obj, collect(values(constraints)))
    if Sys.isunix()
        stream_null = open("/dev/null", "w")
    else
        # assume windows
        stream_null = open("nul", "w")
    end
    redirect_stdout(stream_null) do
        redirect_stderr(stream_null) do
            solve!(problem, optimizer; kwargs...)
        end
    end
    close(stream_null)
    @assert problem.status == Convex.MOI.OPTIMAL
    x_ = isa(x.value, AbstractArray) ? dropdims(x.value, dims = 2) : [x.value]
    if haskey(constraints, "InEq")
        cond_ineq = constraints["InEq"]
        λs_ = isa(cond_ineq.dual, AbstractArray) ? dropdims(cond_ineq.dual, dims = 2) :
            [cond_ineq.dual]
    else
        λs_ = similar(G, 0)
    end
    if haskey(constraints, "Eq")
        cond_eq = constraints["Eq"]
        νs_ = isa(cond_eq.dual, AbstractArray) ? dropdims(abs.(cond_eq.dual), dims = 2) :
            [abs(cond_eq.dual)]
    else
        νs_ = similar(A, 0)
    end
    return x_, λs_, νs_
end

"""
Solving QP in batch

# Parameters

- `Q`: array of size (`dim_x`, `dim_x`, `num_batch`)
- `p`: array of size (`dim_x`, `num_batch`)
- `G`: array of size (`num_ineq`, `dim_x`, `num_batch`)
- `h`: array of size (`num_ineq`, `num_batch`)
- `A`: array of size (`num_eq`, `dim_x`, `num_batch`)
- `b`: array of size (`num_eq`, `num_batch`)

# Return

- `X`: array of size (`dim_x`, `num_batch`) which are the solutions to the given QP batch
- `λs`: array of size (`dim_ineq`, `num_batch`) which are the adjoint values of inequality constraints
- `νs`: array of size (`dim_eq`, `num_batch`) which are the adjoint values of equality constraints
"""
function solve_qp_batch(
    Q::AbstractArray{T,3},
    p::AbstractArray{T,2},
    G::AbstractArray{T,3},
    h::AbstractArray{T,2},
    A::AbstractArray{T,3},
    b::AbstractArray{T,2},
    optimizer = SCS.Optimizer,
) where {T<:Real}
    num_batch = size(Q, 3)
    dim_x = size(Q, 1)
    num_eq = size(A, 1)
    num_ineq = size(G, 1)
    X = similar(p)
    λs = similar(h)
    νs = similar(b)
    @inbounds for i = 1:num_batch
        Q_ = @view Q[:, :, i]
        G_ = @view G[:, :, i]
        p_ = @view p[:, i]
        h_ = @view h[:, i]
        A_ = @view A[:, :, i]
        b_ = @view b[:, i]
        x, λs_, νs_ = solve_qp(Q_, p_, G_, h_, A_, b_, optimizer)
        X[:, i] .= x
        if !isempty(λs_)
            λs[:, i] .= λs_
        end
        if !isempty(νs_)
            νs[:, i] .= νs_
        end
    end
    return X, λs, νs
end

@adjoint function solve_qp_batch(
    Q::AbstractArray{T,3},
    p::AbstractArray{T,2},
    G::AbstractArray{T,3},
    h::AbstractArray{T,2},
    A::AbstractArray{T,3},
    b::AbstractArray{T,2},
    optimizer = SCS.Optimizer,
) where {T<:Real}
    return solve_qp_batch_with_back(Q, p, G, h, A, b, optimizer)
end

function solve_qp_batch_with_back(
    Q::AbstractArray{T,3},
    p::AbstractArray{T,2},
    G::AbstractArray{T,3},
    h::AbstractArray{T,2},
    A::AbstractArray{T,3},
    b::AbstractArray{T,2},
    optimizer = SCS.Optimizer,
) where {T<:Real}
    X, λs, νs = solve_qp_batch(Q, p, G, h, A, b, optimizer)
    num_batch = size(X, 2)
    function backward(Δ)
        ΔX, Δλs, Δνs = Δ
        ΔQ = similar(Q)
        Δp = similar(p)
        ΔG = similar(G)
        Δh = similar(h)
        ΔA = similar(A)
        Δb = similar(b)
        # solve dQ, dp, dG, dh, dA, db
        @inbounds for i = 1:num_batch
            ν_ = @view νs[:, i]
            Q_ = @view Q[:, :, i]
            G_ = @view G[:, :, i]
            h_ = @view h[:, i]
            A_ = @view A[:, :, i]
            ΔX_ = @view ΔX[:, i]
            λ_ = @view λs[:, i]
            x_ = @view X[:, i]
            dx, dλ, dν = solve_kkt(Q_, G_, h_, A_, λ_, x_, ΔX_)
            ΔQ[:, :, i] .= 0.5 * (dx * x_' + x_ * dx') # dQ
            Δp[:, i] .= dx # dp
            ΔG[:, :, i] .= Diagonal(λ_) * (dλ * x_' + λ_ * dx') # dG
            Δh[:, i] .= -Diagonal(λ_) * dλ # dh
            ΔA[:, :, i] .= (dν * x_' + ν_ * dx') # dA
            Δb[:, i] .= -dν # db
        end
        Δoptimizer = nothing
        (ΔQ, Δp, ΔG, Δh, ΔA, Δb, Δoptimizer)
    end
    return (X, λs, νs), backward
end

function solve_kkt(Q, G, h, A, λ, x, Δx)
    # TODO: G, h or A are empty
    # TODO: implement primal-dual interior point method (PDIPM) for solving kkt, better GPU utilization
    M = kkt_matrix(Q, G, h, A, λ, x)
    ndim = size(M, 1)
    Δl = [Δx; zeros(ndim - size(Δx, 1))]
    sol = pinv(-M) * Δl
    dx = sol[1:length(x)]
    dλ = sol[length(x)+1:length(x)+length(λ)]
    dν = sol[length(x)+length(λ)+1:end]
    return dx, dλ, dν
end

function kkt_matrix(Q, G, h, A, λ, x)
    # TODO: G, h or A are empty
    return [
        Q G' * Diagonal(λ) A'
        G Diagonal(G * x - h) zeros(size(G, 1), size(A, 1))
        A zeros(size(A, 1), size(h, 1)) zeros(size(A, 1), size(A, 1))
    ]
end


text:
import Base: +, -, *, /, copy!

for op in (:*, :/)
    @eval begin
        function $op(a::T, s::Spectrum) where T <: Number
            t = deepcopy(s)
            broadcast!($op, t.fid, t.fid, a)
            for (_, p) in t.procs
                broadcast!($op, p.re_ft, p.re_ft, a)
                broadcast!($op, p.im_ft, p.im_ft, a)
            end
            t
        end
    end
end

for op in (:+, :-, :*)
    @eval begin
        function $op(args::Vararg{Spectrum, N}) where N
            shifts = union_shifts(collect(args))
            vals = broadcast($op, (resample(s, shifts) for s in args)...)
            # TODO: make a proper processed Spectrum
            (shifts, vals)
        end
    end
end

function copy!(src::Spectrum, dst::Spectrum, shift_rng::Tuple{Float64,Float64})
    d = resample(src, chemical_shifts(dst))
    rng = ppmtoindex(dst, shift_rng)
    dst[rng] = d[rng]
    dst
end