text:
using Parameters
using NLsolve
using ForwardDiff
using LinearAlgebra
using Plots
using QuantEcon

struct solution
    resultmessage::String
    F::Matrix{Float64}
    Q::Any
    xss::Vector{Float64}

    equations
    A::Matrix{Float64}
    B::Matrix{Float64}
    C::Matrix{Float64}
    E::Any
end

# useful helper functions
atleast_2d(a) = fill(a,1,1)
atleast_2d(a::AbstractArray) = ndims(a) == 1 ? reshape(a, :, 1) : a

function ss_residual(x, equations, ϵ_sd, fargs)
    resid = equations(x, x, x, zero(ϵ_sd), ϵ_sd, fargs...)
    return resid
end

function get_ss(equations, xss_init, ϵ_sd, fargs)

    function ss_resid_nlsolve!(x)
        x = ss_residual(x, equations, ϵ_sd, fargs)
    end

    solver_result = nlsolve(ss_resid_nlsolve!, xss_init)
    return solver_result.zero

end

function rendahl_coeffs(equations, xss, shocks_sd, fargs)

    shocks_sd = atleast_2d(shocks_sd)
    shocks_ss = zero(shocks_sd)

    A = ForwardDiff.jacobian(t -> equations(t, xss, xss, shocks_ss, shocks_sd, fargs...), xss)
    B = ForwardDiff.jacobian(t -> equations(xss, t, xss, shocks_ss, shocks_sd, fargs...), xss)
    C = ForwardDiff.jacobian(t -> equations(xss, xss, t, shocks_ss, shocks_sd, fargs...), xss)
    E = ForwardDiff.jacobian(t -> equations(xss, xss, xss, t, shocks_sd, fargs...), shocks_ss)

    return A, B, C, E

end

function solve_system(A, B, C, E, maxiter=1000, tol=1e-6)

    #==
    Solves for P and Q using Rehndal's Algorithm
    ==#

    F0 = zero(A)
    S0 = zero(A)
    error = one(tol) + tol
    iter = 0

    while error > tol && iter <= maxiter

        F1 = -(C * F0 + B) \ A
        S1 = -(A * S0 + B) \ C

        error = maximum(C * F1 * F1  + B * F1 + A)

        F0 .= F1
        S0 .= S1

        iter += 1

    end

    XP = LinearAlgebra.eigen(F0)
    XS = LinearAlgebra.eigen(S0)

    if iter == maxiter
        outmessage = "Convergence Failed. Max Iterations Reached. Error: $error"
    elseif maximum(abs.(XP.values)) > 1.0
        outmessage = "No Stable Solution Exists!"
    elseif maximum(abs.(XS.values)) > 1.0
        outmessage = "Multiple Solutions Exist!"
    else
        outmessage = "Convergence Successful!"
    end

    Q = -(C * F0 + B) \ E

    println(outmessage)

    return F0, Q, outmessage

end

function solve(equations, fargs, xinit; shocks_sd = 1.0)

    xss = get_ss(equations, xinit, shocks_sd, fargs)
    A, B, C, E = rendahl_coeffs(equations, xss, shocks_sd, fargs)
    F0, Q, outmessage = solve_system(A, B, C, E)

    out = solution(outmessage, F0, Q, xss, equations, A, B, C, E)

    return out

end

function simdata(sol::solution, timeperiods::Int64)

    F = sol.F
    Q = sol.Q
    G = Matrix(I, size(F, 1), size(F, 1))

    lss = LSS(F, Q, G)

    X_simul, _ = simulate(lss, timeperiods)

    return X_simul

end

function compute_irfs(sol::solution, timeperiods = 40, shocks = nothing)

    IRF = zeros(size(sol.F, 1), timeperiods)

    if ~isnothing(shocks)
        shocks = atleast_2d(shocks)
        shocks_ss = zero(shocks)
        E = ForwardDiff.jacobian(t -> equations(sol.xss, sol.xss, sol.xss, t, shocks, fargs...), shocks_ss)
        Q = -(sol.C * sol.F + sol.B) \ E
    else
        Q = sol.Q
    end

    IRF[:, 1] = Q
    for t in 2:timeperiods
        IRF[:, t] = sol.F * IRF[:, t-1]
    end

    IRF = IRF .+ sol.xss

    return IRF

end


text:
# Constants

# Physical parameters

## Materials
const n₁ = 1.49  # acrylate
const n₂ = 1.00  # air

## The sides of the block to be carved
const Block_Side = 0.15 # meters
const Block_Height = Block_Side
const Block_Width = Block_Side

## Depth parameters. If we imagine a reference place, the top face (facing the caustics) is at
## Top_Offset above it. The flat face facing the light source is Bottom below it.
const Top_Offset = 0.01 # 1 centimeters
const Bottom_Offset = 0.02 # 2 centimeters

## Optics
const Focal_Length = 1.0  # meters

const Caustics_Long_Side = 0.20 # m
global Caustics_Height = Caustics_Long_Side
global Caustics_Width = Caustics_Long_Side

# Caustics picture
global N_Pixel_Side = 512
global N_Pixel_Height = N_Pixel_Side
global N_Pixel_Width = N_Pixel_Side

global Meters_Per_Pixel = Caustics_Height / N_Pixel_Height

# calculation
const ω = 1.99
# const ω = 2 / (1 + π / sqrt(N_Pixel_Height * N_Pixel_Width))
const N_Iterations_Convergence = 10_000        # CHECK: What is a reasonable number of iterations?

# # Global allocations to avoid excessive re-allocations. USEFUL ???
# global container = zeros(Float64, 1_024, 1_024)
# global divergence_intensity = zeros(Float64, 1_024, 1_024)
# global target_map = zeros(Float64, 1_024, 1_024)

ext: jl


text:
## characters defined by actions/homomorphisms

function _action_class_fun(
    conjugacy_cls::AbstractVector{CCl},
) where {CCl <: AbstractOrbit{<:PermutationGroups.AbstractPerm}}
    vals = Int[PermutationGroups.nfixedpoints(first(cc)) for cc in conjugacy_cls]
    # in general:
    # vals = [tr(matrix_representative(first(cc))) for cc in conjugacy_cls]
    return Characters.ClassFunction(vals, conjugacy_cls)
end

function _action_class_fun(
    conjugacy_cls::AbstractVector{CCl},
) where {CCl <: AbstractOrbit{<:AbstractMatrix}}
    vals = [tr(first(cc)) for cc in conjugacy_cls]
    return Characters.ClassFunction(vals, conjugacy_cls)
end

function _action_class_fun(
    hom::InducedActionHomomorphism{<:ByPermutations},
    conjugacy_cls
)
    vals = Int[PermutationGroups.nfixedpoints(induce(hom, first(cc))) for cc in conjugacy_cls]
    # in general:
    # vals = [tr(matrix_representative(first(cc))) for cc in conjugacy_cls]
    return Characters.ClassFunction(vals, conjugacy_cls)
end

function _action_class_fun(
    hom::InducedActionHomomorphism{<:ByLinearTransformation},
    conjugacy_cls,
)
    vals = [tr(induce(hom, first(cc))) for cc in conjugacy_cls]
    return Characters.ClassFunction(vals, conjugacy_cls)
end

"""
    action_character([::Type{T}, ]conjugacy_cls, tbl::CharacterTable)
Return the character of the representation given by the elements in the
conjugacy classes `conjugacy_cls`.

This corresponds to the classical definition of characters as a traces of the
representation matrices.
"""
action_character(conjugacy_clss, tbl::CharacterTable) =
    action_character(eltype(tbl), conjugacy_clss, tbl)

function action_character(::Type{T}, conjugacy_clss, tbl::CharacterTable) where T
    ac_char = _action_class_fun(conjugacy_clss)
    constituents = Int[dot(ac_char, χ) for χ in irreducible_characters(tbl)]
    return Character{T}(tbl, constituents)
end

function action_character(::Type{T}, conjugacy_cls, tbl::CharacterTable) where T <: Union{AbstractFloat, ComplexF64}
    ac_char = _action_class_fun(conjugacy_cls)

    constituents = [dot(ac_char, χ) for χ in irreducible_characters(tbl)]
    all(constituents) do c
        ac = abs(c)
        abs(ac - round(ac)) < eps(real(T))*length(conjugacy_cls)
    end

    return Character{T}(tbl, round.(Int, abs.(constituents)))
end

"""
    action_character([::Type{T}, ]hom::InducedActionHomomorphism, tbl::CharacterTable)
Return the character of the representation given by the images under `hom` of
elements in `conjugacy_classes(tbl)`.

This corresponds to the classical definition of characters as a traces of the
representation matrices.
"""
action_character(hom::InducedActionHomomorphism, tbl::CharacterTable) =
    action_character(eltype(tbl), hom, tbl)

function action_character(
    ::Type{T},
    hom::InducedActionHomomorphism,
    tbl::CharacterTable
) where T
    ac_char = _action_class_fun(hom, conjugacy_classes(tbl))
    constituents = Int[dot(ac_char, χ) for χ in irreducible_characters(T, tbl)]
    return Character{T}(tbl, constituents)
end


text:
using Documenter
using RecipesPipeline
using Literate

# create literate versions of the source files
filepath = joinpath(@__DIR__, "..", "src")
files = joinpath.(filepath, readdir(filepath))

Literate.markdown.(files, joinpath(@__DIR__, "src", "generated"); documenter = false)

makedocs(
    sitename = "RecipesPipeline.jl",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages = [
        "index.md",
        "Developer manual" => [
            "Public API" => "api.md",
            "Recipes" => "recipes.md"
            ],
        "Reference" => "reference.md",
        "Source code" => joinpath.("generated",
                [
                "RecipesPipeline.md",
                "api.md",
                "user_recipe.md",
                "plot_recipe.md",
                "type_recipe.md",
                "series_recipe.md",
                "group.md",
                "recipes.md",
                "series.md",
                "group.md",
                "utils.md"
            ]
        ),
    ],
    modules = [RecipesPipeline]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/JuliaPlots/RecipesPipeline.jl",
    push_preview = true
)


text:
@testset "HW3" begin
    using DMUStudent
    using DMUStudent.HW3
    using POMDPs
    using DiscreteValueIteration
    using POMDPTesting
    using POMDPPolicies
    using POMDPModelTools
    using Compose

    @time m = DenseGridWorld()
    # @test has_consistent_distributions(m)

    @time m = rand(DenseGridWorld)
    # @test has_consistent_distributions(m)

    @warn("Skipping HW3 evaluate testing for speed")
    # @test HW3.evaluate(RandomSolver(), "hw3").score <= 100.0
    # @test HW3.evaluate((m,s)->rand(actions(m)), "hw3").score <= 100.0

    r = render(DenseGridWorld(size=(60,60)), (s=GWPos(30,30),))
    draw(SVG(tempname()*".svg"), r)

    q = Dict((1,1)=>0.0, (2,1)=>0.0)
    n = Dict((1,1)=>1, (2,1)=>0)
    t = Dict((1,1,2)=>1)
    tv = visualize_tree(q, n, t, 1)
    sprint((io, x)->show(io, MIME("text/html"), x), tv)

    q = Dict((1,1)=>0.0)
    n = Dict((1,1)=>1)
    t = Dict((1,1,2)=>1)
    tv = visualize_tree(q, n, t, 1)
    sprint((io, x)->show(io, MIME("text/html"), x), tv)
end


text:
@with_kw struct FastTree
    tolerance::Float64 = 1.0
end

# This is the main function

function solve(solver::FastTree, problem::Problem)
    center = problem.input.center
    radius = problem.input.radius[1]
    (W, b) = (problem.network.layers[1].weights, problem.network.layers[1].bias)

    input = forward_affine_map(solver, W, b, problem.input)
    lower, upper = low(input), high(input)
    local_lower, local_upper = similar(lower), similar(lower)

    k_1 = size(W, 1)

    C = vcat(W, -W)
    d = zeros(2k_1)

    stack = Vector{Hyperrectangle}(undef, 0)
    push!(stack, input)
    #count = 0
    while !isempty(stack)
        interval = popfirst!(stack)
        reach = forward_network(solver, problem.network, interval)
        #count += 1
        if issubset(reach, problem.output)
            continue
        else
            if get_largest_width(interval) > solver.tolerance
                sections = bisect(interval)
                for i in 1:2
                    local_lower = low(sections[i])
                    local_upper = high(sections[i])
                    d = vcat(local_upper - b, b - local_lower)

                    inter = intersection(problem.input, HPolyhedron(C, d))

                    if isempty(inter) == false
                        push!(stack, sections[i])
                    end
                end
            else
                return BasicResult(:violated)
            end
        end
    end
    #println("DimTree($(solver.tolerance)): $(count)")
    return BasicResult(:holds)
end

function forward_network(solver::FastTree, nnet::Network, input::Hyperrectangle)
    layers = nnet.layers
    act = layers[1].activation
    reach = Hyperrectangle(low = act.(low(input)), high = act.(high(input)))

    for i in 2:length(layers)
        reach = forward_layer(solver, layers[i], reach)
    end
    return reach
end

function forward_layer(solver::FastTree, L::Layer, input::Hyperrectangle)
    (W, b, act) = (L.weights, L.bias, L.activation)
    center = zeros(size(W, 1))
    gamma  = zeros(size(W, 1))
    for j in 1:size(W, 1)
        node = Node(W[j,:], b[j], act)
        center[j], gamma[j] = forward_node(solver, node, input)
    end
    return Hyperrectangle(center, gamma)
end

function forward_node(solver::FastTree, node::Node, input::Hyperrectangle)
    output    = node.w' * input.center + node.b
    deviation = sum(abs.(node.w) .* input.radius)
    βmax = node.act(output + deviation)
    βmin = node.act(output - deviation)
    return ((βmax + βmin)/2, (βmax - βmin)/2)
end

function forward_affine_map(solver::FastTree, W::Matrix, b::Vector, input::Hyperrectangle)
    center = W * input.center + b
    radius = abs.(W) * input.radius
    return Hyperrectangle(center, radius)
end


text:
using BinaryBuilder

name = "PixhawkData"
version = v"1.0.2"
# Collection of sources required to build PixhawkData
sources = [
    "git://github.com/akshayhiregoudar/PixhawkData.git" =>
    "e57a681e234f8feba7cd4eaa006112b74ba593e222fdc37361e8002d7198aeb1"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/
cd PixhawkData

make
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libPixhawkData", :libPixhawkData),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: akshayhiregoudar/PixhawkDataBuilder

max_issues_repo_name: akshayhiregoudar/PixhawkDataBuilder

max_forks_repo_name: akshayhiregoudar/PixhawkDataBuilder


text:
export AbstractFeatureConstructor, create_features, feature_size

"""
    AbstractFeatureCreator
An abstract feature creator, for feature transformation from the states.
"""
abstract type AbstractFeatureConstructor end

"""
    create_features
Actually create the features
"""
function create_features end

@deprecate feature_size Base.size

"""
    feature_size(fc::AbstractFeatureCreator)

Depricated! Use Base.size instead. Get size of feature vector the features assume exists.
"""
function feature_size(fc::AbstractFeatureConstructor)
    throw("Implement feature size for $(typeof(fc))")
end


text:
using Pkg
using GLMakie, Test
using GLMakie.FileIO
using ImageMagick
# ImageIO seems broken on 1.6 ... and there doesn't
# seem to be a clean way anymore to force not to use a loader library?
filter!(x-> x !== :ImageIO, FileIO.sym2saver[:PNG])
filter!(x-> x !== :ImageIO, FileIO.sym2loader[:PNG])

path = normpath(joinpath(dirname(pathof(AbstractPlotting)), "..", "test", "ReferenceTests"))
Pkg.develop(PackageSpec(path = path))
using ReferenceTests
recorded = joinpath(@__DIR__, "recorded")
rm(recorded; force=true, recursive=true); mkdir(recorded)
ReferenceTests.record_tests(recording_dir=recorded)
ReferenceTests.reference_tests(recorded)
# needs GITHUB_TOKEN to be defined
# ReferenceTests.upload_reference_images()
# Run the below, to generate a html to view all differences:
# recorded, ref_images, scores = ReferenceTests.reference_tests(recorded)
# ReferenceTests.generate_test_summary("preview.html", recorded, ref_images, scores)


text:
# Inflación de Percentiles Ponderados

# InflationPercentileWeighted.jl - Definición de percentiles ponderados

"""
    InflationPercentileWeighted <: InflationFunction
Función de inflación para computar el percentil ponderado k.

## Utilización

    function (inflfn::InflationPercentileWeighted)(base::VarCPIBase{T}) where T

Define cómo opera `InflationPercentileWeighted` sobre un objeto de tipo VarCPIBase.

    function (inflfn::InflationPercentileWeighted)(cs::CountryStructure)

### Ejemplo
Cálculo del percentil 70 de la distribución de variaciones intermensuales ponderadas

```julia-repl
julia> percfn = InflationPercentileWeighted(70)
(::InflationPercentileWeighted) (generic function with 5 methods)
julia>percfn(gtdata) #gtdata es de tipo UniformCountryStructure
231-element Vector{Float32}:
 11.189365
 11.571873
 11.738467
 11.552155
 11.763763
  ⋮
  1.8916845
  2.074194
  2.0474315
  2.1219969
  2.2268414
```
"""
Base.@kwdef struct InflationPercentileWeighted <: InflationFunction
    k::Float32
end

# Métodos para crear funciones de inflación a partir de enteros y flotantes
"""
    InflationPercentileWeighted(k::Int) = InflationPercentileWeighted(k = Float32(k) / 100)
    InflationPercentileWeighted(q::T) where {T <: AbstractFloat}
Permite utilizar valores de `k` que no necesariamente son de tipo de punto flotante, como por ejemplo: enteros, fracciones, Float32 y Float64.

# Ejemplo
Se usa 0.70 en lugar de 70

```julia-repl
julia> percfn = InflationPercentileWeighted(0.70)
(::InflationPercentileWeighted) (generic function with 5 methods)

julia>percfn(gtdata)
231-element Vector{Float32}:
 11.189365
 11.571873
 11.738467
 11.552155
 11.763763
  ⋮
  1.8916845
  2.074194
  2.0474315
  2.1219969
  2.2268414
```
Se obtienen los mismos resultados
"""
InflationPercentileWeighted(k::Int) = InflationPercentileWeighted(k = Float32(k) / 100)
function InflationPercentileWeighted(q::T) where {T <: AbstractFloat}
    q < 1.0 && return InflationPercentileWeighted(convert(Float32, q))
    InflationPercentileWeighted(convert(Float32, q) / 100)
end

"""
    measure_name(inflfn::InflationPercentileWeighted)

Indica qué medida se utiliza para una instancia de una función de inflación.

# Ejemplo
```julia-repl
julia> percfn = InflationPercentileWeighted(0.70)
julia> measure_name(percfn)
"Percentil ponderado 70.0"
```
"""
measure_name(inflfn::InflationPercentileWeighted) = "Percentil ponderado " * string(round(100inflfn.k, digits=2))

# Parámetro de la función de inflación
params(inflfn::InflationPercentileWeighted) = (inflfn.k, )

# Las funciones sobre VarCPIBase se resumen en variaciones intermensuales
function (inflfn::InflationPercentileWeighted)(base::VarCPIBase{T}) where T
    # InflationPercentileWeighted k de la distribución de variaciones intermensuales
    k = inflfn.k

    # Obtener el percentil k de la distribución intermensual
    rows = size(base.v, 1)
    w = StatsBase.aweights(base.w)

    k_interm = Vector{T}(undef, rows)
    Threads.@threads for r in 1:rows
        row = @view base.v[r, :]
        k_interm[r] = StatsBase.quantile(row, w, k)

    end

    k_interm
end

function InflationPercentileWeighted(vec::Vector{<:Real})
    x = vec[1]
    return InflationPercentileWeighted(x)
end


text:
"""
  Three-dimensional panel method for high-Reynolds aerodynamics.

  # AUTHORSHIP
    * Author    : Eduardo J. Alvarez
    * Email     : Edo.AlvarezR@gmail.com
    * Created   : Jun 2018 originally as MyPanel.jl
    * License   : MIT License
"""
module FLOWPanel

# ------------ GENERIC MODULES -------------------------------------------------
import Dierckx
import PyPlot
const plt = PyPlot
import ForwardDiff
using Statistics
using LinearAlgebra

# ------------ FLOW LAB MODULES ------------------------------------------------
# GeometricTools from https://github.com/byuflowlab/GeometricTools.jl
import GeometricTools
const gt = GeometricTools

# ------------ GLOBAL VARIABLES AND DATA STRUCTURES ----------------------------
const module_path = splitdir(@__FILE__)[1]      # Path to this module
                                                # Default path to data files
const def_data_path = joinpath(module_path, "..", "docs", "data")
                                            # Default path to airfoil data files
const def_rfl_path = joinpath(def_data_path, "airfoils")

const RType = Union{Float64,                    # Concrete real types
                    Int64,
                    ForwardDiff.Dual{Nothing,Float64,3},
                    ForwardDiff.Dual{Nothing,Int64,3}
                    }

# Structure of implemented fields
const FIELDS = Dict(
  # Declare fields as follow:
  # "MyField" => Dict( "field_type"  => "scalar" or "vector",
  #                     "entry_type" => "node", "cell", or "system")

  # Flag indicating whether the paneled body has been solved
  "solved"    => Dict(  "field_type"  => "scalar",
                        "entry_type"  => "system"),

  # Freestream velocity at every control point
  "Vinf"      => Dict(  "field_type"  => "vector",
                        "entry_type"  => "cell"),

  # Constant source strength at every panel
  "sigma"     => Dict(  "field_type"  => "scalar",
                        "entry_type"  => "cell"),

  # Constant doublet strength at every panel
  "mu"        => Dict(  "field_type"  => "scalar",
                        "entry_type"  => "cell"),

  # Vortex ring strength at every panel
  "Gamma"     => Dict(  "field_type"  => "scalar",
                        "entry_type"  => "cell"),

  # Freestream direction at every trailing edge point
  "D"         => Dict(  "field_type"  => "vector",
                        "entry_type"  => "system"),

  # Strength of wake elements
  "Gammawake" => Dict(  "field_type"  => "scalar",
                        "entry_type"  => "system"),
)

# Discretization parameter type
const ndivstype = Union{Float64, gt.multidisctype, Nothing}

# ------------ HEADERS ---------------------------------------------------------
for header_name in ["solver", "abstractbody", "utils"]
  include("FLOWPanel_"*header_name*".jl")
end

end # END OF MODULE


text:
# This file is pat of RLinearAlgebra.jl
# 1. Specifies type
# 2. Implements sample function
# 3. Exports Type

# using Random, LinearAlgebra

"""
    LinSysVecColDetermCyclic <: LinSysVecColSelect

An immutable structure without any fields. Specifies deterministic cycling through the
columns of a linear system.
"""
struct LinSysVecColDetermCyclic <: LinSysVecColSelect end

# Common sample interface for linear systems
function sample(
    type::LinSysVecColDetermCyclic,
    A::AbstractArray,
    b::AbstractVector,
    x::AbstractVector,
    iter::Int64
)

    col_ind = mod(iter, 1:length(x))

    # Search direction
    v = zeros(length(x))
    v[col_ind] = 1.0

    # Normal equation residual
    res = dot(A[:,col_ind], A * x - b)

    return v, A, res
end
#export LinSysVecColDetermCyclic


text:
@kernel function sor_car_gpu!(
        potential::AbstractArray{T, 4},
        point_types::AbstractArray{PointType, 4},
        volume_weights::AbstractArray{T, 4},
        q_eff_imp::AbstractArray{T, 4},
        q_eff_fix::AbstractArray{T, 4},
        ϵ_r::AbstractArray{T, 3},
        geom_weights::NTuple{3, <:AbstractArray{T, 2}},
        sor_const::AbstractArray{T, 1},
        update_even_points::Bool,
        depletion_handling_enabled::Bool,
        is_weighting_potential::Bool,
        only2d::Bool
    ) where {T}
    # # GPU Indexing
    # linear_idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x # CUDA.jl
    linear_idx = @index(Global) # KernelAbstractions.jl

    eff_size = broadcast(idim -> 2:size(potential, idim)-1, (1, 2, 3))
    if linear_idx <= prod(length.(eff_size))
        i1, i2, i3 = Tuple(CartesianIndices(eff_size)[linear_idx])
        # Comparison to CPU indecies: (Cyl / Car)
        # i3 <-> idx3 / ir / iz
        # i2 <-> idx2 / iφ / iy
        # i1 <-> idx1 / iz / ix
        in3 = i3 - 1
        in2 = i2 - 1
        in1 = nidx(i1, update_even_points, iseven(i2 + i3))
        iz = i3
        inz = in3
        iy = i2
        iny = in2
        ix = i1
        inx = in1
        ixr = get_rbidx_right_neighbour(ix, update_even_points, iseven(i2 + i3))

        rb_tar_idx, rb_src_idx = update_even_points ? (rb_even::Int, rb_odd::Int) : (rb_odd::Int, rb_even::Int)

        pwwzr::T        = geom_weights[3][1, inz]
        pwwzl::T        = geom_weights[3][2, inz]
        pwΔmpz::T       = geom_weights[3][3, inz]
        Δz_ext_inv_r::T = geom_weights[3][4, inz + 1]
        Δz_ext_inv_l::T = geom_weights[3][4, inz]

        pwwyr::T  = geom_weights[2][1, iny]
        pwwyl::T  = geom_weights[2][2, iny]
        pwΔmpy::T = geom_weights[2][3, iny]
        pwΔmpy_pwΔmpz::T = pwΔmpy * pwΔmpz
        Δy_ext_inv_r_pwΔmpz::T  = geom_weights[2][4, iny + 1] * pwΔmpz
        Δy_ext_inv_l_pwΔmpz::T  = geom_weights[2][4, iny]     * pwΔmpz
        Δz_ext_inv_r_pwΔmpy::T = Δz_ext_inv_r * pwΔmpy
        Δz_ext_inv_l_pwΔmpy::T = Δz_ext_inv_l * pwΔmpy

        pwwyr_pwwzr::T = pwwyr * pwwzr
        pwwyr_pwwzl::T = pwwyr * pwwzl
        pwwyl_pwwzr::T = pwwyl * pwwzr
        pwwyl_pwwzl::T = pwwyl * pwwzl

        pwwxr::T        = geom_weights[1][1, inx]
        pwwxl::T        = geom_weights[1][2, inx]
        pwΔmpx::T       = geom_weights[1][3, inx]
        Δx_ext_inv_r::T = geom_weights[1][4, inx + 1]
        Δx_ext_inv_l::T = geom_weights[1][4, inx]

        pwwxr_pwwzr::T = pwwxr * pwwzr
        pwwxl_pwwzr::T = pwwxl * pwwzr
        pwwxr_pwwzl::T = pwwxr * pwwzl
        pwwxl_pwwzl::T = pwwxl * pwwzl
        pwwxl_pwwyr::T = pwwxl * pwwyr
        pwwxr_pwwyr::T = pwwxr * pwwyr
        pwwxl_pwwyl::T = pwwxl * pwwyl
        pwwxr_pwwyl::T = pwwxr * pwwyl

        ϵ_rrr::T = ϵ_r[ inx + 1,  iy, iz]
        ϵ_rlr::T = ϵ_r[ inx + 1, iny, iz]
        ϵ_rrl::T = ϵ_r[ inx + 1,  iy, inz ]
        ϵ_rll::T = ϵ_r[ inx + 1, iny, inz ]
        ϵ_lrr::T = ϵ_r[ inx,  iy,  iz ]
        ϵ_llr::T = ϵ_r[ inx, iny,  iz ]
        ϵ_lrl::T = ϵ_r[ inx,  iy, inz ]
        ϵ_lll::T = ϵ_r[ inx, iny, inz ]

        # right weight in r: wrr
        wxr::T = ϵ_rrr * pwwyr_pwwzr
        wxr    = muladd(ϵ_rlr, pwwyl_pwwzr, wxr)
        wxr    = muladd(ϵ_rrl, pwwyr_pwwzl, wxr)
        wxr    = muladd(ϵ_rll, pwwyl_pwwzl, wxr)
        # left weight in r: wrr
        wxl::T = ϵ_lrr * pwwyr_pwwzr
        wxl    = muladd(ϵ_llr, pwwyl_pwwzr, wxl)
        wxl    = muladd(ϵ_lrl, pwwyr_pwwzl, wxl)
        wxl    = muladd(ϵ_lll, pwwyl_pwwzl, wxl)
        # right weight in φ: wφr
        wyr::T = ϵ_lrr * pwwxl_pwwzr
        wyr    = muladd(ϵ_rrr, pwwxr_pwwzr, wyr)
        wyr    = muladd(ϵ_lrl, pwwxl_pwwzl, wyr)
        wyr    = muladd(ϵ_rrl, pwwxr_pwwzl, wyr)
        # left weight in φ: wφl
        wyl::T = ϵ_llr * pwwxl_pwwzr
        wyl    = muladd(ϵ_rlr, pwwxr_pwwzr, wyl)
        wyl    = muladd(ϵ_lll, pwwxl_pwwzl, wyl)
        wyl    = muladd(ϵ_rll, pwwxr_pwwzl, wyl)
        # right weight in z: wzr
        wzr::T = ϵ_rrr * pwwxr_pwwyr
        wzr    = muladd(ϵ_rlr, pwwxr_pwwyl, wzr)
        wzr    = muladd(ϵ_lrr, pwwxl_pwwyr, wzr)
        wzr    = muladd(ϵ_llr, pwwxl_pwwyl, wzr)
        # left weight in z: wzr
        wzl::T = ϵ_rrl * pwwxr_pwwyr
        wzl    = muladd(ϵ_rll, pwwxr_pwwyl, wzl)
        wzl    = muladd(ϵ_lrl, pwwxl_pwwyr, wzl)
        wzl    = muladd(ϵ_lll, pwwxl_pwwyl, wzl)

        wxr *= pwΔmpy_pwΔmpz * Δx_ext_inv_r
        wxl *= pwΔmpy_pwΔmpz * Δx_ext_inv_l
        wyr *= Δy_ext_inv_r_pwΔmpz * pwΔmpx
        wyl *= Δy_ext_inv_l_pwΔmpz * pwΔmpx
        wzr *= Δz_ext_inv_r_pwΔmpy * pwΔmpx
        wzl *= Δz_ext_inv_l_pwΔmpy * pwΔmpx

        vxr::T = potential[    ixr,     iy,     iz, rb_src_idx]
        vxl::T = potential[ixr - 1,     iy,     iz, rb_src_idx]
        vyr::T = potential[     ix, iy + 1,     iz, rb_src_idx]
        vyl::T = potential[     ix,    iny,     iz, rb_src_idx]
        vzr::T = potential[     ix,     iy, iz + 1, rb_src_idx]
        vzl::T = potential[     ix,     iy,    inz, rb_src_idx]

        new_potential::T = is_weighting_potential ? zero(T) : (q_eff_imp[ix, iy, iz, rb_tar_idx] + q_eff_fix[ix, iy, iz, rb_tar_idx])
        new_potential = muladd( wxr, vxr, new_potential)
        new_potential = muladd( wxl, vxl, new_potential)
        new_potential = muladd( wyr, vyr, new_potential)
        new_potential = muladd( wyl, vyl, new_potential)
        new_potential = muladd( wzr, vzr, new_potential)
        new_potential = muladd( wzl, vzl, new_potential)

        new_potential *= volume_weights[ix, iy, iz, rb_tar_idx]

        old_potential::T = potential[ix, iy, iz, rb_tar_idx]

        new_potential -= old_potential
        new_potential = muladd(new_potential, sor_const[1], old_potential)

        if depletion_handling_enabled
            vmin::T = min(vxr, vxl, vyr, vyl, vzr, vzl)
            vmax::T = max(vxr, vxl, vyr, vyl, vzr, vzl)

            new_point_type = point_types[ix, iy, iz, rb_tar_idx]
            if new_potential < vmin || new_potential > vmax
                new_potential -= q_eff_imp[ix, iy, iz, rb_tar_idx] * volume_weights[ix, iy, iz, rb_tar_idx] * sor_const[1]
                if (point_types[ix, iy, iz, rb_tar_idx] & undepleted_bit == 0 &&
                    point_types[ix, iy, iz, rb_tar_idx] & pn_junction_bit > 0)
                    # point_types[ix, iy, iz, rb_tar_idx] += undepleted_bit
                    new_point_type += undepleted_bit
                end # mark this point as undepleted
            elseif point_types[ix, iy, iz, rb_tar_idx] & undepleted_bit > 0
                new_point_type -= undepleted_bit
            end
            point_types[ix, iy, iz, rb_tar_idx] = new_point_type
        end

        potential[ix, iy, iz, rb_tar_idx]::T = ifelse(point_types[ix, iy, iz, rb_tar_idx] & update_bit > 0, new_potential, old_potential)
    end
end


text:
using Compat, DataFrames, StatsBase
using Compat.Test
using GLM
using Alistair

# PERFORMANCE COMPARISON OPTIONS
const TEST_SIZE_1 = 100000

Profile.clear()

beta = [0.7; 0.45]
X = hcat(ones(TEST_SIZE_1,1), collect(1.0:TEST_SIZE_1))
Y = X * beta + 0.1*randn(TEST_SIZE_1, 1)

@time data = DataFrame(X1=X[:,2], Y1=Y[:])
@time OLS_result = glm(@formula(Y1 ~ X1), data, Normal(), IdentityLink())

@time OLSRegression = OLS(intercept=true, robust=false)
@time result = linregress(X, Y, OLSRegression)


text:
using DataArrays, DataFrames, FixedEffectModels, Base.Test

df = DataFrame(v1 = @data(["ok", NA, "ok", "first"]), v2 = @data([1, 3, 2, 1]), v3 = @data([NA, 1, 1, NA]))

@test group(df[[:v1]]).pool == [1, 2]
@test group(df[[:v1]]).refs == [2, 0, 2, 1]

@test group(df[[:v2]]) == pool(df[:v2])
@test group(df[[:v2]]) == pool(df[:v2])

@test group(df[[:v3]]).pool ==  [1]
@test group(df[[:v3]]).refs ==  [0, 1, 1, 0]

@test group(df).pool == [1]
@test group(df).refs == [0, 0, 1, 0]

# test different syntaxes
df = DataFrame(v1 = pool(collect(1:1000)), v2 = pool(fill(1, 1000)))
@test group(df, [:v1, :v2]) == collect(1:1000)
@test group(df, :v1, :v2) == collect(1:1000)
@test group(df[:v1], df[:v2]) == collect(1:1000)


text:
# Macro for checking arguments
macro check_args(K, param, cond, desc=string(cond))
    quote
        if !($(esc(cond)))
            throw(
                ArgumentError(
                    string(
                        $(string(K)),
                        ": ",
                        $(string(param)),
                        " = ",
                        $(esc(param)),
                        " does not ",
                        "satisfy the constraint ",
                        $(string(desc)),
                        ".",
                    ),
                ),
            )
        end
    end
end

function vec_of_vecs(X::AbstractMatrix; obsdim::Int=2)
    @assert obsdim ∈ (1, 2) "obsdim should be 1 or 2, see docs of kernelmatrix"
    if obsdim == 1
        RowVecs(X)
    else
        ColVecs(X)
    end
end

"""
    ColVecs(X::AbstractMatrix)

A lightweight wrapper for an `AbstractMatrix` which interprets it as a vector-of-vectors, in
which each _column_ of `X` represents a single vector.

That is, by writing `x = ColVecs(X)`, you are saying "`x` is a vector-of-vectors, each of
which has length `size(X, 1)`. The total number of vectors is `size(X, 2)`."

Phrased differently, `ColVecs(X)` says that `X` should be interpreted as a vector
of horizontally-concatenated column-vectors, hence the name `ColVecs`.

```jldoctest
julia> X = randn(2, 5);

julia> x = ColVecs(X);

julia> length(x) == 5
true

julia> X[:, 3] == x[3]
true
```

`ColVecs` is related to [`RowVecs`](@ref) via transposition:
```jldoctest
julia> X = randn(2, 5);

julia> ColVecs(X) == RowVecs(X')
true
```
"""
struct ColVecs{T,TX<:AbstractMatrix{T},S} <: AbstractVector{S}
    X::TX
    function ColVecs(X::TX) where {T,TX<:AbstractMatrix{T}}
        S = typeof(view(X, :, 1))
        return new{T,TX,S}(X)
    end
end

Base.size(D::ColVecs) = (size(D.X, 2),)
Base.getindex(D::ColVecs, i::Int) = view(D.X, :, i)
Base.getindex(D::ColVecs, i::CartesianIndex{1}) = view(D.X, :, i)
Base.getindex(D::ColVecs, i) = ColVecs(view(D.X, :, i))
Base.setindex!(D::ColVecs, v::AbstractVector, i) = setindex!(D.X, v, :, i)

Base.vcat(a::ColVecs, b::ColVecs) = ColVecs(hcat(a.X, b.X))

dim(x::ColVecs) = size(x.X, 1)

pairwise(d::PreMetric, x::ColVecs) = Distances_pairwise(d, x.X; dims=2)
pairwise(d::PreMetric, x::ColVecs, y::ColVecs) = Distances_pairwise(d, x.X, y.X; dims=2)
function pairwise(d::PreMetric, x::AbstractVector, y::ColVecs)
    return Distances_pairwise(d, reduce(hcat, x), y.X; dims=2)
end
function pairwise(d::PreMetric, x::ColVecs, y::AbstractVector)
    return Distances_pairwise(d, x.X, reduce(hcat, y); dims=2)
end
function pairwise!(out::AbstractMatrix, d::PreMetric, x::ColVecs)
    return Distances.pairwise!(out, d, x.X; dims=2)
end
function pairwise!(out::AbstractMatrix, d::PreMetric, x::ColVecs, y::ColVecs)
    return Distances.pairwise!(out, d, x.X, y.X; dims=2)
end

"""
    RowVecs(X::AbstractMatrix)

A lightweight wrapper for an `AbstractMatrix` which interprets it as a vector-of-vectors, in
which each _row_ of `X` represents a single vector.

That is, by writing `x = RowVecs(X)`, you are saying "`x` is a vector-of-vectors, each of
which has length `size(X, 2)`. The total number of vectors is `size(X, 1)`."

Phrased differently, `RowVecs(X)` says that `X` should be interpreted as a vector
of vertically-concatenated row-vectors, hence the name `RowVecs`.

Internally, the data continues to be represented as an `AbstractMatrix`, so using this type
does not introduce any kind of performance penalty.

```jldoctest
julia> X = randn(5, 2);

julia> x = RowVecs(X);

julia> length(x) == 5
true

julia> X[3, :] == x[3]
true
```

`RowVecs` is related to [`ColVecs`](@ref) via transposition:
```jldoctest
julia> X = randn(5, 2);

julia> RowVecs(X) == ColVecs(X')
true
```
"""
struct RowVecs{T,TX<:AbstractMatrix{T},S} <: AbstractVector{S}
    X::TX
    function RowVecs(X::TX) where {T,TX<:AbstractMatrix{T}}
        S = typeof(view(X, 1, :))
        return new{T,TX,S}(X)
    end
end

RowVecs(x::AbstractVector) = RowVecs(reshape(x, :, 1))

Base.size(D::RowVecs) = (size(D.X, 1),)
Base.getindex(D::RowVecs, i::Int) = view(D.X, i, :)
Base.getindex(D::RowVecs, i::CartesianIndex{1}) = view(D.X, i, :)
Base.getindex(D::RowVecs, i) = RowVecs(view(D.X, i, :))
Base.setindex!(D::RowVecs, v::AbstractVector, i) = setindex!(D.X, v, i, :)

Base.vcat(a::RowVecs, b::RowVecs) = RowVecs(vcat(a.X, b.X))

dim(x::RowVecs) = size(x.X, 2)

pairwise(d::PreMetric, x::RowVecs) = Distances_pairwise(d, x.X; dims=1)
pairwise(d::PreMetric, x::RowVecs, y::RowVecs) = Distances_pairwise(d, x.X, y.X; dims=1)
function pairwise(d::PreMetric, x::AbstractVector, y::RowVecs)
    return Distances_pairwise(d, permutedims(reduce(hcat, x)), y.X; dims=1)
end
function pairwise(d::PreMetric, x::RowVecs, y::AbstractVector)
    return Distances_pairwise(d, x.X, permutedims(reduce(hcat, y)); dims=1)
end
function pairwise!(out::AbstractMatrix, d::PreMetric, x::RowVecs)
    return Distances.pairwise!(out, d, x.X; dims=1)
end
function pairwise!(out::AbstractMatrix, d::PreMetric, x::RowVecs, y::RowVecs)
    return Distances.pairwise!(out, d, x.X, y.X; dims=1)
end

# Resolve ambiguity error for ColVecs vs RowVecs. #346
pairwise(d::PreMetric, x::ColVecs, y::RowVecs) = pairwise(d, x, ColVecs(permutedims(y.X)))
pairwise(d::PreMetric, x::RowVecs, y::ColVecs) = pairwise(d, ColVecs(permutedims(x.X)), y)

dim(x) = 0 # This is the passes-by-default choice. For a proper check, implement `KernelFunctions.dim` for your datatype.
dim(x::AbstractVector) = dim(first(x))
dim(x::AbstractVector{<:AbstractVector{<:Real}}) = length(first(x))
dim(x::AbstractVector{<:Real}) = 1

function validate_inputs(x, y)
    if dim(x) != dim(y) # Passes by default if `dim` is not defined
        throw(
            DimensionMismatch(
                "Dimensionality of x ($(dim(x))) not equality to that of y ($(dim(y)))"
            ),
        )
    end
    return nothing
end

function validate_inplace_dims(K::AbstractMatrix, x::AbstractVector, y::AbstractVector)
    validate_inputs(x, y)
    if size(K) != (length(x), length(y))
        throw(
            DimensionMismatch(
                "Size of the target matrix K ($(size(K))) not consistent with lengths of " *
                "inputs x ($(length(x))) and y ($(length(y)))",
            ),
        )
    end
end

function validate_inplace_dims(K::AbstractVector, x::AbstractVector, y::AbstractVector)
    validate_inputs(x, y)
    n = length(x)
    if length(y) != n
        throw(
            DimensionMismatch(
                "Length of input x ($n) not consistent with length of input y " *
                "($(length(y))",
            ),
        )
    end
    if length(K) != n
        throw(
            DimensionMismatch(
                "Length of target vector K ($(length(K))) not consistent with length of " *
                "inputs ($n)",
            ),
        )
    end
end

function validate_inplace_dims(K::AbstractVecOrMat, x::AbstractVector)
    return validate_inplace_dims(K, x, x)
end


text:
using Base.Test

include("fitting.jl")

function zerodcphase!{T}(data::Array{Complex{T},3})
  nchan, nro, npe = size(data)
  mid = ifloor(nro/2) + 1
  for c in 1:nchan, p in 1:npe
    ϕ = angle(data[c, mid, p])
    data[c,:,p] *= exp(-1im*ϕ)
  end
  data
end

function test_zerodcphase!(nro=32, npe=4)
   print("testing zerodcphase! on 1 x $nro x $npe data ... ")
   nchan = 1
   ϕ = Float64[0:nro-1] - floor(nro/2)   # midpoint at int(nro/2)+1
   mid = ifloor(nro/2) + 1
   Δ = rand(npe) # each line gets a different constant shift
   f = Array(Complex128, nchan, nro, npe)
   g = Array(Complex128, nchan, nro, npe)
   for p in 1:npe
     f[1,:,p] = rand(nro).*exp(1im*ϕ)
     g[1,:,p] = f[1,:,p]*exp(1im*Δ[p])
   end
   zerodcphase!(g)
   for p in 1:npe, r in 1:nro
     @test_approx_eq f[1,r,p] g[1,r,p]
   end
   println("pass")
end

function unshift{T}(data::Array{Complex{T},3}; golden=false, ncalib=32, verbose=false)
  # data assumed to be nchan x nro x npe
  # K assumed to be ndim x nro x npe

  if verbose
    println("correcting for radial shifts")
  end
  # calculate sinogram
  sino = fftshift(ifft(fftshift(data, 2), 2), 2)
  nchan, nro, npe = size(data)
  data_shifted = Array(Complex{T}, nchan, nro, npe)
  slopes = zeros(nchan, npe)
  intercepts = zeros(nchan, npe)
  Δk = golden ? 21 : 1
  @assert Δk < npe "not enough projections for meaningful correction"
  if verbose
    println("using Δk of ±$(Δk)")
  end
  ncalib = min(ncalib, div(nro,2)-1)
  roidxs = [div(nro,2)-ncalib:div(nro,2)+ncalib]
  r = [-int(nro/2): int(nro/2)-1]
  for p in 1:npe, j in 1:nchan
    B = vec(data[j,roidxs,p])
    b = vec(sino[j,roidxs,p])
    if golden
      a = p > Δk        ? sino[j,roidxs,p-Δk] : nothing
      c = p <= npe - Δk ? sino[j,roidxs,p+Δk] : nothing
    else  # linear
      a = p == 1   ? sino[j,roidxs,npe] : sino[j,roidxs,p-1]
      c = p == npe ? sino[j,roidxs,1]   : sino[j,roidxs,p+1]
    end
    # average spatial projections a and c
    if a == nothing
      a = c
      ac = c
    elseif c == nothing
      c = a
      ac = a
    else
      ac = 0.5*(a + c)
    end
    a = vec(a)
    b = vec(b)
    c = vec(c)
    ac = vec(ac)

    # create a mask for the spatial domain projections
    abcmag = abs(a) + abs(b) + abs(c)
    abcmag = (abcmag - minimum(abcmag)) / (maximum(abcmag) - minimum(abcmag))
    #abcmag[abcmag .< 0.25] = 0.0  # TODO threshold_otsu(abcmag)

    # calculate masked phase difference
    Δφ = golden ? angle(ac[end:-1:1] ./ b) : angle(ac ./ b)
    Δφ = Δφ .* abcmag

    # perform a weighted least squares fit
    W = diagm(abcmag[:])
    M = hcat(r[roidxs], ones(length(roidxs)))
    u = M'*W*M \ M'*W*Δφ
    slopes[j,p] = u[1]
    intercepts[j,p] = u[2]

    # apply phase compensation
    phase_ramp = 0.5*(u[1]*r + u[2])
    b_phase_compensated = vec(sino[j,:,p]) .* exp(1im*phase_ramp)
    B_shifted = fftshift(fft(fftshift(b_phase_compensated)))
    data_shifted[j,:,p] = B_shifted[:]
  end
  (data_shifted, slopes, intercepts)
end

function test_unshift(nro=32, npe=32)
  print("testing unshift() on 1 x $nro x $npe data ... ")
  nchan = 1
  ϕ = Float64[0:nro-1] - floor(nro/2)   # midpoint at int(nro/2)+1
  mid = ifloor(nro/2) + 1
  Δ = rand(npe) # each line gets a different constant shift
  f = Array(Complex128, nchan, nro, npe)
  f[1,:,:] = repmat(rand(nro).*exp(1im*ϕ), 1, npe)
  g, m, b = unshift(f)
  for p in 1:npe
    @test_approx_eq_eps m[p] 0.0 eps(Float64)
    @test_approx_eq_eps b[p] 0.0 eps(Float64)
    for r in 1:nro
      @test_approx_eq f[1,r,p] g[1,r,p]
    end
  end
  println("pass")
end

function delaycorrect{T,N}(data::Array{Complex{T},3}, K::Array{T,N}; validate=false)
    # radial gradient delay correction

    outdir = "validation"
    if validate
        println("validation enabled: run time will be much longer")
        isdir(outdir) || mkdir(outdir)
    end

    nrep, nro, npe = size(data)
    ndim = ndims(data)

    nparam = size(K,1)  # number of gradient channels
    @assert nparam == 2  "only 2D implemented for now"
    println("npe: $npe\nnro: $nro\nndim: $ndim")
    println("fitting $nparam parameters")

    # Column $k$ of the $U$ matrix contains a unit vector that points along
    # the $k$-th profile This matrix is then squared on an element-by-element
    # basis to yield the relative contributions of each gradient axis to the
    # delay of the profile parallel to the readout direction.

    U = zeros(nparam, npe)   # unit vectors pointing along each radial profile in readout direction
    angles = ((π / φ) * [0:npe-1]) % 2π
    for i in 1:npe
        # spherical polar angles corresponding to radial readout orientations
        θ = (π / φ) * (i - 1)
        U[1,i] = cos(θ)
        U[2,i] = sin(θ)
    end
    U = U.^2

    # The uncorrected r coordinates
    r = [0:nro-1]/nro - 0.5

    Δp = int(npe / 100)  # number of profiles to skip
    println("taking every $(Δp)-th profile")

    shape{S}(x::Vector{S}, p::Vector{S}) = p[1]*abs(sinc(p[2]*(x - p[3]))) + p[4]
    #shape{S}(x::Vector{S}, p::Vector{S}) = p[1] ./ (p[2] + abs(x - p[3]).^p[4])

    datamin, datamax = extrema(abs(data[1,:,:]))
    p0 = [datamax,1.0, 0.0, 1.5]
    yfit = similar(data)
    params, resid = nlsfit(shape, float64(abs(squeeze(data[1,:,:],1))), [1:Δp:npe], r, p0)
    rcenter = float32(squeeze(params[3,:],1))
    println("residual $(norm(resid))")
    for i in 1:Δp:npe
        yfit[1,:,i] = shape(r, params[:,i])
        if validate   # plot profiles and their fits
             clf()
             plot(r,abs(data[1,i,:]),"b.", r, abs(yfit[:,i]), "g-")
             ylim(datamin, datamax)
             xlabel("spatial frequency")
             ylabel("power")
             savefig("$outdir/fit$i.png")
        end
    end
    if false
        clf()
        plot(angles,rcenter,"o")
        ylabel("fitting center freq")
        xlabel("angle (rad)")
        savefig("$outdir/rcenter.png")
        clf()
        plot([1:npe],vec(params[4,:]),"o")
        ylabel("profile fit power law index")
        xlabel("profile")
        savefig("$outdir/powerindex.png")
    end

    # Finally, we formulate and solve the linear inverse problem $U^2 \delta =
    # \Delta r_{||} $ , where $d$ are the deduced three-axis gradient delays,
    # oriented along the coordinate axes, $\Delta r_{||}$ is the vector of center
    # offsets calculated from the Cauchy distribution fits, and $U$ is the matrix of
    # unit vectors pointing along the radial profiles.  Matrix $U$ is squared
    # because one factor computes the relative contribution from each axis, and the
    # other factor accounts for the portion of the coordinate offset that is
    # parallel to the readout direction, which is the only information available to
    # us.

    delta = U[:,1:Δp:end]' \ rcenter[1:Δp:end]   # <-- THE MAGIC
    println("derived delays: $delta")

    # These three values in delta should be the $x$, $y$, and $z$ gradient delays,
    # respectively, in units of spatial frequency. Since this is a massively
    # overdetermined problem, we have the luxury of re-computing these from subsets
    # of the full data set.  One prudent approach would be to solve the inverse
    # problem on a sliding window to examine temporal consistency.  Since we only
    # need at minimum three profiles to invert to obtain the three gradient delays,
    # we will set our sliding window width to 3.

    # To correct $k$-space, we simply apply the calculated delays from $\delta$ in
    # the appropriate direction for each $(k_x,k_y,k_z)$ triplet.
    for j in 1:nro, i in 1:npe
      K[1,j,i] += U[1,i] * delta[1]
      K[2,j,i] += U[2,i] * delta[2]
    end
    K = clamp(K, -0.5, 0.5)
    convert(Array{T,3}, K)
end

function window!{T}(x::Vector{T})
  n = length(x)
  x = x .* sinc([0:n-1]/n)  # Lanczos factors
end
function window!{T}(x::Array{T,2})
  m, n = size(x)
  x = x .* window(m,n)
end
function window(m::Int, n::Int)
  nx = int(m/2)
  ny = int(n/2)
  fx = sinc([0:nx-1]/nx)
  fy = sinc([0:ny-1]/ny)
  vcat(fx[end:-1:1], fx) * vcat(fy[end:-1:1], fy)'
end

function prewhiten!{T,N}(data::Array{T,N}, noise::Array{T,2})
  # Prewhitens data by removing noise correlation between coils.
  # Assumes:
  #   channel dimensions is the first dimension of 'data'
  #   noise is an nchannel x nsamples array of the same type as 'data'
  insize = size(data) # save input size for reshaping upon output
  nsamp = size(noise, 2)
  μ = mean(noise, 2)
  σ = std(noise, 2)
  for j in 1:nsamp  # for each noise sample
    noise[:,j] -= μ # zero mean
    noise[:,j] ./= σ #  unit std dev
  end
  C = noise*noise' / (nsamp-1)  # C = noise covariance matrix
  L = chol(C, :L)  # Cholesky decomposition => lower triangular result
  data = full(inv(L)) * data[:,:] # decorrelate the data using the noise correlations
       # ^ full() needed because * not yet defined for Triangular and Complex arrays
  reshape(data, insize)
end

function test_prewhiten!(m=32, n=16)
  srand(1337)
  x = 1000*rand(n,m)+1000im*rand(n,m)
  y = copy(x)
  prewhiten!(y, x)
  μ = mean(y,2)
  σ = std(y,2)
  for j in 1:m
    y[:,j] -= μ
    y[:,j] ./= σ
  end
  Σ = abs(y*y') / (m-1)
  @test_approx_eq_eps Σ eye(n) 1e-1
end

function runtests()
  test_zerodcphase!()
  test_unshift()
  test_prewhiten()
end

max_stars_repo_path: code/Julia Recon Scripts/correct.jl

max_issues_repo_path: code/Julia Recon Scripts/correct.jl

max_forks_repo_path: code/Julia Recon Scripts/correct.jl


text:
@inline function improved_estimator(model, param::Parameter, extra)
    p = convert_parameter(model, param)
    return improved_estimator(model, p..., extra)
end

@doc """
    improved_estimator(model::Ising, T::Real, Js::AbstractArray, sw::SWInfo)

Returns the following observables as `Dict{String, Any}` using cluster information `sw`

# Observables
- `"Energy"`
    - Energy per spin (site)
- `"Energy^2"`
- `"Magnetization"`
    - Total magnetization per spin (site)
- `"|Magnetization|"`
- `"|Magnetization|^2"`
- `"|Magnetization|^4"`
"""
function improved_estimator(model::Ising, T::Real, Js::AbstractArray, sw::SWInfo)
    nsites = numsites(model)
    nbonds = numbonds(model)
    nc = numclusters(sw)
    invV = 1.0/nsites

    ## magnetization
    M = 0.0
    M2 = 0.0
    M4 = 0.0
    for (m,s) in zip(sw.clustersize, sw.clusterspin)
        M += m*invV*s
        m2 = (m*invV)^2
        M4 += m2*m2 + 6M2*m2
        M2 += m2
    end

    # energy
    aJ = 2.0*abs.(Js)
    mbeta = -1.0/T
    ns = sw.activated_bonds
    As = -aJ ./ expm1.(mbeta.*aJ)
    Ans = ns.*As
    E0 = 0.0
    for b in 1:numbondtypes(model)
        E0 += Js[b] * numbonds(model,b)
    end
    E = 0.0
    E2 = 0.0
    for b in 1:numbondtypes(model)
        E2 += (aJ[b]-2.0*E0)*Ans[b]
        E2 += Ans[b] * As[b]*(ns[b]-1)
        E2 += 2.0*Ans[b]*E
        E += Ans[b]
    end
    E -= E0
    E2 += E0^2

    E *= -invV
    E2 *= invV*invV

    res = Measurement()
    res["Magnetization"] = M
    res["|Magnetization|"] = abs(M)
    res["Magnetization^2"] = M2
    res["Magnetization^4"] = M4
    res["Energy"] = E
    res["Energy^2"] = E2

    return res
end

@doc """
    improved_estimator(model::Potts, T::Real, Js::AbstractArray, sw::SWInfo)

Returns the following observables as `Dict{String, Any}` using cluster information `sw`

# Observables
- `"Energy"`
- `"Energy^2"`
- `"Magnetization"`
- `"|Magnetization|"`
- `"|Magnetization|^2"`
- `"|Magnetization|^4"`
"""
function improved_estimator(model::Potts, T::Real, Js::AbstractArray, sw::SWInfo)
    nsites = numsites(model)
    Q = model.Q
    nc = numclusters(sw)
    invV = 1.0/nsites

    # magnetization
    I2 = (Q-1)/(Q*Q)
    I4 = (Q-1)*((Q-1)^3+1)/(Q^5)
    M = 0.0
    M2 = 0.0
    M4 = 0.0
    s2 = -1.0/Q
    s1 = 1.0+s2
    for (m,s) in zip(sw.clustersize, sw.clusterspin)
        M += (m*invV) * ifelse(s==1, s1, s2)
        m2 = (m*invV)^2
        M4 += I4*m2*m2 + 6*I2*M2*m2
        M2 += I2*m2
    end

    # energy
    aJ = abs.(Js)
    mbeta = -1.0/T
    ns = sw.activated_bonds
    As = -aJ ./ expm1.(mbeta.*aJ)
    Ans = ns.*As
    E = 0.0
    E2 = 0.0
    for b in 1:numbondtypes(model)
        E2 += aJ[b]*Ans[b]
        E2 += Ans[b] * As[b]*(ns[b]-1)
        E2 += 2.0*Ans[b]*E
        E += Ans[b]
    end

    E *= -invV
    E2 *= invV*invV

    res = Measurement()
    res["Magnetization"] = M
    res["|Magnetization|"] = abs(M)
    res["Magnetization^2"] = M2
    res["Magnetization^4"] = M4
    res["Energy"] = E
    res["Energy^2"] = E2

    return res
end

@doc """
    improved_estimator(model::QuantumXXZ, T::Real, Js::AbstractArray, uf::UnionFind)

Returns the following observables as `Dict{String, Any}` using loop information `uf`

# Observables
- `"Sign"`
    - Sign of the weight function
- `"Sign * Energy"`
    - Energy per spin (site)
- `"Sign * Energy^2"`
- `"Sign * Magnetization"`
    - Total magnetization (Sz) per spin (site)
- `"Sign * |Magnetization|"`
- `"Sign * Magnetization^2"`
- `"Sign * Magnetization^4"`
"""
function improved_estimator(model::QuantumXXZ, T::Real, Jzs::AbstractArray, Jxys::AbstractArray, Gs::AbstractArray, uf::UnionFind)
    S2 = model.S2
    nsites = numsites(model)
    nspins = nsites*S2
    nbonds = numbonds(model)
    nc = numclusters(uf)

    ms = zeros(nc)
    for s in 1:nspins
        i = clusterid(uf, s)
        ms[i] += model.spins[s]
    end
    ms .*= 0.5/nsites

    E0 = 0.0
    for st in 1:numsitetypes(model)
        E0 += 0.5 * numsites(model, st) * Gs[st] * S2
    end
    for bt in 1:numbondtypes(model)
        nb = numbonds(model,bt)*S2*S2
        z = nb*Jzs[bt]
        x = nb*abs(Jxys[bt])
        if z > x
            ## AntiFerroIsing like
            E0 += 0.25z
        elseif z < -x
            ## FerroIsing like
            E0 -= 0.25z
        else
            ## XY like
            E0 += 0.25x
        end
    end
    nops = length(model.ops)
    E = E0-nops*T
    E2 = nops*(nops-1)*T^2 - 2*E0*T*nops + E0^2
    E /= nsites
    E2 /= nsites^2

    M = 0.0
    M2 = 0.0
    M4 = 0.0
    for m in ms
        M += m
        m2 = m*m
        M4 += m2*m2 + 6M2*m2
        M2 += m2
    end

    sgn = 1.0
    for op in model.ops
        if !op.isdiagonal
            if op.let_type == LET_Vertex || op.let_type == LET_Cross
                subbond = op.space
                b,_,_ = subbond2bond(subbond,S2)
                bt = bondtype(model, b)
                sgn *= ifelse(Jxys[bt] > 0.0, -1.0, 1.0)
            end
        end
    end

    res = Measurement()
    res["Sign * Magnetization"] = sgn*M
    res["Sign * |Magnetization|"] = sgn*abs(M)
    res["Sign * Magnetization^2"] = sgn*M2
    res["Sign * Magnetization^4"] = sgn*M4
    res["Sign * Energy"] = sgn*E
    res["Sign * Energy^2"] = sgn*E2
    res["Sign"] = sgn

    return res
end


text:
"""
    r_x(r, vl, vr)

Transfer coordinate r -> x from standard to real elements

"""
r_x(r, vl, vr) = ((1.0 - r) / 2.0) * vl + ((1.0 + r) / 2.0) * vr

"""
Transfer coordinate r -> x from standard to real elements

Triangle elements:

    rs_xy(r, s, v1, v2, v3)
    rs_xy(v, v1, v2, v3)

"""
function rs_xy(
    r,
    s,
    v1::AbstractVector{T},
    v2::AbstractVector{T},
    v3::AbstractVector{T},
) where {T<:Real}
    return @. -(r + s) / 2 * v1 + (r + 1) / 2 * v2 + (s + 1) / 2 * v3
end

function rs_xy(
    v::AbstractVector{T},
    v1::AbstractVector{T},
    v2::AbstractVector{T},
    v3::AbstractVector{T},
) where {T<:Real}
    r, s = v
    return rs_xy(r, s, v1, v2, v3)
end

function rs_xy(
    r,
    s,
    v1::AbstractVector{T},
    v2::AbstractVector{T},
    v3::AbstractVector{T},
    v4::AbstractVector{T},
) where {T<:Real}
    return @. (r - 1.0) * (s - 1.0) / 4.0 * v1 +
        (r + 1.0) * (1.0 - s) / 4.0 * v2 +
        (r + 1.0) * (s + 1.0) / 4.0 * v3 +
        (1.0 - r) * (s + 1.0) / 4.0 * v4
end

"""
    xy_rs(x, y)
    xy_rs(coords)

Transfer coordinates (x, y) -> (r,s) from equilateral to right triangle

"""
function xy_rs(x::T, y::T) where {T<:Real}
    L1 = (sqrt(3.0) * y + 1.0) / 3.0
    L2 = (-3.0 * x - sqrt(3.0) * y + 2.0) / 6.0
    L3 = (3.0 * x - sqrt(3.0) * y + 2.0) / 6.0
    r = -L2 + L3 - L1
    s = -L2 - L3 + L1

    return r, s
end

function xy_rs(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:Real}
    Np = length(x)
    r = zeros(Np)
    s = zeros(Np)

    for n = 1:Np
        r[n], s[n] = xy_rs(x[n], y[n])
    end

    return r, s
end

function xy_rs(coords::AbstractMatrix{T}) where {T<:Real}
    return xy_rs(coords[:, 1], coords[:, 2])
end

"""
    rs_ab(r, s)
    rs_ab(coords)

Transfer coordinates (r,s) -> (a,b) in a triangle

"""
function rs_ab(r::T, s::T) where {T<:Real}
    a = ifelse(s != 1.0, 2.0 * (1.0 + r) / (1.0 - s) - 1.0, -1.0)
    b = 1.0 * s

    return a, b
end

function rs_ab(r::AbstractVector{T}, s::AbstractVector{T}) where {T<:Real}
    a = zero(r)
    b = zero(s)

    for n in eachindex(a)
        a[n], b[n] = rs_ab(r[n], s[n])
    end

    return a, b
end

function rs_ab(coords::AbstractMatrix{T}) where {T<:Real}
    return rs_ab(coords[:, 1], coords[:, 2])
end


text:
module CUDAdrv

using CUDAapi

using CEnum

using Printf

## source code includes

# essential functionality
include("pointer.jl")
const CUdeviceptr = CuPtr{Cvoid}

# low-level wrappers
const libcuda = Sys.iswindows() ? :nvcuda : :libcuda
include("libcuda_common.jl")
include("error.jl")
include("libcuda.jl")
include("libcuda_aliases.jl")

include("util.jl")

# high-level wrappers
include("version.jl")
include("devices.jl")
include("context.jl")
include(joinpath("context", "primary.jl"))
include("stream.jl")
include("memory.jl")
include("module.jl")
include("events.jl")
include("execution.jl")
include("profile.jl")
include("occupancy.jl")

include("deprecated.jl")

## initialization

const __initialized__ = Ref(false)
functional() = __initialized__[]

function __init__()
    precompiling = ccall(:jl_generating_output, Cint, ()) != 0
    silent = parse(Bool, get(ENV, "JULIA_CUDA_SILENT", "false")) || precompiling
    verbose = parse(Bool, get(ENV, "JULIA_CUDA_VERBOSE", "false"))

    try
        if haskey(ENV, "_") && basename(ENV["_"]) == "rr"
            error("Running under rr, which is incompatible with CUDA")
        end

        cuInit(0)

        if version() <= v"9"
            silent || @warn "CUDAdrv.jl only supports NVIDIA drivers for CUDA 9.0 or higher (yours is for CUDA $(version()))"
        end

        __initialized__[] = true
    catch ex
        # don't actually fail to keep the package loadable
        if !silent
            if verbose
                @error "CUDAdrv.jl failed to initialize" exception=(ex, catch_backtrace())
            else
                @info "CUDAdrv.jl failed to initialize, GPU functionality unavailable (set JULIA_CUDA_SILENT or JULIA_CUDA_VERBOSE to silence or expand this message)"
            end
        end
    end
end

end


text:
# Imports
using Random, Statistics, Plots, Printf
Random.seed!(430);

# Part (a)
N = 10000;
Xₙ = randn(2, 1);
μₓ = [0.5; 0.3];
Cₓ = [1 0.6; 0.6 sqrt(2)];

# 2-dimensional Gaussian vector
X = Cₓ * Xₙ + μₓ;
println("Answer of part (a)")
println("2-dimensional Gaussian vector:")
println(X)

# Part (b)
# Generate 10000 vector with related parameters
function generate(matrix, μₓ, n)
    samples = zeros(N, 2)
    for i in 1:n
        X = randn(2, 1);
        X = sqrt(matrix) * X + μₓ;
        samples[i, :] = X;
    end
    return samples
end

# Plot generated samples
samples = generate(Cₓ, μₓ, N)
X₁, X₂ = samples[:, 1], samples[:, 2];
display(plot(X₁, X₂, seriestype = :scatter, title = "Scatter Plot"))

# Part (c)
#∑ = zeros(2,2)
∑ = []
for i in 1:N
    #∑ += samples[i, :]*samples[i, :]'
    push!(∑, samples[i, :]*samples[i, :]')
end
#Sₓ = ∑/N
Sₓ = sum(∑) / N;
println("\nAnswer of part (c) for $N samples: ")
println(Sₓ)

# When we increase "N", convergence will be [1.25 0.75; 0.75 1.504] matrix
println("Convergence limit matrix:")
convergence = Cₓ + μₓ*μₓ'


text:
version = v"14.0.2"

include("../common.jl")

build_tarballs(ARGS, configure_build(ARGS, version; experimental_platforms=true)...;
               preferred_gcc_version=v"7", preferred_llvm_version=v"12", julia_compat="1.9")

max_stars_repo_name: staticfloat/BinaryBuilderRepos

max_issues_repo_name: staticfloat/BinaryBuilderRepos

max_forks_repo_name: staticfloat/BinaryBuilderRepos


text:
using SnoopCompile, Test

cd(@__DIR__)
@testset "bot" begin
    @testset "precompileInclude" begin
        @testset "precompile_pather" begin
            precompilePath, precompileFolder = precompile_pather("TestPackage")
            @test  precompilePath == "\"../deps/SnoopCompile/precompile/precompile_TestPackage.jl\""
            @test precompileFolder == "$(pwd())/deps/SnoopCompile/precompile/"
        end

        @testset "precompile_activator" begin
            Base.write("activated.jl", """
            include("../deps/SnoopCompile/precompile/precompile_TestPackage.jl")
            _precompile_()
            """)

            precompilePath, precompileFolder = precompile_pather("TestPackage")
            @test_logs (:info,"precompile is already activated") precompile_activator("activated.jl", precompilePath)
        end

        @testset "precompile_deactivator" begin
            Base.write("deactivated.jl", """
            # include("../deps/SnoopCompile/precompile/precompile_TestPackage.jl")
            # _precompile_()
            """)

            precompilePath, precompileFolder = precompile_pather("TestPackage")
            @test_logs (:info, "precompile is already deactivated") precompile_deactivator("deactivated.jl", precompilePath)
        end
    end

    # Fails. # snoopi_bot pwd() differs from what it is cd to
    # instead test the code in MatLang actions: https://github.com/juliamatlab/MatLang/actions?query=workflow%3ASnoopCompile
    #=
    @testset "snoopi_bot" begin
        using Pkg; Pkg.develop("MatLang")

        examplePath = Base.read(`cmd /c julia -e 'import MatLang; print(pathof(MatLang))'`, String)
        cd(dirname(dirname(examplePath)))

        @snoopi_bot "MatLang" begin

          using MatLang
          examplePath = joinpath(dirname(dirname(pathof(MatLang))), "examples")
          include(joinpath(examplePath,"Language_Fundamentals", "usage_Matrices_and_Arrays.jl"))
          include(joinpath(examplePath,"Language_Fundamentals", "Data_Types", "usage_Numeric_Types.jl"))
        end
    end

    @testset "snoopi_bench" begin
        using Pkg; Pkg.develop("MatLang")

        examplePath = Base.read(`cmd /c julia -e 'import MatLang; print(pathof(MatLang))'`, String)
        cd(dirname(dirname(examplePath)))

        @snoopi_bench "MatLang" begin
            using MatLang
            examplePath = joinpath(dirname(dirname(pathof(MatLang))), "examples")
            include(joinpath(examplePath,"Language_Fundamentals", "usage_Matrices_and_Arrays.jl"))
            include(joinpath(examplePath,"Language_Fundamentals", "Data_Types", "usage_Numeric_Types.jl"))
        end
    end
    =#

end


text:
@parameters t
D = Differential(t)

function RLC()
    r = Component(:R)
    l = Component(:I)
    c = Component(:C)
    kvl = EqualEffort(name = :kvl)

    bg = BondGraph()
    add_node!(bg, [c, l, kvl, r])

    connect!(bg, r, kvl)
    connect!(bg, l, kvl)
    connect!(bg, c, kvl)
    return bg
end

@testset "Equations" begin
    c = Component(:C)
    @parameters C
    @variables E[1](t) F[1](t) q(t)
    cr = [
        0 ~ q / C - E[1],
        D(q) ~ F[1]
    ]

    @test isequal(equations(c), cr)
    @test isequal(constitutive_relations(c), cr)

    j = EqualEffort()
    @test isequal(equations(j), Equation[])
end

@testset "Parameters" begin
    tf = Component(:TF)
    @parameters n
    @test iszero(parameters(tf) - [n])

    Ce = Component(:Ce)
    @parameters k R T
    @test iszero(parameters(Ce) - [k, R, T])

    Re = Component(:Re)
    @parameters r R T
    @test iszero(parameters(Re) - [r, R, T])

    bg = RLC()
    @parameters C L R
    @test iszero(parameters(bg) - [C, L, R])
end

@testset "State variables" begin
    r = Component(:R)
    @test isempty(states(r))

    @variables q(t)
    c = Component(:C)
    @test isequal(states(c), [q])

    ce = Component(:ce)
    @test isequal(states(ce), [q])

    bg = RLC()
    @variables q(t) p(t)
    @test iszero(states(bg) - [q, p])
end

@testset "Constitutive relations" begin
    eqE = EqualEffort()
    eqF = EqualFlow()
    @test constitutive_relations(eqE) == Equation[]
    @test constitutive_relations(eqF) == Equation[]

    bg = RLC()
    cr_bg = constitutive_relations(bg)
    sys = ODESystem(bg)

    C, L, R = sys.ps
    q, p = sys.states
    cr1 = D(q) ~ -(q / C) / R + (-p) / L
    cr2 = D(p) ~ q / C

    @test isequal(cr_bg[1].lhs,cr1.lhs)
    @test isequal(simplify(cr_bg[1].rhs-cr1.rhs), 0)
    @test isequal(cr_bg[2], cr2)

    cr_bgn = constitutive_relations(BondGraphNode(bg))
    @test isequal(cr_bgn[1].lhs,cr1.lhs)
    @test isequal(simplify(cr_bgn[1].rhs-cr1.rhs), 0)
    @test isequal(cr_bgn[2], cr2)
end

@testset "0-junction equations" begin
    model = BondGraph(:RC)
    C = Component(:C)
    R = Component(:R)
    zero_law = EqualEffort()

    add_node!(model, [R, C, zero_law])
    connect!(model, R, zero_law)
    connect!(model, zero_law, C)

    @test numports(zero_law) == 2

    @variables E[1:2](t) F[1:2](t)
    @test isequal(constitutive_relations(zero_law), [
        0 ~ F[1] + F[2],
        0 ~ E[1] - E[2]
    ])
end

@testset "1-junction equations" begin
    c1 = Component(:C, :C1)
    c2 = Component(:R, :R1)
    c3 = Component(:I, :I1)
    j = EqualFlow()

    bg = BondGraph()
    add_node!(bg, [c1, c2, c3, j])
    connect!(bg, c1, j)
    connect!(bg, j, c2)
    connect!(bg, j, c3)

    @test numports(j) == 3
    @test length(j.weights) == 3
    @test j.weights == [1, -1, -1]

    @variables E[1:3](t) F[1:3](t)
    @test isequal(constitutive_relations(j), [
        0 ~ E[1] - E[2] - E[3],
        0 ~ F[1] + F[2],
        0 ~ F[1] + F[3],
    ])
end

@testset "RC circuit" begin
    r = Component(:R)
    c = Component(:C)
    bg = BondGraph(:RC)
    add_node!(bg, [c, r])
    connect!(bg, r, c)

    sys = ODESystem(bg)
    eqs = constitutive_relations(bg)
    @test length(eqs) == 1

    (C, R) = sys.ps
    x = sys.states[1]
    e1 = eqs[1]
    e2 = D(x) ~ -x / C / R

    @test isequal(e1.lhs, e2.lhs)
    @test isequal(expand(e1.rhs), e2.rhs)
end

@testset "RL circuit" begin
    r = Component(:R)
    l = Component(:I)
    bg = BondGraph(:RL)
    add_node!(bg, [r, l])
    connect!(bg, l, r)

    eqs = constitutive_relations(bg)
    sys = ODESystem(bg)
    x = sys.states[1]
    (R,L) = sys.ps
    @test eqs == [D(x) ~ -R*x/L]
end

@testset "RLC circuit" begin
    bg = RLC()
    eqs = constitutive_relations(bg)
    @test length(eqs) == 2

    sys = ODESystem(bg)
    (C, L, R) = sys.ps
    (qC, pL) = sys.states
    e1 = D(qC) ~ -pL / L + (-qC / C / R)
    e2 = D(pL) ~ qC / C

    @test isequal(simplify(eqs[1].rhs-e1.rhs), 0)
    @test isequal(eqs[2].rhs, e2.rhs)
end

@testset "Chemical reaction A ⇌ B" begin
    A = Component(:ce, :A)
    B = Component(:ce, :B)
    re = Component(:re, :r)
    bg = BondGraph()

    add_node!(bg, [A, B, re])
    connect!(bg, A, re; dstportindex = 1)
    connect!(bg, re, B; srcportindex = 2)
    sys = ODESystem(bg)
    eqs = constitutive_relations(bg)

    (xA, xB) = sys.states
    (KA, KB, r) = sys.ps
    e1 = D(xA) ~ -r * (KA * xA - KB * xB)
    e2 = D(xB) ~ r * (KA * xA - KB * xB)

    @test isequal(eqs[1].rhs, e1.rhs)
    @test isequal(eqs[2].rhs, e2.rhs)
end

@testset "Chemical reaction A ⇌ B + C, C ⇌ D" begin
    C_A = Component(:ce, :A)
    C_B = Component(:ce, :B)
    C_C = Component(:ce, :C)
    C_D = Component(:ce, :D)
    re1 = Component(:re, :r1)
    re2 = Component(:re, :r2)
    common_C = EqualEffort()
    BC = EqualFlow()

    bg = BondGraph()
    add_node!(bg, [C_A, C_B, C_C, C_D, re1, re2, common_C, BC])
    connect!(bg, C_A, re1; dstportindex = 1)
    connect!(bg, re1, BC; srcportindex = 2)
    connect!(bg, BC, C_B)
    connect!(bg, BC, common_C)
    connect!(bg, common_C, C_C)
    connect!(bg, common_C, re2; dstportindex = 1)
    connect!(bg, re2, C_D; srcportindex = 2)

    sys = ODESystem(bg)
    eqs = constitutive_relations(bg)

    (xA, xB, xC, xD) = sys.states
    (KA, KB, KC, KD, r1, r2) = sys.ps
    e1 = D(xA) ~ -r1 * (KA * xA - KB * xB * KC * xC)
    e2 = D(xB) ~ r1 * (KA * xA - KB * xB * KC * xC)
    e3 = D(xC) ~ r1 * (KA * xA - KB * xB * KC * xC) - r2 * (KC * xC - KD * xD)
    e4 = D(xD) ~ r2 * (KC * xC - KD * xD)

    # equations are not simplifying with exp/log rules
    @test isequal(eqs[1].rhs, e1.rhs)
    @test isequal(eqs[2].rhs, e2.rhs)
    @test isequal(eqs[3].rhs, e3.rhs)
    @test isequal(eqs[4].rhs, e4.rhs)
end

max_stars_repo_name: mic-pan/BondGraphs

max_issues_repo_name: mic-pan/BondGraphs

max_forks_repo_name: mic-pan/BondGraphs


text:
# Most of the codes here are copied and modified from PkgServer.jl

const uuid_re = raw"[0-9a-z]{8}-[0-9a-z]{4}-[0-9a-z]{4}-[0-9a-z]{4}-[0-9a-z]{12}(?-i)"
const hash_re = raw"[0-9a-f]{40}"
const registry_re = Regex("^/registry/($uuid_re)/($hash_re)\$")
const resource_re = Regex("""
    ^/registry/$uuid_re/($hash_re)\$
  | ^/package/$uuid_re/($hash_re)\$
  | ^/artifact/($hash_re)\$
""", "x")

function gen_curl_cmd(url, download_path)
    options = [
        "-sSfL",
        "-H", "Content-Type: application/tar+gzip",
        "-H", "Content-Encoding: identity",
    ]
    if haskey(ENV, "BIND_ADDRESS")
        addr = ENV["BIND_ADDRESS"]
        pushfirst!(options, "--interface", addr)
    end

    return `curl $options $url -o $download_path`
end

"""
    query_latest_hash(registry, server)

Interrogate a storage server for a list of registries, match the response against the
registries we are paying attention to, and return the matching hash.

If registry is not avilable in server, then return `nothing`.
"""
function query_latest_hash(registry::RegistryMeta, server::AbstractString)
    resource = "/registries"

    try
        response = timeout_call(30) do
            HTTP.get(server * resource)
        end
        @debug "succeed to fetch resource" resource=server*resource status=response.status
        return get_hash(IOBuffer(response.body), registry.uuid)
    catch err
        @warn "failed to fetch resource" error=err resource=server*resource
        return nothing
    end
end

get_hash(contents::AbstractString, uuid::AbstractString) = get_hash(IOBuffer(contents), uuid)
function get_hash(io::IO, uuid::AbstractString)
    for line in eachline(io)
        m = match(registry_re, strip(line))
        if m !== nothing
            matched_uuid, matched_hash = m.captures
            matched_uuid == uuid && return matched_hash
        end
    end
    return nothing
end

"""
    query_latest_hash(registry::RegistryMeta, upstreams::AbstractVector)

Query `upstreams` for the latest registry hash. Return `nothing` if given registry isn't available
in all upstreams.
"""
function query_latest_hash(registry::RegistryMeta, upstreams::AbstractVector{<:AbstractString})
    upstreams = normalize_upstream.(upstreams)

    # collect current registry hashes from servers
    uuid = registry.uuid
    hash_info = Dict{String, Vector{String}}() # Dict(hashA => [serverA, serverB], ...)
    servers = String[] # [serverA, serverB]
    @sync for server in upstreams
        @async begin
            hash = query_latest_hash(registry, server)
            if !isnothing(hash)
                push!(get!(hash_info, hash, String[]), server)
                push!(servers, server)
            end
        end
    end

    # for each hash check what other servers know about it
    if isempty(hash_info)
        # reach here if none of the upstreams contains the registry we want to mirror
        @error "failed to find available registry" registry=registry.name uuid=registry.uuid upstreams=join(upstreams, ", ")
        return nothing
    end

    # a hash might be known to many upstreams
    for (hash, hash_servers) in hash_info
        @sync for server in servers
            server in hash_servers && continue
            @async url_exists("$server/registry/$uuid/$hash") && push!(hash_servers, server)
        end
    end

    # Ideally, there is an upstream server that knows all hashes, and we set hash in that server
    # as the latest hash.
    # In practice, we set the first non-malicious hash known to fewest servers as the latest hash.
    hashes = sort!(collect(keys(hash_info)))
    sort!(hashes, by = hash -> length(hash_info[hash]))
    hashes[findfirst(x->verify_registry_hash(registry.source_url, x), hashes)]
end

"""
    normalize_upstream(upstream::AbstractString) -> String

Normalize server url by adding necessary prefix and removing trailing `/`
"""
function normalize_upstream(upstream::AbstractString; prefix = "https")
    startswith(upstream, r"\w+://") || (upstream = "$prefix://$upstream")
    return String(rstrip(upstream, '/'))
end

"""
    url_exists(url; timeout=30_000)

Send a `HEAD` request to the specified URL, returns `true` if the response is HTTP 200.

Set `timeout=0` millseconds to disable timeout.
"""
function url_exists(url::AbstractString; timeout::Integer=30_000, throw_warnings=true)
    startswith(url, r"https?://") || throw(ArgumentError("invalid url $url, should be HTTP(S) protocol."))

    f() = HTTP.request("HEAD", url, status_exception=false)
    try
        response = timeout==0 ? f() : timeout_call(f, timeout//1000)
        throw_warnings && @debug "succeed to send HEAD request" response=response url=url
        # FIXME: unsure why could response be nothing
        return !isnothing(response) && response.status == 200
    catch err
        # FIXME: here we temporarily disable unlock warnings as the tasks are interrupted manually
        err isa ErrorException && occursin("unlock from wrong thread", err.msg) && return false

        throw_warnings && @warn "failed to send HEAD request" error=err url=url
        return false
    end
end

"""
    verify_registry_hash(source_url, hash)

Verify that the origin git repository knows about the given registry tree hash.
"""
function verify_registry_hash(source_url::AbstractString, hash::AbstractString)
    url = Pkg.Operations.get_archive_url_for_version(source_url, hash)
    return url === nothing || url_exists(url)
end

"""
    write_atomic(f::Function, path::String)
Performs an atomic filesystem write by writing out to a file on the same
filesystem as the given `path`, then `move()`'ing the file to its eventual
destination.  Requires write access to the file and the containing folder.
Currently stages changes at "<path>.tmp.<randstring>".  If the return value
of `f()` is `false` or an exception is raised, the write will be aborted.
"""
function write_atomic(f::Function, path::String)
    isdir(dirname(path)) || mkpath(dirname(path))
    temp_file = path * ".tmp." * randstring()
    try
        retval = f(temp_file)
        if retval !== false
            mv(temp_file, path; force=true)
        end
    catch e
        rethrow(e)
    finally
        rm(temp_file; force=true)
    end
end

function download_and_verify(
        server::String,
        resource::String,
        tarball::String;
        http_parameters::Dict{Symbol, Any} = Dict{Symbol, Any}(),
)
    http_parameters = merge(default_http_parameters, http_parameters)
    timeout = http_parameters[:timeout]
    delete!(http_parameters, :timeout)

    hash = let m = match(resource_re, resource)
        m !== nothing ? first(filter(x->!isnothing(x), m.captures)) : nothing
    end

    if isnothing(hash)
        @warn "bad resource: valid hash not found" resource=server*resource tarball=tarball
        return false
    end

    # Verifying hash requires a lot of IO reads. A faster way is to only check if the file
    # exists. This is okay if we can make sure the file isn't created when download/creation
    # of tarball fails
    isfile(tarball) && return true

    write_atomic(tarball) do temp_file
        try
            # HTTP.get is problematic and occasionally hangs, use `curl` as a replacement
            # Set Content-Encoding because we want exactly the same content as the upstream serves
            # without any processing
            run(gen_curl_cmd(server*resource, temp_file))
        catch err
            @warn "failed to fetch resource" error=err resource=server*resource
            return false
        end

        # If we're given a hash, then check tarball git hash
        if hash !== nothing
            tree_hash = open(temp_file, "r") do io
                Tar.tree_hash(decompress(io), algorithm="git-sha1")
            end

            # Raise warnings about resource hash mismatches
            if hash != tree_hash
                # julia has changed how hash is calculated, which affects empty directory
                # we specially check whethere hash mismatch is caused by this, and if so,
                # skip the warning and make a copy
                ytree_hash = open(temp_file, "r") do io
                    Tar.tree_hash(decompress(io), algorithm="git-sha1", skip_empty=true)
                end

                if hash == ytree_hash
                    yskip_tarball = joinpath(dirname(tarball), ytree_hash)
                    if !isfile(yskip_tarball)
                        temp_yskip_tarball = yskip_tarball * ".tmp." * randstring()
                        cp(temp_file, temp_yskip_tarball)
                        mv(temp_yskip_tarball, yskip_tarball)
                    end
                    return true
                end

                @warn "resource hash mismatch" resource=server*resource reference_hash=hash actual_hash=tree_hash
                return false
            end
        end

        return true
    end

    return isfile(tarball)
end

function download_and_verify(
        servers::AbstractVector{String},
        resource::String,
        tarball::String;
        throw_warnings = true,
        kwargs...
)
    isfile(tarball) && return true

    race_lock = ReentrantLock()
    task_pool = []

    try
        if length(servers) == 1
            download_and_verify(servers[1], resource, tarball; kwargs...)
        else
            for server in servers
                task = @async begin
                    if url_exists(server*resource; timeout=30_000, throw_warnings=throw_warnings)
                        # the first that hits here start downloading
                        if trylock(race_lock)
                            retval = download_and_verify(server, resource, tarball; kwargs...)
                            unlock(race_lock)
                            retval
                        end
                    else
                        throw_warnings && @warn "resource does not exist" resource=server*resource
                    end
                end
                push!(task_pool, task)
            end

            try
                timeout_call(default_http_parameters[:timeout]) do
                    while true
                        sleep(0.1)

                        if any(t->istaskdone(t) && t.result === true, task_pool)
                            interrupt_task(task_pool)
                            return true
                        end

                        if all(istaskdone, task_pool)
                            @warn "fail to download resource" resource=resource tarball=tarball upstreams=join(servers, ", ")
                            return false
                        end
                    end
                end
            catch err
                @warn err
                return false
            end
        end
    catch e
        throw_warnings && @warn "fail to download resource" err=e resource=resource tarball=tarball upstreams=join(servers, ", ")
        interrupt_task(task_pool)
        return false
    end

    if isfile(tarball)
        return true
    else
        throw_warnings && @warn "unable to download resource from all given upstreams" resource=resource tarball=tarball upstreams=join(servers, ", ")
        return false
    end
end

function interrupt_task(task_pool)
    isempty(task_pool) && return nothing
    for t in task_pool
        @eval :(Base.throwto(t, InterruptException()))
    end
end


text:
# Autogenerated wrapper script for RTKLib_jll for x86_64-linux-gnu
export convbin, librtklib, pos2kml, rtkrcv, str2str, urllist

JLLWrappers.@generate_wrapper_header("RTKLib")
JLLWrappers.@declare_executable_product(convbin)
JLLWrappers.@declare_library_product(librtklib, "librtklib.so")
JLLWrappers.@declare_executable_product(pos2kml)
JLLWrappers.@declare_executable_product(rtkrcv)
JLLWrappers.@declare_executable_product(str2str)
JLLWrappers.@declare_file_product(urllist)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        convbin,
        "bin/convbin",
    )

    JLLWrappers.@init_library_product(
        librtklib,
        "lib/librtklib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        pos2kml,
        "bin/pos2kml",
    )

    JLLWrappers.@init_executable_product(
        rtkrcv,
        "bin/rtkrcv",
    )

    JLLWrappers.@init_executable_product(
        str2str,
        "bin/str2str",
    )

    JLLWrappers.@init_file_product(
        urllist,
        "share/RTKLib/URL_LIST.txt",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ProblemParser
using Utils

using Accessors
using CompositeStructs
using Underscores

const file = find_input(@__FILE__)
const input = @_ file |> slurp |> parse(Lines(FirstRest(Split(), Apply(Symbol), Convert())), __)

@Base.kwdef struct Position
    position::Int = 0
    depth::Int = 0
end

function Base.:+(pos::Position, directions::Tuple{Symbol, Int})
    command, amount = directions
    if command == :forward
        @set pos.position += amount
    elseif command == :down
        @set pos.depth += amount
    elseif command == :up
        @set pos.depth -= amount
    end
end

@composite @Base.kwdef struct AimingPosition
    Position...
    aim::Int = 0
end

function Base.:+(pos::AimingPosition, directions::Tuple{Symbol, Int})
    command, amount = directions
    if command == :forward
        AimingPosition(pos.position + amount, pos.depth + pos.aim * amount, pos.aim)
    elseif command == :down
        @set pos.aim += amount
    elseif command == :up
        @set pos.aim -= amount
    end
end

function solve(input, type)
    final = sum(input; init=type())
    final.position * final.depth
end

part1(input) = solve(input, Position)
assertequal(part1(input), 2039912)

part2(input) = solve(input, AimingPosition)
assertequal(part2(input), 1942068080)

max_stars_repo_name: CmdQ/AoC


text:
using KrylovMethods
using Base.Test
using LinearOperators

println("=== Testing gmres for real matrix === ")
A  = sprandn(100,100,.1) + 10*speye(100)
n  = size(A,2)
D  = diag(A)
M2 = x -> D.\x
rhs = randn(100)
tol = 1e-6;

# test printing and behaviour for early stopping
xtt = gmres(A,rhs ,3,tol=1e-10,maxIter=3,out=2,storeInterm=true)
@test xtt[2]==-1

# test behaviour for zero rhs
xtt = gmres(A,0*rhs,5,tol=tol,maxIter=10,out=2)
@test xtt[2]==-9
@test all(xtt[1].==0)
@test length(xtt[1])==100
@test eltype(xtt[1])==eltype(rhs)

x1 = gmres(A,rhs ,5,tol=tol,maxIter=100,out=1)
x3 = gmres(A,rhs,5,tol=tol,maxIter=100,x=randn(size(rhs)))
x4 = gmres(A,rhs,5,tol=tol,maxIter=100,M=M2)

@test norm(A*x1[1]-rhs)/norm(rhs) < tol
@test norm(A*x3[1]-rhs)/norm(rhs) < tol
# @test norm(A*x4[1]-rhs)/norm(rhs) < tol
@test norm(x3[1]-x1[1])/norm(x1[1]) < 1e-5

println("=== Testing gmres for complex matrix === ")
A  = sprandn(100,100,.1) + 10*speye(100) + im*(sprandn(100,100,.1) + 10*speye(100) )
D  = diag(A)
M3 = x -> D.\x
rhs = complex(randn(100))

# test behaviour for zero rhs
xtt = gmres(A,0*rhs,5,tol=tol,maxIter=10,out=2)
@test xtt[2]==-9
@test all(xtt[1].==0)
@test length(xtt[1])==100
@test eltype(xtt[1])==eltype(rhs)

x1 = gmres(A,rhs ,5,tol=tol,maxIter=100)
x3 = gmres(A,rhs,5,tol=tol,maxIter=100,x=randn(size(rhs)))
x4 = gmres(A,rhs,5,tol=tol,maxIter=100,M=M3)

@test norm(A*x1[1]-rhs)/norm(rhs) < tol
@test norm(A*x3[1]-rhs)/norm(rhs) < tol
@test norm(x3[1]-x1[1])/norm(x1[1]) < 1e-5

println("=== GMRES: All tests passed. ===")


text:
try
    using FoldsPreviewTests
    true
catch
    false
end || begin
    push!(LOAD_PATH, @__DIR__)
    using FoldsPreviewTests
end


text:
export draw

using PyCall
using Conda
Conda.add("python-graphviz", channel="conda-forge")
Conda.add("pydotplus", channel="conda-forge")

const PydotPlus = PyNULL()

function __init__()
    copy!(PydotPlus, pyimport_conda("pydotplus", "pydotplus"))
end

function draw(dot_data)
    graph = PydotPlus.graph_from_dot_data(dot_data)
    graph.progs = Dict("dot" => "$(Conda.BINDIR)/dot")
    display("image/png", Vector{UInt8}(graph.create_png()))
end


text:
module DateUtils

# ARTICLES
#  - http://old.sevierriver.org/et/t110.html
#
#

export date2julian,
    julian_year;

@doc """
Date to Julian Day converter
"""
date2julian(date::Date) = Dates.datetime2julian(convert(DateTime, date))
date2julian(date::Date, since::Date) = Int(date2julian(date) - date2julian(since))

@doc """
Julian Day Of Year
Date to Julian Day converter
"""
julian_year(date::Date) = date2julian(date, Date(Dates.year(date)))

end # module


text:
@testset "Postprocessing" begin
    @test LatexSVG._adjust_web_svg_dim("16.180339pt") == "1.618034rem"

    @test begin
        svg_test = """
        <?xml version="1.0" encoding="utf-8"?>
        <svg width="21.923673pt" height="8.966376pt" viewBox="67.695921 70.081196 21.923673 8.966376">
        <defs/>
        </svg>
        """
        svg_adjusted = LatexSVG._adjust_web_svg_display(svg_test)
        svg_adjusted == """
        <?xml version="1.0" encoding="UTF-8"?>
        <svg height="0.896638rem" viewBox="67.695921 70.081196 21.923673 8.966376" class="latexsvg-display" style="display: block; margin: auto; max-width: 80%">
        <defs/>
        </svg>
        """
    end

    @test begin
        svg_test = """
        <?xml version="1.0" encoding="utf-8"?>
        <svg width="21.923673pt" height="8.966376pt" viewBox="67.695921 70.081196 21.923673 8.966376">
        <defs/>
        </svg>
        """
        svg_adjusted = LatexSVG._adjust_web_svg_inline(svg_test)
        svg_adjusted == """
        <?xml version="1.0" encoding="UTF-8"?>
        <svg height="0.896638rem" viewBox="67.695921 70.081196 21.923673 8.966376" class="latexsvg-inline" style="vertical-align: middle">
        <defs/>
        </svg>
        """
    end
end


text:
using Zygote, Test

function tasks1(x)
  ch = Channel(Inf)
  put!(ch, x^2) + take!(ch)
end

@test gradient(tasks1, 5) == (20,)

function tasks2(x)
  ch = Channel(0)
  t = @async put!(ch, x^2)
  y = take!(ch)
  wait(t)
  return y
end

@test gradient(tasks2, 5) == (10,)

function tasks3(x)
  ch = Channel(0)
  @sync begin
    @async put!(ch, x^2)
    take!(ch)
  end
end

@test gradient(tasks3, 5) == (10,)

tasks4(x) = fetch(@async x^2)

@test gradient(tasks4, 5) == (10,)

VERSION > v"1.3-" && include("threads.jl")


text:
using Statistics
using Plots
using FFTW
using Statistics
using Optim
using Images, FileIO, ImageIO
using Printf
using Revise
using Profile
using LinearAlgebra
using JLD2
using Random
using Distributions
using FITSIO
using StatsBase
using LineSearches
using SparseArrays

push!(LOAD_PATH, pwd()*"/main")
using DHC_2DUtils
push!(LOAD_PATH, pwd()*"/scratch_NM")
using Deriv_Utils_New
using Data_Utils
using Visualization
using ReconFuncs

logbool=false
apdbool = true
isobool=false
invcovhandling = "Full"

datfile = "scratch_NM/StandardizedExp/Nx64/Data_1000.jld2" #Replace w SLURM array
loaddf = load(datfile)
true_img = loaddf["true_img"]
init = loaddf["init"]
Nx=size(true_img)[1]
filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
Nf = size(filter_hash["filt_index"])[1]
dhc_args = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>apdbool, :iso=>isobool)
if dhc_args[:iso]
    coeff_mask = falses(2+S1iso+S2iso)
    coeff_mask[S1iso+3:end] .= true
else #Not iso
    coeff_mask = falses(2+Nf+Nf^2)
    coeff_mask[Nf+3:end] .= triu(trues(Nf, Nf))[:]
end
recon_settings = Dict([("log", logbool), ("covar_type", "sfd_dbn"), ("target_type", "sfd_dbn"), ("Invcov_matrix", invcovhandling), ("epsvalue", 1e-10)])
println("Apdbool", apdbool)
println("Isobool", isobool)
println(recon_settings)

s2mean, s2icov = meancov_generator(true_img, filter_hash, dhc_args, coeff_mask, recon_settings, safety=nothing)
invvar1 = diag(s2icov)

#Check wrt calculation here
dbnimg = readsfd(Nx, logbool=logbool)
dbn = get_dbn_coeffs(dbnimg, filter_hash, dhc_args, coeff_mask = nothing)
dbntriu = dbn[:, coeff_mask]
fmean = mean(dbntriu, dims=1)
fcov = (dbntriu .- fmean)' * (dbntriu .- fmean) ./ 9999
println("CovMat Dynamic Range: ")
println("Max: ", maximum(fcov))
println("Min: ", minimum(fcov))
emax = eigmax(fcov)
emin = eigmin(fcov)
evals = eigvals(fcov)
findall(evals.>1e-6)
svd_out = svd(fcov)
lowrank_approx = svd_out.U[1:end, 1:497] * Diagonal(svd_out.S[1:497]) * svd_out.V'[1:497, 1:end]
lowrank_approx
res = sum((lowrank_approx .- fcov).^2)
sum(svd_out.S[498:end].^2)
cond(fcov)

fullrank_approx = svd_out.U * Diagonal(svd_out.S) * svd_out.V
res = sum((fullrank_approx .- fcov).^2)

fcov -> add eps to the diag -> invcov -> diag(invcov) ?= std(dbncoeffs).^(-2)

end

randmat = randn(1000, 595)
prod = randmat' * randmat ./999
svdrand = svd(prod)
frrand = svdrand.U * Diagonal(svdrand.S) * svdrand.V'
sum((frrand .- prod).^2)

dbntriu_scaled = standardize(UnitRangeTransform, dbntriu, dims=1)
fmean_sc = mean(dbntriu_scaled, dims=1)
fcovsc = (dbntriu_scaled .- fmean_sc)' * (dbntriu_scaled .- fmean_sc) ./ 9999
findmax(dbntriu_scaled[1:end, 1])
findmax(dbntriu_scaled[1:end, 2])
maxcoeff_imgs = (x->findmax(dbntriu_scaled[1:end, x])).(collect(1:595))
badimg  =hcat((x->[x[1], x[2]]).(maxcoeff_imgs)...)[2, 1:end]
badimg = convert(Array{Int16}, badimg)
for img in unique(badimg)
    p= heatmap(dbnimg[:, :, img])
    savefig(p, "scratch_NM/NewWrapper/4-12/Img_"*string(img))
end

#Remove outlier images
imgmask = trues(10000)
imgmask[unique(badimg)] .= false
dbn = get_dbn_coeffs(dbnimg[:, :, imgmask], filter_hash, dhc_args, coeff_mask = nothing)
dbntriu = dbn[:, coeff_mask]
fmean = mean(dbntriu, dims=1)
fcov = (dbntriu .- fmean)' * (dbntriu .- fmean) ./ 9984
svd_out = svd(fcov)
fullrank_approx = svd_out.U * Diagonal(svd_out.S) * svd_out.V
res = sum((fullrank_approx .- fcov).^2)

invvar2 = std(dbn, dims=1)[:]
println(size(invvar2))
invvar2 = (invvar2.^(-2))[coeff_mask]
println("InvVar1", invvar1[:10])
println("InvVar2", invvar2[:10])
println("Mean Abs Error = ", mean(abs.(invvar1 .- invvar2)))
println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")

dhc_argsnorm = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>apdbool, :iso=>isobool, :norm=>true)
(N1iso,Nf) = size(filter_hash["S1_iso_mat"])
(N2iso,Nfsq) = size(filter_hash["S2_iso_mat"])
Ncovsamp = size(dbnimg)[3]
s20_dbn = zeros(Float64, Ncovsamp, 2+ (dhc_args[:iso] ? N1iso+N2iso : Nf+Nf^2))
for idx=1:Ncovsamp
    s20_dbn[idx, :] = DHC_compute_wrapper(dbnimg[:, :, idx], filter_hash; dhc_argsnorm...)
end
if coeff_mask!=nothing
    return s20_dbn[:, coeff_mask]
else
    return s20_dbn
end

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC