text:
module Showoff

using Dates

if isdefined(Base, :Ryu)
    include("ryu.jl")
else
    include("grisu.jl")
end

export showoff

# suppress compile errors when there isn't a grisu_ccall macro
macro grisu_ccall(x, mode, ndigits)
    quote end
end

# Fallback
function showoff(xs::AbstractArray, style=:none)
    result = Vector{String}(undef, length(xs))
    buf = IOBuffer()
    for (i, x) in enumerate(xs)
        show(buf, x)
        result[i] = String(take!(buf))
    end

    return result
end

# Floating-point

function concrete_minimum(xs)
    if isempty(xs)
        throw(ArgumentError("argument must not be empty"))
    end

    x_min = first(xs)
    for x in xs
        if isa(x, AbstractFloat) && isfinite(x)
            x_min = x
            break
        end
    end

    for x in xs
        if isa(x, AbstractFloat) && isfinite(x) && x < x_min
            x_min = x
        end
    end
    return x_min
end

function concrete_maximum(xs)
    if isempty(xs)
        throw(ArgumentError("argument must not be empty"))
    end

    x_max = first(xs)
    for x in xs
        if isa(x, AbstractFloat) && isfinite(x)
            x_max = x
            break
        end
    end

    for x in xs
        if isa(x, AbstractFloat) && isfinite(x) && x > x_max
            x_max = x
        end
    end
    return x_max
end

function scientific_precision_heuristic(xs::AbstractArray{<:AbstractFloat})
    ys = [x == 0.0 ? 0.0 : round(10.0 ^ (z = log10(abs(Float64(x))); z - floor(z)); sigdigits=15)
          for x in xs if isfinite(x)]
    return plain_precision_heuristic(ys) + 1
end

function showoff(xs::AbstractArray{<:AbstractFloat}, style=:auto)
    x_min = concrete_minimum(xs)
    x_max = concrete_maximum(xs)
    x_min = Float64(x_min)
    x_max = Float64(x_max)

    if !isfinite(x_min) || !isfinite(x_max)
        return invoke(showoff,Tuple{AbstractArray,Symbol},xs,:none)
    end

    if style == :auto
        if x_max != x_min && abs(log10(x_max - x_min)) > 4
            style = :scientific
        else
            style = :plain
        end
    end

    if style == :plain
        precision = plain_precision_heuristic(xs)
        return String[format_fixed(x, precision) for x in xs]
    elseif style == :scientific
        precision = scientific_precision_heuristic(xs)
        return String[format_fixed_scientific(x, precision, false)
                      for x in xs]
    elseif style == :engineering
        precision = scientific_precision_heuristic(xs)
        return String[format_fixed_scientific(x, precision, true)
                      for x in xs]
    else
        throw(ArgumentError("$(style) is not a recongnized number format"))
    end
end

const superscript_numerals = ['⁰', '¹', '²', '³', '⁴', '⁵', '⁶', '⁷', '⁸', '⁹']

function showoff(ds::AbstractArray{T}, style=:none) where T<:Union{Date,DateTime}
    years = Set()
    months = Set()
    days = Set()
    hours = Set()
    minutes = Set()
    seconds = Set()
    for d in ds
        push!(years, Dates.year(d))
        push!(months, Dates.month(d))
        push!(days, Dates.day(d))
        push!(hours, Dates.hour(d))
        push!(minutes, Dates.minute(d))
        push!(seconds, Dates.second(d))
    end
    all_same_year         = length(years)   == 1
    all_one_month         = length(months)  == 1 && 1 in months
    all_one_day           = length(days)    == 1 && 1 in days
    all_zero_hour         = length(hours)   == 1 && 0 in hours
    all_zero_minute       = length(minutes) == 1 && 0 in minutes
    all_zero_seconds      = length(minutes) == 1 && 0 in minutes
    all_zero_milliseconds = length(minutes) == 1 && 0 in minutes

    # first label format
    label_months = false
    label_days = false
    f1 = "u d, yyyy"
    f2 = ""
    if !all_zero_seconds
        f2 = "HH:MM:SS.sss"
    elseif !all_zero_seconds
        f2 = "HH:MM:SS"
    elseif !all_zero_hour || !all_zero_minute
        f2 = "HH:MM"
    else
        if !all_one_day
            first_label_format = "u d yyyy"
        elseif !all_one_month
            first_label_format = "u yyyy"
        elseif !all_one_day
            first_label_format = "yyyy"
        end
    end
    if f2 != ""
        first_label_format = string(f1, " ", f2)
    else
        first_label_format = f1
    end

    labels = Vector{String}(undef, length(ds))
    labels[1] = Dates.format(ds[1], first_label_format)
    d_last = ds[1]
    for (i, d) in enumerate(ds[2:end])
        if Dates.year(d) != Dates.year(d_last)
            if all_one_day && all_one_month
                f1 = "yyyy"
            elseif all_one_day && !all_one_month
                f1 = "u yyyy"
            else
                f1 = "u d, yyyy"
            end
        elseif Dates.month(d) != Dates.month(d_last)
            f1 = all_one_day ? "u" : "u d"
        elseif Dates.day(d) != Dates.day(d_last)
            f1 = "d"
        else
            f1 = ""
        end

        if f2 != ""
            f = string(f1, " ", f2)
        elseif f1 != ""
            f = f1
        else
            f = first_label_format
        end

        labels[i+1] = Dates.format(d, f)
        d_last = d
    end

    return labels
end

end # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ipx_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ipx")
JLLWrappers.@generate_main_file("ipx", UUID("b964385b-8358-5f44-856b-d218fe69d290"))
end  # module ipx_jll

ext: jl


text:
using StudentHousing, StatsBase, Distributions, JuMP, Gurobi

srand(32)

# =============================================================================
# Data
# =============================================================================
nbedrooms_range = collect(1:2)
nbedrooms_frequency = Weights([0.5, 0.5])
nbathrooms_range = collect(1:2)
nbathrooms_frequency = Weights([0.5, 0.5])
prices_range_pp = [800.0, 1000.0]
area_ranges = [0.0]
market_data = StudentHousing.MarketData(nbedrooms_range, nbedrooms_frequency,
    nbathrooms_range, nbathrooms_frequency, prices_range_pp, area_ranges)
budget = 1e6
problem_data = StudentHousingData(market_data, nhouses = 2, budget = budget,
                  demand_distribution = Uniform(0.9, 1.1))

d = problem_data
nhouses = length(d.houses)
npatterns = length(d.patterns_allow)

# =============================================================================
# Solve problem
# =============================================================================
m, U_generated, λ_generated, pattern_choice = solve_pattern_generation(d)

# =============================================================================
# Look at solution
# =============================================================================
getobjectivevalue(m)
# Recover the solution
# In case one of our initial columns was in the optimal basis
results = zeros(nhouses, length(d.patterns_allow))
for k = 1:npatterns
    if getvalue(m[:λ][k]) ≈ 1.0
        for i = 1:nhouses
            if house_fits_pattern(i, k, d)
                results[i, k] .= 1.0
                break
            end
        end
    end
end
# Check all the other columns we added
for k = 1:length(U_generated)
    if getvalue(λ_generated[k]) ≈ 1.0
        # Find the assignment that we picked for our house
        houses = find(U_generated[k] .≈ 1.0)
        # Check the pattern index
        p = pattern_choice[k]
        # Mark the assignment
        results[houses, p] .= 1.0
    end
end

println("House 1 assigned to:", find(results[1,:] .≈ 1.0))
# House 1 assigned to:[6, 7]

# sanity check
StudentHousing.house_allowedby(1, d)
# 5-element Array{Int64,1}:
#   6
#   7
#  11
#  12
#  13

println("House 2 assigned to:", find(results[2,:] .≈ 1.0))
# House 2 assigned to:[10]
# sanity check
StudentHousing.house_allowedby(2, d)
# 11-element Array{Int64,1}:
#   3
#   4
#   5
#   6
#   7
#   8
#   9
#  10
#  11
#  12
#  13


text:
# Borrowed from http://www.mathworks.com/help/signal/ref/upfirdn.html

import DSP: kaiser
import Multirate: firdes, filt
import PyPlot

sampleRate    = 48000
interpolation = 147
decimation    = 160
ratio         = interpolation  // decimation
Nt            = 24*interpolation

h = firdes( Nt, 0.5/interpolation, kaiser, beta = 7.8562  ) .* interpolation

n = 0:10239
x = sin.(2*pi*1e3/sampleRate*n)
y = filt( h, x, ratio )

tOriginal = n[1:49]./sampleRate
xOriginal = x[1:49]
tResamp   = n[1:45]./(sampleRate*ratio)
xResamp   = y[12:56]

PyPlot.stem( tOriginal, xOriginal, linefmt = "b-", markerfmt = "b." )
PyPlot.stem( tResamp, xResamp, linefmt = "r-", markerfmt = "r." )
PyPlot.show()

max_stars_repo_path: examples/DAT to CD.jl

max_issues_repo_path: examples/DAT to CD.jl

max_forks_repo_path: examples/DAT to CD.jl


text:
using RecipesBase
using DiffEqBase
import SciMLBase
import StaticArrays: SVector
using RecursiveArrayTools

abstract type AbstractTrajectory{U,T,A} <: AbstractVectorOfArray{U,2,A} end

"""
    Trajectory(u, t[, i])

A trajectory represented as a list of events.

# Fields

- `t`: Vector of event times.
- `u`: Vector of the same length as `t`. `u[i]` is the vector of copy numbers in the interval from `t[i-1]` to `t[i]`.
- `i`: A vector of event/reaction indices corresponding to the event times `t`. If empty, no event indices were recorded.

# Indexing

Indexing a trajectory with an integer index yields a tuple `(u, t, i)`. `u` is the trajectory value **before** the event at time `t`!
"""
struct Trajectory{U,T<:Real,A,B<:AbstractVector{T},C,D,E} <: AbstractTrajectory{U,T,A}
    u::VectorOfArray{U,2,A} # vector of length N. u[i] is the vector of copy numbers in the interval from t[i-1] to t[i]
    t::B # vector of length N, recording the reaction times and the end time of the trajectory
    i::C # vector of length N recording all reaction indices, or `nothing` if no reaction indices
    syms::D
    indepsym::E
end

Trajectory(u::VectorOfArray, t::Vector, i=nothing; syms=nothing, indepsym=nothing) = Trajectory(u, t, i, syms, indepsym)
Trajectory(u::Vector{<:AbstractVector}, t::Vector, i=nothing; syms=nothing, indepsym=nothing) = Trajectory(VectorOfArray(u), t, i; syms, indepsym)

function Base.copy(traj::Trajectory)
    Trajectory(
        copy(traj.u),
        deepcopy(traj.t),
        isnothing(traj.i) ? nothing : copy(traj.i),
        isnothing(traj.syms) ? nothing : copy(traj.syms),
        traj.indepsym
    )
end

Base.@propagate_inbounds Base.getindex(traj::AbstractTrajectory{U,T,A}, I::Int...) where {U,T,A} = traj.u[I[end]][Base.front(I)...]
Base.@propagate_inbounds Base.getindex(traj::AbstractTrajectory{U,T,A}, I::Int) where {U,T,A} = traj.u[I]
Base.@propagate_inbounds Base.getindex(traj::AbstractTrajectory{U,T,A}, I::Colon) where {U,T,A} = traj.u[I]
Base.@propagate_inbounds Base.getindex(traj::AbstractTrajectory{U,T,A}, I::AbstractArray{Int}) where {U,T,A} = Trajectory(traj.u[I], traj.t[I], isempty(traj.t) ? Int[] : traj.i[I])
# Base.@propagate_inbounds Base.getindex(traj::Trajectory, idx::Int) = (traj.u[idx], traj.t[idx], idx > length(traj.i) ? 0 : traj.i[idx])
Base.@propagate_inbounds function Base.getindex(traj::AbstractTrajectory{U,T,A},
    I::Union{Int,AbstractArray{Int},CartesianIndex,Colon,BitArray,AbstractArray{Bool}}...) where {U,T,A}
    traj.u[I...]
end

Base.@propagate_inbounds Base.getindex(t::AbstractTrajectory{U,T,A}, i::Int, ::Colon) where {U,T,A} = [t.u[j][i] for j in 1:length(t)]
Base.@propagate_inbounds Base.getindex(t::AbstractTrajectory{U,T,A}, ::Colon, i::Int) where {U,T,A} = t.u[i]
Base.@propagate_inbounds Base.getindex(t::AbstractTrajectory{U,T,A}, i::Int, II::AbstractArray{Int}) where {U,T,A} = [t.u[j][i] for j in II]

Base.@propagate_inbounds function Base.getindex(traj::AbstractTrajectory{U,T,A}, ii::CartesianIndex) where {U,T,A}
    ti = Tuple(ii)
    i = last(ti)
    jj = CartesianIndex(Base.front(ti))
    return traj.u[i][jj]
end

function equal_i(i1, i2)
    (i1 === nothing || isempty(i1)) && (i2 === nothing || isempty(i2)) || i1 == i2
end
Base.:(==)(traj1::AbstractTrajectory, traj2::AbstractTrajectory) = (traj1.t == traj2.t) && (traj1.u == traj2.u) && equal_i(traj1.i, traj2.i)

function (t::Trajectory)(time::Real)
    index = searchsortedfirst(t.t, time)
    if t.t[index] == time && index != lastindex(t.t)
        index += 1
    end
    if index > lastindex(t.t)
        error("Can't access trajectory that ends at t=$(last(t.t)) at time $time.")
    end
    t.u[index]
end

function (t::Trajectory)(times::AbstractArray{<:Real})
    dims = size(t.u)
    U = eltype(t.u)

    result = Array{U,2}(undef, (dims[1], length(times)))

    i = firstindex(t.u)
    for (j, dtime) in enumerate(times)
        while i < length(t.t) && dtime > t.t[i]
            i += 1
        end
        if i > length(t.u)
            result[:, j] .= t.u[end]
        else
            result[:, j] .= t.u[i]
        end
    end

    result
end

function Trajectory(u::AbstractMatrix{T}, t::AbstractVector{tType}, i=nothing) where {tType<:Real,T<:Real}
    num_components = size(u, 1)
    if num_components > 0
        u_vec = [SVector{num_components}(c) for c in eachcol(u)]
    else
        u_vec = SVector{0,T}[]
    end

    Trajectory(u_vec, convert(Vector{tType}, t), i)
end

function Trajectory(t::Trajectory)
    t
end

function Trajectory(sol::ODESolution)
    Trajectory(sol.u[begin:end-1], sol.t[begin+1:end])
end

Base.length(traj::AbstractTrajectory) = length(traj.t)
Base.firstindex(traj::AbstractTrajectory) = 1
Base.lastindex(traj::AbstractTrajectory) = length(traj)

function duration(traj::AbstractTrajectory)
    if length(traj) == 0
        return 0.0
    end
    traj.t[end] - traj.t[begin]
end

SciMLBase.getsyms(traj::AbstractTrajectory) = traj.syms
SciMLBase.getindepsym(traj::AbstractTrajectory) = traj.indepsym

struct TrajectoryTuplesIter{T}
    traj::T
end

RecursiveArrayTools.tuples(traj::AbstractTrajectory) = TrajectoryTuplesIter(traj)

function Base.iterate(iter::TrajectoryTuplesIter, index=1)
    traj = iter.traj
    index > length(traj) && return nothing
    (traj.u[index], traj.t[index], traj.i === nothing ? 0 : traj.i[index]), index + 1
end

Base.eltype(::Type{TrajectoryTuplesIter{T}}) where {uType,tType,A,T<:AbstractTrajectory{uType,tType,A}} = Tuple{eltype(A),tType,Int}
Base.length(iter::TrajectoryTuplesIter) = length(iter.traj)

SciMLBase.getsyms(iter::TrajectoryTuplesIter) = iter.traj.syms
SciMLBase.getindepsym(iter::TrajectoryTuplesIter) = iter.traj.indepsym

"""
    collect_trajectory(iter, nocopy=false)

Create a `Trajectory` from an iterator that yields a sequence of 3-element tuples `(u, t, i)`.

If `nocopy` is true, the `u` vectors of the tuple will not be copied before adding them to the
trajectory.
"""
function collect_trajectory(iter; nocopy=false)
    (u, t, i), state = iterate(iter)

    indepsym = SciMLBase.getindepsym(iter)
    syms = SciMLBase.getsyms(iter)

    traj = Trajectory([nocopy ? u : copy(u)], [t], Int[i], syms=syms, indepsym=indepsym)

    for (u, t, i) in Iterators.rest(iter, state)
        push!(traj.t, t)
        push!(traj.u, nocopy ? u : copy(u))
        push!(traj.i, i)
    end

    traj
end

function Base.show(io::IO, m::MIME"text/plain", x::AbstractTrajectory)
    print(io, "t: ")
    show(io, m, x.t)
    println(io)
    print(io, "u: ")
    show(io, m, x.u)
    println(io)
    if x.i !== nothing
        print(io, "i:")
        show(io, m, x.i)
    end
end

@recipe function f(traj::AbstractTrajectory{U,T,A}) where {U,T,A}
    seriestype --> :steppre
    xguide --> ((traj.indepsym !== nothing) ? string(traj.indepsym) : "")
    label --> ((traj.syms !== nothing) ? reshape(string.(traj.syms), 1, :) : "")
    traj.t, traj'
end

struct Chain{V,T<:AbstractVector,U<:AbstractVector} <: AbstractVector{V}
    head::T
    tail::U

    function Chain(head::AbstractVector{X}, tail::AbstractVector{Y}) where {X,Y}
        new{promote_type(X, Y),typeof(head),typeof(tail)}(head, tail)
    end
end

Base.IndexStyle(::Type{<:Chain}) = IndexLinear()
Base.copy(ch::Chain) = Chain(copy(ch.head), copy(ch.tail))
Base.size(ch::Chain) = size(ch.head) .+ size(ch.tail)
@inline Base.firstindex(ch::Chain) = 1
@inline Base.lastindex(ch::Chain) = length(ch)
Base.@propagate_inbounds function Base.getindex(ch::Chain{V}, i::Int) where {V}
    if i > length(ch.head)
        val = ch.tail[i-length(ch.head)]
        convert(V, val)
    else
        val = ch.head[i]
        convert(V, val)
    end
end
Base.setindex!(ch::Chain, v, i::Int) = i > length(ch.head) ? ch.tail[i-length(ch.head)] = v : ch.head[i] = v

struct MergeTrajectory{T1,T2}
    traj1::T1
    traj2::T2
end

Base.IteratorSize(::Type{MergeTrajectory{T1,T2}}) where {T1,T2} = Base.SizeUnknown()
function Base.eltype(::Type{MergeTrajectory{T1,T2}}) where {T1,T2}
    E1 = eltype(T1)
    E2 = eltype(T2)
    U1 = fieldtypes(E1)[1]
    U2 = fieldtypes(E2)[1]
    T = fieldtypes(E1)[2]
    Tuple{Chain{promote_type(eltype(U1), eltype(U2)),U1,U2},T,Int}
end

function merge_next(((u1, t1, i1), state1)::Tuple, ((u2, t2, i2), state2)::Tuple)
    t, i = if t1 <= t2
        t1, i1
    else
        t2, i2
    end
    Chain(copy(u1), copy(u2)), t, i
end

function advance_next(mtraj::MergeTrajectory, ((u1, t1, i1), state1)::A, ((u2, t2, i2), state2)::B)::Union{Nothing,Tuple{A,B}} where {A<:Tuple,B<:Tuple}
    if t1 < t2
        next1 = iterate(mtraj.traj1, state1)
        if next1 === nothing
            return nothing
        end
        next1, ((u2, t2, i2), state2)
    elseif t1 > t2
        next2 = iterate(mtraj.traj2, state2)
        if next2 === nothing
            return nothing
        end
        ((u1, t1, i1), state1), next2
    elseif t1 == t2
        next1 = iterate(mtraj.traj1, state1)
        next2 = iterate(mtraj.traj2, state2)
        (next1 === nothing || next2 === nothing) && return nothing
        next1, next2
    else # t1 or t2 is NaN
        nothing
    end
end

function Base.iterate(mtraj::MergeTrajectory)
    s1 = iterate(mtraj.traj1)
    s2 = iterate(mtraj.traj2)
    if s1 === nothing || s2 === nothing
        return nothing
    end
    m = merge_next(s1, s2)
    (m, (s1, s2))
end

function Base.iterate(mtraj::MergeTrajectory, (s1, s2)::Tuple)
    s = advance_next(mtraj, s1, s2)
    if s === nothing
        return nothing
    end
    s1, s2 = s
    m = merge_next(s1, s2)
    (m, s)
end

SciMLBase.getsyms(iter::MergeTrajectory) = nothing
SciMLBase.getindepsym(iter::MergeTrajectory) = nothing

trajectory_iterator(traj::AbstractTrajectory) = tuples(traj)
trajectory_iterator(traj) = traj

# merge more than 2 trajectories using recursion
merge_trajectories(traj) = trajectory_iterator(traj)
function merge_trajectories(traj1, traj2)
    MergeTrajectory(merge_trajectories(traj1), merge_trajectories(traj2))
end
function merge_trajectories(traj1, traj2, other_trajs...)
    merge12 = merge_trajectories(traj1, traj2)
    merge_trajectories(merge12, other_trajs...)
end

function start_collect(::Nothing, indices, syms, indepsym)
    # return empty trajectory
    Trajectory(Vector{Int16}[], Float64[], Int[], syms=syms, indepsym=indepsym)
end

function start_collect(((u, t, i), state)::Tuple, indices, syms, indepsym)
    Trajectory([u[indices]], typeof(t)[], typeof(i)[], syms=syms, indepsym=indepsym)
end

function step_collect!(result::Trajectory, ::Nothing, (uprev, tprev, iprev)::Tuple, indices)
    # end trajectory
    push!(result.t, tprev)
    push!(result.i, iprev)
    result
end

function step_collect!(
    result::Trajectory,
    ((u, t, i), state)::Tuple,
    (uprev, tprev, iprev)::Tuple,
    indices
)
    if (@view u[indices]) != result.u[end]
        push!(result.u, u[indices])
        push!(result.t, tprev)
        push!(result.i, iprev)
    end
    result
end

function collect_sub_trajectory(iter, indices)
    indepsym = SciMLBase.getindepsym(iter)
    syms = SciMLBase.getsyms(iter)

    sub_syms = isnothing(syms) ? nothing : syms[indices]

    f = iterate(iter)
    traj = start_collect(f, indices, sub_syms, indepsym)
    while f !== nothing
        val, state = f
        f = iterate(iter, state)
        traj = step_collect!(traj, f, val, indices)
    end
    traj
end

function collect_sub_trajectories(iter, indices_list...)
    indepsym = SciMLBase.getindepsym(iter)
    syms = SciMLBase.getsyms(iter)

    f = iterate(iter)
    trajs = map(indices -> start_collect(f, indices, isnothing(syms) ? nothing : syms[indices], indepsym), indices_list)
    while f !== nothing
        val, state = f
        f = iterate(iter, state)
        for (traj, indices) in zip(trajs, indices_list)
            step_collect!(traj, f, val, indices)
        end
    end
    trajs
end


text:
using OAuth, Test

@test typeof(oauth_timestamp()) <: AbstractString

@test parse(Int, oauth_timestamp()) > 1422235471

@test typeof(oauth_nonce(15)) <: AbstractString

@test length(oauth_nonce(15)) == 15
@test length(oauth_nonce(20)) == 20
@test length(oauth_nonce(25)) == 25
@test length(oauth_nonce(30)) == 30
@test length(oauth_nonce(32)) == 32

@test oauth_sign_hmac_sha1("randy", "zwitch") == "WqKCG5iyhFiES3fWYVdWJWwinaY="

@test oauth_signing_key("9djdj82h48djs9d2", "kkk9d7dh3k39sjv7") == "9djdj82h48djs9d2&kkk9d7dh3k39sjv7"

@test oauth_signature_base_string("POST", "http://example.com", "?julia=fast&lang=elegant") ==
        "POST&http%3A%2F%2Fexample.com&%3Fjulia%3Dfast%26lang%3Delegant"

@test oauth_percent_encode_keys!(Dict("badkey!"   => "value", "goodkey" => "value")) == Dict("badkey%21" => "value", "goodkey" => "value")

@test oauth_serialize_url_parameters(Dict("language" => "julia", "result" => "awesome")) == "language=julia&result=awesome"

@test encodeURI!(Dict("iv" => 10, "s" => "value!")) == Dict("iv" => 10, "s" => "value%21")

@test oauth_body_hash_encode("Hello, World!") == "CgqfKmdylCVXq1NV12r0Qvj2XgE="

@test oauth_body_hash_data("Hello, World!") == "oauth_body_hash=CgqfKmdylCVXq1NV12r0Qvj2XgE="

test_file = joinpath(dirname(@__FILE__), "auth_body_hash_file.txt")
@test oauth_body_hash_file(test_file) == "oauth_body_hash=CgqfKmdylCVXq1NV12r0Qvj2XgE="

#TODO
#oauth_header()

#TODO
#oauth_request_resource()


text:
using Snopt
using DelimitedFiles
using PyPlot
import ForwardDiff

using CSV
using DataFrames

# set up objective wrapper function
function aep_wrapper(x, params)
    # include relevant globals
    turbine_z = params.turbine_z
    turbine_x = params.turbine_x
    turbine_y = params.turbine_y
    rotor_diameter = params.rotor_diameter
    hub_height = params.hub_height
    ct_models = params.ct_models
    generator_efficiency = params.generator_efficiency
    cut_in_speed = params.cut_in_speed
    cut_out_speed = params.cut_out_speed
    rated_speed = params.rated_speed
    rated_power = params.rated_power
    windresource = params.windresource
    power_models = params.power_models
    model_set = params.model_set
    rotor_points_y = params.rotor_points_y
    rotor_points_z = params.rotor_points_z
    obj_scale = params.obj_scale

    # get number of turbines
    nturbines = length(x)

    # extract yaw of turbines from design variables vector
    turbine_yaw = x

    # calculate AEP
    AEP = obj_scale*ff.calculate_aep(turbine_x, turbine_y, turbine_z, rotor_diameter,
                hub_height, turbine_yaw, ct_models, generator_efficiency, cut_in_speed,
                cut_out_speed, rated_speed, rated_power, windresource, power_models, model_set,
                rotor_sample_points_y=rotor_points_y,rotor_sample_points_z=rotor_points_z, hours_per_year=365.0*24.0)

    # return the objective as an array
    return [AEP]
end

# set up optimization problem wrapper function
function wind_farm_opt(x)

    # calculate the objective function and jacobian (negative sign in order to maximize AEP)
    AEP = -aep_wrapper(x)[1]
    dAEP_dx = -ForwardDiff.jacobian(aep_wrapper,x)

    # set fail flag to false
    fail = false

    c = []
    dcdx = []

    # return objective, constraint, and jacobian values
    return AEP, c, dAEP_dx, dcdx, fail
end

# import model set with wind farm and related details
include("./model_sets/model_set_8_shiloh.jl")

# scale objective to be between 0 and 1
obj_scale = 1E-9

# set globals for use in wrapper functions
struct params_struct4{MS, AF, F, ACTM, WR, APM}
    model_set::MS
    rotor_points_y::AF
    rotor_points_z::AF
    turbine_x::AF
    turbine_y::AF
    turbine_z::AF
    rotor_diameter::AF
    obj_scale::F
    hub_height::AF
    ct_models::ACTM
    generator_efficiency::AF
    cut_in_speed::AF
    cut_out_speed::AF
    rated_speed::AF
    rated_power::AF
    windresource::WR
    power_models::APM
end

params = params_struct4(model_set, rotor_points_y, rotor_points_z, turbine_x, turbine_y, turbine_z,
    rotor_diameter, obj_scale, hub_height,
    ct_models, generator_efficiency, cut_in_speed, cut_out_speed, rated_speed, rated_power,
    windresource, power_models)

# initialize design variable array
x = zeros(nturbines)

# report initial objective value
println("Nturbines: ", nturbines)
println("Nstates: ", nstates)
println("Rotor diameter: ", rotor_diameter[1])
t1 = time()
println("Starting AEP value (GWh): ", aep_wrapper(x, params)*1e-9/obj_scale)
t2 = time()
println("AEP calc took ", t2-t1, "sec.")
# add initial turbine location to plot
for i = 1:length(turbine_x)
    plt.gcf().gca().add_artist(plt.Circle((turbine_x[i],turbine_y[i]), rotor_diameter[1]/2.0, fill=false,color="C0"))
end

# set up options for SNOPT
options = Dict{String, Any}()
options["Derivative option"] = 1
options["Verify level"] = 0
options["Major optimality tolerance"] = 1e-5
options["Major iteration limit"] = 1e6
options["Summary file"] = "snopt_yaw_summary.out"
options["Print file"] = "snopt_yaw_print.out"
lb = zeros(nturbines) .- 30.0*pi/180.0
ub = zeros(nturbines) .+ 30.0*pi/180.0

# generate wrapper function surrogate
aep_wrapper(x) = aep_wrapper(x, params)

initial_aep = aep_wrapper(x)

# set up holders for optimized yaw values for each state
optyaw = zeros((nstates,nturbines))
diraep = zeros(nstates)
diraepinit = zeros(nstates)
# run and time optimization
t1 = time()
for i in 1:nstates
    println("Optimizing for state: ", i, "/", nstates)
    println("Direction: ", winddirections[i])
    println("Wind speed: ", windspeeds[i])
    println("Probability: ", windprobabilities[i])
    options["Summary file"] = String("./yaw/snopt_yaw_summary_$i.out")
    options["Print file"] = String("./yaw/snopt_yaw_print_$i.out")
    params.windresource.wind_directions[1] = winddirections[i]
    params.windresource.wind_speeds[1] = windspeeds[i]
    params.windresource.wind_probabilities[1] = windprobabilities[i]
    params.windresource.measurement_heights[1] = measurementheight[i]
    params.windresource.ambient_tis[1] = ambient_tis[i]
    x = zeros(nturbines)
    diraepinit[i] = aep_wrapper(x)[1]
    t11 = time()
    xopt, fopt, info = snopt(wind_farm_opt, x, lb, ub, options)
    t22 = time()
    println("Finished in : ", t22-t11, " (s)")
    println("info: ", info)
    println("start AEP value for dir: ", diraepinit[i])
    println("end AEP value for dir: ", -fopt)
    println("AEP improvement for dir: ", -fopt/diraepinit[i])
    println()
    diraep[i] = deepcopy(fopt)
    for j in 1:nturbines
        optyaw[i,j] = deepcopy(xopt[j])
    end
end

t2 = time()
clkt = t1-t2

# print optimization results
println("Finished in : ", clkt, " (s)")
println("info: ", info)
println("start AEP value: ", sum(diraepinit))
println("end AEP value: ", sum(diraep))
println("optimized yaw", optyaw)
println("Total AEP improvement: ", sum(diraep)/sum(diraepinit))
# write results to csv files
dataforcsv_diraep = DataFrame(diraep = diraep)
CSV.write("./yaw/optyawdiraep.csv", dataforcsv_diraep)

dataforcsv_optyaw_dir = DataFrame(optyaw)
CSV.write("./yaw/optyaw.csv", dataforcsv_optyaw_dir)

# add final turbine locations to plot
for i = 1:length(turbine_x)
    plt.gcf().gca().add_artist(plt.Circle((turbine_x[i],turbine_y[i]), rotor_diameter[1]/2.0, fill=false,color="C1", linestyle="--"))
end

# set up and show plot
axis("square")
plt.show()


text:
"""
    Evolution

Evolutions define
```
    evolve(iter, value::T)::T
```
and possibly
```
    evolve(iter, key=>value)
```

They guarantee `HasEltype()` and `eltype(iter) == T`.
"""
abstract type Evolution <: DynamicIterator
end
const GEvolution = Union{Evolution, UnitRange, StepRange}
const GDynamicIterator = Union{DynamicIterator, UnitRange, StepRange}

evolve(r::UnitRange, i::Integer) = i < last(r) ?  i + 1 : nothing
evolve(r::UnitRange, ::Nothing) = r.start
function evolve(r::StepRange, i) # Fixme
    i = i + step(r)
    i <= last(r) ?  i : nothing
end
evolve(r::StepRange, ::Nothing) = r.start

@inline dyniterate(r::Union{UnitRange, StepRange}, ::Nothing) = iterate(r)
@inline dyniterate(r::Union{UnitRange, StepRange}, start::Start) = iterate(r, start.value)

dyniterate(E::Evolution, (start,)::Control{<:Union{Start,Value}}, control) = dubwith(evolve(E, start.value, control), Control)
dyniterate(E::GEvolution, value::Pair) = dub(evolve(E, value))
dyniterate(E::Evolution, (value,)::Control, control) = dubwith(evolve(E, value, control), Control)

IteratorSize(::Evolution) = SizeUnknown()

dyniterate(E::Evolution, start::Union{Value,Start}) = dub(evolve(E, start.value))
dyniterate(E::GEvolution, state) = dub(evolve(E, state))

"""
    evolve(f)

Create the DynamicIterator corresponding to the evolution
```
    x = f(x)
```

Integer keys default to increments.
Integer control default to keys (and repetitions).

```
julia> collect(take(from(Evolve(x->x + 1), 10), 5))
5-element Array{Any,1}:
 10
 11
 12
 13
 14
```
"""
struct Evolve{T} <: Evolution
    f::T
end

evolve(F::Evolve, x) = F.f(x)

evolve(F::Evolve, u::Pair, args...) = timelift_evolve(F, u, args...)

dyniterate(E::Evolve, (value,)::Control{<:Pair}, nextkey) = dubwith(evolve(E, value, nextkey), Control)

timelift_evolve(E, (i,x)::Pair) = i+1 => evolve(E, x)
function timelift_evolve(E, (i,x)::Pair{T}, j::T) where {T}
    @assert j ≥ i
    for k in 1:j-i
        x = evolve(E, x)
        x === nothing && return nothing
    end
    j => x
end


text:
using FileIO

struct LaxHeader
    version::UInt32
end

function Base.read(io::IO, ::Type{LaxHeader})
    header = LaxHeader(
        read(io, UInt32),
        )
end

struct LaxQuadtreeHeader
    ssignature::UInt32  # LASS
    stype::UInt32
    qsignature::UInt32  # LASQ
    version::UInt32
    levels::UInt32
    level_index::UInt32
    implicit_levels::UInt32
    min_x::Float64
    max_x::Float64
    min_y::Float64
    max_y::Float64
end

function Base.read(io::IO, ::Type{LaxQuadtreeHeader})
    header = LaxQuadtreeHeader(
        read(io, UInt32),
        read(io, UInt32),
        read(io, UInt32),
        read(io, UInt32),
        read(io, UInt32),
        read(io, UInt32),
        read(io, UInt32),
        Float64(read(io, Float32)),
        Float64(read(io, Float32)),
        Float64(read(io, Float32)),
        Float64(read(io, Float32))
    )
end

struct LaxIntervalHeader
    signature::UInt32
    version::UInt32
    number_cells::UInt32
end

function Base.read(io::IO, ::Type{LaxIntervalHeader})
    header = LaxIntervalHeader(
        read(io, UInt32),
        read(io, UInt32),
        read(io, UInt32),
    )
end

struct LaxIntervalCell
    cell_index::Int32
    number_intervals::UInt32
    number_points::UInt32
end

function Base.read(io::IO, ::Type{LaxIntervalCell})
    header = LaxIntervalCell(
        read(io, UInt32),
        read(io, UInt32),
        read(io, UInt32)
    )
end

struct LaxIntervalCellInterval
    _start::UInt32
    _end::UInt32
end

function Base.read(io::IO, ::Type{UnitRange{T}}) where T <: Integer
    range = UnitRange(
        Int64(read(io, UInt32))+1,  # .lax uses 0 based indices
        Int64(read(io, UInt32))+1   # Int32 could theoretically overflow
    )
end

function load(f::File{format"LAX"})
    open(f) do s
        skipmagic(s)
        load(s)
    end
end

function load(s::Stream{format"LAX"})
    # magic bytes are skipped
    header = read(s, LaxHeader)

    # read quadtree
    qheader = read(s, LaxQuadtreeHeader)
    # if not, qsignature and version are missing
    # but we ignore this other lax type for now
    qheader.level_index != 0 && error("Level index other than 0 is not supported yet.")

    # create quadtree
    qt = quadtree(qheader)

    # read cell header
    iheader = read(s, LaxIntervalHeader)
    ncells = iheader.number_cells

    # total number of points
    total_points = 0

    for i = 1:ncells
        # Read cell and its intervals
        qcell = read(s, LaxIntervalCell)
        total_points += qcell.number_points
        intervals = Vector{UnitRange{Integer}}(undef, qcell.number_intervals)

        for j = 1:qcell.number_intervals
            intervals[j] = read(s, UnitRange{Integer})
        end

        # Create cell in qt with intervals
        quadtree!(qt, qcell.cell_index, intervals)

    end

    # assert we are at end of file
    @assert eof(s)

    @info("Processed $(s.filename) with $(total_points) points.")

    return qt
end


text:
logistic(x) = 1 ./ (1 + exp.(-x))
@diffrule logistic(x::Number) x (logistic(x) .* (1 .- logistic(x))) .* ds

function autoencoder_cost(We1, We2, Wd, b1, b2, x)
    firstLayer = logistic.(We1 * x .+ b1)
    encodedInput = logistic.(We2 * firstLayer .+ b2)
    reconstructedInput = logistic.(Wd * encodedInput)
    cost = sum((reconstructedInput .- x) .^ 2.0)
    return cost
end

function mlp1(w1, w2, w3, x1)
    xx2 = w1 * x1
    x2 = log.(1. + exp.(xx2))
    xx3 = w2 * x2
    x3 = log.(1. + exp.(xx3))
    x4 = w3 * x3
    sum(1. ./ (1. + exp.(-x4)))
end

function mlp2(w1, w2, w3, b1, b2,  b3, x1)
    x2 = logistic.(w1 * x1 .+ b1)
    x3 = logistic.(w2 * x2 .+ b2)
    x4 = logistic.(w3 * x3 .+ b3)
    sum(x4)
end

function rnn(Wxh, Whh, Wy, hprev, x, y)
    h = tanh.(Whh * hprev + Wxh * x)
    yhat = Why * h
    cost = sum((yhat .- y) .^ 2.0)
    return h, cost
end


text:
using JuMP

function test(n::Int)
    m = Model()
    a = rand(n,n)
    @variable(m,x[1:n,1:n])
    b = rand(n,n,n)
    @variable(m,y[1:n,1:n,1:n])
    c = rand(n)
    @variable(m,z[1:n])

    #initialize
    @constraint(m,sum{c[i]*z[i],i=1:n}<=0)

    #Vector
    tic()
    @constraint(m,sum{c[i]*z[i],i=1:n}<=1)
    println("Vector with sum{}: $(toq())")
    tic()
    @constraint(m,vecdot(c,z) <= 1)
    println("Vector with vecdot() : $(toq())")

    #2D Matrix
    tic()
    @constraint(m,sum{a[i,j]*x[i,j],i=1:n,j=1:n}<=1)
    println("2D Matrix with sum{}: $(toq())")
    tic()
    @constraint(m,vecdot(a,x)<=1)
    println("2D Matrix with bigvecdot(): $(toq())")

    #3D Matrix
    tic()
    @constraint(m,sum{b[i,j,k]*y[i,j,k],i=1:n,j=1:n,k=1:n}<=1)
    println("3D Matrix with sum{}: $(toq())")
    tic()
    @constraint(m,vecdot(b,y)<=1)
    println("3D Matrix with vecdot(): $(toq())")
    return 0
end

for n in [10, 50, 100, 200, 300, 400]
    println("n = $n")
    test(n)
end


text:
"""Coefficients of the generalised logistic function to describe the vertical coordinate.
Default coefficients A,K,C,Q,B,M,ν are fitted to the old L31 configuration at ECMWF.
See geometry.jl and function vertical_coordinate for more informaiton.

Following the notation of https://en.wikipedia.org/wiki/Generalised_logistic_function (Dec 15 2021).

Change default parameters for more/fewer levels in the stratosphere vs troposphere vs boundary layer."""
@with_kw struct GenLogisticCoefs
    A::Real=-0.283     # obtained from a fit in /input_date/vertical_coordinate/vertical_resolution.ipynb
    K::Real= 0.871
    C::Real= 0.414
    Q::Real= 6.695
    B::Real=10.336
    M::Real= 0.602
    ν::Real= 5.812
end


text:
push!(LOAD_PATH,pwd())

  import HtmlConstructor

# --------- sample final html text
sample = """

  <meta charset = "utf-8" lang = "en-us" />

  <link rel="stylesheet" type="text/css" href="//code.jquery.com/ui/1.12.1/themes/base/jquery-ui.css" />

Testig HTML Page

<form style="width:80%;margin-left:10%;">

        Please answer <input name="q_1" type="text" required />

        Yes <input name="q_2" type="radio" value="true" checked />
        No <input name="q_2" type="radio" value="true" />

      <input type="submit" onclick="alert('test done!')" />

</form>

test 1

test 2

test 3

"""

# --------- init an html page
page = HtmlConstructor.HtmlPage("index.html","test page")

# ---------
# I.  part
# new an empty  tag
p_head = HtmlConstructor.quickEmptyPairedTag("head","")
# add members to p_head.content
  # 1. <meta> tag
  HtmlConstructor.add!(p_head, HtmlConstructor.tag_meta(charset="utf-8", lang="en-us") )
  # 2. <link> tag
  HtmlConstructor.add!(p_head, HtmlConstructor.link_css(
    "//code.jquery.com/ui/1.12.1/themes/base/jquery-ui.css"
  ))
  # 3.

max_stars_repo_name: Clpr/HtmlConstructor

max_issues_repo_name: Clpr/HtmlConstructor

max_forks_repo_name: Clpr/HtmlConstructor


text:
#=
Here are two approaches to deal with convolutions in the inverse fitting approach.
Taking the convolution of the IRF and the model appears to be computationally expensive.
Each iteration in the least squares fit has to compare the new fit with the data after convolving the model w/IRF.

Here I am trying to take advantage of the fact that the fft(IRF) is fixed for each step in the iteration and there is no reason to repeatedly calculate it.
The challenge is that the conv function from the FFTW library is very optimzed as it is a wrapper of a C/C++ library that is highly optimized.
Even still, there is an opportunity to plan out the FFT beforehand using the fact that the u and v in conv(u, v) are of the same length for our problem.
This allows to calculate the linear operator by the FFT needed to compute fft for a given set of dimensions.
In other words we can calculate a plan_fft (call this p) operator and then calculate the fft by u*p and v*p. If u is the IRF this value is constant for each iteration in the least squares.

=#

model(t, p) = p[1] * exp.(-p[2] *t)

tdata = 1:0.01:10

ydata = model(tdata, [1.0, 2.0]) + randn(length(tdata))

IRF = [Vector(1:500); Vector(500:-1:100)]

ydata = conv(ydata, IRF)

function convFIT(tdata, p, IRF)

    Rt = model(tdata, p)

    Rt = conv(Rt, IRF)
end

function fit(tdata, ydata, IRF)
    curve_fit((tdata, p) -> convFIT(tdata, p, IRF), tdata, ydata, [0.5, 0.5])
end

function convFIT1(tdata, p, pl, IRF_fft)
    Rt = model(tdata, p)
    dims = 2*length(Rt) -1

    Rt_fft = convert(Array{Complex{Float64},1}, [Rt ; zeros(Float64, dims-length(Rt))])
    r = similar(Rt_fft)
    i = similar(Rt_fft)
    c = similar(Rt_fft)

    mul!(r, pl, Rt_fft)
    mul!(i, pl, IRF_fft)

   ldiv!(c, pl, (r.*i))
   return real(c)
end

function fit2(tdata, ydata, IRF)

    dims = 2*length(IRF) -1

    IRF_fft = convert(Array{Complex{Float64},1}, [IRF ; zeros(Float64, dims-length(IRF))])

    pl = plan_fft(IRF_fft)

    curve_fit((tdata, p) -> convFIT1(tdata, p, pl, IRF_fft), tdata, ydata, [0.5, 0.5])
end

#=

############# Run these commands after running previous lines ################
julia> @benchmark f = fit(tdata,ydata,IRF)
BenchmarkTools.Trial:
  memory estimate:  132.05 MiB
  allocs estimate:  837915
  --------------
  minimum time:     792.865 ms (0.64% GC)
  median time:      801.574 ms (0.64% GC)
  mean time:        801.359 ms (0.64% GC)
  maximum time:     807.842 ms (0.63% GC)
  --------------
  samples:          7
  evals/sample:     1

 julia> @benchmark f2 = fit2(tdata,ydata,IRF)
  BenchmarkTools.Trial:
    memory estimate:  255.31 MiB
    allocs estimate:  1578872
    --------------
    minimum time:     212.654 ms (4.39% GC)
    median time:      214.032 ms (4.38% GC)
    mean time:        218.092 ms (4.41% GC)
    maximum time:     242.318 ms (3.91% GC)
    --------------
    samples:          23
    evals/sample:     1
  =#

  su = size(IRF)
  outsize = 2 .* su .-1
 # nffts = nextfastfft(outsize)

  upad = similar(IRF, outsize)
  vpad = similar(IRF, outsize)
  fill!(upad, zero(eltype(upad)))
  fill!(vpad, zero(eltype(vpad)))

  copyto!(upad, IRF)

  p! = plan_fft!(convert(Array{Complex{Float64}}, upad))
  upad = p! * upad

function conv_DT(t, data::input_data, upad, p!, vpad)

	RtDT = Array{Float64}(undef, length(0:0.01:10))
    RtDT = DT_model(0:0.01:10, [0.01,10.0], 1., 1., 1.)

    copyto!(vpad, RtDT)

    vpad = p! * vpad
    vpad .*= upad
    real(ifft!(vpad))

	Rt = Rt./maximum(Rt)

	tidx = findfirst(x -> x == t[1], data.t)

	convDT = convDT[tidx:tidx+length(t)-1]
	return log.(convDT)
end

function pad_fft(u, v)

    su = size(u)
    sv = size(v)
    outsize = su .+ sv .-1
    nffts = nextfastfft(outsize)

    upad = similar(u, nffts)
    vpad = similar(v, nffts)
    fill!(upad, zero(eltype(upad)))
    fill!(vpad, zero(eltype(vpad)))

    copyto!(upad, u)
    copyto!(vpad, v)

    p! = plan_fft!(convert(Array{Complex{Float64}}, upad))
    upad = p! * upad

function testing(upad, vpad, p!)

    vpad = p! * vpad
    vpad .*= upad
    real(ifft!(vpad))
end

outsize = 2 * length(data.IRF) .- 1
nffts = nextfastfft(outsize)

upad = zeros(Float64, nffts)
vpad = zeros(Float64, nffts)
convpad = zeros(Float64, nffts)

copyto!(upad, data.IRF)
pl = plan_rfft(upad)

upadRfft = zeros(Complex{Float64}, 1 .+ div(length(upad), 2, RoundDown))
upadRfft = pl * upad

RtDT = Array{Float64}(undef, length(data.IRF))

vpadRfft = zeros(Complex{Float64}, 1 .+ div(length(upad), 2, RoundDown))
u_vpadRfft = zeros(Complex{Float64}, 1 .+ div(length(upad), 2, RoundDown))

function conv_DT(t, β::Array{Float64,1}, data::input_data,
    vpad, RtDT, pl, vpadRfft, upadRfft, u_vpadRfft,
    convpad
    )

RtDT = DT_model(data.t, β, data.ρ, data.nmed, data.ndet)
copyto!(vpad, RtDT)

vpadRfft = pl * vpad
u_vpadRfft = upadRfft .* vpadRfft
convpad = irfft(u_vpadRfft, length(upad))
convpad = convpad./maximum(convpad)

tidx = findfirst(x -> x == t[1], data.t)

convpad = convpad[tidx:tidx+length(t)-1]
return log.(convpad)
end

#=@benchmark conv_DT(data.t[ind1:ind2], [0.1,10.], data, vpad, RtDT, pl, vpadRfft, upadRfft, u_vpadRfft, convpad)
BenchmarkTools.Trial:
memory estimate:  401.22 KiB
allocs estimate:  134
--------------
minimum time:     194.223 μs (0.00% GC)
median time:      199.783 μs (0.00% GC)
mean time:        241.692 μs (2.87% GC)
maximum time:     5.035 ms (0.00% GC)
--------------
samples:          10000
evals/sample:     1

old version
@benchmark conv_DT(data.t[ind1:ind2], [0.1,10.], data)
BenchmarkTools.Trial:
  memory estimate:  628.81 KiB
  allocs estimate:  184
  --------------
  minimum time:     296.018 μs (0.00% GC)
  median time:      308.628 μs (0.00% GC)
  mean time:        423.330 μs (2.68% GC)
  maximum time:     7.641 ms (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1
=#

function conv_DT(t, β::Array{Float64,1}, data::input_data,
    vpad, RtDT, pl, vpadRfft, upadRfft, u_vpadRfft,
    convpad
    )

RtDT = DT_model(data.t, β, data.ρ, data.nmed, data.ndet)
copyto!(vpad, RtDT)

vpadRfft = pl * vpad
u_vpadRfft = upadRfft .* vpadRfft
convpad = irfft(u_vpadRfft, length(upad))
convpad = convpad./maximum(convpad)

tidx = findfirst(x -> x == t[1], data.t)

convpad = convpad[tidx:tidx+length(t)-1]
return log.(convpad)
end

function getfit(input_data, model_params; alpha=0.1)

    ## get FFT padding and plans
    outsize = 2 * length(data.IRF) .- 1
    nffts = nextfastfft(outsize)

    upad = zeros(Float64, nffts)
    vpad = zeros(Float64, nffts)
    convpad = zeros(Float64, nffts)

    copyto!(upad, data.IRF)
    pl = plan_rfft(upad)

    upadRfft = zeros(Complex{Float64}, 1 .+ div(length(upad), 2, RoundDown))
    upadRfft = pl * upad

    RtDT = Array{Float64}(undef, length(data.IRF))

    vpadRfft = zeros(Complex{Float64}, 1 .+ div(length(upad), 2, RoundDown))
    u_vpadRfft = zeros(Complex{Float64}, 1 .+ div(length(upad), 2, RoundDown))

    #get fitting window
    counts = input_data.DTOF
    maxvalue, maxindex = findmax(counts)

    ind1 = findfirst(x -> x > maxvalue*model_params.risefactor, counts)
    ind2 =  findlast(x -> x >maxvalue*model_params.tailfactor, counts)

    weights = 1 ./ input_data.DTOFerrors

    fit = curve_fit((t, β) -> conv_DT(t, β, input_data, vpad, RtDT, pl, vpadRfft, upadRfft, u_vpadRfft, convpad), input_data.t[ind1:ind2], log.(input_data.DTOF[ind1:ind2]),
    model_params.initparams, lower=model_params.lb, upper=model_params.ub)

    perrors = estimate_errors(fit,alpha)
    marginerror = margin_error(fit,alpha)

    chisq = sum(fit.resid.^2)
    chisqdist = Distributions.Chisq(2)
    probability = Distributions.ccdf(chisqdist,chisq)
    return fitresult(input_data.t[ind1:ind2],log.(input_data.DTOF[ind1:ind2]),input_data.DTOFerrors, DT_model, fit.resid,fit.param,marginerror,perrors,2,chisq,probability)
end

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
# Deterministic method for 2D
function solveNFE(prob::ProbOutput2D,saveat)
    P      = prob.Plan
    Pinv   = prob.PlanInv
    krings = prob.krings
    rings  = prob.Ω.rings2D
    t      = prob.Ω.t
    x      = prob.Ω.x
    y      = prob.Ω.y
    N      = prob.Ω.N
    dt     = prob.Ω.dt
    α      = prob.α
    V      = prob.V0 # Initial condition in natural space (needed if saveat contains t[1])
    v      = prob.v0 # Initial condition in Fourier space
    sv     = prob.sv # Firing Rate in Fourier space (all delays)
    hN     = N÷2+1   # Half of dim N (due to the output of rfft)

    # Pre-allocate matrices
    Vj  = Matrix{Float64}(undef,N*length(saveat),N) # Stores the solution at 'saveat' instants
    svu = Matrix{ComplexF64}(undef,hN,N) # Fourier coefficients of the Firing Rate at delay u
    I   = similar(V)
    î   = similar(svu) # Fourier coefficients of the external input
    a   = similar(svu) # Fourier coefficients of the integral operator A(X,t)
    j   = 1 # Set index for solution saved at saveat instants

    @inbounds for i = 1:length(t) # Time loop
        ti = t[i]
        # Store solution V in saveat instants
        if ti in saveat
            Vj[(j-1)*N+1:N*j,:] .= V
            j += 1
        end

        I.= [prob.I(i,j,ti) for j in y, i in x] # Discretise I
        mul!(î,P,I)                              # rfft of I

        # Compute the integral operator's fourier coef at t_i
        @. a = @views krings[1:hN,:]*sv[1:hN,:] # Init a
        @inbounds for u = 2:rings # Rings loop
            @. a += @views krings[1+hN*(u-1):hN*u,:]*sv[1+hN*(u-1):hN*u,:]
        end

        # Euler explicit scheme: v_i+1 = v_i + (Δt/α)*(î_i - v_i + a_i)
        @. v = v + (dt/α)*(î - v + a)

        V .= Pinv * v # Perform inverse Fourier transform to compute V in natural domain

        # Update sv
        # Every delay block moves one block down. The last one is deleted
        sv[hN+1:end,:].= @view sv[1:hN*(rings-1),:]
        mul!(svu,P,prob.S.(V)) # Apply the Real Fourier Transform to V
        sv[1:hN,:].= svu # Update sv first delay ring
    end #end time loop

    # Build our solution structure output
    sol = SolveOutDet2D(Vj,prob.Ω.x,prob.Ω.y,prob.Ω.t,saveat,N)

    return sol

end

# Stochastic method for 2D
function solveNFE(prob::ProbOutput2D,saveat,ϵ,np,ξ=0.1)
    P      = prob.Plan
    Pinv   = prob.PlanInv
    krings = prob.krings
    rings  = prob.Ω.rings2D
    t      = prob.Ω.t
    x      = prob.Ω.x
    y      = prob.Ω.y
    N      = prob.Ω.N
    dt     = prob.Ω.dt
    α      = prob.α
    hN     = N÷2+1  # Half of dim N (due to the output of rfft)

    # Pre-allocate matrices
    meanVj = Matrix{Float64}(undef,N*length(saveat),N) # Store the mean solution at saveat
    Vj  = Array{Float64,3}(undef,N*length(saveat),N,np) # Store the solution at saveat at path p
    V   = Matrix{Float64}(undef,N,N)
    v   = Matrix{ComplexF64}(undef,hN,N)
    sv  = Matrix{ComplexF64}(undef,hN*rings,N) # Firing rate history in Fourier Space (matrix)
    svu = similar(v) # Fourier coefficients of the Firing Rate at delay u
    a   = similar(v) # Fourier coefficients of the integral operator A(X,t)
    î   = similar(v) # Fourier coefficients of the external input
    I   = similar(V)
    λ   = [exp(-((i^2+j^2)*ξ^2)/(8*pi)) for j = 1:hN, i = 1:N] # Correlation matrix

    # Pre-compute constant term of the stochastic part
    # fft gives the fourier coeff mult by N. Need to scale up noise, to have same magnitude
    # √dt due to the explicit Euler-Maruyama
    # Shift λ due to the also shifted output of fft
    c = (N^2)*sqrt(dt)*ϵ*fftshift(λ) # ϵ is the level of additive noise

    # Trajectories loop
    @inbounds for p = 1:np
        copy!(V,prob.V0) # Initial condition V(X,0) = V0 (Natural domain)
        copy!(v,prob.v0) # Initial condition v(X,0) = v0 (Fourier domain)
        copy!(sv,prob.sv)  # Initialise with the initial values of S

        j = 1 # Set index for tj
        @inbounds for i = 1:length(t) # Time loop
            # w ~ N(0,1) + i*N(0,1)  stochastic term
            w = randn(Float64,hN,N) + im*randn(Float64,hN,N)
            ti = t[i]
            # Store solution V in t[j] instant
            if ti in saveat
                Vj[(j-1)*N+1:N*j,:,p] .= V
                j += 1
            end

            I.= [prob.I(i,j,ti) for j in y, i in x] # Discretise I
            mul!(î,P,I)                             # rfft of I

            # Compute the integral operator's fourier coef at t_i
            @. a = @views krings[1:hN,:]*sv[1:hN,:] # Init a
            @inbounds for u = 2:rings # Ring loop
                @. a += @views krings[1+hN*(u-1):hN*u,:]*sv[1+hN*(u-1):hN*u,:]
            end

            # Euler explicit scheme: v_i+1 = v_i + (Δt/α)*(î_i - v_i + a_i) + √(Δt)*ϵ*λ*w
            @. v = v + (dt/α)*(î - v + a) + c*w

            # Apply inverse Fourier transform to compute V in natural domain
            V .= Pinv * v

            # Update sv
            # Every delay block moves one block down. The last one is deleted
            sv[hN+1:end,:] .= @view sv[1:hN*(rings-1),:]
            mul!(svu,P,prob.S.(V)) # Apply the Real Fourier Transform to V
            sv[1:hN,:].= svu # Update sv first delay ring
        end #end time loop
    end #end trajectories loop

    # Compute the mean solution across all trajectories
    @inbounds for j = 1:N
        @inbounds for i = 1:N*length(saveat)
            meanVj[i,j] = mean(@view Vj[i,j,:])
        end
    end

    # Build our solution structure output
    sol = SolveOutSto2D(Vj,meanVj,prob.Ω.x,prob.Ω.y,prob.Ω.t,saveat,N)

    return sol

end


text:
using PowerSimulationsDynamics
using Sundials

"""
Case 8:
This case study a 19-state virtual synchronous machine against an infinite bus located at bus 1, with VSM located at bus 2.
The perturbation increase the reference power (analogy for mechanical power) from 0.5 to 0.7.
"""

##################################################
############### LOAD DATA ########################
##################################################

include(joinpath(dirname(@__FILE__), "data_tests/test08.jl"))

##################################################
############### SOLVE PROBLEM ####################
##################################################

#time span
tspan = (0.0, 4.0);

#PSCAD benchmark data
csv_file = joinpath(dirname(@__FILE__), "benchmarks/pscad/Test08/Test08_omega.csv")
t_offset = 9.0

#Define Fault using Callbacks
Pref_change = ControlReferenceChange(1.0, case_inv, PSID.P_ref_index, 0.7)

path = (joinpath(pwd(), "test-08"))
!isdir(path) && mkdir(path)
try
    #Define Simulation Problem
    sim = Simulation!(
        ImplicitModel,
        omib_sys, # system
        path,
        tspan,
        Pref_change,
    )

    small_sig = small_signal_analysis(sim)

    #Solve problem in equilibrium
    execute!(sim, Sundials.IDA(), dtmax = 0.005, saveat = 0.005)

    #Obtain frequency data
    series = get_state_series(sim, ("generator-102-1", :ω_oc))
    t = series[1]
    ω = series[2]

    #Obtain PSCAD benchmark data
    M = get_csv_data(csv_file)
    t_pscad = M[:, 1] .- t_offset
    ω_pscad = M[:, 2]

    diff = [0.0]
    res = get_init_values_for_comparison(sim)
    for (k, v) in test08_x0_init
        diff[1] += LinearAlgebra.norm(res[k] - v)
    end
    @test (diff[1] < 1e-3)
    @test sim.solution.retcode == :Success
    @test small_sig.stable

    power = PSID.get_activepower_series(sim, "generator-102-1")
    rpower = PSID.get_reactivepower_series(sim, "generator-102-1")
    @test isa(power, Tuple{Vector{Float64}, Vector{Float64}})
    @test isa(rpower, Tuple{Vector{Float64}, Vector{Float64}})
    @test LinearAlgebra.norm(ω - ω_pscad) <= 1e-4
    @test LinearAlgebra.norm(t - round.(t_pscad, digits = 3)) == 0.0

finally
    @info("removing test files")
    rm(path, force = true, recursive = true)
end


text:
using Test
using TightBindingLattice

@testset "FractCoord" begin
    @testset "Explicit Construction" begin
        fractcoord = FractCoord([-2, 1], [0.8, 0.5])
        @test fractcoord.whole == [-2, 1]
        @test fractcoord.fraction == [0.8, 0.5]
    end

    @testset "Construction by Raw" begin
        fractcoord = FractCoord([-1.2, 1.5])
        @test fractcoord.whole == [-2, 1]
        @test isapprox(fractcoord.fraction, [0.8, 0.5])
    end

    @testset "Zero Construction" begin
        fractcoord = FractCoord(2)
        @test fractcoord.whole == [0,0]
        @test fractcoord.fraction == [0.0, 0.0]
    end

    @testset "methods" begin
        fractcoord = FractCoord([-2, 1], [0.8, 0.5])
        @test dimension(fractcoord) == 2
    end

    @testset "failures" begin
        @test_throws ArgumentError FractCoord(-1)
        @test_throws ArgumentError FractCoord(0)
        @test_throws ArgumentError FractCoord([-1, 1], [0.0])
        @test_throws ArgumentError FractCoord([-1, 1], [1.0, 2.0])
        @test_throws ArgumentError FractCoord([-1, 1], [0.0, 0.0]; tol=-1)
    end

    @testset "isapprox" begin
        fractcoord1 = FractCoord([-2, 1], [0.8, 0.5])
        fractcoord2 = FractCoord([-2, 1], [0.8, 0.5 + 1E-10])
        fractcoord3 = FractCoord([-1, 1], [0.8, 0.5])
        fractcoord4 = FractCoord([-1, 1], [0.8, 0.5+1E-2])

        @test isapprox(fractcoord1, fractcoord2)
        @test !isapprox(fractcoord1, fractcoord3)
        @test !isapprox(fractcoord1, fractcoord4)
    end

    @testset "round" begin
        fractcoord1 = FractCoord([1, 0], [0.0, 0.0])
        fractcoord2 = FractCoord([0, 0], [1.0-1E-12, 0.0]; tol=1E-6)
        @test isapprox(fractcoord1, fractcoord2)
    end

    @testset "Operators" begin
        lvs = [0.5 0.0; 0.0 1.0]
        fc1 = FractCoord([-2, 1], [0.8, 0.5])
        fc2 = FractCoord([1, 2], [0.1, 0.8])
        w = [3, 1]
        @test (+fc1) == fc1
        @test isapprox(+fc1, fc1)
        @test isapprox(-fc1, FractCoord([1, -2], [ 0.2, 0.5]))
        @test isapprox(fc1 + fc2, FractCoord([-1, 4], [0.9, 0.3]))
        @test isapprox(fc1 - fc2, FractCoord([-3, -2], [0.7, 0.7]))
        @test isapprox(fc1 + w, FractCoord([1, 2], [0.8, 0.5]))
        @test isapprox(fc1 - w, FractCoord([-5, 0], [0.8, 0.5]))
    end

    @testset "show" begin
        fc1 = FractCoord([-2, 1], [0.8, 0.5])
        s = string(fc1)
        @test string(fc1) == "FractCoord([-2, 1] + [0.8, 0.5])"
    end

    @testset "fract2carte exceptions" begin
        fc = FractCoord([-2, 1], [0.8, 0.5])

        @test_throws ArgumentError fract2carte([0.5 0.0;], fc)
        @test_throws ArgumentError fract2carte([0.5 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0], fc)
    end

    @testset "fract2carte/carte2fract" begin
        latticevectors = [0.5 0.0; 0.0 1.0]
        rawfractcoord = [-1.2, 1.5]

        correctfractcoord = FractCoord([-2, 1], [0.8, 0.5])
        correctcartecoord = [-0.6, 1.5]

        fractcoord = FractCoord(rawfractcoord)
        cartecoord = fract2carte(latticevectors, fractcoord)
        newfractcoord = carte2fract(latticevectors, cartecoord)
        @test isapprox(fractcoord, correctfractcoord)
        @test isapprox(cartecoord, correctcartecoord)
        @test isapprox(newfractcoord, correctfractcoord)

        @testset "tolerance" begin
            fc1 = carte2fract(latticevectors, [0.5 - 1E-9, 1.0])
            fc2 = carte2fract(latticevectors, [0.5 - 1E-9, 1.0]; tol=0.0)
            @test fc1.whole == [1, 1]
            @test isapprox(fc1.fraction, [0.0, 0.0])
            @test fc2.whole == [0, 1]
            @test isapprox(fc2.fraction, [1.0, 0.0])
        end
    end

    @testset "multiplication" begin
        fc = FractCoord([-2, 1], [0.8, 0.5])
        @test isapprox(FractCoord([1, -2], [0.5, 0.8]), [0 1; 1 0] * fc)
    end

end


text:
@testset "show" begin
    @testset "balance_widths" begin
        @test Tabulars.balance_widths(10, 10, 40) === (10, 10)
        @test Tabulars.balance_widths(20, 20, 40) === (20, 20)
        @test Tabulars.balance_widths(30, 20, 40) === (20, 20)
        @test Tabulars.balance_widths(30, 30, 40) === (20, 20)
        @test Tabulars.balance_widths(30, 30, 40) === (20, 20)
        @test Tabulars.balance_widths(10, 30, 40) === (10, 30)
        @test Tabulars.balance_widths(30, 10, 40) === (30, 10)
    end

    @testset "balance_widths!" begin
        widths = [22, 22, 22, 22]
        @test (n = Tabulars.balance_widths!(widths, 4, 100); (n === 4) & (widths == [22, 22, 22, 22]))

        widths = [23, 23, 23, 23]
        @test (n = Tabulars.balance_widths!(widths, 4, 100); (n === 4) & (widths == [23, 23, 23, 23]))

        widths = [24, 24, 24, 24]
        @test (n = Tabulars.balance_widths!(widths, 4, 100); (n === 4) & (widths == [23, 23, 23, 23]))

        widths = [100, 10, 10, 10]
        @test (n = Tabulars.balance_widths!(widths, 4, 80); (n === 4) & (widths == [42, 10, 10, 10]))

        widths = [10, 100, 10, 10]
        @test (n = Tabulars.balance_widths!(widths, 4, 80); (n === 4) & (widths == [10, 42, 10, 10]))

        widths = [10, 10, 100, 10]
        @test (n = Tabulars.balance_widths!(widths, 4, 80); (n === 4) & (widths == [10, 10, 42, 10]))

        widths = [10, 10, 10, 100]
        @test (n = Tabulars.balance_widths!(widths, 4, 80); (n === 4) & (widths == [10, 10, 10, 42]))

        widths = fill(10, 10)
        @test (n = Tabulars.balance_widths!(widths, 10, 80); (n === 6) & (widths == [10, 10, 10, 10, 10, 10, 10, 10, 10, 10]))
    end

    @testset "show Series" begin
        io = IOBuffer()
        show(io, MIME"text/plain"(), Series(l"a"=>1, l"b"=>2))
        @test String(io) == """
            2-element Series:
             a │ 1
             b │ 2"""

        io = IOBuffer()
        show(io, Series(zeros(100)))
        @test String(io) == """
            100-element Series:
             1  │ 0.0
             2  │ 0.0
             3  │ 0.0
             4  │ 0.0
             5  │ 0.0
             6  │ 0.0
             7  │ 0.0
             8  │ 0.0
             9  │ 0.0
             10 │ 0.0
             11 │ 0.0
             12 │ 0.0
             13 │ 0.0
             14 │ 0.0
             15 │ 0.0
             16 │ 0.0
             17 │ 0.0
             18 │ 0.0
             19 │ 0.0
             ⋮  │  ⋮"""

        io = IOBuffer()
        show(io, Series("This index is really, really long, in fact it's rather too long" => "The value is shorter by far"))
        @test String(io) == """
            1-element Series:
             This index is really, really long, in fact it's … │ The value is shorter by far"""
    end

    @testset "show Table" begin
        io = IOBuffer()
        show(io, MIME"text/plain"(), Table(l"a"=>[1,2,], l"b"=>[3,4]))
        @test String(io) == """
            2×2 Table:
                 a  b
               ┌─────
             1 │ 1  3
             2 │ 2  4"""

        io = IOBuffer()
        show(io, Table(zeros(100,100)))
        @test String(io) == """
            100×100 Table:
                  1    2    3    4    5    6    7    8    9    10   11   12   13   14   ⋯
                ┌────────────────────────────────────────────────────────────────────────
             1  │ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
             2  │ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
             3  │ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
             4  │ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
             5  │ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
             6  │ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
             7  │ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
             8  │ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
             9  │ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
             10 │ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
             11 │ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
             12 │ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
             13 │ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
             14 │ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
             15 │ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
             16 │ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
             17 │ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
             ⋮  │  ⋮    ⋮    ⋮    ⋮    ⋮    ⋮    ⋮    ⋮    ⋮    ⋮    ⋮    ⋮    ⋮    ⋮   ⋱"""

        io = IOBuffer()
        show(io, Table(l"Another cute string" => ["Another cute value"], l"This index is really, really long, in fact it's rather too much too long"=>["The value is shorter by far"]))
        @test String(io) == """
            1×2 Table:
                 Another cute string  This index is really, really long, in fact it's rathe…
               ┌────────────────────────────────────────────────────────────────────────────
             1 │ Another cute value   The value is shorter by far"""

    end
end


text:
using MathepiaInference
using Documenter

DocMeta.setdocmeta!(MathepiaInference, :DocTestSetup, :(using MathepiaInference); recursive=true)

makedocs(;
    modules=[MathepiaInference],
    authors="Pengfei Song",
    repo="https://github.com/Song921012/MathepiaInference.jl/blob/{commit}{path}#{line}",
    sitename="MathepiaInference.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Song921012.github.io/MathepiaInference.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Song921012/MathepiaInference.jl",
    devbranch="main",
)


text:
module EMLpi
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...types.EMLtypes
	import EMLtypes.length
	include("pi/PI_simple.jl")
	include("pi/PI.jl")
end


text:
# clique state machine for tree based initialization and inference

"""
    $SIGNATURES

Specialized info logger print function to show clique state machine information
in a standardized form.
"""
function infocsm(csmc::CliqStateMachineContainer, str::A) where {A <: AbstractString}

  tm = string(Dates.now())
  tmt = split(tm, 'T')[end]

  lbl = getLabel(csmc.cliq)
  lbl1 = split(lbl,',')[1]
  cliqst = getCliqStatus(csmc.cliq)

  with_logger(csmc.logger) do
    @info "$tmt | $(csmc.cliq.index)---$lbl1 @ $(cliqst) | "*str
  end
  flush(csmc.logger.stream)
  nothing
end

"""
    $SIGNATURES

Do cliq downward inference

Notes:
- State machine function nr. 11
"""
function doCliqDownSolve_StateMachine(csmc::CliqStateMachineContainer)
  infocsm(csmc, "11, doCliqDownSolve_StateMachine")
  setCliqDrawColor(csmc.cliq, "red")
  opts = getSolverParams(csmc.dfg)
  csmc.drawtree ? drawTree(csmc.tree, show=false, filepath=joinpath(opts.logpath,"bt.pdf")) : nothing

  # get down msg from parent (assing root clique CSM wont make it here)
  prnt = getParent(csmc.tree, csmc.cliq)
  dwnmsgs = getDwnMsgs(prnt[1])
  infocsm(csmc, "11, doCliqDownSolve_StateMachine -- dwnmsgs=$(collect(keys(dwnmsgs)))")

  # maybe cycle through separators (or better yet, just use values directly -- see next line)
  msgfcts = addMsgFactors!(csmc.cliqSubFg, dwnmsgs)
  # force separator variables to adopt down message values
  updateSubFgFromDownMsgs!(csmc.cliqSubFg, dwnmsgs, getCliqSeparatorVarIds(csmc.cliq))

  # add required all frontal connected factors
  newvars, newfcts = addDownVariableFactors!(csmc.dfg, csmc.cliqSubFg, csmc.cliq, csmc.logger, solvable=1)

  # store the cliqSubFg for later debugging
  if opts.dbg
    DFG.saveDFG(csmc.cliqSubFg, joinpath(opts.logpath,"logs/cliq$(csmc.cliq.index)/fg_beforedownsolve"))
    drawGraph(csmc.cliqSubFg, show=false, filepath=joinpath(opts.logpath,"logs/cliq$(csmc.cliq.index)/fg_beforedownsolve.pdf"))
  end

  ## new way
  # calculate belief on each of the frontal variables and iterate if required
  solveCliqDownFrontalProducts!(csmc.cliqSubFg, csmc.cliq, opts, csmc.logger)

  # compute new down messages
  infocsm(csmc, "11, doCliqDownSolve_StateMachine -- going to set new down msgs.")
  getSetDownMessagesComplete!(csmc.cliqSubFg, csmc.cliq, dwnmsgs, csmc.logger)
  # setDwnMsg!(cliq, drt.keepdwnmsgs)

      # frsyms = getCliqFrontalVarIds(csmc.cliq)
      # # use new localproduct approach
      # # TODO refactor with @spawn jl 1.3
      # if csmc.opts.multiproc
      #   downresult = Dict{Symbol, Tuple{BallTreeDensity, Float64, Vector{Symbol}}}()
      #   @sync for i in 1:length(frsyms)
      #     @async begin
      #       downresult[frsyms[i]] = remotecall_fetch(localProductAndUpdate!, upp2(), csmc.cliqSubFg, frsyms[i])
      #     end
      #   end
      #   infocsm(csmc, "doCliqDownSolve_StateMachine, multiproc keys=$(keys(downresult))")
      #   for fr in frsyms
      #     infocsm(csmc, "doCliqDownSolve_StateMachine, multiproc, setValKDE! $fr, inferdim=$(downresult[fr][2]), lbls=$(downresult[fr][3])")
      #     setValKDE!(csmc.cliqSubFg, fr, downresult[fr][1], false, downresult[fr][2])
      #   end
      # else
      #   for fr in frsyms
      #     localProductAndUpdate!(csmc.cliqSubFg, fr, csmc.logger)
      #   end
      # end

    ## old approach
      # # multiproc = false
      # # call down inference, TODO multiproc
      # if csmc.opts.multiproc
      #   cliqc = deepcopy(csmc.cliq)
      #   cliqcd = getData(cliqc)
      #   # redirect to new unused so that CAN be serialized
      #   cliqcd.initUpChannel = Channel{Symbol}(1)
      #   cliqcd.initDownChannel = Channel{Symbol}(1)
      #   cliqcd.solveCondition = Condition()
      #   cliqcd.statehistory = Vector{Tuple{DateTime, Int, Function, CliqStateMachineContainer}}()
      #   infocsm(csmc, "11, doCliqDownSolve_StateMachine -- MULTIPROC downGibbsCliqueDensity")
      #   # NOTE Cannot send logger to separate process memory
      #   drt = remotecall_fetch(downGibbsCliqueDensity, upp2(), csmc.cliqSubFg, cliqc, dwnmsgs, 100, 3, false, true)
      # else
      #   infocsm(csmc, "11, doCliqDownSolve_StateMachine -- SINGLEPROC downGibbsCliqueDensity")
      #   # usemsgpriors = true so that old plumbing is used properly
      #   drt = downGibbsCliqueDensity(csmc.cliqSubFg, csmc.cliq, dwnmsgs, 100, 3, false, true, csmc.logger)
      # end

      # # update clique subgraph with new status
      # updateFGBT!(csmc.cliqSubFg, csmc.tree, csmc.cliq.index, drt, dbg=false, fillcolor="lightblue", logger=csmc.logger)
      setCliqDrawColor(csmc.cliq, "lightblue")

  csmc.dodownsolve = false
  infocsm(csmc, "11, doCliqDownSolve_StateMachine -- finished with downGibbsCliqueDensity, now update csmc")

  # set MAP est
  for sym in getCliqFrontalVarIds(csmc.cliq)
    setVariablePosteriorEstimates!(csmc.cliqSubFg, sym)
  end

  # store the cliqSubFg for later debugging
  if opts.dbg
    DFG.saveDFG(csmc.cliqSubFg, joinpath(opts.logpath,"logs/cliq$(csmc.cliq.index)/fg_afterdownsolve"))
    drawGraph(csmc.cliqSubFg, show=false, filepath=joinpath(opts.logpath,"logs/cliq$(csmc.cliq.index)/fg_afterdownsolve.pdf"))
  end

  # transfer results to main factor graph
  frsyms = getCliqFrontalVarIds(csmc.cliq)
  infocsm(csmc, "11, finishingCliq -- going for transferUpdateSubGraph! on $frsyms")
  transferUpdateSubGraph!(csmc.dfg, csmc.cliqSubFg, frsyms, csmc.logger)

  # setCliqStatus!(csmc.cliq, :downsolved) # should be a notify
  infocsm(csmc, "11, doCliqDownSolve_StateMachine -- before notifyCliqDownInitStatus!")
  notifyCliqDownInitStatus!(csmc.cliq, :downsolved, logger=csmc.logger)
  infocsm(csmc, "11, doCliqDownSolve_StateMachine -- just notified notifyCliqDownInitStatus!")

  # remove msg factors that were added to the subfg
  infocsm(csmc, "11, doCliqDownSolve_StateMachine -- removing up message factors, length=$(length(msgfcts))")
  deleteMsgFactors!(csmc.cliqSubFg, msgfcts)

  # show tree if required
  csmc.drawtree ? drawTree(csmc.tree, show=false, filepath=joinpath(opts.logpath,"bt.pdf")) : nothing

  infocsm(csmc, "11, doCliqDownSolve_StateMachine -- finished, exiting CSM on clique=$(csmc.cliq.index)")

  # and finished
  return IncrementalInference.exitStateMachine
end

"""
    $SIGNATURES

Direct state machine to continue with downward solve or exit.

Notes
- State machine function nr. 10
"""
function determineCliqIfDownSolve_StateMachine(csmc::CliqStateMachineContainer)
  infocsm(csmc, "10, determineCliqIfDownSolve_StateMachine, csmc.dodownsolve=$(csmc.dodownsolve).")
  # finished and exit downsolve
  if !csmc.dodownsolve
    infocsm(csmc, "10, determineCliqIfDownSolve_StateMachine -- shortcut exit since downsolve not required.")
    return IncrementalInference.exitStateMachine
  end

  # yes, continue with downsolve
  setCliqDrawColor(csmc.cliq, "turquoise")
  csmc.drawtree ? drawTree(csmc.tree, show=false, filepath=joinpath(getSolverParams(csmc.dfg).logpath,"bt.pdf")) : nothing

  # assume separate down solve via solveCliq! call, but need a csmc.cliqSubFg this late in CSM anyway -- so just go copy one
  if length(ls(csmc.cliqSubFg)) == 0
    # first need to fetch cliq sub graph
    # go to 2b
    return buildCliqSubgraphForDown_StateMachine
  end

  # block here until parent is downsolved
  prnt = getParent(csmc.tree, csmc.cliq)
  if length(prnt) > 0
    infocsm(csmc, "10, determineCliqIfDownSolve_StateMachine, going to block on parent.")
    # TODO -- some cleanup
    blockCliqUntilParentDownSolved(prnt[1], logger=csmc.logger)
    prntst = getCliqStatus(prnt[1])
    infocsm(csmc, "10, determineCliqIfDownSolve_StateMachine, parent status=$prntst.")
    if prntst != :downsolved
      infocsm(csmc, "10, determineCliqIfDownSolve_StateMachine, going around again.")
      return determineCliqIfDownSolve_StateMachine
  	end
  else
    # special case for down solve on root clique.  When using solveCliq! following an up pass.

    # this is the root clique, so assume already downsolved -- only special case
    dwnmsgs = getCliqDownMsgsAfterDownSolve(csmc.cliqSubFg, csmc.cliq)
    setCliqDrawColor(csmc.cliq, "lightblue")
    setDwnMsg!(csmc.cliq, dwnmsgs)
    setCliqStatus!(csmc.cliq, :downsolved)
	csmc.dodownsolve = false

	# Update estimates and transfer back to the graph
	frsyms = getCliqFrontalVarIds(csmc.cliq)
	# Calculate estimates
	map(sym -> setVariablePosteriorEstimates!(csmc.cliqSubFg, sym), frsyms)
	# Transfer to parent graph
	transferUpdateSubGraph!(csmc.dfg, csmc.cliqSubFg, frsyms)

    notifyCliqDownInitStatus!(csmc.cliq, :downsolved, logger=csmc.logger)

    return IncrementalInference.exitStateMachine
  end

  infocsm(csmc, "10, going for down solve.")
  # go to 11
  return doCliqDownSolve_StateMachine
end

"""
    $SIGNATURES

Notes
- State machine function nr.9
"""
function finishCliqSolveCheck_StateMachine(csmc::CliqStateMachineContainer)
  csmc.drawtree ? drawTree(csmc.tree, show=false, filepath=joinpath(getSolverParams(csmc.dfg).logpath,"bt.pdf")) : nothing
  cliqst = getCliqStatus(csmc.cliq)
  infocsm(csmc, "9, finishingCliq")
  if cliqst == :upsolved
      frsyms = getCliqFrontalVarIds(csmc.cliq)
    infocsm(csmc, "9, finishingCliq -- going for transferUpdateSubGraph! on $frsyms")
    # TODO what about down solve??
    transferUpdateSubGraph!(csmc.dfg, csmc.cliqSubFg, frsyms, csmc.logger)

    # remove any solvable upward cached data -- TODO will have to be changed for long down partial chains
    # assuming maximally complte up solved cliq at this point
    lockUpStatus!(csmc.cliq)
    sdims = Dict{Symbol,Float64}()
    for varid in getCliqAllVarIds(csmc.cliq)
      sdims[varid] = 0.0
    end
    updateCliqSolvableDims!(csmc.cliq, sdims, csmc.logger)
    unlockUpStatus!(csmc.cliq)

    # go to 10
    return determineCliqIfDownSolve_StateMachine # IncrementalInference.exitStateMachine
  elseif cliqst == :initialized
    setCliqDrawColor(csmc.cliq, "sienna")
    csmc.drawtree ? drawTree(csmc.tree, show=false, filepath=joinpath(getSolverParams(csmc.dfg).logpath,"bt.pdf")) : nothing
    # go to 7
    return determineCliqNeedDownMsg_StateMachine
  else
    infocsm(csmc, "9, finishingCliq -- init not complete and should wait on init down message.")
    setCliqDrawColor(csmc.cliq, "coral")
    # TODO, potential problem with trying to downsolve
    # return isCliqNull_StateMachine # doesCliqNeeddownmsg_StateMachine
  end

  # go to 4
  return isCliqNull_StateMachine # whileCliqNotSolved_StateMachine
end

"""
    $SIGNATURES

Notes
- State machine function nr. 8c
"""
function waitChangeOnParentCondition_StateMachine(csmc::CliqStateMachineContainer)
  prnt = getParent(csmc.tree, csmc.cliq)
  if length(prnt) > 0
    infocsm(csmc, "8c, waitChangeOnParentCondition_StateMachine, wait on parent=$(prnt[1].index) for condition notify.")
    @sync begin
      @async begin
        sleep(1)
        notify(getSolveCondition(prnt[1]))
      end
      wait(getSolveCondition(prnt[1]))
    end
  else
    infocsm(csmc, "8c, waitChangeOnParentCondition_StateMachine, cannot wait on parent for condition notify.")
    @warn "no parent!"
  end

  # go to 4
  return isCliqNull_StateMachine
end

"""
    $SIGNATURES

Do up initialization calculations, loosely translates to solving Chapman-Kolmogorov
transit integral in upward direction.

Notes
- State machine function nr. 8b
- Includes initialization routines.
- TODO: Make multi-core
"""
function attemptCliqInitUp_StateMachine(csmc::CliqStateMachineContainer)
  # get a few pointers to relevant data
  cliqst = getCliqStatus(csmc.cliq)
  opts = getSolverParams(csmc.dfg)

  infocsm(csmc, "8b, attemptCliqInitUp, !areCliqChildrenNeedDownMsg()=$(!areCliqChildrenNeedDownMsg(csmc.tree, csmc.cliq))" )
  if cliqst in [:initialized; :null; :needdownmsg] && !areCliqChildrenNeedDownMsg(csmc.tree, csmc.cliq)
    setCliqDrawColor(csmc.cliq, "red")
    csmc.drawtree ? drawTree(csmc.tree, show=false, filepath=joinpath(opts.logpath,"bt.pdf")) : nothing

    # check if init is required and possible
    infocsm(csmc, "8b, attemptCliqInitUp, going for doCliqAutoInitUpPart1!.")
    # get incoming clique up messages
    upmsgs = getCliqInitUpMsgs(csmc.cliq)
    # add incoming up messages as priors to subfg
    infocsm(csmc, "8b, doCliqAutoInitUpPart1! -- adding up message factors")
    msgfcts = addMsgFactors!(csmc.cliqSubFg, upmsgs)

    # store the cliqSubFg for later debugging
    if opts.dbg
      DFG.saveDFG(csmc.cliqSubFg, joinpath(opts.logpath,"logs/cliq$(csmc.cliq.index)/fg_beforeupsolve"))
      drawGraph(csmc.cliqSubFg, show=false, filepath=joinpath(opts.logpath,"logs/cliq$(csmc.cliq.index)/fg_beforeupsolve.pdf"))
    end

    doCliqAutoInitUpPart1!(csmc.cliqSubFg, csmc.tree, csmc.cliq, logger=csmc.logger)
    infocsm(csmc, "8b, attemptCliqInitUp, areCliqVariablesAllInitialized(subfg, cliq)=$(areCliqVariablesAllInitialized(csmc.cliqSubFg, csmc.cliq))")

    # do actual up solve
    retstatus = doCliqAutoInitUpPart2!(csmc.cliqSubFg, csmc.tree, csmc.cliq, multiproc=csmc.opts.multiproc, logger=csmc.logger)

    # remove msg factors that were added to the subfg
    infocsm(csmc, "8b, doCliqAutoInitUpPart2! -- removing up message factors, length=$(length(msgfcts))")
    deleteMsgFactors!(csmc.cliqSubFg, msgfcts)

    # store the cliqSubFg for later debugging
    if opts.dbg
      DFG.saveDFG(csmc.cliqSubFg, joinpath(opts.logpath,"logs/cliq$(csmc.cliq.index)/fg_afterupsolve"))
      drawGraph(csmc.cliqSubFg, show=false, filepath=joinpath(opts.logpath,"logs/cliq$(csmc.cliq.index)/fg_afterupsolve.pdf"))
    end

    # notify of results
    if cliqst != retstatus
      infocsm(csmc, "8b, attemptCliqInitUp, post-doCliqAu. -- notification retstatus=$retstatus")
      notifyCliqUpInitStatus!(csmc.cliq, retstatus, logger=csmc.logger)
    else
      infocsm(csmc, "8b, attemptCliqInitUp, post-doCliqAu. -- no notification required $cliqst=$retstatus")
    end
  end

  # go to 9
  return finishCliqSolveCheck_StateMachine
end

"""
    $SIGNATURES

Do down initialization calculations, loosely translates to solving Chapman-Kolmogorov
transit integral in downward direction.

Notes
- State machine function nr. 8a
- Includes initialization routines.
- TODO: Make multi-core
"""
function attemptCliqInitDown_StateMachine(csmc::CliqStateMachineContainer)
  #
  # should never happen to
  setCliqDrawColor(csmc.cliq, "green")
  csmc.drawtree ? drawTree(csmc.tree, show=false, filepath=joinpath(getSolverParams(csmc.dfg).logpath,"bt.pdf")) : nothing

  # initialize clique in downward direction
  # not if parent also needs downward init message
  infocsm(csmc, "8a, needs down message -- attempt down init")
  prnt = getParent(csmc.tree, csmc.cliq)[1]

  # take atomic lock when waiting for down ward information
  lockUpStatus!(getData(prnt))

  dwinmsgs = prepCliqInitMsgsDown!(csmc.dfg, csmc.tree, prnt, csmc.cliq, logger=csmc.logger, dbgnew=!haskey(getSolverParams(csmc.dfg).devParams,:dontUseParentFactorsInitDown)) # csmc.cliqSubFg
  dwnkeys = collect(keys(dwinmsgs))

  infocsm(csmc, "8a, attemptCliqInitD., dwinmsgs=$(dwnkeys), adding msg factors")

  ## DEVIdea
  msgfcts = addMsgFactors!(csmc.cliqSubFg, dwinmsgs)

  # determine if more info is needed for partial

  # DEVidea
  sdims = getCliqVariableMoreInitDims(csmc.cliqSubFg, csmc.cliq)
  updateCliqSolvableDims!(csmc.cliq, sdims, csmc.logger)

  infocsm(csmc, "8a, attemptCliqInitD., updated clique solvable dims")

  # priorize solve order for mustinitdown with lowest dependency first
  # follow example from issue #344
  mustwait = false
  if length(intersect(dwnkeys, getCliqSeparatorVarIds(csmc.cliq))) == 0 # length(dwinmsgs) == 0 ||
    infocsm(csmc, "8a, attemptCliqInitDown_StateMachine, no can do, must wait for siblings to update parent first.")
    mustwait = true
  elseif getSiblingsDelayOrder(csmc.tree, csmc.cliq, prnt, dwinmsgs, logger=csmc.logger)
    infocsm(csmc, "8a, attemptCliqInitD., prioritize")
    mustwait = true
  elseif getCliqSiblingsPartialNeeds(csmc.tree, csmc.cliq, prnt, dwinmsgs, logger=csmc.logger)
    infocsm(csmc, "8a, attemptCliqInitD., partialneedsmore")
    mustwait = true
  end

  # remove the downward messages too
  deleteMsgFactors!(csmc.cliqSubFg, msgfcts)

  infocsm(csmc, "8a, attemptCliqInitD., deleted msg factors and unlockUpStatus!")
  # unlock
  unlockUpStatus!(getData(prnt))
  infocsm(csmc, "8a, attemptCliqInitD., unlocked")

  solord = getCliqSiblingsPriorityInitOrder( csmc.tree, prnt, csmc.logger )
  infocsm(csmc, "8a, attemptCliqInitD., before are Siblings")
  noOneElse = areSiblingsRemaingNeedDownOnly(csmc.tree, csmc.cliq)
  infocsm(csmc, "8a, attemptCliqInitD., $mustwait, $noOneElse, solord = $solord")

  if mustwait && csmc.cliq.index!=solord[1] # && !noOneElse
    # go to 8c
    infocsm(csmc, "8a, attemptCliqInitD., mustwait, so wait on change.")
    return waitChangeOnParentCondition_StateMachine
  end

  ## TODO deal with partial inits only, either delay or continue at end...
  # find intersect between downinitmsgs and local clique variables
  # if only partials available, then

  infocsm(csmc, "8a, attemptCliqInitD.,do cliq init down dwinmsgs=$(keys(dwinmsgs))")
  cliqst = doCliqInitDown!(csmc.cliqSubFg, csmc.cliq, dwinmsgs, dbg=getSolverParams(csmc.dfg).dbg, logger=csmc.logger, logpath=getSolverParams(csmc.dfg).logpath )
  # TODO: transfer values changed in the cliques should be transfered to the tree in proc 1 here.

  # # TODO: is status of notify required here?
  setCliqStatus!(csmc.cliq, cliqst)
  # notifyCliqUpInitStatus!(csmc.cliq, cliqst)

  # TODO move out
  children = getChildren(csmc.tree, csmc.cliq)
  if areCliqChildrenNeedDownMsg(children)
    # set messages if children :needdownmsg
    infocsm(csmc, "8a, doCliqInitDown! -- must set messages for future down init")
    # construct init's up msg to place in parent from initialized separator variables
    msg = prepCliqInitMsgsUp(csmc.cliqSubFg, csmc.cliq, csmc.logger) # , tree,

    infocsm(csmc, "8a, putting fake upinitmsg in this cliq, msgs labels $(collect(keys(msg)))")
    # set fake up and notify down status -- repeat change status to same as notifyUp above
    setCliqUpInitMsgs!(csmc.cliq, csmc.cliq.index, msg)
    # setCliqStatus!(csmc.cliq, cliqst)
    setCliqDrawColor(csmc.cliq, "sienna")
    csmc.drawtree ? drawTree(csmc.tree, show=false, filepath=joinpath(getSolverParams(csmc.dfg).logpath,"bt.pdf")) : nothing

    notifyCliqDownInitStatus!(csmc.cliq, cliqst, logger=csmc.logger)

    infocsm(csmc, "8a, near-end down init attempt, $cliqst.")
  end

  return attemptCliqInitUp_StateMachine
end

"""
    $SIGNATURES

Delay loop if waiting on upsolves to complete.

Notes
- State machine 7b
"""
function slowCliqIfChildrenNotUpsolved_StateMachine(csmc::CliqStateMachineContainer)

  # special case short cut
  cliqst = getCliqStatus(csmc.cliq)
  if cliqst == :needdownmsg
    infocsm(csmc, "7b, shortcut on cliq is a needdownmsg status.")
    return isCliqNull_StateMachine
  end
  childs = getChildren(csmc.tree, csmc.cliq)
  len = length(childs)
  @inbounds for i in 1:len
    if !(getCliqStatus(childs[i]) in [:upsolved;:uprecycled;:marginalized])
      infocsm(csmc, "7b, wait condition on upsolve, i=$i, ch_lbl=$(getCliqFrontalVarIds(childs[i])[1]).")
      wait(getSolveCondition(childs[i]))
      break
    end
  end

  # go to 4
  return isCliqNull_StateMachine
end

"""
    $SIGNATURES

Notes
- State machine function nr. 7
"""
function determineCliqNeedDownMsg_StateMachine(csmc::CliqStateMachineContainer)

  infocsm(csmc, "7, start, forceproceed=$(csmc.forceproceed)")

  # fetch children status
  stdict = blockCliqUntilChildrenHaveUpStatus(csmc.tree, csmc.cliq, csmc.logger)

  # hard assumption here on upsolve from leaves to root
  proceed = true
  # fetch status from children (should already be available -- i.e. should not block)
  for (clid, clst) in stdict
    infocsm(csmc, "7, check stdict children: clid=$(clid), clst=$(clst)")
    # :needdownmsg # 'send' downward init msg direction
    !(clst in [:initialized;:upsolved;:marginalized;:downsolved;:uprecycled]) ? (proceed = false) : nothing
  end
  infocsm(csmc, "7, proceed=$(proceed)")

  if proceed || csmc.forceproceed
    # TODO, remove csmc.forceproceed
    csmc.forceproceed = false
    # return doCliqInferAttempt_StateMachine
    cliqst = getCliqStatus(csmc.cliq)
    infocsm(csmc, "7, status=$(cliqst), before attemptCliqInitDown_StateMachine")
    # d1,d2,cliqst = doCliqInitUpOrDown!(csmc.cliqSubFg, csmc.tree, csmc.cliq, isprntnddw)
    if cliqst == :needdownmsg && !isCliqParentNeedDownMsg(csmc.tree, csmc.cliq, csmc.logger)
      # go to 8a
      return attemptCliqInitDown_StateMachine
    # HALF DUPLICATED IN STEP 4
    elseif cliqst == :marginalized
      # go to 1
      return isCliqUpSolved_StateMachine
      ## NOTE -- what about notifyCliqUpInitStatus! ??
      # go to 10
      # return determineCliqIfDownSolve_StateMachine
    end

    # go to 8b
    return attemptCliqInitUp_StateMachine
  else
    # go to 7b
    return slowCliqIfChildrenNotUpsolved_StateMachine
  end
end

"""
    $SIGNATURES

Notes
- State machine function nr. 6c
"""
function blockCliqSiblingsParentChildrenNeedDown_StateMachine(csmc::CliqStateMachineContainer)
  # add blocking case when all siblings and parent :needdownmsg -- until parent :initialized
  infocsm(csmc, "6c, check/block sibl&prnt :needdownmsg")
  blockCliqSiblingsParentNeedDown(csmc.tree, csmc.cliq, logger=csmc.logger)

  # go to 7
  return determineCliqNeedDownMsg_StateMachine
end

"""
    $SIGNATURES

Notes
- State machine function nr. 5
"""
function blockUntilSiblingsStatus_StateMachine(csmc::CliqStateMachineContainer)
  infocsm(csmc, "5, blocking on parent until all sibling cliques have valid status")
  setCliqDrawColor(csmc.cliq, "blueviolet")
  csmc.drawtree ? drawTree(csmc.tree, show=false, filepath=joinpath(getSolverParams(csmc.dfg).logpath,"bt.pdf")) : nothing

  cliqst = getCliqStatus(csmc.cliq)
  infocsm(csmc, "5, block on siblings")
  prnt = getParent(csmc.tree, csmc.cliq)
  if length(prnt) > 0
    infocsm(csmc, "5, has parent clique=$(prnt[1].index)")
    blockCliqUntilChildrenHaveUpStatus(csmc.tree, prnt[1], csmc.logger)
  end

  # go to 6c
  return blockCliqSiblingsParentChildrenNeedDown_StateMachine
end

"""
    $SIGNATURES

Notes
- State machine function nr.4
"""
function isCliqNull_StateMachine(csmc::CliqStateMachineContainer)

  prnt = getParent(csmc.tree, csmc.cliq)
  cliqst = getCliqStatus(csmc.oldcliqdata)
  infocsm(csmc, "4, isCliqNull_StateMachine, $cliqst, len(prnt)=$(length(prnt)), csmc.incremental=$(csmc.incremental)")

  if cliqst in [:marginalized;]
      # if cliqst == :marginalized || cliqst in [:downsolved; :uprecycled] && ( length(prnt) == 0 ) # ||
      # 0 < length(prnt) && getCliqStatus(prnt) in [:downsolved; :uprecycled; :marginalized;] )
      # go to 10 -- Add case for IIF issue #474
      return determineCliqIfDownSolve_StateMachine
  end

  #must happen before if :null
  stdict = blockCliqUntilChildrenHaveUpStatus(csmc.tree, csmc.cliq, csmc.logger)
  csmc.forceproceed = false

  # if clique is marginalized, then no reason to continue here
  # if no parent or parent will not update

  # for recycle computed clique values case
  if csmc.incremental && cliqst == :downsolved
    csmc.incremental = false
    # might be able to recycle the previous clique solve, go to 0b
    return checkChildrenAllUpRecycled_StateMachine
  end

  if 0 == length(prnt)
    # go to 7
    return determineCliqNeedDownMsg_StateMachine
  end

  # go to 4b
  return doesCliqNeeddownmsg_StateMachine
end

"""
    $SIGNATURES

Determine if any down messages are required.

Notes
- State machine function nr.4b
"""
function doesCliqNeeddownmsg_StateMachine(csmc::CliqStateMachineContainer)

  cliqst = getCliqStatus(csmc.cliq)
  infocsm(csmc, "4b, doesCliqNeeddownmsg_StateMachine, cliqst=$cliqst")

  if cliqst != :null
    if cliqst != :needdownmsg
      # go to 6c
      return blockCliqSiblingsParentChildrenNeedDown_StateMachine
    end
  else
    # fetch (should not block)
    stdict = blockCliqUntilChildrenHaveUpStatus(csmc.tree, csmc.cliq, csmc.logger)
    chstatus = collect(values(stdict))
    len = length(chstatus)

    # if all children needdownmsg
    if len > 0 && sum(chstatus .== :needdownmsg) == len
      # TODO maybe can happen where some children need more information?
      infocsm(csmc, "4b, escalating to :needdownmsg since all children :needdownmsg")
      notifyCliqUpInitStatus!(csmc.cliq, :needdownmsg, logger=csmc.logger)
      setCliqDrawColor(csmc.cliq, "yellowgreen")
      csmc.drawtree ? drawTree(csmc.tree, show=false, filepath=joinpath(getSolverParams(csmc.dfg).logpath,"bt.pdf")) : nothing

      # go to 5
      return blockUntilSiblingsStatus_StateMachine
    end

    # go to 6c
    return blockCliqSiblingsParentChildrenNeedDown_StateMachine
  end # != :null

  areChildDown = areCliqChildrenNeedDownMsg(csmc.tree, csmc.cliq)
  infocsm(csmc, "4b, areCliqChildrenNeedDownMsg(csmc.tree, csmc.cliq)=$(areChildDown)")
  # if cliqst == :needdownmsg
    if areChildDown
      infocsm(csmc, "4b, must deal with child :needdownmsg")
      csmc.forceproceed = true
    else
      # go to 5
      return blockUntilSiblingsStatus_StateMachine
    end
  # end

  # go to 6c
  return blockCliqSiblingsParentChildrenNeedDown_StateMachine
end

"""
    $SIGNATURES

Build a sub factor graph for clique variables from the larger factor graph.

Notes
- State machine function nr.2
"""
function buildCliqSubgraph_StateMachine(csmc::CliqStateMachineContainer)
  # build a local subgraph for inference operations
  infocsm(csmc, "2, build subgraph syms=$(getCliqAllVarIds(csmc.cliq))")
  buildCliqSubgraph(csmc.dfg, csmc.tree, csmc.cliq, csmc.cliqSubFg)

  # if dfg, store the cliqSubFg for later debugging
  dbgSaveDFG(csmc.cliqSubFg, "cliq$(csmc.cliq.index)/fg_build")

  # go to 4
  return isCliqNull_StateMachine
end

"""
    $SIGNATURES

Build a sub factor graph for clique variables from the larger factor graph.

Notes
- State machine function nr.2r
"""
function buildCliqSubgraphForDown_StateMachine(csmc::CliqStateMachineContainer)
  # build a local subgraph for inference operations
  syms = getCliqAllVarIds(csmc.cliq)
  infocsm(csmc, "2r, build subgraph syms=$(syms)")
  csmc.cliqSubFg = buildSubgraphFromLabels(csmc.dfg, syms)

  opts = getSolverParams(csmc.dfg)
  # store the cliqSubFg for later debugging
  if opts.dbg
    mkpath(joinpath(opts.logpath,"logs/cliq$(csmc.cliq.index)"))
    DFG.saveDFG(csmc.cliqSubFg, joinpath(opts.logpath,"logs/cliq$(csmc.cliq.index)/fg_build_down"))
    drawGraph(csmc.cliqSubFg, show=false, filepath=joinpath(opts.logpath,"logs/cliq$(csmc.cliq.index)/fg_build_down.pdf"))
  end

  # go to 10
  return determineCliqIfDownSolve_StateMachine
end

"""
    $SIGNATURES

Notify possible parent if clique is upsolved and exit the state machine.

Notes
- State machine function nr.1
"""
function isCliqUpSolved_StateMachine(csmc::CliqStateMachineContainer)

  infocsm(csmc, "1, isCliqUpSolved_StateMachine")
  cliqst = getCliqStatus(csmc.cliq)

  if cliqst in [:upsolved; :downsolved; :marginalized; :uprecycled]  #moved to 4 --- csmc.incremental &&
    # prep and send upward message
    prnt = getParent(csmc.tree, csmc.cliq)
    if length(prnt) > 0
      # not a root clique
      # construct init's up msg to place in parent from initialized separator variables
      msg = prepCliqInitMsgsUp(csmc.dfg, csmc.tree, csmc.cliq)
      setCliqUpInitMsgs!(prnt[1], csmc.cliq.index, msg)
      notifyCliqUpInitStatus!(csmc.cliq, cliqst, logger=csmc.logger)
    end
    #go to 10
    return determineCliqIfDownSolve_StateMachine
  end
  # go to 2
  return buildCliqSubgraph_StateMachine
end

"""
    $SIGNATURES

Final determination on whether can promote clique to `:uprecycled`.

Notes
- State machine function nr.0b
- Assume children clique status is available
- Will return to regular init-solve if new information in children -- ie not uprecycle or marginalized
"""
function checkChildrenAllUpRecycled_StateMachine(csmc::CliqStateMachineContainer)
  count = Int[]
  chldr = getChildren(csmc.tree, csmc.cliq)
  for ch in chldr
    chst = getCliqStatus(ch)
    if chst in [:uprecycled; :marginalized]
      push!(count, 1)
    end
  end
  infocsm(csmc, "0b, checkChildrenAllUpRecycled_StateMachine -- length(chldr)=$(length(chldr)), sum(count)=$(sum(count))")

  # all children can be used for uprecycled -- i.e. no children have new information
  if sum(count) == length(chldr)
    # set up msg and exit go to 1
    sdims = Dict{Symbol,Float64}()
    for varid in getCliqAllVarIds(csmc.cliq)
      sdims[varid] = 0.0
    end
    updateCliqSolvableDims!(csmc.cliq, sdims, csmc.logger)
    setCliqStatus!(csmc.cliq, :uprecycled)
    setCliqDrawColor(csmc.cliq, "orange")
    csmc.drawtree ? drawTree(csmc.tree, show=false, filepath=joinpath(getSolverParams(csmc.dfg).logpath,"bt.pdf")) : nothing
    # go to 1
    return isCliqUpSolved_StateMachine
  end

  # return to regular solve, go to 2
  return buildCliqSubgraph_StateMachine
end

"""
    $SIGNATURES

Notify possible parent if clique is upsolved and exit the state machine.

Notes
- State machine function nr.0
- can recycle if two checks:
  - previous clique was identically downsolved
  - all children are also :uprecycled
"""
function testCliqCanRecycled_StateMachine(csmc::CliqStateMachineContainer)
  # @show getCliqFrontalVarIds(csmc.oldcliqdata), getCliqStatus(csmc.oldcliqdata)
  infocsm(csmc, "0., $(csmc.incremental) ? :uprecycled => getCliqStatus(csmc.oldcliqdata)=$(getCliqStatus(csmc.oldcliqdata))")

  # check if should be trying and can recycle clique computations
  if csmc.incremental && getCliqStatus(csmc.oldcliqdata) == :downsolved
    # check if a subgraph will be needed later
    if csmc.dodownsolve
      # yes need subgraph and need more checks, so go to 2
      return buildCliqSubgraph_StateMachine
    else
       # one or two checks say yes, so go to 4
      return isCliqNull_StateMachine
    end
  end

  # nope, regular clique init-solve, go to 1
  return isCliqUpSolved_StateMachine
end

"""
    $SIGNATURES

EXPERIMENTAL: perform upward inference using a state machine solution approach.

Notes:
- will call on values from children or parent cliques
- can be called multiple times
- Assumes all cliques in tree are being solved simultaneously and in similar manner.
- State machine rev.1 -- copied from first TreeBasedInitialization.jl.
- Doesn't do partial initialized state properly yet.
"""
function cliqInitSolveUpByStateMachine!(dfg::G,
                                        tree::BayesTree,
                                        cliq::TreeClique;
                                        N::Int=100,
										oldcliqdata::BayesTreeNodeData=emptyBTNodeData(),
                                        drawtree::Bool=false,
                                        show::Bool=false,
                                        incremental::Bool=true,
                                        limititers::Int=-1,
                                        upsolve::Bool=true,
                                        downsolve::Bool=true,
                                        recordhistory::Bool=false,
                                        delay::Bool=false,
                                        logger::SimpleLogger=SimpleLogger(Base.stdout)) where {G <: AbstractDFG, AL <: AbstractLogger}
  #
  children = TreeClique[]
  for ch in Graphs.out_neighbors(cliq, tree.bt)
    push!(children, ch)
  end
  prnt = getParent(tree, cliq)

  destType = (G <: InMemoryDFGTypes) ? G : InMemDFGType#GraphsDFG{SolverParams}

  #csmc = CliqStateMachineContainer(dfg, initfg(destType), tree, cliq, prnt, children, false, incremental, drawtree, downsolve, delay, getSolverParams(dfg), oldcliqdata, logger)
  csmc = CliqStateMachineContainer(dfg, initfg(destType, params=getSolverParams(dfg)), tree, cliq, prnt, children, false, incremental, drawtree, downsolve, delay, getSolverParams(dfg), oldcliqdata, logger)

  nxt = upsolve ? testCliqCanRecycled_StateMachine : (downsolve ? testCliqCanRecycled_StateMachine : error("must attempt either up or down solve"))

  statemachine = StateMachine{CliqStateMachineContainer}(next=nxt)
  while statemachine(csmc, verbose=true, iterlimit=limititers, recordhistory=recordhistory); end
  statemachine.history
end

#