text:
using MultivariateStats

export iimg, iimg!
export interp3, interp3with01coords, resize, resizeminmax
export meshgrid, meshgrid3, centeredgrid, centeredmeshgrid, overlaygradient, toranges, tosize, tosize3
export imregionalmin, imregionalmax, monogen, bwlabel, bwlabel!, monoslic, border, sortcoords, bwdist
export blocks, blocks!!, cutbox
export rle, unrle
export stridedblocksub
export inpolygon, inpointcloud
export medfilt
export gausspos, gaussmask
export warp, rotate

#######################################
## iimg

iimg(a) = (r = float(a); iimg!(r); return r)

function iimg!{T}(a::AbstractArray{T,2})
    for x = 1:size(a,2)
        for y = 2:size(a,1)
            a[y,x] += a[y-1,x]
        end
    end

    for y = 1:size(a,1)
        for x = 2:size(a,2)
            a[y,x] += a[y,x-1]
        end
    end
end

if VERSION < v"0.5-"
    view_ = sub
else
    view_ = view
end

function iimg!{T}(a::AbstractArray{T,3})
    for z = 1:size(a,3)
        iimg!(view_(a,1:size(a,1),1:size(a,2),z))
    end

    for x = 1:size(a,2)
        for y = 1:size(a,1)
            for z = 2:size(a,3)
                a[y,x,z] += a[y,x,z-1]
            end
        end
    end
end

function interp3{T,NN,T2<:AbstractFloat}(a::Array{T,NN}, m_::T2, n_::T2, o_::T2)
  # if m_<1 || n_<1 || o_<1 || m_>size(a,1) || n_>size(a,2) || o_>size(a,3)
  #   # @show m_ n_ o_ size(a)
  #   error("interp3: index out of bounds: [$m_, $n_, $o_], size(a) == $(size(a)::Tuple{Int,Int,Int})")
  # end

  # if isnan(m_) || isnan(n_) || isnan(o_)
  #   # @show m_ n_ o_ size(a)
  #   error("interp3: index is NaN")
  # end

  one_ = one(T2)
  wm = one_-rem(m_, one_)::T2
  wn = one_-rem(n_, one_)::T2
  wo = one_-rem(o_, one_)::T2
  m = floor(Int,m_)
  M = ceil(Int,m_)
  n = floor(Int,n_)
  N = ceil(Int,n_)
  o = floor(Int,o_)
  O = ceil(Int,o_)

  #  @show wm wn wo m n o M N O
  # @show "got here" wm wn wo m M n N o O size(a)
  meanm!no = wm*a[m,n,o]+(1-wm)*a[M,n,o]
  meanm!No = wm*a[m,N,o]+(1-wm)*a[M,N,o]
  meanm!nO = wm*a[m,n,O]+(1-wm)*a[M,n,O]
  meanm!NO = wm*a[m,N,O]+(1-wm)*a[M,N,O]

  meanm!nOm!NO = wn*meanm!nO+(1-wn)*meanm!NO
  meanm!nom!No = wn*meanm!no+(1-wn)*meanm!No

  wo*meanm!nom!No + (1-wo)*meanm!nOm!NO
end

interp3with01coords(a,m,n,o) = interp3(a, 1+(size(a,1)-1)*m, 1+(size(a,2)-1)*n, 1+(size(a,3)-1)*o)

#####################################################
##   resize

resize(a, factor::Int; kargs...) = resize(a, asint(factor*siz(a)); kargs...)
resize(a, factor::Number; kargs...) = resize(a, asint(factor*siz(a)); kargs...)

# method: :nearest or :interp
function resize{T<:Real}(a::Array{T}, s::AbstractArray; method = :interp)
    ndims(a) == 3 && sizeo(a) == 3 && length(s) == 2 && return map(a, x->resize(x,s; method = method))

    if size(a) == tuple(s...)
        return Base.copy(a)
    end
    r = Array(method == :nearest ? eltype(a) : Float32, s...)
    mi = linspace(1., size(a,1), size(r,1))
    ni = linspace(1., size(a,2), size(r,2))
    oi = linspace(1., size(a,3), size(r,3))
    if method == :nearest
        mi = asint(mi)
        ni = asint(ni)
        oi = asint(oi)
        resize_kernel_nearest(r, a, mi, ni, oi)
    else
        resize_kernel_interp3(r, a, mi, ni, oi)
    end
    r
end
function resize_kernel_nearest(r, a, mi, ni, oi)
    for o = 1:size(r,3), n = 1:size(r,2), m = 1:size(r,1)
        r[m,n,o] = a[mi[m], ni[n], oi[o]]
    end
end
function resize_kernel_interp3(r, a, mi, ni, oi)
    for o = 1:size(r,3), n = 1:size(r,2), m = 1:size(r,1)
        r[m,n,o] = interp3(a, mi[m], ni[n], oi[o])
    end
end

#####################################################
##   resizeminmax

function resizeminmax(a, mins, maxs; kargs...)

    assert(ndims(a)==length(mins))
    assert(ndims(a)==length(maxs))

    s = [size(a)...]

    if Base.any(s .< mins)
        mi = minimum(s ./ mins)
        minewsiz = round(s / mi)
    else
        minewsiz = s
    end

    if Base.any(s .> maxs)
        ma = maximum(s ./ maxs)
        manewsiz = round(s / ma)
    else
        manewsiz = s
    end
    newsiz = round(Int,FunctionalData.clamp((minewsiz + manewsiz) / 2, mins, maxs))
    #@show size(a) newsiz
    resize(a, newsiz; kargs...)
end

function overlaygradient(img, sp, sp2 = sp*0, sp3 = sp*0)
    if size(img,1)!=size(sp,1) || size(img,2)!=size(sp,2)
        error("overlaygradient: sizes are not equal: size(img) == $(size(img)) and size(sp)==$(size(sp))")
    end
    sp = float(sp)
    sp2 = float(sp2)
    sp3 = float(sp3)
    if size(img,3)>1
        img = mean(img,3)
    end
    img = norm01(img)
    r = cat(3,img,img,img)

    function gradient(a)
        g = falses(size(a))
        for n = 1:size(a,2)-1, m = 1:size(a,1)
            g[m,n] |= a[m,n+1] != a[m,n]
        end
        for n = 1:size(a,2), m = 1:size(a,1)-1
            g[m,n] |= a[m+1,n] !=a[m,n]
        end
        g
    end

    gs = Base.map(gradient, Any[sp, sp2, sp3])

    for n = 1:size(img,2), m = 1:size(img,1)
    r[m,n,1] = gs[1][m,n] ? 1 : r[m,n,1]
    r[m,n,2] = gs[2][m,n] ? 1 : r[m,n,2]
    r[m,n,3] = gs[3][m,n] ? 1 : r[m,n,3]
    end
    r
end

toranges(a) = [a[1,i]:a[2,i] for i in 1:len(a)]
tosize(a) = tuple((a[2,:]-a[1,:].+1)...)
tosize3(a) = (size(a,2)<3 ? a = cat(2,a,[1 2]') : nothing; tosize(a))

function monogrid(rm,rn)
  m = [m for m in rm, n in rn]
  n = [n for m in rm, n in rn]
  return (m,n)
end
function monogrid(rm,rn,ro)
  m = [m for m in rm, n in rn, o in ro]
  n = [n for m in rm, n in rn, o in ro]
  o = [o for m in rm, n in rn, o in ro]
  return (m,n,o)
end
monogrid{T}(a::Array{T,2}) = monogrid(1:size(a,1), 1:size(a,2))
monogrid{T}(a::Array{T,3}) = monogrid(1:size(a,1), 1:size(a,2), 1:size(a,3))

meshgrid(a::Tuple) = meshgrid(a...)
meshgrid(a::Array) = meshgrid(size(a)...)
meshgrid(a::Int, v::Array) = meshgrid(repeat([a], ndims(v))...)
meshgrid(sm::Int, sn::Int) = meshgrid(1:sm, 1:sn)
meshgrid(sm::Int, sn::Int, so::Int) = meshgrid(1:sm, 1:sn, 1:so)
meshgrid(rm::AbstractArray, rn::AbstractArray) = [row([m for m in rm, n in rn]); row([n for m in rm, n in rn])]
meshgrid(rm::AbstractArray, rn::AbstractArray, ro::AbstractArray) = [row([m for m in rm, n in rn, o in ro]); row([n for m in rm, n in rn, o in ro]); row([o for m in rm, n in rn, o in ro])]

meshgrid3{T}(a::Array{T,2}) = meshgrid(reshape(a, size(a,1), size(a,2), 1))
meshgrid3{T}(a::Array{T,3}) = meshgrid(a)

centeredmeshgrid(a...) = let r = meshgrid(a...); ceil(Int, r .- mean(r,2)) end

# centeredgrid{T}(v::Array{T,2}) = grid((1:size(v,1))-size(v,1)/2, (1:size(v,2))-size(v,2)/2)
# centeredgrid{T}(v::Array{T,3}) = grid((1:size(v,1))-size(v,1)/2, (1:size(v,2))-size(v,2)/2, (1:size(v,3))-size(v,3)/2)
# centeredmeshgrid{T}(v::Array{T,2}) = meshgrid((1:size(v,1))-size(v,1)/2, (1:size(v,2))-size(v,2)/2)
# centeredmeshgrid{T}(v::Array{T,3}) = meshgrid((1:size(v,1))-size(v,1)/2, (1:size(v,2))-size(v,2)/2, (1:size(v,3))-size(v,3)/2)

function imregionalmin{T}(img::Array{T,2})
  r = falses(size(img))
  for n = 2:size(img,2)-1, m = 2:size(img,1)-1
    v = img[m,n]
    if v<img[m-1,n] && v<img[m+1,n] && v<img[m,n-1] && v<img[m,n+1]
      r[m,n] = true
    end
  end
  r
end
function imregionalmin{T}(img::Array{T,3})
  r = falses(size(img))
  for o = 2:size(img,3)-1, n = 2:size(img,2)-1, m = 2:size(img,1)-1
    v = img[m,n,o]
    if v<img[m-1,n,o] && v<img[m+1,n,o] && v<img[m,n-1,o] && v<img[m,n+1,o] && v<img[m,n,o-1] && v<img[m,n,o+1]
      r[m,n,o] = true
    end
  end
  r
end;
imregionalmax(img) = imregionmin(-float(img))

monogen(img, spacing) = monogen_(asfloat32(img), asfloat32(spacing*2))

function monogen_(img::Array{Float32,2}, wavelength::Float32)
  Base.FFTW.set_num_threads(nphysicalcores())

  myfft = fft
  myifft(x) = real(ifft(x))

  img = myfft(img)

  rows = size(img,1)
  cols = size(img,2)

  # Generate horizontal and vertical frequency grids that vary from
  (u2, u1) = monogrid([linspace(0.,0.5,floor(rows/2)); linspace(-0.5,0.,ceil(Int,rows/2))],
  [linspace(0.,0.5,floor(cols/2)); linspace(-0.5,.0,ceil(Int,cols/2))])

  radius = sqrt(u1.^2 + u2.^2)    # Matrix values contain frequency
  radius[1,1] = 1
  radius[end,1] = 1
  radius[1,end] = 1
  radius[end,end] = 1

  H1 = Complex64(0,1)*u1./radius   # The two monogenic filters in the frequency domain
  u1 = nothing
  H2 = Complex64(0,1)*u2./radius
  u2 = nothing

  logGabor = exp((-(log(radius*wavelength)).^2) / (2 * log(asfloat32(0.65)).^2))
  #radius = nothing
  logGabor[1,1] = 0                    #% undo the radius fudge.
  logGabor[end,1] = 0
  logGabor[1,end] = 0
  logGabor[end,end] = 0

  img = img.*logGabor
  logGabor=nothing
  h = myifft(img.*H1).^2
  H1=nothing
  psi = atan2(myifft(img),sqrt(h+myifft(img.*H2).^2))
  psi = atan2(myifft(img),sqrt(h+myifft(img.*H2).^2))
  h = nothing
  H2 = nothing
  psi = sign(psi).*exp(-abs(psi))
end

function monogen_(img::Array{Float32,3}, wavelength::Float32)
  if size(img,3)==1
    return monogen_(img[:,:,1], wavelength)
  end

  Base.FFTW.set_num_threads(nphysicalcores())

  myfft = fft
  myifft(x) = real(ifft(x))

  img = myfft(img)
  (rows,cols,channels) = size(img)
  (u2, u1, u3) = monogrid([linspace(0.,0.5,floor(rows/2)), linspace(-0.5,0.,ceil(Int,rows/2))],[linspace(0.,0.5,floor(cols/2)), linspace(-0.5,0.,ceil(Int,cols/2))],
  [linspace(0.,0.5,floor(channels/2)), linspace(-0.5,0.,ceil(Int,channels/2))]);

   radius = sqrt(u1.^2 + u2.^2 + u3.^2);    # Matrix values contain frequency
    # values as a radius from centre
    # (but quadrant shifted)

    # Get rid of the 0 radius value in the middle (at top left corner after
    # fftshifting) so that taking the log of the radius, or dividing by the
    # radius, will not cause trouble.
  radius[1,1,1] = 1;
    radius[1,1,end] = 1;
    radius[1,end,1] = 1;
    radius[1,end,end] = 1;
    radius[end,1,1] = 1;
    radius[end,1,end] = 1;
    radius[end,end,1] = 1;
    radius[end,end,end] = 1;
    H1 = u1./radius;   # The three monogenic filters in the frequency domain
  u1=nothing
    H2 = u2./radius;
  u2=nothing
    H3 = u3./radius;
  u3=nothing
    logGabor = exp((-(log(radius*wavelength)).^2) / (2 * asfloat32(log(0.65))^2));
  radius=nothing
  logGabor[1,1,1] = 0;                    # undo the radius fudge.
  logGabor[1,1,end] = 0;
  logGabor[1,end,1] = 0;
  logGabor[1,end,end] = 0;
  logGabor[end,1,1] = 0;
  logGabor[end,1,end] = 0;
  logGabor[end,end,1] = 0;
  logGabor[end,end,end] = 0;
  img = img.*logGabor;
  logGabor=nothing
  h = myifft(img.*(Complex64(0,1)*H1)).^2;
  H1=nothing
  h += myifft(img.*(Complex64(0,1)*H2)).^2;
  H2=nothing
  h += myifft(img.*(Complex64(0,1)*H3)).^2;
  H3=nothing

  psi = atan2(myifft(img), sqrt(h));
  img = nothing
    psi = sign(psi).*exp(-abs(psi));
end

bwlabel(img, args...) = (a = Base.copy(img); bwlabel!(a, args...); a)
bwlabel!(img) = bwlabel!(img, 1)
function bwlabel!(img, startlabel)
    seeds = Any[]
    labelind = round(Int, startlabel)-1
    for o = 1:size(img,3), n = 1:size(img,2), m = 1:size(img,1)
        if img[m,n,o]==0
            push!(seeds,[m,n,o])
            labelind += 1
            #@show seeds
            while length(seeds)>0
                sm,sn,so = pop!(seeds)
                if img[sm,sn,so] == 0
                    img[sm,sn,so] = labelind
                    for oo = -1:1, on = -1:1, om = -1:1
                        m2 = clamp(sm+om, 1, size(img,1))
                        n2 = clamp(sn+on, 1, size(img,2))
                        o2 = clamp(so+oo, 1, size(img,3))
                        if img[m2,n2,o2] == 0
                            push!(seeds,[m2,n2,o2])
                        end
                    end
                end
            end
        end
    end
end

@inbounds function monoslic_kernel(m,n,o,oi2s, ni2s, mi2s, lookup, newcenters, sv, labels)
  bestd = typemax(Float32)
    for oi2 = oi2s
      for ni2 = ni2s
          for mi2 = mi2s
              if lookup[mi2,ni2,oi2] > 0f0 #
                  diffm = m-newcenters[1,mi2,ni2,oi2]::Float32
                  diffn = n-newcenters[2,mi2,ni2,oi2]::Float32
                  diffo = o-newcenters[3,mi2,ni2,oi2]::Float32
                  d = (diffm*diffm + diffn*diffn + diffo*diffo)::Float32
                  if d < bestd
                      bestd = d
                      sv[m,n,o] = labels[mi2,ni2,oi2]
                  end
              end
          end
      end
  end
  nothing
end

@inbounds function processslice(o, B, clampedois, clampednis, clampedmis, newcenters, sv, labels, signs, notsigns)
    Bind = size(B,1)*size(B,2)*(o-1)::Int
    for n = 1:size(B,2)
        for m = 1:size(B,1)
            Bind += 1
            @inbounds b = B[Bind]
            lookup = b ? signs : notsigns

            monoslic_kernel(m,n,o, clampedois[o], clampednis[n], clampedmis[m], lookup, newcenters, sv, labels)
        end
    end
end

#monoslic(img, spacing; kargs...) = monoslic(asfloat32(img), asfloat32(spacing); kargs...)
#function monoslic(img::Array{Float32}, spacing::Float32; workers = Base.workers())
monoslic(img, spacing; kargs...) = monoslic(asfloat32(img), asfloat32(spacing))
function monoslic(img::Array{Float32}, spacing::Float32)
    println("computing monogen")
    @time M = monogen(img,spacing)

    if size(M,3) > 1
        oind = round(spacing/2):spacing:size(M,3)
    else
        oind = 1
    end

    mind = round(spacing/2):spacing:size(M,1)
    nind = round(spacing/2):spacing:size(M,2)
    centers  = zeros(Float32, 3, length(mind), length(nind), length(oind))
    labels   = zeros(Float32, length(mind), length(nind), length(oind))
    signs    = zeros(Float32, length(mind), length(nind), length(oind))
    bestdist = Inf32*ones(Float32, length(mind), length(nind), length(oind))
    bestind  = zeros(Float32, length(mind), length(nind), length(oind))

    label = 1;
    for o = 1:length(oind), n = 1:length(nind), m = 1:length(mind)
        centers[:,m,n,o] = [mind[m];nind[n];oind[o]]
        labels[m,n,o] = label
        label += 1
    end

    labels[:] = randperm(row(labels))

    cand = imregionalmin(abs(M))
    assert(size(cand)==size(M))
    cand = find(cand)

    newcenters = zeros(Float32, size(centers))
    (_, csm, csn, cso) = size(newcenters)

    myind2sub{T}(a::Array{T,2},i) = ((m,n)=ind2sub(size(a),i);(m,n,1))
    myind2sub(a,i) = ind2sub(size(a),i)
    for c = cand
        (m,n,o) = myind2sub(M,c)
        mi = round(Int, m/spacing+0.5)
        ni = round(Int, n/spacing+0.5)
        oi = round(Int, o/spacing+0.5)
        mi = min(mi, csm)
        ni = min(ni, csn)
        oi = min(oi, cso)
        d = (centers[1,mi,ni,oi]-m).^2 +
        (centers[2,mi,ni,oi]-n).^2 + (centers[3,mi,ni,oi]-o).^2
        if d<bestdist[mi,ni,oi]
            bestdist[mi,ni,oi] = d
            bestind[mi,ni,oi] = c
            newcenters[:,mi,ni,oi] = [m;n;o]
            signs[mi,ni,oi] = asfloat32(M[m,n,o] > 0)
        end
    end

    B = M.>0
    M = nothing
    notsigns = abs(signs.-1f0)

    maxd = typemax(eltype(newcenters))::Float32

    sv = zeros(Float32, size(B))
    oi2s = Array(Int,3)
    ni2s = Array(Int,3)
    mi2s = Array(Int,3)
    println("B loop")
    lookup = signs
    mis = asint(collect(1:size(B,1))./spacing.+0.5f0)
    nis = asint(collect(1:size(B,2))./spacing.+0.5f0)
    ois = asint(collect(1:size(B,3))./spacing.+0.5f0)
    clampedmis = [[clamp(mi-1,1,csm), clamp(mi,1,csm), clamp(mi+1,1,csm)] for mi in mis]
    clampednis = [[clamp(ni-1,1,csn), clamp(ni,1,csn), clamp(ni+1,1,csn)] for ni in nis]
    clampedois = [[clamp(oi-1,1,cso), clamp(oi,1,cso), clamp(oi+1,1,cso)] for oi in ois]
    Bind = 0

    @time for o = 1:size(B,3)
        processslice(o, B, clampedois, clampednis, clampedmis, newcenters, sv, labels, signs, notsigns)
    end
    bwlabel!(sv, maximum(sv)+1)
    uind = unique(asint(sv))
    labels = zeros(1,asint(maximum(uind)))
    labels[uind] = randperm(1:length(uind))
    asint(labels[asint(sv)])
end

border{T}(a::Array{T,2}) = @p border cat(3,a,a,a) | snd
function border{T}(a::Array{T,3})
    offsets = @p row [[m n o]'+1 for m in -1:1, n in -1:1, o in -1:1] | flatten | map subtoind a | (-) 1
    r = zeros(a)
    for ind = subtoind(2*ones(siz(a)),a):subtoind(siz(a)-1,a)
        a[ind] == 0 && continue
        isborder = false
        for i = 1:len(offsets)
            if a[ind+offsets[i]] == 0
                isborder = true
                break
            end
        end
        r[ind] = isborder
    end
    r
end

function sortcoords(coords)
    if sizem(coords)!=2
        error("FunctionalDataUtils.sortcoords: only works for 2D contours")
    end
    r = zeros(coords)
    setat!(r,1,fst(coords))
    nope = fill(typemax(eltype(coords)),size(fst(r)))
    setat!(coords, 1, nope)

    for i = 1:len(coords)-1
        ind = @p distance at(r,i) coords | indmin
        setat!(r,i+1, at(coords,ind))
        setat!(coords, ind, nope)
    end
    r
end

bwdist(a::Array) = @p bwdist findsub(a) meshgrid(a) | reshape siz(a)
function bwdist(a, pos)
    if size(a,1) == size(pos,1)
        @p map pos x->minimum(distance(x,a))
    else
        bwdist(findsub(a), pos)
    end
end

function blocks{T1,T2,N}(pos::Array{T1,2}, a::Array{T2,N}; blocksize = 32,
    scale = 1,
    grid = round(Int, scale .* centeredmeshgrid(repeat(blocksize, ndims(a))...)),
    borderstyle = :staysinside, precompute = false)

    if sizem(pos) == 2 && ndims(a) == 3 && sizeo(a) < 5
        r = @p map a (x->blocks(pos, x; blocksize = blocksize, scale = scale,
            borderstyle = borderstyle, precompute = precompute))
        return ndims(r) == 3 ? permutedims(r,[1,3,2]) : r
    end

    r = similar(a, blocksize^sizem(pos), len(pos))
    if sizem(pos) != ndims(a)
        error("sizem(pos)==$(sizem(pos)) does not match ndims(a)==$(ndims(a))")
    end

    inds = @p map grid+1 subtoind a | minus 1 | vec

    maxoffset = maximum(abs(grid))
    mi = ones(sizem(pos),1)*maxoffset+1
    ma = siz(a)-maxoffset
    blocks!!(r, copy(pos), a, inds, mi, ma, borderstyle)
    precompute ? (r, (inds, mi, ma, borderstyle)) : r
end

function blocks!!{T}(r::Matrix{T}, pos::Matrix{Int}, a::Matrix{T}, inds::Array{Int}, mi, ma, borderstyle)
    if borderstyle == :staysinside
        pos = @p clamp! pos mi ma
        blockssample_internal!(r,pos,a,inds)
    else
        error("unknown borderstyle '$borderstyle'")
    end
    r
end

function blockssample_internal!{T}(r, pos::Array{Int,2}, a::Array{T,2}, inds::Array{Int,1})
    for n = 1:len(pos)
        posind = sizem(a)*(pos[2,n]-1) + pos[1,n]
        for m = 1:length(inds)
            r[m,n] = a[inds[m]+posind]
        end
    end
end

cutbox{T}(ind, a::Array{T,2}) = a[ind[1,1]:ind[1,2],ind[2,1]:ind[2,2]]
cutbox{T}(ind, a::Array{T,3}) = a[ind[1,1]:ind[1,2],ind[2,1]:ind[2,2],:]

function rle(a)
    r = Dict()
    r[:size] = size(a)
    counters = zeros(Int, length(a))
    data = zeros(eltype(a), length(a))

    value = a[1]
    counter = 0
    dcounter = 1

    for i = 1:length(a)
        if a[i]!=value
            data[dcounter] = value
            counters[dcounter] = counter

            value = a[i]
            counter = 1
            dcounter += 1
        else
            counter += 1
        end
    end
    data[dcounter] = value
    counters[dcounter] = counter

    r[:counters] = counters[1:dcounter]
    r[:data] = data[1:dcounter]
    r
end

function unrle(a::Dict)
    data = a[:data]
    counters = a[:counters]
    r = zeros(eltype(data), a[:size])

    rind = 1
    for ind = 1:length(counters)
        value = data[ind]
        for i = 1:counters[ind]
            r[rind] = value
            rind += 1
        end
    end
    r
end

stridedblocksub{T}(a::Array{T,2}, blocksize::Number, stride = blocksize; kargs...) = stridedblocksub(a, blocksize*ones(Int,2,1), stride; kargs...)
function stridedblocksub{T}(a::Array{T,2}, blocksiz, stride = blocksiz; keepshape = false)
    if length(stride) == 1
        stride = asint(stride) * ones(Int,2,1)
    end
    sm, sn = size(a)
    rm = 0:blocksiz[1]-1
    rn = 0:blocksiz[2]-1
    r = [(m+rm,n+rn) for m in 1:stride[1]:sm-blocksiz[1]+1, n in 1:stride[2]:sn-blocksiz[2]+1]
    if length(r) == 0
        @show size(a) blocksiz stride
        error("length(r) == 0")
    end
    r :: Array{Tuple{UnitRange{Int64},UnitRange{Int64}},2}
    keepshape ? r : row(r)
end
stridedblocksub{T}(a::Array{T,3}, blocksize::Number, stride = blocksize; kargs...) = stridedblocksub(a, blocksize*ones(Int,3,1); kargs...)
function stridedblocksub{T}(a::Array{T,3}, blocksiz, stride = blocksiz; keepshape = false)
    if length(stride) == 1
        stride = asint(stride) * ones(Int,3,1)
    end
    sm, sn, so = size(a)
    rm = 0:blocksiz[1]-1
    rn = 0:blocksiz[2]-1
    ro = 0:blocksiz[3]-1
    r = [(m+rm,n+rn,o+ro) for m in 1:stride[1]:sm-blocksiz[1]+1, n in 1:stride[2]:sn-blocksiz[2]+1, o in 1:stride[3]:so-blocksiz[3]+1]
    assert(length(r)>0)
    r :: Array{Tuple{UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}},3}
    keepshape ? r : row(r)
end

inpolygon(point, polygon) = inpolygon(point[1], point[2], polygon)
function inpolygon(m::Int, n::Int, polygon)
    j = len(polygon)
    oddnodes = false
    M = polygon[1,:]
    N = polygon[2,:]

    for i in 1:len(polygon)
        if M[i] < m && M[j] >= m || M[j] < m && M[i] >= m
            if N[i] + (m-M[i]) / (M[j]-M[i]) * (N[j]-N[i]) < n
                oddnodes = !oddnodes
            end
        end
        j = i
    end

    oddnodes
end
poly2mask{T}(polygon::Array{T,2},m::Int,n::Int) = poly2mask(polygon, 1:m, 1:n)
poly2mask{T}(polygon::Array{T,2},img::Array{T,2}) = poly2mask(polygon, 1:size(img,1), 1:size(img,2))
poly2mask{T}(polygon::Array{T,2},M::AbstractArray,N::AbstractArray) = Float32[inpolygon(m,n,polygon) for m in M, n in N]

function inpointcloud(point, cloud, n = 8)
    nearest = at(cloud, indmin(distance(point,cloud)))
    #return distance(nearest,point)
    around = @p sortperm vec(distance(nearest,cloud)) | part cloud _ | take n
    pca = fit(PCA, around.+0.0001*randn(size(around)), pratio = 1., method = :svd)
    line = @p last projection(pca)
    #@show projection(pca)
    dists = distance(around, around)
    diameter = maximum(dists)
    relative = cloud.-nearest
    proj = line * (line' * relative)
    dist2line = sqrt(sum((proj-relative).^2,1))
    inpipe = @p part cloud find(dist2line.<=diameter/2)
    cloud2point = point - nearest
    a = sum(distance(nearest + cloud2point, inpipe))
    b = sum(distance(nearest - cloud2point, inpipe))
    #@show nearest cloud2point around dist2line line diameter inpipe a b size(inpipe)
    sum(distance(nearest + cloud2point, inpipe)) < sum(distance(nearest - cloud2point, inpipe))
end

function medfilt(a, b = 1)
    sm, sn = size(a)
    for n = b+1:sn-b-1, m = 1+b:sm-b-1
        a[m,n] = median(a[m:m,n-b:n+b])
    end
    a
end

function gausspos(a, n, std = 2)
    if length(a)<=3
        ndim = length(a)
        s = a
    else
        ndim = ndims(a)
        s = siz(a)
    end
    mi = ones(s)
    ma = s
    s2 = s/2
    S = s2/std
    @p randn ndim 2*n | times S | plus s2 | reject (x->x!=clamp(x,mi,ma)) | round Int _ | take n
end

gaussmask(a, stdstoedge = 1) = gaussmask(sizem(a), sizen(a), stdstoedge)
function gaussmask(sm::Number, sn::Number, stdstoedge = 1)
    m = mean([sm; sn])
    s = m/stdstoedge
    a = meshgrid(linspace(-m,m,sm), linspace(-m,m,sn))
    d2 = sum(a.*a,1)
    g = exp(-d2 / (2*s^2))
    g = norm01(g)
    r = reshape(g, sm, sn)
end

wf(v,p,alpha) = asfloat32(1./((1+distance(v,p)).^(2.*alpha)+0.0001))
function mlstransform(v::Array{Float32,2},source::Array{Float32,2},target::Array{Float32,2}, alpha = 1)
    q = source
    p = target
    w = wf(v,p,alpha)::Array{Float32,2}
    wnorm = normsum(w)
    pstar = sum(wnorm.*p,2)
    qstar = sum(wnorm.*q,2)
    phat = p.-pstar
    qhat = q.-qstar
    term1 = v-pstar
    term2 = inv(sum([phat[:,i]*w[i]*phat[:,i]' for i in 1:len(phat)]))
    term3 = sum([w[j]*phat[:,j]*qhat[:,j]' for j in 1:len(phat)])
    (term1'*term2*term3)'+qstar
end

function warp(img, from, to, targetsize = siz(img)[1:2]; alpha = 2)
    imgslice = img[:,:,1]
    targetsize = tuple(targetsize[1:2]...)
    coords = @p meshgrid targetsize | asfloat32
    from = asfloat32(from)
    to = asfloat32(to)
    from_m = mean(from,2)
    to_m = mean(to,2)
    from_ = from .- from_m
    to_ = to .- to_m

    coords = @p minus coords to_m | map (x->mlstransform(x,from_,to_,alpha)) | plus from_m
    ind = @p asint coords | clamp imgslice | map subtoind imgslice
    f(a) = @p getindex a ind | reshape targetsize
    ndims(img) == 2 ? f(img) : map(img,f)
end

function rotate{T}(a::Array{T,3}, args...; kargs...)
    @p map a x->rotate(x,args...; kargs...)
end

function rotate{T}(a::Array{T,2}, alpha; background = zero(T))
    r = background*ones(T,sizem(a)+2,sizen(a)+2)
    r[2:end-1,2:end-1] = a
    c = centeredmeshgrid(r)
    alpha = -alpha
    rotmat = [cos(alpha) -sin(alpha); sin(alpha) cos(alpha)]
    r = @p asint rotmat*c | plus div(siz(r),2)+1 | clamp r | map subtoind r | getindex r _ | reshape size(r)
    r[2:end-1,2:end-1]
end


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

using NodeNumbering

P = 9

Q = 12

adjacency = Dict(
                   1 => [3, 8, 9],
                   2 => [6, 7, 3, 8],
                   3 => [1, 2],
                   4 => [12, 13],
                   5 => [13, 16],
                   6 => [2, 7],
                   7 => [2, 6, 10],
                   8 => [1, 2, 10, 15],
                   9 => [1, 15],
                  10 => [8, 7],
                  11 => [12, 13, 14],
                  12 => [4, 11],
                  13 => [5, 16, 4, 11],
                  14 => [11, 16],
                  15 => [9, 8],
                  16 => [13, 5, 14]);

degrees = Dict(
                 1 => 3,
                 2 => 4,
                 3 => 2,
                 4 => 2,
                 5 => 2,
                 6 => 2,
                 7 => 3,
                 8 => 4,
                 9 => 2,
                10 => 2,
                11 => 3,
                12 => 2,
                13 => 4,
                14 => 2,
                15 => 2,
                16 => 3);

@testset "The Reverse Cuthill-McKee algorithm" begin
    expected = [16,5,14,13,11,4,12,6,7,2,10,3,8,1,15,9]
    result = NodeNumbering.RCM(adjacency, degrees, P, Q)
    @test result == expected
end

ext: jl


text:
"""Abstract type generalizes the embedding modules"""
abstract type Embedding <: KnetModule end

"""
`EmbeddingMul`: Basic linear embedding layer using matrix multiplication
for embedding.

# Constructors
    `EmbeddingMul(esize::Int, vsize::Int; winit)`
    esize is the embedding size,
    vsize is the vocabulary size

# Fields
    `l::Linear`: Linear embedding operator.

# Forward execution
    `forward(ctx, e::EmbeddingMul, seq)
    `@mc e(x)`
    `@run e(x)`

If `x` is 3d, `(size, batchsize, time)` is assumed.

"""
type EmbeddingMul <: Embedding
    l::Linear
end

EmbeddingMul(esize::Int, vsize::Int; winit=randn, dtype=Float32) =
    EmbeddingMul(Linear(esize, vsize;
                        winit=winit, dtype=Float32, bias=false))

function (e::EmbeddingMul)(ctx, x)
    if ndims(x) == 3
        dims = size(x, 2, 3)
        x = reshape(x, size(x,1), size(x, 2) * size(x, 3))
    else
        dims = size(x, 2)
    end
    emb = e.l(ctx, x)
    return reshape(emb, (size(emb,1), dims...))
end

"""
`EmbeddingLookup`: Basic linear embedding layer using index access.

# Constructors
    `EmbeddingLookup(esize::Int, vsize::Int; dtype=Float32, winit=randn)`
    esize is the embedding size,
    vsize is the vocabulary size

# Fields
    `w`: Embedding vector

# Forward execution
    `forward(ctx, e::EmbeddingLookup, indices::Array{Int})
    `@mc e(indices)`
    `@run e(indices)`

  output size will be [esize, size(indices)...]
"""
type EmbeddingLookup <: Embedding
    w::Param
end

EmbeddingLookup(emb::Int, vocab::Int; dtype=Float32, winit=randn) =
    EmbeddingLookup(Param(winit(dtype, emb, vocab)))

function (el::EmbeddingLookup)(ctx, indices)#forward(ctx, el::EmbeddingLookup, indices)
    w = val(ctx, el.w)
    emb = w[:, indices[:]]
    if ndims(indices) > 1
        emb = reshape(emb, (size(emb, 1), size(indices)...))
    end
    return emb
end


text:
using RegressionDiscontinuity
using Hypatia
using Roots, Distributions
using GLM, DataFrames
using Test

data = RDData(RegressionDiscontinuity.Lee08())

result = fit(
    ImbensWagerOptRD(
        B = 14.28,
        solver = Hypatia.Optimizer,
        num_buckets = 2000,
        variance = Homoskedastic(),
    ),
    data,
)

#test that the result is close to the R package result with the same settings
@test result.tau_est ≈ 0.05936146 atol = 0.001
@test result.se_est ≈ 0.01975009 atol = 0.0001
@test result.maxbias ≈ 0.01019907 atol = 0.0001

level = 0.95
rel_bias = result.maxbias / result.se_est
zz = fzero(
    z -> cdf(Normal(), rel_bias - z) + cdf(Normal(), -rel_bias - z) + level - 1,
    0,
    rel_bias - quantile(Normal(), (1 - level) / 2.1),
)

result_plusminus = result.tau_est - result.ci[1]
est_plusminus = zz * result.se_est

@test result_plusminus ≈ est_plusminus atol = 0.00001


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

using ACE, BenchmarkTools
using ACE.JacobiPolys: Jacobi
using JuLIP: evaluate

x = 2*rand() - 1
α, β = rand(), rand()
N = 30
P = zeros(N+1)
dP = zeros(N+1)
J = Jacobi(α, β, N)

evaluate!(P, nothing, J, x)
evaluate_d!(P, dP, nothing, J, x)
@info("Timing for evaluate!")
@btime evaluate!($P, nothing, $J, $x)
@btime evaluate!($P, nothing, $J, $x)
@info("Timing for evaluate_d!")
@btime evaluate_d!($P, $dP, nothing, $J, $x)
@btime evaluate_d!($P, $dP, nothing, $J, $x)

# # Julia Bug? => not a bug but a current limitation
# @info("Looking at that strange allocation?")
# using BenchmarkTools
# f(A, B) = A, B
# g(A, B) = A
# A = rand(5)
# B = rand(5)
# @btime f($A, $B);
# @btime f(1, 2);
# @btime g($A, $B);
#
# function runn(P, dP, J, x, N)
#    for n = 1:1000
#       evaluate_d!(P, dP, J, rand(), N)
#    end
#    return nothing
# end
#
# @info("Try again inside a function - no allocation!")
# @btime runn($P, $dP, $J, $x, 15)


text:
import Base: push!, eltype, close
export Signal, push!, value, preserve, unpreserve, close, rename!

using DataStructures

##### Signal #####

const debug_memory = false # Set this to true to debug gc of nodes

const nodeset = WeakKeyDict()
const io_lock = ReentrantLock()

const runner_task = Ref{Task}()

const run_remove_dead_nodes = Ref(false)

const CHANNEL_SIZE = 1024

#run in asynchronous mode by default
const async_mode = Ref(true)

const debug_mode = Ref(false)

const nodes = WeakRef[] #stores the nodes in order of creation (which is a topological order for execution of the nodes' actions)

const edges = Vector{Int}[] #parents to children, useful for plotting graphs

const node_count = Dict{String, Int}() #counts of different signals for naming

if !debug_memory
    type Signal{T}
        id::Int # also its index into `nodes`, and `edges`
        value::T
        parents::Tuple
        active::Bool
        actions::Vector{Function}
        preservers::Dict
        name::String
        function (::Type{Signal{T}}){T}(v, parents, pres, name)
            id = length(nodes) + 1
            n=new{T}(id, v, parents, false, Function[], pres, name)
            push!(nodes, WeakRef(n))
            push!(edges, Int[])
            foreach(p->push!(edges[p.id], id), parents)
            finalizer(n, schedule_node_cleanup)
            n
        end
    end
else
    type Signal{T}
        id::Int
        value::T
        parents::Tuple
        active::Bool
        actions::Vector
        preservers::Dict
        name::String
        bt
        function (::Type{Signal{T}}){T}(v, parents, actions, pres, name)
            id = length(nodes) + 1
            n=new{T}(id, v, parents, false, Function[], pres, name, backtrace())
            push!(nodes, WeakRef(n))
            push!(edges, Int[])
            foreach(p->push!(edges[p.id], id), parents)
            nodeset[n] = nothing
            finalizer(n, log_gc)
            n
        end
    end
end

Signal{T}(x::T, parents=(); name::String = auto_name!("input")) =
    Signal{T}(x, parents, Dict{Signal, Int}(), name)
Signal{T}(::Type{T}, x, parents=(); name::String = auto_name!("input")) =
    Signal{T}(x, parents, Dict{Signal, Int}(), name)
# A signal of types
Signal{T}(t::Type{T}; name::String = auto_name!("input")) = Signal(Type{T}, t, name=name)

function log_gc(n)
    @async begin
        lock(io_lock)
        print(STDERR, "Signal got gc'd. Creation backtrace:")
        Base.show_backtrace(STDERR, n.bt)
        println(STDOUT)
        unlock(io_lock)
    end
end

function auto_name!(name, parents...)
    parents_str = join(map(s->s.name, parents), ", ")
    isempty(parents_str) || (name *= "($parents_str)")
    count = get!(node_count, name, 0) + 1
    node_count[name] = count
    countstr = count > 1 ? "-$(node_count[name])" : ""
    "$name$countstr"
end

"""
`rename!(s::Signal, name::String)`

Change a Signal's name
"""
function rename!(s::Signal, name::String)
    s.name = name
end

# preserve/unpreserve nodes from gc
"""
    preserve(signal::Signal)

prevents `signal` from being garbage collected as long as any of its parents are around. Useful for when you want to do some side effects in a signal.
e.g. `preserve(map(println, x))` - this will continue to print updates to x, until x goes out of scope. `foreach` is a shorthand for `map` with `preserve`.
"""
function preserve(x::Signal)
    for p in x.parents
        p.preservers[x] = get(p.preservers, x, 0)+1
        preserve(p)
    end
    x
end

"""
    unpreserve(signal::Signal)

allow `signal` to be garbage collected. See also `preserve`.
"""
function unpreserve(x::Signal)
    for p in x.parents
        n = get(p.preservers, x, 0)-1
        if n <= 0
            delete!(p.preservers, x)
        else
            p.preservers[x] = n
        end
        unpreserve(p)
    end
    x
end

function Base.show(io::IO, n::Signal)
    active_str = isactive(n) ? "(active)" : ""
    write(io, "$(n.id): \"$(n.name)\" = $(n.value) $(eltype(n)) $active_str")
end

value(n) = n
value(n::Signal) = n.value
value(::Void) = false

eltype{T}(::Signal{T}) = T
eltype{T}(::Type{Signal{T}}) = T

##### Connections #####
const restart_queue = Ref(false)
function add_action!(f, node)
    push!(node.actions, f)
    if current_task() !== runner_task[]
        maybe_restart_queue()
    else
        restart_queue[] = true
    end
    f
end

"""
Removes `action` from `node.actions`
"""
function remove_action!(node, action)
    filter!(a -> a != action, node.actions)
    nothing
end

"""
Schedule a cleanup of dead nodes - called as a finalizer on each GC'd node
"""
function schedule_node_cleanup(n)
    run_remove_dead_nodes[] = true
end

"""
Remove GC'd nodes from `nodes`, is run before push! when scheduled.
Not thread-safe in the sense that if it is run while other code is iterating
through `nodes`, e.g. in run_push, iteration could skip nodes.
"""
function remove_dead_nodes!()
    allnodes = [] # keep a ref to all nodes to avoid GC whilst in this function
    filter!(nodes) do noderef
        node = noderef.value
        node != nothing && push!(allnodes, node)
        node != nothing
    end
    foreach((i_n)->((i,n) = i_n; n.id = i), enumerate(allnodes)) # renumber nodes
    reinit_edges!(allnodes)
    nothing
end

function reinit_edges!(allnodes)
    empty!(edges)
    foreach(n->push!(edges, []), allnodes)
    foreach(allnodes) do n
        foreach(p-> push!(edges[p.id], n.id), n.parents)
    end
    nothing
end

function close(n::Signal, warn_nonleaf=true)
    for p in n.parents
        delete!(p.preservers, n)
    end
    finalize(n) # stop timer, schedule_node_cleanup, etc.
end

set_value!(node::Signal, x) = (node.value = x)

##### Messaging #####

immutable Message
    node
    value
    onerror::Function
end

# Global channel for signal updates
const _messages = Channel{Nullable{Message}}(CHANNEL_SIZE)

run_async(async::Bool) = (async_mode[] = async)

"""
`push!(signal, value, onerror=Reactive.print_error)`

Queue an update to a signal. The update will be propagated when all currently
queued updates are done processing.

The third (optional) argument, `onerror`, is a callback triggered when
the update ends in an error. The callback receives 4 arguments,
`onerror(sig, val, node, capex)`, where `sig` and `val` are the Signal
and value that `push!` was called with, respectively, `node` is the
Signal whose action triggered the error, and `capex` is a
`CapturedException` with the fields `ex` which is the original
exception object, and `processed_bt` which is the backtrace of the
exception.

The default error callback will print the error and backtrace to STDERR.
"""
function Base.push!(n::Signal, x, onerror = print_error)
    if async_mode[]
        async_push!(n, x, onerror)
    else
        run_push(n, x, onerror)
    end
end

function async_push!(n, x, onerror=print_error)
    taken = Base.n_avail(_messages)
    if taken >= CHANNEL_SIZE
        warn("Message queue is full. Ordering may be incorrect.")
        @async put!(_messages, Message(n, x, onerror))
    else
        put!(_messages, Message(n, x, onerror))
    end
    nothing
end

function break_loop()
    put!(_messages, Nullable{Message}())
end

function stop()
    run_till_now() # process all remaining events
    break_loop()
    # it seems to take a long time until the runner_task is actually finished
    # which can be enough to run into maybe_restart_queue with !isdone(runner_task)
    # see #144
    wait(runner_task[])
end

"""
Processes `n` messages from the Reactive event queue.
"""
function run(n::Int=typemax(Int))
    for i=1:n
        message = take!(_messages)

        isnull(message) && break # break on null messages

        msgval = get(message)
        run_push(msgval.node, msgval.value, msgval.onerror)
        restart_queue[] && maybe_restart_queue()
    end
end

set_debug_mode(val = true) = (debug_mode[] = val)

runaction(f) = f()

activate!(node::Signal) = (node.active = true)
deactivate!(node::Signal) = (node.active = false)
isactive(node::Signal) = node.active

isactive(deadnode::Void) = false

activate!(noderef::WeakRef) = (noderef.value != nothing &&
                                (noderef.value.active = true))

deactivate!(noderef::WeakRef) = (noderef.value != nothing &&
                                (noderef.value.active = false))

isactive(noderef::WeakRef) = (noderef.value != nothing && noderef.value.active)

"""
A node's actions should be run if any of its parents are active, since that
generally means one (or more) of the parent nodes' values have changed. If the
node doesn't have actions, don't set it to active, since its value won't be
updated, meaning the update propagation can stop. N.b. The non-active when node
has no actions mechanism is relied on for correct behaviour by fpswhen, and
possibly other operators, i.e. it is not just an optimisation.
"""
function run_node(node::Signal)
    if any(isactive, node.parents) && length(node.actions) > 0
        activate!(node)
        foreach(runaction, node.actions)
    end
end

function run_push(pushnode::Signal, val, onerror, dont_remove_dead = false)
    node = pushnode # ensure node is set for error reporting - see onerror below
    try
        if run_remove_dead_nodes[] && !dont_remove_dead
            run_remove_dead_nodes[] = false
            remove_dead_nodes!()
        end
        set_value!(pushnode, val)
        activate!(pushnode)

        # run the actions for all downstream nodes
        for noderef in nodes[pushnode.id:end]
            noderef.value == nothing && continue
            node = noderef.value # node must be set here - see onerror call below
            run_node(node)
        end

        # reset active status to false for pushnode and all downstream nodes
        for noderef in nodes[pushnode.id:end]
            noderef.value == nothing && continue
            deactivate!(noderef.value)
        end
    catch err
        if isa(err, InterruptException)
            info("Reactive event loop was inturrupted.")
            rethrow()
        else
            bt = catch_backtrace()
            try
                onerror(pushnode, val, node, CapturedException(err, bt))
            catch err_onerror
                if isa(err_onerror, MethodError)
                    println(STDERR, "The syntax for `onerror` has changed, see ?push!")
                end
                rethrow()
            end
        end
    end
end

# Default error handler function
function print_error(pushnode, val, error_node, ex)
    lock(io_lock)
    io = STDERR
    println(io, "Failed to push!")
    print(io, "    ")
    show(io, val)
    println(io)
    println(io, "to node")
    print(io, "    ")
    show(io, pushnode)
    println(io)
    println(io)
    println(io, "error at node: $error_node")
    showerror(io, ex)
    println(io)
    unlock(io_lock)
end

function onerror_rethrow(pushnode, val, error_node, ex)
    rethrow()
end

# Run everything queued up till the instant of calling this function
run_till_now() = run(Base.n_avail(_messages))

# Works around world age problems (see issue #131)
function maybe_restart_queue()
    if !istaskdone(runner_task[])
        break_loop()
        if current_task() === runner_task[]
            # will happen if `add_action!` is called while processing a push!
            prev_runner = current_task()
            @async begin
                # new runner should wait for current runner to process the
                # break_loop (null) message
                wait(prev_runner)
                runner_task[] = current_task()
                run()
            end
        else
            wait(runner_task[])
            runner_task[] = @async run()
        end
    end
    restart_queue[] = false
end

function __init__()
    runner_task[] = @async run()
end


text:
Random.seed!(ω::Ω, T) = Random.seed!(rng(ω), T)
Random.seed!(ω::Ω, ::Type{T}) where T = Random.seed!(rng(ω), T)


text:
######## FANNO FLOW ########

function FLStarOverD(M::MachNumber,::Type{FannoFlow};gas::PerfectGas=DefaultPerfectGas)
  γ = SpecificHeatRatio(gas)
  return FLOverD((1-M^2)/(γ*M^2) + (γ+1)/(2γ)*log((γ+1)*M^2/(2+(γ-1)*M^2)))
end

function SubsonicMachNumber(fL_over_D::FLOverD,::Type{FannoFlow};gas::PerfectGas=DefaultPerfectGas)
  value(fL_over_D) >= 0.0 || error("fL*/D must be positive")
  MachNumber(_subsonic_mach_number_fanno(fL_over_D,gas))
end

function SupersonicMachNumber(fL_over_D::FLOverD,::Type{FannoFlow};gas::PerfectGas=DefaultPerfectGas)
  max_fLD = FLStarOverD(MachNumber(1e10),FannoFlow,gas=gas)
  value(fL_over_D) >= 0.0 || error("fL*/D must be positive")
  value(fL_over_D) > value(max_fLD) ? MachNumber(NaN) : MachNumber(_supersonic_mach_number_fanno(fL_over_D,gas))
end

function MachNumber(fL_over_D::FLOverD,::Type{FannoFlow};gas::PerfectGas=DefaultPerfectGas)
    SubsonicMachNumber(fL_over_D,FannoFlow,gas=gas), SupersonicMachNumber(fL_over_D,FannoFlow,gas=gas)
end

_subsonic_mach_number_fanno(fL_over_D::FLOverD,gas::PerfectGas) =
      find_zero(x -> FLStarOverD(MachNumber(x),FannoFlow,gas=gas)-fL_over_D,(0.001,1),order=16)

_supersonic_mach_number_fanno(fL_over_D::FLOverD,gas::PerfectGas) =
      find_zero(x -> FLStarOverD(MachNumber(x),FannoFlow,gas=gas)-fL_over_D,(1,1e10),order=16)

FLOverD(f::FrictionFactor,L::Length,D::Diameter) = FLOverD(f*L/D)
Length(fLoverD::FLOverD,D::Diameter,f::FrictionFactor) = Length(fLoverD*D/f)

function POverPStar(M::MachNumber,::Type{FannoFlow};gas::PerfectGas=DefaultPerfectGas)
  γ = SpecificHeatRatio(gas)
  return PressureRatio(1/M*sqrt((1+γ)/(2+(γ-1)*M^2)))
end

function ρOverρStar(M::MachNumber,::Type{FannoFlow};gas::PerfectGas=DefaultPerfectGas)
  γ = SpecificHeatRatio(gas)
  return DensityRatio(1/M*sqrt((2+(γ-1)*M^2)/(γ+1)))
end

function TOverTStar(M::MachNumber,::Type{FannoFlow};gas::PerfectGas=DefaultPerfectGas)
  γ = SpecificHeatRatio(gas)
  return TemperatureRatio((γ+1)/(2+(γ-1)*M^2))
end

function MachNumber(T_over_Tstar::TemperatureRatio,::Type{FannoFlow};gas::PerfectGas=DefaultPerfectGas)
    M = find_zero(x -> TOverTStar(MachNumber(x),FannoFlow,gas=gas)-T_over_Tstar,(1e-10,1e10),order=16)
    return MachNumber(M)

end

function P0OverP0Star(M::MachNumber,::Type{FannoFlow};gas::PerfectGas=DefaultPerfectGas)
  γ = SpecificHeatRatio(gas)
  return StagnationPressureRatio(1/M*((2+(γ-1)*M^2)/(γ+1))^((γ+1)/(2(γ-1))))
end

function MachNumber(p0_over_p0star::StagnationPressureRatio,::Type{FannoFlow};gas::PerfectGas=DefaultPerfectGas)
    Msub = find_zero(x -> P0OverP0Star(MachNumber(x),FannoFlow,gas=gas)-p0_over_p0star,(0.001,1),order=16)
    max_p0_over_p0star = P0OverP0Star(MachNumber(1e10),FannoFlow,gas=gas)

    if value(p0_over_p0star) > value(max_p0_over_p0star)
        return MachNumber(Msub)
    else
        Msup = find_zero(x -> P0OverP0Star(MachNumber(x),FannoFlow,gas=gas)-p0_over_p0star,(1,1e10),order=16)
        return MachNumber(Msub), MachNumber(Msup)
    end
end


text:
@testset "basic store-load" begin
    p = Properties()
    setprop(p, "key1", "value1")
    setprop(p, "key2", 157)
    setprop(p, "key3", 157.876341)
    setprop(p, "key4", raw"C:\Users\igor\elefant\.jl")

    fname = "p.jlprop"
    store(p, fname)
    @test isfile(fname)
    pl = load(fname)
    rm(fname)
    @test length(pairs(p)) == length(pairs(pl))
    pkeys = keys(p)
    for pk in pkeys
        @test getprop(p, pk) == getprop(pl, pk)
    end
end

@testset "multiline comment" begin
    p = Properties()
    setprop(p, "key1", "value1")
    setprop(p, "key2", 157)
    setprop(p, "key3", 157.876341)
    setprop(p, "key4", raw"C:\Users\igor\elefant\.jl")

    fname = "p2.jlprop"
    store(p, fname, "this\nis a multiline\ncomment \n and another one")
    @test isfile(fname)

    # test if the comments are written correctly
    lines = readlines(fname)
    commlines = ["#this", "#is a multiline", "#comment ", "# and another one"]
    for i in 1:4
        @test lines[i] == commlines[i]
    end

    pl = load(fname)
    rm(fname)
    @test length(pairs(p)) == length(pairs(pl))
    pkeys = keys(p)
    for pk in pkeys
        @test getprop(p, pk) == getprop(pl, pk)
    end
end


text:
""" Dose response curve (ODE model), Figure 1. """

using DataFrames;
using gcSolver;
using Gadfly;
using Statistics;
using StatsFuns;
gdf = Gadfly;

# Plot of dose response curves
function doseResPlot(ligandName, cellType, date, unkVec, biv = false)
    responseDF = gcSolver.importData(false)
    time = [0.5, 1, 2, 4] .* 60
    doseVec = unique(responseDF, "Dose")
    doseVec = doseVec[!, :Dose]

    predictDF = DataFrame(Dose = Float64[], time = Float64[], pSTAT = Float64[])

    filtFrame = filter(row -> row["Ligand"] .== ligandName, responseDF)
    filter!(row -> row["Cell"] .== cellType, filtFrame)
    filter!(row -> row["Date"] .== date, filtFrame)
    filter!(row -> row["Bivalent"] .== biv, filtFrame)

    realDataDF = filtFrame[!, [:Dose, :Time, :Mean]]
    realDataDF = groupby(realDataDF, [:Time, :Dose])
    realDataDF = combine(realDataDF, :Mean => mean)

    for (i, dose) in enumerate(doseVec)

        # Place ligand dose
        if ligandName == "IL15"
            doseLevel = [0, dose, 0]
        else
            doseLevel = [dose, 0, 0]
        end

        #Gives back 36 parameter long
        idxx = findfirst(responseDF.Cell .== cellType)
        iterParams =
            gcSolver.fitParams(doseLevel, unkVec, 10.0 .^ Vector{Float64}(responseDF[idxx, [:IL2Ra, :IL2Rb, :gc, :IL15Ra, :IL7Ra]]), cellType)
        if ligandName != "IL2" && ligandName != "IL15"
            iterParams = gcSolver.mutAffAdjust(iterParams, responseDF[findfirst(responseDF.Ligand .== ligandName), [:IL2RaKD, :IL2RBGKD]])
        end
        #gives you pstat results
        pstatResults = gcSolver.runCkine(time, iterParams, pSTAT5 = true)
        for indx = 1:length(time)
            #use dataframe and push row into it - enter data into data frame
            push!(predictDF, (dose, time[indx] / 60, pstatResults[indx]))
        end
    end

    DateFrame = gcSolver.importConvFrame()
    predictDF.pSTAT .*= filter(row -> row.Date ∈ [date], DateFrame).Conv

    pl1 = gdf.plot(
        layer(realDataDF, x = :Dose, y = :Mean_mean, color = :Time, Geom.point),
        layer(predictDF, x = :Dose, y = :pSTAT, color = :time, Geom.line),
        Scale.x_log10,
        Guide.title(string(cellType, " Response to ", ligandName)),
        Guide.xlabel("Dose"),
        Guide.ylabel("pSTAT Level"),
        Scale.color_discrete(),
        Guide.colorkey(title = "Time (hr)", labels = ["4", "2", "1", "0.5"]),
    )

    return pl1
end

"""This function compares the fit of Farhat and our Jak/STAT pathway model"""
function fitCompare(fitVec)
    numParams = 18
    farhatParams = gcSolver.getFarhatVec()[1:numParams]
    fitVecComp = fitVec[1:numParams]
    paramNames = gcSolver.getParamNames()[1:numParams]
    compDF = DataFrame(
        Parameter = repeat(paramNames, 2),
        Value = append!(farhatParams, fitVecComp),
        Source = append!(fill("Farhat", numParams), fill("Jak/STAT", numParams)),
    )
    pl1 = gdf.plot(compDF, y = :Value, x = :Parameter, color = :Source, Scale.y_log10)
    return pl1
end

"""Use this if you want to change the parameters here and not input any in the command line"""
function figureJ1()
    fitVec = gcSolver.importFit()
    fitVec = convert(Vector{Float64}, fitVec[!, :Fit])
    fitVec = softplus.(fitVec)

    pParam = fitCompare(fitVec)
    p1 = doseResPlot("IL2", "Treg", "4/18/2019", fitVec, 0)
    p2 = doseResPlot("IL2", "Thelper", "4/18/2019", fitVec, 0)
    p3 = doseResPlot("IL2", "NK", "4/18/2019", fitVec, 0)
    p4 = doseResPlot("IL2", "CD8", "4/18/2019", fitVec, 0)
    p5 = doseResPlot("N88D C-term", "Treg", "12/5/2019", fitVec, 0)
    p6 = doseResPlot("N88D C-term", "Thelper", "12/5/2019", fitVec, 0)
    p7 = doseResPlot("N88D C-term", "NK", "12/5/2019", fitVec, 0)
    p8 = doseResPlot("N88D C-term", "CD8", "12/5/2019", fitVec, 0)
    p9 = doseResPlot("WT N-term", "Treg", "11/8/2019", fitVec, 0)
    p10 = doseResPlot("WT N-term", "Thelper", "11/8/2019", fitVec, 0)
    p11 = doseResPlot("WT N-term", "NK", "11/8/2019", fitVec, 0)
    p12 = doseResPlot("WT N-term", "CD8", "11/8/2019", fitVec, 0)
    p13 = doseResPlot("V91K C-term", "Treg", "11/27/2019", fitVec, 0)
    p14 = doseResPlot("V91K C-term", "Thelper", "11/27/2019", fitVec, 0)
    p15 = doseResPlot("V91K C-term", "NK", "11/27/2019", fitVec, 0)
    p16 = doseResPlot("V91K C-term", "CD8", "11/27/2019", fitVec, 0)
    p17 = doseResPlot("R38Q N-term", "Treg", "12/5/2019", fitVec, 0)
    p18 = doseResPlot("R38Q N-term", "Thelper", "12/5/2019", fitVec, 0)
    p19 = doseResPlot("R38Q N-term", "NK", "12/5/2019", fitVec, 0)
    p20 = doseResPlot("R38Q N-term", "CD8", "12/5/2019", fitVec, 0)
    p21 = doseResPlot("F42Q N-Term", "Treg", "12/5/2019", fitVec, 0)
    p22 = doseResPlot("F42Q N-Term", "Thelper", "12/5/2019", fitVec, 0)
    p23 = doseResPlot("F42Q N-Term", "NK", "12/5/2019", fitVec, 0)
    p24 = doseResPlot("F42Q N-Term", "CD8", "12/5/2019", fitVec, 0)
    p25 = doseResPlot("IL15", "Treg", "4/18/2019", fitVec, 0)
    p26 = doseResPlot("IL15", "Thelper", "4/18/2019", fitVec, 0)
    p27 = doseResPlot("IL15", "NK", "4/18/2019", fitVec, 0)
    p28 = doseResPlot("IL15", "CD8", "4/18/2019", fitVec, 0)

    #draw(SVG("figureJ2.svg", 1000px, 800px), p1)
    draw(
        SVG("figureJ1.svg", 4000px, 2400px),
        gridstack(
            [
                pParam pParam pParam pParam
                p1 p2 p3 p4
                p5 p6 p7 p8
                p9 p10 p11 p12
                p13 p14 p15 p16
                p17 p18 p19 p20
                p21 p22 p23 p24
                p25 p26 p27 p28
            ],
        ),
    )
end


text:
function add_event_handler(obj::AbstractFieldIntegrator, arg0::FieldODEEventHandler, arg1::jdouble, arg2::jdouble, arg3::jint)
    return jcall(obj, "addEventHandler", void, (FieldODEEventHandler, jdouble, jdouble, jint), arg0, arg1, arg2, arg3)
end

function add_event_handler(obj::AbstractFieldIntegrator, arg0::FieldODEEventHandler, arg1::jdouble, arg2::jdouble, arg3::jint, arg4::BracketedRealFieldUnivariateSolver)
    return jcall(obj, "addEventHandler", void, (FieldODEEventHandler, jdouble, jdouble, jint, BracketedRealFieldUnivariateSolver), arg0, arg1, arg2, arg3, arg4)
end

function add_step_handler(obj::AbstractFieldIntegrator, arg0::FieldODEStepHandler)
    return jcall(obj, "addStepHandler", void, (FieldODEStepHandler,), arg0)
end

function clear_event_handlers(obj::AbstractFieldIntegrator)
    return jcall(obj, "clearEventHandlers", void, ())
end

function clear_step_handlers(obj::AbstractFieldIntegrator)
    return jcall(obj, "clearStepHandlers", void, ())
end

function compute_derivatives(obj::AbstractFieldIntegrator, arg0::RealFieldElement, arg1::Vector{RealFieldElement})
    return jcall(obj, "computeDerivatives", Vector{RealFieldElement}, (RealFieldElement, Vector{RealFieldElement}), arg0, arg1)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_current_signed_stepsize(obj::AbstractFieldIntegrator)
    return jcall(obj, "getCurrentSignedStepsize", RealFieldElement, ())
end

function get_evaluations(obj::AbstractFieldIntegrator)
    return jcall(obj, "getEvaluations", jint, ())
end

function get_event_handlers(obj::AbstractFieldIntegrator)
    return jcall(obj, "getEventHandlers", Collection, ())
end

function get_field(obj::AbstractFieldIntegrator)
    return jcall(obj, "getField", Field, ())
end

function get_max_evaluations(obj::AbstractFieldIntegrator)
    return jcall(obj, "getMaxEvaluations", jint, ())
end

function get_name(obj::AbstractFieldIntegrator)
    return jcall(obj, "getName", JString, ())
end

function get_step_handlers(obj::AbstractFieldIntegrator)
    return jcall(obj, "getStepHandlers", Collection, ())
end

function get_step_start(obj::AbstractFieldIntegrator)
    return jcall(obj, "getStepStart", FieldODEStateAndDerivative, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function integrate(obj::FieldODEIntegrator, arg0::FieldExpandableODE, arg1::FieldODEState, arg2::RealFieldElement)
    return jcall(obj, "integrate", FieldODEStateAndDerivative, (FieldExpandableODE, FieldODEState, RealFieldElement), arg0, arg1, arg2)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function set_max_evaluations(obj::AbstractFieldIntegrator, arg0::jint)
    return jcall(obj, "setMaxEvaluations", void, (jint,), arg0)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
"""
Functionality for PyTorch optimisers.
"""

struct PyTorchOptimiser
    optim
end

PyTorchOptimiser(opt::ADAM, ps::PyTorchParams) =
    PyTorchOptimiser(torch.optim.Adam(ps.params, lr=opt.eta, betas=opt.beta))

convertoptimiser(opt::Flux.Optimise.AbstractOptimiser, params) = PyTorchOptimiser(opt, params)
convertoptimiser(opt::PyTorchOptimiser, _) = opt


text:
const REQUIRED_KEYS = Set(["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"])
const PASSPORTS = split(read("day04.txt", String), "\n\n")

const PATTERN1 = r"([a-z]{3}):"

function validate_part1(passport::AbstractString; required_keys=REQUIRED_KEYS, pattern=PATTERN1)
    found_keys = Set(m.captures[1] for m in eachmatch(pattern, passport))
    return found_keys >= required_keys
end

println("Part 1: $(sum(validate_part1(passport) for passport in PASSPORTS))")

const PATTERN2 = r"([a-z]{3}):([a-z0-9#]+)"

function validate_number_between(value::AbstractString, min::Int, max::Int)
    val = tryparse(Float64, value)
    if val === nothing
        return false
    end
    return min <= val <= max
end

function validate_height(value::AbstractString)
    if endswith(value, "cm")
        return validate_number_between(value[begin:end - 2], 150, 193)
    elseif endswith(value, "in")
        return validate_number_between(value[begin:end - 2], 59, 76)
    else
        return false
    end
end

function validate_hair_color(value::AbstractString)
    if startswith(value, '#')
        return occursin(r"[a-f0-9]{6}", value[begin + 1:end])
    else
        return false
    end
end

function validate_eye_color(value::AbstractString)
    return value in ["amb", "blu", "brn", "gry", "grn", "hzl", "oth"]
end

function validate_pid(value::AbstractString)
    return occursin(r"^[0-9]{9}$", value)
end

function validate_passport(passport::AbstractString; required_keys=REQUIRED_KEYS, pattern=PATTERN2)
    data = Dict(m.captures for m in eachmatch(pattern, passport))
    if !validate_part1(passport)
        return false
    end
    return (
        validate_number_between(data["byr"], 1920, 2002) &&
        validate_number_between(data["iyr"], 2010, 2020) &&
        validate_number_between(data["eyr"], 2020, 2030) &&
        validate_height(data["hgt"]) &&
        validate_hair_color(data["hcl"]) &&
        validate_eye_color(data["ecl"]) &&
        validate_pid(data["pid"])
    )
end

println("Part 2: $(sum(validate_passport(passport) for passport in PASSPORTS))")


text:
module ProperOrthogonalDecomposition
using Statistics
using LinearAlgebra

export  PODBasis,
        POD,
        PODsvd,
        PODsvd!,
        PODeigen,
        PODeigen!,
        modeConvergence,
        modeConvergence!

include("types.jl")
include("convergence.jl")
include("PODsvd.jl")
include("PODeigen.jl")

"""
    POD(X)

Computes the Proper Orthogonal Decomposition of `X`. Returns a PODBasis.

This method is just a thin wrapper around `PODsvd(X)`.
"""
function POD(X; subtractmean::Bool = false)
    PODsvd(X, subtractmean = subtractmean)
end

"""
    POD(X, W)

Computes the Proper Orthogonal Decomposition of `X` with weights `W`. Returns a PODBasis.

This method is just a thin wrapper around `PODsvd(X)`.
"""
function POD(X, W; subtractmean::Bool = false)
    PODsvd(X, W, subtractmean = subtractmean)
end

end # module


text:
# Autogenerated wrapper script for Qt5Base_jll for aarch64-apple-darwin
export libqt5concurrent, libqt5core, libqt5dbus, libqt5gui, libqt5network, libqt5opengl, libqt5printsupport, libqt5sql, libqt5test, libqt5widgets, libqt5xml

using Xorg_libXext_jll
using Xorg_libxcb_jll
using Xorg_xcb_util_wm_jll
using Xorg_xcb_util_image_jll
using Xorg_xcb_util_keysyms_jll
using Xorg_xcb_util_renderutil_jll
using xkbcommon_jll
using Libglvnd_jll
using Fontconfig_jll
using Glib_jll
using Zlib_jll
using CompilerSupportLibraries_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("Qt5Base")
JLLWrappers.@declare_library_product(libqt5concurrent, "@rpath/QtConcurrent")
JLLWrappers.@declare_library_product(libqt5core, "@rpath/QtCore")
JLLWrappers.@declare_library_product(libqt5dbus, "@rpath/QtDBus")
JLLWrappers.@declare_library_product(libqt5gui, "@rpath/QtGui")
JLLWrappers.@declare_library_product(libqt5network, "@rpath/QtNetwork")
JLLWrappers.@declare_library_product(libqt5opengl, "@rpath/QtOpenGL")
JLLWrappers.@declare_library_product(libqt5printsupport, "@rpath/QtPrintSupport")
JLLWrappers.@declare_library_product(libqt5sql, "@rpath/QtSql")
JLLWrappers.@declare_library_product(libqt5test, "@rpath/QtTest")
JLLWrappers.@declare_library_product(libqt5widgets, "@rpath/QtWidgets")
JLLWrappers.@declare_library_product(libqt5xml, "@rpath/QtXml")
function __init__()
    JLLWrappers.@generate_init_header(Xorg_libXext_jll, Xorg_libxcb_jll, Xorg_xcb_util_wm_jll, Xorg_xcb_util_image_jll, Xorg_xcb_util_keysyms_jll, Xorg_xcb_util_renderutil_jll, xkbcommon_jll, Libglvnd_jll, Fontconfig_jll, Glib_jll, Zlib_jll, CompilerSupportLibraries_jll, OpenSSL_jll)
    JLLWrappers.@init_library_product(
        libqt5concurrent,
        "lib/QtConcurrent.framework/Versions/5/QtConcurrent",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5core,
        "lib/QtCore.framework/Versions/5/QtCore",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5dbus,
        "lib/QtDBus.framework/Versions/5/QtDBus",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5gui,
        "lib/QtGui.framework/Versions/5/QtGui",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5network,
        "lib/QtNetwork.framework/Versions/5/QtNetwork",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5opengl,
        "lib/QtOpenGL.framework/Versions/5/QtOpenGL",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5printsupport,
        "lib/QtPrintSupport.framework/Versions/5/QtPrintSupport",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5sql,
        "lib/QtSql.framework/Versions/5/QtSql",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5test,
        "lib/QtTest.framework/Versions/5/QtTest",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5widgets,
        "lib/QtWidgets.framework/Versions/5/QtWidgets",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5xml,
        "lib/QtXml.framework/Versions/5/QtXml",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using NumCME
using BenchmarkTools
using Catalyst
using Sundials: CVODE_BDF
using StaticArrays: @MVector

fspalgorithm = AdaptiveFspSparse(
    ode_method=CVODE_BDF(linear_solver=:GMRES),
    space_adapter=RStepAdapter(5, 10, true)
)

# Bursting gene model definition using basic Julia
𝕊 = [[-1, 1, 0] [1, -1, 0] [0, 0, 1] [0, 0, -1]]
x₀ = [1, 0, 0]

a1 = propensity() do x, p
    p[1] * x[1]
end
a2 = propensity() do x, p
    p[2] * x[2]
end
a3 = propensity() do x, p
    p[3] * x[2]
end
a4 = propensity() do x, p
    p[4] * x[3]
end

k₀₁ = 0.05
k₁₀ = 0.1
λ = 5.0
γ = 0.5
θ = [k₀₁, k₁₀, λ, γ]

model = CmeModel(𝕊, [a1, a2, a3, a4], θ)

p0 = FspVectorSparse([@MVector [1, 0, 0]], [1.0])
tspan = (0.0, 300.0)

fspsol1 = solve(model, p0, tspan, fspalgorithm)

# Bursting model definition using Catalyst
@parameters k₀₁ k₁₀ λ γ
bursting_rn = @reaction_network begin
    k₀₁, G0 --> G1
    k₁₀, G1 --> G0
    λ, G1 --> G1 + mRNA
    γ, mRNA --> ∅
end k₀₁ k₁₀ λ γ

parameter_values = [k₀₁ => 0.05, k₁₀ => 0.1, λ => 5.0, γ => 0.5]
model_from_catalyst = CmeModel(bursting_rn, parameter_values)

fspsol2 = solve(model_from_catalyst, p0, tspan, fspalgorithm)

# Check that the two ways to code the model lead to the same numerical outputs
@assert length(fspsol1) == length(fspsol2)
for i ∈ 1:length(fspsol1)
    @assert get_states(fspsol1[i].p) == get_states(fspsol2[i].p)
    @assert get_values(fspsol1[i].p) == get_values(fspsol2[i].p)
end

# Which method is faster?
@btime solve(model, p0, tspan, fspalgorithm);
@btime solve(model_from_catalyst, p0, tspan, fspalgorithm);


text:
# evaluate right bound of scan_func
function get_right_endpoint(
    theta_init::Vector{Float64}, # initial point of parameters
    scan_func::Function, # h(theta) function for predictions or parameters
    loss_func::Function, # lambda(theta) - labmbda_min - delta_lambda
    method::Val{:CICO_ONE_PASS}; # function works only for method ONE_PASS;

    theta_bounds::Vector{Tuple{Float64,Float64}} = fill(
        (-Inf, Inf), length(theta_init)
        ),
    scan_bound::Float64 = 9.0,
    scan_tol::Float64 = 1e-3,
    scan_rtol::Float64 = 0.,
    loss_tol::Float64 = 0., # we have no idea how to implement loss_tol
    # better results for LN_NELDERMEAD, :LD_MMA, :LD_SLSQP, :LD_CCSAQ
    # worse results for :LN_COBYLA, :LN_PRAXIS
    # errors for :LN_BOBYQA, :LN_SBPLX, :LN_NEWUOA
    local_alg::Symbol = :LN_NELDERMEAD,
    # options for local fitter :max_iter
    max_iter::Int = 10^5,
    scan_grad::Union{Function, Symbol} = :EMPTY,
    loss_grad::Union{Function, Symbol} = :EMPTY,
    #kwargs...
)
    # dim of the theta vector
    n_theta = length(theta_init)

    # checking arguments
    # when using :LN_NELDERMEAD initial parameters should not be zero
    if local_alg == :LN_NELDERMEAD
        zeroParameter = [ isapprox(theta_init[i], 0., atol=1e-2) for i in 1:n_theta]
        if any(zeroParameter)
            @warn "Close-to-zero parameters found when using :LN_NELDERMEAD."
            show(findall(zeroParameter))
        end
    end

    # checking scan_grad, loss_grad
    is_gradient = occursin(r"^LD_", String(local_alg))
    if scan_grad == :EMPTY && is_gradient
        throw(ArgumentError("`scan_grad` must be set for gradient local fitter `$(local_alg)`"))
    end
    if loss_grad == :EMPTY && is_gradient
        throw(ArgumentError("`loss_grad` must be set for gradient local fitter `$(local_alg)`"))
    end

    # optimizer
    local_opt = Opt(local_alg, n_theta)
    ftol_abs!(local_opt, scan_tol)
    ftol_rel!(local_opt, scan_rtol)
    # local_opt.initial_step = fill(1., n_theta) # changes fitting but still unstable

    # flags to analyze fitting stop
    out_of_bound::Bool = false
    #constraints_func_error = false

    function constraints_func(x, g) # testing grad methods
        # function constraints_func(x) # testing grad methods
        # this part is necessary to understand the difference between
        # "stop out of bounds" and "stop because of function call error"
        loss_value = try
            loss_func(x)
        catch e
            @warn "Error when call loss_func($x)"
            throw(e)
        end

        if (loss_value < 0.) && (scan_func(x) > scan_bound)
            out_of_bound = true
            throw(NLopt.ForcedStop())
        elseif length(g) > 0
          if isa(loss_grad, Function)
            g .= loss_grad(x)
          elseif loss_grad == :AUTODIFF
            ForwardDiff.gradient!(g, loss_func, x)
          elseif loss_grad == :FINITE
            Calculus.finite_difference!(loss_func,x,g,:central)
          end
        end

        return loss_value
    end

    function obj_func(x, g)
        if length(g) > 0
          if isa(scan_grad, Function)
            g .= scan_grad(x)
          elseif scan_grad == :AUTODIFF
            ForwardDiff.gradient!(g, scan_func, x)
          elseif scan_grad == :FINITE
            Calculus.finite_difference!(scan_func, x, g, :central)
          end
        end
        return scan_func(x)
    end

    # constrain optimizer
    opt = Opt(:LN_AUGLAG, n_theta)
    ftol_abs!(opt, scan_tol)

    max_objective!(
        opt,
        obj_func
        )

    local_optimizer!(opt, local_opt)
    maxeval!(opt, max_iter)

    # inequality constraints
    inequality_constraint!(
        opt,
        constraints_func,
        1e-3 # XXX: magic number, loss_tol
    )

    # version 1: internal :LN_AUGLAG box constrains
    opt.lower_bounds = [tb[1] for tb in theta_bounds]
    opt.upper_bounds = [tb[2] for tb in theta_bounds]

    # version 2: creating lower and upper bounds as general constraints
    #=
    [ inequality_constraint!(
        opt,
        (x, g) -> x[i] - theta_bounds[i][2],
        0.
    ) for i in 1:n_theta ]
    [ inequality_constraint!(
        opt,
        (x, g) -> theta_bounds[i][1] - x[i],
        0.
    ) for i in 1:n_theta ]
    =#

    # start optimization
    (optf, optx, ret) = optimize(opt, theta_init)

    if ret == :FORCED_STOP && !out_of_bound
        pp = ProfilePoint[]
        res = (nothing, pp, :LOSS_ERROR_STOP)
    elseif ret == :MAXEVAL_REACHED
        pp = ProfilePoint[]
        res = (nothing, pp, :MAX_ITER_STOP)
    elseif (ret == :FORCED_STOP || ret == :FAILURE) && out_of_bound # successfull result
        pp = ProfilePoint[]
        res = (nothing, pp, :SCAN_BOUND_REACHED)
    elseif ret == :FTOL_REACHED # successfull result
        loss = loss_func(optx)
        pp = [ ProfilePoint(optf, loss, optx, ret, nothing) ]
        res = (optf, pp, :BORDER_FOUND_BY_SCAN_TOL)
    else
        # this part is reached in case of :FAILURE return code
        # we have seen it in case of gradient methods
        pp = ProfilePoint[]
        res = (nothing, pp, :UNKNOWN_STOP)
    end

    return res
end

function get_right_endpoint(
    theta_init::Vector{Float64}, # initial point of parameters
    theta_num::Int, # number of parameter to scan
    loss_func::Function, # lambda(theta) - labmbda_min - delta_lambda
    method::Val{:CICO_ONE_PASS}; # function works only for method ONE_PASS;

    theta_bounds::Vector{Tuple{Float64,Float64}} = fill(
        (-Inf, Inf), length(theta_init)
        ),
    scan_bound::Float64 = 9.0,
    scan_tol::Float64 = 1e-3,
    scan_rtol::Float64 = 0.,
    loss_tol::Float64 = 0.,
    local_alg::Symbol = :LN_NELDERMEAD,
    max_iter::Int = 10^5,
    loss_grad::Union{Function, Symbol} = :EMPTY,
    #kwargs...
)
    # checking arguments
    if theta_num > length(theta_init)
        throw(DomainError(theta_num, "theta_num exceed theta dimention"))
    end

    scan_func(theta::Vector) = theta[theta_num]
    scan_grad(theta::Vector) = begin
        res = zeros(length(theta_init))
        res[theta_num] = 1.
        return res
    end

    get_right_endpoint(
        theta_init,
        scan_func,
        loss_func,
        method;

        theta_bounds,
        scan_bound,
        scan_tol,
        scan_rtol,
        loss_tol,
        local_alg,
        max_iter,
        scan_grad,
        loss_grad
        #kwargs...
    )
end


text:
using Documenter, MicrobiomePlots

makedocs(
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "MicrobiomePlots.jl",
    pages = [
        "Home" => "index.md",
        "Recipes" => "recipes.md",
        "Contributing" => "contributing.md"
    ],
    authors = "Kevin Bonham, PhD"
)

deploydocs(
    repo = "github.com/BioJulia/MicrobiomePlots.jl.git",
    push_preview=true
    )

max_issues_repo_name: BioJulia/MicrobiomePlots

max_forks_repo_name: BioJulia/MicrobiomePlots


text:
using MultivariateStats: cov_whitening
using NeXLUncertainties

"""
    olssvd(y::AbstractVector{N}, a::AbstractMatrix{N}, sigma::N, xLabels::Vector{<:Label}, tol::N=convert(N,1.0e-10))::UncertainValues where N <: AbstractFloat

Solves the ordinary least squares problem a⋅x = y for x using singular value decomposition for AbstractFloat-based types.
"""
function olssvd(
    y::AbstractVector{N},
    a::AbstractMatrix{N},
    sigma::N,
    xLabels::Vector{<:Label},
    tol::N = convert(N, 1.0e-10)
)::UncertainValues where {N<:AbstractFloat}
    f = svd(a)
    mins = tol * maximum(f.S)
    fs = [s > mins ? one(N) / s : zero(N) for s in f.S]
    # Note: w*a = f.U * Diagonal(f.S) * f.Vt
    genInv = f.V * Diagonal(fs) * transpose(f.U)
    #cov = sigma^2*f.V*(Diagonal(fs)*Diagonal(fs))*Transpose(f.V)
    cov = sigma^2 * genInv * transpose(genInv)
    return uvs(xLabels, genInv * y, cov)
end

"""
    olspinv(y::AbstractVector{N}, a::AbstractMatrix{N}, v::Matrix{N}, xlabels::Vector{<:Label}, tol::N=convert(N,1.0e-10))::UncertainValues

Solves the ordinary least squares problem a⋅x = y for x using the pseudo-inverse for AbstractFloat-based types.
"""
function olspinv(
    y::AbstractVector{N},
    a::AbstractMatrix{N},
    sigma::N,
    xLabels::Vector{<:Label},
    tol::N = convert(N, 1.0e-10),
)::UncertainValues where {N<:AbstractFloat}
    genInv = pinv(a, rtol = 1.0e-6)
    return uvs(xLabels, genInv * y, sigma * genInv * transpose(genInv))
end

"""
    glspinv(y::AbstractVector{N}, a::AbstractMatrix{N}, v::Matrix{N}, xlabels::Vector{<:Label}, tol::N=convert(N,1.0e-10))::UncertainValues

Solve the generalized least squares problem y = x β + ϵ for β where ϵ ~ v σ using the pseudo-inverse.

β = (xᵀv⁻¹x)⁻¹xᵀv⁻¹y where a⁻¹ => pinv(a)

cov[β] = (xᵀv⁻¹x)⁻¹
"""
function glspinv(
    y::AbstractVector{N},
    x::AbstractMatrix{N},
    v::AbstractMatrix{N},
    xLabels::Vector{<:Label},
    tol::N = convert(N, 1.0e-10),
)::UncertainValues where {N<:AbstractFloat}
    txiv = transpose(x) * pinv(v)
    yp, ixp = txiv * y, pinv(txiv * x)
    return uvs(xLabels, ixp * yp, ixp)
end

"""
    glsinv(y::AbstractVector{N}, a::AbstractMatrix{N}, v::Matrix{N}, xlabels::Vector{<:Label}, tol::N=convert(N,1.0e-10))::UncertainValues

Solve the generalized least squares problem y = x β + ϵ for β where ϵ ~ v σ using the inverse.
"""
function glsinv(
    y::AbstractVector{N},
    x::AbstractMatrix{N},
    v::AbstractMatrix{N},
    xLabels::Vector{<:Label},
    tol::N = convert(N, 1.0e-10),
)::UncertainValues where {N<:AbstractFloat}
    txiv = transpose(x) * inv(v)
    yp, ixp = txiv * y, inv(txiv * x)
    return uvs(xLabels, ixp * yp, ixp)
end

"""
    glssvd(y::AbstractVector{N}, a::AbstractMatrix{N}, cov::Matrix{N}, xlabels::Vector{<:Label}, tol::N=convert(N,1.0e-10))::UncertainValues

Solves the generalized least squares problem y = x β + ϵ for β using covariance whitening and the ordinary least squares pseudo-inverse.
"""
function glssvd(
    y::AbstractVector{N},
    x::AbstractMatrix{N},
    cov::AbstractMatrix{N},
    xLabels::Vector{<:Label},
    tol::N = convert(N, 1.0e-10),
)::UncertainValues where {N<:AbstractFloat}
    checkcovariance!(cov)
    w = cov_whitening(Matrix(cov))
    #olssvd(w*y, w*a, 1.0, xLabels, tol)
    return olspinv(w * y, w * x, 1.0, xLabels, tol)
end

"""
    glschol(y::AbstractVector{N}, a::AbstractMatrix{N}, cov::Matrix{N}, xlabels::Vector{<:Label}, tol::N=convert(N,1.0e-10))::UncertainValues

Solves the generalized least squares problem y = x β + ϵ for β using covariance whitening and the ordinary least squares pseudo-inverse.
"""
function glschol(
    y::AbstractVector{N},
    x::AbstractMatrix{N},
    cov::AbstractMatrix{N},
    xLabels::Vector{<:Label},
    tol::N = convert(N, 1.0e-10),
)::UncertainValues where {N<:AbstractFloat}
    checkcovariance!(cov)
    w = inv(cholesky(cov).L)
    return olspinv(w * y, w * x, 1.0, xLabels, tol)
end

"""
    wlssvd(y::AbstractVector{N}, a::AbstractMatrix{N}, cov::AbstractVector{N}, xlabels::Vector{<:Label}, tol::N=convert(N,1.0e-10))::UncertainValues

Solves the weighted least squares problem y = x β + ϵ  for β using singular value decomposition for AbstractFloat-based types.
"""
function wlssvd(
    y::AbstractVector{N},
    x::AbstractMatrix{N},
    cov::AbstractVector{N},
    xLabels::Vector{<:Label},
    tol::N = convert(N, 1.0e-10),
)::UncertainValues where {N<:AbstractFloat}
    w = Diagonal([sqrt(1.0 / cv) for cv in cov])
    return olssvd(w * y, w * x, 1.0, xLabels, tol)
end

"""
    wlspinv(y::AbstractVector{N}, a::AbstractMatrix{N}, cov::AbstractVector{N}, xlabels::Vector{<:Label}, tol::N=convert(N,1.0e-10))::UncertainValues

Solves the weighted least squares problem a⋅x = y for x using singular value decomposition for AbstractFloat-based types.
"""
function wlspinv(
    y::AbstractVector{N},
    a::AbstractMatrix{N},
    cov::AbstractVector{N},
    xLabels::Vector{<:Label},
    tol::N = convert(N, 1.0e-10),
)::UncertainValues where {N<:AbstractFloat}
    w = Diagonal([sqrt(1.0 / cv) for cv in cov])
    return olspinv(w * y, w * a, 1.0, xLabels, tol)
end

"""
    wlspinv2(y::AbstractVector{N}, a::AbstractMatrix{N}, cov::AbstractVector{N}, covscales::AbstractVector{N}, xLabels::Vector{<:Label}, tol::N=convert(N,1.0e-10))::UncertainValues where N <: AbstractFloat

Solves the weighted least squares problem a⋅x = y for x using singular value decomposition for AbstractFloat-based types.  Rescales the rescaleCovariances
according to `covscales`.
"""
function wlspinv2(
    y::AbstractVector{N},
    a::AbstractMatrix{N},
    cov::AbstractVector{N},
    covscales::AbstractVector{N},
    xLabels::Vector{<:Label},
    tol::N = convert(N, 1.0e-10),
)::UncertainValues where {N<:AbstractFloat}
    function rescaleCovariances(
        uvs::UncertainValues,
        covscales::AbstractVector{N},
    )::UncertainValues
        cov = copy(uvs.covariance)
        # rwgts = map(sqrt,wgts)
        a = Base.axes(cov)
        for r in a[1], c in a[2]
            cov[r, c] *= covscales[r] * covscales[c]
        end
        return UncertainValues(uvs.labels, uvs.values, cov)
    end
    w = Diagonal([sqrt(1.0 / cv) for cv in cov])
    return rescaleCovariances(olspinv(w * y, w * a, 1.0, xLabels, tol), covscales)
end

"""
    simple_linear_regression(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})::Tuple{Real, Real}

Peforms simple linear regression. Returns the `( slope, intercept, r, t )` of the unweighted best fit line through
the data `x` and `y`.  ( y = slope*x + intercept ), r is the correlation coefficient and t is the t-statistic.
"""
function simple_linear_regression(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
    sx, sy, n = sum(x), sum(y), length(x)
    ssxy, ssxx, ssyy = (dot(x, y) - sx*sy/n), (dot(x, x) - sx*sx/n), (dot(y, y) - sy*sy/n)
    slope, r = ssxy / ssxx, ssxy/sqrt(ssxx*ssyy)
    return ( slope, (sy - slope*sx)/n, r, r/sqrt((1.0-r^2)/(n-2)) )
end


text:
using LsqFit
using Optim
using MicroLogging
using PyPlot
using Roots

#Fitting a biased Lorentzian to amplitude data.
"""
Returns a five parameter, biased Lorentizian function

Args:
  x: domain of the returned function
  p: a vector of length 5 with function params --
      (amplitude, frequency, line width, linear skew, offset)
Returns:
    The computed function.

See appendix E.4 of Gao\'s thesis (2008)
"""
function bias_lorentzian(x, p)
  return p[1] ./ ((x .- p[2]).^2 .+ (p[3]/2)^2) .+ p[4]*(x .- p[2]) .+ p[5];
end

"""
Fit frequency-amplitude data to a skewed Lorentzian using nonlinear least-squares.
Model is y = a / ((x - b)^2 + (c/2)^2) + d*(x - b) + e

Args:
  freq: Frequency data.
  data: Data to be fitted.
  vars: Variance of data. Defaults to empty.
Returns:
  fit: LsqFitResult object.
  fit_function(x): Biased Lorentzian at fitted parameters.
"""
function fit_biased_lorentzian(freq, data, vars=[])

  BL_fit_dict(p) = Dict("a" => p[1], "b" => p[2], "c" => p[3],
                                "d" => p[4], "e" => p[5])
  p0 = initial_guess_blorentz(freq, data)
  return generic_fit(freq, data, bias_lorentzian, p0, BL_fit_dict,
    "y = a / ((x - b)^2 + (c/2)^2) + d*(x - b) + e", yvars=vars)
end

function initial_guess_blorentz(xpts, ypts)
  """
  Returns a length five vector with an initial guess of
  Lorentizian function parameters.
  """
  # find offset
  e = median(ypts)
  if abs(maximum(ypts) - median(ypts)) <= abs(minimum(ypts) - median(ypts))
        idx = findmin(ypts)[2]
        direction = -1
    else
        idx = argmax(ypts)
        direction = 1
  end
  # center frequency
  b = xpts[idx]
  half = abs(median(ypts) + ypts[idx]) / 2.
  if direction == 1
    idx_l = findfirst(x -> x > half, ypts)
    idx_r = findlast(x -> x > half, ypts)
  else
    idx_l = findfirst(x -> x < half, ypts)
    idx_r = findlast(x -> x < half, ypts)
  end
  c = abs(xpts[idx_l] - xpts[idx_r])
  a =  direction * c^2 * abs(ypts[idx] - e) / 4
  # slope of linear skew
  d = (ypts[end] - ypts[1])/(xpts[end] - xpts[1])
  return [a, b, c, d, e];
end

struct CircleFitParams
  """Container for the result of a circle fit.
    f0: Resonator center frequency.
    Qi: Resonator internal quality factor.
    Qc: Resonator coupling quality factor.
    ϕ: Impedance mismatch angle.
    τ: System phase delay.
    α: System gain angle.
    A: System gain amplitude.
  """
  f0
  Qi
  Qc
  ϕ
  τ
  α
  A
end

struct CircleFitResult
  fit_params::CircleFitParams
  sq_error::Float64
  SNR::Float64
  errors::CircleFitParams
  fit_curve::Function
end

#Fitting to the resonance circle of a quarter-wave resonator
function fit_resonance_circle(freq::Vector{T}, data::Vector{Complex{T}}; kwargs...) where T <: AbstractFloat
  """
  Fits complex-valued data

  Args:
    freq: Frequency data.
    data: Complex S21 data.
    kwargs: Use these if you already know something about the system and would like to override the
      following properties
      τ: Cable delay
      α: System phase shifted
      A: Overall system gain
  Returns:
    Result: The result of the fit.
  """
  @assert length(freq) == length(data) "Frequency and Data vectors must have same length!"
  @assert length(data) > 20 "Too few points."

  kwdict = Dict(kwargs)

  function χ2(data, R, xc, yc)
    x = real.(data)
    y = imag.(data)
    d = sqrt.( (x-xc).^2 + (y-yc).^2 ) - R
    return sum(d.^2)
  end

  #Fit to cable delay
  haskey(kwdict, :τ) ? τ = kwdict[:τ] : τ = fit_delay(freq, data)

  Sp = exp.(1im * 2. * π * freq * τ) .* data
  #Get best-fit circle and translate to origin.
  R, xc, yc = fit_circle(real.(Sp), imag.(Sp))

  #@assert χ2 < 10 "Could not fit circle to delay-corrected data. χ² = $(χ2)."
  @debug "Fit circle to delay corrected data with χ² = $(χ2(Sp, R, xc, yc))."

  #Translate circle to origin and fit overall phase delay and scaling
  St = Sp - (xc + 1im * yc)
  _, _, θ = fit_phase(freq, St)
  #find f → ∞ point -- this part seems fragile...
  β = mod(θ + π, π)
  P = xc + R*cos(β) + 1im*(yc + R*sin(β))

  @debug "Found scale parameters A = $(abs(P)), α = $(angle(P))."

  haskey(kwdict, :A) ? A = kwdict[:A] : A = abs(P)
  haskey(kwdict, :α) ? α = kwdict[:α] : α = angle(P)

  #Calibrate out α, A and get impedance mismatch angle
  Sc = Sp .* exp.(-1im * α) / A
  Rc, xcc, ycc = fit_circle(real.(Sc), imag.(Sc))

  #@assert χ2 < 10 "Could not fit circle to calibrated data. χ² = $(χ2)."
  @debug "Fit circle to calibrated data with χ² = $(χ2(Sc, Rc, xcc, ycc))."

  N = length(data)
  r = sqrt.((real(Sc) - xcc).^2 + (imag(Sc) - ycc).^2)
  SNR = Rc / sqrt(sum((r - Rc).^2)/(N-1))

  @debug "Data has an SNR of $(SNR)."

  ϕ = -asin(ycc/Rc)
  @debug "Found impedance mismatch angle: ϕ = $(ϕ)"
  #Final fit to phase to extract resonant frequency and Q
  Sct = Sc - (xcc + 1im * ycc)
  f0, Qr, _ = fit_phase(freq, Sct)

  Qc_cplx = Qr/(2 * Rc * exp(-1im * ϕ))
  Qi = 1/(1/Qr - real(1/Qc_cplx))
  Qc = 1/abs(1/Qc_cplx)

  @debug "Found quality factors: Qᵢ = $(Qi), Qc = $(Qc)."

  fit_params = CircleFitParams(f0, Qi, Qc, ϕ, τ, α, A)

  refined_params, errors = refine_circle_fit(freq, data, fit_params)

  χ0 = sum(abs.(data - lorentzian_resonance(fit_params, freq)).^2)
  χ1 = sum(abs.(data - lorentzian_resonance(refined_params, freq)).^2)
  @debug "Circle fit found χ² = $(χ0), refined to χ² = $(χ1)"

  if χ0 > χ1
    params = refined_params
    χ = χ1
  else
    params = fit_params
    χ = χ0
  end

  return CircleFitResult(params,
                          χ,
                          SNR,
                          errors,
                          x->lorentzian_resonance(fit_params, x))
end

function refine_circle_fit(freq, data, fit_params)
  """Refine circle fit and get errors using LsqFit.jl"""

  model(x, p) = abs.(lorentzian_resonance(p, x))

  initial_guess = [fit_params.f0, fit_params.Qi, fit_params.Qc,
                  fit_params.ϕ, fit_params.τ, fit_params.α, fit_params.A]

  fit = curve_fit(model, freq, abs.(data), initial_guess)
  errors = margin_error(fit)

  refined_params = CircleFitParams(fit.param[1], fit.param[2], fit.param[3],
                fit.param[4], fit.param[5], fit.param[6], fit.param[7])

  errs = CircleFitParams(errors[1], errors[2], errors[3],
                errors[4], errors[5], errors[6], errors[7])

  return (refined_params, errs)
end

function fit_phase(freq, data)
  """Fit phase of resonance.
  Args:
    freq: Frequency data.
    data: Complex S21 data.
  Returns:
    f0: Center frequency.
    Q: Quality factor.
    Θ0: Offset phase angle.
  """
  model(x, p) = p[1] + 2. * slope * atan.(2 * p[2] *(1. - x / p[3]))
  ϕ = unwrap(angle.(data))
  #guesses for initial parameters
  idx = findmin(abs.(ϕ - mean(ϕ)))[2]
  if mean(ϕ[1:9]) > mean(ϕ[end-9:end])
    j = findfirst(x -> x - ϕ[idx] < π/2., ϕ)
    k = findfirst(x -> x - ϕ[idx] < -π/2., ϕ)
    slope = 1
  else
    j = findfirst(x -> x - ϕ[idx] > π/2., ϕ)
    k = findfirst(x -> x - ϕ[idx] > -π/2., ϕ)
    slope = -1
  end
  Qguess = freq[idx]/abs.(freq[j] - freq[k])
  fit = curve_fit(model, freq, ϕ, [ϕ[idx], Qguess, freq[idx]])
  χ = sum(fit.resid.^2)

  @debug "Frequency vs. Phase fit found: f₀ = $(fit.param[3]), Q = $(fit.param[2]), θ₀ = $(fit.param[3])"
  @debug "Phase fit χ² = $(χ)"
  #@assert χ < 10 "Could not fit phase: χ² = $(χ)"

  return fit.param[3], fit.param[2], fit.param[1]
end

function fit_delay(freq, data)
  """Fit overall phase delay imposed by system transfer function.

  Args:
    freq: Frequency data.
    data: Complex S21 data.
  Returns:
    τ: Linear phase delay
  """
  function delay_model(x)
    ddata = data .* exp.(2. * π * 1im * freq * x)
    R, xc, yc = fit_circle(real.(ddata), imag.(ddata))
    d = sqrt.((real.(ddata) - xc).^2 + (imag.(ddata) - yc).^2) - R
    return sum(d.^2)
  end

  ϕ = unwrap(angle.(data))
  linfit(x,p) = p[1] + x * p[2]
  fit = curve_fit(linfit, freq, ϕ, [mean(ϕ), 0])
  ϕ0 = fit.param[2]
  result = optimize(delay_model, -1.5*abs(ϕ0), 1.5*abs(ϕ0)) #would prefer 1D gradient descent
  τ = Optim.minimizer(result)

  χ2 = delay_model(τ)
  @debug "Cable delay fit found: $(τ) with χ² = $(χ2)."
 #@assert χ2 < 10 "Could not calibrate out cable delay: χ² = $(χ2)."

  return τ
end

function lorentzian_resonance(p::CircleFitParams, f)
  """
  Return a resonance model in S21 amplitude over the range [x] and with
  parameters [p].

  See appendix E of Gao 2008 p. 155 Eq E.1
  """
  Qc_cplx = 1 / ((1/p.Qc)*exp(1im*p.ϕ))
  Q = 1/(1/p.Qi + real(1/Qc_cplx))
  return p.A*exp(1im*p.α).*exp.(-2π*1im*f*p.τ).*(1 - (Q/p.Qc)*exp(1im*p.ϕ)./(1 + 2*1im*Q*(f/p.f0 - 1)));
end

function lorentzian_resonance(p::Array, f)
  return lorentzian_resonance(CircleFitParams(p[1], p[2], p[3], p[4], p[5], p[6], p[7]), f)
end

function simulate_resonance(kwargs...)
  """
  Return synthetic data to test fit functions

  Returns: xpts, ypts, curve_params
  """
  f0 = 6.5 + randn()
  Qc = 120000 + randn()*2000
  Qi = 6e5 + randn()*1e4
  τ = 1.7
  ϕ = 0.17 * π
  Q = 1 ./ (1/Qi + real(1 ./ Qc*exp.(-1im*ϕ)))
  α = 0.35 * π
  A = 0.23
  df = f0/Q;
  freqs = range(f0 - 6*df, stop=f0 + 6*df, length=401)
  p = CircleFitParams(f0, Qi, Qc, ϕ, τ, α, A)
  data = lorentzian_resonance(p, freqs);
  # add some noise
  rand_phase = 2π*0.001*randn(length(freqs));
  rand_amp = 0.002*randn(length(freqs));
  data = data.*exp(-1im*rand_phase).*(1+rand_amp);

  return freqs, data, p
end

function fit_circle(x, y; refine=false)
  """Algebraic fit of (x,y) points to a circle. See:
      N. Chernov and C. J. Lesort, Least Squares Fitting of Circles,
        Journal of Mathematical Imaging and Vision, 23: 239-252, 2005.

      Args:
        x, y: X and Y points to be fitted to a circle.
      Returns:
        R: Circle radius.
        xc: x-coordinate of circle center.
        yc: y-coordinate of circle center.
  """
  @assert length(x) == length(y) "X and Y vector lengths must be equal for circle fit!"
  n = length(x)
  z = x.^2 + y.^2
  Mx = sum(x)
  My = sum(y)
  Mz = sum(z)
  Mxx = sum(x.^2)
  Myy = sum(y.^2)
  Mzz = sum(z.^2)
  Mxy = sum(x .* y)
  Mxz = sum(x .* z)
  Myz = sum(y .* z)
  M = [Mzz Mxz Myz Mz; Mxz Mxx Mxy Mx; Myz Mxy Myy My; Mz Mx My n];
  B = [0 0 0 -2; 0 1 0 0; 0 0 1 0; -2 0 0 0]

  F(η) = det(M - η*B)
  #We use Newton's method to guarantee convergence to the smallest positive eigenvalue
  ηs = newton(F, 0.0)
  ev = nullspace(M - ηs*B)

  xc = -ev[2]/2/ev[1]
  yc = -ev[3]/2/ev[1]
  R = sqrt(ev[2]^2 + ev[3]^2 - 4*ev[1]*ev[4]) / 2 ./ abs(ev[1])

  return R, xc, yc
end

function fit_circle_LM(x, y, initial_guess)
    """Iterative algorithm for refining a circle fit, using radial weighting."""

    function resid(p)
        C = sqrt.((x-p[2]).^2 + (y-p[3]).^2)
        W = 1 ./ sqrt.((p[2]-x).^2 + (p[3]-y).^2)
        return (p[1] - C) .* W
    end

    T = eltype(x)

    result = LsqFit.lmfit(resid, initial_guess, T[])

    return result.param[1], result.param[2], result.param[3]
end