text:
using ChainRulesTestUtils
using FiniteDifferences
using ITensors
using Random

using ChainRulesCore: NoTangent

# TODO: maybe move this into ITensorChainRules module?
# These are useful definitions for testing code. The downside
# is adding dependencies on FiniteDifferences and ChainRulesTestUtils

#
# For ITensor compatibility with FiniteDifferences
#

function FiniteDifferences.to_vec(A::ITensor)
  # TODO: generalize to sparse tensors
  # TODO: define `itensor([1.0])` as well
  # as `itensor([1.0], ())` to help with generic code.
  function vec_to_ITensor(x)
    return isempty(inds(A)) ? ITensor(x[]) : itensor(x, inds(A))
  end
  return vec(array(A)), vec_to_ITensor
end

function FiniteDifferences.to_vec(x::Index)
  return (Bool[], _ -> x)
end

function FiniteDifferences.to_vec(x::Tuple{Vararg{Index}})
  return (Bool[], _ -> x)
end

function FiniteDifferences.to_vec(x::Pair{<:Tuple{Vararg{Index}},<:Tuple{Vararg{Index}}})
  return (Bool[], _ -> x)
end

function FiniteDifferences.rand_tangent(rng::AbstractRNG, A::ITensor)
  # TODO: generalize to sparse tensors
  return isempty(inds(A)) ? ITensor(randn(eltype(A))) : randomITensor(eltype(A), inds(A))
end

function FiniteDifferences.rand_tangent(rng::AbstractRNG, x::Index)
  return NoTangent()
end

function FiniteDifferences.rand_tangent(rng::AbstractRNG, x::Tuple{Vararg{Index}})
  return NoTangent()
end

function FiniteDifferences.rand_tangent(
  rng::AbstractRNG, x::Pair{<:Tuple{Vararg{Index}},<:Tuple{Vararg{Index}}}
)
  return NoTangent()
end

#
# For ITensor compatibility with ChainRulesTestUtils
#

function ChainRulesTestUtils.test_approx(
  actual::ITensor, expected::ITensor, msg=""; kwargs...
)
  ChainRulesTestUtils.@test_msg msg isapprox(actual, expected; kwargs...)
end

function ChainRulesTestUtils.test_approx(
  actual::ITensor, expected::Number, msg=""; kwargs...
)
  ChainRulesTestUtils.@test_msg msg isapprox(actual[], expected; kwargs...)
end

function ChainRulesTestUtils.test_approx(
  actual::Number, expected::ITensor, msg=""; kwargs...
)
  ChainRulesTestUtils.@test_msg msg isapprox(actual, expected[]; kwargs...)
end


text:
module Symplectic

  import Base.size
  import Base.getindex
  import Base.convert
  import Base.ndims
  import Base.size
  import Base.repeat

	import Base.*
	import Base./
	import Base.==

  # New symplectic matrix type
  type Symp
      x::Array{Bool, 2}
      z::Array{Bool, 2}
  end

  # Overload the basic matrix methods
  function size(s::Symp)
    return size(s.x)
  end

  function getindex(s::Symp, i, j)
    return Symp(convert(Array{Bool, 2}, s.x[i,j]), convert(Array{Bool, 2}, s.z[i,j]))
  end

  function convert(::Type{Array{Bool, 2}}, s::Bool)
    return fill(s, (1, 1))::Array{Bool, 2}
  end

  function ndims(s::Symp)
    return 2
  end

  function size(s::Symp, t::Int64)
    return size(s.x::Array{Bool, 2}, t::Int64)
  end

  function repeat(s::Symp; inner=[1, 1], outer=[1, 1])
    return Symp(repeat(s.x, inner=inner, outer=outer), repeat(s.z, inner=inner, outer=outer))
  end

  # Commutation
  function *(a::Symp, b::Symp)
    n_a = size(b)[1];
    n_b = size(a)[1];
    r_a = repeat(a, inner=[n_a, 1]);
    r_b = repeat(b, outer=[n_b, 1]);
    res = reducedim($, r_a.x & r_b.z, [2], false) $ reducedim($, r_a.z & r_b.x, [2], false);
    return reshape(res, (n_a, n_b))
  end

  # Anti-commutation
  function /(a::Symp, b::Symp)
    return ~(*(a::Symp, b::Symp))
  end

  # Equivalence
  function ==(a::Symp, b::Symp)
    eq_x = reducedim(&, a.x .== b.x, [2], true)
    eq_z = reducedim(&, a.z .== b.z, [2], true)
    return (eq_x & eq_z)[1]
  end

  export Symp, size, getindex, convert, ndims, size, repeat, *, /
end # module

ext: jl


text:
function verlet(pos::Float64, acc::Float64, dt::Float64)
    prev_pos = pos
    time = 0.0

    while (pos > 0)
        time += dt
        temp_pos = pos
        pos = pos * 2 - prev_pos + acc * dt * dt
        prev_pos = temp_pos
    end

    return time
end

function stormer_verlet(pos::Float64, acc::Float64, dt::Float64)
    prev_pos = pos
    time = 0.0
    vel = 0.0

    while (pos > 0.0)
        time += dt
        temp_pos = pos
        pos = pos * 2 - prev_pos + acc * dt * dt
        prev_pos = temp_pos

        # Because acceleration is constant, velocity is straightforward
        vel += acc * dt
    end

    return time, vel
end

function velocity_verlet(pos::Float64, acc::Float64, dt::Float64)
    prev_pos = pos
    time = 0.0
    vel = 0.0

    while (pos > 0.0)
        time += dt
        pos += vel * dt + 0.5 * acc * dt * dt;
        vel += acc * dt;
    end

    return time, vel
end

function main()
    time = verlet(5.0, -10.0, 0.01);
    println("Time for Verlet integration is: $(time)\n")

    time, vel = stormer_verlet(5.0, -10.0, 0.01);
    println("Time for Stormer Verlet integration is: $(time)")
    println("Velocity for Stormer Verlet integration is: $(vel)\n")

    time, vel = velocity_verlet(5.0, -10.0, 0.01);
    println("Time for velocity Verlet integration is: $(time)")
    println("Velocity for velocity Verlet integration is: $(vel)\n")
end

main()

max_stars_repo_name: atocil/algorithm-archive

max_issues_repo_name: atocil/algorithm-archive

max_forks_repo_name: atocil/algorithm-archive


text:
using CSV, DataFrames, Distributions, LinearAlgebra, Random

include("common.jl")

fp = open(result_path, "w")
println(fp, join(["experiment", "iter", "algorithm", "solve_time", "is_reward", "is_proportion_sold", "os_reward", "os_proportion_sold"], ","))

for hp_name in ORDERED_KEYS, iter in 1:num_instances
    hp_params = HYPERPARAMS[hp_name]
    dataset = OnlineAdLearning.normal_log_normal_data(num_features, num_samples, ρ=hp_params.ρ, σ=hp_params.σ, α=hp_params.α)
    samples = dataset.samples
    train_data      = OnlineAdLearning.DataSet(samples[1:num_train_points])
    validation_data = OnlineAdLearning.DataSet(samples[(num_train_points)+1:(num_train_points+num_validation_points)])
    test_data       = OnlineAdLearning.DataSet(samples[(num_train_points+num_validation_points+1):(num_train_points+num_validation_points+num_test_points)])
    @assert length(test_data.samples) == num_test_points

    for (algo_name, algo_obj) in ALGORITHMS
        learner = OnlineAdLearning.Learner(train_data, algo_obj; test_data=test_data, validation_data=validation_data)
        solve_time = @elapsed(results = OnlineAdLearning.train(learner))
        println(fp, join([hp_name,
                          iter,
                          algo_name,
                          solve_time,
                          results.final_is_eval.reward,
                          results.final_is_eval.proportion_sold,
                          results.final_os_eval.reward,
                          results.final_os_eval.proportion_sold], ","))
        flush(fp)
    end
end
close(fp)

max_stars_repo_name: joehuchette/reserve-price-optimization

max_issues_repo_name: joehuchette/reserve-price-optimization

max_forks_repo_name: joehuchette/reserve-price-optimization


text:
"""
    L0_reg(x, xbm, z, y, J, k, d, l)

Runs Iterative Hard Thresholding for GWAS data `x`, response `y`, and non-genetic
covariates `z` on a specific sparsity parameter `k`.

One needs to construct a SnpBitMatrix type (`xbm`) before running this function.

# Arguments:
+ `x`: A SnpArray, which can be memory mapped to a file. Does not engage in any linear algebra
+ `xbm`: The bitarray representation of `x`. This matrix is loaded in RAM and performs linear algebra. It's possible to set scale=false for xbm, especially when rare SNPs exist
+ `z`: Matrix of non-genetic covariates. The first column usually denotes the intercept.
+ `y`: Response vector
+ `J`: The number of maximum groups (set as 1 if no group infomation available)
+ `k`: Number of non-zero predictors in each group
+ `d`: A distribution (e.g. Normal, Poisson)
+ `l`: A link function (e.g. Loglink, ProbitLink)

# Optional Arguments:
+ `group` vector storing group membership
+ `weight` vector storing vector of weights containing prior knowledge on each SNP
+ `use_maf` indicates whether we want to scale the projection with minor allele frequencies (see paper)
+ `debias` is boolean indicating whether we debias at each iteration (see paper)
+ `show_info` boolean indicating whether we want to print results if model does not converge. Should set to false for multithread/multicore computing
+ `init` boolean indicating whether we want to initialize β to sensible values through fitting. This is not efficient yet.
+ `tol` is used to track convergence
+ `max_iter` is the maximum IHT iteration for a model to converge. Defaults to 200, or 100 for cross validation
+ `max_step` is the maximum number of backtracking. Since l0 norm is not convex, we have no ascent guarantee
"""
function L0_reg(
    x         :: SnpArray,
    xbm       :: SnpBitMatrix,
    z         :: AbstractMatrix{T},
    y         :: AbstractVector{T},
    J         :: Int,
    k         :: Int,
    d         :: UnivariateDistribution,
    l         :: Link;
    group     :: AbstractVector{Int} = Int[],
    weight    :: AbstractVector{T} = T[],
    use_maf   :: Bool = false,
    debias    :: Bool = false,
    show_info :: Bool = true,          # print things when model didn't converge
    init      :: Bool = false,         # not efficient. whether to initialize β to sensible values
    est_r     :: Int = 0,              # If 1, estimate r using MM algorithm. If 2, estimate r using Newton's Method
    tol       :: T = convert(T, 1e-4), # tolerance for tracking convergence
    max_iter  :: Int = 100,            # maximum IHT iterations
    max_step  :: Int = 5,              # maximum backtracking for each iteration
) where {T <: Float}

    #start timer
    start_time = time()

    # first handle errors
    @assert J >= 0        "Value of J (max number of groups) must be nonnegative!\n"
    @assert k >= 0        "Value of k (max predictors per group) must be nonnegative!\n"
    @assert max_iter >= 0 "Value of max_iter must be nonnegative!\n"
    @assert max_step >= 0 "Value of max_step must be nonnegative!\n"
    @assert tol > eps(T)  "Value of global tol must exceed machine precision!\n"
    checky(y, d) # make sure response data y is in the form compatible with specified GLM

    # initialize constants
    mm_iter     = 0                 # number of iterations
    tot_time    = 0.0               # compute time *within* L0_reg
    next_logl   = oftype(tol,-Inf)  # loglikelihood
    the_norm    = 0.0               # norm(b - b0)
    scaled_norm = 0.0               # the_norm / (norm(b0) + 1)
    η_step      = 0                 # counts number of backtracking steps for η
    converged   = false             # scaled_norm < tol?

    # Initialize variables.
    v = IHTVariables(x, z, y, J, k, group, weight)             # Placeholder variable for cleaner code
    full_grad = zeros(T, size(x, 2) + size(z, 2))              # Preallocated vector for efficiency
    init_iht_indices!(v, xbm, z, y, d, l, J, k)                # initialize non-zero indices
    copyto!(v.xk, @view(x[:, v.idx]), center=true, scale=true) # store relevant components of x for first iteration
    debias && (temp_glm = initialize_glm_object())             # Preallocated GLM variable for debiasing

    # If requested, compute initial guess for model b
    if init
        initialize_beta!(v, y, x, d, l)
        A_mul_B!(v.xb, v.zc, xbm, z, v.b, v.c)
    end

    # Begin 'iterative' hard thresholding algorithm
    for iter = 1:max_iter

        # notify and return current model if maximum iteration exceeded
        if iter >= max_iter
            mm_iter  = iter
            tot_time = time() - start_time
            show_info && printstyled("Did not converge after $max_iter iterations! The run time for IHT was " * string(tot_time) * " seconds and model size was" * string(k) * "\n", color=:red)
            break
        end

        # save values from previous iterate and update loglikelihood
        save_prev!(v)
        logl = next_logl

        # take one IHT step in positive score direction
        (η, η_step, next_logl, d) = iht_one_step!(v, x, xbm, z, y, J, k, d, l, logl, full_grad, iter, max_step, use_maf, est_r)

        # perform debiasing if requested
        if debias && sum(v.idx) == size(v.xk, 2)
            temp_glm = fit(GeneralizedLinearModel, v.xk, y, d, l)
            view(v.b, v.idx) .= temp_glm.pp.beta0
        end

        # track convergence
        the_norm    = max(chebyshev(v.b, v.b0), chebyshev(v.c, v.c0)) #max(abs(x - y))
        scaled_norm = the_norm / (max(norm(v.b0, Inf), norm(v.c0, Inf)) + 1.0)
        converged   = scaled_norm < tol
        if converged
        # if iter > 1 && abs(next_logl - logl) < tol * (abs(logl) + 1.0)
            tot_time = time() - start_time
            mm_iter  = iter
            break
        end
    end

    return ggIHTResults(tot_time, next_logl, mm_iter, v.b, v.c, J, k, v.group, d)
end #function L0_reg

"""
    L0_reg(x, z, y, J, k, d, l)

IHT algorithm that works on general matrices.

# Arguments
+ `x`: A general matrix. User should standardize it first.
+ `z`: Other covariates. The column of intercept should go here.
+ `y`: Response vector
+ `J`: The number of maximum groups (set as 1 if no group infomation available)
+ `k`: Number of non-zero predictors in each group
+ `d`: A distribution (e.g. Normal, Poisson)
+ `l`: A link function (e.g. Loglink, ProbitLink)

# Optional Arguments
Optional arguments are the same as the `L0_reg` that works on SnpArrays.
"""
function L0_reg(
    x         :: AbstractMatrix{T},
    z         :: AbstractMatrix{T},
    y         :: AbstractVector{T},
    J         :: Int,
    k         :: Int,
    d         :: UnivariateDistribution,
    l         :: Link;
    group     :: AbstractVector{Int} = Int[],
    weight    :: AbstractVector{T} = T[],
    use_maf   :: Bool = false,
    debias    :: Bool = false,
    show_info :: Bool = true,            # print things when model didn't converge
    init      :: Bool = false,           # not efficient. initializes β to sensible values
    est_r     :: Int  = 0,               # If 1, estimate r using MM algorithm. If 2, estimate r using Newton's Method
    tol       :: T    = convert(T, 1e-4),# tolerance for tracking convergence
    max_iter  :: Int  = 100,             # maximum IHT iterations
    max_step  :: Int  = 5,               # maximum backtracking for each iteration
) where {T <: Float}

    #start timer
    start_time = time()

    # first handle errors
    @assert J >= 0        "Value of J (max number of groups) must be nonnegative!\n"
    @assert k >= 0        "Value of k (max predictors per group) must be nonnegative!\n"
    @assert max_iter >= 0 "Value of max_iter must be nonnegative!\n"
    @assert max_step >= 0 "Value of max_step must be nonnegative!\n"
    @assert tol > eps(T)  "Value of global tol must exceed machine precision!\n"
    checky(y, d) # make sure response data y is in the form compatible with specified GLM

    # initialize constants
    mm_iter     = 0                 # number of iterations
    tot_time    = zero(T)           # compute time *within* L0_reg
    next_logl   = oftype(tol,-Inf)  # loglikelihood
    the_norm    = zero(T)           # norm(b - b0)
    scaled_norm = zero(T)           # the_norm / (norm(b0) + 1)
    η_step      = 0                 # counts number of backtracking steps for η
    converged   = false             # scaled_norm < tol?

    # Initialize variables.
    v = IHTVariables(x, z, y, J, k, group, weight) # Placeholder variable for cleaner code
    full_grad = zeros(T, size(x, 2) + size(z, 2))  # Preallocated vector for efficiency
    init_iht_indices!(v, x, z, y, d, l, J, k)      # initialize non-zero indices
    copyto!(v.xk, @view(x[:, v.idx]))              # store relevant components of x for first iteration
    debias && (temp_glm = initialize_glm_object()) # Preallocated GLM variable for debiasing

    # If requested, compute initial guess for model b
    if init
        initialize_beta!(v, y, x, d, l)
        A_mul_B!(v.xb, v.zc, x, z, v.b, v.c)
    end

    # Begin 'iterative' hard thresholding algorithm
    for iter = 1:max_iter

        # notify and return current model if maximum iteration exceeded
        if iter >= max_iter
            mm_iter  = iter
            tot_time = time() - start_time
            show_info && printstyled("Did not converge after $max_iter iterations! The run time for IHT was " * string(tot_time) * " seconds and model size was" * string(k) * "\n", color=:red)
            break
        end

        # save values from previous iterate and update loglikelihood
        save_prev!(v)
        logl = next_logl

        # take one IHT step in positive score direction
        (η, η_step, next_logl, d) = iht_one_step!(v, x, z, y, J, k, d, l, logl, full_grad, iter, max_step, use_maf, est_r)

        # perform debiasing if requested
        if debias && sum(v.idx) == size(v.xk, 2)
            temp_glm = fit(GeneralizedLinearModel, v.xk, y, d, l)
            view(v.b, v.idx) .= temp_glm.pp.beta0
        end

        # track convergence
        the_norm    = max(chebyshev(v.b, v.b0), chebyshev(v.c, v.c0)) #max(abs(x - y))
        scaled_norm = the_norm / (max(norm(v.b0, Inf), norm(v.c0, Inf)) + 1.0)
        converged   = scaled_norm < tol
        if converged
        # if iter > 1 && abs(next_logl - logl) < tol * (abs(logl) + 1.0)
            tot_time = time() - start_time
            mm_iter  = iter
            break
        end
    end

    return ggIHTResults(tot_time, next_logl, mm_iter, v.b, v.c, J, k, v.group, d)
end #function L0_reg

"""
This function performs 1 iteration of the IHT algorithm, backtracking a maximum of 5 times.
While IHT can strictly increase loglikelihood, we still allow it to potentially decrease
to avoid bad boundary cases.
"""
function iht_one_step!(v::IHTVariable{T}, x::SnpArray, xbm::SnpBitMatrix, z::AbstractMatrix{T},
    y::AbstractVector{T}, J::Int, k::Int, d::UnivariateDistribution, l::Link, old_logl::T,
    full_grad::AbstractVector{T}, iter::Int, nstep::Int, use_maf::Bool, est_r::Int) where {T <: Float}

    # first calculate step size
    η = iht_stepsize(v, z, d, l)

    # update b and c by taking gradient step v.b = P_k(β + ηv) where v is the score direction
    _iht_gradstep(v, η, J, k, full_grad)

    # update the linear predictors `xb` with the new proposed b, and use that to compute the mean
    update_xb!(v, x, z)
    update_μ!(v.μ, v.xb + v.zc, l)

    # if Negative Binomial, save previous r and update
    if typeof(d) == NegativeBinomial{Float64}
        old_r = d.r

        # Using MM algorithm
        if est_r == 1
            d = update_r!(d, y, v.μ)

        # Using Newton's method
        elseif est_r == 2
            new_r = mle_for_θ(y, v.μ, θ=d.r)
            d = NegativeBinomial(new_r, 0.5)
        end
    end

    # calculate current loglikelihood with the new computed xb and zc
    new_logl = loglikelihood(d, y, v.μ)

    η_step = 0
    while _iht_backtrack_(new_logl, old_logl, η_step, nstep)

        # stephalving
        η /= 2

        # if NegativeBinomial, reset r to previous r
        if typeof(d) == NegativeBinomial{Float64}
            d = NegativeBinomial(old_r, 0.5)
        end

        # recompute gradient step
        copyto!(v.b, v.b0)
        copyto!(v.c, v.c0)
        _iht_gradstep(v, η, J, k, full_grad)

        # recompute η = xb, μ = g(η), and loglikelihood to see if we're now increasing
        update_xb!(v, x, z)
        update_μ!(v.μ, v.xb + v.zc, l)

        # if Negative Binomial, update r
        if typeof(d) == NegativeBinomial{Float64}
            if est_r == 1
                d = update_r!(d, y, v.μ)
            elseif est_r == 2
                new_r = mle_for_θ(y, v.μ, θ=d.r)
                d = NegativeBinomial(new_r, 0.5)
            end
        end

        new_logl = loglikelihood(d, y, v.μ)

        # increment the counter
        η_step += 1
    end

    # compute score with the new mean μ
    score!(d, l, v, xbm, z, y)

    # check for finiteness before moving to the next iteration
    isnan(new_logl) && throw(error("Loglikelihood function is NaN, aborting..."))
    isinf(new_logl) && throw(error("Loglikelihood function is Inf, aborting..."))
    isinf(η) && throw(error("step size not finite! it is $η and max gradient is " * string(maximum(v.gk)) * "!!\n"))

    return η::T, η_step::Int, new_logl::T, d::UnivariateDistribution
end #function iht_one_step

"""
Performs 1 iteration of the IHT algorithm given a general matrix of floating point numbers.
"""
function iht_one_step!(v::IHTVariable{T}, x::AbstractMatrix, z::AbstractMatrix,
    y::AbstractVector{T}, J::Int, k::Int, d::UnivariateDistribution, l::Link, old_logl::T,
    full_grad::AbstractVector{T}, iter::Int, nstep::Int, use_maf::Bool, est_r::Int) where {T <: Float}

    # first calculate step size
    η = iht_stepsize(v, z, d, l)

    # update b and c by taking gradient step v.b = P_k(β + ηv) where v is the score direction
    _iht_gradstep(v, η, J, k, full_grad)

    # update the linear predictors `xb` with the new proposed b, and use that to compute the mean
    update_xb!(v, x, z)
    update_μ!(v.μ, v.xb + v.zc, l)

    # if Negative Binomial, save previous r and update
    if typeof(d) == NegativeBinomial{Float64}
        old_r = d.r

        # Using MM algorithm
        if est_r == 1
            d = update_r!(d, y, v.μ)

        # Using Newton's method
        elseif est_r == 2
            new_r = mle_for_θ(y, v.μ, θ=d.r)
            d = NegativeBinomial(new_r, 0.5)
        end
    end

    # calculate current loglikelihood with the new computed xb and zc
    new_logl = loglikelihood(d, y, v.μ)

    η_step = 0
    while _iht_backtrack_(new_logl, old_logl, η_step, nstep)

        # stephalving
        η /= 2

        # If NegativeBinomial, reset r to previous r
        if typeof(d) == NegativeBinomial{Float64}
            d = NegativeBinomial(old_r, 0.5)
        end

        # recompute gradient step
        copyto!(v.b, v.b0)
        copyto!(v.c, v.c0)
        _iht_gradstep(v, η, J, k, full_grad)

        # recompute η = xb, μ = g(η), and loglikelihood to see if we're now increasing
        update_xb!(v, x, z)
        update_μ!(v.μ, v.xb + v.zc, l)

        # If NegativeBinomial, update r
        if typeof(d) == NegativeBinomial{Float64}
            if est_r == 1
                d = update_r!(d, y, v.μ)
            elseif est_r == 2
                new_r = mle_for_θ(y, v.μ, θ=d.r)
                d = NegativeBinomial(new_r, 0.5)
            end
        end

        new_logl = loglikelihood(d, y, v.μ)

        # increment the counter
        η_step += 1
    end

    # compute score with the new mean μ
    score!(d, l, v, x, z, y)

    # check for finiteness before moving to the next iteration
    isnan(new_logl) && throw(error("Loglikelihood function is NaN, aborting..."))
    isinf(new_logl) && throw(error("Loglikelihood function is Inf, aborting..."))
    isinf(η) && throw(error("step size not finite! it is $η and max gradient is " * string(maximum(v.gk)) * "!!\n"))

    return η::T, η_step::Int, new_logl::T, d::UnivariateDistribution
end #function iht_one_step


text:
module SBTestBoundaries

using SpatialBoundaries
using Test

A = fill(1.0, 10, 12)
W = wombling(A)
@test all(W.m .== zero(eltype(W.m)))
@test all(W.θ .== zero(eltype(W.θ)))
B = boundaries(W)
@test length(B) == prod(size(A).-1)

A[4, 8] = 5.0
A[5, 8] = 5.0
W = wombling(A)
B = boundaries(W; threshold=0.03)
@test length(B) == 5

A = rand(Float64, 10, 10)
n = rand(eachindex(A), 10)
A[n] .= NaN
W = wombling(A)
B = boundaries(W)
@test !any(map(isnan, W.m[B]))

end


text:
# Array subindexing:
function Fgetrow(x::Array,inrow)
    row=Int(inrow[1])
    if row>0
       tmp=zeros(1,size(x,2))
        tmp[:]=x[row,:]
        return (tmp, zeros(size(x)))
    else
        return (zeros(size(x,2)), zeros(size(x)))
    end
end

function Fgetrow_inplace(value,auxvalue,x,inrow)
    row=Int(inrow[1])
    if row>0
        tmp=zeros(1,size(x,2))
        tmp[:]=x[row,:]
        copy!(value,tmp)
        copy!(auxvalue,zeros(size(x)))
    else
        return
        copy!(value,zeros(size(x,2)))
        copy!(auxvalue,zeros(size(x)))
    end
end

function Dgetrow(derivativeIDX,f_c,faux_c,grad_c,grad_n,x,inrow)
    row=Int(inrow[1]) # this is ugly since currently we force all values stored on the graph to be real arrays.
    if derivativeIDX==1
        if row>0
            faux_c[row,:]=grad_c
            axpy!(1.0,faux_c,grad_n)
        end
    end
end

#TODO: gpu version

export Fgetrow

Derivative[Fgetrow]=Dgetrow
Inplace[Fgetrow]=Fgetrow_inplace

import Base.*
#getindex(A::ADnode,:,i::Int)=ADnode(Fgetcol,[A ADint(i)])
getindex(A::ADnode,i::Int,::Colon)=ADnode(Fgetrow,[A ADint(i)])
getrow(A::ADnode,i)=ADnode(Fgetrow,[A ADint(i)])

export getindex
export getrow

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
#=
  Furniture moving (scheduling) in Julia ConstraintSolver.jl

  From Marriott & Stuckey: "Programming with constraints", page  112f

  One result:
     Sp Sc Sb  St
    [0, 0, 30, 45]

  Where the values are the start time for each task:
   Starts with piano time 0  : 3 persons  (30 min)
               chair time 0  : 1 person   (10 min)
               bed   time 30 : 3 persons  (15 min)
               table time 45 : 2 persons  (15 min)

   0       10   15    30      45      60

   piano --------------|bed---|
   piano --------------|       table----|
   piano --------------|bed---|
   chair --|            bed---|table----|

  There are many other solutions...

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP, Printf
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

#
# The cumulative constraint is in constraints_utils.jl
#
function furniture_moving1(print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>16,

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    # Furniture moving
    n = 4
    # [piano, chair, bed, table]
    durations = [30,10,15,15]
    resources = [3,1,3,2] # resource needed per task
    @variable(model, 0 <= start_times[1:n] <= 60, Int)
    @variable(model, minimum(durations) <= end_times[1:n]   <= 60, Int)
    @variable(model, 1 <= limit <= 3, Int)
    @variable(model, 0 <= max_time <= 60+maximum(durations),Int)

    for i in 1:n
        @constraint(model,end_times[i] == start_times[i] + durations[i])
    end
    @constraint(model, end_times .<= max_time)
    cumulative(model, start_times, durations, resources, limit)

    # This is very slow
    # @objective(model, Min, max_time)
    # @objective(model,Min,limit)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                start_timesx = convert.(Integer,JuMP.value.(start_times; result=sol))
                end_timesx = convert.(Integer,JuMP.value.(end_times; result=sol))
                max_timex = convert.(Integer,JuMP.value.(max_time; result=sol))
                limitx = convert.(Integer,JuMP.value.(limit; result=sol))
                println("start_times:$start_timesx")
                println("durations  :$durations")
                println("end_times  :$end_timesx")
                println("resources  :$resources")
                # println("max_time   :$max_timex")
                println("limit      :$limitx")

                for i in 1:n
                    println("Task $i: $(start_timesx[i])..$(end_timesx[i]) (duration:$(durations[i]) resource:$(resources[i])")
                end

                time_line = minimum(start_timesx):maximum(end_timesx)
                println("\nTime line: (tasks active (resource))")
                for t in time_line
                    @printf "time %3d: " t
                    total_resources = 0
                    for i in 1:length(start_timesx)
                        if start_timesx[i] <= t && t < start_timesx[i]+durations[i]
                            total_resources += resources[i]
                            @printf " %3d (%2d) " i resources[i]
                        else
                            @printf " %3s  %3s " "" ""
                        end
                    end
                    println("  total_resources: $total_resources")
                end
            end
        end
    else
        println("status:$status")
    end

    return status
end

print_solutions=true
all_solutions=false
@time furniture_moving1(print_solutions,all_solutions)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Statistics: mean

using FFTW
using GPUifyLoops

using Oceananigans.Operators

function ∇²_ppn!(grid::RegularCartesianGrid, f, ∇²f)
    @loop for k in (1:grid.Nz; blockIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds ∇²f[i, j, k] = ∇²_ppn(grid, f, i, j, k)
            end
        end
    end

    @synchronize
end

function mixed_fft_commutes(N)
    A = rand(N, N, N)
    Ã1 = FFTW.dct(FFTW.rfft(A, [1, 2]), 3)
    Ã2 = FFTW.rfft(FFTW.dct(A, 3), [1, 2])
    Ã1 ≈ Ã2
end

function mixed_ifft_commutes(N)
    A = rand(N, N, N)

    Ã1 = FFTW.dct(FFTW.rfft(A, [1, 2]), 3)
    Ã2 = FFTW.rfft(FFTW.dct(A, 3), [1, 2])

    A11 = FFTW.irfft(FFTW.idct(Ã1, 3), N, [1, 2])
    A12 = FFTW.idct(FFTW.irfft(Ã1, N, [1, 2]), 3)
    A21 = FFTW.irfft(FFTW.idct(Ã2, 3), N, [1, 2])
    A22 = FFTW.idct(FFTW.irfft(Ã2, N, [1, 2]), 3)
    A ≈ A11 && A ≈ A12 && A ≈ A21 && A ≈ A22
end

function fftw_planner_works(ft, Nx, Ny, Nz, planner_flag)
    g = RegularCartesianGrid(ft, (Nx, Ny, Nz), (100, 100, 100))
    RHS = CellField(Complex{ft}, CPU(), g)
    solver = PoissonSolver(g, RHS, FFTW.ESTIMATE)
    true  # Just making sure our PoissonSolver does not error/crash.
end

function poisson_ppn_planned_div_free_cpu(ft, Nx, Ny, Nz, planner_flag)
    g = RegularCartesianGrid(ft, (Nx, Ny, Nz), (100, 100, 100))

    RHS = CellField(Complex{ft}, CPU(), g)
    RHS_orig = CellField(Complex{ft}, CPU(), g)
    ϕ = CellField(Complex{ft}, CPU(), g)
    ∇²ϕ = CellField(Complex{ft}, CPU(), g)

    RHS.data .= rand(Nx, Ny, Nz)
    RHS.data .= RHS.data .- mean(RHS.data)

    RHS_orig.data .= copy(RHS.data)

    solver = PoissonSolver(g, RHS, planner_flag)

    solve_poisson_3d_ppn_planned!(solver, g, RHS, ϕ)
    ∇²_ppn!(g, ϕ, ∇²ϕ)

    ∇²ϕ.data ≈ RHS_orig.data
end


text:
import Term: make_logo

@testset "logo" begin
    @test_nothrow make_logo()

    logo = make_logo()
    @test logo.measure.w == 84
    @test logo.measure.h == 28
end


text:
using Plots, Dates; gr();

export qwckplot

##  qwckplot function
#
#examples:
#  GCMGridLoad()
#  qwckplot(MeshArrays.Depth)
#  qwckplot(MeshArrays.Depth,"Ocean Bottom Depth")

"""
    qwckplot(fld::gcmfaces)

Call qwckplot(fld::gcmfaces,ttl::String) with current date for title
"""
function qwckplot(fld::gcmfaces)
    tmp1=Dates.now()
    tmp1=Dates.format(tmp1, "yyyy-mm-dd HH:MM:SS")
    qwckplot(fld,"Plotted at time "*tmp1)
end

"""
    qwckplot(fld::gcmfaces,ttl::String)

Plot input using convert2array and heatmap + add title
"""
function qwckplot(fld::gcmfaces,ttl::String)
    arr=convert2array(fld)
    arr=permutedims(arr,[2 1])
    #This uses Plots.jl:
    p=heatmap(arr,title=ttl)
end


text:
using DynamicAxisWarping, Plots

ts = range(0, stop=4π, length=128)
x = LinearInterpolation(sin.(ts) .+ 0.1 .* randn.())
y = LinearInterpolation(sin.(1.1 .* ts))

t = range(0, stop = 1, length=100)
orig_plot = plot(t, x, label="x(t)")
plot!(t, y, label="y(t)")
xlabel!("Time t")
ylabel!("Amplitude")
title!("Original signals")

cost, ϕ, ψ = gdtw(x, y)
gdtw_plot = plot(t, x ∘ ϕ, label="x∘ϕ(t)")
plot!(t, y ∘ ψ, label="y∘ψ(t)")
xlabel!("Time t")
ylabel!("Amplitude")
title!("GDTW")

dtw_cost, t_dtw_x, t_dtw_y = dtw(x.x, y.x)
dtw_plot = plot(eachindex(t_dtw_x) ./ length(t_dtw_x), x.x[t_dtw_x], label="x(t_x)")
plot!(eachindex(t_dtw_x) ./ length(t_dtw_x), y.x[t_dtw_y], label="y(t_y)")
xlabel!("Time")
ylabel!("Amplitude")
title!("DTW")

warp_plot = plot(t, ϕ, label = "ϕ(t)")
plot!(t, ψ, label = "ψ(t)")
rescale_t_x = (t_dtw_x .- first(t_dtw_x)) ./ last(t_dtw_x)
rescale_t_y = (t_dtw_y .- first(t_dtw_y)) ./ last(t_dtw_y)
xlabel!("Time")
ylabel!("Rescaled time")
title!("Comparison of warpings")
plot!(t, LinearInterpolation(rescale_t_x), label="t_x")
plot!(t, LinearInterpolation(rescale_t_y), label="t_y")

plot(gdtw_plot, dtw_plot, orig_plot, warp_plot, layout = (2, 2), legend = :outertopright)


text:
############################################################
# joAbstractSAparallelLinearOperator #################################
############################################################

# includes
include("joAbstractSAparallelLinearOperator/constructors.jl")
include("joAbstractSAparallelLinearOperator/base_functions.jl")
include("joAbstractSAparallelLinearOperator/extra_functions.jl")


text:
@testset "linear interpolation" begin
    for i in 1:10
        len = rand(10:80)
        x = Array{Union{Float64, Missing}}(rand(1:100.0,len))
        x[rand(1:len-5)] = missing
        @test length(linear_interpolation(x)) == len
    end
end


text:
module MKLSparse

function __init__()
    ccall((:MKL_Set_Interface_Layer, :libmkl_rt), Cint, (Cint,), Base.USE_BLAS64 ? 1 : 0)
end

include(joinpath("BLAS", "BLAS.jl"))

end # module


text:
# These are a workaround for splatting penalty on 0.4.5

for i = 1:32
    sig = [Symbol("_x_$j") for j = 1:i]

    expr =  Expr(:stagedfunction, Expr(:call, :(Base.call{Names}), :(::Type{Row{Names}}), sig...),
        Expr(:block,
            :(if !isa(Names, Tuple)
                str = "Row parameter 1 (Names) is expected to be a tuple of unique symbols, got $Names" # TODO: reinsert colons in symbols?
                return :(error($str))
            end),
            :(newtype = Row{Names, $(Expr(:curly, :Tuple, sig...))}),
            Expr(:return, Expr(:quote, Expr(:block, Expr(:meta, :inline), Expr(:call, Expr(:$,:newtype), sig...))))))

    eval(expr)

    expr =  Expr(:stagedfunction, Expr(:call, :(Base.call{Names}), :(::Type{Table{Names}}), sig...),
        Expr(:block,
            :(if !isa(Names, Tuple)
                str = "Table parameter 1 (Names) is expected to be a tuple of unique symbols, got $Names" # TODO: reinsert colons in symbols?
                return :(error($str))
            end),
            :(newtype = Table{Names, $(Expr(:curly, :Tuple, sig...))}),
            Expr(:return, Expr(:quote, Expr(:block, Expr(:meta, :inline), Expr(:call, Expr(:$,:newtype), sig...))))))

    eval(expr)
end


text:
# ---------------------------------------------------------------------------- #
#
#   mesh2D.jl
#
#   Type for 2D meshes
#   Inherits from the abstract "mesh" type
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    Mesh2D

Mesh2D type:
Type for 2D meshes consisting of triangles. It holds the nodes and connectivity
information.
"""
struct Mesh2D <: Mesh

  dim::Int64              # Dimension of the problem

  porder::Int64           # Polynomial order

  n::Int64                # Number of nodes

  p::Matrix{Float64}      # Nodal locations
  ploc::Matrix{Float64}   # Local nodal locations
  tloc::Matrix{Int64}     # Local triangles
  t::Matrix{Int64}        # Triangle - node connectivity
  t2f::Matrix{Int64}      # Triangle - face connectivity
  f::Matrix{Int64}        # Face - node/triangle connectivity
  fb::Matrix{Int64}       # Boundary face info
  nodes::Array{Float64,3} # Nodes on which solution is evaluated

  boundtags::Vector{String} # Boundary tags

end

"""
    Mesh2D( name::String, porder_::Int64; N = 5::Int64 )

Constructor for one of the default meshes. Currently "square" is implemented.
"""
function Mesh2D( name::String, porder_::Int64; N = 5, M = N, scl = 1. )

  if name == "square"
    (p_, t_, bel_, tags_) = makesquare( N, M )
  elseif name == "single"
     p_   = [ 0.0 0.0;
              1.0 0.0;
              0.0 1.0 ]
     t_   = [ 1 2 3 ]
     bel_ = [ 2 3 1;
              3 1 2;
              1 2 3 ]
     tags_ = [ "diag", "vert", "hort" ]
  elseif name[end-3:end] == ".su2"
    (p_, t_, bel_, tags_) = readSU2_2D( name )
  elseif name[end-3:end] == ".msh" # BAMG
    (p_, t_, bel_, tags_ ) = readBAMG( name )
  elseif name[end-4:end] == ".mesh" # FEFLOA
    (p_, t_, bel_, tags_ ) = readFEFLOA_2D( name )
  else
    error("Mesh2D: Unknown mesh type")
  end

  p_ .*= scl

  f_, t2f_, nodes_, ploc_, tloc_, fb_ = genmesh( porder_, p_, t_, bel_ )

  n_ = size( p_, 1 )

  Mesh2D( 2, porder_, n_, p_, ploc_, tloc_, t_, t2f_, f_, fb_, nodes_, tags_ )

end

"""
    Mesh2D( porder_::Int64, p::Matrix{Float64}, t::Matrix{Int64}, bel::Matrix{Int64} )

Constructor for nodes locations (`p`), triangle connectivity (`t`), and
boundary element information (`bel`) given.
"""
function Mesh2D( porder_::Int64, p_::Matrix{Float64}, t_::Matrix{Int64}, bel_::Matrix{Int64} )

  (f_, t2f_, nodes_, ploc_, tloc_, fb_) = genmesh( porder_, p_, t_, bel_ )

  n_ = size( p_, 1 )

  Mesh2D( 2, porder_, n_, p_, ploc_, tloc_, t_, t2f_, f_, fb_, nodes_)

end

"""
    genmesh( porder::Int64, p::Matrix{Float64}, t::Matrix{Int64}, bel::Matrix{Int64} )

Generates all necessary connectivity information given nodes locations (`p`),
triangle connectivity (`t`), and boundary element information (`bel`).
"""
function genmesh( porder::Int64, p::Matrix{Float64}, t::Matrix{Int64}, bel_::Matrix{Int64} )

  (f_, t2f_, fb_)           = genFaces2D( t, bel_ )
  nodes_, ploc_, tloc_ = genNodes2D( porder, p, t )

  return f_, t2f_, nodes_, ploc_, tloc_, fb_

end

"""
    genFaces2D( t::Matrix{Int64}, bel::Matrix{Int64} )

Generates face connectivity given triangle connectivity (`t`), and
boundary element information (`bel`).
"""
function genFaces2D( t::Matrix{Int64}, bel::Matrix{Int64} )

  nt = size(t,1)

  edges  = vcat( t[:,[2,3]], t[:,[3,1]], t[:,[1,2]] ) # This holds all edges (but multiple copies)
  triang = 1:nt
  triang = vcat( triang, triang, triang )

  ne = size(edges,1)

  tflip = Array{Int}( ne, 1 )

  for ii = 1:ne
    tflip[ii] = 1
    if edges[ii,1] > edges[ii,2]
      tflip[ii] = -1
      temp = edges[ii,1]
      edges[ii,1] = edges[ii,2]
      edges[ii,2] = temp
    end
  end

  # This index links to the first unique index in the array edges
  ix  = groupslices( edges, 1 )

  t2f = fill(0::Int64, nt, 3)

  # Find unique vector index
  # Need to find an index that links from edges to the unique edges (without gaps)
  jx = Array{Int64}( size(ix) )
  jx[1] = 0
  lind = 0 # maximum unique index
  mind = 0 # maximum index in edges that had unique index
  ne = 0  # number of unique elements

  indUni = fill( false, size(ix,1) ) # index of unique elements

  for ii = 1:length(ix)
    temp    = ix[ii] - lind - 1
    temp2   = ix[ii] - mind - 1
    jx[ii]  = temp
    mind    = max(mind,ix[ii])

    if temp2 >= 0 # Found a unique index
      lind = lind + 1
      ne  = ne + 1
      indUni[ii] = true
    end

    if temp2 < 0 # Need to link to the unique index
      jx[ii] = jx[ix[ii]]
    end

  end

  f = fill( 0::Int64, ne, 4 )

  f[:,1:2] = edges[indUni,:]

  # make t2f and complete f
  for jj = 1:length(tflip)

    ind  = triang[jj]
    indt = Int64(ceil(jj / nt))

    indf = ix[jj] - jx[jj]

    if tflip[jj] == 1
      f[indf,3] = ind
    else
      f[indf,4] = ind
    end

    t2f[ ind, indt ] = tflip[jj] * indf

  end

  ### Boundary Elements
  # include information from 'bel' into 'f': the index of the boundaries

  bel[:,1:2] = sort(bel[:,1:2],2) # sort such that we compare against f
  fb = fill( 0, size(bel,1), 2 )
  nfb = 1
  # Look in third array
  indb = find( f[:,3] .== 0 )
  for jj = 1:length(indb)
    indf = indb[jj]

    # Find corresponding element in boundary element vector
    indel = find( all( bel[:,1:2] .== f[indf,1:2]', 2) )

    # Ensure boundary element is oriented counter-clockwise
    #   flip orientation of face in t2f
    indt = f[indf,4]
    indfint = find( abs.(t2f[indt,:]) .== indf )
    t2f[ f[indf,4], indfint ] = - t2f[ f[indf,4], indfint ]
    #   flip nodes in f
    temp = f[indf,1]
    f[indf,1] = f[indf,2]
    f[indf,2] = temp
    f[indf,3] = f[indf,4]

    # Mark correct boundary
    f[indf,4] = - bel[indel[1],3]
    fb[nfb,1] = indf
    fb[nfb,2] = bel[indel[1],3]
    nfb += 1
  end

  # Look in fourth array
  indb = find( f[:,4] .== 0 )
  for jj = 1:length(indb)
    indf = indb[jj]

    # Find corresponding element in boundary element vector
    indel = find( all( bel[:,1:2] .== f[indf,1:2]', 2) )

    # Mark correct boundary
    f[indf,4] = - bel[indel[1],3]
    fb[nfb,1] = indf
    fb[nfb,2] = bel[indel[1],3]
    nfb += 1
  end

  return f, t2f, fb

end

"""
    genNodes2D( porder::Int64, p::Matrix{Float64}, t::Matrix{Int64} )

Generates the higher order nodes inside each element.
"""
function genNodes2D( porder::Int64, p::Matrix{Float64}, t::Matrix{Int64} )

  nt     = size(t,1)
  (plocal,tlocal) = genlocal( porder )

  npl    = size(plocal,1)
  nodes  = fill( 0.0, npl, 2, nt )

  # We use barycentric coordinates for this
  # for interpolation with barycentric coordinates at the point (xi,yi) with
  # barycentric coordinates (li1, li2, li3), the following holds:
  #   f(xi,yi) = li1 * f(x1,y1) + li2 * f(x2,y2) + li3 * f(x3,y3)

  for ii in 1:nt, jj in 1:npl
    for kk in 1:3
      nodes[jj,1,ii] += plocal[jj,kk] * p[ t[ii,kk], 1 ]
      nodes[jj,2,ii] += plocal[jj,kk] * p[ t[ii,kk], 2 ]
    end
  end

  return nodes, plocal, tlocal

end

"""
    genlocal( porder::Int64 )

Generates the barycentric coordinates for the higher order nodes inside
the master element (triangle).
"""
function genlocal( porder::Int64 )

  # Generate the barycentric coordinates
  # here these are implemented as (1 - x - y, x, y)

  n = Int64( round( 0.5 * (porder+1) * (porder+2) ) )

  plocal = Array{Float64}( n, 2 )
  kk = 1
  for jj = 1:porder+1
    for ii = 1:(porder + 2 - jj)
      plocal[kk,:] = [ (ii-1)/porder, (jj-1)/porder ]
      kk = kk + 1
    end
  end

  temp = 1 - sum(plocal,2)

  plocal = hcat( temp, plocal )

  # Switch to have corners first, then edges
  if porder > 1
    plocal2 = fill( 0.0, size(plocal) )

    #   corners
    plocal2[ 1:3, :] = plocal[ [1 porder+1 n], :]
    #   long edge
    kk = 2*porder + 1
    for jj in 1:(porder-1)
      plocal2[ (3+jj), :] = plocal[ kk, :]
      kk += porder - jj
    end
    #   west edge
    kk = n - 2
    for jj in 1:(porder-1)
      plocal2[ (2+porder+jj), :] = plocal[ kk, :]
      kk -= 2 + jj
    end
    #   south edge
    kk = 2
    for jj in 1:(porder-1)
      plocal2[ (1+2*porder+jj), :] = plocal[ kk, :]
      kk += 1
    end
    #   center
    kk = porder + 3
    pp = 3*(porder-1) + 4
    for jj in 1:(porder-2)
      for ii in 1:(porder-1-jj)
        plocal2[ pp, :] = plocal[ kk, :]
        pp += 1
        kk += 1
      end
      kk += 2
    end
  else
    plocal2 = plocal
  end

  # generate local triangles -- used for plotting
  tlocal = Matrix{Int64}( porder^2, 3 )
  if porder == 1
    tlocal[1,:]  = [1 2 3]
  elseif porder == 2
    tlocal[1,:]  = [1 2 4]
    tlocal[2,:]  = [2 5 4]
    tlocal[3,:]  = [2 3 5]
    tlocal[4,:]  = [4 5 6]
  elseif porder == 3
    tlocal[1,:]  = [1 8 7]
    tlocal[2,:]  = [8 10 7]
    tlocal[3,:]  = [8 9 10]
    tlocal[4,:]  = [9 4 10]
    tlocal[5,:]  = [9 2 4]
    tlocal[6,:]  = [7 10 6]
    tlocal[7,:]  = [10 5 6]
    tlocal[8,:]  = [10 4 5]
    tlocal[9,:]  = [6 5 3]
  elseif porder == 4
    tlocal[1,:]  = [1 10 9]
    tlocal[2,:]  = [10 13 9]
    tlocal[3,:]  = [10 11 13]
    tlocal[4,:]  = [11 14 13]
    tlocal[5,:]  = [11 12 14]
    tlocal[6,:]  = [12 4 14]
    tlocal[7,:]  = [12 2 4]
    tlocal[8,:]  = [9 13 8]
    tlocal[9,:]  = [13 15 8]
    tlocal[10,:] = [13 14 15]
    tlocal[11,:] = [14 5 15]
    tlocal[12,:] = [14 4 5]
    tlocal[13,:] = [8 15 7]
    tlocal[14,:] = [15 6 7]
    tlocal[15,:] = [15 5 6]
    tlocal[16,:] = [7 6 3]
  end

  return (plocal2, tlocal)

end

"""
    makesquare( n::Int64, m::Int64 )

Generates the nodes locations, triangle connectivity and boundary element
information for a square [0,1] x [0,1].
"""
function makesquare( n::Int64, m::Int64 )

  # boundary 1 is South
  # boundary 2 is East
  # boundary 3 is North
  # boundary 4 is West

  p   = Array{Float64}( n*m, 2 )
  t   = Array{Int64}( 2*(n-1)*(m-1), 3 )
  bel = Array{Int64}( 2*(n-1) + 2*(m-1), 3 )

  # nodes
  for ii = 1:n, jj = 1:m

    p[ii + (jj-1) * n, :] = [ 1/(n-1) * (ii - 1), 1/(m-1) * (jj - 1) ]

  end

  # triangles
  for ii = 1:(n-1), jj = 1:(m-1)

    # Lower triangles
    t[2*(ii + (jj-1) * (n-1)) - 1, :] = [ ii + (jj-1) * n,
              ii + 1 + (jj-1) * n,
              ii + 1 + jj     * n ]

    # Upper triangles
    t[2*(ii + (jj-1) * (n-1)),     :] = [ ii + (jj-1) * n,
              ii + 1 + jj     * n,
              ii     + jj     * n ]

  end

  # boundary elements
  kk = 1
  for ii in 1:n-1 # South
    bel[kk,:] = [ ii, ii + 1, 1 ]
    kk = kk + 1
  end
  for jj in 1:m-1 # East
    ii = n - 1
    bel[kk,:] = [ ii + 1 + (jj-1) * n, ii + 1 + jj * n, 2 ]
    kk = kk + 1
  end
  for ii in n-1:-1:1 # North
    jj = m - 1
    bel[kk,:] = [ ii + 1 + jj * n, ii + jj * n, 3 ]
    kk = kk + 1
  end
  for jj in m-1:-1:1 # East
    ii = 1
    bel[kk,:] = [ ii + jj * n, ii + (jj-1) * n, 4 ]
    kk = kk + 1
  end

  tags = [ "south", "east", "north", "west" ]

  return p, t, bel, tags

end


text:
# This file is a part of ElasticArrays.jl, licensed under the MIT License (MIT).

__precompile__(true)

module ElasticArrays

using Compat
using Compat.Markdown

include("util.jl")
include("elasticarray.jl")

end # module


text:
function BiasDict13(; b1, bη, bK², bδη, bη², bKKpara, btd, bδΠ2para, bηΠ2para, bΠ2Kpara, bΠ2para, bΠ3para)
    BiasOperatorDict13 = Dict([("δ3", b1), ("η3", bη) , ("2tr[KK2]", bK²), ("δη2",bδη),
                           ("2ηη2", bη²), ("2KK2∥", bKKpara), ("Otd", btd), ("δΠ2∥", bδΠ2para),
                           ("ηΠ2∥", bηΠ2para), ("Π2K∥", bΠ2Kpara), ("s^k∂kΠ2∥", -bΠ2para),
                           ("u2∥∂∥δ", -b1), ("u2∥∂∥η", -bη), ("u∥∂∥η2", -bη),
                           ("u∥∂∥Π2∥", -bΠ2para), ("Π3∥", bΠ3para+2bΠ2para)])
end

function BiasDictlhd(; b1, bη, P0_ϵ, b∇2δ, β∇2v, β∂2parav, P2_ϵ, P2_ϵε_η)
    BiasOperatorDictlhd = Dict([("b1", b1),("bη", bη), ("P0ϵ", P0_ϵ), ("b∇2δ", b∇2δ), ("β∇2v", β∇2v),
                            ("β∂2∥v", β∂2parav), ("P2ϵ", P2_ϵ), ("P2ϵεη", P2_ϵε_η)])
end

function BiasDict22(; b1, bη, b2, bK², bδη, bη², bKKpara, bΠ2para)
    BiasOperatorDict22 = Dict([("b1",b1), ("bη", bη), ("b2", b2), ("bK2", bK²), ("bδη", bδη),
                           ("bη2", bη²), ("bKK∥", bKKpara), ("bΠ2∥", bΠ2para)])
end

max_stars_repo_name: JosephTomlinson/GeneralBiasPk

max_issues_repo_name: JosephTomlinson/GeneralBiasPk

max_forks_repo_name: JosephTomlinson/GeneralBiasPk


text:
"""
    XfoilGlobals

Holds XFOIL global variables
"""
struct XfoilGlobals{T}
    #cr01
    version::Array{T,1}
    #cr03
    aij::Array{T,2}
    dij::Array{T,2}
    #cr04
    qinv::Array{T,1}
    qvis::Array{T,1}
    cpi::Array{T,1}
    cpv::Array{T,1}
    qinvu::Array{T,2}
    qinv_a::Array{T,1}
    #cr05
    x::Array{T,1}
    y::Array{T,1}
    xp::Array{T,1}
    yp::Array{T,1}
    s::Array{T,1}
    sle::Array{T,1}
    xle::Array{T,1}
    yle::Array{T,1}
    xte::Array{T,1}
    yte::Array{T,1}
    chord::Array{T,1}
    yimage::Array{T,1}
    wgap::Array{T,1}
    waklen::Array{T,1}
    #cr06
    gam::Array{T,1}
    gamu::Array{T,2}
    gam_a::Array{T,1}
    sig::Array{T,1}
    nx::Array{T,1}
    ny::Array{T,1}
    apanel::Array{T,1}
    sst::Array{T,1}
    sst_go::Array{T,1}
    sst_gp::Array{T,1}
    gamte::Array{T,1}
    gamte_a::Array{T,1}
    sigte::Array{T,1}
    sigte_a::Array{T,1}
    dste::Array{T,1}
    ante::Array{T,1}
    aste::Array{T,1}
    #cr07
    ssple::Array{T,1}
    sspec::Array{T,1}
    xspoc::Array{T,1}
    yspoc::Array{T,1}
    qgamm::Array{T,1}
    qspec::Array{T,2}
    qspecp::Array{T,2}
    algam::Array{T,1}
    clgam::Array{T,1}
    cmgam::Array{T,1}
    alqsp::Array{T,1}
    clqsp::Array{T,1}
    cmqsp::Array{T,1}
    qf0::Array{T,1}
    qf1::Array{T,1}
    qf2::Array{T,1}
    qf3::Array{T,1}
    qdof0::Array{T,1}
    qdof1::Array{T,1}
    qdof2::Array{T,1}
    qdof3::Array{T,1}
    clspec::Array{T,1}
    ffilt::Array{T,1}
    #cr09
    adeg::Array{T,1}
    alfa::Array{T,1}
    awake::Array{T,1}
    mvisc::Array{T,1}
    avisc::Array{T,1}
    xcmref::Array{T,1}
    ycmref::Array{T,1}
    cl::Array{T,1}
    cm::Array{T,1}
    cd::Array{T,1}
    cdp::Array{T,1}
    cdf::Array{T,1}
    cl_alf::Array{T,1}
    cl_msq::Array{T,1}
    psio::Array{T,1}
    circ::Array{T,1}
    cosa::Array{T,1}
    sina::Array{T,1}
    qinf::Array{T,1}
    gamma::Array{T,1}
    gamm1::Array{T,1}
    minf1::Array{T,1}
    minf::Array{T,1}
    minf_cl::Array{T,1}
    tklam::Array{T,1}
    tkl_msq::Array{T,1}
    cpstar::Array{T,1}
    qstar::Array{T,1}
    cpmn::Array{T,1}
    cpmni::Array{T,1}
    cpmnv::Array{T,1}
    xcpmni::Array{T,1}
    xcpmnv::Array{T,1}
    #cr10
    xpref::Array{T,1}
    cpref::Array{T,1}
    verspol::Array{T,1}
    cpolxy::Array{T,3}
    machp1::Array{T,1}
    reynp1::Array{T,1}
    acritp::Array{T,1}
    xstripp::Array{T,2}
    #cr11
    pi::Array{T,1}
    hopi::Array{T,1}
    qopi::Array{T,1}
    dtor::Array{T,1}
    #cr12
    cvpar::Array{T,1}
    cterat::Array{T,1}
    ctrrat::Array{T,1}
    xsref1::Array{T,1}
    xsref2::Array{T,1}
    xpref1::Array{T,1}
    xpref2::Array{T,1}
    #cr13
    size::Array{T,1}
    scrnfr::Array{T,1}
    plotar::Array{T,1}
    pfac::Array{T,1}
    qfac::Array{T,1}
    vfac::Array{T,1}
    xwind::Array{T,1}
    ywind::Array{T,1}
    xpage::Array{T,1}
    ypage::Array{T,1}
    xmarg::Array{T,1}
    ymarg::Array{T,1}
    ch::Array{T,1}
    chg::Array{T,1}
    chq::Array{T,1}
    xofair::Array{T,1}
    yofair::Array{T,1}
    facair::Array{T,1}
    xofa::Array{T,1}
    yofa::Array{T,1}
    faca::Array{T,1}
    uprwt::Array{T,1}
    cpmin::Array{T,1}
    cpmax::Array{T,1}
    cpdel::Array{T,1}
    cpolplf::Array{T,2}
    xcdwid::Array{T,1}
    xalwid::Array{T,1}
    xocwid::Array{T,1}
    #cr14
    xb::Array{T,1}
    yb::Array{T,1}
    xbp::Array{T,1}
    ybp::Array{T,1}
    sb::Array{T,1}
    snew::Array{T,1}
    xbf::Array{T,1}
    ybf::Array{T,1}
    xof::Array{T,1}
    yof::Array{T,1}
    hmom::Array{T,1}
    hfx::Array{T,1}
    hfy::Array{T,1}
    xbmin::Array{T,1}
    xbmax::Array{T,1}
    ybmin::Array{T,1}
    ybmax::Array{T,1}
    sble::Array{T,1}
    chordb::Array{T,1}
    areab::Array{T,1}
    radble::Array{T,1}
    angbte::Array{T,1}
    ei11ba::Array{T,1}
    ei22ba::Array{T,1}
    apx1ba::Array{T,1}
    apx2ba::Array{T,1}
    ei11bt::Array{T,1}
    ei22bt::Array{T,1}
    apx1bt::Array{T,1}
    apx2bt::Array{T,1}
    thickb::Array{T,1}
    cambrb::Array{T,1}
    xcm::Array{T,1}
    ycm::Array{T,1}
    scm::Array{T,1}
    xcmp::Array{T,1}
    ycmp::Array{T,1}
    xtk::Array{T,1}
    ytk::Array{T,1}
    stk::Array{T,1}
    xtkp::Array{T,1}
    ytkp::Array{T,1}
    #cr15
    xssi::Array{T,2}
    uedg::Array{T,2}
    uinv::Array{T,2}
    mass::Array{T,2}
    thet::Array{T,2}
    dstr::Array{T,2}
    ctau::Array{T,2}
    delt::Array{T,2}
    uslp::Array{T,2}
    guxq::Array{T,2}
    guxd::Array{T,2}
    tau::Array{T,2}
    dis::Array{T,2}
    ctq::Array{T,2}
    vti::Array{T,2}
    reinf1::Array{T,1}
    reinf::Array{T,1}
    reinf_cl::Array{T,1}
    acrit::Array{T,1}
    xstrip::Array{T,1}
    xoctr::Array{T,1}
    yoctr::Array{T,1}
    xssitr::Array{T,1}
    uinv_a::Array{T,2}
    #cr17
    rmsbl::Array{T,1}
    rmxbl::Array{T,1}
    rlx::Array{T,1}
    vaccel::Array{T,1}
    #cr18
    xsf::Array{T,1}
    ysf::Array{T,1}
    xoff::Array{T,1}
    yoff::Array{T,1}
    xgmin::Array{T,1}
    xgmax::Array{T,1}
    ygmin::Array{T,1}
    ygmax::Array{T,1}
    dxyg::Array{T,1}
    xcmin::Array{T,1}
    xcmax::Array{T,1}
    ycmin::Array{T,1}
    ycmax::Array{T,1}
    dxyc::Array{T,1}
    dyoffc::Array{T,1}
    xpmin::Array{T,1}
    xpmax::Array{T,1}
    ypmin::Array{T,1}
    ypmax::Array{T,1}
    dxyp::Array{T,1}
    dyoffp::Array{T,1}
    ysfp::Array{T,1}
    gtick::Array{T,1}
    #qmat
    q::Array{T,2}
    dq::Array{T,1}
    dzdg::Array{T,1}
    dzdn::Array{T,1}
    dzdm::Array{T,1}
    dqdg::Array{T,1}
    dqdm::Array{T,1}
    qtan1::Array{T,1}
    qtan2::Array{T,1}
    z_qinf::Array{T,1}
    z_alfa::Array{T,1}
    z_qdof0::Array{T,1}
    z_qdof1::Array{T,1}
    z_qdof2::Array{T,1}
    z_qdof3::Array{T,1}
    #vmat
    va::Array{T,3}
    vb::Array{T,3}
    vdel::Array{T,3}
    vm::Array{T,3}
    vz::Array{T,2}
    #ci01
    iq1::Array{Int32,1}
    iq2::Array{Int32,1}
    nsp::Array{Int32,1}
    nqsp::Array{Int32,1}
    kqtarg::Array{Int32,1}
    iacqsp::Array{Int32,1}
    nc1::Array{Int32,1}
    nname::Array{Int32,1}
    nprefix::Array{Int32,1}
    #ci03
    ncpref::Array{Int32,1}
    napol::Array{Int32,1}
    npol::Array{Int32,1}
    ipact::Array{Int32,1}
    nlref::Array{Int32,1}
    icolp::Array{Int32,1}
    icolr::Array{Int32,1}
    imatyp::Array{Int32,1}
    iretyp::Array{Int32,1}
    nxypol::Array{Int32,1}
    npolref::Array{Int32,1}
    ndref::Array{Int32,2}
    #ci04
    n::Array{Int32,1}
    nb::Array{Int32,1}
    nw::Array{Int32,1}
    npan::Array{Int32,1}
    ist::Array{Int32,1}
    kimage::Array{Int32,1}
    itmax::Array{Int32,1}
    nseqex::Array{Int32,1}
    retyp::Array{Int32,1}
    matyp::Array{Int32,1}
    aijpix::Array{Int32,1}
    idev::Array{Int32,1}
    idevrp::Array{Int32,1}
    ipslu::Array{Int32,1}
    ncolor::Array{Int32,1}
    icols::Array{Int32,1}
    nover::Array{Int32,1}
    ncm::Array{Int32,1}
    ntk::Array{Int32,1}
    #ci05
    iblte::Array{Int32,1}
    nbl::Array{Int32,1}
    ipan::Array{Int32,2}
    isys::Array{Int32,2}
    nsys::Array{Int32,1}
    itran::Array{Int32,1}
    #ci06
    imxbl::Array{Int32,1}
    ismxbl::Array{Int32,1}
    #cl01
    ok::Array{Int32,1}
    limage::Array{Int32,1}
    sharp::Array{Int32,1}
    lgamu::Array{Int32,1}
    lqinu::Array{Int32,1}
    lvisc::Array{Int32,1}
    lalfa::Array{Int32,1}
    lwake::Array{Int32,1}
    lpacc::Array{Int32,1}
    lblini::Array{Int32,1}
    lipan::Array{Int32,1}
    lqaij::Array{Int32,1}
    ladij::Array{Int32,1}
    lwdij::Array{Int32,1}
    lcpxx::Array{Int32,1}
    lqvdes::Array{Int32,1}
    lqrefl::Array{Int32,1}
    lqspec::Array{Int32,1}
    lvconv::Array{Int32,1}
    lcpref::Array{Int32,1}
    lclock::Array{Int32,1}
    lpfile::Array{Int32,1}
    lpfilx::Array{Int32,1}
    lppsho::Array{Int32,1}
    lbflap::Array{Int32,1}
    lflap::Array{Int32,1}
    leiw::Array{Int32,1}
    lscini::Array{Int32,1}
    lforef::Array{Int32,1}
    lnorm::Array{Int32,1}
    lgsame::Array{Int32,1}
    lplcam::Array{Int32,1}
    lqsym::Array{Int32,1}
    lgsym::Array{Int32,1}
    lqgrid::Array{Int32,1}
    lggrid::Array{Int32,1}
    lgtick::Array{Int32,1}
    lqslop::Array{Int32,1}
    lgslop::Array{Int32,1}
    lcslop::Array{Int32,1}
    lqsppl::Array{Int32,1}
    lgeopl::Array{Int32,1}
    lgparm::Array{Int32,1}
    lcpgrd::Array{Int32,1}
    lblgrd::Array{Int32,1}
    lblsym::Array{Int32,1}
    lplot::Array{Int32,1}
    lsym::Array{Int32,1}
    liqset::Array{Int32,1}
    lclip::Array{Int32,1}
    lvlab::Array{Int32,1}
    lcurs::Array{Int32,1}
    lland::Array{Int32,1}
    lpgrid::Array{Int32,1}
    lpcdw::Array{Int32,1}
    lplist::Array{Int32,1}
    lplegn::Array{Int32,1}
    lcminp::Array{Int32,1}
    lhmomp::Array{Int32,1}
    lexitflag::Array{Int32,1}
    #cl02
    tforce::Array{Int32,1}
    #cc01
    fname::String
    name::String
    ispars::String
    oname::String
    prefix::String
    pfname::String
    pfnamx::String
    namepol::String
    nameref::String
    #cc02
    labref::String
    #cc03
    vmxbl::String
end

"""
    get_globals()

Wraps the global variables in XFOIL in a struct
"""
get_globals

"""
    get_globals_cs()

`get_globals` for the complex step enabled version of XFOIL.
"""
get_globals_cs

# definition of get_globals
for (T, name, library) in
    ((:Float64, :get_globals, :libxfoil_light),
    (:ComplexF64, :get_globals_cs, :libxfoil_light_cs))

    @eval begin

        function $(name)()

            cr01 = cglobal((:cr01_, $(library)),$T)
            version = unsafe_wrap(Array,cr01,1)

            cr03 = cglobal((:cr03_, $(library)),$T)
            aij = unsafe_wrap(Array,cr03,(IQX,IQX))
            dij = unsafe_wrap(Array,cr03+sizeof($T)*IQX^2,(IZX,IZX))

            cr04 = cglobal((:cr04_, $(library)),$T)
            qinv = unsafe_wrap(Array,cr04,IZX)
            qvis = unsafe_wrap(Array,cr04+sizeof($T)*IZX,IZX)
            cpi = unsafe_wrap(Array,cr04+sizeof($T)*(2*IZX),IZX)
            cpv = unsafe_wrap(Array,cr04+sizeof($T)*(3*IZX),IZX)
            qinvu = unsafe_wrap(Array,cr04+sizeof($T)*(4*IZX),(IZX,2))
            qinv_a = unsafe_wrap(Array,cr04+sizeof($T)*(6*IZX),IZX)

            cr05 = cglobal((:cr05_, $(library)),$T)
            x = unsafe_wrap(Array,cr05,IZX)
            y = unsafe_wrap(Array,cr05+sizeof($T)*IZX,IZX)
            xp = unsafe_wrap(Array,cr05+sizeof($T)*(2*IZX),IZX)
            yp = unsafe_wrap(Array,cr05+sizeof($T)*(3*IZX),IZX)
            s = unsafe_wrap(Array,cr05+sizeof($T)*(4*IZX),IZX)
            sle = unsafe_wrap(Array,cr05+sizeof($T)*(5*IZX),1)
            xle = unsafe_wrap(Array,cr05+sizeof($T)*(5*IZX+1),1)
            yle = unsafe_wrap(Array,cr05+sizeof($T)*(5*IZX+2),1)
            xte = unsafe_wrap(Array,cr05+sizeof($T)*(5*IZX+3),1)
            yte = unsafe_wrap(Array,cr05+sizeof($T)*(5*IZX+4),1)
            chord = unsafe_wrap(Array,cr05+sizeof($T)*(5*IZX+5),1)
            yimage = unsafe_wrap(Array,cr05+sizeof($T)*(5*IZX+6),1)
            wgap = unsafe_wrap(Array,cr05+sizeof($T)*(5*IZX+7),IWK)
            waklen = unsafe_wrap(Array,cr05+sizeof($T)*(5*IZX+7+IWK),1)

            cr06 = cglobal((:cr06_, $(library)),$T)
            gam = unsafe_wrap(Array,cr06,IQX)
            gamu = unsafe_wrap(Array,cr06+sizeof($T)*IQX,(IQX,2))
            gam_a = unsafe_wrap(Array,cr06+sizeof($T)*(IQX*3),IQX)
            sig = unsafe_wrap(Array,cr06+sizeof($T)*(IQX*4),IZX)
            nx = unsafe_wrap(Array,cr06+sizeof($T)*(IQX*4+IZX),IZX)
            ny = unsafe_wrap(Array,cr06+sizeof($T)*(IQX*4+IZX*2),IZX)
            apanel = unsafe_wrap(Array,cr06+sizeof($T)*(IQX*4+IZX*3),IZX)
            sst = unsafe_wrap(Array,cr06+sizeof($T)*(IQX*4+IZX*4),1)
            sst_go = unsafe_wrap(Array,cr06+sizeof($T)*(IQX*4+IZX*4+1),1)
            sst_gp = unsafe_wrap(Array,cr06+sizeof($T)*(IQX*4+IZX*4+2),1)
            gamte = unsafe_wrap(Array,cr06+sizeof($T)*(IQX*4+IZX*4+3),1)
            gamte_a = unsafe_wrap(Array,cr06+sizeof($T)*(IQX*4+IZX*4+4),1)
            sigte = unsafe_wrap(Array,cr06+sizeof($T)*(IQX*4+IZX*4+5),1)
            sigte_a = unsafe_wrap(Array,cr06+sizeof($T)*(IQX*4+IZX*4+6),1)
            dste = unsafe_wrap(Array,cr06+sizeof($T)*(IQX*4+IZX*4+7),1)
            ante = unsafe_wrap(Array,cr06+sizeof($T)*(IQX*4+IZX*4+8),1)
            aste = unsafe_wrap(Array,cr06+sizeof($T)*(IQX*4+IZX*4+9),1)

            cr07 = cglobal((:cr07_, $(library)),$T)
            ssple = unsafe_wrap(Array,cr07,1)
            sspec = unsafe_wrap(Array,cr07+sizeof($T),IBX)
            xspoc = unsafe_wrap(Array,cr07+sizeof($T)*(1+IBX),IBX)
            yspoc = unsafe_wrap(Array,cr07+sizeof($T)*(1+2*IBX),IBX)
            qgamm = unsafe_wrap(Array,cr07+sizeof($T)*(1+3*IBX),IBX)
            qspec = unsafe_wrap(Array,cr07+sizeof($T)*(1+4*IBX),(IBX,IPX))
            qspecp = unsafe_wrap(Array,cr07+sizeof($T)*(1+4*IBX+IBX*IPX),(IBX,IPX))
            algam = unsafe_wrap(Array,cr07+sizeof($T)*(1+4*IBX+2*IBX*IPX),1)
            clgam = unsafe_wrap(Array,cr07+sizeof($T)*(2+4*IBX+2*IBX*IPX),1)
            cmgam = unsafe_wrap(Array,cr07+sizeof($T)*(3+4*IBX+2*IBX*IPX),1)
            alqsp = unsafe_wrap(Array,cr07+sizeof($T)*(4+4*IBX+2*IBX*IPX),IPX)
            clqsp = unsafe_wrap(Array,cr07+sizeof($T)*(4+4*IBX+2*IBX*IPX+IPX),IPX)
            cmqsp = unsafe_wrap(Array,cr07+sizeof($T)*(4+4*IBX+2*IBX*IPX+2*IPX),IPX)
            qf0 = unsafe_wrap(Array,cr07+sizeof($T)*(4+4*IBX+2*IBX*IPX+3*IPX),IQX)
            qf1 = unsafe_wrap(Array,cr07+sizeof($T)*(4+4*IBX+2*IBX*IPX+3*IPX+IQX),IQX)
            qf2 = unsafe_wrap(Array,cr07+sizeof($T)*(4+4*IBX+2*IBX*IPX+3*IPX+2*IQX),IQX)
            qf3 = unsafe_wrap(Array,cr07+sizeof($T)*(4+4*IBX+2*IBX*IPX+3*IPX+3*IQX),IQX)
            qdof0 = unsafe_wrap(Array,cr07+sizeof($T)*(4+4*IBX+2*IBX*IPX+3*IPX+4*IQX),1)
            qdof1 = unsafe_wrap(Array,cr07+sizeof($T)*(5+4*IBX+2*IBX*IPX+3*IPX+4*IQX),1)
            qdof2 = unsafe_wrap(Array,cr07+sizeof($T)*(6+4*IBX+2*IBX*IPX+3*IPX+4*IQX),1)
            qdof3 = unsafe_wrap(Array,cr07+sizeof($T)*(7+4*IBX+2*IBX*IPX+3*IPX+4*IQX),1)
            clspec = unsafe_wrap(Array,cr07+sizeof($T)*(8+4*IBX+2*IBX*IPX+3*IPX+4*IQX),1)
            ffilt = unsafe_wrap(Array,cr07+sizeof($T)*(9+4*IBX+2*IBX*IPX+3*IPX+4*IQX),1)

            cr09 = cglobal((:cr09_, $(library)),$T)
            adeg = unsafe_wrap(Array,cr09+0*sizeof($T),1)
            alfa = unsafe_wrap(Array,cr09+1*sizeof($T),1)
            awake = unsafe_wrap(Array,cr09+2*sizeof($T),1)
            mvisc = unsafe_wrap(Array,cr09+3*sizeof($T),1)
            avisc = unsafe_wrap(Array,cr09+4*sizeof($T),1)
            xcmref = unsafe_wrap(Array,cr09+5*sizeof($T),1)
            ycmref = unsafe_wrap(Array,cr09+6*sizeof($T),1)
            cl = unsafe_wrap(Array,cr09+7*sizeof($T),1)
            cm = unsafe_wrap(Array,cr09+8*sizeof($T),1)
            cd = unsafe_wrap(Array,cr09+9*sizeof($T),1)
            cdp = unsafe_wrap(Array,cr09+10*sizeof($T),1)
            cdf = unsafe_wrap(Array,cr09+11*sizeof($T),1)
            cl_alf = unsafe_wrap(Array,cr09+12*sizeof($T),1)
            cl_msq = unsafe_wrap(Array,cr09+13*sizeof($T),1)
            psio = unsafe_wrap(Array,cr09+14*sizeof($T),1)
            circ = unsafe_wrap(Array,cr09+15*sizeof($T),1)
            cosa = unsafe_wrap(Array,cr09+16*sizeof($T),1)
            sina = unsafe_wrap(Array,cr09+17*sizeof($T),1)
            qinf = unsafe_wrap(Array,cr09+18*sizeof($T),1)
            gamma = unsafe_wrap(Array,cr09+19*sizeof($T),1)
            gamm1 = unsafe_wrap(Array,cr09+20*sizeof($T),1)
            minf1 = unsafe_wrap(Array,cr09+21*sizeof($T),1)
            minf = unsafe_wrap(Array,cr09+22*sizeof($T),1)
            minf_cl = unsafe_wrap(Array,cr09+23*sizeof($T),1)
            tklam = unsafe_wrap(Array,cr09+24*sizeof($T),1)
            tkl_msq = unsafe_wrap(Array,cr09+25*sizeof($T),1)
            cpstar = unsafe_wrap(Array,cr09+26*sizeof($T),1)
            qstar = unsafe_wrap(Array,cr09+27*sizeof($T),1)
            cpmn = unsafe_wrap(Array,cr09+28*sizeof($T),1)
            cpmni = unsafe_wrap(Array,cr09+29*sizeof($T),1)
            cpmnv = unsafe_wrap(Array,cr09+30*sizeof($T),1)
            xcpmni = unsafe_wrap(Array,cr09+31*sizeof($T),1)
            xcpmnv = unsafe_wrap(Array,cr09+32*sizeof($T),1)

            cr10 = cglobal((:cr10_, $(library)),$T)
            xpref = unsafe_wrap(Array,cr10,IQX)
            cpref = unsafe_wrap(Array,cr10+sizeof($T)*IQX,IQX)
            verspol = unsafe_wrap(Array,cr10+sizeof($T)*(2*IQX),NPX)
            cpolxy = unsafe_wrap(Array,cr10+sizeof($T)*(2*IQX+NPX),(IQX,2,NPX))
            machp1 = unsafe_wrap(Array,cr10+sizeof($T)*(2*IQX+NPX+IQX*2*NPX),NPX)
            reynp1 = unsafe_wrap(Array,cr10+sizeof($T)*(2*IQX+2*NPX+IQX*2*NPX),NPX)
            acritp = unsafe_wrap(Array,cr10+sizeof($T)*(2*IQX+3*NPX+IQX*2*NPX),NPX)
            xstripp = unsafe_wrap(Array,cr10+sizeof($T)*(2*IQX+4*NPX+IQX*2*NPX),(ISX,NPX))

            cr11= cglobal((:cr11_, $(library)),$T)
            pi = unsafe_wrap(Array,cr11+sizeof($T)*0,1)
            hopi = unsafe_wrap(Array,cr11+sizeof($T)*1,1)
            qopi = unsafe_wrap(Array,cr11+sizeof($T)*2,1)
            dtor = unsafe_wrap(Array,cr11+sizeof($T)*3,1)

            cr12= cglobal((:cr12_, $(library)),$T)
            cvpar = unsafe_wrap(Array,cr12+sizeof($T)*0,1)
            cterat = unsafe_wrap(Array,cr12+sizeof($T)*1,1)
            ctrrat = unsafe_wrap(Array,cr12+sizeof($T)*2,1)
            xsref1 = unsafe_wrap(Array,cr12+sizeof($T)*3,1)
            xsref2 = unsafe_wrap(Array,cr12+sizeof($T)*4,1)
            xpref1 = unsafe_wrap(Array,cr12+sizeof($T)*5,1)
            xpref2 = unsafe_wrap(Array,cr12+sizeof($T)*6,1)

            cr13 = cglobal((:cr12_, $(library)),$T)
            size = unsafe_wrap(Array,cr13+sizeof($T)*0,1)
            scrnfr = unsafe_wrap(Array,cr13+sizeof($T)*1,1)
            plotar = unsafe_wrap(Array,cr13+sizeof($T)*2,1)
            pfac = unsafe_wrap(Array,cr13+sizeof($T)*3,1)
            qfac = unsafe_wrap(Array,cr13+sizeof($T)*4,1)
            vfac = unsafe_wrap(Array,cr13+sizeof($T)*5,1)
            xwind = unsafe_wrap(Array,cr13+sizeof($T)*6,1)
            ywind = unsafe_wrap(Array,cr13+sizeof($T)*7,1)
            xpage = unsafe_wrap(Array,cr13+sizeof($T)*8,1)
            ypage = unsafe_wrap(Array,cr13+sizeof($T)*9,1)
            xmarg = unsafe_wrap(Array,cr13+sizeof($T)*10,1)
            ymarg = unsafe_wrap(Array,cr13+sizeof($T)*11,1)
            ch = unsafe_wrap(Array,cr13+sizeof($T)*12,1)
            chg = unsafe_wrap(Array,cr13+sizeof($T)*13,1)
            chq = unsafe_wrap(Array,cr13+sizeof($T)*14,1)
            xofair = unsafe_wrap(Array,cr13+sizeof($T)*15,1)
            yofair = unsafe_wrap(Array,cr13+sizeof($T)*16,1)
            facair = unsafe_wrap(Array,cr13+sizeof($T)*17,1)
            xofa = unsafe_wrap(Array,cr13+sizeof($T)*18,1)
            yofa = unsafe_wrap(Array,cr13+sizeof($T)*19,1)
            faca = unsafe_wrap(Array,cr13+sizeof($T)*20,1)
            uprwt = unsafe_wrap(Array,cr13+sizeof($T)*21,1)
            cpmin = unsafe_wrap(Array,cr13+sizeof($T)*22,1)
            cpmax = unsafe_wrap(Array,cr13+sizeof($T)*23,1)
            cpdel = unsafe_wrap(Array,cr13+sizeof($T)*24,1)
            cpolplf = unsafe_wrap(Array,cr13+sizeof($T)*25,(3,4))
            xcdwid = unsafe_wrap(Array,cr13+sizeof($T)*(25+3*4),1)
            xalwid = unsafe_wrap(Array,cr13+sizeof($T)*(26+3*4),1)
            xocwid = unsafe_wrap(Array,cr13+sizeof($T)*(27+3*4),1)

            cr14 = cglobal((:cr14_, $(library)),$T)
            xb = unsafe_wrap(Array,cr14,IBX)
            yb = unsafe_wrap(Array,cr14+sizeof($T)*IBX,IBX)
            xbp = unsafe_wrap(Array,cr14+sizeof($T)*(2*IBX),IBX)
            ybp = unsafe_wrap(Array,cr14+sizeof($T)*(3*IBX),IBX)
            sb = unsafe_wrap(Array,cr14+sizeof($T)*(4*IBX),IBX)
            snew = unsafe_wrap(Array,cr14+sizeof($T)*(5*IBX),4*IBX)
            xbf = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX),1)
            ybf = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+1),1)
            xof = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+2),1)
            yof = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+3),1)
            hmom = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+4),1)
            hfx = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+5),1)
            hfy = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+6),1)
            xbmin = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+7),1)
            xbmax = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+8),1)
            ybmin = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+9),1)
            ybmax = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+10),1)
            sble = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+11),1)
            chordb = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+12),1)
            areab = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+13),1)
            radble = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+14),1)
            angbte = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+15),1)
            ei11ba = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+16),1)
            ei22ba = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+17),1)
            apx1ba = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+18),1)
            apx2ba = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+19),1)
            ei11bt = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+20),1)
            ei22bt = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+21),1)
            apx1bt = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+22),1)
            apx2bt = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+23),1)
            thickb = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+24),1)
            cambrb = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+25),1)
            xcm = unsafe_wrap(Array,cr14+sizeof($T)*(9*IBX+26),2*IBX)
            ycm = unsafe_wrap(Array,cr14+sizeof($T)*(11*IBX+26),2*IBX)
            scm = unsafe_wrap(Array,cr14+sizeof($T)*(13*IBX+26),2*IBX)
            xcmp = unsafe_wrap(Array,cr14+sizeof($T)*(15*IBX+26),2*IBX)
            ycmp = unsafe_wrap(Array,cr14+sizeof($T)*(17*IBX+36),2*IBX)
            xtk = unsafe_wrap(Array,cr14+sizeof($T)*(19*IBX+36),2*IBX)
            ytk = unsafe_wrap(Array,cr14+sizeof($T)*(21*IBX+36),2*IBX)
            stk = unsafe_wrap(Array,cr14+sizeof($T)*(23*IBX+36),2*IBX)
            xtkp = unsafe_wrap(Array,cr14+sizeof($T)*(25*IBX+36),2*IBX)
            ytkp = unsafe_wrap(Array,cr14+sizeof($T)*(27*IBX+36),2*IBX)

            cr15 = cglobal((:cr15_, $(library)),$T)
            xssi = unsafe_wrap(Array,cr15,(IVX,ISX))
            uedg = unsafe_wrap(Array,cr15+sizeof($T)*(IVX*ISX),(IVX,ISX))
            uinv = unsafe_wrap(Array,cr15+sizeof($T)*(2*IVX*ISX),(IVX,ISX))
            mass = unsafe_wrap(Array,cr15+sizeof($T)*(3*IVX*ISX),(IVX,ISX))
            thet = unsafe_wrap(Array,cr15+sizeof($T)*(4*IVX*ISX),(IVX,ISX))
            dstr = unsafe_wrap(Array,cr15+sizeof($T)*(5*IVX*ISX),(IVX,ISX))
            ctau = unsafe_wrap(Array,cr15+sizeof($T)*(6*IVX*ISX),(IVX,ISX))
            delt = unsafe_wrap(Array,cr15+sizeof($T)*(7*IVX*ISX),(IVX,ISX))
            uslp = unsafe_wrap(Array,cr15+sizeof($T)*(8*IVX*ISX),(IVX,ISX))
            guxq = unsafe_wrap(Array,cr15+sizeof($T)*(9*IVX*ISX),(IVX,ISX))
            guxd = unsafe_wrap(Array,cr15+sizeof($T)*(10*IVX*ISX),(IVX,ISX))
            tau = unsafe_wrap(Array,cr15+sizeof($T)*(11*IVX*ISX),(IVX,ISX))
            dis = unsafe_wrap(Array,cr15+sizeof($T)*(12*IVX*ISX),(IVX,ISX))
            ctq = unsafe_wrap(Array,cr15+sizeof($T)*(13*IVX*ISX),(IVX,ISX))
            vti = unsafe_wrap(Array,cr15+sizeof($T)*(14*IVX*ISX),(IVX,ISX))
            reinf1 = unsafe_wrap(Array,cr15+sizeof($T)*(15*IVX*ISX),1)
            reinf = unsafe_wrap(Array,cr15+sizeof($T)*(15*IVX*ISX+1),1)
            reinf_cl = unsafe_wrap(Array,cr15+sizeof($T)*(15*IVX*ISX+2),1)
            acrit = unsafe_wrap(Array,cr15+sizeof($T)*(15*IVX*ISX+3),1)
            xstrip = unsafe_wrap(Array,cr15+sizeof($T)*(15*IVX*ISX+4),ISX)
            xoctr = unsafe_wrap(Array,cr15+sizeof($T)*(15*IVX*ISX+4+ISX),ISX)
            yoctr = unsafe_wrap(Array,cr15+sizeof($T)*(15*IVX*ISX+4+2*ISX),ISX)
            xssitr = unsafe_wrap(Array,cr15+sizeof($T)*(15*IVX*ISX+4+3*ISX),ISX)
            uinv_a = unsafe_wrap(Array,cr15+sizeof($T)*(15*IVX*ISX+4+4*ISX),(IVX,ISX))

            cr17 = cglobal((:cr17_, $(library)),$T)
            rmsbl = unsafe_wrap(Array,cr17,1)
            rmxbl = unsafe_wrap(Array,cr17+sizeof($T),1)
            rlx = unsafe_wrap(Array,cr17+2*sizeof($T),1)
            vaccel = unsafe_wrap(Array,cr17+3*sizeof($T),1)

            cr18 = cglobal((:cr18_, $(library)),$T)
            xsf = unsafe_wrap(Array,cr18,1)
            ysf = unsafe_wrap(Array,cr18+sizeof($T),1)
            xoff = unsafe_wrap(Array,cr18+sizeof($T)*2,1)
            yoff = unsafe_wrap(Array,cr18+sizeof($T)*3,1)
            xgmin = unsafe_wrap(Array,cr18+sizeof($T)*4,1)
            xgmax = unsafe_wrap(Array,cr18+sizeof($T)*5,1)
            ygmin = unsafe_wrap(Array,cr18+sizeof($T)*6,1)
            ygmax = unsafe_wrap(Array,cr18+sizeof($T)*7,1)
            dxyg = unsafe_wrap(Array,cr18+sizeof($T)*8,1)
            xcmin = unsafe_wrap(Array,cr18+sizeof($T)*9,1)
            xcmax = unsafe_wrap(Array,cr18+sizeof($T)*10,1)
            ycmin = unsafe_wrap(Array,cr18+sizeof($T)*11,1)
            ycmax = unsafe_wrap(Array,cr18+sizeof($T)*12,1)
            dxyc = unsafe_wrap(Array,cr18+sizeof($T)*13,1)
            dyoffc = unsafe_wrap(Array,cr18+sizeof($T)*14,1)
            xpmin = unsafe_wrap(Array,cr18+sizeof($T)*15,1)
            xpmax = unsafe_wrap(Array,cr18+sizeof($T)*16,1)
            ypmin = unsafe_wrap(Array,cr18+sizeof($T)*17,1)
            ypmax = unsafe_wrap(Array,cr18+sizeof($T)*18,1)
            dxyp = unsafe_wrap(Array,cr18+sizeof($T)*19,1)
            dyoffp = unsafe_wrap(Array,cr18+sizeof($T)*20,1)
            ysfp = unsafe_wrap(Array,cr18+sizeof($T)*21,1)
            gtick = unsafe_wrap(Array,cr18+sizeof($T)*22,1)

            qmat = cglobal((:qmat_, $(library)),$T)
            q = unsafe_wrap(Array,qmat,(IQX,IQX))
            dq = unsafe_wrap(Array,qmat+sizeof($T)*(IQX*IQX),IQX)
            dzdg = unsafe_wrap(Array,qmat+sizeof($T)*(IQX^2+IQX),IQX)
            dzdn = unsafe_wrap(Array,qmat+sizeof($T)*(IQX^2+2*IQX),IQX)
            dzdm = unsafe_wrap(Array,qmat+sizeof($T)*(IQX^2+3*IQX),IZX)
            dqdg = unsafe_wrap(Array,qmat+sizeof($T)*(IQX^2+3*IQX+IZX),IQX)
            dqdm = unsafe_wrap(Array,qmat+sizeof($T)*(IQX^2+4*IQX+IZX),IZX)
            qtan1 = unsafe_wrap(Array,qmat+sizeof($T)*(IQX^2+4*IQX+2*IZX),1)
            qtan2 = unsafe_wrap(Array,qmat+sizeof($T)*(IQX^2+4*IQX+2*IZX+1),1)
            z_qinf = unsafe_wrap(Array,qmat+sizeof($T)*(IQX^2+4*IQX+2*IZX+2),1)
            z_alfa = unsafe_wrap(Array,qmat+sizeof($T)*(IQX^2+4*IQX+2*IZX+3),1)
            z_qdof0 = unsafe_wrap(Array,qmat+sizeof($T)*(IQX^2+4*IQX+2*IZX+4),1)
            z_qdof1 = unsafe_wrap(Array,qmat+sizeof($T)*(IQX^2+4*IQX+2*IZX+5),1)
            z_qdof2 = unsafe_wrap(Array,qmat+sizeof($T)*(IQX^2+4*IQX+2*IZX+6),1)
            z_qdof3 = unsafe_wrap(Array,qmat+sizeof($T)*(IQX^2+4*IQX+2*IZX+7),1)

            vmat = cglobal((:vmat_, $(library)),$T)
            va = unsafe_wrap(Array,vmat,(3,2,IZX))
            vb = unsafe_wrap(Array,vmat+sizeof($T)*(3*2*IZX),(3,2,IZX))
            vdel = unsafe_wrap(Array,vmat+sizeof($T)*(2*3*2*IZX),(3,2,IZX))
            vm = unsafe_wrap(Array,vmat+sizeof($T)*(3*3*2*IZX),(3,IZX,IZX))
            vz = unsafe_wrap(Array,vmat+sizeof($T)*(3*3*2*IZX+IZX*IZX),(3,2))

            ci01 = cglobal((:ci01_, $(library)),Int32)
            iq1 = unsafe_wrap(Array,ci01,1)
            iq2 = unsafe_wrap(Array,ci01+sizeof(Int32),1)
            nsp = unsafe_wrap(Array,ci01+sizeof(Int32)*2,1)
            nqsp = unsafe_wrap(Array,ci01+sizeof(Int32)*3,1)
            kqtarg = unsafe_wrap(Array,ci01+sizeof(Int32)*4,1)
            iacqsp = unsafe_wrap(Array,ci01+sizeof(Int32)*5,1)
            nc1 = unsafe_wrap(Array,ci01+sizeof(Int32)*6,1)
            nname = unsafe_wrap(Array,ci01+sizeof(Int32)*7,1)
            nprefix = unsafe_wrap(Array,ci01+sizeof(Int32)*8,1)

            ci03 = cglobal((:ci03_, $(library)),Int32)
            ncpref = unsafe_wrap(Array,ci03,1)
            napol = unsafe_wrap(Array,ci03+sizeof(Int32),NPX)
            npol = unsafe_wrap(Array,ci03+sizeof(Int32)*(1+NPX),1)
            ipact = unsafe_wrap(Array,ci03+sizeof(Int32)*(2+NPX),1)
            nlref = unsafe_wrap(Array,ci03+sizeof(Int32)*(3+NPX),1)
            icolp = unsafe_wrap(Array,ci03+sizeof(Int32)*(4+NPX),NPX)
            icolr = unsafe_wrap(Array,ci03+sizeof(Int32)*(4+2*NPX),NPX)
            imatyp = unsafe_wrap(Array,ci03+sizeof(Int32)*(4+3*NPX),NPX)
            iretyp = unsafe_wrap(Array,ci03+sizeof(Int32)*(4+4*NPX),NPX)
            nxypol = unsafe_wrap(Array,ci03+sizeof(Int32)*(4+5*NPX),NPX)
            npolref = unsafe_wrap(Array,ci03+sizeof(Int32)*(4+6*NPX),1)
            ndref = unsafe_wrap(Array,ci03+sizeof(Int32)*(5+6*NPX),(4,NPX))

            ci04 = cglobal((:ci04_, $(library)),Int32)
            n = unsafe_wrap(Array,ci04,1)
            nb = unsafe_wrap(Array,ci04+sizeof(Int32)*1,1)
            nw = unsafe_wrap(Array,ci04+sizeof(Int32)*2,1)
            npan = unsafe_wrap(Array,ci04+sizeof(Int32)*3,1)
            ist = unsafe_wrap(Array,ci04+sizeof(Int32)*4,1)
            kimage = unsafe_wrap(Array,ci04+sizeof(Int32)*5,1)
            itmax = unsafe_wrap(Array,ci04+sizeof(Int32)*6,1)
            nseqex = unsafe_wrap(Array,ci04+sizeof(Int32)*7,1)
            retyp = unsafe_wrap(Array,ci04+sizeof(Int32)*8,1)
            matyp = unsafe_wrap(Array,ci04+sizeof(Int32)*9,1)
            aijpix = unsafe_wrap(Array,ci04+sizeof(Int32)*10,ISX)
            idev = unsafe_wrap(Array,ci04+sizeof(Int32)*(10+ISX),1)
            idevrp = unsafe_wrap(Array,ci04+sizeof(Int32)*(11+ISX),1)
            ipslu = unsafe_wrap(Array,ci04+sizeof(Int32)*(12+ISX),1)
            ncolor = unsafe_wrap(Array,ci04+sizeof(Int32)*(13+ISX),1)
            icols = unsafe_wrap(Array,ci04+sizeof(Int32)*(14+ISX),ISX)
            nover = unsafe_wrap(Array,ci04+sizeof(Int32)*(14+2*ISX),1)
            ncm = unsafe_wrap(Array,ci04+sizeof(Int32)*(15+2*ISX),1)
            ntk = unsafe_wrap(Array,ci04+sizeof(Int32)*(16+2*ISX),1)

            ci05 = cglobal((:ci05_, $(library)),Int32)
            iblte = unsafe_wrap(Array,ci05,ISX)
            nbl = unsafe_wrap(Array,ci05+sizeof(Int32)*ISX,ISX)
            ipan = unsafe_wrap(Array,ci05+sizeof(Int32)*(2*ISX),(IVX,ISX))
            isys = unsafe_wrap(Array,ci05+sizeof(Int32)*(2*ISX+IVX*ISX),(IVX,ISX))
            nsys = unsafe_wrap(Array,ci05+sizeof(Int32)*(2*ISX+2*IVX*ISX),1)
            itran = unsafe_wrap(Array,ci05+sizeof(Int32)*(2*ISX+2*IVX*ISX+1),ISX)

            ci06 = cglobal((:ci06_, $(library)),Int32)
            imxbl = unsafe_wrap(Array,ci06,1)
            ismxbl = unsafe_wrap(Array,ci06+sizeof(Int32),1)

            cl01 = cglobal((:cl01_, $(library)),Int32)
            ok = unsafe_wrap(Array,cl01+sizeof(Int32)*0,1)
            limage = unsafe_wrap(Array,cl01+sizeof(Int32)*1,1)
            sharp = unsafe_wrap(Array,cl01+sizeof(Int32)*2,1)
            lgamu = unsafe_wrap(Array,cl01+sizeof(Int32)*3,1)
            lqinu = unsafe_wrap(Array,cl01+sizeof(Int32)*4,1)
            lvisc = unsafe_wrap(Array,cl01+sizeof(Int32)*5,1)
            lalfa = unsafe_wrap(Array,cl01+sizeof(Int32)*6,1)
            lwake = unsafe_wrap(Array,cl01+sizeof(Int32)*7,1)
            lpacc = unsafe_wrap(Array,cl01+sizeof(Int32)*8,1)
            lblini = unsafe_wrap(Array,cl01+sizeof(Int32)*9,1)
            lipan = unsafe_wrap(Array,cl01+sizeof(Int32)*10,1)
            lqaij = unsafe_wrap(Array,cl01+sizeof(Int32)*11,1)
            ladij = unsafe_wrap(Array,cl01+sizeof(Int32)*12,1)
            lwdij = unsafe_wrap(Array,cl01+sizeof(Int32)*13,1)
            lcpxx = unsafe_wrap(Array,cl01+sizeof(Int32)*14,1)
            lqvdes = unsafe_wrap(Array,cl01+sizeof(Int32)*15,1)
            lqrefl = unsafe_wrap(Array,cl01+sizeof(Int32)*16,1)
            lqspec = unsafe_wrap(Array,cl01+sizeof(Int32)*17,1)
            lvconv = unsafe_wrap(Array,cl01+sizeof(Int32)*18,1)
            lcpref = unsafe_wrap(Array,cl01+sizeof(Int32)*19,1)
            lclock = unsafe_wrap(Array,cl01+sizeof(Int32)*20,1)
            lpfile = unsafe_wrap(Array,cl01+sizeof(Int32)*21,1)
            lpfilx = unsafe_wrap(Array,cl01+sizeof(Int32)*22,1)
            lppsho = unsafe_wrap(Array,cl01+sizeof(Int32)*23,1)
            lbflap = unsafe_wrap(Array,cl01+sizeof(Int32)*24,1)
            lflap = unsafe_wrap(Array,cl01+sizeof(Int32)*25,1)
            leiw = unsafe_wrap(Array,cl01+sizeof(Int32)*26,1)
            lscini = unsafe_wrap(Array,cl01+sizeof(Int32)*27,1)
            lforef = unsafe_wrap(Array,cl01+sizeof(Int32)*28,1)
            lnorm = unsafe_wrap(Array,cl01+sizeof(Int32)*29,1)
            lgsame = unsafe_wrap(Array,cl01+sizeof(Int32)*30,1)
            lplcam = unsafe_wrap(Array,cl01+sizeof(Int32)*31,1)
            lqsym = unsafe_wrap(Array,cl01+sizeof(Int32)*32,1)
            lgsym = unsafe_wrap(Array,cl01+sizeof(Int32)*33,1)
            lqgrid = unsafe_wrap(Array,cl01+sizeof(Int32)*34,1)
            lggrid = unsafe_wrap(Array,cl01+sizeof(Int32)*35,1)
            lgtick = unsafe_wrap(Array,cl01+sizeof(Int32)*36,1)
            lqslop = unsafe_wrap(Array,cl01+sizeof(Int32)*37,1)
            lgslop = unsafe_wrap(Array,cl01+sizeof(Int32)*38,1)
            lcslop = unsafe_wrap(Array,cl01+sizeof(Int32)*39,1)
            lqsppl = unsafe_wrap(Array,cl01+sizeof(Int32)*40,1)
            lgeopl = unsafe_wrap(Array,cl01+sizeof(Int32)*41,1)
            lgparm = unsafe_wrap(Array,cl01+sizeof(Int32)*42,1)
            lcpgrd = unsafe_wrap(Array,cl01+sizeof(Int32)*43,1)
            lblgrd = unsafe_wrap(Array,cl01+sizeof(Int32)*44,1)
            lblsym = unsafe_wrap(Array,cl01+sizeof(Int32)*45,1)
            lplot = unsafe_wrap(Array,cl01+sizeof(Int32)*46,1)
            lsym = unsafe_wrap(Array,cl01+sizeof(Int32)*47,1)
            liqset = unsafe_wrap(Array,cl01+sizeof(Int32)*48,1)
            lclip = unsafe_wrap(Array,cl01+sizeof(Int32)*49,1)
            lvlab = unsafe_wrap(Array,cl01+sizeof(Int32)*50,1)
            lcurs = unsafe_wrap(Array,cl01+sizeof(Int32)*51,1)
            lland = unsafe_wrap(Array,cl01+sizeof(Int32)*52,1)
            lpgrid = unsafe_wrap(Array,cl01+sizeof(Int32)*53,1)
            lpcdw = unsafe_wrap(Array,cl01+sizeof(Int32)*54,1)
            lplist = unsafe_wrap(Array,cl01+sizeof(Int32)*55,1)
            lplegn = unsafe_wrap(Array,cl01+sizeof(Int32)*56,1)
            lcminp = unsafe_wrap(Array,cl01+sizeof(Int32)*57,1)
            lhmomp = unsafe_wrap(Array,cl01+sizeof(Int32)*58,1)
            lexitflag = unsafe_wrap(Array,cl01+sizeof(Int32)*59,1)

            cl02 = cglobal((:cl02_, $(library)),Int32)
            tforce = unsafe_wrap(Array,cl02,ISX)

            cc01 = cglobal((:cc01_, $(library)),UInt8)
            fname = unsafe_string(cc01+sizeof(UInt8)*0,64)
            name = unsafe_string(cc01+sizeof(UInt8)*64,48)
            ispars = unsafe_string(cc01+sizeof(UInt8)*(64+48),80)
            oname = unsafe_string(cc01+sizeof(UInt8)*(64+48+80),64)
            prefix = unsafe_string(cc01+sizeof(UInt8)*(2*64+48+80),64)
            pfname = unsafe_string(cc01+sizeof(UInt8)*(3*64+48+80),64)
            pfnamx = unsafe_string(cc01+sizeof(UInt8)*(4*64+48+80),64)
            namepol = unsafe_string(cc01+sizeof(UInt8)*(5*64+48+80),48)
            nameref = unsafe_string(cc01+sizeof(UInt8)*(5*64+2*48+80),48)

            cc02 = cglobal((:cc02_, $(library)),UInt8)
            labref = unsafe_string(cc02,32)

            cc03 = cglobal((:cc03_, $(library)),UInt8)
            vmxbl = unsafe_string(cc03,1)

            return XfoilGlobals{$T}(
            version,
            aij, dij,
            qinv, qvis, cpi, cpv, qinvu, qinv_a,
            x, y, xp, yp, s, sle, xle, yle, xte, yte, chord, yimage, wgap, waklen,
            gam, gamu, gam_a, sig, nx, ny, apanel, sst, sst_go, sst_gp, gamte, gamte_a,
            sigte, sigte_a, dste, ante, aste,
            ssple, sspec, xspoc, yspoc, qgamm, qspec, qspecp, algam, clgam, cmgam, alqsp,
            clqsp, cmqsp, qf0, qf1, qf2, qf3, qdof0, qdof1, qdof2, qdof3, clspec, ffilt,
            adeg, alfa, awake, mvisc, avisc, xcmref, ycmref, cl, cm, cd, cdp, cdf,
            cl_alf, cl_msq, psio, circ, cosa, sina, qinf, gamma, gamm1, minf1, minf,
            minf_cl, tklam, tkl_msq, cpstar, qstar, cpmn, cpmni, cpmnv, xcpmni, xcpmnv,
            xpref, cpref, verspol, cpolxy, machp1, reynp1, acritp, xstripp,
            pi, hopi, qopi, dtor,
            cvpar, cterat, ctrrat, xsref1, xsref2, xpref1, xpref2,
            size, scrnfr, plotar, pfac, qfac, vfac, xwind, ywind, xpage, ypage, xmarg,
            ymarg, ch, chg, chq, xofair, yofair, facair, xofa, yofa, faca, uprwt, cpmin,
            cpmax, cpdel, cpolplf, xcdwid, xalwid, xocwid,
            xb, yb, xbp, ybp, sb, snew, xbf, ybf, xof, yof, hmom, hfx, hfy, xbmin,
            xbmax, ybmin, ybmax, sble, chordb, areab, radble, angbte, ei11ba, ei22ba,
            apx1ba, apx2ba, ei11bt, ei22bt, apx1bt, apx2bt, thickb, cambrb, xcm, ycm,
            scm, xcmp, ycmp, xtk, ytk, stk, xtkp, ytkp,
            xssi, uedg, uinv, mass, thet, dstr, ctau, delt, uslp, guxq, guxd, tau, dis,
            ctq, vti, reinf1, reinf, reinf_cl, acrit, xstrip, xoctr, yoctr, xssitr,
            uinv_a,
            rmsbl, rmxbl, rlx, vaccel,
            xsf, ysf, xoff, yoff, xgmin, xgmax, ygmin, ygmax, dxyg, xcmin, xcmax, ycmin,
            ycmax, dxyc, dyoffc, xpmin, xpmax, ypmin, ypmax, dxyp, dyoffp, ysfp, gtick,
            q, dq, dzdg, dzdn, dzdm, dqdg, dqdm, qtan1, qtan2, z_qinf, z_alfa, z_qdof0,
            z_qdof1, z_qdof2, z_qdof3,
            va, vb, vdel, vm, vz,
            iq1, iq2, nsp, nqsp, kqtarg, iacqsp, nc1, nname, nprefix,
            ncpref, napol, npol, ipact, nlref, icolp, icolr, imatyp, iretyp, nxypol,
            npolref, ndref,
            n, nb, nw, npan, ist, kimage, itmax, nseqex, retyp, matyp, aijpix, idev,
            idevrp, ipslu, ncolor, icols, nover, ncm, ntk,
            iblte, nbl, ipan, isys, nsys, itran,
            imxbl, ismxbl,
            ok, limage, sharp, lgamu, lqinu, lvisc, lalfa, lwake, lpacc, lblini, lipan, lqaij,
            ladij, lwdij, lcpxx, lqvdes, lqrefl, lqspec, lvconv, lcpref, lclock, lpfile, lpfilx,
            lppsho, lbflap, lflap, leiw, lscini, lforef, lnorm, lgsame, lplcam, lqsym, lgsym,
            lqgrid, lggrid, lgtick, lqslop, lgslop, lcslop, lqsppl, lgeopl, lgparm, lcpgrd,
            lblgrd, lblsym, lplot, lsym, liqset, lclip, lvlab, lcurs, lland, lpgrid, lpcdw,
            lplist, lplegn, lcminp, lhmomp, lexitflag,
            tforce,
            fname, name, ispars, oname, prefix, pfname, pfnamx, namepol, nameref,
            labref,
            vmxbl)
        end

    end

end


text:
### Length

@sjdoc Length """
    Length(expr)

print the length of `expr`.

For Symata expressions, the length is the number or arguments. For scalar Julia
types, the length is zero. For `Array`'s and `Dict`'s the length is the same as
Julia `length`.
"""
@mkapprule Length :nargs => 1
@doap Length(x) = symlength(x)

symlength(mx::Mxpr) = length(mx)
symlength(s::String) = 1
symlength(s::Symbol) = 0
symlength(s::Number) = 0
function symlength(x)
    try
        length(x)
    catch
        0
    end
end

### LeafCount

## This is not the same as Mma.
## 1/2 + I  -> Complex[Rationa[1,2],1]  : LeafCount is 5
## We have Complex{:Rational}(Rational(1,2),Rational(1,1)), LeafCount is 7

"""
    leaf_count(x)

return the number of nodes in `x`. This views `Mxpr`'s as trees and all other objects as nodes.
This is `1` if `x` is a node and the number of nodes in the tree if `x` is an `Mxpr`.
An empty `Mxpr` is an empty tree, i.e. a node.

`LeafCount` is Mma's term. But better might be `NodeCount`,
because it counts all nodes in the tree, not only terminal nodes.
"""
leaf_count(x) = 1
leaf_count(x::Complex) = 3
leaf_count(x::Complex{T}) where {T<:Rational} = 7
leaf_count(x::Rational) = 3
leaf_count(mx::Mxpr) = sum(leaf_count,margs(mx)) + 1  #   1 for the Head

### ByteCount

## Try to count bytes allocated for everything in an object.
## This relies in part on Julia's ability to do this with sizeof.
## But the result is not quite what we want, so we define jssizeof.
## For instance, sizeof a BigInt always returns 16. The amount of data allocated
## is it is probably actually  8 * x.size or x.alloc.
## We also guess what to do for symbol, and we don't yet handle BigFloat.

jssizeof(x) = sizeof(x)
# I think they are 64 bit chunks
jssizeof(x::T) where {T<:BigInt} = 8 * x.alloc
jssizeof(a::Symbol) = Int(ccall(:strlen, Int32, (Ptr{UInt8},), a))
byte_count(x) = jssizeof(x)
function byte_count(mx::Mxpr)
    count::Int = sum(byte_count, margs(mx))
    count += jssizeof(mx)
    return count
end

### Depth

## This is the maximum depth of the tree.
## We do not descend into Complex and Rational.
## Different than tree for LeafCount.
depth(x) = 1
function depth(mx::Mxpr)
    d::Int = 1
    for i in 1:length(mx)
        if is_Mxpr(mx[i])
            nd = depth(mx[i])
            if nd > d
                d = nd
            end
        end
    end
    return d + 1
end

### LeafCount

@sjdoc LeafCount """
    LeafCount(expr)

gives the number of indivisible (`Part` can't be taken) elements in `expr`.

This amounts to counting all the `Heads` and all of the arguments that are not of type `Mxpr`,
that is compound expressions.

A more accurate name is `NodeCount`.
"""
@mkapprule LeafCount :nargs => 1
@doap LeafCount(x) = leaf_count(x)

#### ByteCount

@sjdoc ByteCount """
    ByteCount(expr)

gives number of bytes in `expr`.
"""
@mkapprule ByteCount :nargs => 1
@doap ByteCount(x) = byte_count(x)

#### Depth

@sjdoc Depth """
    Depth(expr)

gives the maximum number of indices required to specify
any part of `expr`, plus `1`.
"""
@mkapprule Depth :nargs => 1
@doap Depth(x) = depth(x)

### Dimensions

@sjdoc Dimensions """
    Dimensions(expr)

returns a list of the dimensions of `expr`.
"""

@mkapprule Dimensions :nargs => 1:2
@doap Dimensions(x) = 0

mutable struct DimensionsData
    dims::Array{Int,1}
    level::Int
    head
    done::Bool
end

@doap function Dimensions(x::Mxpr)
    data = DimensionsData(Array{Int}(undef, 0), 1, mhead(x), false)
    dimensions(x,data)
    MList(data.dims)
end

dimensions(x,data) = 0

function dimensions(x::Mxpr,data)
    length(x) ==  0 && return 0
    if mhead(x[1]) !=  data.head
        data.done = true
    else
        thelength = symlength(x[1])
        for i in 2:length(x)
            if symlength(x[i]) != thelength || mhead(x[i]) != data.head
                data.done = true
                break
            end
        end
    end
    if length(data.dims) < data.level
        push!(data.dims, symlength(x))
    end
    data.done && return
    data.level +=1
    for i in 1:symlength(x)
        dimensions(x[i],data)
        data.done && return
    end
    data.level -=1
end

## Maybe Allocated belongs with Time and Timing. They are all in different places.

### Allocated

@sjdoc Allocated """
    Allocated(expr)

evaluate `expr` and return a list of the memory allocated
and the result of the evaluation.
"""

function apprules(mxt::Mxpr{:Allocated})
    local mx
    a = @allocated begin
        reset_meval_count()
        mx = doeval(mxt[1])
        setsymval(:ans,mx)  ## why here ?
    end
    mxpr(:List,a,mx)
end


text:
#=
The decimal number, 585 = 1001001001 (binary), is palindromic in both bases.

Find the sum of all numbers, less than one million, which are palindromic in base 10 and base 2.

(Please note that the palindromic number, in either base, may not include leading zeros.)
=#

function isdecpalindrome(a)
  digits(a) == reverse(digits(a))
end

function isbinarypalindrome(a)
  digits(a, base=2) == reverse(digits(a, base=2))
end

function calc()
  total = 0
  for i in 1:999_999
    if isdecpalindrome(i) && isbinarypalindrome(i)
      total += i
    end
  end
  total
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
using DataFrames
using CSV

#Function that reads in counts matrix
function ReadCountsMatrix(filename)
    matrix = CSV.read(filename, header = true, delim = '\t')
    return matrix
end

#Function that reads in gene-isoform table
function ReadGeneIsoformRelationships(filename)
    relationships = CSV.read(filename, header = false, delim = ' ')
    return relationships
end

#Function that checks gene-isoform table matches counts matrix
function CheckGenesMatchIsoforms(countsMatrix, GeneIsoformRelationship)
    arr1 = sort(countsMatrix[:Column1])
    arr2 = sort(GeneIsoformRelationship[:Column2])
    if arr1 == arr2
        return true
    end
    return false
end

#Function that adds gene information to isoform table
function AddGeneInfo(countsMatrix, GeneIsoformRelationship)

    #Call the columns something sensible
    rename!(countsMatrix, :Column1 => :Isoforms)
    rename!(GeneIsoformRelationship, :Column1 => :Genes)
    rename!(GeneIsoformRelationship, :Column2 => :Isoforms)

    return join(GeneIsoformRelationship, countsMatrix, on = :Isoforms)
end

#Function that converts dataframe to dict, removes unexpressed isoforms and genes,
#and returns a dict of how many isoforms are expressed per gene.
function KeepGenesWithNIsoforms(countsMatrix, N)

    #Check N is set to a sensible value
    if N<1
        throw(ArgumentError("You can't simulate the expression of less than one
        isoform!\n"))
    end

    if N>4
        println("You have asked to simulate more than 4 isoforms. Not many genes
        have more than 4 isoforms, so this might violate the assumptions of some
         of iPredict's statistical models.")
    end

    #sort countsMatrix. The rationale behind this is that in theory we can make
    #this closer to O(n) instead of O(n^2)
    sort!(countsMatrix, (:Genes))

    #Set first gene
    gene = ""

    #Make empty dicts
    local dict_to_add = Dict{String, Array{Float64, 1}}()
    local output_dict = Dict{String, Any}()
    local unexpr_dict = Dict{String, Int64}()
    local NumIsos = 0
    local MaxNumIsos = 0
    local UnexprIsos = 0

    #Iterate once over countsMatrix
    for i in 1:nrow(countsMatrix)

        row = countsMatrix[i,:]

        #Check gene is still the same
        if row.Genes == gene

            #if N > numisos or N < numalreadyexpressed, continue
            if N > NumIsos
                continue
            end
            if N < length(dict_to_add)
                continue
            end

            #See whether row passes cutoff
            cutoff = sum(row[3:ncol(countsMatrix)].>=5)

            if cutoff >= 2
                #Add isoform to dictionary
                dict_to_add[row.Isoforms] = Vector(row[3:ncol(countsMatrix)])
            end

        else #new gene

            if length(dict_to_add) == N #See whether the previous gene passes
                output_dict[gene] = dict_to_add
                unexpr_dict[gene] = NumIsos - N

                if NumIsos > MaxNumIsos
                    MaxNumIsos = NumIsos
                end
            end

            #Reset gene, NumExpressed and dict_to_add
            gene = countsMatrix.Genes[i]
            dict_to_add = Dict{String, Array{Float64, 1}}()

            #Find NumIsos
            NumIsos = 0
            j = deepcopy(i)
            while countsMatrix.Genes[j] == gene
                j += 1
                NumIsos += 1
                #Make sure we don't go beyond the end of the matrix
                if j > nrow(countsMatrix)
                    break
                end
            end

            if NumIsos < N
                continue
            end

            #See whether row passes cutoff
            cutoff = sum(row[3:ncol(countsMatrix)].>=5)

            if cutoff >= 2
                #Add isoform to dictionary
                dict_to_add[row.Isoforms] = Vector(row[3:ncol(countsMatrix)])
            end
        end

        #Catch last gene in table
        if i == nrow(countsMatrix)
            if length(dict_to_add) == N
                output_dict[gene] = dict_to_add
            end
        end
    end

    #Return dictionary
    return (output_dict, unexpr_dict, MaxNumIsos)

end

#Function that converts expressed isoforms for a single gene to a dict. Used by
#oneGenePrediction()
function ConvertExpressedIsosToDict(countsMatrix)

    #Make empty dicts
    local dict_to_add = Dict{String, Array{Float64, 1}}()
    local output_dict = Dict{String, Any}()

    for row in eachrow(countsMatrix)

        #See whether row passes cutoff
        cutoff = sum(row[3:ncol(countsMatrix)].>=5)

        if cutoff >= 2
            #Add isoform to dictionary
            dict_to_add[row.Isoforms] = Vector(row[3:ncol(countsMatrix)])
        end
    end

    #Add gene info and return dictionary
    output_dict[countsMatrix.Genes[1]] = dict_to_add
    return output_dict
end

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles