text:
using Test
using BioPredictions
using BioSequences

seq = LongAminoAcidSeq([aa for aa in alphabet(AminoAcid)])

@testset "Gravy index" begin
    @test gravy(seq) == -0.35
end

@testset "Molecular weight" begin
    @test isapprox(weight(seq),
                   weight(seq; ignoreunknown = true) + sum(map(aa -> BioPredictions.aminoacidweights[aa], [AA_X, AA_Z, AA_B, AA_J])))
end

@testset "Isoelectric point" begin
    @test pI(seq) == 6.6
    @test pI(seq, "Solomon") == 7.14
    @test pI(seq, collect(1:9)) == 4.5
end

@testset "Cleavage sites" begin
    @test cleave(seq[1:20]) == (positions = [2, 12], peptides = LongAminoAcidSeq[aa"AR", aa"NDCQEGHILK", aa"MFPSTWYV"])
    @test cleave(seq[1:20], BioPredictions.proteinaseK) == (positions = [1, 7, 10, 11, 14, 17, 18, 19], peptides = LongAminoAcidSeq[aa"A", aa"RNDCQE", aa"GHI", aa"L", aa"KMF", aa"PST", aa"W", aa"Y", aa"V"])
    @test cleave(seq[1:20], BioPredictions.lysC) == (positions = [12], peptides = LongAminoAcidSeq[aa"ARNDCQEGHILK", aa"MFPSTWYV"])
    @test cleave(seq[1:20], BioPredictions.proline_endopeptidase) == (positions = Int[], peptides = LongAminoAcidSeq[seq[1:20]])
end


text:
using Dash, DashVtk, DashHtmlComponents
using Random

Random.seed!(42)

app = dash()

volume_view = vtk_view(
   children=vtk_volumedatarepresentation(
     spacing=[1, 1, 1],
     dimensions=[10, 10, 10],
     origin=[0, 0, 0],
     scalars = [rand() for _ in 1:1000],
     rescaleColorMap=false,
   )
)

app.layout = html_div() do
    html_div(
    style = Dict(
          "height" => "calc(100vh - 16px)",
          ),
   children = [html_div(
          children=volume_view,
          style = Dict(
          "height" => "100%",
          "width" =>  "100%",
        ),
      )]
    )
end

run_server(app, "0.0.0.0", debug = true)

ext: jl

max_stars_repo_path: src/Randomely Generated Volume/app.jl

max_issues_repo_path: src/Randomely Generated Volume/app.jl

max_forks_repo_path: src/Randomely Generated Volume/app.jl


text:
#
# Numeric vector sort
#
using Printf, Random

function qsort!(a, lo, hi)
    i, j = lo, hi
    while i < hi
        pivot = a[(lo + hi) >>> 1]
        while i <= j
            while a[i] < pivot
                i += 1
            end
            while a[j] > pivot
                j -= 1
            end
            if i <= j
                a[i], a[j] = a[j], a[i]
                i, j = i+1, j-1
            end
        end
        if lo < j
            qsort!(a, lo, j)
        end
        lo, j = i, hi
    end
    return a
end

# seed rng for more consistent timings
Random.seed!(1776)

n = 10000000
data = rand(n)
t0 = time_ns()
sorted = qsort!(data, 1, n)
t1 = time_ns()

@printf("size = %d, msec = %.1f\n", n, (t1 - t0) / 1e6)

max_stars_repo_name: sergev/vak-opensource

max_issues_repo_name: sergev/vak-opensource

max_forks_repo_name: sergev/vak-opensource


text:
# Viz environment
using Pkg
Pkg.add("Plots")

max_stars_repo_name: charleskawczynski/PkgTemplates

max_issues_repo_name: charleskawczynski/PkgTemplates

max_forks_repo_name: charleskawczynski/PkgTemplates


text:
module ClimateMachine

using Pkg.TOML

const CLIMATEMACHINE_VERSION =
    VersionNumber(TOML.parsefile(joinpath(dirname(@__DIR__), "Project.toml"))["version"])

include(joinpath("Utilities", "TicToc", "TicToc.jl"))
include(joinpath("InputOutput", "Writers", "Writers.jl"))
include(joinpath("Driver", "ConfigTypes", "ConfigTypes.jl"))
include(joinpath("Utilities", "VariableTemplates", "VariableTemplates.jl"))
include(joinpath("Arrays", "MPIStateArrays.jl"))
include(joinpath("Numerics", "Mesh", "Mesh.jl"))
include(joinpath("Common", "CartesianDomains", "CartesianDomains.jl"))
include(joinpath("Common", "CartesianFields", "CartesianFields.jl"))
include(joinpath("Numerics", "DGMethods", "Courant.jl"))
include(joinpath("BalanceLaws", "Problems.jl"))
include(joinpath("BalanceLaws", "BalanceLaws.jl"))
include(joinpath("Numerics", "DGMethods", "DGMethods.jl"))
include(joinpath("Common", "Orientations", "Orientations.jl"))
include(joinpath("Utilities", "SingleStackUtils", "SingleStackUtils.jl"))
include(joinpath("Numerics", "SystemSolvers", "SystemSolvers.jl"))
include(joinpath("Numerics", "ODESolvers", "GenericCallbacks.jl"))
include(joinpath("Numerics", "ODESolvers", "ODESolvers.jl"))
include(joinpath("Land", "Model", "LandModel.jl"))
include(joinpath("InputOutput", "VTK", "VTK.jl"))
include(joinpath("Common", "TurbulenceClosures", "TurbulenceClosures.jl"))
include(joinpath("Common", "TurbulenceConvection", "TurbulenceConvection.jl"))
include(joinpath("Atmos", "Model", "AtmosModel.jl"))
include(joinpath("Common", "Spectra", "Spectra.jl"))
include(joinpath("Diagnostics", "Diagnostics.jl"))
include(joinpath("Diagnostics", "DiagnosticsMachine", "DiagnosticsMachine.jl"))
include(joinpath("Diagnostics", "StdDiagnostics", "StdDiagnostics.jl"))
include(joinpath("Diagnostics", "Debug", "StateCheck.jl"))
include(joinpath("Driver", "Checkpoint", "Checkpoint.jl"))
include(joinpath("Driver", "Callbacks", "Callbacks.jl"))
include(joinpath("Driver", "Driver.jl"))

include(joinpath("Ocean", "Ocean.jl"))

end


text:
# Interface
# single input
function forward(op::Linear{T, Matrix{T}, Vector{T}, Vector{T}}, input::Vector{T}) where T
    op.input = input # NOTE: this is only a reference, no extra memmory alloc here
    BLAS.gemv('N', op.weight, input) + op.bias
end

# batched input
function forward(op::Linear{T, Matrix{T}, Vector{T}, Matrix{T}}, input::Matrix{T}) where T
    # NOTE: we use colum as a sample of data
    # PyTorch use the first dimension
    op.input = input
    return BLAS.gemm('N', 'N', op.weight, input) .+ op.bias # Be careful, here has a broadcast .+
end

# static array (immutable & mutable arrays)
# single input
function forward(op::Linear{T, M, B, I}, input::I) where {
        T,
        M<:StaticMatrix,
        B<:StaticVector,
        I<:StaticVector
    }

    op.input = input
    return op.weight * input + op.bias
end

# batched input
function forward(op::Linear{T, M, B, I}, input::I) where {
        T,
        M<:StaticMatrix,
        B<:StaticVector,
        I<:StaticMatrix
    }

    op.input = input
    return op.weight * input .+ op.bias
end


text:
@testset "Utilities" begin
  a, b, c = P2(0,0), P2(1,0), P2(0,1)
  @test signarea(a, b, c) == T(0.5)
  a, b, c = P2(0,0), P2(0,1), P2(1,0)
  @test signarea(a, b, c) == T(-0.5)

  p1, p2, p3 = P2(0,0), P2(1,1), P2(0.25,0.5)
  s = Segment(P2(0.5,0.0), P2(0.0,1.0))
  @test sideof(p1, s) == :LEFT
  @test sideof(p2, s) == :RIGHT
  @test sideof(p3, s) == :ON
end


text:
# The following code is a modified version of
#
#    function _threadsfor(iter, lbody, schedule)
#    macro threads(args...)
#
# in https://github.com/JuliaLang/julia/blob/9f3265399227fbfc4f0160ec3592a9262bd3eb5f/base/threadingconstructs.jl
#
# Its license is MIT: https://julialang.org/license

using Base.Threads
using Base.Threads: threading_run

function _my_threadsfor(iter, lbody, prebody, postbody, schedule)
    lidx = iter.args[1]         # index
    range = iter.args[2]
    quote
        local threadsfor_fun
        let range = $(esc(range))
        function threadsfor_fun(onethread=false)
            r = range # Load into local variable
            lenr = length(r)
            # divide loop iterations among threads
            if onethread
                tid = 1
                len, rem = lenr, 0
            else
                tid = threadid()
                len, rem = divrem(lenr, nthreads())
            end
            # not enough iterations for all the threads?
            if len == 0
                if tid > rem
                    return
                end
                len, rem = 1, 0
            end
            # compute this thread's iterations
            f = firstindex(r) + ((tid-1) * len)
            l = f + len - 1
            # distribute remaining iterations evenly
            if rem > 0
                if tid <= rem
                    f = f + (tid-1)
                    l = l + tid
                else
                    f = f + rem
                    l = l + rem
                end
            end
            # run prebody
            $(esc(prebody))
            # run this thread's iterations
            for i = f:l
                local $(esc(lidx)) = @inbounds r[i]
                $(esc(lbody))
            end
            # run postbody
            $(esc(postbody))
        end
        end
        if threadid() != 1 || ccall(:jl_in_threaded_region, Cint, ()) != 0
            $(if schedule === :static
              :(error("`@my_threads :static` can only be used from thread 1 and not nested"))
              else
              # only use threads when called from thread 1, outside @threads
              :(Base.invokelatest(threadsfor_fun, true))
              end)
        else
            threading_run(threadsfor_fun)
        end
        nothing
    end
end

"""
    @my_threads

A macro to parallelize a `for` loop to run with multiple threads.
It splits the iteration space among multiple tasks with `prebody` and `postbody`.
It runs those tasks on threads according to a scheduling policy.

Usage:

```julia
@my_threads [schedule] begin
    prebody
end for ...
    ...
end begin
    postbody
end
```
"""
macro my_threads(args...)
    na = length(args)
    if na == 4
        sched, prebody, ex, bostbody = args
        if sched isa QuoteNode
            sched = sched.value
        elseif sched isa Symbol
            # for now only allow quoted symbols
            sched = nothing
        end
        if sched !== :static
            throw(ArgumentError("unsupported schedule argument in @threads"))
        end
    elseif na == 3
        sched = :default
        prebody, ex, postbody = args
    else
        throw(ArgumentError("wrong number of arguments in @my_threads"))
    end
    if !(isa(ex, Expr) && ex.head === :for)
        throw(ArgumentError("@my_threads requires a `for` loop expression"))
    end
    if !(ex.args[1] isa Expr && ex.args[1].head === :(=))
        throw(ArgumentError("nested outer loops are not currently supported by @my_threads"))
    end
    return _my_threadsfor(ex.args[1], ex.args[2], prebody, postbody, sched)
end


text:
@testset "ch04_sim02" begin
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    world = World(xlim, ylim)
    circlings = Array{Agent,1}(undef, 0)
    robots = Array{RealRobot,1}(undef, 0)
    for i = 1:10
        circling = Agent(0.2, 10.0 / 180 * pi)
        robot =
            RealRobot([0.0, 0.0, 0.0], circling, nothing; radius = 0.05, color = "black")
        push!(circlings, circling)
        push!(robots, robot)
        push!(world, robot)
    end

    dt = 0.1
    anim = @animate for i = 1:300
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        p = draw(world, annota)
        for j = 1:10
            obsv = observations(robots[j].sensor_, robots[j].pose_)
            @assert obsv == nothing
            v, ω = decision(circlings[j], obsv)
            state_transition(robots[j], v, ω, dt; move_noise = true)
        end
    end
    if GUI
        gif(anim, "ch04_sim02.gif", fps = 20)
    end
end

@testset "ch04_sim03" begin
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    world = World(xlim, ylim)
    circling_agent = Agent(0.2, 10.0 / 180 * pi)
    nobias_robot = IdealRobot([0.0, 0.0, 0.0], circling_agent, nothing, 0.05, "gray")
    push!(world, nobias_robot)
    biased_robot = RealRobot(
        [0.0, 0.0, 0.0],
        circling_agent,
        nothing;
        radius = 0.05,
        color = "red",
        bias_rate_stds = (0.2, 0.2),
    )
    push!(world, biased_robot)
    dt = 0.1
    anim = @animate for i = 1:300
        annota = "t = $(round(dt * i, sigdigits=3))[s]"
        obsv1 = observations(nobias_robot.sensor_, nobias_robot.pose_)
        v1, ω1 = decision(circling_agent, obsv1)
        obsv2 = observations(biased_robot.sensor_, biased_robot.pose_)
        v2, ω2 = decision(circling_agent, obsv2)
        p = draw(world, annota)
        state_transition(nobias_robot, v1, ω1, dt)
        state_transition(biased_robot, v2, ω2, dt; vel_bias_noise = true)
    end
    if GUI
        gif(anim, "ch04_sim03.gif", fps = 20)
    end
end

@testset "ch04_sim04" begin
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    world = World(xlim, ylim)

    circling_agent = Agent(0.2, 10.0 / 180 * pi)
    robots = Array{RealRobot,1}(undef, 0)
    for i = 1:10
        robot = RealRobot(
            [0.0, 0.0, 0.0],
            circling_agent,
            nothing;
            radius = 0.05,
            color = "gray",
            expected_stuck_time = 60.0,
            expected_escape_time = 60.0,
        )
        push!(robots, robot)
        push!(world, robot)
    end
    ideal_robot = IdealRobot([0.0, 0.0, 0.0], circling_agent, nothing, 0.05, "red")
    push!(world, ideal_robot)

    dt = 0.1
    anim = @animate for i = 1:300
        annota = "t = $(round(dt * i, sigdigits=3))[s]"
        zs = []
        for j = 1:10
            obsv = observations(robots[j].sensor_, robots[j].pose_)
            if obsv == nothing
                push!(zs, nothing)
            else
                push!(zs, copy(obsv))
            end
        end
        obsv1 = observations(ideal_robot.sensor_, ideal_robot.pose_)
        p = draw(world, annota)
        for j = 1:10
            v, ω = decision(circling_agent, zs[j])
            state_transition(robots[j], v, ω, dt; stuck_noise = true)
        end
        v, ω = decision(circling_agent, obsv1)
        state_transition(ideal_robot, v, ω, dt)
    end
    if GUI
        gif(anim, "ch04_sim04.gif", fps = 20)
    end
end

@testset "ch04_sim05" begin
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    world = World(xlim, ylim)

    circling_agent = Agent(0.2, 10.0 / 180 * pi)
    robots = Array{RealRobot,1}(undef, 0)
    for i = 1:10
        robot = RealRobot(
            [0.0, 0.0, 0.0],
            circling_agent,
            nothing;
            radius = 0.05,
            color = "gray",
            expected_kidnap_time = 5.0,
        )
        push!(robots, robot)
        push!(world, robot)
    end
    ideal_robot = IdealRobot([0.0, 0.0, 0.0], circling_agent, nothing, 0.05, "red")
    push!(world, ideal_robot)

    dt = 0.1
    anim = @animate for i = 1:300
        annota = "t = $(round(dt * i, sigdigits=3))[s]"
        zs = []
        for j = 1:10
            obsv = observations(robots[j].sensor_, robots[j].pose_)
            if obsv == nothing
                push!(zs, nothing)
            else
                push!(zs, copy(obsv))
            end
        end
        obsv1 = observations(ideal_robot.sensor_, ideal_robot.pose_)
        p = draw(world, annota)
        for j = 1:10
            v, ω = decision(circling_agent, zs[j])
            state_transition(robots[j], v, ω, dt; kidnap = true)
        end
        v, ω = decision(circling_agent, obsv1)
        state_transition(ideal_robot, v, ω, dt)
    end
    if GUI
        gif(anim, "ch04_sim05.gif", fps = 20)
    end
end

@testset "ch04_sim07" begin
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    world = World(xlim, ylim)
    circling_agent = Agent(0.2, 10.0 / 180 * pi)
    landmarks =
        [Landmark([-4.0, 2.0], 0), Landmark([2.0, -3.0], 1), Landmark([3.0, 3.0], 1)]
    m = Map()
    push!(m, landmarks)
    robot =
        RealRobot([0.0, 0.0, pi / 6], circling_agent, RealCamera(landmarks); color = "red")
    push!(world, robot)
    push!(world, m)
    dt = 0.1
    anim = @animate for i = 1:300
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        obsv = observations(robot.sensor_, robot.pose_; noise = true)
        p = draw(world, annota)
        v, ω = decision(circling_agent, obsv)
        state_transition(robot, v, ω, dt)
    end
    if GUI
        gif(anim, "ch04_sim07.gif", fps = 20)
    end
end

@testset "ch04_sim08" begin
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    world = World(xlim, ylim)
    straight_agent = Agent(0.2, 0.0)
    landmarks = [
        Landmark([-4.0, 2.0], 0),
        Landmark([3.0, -3.0], 1),
        Landmark([3.0, 3.0], 2),
        Landmark([3.0, -2.0], 3),
        Landmark([3.0, 0.0], 4),
        Landmark([3.0, 1.0], 5),
    ]
    m = Map()
    push!(m, landmarks)
    robot = RealRobot([0.0, 0.0, 0.0], straight_agent, RealCamera(landmarks); color = "red")
    push!(world, robot)
    push!(world, m)
    dt = 0.1
    anim = @animate for i = 1:300
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        obsv = observations(robot.sensor_, robot.pose_; noise = true, bias = true)
        p = draw(world, annota)
        v, ω = decision(straight_agent, obsv)
        state_transition(robot, v, ω, dt)
    end
    if GUI
        gif(anim, "ch04_sim08.gif", fps = 20)
    end
end

@testset "ch04_sim09" begin
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    world = World(xlim, ylim)
    landmarks =
        [Landmark([-4.0, 2.0], 0), Landmark([2.0, -3.0], 1), Landmark([3.0, 3.0], 2)]
    m = Map()
    push!(m, landmarks)
    circling = Agent(0.2, 10.0 / 180 * pi)
    robot = RealRobot(
        [0.0, 0.0, 0.0],
        circling,
        RealCamera(landmarks; phantom_prob = 0.2);
        color = "red",
    )
    push!(world, robot)
    push!(world, m)
    dt = 0.1
    anim = @animate for i = 1:300
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        z = observations(
            robot.sensor_,
            robot.pose_;
            noise = true,
            bias = true,
            phantom = true,
        )
        p = draw(world, annota)
        v, ω = decision(circling, z)
        state_transition(robot, v, ω, dt; move_noise = true, vel_bias_noise = true)
    end
    if GUI
        gif(anim, "ch04_sim09.gif", fps = 20)
    end
end

@testset "ch04_sim10" begin
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    world = World(xlim, ylim)
    landmarks =
        [Landmark([-4.0, 2.0], 0), Landmark([2.0, -3.0], 1), Landmark([3.0, 3.0], 2)]
    m = Map()
    push!(m, landmarks)
    circling = Agent(0.2, 10.0 / 180 * pi)
    robot = RealRobot(
        [0.0, 0.0, 0.0],
        circling,
        RealCamera(landmarks; overlook_prob = 0.5);
        color = "red",
    )
    push!(world, robot)
    push!(world, m)
    dt = 0.1
    anim = @animate for i = 1:300
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        z = observations(
            robot.sensor_,
            robot.pose_;
            noise = true,
            bias = true,
            overlook = true,
        )
        p = draw(world, annota)
        v, ω = decision(circling, z)
        state_transition(robot, v, ω, dt; move_noise = true, vel_bias_noise = true)
    end
    if GUI
        gif(anim, "ch04_sim10.gif", fps = 20)
    end
end

@testset "ch04_sim11" begin
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    world = World(xlim, ylim)
    landmarks =
        [Landmark([-4.0, 2.0], 0), Landmark([2.0, -3.0], 1), Landmark([3.0, 3.0], 2)]
    m = Map()
    push!(m, landmarks)
    circling = Agent(0.2, 10.0 / 180 * pi)
    robot = RealRobot(
        [2.0, 2.0, pi / 6],
        circling,
        RealCamera(landmarks; occlusion_prob = 0.1);
        color = "red",
    )
    push!(world, robot)
    push!(world, m)
    dt = 0.1
    anim = @animate for i = 1:300
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        z = observations(
            robot.sensor_,
            robot.pose_;
            noise = true,
            bias = true,
            occlusion = true,
        )
        p = draw(world, annota)
        v, ω = decision(circling, z)
        state_transition(robot, v, ω, dt; move_noise = true, vel_bias_noise = true)
    end
    if GUI
        gif(anim, "ch04_sim11.gif", fps = 20)
    end
end

max_stars_repo_name: soblin/JuliaProbo

max_issues_repo_name: soblin/JuliaProbo

max_forks_repo_name: soblin/JuliaProbo


text:
@inline function tnorm2(x::NTuple{N,T},type::Type{T}) where {N,T}
    x[1]*x[1]+tnorm2(Base.tail(x),T)
end

@inline function tnorm2(x::NTuple{0},type::Type{T}) where {T}
    zero(T)
end

@inline function tnorm2(x::NTuple{N,T}) where {N,T}
    tnorm2(x,T)
end

function tnorm2(x)
    tnorm2(promote(x...))
end

function select_seeds(ima,P)
    S=size(ima)
    l=@. S>>1
    h=@. (S+l)>>1
    seeds=Tuple{ntuple(i->Int,length(S))...}[]
    for p in P
        all(l.<=p.<=h) || continue;
        push!(seeds,p)
    end
    seeds
end

function integrate_disk(ima,P,r)
    integrate_ring(ima,P,oftype(r,0),r)
end

function integrate_ring(ima,P,rs,rh=rs*sqrt(2))
    if rs>rh
        rs,rh=rh,rs
    end
    r_half=rh+1
    op_size=ntuple(i->round(Int,rh+1),ndims(ima))
    kop=KernelOp(ima,op_size) do A,Is,ci_I
        Iim=zero(eltype(A))
        I=Tuple(ci_I)
        N=0
        for ci_Itrasl in Is
            Itrasl=Tuple(ci_Itrasl)
            dist=tnorm2(Itrasl.-I)
            if rs*rs <= dist <= rh*rh
                Iim += A[ci_Itrasl]
                N+=1
            end
        end
        Iim,N
    end
    [kop[round.(Int,p)...] for p in P]
end

function scan_radiuses(ima,P,rs::Int)
    min_r=max(rs>>2,1)
    max_r=(rs<<2)+min_r
    T=eltype(ima)
    r_range=min_r:max_r
    sums=zeros(T,length(r_range))
    i=1
    for r in r_range
        sums[i] = sum(x->x[1],integrate_ring(ima,P,r))
        i+=1
    end
    r_range,sums/length(P)
end

function kern_find_best_radius(ima,P,rs::Int)
    r_range,intensity=scan_radiuses(ima,P,rs)
    r_pos=rs-minimum(r_range)+1
    Δ=intensity[r_pos]*0.2 # imagine an error of 20% on the intensity
    Ns=similar(intensity)
    for i in eachindex(intensity)
        I=intensity[i]
        N=sum(@. I-Δ <= intensity <= I+Δ )
        Ns[i]=N
    end
    Nmax_greedy=maximum(Ns)
    mask=(Ns).>=(Nmax_greedy-1) # there is a noise in Ns hence the -1
    best_radiuses=r_range[mask]
    round(Int,sum(best_radiuses)/length(best_radiuses))
end

function find_best_radius(ima,rawP,rs::Int)
    P=select_seeds(ima,rawP)
    r=rs
    for i in 1:10
        nr=kern_find_best_radius(ima,P,r)
        nr==r && break
        r=nr
    end
    r
end


text:
# Euclidean distance from a set

export DistL2

"""
**Distance from a convex set**

    DistL2(ind_S)

Given `ind_S` the indicator function of a convex set ``S``, and an optional positive parameter `λ`, returns the (weighted) Euclidean distance from ``S``, that is function
```math
g(x) = λ\\mathrm{dist}_S(x) = \\min \\{ λ\\|y - x\\| : y \\in S \\}.
```
"""
struct DistL2{R <: Real, T <: ProximableFunction} <: ProximableFunction
    ind::T
    lambda::R
    function DistL2{R, T}(ind::T, lambda::R) where {R <: Real, T <: ProximableFunction}
        if !is_set(ind) || !is_convex(ind)
            error("`ind` must be a convex set")
        end
        if lambda < 0
            error("parameter `λ` must be nonnegative")
        else
            new(ind, lambda)
        end
    end
end

is_prox_accurate(f::DistL2) = is_prox_accurate(f.ind)
is_convex(f::DistL2) = is_convex(f.ind)

DistL2(ind::T, lambda::R=1) where {R <: Real, T <: ProximableFunction} = DistL2{R, T}(ind, lambda)

function (f::DistL2)(x::AbstractArray{T}) where {R <: Real, T <: RealOrComplex{R}}
    p, = prox(f.ind, x)
    return f.lambda * normdiff(x, p)
end

function prox!(y::AbstractArray{T}, f::DistL2, x::AbstractArray{T}, gamma::R=R(1)) where {R <: Real, T <: RealOrComplex{R}}
    prox!(y, f.ind, x)
    d = normdiff(x, y)
    gamlam = gamma * f.lambda
    if gamlam < d
        gamlamd = gamlam/d
        y .= (1 - gamlamd) .*x .+ gamlamd .* y
        return f.lambda * (d - gamlam)
    end
    return R(0)
end

function gradient!(y::AbstractArray{T}, f::DistL2, x::AbstractArray{T}) where {R <: Real, T <: RealOrComplex{R}}
    prox!(y, f.ind, x) # Use y as temporary storage
    dist = normdiff(x, y)
    if dist > 0
        y .= (f.lambda / dist) .* (x .- y)
    else
        y .= 0
    end
    return f.lambda * dist
end

fun_name(f::DistL2) = "Euclidean distance from a convex set"
fun_dom(f::DistL2) = fun_dom(f.ind)
fun_expr(f::DistL2) = "x ↦ λ inf { ||x-y|| : y ∈ S} "
fun_params(f::DistL2) = string("λ = $(f.lambda), S = ", typeof(f.ind))

function prox_naive(f::DistL2, x::AbstractArray{T}, gamma::R=R(1)) where {R <: Real, T <: RealOrComplex{R}}
    p, = prox(f.ind, x)
    d = norm(x - p)
    gamlam = gamma * f.lambda
    if d > gamlam
        return x + gamlam/d * (p - x), f.lambda * (d - gamlam)
    end
    return p, R(0)
end


text:
using OhMyREPL
using BenchmarkTools
using Pkg

enable_autocomplete_brackets(false)

"""
Display the entire type hierarchy starting from the specified `roottype`.
"""
function subtypetree(roottype, level=1, indent=4)
    level==1 && println(roottype)
    for s in subtypes(roottype)
        println(join(fill(" ", level*indent))*string(s))
        subtypetree(s, level+1, indent)
    end
end # function

"""
Get the subtypes of a Union.
"""
Base.collect(t::Union{Type, DataType, Union{}}) = __startup__collect(t, [])
__startup__collect(t::Type, list) = t<:Union{} ? push!(list, t) :
    __startup__collect(t.b, push!(list, t.a))
__startup__collect(t::Union{DataType,Core.TypeofBottom}, list) = push!(list, t)

# Activate package if REPL started in a package folder.
if isfile("Project.toml") && isfile("Manifest.toml")
    Pkg.activate(".")
elseif isfile("../Project.toml") && isfile("../Manifest.toml")
    Pkg.activate("..")
end

max_stars_repo_name: danyloM/dotfiles

max_issues_repo_name: dmalyuta/dotfiles

max_forks_repo_name: dmalyuta/dotfiles


text:
include("../utils.jl")
include("group_utils.jl")

using Manifolds: connection

@testset "Cartan-Schouten connections" begin
    SO3 = SpecialOrthogonal(3)
    SO3minus = ConnectionManifold(SO3, CartanSchoutenMinus())
    SO3plus = ConnectionManifold(SO3, CartanSchoutenPlus())
    SO3zero = ConnectionManifold(SO3, CartanSchoutenZero())

    e = Matrix{Float64}(I, 3, 3)
    p = exp(hat(SO3, Identity(SO3), [1.0, 2.0, 3.0]))
    q = exp(hat(SO3, Identity(SO3), [3.0, 4.0, 1.0]))
    X = hat(SO3, Identity(SO3), [2.0, 3.0, 4.0])
    SO3e = Identity(SO3)

    @testset "connection" begin
        @test connection(SO3minus) === CartanSchoutenMinus()
        @test connection(SO3plus) === CartanSchoutenPlus()
        @test connection(SO3zero) === CartanSchoutenZero()
    end

    @testset "log/exp" begin
        for CSO3 in [SO3minus, SO3plus, SO3zero]
            @test isapprox(SO3, exp(CSO3, p, X), exp(SO3, p, X))
            @test isapprox(SO3, p, log(CSO3, p, q), log(SO3, p, q); atol=1e-6)
        end
    end

    @testset "Parallel transport" begin
        @test isapprox(SO3, q, X, vector_transport_to(SO3minus, SO3e, X, q))
        @test isapprox(SO3, q, q * X / q, vector_transport_to(SO3plus, SO3e, X, q))
        @test isapprox(
            SO3,
            q,
            vector_transport_to(SO3, e, X, q),
            vector_transport_to(SO3zero, SO3e, X, q),
        )
    end
end


text:
code_0 = 0
code_1 = 1
code_a = "a"
code_b = "b"
l_b = Set([code_a])
l_c = Dict(code_b => code_0)
if abspath(PROGRAM_FILE) == @__FILE__
    @assert("a" ∈ l_b)
    println("OK")
end


text:
function var(a)
    avg = a/length(a)
    x = sum(a-avg)
    s_sq = x/(length(a)-1)
    return(s_sq)
end

max_stars_repo_name: arubhardwaj/Julia

max_issues_repo_name: arubhardwaj/Julia

max_forks_repo_name: arubhardwaj/Julia


text:
function ipopt_solve_run_and_store(snlp::AbstractNLPModel, problem_name::String, test_name::String, my_par::OnePhase.Class_parameters)
    solver = IpoptSolver(print_level=3, tol=my_par.term.tol_opt, max_iter=3000, max_cpu_time=60.0^2, nlp_scaling_method="none", bound_relax_factor = 0.0, acceptable_iter=999999)

    println("RUNNING $problem_name")

    ORG_STDOUT = STDOUT
    file = open("../results/$(test_name)/log/$(problem_name).txt", "w")
    redirect_stdout(file)
    summary = problem_summary2()
    start_time = time()

    try
      t = 0
      function intermediate_ipopt(alg_mod::Int,iter_count::Int,obj_value::Float64,inf_pr::Float64, inf_du::Float64,  mu::Float64, d_norm::Float64, regularization_size::Float64, alpha_du::Float64, alpha_pr::Float64, ls_trials::Int)
        t = iter_count
        return true  # Keep going
      end

      mp = NLPModels.NLPtoMPB(snlp, solver)
      setIntermediateCallback(mp.inner, intermediate_ipopt)
      MathProgBase.optimize!(mp)
      #x = MathProgBase.getsolution(mp)
      #solver = MathProgBase.getrawsolver(mp)
      status = MathProgBase.status(mp)
      summary.status = status;
      summary.it_count = t;

      x = MathProgBase.getsolution(mp)
      set_cutest_info_ipopt!(summary, mp.inner, snlp, x)
    catch(e)
      println("Uncaught error in algorithm!!!")
      @show e;
      summary.status = :ERR
      summary.it_count = -1;
    end

    summary.total_time = time() - start_time;

    redirect_stdout(ORG_STDOUT)
    finalize(snlp)
    close(file)

    println("it count = ", summary.it_count)
    println("status = ", summary.status)

    save("../results/$(test_name)/summary.jld","summary",summary)

    return summary
end

function one_phase_run_and_store(snlp::AbstractNLPModel, problem_name::String, test_name::String, my_par::OnePhase.Class_parameters)
    println("RUNNING $problem_name")
    ORG_STDOUT = STDOUT
    file = open("../results/$(test_name)/log/$(problem_name).txt", "w")
    redirect_stdout(file)
    summary = problem_summary2()
    start_time = time()

    try
        #nlp = Class_CUTEst(snlp)

        iter, status, history, t, err, timer = OnePhase.one_phase_solve(snlp, my_par)

        #master_timer = merge_timers(timer, master_timer)

        save("../results/$(test_name)/jld/$(problem_name).jld","history",history, "timer", timer)

        set_info_me!(summary, history, status)
        #.it_count = t;
    catch(e)
        println("Uncaught error in algorithm!!!")
        @show e;

        if isa(e, Eval_NaN_error)
          summary.status = :NaN_ERR
          summary.it_count = -1;
        else
          summary.status = :ERR
          summary.it_count = -1;
        end
    end
    summary.total_time = time() - start_time;

    redirect_stdout(ORG_STDOUT)
    close(file)

    println("it count = ", summary.it_count)
    println("status = ", summary.status)

    timer_file = open("../results/$(test_name)/timer.txt", "w")
    #print_timer_stats(timer_file, master_timer)
    close(timer_file)

    return summary
end


text:
module LibASICamera

import Libdl
using CEnum

if Sys.isunix()
   const libASICamera2 = joinpath(@__DIR__, "../lib/libASICamera2")

elseif Sys.isapple()  # don't know if apple counts as unix
   const libASICamera2 = joinpath(@__DIR__, "../lib/libASICamera2")

elseif Sys.iswindows()
   const libASICamera2 = joinpath(@__DIR__, "../lib/ASICamera2")
end

include(joinpath(@__DIR__, "LibASICamera_highlevel.jl"))

export
   ASICamera,
   open_camera,
   close_camera,
   init_camera,
   start_exposure,
   stop_exposure,
   start_video,
   stop_video,
   allocate_buffer,
   enable_dark_subtract,
   disable_dark_subtract,
   pulse_guide_on,
   pulse_guide_off,
   send_soft_trigger,
   capture_still
   #capture_video

# export the rest
foreach(names(@__MODULE__, all=true)) do s
   if startswith(string(s), "ASI") || startswith(string(s), "set") ||startswith(string(s), "get")
       @eval export $s
   end
end

end # module


text:
mutable struct vector_3d_t <: LCMType
    x::Float64
    y::Float64
    z::Float64
end

@lcmtypesetup(vector_3d_t)


text:
####################################################

"""
$(TYPEDSIGNATURES)

Calculate the concentration of free SDS (freeSDS) and binding number (NBound).
"""
function calcfreeSDSNbound(protcons::Vector{<:Real},
                           statpts::Array{Vector{Float64}},
                           inflpts::Array{Vector{Float64}},
                           selectstatpts::Vector{<:Array{<:Real,2}},
                           selectinflpts::Vector{<:Array{<:Real,2}},
                           outdir::String,protein::String) #Vector{Array{<:Real,2}})

    lstlines = [selectstatpts..., selectinflpts...]
    nlines = length(lstlines)

    intercept = Vector{Float64}(undef,nlines)
    angcoe = Vector{Float64}(undef,nlines)

    # variance of angcoe and intercept
    # varangcoe = Vector{Float64}(undef,nlines)
    # varintercept = Vector{Float64}(undef,nlines)

    # residuals standard deviation
    resstdev = Vector{Float64}(undef,nlines)

    for i=1:nlines
        # angcoe[i],intercept[i] = lssqregr(lstlines[i])

        # # setup covariance matrices
        # Cd = Diagonal(stdinputdata.^2 .*ones(size(lstlines[i],1)))

        # perform least squares regression
        angcoe[i],intercept[i],resstdev[i] = lssqregr(lstlines[i]) #,Cd)
    end

    ##
    ## totSDS = Nbound * protc + freeSDS
    ##
    Nbound = angcoe
    freeSDS = intercept

    # # standard dev. from variance
    # stdNbound = sqrt.(varangcoe)
    # stdfreeSDS = sqrt.(varintercept)

    # @show Nbound
    # @show freeSDS
    # @show resstdev

    ## Old stuff...
    # y = angcoe.*x .+ intercept
    # x = (-intercept .+ 0)./angcoe
    # Nbound = 1.0./angcoe
    # freeSDS = -intercept./angcoe

    # ## Propagation of uncertainty on Nbound using
    # ##  the formula  δy = | dy/dx | δx
    # stdNbound = (1.0 ./ angcoe.^2) .* stdangcoe

    # ## Propagation of uncertainty on freeSDS using
    # ##  the formula  δz = sqrt( (|dz/dx| δx)^2 + (|dz/dy| δy)^2 )
    # ## [ freeSDS = -intercept./angcoe  ]
    # stdfreeSDS = sqrt.( ((intercept/angcoe.^2).*stdangcoe)^2 .+
    #                     ((-1.0./intercept).*stdintercept).^2 )

    ## Sort in ascending order of SDS concentration
    idxsort = sortperm(freeSDS)
    freeSDS = freeSDS[idxsort]
    Nbound  = Nbound[idxsort]
    resstdev = resstdev[idxsort]
    # stdNbound = stdNbound[idxsort]
    # stdfreeSDS = stdfreeSDS[idxsort]

    ##-------------------------------------------------------
    # save in JLD2 format
    outfile1 = joinpath(outdir,protein*"_ITCbindingisoth.jld2")
    println("Saving binding isotherm data in JL2 to $outfile1\n")
    jldopen(outfile1, "w") do fl
        fl["protcons"] = protcons
        fl["statpts"] = statpts
        fl["inflpts"] = inflpts
        fl["selectedstatpts"] = selectstatpts
        fl["selectedinflpts"] = selectinflpts
        fl["freeSDS"] = freeSDS
        fl["Nbound"] = Nbound
        fl["resstdev"] = resstdev
        # fl["stdfreeSDS"] = stdfreeSDS
        # fl["stdevNbound"] = stdNbound
        # fl["stdinputdata"] = stdinputdata
    end

    ##-------------------------------------------------------
    # save in plain text format
    outdata = [freeSDS Nbound resstdev]
    outfile2 = joinpath(outdir,protein*"_ITCbindingisoth.dat")
    println("Saving results in text file to $outfile2\n")
    header = "# Binding isotherm results for protein $protein, columns: freeSDS, Nbound, stand_dev_residuals"
    ofl = open(outfile2,"w")
    writedlm(ofl,[header])
    writedlm(ofl,outdata)
    close(ofl)

    return freeSDS,Nbound,resstdev #,stdfreeSDS,stdNbound
end

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

"""
$(TYPEDSIGNATURES)

Define a set of features on the Beta mix to subsequently compute the binding isotherm.
It uses the stationary and inflection points at given protein concentrations (protcon).
"""
function findcurvefeatures(betamix::BetaMix2D,protcon::Vector{<:Real})

    function fwd1ptxy(xpt::Real,ypt::Real)
        xy = hcat(xpt,ypt)
        out = forwmod2D(betamix.betpar,xy,betamix.modkonamp)
        return out[1]
    end

    xpt,ypt = 0.0,0.0
    deriv1f(xpt::Real) = ForwardDiff.derivative(xpt->fwd1ptxy(xpt,ypt),xpt)
    deriv2f(xpt::Real) = ForwardDiff.derivative(xpt->deriv1f(xpt),xpt)

    ##==================================================
    ## get stationary and inflection points using
    ##   first and second derivatives
    lowlim = betamix.betpar.a
    uplim  = betamix.betpar.b
    ny = length(protcon)
    statpts = Array{Vector{Float64}}(undef,ny)
    inflpts = Array{Vector{Float64}}(undef,ny)
    for i=1:ny
        ypt = protcon[i] ## define ypt for closure derivf(xpt)
        statpts[i] = find_zeros(deriv1f,lowlim,uplim)
        inflpts[i] = find_zeros(deriv2f,lowlim,uplim)
    end

    return statpts,inflpts
end

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

# """
# $(TYPEDSIGNATURES)

# Define a set of features on a single individual Beta function to subsequently compute the binding isotherm.
# It uses a ratio of amplitude compared to maximum amplitude (ampratios).
# """
# function freeboundsds_singlebeta(betpar::ScaledBeta2DParams,mcur::Vector{<:Real} ;
#                                  ampratios::Union{Vector{<:Real},Nothing}=nothing,
#                                  nptsprot::Integer=10)

#     # init
#     sdsfreeNbound = Array{Float64,2}(undef,0,2) #npoints,2)

#     ## selectd values of amplitudes to trace points
#     if ampratios==nothing
#         ampratios = collect(LinRange(0.1,1.0,6))
#     end

#     ## sample points for protein concentration (vector)
#     ypos = collect(LinRange(betpar.ymin,betpar.ymax,nptsprot))

#     namppoints = length(ampratios)
#     for p=1:namppoints
#         ## find the straight line(s) for a given amplitude
#         fsds,nbou = linefreesdsNbound(betpar,mcur,ypos,ampratios[p])
#         # store results
#         sdsfreeNbound = vcat(sdsfreeNbound,[fsds nbou])
#     end

#     ## sort data
#     sdsfreeNbound = sortslices(sdsfreeNbound,dims=1)

#     return sdsfreeNbound #,sdsfNb_mode
# end

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

# function linefreesdsNbound(betpar::ScaledBeta2DParams,mcur::Vector{<:Real},
#                            ypos::Vector{<:Real},ampratio::Real)

#     ## find the point where ampl/max(ampl)==selectedampl
#     nypts = length(ypos)
#     ## swapping for least squares, y first and x last, [protein] and [SDS]...
#     sdsc = Dict()
#     protc = Dict()
#     ## left, right with respect to the mode!!
#     lssides = ["left","right"]

#     ##------------------------------------------------
#     ## find the x,y values for given amplitudes
#     ## loop on all y values
#     for i=1:nypts
#         ## the following returns 1 or 2-element array
#         xpos,ytmp,sidevec = findxbeta(betpar,mcur,ampratio,ypos[i])
#         ## y first and x last, [protein] and [SDS]...
#         for s=1:length(xpos)
#             try
#                 push!(sdsc[sidevec[s]],xpos[s])
#                 push!(protc[sidevec[s]],ytmp[s])
#             catch
#                 sdsc[sidevec[s]]  = Float64[]
#                 protc[sidevec[s]] = Float64[]
#                 push!(sdsc[sidevec[s]],xpos[s])
#                 push!(protc[sidevec[s]],ytmp[s])
#             end
#         end
#     end

#     nkeys = length(keys(sdsc))

#     ## find the straight line by least squares
#     intercept = Vector{Float64}(undef,nkeys)
#     angcoe = Vector{Float64}(undef,nkeys)
#     for (i,sd) in enumerate(keys(sdsc))
#         ## least squares on swapped x and y, i.e., [SDS] and [protein]
#         G = [protc[sd] ones(length(protc[sd]))]
#         dobs = sdsc[sd]
#         ## overdetermined least squares
#         ## mpost = (G'*G)*G'*dobs
#         ## (G'*G) * h = G'
#         lstmp = (G'*G) \ G'
#         mpost = lstmp * dobs

#         angcoe[i] = mpost[1]
#         intercept[i] = mpost[2]
#     end

#     return intercept,angcoe
# end

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

# function findxbeta(betpar::ScaledBeta2DParams,mcur::Vector{<:Real},
#                    ampratio::Real,yval::Real)

#     ## yval -> protein concentration

#     ##------------------------------------------------------
#     function misfroot(x::Real)
#         bval = scaledbeta(mode,kon,betpar.a,betpar.b,amp,x)
#         ## abs(f(x)) because we have also negative amplitudes,
#         ##  however, the beta is always either all positive or
#         ##  all negative
#         msr = bval-ampval
#         return msr
#     end
#     ##------------------------------------------------------

#     @assert 0.0<ampratio<=1.0

#     ## mode,kon,amp
#     ## get the values of model parameters at y=ycur
#     mode,kon,amp = getmodparbeta(betpar,mcur,yval)

#     ## left, right with respect to the mode!!
#     maxamp = scaledbeta(mode,kon,betpar.a,betpar.b,amp,mode)
#     ## actual amplitude value
#     ampval = ampratio*maxamp

#     lssides = ["left","right"]
#     if ampval==maxamp
#         xpos,ypos,sides = [mode],[yval],["left"]
#     else
#         xpos = Vector{Float64}(undef,2)
#         ypos = Vector{Float64}(undef,2)
#         sides = Vector{String}(undef,2)
#         for (i,sd) in enumerate(lssides)
#             sd=="left" ? (abracket=betpar.a) : (abracket=mode)
#             sd=="left" ? (bbracket=mode) : (bbracket=betpar.b)
#             ## https://github.com/JuliaMath/Roots.jl
#             xpoint = fzero(misfroot,abracket,bbracket)
#  #@show i,sd,abracket,bbracket,xpoint,sd
#             xpos[i],ypos[i],sides[i] = xpoint,yval,sd
#         end
#     end
# #@show xpos,ypos,sides
#     return xpos,ypos,sides
# end

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


text:
using MacroTools: @q

macro nograd(ex)
  isexpr(ex, :tuple) || (ex = Expr(:tuple, ex))
  blk = @q begin end
  for f in ex.args
    back = MacroTools.@q _ -> ($__source__; nothing)
    push!(blk.args, :(@inline Zygote._pullback(::Context, ::Core.Typeof($(esc(f))), args...) = $(esc(f))(args...), $back))
  end
  return blk
end

macro which(ex)
  @capture(ex, f_(args__)) || error("Zygote.@which f(args...)")
  :(InteractiveUtils.@which adjoint(Context(), $(esc(f)), $(esc.(args)...)))
end


text:
function blink!(w::TestWindow, x)
    repr(MIME("text/html"), x)
end


text:
module Ex5_16c
using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

# Example 5.16c - Solution of an H∞ AMMP
println("Example 5.16c")

# define system with control, noise and actuator fault inputs
A = [-.8 0 0; 0 -.5 .6; 0 -.6 -.5];
Bu = [1 1;1 0;0 1]; Bw = 0.25*[0 0;0 1;1 0]; Bf = Bu;
C = [0 1 1; 1 1 0]; Du = zeros(2,2); Dw = zeros(2,2);
p, mu = size(Du); mf = mu; mw = size(Dw,2);

# setup synthesis model with additive actuator faults with Bf = Bu, Df = Du
sysf = fdimodset(dss(A,[Bu Bw],C,[Du Dw]), c = 1:mu, n = mu .+ (1:mw),f = 1:mu);

# define Mr(s) = I
Mr = FDFilterIF(dss(eye(mf)),0,0,mf)

Q, R, info = ammsyn(sysf,Mr; atol = 1.e-7, nullspace = false, reltol = 5.e-4, sdeg = -10, normalize = "dcgain");

gamma_opt0 = info.gammaopt0  # optimal performance for the initial problem
gamma_opt  = info.gammaopt   # optimal performance
gamma_sub  = info.gammasub   # suboptimal performance

# check suboptimal solution
Ge = [sysf.sys[:,[sysf.faults;sysf.noise]]; zeros(mu,mf+mw)]; Me = [info.M*Mr.sys zeros(mf,mw)];
@test iszero(R.sys-Q.sys*Ge, atol = 1.e-7) &&
      fdimmperf(R,info.M*Mr) ≈ gamma_sub && glinfnorm(Me-Q.sys*Ge)[1] ≈ gamma_sub

end # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-3//2*(A+C*sin(e+f*x)^2), x) == :(1//4*a^-1*f^-1*(a+a*sin(e+f*x))^(1/2)*(c+-1*c*sin(e+f*x))^-3//2*(A+C)*cos(e+f*x)+-1//4*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*(A+-3C)*cos(e+f*x)*log(1+-1*sin(e+f*x))+1//4*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*(A+C)*cos(e+f*x)*log(1+sin(e+f*x)))
@test integrate((a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-3//2*(A+B*sin(e+f*x)+C*sin(e+f*x)^2), x) == :(1//4*a^-1*f^-1*(a+a*sin(e+f*x))^(1/2)*(c+-1*c*sin(e+f*x))^-3//2*(A+B+C)*cos(e+f*x)+-1//4*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*(A+-1B+-3C)*cos(e+f*x)*log(1+-1*sin(e+f*x))+1//4*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*(A+C+-1B)*cos(e+f*x)*log(1+sin(e+f*x)))
@test integrate(sin(e+f*x)^-3//2*(a+b*sin(e+f*x))*(A+B*sin(e+f*x)+C*sin(e+f*x)^2), x) == :(f^-1*(-2*a*(A+-1C)+2*B*b)*Elliptic.E((1/2)*e+-1//4*pi+(1/2)*f*x,2)+1//3*f^-1*(2*C*b+6*A*b+6*B*a)*Elliptic.F((1/2)*e+-1//4*pi+(1/2)*f*x,2)+-2*A*a*f^-1*sin(e+f*x)^-1//2*cos(e+f*x)+-2//3*C*b*f^-1*sin(e+f*x)^(1/2)*cos(e+f*x))
@test integrate((a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*sin(e+f*x)^3, x) == :((1/2)*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*cos(e+f*x)^3+-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*cos(e+f*x)*log(cos(e+f*x)))
@test integrate((a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*sin(e+f*x), x) == :(-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*cos(e+f*x)*log(cos(e+f*x)))
@test integrate((a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*csc(e+f*x), x) == :(f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*cos(e+f*x)*log(tan(e+f*x)))
@test integrate((a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*csc(e+f*x)^2, x) == :(-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*cot(e+f*x)+f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*arctanh(sin(e+f*x))*cos(e+f*x))
@test integrate((a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*csc(e+f*x)^3, x) == :(f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*cos(e+f*x)*log(tan(e+f*x))+-1//2*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*cot(e+f*x)^2*cos(e+f*x))
@test integrate((a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*csc(e+f*x)^4, x) == :(-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*cot(e+f*x)+f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*arctanh(sin(e+f*x))*cos(e+f*x)+-1//3*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*csc(e+f*x)^2*cot(e+f*x))
@test integrate((a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-3//2*sin(e+f*x)^3, x) == :(c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*cos(e+f*x)*log(cos(e+f*x))+(1/2)*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*sin(e+f*x)^2*tan(e+f*x)+(1/2)*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*sin(e+f*x)*tan(e+f*x)+-3//2*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*arctanh(sin(e+f*x))*cos(e+f*x)+3//2*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*cos(e+f*x)*sin(e+f*x))
@test integrate((a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-3//2*csc(e+f*x), x) == :((1/2)*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*tan(e+f*x)+c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*cos(e+f*x)*log(tan(e+f*x))+(1/2)*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*arctanh(sin(e+f*x))*cos(e+f*x)+(1/2)*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*sin(e+f*x)*tan(e+f*x))
@test integrate((a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-3//2*csc(e+f*x)^2, x) == :(-3//2*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*cot(e+f*x)+c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*cos(e+f*x)*log(tan(e+f*x))+(1/2)*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*csc(e+f*x)*sec(e+f*x)+(1/2)*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*sin(e+f*x)*tan(e+f*x)+3//2*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*arctanh(sin(e+f*x))*cos(e+f*x))
@test integrate((a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-3//2*csc(e+f*x)^3, x) == :(-3//2*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*cot(e+f*x)+(1/2)*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*csc(e+f*x)*sec(e+f*x)+(1/2)*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*sin(e+f*x)*tan(e+f*x)+2*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*cos(e+f*x)*log(tan(e+f*x))+-1//2*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*cot(e+f*x)^2*cos(e+f*x)+3//2*c^-1*f^-1*(a+a*sin(e+f*x))^-1//2*(c+-1*c*sin(e+f*x))^-1//2*arctanh(sin(e+f*x))*cos(e+f*x))

max_stars_repo_path: test/4 Trig functions/4.1 Sine/4.1.4.2 (a+b sin)^m (c+d sin)^n (A+B sin+C sin^2).jl


text:
function commandQuit(ed::Editor, args::String)
    editorQuit(ed, force=strip(args) == "!")
end

addCommand(:quit, commandQuit,
           help="quit [!]: quit acorn. run 'quit !' to force quit")


text:
# We have to keep the original rulset pointer as it may be modified
# elsewhere like in an Interact.jl interface. `Ruleset` is mutable,
# and rules has an abstract field type.
function precalcrules(sd::SimData, rules::Tuple)
    @set sd.ruleset = ModelParameters.setparent(
        ruleset(sd),
        precalcrules(_proc_setup(proc(sd), ModelParameters.stripparams(rules)), sd)
    )
end

precalcrules(rules::Tuple, simdata) =
    (precalcrule(rules[1], simdata), precalcrules(tail(rules), simdata)...)
precalcrules(rules::Tuple{}, simdata) = ()

# Interface method
precalcrule(chain::Chain{R,W}, simdata) where {R,W} =
    Chain{R,W}(precalcrules(rules(chain), simdata))
# Support for legacy pluralised version
precalcrule(rule, simdata) = precalcrules(rule, simdata)
# The default is to return a rule unchanged
precalcrules(rule, simdata) = rule


text:
"This module contains all inner parts."
module Internal

export Problem, Options, solve

using OffsetArrays
using UnPack

"""
```
Problem(
    p::Function,
    q::Function,
    f::Function,
    a::Float64,
    α::Float64,
    A::Float64,
    b::Float64,
    β::Float64,
    B::Float64,
)
```

A linear boundary value problem.

Expects the problem in the form:

```
y''(x) + p(x) y'(x) + q(x) y(x) = f(x)
y'(a) = α y(a) + A
y'(b) = β y(b) + B
```

"""
struct Problem
    p::Function
    q::Function
    f::Function

    a::Float64
    α::Float64
    A::Float64
    b::Float64
    β::Float64
    B::Float64
end

"""
`Options(n::Integer)`

Options to provide to the solver. Includes the number of nodes in the grid `n`.
"""
struct Options
    n::Integer
end

"The default set of options. Sets `n` to 1000."
DefaultOptions = Options(1000)

include("Functions/solve.jl")

end


text:
using Test, OandaTrade

# Presently only tests that it builds.

@test true

```
@testset "account" begin
    include("account_test.jl")
end

@testset "instrument" begin
    include("instrument_test.jl")
end

@testset "order" begin
    include("order_test.jl")
end

@testset "position" begin
    include("position_test.jl")
end

@testset "transaction" begin
    include("transaction_test.jl")
end

@testset "pricing" begin
    include("pricing_test.jl")
end
```

max_stars_repo_name: CalebDepatie/Julianda

max_forks_repo_name: CalebDepatie/Julianda


text:
##
# New environments

struct BoxND
	# create one with e.g. BoxND(2, [0.0, 0.0], [1.0, 1.0])
	dim::Int64
    min::Vector{Float64}
    max::Vector{Float64}
end

function Base.rand(rng::AbstractRNG, s::Random.SamplerTrivial{BoxND})
	box = s[]
	v_range = box.max - box.min

	return v_range .* rand(rng, box.dim) + box.min
end

##
# Different definitions for action distances and proposals

# Box actions
function action_dist(a1::Vector{Float64}, a2::Vector{Float64})
	return euclidean(a1, a2)
end

function action_sample_centered(best_a::Vector{Float64}, sampler, a_space::BoxND)
	in_bounds = false
	a_new = Vector{Float64}

	# make sure new sample is in bounds
	while !in_bounds
		a_new = rand(sampler) + best_a
		in_bounds = minimum(a_new .> a_space.min) && minimum(a_new .< a_space.max)
	end

	return a_new
end

# TagAction hybrid actions
function action_dist(a1::TagAction, a2::TagAction)
	if a1.look == a2.look
		angle_diff = euclidean(a1.angle, a2.angle)
		return minimum([angle_diff, 2 * pi - angle_diff])
	else
		return Inf
	end
end

function action_sample_centered(best_a::TagAction, sampler, a_space)
	look_temp = best_a.look
	angle_temp = mod(rand(sampler)[1] + best_a.angle, 2 * pi)
	return TagAction(look_temp, angle_temp)
end


text:
module RandUtils

using Random
using StableRNGs
using SparseArrays

samerand(args...) = rand(StableRNG(1), args...)
samesprand(args...) = sprand(StableRNG(1), args...)

samesprandbool(args...) = sprand(StableRNG(1), Bool, args...)

samerandstring(n) = randstring(StableRNG(1), n)

randvec(T, n) = samerand(T, n)
randvec(n) = samerand(n)

randmat(T, n) = samerand(T, n, n)
randmat(n) = samerand(n, n)

export samerand, samesprand, samesprandbool, samerandstring, randvec, randmat, StableRNGs

end


text:
# New version ?
# The fact that we create an aligignment for each exon mess up modeller and does weird things (the hub of CAMK2B is awful)
# but if we add everything on the same line everything is nice.
# the could do this:
# if df.nb_df is the same as the last iteration, append the line, else, start a new one

using BioAlignments
using CSV
using MIToS.MSA
using BioStructures

transcript_name = ARGS[1]

function get_seq_id(d, chain)
    seq_id = Int64[]
    for (i,v) in enumerate(d["_atom_site.label_atom_id"])
        if v =="CA"
            if d["_atom_site.label_asym_id"][i] == chain
                push!(seq_id, parse(Int64, d["_atom_site.label_seq_id"][i]))
            end
        end
    end
    # unique is there to prevent any issues caused by the field ["_atom_site.label_alt_id"]
    # sometimes, the CA is present 2 times for the same number, but the field ["_atom_site.label_alt_id"] differs
    # causing the number to be pushed 2 times
    unique(seq_id)
end

function get_index(seq_id, value)
    for (i,v) in enumerate(seq_id)
        if v == value
            return i
        end
    end
end

open("$(transcript_name)_reconstructed.pir", "w") do f
    df_number = -1

    df = CSV.read("$(transcript_name)_reconstructed.csv", copycols=true)
# println(df)

    open("$(transcript_name)_annotated.pir","r") do e
        global seq_lines = readlines(e)
    end

# we may need 2 loops :
# 1 that appends the query sequence
# 1 that writes each exons from the pir files on each lines

    global queryseq = ""
    global length_gaps = 0
    for i in 1:size(df, 1)
        # if there is no templates found (first_pos_pir == -1; then just add the fasta sequence)
        # add the size of the exon to last_pos_pir
        if df.first_pos_pir[i] ==-1
#         println("not found")
            fasta_start = findfirst(df.symbol[i],seq_lines[2]).start
            fasta_end = findlast(df.symbol[i], seq_lines[2]).stop
#         println(seq_lines[3][fasta_start:fasta_end])
            df.last_pos_pir[i] = fasta_end-fasta_start
            global queryseq = queryseq * seq_lines[3][fasta_start:fasta_end]
        else
            global pir_msa = read("$(transcript_name)_$(df.nb_df[i]).pir", PIR, deletefullgaps=false)
        # the position from the pir are the good ones
        # !!! keep the gaps, we need them, do no do adjustreference!
            exon_qseq = pir_msa[1:1,df.first_pos_pir[i]:df.last_pos_pir[i]]
            global queryseq = queryseq * (join(exon_qseq))
        end
    end
#here, get the query_sequence name + change the header to get the good positions
    query_name = sequencenames(pir_msa)[1]
    query_length = length([ i for i in queryseq  if i != '-'])
    query_header = split(getannotsequence(pir_msa, query_name, "Title"), ":")
    query_header[3] = "1"
    query_header[5] = string(query_length)

    write(f, ">P1;" * split(query_name, "_")[1] * "_" * split(query_name, "_")[2] * "\n")
    write(f, join(query_header, ":") * "\n")
    write(f, queryseq * '*')
    write(f, "\n")

    global templateseq = ""
    global name = ""
    global length_gaps = 0
# global template_header = ""
    first_pos_array = Int64[]
    last_pos_array = Int64[]
    for i in 1:size(df,1)
        if df.first_pos_pir[i] == -1
            global length_gaps = length_gaps + df.last_pos_pir[i]
        elseif df.nb_df[i] != df_number
            if !isempty(templateseq)
#             println("should print the sequence and start a new one")
                template_name = split(name[1], "_")[1]
                template_chain = split(name[1], "_")[2]
                d = MMCIFDict("$(template_name).cif")
                seq_id = get_seq_id(d, template_chain)
                if length(replace(templateseq, "-" => "")) == 0
                    println("template $(template_name) should not be the one taken, it is only gaps")
                else
#             println(seq_id)
#             println(get_index(seq_id, minimum(first_pos_array)))
#             println(length(replace(templateseq, "-" => "")))
                    template_header[5] = string(seq_id[get_index(seq_id, minimum(first_pos_array)) + length(replace(templateseq, "-"=> ""))-1])
                    end_gaps = length(queryseq) - (length_gaps + length(templateseq))
                    template_header[3] = string(minimum(first_pos_array))
                    write(f,">P1;" * name[1] * "\n")
                    write(f,join(template_header, ":")* "\n")
                    write(f,"-" ^ length_gaps * templateseq * "-" ^ end_gaps * "*")
                    write(f, "\n")
                end
                global length_gaps = length_gaps + length(templateseq)
                global templateseq = ""
                global pir_msa = read("$(transcript_name)_$(df.nb_df[i]).pir", PIR, deletefullgaps=false)
                fragment = pir_msa[2:2,df.first_pos_pir[i]:df.last_pos_pir[i]]
                full_seq = replace(join(getsequence(pir_msa, 2)), "-" => "")
                fragment_seq = replace(join(getsequence(fragment, 1)), "-" => "")
                name = sequencenames(fragment)
                first_pos_array = []
                last_pos_array = []
                pos_start = findfirst(fragment_seq, full_seq).start
                pos_end = findfirst(fragment_seq, full_seq).stop
                global template_header = split(getannotsequence(pir_msa, name[1], "Title"), ":")
                starting_position = parse(Int64, template_header[3])
                push!(first_pos_array, starting_position + pos_start-1)
                push!(last_pos_array, starting_position + pos_end-1)
                global templateseq = templateseq * join(fragment)
            else
#             println("get the right exons etc.")
                global pir_msa = read("$(transcript_name)_$(df.nb_df[i]).pir", PIR, deletefullgaps=false)
                fragment = pir_msa[2:2,df.first_pos_pir[i]:df.last_pos_pir[i]]
                full_seq = replace(join(getsequence(pir_msa, 2)), "-" => "")
                fragment_seq = replace(join(getsequence(fragment, 1)), "-" => "")
                name = sequencenames(fragment)
                pos_start = findfirst(fragment_seq, full_seq).start
                pos_end = findfirst(fragment_seq, full_seq).stop
                template_header = split(getannotsequence(pir_msa, name[1], "Title"), ":")
                starting_position = parse(Int64, template_header[3])
                push!(first_pos_array, starting_position + pos_start-1)
                push!(last_pos_array, starting_position + pos_end-1)
                global templateseq = templateseq * join(fragment)
            end
            df_number = df.nb_df[i]
        else
#         println("get the right exons etc.")
            global pir_msa = read("$(transcript_name)_$(df.nb_df[i]).pir", PIR, deletefullgaps=false)
            fragment = pir_msa[2:2,df.first_pos_pir[i]:df.last_pos_pir[i]]
            full_seq = replace(join(getsequence(pir_msa, 2)), "-" => "")
            fragment_seq = replace(join(getsequence(fragment, 1)), "-" => "")
            name = sequencenames(fragment)
            pos_start = findfirst(fragment_seq, full_seq).start
            pos_end = findfirst(fragment_seq, full_seq).stop
            template_header = split(getannotsequence(pir_msa, name[1], "Title"), ":")
            starting_position = parse(Int64, template_header[3])
            push!(first_pos_array, starting_position + pos_start-1)
            push!(last_pos_array, starting_position + pos_end-1)
            global templateseq = templateseq * join(fragment)
        end
    end
    template_name = split(name[1], "_")[1]
    template_chain = split(name[1], "_")[2]
    d = MMCIFDict("$(template_name).cif")
    seq_id = get_seq_id(d, template_chain)
# println(seq_id)
# println(get_index(seq_id, minimum(first_pos_array)))
# println(length(replace(templateseq, "-" => "")))
    if length(replace(templateseq, "-" => "")) == 0
            println("template $(template_name) should not be the one taken, it is only gaps")
    else
        write(f,">P1;" * name[1]* "\n")
        template_header[3] = string(minimum(first_pos_array))
        template_header[5] = string(seq_id[get_index(seq_id, minimum(first_pos_array)) + length(replace(templateseq, "-"=> ""))-1])
        end_gaps = length(queryseq) - (length_gaps + length(templateseq))
        write(f,join(template_header, ":")* "\n")
        write(f,"-" ^ length_gaps * templateseq * "-" ^ end_gaps * "*")
    end
end

max_stars_repo_name: elolaine/PhyloSofS

max_issues_repo_name: elolaine/Phylosofs

max_forks_repo_name: PhyloSofS-Team/PhyloSofS


text:
@testset "report no method matching" begin
    # if there is no method matching case, it should be reported
    let
        # NOTE: we can't just wrap them into `let`, closures can't be inferred correctly
        m = gen_virtual_module()
        interp, frame = Core.eval(m, quote
            foo(a::Integer) = :Integer
            $(profile_call)((AbstractString,)) do a
                foo(a)
            end
        end)
        @test length(interp.reports) === 1
        report = first(interp.reports)
        @test report isa NoMethodErrorReport &&
            report.atype === Tuple{typeof(m.foo), AbstractString}
    end

    # we want to get report on `zero(Any)` for this case, but `Any`-typed statement can't
    # propagate to the error points ...
    let
        interp, report = profile_call(()->sum([]))
        @test_broken !isempty(interp.reports)
    end

    # if there is no method matching case in union-split, it should be reported
    let
        m = gen_virtual_module()
        interp, frame = Core.eval(m, quote
            foo(a::Integer) = :Integer
            foo(a::AbstractString) = "AbstractString"

            $(profile_call)(a->foo(a), (Union{Nothing,Int},))
        end)

        @test length(interp.reports) === 1
        report = first(interp.reports)
        @test report isa NoMethodErrorReport && report.atype === Tuple{typeof(m.foo), Union{Nothing,Int}}
    end
end

@testset "report undefined slots" begin
    let
        interp, frame = profile_call((Bool,)) do b
            if b
                bar = rand(Int)
                return bar
            end
            return bar # undefined in this pass
        end
        @test length(interp.reports) === 1
        @test first(interp.reports) isa LocalUndefVarErrorReport
        @test first(interp.reports).name === :bar
    end

    # deeper level
    let
        m = @def begin
            function foo(b)
                if b
                    bar = rand(Int)
                    return bar
                end
                return bar # undefined in this pass
            end
            baz(a) = foo(a)
        end

        interp, frame = Core.eval(m, :($(profile_call)(baz, (Bool,))))
        @test length(interp.reports) === 1
        @test first(interp.reports) isa LocalUndefVarErrorReport
        @test first(interp.reports).name === :bar

        # works when cached
        interp, frame = Core.eval(m, :($(profile_call)(baz, (Bool,))))
        @test length(interp.reports) === 1
        @test first(interp.reports) isa LocalUndefVarErrorReport
        @test first(interp.reports).name === :bar
    end

    # try to exclude false negatives as possible (by collecting reports in after-optimization pass)
    let
        interp, frame = profile_call((Bool,)) do b
            if b
                bar = rand()
            end

            return if b
                return bar # this shouldn't be reported
            else
                return nothing
            end
        end
        @test isempty(interp.reports)
    end

    let
        interp, frame = profile_call((Bool,)) do b
            if b
                bar = rand()
            end

            return if b
                return nothing
            else
                # ideally we want to have report for this pass, but tons of work will be
                # needed to report this pass
                return bar
            end
        end
        @test_broken length(interp.reports) === 1 &&
            first(interp.reports) isa LocalUndefVarErrorReport &&
            first(interp.reports).name === :bar
    end
end

@testset "report undefined (global) variables" begin
    let
        interp, frame = profile_call(()->foo)
        @test length(interp.reports) === 1
        @test first(interp.reports) isa GlobalUndefVarErrorReport
        @test first(interp.reports).name === :foo
    end

    # deeper level
    let
        m = @def begin
            foo(bar) = bar + baz
            qux(a) = foo(a)
        end

        interp, frame = Core.eval(m, :($(profile_call)(qux, (Int,))))
        @test length(interp.reports) === 1
        @test first(interp.reports) isa GlobalUndefVarErrorReport
        @test first(interp.reports).name === :baz

        # works when cached
        interp, frame = Core.eval(m, :($(profile_call)(qux, (Int,))))
        @test length(interp.reports) === 1
        @test first(interp.reports) isa GlobalUndefVarErrorReport
        @test first(interp.reports).name === :baz
    end
end

@testset "report non-boolean condition error" begin
    let
        interp, frame = profile_call((Int,)) do a
            a ? a : nothing
        end
        @test length(interp.reports) === 1
        er = first(interp.reports)
        @test er isa NonBooleanCondErrorReport
        @test er.t === Int
    end

    let
        interp, frame = profile_call((Any,)) do a
            a ? a : nothing
        end
        @test isempty(interp.reports)
    end

    let
        interp, frame = profile_call() do
            anyary = Any[1,2,3]
            first(anyary) ? first(anyary) : nothing
        end
        @test isempty(interp.reports) # very untyped, we can't report on this ...
    end
end

@testset "inference with virtual global variable" begin
    let
        vmod = gen_virtual_module()
        res, interp = @profile_toplevel vmod begin
            s = "julia"
            sum(s)
        end

        @test is_abstract(vmod, :s)
        @test isa_abstract(vmod.s, String)
        test_sum_over_string(res)
    end

    @testset "union assignment" begin
        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod  begin
                global globalvar
                if rand(Bool)
                    globalvar = "String"
                else
                    globalvar = :Symbol
                end
            end

            @test is_abstract(vmod, :globalvar)
            @test vmod.globalvar.t === Union{String,Symbol}
        end

        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                if rand(Bool)
                    globalvar = "String"
                else
                    globalvar = :Symbol
                end

                foo(s::AbstractString) = length(s)
                foo(globalvar) # union-split no method matching error should be reported
            end

            @test is_abstract(vmod, :globalvar)
            @test vmod.globalvar.t === Union{String,Symbol}
            @test length(res.inference_error_reports) === 1
            er = first(res.inference_error_reports)
            @test er isa NoMethodErrorReport &&
                er.unionsplit # should be true
        end

        # sequential
        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                if rand(Bool)
                    globalvar = "String"
                else
                    globalvar = :Symbol
                end

                foo(s::AbstractString) = length(s)
                foo(globalvar) # union-split no method matching error should be reported

                globalvar = 10
                foo(globalvar) # no method matching error should be reported
            end

            @test is_abstract(vmod, :globalvar)
            @test isa_abstract(vmod.globalvar, Int)
            @test length(res.inference_error_reports) === 2
            let er = first(res.inference_error_reports)
                @test er isa NoMethodErrorReport &&
                er.unionsplit
            end
            let er = last(res.inference_error_reports)
                @test er isa NoMethodErrorReport &&
                !er.unionsplit
            end
        end
    end

    @testset "invalidate code cache" begin
        let
            res, interp = @profile_toplevel begin
                foo(::Integer) = "good call, pal"
                bar() = a

                a = 1
                foo(bar()) # no method error should NOT be reported

                a = '1'
                foo(bar()) # no method error should be reported

                a = 1
                foo(bar()) # no method error should NOT be reported
            end
            @test length(res.inference_error_reports) === 1
            er = first(res.inference_error_reports)
            @test er isa NoMethodErrorReport &&
                first(er.st).file === Symbol(@__FILE__) &&
                first(er.st).line === (@__LINE__) - 9 &&
                er.atype <: Tuple{Any,Char}
        end
    end
end

@testset "UndefKeywordError" begin
    let
        m = gen_virtual_module()
        interp, frame = Core.eval(m, quote
            foo(a; #= can be undef =# kw) =  a, kw
            $(profile_call)(foo, (Any,))
        end)
        @test !isempty(interp.reports)
        @test any(interp.reports) do r
            r isa UndefKeywordErrorReport
            r.err.var === :kw
        end
        # there shouldn't be duplicated report for the `throw` call
        @test !any(Fix2(isa, UncaughtExceptionReport), interp.reports)
    end
end

@testset "report `throw` calls" begin
    # simplest case
    let
        interp, frame = profile_call(()->throw("foo"))
        @test !isempty(interp.reports)
        @test first(interp.reports) isa UncaughtExceptionReport
    end

    # throws in deep level
    let
        foo(a) = throw(a)
        interp, frame = profile_call(()->foo("foo"))
        @test !isempty(interp.reports)
        @test first(interp.reports) isa UncaughtExceptionReport
    end

    # don't report possibly false negative `throw`s
    let
        foo(a) = a ≤ 0 ? throw("a is $(a)") : a
        interp, frame = profile_call(foo, (Int,))
        @test isempty(interp.reports)
    end

    # constant prop sometimes helps exclude false negatives
    let
        foo(a) = a ≤ 0 ? throw("a is $(a)") : a
        interp, frame = profile_call(()->foo(0))
        @test !isempty(interp.reports)
        @test first(interp.reports) isa UncaughtExceptionReport
    end

    # report even if there're other "critical" error exist
    let
        m = gen_virtual_module()
        interp, frame = Core.eval(m, quote
            foo(a) = sum(a) # should be reported
            bar(a) = throw(a) # shouldn't be reported first
            $(profile_call)((Bool, String)) do b, s
                b && foo(s)
                bar(s)
            end
        end)
        @test length(interp.reports) === 3
        test_sum_over_string(interp.reports)
    end

    # end to end
    let
        # this should report `throw(ArgumentError("Sampler for this object is not defined")`
        interp, frame = profile_call(rand, (Char,))
        @test !isempty(interp.reports)
        @test first(interp.reports) isa UncaughtExceptionReport

        # this should not report `throw(DomainError(x, "sin(x) is only defined for finite x."))`
        interp, frame = profile_call(sin, (Int,))
        @test isempty(interp.reports)

        # again, constant prop sometimes can exclude false negatives
        interp, frame = profile_call(()->sin(Inf))
        @test !isempty(interp.reports)
        @test first(interp.reports) isa UncaughtExceptionReport
    end
end

@testset "constant analysis" begin
    # constant prop should limit false positive union-split no method reports
    let
        m = @def begin
            mutable struct P
                i::Int
                s::String
            end
            foo(p, i) = p.i = i
        end

        # "for one of the union split cases, no matching method found for signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Union{Type{Int64}, Type{String}}, v::Int64)" should be threw away
        interp, frame = Core.eval(m, :($(profile_call)(foo, (P, Int))))
        @test isempty(interp.reports)

        # works for cache
        interp, frame = Core.eval(m, :($(profile_call)(foo, (P, Int))))
        @test isempty(interp.reports)
    end

    # more cache test, constant prop should re-run in deeper level
    let
        m = @def begin
            mutable struct P
                i::Int
                s::String
            end
            foo(p, i) = p.i = i
            bar(args...) = foo(args...)
        end

        # "for one of the union split cases, no matching method found for signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Union{Type{Int64}, Type{String}}, v::Int64)" should be threw away
        interp, frame = Core.eval(m, :($(profile_call)(bar, (P, Int))))
        @test isempty(interp.reports)

        # works for cache
        interp, frame = Core.eval(m, :($(profile_call)(bar, (P, Int))))
        @test isempty(interp.reports)
    end

    # constant prop should not exclude those are not related
    let
        m = gen_virtual_module()
        interp, frame = Core.eval(m, quote
            mutable struct P
                i::Int
                s::String
            end
            function foo(p, i, s)
                p.i = i
                p.s = s
            end

            $(profile_call)(foo, (P, Int, #= invalid =# Int))
        end)

        # "for one of the union split cases, no matching method found for signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Union{Type{Int64}, Type{String}}, v::Int64)" should be threw away, while
        # "no matching method found for call signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Type{String}, v::Int64)" should be kept
        @test length(interp.reports) === 1
        er = first(interp.reports)
        @test er isa NoMethodErrorReport &&
            er.atype === Tuple{typeof(convert), Type{String}, Int}
    end

    # constant prop should narrow down union-split no method error to single no method matching error
    let
        m = gen_virtual_module()
        interp, frame = Core.eval(m, quote
            mutable struct P
                i::Int
                s::String
            end
            function foo(p, i, s)
                p.i = i
                p.s = s
            end

            $(profile_call)(foo, (P, String, Int))
        end)

        # "for one of the union split cases, no matching method found for signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Union{Type{Int64}, Type{String}}, v::String)" should be narrowed down to "no matching method found for call signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Type{Int}, v::String)"
        @test !isempty(interp.reports)
        @test any(interp.reports) do report
            return report isa NoMethodErrorReport &&
                report.atype === Tuple{typeof(convert), Type{Int}, String}
        end
        # "no matching method found for call signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Type{String}, v::Int)"
        # won't be reported since `typeinf` early escapes on `Bottom`-annotated statement
    end

    # report-throw away with constant analysis shouldn't throw away reports from the same
    # frame but with the other constants
    let
        m = gen_virtual_module()
        interp, frame = Core.eval(m, quote
            foo(a) = a<0 ? a+string(a) : a
            bar() = foo(-1), foo(1) # constant analysis on `foo(1)` shouldn't throw away reports from `foo(-1)`
            $(profile_call)(bar)
        end)
        @test !isempty(interp.reports)
        @test any(r->isa(r,NoMethodErrorReport), interp.reports)
    end

    let
        m = gen_virtual_module()
        interp, frame = Core.eval(m, quote
            foo(a) = a<0 ? a+string(a) : a
            function bar(b)
                a = b ? foo(-1) : foo(1)
                b = foo(-1)
                return a, b
            end
            $(profile_call)(bar, (Bool,))
        end)
        @test !isempty(interp.reports)
        @test_broken count(isa(report, NoMethodErrorReport) for report in interp.reports) == 2 # FIXME
    end

    @testset "constant analysis throws away false positive reports" begin
        let
            m = @def begin
                foo(a) = a > 0 ? a : "minus"
                bar(a) = foo(a) + 1
            end

            # constant propagation can reveal the error pass can't happen
            interp, frame = Core.eval(m, :($(profile_call)(()->bar(10))))
            @test isempty(interp.reports)

            # for this case, no constant prop' doesn't happen, we can't throw away error pass
            interp, frame = Core.eval(m, :($(profile_call)(bar, (Int,))))
            @test length(interp.reports) === 1
            er = first(interp.reports)
            @test er isa NoMethodErrorReport &&
                er.unionsplit &&
                er.atype ⊑ Tuple{Any,Union{Int,String},Int}

            # if we run constant prop' that leads to the error pass, we should get the reports
            interp, frame = Core.eval(m, :($(profile_call)(()->bar(0))))
            @test length(interp.reports) === 1
            er = first(interp.reports)
            @test er isa NoMethodErrorReport &&
                !er.unionsplit &&
                er.atype ⊑ Tuple{Any,String,Int}
        end

        # we should throw-away reports collected from frames that are revealed as "unreachable"
        # by constant prop'
        let
            m = @def begin
                foo(a) = bar(a)
                function bar(a)
                    return if a < 1
                        baz1(a, "0")
                    else
                        baz2(a, a)
                    end
                end
                baz1(a, b) = a ? b : b
                baz2(a, b) = a + b
            end

            # no constant prop, just report everything
            interp, frame = Core.eval(m, :($(profile_call)(foo, (Int,))))
            @test length(interp.reports) === 1
            er = first(interp.reports)
            @test er isa NonBooleanCondErrorReport &&
                er.t === Int

            # constant prop should throw away the non-boolean condition report from `baz1`
            interp, frame = Core.eval(m, quote
                $(profile_call)() do
                    foo(1)
                end
            end)
            @test isempty(interp.reports)

            # constant prop'ed, still we want to have the non-boolean condition report from `baz1`
            interp, frame = Core.eval(m, quote
                $(profile_call)() do
                    foo(0)
                end
            end)
            @test length(interp.reports) === 1
            er = first(interp.reports)
            @test er isa NonBooleanCondErrorReport &&
                er.t === Int

            # so `Bool` is good for `foo` after all
            interp, frame = Core.eval(m, :($(profile_call)(foo, (Bool,))))
            @test isempty(interp.reports)
        end

        # end to end
        let
            res, interp = @profile_toplevel begin
                function foo(n)
                    if n < 10
                        return n
                    else
                        return "over 10"
                    end
                end

                function bar(n)
                    if n < 10
                        return foo(n) + 1
                    else
                        return foo(n) * "+1"
                    end
                end

                bar(1)
                bar(10)
            end

            @test isempty(res.inference_error_reports)
        end
    end
end

@testset "keyword argument methods" begin
    interp, frame = Core.eval(gen_virtual_module(), quote
        f(a; b = nothing, c = nothing) = return
        $(profile_call)((Any,)) do b
            f(1; b)
        end
    end)
    @test isempty(interp.reports)
end

@testset "don't early escape if type grows up to `Any`" begin
    vmod = gen_virtual_module()
    res, interp = @profile_toplevel vmod begin
        abstract type Foo end
        struct Foo1 <: Foo
            bar
        end
        struct Foo2 <: Foo
            baz # typo
        end
        struct Foo3 <: Foo
            qux # typo
        end
        bar(foo::Foo) = foo.bar

        f = rand(Bool) ? Foo1(1) : rand(Bool) ? Foo2(1) : Foo3(1)
        bar(f)
    end

    # `bar(::Foo1)` is valid but its return type is `Any`, and so we can't collect possible
    # error points for `bar(::Foo2)` and `bar(::Foo3)` if we bail out on `Any`-grew return type
    @test length(res.inference_error_reports) === 2
    @test any(res.inference_error_reports) do er
        return er isa NoFieldErrorReport &&
            er.typ === vmod.Foo2 &&
            er.name === :bar
    end
    @test any(res.inference_error_reports) do er
        return er isa NoFieldErrorReport &&
            er.typ === vmod.Foo3 &&
            er.name === :bar
    end
end


text:
import CartesianGrids: convective_derivative!, divergence!, grad!, curl!

for f in [:divergence!, :grad!, :curl!]
  @eval $f(a,b,sys::ILMSystem) = $f(a,b,sys.base_cache)
end

@inline _scale_derivative!(w,cache::BasicILMCache{N,IndexScaling}) where {N} = w
@inline _scale_derivative!(w,cache::BasicILMCache{N,GridScaling}) where {N} = w ./= cellsize(cache)

_scale_laplacian!(w,cache::BasicILMCache{N,IndexScaling}) where {N} = w
_scale_laplacian!(w,cache::BasicILMCache{N,GridScaling}) where {N} = w # ./= cellsize(cache)^2

_scale_inverse_laplacian!(w,cache::BasicILMCache{N,IndexScaling}) where {N} = w
_scale_inverse_laplacian!(w,cache::BasicILMCache{N,GridScaling}) where {N} = w  #.*= cellsize(cache)^2

"""
    divergence!(p::Nodes{Primal/Dual},v::Edges{Primal/Dual},cache::BasicILMCache)
    divergence!(p::Nodes{Primal/Dual},v::Edges{Primal/Dual},sys::ILMSystem)

Compute the discrete divergence of `v` and return it in `p`, scaling it
by the grid spacing if `cache` (or `sys`) is of `GridScaling` type, or leaving it
as a simple differencing if `cache` (or `sys`) is of `IndexScaling` type.
"""
function divergence!(p::Nodes{C,NX,NY},q::Edges{C,NX,NY},cache::BasicILMCache) where {C,NX,NY}
    fill!(p,0.0)
    if !iszero(q)
      _unscaled_divergence!(p,q,cache)
      _scale_derivative!(p,cache)
    end
end

_unscaled_divergence!(p,q,cache::BasicILMCache) = divergence!(p,q)

"""
    grad!(v::Edges{Primal/Dual},p::Nodes{Primal/Dual},cache::BasicILMCache)
    grad!(v::Edges{Primal/Dual},p::Nodes{Primal/Dual},sys::ILMSystem)

Compute the discrete gradient of `p` and return it in `v`, scaling it
by the grid spacing if `cache` (or `sys`) is of `GridScaling` type, or leaving it
as a simple differencing if `cache` (or `sys`) is of `IndexScaling` type.
"""
function grad!(q::Edges{C,NX,NY},p::Nodes{C,NX,NY},cache::BasicILMCache) where {C,NX,NY}
    fill!(q,0.0)
    if !iszero(p)
      _unscaled_grad!(q,p,cache)
      _scale_derivative!(q,cache)
    end
end

_unscaled_grad!(q,p,cache::BasicILMCache) = grad!(q,p)

"""
    curl!(v::Edges{Primal},s::Nodes{Dual},cache::BasicILMCache)
    curl!(v::Edges{Primal},s::Nodes{Dual},sys::ILMSystem)

Compute the discrete curl of `s` and return it in `v`, scaling it
by the grid spacing if `cache` (or `sys`) is of `GridScaling` type, or leaving it
as a simple differencing if `cache` (or `sys`) is of `IndexScaling` type.
"""
function curl!(q::Edges{Primal,NX,NY},s::Nodes{Dual,NX,NY},cache::BasicILMCache) where {NX,NY}
    fill!(q,0.0)
    if !iszero(s)
      _unscaled_curl!(q,s,cache)
      _scale_derivative!(q,cache)
    end
end

_unscaled_curl!(q::Edges,s::Nodes,cache::BasicILMCache) = curl!(q,s)

"""
    curl!(w::Nodes{Dual},v::Edges{Primal},cache::BasicILMCache)
    curl!(w::Nodes{Dual},v::Edges{Primal},sys::ILMSystem)

Compute the discrete curl of `v` and return it in `w`, scaling it
by the grid spacing if `cache` (or `sys`) is of `GridScaling` type, or leaving it
as a simple differencing if `cache` (or `sys`) is of `IndexScaling` type.
"""
function curl!(w::Nodes{Dual,NX,NY},q::Edges{Primal,NX,NY},cache::BasicILMCache) where {NX,NY}
    fill!(w,0.0)
    if !iszero(q)
      _unscaled_curl!(w,q,cache)
      _scale_derivative!(w,cache)
    end
end

_unscaled_curl!(w::Nodes,q::Edges,cache::BasicILMCache) = curl!(w,q)

"""
    laplacian!(w::GridData,s::GridData,sys::ILMSystem)

Compute the Laplacian of grid data `s`, and divide the result
by unity or by the grid cell size, depending on whether `sys` has `IndexScaling` or `GridScaling`,
respectively, and return the result as `w`.
"""
laplacian!(w,s,sys::ILMSystem) = laplacian!(w,s,sys.base_cache)

"""
    laplacian!(w::GridData,s::GridData,cache::BasicILMCache)

Compute the Laplacian of grid data `s`, and divide the result
by unity or by the grid cell size, depending on whether `cache` has `IndexScaling` or `GridScaling`,
respectively, and return the result as `w`.
"""
function laplacian!(w::GridData,s::GridData,cache::BasicILMCache)
    @unpack L = cache
    _unscaled_laplacian!(w,s,L)
    _scale_laplacian!(w,cache)
end

function _unscaled_laplacian!(w::GridData,s::GridData,L::CartesianGrids.Laplacian)
  fill!(w,0.0)
  if !iszero(s)
    mul!(w,L,s)
  end
end

"""
    inverse_laplacian!(w::GridData,sys::ILMSystem)

Compute the in-place inverse Laplacian of grid data `w`, and multiply the result
by unity or by the grid cell size, depending on whether `sys` has `IndexScaling` or `GridScaling`,
respectively.
"""
inverse_laplacian!(w,sys::ILMSystem) = inverse_laplacian!(w,sys.base_cache)

"""
    inverse_laplacian!(w::GridData,cache::BasicILMCache)

Compute the in-place inverse Laplacian of grid data `w`, and multiply the result
by unity or by the grid cell size, depending on whether `cache` has `IndexScaling` or `GridScaling`,
respectively.
"""
function inverse_laplacian!(w::GridData,cache::BasicILMCache)
    @unpack L = cache
    if !iszero(w)
      _unscaled_inverse_laplacian!(w,L)
      _scale_inverse_laplacian!(w,cache)
    end
end

"""
    inverse_laplacian!(sol::GridData,rhs::GridData,cache::BasicILMCache)

Compute the iinverse Laplacian of grid data `rhs`, multiply the result
by unity or by the grid cell size, depending on whether `cache` has `IndexScaling` or `GridScaling`,
respectively, and return the result as `sol`.
"""
function inverse_laplacian!(sol::GridData{NX,NY},rhs::GridData{NX,NY},cache::BasicILMCache) where {NX,NY}
    @unpack L = cache
    fill!(sol,0.0)
    if !iszero(rhs)
      _unscaled_inverse_laplacian!(sol,L,rhs)
      _scale_inverse_laplacian!(sol,cache)
    end
end

_unscaled_inverse_laplacian!(w::GridData,L::CartesianGrids.Laplacian) = w .= L\w

_unscaled_inverse_laplacian!(sol::GridData,L::CartesianGrids.Laplacian,rhs::GridData) = ldiv!(sol,L,rhs)

###  Convective derivative of velocity-like data by itself ###

"""
    ConvectiveDerivativeCache(dv::EdgeGradient)

Create a cache (a subtype of [`AbstractExtraILMCache`](@ref)) for computing
the convective derivative, using `dv` to define the cache data.
"""
struct ConvectiveDerivativeCache{VTT} <: AbstractExtraILMCache
   vt1_cache :: VTT
   vt2_cache :: VTT
   vt3_cache :: VTT
   ConvectiveDerivativeCache(vt::GridData) = new{typeof(vt)}(similar(vt),similar(vt),similar(vt))
end

"""
    ConvectiveDerivativeCache(cache::BasicILMCache)

Create a cache for computing the convective derivative, based on
the basic ILM cache `cache`.
"""
ConvectiveDerivativeCache(cache::BasicILMCache) = ConvectiveDerivativeCache(similar_gridgrad(cache))

"""
    convective_derivative(v::Edges{Primal},p::Nodes{Primal},base_cache::BasicILMCache)

Compute the convective derivative of `p` with velocity `v`, i.e., ``v\\cdot \\nabla p``. The result is either divided by unity or
the grid cell size depending on whether `base_cache` is of type `IndexScaling`
or `GridScaling`.
"""
function convective_derivative(u::Edges{Primal},p::Nodes{Primal},base_cache::BasicILMCache)
    extra_cache = ConvectiveDerivativeCache(similar_gridgrad(cache))
    udp = similar(p)
    convective_derivative!(udp,u,p,base_cache,extra_cache)
end

"""
    convective_derivative(v::Edges{Dual},w::Nodes{Dual},base_cache::BasicILMCache)

Compute the convective derivative of `w` with velocity `v`, i.e., ``v\\cdot \\nabla w``. The result is either divided by unity or
the grid cell size depending on whether `base_cache` is of type `IndexScaling`
or `GridScaling`.
"""
function convective_derivative(u::Edges{Dual},w::Nodes{Dual},base_cache::BasicILMCache)
    extra_cache = ConvectiveDerivativeCache(similar_gridgradcurl(cache))
    udw = similar(w)
    convective_derivative!(udw,u,w,base_cache,extra_cache)
end

"""
    convective_derivative(v::Edges,base_cache::BasicILMCache)

Compute the convective derivative of `v`, i.e., ``v\\cdot \\nabla v``. The result is either divided by unity or
the grid cell size depending on whether `base_cache` is of type `IndexScaling`
or `GridScaling`.
"""
function convective_derivative(u::Edges{Primal},base_cache::BasicILMCache)
    extra_cache = ConvectiveDerivativeCache(similar_gridgrad(base_cache))
    udu = similar(u)
    convective_derivative!(udu,u,base_cache,extra_cache)
end

"""
    convective_derivative!(vdp::Nodes{Primal},v::Edges,p::Nodes{Primal},base_cache::BasicILMCache,extra_cache::ConvectiveDerivativeCache)

Compute the convective derivative of `p` with velocity `v`, i.e., ``v\\cdot \\nabla p``, and
return the result in `vdp`. The result is either divided by unity or
the grid cell size depending on whether `base_cache` is of type `IndexScaling`
or `GridScaling`. This version of the method uses `extra_cache` of type
[`ConvectiveDerivativeCache`](@ref).
"""
function convective_derivative!(udp::Nodes{Primal},u::Edges{Primal},p::Nodes{Primal},base_cache::BasicILMCache,extra_cache::ConvectiveDerivativeCache)
    fill!(udp,0.0)
    if !iszero(p) && !iszero(u)
      _unscaled_convective_derivative!(udp,u,p,extra_cache)
      _scale_derivative!(udp,base_cache)
    end
end

"""
    convective_derivative!(vdw::Nodes{Dual},v::Edges,w::Nodes{Dual},base_cache::BasicILMCache,extra_cache::ConvectiveDerivativeCache)

Compute the convective derivative of `w` with velocity `v`, i.e., ``v\\cdot \\nabla w``, and
return the result in `vdw`. The result is either divided by unity or
the grid cell size depending on whether `base_cache` is of type `IndexScaling`
or `GridScaling`. This version of the method uses `extra_cache` of type
[`ConvectiveDerivativeCache`](@ref).
"""
function convective_derivative!(udw::Nodes{Dual},u::Edges{Primal},w::Nodes{Dual},base_cache::BasicILMCache,extra_cache::ConvectiveDerivativeCache)
    fill!(udw,0.0)
    if !iszero(w) && !iszero(u)
      _unscaled_convective_derivative!(udw,u,w,extra_cache)
      _scale_derivative!(udw,base_cache)
    end
end

"""
    convective_derivative!(vdu::Edges,v::Edges,u::Edges,base_cache::BasicILMCache,extra_cache::ConvectiveDerivativeCache)

Compute the convective derivative of `u`, i.e., ``v\\cdot \\nabla u``, and
return the result in `vdu`. The result is either divided by unity or
the grid cell size depending on whether `base_cache` is of type `IndexScaling`
or `GridScaling`. This version of the method uses `extra_cache` of type
[`ConvectiveDerivativeCache`](@ref).
"""
function convective_derivative!(vdu::Edges{Primal},v::Edges{Primal},u::Edges{Primal},base_cache::BasicILMCache,extra_cache::ConvectiveDerivativeCache)
    fill!(vdu,0.0)
    if !iszero(u) && !iszero(v)
      _unscaled_convective_derivative!(vdu,v,u,extra_cache)
      _scale_derivative!(vdu,base_cache)
    end
end

"""
    convective_derivative!(vdv::Edges,v::Edges,base_cache::BasicILMCache,extra_cache::ConvectiveDerivativeCache)

Compute the convective derivative of `v`, i.e., ``v\\cdot \\nabla v``, and
return the result in `vdv`. The result is either divided by unity or
the grid cell size depending on whether `base_cache` is of type `IndexScaling`
or `GridScaling`. This version of the method uses `extra_cache` of type
[`ConvectiveDerivativeCache`](@ref).
"""
function convective_derivative!(udu::Edges{Primal},u::Edges{Primal},base_cache::BasicILMCache,extra_cache::ConvectiveDerivativeCache)
    fill!(udu,0.0)
    if !iszero(u)
      _unscaled_convective_derivative!(udu,u,extra_cache)
      _scale_derivative!(udu,base_cache)
    end
end

function _unscaled_convective_derivative!(udp::Nodes{Primal},u::Edges{Primal},p::Nodes{Primal},extra_cache::ConvectiveDerivativeCache)
    @unpack vt1_cache, vt3_cache = extra_cache

    fill!(vt1_cache,0.0)
    grad!(vt1_cache,p)
    product!(vt3_cache,u,vt1_cache)
    fill!(udp,0.0)
    grid_interpolate!(udp,vt3_cache)
    udp
end

function _unscaled_convective_derivative!(udw::Nodes{Dual},u::Edges{Primal},w::Nodes{Dual},extra_cache::ConvectiveDerivativeCache{VT}) where {VT<:Edges{Dual}}
    @unpack vt1_cache, vt2_cache, vt3_cache = extra_cache

    fill!(vt1_cache,0.0)
    grad!(vt1_cache,w)
    fill!(vt2_cache,0.0)
    grid_interpolate!(vt2_cache,u)
    product!(vt3_cache,vt2_cache,vt1_cache)
    fill!(udw,0.0)
    grid_interpolate!(udw,vt3_cache)
    udw
end

function _unscaled_convective_derivative!(udu::Edges{Primal},u::Edges{Primal},extra_cache::ConvectiveDerivativeCache)
    @unpack vt1_cache, vt2_cache, vt3_cache = extra_cache

    fill!(vt1_cache,0.0)
    grid_interpolate!(vt1_cache,u)
    CartesianGrids.transpose!(vt2_cache,vt1_cache)
    fill!(vt1_cache,0.0)
    grad!(vt1_cache,u)
    product!(vt3_cache,vt2_cache,vt1_cache)
    fill!(udu,0.0)
    grid_interpolate!(udu,vt3_cache)
    udu
end

function _unscaled_convective_derivative!(vdu::Edges{Primal},v::Edges{Primal},u::Edges{Primal},extra_cache::ConvectiveDerivativeCache)
    @unpack vt1_cache, vt2_cache, vt3_cache = extra_cache

    fill!(vt1_cache,0.0)
    grid_interpolate!(vt1_cache,v)
    CartesianGrids.transpose!(vt2_cache,vt1_cache)
    fill!(vt1_cache,0.0)
    grad!(vt1_cache,u)
    product!(vt3_cache,vt2_cache,vt1_cache)
    fill!(udv,0.0)
    grid_interpolate!(udv,vt3_cache)
    udu
end


text:
println("*******************************************************")
@testset "qmr" begin
	@testset "real matrix" begin
		A  = sprandn(100,100,.1)
		A  = A + SparseMatrixCSC(100.0I, 100, 100)
		n  = size(A,2)
		D  = diag(A)
		M2 = x -> Vector(D.\x)
		rhs = randn(100)
		tol = 1e-2;

		# test printing and behaviour for early stopping
		xtt = qmr(A,rhs,tol=1e-12,maxIter=3,out=2)
		@test xtt[2]==-1

		# test behaviour for zero rhs
		xtt = qmr(A,0*rhs,tol=tol,maxIter=10,out=2)
		@test xtt[2]==-9
		@test all(xtt[1].==0)
		@test length(xtt[1])==100
		@test eltype(xtt[1])==eltype(rhs)

		x1 = qmr(A,rhs ,tol=tol,maxIter=100,out=1)
		x3 = qmr(A,rhs,tol=tol,maxIter=100,x=randn(size(rhs)))
		x4 = qmr(A,rhs,tol=tol,maxIter=100,M=M2)

		@test norm(A*x1[1]-rhs)/norm(rhs) < tol
		@test norm(A*x3[1]-rhs)/norm(rhs) < tol
		@test norm(A*x4[1]-rhs)/norm(rhs) < tol
		@test norm(x3[1]-x1[1])/norm(x1[1]) < 1e-2
	end

	@testset "complex matrix" begin
		A = sprandn(100,100,.1) + 1im*(sprandn(100,100,.1));
		A  = A + SparseMatrixCSC(100.0I, 100, 100)
		D  = Vector(diag(A))
		M3 = x -> x./D
		rhs = complex(randn(100))
		tol = 1e-3;

		# test behaviour for zero rhs

		xtt = qmr(A,0.0*rhs,tol=tol,maxIter=10,out=2)
		# @test xtt[2]==-9
		# @test all(xtt[1].==0)
		# @test length(xtt[1])==100
		# @test eltype(xtt[1])==eltype(rhs)

		x1 = qmr(A,rhs,tol=tol,maxIter=100)
		x3 = qmr(A,rhs,tol=tol,maxIter=100,x=randn(size(rhs)))
		x4 = qmr(A,rhs,tol=tol,maxIter=100,M=M3)

		@test norm(A*x1[1]-rhs)/norm(rhs) < tol
		@test norm(A*x3[1]-rhs)/norm(rhs) < tol
		@test norm(x3[1]-x1[1])/norm(x1[1]) < 1e-2
	end
end


text: touch(joinpath(@__DIR__, "buildartifact"))

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork