text:
# Kolmogorov distribution
# defined as the sup_{t \in [0,1]} |B(t)|, where B(t) is a Brownian bridge
# used in the Kolmogorov--Smirnov test for large n.

immutable Kolmogorov <: ContinuousUnivariateDistribution
end

insupport(::Kolmogorov, x::Real) = zero(x) <= x < Inf
insupport(::Type{Kolmogorov}, x::Real) = zero(x) <= x < Inf

mean(d::Kolmogorov) = 0.5*√2π*log(2.0)
var(d::Kolmogorov) = pi*pi/12.0 - 0.5*pi*log(2.0)^2
# TODO: higher-order moments also exist, can be obtained by differentiating series

# cdf and ccdf are based on series truncation.
# two different series are available, e.g. see:
#   N. Smirnov, "Table for Estimating the Goodness of Fit of Empirical Distributions",
#   The Annals of Mathematical Statistics , Vol. 19, No. 2 (Jun., 1948), pp. 279-281
#   http://projecteuclid.org/euclid.aoms/1177730256
# use one series for small x, one for large x
# 5 terms seems to be sufficient for Float64 accuracy
# some divergence from Smirnov's table in 6th decimal near 1.0 (e.g. 1.04): occurs in
# both series so assume error in table.

function cdf(d::Kolmogorov,x::Real)
    if x <= 0.0
        return 0.0
    elseif x > 1.0
        return 1.0-ccdf(d,x)
    end
    s = 0.0
    for i = 1:5
        s += exp(-((2*i-1)*π/x)^2/8.0)
    end
    √2π*s/x
end

function ccdf(d::Kolmogorov,x::Real)
    if x <= 1.0
        return 1.0-cdf(d,x)
    end
    s = 0.0
    for i = 1:5
        s += (iseven(i) ? -1 : 1)*exp(-2.0*(i*x)^2)
    end
    2.0*s
end

# TODO: figure out how best to truncate series
function pdf(d::Kolmogorov,x::Real)
    if x <= 0.0
        return 0.0
    elseif x <= 1.0
        c = π/(2.0*x)
        s = 0.0
        for i = 1:20
            k = ((2*i-1)*c)^2
            s += (k-1.0)*exp(-k/2.0)
        end
        return √2π*s/x^2
    else
        s = 0.0
        for i = 1:20
            s += (iseven(i) ? -1 : 1)*i^2*exp(-2.0*(i*x)^2)
        end
        return 8.0*x*s
    end
end

# Devroye IV.5

function rand(d::Kolmogorov)
    # p = cdf(d,t)
    if rand() < p
        # left interval
        while true
            # generate truncated gamma variate
            # g = rand(Truncated(Gamma(1.5,1.0),0.0,t))
            tp = pi*pi/(8*t*t)
            acc = false
            while !acc
                e0, e1 = Base.Random.randmtzig_exprnd(), Base.Random.randmtzig_exprnd()
                e0 = e= / (1- 1/(2*tp))
                e1 = 2*e1
                g = tp + e0
                acc = (e0*e0 <= tp*e1*(g+tp)) || (g/tp - 1 - log(g/tp) <= e1)
            end

            x = pi/sqrt(8g)
            w = 0
            z = 1/(2*g)
            p = exp(-g)
            n = 1
            q = 1
            u = rand()
            while u >= w
                w += z*q
                if u >= w
                    return w
                end
                n += 2
                q = p^(n*n-1)
                w -= n*n*q
            end
        end
    else
        e = Base.Random.randmtzig_exprnd()
        u = rand()
        x = sqrt(t*t+e/2)
        w = 0
        n = 1
        z = exp(-2*x*x)
        while u > w
            n += 1
            w += n*n*z^(n*n-1)
            if u >= w
                return x
            end
            n += 1
            w -= n*n*z^(n*n-1)
        end
    end
end


text:
push!(LOAD_PATH,"../src/")

using Documenter, DocumenterTools
using EAGO, IntervalArithmetic, MathOptInterface, McCormick, JuMP

import McCormick: final_cut, mid3v, precond_and_contract!, AbstractMCCallback, populate_affine!,
                  dline_seg, correct_exp!, cut, mid_grad, preconditioner_storage, newton, secant,
                  MCCallback, contract!, affine_exp!

import EAGO: ExtensionType, Evaluator, variable_dbbt!, set_current_node!,
             VariableInfo, Log, aggressive_filtering!,
             bool_indx_diff, trivial_filtering!, SIPResult, SIPProblem,
             GlobalOptimizer, InputProblem, ParsedProblem, is_integer_feasible_relaxed,
             local_problem_status, default_upper_heurestic, label_branch_variables!,
             label_fixed_variables!, AbstractDirectedGraph, AbstractCache,
             AbstractCacheAttribute, initialize!, f_init!, r_init!, fprop!, rprop!,
             Variable, Subexpression, Expression, Constant, Parameter

import EAGO.Script: dag_flattening!, register_substitution!, Template_Graph,
                    Template_Node, scrub, scrub!, flatten_expression!

const MOI = MathOptInterface

@info "Making documentation..."
makedocs(modules = [EAGO, McCormick],
         doctest = false,
         format = Documenter.HTML(
                prettyurls = get(ENV, "CI", nothing) == "true",
                canonical = "https://PSORLab.github.io/EAGO.jl/stable/",
                collapselevel = 1,
         ),
         authors = "Matthew E. Wilhelm",
         #repo = "https://github.com/PSORLab/EAGO.jl/blob/{commit}{path}#L{line}",
         sitename = "EAGO.jl: Easy Advanced Global Optimization",
         pages = Any["Introduction" => "index.md",
                     "McCormick Operator Library" => Any["McCormick/overview.md",
                                                         "McCormick/usage.md",
                                                         "McCormick/operators.md",
                                                         "McCormick/type.md",
                                                         "McCormick/implicit.md"
                                                        ],
                     "Global Optimizer" => Any["Optimizer/starting.md",
                                               "Optimizer/optimizer.md",
                                               "Optimizer/bnb_back.md",
                                               "Optimizer/relax_back.md",
                                               "Optimizer/domain_reduction.md",
                                               "Optimizer/high_performance.md",
                                               "Optimizer/udf_utilities.md"
                                               ],
                     "Semi-Infinite Programming" => "SemiInfinite/semiinfinite.md",
                     "Contributing to EAGO"      => Any["Dev/contributing.md",
                                                        "Dev/future.md"
                                                        ],
                     "References"  => "ref.md",
                     "Citing EAGO" => "cite.md"]
)

@info "Deploying documentation..."
deploydocs(repo = "github.com/PSORLab/EAGO.jl.git",
           push_preview  = true)

ext: jl

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


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

@doc raw"""StatefulSetStatus represents the current state of a StatefulSet.

    IoK8sApiAppsV1beta1StatefulSetStatus(;
        collisionCount=nothing,
        conditions=nothing,
        currentReplicas=nothing,
        currentRevision=nothing,
        observedGeneration=nothing,
        readyReplicas=nothing,
        replicas=nothing,
        updateRevision=nothing,
        updatedReplicas=nothing,
    )

    - collisionCount::Int32 : collisionCount is the count of hash collisions for the StatefulSet. The StatefulSet controller uses this field as a collision avoidance mechanism when it needs to create the name for the newest ControllerRevision.
    - conditions::Vector{IoK8sApiAppsV1beta1StatefulSetCondition} : Represents the latest available observations of a statefulset's current state.
    - currentReplicas::Int32 : currentReplicas is the number of Pods created by the StatefulSet controller from the StatefulSet version indicated by currentRevision.
    - currentRevision::String : currentRevision, if not empty, indicates the version of the StatefulSet used to generate Pods in the sequence [0,currentReplicas).
    - observedGeneration::Int64 : observedGeneration is the most recent generation observed for this StatefulSet. It corresponds to the StatefulSet's generation, which is updated on mutation by the API Server.
    - readyReplicas::Int32 : readyReplicas is the number of Pods created by the StatefulSet controller that have a Ready Condition.
    - replicas::Int32 : replicas is the number of Pods created by the StatefulSet controller.
    - updateRevision::String : updateRevision, if not empty, indicates the version of the StatefulSet used to generate Pods in the sequence [replicas-updatedReplicas,replicas)
    - updatedReplicas::Int32 : updatedReplicas is the number of Pods created by the StatefulSet controller from the StatefulSet version indicated by updateRevision.
"""
mutable struct IoK8sApiAppsV1beta1StatefulSetStatus <: SwaggerModel
    collisionCount::Any # spec type: Union{ Nothing, Int32 } # spec name: collisionCount
    conditions::Any # spec type: Union{ Nothing, Vector{IoK8sApiAppsV1beta1StatefulSetCondition} } # spec name: conditions
    currentReplicas::Any # spec type: Union{ Nothing, Int32 } # spec name: currentReplicas
    currentRevision::Any # spec type: Union{ Nothing, String } # spec name: currentRevision
    observedGeneration::Any # spec type: Union{ Nothing, Int64 } # spec name: observedGeneration
    readyReplicas::Any # spec type: Union{ Nothing, Int32 } # spec name: readyReplicas
    replicas::Any # spec type: Union{ Nothing, Int32 } # spec name: replicas
    updateRevision::Any # spec type: Union{ Nothing, String } # spec name: updateRevision
    updatedReplicas::Any # spec type: Union{ Nothing, Int32 } # spec name: updatedReplicas

    function IoK8sApiAppsV1beta1StatefulSetStatus(;collisionCount=nothing, conditions=nothing, currentReplicas=nothing, currentRevision=nothing, observedGeneration=nothing, readyReplicas=nothing, replicas=nothing, updateRevision=nothing, updatedReplicas=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1beta1StatefulSetStatus, Symbol("collisionCount"), collisionCount)
        setfield!(o, Symbol("collisionCount"), collisionCount)
        validate_property(IoK8sApiAppsV1beta1StatefulSetStatus, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        validate_property(IoK8sApiAppsV1beta1StatefulSetStatus, Symbol("currentReplicas"), currentReplicas)
        setfield!(o, Symbol("currentReplicas"), currentReplicas)
        validate_property(IoK8sApiAppsV1beta1StatefulSetStatus, Symbol("currentRevision"), currentRevision)
        setfield!(o, Symbol("currentRevision"), currentRevision)
        validate_property(IoK8sApiAppsV1beta1StatefulSetStatus, Symbol("observedGeneration"), observedGeneration)
        setfield!(o, Symbol("observedGeneration"), observedGeneration)
        validate_property(IoK8sApiAppsV1beta1StatefulSetStatus, Symbol("readyReplicas"), readyReplicas)
        setfield!(o, Symbol("readyReplicas"), readyReplicas)
        validate_property(IoK8sApiAppsV1beta1StatefulSetStatus, Symbol("replicas"), replicas)
        setfield!(o, Symbol("replicas"), replicas)
        validate_property(IoK8sApiAppsV1beta1StatefulSetStatus, Symbol("updateRevision"), updateRevision)
        setfield!(o, Symbol("updateRevision"), updateRevision)
        validate_property(IoK8sApiAppsV1beta1StatefulSetStatus, Symbol("updatedReplicas"), updatedReplicas)
        setfield!(o, Symbol("updatedReplicas"), updatedReplicas)
        o
    end
end # type IoK8sApiAppsV1beta1StatefulSetStatus

const _property_map_IoK8sApiAppsV1beta1StatefulSetStatus = Dict{Symbol,Symbol}(Symbol("collisionCount")=>Symbol("collisionCount"), Symbol("conditions")=>Symbol("conditions"), Symbol("currentReplicas")=>Symbol("currentReplicas"), Symbol("currentRevision")=>Symbol("currentRevision"), Symbol("observedGeneration")=>Symbol("observedGeneration"), Symbol("readyReplicas")=>Symbol("readyReplicas"), Symbol("replicas")=>Symbol("replicas"), Symbol("updateRevision")=>Symbol("updateRevision"), Symbol("updatedReplicas")=>Symbol("updatedReplicas"))
const _property_types_IoK8sApiAppsV1beta1StatefulSetStatus = Dict{Symbol,String}(Symbol("collisionCount")=>"Int32", Symbol("conditions")=>"Vector{IoK8sApiAppsV1beta1StatefulSetCondition}", Symbol("currentReplicas")=>"Int32", Symbol("currentRevision")=>"String", Symbol("observedGeneration")=>"Int64", Symbol("readyReplicas")=>"Int32", Symbol("replicas")=>"Int32", Symbol("updateRevision")=>"String", Symbol("updatedReplicas")=>"Int32")
Base.propertynames(::Type{ IoK8sApiAppsV1beta1StatefulSetStatus }) = collect(keys(_property_map_IoK8sApiAppsV1beta1StatefulSetStatus))
Swagger.property_type(::Type{ IoK8sApiAppsV1beta1StatefulSetStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1beta1StatefulSetStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1beta1StatefulSetStatus }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1beta1StatefulSetStatus[property_name]

function check_required(o::IoK8sApiAppsV1beta1StatefulSetStatus)
    (getproperty(o, Symbol("replicas")) === nothing) && (return false)
    true
end

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


text:
relations = Dict{String, Bool}(
  "contains"           => true,
  "contains_properly"  => true,
  "covered_by"         => true,
  "covers"             => true,
  "crosses"            => true,
  "disjoint"           => true,
  "intersects"         => true,
  "overlaps"           => true,
  "touches"            => true,
  "within"             => true
)

formats = Dict{String, Bool}(
  "legacy"                   => true,
  "geojson"                  => true,
  "view"                     => true,
  "application/vnd.geo+json" => true
)

# WKT = Dict{String, Bool}(
#   "point"              => true,
#   "linestring"         => true,
#   "polygon"            => true,
#   "multipoint"         => true,
#   "multilinestring"    => true,
#   "multipolygon"       => true,
#   "geometrycollection" => true
# )

"""
    result = geo_index(db::Database, ddoc::AbstractString, name::AbstractString, index::AbstractString)

Create a geospatial index.

The `index` parameter is a string containing an index function in Javascript.

### Examples

    result = geo_index(db, "geodd", "geoidx",
      "function(doc){if(doc.geometry&&doc.geometry.coordinates){st_index(doc.geometry);}}"
    )

### Returns

Returns a `Dict(...)` from the CouchDB response, of the type

    Dict(
      "ok"  => true,
      "rev" => "1-b950984b19bb1b8bb43513c9d5b235bc",
      "id"  => "_design/geodd"
    )

[API reference](https://docs.cloudant.com/geo.html)
"""
function geo_index(db::Database, ddoc::AbstractString, name::AbstractString, index::AbstractString)
  data = Dict("st_indexes" => Dict(name => Dict("index" => index)))
  relax(HTTP.put, endpoint(db.url, "_design/$ddoc"); json=data, cookies=db.client.cookies)
end

"""
    result = geo_indexinfo(db::Database, ddoc::AbstractString, name::AbstractString)

Retrieve stats for a geospatial index.

### Examples

    result = geo_indexinfo(db, "geodd", "geoidx")

### Returns

Returns a `Dict(...)` from the CouchDB response, of the type

    Dict(
      "name" => "_design/geodd/geoidx",
      "geo_index" => Dict(
        "doc_count" => 269,
        "disk_size" => 33416,
        "data_size" => 26974
      )
    )

[API reference](https://docs.cloudant.com/geo.html#obtaining-information-about-a-cloudant-geo-index)
"""
function geo_indexinfo(db::Database, ddoc::AbstractString, name::AbstractString)
  relax(HTTP.get, endpoint(db.url, "_design/$ddoc/_geo_info/$name"), cookies=db.client.cookies)
end

"""
    geo_query(db::Database, ddoc::AbstractString, name::AbstractString;
      lat::Float64    = -360.0,
      lon::Float64    = -360.0,
      rangex::Float64 = 0.0,
      rangey::Float64 = 0.0,
      radius::Float64 = 0.0,
      bbox::Vector{Float64}  = Vector{Float64}(),
      relation::AbstractString = "intersects",
      nearest = false,
      bookmark::AbstractString = "",
      format::AbstractString = "view",
      skip = 0,
      limit = 0,
      stale = false,
      g::AbstractString = "")

Query a geospatial index. This quickly becomes complicated. See the references below.

The "g" parameter is a string representing a [`Well Known Text`](https://en.wikipedia.org/wiki/Well-known_text)
object (`WKT`). It can be used to describe various geometries, such as lines and polygons. Currently supported
geometric objects are

* point
* linestring
* polygon
* multipoint
* multilinestring
* multipolygon
* geometrycollection

Geo queries can be configured to return its results in a number of different formats
using the `format` parameter. The accepted values are:

* legacy
* geojson
* view (default)
* application/vnd.geo+json

The `relation` parameter follows the [DE-9IM](https://en.wikipedia.org/wiki/DE-9IM)
spec for geometric relationships. Acceptable values are:

* contains
* contains_properly
* covered_by
* covers
* crosses
* disjoint
* intersects (default)
* overlaps
* touches
* within

### Examples

Radial query

    result = geo_query(geodb, "geodd", "geoidx";
      lat    = 42.357963,
      lon    = -71.063991,
      radius = 10000.0,
      limit  = 200)

Polygon query

    result = geo_query(geodb, "geodd", "geoidx";
      g="POLYGON ((-71.0537124 42.3681995 0,-71.054399 42.3675178 0,-71.0522962 42.3667409 0,-71.051631 42.3659324 0,-71.051631 42.3621431 0,-71.0502148 42.3618577 0,-71.0505152 42.3660275 0,-71.0511589 42.3670263 0,-71.0537124 42.3681995 0))")

[API reference](https://docs.cloudant.com/geo.html#querying-a-cloudant-geo-index)
"""

function geo_query(db::Database, ddoc::AbstractString, name::AbstractString;
  lat::Float64    = -360.0, # -90:90 valid range (degrees)
  lon::Float64    = -360.0, # -180:180 valid range (degrees)
  rangex::Float64 = 0.0,
  rangey::Float64 = 0.0,
  radius::Float64 = 0.0,
  bbox::Vector{Float64}  = Vector{Float64}(),
  relation::AbstractString = "intersects",
  nearest = false,
  bookmark::AbstractString = "",
  format::AbstractString = "view", # legacy | geojson | view | application/vnd.geo+json
  skip = 0,
  limit = 0,
  stale = false,
  g::AbstractString = "") # We really need GeoJSON to do this properly

  query::Dict{String, Any} = Dict()

  if lat != -360.0
    query["lat"] = lat
  end

  if lon != -360.0
    query["lon"] = lon
  end

  if rangex > 0.0
    query["rangex"] = rangex
  end

  if rangey > 0.0
    query["rangey"] = rangey
  end

  if radius > 0.0
    query["radius"] = radius
  end

  if g != "" # We really need GeoJSON to do this properly
    query["g"] = g
  end

  if length(bbox) == 4
    query["bbox"] = bbox
  end

  if !haskey(relations, lowercase(relation))
    error("Unknown relation: '$relation'")
  end
  if relation != "intersects"
    query["relation"] = relation
  end

  if !haskey(formats, lowercase(format))
    error("Unknown format: '$format'")
  end
  if format != "view"
    query["format"] = format
  end

  if bookmark != ""
    query["bookmark"] = bookmark
  end

  if nearest
    query["nearest"] = true
  end

  if stale
    query["stale"] = "ok" # ffs
  end

  if skip > 0
    query["skip"] = skip
  end

  if limit > 0
    query["limit"] = limit
  end

  relax(HTTP.get, endpoint(db.url, "_design/$ddoc/_geo/$name"); cookies=db.client.cookies, query=query)
end

max_stars_repo_name: xpqz/couchzilla

max_issues_repo_name: xpqz/couchzilla

max_forks_repo_name: xpqz/couchzilla


text:
"""
KitBase.jl: The lightweight module of physical formulations in Kinetic.jl

Copyright (c) 2020-2022 Tianbai Xiao <tianbaixiao@gmail.com>
"""

module KitBase

if VERSION < v"1.3"
    @warn "Kinetic works better with Julia 1.3 or newer versions."
end

export KB

import Base: *, @kwdef
import BSON
import JLD2
using Base.Threads: @threads
using CSV
using CUDA
using Dates
using DiffEqOperators
using Distributed
using Distributions
using FastGaussQuadrature
using FFTW
using FileIO
using ForwardDiff
using LinearAlgebra
using MultivariatePolynomials
using OffsetArrays
using Optim
using Printf
using PyCall
using RecipesBase
using Reexport
using SpecialFunctions
using StaticArrays
using StructArrays
using TypedPolynomials
using WriteVTK
@reexport using FiniteMesh
using FiniteMesh.DocStringExtensions
using FiniteMesh.ProgressMeter

include("Data/data.jl")
include("Macro/macro.jl")
include("Type/type.jl")
include("IO/io.jl")
include("Math/math.jl")
include("Geometry/geometry.jl")
include("Theory/theory.jl")
include("Phase/phase.jl")
include("Reconstruction/reconstruction.jl")
include("Flux/flux.jl")
include("AD/ad.jl")
include("Config/config.jl")
include("Boundary/boundary.jl")
include("Solver/solver.jl")

const KB = KitBase
const itp = PyNULL()

function __init__()
    np = nworkers()
    nt = Threads.nthreads()
    if nt > 1 || np > 1
        @info "Kinetic will run with $np processors and $nt threads"
    else
        @info "Kinetic will run serially"
    end

    if has_cuda()
        @info "Kinetic will run with CUDA"
        for (i, dev) in enumerate(CUDA.devices())
            @info "$i: $(CUDA.name(dev))"
        end
        @info "Scalar operation is disabled in CUDA"
        CUDA.allowscalar(false)
    end

    copy!(itp, pyimport("scipy.interpolate"))
end

end


text:
using Test
using Mimi

@testset "Packagation" begin

    import MimiPAGE2020

    m = MimiPAGE2020.getpage()
    run(m)

    ## Just test that this was successful
    @test typeof(m) <: Mimi.Model

end


text:
using Test, Awabi

@testset "mecabrc" begin
    @test Awabi.find_mecabrc() == "/etc/mecabrc"
    mecabrc_map = Awabi.get_mecabrc_map(Awabi.find_mecabrc())
    @test mecabrc_map["dicdir"] == "/var/lib/mecab/dic/debian"
    @test Awabi.get_dic_path(mecabrc_map, "sys.dic") == "/var/lib/mecab/dic/debian/sys.dic"
end

@testset "dic" begin
    # Matrix
    m = Awabi.get_matrix(
        Awabi.get_dic_path(Awabi.get_mecabrc_map(Awabi.find_mecabrc()), "matrix.bin")
    )
    @test Awabi.get_trans_cost(m, UInt16(555), UInt16(1283)) == 340
    @test Awabi.get_trans_cost(m, UInt16(10), UInt16(1293)) == -1376

    # CharInfo
    cp = Awabi.get_char_propery(
        Awabi.get_dic_path(Awabi.get_mecabrc_map(Awabi.find_mecabrc()), "char.bin")
    )
    @test cp.category_names == [
        b"DEFAULT", b"SPACE", b"KANJI", b"SYMBOL", b"NUMERIC", b"ALPHA",
        b"HIRAGANA", b"KATAKANA", b"KANJINUMERIC", b"GREEK", b"CYRILLIC"
    ]
    @test Awabi.get_char_info(cp, UInt16(0)) == Awabi.CharInfo(0, 1, 0, 1, 0)         # DEFAULT
    @test Awabi.get_char_info(cp, UInt16(0x20)) == Awabi.CharInfo(1, 2, 0, 1, 0)      # SPACE
    @test Awabi.get_char_info(cp, UInt16(0x09)) == Awabi.CharInfo(1, 2, 0, 1, 0)      # SPACE
    @test Awabi.get_char_info(cp, UInt16(0x6f22)) == Awabi.CharInfo(2, 4, 2, 0, 0)    # KANJI 漢
    @test Awabi.get_char_info(cp, UInt16(0x3007)) == Awabi.CharInfo(3, 264, 0, 1, 1)  # SYMBOL
    @test Awabi.get_char_info(cp, UInt16(0x31)) == Awabi.CharInfo(4, 16, 0, 1, 1)     # NUMERIC 1
    @test Awabi.get_char_info(cp, UInt16(0x3042)) == Awabi.CharInfo(6, 64, 2, 1, 0)   # HIRAGANA あ
    @test Awabi.get_char_info(cp, UInt16(0x4e00)) == Awabi.CharInfo(8, 260, 0, 1, 1)  # KANJINUMERIC 一

    # lookup
    sys_dic = Awabi.get_mecabdic(
        Awabi.get_dic_path(Awabi.get_mecabrc_map(Awabi.find_mecabrc()), "sys.dic")
    )
    s = Vector{UInt8}("すもももももももものうち")
    @test length(Awabi.common_prefix_search(sys_dic, s)) == 3
    @test length(Awabi.lookup(sys_dic, s)) == 9
    s = Vector{UInt8}("もももももも")
    @test length(Awabi.common_prefix_search(sys_dic, s)) == 2
    @test length(Awabi.lookup(sys_dic, s)) == 4

    # lookup_unknowns
    unk_dic = Awabi.get_mecabdic(
        Awabi.get_dic_path(Awabi.get_mecabrc_map(Awabi.find_mecabrc()), "unk.dic")
    )
    @test Awabi.exact_match_search(unk_dic, Vector{UInt8}("SPACE")) == Int32(9729)
    entries, invoke = Awabi.lookup_unknowns(unk_dic, Vector{UInt8}("１９６７年"), cp)
    @test entries[1].original == b"１９６７"
end

@testset "tokenize" begin
    results = [
        [
            ("すもも", "名詞,一般,*,*,*,*,すもも,スモモ,スモモ"),
            ("も", "助詞,係助詞,*,*,*,*,も,モ,モ"),
            ("もも", "名詞,一般,*,*,*,*,もも,モモ,モモ"),
            ("も", "助詞,係助詞,*,*,*,*,も,モ,モ"),
            ("もも", "名詞,一般,*,*,*,*,もも,モモ,モモ"),
            ("の", "助詞,連体化,*,*,*,*,の,ノ,ノ"),
            ("うち", "名詞,非自立,副詞可能,*,*,*,うち,ウチ,ウチ"),
        ],
        [
            ("すもも", "名詞,一般,*,*,*,*,すもも,スモモ,スモモ"),
            ("も", "助詞,係助詞,*,*,*,*,も,モ,モ"),
            ("もも", "名詞,一般,*,*,*,*,もも,モモ,モモ"),
            ("もも", "名詞,一般,*,*,*,*,もも,モモ,モモ"),
            ("も", "助詞,係助詞,*,*,*,*,も,モ,モ"),
            ("の", "助詞,連体化,*,*,*,*,の,ノ,ノ"),
            ("うち", "名詞,非自立,副詞可能,*,*,*,うち,ウチ,ウチ"),
        ],
        [
            ("すもも", "名詞,一般,*,*,*,*,すもも,スモモ,スモモ"),
            ("もも", "名詞,一般,*,*,*,*,もも,モモ,モモ"),
            ("も", "助詞,係助詞,*,*,*,*,も,モ,モ"),
            ("もも", "名詞,一般,*,*,*,*,もも,モモ,モモ"),
            ("も", "助詞,係助詞,*,*,*,*,も,モ,モ"),
            ("の", "助詞,連体化,*,*,*,*,の,ノ,ノ"),
            ("うち", "名詞,非自立,副詞可能,*,*,*,うち,ウチ,ウチ"),
        ],
    ]

    tokenizer = Tokenizer()
    @test tokenize(tokenizer, "すもももももももものうち") == results[1]
    @test tokenize_n_best(tokenizer, "すもももももももものうち", 3) == results
end


text:
#=
    * Need to consider Duration
    * Activeted levels for career paths are not updated
=#
requiredTypes = [ "AcademicLevel", "Affiliation", "Assignment", "CareerPath", "CareerStatus", "Job", "Rank", "SubJob" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include(  reqType * ".jl" )
    end  # if !isdefined( Symbol( ...
end

using LightGraphs

export ConstructCareerPaths
function ConstructCareerPaths(g::DiGraph, DefinedLevels::Array{Vector{AbstractLevel}})
    CareerPaths = Vector{CareerPath}()

		if size(DefinedLevels, 1) < 1
			warn("No level defined.")
			return CareerPaths
		end
		InitPaths = CreateInitPaths(g, DefinedLevels[1])

		for path in InitPaths
			push!(CareerPaths, path)
		end

        for i in 2:length(DefinedLevels)
            NextLevelPaths!(CareerPaths, DefinedLevels[i-1], DefinedLevels[i])
        end

        for i in 1:length(CareerPaths)
            for j in 1:length(DefinedLevels)
                ActivateLevel(CareerPaths[i], typeof(DefinedLevels[j][1]))
            end
        end
    return CareerPaths
end

function CreatePathFromGraph(g::DiGraph, n::Int)
	CarPa = Vector{Array{Int}}()
  if length(neighbors(g, n)) == 0
      push!(CarPa, [-1])
  else
      for neighbor in neighbors(g, n)
          tempPaths = CreatePathFromGraph(g, neighbor)
          for path in tempPaths
              push!(CarPa, [n path])
          end
      end
  end
  return CarPa
end

function CreateInitPaths(g::DiGraph, StLevel::Array{AbstractLevel})

	# Traversing the graph
	InitPaths = Vector{Array{Int}}()
	if is_cyclic(g)
			warn("The graph is cyclic.")
			return Array{CareerPaths}()
	else
			for node in neighbors(g, 1)
					tmpCarPat = CreatePathFromGraph(g, node)
					for path in tmpCarPat
							push!(InitPaths, path)
					end
			end
	end
	InitCareerPaths = Array{CareerPath}(length(InitPaths))
	for path in 1:length(InitPaths)
		InitCareerPaths[path] = CareerPath()
		for p in InitPaths[path][1:(end-1)]
				push!(InitCareerPaths[path].Path, CareerStatus(academicLevel = StLevel[p-1]))
		end
	end
	#println("CreateInitPaths ", InitCareerPaths)
	return InitCareerPaths
end

SetterGetter = Dict{Any,Any}(
                    AcademicLevel => (get_AcademicLevel, set_AcademicLevel),
                    Rank => (get_Rank, set_Rank),
                    Assignment => (get_Assignment, set_Assignment),
                    Affiliation => (get_Affiliation, set_Affiliation),
                    Job => (get_Job, set_Job),
                    SubJob => (get_SubJob, set_SubJob))

function NextLevelPaths!(CareerPaths::Vector{CareerPath}, CurrentLevel::Vector{AbstractLevel}, NextLevel::Vector{AbstractLevel})
    # Getter and Setter for the current and next level
    #currentGet, currentSet = SetterGetter[typeof(CurrentLevel[1])]
    #nextGet, nextSet = SetterGetter[typeof(NextLevel[1])]

    if typeof(NextLevel[1]) == Rank
        Limits = sum(getMinStay(FirstLevel[i]) for i in 1:length(FirstLevel)) + sum(getMaxStay(FirstLevel[i]) for i in 1:length(FirstLevel))
        if Limits == 0
            ConstructRankLevelPathsFreeStLvl(CareerPaths, CurrentLevel, NextLevel)
        else
            ConstructRankLevelPathsLimitedDuration(CareerPaths, CurrentLevel, NextLevel)
        end
    else
        warn("$(typeof(NextLevel[1])) is under Construction.")
    end
end

function ConstructRankLevelPathsFreeStLvl(CareerPaths::Vector{CareerPath}, CurrentLevel::Vector{AbstractLevel}, NextLevel::Vector{AbstractLevel})
    currentGet, currentSet = SetterGetter[typeof(CurrentLevel[1])]
    nextGet, nextSet = SetterGetter[typeof(NextLevel[1])]
    # This is the number of paths that we will exploit
    nbInitPaths = length(CareerPaths)
    for i in 1:nbInitPaths
        PathToExploit = shift!(CareerPaths) # This path is used to generate the possible career paths which will be appended to CareerPaths at the end ***---------
        GeneratedPaths = Vector{CareerPath}()

        #for the first node in the career path we generate possible sub paths
        workinglevel = currentGet(PathToExploit.Path[1])
        for k in 1:length(workinglevel.NextLevel)
            tmppath = CareerPath()
            for l in 1:k
                tmpstatus = CareerStatus( academicLevel = PathToExploit.Path[1].AcademicLevel,
                                          rank = PathToExploit.Path[1].Rank,
                                          assignment = PathToExploit.Path[1].Assignment,
                                          affiliation = PathToExploit.Path[1].Affiliation,
                                          job = PathToExploit.Path[1].Job,
                                          subjob = PathToExploit.Path[1].SubJob)

                nextSet(tmpstatus, NextLevel[findIndex(workinglevel.NextLevel[l], NextLevel)])
                push!(tmppath.Path, tmpstatus)
            end
            push!(GeneratedPaths, tmppath)
        end

        # for each of the remaining nodes in the career path we generate the possible sub paths
        for j in 2:length(PathToExploit.Path)
            workinglevel = currentGet(PathToExploit.Path[j])
            nbPaths = length(GeneratedPaths)
            for p in 1:nbPaths
                workingpath = shift!(GeneratedPaths)

                if nextGet(workingpath.Path[end]).RankCode <= NextLevel[findIndex(workinglevel.NextLevel[1], NextLevel)].RankCode
                    println("StartingIndex = 1")
                    StartingIndex = 1
                elseif nextGet(workingpath.Path[end]).RankCode > NextLevel[findIndex(workinglevel.NextLevel[end], NextLevel)].RankCode
                    continue
                else
                    for indx in 1:length(workinglevel.NextLevel)
                        if nextGet(workingpath.Path[end]).Name == workinglevel.NextLevel[indx]
                            StartingIndex = indx
                            break
                        end
                    end
                    println("StartingIndex = $StartingIndex")
                end
                for k in StartingIndex:length(workinglevel.NextLevel)
                    tmppath = CareerPath(workingpath)
                    for l in StartingIndex:k
                        tmpstatus = CareerStatus( academicLevel = PathToExploit.Path[j].AcademicLevel,
                                                  rank = PathToExploit.Path[j].Rank,
                                                  assignment = PathToExploit.Path[j].Assignment,
                                                  affiliation = PathToExploit.Path[j].Affiliation,
                                                  job = PathToExploit.Path[j].Job,
                                                  subjob = PathToExploit.Path[j].SubJob)

                        nextSet(tmpstatus, NextLevel[findIndex(workinglevel.NextLevel[l], NextLevel)])
                        push!(tmppath.Path, tmpstatus)
                    end
                    push!(GeneratedPaths, tmppath)
                end
            end
        end

        append!(CareerPaths, GeneratedPaths) # Appending generated sub career paths ***---------
    end
end

function ConstructRankLevelPathsLimitedDuration(CareerPaths::Vector{CareerPath}, CurrentLevel::Vector{AbstractLevel}, NextLevel::Vector{AbstractLevel})
    warn("$(typeof(NextLevel[1])) with first level restrictions is under Construction.")

end

function findIndex(name::String, Level::Vector{AbstractLevel})
    for i in 1:length(Level)
        (name == Level[i].Name) && return i
    end
    return -1
end

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
function build_optim_function(f; ad_backend=AD.ForwardDiffBackend())
    ∇f(x) = only(AD.gradient(ad_backend, f, x))
    return build_optim_function(f, ∇f; ad_backend)
end
function build_optim_function(f, ∇f; ad_backend=AD.ForwardDiffBackend())
    # because we need explicit access to grad, we generate these ourselves instead of using
    # GalacticOptim's auto-AD feature.
    # TODO: switch to caching API if available, see
    # https://github.com/JuliaDiff/AbstractDifferentiation.jl/issues/41
    function grad(res, x, p...)
        ∇fx = ∇f(x)
        @. res = -∇fx
        return res
    end
    function hess(res, x, p...)
        H = only(AD.hessian(ad_backend, f, x))
        @. res = -H
        return res
    end
    function hv(res, x, v, p...)
        Hv = only(AD.lazy_hessian(ad_backend, f, x) * v)
        @. res = -Hv
        return res
    end
    return SciMLBase.OptimizationFunction((x, p...) -> -f(x); grad, hess, hv)
end

function build_optim_problem(optim_fun, x₀)
    return SciMLBase.OptimizationProblem(optim_fun, x₀, nothing)
end

function optimize_with_trace(
    prob,
    optimizer;
    progress_name="Optimizing",
    progress_id=nothing,
    maxiters=1_000,
    callback=nothing,
    kwargs...,
)
    u0 = prob.u0
    fun = prob.f
    grad! = fun.grad
    # caches for the trace of x and f(x)
    xs = typeof(u0)[]
    fxs = typeof(fun.f(u0, nothing))[]
    ∇fxs = typeof(u0)[]
    iteration = 0
    function _callback(x, nfx, args...)
        ret = callback !== nothing && callback(x, nfx, args...)

        Base.@logmsg ProgressLogging.ProgressLevel progress_name progress =
            iteration / maxiters _id = progress_id

        iteration += 1
        # some backends mutate x, so we must copy it
        push!(xs, copy(x))
        push!(fxs, -nfx)
        # NOTE: GalacticOptim doesn't have an interface for accessing the gradient trace,
        # so we need to recompute it ourselves
        # see https://github.com/SciML/GalacticOptim.jl/issues/149
        push!(∇fxs, rmul!(grad!(similar(x), x, nothing), -1))
        return ret
    end
    sol = GalacticOptim.solve(prob, optimizer; callback=_callback, maxiters, kwargs...)
    return sol, OptimizationTrace(xs, fxs, ∇fxs)
end

function optimize_with_trace(
    prob,
    optimizer::Union{Optim.FirstOrderOptimizer,Optim.SecondOrderOptimizer};
    progress_name="Optimizing",
    progress_id=nothing,
    maxiters=1_000,
    callback=nothing,
    kwargs...,
)
    iteration = 0
    function _callback(x, nfx, args...)
        ret = callback !== nothing && callback(x, nfx, args...)
        Base.@logmsg ProgressLogging.ProgressLevel progress_name progress =
            iteration / maxiters _id = progress_id
        iteration += 1
        return ret
    end
    new_kwargs = merge(NamedTuple(kwargs), (; store_trace=true, extended_trace=true))
    sol = GalacticOptim.solve(prob, optimizer; callback=_callback, maxiters, new_kwargs...)

    u0 = prob.u0
    xs = Vector{typeof(u0)}(undef, iteration)
    ∇fxs = Vector{typeof(u0)}(undef, iteration)
    result = sol.original
    copyto!(xs, Optim.x_trace(result))
    fxs = -Optim.f_trace(result)
    map!(tr -> -tr.metadata["g(x)"], ∇fxs, Optim.trace(result))

    return sol, OptimizationTrace(xs::Vector{typeof(u0)}, fxs, ∇fxs::Vector{typeof(u0)})
end

struct OptimizationTrace{P,L}
    points::P
    log_densities::L
    gradients::P
end

Base.length(trace::OptimizationTrace) = length(trace.points)

function Base.show(io::IO, trace::OptimizationTrace)
    print(io, "OptimizationTrace with $(length(trace) - 1) iterations")
    return nothing
end


text:
module mbaruniax1
using FinEtools
using FinEtoolsDeforLinear
using FinEtools.FENodeSetModule
using FinEtools.MeshExportModule
using FinEtoolsDeforLinear.MatDeforLinearElasticModule: tangentmoduli!, update!, thermalstrain!
using Test
import LinearAlgebra: norm, cholesky, cross
function test()
    E = 30e6*phun("psi") # Young's modulus
    alpha = 5.5e-6*phun("in")/phun("in")/phun("F")
    for nu in [0.0 0.3 0.45]

        MR = DeforModelRed1DStress
        material = MatDeforElastIso(MR,  0.0, E, nu, alpha)
        D = fill(0.0, 1, 1)
        t::FFlt, dt::FFlt, loc::FFltMat, label::FInt = 0.0, 0.0, [0.0 0.0 0.0], 0
        tangentmoduli!(material,  D,  t, dt, loc, label)
        @test D[1, 1] ≈ E
    end

    true
end
end
using .mbaruniax1
mbaruniax1.test()

module mbaruniax2
using FinEtools
using FinEtoolsDeforLinear
using FinEtools.FENodeSetModule
using FinEtools.MeshExportModule
using FinEtoolsDeforLinear.MatDeforLinearElasticModule: tangentmoduli!, update!, thermalstrain!
using Test
import LinearAlgebra: norm, cholesky, cross
function test()
    E = 30e6*phun("psi") # Young's modulus
    alpha = 5.5e-6*phun("in")/phun("in")/phun("F")
    for nu in [0.0 0.3 0.45]

        MR = DeforModelRed1DStrain
        material = MatDeforElastIso(MR,  0.0, E, nu, alpha)
        D = fill(0.0, 1, 1)
        t::FFlt, dt::FFlt, loc::FFltMat, label::FInt = 0.0, 0.0, [0.0 0.0 0.0], 0
        tangentmoduli!(material,  D,  t, dt, loc, label)
        lambda = E * nu / (1 + nu) / (1 - 2*(nu));
        mu = E / 2. / (1+nu);
        @test D[1, 1] ≈ lambda + 2*mu
    end

    true
end
end
using .mbaruniax2
mbaruniax2.test()


text:
function initialize_deep_weights_per_leyer(neurons,vars,a_all,nodes,a_layer1,n_weigths,x,y,obs,layer1,mat1)

    opti_ww=zeros(n_weigths)

    it1=0
    for i=1:neurons
        # for j=2:vars+1
        for j=2:length(a_all[i])
            it1+=1
            # opti_ww[it1]=a_all[i,j]
            opti_ww[it1]=a_all[i][j]
        end
    end
    println(Dates.format(now(), "HH:MM:SS")," ",it1," weights layer 1")

    layer_prev=copy(layer1)
    a_bias=Array{Float64}(undef,0)
    for i=2:length(nodes)
        layer_new=zeros(obs,neurons)

        # i1=layer_prev.>=0.9
        # layer_prev[i1].=0.9
        # i1=layer_prev.<=0.1
        # layer_prev[i1].=0.1
        # a_all_next=[ones(obs) layer_prev]\isigm1.(y[:,1])
        if i==2
            a_all_next=mat1*[layer_prev ones(obs)]'*isigm1.(y[:,1])
            # a_all_next=[layer_prev ones(obs)]\isigm1.(y[:,1])
        end
        if i==3
            # mat1=inv([layer_prev ones(obs)]'*[layer_prev ones(obs)])
            # a_all_next=mat1*[layer_prev ones(obs)]'*isigm1.(y[:,1])

            a_all_next=[layer_prev ones(obs)]\isigm1.(y[:,1])
        end
        layer_new=repeat(sigm1.([layer_prev ones(obs)]*a_all_next),1,neurons)

        for j=1:neurons # or nodes[i]
            push!(a_bias,a_all_next[end])
            for k=1:neurons
                it1+=1
                opti_ww[it1]=a_all_next[k]
            end
            # println("weight layer $(i) of $(length(nodes)), neuron $(j) of $(neurons)")
        end
        layer_prev=copy(layer_new)
        println(Dates.format(now(), "HH:MM:SS")," ",it1," weights layer ", i)
    end

    a_all_next=[layer_prev ones(obs)]\y
    # a_all_next=mat1*[layer_prev ones(obs)]'*y
    for i=1:neurons
        it1+=1
        opti_ww[it1]=a_all_next[i]
    end
    # for i=1:neurons
    #     it1+=1
    #     opti_ww[it1]=a_layer1[i]
    # end
    println(Dates.format(now(), "HH:MM:SS"),it1," weights layer out")

    # bias
    for i=1:neurons
        it1+=1
        # opti_ww[it1]=a_all[i,1]
        opti_ww[it1]=a_all[i][1]
    end
    println(it1," bias layer 1")

    it2=0
    for i=2:length(nodes)
        for j=1:nodes[i]
            it1+=1
            it2+=1
            opti_ww[it1]=a_bias[it2]
        end
        println(it1," bias layer ",i)
    end

    it1+=1
    opti_ww[it1]=a_all_next[end]
    # it1+=1
    # opti_ww[it1]=a_layer1[end]
    println(Dates.format(now(), "HH:MM:SS")," ",it1," bias layer out")

    return opti_ww,layer_prev

end


text:
using MarginalLogDensities
using Test
using Distributions
using Optim
using ForwardDiff
using LinearAlgebra
using HCubature
using Random

N = 3
σ = 1.5
d = MvNormal(zeros(N), σ*I)
logdensity(x) = logpdf(d, x)
im = [1, 3]
ij = [2]
dmarginal = MvNormal(zeros(length(ij)), σ*I)

@testset "Constructors" begin
    for forwarddiff_sparsity in [false, true]
        hp = HessianConfig(zeros(N, N), zeros(N), N, zeros(N, N), zeros(N, N), zeros(N), zeros(N))
        mld1 = MarginalLogDensity(logdensity, N, im, ij, LaplaceApprox(), hp)
        mld2 = MarginalLogDensity(logdensity, N, im)
        mld3 = MarginalLogDensity(logdensity, N, im, LaplaceApprox(), forwarddiff_sparsity)
        mld4 = MarginalLogDensity(logdensity, N, im,
            Cubature(-100ones(N), 100ones(N)))

        for mld in [mld2, mld3, mld4]
            @test dimension(mld) == dimension(mld1) == N
            @test imarginal(mld) == imarginal(mld1) == im
            @test ijoint(mld) == ijoint(mld1) == ij
            @test nmarginal(mld) == nmarginal(mld1) == length(imarginal(mld))
            @test njoint(mld) == njoint(mld1) == length(ijoint(mld))
            @test njoint(mld) + nmarginal(mld) == dimension(mld)
            @test isempty(setdiff(1:N, union(im, ij)))
        end
    end
end

# @testset "Sparse Hessians" begin
#     f(x) = -logdensity(x)
#     hconf = HessianConfig(f, im, ij)
#     H = zeros(length(im), length(im))

# end

@testset "Approximations" begin
    x = 1.0:3.0
    mld_laplace = MarginalLogDensity(logdensity, N, im, LaplaceApprox())
    mld_cubature = MarginalLogDensity(logdensity, N, im,
        Cubature(-100ones(2), 100ones(2)))

    @test mld_laplace(x[im], x[ij]) == logdensity(x)
    @test mld_cubature(x[im], x[ij]) == logdensity(x)
    # analytical: against 1D Gaussian
    @test mld_laplace(x[ij]) ≈ logpdf(dmarginal, x[ij])
    @test mld_cubature(x[ij]) ≈ logpdf(dmarginal, x[ij])
    # test against numerical integral
    int, err = hcubature(x -> exp(logdensity([x[1], 2, x[2]])),
        -100*ones(2), 100*ones(2))
    @test log(int) ≈ mld_laplace(x[ij])
    @test log(int) ≈ mld_cubature(x[ij])
    # marginalized density should be higher than joint density at same point
    @test mld_laplace(x[ij]) >= mld_laplace(x[im], x[ij])
    @test mld_cubature(x[ij]) >= mld_cubature(x[im], x[ij])
end


text:
using ApproxFun, Base.Test

## ODEs

d=Interval(-20000.,20000.)
x=Fun(identity,d)
u=[dirichlet(d);Derivative(d)^2+I]\[1.,0.]
u=[dirichlet(d);Derivative(d)^2+I]\[1.,0.]
@time u=[dirichlet(d);Derivative(d)^2+I]\[1.,0.]
println("Cos/Sin: should be ~0.017")

d=Interval(-1000.,5.)
x=Fun(identity,d)
u=[dirichlet(d);Derivative(d)^2-x]\[1.,0.]
u=[dirichlet(d);Derivative(d)^2-x]\[1.,0.]
@time u=[dirichlet(d);Derivative(d)^2-x]\[1.,0.]
println("Airy: 0.014356 seconds (1.08 k allocations: 8.015 MB)")

M=cache(ApproxFun.InterlaceOperator([dirichlet(d);Derivative(d)^2-x]);padding=true)
@time ApproxFun.resizedata!(M,12500,:)
println("Airy construct op: 0.005417 seconds (81 allocations: 5.279 MB)")

S=Chebyshev()
x=Fun(identity,S)
D=Derivative(S)
L=D^2+(7+2x+6x^2)
B=dirichlet(S)
n=20000
rhs=ones(n+2)
u=[B;L]\rhs
u=[B;L]\rhs
@time u=[B;L]\rhs
println("Poly: should be ~0.025")

S=Chebyshev()
x=Fun(identity,S)
D=Derivative(S)
L=D^2+cos(x)
B=dirichlet(S)
n=2000
rhs=ones(n+2)
u=linsolve([B;L],rhs;maxlength=Inf)
u=linsolve([B;L],rhs;maxlength=Inf)
@time u=linsolve([B;L],rhs;maxlength=Inf)
println("Cos: should be ~0.0075")

S=Chebyshev()
x=Fun(identity,S)
D=Derivative(S)
L=D^2+sin(x)
B=dirichlet(S)
n=2000
rhs=ones(n+2)
u=linsolve([B;L],rhs;maxlength=Inf)
u=linsolve([B;L],rhs;maxlength=Inf)
@time u=linsolve([B;L],rhs;maxlength=Inf)
println("Sin: should be ~0.008663 seconds (660 allocations: 2.987 MB)")

## Piecewise

x=Fun(identity,[-20.,-10.,-5.,0.,1.,15.])
sp=space(x)
D=Derivative(sp)

u=[dirichlet(sp);
    D^2-x]\[airyai(-20.)];
@time u=[dirichlet(sp);
    D^2-x]\[airyai(-20.)];

println("Piecewise Airy: should be ~0.008")

## Vector
d=Interval()
D=Derivative(d);
B=ldirichlet();
Bn=lneumann();

f=Fun(x->[exp(x),cos(x)],d)

A=[B 0;
   Bn 0;
   0 B;
   D^2-I 2.0I;
   0 D+I];

u=A\Any[0.,0.,0.,f]
@time u=A\Any[0.,0.,0.,f]
println("Systems: should be ~0.0008")

d=Interval(-300.,5.)
x=Fun(identity,d)
A=Derivative(d)^2-x
u=nullspace(A)
@test_approx_eq A[1:10,1:10] (A.')[1:10,1:10].'
@time u=nullspace(A)
println("Nullspace Airy: 0.052730 seconds (75.21 k allocations: 56.736 MB)")


text:
"""
    ==(a::RefMesh, b::RefMesh)

Test RefMesh equality.
"""
function Base.:(==)(a::T, b::T) where {T<:RefMesh}
    isequal(a.name, b.name) &&
        isequal(a.mesh, b.mesh) &&
        isequal(a.normals, b.normals) &&
        isequal(a.texture_coords, b.texture_coords) &&
        isequal(a.material, b.material) &&
        isequal(a.taper, b.taper)
end

"""
    ==(a::RefMeshes, b::RefMeshes)

Test RefMeshes equality.
"""
function Base.:(==)(a::T, b::T) where {T<:RefMeshes}
    isequal(names(a), names(b)) &&
        all([isequal(a.meshes[i], b.meshes[i]) for i in 1:length(a.meshes)])
end

"""
    ==(a::geometry, b::geometry)

Test RefMeshes equality.
"""
function Base.:(==)(a::T, b::T) where {T<:geometry}
    isequal(a.ref_mesh, b.ref_mesh) &&
        isequal(a.ref_mesh_index, b.ref_mesh_index) &&
        isequal(a.transformation([1; 1; 1]), b.transformation([1; 1; 1])) &&
        # NB: transform a point here because transformations can't be compared directly
        isequal(a.dUp, b.dUp) &&
        isequal(a.dDwn, b.dDwn) &&
        isequal(a.mesh, b.mesh)
end


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

module GitTools

using ..Pkg
using ..MiniProgressBars
import ..can_fancyprint, ..printpkgstyle
using SHA
import Base: SHA1
import LibGit2
using Printf

use_cli_git() = get(ENV, "JULIA_PKG_USE_CLI_GIT", "") == "true"

function transfer_progress(progress::Ptr{LibGit2.TransferProgress}, p::Any)
    progress = unsafe_load(progress)
    @assert haskey(p, :transfer_progress)
    bar = p[:transfer_progress]
    @assert typeof(bar) == MiniProgressBar
    if progress.total_deltas != 0
        bar.header = "Resolving Deltas:"
        bar.max = progress.total_deltas
        bar.current = progress.indexed_deltas
    else
        bar.max = progress.total_objects
        bar.current = progress.received_objects
    end
    show_progress(stdout, bar)
    return Cint(0)
end

const GIT_REGEX =
    r"^(?:(?<proto>git|ssh|https)://)?(?:[\w\.\+\-:]+@)?(?<hostname>.+?)(?(<proto>)/|:)(?<path>.+?)(?:\.git)?$"
const GIT_PROTOCOLS = Dict{String, Union{Nothing, String}}()
const GIT_USERS = Dict{String, Union{Nothing, String}}()

@deprecate setprotocol!(proto::Union{Nothing, AbstractString}) setprotocol!(protocol = proto) false

function setprotocol!(;
    domain::AbstractString="github.com",
    protocol::Union{Nothing, AbstractString}=nothing,
    user::Union{Nothing, AbstractString}=(protocol == "ssh" ? "git" : nothing)
)
    domain = lowercase(domain)
    GIT_PROTOCOLS[domain] = protocol
    GIT_USERS[domain] = user
end

function normalize_url(url::AbstractString)
    # LibGit2 is fussy about trailing slash. Make sure there is none.
    url = rstrip(url, '/')
    m = match(GIT_REGEX, url)
    m === nothing && return url

    host = m[:hostname]
    path = "$(m[:path]).git"

    proto = get(GIT_PROTOCOLS, lowercase(host), nothing)

    if proto === nothing
        url
    else
        user = get(GIT_USERS, lowercase(host), nothing)
        user = user === nothing ? "" : "$user@"

        "$proto://$user$host/$path"
    end
end

function ensure_clone(io::IO, target_path, url; kwargs...)
    if ispath(target_path)
        return LibGit2.GitRepo(target_path)
    else
        return GitTools.clone(io, url, target_path; kwargs...)
    end
end

function checkout_tree_to_path(repo::LibGit2.GitRepo, tree::LibGit2.GitObject, path::String)
    GC.@preserve path begin
        opts = LibGit2.CheckoutOptions(
            checkout_strategy = LibGit2.Consts.CHECKOUT_FORCE,
            target_directory = Base.unsafe_convert(Cstring, path)
        )
        LibGit2.checkout_tree(repo, tree, options=opts)
    end
end

function clone(io::IO, url, source_path; header=nothing, credentials=nothing, kwargs...)
    @assert !isdir(source_path) || isempty(readdir(source_path))
    url = normalize_url(url)
    printpkgstyle(io, :Cloning, header === nothing ? "git-repo `$url`" : header)
    bar = MiniProgressBar(header = "Fetching:", color = Base.info_color())
    transfer_payload = MiniProgressBar(header = "Fetching:", color = Base.info_color())
    fancyprint = can_fancyprint(io)
    callbacks = if fancyprint
        LibGit2.Callbacks(
            :transfer_progress => (
                @cfunction(transfer_progress, Cint, (Ptr{LibGit2.TransferProgress}, Any)),
                bar,
            )
        )
    else
        LibGit2.Callbacks()
    end
    fancyprint && start_progress(io, bar)
    if credentials === nothing
        credentials = LibGit2.CachedCredentials()
    end
    try
        if use_cli_git()
            run(`git clone --quiet $url $source_path`)
            return LibGit2.GitRepo(source_path)
        else
            mkpath(source_path)
            return LibGit2.clone(url, source_path; callbacks=callbacks, credentials=credentials, kwargs...)
        end
    catch err
        rm(source_path; force=true, recursive=true)
        err isa LibGit2.GitError || err isa InterruptException || rethrow()
        if err isa InterruptException
            Pkg.Types.pkgerror("git clone of `$url` interrupted")
        elseif (err.class == LibGit2.Error.Net && err.code == LibGit2.Error.EINVALIDSPEC) ||
           (err.class == LibGit2.Error.Repository && err.code == LibGit2.Error.ENOTFOUND)
            Pkg.Types.pkgerror("git repository not found at `$(url)`")
        else
            Pkg.Types.pkgerror("failed to clone from $(url), error: $err")
        end
    finally
        Base.shred!(credentials)
        fancyprint && end_progress(io, bar)
    end
end

function fetch(io::IO, repo::LibGit2.GitRepo, remoteurl=nothing; header=nothing, credentials=nothing, refspecs=[""], kwargs...)
    if remoteurl === nothing
        remoteurl = LibGit2.with(LibGit2.get(LibGit2.GitRemote, repo, "origin")) do remote
            LibGit2.url(remote)
        end
    end
    fancyprint = can_fancyprint(io)
    remoteurl = normalize_url(remoteurl)
    printpkgstyle(io, :Updating, header === nothing ? "git-repo `$remoteurl`" : header)
    bar = MiniProgressBar(header = "Fetching:", color = Base.info_color())
    fancyprint = can_fancyprint(io)
    callbacks = if fancyprint
        LibGit2.Callbacks(
            :transfer_progress => (
                @cfunction(transfer_progress, Cint, (Ptr{LibGit2.TransferProgress}, Any)),
                bar,
            )
        )
    else
        LibGit2.Callbacks()
    end
    fancyprint && start_progress(io, bar)
    if credentials === nothing
        credentials = LibGit2.CachedCredentials()
    end
    try
        if use_cli_git()
            cd(LibGit2.path(repo)) do
                run(`git fetch -q $remoteurl $(only(refspecs))`)
            end
        else
            return LibGit2.fetch(repo; remoteurl=remoteurl, callbacks=callbacks, refspecs=refspecs, kwargs...)
        end
    catch err
        err isa LibGit2.GitError || rethrow()
        if (err.class == LibGit2.Error.Repository && err.code == LibGit2.Error.ERROR)
            Pkg.Types.pkgerror("Git repository not found at '$(remoteurl)'")
        else
            Pkg.Types.pkgerror("failed to fetch from $(remoteurl), error: $err")
        end
    finally
        Base.shred!(credentials)
        fancyprint && end_progress(io, bar)
    end
end

# This code gratefully adapted from https://github.com/simonbyrne/GitX.jl
@enum GitMode mode_dir=0o040000 mode_normal=0o100644 mode_executable=0o100755 mode_symlink=0o120000 mode_submodule=0o160000
Base.string(mode::GitMode) = string(UInt32(mode); base=8)
Base.print(io::IO, mode::GitMode) = print(io, string(mode))

function gitmode(path::AbstractString)
    # Windows doens't deal with executable permissions in quite the same way,
    # `stat()` gives a different answer than we actually want, so we use
    # `isexecutable()` which uses `uv_fs_access()` internally.  On other
    # platforms however, we just want to check via `stat()`.
    function isexec(p)
        @static if Sys.iswindows()
            return Sys.isexecutable(p)
        end
        return !iszero(filemode(p) & 0o100)
    end
    if islink(path)
        return mode_symlink
    elseif isdir(path)
        return mode_dir
    elseif isexec(path)
        return mode_executable
    else
        return mode_normal
    end
end

"""
    blob_hash(HashType::Type, path::AbstractString)

Calculate the git blob hash of a given path.
"""
function blob_hash(::Type{HashType}, path::AbstractString) where HashType
    ctx = HashType()
    if islink(path)
        datalen = length(readlink(path))
    else
        datalen = filesize(path)
    end

    # First, the header
    SHA.update!(ctx, Vector{UInt8}("blob $(datalen)\0"))

    # Next, read data in in chunks of 4KB
    buff = Vector{UInt8}(undef, 4*1024)

    try
        if islink(path)
            update!(ctx, Vector{UInt8}(readlink(path)))
        else
            open(path, "r") do io
                while !eof(io)
                    num_read = readbytes!(io, buff)
                    update!(ctx, buff, num_read)
                end
            end
        end
    catch e
        if isa(e, InterruptException)
            rethrow(e)
        end
        @warn("Unable to open $(path) for hashing; git-tree-sha1 likely suspect")
    end

    # Finish it off and return the digest!
    return SHA.digest!(ctx)
end
blob_hash(path::AbstractString) = blob_hash(SHA1_CTX, path)

"""
    contains_files(root::AbstractString)

Helper function to determine whether a directory contains files; e.g. it is a
direct parent of a file or it contains some other directory that itself is a
direct parent of a file. This is used to exclude directories from tree hashing.
"""
function contains_files(path::AbstractString)
    st = lstat(path)
    ispath(st) || throw(ArgumentError("non-existent path: $(repr(path))"))
    isdir(st) || return true
    for p in readdir(path)
        contains_files(joinpath(path, p)) && return true
    end
    return false
end

"""
    tree_hash(HashType::Type, root::AbstractString)

Calculate the git tree hash of a given path.
"""
function tree_hash(::Type{HashType}, root::AbstractString; debug_out::Union{IO,Nothing} = nothing, indent::Int=0) where HashType
    entries = Tuple{String, Vector{UInt8}, GitMode}[]
    for f in sort(readdir(root; join=true); by = f -> isdir(f) ? f*"/" : f)
        # Skip `.git` directories
        if basename(f) == ".git"
            continue
        end

        filepath = abspath(f)
        mode = gitmode(filepath)
        if mode == mode_dir
            # If this directory contains no files, then skip it
            contains_files(filepath) || continue

            # Otherwise, hash it up!
            child_stream = nothing
            if debug_out !== nothing
                child_stream = IOBuffer()
            end
            hash = tree_hash(HashType, filepath; debug_out=child_stream, indent=indent+1)
            if debug_out !== nothing
                indent_str = "| "^indent
                println(debug_out, "$(indent_str)+ [D] $(basename(filepath)) - $(bytes2hex(hash))")
                print(debug_out, String(take!(child_stream)))
                println(debug_out, indent_str)
            end
        else
            hash = blob_hash(HashType, filepath)
            if debug_out !== nothing
                indent_str = "| "^indent
                mode_str = mode == mode_normal ? "F" : "X"
                println(debug_out, "$(indent_str)[$(mode_str)] $(basename(filepath)) - $(bytes2hex(hash))")
            end
        end
        push!(entries, (basename(filepath), hash, mode))
    end

    content_size = 0
    for (n, h, m) in entries
        content_size += ndigits(UInt32(m); base=8) + 1 + sizeof(n) + 1 + sizeof(h)
    end

    # Return the hash of these entries
    ctx = HashType()
    SHA.update!(ctx, Vector{UInt8}("tree $(content_size)\0"))
    for (name, hash, mode) in entries
        SHA.update!(ctx, Vector{UInt8}("$(mode) $(name)\0"))
        SHA.update!(ctx, hash)
    end
    return SHA.digest!(ctx)
end
tree_hash(root::AbstractString; debug_out::Union{IO,Nothing} = nothing) = tree_hash(SHA.SHA1_CTX, root; debug_out)

function check_valid_HEAD(repo)
    try LibGit2.head(repo)
    catch err
        Pkg.Types.pkgerror("invalid git HEAD ($(err.msg))")
    end
end

function git_file_stream(repo::LibGit2.GitRepo, spec::String; fakeit::Bool=false)::IO
    blob = try LibGit2.GitBlob(repo, spec)
    catch err
        err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow()
        fakeit && return devnull
    end
    iob = IOBuffer(LibGit2.content(blob))
    close(blob)
    return iob
end

end # module


text:
using DataStructures,Distributions,StatsBase,Random,Plots, LaTeXStrings; pyplot()

function simMM1Wait(lambda,mu,T)
    tNextArr = rand(Exponential(1/(lambda)))
    tNextDep = Inf
    t = tNextArr

    waitingRoom = Queue{Float64}()
    serverBusy = false
    waitTimes = Array{Float64,1}()

    while t<T
        if t == tNextArr
            if !serverBusy
                tNextDep = t + rand(Exponential(1/mu))
                serverBusy = true
                push!(waitTimes,0.0)
            else
                enqueue!(waitingRoom,t)
            end
            tNextArr = t + rand(Exponential(1/(lambda)))
        else
            if length(waitingRoom) == 0
               tNextDep = Inf
               serverBusy = false
            else
               tArr = dequeue!(waitingRoom)
               waitTime = t - tArr
               push!(waitTimes, waitTime)
               tNextDep = t + rand(Exponential(1/mu))
            end
        end
        t = min(tNextArr,tNextDep)
    end

    return waitTimes
end

Random.seed!(1)
lambda, mu = 0.8, 1.0
T = 10^3

data = simMM1Wait(lambda,mu,T)
empiricalCDF = ecdf(data)

F(x) = 1-(lambda/mu)*MathConstants.e^(-(mu-lambda)x)
xGrid = 0:0.1:20

plot(xGrid, F.(xGrid),
	c=:blue,label="Analytic CDF of waiting time")
plot!(xGrid, empiricalCDF(xGrid),
	c=:red,label="ECDF of waiting times",
	xlabel=L"x", ylabel=L"\Prob(W \leq x)", xlims=(0,20),ylims=(0,1),
	legend=:bottomright)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# Automatically generated using Clang.jl

# const CSFML_AUDIO_API = CSFML_API_IMPORT

struct sfTimeSpan
    offset::sfTime
    length::sfTime
end

const sfSoundRecorderStartCallback = Ptr{Cvoid}
const sfSoundRecorderProcessCallback = Ptr{Cvoid}
const sfSoundRecorderStopCallback = Ptr{Cvoid}

@cenum sfSoundStatus::UInt32 begin
    sfStopped = 0
    sfPaused = 1
    sfPlaying = 2
end

struct sfSoundStreamChunk
    samples::Ptr{sfInt16}
    sampleCount::UInt32
end

const sfSoundStreamGetDataCallback = Ptr{Cvoid}
const sfSoundStreamSeekCallback = Ptr{Cvoid}
const sfMusic = Cvoid
const sfSound = Cvoid
const sfSoundBuffer = Cvoid
const sfSoundBufferRecorder = Cvoid
const sfSoundRecorder = Cvoid
const sfSoundStream = Cvoid


text:
"""
    dimsRange, dimsDomain = size(op)

Return a tuple containg the dimensions of the operator `op`.
"""
size(op::AbstractLinearTransform) = (op.dimsRange, op.dimsDomain)
function size(op::AbstractLinearTransform, space::Symbol)
    if space == :domain
        op.dimsDomain
    elseif space == :range
        op.dimsRange
    else
        error("Unknown option")
    end
end

function size(op::AbstractLinearTransform, dim::Int)
    if dim == 1
        op.dimsRange
    elseif dim == 2
        op.dimsDomain
    else
        error("Unknown dim $dim")
    end
end

"""
    show(io, op)

Print a description of the generalized linear operator.
"""
function show(io::IO, op::AbstractLinearTransform)
    println(io,"Generalized linear transform")
    println(io,"  dim  range: ", op.dimsDomain)
    print(io,"  dim domain: ", op.dimsRange)
end

"""
    drandn(op::AbstractLinearTransform)

Return a pseudorandom vector in the domain of A.

See also [`rrandn`](@ref).
"""
drandn(op::AbstractLinearTransform) = randn(size(op,:domain))

"""
    rrandn(op::AbstractLinearTransform)

Return a pseudorandom vector in the domain of A.

See also [`drandn`](@ref).
"""
rrandn(op::AbstractLinearTransform) = randn(size(op,:range))

"""
    dottest(op::AbstractLinearTransform, [k=100])

Generate random vectors X and Y (from the domain and range of A), and verifies
that (A*X)'*Y = X'*(A'*Y) within a tolerance of 1E-10. This can help detect
errors in the operator; it canot be used to guarantee correctness.

The function returns true when the test succeeded and false if it failed.
"""
function dottest(A::AbstractLinearTransform, k=100; verbose=false)

    tol = √eps()
    err = 0.
    kpass = 0

    for i=1:k
        x   = drandn(A)
        y   = rrandn(A)
        z1  = dot(A*x, y)
        z2  = dot(x, A'*y)
        err = max(err, norm(z1-z2,Inf))
        if err < tol
            kpass += 1
        end
    end

    if verbose
        if kpass < k
            println("FAILED on $(k-kpass) out of $k tests")
            println("  maximum absolute difference of $err")
        else
            println("PASSED!");
        end
    end

    status = kpass != 0
end

"""
    Matrix(op::AbstractLinearTransform)

Return a matrix representation of the operator, if the sizes are
1-dimensional, and not tuples.
"""
function Matrix(op::AbstractLinearTransform)
    m, n = size(op)
    !(length(m) == length(n) == 1) && error("Dims cannot be tuples")
    T = eltype(op)
    mat = Matrix{T}(undef, (m, n))
    v = fill(zero(T), n)
    @inbounds for i = 1:n
        v[i] = one(T)
        mul!(view(mat, :, i), op, v)
        v[i] = zero(T)
    end
    return mat
end


text:
let
    md = CuModuleFile(joinpath(@__DIR__, "ptx/vadd.ptx"))

    vadd = CuFunction(md, "vadd")
end

# comparisons
let
    md = CuModuleFile(joinpath(@__DIR__, "ptx/vectorops.ptx"))
    vadd = CuFunction(md, "vadd")
    @test vadd == vadd
    @test vadd != CuFunction(md, "vdiv")
end

let
    f = open(joinpath(@__DIR__, "ptx/vadd.ptx"))
    ptx = read(f, String)
    close(f)

    md = CuModule(ptx)
    vadd = CuFunction(md, "vadd")

    md2 = CuModuleFile(joinpath(@__DIR__, "ptx/vadd.ptx"))
    @test md != md2
end

memcheck || @test_throws_cuerror CUDA.ERROR_INVALID_IMAGE CuModule("foobar")

@testset "globals" begin
    md = CuModuleFile(joinpath(@__DIR__, "ptx/global.ptx"))

    var = CuGlobal{Int32}(md, "foobar")
    @test eltype(var) == Int32
    @test eltype(typeof(var)) == Int32

    @test_throws ArgumentError CuGlobal{Int64}(md, "foobar")

    var[] = Int32(42)
    @test var[] == Int32(42)
end

@testset "linker" begin
    link = CuLink()
    @test link == link
    @test link != CuLink()

    # PTX code
    open(joinpath(@__DIR__, "ptx/empty.ptx")) do f
        add_data!(link, "vadd_parent", read(f, String))
    end
    @test_throws ArgumentError add_data!(link, "vadd_parent", "\0")

    # object code
    # TODO: test with valid object code
    # NOTE: apparently, on Windows cuLinkAddData! _does_ accept object data containing \0
    if !Sys.iswindows()
        memcheck || @test_throws_cuerror CUDA.ERROR_UNKNOWN add_data!(link, "vadd_parent", UInt8[0])
    end
end

let
    link = CuLink()
    add_file!(link, joinpath(@__DIR__, "ptx/vadd_child.ptx"), CUDA.JIT_INPUT_PTX)
    open(joinpath(@__DIR__, "ptx/vadd_parent.ptx")) do f
        add_data!(link, "vadd_parent", read(f, String))
    end

    obj = complete(link)
    md = CuModule(obj)

    vadd = CuFunction(md, "vadd")

    options = Dict{CUDA.CUjit_option,Any}()
    options[CUDA.JIT_GENERATE_LINE_INFO] = true

    md = CuModule(obj, options)
    vadd = CuFunction(md, "vadd")
end


text:
module Simulate_FourBar

using  Modia3D
import Modia3D.ModiaMath

# Solids
groundMaterial = Modia3D.Material(color="DarkGreen", transparency=0.5)
massProperties = "Aluminium"
vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)
vmat2 = Modia3D.Material(color="Red")

@assembly Bar(;Lx = 0.1, Ly=Lx/5, Lz=Ly) begin
   frame0 = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), nothing, vmat1))
   frame1 = Modia3D.Object3D(frame0; r=[-Lx/2, 0.0, 0.0])
   frame2 = Modia3D.Object3D(frame0; r=[ Lx/2, 0.0, 0.0])
   cyl    = Modia3D.Cylinder(Ly/2,1.2*Ly; material=vmat2)
   cyl1   = Modia3D.Object3D(frame1, cyl)
   cyl2   = Modia3D.Object3D(frame2, cyl)
end

@assembly Fourbar(;Lx = 0.1, Ly=Lx/5, Lz=Ly, groundWidth=Lx, groundHeight=0.1*Lx) begin
   world  = Modia3D.Object3D(Modia3D.CoordinateSystem(0.6))
   frame1 = Modia3D.Object3D(world ; r=[0.5*Lx, 0.0, groundWidth/2])
   frame2 = Modia3D.Object3D(frame1; r=[Lx    , 0.0, 0.0])
   ground = Modia3D.Object3D(world , Modia3D.Box(3*Lx, groundWidth, groundHeight; material=groundMaterial);
                                              r = [1.5*Lx, -groundHeight/2, groundWidth/2],
                                              R = ModiaMath.rot_nxy([-1,0,0], [0,0,1]))
   bar1   = Bar(Lx=Lx, Ly=Ly, Lz=Lz)
   bar2   = Bar(Lx=Lx, Ly=Ly, Lz=Lz)
   bar3   = Bar(Lx=Lx, Ly=Ly, Lz=Lz)

   rev1   = Modia3D.Revolute(frame1     , bar1.frame1; phi_start =  pi/2)
   rev2   = Modia3D.Revolute(bar1.frame2, bar2.frame1; phi_start = -pi/2)
   rev3   = Modia3D.Revolute(frame2     , bar3.frame1; phi_start =  pi/2)
   rev4   = Modia3D.Revolute(bar3.frame2, bar2.frame2; phi_start = -pi/2)
end

#hier wird eine scene erstellt
fourbar = Fourbar(Lx=1.0, sceneOptions=Modia3D.SceneOptions(visualizeFrames=true,
                                                    defaultFrameLength=0.3,
                                                    enableContactDetection=false))
#hier wird wieder eine scene erstellt
model = Modia3D.SimulationModel(fourbar)

#Modia3D.visualizeAssembly!( fourbar )

println("... success of Simulate_FourBar.jl!")
end


text:
# ---
# title: 1480. Running Sum of 1d Array
# id: problem1480
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/running-sum-of-1d-array/description/>
# hidden: true
# ---
#
# Given an array `nums`. We define a running sum of an array as `runningSum[i] =
# sum(nums[0]…nums[i])`.
#
# Return the running sum of `nums`.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3,4]
#     Output: [1,3,6,10]
#     Explanation: Running sum is obtained as follows: [1, 1+2, 1+2+3, 1+2+3+4].
#
# **Example 2:**
#
#
#
#     Input: nums = [1,1,1,1,1]
#     Output: [1,2,3,4,5]
#     Explanation: Running sum is obtained as follows: [1, 1+1, 1+1+1, 1+1+1+1, 1+1+1+1+1].
#
# **Example 3:**
#
#
#
#     Input: nums = [3,1,2,10,1]
#     Output: [3,4,6,16,17]
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 1000`
#   * `-10^6 <= nums[i] <= 10^6`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Test
using Documenter
using DocumenterEpub

pages = [
    "index.md",
    "Footnotes" => "footnotes.md",
    "Subpage" => "sub/subpage.md",
    "Autodocs" => "autodocs.md",
    hide("Hidden section" => "hiddenpage.md"),
]

module EPUBTest
# a complicated docstring
"""
    QRCompactWY <: Factorization

A QR matrix factorization stored in a compact blocked format, typically obtained from
[`qr`]. If ``A`` is an `m`×`n` matrix, then

```math
A = Q R
```

where ``Q`` is an orthogonal/unitary matrix and ``R`` is upper triangular. It is similar
to the [`QR`] format except that the orthogonal/unitary matrix ``Q`` is stored in
*Compact WY* format [^Schreiber1989].  For the block size ``n_b``, it is stored as
a `m`×`n` lower trapezoidal matrix ``V`` and a matrix ``T = (T_1 \\; T_2 \\; ... \\;
T_{b-1} \\; T_b')`` composed of ``b = \\lceil \\min(m,n) / n_b \\rceil`` upper triangular
matrices ``T_j`` of size ``n_b``×``n_b`` (``j = 1, ..., b-1``) and an upper trapezoidal
``n_b``×``\\min(m,n) - (b-1) n_b`` matrix ``T_b'`` (``j=b``) whose upper square part
denoted with ``T_b`` satisfying

```math
Q = \\prod_{i=1}^{\\min(m,n)} (I - \\tau_i v_i v_i^T)
= \\prod_{j=1}^{b} (I - V_j T_j V_j^T)
```

such that ``v_i`` is the ``i``th column of ``V``, ``\\tau_i`` is the ``i``th element
of `[diag(T_1); diag(T_2); …; diag(T_b)]`, and ``(V_1 \\; V_2 \\; ... \\; V_b)``
is the left `m`×`min(m, n)` block of ``V``.  When constructed using [`qr`],
the block size is given by ``n_b = \\min(m, n, 36)``.

Iterating the decomposition produces the components `Q` and `R`.

The object has two fields:

* `factors`, as in the [`QR`] type, is an `m`×`n` matrix.

  - The upper triangular part contains the elements of ``R``, that is `R =
    triu(F.factors)` for a `QR` object `F`.

  - The subdiagonal part contains the reflectors ``v_i`` stored in a packed format such
    that `V = I + tril(F.factors, -1)`.

* `T` is a ``n_b``-by-``\\min(m,n)`` matrix as described above. The subdiagonal elements
  for each triangular matrix ``T_j`` are ignored.

!!! note

    This format should not to be confused with the older *WY* representation
    [^Bischof1987].

[^Bischof1987]: C Bischof and C Van Loan, "The WY representation for products of Householder matrices", SIAM J Sci Stat Comput 8 (1987), s2-s13. [doi:10.1137/0908009](https://doi.org/10.1137/0908009)

[^Schreiber1989]: R Schreiber and C Van Loan, "A storage-efficient WY representation for products of Householder transformations", SIAM J Sci Stat Comput 10 (1989), 53-57. [doi:10.1137/0910005](https://doi.org/10.1137/0910005)
"""
complicated(x) = x

"""
    simple(x)
A very simple docstring. [`complicated`](@ref)
"""
simple(x) = x

"""
    multidoc(x)
Invokation with one param.
"""
multidoc(x) = x

"""
    multidoc(x,y)
Invokation with two params.
"""
multidoc(x,y) = x + y

"""
    multidoc(x,y,z)
Invokation with three params.
"""
multidoc(x,y,z) = x + y + z

export complicated, simple
end

isepubcheck = false
try
    if Sys.iswindows()
        run(`where epubcheck`)
    else
        run(`which epubcheck`)
    end
global isepubcheck = true
catch e
    @warn "Couldn't find `epubcheck` - necessary for running the test suite"
end

testconfigs = [
    (;sitename="color",color=true,pages=pages,version="",authors="",lang="en"),
    (;sitename="gray",color=false,pages=pages,version="",authors="MemeLord3000",lang="en"),
    (;sitename="nopages",color=false,pages=[],version="3.0",authors="",lang="de")
]

# the test relies on the external tool "epubcheck" to see if there are errors...
@testset "epubcheck $(config[:sitename])" for config in testconfigs

    makedocs(;
        sitename=config[:sitename],
        modules=[EPUBTest],
        version=config[:version],
        authors=config[:authors],
        source=abspath(joinpath(@__DIR__, "doc")),
        pages=config[:pages],
        format=EPUB(color=config[:color], lang=config[:lang]),
    )

    expectedfile = joinpath(@__DIR__, "build", config[:sitename] * ".epub")
    @test isfile(expectedfile)
    if isepubcheck
        cmd = `epubcheck -q $expectedfile`
        err = Pipe()
        proc = run(pipeline(ignorestatus(cmd), stdout=devnull, stderr=err))
        close(err.in)
        errmsg = String(read(err))
        # we expect one error for apple metadata and two additional output lines
        @test count(==('\n'), errmsg) == 3
    end
end