text:
"""
```
get_forecast_input_file(m, input_type)
```

Compute the appropriate forecast input filenames for model `m` and
forecast input type `input_type`.

The default input files for each `input_type` can be overriden by adding entries
to the `Dict{Symbol, String}` returned from
`forecast_input_file_overrides(m)`. For example:

```
overrides = forecast_input_file_overrides(m)
overrides[:mode] = \"path/to/input/file.h5\"
```
"""
function get_forecast_input_file(m, input_type)
    overrides = forecast_input_file_overrides(m)
    if haskey(overrides, input_type)
        override_file = overrides[input_type]
        if ispath(override_file)
            return override_file
        else
            error("Input file $override_file does not exist")
        end
    elseif input_type == :subset
        # If input_type == :subset, also check for existence of overrides[:full]
        return get_forecast_input_file(m, :full)
    end

    if input_type == :mode
        return rawpath(m,"estimate","paramsmode.h5")
    elseif input_type == :mean
        return workpath(m,"estimate","paramsmean.h5")
    elseif input_type == :init
        return ""
    elseif input_type in [:full, :subset]
        return rawpath(m,"estimate","mhsave.h5")
    else
        throw(ArgumentError("Invalid input_type: $(input_type)"))
    end
end

"""
```
get_forecast_filename(m, input_type, cond_type, output_var;
    pathfcn = rawpath, forecast_string = "", fileformat = :jld)

get_forecast_filename(directory, filestring_base, input_type, cond_type,
    output_var; forecast_string = "", fileformat = :jld)
```

### Notes

- If `input_type == :subset`, then the `forecast_string` is also appended to the
filenames. If in this case `forecast_string` is empty, `get_forecast_filename` throws
an error.

- In the second method, `directory` should be a string of the form
  `\"\$saveroot/m990/ss2/forecast/raw/\"`. (Note that a `pathfcn` is therefore
  not required.) `filestring_base` should be equivalent to the result of
  `filestring_base(m)`.
"""
function get_forecast_filename(m::AbstractModel, input_type::Symbol,
                               cond_type::Symbol, output_var::Symbol;
                               pathfcn::Function = rawpath,
                               forecast_string::String = "",
                               fileformat::Symbol = :jld)

    # First, we need to make sure we get all of the settings that have been printed to this filestring
    directory = pathfcn(m, "forecast")
    base      = filestring_base(m)

    # Now we can find the filestring we are looking for
    get_forecast_filename(directory, base, input_type, cond_type, output_var;
                          forecast_string = forecast_string, fileformat = fileformat)
end

function get_forecast_filename(directory::String, filestring_base::Vector{String},
                               input_type::Symbol, cond_type::Symbol, output_var::Symbol;
                               forecast_string::String = "", fileformat = :jld)

    # Base file name
    file_name = String("$output_var.$fileformat")

    # Gather all of the file strings into an array
    filestring_addl = get_forecast_filestring_addl(input_type, cond_type,
                                                   forecast_string = forecast_string)

    return savepath(directory, file_name, filestring_base, filestring_addl)
end

function get_forecast_filestring_addl(input_type, cond_type; forecast_string::String = "")

    filestring_addl = Vector{String}()
    push!(filestring_addl, String("para=" * abbrev_symbol(input_type)))
    push!(filestring_addl, String("cond=" * abbrev_symbol(cond_type)))
    if isempty(forecast_string)
        if input_type == :subset
            error("Must supply a nonempty forecast_string if input_type = subset")
        end
    else
        push!(filestring_addl, String("fcid=" * forecast_string))
    end

    filestring_addl
end

"""
```
get_forecast_output_files(m, input_type, cond_type, output_vars;
    forecast_string = "", fileformat = :jld)

get_forecast_output_files(directory, filestring_base, input_type, cond_type,
    output_vars; forecast_string = "", fileformat = :jld)
```

Compute the appropriate output filenames for model `m`, forecast input
type `input_type`, and conditional type `cond_type`, for each output variable in
`output_vars`. Returns a dictionary of file names with one entry for each output_var.

### Arguments

- See `forecast_one` for descriptions of other non-keyword arguments.

### Keyword Arguments

- `forecast_string::String`: subset identifier for when `input_type = :subset`
- `fileformat::Symbol`: file extension, without a period. Defaults to
  `:jld`, though `:h5` is another common option.

### Notes

See `get_forecast_filename` for more information.
"""
function get_forecast_output_files(m::AbstractModel, input_type::Symbol,
                                   cond_type::Symbol, output_vars::Vector{Symbol};
                                   forecast_string::String = "",
                                   fileformat::Symbol = :jld)

    directory = rawpath(m, "forecast")
    base      = filestring_base(m)

    get_forecast_output_files(directory, base, input_type, cond_type, output_vars,
                              forecast_string = forecast_string,
                              fileformat = fileformat)
end

function get_forecast_output_files(directory::String, filestring_base::Vector{String},
                                   input_type::Symbol, cond_type::Symbol, output_vars::Vector{Symbol};
                                   forecast_string::String = "", fileformat::Symbol = :jld)

    output_files    = Dict{Symbol, String}()
    for var in remove_meansbands_only_output_vars(output_vars)
        output_files[var] = get_forecast_filename(directory, filestring_base,
                                                  input_type, cond_type, var,
                                                  forecast_string = forecast_string,
                                                  fileformat = fileformat)
    end

    return output_files
end

"""
```
write_forecast_outputs(m, output_vars, forecast_output_files, forecast_output;
    df = DataFrame(), block_number = Nullable{Int64}(), block_inds = 1:0,
    verbose = :low)
```

Writes the elements of `forecast_output` indexed by `output_vars` to file, given
`forecast_output_files`, which maps `output_vars` to file names.

If `output_vars` contains `:histobs`, data must be passed in as `df`.
"""
function write_forecast_outputs(m::AbstractModel, input_type::Symbol,
                                output_vars::Vector{Symbol},
                                forecast_output_files::Dict{Symbol, String},
                                forecast_output::Dict{Symbol, Array{Float64}};
                                df::DataFrame = DataFrame(),
                                block_number::Nullable{Int64} = Nullable{Int64}(),
                                block_inds::Range{Int64} = 1:0,
                                subset_inds::Range{Int64} = 1:0,
                                verbose::Symbol = :low)

    for var in output_vars
        prod = get_product(var)
        if prod in [:histut, :hist4q, :forecastut, :bddforecastut, :forecast4q, :bddforecast4q]
            # These are computed and saved in means and bands, not
            # during the forecast itself
            continue
        end
        filepath = forecast_output_files[var]
        if isnull(block_number) || get(block_number) == 1
            jldopen(filepath, "w") do file
                write_forecast_metadata(m, file, var)

                if var == :histobs
                    # :histobs just refers to data, so we only write one draw
                    # (as all draws would be the same)
                    @assert !isempty(df) "df cannot be empty if trying to write :histobs"
                    df1 = df[date_mainsample_start(m) .<= df[:date] .<= date_mainsample_end(m), :]
                    data = df_to_matrix(m, df1)
                    write(file, "arr", data)

                else
                    # Otherwise, pre-allocate HDF5 dataset which will contain
                    # all draws
                    if !isnull(block_number) && get(block_number) == 1
                        # Determine forecast output size
                        dims  = get_forecast_output_dims(m, input_type, var; subset_inds = subset_inds)
                        block_size = forecast_block_size(m)
                        chunk_dims = collect(dims)
                        chunk_dims[1] = block_size

                        # Initialize dataset
                        pfile = file.plain
                        HDF5.d_create(pfile, "arr", datatype(Float64), dataspace(dims...), "chunk", chunk_dims)
                    end
                end
            end
        end

        if var != :histobs
            jldopen(filepath, "r+") do file
                if isnull(block_number)
                    write(file, "arr", forecast_output[var])
                else
                    write_forecast_block(file, forecast_output[var], block_inds)
                end
            end
        end

        if VERBOSITY[verbose] >= VERBOSITY[:high]
            println(" * Wrote $(basename(filepath))")
        end
    end
end

"""
```
write_forecast_metadata(m::AbstractModel, file::JldFile, var::Symbol)
```

Write metadata about the saved forecast output `var` to `filepath`.

Specifically, we save dictionaries mapping dates, as well as state, observable,
pseudo-observable, and shock names, to their respective indices in the saved
forecast output array. The saved dictionaries include:

- `date_indices::Dict{Date, Int}`: saved for all forecast outputs except IRFs
- `state_names::Dict{Symbol, Int}`: saved for `var in [:histstates, :forecaststates, :shockdecstates]`
- `observable_names::Dict{Symbol, Int}`: saved for `var in [:forecastobs, :shockdecobs]`
- `observable_revtransforms::Dict{Symbol, Symbol}`: saved identifiers for reverse transforms used for observables
- `pseudoobservable_names::Dict{Symbol, Int}`: saved for `var in [:histpseudo, :forecastpseudo, :shockdecpseudo]`
- `pseudoobservable_revtransforms::Dict{Symbol, Symbol}`: saved identifiers for reverse transforms used for pseudoobservables
- `shock_names::Dict{Symbol, Int}`: saved for `var in [:histshocks, :forecastshocks, :shockdecstates, :shockdecobs, :shockdecpseudo]`

Note that we don't save dates or transformations for impulse response functions.
"""
function write_forecast_metadata(m::AbstractModel, file::JLD.JldFile, var::Symbol)

    prod  = get_product(var)
    class = get_class(var)

    # Write date range
    if prod != :irf
        dates = if prod == :hist
            quarter_range(date_mainsample_start(m), date_mainsample_end(m))
        elseif prod in [:forecast, :bddforecast]
            quarter_range(date_forecast_start(m), date_forecast_end(m))
        elseif prod in [:shockdec, :dettrend, :trend]
            quarter_range(date_shockdec_start(m), date_shockdec_end(m))
        end

        date_indices = Dict(d::Date => i::Int for (i, d) in enumerate(dates))
        write(file, "date_indices", date_indices)
    end

    # Write state names
    if class == :states
        state_indices = merge(m.endogenous_states, m.endogenous_states_augmented)
        @assert length(state_indices) == n_states_augmented(m) # assert no duplicate keys
        write(file, "state_indices", state_indices)
        rev_transforms = Dict{Symbol,Symbol}(x => Symbol("DSGE.identity") for x in keys(state_indices))
        write(file, "state_revtransforms", rev_transforms)
    end

    # Write observable names and transforms
    if class == :obs
        write(file, "observable_indices", m.observables)
        rev_transforms = if prod != :irf
            Dict{Symbol,Symbol}(x => Symbol(m.observable_mappings[x].rev_transform) for x in keys(m.observables))
        else
            Dict{Symbol,Symbol}(x => Symbol("DSGE.identity") for x in keys(m.observables))
        end
        write(file, "observable_revtransforms", rev_transforms)
    end

    # Write pseudo-observable names and transforms
    if class == :pseudo
        write(file, "pseudoobservable_indices", m.pseudo_observables)
        rev_transforms = if prod != :irf
            Dict{Symbol,Symbol}(x => Symbol(m.pseudo_observable_mappings[x].rev_transform)
                                for x in keys(m.pseudo_observables))
        else
            Dict{Symbol,Symbol}(x => Symbol("DSGE.identity") for x in keys(m.pseudo_observables))
        end
        write(file, "pseudoobservable_revtransforms", rev_transforms)
    end

    # Write shock names and transforms
    if class in [:shocks, :stdshocks] || prod in [:shockdec, :irf]
        write(file, "shock_indices", m.exogenous_shocks)
        if class in [:shocks, :stdshocks]
            rev_transforms = Dict{Symbol,Symbol}(x => Symbol("DSGE.identity") for x in keys(m.exogenous_shocks))
            write(file, "shock_revtransforms", rev_transforms)
        end
    end
end

"""
```
write_forecast_block(file::JLD.JldFile, arr::Array, block_number::Int,
    block_inds::Range{Int64})
```

Writes `arr` to the subarray of `file` indicated by `block_inds`.
"""
function write_forecast_block(file::JLD.JldFile, arr::Array,
                              block_inds::Range{Int64})
    pfile = file.plain
    dataset = HDF5.d_open(pfile, "arr")
    ndims = length(size(dataset))
    dataset[block_inds, fill(Colon(), ndims-1)...] = arr
end

"""
```
read_forecast_metadata(file::JLD.JldFile)
```

Read metadata from forecast output files. This includes dictionaries mapping
dates, as well as state, observable, pseudo-observable, and shock names, to
their respective indices in the saved forecast output array. Depending on the
`output_var`, the saved dictionaries might include:

- `date_indices::Dict{Date, Int}`: not saved for IRFs
- `state_indices::Dict{Symbol, Int}`
- `observable_indices::Dict{Symbol, Int}`
- `pseudoobservable_indices::Dict{Symbol, Int}`
- `shock_indices::Dict{Symbol, Int}`
- `state_revtransforms::Dict{Symbol, Symbol}`: states are not transformed, so
  all values are `:identity`
- `observable_revtransforms::Dict{Symbol, Symbol}`
- `pseudoobservable_revtransforms::Dict{Symbol, Symbol}`
- `shock_revtransforms::Dict{Symbol, Symbol}`: shocks are not transformed, so
  all values are `:identity`
"""
function read_forecast_metadata(file::JLD.JldFile)
    metadata = Dict{Symbol, Any}()
    for field in setdiff(names(file), "arr")
        metadata[Symbol(field)] = read(file, field)
    end
    return metadata
end

function read_forecast_output(m::AbstractModel, input_type::Symbol, cond_type::Symbol,
                              output_var::Symbol, var_name::Symbol,
                              shock_name::Nullable{Symbol} = Nullable{Symbol}();
                              forecast_string::String = "")
    # Get filename
    filename = get_meansbands_input_file(m, input_type, cond_type, output_var,
                                         forecast_string = forecast_string)

    # Determine class and product
    product = get_product(output_var)
    class   = get_class(output_var)

    jldopen(filename, "r") do file
        # Read forecast output
        fcast_series = if isnull(shock_name)
            read_forecast_series(file, class, product, var_name)
        else
            read_forecast_series(file, class, product, var_name, get(shock_name))
        end

        # The `fcast_output` for trends only is of size `ndraws` x `nvars`. We
        # need to use `repeat` below because population adjustments will be
        # different in each period. Now we have something of size `ndraws` x
        # `nvars` x `nperiods`
        if product == :trend
            nperiods = length(read(file, "date_indices"))
            fcast_series = repeat(fcast_series, outer = [1, nperiods])
        end

        # Parse transform
        class_long = get_class_longname(class)
        transforms = read(file, string(class_long) * "_revtransforms")
        transform = parse_transform(transforms[var_name])

        fcast_series, transform
    end
end

"""
```
read_forecast_series(file, class, product, var_name[, shock_name])
```

Read only the forecast output for a particular variable (e.g. for a particular
observable) and possibly a particular shock. Result should be a matrix of size
`ndraws` x `nperiods`.
"""
function read_forecast_series(file::JLD.JldFile, class::Symbol, product::Symbol, var_name::Symbol)
    # Get index corresponding to var_name
    class_long = get_class_longname(class)
    indices = read(file, "$(class_long)_indices")
    var_ind = indices[var_name]

    pfile = file.plain
    filename = pfile.filename
    dataset = HDF5.d_open(pfile, "arr")
    ndims = length(size(dataset))

    # Trends are ndraws x nvars
    if product == :trend
        if ndims == 1 # one draw
            arr = h5read(filename, "arr", (var_ind,))
            arr = reshape(arr, (1, 1))
        elseif ndims == 2 # many draws
            arr = h5read(filename, "arr", (Colon(), var_ind))
        end

    # Other products are ndraws x nvars x nperiods
    elseif product in [:hist, :histut, :hist4q, :forecast, :forecastut, :forecast4q,
                       :bddforecast, :bddforecastut, :bddforecast4q, :dettrend]
        inds_to_read = if ndims == 2 # one draw
            arr = h5read(filename, "arr", (var_ind, Colon()))
        elseif ndims == 3 # many draws
            arr = h5read(filename, "arr", (Colon(), var_ind, Colon()))
            arr = squeeze(arr, 2)
        end
    else
        error("Invalid product: $product for this method")
    end

    return arr
end

function read_forecast_series(file::JLD.JldFile, class::Symbol, product::Symbol, var_name::Symbol,
                              shock_name::Symbol)
    # Get indices corresponding to var_name and shock_name
    class_long = get_class_longname(class)
    indices = read(file, "$(class_long)_indices")
    var_ind = indices[var_name]
    shock_indices = read(file, "shock_indices")
    shock_ind = shock_indices[shock_name]

    pfile = file.plain
    filename = pfile.filename
    dataset = HDF5.d_open(pfile, "arr")
    ndims = length(size(dataset))

    if ndims == 3 # one draw
        arr = h5read(filename, "arr", (var_ind, Colon(), shock_ind))
        arr = squeeze(arr, (1, 3))
        arr = reshape(arr, (1, length(arr)))
    elseif ndims == 4 # many draws
        arr = h5read(filename, "arr", (Colon(), var_ind, Colon(), shock_ind))
        arr = squeeze(arr, (2, 4))
    end

    return arr
end


text:
function mm_spec_kj2_dimspec()::MPCCDimSpec
    dimspec = MPCCDimSpec(2, 1, 2, 0, 0, 1, 0)
    return dimspec
end

function mm_spec_kj2_defn(x, pr, ps)

    f = (x[1] - pr[1])^2 + x[2]^3 + x[2]^2

    ce = Vector{Num}()

    ci = Vector{Num}()

    F = Matrix{Num}(undef, 2, 1)
    F[1,1] = x[1]
    F[2,1] = x[2]

    return MPCCDefinition(f, ce, ci, F)
end

function mm_spec_kj2_nzmask()::MPCCModelNZMask
    nz_mask = MPCCModelNZMask(
        # ce
        [],
        # ci
        [   Set{Int64}([1, 2]),
            Set{Int64}([1]) ],
        # F
        [
            # F1
            [   Set{Int64}([1]) ],
            # F2
            [   Set{Int64}([2]) ]
        ] )
    return nz_mask
end

function mm_spec_kj2_testvectors()
    return Vector{MPCCModelTestVector}(undef, 0)
end

function mm_spec_kj2_knownsols()

    function local_kj2_knownsol(pr::Vector{R}, ps::Vector{Int64}) where {R <: Real}
        @assert ( 1 == length(pr) && 0 == length(ps) )
        t = pr[1]

        if ( t < 0 )
            x_sol = [0, 0]
            f_sol = (-t)^2
        else
            x_sol = [t, 0]
            f_sol = 0
        end

        return (x_sol, f_sol)
    end

    return [ local_kj2_knownsol ]
end

ext: jl


text:
#=
Abstract types
abstract type «name» end
abstract type «name» <: «supertype» end

List of some abstract types and their supertypes in Julia
abstract type Number end
abstract type Real     <: Number end
abstract type AbstractFloat <: Real end
abstract type Integer  <: Real end
abstract type Signed   <: Integer end
abstract type Unsigned <: Integer end

Primitive types
primitive type «name» «bits» end
primitive type «name» <: «supertype» «bits» end

List of primitive types and their supertypes in Julia
primitive type Float16 <: AbstractFloat 16 end
primitive type Float32 <: AbstractFloat 32 end
primitive type Float64 <: AbstractFloat 64 end

primitive type Bool <: Integer 8 end
primitive type Char <: AbstractChar 32 end

primitive type Int8    <: Signed   8 end
primitive type UInt8   <: Unsigned 8 end
primitive type Int16   <: Signed   16 end
primitive type UInt16  <: Unsigned 16 end
primitive type Int32   <: Signed   32 end
primitive type UInt32  <: Unsigned 32 end
primitive type Int64   <: Signed   64 end
primitive type UInt64  <: Unsigned 64 end
primitive type Int128  <: Signed   128 end
primitive type UInt128 <: Unsigned 128 end
=#

primitive type Byte 8 end
# The function reinterpret is used to store the bits of an unsigned integer with 8 bits (UInt8) into the byte.
Byte(val::UInt8) = reinterpret(Byte, val)
b = Byte(0x01)
println(b)

# Parametric typing i.e. Generics
struct Point{T<:Real}
    x::T
    y::T
end

println(Point(0.0, 0.0))
println(Point(1, 3))

struct Pair{T1, T2}
    x::T1
    y::T2
end

println(Pair(3, 1.5))

# Unions
IntOrString = Union{Int64, String}
struct Test
    x::IntOrString
end

println(Test("abc"))
println(Test(123))
# This is invalid:
# println(Test('a'))

a = Vector{T} where T<:Real
# In REPL, type ?a
isintpoint(p::Point{T}) where {T} = (T === Int64)
p = Point(1, 2)
println(isintpoint(p))


text:
using Test
using Statistics
using Random
using LinearAlgebra
using ZigZagBoomerang
using ZigZagBoomerang: poisson_time, sep
const CI = true

Random.seed!(1)

include("poisson.jl")
include("test1d.jl")
include("maintest.jl")
include("testiter.jl")
include("testparallel.jl")
include("staticarrays.jl")

include("sticky.jl")

#include("forwarddiff.jl")


text:
const pyjltypetype = pynew()

function pyjltype_getitem(self::Type, k_)
    if pyistuple(k_)
        k = pyconvert(Vector{Any}, k_)
        Py(self{k...})
    else
        k = pyconvert(Any, k_)
        Py(self{k})
    end
end

function init_jlwrap_type()
    jl = pyjuliacallmodule
    filename = "$(@__FILE__):$(1+@__LINE__)"
    pybuiltins.exec(pybuiltins.compile("""
    class TypeValue(AnyValue):
        __slots__ = ()
        __module__ = "juliacall"
        def __getitem__(self, k):
            return self._jl_callmethod($(pyjl_methodnum(pyjltype_getitem)), k)
        def __setitem__(self, k, v):
            raise TypeError("not supported")
        def __delitem__(self, k):
            raise TypeError("not supported")
    """, filename, "exec"), jl.__dict__)
    pycopy!(pyjltypetype, jl.TypeValue)
end

pyjl(v::Type) = pyjl(pyjltypetype, v)


text:
type ScanGreedyHillClimbing <: GraphSearchStrategy
    cache::ScoreComponentCache
    max_n_parents::Int
    max_depth::Int
    prior::DirichletPrior

    function ScanGreedyHillClimbing(
        cache::ScoreComponentCache;
        max_n_parents::Int=3,
        max_depth::Int=1,
        prior::DirichletPrior=UniformPrior(1.0),
        )

        new(cache, max_n_parents, max_depth, prior)
    end
end

function greedy_score(score_components, n, prior_parent_list, datamat, params::ScanGreedyHillClimbing, ncategories::Vector{Int})
    parent_list = prior_parent_list
    while true

        # Compute score
        best_diff = 0.0
        best_parent_list = parent_list
        for i in 1:n

            # 1) add an edge (j->i)
            if length(parent_list[i]) < params.max_n_parents
                for j in deleteat!(collect(1:n), parent_list[i])
                    if adding_edge_preserves_acyclicity(parent_list, j, i)
                        new_parents = sort!(push!(copy(parent_list[i]), j))
                        #println("A")
                        new_component_score = bayesian_score_component(i, new_parents, ncategories, datamat, params.prior, params.cache)
                        #println("B")
                        if new_component_score - score_components[i] > best_diff
                            best_diff = new_component_score - score_components[i]
                            best_parent_list = deepcopy(parent_list)
                            best_parent_list[i] = new_parents
                        end
                    end
                end
            end

            # 2) remove an edge
            for (idx, j) in enumerate(parent_list[i])

                new_parents = deleteat!(copy(parent_list[i]), idx)
                new_component_score = bayesian_score_component(i, new_parents, ncategories, datamat, params.prior, params.cache)
                if new_component_score - score_components[i] > best_diff
                    best_diff = new_component_score - score_components[i]
                    best_parent_list = deepcopy(parent_list)
                    best_parent_list[i] = new_parents
                end

                # 3) flip an edge
                if length(parent_list[j]) < params.max_n_parents
                    deleteat!(parent_list[i], idx)
                    if adding_edge_preserves_acyclicity(parent_list, i, j)
                        sort!(push!(parent_list[j], i))
                        new_diff = bayesian_score_component(i, copy(parent_list[i]), ncategories, datamat, params.prior, params.cache) - score_components[i]
                        new_diff += bayesian_score_component(j, copy(parent_list[j]), ncategories, datamat, params.prior, params.cache) - score_components[j]
                        if new_diff > best_diff
                            best_diff = new_diff
                            best_parent_list = deepcopy(parent_list)
                        end
                        deleteat!(parent_list[j], findin(parent_list[j], i))
                    end
                    sort!(push!(parent_list[i], j))
                 end
            end
        end

        if best_diff > 0.0
            parent_list = best_parent_list
            score_components = bayesian_score_components(parent_list, ncategories, datamat, params.prior, params.cache)
        else
            break
        end
    end
    sum(score_components), parent_list
end

function Distributions.fit(::Type{DiscreteBayesNet}, data::DataFrame, params::ScanGreedyHillClimbing;
    ncategories::Vector{Int} = map!(i->infer_number_of_instantiations(data[i]), Array{Int}(ncol(data)), 1:ncol(data)),
    )

    n = ncol(data)
    parent_list = map!(i->Int[], Array{Vector{Int}}(n), 1:n)
    datamat = convert(Matrix{Int}, data)'
    score_components = bayesian_score_components(parent_list, ncategories, datamat, params.prior, params.cache)

    # 0 depth
    depth = 0
    best, out_score = greedy_score(score_components, n, parent_list, datamat, params, ncategories)

    # > 1 depth
    greedy_parents = out_score
    while depth < params.max_depth
        depth += 1

        # Scan parameters
        best_parent_list = parent_list
        complete = true
        for i in 1:n

            # 1) add an edge (j->i)
            for j in deleteat!(collect(1:n), parent_list[i])
                if adding_edge_preserves_acyclicity(parent_list, j, i)
                    if length(parent_list[i]) < params.max_n_parents
                        new_parent_list = copy(parent_list)
                        new_parents = sort!(push!(new_parent_list[i], j))
                        new_score, out_score = greedy_score(score_components, n, new_parent_list, datamat, params, ncategories)
                        if new_score > best
                            best = new_score
                            complete = false
                            best_parent_list = new_parent_list
                            greedy_parents = out_score
                        end
                    end
                end
            end

            # 2) did this improve our greedy score?
            # Add the best edge and continue
            parent_list = best_parent_list
        end
        if complete
            break
        end
    end

    # compute the greedy solution
    parent_list = greedy_parents

    # construct the BayesNet
    cpds = Array{DiscreteCPD}(n)
    varnames = names(data)
    for i in 1:n
        name = varnames[i]
        parents = varnames[parent_list[i]]
        cpds[i] = fit(DiscreteCPD, data, name, parents, params.prior, parental_ncategories=ncategories[parent_list[i]], target_ncategories=ncategories[i])
    end
    BayesNet(cpds)
end


text:
"""
A Julia module to perform *o*rthogonal *m*atching *p*ursuit.
OMP.jl includes some simple crossvalidation facilities to find the best predictive model of a given sparsity level. """
module OMP

using IterativeSolvers # need LSQR
using RegressionTools
using PLINK

export omp
export cv_omp

# typealias to only allow single or double precision OMP
typealias Float Union{Float32, Float64}

include("common.jl")
include("cv.jl")
include("gwas.jl")
include("mp.jl")

end


text:
module Sum
using DelimitedFiles

values = readdlm("../input/input.txt")
println(sum(values))

end

max_stars_repo_name: tjpalmer/advent

max_issues_repo_name: tjpalmer/advent

max_forks_repo_name: tjpalmer/advent


text:
"""
module FymLogging

Logging module of FymEnvs.
"""

module FymLogging

import FymEnvs: close!, record

using JLD2
using Dates
# using Debugger

export Logger, close!, record, load, set_info!

############ Logger ############
mutable struct Logger
    path
    mode
    info
    buffer
    len
    max_len
    Logger(args...; kwargs...) = init!(new(), args...; kwargs...)
end

function init!(logger::Logger;
               path=nothing, log_dir=nothing, file_name="data.jld2",
               max_len::Int=1000,
               mode="w")
    if path == nothing
        if log_dir == nothing
            log_dir = joinpath("log", Dates.format(now(), "Ymmd-HMS"))
        end
        mkpath(log_dir)
        path = joinpath(log_dir, file_name)
    end
    logger.path = path
    logger.max_len = max_len
    jldopen(path, mode) do jldfile
        initialise(jldfile, ["data", "info"])
    end
    logger.mode = mode
    logger.info = Dict()
    clear!(logger)
    return logger
end

function initialise(jldfile, group_array)
    JLD2.Group(jldfile, "data")
    JLD2.Group(jldfile, "info")
end

function clear!(logger::Logger)
    logger.buffer = Dict()
    logger.len = 0
end

"Record a dictionary or a numeric data preserving the structure."
function record(logger::Logger, info::Dict)
    _rec_update!(logger.buffer, info)
    logger.len += 1
    if logger.len % logger.max_len == 0
        flush!(logger)
    end
end

function append_data(data, buffer)
    for key in keys(buffer)
        if key in keys(data)
            if typeof(data[key]) <: Dict
                _ = append_data(data[key], buffer[key])
            else
                data[key] = cat(data[key], buffer[key], dims=1)
            end
        else
            data[key] = buffer[key]
        end
    end
    return data
end

function flush!(logger::Logger)
    data = load(logger.path; with_info=false)
    # data appending
    append_data(data, logger.buffer)
    # info update
    info = logger.info  # always overwrited
    jldopen(logger.path, "w") do jldfile
        initialise(jldfile, ["data", "info"])
        _rec_save!(jldfile, "/data/", data)
        _info_save!(jldfile, info)
    end
    clear!(logger)
end

"Close `logger`.
You must close manually defined logger after simulation terminated."
function close!(logger::Logger)
    flush!(logger)
end

function set_info!(logger::Logger, info::Dict)
    _rec_update!(logger.info, info, is_info=true)
end

############ others ############
function _info_save!(jldfile, info::Dict=Dict())
    for (key, val) in info
        jldfile["info"][key] = val
    end
end

"Recursively save the `dic` into the JLD2 file."
function _rec_save!(jldfile, path, dic)
    for (key, val) in dic
        if typeof(val) <: Array
            jldfile[path * key] = val
        elseif typeof(val) <: Dict
            _rec_save!(jldfile, path * key * "/", val)
        else
            error("Cannot save $(typeof(val)) type")
        end
    end
end

"Load JLD2 data saved by `Logger`."
function load(path; with_info=false)
    data = Dict()
    info = Dict()
    jldopen(path, "r") do jldfile
        data = _rec_load(jldfile, "data")
        if with_info
            info = _rec_load(jldfile, "info")
        end
    end
    if with_info
        return data, info
    else
        return data
    end
end

function _rec_load(jldfile, path)
    ans = Dict()
    for key in keys(jldfile[path])
        val = jldfile[path * "/" * key]
        if typeof(val) <: JLD2.Group
            ans[key] = _rec_load(jldfile, path * "/" * key)
        else
            ans[key] = val
        end
    end
    return ans
end

"Recursively update `base_dict` with `input_dict`."
function _rec_update!(base_dict, input_dict; is_info=false)
    for (key, val) in input_dict
        if typeof(val) <: Dict
            _rec_update!(get!(base_dict, key, Dict()), val; is_info=is_info)
        else
            if is_info
                get!(base_dict, key, val)
            else
                if !(typeof(val) <: String)
                    get!(base_dict, key, typeof(val[1])[])
                    if length(size(val)) > 0
                        val = reshape(val, 1, size(val)...)
                    end
                    base_dict[key] = cat(base_dict[key], val, dims=1)
                else
                    error("Unsupported data type")
                end
            end
        end
    end
end

end


text:
"""
    sample!( pg, ps, df, mesh)

Sample from a Particle sampler

- `pg`   : Particle group
- `ps`   : Particle sampler
- `df`   : Distribution function
- `xmin` : lower bound of the domain
- `Lx`   : length of the domain.
"""
function sample!(
    pg::ParticleGroup{1,1}, ps::ParticleSampler, df::CosSumGaussianSpin, mesh::AbstractGrid
)
    s = zeros(pg.n_spin)

    theta = 0.0
    phi = 0.0
    n_rnds = 0
    if df.params.n_gaussians > 1
        n_rnds = 1
    end

    δ = zeros(df.params.n_gaussians)
    for i_v in 1:(df.params.n_gaussians)
        δ[i_v] = sum(df.params.δ[1:i_v])
    end

    n_rnds += 2
    rdn = zeros(3)

    if ps.sampling_type == :sobol
        rng_sobol = SobolSeq(1)
    end

    rng_random = MersenneTwister(ps.seed)

    d = Normal()

    for i_part in 1:(pg.n_particles)
        if ps.sampling_type == :sobol
            x = mesh.xmin + Sobol.next!(rng_sobol)[1] * mesh.dimx
        else
            x = mesh.xmin + rand(rng_random) * mesh.dimx
        end

        v = rand(rng_random, d)

        # For multiple Gaussian, draw which one to take
        rnd_no = rdn[3]

        i_gauss = 1
        while (rnd_no > δ[i_gauss])
            i_gauss += 1
        end
        v = v * df.params.σ[i_gauss] + df.params.μ[i_gauss]

        #for a peaked initial condition in the s direction
        s = [0, 0, 1]
        #for a uniformly distributed initial condition on the sphere
        #       for tt = 1:10
        #            s[1] = randn()
        #            s[2] = randn()
        #            s[3] = randn()
        #            if norm(s) > 10^(-4)
        #                break
        #            end
        #        end
        #        s .= s./norm(s)

        # Set weight according to value of perturbation
        w = eval_x_density(df, x) * mesh.dimx

        # Copy the generated numbers to the particle
        set_x!(pg, i_part, x[1])
        set_v!(pg, i_part, v)
        set_spin!(pg, i_part, 1, s[1])
        set_spin!(pg, i_part, 2, s[2])
        set_spin!(pg, i_part, 3, s[3])
        # Set weights.
        set_weights!(pg, i_part, w)
    end
end


text:
module Registrator

using UUIDs, LibGit2, DataStructures

import Base: PkgId
import Pkg: Pkg, TOML, GitTools

const DEFAULT_REGISTRY = "https://github.com/JuliaRegistries/General"
const REGISTRIES = Dict{String,UUID}()

include("slack.jl")
include("builtin_pkgs.jl")
include("register.jl")
include("Server.jl")

end # module


text:
function get_plot_points(plane::Plane{T}; n = 30) where {T <: AbstractFloat}
    plot_points = Vector{CartesianPoint{T}}[]
    vertices = get_vertices(plane)
    for i in 1:length(vertices)
        push!(plot_points, Vector{CartesianPoint{T}}([vertices[i], vertices[i%length(vertices)+1]]))
    end
    plot_points
end

function mesh(tri::Plane{T, CartesianPoint{T}, Nothing}; n = 30) where {T <: AbstractFloat}
    vertices = [tri.p1]
    append!(vertices, get_vertices(tri))
    Mesh(reshape(map(p -> p.x, vertices), (2,2)),reshape(map(p -> p.y, vertices), (2,2)),reshape(map(p -> p.z, vertices), (2,2)))
end

function mesh(quad::Plane{T, CartesianPoint{T}, CartesianPoint{T}}; n = 30) where {T <: AbstractFloat}
    tri1, tri2 = decompose_into_tiangles(quad)
    vertices = get_vertices(tri1)
    push!(vertices, tri2.p1)
    Mesh(reshape(map(p -> p.x, vertices), (2,2)),reshape(map(p -> p.y, vertices), (2,2)),reshape(map(p -> p.z, vertices), (2,2)))
end


text:
@testset "SHIPs Descriptor" begin

@info("-------- TEST 🚢 DESCRIPTOR ---------")
@info("I'm really just checking that the constructor and evaluator run ok")
using SHIPs, JuLIP, Test, ASE
using SHIPs.Descriptors

try
   desc = SHIPDescriptor(:Si, deg = 5, rcut = 5.0)
   at = bulk("Si", cubic=true) * 2
   B1 = descriptors(desc, at)
   println(@test true)
catch
   println(@test false)
end

end

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
export FixedParameterSystem

"""
    FixedParameterSystem(F, p) <: AbstractSystem

Fix a parametrized system `F(x; p)` at `p`, i.e., it is treated as a system without parameters.
"""
mutable struct FixedParameterSystem{S<:AbstractSystem,P<:AbstractVector} <: AbstractSystem
    F::S
    p::P
end

struct FixedParameterSystemCache{SC<:AbstractSystemCache} <: AbstractSystemCache
    cache::SC
end

cache(F::FixedParameterSystem, x) = FixedParameterSystemCache(cache(F.F, x, F.p))

Base.size(F::FixedParameterSystem) = size(F.F)

evaluate!(u, F::FixedParameterSystem, x, c::FixedParameterSystemCache) =
    evaluate!(u, F.F, x, F.p, c.cache)
evaluate(F::FixedParameterSystem, x, c::FixedParameterSystemCache) =
    evaluate(F.F, x, F.p, c.cache)
jacobian!(U, F::FixedParameterSystem, x, c::FixedParameterSystemCache) =
    jacobian!(U, F.F, x, F.p, c.cache)
jacobian(F::FixedParameterSystem, x, c::FixedParameterSystemCache) =
    jacobian(F.F, x, F.p, c.cache)
evaluate_and_jacobian!(u, U, F::FixedParameterSystem, x, c::FixedParameterSystemCache) =
    evaluate_and_jacobian!(u, U, F.F, x, F.p, c.cache)
evaluate_and_jacobian(F::FixedParameterSystem, x, c::FixedParameterSystemCache) =
    evaluate_and_jacobian(F.F, x, F.p, c.cache)

"""
    set_parameters!(F::FixedParameterSystem, p)

Update the parameters `p` of `F`.
"""
set_parameters!(F::FixedParameterSystem, p) = begin
    F.p .= p
    F
end


text:
abstract type AbstractQRLayout <: MemoryLayout end

"""
   QRCompactWYLayout{SLAY,TLAY}()

represents a Compact-WY QR factorization whose
factors are stored with layout SLAY and τ stored with layout TLAY
"""
struct QRCompactWYLayout{SLAY,TLAY} <: AbstractQRLayout end
"""
    QRPackedLayout{SLAY,TLAY}()

represents a Packed QR factorization whose
factors are stored with layout SLAY and τ stored with layout TLAY
"""
struct QRPackedLayout{SLAY,TLAY} <: AbstractQRLayout end

"""
    LULayout{SLAY}()

represents a Packed QR factorization whose
factors are stored with layout SLAY and τ stored with layout TLAY
"""
struct LULayout{SLAY} <: AbstractQRLayout end

MemoryLayout(::Type{<:LinearAlgebra.QRCompactWY{<:Any,MAT}}) where MAT =
    QRCompactWYLayout{typeof(MemoryLayout(MAT)),DenseColumnMajor}()
MemoryLayout(::Type{<:LinearAlgebra.QR{<:Any,MAT}}) where MAT =
    QRPackedLayout{typeof(MemoryLayout(MAT)),DenseColumnMajor}()
MemoryLayout(::Type{<:LinearAlgebra.LU{<:Any,MAT}}) where MAT =
    LULayout{typeof(MemoryLayout(MAT))}()

function materialize!(L::Ldiv{<:QRCompactWYLayout,<:Any,<:Any,<:AbstractVector})
    A,b = L.A, L.B
    ldiv!(UpperTriangular(A.R), view(lmul!(adjoint(A.Q), b), 1:size(A, 2)))
    b
end

function materialize!(L::Ldiv{<:QRCompactWYLayout,<:Any,<:Any,<:AbstractMatrix})
    A,B = L.A, L.B
    ldiv!(UpperTriangular(A.R), view(lmul!(adjoint(A.Q), B), 1:size(A, 2), 1:size(B, 2)))
    B
end

materialize!(L::Ldiv{<:LULayout{<:AbstractColumnMajor},<:AbstractColumnMajor,<:LU{T},<:AbstractVecOrMat{T}}) where {T<:BlasFloat} =
    LAPACK.getrs!('N', L.A.factors, L.A.ipiv, L.B)

function materialize!(L::Ldiv{<:LULayout})
    A,B = L.A,L.B
    _apply_ipiv_rows!(A, B)
    ldiv!(UpperTriangular(A.factors), ldiv!(UnitLowerTriangular(A.factors), B))
end

# Julia implementation similar to xgelsy
function materialize!(Ldv::Ldiv{<:QRPackedLayout,<:Any,<:Any,<:AbstractMatrix{T}}) where T
    A,B = Ldv.A,Ldv.B
    m, n = size(A)
    minmn = min(m,n)
    mB, nB = size(B)
    lmul!(adjoint(A.Q), view(B, 1:m, :))
    mB ≥ n || throw(DimensionMismatch("Number of rows in B must exceed number of columns in A"))
    R = A.R
    @inbounds begin
        if n > m # minimum norm solution
            τ = zeros(T,m)
            for k = m:-1:1 # Trapezoid to triangular by elementary operation
                x = view(R, k, [k; m + 1:n])
                τk = reflector!(x)
                τ[k] = conj(τk)
                for i = 1:k - 1
                    vRi = R[i,k]
                    for j = m + 1:n
                        vRi += R[i,j]*x[j - m + 1]'
                    end
                    vRi *= τk
                    R[i,k] -= vRi
                    for j = m + 1:n
                        R[i,j] -= vRi*x[j - m + 1]
                    end
                end
            end
        end
        ldiv!(UpperTriangular(view(R, :, 1:minmn)), view(B, 1:minmn, :))
        if n > m # Apply elementary transformation to solution
            B[m + 1:mB,1:nB] .= zero(T)
            for j = 1:nB
                for k = 1:m
                    vBj = B[k,j]
                    for i = m + 1:n
                        vBj += B[i,j]*R[k,i]'
                    end
                    vBj *= τ[k]
                    B[k,j] -= vBj
                    for i = m + 1:n
                        B[i,j] -= R[k,i]*vBj
                    end
                end
            end
        end
    end
    return B
end
function materialize!(Ldv::Ldiv{<:QRPackedLayout,<:Any,<:Any,<:AbstractVector{T}}) where T
    ldiv!(Ldv.A, reshape(Ldv.B, length(Ldv.B), 1))
    Ldv.B
end

abstract type AbstractQLayout <: MemoryLayout end
struct QRPackedQLayout{SLAY,TLAY} <: AbstractQLayout end
struct AdjQRPackedQLayout{SLAY,TLAY} <: AbstractQLayout end
struct QRCompactWYQLayout{SLAY,TLAY} <: AbstractQLayout end
struct AdjQRCompactWYQLayout{SLAY,TLAY} <: AbstractQLayout end

MemoryLayout(::Type{<:LinearAlgebra.QRPackedQ{<:Any,S}}) where {S,T} =
    QRPackedQLayout{typeof(MemoryLayout(S)),DenseColumnMajor}()
MemoryLayout(::Type{<:LinearAlgebra.QRCompactWYQ{<:Any,S}}) where {S,T} =
    QRCompactWYQLayout{typeof(MemoryLayout(S)),DenseColumnMajor}()

adjointlayout(::Type, ::QRPackedQLayout{SLAY,TLAY}) where {SLAY,TLAY} = AdjQRPackedQLayout{SLAY,TLAY}()
adjointlayout(::Type, ::QRCompactWYQLayout{SLAY,TLAY}) where {SLAY,TLAY} = AdjQRCompactWYQLayout{SLAY,TLAY}()

copy(M::Lmul{<:AbstractQLayout}) = copyto!(similar(M), M)
mulreduce(M::Mul{<:AbstractQLayout,<:AbstractQLayout}) = MulAdd(M)
mulreduce(M::Mul{<:AbstractQLayout}) = Lmul(M)
mulreduce(M::Mul{<:Any,<:AbstractQLayout}) = Rmul(M)

function copyto!(dest::AbstractArray{T}, M::Lmul{<:AbstractQLayout}) where T
    A,B = M.A,M.B
    if size(dest,1) == size(B,1)
        copyto!(dest, B)
    else
        copyto!(view(dest,1:size(B,1),:), B)
        zero!(@view(dest[size(B,1)+1:end,:]))
    end
    lmul!(A,dest)
end

function copyto!(dest::AbstractArray, M::Ldiv{<:AbstractQLayout})
    A,B = M.A,M.B
    copyto!(dest, B)
    ldiv!(A,dest)
end

materialize!(M::Lmul{LAY}) where LAY<:AbstractQLayout = error("Overload materialize!(::Lmul{$(LAY)})")
materialize!(M::Rmul{LAY}) where LAY<:AbstractQLayout = error("Overload materialize!(::Rmul{$(LAY)})")

materialize!(M::Ldiv{<:AbstractQLayout}) = materialize!(Lmul(M.A',M.B))

materialize!(M::Lmul{<:QRPackedQLayout{<:AbstractColumnMajor,<:AbstractColumnMajor},<:AbstractColumnMajor,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasFloat =
    LAPACK.ormqr!('L','N',M.A.factors,M.A.τ,M.B)

materialize!(M::Lmul{<:QRCompactWYQLayout{<:AbstractColumnMajor,<:AbstractColumnMajor},<:AbstractColumnMajor,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasFloat =
    LAPACK.gemqrt!('L','N',M.A.factors,M.A.T,M.B)

function materialize!(M::Lmul{<:QRCompactWYQLayout{<:AbstractColumnMajor,<:AbstractColumnMajor},<:AbstractRowMajor,<:AbstractMatrix{T},<:AbstractMatrix{T}}) where T<:BlasReal
    LAPACK.gemqrt!('R','T',M.A.factors,M.A.T,transpose(M.B))
    M.B
end
function materialize!(M::Lmul{<:QRCompactWYQLayout{<:AbstractColumnMajor,<:AbstractColumnMajor},<:ConjLayout{<:AbstractRowMajor},<:AbstractMatrix{T},<:AbstractMatrix{T}}) where T<:BlasComplex
    LAPACK.gemqrt!('R','C',M.A.factors,M.A.T,(M.B)')
    M.B
end

function materialize!(M::Lmul{<:QRPackedQLayout})
    A,B = M.A, M.B
    require_one_based_indexing(B)
    mA, nA = size(A.factors)
    mB, nB = size(B,1), size(B,2)
    if mA != mB
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
    end
    Afactors = A.factors
    @inbounds begin
        for k = min(mA,nA):-1:1
            for j = 1:nB
                vBj = B[k,j]
                for i = k+1:mB
                    vBj += conj(Afactors[i,k])*B[i,j]
                end
                vBj = A.τ[k]*vBj
                B[k,j] -= vBj
                for i = k+1:mB
                    B[i,j] -= Afactors[i,k]*vBj
                end
            end
        end
    end
    B
end

### QcB
materialize!(M::Lmul{<:AdjQRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractStridedLayout,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasFloat =
    (A = M.A.parent; LAPACK.ormqr!('L','T',A.factors,A.τ,M.B))
materialize!(M::Lmul{<:AdjQRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractStridedLayout,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasComplex =
    (A = M.A.parent; LAPACK.ormqr!('L','C',A.factors,A.τ,M.B))
materialize!(M::Lmul{<:AdjQRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractStridedLayout,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasFloat =
    (A = M.A.parent; LAPACK.gemqrt!('L','T',A.factors,A.T,M.B))
materialize!(M::Lmul{<:AdjQRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractStridedLayout,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasComplex =
    (A = M.A.parent; LAPACK.gemqrt!('L','C',A.factors,A.T,M.B))
function materialize!(M::Lmul{<:AdjQRPackedQLayout})
    adjA,B = M.A, M.B
    require_one_based_indexing(B)
    A = adjA.parent
    mA, nA = size(A.factors)
    mB, nB = size(B,1), size(B,2)
    if mA != mB
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
    end
    Afactors = A.factors
    @inbounds begin
        for k = 1:min(mA,nA)
            for j = 1:nB
                vBj = B[k,j]
                for i = k+1:mB
                    vBj += conj(Afactors[i,k])*B[i,j]
                end
                vBj = conj(A.τ[k])*vBj
                B[k,j] -= vBj
                for i = k+1:mB
                    B[i,j] -= Afactors[i,k]*vBj
                end
            end
        end
    end
    B
end

## AQ
materialize!(M::Rmul{<:AbstractStridedLayout,<:QRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasFloat =
    LAPACK.ormqr!('R', 'N', M.B.factors, M.B.τ, M.A)
materialize!(M::Rmul{<:AbstractStridedLayout,<:QRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasFloat =
    LAPACK.gemqrt!('R','N', M.B.factors, M.B.T, M.A)
function materialize!(M::Rmul{<:Any,<:QRPackedQLayout})
    A,Q = M.A,M.B
    mQ, nQ = size(Q.factors)
    mA, nA = size(A,1), size(A,2)
    if nA != mQ
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)"))
    end
    Qfactors = Q.factors
    @inbounds begin
        for k = 1:min(mQ,nQ)
            for i = 1:mA
                vAi = A[i,k]
                for j = k+1:mQ
                    vAi += A[i,j]*Qfactors[j,k]
                end
                vAi = vAi*Q.τ[k]
                A[i,k] -= vAi
                for j = k+1:nA
                    A[i,j] -= vAi*conj(Qfactors[j,k])
                end
            end
        end
    end
    A
end

### AQc
materialize!(M::Rmul{<:AbstractStridedLayout,<:AdjQRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasReal =
    (B = M.B.parent; LAPACK.ormqr!('R','T',B.factors,B.τ,M.A))
materialize!(M::Rmul{<:AbstractStridedLayout,<:AdjQRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasComplex =
    (B = M.B.parent; LAPACK.ormqr!('R','C',B.factors,B.τ,M.A))
materialize!(M::Rmul{<:AbstractStridedLayout,<:AdjQRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasReal =
    (B = M.B.parent; LAPACK.gemqrt!('R','T',B.factors,B.T,M.A))
materialize!(M::Rmul{<:AbstractStridedLayout,<:AdjQRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasComplex =
    (B = M.B.parent; LAPACK.gemqrt!('R','C',B.factors,B.T,M.A))
function materialize!(M::Rmul{<:Any,<:AdjQRPackedQLayout})
    A,adjQ = M.A,M.B
    Q = adjQ.parent
    mQ, nQ = size(Q.factors)
    mA, nA = size(A,1), size(A,2)
    if nA != mQ
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)"))
    end
    Qfactors = Q.factors
    @inbounds begin
        for k = min(mQ,nQ):-1:1
            for i = 1:mA
                vAi = A[i,k]
                for j = k+1:mQ
                    vAi += A[i,j]*Qfactors[j,k]
                end
                vAi = vAi*conj(Q.τ[k])
                A[i,k] -= vAi
                for j = k+1:nA
                    A[i,j] -= vAi*conj(Qfactors[j,k])
                end
            end
        end
    end
    A
end

__qr(layout, lengths, A; kwds...) = Base.invoke(qr, Tuple{AbstractMatrix{eltype(A)}}, A; kwds...)
_qr(layout, axes, A; kwds...) = __qr(layout, map(length, axes), A; kwds...)
_qr(layout, axes, A, pivot::P; kwds...) where P = Base.invoke(qr, Tuple{AbstractMatrix{eltype(A)},P}, A, pivot; kwds...)
_qr!(layout, axes, A, args...; kwds...) = error("Overload _qr!(::$(typeof(layout)), axes, A)")
_lu(layout, axes, A; kwds...) = Base.invoke(lu, Tuple{AbstractMatrix{eltype(A)}}, A; kwds...)
_lu(layout, axes, A, pivot::P; kwds...) where P = Base.invoke(lu, Tuple{AbstractMatrix{eltype(A)},P}, A, pivot; kwds...)
_lu!(layout, axes, A, args...; kwds...) = error("Overload _lu!(::$(typeof(layout)), axes, A)")
_cholesky(layout, axes, A, ::Val{false}=Val(false); check::Bool = true) = cholesky!(cholcopy(A); check = check)
_cholesky(layout, axes, A, ::Val{true}; tol = 0.0, check::Bool = true) = cholesky!(cholcopy(A), Val(true); tol = tol, check = check)
_factorize(layout, axes, A) = qr(A) # Default to QR

_factorize(::AbstractStridedLayout, axes, A) = lu(A)
function _lu!(::AbstractColumnMajor, axes, A::AbstractMatrix{T}, pivot::Union{NoPivot, RowMaximum} = RowMaximum();
            check::Bool = true) where T<:BlasFloat
    if pivot === NoPivot()
        return generic_lufact!(A, pivot; check = check)
    end
    lpt = LAPACK.getrf!(A)
    check && checknonsingular(lpt[3])
    return LU{T,typeof(A)}(lpt[1], lpt[2], lpt[3])
end

# for some reason only defined for StridedMatrix in LinearAlgebra
function getproperty(F::LU{T,<:LayoutMatrix}, d::Symbol) where T
    m, n = size(F)
    if d === :L
        L = tril!(getfield(F, :factors)[1:m, 1:min(m,n)])
        for i = 1:min(m,n); L[i,i] = one(T); end
        return L
    elseif d === :U
        return triu!(getfield(F, :factors)[1:min(m,n), 1:n])
    elseif d === :p
        return ipiv2perm(getfield(F, :ipiv), m)
    elseif d === :P
        return Matrix{T}(I, m, m)[:,invperm(F.p)]
    else
        getfield(F, d)
    end
end

# Cholesky factorization without pivoting (copied from stdlib/LinearAlgebra).

# _chol!. Internal methods for calling unpivoted Cholesky
## BLAS/LAPACK element types
function _chol!(::SymmetricLayout{<:AbstractColumnMajor}, A::AbstractMatrix{<:BlasFloat}, ::Type{UpperTriangular})
    C, info = LAPACK.potrf!('U', A)
    return UpperTriangular(C), info
end
function _chol!(::SymmetricLayout{<:AbstractColumnMajor}, A::AbstractMatrix{<:BlasFloat}, ::Type{LowerTriangular})
    C, info = LAPACK.potrf!('L', A)
    return LowerTriangular(C), info
end

_chol!(_, A, UL) = LinearAlgebra._chol!(A, UL)

function _cholesky!(layout, axes, A::RealHermSymComplexHerm, ::Val{false}; check::Bool = true)
    C, info = _chol!(layout, A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
    check && LinearAlgebra.checkpositivedefinite(info)
    return Cholesky(C.data, A.uplo, info)
end

function _cholesky!(::SymmetricLayout{<:AbstractColumnMajor}, axes, A::AbstractMatrix{<:BlasReal},
    ::Val{true}; tol = 0.0, check::Bool = true)
    AA, piv, rank, info = LAPACK.pstrf!(A.uplo, A.data, tol)
    C = CholeskyPivoted{eltype(AA),typeof(AA)}(AA, A.uplo, piv, rank, tol, info)
    check && chkfullrank(C)
    return C
end

_inv_eye(_, ::Type{T}, axs::NTuple{2,Base.OneTo{Int}}) where T = Matrix{T}(I, map(length,axs)...)
function _inv_eye(A, ::Type{T}, (rows,cols)) where T
    dest = zero!(similar(A, T, (cols,rows)))
    dest[diagind(dest)] .= one(T)
    dest
end

function _inv(layout, axes, A)
    T = eltype(A)
    (rows,cols) = axes
    n = checksquare(A)
    S = typeof(zero(T)/one(T))      # dimensionful
    S0 = typeof(zero(T)/oneunit(T)) # dimensionless
    dest = _inv_eye(A, S0, (cols,rows))
    ldiv!(factorize(convert(AbstractMatrix{S}, A)), dest)
end

macro _layoutfactorizations(Typ)
    esc(quote
        LinearAlgebra.cholesky(A::$Typ, args...; kwds...) = ArrayLayouts._cholesky(ArrayLayouts.MemoryLayout(A), axes(A), A, args...; kwds...)
        LinearAlgebra.cholesky!(A::LinearAlgebra.RealHermSymComplexHerm{<:Real,<:$Typ}, v::Val{false}=Val(false); check::Bool = true) = ArrayLayouts._cholesky!(ArrayLayouts.MemoryLayout(A), axes(A), A, v; check=check)
        LinearAlgebra.cholesky!(A::LinearAlgebra.RealHermSymComplexHerm{<:Real,<:$Typ}, v::Val{true}; check::Bool = true, tol = 0.0) = ArrayLayouts._cholesky!(ArrayLayouts.MemoryLayout(A), axes(A), A, v; check=check, tol=tol)
        LinearAlgebra.qr(A::$Typ, args...; kwds...) = ArrayLayouts._qr(ArrayLayouts.MemoryLayout(A), axes(A), A, args...; kwds...)
        LinearAlgebra.qr!(A::$Typ, args...; kwds...) = ArrayLayouts._qr!(ArrayLayouts.MemoryLayout(A), axes(A), A, args...; kwds...)
        LinearAlgebra.lu(A::$Typ, pivot::Union{ArrayLayouts.NoPivot,ArrayLayouts.RowMaximum}; kwds...) = ArrayLayouts._lu(ArrayLayouts.MemoryLayout(A), axes(A), A, pivot; kwds...)
        LinearAlgebra.lu(A::$Typ{T}; kwds...) where T = ArrayLayouts._lu(ArrayLayouts.MemoryLayout(A), axes(A), A; kwds...)
        LinearAlgebra.lu!(A::$Typ, args...; kwds...) = ArrayLayouts._lu!(ArrayLayouts.MemoryLayout(A), axes(A), A, args...; kwds...)
        LinearAlgebra.factorize(A::$Typ) = ArrayLayouts._factorize(ArrayLayouts.MemoryLayout(A), axes(A), A)
        LinearAlgebra.inv(A::$Typ) = ArrayLayouts._inv(ArrayLayouts.MemoryLayout(A), axes(A), A)
        LinearAlgebra.ldiv!(L::LU{<:Any,<:$Typ}, B) = ArrayLayouts.ldiv!(L, B)
        LinearAlgebra.ldiv!(L::LU{<:Any,<:$Typ}, B::$Typ) = ArrayLayouts.ldiv!(L, B)
    end)
end

macro layoutfactorizations(Typ)
    esc(quote
        ArrayLayouts.@_layoutfactorizations $Typ
        ArrayLayouts.@_layoutfactorizations SubArray{<:Any,2,<:$Typ}
        ArrayLayouts.@_layoutfactorizations LinearAlgebra.RealHermSymComplexHerm{<:Real,<:$Typ}
        ArrayLayouts.@_layoutfactorizations LinearAlgebra.RealHermSymComplexHerm{<:Real,<:SubArray{<:Real,2,<:$Typ}}
    end)
end

function ldiv!(C::Cholesky{<:Any,<:AbstractMatrix}, B::LayoutArray)
    if C.uplo == 'L'
        return ldiv!(adjoint(LowerTriangular(C.factors)), ldiv!(LowerTriangular(C.factors), B))
    else
        return ldiv!(UpperTriangular(C.factors), ldiv!(adjoint(UpperTriangular(C.factors)), B))
    end
end

LinearAlgebra.ldiv!(L::LU{<:Any,<:LayoutMatrix}, B::LayoutVector) = ArrayLayouts.ldiv!(L, B)

if VERSION ≥ v"1.7-"
    @deprecate qr(A::LayoutMatrix, ::Val{true}) qr(A, ColumnNorm())
end