text:
using PowerModelsDistribution
using PowerModels
using Ipopt
using JuMP
using DataFrames
using CSV
# using Memento

# perform changes to the PowerModels data structure
function  data_surgery(data,filename,filename_flex)
    # data = PowerModelsDistribution.parse_file("IEEETestCases/123Bus_GMLC/IEEE123Master.dss")
    data_new = deepcopy(data)

    ADC = build_ADC_list(filename)
    ADC_pm = convert_ADC_pm(data,ADC)

    bus_loads = build_bus_loads(data)

    # @show buspairs = branches_to_buspairs(data,branches_to_remove)

    # change voltage bounds for buses in the adc control area
    for (adc,control_area) in ADC_pm
        # control_area_trimmed  = control_area
        control_area_trimmed = setdiff(control_area,[adc])
        for b in control_area_trimmed
            # delete!(data_new["bus"],b)
            data_new["bus"][b]["vmin"] = [-100,-100,-100]
            data_new["bus"][b]["vmax"] = [100,100,100]
        end
    end

    # increase branch limits for branches in the control area
    branches_to_remove = find_branches_to_remove(data,ADC_pm)
    for b in branches_to_remove
        # delete!(data_new["branch"],b)
        data_new["branch"][b]["rate_a"] = [100,100,100]
        data_new["branch"][b]["rate_b"] = [100,100,100]
        data_new["branch"][b]["rate_c"] = [100,100,100]
    end

    # add loads
    adc_sumloads = find_ADC_sumloads(data,ADC_pm,bus_loads)

    for adc in keys(ADC_pm)
        if adc in keys(bus_loads)
            # println("position 1")
            load = bus_loads[adc][1]
            # println("type of load connection = ", data["load"][load]["conn"])
            # println("type of load model = ", data["load"][load]["model"])

            # println("position 2")
            data_new["load"][load]["pd"]  = adc_sumloads[adc][:,1]
            # println("position 3")
            data_new["load"][load]["qd"] = adc_sumloads[adc][:,2]
            # println("position 4")
            data_new["load"][load]["model"] = "constant_power"
        else
            load_dict = deepcopy(data_new["load"]["1"])
            load_dict["name"] = ""
            load_dict["conn"] = "wye"
            load_dict["model"] = "constant_power"
            load_dict["active_phases"] = [1, 2, 3]
            load_dict["load_bus"] = parse(Int64,adc)
            load_dict["pd"] = adc_sumloads[adc][:,1]
            load_dict["qd"] = adc_sumloads[adc][:,2]
            load_dict["source_id"] = ""
            load_indices = collect(parse(Int64,num) for (num,l) in data_new["load"])
            new_index = maximum(load_indices) + 1
            load_dict["index"] = new_index
            data_new["load"][string(new_index)] = load_dict

            bus_loads[adc] = [string(new_index)]
            # @show length(data_new["load"])
            println("added new load")
        end

    end

    # remove loads
    for (adc,control_area) in ADC_pm
        # control_area_trimmed = control_area
        control_area_trimmed = setdiff(control_area,[adc])
        for bus in control_area_trimmed
            if bus in keys(bus_loads)
                for load in bus_loads[bus]
                    delete!(data_new["load"],load)
                end
            end
        end
    end
    for adc in keys(ADC_pm)
        if adc in keys(bus_loads)
            loads = bus_loads[adc]
            if length(loads) > 1
                println("Indeed too many loads at an ADC")
                for i=2:length(loads)
                    load = loads[i]
                    delete!(data_new["load"],load)
                end
            end
        end
    end

    return data_new
end
# end of data surgery

# parse adc flexibility regions
function build_adc_flex_model(filename_flex, data)

    pm_to_source, source_to_pm = build_bus_renaming_dict(data)

    data_source = CSV.read(filename_flex)
    adc_flex = Dict{String,Array{Float64,2}}()
    for i=1:size(data_source,1)
        adc_name = data_source[i,1]
        pos = findfirst(isequal('C'), adc_name)
        adc_name = adc_name[pos+1:end]
        intervals = split(data_source[i,end],"|")
        pmin = parse(Float64,intervals[1][2:end])
        pmax = parse(Float64,intervals[2])
        qmin = parse(Float64,intervals[3])
        qmax = parse(Float64,intervals[4][1:end-1])

        # if parse(Int64,adc_name) == 52
        #     # custom changes to ADC at node 52
        #
        #
        #     # ADC["52"] = ["52"]
        #     # ADC["53"] = ["53"]
        #     # ADC["54"]= ["54","55","56"]
        #     # ADC["57"] = ["57","58","59"]
        #     # ADC["60"] = ["60","61"]
        # else
            # convert to pm index
        adc_name = source_to_pm[adc_name]
        adc_flex[adc_name]  = [pmin pmax; qmin qmax]
        # end

    end

    return adc_flex
end

# auxiliary functions
function build_ADC_list(filename)
    adc_raw_data = CSV.read(filename)
    adc_raw_data = adc_raw_data[2:2:end,:]

    adc_bus_numbers = unique!(collect(adc_raw_data[i,end] for i=1:size(adc_raw_data,1)))
    adc_bus_numbers = adc_bus_numbers[2:end]

    ADC = Dict{String,Array{String,1}}()

    for i in adc_bus_numbers
        ADC[i] = []
    end

    for i=1:size(adc_raw_data,1)
        if adc_raw_data[i,end]!= "NONE"
            adc_bus_num = adc_raw_data[i,end]
            bus_num = string(adc_raw_data[i,1])
            push!(ADC[adc_bus_num],bus_num)
        end
    end
    # @show ADC

    # custom changes to ADC at node 52
    # delete!(ADC,"52")
    #
    # ADC["52"] = ["52"]
    # ADC["53"] = ["53"]
    # ADC["54"]= ["54","55","56"]
    # ADC["57"] = ["57","58","59"]
    # ADC["60"] = ["60","61"]

    return ADC
end

function convert_ADC_pm(data,ADC)
    pm_to_source, source_to_pm = build_bus_renaming_dict(data)

    ADC_pm = Dict{String,Array{String,1}}()
    for (adc,control_area) in ADC
        if adc in keys(source_to_pm)
            ADC_pm[source_to_pm[adc]] = collect(source_to_pm[b] for b in control_area if b in keys(source_to_pm))
        end
    end

    return ADC_pm
end

function build_bus_renaming_dict(data)
    source_to_pm = Dict{String,String}()
    pm_to_source = Dict{String,String}()
    for (i,bus) in data["bus"]
        if "source_id" in keys(bus)
            continue
        else
            pm_to_source[i] = bus["name"]
            source_to_pm[bus["name"]] = i
        end
    end
    return pm_to_source, source_to_pm
end

function build_bus_loads(data)
    load_buses = unique!(collect(string(load["load_bus"]) for (i,load) in data["load"]))
    bus_loads = Dict{String,Array{String,1}}()
    for i in load_buses
        bus_loads[i] = []
    end

    for (i,load) in data["load"]
        push!(bus_loads[string(load["load_bus"])], i)
    end
    # @show bus_loads
    return bus_loads
end

function find_ADC_sumloads(data,ADC,bus_loads)
    adc_sumloads = Dict{String,Array{Float64,2}}()

    for (adc,control_area) in ADC
        adc_sumloads[adc] = [0 0;0 0; 0 0]
        for bus in control_area
            if bus in keys(bus_loads)
                loads = bus_loads[bus]
                for l in loads
                    adc_sumloads[adc] += [data["load"][l]["pd"] data["load"][l]["qd"]]
                end
            end
        end
    end

    return adc_sumloads
end

    # fix bus renaming issues
function find_branches_to_remove(data,ADC)
    branches_to_remove = []
    for (adc,control_area) in ADC
        control_area_trimmed = setdiff(control_area,[adc])
        # branches = []
        for (b,branch) in data["branch"]
            t_bus = data["bus"][string(branch["t_bus"])]["index"]
            f_bus = data["bus"][string(branch["f_bus"])]["index"]
            if string(branch["t_bus"]) in control_area_trimmed
                push!(branches_to_remove,b)
            elseif string(branch["t_bus"]) in control_area_trimmed
                push!(branches_to_remove,b)
            end
        end
    end
    return unique!(branches_to_remove)
end

function branches_to_buspairs(data, branches)
    # fix bus renaming issues

    buspairs = []
    for b in branches
        t_bus = data["branch"][b]["t_bus"]
        f_bus = data["branch"][b]["f_bus"]

        bp = (data["bus"][string(f_bus)]["index"],data["bus"][string(t_bus)]["index"])

        push!(buspairs,bp)
    end
    return buspairs
end

# data = PowerModelsDistribution.parse_file("IEEETestCases/123Bus_GMLC/IEEE123Master.dss")

max_stars_repo_name: areiman/PFO-ADC-DER-Testbed

max_issues_repo_name: areiman/PFO-ADC-DER-Testbed

max_forks_repo_name: areiman/PFO-ADC-DER-Testbed


text:
type CompositeModF{T<:FP} <: NLregModF{T}
    pt::Ptrans
    nm::NLregModF{T}
end

function (*)(vv::Vector{DataType},nm::NLregModF)
    dt = Dtrans(vv)
    length(dt) == npars(nm) || error("DimensionMisMatch")
    CompositeModF(dt,nm)
end

nlmod(cm::CompositeModF) = cm.nm

npars(cm::CompositeModF) = length(cm.pt)

ptran(cm::CompositeModF) = cm.pt

pnames(cm::CompositeModF) = pnames(ptran(cm), pnames(nlmod(cm)))

StatsBase.residuals(cm::CompositeModF) = residuals(cm.nm)

function updtmu!(cm::CompositeModF,pars::Vector)
    length(pars) == npars(cm) || throw(DimensionMismatch(""))
    nm = nlmod(cm)
    tpars, jac = parjac(ptran(cm), pars)
    res = updtmu!(nm,tpars)
    if isa(jac,Diagonal)                # can update tgrad in place
        A_mul_B!(jac,tgrad(nm))
    else
        updttg = jac * (tg = tgrad(nm))
        tg[:] = updttg
    end
    res
end

type CompositePLModF{T<:FP} <: PLregModF{T}
    pt::Ptrans
    plm::PLregModF{T}
end

function (*)(vv::Vector{DataType},plm::NLregModF)
    dt = Dtrans(vv)
    length(dt) == size(plm,1) || error("DimensionMisMatch")
    CompositePLModF(dt,plm)
end

function initpars(cp::CompositePLModF)
    ini = initpars(plmod(cp))
    pt = ptran(cp)
    res = similar(ini) # a comprehension might be better but the result has type Vector{Any}
    for i in 1:length(ini)
        res[i] = invtrans(pt[i],ini[i])
    end
    res
end

mmjac(cp::CompositePLModF) = mmjac(plmod(cp))

StatsBase.model_response(cp::CompositePLModF) = model_response(plmod(cp))

npars(cp::CompositePLModF) = npars(cp.plm)

plmod(cp::CompositePLModF) = cp.plm

function pnames(cp::CompositePLModF)
    plm = plmod(cp)
    nnl,nl,_ = size(plm)
    innms = pnames(plm)
    [innms[1:nl],pnames(ptran(cp),innms[nl + (1:nnl)])]
end

ptran(cp::CompositePLModF) = cp.pt

StatsBase.residuals(cp::CompositePLModF) = residuals(cp.plm)

Base.size(cp::CompositePLModF) = size(plmod(cp))

Base.size(cp::CompositePLModF,args...) = size(plmod(cp),args...)

tgrad(cp::CompositePLModF) = tgrad(cp.plm)

function updtMM!(cp::CompositePLModF,nlpars::Vector)
    plm = plmod(cp)
    nnl,nl,n = size(plm)
    length(nlpars) == nnl || throw(DimensionMismatch(""))
    tpars, jac = parjac(ptran(cp), nlpars)
    res = updtMM!(plm,tpars)
    if isa(jac,Diagonal)                # can update mmjac in place
        A_mul_B!(jac,reshape(mmjac(plm),(nnl,nl*n)))
    else
        mm = mmjac(plm)
        updtmmd = jac * reshape(mm,(nnl,nl*n))
        mm[:] = updtmmd
    end
    res
end

## There should be a version of this that does not need to update the nlpars part.
## See the definition of the deviance method for PLinearLS, Vector in plreg.jl
function updtmu!(cp::CompositePLModF,pars::Vector)
    plm = plmod(cp)
    nnl,nl,n = size(plm)
    length(pars) == nnl + nl || error("DimensionMisMatch")
    tpars, jac = parjac(ptran(cp),pars[nl + (1:nnl)])
    res = updtmu!(plm,[pars[1:nl],tpars])
    if isa(jac,Diagonal)
        A_mul_B!(jac,reshape(mmjac(plm),(nnl,nl*n)))
        scale!([ones(nl),diag(jac)],tgrad(plm))
    else
        error("code not yet written")
    end
    res
end

ext: jl


text:
"""
    chromatic_adaptation(sourcecolor::Color, sourceWP::Color, destWP::Color = D_series_whitepoint(6504), response::Symbol = :cat02)

Adapt source color according to given source white point and destination white point.
Default destination white point is `D65`. Default cone response is `:cat02`.

# Examples
```jldoctest
julia> chromatic_adaptation(XYZ(0.1, 0.7, 0.05), D_series_whitepoint(4000))
XYZ{Float64}(0.04957926027638688,0.7081276413719826,0.09614857166385628)

julia> chromatic_adaptation(XYZ(0.1, 0.7, 0.05), D_series_whitepoint(4000), D_series_whitepoint(7500), :bradford)
XYZ{Float64}(0.0513224358507469,0.7251796437072336,0.06113763051342081)
```
"""
function chromatic_adaptation(sourcecolor::Color, sourceWP::Color, destWP::Color = D_series_whitepoint(6504), response::Symbol = :cat02)
    MA=coneresponse(response)
    ργβ_S=MA*[sourceWP.x,sourceWP.y,sourceWP.z]
    ργβ_D=MA*[destWP.x,destWP.y,destWP.z]
    cd=MA^-1 *[ργβ_D[1]/ργβ_S[1] 0.0               0.0;
                    0.0               ργβ_D[2]/ργβ_S[2] 0.0;
                    0.0               0.0               ργβ_D[3]/ργβ_S[3]] * MA * [sourcecolor.x, sourcecolor.y, sourcecolor.z]
    # destination color
    Colors.XYZ(cd[1],cd[2],cd[3])
end

"""
    coneresponse(response::Symbol)

Return cone response matrix according to chromatic adaptation transform `response`.

Available transforms:

* `:xyzscaling` Simplest (and worst) method.
* `:vonkries`
* `:bradford` Linearized Bradfod transform. Diagonal Matrix Transform (DMT)-type vonKries-Ives transform.
* `:sharp` Sharp Adaptation Transform. First described by Finlayson et al.
* `:fairchild`
* `:cat97s` CIECAM97s adaptation transform
* `:cat02` Default transform. Developed for the CIECAM02 color appearance model.
* `:hpe` Hunt-Pointer-Estevez transform
* `:bestof20` by *Sabine Süsstrunk*: [Computing Chromatic Adaptation](https://core.ac.uk/download/pdf/147990632.pdf), PhD thesis, 2005, page 104. Sharpened linear cone response matrix optimized for white-point independence, gamut coverage, and hue constancy. Delivers better results than `:cat02` in most cases. Added for experimental purposes.

# Examples

```jldoctest
julia> coneresponse(:xyzscaling)
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> coneresponse(:bradford)
3×3 Matrix{Float64}:
  0.8951   0.2664  -0.1614
 -0.7502   1.7135   0.0367
  0.0389  -0.0685   1.0296

julia> coneresponse(:bestof20)
3×3 Matrix{Float64}:
  1.6351  -0.4071  -0.228
 -0.8044   1.7798   0.0246
  0.0      0.0152   0.9848
```
"""
function coneresponse(response::Symbol = :cat02)
    if response == :xyzscaling
        MA = [1.0 0.0 0.0;
              0.0 1.0 0.0;
              0.0 0.0 1.0]
    elseif response == :vonkries
        MA = [ 0.4002400  0.7076000 -0.0808100;
              -0.2263000  1.1653200  0.0457000;
               0.0000000  0.0000000  0.9182200]
    elseif response == :bradford
        MA = [ 0.8951000  0.2664000 -0.1614000;
              -0.7502000  1.7135000  0.0367000;
               0.0389000 -0.0685000  1.0296000]
    elseif response == :sharp
        MA = [1.2694 -0.0988 -0.1706;
             -0.8364  1.8006  0.03570;
              0.0297 -0.0315  1.0018]
    elseif response == :fairchild
        MA = [0.8562 0.3372 -0.1934;
             -0.8360 1.8327  0.00330;
              0.0357 -0.0469 1.0112]
    elseif response == :cat97s
        MA = [0.8562  0.3372 -0.1934;
             -0.8360  1.8327  0.0033;
              0.0357 -0.0469  1.0112]
    elseif response == :cat02
        MA = [0.7328  0.4296 -0.1624;
             -0.7036  1.6974  0.0061;
              0.0030  0.0136  0.9834]
    elseif response == :hpe
        MA = [0.3897  0.6890 -0.0787;
             -0.2298  1.1834  0.0464;
              0.0     0.0     1.0]
    elseif response == :bestof20
        MA = [1.6351 -0.4071 -0.2280;
             -0.8044  1.7798  0.0246;
              0.0000  0.0152  0.9848]
    else
        throw(DomainError(response, "Cone response matrix does not exist."))
    end
    MA
end

"""
    const INGLING_TSOU

Ingling and Tsou’s suprathreshold transformation matrix.
Converts l,m,s cone fundamentals to color opponent space (Vλ, r-g, y-b) / (WS(λ), YB(λ), RG(λ))
"""
const INGLING_TSOU=
    [0.600  0.400  0.000;
     0.240  0.105 -0.700;
     1.200 -1.600 0.400]

"""
    inglingtsou_opponent(conefund::ConeFund)

Converts cone fundamentals to color opponent channels using the Ingling-Tsou transformation matrix.

Signal flow of the two-stage model of human color vision:

        L              M                    S   photoreceptor
        ┃              ┃                    ┃      stage
        ┃              ┃     ┏━━━━━━━━━━━━━━┫      (CMFs)
        ┃  ┏━━━━━━━━━━━╋━━━━━┃━━━━━━━━━┓    ┃
        ┃  ┃           ┃     ┃         ┃    ┃
     ┏━━┻━━┃━━━━━┳━━━━━┃━━━━━┃━━━━┓    ┃    ┃
     ┃     ┃     ┃     ┃     ┃    ┃    ┃    ┃
     ┃     ┃     ┃     ┃     ┃    ┃    ┃    ┃
     ┃     ┃     ┃     ┃     ┃    ┃    ┃    ┃
    0.6   0.4  0.24  0.105  0.7  1.2  1.6  0.4   factors
     ║     ║     ║     ║     ║    ║    ║    ║
     ╚══╦══╝     ╚═════╬═════╝    ╚════╬════╝
        ║              ║               ║
        ║              ║               ║
     L  +  M     L  +  M  -  S    L  - M +  S     opponent
        ║              ║               ║           stage
        Vλ            y-b             r-g

    lightness     yellow-blue       red-green

# Example

```jldoctest
julia> inglingtsou_opponent(conefund(SPECENV.cmf))
VYBRG(Real[390.0, 391.0, 392.0, 393.0, 394.0, 395.0, 396.0, 397.0, 398.0, 399.0  …  821.0, 822.0, 823.0, 824.0, 825.0, 826.0, 827.0, 828.0, 829.0, 830.0], Real[0.0003878623982, 0.0004737047352, 0.000577319691, 0.0007018085474, 0.0008506078488000001, 0.0010274557607999999, 0.001236329446, 0.0014813467259999998, 0.001766625772, 0.0020960976039999996  …  9.411752703999999e-7, 8.876810727999999e-7, 8.372313131999999e-7, 7.896833636e-7, 7.448974792000001e-7, 7.027372568e-7, 6.630699992e-7, 6.257670254799999e-7, 5.9070389244e-7, 5.5776058904e-7], Real[-0.004164411112444999, -0.00504460427438, -0.006110193693075, -0.007397180969524999, -0.00894723637031, -0.010808111330739998, -0.013033902761799998, -0.01568509182405, -0.018828265622099997, -0.02253541766325  …  3.6881237642999995e-7, 3.4782318591e-7, 3.2803082274e-7, 3.0937903046999996e-7, 2.9181270796499996e-7, 2.7527808831e-7, 2.5972288066499997e-7, 2.450963909835e-7, 2.31349607448e-7, 2.18435275848e-7], Real[0.0023730382284000004, 0.0028717459496, 0.003474333342, 0.00420068739, 0.0050737833632, 0.0061199115808, 0.007368828876, 0.008853794675999999, 0.010611444512, 0.012681445520000003  …  1.5481949744e-6, 1.4590315327999997e-6, 1.3750428792e-6, 1.2959784376e-6, 1.2215934872e-6, 1.1516498848e-6, 1.0859166232e-6, 1.02417030368e-6, 9.661954418399999e-7, 9.117847518400001e-7])
```
"""
function inglingtsou_opponent(conefund::ConeFund)
    it = INGLING_TSOU * [conefund.l, conefund.m, conefund.s]
    VYBRG(conefund.λ, it[1], it[2], it[3])
end

#=
0 	1 	2 	3 	4 	5 	6 	7 	8 	9 	A 	B 	C 	D 	E 	F
U+250x 	─ 	━ 	│ 	┃ 	┄ 	┅ 	┆ 	┇ 	┈ 	┉ 	┊ 	┋ 	┌ 	┍ 	┎ 	┏
U+251x 	┐ 	┑ 	┒ 	┓ 	└ 	┕ 	┖ 	┗ 	┘ 	┙ 	┚ 	┛ 	├ 	┝ 	┞ 	┟
U+252x 	┠ 	┡ 	┢ 	┣ 	┤ 	┥ 	┦ 	┧ 	┨ 	┩ 	┪ 	┫ 	┬ 	┭ 	┮ 	┯
U+253x 	┰ 	┱ 	┲ 	┳ 	┴ 	┵ 	┶ 	┷ 	┸ 	┹ 	┺ 	┻ 	┼ 	┽ 	┾ 	┿
U+254x 	╀ 	╁ 	╂ 	╃ 	╄ 	╅ 	╆ 	╇ 	╈ 	╉ 	╊ 	╋ 	╌ 	╍ 	╎ 	╏
U+255x 	═ 	║ 	╒ 	╓ 	╔ 	╕ 	╖ 	╗ 	╘ 	╙ 	╚ 	╛ 	╜ 	╝ 	╞ 	╟
U+256x 	╠ 	╡ 	╢ 	╣ 	╤ 	╥ 	╦ 	╧ 	╨ 	╩ 	╪ 	╫ 	╬ 	╭ 	╮ 	╯
U+257x 	╰ 	╱ 	╲ 	╳ 	╴ 	╵ 	╶ 	╷ 	╸ 	╹ 	╺ 	╻ 	╼ 	╽ 	╾ 	╿
=#

"""
    const OCS

Combination of M_XYZ→YCC from Kotera with M_LMS→XYZ of CIECAM02
"""
const OCS =
    [0.489  0.134  0.378;
     1.004 -0.790 -0.214;
     0.034 -0.472  0.437]

function mlms_ycc()
    OCS
end

const KOTERA =
    [0.593247  -0.0750709   0.482816;
     0.743587  -0.652887   -0.0907655;
    -0.181412  -0.23549     0.415879]

const IPT =
    [0.4000  0.4000  0.2000;
     4.4550 -4.8510  0.3960;
     0.8056  0.3572 -1.1628]

const IPTinv =
    [1.8501 -1.1383  0.2385;
     0.3668  0.6439 -0.0107;
     0.0000  0.0000  1.0889]

function IPTforward(cfund::ConeFund)
    cfund.l ≥ 0 ? (Lp = cfund.l^0.43) :
    cfund.l < 0 ? (Lp = -(-cfund.l)^0.43) :
    cfund.m ≥ 0 ? (Mp = cfund.m^0.43) :
    cfund.m < 0 ? (Mp = -(-cfund.m)^0.43) :
    cfund.s ≥ 0 ? (Sp = cfund.s^0.43) :
    cfund.s < 0 ? (Sp = -(-cfund.s)^0.43) : nothing

    IPT * [Lp Mp Sp]
end

function IPTinverse(iptcol)
    lmsp = IPTinv * iptcol
    lmsp.l ≥ 0 ? l = lmsp.l^2.3256 :
    lmsp.l < 0 ? l = -(-lmsp.l)^2.3256 :
    lmsp.M ≥ 0 ? m = lmsp.m^2.3256 :
    lmsp.m < 0 ? m = -(-lmsp.m)^2.3256 :
    lmsp.s ≥ 0 ? s = lmsp.s^2.3256 :
    lmsp.s < 0 ? s = -(-lmsp.m)^2.3256 : nothing

    xyz = inv(IPT) * [l m s]
    Colors.XYZ(xyz[1], xyz[2], xyz[3])
end


text:
mutable struct BD{T}
	om::ObsModel{T}
	optm::OptimModel{T}
	gx::X{T}
	sx::X{T}
	sxp::Sxparam
	snorm_flag::Bool 	# restrict s along a unit circle during optimization
	snormmat::Matrix{T}            # stored outer product of s
	dsnorm::Vector{T}		# gradient w.r.t. normalized selet
	attrib_inv::Symbol
	verbose::Bool
	err::DataFrames.DataFrame
	opG::LinearMaps.LinearMap{T}
	opS::LinearMaps.LinearMap{T}
	pbandpass::P_bandpass{T}
end

"""
Constructor for the blind deconvolution problem
"""
function BD(ntg, nt, nr, nts;
	slags=[nts-1, 0],
	dlags=[nt-1, 0],
	glags=[ntg-1, 0],
	    gprecon_attrib=:none,
	       sxp=Sxparam(1,:all),
	       gweights=nothing,
	       goptim=nothing,
	       gαvec=nothing,
	       soptim=nothing,
	       sαvec=nothing,
	       sprecon=nothing,
	       snorm_flag=false,
	       fft_threads=false,
	       fftwflag=FFTW.PATIENT,
	       dobs=nothing, gobs=nothing, sobs=nothing, verbose=false, attrib_inv=:g,
	       size_check=false,
	       )
	# determine type of IBD
	if(!(dobs===nothing))
		(size(dobs)≠(nt,nr)) && error("dobs size error")
		T=eltype(dobs)
	else
		T1=eltype(sobs)
		T2=eltype(gobs)
		(T1≠T2) ? error("type difference") : (T=T1)
		(size(sobs)≠(nts,)) && error("sobs size error")
		(size(gobs)≠(ntg,nr)) && error("sobs size error")
	end

	if(size_check)
		if(ntg+nts-1 ≠ nt)
			error("invalid sizes for convolutional model")
		end
	end

	# use maximum threads for fft
	fft_threads &&  (FFTW.set_num_threads(Sys.CPU_THREADS))

	# store observed data
	om=ObsModel(ntg, nt, nr, nts, T, d=dobs, g=gobs, s=sobs,
	slags=slags,
	dlags=dlags,
	glags=glags, )

	# create models depending on mode
	optm=OptimModel(ntg, nt, nr, nts, T, fftwflag=fftwflag,
	slags=slags,
	dlags=dlags,
	glags=glags,
		 )

	# inversion variables allocation
	gx=X(length(optm.cal.g), T)
	sx=X(length(optm.cal.s), T)

	snorm_flag ?	(snormmat=zeros(T,nts, nts)) : (snormmat=zeros(T,1,1))
	snorm_flag ?	(dsnorm=zeros(T,nts)) : (dsnorm=zeros(T,1))

	err=DataFrame(g=[], g_nodecon=[], s=[], d=[], front_load=[], whiteness=[])

	# dummy
	opG=LinearMap{T}(x->0.0, y->0.0,1,1, ismutating=true)
	opS=LinearMap{T}(x->0.0, y->0.0,1,1, ismutating=true)

	pbandpass=P_bandpass(T, fmin=0.2, fmax=0.4, nt=optm.nt)

	pa=BD(
		om,		optm,			gx,		sx,	sxp, snorm_flag,
		snormmat,		dsnorm,		attrib_inv,		verbose,
		err, opG, opS,pbandpass)		# trying to penalize the energy in the correlations of g (not in practice),

	# update operators
	pa.opG=create_operator(pa, G())
	pa.opS=create_operator(pa, S())

	gobs=pa.om.g
	sobs=pa.om.s
	dobs=pa.om.d

	# obs.g <-- gobs
	replace!(pa.optm, gobs, :obs, :g )
	# obs.s <-- sobs
	replace!(pa.optm, sobs, :obs, :s )
	# obs.d <-- dobs
	copyto!(pa.optm.obs.d, dobs) #

	add_precons!(pa, pa.om.g, attrib=gprecon_attrib)

	initialize!(pa)
	#update_func_grad!(pa,goptim=goptim,soptim=soptim,gαvec=gαvec,sαvec=sαvec)

	return pa

end

function update_prepare!(pa::BD, ::S)
end
function update_prepare!(pa::BD, ::G)
end
function update_finalize!(pa::BD, ::S)

	if(STF_FLAG)
		for i in eachindex(pa.optm.cal.s)
			if(pa.optm.cal.s[i] <0.0)
				pa.optm.cal.s[i]=0.0
			end
		end
	end
end
function update_finalize!(pa::BD, ::G)
end

function model_to_x!(x, pa::BD, ::S)
	for i in eachindex(x)
		x[i]=^(pa.optm.cal.s[i],inv(pa.sxp.n))*pa.sx.precon[i]
	end
	return nothing
end
function model_to_x!(x, pa::BD, ::G)
	for i in eachindex(x)
		x[i]=pa.optm.cal.g[i]*pa.gx.precon[i] 		# multiply by gprecon
	end
	return nothing
end

function x_to_model!(x, pa::BD, ::S)
	for i in eachindex(pa.optm.cal.s)
		# put same in all receivers
		pa.optm.cal.s[i]=^(x[i],pa.sxp.n)*pa.sx.preconI[i]
	end
	if(pa.snorm_flag)
		xn=vecnorm(x)
		scale!(pa.optm.cal.s, inv(xn))
	end
	return nothing
end
function x_to_model!(x, pa::BD, ::G)
	for i in eachindex(pa.optm.cal.g)
		pa.optm.cal.g[i]=x[i]*pa.gx.preconI[i]
	end
	return pa
end

"""
Create preconditioners using the observed Green Functions.
* `cflag` : impose causaulity by creating gprecon using gobs
* `max_tfrac_gprecon` : maximum length of precon windows on g
"""
function add_precons!(pa::BD, gobs; αexp=0.0, cflag=true,
		       max_tfrac_gprecon=1.0, attrib=:focus)

	ntg=pa.om.ntg
	nts=pa.om.nts

	ntgprecon=round(Int,max_tfrac_gprecon*ntg);

	nr=size(gobs,2)
	sprecon=ones(nts)
	gprecon=ones(ntg, nr);
	gweights=ones(ntg, nr);
	if(attrib==:windows)
		for ir in 1:nr
			g=normalize(view(gobs,:,ir))
			indz=findfirst(x->abs(x)>1e-6, g)
		#	if(indz > 1)
		#		indz -= 1 # window one sample less than actual
		#	end
			if(!cflag && indz≠0)
				indz=1
			end
			if(indz≠0)
				for i in 1:indz-1
					gprecon[i,ir]=0.0
					gweights[i,ir]=0.0
				end
				for i in indz:indz+ntgprecon
					if(i≤ntg)
						gweights[i,ir]=exp(αexp*(i-indz-1)/ntg)  # exponential weights
						gprecon[i,ir]=exp(αexp*(i-indz-1)/ntg)  # exponential weights
					end
				end
				for i in indz+ntgprecon+1:ntg
					gprecon[i,ir]=0.0
					gweights[i,ir]=0.0
				end
			else
				gprecon[:,ir]=0.0
				gweights[:,ir]=0.0
			end
		end
	elseif(attrib==:focus)
		ir0=inear(gobs)
		# first receiver is a spike
		gprecon[2:end,ir0]=0.0

	end

	add_gprecon!(pa, gprecon)
	add_gweights!(pa, gweights)
	add_sprecon!(pa, sprecon)

	return pa
end

"return index of closest receiver"
function inear(gobs, threshold=1e-6)
	nr=size(gobs,2)
	ir0=argmin([findfirst(x->abs(x)>threshold, vec(gobs[:,ir])) for ir in 1:nr])
	return ir0
end

function bd!(pa::BD, io=stdout; tol=1e-6)

	if(io===nothing)
		logfilename=joinpath(pwd(),string("XBD",Dates.now(),".log"))
		io=open(logfilename, "a+")
	end

	update_all!(pa, io, max_reroundtrips=1, max_roundtrips=100000, roundtrip_tol=tol)

	err!(pa)
end

function F!(pa::BD,x::AbstractVector, attrib::S)
	compute=(x!=pa.sx.last_x)
	if(compute)
		x_to_model!(x, pa, attrib) #
		copyto!(pa.sx.last_x, x)
		Conv.mod!(pa.optm.cal, Conv.D()) # modify pa.optm.cal.d
		return pa
	end
end

function F!(pa::BD,x::AbstractVector, ::G)
	compute=(x!=pa.gx.last_x)
	if(compute)
		x_to_model!(x, pa, G())
		copyto!(pa.gx.last_x, x)
		Conv.mod!(pa.optm.cal, Conv.D()) # modify pa.optm.cal.d
		return pa
	end
end

"""
Apply Fadj to dcal
"""
function Fadj!(pa::BD, storage, x, dcal, ::S)
	fill!(storage, 0.0)
	Conv.mod!(pa.optm.cal, Conv.S(), d=dcal, s=pa.optm.ds)
	for j in 1:size(pa.optm.ds,1)
		storage[j] = pa.optm.ds[j]
	end

	# apply precon
	for i in eachindex(storage)
		if(iszero(pa.sx.precon[i]))
			storage[i]=0.0
		else
			storage[i] = storage[i]*pa.sxp.n*^(x[i],pa.sxp.n-1)*pa.sx.preconI[i]
		end
	end
	# factor, because s was divided by norm of x
	if(pa.snorm_flag)
		copyto!(pa.optm.dsnorm, storage)
		Misfits.derivative_vector_magnitude!(storage,pa.optm.dsnorm,x,pa.snormmat)
	end
	return storage
end

function Fadj!(pa::BD, storage, x, dcal, ::G)
	Conv.mod!(pa.optm.cal, Conv.G(), g=pa.optm.dg, d=dcal)

	for i in eachindex(storage)
		if(iszero(pa.gx.precon[i]))
			storage[i]=0.0
		else
			storage[i]=pa.optm.dg[i]/pa.gx.precon[i]
		end
	end
	return storage
end

function initialize!(pa::BD)
	# starting random models
	for i in eachindex(pa.optm.cal.s)
		x=(pa.sx.precon[i]≠0.0) ? randn() : 0.0
		if(STF_FLAG)
			pa.optm.cal.s[i]=abs(x)
		else
			pa.optm.cal.s[i]=x
		end
	end
	for i in eachindex(pa.optm.cal.g)
		x=(pa.gx.precon[i]≠0.0) ? randn() : 0.0
		pa.optm.cal.g[i]=x
	end
end

"""
Compute front load, independent of translation, assuming that the most front-loaded channel is 1.
"""
function front_load(x)
	nt=size(x,1)
	nr=size(x,2)
	err=[]
	for is in 1:size(x,1)
		x1=circshift(x,(is,0))

		J=0.0
		for it in 1:nt
			J += abs2(x1[it,1]) * abs2(it-1)
		end
		push!(err,J)
	end
	return minimum(err)
end

"""
compute errors
update pa.err
print?
give either cal
"""
function err!(pa::BD, io=stdout; cal=pa.optm.cal)
	#xg_nodecon=hcat(Conv.xcorr(pa.om.d,Conv.P_xcorr(pa.om.nt, pa.om.nr, cglags=[pa.optm.ntg-1, pa.optm.ntg-1]))...)
	#xgobs=hcat(Conv.xcorr(pa.om.g)...) # compute xcorr with reference g
	#fs = Misfits.error_after_normalized_autocor(cal.s, pa.optm.obs.s)
	#xgcal=hcat(Conv.xcorr(cal.g)...) # compute xcorr with reference g
	#fg = Misfits.error_squared_euclidean!(nothing, xgcal, xgobs, nothing, norm_flag=true)
	#fg_nodecon = Misfits.error_squared_euclidean!(nothing, xg_nodecon, xgobs, nothing, norm_flag=true)
	fg_nodecon= 0.0

	f = Misfits.error_squared_euclidean!(nothing, cal.d, pa.optm.obs.d, nothing, norm_flag=true)
	#fs = Misfits.error_after_translation(cal.s, pa.optm.obs.s)[1]
	fs = 0.0
	fg = Misfits.error_after_translation(cal.g, pa.optm.obs.g)[1]

	whiteness=Conv.func_grad!(nothing, cal.g, Conv.P_misfit_weighted_acorr(pa.om.ntg,pa.om.nr))

	fload=front_load(cal.g) # assuming ir ==1!?

	push!(pa.err[:s],fs)
	push!(pa.err[:d],f)
	push!(pa.err[:g],fg)
	push!(pa.err[:whiteness],whiteness)
	push!(pa.err[:front_load],fload)
	push!(pa.err[:g_nodecon],fg_nodecon)
	if(!(io===nothing))
		write(io,"Blind Decon Errors\t\n")
		write(io,"==================\n")
		write(io, string(pa.err))
	end
end

function update_g!(pa::BD, xg)
	pa.attrib_inv=:g
	fg = update!(pa, xg)
	return fg
end

function update_s!(pa::BD, xs)
	pa.attrib_inv=:s
	fs = update!(pa, xs)
	return fs
end

#=
struct BandLimit <: Manifold
end

function retract!(::BandLimit,x)
	#=
	x_to_model!(x, pa) # modify pa.optm.cal.s or pa.optm.cal.g

	pac=pa.optm.cal
	if(pa.attrib_inv==:g)
	Conv.pad_truncate!(pac.g, pac.gpad, pac.glags[1], pac.glags[2], pac.np2, 1)
	A_mul_B!(pac.gfreq, pac.gfftp, pac.gpad)

	nr=pa.om.nr
	for ir in 1:nr
		for i in 1:10
			pac.gfreq[i,ir]=complex(0.0,0.0)
		end
	end
	A_mul_B!(pac.gpad, pac.gifftp, pac.gfreq)
	Conv.pad_truncate!(pac.g, pac.gpad, pac.glags[1], pac.glags[2], pac.np2, -1)
	=#
end
export retract!

function project_tangent!(::BandLimit,g,x)

end
=#


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

"""
    ValueShapes.default_datatype(T::Type)

Return a default specific type U that is more specific than T, with U <: T.

e.g.

    ValueShapes.default_datatype(Real) == Float64
    ValueShapes.default_datatype(Complex) == Complex{Float64}
"""
function default_datatype end

@inline default_datatype(::Type{>:Int}) = Int
@inline default_datatype(::Type{>:Float64}) = Float64
@inline default_datatype(::Type{>:Real}) = Float64
@inline default_datatype(::Type{>:Complex{Float64}}) = Complex{Float64}
@inline default_datatype(T::Type) = T

"""
    abstract type AbstractValueShape

An `AbstractValueShape` combines type and size information.

Subtypes are defined for shapes of scalars (see [`ScalarShape`](@ref)),
arrays (see [`ArrayShape`](@ref)), constant values
(see [`ConstValueShape`](@ref)) and `NamedTuple`s (see
[`NamedTupleShape`](@ref)).

Subtypes of `AbstractValueShape` must support `eltype`, `size` and
[`totalndof`](@ref).

Value shapes can be used as constructors to generate values of the given
shape with undefined content. If the element type of the shape is an abstract
or union type, a suitable concrete type will be chosen automatically, if
possible (see [`ValueShapes.default_datatype`](@ref)):

```julia
shape = ArrayShape{Real}(2,3)
A = shape(undef)
typeof(A) == Array{Float64,2}
size(A) == (2, 3)
valshape(A) == ArrayShape{Float64}(2,3)
```

Use

    (shape::AbstractValueShape)(data::AbstractVector{<:Real})::eltype(shape)

to view a flat vector of anonymous real values
as a value of the given shape:

```julia
data = [1, 2, 3, 4, 5, 6]
shape(data) == [1 3 5; 2 4 6]
```

In return,

    Base.Vector{T}(undef, shape::AbstractValueShape)
    Base.Vector(undef, shape::AbstractValueShape)

will create a suitable uninitialized vector of the right length to hold such
flat data for the given shape. If no type `T` is given, a suitable data
type will be chosen automatically.

When dealing with multiple vectors of flattened data, use

    shape.(data::ArrayOfArrays.AbstractVectorOfSimilarVectors)

ValueShapes supports this via specialized broadcasting.

In return,

    ArraysOfArrays.VectorOfSimilarVectors{T}(shape::AbstractValueShape)
    ArraysOfArrays.VectorOfSimilarVectors(shape::AbstractValueShape)

will create a suitable vector (of length zero) of vectors that can hold
flattened data for the given shape. The result will be a
`VectorOfSimilarVectors` wrapped around a 2-dimensional `ElasticArray`.
This way, all data is stored in a single contiguous chunk of memory.
"""
abstract type AbstractValueShape end

export AbstractValueShape

# Value shapes behave as scalars under broadcasting:
@inline Base.Broadcast.broadcastable(shape::AbstractValueShape) = Ref(shape)

function _valshapeoftype end

"""
    valshape(x)::AbstractValueShape
    valshape(acc::ValueAccessor)::AbstractValueShape
    valshape(d::Distributions.Distribution)::AbstractValueShape

Get the value shape of an arbitrary value, resp. the shape a `ValueAccessor`
is based on, or the shape of the variates for a `Distribution`.
"""
function valshape end
export valshape

@inline valshape(x::T) where T = _valshapeoftype(T)

"""
    totalndof(shape::AbstractValueShape)

Get the total number of degrees of freedom of values of the given shape.

Equivalent to the length of a vector that would result from flattening the
data into a sequence of real numbers, excluding any constant values.
"""
function totalndof end
export totalndof

@inline nonabstract_eltype(shape::AbstractValueShape) = default_datatype(eltype(shape))

function _checkcompat(shape::AbstractValueShape, data::AbstractVector{<:Real})
    n_shape = totalndof(shape)
    n_data = length(eachindex(data))
     if n_shape != length(eachindex(data))
        throw(ArgumentError("Data vector of length $(n_data) incompatible with value shape with $(n_shape) degrees of freedom"))
    end
    nothing
end

@static if VERSION >= v"1.3"
    @inline function (shape::AbstractValueShape)(data::AbstractVector{<:Real})
        @boundscheck _checkcompat(shape, data)
        data[ValueAccessor(shape, 0)]
    end
else
    # Workaround for Julia issue #14919
    @inline function _apply_shape_to_data(shape::AbstractValueShape, data::AbstractVector{<:Real})
        @boundscheck _checkcompat(shape, data)
        data[ValueAccessor(shape, 0)]
    end
end

Base.Vector{T}(::UndefInitializer, shape::AbstractValueShape) where T =
    Vector{T}(undef, totalndof(shape))

Base.Vector(::UndefInitializer, shape::AbstractValueShape) =
    Vector{nonabstract_eltype(shape)}(undef, totalndof(shape))

ArraysOfArrays.VectorOfSimilarVectors{T}(shape::AbstractValueShape) where T =
    VectorOfSimilarVectors(ElasticArray{T}(undef, totalndof(shape), 0))

ArraysOfArrays.VectorOfSimilarVectors(shape::AbstractValueShape) =
    VectorOfSimilarVectors(ElasticArray{nonabstract_eltype(shape)}(undef, totalndof(shape), 0))

const VSBroadcasted1{F,T} = Base.Broadcast.Broadcasted{
    <:Base.Broadcast.AbstractArrayStyle{1},
    <:Any,
    F,
    <:Tuple{T}
}

const VSBroadcasted2{F,T1,T2} = Base.Broadcast.Broadcasted{
    <:Base.Broadcast.AbstractArrayStyle{1},
    <:Any,
    F,
    <:Tuple{T1,T2}
}

# Specialize (::AbstractValueShape).(::AbstractVector{<:AbstractVector}):
Base.copy(instance::VSBroadcasted1{<:AbstractValueShape,AbstractVector{<:AbstractVector}}) =
    broadcast(getindex, instance.args[1], ValueAccessor(instance.f, 0))


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Pango_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Pango")
JLLWrappers.@generate_main_file("Pango", UUID("36c8627f-9965-5494-a995-c6b170f724f3"))
end  # module Pango_jll


text:
function Base.map(f, t::DDataset; select=nothing)
    # TODO: fix when select has a user-supplied vector
    delayedmap(t.chunks) do x
        map(f, x; select=select)
    end |> fromchunks
end

function dropmissing(t::DDataset, select=(colnames(t)...,))
    delayedmap(t.chunks) do x
        dropmissing(x, select)
    end |> fromchunks
end

function convertmissing(t::DDataset, missingtype)
    delayedmap(t.chunks) do x
        convertmissing(x, missingtype)
    end |> fromchunks
end

function Base.filter(f, t::DDataset; select=isa(f, Union{Tuple, Pair}) ? nothing : valuenames(t))
    g = delayed(x -> filter(f, x; select=select))
    cache_thunks(fromchunks(map(g, t.chunks)))
end

function selectkeys(x::DNDSparse, which; kwargs...)
    ndsparse(rows(keys(x), which), values(x); kwargs...)
end

function selectvalues(x::DNDSparse, which; kwargs...)
    ndsparse(keys(x), rows(values(x), which); kwargs...)
end

Base.@deprecate select(x::DNDSparse, conditions::Pair...) filter(conditions, x)
Base.@deprecate select(x::DNDSparse, which::DimName...; kwargs...) selectkeys(x, which; kwargs...)

"""
    convertdim(x::DNDSparse, d::DimName, xlate; agg::Function, name)

Apply function or dictionary `xlate` to each index in the specified dimension.
If the mapping is many-to-one, `agg` is used to aggregate the results.
`name` optionally specifies a name for the new dimension. `xlate` must be a
monotonically increasing function.

See also [`reduce`](@ref)
"""
function convertdim(t::DNDSparse{K,V}, d::DimName, xlat;
                    agg=nothing, vecagg=nothing, name=nothing) where {K,V}

    if isa(d, Symbol)
        dn = findfirst(dimlabels(t), d)
        if dn == 0
            throw(ArgumentError("table has no dimension \"$d\""))
        end
        d = dn
    end

    chunkf(c) = convertdim(c, d, xlat; agg=agg, vecagg=nothing, name=name)
    chunks = map(delayed(chunkf), t.chunks)

    xlatdim(intv, d) = Interval(tuplesetindex(first(intv), xlat(first(intv)[d]), d),
                                tuplesetindex(last(intv),  xlat(last(intv)[d]), d))

    # TODO: handle name kwarg
    # apply xlat to bounding rectangles
    domains = map(t.domains) do space
        nrows = agg === nothing ? space.nrows : Nullable{Int}()
        IndexSpace(xlatdim(space.interval, d), xlatdim(space.boundingrect, d), nrows)
    end

    t1 = DNDSparse{eltype(domains[1]),V}(domains, chunks)

    if agg !== nothing && has_overlaps(domains)
        overlap_merge(x, y) = merge(x, y, agg=agg)
        chunk_merge(ts...)  = _merge(overlap_merge, ts...)
        cache_thunks(rechunk(t1, merge=chunk_merge, closed=true))
    elseif vecagg != nothing
        groupby(vecagg, t1) # already cached
    else
        cache_thunks(t1)
    end
end

keyindex(t::DNDSparse, i::Int) = i
keyindex(t::DNDSparse{K}, i::Symbol) where {K} = findfirst(x->x===i, fieldnames(K))


text:
function generatePermutations(a, n)
  if n == 1
    println(join(a, ""))
  else
    for i = 1:n-1
      generatePermutations(a, n-1)
      local j = iseven(n) ? i : 1
      a[j], a[n] = a[n], a[j]
    end
    generatePermutations(a, n-1)
  end
end

if length(ARGS) != 1
  println(STDERR, "Exactly one argument is required")
  exit(1)
end
generatePermutations(split(ARGS[1], ""), length(ARGS[1]))

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
# LAR input:  is a collection of five random cubes in the same 3D space
# REMARK:  is not a cellular complex; cells intersect out of boundaries

# julia> @show V,FV,EV;
(V, FV, EV) = ([1.01181 0.215639 0.91979 0.123616 1.02252 0.226347 0.930498 0.134324 0.0458309 -0.301827 0.348275 0.0006172 0.579367 0.23171 0.881811 0.534154 -0.0521776 0.627953 -0.190635 0.489496 -0.0233884 0.656742 -0.161846 0.518285 0.27652 -0.0875132 0.52527 0.161237 0.509324 0.145291 0.758074 0.394041 0.27631 0.564484 0.0058279 0.294002 1.01221 1.30039 0.741732 1.02991; 0.160033 0.0680099 0.956278 0.864255 0.160649 0.0686266 0.956895 0.864872 -0.200245 0.102199 0.417839 0.720283 -0.35354 -0.0510965 0.264543 0.566987 0.682359 0.543901 0.0592036 -0.0792537 0.956374 0.817917 0.333219 0.194762 -0.102028 0.146722 0.324834 0.573584 -0.16916 0.0795901 0.257702 0.506452 0.452356 0.181874 1.15396 0.883481 0.816151 0.545669 1.51776 1.24728; 0.196256 0.206963 0.196872 0.20758 0.997729 1.00844 0.998346 1.00905 0.0677451 0.601282 -0.0855504 0.447986 0.502301 1.03584 0.349005 0.882542 0.159301 0.18809 0.433316 0.462105 0.797002 0.825792 1.07102 1.09981 0.1446 0.377404 0.0774682 0.310272 0.580364 0.813168 0.513232 0.746036 0.403805 1.13971 0.767599 1.5035 0.249344 0.985248 0.613139 1.34904], Array{Int64,1}[[1,2,3,4],[5,6,7,8],[1,2,5,6],[3,4,7,8],[1,3,5,7],[2,4,6,8],[9,10,11,12],[13,14,15,16],[9,10,13,14],[11,12,15,16],[9,11,13,15],[10,12,14,16],[17,18,19,20],[21,22,23,24],[17,18,21,22],[19,20,23,24],[17,19,21,23],[18,20,22,24],[25,26,27,28],[29,30,31,32],[25,26,29,30],[27,28,31,32],[25,27,29,31],[26,28,30,32],[33,34,35,36],[37,38,39,40],[33,34,37,38],[35,36,39,40],[33,35,37,39],[34,36,38,40]],Array{Int64,1}[[1,2],[3,4],[5,6],[7,8],[1,3],[2,4],[5,7],[6,8],[1,5],[2,6],[3,7],[4,8],[9,10],[11,12],[13,14],[15,16],[9,11],[10,12],[13,15],[14,16],[9,13],[10,14],[11,15],[12,16],[17,18],[19,20],[21,22],[23,24],[17,19],[18,20],[21,23],[22,24],[17,21],[18,22],[19,23],[20,24],[25,26],[27,28],[29,30],[31,32],[25,27],[26,28],[29,31],[30,32],[25,29],[26,30],[27,31],[28,32],[33,34],[35,36],[37,38],[39,40],[33,35],[34,36],[37,39],[38,40],[33,37],[34,38],[35,39],[36,40]])


text:
# Main lookup methods.

"""
Find all methods defined by an method definition expression.

    "A docstring for the methods ``f(::Any)`` and ``f(::Any, ::Any)``"
    f(x, y = 1) = x + y

"""
function findmethods(state::State, ex::Expr, codesource)
    line, file  = codesource
    fname, mset = funcname(state, ex), Set{Method}()
    Δ = typemax(Int)
    for m in allmethods(fname)
        # Skip methods defined in other modules.
        Utilities.samemodule(state.mod, m) || continue
        # Get the file and line information for the current method ``m``.
        line′, file′ = lineinfo(m)
        # Filter out methods defined in different files.
        file′ == file || continue
        # Calculate the line distance and update ``mset`` accordingly.
        Δ′ = abs(line′ - line)
        if 0 ≤ Δ′ < Δ
            Δ = Δ′
            empty!(mset); push!(mset, m)
        elseif Δ′ == Δ
            push!(mset, m)
        end
    end
    mset
end

"""
Line number and file name pair for a method ``m``.
"""
function lineinfo(m::Method)
    file, line = Utilities.lsdfield(m, :file), Utilities.lsdfield(m, :line)
    (line, Utilities.expandpath(string(file)))
end

"""
Find the ``Method`` objects referenced by ``(...)`` docstring syntax.

    "Shared docstring for all 2 argument methods, first argument an ``Int``."
    (foo, Int, Any)

"""
function findtuples(state::State, expr::Expr)
    fname = exec(state, exec(state, expr.args[1])) # Run twice to get rid of `QuoteNode`s.
    args = [exec(state, exec(state, arg)) for arg in expr.args[2:end]]
    if isa(fname, Function) && all(x -> isa(x, DataType), args)
        Set{Method}(methods(fname, tup(args...)))
    else
        Set(unshift!(args, fname))
    end
end

if VERSION < v"0.4-dev+4319"
    tup(args...) = tuple(args...)
else
    tup(args...) = Tuple{args...}
end

"""
Find ``Function`` and ``Method`` objects referenced by ``[...]`` syntax.

    "Shared docstring for differently named functions."
    [foobar, foobar!]

"""
function findvcats(state::State, expr::Expr)
    funcs, methods = Set{Function}(), Set{Method}()
    for arg in expr.args
        if isa(arg, Symbol) || isa(arg, QuoteNode)
            push!(funcs, exec(state, exec(state, arg))) # Run twice to remove `Quotenode`s.
        elseif isexpr(arg, :tuple)
            union!(methods, findtuples(state, arg))
        else
            error("Invalid '[objects...]' syntax.")
        end
    end
    funcs, methods
end

# End lookup methods.

# Sandboxed execution of a subset of possible Julian expressions.

"""
Evaluate the expression ``expr`` within the context provided by ``state``.
"""
exec(state::State, expr::Expr) = exec(Head{expr.head}(), state, expr)

exec(H".", state, expr) = getfield(exec(state, expr.args[1]), expr.args[2].value)

exec(H"::", state, expr) = exec(state, expr.args[end])
exec(H"kw", state, expr) = exec(state, expr.args[1])

exec(H":", state, expr) = colon(map(a -> indexer(state, a), expr.args)...)

exec(H"...", state, expr) = Vararg{exec(state, expr.args[1])}

exec(H"curly", state, expr) = exec(state, expr.args[1]){exec(state, expr.args[2:end])...}
exec(H"call", state, expr)  = exec(state, expr.args[1])(exec(state, expr.args[2:end])...)

function exec(H"macrocall", state, expr)
    expr.args[1] == symbol("@doc_str") && return expr.args[end]
    exec(state, exec(state, expr.args[1])(exec(state, expr.args[2:end])...))
end

if VERSION < v"0.4-dev+5551"
    exec(H"triple_quoted_string", state, expr) =
        exec(state, Expr(:macrocall, symbol("@mstr"), expr.args...))
else
    exec(H"triple_quoted_string", state, expr) = string(expr.args...)
end

exec(H"string", state, expr) = string(exec(state, expr.args)...)
exec(H"tuple", state, expr)  = tuple(exec(state, expr.args)...)

exec(H"vcat, vect", state, expr) = vcat(exec(state, expr.args)...)
exec(H"hcat", state, expr)       = hcat(exec(state, expr.args)...)

exec(H"quote", state, expr) = expr.args[1]

exec(H"$", state, expr) = exec(state, expr.args[1])

function exec(H"ref", state, expr)
    object = exec(state, expr.args[1])
    withref(state, object) do
        getindex(object, map(a -> indexer(state, a), expr.args[2:end])...)
    end
end

exec(state::State, args::Vector) = [exec(state, a) for a in args]

function exec(state::State, q::Symbol)
    for scope in reverse(state.scopes)
        haskey(scope, q) && return scope[q]
    end
    getfield(state.mod, q)
end

exec(::State, q::QuoteNode) = q.value
exec(::State, constant)     = constant

function indexer(state::State, arg)
    if arg == symbol("end")
        endof(state.refs[end])
    elseif arg == :(:)
        1:endof(state.refs[end])
    else
        exec(state, arg)
    end
end

# End `exec` methods.

# Helper methods for lookup methods `findmethods` and friends.

"""
Extract the expressions from a ``{}`` in a function definition.
"""
gettvars(expr::Expr) = isexpr(expr.args[1], :curly) ? expr.args[1].args[2:end] : Any[]

"""
Extract the expressions representing a method definition's arguments.
"""
getargs(expr::Expr) = expr.args[2:end]

"""
Return the ``Function`` object represented by a method definition expression.
"""
funcname(state::State, expr::Expr) = funcname(Head{expr.head}(), state, expr)

funcname(H"=, function, call, curly", state::State, expr::Expr) = funcname(state, expr.args[1])

funcname(H"$", state::State, expr::Expr) = funcname(state, getvar(state, expr.args[1]))
funcname(H".", state::State, expr::Expr) = getfield(exec(state, expr.args[1]), expr.args[2].value)

funcname(state::State, q::QuoteNode) = funcname(state, q.value)
funcname(state::State, q::Symbol)    = exec(state, exec(state, q))

funcname(::State, other) = other

mostgeneral(T::DataType) = T{[tvar.ub for tvar in T.parameters]...}
mostgeneral(other)       = other

function allmethods(fname)
    out = Set{Method}()
    for m in methods(fname)
        push!(out, m)
    end

    # In 0.5.0-dev+2007 `mostgeneral` fails when `fname == SharedArray`.
    mg = fname
    try
        mg = mostgeneral(fname)
    catch err
        isa(err, TypeError) || rethrow(err)
    end
    for m in methods(mg)
        push!(out, m)
    end

    out
end


text:
using Test
using JSONPointer
using OrderedCollections

@testset "Basic Tests" begin
    doc = Dict("foo" => ["bar", "baz"],
                ""=> 0,
                "a/b"=> 1,
                "c%d"=> 2,
                "e^f"=> 3,
                "g|h"=> 4,
                "i\\j"=> 5,
                "k\"l"=> 6,
                " "=> 7,
                "m~n"=> 8)

        @test doc[j""] == doc
        @test doc[j"/foo"] == ["bar", "baz"]
        @test doc[JSONPointer.Pointer("/foo/0"; shift_index = true)] == "bar"
        @test doc[j"/foo/1"] == "bar"
        @test doc[j"/"] == 0
        @test doc[j"/a~1b"] == 1
        @test doc[j"/c%d"] == 2
        @test doc[j"/e^f"] == 3
        @test doc[j"/g|h"] == 4
        @test doc[JSONPointer.Pointer("/i\\j")] == 5
        @test doc[j"/k\"l"] == 6
        @test doc[j"/ "] == 7
        @test doc[j"/m~0n"] == 8

        for k in (j"",
                  j"/foo",
                  j"/foo/1",
                  j"/",
                  j"/a~1b",
                  j"/c%d",
                  j"/e^f",
                  j"/g|h",
                  JSONPointer.Pointer("/i\\j"),
                  j"/k\"l",
                  j"/ ",
                  j"/m~0n")
            @test haskey(doc, k)
        end
end

@testset "URI Fragment Tests" begin
    doc = Dict("foo" => ["bar", "baz"],
                ""=> 0,
                "a/b"=> 1,
                "c%d"=> 2,
                "e^f"=> 3,
                "g|h"=> 4,
                "i\\j"=> 5,
                "k\"l"=> 6,
                " "=> 7,
                "m~n"=> 8)

    @test doc[j"#"] == doc
    @test doc[j"#/foo"] == ["bar", "baz"]
    @test doc[JSONPointer.Pointer("#/foo/0"; shift_index = true)] == "bar"
    @test doc[j"#/foo/1"] == "bar"
    @test doc[j"#/"] == 0
    @test doc[j"#/a~1b"] == 1
    @test doc[j"#/c%25d"] == 2
    @test doc[j"#/e%5Ef"] == 3
    @test doc[j"#/g%7Ch"] == 4
    @test doc[j"#/i%5Cj"] == 5
    @test doc[j"#/k%22l"] == 6
    @test doc[j"#/%20"] == 7
    @test doc[j"#/m~0n"] == 8
end

@testset "WrongInputTests" begin

    @test_throws ArgumentError JSONPointer.Pointer("some/thing")

    doc = [0, 1, 2]
    @test_throws ArgumentError doc[j"/a"]
    @test_throws BoundsError doc[j"/10"]

end

@testset "JSONPointer Advanced" begin
    a = j"/a/2/d::Vector"
    b = j"/a/2/e::Vector{Int}"
    c = j"/a/2/f::Vector{Float64}"

    @test a.token == ("a", 2, "d")
    @test eltype(a) <: Array
    @test eltype(b) <: Array{Int, 1}
    @test eltype(c) <: Array{Float64, 1}
end

@testset "construct Dict with JSONPointer" begin
    p1 = j"/a/1/b"
    p2 = j"/cd/2/ef"

    data = Dict(p1 =>1, p2 => 2)
    @test data[p1] == 1
    @test data[p2] == 2
    @test haskey(data, p1)
    @test haskey(data, p2)
    @test !haskey(data, j"/x")
    @test !haskey(data, j"/ba/5")

    p1 = j"/ab/1"
    p2 = j"/cd/2/ef"

    data = OrderedDict(p1 => "This", p2 => "Is my Data")
    @test data[p1] == "This"
    @test data[p2] == "Is my Data"

end

@testset "access deep nested object" begin
    data = [Dict("a" => 10)]
    @test data[j"/1/a"] == 10

    p1 = j"/a/b/c/d/e/f/g/1/2/a/b/c"
    data = Dict(p1 => "sooo deep")
    @test data[p1] == "sooo deep"
    @test get(data, p1, missing) == "sooo deep"

    @test haskey(data, j"/a/b/c")
    @test haskey(data, j"/a/b/c/d")
    @test haskey(data, j"/a/b/c/d/e/f/g/1")
    @test !haskey(data, j"/a/b/c/d/e/f/g/2")

    @test isa(data[j"/a/b/c"], AbstractDict)
    @test isa(data[j"/a/b/c/d"], AbstractDict)
    @test isa(data[j"/a/b/c/d/e"], AbstractDict)
    @test isa(data[j"/a/b/c/d/e/f/g"], Array)
    @test ismissing(data[j"/a/b/c/d/e/f/g/1/1"])

    @test get(data, j"/a/f", missing) |> ismissing
    @test get(data, j"/a/b/c/d/e/f/g/5", 10000) == 10000

    @test_throws KeyError data[j"/a/f"]
    @test_throws KeyError data[j"/x"]
    @test_throws BoundsError data[j"/a/b/c/d/e/f/g/5"]

    data = [[10, 20, 30, ["me"]]]
    @test data[j"/1"] == [10, 20, 30, ["me"]]
    @test data[j"/1/2"] == 20
    @test data[j"/1/4"] == ["me"]
    @test data[j"/1/4/1"] == "me"

    # need to add get for Array?
    @test_broken get(data, j"/1", missing) |> ismissing
end

@testset "grow obejct and array" begin
    d = Dict(j"/a" => Dict())
    d[j"/a/b"] = []
    d[j"/a/b/2"] = 1
    d[j"/a/b/5"] = 2
    @test_throws Exception d[j"/a/5"] = "something"
    @test_throws Exception d[j"/a/b/gd"] = "nothing"

    @test isa(d[j"/a/b"], Array)
    @test isa(d[j"/a/b/1"], Missing)
end

@testset "Enforce type on JSONPointer" begin
    p1 = j"/a/1::String"
    p2 = j"/a/2::Int"

    data = Dict(p1 =>"this", p2 => 20)
    @test data[p1] == "this"
    @test data[p2] == 20

    @test_throws ErrorException data[p1] = 20
    @test_throws ErrorException data[p2] = "this"

    p1 = j"/a/1::Vector{Int}"
    p2 = j"/a/2::Float64"
    p3 = j"/b::Vector{String}"

    data = Dict(p1 =>[10], p2 => 20,  p3 => ["this", "is"])
    @test data[p1] == Int[10]
    @test data[p2] == 20.
    @test data[p3] == ["this", "is"]

    @test_throws ErrorException data[p2] = [1000]
    @test_throws ErrorException data[p3] = "this"

    d = Dict(p1 => missing, p3 => missing)
    @test d[p1] == JSONPointer.null_value(p1)
    @test d[p3] == JSONPointer.null_value(p3)

    # TODO User Defined type
    struct Foo
    end
    # @test_broken j"/foo::Foo"

end

@testset "literal string for a Number" begin
    p1 = j"/\5"
    p2 = j"/\559"
    p3 = j"/\900/10"

    d = Dict(p1 => 1, p2 => 2, p3 => 3)
    @test d[p1] == 1
    @test d["5"] == 1
    @test d[p2] == 2
    @test d["559"] == 2

    @test d[p3] == 3
    @test isa(d["900"], Array)
end


text:
module Unityper

@nospecialize
function expr_to_type(mod::Module, typ)
    typ isa Symbol && return getproperty(mod, typ)
    typ isa Expr || error("oof: $typ")
    base = typ.args[1]
    if base isa Symbol
        baset = getproperty(mod, base)
    else
        @assert Meta.isexpr(base, :curly)
        baset = expr_to_type(mod, base)
    end
    curlytypes = Vector{Any}(undef, length(typ.args)-1)
    for i ∈ eachindex(curlytypes)
        c = typ.args[1+i]
        if c isa Symbol
            curlytypes[i] = getproperty(mod, c)
        elseif Meta.isexpr(base, :curly)
            curlytypes[i] = expr_to_type(mod, c)
        else
            @assert isbitstype(c)
            curlytypes[i] = c
        end
    end
    baset{curlytypes...}
end
@specialize

export @compactify, @compactified
include("compactify.jl")
include("compactified.jl")

end


text:
############################################################
# joCoreBlock ##############################################
############################################################

############################################################
## outer constructors

"""
    julia> op = joCoreBlock(op1[,op2][,...];
        [moffsets=...,][noffsets=...,][weights=...,]
        [mextend=...,][nextend=...,][name=...])

Universal (Core) block operator composed from serial JOLI operators

# Signature

    joCoreBlock(ops::joAbstractLinearOperator...;
        moffsets::LocalVector{Integer},noffsets::LocalVector{Integer},
        weights::LocalVector,mextend::Integer,nextend::Integer,name::String)

# Arguments

- `op#`: JOLI operators (subtypes of joAbstractLinearOperator)
- keywords
    - `moffsets`: vector of starting indecies in 1st dimansion
    - `noffsets`: vector of starting indecies in 2nd dimansion
    - `weights`: vector of waights for each operator
    - `mextend`: size of zero extension in 1st dimension
    - `nextend`: size of zero extension in 2nd dimension
    - `name`: custom name

# Notes

- all given operators must have same domain/range types
- the domain/range types of joCoreBlock are equal to domain/range types of the given operators

# Example

define operators

    a=rand(ComplexF64,4,5);
    A=joMatrix(a;DDT=ComplexF32,RDT=ComplexF64,name="A")

    b=rand(ComplexF64,7,8);
    B=joMatrix(b;DDT=ComplexF32,RDT=ComplexF64,name="B")

    c=rand(ComplexF64,6,8);
    C=joMatrix(c;DDT=ComplexF32,RDT=ComplexF64,name="C")

define supporting vectors for more complex blocks

    moff=[0;5;13]
    noff=[0;6;15]
    w=rand(ComplexF64,3)

basic diagonal-corners adjacent blocks

    BD=joCoreBlock(A,B,C)

sparse blocks

    BD=joCoreBlock(A,B,C;moffsets=moff,noffsets=noff)

sparse blocks with zero extansion of (mextend,nextend) size

    BD=joCoreBlock(A,B,C;moffsets=moff,noffsets=noff,mextend=5,nextend=5)

weighted basic diagonal-corners adjacent blocks

    BD=joCoreBlock(A,B,C;weights=w)

"""
function joCoreBlock(ops::joAbstractLinearOperator...;
        moffsets::LocalVector{OT}=zeros(Int,0),noffsets::LocalVector{OT}=zeros(Int,0),
        weights::LocalVector{WT}=zeros(0),mextend::Integer=0,nextend::Integer=0,
        name::String="joCoreBlock") where {OT<:Integer,WT<:Number}

    isempty(ops) && throw(joCoreBlockException("empty argument list"))
    l=length(ops)
    for i=1:l
        deltype(ops[i])==deltype(ops[1]) || throw(joCoreBlockException("domain type mismatch for $i operator"))
        reltype(ops[i])==reltype(ops[1]) || throw(joCoreBlockException("range type mismatch for $i operator"))
    end
    mo=Base.deepcopy(moffsets)
    (length(mo)==l || length(mo)==0) || throw(joCoreBlockException("lenght of moffsets vector does not match number of operators"))
    no=Base.deepcopy(noffsets)
    (length(no)==l || length(no)==0) || throw(joCoreBlockException("lenght of noffsets vector does not match number of operators"))
    ws=Base.deepcopy(weights)
    (length(ws)==l || length(ws)==0) || throw(joCoreBlockException("lenght of weights vector does not match number of operators"))
    mextend>=0 || throw(joCoreBlockException("mextend must be >=0"))
    nextend>=0 || throw(joCoreBlockException("nextend must be >=0"))
    ms=zeros(Int,l)
    ns=zeros(Int,l)
    for i=1:l
        ms[i]=ops[i].m
        ns[i]=ops[i].n
    end
    if (length(mo)==0 && length(no)==0)
        m=sum(ms)
        n=sum(ns)
        mo=zeros(Int,l)
        no=zeros(Int,l)
        mo[1]=0
        no[1]=0
        for i=2:l
            mo[i]=mo[i-1]+ms[i-1]
            no[i]=no[i-1]+ns[i-1]
        end
    elseif (length(mo)==l && length(no)==0)
        m=max((mo+ms)...)
        n=max(ns...)
        no=zeros(Int,l)
    elseif (length(mo)==0 && length(no)==l)
        m=max(ms...)
        n=max((no+ns)...)
        mo=zeros(Int,l)
    else
        m=max((mo+ms)...)
        n=max((no+ns)...)
    end
    m+=mextend
    n+=nextend
    weighted=(length(ws)==l)
    fops=Vector{joAbstractLinearOperator}(undef,0)
    fops_T=Vector{joAbstractLinearOperator}(undef,0)
    fops_A=Vector{joAbstractLinearOperator}(undef,0)
    fops_C=Vector{joAbstractLinearOperator}(undef,0)
    for i=1:l
        if weighted
            push!(fops,ws[i]*ops[i])
            push!(fops_T,ws[i]*transpose(ops[i]))
            push!(fops_A,conj(ws[i])*adjoint(ops[i]))
            push!(fops_C,conj(ws[i])*conj(ops[i]))
        else
            push!(fops,ops[i])
            push!(fops_T,transpose(ops[i]))
            push!(fops_A,adjoint(ops[i]))
            push!(fops_C,conj(ops[i]))
        end
    end
    return joCoreBlock{deltype(fops[1]),reltype(fops[1])}(name*"($l)",m,n,l,ms,ns,mo,no,ws,
                      fops,fops_T,fops_A,fops_C,@joNF,@joNF,@joNF,@joNF)
end

############################
## overloaded Base functions

# show(jo)
function show(A::joCoreBlock)
    println("Type: $(typeof(A).name)")
    println("Name: $(A.name)")
    println("Size: $(size(A))")
    println(" DDT: $(deltype(A))")
    println(" RDT: $(reltype(A))")
    println("# of ops: $(A.l)")
    for i=1:A.l
        a=A.fop[i]
        println("  $i Type: $(typeof(a).name)")
        println("     Name: $(a.name)")
        println("     Size: $(size(a))")
        println(" M-offset: $(A.mo[i])")
        println(" N-offset: $(A.no[i])")
        try
            println("   Weight: $(A.ws[i])")
        catch
            nothing
        end
        println("      DDT: $(deltype(a))")
        println("      RDT: $(reltype(a))")
    end
end

# conj(jo)
conj(A::joCoreBlock{DDT,RDT}) where {DDT,RDT} =
    joCoreBlock{DDT,RDT}("(conj("*A.name*"))",
        A.m,A.n,A.l,A.ms,A.ns,A.mo,A.no,A.ws,
        A.fop_C,A.fop_A,A.fop_T,A.fop,
        A.iop_C,A.iop_A,A.iop_T,A.iop)

# transpose(jo)
transpose(A::joCoreBlock{DDT,RDT}) where {DDT,RDT} =
    joCoreBlock{RDT,DDT}("(transpose("*A.name*"))",
        A.n,A.m,A.l,A.ns,A.ms,A.no,A.mo,A.ws,
        A.fop_T,A.fop,A.fop_C,A.fop_A,
        A.iop_T,A.iop,A.iop_C,A.iop_A)

# adjoint(jo)
adjoint(A::joCoreBlock{DDT,RDT}) where {DDT,RDT} =
    joCoreBlock{RDT,DDT}("(adjoint("*A.name*"))",
        A.n,A.m,A.l,A.ns,A.ms,A.no,A.mo,A.ws,
        A.fop_A,A.fop_C,A.fop,A.fop_T,
        A.iop_A,A.iop_C,A.iop,A.iop_T)

# *(jo,vec)
function *(A::joCoreBlock{ADDT,ARDT},v::LocalVector{ADDT}) where {ADDT,ARDT}
    size(A,2) == size(v,1) || throw(joCoreBlockException("shape mismatch"))
    V=zeros(ARDT,A.m) # must be preallocated with zeros
    for i=1:1:A.l
        sm=A.mo[i]+1
        em=A.mo[i]+A.ms[i]
        sn=A.no[i]+1
        en=A.no[i]+A.ns[i]
        V[sm:em]+=A.fop[i]*v[sn:en]
    end
    return V
end

# *(jo,mvec)
function *(A::joCoreBlock{ADDT,ARDT},mv::LocalMatrix{ADDT}) where {ADDT,ARDT}
    size(A, 2) == size(mv, 1) || throw(joCoreBlockException("shape mismatch"))
    MV=zeros(ARDT,size(A,1),size(mv,2)) # must be preallocated with zeros
    for i=1:1:A.l
        sm=A.mo[i]+1
        em=A.mo[i]+A.ms[i]
        sn=A.no[i]+1
        en=A.no[i]+A.ns[i]
        MV[sm:em,:]+=A.fop[i]*mv[sn:en,:]
    end
    return MV
end

# -(jo)
function -(A::joCoreBlock{DDT,RDT}) where {DDT,RDT}
    fops=Vector{joAbstractLinearOperator}(undef,0)
    fops_T=Vector{joAbstractLinearOperator}(undef,0)
    fops_A=Vector{joAbstractLinearOperator}(undef,0)
    fops_C=Vector{joAbstractLinearOperator}(undef,0)
    for i=1:A.l
        push!(fops,-A.fop[i])
        push!(fops_T,-A.fop_T[i])
        push!(fops_A,-A.fop_A[i])
        push!(fops_C,-A.fop_C[i])
    end
    return joCoreBlock{DDT,RDT}("(-"*A.name*")",
        A.m,A.n,A.l,A.ms,A.ns,A.mo,A.no,A.ws,
        fops,fops_T,fops_A,fops_C,
        A.iop,A.iop_T,A.iop_A,A.iop_C)
end


text:
#!/usr/bin/env julia

module PE17

const ones=["one","two","three","four","five","six","seven","eight","nine",
    "ten","eleven","twelve","thirteen","fourteen","fifteen","sixteen","seventeen","eighteen","nineteen"];
const tens=["","twenty","thirty","forty","fifty","sixty","seventy","eighty","ninety"];

function toEnglish(n::Int)
    @assert 0<n<1000000 "$(n) is out of range!"
    if n<20
        ones[n]
    elseif n<100
        tens[n÷10]*(n%10≠0 ? ones[n%10] : "")
    elseif n<1000
        ones[n÷100]*"hundred"*(n%100≠0 ? "and"*toEnglish(n%100) : "")
    else
        toEnglish(n÷1000)*"thousand"*(n%1000≠0 ? toEnglish(n%1000) : "")
    end
end

end

export pe17

pe17(n::Integer=1000)=join(PE17.toEnglish.(1:n))|>length

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe17()|>print
end

max_stars_repo_path: 17. Number letter counts/17.jl

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_path: 17. Number letter counts/17.jl

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_path: 17. Number letter counts/17.jl

max_forks_repo_name: miRoox/ProjectEuler


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.
# *********************************************************************************
"""
    struct ElectricTariff

- data for electric tariff in reopt model
- can be defined using custom rates or URDB rate
- very similar to the URDB struct but includes export rates and bins
"""
struct ElectricTariff
    energy_rates::AbstractArray{Float64, 2} # gets a second dim with tiers
    energy_tier_limits::AbstractArray{Float64,1}
    n_energy_tiers::Int

    monthly_demand_rates::AbstractArray{Float64, 2} # gets a second dim with tiers
    time_steps_monthly::AbstractArray{AbstractArray{Int64,1},1}  # length = 0 or 12
    monthly_demand_tier_limits::AbstractArray{Float64,1}
    n_monthly_demand_tiers::Int

    tou_demand_rates::AbstractArray{Float64, 2} # gets a second dim with tiers
    tou_demand_ratchet_time_steps::AbstractArray{AbstractArray{Int64,1},1}  # length = n_tou_demand_ratchets
    tou_demand_tier_limits::AbstractArray{Float64,1}
    n_tou_demand_tiers::Int

    demand_lookback_months::AbstractArray{Int,1}
    demand_lookback_percent::Float64
    demand_lookback_range::Int

    fixed_monthly_charge::Float64
    annual_min_charge::Float64
    min_monthly_charge::Float64

    export_rates::Dict{Symbol, AbstractArray}
    export_bins::AbstractArray{Symbol,1}

    coincident_peak_load_active_time_steps::AbstractVector{AbstractVector{Int64}}
    coincident_peak_load_charge_per_kw::AbstractVector{Float64}
    coincpeak_periods::AbstractVector{Int64}
end

"""
    ElectricTariff

ElectricTariff constructor
```julia
function ElectricTariff(;
    urdb_label::String="",
    urdb_response::Dict=Dict(),
    urdb_utility_name::String="",
    urdb_rate_name::String="",
    year::Int=2020,
    time_steps_per_hour::Int=1,
    NEM::Bool=false,
    wholesale_rate::T1=nothing,
    export_rate_beyond_net_metering_limit::T2=nothing,
    monthly_energy_rates::Array=[],
    monthly_demand_rates::Array=[],
    blended_annual_energy_rate::S=nothing,
    blended_annual_demand_rate::R=nothing,
    add_monthly_rates_to_urdb_rate::Bool=false,
    tou_energy_rates_per_kwh::Array=[],
    add_tou_energy_rates_to_urdb_rate::Bool=false,
    remove_tiers::Bool=false,
    demand_lookback_months::AbstractArray{Int64, 1}=Int64[],
    demand_lookback_percent::Float64=0.0,
    demand_lookback_range::Int=0,
    coincident_peak_load_active_time_steps::Vector{Vector{Int64}}=[Int64[]],
    coincident_peak_load_charge_per_kw::AbstractVector{<:Real}=Real[]
    ) where {
        T1 <: Union{Nothing, Int, Float64, Array},
        T2 <: Union{Nothing, Int, Float64, Array},
        S <: Union{Nothing, Int, Float64},
        R <: Union{Nothing, Int, Float64}
    }
```

!!! note
    The `NEM` boolean is determined by the ElectricUtility.net_metering_limit_kw. There is no need to pass in a `NEM`
    value.

"""
function ElectricTariff(;
    urdb_label::String="",
    urdb_response::Dict=Dict(),
    urdb_utility_name::String="",
    urdb_rate_name::String="",
    year::Int=2020,
    time_steps_per_hour::Int=1,
    NEM::Bool=false,
    wholesale_rate::T1=nothing,
    export_rate_beyond_net_metering_limit::T2=nothing,
    monthly_energy_rates::Array=[],
    monthly_demand_rates::Array=[],
    blended_annual_energy_rate::S=nothing,
    blended_annual_demand_rate::R=nothing,
    add_monthly_rates_to_urdb_rate::Bool=false,
    tou_energy_rates_per_kwh::Array=[],
    add_tou_energy_rates_to_urdb_rate::Bool=false,
    remove_tiers::Bool=false,
    demand_lookback_months::AbstractArray{Int64, 1}=Int64[],
    demand_lookback_percent::Float64=0.0,
    demand_lookback_range::Int=0,
    coincident_peak_load_active_time_steps::Vector{Vector{Int64}}=[Int64[]],
    coincident_peak_load_charge_per_kw::AbstractVector{<:Real}=Real[]
    ) where {
        T1 <: Union{Nothing, Int, Float64, Array},
        T2 <: Union{Nothing, Int, Float64, Array},
        S <: Union{Nothing, Int, Float64},
        R <: Union{Nothing, Int, Float64}
    }
    # TODO remove_tiers for multinode models
    nem_rate = Float64[]

    energy_tier_limits = Float64[]
    n_energy_tiers = 1
    monthly_demand_tier_limits = Float64[]
    n_monthly_demand_tiers = 1
    tou_demand_tier_limits = Float64[]
    n_tou_demand_tiers = 1
    time_steps_monthly = get_monthly_time_steps(year, time_steps_per_hour=time_steps_per_hour)

    u = nothing
    if !isempty(urdb_response)

        u = URDBrate(urdb_response, year, time_steps_per_hour=time_steps_per_hour)

    elseif !isempty(urdb_label)

        u = URDBrate(urdb_label, year, time_steps_per_hour=time_steps_per_hour)

    elseif !isempty(urdb_utility_name) && !isempty(urdb_rate_name)

        u = URDBrate(urdb_utility_name, urdb_rate_name, year, time_steps_per_hour=time_steps_per_hour)

    elseif !isempty(tou_energy_rates_per_kwh) && length(tou_energy_rates_per_kwh) == 8760*time_steps_per_hour

        tou_demand_rates = Float64[]
        tou_demand_ratchet_time_steps = []
        energy_rates = tou_energy_rates_per_kwh
        monthly_demand_rates = convert(Array{Float64}, monthly_demand_rates)

        fixed_monthly_charge = 0.0
        annual_min_charge = 0.0
        min_monthly_charge = 0.0

        if NEM
            nem_rate = [-0.999 * x for x in energy_rates]
        end

    elseif !isempty(monthly_energy_rates)

        invalid_args = String[]
        if !(length(monthly_energy_rates) == 12)
            push!(invalid_args, "length(monthly_energy_rates) must equal 12, got length $(length(monthly_energy_rates))")
        end
        if !isempty(monthly_demand_rates) && !(length(monthly_demand_rates) == 12)
            push!(invalid_args, "length(monthly_demand_rates) must equal 12, got length $(length(monthly_demand_rates))")
        end
        if length(invalid_args) > 0
            error("Invalid argument values: $(invalid_args)")
        end

        if isempty(monthly_demand_rates)
            monthly_demand_rates = repeat([0.0], 12)
        end

        tou_demand_rates = Float64[]
        tou_demand_ratchet_time_steps = []
        energy_rates = Real[]
        for m in 1:12
            append!(energy_rates, [monthly_energy_rates[m] for ts in time_steps_monthly[m]])
        end

        fixed_monthly_charge = 0.0
        annual_min_charge = 0.0
        min_monthly_charge = 0.0

        if NEM
            nem_rate = [-0.999 * x for x in energy_rates]
        end

    elseif !isnothing(blended_annual_energy_rate)

        tou_demand_rates = Float64[]
        tou_demand_ratchet_time_steps = []
        energy_rates = repeat(Real[blended_annual_energy_rate], 8760 * time_steps_per_hour)
        if !isnothing(blended_annual_demand_rate)
            monthly_demand_rates = repeat(Real[blended_annual_demand_rate], 12)
        end
        if isempty(monthly_demand_rates)
            monthly_demand_rates = repeat([0.0], 12)
        end

        fixed_monthly_charge = 0.0
        annual_min_charge = 0.0
        min_monthly_charge = 0.0

        if NEM
            nem_rate = [-0.999 * x for x in energy_rates]
        end

    else
        error("Creating ElectricTariff requires at least urdb_label, urdb_response, monthly rates, annual rates, or tou_energy_rates_per_kwh.")
    end

    if !isnothing(u)  # use URDBrate

        if NEM
            t = get_tier_with_lowest_energy_rate(u)
            nem_rate = [-0.999 * x for x in u.energy_rates[t,:]]
        end

        energy_rates = u.energy_rates
        energy_tier_limits = u.energy_tier_limits
        n_energy_tiers = u.n_energy_tiers
        users_monthly_demand_rates = copy(monthly_demand_rates)
        monthly_demand_rates = u.monthly_demand_rates
        monthly_demand_tier_limits = u.monthly_demand_tier_limits
        n_monthly_demand_tiers = u.n_monthly_demand_tiers
        tou_demand_rates = u.tou_demand_rates
        tou_demand_tier_limits = u.tou_demand_tier_limits
        n_tou_demand_tiers = u.n_tou_demand_tiers

        if remove_tiers
            energy_rates, monthly_demand_rates, tou_demand_rates = remove_tiers_from_urdb_rate(u)
            energy_tier_limits, monthly_demand_tier_limits, tou_demand_tier_limits =
                Float64[], Float64[], Float64[]
            n_energy_tiers, n_monthly_demand_tiers, n_tou_demand_tiers = 1, 1, 1
        end

        tou_demand_ratchet_time_steps = u.tou_demand_ratchet_time_steps
        demand_lookback_months = u.demand_lookback_months
        demand_lookback_percent = u.demand_lookback_percent
        demand_lookback_range = u.demand_lookback_range
        fixed_monthly_charge = u.fixed_monthly_charge
        annual_min_charge = u.annual_min_charge
        min_monthly_charge = u.min_monthly_charge

        if add_monthly_rates_to_urdb_rate
            if length(monthly_energy_rates) == 12
                for tier in 1:size(energy_rates, 2), mth in 1:12, ts in time_steps_monthly[mth]
                    energy_rates[ts, tier] += monthly_energy_rates[mth]
                end
            end
            if length(users_monthly_demand_rates) == 12
                for tier in 1:size(monthly_demand_rates, 2), mth in 1:12
                    monthly_demand_rates[mth, tier] += users_monthly_demand_rates[mth]
                end
            end
        end

        if add_tou_energy_rates_to_urdb_rate && length(tou_energy_rates_per_kwh) == size(energy_rates, 1)
            for tier in 1:size(energy_rates, 2)
                energy_rates[1:end, tier] += tou_energy_rates_per_kwh
            end
        end
    else
        # need to reshape cost vectors to arrays (2nd dim is for tiers)
        energy_rates = reshape(energy_rates, :, 1)
        monthly_demand_rates = reshape(monthly_demand_rates, :, 1)
        tou_demand_rates = reshape(tou_demand_rates, :, 1)
    end

    #= export_rates
    There are three Export tiers and their associated export rates (negative values):
    1. NEM (Net Energy Metering)
    2. WHL (Wholesale)
    3. EXC (Excess, beyond NEM)

    Only one of NEM and Wholesale can be exported into due to the binary constraints.
    Excess can be exported into in the same time step as NEM.

    Excess is meant to be combined with NEM: NEM export is limited to the total grid purchased energy in a year and some
    utilities offer a compensation mechanism for export beyond the site load.
    The Excess tier is not available with the Wholesale tier.

    - if NEM then set ExportRate[:Nem, :] to energy_rate[tier_with_lowest_energy_rate, :]
    - user can provide either scalar wholesale rate or vector of time_steps,
    =#
    whl_rate = create_export_rate(wholesale_rate, length(energy_rates), time_steps_per_hour)
    if !isnothing(u) && sum(u.sell_rates) < 0
        whl_rate += u.sell_rates
    end
    exc_rate = create_export_rate(export_rate_beyond_net_metering_limit, length(energy_rates), time_steps_per_hour)

    if !NEM & (sum(whl_rate) >= 0)
        export_rates = Dict{Symbol, AbstractArray}()
        export_bins = Symbol[]
    elseif !NEM
        export_bins = [:WHL]
        export_rates = Dict(:WHL => whl_rate)
    elseif (sum(whl_rate) >= 0)
        export_bins = [:NEM]
        export_rates = Dict(:NEM => nem_rate)
        if sum(exc_rate) < 0
            push!(export_bins, :EXC)
            export_rates[:EXC] = exc_rate
        end
    else
        export_bins = [:NEM, :WHL]
        export_rates = Dict(:NEM => nem_rate, :WHL => whl_rate)
        if sum(exc_rate) < 0
            push!(export_bins, :EXC)
            export_rates[:EXC] = exc_rate
        end
    end

    coincpeak_periods = Int64[]
    if !isempty(coincident_peak_load_charge_per_kw)
        coincpeak_periods = collect(1:length(coincident_peak_load_charge_per_kw))
    end

    ElectricTariff(
        energy_rates,
        energy_tier_limits,
        n_energy_tiers,
        monthly_demand_rates,
        time_steps_monthly,
        monthly_demand_tier_limits,
        n_monthly_demand_tiers,
        tou_demand_rates,
        tou_demand_ratchet_time_steps,
        tou_demand_tier_limits,
        n_tou_demand_tiers,
        demand_lookback_months,
        demand_lookback_percent,
        demand_lookback_range,
        fixed_monthly_charge,
        annual_min_charge,
        min_monthly_charge,
        export_rates,
        export_bins,
        coincident_peak_load_active_time_steps,
        coincident_peak_load_charge_per_kw,
        coincpeak_periods
    )
end

function get_tier_with_lowest_energy_rate(u::URDBrate)
    """
    ExportRate should be lowest energy cost for tiered rates.
    Otherwise, ExportRate can be > FuelRate, which leads REopt to export all PV energy produced.
    """
    tier_with_lowest_energy_cost = 1
    if length(u.energy_tier_limits) > 1
        annual_energy_charge_sums = Float64[]
        for etier in u.energy_rates
            push!(annual_energy_charge_sums, sum(etier))
        end
        tier_with_lowest_energy_cost =
            findall(annual_energy_charge_sums .== minimum(annual_energy_charge_sums))[1]
    end
    return tier_with_lowest_energy_cost
end

"""
    function create_export_rate(e::Nothing, N::Int, ts_per_hour::Int=1)
No export rate provided by user: set to 0 dollars/kWh for all time
"""
function create_export_rate(e::Nothing, N::Int, ts_per_hour::Int=1)
    [0 for _ in range(1, stop=N) for ts in 1:ts_per_hour]
end

"""
    function create_export_rate(e::T, N::Int, ts_per_hour::Int=1) where T<:Real
Case for scaler export rate provided -> convert to array of time_steps
"""
function create_export_rate(e::T, N::Int, ts_per_hour::Int=1) where T<:Real
    repeat(float(-1*e), N * ts_per_hour)
end

"""
    function create_export_rate(e::AbstractArray{<:Real, 1}, N::Int, ts_per_hour::Int=1)

Check length of e and upsample if length(e) != N
"""
function create_export_rate(e::AbstractArray{<:Real, 1}, N::Int, ts_per_hour::Int=1)
    Ne = length(e)
    if Ne != Int(N/ts_per_hour) || Ne != N
        @error "Export rates do not have correct number of entries. Must be $(N) or $(Int(N/ts_per_hour))."
    end
    if Ne != N  # upsample
        export_rates = [-1*x for x in e for ts in 1:ts_per_hour]
    else
        export_rates = -1*e
    end
    return export_rates
end

"""
    get_monthly_time_steps(year::Int; time_steps_per_hour=1)

return Array{Array{Int64,1},1}, size = (12,)
"""
function get_monthly_time_steps(year::Int; time_steps_per_hour=1)
    a = Array[]
    i = 1
    for m in range(1, stop=12)
        n_days = daysinmonth(Date(string(year) * "-" * string(m)))
        stop = n_days * 24 * time_steps_per_hour + i - 1
        if m == 2 && isleapyear(year)
            stop -= 24 * time_steps_per_hour  # TODO support extra day in leap years?
        end
        steps = [step for step in range(i, stop=stop)]
        append!(a, [steps])
        i = stop + 1
    end
    return a
end

# TODO use this function only for URDBrate
function remove_tiers_from_urdb_rate(u::URDBrate)
    # tariff args: have to validate that there are no tiers
    if length(u.energy_tier_limits) > 1
        @warn "Energy rate contains tiers. Using the first tier!"
    end
    elec_rates = vec(u.energy_rates[1,:])

    if u.n_monthly_demand_tiers > 1
        @warn "Monthly demand rate contains tiers. Using the last tier!"
    end
    if u.n_monthly_demand_tiers > 0
        demand_rates_monthly = vec(u.monthly_demand_rates[:,u.n_monthly_demand_tiers])
    else
        demand_rates_monthly = vec(u.monthly_demand_rates)  # 0×0 Array{Float64,2}
    end

    if u.n_tou_demand_tiers > 1
        @warn "TOU demand rate contains tiers. Using the last tier!"
    end
    if u.n_tou_demand_tiers > 0
        demand_rates = vec(u.tou_demand_rates[:,u.n_tou_demand_tiers])
    else
        demand_rates = vec(u.tou_demand_rates)
    end

    return elec_rates, demand_rates_monthly, demand_rates
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
const Cfunction = Ptr{Cvoid}
const IBString = ConstCxxRef{StdString}

struct JLClient
    managedAccounts::Cfunction
    accountSummary::Cfunction
    accountSummaryEnd::Cfunction
end

struct JLComboLeg
    conId::CxxLong
    ratio::CxxLong
    action::String

    exchange::String
    openClose::CxxLong

    shortSaleSlot::CxxLong
    designatedLocation::String
    exemptCode::CxxInt
end


text: base(::Type{UInt64}) = zero(UInt64)

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using MLStyle

export canonical
using Lazy
canonical(x) = x

function canonical(expr :: Expr)
    # @show expr
    r = canonical
    @match expr begin
        :($x |> $f) => begin
            rf = r(f)
            rx = r(x)
            :($rf($rx)) |> r
        end

        Expr(:do, :($g($x)), :($f)) => begin
            rg = r(g)
            rf = r(f)
            rx = r(x)
            :($rg($rf, $rx)) |> r
        end

        :((iid($n))($dist)) => begin
            rn = r(n)
            rdist = r(dist)
            :(iid($rn,$rdist)) |> r
        end

        :($f($(args...))) => begin
            rf = r(f)
            rx = map(r,args)
            :($rf($(rx...)))
        end

        :(iid($n)($dist)) => begin
            rn = r(n)
            rdist = r(dist)
            :(iid($rn,$rdist)) |> r
        end

        x => x
    end
end

function canonical(m :: Model)
    @as x m begin
        convert.(Expr, x.body)
        canonical.(x)
        convert.(Statement, x)
        Model(m.args, x)
    end
end

ex1 = :(map(1:10) do x x^2 end)

# canonical(ex1)
# canonical(linReg1D)


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct PSSSimple <: PSS
        K_ω::Float64
        K_p::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of a PSS that returns a proportional droop voltage to add to the reference for the AVR

# Arguments
- `K_ω::Float64`: Proportional gain for frequency, validation range: (0, nothing)
- `K_p::Float64`: Proportional gain for active power, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct PSSSimple <: PSS
    "Proportional gain for frequency"
    K_ω::Float64
    "Proportional gain for active power"
    K_p::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function PSSSimple(K_ω, K_p, ext=Dict{String, Any}(), )
    PSSSimple(K_ω, K_p, ext, Vector{Symbol}(), 0, InfrastructureSystemsInternal(), )
end

function PSSSimple(; K_ω, K_p, ext=Dict{String, Any}(), )
    PSSSimple(K_ω, K_p, ext, )
end

# Constructor for demo purposes; non-functional.
function PSSSimple(::Nothing)
    PSSSimple(;
        K_ω=0,
        K_p=0,
        ext=Dict{String, Any}(),
    )
end

"""Get PSSSimple K_ω."""
get_K_ω(value::PSSSimple) = value.K_ω
"""Get PSSSimple K_p."""
get_K_p(value::PSSSimple) = value.K_p
"""Get PSSSimple ext."""
get_ext(value::PSSSimple) = value.ext
"""Get PSSSimple states."""
get_states(value::PSSSimple) = value.states
"""Get PSSSimple n_states."""
get_n_states(value::PSSSimple) = value.n_states
"""Get PSSSimple internal."""
get_internal(value::PSSSimple) = value.internal


text:
mvn = Sys.iswindows() ? "mvn.cmd" : "mvn"
which = Sys.iswindows() ? "where" : "which"

try
    run(`$which $mvn`)
catch
    error("Cannot find maven. Is it installed?")
end

SPARK_VERSION = get(ENV, "BUILD_SPARK_VERSION", "2.1.0")

cd(joinpath(dirname(@__DIR__), "jvm/sparkjl")) do
    run(`$mvn clean package -Dspark.version=$SPARK_VERSION`)
end


text:
using KissThreading
using Compat
using Compat.Test
using Compat.Random
using Compat.Random: GLOBAL_RNG
using Compat.Statistics: std, mean

include("bootstrap.jl")
include("bubblesort.jl")
include("sort_batch.jl")
include("summation.jl")
include("mapreduce.jl")

println("========")


text:
################################
#  DistrustSociety Type   #
################################

# struct DistrustSociety{N, K, T <: Real}
#     """
#     ------
#     w0: array with shape (N, K) - opinion vector with
#     dimension K for N agents
#     C0: array with shape (N, K, K) - opinion uncertainty
#     for each agent
#     mu0: array with shape (N, N) - mistrust array for each agent
#     s20: array with shape (N, N) - mistrust uncertainty for each agent
#     """

"""
    type DistrustSociety{N, K, T}

Type representing a generic Bayesian Agent Society

### Fields
* links (Matrix{T}): N×N-matrix that describes the interactions of the agents of the society
* agents (Vector{MoralVector}): N-vector of the MoralVector{K, T}s of the DistrustSociety
* C (Vector{Matrix{Float64}}):
* mu (Matrix{Float64}):
* s2 (Matrix{Float64}):
"""
struct DistrustSociety{N, K, T <: Real} <: Society{N, K, T}
    links::Matrix{Bool}
    agents::Vector{MoralVector{K, T}}
    C::Vector{Matrix{Float64}}
    mu::Matrix{Float64}
    s2::Matrix{Float64}
end

"""
    DistrustSociety()

Construct a random DistrustSociety with default size NSOC
The agents have the default number of components (KMORAL)
"""
DistrustSociety() =
    DistrustSociety{NSOC, KMORAL, Float64}(Matrix{Bool}(1 - eye(NSOC)), MoralVector{KMORAL, Float64}[MoralVector() for i in 1:NSOC],
                                                [λDEF*eye(KMORAL) for i in 1:NSOC], zeros(NSOC, NSOC), s2DEF*ones(NSOC, NSOC))

"""
    DistrustSociety(n::Int)

Construct a random DistrustSociety with size n
The agents have the default number of components (KMORAL)
"""
DistrustSociety(n::Int; λ::Float64 = λDEF, μ::Float64 = μDEF, s2::Float64 = s2DEF) =
    DistrustSociety{n, KMORAL, Float64}(Matrix{Bool}(1 - eye(n)), MoralVector{KMORAL, Float64}[MoralVector() for i in 1:n],
                                             [λ*eye(KMORAL) for i in 1:n], μ*ones(n, n), s2*ones(n, n))

"""
    DistrustSociety(agentArray::Vector{MoralVector{K, T}})

Construct a DistrustSociety with agents agentArray
"""
DistrustSociety{K, T}(agentArray::Vector{MoralVector{K, T}}; λ::Float64 = λDEF, μ::Float64 = μDEF, s2::Float64 = s2DEF) =
    DistrustSociety{length(agentArray), K, T}(Matrix{Bool}(1 - eye(length(agentArray))), deepcopy(agentArray),
                                                  [λ*eye(KMORAL) for i in 1:length(agentArray)], μ*ones(length(agentArray), length(agentArray)), s2*ones(length(agentArray), length(agentArray)))

########################################
#  DistrustSociety Redefinitions  #
########################################

function show(io::IO, soc::DistrustSociety)
    N, K = size(soc)
    print(io, N, "-sized DistrustSociety on a ", K, "-dimensional Moral space")
    #@printf io "ρ: %.4f\t ε: %.4f" soc.ρ soc.ε
end

########################################
#   DistrustSociety Definitions   #
########################################

"""
    epssoc(soc::DistrustSociety, i::Int, j::Int)

Computes `ε` agent `i` in a Society has about agent `j`
"""
epssoc(soc::DistrustSociety, i::Int, j::Int) = phi(soc.mu[i, j] / sqrt(1 + soc.s2[i, j]))

"""
    epssoc(soc::DistrustSociety)

Computes the matrix of distrust `ε` in a Society
"""
epssoc{N, K, T}(soc::DistrustSociety{N, K, T}) = [epssoc(soc, i, j) for i in 1:N, j in 1:N]

"Computes `γ` of agent `i` in a Society given MoralVector `x`"
gammasoc{N, K, T}(soc::DistrustSociety{N, K, T}, i::Int, x::MoralVector{K, T}) = sqrt(x[:]' * soc.C[i] * x[:])

"Computes `ρ` of agent `i` in a Society given MoralVector `x`"
rhosoc{N, K, T}(soc::DistrustSociety{N, K, T}, i::Int, x::MoralVector{K, T}) = rhosoc(gammasoc(soc, i, x))

"""
    cogcost{K,T}(soc::DistrustSociety{N, K, T}, i::Int, j::Int, x::MoralVector{K,T})

Cognitive cost MoralVector `soc[i]` suffers when learning MoralVector `soc[j]`'s opinion about MoralVector `x`
"""
function cogcost{N, K,T}(soc::DistrustSociety{N, K, T}, i::Int, j::Int, x::MoralVector{K,T})
    gam = gammasoc(soc, i, x)
    ε   = epssoc(soc, i, j)
    agi = soc[i]
    agj = soc[j]

    return - log( ε + (1 - 2ε) * phi(sign(agj ⋅ x) * (agi ⋅ x) / gam) )
end

# """
#     cogcost{K,T}(soc::DistrustSociety{N, K, T}, i::Int, j::Int, x::MoralVector{K,T})

# Cognitive cost MoralVector `agi` suffers when learning MoralVector `soc[i]`'s opinion about MoralVector `x`
# """
# function cogcost{N, K,T}(agi::MoralVector{K,T}, soc::DistrustSociety{N, K, T}, i::Int, j::Int, x::MoralVector{K,T})
#     gam = gammasoc(soc, i, x)
#     ε   = epssoc(soc, i, j)
#     agj = soc[j]

#     return - log( ε + (1 - 2ε) * phi(sign(agj ⋅ x) * (agi ⋅ x) / gam) )
# end

"Returns the two modulation functions: Fw and Fϵ"
function modfunc{N, K, T}(soc::DistrustSociety{N,K,T}, i::Int, j::Int, x::MoralVector{K,T})
    σ  = sign(soc[j] ⋅ x)
    h  = soc[i] ⋅ x
    γ  = gammasoc(soc, i, x)
    μ  = soc.mu[i, j]
    s2 = soc.s2[i, j]

    hσ1γ   = h*σ/γ
    μ1sqs2 = μ/sqrt(1 + s2)

    phiϵ = phi(μ1sqs2)
    phiw = phi(hσ1γ)
    Z    = phiϵ + phiw - 2*phiϵ*phiw

    Fw = (1 - 2phiϵ) * G(hσ1γ) / Z
    Fϵ = (1 - 2phiw) * G(μ1sqs2) / Z

    return (Fw, Fϵ)
end


text: print("HELLO")


text:
using CUDArt
using KUnet
using KUnet: accuracy
require(Pkg.dir("KUnet/test/mnist.jl"))

KUnet.ftype(Float32) # mnist has Float32 data
KUnet.atype(Array)
xtrn = MNIST.xtrn
xtst = MNIST.xtst
ytrn = MNIST.ytrn
ytst = MNIST.ytst
w0 = similar(ytst, size(ytst,1), 0)
d0 = 6
c0 = 1

if true
xnet = Op[Poly(c=c0,d=d0,w=w0), PercLoss()]
for i=1:1
    @time train(xnet, xtrn, ytrn; iters=100)
    println((i, size(xnet[1].s),
             accuracy(ytst, predict(xnet, xtst)),
             )) # accuracy(ytrn, predict(xnet, xtrn))))
end
end # if false

if true
strn = sparse(MNIST.xtrn)
stst = sparse(MNIST.xtst)
snet = Op[Poly(c=c0,d=d0,w=w0), PercLoss()]
for i=1:1
    @time train(snet, strn, ytrn; iters=100)
    println((i, size(snet[1].s),
             accuracy(ytst, predict(snet, stst)),
             )) # accuracy(ytrn, predict(snet, strn))))
end
end # if false

if false # mmul, hcat, ctranspose do not work
KUnet.atype(CudaArray)
cnet = Op[Poly(c=c0,d=d0,w=CudaArray(w0)), PercLoss()]
for i=1:1
    @time train(cnet, xtrn, ytrn; iters=100)
    println((i, size(cnet[1].s),
             accuracy(ytst, predict(cnet, xtst)),
             )) # accuracy(ytrn, predict(cnet, xtrn))))
end
KUnet.atype(Array)
end # if false


text:
module GOR

using Reexport
@reexport using Tables

include("abstractgoriter.jl")

include("gorfile.jl")
include("parquetfile.jl")
include("utils.jl")

include("verifyorder.jl")
include("select.jl")
include("filter.jl")
include("mutate.jl")
include("rename.jl")
include("merge.jl")
include("map.jl")
include("aggregate.jl")
include("groupby.jl")
include("join2.jl")

end # module


text:
"""
Vibration analysis of free-floating thin plate

Free-vibration problem is solved for a homogeneous free-floating
(unsupported) square plate. This is the NAFEMS Benchmark, Test No. FV12.

The plate is discretized with shell elements. Because no displacements are
prevented, the structure has six rigid body modes (six zero vibration
frequencies).

The nonzero benchmark frequencies are (in hertz):

1.622, 2.360, 2.922, 4.190, 4.190,  7.356, 7.356, 7.668.
"""
module barrel_w_stiffeners_examples

using LinearAlgebra
using Arpack
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using DataDrop: with_extension
# using KrylovKit
using FinEtools.MeshExportModule.VTKWrite: vtkwrite

function _execute_model(formul, input, visualize = true)
    E = 200e3*phun("MPa")
    nu = 0.3;
    rho= 7850*phun("KG/M^3");
    thickness = 2.0*phun("mm");
    OmegaShift = (100*2*pi)^2

    # Report
    @info "Mesh: $input"

    if !isfile(joinpath(dirname(@__FILE__()), input))
        success(run(`unzip -qq -d $(dirname(@__FILE__())) $(joinpath(dirname(@__FILE__()), "barrel_w_stiffeners-s3-mesh.zip"))`; wait = false))
    end
    output = FinEtools.MeshImportModule.import_H5MESH(joinpath(dirname(@__FILE__()), input))
    fens, fes  = output["fens"], output["fesets"][1]
    fens.xyz .*= phun("mm");

    # output = import_ABAQUS(joinpath(dirname(@__FILE__()), input))
    # fens = output["fens"]
    # fes = output["fesets"][1]

    # connected = findunconnnodes(fens, fes);
    # fens, new_numbering = compactnodes(fens, connected);
    # fes = renumberconn!(fes, new_numbering);

    # fens, fes = mergenodes(fens, fes, thickness/10)

    # FinEtools.MeshExportModule.H5MESH.write_H5MESH(with_extension(input, "h5mesh"), fens, fes)

    @info "Number of nodes and elements: $(count(fens)), $(count(fes))"

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    associate = formul.associategeometry!
    stiffness = formul.stiffness
    mass = formul.mass

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # l1 = collect(1:count(fens))
    # for i in [6]
    #     setebc!(dchi, l1, true, i)
    # end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    tim = @elapsed begin
        associategeometry!(femm, geom0)
    end
    @info "Associate Geometry time $tim [sec]"
    tim = @elapsed begin
        K = stiffness(femm, geom0, u0, Rfield0, dchi);
    end
    @info "Stiffness assembly time $tim [sec]"
    tim = @elapsed begin
        M = mass(femm, geom0, dchi);
    end
    @info "Mass assembly time $tim [sec]"

    # Solve
    neigvs = 40
    tim = @elapsed begin
        evals, evecs, convinfo = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM, explicittransform=:none)
        # evals, evecs, convinfo = geneigsolve((K+OmegaShift*M, M), neigvs,  :SR, krylovdim = 80)
    end
    @show convinfo
    evals[:] = evals .- OmegaShift;
    fs = real(sqrt.(complex(evals)))/(2*pi)
    @info "Frequencies: $(round.(fs[7:15], digits=4))"
    @info "Time $tim [sec]"

    # Visualization
    if visualize
        for ev in 7:30
            U = evecs[:, ev]
            scattersysvec!(dchi, 1.0/maximum(abs.(U)).*U)
            vtkwrite("$(input)-mode-$(ev).vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3]), ("ur", dchi.values[:, 4:6])])
            # update_rotation_field!(Rfield0, dchi)
            # plots = cat(plot_space_box([[0 0 -L/2]; [L/2 L/2 L/2]]),
            #     plot_nodes(fens),
            #     plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
            #     dims = 1)
            # pl = render(plots; title="$(ev)")
        end
    end

    return true
end

function test_convergence()
    input = "barrel_w_stiffeners-s3.h5mesh"
    formul = FEMMShellT3FFModule
    @info "Barrel With Stiffeners, free vibration, formulation=$(formul)"
    _execute_model(formul, input, true)
    return true
end

end # module

using .barrel_w_stiffeners_examples
m = barrel_w_stiffeners_examples
m.test_convergence()


text:
module Operators

export
    Δx, Δy, ΔzF, ΔzC,
    Δxᶜᶜᵃ, Δxᶠᶜᵃ, Δxᶠᶠᵃ, Δxᶜᶠᵃ,
    Δyᶜᶜᵃ, Δyᶠᶜᵃ, Δyᶠᶠᵃ, Δyᶜᶠᵃ,
    Axᵃᵃᶜ, Axᵃᵃᶠ, Axᶜᶜᶜ, Axᶠᶜᶜ, Axᶠᶠᶜ, Axᶠᶜᶠ,
    Ayᵃᵃᶜ, Ayᵃᵃᶠ, Ayᶜᶜᶜ, Ayᶜᶠᶜ, Ayᶠᶠᶜ, Ayᶜᶠᶠ,
    Azᵃᵃᵃ, Azᶠᶠᵃ, Azᶜᶜᵃ, Azᶠᶜᵃ, Azᶜᶠᵃ,
    Vᵃᵃᶜ, Vᵃᵃᶠ, Vᶜᶜᶜ, Vᶠᶜᶜ, Vᶜᶠᶜ, Vᶜᶜᶠ,
    δxᶜᵃᵃ, δxᶠᵃᵃ, δyᵃᶜᵃ, δyᵃᶠᵃ, δzᵃᵃᶜ, δzᵃᵃᶠ,
    Ax_ψᵃᵃᶠ, Ax_ψᵃᵃᶜ, Ay_ψᵃᵃᶠ, Ay_ψᵃᵃᶜ, Az_ψᵃᵃᵃ,
    Ax_uᶠᶜᶜ, Ay_vᶜᶠᶜ, Az_wᶜᶜᵃ,
    Ax_ζᶠᶠᶜ, Ay_ζᶠᶠᶜ, Ax_ηᶠᶜᶠ, Az_ηᶠᶜᵃ, Ay_ξᶜᶠᶠ, Az_ξᶜᶠᵃ,
    Az_wᶜᶜᵃ, Az_ηᶠᶜᵃ, Az_ξᶜᶠᵃ,
    Ax_cᶜᶜᶜ, Ay_cᶜᶜᶜ, Az_cᶜᶜᵃ,
    δᴶxᶜᵃᶜ, δᴶxᶜᵃᶠ, δᴶxᶠᵃᶜ, δᴶxᶠᵃᶠ, δᴶyᵃᶜᶜ, δᴶyᵃᶜᶠ, δᴶyᵃᶠᶜ, δᴶyᵃᶠᶠ, δᴶzᵃᵃᶜ, δᴶzᵃᵃᶠ,
    ∂xᶜᵃᵃ, ∂xᶠᵃᵃ, ∂yᵃᶜᵃ, ∂yᵃᶠᵃ, ∂zᵃᵃᶜ, ∂zᵃᵃᶠ,
    ∂xᶜᶜᵃ, ∂xᶠᶠᵃ, ∂xᶠᶜᵃ, ∂xᶜᶠᵃ, ∂yᶜᶜᵃ, ∂yᶠᶠᵃ, ∂yᶠᶜᵃ, ∂yᶜᶠᵃ,
    ∂²xᶜᵃᵃ, ∂²xᶠᵃᵃ, ∂²yᵃᶜᵃ, ∂²yᵃᶠᵃ, ∂²zᵃᵃᶜ, ∂²zᵃᵃᶠ,
    ∂⁴xᶜᵃᵃ, ∂⁴xᶠᵃᵃ, ∂⁴yᵃᶜᵃ, ∂⁴yᵃᶠᵃ, ∂⁴zᵃᵃᶜ, ∂⁴zᵃᵃᶠ,
    ℑxᶜᵃᵃ, ℑxᶠᵃᵃ, ℑyᵃᶜᵃ, ℑyᵃᶠᵃ, ℑzᵃᵃᶜ, ℑzᵃᵃᶠ,
    ℑᴶxᶜᵃᵃ, ℑᴶxᶠᵃᵃ, ℑᴶyᵃᶜᵃ, ℑᴶyᵃᶠᵃ, ℑᴶzᵃᵃᶜ, ℑᴶzᵃᵃᶠ,
    ℑxyᶜᶜᵃ, ℑxyᶠᶜᵃ, ℑxyᶠᶠᵃ, ℑxyᶜᶠᵃ, ℑxzᶜᵃᶜ, ℑxzᶠᵃᶜ, ℑxzᶠᵃᶠ, ℑxzᶜᵃᶠ, ℑyzᵃᶜᶜ, ℑyzᵃᶠᶜ, ℑyzᵃᶠᶠ, ℑyzᵃᶜᶠ,
    ℑxyzᶜᶜᶠ, ℑxyzᶜᶠᶜ, ℑxyzᶠᶜᶜ, ℑxyzᶜᶠᶠ, ℑxyzᶠᶜᶠ, ℑxyzᶠᶠᶜ,
    divᶜᶜᶜ, div_xyᶜᶜᵃ, div_xzᶜᵃᶜ, div_yzᵃᶜᶜ, ζ₃ᶠᶠᵃ,
    ∇², ∇²hᶜᶜᵃ, ∇²hᶠᶜᵃ, ∇²hᶜᶠᵃ, ∇⁴hᶜᶜᵃ, ∇⁴hᶠᶜᵃ, ∇⁴hᶜᶠᵃ

#####
##### Convinient aliases
#####

using Oceananigans.Grids: RegularRectilinearGrid, VerticallyStretchedRectilinearGrid, RegularLatitudeLongitudeGrid
using Oceananigans.Grids: AbstractGrid, AbstractRectilinearGrid, AbstractCurvilinearGrid, AbstractHorizontallyCurvilinearGrid

const AG  = AbstractGrid
const ARG = AbstractRectilinearGrid
const RCG = RegularRectilinearGrid
const ACG = AbstractCurvilinearGrid
const AHCG = AbstractHorizontallyCurvilinearGrid

include("areas_and_volumes.jl")
include("field_metric_products.jl")
include("difference_operators.jl")
include("derivative_operators.jl")
include("interpolation_operators.jl")
include("divergence_operators.jl")
include("vorticity_operators.jl")
include("laplacian_operators.jl")
include("interpolation_utils.jl")

end


text:
module Metal

using Reexport
using GPUArrays
using Adapt
using GPUCompiler
using LLVM
using LLVM.Interop
using Metal_LLVM_Tools_jll
using ExprTools: splitdef, combinedef

# core library
include("../lib/core/MTL.jl")
@reexport using .MTL

# device functionality
include("device/utils.jl")
include("device/pointer_abstract.jl")
include("device/pointer_ptr.jl")
include("device/pointer_buf.jl")
include("device/array.jl")
include("device/metal.jl")

# compiler
include("execution/gpucompiler.jl")
include("execution/mtl_type_conversion.jl")
include("execution/device_type_conversion.jl")
include("execution/kernel.jl")
include("execution/reflection.jl")

include("state.jl")

# array abstraction
include("utilities.jl")
include("array.jl")
include("memory.jl")
include("broadcast.jl")
include("random.jl")
include("gpuarrays.jl")

end # module


text:
using PGLib
using PowerModels
using Test

silence() #suppress PowerModels output

case = parse_file(joinpath("./pglib_opf_case3_lmbd.m"))
case_api = parse_file(joinpath("./pglib_opf_case3_lmbd.m"))
case_sad = parse_file(joinpath("./pglib_opf_case3_lmbd.m"))

@testset "pglib" begin
    @testset "open pglib case file" begin
        @test case == pglib("pglib_opf_case3_lmbd.m")
        @test case == pglib("pglib_opf_case3_lmbd")
        @test case == pglib("case3_lmbd.m")
        @test case == pglib("case3_lmbd")
        @test case == pglib("lmbd")

        @test case_api == pglib("lmbd","api")
        @test case_sad == pglib("lmbd","sad")
    end

    @testset "find pglib case" begin
        @test ["pglib_opf_case3_lmbd.m"] == find_pglib_case("lmbd")
        @test ["pglib_opf_case3_lmbd__api.m"] == find_pglib_case("lmbd","api")
        @test ["pglib_opf_case3_lmbd__sad.m"] == find_pglib_case("lmbd","sad")
        @test length(find_pglib_case()) == 59
        @test length(find_pglib_case("","sad")) == 59
        @test length(find_pglib_case("","api")) == 59
    end
end