text:
using Plots, DSGE, ModelConstructors, FileIO, Test
path = dirname(@__FILE__)

# Initialize the plotting backend
gr()
GR.inline("pdf")

# Initialize model object
m = AnSchorfheide(testing = true)
m <= Setting(:saveroot, joinpath(dirname(@__FILE__), "..", "reference"))#tempdir())
m <= Setting(:date_forecast_start, quartertodate("2015-Q4"))
m <= Setting(:date_conditional_end, quartertodate("2015-Q4"))
m <= Setting(:use_population_forecast, true)
m <= Setting(:forecast_horizons, 12)
m <= Setting(:impulse_response_horizons, 8)

# Run full-distribution forecast
estroot = normpath(joinpath(dirname(@__FILE__), "..", "reference"))
overrides = forecast_input_file_overrides(m)
overrides[:mode] = joinpath(estroot, "optimize.h5")
overrides[:full] = joinpath(estroot, "metropolis_hastings.h5")
output_vars = add_requisite_output_vars([:histobs, :forecastobs, :shockdecobs, :irfobs])

m <= Setting(:forecast_block_size, 5)
forecast_one(m, :full, :none, output_vars, verbose = :none)
compute_meansbands(m, :full, :none, output_vars; verbose = :none)
println("The following warning is expected test behavior:")

# Plot history and forecast
if haskey(ENV, "FRED_API_KEY")
    plot_history_and_forecast(m, :obs_nominalrate, :obs, :full, :none,
                              bdd_and_unbdd = true,
                              start_date = DSGE.quartertodate("2007-Q1"),
                              verbose = :none)

    plot_history_and_forecast(m, :obs_nominalrate, :obs, :full, :none,
                              bdd_and_unbdd = false,
                              start_date = DSGE.quartertodate("2007-Q1"),
                              verbose = :none)

    plot_history_and_forecast(m, :obs_nominalrate, :obs, :full, :none,
                              bdd_and_unbdd = false, bdd_and_bdd = true,
                              start_date = DSGE.quartertodate("2007-Q1"),
                              verbose = :none)

    # Plot forecast comparison
    plot_forecast_comparison(m, m, :obs_nominalrate, :obs, :full, :none,
                             bdd_and_unbdd = true,
                             start_date = DSGE.quartertodate("2007-Q1"),
                             verbose = :none)

    plot_forecast_comparison(m, m, :obs_nominalrate, :obs, :full, :none,
                             bdd_and_unbdd = false,
                             start_date = DSGE.quartertodate("2007-Q1"),
                             verbose = :none)

    # Plot shock decomposition
    plot_shock_decomposition(m, :obs_nominalrate, :obs, :full, :none, verbose = :none)

    # Plot IRF
    plot_impulse_response(m, :rm_sh, collect(keys(m.observables)), :obs, :full, :none,
                          verbose = :none)
    plot_impulse_response(m, :rm_sh, collect(keys(m.observables)), :obs, :full, :none,
                          flip = true, verbose = :none)
    plot_prior_posterior(m, :τ)
    plot_prior_posterior(m, [:τ, :σ_g])
    cloud = load("$path/../reference/output_data/an_schorfheide/ss0/estimate/raw/smc_cloud_vint=151001.jld2")["cloud"]
    plot_posterior_intervals(m, cloud = cloud)
    plot_posterior_intervals(m, [cloud, cloud])

    # Plot scenario, zeroing out measurement error
    for para in [:e_y, :e_π, :e_R]
        m[para].value = 0
    end
    alt = Scenario(:altscen, "Test Alternative Scenario", [:obs_gdp, :obs_cpi], [:g_sh, :rm_sh], "REF")
    forecast_scenario(m, alt, verbose = :none)
    scenario_means_bands(m, alt, verbose = :none)
    plot_scenario(m, :obs_nominalrate, :obs, alt, untrans = true, verbose = :none)
    plot_scenario(m, :obs_nominalrate, :obs, alt, verbose = :none)
    plot_scenario(m, :obs_nominalrate, :obs, alt, fourquarter = true, verbose = :none)
    @testset "Test scenario exception" begin
        @test_throws ErrorException plot_scenario(m, :obs_nominalrate, :obs, alt,
                                                  untrans = true, fourquarter = true)
    end

    # Hair plot
    realized = load_data(m, verbose = :none)
    hist_mb = read_mb(m, :full, :none, :histobs)
    fcast_mb = read_mb(m, :full, :none, :bddforecastobs)
    hair_plot(:obs_nominalrate, realized, [hist_mb], [fcast_mb];
              plotroot = saveroot(m), verbose = :none)

    # Things requiring estimated files
    tmp_saveroot = saveroot(m)
    m <= Setting(:saveroot, "$path/../reference")
    plot_history_and_forecast(m, :y_t, :pseudo, :full, :none,
                              bdd_and_unbdd = true,
                              start_date = DSGE.quartertodate("2007-Q1"),
                              verbose = :none, plotroot = tmp_saveroot)
    plot_forecast_comparison(m, m, :y_t, :pseudo, :full, :none,
                             bdd_and_unbdd = true,
                             start_date = DSGE.quartertodate("2007-Q1"),
                             verbose = :none, plotroot = tmp_saveroot)
    plot_history_and_forecast(m, :y_t, :pseudo, :full, :none,
                              bdd_and_unbdd = false,
                              start_date = DSGE.quartertodate("2007-Q1"),
                              verbose = :none, plotroot = tmp_saveroot)
    plot_forecast_comparison(m, m, :y_t, :pseudo, :full, :none,
                             bdd_and_unbdd = false,
                             start_date = DSGE.quartertodate("2007-Q1"),
                             verbose = :none, plotroot = tmp_saveroot)
    plot_history_and_forecast(m, :y_t, :pseudo, :full, :none,
                              bdd_and_unbdd = false, bdd_and_bdd = true,
                              start_date = DSGE.quartertodate("2007-Q1"),
                              verbose = :none, plotroot = tmp_saveroot)
    plot_shock_decomposition(m, :y_t, :pseudo, :full, :none, verbose = :none, plotroot = tmp_saveroot)
    plot_impulse_response(m, :rm_sh, collect(keys(m.pseudo_observables)), :pseudo, :full, :none,
                          verbose = :none, plotroot = tmp_saveroot)
    plot_impulse_response(m, :rm_sh, collect(keys(m.pseudo_observables)), :pseudo, :full, :none,
                          flip = true, verbose = :none, plotroot = tmp_saveroot)
    plot_altpolicies([m], :obs_gdp, :obs, :none, plotroot = tmp_saveroot)
    plot_altpolicies([m], :y_t, :pseudo, :none, plotroot = tmp_saveroot)
    plot_forecast_comparison(m, m, :obs_gdp, :obs, :full, :none, plotroot = tmp_saveroot)
    plot_forecast_comparison(m, m, :y_t, :pseudo, :full, :none, plotroot = tmp_saveroot)

    # Don't test plot_forecast_decomposition b/c on Windows, required reference data has names too long for Windows.
    @test_broken plot_forecast_decomposition(m, m, :obs_gdp, :obs, :full, :none, :none, plotroot = tmp_saveroot)
    @test_broken plot_forecast_decomposition(m, m, :y_t, :pseudo, :full, :none, :none, plotroot = tmp_saveroot)

    plot_forecast_sequence([m,m], [:mode, :mode], [:none, :none], m, :mode, :none, :obs, :obs_gdp,
                           plotroot = tmp_saveroot, start_date = DSGE.quartertodate("1960-Q1"),
                           end_date = DSGE.quartertodate("1961-Q1"))
    plot_forecast_sequence([m,m], [:full, :full], [:none, :none], m, :full, :none, :obs, :obs_gdp,
                           plotroot = tmp_saveroot, start_date = DSGE.quartertodate("1960-Q1"),
                           end_date = DSGE.quartertodate("1961-Q1"))
    plot_forecast_sequence([m,m], [:mode, :mode], [:none, :none], m, :mode, :none, :pseudo, :y_t,
                           plotroot = tmp_saveroot, start_date = DSGE.quartertodate("1960-Q1"),
                           end_date = DSGE.quartertodate("1961-Q1"))
    plot_forecast_sequence([m,m], [:full, :full], [:none, :none], m, :full, :none, :pseudo, :y_t,
                           plotroot = tmp_saveroot, start_date = DSGE.quartertodate("1960-Q1"),
                           end_date = DSGE.quartertodate("1961-Q1"))

    for fn_root in ["../reference/", "../reference/output_data/an_schorfheide/ss0/scenarios/figures/",
                    "../reference/output_data/an_schorfheide/ss0/forecast/figures/",
                    "../reference/output_data/an_schorfheide/ss0/estimate/figures/"]
        for fn in readdir(joinpath(dirname(@__FILE__), fn_root))
            if length(fn) > 3
                if fn[end - 3:end] == ".pdf"
                    rm(joinpath(dirname(@__FILE__), fn_root, fn))
                end
            end
        end
    end

    rm(joinpath(dirname(@__FILE__),
                "../reference/output_data/an_schorfheide/ss0/forecast/raw/bddforecastobs_cond=none_para=full_test.jld2"))
    rm(joinpath(dirname(@__FILE__),
                "../reference/output_data/an_schorfheide/ss0/forecast/raw/dettrendobs_cond=none_para=full_test.jld2"))
    rm(joinpath(dirname(@__FILE__),
                "../reference/output_data/an_schorfheide/ss0/forecast/raw/forecastobs_cond=none_para=full_test.jld2"))
    rm(joinpath(dirname(@__FILE__),
                "../reference/output_data/an_schorfheide/ss0/forecast/raw/shockdecobs_cond=none_para=full_test.jld2"))
    rm(joinpath(dirname(@__FILE__),
                "../reference/output_data/an_schorfheide/ss0/forecast/work/mbbddforecastobs_cond=none_para=full_test.jld2"))
    rm(joinpath(dirname(@__FILE__),
                "../reference/output_data/an_schorfheide/ss0/forecast/work/mbdettrendobs_cond=none_para=full_test.jld2"))
    rm(joinpath(dirname(@__FILE__),
                "../reference/output_data/an_schorfheide/ss0/forecast/work/mbforecastobs_cond=none_para=full_test.jld2"))
    rm(joinpath(dirname(@__FILE__),
                "../reference/output_data/an_schorfheide/ss0/forecast/work/mbshockdecobs_cond=none_para=full_test.jld2"))
end


text:
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Imports
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

using Main.JuliaSub: parseJuliaCode, parseJuliaFile
using Main.JuliaSub: lengthUnique, unionMergeWith!

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Aux values and functions
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

const QUOTE = :quote
const BLOCK = :block

const EMPTY_BLOCK = Expr(BLOCK)

# Wraps `expr` in quote expression
wrapQuote(expr) = Expr(QUOTE, expr)
# Wraps `expr` in block expression
wrapBlock(expr) = Expr(BLOCK, expr)

# Checks that after stripping line numbers,
# `testee` is equal to quoted `reference`
parseTest(testee, reference) =
    MacroTools.striplines(testee) == wrapQuote(reference)

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Tests
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

@testset "utils :: parsing Julia" begin
    @test parseTest(parseJuliaCode(""), EMPTY_BLOCK)
    @test parseTest(parseJuliaFile(testFilePath("empty.jl")), EMPTY_BLOCK)

    @test parseTest(parseJuliaCode("5"), wrapBlock(5))
    @test parseTest(parseJuliaFile(testFilePath("const10.jl")),
            wrapBlock(:(const N = 10)))
end

#--------------------------------------------------

@testset "utils :: multiset" begin
    @test length(Multiset()) == 0
    @test lengthUnique(Multiset()) == 0
    @test length(Multiset(5)) == 1
    @test lengthUnique(Multiset(5)) == 1
    @test length(Multiset(4,2,4,4,6)) == 5
    @test lengthUnique(Multiset(4,2,4,4,6)) == 3

    @test unionMergeWith!(Multiset(), Multiset(5,1,3)) ==
            Multiset(1,3,5)
    @test unionMergeWith!(Multiset("d", "b"), Multiset(), Multiset("a")) ==
            Multiset("a", "b", "d")
    @test unionMergeWith!(Multiset(1,2,2,3), Multiset(5,1)) ==
            Multiset(1,1,2,2,3,5)
    @test unionMergeWith!(Multiset(1,2,2,3), Multiset(5,1), Multiset(0)) ==
            Multiset(5,3,2,2,1,1,0)
    @test unionMergeWith!(Multiset("a","a","b"), Multiset("c","b","c","b")) ==
            Multiset("a","a", "b","b","b", "c","c")
end

ext: jl

max_stars_repo_name: prl-julia/julia-sub

max_issues_repo_name: prl-julia/julia-sub

max_forks_repo_name: prl-julia/julia-sub


text:
module Swagger

using Downloads
using URIs
using JSON
using MbedTLS
using Dates
using LibCURL

import Base: convert, show, summary, getindex, keys, length, getproperty, setproperty!, propertynames, iterate
if isdefined(Base, :hasproperty)
    import Base: hasproperty
end
import JSON: lower

abstract type SwaggerModel end
abstract type SwaggerApi end

include("client.jl")
include("json.jl")
include("val.jl")

end # module Swagger


text:
abstract type AbstractContext end

Base.@kwdef mutable struct Context <: AbstractContext
    index::Index = Index(true)
    trans_units::Vector{TranslationUnit} = TranslationUnit[]
    passes::Vector{AbstractPass} = AbstractPass[]
    dag::ExprDAG = ExprDAG()
    options::Dict = Dict()
end

Base.show(io::IO, x::Context) = print(io, "Context(...)")

function parse_header!(ctx::AbstractContext, header::AbstractString, args)
    flags = CXTranslationUnit_DetailedPreprocessingRecord
    flags |= CXTranslationUnit_SkipFunctionBodies
    push!(ctx.trans_units, parse_header(ctx.index, header, args, flags))
    return ctx
end

function parse_headers!(ctx::AbstractContext, headers::Vector{<:AbstractString}, args)
    flags = CXTranslationUnit_DetailedPreprocessingRecord
    flags |= CXTranslationUnit_SkipFunctionBodies
    tus = parse_headers(ctx.index, headers, args, flags)
    append!(ctx.trans_units, tus)
    return ctx
end

"""
    find_dependent_headers(headers::Vector{T}, args::Vector, system_dirs) where {T<:AbstractString}
Return a vector of headers to which those missing dependent headers are added.
"""
function find_dependent_headers(headers::Vector{T}, args::Vector, system_dirs) where {T<:AbstractString}
    flags = CXTranslationUnit_DetailedPreprocessingRecord
    flags |= CXTranslationUnit_SkipFunctionBodies
    all_headers = copy(headers)
    dependent_headers = T[]
    new_headers = T[]
    idx = Index(false)
    while true
        empty!(new_headers)
        for header in all_headers
            try
                tu = parse_header(idx, header, args, flags)
                GC.@preserve tu begin
                    tu_cursor = getTranslationUnitCursor(tu)
                    for cursor in children(tu_cursor)
                        is_inclusion_directive(cursor) || continue
                        file = getIncludedFile(cursor)
                        file_name = get_filename(file) |> normpath
                        (isempty(file_name) || !isfile(file_name)) && continue
                        # skip system headers
                        any(sysdir->startswith(file_name, sysdir), system_dirs) && continue
                        file_name ∈ all_headers && continue
                        file_name ∈ dependent_headers && continue
                        file_name ∈ new_headers && continue
                        push!(new_headers, file_name)
                    end
                end
            catch err
                @warn "failed to parse $header, skip..."
                continue
            end
        end
        isempty(new_headers) && break
        foreach(new_headers) do h
            @info "Found dependent header: $h"
        end
        append!(dependent_headers, new_headers)
        append!(all_headers, new_headers)
    end

    return normpath.(dependent_headers)
end

"""
    create_context(headers::Vector, args::Vector=String[], options::Dict=Dict())
Create a context from a vector of paths of headers, a vector of compiler flags and
a option dict.
"""
function create_context(headers::Vector, args::Vector=String[], options::Dict=Dict())
    ctx = Context(; options)

    push!(args, "-nostdinc")

    system_dirs = map(x->x[9:end], filter(x->startswith(x, "-isystem"), args))
    dependent_headers = find_dependent_headers(headers, args, system_dirs)

    @info "Parsing headers..."
    parse_headers!(ctx, headers, args)

    push!(ctx.passes, CollectTopLevelNode(ctx.trans_units, dependent_headers, system_dirs))
    push!(ctx.passes, LinkTypedefToAnonymousTagType())
    push!(ctx.passes, LinkTypedefToAnonymousTagType(is_system=true))
    push!(ctx.passes, IndexDefinition())
    push!(ctx.passes, CollectDependentSystemNode())
    push!(ctx.passes, IndexDefinition())
    push!(ctx.passes, CollectNestedRecord())
    push!(ctx.passes, FindOpaques())
    push!(ctx.passes, ResolveDependency())
    push!(ctx.passes, RemoveCircularReference())
    push!(ctx.passes, TopologicalSort())
    push!(ctx.passes, IndexDefinition())
    push!(ctx.passes, ResolveDependency())
    push!(ctx.passes, CatchDuplicatedAnonymousTags())
    push!(ctx.passes, CodegenPreprocessing())

    general_options = get(options, "general", Dict())
    if get(general_options, "smart_de_anonymize", true)
        push!(ctx.passes, DeAnonymize())
    end
    push!(ctx.passes, Audit())
    push!(ctx.passes, Codegen())
    push!(ctx.passes, CodegenMacro())
    # push!(ctx.passes, CodegenPostprocessing())
    if get(general_options, "add_fptr_methods", false)
        push!(ctx.passes, AddFPtrMethods())
    end
    if get(general_options, "auto_mutability", false)
        push!(ctx.passes, TweakMutability())
    end

    # support old behavior
    api_file = get(general_options, "output_api_file_path", "")
    common_file = get(general_options, "output_common_file_path", "")

    output_file_path = get(general_options, "output_file_path", "")

    if isempty(api_file) && isempty(common_file)
        if !isempty(output_file_path)
            push!(ctx.passes, ProloguePrinter(output_file_path))
            push!(ctx.passes, GeneralPrinter(output_file_path))
            push!(ctx.passes, EpiloguePrinter(output_file_path))
        else
            # print to stdout if there is no `output_file_path`
            # this is handy when playing in REPL
            push!(ctx.passes, StdPrinter())
        end
    else
        # let the user handle prologue and epilogue on their own
        push!(ctx.passes, FunctionPrinter(api_file))
        push!(ctx.passes, CommonPrinter(common_file))
    end

    return ctx
end

create_context(header::AbstractString, args=String[], ops=Dict()) = create_context([header], args, ops)

@enum BuildStage begin
    BUILDSTAGE_ALL
    BUILDSTAGE_NO_PRINTING
    BUILDSTAGE_PRINTING_ONLY
end

"""
    build!(ctx::Context, skip_printing = true)
Run all passes.
"""
function build!(ctx::Context, stage=BUILDSTAGE_ALL)
    reset_counter()
    for pass in ctx.passes
        pass isa CollectDependentSystemNode && stage != BUILDSTAGE_PRINTING_ONLY && @info "Building the DAG..."
        pass isa Codegen && stage != BUILDSTAGE_PRINTING_ONLY && @info "Emit Julia expressions..."
        if stage == BUILDSTAGE_NO_PRINTING
            pass isa AbstractPrinter && continue
        elseif stage == BUILDSTAGE_PRINTING_ONLY
            !(pass isa AbstractPrinter) && continue
        end
        pass(ctx.dag, ctx.options)
    end
    @info "Done!"
    return ctx
end

"""
    get_identifier_node(ctx::AbstractContext, s::Symbol)
Return an identifer node that has `s` as its node id.
"""
function get_identifier_node(ctx::AbstractContext, s::Symbol)
    if haskey(ctx.dag.ids, s)
        return ctx.dag.nodes[ctx.dag.ids[s]]
    else
        return nothing
    end
end

"""
    get_tagtype_node(ctx::AbstractContext, s::Symbol)
Return a tag-type node that has `s` as its node id.
"""
function get_tagtype_node(ctx::AbstractContext, s::Symbol)
    if haskey(ctx.dag.tags, s)
        return ctx.dag.nodes[ctx.dag.tags[s]]
    else
        return nothing
    end
end

function get_triple()
    is_libc_musl = occursin("musl", Base.BUILD_TRIPLET)
    if Sys.isapple() && Sys.ARCH === :aarch64
        return "aarch64-apple-darwin20"
    elseif Sys.islinux() && Sys.ARCH === :aarch64 && !is_libc_musl
        return "aarch64-linux-gnu"
    elseif Sys.islinux() && Sys.ARCH === :aarch64 && is_libc_musl
        return "aarch64-linux-musl"
    elseif Sys.islinux() && startswith(string(Sys.ARCH), "arm") && !is_libc_musl
        return "armv7l-linux-gnueabihf"
    elseif Sys.islinux() && startswith(string(Sys.ARCH), "arm") && is_libc_musl
        return "armv7l-linux-musleabihf"
    elseif Sys.islinux() && Sys.ARCH === :i686 && !is_libc_musl
        return "i686-linux-gnu"
    elseif Sys.islinux() && Sys.ARCH === :i686 && is_libc_musl
        return "i686-linux-musl"
    elseif Sys.iswindows() && Sys.ARCH === :i686
        return "i686-w64-mingw32"
    elseif Sys.islinux() && Sys.ARCH === :powerpc64le
        return "powerpc64le-linux-gnu"
    elseif Sys.isapple() && Sys.ARCH === :x86_64
        return "x86_64-apple-darwin14"
    elseif Sys.islinux() && Sys.ARCH === :x86_64 && !is_libc_musl
        return "x86_64-linux-gnu"
    elseif Sys.islinux() && Sys.ARCH === :x86_64 && is_libc_musl
        return "x86_64-linux-musl"
    elseif Sys.isbsd() && !Sys.isapple()
        return "x86_64-unknown-freebsd"
    elseif Sys.iswindows() && Sys.ARCH === :x86_64
        return "x86_64-w64-mingw32"
    end
end

function get_default_args(triple=get_triple())
    args = ["-isystem"*dir for dir in get_system_dirs(triple)]
    target = triple2target(triple)
    push!(args, "--target=$target")
    return args
end

"""
    detect_headers(include_dir, args, options::Dict=Dict())
Detect a set of headers which can span the whole directory.
"""
function detect_headers(include_dir, args, options::Dict=Dict())
    system_dirs = filter(x->startswith(x, "-isystem"), args)
    system_dirs = map(x->x[9:end], system_dirs)

    candidates = Set{String}()
    dependencies = Set{String}()
    idx = Index(false)
    flags = CXTranslationUnit_DetailedPreprocessingRecord
    flags |= CXTranslationUnit_SkipFunctionBodies
    for (root, dirs, files) in walkdir(normpath(include_dir))
        for file in files
            header = joinpath(root, file)
            try
                # FIXME: C++ could also be successful parsed, but that's not what we want
                parse_header(idx, header, args, flags)
            catch err
                @info "failed to parse $header, skip..."
                continue
            end
            push!(candidates, header)
            tu = parse_header(idx, header, args, flags)
            GC.@preserve tu begin
                tu_cursor = getTranslationUnitCursor(tu)
                header_name = spelling(tu_cursor)
                header_dir = dirname(header_name)
                @assert startswith(header_dir, include_dir)
                for cursor in children(tu_cursor)
                    is_inclusion_directive(cursor) || continue
                    incfile = getIncludedFile(cursor)
                    file_path = get_filename(incfile) |> normpath
                    (isempty(file_path) || !isfile(file_path)) && continue
                    # skip system headers
                    any(sysdir->startswith(file_path, sysdir), system_dirs) && continue
                    push!(dependencies, file_path)
                end
            end
        end
    end
    return collect(setdiff(candidates, dependencies))
end


text:
export ASHModel

MLJBase.@mlj_model mutable struct ASHModel <: PDFEstimator end

function MLJBase.fit(model::ASHModel, verbosity, X)
    n_vars = size(X, 2)
    ash_model = ash(eachcol(X)...)

    fitresult = (ash_model, n_vars)
    cache = nothing
    report = nothing
    fitresult, cache, report
end

function MLJBase.transform(model::ASHModel, fitresult, Xnew)
    ash_model, n_vars = fitresult

    ynew = broadcast(eachrow(MLJBase.matrix(Xnew))) do xr
        AverageShiftedHistograms.pdf(ash_model, xr...)
    end
    ynew = reshape(ynew, size(ynew, 1), 1)
    ynew = DataFrame(ynew, :auto)
end

function MLJBase.fitted_params(model::ASHModel, fitresult)
    ash_model, n_vars = fitresult

    (
        ash_model=ash_model,
        n_vars=n_vars,
    )
end


text:
using LibSerialPort

abstract type AbstractNexSerial end

const DEFAULT_BAUDRATE = 9600
const DEFAULT_TIMEOUT_MS = 20000

const v_uint8_eoc = [0xff, 0xff, 0xff]  # Array{UInt8,1}
# const v_char_eoc = Char.([0xff, 0xff, 0xff])  # Array{Char,1} ['ÿ', 'ÿ', 'ÿ']
const s_eoc = String(copy(v_uint8_eoc))  # "\xff\xff\xff"

"""
    _format_command(cmd)

Add \xff\xff\xff at the end of String

# Usage

```
 julia> _format_command("sendme")
"sendme\xff\xff\xff"

julia> transcode(UInt8, _format_command("sendme"))
9-element Base.CodeUnits{UInt8,String}:
 0x73
 0x65
 0x6e
 0x64
 0x6d
 0x65
 0xff
 0xff
 0xff
 ```

"""
function _format_command(cmd)
    cmd * s_eoc
end

function hasend(msg::Vector{UInt8})
    #msg[end] == 0xff && msg[end-1] == 0xff && msg[end-2] == 0xff
    msg[end-2:end] == [0xff, 0xff, 0xff]
end

function ensurehasend(msg::Vector{UInt8})
    if !hasend(msg)
        error("Message must end with $v_uint8_eoc")
    end
end

function my_bytesavailable(sp)::UInt32
    n = bytesavailable(sp)
    if n < 0
        @error "n must be positive"
        n = 0
    elseif n > typemax(UInt32)
        @error "n can't be so big"
        n = 0
    end
    n
end

"""
See https://github.com/andrewadare/LibSerialPort.jl/issues/23
https://github.com/andrewadare/LibSerialPort.jl/issues/24
"""
function my_readuntil(sp::SerialPort, delim::AbstractString, timeout_ms::Real)
    return String(readuntil(sp, Vector{Char}(delim), timeout_ms))
    #return String(take!(readuntil(sp, Vector{Char}(delim), timeout_ms)))
 end

function my_readuntil(sp::SerialPort, delim::Vector{T}, timeout_ms::Real) where {T}
    start_time = time_ns()
    out = Vector{T}()
    lastchars = T[0 for i=1:length(delim)]
    while !eof(sp)
        if (time_ns() - start_time)/1e6 > timeout_ms
            break
        end
        if my_bytesavailable(sp) > 0
            c = read(sp, T)
            push!(out, c)
            lastchars = circshift(lastchars,-1)
            lastchars[end] = c
            if lastchars == delim
                break
            end
        end
    end
    return out
end

function readuntil_eoc(sp::SerialPort, timeout_ms::Integer)
    my_readuntil(sp, v_uint8_eoc, timeout_ms)
end

"""
    wait_response(sp, timeout_ms)

"""
function wait_response(sp::SerialPort, timeout_ms::Integer)
    sp_blocking_read(sp.ref, my_bytesavailable(sp), timeout_ms)
    nb, msg = sp_blocking_read(sp.ref, my_bytesavailable(sp), timeout_ms)
    if nb >= 3
        ensurehasend(msg)
    end
    nb, msg
end

"""
    NexSerial(portname, bps=DEFAULT_BAUDRATE, args...; kwargs...)

Nextion Screen serial connexion object.

By default communication speed is set to DEFAULT_BAUDRATE = 9600 baud.
"""
struct NexSerial <: AbstractNexSerial
    _serial::SerialPort
    timeout_ms::UInt32

    function NexSerial(portname, bps=DEFAULT_BAUDRATE, timeout_ms=DEFAULT_TIMEOUT_MS, args...; kwargs...)
        # new(SerialPort(args...; kwargs...))  # with SerialPorts.jl
        sp = open(portname, bps, args...; kwargs...)
        new(sp)  # with LibSerialPort.jl
    end
end

"""
    NexSerialMock(portname, bps=DEFAULT_BAUDRATE, args...; kwargs...)

Nextion Screen serial connexion object.

This is a mock.

Only use it for unit tests that don't need hardware
"""
struct NexSerialMock <: AbstractNexSerial
    timeout_ms::UInt32

    function NexSerialMock(portname="/dev/null", bps=DEFAULT_BAUDRATE, timeout_ms=DEFAULT_TIMEOUT_MS, args...; kwargs...)
        new(timeout_ms)
    end
end

"""
    init(nexSerial)

Initialize Nextion serial communication.

Serial communication is initialized by baudrate=9600
"""
function init(nexSerial::NexSerial)
    sp = nexSerial._serial
    nb, msg = wait_response(sp, nexSerial.timeout_ms)
    @info msg
    true
end

"""
    send(nexSerial, cmd)

Send a command `cmd` to Nextion using serial communication
defined by `nexSerial` and return code.

`cmd` must be lowercased
"""
function send(nexSerial::NexSerial, cmd::String)
    write(nexSerial, cmd)
    @info "reading..."
    r = read(nexSerial)
    @info "end of read."
    r
end

"""
    write(nexSerial, cmd)

Write a command `cmd` to Nextion using serial communication
defined by `nexSerial`.
"""
function write(nexSerial::NexSerial, cmd::String)
    ser = nexSerial._serial
    @info "write '$cmd' to $ser"
    cmd = _format_command(cmd)
    write(ser, cmd)
end

"""
    read(nexSerial)

Read state of a Nextion using serial communication
defined by `nexSerial`.
"""
function read(nexSerial::NexSerial)
    ser = nexSerial._serial
    #r = read(ser, UInt8)
    r = sp_blocking_read(nexSerial._serial.ref, bytesavailable(nexSerial._serial), nexSerial.timeout_ms)
end

"""
    close(nexSerial)

Close serial communication with Nextion display.
"""
function close(nexSerial::NexSerial)
    ser = nexSerial._serial
    close(ser)
end

"""
    reset(nexSerial)

Reset Nextion device.
"""
function reset(nexSerial::NexSerial)
    send(nexSerial, "rest")
    sleep(1)
    sp = nexSerial._serial
    timeout_ms = nexSerial.timeout_ms
    #nb, r = wait_response(sp, timeout_ms)
    #r = checkcommandcomplete(sp, timeout_ms)
    #@info r
    #@assert r == [0x00, 0x00, 0x00, 0xff]
    #sleep(1)
    #r = checkcommandcomplete(sp, timeout_ms)
    r = readuntil_eoc(sp, nexSerial.timeout_ms)
    @info r
    #@assert r[end-1] == 0x88
end

"""
    checkcommandcomplete(nexSerial)

Return true if command returned succesfully
"""
function checkcommandcomplete(sp::SerialPort, timeout_ms)
    nb, r = sp_blocking_read(sp.ref, 4, timeout_ms)
    r
end

function waitcommandcomplete(sp::SerialPort, timeout_ms)
    sleep(1)
end

struct Backlight
    value::Unsigned
    persist::Bool

    function Backlight(value, persist=false)
        r = 0:100
        if !(value in r)
            error("Backlight value must be in $r")
        end
        new(value, persist)
    end
end

"""
    nexSerial.backlight = Backlight(100)  # Backlight(100, true)  # for persistance

Set backlight.
"""
function setproperty!(nexSerial::NexSerial, property::Symbol, new_val::Backlight)  # caution! new_val type
    if property == :backlight
        val = new_val.value
        persist = new_val.persist
        if persist
            send(nexSerial, "dims=$val")
        else
            send(nexSerial, "dim=$val")
        end
    else
        setfield!(obj, property, new_val)
    end
end

#function getproperty(nexSerial::NexSerial, property::Symbol)
#    if property == :backlight
#        send(nexSerial, "get dim")
#        val = receivenumber()
#        val
#    else
#        setfield!(obj, property, new_val)
#    end
#end


text:
############################################################################################
## This file is part of RLinearAlgebra.jl
##
## Overview: abstractions and methods for routines used in solving linear systems of
## equations.
##
## Contents
## - Abstract Types
## - `rsubsolve!` function documentation
## - Vector Row Projection Solvers
## - Vector Column Projection Solvers
## - Block Row Projection Solvers
## - Block Column Projection Solvers
## - Preconditioned Krylov Solvers
## - Export Statements
##
############################################################################################

# Dependencies: LinearAlgebra

#############################################
# Abstract Types
#############################################
"""
    LinSysSolveRoutine

Abstract supertype that specifies the type of linear system solver routine being deployed.
"""
abstract type LinSysSolveRoutine end

"""
    LinSysVecRowProjection <: LinSysSolveRoutine

Abstract supertype for vector row action projection methods.
"""
abstract type LinSysVecRowProjection <: LinSysSolveRoutine end

"""
    LinSysVecColProjection <: LinSysSolveRoutine

Abstract supertype for vector column action projection methods.
"""
abstract type LinSysVecColProjection <: LinSysSolveRoutine end

"""
    LinSysBlkRowProjection <: LinSysSolveRoutine

Abstract supertype for block row action projection methods.
"""
abstract type LinSysBlkRowProjection <: LinSysSolveRoutine end

"""
    LinSysBlkColProjection <: LinSysSolveRoutine

Abstract supertype for block column action projection methods.
"""
abstract type LinSysBlkColProjection <: LinSysSolveRoutine end

"""
    LinSysPreconKrylov <: LinSysSolveRoutine

Abstract supertype for preconditioned Krylov solver.
"""
abstract type LinSysPreconKrylov <: LinSysSolveRoutine end

#############################################
# `rsubsolve!` Function Documentation`
#############################################
"""
    rsubsolve!(
        type::LinSysSolveRoutine,
        x::AbstractVector,
        samp::Tuple,
        iter::Int64)

A common interface for specifying different strategies for solving a subproblem generated
    by a sampling, selecting or sketching operation on a linear system. The `type` argument
    is used to select the appropriately defined strategy. The argument `x` is the current
    iterate value for the solution. The argument `samp` depends on the subtype
    of `LinSysSolveRoutine` that is being deployed and is described below. The `iter`
    argument is the iteration counter.

The value of `samp` depends on the type of subtype of `LinSysSolveRoutine` being deployed.
To describe this, let `A` be the coefficient matrix of the system, and `b` its constant
vector.
- For `T<:LinSysVecRowProjection`, `samp` is a two-dimensional tuple where the first
    entry is a vector in the row space of `A`; and the second entry is a scalar value
    corresponding to a linear combination of the elements of `b`.
- For `T<:LinSysVecColProjection`, `samp` is a three-dimensional tuple where the first
    entry is a vector of `length(x)` corresponding to the search direction; the second
    entry is a matrix with the same number of rows as `A` (usually it is `A`);
    and the third entry is a scalar-valued residual for the normal system corresponding
    to `samp[1]'* A' * (A * x - b)`.

The function `rsubsolve!` updates the quantity `x` and any fields of `type` that must be
    updated.
"""
function rsubsolve!(
    type::Nothing,
    x::AbstractVector,
    samp::Tuple,
    iter::Int64,
)
    return nothing
end

#############################################
# Vector Row Projection Solvers
#############################################
include("linear_solver_routines/vec_row_projection_std.jl")
include("linear_solver_routines/vec_row_projection_portho.jl")
include("linear_solver_routines/vec_row_projection_fortho.jl")

#############################################
# Vector Column Projection Solvers
#############################################
include("linear_solver_routines/vec_col_projection_std.jl")
include("linear_solver_routines/vec_col_projection_portho.jl")
include("linear_solver_routines/vec_col_projection_fortho.jl")

#############################################
# Block Row Projection Solvers
#############################################
# Exact Solvers

# Approximate Solvers

#############################################
# Block Column Projection Solvers
#############################################
# Exact Solvers

# Approximate Solvers

#############################################
# Precondition Krylov Solvers
#############################################
#Blendenpik

#############################################
# Export Statements
#############################################
# See RLinearAlgebra.jl
#export LinSysSolveRoutine, LinSysVecRowProjection, LinSysVecColProjection,
#    LinSysBlkRowProjection, LInSysBlkColProjection, LinSysPreconKrylov


text:
update!(p::Printer, problem, gradient, ξ, t) = p.verbosity ≥ 2 && @info "t = $t"

function update!(writer::VTKWriter, problem, gradient, ξ, t)
    if writer.n % writer.nupdate == 0
        (; pvd, dir, filename, ifile) = writer
        (; model) = problem

        ncompartment = length(model.mesh.points)
        npoint_cmpts = size.(model.mesh.points, 2)
        inds_cmpts = [0; cumsum(npoint_cmpts[1:end])]

        vtmfile = vtk_multiblock("$dir/$(filename)_$ifile")

        for icmpt = 1:ncompartment
            points = model.mesh.points[icmpt]
            elements = model.mesh.elements[icmpt]

            cells = [
                MeshCell(VTKCellTypes.VTK_TETRA, elements[:, i]) for i = 1:size(elements, 2)
            ]
            ξ_cmpt = ξ[1+inds_cmpts[icmpt]:inds_cmpts[icmpt+1]]

            vtkfile = vtk_grid(vtmfile, points, cells)
            vtkfile["Magnetization (real part)", VTKPointData()] = real(ξ_cmpt)
            vtkfile["Magnetization (imaginary part)", VTKPointData()] = imag(ξ_cmpt)
            vtkfile["Gradient", VTKFieldData()] = gradient(t)
        end

        pvd[t] = vtmfile
        writer.ifile += 1
    end
    writer.n += 1
end

function update!(p::Plotter, problem, gradient, ξ, t)
    if p.n % p.nupdate == 0
        femesh = problem.model.mesh
        M = problem.matrices.M
        ncompartment, nboundary = size(femesh.facets)
        ξ_cmpts = split_field(femesh, ξ)
        push!(p.t[], t)
        if gradient isa ScalarGradient
            p.f[] = push!(p.f[], gradient.profile(t))
        else
            grad = Vec3f(gradient(t))
            p.g⃗[] = [grad]
            p.g⃗_hist[] = push!(p.g⃗_hist[], grad)
        end
        p.attenuation[] = push!(p.attenuation[], sum(abs, M * ξ) / p.S₀)
        p.ξ[] = ξ
        for icmpt = 1:ncompartment, iboundary = 1:nboundary
            p.magnitude[icmpt, iboundary][] = abs.(ξ_cmpts[icmpt])
            p.phase[icmpt, iboundary][] = angle.(ξ_cmpts[icmpt])
        end

        sleep(0.01)
    end
    p.n += 1
end


text:
"""
The `Pfam` module, defines functions to measure the protein contact prediction performance of information measure between column pairs from a Pfam MSA.

**Features**

- Read and download Pfam MSAs
- Obtain PDB information from alignment annotations
- Map between sequence/alignment residues/columns and PDB structures
- Measure of AUC (ROC curve) for contact prediction of MI scores

```julia
using MIToS.Pfam
```
"""
module Pfam

using MIToS.Utils
using MIToS.MSA
using MIToS.SIFTS
using MIToS.PDB
using MIToS.Information
using PairwiseListMatrices
using NamedArrays
using DataStructures
using Requires

export  # Download
        downloadpfam,
        # PDB
        Stockholm,
        getseq2pdb,
        msacolumn2pdbresidue,
        hasresidues,
        msacontacts,
        msaresidues,
        getcontactmasks

include("Download.jl")
include("PDB.jl")

function __init__()
    @require ROCAnalysis="f535d66d-59bb-5153-8d2b-ef0a426c6aff" include("rocanalysis.jl")
end

end


text:
# Methods for checking matrix sizes etc.

issquare(x::Matrix) = size(x, 1) == size(x, 2)
ispossemidef(x::Matrix) = issym(x) && (eigmin(x) >= 0)

function confirm_matrix_sizes(F, B, V, G, W)

    nx = size(F, 1)
    nu = size(B, 2)
    ny = size(G, 1)

    @assert size(F) == (nx, nx)
    @assert size(B) == (nx, nu)
    @assert size(V) == (nx, nx)

    @assert size(G) == (ny, nx)
    @assert size(W) == (ny, ny)

    return nx, ny, nu

end #confirm_matrix_sizes

# function confirm_matrix_sizes(A1, A2, A3, B1, B2, G, Q, C1, C2, C3, D1, D2, H, R, x1, P1)

#     @assert size(B2, 2) == size(D2, 2)

#     nx = size(A1, 1)
#     ny = size(C1, 1)
#     nu = size(B2, 2)

#     na1, na2  = size(A1, 2), size(A2, 2)
#     nb        = size(B1, 2)
#     nc1, nc2  = size(C1, 2), size(C2, 2)
#     nd        = size(D1, 2)

#     nq = size(Q, 1)
#     nr = size(R, 1)

#     @assert size(A1) == (nx, na1)
#     @assert size(A2) == (na1, na2)
#     @assert size(A3) == (na2, nx)

#     @assert size(B1) == (nx, nb)
#     @assert size(B2) == (nb, nu)

#     @assert size(G)  == (nx, nq)
#     @assert size(Q)  == (nq, nq)

#     @assert size(C1) == (ny, nc1)
#     @assert size(C2) == (nc1, nc2)
#     @assert size(C3) == (nc2, nx)

#     @assert size(D1) == (ny, nd)
#     @assert size(D2) == (nd, nu)

#     @assert size(H)  == (ny, nr)
#     @assert size(R)  == (nr, nr)

#     @assert length(x1)  == nx
#     @assert size(P1)    == (nx, nx)

#     return nx, ny, nu, nq, nr

# end #confirm_matrix_sizes


text:
#
module SpectralElements

using Reexport

@reexport using LinearAlgebra
@reexport using LinearSolve
@reexport using UnPack: @unpack
@reexport using Setfield: @set!
@reexport using Plots

using SparseArrays
using NNlib

include("Spectral/Spectral.jl")

import FastGaussQuadrature
import Zygote
#--------------------------------------#
# conveniences
Base.:*(op::Function, x::AbstractArray) = op(x)
linspace(zi::Number,ze::Number,n::Integer) = Array(range(zi,stop=ze,length=n))
#--------------------------------------#

# SEM building blocks
include("jac.jl")
include("interp.jl")
include("derivMat.jl")
include("semmesh.jl")
include("semq.jl")
include("ndgrid.jl")

include("ABu.jl")
include("geom.jl")
include("mask.jl")

include("time.jl")
include("mesh.jl")

include("pcg.jl")
include("adjoint.jl")

include("gatherScatter.jl")
include("grad.jl")
include("mass.jl")
include("lapl.jl")
include("hlmz.jl")
include("advect.jl")
include("diver.jl")

include("diffusion.jl")
include("convectionDiffusion.jl")
#include("stokes.jl")

include("plt.jl")
#--------------------------------------#
#function dolinsolve(integrator, linsolve; A = nothing, u = nothing, b = nothing,
#                             Pl = nothing, Pr = nothing,
#                             reltol = integrator === nothing ? nothing : integrator.opts.reltol)
#  A !== nothing && (linsolve = LinearSolve.set_A(linsolve,A))
#  b !== nothing && (linsolve = LinearSolve.set_b(linsolve,b))
#  u !== nothing && (linsolve = LinearSolve.set_u(linsolve,u))
#  (Pl !== nothing || Pr !== nothing) && (linsolve = LinearSolve.set_prec(Pl,Pr))
#
#  linres = if reltol === nothing
#    solve(linsolve;reltol)
#  else
#    solve(linsolve;reltol)
#  end
#end
#--------------------------------------#

end # module


text:
function veryearly(x::SMatrix{N,A,T}) where {N,A,T}
    x[(N - 1) >> 1 + 3,:]
end

function early(x::SMatrix{N,A,T}) where {N,A,T}
    x[(N - 1) >> 1 + 2,:]
end

function prompt(x::SMatrix{N,A,T}) where {N,A,T}
    x[(N - 1) >> 1 + 1,:]
end

function late(x::SMatrix{N,A,T}) where {N,A,T}
    x[(N - 1) >> 1,:]
end

function verylate(x::SMatrix{N,A,T}) where {N,A,T}
    x[(N - 1) >> 1 - 1,:]
end


text:
#=
ellipse_im.jl
Copyright 2019-03-05, Jeff Fessler, University of Michigan
=#

export ellipse_im, ellipse_im_params
#export ellipse_im_aspire

# using MIRT: image_geom, disk_phantom_params, downsample2, rotate2d
using Plots: plot
using Printf: @sprintf
using Test: @test, @test_throws, @inferred

"""
`phantom = ellipse_im(ig, params ;
	rot=0, oversample=1, hu_scale=1, replace=false)`

Generate ellipse phantom image from parameters:

`[x_center y_center x_radius y_radius angle_degrees amplitude]`

in
- `ig`		from `image_geom()`
- `params`	`[ne 6]` ellipse parameters

# Arguments
- `rot`			rotate ellipses by this amount [degrees]
- `oversample`	oversampling factor, for grayscale boundaries
- `hu_scale`	use 1000 to scale shepp-logan to HU
- `replace`		replace ellipse values if true, else add
- `how`			`:fast` is the only option

out
- `phantom`		[nx ny]	image (Float32)

note: `op ellipse` in aspire with `nsub=3` is `oversample=4 = 2^(3-1)` here

"""
function ellipse_im(ig::MIRT_image_geom,
		params::AbstractMatrix{<:Real} ;
		rot::Real = 0,
		oversample::Int = 1,
		hu_scale::Real = 1,
		replace::Bool = false,
		how::Symbol = :fast, # todo
	)

	if size(params,2) != 6
		throw("bad ellipse parameter vector size")
	end

	params[:,6] .*= hu_scale

#=
	if how === :fast && oversample == 1
		warning("ignoring :fast option for oversample=1")
		how = :slow
	end
=#

	if oversample > 1
		ig = ig.over(oversample)
	end
	args = (ig.nx, ig.ny, params, ig.dx, ig.dy, ig.offset_x, ig.offset_y,
		rot, oversample, replace)

	if how === :fast
		phantom = ellipse_im_fast(args...)
#	elseif how === :slow
#		phantom = ellipse_im_slow(args...)
	else
		throw("bad how $how")
	end

	if oversample > 1
		phantom = downsample2(phantom, oversample)
	end

	return phantom
end

#=
%
% ellipse_im_slow()
% brute force fine grid - can use lots of memory
%
function phantom = ellipse_im_slow(nx, ny, params, dx, dy, ...
	#offset_x, offset_y, rot, over, replace)

if size(params,2) ~= 6
	fail('bad ellipse parameter vector size')
end

% optional rotation of ellipse parameters
if rot ~= 0
	th = deg2rad(rot)
	cx = params[:,1)
	cy = params[:,2)
	params[:,1) = cx * cos(th) + cy * sin(th)
	params[:,2) = -cx * sin(th) + cy * cos(th)
	params[:,5) = params[:,5) + rot
	clear cx cy th
end

wx = (nx*over-1)/2 + offset_x * over
wy = (ny*over-1)/2 + offset_y * over
xx = ((0:nx*over-1) - wx) / over * dx
yy = ((0:ny*over-1) - wy) / over * dy
[xx yy] = ndgrid(xx, yy); % fine grid, equally spaced

phantom = zeros(nx*over, ny*over, 'single'); % fine array

ticker reset
ne = nrow(params)
for ie = 1:ne
	ticker(mfilename, ie, ne)

	ell = params[ie, :)
	cx = ell(1);	rx = ell(3)
	cy = ell(2);	ry = ell(4)
	theta = deg2rad(ell(5))
	[xr yr] = rot2(xx-cx, yy-cy, theta)
	tmp = (xr / rx).^2 + (yr / ry).^2 <= 1

	if replace
		phantom(tmp > 0) = ell(6)
	else
		phantom = phantom + ell(6) * tmp
	end
end

phantom = downsample2(phantom, over)
=#

"""
`phantom = ellipse_im_fast()`
"""
function ellipse_im_fast(nx, ny, params_in, dx, dy,
		offset_x, offset_y, rot, over, replace)

	params = copy(params_in)

	# optional rotation
	if rot != 0
		th = rot * 180/pi
		cx = params[:,1]
		cy = params[:,2]
		params[:,1] = cx * cos(th) + cy * sin(th)
		params[:,2] = -cx * sin(th) + cy * cos(th)
		params[:,5] .+= rot
	end

	phantom = zeros(Float32, nx, ny)

	wx = (nx-1)/2 + offset_x
	wy = (ny-1)/2 + offset_y
	x1 = ((0:nx-1) .- wx) * dx
	y1 = ((0:ny-1) .- wy) * dy
	(xx, yy) = ndgrid(x1, y1)

	hx = abs(dx) / 2
	hy = abs(dy) / 2

	for ie = 1:size(params,1)

		ell = Float32.(params[ie, :])
		cx = ell[1]; rx = ell[3]
		cy = ell[2]; ry = ell[4]
		theta = ell[5] * Float32(pi/180)
		value = Float32(ell[6])

		xs = xx .- cx # shift per ellipse center
		ys = yy .- cy

		# coordinates of "outer" corner of each pixel, relative to ellipse center
		xo = xs + sign.(xs) * hx
		yo = ys + sign.(ys) * hy

		# voxels that are entirely inside the ellipse:
		(xr, yr) = rotate2d(xo, yo, theta)
		is_inside = (xr / rx).^2 + (yr / ry).^2 .<= 1

		if replace
			phantom[is_inside] .= value
		else
			phantom += value * is_inside
		end

	end # ie loop

#@show typeof(phantom)
	return phantom
end # ellipse_im_fast()

"""
`phantom = ellipse_im(nx, dx, params ; kwarg...)`

square image of size `nx` by `nx`,
specifying pixel size `dx` and ellipse `params`
"""
function ellipse_im(nx::Int, dx::Real, params ; kwarg...)
	ig = image_geom(nx=nx, dx=1)
	return ellipse_im(ig, params ; kwarg...)
end

"""
`phantom = ellipse_im(nx::Int, params ; kwarg...)`

square image of size `nx` by `nx` with
pixel size `dx=1` and ellipse `params`
"""
function ellipse_im(nx::Int, params ; kwarg...)
	return ellipse_im(nx, 1., params ; kwarg...)
end

"""
`phantom = ellipse_im(nx::Int ; ny::Int=nx, dx::Real=1, kwarg...)`

image of size `nx` by `ny` (default `nx`) with specified `dx` (default 1),
defaults to `:shepplogan_emis`
"""
function ellipse_im(nx::Int ; ny::Int=nx, dx::Real=1, kwarg...)
	ig = image_geom(nx=nx, ny=ny, dx=dx)
	return ellipse_im(ig, :shepplogan_emis ; kwarg...)
end

"""
`phantom = ellipse_im(nx::Int, ny::Int ; kwarg...)`

`:shepplogan_emis` of size `nx` by `ny`
"""
function ellipse_im(nx::Int, ny::Int ; kwarg...)
	return ellipse_im(nx, ny=ny, dx=1. ; kwarg...)
end

"""
`phantom = ellipse_im(ig, code ; kwarg...)`

`code = :shepplogan | :shepplogan_emis | :shepplogan_brainweb | :southpark`
"""
function ellipse_im(ig::MIRT_image_geom, params::Symbol ; kwarg...)
	params = ellipse_im_params(ig, params)
	return ellipse_im(ig, params ; kwarg...)
end

"""
`phantom = ellipse_im(ig ; kwarg...)`

`:shepplogan` (default) for given image geometry `ig`
"""
function ellipse_im(ig::MIRT_image_geom ; kwarg...)
	return ellipse_im(ig, :shepplogan; kwarg...)
end

"""
`params = ellipse_im_params(ig::MIRT_image_geom, params::Symbol)`

`code = :shepplogan | :shepplogan_emis | :shepplogan_brainweb | :southpark`
"""
function ellipse_im_params(ig::MIRT_image_geom, params::Symbol)
	if params === :disks
		params = disk_phantom_params(fov=ig.fovs[1])
	elseif params === :shepplogan || params === :kak
		params = shepp_logan_parameters(ig.fovs..., case=:kak)
	elseif params === :shepplogan_emis || params === :emis
		params = shepp_logan_parameters(ig.fovs..., case=:emis)
	elseif params === :shepplogan_brainweb || params === :brainweb
		params = shepp_logan_parameters(ig.fovs..., case=:brainweb)
	elseif params === :southpark
		params = south_park_parameters(fov=ig.fovs[1])
	else
		throw("bad phantom symbol $params")
	end
	return params
end

"""
`params = shepp_logan_parameters(xfov, yfov ; case::Symbol)`

parameters from Kak and Slaney text, p. 255

the first four columns are unitless "fractions of field of view"
"""
function shepp_logan_parameters(xfov::Real, yfov::Real ; case::Symbol=:kak)
	params = [
	0		0		0.92	0.69	90	2
	0		-0.0184	0.874	0.6624	90	-0.98
	0.22	0		0.31	0.11	72	-0.02
	-0.22	0		0.41	0.16	108	-0.02
	0		0.35	0.25	0.21	90	0.01
	0		0.1		0.046	0.046	0	0.01
	0		-0.1	0.046	0.046	0	0.01
	-0.08	-0.605	0.046	0.023	0	0.01
	0		-0.605	0.023	0.023	0	0.01
	0.06	-0.605	0.046	0.023	90	0.01]

	params[:,[1,3]] .*= xfov/2
	params[:,[2,4]] .*= yfov/2

	if case === :emis
		params[:,6] = [1, 1, -2, 2, 3, 4, 5, 6, 1, 1]
	elseif case === :brainweb
		params[:,6] = [1, 0, 2, 3, 4, 5, 6, 7, 8, 9] # brainweb uses index 1-10
	elseif case != :kak
		throw("bad phantom case $case")
	end

	return Float32.(params)
end

"""
`param = south_park_parameters( ; fov::Real = 100)`
"""
function south_park_parameters( ; fov::Real = 100)
	xell = [
		0. 0 85 115 0 100
		0 -60 30 20 0 -80 # mouth
		30 20 25 35 30 20 # eyes
		-30 20 25 35 -30 20
		35 25 7 7 0 -100 # pupils
		-15 25 7 7 0 -100
		0 75 60 15 0 -50] # hat
	xell[:,1:4] .*= fov/256
	return Float32.(xell)
end

#%
#% ellipse_im_profile()
#%
#function ellipse_im_profile
#ig = image_geom('nx', 2^9, 'ny', 2^9+2', 'fov', 250)
#profile on
#x0 = ellipse_im(ig, [], 'oversample', 3, 'type', 'fast')
#profile off
#profile report

#=
"""
`ellipse_im()`

show docstring(s)
"""
function ellipse_im()
	@doc ellipse_im
end
=#

"""
`ellipse_im_show()`
"""
function ellipse_im_show()
	ig = image_geom(nx=2^8, ny=2^8+2, fov=250)
#	ig.offset_y = 75.6 / ig.dy

	over = 2^2
	x0 = ellipse_im(ig, oversample=over)
	p1 = jim(ig.x, ig.y, x0, title="Shepp Logan", clim=[0.9,1.1])

	x1 = ellipse_im(ig, :shepplogan_emis, oversample=over)
	p4 = jim(x1, title="Shepp Logan Emission")

	x2 = ellipse_im(ig, :southpark)
	p2 = jim(x2, title="South Park")

	x3 = ellipse_im(ig, :shepplogan_brainweb)
	p5 = jim(x3, title="Shepp Logan Brainweb")

	plot(p1,p2,p4,p5)
end

#=

# compare to aspire
function ellipse_im_aspire()
	nx = 2^6
	ig = image_geom(nx=nx, ny=nx+2, fov=2^7)
	ell = [10, 20, 30, 40, 50, 1]

	file = tempname() * "t.fld"
	file = "tt.fld"
	over = 2^2
#	com = @sprintf "echo y | op -chat 99 ellipse %s %d %d  %g %g %g %g %g %g %d" file ig.nx ig.ny ell ./ [ig.dx, ig.dx, ig.dx, ig.dx, 1, 1] log2(over)+1
	tmp = ""
	for i=1:4
		tmp = tmp * " $(ell[i] / ig.dx)"
	end
	for i=5:6
		tmp = tmp * " $(ell[i])"
	end
	op = "/Users/fessler/bin/mi64-gcc/op"
	com = "$op -chat 99 ellipse $file $(ig.nx) $(ig.ny)$tmp $(log2(over)+1)"
	com = "/bin/ls /Users/fessler/bin/mi64-gcc" # todo: why fails?
	@show com
#	cm = @cmd "$com"
	cm = `"$com"`
	@show cm

	run(cm)
return true
	asp = fld_read(file)
	p1 = jim(asp, title="aspire")

	area_asp = sum(asp[:]) * abs(ig.dx * ig.dy)

	jul = ellipse_im(ig, ell, oversample=over) # 'type', types{ii})
	area_jul = sum(jul[:]) * abs(ig.dx * ig.dy)

	p2 = jim(jul, title="jul")
	p3 = jim((jul-asp)*over^2, "difference: (jul-asp)*over^2")

	@show maximum(abs.(jul - mat)) / ell[6] * over^2
	!isapprox(jul, mat) && throw("approximation error")
#	max_percent_diff(mat, asp)

	area_real = pi * ell[3] * ell[4] * ell[6]
	@show area_real, area_asp, area_jul

	true
end
=#

# ellipse_im_test()
function ellipse_im_test()
#@test ellipse_im_aspire() # todo
	fov = 100
	@inferred shepp_logan_parameters(fov, fov)
	shepp_logan_parameters(fov, fov, case=:kak)
	shepp_logan_parameters(fov, fov, case=:emis)
	shepp_logan_parameters(fov, fov, case=:brainweb)

	# test various ways of calling
	ellipse_im(20)
#	@inferred ellipse_im(20,22) # todo-i: why fails?
	ellipse_im(20,22)
	ellipse_im(30, :shepplogan_emis, oversample=2)

	ig = image_geom(nx=80, dx=1)
	ellipse_im_params(ig, :kak)

	params = shepp_logan_parameters(ig.fovs..., case=:brainweb)
	@inferred ellipse_im_fast(ig.nx, ig.ny, params, ig.dx, ig.dy,
		ig.offset_x, ig.offset_y, 0, 2, false)
	ellipse_im(ig, :disks)
	ellipse_im(ig, params, oversample=2)
	ellipse_im(ig, params, how=:fast, replace=true, rot=30)

	@test_throws String ellipse_im(ig, :bad)
	@test_throws String ellipse_im(ig, params, how=:bad)
	@test_throws String ellipse_im(ig, params')
	@test_throws String shepp_logan_parameters(ig.fovs..., case=:bad)

	ellipse_im(ig, params)

#=
	x1 = ellipse_im(100, :shepplogan_emis)
	x2 = ellipse_im(100, :shepplogan_emis, oversample=2)
	plot(jim(x1), jim(x2), jim(x2-x1))
=#

	ellipse_im_show()

#=
	if has_aspire
		ellipse_im_aspire() # todo
	end
=#
	true
end

"""
`ellipse_im(:test)`
self test
"""
function ellipse_im(test::Symbol)
	test != :test && throw(ArgumentError("test $test"))
	@test ellipse_im_test()
	true
end


text:
# for GISSMO L-Phenylalanine.

# if save_BSON_flag
#     save_path = joinpath(project_folder, "alignment_results.bson")
#     BSON.bson(save_path, region_min_dist = region_min_dist,
#     minfs = minfs,
#     minxs = minxs,
#     rets = rets)
# end

function plotregion(P, U, q_fixed_U, q_free_U, P_y, y, P_cost, y_cost, display_threshold_factor, display_reduction_factor,
    save_folder, title_string, file_name;
    save_plot_flag = true,
    display_plot_flag = true,
    canvas_size = (1000, 400))

    # # reduce the plotting positions for low signal regions. Otherwise the plot store size will be too large, and the time to load the plot will be long.
    # inds, _ = NMRSpectraSimulator.prunelowsignalentries(q_U, display_threshold_factor, display_reduction_factor)
    # P_display = P[inds]
    # U_display = U[inds]
    # q_U_display = q_U[inds]

    P_display = P
    U_display = U
    q_free_U_display = q_free_U
    q_fixed_U_display = q_fixed_U

    # plot.
    isdir(save_folder) || mkpath(save_folder)

    plots_save_path = joinpath(save_folder, file_name)
    #title_string = "$(project_name) alignment results, region $(r), real part"

    plot_obj = Plots.plot( P_display,
        real.(q_free_U_display),
        title = title_string,
        label = "free",
        seriestype = :line,
        ticks = :native,
        xlims = (P_display[1],P_display[end]),
        hover = P_display,
        linewidth = 4,
        xlabel = "ppm",
        ylabel = "real part of spectrum",
        xflip = true,
        size = canvas_size)

    Plots.plot!(plot_obj, P_display, real.(q_fixed_U_display), label = "fixed",
        seriestype = :line,
        linestyle = :dot,
        xflip = true,
        linewidth = 4)

    Plots.plot!(plot_obj, P_cost, real.(y_cost), label = "fit data",
        markershape = :circle,
        seriestype = :scatter,
        xflip = true)

    if save_plot_flag
        Plots.savefig(plot_obj, plots_save_path)
    end

    if display_plot_flag
        display(plot_obj)
    end

    return nothing
end

N_viz = 50000
canvas_size = (1000, 400)

U = LinRange(u_min, u_max, N_viz)
P = hz2ppmfunc.(U)
U_rad = U .* (2*π)

#### visualize.
# minxs[1][1] = 0.000
# minxs[2][1] = 0.000
display_reduction_factor = 100
display_threshold_factor = 0.05/10
if "L-Isoleucine" in molecule_names
    display_reduction_factor = 1
    display_threshold_factor = 0.001/10
end
save_plot_flag = true
display_flag = true

#println("objective: $(minf), return status: $(rets)")
println("shift variable:")
display(minx)
println()

q2 = uu->NMRSpectraSimulator.evalitpproxymixture(uu, As, Es; w = w)

obj_func(minx)
q_free_U = q2.(U_rad)

# file_name = "fitted_kappa_alpha.html"
# title_string = "L-Phenylalanine, fitted κ_α, real part"

fill_val = Es[1].κs_α[1][1]
fill!(Es[1].κs_α[1], fill_val)
fill!(Es[1].κs_α[2], fill_val)
q_fixed_U = q2.(U_rad)

file_name = "kappa_alpha.html"
title_string = "$(molecule_names[1]), fixed vs. free κ_α, real part"

save_folder = project_folder

plotregion(P, U, q_fixed_U, q_free_U, P_y, y, P_cost, y_cost,
    display_threshold_factor, display_reduction_factor,
    save_folder, title_string, file_name;
    save_plot_flag = save_plot_flag,
    display_plot_flag = display_flag,
    canvas_size = canvas_size)


text:
module ThermodynamicsBase

include("abstracts.jl")
include("constants.jl")
include("statestructs.jl")
include("eos_api.jl")

export AbstractPhaseState, AbstractSubstance, AbstractThermoBackend,
AbstractMixture, AbstractEquationOfState, AbstractMixingRules, AbstractProperty

export R, thermocalc, thermocalc!

end


text:
using Test

@testset "UWDDynam" begin
  include("uwd_dynam.jl")
end

@testset "DWDDynam" begin
  include("dwd_dynam.jl")
end

@testset "CPGDynam" begin
  include("cpg_dynam.jl")
end

@testset "Trajectories" begin
  include("trajectories.jl")
end


text:
export ZabbixTrigger, get_triggers

struct ZabbixTrigger
  id::Int
  description::String
  expression::String
end
ZabbixTrigger(dict) =
  ZabbixTrigger(
    parse(Int, dict["triggerid"]),
    dict["description"],
    dict["expression"],
  )

showres(t::ZabbixTrigger; offset=0) =
  println(repeat(" ", offset), "ZabbixTrigger: ", t.id, " => ", t.description)

function get_triggers(z, host::ZabbixHost)
  triggers = ZabbixTrigger[]
  _triggers = @checkerr make_request(z, "trigger.get", Dict("hostids"=>host.id))
  for t in _triggers
    push!(triggers, ZabbixTrigger(t))
  end
  triggers
end


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

## Based on "Multi-Threading and One-Sided Communication in Parallel LU Factorization"
## http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.138.4361&rank=7

hpl_seq(A, b) = hpl_seq(A, b, 32)

function hpl_seq(A::Matrix, b::Vector, blocksize::Integer)
    n = size(A,1)
    A = [A b]

    B_rows = collect(range(0, n, step=blocksize))
    B_rows[end] = n
    B_cols = [B_rows; [n+1]]
    nB = length(B_rows)
    depend = zeros(Bool, nB, nB) # In parallel, depend needs to be able to hold futures

    ## Small matrix case
    if nB <= 1
        x = A[1:n, 1:n] \ A[:,n+1]
        return x
    end

    ## Add a ghost row of dependencies to boostrap the computation
    for j=1:nB; depend[1,j] = true; end

    for i=1:(nB-1)
        ## Threads for panel factorizations
        I = (B_rows[i]+1):B_rows[i+1]
        #(depend[i+1,i], panel_p) = spawn(panel_factor_seq, I, depend[i,i])
        (depend[i+1,i], panel_p) = panel_factor_seq(A, I, depend[i,i])

        ## Threads for trailing updates
        for j=(i+1):nB
            J = (B_cols[j]+1):B_cols[j+1]
            #depend[i+1,j] = spawn(trailing_update_seq, I, J, panel_p, depend[i+1,i],depend[i,j])
            depend[i+1,j] = trailing_update_seq(A, I, J, panel_p, depend[i+1,i],depend[i,j])
        end
    end

    ## Completion of the last diagonal block signals termination
    #wait(depend[nB, nB])

    ## Solve the triangular system
    x = triu(A[1:n,1:n]) \ A[:,n+1]

    return x

end ## hpl()

### Panel factorization ###

function panel_factor_seq(A, I, col_dep)
    n = size(A, 1)

    ## Enforce dependencies
    #wait(col_dep)

    ## Factorize a panel
    K = I[1]:n
    panel_p = lu!(view(A, K, I)).p

    ## Panel permutation
    panel_p = K[panel_p]

    return (true, panel_p)

end ## panel_factor_seq()

### Trailing update ###

function trailing_update_seq(A, I, J, panel_p, row_dep, col_dep)
    n = size(A, 1)

    ## Enforce dependencies
    #wait(row_dep, col_dep)

    ## Apply permutation from pivoting
    K = (I[end]+1):n
    A[I[1]:n, J] = A[panel_p, J]

    ## Compute blocks of U
    L = tril(A[I,I],-1) + LinearAlgebra.I
    A[I, J] = L \ A[I, J]

    ## Trailing submatrix update
    if !isempty(K)
        A[K,J] = A[K,J] - A[K,I]*A[I,J]
    end

    return true

end ## trailing_update_seq()


text:
"""
Computes psuedo ranges

$SIGNATURES
<<<<<<< HEAD
´dc´: Decoder containing ephemeris data of satellite
´code_phase´: Code phase of signal at time of measure
´carrier_phase´: Center frequency of carrier signal
=======
´sat_state´: satellite state, combining decoded data, code- and carrierphase
>>>>>>> dev

Computes relative pseudo ranges of satellite vehicles.
The algorithm is based on the common reception method.
"""
function pseudo_ranges(
    sat_states::AbstractVector{SatelliteState{Float64}}
    )

    for i in 1:length(sat_states)
        is_sat_healthy_and_decodable(sat_states[i].decoder_state) || throw(BadData("SV not decoded properly, #decoder: " * string(i)))
    end

    N_sats = length(sat_states)
    c = sat_states[1].decoder_state.constants.c

    times = map(i -> calc_corrected_time(sat_states[i]), 1:N_sats)

    t_ref = maximum(times)

    τ = map(i -> t_ref - times[i], 1:N_sats)
    pseudoranges = τ .* c
    return pseudoranges
end


text:
using DSB2018
using Flux
using ImageView
using Images
using Plots

function form_regressor(this::DSB2018.Data.Image)
    x̂ = Float64.(Gray.(this.image))

    dens = DSB2018.Model.to_density(x̂, 4)
    dens = abs.(dens) / sum(abs.(dens))

    adaptive_threshold = DSB2018.Model.skimage_filters.threshold_local(dens, 11, offset=0.0)
    global_threshold   = DSB2018.Model.skimage_filters.threshold_otsu(dens)

    w_adaptve = 0.1
    threshold = adaptive_threshold * w_adaptve .+ global_threshold * (1.0 - w_adaptve)

    img_t = dens .> threshold
    img_t = DSB2018.Model.skimage_morphology.remove_small_objects(img_t, min_size=4)
    img_t = DSB2018.Model.ndi.binary_fill_holes(img_t)
    img_f = Float64.(img_t)

    grad_y, grad_x, edges, orient = imedge(img_f, KernelFactors.ando3, "replicate")

    dens, img_f, edges
end

function form_regressand(this::DSB2018.Data.Image)
    ŷ = Float64.(this.mask .> 0)

    grad_y, grad_x, edges_hat, orient = imedge(ŷ, KernelFactors.ando3, "replicate")
    ŵ = copy(edges_hat).^0.5 + 1

    ŷ, ŵ
end

function form_prediction(model::Flux.Chain,
                            batch_size::Tuple{Int64, Int64},
                            img::DSB2018.Data.Image;
                            numsamples::Int64 = 100)

    dens, img_f, edges = form_regressor(img)
    batch_indices = DSB2018.Model.make_batches(size(dens), batch_size, 1)

    batch_indices = rand(batch_indices, numsamples)

    dens_batches  = DSB2018.Model.make_batches(dens, batch_indices)
    img_batches   = DSB2018.Model.make_batches(img_f, batch_indices)
    edges_batches = DSB2018.Model.make_batches(edges, batch_indices)

    res = zeros(size(img.image))
    num_visits = zeros(size(img.image))

    x = Array{Float64}((batch_size[1], batch_size[2], 3, 1))
    for k in eachindex(batch_indices)
        x[:, :, 1, 1] .= dens_batches[k]
        x[:, :, 2, 1] .= img_batches[k]
        x[:, :, 3, 1] .= edges_batches[k]
        prediction = model(x)

        idx = batch_indices[k]
        res[idx[1], idx[2]] .+= prediction.data
        num_visits[idx[1], idx[2]] .+= 1.0
    end

    res = res ./ num_visits
    res[num_visits .== 0.0] = 0.0

    N = sum(abs.(res))
    if N == 0.0
        return res
    end

    return res / N
end

function form_data(this::DSB2018.Data.Image,
                    batch_size::Tuple{Int64, Int64},
                    num::Int64)

    dens, img_f, edges = form_regressor(this)
    y, w = form_regressand(this)

    batch_indices = DSB2018.Model.make_batches(size(this.image), batch_size, 10)

    rnd_batch_indices = unique(rand(batch_indices, num))

    dens_batches   = DSB2018.Model.make_batches(dens, rnd_batch_indices)
    img_f_batches  = DSB2018.Model.make_batches(img_f, rnd_batch_indices)
    edges_batches  = DSB2018.Model.make_batches(edges, rnd_batch_indices)
    reg_batches    = DSB2018.Model.make_batches(y, rnd_batch_indices)
    weight_batches = DSB2018.Model.make_batches(w, rnd_batch_indices)

    num_batches = length(rnd_batch_indices)

    target = Array{Float64}((batch_size[1], batch_size[2], 3, num_batches))
    for k = 1 : num_batches
        target[:, :, 1, k] = dens_batches[k]
        target[:, :, 2, k] = img_f_batches[k]
        target[:, :, 3, k] = edges_batches[k]
    end

    train_data = Vector{Tuple{Array{Float64, 4}, Matrix{Float64}, Matrix{Float64}}}(0)
    for k = 1 : num_batches
        t = Array{Float64, 4}((size(target)[1], size(target)[2], size(target)[3], 1))
        t[:, :, :, 1] = target[:, :, :, k]

        push!(train_data, (t, reg_batches[k], weight_batches[k]))
    end

    train_data
end

function get_metric(this::DSB2018.Data.Image, res::Matrix{Float64})

    adaptive_threshold = DSB2018.Model.skimage_filters.threshold_local(res, 11, offset=0.0)
    global_threshold   = DSB2018.Model.skimage_filters.threshold_otsu(res)

    w_adaptve = 0.1
    threshold = adaptive_threshold * w_adaptve .+ global_threshold * (1.0 - w_adaptve)

    img_t = res .> threshold
    img_t = DSB2018.Model.skimage_morphology.remove_small_objects(img_t, min_size=10)
    img_t = DSB2018.Model.ndi.binary_fill_holes(img_t)

    labeled_prediction = Images.label_components(img_t)
    prediction_vector  = DSB2018.Model.split_labeled_matrix(labeled_prediction)
    truth_images = getTrainImageMasks(this)

    return metric(prediction_vector, truth_images)
end

num_images_to_process = 30
batch_size = (128, 128)
num_samples_per_image = 10

ids = DSB2018.Data.getTrainImageIds()
XT = Tuple{Array{Float64,4},Array{Float64,2},Array{Float64,2}}
train_data = Vector{XT}(0)
for k = 1 : num_images_to_process
    img = DSB2018.Data.loadTrainImage(ids[k])
    append!(train_data, form_data(img, batch_size, num_samples_per_image))
end

all_images = DSB2018.Data.loadTrainImage.(ids)

function loss(this::Tuple) #train_data
    return loss(this[1], this[2], this[3])
end

function loss(x, y, w)
    r = nn_model(x)
    δ = (r .- y) .* w
    z = sum(δ.^2)/length(y)
    (σ.(z/10) - 0.5) / 0.005*0.02
end

bx = batch_size[1]
by = batch_size[2]

struct IMResize
    size::Tuple{Int64, Int64}
end

Flux.treelike(IMResize)
function (this::IMResize)(x)
    return Flux.Tracker.track(Images.imresize, x, this.size)
end

import Flux.Tracker.back

function back(::typeof(Images.imresize), Δ, xs::Flux.TrackedArray, _...)
    back(xs, Images.imresize(Δ, size(xs)))
end

rz_down1 = IMResize((88, 88))
rz_down2 = IMResize((56, 56))

rz_up1   = IMResize((88, 88))
rz_up2   = IMResize((bx, by))

nn_model = Chain(Flux.Conv((9, 9), 3 => 1, pad = (4, 4), σ),
                x -> reshape(x, (bx, by)),
                x -> rz_down1(x),
                x -> reshape(x, (88, 88, 1, 1)),
                Flux.Conv((9, 9), 1 => 1, pad = (4, 4), σ),
                x -> reshape(x, (88, 88)),
                x -> rz_down2(x),
                x -> reshape(x, (56, 56, 1, 1)),
                Flux.Conv((9, 9), 1 => 1, pad = (4, 4), σ),
                x -> reshape(x, (56, 56)),
                x -> rz_up1(x),
                x -> reshape(x, (88, 88, 1, 1)),
                Flux.Conv((9, 9), 1 => 1, pad = (4, 4), σ),
                x -> reshape(x, (88, 88)),
                x -> rz_up2(x)
                )

opt = Flux.Nesterov(params(nn_model), 0.2)
train_percent = 0.6

ltest  = Float64[]
ltrain = Float64[]

num_train = Int64(round(length(train_data)*train_percent))
num_test  = length(train_data) - num_train

jtrain = unique(rand(1 : length(train_data), num_train))
jtest  = setdiff(1:length(train_data), jtrain)

check_points = Vector{Flux.Chain}(0)

Juno.@progress "Training" for k = 1 : 3000
    Flux.train!(loss, train_data[jtrain], opt)
    @show k

    train_loss = mean(loss.(train_data[jtrain]))
    test_loss = mean(loss.(train_data[jtest]))

    push!(ltrain, train_loss.tracker.data)
    push!(ltest, test_loss.tracker.data)
    push!(check_points, deepcopy(nn_model))

    if (mod(k, 10) == 0) || (k == 1)
        plt_train = plot(ltrain[1:end],
                        title = @sprintf("Train Loss: %2.6f", train_loss),
                        fillrange = [minimum(ltrain) maximum(ltrain)],
                        fillalpha = 0.3)

        plt_valid = plot(ltest[1:end],
                        title = @sprintf("Test Loss: %2.6f", test_loss),
                        fillrange = [minimum(ltest) maximum(ltest)],
                        fillalpha = 0.3)

        display(plot(plt_train, plt_valid, layout = grid(2, 1)))
    end
end

scores = Vector{Float64}(length(all_images))
Juno.@progress "Metric" for k in eachindex(all_images)
    scores[k] = get_metric(all_images[k],
                        form_prediction(nn_model, batch_size, all_images[k], numsamples = 90))
    display(scatter(scores[1:k], title=@sprintf("%2.4f", mean(scores[1:k]))))
end

scores_prev = Vector{Float64}(length(all_images))
Juno.@progress "Metric2" for k in eachindex(scores_prev)
    x̂ = Float64.(Gray.(all_images[k].image))
    p̂ = DSB2018.Model.to_density(x̂, 4)
    scores_prev[k] = get_metric(all_images[k], p̂)
end


text:
include("boardrobots.jl")

"""
    read_board(boardstr::String; autodetectformat::Bool=false, emptycell::Char=' ', usecolors::Bool=false)

Return a board represented as a `Matrix` from a `String`. Two formats are
handled by default by setting `autodetectformat=true`.
"""
function read_board(boardstr::String; autodetectformat::Bool=false, emptycell::Char=' ', usecolors::Bool=false)
    board = initialize_board()
    if autodetectformat
        if '-' in boardstr
            emptycell = '-'
        else
            emptycell = ' '
        end
        usecolors = ('r' in boardstr)
    end
    boardstr = replace(boardstr, (emptycell => '0'))
    if usecolors
        colors2int = Dict{Char, Char}([
        'r' => '1',
        'o' => '2',
        'g' => '3',
        'p' => '4',
        'y' => '5',
        'b' => '6'
        ])
        for (key, value) in colors2int
            boardstr = replace(boardstr, (key => value))
        end
    end
    boardstr_lines = split(boardstr, "\n")
    for i in 1:min(5,length(boardstr_lines))
        for j in 1:min(5,length(boardstr_lines[i]))
            board[i,j] = parse(Int, boardstr_lines[i][j])
        end
    end
    return board
end

"""
    read_challenge(file::String; autodetectformat::Bool=false, emptycell::Char=' ', usecolors::Bool=false)

Return a board represented as a `Matrix` from a text file. Two formats are
handled by default by setting `autodetectformat=true`.
"""
function read_challenge(file::String; autodetectformat::Bool=false, emptycell::Char=' ', usecolors::Bool = false)
    return read_board(read(file, String), autodetectformat=autodetectformat, emptycell=emptycell, usecolors=usecolors)
end

max_stars_repo_name: remi-garcia/LunarLockout

max_issues_repo_name: remi-garcia/LunarLockout

max_forks_repo_name: remi-garcia/LunarLockout


text:
module LinearRegression
  using MLLibrary
  function linearRegression(X,Y, qrFunc)
    qrFunc(X,Y)
    return MLLibrary.LinearAlgebra.Commons.backsolve(X,Y)

  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
module Grids

export ARCGrid,ARCIO,ARCTask, ARCDiff, ARCDiffGrid, edges_of_grid, colors
import OffsetArrays
using Images

const colors = [
    colorant"#000000",
    colorant"#0074D9",
    colorant"#FF4136",
    colorant"#2ECC40",
    colorant"#FFDC00",
    colorant"#AAAAAA",
    colorant"#F012BE",
    colorant"#FF851B",
    colorant"#7FDBFF",
    colorant"#870C25",
]

const ARCGrid = Matrix{Int8}
const ARCIO = NamedTuple{(:i,:o),Tuple{ARCGrid,ARCGrid}}

struct ARCTask
    ios::Vector{ARCIO}
    test_ios::Vector{ARCIO}
    path::String
end

struct ARCDiff
    A :: ARCGrid # normal indices
    B :: OffsetArrays.OffsetMatrix # index like A! has whatever size B has
    # match :: OffsetArrays.OffsetMatrix{Bool,BitMatrix}
    match :: BitMatrix # index like A (same size as A)
    # edges_match :: OffsetArrays.OffsetArray{Bool,3, BitArray{3}}
    edges_match :: BitArray{3} # index like A (same size as A)
    offsets :: Tuple{Int,Int} # offset of B relative to A (literally ArcDiff.B.offsets)
    padded_offsets :: CartesianIndex{2} # if you collect(paddedviews(fill,A,B)) and need to erase the, adding these offsets to A indices will give you indices in the larger image (basically, if B has negative offsets (up and to the left of A) then we need to shift A indices by the magnitude of that offset in the positive direction)
end
function ARCDiff(A,B,match,edges_match)
    ARCDiff(
        A,
        B,
        match,
        edges_match,
        B.offsets,
        CartesianIndex(.-min.(B.offsets,0))
    )
end

struct ARCDiffGrid
    grid :: Matrix{ARCDiff}
    A :: ARCGrid
    B :: ARCGrid
end

function edges_of_grid(grid::ARCGrid)
    edges = falses(size(grid)...,4) :: BitArray
    # Is = CartesianIndices(grid)
    # Ifirst, Ilast = first(Is), last(Is)
    # for I in Is:
    #     edges[I,0] =

    # up and down calculated both at once, since "down" for each cell is same as "up" for the one below it
    edges[2:end,:,1] .= edges[1:end-1,:,2] .= grid[2:end,:] .!= grid[1:end-1,:]
    # left and right - all the same stuff but with d1 and d2 swapped by symmetry
    edges[:,2:end,3] .= edges[:,1:end-1,4] .= grid[:,2:end] .!= grid[:,1:end-1]

    return edges

#     edges[:,:,1] .= grid != OffsetArray(grid,-1,0) # up
#     edges[:,:,2] .= grid != OffsetArray(grid,1,0) # down
#     edges[:,:,3] .= grid != OffsetArray(grid,0,-1) # left
#     edges[:,:,4] .= grid != OffsetArray(grid,0,1) # right
end

# struct ARCCell
#     color :: Int8
#     up :: Bool
#     down :: Bool
#     left :: Bool
#     right :: Bool
# end
end


text:
module NavNodeTests

using Base.Test
using Compat
import Documenter: Documents, Builder
import Documenter.Documents: NavNode

type FakeDocumentInternal
    pages   :: Dict{Compat.String, Void}
    navlist :: Vector{NavNode}
    FakeDocumentInternal() = new(Dict(), [])
end
type FakeDocument
    internal :: FakeDocumentInternal
    FakeDocument() = new(FakeDocumentInternal())
end
@test fieldtype(FakeDocumentInternal, :navlist) == fieldtype(Documents.Internal, :navlist)

pages = [
    "page1.md",
    "Page2" => "page2.md",
    "Section" => [
        "page3.md",
        "Page4" => "page4.md",
        "Subsection" => [
            "page5.md",
        ],
    ],
    "page6.md",
]
doc = FakeDocument()
doc.internal.pages = Dict(map(i -> "page$i.md" => nothing, 1:8))
navtree = Builder.walk_navpages(pages, nothing, doc)
navlist = doc.internal.navlist

@test length(navlist) == 6
for (i,navnode) in enumerate(navlist)
    @test get(navnode.page) == "page$i.md"
end

@test isa(navtree, Vector{NavNode})
@test length(navtree) == 4
@test navtree[1] === navlist[1]
@test navtree[2] === navlist[2]
@test navtree[4] === navlist[6]

section = navtree[3]
@test get(section.title_override) == "Section"
@test isnull(section.page)
@test length(section.children) == 3

navpath = Documents.navpath(navlist[5])
@test length(navpath) == 3
@test navpath[1] === section
@test navpath[3] === navlist[5]

end # module NavNodeTests


text:
print("convert, ")

n = NamedArray(randn(5, 10))

@test convert(Array, n) == n.array
@test convert(NamedArray, n) == n

for T in [Float32, Float16]
 	for matrix in (map(T, n), convert(NamedArray{T}, n))
 		@test matrix.array == map(T, n.array)
 		@test names(matrix) == names(n)
 		@test dimnames(matrix) == dimnames(n)
 	end
end

# issue 43

x = @inferred NamedArray(Array([1 2; 3 4]), (["a","b"], [10,11]), (:rows,:cols))
@test convert(Array, x) == x.array

# Test convert don't change array type

n_sparray = NamedArray(SparseArrays.sprand(Int, 5, 5, 0.5))
@test isa(n_sparray.array, SparseMatrixCSC{Int,Int})
@test isa(convert(NamedArray{Float64}, n_sparray).array, SparseMatrixCSC{Float64,Int})


text:
using SomaticEvolution
using Test
using Random
using StatsBase

tests = ["initialisation","multilevel","simulations","testio","process_mutations"]

@testset "SomaticEvolution.jl" begin
    for test in tests
        @testset "$test" begin
            include(test*".jl")
        end
    end
end


text:
function epsilon_dominates_clear(u, v, epsilon)
  veps = v .+ epsilon
  all(u .<= veps) && any(u .< veps)
end

function epsilon_dominates_fast(u, v, epsilon)
  pareto_dominates_fast(u, v .+ epsilon)
end

function pareto_dominates_clear(u, v)
  all(u .<= v) && any(u .< v)
end

function pareto_dominates_fast(u, v)
  one_smaller = false
  for i in 1:length(u)
    if u[i] > v[i]
      return false
    elseif u[i] < v[i]
      one_smaller = true
    end
  end
  return one_smaller
end

# Returns two booleans where the first indicates e-box dominance and the 2nd indicates
# epsilon progress
function epsilon_box_dominates_and_epsilon_progress(u, v, epsilon)
  uindex = floor(u ./ epsilon)
  vindex = floor(v ./ epsilon)
  inner_epsilon_box_dominates_and_epsilon_progress(
    u, uindex, (epsilon .* uindex),
    v, vindex, (epsilon .* vindex))
end

epsilon_box_dominates(u, v, epsilon) = epsilon_box_dominates_and_epsilon_progress(u, v, epsilon)[1]
epsilon_box_progress(u, v, epsilon) = epsilon_box_dominates_and_epsilon_progress(u, v, epsilon)[2]

# Inner boolean expression when the intermediate results have already been calculated.
# This is to speed up processing when there is a whole archive of solutions to compare
# a new solution to.
function inner_epsilon_box_dominates_and_epsilon_progress(
  u, uindex, uindextimesepsilon,
  v, vindex, vindextimesepsilon)

  if pareto_dominates_fast(uindex, vindex)
    return (true, true)
  elseif all(uindex .== vindex)
    smaller_distance = norm(u .- uindextimesepsilon) < norm(v .- vindextimesepsilon)
    return (smaller_distance, false)
  else
    return (false, false)
  end

end


text:
"""
ITensor is a library for rapidly creating correct and efficient tensor network algorithms.

An ITensor is a tensor whose interface is independent of its memory layout. ITensor indices are 'intelligent' meaning they carry extra information and 'recognize' each other automatically when contracting or adding ITensors.

The ITensor library includes composable and extensible algorithms for optimizing and transforming tensor networks, such as matrix product state and matrix product operators.

# Example Usage

Define tensor indices i and j

    i = Index(2, "i")
    j = Index(3, "j")

Make an ITensor with these indices

    A = ITensor(i,j)

Set the i==2,j==1 element to -2.6

    A[j=>1,i=>2] = -2.6
    A[i=>2,j=>1] = -2.6 #this has the same effect

Make an ITensor with random elements

    B = randomITensor(j,i)

Add ITensors A and B together (ok that indices in different order)

    C = A + B

# Other Features of ITensor

- Tools for **tensor networks**, such as matrix product states (MPS) / tensor trains (TT)
- **Algorithms** for solving linear equations in MPS form (such as DMRG) or for integrating differential equations ("time evolving MPS")
- ITensors can have **sparse data** internally, such as block sparsity or diagonal sparsity, while having the same interface as dense ITensors
- ITensors can have **symmetry properties** (invariance or equivariance) under group transformations of the indices. In physics terminology such ITensors conserve quantum numbers.

# Documentation and Resources

ITensor website: https://itensor.org/

Documentation: https://itensor.github.io/ITensors.jl/stable/

"""
module ITensors

#####################################
# External packages
#
using BitIntegers
using ChainRulesCore
using Compat
using HDF5
using IsApprox
using KrylovKit
using LinearAlgebra
using NDTensors
using PackageCompiler
using Pkg
using Printf
using Random
using SerializedElementArrays
using StaticArrays
using TimerOutputs
using Zeros

#####################################
# General utility functions
#
include("utils.jl")

#####################################
# ContractionSequenceOptimization
#
include("ContractionSequenceOptimization/ContractionSequenceOptimization.jl")
using .ContractionSequenceOptimization

#####################################
# LazyApply
#
include("LazyApply/LazyApply.jl")
using .LazyApply: Applied, Sum, ∑, Prod, ∏, Scaled, α, coefficient, materialize

#####################################
# Ops
#
include("Ops/Ops.jl")
using .Ops

#####################################
# Exports
#
include("exports.jl")

#####################################
# Imports
#
include("imports.jl")

#####################################
# Global Variables
#
include("global_variables.jl")

#####################################
# Algorithm type for selecting
# different algorithm backends
# (for internal or advanced usage)
#
include("algorithm.jl")

#####################################
# Index and IndexSet
#
include("lastval.jl")
include("smallstring.jl") # Not currently using in TagSet
include("readwrite.jl")
include("not.jl")
include("tagset.jl")
include("arrow.jl")
include("symmetrystyle.jl")
include("index.jl")
include("set_operations.jl")
include("indexset.jl")

#####################################
# ITensor
#
include("itensor.jl")
include("broadcast.jl")
include("decomp.jl")
include("iterativesolvers.jl")

#####################################
# QNs
#
include("qn/qn.jl")
include("qn/qnindex.jl")
include("qn/qnindexset.jl")
include("qn/qnitensor.jl")

#####################################
# Ops to ITensor conversions
#
include("Ops/ops_itensor.jl")

#####################################
# MPS/MPO
#
include("mps/abstractmps.jl")
include("mps/deprecated.jl")
include("mps/mps.jl")
include("mps/mpo.jl")
include("mps/sweeps.jl")
include("mps/abstractprojmpo.jl")
include("mps/projmpo.jl")
include("mps/diskprojmpo.jl")
include("mps/projmposum.jl")
include("mps/projmps.jl")
include("mps/projmpo_mps.jl")
include("mps/observer.jl")
include("mps/dmrg.jl")

#####################################
# Physics
#
include("physics/sitetype.jl")
include("physics/lattices.jl")
include("physics/site_types/generic_sites.jl")
include("physics/site_types/qubit.jl")
include("physics/site_types/spinhalf.jl")
include("physics/site_types/spinone.jl")
include("physics/site_types/fermion.jl")
include("physics/site_types/electron.jl")
include("physics/site_types/tj.jl")
include("physics/site_types/qudit.jl") # EXPERIMENTAL
include("physics/site_types/boson.jl") # EXPERIMENTAL
include("physics/fermions.jl")
include("physics/autompo.jl")

#####################################
# Ops to MPO conversions
#
include("Ops/ops_mpo.jl")

#####################################
# Trotter-Suzuki decomposition
#
include("Ops/trotter.jl")

#####################################
# ITensorChainRules
#
include("ITensorChainRules/ITensorChainRules.jl")

#####################################
# ITensorNetworkMaps
#
include("ITensorNetworkMaps/ITensorNetworkMaps.jl")

#####################################
# ITensorVisualizationCore
#
include("ITensorVisualizationCore/ITensorVisualizationCore.jl")
using .ITensorVisualizationCore

#####################################
# Deprecations
#
include("deprecated.jl")

#####################################
# Argument parsing
#
include("argsdict/argsdict.jl")

#####################################
# Package compilation
#
include("packagecompile/compile.jl")

#####################################
# Developer tools, for internal
# use only
#
include("developer_tools.jl")

function __init__()
  return resize!(empty!(INDEX_ID_RNGs), Threads.nthreads()) # ensures that we didn't save a bad object
end

#####################################
# Precompile certain functions
#
#if Base.VERSION >= v"1.4.2"
#  include("precompile.jl")
#  _precompile_()
#end

end # module ITensors


text:
# Ari Huttunen # HIDE
using Pkg # HIDE
Pkg.activate(@__DIR__) # HIDE # ?confused about global dir or relative dir, this works for now running from the dir.
Pkg.add(["Makie","CairoMakie","PalmerPenguins", "Colors", "Chain", "DataFrames", "GLM"]) # HIDE
using CairoMakie, PalmerPenguins, Colors, Chain, DataFrames, GLM
CairoMakie.activate!() # HIDE
let
    ENV["DATADEPS_ALWAYS_ACCEPT"] = "true" # Don't ask if the datasets should be downloaded
	penguins = dropmissing(DataFrame(PalmerPenguins.load()))
    set_theme!(Theme(
        Axis = (
            leftspinevisible = false,
            rightspinevisible = false,
            bottomspinevisible = false,
            topspinevisible = false,
            xticksvisible = false,
            yticksvisible = false,
        ),
    ))

    # These colors were taken from https://allisonhorst.github.io/palmerpenguins/
	adeliecolor = colorant"#FC7808"
	chinstrapcolor = colorant"#8C00EC"
	gentoocolor = colorant"#107A78"
	(adeliecolor, chinstrapcolor, gentoocolor)

    style = Dict(
        "Adelie" => (
            color=(adeliecolor,0.7),
            strokecolor=adeliecolor,
            label="Adelie",
            marker=:circle,
            markersize=20,
			linewidth = 5
        ),
        "Chinstrap" => (
            color=(chinstrapcolor,0.7),
            strokecolor=chinstrapcolor,
            label="Chinstrap",
            marker=:utriangle,
            markersize=20,
			linewidth = 5
        ),
        "Gentoo" => (
            color=(gentoocolor,0.7),
            strokecolor=gentoocolor,
            label="Gentoo",
            marker=:rect,
            markersize=20,
			linewidth = 5
        )
    )

	fig = Figure(resolution = (1000, 700))
	ax1 = fig[1, 1] = Axis(fig, title = "Flipper and bill length",
		xlabel="Flipper length (mm)", ylabel="Bill length (mm)")
	ax1.xticks = 160:10:240
	ax1.yticks = 30:5:60

	specie = @chain penguins begin
		filter(:species => ==("Adelie"), _)
	end
	X = specie[!,:flipper_length_mm]
	Y = specie[!,:bill_length_mm]
	scatter!(ax1, X, Y; style["Adelie"]...)
	ols = lm(@formula(Y ~ X), DataFrame(X=X, Y=Y))
	y = predict(ols)
	lines!(ax1, X, y; style["Adelie"]...)

	chinstrap = @chain penguins begin
		filter(:species => ==("Chinstrap"), _)
	end
	X = chinstrap[!,:flipper_length_mm]
	Y = chinstrap[!,:bill_length_mm]
	scatter!(ax1, X,Y; style["Chinstrap"]...)
	ols = lm(@formula(Y ~ X), DataFrame(X=X, Y=Y))
	y = predict(ols)
	lines!(ax1, X, y; style["Chinstrap"]...)

	gentoo = @chain penguins begin
		filter(:species => ==("Gentoo"), _)
	end
	X = gentoo[!,:flipper_length_mm]
	Y = gentoo[!,:bill_length_mm]
	scatter!(ax1,X, Y; style["Gentoo"]...)
	ols = lm(@formula(Y ~ X), DataFrame(X=X, Y=Y))
	y = predict(ols)
	lines!(ax1, X, y; style["Gentoo"]...)
	axislegend("Penguin species", position = :rb)

    save(joinpath(@__DIR__, "output", "penguinsDataSetWithChain.png"), fig)
end

using Pkg # HIDE
Pkg.status() # HIDE
Pkg.activate() # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
G = CImGui

using Printf

DoGui() do
    @cstatic f=Cfloat(0.0) counter=Cint(0) one=false two=false clear_color=copy(CImGuiFrontEnd.default_clear_color) begin
        G.Begin("Hello, world!")  # create a window called "Hello, world!" and append into it.

        TreeNode("First section") do
            G.Text("This is some useful text.")  # display some text
            TreeNode("Secret") do
                G.Text("Not really...")
            end
            TreeNode("Even more secret") do
                G.Text("Nothing to see here")
            end
        end
        WithColumns(2, ["A", "B"]) do
            @c G.Checkbox("Column 1 Box", &one)  # edit bools storing our window open/close state
            G.NextColumn()
            @c G.Checkbox("Column 2 Box", &two)  # edit bools storing our window open/close state
        end

        @SameLine begin
            G.Text("testing")
            G.Button("Button") && (counter += 1)
            G.Text("counter = $counter")
            @c G.Checkbox("Something", &two)
            G.Text("$one $two")
        end

        WithStyleColors(G.ImGuiCol_Text => G.HSV(0.0, 0.9, 0.8)) do
            G.Text(@sprintf("Application average %.3f ms/frame (%.1f FPS)", 1000 / G.GetIO().Framerate, G.GetIO().Framerate))
        end

        G.End()

    end

end


text:
import Base
using FillArrays: Fill
# """
# A power measure is a product of a measure with itself. The number of elements in
# the product determines the dimensionality of the resulting support.

# Note that power measures are only well-defined for integer powers.

# The nth power of a measure μ can be written μ^x.
# """
# PowerMeasure{M,N,D} = ProductMeasure{Fill{M,N,D}}

# function Base.show(io::IO, μ::PowerMeasure)
#     io = IOContext(io, :compact => true)
#     print(io, μ.data.value, " ^ ", size(μ.data))
# end

# function Base.show_unquoted(io::IO, μ::PowerMeasure{M,N,D}, indent::Int, prec::Int) where {M,N,D}
#     io = IOContext(io, :compact => true)
#     if Base.operator_precedence(:^) ≤ prec
#         print(io, "(")
#         show(io, μ.data.value)
#         print(io, ")")
#     else
#         show(io, size(μ.data))
#     end
#     return nothing
# end

export PowerMeasure

const PowerMeasure{F,T,N,A} = ProductMeasure{F,Fill{T,N,A}}

function Base.:^(μ::AbstractMeasure, dims::Integer...)
    return μ^dims
end

function Base.:^(μ::M, dims::NTuple{N,I}) where {M<:AbstractMeasure,N,I<:Integer}
    ProductMeasure(identity, Fill(μ, dims))
end

function Base.show(io::IO, d::PowerMeasure)
    io = IOContext(io, :compact => true)
    print(io, d.f(first(d.pars)), " ^ ", size(d.pars))
end

function Base.show(io::IO, d::PowerMeasure{F,T,1}) where {F,T}
    io = IOContext(io, :compact => true)
    print(io, d.f(first(d.pars)), " ^ ", size(d.pars)[1])
end

# sampletype(d::PowerMeasure{M,N}) where {M,N} = @inbounds Array{sampletype(first(marginals(d))), N}

function Base.:^(μ::WeightedMeasure, dims::NTuple{N,I}) where {N,I<:Integer}
    k = prod(dims) * μ.logweight
    return WeightedMeasure(k, μ.base^dims)
end

params(d::PowerMeasure{D}) where {D} = params(first(marginals(d)))
params(::Type{PowerMeasure{D}}) where {D} = params(D)

# basemeasure(μ::PowerMeasure) = @inbounds basemeasure(first(μ.data))^size(μ.data)


text:
function setup_polishing()
    options = Dict(:verbose => false,
                   :polish => true,
                   :eps_abs => 1e-03,
                   :eps_rel => 1e-03,
                   :verbose => false,
                   :max_iter => 5000)
    return options
end

tol = 1e-3

@testset "polishing" begin

    @testset "polishing_problem" begin
        P = sparse(Diagonal([11.; 0.]))
        q = [3.; 4.]
        A = sparse([-1. 0.; 0. -1.; -1. -3; 2. 5.; 3. 4.])
        u = [0.; 0.; -15.; 100.; 80]
        l = -Inf * ones(length(u))
        (m, n) = size(A)

        # Solve problem
        options = setup_polishing()
        model = OSQP.Model()
        OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)
        results = OSQP.solve!(model)

        x_test = [9.90341e-11; 5.0]
        y_test = [1.66667; 0.0; 1.33333; 1.20431e-14; 1.49741e-14]
        obj_test = 20.

        @test isapprox(results.x, x_test, atol=tol)
        @test isapprox(results.y, y_test, atol=tol)
        @test isapprox(results.info.obj_val, obj_test, atol=tol)
        @test results.info.status_polish == 1

    end

    @testset "polishing_unconstrained" begin

        seed!(1)

        n = 10
        m = n
        P = sparse(Diagonal(rand(n)) + .2 * sparse(I, n, n))
        q = randn(n)
        A = sparse(I, n, n)
        l = -100 * ones(m)
        u = 100 * ones(m)

        # Solve problem
        options = setup_polishing()
        model = OSQP.Model()
        OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)
        results = OSQP.solve!(model)

        # Explicit solution
        invP = inv(Array(P))
        x_test = - invP * q
        y_test = zeros(m)
        obj_test = - .5 * q' * invP * q

        @test isapprox(results.x, x_test, atol=tol)
        @test isapprox(results.y, zeros(m), atol=tol)
        @test isapprox(results.info.obj_val, obj_test, atol=tol)
        @test results.info.status_polish == 1
    end

    @testset "polish_random" begin

        # load randomly generated problem with known accurate solution from Mosek
        problem_data = FileIO.load(joinpath(@__DIR__, "problem_data/random_polish_qp.jld2"))
        P = problem_data["P"]; q = problem_data["q"]; A = problem_data["A"]; u = problem_data["u"]; l = problem_data["l"];
        options = setup_polishing()
        model = OSQP.Model()
        OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)
        results = OSQP.solve!(model)

        @test isapprox(results.x, problem_data["x_test"], atol=tol)
        @test isapprox(results.y, problem_data["y_test"], atol=tol)
        @test isapprox(results.info.obj_val, problem_data["obj_test"], atol=tol)
        @test results.info.status_polish == 1

    end
end