text:
############### print_human_readable(X::AbstractMultiScaleArray
function add_separator!(first_element,toprint,level)
    if first_element
        toprint[level][end] = "+"*toprint[level][end]
    end
    return false
end

function take_first_char(s::Union{String,SubString{String}},N::Int64)     s[1:min(length(s),N)]     end

function toprint_AbstractMultiScaleArray!(toprint,X::AbstractMultiScaleArray,fields,levelmax,n_char_per_name,level = 1)
    if length(toprint) < level
        push!(toprint,Vector{String}())
    end
    first_element = true
    if level<levelmax
        for x in X.nodes
            if fieldname(typeof(x),1) == :nodes
                push!(toprint[level],  "|"*take_first_char(split(string(typeof(x)),"{")[1],n_char_per_name)  )
                first_element = add_separator!(first_element,toprint,level)
                toprint_AbstractMultiScaleArray!(toprint,x,fields,levelmax,n_char_per_name,level+1)
            else
                if fields == nothing
                    push!(toprint[level],  take_first_char(split(string(typeof(x)),"{")[1],n_char_per_name)  )
                    first_element = add_separator!(first_element,toprint,level)
                else
                    if !(typeof(x) <: AbstractMultiScaleArrayLeaf)
                        push!(toprint[level],  take_first_char(split(string(typeof(x)),"{")[1],n_char_per_name)  )
                        first_element = add_separator!(first_element,toprint,level)
                    else
                        push!(toprint[level],  join([take_first_char(string(field),n_char_per_name)*": "*string(getfield(x,Symbol(field))) for field in collect(fields)],", ")  )
                        first_element = add_separator!(first_element,toprint,level)
                    end
                end
            end
        end
    end
end

"""
```julia
print_human_readable(embryo)
# +|Tissue;                                                 |Tissue
#  +|Popula;           |Popula;           |Popula;          +|Popula;           |Popula;           |Popula
#   +Cell; Cell; Cell; +Cell; Cell; Cell; +Cell; Cell; Cell; +Cell; Cell; Cell; +Cell; Cell; Cell; +Cell; Cell; Cell

print_human_readable(embryo;NcharPerName=2)
# +|Ti;                                   |Ti
#  +|Po;         |Po;         |Po;        +|Po;         |Po;         |Po
#   +Ce; Ce; Ce; +Ce; Ce; Ce; +Ce; Ce; Ce; +Ce; Ce; Ce; +Ce; Ce; Ce; +Ce; Ce; Ce
```
Here, if the 'AbstractMultiScaleArrayLeaf's contain several fields, you can specify them with fields = [field1,field2,...]
```julia
print_human_readable(embryo;NcharPerName=2,fields=[:values])
# +|Ti;                                                                                |Ti
#  +|Po;                                                       |Po;                                                       |Po;                                                      +|Po;                                                       |Po;                                                       |Po
#   +va: [1.0, 2.0, 3.0]; va: [3.0, 2.0, 5.0]; va: [4.0, 6.0]; +va: [1.0, 2.0, 3.0]; va: [3.0, 2.0, 5.0]; va: [4.0, 6.0]; +va: [1.0, 2.0, 3.0]; va: [3.0, 2.0, 5.0]; va: [4.0, 6.0]; +va: [1.0, 2.0, 3.0]; va: [3.0, 2.0, 5.0]; va: [4.0, 6.0]; +va: [1.0, 2.0, 3.0]; va: [3.0, 2.0, 5.0]; va: [4.0, 6.0]; +va: [1.0, 2.0, 3.0]; va: [3.0, 2.0, 5.0]; va: [4.0, 6.0]
```
if your screen is small, then print a sub-part of the AbstractMultiScaleArray:
```julia
print_human_readable(embryo.nodes[1].nodes[1];fields=[:values])
# +values: [1.0, 2.0, 3.0]; values: [3.0, 2.0, 5.0]; values: [4.0, 6.0]
```
"""
function print_human_readable(X::AbstractMultiScaleArray ; n_char_per_name = 6, fields = nothing, levelmax=Inf, n_item_max_per_levels = Inf) # fields = nothing OR [field1,field2,...]
#     if typeof(X) <: AbstractMultiScaleArrayLeaf
#         println(X)
#     else
        toprint = Vector{Vector{String}}([[]]) # one vector per row, one string per element
        toprint_AbstractMultiScaleArray!(toprint,X,fields,levelmax,n_char_per_name)
        toprint = map(x -> join(x,"; "),toprint)
    #     toprint[1] = join(split(toprint[1]," "),"; ")
        Done = false
        while !Done
            Done = true # temporary
            for level1 in length(toprint):-1:2
                toprint0 = split(toprint[level1-1],"|")
                toprint1 = split(toprint[level1  ],"+")
                toprint0_ , toprint1_ = "", ""
                for (x0, x1) in zip(toprint0,toprint1)
                    if length(x0) > length(x1)
                        x1=join([x1,repeat(" ",length(x0)-length(x1))],"")
                        Done = false
                    elseif length(x0) < length(x1)
                        x0=join([x0,repeat(" ",length(x1)-length(x0))],"")
                        Done = false
                    end
                    toprint0_ = join([toprint0_, x0],"|")
                    toprint1_ = join([toprint1_, x1],"+")
                end
                toprint0_ = toprint0_[2:end]
                toprint1_ = toprint1_[2:end]
                toprint[level1-1] = toprint0_
                toprint[level1  ] = toprint1_
            end
        end
        # handle "ula; +         |Popu" cases
        #              ^^^^^^^^^^^
        for level in length(toprint):-1:1
            while occursin("+ ", toprint[level])
                Move = findfirst("+ ", toprint[level])[1]
                To   = findnext(" |", toprint[level], Move[end])[1]
                toprint[level] = join([toprint[level][1:(Move-1)],  toprint[level][(Move+1):To] , toprint[level][Move]  , toprint[level][(To+1):end]])
                for level_ in 1:(level-1)
                    toprint[level_] = join([toprint[level_][1:(Move-1)],  toprint[level_][(To+1):(To+To-Move)],  toprint[level_][Move:To],  toprint[level_][(To+To-Move+1):end]])
                end
            end
        end
        println(join(toprint,"\n"))
#     end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct StorageAccountCreateParameters <: SwaggerModel
    sku::Any # spec type: Union{ Nothing, Sku } # spec name: sku
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    location::Any # spec type: Union{ Nothing, String } # spec name: location
    tags::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: tags
    identity::Any # spec type: Union{ Nothing, Identity } # spec name: identity
    properties::Any # spec type: Union{ Nothing, StorageAccountPropertiesCreateParameters } # spec name: properties

    function StorageAccountCreateParameters(;sku=nothing, kind=nothing, location=nothing, tags=nothing, identity=nothing, properties=nothing)
        o = new()
        validate_property(StorageAccountCreateParameters, Symbol("sku"), sku)
        setfield!(o, Symbol("sku"), sku)
        validate_property(StorageAccountCreateParameters, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(StorageAccountCreateParameters, Symbol("location"), location)
        setfield!(o, Symbol("location"), location)
        validate_property(StorageAccountCreateParameters, Symbol("tags"), tags)
        setfield!(o, Symbol("tags"), tags)
        validate_property(StorageAccountCreateParameters, Symbol("identity"), identity)
        setfield!(o, Symbol("identity"), identity)
        validate_property(StorageAccountCreateParameters, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type StorageAccountCreateParameters

const _property_map_StorageAccountCreateParameters = Dict{Symbol,Symbol}(Symbol("sku")=>Symbol("sku"), Symbol("kind")=>Symbol("kind"), Symbol("location")=>Symbol("location"), Symbol("tags")=>Symbol("tags"), Symbol("identity")=>Symbol("identity"), Symbol("properties")=>Symbol("properties"))
const _property_types_StorageAccountCreateParameters = Dict{Symbol,String}(Symbol("sku")=>"Sku", Symbol("kind")=>"String", Symbol("location")=>"String", Symbol("tags")=>"Dict{String, String}", Symbol("identity")=>"Identity", Symbol("properties")=>"StorageAccountPropertiesCreateParameters")
Base.propertynames(::Type{ StorageAccountCreateParameters }) = collect(keys(_property_map_StorageAccountCreateParameters))
Swagger.property_type(::Type{ StorageAccountCreateParameters }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_StorageAccountCreateParameters[name]))}
Swagger.field_name(::Type{ StorageAccountCreateParameters }, property_name::Symbol) =  _property_map_StorageAccountCreateParameters[property_name]

const _allowed_StorageAccountCreateParameters_kind = ["Storage", "BlobStorage"]

function check_required(o::StorageAccountCreateParameters)
    (getproperty(o, Symbol("sku")) === nothing) && (return false)
    (getproperty(o, Symbol("location")) === nothing) && (return false)
    true
end

function validate_property(::Type{ StorageAccountCreateParameters }, name::Symbol, val)
    if name === Symbol("kind")
        Swagger.validate_param(name, "StorageAccountCreateParameters", :enum, val, _allowed_StorageAccountCreateParameters_kind)
    end
end

ext: jl


text:
#=
Three distinct points are plotted at random on a Cartesian plane, for which -1000 ≤ x, y ≤ 1000, such that a triangle is formed.

Consider the following two triangles:

A(-340,495), B(-153,-910), C(835,-947)

X(-175,41), Y(-421,-714), Z(574,-645)

It can be verified that triangle ABC contains the origin, whereas triangle XYZ does not.

Using triangles.txt (right click and 'Save Link/Target As...'), a 27K text file containing the co-ordinates of one thousand "random" triangles, find the number of triangles for which the interior contains the origin.

NOTE: The first two examples in the file represent the triangles in the example given above.
=#
using DelimitedFiles

function pointInTriangle(pointX, pointY, v1X, v1Y, v2X, v2Y, v3X, v3Y) # Find point in triangle using barycentric coordinates
  A = (-v2Y * v3X + v1Y * (-v2X + v3X) + v1X * (v2Y - v3Y) + v2X * v3Y) / 2;
  sign = A < 0 ? -1 : 1;
  s = (v1Y * v3X - v1X * v3Y + (v3Y - v1Y) * pointX + (v1X - v3X) * pointY) * sign;
  t = (v1X * v2Y - v1Y * v2X + (v1Y - v2Y) * pointX + (v2X - v1X) * pointY) * sign;
  return s > 0 && t > 0 && s + t < 2 * A * sign;
end

function calc()
  x = readdlm(dirname(@__FILE__()) * "/../Resources/p102.txt", ',', Int64)
  count = 0
  for i in collect(1:1000)
    count = pointInTriangle(0, 0, x[i,1], x[i, 2], x[i, 3], x[i, 4], x[i, 5], x[i, 6]) ? count + 1 : count
  end
  count
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
int_rules_7_2_5 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 7.2.5*u*(a+b*arccosh(c*x))^n =#
    @apply_utils Antiderivative((~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n') / (~(d') + ~(e') * ~x), ~x) => Subst(Antiderivative(((~a + ~b * ~x) ^ ~n * sinh(~x)) / (~c * ~d + ~e * cosh(~x)), ~x), ~x, acosh(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => ((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * acosh(~c * ~x)) ^ ~n) / (~e * (~m + 1)) - ((~b * ~c * ~n) / (~e * (~m + 1))) * Antiderivative(((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * acosh(~c * ~x)) ^ (~n - 1)) / (sqrt(-1 + ~c * ~x) * sqrt(1 + ~c * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (IGtQ(~n, 0) && NeQ(~m, -1))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~n, ~x) => Antiderivative(ExpandIntegrand((~d + ~e * ~x) ^ ~m * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IGtQ(~m, 0) && LtQ(~n, -1))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => (1 / (~c) ^ (~m + 1)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * sinh(~x) * (~c * ~d + ~e * cosh(~x)) ^ ~m, ~x), ~x, acosh(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * acosh(~(c') * ~x)), ~x) => With([u = IntHide(ExpandExpression(~Px, ~x), ~x)], Dist(~a + ~b * acosh(~c * ~x), u, ~x) - ((~b * ~c * sqrt(1 - (~c) ^ 2 * (~x) ^ 2)) / (sqrt(-1 + ~c * ~x) * sqrt(1 + ~c * ~x))) * Antiderivative(SimplifyIntegrand(u / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && PolyQ(~Px, ~x)
    #= Antiderivative((~Px)*((~a')+(~b')*acosh((~c')*(~x)))^(~n'),~x) :=  With([u=IntHide(Px,x)], Dist((a+b*acosh(c*x))^n,u,x) - b*c*n*sqrt(1-c^2*x^2)/(sqrt(-1+c*x)*sqrt(1+c*x))*Antiderivative( SimplifyIntegrand(u*(a+b*acosh(c*x))^(n-1)/sqrt(1-c^2*x^2),x),x)) <--  FreeQ([a,b,c],x) && PolyQ(Px,x) && IGtQ(n,0) =#
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~n, ~x) => Antiderivative(ExpandIntegrand(~Px * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && PolyQ(~Px, ~x)
    @apply_utils Antiderivative(~Px * (~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * acosh(~(c') * ~x)), ~x) => With([u = IntHide(~Px * (~d + ~e * ~x) ^ ~m, ~x)], Dist(~a + ~b * acosh(~c * ~x), u, ~x) - ((~b * ~c * sqrt(1 - (~c) ^ 2 * (~x) ^ 2)) / (sqrt(-1 + ~c * ~x) * sqrt(1 + ~c * ~x))) * Antiderivative(SimplifyIntegrand(u / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && PolyQ(~Px, ~x)
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(p') * (~d + ~(e') * ~x) ^ ~m * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~n, ~x) => With([u = IntHide((~f + ~g * ~x) ^ ~p * (~d + ~e * ~x) ^ ~m, ~x)], Dist((~a + ~b * acosh(~c * ~x)) ^ ~n, u, ~x) - ~b * ~c * ~n * Antiderivative(SimplifyIntegrand((u * (~a + ~b * acosh(~c * ~x)) ^ (~n - 1)) / (sqrt(-1 + ~c * ~x) * sqrt(1 + ~c * ~x)), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (IGtQ(~n, 0) && (IGtQ(~p, 0) && (ILtQ(~m, 0) && LtQ(~m + ~p + 1, 0))))
    @apply_utils Antiderivative(((~(f') + ~(g') * ~x + ~(h') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~n) / (~d + ~(e') * ~x) ^ 2, ~x) => With([u = IntHide((~f + ~g * ~x + ~h * (~x) ^ 2) ^ ~p / (~d + ~e * ~x) ^ 2, ~x)], Dist((~a + ~b * acosh(~c * ~x)) ^ ~n, u, ~x) - ~b * ~c * ~n * Antiderivative(SimplifyIntegrand((u * (~a + ~b * acosh(~c * ~x)) ^ (~n - 1)) / (sqrt(-1 + ~c * ~x) * sqrt(1 + ~c * ~x)), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (IGtQ(~n, 0) && (IGtQ(~p, 0) && EqQ(~e * ~g - 2 * ~d * ~h, 0)))
    @apply_utils Antiderivative(~Px * (~d + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand(~Px * (~d + ~e * ~x) ^ ~m * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (PolyQ(~Px, ~x) && (IGtQ(~n, 0) && IntegerQ(~m)))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => (((-(~d)) ^ IntPart(~p) * (~d + ~e * (~x) ^ 2) ^ FracPart(~p)) / ((-1 + ~c * ~x) ^ FracPart(~p) * (1 + ~c * ~x) ^ FracPart(~p))) * Antiderivative((~f + ~g * ~x) ^ ~m * (-1 + ~c * ~x) ^ ~p * (1 + ~c * ~x) ^ ~p * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~p - 1 / 2) && IntegerQ(~m)))
    @apply_utils Antiderivative(log(~(h') * (~(f') + ~(g') * ~x) ^ ~(m')) * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => (((-(~d)) ^ IntPart(~p) * (~d + ~e * (~x) ^ 2) ^ FracPart(~p)) / ((-1 + ~c * ~x) ^ FracPart(~p) * (1 + ~c * ~x) ^ FracPart(~p))) * Antiderivative(log(~h * (~f + ~g * ~x) ^ ~m) * (-1 + ~c * ~x) ^ ~p * (1 + ~c * ~x) ^ ~p * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && IntegerQ(~p - 1 / 2))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d1 + ~(e1') * ~x) ^ ~p * (~d2 + ~(e2') * ~x) ^ ~p * (~(a') + ~(b') * acosh(~(c') * ~x)), ~x) => With([u = IntHide((~f + ~g * ~x) ^ ~m * (~d1 + ~e1 * ~x) ^ ~p * (~d2 + ~e2 * ~x) ^ ~p, ~x)], Dist(~a + ~b * acosh(~c * ~x), u, ~x) - ~b * ~c * Antiderivative(Dist(1 / (sqrt(-1 + ~c * ~x) * sqrt(1 + ~c * ~x)), u, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~f, ~g], ~x) && (EqQ(~e1 - ~c * ~d1, 0) && (EqQ(~e2 + ~c * ~d2, 0) && (IGtQ(~m, 0) && (ILtQ(~p + 1 / 2, 0) && (GtQ(~d1, 0) && (LtQ(~d2, 0) && (GtQ(~m, 3) || LtQ(~m, -2 * ~p - 1))))))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d1 + ~(e1') * ~x) ^ ~p * (~d2 + ~(e2') * ~x) ^ ~p * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~d1 + ~e1 * ~x) ^ ~p * (~d2 + ~e2 * ~x) ^ ~p * (~a + ~b * acosh(~c * ~x)) ^ ~n, (~f + ~g * ~x) ^ ~m, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~f, ~g], ~x) && (EqQ(~e1 - ~c * ~d1, 0) && (EqQ(~e2 + ~c * ~d2, 0) && (IGtQ(~m, 0) && (IntegerQ(~p + 1 / 2) && (GtQ(~d1, 0) && (LtQ(~d2, 0) && (IGtQ(~n, 0) && (EqQ(~n, 1) && GtQ(~p, -1) || (GtQ(~p, 0) || (EqQ(~m, 1) || EqQ(~m, 2) && LtQ(~p, -2)))))))))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~m * sqrt(~d1 + ~(e1') * ~x) * sqrt(~d2 + ~(e2') * ~x) * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => ((~f + ~g * ~x) ^ ~m * (~d1 * ~d2 + ~e1 * ~e2 * (~x) ^ 2) * (~a + ~b * acosh(~c * ~x)) ^ (~n + 1)) / (~b * ~c * sqrt(-(~d1) * ~d2) * (~n + 1)) - (1 / (~b * ~c * sqrt(-(~d1) * ~d2) * (~n + 1))) * Antiderivative((~d1 * ~d2 * ~g * ~m + 2 * ~e1 * ~e2 * ~f * ~x + ~e1 * ~e2 * ~g * (~m + 2) * (~x) ^ 2) * (~f + ~g * ~x) ^ (~m - 1) * (~a + ~b * acosh(~c * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~f, ~g], ~x) && (EqQ(~e1 - ~c * ~d1, 0) && (EqQ(~e2 + ~c * ~d2, 0) && (ILtQ(~m, 0) && (GtQ(~d1, 0) && (LtQ(~d2, 0) && IGtQ(~n, 0))))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d1 + ~(e1') * ~x) ^ ~p * (~d2 + ~(e2') * ~x) ^ ~p * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand(sqrt(~d1 + ~e1 * ~x) * sqrt(~d2 + ~e2 * ~x) * (~a + ~b * acosh(~c * ~x)) ^ ~n, (~f + ~g * ~x) ^ ~m * (~d1 + ~e1 * ~x) ^ (~p - 1 / 2) * (~d2 + ~e2 * ~x) ^ (~p - 1 / 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~f, ~g], ~x) && (EqQ(~e1 - ~c * ~d1, 0) && (EqQ(~e2 + ~c * ~d2, 0) && (IntegerQ(~m) && (IGtQ(~p + 1 / 2, 0) && (GtQ(~d1, 0) && (LtQ(~d2, 0) && IGtQ(~n, 0)))))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d1 + ~(e1') * ~x) ^ ~p * (~d2 + ~(e2') * ~x) ^ ~p * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => ((~f + ~g * ~x) ^ ~m * (~d1 + ~e1 * ~x) ^ (~p + 1 / 2) * (~d2 + ~e2 * ~x) ^ (~p + 1 / 2) * (~a + ~b * acosh(~c * ~x)) ^ (~n + 1)) / (~b * ~c * sqrt(-(~d1) * ~d2) * (~n + 1)) - (1 / (~b * ~c * sqrt(-(~d1) * ~d2) * (~n + 1))) * Antiderivative(ExpandIntegrand((~f + ~g * ~x) ^ (~m - 1) * (~a + ~b * acosh(~c * ~x)) ^ (~n + 1), (~d1 * ~d2 * ~g * ~m + ~e1 * ~e2 * ~f * (2 * ~p + 1) * ~x + ~e1 * ~e2 * ~g * (~m + 2 * ~p + 1) * (~x) ^ 2) * (~d1 + ~e1 * ~x) ^ (~p - 1 / 2) * (~d2 + ~e2 * ~x) ^ (~p - 1 / 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~f, ~g], ~x) && (EqQ(~e1 - ~c * ~d1, 0) && (EqQ(~e2 + ~c * ~d2, 0) && (ILtQ(~m, 0) && (IGtQ(~p - 1 / 2, 0) && (GtQ(~d1, 0) && (LtQ(~d2, 0) && IGtQ(~n, 0)))))))
    @apply_utils Antiderivative(((~f + ~(g') * ~x) ^ ~(m') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~n) / (sqrt(~d1 + ~(e1') * ~x) * sqrt(~d2 + ~(e2') * ~x)), ~x) => ((~f + ~g * ~x) ^ ~m * (~a + ~b * acosh(~c * ~x)) ^ (~n + 1)) / (~b * ~c * sqrt(-(~d1) * ~d2) * (~n + 1)) - ((~g * ~m) / (~b * ~c * sqrt(-(~d1) * ~d2) * (~n + 1))) * Antiderivative((~f + ~g * ~x) ^ (~m - 1) * (~a + ~b * acosh(~c * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~f, ~g], ~x) && (EqQ(~e1 - ~c * ~d1, 0) && (EqQ(~e2 + ~c * ~d2, 0) && (IGtQ(~m, 0) && (GtQ(~d1, 0) && (LtQ(~d2, 0) && LtQ(~n, -1))))))
    @apply_utils Antiderivative(((~f + ~(g') * ~x) ^ ~(m') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n')) / (sqrt(~d1 + ~(e1') * ~x) * sqrt(~d2 + ~(e2') * ~x)), ~x) => (1 / ((~c) ^ (~m + 1) * sqrt(-(~d1) * ~d2))) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * (~c * ~f + ~g * cosh(~x)) ^ ~m, ~x), ~x, acosh(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~f, ~g, ~n], ~x) && (EqQ(~e1 - ~c * ~d1, 0) && (EqQ(~e2 + ~c * ~d2, 0) && (IntegerQ(~m) && (GtQ(~d1, 0) && (LtQ(~d2, 0) && (GtQ(~m, 0) || IGtQ(~n, 0)))))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d1 + ~(e1') * ~x) ^ ~p * (~d2 + ~(e2') * ~x) ^ ~p * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~a + ~b * acosh(~c * ~x)) ^ ~n / (sqrt(~d1 + ~e1 * ~x) * sqrt(~d2 + ~e2 * ~x)), (~f + ~g * ~x) ^ ~m * (~d1 + ~e1 * ~x) ^ (~p + 1 / 2) * (~d2 + ~e2 * ~x) ^ (~p + 1 / 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~f, ~g], ~x) && (EqQ(~e1 - ~c * ~d1, 0) && (EqQ(~e2 + ~c * ~d2, 0) && (IntegerQ(~m) && (ILtQ(~p + 1 / 2, 0) && (GtQ(~d1, 0) && (LtQ(~d2, 0) && IGtQ(~n, 0)))))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d1 + ~(e1') * ~x) ^ ~p * (~d2 + ~(e2') * ~x) ^ ~p * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => (((-(~d1) * ~d2) ^ IntPart(~p) * (~d1 + ~e1 * ~x) ^ FracPart(~p) * (~d2 + ~e2 * ~x) ^ FracPart(~p)) / ((-1 + ~c * ~x) ^ FracPart(~p) * (1 + ~c * ~x) ^ FracPart(~p))) * Antiderivative((~f + ~g * ~x) ^ ~m * (-1 + ~c * ~x) ^ ~p * (1 + ~c * ~x) ^ ~p * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~f, ~g, ~n], ~x) && (EqQ(~e1 - ~c * ~d1, 0) && (EqQ(~e2 + ~c * ~d2, 0) && (IntegerQ(~m) && (IntegerQ(~p - 1 / 2) && Not(GtQ(~d1, 0) && LtQ(~d2, 0))))))
    @apply_utils Antiderivative((log(~(h') * (~(f') + ~(g') * ~x) ^ ~(m')) * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n')) / (sqrt(~d1 + ~(e1') * ~x) * sqrt(~d2 + ~(e2') * ~x)), ~x) => (log(~h * (~f + ~g * ~x) ^ ~m) * (~a + ~b * acosh(~c * ~x)) ^ (~n + 1)) / (~b * ~c * sqrt(-(~d1) * ~d2) * (~n + 1)) - ((~g * ~m) / (~b * ~c * sqrt(-(~d1) * ~d2) * (~n + 1))) * Antiderivative((~a + ~b * acosh(~c * ~x)) ^ (~n + 1) / (~f + ~g * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~f, ~g, ~h, ~m], ~x) && (EqQ(~e1 - ~c * ~d1, 0) && (EqQ(~e2 + ~c * ~d2, 0) && (GtQ(~d1, 0) && (LtQ(~d2, 0) && IGtQ(~n, 0)))))
    @apply_utils Antiderivative(log(~(h') * (~(f') + ~(g') * ~x) ^ ~(m')) * (~d1 + ~(e1') * ~x) ^ ~p * (~d2 + ~(e2') * ~x) ^ ~p * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => (((-(~d1) * ~d2) ^ IntPart(~p) * (~d1 + ~e1 * ~x) ^ FracPart(~p) * (~d2 + ~e2 * ~x) ^ FracPart(~p)) / ((-1 + ~c * ~x) ^ FracPart(~p) * (1 + ~c * ~x) ^ FracPart(~p))) * Antiderivative(log(~h * (~f + ~g * ~x) ^ ~m) * (-1 + ~c * ~x) ^ ~p * (1 + ~c * ~x) ^ ~p * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~f, ~g, ~h, ~m, ~n], ~x) && (EqQ(~e1 - ~c * ~d1, 0) && (EqQ(~e2 + ~c * ~d2, 0) && (IntegerQ(~p - 1 / 2) && Not(GtQ(~d1, 0) && LtQ(~d2, 0)))))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~m * (~f + ~(g') * ~x) ^ ~m * (~(a') + ~(b') * acosh(~(c') * ~x)), ~x) => With([u = IntHide((~d + ~e * ~x) ^ ~m * (~f + ~g * ~x) ^ ~m, ~x)], Dist(~a + ~b * acosh(~c * ~x), u, ~x) - ~b * ~c * Antiderivative(Dist(1 / (sqrt(-1 + ~c * ~x) * sqrt(1 + ~c * ~x)), u, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && ILtQ(~m + 1 / 2, 0)
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * (~f + ~(g') * ~x) ^ ~(m') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~a + ~b * acosh(~c * ~x)) ^ ~n, (~d + ~e * ~x) ^ ~m * (~f + ~g * ~x) ^ ~m, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && IntegerQ(~m)
    @apply_utils Antiderivative(~u * (~(a') + ~(b') * acosh(~(c') * ~x)), ~x) => With([v = IntHide(~u, ~x)], Dist(~a + ~b * acosh(~c * ~x), v, ~x) - ((~b * ~c * sqrt(1 - (~c) ^ 2 * (~x) ^ 2)) / (sqrt(-1 + ~c * ~x) * sqrt(1 + ~c * ~x))) * Antiderivative(SimplifyIntegrand(v / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x), ~x) <-- InverseFunctionFreeQ(v, ~x)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(~Px * (~d1 + ~(e1') * ~x) ^ ~p * (~d2 + ~(e2') * ~x) ^ ~p * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~n, ~x) => With([u = ExpandIntegrand(~Px * (~d1 + ~e1 * ~x) ^ ~p * (~d2 + ~e2 * ~x) ^ ~p * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~n], ~x) && (PolyQ(~Px, ~x) && (EqQ(~e1 - ~c * ~d1, 0) && (EqQ(~e2 + ~c * ~d2, 0) && IntegerQ(~p - 1 / 2))))
    @apply_utils Antiderivative(~(Px') * (~f + ~(g') * (~d1 + ~(e1') * ~x) ^ ~p * (~d2 + ~(e2') * ~x) ^ ~p) ^ ~(m') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => With([u = ExpandIntegrand(~Px * (~f + ~g * (~d1 + ~e1 * ~x) ^ ~p * (~d2 + ~e2 * ~x) ^ ~p) ^ ~m * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~f, ~g], ~x) && (PolyQ(~Px, ~x) && (EqQ(~e1 - ~c * ~d1, 0) && (EqQ(~e2 + ~c * ~d2, 0) && (IGtQ(~p + 1 / 2, 0) && IntegersQ(~m, ~n)))))
    @apply_utils Antiderivative(~RFx * acosh(~(c') * ~x) ^ ~(n'), ~x) => With([u = ExpandIntegrand(acosh(~c * ~x) ^ ~n, ~RFx, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ(~c, ~x) && (RationalFunctionQ(~RFx, ~x) && IGtQ(~n, 0))
    @apply_utils Antiderivative(~RFx * (~a + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand(~RFx * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (RationalFunctionQ(~RFx, ~x) && IGtQ(~n, 0))
    @apply_utils Antiderivative(~RFx * (~d1 + ~(e1') * ~x) ^ ~p * (~d2 + ~(e2') * ~x) ^ ~p * acosh(~(c') * ~x) ^ ~(n'), ~x) => With([u = ExpandIntegrand((~d1 + ~e1 * ~x) ^ ~p * (~d2 + ~e2 * ~x) ^ ~p * acosh(~c * ~x) ^ ~n, ~RFx, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ([~c, ~d1, ~e1, ~d2, ~e2], ~x) && (RationalFunctionQ(~RFx, ~x) && (IGtQ(~n, 0) && (EqQ(~e1 - ~c * ~d1, 0) && (EqQ(~e2 + ~c * ~d2, 0) && IntegerQ(~p - 1 / 2)))))
    @apply_utils Antiderivative(~RFx * (~d1 + ~(e1') * ~x) ^ ~p * (~d2 + ~(e2') * ~x) ^ ~p * (~a + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~d1 + ~e1 * ~x) ^ ~p * (~d2 + ~e2 * ~x) ^ ~p, ~RFx * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2], ~x) && (RationalFunctionQ(~RFx, ~x) && (IGtQ(~n, 0) && (EqQ(~e1 - ~c * ~d1, 0) && (EqQ(~e2 + ~c * ~d2, 0) && IntegerQ(~p - 1 / 2)))))
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => Unintegrable(~u * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
end

max_stars_repo_path: src/rules/7 Inverse hyperbolic functions/7.2 Inverse hyperbolic cosine/7.2.5 u (a+b arccosh(c x))^n.jl


text:
module GridapDistributedTests

include("sequential/runtests.jl")

include("mpi/runtests.jl")

end # module


text:
# TemplateParameterList
function getNumTemplateParameterLists(x::AbstractTagDecl)
    @check_ptrs x
    return clang_TypeDecl_getNumTemplateParameterLists(x)
end

function getTemplateParameterList(x::AbstractTagDecl, i::Integer)
    @check_ptrs x
    return TemplateParameterList(clang_TypeDecl_getTemplateParameterList(x, i))
end

function getParam(x::TemplateParameterList, i::Integer)
    @check_ptrs x
    return NamedDecl(clang_TemplateParameterList_getParam(x, i))
end

function Base.size(x::TemplateParameterList)
    @check_ptrs x
    return clang_TemplateParameterList_size(x)
end

function getDescribedTemplateParams(x::AbstractDecl)
    @check_ptrs x
    return TemplateParameterList(clang_Decl_getDescribedTemplateParams(x))
end

# TemplateArgumentList
function TemplateArgumentList(ctx::ASTContext, args::Vector{CXTemplateArgument})
    @check_ptrs ctx
    list = clang_TemplateArgumentList_CreateCopy(ctx, args, length(args))
    return TemplateArgumentList(list)
end

function TemplateArgumentList(ctx::ASTContext, args::Vector{TemplateArgument})
    return TemplateArgumentList(ctx, [arg.ptr for arg in args])
end

function Base.size(x::TemplateArgumentList)
    @check_ptrs x
    return clang_TemplateArgumentList_size(x)
end

function data(x::TemplateArgumentList)
    @check_ptrs x
    return clang_TemplateArgumentList_data(x)
end

function get(x::TemplateArgumentList, i::Integer)
    @check_ptrs x
    return clang_TemplateArgumentList_get(x, i)
end

# TemplateDecl
function init(x::AbstractTemplateDecl, nd::NamedDecl, tp::TemplateParameterList)
    @check_ptrs x nd tp
    return clang_TemplateDecl_init(x, nd, tp)
end

getAsTemplateDecl(x::TemplateName) = TemplateDecl(clang_TemplateName_getAsTemplateDecl(x))

# RedeclarableTemplateDecl
function getCanonicalDecl(x::AbstractRedeclarableTemplateDecl)
    @check_ptrs x
    return RedeclarableTemplateDecl(clang_RedeclarableTemplateDecl_getCanonicalDecl(x))
end

function isMemberSpecialization(x::AbstractRedeclarableTemplateDecl)
    @check_ptrs x
    return clang_RedeclarableTemplateDecl_isMemberSpecialization(x)
end

function setMemberSpecialization(x::AbstractRedeclarableTemplateDecl)
    @check_ptrs x
    return clang_RedeclarableTemplateDecl_setMemberSpecialization(x)
end

function getTemplatedDecl(x::AbstractRedeclarableTemplateDecl)
    @check_ptrs x
    return CXXRecordDecl(clang_ClassTemplateDecl_getTemplatedDecl(x))
end

function isThisDeclarationADefinition(x::AbstractRedeclarableTemplateDecl)
    @check_ptrs x
    return clang_ClassTemplateDecl_isThisDeclarationADefinition(x)
end

function findSpecialization(x::AbstractRedeclarableTemplateDecl, list::TemplateArgumentList,
                            insert_pos=C_NULL)
    @check_ptrs x list
    ctsd = clang_ClassTemplateDecl_findSpecialization(x, list, insert_pos)
    return ClassTemplateSpecializationDecl(ctsd)
end

# ClassTemplateDecl
function getCanonicalDecl(x::AbstractClassTemplateDecl)
    @check_ptrs x
    return ClassTemplateDecl(clang_ClassTemplateDecl_getCanonicalDecl(x))
end

function getTemplateArgs(x::AbstractClassTemplateDecl)
    @check_ptrs x
    return TemplateArgumentList(clang_ClassTemplateSpecializationDecl_getTemplateArgs(x))
end

function setTemplateArgs(x::AbstractClassTemplateDecl, list::TemplateArgumentList)
    @check_ptrs x list
    return clang_ClassTemplateSpecializationDecl_setTemplateArgs(x, list)
end

# ClassTemplateSpecializationDecl
function ClassTemplateSpecializationDecl(ctx::ASTContext, tk::CXTagTypeKind,
                                         dc::DeclContext, start_loc::SourceLocation,
                                         id_loc::SourceLocation,
                                         template::ClassTemplateDecl,
                                         args::TemplateArgumentList,
                                         prev_decl::ClassTemplateSpecializationDecl)
    @check_ptrs ctx dc template args
    ctsd = clang_ClassTemplateSpecializationDecl_Create(ctx, tk, dc, start_loc, id_loc,
                                                        template, args, prev_decl)
    return ClassTemplateSpecializationDecl(ctsd)
end

function ClassTemplateSpecializationDecl(ctx::ASTContext, template::ClassTemplateDecl,
                                         args::TemplateArgumentList,
                                         prev_decl::ClassTemplateSpecializationDecl=ClassTemplateSpecializationDecl(C_NULL))
    tdecl = getTemplatedDecl(template)
    tk = getTagKind(tdecl)
    dc_ctx = getDeclContext(template)
    start_loc = getBeginLoc(tdecl)
    id_loc = getLocation(template)
    return ClassTemplateSpecializationDecl(ctx, tk, dc_ctx, start_loc, id_loc, template,
                                           args, prev_decl)
end

function AddSpecialization(x::AbstractRedeclarableTemplateDecl,
                           ctsd::ClassTemplateSpecializationDecl, insert_pos=C_NULL)
    @check_ptrs x ctsd
    return clang_ClassTemplateDecl_AddSpecialization(x, ctsd, insert_pos)
end


text:
using Test
using CLIMA
using CLIMA.RootSolvers

@testset "RootSolvers - compact solution correctness" begin
  f(x) = x^2 - 100^2
  f′(x) = 2x
  for FT in [Float32, Float64]
    sol = find_zero(f, FT(0.0), FT(1000.0), SecantMethod(), CompactSolution())
    @test sol.converged
    @test sol.root isa FT
    @test sol.root ≈ 100

    sol = find_zero(f, FT(0.0), FT(1000.0), RegulaFalsiMethod(), CompactSolution())
    @test sol.converged
    @test sol.root isa FT
    @test sol.root ≈ 100

    sol = find_zero(f, FT(1.0), NewtonsMethodAD(), CompactSolution())
    @test sol.converged
    @test sol.root isa FT
    @test sol.root ≈ 100

    sol = find_zero(f, f′, FT(1.0), NewtonsMethod(), CompactSolution())
    @test sol.converged
    @test sol.root isa FT
    @test sol.root ≈ 100
  end
end

@testset "RootSolvers - verbose solution correctness" begin
  f(x) = x^2 - 100^2
  f′(x) = 2x
  for FT in [Float32, Float64]
    sol = find_zero(f, FT(0.0), FT(1000.0), SecantMethod(), VerboseSolution())
    @test sol.converged
    @test sol.root isa FT
    @test sol.root ≈ 100
    @test sol.err < 1e-3
    @test sol.iter_performed < 20
    @test sol.iter_performed+1 == length(sol.root_history) == length(sol.err_history)

    sol = find_zero(f, FT(0.0), FT(1000.0), RegulaFalsiMethod(), VerboseSolution())
    @test sol.converged
    @test sol.root isa FT
    @test sol.root ≈ 100
    @test sol.err < 1e-3
    @test sol.iter_performed < 20
    @test sol.iter_performed+1 == length(sol.root_history) == length(sol.err_history)

    sol = find_zero(f, FT(1.0), NewtonsMethodAD(), VerboseSolution())
    @test sol.converged
    @test sol.root isa FT
    @test sol.root ≈ 100
    @test sol.err < 1e-3
    @test sol.iter_performed < 20
    @test sol.iter_performed+1 == length(sol.root_history) == length(sol.err_history)

    sol = find_zero(f, f′, FT(1.0), NewtonsMethod(), VerboseSolution())
    @test sol.converged
    @test sol.root isa FT
    @test sol.root ≈ 100
    @test sol.err < 1e-3
    @test sol.iter_performed < 20
    @test sol.iter_performed+1 == length(sol.root_history) == length(sol.err_history)
  end
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
# verify processing of input files works
"""
  This function tests whether the known keys list works correctly
"""
function test_input()
  @testset "---- Testing input processing ------" begin

    extract_path = joinpath(Pkg.dir("PDESolver"), "src/input/extract_keys.jl")

    include(extract_path)
#    include("known_keys.jl")
    @test ( haskey(Input.KNOWN_KEYS, "key1") )== false
    @test ( haskey(Input.KNOWN_KEYS, "smb_name") )== true
    @test ( haskey(Input.KNOWN_KEYS, "var1") )== false

    arg_dict = Input.read_input_file("input_test.jl")

    @test ( Input.checkKeys(arg_dict, Input.KNOWN_KEYS) )== 1

    arg_dict["numBC"] = 1
    arg_dict["BC1"] = [0, 0]

    @test_throws Exception  Inputs.checkBCOptions(arg_dict)

    # pick another random input file
    opts = Input.read_input("input_vals_channel.jl")

    # verify that supplying default values is idempotent
    opts_orig = copy(opts)
    Input.read_input(opts)

    for (key, val) in opts
      @test ( val )== opts_orig[key]
    end

  end
end

#test_input()
add_func1!(EulerTests, test_input, [TAG_INPUT, TAG_SHORTTEST])


text:
using AtomsBase
using Documenter

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

makedocs(;
    modules=[AtomsBase],
    authors="JuliaMolSim community",
    repo="https://github.com/JuliaMolSim/AtomsBase.jl/blob/{commit}{path}#{line}",
    sitename="AtomsBase.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://juliamolsim.github.io/AtomsBase.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "overview.md",
        "atomicsystems.md",
        "apireference.md"
    ],
    checkdocs=:exports,
    strict=true,
)

deploydocs(;
    repo="github.com/JuliaMolSim/AtomsBase.jl",
    devbranch="master",
)


text:
using Documenter
using HASlib

makedocs(
    sitename = "HASlib",
    modules = [HASlib],
    format = Documenter.HTML(prettyurls = false)
)

deploydocs(
    repo = "github.com/feanor12/HASlib.jl.git",
)


text:
# Plot solutions to 2D Fisher-Stefan
# Alex Tam, 09/11/2021

using Plots
using Measures
using LaTeXStrings
using DelimitedFiles

"Plot solutions as 2D heat maps"
function draw_heat(x, y, U, V, ϕ, i, Lx, Ly)
    gr(); plot() # Load GR plotting backend and clear previous plots
    default(titlefont = (18, "Computer Modern"), guidefont = (26, "Computer Modern"), tickfont = (18, "Computer Modern"))
    heatmap(x,y,transpose(U), xlabel = L"$x$", ylabel = L"$y$", margin=3mm, aspect_ratio=:equal, xlims=(0,Lx), ylims =(0,Ly), tick_direction=:out, c=:plasma, clims=(0.0, 1.0))
    savefig("u-$i.pdf")
    heatmap(x,y,transpose(V), xlabel = L"$x$", ylabel = L"$y$", margin=3mm, aspect_ratio=:equal, xlims=(0,Lx), ylims =(0,Ly), tick_direction=:out, c=:plasma)
    savefig("V-$i.pdf")
    heatmap(x,y,transpose(ϕ), xlabel = L"$x$", ylabel = L"$y$", margin=3mm, aspect_ratio=:equal, xlims=(0,Lx), ylims =(0,Ly), tick_direction=:out, c=:plasma)
    savefig("phi-$i.pdf")
end

function draw_heat(x, y, U, ϕ, i, Lx, Ly)
    gr(); plot() # Load GR plotting backend and clear previous plots
    default(titlefont = (18, "Computer Modern"), guidefont = (26, "Computer Modern"), tickfont = (18, "Computer Modern"))
    heatmap(x,y,transpose(U), xlabel = L"$x$", ylabel = L"$y$", margin=3mm, aspect_ratio=:equal, xlims=(0,Lx), ylims =(0,Ly), tick_direction=:out, c=:plasma, clims=(0.0, 1.0))
    savefig("u-$i.pdf")
    heatmap(x,y,transpose(ϕ), xlabel = L"$x$", ylabel = L"$y$", margin=3mm, aspect_ratio=:equal, xlims=(0,Lx), ylims =(0,Ly), tick_direction=:out, c=:plasma)
    savefig("phi-$i.pdf")
end

"Density slice plots"
function draw_slices(x, y, nx, ny, Lx, Ly, plot_times)
    gr(); plot() # Load GR plotting backend and clear previous plots
    default(titlefont = (18, "Computer Modern"), guidefont = (26, "Computer Modern"), tickfont = (18, "Computer Modern"))
    U = readdlm("U-0.csv")
    plot(x, U[:,ny], xlabel = L"$x$", ylabel = L"$u(x,10,t)$", linecolor = :black, linewidth = 2, aspect_ratio = 20.0, grid = false, margin=3mm, legend = false, xlims=(0,Lx), ylims=(0,1))
    for i in plot_times[2:2:end]
        u = readdlm("ux-$i.csv")
        plot!(x, u, linecolor = :black, linestyle = :dash, linewidth = 2)
    end
    # plot!([7.595174442,7.595174442], [0,1], linecolor = :red, linewidth = 2)
    # plot!([12.404825558,12.404825558], [0,1], linecolor = :red, linewidth = 2)
    savefig("u_slice_x.pdf")
    plot(y, U[nx,:], xlabel = L"$y$", ylabel = L"$u(10,y,t)$", linecolor = :black, linewidth = 2, aspect_ratio = 20.0, grid = false, margin=3mm, legend = false, xlims=(0,Ly), ylims=(0,1))
    for i in plot_times[2:2:end]
        u = readdlm("uy-$i.csv")
        plot!(y, u, linecolor = :black, linestyle = :dash, linewidth = 2)
    end
    # plot!([7.595174442,7.595174442], [0,1], linecolor = :red, linewidth = 2)
    # plot!([12.404825558,12.404825558], [0,1], linecolor = :red, linewidth = 2)
    savefig("u_slice_y.pdf")
end

# Control function for plotting
function draw()
    plot_times = convert(Vector{Int}, vec(readdlm("plot_times.csv")))
    x = vec(readdlm("x.csv"))
    y = vec(readdlm("y.csv"))
    nx::Int = (length(x)-1)/2; ny::Int = (length(y)-1)/2
    Lx = maximum(x); Ly = maximum(y)
    U = readdlm("U-0.csv")
    ϕ = readdlm("Phi-0.csv")
    draw_heat(x, y, U, ϕ, 0, Lx, Ly)
    draw_slices(x, y, nx, ny, Lx, Ly, plot_times)
    for i in plot_times
        U = readdlm("U-$i.csv")
        V = readdlm("V-$i.csv")
        ϕ = readdlm("Phi-$i.csv")
        draw_heat(x, y, U, V, ϕ, i, Lx, Ly)
    end
end

@time draw()


text:
"""
    B,T=hkz(H)

Do Hermite-Korkine-Zolotarev (HKZ) reduction of the basis `H`, returning
the reduced basis `B` and the unimodular rotation `T` such that B = H*T.
HKZ reduction is sometimes called "Hermite-Korkine-Zolotareff" or
"Korkine-Zolotareff" reduction.

Based on "Practical HKZ and Minkowski Lattice Reduction Algorithms" by Wen
Zhang, Sanzheng Qiao, and Yimin Wei, 17 Aug. 2011.
http://www.cas.mcmaster.ca/~qiao/publications/ZQW11.pdf

# Examples
```jldoctest
julia> H=[1 2; 3 4];

julia> B,T= hkz(H); B
2×2 Matrix{Int64}:
  1  -1
 -1  -1

julia> N=9; H=rand(0:10,N,N); B,_=hkz(H);

julia> ishkzreduced(B)
true

julia> A = [10.6347 -66.2715  9.3046 17.5349 24.9625 # From
                0     8.6759 -4.7536 -3.9379 -2.3318 # TransIT, v65,
                0        0    0.3876  0.1296 -0.2879 # n3, Mar 2019,
                0        0    0       0.0133 -0.0082 # p 1929.
                0        0    0       0       0.0015];

julia> B,_= hkz(A); # our HKZ doesn't work well on A ...

julia> ishkzreduced(B) # ...because A is ill-conditioned
false

```
"""
function hkz(B::AbstractArray{Td,2}) where {Td<:Number}
    Z,_ = hkz_red(B)
    return B*Z,Z
end

"""
    x=hkz_red(B)

Do Hermite-Korkine-Zolotareff reduction of the basis B.

This is the "HKZ-red" function for HKZ reduction in
"Practical HKZ and Minkowski Lattice Reduction Algorithms" by
Wen Zhang, Sanzheng Qiao, and Yimin Wei, 17 Aug. 2011.
http://www.cas.mcmaster.ca/~qiao/publications/ZQW11.pdf

# Examples
```jldoctest
julia> B=[1 2; 3 4];

julia> Z,_= LLLplus.hkz_red(B);

julia> B*Z
2×2 Matrix{Int64}:
  1  -1
 -1  -1

julia> N=10; Bo=rand(0:100,N,N); Z,_=LLLplus.hkz_red(Bo); B=Bo*Z;

julia> ishkzreduced(B)
true

julia> A = [10.6347 -66.2715  9.3046 17.5349 24.9625 # From
                0     8.6759 -4.7536 -3.9379 -2.3318 # TransIT, v65,
                0        0    0.3876  0.1296 -0.2879 # n3, Mar 2019,
                0        0    0       0.0133 -0.0082 # p 1929.
                0        0    0       0       0.0015];

julia> Z,_= LLLplus.hkz_red(A); # our HKZ doesn't work well on A ...

julia> ishkzreduced(A*Z)  # ...because A is ill-conditioned
false

```
"""
function hkz_red(B::AbstractArray{Td,2}) where {Td<:Number}
    m,n = size(B)
    Ti= getIntType(Td)

    #Q,R = qr(B)  # QR doesn't give positive diagonal elements
    # HKZ-Red doesn't specify positive diagonal, but elsewhere in the
    # paper it does for QR, so I'll use it
    R = Matrix(cholesky(B'*B).U)
    Z = Matrix{Ti}(I, n,n)
    for k=1:n-1
        z = svpu(R[k:n,k:n])
        transform!(R,Z,z,k)
    end
    # size reduce
    for lx=2:n
        # reduce lx-th column of B
        for k=lx-1:-1:1
            rk = R[k,lx]/R[k,k]
            mu = round(rk)
            if abs(mu)>0
                R[1:k,lx] .-= mu .* view(R,1:k,k)
                Z[:,lx]   .-= mu .* view(Z,:,k)
            end
        end
    end
    return Z,R
end

"""
    transform!(R,Z,z,k)

Basis expansion. Given n x n upper triangular R, n x n unimodular Z, n-k+1
vector of integer z, and index k, return updated R and Z such that
R[k,k] = ||R[k:n,k:n]*z||.

This is the "TRANSFORM" basis expansion function from
"Practical HKZ and Minkowski Lattice Reduction Algorithms" by
Wen Zhang, Sanzheng Qiao, and Yimin Wei, 17 Aug. 2011.
http://www.cas.mcmaster.ca/~qiao/publications/ZQW11.pdf

"""
function transform!(R::AbstractArray{Td,2},Z::AbstractArray{<:Integer,2},
                    z::AbstractArray{<:Integer,1},k::Integer) where {Td<:Number}
    n = size(R,1)
    for j=n-k+1:-1:2
        if z[j]≠0
            # d=gcd(z[j-1],z[j]),  a,b are such that a*z[j−1] + b*z[j] = d
            d,a,b=gcdx(z[j-1],z[j])
            M = [z[j-1]÷d -b;
                 z[j]÷d    a]
            z[j-1] = d
            R[1:j+k-1, j+k-2:j+k-1]= R[1:j+k-1, j+k-2:j+k-1]*M
            c,s = givensrotation(R[j+k-2,j+k-2],R[j+k-1,j+k-2])
            G = [c -s; s c]
            R[j+k-2:j+k-1,j+k-2:n] = G'*R[j+k-2:j+k-1,j+k-2:n]
            Z[:,j+k-2:j+k-1] = Z[:,j+k-2:j+k-1]*M
        end
    end
end

"""
    c,s = givensrotation(a,b)
"""
function givensrotation(a,b)
    if b == 0
        c = 1
        s = 0
    else
        if abs(b) > abs(a)
            r = a / b
            s = 1 / sqrt(1 + r^2)
            c = s*r
        else
            r = b / a
            c = 1 / sqrt(1 + r^2)
            s = c*r
        end
    end
    return c,s
end

"""
    ishkzreduced(B)

Determine if the matrix B is Hermite-Korkine-Zolotarev (HKZ) reduced or
not. See the `hkz` function.

# Examples
```jldoctest
julia> H= [1 2; 3 4];ishkzreduced(H)
false

julia> B,_=hkz(H); ishkzreduced(B)
true

```
"""
function ishkzreduced(B)
    if issizereduced(B)
        m,n = size(B)
        R = Matrix(cholesky(B'*B).U)
        for i=1:n
            s = svp(R[i:n,i:n])
            if !(R[i,i]≈sqrt(s'*s))
                return false
            end
        end
        return true
    else
        return false
    end
end


text:
function check_codesigned(path::AbstractString, platform::AbstractPlatform)
    # We only perform ad-hoc codesigning on Apple platforms
    if !Sys.isapple(platform)
        return true
    end

    ur = preferred_runner()(dirname(path); cwd="/workspace/", platform=platform)
    return run(ur, `/usr/local/bin/ldid -d $(basename(path))`)
end

function ensure_codesigned(path::AbstractString, prefix::Prefix, platform::AbstractPlatform;
                           verbose::Bool = false, subdir::AbstractString="")
    # We only perform ad-hoc codesigning on Apple platforms
    if !Sys.isapple(platform)
        return true
    end

    rel_path = relpath(path, prefix.path)
    ur = preferred_runner()(prefix.path; cwd="/workspace/", platform=platform)
    with_logfile(prefix, "ldid_$(basename(rel_path)).log"; subdir) do io
        run(ur, `/usr/local/bin/ldid -S -d $(rel_path)`, io; verbose=verbose)
    end
end


text:
include("imports.jl")

@testset ExtendedTestSet "`methodlearner`" begin
    method = BlockMethod((Label(1:2), Label(1:2)), (OneHot(),))
    data = (rand(1:2, 1000), rand(1:2, 1000))
    @test_nowarn learner = methodlearner(method, data, model=identity)

    @testset ExtendedTestSet "batch sizes" begin
        learner = methodlearner(method, data, model=identity, batchsize=100)
        @test length(learner.data.training) == 8
        @test length(learner.data.validation) == 1

        learner = methodlearner(method, data, model=identity, pctgval=0.4, batchsize=100)
        @test length(learner.data.training) == 6
        @test length(learner.data.validation) == 2

        learner = methodlearner(method, data, model=identity, batchsize=100, validbsfactor=1)
        @test length(learner.data.training) == 8
        @test length(learner.data.validation) == 2
    end

    @testset ExtendedTestSet "callbacks" begin
        learner = methodlearner(
            method, data, model=identity,
            callbacks=[ToGPU(), Checkpointer(mktempdir())])
        @test !isnothing(FluxTraining.getcallback(learner, Checkpointer))

    end
end

@testset ExtendedTestSet "`blockbackbone`" begin
    @test_nowarn FastAI.blockbackbone(FastAI.ImageTensor{2}(3))

    @test_nowarn FastAI.blockbackbone(FastAI.EncodedTableRow((:x,), (:y,), Dict(:x => [1, 2])))
end

@testset ExtendedTestSet "`blockmodel`" begin
    method = ImageClassificationSingle((Image{2}(), Label(1:2)))
    @test_nowarn methodmodel(method)
end


text:
using Test

include("roman-numerals.jl")

samples = Dict(
    1 => "I",
    2 => "II",
    3 => "III",
    4 => "IV",
    5 => "V",
    6 => "VI",
    9 => "IX",
    27 => "XXVII",
    48 => "XLVIII",
    49 => "XLIX",
    59 => "LIX",
    93 => "XCIII",
    141 => "CXLI",
    163 => "CLXIII",
    402 => "CDII",
    575 => "DLXXV",
    911 => "CMXI",
    1024 => "MXXIV",
    1703 => "MDCCIII",
    1991 => "MCMXCI",
    2017 => "MMXVII",
    3000 => "MMM"
)

@testset "convert $sample[1] to roman numeral" for sample in samples
    @test to_roman(sample[1]) == sample[2]
end

@testset "error handling" begin
    @test_throws ErrorException to_roman(0)
    @test_throws ErrorException to_roman(-2017)
end

max_stars_repo_name: tripitakit/julia

max_issues_repo_name: tripitakit/julia

max_forks_repo_name: tripitakit/julia


text:
function findRequiredIdxs(A::LagrangeInterpolation, t)
  idxs = sortperm(A.t,by=x->abs(t-x))
  idxs[1:(A.n+1)]
end

function spline_coefficients(n, d, k, u::Number)
  N = zeros(n)
  if u == k[1]
    N[1] = one(u)
  elseif u == k[end]
    N[end] = one(u)
  else
    i = findfirst(x->x>u,k) - 1
    N[i] = one(u)
    for deg = 1:d
      N[i-deg] = (k[i+1]-u)/(k[i+1]-k[i-deg+1]) * N[i-deg+1]
      for j = (i-deg+1):(i-1)
        N[j] = (u-k[j])/(k[j+deg]-k[j]) * N[j] + (k[j+deg+1]-u)/(k[j+deg+1]-k[j+1]) * N[j+1]
      end
      N[i] = (u-k[i])/(k[i+deg]-k[i]) * N[i]
    end
  end
  N
end

function spline_coefficients(n, d, k, u::AbstractVector)
  N = zeros(eltype(u),n,n)
  for i = 1:n
    N[i,:] .= spline_coefficients(n,d,k,u[i])
  end
  N
end

function weibull_fun(x,p)
  ans = copy(x)
  for i = 1:length(x)
    if x[i] >= 0
      ans[i] = p[1] * (1 - exp(-1 * (x[i]/p[2])^p[3]))
    else
      ans[i] = 0.0
    end
  end
  ans
end

# helper functions to get t_max and c_max (OneCompartmentPK)
function OneCompartmentPK_tmax(A)
  p = A.c_f.param
  2.303 * log(10.0,p[1]/p[4]) / (p[1] - p[4])
end

function OneCompartmentPK_cmax(A)
  p = A.c_f.param
  t = OneCompartmentPK_tmax(A)
  p[1] * p[2] * p[3] * (exp(-p[4]*t) - exp(-p[1]*t)) / (p[5] * (p[1] - p[4]))
end


text:
module LinearOperatorsRevamp

using FastClosures, LinearAlgebra, Printf, SparseArrays

# Basic defitions
include("abstract.jl")
include("constructors.jl")

# Operations
include("operations.jl") # This first
include("adjtrans.jl")
include("cat.jl")
include("kron.jl")

# quasi-Newton operators
include("qn.jl")

# Special operators
include("linalg.jl")
# include("PreallocatedLinearOperators.jl")
include("special-operators.jl")
include("TimedOperators.jl")

# Utilities
include("utilities.jl")

include("deprecated.jl")

end # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Overhead structure represents the resource overhead associated with running a pod.

    IoK8sApiNodeV1beta1Overhead(;
        podFixed=nothing,
    )

    - podFixed::Dict{String, IoK8sApimachineryPkgApiResourceQuantity} : PodFixed represents the fixed resource overhead associated with running a pod.
"""
mutable struct IoK8sApiNodeV1beta1Overhead <: SwaggerModel
    podFixed::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: podFixed

    function IoK8sApiNodeV1beta1Overhead(;podFixed=nothing)
        o = new()
        validate_property(IoK8sApiNodeV1beta1Overhead, Symbol("podFixed"), podFixed)
        setfield!(o, Symbol("podFixed"), podFixed)
        o
    end
end # type IoK8sApiNodeV1beta1Overhead

const _property_map_IoK8sApiNodeV1beta1Overhead = Dict{Symbol,Symbol}(Symbol("podFixed")=>Symbol("podFixed"))
const _property_types_IoK8sApiNodeV1beta1Overhead = Dict{Symbol,String}(Symbol("podFixed")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}")
Base.propertynames(::Type{ IoK8sApiNodeV1beta1Overhead }) = collect(keys(_property_map_IoK8sApiNodeV1beta1Overhead))
Swagger.property_type(::Type{ IoK8sApiNodeV1beta1Overhead }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiNodeV1beta1Overhead[name]))}
Swagger.field_name(::Type{ IoK8sApiNodeV1beta1Overhead }, property_name::Symbol) =  _property_map_IoK8sApiNodeV1beta1Overhead[property_name]

function check_required(o::IoK8sApiNodeV1beta1Overhead)
    true
end

function validate_property(::Type{ IoK8sApiNodeV1beta1Overhead }, name::Symbol, val)
end


text:
module Example1

using Mapper2

# Export many functions for testing
export build_input_primitive,
    build_output_primitive,
    build_general_primitive,
    build_io_tile,
    build_general_tile,
    build_super_tile,
    build_double_general_tile,
    build_routing_tile,
    build_test_arch,
    # Taskgraph
    make_taskgraph,
    make_nodes,
    make_edges,
    # Map
    make_map,
    make_fanout

# Extensions of the base Mapper types
struct TestRuleSet <: RuleSet end

make_map() = Map(TestRuleSet(), build_test_arch(), make_taskgraph())
make_fanout() = Map(TestRuleSet(), build_test_arch(), make_fanout_taskgraph())

include("Architecture.jl")
include("Taskgraph.jl")
include("Placement.jl")
include("Routing.jl")

end


text:
using ApproxFunOrthogonalPolynomials, ApproxFunRational,
 ApproxFunFourier, ApproxFunBase, ApproxFun, AbstractIterativeSolvers, Plots

### Scalar case
L = 1.
g = x -> .2exp(-x^2)
G = Fun(zai(g),OscLaurent(-1.,L),100) + Fun(zai(g),OscLaurent(1.,L),100)
h = x -> .5/(1 + x^2)
H = Fun(zai(h),OscLaurent(-1.,L),100)
𝓒 = Cauchy(-1)
𝓢 = Cauchy(1)

op = x -> condense(x - G*(𝓒*x))
out = GMRES(op,H,⋅,1e-8,100, condense)
u = sum([out[2][i]*out[1][i] for i=1:length(out[2]) ])

1.0 + (𝓢*u)(1.)
1.0 + (𝓒*u)(1.)

jump = x -> (𝓢*u)(x) - (𝓒*u)(x)*G(x) - H(x)
jump(.11)

jump = x -> (𝓢*u)(x) + 1 - ((𝓒*u)(x)+1)*(1 + G(x))
jump(.11)

jump.space

conj(jump).space

jump

𝓕 = FourierTransform(1.0)
U = 𝓕*u

x = 0:.01:10
y = real(map(U,x))
plot(x,y)

res = op(u)-b
sqrt(res⋅res)


text:
"""
$(DocStringExtensions.README)
"""
module Catalyst

using DocStringExtensions
using SparseArrays, DiffEqBase, Reexport
using LaTeXStrings, Latexify, Requires
using DiffEqJump: DiffEqJump, JumpProblem, MassActionJump, ConstantRateJump, VariableRateJump

# ModelingToolkit imports and convenience functions we use
using ModelingToolkit; const MT = ModelingToolkit
@reexport using ModelingToolkit
using Symbolics
using ModelingToolkit: Symbolic, value, istree, get_states, get_ps, get_iv, get_systems,
                       get_eqs, get_defaults, toparam, get_var_to_name, get_observed, getvar

import ModelingToolkit: get_variables, namespace_expr, namespace_equation, get_variables!,
                        modified_states!, validate, namespace_variables, namespace_parameters,
                        rename, renamespace, getname, flatten
# internal but needed ModelingToolkit functions
import ModelingToolkit: check_variables, check_parameters, _iszero, _merge, check_units, get_unit

import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show
import MacroTools, Graphs
import DataStructures: OrderedDict, OrderedSet
import Parameters: @with_kw_noshow

# globals for the modulate
const DEFAULT_IV = (@parameters t)[1]

# as used in Catlab
const USE_GV_JLL = Ref(false)
function __init__()
    @require Graphviz_jll="3c863552-8265-54e4-a6dc-903eb78fde85" begin
      USE_GV_JLL[] = true
      let cfg = joinpath(Graphviz_jll.artifact_dir, "lib", "graphviz", "config6")
        if !isfile(cfg)
          Graphviz_jll.dot(path -> run(`$path -c`))
        end
      end
    end
  end

# base system type and features
include("reactionsystem.jl")
export Reaction, ReactionSystem, ismassaction, oderatelaw, jumpratelaw
export ODEProblem, SDEProblem, JumpProblem, NonlinearProblem, DiscreteProblem, SteadyStateProblem
export get_constraints, has_constraints, get_combinatoric_ratelaws

# reaction_network macro
const ExprValues = Union{Expr,Symbol,Float64,Int}
include("expression_utils.jl")
include("reaction_network.jl")
export @reaction_network, @add_reactions, @reaction

# registers CRN specific functions using Symbolics.jl
include("registered_functions.jl")
export mm, mmr, hill, hillr, hillar

# functions to query network properties
include("networkapi.jl")
export species, reactionparams, reactions, speciesmap, paramsmap, reactionparamsmap
export numspecies, numreactions, numreactionparams, setdefaults!, symmap_to_varmap
export make_empty_network, addspecies!, addparam!, addreaction!
export dependants, dependents, substoichmat, prodstoichmat, netstoichmat
export conservationlaws, conservedquantities, conservedequations, conservationlaw_constants

# depreciated functions to remove in future releases
export params, numparams

# network analysis functions
export reactioncomplexmap, reactioncomplexes, incidencemat, reactionrates, complexstoichmat
export complexoutgoingmat, incidencematgraph, linkageclasses, deficiency, subnetworks
export linkagedeficiencies, isreversible, isweaklyreversible

# for Latex printing of ReactionSystems
include("latexify_recipes.jl")

# for making and saving graphs
include("graphs.jl")
export Graph, savegraph, complexgraph

end # module


text:
function level_bools(level, agg, term)
    bools = [false for _ in 1:length(agg)]'
    for i in 1:level - 1
        bools = bools .| ((agg .== i) .& (term .== 1))
    end
    bools = bools .| (agg .== level)
    return bools
end

max_stars_repo_name: jungangl/forecast-machine-learning

max_issues_repo_name: jungangl/forecast-machine-learning

max_forks_repo_name: jungangl/forecast-machine-learning


text:
# Linear algebra functions, in particular extending the `LinearAlgebra` stdlib

for f in (:adjoint, :eigvecs, :inv, :pinv, :transpose)
    @eval LinearAlgebra.$f(B::BlockDiagonal) = BlockDiagonal(map($f, blocks(B)))
end

LinearAlgebra.diag(B::BlockDiagonal) = mapreduce(diag, vcat, blocks(B))
LinearAlgebra.det(B::BlockDiagonal) = prod(det, blocks(B))
LinearAlgebra.logdet(B::BlockDiagonal) = sum(logdet, blocks(B))
LinearAlgebra.tr(B::BlockDiagonal) = sum(tr, blocks(B))

for f in (:Symmetric, :Hermitian)
    @eval LinearAlgebra.$f(B::BlockDiagonal, uplo::Symbol=:U) = BlockDiagonal([$f(block, uplo) for block in blocks(B)])
end

# Real matrices can have Complex eigenvalues; `eigvals` is not type stable.
@static if VERSION < v"1.2.0-DEV.275"
    # No convention for sorting complex eigenvalues in earlier versions of Julia.
    function LinearAlgebra.eigvals(B::BlockDiagonal, args...; kwargs...)
        vals = mapreduce(b -> eigvals(b, args...; kwargs...), vcat, blocks(B))
        return !isa(vals, Vector{<:Complex}) ? sort(vals) : vals
    end
else
    # Sorting was introduced in Julia v1.2 by https://github.com/JuliaLang/julia/pull/21598
    function LinearAlgebra.eigvals(
        B::BlockDiagonal, args...; sortby::Union{Function, Nothing}=LinearAlgebra.eigsortby, kwargs...
    )
        vals = mapreduce(b -> eigvals(b, args...; kwargs...), vcat, blocks(B))
        return LinearAlgebra.sorteig!(vals, sortby)
    end
end

if VERSION < v"1.3.0-DEV.426"
    # This is copy of the definition for LinearAlgebra, only used to workaround
    # https://github.com/JuliaLang/julia/issues/31843 which was fixed in Julia v1.3
    function LinearAlgebra.eigmax(B::BlockDiagonal; kwargs...)
        v = eigvals(B; kwargs...)
        if eltype(v) <: Complex
            throw(DomainError(A, "`A` cannot have complex eigenvalues."))
        end
        return maximum(v)
    end
end

"""
    eigen_blockwise(B::BlockDiagonal, args...; kwargs...) -> values, vectors

Computes the eigendecomposition for each block separately and keeps the block diagonal
structure in the matrix of eigenvectors. Hence any parameters given are applied to each
eigendecomposition separately, but there is e.g. no global sorting of eigenvalues.
"""
function eigen_blockwise(B::BlockDiagonal, args...; kwargs...)
    eigens = [eigen(b, args...; kwargs...) for b in blocks(B)]
    # promote to common types to avoid vectors of unclear numerical type
    # This may happen because eigen is not type stable!
    # e.g typeof(eigen(randn(2,2))) may yield Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}
    # or Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}}
    values = promote([e.values for e in eigens]...)
    vectors = promote([e.vectors for e in eigens]...)
    vcat(values...), BlockDiagonal([vectors...])
end

## This function never keeps the block diagonal structure
function LinearAlgebra.eigen(B::BlockDiagonal, args...; kwargs...)
    values, vectors = eigen_blockwise(B, args...; kwargs...)
    vectors = Matrix(vectors) # always convert to avoid type instability (also it speeds up the permutation step)
    @static if VERSION > v"1.2.0-DEV.275"
        Eigen(LinearAlgebra.sorteig!(values, vectors, kwargs...)...)
    else
        Eigen(values, vectors)
    end
end

svdvals_blockwise(B::BlockDiagonal) = mapreduce(svdvals, vcat, blocks(B))
LinearAlgebra.svdvals(B::BlockDiagonal) = sort!(svdvals_blockwise(B); rev=true)

# `B = U * Diagonal(S) * Vt` with `U` and `Vt` `BlockDiagonal` (`S` only sorted block-wise).
function svd_blockwise(B::BlockDiagonal{T}; full::Bool=false) where T
    U = Matrix{float(T)}[]
    S = Vector{float(T)}()
    Vt = Matrix{float(T)}[]
    for b in blocks(B)
        F = svd(b, full=full)
        push!(U, F.U)
        append!(S, F.S)
        push!(Vt, F.Vt)
    end
    return BlockDiagonal(U), S, BlockDiagonal(Vt)
end

function LinearAlgebra.svd(B::BlockDiagonal; full::Bool=false)::SVD
    U, S, Vt = svd_blockwise(B, full=full)
    # Sort singular values in descending order by convention.
    # This means `U` and `Vt` will be `Matrix`s, not `BlockDiagonal`s.
    p = sortperm(S, rev=true)
    return SVD(U[:, p], S[p], Vt[p, :])
end

function LinearAlgebra.cholesky(B::BlockDiagonal)::Cholesky
    C = BlockDiagonal(map(b -> cholesky(b).U, blocks(B)))
    return Cholesky(C, 'U', 0)
end

# Make getproperty on a Cholesky factorized BlockDiagonal return another BlockDiagonal
# where each block is an upper or lower triangular matrix. This ensures that optimizations
# for BlockDiagonal matrices are preserved, though it comes at the cost of reallocating
# a vector of triangular wrappers on each call.
function Base.getproperty(C::Cholesky{T, <:BlockDiagonal{T}}, x::Symbol) where T
    B = getfield(C, :factors)
    uplo = getfield(C, :uplo)
    f = if x === :U
        uplo == 'U' ? UpperTriangular : (X -> UpperTriangular(X'))
    elseif x === :L
        uplo == 'L' ? LowerTriangular : (X -> LowerTriangular(X'))
    elseif x === :UL
        uplo == 'U' ? UpperTriangular : LowerTriangular
    else
        return getfield(C, x)
    end
    return BlockDiagonal(map(f, blocks(B)))
end

# 3-Argument mul!
LinearAlgebra.mul!(C::BlockDiagonal, A::BlockDiagonal, B::BlockDiagonal) = _mul!(C, A, B)

if VERSION ≥ v"1.3"
    function LinearAlgebra.mul!(C::BlockDiagonal, A::BlockDiagonal, B::BlockDiagonal, α::Number, β::Number)
        return _mul!(C, A, B, α, β)
    end
end

function _mul!(C::BlockDiagonal, A::BlockDiagonal, B::BlockDiagonal)
    isequal_blocksizes(A, B) || throw(DimensionMismatch("A and B have different block sizes"))
    isequal_blocksizes(C, A) || throw(DimensionMismatch("C has incompatible block sizes"))
    for i in eachindex(blocks(C))
        @inbounds LinearAlgebra.mul!(C.blocks[i], A.blocks[i], B.blocks[i])
    end

    return C
end

function _mul!(C::BlockDiagonal, A::BlockDiagonal, B::BlockDiagonal, α::Number, β::Number)
    isequal_blocksizes(A, B) || throw(DimensionMismatch("A and B have different block sizes"))
    isequal_blocksizes(C, A) || throw(DimensionMismatch("C has incompatible block sizes"))
    for i in eachindex(blocks(C))
        @inbounds LinearAlgebra.mul!(C.blocks[i], A.blocks[i], B.blocks[i], α, β)
    end

    return C
end

function LinearAlgebra.:\(B::BlockDiagonal, vm::AbstractVecOrMat)
    row_i = 1
    # BlockDiagonals with non-square blocks
    if !all(is_square, blocks(B))
        return Matrix(B) \ vm # Fallback on the generic LinearAlgebra method
    end
    result = similar(vm)
    for block in blocks(B)
        nrow = size(block, 1)
        result[row_i:(row_i + nrow - 1), :] = block \ vm[row_i:(row_i + nrow - 1), :]
        row_i += nrow
    end
    return result
end


text:
using IncrementalInference, KernelDensityEstimate

fg = emptyFactorGraph()

N=100

doors = reshape(Float64[-100.0;0.0;100.0;300.0],1,4)
pd = kde!(doors,[3.0])
pd = resample(pd,N);
bws = getBW(pd)[:,1]
doors2 = getPoints(pd);

v1 = addNode!(fg,:x0,doors,N=N)
f1  = addFactor!(fg,[v1],Obsv2( doors2, reshape(bws,1,1), [1.0])) #, samplefnc=getSample

tem = 2.0*randn(1,N)+getVal(v1)+50.0
v2 = addNode!(fg,:x2, tem, N=N)
addFactor!(fg, [:x0; :x2], Odo(50.0*ones(1,1),2.0*ones(1,1),[1.0])) #, samplefnc=getSample
# addFactor!(fg, [v1;v2], Odo(50.0*ones(1,1),[2.0]',[1.0])) #, samplefnc=getSample

# monocular sighting would look something like
#addFactor!(fg, Mono, [:x3,:l1], [14.0], [1.0], [1.0])
#addFactor!(fg, Mono, [:x4,:l1], [11.0], [1.0], [1.0])

v3=addNode!(fg,:x3,4.0*randn(1,N)+getVal(v2)+50.0, N=N)
addFactor!(fg,[v2;v3],Odo(50.0*ones(1,1),4.0*ones(1,1),[1.0])) #, samplefnc=getSample
f2 = addFactor!(fg,[v3], Obsv2(doors2, bws', [1.0]))

v4=addNode!(fg,:x4,2.0*randn(1,N)+getVal(v3)+50.0, N=N)
addFactor!(fg,[v3;v4],Odo(50.0*ones(1,1),2.0*ones(1,1),[1.0])) #, samplefnc=getSample

    l1=addNode!(fg, :l1, 0.5*randn(1,N)+getVal(v3)+64.0, N=N)
    addFactor!(fg, [v3,l1], Ranged([64.0],[0.5],[1.0])) #, samplefnc=getSample
    addFactor!(fg, [v4,l1], Ranged([16.0],[0.5],[1.0])) #, samplefnc=getSample

v5=addNode!(fg,:x5,2.0*randn(1,N)+getVal(v4)+50.0, N=N)
addFactor!(fg,[v4;v5],Odo(50.0*ones(1,1),2.0*ones(1,1),[1.0])) #, samplefnc=getSample

v6=addNode!(fg,:x6,1.25*randn(1,N)+getVal(v5)+40.0, N=N)
addFactor!(fg,[v5;v6],Odo(40.0*ones(1,1),1.20*ones(1,1),[1.0])) #, samplefnc=getSample

v7=addNode!(fg,:x7,2.0*randn(1,N)+getVal(v6) +60.0, N=N)
addFactor!(fg,[v6;v7],Odo(60.0*ones(1,1),2.0*ones(1,1),[1.0])) #, samplefnc=getSample

f3 = addFactor!(fg,[v7], Obsv2(doors, reshape(bws,1,1), [1.0])) #, samplefnc=getSample

# HMM computed ground truth, extended for 7 poses with landmark
gt = Dict{Symbol, Array{Float64,2}}()
gt[:x0]=reshape(Float64[0.0;1.97304 ],2,1) # -0.0342366
gt[:x2]=reshape(Float64[50.0; 2.83153 ],2,1) # 49.8797
gt[:x3]=reshape(Float64[100.0; 1.65557 ],2,1) # 99.8351
gt[:x4]=reshape(Float64[150.0; 1.64945 ],2,1) # 148.637
gt[:x5]=reshape(Float64[200.0; 1.77992 ],2,1) # 198.62
gt[:x6]=reshape(Float64[240.0; 2.20466 ],2,1) # 238.492
gt[:x7]=reshape(Float64[300.0; 2.14353 ],2,1) # 298.467
gt[:l1]=reshape(Float64[165.0; 1.17284 ],2,1) # 164.102

tree = prepBatchTree!(fg, drawpdf=false);

# list vertices in fg
@show xx,ll = ls(fg)

# do belief propagation inference over tree once
# using recursive single core approach (better stack trace for development)
# inferOverTreeR!(fg, tree)
inferOverTreeR!(fg, tree, N=N, dbg=true)
#
# test multi-processor solve (operational fast solving)
inferOverTree!(fg, tree)
# inferOverTree!(fg, tree, dbg=true)

#


text:
using SimpleDiffEq, StaticArrays, OrdinaryDiffEq, Test

function loop(u, p, t)
    @inbounds begin
        σ = p[1]; ρ = p[2]; β = p[3]
        du1 = σ*(u[2]-u[1])
        du2 = u[1]*(ρ-u[3]) - u[2]
        du3 = u[1]*u[2] - β*u[3]
        return SVector{3}(du1, du2, du3)
    end
end
function liip(du, u, p, t)
    σ = p[1]; ρ = p[2]; β = p[3]
    du[1] = σ*(u[2]-u[1])
    du[2] = u[1]*(ρ-u[3]) - u[2]
    du[3] = u[1]*u[2] - β*u[3]
    return nothing
end

u0 = 10ones(3)
dt = 1e-2

odeoop = ODEProblem{false}(loop, SVector{3}(u0), (0.0, 100.0),  [10, 28, 8/3])
sol  = solve(odeoop,SimpleATsit5()   ,dt=dt)
sol2 = solve(odeoop,GPUSimpleATsit5(),dt=dt,abstol=1e-6,reltol=1e-3)

@test sol.u == sol2.u
@test sol.t == sol2.t

sol  = solve(odeoop,Tsit5()          ,dt=dt,saveat=0.0:0.1:100.0)
sol2 = solve(odeoop,GPUSimpleATsit5(),dt=dt,saveat=0.0:0.1:100.0,abstol=1e-6,reltol=1e-3)
sol3 = solve(odeoop,SimpleATsit5()   ,dt=dt,saveat=0.0:0.1:100.0)

@test sol[20] ≈ sol2[20]
@test sol2.u == sol3.u
@test sol.t  == sol2.t

dt = 1e-1

sol  = solve(odeoop,SimpleTsit5()   ,dt=dt)
sol2 = solve(odeoop,GPUSimpleTsit5(),dt=dt)

@test sol.u == sol2.u
@test sol.t == sol2.t


text:
using Images

function lerp(v0, v1, t)
  v0 + t * (v1 - v0)
end

function make_color_scale(
    color1::Tuple{Int, Int, Int},
    color2::Tuple{Int, Int, Int},
    steps::Int,
    step_size::Int,
    height::Int)

  img = Array{RGB}(undef, height, step_size * steps)
  println("generated empty image of size ", size(img))

  for i = 1:steps
    println("step ", i)
    println("calculating color...")
    this_color = RGB(
      lerp(color1[1], color2[1], 1.0/convert(AbstractFloat, step_size)),
      lerp(color1[2], color2[2], 1.0/convert(AbstractFloat, step_size)),
      lerp(color1[3], color2[3], 1.0/convert(AbstractFloat, step_size))
    )

    println("assigning color to section of image")
    # this line causes problems. Perhaps I don't understand array indexing in
    # julia
    img[1:height][(i-1)*step_size+1:i*step_size] .= this_color

    println()
  end
end


text:
"""
ProjectiveCamera can be 'regular' or orthographic depending on the projection.
"""
struct ProjectiveCamera{F<:Film} <: Camera
    film::F
    camera_to_world::Transformation
    raster_to_camera::Transformation
    camera_to_raster::Transformation
    camera_to_screen::Transformation
    screen_to_camera::Transformation
    dxcam::Vector3
    dycam::Vector3
    lensr::Float64
    focald::Float64
end

"""
    ProjectiveCamera(cam2world, projection, window, lensR, focal_depth, film)

Construct a ProjectiveCamera.
"""
function ProjectiveCamera(cam2w::Transformation, cam2screen::Transformation, window::Array{Float64,1}, lensr::Float64, focald::Float64, f::Film)
    screen2raster = (scaling(f.resX, f.resY, 1.0)
                   * scaling(1.0/(window[2] - window[1]), 1.0/(window[4] - window[3]), 1.0)
                   * translation(-window[1], -window[3], 0.0))
    raster2cam = inv(cam2screen) * inv(screen2raster)
    dxcam = raster2cam(Vector3(1,0,0))
    dycam = raster2cam(Vector3(0,1,0))
    ProjectiveCamera(f, cam2w, raster2cam, inv(raster2cam), cam2screen, inv(cam2screen), dxcam, dycam, lensr, focald)
end

"""
    OrthographicCamera(cam2world, window, lensR, focal_depth, film)

Construct an orthographic camera.
"""
function OrthographicCamera(cam2w::Transformation, window::Array{Float64,1}, lensr::Float64, focald::Float64, f::Film)
    ProjectiveCamera(cam2w, orthographic(0.0, 1.0), window, lensr, focald, f)
end

"""
    OrthographicCamera(cam2world, Xwidth, Ywidth, lensR, focal_depth, film)

Construct an orthographic camera.
"""
function OrthographicCamera(cam2w::Transformation, Xwidth::Real, Ywidth::Real, lensr::Float64, focald::Float64, f::Film)
    ProjectiveCamera(cam2w, orthographic(0.0, 1.0), [-Xwidth/2, Xwidth/2, -Ywidth/2, Ywidth/2], lensr, focald, f)
end

"""
    orthographic(znear::Real, zfar::Real)

Computes the orthographic projection transformation given near and far distances.
"""
function orthographic(znear::Real, zfar::Real)
    scaling(1.0, 1.0, 1.0/(zfar - znear)) * translation(0.0, 0.0, -znear)
end

"""
    generate_ray(C::ProjectiveCamera, sample::CameraSample)

Generate a ray leaving the camera based on the sample.
"""
function generate_ray(C::ProjectiveCamera, sample::CameraSample)
    pras = Point3(sample.imgX, sample.imgY, 0.0)
    pcam = C.raster_to_camera(pras)
    R = Ray(pcam, Vector3(0,0,1), 0.0, Inf, 1)
    return 1.0, C.camera_to_world(R)
end

"""
    count_tasks(camera::ProjectiveCamera, nCores::Integer)

A heuristic function to count the number of rendering tasks given the number of processor
cores and the resolution of the output image. This number is either four tasks per core,
or the number of 16x16 pixel blocks in the output, whichever is higher..

"""
function count_tasks(camera::ProjectiveCamera, nCores::Integer)
    npix = camera.film.resX * camera.film.resY
    nTasks = max(4*nCores, div(npix, 256))
    return round_pow2(nTasks)
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
module PolicyViz

using GridInterpolations, Interact, PGFPlots, Colors, ColorBrewer

include("./viz_policy_constants.jl")
include("./viz_policy_drl.jl")
include("./viz_policy.jl")
include("./nnet_calculations.jl")
include("./policy_functions.jl")

end #module PolicyViz


text:
# ------------------------------------------------------------------------------
# Tests of the FP Growth algorithm and the associated functions
#
# TODO Test tree building functions
# TODO Test
# ------------------------------------------------------------------------------

max_stars_repo_name: DopedCow/FPGrowth

max_issues_repo_name: DopedCow/FPGrowth

max_forks_repo_name: DopedCow/FPGrowth


text:
function numerical_integration(x, y)
    sum(diff(x) .* midpoints(y))
end
function cdf_distance(F_X, F_Y, x, p)
    if p == Inf
        return maximum(abs.(F_X .- F_Y))
    end
    diff_p = abs.(F_X .- F_Y).^p
    return numerical_integration(x, diff_p)^(1/p)
end
function cdf_distance(Xdistr::UnivariateDistribution, Ydistr::UnivariateDistribution, x, p)
    F_X = cdf.(Xdistr, x)
    F_Y = cdf.(Ydistr, x)
    return cdf_distance(F_X, F_Y, x, p)
end

function distance_test_statistic(p, method::DeconvolutionMethod, xgrid::AbstractVector)
    function ϕ(X, Y, σ_X, σ_Y)
        ϵ_X = Normal.(0.0, σ_X)
        Xdistr_hat = decon(method, X, ϵ_X)

        ϵ_Y = Normal.(0.0, σ_Y)
        Ydistr_hat = decon(method, Y, ϵ_Y)
        return cdf_distance(Xdistr_hat, Ydistr_hat, xgrid, p)
    end
    return ϕ
end

KS_test_statistic(X, Y, σ_X, σ_Y) = HypothesisTests.ApproximateTwoSampleKSTest(X, Y).δ


text:
module HyperGraphs

using SimpleTraits.SimpleTraits
using LinearAlgebra: Diagonal

include("core.jl")
include("conversions.jl")
include("matrices.jl")
include("operations.jl")
include("properties.jl")
include("queries.jl")
include("utils.jl")

export

# core: types
AbstractHyperGraph, AbstractHyperEdge,
HyperEdge, HyperGraph,
SpeciesSet, ChemicalHyperEdge, ChemicalHyperGraph,

# core: traits
OrientedStructs, WeightedStructs, # constants holding type unions
IsOriented, IsWeighted, # traits
isoriented, isweighted, # trait functions

# core: functions
objects, multiplicities, # incidence functions
vertices, hyperedges, # hyperedges and hypergraphs functions
src, tgt, src_multiplicities, tgt_multiplicities, # oriented hyperedges functions
weight, weights, # weighted hyperedges functions

# core: type functions
species, stoich, src_stoich, tgt_stoich, inputs_stoich, outputs_stoich,
inputs, outputs, rate,

# matrices
degree_matrix, incidence_matrix, cardinality_matrix, weight_matrix,

# operations
add_vertex!, add_vertices!, del_vertex!, del_vertices!,
add_hyperedge!, add_hyperedges!, del_hyperedge!, del_hyperedges!, del_empty_hyperedges!,
replace_vertex!, replace_hyperedge!, replace_hyperedges!,
update_weight!, set_all_weights!,
merge,
subhypergraph,
switch,

# properties
degree, indegree, outdegree, degrees, indegrees, outdegrees, # vertices properties
length, cardinality, cardinalities, nsrcs, ntgts, # hyperedges properties
nv, nhe, delta, Δ, rank, co_rank, order, hypergraph_size, volume, # hypergraphs properties

# queries: Boolean
has_vertex, has_vertices, has_hyperedge, has_hyperedges,
in_src, in_tgt,
isincident,
issimple, isgraph,
issubhyperedge, isuniquesubhyperedge,
isloop, is_positive_loop, has_loops,
iskuniform, iskregular,
is_netstoich_null, iscatalyst,
has_parallel_hyperedges, has_multi_hyperedges,

# queries: content
incident_hyperedges,
empty_hyperedges, num_empty_hyperedges, has_empty_hyperedges,
neighbors, inneighbors, outneighbors, all_neighbors,
loops, positive_loops, num_loops,
catalysts,
parallel_hyperedges, num_parallel_hyperedges,
multi_hyperedges, num_multi_hyperedges

end # module


text:
using BHAPtfem, Test

ProjDir = dirname(@__FILE__)
ProjName = split(ProjDir, "/")[end]

bha[:segs] = [
# Element type,  Material,    Length,     OD,         ID,      fc,    NoOfElements
  (:bit,          :steel,     0.00,   2.75,   9.0,  0.0),
  (:collar,       :steel,    45.00,   2.75,   7.0,  0.0),
  (:stabilizer,   :steel,     0.00,   2.75,   9.0,  0.0),
  (:pipe,         :steel,    30.00,   2.75,   7.0,  0.0),
  (:stabilizer,   :steel,     0.00,   2.75,   9.0,  0.0),
  (:pipe,         :steel,    90.00,   2.75,   7.0,  0.0),
  (:stabilizer,   :steel,     0.00,   2.75,   9.0,  0.0),
  (:pipe,         :steel,   100.00,   2.75,   7.0,  0.0)
]

traj = [
#   Heading,      Diameter
  ( 60.0,      9.0)
]

wobs = 20:10:40
incls = 20:10:40               # Or e.g. incls = [5 10 20 30 40 45 50]

properties, nodes, elements = input!(bha)

bha[:segments] |> display
println("\n")

bha[:properties] |> display
println("\n")

bha[:nodes] |> display
println("\n")

bha[:elements] |> display
println("\n")

pfem = problem(bha, traj, wobs, incls, pdir=ProjDir)

pfem1 = Dict(
  # Frame(nels, nn, ndim, nst, nip, finite_element(nod, nodof))
  :struc_el => Frame(10, 11, 2, 1, 1, Line(2, 3)),
  :properties => [
    1.07207e9 4.53116e9;],
  :x_coords =>0.0:1.0:10.0,
  :y_coords => [0.0 for i in 1:11],
  :g_num => [
    1 2 3 4 5 6 7 8 9 10;
    2 3 4 5 6 7 8 9 10 11],
  :support => [
    (1, [1, 0, 1])
    (11, [0, 0, 1])
    ],
  :loaded_nodes => [
    (1, [10000.0 0.0 0.0]),
    (5, [0.0 10000.0 0.0]),
    ],
  :penalty => 1e19,
  #=
  :eq_nodal_forces_and_moments => [
    (1, [0.0 -60.0 -60.0 0.0 -60.0 60.0]),
    (2, [0.0 -120.0 -140.0 0.0 -120.0 140.0]),
    (3, [0.0 -20.0 -6.67 0.0 -20.0 6.67])
  ]
  =#
)

fem, dis_df, fm_df = solve(pfem1)

dis_df


text:
using AtomicLevels

function nonh()
    info("Executing nonh")
    # Angular integrals
    clean_mpi_tmp(r"fort.[0-9]+",
                  r"yint.lst.[0-9]+",
                  r"c.lst.[0-9]+")
    mpi_run(mpi_cmd("$atsp/nonh_mpi"))
    info("Finished executing nonh")
end

function hf(name, term::Term, Z,
            closed::Config, ref_set::Config,
            vary = "all",
            default_electron_parameters = true,
            rem_parameters = ["y","n","T,1.D-8","y","y"],
            additional_parameters = [])
    closed = join([@sprintf("%2d%s", c[1], ells[c[2]+1]) for c in closed], " ")

    pipe_file_run("$atsp/hf",
                  """$name,$(string(term,false)),$(float(Z)) ! Name, final term, Z (next row: closed orbitals)
 $closed
$(orbital_string(ref_set, true, false; escape_principal_qn = true)) ! Electrons outside closed orbitals
$vary ! Which orbitals to vary
$(y_or_n(default_electron_parameters)) ! Default electron parameters
$(y_or_n(length(rem_parameters)==0)) ! Remaining parameters
$(join(rem_parameters, "\n"))
$(y_or_n(length(additional_parameters)>0)) ! Additional parameters
$(join(additional_parameters, "\n"))
n ! Don't continue along the sequence
""")
end

function mchf(name,Z,
              term_weights,
              active::Config,
              spectroscopic::Config = [];
              iter = 200,
              cfg_tol = 1e-8,
              scf_tol = 1e-8)
    pipe_file_run("$atsp/mchf_mpi",
                  """$name,$(float(Z))
$(join(term_weights, '\n'))
$(orbital_string(active,false,false,","; escape_principal_qn = true))
$(orbital_string(spectroscopic,false,false,"," ; escape_principal_qn = true))
y ! Default electron parameters
n ! Custom other parameters
y ! Default values for NO,REL,STRONG
n ! Custom values for PRINT,CFGTOL,SCFTOL
F,$cfg_tol,$scf_tol
n ! Custom values for NSCF,IC
$iter,0
y ! Default values for ACFG,LD,TRACE
""")
end

function atsp_save_run(config, term)
    cpf("wfn.out", "wfn.inp")
    name = active_file(config, term)
    cpf("wfn.out", "$name.w")
    cpf("summry", "$name.s")
    cpf("clist.out", "$name.c")
    cpf(filter(f -> ismatch(Regex("^$(string(term)).*\\.l\$"), f), readdir())[1], "$name.l")
end

function breit_pauli(config, term;
                     restarting = false,
                     guesses = false,
                     which_op = 2,
                     all_rel = true,
                     all_inter = true,
                     atom_mass = Inf)
    name = active_file(config, term)

    (isfile("$name.l") && isfile("$name.c")) || error("Could not find $name.l/c")
    cpf("$name.l", "$name.l.bak")
    Js = J_range(term)
    j2 = map(Js) do JJ
        round(Int, 2JJ)
    end
    eiv = join([1 for JJ in reverse(Js)], '\n')
    println("$term ⟹ J ∈ $Js")

    atom_mass_s = (atom_mass == Inf) ? "" : "\n$atom_mass"

    pipe_file_run("$atsp/bpci",
                  """$name, y, n ! Atom, relativistic, mass correction
$(y_or_n(restarting)) ! Restarting
$(y_or_n(guesses)) ! Existing $name.l/.j as guesses
$(maximum(j2)),$(minimum(j2)) ! max(2J),min(2J)
$eiv ! Eigenvalues
$which_op ! Non-relativistic operators and selected relativistic
$(y_or_n(all_rel)) ! All relativistic operators
$(y_or_n(all_inter)) ! All interactions
$(y_or_n(atom_mass == Inf)) ! Default Rydberg constant$(atom_mass_s)
""")
    cpf("$name.l.bak", "$name.l")
end

function atsp_cp_wfn(a,b)
    mkpath("$b")
    cpf("$a/wfn.out", "$b/wfn.inp")
end

function read_hf_eng()
    # Read the total energy from hf.log
    open("hf.log") do file
        float(split(readall(file))[end-11])
    end
end

function read_eng_blocks(f::Function, config, term, ending)
    open("$(active_file(config, term)).$ending") do file
        f() do
            while !ismatch(r"^[ ]*2\*J =", readline(file)); end
            readline(file)
            float(split(readline(file))[2])
        end
    end
end

function read_mchf_eng(config, term)
    # Read the total (non-relativistic) energy from $name.l
    read_eng_blocks(config, term, "l") do read_block
        read_block()
    end
end

function read_breit_pauli_eng(config, term)
    Js = J_range(term)
    # Read the total (Breit–Pauli corrected) energies from $name.j
    read_eng_blocks(config, term, "j") do read_block
        [read_block()
         for JJ in reverse(Js)]
    end
end

function atsp_clean()
    # Clean up temporary files that take up a lot of space
    for f in filter(f -> contains(f, ".lst"), readdir())
        rm(f)
    end
end

function hf_mchf_bp(config::Config,
                    term::Term,
                    Z,
                    ncorr,
                    wfn0 = nothing;
                    active::Function = active_set,
                    overwrite = true,
                    atom_mass = Inf,
                    csf_filter = (a...) -> true,
                    hf_core = true,
                    bp = true)
    println(repeat("=", 80))
    conf = active_file(config, term)
    println(config, " ", term)

    overwrite_i = i -> overwrite && isfile("$i/wfn.out") || !isfile("$i/wfn.out")
    dir_run(conf) do
        energies = zeros(ncorr+1)

        Js = J_range(term)
        bp_energies = zeros(ncorr+1, length(Js))

        # Hartree--Fock run
        overwrite_i(0) && dir_run("0") do
            if wfn0 != nothing
                println("Initial guess: $wfn0")
                cpf(wfn0, "wfn.inp")
            end
            csfgenerate(active(config, 0),
                        csf_filter,
                        (term, 0, config))
            nonh()
            hf(conf, term, Z,
               closed(config),
               open(config),
               hf_core ? "all" : "=$(length(open(config)))")
            atsp_clean()
        end

        dir_run("0") do
            energies[1] = read_hf_eng()
            bp_energies[1,:] = energies[1]
        end

        # Multiconfigurational Hartree--Fock, layer by layer, with
        # additional Breit--Pauli calculation at the end of each run
        for i = 1:ncorr
            atsp_cp_wfn(i-1,i)
            overwrite_i(i) && dir_run("$i") do
                act = csfgenerate(active(config, i),
                                  csf_filter,
                                  (term, 2, config))
                nonh()
                mchf(conf, Z, 1,
                     act, config)
                atsp_save_run(config, term)
                bp && breit_pauli(config, term; atom_mass = atom_mass)
                atsp_clean()
            end

            dir_run("$i") do
                energies[i+1] = read_mchf_eng(config, term)
                bp && (bp_energies[i+1,:] = read_breit_pauli_eng(config, term))
            end
        end

        save_eng([energies bp_energies]')
    end
    conf, abspath("$conf/$ncorr/wfn.out")
end

export nonh, hf, mchf,
atsp_save_run, breit_pauli,
atsp_cp_wfn, read_hf_eng, read_mchf_eng, read_breit_pauli_eng,
atsp_clean,
hf_mchf_bp


text:
export
    # cudnnConvolutionMode_t
    CUDNN_CONVOLUTION,
    CUDNN_CROSS_CORRELATION,

    # cudnnConvolutionFwdPreference_t
    CUDNN_CONVOLUTION_FWD_NO_WORKSPACE,
    CUDNN_CONVOLUTION_FWD_PREFER_FASTEST,
    CUDNN_CONVOLUTION_FWD_SPECIFY_WORKSPACE_LIMIT

type ConvolutionDesc
    ptr::Ptr{Void}

    function ConvolutionDesc()
        p = Ptr{Void}[0]
        cudnnCreateConvolutionDescriptor(p)
        desc = new(p[1])
        finalizer(desc, cudnnDestroyConvolutionDescriptor)
        desc
    end
end

Base.unsafe_convert(::Type{Ptr{Void}}, desc::ConvolutionDesc) = desc.ptr

function convolution{T,N}(x::CuArray{T}, w::CuArray{T}, b::CuArray{T},
    pads::NTuple{N,Int}, strides::NTuple{N,Int}; mode=CUDNN_CROSS_CORRELATION)

    convdesc = ConvolutionDesc()
    c_pads = Cint[pads[i] for i=N:-1:1]
    c_strides = Cint[strides[i] for i=N:-1:1]
    c_upscale = fill(Cint(1), N)
    cudnnSetConvolutionNdDescriptor(convdesc, N, c_pads, c_strides, c_upscale, mode, datatype(T))

    outdims = ntuple(N) do i
        (size(x,i) + 2pads[i] - size(w,i)) ÷ strides[i] + 1
    end
    y = similar(x, outdims..., size(w,N+2), size(x,N+2))
    xdesc = TensorDesc(x)
    wdesc = FilterDesc(w)
    ydesc = TensorDesc(y)

    h = handle(x)
    p = cudnnConvolutionFwdAlgo_t[0]
    cudnnGetConvolutionForwardAlgorithm(h, xdesc, wdesc, convdesc, ydesc,
        CUDNN_CONVOLUTION_FWD_PREFER_FASTEST, 0, p)
    algo = p[1]

    p = Cint[0]
    cudnnGetConvolutionForwardWorkspaceSize(h, xdesc, wdesc, convdesc, ydesc, algo, p)
    worksize = p[1]
    workspace = CuArray{Int8}(Int(worksize))

    cudnnConvolutionForward(h, T[1], xdesc, x, wdesc, w, convdesc,
        algo, workspace, worksize, T[0], ydesc, y)
    dims1 = ntuple(_ -> 1, ndims(y)-ndims(b))
    b = reshape(b, size(b)..., dims1...)
    addtensor!(b, y)

    function backward!(gy, gx, gw, gb)
        isvoid(gw) || ∇convolution_filter!(xdesc, x, ydesc, gy, convdesc, wdesc, gw)
        isvoid(gx) || ∇convolution_data!(wdesc, w, ydesc, gy, convdesc, xdesc, gx)
        isvoid(gb) || ∇convolution_bias!(ydesc, gy, bdesc, gb)
    end
    y, backward!
end

function ∇convolution_filter!(xdesc, x, dydesc, dy, convdesc, dwdesc, dw)
    h = handle(x)
    p = cudnnConvolutionBwdFilterAlgo_t[0]
    cudnnGetConvolutionBackwardFilterAlgorithm(h, xdesc, dydesc, convdesc, dwdesc,
        CUDNN_CONVOLUTION_BWD_FILTER_PREFER_FASTEST, 0, p)
    algo = p[1]

    p = Cint[0]
    cudnnGetConvolutionBackwardFilterWorkspaceSize(h, xdesc, dydesc, convdesc, dwdesc, algo, p)
    worksize = p[1]
    workspace = CuArray{Int8}(Int(p[1]))

    T = eltype(x)
    cudnnConvolutionBackwardFilter(h, T[1], xdesc, x, dydesc, dy, convdesc,
        algo, workspace, worksize, T[1], dwdesc, dw)
end

function ∇convolution_data!(wdesc, w, dydesc, dy, convdesc, dxdesc, dx)
    h = handle(dy)
    p = cudnnConvolutionBwdDataAlgo_t[0]
    cudnnGetConvolutionBackwardDataAlgorithm(h, wdesc, dydesc, convdesc, dxdesc,
        CUDNN_CONVOLUTION_BWD_DATA_PREFER_FASTEST, 0, p)
    algo = p[1]

    p = Cint[0]
    cudnnGetConvolutionBackwardDataWorkspaceSize(h, wdesc, dydesc, convdesc,
        dxdesc, algo, p)
    worksize = p[1]
    workspace = CuArray{Int8}(Int(worksize))

    T = eltype(dy)
    cudnnConvolutionBackwardData(h, T[1], wdesc, w, dydesc, dy, convdesc,
        algo, workspace, worksize, T[1], dxdesc, dx)
end

function ∇convolution_bias!(dydesc, dy, dbdesc, db)
    h = handle(dy)
    T = eltype(dy)
    cudnnConvolutionBackwardBias(h, T[1], dydesc, dy, T[1], dbdesc, db)
end


text:
import Base.show

function show(io::IO,pdb::PdbTool.Pdb)
	write(io,"Pdb with $(length(pdb.chain)) chain(s)")
end

function show(io::IO,chain::PdbTool.Chain)
	write(io,"Chain with $(length(chain.residue)) residue(s), $(length(chain.align)) mapped")
end

function show(io::IO,res::PdbTool.Residue)
	write(io,"Residue with $(length(res.atom)) atom(s)")
end

max_stars_repo_name: christophfeinauer/PdbTool

max_issues_repo_name: christophfeinauer/PdbTool

max_forks_repo_name: christophfeinauer/PdbTool


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM.Testing

@testset "test 2d linear elasticity with surface + volume load" begin
    meshfile = "/geometry/2d_block/BLOCK_1elem.med"
    mesh = aster_read_mesh(Pkg.dir("JuliaFEM")*meshfile)

    # field problem
    block = Problem(Elasticity, "BLOCK", 2)
    block.properties.formulation = :plane_strain
    block.properties.finite_strain = false
    block.properties.geometric_stiffness = false
    block.elements = create_elements(mesh, "BLOCK")
    update!(block.elements, "youngs modulus", 288.0)
    update!(block.elements, "poissons ratio", 1/3)

    # traction
    traction = Problem(Elasticity, "TRACTION", 2)
    traction.properties.formulation = :plane_strain
    traction.properties.finite_strain = false
    traction.properties.geometric_stiffness = false
    traction.elements = create_elements(mesh, "TOP")
    update!(traction, "displacement traction force 2", 288.0*9/8)

    # boundary conditions
    bc_sym_23 = Problem(Dirichlet, "symmetry bc 23", 2, "displacement")
    bc_sym_23.elements = create_elements(mesh, "LEFT")
    update!(bc_sym_23, "displacement 1", 0.0)
    bc_sym_13 = Problem(Dirichlet, "symmetry bc 13", 2, "displacement")
    bc_sym_13.elements = create_elements(mesh, "BOTTOM")
    update!(bc_sym_13, "displacement 2", 0.0)

    solver = LinearSolver(block, traction, bc_sym_23, bc_sym_13)
    solver()

    info("u = ", block.assembly.u)
    info("λ = ", block.assembly.la)

end


text:
using Test
using BetweenFlags
import BetweenFlags

@testset "Greedy - length(flag)>length(word bc)" begin
  flag_set = FlagSet([
      FlagPair{GreedyType}(
          StartFlag("{", ["*"], [""]),
          StopFlag( "}", [""], ["*"])
      )
  ])

  text = "Foo, bar..."
  token_stream = TokenStream(text, flag_set)
  @test token_stream("{-}") == ""

  text = "Foo,*{}*, bar..."
  token_stream = TokenStream(text, flag_set; rm_outer_bcs=true)
  @test token_stream("{-}") == "{}"

  text = "Foo,*{a}{b}*, bar..."
  token_stream = TokenStream(text, flag_set; rm_outer_bcs=true)
  @test token_stream("{-}") == "{a}{b}"

  text = "Foo, *{a} {b}*, bar..."
  token_stream = TokenStream(text, flag_set; rm_outer_bcs=true)
  # Consecutive scopes are concatenated (in `get_string`):
  @test token_stream("{-}") == "{a} {b}"

  text = "Foo, *{bar}*, foobar..."
  token_stream = TokenStream(text, flag_set; rm_outer_bcs=true)
  @test token_stream("{-}") == "{bar}"

  text = "Foo, *{{foobar}*}, foobaz..."
  token_stream = TokenStream(text, flag_set; rm_outer_bcs=true)
  @test token_stream("{-}") == "{{foobar}"
end

@testset "Scope - flag subsets" begin
  flag_set = FlagSet([
      FlagPair{ScopeType}(
          StartFlag("do", ["*"], [""]),
          StopFlag( "end do", [""], ["*"])
      )
  ])
  text = "Foo, *do bar; foo end do*, foobar..."
  token_stream = TokenStream(text, flag_set; rm_outer_bcs=true)
  @test token_stream("do-end do") == "do bar; foo end do"
end

@testset "Scope - nested" begin
  flag_set = FlagSet([
      FlagPair{ScopeType}(
          StartFlag("{", ["*"], [""]),
          StopFlag( "}", [""], ["*"])
      )
  ])

  text = "Foo, *{bar}*, foobar..."
  token_stream = TokenStream(text, flag_set; rm_outer_bcs=true)
  @test token_stream("{-}") == "{bar}"

  text = "Foo, *{{bar}}*, foobar..."
  token_stream = TokenStream(text, flag_set; rm_outer_bcs=true)
  @test token_stream("{-}") == "{{bar}}"

  text = "Foo, *{{bar}{baz}}*, foobar..."
  token_stream = TokenStream(text, flag_set; rm_outer_bcs=true)
  @test token_stream("{-}") == "{{bar}{baz}}"

end

@testset "Greedy - length(flag)==length(word bc)" begin
  flag_set = FlagSet([
      FlagPair{GreedyType}(
          StartFlag("ABC_ABC_L", ["STA_WBC_L"], ["STA_WBC_R"]),
          StopFlag( "ABC_ABC_R", ["STO_WBC_L"], ["STO_WBC_R"])
      )
  ])

  text = "foo STA_WBC_LABC_ABC_LSTA_WBC_R bar STO_WBC_LABC_ABC_RSTO_WBC_R foobar"
  token_stream = TokenStream(text, flag_set; rm_outer_bcs=true)
  @test token_stream("ABC_ABC_L-ABC_ABC_R") == "ABC_ABC_LSTA_WBC_R bar STO_WBC_LABC_ABC_R"

end

@testset "Greedy - length(flag)<length(word bc)" begin
  flag_set = FlagSet([
      FlagPair{GreedyType}(
          StartFlag("ABC_L", ["XYZ_XYZ_STA_WBC_L"], ["XYZ_XYZ_STA_WBC_R"]),
          StopFlag( "ABC_R", ["XYZ_XYZ_STO_WBC_L"], ["XYZ_XYZ_STO_WBC_R"])
      )
  ])

  text = "foo XYZ_XYZ_STA_WBC_LABC_LXYZ_XYZ_STA_WBC_R bar XYZ_XYZ_STO_WBC_LABC_RXYZ_XYZ_STO_WBC_R foobar"
  token_stream = TokenStream(text, flag_set; rm_outer_bcs=true)
  @test token_stream("ABC_L-ABC_R") == "ABC_LXYZ_XYZ_STA_WBC_R bar XYZ_XYZ_STO_WBC_LABC_R"

end

@testset "start_flag == stop_flag, same word bcs" begin
  flag_set = FlagSet([
        FlagPair{GreedyType}(
        StartFlag("|", ["*"], [" "]),
        StopFlag( "|", [" "], ["*"])
        )
        ])
  text = "baz*| foo |*baz*| bar |*foobaz"
  token_stream = TokenStream(text, flag_set; rm_outer_bcs=true)
  @test token_stream("|-|") == "| foo || bar |"
end

@testset "Unicode" begin
  flag_set = FlagSet([
        FlagPair{GreedyType}(
        StartFlag("|", ["*"], [" "]),
        StopFlag( "|", [" "], ["*"])
        )
        ])
  text = "baz*| θ |*baz*| ϕ |*foobaz"
  token_stream = TokenStream(text, flag_set; rm_outer_bcs=true)
  @test token_stream("|-|") == "| θ || ϕ |"
end


text:
using CPUTime
using Plots
using StaticArrays
using ArraysOfArrays

function split_vec_of_arrays(u)
    vec.(u) |>
    x -> VectorOfSimilarVectors(x).data |>
    transpose |>
    VectorOfSimilarVectors
end

function jisaku_solve_euler(dx, x₀, t_span, Δt)
    """オイラー法
    ・参考にしました: https://twitter.com/genkuroki/status/1301832571131633665/photo/1
    """

    t = range(t_span..., step = Δt)  # 時間軸
    x = Vector{typeof(x₀)}(undef, length(t))  # 解を格納する1次元配列

    x[1] = x₀  # 初期値
    for i in 1:length(x)-1
        x[i+1] = x[i] .+ dx(t[i], x[i]) .* Δt
    end

    t, x
end

function jisaku_solve_RungeKutta(dx, x₀, t_span, Δt)
    """ルンゲクッタ法（4次）"""

    t = range(t_span..., step = Δt)  # 時間軸
    x = Vector{typeof(x₀)}(undef, length(t))  # 解を格納する1次元配列

    x[1] = x₀  # 初期値
    for i in 1:length(x)-1
        k₁ = dx(t[i], x[i])
        k₂ = dx(t[i]+Δt/2, x[i]+k₁*Δt/2)
        k₃ = dx(t[i]+Δt/2, x[i]+k₂*Δt/2)
        k₄ = dx(t[i]+Δt, x[i]+k₃*Δt)
        x[i+1] = x[i] .+ (k₁ .+ 2k₂ .+ 2k₃ .+k₄) .* Δt/6
    end

    t, x
end

function VanDerPol(t, x)
    K = 1.0
    [
        x[2]
        K * (1 - x[1]^2) * x[2] - x[1]
    ]
end

# 数値シミュレーション実行
x₀ = [0.1, 0.1]  # 初期値
t_span = (0.0, 50.0)  # 時間幅
Δt = 0.01  # 刻み時間
#@time t, x = jisaku_solve_euler(VanDerPol, x₀, t_span, Δt)  # 解く
@time t, x = jisaku_solve_RungeKutta(VanDerPol, x₀, t_span, Δt)  # 解く

# plot
x, v = split_vec_of_arrays(x)  # 解をplotしやすいように分割
fig1 = plot(
    t, x,
    #linewidth=3,
    xaxis="time",
    label="x",
)
plot!(fig1, t, v, label = "ẋ")

fig2 = plot(
    x, v,
    #linewidth=3,
    xaxis="x",
    yaxis="dx",
    label="trajectory",
)

fig3 = plot(
    t, x, v,
    xaxis="time", yaxis="x", zaxis="v",
    label="trajectory"
)

plot(
    fig1, fig2, fig3, layout=(3,1),
    size=(600,900)
)


text:
function find_passage(r::Robot,side::HorizonSide)
    num_steps=0; side=Ost
    while isborder(r,Nord)==true # прохода сверху нет
        num_steps+=1
        move!(r,side,num_steps)
        side=inverse(side)
    end
end

move!(r,side,num_steps)=for _ in 1:num_steps move!(r,side) end

inverse(side::HorizonSide) = HorizonSide(mod(Int(side)+2, 4))

max_stars_repo_name: Anpilogova/-

max_issues_repo_name: Anpilogova/-

max_forks_repo_name: Anpilogova/-


text:
include("common.jl")
@testset "SDIRK integrators" begin
    prob_inplace = prob_dde_1delay
    prob_notinplace = prob_dde_1delay_scalar_notinplace

    # ODE algorithms
    algs = [GenericImplicitEuler(), GenericTrapezoid(),
            ImplicitEuler(), ImplicitMidpoint(), Trapezoid(),
            TRBDF2(), SDIRK2(), SSPSDIRK2(),
            Kvaerno3(), KenCarp3(),
            Cash4(), Hairer4(), Hairer42(), Kvaerno4(), KenCarp4(),
            Kvaerno5(), KenCarp5()]

    @testset for alg in algs
        @show alg
        stepsalg = MethodOfSteps(alg)
        solve(prob_inplace, stepsalg)
        solve(prob_notinplace, stepsalg)
    end
end


text:
using Test

@testset "Zero shaves to zero" begin
    for T in [Float16,Float32,Float64]
        for k in -5:50
            A = zeros(T,2,3)
            Ar = shave(A,k)
            @test A == Ar
            @test zero(T) == round(zero(T),k)
        end
    end
end

@testset "one shaves to one" begin
    for T in [Float16,Float32,Float64]
        for k in 0:50
            A = ones(T,2,3)
            Ar = shave(A,k)
            @test A == Ar
            @test one(T) == round(one(T),k)
        end
    end
end

@testset "minus one shaves to minus one" begin
    for T in [Float16,Float32,Float64]
        for k in 0:50
            A = -ones(T,2,3)
            Ar = shave(A,k)
            @test A == Ar
            @test -one(T) == round(-one(T),k)
        end
    end
end

@testset "No (half)shaving/setting/grooming for keepbits=10,23,52" begin
    for (T,k) in zip([Float16,Float32,Float64],
                        [10,23,52])
        A = rand(T,200,300)
        Ar = shave(A,k)
        @test A == Ar

        Ar = halfshave(A,k)
        @test A == Ar

        Ar = set_one(A,k)
        @test A == Ar

        Ar = groom(A,k)
        @test A == Ar

        # and a single one
        r = rand(T)
        @test r == shave(r,k)
        @test r == halfshave(r,k)
        @test r == set_one(r,k)
    end
end

@testset "Approx equal for keepbits=5,10,25" begin
    for (T,k) in zip([Float16,Float32,Float64],
                        [6,11,26])
        A = rand(T,200,300)
        Ar = shave(A,k)
        @test A ≈ Ar

        Ar = halfshave(A,k)
        @test A ≈ Ar

        Ar = set_one(A,k)
        @test A ≈ Ar

        Ar = groom(A,k)
        @test A ≈ Ar

        # and a single one
        r = rand(T)
        @test r ≈ shave(r,k)
        @test r ≈ set_one(r,k)
        @test r ≈ halfshave(r,k)
    end
end

@testset "Idempotence" begin
    for T in [Float16,Float32,Float64]
        for k in 0:20
            A = rand(T,200,300)
            Ar = shave(A,k)
            Ar2 = shave(Ar,k)
            @test Ar == Ar2

            Ar = halfshave(A,k)
            Ar2 = halfshave(Ar,k)
            @test Ar == Ar2

            Ar = set_one(A,k)
            Ar2 = set_one(Ar,k)
            @test Ar == Ar2

            Ar = groom(A,k)
            Ar2 = groom(Ar,k)
            @test Ar == Ar2
        end
    end
end

@testset "Shave/set = round towards/away from zero?" begin
    N = 1000
    for _ in 1:N
        for (T,UIntT) in zip([Float16,Float32,Float64],
                                [UInt16,UInt32,UInt64])
            for k in 1:9
                x = randn(T)
                xr = shave(x,k)
                @test abs(xr) <= abs(x)

                xr = set_one(x,k)
                @test abs(xr) >= abs(x)
            end
        end
    end
end


text:
module BenchDivideAtMean
using Random
using BenchmarkTools
using TapirBenchmarks

Random.seed!(1234)

const SUITE = BenchmarkGroup()

for n in [2^20]
    SUITE[:n=>n] = s0 = BenchmarkGroup()
    xs = rand(n)
    s0[:impl=>:seq] = @benchmarkable divide_at_mean_seq($xs)
    s0[:impl=>:threads] = @benchmarkable divide_at_mean_threads($xs)
    s0[:impl=>:tapir] = @benchmarkable divide_at_mean_tapir($xs)
end

end
BenchDivideAtMean.SUITE