text:
using IterativeSolvers
using RecipesBase
using Test

@testset "ConvergenceHistory" begin

    KW = Dict{Symbol, Any}

    RecipesBase.is_key_supported(k::Symbol) = k == :sep ? false : true

    begin
        history = ConvergenceHistory(partial = false)
        history.iters = 3
        history.isconverged = true
        @test string(history) == "Converged after 3 iterations."
        history.isconverged = false
        @test string(history) == "Not converged after 3 iterations."
    end

    # No plottables
    begin
        history = ConvergenceHistory(partial = false)
        @test_throws ErrorException RecipesBase.apply_recipe(KW(), history)
    end

    # Without restart markers
    begin
        history = ConvergenceHistory(partial = false)
        history.iters = 3
        history.data[:resnorm] = [10.0, 3.0, 0.1]

        plots = (
            RecipesBase.apply_recipe(KW(), history),
            RecipesBase.apply_recipe(KW(), history, :resnorm)
        )

        for data in plots
            @test length(data) == 1
            @test data[1].d[:label] == "resnorm"
            @test data[1].d[:seriestype] == :line
        end
    end

    # With restart markers
    begin
        history = ConvergenceHistory(partial = false, restart = 2)
        history.iters = 3
        history.data[:resnorm] = [10.0, 3.0, 0.1]

        plots = (
            RecipesBase.apply_recipe(KW(), history),
            RecipesBase.apply_recipe(KW(), history, :resnorm)
        )

        for data in plots
            @test length(data) == 2
            @test data[2].d[:linecolor] == :white
        end
    end

    # Custom color
    begin
        history = ConvergenceHistory(partial = false, restart = 2)
        history.iters = 3
        history.data[:resnorm] = [10.0, 3.0, 0.1]

        plots = (
            RecipesBase.apply_recipe(KW(:sep => :red), history),
            RecipesBase.apply_recipe(KW(:sep => :red), history, :resnorm)
        )

        for data in plots
            @test length(data) == 2
            @test data[2].d[:linecolor] == :red
        end
    end
end


text:
"""
    isdense(M) -> Bool

Return true if a matrix is dense.

Note:
It is not exactly same as !isparse, e.g. Diagonal, IMatrix and PermMatrix are both not isdense and not isparse.
"""
function isdense end

isdense(M)::Bool = !issparse(M)
isdense(::Diagonal) = false
isdense(x::Transpose) = isdense(parent(x))
isdense(x::Adjoint) = isdense(parent(x))

"""faster invperm"""
function fast_invperm(order)
    v = similar(order)
    @inbounds @simd for i=1:length(order)
        v[order[i]] = i
    end
    v
end

dropzeros!(A::Diagonal; trim::Bool=false) = A

ext: jl


text:
using Smg2s, SparseArrays, LinearAlgebra, Test

dim = 5000
nzeros = 10
nilp = Nilp(10, dim)
#nilp.nilpMat
M=sprandn(dim, dim, 0.25)

@test MNilpM(M, nilp) == M*NilpMat(nilp) - NilpMat(nilp)*M

@time MNilpM(M, nilp)
@time M*NilpMat(nilp) - NilpMat(nilp)*M

@info "done"


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
    ElectricLoad(;
        loads_kw::Union{Missing, Array{<:Real,1}} = missing,
        year::Int = 2020,
        doe_reference_name::Union{Missing, String} = missing,
        blended_doe_reference_names::Array{String, 1} = String[],
        blended_doe_reference_percents::Array{<:Float64,1} = Float64[],
        city::Union{Missing, String} = missing,
        annual_kwh::Union{Real, Nothing} = nothing,
        monthly_totals_kwh::Array{<:Real,1} = Real[],
        critical_loads_kw::Union{Missing, Array{Real,1}} = missing,
        loads_kw_is_net::Bool = true,
        critical_loads_kw_is_net::Bool = false,
        critical_load_pct::Real = 0.5
    )

Must provide either `loads_kw` or [`doe_reference_name` and `city`] or `doe_reference_name` or [`blended_doe_reference_names` and `blended_doe_reference_percents`].

When only `doe_reference_name` is provided the `Site.latitude` and `Site.longitude` are used to look up the ASHRAE climate zone, which determines the appropriate DoE Commercial Reference Building profile.

When using the [`doe_reference_name` and `city`] option, choose `city` from one of the
cities used to represent the ASHRAE climate zones:
- Albuquerque
- Atlanta
- Baltimore
- Boulder
- Chicago
- Duluth
- Fairbanks
- Helena
- Houston
- LosAngeles
- Miami
- Minneapolis
- Phoenix
- SanFrancisco
- Seattle
and `doe_reference_name` from:
- FastFoodRest
- FullServiceRest
- Hospital
- LargeHotel
- LargeOffice
- MediumOffice
- MidriseApartment
- Outpatient
- PrimarySchool
- RetailStore
- SecondarySchool
- SmallHotel
- SmallOffice
- StripMall
- Supermarket
- Warehouse
- FlatLoad

Each `city` and `doe_reference_name` combination has a default `annual_kwh`, or you can provide your
own `annual_kwh` or `monthly_totals_kwh` and the reference profile will be scaled appropriately.
"""
mutable struct ElectricLoad  # mutable to adjust (critical_)loads_kw based off of (critical_)loads_kw_is_net
    loads_kw::Array{Real,1}
    year::Int  # used in ElectricTariff to align rate schedule with weekdays/weekends
    critical_loads_kw::Array{Real,1}
    loads_kw_is_net::Bool
    critical_loads_kw_is_net::Bool
    city::String

    function ElectricLoad(;
        loads_kw::Array{<:Real,1} = Real[],
        year::Int = 2020,
        doe_reference_name::String = "",
        blended_doe_reference_names::Array{String, 1} = String[],
        blended_doe_reference_percents::Array{<:Real,1} = Real[],
        city::String = "",
        annual_kwh::Union{Real, Nothing} = nothing,
        monthly_totals_kwh::Array{<:Real,1} = Real[],
        critical_loads_kw::Union{Missing, Array{Real,1}} = missing,
        loads_kw_is_net::Bool = true,
        critical_loads_kw_is_net::Bool = false,
        critical_load_pct::Real = 0.5,
        latitude::Float64,
        longitude::Float64,
        time_steps_per_hour::Int = 1
        )

        if length(loads_kw) > 0

            if !(length(loads_kw) / time_steps_per_hour ≈ 8760)
                @error "Provided electric load does not match the time_steps_per_hour."
            end

            if ismissing(critical_loads_kw)
                critical_loads_kw = critical_load_pct * loads_kw
            end

        elseif !isempty(doe_reference_name)
            # NOTE: must use year that starts on Sunday with DOE reference doe_ref_profiles
            if year != 2017
                @warn "Changing load profile year to 2017 because DOE reference profiles start on a Sunday."
            end
            year = 2017
            loads_kw = BuiltInElectricLoad(city, doe_reference_name, latitude, longitude, year, annual_kwh, monthly_totals_kwh)
            if ismissing(critical_loads_kw)
                critical_loads_kw = critical_load_pct * loads_kw
            end

        elseif length(blended_doe_reference_names) > 1 &&
            length(blended_doe_reference_names) == length(blended_doe_reference_percents)
            loads_kw = blend_and_scale_doe_profiles(BuiltInElectricLoad, latitude, longitude, year,
                                                    blended_doe_reference_names, blended_doe_reference_percents, city,
                                                    annual_kwh, monthly_totals_kwh)

            if ismissing(critical_loads_kw)
                critical_loads_kw = critical_load_pct * loads_kw
            end
        else
            error("Cannot construct ElectricLoad. You must provide either [loads_kw], [doe_reference_name, city],
                  [doe_reference_name, latitude, longitude],
                  or [blended_doe_reference_names, blended_doe_reference_percents] with city or latitude and longitude.")
        end

        if length(loads_kw) < 8760*time_steps_per_hour
            loads_kw = repeat(loads_kw, inner=time_steps_per_hour / (length(loads_kw)/8760))
            @info "Repeating electric loads in each hour to match the time_steps_per_hour."
        end

        new(
            loads_kw,
            year,
            critical_loads_kw,
            loads_kw_is_net,
            critical_loads_kw_is_net
        )
    end
end

function BuiltInElectricLoad(
    city::String,
    buildingtype::String,
    latitude::Float64,
    longitude::Float64,
    year::Int,
    annual_kwh::Union{<:Real, Nothing}=nothing,
    monthly_totals_kwh::Union{<:Real, Vector{<:Real}}=nothing,
    )

    annual_loads = Dict(
        "Albuquerque" => Dict(
            "fastfoodrest" => 193235,
            "fullservicerest" => 367661,
            "hospital" => 8468546,
            "largehotel" => 2407649,
            "largeoffice" => 6303595,
            "mediumoffice" => 884408,
            "midriseapartment" => 269734,
            "outpatient" => 1678720,
            "primaryschool" => 1070908,
            "retailstore" => 505417,
            "secondaryschool" => 2588879,
            "smallhotel" => 755373,
            "smalloffice" => 87008,
            "stripmall" => 497132,
            "supermarket" => 1947654,
            "warehouse" => 228939,
            "flatload" => 500000
        ),
        "Atlanta" => Dict(
            "fastfoodrest" => 197467,
            "fullservicerest" => 353750,
            "hospital" => 9054747,
            "largehotel" => 2649819,
            "largeoffice" => 6995864,
            "mediumoffice" => 929349,
            "midriseapartment" => 285349,
            "outpatient" => 1672434,
            "primaryschool" => 1128702,
            "retailstore" => 543340,
            "secondaryschool" => 2849901,
            "smallhotel" => 795777,
            "smalloffice" => 90162,
            "stripmall" => 529719,
            "supermarket" => 2092966,
            "warehouse" => 223009,
            "flatload" => 500000
        ),
        "Baltimore" => Dict(
            "fastfoodrest" => 192831,
            "fullservicerest" => 341893,
            "hospital" => 8895223,
            "largehotel" => 2534272,
            "largeoffice" => 6836130,
            "mediumoffice" => 945425,
            "midriseapartment" => 273225,
            "outpatient" => 1623103,
            "primaryschool" => 1077312,
            "retailstore" => 510257,
            "secondaryschool" => 2698987,
            "smallhotel" => 767538,
            "smalloffice" => 86112,
            "stripmall" => 504715,
            "supermarket" => 2018760,
            "warehouse" => 229712,
            "flatload" => 500000
        ),
        "Boulder" => Dict(
            "fastfoodrest" => 189092,
            "fullservicerest" => 334005,
            "hospital" => 8281865,
            "largehotel" => 2313151,
            "largeoffice" => 6127030,
            "mediumoffice" => 884726,
            "midriseapartment" => 255428,
            "outpatient" => 1621950,
            "primaryschool" => 1018424,
            "retailstore" => 504256,
            "secondaryschool" => 2441588,
            "smallhotel" => 736174,
            "smalloffice" => 84900,
            "stripmall" => 495018,
            "supermarket" => 1956244,
            "warehouse" => 243615,
            "flatload" => 500000
        ),
        "Chicago" => Dict(
            "fastfoodrest" => 189558,
            "fullservicerest" => 333659,
            "hospital" => 8567087,
            "largehotel" => 2402021,
            "largeoffice" => 6369028,
            "mediumoffice" => 972772,
            "midriseapartment" => 265528,
            "outpatient" => 1587062,
            "primaryschool" => 1045477,
            "retailstore" => 513106,
            "secondaryschool" => 2568086,
            "smallhotel" => 759657,
            "smalloffice" => 86224,
            "stripmall" => 506886,
            "supermarket" => 2025507,
            "warehouse" => 245750,
            "flatload" => 500000
        ),
        "Duluth" => Dict(
            "fastfoodrest" => 183713,
            "fullservicerest" => 318867,
            "hospital" => 8134328,
            "largehotel" => 2231678,
            "largeoffice" => 6036003,
            "mediumoffice" => 1032533,
            "midriseapartment" => 256393,
            "outpatient" => 1534322,
            "primaryschool" => 982163,
            "retailstore" => 532503,
            "secondaryschool" => 2333466,
            "smallhotel" => 752284,
            "smalloffice" => 83759,
            "stripmall" => 500979,
            "supermarket" => 1980986,
            "warehouse" => 256575,
            "flatload" => 500000
        ),
        "Fairbanks" => Dict(
            "fastfoodrest" => 182495,
            "fullservicerest" => 314760,
            "hospital" => 7899166,
            "largehotel" => 2181664,
            "largeoffice" => 5956232,
            "mediumoffice" => 1267132,
            "midriseapartment" => 271840,
            "outpatient" => 1620270,
            "primaryschool" => 986128,
            "retailstore" => 573411,
            "secondaryschool" => 2344790,
            "smallhotel" => 831480,
            "smalloffice" => 86614,
            "stripmall" => 545421,
            "supermarket" => 2033295,
            "warehouse" => 285064,
            "flatload" => 500000
        ),
        "Helena" => Dict(
            "fastfoodrest" => 185877,
            "fullservicerest" => 325263,
            "hospital" => 8068698,
            "largehotel" => 2246239,
            "largeoffice" => 6003137,
            "mediumoffice" => 930630,
            "midriseapartment" => 252659,
            "outpatient" => 1568262,
            "primaryschool" => 994496,
            "retailstore" => 534933,
            "secondaryschool" => 2357548,
            "smallhotel" => 729797,
            "smalloffice" => 84219,
            "stripmall" => 503504,
            "supermarket" => 1969137,
            "warehouse" => 252245,
            "flatload" => 500000
        ),
        "Houston" => Dict(
            "fastfoodrest" => 210283,
            "fullservicerest" => 383987,
            "hospital" => 9634661,
            "largehotel" => 3050370,
            "largeoffice" => 7539308,
            "mediumoffice" => 972535,
            "midriseapartment" => 335063,
            "outpatient" => 1756541,
            "primaryschool" => 1258146,
            "retailstore" => 589419,
            "secondaryschool" => 3421024,
            "smallhotel" => 863952,
            "smalloffice" => 98508,
            "stripmall" => 577987,
            "supermarket" => 2225265,
            "warehouse" => 221593,
            "flatload" => 500000
        ),
        "LosAngeles" => Dict(
            "fastfoodrest" => 188857,
            "fullservicerest" => 352240,
            "hospital" => 8498389,
            "largehotel" => 2458786,
            "largeoffice" => 6642784,
            "mediumoffice" => 846742,
            "midriseapartment" => 248028,
            "outpatient" => 1565008,
            "primaryschool" => 1095263,
            "retailstore" => 486188,
            "secondaryschool" => 2584380,
            "smallhotel" => 751880,
            "smalloffice" => 86655,
            "stripmall" => 491972,
            "supermarket" => 1935886,
            "warehouse" => 182085,
            "flatload" => 500000
        ),
        "Miami" => Dict(
            "fastfoodrest" => 224494,
            "fullservicerest" => 448713,
            "hospital" => 10062043,
            "largehotel" => 3437188,
            "largeoffice" => 8002063,
            "mediumoffice" => 1021224,
            "midriseapartment" => 424956,
            "outpatient" => 1929148,
            "primaryschool" => 1426635,
            "retailstore" => 635086,
            "secondaryschool" => 4074081,
            "smallhotel" => 972090,
            "smalloffice" => 108279,
            "stripmall" => 675793,
            "supermarket" => 2260929,
            "warehouse" => 202082,
            "flatload" => 500000
        ),
        "Minneapolis" => Dict(
            "fastfoodrest" => 188368,
            "fullservicerest" => 330920,
            "hospital" => 8425063,
            "largehotel" => 2378872,
            "largeoffice" => 6306693,
            "mediumoffice" => 1005875,
            "midriseapartment" => 267383,
            "outpatient" => 1582701,
            "primaryschool" => 1022667,
            "retailstore" => 539203,
            "secondaryschool" => 2498647,
            "smallhotel" => 774571,
            "smalloffice" => 85921,
            "stripmall" => 511567,
            "supermarket" => 2034650,
            "warehouse" => 249332,
            "flatload" => 500000
        ),
        "Phoenix" => Dict(
            "fastfoodrest" => 216088,
            "fullservicerest" => 389739,
            "hospital" => 9265786,
            "largehotel" => 2990053,
            "largeoffice" => 7211666,
            "mediumoffice" => 1004988,
            "midriseapartment" => 378378,
            "outpatient" => 1849358,
            "primaryschool" => 1289084,
            "retailstore" => 593924,
            "secondaryschool" => 3503727,
            "smallhotel" => 881843,
            "smalloffice" => 104583,
            "stripmall" => 590954,
            "supermarket" => 2056195,
            "warehouse" => 241585,
            "flatload" => 500000
        ),
        "SanFrancisco" => Dict(
            "fastfoodrest" => 183953,
            "fullservicerest" => 317124,
            "hospital" => 7752817,
            "largehotel" => 2206880,
            "largeoffice" => 6085403,
            "mediumoffice" => 792199,
            "midriseapartment" => 229671,
            "outpatient" => 1394447,
            "primaryschool" => 1009369,
            "retailstore" => 449025,
            "secondaryschool" => 2327074,
            "smallhotel" => 698095,
            "smalloffice" => 78132,
            "stripmall" => 455802,
            "supermarket" => 1841655,
            "warehouse" => 185889,
            "flatload" => 500000
        ),
        "Seattle" => Dict(
            "fastfoodrest" => 184142,
            "fullservicerest" => 318741,
            "hospital" => 7912504,
            "largehotel" => 2212410,
            "largeoffice" => 6019271,
            "mediumoffice" => 878390,
            "midriseapartment" => 237242,
            "outpatient" => 1434195,
            "primaryschool" => 983498,
            "retailstore" => 455854,
            "secondaryschool" => 2282972,
            "smallhotel" => 693921,
            "smalloffice" => 79716,
            "stripmall" => 460449,
            "supermarket" => 1868973,
            "warehouse" => 210300,
            "flatload" => 500000
        ),
    )
    if !(buildingtype in default_buildings)
        error("buildingtype $(buildingtype) not in $(default_buildings).")
    end

    if isempty(city)
        city = find_ashrae_zone_city(latitude, longitude)
    end

    if isnothing(annual_kwh)
        annual_kwh = annual_loads[city][lowercase(buildingtype)]
    end

    built_in_load("electric", city, buildingtype, year, annual_kwh, monthly_totals_kwh)
end

max_stars_repo_name: NREL/REoptLite

max_issues_repo_name: NREL/REoptLite

max_forks_repo_name: NREL/REoptLite


text:
using LightGraphs
using BALUtils
using StaticGraphs

ba = readbal(ARGS[1])
dia = diameter(StaticGraph(Graph(visibility_graph(ba))))
open(ARGS[2], "w") do io
    write(io, "diameter\n$dia\n")
end


text:
# This implements src/snes/examples/tutorials/ex2.c from PETSc using the PETSc.jl package, using SNES
#
# This is the same as SNES_ex2b.j, except that we show how automatic differentiation can be used to
# compute the jacobian.
#
# Newton method to solve u'' + u^{2} = f, sequentially.

using PETSc, MPI, LinearAlgebra, SparseArrays, Plots, ForwardDiff

if ~MPI.Initialized()
    MPI.Init()
end

PETSc.initialize()

```
    Computes initial guess
```
function FormInitialGuess!(x)
    for i=1:length(x)
        x[i] = 0.50;
    end
end

```
    Computes the residual f, given solution vector x
```
function FormResidual!(f,x)
    n       =   length(x);
    xp      =   LinRange(0.0,1.0, n);
    F       =   6.0.*xp .+ (xp .+1.e-12).^6.0;      # define source term function

    dx      =   1.0/(n-1.0);
    f[1]    =   x[1] - 0.0;
    for i=2:n-1
        f[i] = (x[i-1] - 2.0*x[i] + x[i+1])/dx^2 + x[i]*x[i] - F[i]
    end
    f[n]    =   x[n] - 1.0;

end

```
    Wrapper which makes it easier to compute the jacobian using automatic differntiation
```
function  ForwardDiff_res(x)

    f   = zero(x)               # vector of zeros, of same type as e
    FormResidual!(f,x);

    return f;
end

```
    Computes the jacobian, given solution vector x
```
function FormJacobian!(x, args...)

    J        =  args[1];        # preconditioner = args[2], in case we want it to be different from J

    # Use AD to compute jacobian; by transferring x into sparse, the output will be sparse
    J_julia  =  ForwardDiff.jacobian(ForwardDiff_res,sparse(x));

    J       .=  J_julia;   # copy to petsc format
end

# ==========================================
# Main code

# Compute initial solution
n   =   101;
x   =   zeros(n);

FormInitialGuess!(x);

# Compute initial jacobian using a julia structure to obtain the nonzero structure
# Note that we can also obtain this structure in a different manner
Jstruct  = zeros(n,n);
FormJacobian!(x, Jstruct);                              # jacobian in julia form
Jsp      =   sparse(Float64.(abs.(Jstruct) .> 0))       # sparse julia, with 1.0 in nonzero spots
PJ       =   PETSc.MatSeqAIJ(Jsp);                      # transfer to PETSc (initialize matrix with correct nonzero pattern)

# Setup SNES
x_s = PETSc.VecSeq(x);                  # solution vector
res = PETSc.VecSeq(zeros(size(x)));     # residual vector

S = PETSc.SNES{Float64}(MPI.COMM_SELF;
        snes_rtol=1e-12,
        snes_monitor=true,
        snes_converged_reason=true);
PETSc.setfunction!(S, FormResidual!, res)
PETSc.setjacobian!(S, FormJacobian!, PJ, PJ)

# solve
PETSc.solve!(x_s, S);

# Extract & plot solution
x_sol = x_s.array;                  # convert solution to julia format
FormResidual!(res.array,x_sol)      # just for checking, compute residual
@show norm(res.array)

plot(LinRange(0,1,n),x_sol,xlabel="width",ylabel="solution")

#PETSc.finalize()


text:
const graphs = [
:bull,
:chvatal,
:house
]
@testset "Connectivity" begin
     for s in graphs
         gx = smallgraph(s)

         for g in testgraphs(gx)
             a = adjacency_matrix(g)
             x = connected_components(g)
             y = connected_components(SimpleWeightedGraph(a))
             @test x == y
         end
     end
 end


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

using ParallelAccelerator
using DocOpt

#ParallelAccelerator.set_debug_level(3)
#ParallelAccelerator.DomainIR.set_debug_level(3)
#ParallelAccelerator.CGen.set_debug_level(3)

@acc function calcPi(n)
    x = rand(n) * 2.0 - 1.0
    y = rand(n) * 2.0 - 1.0
    return 4.0*sum(x.^2 + y.^2 .< 1.0)/n
end

function main()
    doc = """pi.jl

Estimate Pi using Monte Carlo method.

Usage:
  pi.jl -h | --help
  pi.jl [--points=<points>]

Options:
  -h --help          Show this screen.
  --points=<points>  Specify number of generated random points [default: 10000000].
"""
    arguments = docopt(doc)

    points = parse(Int, arguments["--points"])

    println("points = ", points)

    tic()
    calcPi(100)
    println("SELFPRIMED ", toq())

    tic()
    pi_val = calcPi(points)
    time = toq()
    println("pi = ", pi_val)
    println("SELFTIMED ", time)

end

main()


text:
using Tensorial
using BenchmarkTools
using InteractiveUtils # for `versioninfo`
using Markdown

const run_array = true

const suite = BenchmarkGroup()
suite["Array"] = BenchmarkGroup()
suite["Tensor"] = BenchmarkGroup()

const tensors = let
    dict = Dict{Tuple{Space, DataType}, Tensor}()
    for S in (Space(3),
              Space(3,3),
              Space(Symmetry(3,3)),
              Space(3,3,3),
              Space(Symmetry(3,3,3)),
              Space(3,3,3,3),
              Space(Symmetry(3,3),Symmetry(3,3)))
        for T in (Float32, Float64)
            dict[(S, T)] = rand(Tensorial.tensortype(S){T})
        end
    end
    dict
end

include("ops.jl")

tune!(suite)
const results = run(suite)

let path = "../docs/src/Benchmarks.md"

    path = joinpath(dirname(@__FILE__), path)

    open(path, "w") do file

        function printheader(head)
            println(file, "| **$(head)** | | | |")
        end
        function printrow(op, tt, ta)
            pretty = (t) -> BenchmarkTools.prettytime(BenchmarkTools.time(minimum(t)))
            speedup = (ta, tt) -> round(10*BenchmarkTools.time(minimum(ta))/BenchmarkTools.time(minimum(tt)))/10
            println(file, "| `$(op)` | $(pretty(tt)) | $(pretty(ta)) | ×$(speedup(ta, tt)) |")
        end

        print(file,
              """
              # Benchmarks

              ```julia
              a = rand(Vec{3})
              A = rand(SecondOrderTensor{3})
              S = rand(SymmetricSecondOrderTensor{3})
              B = rand(Tensor{Tuple{3,3,3}})
              AA = rand(FourthOrderTensor{3})
              SS = rand(SymmetricFourthOrderTensor{3})
              ```

              | Operation  | `Tensor` | `Array` | speed-up |
              |:-----------|---------:|--------:|---------:|
              """)

        printheader("Single contraction")
        printrow("a ⋅ a", results["Tensor"]["Space(3) ⋅ Space(3)"],
                          results["Array" ]["Space(3) ⋅ Space(3)"])
        printrow("A ⋅ a", results["Tensor"]["Space(3,3) ⋅ Space(3)"],
                          results["Array" ]["Space(3,3) ⋅ Space(3)"])
        printrow("S ⋅ a", results["Tensor"]["Space(Symmetry(3,3)) ⋅ Space(3)"],
                          results["Array" ]["Space(Symmetry(3,3)) ⋅ Space(3)"])

        printheader("Double contraction")
        printrow("A ⊡ A", results["Tensor"]["Space(3,3) ⊡ Space(3,3)"],
                          results["Array" ]["Space(3,3) ⊡ Space(3,3)"])
        printrow("S ⊡ S", results["Tensor"]["Space(Symmetry(3,3)) ⊡ Space(Symmetry(3,3))"],
                          results["Array" ]["Space(Symmetry(3,3)) ⊡ Space(Symmetry(3,3))"])
        printrow("B ⊡ A", results["Tensor"]["Space(3,3,3) ⊡ Space(3,3)"],
                          results["Array" ]["Space(3,3,3) ⊡ Space(3,3)"])
        printrow("AA ⊡ A", results["Tensor"]["Space(3,3,3,3) ⊡ Space(3,3)"],
                           results["Array" ]["Space(3,3,3,3) ⊡ Space(3,3)"])
        printrow("SS ⊡ S", results["Tensor"]["Space(Symmetry(3,3),Symmetry(3,3)) ⊡ Space(Symmetry(3,3))"],
                           results["Array" ]["Space(Symmetry(3,3),Symmetry(3,3)) ⊡ Space(Symmetry(3,3))"])

        printheader("Tensor product")
        printrow("a ⊗ a", results["Tensor"]["Space(3) ⊗ Space(3)"],
                          results["Array" ]["Space(3) ⊗ Space(3)"])

        printheader("Cross product")
        printrow("a × a", results["Tensor"]["Space(3) ⊗ Space(3)"],
                          results["Array" ]["Space(3) ⊗ Space(3)"])

        printheader("Determinant")
        printrow("det(A)", results["Tensor"]["det(Space(3,3))"],
                           results["Array" ]["det(Space(3,3))"])
        printrow("det(S)", results["Tensor"]["det(Space(Symmetry(3,3)))"],
                           results["Array" ]["det(Space(Symmetry(3,3)))"])

        printheader("Inverse")
        printrow("inv(A)", results["Tensor"]["inv(Space(3,3))"],
                           results["Array" ]["inv(Space(3,3))"])
        printrow("inv(S)", results["Tensor"]["inv(Space(Symmetry(3,3)))"],
                           results["Array" ]["inv(Space(Symmetry(3,3)))"])
        printrow("inv(AA)", results["Tensor"]["inv(Space(3,3,3,3))"],
                            results["Array" ]["inv(Space(3,3,3,3))"])
        printrow("inv(SS)", results["Tensor"]["inv(Space(Symmetry(3,3),Symmetry(3,3)))"],
                            results["Array" ]["inv(Space(Symmetry(3,3),Symmetry(3,3)))"])

        # print versioninfo
        print(file,
              """

              The benchmarks are generated by
              [`runbenchmarks.jl`](https://github.com/KeitaNakamura/Tensorial.jl/blob/master/benchmark/runbenchmarks.jl)
              on the following system:

              ```julia
              julia> versioninfo()
              $((x = IOBuffer(); versioninfo(x); String(take!(x))))
              """)
        println(file, "```")
    end

    # display results
    open(path, "r") do file
        display(Markdown.parse(file))
    end
end


text:
using FinElt
using FinElt.PlanarPoisson

run(`gmsh -2 -o complic_keyhole.msh -clmax 0.05 -format msh22
complic_keyhole.geo`)
mesh = read_msh_file("complic_keyhole.msh")
essential_bc = [ "DarkBlue", "Black" ]
vp = VariationalProblem(mesh, essential_bc)
g(x) = -hypot(x[1],x[2])/2
assign_bdry_vals!(vp, "DarkBlue", g)
add_bilin_form!(vp,  "LightBlue", grad_dot_grad!, 1.0)
add_bilin_form!(vp,  "LightGreen", grad_dot_grad!, 10.0)
add_lin_functnl!(vp, "LightBlue", source_times_func!, x->1.0)
add_lin_functnl!(vp, "LightGreen", source_times_func!, x->4.0)
add_lin_functnl!(vp, "Violet",  bdry_source_times_func!, x->-1.0)

A, b = assembled_linear_system(vp)

ufree = A \ b

u = complete_soln(ufree, vp)
write_pos_file("complic_keyhole.pos") do fid
    save_warp_nodal_scalar_field(u, "u", fid)
end

run(`gmsh complic_keyhole.script`)


text:
dir = joinpath(Pkg.dir("Augmentor"), "test/sampledir")

@testset "listfiles" begin
    lst = Augmentor.listfiles(dir;
            expand=false, hidden=false, recursive=false)
    @test typeof(lst) <: Vector{UTF8String}
    @test length(lst) == 2

    lst = Augmentor.listfiles(dir;
            expand=false, hidden=true, recursive=false)
    @test typeof(lst) <: Vector{UTF8String}
    @test length(lst) == 3

    lst = Augmentor.listfiles(dir;
            expand=false, hidden=false, recursive=true)
    @test typeof(lst) <: Vector{UTF8String}
    @test length(lst) == 3

    lst = Augmentor.listfiles(dir;
            expand=false, hidden=true, recursive=true)
    @test typeof(lst) <: Vector{UTF8String}
    @test length(lst) == 5

    lst = Augmentor.listfiles(dir;
            expand=true, hidden=true, recursive=true)
    @test typeof(lst) <: Vector{UTF8String}
    @test length(lst) == 5
end

@testset "DirImageSource" begin
    @test DirImageSource <: ImageSource
    src = DirImageSource(dir)
    @test length(src) == 3
    @test endof(src) == 3
    @test src[end] == src[3]
    @test eltype(src) <: Image
    show(src)
    println()

    img = rand(src)
    @test typeof(img) <: Image

    imgs = rand(src, 2)
    @test typeof(imgs) <: Vector
    @test eltype(imgs) <: Image

    @testset "iterator interface" begin
        count = 0
        for img in src
            count += 1
            @test typeof(img) <: Image
        end
        @test count == 3
    end
end


text:
using StaticCompiler
using Test
using Libdl
using LinearAlgebra

@testset "Basics" begin

    simple_sum(x) = x + one(typeof(x))

    # This probably needs a macro
    @test ccall(generate_shlib_fptr(simple_sum, (Int,)), Int, (Int,), 1) == Int(2)
    @test ccall(generate_shlib_fptr(simple_sum, (Float64,)), Float64, (Float64 ,), 1) == Float64(2)

    @test ccall(generate_shlib_fptr(simple_sum, (Int32,)), Int32, (Int32,), 1) == Int32(2)
    @test ccall(generate_shlib_fptr(simple_sum, (Float32,)), Float32, (Float32 ,), 1) == Float16(2)

    @test ccall(generate_shlib_fptr(simple_sum, (Int16,)), Int16, (Int16,), 1) == Int16(2)
    @test ccall(generate_shlib_fptr(simple_sum, (Float16,)), Float16, (Float16 ,), 1) == Float16(2)

end

@testset "Recursion" begin
    fib(n) = n <= 1 ? n : fib(n - 1) + fib(n - 2)
    # This works on the REPL but fails here
    @test_skip ccall(generate_shlib_fptr(fib, (Int,)), Int, (Int,), 10) == 55
end

# Call binaries for testing
# @testset "Generate binary" begin
#     fib(n) = n <= 1 ? n : fib(n - 1) + fib(n - 2)
#     libname = tempname()
#     generate_shlib(fib, (Int,), libname)
#     ptr = Libdl.dlopen(libname * "." * Libdl.dlext, Libdl.RTLD_LOCAL)
#     fptr = Libdl.dlsym(ptr, "julia_fib")
#     @assert fptr != C_NULL
#     # This works on REPL
#     @test_skip ccall(fptr, Int, (Int,), 10) == 55
# end

@testset "Loops" begin
    function sum_first_N_int(N)
        s = 0
        for a in 1:N
            s += a
        end
        s
    end
    @test ccall(generate_shlib_fptr(sum_first_N_int, (Int,)), Int, (Int,), 10) == 55

    function sum_first_N_float64(N)
        s = Float64(0)
        for a in 1:N
            s += Float64(a)
        end
        s
    end
    @test ccall(generate_shlib_fptr(sum_first_N_float64, (Int,)), Float64, (Int,), 10) == 55.

    function sum_first_N_int_inbounds(N)
        s = 0
        @inbounds for a in 1:N
            s += a
        end
        s
    end
    @test ccall(generate_shlib_fptr(sum_first_N_int_inbounds, (Int,)), Int, (Int,), 10) == 55

    function sum_first_N_float64_inbounds(N)
        s = Float64(0)
        @inbounds for a in 1:N
            s += Float64(a)
        end
        s
    end
    @test ccall(generate_shlib_fptr(sum_first_N_float64_inbounds, (Int,)), Float64, (Int,), 10) == 55.

end

# Arrays with different input types Int32, Int64, Float32, Float64, Complex?
@testset "Arrays" begin

    arr = collect(1:10)
    function array_sum(n, A)
        s = zero(eltype(A))
        for i in 1:n
            s += A[i]
        end
        s
    end

    #This segfaults, not sure if this is how you pass around arrays
    @test_skip ccall(generate_shlib_fptr(array_sum, (Csize_t, Ptr{Float64})), Int, (Csize_t, Ptr{Float64}), length(arr), arr) == 55

end

# Just to call external libraries
@testset "BLAS" begin
    function mydot(N)
        a = Float64.(1:N)
        BLAS.dot(N, a, 1, a, 1)
    end
    @test_skip ccall(generate_shlib_fptr(mydot, (Int,)), Float64, (Int,), 2) == 5.
end

@testset "Hello World" begin
    function hello(N)
        println("Hello World $N")
        N
    end
    # How do I test this?
    # Also ... this segfaults
    @test_skip ccall(generate_shlib_fptr(hello, (Int,)), Int, (Int,), 1) == 1
end

# data structures, dictionaries, tuples, named tuples
# passing pointers?
# @inbounds LoopVectorization


text:
module jjson

include("./lexer.jl")
include("./parser.jl")

using Mmap

function popnfirst!(a::Array{T,1}, n::Int64) where T
    for i in 1:n
        popfirst!(a)
    end
end

const JSON_COMMA        = ','
const JSON_COLON        = ':'
const JSON_LEFTBRACKET  = '['
const JSON_RIGHTBRACKET = ']'
const JSON_LEFTBRACE    = '{'
const JSON_RIGHTBRACE   = '}'
const JSON_QUOTE        = '"'
const JSON_QUOTE      = '"'
const JSON_WHITESPACE = [' ', '\t', '\b', '\n', '\r']
const JSON_SYNTAX     = [JSON_COMMA, JSON_COLON, JSON_LEFTBRACKET,
                         JSON_RIGHTBRACKET, JSON_LEFTBRACE, JSON_RIGHTBRACE]

function parse_file(filename::String)
    tokens = lexer(filename)
    return parser!(tokens)
end # parse_file

end # module

max_stars_repo_name: Iaggelis/jjson

max_issues_repo_name: Iaggelis/jjson

max_forks_repo_name: Iaggelis/jjson


text:
module HillsEigenvalues
using SparseArrays, LinearAlgebra, Arpack
# Write your package code here.

include("eig_system.jl")

export matrix_construct, eigval_iterator

end

max_stars_repo_name: Mikejmnez/HillsEigenvalues

max_issues_repo_name: Mikejmnez/HillsEigenvalues

max_forks_repo_name: Mikejmnez/HillsEigenvalues


text:
using SMACOF
using Plots
using Random
Random.seed!(2021)
gr()

script_dir = dirname(@__FILE__)

X = rand(400, 2)
X = X[.!((1 / 4 .< X[:,1] .< 3 / 4) .& (1 / 4 .< X[:,2] .< 3 / 4)), :]
# Xhist = wda_smacof(dist(X), η = 0.8, verbose = true, anchors = X, ε = 1e-6, return_history = true)[2]
sm = smacof(dist(X), Xinit = randn(size(X)) * 0.1, anchors = X, return_hist = true)
Xhist = sm.hist
# Make mp4 and then convert to gif
anim_dir = joinpath(dirname(script_dir), "tmp")
animation = Animation(anim_dir, String[])
lst_mse  = []
lst_str = sm.σ / maximum(sm.σ) * 0.35
for i ∈ 1:length(Xhist)
    Xi = align(Xhist[i], X)
    push!(lst_mse, mse(Xi, X))

    p1 = scatter(X[:, 1], X[:, 2], alpha = 0.7, markersize = 8, legend = nothing);
    scatter!(Xi[:,1], Xi[:, 2]);
    xlims!((-0.2, 1.2))
    ylims!((-0.2, 1.2))
    title!("$i")

    p2 = plot(1:i, lst_mse, label = "mse", legend = :best)
    plot!(p2, 1:i, lst_str[1:i], label = "stress")
    xlims!((0, length(Xhist)))
    ylims!((0, 0.35))
    title!("Performance")
    p = plot(p1, p2, layout = grid(1, 2, widths = [.6, .4]))
    plot!(size = (600, 400))
    frame(animation, p)
end

gif_path = joinpath(script_dir, "simple")
run(`ffmpeg -y -r 15 -i $anim_dir/"%06d.png" -vcodec libx264 -crf 25 "$(gif_path).mp4"`)
run(`ffmpeg -y -i $(gif_path).mp4 -vf "fps=10,scale=1000:-1:flags=lanczos,split[s0][s1];[s0]palettegen[p];[s1][p]paletteuse" -loop 0 $(gif_path).gif`)
run(`rm $(gif_path).mp4`)

# Clear image folder
for f in readdir(anim_dir)
    println(f)
    run(`rm $(joinpath(anim_dir, f))`)
end


text:
using InMemoryDatasets, PooledArrays, Random, Test, CategoricalArrays

@testset "Basic usage" begin
    ds = Dataset(x1 = rand(1:1000,100))
    dss = sort(ds, 1)
    @test issorted(dss.x1)
    dss = sort(ds, 1, rev = true)
    @test issorted(dss.x1, rev = true)
    setformat!(ds, 1=>isodd)
    dss = sort(ds, 1)
    @test issorted(dss.x1, by = isodd)
    dss = sort(ds, 1, rev = true)
    @test issorted(dss.x1, by = isodd, rev = true)
    dss = sort(ds, 1, rev = true, mapformats = false)
    @test issorted(dss.x1, rev = true)
    ds = Dataset(x1 = rand(3333)*100)
    dss = sort(ds, 1)
    @test issorted(dss.x1)
    dss = sort(ds, 1, rev = true)
    @test issorted(dss.x1, rev = true)
    setformat!(ds, 1=>round)
    dss = sort(ds, 1)
    @test issorted(dss.x1, by = round)
    dss = sort(ds, 1, rev = true)
    @test issorted(dss.x1, by = round, rev = true)
    dss = sort(ds, 1, rev = true, mapformats = false)
    @test issorted(dss.x1, rev = true)
    x1 = Vector{Int}(undef, 5050)
    y1 = Vector{Float64}(undef, 5050)
    cnt = 1
    for i in 1:100
        for j in 1:i
            x1[cnt] = i
            y1[cnt] = j + rand()
            cnt += 1
        end
    end
    shuffle_row = shuffle(1:5050)
    ds = Dataset(x = x1[shuffle_row], y = y1[shuffle_row])
    sort!(ds, 1:2)
    @test ds.x == x1
    @test ds.y == y1
    sort!(ds,2)
    @test issorted(ds.y)
    ds2 = ds[shuffle_row, :]
    sort!(ds2, 1)
    ds3 = sort(ds2, 1:2)
    @test ds3 == sort(ds, :)
    ds4 = sort(ds, 1:2, rev = true)
    @test ds4.x == reverse(x1)
    @test ds4.y == reverse(y1)
    ds = Dataset(Float64.(rand(1:3, 100, 4)), :auto)
    map!(ds, x->rand()<.3 ? missing : x, :)
    ds2 = copy(ds)
    for j in ncol(ds2):-1:1
        r = sortperm(ds2[!, j].val)
        for i in 1:ncol(ds2)
            ds2[!, i] = ds2[r, i]
        end
    end
    sort!(ds, 1:2)
    @test ds.x1 == ds2.x1
    @test ds.x2 == ds2.x2
    @test sort(ds, 1:4) == ds2

    x1 = rand(1:3, 777)
    x2 = rand(Date(100):Day(1):Date(101), 777)
    x3 = PooledArray(rand([missing, 1.1, 20000.0, 123.0], 777))
    x4 = [randstring(2) for _ in 1:777]
    x5 = rand(1:2, 777)
    ds = Dataset(col1 = x3, col2 = x2, col3 = x1, col4 = x5, col5 = x4)

    ds2 = copy(ds)
    for j in ncol(ds2):-1:1
        r = sortperm(ds2[!, j].val)
        for i in 1:ncol(ds2)
            ds2[!, i] = ds2[r, i]
        end
    end
    @test sort(ds, 1:5) == ds2
    sort!(ds, :)

    @test sort(ds,1, stable = false).col1 == ds2.col1
    sort!(ds, 5)
    @test issorted(ds.col5)
    @test issorted(sort(ds, 5, rev = true).col5, rev = true)

    dv1 = [9, 1, 8, missing, 3, 3, 7, missing]
    dv2 = [9, 1, 8, missing, 3, 3, 7, missing]
    dv3 = Vector{Union{Int, Missing}}(1:8)
    cv1 = CategoricalArray(dv1, ordered=true)

    d = Dataset(dv1 = dv1, dv2 = dv2, dv3 = dv3, cv1 = cv1)

    @test sortperm(d, :) == sortperm(dv1)
    @test sortperm(d[:, [:dv3, :dv1]], :) == sortperm(dv3)
    @test sort(d, :dv1)[!, :dv3] == sort(d, "dv1")[!, "dv3"] == sortperm(dv1)
    @test sort(d, :dv2)[!, :dv3] == sortperm(dv1)
    @test sort(d, :cv1)[!, :dv3] == sortperm(dv1)
    @test sort(d, [:dv1, :cv1])[!, :dv3] == sortperm(dv1)
    @test sort(d, [:dv1, :dv3])[!, :dv3] == sortperm(dv1)

    x = CategoricalArray{Union{Characters{6, UInt8}, String}}(["Old", "Young", "Middle", "Young"])
    levels!(x, ["Young", "Middle", "Old"])
    ds = Dataset(x = x)
    ds_s = sort(ds, :x)
    @test ds_s.x == ["Young", "Young", "Middle", "Old"]
    ordered!(x, true)
    ds = Dataset(x = x)
    ds_s = sort(ds, :x)
    @test ds_s.x == ["Young", "Young", "Middle", "Old"]
    a = rand(1:3, 1000)
    c = rand(1:10, 1000)
    f = rand(1:2, 1000)
    pc = PooledArray(c)
    pa = PooledArray(a)
    pf = categorical(f)
    ds = Dataset(a = a, c = c, f = f)
    ds_pa = Dataset(a = pa, c = pc, f = pf)
    ds2 = copy(ds)
    for j in ncol(ds2):-1:1
        r = sortperm(ds2[!, j].val)
        for i in 1:ncol(ds2)
            ds2[!, i] = ds2[r, i]
        end
    end
    for i in 1:20
        @test sort(ds, :) == sort(ds_pa, :) == ds2
    end
    x1 = -rand(1:1000, 5000)
    x2 = -rand(1:100, 5000)
    dsl = Dataset(x1 = Characters{6, UInt8}.(c"id" .* string.(-x1)), x2 = Characters{5, UInt8}.(c"id" .* string.(-x2)))
    dsr = Dataset(x1 = x1, x2 = x2)
    for i in 1:2
        dsl[!, i] = PooledArray(dsl[!, i])
        dsr[!, i] = PooledArray(dsr[!, i])
    end
    fmtfun3(x) = @views -parse(Int, x[3:end])
    setformat!(dsl, 1:2=>fmtfun3)
    @test sortperm(dsl, :)==sortperm(dsr, :)
    @test sortperm(dsl, :, rev=true)==sortperm(dsr, :, rev=true)

    ds = Dataset(rand(-100:1000000, 10^7, 3), :auto)
    sort!(ds, :)
    @test issorted(ds, :)
    sort!(ds, :, rev = true)
    @test issorted(ds, :, rev = true)
    sort!(ds, :, rev = [true, false, true])
    @test issorted(ds, :, rev = [true, false, true])

    ds = Dataset(x = Int8[-128, 1,2,1,127, -127])
    @test sort(ds, 1) == ds[[1,6,2,4,3,5], :]
    @test sort(ds, 1, rev = true) == ds[[5,3,4,2,6,1], :]
end

@testset "_find_starts_of_groups" begin
    # Suppose there will be at least one row and one column.
    ds = Dataset(x1 = [1])
    colsidx = [1]
    T = nrow(ds) < typemax(Int32) ? Val(Int32) : Val(Int64)
    ranges = [1]
    last_valid_index = 1
    r1, r2, r3 = IMD._find_starts_of_groups(ds, colsidx, T)
    @test r1 == colsidx
    @test r3 == last_valid_index
    @test r2 == ranges[1:r3]

    ds = Dataset(x1 = [1,1,1,3,3,1,1])
    colsidx = [1]
    T = nrow(ds) < typemax(Int32) ? Val(Int32) : Val(Int64)
    ranges = [1, 4, 6]
    last_valid_index = 3
    r1, r2, r3 = IMD._find_starts_of_groups(ds, colsidx, T)
    @test r1 == colsidx
    @test r3 == last_valid_index
    @test r2[1:r3] == ranges

    ds = Dataset(x1 = [1], x2 = [1], x3 = [2], x4 = [4])
    colsidx = 1:3
    T = nrow(ds) < typemax(Int32) ? Val(Int32) : Val(Int64)
    ranges = [1]
    last_valid_index = 1
    r1, r2, r3 = IMD._find_starts_of_groups(ds, colsidx, T)
    @test r1 == colsidx
    @test r3 == last_valid_index
    @test r2[1:r3] == ranges

    # Not only for Dataset, but also GroupBy.
    ds = Dataset(x1 = [1, 1, 1], x2 = [1, 2, 1])
    gb = groupby(ds, 2)
    colsidx = [2]
    T = nrow(ds) < typemax(Int32) ? Val(Int32) : Val(Int64)
    ranges = [1, 3]
    last_valid_index = 2
    r1, r2, r3 = IMD._find_starts_of_groups(gb, colsidx, T)
    @test r1 == colsidx
    @test r3 == last_valid_index
    @test r2[1:r3] == ranges

    # Use formats for some columns.
    format1(x) = isodd(x)
    ds = Dataset(x1 = [1, 1, 1, 1, 3, 3, 3], x2 = [1, 6, 5, 5, 5, 5, 2], x3 = [2, 2, 4, 4, 4, 4, 1], x4 = [4, 1, 1, 4, 4, 4, 1])
    setformat!(ds, [1, 3, 4] => format1)
    colsidx = 1:3
    T = nrow(ds) < typemax(Int32) ? Val(Int32) : Val(Int64)
    ranges = [1, 2, 3, 5, 7]
    last_valid_index = 5
    r1, r2, r3 = IMD._find_starts_of_groups(ds, colsidx, T, mapformats = false) # Do not use formatted values.
    @test r1 == colsidx
    @test r3 == last_valid_index
    @test r2[1:r3] == ranges

    ranges = [1, 2, 3, 7]
    last_valid_index = 4
    r1, r2, r3 = IMD._find_starts_of_groups(ds, colsidx, T) # Use formatted values.
    @test r1 == colsidx
    @test r3 == last_valid_index
    @test r2[1:r3] == ranges

    # Consider whether ds is sorted using formatted values.
    # If ds is sorted using formatted values.
    format2(x) = abs(2.5 - x)
    ds = Dataset(x1 = [1, 1, 1, 1, 3, 3, 3], x2 = [1, 6, 5, 5, 5, 5, 2], x3 = [2, 1, 4, 4, 4, 4, 1], x4 = [4, 1, 1, 4, 4, 4, 1])
    setformat!(ds, 3:4 => format2)
    dss = sort(ds, 3, mapformats = true) # Use formatted values for sorting.
    colsidx = [1, 3, 4]
    T = nrow(dss) < typemax(Int32) ? Val(Int32) : Val(Int64)
    ranges = [1, 2, 3, 4, 5, 7]
    last_valid_index = 6
    r1, r2, r3 = IMD._find_starts_of_groups(dss, colsidx, T, mapformats = false) # Do not use formatted values.
    @test r1 == colsidx
    @test r3 == last_valid_index
    @test r2[1:r3] == ranges

    ranges = [1, 2, 5]
    last_valid_index = 3
    r1, r2, r3 = IMD._find_starts_of_groups(dss, colsidx, T) # Use formatted values.
    @test r1 == colsidx
    @test r3 == last_valid_index
    @test r2[1:r3] == ranges

    # If ds is sorted using unformatted values.
    ds = Dataset(x1 = [1, 1, 1, 1, 3, 3, 3], x2 = [1, 6, 5, 5, 5, 5, 2], x3 = [2, 1, 4, 4, 4, 4, 1], x4 = [4, 1, 1, 4, 4, 4, 1])
    setformat!(ds, 3:4 => format2)
    dss = sort(ds, 3, mapformats = false) # Use unformatted values for sorting.
    colsidx = [1, 3, 4]
    T = nrow(dss) < typemax(Int32) ? Val(Int32) : Val(Int64)
    ranges = [1, 2, 3, 4, 5, 6]
    last_valid_index = 6
    r1, r2, r3 = IMD._find_starts_of_groups(dss, colsidx, T, mapformats = false) # Do not use formatted values.
    @test r1 == colsidx
    @test r3 == last_valid_index
    @test r2[1:r3] == ranges

    ranges = [1, 2, 3, 4, 6]
    last_valid_index = 5
    r1, r2, r3 = IMD._find_starts_of_groups(dss, colsidx, T) # Use formatted values.
    @test r1 == colsidx
    @test r3 == last_valid_index
    @test r2[1:r3] == ranges

    # Consider ds with multiple types.
    c1 = PooledArray(["string", "string", 1.1, 1.1, 1.1, 20000.0, 123.0])
    c2 = PooledArray(["string", missing, 1.1, 1.1, 'a', 'a', 'b'])
    c3 = PooledArray([missing, missing, missing, 1.1, 1.1, 20000.0, 123.0])
    c4 = CategoricalArray{Union{Characters{6, UInt8}, String}}(["Old", "Young", "Young", "Young", "Old", "Young", "Middle"])
    levels!(c4, ["Young", "Middle", "Old"])
    ds = Dataset(x1 = c1, x2 = c2, x3 = c3, x4 = c4)
    colsidx = [1, 2]
    T = nrow(ds) < typemax(Int32) ? Val(Int32) : Val(Int64)
    ranges = [1, 2, 3, 5, 6, 7]
    last_valid_index = 6
    r1, r2, r3 = IMD._find_starts_of_groups(ds, colsidx, T)
    @test r1 == colsidx
    @test r3 == last_valid_index
    @test r2[1:r3] == ranges

    # Simple functions to test large data set.
    function _zero_or_one(ds, colsidx)
        re = zeros(Bool, nrow(ds))
        re[1] = true
        re[2:end] = byrow(ds[2:end, colsidx] .!== ds[1:(end-1), colsidx], any)
        re
    end

    function _get_starts(ds, colsidx, ::Val{T}) where T
        zero_or_one = _zero_or_one(ds, colsidx)
        cols = IMD.index(ds)[colsidx]
        last_valid_index = 1
        ranges = Vector{T}(undef, length(zero_or_one))
        @inbounds for i in 1:length(zero_or_one)
            if zero_or_one[i] == true
                ranges[last_valid_index] = i
                last_valid_index += 1
            end
        end
        return cols, ranges, (last_valid_index - 1)
    end

    # Test for large data set with few levels.
    c1 = rand(1:3, 10^6)
    c2 = PooledArray(rand([missing, 1.1, 20000.0, 123.0], 10^6))
    c3 = PooledArray(rand([missing, 1.1, 20000.0, 123.0], 10^6))
    c4 = PooledArray(rand([missing, 1.1, 20000.0, 123.0], 10^6))
    c5 = rand(1:8, 10^6)
    ds = Dataset(x1 = c1, x2 = c2, x3 = c3, x4 = c4, x5 = c5)
    colsidx = [1, 2, 4, 5]
    T = nrow(ds) < typemax(Int32) ? Val(Int32) : Val(Int64)
    ra1, ra2, ra3 = IMD._find_starts_of_groups(ds, colsidx, T)
    rb1, rb2, rb3 = _get_starts(ds, colsidx, T)
    @test ra1 == rb1
    @test ra3 == rb3
    @test ra2[1:ra3] == rb2[1:rb3]

    # Test for large data set with many levels.
    c1 = rand(1:3, 10^6)
    c2 = rand(Date(100):Day(1):Date(101), 10^6)
    c3 = PooledArray(rand([missing, 1.1, 20000.0, 123.0], 10^6))
    c4 = [randstring(2) for _ in 1:10^6]
    c5 = rand(1:2, 10^6)
    ds = Dataset(x1 = c1, x2 = c2, x3 = c3, x4 = c4, x5 = c5)
    colsidx = [1, 2, 4, 5]
    T = nrow(ds) < typemax(Int32) ? Val(Int32) : Val(Int64)
    ra1, ra2, ra3 = IMD._find_starts_of_groups(ds, colsidx, T)
    rb1, rb2, rb3 = _get_starts(ds, colsidx, T)
    @test ra1 == rb1
    @test ra3 == rb3
    @test ra2[1:ra3] == rb2[1:rb3]

end

@testset "bigInt test set" begin
    ds = Dataset(x = big.([1,4,-1,1,100]), x2 = [45,3,98,100,10])
    @test sortperm(ds, 1) == [3,1,4,2,5]
    @test sortperm(ds, 1, rev = true) == [5,2,1,4,3]
    @test sortperm(ds, 1:2) == [3, 1, 4, 2, 5]
    @test sortperm(ds, 1:2, rev = [false, true]) == [3,4,1,2,5]
    ds[2,1]=missing
    @test sortperm(ds, 1) == [3,1,4,5,2]
    @test sortperm(ds, 1, rev = true) == [2,5,1,4,3]

    x = rand(big(1):big(100), 10000)
    y = rand(1:200, 10000)
    ds = Dataset(x = x, y = y)
    @test sortperm(ds, 1) == sortperm(x)
    @test sortperm(ds, 1, rev=true) == sortperm(x, rev = true)

    x = rand(Int128, 1000)
    y = rand(1:100, 1000)
    ds = Dataset(x = x, y = y)
    @test sortperm(ds, 1) == sortperm(x)
    @test sortperm(ds, 1, rev = true) == sortperm(x, rev = true)
    setformat!(ds, 1=>isodd)
    @test sortperm(ds, 1) == sortperm(x, by = isodd)
    @test sortperm(ds, 1, rev = true) == sortperm(x, by = isodd, rev = true)

    ds = Dataset(x = [0xfffffffffffffff3, 0xfffffffffffffff2, 0xfffffffffffffff4, 0xfffffffffffffff1], y = [1,1,2,2])
    @test sort(ds,1) == ds[[4,2,1,3],:]
    @test sort(ds,1, rev = true) == ds[[3,1,2,4],:]
    setformat!(ds, 1=>isodd)
    @test sort(ds,1) == ds[[2,3,1,4],:]
    @test sort(ds,1, rev=true) == ds[[1,4,2,3],:]
    @test sort(ds, 1:2) == ds[[2,3,1,4], :]
    @test sort(ds, 1:2, rev = [false, true]) == ds[[3,2,4,1], :]
end

@testset "sort views" begin
    ds = Dataset(x = [missing, 1.0, 5.0, 3.2, missing], y = [2, 1, 7, 4, 1])
    sds = view(ds, [2,1,2,2,2,1,3,4,5], [2,1])
    @test sort(ds, 1) == sort(ds[!, [1,2]], 1)
    @test sort(ds, 1) == sort(ds[!, [2,1]], 2)[:, [2,1]]
    @test sortperm(sds, 2) == [1,3,4,5,8, 7, 2,6,9]
    @test sortperm(sds, [2,1]) == [1,3,4,5,8,7, 9, 2,6]
    if !Base.Sys.iswindows()
        for i = 1:100
            ds = Dataset(rand(1:10, 1000, 3), :auto)
            @test sort(ds, :) == sort(ds[!, 1:3], :)
            ds = Dataset(rand(1:1000000, 1000, 3), :auto)
            @test sort(ds, :) == sort(ds[!, 1:3], :)
            ds = Dataset(rand(1000, 2), :auto)
            @test sort(ds, :) == sort(ds[!, 1:2], :)
        end
    end
end

@testset "issorted/issorted!" begin
    dv1 = [9, 1, 8, missing, 3, 3, 7, missing]
    dv2 = [9, 1, 8, missing, 3, 3, 7, missing]
    dv3 = Vector{Union{Int, Missing}}(1:8)
    cv1 = CategoricalArray(dv1, ordered=true)

    d = Dataset(dv1=dv1, dv2=dv2, dv3=dv3, cv1=cv1)

    @test !issorted(d, :cv1)
    @test issorted(d, :dv3)
    @test !issorted(d, :dv1)

    dv1 = [1,3,3,7,8,9, missing, missing]
    dv2 = [9, 1, 8, missing, 3, 3, 7, missing]
    dv3 = Vector{Union{Int, Missing}}(1:8)
    cv1 = CategoricalArray(dv1, ordered=true)

    d = Dataset(dv1=dv1, dv2=dv2, dv3=dv3, cv1=cv1)
    @test issorted(d, :cv1)
    @test issorted(d, :dv1)
    @test !issorted(d, :dv2)

    ds = Dataset(x = [0xfffffffffffffff3, 0xfffffffffffffff2, 0xfffffffffffffff4, 0xfffffffffffffff1], y = [1,1,2,2])
    @test issorted(ds[[4,2,1,3],:],1)
    @test issorted(view(ds, [4,2,1,3], :), 1)
    @test issorted(ds[[3,1,2,4],:],1, rev = true)
    setformat!(ds, 1=>isodd)
    @test issorted(ds[[2,3,1,4],:],1)
    @test issorted(view(ds, [2,3,1,4], :), 1)
    @test issorted(ds[[1,4,2,3],:],1, rev=true)
    @test issorted(ds[[2,3,1,4], :], 1:2)
    @test issorted(view(ds, [2,3,1,4], :), 1:2)
    @test issorted(ds[[3,2,4,1], :], 1:2, rev = [false, true])
    @test issorted(view(ds, [3,2,4,1], :), 1:2, rev = [false, true])

    x = rand(Int128, 1000)
    y = rand(1:100, 1000)
    ds = Dataset(x = x, y = y)
    @test issorted(sort(ds, 1),1)
    @test issorted(sort(ds, 1, rev = true), 1, rev=true)
    setformat!(ds, 1=>isodd)
    @test issorted(sort(ds, 1),1)
    @test issorted(sort(ds, 1, rev = true), 1, rev = true)

    ds = Dataset(x = big.([1,4,-1,1,100]), x2 = [45,3,98,100,10])
    @test !issorted(ds, 1)
    @test issorted(ds[[3,1,4,2,5], 1:1], 1)
    @test issorted(view(ds, [5,2,1,4,3], [2,1]), 2, rev = true)
    @test issorted(ds[[3, 1, 4, 2, 5], :], 1:2)
    @test issorted(ds[[3,4,1,2,5],:], 1:2, rev = [false, true])
    ds[2,1]=missing
    @test !issorted(ds, 1)
    @test issorted(ds[[3,1,4,5,2], :], 1)
    @test issorted(view(ds, [2,5,1,4,3], :), 1, rev = true)
    if !Base.Sys.iswindows()
        for i in 1:100
            ds = Dataset(rand(1:10, 1000, 10), :auto)
            for j in 1:10
                @test issorted(sort(ds, 1:j), 1:j)
                @test issorted(sort(ds, 1:j, rev = true), 1:j, rev = true)
                setformat!(ds, 1:10=>isodd)
                @test issorted(sort(ds, 1:j), 1:j)
                @test issorted(sort(ds, 1:j, rev = true), 1:j, rev = true)
            end
            ds = Dataset(rand(1:10., 1000, 10), :auto)
            map!(ds, x->rand()<.1 ? missing : x, :)
            for j in 1:10
                @test issorted(sort(ds, 1:j), 1:j)
                @test issorted(sort(ds, 1:j, rev = true), 1:j, rev = true)
                setformat!(ds, 1:10=>sign)
                @test issorted(sort(ds, 1:j), 1:j)
                @test issorted(sort(ds, 1:j, rev = true), 1:j, rev = true)
            end
            ds = Dataset(rand(1:10., 1000, 10), :auto)
            map!(ds, x->rand()<.1 ? missing : x, :)
            for j in 1:10
                ds[!, j] = PooledArray(ds[!, j])
            end
            for j in 1:10
                @test issorted(sort(ds, 1:j), 1:j)
                @test issorted(sort(ds, 1:j, rev = true), 1:j, rev = true)
                setformat!(ds, 1:10=>sign)
                @test issorted(sort(ds, 1:j), 1:j)
                @test issorted(sort(ds, 1:j, rev = true), 1:j, rev = true)
            end
        end
        for i in 1:100
            ds = Dataset(rand(1:10, 1000, 10), :auto)
            for j in 1:10
                sort!(ds, 1:j)
                issorted!(ds, 1:j)
                @test IMD._sortedcols(ds) == 1:j
                @test issorted(ds, 1:j)

                setformat!(ds, 1:10=>isodd)
                sort!(ds, 1:j, rev = true)
                issorted!(ds, 1:j, rev = true)
                @test IMD._sortedcols(ds) == 1:j
                @test issorted(ds, 1:j, rev = true)
            end
            ds = Dataset(rand(1:10., 1000, 10), :auto)
            map!(ds, x->rand()<.1 ? missing : x, :)
            for j in 1:10
                sort!(ds, 1:2:j)
                issorted!(ds, 1:2:j)
                @test IMD._sortedcols(ds) == collect(1:2:j)
                @test issorted(ds, 1:2:j)

                setformat!(ds, 1:10=>sign)
                sort!(ds, 1:2:j, rev = true)
                issorted!(ds, 1:2:j, rev = true)
                @test IMD._sortedcols(ds) == collect(1:2:j)
                @test issorted(ds, 1:2:j, rev = true)
            end
            ds = Dataset(rand(1:10., 1000, 10), :auto)
            map!(ds, x->rand()<.1 ? missing : x, :)
            for j in 1:10
                ds[!, j] = PooledArray(ds[!, j])
            end
            for j in 1:10
                sort!(ds, 1:2:j)
                issorted!(ds, 1:2:j)
                @test IMD._sortedcols(ds) == collect(1:2:j)
                @test issorted(ds, 1:2:j)

                setformat!(ds, 1:10=>sign)
                sort!(ds, 1:2:j, rev = true)
                issorted!(ds, 1:2:j, rev = true)
                @test IMD._sortedcols(ds) == collect(1:2:j)
                @test issorted(ds, 1:2:j, rev = true)
            end
        end
    end

end


text:
using Documenter, HoloProcessing
include("../contrib/LaTeXWriter.jl")

const PAGES = [
    "主页" => "index.md",
    "手册" => [
        "manual/processing.md",
        "manual/denoising.md",
        "manual/evaluating.md",
        "manual/function.md",
    ],
]

const render_pdf = "pdf" in ARGS

const format = if render_pdf
    Documenter.LaTeX(platform = "texplatform=docker" in ARGS ? "docker" : "native")
else
    Documenter.HTML(prettyurls=get(ENV, "CI", nothing) == "true")
end

makedocs(
    clean     = true,
    doctest   = ("doctest=fix" in ARGS) ? (:fix) : ("doctest=true" in ARGS) ? true : false,
    linkcheck = "linkcheck=true" in ARGS,
    checkdocs = :none,
    format    = format,
    sitename  = "HoloProcessing中文文档",
    authors   = "Qling",
    pages     = PAGES,
)

deploydocs(
    repo = "github.com/rfhklwt/HoloProcessing.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
    branch = render_pdf ? "pdf" : "gh-pages"
)


text:
# ---
# title: 238. Product of Array Except Self
# id: problem238
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/product-of-array-except-self/description/>
# hidden: true
# ---
#
# Given an array `nums` of _n_ integers where _n_ > 1,  return an array `output`
# such that `output[i]` is equal to the product of all the elements of `nums`
# except `nums[i]`.
#
# **Example:**
#
#
#
#     Input:  [1,2,3,4]
#     Output: [24,12,8,6]
#
#
# **Constraint:**  It's guaranteed that the product of the elements of any
# prefix or suffix of the array (including the whole array) fits in a 32 bit
# integer.
#
# **Note:** Please solve it **without division** and in O( _n_ ).
#
# **Follow up:**
# Could you solve it with constant space complexity? (The output array **does
# not** count as extra space for the purpose of space complexity analysis.)
#
#
## @lc code=start
using LeetCode

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


text:
function intersects_geometry(m::AbstractAccretionGeometry{T}, line_element) where {T}
    if in_nearby_region(m, line_element)
        return has_intersect(m, line_element)
    end
    false
end

in_nearby_region(d::AbstractAccretionGeometry{T}, line_element) where {T} =
    error("Not implemented for $(typeof(d))")
has_intersect(d::AbstractAccretionGeometry{T}, line_element) where {T} =
    error("Not implemented for $(typeof(d))")

# Jiménez, Segura, Feito. Computation Geometry 43 (2010) 474-492
function jsr_algorithm(V₁::T, V₂::T, V₃::T, Q₁::V, Q₂::V; ϵ = 1e-8) where {T,V}
    A = Q₁ .- V₃
    B = V₁ .- V₃
    C = V₂ .- V₃
    W₁ = B × C
    w = A ⋅ W₁
    if w > ϵ
        D = Q₂ .- V₃
        s = D ⋅ W₁
        s > ϵ && return false
        W₂ = A × D
        t = W₂ ⋅ C
        t < -ϵ && return false
        u = -W₂ ⋅ B
        u < -ϵ && return false
        w < s + t + u && return false
    elseif w < -ϵ
        return false
    else # w == 0
        D = Q₂ .- V₃
        s = D ⋅ W₁
        if s > ϵ
            return false
        elseif s < -ϵ
            W₂ = D × A
            t = W₂ ⋅ C
            t > ϵ && return false
            u = -W₂ ⋅ B
            u > ϵ && return false
            -s > t + u && return false
        else
            return false
        end
    end
    true
end


text:
function get_components_to_pm(sys::PSY.System, ::Type{T}) where {T <: PSY.Area}
    areas = PSY.get_components(PSY.Area, sys)

    PM_areas = Dict{String, Any}()

    for (ix, area) in enumerate(areas)
        buses = PSY.get_buses(sys, area)
        ref_bus = [b for b in buses if PSY.get_bustype(b) == PSY.BusTypes.REF]
        if length(ref_bus) == 1
            ref_bus = ref_bus[1]
        elseif PSY.BusTypes.PV ∈ PSY.get_bustype.(buses)
            # get largest gen bus
            pv_buses = [b for b in buses if PSY.get_bustype(b) == PSY.BusTypes.PV]
            gens = PSY.get_components(PSY.Generator, sys, x -> PSY.get_bus(x) ∈ pv_buses)
            ref_bus = PSY.get_bus(
                last(sort(collect(gens), by = x -> PSY.get_max_active_power(x))),
            )
        else
            ref_bus = buses[1]
        end

        PM_areas["$(ix)"] = Dict{String, Any}(
            "source_id" => ["areas", ix],
            "col_1" => ix,
            "col_2" => PSY.get_number(ref_bus),
            "index" => ix,
        )
    end
    return PM_areas
end

max_stars_repo_name: NREL-SIIP/PowerModelsInterface

max_issues_repo_name: NREL-SIIP/PowerModelsInterface

max_forks_repo_name: NREL-SIIP/PowerModelsInterface


text:
export Wigner3j!
@doc raw"""
    w3j, jmin, jmax = Wigner3j!(w3j::AbstractVector{Cdouble}, j2::Integer, j3::Integer,
                                m1::Integer, m2::Integer, m3::Integer,
                                exitstatus::Union{Nothing,Ref{<:Integer}} = nothing)

Compute the Wigner 3j symbol

```math
\left(\begin{array}{ccc}
j & j_{2} & j_{3}\\
m_{1} & m_{2} & m_{3}
\end{array}\right)
```

for all integral values of `j` that satisfy `max(abs(m1), abs(j2 - j3)) <= j <= j2 + j3`.
The values of `jmin` and `jmax` are set to `max(abs(m1), abs(j2 - j3))` and `j2 + j3` respectively.
The vector `w3j` is filled with the computed 3j symbols.
Its first index corresponds to `j = jmin`,
and it has its first `jmax - jmin + 1` indices populated.

!!! note
    `w3j` must be sufficiently long to contain all the elements.
    Pre-existing values in `w3j` may be overwritten.

See also: [`Wigner3j`](@ref)

SHTOOLS page: [`Wigner3j`](https://shtools.oca.eu/shtools/public/wigner3j.html)
"""
function Wigner3j!(w3j::AbstractVector{Cdouble},
    j2::Integer, j3::Integer,
    m1::Integer, m2::Integer, m3::Integer,
    exitstatus::Optional{Ref{<:Integer}}=nothing)

    @assert j2 ≥ 0 "j2 must be >= 0"
    @assert j3 ≥ 0 "j2 must be >= 0"
    @assert abs(m2) ≤ j2 "m2 must satisfy $(-j2) <= m2 <= $j2 for j2 = $j2"
    @assert abs(m3) ≤ j3 "m3 must satisfy $(-j3) <= m3 <= $j3 for j3 = $j3"
    @assert abs(m1) ≤ j2 + j3 "m1 must satisfy $(-(j2 + j3)) <= m1 <= $(j2 + j3)"
    w3j_minlength = j2 + j3 - max(abs(m1), abs(j2 - j3)) + 1
    @assert length(w3j) ≥ w3j_minlength "length of w3j must be >= $w3j_minlength"

    exitstatus′ = Ref{Cint}()
    jmin = Ref{Cint}()
    jmax = Ref{Cint}()
    ccall((:Wigner3j, libSHTOOLS), Cvoid,
          (Ptr{Cdouble}, Ref{Cint}, Ref{Cint},
            Cint, Cint, Cint, Cint, Cint, Ref{Cint}),
          w3j, jmin, jmax, j2, j3, m1, m2, m3, exitstatus′)
    if exitstatus === nothing
        exitstatus′[] ≠ 0 &&
            error("Wigner3j!: Error code $(exitstatus′[])")
    else
        exitstatus[] = exitstatus′[]
    end

    return w3j, jmin[], jmax[]
end

export Wigner3j
@doc raw"""
    w3j, jmin, jmax = Wigner3j(j2::Integer, j3::Integer,
                                m1::Integer, m2::Integer, m3::Integer,
                                exitstatus::Union{Nothing,Ref{<:Integer}} = nothing)

Compute the Wigner 3j symbol

```math
\left(\begin{array}{ccc}
j & j_{2} & j_{3}\\
m_{1} & m_{2} & m_{3}
\end{array}\right)
```

for all integral values of `j` that satisfy `max(abs(m1), abs(j2 - j3)) <= j <= j2 + j3`.
The values of `jmin` and `jmax` are set to `max(abs(m1), abs(j2 - j3))` and `j2 + j3` respectively.
The vector `w3j` contains the computed 3j symbols. Its first index corresponds to `j = jmin`,
and last index corresponds to `j = jmax`.

See also: [`Wigner3j!`](@ref)

SHTOOLS page: [`Wigner3j`](https://shtools.oca.eu/shtools/public/wigner3j.html)
"""
function Wigner3j(j2::Integer, j3::Integer,
    m1::Integer, m2::Integer, m3::Integer, exitstatus::Optional{Ref{<:Integer}}=nothing)

    @assert j2 ≥ 0 "j2 must be >= 0"
    @assert j3 ≥ 0 "j2 must be >= 0"
    @assert abs(m1) ≤ j2 + j3 "m1 must satisfy $(-(j2 + j3)) <= m1 <= $(j2 + j3)"

    w3j_minlength = j2 + j3 - max(abs(m1), abs(j2 - j3)) + 1
    w3j = Vector{Cdouble}(undef, w3j_minlength)
    w3j, jmin, jmax = Wigner3j!(w3j, j2, j3, m1, m2, m3, exitstatus)
    return w3j, jmin, jmax
end


text:
using Redis: RedisConnection

function getdatabase(;host::AbstractString="127.0.0.1", port::Integer=6379, password::AbstractString="", db::Integer=0)
    db_conn = RedisConnection(host=host, port=port, password=password, db=db)
    return db_conn
end


text:
@testset verbose=true "Quaternion" begin
    @test eltype([1.0, imx]) === QuaternionF64
    @test [1.0, imx] == QuaternionF64[QuaternionF64(1.0, 0.0, 0.0, 0.0), QuaternionF64(0.0, 1.0, 0.0, 0.0)]

    @testset "$Q{T}" for Q in [Quaternion, Rotor, QuatVec]
        for T in Types
            @test Q(T) === Q{T}
            @test Q(Q{T}) === Q{T}
            if Q === Rotor
                @test real(Q{T}) === T
            else
                @test real(Q{T}) === T
            end
        end

        for TypeGroup in [FloatTypes, IntTypes]
            for i1 in 1:length(TypeGroup)
                T1 = TypeGroup[i1]
                for i2 in i1+1:length(TypeGroup)
                    T2 = TypeGroup[i2]
                    if all(Ti ∈ [BigFloat, Symbolics.Num] for Ti ∈ [T1, T2])
                        continue
                    end
                    @test promote_rule(T1, T2) === T1
                    if Q === QuatVec
                        @test promote_rule(Q{T1}, T2) === Quaternion{T1}
                    else
                        @test promote_rule(Q{T1}, T2) === Q{T1}
                    end
                    @test promote_rule(Q{T1}, Q{T2}) === Q{T1}
                end
            end
        end

        for T1 in Types
            for T2 in Types
                if all(Ti ∈ [BigFloat, Symbolics.Num] for Ti ∈ [T1, T2])
                    continue
                end
                for Q1 in QTypes
                    for Q2 in QTypes
                        Qexpected = Q1===Q2 ? Q1 : Quaternion
                        Texpected = promote_type(T1, T2)
                        expected = Qexpected{Texpected}
                        observed = promote_rule(Q1{T1}, Q2{T2})
                        @test expected === observed
                    end
                end
            end
        end

        for TypeGroup in [FloatTypes, IntTypes]
            @test widen(Q{TypeGroup[1]}) === Q{TypeGroup[1]}
            for i in 2:length(TypeGroup)
                T1 = TypeGroup[i-1]
                T2 = TypeGroup[i]
                @test widen(Q{T2}) === Q{T1}
            end
        end

        for T in FloatTypes
            @test float(Q{T}) === Q{T}
            @test big(Q{T}) === Q{big(T)}
            if Q !== Rotor
                @test float(Q{T}(1, 2, 3, 4)) == Q(float(T)(1), float(T)(2), float(T)(3), float(T)(4))
                @test big(Q{T}(1, 2, 3, 4)) == Q(big(T)(1), big(T)(2), big(T)(3), big(T)(4))
                @test big(Q{T}(1, 2, 3, 4)) == Q{big(T)}(1, 2, 3, 4)
                @test typeof(Q{T}(1, 2, 3, 4) * rand(T)) === Q{T}
                @test typeof(Q{T}(1, 2, 3, 4) / rand(T)) === Q{T}
                @test typeof(rand(T) * Q{T}(1, 2, 3, 4)) === Q{T}
                @test typeof(rand(T) / Q{T}(1, 2, 3, 4)) === Q{T}
            end
            @test float(Q{T}(1, 2, 3, 4)) == Q{T}(1, 2, 3, 4)
        end

        for T in IntTypes
            @test float(Q{T}) === Q{float(T)}
            if Q !== Rotor
                @test float(Q{T}(1, 2, 3, 4)) == Q(float(T)(1), float(T)(2), float(T)(3), float(T)(4))
            end
            for i in 1:4
                v = zeros(T, 4)
                v[i] = one(T)
                v = SVector(v...)
                @test float(Q{T}(v)) == Q(float(T).(v))
                @test float(Q{T}(v)) == Q{T}(float(T).(v))
                v = Vector(v)
                @test float(Q{T}(v)) == Q{T}(float(T).(v))
            end
        end

        for T in SymbolicTypes
            q = Q{T}(1, 2, 3, 4)
            if Q !== Rotor
                @test typeof(a * q) === Q{T}
                @test typeof(a / q) === Q{T}
                @test typeof(q * a) === Q{T}
                @test typeof(q / a) === Q{T}
            end
        end
    end

end


text:
"""
This function adds the power limits of generators when there are no CommitmentVariables
"""
function activepower(m::JuMP.AbstractModel, devices::Array{T,1}, device_formulation::Type{D}, system_formulation::Type{S}, time_periods::Int64) where {T <: PowerSystems.ThermalGen, D <: AbstractThermalDispatchForm, S <: PM.AbstractPowerFormulation}

    p_th = m[:p_th]
    time_index = m[:p_th].axes[2]
    name_index = m[:p_th].axes[1]

    (length(time_index) != time_periods) ? @error("Length of time dimension inconsistent") : true

    pmax_th = JuMP.JuMPArray(Array{ConstraintRef}(undef,length(name_index), time_periods), name_index, time_index)
    pmin_th = JuMP.JuMPArray(Array{ConstraintRef}(undef,length(name_index), time_periods), name_index, time_index)

    for t in time_index, (ix, name) in enumerate(name_index)

        if name == devices[ix].name

            pmin_th[name, t] = @constraint(m, p_th[name, t] >= devices[ix].tech.activepowerlimits.min)
            pmax_th[name, t] = @constraint(m, p_th[name, t] <= devices[ix].tech.activepowerlimits.max)

        else
            @error "Bus name in Array and variable do not match"
        end

    end

    JuMP.register_object(m, :pmax_th, pmax_th)
    JuMP.register_object(m, :pmin_th, pmin_th)

    return m
end

"""
This function adds the power limits of generators when there are CommitmentVariables
"""
function activepower(m::JuMP.AbstractModel, devices::Array{T,1}, device_formulation::Type{D}, system_formulation::Type{S}, time_periods::Int64) where {T <: PowerSystems.ThermalGen, D <: AbstractThermalCommitmentForm, S <: AbstractDCPowerModel}

    p_th = m[:p_th]
    on_th = m[:on_th]

    time_index = m[:p_th].axes[2]
    name_index = m[:p_th].axes[1]

    (length(time_index) != time_periods) ? @error("Length of time dimension inconsistent") : true

    pmax_th = JuMP.JuMPArray(Array{ConstraintRef}(undef,length(name_index), time_periods), name_index, time_index)
    pmin_th = JuMP.JuMPArray(Array{ConstraintRef}(undef,length(name_index), time_periods), name_index, time_index)

    for t in time_index, (ix, name) in enumerate(name_index)

        if name == devices[ix].name
            pmin_th[name, t] = @constraint(m, p_th[name, t] >= devices[ix].tech.activepowerlimits.min*on_th[name,t])
            pmax_th[name, t] = @constraint(m, p_th[name, t] <= devices[ix].tech.activepowerlimits.max*on_th[name,t])
        else
            @error "Bus name in Array and variable do not match"
        end

    end

    JuMP.register_object(m, :pmax_th, pmax_th)
    JuMP.register_object(m, :pmin_th, pmin_th)

    return m
end

"""
This function adds the power limits of generators when there are no CommitmentVariables
"""
function reactivepower(m::JuMP.AbstractModel, devices::Array{T,1}, device_formulation::Type{D}, system_formulation::Type{S}, time_periods::Int64) where {T <: PowerSystems.ThermalGen, D <: AbstractThermalDispatchForm, S <: AbstractACPowerModel}

    q_th = m[:q_th]
    time_index = m[:q_th].axes[2]
    name_index = m[:q_th].axes[1]

    (length(time_index) != time_periods) ? @error("Length of time dimension inconsistent") : true

    qmax_th = JuMP.JuMPArray(Array{ConstraintRef}(undef,length(name_index), time_periods), name_index, time_index)
    qmin_th = JuMP.JuMPArray(Array{ConstraintRef}(undef,length(name_index), time_periods), name_index, time_index)

    for t in time_index, (ix, name) in enumerate(name_index)

        if name == devices[ix].name

            qmin_th[name, t] = @constraint(m, q_th[name, t] >= devices[ix].tech.reactivepowerlimits.min)
            qmax_th[name, t] = @constraint(m, q_th[name, t] <= devices[ix].tech.reactivepowerlimits.max)

        else
            @error "Bus name in Array and variable do not match"
        end

    end

    JuMP.register_object(m, :qmax_th, qmax_th)
    JuMP.register_object(m, :qmin_th, qmin_th)

    return m
end

"""
This function adds the power limits of generators when there are CommitmentVariables
"""
function reactivepower(m::JuMP.AbstractModel, devices::Array{T,1}, device_formulation::Type{D}, system_formulation::Type{S}, time_periods::Int64) where {T <: PowerSystems.ThermalGen, D <: AbstractThermalCommitmentForm, S <: AbstractACPowerModel}

    q_th = m[:p_th]
    on_th = m[:on_th]

    time_index = m[:q_th].axes[2]
    name_index = m[:q_th].axes[1]

    (length(time_index) != time_periods) ? @error("Length of time dimension inconsistent") : true

    qmax_th = JuMP.JuMPArray(Array{ConstraintRef}(undef,length(name_index), time_periods), name_index, time_index)
    qmin_th = JuMP.JuMPArray(Array{ConstraintRef}(undef,llength(name_index), time_periods), name_index, time_index)

    for t in time_index, (ix, name) in enumerate(name_index)

        if name == devices[ix].name
            qmin_th[name, t] = @constraint(m, q_th[name, t] >= devices[ix].tech.reactivepowerlimits.min*on_th[name,t])
            qmax_th[name, t] = @constraint(m, q_th[name, t] <= devices[ix].tech.reactivepowerlimits.max*on_th[name,t])
        else
            @error "Bus name in Array and variable do not match"
        end

    end

    JuMP.register_object(m, :qmax_th, qmax_th)
    JuMP.register_object(m, :qmin_th, qmin_th)

    return m
end


text:
"""
    AbstractMetaData

Thin wrapper for data to store cache and 'mtime' for the data file
"""
abstract type AbstractMetaData end

Base.getindex(x::AbstractMetaData, i) = getindex(x.data, i)

filepath(x::AbstractMetaData) = x.filepath

function isloaded(x::AbstractMetaData)
    !ismissing(x.data)
end
function ismodified(x::AbstractMetaData)
    if isloaded(x)
        return x.mtime != mtime(filepath(x))
    else
        return true
    end
end
function update!(x::AbstractMetaData)
    if ismodified(x)
        loaddata!(x)
    end
    return x
end

mutable struct SchemaData <: AbstractMetaData
    filepath::AbstractString
    data::Union{Missing, JSONSchema.Schema}
    mtime::Float64
end
function SchemaData(file)
    SchemaData(file, missing, 0.)
end
function loaddata!(x::SchemaData)
    x.data = JSONSchema.Schema(JSON.parsefile(x.filepath; use_mmap=false))
    x.mtime = mtime(x.filepath)
    return x
end


text:
#################################################################
# The policy is a mapping from a belief to an action in a POMDP,
# and it maps states to actions in an MDP.
# The policy is extracted through calls to the action() function.
#################################################################

abstract Policy

# creates an action object (for preallocation purposes)
@pomdp_func create_action(pomdp::POMDP)

# returns a default BeliefUpdater appropriate for a belief type that policy `p` can use
@pomdp_func action(p::Policy, s::State, a::Action)
@pomdp_func action(p::Policy, s::State)
@pomdp_func action(p::Policy, b::Belief, a::Action)
@pomdp_func action(p::Policy, b::Belief)

# returns a default BeliefUpdater appropriate for a belief type that policy `p` can use
@pomdp_func updater(policy::Policy)

# returns the utility value from policy p given the belief
@pomdp_func value(p::Policy, belief::Belief)
# returns the utility value from policy p given the state
@pomdp_func value(p::Policy, state::State)

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
export MatchesIterator, ParseMatch

"""
Iterator type for `match_all` and `parse_all`.
"""
@auto_hash_equals struct MatchesIterator{P,S}
    parser::P
    sequence::S
    from::Int
    till::Int
    MatchesIterator(parser::P,sequence::S,from=1,till=lastindex(sequence)) where {P,S} =
        new{P,S}(parser,sequence,from,till)
end
result_type(::Type{<:MatchesIterator{P}}) where P = result_type(P)
Base.eltype(T::Type{<:MatchesIterator{P,S}}) where {P,S} =
    ParseMatch{P,S,state_type(P)}
Base.IteratorSize(::Type{<:MatchesIterator}) = Base.SizeUnknown()

export ParseMatch
"""
Wrapper type for [`CombinedParsers.Regexp.SequenceWithCaptures`](@ref), providing
`getindex` and `getproperty` behavior like `RegexMatch`.

```jldoctest
julia> m = match(re"(?<a>so)+ (or)", "soso or")
ParseMatch("soso or", a="so", 2="or")

julia> m[:a]
"so"

julia> m[2]
"or"

julia> m.match, m.captures
("soso or", SubString{String}["so", "or"])

```
"""
@auto_hash_equals struct ParseMatch{P,S,State}
    parsings::MatchesIterator{P,S}
    start::Int
    stop::Int
    state::State
    function ParseMatch(p::MatchesIterator{P,S},start::Integer,stop::Integer,state::ST) where {P,S,ST}
        new{P,S,ST}(
            p,
            start,stop,state)
    end
end
function ParseMatch(parser,s,start=1,stop=start,state=nothing)
    ParseMatch(MatchesIterator(parser,s),
               start,stop,state)
end
result_type(::Type{<:ParseMatch{P}}) where P = result_type(P)
Base.eltype(T::Type{<:ParseMatch}) = ParseMatch
Base.IteratorSize(::Type{<:ParseMatch}) = Base.SizeUnknown()
_getproperty(x::ParseMatch,key) =
    if key == :parsings
        getfield(x,1)
    elseif key == :start
        getfield(x,2)
    elseif key == :stop
        getfield(x,3)
    elseif key == :state
        getfield(x,4)
    else
        error("no property $key")
    end

function Base.show(io::IO,m::ParseMatch{<:Any,<:AbstractString,<:Any})
    print(io,"ParseMatch(\"",escape_string(m.match),"\"")
    print(io,")")
end

function Base.getproperty(x::ParseMatch{<:Any,<:AbstractString,<:Any},key::Symbol)
    if key==:captures
        AbstractString[ ]
    elseif key==:match
        SubString(x.parsings.sequence,
                  x.start,
                  _prevind(x.parsings.sequence,x.stop))
    else
        _getproperty(x,key)
    end
end

function Base.show(io::IO,m::ParseMatch{<:Any,<:AbstractString,Nothing})
    print(io,"no match")
end
function Base.show(io::IO,m::ParseMatch{<:Any,<:Any,Nothing})
    print(io,"no match")
end

Base.getproperty(x::ParseMatch{<:Any,<:AbstractString,Nothing},key::Symbol) =
    if key==:captures
        AbstractString[ ]
    elseif key==:match
        ""
    else
        _getproperty(x,key)
    end

"""
    parsematch_tuple(m,start,state)

ParseMatch iteration has the first match as iterator, the last match as a state.
(Turned out to be fastest.)
"""
parsematch_tuple(m,start,state) =
    let r = ParseMatch(m,start,tuple_pos(state),tuple_state(state))
        return tuple(r,r)
    end
parsematch_tuple(m,start,state::Nothing) = nothing

"""
    iterate(m::CombinedParser,s::AbstractString)

Returns `iterate(ParseMatch(m,s,1,1,nothing))`.

See also [`ParseMatch`](@ref).
"""
Base.iterate(m::CombinedParser,s::AbstractString) =
    iterate(ParseMatch(m,s,1,1,nothing))
import Base: iterate
"""
    Base.iterate(x::ParseMatch[, m::ParseMatch=x])

Returns next match at `m.start` after `m.state`, see [`_iterate`](@ref)(m).

See also [`ParseMatch`](@ref).
"""
function Base.iterate(x::ParseMatch, m=x)
    i = _iterate(m)
    parsematch_tuple(m.parsings,m.start,i)
end
_iterate(m::ParseMatch) =
    _iterate(m.parsings,m.start,m.stop,m.state)
export match_all

"""
    match_all(parser::CombinedParser, sequence, idx=1; log=nothing)

Returns an iterator over all matches at all indices in the sequence.
"""
function match_all(parser::CombinedParser, sequence, idx=1; log=nothing)
    p = (log === nothing || log == false ) ? parser : log_names(parser,log)
    MatchesIterator(p,sequence,idx)
end

import Base: iterate
"""
    Base.iterate(x::MatchesIterator[, s::ParseMatch=ParseMatch(x,x.from,x.from,nothing)])

Iterate match `s` at current position.
While no match is found and `s.start<=x.till`, increase s.start.
Return first next [`ParseMatch`](@ref) (as return value and state) or `nothing` when at `m.till`.
"""
@inline Base.iterate(x::MatchesIterator) =
    iterate(x,ParseMatch(x,x.from,x.from,nothing))
function Base.iterate(m::MatchesIterator,
                      s::ParseMatch)
    start,stop = s.start, s.stop
    till = m.till
    state = _iterate(m,start,stop,s.state)
    while start <= till+1 && state===nothing
        # state = iterate(m.parsings,(start,nothing))
        start > till && break
        start = _nextind(m.sequence,start)
        state = _iterate(m,start,start,nothing)
    end
    parsematch_tuple(m,start,state)
end

_iterate(p::CombinedParser,s) =
    _iterate(MatchesIterator(p,s),1,1,nothing)
_iterate(mi::MatchesIterator,a...) =
    _iterate(mi.parser, mi.sequence, mi.till, a...)

"""
    Base.match(parser::CombinedParser,sequence::AbstractString[, idx::Integer]; log=nothing)

Plug-in replacement for `match(::Regex,sequence)`.

If `log!==nothing`, parser is transformed with `log_names(p, log)`.
See also [`log_names`](@ref).
"""
function Base.match(parser::CombinedParser, sequence, idx=1; log=nothing)
    p = (log === nothing || log == false ) ? parser : log_names(parser,log)
    i = iterate(MatchesIterator(p,sequence,idx))
    i === nothing && return nothing
    i[1]
end

import Base: tryparse, parse
"""
    parse(parser::CombinedParser, str::AbstractString; log=nothing)

Parse a string with a CombinedParser as an instance of `result_type(parser)`.

If `log!==nothing`, parser is transformed with `log_names(p, log)`.
See also [`log_names`](@ref).

```@meta
DocTestFilters = r"map\\(.+\\)"
```

```jldoctest
julia> using TextParse

julia> p = ("Number: "*TextParse.Numeric(Int))[2]
🗄 Sequence[2]
├─ Number\\:\\
└─ <Int64>
::Int64

julia> parse(p,"Number: 42")
42

```

"""
function Base.parse(p::AbstractToken, s; kw...)
    i = tryparse(p, s; kw...)
    i === nothing && throw(ArgumentError("no successfull parsing."))
    i
end

"""
    tryparse(parser::CombinedParser, str::AbstractString[, idx=1]; log=nothing)

Like `parse`, but returns either a value of `result_type(parser)` or `nothing` if string does not start with with a match.

If `log!==nothing`, parser is transformed with `log_names(p, log)`.
See also [`log_names`](@ref).
"""
function Base.tryparse(p::AbstractToken, s, idx=1; log=nothing)
    i = iterate(ParseMatch(log !== nothing ? log_names(p,log) : p,s,idx))
    i === nothing && return nothing
    get(tuple_state(i))
end

export tryparse_pos
"""
    tryparse_pos(parser::CombinedParser, str::AbstractString[, idx=1, till=lastindex(s)])

Like `parse`, but returns either a tuple of `result_type(parser)` and the position after the match, or `nothing` if string does not start with with a match.
"""
function tryparse_pos(p,s, idx=1, till=lastindex(s))
    i = _iterate(p,s,till,idx,idx,nothing)
    i === nothing && return nothing
    get(p,s,lastindex(s),tuple_pos(i),1,tuple_state(i)),tuple_pos(i)
end

export parse_all
"""
    parse_all(parser::CombinedParser, sequence, idx=1)

Returns an iterator over all parsings of the sequence starting at `idx`.
"""
function parse_all(parser::CombinedParser, sequence, idx=1)
    ( get(p) for p=ParseMatch(parser,sequence,idx) )
end

"""
    Base.get(x::ParseMatch{<:MatchTuple})

Get the result of a match result.

```jldoctest
julia> m = match(re"(?<a>so)+ (or)", "soso or")
ParseMatch("soso or", a="so", 2="or")

julia> get(m)
([('s', 'o'), ('s', 'o')], ' ', ('o', 'r'))

julia> m[2]
"or"

julia> m.match, m.captures
("soso or", SubString{String}["so", "or"])

```
"""
Base.get(x::ParseMatch)=
    get(x.parsings,x.stop,x.start,x.state)

Base.get(x::MatchesIterator,a...)=
    get(x.parser,x.sequence,x.till,a...)