text:
export SpiralPerturbedTrajectory, spiralPerturbedTrajectoryNodes

function SpiralPerturbedTrajectory(numProfiles, numSamplingPerProfile
                  ; TE::Real=0.0
                  , AQ::Real=1.e-3
                  , windings::Real= 6.25
                  , alpha::Real=2.0
                  , kmax::Real=0.5
                  , angleOffset::String="equispaced"
                  , kargs...)
  nodes = spiralPerturbedTrajectoryNodes(numProfiles, numSamplingPerProfile;
            windings=windings, alpha=alpha, angleOffset=angleOffset, kmax=kmax)
  times = readoutTimes(numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ)
  return  Trajectory("SpiralPerturbed", nodes, times, TE, AQ, numProfiles, numSamplingPerProfile, 1, false, true)
end

# Constructing the amplitude limited case of the variable density spiral Trajectory
function spiralPerturbedTrajectoryNodes(numProfiles::Int64
                        , numSamplingPerProfile::Int64
                        ; windings=6.25
                        , alpha=2.0
                        , kmax::Real=0.5
                        , angleOffset= "equispaced"
                        , kargs...)

  nodes = zeros(2,numSamplingPerProfile, numProfiles)
  angles = collect((0:numProfiles-1)/numProfiles) # Phase offset for multiple Profiles

  if angleOffset == "golden"
      angles = [i*(3-sqrt(5))/2  for i=0:numProfiles-1 ]
  elseif angleOffset == "random"
      angles = sort(rand(numProfiles))
  elseif angleOffset == "equispaced"
      angles = collect((0:numProfiles-1)/numProfiles)
  end

  T_ea = 1/((alpha+1)*1/(2pi*windings))
  times = range(0,stop=T_ea,length=numSamplingPerProfile)
  for l = 1:numProfiles
    for k = 1:numSamplingPerProfile
      tau = (1/(2pi*windings) *(alpha+1)*times[k])^(1/(alpha+1))
      tmp = (1+0.025*cos(k*0.1))*(tau^alpha)*(kmax-0.025)*exp(1im*2*pi*( windings*tau + angles[l] ))
      # Thresholding if the petrubed spiral goes beyond |x_j| > [0.5,0.5]
      if real(tmp) >= 0.5
          nodes[1,k,l] = 0.5 - 1e-6
      elseif real(tmp) < -0.5
          nodes[1,k,l] = -0.5
      else
          nodes[1,k,l] = real(tmp)
      end

      if imag(tmp) >= 0.5
          nodes[2,k,l] = 0.5 - 1e-6
      elseif imag(tmp) < -0.5
          nodes[2,k,l] = -0.5
      else
          nodes[2,k,l] = imag(tmp)
      end
    end
  end
  return reshape(nodes, 2, numSamplingPerProfile*numProfiles)
end


text:
# obtain the policies for different time horizon
addprocs(30);
@everywhere include("loadMod.jl");
@everywhere const GUROBI_ENV = Gurobi.Env();

Δt = 0.25;
N = 5;
iterMax = 20;

caseList = [13,33,123];

policyDict = Dict();
for ci in 1:length(caseList)
    for T in TList
        τ = Int64(1/6*T);
        for j in procs()
            remotecall_fetch(readInData,j,ci,caseList,T);
        end

        startPGT = time();
        cutDictPG = preGen(τ, T, Δt, N, iterMax);
        preGenT = time() - startPGT;

        startT = time();
        cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(τ, T, Δt, N, false, false, max(Int64(round(200/N)),20), max(Int64(round(200/N)),20),cutDictPG);
        elapsedT = time() - startT;

        policyDict[T] = [cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist,elapsedT,preGenT];
        save("policy_trained_$(ci).jld","policyDict",policyDict);
    end
end

ext: jl


text:
# Only 2D for now, 3 features from two poses
# will add multi-modality soon

# using Distributions, TransformUtils

function getUvecScaleFeature2D(wTb, bTc, wTo)
  cTo = (wTb * bTc) \ wTo
  r = cTo[1:2,3]
  α = norm(r)
  return r, r/α, α
end

function getUvecScaleBaseline2D(wTbi, wTbj, bTc)
  ciTcj = (wTbi * bTc) \ (wTbj * bTc)
  b = ciTcj[1:2,3]
  β = norm(b)
  return b, b/β, β
end

# function getEstFeatDirectii(wTbi, wTobjs::Vector, bTc)
#   xirvecs = Vector{Tuple}()
#   for wTo in wTobjs
#     push!(xirvecs, getUvecScaleFeature2D(wTbi, bTc, wTo))
#   end
#   return xirvecs
# end

function residruvec()
  nothing
end

# two sets of three feature sightings and a body to camera lever arm transform
type MultipleFeatures2D <: IncrementalInference.FunctorPairwiseMinimize
  # treat these as angles positive X->Y, X forward, Y left and Z up
  xir1::Normal
  xir2::Normal
  xir3::Normal
  xjr1::Normal
  xjr2::Normal
  xjr3a::Normal
  xjr3b::Normal
  abC::Categorical
  bTc::Array{Float64,2}
  zDim::Tuple{Int,Int,Int,Int,Int}
  # Takes five
  MultipleFeatures2D(a,b,c,d,e,f,g,h,i ) = new(a,b,c,d,e,f,g,h,i, (3,3,2,2,2) )
end

function getSample(mm2d::MultipleFeatures2D, N::Int=1)

  meas = Array{Float64,3}(3,N,2)
  meas[1,:,1] = rand(mm2d.xir1,N)
  meas[2,:,1] = rand(mm2d.xir2,N)
  meas[3,:,1] = rand(mm2d.xir3,N)
  meas[1,:,2] = rand(mm2d.xjr1,N)
  meas[2,:,2] = rand(mm2d.xjr2,N)
  meas[3,:,2] = rand(mm2d.xjr3a,N)

  select = rand(mm2d.abC,N)
  beb = select .== 2
  numb = sum(beb)
  meas[3,beb,2] = rand(mm2d.xjr3b,numb)

  return (meas, beb, numb)
end

# function (mm2d::MultipleFeatures2D)(res::Array{Float64}, idx::Int, meas::Tuple,
#           wAbi::Array{Float64,2},
#           wAbj::Array{Float64,2},
#           wAo1::Array{Float64,2},
#           wAo2::Array{Float64,2},
#           wAo3::Array{Float64,2}  )
#   #
#   wTbi, wTbj = SE2(wAbi[:,idx]), SE2(wAbj[:,idx])
#   wTo = ( SE2([wAo1[:,idx];0.0]),SE2([wAo2[:,idx];0.0]),SE2([wAo3[:,idx];0.0]) )
#
#   # b, bhat, β = getUvecScaleBaseline2D(wTbi, wTbj, mm2d.bTc)
#
#   rhat, resid = zeros(2), zeros(2)
#   # res[:] = -1e-5
#   # res[3] = -3*β
#
#   for i in 1:3
#     # sightings from first pose
#     r1, rhathat1, α1 = getUvecScaleFeature2D(wTbi, mm2d.bTc, wTo[i])
#     rhat[1] = cos(meas[1][i,idx,1])
#     rhat[2] = sin(meas[1][i,idx,1])
#     resid[1:2] = rhat-rhathat1
#     res[1] += dot(resid, resid)
#
#     # sightingings from second pose
#     r2, rhathat2, α2 = getUvecScaleFeature2D(wTbj, mm2d.bTc, wTo[i])
#     rhat[1] = cos(meas[1][i,idx,2])
#     rhat[2] = sin(meas[1][i,idx,2])
#     resid[1:2] = rhat-rhathat2
#     res[2] += dot(resid, resid)
#
#     # baseline constraint
#     # res[3] += dot(r1, bhat)
#     # res[3] += dot(r2, -bhat)
#   end
#
#   nothing
# end

# redo with angles on sightings for minimization
# NOTE -- an inefficient implementation
function (mm2d::MultipleFeatures2D)(res::Array{Float64}, idx::Int, meas::Tuple,
          wAbi::Array{Float64,2},
          wAbj::Array{Float64,2},
          wAo1::Array{Float64,2},
          wAo2::Array{Float64,2},
          wAo3::Array{Float64,2},
          wAo3b::Union{Array{Float64,2},Bool}=false  )
  #
  len3b = wAo3b==false ? 0 : size(wAo3b,2)
  dobimodal = len3b > 0
  meas[3] > 0 && !dobimodal ? error("mm2d::MultipleFeatures2D -- bi-modal mismatch") : nothing

  wTbi = SE2(wAbi[:,idx])
  wTbj = SE2(wAbj[:,idx])
  wTo1 = SE2([wAo1[:,idx];0.0])
  wTo2 = SE2([wAo2[:,idx];0.0])
  wTo3 = SE2([wAo3[:,idx];0.0])
  wTo3b = dobimodal ? SE2([wAo3b[:,idx];0.0]) : nothing

  ciTo1 = (wTbi * mm2d.bTc) \ wTo1
  ciTo2 = (wTbi * mm2d.bTc) \ wTo2
  ciTo3 = (wTbi * mm2d.bTc) \ wTo3

  cjTo1 = (wTbj * mm2d.bTc) \ wTo1
  cjTo2 = (wTbj * mm2d.bTc) \ wTo2
  cjTo3 = (wTbj * mm2d.bTc) \ wTo3
  cjTo3b = dobimodal ? (wTbj * mm2d.bTc) \ wTo3 : nothing

  the1 = atan2(se2vee(ciTo1)[2], se2vee(ciTo1)[1])
  the2 = atan2(se2vee(ciTo2)[2], se2vee(ciTo2)[1])
  the3 = atan2(se2vee(ciTo3)[2], se2vee(ciTo3)[1])

  res[1] = 0.0
  res[1] += ( meas[1][1,idx,1] - the1 )^2
  res[1] += ( meas[1][2,idx,1] - the2 )^2
  res[1] += ( meas[1][3,idx,1] - the3 )^2

  the1 = atan2(se2vee(cjTo1)[2], se2vee(cjTo1)[1])
  the2 = atan2(se2vee(cjTo2)[2], se2vee(cjTo2)[1])
  the3 = atan2(se2vee(cjTo3)[2], se2vee(cjTo3)[1])
  the3b = dobimodal ? atan2(se2vee(cjTo3b)[2], se2vee(cjTo3b)[1]) : nothing

  res[2] = 0.0
  res[2] += ( meas[1][1,idx,2] - the1 )^2
  res[2] += ( meas[1][2,idx,2] - the2 )^2

  if dobimodal
    if !meas[2][idx]
      # 3a case
      res[2] += ( meas[1][3,idx,2] - the3 )^2
    else
      # 3b case
      res[2] += ( meas[1][3,idx,2] - the3b )^2
    end
  else
    res[2] += ( meas[1][3,idx,2] - the3 )^2
  end
  sum(res)
end

#


text:
export findRelatedFromPotential

"""
    $SIGNATURES
Internal method to set which dimensions should be used as the decision variables for later numerical optimization.
"""
function _setCCWDecisionDimsConv!(ccwl::Union{CommonConvWrapper{F},
                                              CommonConvWrapper{Mixture{N_,F,S,T}}} ) where {N_,F<:Union{AbstractRelativeMinimize, AbstractPrior},S,T}
  #
  p = if ccwl.partial
    Int[ccwl.usrfnc!.partial...]
  else
    Int[1:ccwl.xDim...]
  end

  for thrid in 1:Threads.nthreads()
    length(ccwl.cpt[thrid].p) != length(p) ? resize!(ccwl.cpt[thrid].p, length(p)) : nothing
    ccwl.cpt[thrid].p .= p
  end
  nothing
end

function _setCCWDecisionDimsConv!(ccwl::Union{CommonConvWrapper{F},
                                              CommonConvWrapper{Mixture{N_,F,S,T}}} ) where {N_,F<:AbstractRelativeRoots,S,T}
  #
  for thrid in 1:Threads.nthreads()
    length(ccwl.cpt[thrid].p) != ccwl.xDim ? resize!(ccwl.cpt[thrid].p, ccwl.xDim) : nothing
    ccwl.cpt[thrid].p .= Int[1:ccwl.xDim;]
  end
  nothing
end

"""
    $(SIGNATURES)

Perform the nonlinear numerical operations to approximate the convolution with a particular user defined likelihood function (conditional), which as been prepared in the `frl` object.  This function uses root finding to enforce a non-linear function constraint.

Notes:
- remember this is a deepcopy of original sfidx, since we are generating a proposal distribution and not directly replacing the existing variable belief estimate

Future work:
- once Threads.@threads have been optmized JuliaLang/julia#19967, also see area4 branch
- improve handling of n and particleidx, especially considering future multithreading support
"""
function approxConvOnElements!( ccwl::Union{CommonConvWrapper{F},
                                            CommonConvWrapper{Mixture{N_,F,S,T}}},
                                elements::Union{Vector{Int}, UnitRange{Int}}, ::Type{MultiThreaded}  ) where {N_,F<:AbstractRelative,S,T}
  #
  Threads.@threads for n in elements
    # ccwl.thrid_ = Threads.threadid()
    ccwl.cpt[Threads.threadid()].particleidx = n

    # ccall(:jl_, Nothing, (Any,), "starting loop, thrid_=$(Threads.threadid()), partidx=$(ccwl.cpt[Threads.threadid()].particleidx)")
    _solveCCWNumeric!( ccwl )
  end
  nothing
end

function approxConvOnElements!( ccwl::Union{CommonConvWrapper{F},
                                            CommonConvWrapper{Mixture{N_,F,S,T}}},
                                elements::Union{Vector{Int}, UnitRange{Int}}, ::Type{SingleThreaded}) where {N_,F<:AbstractRelative,S,T}
  #
  for n in elements
    ccwl.cpt[Threads.threadid()].particleidx = n
    _solveCCWNumeric!( ccwl )
  end
  nothing
end

function approxConvOnElements!( ccwl::Union{CommonConvWrapper{F},
                                            CommonConvWrapper{Mixture{N_,F,S,T}}},
                                elements::Union{Vector{Int}, UnitRange{Int}} )  where {N_,F<:AbstractRelative,S,T}
  #
  approxConvOnElements!(ccwl, elements, ccwl.threadmodel)
end

"""
    $(SIGNATURES)

Prepare a common functor computation object `prepareCommonConvWrapper{T}` containing
the user factor functor along with additional variables and information using during
approximate convolution computations.
"""
function prepareCommonConvWrapper!( F_::Type{<:AbstractRelative},
                                    ccwl::CommonConvWrapper{F},
                                    Xi::AbstractVector{<:DFGVariable},
                                    solvefor::Symbol,
                                    N::Int;
                                    needFreshMeasurements::Bool=true,
                                    solveKey::Symbol=:default  ) where {F <: FunctorInferenceType}
  #

  # FIXME, order of fmd ccwl cf are a little weird and should be revised.
  ARR = Array{Array{Float64,2},1}()
  # FIXME maxlen should parrot N (barring multi-/nullhypo issues)
  maxlen, sfidx, manis = prepareparamsarray!(ARR, Xi, solvefor, N, solveKey=solveKey)
  # TODO should be selecting for the correct multihypothesis mode
  ccwl.params = ARR
  # get factor metadata -- TODO, populate, also see #784
  fmd = FactorMetadata(Xi, getLabel.(Xi), ARR, solvefor, nothing)

  # TODO consolidate with ccwl??
  # FIXME do not divert Mixture for sampling
  # cf = _buildCalcFactorMixture(ccwl, fmd, 1, ccwl.measurement, ARR) # TODO perhaps 0 is safer
  cf = CalcFactor( ccwl.usrfnc!, fmd, 0, length(ccwl.measurement), ccwl.measurement, ARR)

  #  get variable node data
  vnds = Xi

  # option to disable fresh samples
  if needFreshMeasurements
    # TODO refactor
    ccwl.measurement = sampleFactor(cf, maxlen)
    # sampleFactor!(ccwl, maxlen, fmd, vnds)
  end

  if ccwl.specialzDim
    ccwl.zDim = ccwl.usrfnc!.zDim[sfidx]
  else
    ccwl.zDim = size(ccwl.measurement[1],1) # TODO -- zDim aspect needs to be reviewed
  end
  ccwl.varidx = sfidx

  ccwl.xDim = size(ccwl.params[sfidx],1)
  # ccwl.xDim = size(ccwl.cpt[1].X,1)
  # info("what? sfidx=$(sfidx), ccwl.xDim = size(ccwl.params[sfidx]) = $(ccwl.xDim), size=$(size(ccwl.params[sfidx]))")
  for thrid in 1:Threads.nthreads()
    cpt_ = ccwl.cpt[thrid]
    cpt_.X = ccwl.params[sfidx]
    resize!(cpt_.p, ccwl.xDim)
    cpt_.p .= 1:ccwl.xDim
    # used in ccw functor for AbstractRelativeMinimize
    # TODO JT - Confirm it should be updated here. Testing in prepgenericconvolution
    resize!(cpt_.res, ccwl.zDim)
    fill!(cpt_.res, 0.0)
    # cpt_.res = zeros(ccwl.xDim)
  end

  return sfidx, maxlen, manis
end

function prepareCommonConvWrapper!( ccwl::Union{CommonConvWrapper{F},
                                                CommonConvWrapper{Mixture{N_,F,S,T}}},
                                    Xi::AbstractVector{<:DFGVariable},
                                    solvefor::Symbol,
                                    N::Int;
                                    needFreshMeasurements::Bool=true,
                                    solveKey::Symbol=:default  ) where {N_,F<:AbstractRelative,S,T}
  #
  prepareCommonConvWrapper!(F, ccwl, Xi, solvefor, N, needFreshMeasurements=needFreshMeasurements, solveKey=solveKey)
end

function generateNullhypoEntropy( val::AbstractMatrix{<:Real},
                                  maxlen::Int,
                                  spreadfactor::Real=10  )
  #
  # covD = sqrt.(vec(Statistics.var(val,dims=2))) .+ 1e-3
  # cVar = diagm((spreadfactor*covD).^2)
  len = size(val, 1)
  cVar = diagm((spreadfactor*ones(len)).^2)
  mu = zeros( len )
  MvNormal( mu, cVar )
end

function calcVariableCovarianceBasic(arr::AbstractMatrix)
  # cannot calculate the stdev from uninitialized state
  msst = Statistics.std(arr, dims=2)
  # FIXME use adaptive scale, see #802
  msst_ = 0 < sum(1e-10 .< msst) ? maximum(msst) : 1.0
  return msst_
end

"""
    $SIGNATURES

Control the amount of entropy to add to null-hypothesis in multihypo case.

Notes:
- FIXME, Currently only supports Euclidean domains.
"""
function calcVariableDistanceExpectedFractional(ccwl::CommonConvWrapper,
                                                sfidx::Int,
                                                certainidx::Vector{Int};
                                                kappa::Float64=3.0  )
  #
  if sfidx in certainidx
    msst_ = calcVariableCovarianceBasic(ccwl.params[sfidx])
    return kappa*msst_
  end
  # @assert !(sfidx in certainidx) "null hypo distance does not work for sfidx in certainidx"

  # get mean of all fractional variables
  uncertainidx = setdiff(1:length(ccwl.params), certainidx)
  uncMeans = zeros(size(ccwl.params[sfidx],1), length(uncertainidx))
  dists = zeros(length(uncertainidx)+length(certainidx))
  dims = size(ccwl.params[sfidx],1)
  count = 0
  for i in uncertainidx
    count += 1
    uncMeans[:,count] = Statistics.mean(ccwl.params[i], dims=2)[:]
  end
  count = 0
  refMean = Statistics.mean(ccwl.params[sfidx], dims=2)[:]
  for i in uncertainidx
    count += 1
    dists[count] = norm(refMean - uncMeans[:,count])
  end
  # also check distance to certainidx for general scale reference (workaround heuristic)
  for cidx in certainidx
    count += 1
    cerMean = Statistics.mean(ccwl.params[cidx], dims=2)[:]
    dists[count] = norm(refMean[1:dims] - cerMean[1:dims])
  end

  push!(dists, 1e-2)
  return kappa*maximum(dists)
end

function addEntropyOnManifoldHack!(addEntr::Union{AbstractMatrix{<:Real},SubArray}, maniAddOps, spreadDist::Real)
  # add 1σ "noise" level to max distance as control
  for i in 1:size(addEntr, 1)
    for j in 1:size(addEntr,2)
      addEntr[i,j] = maniAddOps[i](addEntr[i,j], spreadDist*(rand()-0.5))
    end
  end
  nothing
end

"""
    $(SIGNATURES)

Common function to compute across a single user defined multi-hypothesis ambiguity per factor.  This function dispatches both `AbstractRelativeRoots` and `AbstractRelativeMinimize` factors.
"""
function computeAcrossHypothesis!(ccwl::Union{CommonConvWrapper{F},
                                              CommonConvWrapper{Mixture{N_,F,S,T}}},
                                  allelements::AbstractVector,
                                  activehypo,
                                  certainidx::Vector{Int},
                                  sfidx::Int,
                                  maxlen::Int,
                                  maniAddOps::Tuple;
                                  spreadNH::Real=3.0 ) where {N_,F<:AbstractRelative,S,T}
  #
  count = 0

  # setup the partial or complete decision variable dimensions for this ccwl object
  # NOTE perhaps deconv has changed the decision variable list, so placed here during consolidation phase
  _setCCWDecisionDimsConv!(ccwl)

  # TODO remove assert once all GenericWrapParam has been removed
  # @assert norm(ccwl.certainhypo - certainidx) < 1e-6
  for (hypoidx, vars) in activehypo
    count += 1
    if sfidx in certainidx && hypoidx != 0 || hypoidx in certainidx || hypoidx == sfidx
      # hypo case hypoidx, sfidx = $hypoidx, $sfidx
      for i in 1:Threads.nthreads()  ccwl.cpt[i].activehypo = vars; end
      approxConvOnElements!(ccwl, allelements[count])
    elseif hypoidx != sfidx && hypoidx != 0  # sfidx in uncertnidx
      # multihypo, take other value case
      # sfidx=2, hypoidx=3:  2 should take a value from 3
      # sfidx=3, hypoidx=2:  3 should take a value from 2
        # DEBUG sfidx=2, hypoidx=1 -- bad when do something like multihypo=[0.5;0.5] -- issue 424
      ccwl.params[sfidx][:,allelements[count]] = view(ccwl.params[hypoidx],:,allelements[count])
    elseif hypoidx == 0
      # basically do nothing since the factor is not active for these allelements[count]
      # inject lots of entropy in nullhypo case
      addEntr = view(ccwl.params[sfidx], :, allelements[count])
      # make spread (1σ) equal to mean distance of other fractionals
      spreadDist = calcVariableDistanceExpectedFractional(ccwl, sfidx, certainidx, kappa=spreadNH)
      # ENT = generateNullhypoEntropy(addEntr, maxlen, spreadDist) # TODO
      # add 1σ "noise" level to max distance as control
      addEntropyOnManifoldHack!(addEntr, maniAddOps, spreadDist)
      # for i in 1:size(addEntr, 1)
      #   for j in 1:size(addEntr,2)
      #     # FIXME, should be with ENT
      #     addEntr[i,j] = maniAddOps[i](addEntr[i,j], spreadDist*(rand()-0.5))
      #   end
      # end
    else
      error("computeAcrossHypothesis -- not dealing with multi-hypothesis case correctly")
    end
  end
  nothing
end
  # elseif hypoidx == sfidx
  #   # multihypo, do conv case, hypoidx == sfidx
  #   ah = sort(union([sfidx;], certainidx))
  #   @assert norm(ah - vars) < 1e-10
  #   for i in 1:Threads.nthreads()  ccwl.cpt[i].activehypo = ah; end
  #   approxConvOnElements!(ccwl, allelements[count])

"""
    $(SIGNATURES)

Multiple dispatch wrapper for `<:AbstractRelativeRoots` types, to prepare and execute the general approximate convolution with user defined factor residual functions.  This method also supports multihypothesis operations as one mechanism to introduce new modality into the proposal beliefs.

Planned changes will fold null hypothesis in as a standard feature and no longer appear as a separate `InferenceType`.
"""
function evalPotentialSpecific( Xi::AbstractVector{<:DFGVariable},
                                ccwl::CommonConvWrapper{T},
                                solvefor::Symbol,
                                T_::Type{<:AbstractRelative},
                                measurement::Tuple=(zeros(0,100),);
                                needFreshMeasurements::Bool=true,
                                solveKey::Symbol=:default,
                                N::Int=size(measurement[1],2),
                                spreadNH::Real=3.0,
                                dbg::Bool=false  ) where {T <: AbstractFactor}
  #

  # Prep computation variables
  # FIXME #1025, should FMD be built here?
  sfidx, maxlen, manis = prepareCommonConvWrapper!(ccwl, Xi, solvefor, N, needFreshMeasurements=needFreshMeasurements, solveKey=solveKey)
  # check for user desired measurement values
  if 0 < size(measurement[1],1)
    ccwl.measurement = measurement
  end

  # Check which variables have been initialized
  isinit = map(x->isInitialized(x), Xi)

  # get manifold add operations
  # TODO, make better use of dispatch, see JuliaRobotics/RoME.jl#244
  addOps, d1, d2, d3 = buildHybridManifoldCallbacks(manis)

  # assemble how hypotheses should be computed
  _, allelements, activehypo, mhidx = assembleHypothesesElements!(ccwl.hypotheses, maxlen, sfidx, length(Xi), isinit, ccwl.nullhypo )
  certainidx = ccwl.certainhypo

  # perform the numeric solutions on the indicated elements
  # error("ccwl.xDim=$(ccwl.xDim)")
  computeAcrossHypothesis!(ccwl, allelements, activehypo, certainidx, sfidx, maxlen, addOps, spreadNH=spreadNH)

  return ccwl.params[ccwl.varidx]
end

# TODO `measurement` might not be properly wired up yet
function evalPotentialSpecific( Xi::AbstractVector{<:DFGVariable},
                                ccwl::CommonConvWrapper{T},
                                solvefor::Symbol,
                                T_::Type{<:AbstractPrior},
                                measurement::Tuple=(zeros(0,0),);
                                needFreshMeasurements::Bool=true,
                                solveKey::Symbol=:default,
                                N::Int=size(measurement[1],2),
                                dbg::Bool=false,
                                spreadNH::Real=3.0 ) where {T <: AbstractFactor}
  #
  # FIXME, NEEDS TO BE CLEANED UP AND WORK ON MANIFOLDS PROPER
  fnc = ccwl.usrfnc!
  sfidx = 1
  oldVal = getVal(Xi[sfidx], solveKey=solveKey)
  nn = maximum([N; size(measurement[1],2); size(oldVal,2); size(ccwl.params[sfidx],2)]) # (N <= 0 ? size(getVal(Xi[1]),2) : N)
  vnds = Xi # (x->getSolverData(x)).(Xi)
  # FIXME better standardize in-place operations (considering solveKey)
  if needFreshMeasurements
    cf = CalcFactor( ccwl.usrfnc!, _getFMdThread(ccwl), 0, length(ccwl.measurement), ccwl.measurement, ccwl.params)
    ccwl.measurement = sampleFactor(cf, nn)
    # fmd = FactorMetadata(Xi, getLabel.(Xi), ccwl.params, solvefor, nothing)
    # sampleFactor!(ccwl, nn, fmd, vnds)
  end
  # Check which variables have been initialized
  isinit = map(x->isInitialized(x), Xi)
  _, allelements, activehypo, mhidx = assembleHypothesesElements!(ccwl.hypotheses, nn, sfidx, length(Xi), isinit, ccwl.nullhypo )
  # get solvefor manifolds
  manis = getManifolds(Xi[sfidx])
  addOps, d1, d2, d3 = buildHybridManifoldCallbacks(manis)
  # two cases on how to use the measurement
  nhmask = mhidx .== 0
  ahmask = mhidx .== 1
  # generate nullhypo samples
  # inject lots of entropy in nullhypo case
  # make spread (1σ) equal to mean distance of other fractionals
  # FIXME better standardize in-place operations (considering solveKey)
  addEntr = if size(oldVal,2) == nn
    deepcopy(oldVal)  #ccwl.params[sfidx])
  else
    ret = zeros(size(oldVal,1),nn)
    ret[:,1:size(oldVal,2)] .= oldVal #ccwl.params[sfidx]
    ret
  end
  # @show nn, size(addEntr), size(nhmask), size(oldVal)
  addEntrNH = view(addEntr, :, nhmask)
  spreadDist = spreadNH*calcVariableCovarianceBasic(addEntr)
  # ENT = generateNullhypoEntropy(addEntr, nn, spreadDist)
  if !ccwl.partial
      addEntr[:,ahmask] = ccwl.measurement[1][:,ahmask]
      # ongoing part of RoME.jl #244
      addEntropyOnManifoldHack!(addEntrNH, addOps, spreadDist)
    # return ccwl.measurement[1]
  else
    i = 0
    for dimnum in fnc.partial
      i += 1
      addEntr[dimnum,ahmask] = ccwl.measurement[1][i,ahmask]
      addEntrNHp = view(addEntr, dimnum, nhmask)
      # ongoing part of RoME.jl #244
      addEntropyOnManifoldHack!(addEntrNHp, addOps[dimnum:dimnum], spreadDist)
    end
  end
  return addEntr
end

function evalPotentialSpecific( Xi::AbstractVector{<:DFGVariable},
                                ccwl::CommonConvWrapper{Mixture{N_,F,S,T}},
                                solvefor::Symbol,
                                measurement::Tuple=(zeros(0,0),);
                                needFreshMeasurements::Bool=true,
                                solveKey::Symbol=:default,
                                N::Int=size(measurement[1],2),
                                dbg::Bool=false,
                                spreadNH::Real=3.0 ) where {N_,F<:FunctorInferenceType,S,T}
  #
  evalPotentialSpecific(Xi,
                        ccwl,
                        solvefor,
                        F,
                        measurement;
                        needFreshMeasurements=needFreshMeasurements,
                        solveKey=solveKey,
                        N=N,
                        dbg=dbg,
                        spreadNH=spreadNH )
end

function evalPotentialSpecific( Xi::AbstractVector{<:DFGVariable},
                                ccwl::CommonConvWrapper{F},
                                solvefor::Symbol,
                                measurement::Tuple=(zeros(0,0),);
                                needFreshMeasurements::Bool=true,
                                solveKey::Symbol=:default,
                                N::Int=size(measurement[1],2),
                                dbg::Bool=false,
                                spreadNH::Real=3.0 ) where {F <: FunctorInferenceType}
  #
  evalPotentialSpecific(Xi,
                        ccwl,
                        solvefor,
                        F,
                        measurement;
                        needFreshMeasurements=needFreshMeasurements,
                        solveKey=solveKey,
                        N=N,
                        dbg=dbg,
                        spreadNH=spreadNH )
end

"""
    $(SIGNATURES)

Single entry point for evaluating factors from factor graph, using multiple dispatch to locate the correct `evalPotentialSpecific` function.
"""
function evalFactor(dfg::AbstractDFG,
                    fct::DFGFactor,
                    solvefor::Symbol,
                    measurement::Tuple=(zeros(0,100),);
                    needFreshMeasurements::Bool=true,
                    solveKey::Symbol=:default,
                    N::Int=size(measurement[1],2),
                    dbg::Bool=false  )
  #

  ccw = _getCCW(fct)
  # TODO -- this build up of Xi is excessive and could happen at addFactor time
  variablelist = getVariableOrder(fct)
  Xi = getVariable.(dfg, variablelist)

  # setup operational values before compute (likely to be refactored)
  for i in 1:Threads.nthreads()
    ccw.cpt[i].factormetadata.variablelist = variablelist
    ccw.cpt[i].factormetadata.solvefor = solvefor
  end

  return evalPotentialSpecific( Xi, ccw, solvefor, measurement, needFreshMeasurements=needFreshMeasurements,
                                solveKey=solveKey, N=N, dbg=dbg, spreadNH=getSolverParams(dfg).spreadNH )
  #
end

function approxConv(dfg::AbstractDFG,
                    fc::DFGFactor,
                    target::Symbol,
                    measurement::Tuple=(zeros(0,0),);
                    solveKey::Symbol=:default,
                    N::Int=size(measurement[1],2) )
  #
  v1 = getVariable(dfg, target)
  N = N == 0 ? getNumPts(v1) : N
  return evalFactor(dfg, fc, v1.label, measurement, solveKey=solveKey, N=N)
end

"""
    $SIGNATURES

Calculate the sequential series of convolutions in order as listed by `fctLabels`, and starting from the
value already contained in the first variable.

Notes
- The ultimate `target` variable must be given to allow path discovery through n-ary factors.
- Fresh starting point will be used if first element in `fctLabels` is a unary `<:AbstractPrior`.
- This function will not change any values in `dfg`, and might have slightly less speed performance to meet this requirement.
- pass in `tfg` to get a recoverable result of all convolutions in the chain.
- `setPPE` and `setPPEmethod` can be used to store PPE information in temporary `tfg`

DevNotes
- TODO strong requirement that this function is super efficient on single factor/variable case!
- FIXME must consolidate with `accumulateFactorMeans`
- TODO `solveKey` not fully wired up everywhere yet
  - tfg gets all the solveKeys inside the source `dfg` variables
- TODO add a approxConv on PPE option
  - Consolidate with [`accumulateFactorMeans`](@ref), `approxConvBinary`

Related

[`approxDeconv`](@ref), `LightDFG.findShortestPathDijkstra`, [`evalFactor`](@ref)
"""
function approxConv(dfg::AbstractDFG,
                    from::Symbol,
                    target::Symbol,
                    measurement::Tuple=(zeros(0,0),);
                    solveKey::Symbol=:default,
                    N::Int = size(measurement[1],2),
                    tfg::AbstractDFG = initfg(),
                    setPPEmethod::Union{Nothing, Type{<:AbstractPointParametricEst}}=nothing,
                    setPPE::Bool= setPPEmethod !== nothing,
                    path::AbstractVector{Symbol}=Symbol[]  )
  #
  @assert isVariable(dfg, target) "approxConv(dfg, from, target,...) where `target`=$target must be a variable in `dfg`"

  if from in ls(dfg, target)
    # direct request
    # TODO avoid this allocation for direct cases ( dfg, :x1x2f1, :x2[/:x1] )
    path = Symbol[from; target]
    varLbls = Symbol[target;]
  else
    # must first discover shortest factor path in dfg
    # TODO DFG only supports LightDFG.findShortestPathDijkstra at the time of writing (DFG v0.10.9)
    path = 0 == length(path) ? findShortestPathDijkstra(dfg, from, target) : path
    @assert path[1] == from "sanity check that shortest path function is working as expected"

    # list of variables
    fctMsk = isFactor.(dfg, path)
    # which factors in the path
    fctLbls = path[fctMsk]
    # must still add
    varLbls =  union(lsf.(dfg, fctLbls)...)
    neMsk = exists.(tfg, varLbls) .|> x-> xor(x,true)
    # put the non-existing variables into the temporary graph `tfg`
    # bring all the solveKeys too
    addVariable!.(tfg, getVariable.(dfg, varLbls[neMsk]))
    # variables adjacent to the shortest path should be initialized from dfg
    setdiff(varLbls, path[xor.(fctMsk,true)]) .|> x->initManual!(tfg, x, getBelief(dfg, x))
  end

  # find/set the starting point
  idxS = 1
  pts = if varLbls[1] == from
    # starting from a variable
    pts0 = getBelief(dfg, varLbls[1]) |> getPoints
  else
    # chain would start one later
    idxS += 1
    # get the factor
    fct0 = getFactor(dfg,from)
    # get the Matrix{<:Real} of projected points
    pts1 = approxConv(dfg, fct0, path[2], measurement, solveKey=solveKey, N=N)
    length(path) == 2 ? (return pts1) : pts1
  end
  # didn't return early so shift focus to using `tfg` more intensely
  initManual!(tfg, varLbls[1], pts)
  # use in combination with setPPE and setPPEmethod keyword arguments
  ppemethod = setPPEmethod === nothing ? MeanMaxPPE : setPPEmethod
  !setPPE ? nothing : setPPE!(tfg, varLbls[1], solveKey, ppemethod)

  # do chain of convolutions
  for idx in idxS:length(path)
    if fctMsk[idx]
      # this is a factor path[idx]
      fct = getFactor(dfg, path[idx])
      addFactor!(tfg, fct)
      pts = approxConv(tfg, fct, path[idx+1], solveKey=solveKey, N=N)
      initManual!(tfg, path[idx+1], pts)
      !setPPE ? nothing : setPPE!(tfg, path[idx+1], solveKey, ppemethod)
    end
  end

  # return target variable values
  return getBelief(tfg, target) |> getPoints
end

## ====================================================================================
## TODO better consolidate below with existing functions
## ====================================================================================

# TODO should this be consolidated with regular approxConv?
# TODO, perhaps pass Xi::Vector{DFGVariable} instead?
function approxConvBinary(arr::Array{Float64,2},
                          meas::FunctorInferenceType,
                          outdims::Int,
                          fmd::FactorMetadata,
                          measurement::Tuple=(zeros(0,size(arr,2)),);
                          varidx::Int=2,
                          N::Int=size(arr,2),
                          vnds=DFGVariable[] )
  #
  # N = N == 0 ? size(arr,2) : N
  pts = zeros(outdims,N);
  t = Array{Array{Float64,2},1}()
  push!(t,arr)
  push!(t,pts)

  fmd.arrRef = t

  # TODO consolidate with ccwl??
  # FIXME do not divert Mixture for sampling
  # cf = _buildCalcFactorMixture(ccwl, fmd, 1, ccwl.measurement, ARR) # TODO perhaps 0 is safer
  # FIXME 0, 0, ()
  cf = CalcFactor( meas, fmd, 0, 0, (), t)

  measurement = size(measurement[1],2) == 0 ? sampleFactor(cf, N) : measurement
  # measurement = size(measurement[1],2) == 0 ? sampleFactor(meas, N, fmd, vnds) : measurement

  zDim = size(measurement[1],1)
  ccw = CommonConvWrapper(meas, t[varidx], zDim, t, fmd, varidx=varidx, measurement=measurement)  # N=> size(measurement[1],2)

  for n in 1:N
    ccw.cpt[Threads.threadid()].particleidx = n
    _solveCCWNumeric!( ccw )
  end
  return pts
end

"""
    $SIGNATURES

Calculate both measured and predicted relative variable values, starting with `from` at zeros up to `to::Symbol`.

Notes
- assume single variable separators only.
"""
function accumulateFactorChain( dfg::AbstractDFG,
                                from::Symbol,
                                to::Symbol,
                                fsyms::Vector{Symbol}=findFactorsBetweenNaive(dfg, from, to);
                                initval=zeros(size(getVal(dfg, from))))

  # get associated variables
  svars = union(ls.(dfg, fsyms)...)

  # use subgraph copys to do calculations
  tfg_meas = buildSubgraph(dfg, [svars;fsyms])
  tfg_pred = buildSubgraph(dfg, [svars;fsyms])

  # drive variable values manually to ensure no additional stochastics are introduced.
  nextvar = from
  initManual!(tfg_meas, nextvar, initval)
  initManual!(tfg_pred, nextvar, initval)

  # nextfct = fsyms[1] # for debugging
  for nextfct in fsyms
    nextvars = setdiff(ls(tfg_meas,nextfct),[nextvar])
    @assert length(nextvars) == 1 "accumulateFactorChain requires each factor pair to separated by a single variable"
    nextvar = nextvars[1]
    meas, pred = approxDeconv(dfg, nextfct) # solveFactorMeasurements
    pts_meas = approxConv(tfg_meas, nextfct, nextvar, (meas,ones(Int,100),collect(1:100)))
    pts_pred = approxConv(tfg_pred, nextfct, nextvar, (pred,ones(Int,100),collect(1:100)))
    initManual!(tfg_meas, nextvar, pts_meas)
    initManual!(tfg_pred, nextvar, pts_pred)
  end
  return getVal(tfg_meas,nextvar), getVal(tfg_pred,nextvar)
end

"""
    $(SIGNATURES)

Compute proposal belief on `vertid` through `fct` representing some constraint in factor graph.
Always full dimension variable node -- partial constraints will only influence subset of variable dimensions.
The remaining dimensions will keep pre-existing variable values.

Notes
- fulldim is true when "rank-deficient" -- TODO swap to false (or even float)
"""
function findRelatedFromPotential(dfg::AbstractDFG,
                                  fct::DFGFactor,
                                  target::Symbol,
                                  measurement::Tuple=(zeros(0,0),);
                                  N::Int=size(measurement[1],2),
                                  solveKey::Symbol=:default,
                                  dbg::Bool=false  )
  #

  # # assuming it is properly initialized TODO
  pts = evalFactor(dfg, fct, target, solveKey=solveKey, N=N, dbg=dbg);
  # pts = approxConv(dfg, fct, target, measurement, N=N, solveKey=solveKey)

  # # determine if evaluation is "dimension-deficient"
  # solvable dimension
  inferdim = getFactorSolvableDim(dfg, fct, target)
  # zdim = getFactorDim(fct)
  # vdim = getVariableDim(DFG.getVariable(dfg, target))

  # TODO -- better to upsample before the projection
  Ndim = size(pts,1)
  Npoints = size(pts,2)
  # Assume we only have large particle population sizes, thanks to addNode!
  manis = getManifolds(dfg, target)
  # manis = getSofttype(DFG.getVariable(dfg, target)).manifolds # older
  proposal = AMP.manikde!(pts, manis)

  # FIXME consolidate with approxConv method instead
  if Npoints != N # this is where we control the overall particle set size
      proposal = resample(proposal,N)
  end
  return (proposal, inferdim)
end

"""
    $SIGNATURES

Compute the proposals of a destination vertex for each of `factors` and place the result
as belief estimates in both `dens` and `partials` respectively.

Notes
- TODO: also return if proposals were "dimension-deficient" (aka ~rank-deficient).
"""
function proposalbeliefs!(dfg::AbstractDFG,
                          destvertlabel::Symbol,
                          factors::AbstractVector{<:DFGFactor},
                          dens::Vector{BallTreeDensity},
                          partials::Dict{Int, Vector{BallTreeDensity}},
                          measurement::Tuple=(zeros(0,0),);
                          solveKey::Symbol=:default,
                          N::Int=100,
                          dbg::Bool=false  )
  #
  inferddimproposal = Vector{Float64}(undef, length(factors))
  for (count,fct) in enumerate(factors)
    data = getSolverData(fct)
    p, inferd = findRelatedFromPotential(dfg, fct, destvertlabel, measurement, N=N, dbg=dbg, solveKey=solveKey)
    if _getCCW(data).partial   # partial density
      pardims = _getCCW(data).usrfnc!.partial
      for dimnum in pardims
        if haskey(partials, dimnum)
          push!(partials[dimnum], marginal(p,[dimnum]))
        else
          partials[dimnum] = BallTreeDensity[marginal(p,[dimnum])]
        end
      end
    else # add onto full density list
      push!(dens, p)
    end
    inferddimproposal[count] = inferd
  end
  inferddimproposal
end

#


text:
using ModelingToolkitStandardLibrary.Thermal, ModelingToolkit, OrdinaryDiffEq, Test
@parameters t

# Test HeatCapacitor, TemperatureSensor, RelativeTemperatureSensor, FixedTemperature
@testset "Heat systems" begin
    T, C, G = 10.0, 10.0, 10.0
    @variables final_T(t)
    @named mass1         = HeatCapacitor(C=C)
    @named mass2         = HeatCapacitor(C=C)
    @named th_conductor  = ThermalConductor(G=G)
    @named reltem_sensor = RelativeTemperatureSensor()
    @named T_sensor1     = TemperatureSensor()
    @named T_sensor2     = TemperatureSensor()
    @named tem_src       = FixedTemperature(T=T)

    @info "Building a single-body system..."
    eqs = [
        connect(mass1.hp, th_conductor.hp1)
        connect(th_conductor.hp2, reltem_sensor.hp1)
        connect(reltem_sensor.hp2, tem_src.hp)
    ]
    @named h1 = ODESystem(eqs, t, systems=[mass1, reltem_sensor, tem_src, th_conductor])
    sys = structural_simplify(h1)

    u0 = [
        mass1.T        => 2.0
        th_conductor.T => 10.0
    ]
    prob = ODEProblem(sys, u0, (0, 2.0))
    sol = solve(prob, Rosenbrock23())

    temperatures = reduce(hcat, sol.u)
    # Check if Relative temperature sensor reads the temperature of heat capacitor
    # when connected to a thermal conductor and a fixed temperature source
    @test sol[reltem_sensor.T] == temperatures[1, :] - temperatures[2, :] - sol[tem_src.hp.T]

    @info "Building a two-body system..."
    eqs = [
        connect(T_sensor1.hp, mass1.hp, th_conductor.hp1)
        connect(th_conductor.hp2, mass2.hp, T_sensor2.hp)
        final_T ~ (mass1.C * mass1.T + mass2.C * mass2.T) /
        (mass1.C + mass2.C)
    ]
    @named h2 = ODESystem(eqs, t, [final_T], [],
                        systems=[mass1, mass2, T_sensor1, T_sensor2, th_conductor])
    sys = structural_simplify(h2)

    u0 = [
        mass1.T => 1.0
        mass2.T => 10.0
        final_T => 12
    ]
    prob = ODEProblem(sys, u0, (0, 3.0))
    sol = solve(prob, Rosenbrock23())

    m1, m2 = sol.u[end]
    @test m1 ≈ m2 atol=1e-1
    mass_T = reduce(hcat, sol.u)
    @test sol[T_sensor1.T] == mass_T[1, :]
    @test sol[T_sensor2.T] == mass_T[2, :]
end

# Test HeatFlowSensor, FixedHeatFlow, ThermalResistor, ThermalConductor, ThermalGround
@testset "Heat flow system" begin
    C, G, R = 10, 10, 10
    @named flow_src     = FixedHeatFlow(Q_flow=50, α=100)
    @named mass1        = HeatCapacitor(C=C)
    @named hf_sensor1   = HeatFlowSensor()
    @named hf_sensor2   = HeatFlowSensor()
    @named th_conductor = ThermalConductor(G=G)
    @named th_resistor  = ThermalResistor(R=R)
    @named th_ground    = ThermalGround()

    @info "Building a heat-flow system..."
    eqs = [
        connect(mass1.hp, th_resistor.hp1, th_conductor.hp1)
        connect(th_conductor.hp2, flow_src.hp, hf_sensor1.hp1, hf_sensor2.hp1)
        connect(th_resistor.hp2, hf_sensor1.hp2, hf_sensor2.hp2, th_ground.hp)
    ]
    @named h2 = ODESystem(eqs, t,
                          systems=[mass1, hf_sensor1, hf_sensor2,
                                   th_resistor, flow_src, th_ground, th_conductor])
    sys = structural_simplify(h2)

    u0 = [
        mass1.T            => 10.0
        th_resistor.Q_flow => 1.0
    ]
    prob = ODEProblem(sys, u0, (0, 3.0))
    sol  = solve(prob, Rosenbrock23())

    @test sol[th_conductor.T]*G == sol[th_conductor.Q_flow]
    @test sol[th_conductor.Q_flow] ≈ sol[hf_sensor1.Q_flow] + sol[flow_src.hp.Q_flow]

    @test sol[mass1.T] == sol[th_resistor.T]
    @test sol[th_resistor.T]./R ≈ sol[th_resistor.Q_flow]

end

# Tests ConvectiveResistor and includes FixedTemperature and ThermalResistor
@testset "Piston cylinder wall" begin
    Tᵧ, Tᵪ = 1000, 10 # ᵧ -> gas and ᵪ -> coolant
    Rᵧ, Rᵪ = 50e-4, 10e-4 # R = 1/h; h is convection co-efficient
    R_wall = 1.5e-4
    @named coolant     = ConvectiveResistor(R=Rᵪ)
    @named gas         = ConvectiveResistor(R=Rᵧ)
    @named wall        = ThermalResistor(R=R_wall)
    @named gas_tem     = FixedTemperature(T=Tᵧ)
    @named coolant_tem = FixedTemperature(T=Tᵪ)

    @info "Building a piston-cylinder..."
    eqs = [
        connect(gas_tem.hp, gas.solidport)
        connect(gas.fluidport, wall.hp1)
        connect(wall.hp2, coolant.fluidport)
        connect(coolant.solidport, coolant_tem.hp)
    ]
    @named piston = ODESystem(eqs, t, systems=[gas_tem, wall, gas, coolant, coolant_tem])
    sys = structural_simplify(piston)

    u0 = [
        coolant.dT  => 5.0
        wall.Q_flow => 10.0
    ]
    prob = ODEProblem(sys, u0, (0, 3.0))
    sol = solve(prob, Rosenbrock23())

    # Heat-flow-rate is equal in magnitude
    # and opposite in direction
    @test sol[gas.Q_flow] + sol[coolant.Q_flow] == zeros(length(sol))
end

# Test ConvectiveConductor, BodyRadiation
@testset "Radiator system" begin
    Tᵧ, Tᵪ = 1000, 10 # ᵧ -> gas and ᵪ -> coolant
    R_wall = 10
    G = 0.04
    σ   = 5.6703744191844294e-8 # Stefan-Boltzmann constant

    @named base        = ThermalResistor(R=R_wall)
    @named gas_tem     = FixedTemperature(T=Tᵧ)
    @named coolant_tem = FixedTemperature(T=Tᵪ)
    @named radiator    = BodyRadiation(G=G)
    @named ground      = ThermalGround()
    @named dissipator  = ConvectiveConductor(G=10)

    @info "Building a radiator..."
    eqs = [
        connect(gas_tem.hp, radiator.hp1, base.hp1, dissipator.solidport)
        connect(base.hp2, radiator.hp2, coolant_tem.hp, dissipator.fluidport)
    ]
    @named rad = ODESystem(eqs, t, systems=[base, gas_tem, radiator, dissipator, coolant_tem])
    sys = structural_simplify(rad)

    u0 = [
        base.Q_flow => 10
        dissipator.Q_flow => 10
    ]
    prob = ODEProblem(sys, u0, (0, 3.0))
    sol = solve(prob, Rosenbrock23())

    @test sol[dissipator.dT] == sol[radiator.hp1.T] - sol[radiator.hp2.T]
    rad_Q_flow = G*σ*(Tᵧ^4 - Tᵪ^4)
    @test sol[radiator.Q_flow] == fill(rad_Q_flow, length(sol[radiator.Q_flow]))
end

@testset "Thermal Collector" begin
    @named flow_src    = FixedHeatFlow(Q_flow=50, α=100)
    @named hf_sensor   = HeatFlowSensor()
    @named th_ground   = ThermalGround()
    @named collector   = ThermalCollector(N=2)
    @named th_resistor = ThermalResistor(R=10)
    @named tem_src     = FixedTemperature(T=10)

    @info "Building a heat collector..."
    eqs = [
        connect(flow_src.hp, collector.hp1, th_resistor.hp1)
        connect(tem_src.hp, collector.hp2)
        connect(hf_sensor.hp1, collector.collector_port)
        connect(hf_sensor.hp2, th_ground.hp, th_resistor.hp2)
        ]
    @named coll = ODESystem(eqs, t,
                            systems=[hf_sensor,flow_src, tem_src,
                            collector, th_resistor])
    sys = structural_simplify(coll)

    u0 = [
        th_resistor.Q_flow => 1.0
    ]
    prob = ODEProblem(sys, u0, (0, 3.0))
    sol  = solve(prob, Rosenbrock23())

    @test sol[collector.collector_port.Q_flow] + sol[collector.hp1.Q_flow] + sol[collector.hp2.Q_flow] ==
        zeros(length(sol[collector.collector_port.Q_flow]))
    @test sol[collector.collector_port.T] == sol[collector.hp1.T] == sol[collector.hp2.T]
end


text:
using Test, OceanBasins

ENV["DATADEPS_ALWAYS_ACCEPT"] = true
const OCEANS = oceanpolygons()

@testset "Testing Oceanbasins" begin
    @test ispacific(0, -160, OCEANS)
    @test isarctic(89, 0, OCEANS)
    @test isatlantic(0, -30, OCEANS)
    @test isindian(0, 90, OCEANS)
    @test ispacific2(0, -160, OCEANS)
    @test isatlantic2(0, -30, OCEANS)
    @test isindian2(0, 90, OCEANS)
    @test isSOpacific(-80, -160, OCEANS)
    @test isSOatlantic(-80, -30, OCEANS)
    @test isSOindian(-80, 90, OCEANS)
    @test isantarctic(-40, 0, OCEANS)
    @test ismediterranean(31, 15, OCEANS)
    @test isindonesian(0, 120, OCEANS)
    @test ismediterranean2(31, 15, OCEANS)
    @test isNPAC(30, -160, OCEANS)
    @test isNATL(30, -30, OCEANS)
    @test isSPAC(-30, -160, OCEANS)
    @test isSATL(-30, -30, OCEANS)
end


text:
input = readlines("./src/Day3/input.txt")

function count_trees(input, slope, down)
    pattern = copy(input)
    numtrees = 0
    pos = 1
    width = length(pattern[1])
    for row in 1:down:length(pattern)
        pos = mod(pos - 1, width) + 1
        val = pattern[row][pos]
        letters = string.(split(pattern[row],""))
        if val == '#'
            numtrees += 1
            letters[pos] = "X"
        else
            # print(val)
            letters[pos] = "O"
        end
        pattern[row] = join(letters)
        pos += slope
    end
    return numtrees, pattern
end

hill = copy(input)
num11, pat = count_trees(hill, 1, 1)
num31, pat = count_trees(hill, 3, 1)
num51, pat = count_trees(hill, 5, 1)
num71, pat = count_trees(hill, 7, 1)
num12, pat = count_trees(hill, 1, 2)

num11*num31*num51*num71*num12


text:
@testset "Utility and Printing Functions" begin
    @testset "Pretty printing of results" begin
        @testset "Pretty printing of TopsisResult" begin
            t = TopsisResult(
                DataFrame(rand(3, 4), :auto),
                rand(3),
                DataFrame(rand(3, 4), :auto),
                DataFrame(rand(3, 4), :auto),
                1,
                [1, 2, 3],
            )
            io = IOBuffer()
            show(io, t)
            str_expected = String(take!(io))
            @test str_expected == "Scores:\n[1.0, 2.0, 3.0]\nBest indices:\n1\n"
        end

        @testset "Pretty printing of ARAS Result" begin
            t = ARASResult(
                rand(3),
                rand(3, 3),
                rand(3, 3),
                rand(3),
                [1.0, 2.0, 3.0],
                [3, 2, 1],
                2,
            )
            io = IOBuffer()
            show(io, t)
            str_expected = String(take!(io))
            @test str_expected ==
                  "Scores:\n[1.0, 2.0, 3.0]\nOrderings: \n[3, 2, 1]\nBest indices:\n2\n"
        end
    end

    @testset "Utility functions" begin
        @testset "Identity matrix" begin
            @test JMcDM.I(2) == [1.0 0; 0 1]
            @test JMcDM.I(3) == [1.0 0 0; 0 1 0; 0 0 1]
            @test JMcDM.I(4) == [1.0 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
        end

        @testset "mean" begin
            @test JMcDM.mean([1.0, 2.0, 3.0]) == 2.0
            @test JMcDM.mean([1.0, 2.0, 3.0, 4.0]) == 2.5
        end

        @testset "var" begin
            @test JMcDM.var([1.0, 2.0, 3.0]) == 1.0
            @test isapprox(JMcDM.var([1.0, 2.0, 3.0, 4.0]), 1.666667, atol = 0.001)
        end

        @testset "std" begin
            @test JMcDM.std([1.0, 2.0, 3.0]) == 1.0
            @test isapprox(JMcDM.std([1.0, 2.0, 3.0, 4.0]), 1.290994, atol = 0.001)
        end

        @testset "geomean" begin
            @test isapprox(JMcDM.geomean([1.0, 2.0, 3.0]), 1.817121, atol = 0.001)
            @test isapprox(JMcDM.geomean([1.0, 2.0, 3.0, 4.0]), 2.213364, atol = 0.001)
        end

        @testset "cor" begin
            x = [1.0, 2.0, 3.0, 4.0, 5.0]
            y = [5.0, 4.0, 3.0, 2.0, 1.0]
            @test JMcDM.cor(x, x) == 1.0
            @test JMcDM.cor(x, y) == -1.0

            mat = hcat(x, y)
            cormat = JMcDM.cor(mat)
            @test cormat[1, 1] == 1.0
            @test cormat[1, 2] == -1.0
            @test cormat[2, 1] == -1.0
            @test cormat[2, 2] == 1.0
        end

        @testset "Euclidean distance" begin
            @test euclidean([0.0, 1.0, 2.0], [0.0, 1.0, 2.0]) == 0.0
            @test euclidean([0.0, 0.0, 0.0]) == 0.0
            @test euclidean([0.0, 0.0, 1.0]) == 1.0
            @test euclidean([0.0, 0.0, 1.0], [0.0, 0.0, 2.0]) == 1.0
            @test euclidean([0, 0, 0]) == 0.0
            @test euclidean([0, 0, 0, 0], [1, 0, 0, 0]) == 1.0
            @test euclidean([0, 0, 0, 0], [1, 0, 0, 0]) isa Float64
            @test euclidean([0, 0, 0, 0], [1.0, 0, 0, 0]) isa Float64
        end

        @testset "Normalization" begin
            tol = 0.00001
            nz = normalize([1.0, 2.0, 3.0, -1.0, 0.0])
            @test isapprox(nz[1], 0.2581989, atol = tol)
            @test isapprox(nz[2], 0.5163978, atol = tol)
            @test isapprox(nz[3], 0.7745967, atol = tol)
            @test isapprox(nz[4], -0.2581989, atol = tol)
            @test isapprox(nz[5], 0.0000000, atol = tol)

            nzint = normalize([1, 2, 3, -1, 0])
            @test isapprox(nzint[1], 0.2581989, atol = tol)
            @test isapprox(nzint[2], 0.5163978, atol = tol)
            @test isapprox(nzint[3], 0.7745967, atol = tol)
            @test isapprox(nzint[4], -0.2581989, atol = tol)
            @test isapprox(nzint[5], 0.0000000, atol = tol)
        end

        @testset "Column min and max vectors" begin
            df = DataFrame()
            df[:, :x] = [0.0, 1.0, 10.0]
            df[:, :y] = [0.0, -1.0, -10.0]
            @test colmins(df) == [0.0, -10.0]
            @test colmaxs(df) == [10.0, 0.0]
        end

        @testset "Unitize vector" begin
            x = [1.0, 2.0, 3.0, 4.0, 5.0]
            result = x |> unitize |> sum
            @test result == 1.0

            x = [1, 2, 3, 4, 5]
            result = x |> unitize |> sum
            @test result == 1.0

            x = [1, 1, 1, 1, 1]
            result = x |> unitize
            @test result == [0.20, 0.20, 0.20, 0.20, 0.20]

            x = [0.20, 0.20, 0.20, 0.20, 0.20]
            result = x |> unitize
            @test result == [0.20, 0.20, 0.20, 0.20, 0.20]
        end

        @testset "Product weights with DataFrame" begin
            df = DataFrame()
            df[:, :x] = [1.0, 2.0, 4.0, 8.0]
            df[:, :y] = [10.0, 20.0, 30.0, 40.0]
            w = [0.60, 0.40]
            result = w * df
            @test result[:, :x] == [0.6, 1.2, 2.4, 4.8]
            @test result[:, :y] == [4.0, 8.0, 12.0, 16.0]

            dfint = DataFrame()
            dfint[:, :x] = [1, 2, 4, 8]
            dfint[:, :y] = [10, 20, 30, 40]
            w = [0.60, 0.40]
            result = w * dfint
            @test result[:, :x] == [0.6, 1.2, 2.4, 4.8]
            @test result[:, :y] == [4.0, 8.0, 12.0, 16.0]
        end

        @testset "Make Decision Matrix" begin
            m = rand(5, 10)
            df = makeDecisionMatrix(m)

            @test isa(df, DataFrame)
            @test size(df) == (5, 10)
            @test df[:, 1] isa Array{Float64,1}
            @test df[:, 2] isa Array{Float64,1}
            @test df[:, 3] isa Array{Float64,1}
            @test df[:, 4] isa Array{Float64,1}
            @test df[:, 5] isa Array{Float64,1}

            @test names(df)[1] == "Crt1"
            @test names(df)[2] == "Crt2"

            m = rand(3, 5)
            dfwithnames = makeDecisionMatrix(m, names = ["A1", "B", "CD", "EE", "FG"])
            @test names(dfwithnames)[1] == "A1"
            @test names(dfwithnames)[2] == "B"
            @test names(dfwithnames)[3] == "CD"
            @test names(dfwithnames)[4] == "EE"
            @test names(dfwithnames)[5] == "FG"
        end

        @testset "Reverse minimum & maximum array" begin

            fns = [minimum, maximum, maximum, minimum, maximum]
            revfns = [maximum, minimum, minimum, maximum, minimum]

            @test reverseminmax(fns) == revfns

            @test reverseminmax(revfns) == fns
        end

        @testset "Make Array of minimum and maximum" begin
            result1 = makeminmax([maximum, maximum, maximum, maximum])

            @test typeof(result1) == Array{Function,1}

            @test typeof(result1[1]([1.0, 2.0, 3.0])) == Float64

            @test result1[1]([1.0, 2.0, 3.0]) == 3.0

        end
    end
end

max_stars_repo_name: drvinceknight/JMcDM

max_issues_repo_name: drvinceknight/JMcDM

max_forks_repo_name: drvinceknight/JMcDM


text:
abstract type AbstractBitstringProblem <: AbstractProblem end

import Base: length

function fitness(x::Array{Bool,1}, problem)
    fitness(BitArray(x), problem)
end

function fitness(
    x::BitArray,
    problem::T,
    fitnessvalue::Number,
    flipsposition,
) where {T<:AbstractProblem}
    y = copy(x)
    for item ∈ flipsposition
        flip!(y, item)
    end
    fitness(y, problem)
end
"""
    struct OneMax <: AbstractBitstringProblem
        n::Integer
    end

Maximization of the number of `true` in a `BitArray`. The mathematical definition of this problem is
maximizing the number of one-bits in a bitstring, which can be seen in the following formula:

```math
OneMax(x_1,\\dots,x_n)\\coloneqq\\sum_{i=1}^{n}x_i.
```

## Examples
```jldoctest; setup = :(using Pidoh)
julia> problem = OneMax(4);

julia> fitness(BitArray([true,true,false,true]), problem)
3

julia> optimum(problem)
4

```
"""
struct OneMax <: AbstractBitstringProblem
    n::Integer
end
fitness(x::BitArray, problem::OneMax) = sum(x)
optimum(problem::OneMax) = problem.n
length(problem::OneMax) = problem.n

"""
    struct LeadingOnes <: AbstractBitstringProblem
        n::Integer
    end

Maximization of the length of a prefix of `BitArray`  that consits only of `true`. More formally

```math
LeadingOnes(x_1,\\dots,x_n)\\coloneqq\\sum_{i=1}^{n}\\prod_{j=1}^ix_j.
```

## Examples
```jldoctest; setup = :(using Pidoh)
julia> problem = LeadingOnes(4);

julia> fitness(BitArray([true,true,false,true]), problem)
2

julia> optimum(problem)
4

```
"""
struct LeadingOnes <: AbstractBitstringProblem
    n::Integer
end

function fitness(x::BitArray, problem::LeadingOnes)
    lo = 0
    while lo < problem.n && x[lo+1]
        lo += 1
    end
    lo
end

optimum(problem::LeadingOnes) = problem.n
length(problem::LeadingOnes) = problem.n

"""
    struct CliffTwoParameters <: AbstractBitstringProblem
        n::Integer
    end

## Examples
```jldoctest; setup = :(using Pidoh)
julia> problem = CliffTwoParameters(10, 4, 2);

```
"""
struct CliffTwoParameters <: AbstractBitstringProblem
    n::Integer
    m::Integer
    d::Real
end

function fitness(x::BitArray, problem::CliffTwoParameters)
    onemax = sum(x)
    n = problem.n
    m = problem.m
    d = problem.d
    if onemax ≤ n - m
        return onemax * 1.0
    else
        return onemax*1.0 - d - 1
    end
end

optimum(problem::CliffTwoParameters) = problem.n - problem.d - 1
length(problem::CliffTwoParameters) = problem.n

struct ZeroMax <: AbstractBitstringProblem
    n::Integer
end

fitness(x::BitArray, problem::ZeroMax) = problem.n - sum(x)
optimum(problem::ZeroMax) = 0

struct Jump <: AbstractBitstringProblem
    n::Integer
    jumpsize::Integer
end

function fitness(x::BitArray, problem::Jump)
    onemax = sum(x)
    n = problem.n
    k = problem.jumpsize
    jump = -1
    if ((onemax ≤ n - k) || (onemax == n))
        jump = k + onemax
    else
        jump = n - onemax
    end
    jump
end

optimum(problem::Jump) = problem.n + problem.jumpsize
length(problem::Jump) = problem.n

struct GeneralJump <: AbstractBitstringProblem
    n::Integer
    k::Integer
    δ::Integer
end

function fitness(x::BitArray, problem::GeneralJump)
    onemax = sum(x)

    if onemax ≤ problem.n - problem.k || onemax ≥ problem.n-problem.k+problem.δ
        return onemax
    else
        return -onemax
    end
end

optimum(problem::GeneralJump) = problem.n
length(problem::GeneralJump) = problem.n

# function fitness(x::BitArray, problem::OneMax, fitnessvalue:: Number, flipsposition)
#     if length(flipsposition) == 0
#         return fitnessvalue
#     end
#     fit = fitnessvalue
#     for item ∈ flipsposition
#         @inbounds if x[item]
#             fit -= 1
#         else
#             fit += 1
#         end
#     end
#     fit
# end

# function fitness(x::BitArray, problem::Jump, fitnessvalue:: Number, flipsposition)
#     if length(flipsposition) == 0
#         return fitnessvalue
#     end
#     n=problem.n
#     k = problem.jumpsize
#
#     onemax = fitnessvalue
#
#     if fitnessvalue < k
#         onemax = n - fitnessvalue
#     else
#         onemax = fitnessvalue - k
#     end
#
#     if ((onemax ≤ n - k) || (onemax == n))
#         jump =k + onemax
#     else
#         jump = n - onemax
#     end
#
#     for item ∈ flipsposition
#         @inbounds if x[item]
#             onemax -= 1
#         else
#             onemax += 1
#         end
#     end
#
#     jump = -1
#     if ((onemax ≤ n - k) || (onemax == n))
#         jump =k + onemax
#     else
#         jump = n - onemax
#     end
#     jump
# end

max_stars_repo_name: DTUComputeTONIA/SDRLS

max_issues_repo_name: DTUComputeTONIA/SDRLS

max_forks_repo_name: DTUComputeTONIA/SDRLS


text:
using Plots

function binnedboxplot!(
    df::DataFrame,
    metric::Union{Symbol, String},
    var::Union{Symbol, String},
    bins::AbstractVector;
    target = missing;
    size = (800, 800),
    id = :ID,
    α = 3,
    color = :lightblue,
)
    individuals = length(unique(df[:, id]))
    data = PopPKBase.bin(df[:, metric], bins, bins)
    for i in bins
        relevant = filter(x -> (j < x[metrics[i]] < j + 1) && (x[var] !== missing), df)
        data = relevant[:, var]
        boxplot!(
            [j + 0.5],
            isempty(data) ? [0] : data,
            outliers = false,
            color = color,
            subplot = i,
            title = metrics[i],
            α = α * nrow(relevant) / individuals
        )
    end
    hline!(transpose([target for i = 1:length(metrics)*2]), color = :red)
    ylabel!("Count")
    xlabel!("Z-Score")
end


text:
export freqContBasic;
export freqContTraceEstimation;
using jInv.InverseSolve
# using JLD
using KrylovMethods

"""
	function freqContBasic
	Frequency continuation procedure for running FWI.
	This function runs GaussNewton on misfit functions defined by pMis with nfreq frequencies.

	Input:
		mc    		- current model
		pInv		- Inverse param
		pMis 		- misfit params (remote)
		nfreq		- number of frequencies in the problem
		windowSize  - How many frequencies to treat at once at the most.
		dumpFun     - a function for plotting, saving and doing all the things with the intermidiate results.
		resultsFilename - file name to save results
		startFrom   - a start index for the continuation. Usefull when our run broke in the middle of the iterations.

"""
function freqContBasic(mc, pInv::InverseParam, pMis::Array{RemoteChannel},nfreq::Int64, windowSize::Int64,
			dumpFun::Function, resultsFilename::String, startFrom::Int64 = 1)
Dc = 0;
flag = -1;
HIS = [];
for freqIdx = startFrom:nfreq
	println("start freqCont iteration from: ", freqIdx)
	tstart = time_ns();

	reqIdx1 = freqIdx;
	if freqIdx > 1
		reqIdx1 = max(1,freqIdx-windowSize+1);
	end
	reqIdx2 = freqIdx;
	currentProblems = reqIdx1:reqIdx2;
	println("\n======= New Continuation Stage: selecting continuation batches: ",reqIdx1," to ",reqIdx2,"=======\n");
	pMisTemp = pMis[currentProblems];
	pInv.mref = mc[:];

	if resultsFilename == ""
			filename = "";
		else
			Temp = splitext(resultsFilename);
			filename = string(Temp[1],"_FC",freqIdx,"_GN",Temp[2]);
	end
	# Here we set a dump function for GN for this iteracion of FC
	function dumpGN(mc,Dc,iter,pInv,PF)
		dumpFun(mc,Dc,iter,pInv,PF,filename);
	end

	mc,Dc,flag,His = projGNCG(mc,pInv,pMisTemp,dumpResults = dumpGN);

	Dc,F,dF,d2F,pMisNone,times,indDebit = computeMisfit(mc, map(pm -> fetch(pm), pMisTemp));
	println("Misfit for freq: ", freqIdx, " after 5 GN iterations is: ", F);
	clear!(pMisTemp);
	tend = time_ns();
    println("Runtime of freqCont iteration: ");
    println((tend - tstart)/1.0e9);
end
return mc,Dc,flag,HIS;
end

function calculateZ2(misfitCalc::Function, p::Integer, nsrc::Integer,
	numOfCurrentProblems::Integer, Wd::Array, HinvPs::Array,
	pMisCurrent::Array{MisfitParam}, Z1::Matrix, alpha::Float64)

	println("misfit at start:: ", misfitCalc())
	rhs = zeros(ComplexF64, (p, nsrc));
	for i = 1:numOfCurrentProblems
		pm = pMisCurrent[i]
		println("size wd: ", mean(Wd[i]))

		rhs += (mean(Wd[i])^2) .* Z1' * HinvPs[i] * (-HinvPs[i]' * pm.pFor.Sources + pm.dobs[:,:,1]);
	end
	lhs = zeros(ComplexF64, (p,p));
	println("size Wd: ", size(Wd[1]));
	println("size HINV: ", size(HinvPs[1]));
	for i = 1:numOfCurrentProblems
		lhs += (mean(Wd[i])^2) .* Z1' * HinvPs[i] * HinvPs[i]' * Z1;
	end

	lhs += alpha * I;

	Z2 = lhs\rhs;
end

function freqContExtendedSources(mc, pInv::InverseParam, pMis::Array{RemoteChannel},
	nfreq::Int64, windowSize::Int64, Iact,
	mback::Union{Vector,AbstractFloat,AbstractModel, Array{Float64,1}},
	dumpFun::Function, resultsFilename::String, startFrom::Int64 = 1)
Dc = 0;
flag = -1;
HIS = [];

pFor = fetch(pMis[1]).pFor
mSizeMat = pFor.Mesh.n .+ 1;
m = mSizeMat[1];
n = mSizeMat[2];
mSizeVec = mSizeMat[1] * mSizeMat[2];
# Z = copy(pFor.originalSources);
# sizeQ = size(Z);
nrec = size(pFor.Receivers, 2);
sizeH = size(pFor.Ainv[1]);

nsrc = size(pFor.Sources, 2);
pFor = nothing
alpha = 2e-3;
p = 10;

for freqIdx = startFrom:nfreq
	Z1 = rand(ComplexF64,(m*n, p)) .+ 0.01;
	Z2 = rand(ComplexF64, (p, nsrc)) .+ 0.01;
	println("start freqCont Zs iteration from: ", freqIdx)
	tstart = time_ns();
	reqIdx1 = freqIdx;
	if freqIdx > 1
		reqIdx1 = max(1,freqIdx-windowSize+1);
	end
	reqIdx2 = freqIdx;
	currentProblems = reqIdx1:reqIdx2;

	println("\n======= New Continuation Stage: selecting continuation batches: ",reqIdx1," to ",reqIdx2,"=======\n");
	runningProcs = map(x->x.where, pMis[currentProblems]);

	for j=1:5

		pMisCurrent = map(fetch, pMis[currentProblems]);
		pForpCurrent =  map(x->x.pFor, pMisCurrent);
		Dp,pForp = getData(vec(mc), pForpCurrent);

		pForCurrent = map(x->fetch(x), pForp);
		numOfCurrentProblems = size(currentProblems, 1);
		map((pm,pf) -> pm.pFor = pf , pMisCurrent, pForCurrent);
		HinvPs = Vector{Array}(undef, numOfCurrentProblems);
		pMisTemp = Array{RemoteChannel}(undef, length(currentProblems));
		t1 = time_ns();
		for freqs = 1:numOfCurrentProblems
			HinvPs[freqs] = (pForCurrent[freqs].Ainv[1])' \ Matrix(pForCurrent[freqs].Receivers);
			println("HINVP done");
		end
		e1 = time_ns();
		println("runtime of HINVPs");
		println((e1 - t1)/1.0e9);

		if resultsFilename == ""
				filename = "";
			else
				Temp = splitext(resultsFilename);
				filename = string(Temp[1],"_FC",freqIdx,"_",j,"_GN",Temp[2]);
		end

		function dumpGN(mc,Dc,iter,pInv,PF)
			dumpFun(mc,Dc,iter,pInv,PF,filename);
		end

		t2 = time_ns();
		Dc,F,dF,d2F,pMisNone,times,indDebit = computeMisfit(mc, pMisCurrent);
		e2 = time_ns();
		println("Misfit B4 mzs at GN ", j, "frequncy idx: ", freqIdx, " Is: ", F);
		println("runtime of compute misfit");
		println((e2 - t2)/1.0e9);

		t3 = time_ns();

		pMisArr = map(pm -> fetch(pm), pMisCurrent);
		Ap = Vector{Matrix}(undef, numOfCurrentProblems);
		diags = Vector{SparseMatrixCSC}(undef, numOfCurrentProblems);
		B = Vector{Vector}(undef, numOfCurrentProblems);

		toVec(mat) = reshape(mat, mSizeVec);
		toMat(vec) = reshape(vec, tuple(mSizeMat...));

		Wd = map(pm -> pm.Wd[:,:,1], pMisCurrent)
		function misfitCalc()
			sum = 0;
			for i = 1:numOfCurrentProblems
				sum += (mean(Wd[i])^2) .* norm(HinvPs[i]' * (pMisCurrent[i].pFor.Sources + Z1 * Z2) - pMisCurrent[i].dobs[:,:,1])^2;
			end

			sum	+= alpha * norm(Z1)^2 + alpha * norm(Z2)^2;
			return sum;
		end

		Z2 = calculateZ2(misfitCalc, p, nsrc, numOfCurrentProblems, Wd, HinvPs,
			pMisCurrent, Z1, alpha);

		println("misfit at Z2:: ", misfitCalc())

		function multOP(R,HinvP)
			return HinvP' * R * Z2;
		end

		function multOPT(R,HinvP)
			return HinvP * R * Z2';
		end

		function multAll(x)
			sum = zeros(ComplexF64, (mSizeVec, p));

			for i = 1:numOfCurrentProblems
				sum += (mean(Wd[i])^2) .* multOPT(multOP(x, HinvPs[i]), HinvPs[i]);
			end

			sum += alpha * x;
			return sum;
		end

		rhs = zeros(ComplexF64, (mSizeVec, p));
		for i = 1:numOfCurrentProblems
			pm = pMisCurrent[i]
			rhs += (mean(Wd[i])^2).*multOPT(-HinvPs[i]' * pm.pFor.Sources + pm.dobs[:,:,1], HinvPs[i]);
		end

		Z1 = KrylovMethods.blockBiCGSTB(x-> real(multAll(x)), real(rhs))[1];

		for i = 1:numOfCurrentProblems
			pMisCurrent[i].pFor.Sources += Z1 * Z2;
		end

		println("misfit at Z1:: ", misfitCalc())

		Dc,F,dF,d2F,pMisNone,times,indDebit = computeMisfit(mc, pMisCurrent);

		println("Misfit after GN ", j, "frequncy idx: ", freqIdx, " Is: ", F);
		pMispCurrent = Array{RemoteChannel}(undef, numOfCurrentProblems);
		for i=1:numOfCurrentProblems
			pMispCurrent[i] = initRemoteChannel(x->x, runningProcs[i], pMisCurrent[i]);
		end
		pInv.mref = mc[:];
		t4 = time_ns();
		mc,Dc,flag,His = projGNCG(mc,pInv,pMispCurrent,dumpResults = dumpGN);
		e4 = time_ns();
		println((t4 - e4)/1.0e9);
		Dc,F,dF,d2F,pMisNone,times,indDebit = computeMisfit(mc,pMisCurrent);

		println("Misfit after GN ", j, "frequncy idx: ", freqIdx, " Is: ", F);

		pMis[currentProblems] = pMispCurrent;
		clear!(pMispCurrent);
	end

	tend = time_ns();
    println("Runtime of freqCont iteration: ");
    println((tend - tstart)/1.0e9);
	global inx = inx + 1;

end
return mc,Dc,flag,HIS;
end

"""
	function freqContTraceEstimation
	Frequency continuation procedure for running FWI.
	this function first reduces the number of sources by trace estimation,
	then minimizing the misfit using alternating minimization between dsj and m

	Input:
		mc    		- current model
		pInv		- Inverse param
		pMis 		- misfit params (remote)
		nfreq		- number of frequencies in the problem
		windowSize  - How many frequencies to treat at once at the most.
		Iact        - Projector to active cells.
		sigmaBack   - Background model ("frozen" cells).
		dumpFun     - a function for plotting, saving and doing all the things with the intermidiate results.
		resultsFilename - file name to save results
		startFrom   - a start index for the continuation. Usefull when our run broke in the middle of the iterations.

"""
function freqContTraceEstimation(mc, pInv::InverseParam, pMis::Array{RemoteChannel},nfreq::Int64, windowSize::Int64,
			Iact,mback::Union{Vector,AbstractFloat,AbstractModel},
			dumpFun::Function, resultsFilename::String, startFrom::Int64 = 1)
Dc = 0;
flag = -1;
HIS = [];

for freqIdx = startFrom:nfreq
	println("start freqCont iteration from: ", freqIdx)
	tstart = time_ns();

	reqIdx1 = freqIdx;
	if freqIdx > 1
		reqIdx1 = max(1,freqIdx-windowSize+1);
	end
	reqIdx2 = freqIdx;
	currentProblems = reqIdx1:reqIdx2;
	pMisTE = calculateReducedMisfitParams(mc, currentProblems, pMis, Iact, mback);
	pInv.mref = mc[:];

	println("\n======= New Continuation Stage: selecting continuation batches: ",reqIdx1," to ",reqIdx2,"=======\n");

	if resultsFilename == ""
			filename = "";
		else
			Temp = splitext(resultsFilename);
			filename = string(Temp[1],"_FC",freqIdx,"_GN",Temp[2]);
	end
	# Here we set a dump function for GN for this iteracion of FC
	function dumpGN(mc,Dc,iter,pInv,PF)
		dumpFun(mc,Dc,iter,pInv,PF,filename);
	end

	mc,Dc,flag,His = projGNCG(mc,pInv,pMisTE,dumpResults = dumpGN);

	clear!(pMisTE);
	tend = time_ns();
    println("Runtime of freqCont iteration: ");
    println((tend - tstart)/1.0e9);
end

return mc,Dc,flag,HIS;
end

"""
	Function to calculate MistfitParams of new dimensions after trace estimation
"""
function calculateReducedMisfitParams(mc, currentProblems::UnitRange, pMis::Array{RemoteChannel},
			Iact,mback::Union{Vector,AbstractFloat,AbstractModel})
	numOfCurrentProblems = size(currentProblems, 1);

	eta = 20.0;
	newDim = 10;

	runningProcs = map(x->x.where, pMis[currentProblems]);
	pMisCurrent = map(fetch, pMis[currentProblems]);
	pForpCurrent =  map(x->x.pFor, pMisCurrent);

	DobsCurrent = map(x->x.dobs[:,:], pMisCurrent);
	WdCurrent = map(x->x.Wd[:,:], pMisCurrent);
	DobsNew = copy(DobsCurrent[:]);
	nsrc = size(DobsNew[1],2);
	TEmat = rand([-1,1],(nsrc,newDim));
	WdEta = eta.*WdCurrent;
	DpCurrent, = getData(vec(mc),pForpCurrent);
	WdNew = Array{Array}(undef, numOfCurrentProblems);
	for i=1:numOfCurrentProblems
		DobsTemp = fetch(DpCurrent[i]);
		for s=1:nsrc
			etaM = 2 * eta .*diagm(0 => vec(WdCurrent[i][:,s])).*diagm(0 => vec(WdCurrent[i][:,s]));
			A = (2 .*I + etaM)
			b = (etaM)*(A\DobsCurrent[i][:,s])
			DobsNew[i][:,s] = 2 .* (A\DobsTemp[:,s]) + b;
		end
		DobsNew[i] = DobsNew[i][:,:] * TEmat;
		WdNew[i] = 1.0./(sqrt(newDim) .* abs.(real.(DobsNew[i])) .+ 1e-1*mean(abs.(DobsNew[i])));
	end

	pForReduced = Array{RemoteChannel}(undef, numOfCurrentProblems);
	for i=1:numOfCurrentProblems
		pForTemp = pForpCurrent[i];
		pForTemp.Sources = pForTemp.originalSources * TEmat;
		pForReduced[i] = initRemoteChannel(x->x, runningProcs[i], pForTemp);
	end
	return getMisfitParam(pForReduced, WdNew, DobsNew, SSDFun, Iact, mback);
end

function dummy(mc,Dc,iter,pInv,pMis)
# this function does nothing and is used as a default for dumpResults().
end


text:
mutable struct DataPoint
    features::Vector{Float64}
    bid_price::Float64
    pay_price::Float64
end

mutable struct DataSet
    samples::Vector{DataPoint}
end
function DataSet(X::Vector{Vector{T}}, bid_prices::Vector{R}, pay_prices::Vector{S}) where {T, R, S <: Real}
    if !(length(X) == length(bid_prices) == length(pay_prices))
        error("Data set with sample size mismatch.")
    end
    if isempty(X)
        error("Cannot instantiate a problem with no data points.")
    end
    num_features = length(first(X))
    for point in X
        if length(point) != num_features
            error("Data points with different numbers of features.")
        end
    end
    for i in 1:length(bid_prices)
        if bid_prices[i] < pay_prices[i]
            error("Cannot set bid price below pay price.")
        end
        if pay_prices[i] < 0
            error("Cannot have negative pay price.")
        end
    end
    return DataSet([DataPoint(X[i], bid_prices[i], pay_prices[i]) for i in 1:length(X)])
end

function sample_count(dataset::DataSet)
    return length(dataset.samples)
end

function feature_count(dataset::DataSet)
    return length(dataset.samples[1].features)
end

function _max_abs_val(dataset::DataSet)
    return [maximum(abs.(sample.features[i] for sample in dataset.samples)) for i in 1:feature_count(dataset)]
end

function split_data(dataset::DataSet, num_train_points, num_validation_points, num_test_points)
    samples = dataset.samples
    @assert length(samples) == num_train_points + num_validation_points + num_test_points
    return samples[1:num_train_points],
           samples[(num_train_points+1):(num_train_points+num_validation_points)],
           samples[(num_train_points+num_validation_points+1):(num_train_points+num_validation_points+num_test_points)]
end

max_stars_repo_name: joehuchette/reserve-price-optimization

max_issues_repo_name: joehuchette/reserve-price-optimization

max_forks_repo_name: joehuchette/reserve-price-optimization


text:
# Bi-objective linear problem (bilp)
#
# Example 6.3 (from Steuer, 1985), page 154 of
# Multicriteria Optimization (2nd edt), M. Ehrgott, Springer 2005.

# ---- Packages to use
using vOptGeneric, JuMP, GLPK

# ---- setting the model + values
bilp = vModel( GLPK.Optimizer )

@variable( bilp, x1 >= 0 )
@variable( bilp, x2 >= 0 )

@addobjective( bilp, Min, 3*x1 + x2 )
@addobjective( bilp, Min, -x1 - 2*x2 )

@constraint( bilp, cst1, x2 <= 3 )
@constraint( bilp, cst2, 3*x1 - x2 <= 6 )

# ---- Invoking the solver (lexicographic method)
vSolve( bilp, method=:lex )

# ---- Querying the results
Y_N = getY_N( bilp )

# ---- Displaying the results for lex(1,2) and lex(2,1)
for s = 1:2
    print("X = [")
    print(round(getvOptData(bilp).X_E[s][1],digits=5))
    print(" ")
    print(round(getvOptData(bilp).X_E[s][2],digits=5))
    print("] | Z = [")
    print(round(Y_N[s][1], digits=5))
    print(" ")
    print(round(Y_N[s][2], digits=5))
    println("]")
end

max_stars_repo_name: vOptSolver/vOptGeneric

max_issues_repo_name: vOptSolver/vOptGeneric

max_forks_repo_name: vOptSolver/vOptGeneric


text:
"""
lookup : keeps track of code that's main function is to look up external data.
"""

module lookup

using DelimitedFiles, Statistics

function get_energy_params(gene_name, param_dir="energies/") #NOTE: extra dependency
	deg_energy_file = open("$param_dir" * "energy_$(gene_name).txt")
	deg_energy_read = readlines(deg_energy_file)
	close(deg_energy_file)

	#energies = deg_energy[strip(split(der, '\t')[2]) for der in deg_energy_read]
	energies = map(xx -> parse(Float64, xx), deg_energy_read)
	mean_energy = mean(energies)
	std_energy = std(energies)
	return mean_energy, std_energy
end

function mu_sig(gene_name, param_dir="psls/") #NOTE: extra dependency
	for csv_file in readdir(param_dir)
		if csv_file == "psls_$(gene_name).txt"
			#Then it's a match
			raw_data = readdlm(param_dir * "/" * csv_file)
			raw_values = map(Float64, raw_data[1:end, 1])
			mu = mean(raw_values)
			sig = std(raw_values)
			return mu, sig
		end
	end
	return -1, -1 #dummy vars for not successful.
end

end


text:
export Algorithm

abstract type Algorithm end

abstract type Detect <: Algorithm end
abstract type Align <:Algorithm end
abstract type Cluster <:Algorithm end
abstract type Feature <:Algorithm end
abstract type Reduction <:Algorithm end
abstract type Threshold <:Algorithm end

abstract type Sorting end

#Data Structure to store range of a spike and cluster ID
struct Spike
    inds::UnitRange{Int64}
    id::UInt8
end

Spike()=Spike(0:0,0)

struct s_sizes
    win::Int16
    s_end::Int16
    win2::Int16
end

s_sizes(win)=s_sizes(win,win+div(win,2),div(win,2))

function output_buffer(channels::Int64,par=false)

    nums=zeros(UInt16,channels)

    if par==false
        buf=Spike[Spike() for i=1:100,j=1:channels]
    else
        buf=convert(SharedArray{Spike,2},Spike[Spike() for i=1:100,j=1:channels])
        nums=convert(SharedArray{UInt16,1},nums)
    end

    (buf,nums)
end

global sorting_num = 1

function MakeSorting()
end

function gen_sorting(D::Detect,C::Cluster,A::Align,F::Feature,R::Reduction,T::Threshold,in_type=Int16)

    global sorting_num

    f_type=feature_type(F,in_type)

    @eval begin
        mutable struct $(Symbol("Sorting_$sorting_num")) <: Sorting
            d::($(typeof(D)))
            c::($(typeof(C)))
            a::($(typeof(A)))
            f::($(typeof(F)))
            r::($(typeof(R)))
            t::($(typeof(T)))
            id::Int16
            sigend::Array{$(in_type),1}
            index::Int16
            p_temp::Array{$(in_type),1}
            features::Array{$(f_type),1}
            fullfeature::Array{$(f_type),1}
            dims::Array{Int16,1}
            thres::Float64
            #waveform::Array{$(in_type),1}
            cent::Int16
	    s::s_sizes
        end

        function MakeSorting(D::($(typeof(D))),C::($(typeof(C))),A::($(typeof(A))),F::($(typeof(F))),R::($(typeof(R))),T::($(typeof(T))),window::Int64,in_type=Int16)

            #determine size of alignment output
            wavelength=mysize(A,window)

            #determine feature size
            fulllength=mysize(F,wavelength)

            f_type=feature_type(F,in_type)

            if typeof(R)==ReductionNone
                reducedims=fulllength
            else
                R=typeof(R)(fulllength,R.mydims)
                reducedims=R.mydims
            end
            F=typeof(F)(wavelength,reducedims)
            C=typeof(C)(reducedims)
            p_temp=zeros(in_type,window*2)
            $(Symbol("Sorting_$sorting_num"))(D,C,A,F,R,T,
                    1,zeros(in_type,window+div(window,2)),0,
                    p_temp,zeros(f_type,reducedims),zeros(f_type,fulllength),
                    collect(1:reducedims),1.0,1,s_sizes(window))
        end
    end

    sorting_num+=1
    nothing
end

function create_multi(d::Detect,c::Cluster,a::Align,f::Feature,r::Reduction,t::Threshold,num::Int64,window=50,in_type=Int16)

    if method_exists(MakeSorting,(typeof(d),typeof(c),typeof(a),typeof(f),typeof(r),typeof(t),Int64,DataType))
    else
        gen_sorting(d,c,a,f,r,t,in_type)
    end

    if VERSION > v"0.7-"
        st=Array{typeof(Base.invokelatest(MakeSorting,d,c,a,f,r,t,window,in_type))}(undef, num)
    else
        st=Array{typeof(Base.invokelatest(MakeSorting,d,c,a,f,r,t,window,in_type))}(num)
    end

    for i=1:num
        st[i]=Base.invokelatest(MakeSorting,d,c,a,f,r,t,window,in_type)
        st[i].id=i
    end

    st
end

function create_multi(d::Detect,c::Cluster,a::Align,f::Feature,r::Reduction,t::Threshold,num::Int64,cores::UnitRange{Int64},window=50,in_type=Int16)

    st=create_multi(d,c,a,f,r,t,num,window,in_type)
    st=distribute(st,procs=collect(cores))
end


text:
# include the includes -
include("Include.jl")

# mean center -
function mean_center_array(results_array::Array{Float64,2})::Array{Float64,2}

    # get the size -
    (NR,NC) = size(results_array)
    scaled_array = zeros(NR,NC)

    for col_index = 1:NC

        data_col = results_array[:,col_index]
        mu_value = mean(data_col)
        std_value = std(data_col)

        for row_index = 1:NR
            scaled_array[row_index,col_index] = (data_col[row_index] - mu_value)/(std_value)
        end
    end

    return scaled_array
end

# computes the model performance -
function model_performance(parameter_guess_array,index)

    # what is the host_type?
    host_type = :cell_free

    # load the default data_dictionary -
    time_start = 0.0
    time_stop = 16.0
    time_step_size = 0.01

    # what is the host_type?
    host_type = :cell_free

    # path to parameters -
    path_to_biophysical_constants_file = "./CellFree.json"
    #path_to_data_dir = "$(pwd())/data"

    # Load the data dictionary (uses the default biophysical_constants file)
    default_data_dictionary = build_data_dictionary((time_start,time_stop,time_step_size), path_to_biophysical_constants_file, host_type)

    # Set some values for the weights, gene lengths etc -
    model_data_dictionary = customize_data_dictionary(default_data_dictionary,host_type)

   # Phase 1: parameter update =========================================================================== #
   # update the paramaters in the model data dictionary -
   # for now - lets only search over dG's -
   R = model_data_dictionary["R"]
   T_K = model_data_dictionary["T_K"]

   # compute W -
   tmp_W_array = Float64[]
   for index = 1:2
       parameter_guess = parameter_guess_array[index]
       value = exp(-1*parameter_guess/(R*T_K))
       push!(tmp_W_array,value)
   end

   # update the control W's -
   control_parameter_dictionary = model_data_dictionary["control_parameter_dictionary"]
   control_parameter_dictionary["W_deGFP_RNAP"] = tmp_W_array[1]       # 3
   control_parameter_dictionary["W_deGFP_sigma_70"] = tmp_W_array[2]   # 4
   model_data_dictionary["control_parameter_dictionary"] = control_parameter_dictionary

   binding_parameter_dictionary = model_data_dictionary["binding_parameter_dictionary"]
   binding_parameter_dictionary["n_deGFP_sigma_70"] = parameter_guess_array[3]
   binding_parameter_dictionary["K_deGFP_sigma_70"] = parameter_guess_array[4]
   model_data_dictionary["binding_parameter_dictionary"] = binding_parameter_dictionary

   # time constant modifier -
   # time constant modifiers -
   time_constant_modifier_array = [
       0.0	                        ;	# 1	deGFP
       0.0	                        ;	# 2	sigma_70
       parameter_guess_array[5] 	;	# 3	mRNA_deGFP
       1.0	                        ;	# 4	mRNA_sigma_70
       parameter_guess_array[6] 	;	# 5	protein_deGFP
       1.0	                        ;	# 6	protein_sigma_70
   ]
   model_data_dictionary["time_constant_modifier_array"] = time_constant_modifier_array

   # setup degradation_modifier_array -
   degradation_modifier_array = [
       0.0	;	# 1	deGFP
       0.0	;	# 2	sigma_70
       parameter_guess_array[7]	;	# 3	mRNA_deGFP
       1.0	;	# 4	mRNA_sigma_70
       parameter_guess_array[8]	;	# 5	protein_deGFP
       parameter_guess_array[9]	;	# 6	protein_sigma_70
   ]

   # update the translation time -
   model_data_dictionary["half_life_translation_capacity"] = parameter_guess_array[10]

   # lastly, update KL -
   biophysical_constants_dictionary = model_data_dictionary["biophysical_constants_dictionary"]
   biophysical_constants_dictionary["translation_saturation_constant"] = parameter_guess_array[11]
   model_data_dictionary["biophysical_constants_dictionary"] = biophysical_constants_dictionary

   # grab defaults -
   species_symbol_type_array = model_data_dictionary["species_symbol_type_array"]
   protein_coding_length_array = model_data_dictionary["protein_coding_length_array"]
   gene_coding_length_array = model_data_dictionary["gene_coding_length_array"]
   time_constant_modifier_array = model_data_dictionary["time_constant_modifier_array"]
   initial_condition_array = model_data_dictionary["initial_condition_array"]

   # # get gene IC -
   idx_gene = findall(x->x==:gene,species_symbol_type_array)
   gene_abundance_array = initial_condition_array[idx_gene]

   # Precompute the translation parameters -
   translation_parameter_array = precompute_translation_parameter_array(biophysical_constants_dictionary, protein_coding_length_array, time_constant_modifier_array,host_type)
   model_data_dictionary["translation_parameter_array"] = translation_parameter_array

   # Precompute the kinetic limit of transcription -
   transcription_kinetic_limit_array = precompute_transcription_kinetic_limit_array(biophysical_constants_dictionary, gene_coding_length_array, gene_abundance_array, time_constant_modifier_array, host_type)
   model_data_dictionary["transcription_kinetic_limit_array"] = transcription_kinetic_limit_array

   # Dilution degrdation matrix -
   dilution_degradation_matrix = build_dilution_degradation_matrix(biophysical_constants_dictionary, species_symbol_type_array, degradation_modifier_array)
   model_data_dictionary["dilution_degradation_matrix"] = dilution_degradation_matrix
   # ===================================================================================================== #

   # Phase 2:  solve model equations ===================================================================== #
   # solve the balance equations -
   (TSIM,XSIM) = SolveBalances(time_start,time_stop,time_step_size,model_data_dictionary)
   # ===================================================================================================== #

    # Phase 3: compute the model performance metrics ====================================================== #
    p_GFP_AUC = integrate(TSIM,XSIM[:,index])
    # ===================================================================================================== #

    # return the performance_array -
    return p_GFP_AUC
end

function main(path_to_ensemble_file::String,index)

    # setup the sensitivity function -
    SF(P) = model_performance(P,index)

    # setup ranges -
    sample_bounds_array = Array{Tuple,1}()
    ensemble_array = readdlm(path_to_ensemble_file)
    (number_of_parameters,number_of_trials) = size(ensemble_array)
    for parameter_index = 1:(number_of_parameters-1)

        # get row of parameters -
        parameter_row = ensemble_array[parameter_index,:]
        min_value = minimum(parameter_row)
        max_value = maximum(parameter_row)

        # create the tuple -
        tmp_tuple = (min_value,max_value)

        # cache -
        push!(sample_bounds_array,tmp_tuple)
    end

    @show index

    # do the global sensitivity analysis -
    sensitivity_results = gsa(SF,Morris(total_num_trajectory=10000,num_trajectory=1000),sample_bounds_array)

    # return -
    return sensitivity_results
end

# setup paths -
path_to_ensemble_file = "$(pwd())/Ensemble-poets-MVP-T20.dat"

# compute a sensitivity array for the AUC of each species -
species_index_array = [3 5]
number_of_species = length(species_index_array)
number_of_parameters = 11
results_array = zeros(number_of_parameters,1)
for species_index in species_index_array

    global results_array

    # conduct senstivity analysis -
    sensitivity_results = main(path_to_ensemble_file,species_index)

    # get the μ and σ^2
    mu = sensitivity_results.means
    var = sensitivity_results.variances

    #@show mu, var

    results_array = [results_array transpose(mu) transpose(var)]

end

results_array = results_array[:,2:end]

max_stars_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_issues_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_forks_repo_name: varnerlab/Biophysical-TXTL-Model-Code


text:
# exponentiation rapide

function expo(x, n)
    # n >= 1
    if n == 1
        x
    else
        y = expo(x, div(n, 2))
        y = mult(y, y)
        if n % 2 == 0 y else mult(x, y) end
    end
end

function mult(x::Int64, y::Int64) x * y end

println(expo(2, 1))
println(expo(2, 10))

# fibonacci en temps O(log n)

struct Mat22 a11; a12; a21; a22 end

function mult(x::Mat22, y::Mat22)
    Mat22(x.a11 * y.a11 + x.a12 * y.a21, x.a11 * y.a12 + x.a12 * y.a22,
          x.a21 * y.a11 + x.a22 * y.a21, x.a21 * y.a12 + x.a22 * y.a22)
end

fibM = Mat22(1, 1, 1, 0)
function fib(n)
    m = expo(fibM, n)
    m.a12
end

for n = 1:10 println(fib(n)) end

println(fib(1000000000))  # avec un débordement, bien sûr, mais instantané

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using LinearAlgebra

export calcmoit, eigenmoit

function calcmoit(atoms)
    # We center the atoms at their center of mass
    Atoms = Molecules.translate(atoms, Molecules.center_of_mass(atoms))
    I = zeros(Float64, (3,3))
    len = size(Atoms)[1]
    for i = 1:3, j = 1:3
        if i == j
            a = i % 3 + 1
            b = (i+1) % 3 + 1
            suum = 0.0
            for k = 1:len
                suum += Atoms[k].mass * (Atoms[k].xyz[a]^2 + Atoms[k].xyz[b]^2)
            end
            I[i,i] = suum
        else
            suum = 0.0
            for k = 1:len
                suum -= Atoms[k].mass * Atoms[k].xyz[i] * Atoms[k].xyz[j]
            end
            I[i,j] = suum
            I[j,i] = I[i,j]
        end
    end
    return I
end

function eigenmoit(moit)
    A = LinearAlgebra.Symmetric(moit)
    eval, evec = LinearAlgebra.eigen(A)
    return (eval, evec)
end


text:
using Statistics, Memoization

struct PeriodicFiniteDifferenceOperator{T} <: AbstractMatrix{T}
  N::Int
  difforder::Int
  accuracyorder::Int
  stencil::Vector{T}
end
function PeriodicFiniteDifferenceOperator(N::Int, L::T, difforder::Int, accuracyorder::Int=2) where {T}
  @assert iseven(accuracyorder)
  Δ = L / N
  b = zeros(accuracyorder + 1)
  b[difforder + 1] = factorial(difforder)
  x = (collect(1:accuracyorder + 1) .- (accuracy÷2 + 1)) .* Δ
  stencil = hcat([x[i+1].^(0:accuracyorder) for i in 0:accuracyorder]...) \ b ./ Δ^difforder
  return PeriodicFiniteDifferenceOperator{T}(N, difforder, accuracyorder, stencil)
end
Base.eltype(p::PeriodicFiniteDifferenceOperator{T}) where T = T
Base.size(p::PeriodicFiniteDifferenceOperator{T}) where {T} = (p.N, p.N)
function Base.getindex(p::PeriodicFiniteDifferenceOperator{T},
                       I::Vararg{Int, 2}) where {T}
  i, j = I
  A = p.order
  for jj ∈ (-A÷2:A÷2)
    j == mod1(jj + i, p.N) && return p.stencil[jj + A÷2 + 1]
  end
  return zero(T)
end

struct PeriodicFiniteIntegratorOperator{A,T} <: Function
  N::Int
  Δ::T
  function PeriodicFiniteIntegratorOperator(N::Int, L::T, order::Int=2) where {T}
    Δ = L / N
    return new{order,T}(N, Δ)
  end
end
@memoize polymatrix(N) = lu(hcat([(0:N-1).^i for i in 0:N-1]...))
function integratepoly(y, a, b) # TODO calculate coefficients and not do this
  N = length(y)
  N > 8 && throw(error("length of y is too large, gives badly conditioned matrix"))
  coeffs = zeros(N+1)
  coeffs[2:end] = (polymatrix(N) \ y) ./ (1:N)
  return sum(coeffs .* b.^(0:N)) - sum(coeffs .* a.^(0:N))
end
function (p::PeriodicFiniteIntegratorOperator{1})(z, y)
  z[1] = y[1]/2 * p.Δ
  for i ∈ 2:length(z)
    z[i] = z[i-1] + y[i] * p.Δ
  end
  return z
end
function (p::PeriodicFiniteIntegratorOperator{2})(z, y)
  z[1] = (y[end]/4 + 3y[1]/4) * p.Δ
  for i in 2:length(z)
    z[i] = z[i-1] + (y[i-1] + y[i])/2 * p.Δ
  end
  return z
end
# this operator is 4th order accurate despite using a 3 point stencil
function (p::PeriodicFiniteIntegratorOperator{3})(z, y)
  fill!(z, 0)
  z[1] = integratepoly([y[end], y[1], y[2]], 0.5, 1) * p.Δ
  z[2] = integratepoly([y[end], y[1], y[2]], 1, 1.5) * p.Δ
  for i in 2:length(z)-1
    z[i] += z[i-1] + integratepoly(y[i-1:i+1], 0.5, 1) * p.Δ
    z[i+1] += integratepoly(y[i-1:i+1], 1, 1.5) * p.Δ
  end
  z[end] += z[end-1] + integratepoly([y[end-1], y[end], y[1]], 0.5, 1) * p.Δ
  return z
end
function (p::PeriodicFiniteIntegratorOperator{4})(z, y)
  z[1] = integratepoly([y[end-1], y[end], y[1], y[2]], 1.5, 2) * p.Δ
  z[2] = z[1] + integratepoly([y[end], y[1], y[2], y[3]], 1, 2) * p.Δ
  for i in 3:length(z)-1
    z[i] = z[i-1] + integratepoly(y[i-2:i+1], 1, 2) * p.Δ
  end
  z[end] = z[end-1] + integratepoly([y[end-2], y[end-1], y[end], y[1]], 1, 2) * p.Δ
  return z
end
(p::PeriodicFiniteIntegratorOperator{A,T})(y) where {A,T} = p(deepcopy(y), y)

struct FiniteDifferenceField{BC<:PeriodicGridBC, T, A} <: AbstractField{BC}
  chargedensity::EquispacedValueGrid{BC,T}
  electricfield::EquispacedValueGrid{BC,T}
  integrator::PeriodicFiniteIntegratorOperator{A,T}
  function FiniteDifferenceField(chargedensity::EquispacedValueGrid{BC,T},
                                 electricfield::EquispacedValueGrid{BC,T},
                                 integrator::PeriodicFiniteIntegratorOperator{A,T}) where {BC,T,A}

    zero!(electricfield)
    finitediff = new{BC,T,A}(chargedensity, electricfield, integrator)
    solve!(finitediff)
    return finitediff
  end
end
function FiniteDifferenceField(chargedensity::EquispacedValueGrid{PeriodicGridBC};
                               order::Int=2)
  electricfield = deepcopy(chargedensity)
  zero!(electricfield)
  N = first(size(chargedensity))
  integrator = PeriodicFiniteIntegratorOperator(N, chargedensity.L, order)
  return FiniteDifferenceField(chargedensity, electricfield, integrator)
end
function FiniteDifferenceField(N::Int, L::Real, ::Type{BC}=PeriodicGridBC; order::Int=2) where {BC}
  return FiniteDifferenceField(EquispacedValueGrid(N, L, BC), order=order)
end

presolve(g::EquispacedValueGrid{PeriodicGridBC}) = g .- mean(g)
presolve(g::EquispacedValueGrid{<:AbstractBC}) = g
function solve!(f::FiniteDifferenceField{PeriodicGridBC})
  ρ = presolve(f.chargedensity)
  f.electricfield .= demean!(f.integrator(ρ))
  return f
end

zero!(f::FiniteDifferenceField) = (zero!(f.chargedensity); zero!(f.electricfield))


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    HierarchicalDensity <: DistLikeDensity

A hierarchical density, useful for hierarchical models/priors.

Constructors:

```julia
HierarchicalDensity(f::Function, parent_density::DistLikeDensity)
```

with a functon `f` that returns a `DistLikeDensity` for any variate `v` drawn
from `parent_dist`.

`Distribution`s are automatically converted to `DistLikeDensity`s.
`varshape(parent_density)` and `varshape(f(v))` must be a `NamedTupleShape`.

Example:

```julia
hd = HierarchicalDensity(
    v -> NamedTupleDist(
        baz = fill(Normal(v.bar, v.foo), 3)
    ),
    NamedTupleDist(
        foo = Exponential(3.5),
        bar = Normal(2.0, 1.0)
    )
)

varshape(hd) == NamedTupleShape(
    foo = ScalarShape{Real}(),
    bar = ScalarShape{Real}(),
    baz = ArrayShape{Real}(3)
)

v = rand(sampler(hd))
)
```
"""
struct HierarchicalDensity{
    F<:Function,
    D<:DistLikeDensity,
    S<:AbstractValueShape
} <: DistLikeDensity
    f::F
    pd::D
    vs::S
end

export HierarchicalDensity

function HierarchicalDensity(f::Function, pd::Any)
    pd_conv = convert(DistLikeDensity, pd)
    vs_pd = varshape(pd_conv)
    v_pd = rand(bat_determ_rng(), sampler(pd_conv))
    cd = _hd_cd(f, vs_pd, v_pd)
    vs = NamedTupleShape(;vs_pd..., varshape(cd)...)
    HierarchicalDensity(f, pd_conv, vs)
end

function _hd_cd(f::Function, vs_pd::AbstractValueShape, v_pd::AbstractVector{<:Real})
    shaped_v = stripscalar(vs_pd(v_pd))
    convert(AbstractDensity, f(shaped_v))
end

_hd_cd(f::Function, vs_pd::AbstractValueShape, v_pd::Any) = convert(AbstractDensity, f(v_pd))

function _hd_cd(d::HierarchicalDensity, v_pd::Any)
    vs_pd = varshape(d.pd)
    _hd_cd(d.f, vs_pd, v_pd)
end

_nt_type_names(::Type{<:NamedTuple{names}}) where names = names

@inline @generated function _split_nt(nt::NamedTuple, ::Val{names}) where {names}
    all_names = _nt_type_names(nt)
    rest_names = filter(n -> !(n in names), [all_names...])
    expr1 = Expr(:tuple, map(key -> :($key = nt.$key), names)...)
    expr2 = Expr(:tuple, map(key -> :($key = nt.$key), rest_names)...)
    Expr(:tuple, expr1, expr2)
end

@inline _split_nt(nt::NamedTuple, ::NamedTupleShape{names}) where {names} =
    _split_nt(nt, Val(names))

@inline function _split_v(d::HierarchicalDensity, v::NamedTuple)
    vsp = varshape(d.pd)
    _split_nt(v, vsp)
end

function _split_v(d::HierarchicalDensity, v::AbstractVector{<:Real})
    np = totalndof(d.pd)
    idxs = eachindex(v)
    idxs_a = first(idxs):(first(idxs) + np - 1)
    idxs_b = (first(idxs) + np):last(idxs)
    (view(v, idxs_a), view(v, idxs_b))
end

function eval_logval_unchecked(
    density::HierarchicalDensity,
    v::Any
)
    d = density
    v1, v2 = _split_v(d, v)
    logval1 = eval_logval_unchecked(d.pd, v1)
    cd = _hd_cd(d, v1)
    logval2 = eval_logval_unchecked(cd, v2)
    logval1 + logval2
end

ValueShapes.varshape(density::HierarchicalDensity) = density.vs

ValueShapes.totalndof(density::HierarchicalDensity) = totalndof(density.vs)

function Statistics.cov(density::HierarchicalDensity)
    cov(nestedview(rand(bat_determ_rng(), sampler(density), 10^5)))
end

Distributions.sampler(density::HierarchicalDensity) = HierarchicalDensitySampler(density)

var_bounds(density::HierarchicalDensity) = HierarchicalDensityBounds(density)

struct HierarchicalDensitySampler{D<:HierarchicalDensity} <: Sampleable{Multivariate,Continuous}
    d::D
end

function Distributions._rand!(rng::AbstractRNG, s::HierarchicalDensitySampler, v::AbstractVector{<:Real})
    d = s.d
    v1, v2 = _split_v(d, v)
    rand!(rng, sampler(d.pd),  v1)
    cd = _hd_cd(d, v1)
    rand!(rng, sampler(cd),  v2)
    v
end

Base.length(s::HierarchicalDensitySampler) = totalndof(varshape(s.d))

struct HierarchicalDensityBounds{D<:HierarchicalDensity} <: AbstractVarBounds
    d::D
end

ValueShapes.totalndof(bounds::HierarchicalDensityBounds) = totalndof(bounds.d)

function Base.eltype(bounds::HierarchicalDensityBounds)
    d = bounds.d
    v = rand(bat_determ_rng(), sampler(d))
    bounds1 = var_bounds(d.pd)
    v1, v2 = _split_v(d, v)
    cd = _hd_cd(d, v1)
    bounds2 = var_bounds(cd)
    promote_type(eltype(bounds1), eltype(bounds2))
end

function Base.in(v::AbstractVector{<:Real}, bounds::HierarchicalDensityBounds)
    d = bounds.d
    bounds1 = var_bounds(d.pd)
    v1, v2 = _split_v(d, v)
    if v1 in bounds1
        cd = _hd_cd(d, v1)
        bounds2 = var_bounds(cd)
        v2 in bounds2
    else
        false
    end
end

Base.in(v::Any, bounds::HierarchicalDensityBounds) = unshaped(v) in bounds

function renormalize_variate!(v_renorm::AbstractVector{<:Real}, bounds::HierarchicalDensityBounds, v::AbstractVector{<:Real})
    d = bounds.d
    bounds1 = var_bounds(d.pd)
    v1, v2 = _split_v(d, v)
    v1_renorm, v2_renorm = _split_v(d, v_renorm)
    renormalize_variate!(v1_renorm, bounds1, v1)

    new_v = if v1_renorm in bounds1
        cd = _hd_cd(d, v1)
        bounds2 = var_bounds(cd)
        renormalize_variate!(v2_renorm, bounds2, v2)
    end

    v_renorm
end