text:
#############################################################################
#Copyright (c) 2016-2019 Charles Le Losq
#
#The MIT License (MIT)
#
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the #Software without restriction, including without limitation the rights to use, copy, #modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, #and to permit persons to whom the Software is furnished to do so, subject to the #following conditions:
#
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#
#############################################################################

VERSION >= v"1.0.0" && __precompile__()

module Spectra

using Dierckx
using LinearAlgebra
using LsqFit
using Polynomials
using PyCall
using Random
using SparseArrays
using StatsBase
using Statistics

# For PyCall modules
const rampy = PyNULL()

function __init__()
	copy!(rampy, pyimport("rampy"))
end

include("integrale.jl")
include("functions.jl")
include("baseline.jl")
include("bootstrap.jl")
include("tlcorrection.jl")
include("deprecated.jl")
include("ml.jl")
include("peakmeasurement.jl")

#From integrale.jl
export trapz, bandarea

#From functions.jl
export poly, gaussiennes, lorentziennes, pseudovoigts, pearson7, normal_dist, xshift_inversion, xshift_direct,xshift_correction, smooth, flipsp, resample

#From baseline.jl
export baseline

#From bootstrap
export bootsample, bootperf

#From tlcorrection
export tlcorrection

#From ml_regressor
export mlregressor, mlexplorer

# From peakmeasurement
export peakmeas

end # module


text:
using Parameters, SparseArrays, Distributed, CSV, JLD2, DataFrames, ProgressMeter, Quadmath, ParallelUtilities, StatsBase, RCall

# Analytical solutions
import Roots: find_zero
import NLsolve: nlsolve
import SpecialFunctions: polygamma, zeta
import PoissonRandom: pois_rand
import Distributions: Binomial, pdf

# Parse data
import GZip: open
import Parsers: parse
import OrderedCollections: OrderedDict
import FastaIO: readfasta
import Random: randstring

# MK-approaches
import LsqFit: curve_fit, confidence_interval
import HypothesisTests: pvalue,FisherExactTest

include("/home/jmurga/.julia/dev/Analytical/src/parameters.jl")
include("/home/jmurga/.julia/dev/Analytical/src/fixations.jl")
include("/home/jmurga/.julia/dev/Analytical/src/polymorphism.jl")
include("/home/jmurga/.julia/dev/Analytical/src/summaryStatistics.jl")
include("/home/jmurga/.julia/dev/Analytical/src/rates.jl")
include("/home/jmurga/.julia/dev/Analytical/src/inferTools.jl")
include("/home/jmurga/.julia/dev/Analytical/src/readFasta.jl")
include("/home/jmurga/.julia/dev/Analytical/src/methods.jl")

function solve(param::parameters,convoluted_samples::binomial_dict,alTot::Float64,alLow::Float64,gH::Int64,gL::Int64,gamNeg::Int64,afac::Float64,θ::Float64,ρ::Float64)
        # Creating model to solve
    # Γ distribution
    param.al = afac; param.be = abs(afac/gamNeg); param.gamNeg = gamNeg
    # α, αW
    param.alLow = alLow; param.alTot = alTot;
    # Positive selection coefficients
    param.gH    = gH;param.gL = gL
    # Mutation rate and recomb
    param.thetaMidNeutral = θ; param.rho = ρ
    #=0.001.thetaMidNeutral = θ; param.θᵣ .= θᵣ; param.rho = ρ=#
    # Solving θ on non-coding region and probabilites to get α value without BGS
    param.B = 0.999
    setThetaF!(param)
    setPpos!(param)

    return([param.pposL,param.pposH,param.thetaF])

end

function r_old(param::parameters,convoluted_samples::binomial_dict,alTot::Float64,alLow::Float64,gH::Int64,gL::Int64,gamNeg::Int64,afac::Float64,θ::Float64,ρ::Float64)

    # Creating model to solve
    # Γ distribution
    param.al    = afac; param.be = abs(afac/gamNeg); param.gamNeg = gamNeg
    # α, αW
    param.alLow = alLow; param.alTot = alTot;
    # Positive selection coefficients
    param.gH    = gH;param.gL = gL
    # Mutation rate and recomb
    param.thetaMidNeutral = θ; param.rho = ρ
    #=0.001.thetaMidNeutral = θ; param.θᵣ .= θᵣ; param.rho = ρ=#
    # Solving θ on non-coding region and probabilites to get α value without BGS
    param.B = 0.999
    setThetaF!(param)
    setPpos!(param)

    # Allocate array to solve the model for all B values
    r = zeros(size(param.B_bins,1),(size(param.dac,1) * 2) + 14)
    for j in eachindex(param.B_bins)
        # Set B value
        param.B = param.B_bins[j]

        cnvBinom = convoluted_samples.bn[param.B]
        # Solve θ non-coding for the B value.
        setThetaF!(param)
        # Solve model for the B value
        tmp = try
            # Fixation
            fN       = param.B*fixNeut(param)
            fNeg     = param.B*fixNegB(param,0.5*param.pposH+0.5*param.pposL)
            fPosL    = fixPosSim(param,param.gL,0.5*param.pposL)
            fPosH    = fixPosSim(param,param.gH,0.5*param.pposH)

            ds       = fN
            dn       = fNeg + fPosL + fPosH

            # Polymorphism
            neut = DiscSFSNeutDown(param,cnvBinom)
            selH = if isinf(exp(param.gH * 2))
                DiscSFSSelPosDownArb(param,param.gH,param.pposH,cnvBinom)
            else
                DiscSFSSelPosDown(param,param.gH,param.pposH,cnvBinom)
            end
            selL = DiscSFSSelPosDown(param,param.gL,param.pposL,cnvBinom)
            selN = DiscSFSSelNegDown(param,param.pposH+param.pposL,cnvBinom)

            # Cumulative rates
            tmp = cumulative_sfs(hcat(neut,selH,selL,selN),false)
            splitColumns(matrix::Array{Float64,2}) = (view(matrix, :, i) for i in 1:size(matrix, 2));
            neut, selH, selL, selN = splitColumns(tmp)
            sel = (selH+selL)+selN

            ##########
            # Output #
            ##########
           tmp = vcat(param.B,param.alLow,param.alTot,param.gamNeg,param.gL,param.gH,param.al,param.thetaMidNeutral,neut[param.dac],sel[param.dac],ds,dn,fPosL,fPosH,neut[1],sel[1])'
        catch e
            tmp = zeros(32)
        end
        @inbounds r[j,:] = tmp
        #=@inbounds r = tmp=#
    end

    return r
end

function r(param::parameters,convoluted_samples::binomial_dict,alTot::Float64,alLow::Float64,gH::Int64,gL::Int64,gamNeg::Int64,afac::Float64,θ::Float64,ρ::Float64,e::Vector{Float64})

    # Creating model to solve
    # Γ distribution
    param.al    = afac; param.be = abs(afac/gamNeg); param.gamNeg = gamNeg
    # α, αW
    param.alLow = alLow; param.alTot = alTot;
    # Positive selection coefficients
    param.gH    = gH;param.gL = gL
    # Mutation rate and recomb
    param.thetaF = e[3]; param.thetaMidNeutral = θ; param.rho = ρ
    # Fixation probabilites
    param.pposL = e[1]; param.pposH = e[2]

    # Solving θ on non-coding region and probabilites to get α value without BGS
    param.B = e[4]
    cnvBinom = convoluted_samples.bn[e[4]]
    # Solve θ non-coding for the B value.
    setThetaF!(param)
    # Solve model for the B value

    # Fixation
    fN       = param.B*fixNeut(param)
    fNeg     = param.B*fixNegB(param,0.5*param.pposH+0.5*param.pposL)
    fPosL    = fixPosSim(param,param.gL,0.5*param.pposL)
    fPosH    = fixPosSim(param,param.gH,0.5*param.pposH)

    ds       = fN
    dn       = fNeg + fPosL + fPosH

    # Polymorphism
    neut = DiscSFSNeutDown(param,cnvBinom)
    selH = if isinf(exp(param.gH * 2))
        DiscSFSSelPosDownArb(param,param.gH,param.pposH,cnvBinom)
    else
        DiscSFSSelPosDown(param,param.gH,param.pposH,cnvBinom)
    end
    selL = DiscSFSSelPosDown(param,param.gL,param.pposL,cnvBinom)
    selN = DiscSFSSelNegDown(param,param.pposH+param.pposL,cnvBinom)

    # Cumulative rates
    tmp = cumulative_sfs(hcat(neut,selH,selL,selN),false)
    splitColumns(matrix::Array{Float64,2}) = (view(matrix, :, i) for i in 1:size(matrix, 2));
    neut, selH, selL, selN = splitColumns(tmp)
    sel = (selH+selL)+selN

    ##########
    # Output #
    ##########
    r = vcat(param.B,param.alLow,param.alTot,param.gamNeg,param.gL,param.gH,param.al,param.thetaMidNeutral,neut[param.dac],sel[param.dac],ds,dn,fPosL,fPosH,neut[1],sel[1])'

    return r
end

param = parameters()
convoluted_samples = binomial_dict()
binomOp!(param,convoluted_samples.bn)

#gH=200:1000;gL=1:10;gamNeg=-2000:-200;iterations = 10^2;shape=adap.al;theta=0.001;rho=0.001

iterations = 10^3
afac       = fill(param.al,iterations)

# Random α values
nTot       = fill(0.4,iterations)
# Random weak selection coefficients
nLow       = fill(0.2,iterations)
# Creating N models to iter in threads. Set N models (paramerters) and sampling probabilites (binomial_dict)
nParam     = [param for i in 1:iterations];
nBinom     = [convoluted_samples for i in 1:iterations];
# Random strong selection coefficients
ngh        = fill(200,iterations)
ngl        = fill(5,iterations)
# Random negative selection coefficients
ngamNeg    = fill(-457,iterations)
θ          = fill(0.001,iterations)
ρ          = fill(0.001,iterations)

function rates_threads(;param::parameters,
                convoluted_samples::binomial_dict,
                gH::S,
                gL::S,
                gamNeg::S,
                theta::Union{Float64,Nothing}=0.001,
                rho::Union{Float64,Nothing}=0.001,
                shape::Float64=0.184,
                iterations::Int64,
                output::String) where S <: Union{Array{Int64,1},UnitRange{Int64},Nothing}

    # Iterations = models to solve
    # Factor to modify input Γ(shape) parameter. Flexible Γ distribution over negative alleles
    fac     = rand(-2:0.05:2,iterations)
    afac    = @. param.al*(2^fac)

    # Deleting shape > 1. Negative alpha_x values
    idx = findall(afac .> 1)
    if !isempty(idx)
        afac[idx] = rand(afac[afac .< 1],size(idx,1))
    end

    # Random α values
    nTot    = rand(0.1:0.01:0.9,iterations)

    # Defining αW. It is possible to solve non-accounting for weak fixations
    if isnothing(gL)
        # Setting αW to 0 for all estimations
        nLow    = fill(0.0,iterations)
        # Random strong selection coefficients
        ngl     = rand(repeat([1],iterations),iterations);
    else
        # Setting αW as proportion of α
        lfac    = rand(0.0:0.05:0.9,iterations)
        nLow    = @. nTot * lfac
        # Random weak selection coefficients
        ngl     = rand(repeat(gL,iterations),iterations);
    end

    # Creating N models to iter in threads. Set N models (paramerters) and sampling probabilites (binomial_dict)
    nParam  = [param for i in 1:iterations];
    nBinom  = [convoluted_samples for i in 1:iterations];

    # Random strong selection coefficients
    ngh     = rand(repeat(gH,iterations),iterations);
    # Random negative selection coefficients
    ngamNeg = rand(repeat(gamNeg,iterations),iterations);

    # Random θ on coding regions
    #=obsNeut        = sfs.p0./sum(sfs.p0)
    n              = 1 ./collect(1:(param.nn-1))
    expNeut        = n ./ sum(n)=#

    if !isnothing(theta)
        θ = fill(theta,iterations)
        #=θᵣ= fill(theta .* obsNeut ./ expNeut,iterations)=#
    else
        θ = rand(0.0005:0.0005:0.01,iterations)
        #=θᵣ= map(x -> x .* obsNeut ./ expNeut,θ)=#
    end

    # Random ρ on coding regions
    if !isnothing(rho)
        ρ = fill(rho,iterations)
    else
        ρ = rand(0.0005:0.0005:0.05,iterations)
    end

    ######Parallel; 10^5 -> 62.76"
    x = zeros(1,3,iterations);
    @inbounds @sync for i in 1:iterations
        Base.Threads.@spawn begin
            @inbounds x[:,:,i] = solve(param,convoluted_samples, nTot[i], nLow[i], ngh[i], ngl[i], ngamNeg[i], afac[i], θ[i], ρ[i])
        end
    end

    #y = ThreadsX.map(solve, nParam,nBinom, nTot, nLow, ngh, ngl, ngamNeg, afac, θ, ρ)

    x = repeat(x, outer = [size(param.B_bins,1), 1, 1]);
    x = vcat([ @view x[:,:,i] for i=1:iterations]...);
    b = repeat(param.B_bins,iterations);
    e = hcat(x,b);

    out = zeros(iterations,32,size(param.B_bins,1));

    @showprogress for (j,val) in enumerate(reverse(param.B_bins))
        tmp = e[e[:,4] .== val,:]
        for i in 1:iterations
            Base.Threads.@spawn begin
               out[i,:,j] = r(nParam[i],nBinom[i], nTot[i], nLow[i], ngh[i], ngl[i], ngamNeg[i], afac[i], θ[i], ρ[i],tmp[i,:]);
            end
        end
    end

    df = vcat([@view out[:,:,i] for i=1:size(out,3)]...);

    # Saving models and rates
    models = DataFrame(df[:,1:8],[:B,:alLow,:alTot,:gamNeg,:gL,:gH,:al,:ρ]);
    neut   = df[:,9:(8+size(param.dac,1))];
    sel    = df[:,(9+size(param.dac,1)):(8+size(param.dac,1)*2)];
    dsdn   = Array(df[:,(end-5):end-2]);

    sum_pol = df[:,end-1:end];

    # Saving multiple summary statistics
    n = OrderedDict{Int,Array}();
    s = OrderedDict{Int,Array}();
    for i in eachindex(param.dac)
        n[param.dac[i]] = neut[:,i]
        s[param.dac[i]] = sel[:,i]
    end

    #=return models=#
    # Writting HDF5 file
    JLD2.jldopen(output, "a+") do file
        file[string(param.N)* "/" * string(param.n) * "/models"] = models;
        file[string(param.N)* "/" * string(param.n) * "/neut"]   = n;
        file[string(param.N)* "/" * string(param.n) * "/sel"]    = s;
        file[string(param.N)* "/" * string(param.n) * "/dsdn"]   = dsdn;
        file[string(param.N)* "/" * string(param.n) * "/pol"]   = sum_pol;
        file[string(param.N)* "/" * string(param.n) * "/dac"]    = param.dac;
    end;
end

@time df = rates_threads(param = param,convoluted_samples=convoluted_samples,gH=200:2000,gL=1:10,gamNeg=-2000:-200,iterations = 10^3,shape=adap.al,output="/home/jmurga/rates_threads.jld2");rm("/home/jmurga/rates_threads.jld2")

@time df = Analytical.rates(param = param,convoluted_samples=convoluted_samples,gH=200:2000,gL=1:10,gamNeg=-2000:-200,iterations = 10^2,shape=adap.al,output="/home/jmurga/rates_distributed.jld2");rm("/home/jmurga/rates_distributed.jld2")

######Sequential

real = fill(Matrix{Float64}(undef, 0,0),iterations);
@time for i in eachindex(nTot)
    real[i] = r(param,convoluted_samples, nTot[i], nLow[i], ngh[i], ngl[i], ngamNeg[i], afac[i], θ[i], ρ[i])
end

### ThreadsPools
#=@time @sync begin @qbthreads for i in 1:iterations
        x[:,:,i] = solve(nParam[i],nBinom[i], nTot[i], nLow[i], ngh[i], ngl[i], ngamNeg[i], afac[i], θ[i], ρ[i])
    end
end=#
#=@showprogress for (j,val) in enumerate(param.B_bins)

    tmp = e[e[:,4] .== val,:]
    @sync begin @qbthreads for i in 1:iterations
           out[i,:,j] = r(nParam[i],nBinom[i], nTot[i], nLow[i], ngh[i], ngl[i], ngamNeg[i], afac[i], θ[i], ρ[i],tmp[i,:])
        end
    end
end
=#

ext: jl


text:
using JLD
using MAT
using JSON
using Knet

include("tokenizer.jl")

global vocab = Dict{String, Int}()
caption_dicts = Array{Array{Tuple{Tuple{Int64,Array{String,1}},Int64},1},1}()
Tokenizer.tokenize(vocab, caption_dicts)

dataset = open("./data/Flickr30k/karpathy/dataset.json");
data = JSON.parse(dataset);
images = data["images"];
file = matopen("./data/Flickr30k/karpathy/vgg_feats.mat");
features = read(file,"feats");

feats = Dict{Int,Array{Float32}}()

# for image in images
#     feature = features[:,image["imgid"]+1];
#     filename = split(image["filename"],'.')[1];
#     save("./data/Flickr30k/karpathy/features/$(filename).jld", "feature", feature);
# end

for image in images
     feature = features[:,image["imgid"]+1];
     filename = split(image["filename"],'.')[1];
     get!(feats,parse(Int,filename),feature);
end

println(length(feats))
count = 1;
for (el,v) in caption_dicts[1]
     id = el[1];
     filename = "./data/Flickr30k/karpathy/features/$(id).jld";
     if isfile(filename)
         exist = get(feats,id,nothing);
         if exist == nothing
           println(id)
           get!(feats,id,load(filename,"feature"))
         end
    end
    count += 1;
end

println(length(feats))

save("./data/Flickr30k/karpathy/features/feats.jld", "features", feats);

max_issues_repo_name: ekinakyurek/LRCN

max_forks_repo_name: ekinakyurek/LRCN


text:
import ChainRulesCore: frule, rrule
using LinearAlgebra

function compute_grad(ld::QuadLimbDark, b::S, r) where S
    T = float(S)
    bcut = 1e-3
    dfdg1 = zero(T)
    dfdg2 = zero(T)
    dfdg3 = zero(T)

    ## check for trivial cases
    if b ≥ 1 + r || iszero(r)
        # completely unobscured
        return one(T), SA[dfdg1, dfdg2, dfdg3], zero(T), zero(T)
    elseif r ≥ 1 + b
        # completely obscured
        return zero(T), SA[dfdg1, dfdg2, dfdg3], zero(T), zero(T)
    end

    r2 = r^2
    b2 = b^2

    if iszero(b)
        onemr2 = 1 - r2
        sqrt1mr2 = sqrt(onemr2)

        # annular ellipse
        facd = -2 * r
        if ld.n_max == 0
            flux = ld.g_n[begin] * onemr2
            dfdr = ld.g_n[begin] * facd
        else
            flux = ld.g_n[begin] * onemr2 + 2 / 3 * ld.g_n[begin + 1] * sqrt1mr2^3
            dfdr = ld.g_n[begin] * facd + ld.g_n[begin + 1] * facd * sqrt1mr2
        end
        fac = 2 * r2 * onemr2
        if ld.n_max > 1
            flux -= ld.g_n[begin + 2] * fac
            dfdr += ld.g_n[begin + 2] * facd * 2 * (onemr2 - r2)
        end
        dfdr *= π * ld.norm
        dfdg = SA[(onemr2 - flux) * π * ld.norm,
                  2 / 3 * (sqrt1mr2^3 - flux) * π * ld.norm,
                  -fac]
        return flux * π * ld.norm, dfdg * ld.norm, zero(T), dfdr
    end

    # take a moment to calculate values used repeatedly in the
    # numerical routines ahead
    onembmr2 = (r - b + 1) * (1 - r + b)
    onembmr2inv = inv(onembmr2)
    sqonembmr2 = sqrt(onembmr2)
    br = b * r
    fourbr = 4 * br
    fourbrinv = inv(fourbr)
    sqbr = sqrt(br)
    sqbrinv = inv(sqbr)
    onembpr2 = (1 - r - b) * (1 + b + r)
    sqarea = sqarea_triangle(one(T), r, b)
    k2 = max(zero(T), onembmr2 * fourbrinv)
    k = sqrt(k2)
    omemr2 = (1 - r) * (1 + r)
    onemr2mb2 = omemr2 - b2
    onemr2pb2 = omemr2 + b2
    if k2 ≥ 1
        if k2 > 2
            kc2 = 1 - inv(k2)
        else
            kc2 = onembpr2 * onembmr2inv
        end
    else
        if k2 > 0.5
            kc2 = (r - 1 + b) * (b + r + 1) * fourbrinv
        else
            kc2 = 1 - k2
        end
    end
    kc = sqrt(abs(kc2))

    ## Compute uniform term
    (s0, kap0, kite_area2, kck), ∇s0 = compute_uniform_grad(b, r; b2, r2, sqarea, fourbrinv)

    flux = ld.g_n[begin] * s0
    ∇flux = ld.g_n[begin] * ∇s0
    dfdg1 = s0

    if ld.n_max == 0
        dfdg1 -= flux * ld.norm * π
        dfdb, dfdr = ∇flux * ld.norm
        dfdg = SA[dfdg1, dfdg2, dfdg3]
        return flux * ld.norm, dfdg * ld.norm, dfdb, dfdr
    end

    ## compute linear term
    (s1, Eofk, Em1mKdm), ∇s1 = compute_linear_grad(b, r; k2, kc, kc2, r2, b2, br, fourbr, sqbr, onembmr2, onembpr2, onembmr2inv, sqonembmr2)

    flux += ld.g_n[begin + 1] * s1
    ∇flux += ld.g_n[begin + 1] * ∇s1
    dfdg2 = s1

    if ld.n_max == 1
        dfdg1 -= flux * ld.norm * π
        dfdg2 -= flux * ld.norm * π * 2 / 3
        dfdb, dfdr = ∇flux * ld.norm
        dfdg = SA[dfdg1, dfdg2, dfdg3]
        return flux * ld.norm, dfdg * ld.norm, dfdb, dfdr
    end

    ## calculate quadratic term
    s2, ∇s2 = compute_quadratic_grad(b, r; s0, r2, b2, kap0, kite_area2, k2, ∇s0)
    flux += ld.g_n[begin + 2] * s2
    ∇flux += ld.g_n[begin + 2] * ∇s2
    dfdg3 = s2

    dfdg1 -= flux * ld.norm * π
    dfdg2 -= flux * ld.norm * π * 2 / 3
    dfdb, dfdr = ∇flux * ld.norm
    dfdg = SA[dfdg1, dfdg2, dfdg3]
    return flux * ld.norm, dfdg * ld.norm, dfdb, dfdr
end

####

function frule((_, Δld, Δb, Δr), ::typeof(compute), ld::LD, b, r) where {LD <: QuadLimbDark}
    f, dfdg, dfdb, dfdr = compute_grad(ld, b, r)
    ∂g_n = dot(dfdg, Δld.g_n)
    return f, ∂g_n + dfdb * Δb + dfdr * Δr
end

function rrule(::typeof(compute), ld::LD, b, r) where {LD <: QuadLimbDark}
    f, dfdg, dfdb, dfdr = compute_grad(ld, b, r)
    function compute_pullback(Δf)
        ∂ld = Tangent{LD}(g_n=dfdg * Δf)
        ∂b = dfdb * Δf
        ∂r = dfdr * Δf
        return NoTangent(), ∂ld, ∂b, ∂r
    end
    return f, compute_pullback
end

function frule((_, Δu_n), ::Type{<:QuadLimbDark}, u_n::AbstractVector{T}) where T
    Ω = QuadLimbDark(u_n)
    ∇g_n = SA[-one(T) -1.5
               one(T)  2.0
              zero(T) -0.25]

    # calculate flux normalization factor, which only depends on first two terms
    N = length(u_n)
    if N == 0
        ∂g_n = @SVector zeros(T, 3)
    elseif N == 1
        ∂g_n = SA[-Δu_n[begin], Δu_n[begin], zero(T)]
    else
        Δu_n_full = SA[Δu_n[begin], Δu_n[begin + 1]]
        ∂g_n = ∇g_n * Δu_n_full
    end
    ∂Ω = Tangent{typeof(Ω)}(g_n=∂g_n)
    return Ω, ∂Ω
end

function rrule(::Type{<:QuadLimbDark}, u_n::AbstractVector{T}; maxiter=100) where T
    Ω = QuadLimbDark(u_n)
    ∇g_n = SA[-one(T)  one(T) zero(T)
                 -1.5     2.0  -0.25]
    N = length(u_n)
    proj = ProjectTo(u_n)
    function QuadLimbDark_pullback(Δld)
        ∂u = @view(∇g_n[begin:begin+N-1, :]) * Δld.g_n
        return NoTangent(), proj(∂u)
    end
    return Ω, QuadLimbDark_pullback
end


text:
using Test

@testset "Testing Lab 1" begin
   include("test1.jl")
   include("test2.jl")
end;


text:
using SafeTestsets

@safetestset "just some tests" begin include("basic_tests.jl") end


text:
# This file includes portions from JuliaMath/Combinatorics.jl in modified form.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module Setpartitions
using Nemo

export ModuleSetPartitions
export SetPartitions, SetNumber

"""

* SetPartitions(s::AbstractVector)
Return an iterator over all set partitions of the elements of the array ``s``, represented as arrays of arrays.

* SetPartitions(n::Int)
Return an iterator over all set partitions of the elements of the array ``[1,2,...,n]``, represented as arrays of arrays.

* SetPartitions(s::AbstractVector, m::Int)
Return all set partitions of the elements of the array ``s`` into exactly ``m`` subsets, represented as arrays of arrays.

* SetPartitions(n::Int, m::Int)
Return all set partitions of the elements of the array ``{1,2,3,...,n}`` into exactly ``m`` subsets, represented as arrays of arrays.

* SetNumber(n::Int)
Return the numbers of partitions of an ``n``-set into nonempty subsets.

* SetNumber(n::Int, m::Int)
Return the numbers of partitions of an ``n``-set into ``m`` nonempty subsets.
"""
const ModuleSetPartitions = ""

struct SetIterator{T<:AbstractVector}
    s::T
end

Base.length(p::SetIterator) = SetNumber(length(p.s))
Base.eltype(p::SetIterator) = Vector{Vector{eltype(p.s)}}

function Base.iterate(p::SetIterator)
    n = length(p.s)
    iterate(p, (zeros(Int, n), ones(Int, n - 1), n, 1))
end

function Base.iterate(p::SetIterator, s)
    s[1][1] > 0 && return
    NextSetPartition(p.s, s...)
end

# Algorithm H from TAOCP 7.2.1.5
function NextSetPartition(s::AbstractVector, a, b, n, m)
    function makeparts(s, a, m)
        temp = [similar(s, 0) for k = 0:m]
        for i = 1:n
            push!(temp[a[i]+1], s[i])
        end
        filter!(!isempty, temp)
    end

    if isempty(s)
        return ([s], (eltype(a)[1], eltype(b)[], n, 1))
    end

    part = makeparts(s, a, m)

    if a[end] != m
        a[end] += 1
    else
        j = 0
        for jj = n - 1:-1:1
            j = jj
            a[jj] == b[jj] || break
        end
        a[j] += 1
        m = Int(b[j]) + (a[j] == b[j])
        for k = j+1:n-1
            a[k] = 0
            b[k] = m
        end
        a[end] = 0
    end

    return (part, (a, b, n, m))
end

"""

Return an iterator over all set partitions of the elements of an array ``s``, represented as arrays of arrays.
"""
function SetPartitions(s::AbstractVector)
    # The only partition of the empty set is the set that contains the empty set.
    isempty(s) && return ([[]])
    SetIterator(s)
end

"""

Return an iterator over all set partitions of the elements of the array ``[1,2,...,n]``, represented as arrays of arrays.
"""
function SetPartitions(n::Int)
    n < 0 && throw(DomainError(n), "n must be nonnegative")
    n == 0 && return ([[]])
    SetIterator(1:n)
end

#######################################################################

struct KSetIterator{T<:AbstractVector}
    s::T
    m::Int
end

Base.length(p::KSetIterator) = SetNumber(length(p.s), p.m)
Base.eltype(p::KSetIterator) = Vector{Vector{eltype(p.s)}}

"""

Return all set partitions of the elements of an array ``s`` into exactly ``m`` subsets, represented as arrays of arrays.
"""
function SetPartitions(s::AbstractVector, m::Int)
    # The only partition of the empty set is the set that contains the empty set.
    if length(s) == 0 && m == 0 return [[]] end
    if m == 0 || length(s) < m  return [] end
    if (length(s) >= 1 && m >= 1) return KSetIterator(s, m) end
    throw(DomainError((length(s), m), "length(s) and m must be nonnegative"))
end

# State consists of:
# vector a of length n describing to which partition every element of s belongs
# vector b of length n describing the first index b[i] that belongs to partition i
# integer n
function Base.iterate(p::KSetIterator)
    n = length(p.s)
    m = p.m
    state = m <= n ?
            (vcat(ones(Int, n - m), 1:m), vcat(1:1, n-m+2:n), n) :
            (Int[], Int[], n)

    iterate(p, state)
end

function Base.iterate(p::KSetIterator, s)
    (isempty(s[1]) || s[1][1] > 1) && return
    NextKSetPartition(p.s, p.m, s...)
end

function NextKSetPartition(s::AbstractVector, m, a, b, n)
    function makeparts(s, a)
        local part = [similar(s, 0) for k = 1:m]
        for i = 1:n
            push!(part[a[i]], s[i])
        end
        return part
    end

    part = makeparts(s, a)

    if m == 1
        a[1] = 2
        return (part, (a, b, n))
    end

    if a[end] != m
        a[end] += 1
    else
        j = k = 0
        for jj = n - 1:-1:1
            j = jj
            if a[j] < m && b[a[j]+1] < j
                break
            end
        end
        if j > 1
            a[j] += 1
            for p = j+1:n
                if b[a[p]] != p
                    a[p] = 1
                end
            end
        else
            for kk = m:-1:2
                k = kk
                if b[k-1] < b[k] - 1
                    break
                end
            end
            b[k] -= 1
            b[k+1:m] = n - m + k + 1:n
            a[1:n] .= 1
            a[b] = 1:m
        end
    end

    return (part, (a, b, n))
end

"""

Return all set partitions of the elements of the array ``{1,2,3,...,n}`` into exactly ``m`` subsets, represented as arrays of arrays.
"""
function SetPartitions(n::Int, m::Int)
    if (n < 0) || (m < 0)
        throw(DomainError(n, m), "n and m must be nonnegative")
    end
    SetPartitions(1:n, m)
end

"""

Return the numbers of partitions of an ``n``-set into ``m`` nonempty subsets.
"""
function SetNumber(n::Int, m::Int)
    numpart = 0
    for k = 0:m
        numpart += (-1)^(m - k) * binomial(m, k) * (k^n)
    end
    numpart = div(numpart, factorial(m))
    return numpart
end

"""

Return the numbers of partitions of an ``n``-set into nonempty subsets.
"""
SetNumber(n::Int) = Int(Nemo.bell(n))

#START-TEST-########################################################

using Test, SeqUtils

function test()
    @testset "SetPartitions" begin
        @test collect(SetPartitions([])) == [[]]

        @test collect(SetPartitions([], 0)) == [[]]
        @test collect(SetPartitions([], 1)) == []

        @test collect(SetPartitions([1])) == [[[1]]]
        @test collect(SetPartitions([1], 0)) == []
        @test collect(SetPartitions([1], 1)) == [[[1]]]
        @test collect(SetPartitions([1], 2)) == []

        @test collect(SetPartitions([1,2])) == [[[1,2]],[[1],[2]]]
        @test collect(SetPartitions([1,2], 0)) == []
        @test collect(SetPartitions([1,2], 1)) == [[[1, 2]]]
        @test collect(SetPartitions([1,2], 2)) == [[[1], [2]]]
        @test collect(SetPartitions([1,2], 3)) == []

        @test collect(SetPartitions(0)) == [[]]

        @test collect(SetPartitions(0, 0)) == [[]]
        @test collect(SetPartitions(0, 1)) == []

        @test collect(SetPartitions(1)) == [[[1]]]
        @test collect(SetPartitions(1, 0)) == []
        @test collect(SetPartitions(1, 1)) == [[[1]]]
        @test collect(SetPartitions(1, 2)) == []

        @test collect(SetPartitions(2)) == [[[1,2]],[[1],[2]]]
        @test collect(SetPartitions(2, 0)) == []
        @test collect(SetPartitions(2, 1)) == [[[1, 2]]]
        @test collect(SetPartitions(2, 2)) == [[[1], [2]]]
        @test collect(SetPartitions(2, 3)) == []

        SP = [ [[1,2,3]], [[1,2], [3]], [[1,3], [2]], [[1], [2,3]], [[1], [2], [3]] ]

        CP = collect(SetPartitions([1, 2, 3]))
        @test SP == CP

        RP = reduce(vcat, [collect(SetPartitions([1, 2, 3], k)) for k ∈ 1:3])
        @test SP == RP

        for n ∈ 1:7
            s = collect(0:n-1)
            @test length(collect(SetPartitions(s))) == SetNumber(length(s))
        end

        for n ∈ 1:7, k ∈ 1:n
            sn = SetNumber(n, k)
            @test length(collect(SetPartitions(1:n, k))) == sn
        end

        @test_throws DomainError SetPartitions([ ], -1)
        @test_throws DomainError SetPartitions([1], -1)
    end
end

function demo()
    println()
    s = SetPartitions([])
    println("The empty set has a partition, namely: ", reduce(vcat, s))
    println()

    for n ∈ 1:4
        println("-- ", n)
        for k ∈ 1:n
            for p ∈ SetPartitions(n, k)
            # same as: for p ∈ SetPartitions(1:n, k)
                Println(p)
            end
        end
        println()
    end

    for n ∈ 1:4
        println("-- ", n)
        for k ∈ 1:n
            for p ∈ SetPartitions(n, k)
                Println(p)
            end
        end
        println()
    end

    println("The Set numbers (aka Stirling2 numbers):")
    for n ∈ 0:6
        println([SetNumber(n, k) for k ∈ 0:n])
    end
    println()
    println("The Set numbers (aka Bell numbers):")
    Println([SetNumber(n) for n ∈ 0:9])
end

function perf() end

function main()
    test()
    demo()
    perf()
end

main()

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
"""
    abstract type AbstractIntervalGrid{T} <: AbstractGrid1d{T}

An AbstractIntervalGrid is a grid that is defined on an interval.
"""
abstract type AbstractIntervalGrid{T} <: AbstractGrid1d{T} end

isperiodic(::AbstractIntervalGrid) = false

"""
    abstract type AbstractEquispacedGrid{T} <: AbstractIntervalGrid{T}

An equispaced grid has equispaced points, and therefore it has a step.
"""
abstract type AbstractEquispacedGrid{T} <: AbstractIntervalGrid{T} end

==(g1::AbstractEquispacedGrid, g2::AbstractEquispacedGrid) =
    covering(g1) == covering(g2) && size(g1)==size(g2) &&
		(g1[1]==g2[1]) && (g1[end]==g2[end])

# Mainly for technical reasons we introduce two intermediate types
abstract type AbstractEquispacedRangeGrid{T} <: AbstractEquispacedGrid{T} end
abstract type AbstractUnitEquispacedGrid{T} <: AbstractEquispacedGrid{T} end

# the range grids are defined in terms of a range
range(g::AbstractEquispacedRangeGrid) = g.range
size(g::AbstractEquispacedRangeGrid) = size(range(g))
step(g::AbstractEquispacedRangeGrid) = step(range(g))

@inline function unsafe_grid_getindex(grid::AbstractEquispacedRangeGrid, i::Int)
	@inbounds getindex(range(grid), i)
end

# the unit equispaced grids are defined on the unit interval [0,1]
covering(g::AbstractUnitEquispacedGrid{T}) where {T} = UnitInterval{T}()
size(g::AbstractUnitEquispacedGrid) = (g.n,)

"""
    struct EquispacedGrid{T} <: AbstractEquispacedGrid{T}

An equispaced grid with n points on an interval [a,b], including the endpoints.
It has step (b-a)/(n-1).

# Example
```jldocs
julia> EquispacedGrid(4,0,1)
4-element EquispacedGrid{Float64}:
 0.0
 0.3333333333333333
 0.6666666666666666
 1.0
```
"""
struct EquispacedGrid{T} <: AbstractEquispacedRangeGrid{T}
    # Use StepRangeLen for higher precision
    range   :: LinRange{T}

    EquispacedGrid{T}(n::Int, a, b) where {T} = new(LinRange(T(a),T(b),n))
end

covering(g::EquispacedGrid) = g[1]..g[end]

"An equispaced grid on the unit interval `[0,1]`."
struct UnitEquispacedGrid{T} <: AbstractUnitEquispacedGrid{T}
	n	::	Int
end

UnitEquispacedGrid(n::Int) = UnitEquispacedGrid{Float64}(n)
similargrid(g::UnitEquispacedGrid, ::Type{T}, n::Int) where {T} =
	UnitEquispacedGrid{T}(n)
step(g::UnitEquispacedGrid{T}) where {T} = T(1)/(g.n-1)

unsafe_grid_getindex(grid::UnitEquispacedGrid{T}, i::Int) where {T} =
	convert(T, (i-1))/(grid.n-1)

canonicalgrid(g::EquispacedGrid{T}) where {T} = UnitEquispacedGrid{T}(length(g))
mapfrom_canonical(g::EquispacedGrid) = mapto(0..1, covering(g))

show(io::IO, g::UnitEquispacedGrid{T}) where T =
	T == Float64 ? print(io, "UnitEquispacedGrid($(length(g)))") : print(io, "UnitEquispacedGrid{$(T)}($(length(g)))")

"""
    struct PeriodicEquispacedGrid{T} <: AbstractEquispacedGrid{T}

A periodic equispaced grid is an equispaced grid that omits the right endpoint.
It has step size (b-a)/n.

# Example
```jldocs
julia> PeriodicEquispacedGrid(4,0,1)
4-element PeriodicEquispacedGrid{Float64}:
 0.0
 0.25
 0.5
 0.75
```
"""
struct PeriodicEquispacedGrid{T} <: AbstractEquispacedRangeGrid{T}
    range   :: LinRange{T}
    a   ::  T
    b   ::  T

    PeriodicEquispacedGrid{T}(n::Int, a, b) where {T} =
		new(LinRange(T(a),T(b),n+1)[1:end-1], a, b)
end

range(g::PeriodicEquispacedGrid) = g.range

covering(grid::PeriodicEquispacedGrid) = grid.a..grid.b
isperiodic(::PeriodicEquispacedGrid) = true

"A periodic equispaced grid on the unit interval `[0,1]`."
struct UnitPeriodicEquispacedGrid{T} <: AbstractUnitEquispacedGrid{T}
	n	::	Int
end

UnitPeriodicEquispacedGrid(n::Int) = UnitPeriodicEquispacedGrid{Float64}(n)
similargrid(g::UnitPeriodicEquispacedGrid, ::Type{T}, n::Int) where {T} =
	UnitPeriodicEquispacedGrid{T}(n)
isperiodic(::UnitPeriodicEquispacedGrid) = true
step(g::UnitPeriodicEquispacedGrid{T}) where {T} = T(1)/g.n

unsafe_grid_getindex(grid::UnitPeriodicEquispacedGrid{T}, i::Int) where {T} =
	convert(T, (i-1))/grid.n

const FourierGrid = UnitPeriodicEquispacedGrid
@deprecate FourierGrid(n, a, b) rescale(FourierGrid(n), a, b)

canonicalgrid(g::PeriodicEquispacedGrid{T}) where {T} =
	UnitPeriodicEquispacedGrid{T}(length(g))
mapfrom_canonical(g::PeriodicEquispacedGrid) = mapto(0..1, covering(g))

show(io::IO, g::UnitPeriodicEquispacedGrid{T}) where T =
	T == Float64 ? print(io, "FourierGrid($(length(g)))") : print(io, "FourierGrid{$(T)}($(length(g)))")

"""
    struct MidpointEquispacedGrid{T} <: AbstractEquispacedGrid{T}

A MidpointEquispaced grid is an equispaced grid with grid points in the centers of the equispaced
subintervals. In other words, this is a DCT-II grid.
It has step `(b-a)/n`.

# Example
```jldocs
julia> MidpointEquispacedGrid(4,0,1)
4-element MidpointEquispacedGrid{Float64}:
 0.125
 0.375
 0.6249999999999999
 0.875
```
"""
struct MidpointEquispacedGrid{T} <: AbstractEquispacedRangeGrid{T}
    range   ::LinRange{T}
    a   ::  T
    b   ::  T

    MidpointEquispacedGrid{T}(n::Int, a, b) where {T} =
		new(LinRange(T(a),T(b),2n+1)[2:2:end], a, b)
end

range(g::MidpointEquispacedGrid) = g.range

covering(grid::MidpointEquispacedGrid) = grid.a..grid.b
isperiodic(::MidpointEquispacedGrid) = true

"A midpoint equispaced grid on the unit interval `[0,1]`."
struct UnitMidpointEquispacedGrid{T} <: AbstractUnitEquispacedGrid{T}
	n	::	Int
end

UnitMidpointEquispacedGrid(n::Int) = UnitMidpointEquispacedGrid{Float64}(n)
similargrid(g::UnitMidpointEquispacedGrid, ::Type{T}, n::Int) where {T} = UnitMidpointEquispacedGrid{T}(n)
step(g::UnitMidpointEquispacedGrid{T}) where {T} = T(2)/(2*g.n-1)

unsafe_grid_getindex(grid::UnitMidpointEquispacedGrid{T}, i::Int) where {T} =
	convert(T, 2(i-1))/(2grid.n-1)

canonicalgrid(g::MidpointEquispacedGrid{T}) where {T} =
	UnitMidpointEquispacedGrid{T}(length(g))
mapfrom_canonical(g::MidpointEquispacedGrid) = mapto(0..1, covering(g))

# Grids with flexible support
for GRID in (:PeriodicEquispacedGrid, :MidpointEquispacedGrid, :EquispacedGrid)
	@eval $GRID(n::Int, a, b) = $GRID(n, promote(a, b)...)
	@eval $GRID(n::Int, a::T, b::T) where {T} = $GRID{float(T)}(n, a, b)
	@eval $GRID(n::Int, d::AbstractInterval) =
		$GRID(n, infimum(d), supremum(d))
	@eval $GRID{T}(n::Int, d::AbstractInterval) where {T} =
		$GRID{T}(n, infimum(d), supremum(d))
    @eval similargrid(grid::$GRID, ::Type{T}, n::Int) where {T} =
        $GRID{T}(n, covering(grid))
    @eval rescale(grid::$GRID{T}, a::S, b::U) where {U,S,T} =
        $GRID{promote_type(float(S),float(U),T)}(length(grid), a, b)
    # @eval map_grid(grid::$GRID, map::AbstractAffineMap) =
    #     $GRID(length(grid), map_domain(map, covering(grid)))
end

map_grid(map::DomainSets.ScalarAffineMap, g::AbstractEquispacedRangeGrid) =
	rescale(g, endpoints(map_domain(map, covering(g)))...)

# extensible grids
# TODO: deprecate extend (in favour of resize)
extend(grid::AbstractEquispacedGrid, factor::Int) =
	resize(grid, extension_size(grid, length(grid), factor))

hasextension(::Union{PeriodicEquispacedGrid,FourierGrid,EquispacedGrid}) = true

extension_size(::PeriodicEquispacedGrid, n::Int, factor::Int) = factor*n
extension_size(::FourierGrid, n::Int, factor::Int) = factor*n
extension_size(::EquispacedGrid, n::Int, factor::Int) = factor*n-1


text:
module TestDualNumbers
    using Calculus2
    using Base.Test

    x = dual(2, 1)

    y = x^3
    @test_approx_eq(
        real(y),
        2^3
    )
    @test_approx_eq(
        epsilon(y),
        3*2^2
    )

    y = sin(x) + exp(x)
    @test_approx_eq(
        real(y),
        sin(2) + exp(2)
    )
    @test_approx_eq(
        epsilon(y),
        cos(2) + exp(2)
    )

    @test x > 1

    y = abs(-x)
    @test_approx_eq(
        real(y),
        2.0
    )
    @test_approx_eq(
        epsilon(y),
        1.0
    )

    @test isequal(1.0, dual(1.0))

    y = 1 / x
    @test_approx_eq(
        real(y),
        1 / 2
    )
    @test_approx_eq(
        epsilon(y),
        -1 / 2^2
    )

    Q = [
        1.0 0.1;
        0.1 1.0;
    ]
    x = dual([1.0, 2.0])
    x[1] = dual(1.0, 1.0)
    y = 1 / 2 * dot(x, Q * x)
    @test_approx_eq(
        real(y),
        2.7
    )
    @test_approx_eq(
        epsilon(y),
        1.2
    )
end


text:
# Command: Send Alert Text
"""
    alert_text!(session::Session, text::AbstractString)::Nothing

Sends keystrokes to a JavaScript prompt() dialog.

!!! note

    This command uses the JSON Wire Protocol instead of the current W3c WebDriver API.
"""
function alert_text!(session::Session, text::AbstractString)::Nothing
    @unpack addr, id = session
    response = HTTP.post(
        "$addr/session/$id/alert_text",
        [("Content-Type" => "application/json; charset=utf-8")],
        JSON3.write(Dict("text" => text)),
    )
    @assert response.status == 200
end


text:
using Documenter
using ForestBiometrics

makedocs(
        modules =[ForestBiometrics],
        sitename= "ForestBiometrics.jl",
        format = :html,
        authors="Casey Ghilardi",
        doctest=false,
        clean=true,
        pages = Any[
        "Home" => "index.md",
        "Functionality" =>Any["functionality/height_diameter.md",
        "functionality/density.md",
        "functionality/other_functions.md"],
        "Volume Equations" => Any["volume_equations/volume_eqs.md"],
        "Visualizations" =>Any["visualizations/visualizations.md"],
        "Included Datasets"=>Any["data/data.md"]
        ])

 deploydocs(
            deps = nothing,
            branch = "gh-pages",
            latest = "master",
            julia ="0.6",
            repo = "github.com/Crghilardi/ForestBiometrics.jl.git",
            target = "build",
            osname ="linux",
            make = nothing
   )


text:
export bloch_expand

function blochify(ii,jj,mm, naxis,nxbloch,nsector,naxis_ln,nsector_ln,N_points; axis = true)
    #nsector=naxis+nxbloch+nbody+nxsymmetry+nbody
    blochshift=nsector-naxis
    blochshift_ln=nsector_ln-naxis_ln# blochshift is necessary because point dofs are reduced by blochshift
    MM=ComplexF64[]
    MM_plus=ComplexF64[]
    MM_minus=ComplexF64[]
    MM_plus_axis=ComplexF64[]
    MM_minus_axis=ComplexF64[]
    MM_axis=ComplexF64[]

    II=UInt32[]
    JJ=UInt32[]
    II_plus=UInt32[]#TODO: consider preallocation
    JJ_plus=UInt32[]
    II_minus=UInt32[]
    JJ_minus=UInt32[]
    II_axis=UInt32[]
    JJ_axis=UInt32[]
    II_plus_axis=UInt32[]
    JJ_plus_axis=UInt32[]
    II_minus_axis=UInt32[]
    JJ_minus_axis=UInt32[]
    #println("$(length(ii)) $(length(jj)) $(length(mm))")
    for (i,j,m) in zip(ii,jj,mm)
        if i<=N_points #deal with point dof
            i_check = i>nsector
            # map Bloch_ref to Bloch_img
            if i_check
                i-=blochshift
            end
        else #deal with line dof
            i_check = i> nsector_ln
            if i_check
                i-=blochshift_ln
            end
        end
        if j<=N_points #deal with point dof
            j_check = j>nsector
            # map Bloch_ref to Bloch_img
            if j_check
                j-=blochshift
            end
        else #deal with line dof
            j_check = j> nsector_ln
            if j_check
                j-=blochshift_ln
            end
        end

        if axis && ( i<=naxis || j<=naxis || N_points<i<=naxis_ln || N_points<j<=naxis_ln) #&& !(i<=naxis && j<=naxis) && !(N_points>i<=naxis_ln && N_points>j<=naxis_ln) &&!(i<=naxis && N_points>j<=naxis_ln) && !(N_points>i<=naxis_ln && j<=naxis)
            axis_check=true
        else
            axis_check=false
        end

        #account for reduced number of point dofs
        if i>N_points
            i-=nxbloch
        end
        if j>N_points
            j-=nxbloch
        end

        #sort into operators matrices
        if (!i_check && !j_check) || (i_check && j_check)  #no manipulation of matrix entries
            if axis_check
                append!(II_axis,i)
                append!(JJ_axis,j)
                append!(MM_axis,m)
            else
                append!(II,i)
                append!(JJ,j)
                append!(MM,m)
            end

        elseif !i_check && j_check
            if axis_check
                append!(II_plus_axis,i)
                append!(JJ_plus_axis,j)
                append!(MM_plus_axis,m)
            else
                append!(II_plus,i)
                append!(JJ_plus,j)
                append!(MM_plus,m)
            end

        elseif i_check && !j_check
            if axis_check
                append!(II_minus_axis,i)
                append!(JJ_minus_axis,j)
                append!(MM_minus_axis,m)
            else
                append!(II_minus,i)
                append!(JJ_minus,j)
                append!(MM_minus,m)
            end
        else
            println("ERROR: in blochification")
            return nothing
        end
    end

    if naxis==0
        return (II,II_plus,II_minus), (JJ,JJ_plus,JJ_minus), (MM, MM_plus, MM_minus)
    else
        return (II, II_plus, II_minus, II_axis, II_plus_axis, II_minus_axis), (JJ, JJ_plus, JJ_minus, JJ_axis, JJ_plus_axis, JJ_minus_axis), (MM, MM_plus, MM_minus, MM_axis, MM_plus_axis, MM_minus_axis)
    end
end

"""
    v=bloch_expand(mesh::Mesh,sol::Solution,b=:b)

Expand solution vector `sol.v` on mesh `mesh`with Bloch wave number
`sol.params[b]` and return it as `v`.
"""
function bloch_expand(mesh::Mesh,sol::Solution,b=:b)
    naxis=mesh.dos.naxis
    nxsector=mesh.dos.nxsector
    DOS=mesh.dos.DOS
    v=zeros(ComplexF64,naxis+nxsector*DOS)
    v[1:naxis]=sol.v[1:naxis]
    B=sol.params[b]
    for s=0:DOS-1
        v[naxis+1+s*nxsector:naxis+(s+1)*nxsector]=sol.v[naxis+1:naxis+nxsector].*exp(+2.0im*pi/DOS*B*s)
    end
    return v
end

function bloch_expand(mesh::Mesh,vec::Array,b::Real=0)
    naxis=mesh.dos.naxis
    nxsector=mesh.dos.nxsector
    DOS=mesh.dos.DOS
    v=zeros(ComplexF64,naxis+nxsector*DOS)
    v[1:naxis]=vec[1:naxis]
    for s=0:DOS-1
        v[naxis+1+s*nxsector:naxis+(s+1)*nxsector]=vec[naxis+1:naxis+nxsector].*exp(+2.0im*pi/DOS*b*s)
    end
    return v
end