text:
using LinearAlgebra
using Distributions
using SpecialFunctions
using PDMats
using Printf

include("kernel.jl")
include("examples.jl")

"""
Define prediction/inference problem by supplying known parameters, including design matrices,
observation locations, prediction locations, nonlinear transform, prior marginals on theta and lambda,
and anticipated ranges for theta and lambda. Returns Bayesian predictive density function,
namely the integral of p(z_0| theta, lambda, z)*p(theta, lambda | z) over Theta x Lambda
(See Equations 8 and 12)

INPUTS:
X: k x p design matrix
X0: n x p design matrix
s: observed locations
s0: prediction locations
g: nonlinear transform parametrized by lambda. Usage: g(lambda, z)
gprime: derivative of g. Usage: gprime(z, lambda)
ptheta: prior marginal on theta; function handle
plambda: prior marginal on lambda; function handle
theta: structural parameter
lambda: transform parameter
z: observed random vector
z0: unobserved random vector

OUTPUTS:
probability density function z0 -> f(z_0|z)
cumulative density function z0 -> F(z0|z)
"""

function model(setting, g, gprime, pθ, pλ, rangeθ, rangeλ)
    s = setting.s; s0 = setting.s0; X = setting.X; X0 = setting.X0; z = setting.z;
    n = size(X, 1)
    p = size(X, 2)
    k = size(X0, 1)

    function func(θ)
        Eθ = K(s0, s0, θ, rbf)
        Σθ = K(s, s, θ, rbf)
        Bθ = K(s0, s, θ, rbf)
        choleskyΣθ = cholesky(Σθ) #precompute Cholesky decomposition of Sigma
        choleskyXΣX = cholesky(Hermitian(X'*(choleskyΣθ\X))) #precompute Cholesky decomposition of XSigma\X
        Dθ = Eθ - Bθ*(choleskyΣθ\Bθ')
        Hθ = X0 - Bθ*(choleskyΣθ\X)
        Cθ = Dθ + Hθ*(choleskyXΣX\Hθ')
        args = [Eθ, Σθ, Bθ, choleskyΣθ, choleskyXΣX, Dθ, Hθ, Cθ]
        return args
    end
    function density(λ, z0, args)
        Eθ, Σθ, Bθ, choleskyΣθ, choleskyXΣX, Dθ, Hθ, Cθ = args
        βhat = (X'*(choleskyΣθ\X))\(X'*(choleskyΣθ\g(z, λ)))
        qtilde = (expr = g(z, λ)-X*βhat; expr'*(choleskyΣθ\expr))
        m = Bθ*(cholesky(Σθ)\g(z, λ)) + Hθ*βhat
        J = λ -> abs(reduce(*, map(x -> gprime(x, λ), z)))
        post = det(cholesky(Σθ))^(-0.5)*det(choleskyXΣX)^(-0.5)*qtilde[1]^(-(n-p)/2)*J(λ)^(1-p/n)*pθ(θ)*pλ(λ)
        jac = abs(reduce(*, map(x -> gprime(x, λ), z0)))
        t = LocationScale(m[1], sqrt(qtilde[1]*Cθ[1]/(n-p)), TDist(n-p))
        return post*Distributions.pdf(t, g(z0, λ))*jac, post*Distributions.cdf(t, g(z0, λ))*jac
    end
    #simultaneously define PDF and CDF - (reuse integration nodes and weights)
    (z0 ->  int1D(θ -> (args = func(θ); int1D(λ -> density(λ, z0, args)[1], rangeλ)), rangeθ),
    z0 ->  int1D(θ -> (args = func(θ); int1D(λ -> density(λ, z0, args)[2], rangeλ)), rangeθ))
end

if false

end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module TestExamples1DCurved

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "1d")

@testset "Curved Mesh" begin
  @testset "elixir_advection_basic_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic_curved.jl"),
      l2   = [6.0388296447998465e-6],
      linf = [3.217887726258972e-5])
  end

  @testset "elixir_advection_nonperiodic_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonperiodic_curved.jl"),
      l2   = [5.641921365468918e-5],
      linf = [0.00021049780975179733])
  end

  @testset "elixir_euler_source_terms_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_curved.jl"),
      l2   = [2.2527950196212703e-8, 1.8187357193835156e-8, 7.705669939973104e-8],
      linf = [1.6205433861493646e-7, 1.465427772462391e-7, 5.372255111879554e-7])
  end

  @testset "elixir_euler_nonperiodic_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_nonperiodic_curved.jl"),
      l2   = [3.8099996914101204e-6, 1.6745575717106341e-6, 7.732189531480852e-6],
      linf = [1.2971473393186272e-5, 9.270328934274374e-6, 3.092514399671842e-5])
  end
end

end # module

ext: jl


text:
using Statistics
using ImageTransformations: center
using Photometry
using HCIToolbox
using Dierckx
using ProgressLogging
using StaticKernels
using LinearAlgebra: dot

"""
    contrast_curve(alg, cube, angles, psf, args...;
                   fwhm, sigma=5, nbranch=1, theta=0, inner_rad=1,
                   starphot=Metrics.estimate_starphot(cube, fwhm),
                   fc_rad_sep=3, snr=100, k=2, smooth=true,
                   subsample=true, kwargs...)

Calculate the throughput-calibrated contrast. This first processes the algorithmic [`throughput`](@ref) by injecting instances of `psf` into `cube`. These are processed through `alg` and the ratio of the recovered flux to the injected flux is calculated. These companions are injected in resolution elements across the frame, which can be changed via the various keyword arguments.

The throughput can only be calculated for discrete resolution elements, but we typically want a much smoother curve. To accomplish this, we measure the noise (the standard deviation of all resolution elements in an annulus at a given radius) for every pixel in increasing radii. We then interpolate the throughput to this grid and return the subsampled curves.

# Returned Fields
* `distance` - The radial distance (in pixels) for each measurement
* `contrast` - The Gaussian sensitivity
* `contrast_corr` - The Student-t sensitivity
* `noise` - The noise measured for each distance
* `throughput` - The throughput measured for each distance.

# Keyword Arguments
* `sigma` - The confidence limit in terms of Gaussian standard deviations
* `starphot` - The flux of the star. By default calculates the flux in the central core.
**Injection Options** (See also [`throughput`](@ref))
* `nbranch` - number of azimuthal branches to use
* `theta` - position angle of initial branch
* `inner_rad` - position of innermost planet in FWHM
* `fc_rad_sep` - the separation between planets in FWHM for a single reduction
* `snr` - the target signal to noise ratio of the injected planets
**Subsampling Options** (See also [`Metrics.subsample_contrast`](@ref))
* `subsample` - If true, subsamples the throughput measurements to increase density of curve
* `k` - The order of the BSpline used for subsampling the throughput
* `smooth` - If true, will smooth the subsampled noise measurements with a 2nd order Savitzky-Golay filter

!!! tip
    If you prefer a tabular format, simply pipe the output of this function into any type supporting the Tables.jl interface, e.g.
    ```
    contrast_curve(alg, cube, angles, psf; fwhm=fwhm) |> DataFrame
    ```
"""
function contrast_curve(alg, cube, angles, psf, args...;
        fwhm, sigma=5, nbranch=1, theta=0, inner_rad=1,
        starphot=Metrics.estimate_starphot(cube, fwhm),
        fc_rad_sep=3, snr=100, k=2, smooth=true,
        subsample=true, kwargs...)

    # measure the noise and throughput in consecutive resolution elements
    # across azimuthal branches
    @info "Calculating Throughput"
    reduced_empty = alg(cube, angles, args...; kwargs...)

    through, meta = throughput(alg, cube, angles, psf, args...;
                               fwhm=fwhm, inner_rad=inner_rad, fc_rad_sep=fc_rad_sep, theta=theta,
                               nbranch=nbranch, snr=snr, reduced_empty=reduced_empty, kwargs...)

    through_mean = mean(through, dims=2) |> vec

    if subsample
        return subsample_contrast(reduced_empty, meta.distance, through_mean;
                                  fwhm=fwhm, starphot=starphot, sigma=sigma, inner_rad=inner_rad,
                                  theta=theta, smooth=smooth, k=k)
    end

    # calculate common terms once
    unit_contrast = @. meta.noise / (through_mean * starphot)
    # gaussian contrast (invalid values become NaN)
    contrast = calculate_contrast.(sigma, unit_contrast)

    # get correction for small-sample statistics
    sigma_corr = correction_factor.(meta.distance, fwhm, sigma)
    # student-t contrast (invalid values become NaN)
    contrast_corr = calculate_contrast.(sigma_corr, unit_contrast)

    return (distance=meta.distance,
            throughput=through_mean,
            contrast=contrast,
            contrast_corr=contrast_corr,
            noise=meta.noise)
end

function correction_factor(radius, fwhm, sigma)
    n_res_els = 2 * π * radius ÷ fwhm
    ss_corr = sqrt(1 + 1 / (n_res_els - 1))
    return quantile(TDist(n_res_els), cdf(Normal(), sigma)) * ss_corr
end

"""
    Metrics.subsample_contrast(empty_frame, distance, throughput;
                               fwhm, starphot, sigma=5, inner_rad=1,
                               theta=0, smooth=true, k=2)

Helper function to subsample and smooth a contrast curve.

Contrast curves, by definition, are calculated with discrete resolution elements. This can cause contrast curves to have very few points instead of appearing as a continuously measured statistic across the pixels. We alleviate this by sub-sampling the throughput (via BSpline interpolation) across each pixel (instead of each resolution element).

The noise can be found efficiently enough, so rather than interpolate we measure the noise in annuli of width `fwhm` increasing in distance by 1 pixel. We measure this noise in `empty_frame`, which should be a 2D reduced ADI frame.

The noise measurements can be noisy, so a 2nd order Savitzky-Golay filter can be applied via `smooth`. This fits a quadratic polynomial over a window of `fwhm/2` points together to reduce high-frequency jitter.

# Examples

Here is an example which calculates the exact contrast curve in addition to a subsampled curve without re-calculating the throughput.

```julia
cube, angles, psf = # load data

alg = PCA(10)
cc = contrast_curve(alg, cube, angles, psf; fwhm=8.4, subsample=false)
reduced_empty = alg(cube, angles)
cc_sub = Metrics.subsample_contrast(reduced_empty, cc.distance, cc.throughput; fwhm=8.4)
```
"""
function subsample_contrast(empty_frame, distance, throughput;
                            fwhm, starphot, sigma=5, inner_rad=1,
                            theta=0, smooth=true, k=2)
    # measure the noise with high sub-sampling-
    # at every pixel instead of every resolution element
    cy, cx = center(empty_frame)
    radii_subsample = first(distance):1:last(distance) + 1
    through_subsample = Spline1D(distance, throughput, k=k)(radii_subsample)
    noise_subsample = @. annulus_noise((empty_frame,), fwhm, cy, cx, radii_subsample, theta)

    if smooth
        window_size = min(length(noise_subsample) - 2, round(Int, 2 * fwhm))
        iseven(window_size) && (window_size += 1)
        width = window_size ÷ 2
        coeffs = savgol_coeffs(window_size, 2) |> Tuple
        smooth_kernel = Kernel{(-width:width,)}(w -> dot(coeffs, Tuple(w)))
        noise_smoothed = map(smooth_kernel, extend(noise_subsample, StaticKernels.ExtensionReplicate()))
    else
        noise_smoothed = noise_subsample
    end

    # calculate common terms once
    unit_contrast = @. noise_smoothed / (through_subsample * starphot)
    # gaussian contrast (invalid values become NaN)
    contrast = calculate_contrast.(sigma, unit_contrast)

    # get correction for small-sample statistics
    sigma_corr = correction_factor.(radii_subsample, fwhm, sigma)
    # student-t contrast (invalid values become NaN)
    contrast_corr = calculate_contrast.(sigma_corr, unit_contrast)

    return (distance=radii_subsample,
            throughput=through_subsample,
            contrast=contrast,
            contrast_corr=contrast_corr,
            noise=noise_smoothed)
end

# simple function to inline the contrast calculation _and_ clipping
@inline function calculate_contrast(k, unit_contrast)
    contrast = k * unit_contrast
    return 0 ≤ contrast ≤ 1 ? contrast : NaN
end

function savgol_coeffs(window_size, order=2)
    pos = window_size ÷ 2
    # form vandermonde matrix
    x = -pos:window_size - pos - 1
    powers = 0:order
    A = x' .^ powers
    # solve least squares equation Ac = y
    y = zeros(order + 1)
    y[1] = 1
    return A \ y
end

"""
    Metrics.estimate_starphot(cube, fwhm)
    Metrics.estimate_starphot(frame, fwhm)

Simple utility to estimate the stellar photometry by placing a circular aperture with `fwhm` diameter in the center of the `frame`. If a cube is provided, first the median frame will be found.
"""
function estimate_starphot(frame::AbstractMatrix, fwhm)
    ap = CircularAperture(reverse(center(frame)), fwhm/2)
    return photometry(ap, frame).aperture_sum
end

estimate_starphot(cube::AbstractArray{T, 3}, fwhm) where {T} = estimate_starphot(collapse(cube, method=median), fwhm)

"""
    throughput(alg, cube, angles, psf, args...;
               fwhm, nbranch=1, theta=0, inner_rad=1,
               fc_rad_sep=3, snr=100, kwargs...)

Calculate the throughput of `alg` by injecting fake companions into `cube` and measuring the relative photometry of each companion in the reduced frame. Any additional `args` or `kwargs` will be passed to `alg` when it is called.

# Keyword Arguments
* `nbranch` - number of azimuthal branches to use
* `theta` - position angle of initial branch
* `inner_rad` - position of innermost planet in FWHM
* `fc_rad_sep` - the separation between planets in FWHM for a single reduction
* `snr` - the target signal to noise ratio of the injected planets
* `reduced_empty` - the collapsed residual frame for estimating the noise. Will process using `alg` if not provided.
"""
function throughput(alg, cube::AbstractArray{T,3}, angles, psf_model, args...;
                    fwhm, nbranch=1, theta=0, inner_rad=1, fc_rad_sep=3,
                    snr=100, reduced_empty = nothing, kwargs...) where T
    maxfcsep = size(cube, 2) ÷ (2 * fwhm) - 1
    # too large separation between companions in the radial patterns
    3 ≤ fc_rad_sep ≤ maxfcsep || error("`fc_rad_sep` should lie ∈[3, $(maxfcsep)], got $fc_rad_sep")

    # compute noise in concentric annuli on the empty frame
    if isnothing(reduced_empty)
        reduced_empty = alg(cube, angles, args...; kwargs...)
    end

    cy, cx = center(reduced_empty)

    n_annuli = floor(Int, (cy - fwhm) / fwhm) - 1
    radii = fwhm .* (inner_rad:n_annuli)
    δy, δx = sincosd(theta)
    noise = @. annulus_noise((reduced_empty,), fwhm, cy, cx, radii, theta)

    angle_per_branch = 360 / nbranch
    output = similar(cube, length(radii), nbranch)

    fake_comps_full = zero(reduced_empty)
    @progress "branch" for branch in 1:nbranch
        θ = theta + angle_per_branch * (branch - 1)
        @progress "pattern" for init_rad in 1:fc_rad_sep
            slice = init_rad:fc_rad_sep:lastindex(radii)
            fake_comps = zero(reduced_empty)

            cube_fake_comps = copy(cube)

            apertures = map(slice) do ann
                r = radii[ann]
                δy, δx = sincosd(θ)
                x = r * δx + cx
                y = r * δy + cy

                A = snr * noise[ann]

                inject!(fake_comps, psf_model; A=A, r=r, θ=θ)
                fake_comps_full .+= fake_comps
                inject!(cube_fake_comps, psf_model, angles; A=A, r=r, θ=θ)

                return CircularAperture(x, y, fwhm / 2)
            end
            reduced = alg(cube_fake_comps, angles, args...; kwargs...)

            injected_flux = photometry(apertures, fake_comps).aperture_sum
            recovered_flux = photometry(apertures, reduced .- reduced_empty).aperture_sum
            @. output[slice, branch] = max(zero(T), recovered_flux / injected_flux)
        end
    end

    return output, (distance=radii, fake_comps=fake_comps_full, noise=noise)
end

function annulus_noise(frame, fwhm, cy, cx, r, θ=0)
    δy, δx = sincosd(θ)
    x = r * δx + cx
    y = r * δy + cy
    Metrics.noise(frame, (x, y), fwhm)
end


text:
nt = 1000
m = ABCDQR_simula(nt)

y = m["y"]
u = m["u"]
A = m["A"]
B = m["B"]
C = m["C"]
D = m["D"]
Q = m["Q"]
R = m["R"]
x10 = 0*A[:,1]
P10 = 0*A

###########################
print("Testing ABCDQR_kfilter: ")
xtt,Ptt,xtt1,Ptt1,et,St,Kt,loglik = ABCDQR_kfilter(y,u,A,B,C,D,Q,R,x10,P10)
@test sum(xtt.^2) > 1e-6
println("OK")

#######################
print("Testing ABCDQR_kfilter_s: ")
xtt,Ptt,xtt1,Ptt1,et,St,Kt,loglik = ABCDQR_kfilter_s(y,u,A,B,C,D,Q,R,x10)
@test sum(xtt.^2) > 1e-6
println("OK")


text:
# Bland's rule
function pivot_Bland(spx::SpxData)
    spx.enter = -1

    # compute reduced cost
    compute_reduced_cost(spx)
    # @show spx.r
    # @show spx.nonbasic

    min_enter = ncols(spx.lpdata) + 1
    for j = 1:length(spx.r)
        if spx.nonbasic[j] < min_enter
            if spx.basis_status[spx.nonbasic[j]] == BASIS_AT_UPPER
                if spx.r[j] > 1e-10
                    spx.enter = spx.nonbasic[j]
                    spx.enter_pos = j
                    min_enter = spx.enter
                end
            else
                if spx.r[j] < -1e-10
                    spx.enter = spx.nonbasic[j]
                    spx.enter_pos = j
                    min_enter = spx.enter
                end
            end
        end
    end
    # @show spx.enter, spx.r[spx.enter_pos]
end

# Steepest-edge rule
"""
Steepest-edge rule:

This is based on the following references.
1. Bieling et al. An efficient GPU implementation of the revised simplex method. 2010
2. Forest and Goldfarb. Steepest-edge simplex algorithms for linear programming. 1992
3. Goldfarb and Reid. A practicable steepest-edge simplex algorithm. 1977
"""
function pivot_steepest_edge(spx::SpxData)

    compute_reduced_cost(spx)
    # @show spx.r
    # @show spx.nonbasic

    compute_steepest_edge_weights(spx)

    # compute the local slope
    # @show length(spx.r), length(spx.gamma)
    spx.s .= (spx.r).^2 ./ spx.gamma

    spx.enter = -1
    max_s = 0.0
    for j = 1:length(spx.r)
        if spx.s[j] > max_s
            if spx.basis_status[spx.nonbasic[j]] == BASIS_AT_UPPER && sign(spx.r[j]) > 0
                spx.enter = spx.nonbasic[j]
                spx.enter_pos = j
                max_s = spx.s[j]
            elseif spx.basis_status[spx.nonbasic[j]] == BASIS_AT_LOWER && sign(spx.r[j]) < 0
                spx.enter = spx.nonbasic[j]
                spx.enter_pos = j
                max_s = spx.s[j]
            end
        end
    end
    # @show spx.enter, s[spx.enter_pos]
end

function compute_steepest_edge_weights(spx::SpxData)
    @timeit TO "compute steepest edge weights" begin
        if spx.iter == 1
            for j = 1:length(spx.nonbasic)
                spx.g .= spx.invB * spx.lpdata.A[:,spx.nonbasic[j]]
                spx.gamma[j] = spx.g' * spx.g
            end
        else
            # Here, spx.enter and spx.leave were obtained from the previous iteration.
            @timeit TO "gamma_j" begin
                spx.gamma .+= 2 .* (spx.lpdata.A[:,spx.nonbasic]' * spx.invB[:,spx.leave]) .* (spx.lpdata.A[:,spx.nonbasic]' * spx.v) .+ (spx.lpdata.A[:,spx.nonbasic]' * spx.invB[:,spx.leave]).^2 .* (1 + spx.d' * spx.d)
            end

            @timeit TO "gamm_e" begin
                spx.g .= spx.invB * spx.lpdata.A[:,spx.nonbasic[spx.enter_pos]]
                spx.gamma[spx.enter_pos] = spx.g' * spx.g
            end
        end
    end
end

# Dantzig's rule
function pivot_Dantzig(spx::SpxData)
    spx.enter = -1

    # compute reduced cost
    compute_reduced_cost(spx)

    max_r = 0.0
    for j = 1:length(spx.r)
        if abs(spx.r[j]) > max_r
            if spx.basis_status[spx.nonbasic[j]] == BASIS_AT_UPPER
                if spx.r[j] > 1e-6
                    spx.enter = spx.nonbasic[j]
                    spx.enter_pos = j
                    max_r = spx.r[j]
                end
            else
                if spx.r[j] < -1e-6
                    spx.enter = spx.nonbasic[j]
                    spx.enter_pos = j
                    max_r = -spx.r[j]
                end
            end
        end
    end
    # @show spx.enter, spx.r[spx.enter_pos]
end

# Pivot to remove artificial variables from basis (based on Bland's rule)
function pivot_to_remove_artificials(spx::SpxData)
    spx.enter = -1

    # compute reduced cost
    compute_reduced_cost(spx)

    min_enter = spx.start_artvars
    for j = 1:length(spx.r)
        if spx.nonbasic[j] < min_enter
            if spx.basis_status[spx.nonbasic[j]] == BASIS_AT_UPPER
                if spx.r[j] > -1e-10
                    spx.enter = spx.nonbasic[j]
                    spx.enter_pos = j
                    min_enter = spx.enter
                end
            else
                if spx.r[j] < +1e-10
                    spx.enter = spx.nonbasic[j]
                    spx.enter_pos = j
                    min_enter = spx.enter
                end
            end
        end
    end
end


text:
# Formats package and version into a standard string
package_name(package::Package) = package_name(package.name, package.version)
package_name(package::String, version::VersionNumber) = "$package-$version"
# Creates standard work-directory name
builddir(package::Package, workdir::String) = abspath(joinpath(workdir, package_name(package)))
# Directory where the package is located with debian files and such
function packagedir(package::Package, workdir::String)
  joinpath(builddir(package, workdir), package_name(package))
end
# Creates standard tar-file name
function tarfile(package::Package, workdir::String)
  const version = VersionNumber(package.version.major, package.version.minor, package.version.patch)
  joinpath(builddir(package, workdir), "$(package.name)_$(version).orig.tar.gz")
end
function sourcedir(package::Package, workdir::String)
  joinpath(builddir(package, workdir), "source", package_name(package))
end

# Converts from kwargs for list of dictionaries
# Each dictionary specifies a binary package
function get_binaries(name::String; kwargs...)
  # Transform variable keyword arguments into binary packages
  current = Dict{Symbol, Any}()
  binaries = Dict{Symbol, Any}[]

  # Adds to list of packages, making sure a package name is available
  function add_package(current, binaries)
    if !haskey(current, :package); current[:package] = name end
    push!(binaries, current)
  end

  for (key, value) in kwargs
    if key == :package && length(current) > 0
      add_package(current, binaries)
      current = {key => value}
    else
      current[key] = value
    end
  end
  if length(current) > 0; add_package(current, binaries) end
  if length(binaries) == 0; push!(binaries, {:package => name}) end
  binaries
end

max_stars_repo_name: UCL/DebbyPacker

max_issues_repo_name: UCL/DebbyPacker

max_forks_repo_name: UCL/DebbyPacker


text:
# policy file better sanity fewer contacts
adaptivepolicythresholdoff = 0
adaptivepolicythresholdon = 9999999999999

socialmaxyyt = socialmaxyy_o # max no. people met outside firm young_young
socialmaxoot = socialmaxoo_o # max no. people met outside firm old_old
socialmaxoyt = socialmaxoy_o
socialmaxyot = socialmaxyo_o

genhomeoffice = false

ptravel = ptravel_o

pshopt = [[1,1,1],[1,1,1]]

tauadj = 0.05

shorttimeprogram = false


text:
## Util.jl##

if VERSION>=v"0.4.0-rc1"
    tupletypelength(a)=length(a.parameters)
else
    tupletypelength(a)=length(a)
end

macro replaceBase(fs...)
    b = Expr(:block)
    for f in fs
        if f == :./
            fa = :fdiv
        elseif f == :.^
            fa = :pow
        elseif f == :.+
            fa = :vadd
        elseif f == :.-
            fa = :vsub
        elseif f == :.*
            fa = :vmul
        else
            fa = f
        end
        if fa in (:ceil,:floor,:sqrt,:rsqrt,:rec,
                  :exp,:exp2,:expm1,:log,:log1p,:log2,:log10,
                  :sin,:sinpi,:cos,:cospi,:tan,:tanpi,:asin,:acos,:atan,
                  :sinh,:cosh,:tanh,:asinh,:acosh,:atanh,
                  :trunc,:round,:exponent,:abs,:sincos,:cis)
            e = quote
                (Base.$f)(X::Array{Float64}) = ($fa)(X)
                (Base.$f)(X::Array{Float32}) = ($fa)(X)
            end
        elseif fa in (:copysign,:atan2,:pow,:rem,:fdiv, :vadd, :vsub, :vmul)
            e = quote
                (Base.$f)(X::Array{Float64},Y::Array{Float64}) = ($fa)(X,Y)
                (Base.$f)(X::Array{Float32},Y::Array{Float32}) = ($fa)(X,Y)
            end
        else
            error("Function $f not defined by AppleAccelerate.jl")
        end
        push!(b.args,e)
    end
    b
end


text:
# Counting Gaps and Coverage
# --------------------------

"""
It calculates the fraction of gaps on the `Array` (alignment, sequence, column, etc.). This
function can take an extra `dimension` argument for calculation of the gap fraction over the
given dimension.
"""
function gapfraction(x::AbstractArray{Residue})
    counter = 0
    len = 0
    for res in x
        counter += ifelse(res == GAP, 1, 0)
        len += 1
    end
    float(counter) / float(len)
end

function gapfraction(x::AbstractArray{Residue}, dimension::Int)
    mapslices(gapfraction, x, dims=dimension)
end

"""
It calculates the fraction of residues (no gaps) on the `Array`
(alignment, sequence, column, etc.). This function can take an extra `dimension` argument
for calculation of the residue fraction over the given dimension
"""
residuefraction(x::AbstractArray{Residue}) = 1.0 - gapfraction(x)

function residuefraction(x::AbstractArray{Residue}, dimension::Int)
    mapslices(residuefraction, x, dims=dimension)
end

macro keep_names_dimension(functions)
    function_names = functions.args
    n = length(function_names)
    definitions = Array{Any}(undef, n)

    for i in 1:n
        f = esc(function_names[i])
        definitions[i] = quote

            function ($f)(msa::NamedResidueMatrix{T}, dimension::Int) where T
                result = ($f)(getarray(msa), dimension)
                if dimension == 1
                    name_list = names(msa,2)
                    N = length(name_list)
                    NamedArray(result,
                        (OrderedDict{String,Int}(Utils._get_function_name(string($f))=>1),
                         OrderedDict{String,Int}(name_list[i]=>i for i in 1:N)),
                        ("Function","Col"))
                elseif dimension == 2
                    name_list = names(msa,1)
                    N = length(name_list)
                    NamedArray(result,
                        (OrderedDict{String,Int}(name_list[i]=>i for i in 1:N),
                         OrderedDict{String,Int}(Utils._get_function_name(string($f))=>1)),
                        ("Seq","Function"))
                else
                    throw(ArgumentError("Dimension must be 1 or 2."))
                end
            end

            ($f)(a::AbstractAlignedObject, dimension::Int) = ($f)(namedmatrix(a), dimension)
        end
    end

    return Expr(:block, definitions...)
end

@keep_names_dimension([gapfraction, residuefraction])

"Coverage of the sequences with respect of the number of positions on the MSA"
function coverage(msa::AbstractMatrix{Residue})
    result = residuefraction(msa, 2)
    if isa(result, NamedArray) && ndims(result) == 2
        setnames!(result,["coverage"],2)
    end
    result
end

coverage(msa::AbstractAlignedObject) = coverage(namedmatrix(msa))

"Fraction of gaps per column/position on the MSA"
columngapfraction(msa::AbstractMatrix{Residue}) = gapfraction(msa, 1)
columngapfraction(msa::AbstractAlignedObject) = columngapfraction(namedmatrix(msa))


text:
module SymbolicUtils
using SymEngine
using CommonSubexpressions
using MacroTools: postwalk
using ..RocketlandDefns

struct Deriv
	fn::Symbol
	derv::Vector{Int}
end
Base.isequal(a::Deriv, b::Deriv) = (a.fn == b.fn) && (a.derv == b.derv)
Base.:(==)(a::Deriv, b::Deriv) = isequal(a,b)
Base.hash(a::Deriv, h::UInt=0) = hash(a.derv, hash(a.fn, h))
Base.hash(a::Deriv, h::Int=0) = hash(a, convert(UInt, h))

function diff_fun(f::Function, nargs::Int, funs_and_derivs::Dict{Deriv, Symbol})
	args = [symbols("x$i") for i=1:nargs]
	symbolic = f(args...)
	return symbolic
end

fs_rgx = r"(.*?)\("
function get_fnsym_name(ex::Basic)
	return match(fs_rgx, SymEngine.toString(ex))[1]
end

function _convert_deriv(ex::Basic, derivs::Dict{Deriv, Symbol}, substs::Dict{Symbol, Any})
	args = SymEngine.get_args(ex)
	call = args[1]
	fnd = Symbol(get_fnsym_name(call))
	fnargs = SymEngine.get_args(call)
	argo = Dict{Symbol, Int64}(map(p->(Symbol(p[2]),p[1]), enumerate(fnargs)))
	sig = [argo[Symbol(var)] for var in args[2:end]]
	key = Deriv(fnd, sig)
	#dfun = derivs[key]

	return Expr(:ref, Expr(:call, fnd, Val{:jac}, [_convert_diff(a, derivs, substs) for a in fnargs]...), sig[1])
end

SymEngine.N(b::SymEngine.BasicType{Val{:Integer}}) = convert(Float64, convert(BigInt, b))
SymEngine.N(b::SymEngine.BasicType{Val{:Rational}}) = N(numerator(b))/N(denominator(b))
function _convert_diff(ex::Basic, derivs::Dict{Deriv, Symbol}, substs::Dict{Symbol, Any})
    fn = SymEngine.get_symengine_class(ex)
    if fn == :Symbol
    	sym = Symbol(SymEngine.toString(ex))
    	if haskey(substs, sym)
    		return substs[sym]
    	else
    		return sym
    	end
    elseif (fn in SymEngine.number_types) || (fn == :Constant)
        return N(ex)
    elseif fn == :Pow
    	as = SymEngine.get_args(ex)
    	bdy = _convert_diff(as[1], derivs, substs)
    	pow = _convert_diff(as[2], derivs, substs)
    	if pow == 0.5
    		return :(sqrt($bdy))
    	elseif pow == -0.5
    		return :(1.0/sqrt($bdy))
    	elseif pow == 1.0
    		return bdy
    	elseif pow == 2.0
    		return :($bdy ^ 2)
    	else
    		return :($bdy ^ $pow)
    	end
    elseif fn == :FunctionSymbol
    	fn = Symbol(get_fnsym_name(ex))
    elseif fn == :Derivative
    	return _convert_deriv(ex, derivs, substs)
    elseif fn == :Subs
    	nsubsts = copy(substs)
    	args = SymEngine.get_args(ex)
    	nsubsts[_convert_diff(args[2], derivs, substs)] =
    		_convert_diff(args[3], derivs, substs)
    	return _convert_diff(args[1], derivs, nsubsts)
    end

    as = SymEngine.get_args(ex)

    Expr(:call, SymEngine.map_fn(fn, SymEngine.fn_map),
    	 [_convert_diff(a,derivs,substs) for a in as]...)
end

function convert_diff(ex::Basic, derivs::Dict{Deriv, Symbol}, substs::Dict{Symbol, Any})
    fn = SymEngine.get_symengine_class(ex)

    if fn == :Symbol
        return Expr(:call, :*, Symbol(SymEngine.toString(ex)), 1)
    elseif (fn in SymEngine.number_types) || (fn == :Constant)
        return SymEngine.N(ex)
    end

    return _convert_diff(ex, derivs, substs)
end

function convert_diff(exes::AbstractArray{Basic}, derivs::Dict{Deriv, Symbol}, substs::Dict{Symbol, Any})
	return map(expr -> convert_diff(expr, derivs, substs), exes)
end

function make_simplified(name, fun, nargs; postprocess = nothing, expected_args = Set{Symbol}())
	# symbolically execute the body
	args = [symbols("x$i") for i = 1:nargs]
	retv = fun(args)
	exprs = convert_diff(retv, Dict{Deriv,Symbol}(), Dict{Symbol,Any}())
	if exprs isa Array{T, 1} where T
		rexpr = Expr(:vect, exprs...)
	elseif exprs isa Array{T, 2} where T
		rexpr = Expr(:vect, reshape(exprs, length(exprs))...)
	else
		rexpr = Expr(:vect, exprs)
	end
	if !isnothing(postprocess)
		rexpr = postprocess(rexpr)
	end
	body = macroexpand(CommonSubexpressions,:(@cse $rexpr))

	# compute the extra arguments
	need_args = Set{Symbol}()
	map(x -> push!(need_args, x), vcat(map(x -> map(Symbol, SymEngine.free_symbols(x)), retv)...))
	needed_state = intersect(need_args, Set{Symbol}(map(Symbol, args)))
	map(x -> delete!(need_args, Symbol(x)), args)
	ordered_args = sort(collect(union(need_args, expected_args)))

	# generate the header
	header = map((idx,var) -> :($(Symbol(var)) = inp[$idx]),
		parse.(Int, map(x -> String(x)[2:end], collect(needed_state))), needed_state)

	# generate the in-place assignment
	if exprs isa Array{T, 1} where T
		assignment = map(x -> :(J[$(x[1])] = $(x[2])), enumerate(body.args[end].args))
	elseif exprs isa Array{T, 2} where T
		bargs = reshape(body.args[end].args, size(exprs))
		assignment = vcat([[:(J[$rown, $coln] = $(bargs[rown, coln]))
				for (rown, cell) in enumerate(col) if bargs[rown, coln] != 0.0]
				for (coln, col) in enumerate(eachcol(exprs))]...)
	else
		assignment = Expr[:(return $(body.args[end].args...))]
	end

	#output the simplified function
	fbody = (Expr(:block, header...,
				 body.args[1:end-1]...,
				 assignment...))

	if exprs isa Array
		return :($(name)(J, inp, $(ordered_args...)) = $(fbody.args...))
	else
		return :($(name)(inp, $(ordered_args...)) = $(fbody.args...))
	end
end

function make_simplified_array(name, fun, nargs; postprocess = nothing, expected_args = Set{Symbol}())
	# symbolically execute the body
	args = [symbols("x$i") for i = 1:nargs]
	retv = fun(args)
	exprs = convert_diff(retv, Dict{Deriv,Symbol}(), Dict{Symbol,Any}())
	if exprs isa Array{T, 1} where T
		rexpr = Expr(:vect, exprs...)
	elseif exprs isa Array{T, 2} where T
		rexpr = Expr(:vect, reshape(exprs, length(exprs))...)
	else
		rexpr = Expr(:vect, exprs)
	end
	if !isnothing(postprocess)
		rexpr = postprocess(rexpr)
	end
	body = macroexpand(CommonSubexpressions,:(@cse $rexpr))

	# compute the extra arguments
	need_args = Set{Symbol}()
	map(x -> push!(need_args, x), vcat(map(x -> map(Symbol, SymEngine.free_symbols(x)), retv)...))
	needed_state = intersect(need_args, Set{Symbol}(map(Symbol, args)))
	map(x -> delete!(need_args, Symbol(x)), args)
	ordered_args = sort(collect(union(need_args, expected_args)))

	# generate the header
	header = map((idx,var) -> :($(Symbol(var)) = inp[$idx]),
		parse.(Int, map(x -> String(x)[2:end], collect(needed_state))), needed_state)
	result = body.args[end]
	result.head = :call
	pushfirst!(result.args, :SVector)
	#output the simplified function
	fbody = (Expr(:block, header...,
				 body.args[1:end-1]...,
				 result))

	outexp = :($(name)(inp, $(ordered_args...)) = $(fbody.args...))
	#println(outexp)
	return outexp
end

function make_jacobian(fnname, fun, nargs)
	#=
	dfs = Dict{Deriv,Symbol}(
		[Deriv(c.args[2].args[1], [c.args[2].args[2:end]...]) => c.args[3] for c in idfs.args])
	vects = Dict{Symbol,Tuple{Symbol, Int}}(
		[c.args[2] => (c.args[3].args[1], c.args[3].args[2]) for c in ivects.args])
	;
		postprocess = mat -> postwalk(x -> x isa Expr && x.head == :call && haskey(vects, x.args[1]) ?
				begin (fn,idx)=vects[x.args[1]]; op = :($(fn)($(x.args[2:end]...))[$(idx)]); op end : x, mat)
		=#
	return make_simplified(fnname, x -> begin frs = fun(x);
		res = [diff(row, var) for row in frs, var in x];
		return res end, nargs)
end
#=

dp = DescentProblem()
pri = ProbInfo(dp)
@make_simplified(dxs, st -> dx(st[1:14], st[15:17], pri), 17)

@make_jacobian(dx, st -> dx(st[1:14], st[15:17], pri), 17,
	[ilift[1]=>ilift1, ilift[2]=>ilift2,
	 idrag[1]=>idrag1, idrag[2]=>idrag2],
	[ilift => ilift[1], ilift1 => ilift[2], ilift2 => ilift[3],
	 idrag => idrag[1], idrag1 => idrag[2], idrag2 => idrag[3]])
	 =#
end

max_stars_repo_name: BenChung/SuccessiveConvexification

max_issues_repo_name: BenChung/SuccessiveConvexification

max_forks_repo_name: BenChung/SuccessiveConvexification


text:
module MagneticReadHead

using Base: invokelatest
using Cassette
using MacroTools
using InteractiveUtils
using CodeTracking
# We don't use Revise, but if it isn't loaded CodeTracking has issues
using Revise: Revise
using OrderedCollections

export @run, @enter

include("utils.jl")
include("method_utils.jl")

include("breakpoint_rules.jl")

include("core_control.jl")
@nospecialize
include("pass.jl")
@specialize

include("inner_repl.jl")
include("break_action.jl")
include("locate.jl")
include("breakpoints.jl")

struct UserAbortedException <: Exception end

function iron_debug(debugbody)
    try
        ctx = new_debug_ctx()
        return Cassette.recurse(ctx, debugbody)
    catch err
        err isa UserAbortedException || rethrow()
        nothing
    finally
        # Disable any stepping left-over
        GLOBAL_STEPPING_MODE[] =  StepContinue
    end
end

"""
    @run the_code
Run until the_code until a breakpoint is hit.
"""
macro run(body)
    quote
        iron_debug() do
            $(esc(body))
        end
    end
end

"""
    @enter the_code
Begin debugging and break on the start of the_code.
"""
macro enter(body)
    if body isa Expr && body.head==:call && body.args[1] isa Symbol
        body = MacroTools.striplines(body)
        break_target = :(InteractiveUtils.which($(body.args[1]), Base.typesof($(body.args[2:end])...)))
        quote
            set_breakpoint!($(esc(break_target)))
            try
                iron_debug() do
                    $(esc(body))
                end
            finally
                rm_breakpoint!($(esc(break_target)))
            end
        end
    else
        quote
            error("Expression too complex to `@enter`. Please use `@run` with manual breakpoint set")
        end
    end
end

end # module


text:
using HiQGA.transD_GP
include("BarPhysics.jl")
include("BarPhysicsInversion.jl")
using .BarPhysicsInversion
# Forward setttings
# dummy physics struct defined in BarPhysics.jl
F = BarPhysics.Bar(25, rand(10))

# some dummy forward
@info BarPhysics.returnphysics!(F, rand(3))

# struct for data and physics operator defined in BarPhysicsInversion.jl
FOp = BarPhysicsInversion.BarInversion(rand(10), F)

## Inverse settings
# prior bounds and setings
fbounds = [0 1.]
xall = collect(permutedims(0:0.1:1))
xbounds = permutedims([extrema(xall)...])
λ = vec(diff(xbounds, dims=2))
nmin, nmax = 2, 5
demean, sampledc = false, true
# McMC proposals
δ = 0.1*diff(fbounds, dims=2)[1]
sdev_pos = abs.(vec(diff(xbounds, dims=2)))
sdev_prop = abs.(vec(0.05*diff(fbounds, dims=2)))
K = transD_GP.GP.OrstUhn()

# make MCMC options using prior and proposals
opt = transD_GP.OptionsStat(nmin = nmin,
                        nmax = nmax,
                        xbounds = xbounds,
                        fbounds = fbounds,
                        xall = xall,
                        λ = λ,
                        δ = δ,
                        demean = demean,
                        sampledc = sampledc,
                        sdev_prop = sdev_prop,
                        sdev_pos = sdev_pos,
                        quasimultid = false,
                        K = K
                        )

## Put forward and inverse together
# initialise a random model
m = transD_GP.init(opt)

# compute misfit
transD_GP.get_misfit(m, opt, FOp)

# if you've gotten till here you're golden
# you don't really need to initialize a model and compute
# misfit outside a module

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
using Pkg.Artifacts
using HTTP

# This is the path to the Artifacts.toml we will manipulate
artifact_toml = joinpath(@__DIR__, "Artifacts.toml")

# Query the `Artifacts.toml` file for the hash bound to the name "iris"
# (returns `nothing` if no such binding exists)
hitran_hash = artifact_hash("hitran", artifact_toml)

# If the name was not bound, or the hash it was bound to does not exist, create it!
if hitran_hash === nothing || !artifact_exists(hitran_hash)
    # create_artifact() returns the content-hash of the artifact directory once we're finished creating it
    hitran_hash = create_artifact() do artifact_dir
        # We create the artifact by simply downloading a few files into the new artifact directory
        hitran_url_base = "https://hitran.org/lbl/api?numin=0&numax=50000&iso_ids_list="
        for i = 1:130
            println(i)
            try
                download(hitran_url_base * string(i), joinpath(artifact_dir, "hitran_global_id_" * string(i) * ".par"))
            catch
            end
        end
    end

    # Now bind that hash within our `Artifacts.toml`.  `force = true` means that if it already exists,
    # just overwrite with the new content-hash.  Unless the source files change, we do not expect
    # the content hash to change, so this should not cause unnecessary version control churn.
    bind_artifact!(artifact_toml, "hitran", hitran_hash)
end

# Get the path of the iris dataset, either newly created or previously generated.
# this should be something like `~/.julia/artifacts/dbd04e28be047a54fbe9bf67e934be5b5e0d357a`
hitran_dataset_path = artifact_path(hitran_hash)


text:
function testone1()

    M0 = 1
    T1 = 1000
    T2 = 80
    Δω = 10 * 2π
    κ = 0.9
    Tfree = 8
    Tg = 3
    α = deg2rad(15)
    β = deg2rad(13)
    ϕ = deg2rad(30)
    TE = Tfree / 2
    s1 = stfr(M0, T1, T2, Δω, κ, Tfree, Tg, α, β, ϕ, TE)

    spin = Spin(M0, T1, T2, Δω / 2π)
    stfrblochsim! = STFRBlochSim(Tfree, Tg, Tfree / 2, κ*α, κ*β, ϕ)
    stfrblochsim!(spin)

    show(devnull, stfrblochsim!)
    show(devnull, "text/plain", stfrblochsim!)
    return s1 ≈ signal(spin)

end

function testone2()

    M0 = [1.04, 0.8, 1.3]
    T1 = [833, 900, 1100]
    T2 = [76, 90, 100]
    Δω = [10, 0, 30] * 2π
    κ = [0.9, 1, 1.1]
    Tfree = [8, 10]
    Tg = [3, 4]
    α = deg2rad.([15, 14])
    β = deg2rad.([13, 3])
    ϕ = deg2rad.([30, 0])
    TE = Tfree / 2
    s1 = stfr.(M0', T1', T2', Δω', κ', Tfree, Tg, α, β, ϕ, TE)

    s2 = Matrix{ComplexF64}(undef, length(Tfree), length(M0))
    workspace = STFRBlochSimWorkspace(Spin{Float64}, STFRBlochSim{<:Real,InstantaneousRF,InstantaneousRF,IdealSpoiling})
    for j = 1:length(M0), i = 1:length(Tfree)
        spin = Spin(M0[j], T1[j], T2[j], Δω[j] / 2π)
        stfrblochsim! = STFRBlochSim(Tfree[i], Tg[i], Tfree[i] / 2, κ[j] * α[i], κ[j] * β[i], ϕ[i])
        stfrblochsim!(spin, workspace)
        s2[i,j] = signal(spin)
    end

    return s1 ≈ s2

end

function testtwo1()

    M0 = 0.8
    ff = 0.15
    T1f = 400
    T1s = 1000
    T2f = 20
    T2s = 80
    Δωf = 17 * 2π
    Δω = 0 * 2π
    κ = 1
    Tfree = 10
    Tg = 5
    α = deg2rad(11)
    β = deg2rad(1)
    ϕ = deg2rad(0)
    TE = Tfree / 2
    s1 = stfr(M0, ff, T1f, T1s, T2f, T2s, Δωf, Δω, κ, Tfree, Tg, α, β, ϕ, TE)

    frac = (ff, 1-ff)
    T1 = (T1f, T1s)
    T2 = (T2f, T2s)
    Δω2 = (Δωf + Δω, Δω)
    τ = (Inf, Inf)
    spin = SpinMC(M0, frac, T1, T2, Δω2 ./ 2π, τ)
    stfrblochsim! = STFRBlochSim(Tfree, Tg, Tfree / 2, κ*α, κ*β, ϕ)
    stfrblochsim!(spin)

    return s1 ≈ signal(spin)

end

function testtwo2()

    M0 = [1.04, 0.8, 1.3]
    ff = [0.15, 0.10, 0.23]
    T1f = [400, 320, 833]
    T1s = [833, 900, 1100]
    T2f = [20,  16, 24]
    T2s = [76, 90, 100]
    Δωf = [17, 10, 3] * 2π
    Δω = [10, 0, 30] * 2π
    κ = [0.9, 1, 1.1]
    Tfree = [8, 10]
    Tg = [3, 4]
    α = deg2rad.([15, 14])
    β = deg2rad.([13, 3])
    ϕ = deg2rad.([30, 0])
    frac = permutedims([[ff[i], 1 .- ff[i]] for i = 1:length(ff)])
    T1 = permutedims([[T1f[i], T1s[i]] for i = 1:length(T1f)])
    T2 = permutedims([[T2f[i], T2s[i]] for i = 1:length(T2f)])
    Δω2 = permutedims([[Δωf[i], 0] .+ Δω[i] for i = 1:length(Δω)])
    τ = [[Inf, Inf]]
    TE = Tfree / 2
    s1 = stfr.(M0', ff', T1f', T1s', T2f', T2s', Δωf', Δω', κ', Tfree, Tg, α, β, ϕ, TE)

    s2 = Matrix{ComplexF64}(undef, length(Tfree), length(M0))
    workspace = STFRBlochSimWorkspace(SpinMC{Float64,2}, STFRBlochSim{<:Real,InstantaneousRF,InstantaneousRF,IdealSpoiling})
    for j = 1:length(M0), i = 1:length(Tfree)
        spin = SpinMC(M0[j], (ff[j], 1 - ff[j]), (T1f[j], T1s[j]), (T2f[j], T2s[j]), (Δωf[j] + Δω[j], Δω[j]) ./ 2π, (Inf, Inf))
        stfrblochsim! = STFRBlochSim(Tfree[i], Tg[i], Tfree[i] / 2, κ[j] * α[i], κ[j] * β[i], ϕ[i])
        stfrblochsim!(spin, workspace)
        s2[i,j] = signal(spin)
    end

    return s1 ≈ s2

end

function testnonideal1()

    M0 = [1.04, 0.8, 1.3]
    T1 = [833, 900, 1100]
    T2 = [76, 90, 100]
    Δω = [10, 0, 30] * 2π
    κ = [0.9, 1, 1.1]
    Tfree = [8, 10]
    Tg = [3, 4]
    α = deg2rad.([15, 14])
    β = deg2rad.([13, 3])
    ϕ = deg2rad.([30, 0])
    s1 = Matrix{ComplexF64}(undef, length(Tfree), length(M0))
    s2 = similar(s1)
    s3 = similar(s1)
    w1 = STFRBlochSimWorkspace(Spin{Float64}, STFRBlochSim{<:Real,InstantaneousRF,InstantaneousRF,GradientSpoiling})
    w2 = STFRBlochSimWorkspace(Spin{Float64}, STFRBlochSim{<:Real,InstantaneousRF,InstantaneousRF,RFSpoiling})
    w3 = STFRBlochSimWorkspace(Spin{Float64}, STFRBlochSim{<:Real,InstantaneousRF,InstantaneousRF,RFandGradientSpoiling})
    spoil1 = GradientSpoiling(0, 0, 0, 2)
    spoil2 = RFSpoiling(2π)
    spoil3 = RFandGradientSpoiling(spoil1, spoil2)
    nTR = Val(2000)
    for j = 1:length(M0), i = 1:length(Tfree)
        spin = Spin(M0[j], T1[j], T2[j], Δω[j] / 2π, Position(1, 1, 1))
        stfrblochsim! = STFRBlochSim(Tfree[i], Tg[i], Tfree[i] / 2, κ[j] * α[i], κ[j] * β[i], ϕ[i], spoil1)
        stfrblochsim!(spin, w1)
        s1[i,j] = signal(spin)
        spin.M.x = 0; spin.M.y = 0; spin.M.z = M0[j]
        stfrblochsim! = STFRBlochSim(Tfree[i], Tg[i], Tfree[i] / 2, κ[j] * α[i], κ[j] * β[i], ϕ[i], spoil2, nTR)
        stfrblochsim!(spin, w2)
        s2[i,j] = signal(spin)
        spin.M.x = 0; spin.M.y = 0; spin.M.z = M0[j]
        stfrblochsim! = STFRBlochSim(Tfree[i], Tg[i], Tfree[i] / 2, κ[j] * α[i], κ[j] * β[i], ϕ[i], spoil3, nTR, Val(true))
        stfrblochsim!(spin, w3)
        s3[i,j] = signal(spin)
    end

    return s1 ≈ s2 ≈ s3

end

function testnonideal2()

    M0 = [1.04, 0.8, 1.3]
    frac = [(0.1, 0.9), (0.15, 0.85), (0.2, 0.8)]
    T1 = [(300, 833), (400, 900), (800, 1100)]
    T2 = [(10, 76), (15, 90), (20, 100)]
    Δf = [(0, 0), (10, 0), (30, 10)]
    τ = [(20, 100), (30, 130), (70, 70)]
    κ = [0.9, 1, 1.1]
    Tfree = [8, 10]
    Tg = [3, 4]
    α = deg2rad.([15, 14])
    β = deg2rad.([13, 3])
    ϕ = deg2rad.([30, 0])
    s1 = Matrix{ComplexF64}(undef, length(Tfree), length(M0))
    s2 = similar(s1)
    s3 = similar(s1)
    w1 = STFRBlochSimWorkspace(SpinMC{Float64,2}, STFRBlochSim{<:Real,InstantaneousRF,InstantaneousRF,GradientSpoiling})
    w2 = STFRBlochSimWorkspace(SpinMC{Float64,2}, STFRBlochSim{<:Real,InstantaneousRF,InstantaneousRF,RFSpoiling})
    w3 = STFRBlochSimWorkspace(SpinMC{Float64,2}, STFRBlochSim{<:Real,InstantaneousRF,InstantaneousRF,RFandGradientSpoiling})
    spoil1 = GradientSpoiling(0, 0, 0, 2)
    spoil2 = RFSpoiling(2π)
    spoil3 = RFandGradientSpoiling(spoil1, spoil2)
    nTR = Val(2000)
    for j = 1:length(M0), i = 1:length(Tfree)
        spin = SpinMC(M0[j], frac[j], T1[j], T2[j], Δf[j], τ[j], Position(1, 1, 1))
        stfrblochsim! = STFRBlochSim(Tfree[i], Tg[i], Tfree[i] / 2, κ[j] * α[i], κ[j] * β[i], ϕ[i], spoil1)
        stfrblochsim!(spin, w1)
        s1[i,j] = signal(spin)
        for k = 1:2 spin.M[k].x = 0; spin.M[k].y = 0; spin.M[k].z = frac[j][k] * M0[j] end
        stfrblochsim! = STFRBlochSim(Tfree[i], Tg[i], Tfree[i] / 2, κ[j] * α[i], κ[j] * β[i], ϕ[i], spoil2, nTR)
        stfrblochsim!(spin, w2)
        s2[i,j] = signal(spin)
        for k = 1:2 spin.M[k].x = 0; spin.M[k].y = 0; spin.M[k].z = frac[j][k] * M0[j] end
        stfrblochsim! = STFRBlochSim(Tfree[i], Tg[i], Tfree[i] / 2, κ[j] * α[i], κ[j] * β[i], ϕ[i], spoil3, nTR, Val(true))
        stfrblochsim!(spin, w3)
        s3[i,j] = signal(spin)
    end

    return s1 ≈ s2 ≈ s3

end

function testnonideal3()

    Tfree = 10
    Tg = 3
    TE = Tfree / 2
    α = deg2rad(15)
    β = deg2rad(13)
    ϕ = 0
    (M0, T1, T2, Δf) = (1, 1000, 100, 0)
    spin = Spin(M0, T1, T2, Δf)
    stfrblochsim! = STFRBlochSim(Tfree, Tg, TE, α, β, ϕ)
    stfrblochsim!(spin)
    M1 = copy(spin.M)

    grad = Gradient(0, 0, 0.3)
    zmax = 2π / (GAMMA * grad.z * Tg / 1000)
    nz = 100
    z = (1:nz)/nz * zmax
    spoil = RFandGradientSpoiling(grad, Tg)
    stfrblochsim! = STFRBlochSim(Tfree, Tg, TE, α, β, ϕ, spoil, Val(2000))
    workspace = STFRBlochSimWorkspace(spin, stfrblochsim!)
    M2 = Magnetization(0.0, 0.0, 0.0)
    foreach(z) do z
        spin.M.x = 0; spin.M.y = 0; spin.M.z = M0; spin.pos.z = z
        stfrblochsim!(spin, workspace)
        add!(M2, spin.M)
    end
    div!(M2, nz)

    show(devnull, "text/plain", stfrblochsim!)
    return isapprox(M1, M2; atol = 1e-4)

end

function testnoninstant1()

    M0 = [1.04, 0.8, 1.3]
    ff = [0.15, 0.10, 0.23]
    T1f = [400, 320, 833]
    T1s = [833, 900, 1100]
    T2f = [20,  16, 24]
    T2s = [76, 90, 100]
    Δωf = [17, 10, 3] * 2π
    Δω = [10, 0, 30] * 2π
    κ = [0.9, 1, 1.1]
    Tfree = [8, 10]
    Tg = [3, 4]
    α = deg2rad.([15, 14])
    β = deg2rad.([13, 3])
    ϕ = deg2rad.([30, 0])
    frac = [(ff[i], 1 - ff[i]) for i = 1:length(ff)]
    T1 = collect(zip(T1f, T1s))
    T2 = collect(zip(T2f, T2s))
    Δf = [(Δωf[i] + Δω[i], Δω[i]) ./ 2π for i = 1:length(Δωf)]
    τ = (Inf, Inf)
    Δt = 0.0002
    TE = Tfree / 2
    s1 = stfr.(M0', ff', T1f', T1s', T2f', T2s', Δωf', Δω', κ', Tfree, Tg, α, β, ϕ, TE)

    s2 = Matrix{ComplexF64}(undef, length(Tfree), length(M0))
    workspace = STFRBlochSimWorkspace(SpinMC{Float64,2}, STFRBlochSim{<:Real,RF,RF,IdealSpoiling})
    for j = 1:length(M0), i = 1:length(Tfree)
        spin = SpinMC(M0[j], frac[j], T1[j], T2[j], Δf[j], τ)
        rftipdown = RF(fill(κ[j] * α[i] * 1000 / Δt / GAMMA / 3, 3), Δt)
        rftipup = RF(fill(κ[j] * -β[i] * 1000 / Δt / GAMMA / 3, 3), Δt, ϕ[i])
        stfrblochsim! = STFRBlochSim(Tfree[i], Tg[i], Tfree[i] / 2, rftipdown, rftipup)
        stfrblochsim!(spin, workspace)
        s2[i,j] = signal(spin)
    end

    return isapprox(s1, s2, atol = 1e-5)

end

function testnoninstant2()

    M0 = [1.04, 0.8, 1.3]
    frac = [(0.1, 0.9), (0.15, 0.85), (0.2, 0.8)]
    T1 = [(300, 833), (400, 900), (800, 1100)]
    T2 = [(10, 76), (15, 90), (20, 100)]
    Δf = [(0, 0), (10, 0), (30, 10)]
    τ = [(20, 100), (30, 130), (70, 70)]
    κ = [0.9, 1, 1.1]
    Tfree = [8, 10]
    Tg = [3, 4]
    α = deg2rad.([15, 14])
    β = deg2rad.([13, 3])
    ϕ = deg2rad.([30, 0])
    Δt = 0.5
    s1 = Matrix{ComplexF64}(undef, length(Tfree), length(M0))
    s2 = similar(s1)
    w1 = STFRBlochSimWorkspace(SpinMC{Float64,2}, STFRBlochSim{<:Real,RF,RF,GradientSpoiling})
    w2 = STFRBlochSimWorkspace(SpinMC{Float64,2}, STFRBlochSim{<:Real,RF,RF,RFandGradientSpoiling})
    spoil1 = GradientSpoiling(1, 2, 3, 1)
    spoil2 = RFandGradientSpoiling(spoil1, 2π)
    nTR = Val(2000)
    for j = 1:length(M0), i = 1:length(Tfree)
        spin = SpinMC(M0[j], frac[j], T1[j], T2[j], Δf[j], τ[j], Position(1, 1, 1))
        rftipdown = RF(fill(κ[j] * α[i] * 1000 / Δt / GAMMA / 3, 3), Δt)
        rftipup = RF(fill(κ[j] * -β[i] * 1000 / Δt / GAMMA / 3, 3), Δt, ϕ[i])
        stfrblochsim! = STFRBlochSim(Tfree[i], Tg[i], Tfree[i] / 2, rftipdown, rftipup, spoil1)
        stfrblochsim!(spin, w1)
        s1[i,j] = signal(spin)
        for k = 1:2 spin.M[k].x = 0; spin.M[k].y = 0; spin.M[k].z = frac[j][k] * M0[j] end
        stfrblochsim! = STFRBlochSim(Tfree[i], Tg[i], Tfree[i] / 2, rftipdown, rftipup, spoil2, nTR)
        stfrblochsim!(spin, w2)
        s2[i,j] = signal(spin)
    end

    return s1 ≈ s2

end

function testnoninstant3()

    M0 = 1
    T1 = 1000
    T2 = 100
    Δf = 0
    Tfree = 8
    Tg = 2.8
    TE = Tfree / 2
    α = deg2rad(15)
    β = deg2rad(15)
    ϕ = 0
    Δt = 0.001
    rftipdown = RF([α * 1000 / Δt / GAMMA], Δt)
    rftipup = RF([-β * 1000 / Δt / GAMMA], Δt)

    spin = Spin(M0, T1, T2, Δf)
    stfr! = STFRBlochSim(Tfree, Tg, TE, rftipdown, rftipup, Val(2000))
    stfr!(spin)
    s1 = spin.M

    spin = Spin(M0, T1, T2, Δf)
    stfr! = STFRBlochSim(Tfree, Tg, TE, α, β, ϕ, Val(2000))
    stfr!(spin)
    s2 = spin.M

    return s1 ≈ s2

end

@testset "Bloch Simulation" begin

    @testset "Single-Compartment" begin

        @test testone1()
        @test testone2()

    end

    @testset "Two-Compartment" begin

        @test testtwo1()
        @test testtwo2()

    end

    @testset "Nonideal Spoiling" begin

        @test testnonideal1()
        @test testnonideal2()
        @test testnonideal3()

    end

    @testset "Non-instantaneous Excitation" begin

        @test testnoninstant1()
        @test testnoninstant2()
        @test testnoninstant3()

    end

end


text:
@testset "Solutions" begin
  @testset "Estimation" begin
    d = PointSet(rand(2,100))
    m = Dict(:z => 1:100)
    v = Dict(:z => 1:100)
    s = EstimationSolution(d, m, v)
    @test s[:z] == georef((z=1:100, zvar=1:100), d)

    @test sprint(show, s) == "2D EstimationSolution"
    @test sprint(show, MIME"text/plain"(), s) == "2D EstimationSolution\n  domain: 100 PointSet{Float64,2}\n  variables: z"

    d = RegularGrid(10,10)
    s = EstimationSolution(d, m, v)
    @test s[:z] == georef((z=1:100, zvar=1:100), d)

    @test sprint(show, s) == "2D EstimationSolution"
    @test sprint(show, MIME"text/plain"(), s) == "2D EstimationSolution\n  domain: 10×10 RegularGrid{Float64,2}\n  variables: z"

    if visualtests
      @plottest plot(s,size=(800,400)) joinpath(datadir,"estimsol.png") !istravis
    end
  end

  @testset "Simulation" begin
    d = PointSet(rand(2,100))
    r = Dict(:z => [1:100 for i in 1:10])
    s = SimulationSolution(d, r)
    @test s[:z] == r[:z]

    @test sprint(show, s) == "2D SimulationSolution"
    @test sprint(show, MIME"text/plain"(), s) == "2D SimulationSolution\n  domain: 100 PointSet{Float64,2}\n  variables: z\n  N° reals:  10"

    d = RegularGrid(10,10)
    s = SimulationSolution(d, r)

    @test sprint(show, s) == "2D SimulationSolution"
    @test sprint(show, MIME"text/plain"(), s) == "2D SimulationSolution\n  domain: 10×10 RegularGrid{Float64,2}\n  variables: z\n  N° reals:  10"

    if visualtests
      @plottest plot(s,size=(800,300)) joinpath(datadir,"simsol.png") !istravis
    end
  end
end


text:
# Daily numbers across C: Confirmed, R:Recovered , D:Deceased and Tested per state (historical data).
timeseries_url = "https://api.covid19india.org/v4/min/timeseries.min.json"

# Current day numbers across districts and states.
data_today_url = "https://api.covid19india.org/v4/min/data.min.json"

# Per day numbers across districts and states - consider using timeseries in place of this.
# This is a huge file and is a mix of timeseries and data.min.json
data_all_url = "https://api.covid19india.org/v4/min/data-all.min.json"

# Aggregated sheets (CSV files)

# 1. India level timeseries for Confirmed, Recovered and Deceased cases [Reading]
case_time_series_csv_url = "https://api.covid19india.org/csv/latest/case_time_series.csv"

# 2. Statewise timeseries of Confirmed, Recovered and Deceased numbers. [Reading]
states_csv_url = "https://api.covid19india.org/csv/latest/states.csv"

# 3. Districtwise timeseries of Confirmed, Recovered and Deceased numbers. [Reading]
districts_csv_url = "https://api.covid19india.org/csv/latest/districts.csv"

# 4. Statewise per day delta of Confirmed, Recovered and Deceased numbers. [Reading]
state_wise_daily_csv_url = "https://api.covid19india.org/csv/latest/state_wise_daily.csv"

# 5. Statewise cumulative numbers till date. [Reading]
state_wise_csv_url = "https://api.covid19india.org/csv/latest/state_wise.csv"

# 6. Districtwise Cumulative numbers till date. [Reading]
district_wise_csv_url =	"https://api.covid19india.org/csv/latest/district_wise.csv"

# 7. Number of tests conducted in each state, ventilators ,hospital bed occupany reported in state bulletins
statewise_tested_numbers_data_csv_url = "https://api.covid19india.org/csv/latest/statewise_tested_numbers_data.csv"

# 8. Number of tests reported by ICMR
tested_numbers_icmr_data_csv_url = "https://api.covid19india.org/csv/latest/tested_numbers_icmr_data.csv"

# 9. Number of Labs in each state as per ICMR
icmr_labs_statewise_csv_url = "https://api.covid19india.org/csv/latest/icmr_labs_statewise.csv"

# 10. List of sources that we are using.
sources_list_csv_url = "https://api.covid19india.org/csv/latest/sources_list.csv"

# 11. Number of RTPCR samples collected statewise in ICMR Application
rtpcr_samples_collected_csv_url	= "http://api.covid19india.org/csv/latest/icmr_rtpcr_tests_daily.csv"

# 12. Number of vaccine doses administered statewise [Reading]
vaccine_doses_administered_statewise_csv_url = "http://api.covid19india.org/csv/latest/vaccine_doses_statewise.csv"

# 13. Key data points from CoWin database at a state level
cowin_vaccine_data_statewise_csv_url = "http://api.covid19india.org/csv/latest/cowin_vaccine_data_statewise.csv"

# 14. Key data points from CoWin database at a district level
cowin_vaccine_data_districtwise_csv_url = "http://api.covid19india.org/csv/latest/cowin_vaccine_data_districtwise.csv"


text:
using bayesian_nlde
using DataFrames
using Plots
using Turing
# using ReverseDiff
# using DiffEqSensitivity

Turing.setadbackend(:forwarddiff)

# load data and format for inference
dat = import_data(download_data = false)
grouped_dat = [DataFrame(df) for df in groupby(dat, :Country_Region)];
grouped_sir_data = [Matrix((df[:, [:susc_prop, :case_prop, :recov_prop]]))' for df in groupby(dat, :Country_Region)];
country_pops = [df.pop[1] for df in grouped_dat]
country_names = [df.Country_Region[1] for df in grouped_dat]

# look at countries with more than a months worth of reporting
sub_idx = [i for (i,x) in enumerate(grouped_sir_data) if size(x, 2) > 28]
sub_idx = sub_idx[1:10]
sub_grouped_sir_data = grouped_sir_data[sub_idx]
sub_country_pops = country_pops[sub_idx]

model = hierarchical_sir_model(grouped_sir = sub_grouped_sir_data,
    m_countries = length(sub_grouped_sir_data),
    country_pops = sub_country_pops)

varinfo = Turing.VarInfo(hierarchical_sir_model)
spl = Turing.SampleFromPrior()
@code_warntype model.f(varinfo, spl, Turing.DefaultContext(), model)

# infer paramters
@time chains = sample(hierarchical_sir_model(grouped_sir = sub_grouped_sir_data,
    m_countries = length(sub_grouped_sir_data), country_pops = sub_country_pops),
    NUTS(100, 0.8), 500; progress = true, verbose = true)
write(string(dirname(@__DIR__)) * "/data/chains/posterior_samples.jls", chains)

#plotlyjs()

# plot distribution over trajectories
for i in sub_idx
    plt = plot_country_posterior(chains, grouped_dat, i; subsample_size = 200)
    plt_path = string(dirname(@__DIR__)) * "/data/plots/" * string(country_names[i]) * ".svg"
    savefig(plt, plt_path)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test

function get_test_model()
    X = Dict(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 0.5],
        4 => [0.0, 0.5],
        5 => [0.0, 0.6],
        6 => [1.0, 0.6],
        7 => [1.0, 1.1],
        8 => [0.0, 1.1])
    el1 = Element(Quad4, (1, 2, 3, 4))
    el2 = Element(Quad4, (5, 6, 7, 8))
    el3 = Element(Seg2, (1, 2))
    el4 = Element(Seg2, (7, 8))
    el5 = Element(Seg2, (4, 3))
    el6 = Element(Seg2, (5, 6))
    update!((el1, el2, el3, el4, el5, el6), "geometry", X)
    update!((el1, el2), "youngs modulus", 96.0)
    update!((el1, el2), "poissons ratio", 1/3)
    update!(el3, "displacement 1", 0.0)
    update!(el3, "displacement 2", 0.0)
    update!(el4, "displacement 1", 0.0)
    update!(el4, "displacement 2", 0.0)
    update!(el6, "master elements", [el5])
    p1 = Problem(Elasticity, "body1", 2)
    p2 = Problem(Elasticity, "body2", 2)
    p3 = Problem(Dirichlet, "fixed", 2, "displacement")
    p4 = Problem(Mortar2D, "interface", 2, "displacement")
    add_elements!(p1, el1)
    add_elements!(p2, el2)
    add_elements!(p3, el3, el4)
    add_elements!(p4, el5, el6)
    return p1, p2, p3, p4
end

## test adjust setting in 2d tie contact
p1, p2, p3, p4 = get_test_model()
p1.properties.formulation = :plane_stress
p2.properties.formulation = :plane_stress
#p4.properties.adjust = true
#p4.properties.rotate_normals = false
analysis = Analysis(Linear)
add_problems!(analysis, p1, p2, p3, p4)
run!(analysis)
el5 = first(get_elements(p4))
xi, time = (0.0, ), 0.0
u = el5("displacement", xi, time)
@debug("u = $u")
@test_broken isapprox(u, [0.0, 0.05])

## test that interface transfers constant field without error
meshfile = @__DIR__() * "/testdata/block_2d.med"
mesh = aster_read_mesh(meshfile)

upper = Problem(PlaneHeat, "upper", 1)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "thermal conductivity", 1.0)
add_elements!(upper, upper_elements)

lower = Problem(PlaneHeat, "lower", 1)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "thermal conductivity", 1.0)
add_elements!(lower, lower_elements)

bc_upper = Problem(Dirichlet, "upper boundary", 1, "temperature")
bc_upper_elements = create_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "temperature 1", 0.0)
add_elements!(bc_upper, bc_upper_elements)

bc_lower = Problem(Dirichlet, "lower boundary", 1, "temperature")
bc_lower_elements = create_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "temperature 1", 1.0)
add_elements!(bc_lower, bc_lower_elements)

interface = Problem(Mortar2D, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_elements(mesh, "LOWER_TOP")
interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
add_master_elements!(interface, interface_master_elements)
add_slave_elements!(interface, interface_slave_elements)

analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, interface)
run!(analysis)

#interface_norm = norm(interface.assembly)
# for bi-orthogonal:
#interface_norm_expected = [0.0, 0.0, 0.0, 0.0, 0.0, 0.44870723441585775, 0.44870723441585775, 0.0, 0.0, 0.0]
#interface_norm_expected = [0.0, 0.0, 0.0, 0.0, 0.0, 0.39361633468943247, 0.39361633468943247, 0.0, 0.0, 0.0]
#@info("Interface norm:          $interface_norm")
#@info("Interface norm expected: $interface_norm_expected")
#@test isapprox(interface_norm, interface_norm_expected)

xi, time = (0.0, ), 0.0
T_upper = first(bc_upper_elements)("temperature", xi, time)
T_lower = first(bc_lower_elements)("temperature", xi, time)
T_middle = first(interface_slave_elements)("temperature", xi, time)
@debug("T upper: $T_upper, T lower: $T_lower, T interface: $T_middle")

node_ids, temperature = get_nodal_vector(interface_slave_elements, "temperature", 0.0)
T = [t[1] for t in temperature]
minT = minimum(T)
maxT = maximum(T)
@debug("minT = $minT, maxT = $maxT")
@test isapprox(minT, 0.5)
@test isapprox(maxT, 0.5)

## test mesh tie with splitted block and plane stress elasticity

meshfile = @__DIR__() * "/testdata/block_2d.med"
mesh = aster_read_mesh(meshfile)

upper = Problem(Elasticity, "upper", 2)
upper.properties.formulation = :plane_stress
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "youngs modulus", 100.0)
update!(upper_elements, "poissons ratio", 1/3)
add_elements!(upper, upper_elements)

lower = Problem(Elasticity, "lower", 2)
lower.properties.formulation = :plane_stress
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "youngs modulus", 100.0)
update!(lower_elements, "poissons ratio", 1/3)
add_elements!(lower, lower_elements)

bc_upper = Problem(Dirichlet, "upper boundary", 2, "displacement")
bc_upper_elements = create_elements(mesh, "UPPER_TOP")
#   update!(bc_upper.elements, "displacement 1", 0.1)
update!(bc_upper_elements, "displacement 2", -0.1)
add_elements!(bc_upper, bc_upper_elements)

bc_lower = Problem(Dirichlet, "lower boundary", 2, "displacement")
bc_lower_elements = create_elements(mesh, "LOWER_BOTTOM")
#   update!(bc_lower.elements, "displacement 1", 0.0)
update!(bc_lower_elements, "displacement 2", 0.0)
add_elements!(bc_lower, bc_lower_elements)

bc_corner = Problem(Dirichlet, "fix model from lower left corner to prevent singularity", 2, "displacement")
node_ids = find_nearest_nodes(mesh, [0.0, 0.0])
bc_corner_elements = [Element(Poi1, node_ids)]
update!(bc_corner_elements, "geometry", mesh.nodes)
update!(bc_corner_elements, "displacement 1", 0.0)
add_elements!(bc_corner, bc_corner_elements)

interface = Problem(Mortar2D, "interface between upper and lower block", 2, "displacement")
interface_slave_elements = create_elements(mesh, "LOWER_TOP")
interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
add_master_elements!(interface, interface_master_elements)
add_slave_elements!(interface, interface_slave_elements)

analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, interface, bc_corner)
run!(analysis)

slave_elements = get_slave_elements(interface)
node_ids, la = get_nodal_vector(slave_elements, "lambda", 0.0)
for lai in la
    @test isapprox(lai, [0.0, 10.0])
end


text:
export Role

mutable struct Role <: AbstractModel
  _table_name::String
  _id::String

  id::Nullable{SearchLight.DbId}
  name::Symbol

  has_one::Vector{SearchLight.SQLRelation}
  on_dehydration::Function

  Role(;
    id = Nullable{SearchLight.DbId}(),
    name = :user,
    has_one = [SearchLight.SQLRelation(User)],
    on_dehydration = Roles.dehydrate
  ) = new("roles", "id", id, name, has_one, on_dehydration)
end

module Roles

using App

function dehydrate(r::Role, field::Symbol, value::Any)
  return  if field == :name
            value = string(value)
          end
end

end


text:
module CUTENSOR

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, cudaDataType
using ..CUDA: libcutensor,  @retry_reclaim

using CEnum

using Memoize

using DataStructures

const cudaDataType_t = cudaDataType

# core library
include("libcutensor_common.jl")
include("error.jl")
include("libcutensor.jl")

# low-level wrappers
include("tensor.jl")
include("wrappers.jl")

# high-level integrations
include("interfaces.jl")

# thread cache for task-local library handles
const thread_handles = Vector{Union{Nothing,Base.RefValue{cutensorHandle_t}}}()

# cache for created, but unused handles
const old_handles = DefaultDict{CuContext,Vector{Base.RefValue{cutensorHandle_t}}}(()->Base.RefValue{cutensorHandle_t}[])

function handle()
    CUDA.detect_state_changes()
    tid = Threads.threadid()
    if @inbounds thread_handles[tid] === nothing
        ctx = context()
        thread_handles[tid] = get!(task_local_storage(), (:CUTENSOR, ctx)) do
            if isempty(old_handles[ctx])
                handle = Ref{cutensorHandle_t}()
                cutensorInit(handle)
            else
                handle = pop!(old_handles[ctx])
            end

            finalizer(current_task()) do task
                push!(old_handles[ctx], handle)
            end
            # TODO: destroy to preserve memory, or at exit?

            handle
        end
    end
    something(@inbounds thread_handles[tid])
end

@inline function set_stream(stream::CuStream)
    # CUTENSOR uses stream arguments per operation
    return
end

function __init__()
    resize!(thread_handles, Threads.nthreads())
    fill!(thread_handles, nothing)

    CUDA.atdeviceswitch() do
        tid = Threads.threadid()
        thread_handles[tid] = nothing
    end

    CUDA.attaskswitch() do
        tid = Threads.threadid()
        thread_handles[tid] = nothing
    end
end

end


text:
# Code to test whether Outlook is running on Windows.   This is used
# to work around a bug in the Miniconda installer, which crashes Outlook
# if it is running.  See Luthaf/Conda.jl#15

immutable PROCESSENTRY32
    dwSize::Int32
    cntUsage::Int32
    th32ProcessID::Int32
    th32DefaultHeapID::UInt
    th32ModuleID::Int32
    cntThreads::Int32
    th32ParentProcessID::Int32
    pcPriClassBase::Int32
    dwFlags::UInt32
    szExeFile::NTuple{260,UInt8}
    PROCESSENTRY32() = new(sizeof(PROCESSENTRY32),0,0,0,0,0,0,0,0)
end
szExeFile(p::PROCESSENTRY32) = unsafe_string(pointer(collect(p.szExeFile)))
const TH32CS_SNAPPROCESS = 0x00000002
function isrunning(exefile::AbstractString)
    snapshot = ccall(:CreateToolhelp32Snapshot, stdcall, Ptr{Void}, (UInt32,Int32), TH32CS_SNAPPROCESS, 0)
    try
        entry = Ref(PROCESSENTRY32())
        if ccall(:Process32First, stdcall, Cint, (Ptr{Void}, Ref{PROCESSENTRY32}), snapshot, entry) == 1
            szExeFile(entry[]) == exefile && return true
            while ccall(:Process32Next, stdcall, Cint, (Ptr{Void}, Ref{PROCESSENTRY32}), snapshot, entry) == 1
                szExeFile(entry[]) == exefile && return true
            end
        end
    finally
        ccall(:CloseHandle, stdcall, Cint, (Ptr{Void},), snapshot)
    end
    return false
end
is_outlook_running() = isrunning("outlook.exe")


text:
#ENO coefficients for uniform mesh
#Based on: Shu C. Essentially Non-Oscillatory and Weighted Essentially
#Non-Oscillatory Schemes for Hyperbolic Conservation Laws
function unif_crj(k::Int)
  if k == 1
      return([1;1])
  end
  crj = zeros(k+1,k)
  for i = 1:(k+1)
    r = i-2
    for j = 1:k
      psum = 0.0
      for m = j:k
        numer = 0.0
        for l = 0:k
          if (l!=m)
            numer_prod = 1.0
            for q = 0:k
              if (q != m && q != l)
                numer_prod = numer_prod*(r-q+1)
              end
            end
            numer = numer + numer_prod
          end
        end
        denom = 1.0
        for l = 0:k
          if (l!=m);denom = denom * (m-l);end
        end
        psum = psum + numer/denom
      end
      crj[i,j] = psum
    end
  end
  crj
end

# Reconstruction types
struct ENO_Reconstruction{kType} <: AbstractReconstruction
    order::kType
    crj::Matrix{Float64}
end
ENO_Reconstruction(order::Int) = ENO_Reconstruction(order, unif_crj(Int((order + 1)/2)))

struct WENO_Reconstruction{kType,eType} <: AbstractReconstruction
    order::kType
    crj::Matrix{Float64}
    ɛ::eType
end
WENO_Reconstruction(order::Int; ɛ = 1e-6) = WENO_Reconstruction(order, unif_crj(Int((order + 1)/2)), ɛ)

struct MWENO_Reconstruction{kType,eType} <: AbstractReconstruction
    order::kType
    crj::Matrix{Float64}
    ɛ::eType
end
MWENO_Reconstruction(order::Int; ɛ = 1e-6) = MWENO_Reconstruction(order, unif_crj(Int((order + 1)/2)), ɛ)

#Eno reconstruction for uniform mesh
function reconstruct(vloc::AbstractVector, dx::Real, eno_rec::ENO_Reconstruction)
  order = eno_rec.order
  k = Int((order + 1)/2)
  crj = eno_rec.crj
  vdiffs = Vector{typeof(vloc)}(undef,0)
  vl = zero(eltype(vloc))
  vr = zero(eltype(vloc))
  N = size(vloc,1)
  if (N != 2*k-1)
    throw("dimension of vloc is not consistent with order $k ENO")
  end
  # Calculation of divided differences
  push!(vdiffs,copy(vloc))
  for i in 2:k
    push!(vdiffs,(vdiffs[i-1][2:end]-vdiffs[i-1][1:end-1])/dx)
  end
  #Calculation of Stencil
  r = 0
  if k < 2;return(vloc[1],vloc[1]);end
  for m = 2:k
    if (abs(vdiffs[m][k-r-1])<abs(vdiffs[m][k-r]))
      r = r + 1
    end
  end
  for j = 0:(k-1)
    vl = vl + vloc[k-r+j]*crj[r+1,j+1]
    vr = vr + vloc[k-r+j]*crj[r+2,j+1]
  end
  return(vl,vr)
end

#Eno reconstruction for unestructured mesh
function reconstruct(vloc::AbstractVector, xloc::AbstractVector, eno_rec::ENO_Reconstruction)
  order = eno_rec.order
  k = Int((order + 1)/2)
  crj = eno_rec.crj
  vdiffs = Vector{typeof(vloc)}()
  vl = zero(eltype(vloc))
  vr = zero(eltype(vloc))
  N = size(vloc,1)
  if (N != 2*k-1)
    throw("dimension of vloc is not consistent with order $k ENO")
  end
  # Calculation of divided differences
  push!(vdiffs,copy(vloc))
  for i in 2:k
    push!(vdiffs,(vdiffs[i-1][2:end]-vdiffs[i-1][1:end-1])./(xloc[i:N]-xloc[1:N-i+1]))
  end
  #Calculation of Stencil
  r = 0
  if k < 2;return(vloc[1],vloc[1]);end
  for m = 2:k
    if (abs(vdiffs[m][k-r-1])<abs(vdiffs[m][k-r]))
      r = r + 1
    end
  end
  for j = 0:(k-1)
    vl = vl + vloc[k-r+j]*crj[r+1,j+1]
    vr = vr + vloc[k-r+j]*crj[r+2,j+1]
  end
  return(vl,vr)
end

#WEno reconstruction for uniform mesh
#Order available: 1, 3, 5
# Jiang, Shu, Efficient Implementation of Weighted ENO Schemes
@inline function get_βk(vloc::Vector, k)
  βk = zeros(k); dr = zeros(k);
  if (k==2)
      dr = [2//3,1//3]
      βk = [(vloc[3]-vloc[2])^2, (vloc[2]-vloc[1])^2]
  elseif (k==3)
      dr = [3//10, 3//5, 1//10]
      βk = [13//12*(vloc[3]-2*vloc[4]+vloc[5])^2 + 1//4*(3*vloc[3]-4*vloc[4]+vloc[5])^2,
      13//12*(vloc[2]-2*vloc[3]+vloc[4])^2 + 1//4*(vloc[2]-vloc[4])^2,
      13//12*(vloc[1]-2*vloc[2]+vloc[3])^2 + 1//4*(3*vloc[3]-4*vloc[2]+vloc[1])^2]
  else
    throw("WENO reconstruction of order $order is not implemented yet!")
  end
  return βk, dr
end

function reconstruct(vloc::AbstractVector, rec_scheme::WENO_Reconstruction)
  order = rec_scheme.order
  k = Int((order + 1)/2)
  crj = rec_scheme.crj
  ɛ = rech_scheme.ɛ
  vl = zero(eltype(vloc))
  vr = zero(eltype(vloc))
  N = size(vloc,1)
  if (N != order)
    throw("dimension of vloc is not consistent with order $order WENO")
  end

  #Special case k = 1
  if (k == 1)
    vl = vloc[1]; vr = vloc[1]
    return vl, vr
  end

  # Apply WENO procedure
  αl = zeros(k); αr = zeros(k);
  ωl = zeros(k); ωr = zeros(k);
  βk = zeros(k); dr = zeros(k);

  # Compute k values of xl and xr based on different stencils
  ulr = zeros(k); urr = zeros(k);
  for r=0:(k-1)
      for i=0:k-1
          urr[r+1] = urr[r+1] + crj[r+2,i+1]*vloc[k-r+i];
          ulr[r+1] = ulr[r+1] + crj[r+1,i+1]*vloc[k-r+i];
      end
  end

  # Set up WENO coefficients for different orders - 2k-1
  βk, dr = get_βk(vloc, k)

  # Compute α parameters
  for r=1:k
      αr[r] = dr[r]/(ɛ+βk[r])^2;
      αl[r] = dr[k+1-r]/(ɛ+βk[r])^2;
  end

  # Compute wENO weights parameters
  for r=1:k
      ωl[r] = αl[r]/sum(αl);
      ωr[r] = αr[r]/sum(αr);
  end

  # Compute cell interface values
  for r=1:k
      vl = vl + ωl[r]*ulr[r];
      vr = vr + ωr[r]*urr[r];
  end
  return(vl,vr)
end

function reconstruct(vmloc::AbstractVector, vploc::AbstractVector,rec_scheme::WENO_Reconstruction)
  order = rec_scheme.order
  k = Int((order + 1)/2)
  crj = rec_scheme.crj
  ɛ = rec_scheme.ɛ
  vl = zero(eltype(vmloc))
  vr = zero(eltype(vploc))
  N = size(vmloc,1)
  if (N != order)
    throw("dimension of vloc is not consistent with order $order WENO")
  end
  #Special case k = 1
  if (k == 1)
    vl = vmloc[1]; vr = vploc[1]
    return vl, vr
  end

  # Apply WENO procedure
  αl = zeros(k); αr = zeros(k);
  ωl = zeros(k); ωr = zeros(k);
  βk = zeros(k); dr = zeros(k);

  # Compute k values of xl and xr based on different stencils
  ulr = zeros(k); urr = zeros(k);
  for r=0:(k-1)
      for i=0:k-1
          urr[r+1] = urr[r+1] + crj[r+2,i+1]*vploc[k-r+i];
          ulr[r+1] = ulr[r+1] + crj[r+1,i+1]*vmloc[k-r+i];
      end
  end

  # Set up WENO coefficients for different orders - 2k-1
  βrk, dr = get_βk(vploc, k)
  βlk, dr = get_βk(vmloc, k)

  # Compute α parameters
  for r=1:k
      αr[r] = dr[r]/(ɛ+βrk[r])^2;
      αl[r] = dr[k+1-r]/(ɛ+βlk[r])^2;
  end

  # Compute wENO weights parameters
  for r=1:k
      ωl[r] = αl[r]/sum(αl);
      ωr[r] = αr[r]/sum(αr);
  end

  # Compute cell interface values
  for r=1:k
      vl = vl + ωl[r]*ulr[r];
      vr = vr + ωr[r]*urr[r];
  end
  return(vl,vr)
end

#Mapped WEno reconstruction for uniform mesh
#Reference:
# A. Henrick, T. Aslam, J. Powers, Mapped weighted essentially non-oscillatory
# schemes: Achiving optimal order near critical points

#Mapping function
@inline function gk(ω::Vector, dr::Vector)
  g = zero(ω)
  for i in 1:size(ω,1)
    g[i] = ω[i]*(dr[i]+dr[i]^2-3*dr[i]*ω[i]+ω[i]^2)/(dr[i]^2+ω[i]*(1-2*dr[i]))
  end
  g
end

function reconstruct(vloc::AbstractVector, rec_scheme::MWENO_Reconstruction)
  order = rec_scheme.order
  k = Int((order + 1)/2)
  crj = rec_scheme.crj
  ɛ = rech_scheme.ɛ
  vl = zero(eltype(vloc))
  vr = zero(eltype(vloc))
  N = size(vloc,1)
  if (N != order)
    throw("dimension of vloc is not consistent with order $order WENO")
  end
  #Special case k = 1
  if (k == 1)
    vl = vloc[1]; vr = vloc[1]
    return vl, vr
  end

  # Apply WENO procedure
  αl = zeros(k); αr = zeros(k);
  ωl = zeros(k); ωr = zeros(k);
  αml = zeros(k); αmr = zeros(k);
  ωml = zeros(k); ωmr = zeros(k);
  βk = zeros(k); dr = zeros(k);

  # Compute k values of xl and xr based on different stencils
  ulr = zeros(k); urr = zeros(k);
  for r=0:(k-1)
      for i=0:k-1
          urr[r+1] = urr[r+1] + crj[r+2,i+1]*vloc[k-r+i];
          ulr[r+1] = ulr[r+1] + crj[r+1,i+1]*vloc[k-r+i];
      end
  end

  # Set up WENO coefficients for different orders - 2k-1
  if (k==2)
      dr = [2//3,1//3]
      βk = [(vloc[3]-vloc[2])^2, (vloc[2]-vloc[1])^2]
  elseif (k==3)
      dr = [3//10, 3//5, 1//10]
      βk = [13//12*(vloc[3]-2*vloc[4]+vloc[5])^2 + 1//4*(3*vloc[3]-4*vloc[4]+vloc[5])^2,
      13//12*(vloc[2]-2*vloc[3]+vloc[4])^2 + 1//4*(vloc[2]-vloc[4])^2,
      13//12*(vloc[1]-2*vloc[2]+vloc[3])^2 + 1//4*(3*vloc[3]-4*vloc[2]+vloc[1])^2]
  else
    throw("WENO reconstruction of order $order is not implemented yet!")
  end

  # Compute α parameters
  for r=1:k
      αr[r] = dr[r]/(ɛ+βk[r])^2;
      αl[r] = dr[k+1-r]/(ɛ+βk[r])^2;
  end

  # Compute wENO weights parameters
  for r=1:k
      ωl[r] = αl[r]/sum(αl);
      ωr[r] = αr[r]/sum(αr);
  end

  # Compute α mapped parameters
  αmr = gk(ωr)
  αml = gk(ωl)

  # Compute mapped wENO weights parameters
  for r=1:k
      ωml[r] = αml[r]/sum(αml);
      ωmr[r] = αmr[r]/sum(αmr);
  end

  # Compute cell interface values
  for r=1:k
      vl = vl + ωml[r]*ulr[r];
      vr = vr + ωmr[r]*urr[r];
  end
  return(vl,vr)
end

function reconstruct(vmloc::AbstractVector, vploc::AbstractVector,rec_scheme::MWENO_Reconstruction)
  order = rec_scheme.order
  k = Int((order + 1)/2)
  crj = rec_scheme.crj
  ɛ = rec_scheme.ɛ
  vl = zero(eltype(vmloc))
  vr = zero(eltype(vploc))
  N = size(vmloc,1)
  if (N != order)
    throw("dimension of vloc is not consistent with order $order WENO")
  end
  #Special case k = 1
  if (k == 1)
    vl = vmloc[1]; vr = vploc[1]
    return vl, vr
  end

  # Apply WENO procedure
  αl = zeros(k); αr = zeros(k);
  ωl = zeros(k); ωr = zeros(k);
  αml = zeros(k); αmr = zeros(k);
  ωml = zeros(k); ωmr = zeros(k);
  βk = zeros(k); dr = zeros(k);

  # Compute k values of xl and xr based on different stencils
  ulr = zeros(k); urr = zeros(k);
  for r=0:(k-1)
      for i=0:k-1
          urr[r+1] = urr[r+1] + crj[r+2,i+1]*vploc[k-r+i];
          ulr[r+1] = ulr[r+1] + crj[r+1,i+1]*vmloc[k-r+i];
      end
  end

  # Set up WENO coefficients for different orders - 2k-1
  βrk, dr = get_βk(vploc, k)
  βlk, dr = get_βk(vmloc, k)

  # Compute α parameters
  for r=1:k
      αr[r] = dr[r]/(ɛ+βrk[r])^2;
      αl[r] = dr[k+1-r]/(ɛ+βlk[r])^2;
  end

  # Compute wENO weights parameters
  for r=1:k
      ωl[r] = αl[r]/sum(αl);
      ωr[r] = αr[r]/sum(αr);
  end

  # Compute α mapped parameters
  drl = reverse(dr)
  αmr = gk(ωr,dr)
  αml = gk(ωl,drl)

  # Compute mapped wENO weights parameters
  for r=1:k
      ωml[r] = αml[r]/sum(αml);
      ωmr[r] = αmr[r]/sum(αmr);
  end

  # Compute cell interface values
  for r=1:k
      vl = vl + ωml[r]*ulr[r];
      vr = vr + ωmr[r]*urr[r];
  end
  return(vl,vr)
end


text:
include("../Happy-Number/HN.jl")

function addDigit(n::Int) ::Int
    n<0&& return -1
    n in range(0,stop=9) && return n
    n >= 10 && addDigit(sum(eachBit(n)))
end

max_stars_repo_name: ccqpein/Arithmetic-Exercises

max_issues_repo_name: ccqpein/Arithmetic-Exercises

max_forks_repo_name: ccqpein/Arithmetic-Exercises


text:
# Copyright (C) 2021 Ilya Baran.  This program is distributed under the terms of the MIT license.

using ModelBlocks
using Test
using OrdinaryDiffEq
using Plots

println("Bindings Test starting");

variables = Variables([
    Variable("X", 1, (0, 100), "", ""),
    Variable("Y", 1, (0, 100), "", ""),
]);

parameters = Variables([
    Variable("p1", 1.0, (0, 100), "", ""),
    Variable("p2", 2.0, (0, 100), "", ""),
    Variable("p3", 3.0, (0, 100), "", ""),
    Variable("p4", 4.0, (0, 100), "", ""),
]);

reactions = [
    GeneralReaction("Really random stuff", (dvdt, t, v, p, e) -> begin
        dvdt.X += p.p1 * v.Y;
        dvdt.Y -= p.p2 * v.X;
        @test p.p1 == 3;
        @test p.p3 == t / 5.;
    end),
];

block = Block(variables, parameters, reactions);
setTimeRange!(block, 0:100);

boundBlock = BlockWithBindings(block, Dict([
    ("p1", 3),
    ("p3", (t, p, v) -> t / 5.)
]));

@test getVariables(boundBlock).values == getVariables(block).values;
@test getParameters(boundBlock).values == [2.0, 4.0];

@time solution = runBlock(boundBlock);

boundBlock = BlockWithBindings(block, Dict([
    ("p2", 3),
]));

@test getParameters(boundBlock).values == [1.0, 3.0, 4.0];

println("Bindings Test finished");

max_stars_repo_name: ilya-baran-personal/modelBlocks

max_issues_repo_name: ilya-baran-personal/modelBlocks

max_forks_repo_name: ilya-baran-personal/modelBlocks


text:
#=
# 3D Poisson, Dirichlet bc
# Uses Dendro imported as a custom gen target
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("poissondendro");

useLog("poissondendrolog", level=3)

# default values (max_depth=6, wavelet_tol = 0.1, partition_tol = 0.3, solve_tol = 1e-6, max_iters = 100)
generateFor(DENDRO, params=(7, 0.01, 0.3, 0.000001, 1000))

domain(3)
functionSpace(order=4)

u = variable("u")
testSymbol("v")

boundary(u, 1, DIRICHLET, "0")

f = coefficient("f", "-30*pi*pi*sin(5*pi*x)*sin(2*pi*y)*sin(pi*z)")
weakForm(u, "-dot(grad(u),grad(v)) - f*v")

build_octree_with(f);
solve(u);

finalize_finch()

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
function EvaluationRmsChecker(arg0::jdouble)
    return EvaluationRmsChecker((jdouble,), arg0)
end

function EvaluationRmsChecker(arg0::jdouble, arg1::jdouble)
    return EvaluationRmsChecker((jdouble, jdouble), arg0, arg1)
end

function converged(obj::EvaluationRmsChecker, arg0::jint, arg1::LeastSquaresProblem_Evaluation, arg2::LeastSquaresProblem_Evaluation)
    return jcall(obj, "converged", jboolean, (jint, LeastSquaresProblem_Evaluation, LeastSquaresProblem_Evaluation), arg0, arg1, arg2)
end


text:
parse_repo(res::AbstractString) =
    parse(Core.ResRepoTree, joinpath(get_working_dir(), extract_basename(res)));
parse_repo(res::AbstractString, tt) =
    Core.ResRepoTree(get_working_dir(),res,tt);

function extract_basename(res::AbstractString)
    res_name = basename(res);
    res != res_name && warn("Only the file name $(res_name) is taken, the rest is ignored!!!");
    return res_name;
end


text:
function dDrs(d1, d2)
    n, m = map(min, size(d1), size(d2))
    weights = (m:-1:1) ./ m
    @views vd1 = d1[1:n, 1:m]
    @views vd2 = d2[1:n, 1:m]
    abs2.(vd1 - vd2) * weights ./ sum(weights) .|> sqrt
end

function analyze(rs, ϕs; fs = [normK, normT, normH], isD=true)
    label = ["ϕ$i" for i in 1:length(ϕs)]
    for (ϕ, lbl) in zip(ϕs, label)
        scatter(Tuple.(ϕ); title=lbl, label=nothing, markersize=2) |> display
    end

    for f in fs
        vs = [f(rs, ϕ) for ϕ in ϕs]
        plot(rs, vs; label=permutedims(label), title="$f(r)") |> display
    end

    if isD
        ds = [normD(length(ϕ) - 1, rs, ϕ) for ϕ in ϕs]
        [plot(ds[i], label=nothing, title=label[i]) for i in 1:length(ϕs)] .|> display

        d1 = ds[1]
        d2 = ds[end]

        dd = dDrs(d1, d2)
        plot(dd, title="∑Dk(r)/k") |> display
    end
    return
end

lossK(ϕ, p) = mapreduce((k, k0) -> abs2(k - k0), +, normK(p.rs, ϕ), p.k0) / length(p.rs)
lossH(ϕ, p) = mapreduce((h, h0) -> abs2(h - h0), +, normH(p.rs, ϕ), p.h0) / length(p.rs)

function errD(d1, d2)
    ddrs = PointsReconstruction.dDrs(d1, d2)
    sum(ddrs) / length(ddrs)
end

lossD(ϕ, p) = errD(normD(p.maxk, p.rs, ϕ), p.d0)

function default_loss_params(rs, ϕ; isK=false, isH=false, isD=false, isT=false, maxk=length(ϕ) - 1)
    p = (;rs)

    if isK
        k = (; k0 = normK(rs, ϕ))
        p = merge(p, k)
    end

    if isH
        h = (; h0 = normH(rs, ϕ))
        p = merge(p, h)
    end

    if isD
        d = (; maxk, d0 = normD(maxk, rs, ϕ))
        p = merge(p, d)
    end

    p
end


text:
"""
`AbstractRidgePredictor` is supposed to implement the interface
* `update_λs!`
* `trace_XtX`
* `XtXpΛ_ldiv_XtX`
* `LinearAlgebra.ldiv!`
* `Base.\`
Concrete subtypes available are `CholeskyRidgePredictor` and
`WoodburyRidgePredictor`.
"""
abstract type AbstractRidgePredictor end

"""
Used typically for p < n.
"""
struct CholeskyRidgePredictor{M<:AbstractMatrix,SYM<:Symmetric,C<:Cholesky} <:
       AbstractRidgePredictor
    X::M
    XtX::SYM
    XtXpΛ::SYM
    XtXpΛ_chol::C
end

function CholeskyRidgePredictor(X)
    (n, p) = size(X)
    XtX = Symmetric(X' * X ./ n)
    XtXpΛ = XtX + 1.0 * I
    XtXpΛ_chol = cholesky!(XtXpΛ)
    CholeskyRidgePredictor(X, XtX, XtXpΛ, XtXpΛ_chol)
end

function update_λs!(chol::CholeskyRidgePredictor, groups, λs)
    chol.XtXpΛ .= Symmetric(chol.XtX + Diagonal(group_expand(groups, λs)))
    cholesky!(chol.XtXpΛ)
end

function ldiv!(A, chol::CholeskyRidgePredictor, B)
    ldiv!(A, chol.XtXpΛ_chol, B)
end

function \(chol::CholeskyRidgePredictor, B)
    chol.XtXpΛ_chol \ B
end

function XtXpΛ_ldiv_XtX(chol::CholeskyRidgePredictor)
    chol.XtXpΛ_chol \ chol.XtX
end

function trace_XtX(chol::CholeskyRidgePredictor)
    tr(chol.XtX)
end

"""
Used typically for p >> n and n reasonably small
"""
struct WoodburyRidgePredictor{M<:AbstractMatrix,S<:SymWoodbury} <:
       SigmaRidgeRegression.AbstractRidgePredictor
    X::M
    wdb::S
end

function WoodburyRidgePredictor(X)
    (n, p) = size(X)
    wdb = SymWoodbury(1.0 * I(p), X', I(n) / n)
    WoodburyRidgePredictor(X, wdb)
end

#It would be very useful if there could be an implementation of
#```julia
#ldiv!(dest::AbstracMatrix, W::AbstractWoodbury, B::AbstractMatrix)
#```
#Right now I think this only works with `AbstractVector`. Before implementing and filing a pull request, I was wondering whether you think it is an OK approach to

#```julia
# for i=1:ncols
#      ldiv!(view(dest, :, i), A, view(B,:,i))
#end
#```
#------------------------------------------------------------------------
# TODO: Fix the following two things upstream on WoodburyMatrices.jl
#------------------------------------------------------------------------
function _ldiv!(dest, W::SymWoodbury, A::Diagonal, B)
    WoodburyMatrices.myldiv!(W.tmpN1, A, B)
    mul!(W.tmpk1, W.V, W.tmpN1)
    mul!(W.tmpk2, W.Cp, W.tmpk1)
    mul!(W.tmpN2, W.U, W.tmpk2)
    WoodburyMatrices.myldiv!(A, W.tmpN2)
    for i = 1:length(W.tmpN2)
        @inbounds dest[i] = W.tmpN1[i] - W.tmpN2[i]
    end
    return dest
end

#-----------------------------------------------------------------------
function LinearAlgebra.ldiv!(Y::AbstractMatrix, A::SymWoodbury, B::AbstractMatrix)
    ncols = size(B, 2)
    for i = 1:ncols
        ldiv!(view(Y, :, i), A, view(B, :, i))
    end
    Y
end
#-----------------------------------------------------------------------

function update_λs!(wbpred::WoodburyRidgePredictor, groups, λs)
    wdb = wbpred.wdb
    n = size(wdb.D, 1)
    A = Diagonal(group_expand(groups, λs))
    wdb.A .= A
    wdb.Dp .= inv(n * I + wdb.B' * (A \ wdb.B))
end

function ldiv!(A, wbpred::WoodburyRidgePredictor, B)
    ldiv!(A, wbpred.wdb, B)
end

function \(wbpred::WoodburyRidgePredictor, B)
    wbpred.wdb \ B
end

function XtXpΛ_ldiv_XtX(wbpred::WoodburyRidgePredictor)
    n = size(wbpred.X, 1)
    (wbpred.wdb \ wbpred.X') * wbpred.X ./ n
end

function trace_XtX(wbpred::WoodburyRidgePredictor)
    n = size(wbpred.X, 1)
    # recall XtX here really is XtX/n
    tr(wbpred.X' * wbpred.X) / n #make more efficient later.
end

Base.@kwdef mutable struct BasicGroupRidgeWorkspace{
    CP<:AbstractRidgePredictor,
    M<:AbstractMatrix,
    V<:AbstractVector,
}
    X::M
    Y::V
    groups::GroupedFeatures
    n::Integer = size(X, 1)
    p::Integer = size(X, 2)
    λs::V = ones(groups.num_groups)
    XtY::V = X' * Y ./ n
    XtXpΛ_chol::CP = CholeskyRidgePredictor(X)
    XtXpΛ_div_Xt::M = XtXpΛ_chol \ X' .\ n
    β_curr::V = XtXpΛ_chol \ XtY
    leverage_store::V = zeros(n)
    Y_hat::V = X * β_curr
    cache = nothing
end

ngroups(rdg::BasicGroupRidgeWorkspace) = ngroups(rdg.groups)

# StatsBase.jl interace
coef(rdg::BasicGroupRidgeWorkspace) = rdg.β_curr
islinear(rdg::BasicGroupRidgeWorkspace) = true
leverage(rdg::BasicGroupRidgeWorkspace) = rdg.leverage_store
modelmatrix(rdg::BasicGroupRidgeWorkspace) = rdg.X
predict(rdg::BasicGroupRidgeWorkspace, X) = X * coef(rdg)
response(rdg::BasicGroupRidgeWorkspace) = rdg.Y

function loo_error(rdg::BasicGroupRidgeWorkspace)
    mean(abs2.((rdg.Y .- rdg.Y_hat) ./ (1.0 .- rdg.leverage_store)))
end

function mse_ridge(rdg::BasicGroupRidgeWorkspace, X_test, Y_test)
    mean(abs2.(Y_test - X_test * rdg.β_curr))
end

function StatsBase.fit!(rdg::BasicGroupRidgeWorkspace, λs)
    λs = isa(λs, MutableNamedTuple) ? collect(values(λs)) : λs
    rdg.λs .= λs
    update_λs!(rdg.XtXpΛ_chol, rdg.groups, λs)
    #rdg.XtXpΛ .= Symmetric(rdg.XtX + Diagonal(group_expand(rdg.groups, λs)))
    #cholesky!(rdg.XtXpΛ)
    ldiv!(rdg.β_curr, rdg.XtXpΛ_chol, rdg.XtY)
    mul!(rdg.Y_hat, rdg.X, rdg.β_curr)
    ldiv!(rdg.XtXpΛ_div_Xt, rdg.XtXpΛ_chol, rdg.X')
    rdg.XtXpΛ_div_Xt ./= rdg.n
    _prod_diagonals!(rdg.leverage_store, rdg.X, rdg.XtXpΛ_div_Xt)
    loo_error(rdg)
end

"""
    λωλας_λ(rdg; multiplier=0.1)

Implements the Panagiotis Lolas rule of thumb for picking an optimal λ.
"""
function λωλας_λ(rdg; multiplier = 0.1)
    multiplier * rdg.p^2 / rdg.n / trace_XtX(rdg.XtXpΛ_chol) #TODO 2s
end

#function max_σ_squared(rdg)
#   mean(abs2, rdg.Y)
#end

# Tuning through Moment Fitting

Base.@kwdef struct MomentTunerSetup{
    IV<:AbstractVector,
    FV<:AbstractVector,
    FM<:AbstractMatrix,
}
    ps::IV
    n::Integer
    beta_norms_squared::FV
    N_norms_squared::FV
    M_squared::FM
end

function MomentTunerSetup(rdg::BasicGroupRidgeWorkspace)
    grps = rdg.groups
    n = rdg.n
    ps = grps.ps
    ngroups = grps.num_groups
    beta_norms_squared = group_summary(grps, rdg.β_curr, x -> sum(abs2, x))
    N_matrix = rdg.XtXpΛ_div_Xt #sqrt(n)*N from paper
    M_matrix = XtXpΛ_ldiv_XtX(rdg.XtXpΛ_chol) #TODO 1
    N_norms_squared = Vector{eltype(beta_norms_squared)}(undef, ngroups)
    M_squared = Matrix{eltype(beta_norms_squared)}(undef, ngroups, ngroups)

    for g = 1:ngroups
        N_norms_squared[g] = sum(abs2, N_matrix[group_idx(grps, g), :])
        for h = 1:ngroups
            # Mij  is entry (j,i)  and dived by p_i
            M_squared[h, g] = sum(abs2, M_matrix[group_idx(grps, g), group_idx(grps, h)])
        end
    end
    MomentTunerSetup(
        ps = ps,
        n = n,
        beta_norms_squared = beta_norms_squared,
        N_norms_squared = N_norms_squared,
        M_squared = M_squared,
    )

end

function σ_squared_max(mom::MomentTunerSetup)
    u = mom.beta_norms_squared
    v = mom.N_norms_squared
    maximum(u ./ v)
end

function get_αs_squared(mom::MomentTunerSetup, σ_squared)
    rhs = mom.beta_norms_squared .- σ_squared .* mom.N_norms_squared
    α_sq_by_p = vec(nonneg_lsq(mom.M_squared, rhs; alg = :fnnls)) #  mom.M_squared\rhs\
    α_sq_by_p .* mom.ps
end

function get_λs(mom::MomentTunerSetup, σ_squared)
    αs_squared = get_αs_squared(mom, σ_squared)
    γs = mom.ps ./ mom.n
    σ_squared .* γs ./ αs_squared
end

function sigma_squared_path(
    rdg::BasicGroupRidgeWorkspace,
    mom::MomentTunerSetup,
    σs_squared,
)

    n_σs = length(σs_squared)
    n_groups = ngroups(rdg)
    loos_hat = zeros(n_σs)
    λs = zeros(n_σs, n_groups)
    βs = zeros(n_σs, rdg.groups.p)
    for (i, σ_squared) in enumerate(σs_squared)
        λs_tmp = get_λs(mom, σ_squared)
        λs[i, :] = λs_tmp
        loos_hat[i] = fit!(rdg, λs_tmp)
        βs[i, :] = rdg.β_curr
    end
    (λs = λs, loos = loos_hat, βs = βs)
end


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

using Test, Random
import Base64:
    Base64EncodePipe,
    base64encode,
    Base64DecodePipe,
    base64decode

const inputText = "Man is distinguished, not only by his reason, but by this singular passion from other animals, which is a lust of the mind, that by a perseverance of delight in the continued and indefatigable generation of knowledge, exceeds the short vehemence of any carnal pleasure."
const encodedMaxLine76 = """
TWFuIGlzIGRpc3Rpbmd1aXNoZWQsIG5vdCBvbmx5IGJ5IGhpcyByZWFzb24sIGJ1dCBieSB0aGlz
IHNpbmd1bGFyIHBhc3Npb24gZnJvbSBvdGhlciBhbmltYWxzLCB3aGljaCBpcyBhIGx1c3Qgb2Yg
dGhlIG1pbmQsIHRoYXQgYnkgYSBwZXJzZXZlcmFuY2Ugb2YgZGVsaWdodCBpbiB0aGUgY29udGlu
dWVkIGFuZCBpbmRlZmF0aWdhYmxlIGdlbmVyYXRpb24gb2Yga25vd2xlZGdlLCBleGNlZWRzIHRo
ZSBzaG9ydCB2ZWhlbWVuY2Ugb2YgYW55IGNhcm5hbCBwbGVhc3VyZS4="""

@testset "Examples" begin
    # Encode and decode
    fname = tempname()
    open(fname, "w") do f
        opipe = Base64EncodePipe(f)
        write(opipe,inputText)
        @test close(opipe) === nothing
    end

    open(fname, "r") do f
        ipipe = Base64DecodePipe(f)
        @test read(ipipe, String) == inputText
        @test close(ipipe) === nothing
    end
    rm(fname)

    # Byte-by-byte encode and decode.
    buf = IOBuffer()
    pipe = Base64EncodePipe(buf)
    for char in inputText
        write(pipe, UInt8(char))
    end
    close(pipe)
    pipe = Base64DecodePipe(IOBuffer(take!(buf)))
    decoded = UInt8[]
    while !eof(pipe)
        push!(decoded, read(pipe, UInt8))
    end
    @test String(decoded) == inputText

    # Encode to string and decode
    @test String(base64decode(base64encode(inputText))) == inputText

    # Decode with max line chars = 76 and padding
    ipipe = Base64DecodePipe(IOBuffer(encodedMaxLine76))
    @test read(ipipe, String) == inputText

    # Decode with max line chars = 76 and no padding
    #ipipe = Base64DecodePipe(IOBuffer(encodedMaxLine76[1:end-1]))
    #@test read(ipipe, String) == inputText

    # Decode with two padding characters ("==")
    ipipe = Base64DecodePipe(IOBuffer(string(encodedMaxLine76[1:end-2],"==")))
    @test read(ipipe, String) == inputText[1:end-1]

    # Test incorrect format
    ipipe = Base64DecodePipe(IOBuffer(encodedMaxLine76[1:end-3]))
    @test_throws ArgumentError read(ipipe, String)

    # issue #21314
    @test base64decode(chomp("test")) == base64decode("test")
end

@testset "Random data" begin
    mt = MersenneTwister(1234)
    for _ in 1:1000
        data = rand(mt, UInt8, rand(0:300))
        @test hash(base64decode(base64encode(data))) == hash(data)
    end
end

max_stars_repo_name: nreigl/julia

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI


text:
include("./predict.jl")
using Test

@test predict(34.1330949, -117.9143879, 33.8068768, -118.3527671, 3600) isa Number
@test predict(34.1330949, -117.9143879, 33.8068768, -118.3527671, 3600) < 1
@test predict(34.1330949, -117.9143879, 33.8068768, -118.3527671, 3600) >= 0

@test read_response() isa Dict

println("Test Successful")

max_issues_repo_name: MLH-Fellowship/CodeVidLive


text:
using CEnum
using Test

@cenum(ZeroOne, zero, one)
@test zero == 0
@test one == 1

@cenum(Fruit, apple=1, orange=2, kiwi=2)
@test orange == kiwi
@test apple | orange == 3
@test apple & orange == 0
@test apple + 1 == 2
@test kiwi - 1 == 1
@test kiwi ⊻ kiwi == 0
@test ~orange == -3
@test_nowarn print(devnull, Fruit(apple | orange))

@cenum(Boolean::Bool, alternativefact, fact)
@test alternativefact == false

@cenum(Day, Mon=1, Tue, Wed=3, suiyoubi=3, Fri=5, Sat)
@test Mon == 1
@test Tue == 2
@test Wed == 3 == suiyoubi
@test Fri == 5
@test Sat == 6

# issue #11
@test 0 - one == -1


text:
#Reineke SDI calculation
function sdi(tpa,qmd)
    sdi=tpa*(qmd/10.0)^1.605
end

function qmd(ba, tpa, constant)
    qmd=sqrt((ba/tpa)/constant)
end