text:
"""
This module will be deleted shortly, when we beginning using a rasterized PSF.
In the future, PSF code will appear only in `SDSSIO.jl` (to load the eigen-PSF)
and `model/psf.jl` (functions for computing parts of the log probability
relating to the PSF).
"""
module PSF

using Celeste
using ..Model
using ..Transform
using ..SensitiveFloats.SensitiveFloat
import ..SensitiveFloats.clear!
import ..SensitiveFloats

import Optim
import WCS

include("bivariate_normals.jl")

const ID_MAT_2D = eye(Float64, 2)

"""
A data type to store functions related to optimizing a PSF fit. Initialize
using the transform and number of components, and then call fit_psf
to perform a fit to a specificed psf and initial parameters.
"""
type PsfOptimizer
    psf_transform::DataTransform
    ftol::Float64
    grtol::Float64
    num_iters::Int
    raw_psf::Matrix{Float64}
    K::Int

    # Variable that will be allocated in optimization:
    x_mat::Matrix{Float64}
    bvn_derivs::BivariateNormalDerivatives{Float64}

    log_pdf::SensitiveFloat{PsfParams, Float64}
    pdf::SensitiveFloat{PsfParams, Float64}
    pixel_value::SensitiveFloat{PsfParams, Float64}
    squared_error::SensitiveFloat{PsfParams, Float64}
    sf_free::SensitiveFloat{PsfParams, Float64}

    psf_params_free_vec_cache::Vector{Float64}

    # functions
    psf_2df::Optim.TwiceDifferentiableFunction
    fit_psf::Function

    function PsfOptimizer(psf_transform::DataTransform, K::Int;
                          verbose::Bool=false, ftol::Float64=1e-9, grtol::Float64=1e-9)

        num_iters = 50

        bvn_derivs = BivariateNormalDerivatives{Float64}(Float64)

        log_pdf = SensitiveFloats.zero_sensitive_float(PsfParams, Float64, 1)
        pdf = SensitiveFloats.zero_sensitive_float(PsfParams, Float64, 1)
        pixel_value = SensitiveFloats.zero_sensitive_float(PsfParams, Float64, K)
        squared_error = SensitiveFloats.zero_sensitive_float(PsfParams, Float64, K)
        sf_free = SensitiveFloats.zero_sensitive_float(PsfParams, Float64, K)

        x_mat = Array(Float64, 0, 0)
        raw_psf = Array(Float64, 0, 0)
        psf_params_free_vec_cache = fill(NaN, K * length(PsfParams))

        function psf_fit_for_optim{NumType <: Number}(
                psf_params_free_vec::Vector{NumType})

            if psf_params_free_vec == psf_params_free_vec_cache
                return sf_free
            else
                psf_params_free_vec_cache = deepcopy(psf_params_free_vec)
            end
            psf_params_free = unwrap_psf_params(psf_params_free_vec)
            psf_params = constrain_psf_params(psf_params_free, psf_transform)

            # Update squared_error in place.
            evaluate_psf_fit!(
                    psf_params, raw_psf, x_mat, bvn_derivs,
                    log_pdf, pdf, pixel_value, squared_error, true)

            # Update sf_free in place.
            transform_psf_sensitive_float!(
                psf_params, psf_transform, squared_error, sf_free, true)

            sf_free
        end

        function psf_fit_value{NumType <: Number}(psf_params_free_vec::Vector{NumType})
            psf_fit_for_optim(psf_params_free_vec).v[1]
        end

        function psf_fit_grad!(
                psf_params_free_vec::Vector{Float64}, grad::Vector{Float64})
            grad[:] = psf_fit_for_optim(psf_params_free_vec).d[:]
        end

        function psf_fit_hess!(
                psf_params_free_vec::Vector{Float64}, hess::Matrix{Float64})
            hess[:] = psf_fit_for_optim(psf_params_free_vec).h
            hess[:] = 0.5 * (hess + hess')
        end

        psf_2df = Optim.TwiceDifferentiableFunction(
            psf_fit_value, psf_fit_grad!, psf_fit_hess!)

        function fit_psf(psf::Matrix{Float64}, initial_params::Vector{Vector{Float64}})
            raw_psf = psf
            x_mat = get_x_matrix_from_psf(raw_psf)
            psf_params_free = unconstrain_psf_params(initial_params, psf_transform)
            psf_params_free_vec = vec(wrap_psf_params(psf_params_free))

            tr_method =
                    Optim.NewtonTrustRegion(initial_delta=10.0, delta_hat=1e9, eta=0.1,
                                            rho_lower=0.2, rho_upper=0.75)

            options = Optim.OptimizationOptions(;
                    x_tol = 0.0, # Don't allow convergence in params
                    f_tol = ftol, g_tol = grtol,
                    iterations = num_iters, store_trace = verbose,
                    show_trace = false, extended_trace = verbose)

            nm_result = Optim.optimize(psf_fit_value,
                                       psf_fit_grad!,
                                       psf_fit_hess!,
                                       psf_params_free_vec, tr_method, options)
            nm_result
        end

        new(psf_transform, ftol, grtol, num_iters, raw_psf, K,
            x_mat, bvn_derivs, log_pdf, pdf, pixel_value, squared_error, sf_free,
            psf_params_free_vec_cache, psf_2df, fit_psf)
    end
end

"""
Return an image of a Celeste GMM PSF evaluated at rows, cols.

Args:
- psf_array: The PSF to be evaluated as an array of PsfComponent
- rows: The rows in the image (usually in pixel coordinates)
- cols: The column in the image (usually in pixel coordinates)

Returns:
 - The PSF values at rows and cols.    The default size is the same as
     that returned by get_psf_at_point applied to FITS header values.

Note that the point in the image at which the PSF is evaluated --
that is, the center of the image returned by this function -- is
already implicit in the value of psf_array.
"""
function get_psf_at_point(psf_array::Array{PsfComponent, 1};
                          rows=collect(-25:25), cols=collect(-25:25))

    function get_psf_value(psf::PsfComponent, row::Float64, col::Float64)
         x = Float64[row, col] - psf.xiBar
         exp_term = exp(-0.5 * x' * psf.tauBarInv * x - 0.5 * psf.tauBarLd)
         (psf.alphaBar * exp_term / (2 * pi))[1]
    end

    Float64[
     sum([ get_psf_value(psf, float(row), float(col)) for psf in psf_array ])
         for row in rows, col in cols ]
end

"""
Get the PSF located at a particular world location in an image.

Args:
 - world_loc: A location in world coordinates.
 - img: An TiledImage

Returns:
 - An array of PsfComponent objects that represents the PSF as a mixture
     of Gaussians.
"""
function get_source_psf(world_loc::Vector{Float64}, img::TiledImage)
    # Some stamps or simulated data have no raw psf information.    In that case,
    # just use the psf from the image.
    if size(img.raw_psf_comp.rrows) == (0, 0)
        # Also return a vector of empty psf params
        return img.psf, fill(fill(NaN, length(PsfParams)), psf_K)
    else
        pixel_loc = WCS.world_to_pix(img.wcs, world_loc)
        psfstamp = img.raw_psf_comp(pixel_loc[1], pixel_loc[2])
        return PSF.fit_raw_psf_for_celeste(psfstamp)
    end
end

"""
Get the PSF located at a particular world location in an image.

Args:
 - world_loc: A location in world coordinates.
 - img: An TiledImage

Returns:
 - An array of PsfComponent objects that represents the PSF as a mixture
     of Gaussians.
"""
function get_source_psf(world_loc::Vector{Float64},
                        img::TiledImage,
                        psf_optimizer::PSF.PsfOptimizer,
                        initial_psf_params::Vector{Vector{Float64}})
    # Some stamps or simulated data have no raw psf information. In that case,
    # just use the psf from the image.
    if size(img.raw_psf_comp.rrows) == (0, 0)
        return img.psf
    else
        pixel_loc = WCS.world_to_pix(img.wcs, world_loc)
        psfstamp = img.raw_psf_comp(pixel_loc[1], pixel_loc[2])
        return PSF.fit_raw_psf_for_celeste(psfstamp, psf_optimizer, initial_psf_params)
    end
end

"""
Get initial values for a set of PSF parameters.

Args:
    - K: The number of components.
    - for_test: If true, return a set of slightly asymmetric parameters so that
                            derivatives are interesting for testing.

Returns:
    - A vector of PSF parameters, one for each component.
"""
function initialize_psf_params(K::Int; for_test::Bool=false)
    psf_params = Array(Vector{Float64}, K)
    for k=1:K
        if for_test
            # Choose asymmetric values for testing.
            psf_params[k] = zeros(length(PsfParams))
            psf_params[k][psf_ids.mu] = [0.1, 0.2]
            psf_params[k][psf_ids.e_axis] = 0.8
            psf_params[k][psf_ids.e_angle] = pi / 4
            psf_params[k][psf_ids.e_scale] = sqrt(2 * k)
            psf_params[k][psf_ids.weight] = 1 / K + k / 10
        else
            psf_params[k] = zeros(length(PsfParams))
            psf_params[k][psf_ids.mu] = [0.0, 0.0]
            psf_params[k][psf_ids.e_axis] = 0.95
            psf_params[k][psf_ids.e_angle] = 0.0
            psf_params[k][psf_ids.e_scale] = sqrt(2 * k)
            psf_params[k][psf_ids.weight] = 1 / K
        end
    end

    psf_params
end

"""
Given PSF parameters, return a DataTransform to transfrom them to and from
an unconstrained parameterization.

Args:
    - psf_params: The parameters to be transformed
    - scale: A vector as long as PsfParams to define a linear rescaling.

Returns:
    - A DataTransform object.
"""
function get_psf_transform(
        psf_params::Vector{Vector{Float64}};
        scale::Vector{Float64}=ones(length(PsfParams)))

    K = length(psf_params)
    bounds = Array(ParamBounds, length(psf_params))
    # Note that, for numerical reasons, the bounds must be on the scale
    # of reasonably meaningful changes.
    for k in 1:K
        bounds[k] = ParamBounds()
        bounds[k][:mu] = ParamBox[ ParamBox(-5.0, 5.0, scale[psf_ids.mu[1]]),
                                                             ParamBox(-5.0, 5.0, scale[psf_ids.mu[2]]) ]
        bounds[k][:e_axis] = ParamBox[ ParamBox(0.1, 1.0, scale[psf_ids.e_axis] ) ]
        bounds[k][:e_angle] =
            ParamBox[ ParamBox(-4 * pi, 4 * pi, scale[psf_ids.e_angle] ) ]
        bounds[k][:e_scale] =
            ParamBox[ ParamBox(0.05, 10.0, scale[psf_ids.e_scale] ) ]

        # Note that the weights do not need to sum to one.
        bounds[k][:weight] = ParamBox[ ParamBox(0.05, 2.0, scale[psf_ids.weight] ) ]
    end
    DataTransform(bounds, active_sources=collect(1:K), S=K)
end

"""
Transform psf_params to and from psf_params_free in place.

Args:
    - psf_params, psf_params_free: The constrained and unconstrained parameters,
            respectively.
    - psf_transform: The DataTransform to be applied.
    - to_unconstrained: If true, then update psf_params_free with the unconstrained
            representation of psf_params.    If false, update psf_params with the
            constrained version of psf_params_free.

Returns:
    - Updates psf_params or psf_params_free in place.
"""
function transform_psf_params!{NumType <: Number}(
        psf_params::Vector{Vector{NumType}}, psf_params_free::Vector{Vector{NumType}},
        psf_transform::DataTransform, to_unconstrained::Bool)

    for k=1:length(psf_params)
        for (param, constraint_vec) in psf_transform.bounds[k]
            for ind in 1:length(getfield(psf_ids, param))
                param_ind = getfield(psf_ids, param)[ind]
                constraint = constraint_vec[ind]
                to_unconstrained ?
                    psf_params_free[k][param_ind] =
                        Transform.unbox_parameter(psf_params[k][param_ind], constraint):
                    psf_params[k][param_ind] =
                        Transform.box_parameter(psf_params_free[k][param_ind], constraint)
            end
        end
    end

    true # return type
end

"""
Allocate memory for and return a constrained parameter set.
"""
function constrain_psf_params{NumType <: Number}(
        psf_params_free::Vector{Vector{NumType}}, psf_transform::DataTransform)

    K = length(psf_params_free)
    psf_params = Array(Vector{NumType}, K)
    for k=1:K
        psf_params[k] = zeros(NumType, length(PsfParams))
    end

    transform_psf_params!(psf_params, psf_params_free, psf_transform, false)

    psf_params
end

"""
Allocate memory for and return an unconstrained parameter set.
"""
function unconstrain_psf_params{NumType <: Number}(
        psf_params::Vector{Vector{NumType}}, psf_transform::DataTransform)

    K = length(psf_params)
    psf_params_free = Array(Vector{NumType}, K)
    for k=1:K
        psf_params_free[k] = zeros(NumType, length(PsfParams))
    end

    transform_psf_params!(psf_params, psf_params_free, psf_transform, true)

    psf_params_free
end

"""
Convert a single vector of psf parameters to a vector of vectors.
"""
function unwrap_psf_params{NumType <: Number}(psf_param_vec::Vector{NumType})
    @assert length(psf_param_vec) % length(PsfParams) == 0
    K = round(Int, length(psf_param_vec) / length(PsfParams))
    psf_param_mat = reshape(psf_param_vec, length(PsfParams), K)
    psf_params = Array(Vector{NumType}, K)
    for k = 1:K
        psf_params[k] = psf_param_mat[:, k]
    end
    psf_params
end

"""
Convert a vector of vectors of psf parameters to a single vector.
"""
function wrap_psf_params{NumType <: Number}(psf_params::Vector{Vector{NumType}})
    psf_params_mat = zeros(NumType, length(PsfParams), length(psf_params))
    for k=1:length(psf_params)
        psf_params_mat[:, k] = psf_params[k]
    end
    psf_params_mat[:]
end

"""
Return a sensitive float representing the value of the psf at pixel x
with all its associated derivatives (with respect to the constrained
parameterization).    Note that "fit" is a bad name -- it is just the pixel value.
TODO: fix the name.

Args:
    - x: The 2d location at which to evaluate the pdf
    - calculate_derivs: If false, only update the value of pixel_value
    - other values: Pre-allocated memory for intermediate calculations

Returns:
    - Updates pixel_value in place (and all the other placeholder values as well)
"""
function evaluate_psf_pixel_fit!{NumType <: Number}(
        x::Vector{Float64}, psf_params::Vector{Vector{NumType}},
        sigma_vec::Vector{Matrix{NumType}},
        sig_sf_vec::Vector{GalaxySigmaDerivs{NumType}},
        bvn_vec::Vector{BvnComponent{NumType}},
        bvn_derivs::BivariateNormalDerivatives{NumType},
        log_pdf::SensitiveFloat{PsfParams, NumType},
        pdf::SensitiveFloat{PsfParams, NumType},
        pixel_value::SensitiveFloat{PsfParams, NumType},
        calculate_derivs::Bool)

    clear!(pixel_value)

    K = length(psf_params)
    sigma_ids = (psf_ids.e_axis, psf_ids.e_angle, psf_ids.e_scale)
    for k = 1:K
        # I will put in the weights later so that the log pdf sensitive float
        # is accurate.
        bvn = bvn_vec[k]
        eval_bvn_pdf!(bvn_derivs, bvn, x)
        get_bvn_derivs!(bvn_derivs, bvn, true, true)
        transform_bvn_derivs!(bvn_derivs, sig_sf_vec[k], ID_MAT_2D, true)

        clear!(log_pdf)
        clear!(pdf)

        # This is redundant, but it's what eval_bvn_pdf returns.
        log_pdf.v[1] = log(bvn_derivs.f_pre[1])

        if calculate_derivs
            for ind=1:2
                log_pdf.d[psf_ids.mu[ind]] = bvn_derivs.bvn_u_d[ind]
            end
            for ind=1:3
                log_pdf.d[sigma_ids[ind]] = bvn_derivs.bvn_s_d[ind]
            end
            log_pdf.d[psf_ids.weight] = 0

            for ind1 = 1:2, ind2 = 1:2
                log_pdf.h[psf_ids.mu[ind1], psf_ids.mu[ind2]] =
                    bvn_derivs.bvn_uu_h[ind1, ind2]
            end
            for mu_ind = 1:2, sig_ind = 1:3
                log_pdf.h[psf_ids.mu[mu_ind], sigma_ids[sig_ind]] =
                log_pdf.h[sigma_ids[sig_ind], psf_ids.mu[mu_ind]] =
                    bvn_derivs.bvn_us_h[mu_ind, sig_ind]
            end
            for ind1 = 1:3, ind2 = 1:3
                log_pdf.h[sigma_ids[ind1], sigma_ids[ind2]] =
                    bvn_derivs.bvn_ss_h[ind1, ind2]
            end
        end

        pdf_val = exp(log_pdf.v[1])
        pdf.v[1] = pdf_val

        if calculate_derivs
            for ind1 = 1:length(PsfParams)
                if ind1 == psf_ids.weight
                    pdf.d[ind1] = pdf_val
                else
                    pdf.d[ind1] = psf_params[k][psf_ids.weight] * pdf_val * log_pdf.d[ind1]
                end

                for ind2 = 1:ind1
                    pdf.h[ind1, ind2] = pdf.h[ind2, ind1] =
                        psf_params[k][psf_ids.weight] * pdf_val *
                        (log_pdf.h[ind1, ind2] + log_pdf.d[ind1] * log_pdf.d[ind2])
                end
            end

            # Weight hessian terms.
            for ind1 = 1:length(PsfParams)
                pdf.h[psf_ids.weight, ind1] = pdf.h[ind1, psf_ids.weight] =
                    pdf_val * log_pdf.d[ind1]
            end

        end

        pdf.v[1] *= psf_params[k][psf_ids.weight]

        SensitiveFloats.add_sources_sf!(pixel_value, pdf, k, calculate_derivs)
    end

    true # Set return type
end

"""
Convert PSF parameters to covariance matrices and derivatives and BvnComponents.
"""
function get_sigma_from_params{NumType <: Number}(
        psf_params::Vector{Vector{NumType}})

    K = length(psf_params)
    sigma_vec = Array(Matrix{NumType}, K)
    sig_sf_vec = Array(GalaxySigmaDerivs{NumType}, K)
    bvn_vec = Array(BvnComponent{NumType}, K)
    for k = 1:K
        sigma_vec[k] = get_bvn_cov(psf_params[k][psf_ids.e_axis],
                                                                        psf_params[k][psf_ids.e_angle],
                                                                        psf_params[k][psf_ids.e_scale])
        sig_sf_vec[k] = GalaxySigmaDerivs(
            psf_params[k][psf_ids.e_angle],
            psf_params[k][psf_ids.e_axis],
            psf_params[k][psf_ids.e_scale], sigma_vec[k], calculate_tensor=true)

        bvn_vec[k] =
            BvnComponent{NumType}(psf_params[k][psf_ids.mu], sigma_vec[k], 1.0)
    end
    sigma_vec, sig_sf_vec, bvn_vec
end

"""
evaluate_psf_fit but with pre-allocated memory for intermediate calculations.
"""
function evaluate_psf_fit!{NumType <: Number}(
        psf_params::Vector{Vector{NumType}}, raw_psf::Matrix{Float64},
        x_mat::Matrix{Vector{Float64}},
        bvn_derivs::BivariateNormalDerivatives{NumType},
        log_pdf::SensitiveFloat{PsfParams, NumType},
        pdf::SensitiveFloat{PsfParams, NumType},
        pixel_value::SensitiveFloat{PsfParams, NumType},
        squared_error::SensitiveFloat{PsfParams, NumType},
        calculate_derivs::Bool)

    K = length(psf_params)
    sigma_vec, sig_sf_vec, bvn_vec = get_sigma_from_params(psf_params)
    clear!(squared_error)

    for x_ind in 1:length(x_mat)
        clear!(pixel_value)
        evaluate_psf_pixel_fit!(
                x_mat[x_ind], psf_params, sigma_vec, sig_sf_vec, bvn_vec,
                bvn_derivs, log_pdf, pdf, pixel_value, calculate_derivs)

        diff = (pixel_value.v[1] - raw_psf[x_ind])
        squared_error.v[1] +=    diff ^ 2
        if calculate_derivs
            for ind1 = 1:length(squared_error.d)
                squared_error.d[ind1] += 2 * diff * pixel_value.d[ind1]
                for ind2 = 1:ind1
                    squared_error.h[ind1, ind2] +=
                        2 * (diff * pixel_value.h[ind1, ind2] +
                                 pixel_value.d[ind1] * pixel_value.d[ind2]')
                    squared_error.h[ind2, ind1] = squared_error.h[ind1, ind2]
                end
            end
        end # if calculate_derivs
    end

    squared_error
end

"""
Transform a sf, a SensitiveFloat that has derivatives with respect to
constrianed parameters, into sf_free, a SensitiveFloat that has derivatives
with respect to unconstrained parameteres.

Args:
    - psf_params: The constrianed parameters
    - psf_transform: The data transform to be applied
    - sf: The SensitiveFloat with derivatives with respect to the constrained params
    - sf_free: Updated in place.    The SensitiveFloat with derivatives with respect
                         to the unconstrained parameters.
    - calculate_derivs: If false, only calculate the value.

Returns:
    - Updates sf_free in place.
"""
function transform_psf_sensitive_float!{NumType <: Number}(
        psf_params::Vector{Vector{NumType}}, psf_transform::Transform.DataTransform,
        sf::SensitiveFloat{PsfParams, NumType}, sf_free::SensitiveFloat{PsfParams, NumType},
        calculate_derivs::Bool)

    sf_free.v[1] = sf.v[1]
    if calculate_derivs
        K = length(psf_params)

        # These are the hessians of each individual parameter's transform.    We
        # can represent it this way since each parameter's transform only depends on
        # its own value and not on others.

        # This is the diagonal of the Jacobian transform.
        jacobian_diag = zeros(length(PsfParams) * K)

        # These are hte Hessians of the each parameter's transform.
        hessian_values = zeros(length(PsfParams) * K)

        for k = 1:K
            offset = length(PsfParams) * (k - 1)
            for ind = 1:2
                mu_ind = psf_ids.mu[ind]
                jac, hess = Transform.box_derivatives(
                    psf_params[k][mu_ind], psf_transform.bounds[k][:mu][ind])
                jacobian_diag[offset + mu_ind] = jac
                hessian_values[offset + mu_ind] = hess
            end

            # The rest are one-dimensional.
            for field in setdiff(fieldnames(PsfParams), [ :mu ])
                ind = getfield(psf_ids, field)
                jac, hess = Transform.box_derivatives(
                    psf_params[k][ind], psf_transform.bounds[k][field][1])
                jacobian_diag[offset + ind] = jac
                hessian_values[offset + ind] = hess
            end
        end

        # Apply the transformations.
        for ind1 = 1:length(sf_free.d)
            sf_free.d[ind1] = jacobian_diag[ind1] * sf.d[ind1]

            for ind2 = 1:ind1
                # Calculate the Hessian
                sf_free.h[ind1, ind2] = sf_free.h[ind2, ind1] =
                    (jacobian_diag[ind1] * jacobian_diag[ind2]) * sf.h[ind1, ind2]
                if ind1 == ind2
                    sf_free.h[ind1, ind2] +=    hessian_values[ind1] * sf.d[ind1]
                    diagm(hessian_values .* sf.d[:])
                end
            end
        end
    end

    true # return type
end

"""
Given a psf matrix, return a matrix of the same size, where each element of
the matrix is a 2-length vector of the [x1, x2] location of that matrix cell.
The locations are chosen so that the scale is pixels, but the center of the
psf is [0, 0].
"""
function get_x_matrix_from_psf(psf::Matrix{Float64})
    psf_center = Float64[ (size(psf, i) - 1) / 2 + 1 for i=1:2 ]
    Vector{Float64}[ Float64[i, j] - psf_center for i=1:size(psf, 1), j=1:size(psf, 2) ]
end

"""
Fit a mixture of 2d Gaussians to a PSF image (evaluated at a single point).

Args:
    - raw_psf: An matrix image of the point spread function.
    - psf_optimizer: A PsfOptimizer object for the fit.

Returns:
    - A vector of PsfComponents fit to the raw_psf.
"""
function fit_raw_psf_for_celeste(
        raw_psf::Array{Float64, 2}, psf_optimizer::PsfOptimizer,
        initial_psf_params::Vector{Vector{Float64}})

    K = length(initial_psf_params)
    @assert K == psf_optimizer.K
    optim_result = psf_optimizer.fit_psf(raw_psf, initial_psf_params)
    psf_params_fit =
        constrain_psf_params(
            unwrap_psf_params(optim_result.minimum), psf_optimizer.psf_transform)

    sigma_vec = get_sigma_from_params(psf_params_fit)[1]
    celeste_psf = Array(PsfComponent, K)
    for k=1:K
        mu = psf_params_fit[k][psf_ids.mu]
        weight = psf_params_fit[k][psf_ids.weight]
        celeste_psf[k] = PsfComponent(weight, mu, sigma_vec[k])
    end

    celeste_psf, psf_params_fit
end

"""
Fit a mixture of 2d Gaussians to a PSF image (evaluated at a single point).

Args:
 - raw_psf: An matrix image of the point spread function.
 - K: The number of components to fit.

 Returns:
     - A vector of PsfComponents fit to the raw_psf.
"""
function fit_raw_psf_for_celeste(raw_psf::Array{Float64, 2}; K=psf_K, ftol=1e-9)
    psf_params = PSF.initialize_psf_params(K, for_test=false)
    psf_transform = PSF.get_psf_transform(psf_params)
    psf_optimizer = PsfOptimizer(psf_transform, K)
    psf_optimizer.ftol = ftol
    fit_raw_psf_for_celeste(raw_psf, psf_optimizer, psf_params)
end

end # module


text:
# example "Case f" (plot in ipad)
# this is the bad diamond case
# Claudia August 2014
#
# in julia: include("case_f_example.jl")

ed1=Edge(1,0.6,true,0.7);
ed2=Edge(2,0.7,true,0.3);
ed3=Edge(3,0.9);
ed4=Edge(4,0.1);
ed5=Edge(5,0.2);
ed6=Edge(6,0.1);
ed7=Edge(7,0.1);
ed8=Edge(8,0.1);
ed9=Edge(9,0.1);
ed10=Edge(10,0.1);

n1=Node(1,false,false,[ed1,ed5,ed6]);
n2=Node(2,false,true,[ed1,ed2,ed3]);
n3=Node(3,false,false,[ed2,ed4,ed7]);
n4=Node(4,true,false,[ed3]);
n5=Node(5,false,false,[ed4,ed5,ed9]);
n6=Node(6,true, false,[ed6]);
n7=Node(7,true, false,[ed7]);
n8=Node(8,true, false,[ed8]);
n9=Node(9,false,false,[ed8,ed9,ed10]);
n10=Node(10,true, false,[ed10]);

setNode!(ed1,n1)
setNode!(ed1,n2)
setNode!(ed2,n2);
setNode!(ed2,n3);
setNode!(ed5,[n5,n1]);
setNode!(ed4,[n3,n5]);
setNode!(ed3,n2);
setNode!(ed3,n4);
setNode!(ed6,[n1,n6]);
setNode!(ed7,[n3,n7]);
setNode!(ed8,[n8,n9]);
setNode!(ed9,[n5,n9]);
setNode!(ed10,[n9,n10]);

#ed1.inCycle=2;
#ed2.inCycle=2;
#ed4.inCycle=2;
#ed5.inCycle=2;

net=HybridNetwork([n1,n2,n3,n4,n5,n6,n7,n8,n9,n10],[ed1,ed2,ed3,ed4,ed5,ed6,ed7,ed8,ed9,ed10]);
node=searchHybridNode(net);
n2.name = "H1"; n4.name = "4"; n6.name = "6"; n7.name = "7";
n8.name = "8"; n10.name = "10"
net.names=["1","2","3","4","5","6","7","8","9","10"]

flag, nocycle,edges, nodes = updateInCycle!(net,node[1]);
flag2, edges2 = updateContainRoot!(net,node[1]);
flag3, edges3 = updateGammaz!(net,node[1]);

#printEdges(net)

#deleteHybrid!(node[1],net,false)

#changeDirectionUpdate!(node[1],net);

ext: jl


text:
using FractionalDiffEq, Plots
plotlyjs()
rightfun(x, t) = -(1+x).*exp.(-t).*x.^3
T=1
M=100
N=M
h=1/M
tau=T/N
x=collect(0:h:1)
t=collect(0:tau:T)
alph=1.8
initial_condition(x) = x.^3
left_boundry(t) = 0
right_boundry(t) = exp.(-t)
d(x)  = gamma(2.2)*x.^2.8/6
(X, Y)=meshgrid(x, t)
sol=solve(alph, d, rightfun, M, N, initial_condition, left_boundry, right_boundry, GLDiff())
plot(X, Y, sol, st=:surface)


text:
module SixelTerm

using ImageMagick
using FileIO

struct SixelDisplay <: AbstractDisplay end

function Base.display(d::SixelDisplay, m::MIME{Symbol("image/png")}, x)
    a = repr("image/png", x)
    im = ImageMagick.load_(a)
    io = IOBuffer()
    s = Stream(format"six", io, "t.six")
    ImageMagick.save(s, im)
    write(stdout, take!(io))
end

function Base.display(d::SixelDisplay, x)
    if showable("image/png", x)
        display(d, "image/png", x)
    else
        # fall through to the Displays lower in the display stack
        throw(MethodError(display, "nope"))
    end
end

function __init__()
    Base.Multimedia.pushdisplay(SixelDisplay())
end

end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm, SparseArrays, DataStructures
using Distributed

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

function splitedges(V,copEV,sigma,return_edge_map,multiproc)
@show multiproc

    edgenum = size(copEV, 1)
    edge_map = Array{Array{Int, 1}, 1}(undef,edgenum)
    rV = Lar.Points(zeros(0, 2))
    rEV = SparseArrays.spzeros(Int8, 0, 0)
    finalcells_num = 0

    if (multiproc == true)
        in_chan = Distributed.RemoteChannel(()->Channel{Int64}(0))
        out_chan = Distributed.RemoteChannel(()->Channel{Tuple}(0))
        ordered_dict = SortedDict{Int64,Tuple}()
        @async begin
            for i in 1:edgenum
                put!(in_chan,i)
            end
            for p in Distributed.workers()
                put!(in_chan,-1)
            end
        end
        for p in Distributed.workers()
            @async Base.remote_do(frag_edge_channel, p, in_chan, out_chan, V, copEV)
        end
        for i in 1:edgenum
            frag_done_job = take!(out_chan)
            ordered_dict[frag_done_job[1]] = frag_done_job[2]
        end
        for (dkey, dval) in ordered_dict
            i = dkey
            v, ev = dval
            newedges_nums = map(x->x+finalcells_num, collect(1:size(ev, 1)))

            edge_map[i] = newedges_nums

            finalcells_num += size(ev, 1)

            rV, rEV = Lar.skel_merge(rV, rEV, v, ev)
        end
    else
        for i in 1:edgenum
            v, ev = Lar.Arrangement.frag_edge(V, copEV, i)
            newedges_nums = map(x->x+finalcells_num, collect(1:size(ev, 1)))
            edge_map[i] = newedges_nums
            finalcells_num += size(ev, 1)
            rV, rEV = Lar.skel_merge(rV, rEV, v, ev)
        end
    end
    return V,copEV,sigma,return_edge_map,edge_map,rV, rEV
end

function deleteoutsidesigma(V,copEV,sigma,return_edge_map,edge_map)
    # Deletes edges outside sigma area
    if sigma.n > 0
        todel = []

        new_edges = []
        map(i->new_edges=union(new_edges, edge_map[i]), sigma.nzind)
        ev = copEV[new_edges, :]

        for e in 1:copEV.m
            if !(e in new_edges)

                vidxs = copEV[e, :].nzind
                v1, v2 = map(i->V[vidxs[i], :], [1,2])
                centroid = .5*(v1 + v2)

                if ! Lar.point_in_face(centroid, V, ev)
                    push!(todel, e)
                end
            end
        end

        for i in reverse(todel)
            for row in edge_map

                filter!(x->x!=i, row)

                for j in 1:length(row)
                    if row[j] > i
                        row[j] -= 1
                    end
                end
            end
        end

        V, copEV = Lar.delete_edges(todel, V, copEV)
    end
    return V,copEV,sigma,return_edge_map,edge_map
end

function minimal_cycles(angles_fn::Function, verbose=false)

    function _minimal_cycles(V::Lar.Points,
    ld_bounds::Lar.ChainOp)

        lld_cellsnum, ld_cellsnum = size(ld_bounds)
        count_marks = zeros(Int8, ld_cellsnum)
        dir_marks = zeros(Int8, ld_cellsnum)
        d_bounds = spzeros(Int8, ld_cellsnum, 0)

        angles = Array{Array{Int64, 1}, 1}(undef, lld_cellsnum)

        function get_seed_cell()
            s = -1
            for i in 1:ld_cellsnum
                if count_marks[i] == 0
                    return i
                elseif count_marks[i] == 1 && s < 0
                    s = i
                end
            end
            return s
        end
        for lld in 1:lld_cellsnum
            as = []
            for ld in ld_bounds[lld, :].nzind
                push!(as, (ld, angles_fn(lld, ld)))
            end
            sort!(as, lt=(a,b)->a[2]<b[2])
            as = map(a->a[1], as)
            angles[lld] = as
        end

        function nextprev(lld::Int64, ld::Int64, norp)
            as = angles[lld]
            #ne = findfirst(as, ld)  (findfirst(isequal(v), A), 0)[1]
            ne = (findfirst(isequal(ld), as), 0)[1]
            while true
                ne += norp
                if ne > length(as)
                    ne = 1
                elseif ne < 1
                    ne = length(as)
                end

                if count_marks[as[ne]] < 2
                    break
                end
            end
            as[ne]
        end

        while (sigma = get_seed_cell()) > 0

            if verbose
                print(Int(floor(50 * sum(count_marks) / ld_cellsnum)), "%\r")
            end

            c_ld = spzeros(Int8, ld_cellsnum)
            if count_marks[sigma] == 0
                c_ld[sigma] = 1
            else
                c_ld[sigma] = -dir_marks[sigma]
            end
            c_lld = ld_bounds*c_ld
            while c_lld.nzind != []
                corolla = spzeros(Int8, ld_cellsnum)
                for tau in c_lld.nzind
                    b_ld = ld_bounds[tau, :]
                    pivot = intersect(c_ld.nzind, b_ld.nzind)[1]
                    adj = nextprev(tau, pivot, sign(-c_lld[tau]))
                    corolla[adj] = c_ld[pivot]
                    if b_ld[adj] == b_ld[pivot]
                        corolla[adj] *= -1
                    end
                end
                c_ld += corolla
                c_lld = ld_bounds*c_ld
            end
            map(s->count_marks[s] += 1, c_ld.nzind)
            map(s->dir_marks[s] = c_ld[s], c_ld.nzind)
            d_bounds = [d_bounds c_ld]

        end

        return d_bounds

    end

    return _minimal_cycles
end

function minimal_2cycles(V::Lar.Points, EV::Lar.ChainOp)

    function edge_angle(v::Int, e::Int)
        edge = EV[e, :]
        v2 = setdiff(edge.nzind, [v])[1]
        x, y = V[v2, :] - V[v, :]
        return atan(y, x)
    end

    for i in 1:EV.m
        j = min(EV[i,:].nzind...)
        EV[i, j] = -1
    end
    VE = convert(Lar.ChainOp, SparseArrays.transpose(EV))

    EF = minimal_cycles(edge_angle)(V, VE)

    return convert(Lar.ChainOp, SparseArrays.transpose(EF))
end

"""
    test_planar_arrangement(V::Points, EV::ChainOp, [sigma::Chain], [return_edge_map::Bool], [multiproc::Bool])

Compute the arrangement on the given cellular complex 1-skeleton in 2D.

A cellular complex is arranged when the intersection of every possible pair of cell
of the complex is empty and the union of all the cells is the whole Euclidean space.
The basic method of the function without the `sigma`, `return_edge_map` and `multiproc` arguments
returns the full arranged complex `V`, `EV` and `FE`.

## Additional arguments:
- `sigma::Chain`: if specified, `planar_arrangement` will delete from the output every edge and face outside this cell. Defaults to an empty cell.
- `return_edge_map::Bool`: makes the function return also an `edge_map` which maps the edges of the imput to the one of the output. Defaults to `false`.
- `multiproc::Bool`: Runs the computation in parallel mode. Defaults to `false`.
"""
function test_planar_arrangement(
	V::Lar.Points,
	copEV::Lar.ChainOp;
	sigma::Lar.Chain=spzeros(Int8, 0),
	return_edge_map::Bool=false,
	multiproc::Bool=false)

	V,copEV,sigma,return_edge_map,edge_map,rV, rEV = splitedges(
							V,copEV,sigma,return_edge_map,multiproc)
    V, copEV = rV, rEV
    V, copEV = Lar.Arrangement.merge_vertices!(V, copEV, edge_map)

	V,copEV,sigma,return_edge_map,edge_map = deleteoutsidesigma(
							V,copEV,sigma,return_edge_map,edge_map )

    bicon_comps = Lar.Arrangement.biconnected_components(copEV)

    if isempty(bicon_comps)
        println("No biconnected components found.")
        if (return_edge_map)
            return (nothing, nothing, nothing, nothing)
        else
            return (nothing, nothing, nothing)
        end
    end

    edges = sort(union(bicon_comps...))
    todel = sort(setdiff(collect(1:size(copEV,1)), edges))

    for i in reverse(todel)
        for row in edge_map

            filter!(x->x!=i, row)

            for j in 1:length(row)
                if row[j] > i
                    row[j] -= 1
                end
            end
        end
    end

    V, copEV = Lar.delete_edges(todel, V, copEV)

    bicon_comps = Lar.Arrangement.biconnected_components(copEV)

    n = size(bicon_comps, 1)
    shells = Array{Lar.Chain, 1}(undef, n)
    boundaries = Array{Lar.ChainOp, 1}(undef, n)
    EVs = Array{Lar.ChainOp, 1}(undef, n)
    for p in 1:n
        ev = copEV[sort(bicon_comps[p]), :]
        fe = minimal_2cycles(V, ev)
        shell_num = Lar.Arrangement.get_external_cycle(
        	V, ev, fe)

        EVs[p] = ev
        tokeep = setdiff(1:fe.m, shell_num)
        boundaries[p] = fe[tokeep, :]
        shells[p] = fe[shell_num, :]
    end

    shell_bboxes = []
    for i in 1:n
        vs_indexes = (abs.(EVs[i]')*abs.(shells[i])).nzind
        push!(shell_bboxes, Lar.bbox(V[vs_indexes, :]))
    end

    containment_graph = Lar.Arrangement.pre_containment_test(shell_bboxes)
    containment_graph =  Lar.Arrangement.prune_containment_graph(n, V, EVs, shells, containment_graph)

     Lar.Arrangement.transitive_reduction!(containment_graph)

    copEV, FE = Lar.Arrangement.cell_merging(
    	n, containment_graph, V, EVs, boundaries, shells, shell_bboxes)

    if (return_edge_map)
        return V, copEV, FE, edge_map
    else
        return V, copEV, FE
    end
end

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

V,EV = Lar.randomcuboids(10, .75)
V = Plasm.normalize(V,flag=true)
Plasm.view(Plasm.numbering(.05)((V,[[[k] for k=1:size(V,2)], EV])))

W = convert(Lar.Points, V')
cop_EV = Lar.coboundary_0(EV::Lar.Cells)
cop_EW = convert(Lar.ChainOp, cop_EV)
V, copEV, copFE = test_planar_arrangement(W::Lar.Points, cop_EW::Lar.ChainOp;
	multiproc=false)

triangulated_faces = Lar.triangulate2D(V, [copEV, copFE])
V = convert(Lar.Points, V')
FVs = convert(Array{Lar.Cells}, triangulated_faces)
Plasm.viewcolor(V::Lar.Points, FVs::Array{Lar.Cells})

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


text:
# Poisson distribution

export Poisson
import Base
using SpecialFunctions: logfactorial

@parameterized Poisson(λ) ≪ CountingMeasure(ℤ[0:∞])

Base.eltype(::Type{P}) where {P<:Poisson} = Int

function logdensity(d::Poisson{(:λ,)}, y)
    λ = d.λ
    return y * log(λ) - λ - logfactorial(y)
end

function logdensity(d::Poisson{(:logλ,)}, y)
    return y * d.logλ - exp(d.logλ) - logfactorial(y)
end

asparams(::Type{<:Poisson}, ::Val{:λ}) = asℝ₊
asparams(::Type{<:Poisson}, ::Val{:logλ}) = asℝ

sampletype(::Poisson) = Int

Base.rand(rng::AbstractRNG, T::Type, d::Poisson{(:λ,)}) = rand(rng, Dists.Poisson(d.λ))
Base.rand(rng::AbstractRNG, T::Type, d::Poisson{(:logλ,)}) = rand(rng, Dists.Poisson(exp(d.logλ)))

≪(::Poisson, ::IntegerRange{lo,hi}) where {lo, hi} = lo ≤ 0 && isinf(hi)


text:
using FvCFD

# Creates a simple sample cell from Moukalled with one fully-defined cell and a centroid-only neighbour
# For laplacian test cases
function meshMoukalled813()
    f1Pts = [
        [ 2.5, 4.0, 1.0],
        [ 3.5, 1.5, 1.0],
        [ 3.5, 1.5, 0.0],
        [ 2.5, 4.0, 0.0]
    ]
    f2Pts = [
        [ 0.5, 3.0, 1.0],
        [ 2.5, 4.0, 1.0],
        [ 2.5, 4.0, 0.0],
        [ 0.5, 3.0, 0.0]
    ]
    f3Pts = [
        [ 0.0, 1.0, 1.0],
        [ 0.5, 3.0, 1.0],
        [ 0.5, 3.0, 0.0],
        [ 0.0, 1.0, 0.0]
    ]
    f4Pts = [
        [ 2.0, 0.0, 1.0],
        [ 0.0, 1.0, 1.0],
        [ 0.0, 1.0, 0.0],
        [ 2.0, 0.0, 0.0]
    ]
    f5Pts = [
        [ 2.0, 0.0, 1.0],
        [ 2.0, 0.0, 0.0],
        [ 3.5, 1.5, 0.0],
        [ 3.5, 1.5, 1.0]
    ]
    f6Pts = [
        [ 2.5, 4.0, 1.0],
        [ 0.5, 3.0, 1.0],
        [ 0.0, 1.0, 1.0],
        [ 2.0, 0.0, 1.0],
        [ 3.5, 1.5, 1.0],
    ]
    f7Pts = [
        [ 2.5, 4.0, 0.0],
        [ 3.5, 1.5, 0.0],
        [ 2.0, 0.0, 0.0],
        [ 0.0, 1.0, 0.0],
        [ 0.5, 3.0, 0.0]
    ]
    points = [
        [ 2.5, 4.0, 0.0],
        [ 3.5, 1.5, 0.0],
        [ 2.5, 4.0, 1.0],
        [ 3.5, 1.5, 1.0],
        [ 2.0, 0.0, 0.0],
        [ 0.0, 1.0, 0.0],
        [ 0.0, 1.0, 1.0],
        [ 2.0, 0.0, 1.0],
        [ 0.5, 3.0, 0.0],
        [ 0.5, 3.0, 1.0]
    ]
    cells = []

    facePoints = [ f1Pts, f2Pts, f3Pts, f4Pts, f5Pts, f6Pts, f7Pts ]
    fAVecs = Array{Array{Float64, 1}}(undef, 7)
    fCenters = Array{Array{Float64, 1}}(undef, 7)
    for i in 1:7
        fAVec, fCenter = FvCFD.faceAreaCentroid(facePoints[i])
        fAVecs[i] = fAVec
        fCenters[i] = fCenter
    end
    faces = [ [ 1,2 ] ]
    boundaryFaces = []

    cellVol, cellCentroid = FvCFD.cellVolCentroid(points, fAVecs, fCenters)
    cVols = [ cellVol, 1 ]
    cCenters = [ cellCentroid, [4.25, 3.5, 0.5 ]]

    return [ cells, cVols, cCenters, faces, fAVecs, fCenters, boundaryFaces ]
end

max_stars_repo_name: henrystoldt/fvCFD


text:
include("DOM.jl")
using .DOM

include("commands.jl")
include("objects.jl")

# Utils

function tohtml(m::MIME"text/html", x)
  HTML() do io
    writemime(io, m, x)
  end
end

function tohtml(m::MIME"text/plain", x)
  rep = stringmime(m, x)
  lines = split(rep, "\n")
  html = span(".text", rep)
  length(lines) == 1 && length(lines[1]) ≤ 50 ?
    Collapsible(html) :
    Collapsible(strong(summary(x)), html)
end

function tohtml(m::MIME"image/png", img)
  HTML() do io
    print(io, """<img src="data:image/png;base64,""")
    print(io, stringmime(m, img))
    print(io, "\" />")
  end
end

function tohtml(m::MIME"image/svg+xml", img)
  HTML() do io
     writemime(io, m, img)
  end
end

# Display infrastructure

function bestmime(val)
  for mime in ("text/html", "image/svg+xml", "image/png", "text/plain")
    mimewritable(mime, val) && return MIME(symbol(mime))
  end
  error("Cannot display $val.")
end

displayinline(x::HTML) = x

# Catch-all fallback
displayinline(x) = tohtml(bestmime(x), x)

function applydisplayinline(x)
  while (x′ = displayinline(x)) !== x
    x = x′
  end
  return x
end

displayinline!(x, opts) =
  displayinline!(displayinline(x), opts)

displayinline!(html::HTML, opts) =
  showresult(stringmime("text/html", html), opts, html=true, under=true)

# Light Table's Console as a display

type LTConsole <: Display end

import Base: display, writemime

function display(d::LTConsole, m::MIME"text/plain", x)
  console(stringmime(m, x))
end

function display(d::LTConsole, m::MIME"text/html", x)
  console(stringmime(m, x), html = true)
end

display(d::LTConsole, m::MIME"image/png", x) = display(d, tohtml(m, x))

display(d::LTConsole, x) = display(d, bestmime(x), x)


text:
module Schwarz

include("constants.jl")
include("potentials.jl")
include("rotations.jl")
include("orbits.jl")
include("observables.jl")
# include("fit.jl")

using .Constants
c = Constants
using .Observables: KinematicData
using .Potentials: PowerLawPotential, mass, density, potential
using .Orbits: Orbit, BinnedOrbit, OrbitLibrary, BinnedLibrary, integrals
# using .Optimize: fit

export c, mass, density, potential, PowerLawPotential, Orbit, OrbitLibrary,
    BinnedOrbit, BinnedLibrary, KinematicData

end


text:
# Copyright 2017 Gandalf Software, Inc., Scott P. Jones
# Licensed under MIT License, see LICENSE.md

# Todo, rewrite this as a parameterized type, with basetype, swapped, aligned
# Also, add UInt64 and UInt128

abstract type AccessType end
primitive type UInt16_U  <: AccessType 16 end # Unaligned 16-bit word
primitive type UInt32_U  <: AccessType 32 end # Unaligned 32-bit word
primitive type UInt16_S  <: AccessType 16 end # Swapped 16-bit word
primitive type UInt32_S  <: AccessType 32 end # Swapped 32-bit word
primitive type UInt16_US <: AccessType 16 end # Unaligned & swapped 16-bit word
primitive type UInt32_US <: AccessType 32 end # Unaligned & swapped 32-bit word

basetype(::Type{UInt16})    = UInt16
basetype(::Type{UInt16_U})  = UInt16
basetype(::Type{UInt16_S})  = UInt16
basetype(::Type{UInt16_US}) = UInt16

basetype(::Type{UInt32})    = UInt32
basetype(::Type{UInt32_U})  = UInt32
basetype(::Type{UInt32_S})  = UInt32
basetype(::Type{UInt32_US}) = UInt32

# Provide the correct type to access as swapped
swappedtype(::Type{UInt16})    = UInt16_S
swappedtype(::Type{UInt16_U})  = UInt16_US
swappedtype(::Type{UInt16_S})  = UInt16
swappedtype(::Type{UInt16_US}) = UInt16_U

swappedtype(::Type{UInt32})    = UInt32_S
swappedtype(::Type{UInt32_U})  = UInt32_US
swappedtype(::Type{UInt32_S})  = UInt32
swappedtype(::Type{UInt32_US}) = UInt32_U

# Provide the correct type to access if aligned
alignedtype(::Type{UInt16})    = UInt16
alignedtype(::Type{UInt16_U})  = UInt16
alignedtype(::Type{UInt16_S})  = UInt16_S
alignedtype(::Type{UInt16_US}) = UInt16_S

alignedtype(::Type{UInt32})    = UInt32
alignedtype(::Type{UInt32_U})  = UInt32
alignedtype(::Type{UInt32_S})  = UInt32_S
alignedtype(::Type{UInt32_US}) = UInt32_S

@inline _ul(pnt::Ptr{T}, shift) where {T} =
    (unsafe_load(reinterpret(Ptr{UInt8}, pnt))%basetype(T)) << shift

@inline _load(pnt::Ptr{T}) where {T<:Union{UInt16_US,UInt32_US}} =
    bswap(unsafe_load(reinterpret(Ptr{basetype(T)}, pnt)))

@inline _load(pnt::Ptr{UInt16_U})  = _ul(pnt, 0) | _ul(pnt+1, 8)
@inline _load(pnt::Ptr{UInt16_US}) = _ul(pnt, 8) | _ul(pnt+1, 0)
@inline _load(pnt::Ptr{UInt32_U})  = _ul(pnt, 0) | _ul(pnt+1, 8) | _ul(pnt+2, 16) | _ul(pnt+3, 24)
@inline _load(pnt::Ptr{UInt32_US}) = _ul(pnt, 24) | _ul(pnt+1, 16) | _ul(pnt+2, 8) | _ul(pnt+3, 0)

@inline unsafe_load(pnt::Ptr{T}) where {T<:AccessType} = _load(pnt)
@inline unsafe_load(pnt::Ptr{T}, i) where {T<:AccessType} = _load(pnt + (i-1)*sizeof(basetype(T)))


text:
module FluxTools

using Functors: @functor

export Model, @model

include("model.jl")
include("macros.jl")

end


text:
# * Orbital integrals

"""
    OrbitalIntegral{N}

Abstract type for integrals of rank `N` of orbitals, whose values need
to be recomputed every time the orbitals are updated. Rank 0
corresponds to a scalar value, rank 1 to a diagonal matrix, etc.
"""
abstract type OrbitalIntegral{N,aO<:AbstractOrbital,bO<:AbstractOrbital,T} end

Base.iszero(::OrbitalIntegral) = false

# ** Zero integral
struct ZeroIntegral{aO,bO,T} <: OrbitalIntegral{0,aO,bO,T} end

integral_value(::ZeroIntegral{aO,bO,T}) where {aO,bO,T} = zero(T)

SCF.update!(::ZeroIntegral, ::Atom; kwargs...) = nothing

# ** Orbital overlap integral

"""
    OrbitalOverlapIntegral(a, b, av, bv, value)

Represents the orbital overlap integral `⟨a|b⟩`, for orbitals `a` and
`b`, along with `view`s of their radial orbitals `av` and `bv` and the
current `value` of the integral.
"""
mutable struct OrbitalOverlapIntegral{aO,bO,T,OV} <: OrbitalIntegral{0,aO,bO,T}
    a::aO
    b::bO
    av::OV
    bv::OV
    value::T
end

function OrbitalOverlapIntegral(a::aO, b::bO, atom::Atom{T}) where {aO,bO,T}
    oo = OrbitalOverlapIntegral(a, b, view(atom, a), view(atom, b), zero(T))
    SCF.update!(oo)
    oo
end

Base.show(io::IO, oo::OrbitalOverlapIntegral) =
    write(io, "⟨$(oo.a)|$(oo.b)⟩")

"""
    SCF.update!(oo::OrbitalOverlapIntegral)

Update the value of the integral `oo`.
"""
function SCF.update!(oo::OrbitalOverlapIntegral; kwargs...)
    oo.value = materialize(applied(*, oo.av', oo.bv))[1]
end

"""
    SCF.update!(oo::OrbitalOverlapIntegral)

Update the value of the integral `oo` with respect to `atom`.
"""
function SCF.update!(oo::OrbitalOverlapIntegral, atom::Atom; kwargs...)
    oo.av = view(atom, oo.a)
    oo.bv = view(atom, oo.b)
    SCF.update!(oo; kwargs...)
end

integral_value(oo::OrbitalOverlapIntegral) = oo.value

# ** Operator matrix element

"""
    OperatorMatrixElement(a, b, Â, coeff, value)

Represents the matrix element `coeff*⟨a|Â|b⟩`, for the operator `Â`
and the orbitals `a` and `b`, along with the current `value` of the
integral. Typically, `Â` is a radial part of an operator and `coeff`
is the associated angular coefficient; `coeff` can be of any type
convertible to a scalar.
"""
mutable struct OperatorMatrixElement{aO,bO,OV,RO,T,Coeff} <: OrbitalIntegral{0,aO,bO,T}
    a::aO
    b::bO
    av::OV
    bv::OV
    Â::RO
    coeff::Coeff
    value::T
end

function OperatorMatrixElement(a::aO, b::bO, Â::RO, atom::Atom{T}, coeff) where {aO,bO,T,RO}
    ome = OperatorMatrixElement(a, b, view(atom, a), view(atom, b),
                                Â, coeff, zero(T))
    SCF.update!(ome)
    ome
end

function OperatorMatrixElement(a::aO, b::bO, A::M, atom::Atom{T}, coeff) where {aO,bO,T,
                                                                                M<:AbstractMatrix}
    R = radial_basis(atom)
    OperatorMatrixElement(a, b, applied(*,R,A,R'), atom, coeff)
end

function Base.show(io::IO, ome::OperatorMatrixElement)
    write(io, "⟨$(ome.a)|")
    show(io, ome.Â)
    write(io, "|$(ome.b)⟩")
end

function SCF.update!(ome::OperatorMatrixElement{aO,bO,OV,RO,T}; kwargs...) where {aO,bO,OV,RO,T}
    # NB: It is assumed that ome.Â, if necessary, is updated /before/
    # SCF.update!(ome) is called.
    #
    # TODO: This is not particularly efficient, since it allocates a
    # temporary vector; rewrite using applied(*, ...).
    b = ome.bv
    tmp = similar(b)
    tmp.args[2] .= zero(T)
    materialize!(MulAdd(one(T), ome.Â, b, zero(T), tmp))
    ome.value = convert(T,ome.coeff)*materialize(applied(*, ome.av', tmp))
end

function SCF.update!(ome::OperatorMatrixElement, atom::Atom; kwargs...)
    ome.av = view(atom, ome.a)
    ome.bv = view(atom, ome.b)
    SCF.update!(ome; kwargs...)
end

integral_value(ome::OperatorMatrixElement) = ome.value

# ** Hartree–Fock potentials

const HFPotentialOperator{T,B<:Basis} = RadialOperator{T,B,Diagonal{T,<:AbstractVector{T}}}

"""
    HFPotential(k, a, b, av, bv, V̂, poisson)

Represents the `k`:th multipole exansion of the Hartree–Fock potential
formed by orbitals `a` and `b` (`av` and `bv` being `view`s of their
corresponding radial orbitals). `V̂` is the resultant one-body
potential formed, which can act on a third orbital and `poisson`
computes the potential by solving Poisson's problem.
"""
mutable struct HFPotential{kind,aO,bO,T,
                           OV<:RadialOrbital,
                           RO<:HFPotentialOperator{T},P<:AbstractPoissonProblem} <: OrbitalIntegral{1,aO,bO,T}
    k::Int
    a::aO
    b::bO
    av::OV
    bv::OV
    V̂::RO
    poisson::P
end
HFPotential(kind::Symbol, k::Int, a::aO, b::bO, av::OV, bv::OV, V̂::RO, poisson::P) where {aO,bO,T,OV,RO<:RadialOperator{T},P} =
    HFPotential{kind,aO,bO,T,OV,RO,P}(k, a, b, av, bv, V̂, poisson)

get_poisson(::CoulombInteraction{Nothing}, args...; kwargs...) =
    PoissonProblem(args...; kwargs...)

get_poisson(g::CoulombInteraction{<:AbstractQuasiMatrix}, args...; kwargs...) =
    AsymptoticPoissonProblem(args..., g.o; kwargs...)

function HFPotential(kind::Symbol, k::Int, a::aO, b::bO, atom::Atom{T}, g::CoulombInteraction; kwargs...) where {aO,bO,T}
    av, bv = view(atom, a), view(atom, b)
    R = av.args[1]
    D = Diagonal(Vector{T}(undef, size(R,2)))
    D.diag .= zero(T)
    V̂ = applied(*, R, D, R')
    poisson = get_poisson(g, k, av, bv; w′=applied(*, R, D.diag), kwargs...)
    update!(HFPotential(kind, k, a, b, av, bv, V̂, poisson), atom)
end

Base.convert(::Type{HFPotential{kind,aO₁,bO₁,T,OV,RO,P}},
             hfpotential::HFPotential{kind,aO₂,bO₂,T,OV,RO,P}) where {kind,aO₁,bO₁,aO₂,bO₂,T,OV,RO,P} =
                 HFPotential{kind,aO₁,bO₁,T,OV,RO,P}(hfpotential.k,
                                                     hfpotential.a, hfpotential.b,
                                                     hfpotential.av, hfpotential.bv,
                                                     hfpotential.V̂, hfpotential.poisson)

# *** Direct potential

"""
    DirectPotential

Special case of [`HFPotential`](@ref) for the direct interaction, in
which case the potential formed from two orbitals can be precomputed
before acting on a third orbital.
"""
const DirectPotential{aO,bO,T,OV,RO,P} = HFPotential{:direct,aO,bO,T,OV,RO,P}

Base.show(io::IO, Y::DirectPotential) =
    write(io, "r⁻¹×Y", to_superscript(Y.k), "($(Y.a), $(Y.b))")

"""
    SCF.update!(p::DirectPotential)

Update the direct potential `p` by solving the Poisson problem with
the current values of the orbitals forming the mutual density.
"""
function SCF.update!(p::DirectPotential{aO,bO,T,OV,RO,P}; kwargs...) where {aO,bO,T,OV,RO,P}
    p.poisson(p.av .⋆ p.bv; kwargs...)
    p
end

"""
    SCF.update!(p::DirectPotential, atom::Atom)

Update the direct potential `p` by solving the Poisson problem with
the current values of the orbitals of `atom` forming the mutual
density.
"""
function SCF.update!(p::DirectPotential{aO,bO,T,OV,RO,P}, atom::Atom; kwargs...) where {aO,bO,T,OV,RO,P}
    p.av = view(atom, p.a)
    p.bv = view(atom, p.b)
    SCF.update!(p; kwargs...)
end

"""
    materialize!(ma::MulAdd{<:Any, <:Any, <:Any, T, <:DirectPotential, Source, Dest})

Materialize the lazy multiplication–addition of the type `y ←
α*V̂*x + β*y` where `V̂` is a [`DirectPotential`](@ref) (with a
precomputed direct potential computed via `SCF.update!`) and `x` and
`y` are [`RadialOrbital`](@ref)s.
"""
LazyArrays.materialize!(ma::MulAdd{<:Any, <:Any, <:Any, T, <:DirectPotential, Source, Dest}) where {T,Source,Dest} =
    materialize!(MulAdd(ma.α, ma.A.V̂, ma.B, ma.β, ma.C))

# *** Exchange potential

"""
    ExchangePotential

Special case of [`HFPotential`](@ref) for the exchange interaction, in
which case the potential is formed from the orbital acted upon, along
with another orbital, and then applied to a third orbital. Thus this
potential *cannot* be precomputed, but must be recomputed every time
the operator is applied. This makes this potential expensive to handle
and the number of times it is applied should be minimized, if possible.
"""
const ExchangePotential{aO,bO,T,OV,RO,P} = HFPotential{:exchange,aO,bO,T,OV,RO,P}

Base.show(io::IO, Y::ExchangePotential) =
    write(io, "|$(Y.b)⟩r⁻¹×Y", to_superscript(Y.k), "($(Y.a), ●)")

SCF.update!(p::ExchangePotential; kwargs...) = p

function SCF.update!(p::ExchangePotential, atom::Atom; kwargs...)
    p.av = view(atom, p.a)
    p.bv = view(atom, p.b)
    p
end

"""
    materialize!(ma::MulAdd{<:Any, <:Any, <:Any, T, <:ExchangePotential, Source, Dest})

Materialize the lazy multiplication–addition of the type `y ← α*V̂*x +
β*y` where `V̂` is a [`ExchangePotential`](@ref) (by solving the
Poisson problem with `x` as one of the constituent source orbitals in
the mutual density) and `x` and `y` are [`RadialOrbital`](@ref)s.
"""
function LazyArrays.materialize!(ma::MulAdd{<:Any, <:Any, <:Any, T, <:ExchangePotential, Source, Dest}) where {T,Source,Dest}
    p = ma.A
    p.poisson(p.av .⋆ ma.B) # Form exchange potential from conj(p.a)*b
    # Act with the exchange potential on p.bv
    materialize!(MulAdd(ma.α, p.V̂, p.bv, ma.β, ma.C))
end

# * Source terms

"""
    SourceTerm(operator, source_orbital, ov)

The point of `SourceTerm` is to implement inhomogeneous terms that
contribute to the equation for an orbital, and whose input is some
other `source_orbital`. This kind of term appears in
multi-configurational problems.
"""
mutable struct SourceTerm{QO,O,OV}
    operator::QO
    source_orbital::O
    ov::OV
end

Base.show(io::IO, st::SourceTerm) = write(io, "SourceTerm($(st.operator)|$(st.source_orbital)⟩)")

Base.iszero(::SourceTerm) = false
Base.similar(st::SourceTerm) = similar(st.ov)

function Base.copyto!(dest::Mul{<:Any,<:Tuple{<:AbstractQuasiMatrix,<:AbstractArray{<:Any,N}}},
                      src::Mul{<:Any,<:Tuple{<:AbstractQuasiMatrix,<:AbstractArray{<:Any,N}}}) where N
    d = last(dest.args)
    s = last(src.args)
    copyto!(IndexStyle(d), d, IndexStyle(s), s)
    dest
end

update!(st::SourceTerm; kwargs...) = nothing

function update!(st::SourceTerm, atom::Atom; kwargs...)
    st.ov = view(atom, st.source_orbital)
end

# These are source terms that do not depend on the atom (e.g. external
# source term, or a constant orbital).
update!(::SourceTerm{<:IdentityOperator,<:AbstractString}, ::Atom) = nothing

LazyArrays.materialize!(ma::MulAdd{<:Any, <:Any, <:Any, T, <:SourceTerm, Source, Dest}) where {T,Source,Dest} =
    materialize!(MulAdd(ma.α, ma.A.operator, ma.A.ov, ma.β, ma.C))

function LazyArrays.materialize!(ma::MulAdd{<:Any, <:Any, <:Any, T, <:IdentityOperator{1}, Source, Dest}) where {T,Source,Dest}
    isone(ma.β) || lmul!(ma.β, ma.C.args[2])
    BLAS.axpy!(ma.α, ma.B.args[2], ma.C.args[2])
end

# * Shift terms

"""
    ShiftTerm(λ)

The point of `ShiftTerm` is to implement an overall energy shift of
the Hamiltonian.
"""
struct ShiftTerm{T}
    shift::UniformScaling{T}
end

Base.iszero(::ShiftTerm) = false

LazyArrays.materialize!(ma::MulAdd{<:Any, <:Any, <:Any, T, <:ShiftTerm, Source, Dest}) where {T,Source,Dest} =
    BLAS.axpy!(ma.α*ma.A.shift.λ, ma.B.args[2], ma.C.args[2])


text:
# BinaryEntropy loss: f(x,y)=sum(x.*log(x./y)+(1.-x).*log((1.-x)./(1.-y)))/prod(size(x))
# It's often a bad idea to use this function since it is numerically unstable.
# BinaryKullbackLeiblerLossXsigmaY is usually better.

#FBinaryEntropyLoss(x::Array,y::Array)=([sum(x.*log(x./y)+(1.-x).*log((1.-x)./(1.-y)))/length(x)],[])

FBinaryEntropyLoss(x::Array,y::Array)=([sum(x.*log(x./y)+(1.-x).*log((1.-x)./(1.-y)))/length(x)],((y-x)./(y.*(1.-y))/length(x)))

#FBinaryEntropyLoss_inplace(value,aux,x::Array,y::Array)=copy!(value,sum(x.*log(x./y)+(1.-x).*log((1.-x)./(1.-y)))/length(x))
function FBinaryEntropyLoss_inplace(value,aux,x::Array,y::Array)
    copy!(value,sum(x.*log(x./y)+(1.-x).*log((1.-x)./(1.-y)))/length(x))
    copy!(aux,((y-x)./(y.*(1.-y)))/length(x))
end

function DBinaryEntropyLoss(derivativeIDX,f_c,faux_c,grad_c,grad_n,x,y)
    if derivativeIDX==1
        axpy!(grad_c[1],log(x.*(1.-y)./(y.*(1.-x)))/length(x),grad_n)
    elseif derivativeIDX==2
        #axpy!(grad_c[1],((y-x)./(y.*(1.-y)))/length(x),grad_n)
        axpy!(grad_c[1],faux_c,grad_n)
    end
end

if PROC=="GPU"
    FBinaryEntropyLoss(x::CudaArray,y::CudaArray)=(binaryentropy(x,y),nothing)
    FBinaryEntropyLoss_inplace(value,aux,x::CudaArray,y::CudaArray)=copy!(value,binaryentropy(x,y))

    function DBinaryEntropyLoss(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray,y::CudaArray)
        if derivativeIDX==1
            DXbinaryentropy!(x,y,grad_c,grad_n)
        elseif derivativeIDX==2
            DYbinaryentropy!(x,y,grad_c,grad_n)
        end
    end
end

Derivative[FBinaryEntropyLoss]=DBinaryEntropyLoss
Inplace[FBinaryEntropyLoss]=FBinaryEntropyLoss_inplace

BinaryKullbackLeiblerLoss(nx,ny)=ADnode(FBinaryEntropyLoss,[nx ny])
export BinaryKullbackLeiblerLoss

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
function set_wdot(Qbase, cellxmax, cellymax, Rd, nch, nre, nval)

    wdot = zeros(cellxmax, cellymax, nval)
    npre = nval - nch

    mw_N2_N = [28e-3, 14e-3]
    avgdro  = 6.022*10^23

    # N2 + N2 --> N + N + N2
    # N2 + N  --> N + N + N
    fmol = [[2, 0],
            [1, 1]]
    # N2 + N2 <-- N + N + N2
    # N2 + N  <-- N + N + N
    bmol = [[1, 2],
            [0, 3]]

    for i in 2:cellxmax-1
        for j in 2:cellymax-1

            # Initial setting
            total = 0.0
            rho_mw = zeros(nch)
            for s in 1:nch
                rho_mw[s] = Qbase[i,j,npre+s] / mw_N2_N[s]   # [mol/m3] = [kg/m3] / [kg/mol]
                total     = total + rho_mw[s]
            end
            total = log10(total*avgdro)

            if total <= 20.0
                nsb = 1
            elseif total <= 21.0
                nsb = 2
            elseif total <= 22.0
                nsb = 3
            elseif total <= 23.0
                nsb = 4
            elseif total <= 24.0
                nsb = 5
            else
                nsb = 6
            end

            T   = Qbase[i,j,npre]/(Qbase[i,j,1]*Rd)

            keq = k_eq(T, nsb, nre)                 # [mol/(cm^3 s)]
            kf  = freaction_rate(T, nre)            # [mol/(cm^3 s)]
            kb  = breaction_rate(kf, keq, nre)      # [mol/(cm^3 s)]

            # Forward and backward rate for reaction ir
            lf = ones(nre)
            lb = ones(nre)
            for r in 1:nre
                for s in 1:nch
                    # 並行定数の単位に合わせるため，モル濃度を[mol/cm^3]の単位で計算している
                    lf[r] = lf[r]*1.0e-6*rho_mw[s]^fmol[r][s]
                    lb[r] = lb[r]*1.0e-6*rho_mw[s]^bmol[r][s]
                end
                lf[r] = lf[r]*kf[r]
                lb[r] = lb[r]*kb[r]
            end

            for s in 1:nch
                wdot[i,j,npre+s] = 0.0
                for r in 1:nre
                    temp = (fmol[r][s] - bmol[r][s])*(lf[r] - lb[r])
                    wdot[i,j,npre+s] += temp
                end
                # モル濃度[mol/cm^3]から[mol/m^3]に変換するため，1.0e6を乗じている，
                wdot[i,j,npre+s] = 1.0e6 * wdot[i,j,npre+s] * mw_N2_N[s]
            end
        end
    end
    return wdot
end

function chemical_value(Qbase, cellxmax, cellymax, Rd, nval, nch)

    chmu        = zeros(cellxmax, cellymax)
    chlambda_tr = zeros(cellxmax, cellymax)
    chD         = zeros(cellxmax, cellymax, nch)
    chi         = zeros(cellxmax, cellymax, nch)
    # N2 - N2
    # N2 - N
    # N  - N
    num_nncross = binomial(BigInt(2*nch-1), BigInt(nch)) # 重複組み合わせ nHr = n+r-1Cr
    npre = nval - nch

    mw = [28.0e-3, 14.0e-3]  # 原子量[kg/mol]
    avgdro  = 6.022*10^23    # アボガドロ定数
    kb = 1.380/10^(23)      # ボルツマン定数
    pi = 3.1415              # pi

    ms = zeros(nch)          # 原子一個当たりの質量[kg]
    for s in 1:nch
        ms[s] = mw[s] / avgdro
    end

    for i in 1:cellxmax
        for j in 1:cellymax

            total = 0.0
            for s in 1:nch
                total += Qbase[i,j,npre+s]/mw[s]
            end
            for s in 1:nch
                chi[i,j,s] = (Qbase[i,j,npre+s]/mw[s]) /total
            end

            T = Qbase[i,j,npre]/(Qbase[i,j,1]*Rd)
            p = Qbase[i,j,npre]

            # 対象行列に格納
            # [ N2-N2   N2-N ]
            # [   0.0    N-N ]
            #
            ps11 = collision_cross_section_pisigma11(T, nch)
            ps22 = collision_cross_section_pisigma22(T, nch)

            deltaij1 = zeros(nch, nch)
            deltaij2 = zeros(nch, nch)

            # 以下では行列のすべてを計算しているが、下三角は0なのでキャンセルされる
            for si in 1:nch
                for sj in 1:nch
                    mi = ms[si]
                    mj = ms[sj]

                    temp = 2*mi*mj / (pi*kb*T*(mi+mj))
                    deltaij1[si,sj] = 8/3 * temp^0.5 * ps11[si,sj]
                    deltaij2[si,sj] = 16/5 * temp^0.5 * ps22[si,sj]
                end
            end

            chmu[i,j] = 0.0
            temp      = 0.0
            for si in 1:nch
                for sj in 1:nch
                    temp += chi[i,j,sj] * deltaij2[si,sj]
                end
                chmu[i,j] += ms[si]*chi[i,j,si]/temp
            end

            chlambda_tr[i,j] = 0.0
            temp =0.0
            for si in 1:nch
                for sj in 1:nch
                    alphaij = 1 + (1-ms[si]/ms[sj]) * (0.45-2.54*ms[si]/ms[sj]) / (1+ms[si]/ms[sj])^2
                    temp += alphaij * chi[i,j,sj] * deltaij2[si, sj]
                end
                chlambda_tr[i,j] += chi[i,j,si]/temp
            end
            chlambda_tr[i,j] = 15/4 * kb * chlambda_tr[i,j]

            # 対象の行と列を足して(si,si)を引く
            # [ N2-N2   N2-N ]
            # [   0.0    N-N ]
            #
            #=
            (例) N2
                sum[1,si] + sum[si,1] - [si,si]*2 = [ N2-N2  N2-N ] + [ N2-N2 0.0 ] - [N2-N2]*2
                                                  = [ N2-N ]

            Dijを導入すると0割りが発生するため、inverse Dijを計算している
            =#
            inv_Dij = zeros(nch,nch)
            for si in 1:nch
                for sj in 1:nch
                    inv_Dij[si,sj] = (p*deltaij1[si,sj])/kb*T
                end
            end

            for si in 1:nch
                tempsum = 0
                for sj in 1:nch
                    tempsum += chi[i,j,sj]*inv_Dij[si,sj]
                    tempsum += chi[i,j,sj]*inv_Dij[sj,si]
                end
                tempsum -= chi[i,j,si]*inv_Dij[si,si]*2
                rho  = Qbase[i,j,1]
                rhos = Qbase[i,j,npre+si]
                Cs   = rhos / rho
                if chi[i,j,si] == 0.0 || tempsum == 0.0
                    chD[i,j,si] = 0.0
                else
                    chD[i,j,si] = (1-Cs)/tempsum * Cs/chi[i,j,si]
                end
            end
        end
    end

    return chmu, chlambda_tr, chD, chi
end


text:
@safetestset "MINERVA" begin
    using GlobalMatchingModels, Statistics, Test, Random
    M = [1 1 0; 0 -1 0]
    p = [0, 1, 0]

    s = compute_similarity(M, p)
    @test s == [.5,-1]

    a = compute_activation(M, p)
    @test a == [.125,-1]

    x = compute_intensity(M, p)
    @test x == -0.875

    include("../Examples/MINERVA_Functions.jl")
    Random.seed!(1125435)
    # initialize model
    model = MINERVA(n_features=20, n_traces=60, L=.50)
    # simulate experiment
    sim_data = mapreduce(x->simulate(model, 4, 5), hcat, 1:20000)'
    means = mean(sim_data, dims=1)
    @test means ≈ [.0 0.153 0.31  0.46  0.61  0.77] atol = .01
end


text:
using Documenter, ParquetFiles

makedocs(
	modules=[ParquetFiles],
	sitename="ParquetFiles.jl",
	analytics="UA-132838790-1",
	pages=[
        "Introduction" => "index.md"
    ]
)

deploydocs(
    repo="github.com/queryverse/ParquetFiles.jl.git"
)


text:
using Test

include("dataPreProcessing.jl")
include("testfunctions.jl")

@testset "Test readIBW                   " begin
    # TO DO:
    # ERROR HANDLING
    @test  begin # "Read back written data: small array"
        filename = "Test1.ibw"
        inputData = Float32.(rand(100,3))
        writeIBW(filename,inputData)
        all(reshape(inputData, : , 1) .== readIBW(filename))
    end
    rm("Test1.ibw")
    @test  begin # "Read back written data: big array"
        filename = "Test2.ibw"
        inputData = Float32.(rand(10_000,3))
        writeIBW(filename,inputData)
        all(reshape(inputData, : , 1) .== readIBW(filename))
    end
    rm("Test2.ibw")
    @test begin # "Data supplied as a Float64 matrix"
        filename = "Test3.ibw"
        inputData = Float64.(rand(100,3))
        writeIBW(filename,inputData)
        all(Float32.(reshape(inputData, : , 1)) .== readIBW(filename))
    end
    rm("Test3.ibw")

    ### Here comes a series of regression tests designed on real data files.
    ### Obviously any failing test means only that behavior has changed, not
    ### necessarily that the previous implementation was correct. This is
    ### because the real .IBW file format is "black-box" to me.
    @test begin
        # Regression test 1
        inputVec = readIBW("sampleRealData/MOF_Mg_2_1_0009.ibw")
        length(inputVec) == 254541
    end
    @test begin
        # Regression test 2
        inputVec = readIBW("sampleRealData/MOF_Mg_2_1_0009.ibw")
        inputVec[12345] == -3.404614f-6
    end
    @test begin
        # Regression test 3
        inputVec = readIBW("sampleRealData/MOF_Zn_2_1_0008.ibw")
        length(inputVec) == 254466
    end
    @test begin
        # Regression test 4
        inputVec = readIBW("sampleRealData/MOF_Zn_2_1_0008.ibw")
        inputVec[end] == -2.6649975f-6
    end
    @test begin
        # Regression test 5
        inputVec = readIBW("sampleRealData/PMMA_0013.ibw")
        length(inputVec) == 254961
    end
    @test begin
        # Regression test 6
        inputVec = readIBW("sampleRealData/PMMA_0013.ibw")
        inputVec[1] == -2.3641041f-6
    end
    @test begin
        # Regression test 7
        inputVec = readIBW("sampleRealData/RH45_1_2_0001.ibw")
        length(inputVec) == 4299111
    end
    @test begin
        # Regression test 8
        inputVec = readIBW("sampleRealData/RH45_1_2_0001.ibw")
        inputVec[1337] == 1.5439662f-7
    end
end

@testset "Test IBWtoTXT                  " begin
    @test begin
        filename = "Test1.ibw"
        inputData = Float32.([collect(1:100) collect(101:200) collect(201:300)])
        writeIBW(filename,inputData)
        all(inputData[:,3:-1:2] .== IBWtoTXT(filename))
    end
    rm("Test1.ibw")
    @test begin
        filename = "Test1.ibw"
        inputData = Float32.(rand(100,3))
        writeIBW(filename,inputData)
        all(inputData[:,3:-1:2] .== IBWtoTXT(filename))
    end
    rm("Test1.ibw")
end

@testset "Test offsetAndDriftCompensation" begin
    ############################################
    # First data set:
    #   - No noise
    #   - Linear interpolation between keypoints
    #   -
    sampleRate = 1000.0 # Hz
    #                 TIME   FORCE    DISPLACEMENT
        matOfKeyPoints = [ 0.0    100.0    -200.0        # Starting position
                           5.0    100.0     100.0        # Position of contact
                           6.0   1500.0     110.0        # End of ramp / Start of hold
                          16.0   1500.0     125.0        # End of hold / Start of unload
                          17.0    175.0     105.0        # End of unload / Start of thermal hold
                          47.0    175.0     101.0        # End of thermal hold / Start of retraction
                          48.0    100.0     100.0        # End of retraction / Start of exit
                          50.0    100.0    -100.0]       # End of exit & End of experiment
    xy = generateSyntheticForceDisplacementSignal(matOfKeyPoints,sampleRate)
    xy , basefit , zeroLine , rampStartIdx, xZero = offsetAndDriftCompensation(Float32.(xy))
    @test begin
        isapprox(zeroLine[1] , (matOfKeyPoints[2,2]-matOfKeyPoints[1,2])/(matOfKeyPoints[2,3]-matOfKeyPoints[1,3]), atol = 1e-6)
    end
    @test begin
        isapprox(zeroLine[2] , 0.5*(matOfKeyPoints[1,2]+matOfKeyPoints[2,2]), atol = 1e-6)
    end
    @test begin
        isapprox(rampStartIdx , sampleRate.*matOfKeyPoints[2,1], rtol = 0.01)
    end
    @test begin
        isapprox(Float32(xZero), matOfKeyPoints[2,3], rtol = 0.01)
    end
    ############################################
    # Second data set:
    #   - Noise
    #   - Linear interpolation between keypoints
    #   -
    sampleRate = 1000.0 # Hz
    #                 TIME   FORCE    DISPLACEMENT
        matOfKeyPoints = [ 0.0    100.0    -200.0        # Starting position
                           5.0    100.0     100.0        # Position of contact
                           6.0   1500.0     110.0        # End of ramp / Start of hold
                          16.0   1500.0     125.0        # End of hold / Start of unload
                          17.0    175.0     105.0        # End of unload / Start of thermal hold
                          47.0    175.0     101.0        # End of thermal hold / Start of retraction
                          48.0    100.0     100.0        # End of retraction / Start of exit
                          50.0    100.0    -100.0]       # End of exit & End of experiment
    # xy = generateSyntheticForceDisplacementSignal(matOfKeyPoints,sampleRate, stdOfNoise = 2.0,sizeOfNoise =[0.0 0.4 0.4])
    # xy , basefit , zeroLine , rampStartIdx, xZero = offsetAndDriftCompensation(Float32.(xy))
    # @test begin
    #     isapprox(zeroLine[1] , (matOfKeyPoints[2,2]-matOfKeyPoints[1,2])/(matOfKeyPoints[2,3]-matOfKeyPoints[1,3]), atol = 0.02)
    # end
    # @test begin
    #     isapprox(zeroLine[2] , 0.5*(matOfKeyPoints[1,2]+matOfKeyPoints[2,2]), atol = 10.0)
    # end
    # @test begin
    #     #println(rampStartIdx)
    #     #println(sampleRate.*matOfKeyPoints[2,1])
    #     isapprox(rampStartIdx, sampleRate.*matOfKeyPoints[2,1], atol = 25)
    # end
    # @test begin
    #     isapprox(Float32(xZero), matOfKeyPoints[2,3], rtol = 0.05)
    # end

end

@testset "Test findStartOfHold           " begin
    sampleRate = 1000.0 # Hz
    #                 TIME   FORCE    DISPLACEMENT
    matOfKeyPoints  = [ 0.0    100.0    -200.0        # Starting position
                        5.0    100.0     100.0        # Position of contact
                        6.0   1500.0     110.0        # End of ramp / Start of hold
                        16.0   1500.0     125.0        # End of hold / Start of unload
                        17.0    175.0     105.0        # End of unload / Start of thermal hold
                        47.0    175.0     101.0        # End of thermal hold / Start of retraction
                        48.0    100.0     100.0        # End of retraction / Start of exit
                        50.0    100.0    -100.0]       # End of exit & End of experiment
    xy = generateSyntheticForceDisplacementSignal(matOfKeyPoints,sampleRate)
    xy , ~ , ~ , ~, ~ = offsetAndDriftCompensation(Float32.(xy))
    @test begin
        isapprox(findStartOfHold(xy, "first") , sampleRate*matOfKeyPoints[3,1], atol = 10)
    end
    @test begin
        isapprox(findStartOfHold(xy, "last") , sampleRate*matOfKeyPoints[4,1], atol = 10)
    end
    xy = generateSyntheticForceDisplacementSignal(matOfKeyPoints,sampleRate, stdOfNoise = 2.0, sizeOfNoise =[0.0 0.4 0.4])
    xy , ~ , ~ , ~, ~ = offsetAndDriftCompensation(Float32.(xy))
    @test begin
        isapprox(findStartOfHold(xy, "first") , sampleRate*matOfKeyPoints[3,1], atol = 10)
    end
    @test begin
        isapprox(findStartOfHold(xy, "last") , sampleRate*matOfKeyPoints[4,1], atol = 10)
    end
end

@testset "determineThermalCreep          " begin
    sampleRate = 1000
    holdTimeVals = collect(1:10000) ./ sampleRate
    thermalCreepGenerator(x,y,z) = x .+ y.*holdTimeVals.^z
    @test begin
        xy_hold = thermalCreepGenerator(44.0,4.0,1.1)
        isapprox(4.0 * 1.1*holdTimeVals[end]^(1.1 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),rtol = 0.01)
    end
    @test begin
        xy_hold = thermalCreepGenerator(123.0,4.0,1.1)
        isapprox(4.0 * 1.1*holdTimeVals[end]^(1.1 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),rtol = 0.01)
    end
    @test begin
        xy_hold = thermalCreepGenerator(123.0,0.20,1.1)
        isapprox(0.2 * 1.1*holdTimeVals[end]^(1.1 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),rtol = 0.01)
    end
    @test begin
        xy_hold = thermalCreepGenerator(123.0,4.00,0.8)
        isapprox(4.0 * 0.8*holdTimeVals[end]^(0.8 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),rtol = 0.01)
    end
    @test begin
        xy_hold = thermalCreepGenerator(0.0,4.00,0.8)
        isapprox(4.0 * 0.8*holdTimeVals[end]^(0.8 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),rtol = 0.02)
    end
    @test begin
        xy_hold = thermalCreepGenerator(123.0,0.00,0.8)
        isapprox(0.0 * 0.8*holdTimeVals[end]^(0.8 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),atol = 0.01)
    end
    @test begin
        xy_hold = thermalCreepGenerator(123.0,4.00,0.0)
        isapprox(4.0 * 0.0*holdTimeVals[end]^(0.0 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),atol = 0.01)
    end
end

@testset "determineCreepDuringHold       " begin
    sampleRate = 1000
    holdTimeVals = collect(1:10000) ./ sampleRate
    holdCreepGenerator(x,y,z) = x .+ y.*holdTimeVals.^z
    @test begin
        xy_hold = holdCreepGenerator(44.0,4.0,1.1)
        isapprox(4.0 * 1.1*holdTimeVals[end]^(1.1 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),rtol = 0.01)
    end
    @test begin
        xy_hold = holdCreepGenerator(123.0,4.0,1.1)
        isapprox(4.0 * 1.1*holdTimeVals[end]^(1.1 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),rtol = 0.01)
    end
    @test begin
        xy_hold = holdCreepGenerator(123.0,0.20,1.1)
        isapprox(0.2 * 1.1*holdTimeVals[end]^(1.1 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),rtol = 0.01)
    end
    @test begin
        xy_hold = holdCreepGenerator(123.0,4.00,0.8)
        isapprox(4.0 * 0.8*holdTimeVals[end]^(0.8 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),rtol = 0.01)
    end
    @test begin
        xy_hold = holdCreepGenerator(0.0,4.00,0.8)
        isapprox(4.0 * 0.8*holdTimeVals[end]^(0.8 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),rtol = 0.02)
    end
    @test begin
        xy_hold = holdCreepGenerator(123.0,0.00,0.8)
        isapprox(0.0 * 0.8*holdTimeVals[end]^(0.8 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),atol = 0.01)
    end
    @test begin
        xy_hold = holdCreepGenerator(123.0,4.00,0.0)
        isapprox(4.0 * 0.0*holdTimeVals[end]^(0.0 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),atol = 0.01)
    end
end

@testset "modulusfitter                  " begin
    ctrl = control(false, false)
    indentSet = metaInfoExperimentalSeries("Test2",
                                            50.0,
                                        "hemisphere",
                                        "X",
                                        217.51,
                                        "sampleRealData/LRH250_1_0001_arfun.txt",
                                        "sampleRealData/",
                                        50000)

    # Regression tests
    @test begin
        hp =  hyperParameters( 2000, 1400, "Oliver-Pharr", false, 0, 0, 0.0)
        isapprox(1.8696, modulusfitter( indentSet, hp, ctrl, "MOF_Mg_2_1_0009.ibw"),rtol = 0.05)
    end
    @test begin
        hp =  hyperParameters( 2000, 1400, "Oliver-Pharr", true, 0, 0, 0.0)
        isapprox(1.9, modulusfitter( indentSet, hp, ctrl, "MOF_Mg_2_1_0009.ibw"),rtol = 0.05)
    end
    @test begin
        hp =  hyperParameters( 2000, 1400, "Feng", false, 0, 0, 0.0)
        isapprox( 1.9289, modulusfitter( indentSet, hp, ctrl, "MOF_Mg_2_1_0009.ibw"),rtol = 0.05)
    end
    @test begin
        hp =  hyperParameters( 2000, 1400, "Feng", true, 0, 0, 0.0)
        isapprox( 1.92197, modulusfitter( indentSet, hp, ctrl, "MOF_Mg_2_1_0009.ibw"), rtol = 0.05)
    end
end

max_stars_repo_name: abrandberg/afmIndentation

max_issues_repo_name: abrandberg/afmIndentation

max_forks_repo_name: abrandberg/afmIndentation


text:
function eye(n::Int64)
    res = zeros(n,n)
    for i = 1:n res[i,i] = 1.0 end
    return res
end
#Pkg FillArrays (Eye) ?

function eye(dtype::DataType, n::Int64)
    res = zeros(dtype, n, n)
    I = 1.0 |> dtype
    for i = 1:n res[i,i] = I end
    return res
end

function delta(x::Real; c::Real = 0., η::Real = 0.05)
    num = η
    den = ((x-c)^2 + η^2) * π
    return num/den
end

function gaussian(x::Real, fac::AbstractVector)
    # f(x) = 1/(σ√(2π)) exp(-((x-μ)/σ)^2 /2 )
    μ = fac[1]; σ2 = fac[2]
    σ2 <= 0. ? (@error "σ2 should be positive.") :
    norm = √(2π * σ2)
    m = (x-μ)^2 / σ2
    return 1/norm*exp(-m/2)
end

function multi_gaussian(N::Int64; ms::Vector{Float64}=[-3.,3.], smax::Real=3.)
    fac = rand(N,2)
    fac[:,1] = rand(N) .* (ms[2] - ms[1]) .+ ms[1]
    fac[:,2] = rand(N) .* smax
    function f(x::Real)
        res = 0.
        for i = 1: N
            res += gaussian(x, fac[i,:])
        end
        res/N
    end
    return f, fac
end

function Masubara_freq(n::Int64, β::Real; type::Symbol= :f)
    if type == :b  N = 2n
    elseif type == :f  N = 2n + 1
    else @error "type should be :b for bosons and :f for fermions"
    end
    return N*π/β
end

function Masubara_GF(n::Int64, A::Function, β::Real;
    type::Symbol = :f, Λ::Float64=100., L::Int64=1000000)
    # G(iωn) = 1/2π ∫dΩ A(Ω)/(iωn - Ω)
    ωn = Masubara_freq(n, β, type=type)
    dL = 2*Λ/L; lowlim = -Λ
    Gn = 0.
    for n = 0:L
        ω = lowlim + dL*n
        Gn += A(ω)/(1.0im*ωn-ω) * dL
    end
    return Gn/(2π)
end

function h(z::Number)
    #Conformal mapping(Mobius transform)
    #\bar{C+} -> \bar{D}
    return (z - 1.0im)/(z + 1.0im)
end

function invh(z::Number)
    return 1.0im*(1 + z)/(1 - z)
end

function h1(z::Number, Y::Number)
    #Conformal mapping(Mobius transform)
    #C+ -> D
    return (z - Y)/(z - Y')
end

function invh1(z::Complex, Y::Number)
    return (z*Y'-Y)/(z-1)
end

function ispossemidef(A::Matrix)
    evals = eigvals(A)
    return all(evals .>= 0)
end


text:
using MacroTools: @capture, prewalk

using Calculus: differentiate

using CommonSubexpressions

import FunctionWrappers
import FunctionWrappers: FunctionWrapper

export AnalyticFunction, @analytic

const ScalarFun{T} = FunctionWrapper{T, Tuple{T}}

"""
`F64fun`: `FunctionWrapper` to wrap many different potentials within a single
type. Can be used as  `F64fun(::Function)` or as
`F64fun(::AnalyticFunction)`
"""
const F64fun = ScalarFun{Float64}

"""
take an expression of the form `r -> f(r)` and return the expression
`r -> f'(r)`
"""
function fdiff( ex, ndiff )
   @assert 1 <= ndiff <= 2
   @assert @capture(ex, var_ -> expr_)
   d_ex = differentiate(expr.args[2], var)
   if ndiff == 2
      d_ex = differentiate(d_ex, var)
   end
   d_ex = CommonSubexpressions.cse(d_ex)
   return :(  $var -> $d_ex )
end

"""
auxiliary function to allow expression substitution from `@Analytic`;
see `?@Analytic` for more detail
"""
function substitute(args)
    subs = Dict{Symbol,Union{Expr, Symbol}}()
    for i = 2:length(args)
        @assert @capture(args[i], var_ = sub_)
        subs[var] = sub
    end
    prewalk(expr -> get(subs, expr, expr), args[1])
end

# ================== Analytical Potentials ==========================

"""
`struct AnalyticFunction`: described an analytic function, allowing to
evaluate at least 2 derivatives.

Formally, `AnalyticFunction <: PairPotential`, which simplifies dispatch,
but it an `AnalyticFunction` should normally **not be used as a Calculator**!
This type hierarchy may need to be revisited.

### Usage:
```julia
lj = @analytic r -> r^(-12) - 2.0*r^(-6)
morse = let A = 4.0, r0 = 1.234
   @analytic r -> exp(-2.0*A*(r/r0-1.0)) - 2.0*exp(-A*(r/r0-1.0))
end

To create a "wrapped" AnalyticFunction{F64fun, ...} , use
```
lj_wrapped = F64fun(lj)
```

if an formula for an analytic function contains a sub-expression that
occurs multiple time, then this can be constructed as follows:
```
V = @analytic( r -> exp(s) * s, s = r^2 )
# is the same as
V = @analytic r -> exp(r^2) * r^2
```
"""
struct AnalyticFunction{F0,F1,F2} <: ExplicitPairPotential
   f::F0
   f_d::F1
   f_dd::F2
end

@pot AnalyticFunction

F64fun(p::AnalyticFunction) =
   AnalyticFunction(F64fun(p.f), F64fun(p.f_d), F64fun(p.f_dd))

cutoff(V::AnalyticFunction) = Inf

const WrappedAnalyticFunction = AnalyticFunction{F64fun,F64fun,F64fun}

"""
`@analytic`: generate C2 function from symbol
"""
macro analytic(args...)
   fexpr = substitute(args)
   quote
      AnalyticFunction(
        $(Base.FastMath.make_fastmath(esc(fexpr))),
        $(Base.FastMath.make_fastmath(esc(fdiff(fexpr, 1)))),
        $(Base.FastMath.make_fastmath(esc(fdiff(fexpr, 2))))
      )
   end
end


text:
import Dates: now
import SeisIO: InstrumentPosition,
  code2typ,
  fctoresp,
  getbandcode,
  note!,
  sμ,
  μs


text:
# Automatically generated using Clang.jl

using CEnum

# Skipping MacroDefinition: H3_EXPORT ( name ) name

const H3_VERSION_MAJOR = 3
const H3_VERSION_MINOR = 4
const H3_VERSION_PATCH = 2
const MAX_CELL_BNDRY_VERTS = 10
const H3Index = UInt64

struct GeoCoord
    lat::Cdouble
    lon::Cdouble
end

struct GeoBoundary
    numVerts::Cint
    verts::NTuple{10, GeoCoord}
end

struct Geofence
    numVerts::Cint
    verts::Ptr{GeoCoord}
end

struct GeoPolygon
    geofence::Geofence
    numHoles::Cint
    holes::Ptr{Geofence}
end

struct GeoMultiPolygon
    numPolygons::Cint
    polygons::Ptr{GeoPolygon}
end

struct LinkedGeoCoord
    vertex::GeoCoord
    next::Ptr{LinkedGeoCoord}
end

struct LinkedGeoLoop
    first::Ptr{LinkedGeoCoord}
    last::Ptr{LinkedGeoCoord}
    next::Ptr{LinkedGeoLoop}
end

struct LinkedGeoPolygon
    first::Ptr{LinkedGeoLoop}
    last::Ptr{LinkedGeoLoop}
    next::Ptr{LinkedGeoPolygon}
end

struct CoordIJ
    i::Cint
    j::Cint
end


text:
################################################################################
# QuickTest
#
# A Julia implementation of QuickCheck, a specification-based tester
#
# QuickCheck was originally written for Haskell by Koen Claessen and John Hughes
# http://www.cse.chalmers.se/~rjmh/QuickCheck/
#
################################################################################
module QuickTest

include("GenerateTestValue.jl")
include("TestProp.jl")

end


text:
using Revise
using TensorKit
using statmech_tm_solver
using Printf

Δ = -0.5
T = cmpo_xxz(Δ)

println("========================================================================")

psi = qbimps(rand, 4, 2, 3)
for chi in 8:8:80
    psi = expand(psi, chi, 1e-9)
    println("===========================  chi = $chi   ==============================")
    ix, δ = 0, 42.
    while ix < 250 && δ > 1e-12
        ΛA, psi = A_canonical(T, psi)
        ΛB, psi = B_canonical(T, psi)
        ix += 1
        δ = norm(ΛA - ΛB) / norm(ΛA + ΛB)
        println("$ix, ΛA: $ΛA, ΛB: $ΛB, δ: $δ")
        println("energy density: ", energy_quantum_xxz(psi.A, Δ))
    end

    sA = entanglement_spectrum(psi.A)
    sB = entanglement_spectrum(psi.B)

    io = open(@sprintf("result_qbimps_chi%d_xxz_Delta%.2f.txt", chi, Δ), "w")
    for s in sA print(io, "$s ") end
    println(io, " ")
    for s in sB print(io, "$s ") end
    println(io, " ")
    close(io)
end


text:
@i function i_fr_step!(r::AbstractVector{V3{T}}, v::AbstractVector{V3{T}}, m::AbstractVector; G=G_year_AU, dt) where T
    @routine @invcheckoff begin
        nplanets ← length(m)
        θ ← 1/(2-2^(1/3))
        @zeros typeof(dt) halfθdt θdt halfdt_minus_halfθdt dt_minus_twoθdt halfdt
        halfdt += dt / 2
        θdt += dt * θ
        halfθdt += θdt / 2
        halfdt_minus_halfθdt += halfdt - halfθdt
        dt_minus_twoθdt += dt
        dt_minus_twoθdt -= 2 * θdt
    end
    @inbounds @invcheckoff begin
        a ← zeros(V3{T}, nplanets)
        for j=1:nplanets
            r[j] += halfθdt*v[j]
        end
        @routine i_update_acceleration!(a, r, m, G)
        for j=1:nplanets
            v[j] += θdt*a[j]
        end
        ~@routine
        for j=1:nplanets
            r[j] += halfdt_minus_halfθdt*v[j]
        end
        @routine i_update_acceleration!(a, r, m, G)
        for j=1:nplanets
            v[j] += dt_minus_twoθdt*a[j]
        end
        ~@routine
        for j=1:nplanets
            r[j] += halfdt_minus_halfθdt*v[j]
        end
        @routine i_update_acceleration!(a, r, m, G)
        for j=1:nplanets
            v[j] += θdt*a[j]
        end
        ~@routine
        for j=1:nplanets
            r[j] += halfθdt*v[j]
        end
        a → zeros(V3{T}, nplanets)
    end

    ~@routine
end

@i function i_pefrl_step!(r::AbstractVector{V3{T}}, v::AbstractVector{V3{T}}, m::AbstractVector; G=G_year_AU, dt) where T
    @routine @invcheckoff begin
        nplanets ← length(m)
        ξ ← +0.1786178958448091E+00
        λ ← -0.2123418310626054E+00
        χ ← -0.6626458266981849E-01
        @zeros typeof(dt) ξh χh λh halfh_minus_λh h_minus_2χh_minus_2ξh halfh
        halfh += dt / 2
        ξh += ξ * dt
        λh += λ * dt
        χh += χ * dt
        halfh_minus_λh += halfh - λh
        h_minus_2χh_minus_2ξh += dt
        h_minus_2χh_minus_2ξh -= 2 * χh
        h_minus_2χh_minus_2ξh -= 2 * ξh
    end
    @invcheckoff @inbounds begin
        a ← zeros(V3{T}, nplanets)
        for j=1:nplanets
            r[j] += ξh*v[j]
        end
        @routine i_update_acceleration!(a, r, m, G)
        for j=1:nplanets
            v[j] += halfh_minus_λh*a[j]
        end
        ~@routine
        for j=1:nplanets
            r[j] += χh*v[j]
        end
        @routine i_update_acceleration!(a, r, m, G)
        for j=1:nplanets
            v[j] += λh*a[j]
        end
        ~@routine
        for j=1:nplanets
            r[j] += h_minus_2χh_minus_2ξh*v[j]
        end
        @routine i_update_acceleration!(a, r, m, G)
        for j=1:nplanets
            v[j] += λh*a[j]
        end
        ~@routine
        for j=1:nplanets
            r[j] += χh*v[j]
        end
        @routine i_update_acceleration!(a, r, m, G)
        for j=1:nplanets
            v[j] += halfh_minus_λh*a[j]
        end
        ~@routine
        for j=1:nplanets
            r[j] += ξh*v[j]
        end
        a → zeros(V3{T}, nplanets)
    end
    ~@routine
end

@i function i_pefrl!(r::AbstractVector{V3{T}}, v::AbstractVector{V3{T}}, m::AbstractVector; G=G_year_AU, n, dt) where T
    @invcheckoff for i=1:n
        i_pefrl_step!(r, v, m; dt=dt, G=G)
    end
end

@i function i_fr!(r::AbstractVector{V3{T}}, v::AbstractVector{V3{T}}, m::AbstractVector; G=G_year_AU, n, dt) where T
    @invcheckoff for i=1:n
        i_fr_step!(r, v, m; dt=dt, G=G)
    end
end


text:
@withkw @mixin immutable IF2Parameter <: IFParameter
    Ee::Float = 0mV
    Ei::Float = 0mV
end

@withkw @mixin type IF2 <: IF
    param::IF2Parameter = IF2Parameter()
end

@replace function integrate!(p::IF2, param::IF2Parameter, dt::Float)
    @inbounds for i = 1:N
        v[i] += dt * (ge[i] * (Ee - v[i]) + gi[i] * (Ei - v[i]) - (v[i] - El)) / τm
        ge[i] += dt * -ge[i] / τe
        gi[i] += dt * -gi[i] / τi
    end
    @inbounds for i = 1:N
        fire[i] = v[i] > Vt
        v[i] = ifelse(fire[i], Vr, v[i])
    end
end


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

function svdtest(n, iter)
    A = rand(n,n)
    B = svdfact(A)
    for i = 1:iter-1
        B = svdfact(A)
    end
    B
end

function schurtest(n, iter)
    A = rand(n,n)
    B = schurfact(A)
    for i = 1:iter-1
        B = schurfact(A)
    end
    B
end

function choleskytest(n, iter)
    A = rand(n,n)
    A = A'*A
    B = cholfact(A)
    for i = 1:iter-1
        B = cholfact(A)
    end
    B
end

function qrtest(n, iter)
    A = rand(n,n)
    B = qrfact(A)
    for i = 1:iter-1
        B = qrfact(A)
    end
    B
end

function lutest(n, iter)
    A = rand(n,n)
    B = lufact(A)
    for i = 1:iter-1
        B = lufact(A)
    end
    B
end

problemsizes = [(2, 20000, "tiny"), (2^4, 1000, "small"), (2^6, 100, "medium"), (2^8, 5, "large"), (2^10, 1, "huge")]
qr_problemsizes = [(2, 100000, "tiny"), (2^4, 5000, "small"), (2^6, 500, "medium"), (2^8, 5, "large"), (2^10, 1, "huge")]
lu_problemsizes = [(2, 100000, "tiny"), (2^4, 10000, "small"), (2^6, 1000, "medium"), (2^8, 100, "large"), (2^10, 2, "huge")]
testdata = [(choleskytest, "choleskytest", "Cholesky factorization", problemsizes),
            (schurtest, "schurtest", "Schur factorization", problemsizes),
            (svdtest, "svdtest", "Singular value decomposition", problemsizes),
            (qrtest, "qrtest", "QR factorization", qr_problemsizes),
            (lutest, "lutest", "LU factorization", lu_problemsizes)]
include("../perfgeneric.jl")


text:
# benchmark the motif finding algorithms
include("utils.jl")
include("search.jl")
include("data.jl")

"""
Print benchmarks for motifsearch_greedy
"""
function bench_greedy()

    count = 16
    len = 300
    motiflen = 25

    scores = []
    times = []

    for i = 1:100
        seqs, motif, inds = gendata(count, len, motiflen, 2)
        truescore = score(seqs, inds, motiflen)
        ret = @timed score(seqs, motifsearch_greedy(seqs, motiflen), motiflen) / truescore
        push!(scores, ret[1])
        push!(times, ret[2])
        if i % 10 == 0
            @printf("Step %d\tMean time per execution: %gs\n", i, mean(times))
        end
    end

    @printf("Mean score for motifsearch_greedy: %g\n", mean(scores))

end

max_stars_repo_name: dillondaudert/JuliaMotifs

max_issues_repo_name: dillondaudert/JuliaMotifs

max_forks_repo_name: dillondaudert/JuliaMotifs


text:
#####
##### Profiling MCTS
#####

using AlphaZero
using Base: @timed

# using Revise; Revise.includet("scripts/profile/mcts.jl")
# profile_mcts(nrep=64, nworkers=64, batched=true)

# Cost of simulation on one thread: 11μs
# Cost of simulation on all threads: 3μs
# To be compared to the cost of inference: 15μs to 30μs
# Also, by comparing `profile_mcts` with `batched=true` and `batched=false`, we conclude
# that the overhead of spawning tasks is not too big.

"""
Return the search time per inference request sent in μs/request.

!!! note
    For both the GPU and CPU to be maximally utilized, this number must be comparable to
    the result of `profile_inference`.
"""
function profile_mcts(
  exp::Experiment = Examples.experiments["connect-four"];
  batched=false,
  nrep=30,
  nworkers=1)

  oracle = MCTS.RandomOracle(exp.gspec)
  if batched
    spawn_oracle, done! =
      AlphaZero.batchify_oracles(
        oracle;
        num_workers=nworkers,
        batch_size=nworkers,
        fill_batches=false)
    oracles = [spawn_oracle() for _ in 1:nworkers]
  else
    oracles = [oracle for _ in 1:nworkers]
    done!() = nothing
  end

  players = [
    MctsPlayer(exp.gspec, oracle, exp.params.self_play.mcts)
    for oracle in oracles ]

  think(players[1], GI.init(exp.gspec)) # compile everything
  AlphaZero.reset_player!(players[1])

  info = @timed Threads.@threads for player in players
    for i in 1:nrep
      think(player, GI.init(exp.gspec))
    end
    done!()
  end
  traversed = sum(p.mcts.total_simulations for p in players)
  return info.time * 1_000_000 / traversed
end

# profile_mcts(nrep=64, nworkers=1)
# profile_mcts(nrep=64, nworkers=64)
# profile_mcts(nrep=64, nworkers=64, batched=true)


text:
module Collections

import Base: setindex!, done, get, haskey, isempty, length, next, getindex, start
import ..Order: Forward, Ordering, lt

export
    PriorityQueue,
    dequeue!,
    enqueue!,
    heapify!,
    heapify,
    heappop!,
    heappush!,
    isheap,
    peek

# Heap operations on flat arrays
# ------------------------------

# Binary heap indexing
heapleft(i::Integer) = 2i
heapright(i::Integer) = 2i + 1
heapparent(i::Integer) = div(i, 2)

# Binary min-heap percolate down.
function percolate_down!(xs::AbstractArray, i::Integer, x=xs[i], o::Ordering=Forward, len::Integer=length(xs))
    @inbounds while (l = heapleft(i)) <= len
        r = heapright(i)
        j = r > len || lt(o, xs[l], xs[r]) ? l : r
        if lt(o, xs[j], x)
            xs[i] = xs[j]
            i = j
        else
            break
        end
    end
    xs[i] = x
end

percolate_down!(xs::AbstractArray, i::Integer, o::Ordering, len::Integer=length(xs)) = percolate_down!(xs, i, xs[i], o, len)

# Binary min-heap percolate up.
function percolate_up!(xs::AbstractArray, i::Integer, x=xs[i], o::Ordering=Forward)
    @inbounds while (j = heapparent(i)) >= 1
        if lt(o, x, xs[j])
            xs[i] = xs[j]
            i = j
        else
            break
        end
    end
    xs[i] = x
end

percolate_up!{T}(xs::AbstractArray{T}, i::Integer, o::Ordering) = percolate_up!(xs, i, xs[i], o)

# Binary min-heap pop.
function heappop!(xs::AbstractArray, o::Ordering=Forward)
    x = xs[1]
    y = pop!(xs)
    if !isempty(xs)
        percolate_down!(xs, 1, y, o)
    end
    x
end

# Binary min-heap push.
function heappush!(xs::AbstractArray, x, o::Ordering=Forward)
    push!(xs, x)
    percolate_up!(xs, length(xs), x, o)
    xs
end

# Turn an arbitrary array into a binary min-heap in linear time.
function heapify!(xs::AbstractArray, o::Ordering=Forward)
    for i in heapparent(length(xs)):-1:1
        percolate_down!(xs, i, o)
    end
    xs
end

heapify(xs::AbstractArray, o::Ordering=Forward) = heapify!(copy(xs), o)

# Is an arbitrary array heap ordered?
function isheap(xs::AbstractArray, o::Ordering=Forward)
    for i in 1:div(length(xs), 2)
        if lt(o, xs[heapleft(i)], xs[i]) ||
           (heapright(i) <= length(xs) && lt(o, xs[heapright(i)], xs[i]))
            return false
        end
    end
    true
end

# PriorityQueue
# -------------

# A PriorityQueue that acts like a Dict, mapping values to their priorities,
# with the addition of a dequeue! function to remove the lowest priority
# element.
type PriorityQueue{K,V} <: Associative{K,V}
    # Binary heap of (element, priority) pairs.
    xs::Array{(K, V), 1}
    o::Ordering

    # Map elements to their index in xs
    index::Dict{K, Int}

    function PriorityQueue(o::Ordering)
        new(Array((K, V), 0), o, Dict{K, Int}())
    end

    PriorityQueue() = PriorityQueue{K,V}(Forward)

    function PriorityQueue(ks::AbstractArray{K}, vs::AbstractArray{V},
                           o::Ordering)
        # TODO: maybe deprecate
        if length(ks) != length(vs)
            error("key and value arrays must have equal lengths")
        end
        PriorityQueue{K,V}(zip(ks, vs))
    end

    function PriorityQueue(itr, o::Ordering)
        xs = Array((K, V), length(itr))
        index = Dict{K, Int}()
        for (i, (k, v)) in enumerate(itr)
            xs[i] = (k, v)
            if haskey(index, k)
                error("PriorityQueue keys must be unique")
            end
            index[k] = i
        end
        pq = new(xs, o, index)

        # heapify
        for i in heapparent(length(pq.xs)):-1:1
            percolate_down!(pq, i)
        end

        pq
    end
end

PriorityQueue(o::Ordering=Forward) = PriorityQueue{Any,Any}(o)

# TODO: maybe deprecate
PriorityQueue{K,V}(ks::AbstractArray{K}, vs::AbstractArray{V},
                   o::Ordering=Forward) = PriorityQueue{K,V}(ks, vs, o)

PriorityQueue{K,V}(kvs::Associative{K,V}, o::Ordering=Forward) = PriorityQueue{K,V}(kvs, o)

PriorityQueue{K,V}(a::AbstractArray{(K,V)}, o::Ordering=Forward) = PriorityQueue{K,V}(a, o)

length(pq::PriorityQueue) = length(pq.xs)
isempty(pq::PriorityQueue) = isempty(pq.xs)
haskey(pq::PriorityQueue, key) = haskey(pq.index, key)
peek(pq::PriorityQueue) = pq.xs[1]

function percolate_down!(pq::PriorityQueue, i::Integer)
    x = pq.xs[i]
    @inbounds while (l = heapleft(i)) <= length(pq)
        r = heapright(i)
        j = r > length(pq) || lt(pq.o, pq.xs[l][2], pq.xs[r][2]) ? l : r
        if lt(pq.o, pq.xs[j][2], x[2])
            pq.index[pq.xs[j][1]] = i
            pq.xs[i] = pq.xs[j]
            i = j
        else
            break
        end
    end
    pq.index[x[1]] = i
    pq.xs[i] = x
end

function percolate_up!(pq::PriorityQueue, i::Integer)
    x = pq.xs[i]
    @inbounds while i > 1
        j = heapparent(i)
        if lt(pq.o, x[2], pq.xs[j][2])
            pq.index[pq.xs[j][1]] = i
            pq.xs[i] = pq.xs[j]
            i = j
        else
            break
        end
    end
    pq.index[x[1]] = i
    pq.xs[i] = x
end

function getindex{K,V}(pq::PriorityQueue{K,V}, key)
    pq.xs[pq.index[key]][2]
end

function get{K,V}(pq::PriorityQueue{K,V}, key, deflt)
    i = get(pq.index, key, 0)
    i == 0 ? deflt : pq.xs[i][2]
end

# Change the priority of an existing element, or equeue it if it isn't present.
function setindex!{K,V}(pq::PriorityQueue{K, V}, value, key)
    if haskey(pq, key)
        i = pq.index[key]
        _, oldvalue = pq.xs[i]
        pq.xs[i] = (key, value)
        if lt(pq.o, oldvalue, value)
            percolate_down!(pq, i)
        else
            percolate_up!(pq, i)
        end
    else
        enqueue!(pq, key, value)
    end
end

function enqueue!{K,V}(pq::PriorityQueue{K,V}, key, value)
    if haskey(pq, key)
        error("PriorityQueue keys must be unique")
    end

    push!(pq.xs, (key, value))
    pq.index[key] = length(pq)
    percolate_up!(pq, length(pq))
    pq
end

function dequeue!(pq::PriorityQueue)
    x = pq.xs[1]
    y = pop!(pq.xs)
    if !isempty(pq)
        pq.xs[1] = y
        pq.index[pq.xs[1][1]] = 1
        percolate_down!(pq, 1)
    end
    delete!(pq.index, x[1])
    x[1]
end

function dequeue!(pq::PriorityQueue, key)
    idx = pop!(pq.index, key)  # throws key error if missing
    splice!(pq.xs, idx)
    for (k,v) in pq.index
        (v >= idx) && (pq.index[k] = (v-1))
    end
    key
end

# Unordered iteration through key value pairs in a PriorityQueue
start(pq::PriorityQueue) = start(pq.index)

done(pq::PriorityQueue, i) = done(pq.index, i)

function next(pq::PriorityQueue, i)
    (k, idx), i = next(pq.index, i)
    return ((k, pq.xs[idx][2]), i)
end

end # module Collections

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia