text:
module emodel

# Uses an eccentric disk

export write_grid, write_model, write_lambda, write_dust, Parameters, Grid

using ..constants

# Write the wavelength sampling file. Only run on setup
function write_lambda(lams::Array{Float64, 1}, basedir::AbstractString)
    fcam = open(basedir * "camera_wavelength_micron.inp", "w")
    nlam = length(lams)
    @printf(fcam, "%d\n", nlam)
    for lam in lams
        @printf(fcam, "%.9e\n", lam) # [microns]
    end

    close(fcam)
end

const eqmirror = true # mirror the grid about the z=0 midplane ?
# if we decide to mirror, then ncells = 1/2 of the true value

# Define a grid object which stores all of these variables
# This will not change for the duration of the run
immutable Grid
    nr::Int
    ntheta::Int
    nphi::Int
    ncells::Int
    # cell edges
    Rs::Vector{Float64}
    Thetas::Vector{Float64}
    Phis::Vector{Float64}
    As::Vector{Float64}
    # cell centers
    rs::Vector{Float64}
    thetas::Vector{Float64}
    phis::Vector{Float64}
    as::Vector{Float64}
    a_widths::Vector{Float64}
end

function Grid(nr::Int, ntheta::Int, nphi::Int, r_in::Real, r_out::Real, na::Int, eqmirror::Bool)
    # Specify a 2D axisymmetric *separable* grid in spherical coordinates:
    # {r, theta, phi}, where theta is angle from zenith, phi is azimuth

    # na is the number of eccentric bins that we want to have. Generally, I
    # think this should be less than nr.

    ncells = nr * ntheta * nphi
    r_in = convert(Float64, r_in) * AU # Inner extent of disk
    r_out = convert(Float64, r_out) * AU # Outer extent of disk

    #Define the cell *walls*
    Rs = logspace(log10(r_in), log10(r_out), nr+1) # [cm] logarithmically spaced

    if eqmirror
        ped = 0.1
        #Thetas = linspace(0, pi/2., ntheta+1)
        # [rad] Angles are internally defined in radians, not degrees
        Thetas = pi/2. - (logspace(log10(ped), log10(pi/2. + ped), ntheta+1) - ped)[end:-1:1]
        #Spaced closer near the z=0
    else
        Thetas = linspace(0, pi, ntheta+1)
        # [rad] Angles are internally defined in radians, not degrees
    end

    Phis = linspace(0, 2pi, nphi + 1) # [rad] cell walls for inactive coordinate

    # Define the semi-major axis bin edges such that we fully encapsulate the full range of radii
    # Assume 0 <= e <= 0.7
    eps = 0.1 * AU # [cm], tiny extra bit of comfort
    e = 0.7 # Maximum e we'll tolerate
    A_1 = (Rs[1] - eps)/(1 + e)
    A_N = (Rs[end] + eps)/(1 - e)

    As = logspace(log10(A_1), log10(A_N), na+1)

    #Define the cell centers as the average between walls
    rs = 0.5 * (Rs[1:end-1] + Rs[2:end])
    thetas = 0.5 * (Thetas[1:end-1] + Thetas[2:end])
    phis = 0.5 * (Phis[1:end-1] + Phis[2:end])

    as = 0.5 * (As[1:end-1] + As[2:end])
    a_widths = As[2:end] - As[1:end-1]

    return Grid(nr, ntheta, nphi, ncells, Rs, Thetas, Phis, As, rs, thetas, phis, as, a_widths)

end

#This function only needs to be run once, upon setup.
function write_grid(basedir::AbstractString, grid::Grid)
    #amr_grid.inp
    f = open(basedir * "amr_grid.inp", "w")

    #Write the header
    @printf(f, "%d\n", 1) #iformat
    @printf(f, "%d\n", 0) #regular grid (no AMR or Oct-tree)
    @printf(f, "%d\n", 100) #spherical coordiantes
    @printf(f, "%d\n", 0) #gridinfo (none needed for now)
    #incl_r incl_phi incl_z #use this axis?
    @printf(f, "%d %d %d \n", 1, 1, 1) # Full 3D grid
    #n_r    n_phi   n_z #number of cells in this dimension
    @printf(f, "%d %d %d \n", grid.nr, grid.ntheta, grid.nphi)

    # Write out the cell walls
    for R in grid.Rs
        @printf(f, "%.9e\n", R)
    end

    for Theta in grid.Thetas
        @printf(f, "%.9e\n", Theta)
    end

    for Phi in grid.Phis
        @printf(f, "%.9e\n", Phi)
    end

    close(f)
end

#Let's try defining a parameters type, an object of which gets passed around.

type Parameters
    M_star::Float64 # [M_sun] stellar mass
    a_c::Float64 # [AU] characteristic radius
    T_10::Float64 # [K] temperature at 10 AU
    q::Float64 # temperature gradient exponent
    gamma::Float64 # surface temperature gradient exponent
    M_CO::Float64 # [M_earth] disk mass of CO
    ksi::Float64 # [cm s^{-1}] microturbulence
    dpc::Float64 # [pc] distance to system
    incl::Float64 # [degrees] inclination 0 deg = face on, 90 = edge on.
    PA::Float64 # [degrees] position angle (East of North)
    e::Float64 # eccentricity
    w::Float64 # [degrees] argument of periastron
    vel::Float64 # [km/s] systemic velocity (positive is redshift/receeding)
    mu_RA::Float64 # [arcsec] central offset in RA
    mu_DEC::Float64 # [arcsec] central offset in DEC
end

# Given an array of bin edges (N) in radius from the star, and an array of radii (xx), return an
# array of indices of which bins (N-1) the radii fall into.
# Assumes both bins and xx are sorted in increasing order.
function get_bins(bin_edges, xx)
    inds = Array(Int64, length(xx))
    b::Int64 = 1
    nb = length(bin_edges)
    for (i,x) in enumerate(xx)
        while true
            # Greater than the lower bin
            if x >= bin_edges[b]

                # Smaller than the upper bin
                if x < bin_edges[b+1]
                    inds[i] = b
                    break

                # larger than the upper bin, means we need to try a larger bin index
                else
                    b += 1
                    if b + 1 > nb
                        throw(error(@sprintf("x value % s out of bin range.", x)))
                    end
                end

            # Smaller than the lower bin, something went wrong?
            else
                throw(error(@sprintf("x value %s out of bin range.", x)))

            end
        end

    end
    return inds
end

# Assume all inputs to these functions are in CGS units and in *cylindrical* coordinates.

# Assume all angles are in *radians*, including argument of periapse

# phi is azimuthal angle, a is semi-major axis, e is eccentricity
function radius{T}(phi::T, a::Float64, e::Float64, w::Float64)
    a * (1 - e^2)./(1 + e .* cos(phi - w))
end

radius{T}(phi::T, a::Float64, pars::Parameters) = radius(phi, a, pars.e, pars.w * deg)

# Linear density of the ring for the mass-on-a-wire approximation
# m is the total mass of the ring
function lambda(phi::Float64, a::Float64, e::Float64, w::Float64, m::Float64)
    m * sqrt(1 - e^2) / (2pi * a * sqrt(1 + 2e * cos(phi - w) + e^2))
end

# Parametric type T allows passing individual Float64 or Vectors.
# Alternate functions accept pars passed around, where pars is in M_star, AU, etc...
function velocity(phi::Float64, a::Float64, e::Float64, w::Float64, M_star::Float64)
    v_r = sqrt(G * M_star / a) * e * sin(phi - w) / sqrt(1 - e^2)
    v_theta = 0.0
    v_phi = sqrt(G * M_star / a) * (1 + e * cos(phi - w)) / sqrt(1 - e^2)
    return Float64[v_r, v_theta, v_phi]
end
velocity(phi::Float64, a::Float64, pars::Parameters) = velocity(phi, a, pars.e, pars.w * deg, pars.M_star * M_sun)

function temperature{T}(r::T, T_10::Float64, q::Float64)
    T_10 * (r ./ (10. * AU)).^(-q)
end
temperature{T}(r::T, pars::Parameters) = temperature(r, pars.T_10, pars.q)

function Hp{T}(r::T, M_star::Float64, T_10::Float64, q::Float64)
    temp = temperature(r, T_10, q)
    sqrt(kB * temp .* r.^3./(mu_gas * m_H * G * M_star))
end
Hp{T}(r::T,  pars::Parameters) = Hp(r, pars.M_star * M_sun, pars.T_10, pars.q)

# No parametric type for number density, because it is a 2D function.
function n_CO(r::Float64, z::Float64, M_star::Float64, T_10::Float64, q::Float64)
    H = Hp(r, M_star, T_10, q)
    return 1./ (m_CO * sqrt(2pi) * H) * exp(-0.5 * (z/ H)^2)
end

n_CO(r::Float64, z::Float64, pars::Parameters) = n_CO(r, z, pars.M_star * M_sun, pars.T_10, pars.q)

function n_nobin(r::Float64, phi::Float64, z::Float64, M_star::Float64, a_c::Float64, T_10::Float64, q::Float64, gamma::Float64, M_CO::Float64, e::Float64, w::Float64)
    n = n_CO(r, z, M_star, T_10, q)
    a = r * (1 + e * cos(phi - w))/(1 - e^2)
    Sigma = (2 - gamma) * M_CO / (2pi * a * a_c) * (a/a_c)^(1 - gamma) * exp(- (a/a_c)^(2 - gamma)) * (1 + e * cos(phi - w))/(sqrt(1 + 2 * e * cos(phi - w)) * sqrt(1 - e^2))
    return Sigma * n
end

n_nobin(r::Float64, phi::Float64, z::Float64, pars::Parameters) = n_nobin(r, phi, z, pars.M_star * M_sun, pars.a_c * AU, pars.T_10, pars.q, pars.gamma, pars.M_CO * M_earth, pars.e, pars.w * deg)

# Calculate the mass in each ring using a_c and M_co
function ring_mass(as::Vector{Float64}, a_c::Float64, M_CO::Float64, gamma::Float64)
    a = as ./ a_c
    ms = M_CO .* a.^(-gamma) .* exp(-a.^(2 - gamma)) / sum(a.^(-gamma) .* exp(-a.^(2 - gamma)))
    return ms
end

ring_mass(grid, pars) = ring_mass(grid.as, pars.a_c * AU, pars.M_CO * M_earth, pars.gamma)

function rho_dust(r::Float64, z::Float64, pars::Parameters)
    nCO = n_CO(r, z, pars) # number of CO molecules per cm^3

    # Convert from nCO to nH2
    nH2 = nCO / 7.e-5 # number density ratio

    # Convert from nH2 (assuming nH2 ~ nGas ) to mGas
    mGas = constants.m0 * nH2 # [g]

    # Convert from mGas to mDust using Gas/Dust ratio of 100
    mDust = mGas * 0.01 # [g]

    return mDust
end

# Ksi is microturbulent broadining width in units of km/s. Output of this function
# is in cm/s for RADMC (RADMC manual, eqn 7.12)
function microturbulence(ksi::Float64)
    return ksi * 1.e5 # convert from km/s to cm/s
end

microturbulence(pars::Parameters) = microturbulence(pars.ksi)

function write_model(pars::Parameters, basedir::AbstractString, grid::Grid)
    # numberdens_co.inp
    fdens = open(basedir * "numberdens_co.inp", "w")
    @printf(fdens, "%d\n", 1) #iformat
    @printf(fdens, "%d\n", grid.ncells)

    # gas_velocity.inp
    fvel = open(basedir * "gas_velocity.inp", "w")
    @printf(fvel, "%d\n", 1) #iformat
    @printf(fvel, "%d\n", grid.ncells)

    # gas_temperature.inp
    ftemp = open(basedir * "gas_temperature.inp", "w")
    @printf(ftemp, "%d\n", 1) #iformat
    @printf(ftemp, "%d\n", grid.ncells)

    # microturbulence.inp
    fmicro = open(basedir * "microturbulence.inp", "w")
    @printf(fmicro, "%d\n", 1) #iformat
    @printf(fmicro, "%d\n", grid.ncells)

    # Now, we will need to write the three other files as a function of grid position.
    # Therefore we will do *one* loop over these indices, calculate the required value,
    # and write it to the appropriate file.

    # Calculate the mass in each ring using a_c and M_co

    #Looping over the cell centers
    for phi in grid.phis
        # Calculate the ellipse bin-edges
        # R_edges = grid.As .* (1 - pars.e^2)/(1 + pars.e * cos(phi - pars.w * deg))
        #
        # # println("grid.As", grid.As)
        # # println("R_edges", R_edges)
        # # println("grid.rs", grid.rs)
        #
        # # Determine the array of bin indices which correspond to each point in grid.rs
        # inds = get_bins(R_edges, grid.rs)
        #
        # # Calculate all following phi-dependent quantities in vectorized form:
        # # m_i, a_i, dr_i, Sigma_CO_i
        #
        # m_rings = ring_mass(grid, pars)
        #
        # # Calculate arrays of these values for each radius
        # ms = [m_rings[ind] for ind in inds]
        # as = [grid.as[ind] for ind in inds]
        # das = [grid.a_widths[ind] for ind in inds]
        #
        # drs = das .* (1 - pars.e^2)/(1 + pars.e * cos(phi - pars.w * deg))
        #
        # lambdas = ms ./ (2pi .* as) .* (sqrt(1 - pars.e^2) ./ sqrt(1 + 2 * pars.e .* cos(phi - pars.w * deg) + pars.e^2))
        #
        # Sigmas = lambdas ./ drs

        for theta in grid.thetas

            for (i,r) in enumerate(grid.rs)
                #Convert from spherical to cylindrical coordinates
                z = r * cos(theta)
                r_cyl = r * sin(theta)

                a = r_cyl * (1 + pars.e * cos(phi - pars.w * deg))/(1 - pars.e^2)
                #Calculate unscaled density as function of (r_cyl, z)
                # n = Sigmas[i] * n_CO(r_cyl, z, pars)

                v = velocity(phi, a, pars)

                # @printf(fdens, "%.9e\n", n)
                @printf(fdens, "%.9e\n", n_nobin(r_cyl, phi, z, pars))
                @printf(fvel, "%.9e %.9e %.9e\n", v[1], v[2], v[3])
                @printf(ftemp, "%.9e\n", temperature(r_cyl, pars))
                @printf(fmicro, "%.9e\n", microturbulence(pars))
            end
        end
    end

    close(fdens)
    close(fvel)
    close(ftemp)
    close(fmicro)

end

function write_dust(pars::Parameters, basedir::AbstractString, grid::Grid)
    fdens = open(basedir * "dust_density.inp", "w")
    @printf(fdens, "%d\n", 1) #iformat
    @printf(fdens, "%d\n", grid.ncells)
    @printf(fdens, "%d\n", 1) # number of dust species

    for phi in grid.phis
        for theta in grid.thetas
            for r in grid.rs
                #Convert from spherical to cylindrical coordinates
                z = r * cos(theta)
                r_cyl = r * sin(theta)

                @printf(fdens, "%.9e\n", rho_dust(r_cyl, z, pars))
            end
        end
    end

    close(fdens)
end

end

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
using ModelConstructors, HDF5, Random, JLD2, FileIO, SMC, Test
include("modelsetup.jl")

path = dirname(@__FILE__)

if VERSION < v"1.5"
    ver = "111"
else
    ver = "150"
end

m    = setup_linear_model(; regime_switching = false)
m_rs = setup_linear_model(; regime_switching = true)
save = normpath(joinpath(dirname(@__FILE__),"save"))
m    <= Setting(:saveroot, save)
m_rs <= Setting(:saveroot, save)
savepath = rawpath(m, "estimate", "smc_cloud.jld2")
particle_store_path = rawpath(m, "estimate", "smcsave.h5")

data    = h5read("reference/test_data.h5", "data")
data_rs = h5read("reference/test_data.h5", "rsdata")

Random.seed!(1793)

# mean_para = mean(SMC.get_vals(test_cloud), dims = 2)
true_para = [1., 1., 1., # α1, β1, σ1 (regime 1)
             2., 2., 1., # α2, β2, σ2 (regime 1)
             3., 3., 1., # α3, β3, σ3 (regime 1)
             1., 1.,     # α1 regimes = 2-3
             2., 3.,     # β1 regimes = 2-3
             2., 2.,     # α2 regimes = 2-3
             3., 4.,     # β2 regimes = 2-3
             3., 3.,     # α3 regimes = 2-3
             4., 5.]     # β3 regimes = 2-3

saved_file  = JLD2.jldopen(string("reference/smc_cloud_fix=true_version=", ver, ".jld2"), "r")
saved_cloud = saved_file["cloud"]
saved_w     = saved_file["w"]
saved_W     = saved_file["W"]

saved_file_rs  = JLD2.jldopen(string("reference/smc_cloud_fix=true_rs=true_version=", ver, ".jld2"), "r")
saved_cloud_rs = saved_file_rs["cloud"]
saved_w_rs     = saved_file_rs["w"]
saved_W_rs     = saved_file_rs["W"]

n_para    = length(m.parameters)
n_para_rs = n_parameters_regime_switching(m_rs)
old_para_inds = vcat(trues(n_para), falses(n_para_rs - n_para))

m_rs.parameters[7].fixed = false # need to unfix this or the log-likelihoods won't match
m_rs.parameters[7].valuebounds = m_rs.parameters[4].valuebounds
set_regime_fixed!(m_rs.parameters[7], 1, false) # update! doesn't seem to be working properly???
set_regime_valuebounds!(m_rs.parameters[7], 1, m_rs.parameters[4].valuebounds)
new_cloud = SMC.add_parameters_to_cloud(saved_cloud, m_rs.parameters,
                                        old_para_inds; regime_switching = true)
new_cloud2 = SMC.add_parameters_to_cloud(joinpath("reference", "smc_cloud_fix=true_version=" *
                                                  ver * ".jld2"), m_rs.parameters,
                                         old_para_inds; regime_switching = true)

@test saved_cloud.particles[:, 1:9] ≈ new_cloud.particles[:, 1:9]
@test saved_cloud.particles[:, end - 2:end] ≈ new_cloud.particles[:, end - 2:end] # old_loglh, accept, and weight should be the same
@test saved_cloud.particles[:, end - 4] ≈ new_cloud.particles[:, end - 4] # loglh should be the same
@test !(saved_cloud.particles[:, end - 3] ≈ new_cloud.particles[:, end - 3]) # prior should not be the same b/c new parameters & resampling from prior
@test saved_cloud.particles[:, 1:9] ≈ new_cloud2.particles[:, 1:9]
@test saved_cloud.particles[:, end - 2:end] ≈ new_cloud2.particles[:, end - 2:end] # old_loglh, accept, and weight should be the same
@test saved_cloud.particles[:, end - 4] ≈ new_cloud2.particles[:, end - 4] # loglh should be the same
@test !(saved_cloud.particles[:, end - 3] ≈ new_cloud2.particles[:, end - 3]) # prior should not be the same b/c new parameters & resampling from prior

ext: jl


text:
using Test
using OverflowContexts

for i ∈ (1, 2)
    if i == 1
        testname = "checked default"
        @default_checked

        @testset "$testname" begin
            @test_throws OverflowError -typemin(Int)
            @test_throws OverflowError -UInt(1)

            @test_throws OverflowError typemax(Int) + 1
            @test_throws OverflowError typemax(UInt) + 1

            @test_throws OverflowError typemin(Int) - 1
            @test_throws OverflowError typemin(UInt) - 1

            @test_throws OverflowError typemax(Int) * 2
            @test_throws OverflowError typemin(Int) * 2
            @test_throws OverflowError typemax(UInt) * 2

            @test_throws OverflowError typemax(Int) ^ 2
            @test_throws OverflowError typemin(Int) ^ 2
            @test_throws OverflowError typemax(UInt) ^ 2

            @test_throws OverflowError abs(typemin(Int))
        end
    else
        testname = "unchecked default"
        @default_unchecked

        @testset "$testname" begin
            @test -typemin(Int) == typemin(Int)
            @test -UInt(1) == typemax(UInt)

            @test typemax(Int) + 1 == typemin(Int)
            @test typemax(UInt) + 1 == typemin(UInt)

            @test typemin(Int) - 1 == typemax(Int)
            @test typemin(UInt) - 1 == typemax(UInt)

            @test typemax(Int) * 2 == -2
            @test typemin(Int) * 2 == 0
            @test typemax(UInt) * 2 == typemax(UInt) - 1

            @test typemax(Int) ^ 2 == 1
            @test typemin(Int) ^ 2 == 0
            @test typemax(UInt) ^ 2 == UInt(1)

            @test abs(typemin(Int)) == typemin(Int)
        end
    end

    @testset "$testname with expressions" begin
        @testset "checked expressions" begin
            @test_throws OverflowError @checked -typemin(Int)
            @test_throws OverflowError @checked -UInt(1)

            @test_throws OverflowError @checked typemax(Int) + 1
            @test_throws OverflowError @checked typemax(UInt) + 1

            @test_throws OverflowError @checked typemin(Int) - 1
            @test_throws OverflowError @checked typemin(UInt) - 1

            @test_throws OverflowError @checked typemax(Int) * 2
            @test_throws OverflowError @checked typemin(Int) * 2
            @test_throws OverflowError @checked typemax(UInt) * 2

            @test_throws OverflowError @checked typemax(Int) ^ 2
            @test_throws OverflowError @checked typemin(Int) ^ 2
            @test_throws OverflowError @checked typemax(UInt) ^ 2

            @test_throws OverflowError @checked abs(typemin(Int))
        end

        @testset "unchecked expressions" begin
            @test @unchecked -typemin(Int) == typemin(Int)
            @test @unchecked -UInt(1) == typemax(UInt)

            @test @unchecked typemax(Int) + 1 == typemin(Int)
            @test @unchecked typemax(UInt) + 1 == typemin(UInt)

            @test @unchecked typemin(Int) - 1 == typemax(Int)
            @test @unchecked typemin(UInt) - 1 == typemax(UInt)

            @test @unchecked typemax(Int) * 2 == -2
            @test @unchecked typemin(Int) * 2 == 0
            @test @unchecked typemax(UInt) * 2 == typemax(UInt) - 1

            @test @unchecked typemax(Int) ^ 2 == 1
            @test @unchecked typemin(Int) ^ 2 == 0
            @test @unchecked typemax(UInt) ^ 2 == UInt(1)

            @test @unchecked abs(typemin(Int)) == typemin(Int)
        end

        @testset "lowest-level macro takes priority" begin
            @checked begin
                @test @unchecked typemax(Int) + 1 == typemin(Int)
            end
            @unchecked begin
                @test_throws OverflowError @checked typemax(Int) + 1
            end
        end

        @testset "non-integer math still works" begin
            @test @checked -1.0 == -1
            @test @unchecked -1.0 == -1
            @test @checked 1.0 + 3.0 == 4.0
            @test @unchecked 1.0 + 3.0 == 4.0
            @test @checked 1 + 3.0 == 4.0
            @test @unchecked 1 + 3.0 == 4.0
            @test @checked 1.0 - 3.0 == -2.0
            @test @unchecked 1.0 - 3.0 == -2.0
            @test @checked 1 - 3.0 == -2.0
            @test @unchecked 1 - 3.0 == -2.0
            @test @checked 1.0 * 3.0 == 3.0
            @test @unchecked 1.0 * 3.0 == 3.0
            @test @checked 1 * 3.0 == 3.0
            @test @unchecked 1 * 3.0 == 3.0
            @test @checked 1.0 ^ 3.0 == 1.0
            @test @unchecked 1.0 ^ 3.0 == 1.0
            @test @checked 1 ^ 3.0 == 1.0
            @test @unchecked 1 ^ 3.0 == 1.0
            @test @checked abs(-1.0) == 1.0
            @test @unchecked abs(-1.0) == 1.0
        end
    end
end

@default_unchecked # reset to Julia default

@testset "symbol replacement" begin
    expr = @macroexpand @checked foldl(+, [])
    @test expr.args[2] == :checked_add

    expr = @macroexpand @unchecked foldl(+, [])
    @test expr.args[2] == :unchecked_add

    expr = @macroexpand @checked foldl(-, [])
    @test expr.args[2] == :checked_negsub

    expr = @macroexpand @unchecked foldl(-, [])
    @test expr.args[2] == :unchecked_negsub

    expr = @macroexpand @checked foldl(*, [])
    @test expr.args[2] == :checked_mul

    expr = @macroexpand @unchecked foldl(*, [])
    @test expr.args[2] == :unchecked_mul

    expr = @macroexpand @checked foldl(^, [])
    @test expr.args[2] == :checked_pow

    expr = @macroexpand @unchecked foldl(^, [])
    @test expr.args[2] == :unchecked_pow

    expr = @macroexpand @checked foldl(:abs, [])
    @test expr.args[2] == :checked_abs

    expr = @macroexpand @unchecked foldl(:abs, [])
    @test expr.args[2] == :unchecked_abs
end

@testset "negsub helper methods dispatch correctly" begin
    @test unchecked_negsub(1) == -1
    @test unchecked_negsub(1, 2) == 1 - 2
end

@testset "in-place assignement" begin
    a = typemax(Int)
    @test_throws OverflowError @checked a += 1
    @unchecked a += 1
    @test a == typemin(Int)

    a = typemin(Int)
    @test_throws OverflowError @checked a -= 1
    @unchecked a -= 1
    @test a == typemax(Int)

    a = typemax(Int)
    @test_throws OverflowError @checked a *= 2
    @unchecked a *= 2
    @test a == -2

    a = typemax(Int)
    @test_throws OverflowError @checked a ^= 2
    @unchecked a ^= 2
    @test a == 1
end

module CheckedModule
    using OverflowContexts, Test
    @default_checked
    testfunc() = @test_throws OverflowError typemax(Int) + 1

    module NestedUncheckedModule
        using OverflowContexts, Test
        @default_unchecked
        testfunc() = @test typemax(Int) + 1 == typemin(Int)
    end
end

module UncheckedModule
    using OverflowContexts, Test
    @default_unchecked
    testfunc() = @test typemax(Int) + 1 == typemin(Int)
    module NestedCheckedModule
        using OverflowContexts, Test
        @default_checked
        testfunc() = @test_throws OverflowError typemax(Int) + 1
    end
end

@testset "module-specific contexts" begin
    CheckedModule.testfunc()
    CheckedModule.NestedUncheckedModule.testfunc()
    UncheckedModule.testfunc()
    UncheckedModule.NestedCheckedModule.testfunc()
end

@testset "ensure pow methods don't promote on the power" begin
    @test typeof(@unchecked 3 ^ UInt(4)) == Int
    @test typeof(@checked 3 ^ UInt(4)) == Int
end


text:
using DFTK: bzmesh_uniform, bzmesh_ir_wedge, ElementCoulomb, Vec3, Mat3
using DFTK: pymatgen_structure, load_lattice, standardize_atoms
using LinearAlgebra
using PyCall
using Test
include("testcases.jl")

@testset "bzmesh_uniform agrees with spglib" begin
    spglib = pyimport_conda("spglib", "spglib")

    function test_against_spglib(kgrid_size; kshift=[0, 0, 0])
        kgrid_size = Vec3(kgrid_size)
        identity = [reshape(Mat3{Int}(I), 1, 3, 3)]
        is_shift = ifelse.(kshift .== 0, 0, 1)
        _, grid = spglib.get_stabilized_reciprocal_mesh(kgrid_size, identity,
                                                        is_shift=is_shift)
        kcoords_spglib = [(kshift .+ Vec3{Int}(grid[ik, :])) .// kgrid_size
                          for ik in 1:size(grid, 1)]
        kcoords_spglib = DFTK.normalize_kpoint_coordinate.(kcoords_spglib)
        sort!(kcoords_spglib)

        kcoords, _ = bzmesh_uniform(kgrid_size, kshift=kshift)
        sort!(kcoords)

        @test kcoords == kcoords_spglib
    end

    test_against_spglib([ 2,  3,  2])
    test_against_spglib([ 3,  3,  3])
    test_against_spglib([ 3,  3,  3], kshift=[1//2, 0, 0])
    test_against_spglib([ 2,  3,  4])
    test_against_spglib([ 9, 11, 13])
end

@testset "bzmesh_ir_wedge is correct reduction" begin
    function test_reduction(system, kgrid_size; supercell=[1, 1, 1], kshift=[0, 0, 0])
        lattice = system.lattice
        atoms = [ElementCoulomb(system.atnum) => system.positions]
        if supercell != [1, 1, 1]  # Make a supercell
            pystruct = pymatgen_structure(lattice, atoms)
            pystruct.make_supercell(supercell)
            lattice = load_lattice(pystruct)
            el = ElementCoulomb(system.atnum)
            atoms = [el => [s.frac_coords for s in pystruct.sites]]
        end

        red_kcoords, _ = bzmesh_uniform(kgrid_size, kshift=kshift)
        irred_kcoords, ksymops = bzmesh_ir_wedge(kgrid_size, DFTK.symmetry_operations(lattice, atoms); kshift=kshift)

        # Try to reproduce all kcoords from irred_kcoords
        all_kcoords = Vector{Vec3{Rational{Int}}}()
        for (ik, k) in enumerate(irred_kcoords)
            append!(all_kcoords, [S * k for (S, τ) in ksymops[ik]])
        end

        # Normalize the obtained k-Points and test for equality
        red_kcoords = sort([mod.(k .* kgrid_size, kgrid_size) for k in red_kcoords])
        all_kcoords = sort([mod.(k .* kgrid_size, kgrid_size) for k in all_kcoords])
        @test all_kcoords == red_kcoords
    end

    test_reduction(silicon, [ 2,  3,  2])
    test_reduction(silicon, [ 3,  3,  3])
    test_reduction(silicon, [ 2,  3,  4])
    test_reduction(silicon, [ 9, 11, 13])

    test_reduction(silicon, [ 3,  3,  3], kshift=[1//2, 1//2, 1//2])
    test_reduction(silicon, [ 3,  3,  3], kshift=[1//2, 0, 1//2])
    test_reduction(silicon, [ 3,  3,  3], kshift=[0, 1//2, 0])

    test_reduction(silicon, [ 1,  4,  4], supercell=[2, 1, 1])
    test_reduction(silicon, [ 1,  16,  16], supercell=[4, 1, 1])

    test_reduction(magnesium, [ 2,  3,  2])
    test_reduction(magnesium, [ 3,  3,  3])
    test_reduction(magnesium, [ 2,  3,  4])
    test_reduction(magnesium, [ 9, 11, 13])
end

@testset "standardize_atoms" begin
    # Test unperturbed structure
    atoms = [ElementCoulomb(:Si) => silicon.positions]
    slattice, satoms = standardize_atoms(silicon.lattice, atoms, primitive=true)
    @test length(atoms) == 1
    @test atoms[1][1] == ElementCoulomb(:Si)
    @test length(atoms[1][2]) == 2
    @test atoms[1][2][1] == ones(3) ./ 8
    @test atoms[1][2][2] == -ones(3) ./ 8

    # Perturb structure
    plattice = silicon.lattice .+ 1e-8rand(3, 3)
    patoms = [ElementCoulomb(:Si) => [p + 1e-8rand(3) for p in silicon.positions]]
    plattice, patoms = standardize_atoms(plattice, patoms, primitive=true)

    # And check we get the usual silicon primitive cell back:
    a = plattice[1, 2]
    @test plattice == [0  a  a; a 0 a; a a 0]
    @test length(atoms) == 1
    @test atoms[1][1] == ElementCoulomb(:Si)
    @test length(atoms[1][2]) == 2
    @test atoms[1][2][1] - atoms[1][2][2] == ones(3) ./ 4
end


text:
module fkMigration
	using MAT
	using FFTW

	export demo, loadDataset, calibrate!, downsampleAndCrop, reconstruct

	const c = 299792458  # speed of light (m/s)

	demo(dataset = "teaser", resXY = 512, cropT = 1024) = reconstruct(downsampleAndCrop(calibrate!(loadDataset(dataset)...), resXY, cropT))

	function loadDataset(dataset, filename = "meas_180min.mat")
		data = matread(joinpath(dataset, filename))
		tau = data["meas"]

		data = matread(joinpath(dataset, "tof.mat"))
		calib = data["tofgrid"]

		return (tau, calib)
	end

	function calibrate!(tau, calib; wallSize = 2, resT = 32e-12)
		# calibrate so that t=0 is when light reaches the scan surface
		for x in 1:size(tau, 1), y in 1:size(tau, 2)
			delay = floor(Int, calib[x, y]*1e-12 / resT)
			tau[x, y, :] .= circshift(tau[x, y, :], -delay)
		end
		return tau
	end

	function downsampleAndCrop(tau, resXY, cropT)
		function downsample(A)
			Isize = (size(A, 1)+1)÷2
			Jsize = (size(A, 2)+1)÷2
			for i in 1:Isize, j in 1:Jsize
				I = 2*(i-1) + 1
				J = 2*(j-1) + 1
				A[i, j, :] .= 0.25 .* (
					A[min(I,   size(A, 1)), min(J,   size(A, 2)), :] .+
					A[min(I+1, size(A, 1)), min(J,   size(A, 2)), :] .+
					A[min(I+1, size(A, 1)), min(J+1, size(A, 2)), :] .+
					A[min(I,   size(A, 1)), min(J+1, size(A, 2)), :]
				)
			end
			return A[1:Isize, 1:Jsize, :]
		end

		tau = tau[:, :, 1:min(cropT, size(tau, 3))]
		while size(tau, 1) > resXY || size(tau, 2) > resXY
			tau = downsample(tau)
		end
		return tau
	end

	function reconstruct(tau; wallSize = 2, resT = 32e-12)
		(X, Y, T) = size(tau)
		Psi = zeros(eltype(tau), map(s -> 2*s, size(tau)))
		for t in 1:T
			Psi[1:X, 1:Y, t] .= sqrt.(tau[1:X, 1:Y, t])
		end

		PhiBar = fftshift(fft(Psi))

		Phi = zero(PhiBar)
		for K_x in 1:size(Phi, 1), K_y in 1:size(Phi, 2), K_z in 1:size(Phi, 3)
			Phi[K_x, K_y, K_z] = stoltInterp(PhiBar, K_x, K_y, K_z, physWidth = wallSize, physHeight = wallSize, resT = resT)
		end

		Psi = ifft(ifftshift(Phi))
		Psi  = abs2.(Psi[1:X, 1:Y, 1:T])

		return Psi
	end

	function stoltLinearInterpZ(A, x, y, z)
		zLow  = floor(Int, z)
		zHigh = ceil(Int, z)
		low   = 1 <= zLow  <= size(A, 3) ? A[x, y, zLow] : 0
		high  = 1 <= zHigh <= size(A, 3) ? A[x, y, zHigh] : 0
		return low*(zHigh-z) + high*(z-zLow)
	end

	function stoltInterp(A, x, y, z; physWidth, physHeight, resT)
		# normalize indices and convert i and j to units of k
		i = 2*(x-1)/size(A, 1) - 1
		j = 2*(y-1)/size(A, 2) - 1
		k = 2*(z-1)/size(A, 3) - 1

		i *= (size(A, 1)/2*c*resT) / (4*(physWidth/2))
		j *= (size(A, 2)/2*c*resT) / (4*(physHeight/2))

		l = sqrt(i^2 + j^2 + k^2)
		f = (l + 1)*size(A, 3)/2 + 1
		S = max(c/2 * k/max(l, eps()), 0)
		return k > 0 ? S*stoltLinearInterpZ(A, x, y, f) : 0
	end
end


text:
# tensor for classical 2-d model
const isingβc = log(1+sqrt(2))/2

"""
    tensorfromclassical(h::Matrix)

given a classical 2-body hamiltonian `h`, return the corresponding tensor
for use in e.g. `trg` for a two-dimensional square-lattice.

# Example
```julia
julia> model_tensor(Ising(),β) ≈ tensorfromclassical([β -β; -β β])

true
```
"""
function tensorfromclassical(ham::Matrix)
    wboltzmann = exp.(ham)
    q = sqrt(wboltzmann)
    ein"ij,ik,il,im -> jklm"(q,q,q,q)
end

"""
    model_tensor(::Ising,β)

return the isingtensor at inverse temperature `β` for a two-dimensional
square lattice tensor-network.
"""
function model_tensor(::Ising, β::Real)
    a = reshape(Float64[1 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 1], 2,2,2,2)
    cβ, sβ = sqrt(cosh(β)), sqrt(sinh(β))
    q = 1/sqrt(2) * [cβ+sβ cβ-sβ; cβ-sβ cβ+sβ]
    ein"abcd,ai,bj,ck,dl -> ijkl"(a,q,q,q,q)
end

"""
    mag_tensor(::Ising,β)

return the operator for the magnetisation at inverse temperature `β`
at a site in the two-dimensional ising model on a square lattice in tensor-network form.
"""
function mag_tensor(::Ising, β)
    a = reshape(Float64[1 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 -1] , 2,2,2,2)
    cβ, sβ = sqrt(cosh(β)), sqrt(sinh(β))
    q = 1/sqrt(2) * [cβ+sβ cβ-sβ; cβ-sβ cβ+sβ]
    ein"abcd,ai,bj,ck,dl -> ijkl"(a,q,q,q,q)
end

"""
    magnetisation(model<:HamiltonianModel, β, χ)

return the magnetisation of the `model` as a function of the inverse
temperature `β` and the environment bonddimension `χ` as calculated with
ctmrg. Requires that `mag_tensor` and `model_tensor` are defined for `model`.
"""
function magnetisation(model::MT, β, χ) where {MT <: HamiltonianModel}
    a = model_tensor(model,β)
    m = mag_tensor(model, β)
    rt = SquareCTMRGRuntime(a, Val(:random), χ)
    env = ctmrg(rt; tol=1e-6, maxit=100)
    corner, edge = env.corner, env.edge
    ctc  = ein"ia,ajb,bk -> ijk"(corner,edge,corner)
    env  = ein"alc,ckd,bjd,bia -> ijkl"(ctc,edge,ctc,edge)
    mag  = ein"ijkl,ijkl ->"(env,m)[]
    norm = ein"ijkl,ijkl ->"(env,a)[]

    return abs(mag/norm)
end

"""
    magofβ(::Ising,β)

return the analytical result for the magnetisation at inverse temperature
`β` for the 2d classical ising model.
"""
magofβ(::Ising, β) = β > isingβc ? (1-sinh(2*β)^-4)^(1/8) : 0.


text:
"""

# CmdStanModels

Base model specialized in:

### Method
```julia
*  SampleModel                  : StanSample.jl
*  OptimizeModel                : StanOptimize.jl
*  VariationalModel             : StanVariational.jl
*  DiagnoseModel                : StanDiagnose.jl
```
"""
abstract type CmdStanModels end

make_string = "make"

if Sys.iswindows()
    make_command =  "mingw32-" * make
end

"""
set_make_string!()

Update the name of the make command, e.g. "make" or "minw32-make".
By default make_string == "make" and on Windows make_string == "ming32-make".

Display value as `StanBase.make_string`.

$(SIGNATURES)

# Extended help

Internal, not exported.
"""
function set_make_string!(str = make_string)
    make_string = str
end

"""
make command.

$(SIGNATURES)

# Extended help

Internal, not exported.
"""
function make_command(str = make_string)
    str
end

"""

Executable path corresponding to a source file, or a model.

$(SIGNATURES)

# Extended help

Internal, not exported.
"""
function executable_path(source_path::AbstractString)
    if Sys.iswindows()
        source_path =  source_path * ".exe"
    end
    source_path = source_path
end

executable_path(m::T) where T <: CmdStanModels =
    executable_path(m.tmpdir)

"""

Error thrown when a Stan model fails to compile.

$(TYPEDEF)

# Extended help

Accessing fields directly is part of the API.

$(FIELDS)
"""
struct StanModelError <: Exception
    name::String
    message::String
end

function Base.showerror(io::IO, e::StanModelError)
    print(io, "\nError when compiling SampleModel ", e.name, ":\n",
          e.message)
end

"""

Ensure that a compiled model executable exists, and return its path.

$(SIGNATURES)

# Extended help

If compilation fails, a `StanModelError` is returned instead.

Internal, not exported.
"""
function ensure_executable(m::T) where T <: CmdStanModels
    @unpack cmdstan_home, exec_path = m
    error_output = IOBuffer()
    is_ok = cd(cmdstan_home) do
        success(pipeline(`$(make_command()) -f $(cmdstan_home)/makefile -C $(cmdstan_home) $(exec_path)`;
                         stderr = error_output))
    end
    if is_ok
        exec_path
    else
        throw(StanModelError(m.name, String(take!(error_output))))
    end
end

"""

Compile a model, throwing an error if it failed.

$(SIGNATURES)
"""
function stan_compile(m::T) where T <: CmdStanModels
    ensure_executable(m)
    nothing
end

"""

Default `output_base` data files, in tmpdir.

$(SIGNATURES)

# Extended help

Internal, not exported.
"""
data_file_path(output_base::AbstractString, id::Int) = output_base * "_data_$(id).R"

"""

Default `output_base` init files, in tmpdir.

$(SIGNATURES)

# Extended help

Internal, not exported.
"""
init_file_path(output_base::AbstractString, id::Int) =
  output_base * "_init_$(id).R"

"""

Default `output_base` chain files, in tmpdir.

$(SIGNATURES)

# Extended help

Internal, not exported.
"""
sample_file_path(output_base::AbstractString, id::Int) =
  output_base * "_chain_$(id).csv"

"""

Default `output_base` for generated_quatities files, in tmpdir.

$(SIGNATURES)

# Extended help

Internal, not exported.
"""
generated_quantities_file_path(output_base::AbstractString, id::Int) =
  output_base * "_generated_quantities_$(id).csv"

"""

Default `output_base` log files, in tmpdir.

$(SIGNATURES)

# Extended help

Internal, not exported.
"""
log_file_path(output_base::AbstractString, id::Int) =
  output_base * "_log_$(id).log"

"""

Default `output_base` diagnostic files, in tmpdir.

$(SIGNATURES)

# Extended help

Internal, not exported.
"""
diagnostic_file_path(output_base::AbstractString, id::Int) =
  output_base * "_diagnostic_$(id).csv"

"""
Helper function for the (deprecated) diagnostics file generation.

$(SIGNATURES)

# Extended help

I am not aware the diagnostic files contain other info then the regular .csv files.
Currently I have not disabled this functionality. Please let me know if this
feature should be included/enabled.
"""
function setup_diagnostics(m, num_chains)
  for i in 1:num_chains
    append!(m.diagnostic_file, [m.output_base*"_diagnostic_$i.log"])
  end
end


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

# Eigendecomposition
"""
    Eigen <: Factorization

Matrix factorization type of the eigenvalue/spectral decomposition of a square
matrix `A`. This is the return type of [`eigen`](@ref), the corresponding matrix
factorization function.

If `F::Eigen` is the factorization object, the eigenvalues can be obtained via
`F.values` and the eigenvectors as the columns of the matrix `F.vectors`.
(The `k`th eigenvector can be obtained from the slice `F.vectors[:, k]`.)

Iterating the decomposition produces the components `F.values` and `F.vectors`.

# Examples
```jldoctest
julia> F = eigen([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
3-element Vector{Float64}:
  1.0
  3.0
 18.0
vectors:
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> F.values
3-element Vector{Float64}:
  1.0
  3.0
 18.0

julia> F.vectors
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> vals, vecs = F; # destructuring via iteration

julia> vals == F.values && vecs == F.vectors
true
```
"""
struct Eigen{T,V,S<:AbstractMatrix,U<:AbstractVector} <: Factorization{T}
    values::U
    vectors::S
    Eigen{T,V,S,U}(values::AbstractVector{V}, vectors::AbstractMatrix{T}) where {T,V,S,U} =
        new(values, vectors)
end
Eigen(values::AbstractVector{V}, vectors::AbstractMatrix{T}) where {T,V} =
    Eigen{T,V,typeof(vectors),typeof(values)}(values, vectors)

# Generalized eigenvalue problem.
"""
    GeneralizedEigen <: Factorization

Matrix factorization type of the generalized eigenvalue/spectral decomposition of
`A` and `B`. This is the return type of [`eigen`](@ref), the corresponding
matrix factorization function, when called with two matrix arguments.

If `F::GeneralizedEigen` is the factorization object, the eigenvalues can be obtained via
`F.values` and the eigenvectors as the columns of the matrix `F.vectors`.
(The `k`th eigenvector can be obtained from the slice `F.vectors[:, k]`.)

Iterating the decomposition produces the components `F.values` and `F.vectors`.

# Examples
```jldoctest
julia> A = [1 0; 0 -1]
2×2 Matrix{Int64}:
 1   0
 0  -1

julia> B = [0 1; 1 0]
2×2 Matrix{Int64}:
 0  1
 1  0

julia> F = eigen(A, B)
GeneralizedEigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}
values:
2-element Vector{ComplexF64}:
 0.0 - 1.0im
 0.0 + 1.0im
vectors:
2×2 Matrix{ComplexF64}:
  0.0+1.0im   0.0-1.0im
 -1.0+0.0im  -1.0-0.0im

julia> F.values
2-element Vector{ComplexF64}:
 0.0 - 1.0im
 0.0 + 1.0im

julia> F.vectors
2×2 Matrix{ComplexF64}:
  0.0+1.0im   0.0-1.0im
 -1.0+0.0im  -1.0-0.0im

julia> vals, vecs = F; # destructuring via iteration

julia> vals == F.values && vecs == F.vectors
true
```
"""
struct GeneralizedEigen{T,V,S<:AbstractMatrix,U<:AbstractVector} <: Factorization{T}
    values::U
    vectors::S
    GeneralizedEigen{T,V,S,U}(values::AbstractVector{V}, vectors::AbstractMatrix{T}) where {T,V,S,U} =
        new(values, vectors)
end
GeneralizedEigen(values::AbstractVector{V}, vectors::AbstractMatrix{T}) where {T,V} =
    GeneralizedEigen{T,V,typeof(vectors),typeof(values)}(values, vectors)

# iteration for destructuring into components
Base.iterate(S::Union{Eigen,GeneralizedEigen}) = (S.values, Val(:vectors))
Base.iterate(S::Union{Eigen,GeneralizedEigen}, ::Val{:vectors}) = (S.vectors, Val(:done))
Base.iterate(S::Union{Eigen,GeneralizedEigen}, ::Val{:done}) = nothing

isposdef(A::Union{Eigen,GeneralizedEigen}) = isreal(A.values) && all(x -> x > 0, A.values)

# pick a canonical ordering to avoid returning eigenvalues in "random" order
# as is the LAPACK default (for complex λ — LAPACK sorts by λ for the Hermitian/Symmetric case)
eigsortby(λ::Real) = λ
eigsortby(λ::Complex) = (real(λ),imag(λ))
function sorteig!(λ::AbstractVector, X::AbstractMatrix, sortby::Union{Function,Nothing}=eigsortby)
    if sortby !== nothing && !issorted(λ, by=sortby)
        p = sortperm(λ; alg=QuickSort, by=sortby)
        permute!(λ, p)
        Base.permutecols!!(X, p)
    end
    return λ, X
end
sorteig!(λ::AbstractVector, sortby::Union{Function,Nothing}=eigsortby) = sortby === nothing ? λ : sort!(λ, by=sortby)

"""
    eigen!(A, [B]; permute, scale, sortby)

Same as [`eigen`](@ref), but saves space by overwriting the input `A` (and
`B`), instead of creating a copy.
"""
function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T<:BlasReal
    n = size(A, 2)
    n == 0 && return Eigen(zeros(T, 0), zeros(T, 0, 0))
    issymmetric(A) && return eigen!(Symmetric(A), sortby=sortby)
    A, WR, WI, VL, VR, _ = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)
    iszero(WI) && return Eigen(sorteig!(WR, VR, sortby)...)
    evec = zeros(Complex{T}, n, n)
    j = 1
    while j <= n
        if WI[j] == 0
            evec[:,j] = view(VR, :, j)
        else
            for i = 1:n
                evec[i,j]   = VR[i,j] + im*VR[i,j+1]
                evec[i,j+1] = VR[i,j] - im*VR[i,j+1]
            end
            j += 1
        end
        j += 1
    end
    return Eigen(sorteig!(complex.(WR, WI), evec, sortby)...)
end

function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T<:BlasComplex
    n = size(A, 2)
    n == 0 && return Eigen(zeros(T, 0), zeros(T, 0, 0))
    ishermitian(A) && return eigen!(Hermitian(A), sortby=sortby)
    eval, evec = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)[[2,4]]
    return Eigen(sorteig!(eval, evec, sortby)...)
end

"""
    eigen(A; permute::Bool=true, scale::Bool=true, sortby) -> Eigen

Computes the eigenvalue decomposition of `A`, returning an [`Eigen`](@ref) factorization object `F`
which contains the eigenvalues in `F.values` and the eigenvectors in the columns of the
matrix `F.vectors`. (The `k`th eigenvector can be obtained from the slice `F.vectors[:, k]`.)

Iterating the decomposition produces the components `F.values` and `F.vectors`.

The following functions are available for `Eigen` objects: [`inv`](@ref), [`det`](@ref), and [`isposdef`](@ref).

For general nonsymmetric matrices it is possible to specify how the matrix is balanced
before the eigenvector calculation. The option `permute=true` permutes the matrix to become
closer to upper triangular, and `scale=true` scales the matrix by its diagonal elements to
make rows and columns more equal in norm. The default is `true` for both options.

By default, the eigenvalues and vectors are sorted lexicographically by `(real(λ),imag(λ))`.
A different comparison function `by(λ)` can be passed to `sortby`, or you can pass
`sortby=nothing` to leave the eigenvalues in an arbitrary order.   Some special matrix types
(e.g. [`Diagonal`](@ref) or [`SymTridiagonal`](@ref)) may implement their own sorting convention and not
accept a `sortby` keyword.

# Examples
```jldoctest
julia> F = eigen([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
3-element Vector{Float64}:
  1.0
  3.0
 18.0
vectors:
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> F.values
3-element Vector{Float64}:
  1.0
  3.0
 18.0

julia> F.vectors
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> vals, vecs = F; # destructuring via iteration

julia> vals == F.values && vecs == F.vectors
true
```
"""
function eigen(A::AbstractMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T
    AA = copy_oftype(A, eigtype(T))
    isdiag(AA) && return eigen(Diagonal(AA); permute=permute, scale=scale, sortby=sortby)
    return eigen!(AA; permute=permute, scale=scale, sortby=sortby)
end
eigen(x::Number) = Eigen([x], fill(one(x), 1, 1))

"""
    eigvecs(A; permute::Bool=true, scale::Bool=true, `sortby`) -> Matrix

Return a matrix `M` whose columns are the eigenvectors of `A`. (The `k`th eigenvector can
be obtained from the slice `M[:, k]`.) The `permute`, `scale`, and `sortby` keywords are the same as
for [`eigen`](@ref).

# Examples
```jldoctest
julia> eigvecs([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0
```
"""
eigvecs(A::Union{Number, AbstractMatrix}; kws...) =
    eigvecs(eigen(A; kws...))
eigvecs(F::Union{Eigen, GeneralizedEigen}) = F.vectors

eigvals(F::Union{Eigen, GeneralizedEigen}) = F.values

"""
    eigvals!(A; permute::Bool=true, scale::Bool=true, sortby) -> values

Same as [`eigvals`](@ref), but saves space by overwriting the input `A`, instead of creating a copy.
The `permute`, `scale`, and `sortby` keywords are the same as for [`eigen`](@ref).

!!! note
    The input matrix `A` will not contain its eigenvalues after `eigvals!` is
    called on it - `A` is used as a workspace.

# Examples
```jldoctest
julia> A = [1. 2.; 3. 4.]
2×2 Matrix{Float64}:
 1.0  2.0
 3.0  4.0

julia> eigvals!(A)
2-element Vector{Float64}:
 -0.3722813232690143
  5.372281323269014

julia> A
2×2 Matrix{Float64}:
 -0.372281  -1.0
  0.0        5.37228
```
"""
function eigvals!(A::StridedMatrix{<:BlasReal}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby)
    issymmetric(A) && return sorteig!(eigvals!(Symmetric(A)), sortby)
    _, valsre, valsim, _ = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'N', 'N', A)
    return sorteig!(iszero(valsim) ? valsre : complex.(valsre, valsim), sortby)
end
function eigvals!(A::StridedMatrix{<:BlasComplex}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby)
    ishermitian(A) && return sorteig!(eigvals(Hermitian(A)), sortby)
    return sorteig!(LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'N', 'N', A)[2], sortby)
end

# promotion type to use for eigenvalues of a Matrix{T}
eigtype(T) = promote_type(Float32, typeof(zero(T)/sqrt(abs2(one(T)))))

"""
    eigvals(A; permute::Bool=true, scale::Bool=true, sortby) -> values

Return the eigenvalues of `A`.

For general non-symmetric matrices it is possible to specify how the matrix is balanced
before the eigenvalue calculation. The `permute`, `scale`, and `sortby` keywords are
the same as for [`eigen!`](@ref).

# Examples
```jldoctest
julia> diag_matrix = [1 0; 0 4]
2×2 Matrix{Int64}:
 1  0
 0  4

julia> eigvals(diag_matrix)
2-element Vector{Float64}:
 1.0
 4.0
```
"""
eigvals(A::AbstractMatrix{T}; kws...) where T =
    eigvals!(copy_oftype(A, eigtype(T)); kws...)

"""
For a scalar input, `eigvals` will return a scalar.

# Example
```jldoctest
julia> eigvals(-2)
-2
```
"""
eigvals(x::Number; kwargs...) = imag(x) == 0 ? real(x) : x

"""
    eigmax(A; permute::Bool=true, scale::Bool=true)

Return the largest eigenvalue of `A`.
The option `permute=true` permutes the matrix to become
closer to upper triangular, and `scale=true` scales the matrix by its diagonal elements to
make rows and columns more equal in norm.
Note that if the eigenvalues of `A` are complex,
this method will fail, since complex numbers cannot
be sorted.

# Examples
```jldoctest
julia> A = [0 im; -im 0]
2×2 Matrix{Complex{Int64}}:
 0+0im  0+1im
 0-1im  0+0im

julia> eigmax(A)
1.0

julia> A = [0 im; -1 0]
2×2 Matrix{Complex{Int64}}:
  0+0im  0+1im
 -1+0im  0+0im

julia> eigmax(A)
ERROR: DomainError with Complex{Int64}[0+0im 0+1im; -1+0im 0+0im]:
`A` cannot have complex eigenvalues.
Stacktrace:
[...]
```
"""
function eigmax(A::Union{Number, AbstractMatrix}; permute::Bool=true, scale::Bool=true)
    v = eigvals(A, permute = permute, scale = scale)
    if eltype(v)<:Complex
        throw(DomainError(A, "`A` cannot have complex eigenvalues."))
    end
    maximum(v)
end

"""
    eigmin(A; permute::Bool=true, scale::Bool=true)

Return the smallest eigenvalue of `A`.
The option `permute=true` permutes the matrix to become
closer to upper triangular, and `scale=true` scales the matrix by its diagonal elements to
make rows and columns more equal in norm.
Note that if the eigenvalues of `A` are complex,
this method will fail, since complex numbers cannot
be sorted.

# Examples
```jldoctest
julia> A = [0 im; -im 0]
2×2 Matrix{Complex{Int64}}:
 0+0im  0+1im
 0-1im  0+0im

julia> eigmin(A)
-1.0

julia> A = [0 im; -1 0]
2×2 Matrix{Complex{Int64}}:
  0+0im  0+1im
 -1+0im  0+0im

julia> eigmin(A)
ERROR: DomainError with Complex{Int64}[0+0im 0+1im; -1+0im 0+0im]:
`A` cannot have complex eigenvalues.
Stacktrace:
[...]
```
"""
function eigmin(A::Union{Number, AbstractMatrix};
                permute::Bool=true, scale::Bool=true)
    v = eigvals(A, permute = permute, scale = scale)
    if eltype(v)<:Complex
        throw(DomainError(A, "`A` cannot have complex eigenvalues."))
    end
    minimum(v)
end

inv(A::Eigen) = A.vectors * inv(Diagonal(A.values)) / A.vectors
det(A::Eigen) = prod(A.values)

# Generalized eigenproblem
function eigen!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function,Nothing}=eigsortby) where T<:BlasReal
    issymmetric(A) && isposdef(B) && return eigen!(Symmetric(A), Symmetric(B), sortby=sortby)
    n = size(A, 1)
    alphar, alphai, beta, _, vr = LAPACK.ggev!('N', 'V', A, B)
    iszero(alphai) && return GeneralizedEigen(sorteig!(alphar ./ beta, vr, sortby)...)

    vecs = zeros(Complex{T}, n, n)
    j = 1
    while j <= n
        if alphai[j] == 0
            vecs[:,j] = view(vr, :, j)
        else
            for i = 1:n
                vecs[i,j  ] = vr[i,j] + im*vr[i,j+1]
                vecs[i,j+1] = vr[i,j] - im*vr[i,j+1]
            end
            j += 1
        end
        j += 1
    end
    return GeneralizedEigen(sorteig!(complex.(alphar, alphai)./beta, vecs, sortby)...)
end

function eigen!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function,Nothing}=eigsortby) where T<:BlasComplex
    ishermitian(A) && isposdef(B) && return eigen!(Hermitian(A), Hermitian(B), sortby=sortby)
    alpha, beta, _, vr = LAPACK.ggev!('N', 'V', A, B)
    return GeneralizedEigen(sorteig!(alpha./beta, vr, sortby)...)
end

"""
    eigen(A, B) -> GeneralizedEigen

Computes the generalized eigenvalue decomposition of `A` and `B`, returning a
[`GeneralizedEigen`](@ref) factorization object `F` which contains the generalized eigenvalues in
`F.values` and the generalized eigenvectors in the columns of the matrix `F.vectors`.
(The `k`th generalized eigenvector can be obtained from the slice `F.vectors[:, k]`.)

Iterating the decomposition produces the components `F.values` and `F.vectors`.

Any keyword arguments passed to `eigen` are passed through to the lower-level
[`eigen!`](@ref) function.

# Examples
```jldoctest
julia> A = [1 0; 0 -1]
2×2 Matrix{Int64}:
 1   0
 0  -1

julia> B = [0 1; 1 0]
2×2 Matrix{Int64}:
 0  1
 1  0

julia> F = eigen(A, B);

julia> F.values
2-element Vector{ComplexF64}:
 0.0 - 1.0im
 0.0 + 1.0im

julia> F.vectors
2×2 Matrix{ComplexF64}:
  0.0+1.0im   0.0-1.0im
 -1.0+0.0im  -1.0-0.0im

julia> vals, vecs = F; # destructuring via iteration

julia> vals == F.values && vecs == F.vectors
true
```
"""
function eigen(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}; kws...) where {TA,TB}
    S = promote_type(eigtype(TA),TB)
    eigen!(copy_oftype(A, S), copy_oftype(B, S); kws...)
end

eigen(A::Number, B::Number) = eigen(fill(A,1,1), fill(B,1,1))

"""
    eigvals!(A, B; sortby) -> values

Same as [`eigvals`](@ref), but saves space by overwriting the input `A` (and `B`),
instead of creating copies.

!!! note
    The input matrices `A` and `B` will not contain their eigenvalues after
    `eigvals!` is called. They are used as workspaces.

# Examples
```jldoctest
julia> A = [1. 0.; 0. -1.]
2×2 Matrix{Float64}:
 1.0   0.0
 0.0  -1.0

julia> B = [0. 1.; 1. 0.]
2×2 Matrix{Float64}:
 0.0  1.0
 1.0  0.0

julia> eigvals!(A, B)
2-element Vector{ComplexF64}:
 0.0 - 1.0im
 0.0 + 1.0im

julia> A
2×2 Matrix{Float64}:
 -0.0  -1.0
  1.0  -0.0

julia> B
2×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0
```
"""
function eigvals!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function,Nothing}=eigsortby) where T<:BlasReal
    issymmetric(A) && isposdef(B) && return sorteig!(eigvals!(Symmetric(A), Symmetric(B)), sortby)
    alphar, alphai, beta, vl, vr = LAPACK.ggev!('N', 'N', A, B)
    return sorteig!((iszero(alphai) ? alphar : complex.(alphar, alphai))./beta, sortby)
end
function eigvals!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function,Nothing}=eigsortby) where T<:BlasComplex
    ishermitian(A) && isposdef(B) && return sorteig!(eigvals!(Hermitian(A), Hermitian(B)), sortby)
    alpha, beta, vl, vr = LAPACK.ggev!('N', 'N', A, B)
    return sorteig!(alpha./beta, sortby)
end

"""
    eigvals(A, B) -> values

Computes the generalized eigenvalues of `A` and `B`.

# Examples
```jldoctest
julia> A = [1 0; 0 -1]
2×2 Matrix{Int64}:
 1   0
 0  -1

julia> B = [0 1; 1 0]
2×2 Matrix{Int64}:
 0  1
 1  0

julia> eigvals(A,B)
2-element Vector{ComplexF64}:
 0.0 - 1.0im
 0.0 + 1.0im
```
"""
function eigvals(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}; kws...) where {TA,TB}
    S = promote_type(eigtype(TA),TB)
    return eigvals!(copy_oftype(A, S), copy_oftype(B, S); kws...)
end

"""
    eigvecs(A, B) -> Matrix

Return a matrix `M` whose columns are the generalized eigenvectors of `A` and `B`. (The `k`th eigenvector can
be obtained from the slice `M[:, k]`.)

# Examples
```jldoctest
julia> A = [1 0; 0 -1]
2×2 Matrix{Int64}:
 1   0
 0  -1

julia> B = [0 1; 1 0]
2×2 Matrix{Int64}:
 0  1
 1  0

julia> eigvecs(A, B)
2×2 Matrix{ComplexF64}:
  0.0+1.0im   0.0-1.0im
 -1.0+0.0im  -1.0-0.0im
```
"""
eigvecs(A::AbstractMatrix, B::AbstractMatrix; kws...) = eigvecs(eigen(A, B; kws...))

function show(io::IO, mime::MIME{Symbol("text/plain")}, F::Union{Eigen,GeneralizedEigen})
    summary(io, F); println(io)
    println(io, "values:")
    show(io, mime, F.values)
    println(io, "\nvectors:")
    show(io, mime, F.vectors)
end

# Conversion methods

## Can we determine the source/result is Real?  This is not stored in the type Eigen
AbstractMatrix(F::Eigen) = F.vectors * Diagonal(F.values) / F.vectors
AbstractArray(F::Eigen) = AbstractMatrix(F)
Matrix(F::Eigen) = Array(AbstractArray(F))
Array(F::Eigen) = Matrix(F)

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using NumericalMethodsforEngineers, DataFrames, Plots

@sym begin
  LagrangePolynomial(xi_, yi_) := Module([N, sum, num, den],
    begin
      N = Length(xi)
      sum = 0
      For( i=1, i <= N, Increment(i),
        begin
          num = 1
          den = 1
          For( j=1, j <= N, Increment(j),
            begin
              If(j != i, num *= (x-xi[j]))
              If(j != i, den *= (xi[i]-xi[j]))
            end
          )
          sum += yi[i] * num/den
        end
      ),
      Return(Simplify(sum))
    end
  )
end

println()
@sym begin
  xi = [1, 3, 6]
  yi = [1, 5, 10]
  # LagrangePolynomial defined in ../../src/ch07/sym_fcns.jl
  Q(x_) = LagrangePolynomial(xi, yi)
  SetJ(qs, ToString(Q(x)))
  Print("Q(x) = $(Q(x))\n")
end

@sym Print("Q(4.5) = $(N(Q(4.5)))\n");
println("\n")

ProjDir = dirname(@__FILE__)
cd(ProjDir) #do

  x = [1.0, 3.0, 6.0]
  y = [1.0, 5.0, 10.0]
  xi = [2.0, 4.5]
  (dfin, dfxi) = lagrangianpolynomial(length(x), x, y, xi)

  dfin |> display
  println()

  dfxi |> display
  println()

  @eval f(x) = $(Meta.parse(qs))

  xint = 1:1:10
  p = plot(xint, f.(xint), line=(:path, 1), label="interpolated curve")
  scatter!(p, dfin[:x], dfin[:y], marker=(:circle, 4), label="input points", color=:blue)
  scatter!(p, dfxi[:xi], dfxi[:yi], marker=(:star, 8), color=:red, label="interpolated points")

  plot(p)
  savefig("ex.5.1_symata.png")
  #gui()

  #end


text:
using Base: splat

test = """
be cfbegad cbdgef fgaecd cgeb fdcge agebfd fecdb fabcd edb | fdgacbe cefdb cefbgd gcbe
edbfga begcd cbg gc gcadebf fbgde acbgfd abcde gfcbed gfec | fcgedb cgb dgebacf gc
fgaebd cg bdaec gdafb agbcfd gdcbef bgcad gfac gcb cdgabef | cg cg fdcagb cbg
fbegcd cbd adcefb dageb afcb bc aefdc ecdab fgdeca fcdbega | efabcd cedba gadfec cb
aecbfdg fbg gf bafeg dbefa fcge gcbea fcaegb dgceab fcbdga | gecf egdcabf bgf bfgea
fgeab ca afcebg bdacfeg cfaedg gcfdb baec bfadeg bafgc acf | gebdcfa ecba ca fadegcb
dbcfg fgd bdegcaf fgec aegbdf ecdfab fbedc dacgb gdcebf gf | cefg dcbef fcge gbcadfe
bdfegc cbegaf gecbf dfcage bdacg ed bedf ced adcbefg gebcd | ed bcgafe cdgba cbgef
egadfb cdbfeg cegd fecab cgb gbdefca cg fgcdab egfdb bfceg | gbdfcae bgc cg cgb
gcafb gcf dcaebfg ecagb gf abcdeg gaef cafbge fdbac fegbdc | fgae cfgab fg bagce
"""

function parse_input(input::AbstractString)::Vector{Tuple}
    return map(split(input, "\n"; keepempty = false)) do line
        return split.(split(line, " | "), ' '; keepempty = false) |> Tuple
    end
end

function q1(specs::AbstractVector)::Int
    return sum(specs) do (_, display)
        return count(∈([2, 3, 4, 7]) ∘ length, display)
    end
end

function decode(digits::AbstractVector, display::AbstractVector)::Int
    all_digits = [digits; display]
    l(n...) = s -> length(s) ∈ n
    d(n...) = filter(l(n...), all_digits)
    i(v...) = intersect(d(v...)...)
    u(v...) = union(d(v...)...)
    Δ(vs...) = setdiff(vs...)

    # 1: 2segs, 7: 3segs, 4: 4segs, 2, 3, 5: 5segs, 0, 6, 9: 6segs, 8: 7segs
    segments = [Δ(i(3, 5, 6, 7), u(2, 4)),    # in all but 1, 4
                Δ(i(4, 6), u(2, 3), i(5, 7)), # in 0, 6, 9, not in 1, 7
                Δ(i(2, 3, 4, 7), i(5, 6)),    # in all but 5, 6
                Δ(i(4, 5, 7), u(2, 3), i(6)), # in 2, 3, 4, 5, 8, not in 1, 7
                Δ(i(7), u(2, 3, 4), i(5, 6)), # in 8, not in 1, 4, 7
                Δ(i(2, 3, 4, 6, 7), i(5)),    # in all but 2
                Δ(i(5, 6, 7), i(2, 3, 4))]    # in all but 1, 4, 7
    # resolve remaining unambiguities
    segments = map(segments) do seg
        length(seg) == 1 && return only(seg)
        unambiguous = union(filter(l(1), segments)...)
        return Δ(seg, unambiguous) |> only
    end

    s(v...) = segments[[v...]] |> sort! |> v -> join(v, "")
    mapping = Dict(s(1, 2, 3, 5, 6, 7)    => 0,
                   s(3, 6)                => 1,
                   s(1, 3, 4, 5, 7)       => 2,
                   s(1, 3, 4, 6, 7)       => 3,
                   s(2, 3, 4, 6)          => 4,
                   s(1, 2, 4, 6, 7)       => 5,
                   s(1, 2, 4, 5, 6, 7)    => 6,
                   s(1, 3, 6)             => 7,
                   s(1, 2, 3, 4, 5, 6, 7) => 8,
                   s(1, 2, 3, 4, 6, 7)    => 9)

    v(d) = mapping[join(sort(split(d, "")), "")]
    return only([1000 100 10 1] * [v.(display)...])
end

q2(specs::AbstractVector)::Int = sum(splat(decode), specs)

let specs = parse_input(test)
    @assert q1(specs) == 26
    @assert q2(specs) == 61229
end

specs = parse_input(read("day08.in", String))
println("Q1: ", q1(specs))
println("Q2: ", q2(specs))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
# R2_IRI_assRCS_kfold_include.jl

include("../tcm_code/tcm_estimateRanking_CS.jl")
include("../tcm_code/tcm_predictRanking.jl")

include("../tcm_code/tcm_estimateRanking_kcv.jl")

include("../tcm_code/tcm_evaluateKL.jl")

include("../tcm_code/tcm_rankingStats.jl")

include("R2_IRI_assRCS_kfold_helperfn.jl")

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
function covcor_axes(old_axes::NTuple{2,Any}, new_indices::NTuple{2,Any}, dim::Int)
    if dim === 1
        return (
            unsafe_reconstruct(last(old_axes), first(new_indices)),
            StaticRanges.resize_last(last(old_axes), last(new_indices))
        )
    elseif dim === 2
        return (
            StaticRanges.resize_last(first(old_axes), first(new_indices)),
            unsafe_reconstruct(first(old_axes), last(new_indices))
        )
    else
        return (
            StaticRanges.resize_last(first(old_axes), first(new_indices)),
            StaticRanges.resize_last(last(old_axes), last(new_indices))
        )
    end
end

for fun in (:cor, :cov)
    fun_doc = """
        $fun(x::AxisMatrix; dims=1, kwargs...)

    Performs `$fun` on the parent matrix of `x` and reconstructs a similar type
    with the appropriate axes.

    ## Examples
    ```jldoctest
    julia> using AxisIndices, Statistics

    julia> A = AxisArray([1 2 3; 4 5 6; 7 8 9], ["a", "b", "c"], [:one, :two, :three]);

    julia> keys.(axes($fun(A, dims = 2)))
    (["a", "b", "c"], ["a", "b", "c"])

    julia> keys.(axes($fun(A, dims = 1)))
    ([:one, :two, :three], [:one, :two, :three])

    ```
    """
    @eval begin
        @doc $fun_doc
        function Statistics.$fun(x::AxisArray{T,2}; dims=1, kwargs...) where {T}
            p = Statistics.$fun(parent(x); dims=dims, kwargs...)
            return AxisArray(p, covcor_axes(axes(x), axes(p), dims); checks=NoChecks)
        end
    end
end

###
### 🐇/*
###
#
#    matmul_axes(a, b, p) -> Tuple
#
# The arrays from `a * b = p`, where `p` will be the new parent array of an AxisArray.

matmul_axes(a::Tuple, b::Tuple, p::Tuple) = _matmul_axes(a, b, p)
matmul_axes(a::Tuple, b::Tuple, p::Tuple{}) = ()

function _matmul_axes(a::Tuple{Any}, b::Tuple{Any,Any}, p::Tuple{Any,Any})
    return (
        _matmul_unsafe_reconstruct(first(a), first(p)),
        _matmul_unsafe_reconstruct(last(b), last(p))
    )
end

function _matmul_axes(a::Tuple{Any,Any}, b::Tuple{Any,Any}, p::Tuple{Any,Any})
    return (
        _matmul_unsafe_reconstruct(first(a), first(p)),
        _matmul_unsafe_reconstruct(last(b), last(p))
    )
end

function _matmul_axes(a::Tuple{Any,Any}, b::Tuple{Any}, p::Tuple{Any})
    return (_matmul_unsafe_reconstruct(first(a), first(p)),)
end

_matmul_unsafe_reconstruct(axis::AbstractAxis, inds) = unsafe_reconstruct(axis, inds; keys=keys(axis))
_matmul_unsafe_reconstruct(axis, inds) = SimpleAxis(inds)

_matmul(p, axs::Tuple) = AxisArray(p, axs)
_matmul(p, axs::Tuple{}) = p

function Base.:*(a::AxisMatrix, b::AxisMatrix)
    p = *(parent(a), parent(b))
    return _matmul(p, matmul_axes(axes(a), axes(b), axes(p)))
end
function Base.:*(a::AxisMatrix, b::AxisVector)
    p = *(parent(a), parent(b))
    return _matmul(p, matmul_axes(axes(a), axes(b), axes(p)))
end
function Base.:*(a::AxisVector, b::AxisMatrix)
    p = *(parent(a), parent(b))
    return _matmul(p, matmul_axes(axes(a), axes(b), axes(p)))
end

#Base.:*(x::AxisVector, y::AxisVector) = throw(MethodError(*, (x, y)))
function Base.:*(a::AxisVector, b::AxisVector)
    p = *(parent(a), parent(b))
    return _matmul(p, matmul_axes(axes(a), axes(b), axes(p)))
end

macro declare_axis_array_matmul_left(X, Y)
    esc(quote
        function Base.:*(x::$X, y::$Y)
            p = *(parent(x), y)
            return _matmul(p, matmul_axes(axes(x), axes(y), axes(p)))
        end
    end)
end

macro declare_axis_array_matmul_right(X, Y)
    esc(quote
        function Base.:*(x::$X, y::$Y)
            p = *(x, parent(y))
            return _matmul(p, matmul_axes(axes(x), axes(y), axes(p)))
        end
    end)
end

#     @declare_axisarray_matmul(T)
#
# This is used to quickly resolve ambiguities with `*`.
# It simply defines a method for `*(::T, ::AxisArray)` and
# `*(::AxisArray, ::T)`, which unwrap the `AxisArray` for subsequent
# `*` methods and than rewraps the appropriate axes as an axis array. It isn't
# exported and doesn't have an official docstring b/c it's unlikely to be useful
# outside the internals of this package.
macro declare_axisarray_matmul(T)
    esc(quote
        @declare_axis_array_matmul_left(AxisMatrix, $T)
        @declare_axis_array_matmul_left(AxisVector, $T)
        @declare_axis_array_matmul_right($T, AxisMatrix)
        @declare_axis_array_matmul_right($T, AxisVector)
    end)
end

@declare_axisarray_matmul AbstractVector

@declare_axisarray_matmul AbstractMatrix

@declare_axisarray_matmul Diagonal

@declare_axis_array_matmul_left(AxisMatrix, Transpose{<:Any,<:LinearAlgebra.RealHermSymComplexSym})
@declare_axis_array_matmul_left(AxisMatrix, Transpose{<:Any,<:LinearAlgebra.AbstractTriangular})
@declare_axis_array_matmul_left(AxisMatrix, Adjoint{<:Any,<:LinearAlgebra.RealHermSymComplexHerm})
@declare_axis_array_matmul_left(AxisMatrix, LinearAlgebra.AbstractTriangular)
@declare_axis_array_matmul_left(AxisVector, LinearAlgebra.AdjOrTransAbsVec)
@declare_axis_array_matmul_left(AxisVector, Transpose{<:Any,<:AbstractMatrix})
@declare_axis_array_matmul_left(AxisVector, Adjoint{<:Any,<:AbstractMatrix})

@declare_axis_array_matmul_right(LinearAlgebra.AdjointAbsVec{<:Number}, AxisVector{<:Number})
@declare_axis_array_matmul_right(LinearAlgebra.AdjOrTransAbsVec, AxisVector)
@declare_axis_array_matmul_right(LinearAlgebra.TransposeAbsVec, AxisMatrix)

for T in (AxisVector, AxisMatrix)
    @eval begin
        @declare_axis_array_matmul_right(LinearAlgebra.AbstractTriangular, $T)
        @declare_axis_array_matmul_right(Adjoint{<:Any,<:LinearAlgebra.AbstractTriangular}, $T)
        @declare_axis_array_matmul_right(Transpose{<:Any,<:LinearAlgebra.AbstractTriangular}, $T)
        @declare_axis_array_matmul_right(Adjoint{<:Any,<:LinearAlgebra.RealHermSymComplexHerm}, $T)
        @declare_axis_array_matmul_right(LinearAlgebra.AdjointAbsVec, $T)
        @declare_axis_array_matmul_right(Transpose{<:Any,<:LinearAlgebra.RealHermSymComplexSym}, $T)

        @declare_axis_array_matmul_left($T, Adjoint{<:Any,<:LinearAlgebra.AbstractTriangular})
        @declare_axis_array_matmul_left($T, Adjoint{<:Any,<:LinearAlgebra.AbstractRotation})
    end
end
function Base.:*(a::Transpose{<:Any,<:AbstractMatrix{T}}, b::AxisVector{S}) where {T,S}
    p = *(a, parent(b))
    return _matmul(p, matmul_axes(axes(a), axes(b), axes(p)))
end

function Base.:*(x::LinearAlgebra.TransposeAbsVec{T}, y::AxisVector{T}) where {T<:Real}
    p = *(x, parent(y))
    return _matmul(p, matmul_axes(axes(a), axes(b), axes(p)))
end

function Base.:*(a::Adjoint{<:Any,<:AbstractMatrix{T}}, b::AxisVector{S}) where {T,S}
    p = *(a, parent(b))
    return _matmul(p, matmul_axes(axes(a), axes(b), axes(p)))
end
function Base.:*(a::Adjoint{T,<:AbstractArray{T,1}}, b::AxisVector{T}) where {T}
    p = *(a, parent(b))
    return _matmul(p, matmul_axes(axes(a), axes(b), axes(p)))
end
function Base.:*(a::Adjoint{T,<:AbstractArray{T,1}}, b::AxisVector{T}) where {T<:Number}
    p = *(a, parent(b))
    return _matmul(p, matmul_axes(axes(a), axes(b), axes(p)))
end
function Base.:*(a::Transpose{T,<:AbstractArray{T,1}}, b::AxisVector{T}) where {T<:Real}
    p = *(a, parent(b))
    return _matmul(p, matmul_axes(axes(a), axes(b), axes(p)))
end

#=
    get_factorization(F::Factorization, A::AbstractArray, d::Symbol)

Used internally to compose an `AxisArray` for each component of a factor
decomposition. `F` is the result of decomposition, `A` is an arry (likely
a subtype of `AxisArray`), and `d` is a symbol referring to a component
of the factorization.
=#
function get_factorization end

@doc """
    lu(A::AxisArray, args...; kwargs...)

Compute the LU factorization of an `AxisArray` `A`.

## Examples
```jldoctest
julia> using AxisIndices, LinearAlgebra

julia> m = AxisArray([1.0 2; 3 4], (2:3, 3:4));

julia> F = lu(m);

julia> axes(F.L)
(offset(1)(SimpleAxis(1:2)), SimpleAxis(1:2))

julia> axes(F.U)
(SimpleAxis(1:2), offset(2)(SimpleAxis(1:2)))

julia> F.p
2-element Array{Int64,1}:
 3
 2

julia> axes(F.P)
(offset(1)(SimpleAxis(1:2)), offset(1)(SimpleAxis(1:2)))

julia> axes(F.P * m)
(offset(1)(SimpleAxis(1:2)), offset(2)(SimpleAxis(1:2)))

julia> axes(F.L * F.U)
(offset(1)(SimpleAxis(1:2)), offset(2)(SimpleAxis(1:2)))

```
""" lu

function LinearAlgebra.lu!(A::AxisArray, args...; kwargs...)
    inner_lu = lu!(parent(A), args...; kwargs...)
    return LU(
        unsafe_reconstruct(A, getfield(inner_lu, :factors); axes=axes(A)),
        getfield(inner_lu, :ipiv),
        getfield(inner_lu, :info)
    )
end

function Base.parent(F::LU{T,<:AxisArray}) where {T}
    return LU(parent(getfield(F, :factors)), getfield(F, :ipiv), getfield(F, :info))
end

@inline function Base.getproperty(F::LU{T,<:AxisArray}, d::Symbol) where {T}
    return get_factorization(parent(F), getfield(F, :factors), d)
end

function get_factorization(F::LU, A::AbstractArray, d::Symbol)
    inner = getproperty(F, d)
    if d === :L
        return unsafe_reconstruct(A, inner, axes=(axes(A, 1), SimpleAxis(OneTo(size(inner, 2)))))
    elseif d === :U
        return unsafe_reconstruct(A, inner, axes=(SimpleAxis(OneTo(size(inner, 1))), axes(A, 2)))
    elseif d === :P
        return unsafe_reconstruct(A, inner, axes=(axes(A, 1), axes(A, 1)))
    elseif d === :p
        return _add_offset(axes(A, 1), inner)
        #=
        if is_offset_axis(axis)
            return inner .+ offsets(axis, 1)
        else
            return inner
        end
        =#
    else
        return inner
    end
end

"""
    lq(A::AxisArray, args...; kwargs...)

Compute the LQ factorization of an `AxisArray` `A`.

## Examples
```jldoctest
julia> using AxisIndices, LinearAlgebra

julia> m = AxisArray([1.0 2; 3 4], (2:3, 3:4));

julia> F = lq(m);

julia> axes(F.L)
(offset(1)(SimpleAxis(1:2)), SimpleAxis(1:2))

julia> axes(F.Q)
(SimpleAxis(1:2), offset(2)(SimpleAxis(1:2)))

julia> axes(F.L * F.Q)
(offset(1)(SimpleAxis(1:2)), offset(2)(SimpleAxis(1:2)))
```
"""
LinearAlgebra.lq(A::AxisArray, args...; kws...) = lq!(copy(A), args...; kws...)
function LinearAlgebra.lq!(A::AxisArray, args...; kwargs...)
    F = lq!(parent(A), args...; kwargs...)
    inner = getfield(F, :factors)
    return LQ(unsafe_reconstruct(A, inner, axes=axes(A)), getfield(F, :τ))
end
function Base.parent(F::LQ{T,<:AxisArray}) where {T}
    return LQ(parent(getfield(F, :factors)), getfield(F, :τ))
end

@inline function Base.getproperty(F::LQ{T,<:AxisArray}, d::Symbol) where {T}
    return get_factorization(parent(F), getfield(F, :factors), d)
end

function get_factorization(F::LQ, A::AbstractArray, d::Symbol)
    inner = getproperty(F, d)
    if d === :L
        return unsafe_reconstruct(A, inner, axes=(axes(A, 1), SimpleAxis(OneTo(size(inner, 2)))))
    elseif d === :Q
        return unsafe_reconstruct(A, inner, axes=(SimpleAxis(OneTo(size(inner, 1))), axes(A, 2)))
    else
        return inner
    end
end

const AIQRUnion{T} = Union{LinearAlgebra.QRCompactWY{T,<:AxisArray},
                                     QRPivoted{T,<:AxisArray},
                                     QR{T,<:AxisArray}}

"""
    qr(F::AxisArray, args...; kwargs...)

Compute the QR factorization of an `AxisArray` `A`.

## Examples
```jldoctest
julia> using AxisIndices, LinearAlgebra

julia> m = AxisArray([1.0 2; 3 4], (2:3, 3:4));

julia> F = qr(m, Val(true));

julia> axes(F.Q)
(offset(1)(SimpleAxis(1:2)), SimpleAxis(1:2))

julia> axes(F.R)
(SimpleAxis(1:2), offset(2)(SimpleAxis(1:2)))

julia> axes(F.Q * F.R)
(offset(1)(SimpleAxis(1:2)), offset(2)(SimpleAxis(1:2)))

julia> axes(F.p)
(Base.OneTo(2),)

julia> axes(F.P)
(offset(1)(SimpleAxis(1:2)), offset(1)(SimpleAxis(1:2)))

julia> axes(F.P * AxisArray([1.0 2; 3 4], (2:3, 3:4)))
(offset(1)(SimpleAxis(1:2)), offset(2)(SimpleAxis(1:2)))

```

"""
function LinearAlgebra.qr(A::AxisArray{T,2}, args...; kwargs...) where T
    Base.require_one_based_indexing(A)
    AA = similar(A, LinearAlgebra._qreltype(T), axes(A))
    copyto!(AA, A)
    return qr!(AA, args...; kwargs...)
end

function LinearAlgebra.qr!(a::AxisArray, args...; kwargs...)
    return _qr(a, qr!(parent(a), args...; kwargs...), axes(a))
end

function _qr(a::AxisArray, F::QR, axs::Tuple)
    return QR(unsafe_reconstruct(a, getfield(F, :factors); axes=axs), F.τ)
end
function Base.parent(F::QR{<:Any,<:AxisArray})
    return QR(parent(getfield(F, :factors)), getfield(F, :τ))
end

function _qr(a::AxisArray, F::LinearAlgebra.QRCompactWY, axs::Tuple)
    return LinearAlgebra.QRCompactWY(
        unsafe_reconstruct(a, getfield(F, :factors); axes=axs),
        F.T
    )
end
function Base.parent(F::LinearAlgebra.QRCompactWY{<:Any, <:AxisArray})
    return LinearAlgebra.QRCompactWY(parent(getfield(F, :factors)), getfield(F, :T))
end

function _qr(a::AxisArray, F::QRPivoted, axs::Tuple)
    return QRPivoted(
        unsafe_reconstruct(a, getfield(F, :factors), axes=axs),
        getfield(F, :τ),
        getfield(F, :jpvt)
    )
end
function Base.parent(F::QRPivoted{<:Any, <:AxisArray})
    return QRPivoted(parent(getfield(F, :factors)), getfield(F, :τ), getfield(F, :jpvt))
end

@inline function Base.getproperty(F::AIQRUnion, d::Symbol)
    return get_factorization(parent(F), getfield(F, :factors), d)
end

function get_factorization(F::Q, A::AbstractArray, d::Symbol) where {Q<:Union{LinearAlgebra.QRCompactWY,QRPivoted,QR}}
    inner = getproperty(F, d)
    if d === :Q
        return unsafe_reconstruct(A, inner; axes=(axes(A, 1), SimpleAxis(OneTo(size(inner, 2)))))
    elseif d === :R
        return unsafe_reconstruct(A, inner; axes=(SimpleAxis(OneTo(size(inner, 1))), axes(A, 2)))
    elseif F isa QRPivoted && d === :P
        return unsafe_reconstruct(A, inner; axes=(axes(A, 1), axes(A, 1)))
    elseif F isa QRPivoted && d === :p
        return _add_offset(axes(A, 1), inner)
        #=
        axis =
        if known_offset(axis) === nothing || known_offset(axis) !== 1
            return inner .+ offsets(axis, 1)
        else
            return inner
        end
        =#
    else
        return inner
    end
end

struct AxisSVD{T,F<:SVD{T},A<:AxisArray} <: Factorization{T}
    factor::F
    array::A
end

"""
    svd(F::AxisArray, args...; kwargs...)

Compute the singular value decomposition (SVD) of an `AxisArray` `A`.

## Examples
```jldoctest
julia> using AxisIndices, LinearAlgebra

julia> m = AxisArray([1.0 2; 3 4], (2:3, 3:4));

julia> F = svd(m);

julia> axes(F.U)
(offset(1)(SimpleAxis(1:2)), SimpleAxis(1:2))

julia> axes(F.V)
(offset(2)(SimpleAxis(1:2)), SimpleAxis(1:2))

julia> axes(F.Vt)
(SimpleAxis(1:2), offset(2)(SimpleAxis(1:2)))

julia> axes(F.U * Diagonal(F.S) * F.Vt)
(offset(1)(SimpleAxis(1:2)), offset(2)(SimpleAxis(1:2)))

```
"""
function LinearAlgebra.svd(A::AxisArray, args...; kwargs...)
    return AxisSVD(svd(parent(A), args...; kwargs...), A)
end

function LinearAlgebra.svd!(A::AxisArray, args...; kwargs...)
    return AxisSVD(svd!(parent(A), args...; kwargs...), A)
end

Base.parent(F::AxisSVD) = getfield(F, :factor)

Base.size(F::AxisSVD) = size(parent(F))

Base.size(F::AxisSVD, i) = size(parent(F), i)

function Base.propertynames(F::AxisSVD, private::Bool=false)
    return private ? (:V, fieldnames(typeof(parent(F)))...) : (:U, :S, :V, :Vt)
end

Base.show(io::IO, F::AxisSVD) = show(io, MIME"text/plain"(), F)
function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, F::AxisSVD{T}) where {T}
    print(io, "AxisSVD{$T}\n")
    println(io, "U factor:")
    show(io, mime, F.U)
    println(io, "\nsingular values:")
    show(io, mime, F.S)
    println(io, "\nVt factor:")
    show(io, mime, F.Vt)
end

LinearAlgebra.svdvals(A::AxisArray) = svdvals(parent(A))

# iteration for destructuring into components
Base.iterate(S::AxisSVD) = (S.U, Val(:S))
Base.iterate(S::AxisSVD, ::Val{:S}) = (S.S, Val(:V))
Base.iterate(S::AxisSVD, ::Val{:V}) = (S.V, Val(:done))
Base.iterate(S::AxisSVD, ::Val{:done}) = nothing
# TODO GeneralizedSVD

@inline function Base.getproperty(F::AxisSVD, d::Symbol)
    return get_factorization(parent(F), getfield(F, :array), d)
end

function get_factorization(F::SVD, A::AbstractArray, d::Symbol)
    inner = getproperty(F, d)
    if d === :U
        return unsafe_reconstruct(A, inner, axes=(axes(A, 1), SimpleAxis(OneTo(size(inner, 2)))))
    elseif d === :V
        return unsafe_reconstruct(A, inner, axes=(axes(A, 2), SimpleAxis(OneTo(size(inner, 2)))))
    elseif d === :Vt
        return unsafe_reconstruct(A, inner, axes=(SimpleAxis(OneTo(size(inner, 1))), axes(A, 2)))
    else  # d === :S
        return inner
    end
end

function LinearAlgebra.eigen(A::AxisArray{T,N,P,AI}; kwargs...) where {T,N,P,AI}
    vals, vecs = LinearAlgebra.eigen(parent(A); kwargs...)
    return Eigen(vals, unsafe_reconstruct(A, vecs; axes=axes(A)))
end

function LinearAlgebra.eigvals(A::AxisArray; kwargs...)
    return LinearAlgebra.eigvals(parent(A); kwargs...)
end

# TODO test this
function LinearAlgebra.eigen!(A::AxisArray{T,N,P,AI}; kwargs...) where {T,N,P,AI}
    vals, vecs = LinearAlgebra.eigen!(parent(A); kwargs...)
    return Eigen(vals, unsafe_reconstruct(A, vecs; axes=axes(A)))
end

function LinearAlgebra.eigvals!(A::AxisArray; kwargs...)
    return LinearAlgebra.eigvals!(parent(A); kwargs...)
end

#TODO eigen!(::AbstractArray, ::AbstractArray)


text:
__precompile__()

module Preprocessing

using Types, Master, SymPy

export createhighordermesh, preprocessing, initializeexasim, findexec

include("xiny.jl");
include("createdgnodes.jl");
include("getelemface.jl");
include("localbasis.jl");
include("sortcolumns.jl");
include("sortrows.jl");
include("uniquerows.jl");
include("mkf2e.jl");
include("mkf2t.jl");
include("mkv2t.jl");
include("node2elem.jl");
include("mkent2elem.jl");
include("mkelconcg.jl");
include("facenumbering.jl");
include("createhighordermesh.jl");
include("partition.jl");
include("neighboringelements.jl");
include("permindex.jl");
include("elementpartition2.jl");
include("faceconnectivity2.jl");
include("facepartition2.jl");
include("meshpartition2.jl");
# include("elementpartition.jl");
# include("faceconnectivity.jl");
# include("facepartition.jl");
# include("meshpartition.jl");
include("mkcgent2dgent.jl");
include("mkdge2dgf.jl");
include("mkelemblocks.jl");
include("mkfaceblocks.jl");
include("writesol.jl");
include("writeapp.jl");
include("writemaster.jl");
include("initializepde.jl");
include("initializemesh.jl");
include("initializeexasim.jl");
include("findexec.jl");

function preprocessing(app,mesh)

if !isdir("datain")
    mkdir("datain");
end
if !isdir("dataout")
    mkdir("dataout");
end

filename = "datain/";
fileapp = filename * "app.bin";
filemaster = filename * "master.bin";

if app.preprocessmode==0
    # update app structure
    app = writeapp(app,fileapp);
    return app;
end

app.nd  = size(mesh.p,1);
app.ncx = app.nd; #size(mesh.dgnodes,2);
nve,ne = size(mesh.t);

app.elemtype = 0;
if (app.nd==2) && (nve==4)
    app.elemtype=1;
end
if (app.nd==3) && (nve==8)
    app.elemtype=1;
end
app.pgauss = 2*app.porder;

master = Master.mkmaster(app.nd,app.porder,app.pgauss,app.elemtype,app.nodetype);
writemaster(master,filemaster);

app.boundaryconditions = mesh.boundarycondition;
app.uinf = app.externalparam;
nuinf = length(app.uinf);
nparam = length(app.physicsparam);
xdgsym = [SymPy.symbols("xdg$i") for i=1:app.ncx];
uinfsym = [SymPy.symbols("uinf$i") for i=1:nuinf];
paramsym = [SymPy.symbols("param$i") for i=1:nparam];
if isdefined(Main, Symbol("initu"))
    udgsym = Main.initu(xdgsym, paramsym, uinfsym);
    app.ncu = length(udgsym);
else
    error("pde.initu is not defined");
end
if isdefined(Main, Symbol("initv"))
    vdgsym = Main.initv(xdgsym, paramsym, uinfsym);
    app.nco = length(vdgsym);
else
    app.nco = 0;
end
if isdefined(Main, Symbol("initw"))
    wdgsym = Main.initw(xdgsym, paramsym, uinfsym);
    app.ncw = length(wdgsym);
else
    app.ncw = 0;
end

if app.model=="ModelC" || app.model=="modelC"
    app.wave = 0;
    app.nc = app.ncu;
elseif app.model=="ModelD" || app.model=="modelD"
    app.wave = 0;
    app.nc = round((app.ncu)*(app.nd+1));
elseif app.model=="ModelW" || app.model=="modelW"
    app.tdep = 1;
    app.wave = 1;
    app.nc = round((app.ncu)*(app.nd+1));
end
app.ncq = app.nc - app.ncu;
app.nch  = app.ncu;

if maximum(app.dt)>0.0
    app.tdep = 1;
else
    app.tdep = 0;
end

print("run facenumbering...\n");
mesh.f, mesh.tprd, t2t = facenumbering(mesh.p,mesh.t,app.elemtype,mesh.boundaryexpr,mesh.periodicexpr);

mpiprocs = app.mpiprocs;
dmd = Array{DMDStruct, 1}(undef, mpiprocs);
for i = 1:mpiprocs
    dmd[i] = DMDStruct();
end

# app.neb = 50000;
# app.nfb = 100000;

#dmd = meshpartition(dmd,mesh.p,mesh.t,mesh.f,t2t,mesh.tprd,app.elemtype,app.boundaryconditions,mesh.boundaryexpr,mesh.periodicexpr,app.porder,mpiprocs,app.metis);
dmd = meshpartition2(dmd,mesh.tprd,mesh.f,t2t,app.boundaryconditions,app.nd,app.elemtype,app.porder,mpiprocs,app.metis);
for i = 1:mpiprocs
    #xdg = mesh.dgnodes[:,:,dmd[i].elempart[:]];
    # create DG nodes
    xdg = createdgnodes(mesh.p,mesh.t[:,dmd[i].elempart[:]],mesh.f[:,dmd[i].elempart[:]],mesh.curvedboundary,mesh.curvedboundaryexpr,app.porder);

    cgelcon,rowent2elem,colent2elem,cgent2dgent,~ = mkcgent2dgent(xdg,1e-8);

    if mpiprocs==1
        writesol(filename  * "/sol",0,xdg);
        ne = length(dmd[i].elempart);
        eblks,nbe = mkelemblocks(ne,app.neb);
        eblks = vcat(eblks, zeros(Int, 1, size(eblks,2)));
        mf = cumsum([0; dmd[i].facepartpts[:]]);
        fblks,nbf = mkfaceblocks(mf,dmd[i].facepartbnd[:],app.nfb);
        neb = maximum(eblks[2,:]-eblks[1,:])+1;
        nfb = maximum(fblks[2,:]-fblks[1,:])+1;
    else
        writesol(filename  * "/sol",i,xdg);
        me = cumsum([0; dmd[i].elempartpts[1]; dmd[i].elempartpts[2]; dmd[i].elempartpts[3]]);
        eblks,nbe = mkfaceblocks(me,[0 1 2],app.neb);
        mf = cumsum([0; dmd[i].facepartpts[:]]);
        fblks,nbf = mkfaceblocks(mf,dmd[i].facepartbnd[:],app.nfb);
        neb = maximum(eblks[2,:]-eblks[1,:])+1;
        nfb = maximum(fblks[2,:]-fblks[1,:])+1;
    end

    npe = master.npe;
    nfe = size(master.perm,2);
    facecon1 = reshape(dmd[i].facecon[:,1,:],size(dmd[i].facecon,1), size(dmd[i].facecon,3));
    facecon2 = reshape(dmd[i].facecon[:,2,:],size(dmd[i].facecon,1), size(dmd[i].facecon,3));
    ind = [];
    for ii = 1:size(fblks,2)
        if fblks[3,ii]>0
            ind = vcat(ind, Vector(fblks[1,ii]:fblks[2,ii]));
        end
    end
    if length(ind)>0
        facecon2 = facecon2[:, setdiff(1:end, ind)];
    end
    rowe2f1,cole2f1,ent2ind1 = mkdge2dgf(facecon1,master.npe*length(dmd[i].elempart));
    rowe2f2,cole2f2,ent2ind2 = mkdge2dgf(facecon2,master.npe*length(dmd[i].elempart));
    ndims = zeros(Int,20,1);
    ndims[1] = size(mesh.p,1);
    ndims[2] = length(dmd[i].elempart);
    ndims[3] = sum(dmd[i].facepartpts);
    ndims[4] = maximum(mesh.t[:]);
    ndims[5] = nfe;
    ndims[6] = nbe;
    ndims[7] = neb;
    ndims[8] = nbf;
    ndims[9] = nfb;

    nsize = zeros(Int,21,1);
    nsize[1] = length(ndims);
    nsize[2] = length(dmd[i].facecon);
    nsize[3] = length(eblks);
    nsize[4] = length(fblks);
    nsize[5] = length(dmd[i].nbsd);
    nsize[6] = length(dmd[i].elemsend);
    nsize[7] = length(dmd[i].elemrecv);
    nsize[8] = length(dmd[i].elemsendpts);
    nsize[9] = length(dmd[i].elemrecvpts);
    nsize[10] = length(dmd[i].elempart);
    nsize[11] = length(dmd[i].elempartpts);
    nsize[12] = length(cgelcon);
    nsize[13] = length(rowent2elem);
    nsize[14] = length(cgent2dgent);
    nsize[15] = length(colent2elem);
    nsize[16] = length(rowe2f1);
    nsize[17] = length(cole2f1);
    nsize[18] = length(ent2ind1);
    nsize[19] = length(rowe2f2);
    nsize[20] = length(cole2f2);
    nsize[21] = length(ent2ind2);

    if (mpiprocs>1)
        print("Writing mesh into file " * string(i) * "\n");
        fileID = open(string(filename  * "/mesh", string(string(i), ".bin")),"w");
    else
        print("Writing mesh into file \n");
        fileID = open(string(filename  * "/mesh", ".bin"),"w");
    end

    write(fileID,Float64(length(nsize[:])));
    write(fileID,Float64.(nsize[:]));
    write(fileID,Float64.(ndims[:]));
    write(fileID,Float64.(reshape(permutedims(dmd[i].facecon,[2,1,3]),length(dmd[i].facecon),1)));
    write(fileID,Float64.(eblks[:]));
    write(fileID,Float64.(fblks[:]));
    write(fileID,Float64.(dmd[i].nbsd[:]));
    write(fileID,Float64.(dmd[i].elemsend[:]));
    write(fileID,Float64.(dmd[i].elemrecv[:]));
    write(fileID,Float64.(dmd[i].elemsendpts[:]));
    write(fileID,Float64.(dmd[i].elemrecvpts[:]));
    write(fileID,Float64.(dmd[i].elempart[:]));
    write(fileID,Float64.(dmd[i].elempartpts[:]));
    write(fileID,Float64.(cgelcon[:].-1));
    write(fileID,Float64.(rowent2elem[:]));
    write(fileID,Float64.(cgent2dgent[:].-1));
    write(fileID,Float64.(colent2elem[:].-1));
    write(fileID,Float64.(rowe2f1[:]));
    write(fileID,Float64.(cole2f1[:].-1));
    write(fileID,Float64.(ent2ind1[:].-1));
    write(fileID,Float64.(rowe2f2[:]));
    write(fileID,Float64.(cole2f2[:].-1));
    write(fileID,Float64.(ent2ind2[:].-1));
    close(fileID);
end

app = writeapp(app,fileapp);

mesh.telem = master.telem;
mesh.tface = master.telem;
mesh.xpe = master.xpe;
mesh.xpf = master.xpf;
for i = 1:mpiprocs
    dmd[i].nbsd = reshape([0],1,1);
    #dmd[i].elem2cpu = reshape([0],1,1);
    dmd[i].elemrecv = reshape([0],1,1);
    dmd[i].elemsend = reshape([0],1,1);
    dmd[i].elemrecvpts = reshape([0],1,1);
    dmd[i].elemsendpts = reshape([0],1,1);
    dmd[i].facepartpts = reshape([0],1,1);
    dmd[i].facepartbnd = reshape([0],1,1);
    dmd[i].facecon = reshape([0],1,1,1);
end

return app, mesh, master, dmd

end

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
using FibonacciOptimization, Test

@testset "fibonacci with accuracy" begin
    lower_bound, upper_bound = -1, 1
    func, accuracy = x -> x^2, 0.1
    true_x, true_f = 0, 0
    x, f = fibonacci_search(func, lower_bound, upper_bound, accuracy=accuracy)
    @test abs(x - true_x) < accuracy && abs(f - true_f) < accuracy

    lower_bound, upper_bound = 0, 1
    func, accuracy = x -> x^4 - 14x^3 + 60x^2 - 70x, 0.001
    true_x, true_f = 0.780884, -24.3696
    x, f = fibonacci_search(func, lower_bound, upper_bound, accuracy=accuracy)
    @test abs(x - true_x) < accuracy && abs(f - true_f) < accuracy

    lower_bound, upper_bound = pi, 2 * pi
    func, accuracy = x -> sin(x), 1e-4
    true_x, true_f = 1.5 * pi, -1
    x, f = fibonacci_search(func, lower_bound, upper_bound) # test default accuracy
    @test abs(x - true_x) < accuracy && abs(f - true_f) < accuracy
end

@testset "fibonacci with iteration number" begin
    lower_bound, upper_bound = -1, 1
    func, iter_num = x -> x^2, 5
    expected_x, expected_f = 0, 0
    @test fibonacci_search(func, lower_bound, upper_bound, iterations=iter_num) == (expected_x, expected_f)

    lower_bound, upper_bound = -2, 2
    func, iter_num = x -> x^4 - 14x^3 + 60x^2 - 70x, 12
    expected_x, expected_f = 0.781115879828326, -24.369599908878634
    @test fibonacci_search(func, lower_bound, upper_bound, iterations=iter_num) == (expected_x, expected_f)

    lower_bound, upper_bound = -5, -3
    func, iter_num = x -> sin(x) * cos(x), 25
    expected_x, expected_f = -3.926989200365754, -0.49999999999738654
    @test fibonacci_search(func, lower_bound, upper_bound, iterations=iter_num) == (expected_x, expected_f)

end