text:
module GlobalEpistasis

export prepdata, nonepistatic_model, fit, boot, boot_stats, predict, cvpredict, lrt

using NLopt
using DataFrames

function prepdata(df, seqname, kind, wt, yname; delim = '-', cname = nothing, vname = nothing, condition_type = :categorical)
    n = nrow(df)
    I = collect(1:n)
    J = ones(Int64, n)
    j = 2
    aa = Char[]
    pos = Int64[]
    code = Dict{String, Int64}()
    if kind == :listofmuts
        muts = df[seqname]
        for i = 1:n
            if muts[i] != wt
                for k in split(muts[i], delim)
                    if !haskey(code, k)
                        code[k] = j
                        push!(aa, k[end])
                        push!(pos, parse(k[2:end-1]))
                        j += 1
                    end
                    push!(J, code[k])
                    push!(I, i)
                end
            end
        end
    elseif kind == :sequence
        seq = df[seqname]
        ns = length(wt)
        for i = 1:n
            for p = 1:ns
				a = seq[i][p]
	            if a != wt[p]
                    k = string(wt[p], p, a)
                    if !haskey(code, k)
                        code[k] = j
                        push!(aa, a)
                        push!(pos, p)
                        j += 1
                    end
                    push!(J, code[k])
                    push!(I, i)
                end
            end
        end
    end
    x = sparse(I, J, 1.0)
    g = trues(size(x,2))
    g[1] = false
    ham = ceil.(Int64, vec(sum(x,2))-1)

    if any(ismissing.(df[yname]))
		error("missing values in phenotype")
	end
	y = collect(Missings.skipmissing(df[yname]))
    dout = Dict(:x => x, :y => y, :code => code, :g => g,
                :ham=> ham, :pos => pos, :aa => aa, :wt => wt)
    if vname != nothing
        dout[:v] = collect(Missings.replace(df[vname], 0.0))
		if any(dout[:v] .< 0.0)
			error("negative variances")
		end
    end
	if cname != nothing
        if any(ismissing.(df[cname]))
			error("missing values in conditions")
		end
		c = collect(Missings.skipmissing(df[cname]))
		if condition_type == :categorical
			lev1 = levels(c)[1]
			for lev in levels(c)[2:end]
				x = hcat(x, sparse(c .== lev))
				push!(g, false)
			end
		elseif condition_type == :continuous
			c = convert(Vector{Float64}, c)
			x = hcat(x, c)
			push!(g, false)
		end
		dout[:c] = c
		dout[:x] = x
		dout[:g] = g
	end
    return dout
end

function nonepistatic_model(data; kwargs...)
    x = data[:x]
    y = data[:y]
    b = x \ y
    yhat = x * b
    mse = mean((y-yhat).^2)
    r2 = cor(y, yhat)^2
    ll = -length(y)/2 - 1/2*log(mse)
    m = Dict(:phi => yhat, :yhat => yhat, :b => b, :sigma2=>mse, :sigma2p => 1e-6,
        :method => :leftdiv, :data => data, :r2 => r2, :ll => ll)
    if haskey(data, :c)
	    g = data[:g]
	    m[:phiG] = view(x, :, g) * view(b, g)
    end
    if haskey(data, :v)
        return spmopt(m, false, true; kwargs...)
    else
        return m
    end
end

function monosplinebasis1(x, knots, k)
  @assert all(diff(knots) .> 0) "knots must be increasing"
  @assert length(knots) > 2 "must have at least 3 knots"
  nk = length(knots)
  na = nk+k-1
  I = zeros(na)
  I[1] = 1.0
  M = zeros(na)
  t = [[knots[1] for j in 1:k-1]; knots; [knots[end] for j in 1:k-1]]
  slope1, slope2 = monosplinebasisslopes(t, na-1, k)

  Mall = zeros(length(x), na)
  Iall = zeros(length(x), na)
  Iall[:,1] = 1.0
  # monosplinebasis!(Mall, Iall, slope1, slope2, x, t, k)
  for (i,xi) in enumerate(x)
    monosplinebasis1!(M, I, slope1, slope2, xi, t, k)
    Mall[i,:] = M
    Iall[i,:] = I
  end
  return (Mall, Iall, t, slope1, slope2)
end

function monosplinebasisslopes(t, na, kk)
  slope = zeros(2, na)
  for i in 1:na
    if t[i] < t[i+1]
      slope[1, i] = 1.0/(t[i+1]-t[i])
      slope[2, i] = 1.0/(t[i+1]-t[i])
    end
  end
  for k in 2:kk
    for i in 1:na
      if t[i] < t[i+k]
        slope[1, i] = (t[1]-t[i])*slope[1, i]*k/(k-1)/(t[i+k]-t[i])
        slope[2, i] = (t[end]-t[i])*slope[2, i]*k/(k-1)/(t[i+k]-t[i])
        if i+1 < na
          slope[1, i] += (t[i+k]-t[1])*slope[1, i+1]*k/(k-1)/(t[i+k]-t[i])
          slope[2, i] += (t[i+k]-t[end])*slope[2, i+1]*k/(k-1)/(t[i+k]-t[i])
        end
      end
    end
  end
  s1 = slope[1,1]
  s2 = slope[2,end]
  return (s1, s2)
end

function monosplinebasis1!(M::Vector{Float64}, I::Vector{Float64}, slope1::Float64, slope2::Float64, x::Float64, t::Vector{Float64}, kk::Int64)
  # M spline
  na = length(M) - 1
  M[2:na+1] = 0.0
  # extend basis to extrapolate linearly beyond knots
  if x <= t[1]
    # M[3:end] = 0.0
    M[2] = slope1
    I[3:end] = 0.0
    I[2] = slope1 * (x-t[1])
  elseif x >= t[end]
    # M[2:na] = 0.0
    M[na+1] = slope2
    I[2:na] = 1.0
    I[na+1] = 1.0 + slope2 * (x-t[end])
  else
    for i in 1:na
      ti1 = 1.0/(t[i+1]-t[i])
      if t[i] < x <= t[i+1]# || x[j] == t[1]
        M[i+1] = ti1
      # else
      #   M[i+1] = 0.0
      end
    end
    for k in 2:kk
      for i in 1:na
        tikk = k/(k-1)/(t[i+k]-t[i])
        m = 0.0
        if t[i] < x <= t[i+k]
          if x <= t[i+k-1]
            m += (x-t[i])*M[i+1]
          end
          if t[i+1] < x
            m += (t[i+k]-x)*M[i+2]
          end
          m *= tikk
        end
        M[i+1] = m
      end
    end
    # I spline
    nt = length(t)
    for i in 1:na
      j = 1
      for jj in 2:nt
        if x <= t[jj]
          j = jj-1
          break
        end
      end
      Itemp = 0.0
      if i > j
        Itemp = 0.0
      elseif i < j - kk + 1
        Itemp = 1.0
      else
        for m in i:j
          if m+1 <= na
            Itemp += 1/kk*((x-t[m])*M[m+1] + (t[m+kk+1]-x)*M[m+2])
          else
            Itemp += 1/kk*((x-t[m])*M[m+1])
          end
        end
      end
      I[i+1] = Itemp
    end
  end
end

function monosplinebasis!(M, I, slope1, slope2, xv, t, kk)
  # M spline
  na = size(M,2)-1
  M[:,2:na+1] = 0.0
  for (p, xp) in enumerate(xv)
#     Threads.@threads for p = 1:length(xv)
#     xp = xv[p]
    # extend basis to extrapolate linearly beyond knots
    if xp <= t[1]
      # M[3:end] = 0.0
      M[p,2] = slope1
      I[p,3:end] = 0.0
      I[p,2] = slope1 * (xp-t[1])
    elseif xp >= t[end]
      # M[2:na] = 0.0
      M[p,na+1] = slope2
      I[p,2:na] = 1.0
      I[p,na+1] = 1.0 + slope2 * (xp-t[end])
    else
      for i in 1:na
        ti1 = 1.0/(t[i+1]-t[i])
        if t[i] < xp <= t[i+1]# || x[j] == t[1]
          M[p,i+1] = ti1
        # else
        #   M[i+1] = 0.0
        end
      end
      for k in 2:kk
        for i in 1:na
          tikk = k/(k-1)/(t[i+k]-t[i])
          m = 0.0
          if t[i] < xp <= t[i+k]
            if xp <= t[i+k-1]
              m += (xp-t[i])*M[p,i+1]
            end
            if t[i+1] < xp
              m += (t[i+k]-xp)*M[p,i+2]
            end
            m *= tikk
          end
          M[p,i+1] = m
        end
      end
      # I spline
      nt = length(t)
      for i in 1:na
        j = 1
        for jj in 2:nt
          if xp <= t[jj]
            j = jj-1
            break
          end
        end
        Itemp = 0.0
        if i > j
          Itemp = 0.0
        elseif i < j - kk + 1
          Itemp = 1.0
        else
          for m in i:j
            if m+1 <= na
              Itemp += 1/kk*((xp-t[m])*M[p,m+1] + (t[m+kk+1]-xp)*M[p,m+2])
            else
              Itemp += 1/kk*((xp-t[m])*M[p,m+1])
            end
          end
        end
        I[p,i+1] = Itemp
      end
    end
  end
end

function spmopt(mi, estimate_alpha = haskey(mi, :a), estimate_beta = haskey(mi, :b);
        nk = 4, knots = :linear,
        a_upper_bound = get(mi, :a_upper_bound, [Inf, Inf]),
		a_lower_bound = get(mi, :a_lower_bound, [0.0, 0.0]),
        maxit = 1000000, alg = :LD_LBFGS, tol=1e-14, constrainbeta = true)

    data = mi[:data]
    x = data[:x]
    y = data[:y]

    n, nb = size(x)
    pi = zeros(0)
    lbounds = zeros(0)
    ubounds = zeros(0)

    estimate_sigma2 = haskey(data, :v)
    estimate_sigma2p = false
    if estimate_sigma2
        v = data[:v]
        estimate_sigma2p = any(v .== 0.0)
        push!(pi, log(mi[:sigma2]))
        push!(lbounds, -Inf)
        push!(ubounds, Inf)
        if estimate_sigma2p
            push!(pi, log(mi[:sigma2p]))
            push!(lbounds, -Inf)
            push!(ubounds, Inf)
        end
        iv0 = find(v .== 0.0)
        iv = find(v .!= 0.0)
    else
        iv0 = nothing
        iv = nothing
	    v = nothing
    end

    phi = deepcopy(mi[:phi])
    b = deepcopy(mi[:b])
    if estimate_alpha && !estimate_beta
        minphi, maxphi = extrema(phi)
        # rescale phi to have range 0..1
        b = b/(maxphi-minphi)
        b[1] = b[1] - minphi/(maxphi-minphi)
        phi = (phi-minphi)/(maxphi-minphi)
    end
    if estimate_alpha
        if haskey(mi, :knots)
            knots = mi[:knots]
        elseif knots == :linear
            knots = linspace(0, 1, nk)
        elseif knots == :quantile
            knots = unique(quantile(phi,  linspace(0, 1, nk)))
            if length(knots) < nk
                warn("knots removed  from quantile redundancy")
            end
        end
        M, I, t, slope1, slope2 = monosplinebasis1(phi, knots, 3)
        na = size(I,2)
        a = zeros(na)
        arange = length(pi) + (1:na)
        if haskey(mi, :a)
            #push!(pi, mi[:a][1])
            #append!(pi, log.(mi[:a][2:end]))
            a0 = copy(mi[:a])
            if a_upper_bound[1] < a0[2]
                a0[2] = a_upper_bound[1]
            end
            if a_upper_bound[2] < a0[end]
                a0[end] = a_upper_bound[2]
            end
            append!(pi, a0)
        else
            #push!(pi, 0.0)
            #append!(pi, log.(0.001*ones(na-1)))
            append!(pi, zeros(na))
        end
        append!(lbounds, [-Inf; a_lower_bound[1]; zeros(na-3); a_lower_bound[2]])
        append!(ubounds, [Inf; a_upper_bound[1]; Inf*ones(na-3); a_upper_bound[2]])
    else
        arange = nothing
        knots = nothing
        M = nothing
        I = nothing
        t = nothing
        slope1 = nothing
        slope2 = nothing
        a = nothing
    end
    if estimate_beta
        brange = length(pi) + (1:nb)
        append!(pi, b)
        append!(lbounds, -Inf*ones(nb))
        append!(ubounds, Inf*ones(nb))
    else
        brange = nothing
    end
    opt = Opt(alg, length(pi))
    lower_bounds!(opt, lbounds)
    upper_bounds!(opt, ubounds)
    maxeval!(opt, maxit)
    ftol_rel!(opt, tol)

    yhatp = ones(y)
    llmem = similar(y)
    gs2 = similar(y)
    gs2p = similar(y)
    rsy = similar(y)
    rsvi = similar(y)

    cb = (p, g) ->  spm_objective(p, g, x, y, v,
            estimate_sigma2, estimate_sigma2p, estimate_alpha, estimate_beta,
            arange, brange, M, I, t, slope1, slope2,
            phi, yhatp, llmem, gs2, gs2p, rsy, rsvi, iv0, iv)
    max_objective!(opt, cb)
	(ll, p, ret) = optimize!(opt, pi)
    #if estimate_alpha && estimate_beta
    #    a = p[arange]
    #    if a[2] < 0.01 || a[end] < 0.01
    #        display("reoptimizing with lower a bound 0.01")
    #        pi = p
    #        ai = [arange[2], arange[end]]
    #        lbounds[ai] = 0.01
    #        pi[ai] = 0.01
    #        lower_bounds!(opt, lbounds)
    #        (ll, p, ret) = optimize!(opt, pi)
    #        display("reoptimizing with lower a bound 0.0")
    #        pi = p
    #        lbounds[ai] = 0.0
    #        lower_bounds!(opt, lbounds)
    #        (ll, p, ret) = optimize!(opt, pi)
    #    end
    #end
    g = data[:g]
	if constrainbeta && estimate_beta && estimate_alpha && all(data[:ham] .<= 1) && haskey(data, :c) # single mutants multi conditions
		betabad = p[brange]
		phiE = x[:, .!g] * betabad[.!g]
	    constrain_betas!(p[arange], betabad, g, phiE, 0.1)
		p[brange] = betabad
		display("constraining betas")
    	#(ll, p, ret) = optimize!(opt, p)
		#constrain_betas!(p[arange], betabad, g, phiE, 0.1)
		#p[brange] = betabad
	end

    m = Dict(:data => data, :nlopt_return => ret, :ll => ll*n)
    if estimate_beta
        m[:b] = p[brange]
    else
        m[:b] = b
    end

    m[:phi] = x*m[:b]

    if estimate_alpha
        m[:a] = p[arange]
        m[:yhat] = monosplinebasis1(m[:phi], knots, 3)[2]*m[:a]
        m[:knots] = knots
        m[:a_upper_bound] = a_upper_bound
    else
        m[:yhat] = m[:phi]
    end

    m[:r2] = cor(y, m[:yhat])^2
	m[:rmse] = sqrt(mean((m[:yhat]-data[:y]).^2))
    if estimate_sigma2
        m[:sigma2] = exp(p[1])
        if estimate_sigma2p
            m[:sigma2p] = exp(p[2])
        end
    else
        m[:sigma2] = sum_kbn((y.-m[:yhat]).^2)/n
        m[:ll] = -n/2*log(m[:sigma2])
    end
	m[:beta] = DataFrame(:pos => data[:pos], :aa => data[:aa])
	m[:prediction] = DataFrame(:y => m[:data][:y], :yhat => m[:yhat])
	beta = m[:b][g]
    if estimate_alpha
		m[:beta][:b] = beta/mean(abs.(beta))
		m[:prediction][:phiG] = x[:, g] * m[:beta][:b]
	else
		m[:beta][:b] = beta
	end
	#sort!(m[:beta], :pos)
	if haskey(data, :c)
		m[:prediction][:c] = data[:c]
        m[:be] = m[:b][.!g]
	    m[:prediction][:phiE] = x[:, .!g] * m[:be]
    end

    return m
end

function constrain_betas!(a, b, g, phiE, delta = 0.0)
  if a[end] == 0.0 # curve is flat on right
    minbc = minimum(phiE)
    for i in find(g)
      if minbc + b[i] > 1
        b[i] = 1 - minbc - delta
      end
    end
  end
  if a[2] == 0.0 # curve is flat on left
    maxbc = maximum(phiE)
    for i in find(g)
      if maxbc + b[i] < 0
        b[i] = -maxbc + delta
      end
    end
  end
end

function spm_objective(p, g, x, y, v,
            estimate_sigma2, estimate_sigma2p, estimate_alpha, estimate_beta,
            arange, brange, M, I, t, slope1, slope2,
            yhat, yhatp, ll::Vector{Float64}, gs2, gs2p, rsy, rsvi::Vector{Float64}, iv0, iv)
    if estimate_beta
        A_mul_B!(yhat, x, view(p, brange))
    end
    if estimate_alpha
        #a .= exp.(view(p, arange))
        #a[1] = p[arange[1]]
        #for i = 2:length(a)
        #    ai = p[arange[i]]
        #    if ai > 20
        #        ai = 20
        #    end
        #    a[i] = exp(ai)
        #end
        a = view(p, arange)
        if estimate_beta
            monosplinebasis!(M, I, slope1, slope2, yhat, t, 3)
        end
        A_mul_B!(yhat, I, a)
        A_mul_B!(yhatp, M, a)
    end

    n, nb = size(x)
#     Threads.@threads for j = 1:n
    if estimate_sigma2
        s = exp(p[1])
        for j in iv
	        r = y[j]-yhat[j]
	        svi = 1.0/(s+v[j])
	        rsvi[j] = r*svi
	        gs2p[j] = 0.0
	        gs2[j] = (rsvi[j]*r-1)*svi*s/2
	        ll[j] = -rsvi[j]*r-log(s+v[j])
        end
        if estimate_sigma2p
            s2 = exp(p[2])
            svi2 = 1/(s+s2)
            lsvi2 = log(s+s2)
            for j in iv0
                r = y[j]-yhat[j]
                rsvi[j] = r*svi2
                gsvi = (rsvi[j]*r-1)*svi2/2
                gs2p[j] = gsvi*s2
                gs2[j] = gsvi*s
                ll[j] = -rsvi[j]*r-lsvi2
            end
            g[2] = sum_kbn(gs2p)
        end
        g[1] = sum_kbn(gs2)
    else
        for j = 1:n
            r = y[j]-yhat[j]
            rsvi[j] = r
            ll[j] = -r^2
        end
    end
    sll = sum_kbn(ll)/2/n
   # display(any(isnan.(yhat)))
#     if isnan(sll)
#     #    display("$s $s2")
#         display(p[brange])
#         display(p[arange])
#         display(find(isnan.(yhat)))
#         display(find(isnan.(p[brange])))
#         #display(find(isnan.(p[arange])))
#         error("nan ll")
#     end
    if estimate_beta
        rsy .= rsvi.*yhatp
        At_mul_B!(view(g, brange), x, rsy)
    end
    if estimate_alpha
        At_mul_B!(view(g, arange), I, rsvi)
        #for i = 2:length(a)
        #    g[arange[i]] = g[arange[i]]*a[i]
        #end
    end
    g .= g/n
	return sll

end

function fit(m0::Dict; kwargs...)
	if !haskey(m0, :data)
		data = m0
		m0 = nonepistatic_model(data; kwargs...)
	end
	if !haskey(m0, :a)
	    m0 = spmopt(m0, true, false; kwargs...)
	end
    spmopt(m0, true, true; kwargs...)
end

function boot(m; kwargs...)
    mboot = copy(m)
    mboot[:data] = copy(m[:data])
    if haskey(m[:data], :v)
        v = deepcopy(mboot[:data][:v])
        if haskey(m, :sigma2p)
            v[v.==0.0] = m[:sigma2p]
        end
        v .= sqrt.(v .+ m[:sigma2])
    else
        v = sqrt(m[:sigma2])
    end
    mboot[:data][:y] = randn(length(m[:yhat])).*v .+ m[:yhat]
    wt = m[:data][:ham] .== 0
    if all(m[:data][:y][wt] .== 0.0)
        mboot[:data][:y][wt] = 0.0
    end
    spmopt(mboot; kwargs...)
end

function boot(mi, i, search = false, tol=1e-4; kwargs...)
    mdls = []
    lli = mi[:ll]
    while length(mdls) < i
        mb = boot(mi; kwargs...)
        mb[:data] = mi[:data]
        if search
            m = spmopt(mb; kwargs...)
            dll = m[:ll]-mi[:ll]
            if dll > tol
                mi = m
                mdls = []
            else
                delete!(mb, :yhat)
                delete!(mb, :phi)
				delete!(mb, :prediction)
				if haskey(mi[:data], :c)
					delete!(mb, :phiE)
					delete!(mb, :be)
				end
				delete!(mb, :beta)
                push!(mdls, mb)
            end
            display("deltaLL $dll, i $(length(mdls))")
        else
            delete!(mb, :yhat)
            delete!(mb, :phi)
			delete!(mb, :prediction)
			if haskey(mi[:data], :c)
				delete!(mb, :phiE)
				delete!(mb, :be)
			end
			delete!(mb, :beta)
            push!(mdls, mb)
            display(length(mdls))
        end
    end
    display("total change in LL $(mi[:ll]-lli)")
    return mdls
end

#using GLM
#using DataFrames
function boot_stats(m, mb)
  bboot = DataFrame(b = Float64[], name = String[])
  data = m[:data]
  g = data[:g]
  names = ["$p$a" for (p,a) in zip(data[:pos],data[:aa])]
  bnorm = mean(abs.(m[:b][g]))
  aboot = DataFrame(a = Float64[], i = Int64[])
  i = collect(1:length(m[:a]))
	phi = [linspace(minimum(m[:phi]), 0, 101); 0:.01:1; linspace(1,maximum(m[:phi]), 101)]
    bootCI = DataFrame(phi = Float64[], yhat = Float64[], ddy = Float64[])
	phi1 = hcat(ones(length(m[:phi])), m[:phi])
  for mm in mb
	c = phi1\(data[:x] * mm[:b])

	append!(aboot, DataFrame(a = mm[:a]*bnorm, i = i))
	append!(bboot, DataFrame(b = mm[:b][g]/c[2]/bnorm, name = names))

	sp = monosplinebasis1(phi*c[2]+c[1], m[:knots], 3)
	dy = sp[1]*mm[:a]
    ddy = (dy[2:end].-dy[1:end-1])/(phi[2]-phi[1])
    append!(bootCI, DataFrame(phi = phi[1:end-1], yhat = (sp[2]*mm[:a])[1:end-1], ddy = ddy))
  end
  bbs = by(bboot, :name, d -> DataFrame(med = median(d[:b]),
										upper = quantile(d[:b], .975),
										lower = quantile(d[:b], .025),
										se = sqrt(var(d[:b])/length(d[:b]))))
  bbs = join(bbs, DataFrame(b = m[:beta][:b], name = names), on = :name)
  ben = bbs[:lower] .> 0
  del = bbs[:upper] .< 0
  neut = .!ben .& .!del
  bbs[:effect] = ["" for i = 1:nrow(bbs)]
  bbs[:effect][ben] = "beneficial"
  bbs[:effect][del] = "deleterious"
  bbs[:effect][neut] = "neutral"
  bbs[:pos] = data[:pos]
  bbs[:aa] = data[:aa]
  bben = sum(bbs[:lower].>0)
  bdel = sum(bbs[:upper].<0)
  l = sum(g)
  display("beta $(l), positive $bben, deleterious $bdel, neutral $(l-bben-bdel) ")
  ciw = mean(bbs[:upper].-bbs[:lower])
  display("beta average 95% CI width $ciw")
  display("beta average SE, $(sqrt(mean(bbs[:se].^2)))")

  aboot = by(aboot, :i, d -> DataFrame(med = median(d[:a]), upper = quantile(d[:a], .975), lower = quantile(d[:a], .025)))
  aboot = join(aboot, DataFrame(a = m[:a]*bnorm, i = i), on = :i)
	lower(x) = quantile(x, 0.025)
	upper(x) = quantile(x, 0.975)
	bootCI = by(bootCI, :phi, d -> DataFrame(yhat_upper = upper(d[:yhat]), yhat_lower = lower(d[:yhat]),
					ddy_upper = upper(d[:ddy]), ddy_lower = lower(d[:ddy])))
#    bootCI = aggregate(bootCI, :phi, [lower, upper])
	bootCI[:phiG] = (bootCI[:phi]-m[:b][1])/mean(abs.(m[:b][g]))
	#display(bootCI)
	ciw = mean(bootCI[:yhat_upper].-bootCI[:yhat_lower])
	display("g(phi) average CI $ciw")
	g2k = ([0,1]-m[:b][1])/mean(abs.(m[:b][g]))
	display("g(phi) boundary $g2k")
	if any(bootCI[:ddy_lower] .> 0.0)
		g2p = extrema(bootCI[:phiG][bootCI[:ddy_lower] .> 0.0])
		display("g(phi) second deriv+ CI $g2p")
	end
	if any(bootCI[:ddy_upper] .< 0.0)
		g2m = extrema(bootCI[:phiG][bootCI[:ddy_upper] .< 0.0])
		display("g(phi) second deriv- CI $g2m")
	end
	#g2p = extrema(bootCI[:phi][bootCI[:ddy_lower] .> 0.0])
	#display("g(phi) unscaled second deriv+ CI $g2p")
	#g2m = extrema(bootCI[:phi][bootCI[:ddy_upper] .< 0.0])
	#display("g(phi) unscaled second deriv- CI $g2m")
	if haskey(data, :c)
		lc = levels(data[:c])
		bootCI[:c] = lc[1]
		c = m[:b][.!g]
		gboot1 = deepcopy(bootCI)
		for i in 2:length(lc)
			gboot1[:phiG] = (gboot1[:phi]-c[1]-c[i])/bnorm
			gboot1[:c] = lc[i]
			append!(bootCI, gboot1)
		end
		minphiG, maxphiG = extrema(m[:prediction][:phiG])
		phiGrange = (bootCI[:phiG] .> minphiG) .& (bootCI[:phiG] .< maxphiG)
		bootCI = bootCI[phiGrange,:]
	end
	bootCI[:yhat] = monosplinebasis1(bootCI[:phi], m[:knots], 3)[2]*m[:a]

  return (bbs, aboot, bootCI)
end

import StatsBase.predict

function designmatrix_listofmuts(wt, seq, delim, code = Dict{String, Int64}())
    n = length(seq)
	I = collect(1:n)
    J = ones(Int64, n)
    j = 2
    for (i,s) in enumerate(seq)
        if s != wt
            for k in split(s, delim)
                if !haskey(code, k)
                    code[k] = j
                    j += 1
                end
                push!(I, i)
                push!(J, code[k])
            end
        end
    end
	return sparse(I, J, 1.0)
end

function designmatrix_sequence(wt, seq, pos, code = Dict{String, Int64}())
    n = length(seq)
	I = collect(1:n)
    J = ones(Int64, n)
    j = 2
    ns = length(seq[1])
    for (i,s) in enumerate(seq)
        for p = 1:ns
            if s[p] != wt[p]
                k = string(wt[p], pos[p], s[p])
                push!(I, i)
                push!(J, code[k])
            end
        end
    end
	return sparse(I, J, 1.0)
end

function designmatrix(data, wt, seq, kind, delim, pos)
    code = data[:code]
    n = length(seq)
	I = collect(1:n)
    J = ones(Int64, n)
    if kind == :listofmuts
        for (i,s) in enumerate(seq)
            if s != wt
                for k in split(s, delim)
                    push!(I, i)
                    push!(J, code[k])
                end
            end
        end
    elseif kind == :sequence
        ns = length(seq[1])
        for (i,s) in enumerate(seq)
            for p = 1:ns
                if s[p] != wt[p]
                    k = string(wt[p], pos[p], s[p])
                    push!(I, i)
                    push!(J, code[k])
                end
            end
		end
    else
        error("must be :listofmuts or :sequence")
    end
	return sparse(I, J, 1.0, n, size(data[:x], 2))
end

function predictx(m, x::SparseMatrixCSC{Float64,Int64})
    g = m[:data][:g]
    phi = x*m[:b]
    if haskey(m, :a)
        yhat = monosplinebasis1(phi, m[:knots], 3)[2]*m[:a]
	    phiG = x[:, g] * m[:beta][:b]
		return DataFrame(phi = phi, phiG = phiG, yhat = yhat)
    else
        return DataFrame(yhat = phi)
    end
end

function predict(m, seq, kind; wt = m[:data][:wt], delim = '-', pos = 1:length(seq[1]))
    x = designmatrix(m[:data], wt, seq, kind, delim, pos)
	predictx(m, x)
end

function cvpredict(data, fun, nfold; eachcondition = false, kwargs...)
  x = data[:x]
  g = data[:g]
	if eachcondition
	  nfold = ceil(Int64, mean(sum(x[:, g], 1)))
	  p = zeros(Int64, size(x,1))
	  for i in find(g)
		xi = find(view(x, :, i))
		l = length(xi)
		p[xi] = shuffle!(collect(1:l))
	  end
	else
		n = length(data[:y])
		p = ceil.(Int, collect(1:n)/n*nfold)
		shuffle!(p)
	end
    yhat = zeros(data[:y])
	function pred1(i)
		display(i)
		data1 = copy(data)
		traini = p .!= i
		data1[:x] = data[:x][traini,:]
		data1[:y] = data[:y][traini]
        if haskey(data, :v)
            data1[:v] = data[:v][traini]
        end
		if haskey(data, :c)
			data1[:c] = data[:c][traini]
		end
		mcv = fun(data1; kwargs...)
        testi = p .== i
		yhat[testi] = predictx(mcv, data[:x][testi, :])[:yhat]
	end
	map(pred1, 1:nfold)
	display("$nfold fold cross-validated correlation $(cor(yhat,data[:y]))")
	return yhat
end

function lrt(ii, m0, f1, f2)
    m1 = f2(m0)
	dl = m1[:ll]-m0[:ll]
    display("delta LL $dl")
	if dl == 0.0
		error("no change in likelihood")
	end
    m0b = pmap(f1, repeat([m0], inner=ii))
    m1b = pmap(f2, m0b)
    dl = map((m0,m1) -> m1[:ll] - m0[:ll], m0b, m1b)
    display("p = $(mean(dl .> m1[:ll]-m0[:ll]))")
    return dl
end

end # module


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    CrossValidation(k; shuffle=true, loss=Dict())

`k`-fold cross-validation. Optionally, `shuffle` the
data, and specify `loss` function  from `LossFunctions.jl`
for some of the variables.

    CrossValidation(partitioner; loss=Dict())

Generalization of k-fold cross-validation in which
the data is split using `partitioner`.

## References

* Geisser, S. 1975. [The predictive sample reuse method with applications]
  (https://www.jstor.org/stable/2285815)
* Burman, P. 1989. [A comparative study of ordinary cross-validation, v-fold
  cross-validation and the repeated learning-testing methods]
  (https://www.jstor.org/stable/2336116)
"""
struct CrossValidation{P<:AbstractPartitioner} <: AbstractErrorEstimator
  partitioner::P
  loss::Dict{Symbol,SupervisedLoss}
end

CrossValidation(partitioner::AbstractPartitioner; loss=Dict()) =
  CrossValidation{typeof(partitioner)}(partitioner, loss)

CrossValidation(k::Int; shuffle=true, loss=Dict()) =
  CrossValidation(UniformPartitioner(k, shuffle), loss=loss)

function Base.error(solver::AbstractLearningSolver,
                    problem::LearningProblem,
                    eestimator::CrossValidation)
  # retrieve problem info
  sdata = sourcedata(problem)
  ovars = outputvars(task(problem))
  partitioner = eestimator.partitioner
  loss  = eestimator.loss
  for var in ovars
    if var ∉ keys(loss)
      loss[var] = defaultloss(sdata[1,var])
    end
  end

  # folds for cross-validation
  folds  = subsets(partition(sdata, partitioner))
  nfolds = length(folds)

  solutions = pmap(1:nfolds) do k
    # source and target indices
    sinds = [ind for i in vcat(1:k-1, k+1:nfolds) for ind in folds[i]]
    tinds = folds[k]

    # source and target data
    train = view(sdata, sinds)
    hold  = view(sdata, tinds)

    # setup and solve sub-problem
    subproblem = LearningProblem(train, hold, task(problem))
    solve(subproblem, solver)
  end

  result = pmap(ovars) do var
    losses = map(1:nfolds) do k
      hold = view(sdata, folds[k])
      y = hold[var]
      ŷ = solutions[k][var]
      value(loss[var], y, ŷ, AggMode.Mean())
    end
    var => mean(losses)
  end

  Dict(result)
end

function Base.error(solver::AbstractEstimationSolver,
                    problem::EstimationProblem,
                    eestimator::CrossValidation)
  # retrieve problem info
  pdata = data(problem)
  pdomain = domain(problem)
  partitioner = eestimator.partitioner

  # folds for cross-validation
  folds  = subsets(partition(pdata, partitioner))
  nfolds = length(folds)

  result = []
  for (var, V) in variables(problem)
    # mappings from data to domain locations
    varmap = Dict(datloc => domloc for (domloc, datloc) in datamap(problem, var))

    # k-fold validation loop
    losses = pmap(1:nfolds) do k
      # training and holdout set
      train = [ind for i in vcat(1:k-1, k+1:nfolds) for ind in folds[i]]
      hold  = folds[k]

      # discard indices that were filtered out by mapping strategy (e.g. missing values)
      train = filter(in(keys(varmap)), train)
      hold  = filter(in(keys(varmap)), hold)

      if isempty(train)
        missing
      else
        # copy data to their locations in domain
        mapper = CopyMapper([varmap[ind] for ind in train])

        # setup and solve sub-problem
        subproblem = EstimationProblem(view(pdata, train), pdomain,
                                       var, mapper=mapper)
        solution = solve(subproblem, solver)

        # get true holdout values
        y = pdata[hold,var]

        # get solver estimate at holdout locations
        μ, _ = solution[var]
        ŷ = [μ[varmap[loc]] for loc in hold]

        # loss for fold
        mean((ŷ .- y).^2)
      end
    end

    push!(result, var => mean(skipmissing(losses)))
  end

  Dict(result)
end

ext: jl


text:
using Documenter, BioMedQuery
using Literate

# compile all examples in BioMedQuery/examples/literate_src into markdown and jupyter notebooks for documentation
for (root, dirs, files) in walkdir("examples/literate_src")
    for file in files
        Literate.notebook(joinpath(root,file), joinpath(@__DIR__, "src", "notebooks"))
    end
end

for (root, dirs, files) in walkdir("examples/literate_src")
    for file in files
        Literate.markdown(joinpath(root,file), joinpath(@__DIR__, "src", "examples"))
    end
end

makedocs()

deploydocs(
    deps   = Deps.pip("mkdocs==0.17.5", "mkdocs-material==2.9.4"),
    repo = "github.com/bcbi/BioMedQuery.jl.git",
    julia  = "0.7",
    osname = "linux"
)


text:
type Distribution
    tG::Tuple{Vararg{Array{Float64}}}
    td::Array{Float64}
    T::SparseMatrixCSC
    d::Array{Float64}
end

type temp
    error::Array{Float64,2}
    J::Array{Float64,2}
    foc::Expr
    ploc::Vector{Int}
    updateD::Function
    f2::Expr
end

type Model
    G::NGrid
    S::Array{Float64,2}
    SP::Array{Float64,2}
    X::Array{Float64,2}
    XP::Array{Float64,2}
    ProbWeights::Array{Float64,2}
    variables::Vector{ModVar}
    F::Function
    J::Function
    distribution::Distribution
    temp::temp
end

(M::Model)(s::Symbol,x::Array{Float64,2}) = M.G(M[s,0],x)

function display(M::Model)
    [println("\t",x) for x in M.temp.foc.args]
    println("S: ",names(State(M.variables)))
    println("P: ",names(Policy(M.variables)))
end

function Model(foc,states,policy,vars,params)
    slist,G,S,SP,ProbWeights    = parsestate(states)
    vlist,Dlist,dlist,plist     = parsevars(slist,policy,vars,params)
    f1,f2,j2,expects = parsefoc(foc,vlist,dlist,plist)

    vlist = vcat(vlist,Dlist)
    vlist = vcat(vlist,Future(f2,slist,vlist))

    X = hcat([((length(v.def.args)==4 && vtype(v)==Policy) ? S[:,find(slist,v.def.args[3])]*v.init : ones(length(G))*v.init[1]) for v in Future(vlist,false)]...)
    XP = zeros(length(G)*size(ProbWeights,2),length(Future(vlist)))

    tG = ndgrid(Vector{Float64}[s.val.x for s in slist]...)
    td = zeros(size(tG[1]))+1/length(tG[1])
    T   = spzeros(length(td),length(td))
    d  = zeros(length(G))
    D = Distribution(tG,td,T,d)

    vlist = vcat(slist,vlist)

    return Model(G,S,SP,X,XP,ProbWeights,vlist,
        eval(:( $(gensym(:F))(M::Model) = @fastmath $(buildfunc(addpweights!(inserthardloc(f2,vlist),size(ProbWeights,2)),:(M.temp.error))))),
        eval(:( $(gensym(:J))(M::Model,i::Int) = @fastmath $(buildJ(vec(addpweights!(inserthardloc(j2,vlist),size(ProbWeights,2))))))),
        D,
        temp(ones(length(G),length(Policy(vlist)))
            ,zeros(length(Policy(vlist)),length(Policy(vlist)))
            ,foc,
            Int[findfirst(names(State(vlist,false)),p)+size(S,2) for p in names(Policy(vlist))],
            eval(:($(gensym(:updateD))(M) = $(buildfunc(inserthardloc(addindex(subs(Expr(:vect,[v.def for v in Dlist]...),plist)),vlist),:(M.X),length(vlist)-length(slist)-length(Dlist)-length(Future(vlist)))))),
            :( $(gensym(:F))(M) = @fastmath $(buildfunc(addpweights!(inserthardloc(f2,vlist),size(ProbWeights,2)),:(M.temp.error))))
            )
        )
end

function getindex(M::Model,x::Symbol)
    if in(x,names(State(M.variables)))
        return M.variables[findfirst(names(State(M.variables)),x)].val
    else
        error("$x not found")
    end
end

function getindex(M::Model,x::Symbol,t::Int)
    ns = length(State(M.variables))
    for i ∈ 1:length((M.variables))
        if  t==-1 && vtype(M.variables[i])==Endogenous && M.variables[i].name==x
            return M.S[:,i]
        elseif vtype(M.variables[i])==Aggregate && M.variables[i].name==x && M.variables[i].val.target[2]==t
            return M.X[:,i-ns]
        elseif  t==0 && vtype(M.variables[i]) <: Stochastic && M.variables[i].name==x
           return  M.S[:,i]
        elseif  t==0 && vtype(M.variables[i]) <: Union{Policy,Static}  && M.variables[i].name==x
           return  M.X[:,i-ns]
        end
   end
end

function setindex!(M::Model,X::Vector{Float64},x::Symbol,t::Int)
    ns = length(State(M.variables))
    for i ∈ 1:length((M.variables))
        if  t==-1 && vtype(M.variables[i])==Endogenous && M.variables[i].name==x
            return M.S[:,i] = X
        elseif vtype(M.variables[i])==Aggregate && M.variables[i].name==x && M.variables[i].val.target[2]==t
            return M.X[:,i-ns] = X
        elseif  t==0 && vtype(M.variables[i]) <: Stochastic && M.variables[i].name==x
           return  M.S[:,i] = X
        elseif  t==0 && vtype(M.variables[i]) <: Union{Policy,Static}  && M.variables[i].name==x
           return  M.X[:,i-ns] = X
        end
   end
end


text:
function keep_above_threshold(Λ_of_t, wms_of_t, ε)
    # Might be a good idea to write a specialised function for generators
    # Anyway this has two passes on the weight vector, could be done in one pass
    wms_of_t_c = collect(wms_of_t)
    Λ_of_t_c = collect(Λ_of_t)
    Λ_of_t_kept = [Λ_of_t_c[i] for i in eachindex(Λ_of_t_c) if wms_of_t_c[i] >= ε]
    wms_of_t_kept = [wms_of_t_c[i] for i in eachindex(Λ_of_t_c) if wms_of_t_c[i] >= ε]
    return Λ_of_t_kept, wms_of_t_kept
end

function keep_fixed_number_of_weights(Λ_of_t, wms_of_t, k)
    if length(wms_of_t) <= k
        return Λ_of_t, wms_of_t
    else
        #This is intended to work also with iterators
        #Some more work might be needed, but it will only be implemented if this part of the code is critical for performance.
        wms_of_t_c = collect(wms_of_t)
        last_w = kmax(wms_of_t_c, k) |> minimum #smallest weight kept
        return keep_above_threshold(collect(Λ_of_t), wms_of_t_c, last_w)
    end
end

function keep_fixed_fraction(Λ_of_t, wms_of_t, fraction; total = 1.)
    #This is not guaranteed to keep exactly the correct fraction, should be called "keep at least fraction". This keeps more when there are several weights equal to the threshold
    sorted_wms_of_t = sort(wms_of_t, rev = true)
    cumulative_sum::Float64 = 0.
    i::Int64 = 0
    max_i = length(sorted_wms_of_t)
    # @show max_i
    frac_times_total::Float64 = fraction * total
    @inbounds while i <= max_i && cumulative_sum < frac_times_total
        i +=1
        # @show i
        cumulative_sum += sorted_wms_of_t[i]
    end
    threshold = sorted_wms_of_t[i]
    return keep_above_threshold(Λ_of_t, wms_of_t, threshold)
end

function keep_fixed_fraction_logw(Λ_of_t, logwms_of_t, fraction; logtotal = 0.)
    #This is not guaranteed to keep exactly the correct fraction, should be called "keep at least fraction". This keeps more when there are several weights equal to the threshold
    sorted_logwms_of_t = sort(logwms_of_t, rev = true)
    cumulative_logsum::Float64 = -Inf
    i::Int64 = 0
    max_i = length(sorted_logwms_of_t)
    # @show max_i
    logfrac_times_total::Float64 = log(fraction) + logtotal
    @inbounds while i <= max_i && cumulative_logsum < logfrac_times_total
        i +=1
        # @show i
        cumulative_logsum = logaddexp(cumulative_logsum, sorted_logwms_of_t[i])
    end
    logthreshold = sorted_logwms_of_t[i]
    return keep_above_threshold(Λ_of_t, logwms_of_t, logthreshold)
end

function prune_keeping_fixed_fraction_log_wms(Λ_of_t, log_wms_of_t, fraction::Real)
    Λ_of_t_kept, wms_of_t_kept = keep_fixed_fraction(Λ_of_t, exp.(log_wms_of_t), fraction)
    return Λ_of_t_kept, log.(normalise(wms_of_t_kept))
end

function prune_keeping_fixed_number_log_wms(Λ_of_t, log_wms_of_t, fixed_number::Integer)
    Λ_of_t_kept, log_wms_of_t_kept = keep_fixed_number_of_weights(Λ_of_t, log_wms_of_t, fixed_number)
    return Λ_of_t_kept, log_wms_of_t_kept .- logsumexp(log_wms_of_t_kept)
end

function prune_keeping_above_threshold_log_wms(Λ_of_t, log_wms_of_t, ε::Real)
    Λ_of_t_kept, log_wms_of_t_kept = keep_above_threshold(Λ_of_t, log_wms_of_t, log(ε))
    return Λ_of_t_kept, log_wms_of_t_kept .- logsumexp(log_wms_of_t_kept)
end

function prune_all_dicts(Λ_of_ts, wms_of_ts, pruning_function)
    Λ_of_ts_pruned = Dict()
    wms_of_ts_pruned = Dict()
    for k in keys(Λ_of_ts)
        Λ_of_ts_pruned[k], wms_of_ts_pruned[k] = pruning_function(Λ_of_ts[k], wms_of_ts[k])
    end
    return Λ_of_ts_pruned, wms_of_ts_pruned
end
function prune_all_dicts(Λ_of_ts::Dict{Float64, Array{Int64,1}}, wms_of_ts::Dict{Float64, Array{Float64,1}}, pruning_function::Function)
    Λ_of_ts_pruned = Dict{Float64, Array{Int64,1}}()
    wms_of_ts_pruned = Dict{Float64, Array{Float64,1}}()
    for k in keys(Λ_of_ts)
        Λ_of_ts_pruned[k], wms_of_ts_pruned[k] = pruning_function(Λ_of_ts[k], wms_of_ts[k])
    end
    return Λ_of_ts_pruned, wms_of_ts_pruned
end

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
################################################################################
#
#   debug.jl : debugging functions
#
################################################################################

# Internal debugging functions, not documented

function checkeq(P, M, m, K)
    i = 1
    R = parent(P)
    F = base_ring(P)
    defPol = K.definingPolynomial
    tmp2=R(1)
    for y in F
        if M[i, 1] == 1
            tmp2 *= P-y
        end
        i += 1
    end
    tmp2 %= defPol

    tmp1 = powmod(gcdinv(K.h1, defPol)[2], degree(P), defPol)
    tmp1 = mulmod(makeEquation(m, P, K.h0, K.h1), tmp1, defPol)

    return tmp1,tmp2
end

function checklog(L, defPol)
    l = length(L.factors)
    R = parent(L.factors[1])
    res = R(1)
    for i in 1:l
        tmp = powmod(L.factors[i], L.coefs[i], defPol)
        res = mulmod(res, tmp, defPol)
    end
    return res*L.unit
end

function checknum(num, M, P, h1, Q)
    res1 = Q(1)
    res2 = Q(1)
    N, det = rref(M)
    m, j = size(M)
    sol = fmpz[N[i, j] for i in 1:m]
    d = degree(P)

    piv = pivots(N, m-1)
    for i in 1:(m-1)
        loc = num[piv[i]]
        if sol[i] < 0
            exp = -BigInt(sol[i])
            res1 *= inv(Q(loc))^(exp)*Q(h1)^(d*exp)
            for f in factor(loc)
                res2 *= inv(Q(f[1]))^(f[2]*exp)
            end
            res2 *= Q(h1)^(d*exp)*inv(coeff(loc, degree(loc)))^exp
        else
            exp = BigInt(sol[i])
            for f in factor(loc)
                res2 *= Q(f[1])^(f[2]*exp)
            end
            res1 *= Q(loc)^exp*inv(Q(h1))^(d*exp)
            res2 *= inv(Q(h1))^(d*exp)*coeff(loc, degree(loc))^exp
        end
    end
    return res1, res2
end

function checkcol(M, j, P, defPol)
    i = 0
    F = base_ring(P)
    ζ = parent(P)(1)
    for y in F
        i += 1
        if M[i, j] == 1
            ζ *= P-y
        end
    end
    return ζ%defPol
end

function checkmat(M::MatElem, K, T::PolyElem)
    g = K.gen
    R = parent(g)
    F = base_ring(g)
    defPol = K.definingPolynomial
    r, c = size(M)
    for j in 1:c
        ζ = R(1)
        i = 0
        for y in F
            i += 1
            if M[i, j] != 0
                exp = Int(M[i, j])
                ζ *= powmod(T-y, exp, defPol)
            end
        end
        ζ *= K.h1
        ζ %= defPol
        ξ = powmod(g, M[r, j], defPol)
        if ζ != ξ
            break
        end
    end
end

function isConjugate(A, B)
    C = A*inv(B)
    q::Int = sqrt(length(parent(A[1, 1])))
    if C[1, 1] != 0
        C = inv(C[1, 1])*C
        if C[1, 2]^q == C[1, 2] && C[2, 1]^q == C[2, 1] && C[2, 2]^q == C[2, 2]
            return true
        end
    else
        C = inv(C[1, 2])*C
        if C[2, 1]^q == C[2, 1] && C[2, 2]^q == C[2, 2]
            return true
        end
    end
    return false
end

function truePgl(pgl)
    new = Array{typeof(pgl[1]), 1}()
    for M in pgl
        b = true
        for N in new
            if isConjugate(M, N)
                b = false
                break
            end
        end
        if b
            push!(new, M)
        end
    end
    return new
end

function isSplit(P::PolyElem)
    for f in factor(P)
        if degree(f[1]) != 1
            return false
        end
    end
    return true
end

function substitute(P::PolyElem, S::PolyElem)
    R = parent(S)()
    for i in 0:degree(P)
        R += coeff(P, i)*S^i
    end
    return R
end

function checkgkz(L, K, Q)
    l = length(L)
    R = parent(L[1][1])
    X = gen(R)
    res = one(R)
    q = length(base_ring(K.h0))
#    Q = K.h1*X^q-K.h0
    for j in 1:l
        poly = L[j][1]
        coef = L[j][2]
        if coef < 0
            poly = gcdinv(poly, Q)[2]
            coef = -coef
        end
        res = mulmod(res, poly^coef, Q)
    end
    return res
end

function checksplit(R, q)
    F = base_ring(R)
    X = gen(R)
    κ = 0
    l = length(F)
    ι = 0
    for B in F
        ι += 1
        P = X^(q+1)-B*X+B
        if isSplit(P)
            κ += 1
        end
        println((ι/l, κ))
    end
    return κ
end


text:
##==============================================================================
## AbstractDataEntry - Defined in src/entities/AbstractDFG.jl
##==============================================================================
# Fields to be implemented
# label
# id

getLabel(entry::AbstractDataEntry) = entry.label
getId(entry::AbstractDataEntry) = entry.id
getHash(entry::AbstractDataEntry) = hex2bytes(entry.hash)
getCreatedTimestamp(entry::AbstractDataEntry) = entry.createdTimestamp

##==============================================================================
## BlobStoreEntry
##==============================================================================
export BlobStoreEntry

"""
    $(TYPEDEF)
Genaral Data Store Entry.
"""
struct BlobStoreEntry <: AbstractDataEntry
    label::Symbol
    id::UUID
    blobstore::Symbol
    hash::String # Probably https://docs.julialang.org/en/v1/stdlib/SHA
    origin::String # E.g. user|robot|session|varlabel
    description::String
    mimeType::String
    createdTimestamp::ZonedDateTime # of when the entry was created
end

# TODO
"""
    $(TYPEDEF)
Data Entry in MongoDB.
"""
struct MongodbDataEntry <: AbstractDataEntry
    label::Symbol
    id::UUID
    oid::NTuple{12, UInt8} #mongodb object id - TODO Not needed with id::UUID unique, but perhaps usefull
    hash::String
    createdTimestamp::ZonedDateTime
    # mongodb
    # mongoConfig::MongoConfig
    #maybe other fields such as:
    #flags::Bool ready, valid, locked, permissions
    #MIMEType::String
end


text:
export PoolDims

"""
    PoolDims(x_size::NTuple{M}, k::Union{NTuple{L, Int}, Int};
            stride=k, padding=0, dilation=1)  where {M, L}

Dimensions for a "pooling" operation that can have an arbitrary input size, kernel size,
stride, dilation, and channel count.  Used to dispatch onto efficient implementations at
compile-time.
"""
struct PoolDims{N,K,S,P,D} <: ConvDims{N, S, P, D, false}
    I::NTuple{N,Int}
    C_in::Int
end

# Getters for both type parameters and fields
kernel_size(c::PoolDims{N,K,S,P,D}) where {N, K, S, P, D} = K
input_size(c::PoolDims) = c.I
channels_in(c::PoolDims) = c.C_in
channels_out(c::PoolDims) = c.C_in

# Convenience wrapper to create DenseConvDims objects
function PoolDims(x_size::NTuple{M}, k::Union{NTuple{L, Int}, Int};
                  stride=k, padding=0, dilation=1) where {M, L}
    # Expand `k` up to a tuple
    if typeof(k) <: Number
        k = ntuple(_ -> k, M - 2)
    end

    # Do common parameter validation
    stride, padding, dilation = check_spdf(x_size, (k..., 1, 1), stride, padding, dilation)

    # Build it
    return PoolDims{
        M - 2,
        k,
        stride,
        padding,
        dilation
    }(
        # Image spatial size
        x_size[1:end-2],

        # Input channels
        x_size[end-1],
    )
end

# Auto-take `size(x)` when `x` is an array.
function PoolDims(x::AbstractArray, k; kwargs...)
    return PoolDims(size(x), k; kwargs...)
end

# Useful for constructing a new PoolDims that has only a few elements different
# from the original progenitor object that it inherits shapes from.
function PoolDims(c::ConvDims; N=spatial_dims(c), I=input_size(c), K=kernel_size(c),
                       C_in=channels_in(c), S=stride(c), P=padding(c), D=dilation(c))
    return PoolDims{N, K, S, P, D}(I, C_in)
end

function check_dims(x::NTuple{M}, y::NTuple{M}, pdims::PoolDims) where {M}
    # First, check that channel counts are all correct:
    @assert x[end-1] == channels_in(pdims) DimensionMismatch("Data input channel count ($(x[end-1]) vs. $(channels_in(pdims)))")
    @assert y[end-1] == channels_out(pdims) DimensionMismatch("Data output channel count ($(y[end-1]) vs. $(channels_out(pdims)))")

    # Next, check that the spatial dimensions match up
    @assert x[1:end-2] == input_size(pdims) DimensionMismatch("Data input spatial size ($(x[1:end-2]) vs. $(input_size(pdims)))")
    @assert y[1:end-2] == output_size(pdims) DimensionMismatch("Data output spatial size ($(y[1:end-2]) vs. $(output_size(pdims)))")

    # Finally, check that the batch size matches
    @assert x[end] == y[end] DimensionMismatch("Batch size ($(x[end]) vs. $(y[end]))")
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# 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.
# *****************************************************************************

"""
GrammaticalEvolution
"""
module GE

export GEESParams, GEESResult, ge_search, exprsearch, SearchParams, SearchResult, get_derivtree

using Reexport
using ExprSearch
using RLESUtils, GitUtils, CPUTimeUtils, Observers, LogSystems
using GrammaticalEvolution
@reexport using LinearDerivTrees  #for pretty strings
using CPUTime
using RLESTypes.SymbolTable

import ..ExprSearch: exprsearch, get_fitness, get_derivtree, get_expr

include("logdefs.jl")

type GEESParams <: SearchParams
    #GrammaticalEvolution params
    genome_size::Int64
    pop_size::Int64
    maxwraps::Int64
    top_keep::Float64
    top_seed::Float64
    rand_frac::Float64
    prob_mutation::Float64
    mutation_rate::Float64
    default_code::Any
    max_iters::Int64
    logsys::LogSystem
    userargs::SymbolTable
end
GEESParams(genome_size::Int64, pop_size::Int64, maxwraps::Int64, top_keep::Float64,
    top_seed::Float64, rand_frac::Float64, prob_mutation::Float64, mutation_rate::Float64,
    default_code::Any, max_iters::Int64, logsys::LogSystem=logsystem();
    userargs::SymbolTable=SymbolTable()) =
        GEESParams(genome_size, pop_size, maxwraps, top_keep, top_seed, rand_frac,
        prob_mutation, mutation_rate, default_code, max_iters, logsys, userargs)

type GEESResult <: SearchResult
    tree::LinearDerivTree
    genome::Vector{Int64}
    fitness::Float64
    expr
    best_at_eval::Int64
    totalevals::Int64
end

exprsearch(p::GEESParams, problem::ExprProblem) = ge_search(p, problem::ExprProblem)

get_derivtree(result::GEESResult) = get_derivtree(result.tree)
get_expr(result::GEESResult) = result.expr
get_fitness(result::GEESResult) = result.fitness

function ge_search(p::GEESParams, problem::ExprProblem)
    @notify_observer(p.logsys.observer, "verbose1", ["Starting GE search"])
    @notify_observer(p.logsys.observer, "computeinfo", ["starttime", string(now())])

    initialize!(problem)
    grammar = get_grammar(problem)

    tree_params = LDTParams(grammar, p.genome_size)
    tree = LinearDerivTree(tree_params)

    tstart = CPUtime_start()

    #first iter is random init
    iter = 1
    pop = ExamplePopulation(p.pop_size, p.genome_size)
    evaluate!(grammar, pop, p, problem, tree)
    sort!(pop)
    fitness = pop[1].fitness
    code = pop[1].code
    nevals = iter * p.pop_size
    @notify_observer(p.logsys.observer, "elapsed_cpu_s", [nevals, CPUtime_elapsed_s(tstart)])
    @notify_observer(p.logsys.observer, "fitness", Any[iter, fitness])
    @notify_observer(p.logsys.observer, "fitness5", Any[iter, [pop[i].fitness for i=1:5]...])
    @notify_observer(p.logsys.observer, "code", Any[iter, string(code)])
    @notify_observer(p.logsys.observer, "population", Any[iter, pop])
    @notify_observer(p.logsys.observer, "current_best", [nevals, fitness, string(code)])
    iter += 1

    while iter <= p.max_iters
        pop = generate(grammar, pop, p.top_keep, p.top_seed, p.rand_frac, p.prob_mutation,
            p.mutation_rate, p, problem::ExprProblem, tree)
        fitness = pop[1].fitness #population is sorted, so first entry is the best
        code = pop[1].code
        nevals = iter * p.pop_size
        @notify_observer(p.logsys.observer, "elapsed_cpu_s", [nevals, CPUtime_elapsed_s(tstart)])
        @notify_observer(p.logsys.observer, "fitness", Any[iter, fitness])
        @notify_observer(p.logsys.observer, "fitness5", Any[iter, [pop[i].fitness for i=1:5]...])
        @notify_observer(p.logsys.observer, "code", Any[iter, string(code)])
        @notify_observer(p.logsys.observer, "population", Any[iter, pop])
        @notify_observer(p.logsys.observer, "current_best", [nevals, fitness, string(code)])
        iter += 1
    end
    @assert pop.best_ind.fitness == pop.best_fitness <= pop[1].fitness

    fitness = pop.best_fitness
    ind = pop.best_ind
    genome = ind.genome
    expr = ind.code
    best_at_eval = pop.best_at_eval
    totalevals = pop.totalevals

    play!(tree, ind)

    @assert expr == get_expr(tree) "expr=$expr, get_expr(tree)=$(get_expr(tree))"

    @notify_observer(p.logsys.observer, "result", [fitness, string(expr), best_at_eval, totalevals])

    #meta info
    @notify_observer(p.logsys.observer, "computeinfo", ["endtime",  string(now())])
    @notify_observer(p.logsys.observer, "computeinfo", ["hostname", gethostname()])
    @notify_observer(p.logsys.observer, "computeinfo", ["gitSHA",  get_SHA(dirname(@__FILE__))])
    @notify_observer(p.logsys.observer, "computeinfo", ["cpu_time", CPUtime_elapsed_s(tstart)])
    @notify_observer(p.logsys.observer, "parameters", ["genome_size", p.genome_size])
    @notify_observer(p.logsys.observer, "parameters", ["pop_size", p.pop_size])
    @notify_observer(p.logsys.observer, "parameters", ["maxwraps", p.maxwraps])
    @notify_observer(p.logsys.observer, "parameters", ["top_keep", p.top_keep])
    @notify_observer(p.logsys.observer, "parameters", ["top_seed", p.top_seed])
    @notify_observer(p.logsys.observer, "parameters", ["rand_frac", p.rand_frac])
    @notify_observer(p.logsys.observer, "parameters", ["prob_mutation", p.prob_mutation])
    @notify_observer(p.logsys.observer, "parameters", ["mutation_rate", p.mutation_rate])
    @notify_observer(p.logsys.observer, "parameters", ["default_code", string(p.default_code)])
    @notify_observer(p.logsys.observer, "parameters", ["max_iters", p.max_iters])

    GEESResult(tree, genome, fitness, expr, best_at_eval, totalevals)
end

function GrammaticalEvolution.evaluate!(grammar::Grammar, ind::ExampleIndividual,
    pop::ExamplePopulation, p::GEESParams, problem::ExprProblem, tree::LinearDerivTree)
    try
        #ind.code = transform(grammar, ind, maxwraps=p.maxwraps)
        play!(tree, ind)
        ind.code = get_expr(tree)
        ind.fitness = get_fitness(problem, get_derivtree(tree), p.userargs)
        pop.totalevals += 1
        if ind.fitness < pop.best_fitness
            pop.best_fitness = ind.fitness
            pop.best_ind = ind
            pop.best_at_eval = pop.totalevals
        end
    catch e
        #if !isa(e, MaxWrapException)
        if !isa(e, IncompleteException)
            rethrow(e)
        end
        ind.code = p.default_code
        ind.fitness = realmax(Float64)
    end
end

end #module

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
#############################################################################
#Copyright (c) 2016 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.
#
#############################################################################

"""
	function mlregressor(x::Array{Float64},y::Array{Float64},algorithm::AbstractString;X_test::Array{Float64}=[0.0],y_test::Array{Float64}=[0.0],test_sz=0.3,scaler="MinMaxScaler",rand_state=42,param_grid_kr = Dict("alpha"=> [1e1, 1e0, 0.5, 0.1, 5e-2, 1e-2, 5e-3, 1e-3],"gamma"=> logspace(-4, 4, 9)),param_grid_svm=Dict("C"=> [1e0, 2e0, 5e0, 1e1, 5e1, 1e2, 5e2, 1e3, 5e3, 1e4, 5e4, 1e5],"gamma"=> logspace(-4, 4, 9)),user_kernel="rbf")

Working on that, please be careful when using it.

INPUTS

	x: Array{Float64}, the spectra organised in rows (1 row = one spectrum). The spectra should share a common X axis;

	y::Array{Float64}, the target. Only a sigle target is possible for now.

	algorithm: AbstractString, "KernelRidge" or "SVM" or "LinearRegression"

OPTIONS

	X_test: Array{Float64}, spectra organised in rows (1 row = one spectrum) that you want to use as a testing dataset. THose spectra should not be present in the x (training) dataset. The spectra should share a common X axis;

	y_test: Array{Float64}, the target that you want to use as a testing dataset. Those targets should not be present in the y (training) dataset;

	scaler: String, the type of scaling performed. Choose between MinMaxScaler or StandardScaler, see http://scikit-learn.org/stable/modules/preprocessing.html for details. Default = "MinMaxScaler";

	rand_state: Float64, the random seed that is used for reproductibility of the results. Default = 42;

	param_grid_kr: Dictionary, containg the values of the hyperparameters that should be checked by gridsearch for the Kernel Ridge regression algorithm;

	param_grid_svm: Dictionary, containg the values of the hyperparameters that should be checked by gridsearch for the Support Vector regression algorithm.

For the last two parameters, the user is refered to the documentation of SciKit Learn. See the pages:

http://scikit-learn.org/stable/modules/kernel_ridge.html

http://scikit-learn.org/stable/modules/generated/sklearn.svm.SVR.html

OUTPUTS

	prediction_train: Array{Float64}, the predicted target values for the training y dataset;

	prediction_test: Array{Float64}, the predicted target values for the testing y_test dataset;

	model: A Scikit Learn object model, see the above link for details;

	X_scaler: A Scikit Learn scaler object for the x values;

	Y_scaler: A Scikit Learn scaler object for the y values;

NOTES

	For Support Vector and Kernel Ridge regressions, mlregressor performs a cross_validation search with using 5 KFold cross validators.

	If the results are poor with Support Vector and Kernel Ridge regressions, you will have to tune the param_grid_kr or param_grid_svm dictionnary that records the hyperparameter space to investigate during the cross validation.

"""
function mlregressor(x::Array{Float64},y::Array{Float64},algorithm::AbstractString;X_test::Array{Float64}=[0.0],y_test::Array{Float64}=[0.0],test_sz=0.3,scaling="yes",scaler="MinMaxScaler",rand_state=42,param_grid_kr = Dict("alpha"=> [1e1, 1e0, 0.5, 0.1, 5e-2, 1e-2, 5e-3, 1e-3],"gamma"=> logspace(-4, 4, 9)),param_grid_svm=Dict("C"=> [1e0, 2e0, 5e0, 1e1, 5e1, 1e2, 5e2, 1e3, 5e3, 1e4, 5e4, 1e5],"gamma"=> logspace(-4, 4, 9)),user_kernel="rbf")

	if size(X_test,1) == 1
		X_train, X_test, y_train, y_test = model_selection[:train_test_split](x, reshape(y,size(y,1),1), test_size=test_sz, random_state=rand_state)
	elseif size(X_test,2) == size(x,2)
		X_train = x[:,:]
		y_train = y[:,:]
	else
		error("You tried to provide a testing dataset that has a different number of features (in columns) than the training set. Please correct this.")
	end

	# to avoid any problem with SciKit Learn quite annoying demands for the shape of arrays...
	y_train =reshape(y_train,size(y_train,1),1)
	y_test = reshape(y_test,size(y_test,1),1)

	# initialising the preprocessor scaler
	if scaler == "StandardScaler"
		X_scaler = preprocessing[:StandardScaler]()
		Y_scaler = preprocessing[:StandardScaler]()
	elseif scaler == "MinMaxScaler"
		X_scaler = preprocessing[:MinMaxScaler]()
		Y_scaler = preprocessing[:MinMaxScaler]()
	else
		error("Choose the scaler between MinMaxScaler and StandardScaler")
	end

	X_scaler[:fit](X_train)
	Y_scaler[:fit](y_train)

	#scaling the data
	X_train_sc = X_scaler[:transform](X_train)
	y_train_sc = Y_scaler[:transform](y_train)

	X_test_sc = X_scaler[:transform](X_test)
	y_test_sc = Y_scaler[:transform](y_test)

	if algorithm=="KernelRidge"
		clf_kr = kernel_ridge[:KernelRidge](kernel=user_kernel, gamma=0.1)
		model = model_selection[:GridSearchCV](clf_kr,cv=5,param_grid=param_grid_kr)
	elseif algorithm=="SVM"
		clf_svm = svm[:SVR](kernel=user_kernel, gamma=0.1)
		model = model_selection[:GridSearchCV](clf_svm,cv=5,param_grid=param_grid_svm)
	elseif algorithm=="Lasso"
		clf_lasso = linear_model[:Lasso](alpha=0.1)
		model = model_selection[:GridSearchCV](clf_lasso,cv=5,param_grid=Dict("alpha"=> [1e-3, 1e-2, 1e-1, 1., 1e1, 1e2, 1e3, 1e4]))
	elseif algorithm=="ElasticNet"
		clf_ElasticNet = linear_model[:ElasticNet](alpha=0.1,l1_ratio=0.5)
		model = model_selection[:GridSearchCV](clf_ElasticNet,cv=5,param_grid=Dict("alpha"=> [1e-3, 1e-2, 1e-1, 1., 1e1, 1e2, 1e3, 1e4]))
	elseif algorithm=="LinearRegression"
		model = linear_model[:LinearRegression]()
	end

	if scaling == "yes"
		model[:fit](X_train_sc, vec(y_train_sc))
		predict_train_sc = model[:predict](X_train_sc)
		prediction_train = Y_scaler[:inverse_transform](predict_train_sc)
		predict_test_sc = model[:predict](X_test_sc)
		prediction_test = Y_scaler[:inverse_transform](predict_test_sc)
		return prediction_train, prediction_test, model, X_scaler, Y_scaler
	else
		model[:fit](X_train, vec(y_train))
		prediction_train = model[:predict](X_train)
		prediction_test = model[:predict](X_test)
		return prediction_train, prediction_test, model
	end
end


text:
"""
**Template Likelihood**

Template file for likelihood creation

```julia
TemplateLikelihood()
```
See all functions you need to implement
---

"""
struct TemplateLikelihood{T<:Real} <: Likelihood{T}
    θ::LatentArray{AbstractVector{T}}
    function TemplateLikelihood{T}() where {T<:Real}
        new{T}()
    end
    function TemplateLikelihood{T}(θ::AbstractVector{<:AbstractVector{<:Real}}) where {T<:Real}
        new{T}(θ)
    end
end

function TemplateLikelihood()
    TemplateLikelihood{Float64}()
end

function init_likelihood(likelihood::TemplateLikelihood{T},inference::Inference{T},nLatent::Integer,nSamplesUsed::Integer) where T
    if inference isa AnalyticVI || inference isa GibbsSampling
        TemplateLikelihood{T}([zeros(T,nSamplesUsed) for _ in 1:nLatent])
    else
        TemplateLikelihood{T}()
    end
end

function pdf(l::TemplateLikelihood,y::Real,f::Real)

end

function Base.show(io::IO,model::TemplateLikelihood{T}) where T
    print(io,"Template Likelihood")
end

function compute_proba(l::TemplateLikelihood{T},μ::AbstractVector{T},σ²::AbstractVector{T}) where {T<:Real}
    N = length(μ)
    pred = zeros(T,N)
    for i in 1:N
        pred[i]  = 0.0
    end
    return pred
end

### Local Updates Section ###

function local_updates!(model::VGP{<:TemplateLikelihood,<:AnalyticVI})
end

function local_updates!(model::SVGP{<:LogisticLikelihood,<:AnalyticVI})
end

function sample_local!(model::VGP{<:LogisticLikelihood,<:GibbsSampling})
    return nothing
end

### Natural Gradient Section ###

function expec_μ(model::VGP{<:LogisticLikelihood,<:AnalyticVI},index::Integer)
end

function ∇μ(model::VGP{<:LogisticLikelihood})
end

function expec_μ(model::SVGP{<:LogisticLikelihood,<:AnalyticVI},index::Integer)
end

function ∇μ(model::SVGP{<:LogisticLikelihood})
end

function expec_Σ(model::AbstractGP{<:LogisticLikelihood,<:AnalyticVI},index::Integer)
    return model.likelihood.θ[index]
end

function ∇Σ(model::AbstractGP{<:LogisticLikelihood})
    return model.likelihood.θ
end

### ELBO Section ###

function ELBO(model::AbstractGP{<:LogisticLikelihood,<:AnalyticVI})
    return expecLogLikelihood(model) - GaussianKL(model) - PolyaGammaKL(model)
end

function expecLogLikelihood(model::VGP{<:LogisticLikelihood,<:AnalyticVI})
    tot = -model.nLatent*(0.5*model.nSample*logtwo)
    tot += sum(broadcast((μ,y,θ,Σ)->0.5.*(sum(μ.*y)-dot(θ,Σ+abs2.(μ))),
                        model.μ,model.y,model.likelihood.θ,diag.(model.Σ)))
    return tot
end

function expecLogLikelihood(model::SVGP{<:LogisticLikelihood,<:AnalyticVI})
    tot = -model.nLatent*(0.5*model.inference.nSamplesUsed*logtwo)
    tot += sum(broadcast((κμ,y,θ,κΣκ,K̃)->0.5.*(sum(κμ.*y[model.inference.MBIndices])-dot(θ,K̃+κΣκ+abs2.(κμ))),
                        model.κ.*model.μ,model.y,model.likelihood.θ,opt_diag.(model.κ.*model.Σ,model.κ),model.K̃))
    return model.inference.ρ*tot
end

function PolyaGammaKL(model::AbstractGP{<:LogisticLikelihood})
    model.inference.ρ*sum(broadcast(PolyaGammaKL,[ones(length(model.likelihood.c[1]))],model.likelihood.c,model.likelihood.θ))
end

### Gradient Section ###

function gradpdf(::LogisticLikelihood,y::Int,f::T) where {T<:Real}
    σ=y*f
    σ*(one(T)-σ)
end

function hessiandiagpdf(::LogisticLikelihood,y::Int,f::T) where {T<:Real}
    σ=y*f
    σ*(one(T)-2σ + abs2(σ))
end


text:
type StaticLaneFollowingDriver <: LaneFollowingDriver
    a::LaneFollowingAccel
    StaticLaneFollowingDriver(a::LaneFollowingAccel=LaneFollowingAccel(0.0)) = new(a)
    StaticLaneFollowingDriver(a::Float64) = new(LaneFollowingAccel(a))
end
get_name(::StaticLaneFollowingDriver) = "ProportionalSpeedTracker"
Base.rand(model::StaticLaneFollowingDriver) = model.a
Distributions.pdf(model::StaticLaneFollowingDriver, a::LaneFollowingAccel) = isapprox(a.a, model.a.a) ? Inf : 0.0
Distributions.logpdf(model::StaticLaneFollowingDriver, a::LaneFollowingAccel) = isapprox(a.a, model.a.a) ? Inf : -Inf


text:
# include("D:\\Programs\\TSSOS\\src\\TSSOS.jl")
using TSSOS
# You need to add the package PowerModels.
include("D:\\Programs\\TSSOS\\example\\modelopf.jl")

# You need to download the problem data from PGLiB (https://github.com/power-grid-lib/pglib-opf) first.
cd("D:\\Programs\\PolyOPF\\pglib")

silence()

case = "pglib_opf_case3_lmbd"
AC = 5812.6
opfdata = parse_file(case * ".m")

# the first order relaxation
model,_ = pop_opf_real(opfdata, normal=true, AngleCons=true, LineLimit="relax")
n = model.n
m = model.m
numeq = model.numeq
supp = model.supp
coe = model.coe
mc = maximum(abs.(coe[1]))
coe[1]=coe[1]./mc

time = @elapsed begin
opt,sol,popd = cs_tssos_first(supp, coe, n, 1, numeq=numeq, CS=false, TS="MF", MomentOne=false)
end
opt *= mc
mb = maximum(maximum.(popd.sb)) # maximal block size
gap = (AC-opt)*100/AC # optimality gap
println("n = $n, m = $m")
println("opt = $opt, time = $time, mb = $mb, gap = $gap%")

# the minimum order relaxation
model,_ = pop_opf_real(opfdata, normal=true, AngleCons=true, LineLimit=true)
n = model.n
m = model.m
numeq = model.numeq
supp = model.supp
coe = model.coe
mc = maximum(abs.(coe[1]))
coe[1] = coe[1]./mc

time = @elapsed begin
opt,sol,popd = cs_tssos_first(supp, coe, n, "min", numeq=numeq, CS="MF", TS="block", MomentOne=true)
end
opt *= mc
maxc = maximum(popd.cliquesize) # maximal clique size
mb = maximum(maximum.(popd.sb)) # maximal block size
gap = 100*(AC-opt)/AC # optimality gap
println("n = $n, m = $m")
println("mc = $maxc, opt = $opt, time = $time, mb = $mb, gap = $gap%")


text:
##
cd(@__DIR__)
⊕(x, y) = x * y
⊗(x, y) = x + y
function parse_line(line)
    Meta.parse(replace(replace(line, r"\*" => s"⊕"), r"\+" => s"⊗"))
end
input = parse_line.(eachline("input.txt"))

##
sum(eval.(input))


text:
module Atmos

export AtmosModel,
       AtmosAcousticLinearModel, AtmosAcousticGravityLinearModel,
       RemainderModel

using LinearAlgebra, StaticArrays
using ..VariableTemplates
using ..MoistThermodynamics
using ..PlanetParameters
import ..MoistThermodynamics: internal_energy
using ..SubgridScaleParameters
using GPUifyLoops
using ..MPIStateArrays: MPIStateArray

import CLIMA.DGmethods: BalanceLaw, vars_aux, vars_state, vars_gradient,
                        vars_diffusive, vars_integrals, flux_nondiffusive!,
                        flux_diffusive!, source!, wavespeed, boundary_state!,
                        gradvariables!, diffusive!, init_aux!, init_state!,
                        update_aux!, integrate_aux!, LocalGeometry, lengthscale,
                        resolutionmetric, DGModel, num_integrals,
                        nodal_update_aux!, indefinite_stack_integral!,
                        reverse_indefinite_stack_integral!
using ..DGmethods.NumericalFluxes

"""
    AtmosModel <: BalanceLaw

A `BalanceLaw` for atmosphere modeling.

# Usage

    AtmosModel(orientation, ref_state, turbulence, moisture, radiation, source,
               boundarycondition, init_state)

"""
struct AtmosModel{O,RS,T,M,R,SU,S,BC,IS} <: BalanceLaw
  orientation::O
  ref_state::RS
  turbulence::T
  moisture::M
  radiation::R
  subsidence::SU
  source::S
  # TODO: Probably want to have different bc for state and diffusion...
  boundarycondition::BC
  init_state::IS
end

function vars_state(m::AtmosModel, FT)
  @vars begin
    ρ::FT
    ρu::SVector{3,FT}
    ρe::FT
    turbulence::vars_state(m.turbulence, FT)
    moisture::vars_state(m.moisture, FT)
    radiation::vars_state(m.radiation, FT)
  end
end
function vars_gradient(m::AtmosModel, FT)
  @vars begin
    u::SVector{3,FT}
    h_tot::FT
    turbulence::vars_gradient(m.turbulence,FT)
    moisture::vars_gradient(m.moisture,FT)
  end
end
function vars_diffusive(m::AtmosModel, FT)
  @vars begin
    ρτ::SHermitianCompact{3,FT,6}
    ρd_h_tot::SVector{3,FT}
    turbulence::vars_diffusive(m.turbulence,FT)
    moisture::vars_diffusive(m.moisture,FT)
  end
end

function vars_aux(m::AtmosModel, FT)
  @vars begin
    ∫dz::vars_integrals(m, FT)
    ∫dnz::vars_integrals(m, FT)
    coord::SVector{3,FT}
    orientation::vars_aux(m.orientation, FT)
    ref_state::vars_aux(m.ref_state,FT)
    turbulence::vars_aux(m.turbulence,FT)
    moisture::vars_aux(m.moisture,FT)
    radiation::vars_aux(m.radiation,FT)
  end
end
function vars_integrals(m::AtmosModel,FT)
  @vars begin
    radiation::vars_integrals(m.radiation,FT)
  end
end

include("orientation.jl")
include("ref_state.jl")
include("turbulence.jl")
include("moisture.jl")
include("subsidence.jl")
include("radiation.jl")
include("source.jl")
include("boundaryconditions.jl")
include("linear.jl")
include("remainder.jl")

"""
    flux_nondiffusive!(m::AtmosModel, flux::Grad, state::Vars, aux::Vars,
                       t::Real)

Computes flux non-diffusive flux portion of `F` in:

```
∂Y
-- = - ∇ • (F_{adv} + F_{press} + F_{nondiff} + F_{diff}) + S(Y)
∂t
```
Where

 - `F_{adv}`      Advective flux             ; see [`flux_advective!`]@ref()
 - `F_{press}`    Pressure flux              ; see [`flux_pressure!`]@ref()
 - `F_{diff}`     Fluxes that state gradients; see [`flux_diffusive!`]@ref()
"""
@inline function flux_nondiffusive!(m::AtmosModel, flux::Grad, state::Vars,
                                    aux::Vars, t::Real)
  ρ = state.ρ
  ρinv = 1/ρ
  ρu = state.ρu
  u = ρinv * ρu
  z  = altitude(m.orientation, aux)

  # advective terms
  usub = subsidence_velocity(m.subsidence, z)
  ẑ = vertical_unit_vector(m.orientation, aux)
  u_tot = u .- usub * ẑ
  flux.ρ   = ρ * u_tot
  flux.ρu  = ρ * u_tot .* u_tot'
  flux.ρe  = u_tot * state.ρe

  # pressure terms
  p = pressure(m.moisture, m.orientation, state, aux)

  if m.ref_state isa HydrostaticState
    flux.ρu += (p-aux.ref_state.p)*I
  else
    flux.ρu += p*I
  end
  flux.ρe += u*p
  flux_radiation!(m.radiation, m, flux, state, aux, t)
  flux_moisture!(m.moisture, m, flux, state, aux, t)
end

@inline function flux_diffusive!(m::AtmosModel, flux::Grad, state::Vars,
                                 diffusive::Vars, aux::Vars, t::Real)
  ρinv = 1/state.ρ
  u = ρinv * state.ρu

  # diffusive
  ρτ = diffusive.ρτ
  ρd_h_tot = diffusive.ρd_h_tot
  flux.ρu += ρτ
  flux.ρe += ρτ*u
  flux.ρe += ρd_h_tot
  flux_diffusive!(m.moisture, flux, state, diffusive, aux, t)
end

@inline function wavespeed(m::AtmosModel, nM, state::Vars, aux::Vars, t::Real)
  ρinv = 1/state.ρ
  u = ρinv * state.ρu
  return abs(dot(nM, u)) + soundspeed(m.moisture, m.orientation, state, aux)
end

function gradvariables!(atmos::AtmosModel, transform::Vars, state::Vars, aux::Vars, t::Real)
  ρinv = 1/state.ρ
  transform.u = ρinv * state.ρu
  transform.h_tot = total_specific_enthalpy(atmos.moisture, atmos.orientation, state, aux)

  gradvariables!(atmos.moisture, transform, state, aux, t)
  gradvariables!(atmos.turbulence, transform, state, aux, t)
end

function symmetrize(X::StaticArray{Tuple{3,3}})
  SHermitianCompact(SVector(X[1,1], (X[2,1] + X[1,2])/2, (X[3,1] + X[1,3])/2, X[2,2], (X[3,2] + X[2,3])/2, X[3,3]))
end

function diffusive!(m::AtmosModel, diffusive::Vars, ∇transform::Grad, state::Vars, aux::Vars, t::Real)
  ∇u = ∇transform.u
  # strain rate tensor
  S = symmetrize(∇u)
  # kinematic viscosity tensor
  ρν = dynamic_viscosity_tensor(m.turbulence, S, state, diffusive, ∇transform, aux, t)
  # momentum flux tensor
  diffusive.ρτ = scaled_momentum_flux_tensor(m.turbulence, ρν, S)

  ∇h_tot = ∇transform.h_tot
  # turbulent Prandtl number
  diag_ρν = ρν isa Real ? ρν : diag(ρν) # either a scalar or matrix
  # Diffusivity ρD_t = ρν/Prandtl_turb
  ρD_t = diag_ρν * inv_Pr_turb
  # diffusive flux of total energy
  diffusive.ρd_h_tot = -ρD_t .* ∇transform.h_tot

  # diffusivity of moisture components
  diffusive!(m.moisture, diffusive, ∇transform, state, aux, t, ρD_t)
  # diffusion terms required for SGS turbulence computations
  diffusive!(m.turbulence, diffusive, ∇transform, state, aux, t, ρD_t)
end

function update_aux!(dg::DGModel, m::AtmosModel, Q::MPIStateArray, t::Real)
  FT = eltype(Q)
  auxstate = dg.auxstate

  if num_integrals(m, FT) > 0
    indefinite_stack_integral!(dg, m, Q, auxstate, t)
    reverse_indefinite_stack_integral!(dg, m, auxstate, t)
  end

  nodal_update_aux!(atmos_nodal_update_aux!, dg, m, Q, t)

  return true
end

function atmos_nodal_update_aux!(m::AtmosModel, state::Vars, aux::Vars,
                                 t::Real)
  atmos_nodal_update_aux!(m.moisture, m, state, aux, t)
  atmos_nodal_update_aux!(m.radiation, m, state, aux, t)
  atmos_nodal_update_aux!(m.turbulence, m, state, aux, t)
end

function integrate_aux!(m::AtmosModel, integ::Vars, state::Vars, aux::Vars)
  integrate_aux!(m.radiation, integ, state, aux)
end

# TODO: figure out a nice way to handle this
function init_aux!(m::AtmosModel, aux::Vars, geom::LocalGeometry)
  aux.coord = geom.coord
  atmos_init_aux!(m.orientation, m, aux, geom)
  atmos_init_aux!(m.ref_state, m, aux, geom)
  atmos_init_aux!(m.turbulence, m, aux, geom)
end

"""
    source!(m::AtmosModel, source::Vars, state::Vars, aux::Vars, t::Real)
Computes flux `S(Y)` in:
```
∂Y
-- = - ∇ • F + S(Y)
∂t
```
"""
function source!(m::AtmosModel, source::Vars, state::Vars, aux::Vars, t::Real)
  atmos_source!(m.source, m, source, state, aux, t)
end

boundary_state!(nf, m::AtmosModel, x...) =
  atmos_boundary_state!(nf, m.boundarycondition, m, x...)

# FIXME: This is probably not right....
boundary_state!(::CentralGradPenalty, bl::AtmosModel, _...) = nothing

function init_state!(m::AtmosModel, state::Vars, aux::Vars, coords, t, args...)
  m.init_state(state, aux, coords, t, args...)
end

end # module

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
include("../utils.jl")
include("group_utils.jl")

@testset "Special Orthogonal group" begin
    G = SpecialOrthogonal(3)
    @test repr(G) == "SpecialOrthogonal(3)"
    M = base_manifold(G)
    @test M === Rotations(3)
    x = Matrix(I, 3, 3)
    @test is_default_metric(MetricManifold(G, EuclideanMetric()))

    types = [Matrix{Float64}]
    ω = [[1.0, 2.0, 3.0], [3.0, 2.0, 1.0], [1.0, 3.0, 2.0]]
    pts = [exp(M, x, hat(M, x, ωi)) for ωi in ω]
    vpts = [hat(M, x, [-1.0, 2.0, 0.5]), hat(M, x, [1.0, 0.0, 0.5])]

    ge = allocate(pts[1])
    identity_element!(G, ge)
    @test isapprox(ge, I; atol=1e-10)

    gI = Identity(G)
    gT = allocate_result(G, exp, gI, log(G, pts[1], pts[2]))
    @test size(gT) == size(ge)
    @test eltype(gT) == eltype(ge)
    gT = allocate_result(G, log, gI, pts[1])
    @test size(gT) == size(ge)
    @test eltype(gT) == eltype(ge)

    for T in types
        gpts = convert.(T, pts)
        vgpts = convert.(T, vpts)
        @test compose(G, gpts[1], gpts[2]) ≈ gpts[1] * gpts[2]
        @test translate_diff(G, gpts[2], gpts[1], vgpts[1], LeftAction()) ≈ vgpts[1]
        @test translate_diff(G, gpts[2], gpts[1], vgpts[1], RightAction()) ≈
              transpose(gpts[2]) * vgpts[1] * gpts[2]
        test_group(
            G,
            gpts,
            vgpts,
            vgpts;
            test_diff=true,
            test_invariance=true,
            test_lie_bracket=true,
            test_adjoint_action=true,
        )
    end

    @testset "Decorator forwards to group" begin
        DM = NotImplementedGroupDecorator(G)
        test_group(DM, pts, vpts, vpts; test_diff=true)
    end

    @testset "Group forwards to decorated" begin
        retraction_methods = [
            Manifolds.PolarRetraction(),
            Manifolds.QRRetraction(),
            Manifolds.GroupExponentialRetraction(LeftAction()),
            Manifolds.GroupExponentialRetraction(RightAction()),
        ]

        inverse_retraction_methods = [
            Manifolds.PolarInverseRetraction(),
            Manifolds.QRInverseRetraction(),
            Manifolds.GroupLogarithmicInverseRetraction(LeftAction()),
            Manifolds.GroupLogarithmicInverseRetraction(RightAction()),
        ]

        test_manifold(
            G,
            pts;
            test_injectivity_radius=false,
            test_project_tangent=true,
            test_musical_isomorphisms=false,
            retraction_methods=retraction_methods,
            inverse_retraction_methods=inverse_retraction_methods,
            exp_log_atol_multiplier=20,
            retraction_atol_multiplier=12,
            is_tangent_atol_multiplier=1.2,
            test_atlases=(Manifolds.RetractionAtlas(),),
        )

        @test injectivity_radius(G) == injectivity_radius(M)
        @test injectivity_radius(G, pts[1]) == injectivity_radius(M, pts[1])
        @test injectivity_radius(G, pts[1], PolarRetraction()) ==
              injectivity_radius(M, pts[1], PolarRetraction())

        y = allocate(pts[1])
        exp!(G, y, pts[1], vpts[1])
        @test isapprox(M, y, exp(M, pts[1], vpts[1]))

        y = allocate(pts[1])
        retract!(G, y, pts[1], vpts[1])
        @test isapprox(M, y, retract(M, pts[1], vpts[1]))

        w = allocate(vpts[1])
        inverse_retract!(G, w, pts[1], pts[2])
        @test isapprox(M, pts[1], w, inverse_retract(M, pts[1], pts[2]))

        w = allocate(vpts[1])
        inverse_retract!(G, w, pts[1], pts[2], QRInverseRetraction())
        @test isapprox(
            M,
            pts[1],
            w,
            inverse_retract(M, pts[1], pts[2], QRInverseRetraction()),
        )

        z = collect(reshape(1.0:9.0, 3, 3))
        @test isapprox(M, project(G, z), project(M, z))
        z2 = allocate(pts[1])
        project!(G, z2, z)
        @test isapprox(M, z2, project(M, z))
    end

    @testset "vee/hat" begin
        X = vpts[1]
        pe = Identity(G)

        Xⁱ = vee(G, identity_element(G), X)
        @test Xⁱ ≈ vee(G, pe, X)

        X2 = hat(G, pts[1], Xⁱ)
        @test isapprox(M, pe, X2, hat(G, pe, Xⁱ); atol=1e-6)
    end
    @testset "Identity and get_vector/get_coordinates" begin
        e = Identity(G)

        eF = Identity(SpecialEuclidean(3))
        c = [1.0, 0.0, 0.0]
        Y = zeros(representation_size(G))
        get_vector_lie!(G, Y, c, Manifolds.VeeOrthogonalBasis())
        @test Y ≈ get_vector_lie(G, c, Manifolds.VeeOrthogonalBasis())
        get_vector!(M, Y, identity_element(G), c, Manifolds.VeeOrthogonalBasis())
        @test Y ≈ get_vector(M, identity_element(G), c, Manifolds.VeeOrthogonalBasis())

        @test get_coordinates(M, identity_element(G), Y, Manifolds.VeeOrthogonalBasis()) ==
              c
        c2 = similar(c)
        get_coordinates_lie!(G, c2, Y, Manifolds.VeeOrthogonalBasis())
        @test c == c2
        get_coordinates!(M, c2, identity_element(G), Y, Manifolds.VeeOrthogonalBasis())
        @test c == c2

        q = zeros(3, 3)
        mul!(q, e, e)
        @test isone(q)
        e2 = Identity(G)
        @test mul!(e2, e, e) === e2
    end
end

max_stars_repo_name: JuliaNLSolvers/ManifoldMuseum

max_issues_repo_name: JuliaNLSolvers/ManifoldMuseum

max_forks_repo_name: JuliaNLSolvers/ManifoldMuseum


text:
using MuZero
import Zygote # @todo figure out how to remove

# User defined state and game implementations
include("state.jl")
include("game.jl")

# Configure Learner
learner = Learner(TicTacToe)

# Let's make sure our policy test cases look reasonable.
learner |> tests |> svg

# The predictions from this untrained model should be random
# and mostly incorrect.
learner |> test_results |> svg

#learner |> policy_tests |> first |> svg
#pt = learner |> policy_tests |> first
#policy_test_result = PolicyTestResult(pt, Flux.softmax(collect(rand(9))))

# Self-play 10 sample games to ensure proper game implementation.
selfplay!(learner; limit=10)

# Show a random self-played game
learner |> randgame |> svg

# Get a random training sample
learner |> randtrainingsample |> svg

# Get a random batch of training samples
learner |> randtrainingbatch |> show

# If game play seems to be working properly, then start indefinite self play on multiple threads.
using Base.Threads: @spawn
for i = 1:num_actors(learner)
    @spawn selfplay!(learner)
end

# Let's make sure threading is working propery and the games are being recorded
length(learner.replay_buffer)

# Train three steps to make sure there are no errors. This 1st step
# may take a several minutes if your machine has no usable
# Nvidia GPU.
@time fit!(learner, steps=3)

# Let's ensure the initial loss was stored, which means the model is training
learner |> losses |> svg

# Let's see how close our predictions are to their targets
learner |> rand_training_result |> svg
#learner |> training_results |> svg

# Visualize test set errors
learner |> test_results |> svg

# Train network for a few more steps
fit!(learner, steps=20)

# We should see steady decreases in the loss
learner |> losses |> svg

# We should see steady increase in the accuracy
learner |> accuracy |> svg

# Let's train indefinitely on a separate thread and periodically check the
# progress.
res = @spawn fit!(learner)

# function losses_svg()
#     function Base.show(io::IO, m::MIME"image/svg+xml", g::GridSVG)
#         draw(g)
#         print(io, read(output_file(g), String))
#     end
# end

# # Visualize any game from batch using the encoded state key
# df = get_df(batch)
# state = TicTacToeState(df."State Key"[1])
# state |> show

# # Or, use the following convenience function to get the state
# # @todo Stopped here!
# get_state(batch, 1) |> show

# Learning rate finder
# @todo implement

# Open up "losses.png" and "accuracy.png" in separate windows
# to view training loss and testing accuracy graphs in real-time.

# Alternatively, print the "losses" and/or "accuracy" arrays
# in the REPL to view progress.
losses
accuracy

# If losses is empty and/or not increasing, see status.txt
# to check if there was an error during training.
`cat ./status.txt` |> run

# Also, you can check "fetch(result)" to see if there was an exception.
# If no exception prints (and the call to fetch hangs), then the training is still
# starting, and you can interrupt the call to fetch
fetch(result)

# Export Flux model
MuZero.export(learner)


text:
module Quadrature

typealias Rule1d Tuple{Vector{Float64}, Vector{Float64}}

include("gauss_quad.jl")
include("quadgrid.jl")

export Rule1d
end


text:
using Terra

f = "/mnt/z/ChinaHW/Observed/clusterId_ncell=16/clusterId_HI_Tmax_(1961-2016).nc"

@time r = brick(f);


text:
include("1.2.2.1 (a+b x^2+c x^4)^p.jl")
include("1.2.2.2 (d x)^m (a+b x^2+c x^4)^p.jl")
include("1.2.2.3 (d+e x^2)^q (a+b x^2+c x^4)^p.jl")
include("1.2.2.4 (f x)^m (d+e x^2)^q (a+b x^2+c x^4)^p.jl")
include("1.2.2.5 P(x) (a+b x^2+c x^4)^p.jl")
include("1.2.2.6 P(x) (d x)^m (a+b x^2+c x^4)^p.jl")
include("1.2.2.7 P(x) (d+e x^2)^q (a+b x^2+c x^4)^p.jl")
include("1.2.2.8 P(x) (d+e x)^q (a+b x^2+c x^4)^p.jl")

max_stars_repo_path: src/rules/1 Algebraic functions/1.2 Trinomial products/1.2.2 Quartic/.jl


text:
function getSubParams(Hparam::RegularMesh, M::RegularMesh,i::Array{Int64},NumCells::Array{Int64},overlap::Array{Int64})
	subMesh   = getSubMeshOfCell(NumCells,overlap,i,Hparam);
	return subMesh;
end

function getLap(M::RegularMesh)
	G = getNodalGradientMatrix(M);
	Ar = G'*G;
	Ar = Ar + 1e-5*opnorm(Ar,1)*sparse(1.0I,size(Ar,2),size(Ar,2));
	return Ar;
end

function getDirichletMassNodalMesh(DDparam::DomainDecompositionParam,problem_param::RegularMesh,i::Array{Int64})
	d = getDirichletMassNodal(DDparam.numDomains,DDparam.overlap,i,DDparam.Mesh.n);
	d.*=(0.1*4.0)/prod(DDparam.Mesh.h);
	return d;
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

include("model_Access.jl")
include("model_AuthenticationMethod.jl")
include("model_AzureAsyncOperationResult.jl")
include("model_CloudError.jl")
include("model_CloudErrorBody.jl")
include("model_Error.jl")
include("model_ErrorDetails.jl")
include("model_IPAllocationMethod.jl")
include("model_IPVersion.jl")
include("model_ManagedServiceIdentity.jl")
include("model_ManagedServiceIdentityUserAssignedIdentities.jl")
include("model_ProvisioningState.jl")
include("model_Resource.jl")
include("model_SubResource.jl")
include("model_TagsObject.jl")


text:
function ftriple!(f::Vector{Float64},x::Vector{Float64},beta::Vector{Float64},A0::Vector{Float64},
    rho::Float64,c0::Vector{Float64},W::Vector{Float64})
    # f = zeros(8);
    # conservation of mass
    f[1] = x[1] .- x[3] .- x[5] .- x[7];

    # total P
    f[2] = (beta[1].*(sqrt.(x[2]) .- sqrt.(A0[1])) + 0.5.*rho.*(x[1]/x[2]).^2) .-
        (beta[2].*(sqrt.(x[4]) .- sqrt.(A0[2])) + 0.5.*rho.*(x[3]/x[4]).^2);
    f[3] = (beta[1].*(sqrt.(x[2]) .- sqrt.(A0[1])) + 0.5.*rho.*(x[1]/x[2]).^2) .-
        (beta[3].*(sqrt.(x[6]) .- sqrt.(A0[3])) + 0.5.*rho.*(x[5]/x[6]).^2);
    f[4] = (beta[1].*(sqrt.(x[2]) .- sqrt.(A0[1])) + 0.5.*rho.*(x[1]/x[2]).^2) .-
        (beta[4].*(sqrt.(x[8]) .- sqrt.(A0[4])) + 0.5.*rho.*(x[7]/x[8]).^2);

    # Riemann invariants
    f[5] = (x[1]./x[2]) .+ 4.*(sqrt.(0.5.*beta[1]./rho).*x[2].^0.25 .- c0[1]) .- W[1];
    f[6] = (x[3]./x[4]) .- 4.*(sqrt.(0.5.*beta[2]./rho).*x[4].^0.25 .- c0[2]) .- W[2];
    f[7] = (x[5]./x[6]) .- 4.*(sqrt.(0.5.*beta[3]./rho).*x[6].^0.25 .- c0[3]) .- W[3];
    f[8] = (x[7]./x[8]) .- 4.*(sqrt.(0.5.*beta[4]./rho).*x[8].^0.25 .- c0[4]) .- W[4];

    # return f
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
module Propagators

include("utils.jl")
include("impl.jl")
include("lib.jl")
include("integrations.jl")

end


text:
abstract ValueHistory
abstract UnivalueHistory{I} <: ValueHistory
abstract MultivalueHistory <: ValueHistory

Base.push!(history::UnivalueHistory, iteration, value) =
    throw(ArgumentError("The specified arguments are of incompatible type"))

# length(history::ValueHistory) =  error()
# enumerate(history::ValueHistory) = error()
# get(history::ValueHistory) = error()
# first(history::ValueHistory) = error()
# last(history::ValueHistory) = error()


text:
using FLICameras
using Test

@testset "FLICameras.jl" begin
    # Write your tests here.
end


text:
"""
    naive_lpd(log_likelihood::AbstractArray{<:Real}[, chain_index])

Calculate the naive (in-sample) estimate of the expected log probability density, otherwise
known as the in-sample Bayes score. This method yields heavily biased results, and we advise
against using it; it is included only for pedagogical purposes.

This method is unexported and can only be accessed by calling `ParetoSmooth.naive_lpd`.

# Arguments
  - $LOG_LIK_ARR
  - $CHAIN_INDEX_DOC
"""
function naive_lpd(log_likelihood::AbstractArray{<:Real, 3})
    @info "We advise against using `naive_lpd`, as it gives inconsistent and strongly " *
    "biased estimates. Use `psis_loo` instead."
    return _naive_lpd(log_likelihood)
end

function naive_lpd(
    log_likelihood::AbstractMatrix{<:Real},
    chain_index::AbstractVector = _assume_one_chain(log_likelihood)
)
    chain_index = Int.(chain_index)
    log_likelihood = _convert_to_array(log_likelihood, chain_index)
    return naive_lpd(log_likelihood)
end

function _naive_lpd(log_likelihood::AbstractArray{<:Real, 3})
    dims = size(log_likelihood)
    data_size = dims[1]
    mcmc_count = dims[2] * dims[3]
    log_count = log(mcmc_count)

    pointwise_naive = similar(log_likelihood, data_size)
    @tullio pointwise_naive[i] = exp(log_likelihood[i, j, k] - log_count) |> log
    return @tullio naive := pointwise_naive[i]
end

max_stars_repo_name: ParadaCarleton/JuLOOa

max_issues_repo_name: ParadaCarleton/JuLOOa

max_forks_repo_name: ParadaCarleton/JuLOOa


text:
always_print(::Matcher) = false

print_field{N}(m::Matcher, n::Type{Val{N}}) = "$(N)"

print_matcher(m::Matcher) = print_matcher(m, Set{Matcher}())

# this is optimized to be compact.  it could be prettier with more
# spaces, but then it's not as useful for large grammars.

function print_matcher(m::Matcher, known::Set{Matcher})
    function producer()
        if m in known
            produce("$(m.name)...")
        else
            produce("$(m.name)")
            if !always_print(m)
                push!(known, m)
            end
            names = filter(n -> n != :name, fieldnames(m))
            for name in names
                if isa(getfield(m, name), Matcher)
                    for (i, line) = enumerate(print_matcher(getfield(m, name), known))
                        if name == names[end]
                            produce(i == 1 ? "`-$(line)" : "  $(line)")
                        else
                            produce(i == 1 ? "+-$(line)" : "| $(line)")
                        end
                    end
                elseif isa(getfield(m, name), Array{Matcher,1})
                    for (j, x) in enumerate(getfield(m, name))
                        tag = name == :matchers ? "[$j]" : "$(name)[$j]"
                        for (i, line) = enumerate(print_matcher(getfield(m, name)[j], known))
                            if name == names[end] && j == length(getfield(m, name))
                                produce(i == 1 ? "`-$(tag):$(line)" : "  $(line)")
                            else
                                produce(i == 1 ? "+-$(tag):$(line)" : "| $(line)")
                            end
                        end
                    end
                else
                    if name == names[end]
                        produce("`-$(print_field(m, Val{name}))")
                    else
                        produce("+-$(print_field(m, Val{name}))")
                    end
                end
            end
        end
    end
    Task(producer)
end

function Base.print(io::Base.IO, m::Matcher)
    print(io, join(print_matcher(m), "\n"))
end