text:
## Copyright (c) 2015, Júlio Hoffimann Mendes <juliohm@stanford.edu>
##
## Permission to use, copy, modify, and/or distribute this software for any
## purpose with or without fee is hereby granted, provided that the above
## copyright notice and this permission notice appear in all copies.
##
## THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
## WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
## MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
## ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
## WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
## ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
## OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

function dijkstra_cut(A1::AbstractArray, A2::AbstractArray, dir::Symbol)
  # permute dimensions so that the algorithm is
  # the same for cuts in x, y and z directions
  B = abs.(A1 - A2)
  if dir == :x
    B = permutedims(B, [1,2,3])
  elseif dir == :y
    B = permutedims(B, [2,1,3])
  elseif dir == :z
    B = permutedims(B, [3,2,1])
  end

  mx, my, mz = size(B)

  # accumulation cube and overlap mask
  E = zeros(B); M = falses(B)

  # pad accumulation cube with +inf
  Epad(i,j,k) = all(0 .< [i,j,k] .≤ [mx,my,mz]) ? E[i,j,k] : Inf

  # forward accumulation along 3D cube
  E[:,:,1] = B[:,:,1]
  for k=2:mz, i=1:mx, j=1:my
    square = [Epad(i-1 , j-1 , k-1) , Epad(i-1 , j , k-1) , Epad(i-1 , j+1 , k-1),
              Epad(i   , j-1 , k-1) , Epad(i   , j , k-1) , Epad(i   , j+1 , k-1),
              Epad(i+1 , j-1 , k-1) , Epad(i+1 , j , k-1) , Epad(i+1 , j+1 , k-1)]
    E[i,j,k] = B[i,j,k] + minimum(square)
  end

  # forward accumulation along 2D square (last slice of the cube)
  zslice = view(E,:,:,mz)
  for j=2:my
    zslice[1,j] += minimum(zslice[1:2,j-1])
    for i=2:mx-1
      zslice[i,j] += minimum(zslice[i-1:i+1,j-1])
    end
    zslice[mx,j] += minimum(zslice[mx-1:mx,j-1])
  end

  # backward search along last slice
  idx = indmin(zslice[:,my])
  mslice = view(M,:,:,mz)
  mslice[1:idx,my] = trues(idx)
  idxvec = zeros(Int, my); idxvec[my] = idx # keep track of indexes
  width = isodd(mx) ? mx+1 : mx # avoid zig-zag artifact
  for j=my-1:-1:1
    for i=1:width
      if idx < mx && minimum(zslice[max(idx-1,1):idx+1,j]) == zslice[idx+1,j]
        idx += 1
      elseif idx > 1 && zslice[idx-1,j] ≤ zslice[idx,j]
        idx -= 1
      end
    end
    mslice[1:idx,j] = true
    idxvec[j] = idx
  end

  # backward search along cube
  for j=1:my
    yslice = view(E,:,j,:)
    idx = idxvec[j]
    for k=mz-1:-1:1
      for i=1:width
        if idx < mx && minimum(yslice[max(idx-1,1):idx+1,k]) == yslice[idx+1,k]
          idx += 1
        elseif idx > 1 && yslice[idx-1,k] ≤ yslice[idx,k]
          idx -= 1
        end
      end
      M[1:idx,j,k] = true
    end
  end

  # permute back to original shape
  if dir == :x
    M = permutedims(M, [1,2,3])
  elseif dir == :y
    M = permutedims(M, [2,1,3])
  elseif dir == :z
    M = permutedims(M, [3,2,1])
  end

  M
end


text:
#####
##### Nothing happens when your boundary condition is nothing
#####

fill_west_halo!(c, ::Nothing, arch, barrier, grid, args...) = nothing
fill_east_halo!(c, ::Nothing, arch, barrier, grid, args...) = nothing

fill_north_halo!(c, ::Nothing, arch, barrier, grid, args...) = nothing
fill_south_halo!(c, ::Nothing, arch, barrier, grid, args...) = nothing

fill_top_halo!(c, ::Nothing, arch, barrier, grid, args...) = nothing
fill_bottom_halo!(c, ::Nothing, arch, barrier, grid, args...) = nothing

ext: jl


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

"""
    AbstractField

Abstract supertype for all fields in JuliaFEM.
"""
abstract type AbstractField end

function length(f::F) where F<:AbstractField
    return length(f.data)
end

function size(f::F) where F<:AbstractField
    return size(f.data)
end

function ==(x::F, y) where F<:AbstractField
    return ==(x.data, y)
end

function ==(x, y::F) where F<:AbstractField
    return ==(x, y.data)
end

function ==(x::F, y::F) where F<:AbstractField
    return ==(x.data, y.data)
end

function getindex(f::F, i::Int) where F<:AbstractField
    return getindex(f.data, i)
end

function interpolate_field(field::AbstractField, ::Any)
    return field.data
end

function update_field!(field::AbstractField, data)
   field.data = data
end

"""
    DCTI{T} <: AbstractField

Discrete, constant, time-invariant field.

This field is constant in both spatial direction and time direction,
i.e. df/dX = 0 and df/dt = 0.

# Example

```jldoctest
julia> DCTI(1)
FEMBase.DCTI{Int64}(1)
```
"""
mutable struct DCTI{T} <: AbstractField
    data :: T
end

function getindex(field::DCTI, ::Int)
    return field.data
end

"""
    DVTI{N,T} <: AbstractField

Discrete, variable, time-invariant field.

This is constant in time direction, but not in spatial direction, i.e. df/dt = 0
but df/dX != 0. The basic structure of data is `Tuple`, and it is implicitly
assumed that length of field matches to the number of shape functions, so that
interpolation in spatial direction works.

# Example

```jldoctest
julia> DVTI(1, 2, 3)
FEMBase.DVTI{3,Int64}((1, 2, 3))
```
"""
mutable struct DVTI{N,T} <: AbstractField
    data :: NTuple{N,T}
end

function DVTI(data...)
    return DVTI(data)
end

"""
    DCTV{T} <: AbstractField

Discrete, constant, time variant field. This type of field can change in time
direction but not in spatial direction.

# Example

Field having value 5 at time 0.0 and value 10 at time 1.0:

```jldoctest
julia> DCTV(0.0 => 5, 1.0 => 10)
FEMBase.DCTV{Int64}(Pair{Float64,Int64}[0.0=>5, 1.0=>10])
```

"""
mutable struct DCTV{T} <: AbstractField
    data :: Vector{Pair{Float64,T}}
end

function DCTV(data::Pair{Float64,T}...) where T
    return DCTV(collect(data))
end

function update_field!(f::DCTV, data::Pair{Float64, T}) where T
    if isapprox(last(f.data).first, data.first)
        f.data[end] = data
    else
        push!(f.data, data)
    end
end

function interpolate_field(field::DCTV, time)
    time < first(field.data).first && return first(field.data).second
    time > last(field.data).first && return last(field.data).second
    for i=reverse(1:length(field))
        isapprox(field.data[i].first, time) && return field.data[i].second
    end
    for i=length(field.data):-1:2
        t0 = field.data[i-1].first
        t1 = field.data[i].first
        if t0 < time < t1
            y0 = field.data[i-1].second
            y1 = field.data[i].second
            dy = y1-y0
            dt = t1-t0
            return y0 + (time-t0)*dy/dt
        end
    end
end

"""
    DVTV{N,T} <: AbstractField

Discrete, variable, time variant field. The most general discrete field can
change in both temporal and spatial direction.

# Example

```jldoctest
julia> DVTV(0.0 => (1, 2), 1.0 => (2, 3))
FEMBase.DVTV{2,Int64}(Pair{Float64,Tuple{Int64,Int64}}[0.0=>(1, 2), 1.0=>(2, 3)])
```
"""
mutable struct DVTV{N,T} <: AbstractField
    data :: Vector{Pair{Float64,NTuple{N,T}}}
end

function DVTV(data::Pair{Float64,NTuple{N,T}}...) where {N,T}
    return DVTV(collect(data))
end

function update_field!(f::DVTV, data::Pair{Float64, NTuple{N,T}}) where {N,T}
    if isapprox(last(f.data).first, data.first)
        f.data[end] = data
    else
        push!(f.data, data)
    end
end

function interpolate_field(field::DVTV{N,T}, time) where {N,T}
    time < first(field.data).first && return first(field.data).second
    time > last(field.data).first && return last(field.data).second
    for i=reverse(1:length(field))
        isapprox(field.data[i].first, time) && return field.data[i].second
    end
    for i=length(field.data):-1:2
        t0 = field.data[i-1].first
        t1 = field.data[i].first
        if t0 < time < t1
            y0 = field.data[i-1].second
            y1 = field.data[i].second
            dt = t1-t0
            return map((a,b) -> a + (time-t0)*(b-a)/dt, y0, y1)
        end
    end
end

"""
    CVTV <: AbstractField

Continuous, variable, time variant field.

# Example

```jldoctest
julia> f = CVTV((xi,t) -> xi*t)
FEMBase.CVTV(#1)
```
"""
mutable struct CVTV <: AbstractField
    data :: Function
end

function (f::CVTV)(xi, time)
    return f.data(xi, time)
end

"""
    DVTId(X::Dict)

Discrete, variable, time invariant dictionary field.
"""
mutable struct DVTId{T} <: AbstractField
    data :: Dict{Int, T}
end

function update_field!(field::DVTId{T}, data::Dict{Int, T}) where T
    merge!(field.data, data)
end

"""
    DVTVd(time => data::Dict)

Discrete, variable, time variant dictionary field.
"""
mutable struct DVTVd{T} <: AbstractField
    data :: Vector{Pair{Float64,Dict{Int,T}}}
end

function DVTVd(data::Pair{Float64,Dict{Int,T}}...) where T
    return DVTVd(collect(data))
end

function interpolate_field(field::DVTVd{T}, time) where T
    time >= last(field.data).first && return last(field.data).second
    time <= first(field.data).first && return first(field.data).second
    for i=reverse(1:length(field))
        isapprox(field.data[i].first, time) && return field.data[i].second
    end
    for i=length(field.data):-1:2
        t0 = field.data[i-1].first
        t1 = field.data[i].first
        if t0 < time < t1
            y0 = field.data[i-1].second
            y1 = field.data[i].second
            f = (time-t0)/(t1-t0)
            new_data = empty(y0)
            for i in keys(y0)
                new_data[i] = f*y0[i] + (1-f)*y1[i]
            end
            return new_data
        end
    end
end

function update_field!(f::DVTVd, data::Pair{Float64,Dict{Int,T}}) where T
    if isapprox(last(f.data).first, data.first)
        f.data[end] = data
    else
        push!(f.data, data)
    end
end

function new_field(data)
    return DCTI(data)
end

function new_field(data...)
    return DVTI(data)
end

function new_field(data::NTuple{N,T}) where {N,T}
    return DVTI(data)
end

function new_field(data::Pair{Float64,T}...) where T
    return DCTV(collect(data))
end

function new_field(data::Pair{Float64,NTuple{N,T}}...) where {N,T}
    return DVTV(collect(data))
end

function new_field(data::Function)
    return CVTV(data)
end

function new_field(data::Pair{Int, T}...) where T
    return DVTId(Dict(data))
end

function new_field(data::Pair{Float64, NTuple{N, Pair{Int, T}}}...) where {N,T}
    return DVTVd(collect(t => Dict(d) for (t, d) in data))
end

function new_field(data::Dict{Int,T}) where T
    return DVTId(data)
end

function new_field(data::Pair{Float64, Dict{Int, T}}...) where T
    return DVTVd(collect(data))
end

"""
    field(x)

Create new field. Field type is deduced from data type.
"""
function field(data...)
    return new_field(data...)
end

"""
    interpolate(field, time)

Interpolate field in time direction.

# Examples

For time invariant fields [`DCTI`](@ref), [`DVTI`](@ref), [`DVTId`](@ref)
solution is trivially the data inside field as fields does not depend from
the time:

```jldoctest
julia> a = field(1.0)
FEMBase.DCTI{Float64}(1.0)

julia> interpolate(a, 0.0)
1.0
```

```jldoctest
julia> a = field((1.0, 2.0))
FEMBase.DVTI{2,Float64}((1.0, 2.0))

julia> interpolate(a, 0.0)
(1.0, 2.0)
```

```jldoctest
julia> a = field(1=>1.0, 2=>2.0)
FEMBase.DVTId{Float64}(Dict(2=>2.0,1=>1.0))

julia> interpolate(a, 0.0)
Dict{Int64,Float64} with 2 entries:
  2 => 2.0
  1 => 1.0
```

DVTId trivial solution is returned. For time variant fields DCTV, DVTV, DVTVd
linear interpolation is performed.

# Other notes

First algorithm checks that is time out of range, i.e. time is smaller than
time of first frame or larger than last frame. If that is the case, return
first or last frame. Secondly algorithm finds is given time exact match to
time of some frame and return that frame. At last, we find correct bin so
that t0 < time < t1 and use linear interpolation.

"""
function interpolate(field::AbstractField, time)
    return interpolate_field(field, time)
end

"""
    interpolate(a, b)

A helper function for interpolate routines. Given iterables `a` and `b`,
calculate c = aᵢbᵢ. Length of `a` can be less than `b`, but not vice versa.
"""
function interpolate(a, b)
    @assert length(a) <= length(b)
    return sum(a[i]*b[i] for i=1:length(a))
end


text:
function svdbyeigen(A; nsv::Integer=3)
	P,N = size(A)
	K = Symmetric(N<=P ? A'A : A*A')
	M = size(K,1)
	F = eigen(K, M-nsv+1:M)
	S = sqrt.(max.(0.,reverse(F.values)))

	V = F.vectors[:,end:-1:1]
	N<=P ? SVD(A*V./S',S,V') : SVD(V,S,V'A./S)
end


text:
using HetaSimulator, Plots

model = load_jlmodel("../dist/julia/model.jl")

res = sim(model; tspan = (1,340))

plotd = plot(res; vars = [:P_p, :A_p, :O_p])
savefig(plotd, "julia-plot-1.png")

plotd = plot(res; vars = [:PFM_p])
savefig(plotd, "julia-plot-2.png")

max_stars_repo_name: insysbio/faah-inhibitor

max_issues_repo_name: insysbio/faah-inhibitor

max_forks_repo_name: insysbio/faah-inhibitor


text:
@testset "LAPACK" begin
    @testset "trsyl!" begin
        @testset "T=$T, m=$m, n=$n, transa='$transa', transb='$transb', isgn=$isgn" for T in
                                                                                (
                Float64, ComplexF64
            ),
            transa in (T <: Real ? ('N', 'C', 'T') : ('N', 'C')),
            transb in (T <: Real ? ('N', 'C', 'T') : ('N', 'C')),
            m in (2, 3),
            n in (1, 3),
            isgn in (1, -1)

            # make A and B quasi upper-triangular (or upper-triangular for complex)
            # and their tangents have the same sparsity pattern
            A = schur(randn(T, m, m)).T
            B = schur(randn(T, n, n)).T
            C = randn(T, m, n)
            test_frule(
                LAPACK.trsyl!,
                transa,
                transb,
                A ⊢ rand_tangent(A) .* (!iszero).(A),  # Match sparsity pattern
                B ⊢ rand_tangent(B) .* (!iszero).(B),
                C,
                isgn,
            )
        end
    end
end


text:
function config_solver(solver::MathOptInterface.AbstractOptimizer,
    solver_config::Dict)
    for (k, v) in solver_config
        MOI.set(solver,
            MOI.RawOptimizerAttribute(k), v)
    end
end

"""
    ordinalgwas(nullformula, covfile, geneticfile; kwargs...)
    ordinalgwas(nullformula, df, geneticfile; kwargs...)
    ordinalgwas(fittednullmodel, geneticfile; kwargs...)
    ordinalgwas(fittednullmodel, bedfile, bimfile, bedn; kwargs...)

# Positional arguments
- `nullformula::FormulaTerm`: formula for the null model.
- `covfile::AbstractString`: covariate file (csv) with one header line. One column
    should be the ordinal phenotype coded as integers starting from 1.  For example,
    ordinal phenotypes can be coded as 1, 2, 3, 4 but not 0, 1, 2, 3.
- `df::DataFrame`: DataFrame containing response and regressors for null model.
- `geneticfile::Union{Nothing, AbstractString}`: File containing genetic information for GWAS.
    This includes a PLINK file name without the .bed, .fam, or .bim
    extensions or a VCF file without the .vcf extension. If `geneticfile==nothing`,
    only null model is fitted. If `geneticfile` is provided, bed, bim, and fam file (or vcf) with
    the same `geneticfile` prefix need to exist. Compressed file formats such as gz and bz2
    are allowed. Check all allowed formats by `SnpArrays.ALLOWED_FORMAT`. If you're using a VCF file,
    make sure to use the `geneticformat = "VCF"` keyword option, and specificy dosage (:DS) or
    genotype (:GT) data with the `vcftype` command.
- `fittednullmodel::StatsModels.TableRegressionModel`: the fitted null model
    output from `ordinalgwas(nullformula, covfile)` or `ordinalgwas(nullformula, df)`.
- `bedfile::Union{AbstractString,IOStream}`: path to Plink bed file with full file name.
- `bimfile::Union{AbstractString,IOStream}`: path to Plink bim file with full file name.
- `bedn::Integer`: number of samples in bed/vcf file.

# Keyword arguments
- `analysistype`::AbstractString: Type of analysis to conduct. Default is `singlesnp`. Other options are `snpset` and `gxe`.
- `geneticformat`::AbstractString: Type of file used for the genetic analysis. `"PLINK"`, `"VCF"`, and `"BGEN"` are currently supported. Default is PLINK.
- `vcftype`::Union{Symbol, Nothing}: Data to extract from the VCF file for the GWAS analysis. `:DS` for dosage or `:GT` for genotypes. Default is nothing.
- `nullfile::Union{AbstractString, IOStream}`: output file for the fitted null model;
    default is `ordinalgwas.null.txt`.
- `pvalfile::Union{AbstractString, IOStream}`: output file for the gwas p-values; default is
    `ordinalgwas.pval.txt`.
- `covtype::Vector{DataType}`: type information for `covfile`. This is useful
    when `CSV.read(covarfile)` has parsing errors.
- `covrowinds::Union{Nothing,AbstractVector{<:Integer}}`: sample indices for covariate file.
- `testformula::FormulaTerm`: formula for test unit. Default is `@formula(trait ~ 0 + snp)`.
- `test::Symbol`: `:score` (default) or `:lrt`.
- `link::GLM.Link`: `LogitLink()` (default), `ProbitLink()`, `CauchitLink()`,
    or `CloglogLink()`.
- `snpmodel`: `ADDITIVE_MODEL` (default), `DOMINANT_MODEL`, or `RECESSIVE_MODEL`.
- `snpinds::Union{Nothing,AbstractVector{<:Integer}}`: SNP indices for bed/vcf file.
- `geneticrowinds::Union{Nothing,AbstractVector{<:Integer}}`: sample indices for bed/vcf file.
- `solver`: an optimizer supported by MathOptInterface. Default is
    `NLopt.Optimizer()` with `algorithm=:LD_SLSQP` and `maxeval=4000`. Another common choice is
    `Ipopt.Optimizer()`.
- `verbose::Bool`: default is `false`.
- `snpset::Union{Nothing, Integer, AbstractString, AbstractVector{<:Integer}}`: Only include
    if you are conducting a snpset analysis. An integer indicates a window of SNPs
    (i.e. every 500 snps). An abstract string allows you to specify an input file,
    with no header and two columns separated by a space. The first column must contain the snpset ID
    and the second column must contain the snpid's identical to the bimfile. An AbstractVector
    allows you to specify the snps you want to perform one joint snpset test for.
- `e::Union{AbstractString,Symbol}`: Only include if you are conducting a GxE analysis.
    Enviromental variable to be used to test the GxE interaction.
    For instance, for testing `sex & snp` interaction, use `:sex` or `"sex"`.

# Examples
The following is an example of basic GWAS with PLINK files:
```julia
plinkfile = "plinkexample"
covfile = "covexample"
ordinalgwas(@formula(trait ~ sex), covfile, plkfile)
```

The following is an example of basic GWAS with a VCF file using dosages then genotypes:
```julia
vcffile = "vcfexample"
covfile = "covexample"
ordinalgwas(@formula(trait ~ sex), covfile, vcffile;
    geneticfile = "VCF", vcftype = :DS)

ordinalgwas(@formula(trait ~ sex), covfile, vcffile;
    geneticfile = "VCF", vcftype = :GT)
```

The following is an example of snpset GWAS (every 50 snps). For more types of snpset analyses see documentation:
```julia
ordinalgwas(@formula(trait ~ sex), covfile, plkfile;
    analysistype = "snpset", snpset = 50)
```

The following is an example of GxE GWAS testing the interaction effect:
```julia
ordinalgwas(@formula(trait ~ sex), covfile, plkfile;
    analysistype = "gxe", e = :sex)
```
"""
function ordinalgwas(
    # positional arguments
    nullformula::FormulaTerm,
    covfile::AbstractString,
    geneticfile::Union{Nothing, AbstractString} = nothing;
    # keyword arguments
    covtype::Union{Nothing, Vector{DataType}} = nothing,
    covrowinds::Union{Nothing, AbstractVector{<:Integer}} = nothing,
    kwargs...
    )
    covdf = SnpArrays.makestream(covfile) do io
        CSV.read(io, DataFrame; types=covtype)
    end
    ordinalgwas(nullformula, covrowinds === nothing ? covdf : covdf[covrowinds, :],
        geneticfile; kwargs...)
end

function ordinalgwas(
    nullformula::FormulaTerm,
    nulldf::DataFrame,
    geneticfile::Union{Nothing, AbstractString} = nothing;
    nullfile::Union{AbstractString, IOStream} = "ordinalgwas.null.txt",
    link::GLM.Link = LogitLink(),
    solver::MOI.AbstractOptimizer = NLopt.Optimizer(),
    solver_config = Dict("algorithm" => :LD_SLSQP, "max_iter" => 4000),
    verbose::Bool = false,
    kwargs...
    )
    # fit and output null model
    nm = polr(nullformula, nulldf, link, solver)
    verbose && show(nm)
    SnpArrays.makestream(nullfile, "w") do io
        show(io, nm)
    end
    geneticfile === nothing && (return nm)
    ordinalgwas(nm, geneticfile; solver=solver,
        solver_config=solver_config, verbose=verbose, kwargs...)
end

function ordinalgwas(
    # positional arguments
    fittednullmodel::StatsModels.TableRegressionModel,
    geneticfile::AbstractString;
    # keyword arguments
    analysistype::AbstractString = "singlesnp",
    geneticformat::AbstractString = "PLINK",
    vcftype::Union{Symbol, Nothing} = nothing,
    testformula::FormulaTerm = fittednullmodel.mf.f.lhs ~ Term(:snp),
    test::Symbol = :score,
    pvalfile::Union{AbstractString, IOStream} = "ordinalgwas.pval.txt",
    snpmodel::Union{Val{1}, Val{2}, Val{3}} = ADDITIVE_MODEL,
    snpinds::Union{Nothing, AbstractVector{<:Integer}} = nothing,
    geneticrowinds::Union{Nothing, AbstractVector{<:Integer}} = nothing,
    solver::MOI.AbstractOptimizer = NLopt.Optimizer(),
    solver_config = Dict("algorithm" => :LD_SLSQP, "max_iter" => 4000),
    verbose::Bool = false,
    snpset::Union{Nothing, Integer, AbstractString, #for snpset analysis
        AbstractVector{<:Integer}} = nothing,
    e::Union{Nothing, AbstractString, Symbol} = nothing # for GxE analysis
    )
    # locate plink bed, fam, bim files or VCF file
    lowercase(geneticformat) in ["plink", "vcf", "bgen"] || error("`geneticformat` $geneticformat not valid. Please use 'VCF' or 'PLINK'.")
    isplink = "plink" == lowercase(geneticformat)
    if lowercase(geneticformat) == "plink"
        if isfile(geneticfile * ".bed")
            bedfile = geneticfile * ".bed"
        else
            fmt = findfirst(isfile, geneticfile * ".bed." .* SnpArrays.ALLOWED_FORMAT)
            fmt === nothing && throw(ArgumentError("bed file not found"))
            bedfile = geneticfile * ".bed." * SnpArrays.ALLOWED_FORMAT[fmt]
        end
        famfile = replace(bedfile, ".bed" => ".fam")
        isfile(famfile) || throw(ArgumentError("fam file not found"))
        bimfile = replace(bedfile, ".bed" => ".bim")
        isfile(bimfile) || throw(ArgumentError("bim file not found"))
        # selected rows should match nobs in null model
        bedn = SnpArrays.makestream(countlines, famfile)
    elseif lowercase(geneticformat) == "vcf"
        vcftype in [:GT, :DS] || throw(ArgumentError("vcftype not specified. Allowable types are :GT for genotypes and :DS for dosages."))
        if isfile(geneticfile * ".vcf")
            vcffile = geneticfile * ".vcf"
        else
            fmt = findfirst(isfile, geneticfile * ".vcf." .* SnpArrays.ALLOWED_FORMAT)
            fmt === nothing && throw(ArgumentError("VCF file not found"))
            vcffile = geneticfile * ".vcf." * SnpArrays.ALLOWED_FORMAT[fmt]
        end
        bedn = VCFTools.nsamples(vcffile)
    elseif lowercase(geneticformat) == "bgen"
        if isfile(geneticfile * ".bgen")
            bgenfile = geneticfile * ".bgen"
        else
            fmt = findfirst(isfile, geneticfile * ".bgen." .* SnpArrays.ALLOWED_FORMAT)
            fmt === nothing && throw(ArgumentError("BGEN file not found"))
            bgenfile = geneticfile * ".bgen." * SnpArrays.ALLOWED_FORMAT[fmt]
        end
        b = Bgen(bgenfile)
        bedn = n_samples(b)
    end
    if geneticrowinds === nothing
        nbedrows = bedn
        rowinds = 1:bedn
    else
        nbedrows = eltype(geneticrowinds) == Bool ? count(geneticrowinds) : length(geneticrowinds)
        rowinds = geneticrowinds
    end

    nbedrows == nobs(fittednullmodel) ||
        throw(ArgumentError("number of samples in geneticrowinds does not match null model"))

    # validate testing method
    test = Symbol(lowercase(string(test)))
    test == :score || test == :lrt || throw(ArgumentError("unrecognized test $test"))

    # gwas
    if lowercase(geneticformat) == "plink" #plink
        ordinalgwas(fittednullmodel, bedfile, bimfile, bedn;
            analysistype = analysistype,
            testformula = testformula,
            test = test,
            pvalfile = pvalfile,
            snpmodel = snpmodel,
            snpinds = snpinds,
            bedrowinds = rowinds,
            solver = solver,
            solver_config = solver_config,
            verbose = verbose,
            snpset = snpset,
            e = e)
    elseif lowercase(geneticformat) == "vcf" #vcf
        ordinalgwas(fittednullmodel, vcffile, bedn, vcftype;
            analysistype = analysistype,
            testformula = testformula,
            test = test,
            pvalfile = pvalfile,
            snpmodel = snpmodel,
            snpinds = snpinds,
            vcfrowinds = rowinds,
            solver = solver,
            solver_config = solver_config,
            verbose = verbose,
            snpset = snpset,
            e = e)
    else #bgen
        ordinalgwas(fittednullmodel, bgenfile, bedn;
            analysistype = analysistype,
            testformula = testformula,
            test = test,
            pvalfile = pvalfile,
            snpmodel = snpmodel,
            snpinds = snpinds,
            bgenrowinds = rowinds,
            solver = solver,
            solver_config = solver_config,
            verbose = verbose,
            snpset = snpset,
            e = e)
    end
end

# For PLINK Analysis
function ordinalgwas(
    fittednullmodel::StatsModels.TableRegressionModel,
    bedfile::Union{AbstractString, IOStream}, # full path and bed file name
    bimfile::Union{AbstractString, IOStream}, # full path and bim file name
    bedn::Integer;           # number of samples in bed file
    analysistype::AbstractString = "singlesnp",
    testformula::FormulaTerm = fittednullmodel.mf.f.lhs ~ Term(:snp),
    test::Symbol = :score,
    pvalfile::Union{AbstractString, IOStream} = "ordinalgwas.pval.txt",
    snpmodel::Union{Val{1}, Val{2}, Val{3}} = ADDITIVE_MODEL,
    snpinds::Union{Nothing, AbstractVector{<:Integer}} = nothing,
    bedrowinds::AbstractVector{<:Integer} = 1:bedn, # row indices for SnpArray
    solver::MOI.AbstractOptimizer = NLopt.Optimizer(),
    solver_config = Dict("algorithm" => :LD_SLSQP, "max_iter" => 4000),
    verbose::Bool = false,
    snpset::Union{Nothing, Integer, AbstractString, #for snpset analysis
        AbstractVector{<:Integer}} = nothing,
    e::Union{Nothing, AbstractString, Symbol} = nothing # for GxE analysis
    )
    config_solver(solver, solver_config)
    # create SnpArray
    genomat = SnpArrays.SnpArray(bedfile, bedn)
    cc = SnpArrays.counts(genomat, dims=1) # column counts of genomat
    mafs = SnpArrays.maf(genomat)

    # create SNP mask vector
    if snpinds === nothing
        snpmask = trues(SnpArrays.makestream(countlines, bimfile))
    elseif eltype(snpinds) == Bool
        snpmask = snpinds
    else
        snpmask = falses(SnpArrays.makestream(countlines, bimfile))
        snpmask[snpinds] .= true
    end

    analysistype = lowercase(analysistype)
    analysistype in ["singlesnp", "snpset", "gxe"] || error("Analysis type $analysis invalid option.
    Available options are 'singlesnp', 'snpset' and 'gxe'.")

    # determine analysis type
    if analysistype == "singlesnp"

        # carry out score or LRT test SNP by SNP
        snponly = testformula.rhs == Term(:snp)

        # extra columns in design matrix to be tested
        testdf = DataFrame(fittednullmodel.mf.data) # TODO: not type stable here
        testdf[!, :snp] = zeros(size(fittednullmodel.mm, 1))
        Z = similar(modelmatrix(testformula, testdf))
        SnpArrays.makestream(pvalfile, "w") do io
            if test == :score
                println(io, "chr,pos,snpid,maf,hwepval,pval")
                ts = OrdinalMultinomialScoreTest(fittednullmodel.model, Z)
            else # lrt
                nulldev = deviance(fittednullmodel.model)
                Xaug = [fittednullmodel.model.X Z]
                q = size(Z, 2)
                γ̂ = Vector{Float64}(undef, q) # effect size for columns being tested
                if snponly
                    println(io, "chr,pos,snpid,maf,hwepval,effect,stder,pval")
                else
                    print(io, "chr,pos,snpid,maf,hwepval,")
                    for j in 1:q
                        print(io, "effect$j,")
                    end
                    println(io, "pval")
                end
            end
            SnpArrays.makestream(bimfile) do bimio
                for j in eachindex(snpmask)
                    row = readline(bimio)
                    snpmask[j] || continue
                    snpj = split(row)
                    hwepval = SnpArrays.hwe(cc[1, j], cc[3, j], cc[4, j])
                    maf = mafs[j]
                    if test == :score
                        if maf == 0 # mono-allelic
                            pval = 1.0
                        else
                            if snponly
                                copyto!(ts.Z, @view(genomat[bedrowinds, j]),
                                impute=true, model=snpmodel)
                            else # snp + other terms
                                copyto!(testdf[!, :snp], @view(genomat[bedrowinds, j]),
                                    impute = true, model = snpmodel)
                                ts.Z[:] = modelmatrix(testformula, testdf)
                            end
                            pval = polrtest(ts)
                        end
                        println(io, "$(snpj[1]),$(snpj[4]),$(snpj[2]),",
                            "$maf,$hwepval,$pval")
                    elseif test == :lrt
                        if maf == 0 # mono-allelic
                            fill!(γ̂, 0)
                            stderr = -1.0
                            pval = 1.0
                        else
                            if snponly
                                copyto!(@view(Xaug[:, fittednullmodel.model.p+1]),
                                    @view(genomat[bedrowinds, j]),
                                    impute=true, model=snpmodel)
                            else # snp + other terms
                                copyto!(testdf[!, :snp], @view(genomat[bedrowinds, j]),
                                    impute=true, model=snpmodel)
                                Xaug[:, fittednullmodel.model.p+1:end] = modelmatrix(testformula, testdf)
                            end
                            altmodel = polr(Xaug, fittednullmodel.model.Y,
                                fittednullmodel.model.link, solver,
                                wts = fittednullmodel.model.wts)
                            copyto!(γ̂, 1, altmodel.β, fittednullmodel.model.p + 1, q)
                            pval = ccdf(Chisq(q), nulldev - deviance(altmodel))
                            stderr = stderror(altmodel)[fittednullmodel.model.npar + 1]
                        end
                        if snponly
                            println(io, "$(snpj[1]),$(snpj[4]),$(snpj[2]),$maf,$hwepval,",
                                "$(γ̂[1]),$stderr,$pval")
                        else
                            print(io, "$(snpj[1]),$(snpj[4]),$(snpj[2]),$maf,$hwepval,")
                            for j in 1:q
                                print(io, "$(γ̂[j]),")
                            end
                            println(io, pval)
                        end
                    end
                end
            end
        end
    elseif analysistype == "snpset"
        # determine type of snpset analysis
        if isa(snpset, Nothing)
            @warn("Nothing set for `snpset`.
            This will default to `singlesnp` analysis (windowsize = 1).")
            setlength = 1
        elseif isa(snpset, AbstractString)
            isfile(snpset) || throw(ArgumentError("snpset file not found,
            to specify a window replace snpset string with a window size"))
            #first column SNPset ID, second column SNP ID
            snpsetFile = CSV.read(snpset, DataFrame, header = [:snpset_id, :snp_id], delim = " ")
            #make sure it matches bim file
            biminfo = CSV.read(bimfile, DataFrame, header = [:chr, :snp_id, :c3, :bp, :c5, :c6], delim = "\t")
            snpsetFile[!, :snp_id] == biminfo[!, :snp_id] || throw(ArgumentError("snp order in snpset file
            must match (in the same order) bimfile"))
            snpset_ids = unique(snpsetFile[!, :snpset_id])
            nSets = length(snpset_ids)
            setlength = 0
        elseif isa(snpset, Integer)
            setlength = snpset
        else #abstract vector (boolean of true at indicies or range or indicies)
            setlength = -1
        end

        # conduct analysis based on type
        if setlength > 0 #single snp analysis or window
            Z = zeros(size(fittednullmodel.mm, 1), setlength) # column counts of genomat
            totalsnps = SnpArrays.makestream(countlines, bimfile)
            SnpArrays.makestream(pvalfile, "w") do io
                if test == :score
                    println(io, "startchr,startpos,startsnpid,endchr,endpos,endsnpid,pval")
                    ts = OrdinalMultinomialScoreTest(fittednullmodel.model, Z)
                else # lrt
                    println(io, "startchr,startpos,startsnpid,endchr,",
                    "endpos,endsnpid,l2normeffect,pval")
                    nulldev = deviance(fittednullmodel.model)
                    Xaug = [fittednullmodel.model.X Z]
                    γ̂ = Vector{Float64}(undef, setlength) # effect size for columns being tested
                end
                SnpArrays.makestream(bimfile) do bimio
                    q = setlength
                    for j in 1:q:totalsnps
                        endj = j + q - 1
                        rowj = readline(bimio)
                        if endj >= totalsnps
                            endj = totalsnps
                            q = totalsnps - j + 1
                            #length of Z will be different
                            if test == :score
                                Z = zeros(size(fittednullmodel.mm, 1), q)
                                ts = OrdinalMultinomialScoreTest(fittednullmodel.model, Z)
                            else
                                Xaug = [fittednullmodel.model.X zeros(size(
                                    fittednullmodel.mm, 1), q)]
                            end
                        end
                        for i in 1:(q - 2) #
                            readline(bimio)
                        end
                        endj == totalsnps ? rowj_s = rowj : rowj_s = readline(bimio)
                        snpj = split(rowj)
                        snpj_s = split(rowj_s)
                        if test == :score
                            if all(@view(mafs[j:endj]) .== 0) # all mono-allelic, unlikely but just in case
                                pval = 1.0
                            else
                                copyto!(ts.Z, @view(genomat[bedrowinds, j:endj]), impute=true, model=snpmodel)
                                pval = polrtest(ts)
                            end
                            println(io, "$(snpj[1]),$(snpj[4]),$(snpj[2]),$(snpj_s[1]),",
                                "$(snpj_s[4]),$(snpj_s[2]),$pval")
                        elseif test == :lrt
                            if all(@view(mafs[j:endj]) .== 0) # all mono-allelic, unlikely but just in case
                                fill!(γ̂, 0)
                                pval = 1.0
                            else
                                copyto!(@view(Xaug[:, (fittednullmodel.model.p+1):end]),
                                        @view(genomat[bedrowinds, j:endj]),
                                        impute=true, model=snpmodel)
                                altmodel = polr(Xaug, fittednullmodel.model.Y,
                                    fittednullmodel.model.link, solver,
                                    wts = fittednullmodel.model.wts)
                                copyto!(γ̂, @view(altmodel.β[(fittednullmodel.model.p+1):end]))#, fittednullmodel.model.p + 1, setlength)
                                l2normeffect = norm(γ̂)
                                pval = ccdf(Chisq(q), nulldev - deviance(altmodel))
                            end
                            println(io, "$(snpj[1]),$(snpj[4]),$(snpj[2]),",
                                "$(snpj_s[1]),$(snpj_s[4]),$(snpj_s[2]),$l2normeffect,$pval")
                        end
                    end
                end
            end
        elseif setlength == 0 #snpset is defined by snpset file
            SnpArrays.makestream(pvalfile, "w") do io
                test == :score ? println(io, "snpsetid,nsnps,pval") : println(io,
                    "snpsetid,nsnps,l2normeffect,pval")
                for j in eachindex(snpset_ids)
                    snpset_id = snpset_ids[j]
                    snpinds = findall(snpsetFile[!, :snpset_id] .== snpset_id)
                    q = length(snpinds)
                    Z = zeros(size(fittednullmodel.mm, 1), q)
                    γ̂ = Vector{Float64}(undef, q)
                    Xaug = [fittednullmodel.model.X Z]
                    if all(@view(mafs[snpinds]) .== 0) # all mono-allelic, unlikely but just in case
                        l2normeffect = 0.0
                        pval = 1.0
                        test == :score ? println(io, "$(snpset_id),$q,$pval") :
                        println(io, "$(snpset_id),$q,$l2normeffect,$pval")
                    elseif test == :score
                        ts = OrdinalMultinomialScoreTest(fittednullmodel.model, Z)
                        copyto!(ts.Z, @view(genomat[bedrowinds, snpinds]), impute=true,
                            model=snpmodel)
                        pval = polrtest(ts)
                        println(io, "$(snpset_id),$q,$pval")
                    elseif test == :lrt
                        nulldev = deviance(fittednullmodel.model)
                        copyto!(@view(Xaug[:, fittednullmodel.model.p+1:end]),
                                @view(genomat[bedrowinds, snpinds]),
                                impute=true, model=snpmodel)
                        altmodel = polr(Xaug, fittednullmodel.model.Y,
                            fittednullmodel.model.link, solver,
                            wts = fittednullmodel.model.wts)
                        copyto!(γ̂, 1, altmodel.β, fittednullmodel.model.p + 1, q)
                        l2normeffect = norm(γ̂)
                        pval = ccdf(Chisq(q), nulldev - deviance(altmodel))
                        println(io, "$(snpset_id),$q,$l2normeffect,$pval")
                    end
                end
            end
        else #setlength == -1 (testing just one set with specified snps in snpset)
            snpset = eltype(snpset) == Bool ? findall(snpset) : snpset
            SnpArrays.makestream(pvalfile, "w") do io
                if all(@view(mafs[snpset]) .== 0) # all mono-allelic, unlikely but just in case
                    l2normeffect = 0.0
                    pval = 1.0
                    test == :score ? println(io, "The joint pvalue of snps indexed",
                    " at $(snpset) is $pval") : println(io, "The l2norm of the effect size vector",
                    " is $l2normeffect and joint pvalue of snps indexed",
                    " at $(snpset) is $pval")
                else
                    q = length(snpset)
                    γ̂ = Vector{Float64}(undef, q)
                    Z = zeros(size(fittednullmodel.mm, 1), q)
                    if test == :score
                        ts = OrdinalMultinomialScoreTest(fittednullmodel.model, Z)
                        copyto!(ts.Z, @view(genomat[bedrowinds, snpset]), impute=true, model=snpmodel)
                        pval = polrtest(ts)
                        println(io, "The joint pvalue of snps indexed",
                         " at $(snpset) is $pval")
                    elseif test == :lrt
                        nulldev = deviance(fittednullmodel.model)
                        Xaug = [fittednullmodel.model.X Z]
                        copyto!(@view(Xaug[:, fittednullmodel.model.p+1:end]),
                                @view(genomat[bedrowinds, snpset]),
                                impute=true, model=snpmodel)
                        altmodel = polr(Xaug, fittednullmodel.model.Y,
                            fittednullmodel.model.link, solver,
                            wts = fittednullmodel.model.wts)
                        copyto!(γ̂, 1, altmodel.β, fittednullmodel.model.p + 1, q)
                        l2normeffect = norm(γ̂)
                        pval = ccdf(Chisq(q), nulldev - deviance(altmodel))
                        println(io, "The l2norm of the effect size vector",
                        " is $l2normeffect and joint pvalue of snps indexed",
                        " at $(snpset) is $pval")
                    end
                end
            end
        end
    else #analysistype == "gxe"
        isnothing(e) &&
            @error("GxE analysis indicated but not environmental variable keyword argument: `e` set.")
        Xaug = [fittednullmodel.model.X zeros(size(fittednullmodel.mm, 1))]
        Xaug2 = [fittednullmodel.model.X zeros(size(fittednullmodel.mm, 1), 2)]
        envvar = modelmatrix(FormulaTerm(fittednullmodel.mf.f.lhs, Term(Symbol(e))),
                 DataFrame(fittednullmodel.mf.data))
        testvec = Matrix{Float64}(undef, size(envvar))
        snpeffectnull = 0.0
        SnpArrays.makestream(pvalfile, "w") do io
            if test == :score
                println(io, "chr,pos,snpid,maf,hwepval,snpeffectnull,pval")
            else
                γ̂ = 0.0 # effect size for columns being tested
                println(io, "chr,pos,snpid,maf,hwepval,snpeffectnull,snpeffectfull,GxEeffect,pval")
            end
            SnpArrays.makestream(bimfile) do bimio
                for j in eachindex(snpmask)
                    row = readline(bimio)
                    snpmask[j] || continue
                    hwepval = SnpArrays.hwe(cc[1, j], cc[3, j], cc[4, j])
                    maf = mafs[j]
                    snpj = split(row)
                    if maf == 0 # mono-allelic
                        γ̂ = 0.0
                        pval = 1.0
                        snpeffectfull = 0.0
                        snpeffectnull = 0.0
                    elseif test == :score
                        copyto!(@view(Xaug[:, end]), @view(genomat[bedrowinds,
                        j]), impute=true, model=snpmodel)
                        copyto!(testvec, @view(Xaug[:, end]) .* envvar)
                        nm = polr(Xaug, fittednullmodel.model.Y,
                            fittednullmodel.model.link, solver, wts = fittednullmodel.model.wts)
                        if Inf in nm.vcov #singular design matrix
                            snpeffectnull = 0.0
                            pval = 1.0
                        else
                            snpeffectnull = nm.β[end]
                            ts = OrdinalMultinomialScoreTest(nm, testvec)
                            pval = polrtest(ts)
                        end
                    elseif test == :lrt
                        copyto!(@view(Xaug[:, end]), @view(genomat[bedrowinds,
                            j]), impute=true, model=snpmodel)
                        copyto!(@view(Xaug2[:, end - 1]), @view(Xaug[:, end]))
                        copyto!(@view(Xaug2[:, end]), @view(Xaug[:, end]) .*
                            envvar)
                        nm = polr(Xaug, fittednullmodel.model.Y,
                        fittednullmodel.model.link, solver, wts = fittednullmodel.model.wts)
                        if Inf in nm.vcov #singular design matrix
                            snpeffectnull = 0.0
                            snpeffectfull = 0.0
                            γ̂ = 0.0
                            pval = 1.0
                        else
                            snpeffectnull = nm.β[end]
                            nulldev = deviance(nm)
                            altmodel = polr(Xaug2, fittednullmodel.model.Y,
                                fittednullmodel.model.link, solver,
                                wts = fittednullmodel.model.wts)
                            γ̂ = altmodel.β[end]
                            snpeffectfull = altmodel.β[end-1]
                            pval = ccdf(Chisq(1), nulldev - deviance(altmodel))
                        end
                    end
                    if test == :score
                        println(io, "$(snpj[1]),$(snpj[4]),$(snpj[2]),$maf,",
                        "$hwepval,$snpeffectnull,$pval")
                    else
                        println(io, "$(snpj[1]),$(snpj[4]),$(snpj[2]),$maf,",
                        "$hwepval,$snpeffectnull,$snpeffectfull,$γ̂,$pval")
                    end
                end
            end
        end
    end
    return fittednullmodel
end

# For VCF Analysis
function ordinalgwas(
    fittednullmodel::StatsModels.TableRegressionModel,
    vcffile::Union{AbstractString, IOStream}, # full path and vcf file name
    nsamples::Integer,          # number of samples in bed file
    vcftype::Symbol;           # :GT = genotype, :DS = dosage
    analysistype::AbstractString = "singlesnp",
    testformula::FormulaTerm = fittednullmodel.mf.f.lhs ~ Term(:snp),
    test::Symbol = :score,
    pvalfile::Union{AbstractString, IOStream} = "ordinalgwas.pval.txt",
    snpmodel::Union{Val{1}, Val{2}, Val{3}} = ADDITIVE_MODEL,
    snpinds::Union{Nothing, AbstractVector{<:Integer}} = nothing,
    vcfrowinds::AbstractVector{<:Integer} = 1:nsamples, # row indices for VCF array
    solver::MOI.AbstractOptimizer = NLopt.Optimizer(),
    solver_config = Dict("algorithm" => :LD_SLSQP, "max_iter" => 4000),
    verbose::Bool = false,
    snpset::Union{Nothing, Integer, AbstractString, #for snpset analysis
        AbstractVector{<:Integer}} = nothing,
    e::Union{Nothing, AbstractString, Symbol} = nothing # for GxE analysis
    )
    config_solver(solver, solver_config)
    # get number of SNPs in file
    nsnps = nrecords(vcffile)
    # these are currently only based on genotype data -- not dosage. Comment out.
    # nsnps, _, _, _, _, mafs, _, hwes = gtstats(vcffile)

    # for VCFTools, snpmodel is coded differently
    snpmodel = modelingdict[snpmodel]

    # create SNP mask vector
    if snpinds === nothing
        snpmask = trues(nsnps)
    elseif eltype(snpinds) == Bool
        snpmask = snpinds
    else
        snpmask = falses(nsnps)
        snpmask[snpinds] .= true
    end

    analysistype = lowercase(analysistype)
    analysistype in ["singlesnp", "snpset", "gxe"] || error("Analysis type $analysis invalid option.
    Available options are 'singlesnp', 'snpset' and 'gxe'.")
    # open VCF File
    reader = VCF.Reader(openvcf(vcffile))

    # determine analysis type
    if analysistype == "singlesnp"
        # extra columns in design matrix to be tested
        testdf = DataFrame(fittednullmodel.mf.data) # TODO: not type stable here
        testdf[!, :snp] = zeros(size(fittednullmodel.mm, 1))
        Z = similar(modelmatrix(testformula, testdf))

        # create holders for chromome, position, id, dosage/snps
        rec_chr = Array{Any, 1}(undef, 1)
        rec_pos = Array{Any, 1}(undef, 1)
        rec_ids = Array{Any, 1}(undef, 1)
        gholder = zeros(Union{Missing, Float64}, nsamples)

        # carry out score or LRT test SNP by SNP
        snponly = testformula.rhs == Term(:snp)
        SnpArrays.makestream(pvalfile, "w") do io
            if test == :score
                println(io, "chr,pos,snpid,pval")
                ts = OrdinalMultinomialScoreTest(fittednullmodel.model, Z)
            else
                nulldev = deviance(fittednullmodel.model)
                Xaug = [fittednullmodel.model.X Z]
                q = size(Z, 2)
                γ̂ = Vector{Float64}(undef, q) # effect size for columns being tested
                if snponly
                    println(io, "chr,pos,snpid,effect,stder,pval")
                else
                    print(io, "chr,pos,snpid,")
                    for j in 1:q
                        print(io, "effect$j,")
                    end
                    println(io, "pval")
                end
            end
            for j in eachindex(snpmask)
                if !snpmask[j] #skip snp, must read marker still.
                    copy_gt!(gholder, reader; model = snpmodel, impute = true,
                    record_chr = rec_chr, record_pos = rec_pos, record_ids = rec_ids)
                    continue
                end
                if vcftype == :GT #genotype
                    copy_gt!(gholder, reader; model = snpmodel, impute = true,
                    record_chr = rec_chr, record_pos = rec_pos, record_ids = rec_ids)
                else #dosage
                    copy_ds!(gholder, reader; model = snpmodel, impute = true,
                    record_chr = rec_chr, record_pos = rec_pos, record_ids = rec_ids)
                end
                if test == :score
                    if snponly
                        copyto!(ts.Z, @view(gholder[vcfrowinds]))
                    else # snp + other terms
                        copyto!(testdf[!, :snp], @view(gholder[vcfrowinds]))
                        ts.Z[:] = modelmatrix(testformula, testdf)
                    end
                    if var(@view(gholder[vcfrowinds])) == 0.0 #mafs[j] == 0.0
                        pval = 1.0
                    else
                        pval = polrtest(ts)
                    end
                    println(io, "$(rec_chr[1]),$(rec_pos[1]),$(rec_ids[1][1]),",
                    # "$(mafs[j]),$(hwes[j]),"
                    "$pval")
                elseif test == :lrt
                    if snponly
                        copyto!(@view(Xaug[:, fittednullmodel.model.p+1]),
                            @view(gholder[vcfrowinds]))
                    else # snp + other terms
                        copyto!(testdf[!, :snp], @view(gholder[vcfrowinds]))
                        Xaug[:, fittednullmodel.model.p+1:end] = modelmatrix(testformula, testdf)
                    end
                    if var(@view(gholder[vcfrowinds])) == 0.0 # mafs[j] == 0.0
                        pval = 1.0
                        stderr = NaN
                    else
                        altmodel = polr(Xaug, fittednullmodel.model.Y,
                            fittednullmodel.model.link, solver,
                            wts = fittednullmodel.model.wts)
                        copyto!(γ̂, 1, altmodel.β, fittednullmodel.model.p + 1, q)
                        pval = ccdf(Chisq(q), nulldev - deviance(altmodel))
                        stderr = stderror(altmodel)[fittednullmodel.model.npar + 1]
                    end
                    if snponly
                        println(io, "$(rec_chr[1]),$(rec_pos[1]),$(rec_ids[1][1]),",
                        # "$(mafs[j]),$(hwes[j]),",
                        "$(γ̂[1]),$stderr,$pval")
                    else
                        print(io, "$(rec_chr[1]),$(rec_pos[1]),$(rec_ids[1][1]),",
                        # "$(mafs[j]),$(hwes[j]),"
                        )
                        for j in 1:q
                            print(io, "$(γ̂[j]),")
                        end
                        println(io, pval)
                    end
                end
            end
        end
    elseif analysistype == "snpset"
        # max size of a snpset length
        maxsnpset = 1

        #determine snpset
        if isa(snpset, Nothing)
            setlength = 1
            maxsnpset = 1
        elseif isa(snpset, AbstractString)
            isfile(snpset) || throw(ArgumentError("snpset file not found,
            to specify a window replace snpset string with a window size"))
            #first column SNPset ID, second column SNP ID
            snpsetFile = CSV.read(snpset, DataFrame, header = [:snpset_id, :snp_id], delim = " ")
            maxsnpset = combine(groupby(snpsetFile, :snpset_id), :snp_id => length => :snpset_length) |>
                x -> maximum(x.snpset_length)
            snpset_ids = unique(snpsetFile[!, :snpset_id])
            nSets = length(snpset_ids)
            setlength = 0
        elseif isa(snpset, Integer)
            setlength = snpset
            maxsnpset = snpset
        else #abstract vector (boolean of true at indicies or range or indicies)
            setlength = -1
            maxsnpset = count(snpset .!= 0)
        end

        # create holders for chromome, position, id, dosage/gt
        rec_chr = Array{Any, 1}(undef, maxsnpset)
        rec_pos = Array{Any, 1}(undef, maxsnpset)
        rec_ids = Array{Any, 1}(undef, maxsnpset)
        gholder = zeros(Union{Missing, Float64}, nsamples, maxsnpset)

        if setlength > 0 #single snp analysis or window
            Z = zeros(size(fittednullmodel.mm, 1), setlength) #
            q = setlength
            SnpArrays.makestream(pvalfile, "w") do io
                if test == :score
                    println(io, "startchr,startpos,startsnpid,endchr,",
                    "endpos,endsnpid,pval")
                    ts = OrdinalMultinomialScoreTest(fittednullmodel.model, Z)
                elseif test == :lrt
                    println(io, "startchr,startpos,startsnpid,endchr,",
                    "endpos,endsnpid,l2normeffect,pval")
                    nulldev = deviance(fittednullmodel.model)
                    Xaug = [fittednullmodel.model.X Z]
                    γ̂ = Vector{Float64}(undef, setlength) # effect size for columns being tested
                end
                for j in 1:q:nsnps
                    endj = j + q - 1
                    if endj >= nsnps
                        endj = nsnps
                        q = nsnps - j + 1
                        #length of Z will be different
                        gholder = zeros(Union{Missing, Float64}, nsamples, q)
                        rec_chr = Array{Any, 1}(undef, q)
                        rec_pos = Array{Any, 1}(undef, q)
                        rec_ids = Array{Any, 1}(undef, q)
                        if test == :score
                            Z = zeros(size(fittednullmodel.mm, 1), q)
                            ts = OrdinalMultinomialScoreTest(fittednullmodel.model, Z)
                        elseif test == :lrt
                            Xaug = [fittednullmodel.model.X zeros(size(
                            fittednullmodel.mm, 1), q)]
                        end
                    end
                    if vcftype == :GT #genotype
                        copy_gt!(gholder, reader; model = snpmodel, impute = true,
                        record_chr = rec_chr, record_pos = rec_pos, record_ids = rec_ids)
                    else #dosage
                        copy_ds!(gholder, reader; model = snpmodel, impute = true,
                        record_chr = rec_chr, record_pos = rec_pos, record_ids = rec_ids)
                    end
                    if test == :score
                        # if all(@view(mafs[j:endj]) .== 0) # all mono-allelic, unlikely but just in case
                        if all(var(@view(gholder[vcfrowinds, :]), dims = [1]) == 0.0)
                            pval = 1.0
                        else
                            copyto!(ts.Z, @view(gholder[vcfrowinds, :]))
                            pval = polrtest(ts)
                        end
                        println(io, "$(rec_chr[1]),$(rec_pos[1]),$(rec_ids[1][1]),",
                        "$(rec_chr[end]),$(rec_pos[end]),$(rec_ids[end][end]),$pval")
                    elseif test == :lrt
                        # if all(@view(mafs[j:endj]) .== 0) # all mono-allelic, unlikely but just in case
                        if all(var(@view(gholder[vcfrowinds, :]), dims = [1]) == 0.0)
                            l2normeffect = 0.0
                            pval = 1.0
                        else
                            copyto!(@view(Xaug[:, (fittednullmodel.model.p+1):end]),
                            @view(gholder[vcfrowinds]))
                            altmodel = polr(Xaug, fittednullmodel.model.Y,
                            fittednullmodel.model.link, solver,
                            wts = fittednullmodel.model.wts)
                            copyto!(γ̂, @view(altmodel.β[(fittednullmodel.model.p+1):end]))#, fittednullmodel.model.p + 1, setlength)
                            l2normeffect = norm(γ̂)
                            pval = ccdf(Chisq(q), nulldev - deviance(altmodel))
                        end
                        println(io, "$(rec_chr[1]),$(rec_pos[1]),$(rec_ids[1][1]),",
                        "$(rec_chr[end]),$(rec_pos[end]),$(rec_ids[end][end]),",
                        "$l2normeffect,$pval")
                    end
                end
            end
        elseif setlength == 0 #snpset is defined by snpset file
            @warn("This method requires reading in the entire VCF File.
             This can take a lot of memory for large files, as they must be brought into memory.")
            if vcftype == :GT #genotype
                genomat = convert_gt(Float64, vcffile;
                model = snpmodel, impute = true,
                center = false, scale = false)
            else #dosage
                genomat = convert_ds(Float64, vcffile; model = snpmodel,
                key="DS", impute = true, center = false, scale = false)
            end
            SnpArrays.makestream(pvalfile, "w") do io
                test == :score ? println(io, "snpsetid,nsnps,pval") : println(io,
                    "snpsetid,nsnps,l2normeffect,pval")
                for j in eachindex(snpset_ids)
                    snpset_id = snpset_ids[j]
                    snpinds = findall(snpsetFile[!, :snpset_id] .== snpset_id)
                    q = length(snpinds)
                    Z = zeros(size(fittednullmodel.mm, 1), q)
                    if test == :score
                        # if all(@view(mafs[snpinds]) .== 0)
                        if all(var(@view(genomat[vcfrowinds, snpinds]), dims = [1]) == 0.0)
                            pval = 1.0
                        else
                            ts = OrdinalMultinomialScoreTest(fittednullmodel.model, Z)
                            copyto!(ts.Z, @view(genomat[vcfrowinds, snpinds]))
                            pval = polrtest(ts)
                        end
                        println(io, "$(snpset_id),$q,$pval")
                    elseif test == :lrt
                        # if all(@view(mafs[snpinds]) .== 0)
                        if all(var(@view(genomat[vcfrowinds, snpinds]), dims = [1]) == 0.0)
                            l2normeffect = 0.0
                            pval = 1.0
                        else
                            γ̂ = Vector{Float64}(undef, q)
                            Xaug = [fittednullmodel.model.X Z]
                            nulldev = deviance(fittednullmodel.model)
                            copyto!(@view(Xaug[:, fittednullmodel.model.p+1:end]),
                                    @view(genomat[vcfrowinds, snpinds]))
                            altmodel = polr(Xaug, fittednullmodel.model.Y,
                                fittednullmodel.model.link, solver,
                                wts = fittednullmodel.model.wts)
                            copyto!(γ̂, 1, altmodel.β, fittednullmodel.model.p + 1, q)
                            l2normeffect = norm(γ̂)
                            pval = ccdf(Chisq(q), nulldev - deviance(altmodel))
                        end
                        println(io, "$(snpset_id),$q,$l2normeffect,$pval")
                    end
                end
            end
        else #setlength == -1 (testing just one set with specified snps in snpset)
            @warn("This method requires reading in the entire VCF File.
            This can take a lot of memory for large files, as they must be brought into memory.")
            snpset = eltype(snpset) == Bool ? findall(snpset) : snpset
            if vcftype == :GT #genotype
                genomat = convert_gt(Float64, vcffile;
                model = snpmodel, impute = true,
                center = false, scale = false)
            else #dosage
                genomat = convert_ds(Float64, vcffile; model = snpmodel,
                key="DS", impute=true, center=false, scale=false)
            end
            SnpArrays.makestream(pvalfile, "w") do io
                q = length(snpset)
                γ̂ = Vector{Float64}(undef, q)
                Z = zeros(size(fittednullmodel.mm, 1), q)
                if all(@view(genomat[vcfrowinds, snpset]) .== 0)
                    println(io, "All SNPs had 0 variation. No results to report.")
                elseif test == :score
                    ts = OrdinalMultinomialScoreTest(fittednullmodel.model, Z)
                    copyto!(ts.Z, @view(genomat[vcfrowinds, snpset]))
                    pval = polrtest(ts)
                    println(io, "The joint pvalue of snps indexed",
                        " at $(snpset) is $pval")
                elseif test == :lrt
                    nulldev = deviance(fittednullmodel.model)
                    Xaug = [fittednullmodel.model.X Z]
                    copyto!(@view(Xaug[:, fittednullmodel.model.p+1:end]),
                            @view(genomat[vcfrowinds, snpset]))
                    altmodel = polr(Xaug, fittednullmodel.model.Y,
                        fittednullmodel.model.link, solver,
                        wts = fittednullmodel.model.wts)
                    copyto!(γ̂, 1, altmodel.β, fittednullmodel.model.p + 1, q)
                    l2normeffect = norm(γ̂)
                    pval = ccdf(Chisq(q), nulldev - deviance(altmodel))
                    println(io, "The l2norm of the effect size vector",
                    " is $l2normeffect and joint pvalue of snps indexed",
                    " at $(snpset) is $pval")
                end
            end
        end
    else #analysistype == "gxe"
        isnothing(e) &&
            @error("GxE analysis indicated but not environmental variable keyword argument: `e` set.")
        # create holders for chromome, position, id
        rec_chr = Array{Any, 1}(undef, 1)
        rec_pos = Array{Any, 1}(undef, 1)
        rec_ids = Array{Any, 1}(undef, 1)
        gholder = zeros(Union{Missing, Float64}, nsamples)
        Xaug = [fittednullmodel.model.X zeros(size(fittednullmodel.mm, 1))]
        Xaug2 = [fittednullmodel.model.X zeros(size(fittednullmodel.mm, 1), 2)] #or get Xaug to point to part of it

        # create array for environmental variable and testing
        envvar = modelmatrix(FormulaTerm(fittednullmodel.mf.f.lhs, Term(Symbol(e))),
                 DataFrame(fittednullmodel.mf.data))
        testvec = Matrix{Float64}(undef, size(envvar))
        snpeffectnull = 0.0
        SnpArrays.makestream(pvalfile, "w") do io
            if test == :score
                println(io, "chr,pos,snpid,snpeffectnull,pval")
            else
                println(io, "chr,pos,snpid,snpeffectnull,",
                "snpeffectfull,GxEeffect,pval")
            end
            for j in eachindex(snpmask)
                if !snpmask[j] #skip snp, must read marker still.
                    copy_gt!(gholder, reader; model = snpmodel, impute = true,
                    record_chr = rec_chr, record_pos = rec_pos, record_ids = rec_ids)
                    continue
                end
                if vcftype == :GT #genotype
                    copy_gt!(gholder, reader; model = snpmodel, impute = true,
                    record_chr = rec_chr, record_pos = rec_pos, record_ids = rec_ids)
                else #dosage
                    copy_ds!(gholder, reader; model = snpmodel, impute = true,
                    record_chr = rec_chr, record_pos = rec_pos, record_ids = rec_ids)
                end
                copyto!(@view(Xaug[:, end]), @view(gholder[vcfrowinds]))
                zeromaf = var(@view(gholder[vcfrowinds])) == 0
                if test == :score
                    copyto!(testvec, @view(Xaug[:, end]) .* envvar)
                    nm = polr(Xaug, fittednullmodel.model.Y,
                        fittednullmodel.model.link, solver, wts = fittednullmodel.model.wts)
                    if (Inf in nm.vcov) || zeromaf #(mafs[j] == 0.0) #singular design matrix
                        snpeffectnull = 0.0
                        pval = 1.0
                    else
                        snpeffectnull = nm.β[end]
                        ts = OrdinalMultinomialScoreTest(nm, testvec)
                        pval = polrtest(ts)
                    end
                    println(io, "$(rec_chr[1]),$(rec_pos[1]),$(rec_ids[1][1]),",
                    # "$(mafs[j]),$(hwes[j]),",
                    "$snpeffectnull,$pval")
                elseif test == :lrt
                    γ̂ = 0.0 # effect size for columns being tested
                    copyto!(@view(Xaug2[:, end - 1]), @view(Xaug[:, end]))
                    copyto!(@view(Xaug2[:, end]), @view(Xaug[:, end]) .*
                        envvar)
                    nm = polr(Xaug, fittednullmodel.model.Y,
                    fittednullmodel.model.link, solver, wts = fittednullmodel.model.wts)
                    if (Inf in nm.vcov) || zeromaf #(mafs[j] == 0.0) #singular design matrix
                        snpeffectnull = 0.0
                        snpeffectfull = 0.0
                        γ̂ = 0.0
                        pval = 1.0
                    else
                        snpeffectnull = nm.β[end]
                        nulldev = deviance(nm)
                        altmodel = polr(Xaug2, fittednullmodel.model.Y,
                            fittednullmodel.model.link, solver,
                            wts = fittednullmodel.model.wts)
                        γ̂ = altmodel.β[end]
                        snpeffectfull = altmodel.β[end-1]
                        pval = ccdf(Chisq(1), nulldev - deviance(altmodel))
                    end
                    println(io, "$(rec_chr[1]),$(rec_pos[1]),$(rec_ids[1][1]),",
                        #"$(mafs[j]),$(hwes[j]),",
                        "$snpeffectnull,$snpeffectfull,$γ̂,$pval")
                end
            end
        end
    end
    close(reader)
    return fittednullmodel
end

# For BGEN Analysis
function ordinalgwas(
    fittednullmodel::StatsModels.TableRegressionModel,
    bgenfile::Union{AbstractString, IOStream}, # full path and bgen file name
    nsamples::Integer;          # number of samples in bed file
    analysistype::AbstractString = "singlesnp",
    testformula::FormulaTerm = fittednullmodel.mf.f.lhs ~ Term(:snp),
    test::Symbol = :score,
    pvalfile::Union{AbstractString, IOStream} = "ordinalgwas.pval.txt",
    snpmodel::Union{Val{1}, Val{2}, Val{3}} = ADDITIVE_MODEL,
    snpinds::Union{Nothing, AbstractVector{<:Integer}} = nothing,
    bgenrowinds::AbstractVector{<:Integer} = 1:nsamples, # row indices for VCF array
    solver::MOI.AbstractOptimizer = NLopt.Optimizer(),
    solver_config = Dict("algorithm" => :LD_SLSQP, "max_iter" => 4000),
    verbose::Bool = false,
    snpset::Union{Nothing, Integer, AbstractString, #for snpset analysis
        AbstractVector{<:Integer}} = nothing,
    e::Union{Nothing, AbstractString, Symbol} = nothing # for GxE analysis
    )
    config_solver(solver, solver_config)
    # open BGEN file and get number of SNPs in file
    bgendata = Bgen(bgenfile)
    nsnps = n_variants(bgendata)
    bgen_iterator = iterator(bgendata)

    dosageholder = Vector{Float32}(undef, n_samples(bgendata))
    decompressed_length, _ = BGEN.check_decompressed_length(
        bgendata.io, first(bgen_iterator), bgendata.header)
    decompressed = Vector{UInt8}(undef, decompressed_length)
    bgenrowmask_UInt16 = zeros(UInt16, n_samples(bgendata))
    bgenrowmask_UInt16[bgenrowinds] .= 1

    # create SNP mask vector
    if snpinds === nothing
        snpmask = trues(nsnps)
    elseif eltype(snpinds) == Bool
        snpmask = snpinds
    else
        snpmask = falses(nsnps)
        snpmask[snpinds] .= true
    end

    analysistype = lowercase(analysistype)
    analysistype in ["singlesnp", "snpset", "gxe"] || error("Analysis type $analysis invalid option.
    Available options are 'singlesnp', 'snpset' and 'gxe'.")

    # determine analysis type
    if analysistype == "singlesnp"
        # extra columns in design matrix to be tested
        testdf = DataFrame(fittednullmodel.mf.data) # TODO: not type stable here
        testdf[!, :snp] = zeros(size(fittednullmodel.mm, 1))
        Z = similar(modelmatrix(testformula, testdf))

        # create holder for dosage/snps
        snpholder = zeros(Union{Missing, Float64}, size(fittednullmodel.mm, 1))

        # carry out score or LRT test SNP by SNP
        snponly = testformula.rhs == Term(:snp)
        SnpArrays.makestream(pvalfile, "w") do io
            if test == :score
                println(io, "chr,pos,snpid,varid,maf,hwepval,infoscore,pval")
                ts = OrdinalMultinomialScoreTest(fittednullmodel.model, Z)
            else
                nulldev = deviance(fittednullmodel.model)
                Xaug = [fittednullmodel.model.X Z]
                q = size(Z, 2)
                γ̂ = Vector{Float64}(undef, q) # effect size for columns being tested
                if snponly
                    println(io, "chr,pos,snpid,varid,maf,hwepval,",
                    "infoscore,effect,stder,pval")
                else
                    print(io, "chr,pos,snpid,varid,maf,hwepval,infoscore,")
                    for j in 1:q
                        print(io, "effect$j,")
                    end
                    println(io, "pval")
                end
            end
            for (j, variant) in enumerate(bgen_iterator)
                if !snpmask[j] #skip snp
                    continue
                end
                minor_allele_dosage!(bgendata, variant;
                    T = Float64, mean_impute = true, data = dosageholder,
                    decompressed = decompressed)
                @views copyto!(snpholder, dosageholder[bgenrowinds])
                hwepval = 9.0
                maf = 9.0
                infoscore = 9.0
                try
                    hwepval = BGEN.hwe(bgendata, variant; rmask = bgenrowmask_UInt16)
                catch nothing
                end
                try
                    maf = BGEN.maf(bgendata, variant; rmask = bgenrowmask_UInt16)
                catch nothing
                end
                try
                    infoscore = BGEN.info_score(bgendata, variant; rmask = bgenrowmask_UInt16)
                catch nothing
                end
                if test == :score
                    if snponly
                        copyto!(ts.Z, snpholder)
                    else # snp + other terms
                        copyto!(testdf[!, :snp], snpholder)
                        ts.Z[:] = modelmatrix(testformula, testdf)
                    end
                    if maf == 0.0
                        pval = 1.0
                    else
                        pval = polrtest(ts)
                    end
                    println(io, "$(variant.chrom),$(variant.pos),$(variant.rsid),",
                    "$(variant.varid),$(maf),$(hwepval),$(infoscore),",
                    "$pval")
                elseif test == :lrt
                    if snponly
                        copyto!(@view(Xaug[:, fittednullmodel.model.p+1]),
                        snpholder)
                    else # snp + other terms
                        copyto!(testdf[!, :snp], snpholder)
                        Xaug[:, fittednullmodel.model.p+1:end] = modelmatrix(testformula, testdf)
                    end
                    if maf == 0.0
                        pval = 1.0
                        stderr = NaN
                    else
                        altmodel = polr(Xaug, fittednullmodel.model.Y,
                        fittednullmodel.model.link, solver,
                        wts = fittednullmodel.model.wts)
                        copyto!(γ̂, 1, altmodel.β, fittednullmodel.model.p + 1, q)
                        pval = ccdf(Chisq(q), nulldev - deviance(altmodel))
                        stderr = stderror(altmodel)[fittednullmodel.model.npar + 1]
                    end
                    if snponly
                        println(io, "$(variant.chrom),$(variant.pos),$(variant.rsid),",
                        "$(variant.varid),$(maf),$(hwepval),$(infoscore),",
                        "$(γ̂[1]),$stderr,$pval")
                    else
                        print(io, "$(variant.chrom),$(variant.pos),$(variant.rsid),",
                        "$(variant.varid),$(maf),$(hwepval),$(infoscore),")
                        for j in 1:q
                            print(io, "$(γ̂[j]),")
                        end
                        println(io, pval)
                    end
                end
            end
        end
    elseif analysistype == "snpset"
        # max size of a snpset length
        maxsnpset = 1

        #determine snpset
        if isa(snpset, Nothing)
            setlength = 1
            maxsnpset = 1
        elseif isa(snpset, AbstractString)
            isfile(snpset) || throw(ArgumentError("snpset file not found,
            to specify a window replace snpset string with a window size"))
            #first column SNPset ID, second column SNP ID
            snpsetFile = CSV.read(snpset, DataFrame, header = [:snpset_id, :snp_id], delim = " ")
            maxsnpset = combine(groupby(snpsetFile, :snpset_id), :snp_id => length => :snpset_length) |>
                x -> maximum(x.snpset_length)
            snpset_ids = unique(snpsetFile[!, :snpset_id])
            nSets = length(snpset_ids)
            setlength = 0
        elseif isa(snpset, Integer)
            setlength = snpset
            maxsnpset = snpset
        else #abstract vector (boolean of true at indicies or range or indicies)
            setlength = -1
            maxsnpset = count(snpset .!= 0)
        end

        # create holders for chromome, position, id, dosage/gt
        snpholder = zeros(Union{Missing, Float64},
            size(fittednullmodel.mm, 1), maxsnpset)
        dosageholder = Vector{Float32}(undef, n_samples(bgendata))
        decompressed = Vector{UInt8}(undef, 3 * n_samples(bgendata) + 10)

        if setlength > 0 #single snp analysis or window
            Z = zeros(size(fittednullmodel.mm, 1), setlength) #
            q = setlength
            SnpArrays.makestream(pvalfile, "w") do io
                if test == :score
                    println(io, "startchr,startpos,startsnpid,startvarid,",
                        "endchr,endpos,endsnpid,endvarid,pval")
                    ts = OrdinalMultinomialScoreTest(fittednullmodel.model, Z)
                elseif test == :lrt
                    println(io, "startchr,startpos,startsnpid,startvarid,",
                    "endchr,endpos,endsnpid,endvarid,l2normeffect,pval")
                    nulldev = deviance(fittednullmodel.model)
                    Xaug = [fittednullmodel.model.X Z]
                    γ̂ = Vector{Float64}(undef, setlength) # effect size for columns being tested
                end
                chrstart, posstart, rsidstart, varidstart = "", "", "", ""
                chrend, posend, rsidend, varidend = "", "", "", ""
                for j in 1:q:nsnps
                    endj = j + q - 1
                    if endj >= nsnps
                        q = nsnps - j + 1
                        #length of Z will be different
                        snpholder = zeros(size(fittednullmodel.mm, 1), q)
                        if test == :score
                            Z = zeros(size(fittednullmodel.mm, 1), q)
                            ts = OrdinalMultinomialScoreTest(fittednullmodel.model, Z)
                        elseif test == :lrt
                            Xaug = [fittednullmodel.model.X zeros(size(
                            fittednullmodel.mm, 1), q)]
                        end
                    end
                    for i in 1:q
                        variant = variant_by_index(bgendata, j + i - 1)
                        minor_allele_dosage!(bgendata, variant;
                        T = Float64, mean_impute = true, data = dosageholder,
                        decompressed = decompressed)
                        @views copyto!(snpholder[:, i], dosageholder[bgenrowinds])
                        if i == 1
                            chrstart = variant.chrom
                            posstart = variant.pos
                            rsidstart = variant.rsid
                            varidstart = variant.rsid
                        end
                        if i == q
                            chrend = variant.chrom
                            posend = variant.pos
                            rsidend = variant.rsid
                            varidend = variant.rsid
                        end
                    end
                    if test == :score
                        if all(var(snpholder, dims = [1]) .== 0)
                            pval = 1.0
                        else
                            copyto!(ts.Z, snpholder)
                            pval = polrtest(ts)
                        end
                        println(io, "$chrstart,$posstart,$rsidstart,$varidstart,",
                            "$chrend,$posend,$rsidend,$varidend,$pval")
                    elseif test == :lrt
                        if all(var(snpholder, dims = [1]) .== 0)
                            l2normeffect = 0.0
                            pval = 1.0
                        else
                            copyto!(@view(Xaug[:, (fittednullmodel.model.p+1):end]),
                            snpholder)
                            altmodel = polr(Xaug, fittednullmodel.model.Y,
                            fittednullmodel.model.link, solver,
                            wts = fittednullmodel.model.wts)
                            copyto!(γ̂, @view(altmodel.β[(fittednullmodel.model.p+1):end]))#, fittednullmodel.model.p + 1, setlength)
                            l2normeffect = norm(γ̂)
                            pval = ccdf(Chisq(q), nulldev - deviance(altmodel))
                        end
                        println(io, "$chrstart,$posstart,$rsidstart,$varidstart,",
                        "$chrend,$posend,$rsidend,$varidend,$l2normeffect,$pval")
                    end
                end
            end
        elseif setlength == 0 #snpset is defined by snpset file
            SnpArrays.makestream(pvalfile, "w") do io
                test == :score ? println(io, "snpsetid,nsnps,pval") : println(io,
                    "snpsetid,nsnps,l2normeffect,pval")
                for j in eachindex(snpset_ids)
                    snpset_id = snpset_ids[j]
                    snpinds = findall(snpsetFile[!, :snpset_id] .== snpset_id)
                    q = length(snpinds)
                    Z = zeros(size(fittednullmodel.mm, 1), q)
                    for i in 1:q
                        variant = variant_by_index(bgendata, snpinds[i])
                        minor_allele_dosage!(bgendata, variant;
                            T = Float64, mean_impute = true, data = dosageholder,
                            decompressed = decompressed)
                        @views copyto!(Z[:, i], dosageholder[bgenrowinds])
                    end
                    if test == :score
                        if all(var(Z, dims = [1]) .== 0)
                            pval = 1.0
                        else
                            ts = OrdinalMultinomialScoreTest(fittednullmodel.model, Z)
                            pval = polrtest(ts)
                        end
                        println(io, "$(snpset_id),$q,$pval")
                    elseif test == :lrt
                        if all(var(Z, dims = [1]) .== 0)
                            l2normeffect = 0.0
                            pval = 1.0
                        else
                            γ̂ = Vector{Float64}(undef, q)
                            Xaug = [fittednullmodel.model.X Z]
                            nulldev = deviance(fittednullmodel.model)
                            altmodel = polr(Xaug, fittednullmodel.model.Y,
                                fittednullmodel.model.link, solver,
                                wts = fittednullmodel.model.wts)
                            copyto!(γ̂, 1, altmodel.β, fittednullmodel.model.p + 1, q)
                            l2normeffect = norm(γ̂)
                            pval = ccdf(Chisq(q), nulldev - deviance(altmodel))
                        end
                        println(io, "$(snpset_id),$q,$l2normeffect,$pval")
                    end
                end
            end
        else #setlength == -1 (testing just one set with specified snps in snpset)
            snpset = eltype(snpset) == Bool ? findall(snpset) : snpset
            q = length(snpset)
            γ̂ = Vector{Float64}(undef, q)
            Z = zeros(size(fittednullmodel.mm, 1), q)
            for i in 1:length(snpset)
                variant = variant_by_index(bgendata, snpset[i])
                minor_allele_dosage!(bgendata, variant;
                    T = Float64, mean_impute = true, data = dosageholder,
                    decompressed = decompressed)
                @views copyto!(Z[:, i], dosageholder[bgenrowinds])
            end
            SnpArrays.makestream(pvalfile, "w") do io
                if all(var(Z, dims = [1]) .== 0)
                    println(io, "All SNPs had 0 variation. No results to report.")
                elseif test == :score
                    ts = OrdinalMultinomialScoreTest(fittednullmodel.model, Z)
                    pval = polrtest(ts)
                    println(io, "The joint pvalue of snps indexed",
                        " at $(snpset) is $pval")
                elseif test == :lrt
                    nulldev = deviance(fittednullmodel.model)
                    Xaug = [fittednullmodel.model.X Z]
                    altmodel = polr(Xaug, fittednullmodel.model.Y,
                        fittednullmodel.model.link, solver,
                        wts = fittednullmodel.model.wts)
                    copyto!(γ̂, 1, altmodel.β, fittednullmodel.model.p + 1, q)
                    l2normeffect = norm(γ̂)
                    pval = ccdf(Chisq(q), nulldev - deviance(altmodel))
                    println(io, "The l2norm of the effect size vector",
                    " is $l2normeffect and joint pvalue of snps indexed",
                    " at $(snpset) is $pval")
                end
            end
        end
    else #analysistype == "gxe"
        isnothing(e) &&
            @error("GxE analysis indicated but not environmental variable keyword argument: `e` set.")

        Xaug = [fittednullmodel.model.X zeros(size(fittednullmodel.mm, 1))]
        Xaug2 = [fittednullmodel.model.X zeros(size(fittednullmodel.mm, 1), 2)] #or get Xaug to point to part of it

        # create array for environmental variable and testing
        envvar = modelmatrix(FormulaTerm(fittednullmodel.mf.f.lhs, Term(Symbol(e))),
                 DataFrame(fittednullmodel.mf.data))
        testvec = Matrix{Float64}(undef, size(envvar))
        snpeffectnull = 0.0
        SnpArrays.makestream(pvalfile, "w") do io
            if test == :score
                println(io, "chr,pos,snpid,varid,maf,hwepval,infoscore,snpeffectnull,pval")
            else
                println(io, "chr,pos,snpid,varid,maf,hwepval,infoscore,",
                    "snpeffectnull,snpeffectfull,GxEeffect,pval")
            end
            for (j, variant) in enumerate(bgen_iterator)
                if !snpmask[j] #skip snp, must read marker still.
                    continue
                end
                minor_allele_dosage!(bgendata, variant;
                    T = Float64, mean_impute = true, data = dosageholder,
                    decompressed = decompressed)
                copyto!(@view(Xaug[:, end]), dosageholder[bgenrowinds])
                hwepval = 9.0
                maf = 9.0
                infoscore = 9.0
                try
                    hwepval = BGEN.hwe(bgendata, variant; rmask = bgenrowmask_UInt16)
                catch nothing
                end
                try
                    maf = BGEN.maf(bgendata, variant; rmask = bgenrowmask_UInt16)
                catch nothing
                end
                try
                    infoscore = BGEN.info_score(bgendata, variant; rmask = bgenrowmask_UInt16)
                catch nothing
                end
                if test == :score
                    copyto!(testvec, @view(Xaug[:, end]) .* envvar)
                    nm = polr(Xaug, fittednullmodel.model.Y,
                        fittednullmodel.model.link, solver, wts = fittednullmodel.model.wts)
                    if (Inf in nm.vcov) || (maf == 0.0) #singular design matrix or 0 MAF
                        snpeffectnull = 0.0
                        pval = 1.0
                    else
                        snpeffectnull = nm.β[end]
                        ts = OrdinalMultinomialScoreTest(nm, testvec)
                        pval = polrtest(ts)
                    end
                    println(io, "$(variant.chrom),$(variant.pos),$(variant.rsid),",
                    "$(variant.varid),$(maf),$(hwepval),$(infoscore),",
                    "$snpeffectnull,$pval")
                elseif test == :lrt
                    γ̂ = 0.0 # effect size for columns being tested
                    copyto!(@view(Xaug2[:, end - 1]), @view(Xaug[:, end]))
                    copyto!(@view(Xaug2[:, end]), @view(Xaug[:, end]) .*
                        envvar)
                    nm = polr(Xaug, fittednullmodel.model.Y,
                    fittednullmodel.model.link, solver, wts = fittednullmodel.model.wts)
                    if (Inf in nm.vcov) || (maf == 0.0) #singular design matrix
                        snpeffectnull = 0.0
                        snpeffectfull = 0.0
                        γ̂ = 0.0
                        pval = 1.0
                    else
                        snpeffectnull = nm.β[end]
                        nulldev = deviance(nm)
                        altmodel = polr(Xaug2, fittednullmodel.model.Y,
                            fittednullmodel.model.link, solver,
                            wts = fittednullmodel.model.wts)
                        γ̂ = altmodel.β[end]
                        snpeffectfull = altmodel.β[end-1]
                        pval = ccdf(Chisq(1), nulldev - deviance(altmodel))
                    end
                    println(io, "$(variant.chrom),$(variant.pos),$(variant.rsid),",
                        "$(variant.varid),$(maf),$(hwepval),$(infoscore),",
                        "$snpeffectnull,$snpeffectfull,$γ̂,$pval")
                end
            end
        end
    end
    return fittednullmodel
end

# VCFTools uses different coding for additive, dominant, recessive models than SnpArrays
modelingdict = Dict(
    Val{1}() => :additive,
    Val{2}() => :dominant,
    Val{3}() => :recessive
    )

max_stars_repo_name: Hua-Zhou/PolrGWAS

max_issues_repo_name: Hua-Zhou/PolrGWAS


text:
isdefined(Base, :__precompile__) && __precompile__()

"""
    NonNegLeastSquares
Nonnegative least squares
"""
module NonNegLeastSquares
    using Distributed, LinearAlgebra
    import SparseArrays

    export nonneg_lsq

    ## Algorithms
    include("nnls.jl")
    using .NNLS: nnls
    include("fnnls.jl")
    include("pivot.jl")
    include("pivot_comb.jl")
    include("pivot_cache.jl")
    include("admm.jl")
    include("coorddesc.jl")

    ## Common interface to algorithms
    include("interface.jl")

    ## Helper functions
    include("cssls.jl") # combinatorial subspace least squares (CSSLS)
end


text:
# integration with LinearAlgebra.jl

CuMatOrAdj{T} = Union{CuMatrix, LinearAlgebra.Adjoint{T, <:CuMatrix{T}}, LinearAlgebra.Transpose{T, <:CuMatrix{T}}}
CuOrAdj{T} = Union{CuVecOrMat, LinearAlgebra.Adjoint{T, <:CuVecOrMat{T}}, LinearAlgebra.Transpose{T, <:CuVecOrMat{T}}}

# matrix division

function Base.:\(_A::CuMatOrAdj, _B::CuOrAdj)
    A, B = copy(_A), copy(_B)
    A, ipiv = CuArrays.CUSOLVER.getrf!(A)
    return CuArrays.CUSOLVER.getrs!('N', A, ipiv, B)
end


text:
module GLHelper

using ModernGL

glUniform(index::GLint, value::Float32) = glUniform1f(index, value)
glUniform(index::GLint, value::Tuple{Float32, Float32}) = glUniform2f(index, value[1], value[2])
glUniform(index::GLint, value::Tuple{Float32, Float32, Float32}) = glUniform3f(index, value[1], value[2], value[3])
glUniform(index::GLint, value::Tuple{Float32, Float32, Float32, Float32}) = glUniform4f(index, value[1], value[2], value[3], value[4])

glUniform(index::GLint, value::Int32) = glUniform1i(index, value)
glUniform(index::GLint, value::Tuple{Int32, Int32}) = glUniform2i(index, value[1], value[2])
glUniform(index::GLint, value::Tuple{Int32, Int32, Int32}) = glUniform3i(index, value[1], value[2], value[3])
glUniform(index::GLint, value::Tuple{Int32, Int32, Int32, Int32}) = glUniform4i(index, value[1], value[2], value[3], value[4])

glUniform(index::GLint, value::UInt32) = glUniform1ui(index, value)
glUniform(index::GLint, value::Tuple{UInt32, UInt32}) = glUniform2ui(index, value[1], value[2])
glUniform(index::GLint, value::Tuple{UInt32, UInt32, UInt32}) = glUniform3ui(index, value[1], value[2], value[3])
glUniform(index::GLint, value::Tuple{UInt32, UInt32, UInt32, UInt32}) = glUniform4ui(index, value[1], value[2], value[3], value[4])

nofunc(x...) = nothing

macro notranspose(func)
	quote
		function (location::GLint, count::GLsizei, value::Ptr{GLfloat})
			$func(location, count, GL_FALSE, value)
		end
	end
end

const uniformFuncTable =
	[glUniform1fv		 nofunc							 nofunc							 nofunc							;
	 glUniform2fv		 @notranspose(glUniformMatrix2fv)   @notranspose(glUniformMatrix3x2fv) @notranspose(glUniformMatrix4x2fv);
	 glUniform3fv		 @notranspose(glUniformMatrix2x3fv) @notranspose(glUniformMatrix3fv)   @notranspose(glUniformMatrix4x3fv);
	 glUniform4fv		 @notranspose(glUniformMatrix2x4fv) @notranspose(glUniformMatrix3x4fv) @notranspose(glUniformMatrix4fv)  ]

function get_uniform_func_float(rows::Int, cols::Int)
	func = uniformFuncTable[rows, cols]
	@assert func != nofunc
	return func
end

function glUniform(index::GLint, value::Vector{Float32})
	glUniformFunc = get_uniform_func_float(length(value), 1)
	glUniformFunc(index, one(GLsizei), convert(Ptr{Float32}, value))
end

function glUniform(index::GLint, value::Matrix{Float32})
	rows, cols = size(value)
	glUniformMatrixFunc = get_uniform_func_float(rows, cols)
	glUniformMatrixFunc(index, one(GLsizei), convert(Ptr{Float32}, value))
end

function glUniform(index::GLint, value::Array{Float32, 3})
	rows, cols, elements = size(value)
	local glUniformFunc
	glUniformFunc = get_uniform_func_float(rows, cols)
	glUniformFunc(index, convert(GLsizei, elements), convert(Ptr{Float32}, value))
end

# todo: add glUniform for Vector{Int32 & UInt32}

function gl_get_current_program()
	currentProgram = GLint[0]
	glGetIntegerv(GL_CURRENT_PROGRAM, currentProgram)
	return currentProgram[1]
end

function gl_clear_buffers(color, depth, stencil)
	mask = 0
	if color !== nothing
		glClearColor(color...)
		mask |= GL_COLOR_BUFFER_BIT
	end
	if depth !== nothing
		glClearDepthf(depth)
		mask |= GL_DEPTH_BUFFER_BIT
	end
	if stencil !== nothing
		glClearStencil(stencil)
		mask |= GL_STENCIL_BUFFER_BIT
	end
	glClear(mask)
end

function gl_info()
	vendor = unsafe_string(glGetString(GL_VENDOR))
	version = unsafe_string(glGetString(GL_VERSION))
	renderer = unsafe_string(glGetString(GL_RENDERER))
	@info("OpenGL Vendor: $vendor, Version: $version, Renderer: $renderer")
end

export glUniform, gl_get_current_program, gl_clear_buffers, gl_info

end


text:
using Distributions

mutable struct Lineage
  fitness::Float64
  size::Int64
  state::Vector{Float64}
  mut_rate::Float64
end

function average_fitness(population::Dict{Array{Float64,1}, Lineage})
#calculate average fitness of the population. also outputs population size
  popN = Int64
  popN = 0
  popW = Float64
  popW = 0
  for line in values(population)
    popN+=line.size
    popW+=line.size * line.fitness
  end
  popW = popW/popN
  return popW, popN
end

function assay_mutation_rate(population::Dict{Array{Float64,1}, Lineage})
#calculates mutator frequency
  popN = Int64
  popN = 0
  mutN = Int64
  mutN = 0
  for line in values(population)
    popN+=line.size
    if line.mut_rate>1 mutN+=line.size end
  end
  mutF = Float64
  mutF = mutN/popN
  return mutF
end

function avail_sites(state::Array{Float64,1})
  sites = Int64[]
  for i = 1:length(state)
    if state[i] == 0 push!(sites,i) end
  end
  return sites
end

function mutate_population(population::Dict{Array{Float64,1}, Lineage}, Ub::Float64, sb::Float64, Ud::Float64, sd::Float64, mut_multiplier::Float64)
  new_population = Dict{Array{Float64,1}, Lineage}()
  #initialize new population dictionary
  for (key, line) in population
    mutations = rand(Poisson(line.mut_rate * (Ud + Ub) * line.size))
    if mutations > line.size mutations = line.size end
    bmutations = rand(Binomial(mutations, Ub/(Ub+Ud)))
    dmutations = mutations - bmutations

    open_sites =   avail_sites(line.state)
    #generates a list of indexes of all loci that are available for mutation
    mutation_positions = open_sites[rand(1:end,(dmutations+bmutations))] #generates positions of new mutations by randomly picking indexes from open_sites

      #for each new mutation, copies the state, adds mutation and generates the key.
      #check if it is alrady in the new_population. if it is, add 1 individual to it, if it's not - make another lineage with
      #the key and size 1
    c=1
    for i = 1:bmutations
      #assign beneficial mutations
      newstate = copy(line.state)#
      newstate[mutation_positions[c]] = sb
      c+=1
      if newstate in keys(new_population) new_population[newstate].size+=1
      else
        newfit::Float64 = 1 + sum(newstate[2:end])
        new_population[newstate] = Lineage(newfit, 1, newstate, newstate[1])
      end
    end

    for i = 1:dmutations
      #assign deleterious mutations
      newstate = copy(line.state)#
      newstate[mutation_positions[c]] = -sd
      c+=1
      if newstate in keys(new_population) new_population[newstate].size+=1
      else
        newfit::Float64 = 1 + sum(newstate[2:end])
        new_population[newstate] = Lineage(newfit, 1, newstate, newstate[1])
      end
    end

    new_size = line.size - bmutations - dmutations
    if  new_size > 0
      if key in keys(new_population) new_population[key].size+=new_size
      else new_population[key] = Lineage(line.fitness, new_size, line.state, line.mut_rate)
      end
    end
  end
  return new_population
end

function wright_fisher_reproduction(population::Dict{Array{Float64,1}, Lineage}, N0::Int64, Nnew::Int64, popw::Float64)
  proby_list = Float64[] #array of probabilities
  lineage_list = Lineage[] #initialize lineage liste
  for (key,line) in population
    push!(proby_list, max(line.size/N0 * line.fitness/popw, 0.0))
    #representation of a lineage in the next generations depends on its frequency and relative fitness
    push!(lineage_list, line)
  end
  new_counts_list = rand(Multinomial(Nnew, proby_list))
  new_population = Dict{Array{Float64,1}, Lineage}()
  for i = 1:length(new_counts_list)
    if new_counts_list[i] > 0
      new_population[lineage_list[i].state] = Lineage(lineage_list[i].fitness, new_counts_list[i], lineage_list[i].state, lineage_list[i].mut_rate)
    end
  end
  return new_population
end

function simulate(n0,job_id)
  pop_Ni = 10000000
  init_mut_N = n0
  sb = 0.1
  sd = 0.1
  mutator_strength = 100.0
  Ub = 0.000001
  Ud = 0.0001

  population = Dict{Array{Float64,1}, Lineage}()
  wt_state = zeros(Float64,100)
  wt_state[1] = 1.0
  population[wt_state] = Lineage(1,pop_Ni-init_mut_N,wt_state,1.0)
  mut_state = copy(wt_state)
  mut_state[1] = mutator_strength
  population[mut_state] = Lineage(1,init_mut_N,mut_state,mutator_strength)
  mut_f = assay_mutation_rate(population)
  popw, popN = average_fitness(population)
  generations = Int64
  generations = 0
  #initialize generations counter
  mut_trajectory = Float64[]
  push!(mut_trajectory, mut_f)

  while 0.0<mut_f<1.0
    popw, pop_N = average_fitness(population)
    population = wright_fisher_reproduction(population, pop_N, pop_Ni, popw)
    population = mutate_population(population, Ub, sb, Ud, sd, mutator_strength)
    mut_f = assay_mutation_rate(population)
    popw, pop_N = average_fitness(population)
    push!(mut_trajectory, mut_f)
    generations+=1
    #println("At the end of generation: ", generations, " Mut frequency: ", mut_f, " Fitness: ", popw, " Population size: ", pop_N, " Number of lines ", length(population))

  end
  outfile = open(string(job_id,"mutator_traces.csv"), "a")
  write(outfile, join(mut_trajectory, ","), "\n")
  close(outfile)
  return mut_f ==1.0, generations, popw
end

job_id = ARGS[1]
for n0 in [1, 3, 10, 30, 100, 300, 1000, 3000, 10000, 30000]
  time_to_mut= Float64[]
  #time_to_mut list records all times of successful mutator hitchhiking
  for run = 1:25
    output = simulate(n0, job_id)
    if output[1]
      push!(time_to_mut, output[2])
    else
      push!(time_to_mut, 0)
    end
  end
  outfile = open(string(job_id,"time_to_mut.csv"), "a")
  write(outfile, join(time_to_mut, ","), "\n")
  close(outfile)
end