text:
const DOCUMENTER_DEP = PackageSpec(;
    name="Documenter",
    uuid="e30172f5-a6a5-5a46-863b-614d45cd2de4",
)

struct NoDeploy end
const DeployStyle = Union{TravisCI, GitHubActions, GitLabCI, NoDeploy}
const YesDeploy = Union{TravisCI, GitHubActions, GitLabCI}
const GitHubPagesStyle = Union{TravisCI, GitHubActions}

"""
    Logo(; light=nothing, dark=nothing)

Logo information for documentation.

## Keyword Arguments
- `light::AbstractString`: Path to a logo file for the light (default) theme.
- `dark::AbstractString`: Path to a logo file for the dark theme.
"""
@with_kw_noshow struct Logo
    light::Union{String, Nothing} = nothing
    dark::Union{String, Nothing} = nothing
end

"""
    Documenter{T<:Union{TravisCI, GitLabCI, GitHubActions, NoDeploy}}(;
        make_jl="$(contractuser(default_file("docs", "make.jl")))",
        index_md="$(contractuser(default_file("docs", "src", "index.md")))",
        assets=String[],
        logo=Logo(),
        canonical_url=make_canonical(T),
        makedocs_kwargs=Dict{Symbol, Any}(),
    )

Sets up documentation generation via [Documenter.jl](https://github.com/JuliaDocs/Documenter.jl).
Documentation deployment depends on `T`, where `T` is some supported CI plugin,
or `Nothing` to only support local documentation builds.

## Supported Type Parameters
- `GitHubActions`: Deploys documentation to [GitHub Pages](https://pages.github.com)
  with the help of [`GitHubActions`](@ref).
- `TravisCI`: Deploys documentation to [GitHub Pages](https://pages.github.com)
  with the help of [`TravisCI`](@ref).
- `GitLabCI`: Deploys documentation to [GitLab Pages](https://pages.gitlab.com)
  with the help of [`GitLabCI`](@ref).
- `NoDeploy` (default): Does not set up documentation deployment.

## Keyword Arguments
- `make_jl::AbstractString`: Template file for `make.jl`.
- `index_md::AbstractString`: Template file for `index.md`.
- `assets::Vector{<:AbstractString}`: Extra assets for the generated site.
- `logo::Logo`: A [`Logo`](@ref) containing documentation logo information.
- `canonical_url::Union{Function, Nothing}`: A function to generate the site's canonical URL.
  The default value will compute GitHub Pages and GitLab Pages URLs
  for [`TravisCI`](@ref) and [`GitLabCI`](@ref), respectively.
  If set to `nothing`, no canonical URL is set.
- `makedocs_kwargs::Dict{Symbol}`: Extra keyword arguments to be inserted into `makedocs`.
- `devbranch::Union{AbstractString, Nothing}`: Branch that will trigger docs deployment.

!!! note
    If deploying documentation with Travis CI, don't forget to complete
    [the required configuration](https://juliadocs.github.io/Documenter.jl/stable/man/hosting/#SSH-Deploy-Keys-1).
"""
struct Documenter{T<:DeployStyle} <: Plugin
    assets::Vector{String}
    logo::Logo
    makedocs_kwargs::Dict{Symbol}
    canonical_url::Union{Function, Nothing}
    make_jl::String
    index_md::String
    devbranch::Union{String, Nothing}
end

# Can't use @plugin because we're implementing our own no-arguments constructor.
function Documenter{T}(;
    assets::Vector{<:AbstractString}=String[],
    logo::Logo=Logo(),
    makedocs_kwargs::Dict{Symbol}=Dict{Symbol, Any}(),
    canonical_url::Union{Function, Nothing}=make_canonical(T),
    make_jl::AbstractString=default_file("docs", "make.jl"),
    index_md::AbstractString=default_file("docs", "src", "index.md"),
    devbranch::Union{String, Nothing}=nothing,
) where T <: DeployStyle
    return Documenter{T}(
        assets,
        logo,
        makedocs_kwargs,
        canonical_url,
        make_jl,
        index_md,
        devbranch,
    )
end

Documenter(; kwargs...) = Documenter{NoDeploy}(; kwargs...)

# We have to define these manually because we didn't use @plugin.
defaultkw(::Type{<:Documenter}, ::Val{:assets}) = String[]
defaultkw(::Type{<:Documenter}, ::Val{:logo}) = Logo()
defaultkw(::Type{<:Documenter}, ::Val{:make_jl}) = default_file("docs", "make.jl")
defaultkw(::Type{<:Documenter}, ::Val{:index_md}) = default_file("docs", "src", "index.md")
defaultkw(::Type{<:Documenter}, ::Val{:devbranch}) = nothing

gitignore(::Documenter) = ["/docs/build/"]
priority(::Documenter, ::Function) = DEFAULT_PRIORITY - 1  # We need SrcDir to go first.

badges(::Documenter) = Badge[]
badges(::Documenter{<:GitHubPagesStyle}) = [
    Badge(
        "Stable",
        "https://img.shields.io/badge/docs-stable-blue.svg",
        "https://{{{USER}}}.github.io/{{{PKG}}}.jl/stable",
    ),
    Badge(
        "Dev",
        "https://img.shields.io/badge/docs-dev-blue.svg",
        "https://{{{USER}}}.github.io/{{{PKG}}}.jl/dev",
    ),
]
badges(::Documenter{GitLabCI}) = Badge(
    "Dev",
    "https://img.shields.io/badge/docs-dev-blue.svg",
    # TODO: Support custom domain here.
    "https://{{{USER}}}.gitlab.io/{{{PKG}}}.jl/dev",
)

view(p::Documenter, t::Template, pkg::AbstractString) = Dict(
    "ASSETS" => map(basename, p.assets),
    "AUTHORS" => join(t.authors, ", "),
    "CANONICAL" => p.canonical_url === nothing ? nothing : p.canonical_url(t, pkg),
    "HAS_ASSETS" => !isempty(p.assets),
    "MAKEDOCS_KWARGS" => map(((k, v),) -> k => repr(v), sort(collect(p.makedocs_kwargs))),
    "PKG" => pkg,
    "REPO" => "$(t.host)/$(t.user)/$pkg.jl",
    "USER" => t.user,
    "BRANCH" => p.devbranch === nothing ? default_branch(t) : p.devbranch,
)

function view(p::Documenter{<:GitHubPagesStyle}, t::Template, pkg::AbstractString)
    base = invoke(view, Tuple{Documenter, Template, AbstractString}, p, t, pkg)
    return merge(base, Dict("HAS_DEPLOY" => true))
end

function validate(p::Documenter, ::Template)
    foreach(p.assets) do a
        isfile(a) || throw(ArgumentError("Asset file $a does not exist"))
    end
    foreach((:light, :dark)) do k
        logo = getfield(p.logo, k)
        if logo !== nothing && !isfile(logo)
            throw(ArgumentError("Logo file $logo does not exist"))
        end
    end
end

function validate(p::Documenter{T}, t::Template) where T <: YesDeploy
    invoke(validate, Tuple{Documenter, Template}, p, t)
    if !hasplugin(t, T)
        name = nameof(T)
        s = "Documenter: The $name plugin must be included for docs deployment to be set up"
        throw(ArgumentError(s))
    end
end

function hook(p::Documenter, t::Template, pkg_dir::AbstractString)
    pkg = basename(pkg_dir)
    docs_dir = joinpath(pkg_dir, "docs")

    # Generate files.
    make = render_file(p.make_jl, combined_view(p, t, pkg), tags(p))
    index = render_file(p.index_md, combined_view(p, t, pkg), tags(p))
    gen_file(joinpath(docs_dir, "make.jl"), make)
    gen_file(joinpath(docs_dir, "src", "index.md"), index)

    # Copy over any assets.
    assets_dir = joinpath(docs_dir, "src", "assets")
    mkpath(assets_dir)
    foreach(a -> cp(a, joinpath(assets_dir, basename(a))), p.assets)
    foreach((:light => "logo", :dark => "logo-dark")) do (k, f)
        logo = getfield(p.logo, k)
        if logo !== nothing
            _, ext = splitext(logo)
            cp(logo, joinpath(assets_dir, "$f$ext"))
        end
    end
    isempty(readdir(assets_dir)) && rm(assets_dir)

    # Create the documentation project.
    with_project(docs_dir) do
        Pkg.add(DOCUMENTER_DEP)
        cd(() -> Pkg.develop(PackageSpec(; path="..")), docs_dir)
    end
end

github_pages_url(t::Template, pkg::AbstractString) = "https://$(t.user).github.io/$pkg.jl"
gitlab_pages_url(t::Template, pkg::AbstractString) = "https://$(t.user).gitlab.io/$pkg.jl"

make_canonical(::Type{<:GitHubPagesStyle}) = github_pages_url
make_canonical(::Type{GitLabCI}) = gitlab_pages_url
make_canonical(::Type{NoDeploy}) = nothing

needs_username(::Documenter) = true

function customizable(::Type{<:Documenter})
    return (:canonical_url => NotCustomizable, :makedocs_kwargs => NotCustomizable)
end

function interactive(::Type{Documenter})
    styles = [NoDeploy, TravisCI, GitLabCI, GitHubActions]
    menu = RadioMenu(map(string, styles); pagesize=length(styles))
    println("Documenter deploy style:")
    idx = request(menu)
    return interactive(Documenter{styles[idx]})
end

function prompt(::Type{<:Documenter}, ::Type{Logo}, ::Val{:logo})
    light = Base.prompt("Enter value for 'logo.light' (String, default=nothing)")
    dark = Base.prompt("Enter value for 'logo.dark' (String, default=nothing)")
    return Logo(; light=light, dark=dark)
end


text:
# p9.jl - polynomial interpolation in equispaced and Chebyshev pts

N = 16;
xx = -1.01:.005:1.01; clf();
for i = 1:2
    i==1 && ( (s,x) = ("equispaced points", -1 + 2*(0:N)/N) );
    i==2 && ( (s,x) = ("Chebyshev points", cos.(pi*(0:N)/N)) );
    subplot(2,2,i)
    u = 1./(1+16*x.^2);
    uu = 1./(1+16*xx.^2);
    p = polyfit(x,u);              # interpolation
    pp = p(xx);                    # evaluation of interpolant
    plot(x,u,".",markersize=6)
    plot(xx,pp)
    axis([-1.1,1.1,-1,1.5]); title(s);
    error = signif( norm(uu-pp,Inf), 5);
    text(-.5,-.5,"max error = $error",fontsize=8);
end

ext: jl


text:
println("Hello world!");

@everywhere using Printf;

@show workdir = if length(ARGS) > 0
  ARGS[1];
else
  ".";
end

@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x));
@everywhere fmt_int(x) = @sprintf("%03d", x);

@everywhere function prefix(case)
  "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))";
end

cases = Any[];
E0s = [10.0; 1.0; 1e-3];
Ks = zip([1.0; 0.0; 0.5; 1.0], [0.0; 1.0; 1.0; 0.5])
kTs = [100.0; 1.0; 1e-3];
Fzs = [0.0; 1e-3; 1.0; 10.0];
Fxs = [0.0; 1e-3; 1.0; 10.0];
ns = Int[10; 100; 250; 1000];
bs = [1.0; 0.1; 10.0];
for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, Kpair in Ks,
    E0 in E0s
    push!(cases, Dict(:E0 => E0, :K1 => Kpair[1], :K2 => Kpair[2],
                      :kT => kT, :Fz => Fz,
                      :Fx => Fx, :n => n, :b => b));
end

@info "total number of cases to run: $(length(cases))";

pmap(case -> begin;
  for iter=1:10
    outfile = joinpath(workdir, "$(prefix(case))_iter-$iter.out");
    if !isfile(outfile)
      println("Running case: $case.");
      command = `julia -O 3 interacting_chain.jl -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 100000 -v 2 --prefix $(joinpath(workdir, prefix(case)*"_iter-$iter"))`;
      output = read(command, String);
      write(outfile, output);
    else
      println("Case: $case, $iter has already been run.");
    end
  end
end, cases);

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats


text:
"""
```
plotNodes(network::DislocationNetwork, args...; kw...)
```
Plots [`DislocationNetwork`](@ref).
"""
function plotNodes(network::DislocationNetwork, args...; kw...)
    fig = plot()
    plotNodes!(fig, network, args...; kw...)
end
"""
```
plotNodes!(fig, network::DislocationNetwork, args...; kw...)
```
In-place plots [`DislocationNetwork`](@ref).
"""
function plotNodes!(fig, network::DislocationNetwork, args...; kw...)
    label = network.label
    coord = network.coord
    links = network.links
    numNode = network.numNode[1]
    elemT = eltype(links)
    @inbounds for i in 1:numNode
        n1n2 = SVector{2, elemT}(links[1, i], links[2, i])
        label[n1n2[1]] == extDln || label[n1n2[2]] == extDln && continue
        plot!(fig, coord[1, n1n2], coord[2, n1n2], coord[3, n1n2], args...; kw...)
        #=
        # quiver needs to be implemented in Plots.jl but we can use python.
        lVec = coord[n2, :] - coord[n1, :]
        quiver!([coord[n1,1]], [coord[n1,2]], [coord[n1,3]], args...; quiver=([lVec[1]], [lVec[2]], [lVec[3]]), kw...) =#
    end
    plot!(fig; xlabel = "x", ylabel = "y", zlabel = "z")
    return fig
end

"""
```
plotNodes(mesh::AbstractMesh, network::DislocationNetwork, args...; kw...)
```
Plots [`DislocationNetwork`](@ref) inside [`AbstractMesh`](@ref).
"""
function plotNodes(mesh::AbstractMesh, network::DislocationNetwork, args...; kw...)
    fig = plot()
    plotNodes!(fig, mesh, network, args...; kw...)
end

"""
```
plotNodes!(fig, mesh::AbstractMesh, network::DislocationNetwork, args...; kw...)
```
In-place plots [`DislocationNetwork`](@ref) inside [`AbstractMesh`](@ref).
"""
function plotNodes!(
    fig,
    mesh::T1,
    network::DislocationNetwork,
    args...;
    kw...,
) where {T1 <: AbstractMesh}
    label = network.label
    coord = network.coord
    links = network.links
    numNode = network.numNode[1]
    elemT = eltype(links)
    @inbounds for i in 1:numNode
        n1n2 = SVector{2, elemT}(links[1, i], links[2, i])
        label[n1n2[1]] == extDln || label[n1n2[2]] == extDln ? continue : nothing
        plot!(fig, coord[1, n1n2], coord[2, n1n2], coord[3, n1n2], args...; kw...)
        # quiver needs to be implemented in Plots.jl but we can use python.
        #=
        lVec = coord[n2, :] - coord[n1, :]
        quiver!([coord[n1,1]], [coord[n1,2]], [coord[n1,3]], args...; quiver=([lVec[1]], [lVec[2]], [lVec[3]]), kw...) =#
    end

    # Plot FEM domain.
    vertices = reshape(collect(Iterators.flatten(mesh.vertices.vertices)), 3, :)
    vertices = reshape(collect(Iterators.flatten(mesh.vertices.vertices)), 3, :)
    face1 = [mesh.faces[:, 5]; mesh.faces[1, 5]]
    face2 = [mesh.faces[:, 6]; mesh.faces[1, 6]]
    surf1 = vertices[:, face1]
    surf2 = vertices[:, face2]
    plot!(fig, surf1[1, :], surf1[2, :], surf1[3, :], linecolor = :black, linewidth = 2)
    plot!(fig, surf2[1, :], surf2[2, :], surf2[3, :], linecolor = :black, linewidth = 2)
    side = vertices[:, [1, 5]]
    plot!(fig, side[1, :], side[2, :], side[3, :], linecolor = :black, linewidth = 2)
    side = vertices[:, [2, 6]]
    plot!(fig, side[1, :], side[2, :], side[3, :], linecolor = :black, linewidth = 2)
    side = vertices[:, [3, 7]]
    plot!(fig, side[1, :], side[2, :], side[3, :], linecolor = :black, linewidth = 2)
    side = vertices[:, [4, 8]]
    plot!(fig, side[1, :], side[2, :], side[3, :], linecolor = :black, linewidth = 2)

    xlims = (
        minimum(vertices[1, :]) - 0.2 * maximum(vertices[1, :]),
        maximum(vertices[1, :]) + 0.2 * maximum(vertices[1, :]),
    )
    ylims = (
        minimum(vertices[2, :]) - 0.2 * maximum(vertices[2, :]),
        maximum(vertices[2, :]) + 0.2 * maximum(vertices[2, :]),
    )
    zlims = (
        minimum(vertices[3, :]) - 0.2 * maximum(vertices[3, :]),
        maximum(vertices[3, :]) + 0.2 * maximum(vertices[1, :]),
    )
    plot!(fig, xlims = xlims, ylims = ylims, zlims = zlims)
    plot!(fig; xlabel = "x", ylabel = "y", zlabel = "z")

    return fig
end

"""
```
plotNodes(loop::DislocationLoop, args...; kw...)
```
Plots [`DislocationLoop`](@ref).
"""
function plotNodes(loop::DislocationLoop, args...; kw...)
    fig = plot()
    plotNodes!(fig, loop::DislocationLoop, args...; kw...)
end
"""
```
plotNodes!(fig, loop::DislocationLoop, args...; kw...)
```
In-place plots [`DislocationLoop`](@ref).
"""
function plotNodes!(fig, loop::DislocationLoop, args...; kw...)
    idx = findall(x -> x != 0, loop.label)
    coord = loop.coord
    links = loop.links
    elemT = eltype(links)
    @inbounds for i in idx
        n1n2 = SVector{2, elemT}(links[1, i], links[2, i])
        plot!(fig, coord[1, n1n2], coord[2, n1n2], coord[3, n1n2], args...; kw...)
        #=
        # quiver needs to be implemented in Plots.jl but we can use python.
        lVec = coord[n2, :] - coord[n1, :]
        quiver!([coord[n1,1]], [coord[n1,2]], [coord[n1,3]], args...; quiver=([lVec[1]], [lVec[2]], [lVec[3]]), kw...) =#
    end
    plot!(fig; xlabel = "x", ylabel = "y", zlabel = "z")
    return fig
end

"""
```
plotFEDomain(mesh::AbstractMesh)
```
Plots corners, edges and surfaces of [`AbstractMesh`](@ref).
"""
function plotFEDomain(mesh::AbstractMesh, args...; kw...)
    surfNode = mesh.surfNode
    cornerNode = mesh.cornerNode
    edgeNode = mesh.edgeNode
    faceNode = mesh.faceNode
    coord = mesh.coord

    cornerIdx =
        collect(Iterators.flatten([surfNode[cornerNode[i]] for i in 1:length(cornerNode)]))
    fig = scatter(
        coord[1, cornerIdx],
        coord[2, cornerIdx],
        coord[3, cornerIdx],
        markershape = :diamond,
        markersize = 3,
        label = "Corners",
        args...;
        kw...,
    )

    @inbounds for i in 1:length(edgeNode)
        scatter!(
            fig,
            coord[1, surfNode[edgeNode[i]]],
            coord[2, surfNode[edgeNode[i]]],
            coord[3, surfNode[edgeNode[i]]],
            markershape = :circle,
            markersize = 3,
            label = "Edge $i",
            args...;
            kw...,
        )
    end
    @inbounds for i in 1:length(faceNode)
        scatter!(
            fig,
            coord[1, surfNode[faceNode[i]]],
            coord[2, surfNode[faceNode[i]]],
            coord[3, surfNode[faceNode[i]]],
            markershape = :square,
            markersize = 3,
            label = "Face $i",
            args...;
            kw...,
        )
    end
    plot!(fig; xlabel = "x", ylabel = "y", zlabel = "z")
    return fig
end

function plotBoundaries(boundaries::Boundaries, mesh::RegularCuboidMesh, args...; kw...)
    uGamma = boundaries.uGamma.node
    mGamma = boundaries.mGamma.node
    tGamma = boundaries.tGamma.node
    feCoord = mesh.coord

    fig = scatter(
        feCoord[1, uGamma],
        feCoord[2, uGamma],
        feCoord[3, uGamma];
        markershape = :square,
        markercolor = :red,
        markersize = 3,
        label = "Fixed",
    )
    scatter!(
        fig,
        feCoord[1, mGamma],
        feCoord[2, mGamma],
        feCoord[3, mGamma];
        markershape = :diamond,
        markercolor = :cyan,
        markersize = 3,
        label = "Load",
    )
    scatter!(
        fig,
        feCoord[1, tGamma],
        feCoord[2, tGamma],
        feCoord[3, tGamma];
        markershape = :circle,
        markercolor = :black,
        markersize = 3,
        label = "Traction",
    )
    plot!(fig, args...; xlabel = "x", ylabel = "y", zlabel = "z", kw...)
    return fig
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


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

mutable struct LoadBalancerFrontendIPConfigurationListResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{FrontendIPConfiguration} } # spec name: value
    nextLink::Any # spec type: Union{ Nothing, String } # spec name: nextLink

    function LoadBalancerFrontendIPConfigurationListResult(;value=nothing, nextLink=nothing)
        o = new()
        validate_property(LoadBalancerFrontendIPConfigurationListResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(LoadBalancerFrontendIPConfigurationListResult, Symbol("nextLink"), nextLink)
        setfield!(o, Symbol("nextLink"), nextLink)
        o
    end
end # type LoadBalancerFrontendIPConfigurationListResult

const _property_map_LoadBalancerFrontendIPConfigurationListResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"), Symbol("nextLink")=>Symbol("nextLink"))
const _property_types_LoadBalancerFrontendIPConfigurationListResult = Dict{Symbol,String}(Symbol("value")=>"Vector{FrontendIPConfiguration}", Symbol("nextLink")=>"String")
Base.propertynames(::Type{ LoadBalancerFrontendIPConfigurationListResult }) = collect(keys(_property_map_LoadBalancerFrontendIPConfigurationListResult))
Swagger.property_type(::Type{ LoadBalancerFrontendIPConfigurationListResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_LoadBalancerFrontendIPConfigurationListResult[name]))}
Swagger.field_name(::Type{ LoadBalancerFrontendIPConfigurationListResult }, property_name::Symbol) =  _property_map_LoadBalancerFrontendIPConfigurationListResult[property_name]

function check_required(o::LoadBalancerFrontendIPConfigurationListResult)
    true
end

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


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dcc_interval

"""
    dcc_interval(;kwargs...)

An Interval component.
A component that repeatedly increments a counter `n_intervals`
with a fixed time delay between each increment.
Interval is good for triggering a component on a recurring basis.
The time delay is set with the property "interval" in milliseconds.
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `interval` (Real; optional): This component will increment the counter `n_intervals` every
`interval` milliseconds
- `disabled` (Bool; optional): If True, the counter will no longer update
- `n_intervals` (Real; optional): Number of times the interval has passed
- `max_intervals` (Real; optional): Number of times the interval will be fired.
If -1, then the interval has no limit (the default)
and if 0 then the interval stops running.
"""
function dcc_interval(; kwargs...)
        available_props = Symbol[:id, :interval, :disabled, :n_intervals, :max_intervals]
        wild_props = Symbol[]
        return Component("dcc_interval", "Interval", "dash_core_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: asc-csa/dash-core-components

max_issues_repo_name: asc-csa/dash-core-components

max_forks_repo_name: asc-csa/dash-core-components


text:
# linearfitxy.jl

struct stfitxy
    X   :: Vector{Float64}
    Y   :: Vector{Float64}
    σX  :: Union{Float64, Vector{Float64}}
    σY  :: Union{Float64, Vector{Float64}}
    r   :: Union{Float64, Vector{Float64}}
    a   :: Float64
    b   :: Float64
    σa  :: Float64
    σb  :: Float64
    σa95:: Float64
    σb95:: Float64
    S   :: Float64
    ρ   :: Float64
    bᵢ  :: Vector{Float64}
    i   :: Int64
  end

"""
linearfitxy(X, Y; σX=0, σY=0, r=0, isplot=false, ratio=1)

Performs 1D linear fitting of experimental data with uncertainties in  X and Y:
- Linear fit:             `Y = a + b*X`                               [1]
- Errors:                 ``X ± σX;  Y ± σY``                         [2]
- Errors' correlation:    ``r =  = cov(σX, σY) / (σX * σY)``          [3]

# Arguments:
- `X` and `Y` are input data vectors with length ≥ 3
- Optional standard deviation errors ``σX`` and ``σY`` are vectors or scalars
- Optional `r` is the correlation between the ``σX`` and ``σY`` errors
`r` can be a vector or scalar

``σX`` and ``σY`` errors (error ellipses) with bivariate Gaussian distribution assumed.
If no errors, or if only ``σX`` or ``σY`` are provided, then the results are equivalent
to those from the [LsqFit.jl](https://github.com/JuliaNLSolvers/LsqFit.jl) package.

Based on York et al. (2004) with extensions (confidence intervals, diluted corr. coeff.).

# Examples:
```julia-repl
st = linearfitxy(X, Y)    # no errors in X and Y, no plot displayed

st = linearfitxy(X, Y; σX, σY, isplot=true) # XY errors not correlated (r=0); plot ratio=1

linearfitxy(X, Y; σX,σY,r=0,isplot=true,ratio=:auto) # XY errors not correlated (r=0); plot
```

The results are in the fields of the returned st::stfitxy structure:
- The intercept `a`, the slope `b` and their uncertainties `σa` and `σb`
- ``σa95`` and ``σb95``: 95%-confidence interval using two-tailed t-Student distribution,
    e.g.: ``b ± σb95 = b ± t(0.975,N-2)*σb``
- Goodness of fit `S` (reduced ``Χ²`` test): quantity with ``Χ²`` N-2 degrees of freedom
  `S ~ 1`: fit consistent with errors, `S > 1`: poor fit, `S >> 1`: errors underestimated,
  `S < 1`: overfitting or errors overestimated
- Pearson's correlation coefficient ``ρ`` that accounts for data errors
- Optional display of fit results with error ellipses and confidence intervals

The default argument `isplot=false` can be turned on to plot the results.
Currently Plots.jl's gr() is used

"""
function linearfitxy(X, Y; σX=0, σY=0, r=0, isplot=false, ratio=1)

    N =   length(X)
    d = TDist(N-2)     # t-Student distribution with N-2 degrees of freedom
    cf = quantile(d, 0.975)  # correction factor for 95% confidence intervals (two-tailed distribution)

    if (σX == 0) && (σY == 0)
        σX = σY = r = zeros(N)

        M = [ones(N) X]
        (a, b) = M \ Y
        X̄ = mean(X)    # Ȳ = mean(Y)

        # standard error * Ŝ factor defined in Cantrell (2008)
        Ŝ = sqrt(sum((Y - b*X .- a).^2)/(N-2))  # goodness of fit
        σb = sqrt(1/sum((X .- X̄ ).^2))
        σa = Ŝ * sqrt(1/N + X̄^2 * σb^2)
        σb *= Ŝ

        # Pearson's correlation coefficient:
        ρ = cov(X,Y)/sqrt(var(X) * var(Y))

        st = stfitxy(X, Y, σX, σY, r, a, b, σa, σb, cf*σa, cf*σb, Ŝ, ρ, [b], 1)

    else
        Nmax = 50;                  # maximum number of iterations
        tol = 1e-15;                #relative tolerance to stop at

        σX == 0 && (σX = 1e-16)
        σY == 0 && (σY = 1e-16)
        length(σX) == 1 &&  (σX = σX*ones(N))
        length(σY) == 1 &&  (σY = σY*ones(N))
        length(r) == 1  &&  (r = r*ones(N))

        M = [ones(N) X]
        (a, b) = M \ Y              # initial guess for b
        bᵢ = zeros(Nmax+1);         #vector to save b iterations in
        bᵢ[1] = b
        X̄ = 0; Ȳ = 0;
        β = zeros(N)
        W = zeros(N)
        local i
        for outer i = 1:Nmax        # KEY TRICK
            ωX = 1 ./σX.^2
            ωY = 1 ./σY.^2
            α = sqrt.(ωX .* ωY)
            W .= ωX .* ωY ./(ωX + b^2 * ωY - 2*b * r .* α)
            X̄ = sum(W .* X)/sum(W)
            Ȳ = sum(W .* Y)/sum(W)
            U = X .- X̄
            V = Y .- Ȳ
            β .= W .*(U ./ ωY + b * V ./ ωX - (b * U + V) .* r ./ α)
            b = sum(W .* β .* V)/sum(W .* β .* U)
            bᵢ[i+1] = b
            abs((bᵢ[i+1] - bᵢ[i])/bᵢ[i+1]) < tol && break
        end
        a = Ȳ - b*X̄
        x = X̄ .+ β                  # y = Ȳ + b*β
        X̄ = sum(W .* x)/sum(W)      # Ȳ = sum(W .* y)/sum(W)

        # compare to: Ŝ = sqrt(sum((Y - b*X .- a).^2)/(N-2))
        Ŝ = sqrt(sum(W .* (Y - b*X .- a).^2) /(N-2))  # goodness of fit (York, 2004), (Cantrell, 2008)

        σb = sqrt(1/sum(W .* (x .- X̄).^2))
        σa = Ŝ * sqrt(1/sum(W) + X̄^2 * σb^2)
        σb *= Ŝ             # standard error * Ŝ factor (Cantrell, 2008)

        # See Wikipedia Regression dilution (Pearson's correlation coefficient with errors in variables)
        vX = var(X); vY = var(Y)
        ρ = cov(X,Y)/sqrt(vX*vY) * sqrt(vX/(vX + var(σX))) * sqrt(vY/(vY + var(σY)))

        st = stfitxy(X, Y, σX, σY, r, a, b, σa, σb, cf*σa, cf*σb, Ŝ, ρ, bᵢ, i)
    end

    if isplot
        plot_linfitxy(st; ratio=ratio)
    end

    @printf("\n>>> [± σ]  Y = (%.4f ± %.4f) + (%.4f ± %.4f)*X \n", a,  σa, b, σb)
    @printf(">>> [± 95%% CI]  Y = (%.4f ± %.4f) + (%.4f ± %.4f)*X \n", a, cf*σa, b, cf*σb)
    @printf(">>> Pearson ρ = %.3f;  Goodness of fit = %.3f \n\n", ρ, Ŝ)

    return st
end

function plot_covariance_ellipses!(X, Y, a,  b, c; lc=:lightgrey, lw=0.3)
    # https://cookierobotics.com/007/
    # covariance matrix = [a b; b c]
    t = LinRange(0, 2π, 72)
    st = sin.(t)
    ct = cos.(t)
    for (X, Y, a, b, c) in zip(X, Y, a, b, c)
        r = sqrt(((a - c)/2)^2 + b^2)
        λ₁ = sqrt(abs((a + c)/2 + r))
        λ₂ = sqrt(abs((a + c)/2 - r))

        (b == 0 && a >= c) ? (θ = 0.) :
        (b == 0 && a < c)  ? (θ = π/2) :
        (θ = atan(λ₁^2 - a, b))

        xₜ = @. X + λ₁*cos(θ)*ct - λ₂*sin(θ)*st
        yₜ = @. Y + λ₁*sin(θ)*ct + λ₂*cos(θ)*st
        plot!(xₜ, yₜ, lc=lc, lw=0.3, label=false)
    end
    display(plot!())
end

function  plot_linfitxy(st::stfitxy; ratio=1)
# set plotting defaults
    plot_font = "Computer Modern";
    default(fontfamily=plot_font,framestyle=:axes,yminorgrid=true, legendtitlefontsize=6,fg_color_legend=nothing,
            legendfontsize=6, guidefont=(7,:black),tickfont=(6,:black),size=(600,400), dpi=300,
            titlefont = (6, plot_font))

    X = st.X;  Y = st.Y;  σX = st.σX;  σY = st.σY;  r = st.r;
    covXY = σX .* σY .* r
    # sort data because of ribbon:
    I = sortperm(X);
    X = X[I];  Y = Y[I]; σX = σX[I];  σY = σY[I]; covXY = covXY[I];

    str1 = @sprintf("[± σ]  Y = (%.4f ± %.4f) + (%.4f ± %.4f)*X", st.a, st.σa, st.b, st.σb)
    str2 = @sprintf("\n[± 95%% CI]  Y = (%.4f ± %.4f) + (%.4f ± %.4f)*X", st.a, st.σa95, st.b, st.σb95)
    str3 = @sprintf("\nPearson ρ = %.2f; Goodness of fit = %.2f", st.ρ, st.S)

    dX = diff([extrema(X)...])[1]/7    # extends plot x-axis by 1/7 each side
    x1, x2 = (-dX, dX) .+ extrema(X)
    xx = [x1; X; x2]
    tl, bl = (st.a - st.σa) .+ (st.b + st.σb)*xx,  (st.a + st.σa) .+ (st.b - st.σb)*xx
    σp, σm = maximum([tl bl], dims=2) .-  (st.a .+ st.b*xx),  (st.a .+ st.b*xx) .- minimum([tl bl], dims=2)

    #TODO: 95%-confidence interval computation for X-Y errors
    N = length(X)
    d = TDist(N-2)     # t-Student distribution with N-2 degrees of freedom
    cf = quantile(d, 0.975)  # correction factor for 95% confidence intervals (two-tailed distribution)
    σx95 = cf * st.S * sqrt.(1/N .+ (X .- mean(X)).^2 / var(X) /(N-1))

    p = plot(xlims=(x1,x2), title=str1*str2*str3, ratio=ratio, legend=:outerbottomright)
    scatter!([x1-dX],[Y[1]], mc=:lightblue, label="Y=(a±σa)+(b±σb)X", marker=:rect)

    #TODO: 95%-confidence interval computation for X-Y errors
    if all(≈(0), σX) && all(≈(0), σY)
        scatter!([x1-dX],[Y[1]], mc=:reds, label="95% confidence interval", marker=:rect)
        plot!(X, st.a .+ st.b*X, color=:reds, ribbon=(σx95,σx95), label=false)
        plot!(xlims = (-dX, dX)./5 .+ extrema(X))
    end
    plot!([x1-dX],[Y[1]], lc=:white, label=" ")   # adds space in legend
    plot!(xx, st.a .+ st.b*xx, color=:lightblue, ribbon=(σp,σm), label=false)
    plot!(xx, st.a .+ st.b*xx, color=:blue, lw=0.5, xlabel="X", ylabel="Y", label="LinearFitXY")

    scatter!(X, Y, msw=0.1, ms=1., msc=:lightgrey, xerror= σX, yerror= σY, label=false)
    scatter!(X, Y, msw=0.1, ms=1.5, mc=:blue, label=false)

    if !all(≈(0), σX) || !all(≈(0), σY)
        plot_covariance_ellipses!(X, Y, σX.^2, covXY, σY.^2; lc=:grey)
    end
    display(p)
end


text:
struct AsymptoticVars{T, U, V}
    f ::T
    fˢ::U
    ζ∞::V
end

function initialize_asymptotics(structure)
    @unpack f, grid = structure
    liquid = f.liquid

    D₀ = diffusion_coeff(liquid)
    K  = nodes(grid)
    S  = project.(f.(K))
    Sˢ = onesof(S)
    d  = dimensionality(liquid)
    w  = weights(Asymptotics, K, S, d, grid)
    υ  = weight(one(D₀), volume_fraction(liquid), d, grid)
    Λ  = lambdaof(liquid).(K)
    B  = K.^2 ./ gett.(Λ)
    Bˢ = nothing

    svars = StaticVars(S, Sˢ, B, Bˢ, w, υ)
    kvars = KernelStaticVars(svars, Λ)

    m  = length(S)
    fˢ = zeros(eltype(Sˢ), m)
    f  = zeros(eltype(S), m)
    ζ∞ = Box(υ * sum(w))

    avars = AsymptoticVars(f, fˢ, ζ∞)

    return avars, kvars, D₀
end


text:
@testset "$TEST $G" begin

g5 = DG(4)
add_edge!(g5,1,2); add_edge!(g5,2,3); add_edge!(g5,1,3); add_edge!(g5,3,4)

@test degree_centrality(g5) == [0.6666666666666666, 0.6666666666666666, 1.0, 0.3333333333333333]
@test in_degree_centrality(g5, normalize=false) == [0.0, 1.0, 2.0, 1.0]
@test out_degree_centrality(g5; normalize=false) == [2.0, 1.0, 1.0, 0.0]

end # testset


text:
import ImageView
using Images, OffsetArrays

# Grayscale
ImageView.imshow(rand(Gray{N0f8}, 10, 10))
# clamping checks
ImageView.imshow(rand(UInt8, 10, 10))  # most/all pixels should be white
A = randn(10, 10)
ImageView.imshow(A)
A[2,2] = NaN
A[3,3] = -Inf
A[4,4] = Inf
ImageView.imshow(A)

# RGB
ImageView.imshow(rand(RGB{N0f8}, 10, 10))
# clamping checks
A = randn(3, 10, 10)
ImageView.imshow(colorview(RGB, A))
A[1,2,2] = NaN
A[1,3,3] = -Inf
A[1,4,4] = Inf
ImageView.imshow(colorview(RGB, A))

# Non-1 indices
A = OffsetArray(rand(11, 10), -5:5, 0:9)
ret = ImageView.imshow(A)


text:
"""
    gaborkernel(;
        sqexponential_transform=IdentityTransform(), cosine_tranform=IdentityTransform()
    )

Construct a Gabor kernel with transformations `sqexponential_transform` and
`cosine_transform` of the inputs of the underlying squared exponential and cosine kernel,
respectively.

# Definition

For inputs ``x, x' \\in \\mathbb{R}^d``, the Gabor kernel with transformations ``f``
and ``g`` of the inputs to the squared exponential and cosine kernel, respectively,
is defined as
```math
k(x, x'; f, g) = \\exp\\bigg(- \\frac{\\| f(x) - f(x')\\|_2^2}{2}\\bigg)
                 \\cos\\big(\\pi \\|g(x) - g(x')\\|_2 \\big).
```
"""
function gaborkernel(;
    sqexponential_transform=IdentityTransform(), cosine_transform=IdentityTransform()
)
    return (SqExponentialKernel() ∘ sqexponential_transform) *
           (CosineKernel() ∘ cosine_transform)
end