text:
"""
Update positions and velocities with a simple leapfrog langevin dynamics BBAOA integrator
"""
function refcpukernel!(pos,
    vel,
    specificforces,
    specificinteractions,
    specificparams,
    chain_force,
    perbead_params,
    global_params,
    invmass,
    γ,
    invβ,
    VNbeads::Val{Nbeads},
    step0::UInt64,
    time0,
    Vrngkey::Val{rngkey},
    Δt,
    VNdims::Val{Ndims},
    VChainbounds::Val{Chainbounds},
    VNsteps::Val{Nsteps}) where {Ndims, Nbeads, rngkey, Chainbounds, Nsteps}
    force= zero(vel)
    for step in 1:Nsteps
        time= (step-1)*Δt + time0
        force .= 0
        #add specific forces
        for (sf, si, sp) in zip(specificforces,specificinteractions,specificparams)
            for interaction in si
                rs= map(beadid->pos[beadid,:], interaction)
                fs= sf(rs,time,sp)
                for i in 1:length(interaction)
                    force[interaction[i],:] .+= fs[i]
                end
            end
        end
        #add chain forces and do kick/ drift BAOAB stuff
        for beadid in 1:Nbeads
            previd = prev_beadid(beadid,Chainbounds)
            nextid = next_beadid(beadid,Chainbounds)
            p= pos[previd,:]
            b= pos[beadid,:]
            n= pos[nextid,:]
            pf, bf, nf = chain_force(p,b,n,perbead_params[beadid], global_params)
            force[previd,:] .+= pf
            force[beadid,:] .+= bf
            force[nextid,:] .+= nf
        end
        for beadid in 1:Nbeads
            vel[beadid,:] .+= (invmass*Δt) .* force[beadid,:]
            pos[beadid,:] .+= (Δt/2) .* vel[beadid,:]
            vel[beadid,:] .*= exp(-Δt*γ)
            randns= randn_32(rngkey, step0+step, (beadid<<32)%UInt64,VNdims)
            vel[beadid,:] .+= (√(1-exp(-2γ*Δt))*√(invβ*invmass)) .* randns
            pos[beadid,:] .+= (Δt/2) .* vel[beadid,:]
        end
    end
end


text:
function prepareFWIDataFiles(m,Minv::RegularMesh,mref,boundsHigh,boundsLow,
							filenamePrefix::String,omega::Array{Float64,1},waveCoef::Array{Complex128,1},
							pad::Int64,ABLpad::Int64,jump::Int64,offset::Int64=prod(Minv.n+1),workerList = workers(),
							maxBatchSize::Int64=48, Ainv::AbstractSolver = getMUMPSsolver([],0,0,2),useFilesForFields::Bool = false)

########################## m is in Velocity here. ###################################
RCVfile = string(filenamePrefix,"_rcvMap.dat");
SRCfile = string(filenamePrefix,"_srcMap.dat");
writeSrcRcvLocFile(SRCfile,Minv,ABLpad,jump);
writeSrcRcvLocFile(RCVfile,Minv,ABLpad,1);

dataFullFilenamePrefix = string(filenamePrefix,"_freq");
gamma = prepareFWIDataFiles2(m, Minv, filenamePrefix,dataFullFilenamePrefix,omega,waveCoef,pad,ABLpad,offset,workerList,maxBatchSize,Ainv);

file = matopen(string(filenamePrefix,"_PARAM.mat"), "w");
write(file,"boundsLow",boundsLow);
write(file,"boundsHigh",boundsHigh);
write(file,"mref",mref);
write(file,"domain",Minv.domain);
write(file,"n",Minv.n);
write(file,"gamma",gamma);
write(file,"pad",pad);
write(file,"omega",omega);
write(file,"waveCoef",waveCoef);
close(file);
return gamma;
end

function prepareJointTravelTimeAndFWIDataFiles(m,Minv::RegularMesh,mref,boundsHigh,boundsLow,
							filenamePrefix::String,omega::Array{Float64,1},waveCoef::Array{Complex128,1},
							pad::Int64,ABLpad::Int64,jump::Int64,offset::Int64=prod(Minv.n+1),workerList = workers(),
							maxBatchSize::Int64=48, Ainv::AbstractSolver = getMUMPSsolver([],0,0,2),useFilesForFields::Bool = false)

########################## m is in Velocity here. ###################################
RCVfile = string(filenamePrefix,"_rcvMap.dat");
SRCfile = string(filenamePrefix,"_srcMap.dat");
writeSrcRcvLocFile(SRCfile,Minv,ABLpad,jump);
writeSrcRcvLocFile(RCVfile,Minv,ABLpad,1);

dataFullFilenamePrefix = string(filenamePrefix,"_freq");
gamma = prepareFWIDataFiles2(m, Minv, filenamePrefix,dataFullFilenamePrefix,omega,waveCoef,pad,ABLpad,offset,workerList,maxBatchSize,Ainv,useFilesForFields);
dataFullFilenamePrefix = string(filenamePrefix,"_travelTime");
HO = false;
prepareTravelTimeDataFiles(m, Minv, filenamePrefix,dataFullFilenamePrefix,offset,HO,useFilesForFields);

file = matopen(string(filenamePrefix,"_PARAM.mat"), "w");
write(file,"boundsLow",boundsLow);
write(file,"boundsHigh",boundsHigh);
write(file,"mref",mref);
write(file,"MinvOmega",Minv.domain);
write(file,"MinvN",Minv.n);
write(file,"gamma",gamma);
write(file,"pad",pad);
write(file,"omega",omega);
write(file,"waveCoef",waveCoef);
write(file,"HO",HO);
close(file);
end

function prepareFWIDataFiles2(m, Minv::RegularMesh, filenamePrefix::String,dataFullFilenamePrefix::String, omega::Array{Float64,1},
								waveCoef::Array{Complex128,1}, pad::Int64,ABLpad::Int64,offset::Int64,workerList::Array{Int64,1},maxBatchSize::Int64,
								Ainv::AbstractSolver,useFilesForFields::Bool = false)

println("maxOmega*maxKappaSq*h: should be below 0.6");
println(omega[end]*maximum(Minv.h)*(1./minimum(m)));

RCVfile = string(filenamePrefix,"_rcvMap.dat");
SRCfile = string(filenamePrefix,"_srcMap.dat");
srcNodeMap = readSrcRcvLocationFile(SRCfile,Minv);
rcvNodeMap = readSrcRcvLocationFile(RCVfile,Minv);

Q = generateSrcRcvProjOperators(Minv.n+1,srcNodeMap);
P = generateSrcRcvProjOperators(Minv.n+1,rcvNodeMap);
Q = Q.*1/(norm(Minv.h)^2);
println("We have ",size(Q,2)," sources");
# compute observed data

ABLamp = getMaximalFrequency(1./(minimum(m).^2),Minv);
gamma = getABL(Minv,true,ones(Int64,Minv.dim)*ABLpad,ABLamp);
attenuation = 0.01*ABLamp;
gamma += attenuation; # adding Attenuation.

println("~~~~~~~ Getting data FWI: ~~~~~~~");

# Solve forward problem (should be relaced in reading data from file)

batch = min(size(Q,2),maxBatchSize);
(pFor,contDiv,SourcesSubInd) = getFWIparam(omega,waveCoef,vec(gamma),Q,P,Minv,Ainv,workerList,batch,useFilesForFields);

(D,pFor) = getData(velocityToSlowSquared(m[:])[1],pFor,ones(length(pFor)),true);

nsrc = size(Q,2);
nrcv = size(P,2);

for k = 1:length(omega)
	I = contDiv[k] : contDiv[k+1] - 1
	Dobsk = Array{Array{Complex128,2}}(length(I));
	for i = 1:length(I)
		Dobsk[i] = fetch(D[I[i]]);
	end
	Dobsk = arrangeRemoteCallDataIntoLocalData(Dobsk);
	omRound = string(round((omega[k]/(2*pi))*100.0)/100.0);
	Wd_k = (1./(abs.(real(Dobsk))+0.1*mean(abs.(Dobsk)))) + 1im*(1./(abs.(imag(Dobsk))+0.1*mean(abs.(Dobsk))));
	# Wd_k = (1./(0.0*abs(real(Dobsk))+1.0*mean(abs(Dobsk)))) + 1im*(1./(0.0*abs(imag(Dobsk))+1.0*mean(abs(Dobsk))));
	Wd_k = limitDataToOffset(Wd_k,srcNodeMap,rcvNodeMap,offset);
	filename = string(dataFullFilenamePrefix,omRound,".dat");
	writeDataFile(filename,Dobsk,Wd_k,srcNodeMap,rcvNodeMap);
end
return gamma;
end

ext: jl


text:
mutable struct Event
    time::Int
    args::Dict{Symbol, Any}
end

Event(args::Dict{Symbol, Any}) = Event(0, args)


text:
# BayesTree declarations
type BayesTree
  bt
  btid::Int
  cliques::Dict{Int,Graphs.ExVertex}
  #edges
  frontals::Dict{ASCIIString,Int64}
end

function emptyBayesTree()
    bt =   BayesTree(Graphs.inclist(Graphs.ExVertex,is_directed=true),
                     0,
                     Dict{Int,Graphs.ExVertex}(),
                     #[],
                     Dict{AbstractString, Int}())
    return bt
end

# create a new clique
function addClique!(bt::BayesTree, fg::FactorGraph, varID::Int, condIDs::Array{Int}=Int[])
  bt.btid += 1
  clq = Graphs.add_vertex!(bt.bt, ExVertex(bt.btid,string("Clique",bt.btid)))
  bt.cliques[bt.btid] = clq
  clq.attributes["frontalIDs"] = Int[]
  clq.attributes["conditIDs"] = Int[]

  clq.attributes["label"] = ""
  clq.attributes["potentials"] = []

  appendClique!(bt, bt.btid, fg, varID, condIDs)
  return clq
end

# generate the label for particular clique -- graphviz drawing
function makeCliqueLabel(fgl::FactorGraph, bt::BayesTree, clqID::Int)
  clq = bt.cliques[clqID]
  flbl = ""
  clbl = ""
  for fr in clq.attributes["frontalIDs"]
    flbl = string(flbl,dlapi.getvertex(fgl,fr).attributes["label"], ",") #fgl.v[fr].
  end
  for cond in clq.attributes["conditIDs"]
    clbl = string(clbl, dlapi.getvertex(fgl,cond).attributes["label"], ",") # fgl.v[cond].
  end
  clq.attributes["label"] = string(flbl, ": ", clbl)
end

# add a conditional ID to clique
function appendConditional(bt::BayesTree, clqID::Int, condIDs::Array{Int,1})
  #@show "adding conditionals", condIDs
  clq = bt.cliques[clqID]
  clq.attributes["conditIDs"] = union(clq.attributes["conditIDs"], condIDs)
end

# Add a new frontal variable to clique
function appendClique!(bt::BayesTree, clqID::Int, fg::FactorGraph, varID::Int, condIDs::Array{Int,1}=Int[])
  clq = bt.cliques[clqID]
  var = dlapi.getvertex(fg, varID) # fg.v[varID]
  # add frontal variable
  push!(clq.attributes["frontalIDs"],varID)
  # total dictionary of frontals for easy access
  bt.frontals[var.attributes["label"]] = clqID#bt.btid

  appendConditional(bt, clqID, condIDs)
  makeCliqueLabel(fg, bt, clqID)
  Union{}
end

# instantiate a new child clique in the tree
function newChildClique!(bt::BayesTree, fg::FactorGraph, CpID::Int, varID::Int, Sepj::Array{Int,1})
  chclq = addClique!(bt, fg, varID, Sepj)
  parent = bt.cliques[CpID]
  # Staying with Graphs.jl for tree in first stage
  edge = Graphs.make_edge(bt.bt, parent, chclq)
  Graphs.add_edge!(bt.bt, edge)
  # addEdge!(bt.bt, parent, chclq)

  return chclq
end

# post order tree traversal and build potential functions
function findCliqueFromFrontal(bt::BayesTree, frtlID::Int)
  for cliqPair in bt.cliques
    id = cliqPair[1]
    cliq = cliqPair[2]
    for frtl in cliq.attributes["frontalIDs"]
      if frtl == frtlID
        return cliq
      end
    end
  end
  error("Clique with desired frontal ID not found")
end

# eliminate a variable for new
function newPotential(tree::BayesTree, fg::FactorGraph, var::Int, prevVar::Int, p::Array{Int,1})
    #@show fg.v[var].attributes["label"]
    firvert = dlapi.getvertex(fg,var)
    if (length(firvert.attributes["data"].separator) == 0) #fg.v[var].
      warn("newPotential -- sep length is 0")
      if (length(tree.cliques) == 0)
        addClique!(tree, fg, var)
      else
        appendClique!(tree, 1, fg, var) # add to root
      end
    else
      Sj = firvert.attributes["data"].separator  # fg.v[var].
      # Sj = fg.v[var].attributes["separator"]
      # find parent clique Cp that containts the first eliminated variable of Sj as frontal
      firstelim = 99999999999
      for s in Sj
        temp = findfirst(p, s)
        if (temp < firstelim)
          firstelim = temp
        end
      end
      felbl = dlapi.getvertex(fg,p[firstelim]).attributes["label"] # fg.v[p[firstelim]].
      CpID = tree.frontals[felbl]
      # look to add this conditional to the tree
      unFC = union(tree.cliques[CpID].attributes["frontalIDs"], tree.cliques[CpID].attributes["conditIDs"])
      if (sort(unFC) == sort(Sj))
        appendClique!(tree, CpID, fg, var)
      else
        newChildClique!(tree, fg, CpID, var, Sj)
      end
    end
end

# build the whole tree in batch format
function buildTree!(tree::BayesTree, fg::FactorGraph, p::Array{Int,1})
  rp = flipdim(p,1)
  prevVar = 0
  for var in rp
    newPotential(tree, fg, var, prevVar, p)
    prevVar = var
  end
end

# marg(x...) = -1, [0.0] #

## Find batch belief propagation solution
function prepBatchTree!(fg::FactorGraph; ordering::Symbol=:qr,drawpdf::Bool=false)
  p = getEliminationOrder(fg, ordering=ordering)
  #@show p=[1, 3, 8, 10, 5]
  # for v in p
  #     print("$(fg.v[v].label),")
  # end
  println()

  fge = deepcopy(fg)
  println("Building Bayes net...")
  buildBayesNet!(fge, p)

  tree = emptyBayesTree()
  buildTree!(tree, fge, p)

  println("Bayes Net")
  sleep(0.1)
  #fid = open("bn.dot","w+")
  #write(fid,to_dot(fge.bn))
  #close(fid)

  #GraphViz.Graph(to_dot(fge.bn))
  # Michael reference -- x2->x1, x2->x3, x2->x4, x2->l1, x4->x3, l1->x3, l1->x4
  println("Bayes Tree")
  if drawpdf
    fid = open("bt.dot","w+")
    write(fid,to_dot(tree.bt))
    close(fid)
    run(`dot bt.dot -Tpdf -o bt.pdf`)
  end

  # GraphViz.Graph(to_dot(tree.bt))
  #Michael reference 3sig -- x2l1x4x3    x1|x2

  println("Find potential functions for each clique")
  cliq = tree.cliques[1] # start at the root
  buildCliquePotentials(fg, tree, cliq); # fg does not have the marginals as fge does

  return tree
end

function resetData!(vdata::VariableNodeData)
  vdata.eliminated = false
  vdata.BayesNetOutVertIDs = Int64[]
  vdata.BayesNetVertID = 0
  vdata.separator = Int64[]
  nothing
end

function resetData!(vdata::FunctionNodeData)
  vdata.eliminated = false #f[2].attributes["eliminated"] = false
  vdata.potentialused = false #f[2].attributes["potentialused"] = false
  nothing
end

function resetFactorGraphNewTree!(fg::FactorGraph)
  for v in fg.g.vertices
    resetData!(getData(v))
    dlapi.updatevertex!(fg, v)
  end

  # for v in fg.v
  #   # v[2].attributes["data"].eliminated = false
  #   # # v[2].attributes["eliminated"] = false
  #   # v[2].attributes["data"].BayesNetOutVertIDs = Int64[]
  #   # v[2].attributes["data"].BayesNetVertID = 0
  #   # v[2].attributes["data"].separator = Int64[]
  #   # # v[2].attributes["BayesNetVert"] = Union{}
  #   # # v[2].attributes["separator"] = Int64[]
  #   resetData!(v[2].attributes["data"])
  #   dlapi.updatevertex!(fg, v[2])
  # end
  # for f in fg.f
  #   # f[2].attributes["data"].eliminated = false #f[2].attributes["eliminated"] = false
  #   # f[2].attributes["data"].potentialused = false #f[2].attributes["potentialused"] = false
  #   resetData!(f[2].attributes["data"])
  #   dlapi.updatevertex!(fg, f[2])
  # end

  nothing
end

function wipeBuildNewTree!(fg::FactorGraph; ordering=:qr)
  resetFactorGraphNewTree!(fg);
  return prepBatchTree!(fg, ordering=ordering);
end

function whichCliq(bt::BayesTree, frt::ASCIIString)
    bt.cliques[bt.frontals[frt]]
end

function getCliquePotentials!(fg::FactorGraph, bt::BayesTree, cliq::Graphs.ExVertex)
    frtl = cliq.attributes["frontalIDs"]
    cond = cliq.attributes["conditIDs"]
    allids = [frtl;cond]
    alldimIDs = Int[]
    for fid in frtl
      alldimIDs = [alldimIDs;dlapi.getvertex(fg,fid).attributes["data"].dimIDs] # ;fg.v[fid].
      # alldimIDs = [alldimIDs;fg.v[fid].attributes["dimIDs"]]
    end
    for cid in cond
      alldimIDs = [alldimIDs;dlapi.getvertex(fg,cid).attributes["data"].dimIDs] # ;fg.v[cid].
      # alldimIDs = [alldimIDs;fg.v[cid].attributes["dimIDs"]]
    end

    for fid in frtl
        usefcts = []
        for fct in dlapi.outneighbors(fg, dlapi.getvertex(fg,fid)) #out_neighbors(dlapi.getvertex(fg,fid),fg.g) # (fg.v[fid],
            #println("Checking fct=$(fct.label)")
            if fct.attributes["data"].potentialused!=true #fct.attributes["potentialused"]!=true ## USED TO HAVE == AND continue end here
                loutn = dlapi.outneighbors(fg, fct) #out_neighbors(fct, fg.g)
                if length(loutn)==1
                    appendUseFcts!(usefcts, fg.IDs[loutn[1].label], fct, fid) #out_neighbors(fct, fg.g)
                    # TODO -- make update vertex call
                    fct.attributes["data"].potentialused = true #fct.attributes["potentialused"] = true
                    dlapi.updatevertex!(fg, fct)
                end
                for sepSearch in loutn # out_neighbors(fct, fg.g)
                    if (fg.IDs[sepSearch.label] == fid)
                        continue # skip the fid itself
                    end
                    sea = findmin(abs(allids-fg.IDs[sepSearch.label]))
                    if sea[1]==0.0
                        appendUseFcts!(usefcts, fg.IDs[sepSearch.label], fct, fid)
                        # usefcts = [usefcts;(fg.IDs[sepSearch.label], fct, fid)]
                        fct.attributes["data"].potentialused = true #fct.attributes["potentialused"] = true
                        dlapi.updatevertex!(fg, fct)
                    end
                end
            end
        end
        cliq.attributes["potentials"]=[cliq.attributes["potentials"];usefcts]
    end
    return nothing
end

function getCliquePotentials!(fg::FactorGraph, bt::BayesTree, chkcliq::Int64)
    getCliquePotentials!(fg, bt.cliques[chkcliq])
end

function cliqPotentialIDs(cliq::Graphs.ExVertex)
  potIDs = Int[]
  for idfct in cliq.attributes["potentials"]
    push!(potIDs,idfct[1])
  end
  return potIDs
end

function collectSeparators(bt::BayesTree, cliq::Graphs.ExVertex)
  allseps = Int[]
  for child in out_neighbors(cliq, bt.bt)#tree
    # if length(child.attributes["conditIDs"]) > 0
      allseps = [allseps; child.attributes["conditIDs"]]
    # end
  end
  return allseps
end

function compCliqAssocMatrices!(bt::BayesTree, cliq::Graphs.ExVertex)
  frtl = cliq.attributes["frontalIDs"]
  cond = cliq.attributes["conditIDs"]
  inmsgIDs = collectSeparators(bt, cliq)
  potIDs = cliqPotentialIDs(cliq)
  # Construct associations matrix here
  # matrix has variables are columns, and messages/constraints as rows
  cols = [frtl;cond]
  cliq.attributes["inmsgIDs"] = inmsgIDs
  cliq.attributes["potIDs"] = potIDs
  cliqAssocMat = Array{Bool,2}(length(potIDs), length(cols))
  cliqMsgMat = Array{Bool,2}(length(inmsgIDs), length(cols))
  fill!(cliqAssocMat, false)
  fill!(cliqMsgMat, false)
  for j in 1:length(cols)
    for i in 1:length(inmsgIDs)
      if cols[j] == inmsgIDs[i]
        cliqMsgMat[i,j] = true
      end
    end
    for i in 1:length(potIDs)
      idfct = cliq.attributes["potentials"][i]
      if idfct[2].index == potIDs[i] # sanity check on clique potentials ordering
        for vertidx in idfct[2].attributes["data"].fncargvID #for vert in idfct[2].attributes["fnc"].Xi
          if vertidx == cols[j]
            cliqAssocMat[i,j] = true
          end
        end
        # for vert in idfct[2].attributes["data"].fnc.Xi #for vert in idfct[2].attributes["fnc"].Xi
        #   if vert.index == cols[j]
        #     cliqAssocMat[i,j] = true
        #   end
        # end
      else
        prtslperr("compCliqAssocMatrices! -- potential ID ordering was lost")
      end
    end
  end
  cliq.attributes["cliqAssocMat"] = cliqAssocMat
  cliq.attributes["cliqMsgMat"] = cliqMsgMat
  nothing
end

function getCliqMat(cliq::Graphs.ExVertex; showmsg=true)
  assocMat = cliq.attributes["cliqAssocMat"]
  msgMat = cliq.attributes["cliqMsgMat"]
  mat = showmsg ? [assocMat;msgMat] : assocMat
  # @show size(mat), size(assocMat), size(msgMat)
  return mat
end

function spyCliqMat(cliq::Graphs.ExVertex; showmsg=true)
  mat = getCliqMat(cliq, showmsg=showmsg)
  Gadfly.spy(mat)
end

function countSkips(bt::BayesTree)
  skps = 0
  for cliq in bt.cliques
    m = getCliqMat(cliq[2])
    mi = map(Int,m)
    skps += sum(map(Int,sum(mi,1) .== 1))
  end
  return skps
end

function skipThroughMsgsIDs(cliq::Graphs.ExVertex)
  numfrtl1 = floor(Int,length(cliq.attributes["frontalIDs"])+1)
  condAssocMat = cliq.attributes["cliqAssocMat"][:,numfrtl1:end]
  condMsgMat = cliq.attributes["cliqMsgMat"][:,numfrtl1:end]
  mat = [condAssocMat;condMsgMat];
  mab = sum(map(Int,mat),1) .== 1
  mabM = sum(map(Int,condMsgMat),1) .== 1
  mab = mab & mabM
  # rang = 1:size(condMsgMat,2)
  msgidx = cliq.attributes["conditIDs"][collect(mab)]
  return msgidx
end

function directAssignmentIDs(cliq::Graphs.ExVertex)
  assocMat = cliq.attributes["cliqAssocMat"]
  msgMat = cliq.attributes["cliqMsgMat"]
  mat = [assocMat;msgMat];
  mab = sum(map(Int,mat),1) .== 1
  mabA = sum(map(Int,assocMat),1) .== 1
  mab = mab & mabA
  frtl = cliq.attributes["frontalIDs"]
  cond = cliq.attributes["conditIDs"]
  cols = [frtl;cond]
  return cols[collect(mab)]
end

function directFrtlMsgIDs(cliq::Graphs.ExVertex)
  numfrtl = length(cliq.attributes["frontalIDs"])
  frntAssocMat = cliq.attributes["cliqAssocMat"][:,1:numfrtl]
  frtlMsgMat = cliq.attributes["cliqMsgMat"][:,1:numfrtl]
  mat = [frntAssocMat; frtlMsgMat];
  mab = sum(map(Int,mat),1) .== 1
  mabM = sum(map(Int,frtlMsgMat),1) .== 1
  mab = mab & mabM
  return cliq.attributes["frontalIDs"][collect(mab)]
end

function mcmcIterationIDs(cliq::Graphs.ExVertex)
  assocMat = cliq.attributes["cliqAssocMat"]
  msgMat = cliq.attributes["cliqMsgMat"]
  mat = [assocMat;msgMat];
  sum(sum(map(Int,mat),1)) == 0 ? error("mcmcIterationIDs -- unaccounted variables") : nothing
  mab = 1 .< sum(map(Int,mat),1)
  frtl = cliq.attributes["frontalIDs"]
  cond = cliq.attributes["conditIDs"]
  cols = [frtl;cond]
  return cols[collect(mab)]
end

function setCliqMCIDs!(cliq::Graphs.ExVertex)
  cliq.attributes["itervarIDs"] = mcmcIterationIDs(cliq)
  cliq.attributes["msgskipIDs"] = skipThroughMsgsIDs(cliq)
  cliq.attributes["directvarIDs"] = directAssignmentIDs(cliq)
  cliq.attributes["directFrtlMsgIDs"] = directFrtlMsgIDs(cliq)
  nothing
end

# function drawCliqAssocMat(fgl::FactorGrapha, cliq::Graphs.ExVertex)
#   cliqAssocMat = cliq.attributes["cliqAssocMat"]
#   ro, co = size(cliqAssocMat)
#   [@show fg.f[i].label for i in tree.cliques[1].attributes["potIDs"]];
#
#   showmat = Array{ASCIIString
#
# end

# post order tree traversal and build potential functions
function buildCliquePotentials(fg::FactorGraph, bt::BayesTree, cliq::Graphs.ExVertex)
    for child in out_neighbors(cliq, bt.bt)#tree
        buildCliquePotentials(fg, bt, child)
    end
    println("Get potentials $(cliq.attributes["label"])");
    getCliquePotentials!(fg, bt, cliq);

    compCliqAssocMatrices!(bt, cliq);
    setCliqMCIDs!(cliq);

    return Union{}
end


text:
# phylo/mutation_counting.jl
# ==========================
#
# Types and methods for counting different kinds of mutations between two
# DNA or RNA sequences.
#
# Part of the Bio.Phylo module.
#
# This file is a part of BioJulia. License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# Counting mutations
# ------------------

# Mutation types

abstract MutationType

"""
`AnyMutation` describes a site where two aligned nucleotides are not the
same.

Every kind of difference is counted.
"""
immutable AnyMutation <: MutationType end

"""
`TransitionMutation` describes a situation with two aligned nucleotides, where a
purine has mutated into another purine, or a pyrimadine has mutated into another
pyrimadine.

Possible transition mutations are:
A <-> G
C <-> T
"""
immutable TransitionMutation <: MutationType end

"""
`TransversionMutation` describes a situation with two aligned nucleotides,
where a purine has mutated into a pyrimadine or vice versa.

Possible transversion mutations are:
A <-> T
C <-> G
A <-> C
T <-> G
"""
immutable TransversionMutation <: MutationType end

"""
    is_ambiguous_strict{T<:Nucleotide}(a::T, b::T)

The strictest test of ambiguity, if either nucleotide is not
A, T, G, or C, then this function returns true.
"""
@inline function is_ambiguous_strict{T<:Nucleotide}(a::T, b::T)
    return isambiguous(a) | isambiguous(b)
end

"""
    is_mutation{T<:Nucleotide}(::Type{AnyMutation}, a::T, b::T)

Test if two nucleotides constitute a `AnyMutation`.
"""
@inline function is_mutation{T<:Nucleotide}(::Type{AnyMutation}, a::T, b::T)
    return a != b
end

"""
    is_mutation{T<:Nucleotide}(::Type{TransitionMutation}, a::T, b::T)

Test if two nucleotides constitute a `TransitionMutation`.
"""
@inline function is_mutation{T<:Nucleotide}(::Type{TransitionMutation}, a::T, b::T)
    return (a != b) & ((ispurine(a) & ispurine(b)) | (ispyrimidine(a) & ispyrimidine(b)))
end

"""
    is_mutation{T<:Nucleotide}(::Type{TransversionMutation}, a::T, b::T)

Test if two nucleotides constitute a `TransversionMutation`.
"""
@inline function is_mutation{T<:Nucleotide}(::Type{TransversionMutation}, a::T, b::T)
    return (a != b) & ((ispurine(a) & ispyrimidine(b)) | (ispyrimidine(a) & ispurine(b)))
end

"""
    flagmutations{M<:MutationType,N<:Nucleotide}(::Type{M}, seqs::Matrix{N})

For every pair of sequences, flag which base positions are mutations of type `T`.

This function returns a tuple of two matrices. Both matrices contain logical
values, and contain a column for every pairwise combination of sequences,
and a number of rows equal to the length of the aligned sequences.

For the first matrix, for a cell of a given column,
a value of `true` indicates that there is a mutation of type `T` between two
sequences at that base position.

For the second matrix, for a cell of a given column,
a value of `true` indicates that there is an ambiguous base in either of the two
sequences at that base position, consequently it is not possible to conclude
whether there is a mutation at that base position or not.

**Note: This method assumes that the sequences are stored in the `Matrix{N}`
provided as `seqs` in sequence major order i.e. each column of the matrix is one
complete nucleotide sequence.**

# Example

```julia
dnas = [dna"ATTG-ACCTGGNTTTCCGAA", dna"A-ACAGAGTATACRGTCGTC"]
2-element Array{Bio.Seq.BioSequence{Bio.Seq.DNAAlphabet{4}},1}:
 ATTG-ACCTGGNTTTCCGAA
 A-ACAGAGTATACRGTCGTC

julia> m = seqmatrix(dnas, :seq)
20×2 Array{Bio.Seq.DNANucleotide,2}:
  DNA_A    DNA_A
  DNA_T    DNA_Gap
  DNA_T    DNA_A
  DNA_G    DNA_C
  DNA_Gap  DNA_A
  DNA_A    DNA_G
  DNA_C    DNA_A
  DNA_C    DNA_G
  DNA_T    DNA_T
  DNA_G    DNA_A
  DNA_G    DNA_T
  DNA_N    DNA_A
  DNA_T    DNA_C
  DNA_T    DNA_R
  DNA_T    DNA_G
  DNA_C    DNA_T
  DNA_C    DNA_C
  DNA_G    DNA_G
  DNA_A    DNA_T
  DNA_A    DNA_C

julia> r = flagmutations(AnyMutation, m)
(
Bool[false; false; … ; true; true],

Bool[false; true; … ; false; false])

julia> r[1]
20×1 Array{Bool,2}:
 false
 false
  true
  true
 false
  true
  true
  true
 false
  true
  true
 false
  true
 false
  true
  true
 false
 false
  true
  true
```

In the above example of two sequences, positions 3:4, 7:8 are mutations and
positions 1:2, 5, 19:20 are not mutations.

"""
function flagmutations{M<:MutationType,N<:Nucleotide}(::Type{M}, seqs::Matrix{N})
    seqsize, nseqs = size(seqs)
    ismutant = Matrix{Bool}(seqsize, binomial(nseqs, 2))
    isambiguous = Matrix{Bool}(seqsize, binomial(nseqs, 2))
    col = 1
    @inbounds for i1 in 1:nseqs
        for i2 in i1+1:nseqs
            for s in 1:seqsize
                s1 = seqs[(i1 - 1) * seqsize + s]
                s2 = seqs[(i2 - 1) * seqsize + s]
                isamb = is_ambiguous_strict(s1, s2)
                ismut = is_mutation(M, s1, s2)
                isambiguous[(col - 1) * seqsize + s] = isamb
                ismutant[(col - 1) * seqsize + s] = !isamb & ismut
            end
            col += 1
        end
    end
    return ismutant, isambiguous
end

function flagmutations{M<:MutationType,A<:NucleotideAlphabet}(::Type{M}, seqs::Vector{BioSequence{A}})
    return flagmutations(M, seqmatrix(seqs, :seq))
end

"""
    count_mutations{T<:MutationType,N<:Nucleotide}(::Type{T}, seqs::Matrix{N})

Count the number of mutations between DNA sequences in a pairwise manner.

Different types of mutation can be counted:
`AnyMutation`, `TransitionMutation`, `TransversionMutation`.

Returns a tuple of: 1. A vector containing the number of mutations between each,
possible pair of sequences, and 2. a vector containing the number of sites
considered (sites with any ambiguity characters are not considered) for each
possible pair of sequences.

**Note: This method assumes that the sequences are stored in the `Matrix{N}`
provided as `seqs` in sequence major order i.e. each column of the matrix is one
complete nucleotide sequence.**
"""
function count_mutations{T<:MutationType,N<:Nucleotide}(::Type{T}, seqs::Matrix{N})
    # This method has been written with the aim of improving performance by taking
    # advantage of the memory layout of matrices of nucleotides, as well as
    # getting julia to emit simd code for the innermost loop.
    seqsize, nseqs = size(seqs)
    c = binomial(nseqs, 2)
    lengths = Vector{Int}(c)
    nmutations = Vector{Int}(c)
    target = 1
    @inbounds for i1 in 1:nseqs
        for i2 in i1+1:nseqs
            L = seqsize
            Nd = 0
            for s in 1:seqsize
                s1 = seqs[(i1 - 1) * seqsize + s]
                s2 = seqs[(i2 - 1) * seqsize + s]
                isamb = is_ambiguous_strict(s1, s2)
                ismut = is_mutation(T, s1, s2)
                L -= isamb
                Nd += !isamb & ismut
            end
            lengths[target] = L
            nmutations[target] = Nd
            target += 1
        end
    end
    return nmutations, lengths
end

"""
    count_mutations{T<:MutationType,A<:NucleotideAlphabet}(::Type{T}, sequences::Vector{BioSequence{A}})

Count the number of mutations between DNA sequences in a pairwise manner.

Different types of mutation can be counted:
`AnyMutation`, `TransitionMutation`, `TransversionMutation`.

Returns a tuple of: 1. A vector containing the number of mutations between each,
possible pair of sequences, and 2. a vector containing the number of sites
considered (sites with any ambiguity characters are not considered) for each
possible pair of sequences.
"""
function count_mutations{T<:MutationType,A<:NucleotideAlphabet}(::Type{T}, sequences::Vector{BioSequence{A}})
    seqs = seqmatrix(sequences, :seq)
    return count_mutations(T, seqs)
end

function flagmutations{N<:Nucleotide}(::Type{TransitionMutation}, ::Type{TransversionMutation}, seqs::Matrix{N})
    seqsize, nseqs = size(seqs)
    istransition = Matrix{Bool}(seqsize, binomial(nseqs, 2))
    istransversion = Matrix{Bool}(seqsize, binomial(nseqs, 2))
    isambiguous = Matrix{Bool}(seqsize, binomial(nseqs, 2))
    col = 1
    @inbounds for i1 in 1:nseqs
        s1offset = (i1 - 1) * seqsize
        for i2 in i1+1:nseqs
            s2offset = (i2 - 1) * seqsize
            resoffset = (col - 1) * seqsize
            for s in 1:seqsize
                s1 = seqs[s1offset + s]
                s2 = seqs[s2offset + s]
                isamb = is_ambiguous_strict(s1, s2)
                isdiff = s1 != s2
                ists = is_mutation(TransitionMutation, s1, s2)
                isambiguous[resoffset + s] = isamb
                istransition[resoffset + s] = !isamb & isdiff & ists
                istransversion[resoffset + s] = !isamb & isdiff & !ists
            end
            col += 1
        end
    end
    return istransition, istransversion, isambiguous
end

function flagmutations{A<:NucleotideAlphabet}(::Type{TransitionMutation}, ::Type{TransversionMutation}, seqs::Vector{BioSequence{A}})
    return flagmutations(TransitionMutation, TransversionMutation, seqmatrix(seqs, :seq))
end

"""
    count_mutations{N<:Nucleotide}(::Type{TransitionMutation}, ::Type{TransversionMutation}, sequences::Matrix{N})

Count the number of `TransitionMutation`s and `TransversionMutation`s in a
pairwise manner, between each possible pair of sequences.

Returns a tuple of: 1. A vector containing the number of transitions between
each, possible pair of sequences, and 2. a vector containing the number of
transversions between each, possible pair of sequences, and 3. a vector
containing the number of sites considered (sites with any ambiguity characters
are not considered) for each possible pair of sequences.

**Note: This method assumes that the sequences are stored in the `Matrix{N}`
provided as `seqs` in sequence major order i.e. each column of the matrix is one
complete nucleotide sequence.**
"""
function count_mutations{A<:Nucleotide}(::Type{TransitionMutation}, ::Type{TransversionMutation}, seqs::Matrix{A})
    # This method has been written with the aim of improving performance by taking
    # advantage of the memory layout of matrices of nucleotides, as well as
    # getting julia to emit simd code for the innermost loop.
    seqsize, nseqs = size(seqs)
    c = binomial(nseqs, 2)
    lengths = Vector{Int}(c)
    ntransition = Vector{Int}(c)
    ntransversion = Vector{Int}(c)
    target = 1
    @inbounds for i1 in 1:nseqs
        for i2 in i1+1:nseqs
            L = seqsize
            Ns = 0
            Nv = 0
            for s in 1:seqsize
                s1 = seqs[(i1 - 1) * seqsize + s]
                s2 = seqs[(i2 - 1) * seqsize + s]
                isamb = is_ambiguous_strict(s1, s2)
                isdiff = s1 != s2
                istrans = (ispurine(s1) & ispurine(s2)) | (ispyrimidine(s1) & ispyrimidine(s2))
                L -= isamb
                Ns += !isamb & isdiff & istrans
                Nv += !isamb & isdiff & !istrans
            end
            lengths[target] = L
            ntransition[target] = Ns
            ntransversion[target] = Nv
            target += 1
        end
    end
    return ntransition, ntransversion, lengths
end

function count_mutations{A<:Nucleotide}(::Type{TransversionMutation}, ::Type{TransitionMutation}, seqs::Matrix{A})
    return count_mutations(TransitionMutation, TransversionMutation, seqs)
end

"""
    count_mutations{A<:NucleotideAlphabet}(::Type{TransitionMutation}, ::Type{TransversionMutation}, sequences::Vector{BioSequence{A}})

Count the number of `TransitionMutation`s and `TransversionMutation`s in a
pairwise manner, between each possible pair of sequences.

Returns a tuple of: 1. A vector containing the number of transitions between
each, possible pair of sequences, and 2. a vector containing the number of
transversions between each, possible pair of sequences, and 3. a vector
containing the number of sites considered (sites with any ambiguity characters
are not considered) for each possible pair of sequences.
"""
function count_mutations{A<:NucleotideAlphabet}(::Type{TransitionMutation}, ::Type{TransversionMutation}, sequences::Vector{BioSequence{A}})
    seqs = seqmatrix(sequences, :seq)
    return count_mutations(TransitionMutation, TransversionMutation, seqs)
end

function count_mutations{A<:NucleotideAlphabet}(::Type{TransversionMutation}, ::Type{TransitionMutation}, sequences::Vector{BioSequence{A}})
    return count_mutations(TransitionMutation, TransversionMutation, sequences)
end


text:
export ValueFunctionArrayOnly,
    ValueFunction,
    EV_scaled_itp,
   dEV_scaled_itp

function check_a_eq_b_sans_2nd_to_last(a::Tuple, b::Tuple)
   length(b) == length(a)+1 || return false
   b[1:end-2] == a[1:end-1] || return false
   b[end] == a[end]         || return false
   return true
end

abstract type AbstractValueFunction end

# -----------------------------------------
# Value Function Arrays
# -----------------------------------------

"Value function arrays"
struct ValueFunctionArrayOnly{T<:Real} <: AbstractValueFunction
   EV::Array{T,3}
   dEV::Array{T,4}
   function ValueFunctionArrayOnly(EV, dEV)
      check_a_eq_b_sans_2nd_to_last(size(EV), size(dEV))  ||
        throw(DimensionMismatch("EV is $(size(EV)) but dEV is $(size(dEV))"))
      new{eltype(EV)}(EV, dEV)
   end
end

"Array suitable for VF iterations"
EV(  x::ValueFunctionArrayOnly) = x.EV
dEV( x::ValueFunctionArrayOnly) = x.dEV
size(x::ValueFunctionArrayOnly) = size(dEV(x))

function fill!(x::ValueFunctionArrayOnly, y)
    fill!(EV(x), y)
    fill!(dEV(x), y)
end

# -----------------------------------------
# Value Function Array w/ interpolation info
# -----------------------------------------

"Value Function arrays with ability to interpolate"
struct ValueFunction{T<:Real, TC<:Real, N, NP1,
    EV_ITP <:InPlaceInterp{T, TC, N,   Array{T,N}  },
    DEV_ITP<:InPlaceInterp{T, TC, NP1, Array{T,NP1}}
} <: AbstractValueFunction

    EV::EV_ITP
    dEV::DEV_ITP

    function ValueFunction(ev::EV_ITP, dev::DEV_ITP) where {
        T<:Real, TC<:Real, N, NP1,
        EV_ITP <:InPlaceInterp{T, TC, N,   Array{T,N}  },
        DEV_ITP<:InPlaceInterp{T, TC, NP1, Array{T,NP1}}
    }

        check_a_eq_b_sans_2nd_to_last(size(ev), size(dev))       || throw(DimensionMismatch())
        check_a_eq_b_sans_2nd_to_last(itpflag(ev), itpflag(dev)) || throw(error("itpflags are different"))
        check_a_eq_b_sans_2nd_to_last(ranges(ev), ranges(dev))   || throw(error("ranges are different"))

        return new{T,TC,N,NP1,EV_ITP,DEV_ITP}(ev, dev)
    end
end

EVobj( x::ValueFunction) = x.EV
dEVobj(x::ValueFunction) = x.dEV

function EV(x::ValueFunction)
    v = data(EVobj(x))
    nd = ndims(v)
    return reshape(v, :, size(v, nd-1), size(v, nd))
end

function dEV(x::ValueFunction)
    v = data(dEVobj(x))
    nd = ndims(v)
    return reshape(v, :, size(v, nd-2), size(v, nd-1), size(v, nd))
end

ValueFunctionArrayOnly(x::ValueFunction) = ValueFunctionArrayOnly(EV(x), dEV(x))

function fill!(x::ValueFunction, y)
    fill!(data( EVobj(x)), y)
    fill!(data(dEVobj(x)), y)
end

function update_interpolation!(x::ValueFunction, dograd)
    dograd && update_interpolation!(dEVobj(x))
    update_interpolation!(EVobj(x))
end

update_interpolation!(x...) = nothing
value_function(x...) = nothing

 EV_scaled_itp(x::ValueFunction) = scaled_interpolation( EVobj(x))
dEV_scaled_itp(x::ValueFunction) = scaled_interpolation(dEVobj(x))

# -----------------------------------------
# Outer constructors for VF
# -----------------------------------------

NoValueFunction(args...) = nothing

function ValueFunctionArrayOnly(reward, discount::T, statespace, zspace, ztransition, psispace, args...) where {T<:Real}
    nz = size(ztransition, 1)
    nψ = length(psispace)
    nK = _nparm(reward)
    nS = length(statespace)

    dev = zeros(T, nz, nψ, nK, nS)
     ev = zeros(T, nz, nψ,     nS)
    return ValueFunctionArrayOnly(ev, dev)
end

const DEFAULT_SPLINE = BSpline(Quadratic(Free(OnGrid())))
splinetype(r::AbstractRange, spline::Union{NoInterp,<:BSpline}=DEFAULT_SPLINE) = spline
splinetype(r::UnitRange, spline... ) = BSpline(Constant())
splinetype(r::AbstractUnitRange, spline...) = NoInterp()

function ValueFunction(ev, dev, reward, discount, statespace, zspace, ztransition, psispace, args...)
    nK = _nparm(reward)
    nS = length(statespace)
    nψ = length(psispace)
    nz = length.(zspace)

    ev_ranges  = (zspace..., psispace,            OneTo(nS))
    dev_ranges = (zspace..., psispace, OneTo(nK), OneTo(nS))

    ev_it  = splinetype.( ev_ranges, args...)
    dev_it = splinetype.(dev_ranges, args...)

    ev  = reshape(ev,  nz..., nψ,    nS)
    dev = reshape(dev, nz..., nψ, nK, nS)

    EV  = InPlaceInterp( ev,  ev_it,  ev_ranges)
    dEV = InPlaceInterp(dev, dev_it, dev_ranges)

    return ValueFunction(EV, dEV)
end

function ValueFunction(reward, discount::T, statespace, zspace, ztransition, psispace, args...) where {T<:Real}
    nK = _nparm(reward)
    nS = length(statespace)
    nψ = length(psispace)
    nz = length.(zspace)
    ev  = zeros(T, nz..., nψ,    nS)
    dev = zeros(T, nz..., nψ, nK, nS)
    return ValueFunction(ev, dev, reward, discount, statespace, zspace, ztransition, psispace, args...)
end


text:
using Statistics

function aux_stat(y)
    α = sqrt(mean(y.^2.0))
    y = abs.(y)
    m = mean(y)
    s = std(y)
    k = std(y.^2.0)
    # look for evidence of volatility clusters, for ρ
    mm = ma(y,5)
    mm = mm[5:end]
    clusters1 = quantile(mm,0.75)-quantile(mm, 0.25)
    mm = ma(y,10)
    mm = mm[10:end]
    clusters2 = quantile(mm,0.75)-quantile(mm, 0.25)
    # HAR model, for all params
    vcat(α, m, s, k, clusters1, clusters2, HAR(y))
end

#=
function RawMoment(θ, m)
    α = 2.0*log(θ[1])
    ρ = θ[2]
    σ = θ[3]
    exp(m*α/2.0 + (m^2.0)*σ/(1.0 - ρ^2.0)/4.0) # Shepphard GMM notes, page 384
end

function UnconditionalMoments(θ)
    α = 2.0*log(θ[1])
    ρ = θ[2]
    σ = θ[3]
    vcat(
        sqrt(2.0/pi)*RawMoment(θ,1),   # mean (abs(y))
        RawMoment(θ,2) # mean y^2
        )
        #2.0*sqrt(2.0/pi)*RawMoment(θ,3), # mean abs(y^3))
        #3.0*RawMoment(θ,4)               # mean y^4
        #)
end
=##=
function aux_stat(y, θ)
    α = sqrt(mean(y.^2.0))
    y = abs.(y)
    # look for evidence of volatility clusters, for ρ
    mm = ma(y,5)
    mm = mm[5:end]
    clusters = 0.0
    try
        clusters = quantile(mm,0.75)-quantile(mm, 0.25)
    catch
        clusters = 1.0
    end
    # HAR model, for all params
    vcat(α, clusters, HAR(y), UnconditionalMoments(θ))
end
=#

max_stars_repo_name: mcreel/SV

max_issues_repo_name: mcreel/SV

max_forks_repo_name: mcreel/SV


text:
using Documenter: deploydocs, makedocs
using QuerySQLite

makedocs(sitename="QuerySQLite.jl", modules=[QuerySQLite], doctest=false)
deploydocs(repo="github.com/queryverse/QuerySQLite.jl.git")


text:
#using PolyhedralOmega, LinearAlgebra, Test
using LinearAlgebra, Test, PolyhedralOmega

include("../src/SymbolicCone.jl")
@testset "Symbolic Cone Test Set" begin
    @test true
end

max_stars_repo_name: kevseryolcu/PolyhedralOmega

max_issues_repo_name: kevseryolcu/PolyhedralOmega

max_forks_repo_name: kevseryolcu/PolyhedralOmega


text:
function removeBubble(bubble, lc, rc)
    Lver, Rver = bubble.lver, bubble.rver
    para = bubble.parent.para
    chan = bubble.channel

    if NoBubble in para.filter
        if Lver.para.innerLoopNum == 0 && Rver.para.innerLoopNum == 0
            if chan == T || chan == U
                if lc == DI && rc == DI
                    return true
                end
            end
        end
    end

    return false
end

function notProper(para, K)
    if Proper in para.filter
        transferLoop = para.transferLoop
        @assert isempty(transferLoop) == false "Please initialize para.transferLoop to check proper diagrams."
        if transferLoop ≈ K
            return true
        end
    end
    return false
end

# check if G exist without creating objects in the pool
function isValidG(filter, innerLoopNum::Int)
    if (NoFock in filter) && (innerLoopNum == 1)
        return false
    end

    if (Girreducible in filter) && (innerLoopNum > 0)
        return false
    end

    return true
end

function isValidG(para::GenericPara)
    @assert para.diagType == GreenDiag
    return isValidG(para.filter, para.innerLoopNum)
end

function isValidSigma(filter, innerLoopNum::Int, subdiagram::Bool)
    @assert innerLoopNum >= 0
    if innerLoopNum == 0
        return false
    end
    if subdiagram && (Girreducible in filter)
        return false
    end

    if subdiagram && (NoFock in filter) && innerLoopNum == 1
        return false
    end

    return true
end

function isValidPolarization(filter, innerLoopNum::Int, subdiagram::Bool)
    @assert innerLoopNum >= 0
    if innerLoopNum == 0
        return false
    end
    if subdiagram && (Wirreducible in filter)
        return false
    end
    if subdiagram && (NoBubble in filer) && innerLoopNum == 1
        return false
    end
    return true
end


text:
module MITgcmTools

using Dates, DataFrames, NetCDF, Printf, MeshArrays, SparseArrays, Pkg.Artifacts

include("ReadFiles.jl")
include("FormatConversions.jl")
include("PhysicalOceanography.jl")
include("moreMITgcmTools.jl")

export MatrixInterp, convert2array, convert2gcmfaces
export read_bin, read_flt, read_nctiles, findtiles
export read_mdsio, read_meta, read_available_diagnostics
export SeaWaterDensity, MixedLayerDepth
export MITgcm_path, testreport
export extract_timeseries, matmul, nancount_gcmarray
export read_netcdf

p=dirname(pathof(MITgcmTools))
artifact_toml = joinpath(p, "../Artifacts.toml")
MITgcm_hash = artifact_hash("MITgcm", artifact_toml)
MITgcm_path = joinpath(artifact_path(MITgcm_hash)*"/","MITgcm-checkpoint67s/")

"""
    testreport(nam::String)

```
tmp=testreport("front_relax");
```
"""
function testreport(nm::String)
    c=`$(MITgcm_path)/verification/testreport -t $(MITgcm_path)/verification/$(nm)`
    run(c)
end

#more:
#
#using Plots; include("recipes_plots.jl"); export qwckplot

#deprecated:
#
#include("deprecated.jl"); export prep_MTRX, read_SPM

end # module


text:
using Knet
using CUDA

struct Conv3x3; w; act; f; end
(c::Conv3x3)(x) = if c.act c.f.(conv4(c.w, x, padding = 1, stride = 1)) else conv4(c.w, x, padding = 1, stride = 1) end
Conv3x3(cx::Int,cy::Int,act,f=relu) = Conv3x3(param(3,3,cx,cy),act,f)

struct Conv1x1; w; act; f; end
(c::Conv1x1)(x) = if c.act c.f.(conv4(c.w, x, padding = 0, stride = 1)) else conv4(c.w, x, padding = 0, stride = 1) end
Conv1x1(cx::Int,cy::Int,act,f=relu) = Conv1x1(param(1,1,cx,cy),act,f)

struct Dense; w; b; p; end
(d::Dense)(x) = d.w * mat(dropout(x,d.p)) .+ d.b # mat reshapes 4-D tensor to 2-D matrix so we can use matmul
Dense(i::Int,o::Int;pdrop=0) = Dense(param(o,i), param0(o), pdrop)

struct BNorm
    C
    moments
    params
    act
end
(bn::BNorm)(x) = if bn.act relu.(batchnorm(x, bn.moments, bn.params)) else batchnorm(x, bn.moments, bn.params) end

struct Upsample; stride; end
(u::Upsample)(x) = unpool(x, stride=u.stride)

struct Chain
    layers
    Chain(layers...) = new(layers)
end
(c::Chain)(x) = (for l in c.layers; x = l(x); end; x)

struct Conv; w; f; end
(c::Conv)(x) = c.f.(conv4(c.w, x, padding = 1, stride = 2))
Conv(w1::Int,w2::Int,cx::Int,cy::Int,f=relu) = Conv(param(w1,w2,cx,cy),f)

struct Dense1; w; b; f; p; end
(d::Dense1)(x) = d.f.(d.w * mat(dropout(x,d.p)) .+ d.b) # mat reshapes 4-D tensor to 2-D matrix so we can use matmul
Dense1(i::Int,o::Int,f=sigm;pdrop=0) = Dense1(param(o,i), param0(o), f, pdrop)


text:
module Fmincon

using MATLAB

function fmincon(fundef::String, fun::String, x0::AbstractArray{<:Real,1},
    lb::AbstractArray{<:Real,1}, ub::AbstractArray{<:Real,1};
    options::AbstractDict=Dict(),
    A::AbstractArray{<:Real,2}=zeros(0,0), b::AbstractArray{<:Real,1}=zeros(0),
    Aeq::AbstractArray{<:Real,2}=zeros(0,0), beq::AbstractArray{<:Real,1}=zeros(0),
    gradients::Bool=false, printfile::String="")

    # get the default matlab session
    msession = MATLAB.get_default_msession()

    # add matlab optimization function to the MATLAB path
    moptimpath = dirname(pathof(Fmincon))
    eval_string(msession, "addpath('$moptimpath');")

    # make sure inputs are compatible with MATLAB
    x0_m = Float64.(x0)
    A_m = Float64.(A)
    b_m = Float64.(b)
    Aeq_m = Float64.(Aeq)
    beq_m = Float64.(beq)
    lb_m = Float64.(lb)
    ub_m = Float64.(ub)
    options_m = Dict{String, Any}()
    for key in keys(options)
        if isa(options[key], Number)
            options_m[key] = Float64.(options[key])
        else
            options_m[key] = options[key]
        end
    end
    gradients_m = Bool(gradients)

    # # copy input variables to MATLAB session
    put_variable(msession, :x0, x0_m)
    put_variable(msession, :A, A_m)
    put_variable(msession, :b, b_m)
    put_variable(msession, :Aeq, Aeq_m)
    put_variable(msession, :beq, beq_m)
    put_variable(msession, :lb, lb_m)
    put_variable(msession, :ub, ub_m)
    put_variable(msession, :options, options_m)
    put_variable(msession, :gradients, gradients_m)

    # load objective/constraint function into embedded Julia
    if isfile(fundef)
        eval_string(msession, "jl.include('$fundef');")
    else
        eval_string(msession, "jl.eval('$fundef');")
    end

    nargout = gradients ? 6 : 3

    # make user function MEX-like
    mexusrfun =
        """
        jl.eval(...
            ['function mexusrfun(prhs)' newline...
             '    args = vec(jarray(prhs[1]))' newline...
             '    outs = $fun(args)' newline...
             '    return [Float64.(outs[i]) for i=1:$nargout]' newline...
             'end']);
        """
    eval_string(msession, mexusrfun)

    # create a function handle for the user function
    eval_string(msession, "mexusrfun = jl.wrapmex('mexusrfun', $nargout);")

    # clear outputs
    eval_string(msession, "clear xopt fopt exitflag output")

    # perform optimization
    eval_string(msession, "[xopt, fopt, exitflag, output] = optimize(mexusrfun, x0, A, b, Aeq, beq, lb, ub, options, gradients, '$printfile');")

    # check if optimization was performed succesfully
    eval_string(msession, "fail = exist('xopt') == 0;")
    fail = jvalue(get_mvariable(msession, :fail))

    # throw error if it didn't, and redirect user to other error messages
    if fail
        error("An error occured during call to fmincon, see error messages from MATLAB and embedded Julia")
    end

    # copy output variables from MATLAB session
    xopt = vec(jarray(get_mvariable(msession, :xopt)))
    fopt = jvalue(get_mvariable(msession, :fopt))
    exitflag = jvalue(get_mvariable(msession, :exitflag))
    output = jvalue(get_mvariable(msession, :output))

    return xopt, fopt, exitflag, output

end

end


text:
using ProgressMeter
using LinearAlgebra
#using MDToolbox

println(nprocs())

M = Matrix{Float64}[rand(1000,1000) for i = 1:50];

@time mvalues = @showprogress pmap(svdvals, M);

@everywhere f(x, k) = x .+ k

@time fvalues = progress_pmap(x -> f(x, 1.0), M);

max_stars_repo_name: matsunagalab/notebooks

max_issues_repo_name: matsunagalab/notebooks

max_forks_repo_name: matsunagalab/notebooks


text:
function writeresponse(response::HTTP.Messages.Response)
#TODO

end

function hasvalue(response)
    return haskey(response, "obj")
end

function parseresponse(response::HTTP.Messages.Response)
    bod = response.body
    code = response.status
    if code ∉ 200:399
        ret_dict = Dict{String, Any}()
        ret_dict["code"] = code
        ret_dict["url"] = response.request.target
        try
            ret_dict["return"] = JSON.parse(String(bod))
        catch e
            ret_dict["return"] = missing
        end
        return ret_dict
    end
    ret_dict = JSON.parse(String(bod))
    return ret_dict
end

function flattenlists!(responsedict)
    for dict in responsedict["obj"]
        for (key, value) in dict
            dict[key] = typeof(value) <: Array ? join(value, ", ") : value
        end
    end
end

function buildrequesturl(path)
    base = "https://api.chartmetric.com/api/"
    urls = join(path, '/')
    url_req =  base * urls
    return url_req
end

# Not exported
function ratelimitprotect(response::HTTP.Messages.Response; verbose = false)
    headers = Dict(response.headers)
    remaining = haskey(headers, "X-RateLimit-Remaining") ? parse(Int, headers["X-RateLimit-Remaining"]) : 0
    if remaining < 1
        resetin =  haskey(headers, "X-RateLimit-Reset") ? parse(Int, headers["X-RateLimit-Reset"]) - time() : 60
        sleepfor = maximum([0.1, resetin])
        verbose && println("Protecting rate limit")
        sleep((sleepfor + 5))
    end
    return nothing
end

function parsecmdate(date::String)
    date_tmp = string(SubString(date, 1:25))
    return parse(Dates.DateTime, date_tmp, Dates.RFC1123Format)
end


text:
using AROW
using Base.Test

using SVMLightLoader

function test_iris()
    iris = Array{Float64}[
        [5.1,3.5,1.4,0.2,-1],
        [4.9,3.0,1.4,0.2,-1],
        [4.7,3.2,1.3,0.2,-1],
        [4.6,3.1,1.5,0.2,-1],
        [5.0,3.6,1.4,0.2,-1],
        [5.4,3.9,1.7,0.4,-1],
        [4.6,3.4,1.4,0.3,-1],
        [5.0,3.4,1.5,0.2,-1],
        [4.4,2.9,1.4,0.2,-1],
        [4.9,3.1,1.5,0.1,-1],
        [5.4,3.7,1.5,0.2,-1],
        [4.8,3.4,1.6,0.2,-1],
        [4.8,3.0,1.4,0.1,-1],
        [4.3,3.0,1.1,0.1,-1],
        [5.8,4.0,1.2,0.2,-1],
        [5.7,4.4,1.5,0.4,-1],
        [5.4,3.9,1.3,0.4,-1],
        [5.1,3.5,1.4,0.3,-1],
        [5.7,3.8,1.7,0.3,-1],
        [5.1,3.8,1.5,0.3,-1],
        [5.4,3.4,1.7,0.2,-1],
        [5.1,3.7,1.5,0.4,-1],
        [4.6,3.6,1.0,0.2,-1],
        [5.1,3.3,1.7,0.5,-1],
        [4.8,3.4,1.9,0.2,-1],
        [5.0,3.0,1.6,0.2,-1],
        [5.0,3.4,1.6,0.4,-1],
        [5.2,3.5,1.5,0.2,-1],
        [5.2,3.4,1.4,0.2,-1],
        [4.7,3.2,1.6,0.2,-1],
        [4.8,3.1,1.6,0.2,-1],
        [5.4,3.4,1.5,0.4,-1],
        [5.2,4.1,1.5,0.1,-1],
        [5.5,4.2,1.4,0.2,-1],
        [4.9,3.1,1.5,0.1,-1],
        [5.0,3.2,1.2,0.2,-1],
        [5.5,3.5,1.3,0.2,-1],
        [4.9,3.1,1.5,0.1,-1],
        [4.4,3.0,1.3,0.2,-1],
        [5.1,3.4,1.5,0.2,-1],
        [5.0,3.5,1.3,0.3,-1],
        [4.5,2.3,1.3,0.3,-1],
        [4.4,3.2,1.3,0.2,-1],
        [5.0,3.5,1.6,0.6,-1],
        [5.1,3.8,1.9,0.4,-1],

        [5.7,2.8,4.5,1.3,1],
        [6.3,3.3,4.7,1.6,1],
        [4.9,2.4,3.3,1.0,1],
        [6.6,2.9,4.6,1.3,1],
        [5.2,2.7,3.9,1.4,1],
        [5.0,2.0,3.5,1.0,1],
        [5.9,3.0,4.2,1.5,1],
        [6.0,2.2,4.0,1.0,1],
        [6.1,2.9,4.7,1.4,1],
        [5.6,2.9,3.6,1.3,1],
        [6.7,3.1,4.4,1.4,1],
        [5.6,3.0,4.5,1.5,1],
        [5.8,2.7,4.1,1.0,1],
        [6.2,2.2,4.5,1.5,1],
        [5.6,2.5,3.9,1.1,1],
        [5.9,3.2,4.8,1.8,1],
        [6.1,2.8,4.0,1.3,1],
        [6.3,2.5,4.9,1.5,1],
        [6.1,2.8,4.7,1.2,1],
        [6.4,2.9,4.3,1.3,1],
        [6.6,3.0,4.4,1.4,1],
        [6.8,2.8,4.8,1.4,1],
        [6.7,3.0,5.0,1.7,1],
        [6.0,2.9,4.5,1.5,1],
        [5.7,2.6,3.5,1.0,1],
        [5.5,2.4,3.8,1.1,1],
        [5.5,2.4,3.7,1.0,1],
        [5.8,2.7,3.9,1.2,1],
        [6.0,2.7,5.1,1.6,1],
        [5.4,3.0,4.5,1.5,1],
        [6.0,3.4,4.5,1.6,1],
        [6.7,3.1,4.7,1.5,1],
        [6.3,2.3,4.4,1.3,1],
        [5.6,3.0,4.1,1.3,1],
        [5.5,2.5,4.0,1.3,1],
        [5.5,2.6,4.4,1.2,1],
        [6.1,3.0,4.6,1.4,1],
        [5.8,2.6,4.0,1.2,1],
        [5.0,2.3,3.3,1.0,1],
        [5.6,2.7,4.2,1.3,1],
        [5.7,3.0,4.2,1.2,1],
        [5.7,2.9,4.2,1.3,1],
        [6.2,2.9,4.3,1.3,1],
        [5.1,2.5,3.0,1.1,1],
        [5.7,2.8,4.1,1.3,1],
    ]

    iris_test = Array{Float64}[
        [4.8,3.0,1.4,0.3,-1],
        [5.1,3.8,1.6,0.2,-1],
        [4.6,3.2,1.4,0.2,-1],
        [5.3,3.7,1.5,0.2,-1],
        [5.0,3.3,1.4,0.2,-1],

        [7.0,3.2,4.7,1.4,1],
        [6.4,3.2,4.5,1.5,1],
        [6.9,3.1,4.9,1.5,1],
        [5.5,2.3,4.0,1.3,1],
        [6.5,2.8,4.6,1.5,1],
    ]

    arow = AROW.Classifier(4)

    shuffle!(iris)
    for v in iris
        fit(arow, v[1:4], Int(v[5]))
    end

    for v in iris_test
        @test predict(arow, v[1:4]) == Int(v[5])
    end

    arow1 = AROW.Classifier(4)
    arow2 = AROW.Classifier(4)
    shuffle!(iris)
    for v = iris[1:22]
        fit(arow1, v[1:4], Int(v[5]))
    end
    for v = iris[23:45]
        fit(arow2, v[1:4], Int(v[5]))
    end
    composition!(arow1, arow2)
    for v in iris_test
        @test predict(arow1, v[1:4]) == Int(v[5])
    end
end

function test_news20()
    dim = 62061
    arow = AROW.MultiClassifier(dim, 20, 0.1)

    for (v, label) in SVMLightFile("data/news20.scale", dim)
        fit(arow, vec(full(v)), Int(label))
    end

    num = 0
    correct = 0
    for (v, label) in SVMLightFile("data/news20.t.scale", dim)
        result = predict(arow, vec(full(v)))

        if(result == Int(label))
            correct += 1
        end
        num += 1
    end

    @test correct / num >= 0.7
end

test_iris()
test_news20()


text:
#softmax function: f(x)=exp(x)/sum(exp(x))

Fsoftmax(x::Array{Float64,2})=(exp(x)./sum(exp(x),1),[]) # TODO: better to subtract the max of x to make it numerically more stable

function Fsoftmax_inplace(value,auxvalue,x::Array{Float64,2})
    copy!(value,exp(x)./sum(exp(x),1)) # TODO: better to subtract the max of x to make it numerically more stable
end

function Dsoftmax(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::Array{Float64,2})
    axpy!(1.0,f_c.*(grad_c-repmat(sum(f_c.*grad_c,1),size(f_c,1),1)),grad_n)
end

if PROC=="GPU"

    function Fsoftmax(X::CudaArray)
        expX=CudaArray(Float64,size(X))
        exp!(X,expX)
        onr=CudaArray(Float64,(1,size(X,1))); fill!(onr,1.0);
        colsum=onr*expX
        onc=CudaArray(Float64,(size(X,1),1)); fill!(onc,1.0);
        out=CudaArray(Float64,size(X))
        vdiv!(1.0,expX,onc*colsum,out)
        free(expX); free(onr); free(colsum); free(onc)
        return (out,[]) # memory leak -- how to free out
    end

    function softmax!(X::CudaArray,out::CudaArray)
        expX=CudaArray(Float64,size(X))
        exp!(X,expX)
        onr=CudaArray(Float64,(1,size(X,1))); fill!(onr,1.0);
        colsum=onr*expX
        onc=CudaArray(Float64,(size(X,1),1)); fill!(onc,1.0);
        vdiv!(1.0,expX,onc*colsum,out)
        free(expX); free(onr); free(colsum); free(onc)
    end
    export softmax!

    function Fsoftmax_inplace(value,auxvalue,X::CudaArray)
        expX=CudaArray(Float64,size(X))
        exp!(X,expX)
        onr=CudaArray(Float64,(1,size(X,1))); fill!(onr,1.0);
        colsum=onr*expX
        onc=CudaArray(Float64,(size(X,1),1)); fill!(onc,1.0);
        vdiv!(1.0,expX,onc*colsum,value)
        free(expX); free(onr); free(colsum); free(onc)
    end

    function Dsoftmax(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray)
        tmp=CudaArray(Float64,size(f_c))
        vmult!(1.0,f_c,grad_c,tmp)
        onr=CudaArray(Float64,(1,size(tmp,1))); fill!(onr,1.0);
        colsum=onr*tmp
        onc=CudaArray(Float64,(size(tmp,1),1)); fill!(onc,1.0);
        copy!(tmp,grad_c)
        axpy!(-1.0,onc*colsum,tmp)
        vmult!(1.0,f_c,tmp,tmp)
        axpy!(1.0,tmp,grad_n)
        free(tmp); free(onr); free(onc);
    end

end

Derivative[Fsoftmax]=Dsoftmax
Inplace[Fsoftmax]=Fsoftmax_inplace

ADsoftmax(n)=ADnode(Fsoftmax,n)

softmax(n::ADnode)=ADnode(Fsoftmax,n)
softmax(x::Array)=exp(x)./sum(exp(x),1) # TODO: better to subtract the max of x to make it numerically more stableFsoftmax(x)

export softmax

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
module ModBC

using ..ModLattice: LatticeD2Q4,LatticeD2Q5,LatticeD2Q9, opp_lat_vec
using ..ModParameters: Param, BCDirichlet, BCRiemann
using ..ModMesh: Mesh, East_,West_,North_,South_

export setBCs!

"""
    setBCs!(lat,param,mesh)

Set the boundary conditions. Note that BCDir: right 1, left 2, top 3, bottom 4
This is consistent with the lattice direction, but different from the
connectivity orientation.
"""
function setBCs!(lat::LatticeD2Q4,param::Param,mesh::Mesh)
   Temp = param.Temp
   BCtype = param.BCtype
   BCDir  = mesh.BCDir
   BCnode = mesh.BCnode
   Conn   = mesh.Connectivity
   w,nw = lat.w, lat.nw

   for iBC = 1:length(BCtype)
      if BCtype[iBC] == BCDirichlet
         i = BCDir[iBC]
         j = opp_lat_vec(lat,i)
         lat.f[j,BCnode[iBC]] = (w[i]+w[j])*Temp[iBC] .- lat.f[i,BCnode[iBC]]

      elseif BCtype[iBC] == BCRiemann
         i = BCDir[iBC]
         for iBCnode in BCnode[iBC]
            lat.f[:,iBCnode] = lat.f[:,Conn[i,iBCnode]]
         end
      end
   end

end

function setBCs!(lat::LatticeD2Q5,param::Param)
   T = param.T
   nI,nJ = param.nI,param.nJ
   BCtype= param.BCtype
   w = lat.w

   if BCtype == BCDirichlet

   elseif BCtype == BCRiemann

   end
end

function setBCs!(lat::LatticeD2Q9,param::Param)
   Tleft,Tright,Ttop,Tbot = param.Tleft, param.Tright, param.Ttop, param.Tbot
   nI,nJ = param.nI,param.nJ
   BCtype= param.BCtype
   w = lat.w

   if BCtype == BCDirichlet

   elseif BCtype == BCRiemann

   end
end

end


text:
using ADCME
using PyPlot
using AdFem

Random.seed!(233)
idx = fem_randidx(100, m, n, h)

function poisson1d(c)
    m = 40
    n = 40
    h = 0.1
    bdnode = bcnode("all", m, n, h)
    c = constant(c)
    xy = gauss_nodes(m, n, h)
    κ = exp(c[1] + c[2] * xy[:,1] + c[3]*xy[:,2])
    κ = compute_space_varying_tangent_elasticity_matrix(κ, m, n, h)
    K = compute_fem_stiffness_matrix1(κ, m, n, h)
    K, _ = fem_impose_Dirichlet_boundary_condition1(K, bdnode, m, n, h)
    rhs = compute_fem_source_term1(ones(4m*n), m, n, h)
    rhs[bdnode] .= 0.0
    sol = K\rhs
    sol[idx]
end

c = Variable(rand(3))
y = poisson1d(c)
Γ = gradients(y, c)
Γ = reshape(Γ, (100, 3))

# generate data
sess = Session(); init(sess)
run(sess, assign(c, [1.0;2.0;3.0]))
obs = run(sess, y) + 1e-3 * randn(100)

# Inverse modeling
loss = sum((y - obs)^2)
init(sess)
BFGS!(sess, loss)

y = obs
H = run(sess, Γ)
R = (2e-3)^2 * diagm(0=>ones(100))
X = run(sess, c)
Q = diagm(0=>ones(3))
m, V = uqlin(y, H, R, X, Q)
plot([1;2;3], [1.;2.;3.], "o", label="Reference")
errorbar([1;2;3],m + run(sess, c), yerr=2diag(V), label="Estimated")
legend()

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
"""
Rational quardratic kernel
"""
struct ArdRQKernel{T, VT<:AbstractVector{T}} <: AbstractKernel
	ll::VT
	lα::VT
	lσ::VT
end
@functor ArdRQKernel
function (RQK::ArdRQKernel{T, VT})(x; λ=T(1e-6)) where {T, VT}
	length(RQK.ll) == size(x, 1) || throw(DimensionMismatch("size of length scale parameter should be the same as number of features"))
	n = size(x, 2)
	scaled_x = @. x*exp(-RQK.ll)
	d = square_metric(scaled_x)
	σ_square = exp(2*RQK.lσ[1])
	σ_square*power(1.0 .+ 0.5*exp(-RQK.lα[1])*d, -exp(RQK.lα[1])) + Diagonal(λ*ones(n))
end
function (RQK::ArdRQKernel{T, VT})(x, xo) where {T, VT}
	σ_square = exp(2*RQK.lσ[1])
	scaled_x = @. x*exp(-RQK.ll)
	scaled_xo = @. xo*exp(-RQK.ll)
	d = square_metric(scaled_x, scaled_xo)
	σ_square*power(1.0 .+ 0.5*exp(-RQK.lα[1])*d, -exp(RQK.lα[1]))
end

function ArdRQKernel(n_features::Int)
	l = rand(n_features)
	lα = rand(1)
	lσ = rand(1)
	ArdRQKernel(ll, lα, lσ)
end

reset(RQK::ArdRQKernel) = ArdRQKernel(length(RQK.ll))

function Base.show(io::IO, kernel::ArdRQKernel)
	print(io, "ArdRQKernel(")
	print(io, "ll=", kernel.ll, ", ", "lα=", kernel.lα, ", ", "lσ=", kernel.lσ)
	print(")")
end

## derived kernel
struct IsoRQKernel{T, VT<:AbstractVector{T}} <: AbstractKernel
	ll::VT
	lα::VT
	lσ::VT
end
@functor IsoRQKernel
function (IsoRQK::IsoRQKernel{T, VT})(x; λ=T(1e-6)) where {T, VT}
	n = size(x, 2)
	scaled_x = x .* exp(-IsoRQK.ll[1])
	d = square_metric(scaled_x)
	σ_square = exp(2*IsoRQK.lσ[1])
	σ_square*power(1.0 .+ 0.5*exp(-IsoRQK.lα[1])*d, -exp(IsoRQK.lα[1])) + Diagonal(λ*ones(n))
end
function (IsoRQK::IsoRQKernel{T, VT})(x, xo) where {T, VT}
	σ_square = exp(2*IsoRQK.lσ[1])
	scaled_x = x .* exp(-IsoRQK.ll[1])
	scaled_xo = xo .* exp(-IsoRQK.ll[1])
	d = square_metric(scaled_x, scaled_xo)
	σ_square*power(1.0 .+ 0.5*exp(-IsoRQK.lα[1])*d, -exp(IsoRQK.lα[1]))
end

function IsoRQKernel()
	ll = rand(1)
	lα = rand(1)
	lσ = rand(1)
	IsoRQKernel(ll, lα, lσ)
end

reset(IsoRQK::IsoRQKernel) = IsoRQKernel()

function Base.show(io::IO, kernel::IsoRQKernel)
	print(io, "IsoRQKernel(")
	print(io, "ll=", kernel.ll[1], ", ", "lα=", kernel.lα[1], ", ", "lσ=", kernel.lσ[1])
	print(")")
end

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
time_namespace_modify = Namespace(
    "time" => AetherlangObject((line::Ref{UInt}, ns::Namespace, dummy::AetherlangObject)->AetherlangObject(time())),
    "c.time" => AetherlangObject((line::Ref{UInt}, ns::Namespace, n::AetherlangObject)->AetherlangObject(ccall(:time, (UInt64), (Int,), n.dataref)))
)

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
"""
    block_spec(λmin, λmax, Δλ, λs, λe, str = 1.0, mode = 0)

Regular `block spectrum` initialized to `0.0`, with blocks `λs..λe = 1.0`

    block_spec(λmin, λmax, Δλ, λs, λe, str = `1.0`, 1)

Inversed `block spectrum` initialized to `1.0`, with blocks `λs..λe = 0.0`
"""
function block_spec(λmin, λmax, Δλ, λs, λe, str = 1.0, mode = 0)
    r = mode == 0 ? (zeros(Float64, Int(1.0 + (λmax - λmin) / Δλ))) : (ones(Float64, Int(1.0 + (λmax - λmin) / Δλ)))
    for n = λs:Δλ:λe
        r[Int((n-λmin)/Δλ)+1] = mode==0 ? str : 0.0
    end
    λ = collect(λmin:Δλ:λmax)
    reflectance_spec(λ, r)
end

"""
    block_spec(env::SpecEnvironment, λ)

Zero block spectrum with only wavelength `λ` set to `1.0`.

Fast version to create highest purity optimal colors. Linear interpolation for intermediate `λ` values.
"""
function block_spec(env::SpecEnvironment,λ::Real)
    env.λmin ≤ λ ≤ env.λmax ? nothing : DomainError(λ, "Outside of boundaries defined by spectral environment.")
    lam = collect(env.λmin:env.Δλ:env.λmax)
    lo = indexin(floor(Float64(λ)), lam)[1]
    s = zeros(Float64, length(lam))
    hi = lo + 1
    dh = lam[hi] - λ
    dl = λ - lam[lo]
    il = (dl / dh + 1)^-1
    s[lo] = il
    s[hi] = 1 - il
    reflectance_spec(lam, s)
end

"""
    block_spec(λ)

Zero block spectrum with only wavelength `λ` set to `1.0`
Fast version to create highest purity optimal colors.

Spectral limits are taken from the spectral environment.
"""
block_spec(λ::Real) = block_spec(SPECENV, λ)

"""
    led_spec(λmin, λmax, Δλ, λ0, Δλ1_2)

LED spectrum simulation

- `λ0`: peak wavelength
- `Δλ1_2`: half spectral width, width at 50% peak intensity

Result: `LSpec`

Source: Yoshi Ohno (2004), Simulation Analysis of White LED Spectra and Color Rendering,
        Proceedings for CIE Expert Sumposium on LED Light Sources, Tokyo, JA (Accessed May 27, 2021)
"""
function led_spec(λmin, λmax, Δλ, λ0, Δλ1_2)
    λ = collect(λmin:Δλ:λmax)
    n = length(λ)
    l = zeros(n)
    for i in 1:n
        l[i] = led(λ[i], λ0, Δλ1_2)
    end
    luminance_spec(λ,l)
end

"""
LED helper function
"""
function led(λ,λ0, Δλ1_2)
    g = exp(-((λ - λ0) / Δλ1_2)^2)
    (g + 2 * g^5) / 3
end

"""
    D_series_generator(T::Real)

The relative spectral power distribution SPD(λ) of a D series illuminant at CCT `T`, from `300` to `830 nm` in `5 nm` steps. By definition the output spectra are normalized in a way that their luminance is always `100.0` at a wavelength of `560 nm`.

Allowed CCT range: `4000 K` ≤ `CCT` ≤ `25000 K`
"""
function D_series_generator(T::Real)
    #xD,yD = CIE Daylight Locus
    if T < 4000 || T > 25000
        throw(DomainError(T, "CCT must be between 4000 and 25000 K."))
    elseif 4000 ≤ T ≤ 7000
        xD = 0.244063 + 0.09911 * 1000.0 / T + 2.9678 * 1_000_000.0 / T^2 - 4.6070 * 1_000_000_000.0 / T^3
    elseif 7000 < T ≤ 25000
        xD = 0.237040 + 0.24748 * 1000.0 / T + 1.9018 * 1_000_000.0 / T^2 - 2.0064 * 1000_000_000.0 / T^3
    end

    yD = -3.000 * xD^2 + 2.870 * xD - 0.275
    M  = 0.0241 + 0.2562 * xD - 0.7341 * yD
    M1 = (-1.3515 - 1.7703 * xD + 5.9114 * yD) / M
    M2 = (0.03000 - 31.4424 * xD + 30.0717 * yD) / M
    S0 = daylight_generator_table[:, 2]
    S1 = daylight_generator_table[:, 3]
    S2 = daylight_generator_table[:, 4]
    SD = S0 .+ M1 * S1 .+ M2 * S2
    luminance_spec(daylight_generator_table[:, 1], SD)
end

"""
    D_series_whitepoint(env::SpecEnv,T::Real)

produces the white point according to the CIE Daylight Series.

result: `::XYZ`
"""
function D_series_whitepoint(env::SpecEnv, T::Real)
    ×(D_series_illuminant(env,T), env.cmf)
end

"""
    D_series_whitepoint(T::Real)

produces the white point according to the CIE Daylight Series.

result: `::XYZ`
"""
function D_series_whitepoint(T::Real)
    D_series_whitepoint(SPECENV,T)
end

"""
    D_series_illuminant(env::SpecEnv, T::Real)

produces CIE Daylight spectrum according to desired CCT.

result: `::LSpec`
"""
function D_series_illuminant(env::SpecEnv, T::Real)
    dlum = D_series_generator(T)
    if dlum.λ[1] < env.λmin || dlum.λ[end] > env.λmax
        extr = extrap(linearspline(dlum.λ,dlum.l),collect(env.λmin:env.Δλ:env.λmax), :linear)
        intr = interp(extr,collect(env.λmin:env.Δλ:env.λmax))
    else
        intr = interp(dlum,collect(env.λmin:env.Δλ:env.λmax))
    end
    luminance_spec(intr[1], intr[2])
end

"""
    D_series_illuminant(T::Real)

produces CIE Daylight spectrum according to desired CCT.

result: `::LSpec`
"""
function D_series_illuminant(T::Real)
    D_series_illuminant(SPECENV, T)
end

# Selection of standard daytime illuminants using proper CCTs

"""
    D40_illuminant()

Result: `::LSpec` of D40 illuminant, with a CCT of  4003 K
"""
D40_illuminant() = D_series_proper_illuminant(4000)

"""
    D50_illuminant()

Result: `::LSpec` of D50 illuminant, with a CCT of  5003 K
"""
D50_illuminant() = D_series_proper_illuminant(5000)

"""
    D55_illuminant()

Result: `::LSpec` of D55 illuminant, with a CCT of  5503 K
"""
D55_illuminant() = D_series_proper_illuminant(5500)

"""
    D65_illuminant()

Result: `::LSpec` of D55 illuminant, with a CCT of  6504 K
"""
D65_illuminant() = D_series_proper_illuminant(6500)

"""
    D75_illuminant()

Result: `::LSpec` of D75 illuminant, with a CCT of 7504 K
"""
D75_illuminant() = D_series_proper_illuminant(7500)

"""
    D93_illuminant()

Result: `::LSpec` of D93 illuminant, with a CCT of 9304 K
"""
D93_illuminant() = D_series_proper_illuminant(9300)

# General standard daytime illuminants with proper CCTs

"""
    D_series_proper_illuminant(T::Real)

Generates a D-series illuminant with the proper CCT. For example, a D65 illuminant has a proper CCT of 6504 K. CCT according to the International Practical Temperature Scale from 1968.
"""
function D_series_proper_illuminant(T::Real)
    D_series_illuminant(T * (1.4388 / 1.4380))
end

"""
    normalized_D_series(env::SpecEnv = SPECENV, T::Real)

produces CIE Daylight spectrum according to desired CCT, normalized to l=100 at λ=560 nm.

# Examples
```jldoctest
julia> normalized_D_series(SPECENV, 6500)
LSpec(Real[390.0, 391.0, 392.0, 393.0, 394.0, 395.0, 396.0, 397.0, 398.0, 399.0  …  821.0, 822.0, 823.0, 824.0, 825.0, 826.0, 827.0, 828.0, 829.0, 830.0], Real[0.5460212796545921, 0.5650827295052294, 0.590288902296674, 0.620011293254124, 0.6526213976027774, 0.6864907105678326, 0.7200645866061678, 0.7520838171013811, 0.7813630526687508, 0.8067169439235554  …  0.578028743266789, 0.5810475635333059, 0.5837744258575818, 0.5863553092107371, 0.5889361925638924, 0.5916338590939442, 0.5944483088008923, 0.5973503458905127, 0.6003107745685814, 0.6033003990408742])

normalized_D_series(9300)
LSpec(Real[390.0, 391.0, 392.0, 393.0, 394.0, 395.0, 396.0, 397.0, 398.0, 399.0  …  821.0, 822.0, 823.0, 824.0, 825.0, 826.0, 827.0, 828.0, 829.0, 830.0], Real[0.9296442162075731, 0.9527892027932645, 0.9854018656135435, 1.0248615127313507, 1.0685474522096268, 1.1138389921113123, 1.158318458802078, 1.2003802518585145, 1.2386217891599418, 1.2716404885856807  …  0.4869584146523383, 0.4894740470571454, 0.4917569270423351, 0.4939234308177161, 0.4960899345930971, 0.49834953933632503, 0.5007022450474, 0.5031247764843602, 0.5055938584052438, 0.5080862155680892])
```
"""
function normalized_D_series(env::SpecEnv, T::Real)
    nT = normalize_D_series(D_series_generator(T))
    sT = SSpline.cubicspline(nT.λ, nT.l,:natural)
    iT = SSpline.interp(sT, collect(env.λmin:env.Δλ:env.λmax))
    luminance_spec(iT[1], iT[2])
end

normalized_D_series(T::Real) = normalized_D_series(SPECENV, T)

"""
    blackbody_illuminant(env::SpecEnv,T::Real)

produces blackbody spectrum according to desired blackbody temperature in K.

# Examples
```jldoctest
```
"""
function blackbody_illuminant(env::SpecEnv,T::Real)
    h = 6.62607015e-34 # [J*s] Planck’s constant
    k = 1.380649e-23   # [J/K] Boltzmann constant
    c = 299792458      # [m/s] Light speed
    wl=collect(env.λmin:env.Δλ:env.λmax)
    il=zeros(length(wl))
    @inbounds for n in 1:length(wl)
        λ = wl * 1e-9
        il[n] = (2 * π * h * c^2) / (λ[n]^5 * (exp(h * c / (λ[n] * k * T)) - 1))
    end
    luminance_spec(wl, il)
end

blackbody_illuminant(T::Real) = blackbody_illuminant(SPECENV, T)

"""
`blackbody_whitepoint(env::SpecEnv, T::Real)`

produces the white point according to desired blackbody temperature in K.

# Examples
```jldoctest
```
"""
function blackbody_whitepoint(env::SpecEnv, T::Real)
    ×(normalize_spec(blackbody_illuminant(env, T)), env.cmf)
end

# Illuminant A spectrum
# Source: Dr. H. W. G. Hunt, “Measuring Colour”, third edition
# Appendix 5, “Relative spectral power distributions of illuminants”
"""
    illuminant_A(env::SpecEnv = SPECENV)

Produces the spectrum of illuminant A, which represents tungsten-filament lighting, depending on the limits defined by SPECENV.
The relative spectral power distribution is equal to a Planckian spectrum at T ≊ 2856 K, which is also its CCT.

CCT according to the *International Practical Temperature Scale* from 1968.

This is the illuminant that should always be used when incandescent lighting is involved.

# Example
```jldoctest
julia> illuminant_A()
LSpec(Real[390.0, 391.0, 392.0, 393.0, 394.0, 395.0, 396.0, 397.0, 398.0, 399.0  …  821.0, 822.0, 823.0, 824.0, 825.0, 826.0, 827.0, 828.0, 829.0, 830.0], Real[1.0161486323473047e11, 1.0369227145593657e11, 1.05797774054285e11, 1.0793148447599527e11, 1.1009351284197664e11, 1.1228396592637178e11, 1.1450294713618913e11, 1.1675055649202544e11, 1.1902689060987242e11, 1.2133204268399858e11  …  2.17270715240179e12, 2.1757420157602646e12, 2.178757889232851e12, 2.1817547874185107e12, 2.1847327254181292e12, 2.1876917188295703e12, 2.190631783742765e12, 2.1935529367348203e12, 2.196455194865146e12, 2.199338575670583e12])
```
"""
function illuminant_A(env::SpecEnv = SPECENV)
    blackbody_illuminant(env, 1.4388 / 1.4350 * 2848)
end

# CIE Standard Illuminant B
# Source: R. G. W. Hunt: “Measuring Colour”, third edition
# Chapter 4.14, “Standard Illuminants B and C”
"""
    illuminant_B(env::SpecEnv = SPECENV)

Produces the spectrum of CIE Standard Illuminant B, with a CCT of ≈4874 K.
Originally intended to represent sunlight.
Now obsolete, but according to Hunt the European brewing industry still uses an approximation of this illuminant.

# Example
```jldoctest
julia> illuminant_B()
LSpec(Real[390.0, 391.0, 392.0, 393.0, 394.0, 395.0, 396.0, 397.0, 398.0, 399.0  …  821.0, 822.0, 823.0, 824.0, 825.0, 826.0, 827.0, 828.0, 829.0, 830.0], Real[0.313, 0.3223587343061588, 0.3319703445293299, 0.3417825875994216, 0.3517432204463422, 0.3618, 0.3719111024629742, 0.38207638112852754, 0.3923061085625939, 0.4026105573311068  …  7.557164594645486e-6, 9.446455743306857e-6, 7.557164594645484e-6, 3.778582297322749e-6, 1.734723475976807e-20, -2.267149378393646e-6, -3.0228658378581947e-6, -2.6450076081259192e-6, -1.5114329189290965e-6, -2.168404344971009e-21])
```
"""
function illuminant_B(env::SpecEnv = SPECENV)
    b = b_c_series_table
    bsp = SSpline.cubicspline(b[:,1], b[:,2], :natural)
    bin = SSpline.interp(bsp, collect(env.λmin:env.Δλ:env.λmax))
    luminance_spec(bin[1], bin[2] / 100.0)
end

# CIE Standard Illuminant C
# Source: R. G. W. Hunt: “Measuring Colour”, third edition
# Chapter 4.14, “Standard Illuminants B and C”
"""
    illuminant_C(env::SpecEnv)

Produces the spectrum of CIE Standard Illuminant C, with a CCT of ≈ 6774 K.
Originally intended to represent average daylight.

# Example
```jldoctest
julia> illuminant_C()
LSpec(Real[390.0, 391.0, 392.0, 393.0, 394.0, 395.0, 396.0, 397.0, 398.0, 399.0  …  821.0, 822.0, 823.0, 824.0, 825.0, 826.0, 827.0, 828.0, 829.0, 830.0], Real[0.474, 0.489271743700384, 0.5046749637120697, 0.5202123118735634, 0.5358864400233714, 0.5517, 0.5676560269088902, 0.583759088923221, 0.6000141374831071, 0.616426124028662  …  0.7121636730114121, 0.7159995912642652, 0.7194836730114121, 0.722791836505706, 0.7261, 0.7295488980965763, 0.7331385307954352, 0.7368337144460056, 0.7405992653977175, 0.7444])
```

How the CIE Standard Illuminant C is produced:

Standard Source A is combined with a filter made of two 1 cm thick layers, containing two solutions, C₁ and C₂.
The solutions are contained in a double cell that is made of a colorless optical glass.

Solution C₁:
- 3.412 g copper sulphate (``CuSO_4·5H_2O``)
- 3.412 g mannite [``C_6H_8(OH)_6``]
- 30.0 ml pyridine (``C_5H_5N``)
- 1000.0 ml distilled water

Solution C₂:
- 30.58 g cobalt ammonium sulphate [``CoSO_4·(NH_4)_2SO_4·6H_2O``]
- 22.52 g copper sulphate (``CuSO_4·5H_2O``)
- 10.0 ml sulphuric acid (density 1.835 g/ml)
- 1000.0 ml distilled water
"""
function illuminant_C(env::SpecEnv = SPECENV)
    c = b_c_series_table
    csp = SSpline.cubicspline(c[:,1], c[:,3], :natural)
    cin = SSpline.interp(csp, collect(env.λmin:env.Δλ:env.λmax))
    luminance_spec(cin[1], cin[2] / 100.0)
end

"""
    F_series_illuminant(type::Symbol)

Spectral power distribution of representative fluorescent lamps from 380 to 780 nm in 5 nm steps.

Available types:

`:f1` ... `:f12`

`F2`, `F7`, and `F11` should be preferred over the other distributions
as long as the choice in each of these three groups is not important or critical.
The three recommendend illuminants are also marked in the CRI lists below.

*For wavelength outside the `380 nm ... 780 nm` range the SPDs should be set to zero.*

Result: `::LSpec`

These are not CIE Standard Illuminants. They were compiled by the CIE as a representative
collection of practical fluorescent lamps. The 12 lamps in this collection can be grouped
in 3 groups:

`normal`:       two semi-broad band emmisions of antimony and manganese activations in
calcium halo-phosphate phosphor.

**CRI table**
* F1: 76
* `F2`: 64
* F3: 57
* F4: 51
* F5: 72
* F6: 59

`broad-band`:   enhanced color rendering properties in comparison to the `normal` group, mostly by using multiple phosphors. Flatter spectral distributions with a wider spectral range.

**CRI table**
* `F7:` 90
* F8: 95
* F9: 90

`three-band`:   three narrow-band emissions in the red, green, and blue spectral regions.
The narrow-band emissions are produced by ternary compositions of rare-earth phosphors.

**CRI table**
* F10: 81
* `F11`: 83
* F12: 83

# Examples
```jldoctest
julia> F_series_illuminant(:f11)
LSpec(Real[390.0, 391.0, 392.0, 393.0, 394.0, 395.0, 396.0, 397.0, 398.0, 399.0  …  821.0, 822.0, 823.0, 824.0, 825.0, 826.0, 827.0, 828.0, 829.0, 830.0], Real[0.46, 0.442, 0.424, 0.406, 0.388, 0.37, 0.554, 0.738, 0.922, 1.1059999999999999  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])

julia> F_series_illuminant(:f0)
ERROR: DomainError with f0:
F-series illuminant type does not exist.
```
"""
function F_series_illuminant(type::Symbol)
    if type in (:f1, :f2, :f3, :f4, :f5, :f6, :f7, :f8, :f9, :f10, :f11, :f12)
        n = parse(Int,string(type)[2:end])
        F_series_illuminant(SPECENV, n)
    else throw(DomainError(type, "F-series illuminant type does not exist."))
    end
end

"""
    F_series_illuminant(n::Int)

Return F series illuminant **Fn** for `n` from 1 to 12.

# Examples
```jldoctest
julia> F_series_illuminant(1)
LSpec(Real[390.0, 391.0, 392.0, 393.0, 394.0, 395.0, 396.0, 397.0, 398.0, 399.0  …  821.0, 822.0, 823.0, 824.0, 825.0, 826.0, 827.0, 828.0, 829.0, 830.0], Real[2.94, 3.046, 3.152, 3.258, 3.3640000000000003, 3.47, 3.81, 4.15, 4.49, 4.83  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])

julia> F_series_illuminant(13)
ERROR: DomainError with 13:
F series illuminant outside the range 1 ≤ n ≤ 12.
```
"""
function F_series_illuminant(env::SpecEnv, n::Int)
    if 0 < n < 13
        c = luminance_spec(f_series_table[:,1], f_series_table[:,n+1])
        csp = SSpline.spline1(c.λ, c.l)
        cin = SSpline.interp(csp, collect(env.λmin:env.Δλ:env.λmax))
        luminance_spec(cin[1], cin[2])
    else
        throw(DomainError(n, "F series illuminant outside the range 1 ≤ n ≤ 12."))
    end
end

F_series_illuminant(n::Int) = F_series_illuminant(SPECENV, n)

"""
    gas_discharge_illuminant(type::Symbol)

SPDs of common gas discharge lamps.

available types:

* `:lps`:   Low pressure sodium
* `:hps`:   High pressure sodium
* `:mb`:    High pressure mercury
* `:mbf`:   High pressure mercury with red-emitting phosphor coating, better CRI than MB
* `:mbtf`:  High pressure mercury with red-emitting phopsphor coating and tungsten filament ballast, better CRI
* `:hmi`:   High pressure mercury with metal halides. Medium arc, iodides, typical stadium lighting, was used in television to supplement daylight
* `:xenon`: Xenon, used in film projectors in cinemas, floodlighting, in lighthouses, for accelerated fading tests, flash photography, lighting fluorescent materials

# Example
```jldoctest
julia> gas_discharge_illuminant(:lps)
LSpec(Real[390.0, 391.0, 392.0, 393.0, 394.0, 395.0, 396.0, 397.0, 398.0, 399.0  …  821.0, 822.0, 823.0, 824.0, 825.0, 826.0, 827.0, 828.0, 829.0, 830.0], Real[0.0, 0.02, 0.04, 0.06, 0.08, 0.1, 0.08, 0.060000000000000005, 0.04000000000000001, 0.020000000000000004  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])

julia> gas_discharge_illuminant(:xenon)
LSpec(Real[390.0, 391.0, 392.0, 393.0, 394.0, 395.0, 396.0, 397.0, 398.0, 399.0  …  821.0, 822.0, 823.0, 824.0, 825.0, 826.0, 827.0, 828.0, 829.0, 830.0], Real[96.33, 97.176, 98.022, 98.868, 99.714, 100.56, 101.01, 101.46000000000001, 101.91, 102.36  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
```
"""
function gas_discharge_illuminant(type::Symbol)
    type == :lps  ? gas_discharge_illuminant(SPECENV, 1) :
    type == :hps  ? gas_discharge_illuminant(SPECENV, 2) :
    type == :mb   ? gas_discharge_illuminant(SPECENV, 3) :
    type == :mbf  ? gas_discharge_illuminant(SPECENV, 4) :
    type == :mbtf ? gas_discharge_illuminant(SPECENV, 5) :
    type == :hmi  ? gas_discharge_illuminant(SPECENV, 6) :
    type == :xenon ? gas_discharge_illuminant(SPECENV, 7) :
    throw(DomainError(type, "Lamp type does not exist!"))
end

gas_discharge_illuminant(n::Int) = gas_discharge_illuminant(SPECENV, n)

"""
`gas_discharge_illuminant(n::Int)`

Returns gas discharge type illuminant of index `n` for n from 1 to 7.

Result: `::LSpec`
"""
function gas_discharge_illuminant(env::SpecEnv, n::Int)
    if 0 < n < 8
        c = luminance_spec(gas_discharge_table[:,1], gas_discharge_table[:,n+1])
        csp = SSpline.linearspline(c.λ, c.l)
        cin = SSpline.interp(csp, collect(env.λmin:env.Δλ:env.λmax))
        luminance_spec(cin[1], cin[2])
    else
        throw(DomainError(n, "Gas discharge illuminant index outside the range 1 ≤ n ≤ 7."))
    end
end

"""
    sinusoidal_spd(freq::Real, modfac::Real = 1.0, ampl::Real = 1.0, phase::Real, r = (SPECENV.λmin:SPECENV.Δλ:SPECENV.λmax))

Generates a sinusoidal spectral power distribution. Useful for computational gamut and spectral calculations.
"""
function sinusoidal_spd(stype::Symbol, freq::Real, modfac::Real = 1.0, ampl::Real = 1.0, phase::Real = 0.0, r = (SPECENV.λmin:SPECENV.Δλ:SPECENV.λmax))
    if 0 ≤ ampl ≤ 1.0
        if 0 ≤ modfac ≤ 1.0
            c = 0.5 .* ampl .* (1 .+ modfac * sin.(2π * freq .* (collect(r) .- r[1]) .+ phase))
            if stype == :l
                return luminance_spec(collect(r), c)
            elseif stype == :r
                return reflectance_spec(collect(r),c)
            else
                throw(DomainError(stype, "Only :l and :r allowed."))
            end
        else
            throw(DomainError(modfac, "Expecting a value 0 ≤ modfac ≤ 1.0."))
        end
    else
        throw(DomainError(ampl, "Expecting a value 0 ≤ ampl ≤ 1.0."))
    end
end

#=
    Sources:
    George M. Hale and Marvin R. Querry, “Optical Constants of Water in the 200-nm to 200-μm Wavelength Region”, Appl. Opt. 12/1973, https://doi.org/10.1364/AO.12.000555

    W.G. Zijlstra, A. Buursma, W.P. Meeuwsen-van der Roest, “Absorption spectra of human fetal and adult oxyhemoglobin, de-oxyhemoglobin, carboxyhemoglobin, and methemoglobin”, Clin. Chem. 37 (1991) 1633/1638

    I.S. Saidi, “Transcutaneous optical measurement of hyperbilirubinemia in neonates”, PhD thesis, Rice University, Houston, TX, 1992

=#
"""
`skin_simulation()`

Skin optics Simulation
"""
function skin_simulation()

end

function μa(layer,λ)
    #=  FHB: volume fraction of hemoglobin
        FRBC: volume fratcion of erythrocytes in total volume of all blood cells
        Ht: hemotocrit
        assumption: hemoglobin is contained in the erythrocytes only

        absorption coefficients:
        μaHb

    =#
    μS = abscoeff(layer)[1]
    g = abscoeff(layer)[2]
    n = abscoeff(layer)[3]
    CBlood = abscoeff(layer)[4]
    S = abscoeff(layer)[5]
    Ht = abscoeff(layer)[6]
    FHb = abscoeff(layer)[7]
    FRBC = abscoeff(layer)[8]
    CH₂O = abscoeff(layer)[9]

    γ = FHb * FRBC * Ht

    μalayer = (1-S) * γ * Cblood * μaHb + S * γ * CBlood * μaHbO₂(λ) +
              (1- γ * CBlood) * CH₂O * μaH₂O(λ) +
              (1- γ * CBlood) * (1-CH₂O) * μaOther(λ)
    μalayer
end

function abscoeff(layer::Int)
    ABSORPTIONCOEFF[layer,:]
end

# absorption coefficient matrix for the 7 skin layers
#   μs    g     n   CBlood  S    Ht   FHb   FRBC  CH₂O
const ABSORPTIONCOEFF =
   [100  0.86  1.50  0.00  0.0  0.00  0.00  0.00  0.05;
     45  0.80  1.34  0.00  0.0  0.00  0.00  0.00  0.20;
     30  0.90  1.40  0.04  0.6  0.45  0.99  0.25  0.50;
     35  0.95  1.39  0.30  0.6  0.45  0.99  0.25  0.60;
     25  0.80  1.40  0.04  0.6  0.45  0.99  0.25  0.70;
     30  0.95  1.38  0.10  0.6  0.45  0.99  0.25  0.70;
      5  0.75  1.44  0.05  0.6  0.45  0.99  0.25  0.70]

function μaOther(λ)
    7.87e7 * λ^(-3.255)
end

function μaStratum(λ)
    (0.1 - 8.3e-4 * λ) + 0.125 * μaOther(λ)
end

function μalivingepidermis(λ)
    (0.5e10 * λ^(-3.33)) * (1 - CH₂O) + CH₂O * μaH₂O
end

function μaH₂O(λ)
    WATERCONST
    absorption(k,λ)
end

# if λ in nm, multiply by 1e7 to get absorption in cm^-1
function absorption(k,λ)
    4 * π * k / λ
end


text:
#=
    OBSOLETE - SUBSUMED BY SCORE

    evidence_ops.jl : Support for operators that process evidence and interventions
=#

export
    evidence_message,
    evidence_entry,
    intervention_message,
    intervention_entry,
    NoiseEvidence
import LinearAlgebra

function evidence_entry(hard_evidence::T, value::T) where {T}
    return hard_evidence == value ? 1.0 : 0.0
end

function evidence_entry(soft_evidence::Dict{T, Float64}, value::T) where {T}
    return get(soft_evidence, value, 0.0)
end

struct NoiseEvidence
    mean ::Union{Array{Float64,1}, Float64}
    std ::Union{Array{Float64,2}, Float64}
end

function evidence_entry(noise_evid::NoiseEvidence, value::Float64)
    val = 1.0/(noise_evid.std * sqrt(2*pi)) * exp(-0.5 * ((noise_evid.mean - value)/noise_evid.std)^2)
    return val
end

function evidence_entry(noise_evid::NoiseEvidence, value::Array{Float64,1})
    k = length(value)
    val = exp.(-0.5 * transpose(value - noise_evid.mean) * LinearAlgebra.inv(noise_evid.std) * (value - noise_evid.mean))/sqrt((2*pi)^k * LinearAlgebra.det(noise_evid.std))
    return val
end

function evidence_entry(f::Function , value) ::Float64 # return f(v), has to be Float64
    return f(value)
end

function intervention_entry(hard_evidence::T, value::T) where {T}
    return hard_evidence == value ? 1.0 : 0.0
end

function evidence_message(::SFunc{I,O,P}, range::Vector{O}, evidence::Score{O})::Score{O} where {I,O,P}
    n = length(range)
    ps = Array{Float64, 1}(undef, n)
    for i = 1:n
        ps[i] = get_score(evidence, range[i])
    end
    return SoftScore{O}(range, ps)
end

function intervention_message(SFunc, range, intervention)
    n = length(range)
    ps = Array{Float64, 1}(undef, n)
    for i = 1:n
        ps[i] = intervention_entry(intervention, range[i])
    end
    return Cat(range, ps)
end