text:
### Instruction:
### This file runs post estimation after mainMulti.jl is complete.
# The CSV files exported from mainMulti.jl must be in the same output folder.
# For definition of function inputs, see UnitTest.jl and PostEstim.jl

push!(LOAD_PATH,pwd())

using Revise
using TypeDef, Helper, Estim, Gmm, GMMInterface, UserFileIO, UnitTest
using PostEstim

const ver   = "824"
const fCSV1 = "gmmStage1.csv"
const fCSV2 = "gmmStage2.csv"
const fJLD1 = "gmmStage1.jld"
const fJLD2 = "gmmStage2.jld"
const fINT1 = "gmmStage1Interim.csv"
const fINT2 = "gmmStage2Interim.csv"
const fileName = fileNameType(fCSV1,fCSV2,fJLD1,fJLD2,fINT1,fINT2)

### Export e(δ) for weak IV test
testIV!( 0, 0,20,"opt",ver,"base")
testIV!(15, 0,10,"opt",ver,"base")
testIV!( 0, 2, 1,"blp",ver,"base")
testIV!( 0, 2, 1,"diff-quad",ver,"base")
testIV!( 0, 2, 1,"diff-local",ver,"base")
testIV!(15, 2, 1,"blp",ver,"base")
testIV!(15, 2, 1,"diff-quad",ver,"base")
testIV!(15, 2, 1,"diff-local",ver,"base")

### Print Tables 4, A.1 & A.15
printEstim!( 0, 0,20,"opt",ver,"base")
printEstim!(15, 0,10,"opt",ver,"base")
printEstim!(27, 0,17,"opt",ver,"base",true,0,3000)
printEstim!(15, 0, 2,"opt",ver,"ms15")
printEstim!(15, 0,16,"opt",ver,"noAllow")
printEstim!(15, 0, 2,"opt",ver,"extra")

### Print Table A.17: Comparison of alternative IVs
printEstim!( 0, 2, 1,"blp",ver,"base",false)
printEstim!( 0, 2, 1,"diff-quad",ver,"base",false)
printEstim!( 0, 2, 1,"diff-local",ver,"base",false)
printEstim!(15, 2, 1,"blp",ver,"base",false)
printEstim!(15, 2, 1,"diff-quad",ver,"base",false)
printEstim!(15, 2, 1,"diff-local",ver,"base",false)

### Table A.2: Monte Carlo Average income for subscribers (yAvg)
yAvg, y025, y975 = MonteCarloIncome(15,10,200,"opt",ver,"base",true;lowB=false)

### Bootstrap samples export for Monte Carlo simulation
yAvg, y025, y975 = MonteCarloIncome( 0,20,200,"opt",ver,"base",true)
yAvg, y025, y975 = MonteCarloIncome(15,10,200,"opt",ver,"base",true)

max_stars_repo_name: yutec/FightBrand

max_issues_repo_name: yutec/FightBrand

max_forks_repo_name: yutec/FightBrand


text:
function cleanup()
    tmp = SwaggerMarkdown.TMP
    rm(tmp, force=true, recursive=true)
end

ext: jl


text:
struct PointerMatrix{T,P <: AbstractStridedPointer} <: AbstractMatrix{T}
    ptr::P
    size::Tuple{Int,Int}
    PointerMatrix(ptr::P, size::Tuple{Int,Int}) where {T, P <: AbstractStridedPointer{T}} = new{T,P}(ptr, size)
end
@inline PtrMatrix(A::AbstractMatrix) = PointerMatrix(stridedpointer(A), size(A))
@inline Base.pointer(A::PointerMatrix) = pointer(A.ptr)
@inline Base.size(A::PointerMatrix) = A.size
@inline Base.strides(A::PointerMatrix) = strides(A.ptr)
@inline stridedpointer(A::PointerMatrix) = A.ptr

@inline Base.maybeview(A::PointerMatrix, r::UnitRange, c::UnitRange) = PointerMatrix(gesp(A.ptr, (first(r) - 1, (first(c) - 1))), (length(r), length(c)))
# getindex is important for the sake of printing the AbstractPointerMatrix. If we call something a Matrix, it's nice to support the interface if possible.
Base.@propagate_inbounds function Base.getindex(A::PointerMatrix, i::Integer, j::Integer)
    @boundscheck begin
        M, N = size(A)
        (M < i || N < j) && throw(BoundsError(A, (i,j)))
    end
    vload(stridedpointer(A), (i-1, j-1))
end
Base.@propagate_inbounds function Base.getindex(A::PointerMatrix, v, i::Integer, j::Integer)
    @boundscheck begin
        M, N = size(A)
        (M < i || N < j) && throw(BoundsError(A, (i,j)))
    end
    vstore!(stridedpointer(A), v, (i-1, j-1))
end
Base.IndexStyle(::Type{<:PointerMatrix}) = IndexCartesian()


text:
using StochasticDiffEq, Test
u0=1/2
f(u,p,t) = u
g(u,p,t) = u
dt = 1//2^(4)
tspan = (0.0,1.0)
prob = SDEProblem(f,g,u0,(0.0,1.0))
sol = solve(prob,EM(),dt=dt)
@test_broken @inferred solve(prob,EM(),dt=dt)


text:
function createsegmentdf(segments)
  segment_types = [Symbol, Symbol, Float64, Float64, Float64, Float64]
  segment_names = [:eltype, :material, :length, :id, :od, :frictioncoefficient]
  #segs = reshape(segments, (length(segment_types), :))

  segmentdf = DataFrame(segment_types, segment_names, 0)
  for i in 1:size(segments, 1)
    seg = segments[i]
      append!(segmentdf,
        DataFrame(eltype=seg[1], material=seg[2], length=seg[3], id=seg[4], od=seg[5],
          frictioncoefficient=seg[6])
      )
  end

  segmentdf

end


text:
module Methods
import Base.push!, Base.>=, Base.&
import ..Dispatch.domainof, ..Dispatch.signatureof, ..Dispatch.make_namer
import ..Dispatch.is_empty_domain, ..Dispatch.hullof
import ..Nodes
using ..Meta, ..PartialOrder, ..Patterns, ..Dispatch, ..Encode, ..Recode

export @pattern, show_dispatch

# ---- Method -----------------------------------------------------------------

type Method
    sig::Pattern
    bindings::Vector{Node}
    body::Union{Function,Void}
    body_ex
    hullT::Union{Tuple,Void}
    id::Int

    function Method(sig::Pattern, bs, body, body_ex)
        new(sig, bs, body, body_ex, Nodes.julia_signature_of(sig))
    end
end

const nomethod = Method(Pattern(anything), Node[], nothing, nothing)

>=(x::Method, y::Method) = x.sig.intent >= y.sig.intent
# ==(x::Method, y::Method) = x.sig.intent == y.sig.intent
(&)(m::Method,  i::Intension) = m.sig.intent & i

domainof(m::Method) = m.sig.intent
hullof(m::Method)   = m.hullT === nothing ? anything : intension(m.hullT)
signatureof(m::Method)         = m.sig
signatureof(m::Method, suffix) = suffix_bindings(m.sig, suffix)

function make_namer(methods::Vector{Method})
    (node::Node)->begin
        for method in methods
            rb = method.sig.rev_bindings
            if haskey(rb, node)
                return symbol(string(rb[node], '_', method.id))
            end
        end
        nothing
    end
end

is_empty_domain(domain::Intension) = (domain === naught)

# ---- MethodTable ------------------------------------------------------------

typealias MethodNode PartialOrder.Node{Method}

type MethodTable
    name::Symbol
    top::MethodNode

    compiled::Bool
    julia_methods::Dict{Tuple,Any}
    method_counter::Int
    f::Function

    function MethodTable(name::Symbol)
        mt = new(name, MethodNode(nomethod), false, Dict{Tuple,Any}(), 0)
        mt.f = eval(:(let
                function $name(args...)
                    compile!($(quot(mt)))
                    $(quot(mt)).f(args...)
                end
                $name
            end))
        mt
    end
end

methodsof(mt::MethodTable) = [m.value for m in ordered_subDAGof(mt.top)]

function push!(mt::MethodTable, m::Method)
    m.id = (mt.method_counter += 1)
    addmethod!(mt.top, mt.name, m)

    # todo: recompile only when necessary (when is that?)
    if mt.compiled;  compile!(mt)  end
end

function compile!(mt::MethodTable)
    mt.compiled = true
    eval(:(let
            const f = $(quot(mt.f))
            f(args...) = error("No matching pattern method found")
        end))

    # NB: Lists methods in topological order;
    # avoids ambiguity warnings from julia as long
    # as there is no ambiguity among the patterns.
    methods = methodsof(mt)
#    hullTs = Tuple[m.hullT for m in filter(m->(m != nomethod), methods)]
    # Each element is a Tuple, but how to express that to the comprehension?
    hullTs = [m.hullT for m in filter(m->(m != nomethod), methods)]

    compiled = Set{Tuple}()
    for hullT in reverse(hullTs)
        if hullT in compiled; continue end
        push!(compiled, hullT)
        compile!(mt, methods, hullT)
    end
end

function compile!(mt::MethodTable, methods::Vector{Method}, hullT::Tuple)
    hull  = intension(hullT)
    top   = simplify(mt.top, hull)
    dtree = build_dtree(top)

    argsyms = seq_dispatch!(dtree, methods, hullT)
    code    = code_dispatch(dtree)

    args = Any[:($argsym::$(quot(T))) for (argsym,T) in zip(argsyms,hullT)]
    fdef = Expr(:function, :( $(mt.name)($(args...)) ), code)
    mt.julia_methods[hullT] = Expr(:function, :( dispatch($(args...)) ), code)

    eval(:(let
            const f = $(quot(mt.f))
            f($(args...)) = $code
        end))
end

show_dispatch(mt::MethodTable, args...) = show_dispatch(STDOUT,
                                                        mt, args...)
show_dispatch(io::IO, mt::MethodTable) = show_dispatch(io, mt, nothing)
function show_dispatch(io::IO, mt::MethodTable, Ts::Union{Tuple,Void})
    if !mt.compiled;  compile!(mt)  end

    println("const ", mt.name, " = (args...)->dispatch(args...)")

    println("\n# ---- Pattern methods: ----")
    methods = methodsof(mt)
    mnames = Dict{Function,Symbol}()
    for (id, method) in sort([(m.id, m) for m in methods])
        if method === nomethod; continue end

        println(io, "# ", mt.name, method.sig)
        mname = symbol(string("match", method.id))
        mnames[method.body] = mname

        args = keys(method.sig.bindings) # Right order? Does it matter?
        Base.show_unquoted(io, Expr(:function, :($mname($(args...))),
                                    method.body_ex))
        print(io,"\n\n")
    end

    subs_invocation(ex) = begin
        if is_expr(ex, :quote) && haskey(mnames, ex.args[1]); mnames[ex.args[1]]
        else; nothing
        end
    end

    println("# ---- Dispatch methods: ----")
    for (f_Ts, fdef) in mt.julia_methods
        if (Ts === nothing) || (Tuple{f_Ts...} <: Tuple{Ts...})
            Base.show_unquoted(io, subs_ex(subs_invocation, fdef))
            print(io, "\n\n")
        end
    end
end

# ---- @pattern etc. ----------------------------------------------------------

const method_tables = Dict{Function, MethodTable}()

function show_dispatch(f::Function, args...)
    if !haskey(method_tables, f);  error("not a pattern function: $f")  end
    show_dispatch(method_tables[f], args...)
end

macro pattern(ex)
    code_pattern(ex)
end
function code_pattern(ex)
    sig, body = try
        split_fdef(ex)
    catch e
        error("@pattern: not a function definition ($e)")
    end
    @expect is_expr(sig, :call)
    fname, args = sig.args[1], sig.args[2:end]
    if is_expr(fname, :curly)
        error("@pattern: type parameters are not implemented")
    end
    psig = :($(args...),)
    p_ex, bodyargs = recode(psig)

    f = esc(fname::Symbol)
    quote
        local p = $p_ex
        local bindings = Node[p.bindings[name] for name in $(quot(bodyargs))]
        local bodyfun = $(esc(:(($(bodyargs...),)->$body)))
        local method = Method(p, bindings, bodyfun, $(quot(body)))

        local wasbound = try
            f = $f
            true
        catch e
            false
        end

        local mt
        if !wasbound
            mt = MethodTable($(quot(fname)))
            local const f = mt.f
            const $f = (args...)->f(args...)
            method_tables[$f] = mt
            push!(mt, method)
        else
            if !haskey(method_tables, $f)
                error($("$fname is not a pattern function"))
            end
            push!(method_tables[$f], method)
        end
    end
end

end # module


text:
using JavaCall
DBtableTypeorString = Union{DBtableType, AbstractString}

function ktrussadj(A::DBtableTypeorString, Rname::AbstractString, k::Number; filterRowCol::AbstractString="", forceDelete::Bool=true, maxiter = 2^31 - 1::Number)
# /**
# * From input unweighted, undirected adjacency table Aorig, put the k-Truss
# * of Aorig in Rfinal.
# * @param Aorig Unweighted, undirected adjacency table.
# * @param Rfinal Does not have to previously exist. Writes the kTruss into Rfinal if it already exists.
# *               Use a combiner if you want to sum it in.
# * @param k Trivial if k <= 2.
# * @param filterRowCol Filter applied to rows and columns of Aorig
# *                     (must apply to both rows and cols because A is undirected Adjacency table).
# * @param forceDelete False means throws exception if the temporary tables used inside the algorithm already exist.
# *                    True means delete them if they exist.
# * @param Aauthorizations Authorizations for scanning Atable. Null means use default: Authorizations.EMPTY
# * @param RNewVisibility Visibility label for new entries created. Null means no visibility label.
# * @param maxiter A bound on the number of iterations. The algorithm will halt
# *                either at convergence or after reaching the maximum number of iterations.
# *                Notemitimmy that if the algorithm stops before convergence, the result may not be correct.
# * @return nnz of the kTruss subgraph, which is 2* the number of edges in the kTruss subgraph.
# *          Returns -1 if k < 2 since there is no point in counting the number of edges.
# */
# public long kTrussAdj(String Aorig, String Rfinal, int k,
#                      String filterRowCol, boolean forceDelete,
#                      Authorizations Aauthorizations, String RNewVisibility,
#                      int maxiter) {

    if isa(A, DBtablePair)
        Aname = A.name1
    elseif isa(A, DBtable)
        Aname = A.name
    else
        Aname = A
    end

    JAuthorizations = @jimport "org.apache.accumulo.core.security.Authorizations"
    emptyAuth = JAuthorizations((), )
    newVisibility = ""

    return jcall(DB.Graphulo,"kTrussAdj", jlong,
        (JString, JString, jint, JString, jboolean, JAuthorizations, JString, jint),
        Aname, Rname, k, toDBstring(filterRowCol), forceDelete, emptyAuth, newVisibility, maxiter)

end

function ktrussedge(E::DBtableTypeorString, ET::DBtableTypeorString, Rname::AbstractString, RTname::AbstractString, k::Number; edgeFilter::AbstractString="", forceDelete::Bool=true)
    # /**
    # * From input unweighted, undirected incidence table Eorig, put the k-Truss
    # * of Eorig in Rfinal.  Needs transpose ETorig, and can output transpose of k-Truss subgraph too.
    # * @param Eorig Unweighted, undirected incidence table.
    # * @param ETorig Transpose of input incidence table.
    # *               Optional.  Created on the fly if not given.
    # * @param Rfinal Does not have to previously exist. Writes the kTruss into Rfinal if it already exists.
    # *               Use a combiner if you want to sum it in.
    # * @param RTfinal Does not have to previously exist. Writes in the transpose of the kTruss subgraph.
    # * @param k Trivial if k <= 2.
    # * @param edgeFilter Filter on rows of Eorig, i.e., the edges in an incidence table.
    # * @param forceDelete False means throws exception if the temporary tables used inside the algorithm already exist.
    # *                    True means delete them if they exist.
    # * @param Eauthorizations Authorizations for scanning Atable. Null means use default: Authorizations.EMPTY
    # * @return  nnz of the kTruss subgraph, which is 2* the number of edges in the kTruss subgraph.
    # *          Returns -1 if k < 2 since there is no point in counting the number of edges.
    # *  public long kTrussEdge(String Eorig, String ETorig, String Rfinal, String RTfinal, int k,
    # String edgeFilter, boolean forceDelete, Authorizations Eauthorizations)

    if isa(E, DBtablePair)
        Ename = E.name1
    elseif isa(E, DBtable)
        Ename = E.name
    else
        Ename = E
    end

    if isa(ET, DBtablePair)
        ETname = ET.name1
    elseif isa(E, DBtable)
        ETname = ET.name
    else
        ETname = ET
    end

    JAuthorizations = @jimport "org.apache.accumulo.core.security.Authorizations"
    emptyAuth = JAuthorizations((), )
    newVisibility = ""

    return jcall(DB.Graphulo,"kTrussEdge", jlong,
        (JString, JString, JString, JString, jint, JString, jboolean, JAuthorizations, JString, jint),
        Ename, ETname, Rname, RTname, k, toDBstring(edgeFilter), forceDelete, emptyAuth)

end


text:
module simulate_quantum_memory
  using Symplectic
  using css_code
  using noise_models
  export SingleAngleModel
  using Distributions: Geometric, Categorical, Bernoulli, Normal

  # Utility functions
  function changelocs(arr)
    n = size(arr)[1];
    locs = collect(1:n-1)[(arr[2:n]-arr[1:n-1] .!= 0)];
    return locs
  end
  function nonzerolocs(arr)
    n = size(arr)[1];
    locs = collect(1:n)[(arr .!= 0)];
    return locs
  end
  function sum_to_right_edge(arr, edges)
    sums = zeros(length(edges));
    edges = [1; edges];
    for k=1:length(sums)
        sums[k] = sum(arr[edges[k]:edges[k+1]]);
    end
    return sums
  end
  export changelocs, nonzerolocs, sum_to_right_edge

################################################################################
##                     [[15, 1, 3]] Asymmetric CSS Code                       ##
################################################################################

    function RMCode()
        # The bitstrings that describe each stabilizer operator
        x_S=[[0 0 0 0 0 0 0 1 1 1 1 1 1 1 1];
             [0 0 0 1 1 1 1 0 0 0 0 1 1 1 1];
             [0 1 1 0 0 1 1 0 0 1 1 0 0 1 1];
             [1 0 1 0 1 0 1 0 1 0 1 0 1 0 1]];
        z_S=[[0 0 0 0 0 0 0 1 1 1 1 1 1 1 1];
             [0 0 0 1 1 1 1 0 0 0 0 1 1 1 1];
             [0 1 1 0 0 1 1 0 0 1 1 0 0 1 1];
             [1 0 1 0 1 0 1 0 1 0 1 0 1 0 1];
             [1 0 0 0 0 0 1 0 0 1 0 1 0 0 0];
             [0 1 0 0 0 0 1 0 1 0 0 1 0 0 0];
             [0 0 1 0 0 0 1 1 0 0 0 1 0 0 0];
             [0 0 0 1 0 0 1 0 1 1 0 0 0 0 0];
             [0 0 0 0 1 0 1 1 0 1 0 0 0 0 0];
             [0 0 0 0 0 1 1 1 1 0 0 0 0 0 0]];

        # The same operators, now in symplectic form
        X_S = Symp(zeros(size(x_S)), (x_S .== 1));
        Z_S = Symp((z_S .== 1), zeros(size(z_S)));
        code = Code(1, 3, X_S, Z_S)
        return code
    end
    export RMCode

################################################################################
##                 Single-angle dephasing noise simulation                    ##
################################################################################

    function single_angle_full_simulation(model::SingleAngleModel,
                                          code::Code,
                                          MAX_T::Number;
                                          simulate_decoder=simulate_decoder::Bool,
                                          report=report::Bool)

     #########################
     #         Setup         #
     #########################

     theta_true = (0.05+0.9*rand())*pi/2;
     theta_est = pi/4;
     time = 0;
     if report == true
       W_X = [0.0];
       W_Z = [0.0];
       T0s = [0.0];
       TEs = [0.0];
       TIMEs = [0.0];
       report_idx = 1;
     end

     #########################
     #       Simulation      #
     #########################

     for T=1:MAX_T
         original_error = rand_error(model, code.n, theta_true - theta_est);
         w_x = sum(original_error.x); w_z = sum(original_error.z);

         if simulate_decoder == true
           (s_x, s_z) = code.measure(original_error);
           decoded_error, flag = code.decode(s_x, s_z);
         else
           decoded_error, flag = original_error, convert(Int64, (w_x <= 1) & (w_z <=3))
         end

         if flag != 1
             break
         else
             time = time + 1
         end

         model = update(model, decoded_error, theta_est);
         theta_est = deepcopy(model.theta[indmax(model.p_theta)]);

         ## REPORTING
         if report == true
           push!(W_X, w_x);
           push!(W_Z, w_z);
           push!(T0s, theta_true);
           push!(TEs, theta_est);
           push!(TIMEs, time);
           report_idx = report_idx + 1;
         end
     end

     #########################
     #         Report        #
     #########################

     if report == false
       return (time, abs(theta_true - theta_est))
     else
       return (TIMEs, T0s, TEs, W_X, W_Z)
     end

    end
    export single_angle_full_simulation

################################################################################
##           Single-angle dephasing noise fast-forward simulation             ##
################################################################################

    function single_angle_fast_forward(model::SingleAngleModel,
                                       MAXTIME::Number;
                                       report=report::Bool)

        #########################
        #         Setup         #
        #########################
        time = 0;
        EXPORT = Dict("time" => 0.0, "exit_status" => 0.0, "p_uncorrectable" => 0.0)
        if report == true
            REPORT = Dict("W_X" => [0.0],
                          "W_Z" => [0.0],
                          "T0s" => [0.0],
                          "TEs" => [0.0],
                          "TIMEs" => [0.0],
                          "CYCLE_TYPEs" => [0],
                          "PROFILEs" => [0.0],
                          "M_HATs" => [0],
                          "M_TRUEs" => [0]);
        end

        #########################
        #       Simulation      #
        #########################
        uninformative_continue = true;
        informative_continue = true;
        exit_status = 0;
        while (time < MAXTIME) & uninformative_continue & informative_continue

            idx = convert(Int64, 2*model.grains + model.M_hat);

            # Keep simulating plateaus

            ## SAMPLE TIME
            # Sample the amount time until the next informative error
            p_inf = model.sampling_rates["W_Z>0"][idx]
            T_inf = Geometric(p_inf);
            t_inf = rand(T_inf);

            # Sample uncorrectable uninformative errors that occured during t
            f_uninf = -randexp() / log1p(-model.sampling_rates["W_X>d_x"][idx]);

            # Test for failure
            uninformative_continue = (t_inf < f_uninf);

            ## SAMPLE INFORMATIVE ERROR
            # Sample the weight of the informative error
            M_inf = Bernoulli(model.sampling_rates["W_Z>d_z|W_Z>0"][idx]);
            m_inf = rand(M_inf)+1;

            # Test for failure
            informative_continue = (m_inf <= model.d_z);

            ## SAMPLE UNINFORMATIVE ERROR
            # Sample the number of uninformative errors that have occured up to t
            n_uninf = 0;
            p_uninf = model.sampling_rates["W_X>0"][idx];
            p_cat = model.sampling_rates["W_X=k|d_x>=W_X>0"][:, idx]
            if t_inf * p_uninf < 5
                tic();
                cycle_type=0;
                t_skip = 0;
                Et_uninf = 1/p_uninf;
                T_uninf = Geometric(p_uninf);
                while (t_skip <= t_inf)
                    t_uninf = rand(T_uninf);
                    n_uninf = n_uninf + 1;
                    t_skip = t_skip + t_uninf;
                end
                M_uninf = Categorical(p_cat);
                m_uninf = rand(M_uninf, n_uninf);
                m_uninfs = hist(m_uninf,0.5:model.d_x+0.5)[2];
                profile = toq();
            else
                tic();
                N_uninf = Normal(t_inf*p_uninf, t_inf*p_uninf*(1-p_uninf));
                n_uninf = -1;
                while n_uninf < 0
                    n_uninf = convert(Int64, floor(rand(N_uninf)));
                end
                if n_uninf > 1 / p_cat[2]
                    cycle_type=2;
                    m_uninfs = convert(Array{Int64, 1}, floor(n_uninf .* p_cat));
                else
                    cycle_type=1;
                    M_uninf = Categorical(p_cat);
                    m_uninf = rand(M_uninf, n_uninf);
                    m_uninfs = hist(m_uninf,0.5:model.d_x+0.5)[2];
                end
                profile = toq();
            end

            ## UPDATE MODEL
            # Informative update
            time = copy(time + t_inf)

            # Update the probability distribution based on the informative error
            idx = convert(Int64, 2*model.grains-model.M_hat);
            g = convert(Int64, model.grains-1);
            update = model.update_functions["Z=1, X=0"][idx:idx+g]';

            # Uninformative update
            for k=1:model.d_x
                update = update + m_uninfs[k] * model.update_functions["Z=0, X=1"][idx:idx+g]'
            end

            model.p_theta = model.p_theta + update[:];
            model.p_theta = model.p_theta - mean(model.p_theta);

            # Update the MLE
            model.M_hat = indmax(model.p_theta);
            idx_hat = convert(Int64, 2*model.grains + model.M_hat);
            model.theta_hat = copy(model.M[idx_hat]);

            M_d = model.M_hat - model.M_true;
            p_uncorr = model.sampling_rates["uncorrectable"][idx_hat]

            ## REPORTING
            if report == true
                print("$M_d, $p_inf, $p_uncorr, $t_inf\n");
                push!(REPORT["W_X"], sum(m_uninfs));
                push!(REPORT["W_Z"], m_inf);
                push!(REPORT["T0s"], model.theta_true);
                push!(REPORT["TEs"], model.theta_hat);
                push!(REPORT["TIMEs"], time);
                push!(REPORT["CYCLE_TYPEs"], cycle_type);
                push!(REPORT["PROFILEs"], profile);
                push!(REPORT["M_HATs"], model.M_hat);
                push!(REPORT["M_TRUEs"], model.M_true);
            end

            # EXITING THE LOOP
            if !informative_continue
                exit_status = 0 # informative_overweight
            elseif !uninformative_continue
                exit_status = 1 # uninformative_overweight
            elseif !(time < MAXTIME)
                exit_status = 2 # time_out
            end
        end

        #########################
        #         Report        #
        #########################

        EXPORT["time"] = time;
        EXPORT["exit_status"] = exit_status;
        idx_hat = convert(Int64, 2*model.grains + model.M_hat);
        EXPORT["p_uncorrectable"] = model.sampling_rates["uncorrectable"][idx_hat];

        if report == false
            return EXPORT
        else
            return (EXPORT, REPORT)
        end
  end
  export single_angle_fast_forward
end


text:
# ------------------------------------
# Tests adapted from Git2Go Library
# ------------------------------------
context("test Signature constructor") do
    test_name   = "test test"
    test_email  = "test@test.com"
    test_time   = 0
    test_offset = 0

    # throw on bad name
    @test_throws LibGitError{:Invalid,:Error} Signature("", "")

    let sig = Signature(test_name, test_email, test_time, test_offset)
        @test sig.time == time(sig)
        @test sig.time_offset == time_offset(sig)
        @test time(sig) == test_time
        @test time_offset(sig) == test_offset
    end

    let now = floor(Int, time()), sig = Signature(test_name, test_email)
        @test time(sig) == now
    end
end

# TODO: Win32 Failure
context("test round-trip") do
    sig1 = Signature("test1", "test1@test.com")
    sig1_ptr = convert(Ptr{LibGit2.SignatureStruct}, sig1)
    sig2 = Signature(sig1_ptr)
    @test sig1 == sig2
end

# -----------------------------------------
# Tests adapted from Ruby's Rugged Library
# -----------------------------------------


text:
# unit tests of extract_string
x = "12345678"
@assert DataFrames.extract_string(x, 3, 6) == "3456"
@assert DataFrames.extract_string(x, 3, 3) == "3"
@assert DataFrames.extract_string(x, 3, 6, Set(3)) == "456"
@assert DataFrames.extract_string(x, 3, 6, Set(5, 3)) == "46"
x = "\"Güerín\",\"Sí\",\"No\""
@assert DataFrames.extract_string(x, 2, 7, Set(3)) == "Gerí"
@assert DataFrames.extract_string("", 0,0,Set()) == ""

# Handle the empty string properly
test0 = IOString("")
res0 = DataFrames.read_separated_line(test0, ',', '"')
@assert isempty(res0)

test1 = IOString("I'm A,I'm B,I'm C,-0.3932755625236671,20.157657978753534")
res1 = DataFrames.read_separated_line(test1, ',', '"')
@assert res1[2] == "I'm B"
@assert res1[4] == "-0.3932755625236671"

test2 = IOString("123, 456 , \"789\",TRUE")
res2 = DataFrames.read_separated_line(test2, ',', '"')
@assert res2[2] == "456"
@assert res2[3] == "789"

test3 = IOString("123 ,456 , \"789\" , \"TRUE\"")
res3 = DataFrames.read_separated_line(test3, ',', '"')
@assert res3[2] == "456"
@assert res3[4] == "TRUE"

test4 = IOString("123 ,456 ,  , \"TRUE\"")
res4 = DataFrames.read_separated_line(test4, ',', '"')
@assert res4[3] == ""
@assert res4[4] == "TRUE"

test5 = IOString("123 ,456 , \"a\"\"b\" ,\"TRUE\"")
res5 = DataFrames.read_separated_line(test5, ',', '"')
@assert res5[3] == "a\"b"
@assert res5[4] == "TRUE"

test6 = IOString("123 ,456 , \"a
b\" ,\"TRUE\"")
res6 = DataFrames.read_separated_line(test6, ',', '"')
@assert res6[3] == "a\nb"
@assert res6[4] == "TRUE"

# Should this be one NA?
# test7 = IOString("")
# res7 = DataFrames.read_separated_line(test7, ',', '"')
# @assert length(res7) == 1

test8 = IOString("a,\"b\",\"cd\",1.0,1\na,\"b\",\"cd\",1.0,1")
res8 = DataFrames.read_separated_line(test8, ',', '"')
@assert length(res8) == 5
@assert res8[5] == "1"

test9 = IOString("\"Güerín\",\"Sí\",\"No\"")
res9 = DataFrames.read_separated_line(test9, ',', '"')
@assert res9[2] == "Sí"

test10 = IOString("1,2,3,,")
res10 = DataFrames.read_separated_line(test10, ',', '"')
@assert length(res10) == 5

filename = Pkg.dir("DataFrames", "test", "data", "simple_data.csv")
open(filename,"r") do io
    t1 = read_table(io, ',', '"', DataFrames.DEFAULT_MISSINGNESS_INDICATORS, false, ["1","2","3","4","5"], 2)
    @assert nrow(t1) == 2
    @assert t1[1,2] == "b"
end

# Test separated line splitting
#
# TODO: Test minimially-quoted
# TODO: Test only-strings-quoted

separators = [',', '\t', ' ']
quotation_characters = ['\'', '"']

# Test all-entries-quoted for all quote characters and separators
items = {"a", "b", "c,d", "1.0", "1"}
item_buffer = Array(UTF8String, length(items))

# TODO: make this work with new splitting code
# for separator in separators
#   for quotation_character in quotation_characters
#     line = join(map(x -> string(quotation_character, x, quotation_character),
#                     items),
#                 separator)
#     current_item_buffer = Array(Char, strlen(line))
#     split_results = DataFrames.split_separated_line(line, separator, quotation_character, item_buffer, current_item_buffer)
#     @assert all(split_results .== items)
#   end
# end

# Test reading
@assert DataFrames.determine_separator("blah.csv") == ','
@assert DataFrames.determine_separator("blah.tsv") == '\t'
@assert DataFrames.determine_separator("blah.wsv") == ' '
# @assert DataFrames.determine_separator("blah.txt")
# Need to change to use @expects to test that error gets raised

filename = Pkg.dir("DataFrames", "test", "data", "big_data.csv")
separator = DataFrames.determine_separator(filename)
quotation_character = '"'
missingness_indicators = ["", "NA"]
header = true
column_names = UTF8String["A", "B", "C", "D", "E"]
minibatch_size = 10

file = open(filename, "r")
readline(file)
minibatch = read_minibatch(file,
	                       separator,
	                       quotation_character,
	                       missingness_indicators,
	                       column_names,
	                       minibatch_size)
@assert nrow(minibatch) == minibatch_size
@assert ncol(minibatch) == length(column_names)
@assert colnames(minibatch) == column_names
@assert eltype(minibatch[:, 1]) == UTF8String
@assert eltype(minibatch[:, 2]) == UTF8String
@assert eltype(minibatch[:, 3]) == UTF8String
@assert eltype(minibatch[:, 4]) == Float64
@assert eltype(minibatch[:, 5]) == Float64
close(file)

@elapsed df = read_table(filename)
@assert nrow(df) == 10_000
@assert ncol(df) == 5
@assert colnames(df) == column_names
@assert typeof(df[:, 1]) == DataVector{UTF8String}
@assert typeof(df[:, 2]) == DataVector{UTF8String}
@assert typeof(df[:, 3]) == DataVector{UTF8String}
@assert typeof(df[:, 4]) == DataVector{Float64}
@assert typeof(df[:, 5]) == DataVector{Float64}

# TODO: Split apart methods that perform seek() from those that don't
text_data = convert(Array{UTF8String, 2}, (["1" "3" "A"; "2" "3" "NA"; "3" "3.1" "C"]))

true_df = DataFrame(quote
                      x1 = DataArray([1, 2, 3])
                      x2 = DataArray([3, 3, 3.1])
                      x3 = DataArray(UTF8String["A", "", "C"], [false, true, false])
                    end)
df = DataFrames.convert_to_dataframe(text_data,
                                     ["", "NA"],
                                     ["x1", "x2", "x3"])
@assert isequal(df, true_df)
@assert isequal(eltype(df["x1"]), Int64)
@assert isequal(eltype(df["x2"]), Float64)
@assert isequal(eltype(df["x3"]), UTF8String)

filename = Pkg.dir("DataFrames", "test", "data", "big_data.csv")
separator = DataFrames.determine_separator(filename)
quotation_character = '"'
missingness_indicators = ["", "NA"]
header = true

nrows = DataFrames.determine_nrows(filename, header)
@assert nrows == 10_000

io = open(filename, "r")

column_names = DataFrames.determine_column_names(io, separator, quotation_character, header)
@assert column_names == UTF8String["A", "B", "C", "D", "E"]

seek(io, 0)
if header
  readline(io)
end
text_data = DataFrames.read_separated_text(io, nrows, separator, quotation_character)
@assert eltype(text_data) == UTF8String
@assert size(text_data) == (10_000, 5)

df = read_table(io,
                separator,
                quotation_character,
                missingness_indicators,
                header,
                column_names,
                nrows)
@assert nrow(df) == 10_000
@assert ncol(df) == 5
@assert colnames(df) == column_names
@assert eltype(df[:, 1]) == UTF8String
@assert eltype(df[:, 2]) == UTF8String
@assert eltype(df[:, 3]) == UTF8String
@assert eltype(df[:, 4]) == Float64
@assert eltype(df[:, 5]) == Float64

df = read_table(filename)
@assert nrow(df) == 10_000
@assert ncol(df) == 5
@assert colnames(df) == column_names
@assert eltype(df[:, 1]) == UTF8String
@assert eltype(df[:, 2]) == UTF8String
@assert eltype(df[:, 3]) == UTF8String
@assert eltype(df[:, 4]) == Float64
@assert eltype(df[:, 5]) == Float64

# TODO: Add test case in which data file has header, but no rows
# Example "RDatasets/data/Zelig/sna.ex.csv"
# "","Var1","Var2","Var3","Var4","Var5"

# Additional data sets

@elapsed df = read_table("test/data/big_data.csv")
# TODO: Make this faster
@elapsed df = read_table("test/data/movies.csv")
# TODO: Release this data set publicly
#@elapsed df = read_table("test/data/bigrams.tsv")
@elapsed df = read_table("test/data/utf8.csv")
@elapsed df = read_table("test/data/bool.csv")
@elapsed df = read_table("test/data/types.csv")
@elapsed df = read_table("test/data/space_after_delimiter.csv")
@elapsed df = read_table("test/data/space_before_delimiter.csv")
@elapsed df = read_table("test/data/space_around_delimiter.csv")
@elapsed df = read_table("test/data/corrupt_utf8.csv")


text:
using Luxor, Test

function test_polytriangulate(fname)
    width, height = 800, 800
    Drawing(width, height, fname)
    origin()
    background("white")

    rawpts = ngon(O, 300, 8, vertices = true)

    push!(rawpts, midpoint(rawpts[1], rawpts[4]))
    push!(rawpts, midpoint(rawpts[2], rawpts[5]))
    push!(rawpts, midpoint(rawpts[3], rawpts[6]))

    ptriangles = polytriangulate(rawpts)

    @test length(ptriangles) == 12

    fontsize(16)
    # visual inspection... :(
    setline(3)
    for (n, p) in enumerate(ptriangles)
        sethue([Luxor.julia_purple, Luxor.julia_blue,  Luxor.julia_red,  Luxor.julia_green][mod1(n, end)])
        pgon = Point[p[1], p[2], p[3]]
        poly(pgon,  :fillpreserve, close = true)
        sethue("white")
        strokepath()
        text(string(n), polycentroid(pgon), halign=:middle)
    end
    sethue("white")
    circle.(rawpts, 5, :fill)
    @test finish() == true
end

fname = "polytriangulate.png"
test_polytriangulate(fname)
println("...finished test: output in $(fname)")


text:
export covariance, cov2cor!, filter2d, smooth, addConstColumn, weights_masker

##################################################
# Calculates covariance matrix for binary matrix
##################################################
function covariance(bigData::BitArray{2}; chunkSize=100000, quiet=false)
    P,N = size(bigData)
    XtX = zeros(Float64, P, P)
    varSums = zeros(Float64, P, 1)

    # force the chunk size to line up with 64 bit word boundaries,
    # this is most important for loading from a file, but we also use it here.
    # we try and keep the size close to what was requested
    chunkSize = max(round(Int64, chunkSize/64),1)*64
    # build XtX incrementally and also the totals of every variable.
    chunk = Array(Float32, P, chunkSize)
    numChunks = round(Int64, ceil(N/chunkSize))
    for i in 1:numChunks-1
        chunk[:,:] = bigData[:,(i-1)*chunkSize+1:i*chunkSize]
        XtX .+= A_mul_Bt(chunk,chunk) # using a float array is important to get LAPACK speed
        varSums .+= sum(chunk,2)
        #if !quiet println(STDERR, "processed $(i*chunkSize*1000) bp...") end
    end

    # get the last unevenly sized chunk
    chunk = Array(Float32, P, N - (numChunks-1)*chunkSize)
    chunk[:,:] = bigData[:,(numChunks-1)*chunkSize+1:end]
    XtX .+= A_mul_Bt(chunk,chunk)
    varSums .+= sum(chunk,2)

    # convert XtX to a covariance matrix
    XtX .-= varSums*varSums'/N
    XtX ./= (N-1)
end

#####################################
# Converts covariance to correlation
#####################################
function cov2cor!(M)
    for i in 1:size(M)[1]
        val = sqrt(M[i,i])
        if val > 0.0
            M[i,:] ./= val
            M[:,i] ./= val
        end
    end
    M
end

######################################
# Takes in matrix and true/fase masks
######################################
function filter2d(matrix, ymask, xmask)
    temp = filter(x->x[2], [(matrix[i,1:end], ymask[i]) for i in 1:size(matrix)[1]])
    temp = [i[1] for i in temp]
    temp = hcat(temp...)

    temp = filter(x->x[2], [(temp[i,1:end], xmask[i]) for i in 1:size(temp)[1]])
    temp = [i[1] for i in temp]
    temp = hcat(temp...)
end

######################################
# Smoothing function for peak calling
######################################
function smooth(a, width)
    if width > length(a)
        width = length(a)
    end

    ret = copy(a)
    counts = ones(length(a))

    #Aggregating the values
    for i in 1:width
        for j in 1:length(a)-i
            ret[j] += a[j+i]
            ret[end-j+1] += a[end-j-i+1]
            counts[j] += 1
            counts[end-j+1] += 1
        end
    end
    ret./counts
end

######################################
# Adds constant column to data matrix
######################################
function addConstColumn(M::Array{Int64, 2})
    _const = ones(Int64, (size(M)[1], 1))
    hcat(M, _const)
end

##########################################################################
# masks weight based on some conditions
# mask = weights_masker("ENCFF000YRS", 0, JLD.load(ctrllistdata)["ctrls"])
##########################################################################
function weights_masker(expID::String, mode::Int, clist; metadata::String="../data/metadata.csv", constant::Bool=true)
    m = CSV.read(metadata)

    # mode 0 removes matched control
    if mode == 0
        cexp = convert(String, m[m[:ID].==expID, :CTRL1][1])
        ignorelist = m[m[:EXP].==cexp, :ID]

    # mode 1 removes conrol form the same cellline
    elseif mode == 1
        ct = convert(String, m[m[:ID].==expID, :CELLTYPE][1])
        ignorelist = m[(m[:IFCTRL].==true).&(m[:CELLTYPE].==ct), :ID]

    # mode 2 removes conrol form the same lab
    elseif mode == 2
        lab = convert(String, m[m[:ID].==expID, :LAB][1])
        ignorelist = m[(m[:IFCTRL].==true).&(m[:LAB].==lab), :ID]

    # combination of mode 1 and 2
    elseif mode == 3
        ct = convert(String, m[m[:ID].==expID, :CELLTYPE][1])
        ignorelist1 = m[(m[:IFCTRL].==true).&(m[:CELLTYPE].==ct), :ID]
        lab = convert(String, m[m[:ID].==expID, :LAB][1])
        ignorelist2 = m[(m[:IFCTRL].==true).&(m[:LAB].==lab), :ID]
        ignorelist = append!(ignorelist1,ignorelist2)
    end

    ignorelist = collect(Set(ignorelist))

    mask = []
    for c in clist
        if c in ignorelist
            push!(mask, false)
        else
            push!(mask, true)
        end
    end

    # some assertion
    assert(length(ignorelist) >= length(mask)-sum(mask))

    if constant
        push!(mask, true)
    end

    convert(BitArray{1}, mask)
end


text:
"""
    periodic_mean(df, target, pdir)

Filter DataFrame df by time directive and get averages

time directive must be in :hour (24 hour), :day (365 day), :quarter (4 quarter)
"""
function periodic_mean(df::DataFrame, target::Symbol, pdir::Symbol)
    if pdir == :hour
        max_period = 24
    elseif pdir == :week
        max_period = 53
    elseif pdir == :month
        max_period = 12
    else
        error("Period directive must be between :hour, :day, :month, :quarter")
    end

    avgs = zeros(max_period)

    for period=1:max_period
        if pdir == :hour
            _df = @from i in df begin
                @where hour(i.date) == period - 1
                @orderby i.date
                @select i
                @collect DataFrame
            end
        elseif pdir == :week
            _df = @from i in df begin
                @where week(i.date) == period
                @orderby i.date
                @select i
                @collect DataFrame
            end
        elseif pdir == :month
            _df = @from i in df begin
                @where month(i.date) == period
                @orderby i.date
                @select i
                @collect DataFrame
            end
        end

        avgs[period] = mean(_df[!, target])
    end

    avgs
end

"""
    periodic_mean(df, target, pdir)

decompose double seasonality

time directive must be in :hour (24 hour), :day (365 day), :quarter (4 quarter)
"""
function populate_periodic_mean(df::DataFrame, target::Symbol, pdir::Symbol)
    means = periodic_mean(df, target, pdir)

    # adjust index only for hourly (0-23)
    if pdir == :hour
        means_repeated = means[df[!, pdir] .+ 1]
    else
        means_repeated = means[df[!, pdir]]
    end

    means, means_repeated
end

"""
    populate_periodic_mean!(df1, df2, key, target)

populate `df2` column `target` to `df1` by `key` and return joined result

* df1 must have `key` column
* df2 must have `key` and `target` column
"""
#populate_periodic_mean(df1::DF, df2::DF, key::Symbol, target::Symbol) where {DF<:DataFrame} =
#    DataFrames.join(df1[!, [key]], df2[!, [key, target]], on = key)
function populate_periodic_mean(df1::DF, df2::DF, key::Symbol, target::Symbol) where {DF<:DataFrame}
    #@show names(df1[!, [key]]), names(df2[!, [key, target]])
    _joined = DataFrames.join(df1[!, [key]], df2[!, [key, target]], on = key)
    #@show first(_joined, 5)
    DataFrames.join(df1[!, [key]], df2[!, [key, target]], on = key)
end

function smoothing_series(xs, ys; span = 0.05)
    model = Loess.loess(xs, ys; span = span)

    us = collect(range(extrema(xs)...; step = 1))
    vs = Loess.predict(model, us)

    us, vs
end

function compute_annual_mean(df::DataFrame, target::Symbol)
    # to deal with leap year, just make large array for indexing
    # when moving window have some error but it is ignorable
    year_wd_mean = zeros(12*31*24)
    year_wd_count = zeros(Int64, 12*31*24)

    last_date = df[end, :date]
    window = 24

    # mean of window
    mean_vals = map(i -> mean(@view df[i:i+(window-1), ycol]), 1:size(df, 1)-window+1)
    for (i, row) in enumerate(eachrow(df))
        t = getproperty(row, :date)
        if t > last_date - Dates.Hour(window)
            break
        end
        v = mean_vals[i]

        idx = (Dates.month(t) - 1) * 31 * 24 + (Dates.day(t) - 1) * 24 + Dates.hour(t) + 1
        year_wd_count[idx] += 1
        year_wd_mean[idx] = year_wd_mean[idx] + (v - year_wd_mean[idx]) / year_wd_count[idx]
    end

    year_wd_mean
end

"""
    decompose_seasonality(df, target)

decompose seasonality to daily + annual seasonality
"""
function decompose_seasonality(df::DataFrame, target::Symbol)
    # daily mean
    # remove time information
    #df[!, :date_day] = Date.(DateTime.(df[!, :date], Local))
    df[!, :hour] = Dates.hour.(df[!, :date])

    # find out periodic mean (averaged over hours)
    _df_avg_hour = DataFrames.by(df, :hour, avg_hour = target => mean)
    # same as df_avg_hour
    # names -> (:date, :hour, :avg_hour, :raw)
    # row -> (yymmddTHH , HH, yymmdd, ..., ...)
    df_avg_hour = populate_periodic_mean(df, _df_avg_hour, :hour, :avg_hour)

    DataFrames.insertcols!(df_avg_hour, 1, :raw => df[!, target])
    # add date column for periodic mean computation
    # names(df_avg_day) = (:hour, :date, :date_day, :(target)_mean)
    DataFrames.insertcols!(df_avg_hour, 1, :date => DateTime.(df[!, :date], Local))

    # now daily seasonality is df_avg_hour[!, :avg_hour]
    day_sea = df_avg_hour[!, :avg_hour]
    day_res = df[!, target] .- df_avg_hour[!, :avg_hour]

    # now get annual seasonality from day_res1
    DataFrames.insertcols!(df_avg_hour, 3, :day_sea => day_sea)
    DataFrames.insertcols!(df_avg_hour, 3, :day_res => day_res)
    @show first(df_avg_hour, 30)

    df_avg_hour[!, :date_day] = Date.(df_avg_hour[!, :date])
    # Metaprogamming for ycol_mean = ycol => x -> mean
    # https://discourse.julialang.org/t/metaprogramming-function-with-keyword-argument/12551/2?u=appleparan
    z = Expr(:kw, :day_res_mean, :day_res => mean)
    #z2 = Expr(:kw, :date_day, :date => x -> Dates.Date.(DateTime.(x, Local)))
    # size(df_daymean, 1) == size(df, 1) / 24
    # 1 value per YYMMDD
    df_daymean = eval(:(DataFrames.by($(df_avg_hour), :date_day; $z)))
    # daily averaged values are stored
    # size(df_avg_day, 1) == size(df, 1)
    # names -> (:date, :hour, :date_day, target_daymean, :raw)
    # row -> (yymmddTHH , HH, YYMMDD, ..., ...)
    df_avg_day = populate_periodic_mean(df_avg_hour, df_daymean, :date_day, :day_res_mean)
    DataFrames.insertcols!(df_avg_day, 1, :hour => Dates.hour.(df_avg_hour[!, :date]))
    DataFrames.insertcols!(df_avg_day, 1, :date => df_avg_hour[!, :date])
    DataFrames.insertcols!(df_avg_day, 5, :raw => df_avg_hour[!, :raw])
    DataFrames.insertcols!(df_avg_day, 5, :day_res => df_avg_hour[!, :day_res])
    DataFrames.insertcols!(df_avg_day, 5, :day_sea => df_avg_hour[!, :day_sea])
    @show first(df_avg_day, 5)

    # construct keys like MMDDHH
    months = Dates.month.(df_avg_day[!, :date])
    days = Dates.day.(df_avg_day[!, :date])
    hours = Dates.hour.(df_avg_day[!, :date])

    day_keys = Symbol.(lpad.(string.(months), 2, "0") .* lpad.(string.(days), 2, "0") .* lpad.(string.(hours), 2, "0") )
    DataFrames.insertcols!(df_avg_day, 3, :day_keys => day_keys)

    # join and get seasonality
    _df_avg_year = DataFrames.by(df_avg_day, :day_keys, :day_res => mean)
    df_avg_year = populate_periodic_mean(df_avg_day, _df_avg_year, :day_keys, :day_res_mean)

    year_sea = df_avg_year[!, :day_res_mean]
    year_res = df_avg_day[!, :day_res] .- df_avg_year[!, :day_res_mean]

    # insert results (raw residual)
    DataFrames.insertcols!(df_avg_day, 3, :year_sea => year_sea)
    DataFrames.insertcols!(df_avg_day, 3, :year_res => year_res)

    # get also smoothed residual
    year_sea_s_x, year_sea_s = smoothing_series(float.(collect(1:length(df_avg_day[!, :year_sea]))), float.(df_avg_day[!, :year_sea]))
    year_res_s = df_avg_day[!, :day_res] .- year_sea_s
    # insert results (smoothed residual)
    DataFrames.insertcols!(df_avg_day, 3, :year_sea_s => year_sea_s)
    DataFrames.insertcols!(df_avg_day, 3, :year_res_s => year_res_s)
    @show first(df_avg_day, 5)

    df_avg_day
end

"""
    season_adj_lee(df::DataFrame, ycol::Symbol)
simple implementation for Prof. Lee's seasonal adjustment

Assume data as

y = smoothing(S_year) + S_day + residual

S_year : annual seasonality
S_day : daily seasonality
residual : residual

`S_year` and `S_day` is computed by mean
naive but easy to implement
"""
function season_adj_lee(df::DataFrame, ycol::Symbol)
    f_date = df[1, :date]
    t_date = df[end, :date]

    season_adj_lee(df, ycol, f_date, t_date)
end

"""
    season_adj_lee(df::DataFrame, ycol::Symbol, f_date::DateTime, t_date::DateTime)
Prof. Lee's seasonal adjustment with date range
station fluctuation will be averaged

Assume data as

y = smoothing(S_year) + S_day + residual

S_year : annual seasonality
S_day : daily seasonality
residual : residual

`S_year` and `S_day` is computed by mean
naive but easy to implement
"""
season_adj_lee(df::DataFrame, ycol::Symbol,
    f_date::ZonedDateTime, t_date::ZonedDateTime) =
    season_adj_lee(df, ycol, DateTime(f_date, Local), DateTime(t_date, Local))

function season_adj_lee(df::DataFrame, ycol::Symbol,
    f_date::DateTime, t_date::DateTime)

    # filter by given date range
    _df = @from i in df begin
        @where f_date <= DateTime(i.date, Local) <= t_date
        @orderby i.date
        @select i
        @collect DataFrame
    end

    # all station values are averaged by date and filter date, ycol, and mean_ycol only
    ycol_mean = Symbol(ycol, "_mean")

    # Metaprogamming for ycol_mean = ycol => x -> mean
    # https://discourse.julialang.org/t/metaprogramming-function-with-keyword-argument/12551/2?u=appleparan
    z = Expr(:kw, ycol_mean, ycol => mean)
    period_df = eval(:(DataFrames.by($(_df), :date; $z)))
    # insert original columns

    last_date = period_df[end, :date]
    window = 24

    # window mean
    year_wd_mean = compute_annual_mean(period_df, ycol_mean)
    # remove annual seasonality
    # fluctuation of annual average by 24 hour window
    year_res1 = zeros(length(period_df[!, :date]))
    year_sea1 = zeros(length(period_df[!, :date]))

    for (i, row) in enumerate(eachrow(period_df))
        t = getproperty(row, :date)
        if t > last_date - Dates.Hour(window)
            break
        end
        idx = (Dates.month(t) - 1) * 31 * 24 + (Dates.day(t) - 1) * 24 + Dates.hour(t) + 1

        year_res1[i] = getproperty(row, ycol_mean) - year_wd_mean[idx]
        year_sea1[i] = year_wd_mean[idx]
    end

    # smoothing(S_year)
    # The leap year error will be ignorable thanks to smoothing
    year_sea2_x, year_sea2_y = smoothing_series(float.(collect(1:length(year_sea1))), year_sea1; span = 0.1)

    #_model = Loess.loess(float.(collect(1:length(year_sea1))), year_sea1)
    #year_sea2_x = collect(range(extrema(float.(collect(1:length(year_sea1))))...; step = 1))
    #year_sea2_y = Loess.predict(_model, year_sea2_x)

    # Interpolations
    #itp = interpolate((year_sea2_x, year_sea2_y), BSpline(Cubic(Line(OnGrid()))))
    #year_sea2 = map(i -> itp(float(i)), 1:length(period_df[!, :date]))
    # Dierckx
    spl = Spline1D(year_sea2_x, year_sea2_y)
    year_sea2 = evaluate(spl, float.(collect(1:length(period_df[!, :date]))))
    year_res2 = period_df[!, ycol_mean] .- year_sea2
    DataFrames.insertcols!(period_df, 3, :year_sea2 => year_sea2)
    DataFrames.insertcols!(period_df, 3, :year_res2 => year_res2)

    # now calculate S_day
    hour_col = hour.(period_df[!, :date])
    DataFrames.insertcols!(period_df, 3, :hour => hour_col)

    # put mean values to dataframe, adjust index only for hour (0-23)
    day_sea1_summary, day_sea1 = populate_periodic_mean(period_df, :year_res2, :hour)
    DataFrames.insertcols!(period_df, 3, :day_sea1 => day_sea1)
    day_res1 = zeros(length(period_df[!, :date]))

    for (i, row) in enumerate(eachrow(period_df))
        t = getproperty(row, :date)
        if t > last_date - Dates.Hour(window)
            break
        end

        day_res1[i] = getproperty(row, :year_res2) - getproperty(row, :day_sea1)
    end

    # S_year, smoothed(S_year), S_day, R
    year_sea1, year_sea2, year_res2, day_sea1, day_res1
end
#=
"""
    bsts()

Decompose series by Bayesian structural time series models

```math
    yₜ = μₜ + τₜ + βₜ x + ϵₜ

    where ϵₜ ~ N(0, σ²_ϵ  )
```
μₜ : local linear trend
τₜ : seasonality
βₜ x : regression component
ϵₜ : noise

# Arguments
- level::Symbol, optional : Whether or not to include a level component.
- trend::Bool, optional : Whether or not to include a trend component. If true, level must not be irregular
- seaasonal::Integer, optional : The period of the seasonal component
- freq_seasonal::Array{NamedTuple,1} : Whether (and how) to model seasonal component(s) with trig
    form of Array of NamedTuple, like (period = some_value, harmonics = some_value)
    period must be present
- cycle::Bool : Whether or not to include a cycle component
- autoregressive::Integer : The order of the autoregressive component
- exog::Array{Float64,1} : Exogenous variables
- irregular::Bool : Whether or not to include an irregular component
- stochastic_level : Whether or not any level component is stochastic
- stochastic_trend : Whether or not any trend component is stochastic
- stochastic_seasonal : Whether or not any seasonal component is stochastic
- stochastic_freq_seasonal : Whether or not each seasonal component(s) is (are) stochastic
- stochastic_cycle : Whether or not any cycle component is stochastic.
- damped_cycle : Whether or not the cycle component is damped

https://www.statsmodels.org/dev/generated/statsmodels.tsa.statespace.structural.UnobservedComponents.html
"""
function bsts(;level::Symbol=:irregular, trend::Bool=false, seasonal::Integer=0,
    freq_seasonal::Array{NamedTuple,1}=[(period=0., harmonics=0.)],
    cycle::Bool=false, autoregressive::Integer=0,
    exog::Array{Float64,1}=[], irregular::Bool=false,
    stochastic_level::Bool=false, stochastic_seasonal::Bool=false,
    stochastic_freq_seasonal::Array{Bool, 1}, stochastic_cycle::Bool=false,
    damped_cycle::Bool=false, cycle_period_bounds::Tuple{Float64, Float64}=(0.5, Inf),
    use_exact_diffuse::Bool=false)

    # Local linear trend

end
=#


text:
using PyCall

@pyimport dolo

modelpy = dolo.yaml_import("examples/models/consumption_savings_iid.yaml")
modelpy.__exogenous__ = modelpy.exogenous.processes[1]

@pyimport dolo.algos.time_iteration_helpers as tih
Fpy = tih.Euler(modelpy)

using Dolo

model = Dolo.yaml_import("examples/models/consumption_savings_iid.yaml")

@time Dolo.time_iteration(model; maxit=500)

@time Dolo.improved_time_iteration(model)

F = Dolo.Euler(model)

# set equal dprocess

Dolo.n_inodes(F.dprocess,1)
Dolo.inode(F.dprocess,1,1)
Dolo.iweight(F.dprocess,1,1)

Fpy.dprocess.n_inodes(1)
Fpy.dprocess.iweight(1,0)

jpoints = [Dolo.inode(F.dprocess, 1, j) for j=1:Dolo.n_inodes(F.dprocess, 1)]
jweights = [Dolo.iweight(F.dprocess, 1, j) for j=1:Dolo.n_inodes(F.dprocess, 1)]

ppoints  = [Fpy.dprocess.inode(0,j) for j=0:6]
pweights  = [Fpy.dprocess.iweight(0,j) for j=0:6]

F.dprocess.integration_weights = pweights
F.dprocess.integration_nodes = copy(hcat(ppoints...)')

# check derivatives without complementarities

flatten(x0::Dolo.MSM) = cat([e for e in x0.data]...;dims=1)

# same residual
x0_py = Fpy.x0
respy = Fpy(x0_py, x0_py)
resj = F(F.x0, F.x0; ignore_constraints=false)
maximum(abs, (respy.data)[:] - flatten(resj))

# same jacobian
rpy, Jpy = Fpy.d_A(x0_py, x0_py)
Jj = Dolo.df_A(F, F.x0, F.x0)
diff = [Jj.data[n] - (Jpy.data)[n,:,:] for n=1:length(Jj.data)]
ddiff = cat(diff...; dims=1)
maximum(abs, ddiff)

# same step
delta_py = Jpy.solve(rpy)
delta_jl = Jj\resj
maximum(abs, (delta_py.data)[:] - flatten(delta_jl))

x0j = F.x0
x0p = Fpy.x0

x0j = x1j
x0p = x1p

x1j = x0j - delta_jl
x1p = x0p - delta_py
maximum(abs, (x1p.data)[:] - flatten(x1j))

res1j = F(x1j, x0j; ignore_constraints=false)
res1p = Fpy(x1p, x0p; ignore_constraints=false)

J1j = Dolo.df_A(F, x1j, x0j)
r1p, J1p = Fpy.d_A(x1p, x0p)

Dolo.norm(res1j)
res1p.norm()

maximum(abs, (res1p.data)[:] - flatten(res1j))

diff1 = [J1j.data[n] - (J1p.data)[n,:,:] for n=1:length(J1j.data)]
ddiff1 = cat(diff...; dims=1)
maximum(abs, ddiff1)

#####

time_iteration(model, maxit=100)

dolo.time_iteration(modelpy, maxit=100)

sol = Dolo.improved_time_iteration(model)

tab = Dolo.tabulate(model, sol.dr, :w)

using SimplePlots

plot(tab[:w], tab[:w], ylims=(0,2))
plot!(tab[:w], tab[:c])


text:
# test PartialPose3XYYaw

using RoME, CoordinateTransformations, Distributions, TransformUtils
using Base: Test

println("test x translation case")
# trivial cases first, orientation based tests below
wTx = Vector{AffineMap}(2)
wTx[1] = Translation(0.0,0,0) ∘ LinearMap(Quat(1.0, 0, 0, 0))
wTx[2] = Translation(10.0, 0, 0) ∘ LinearMap(Quat(1.0, 0, 0, 0))

## Recalculate XYH
# change toolbox
wTx1 = convert(SE3, wTx[1])
wTx2 = convert(SE3, wTx[2])
# get rotation to world frame for local level
wEx1 = convert(Euler, wTx1.R)
wEx1.Y = 0.0
wRlx1 = SE3(zeros(3), wEx1)
# wRx2 = deepcopy(wTx2); wRx2.t = zeros(3);
# Odometries
x1Tx2 = (wTx1\wTx2)
wRlx1Tx2 = wRlx1 * x1Tx2
vEx1_2 = veeEuler(wRlx1Tx2)
XYH1_2 = [vEx1_2[1], vEx1_2[2], vEx1_2[6]]
prz2 = veeEuler(wTx1)[[3;4;5]]

# test with Pose3Pose3
testpp3 = Pose3Pose3(MvNormal(veeEuler(x1Tx2),0.001*eye(6)))
res = zeros(6)
testpp3(res, 1, (veeEuler(x1Tx2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@test norm(res[1:3]) < 1e-10
@test norm(res[4:6]) < 1e-10

# test with PartialXYH
testppxyh = PartialPose3XYYaw(  MvNormal(XYH1_2,0.001*eye(3))  )
res = zeros(3)
testppxyh(res, 1, (vectoarr2(XYH1_2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@test norm(res[1:2]) < 1e-10
@test abs(res[3]) < 1e-10

println("test z translation case")
# z translation only
wTx = Vector{AffineMap}(2)
wTx[1] = Translation(0.0,0,0) ∘ LinearMap(Quat(1.0, 0, 0, 0))
wTx[2] = Translation(0, 0, 10.0) ∘ LinearMap(Quat(1.0, 0, 0, 0))

## Recalculate XYH
# change toolbox
wTx1 = convert(SE3, wTx[1])
wTx2 = convert(SE3, wTx[2])
# get rotation to world frame for local level
wEx1 = convert(Euler, wTx1.R)
wEx1.Y = 0.0
wRlx1 = SE3(zeros(3), wEx1)
# wRx2 = deepcopy(wTx2); wRx2.t = zeros(3);
# Odometries
x1Tx2 = (wTx1\wTx2)
wRlx1Tx2 = wRlx1 * x1Tx2
vEx1_2 = veeEuler(wRlx1Tx2)
XYH1_2 = [vEx1_2[1], vEx1_2[2], vEx1_2[6]]
prz2 = veeEuler(wTx1)[[3;4;5]]

# test with Pose3Pose3
testpp3 = Pose3Pose3(MvNormal(veeEuler(x1Tx2),0.001*eye(6)))
res = zeros(6)
testpp3(res, 1, (veeEuler(x1Tx2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@test norm(res[1:3]) < 1e-10
@test norm(res[4:6]) < 1e-10

# test with PartialXYH
testppxyh = PartialPose3XYYaw(  MvNormal(XYH1_2,0.001*eye(3))  )
res = zeros(3)
testppxyh(res, 1, (vectoarr2(XYH1_2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@test norm(res[1:2]) < 1e-10
@test abs(res[3]) < 1e-10

println("test roll and translate case 1")
# different orientation, roll
sq2 = 1.0/sqrt(2)
wTx[1] = Translation(0.0,0,0) ∘ LinearMap(Quat(sq2, sq2, 0, 0))
wTx[2] = Translation(10.0, 0, 0) ∘ LinearMap(Quat(sq2, sq2, 0, 0))

## Recalculate XYH
# change toolbox
wTx1 = convert(SE3, wTx[1])
wTx2 = convert(SE3, wTx[2])
# get rotation to world frame for local level
wEx1 = convert(Euler, wTx1.R)
wEx1.Y = 0.0
wRlx1 = SE3(zeros(3), wEx1)
# wRx2 = deepcopy(wTx2); wRx2.t = zeros(3);
# Odometries
x1Tx2 = (wTx1\wTx2)
wRlx1Tx2 = wRlx1 * x1Tx2
vEx1_2 = veeEuler(wRlx1Tx2)
XYH1_2 = [vEx1_2[1], vEx1_2[2], vEx1_2[6]]
prz2 = veeEuler(wTx1)[[3;4;5]]

# test with Pose3Pose3
testpp3 = Pose3Pose3(MvNormal(veeEuler(x1Tx2),0.001*eye(6)))
res = zeros(6)
testpp3(res, 1, (veeEuler(x1Tx2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@show res
@test norm(res[1:3]) < 1e-10
@test norm(res[4:6]) < 1e-10

# test with PartialXYH
testppxyh = PartialPose3XYYaw(  MvNormal(XYH1_2,0.001*eye(3))  )
res = zeros(3)
testppxyh(res, 1, (vectoarr2(XYH1_2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@test norm(res[1:2]) < 1e-10
@test abs(res[3]) < 1e-10

println("test roll and translate case 2")
# different orientation, roll
sq2 = 1.0/sqrt(2)
wTx[1] = Translation(0.0,0,0) ∘ LinearMap(Quat(sq2, sq2, 0, 0))
wTx[2] = Translation(0, 10.0, 0) ∘ LinearMap(Quat(sq2, sq2, 0, 0))

## Recalculate XYH
# change toolbox
wTx1 = convert(SE3, wTx[1])
wTx2 = convert(SE3, wTx[2])
# get rotation to world frame for local level
wEx1 = convert(Euler, wTx1.R)
wEx1.Y = 0.0
wRlx1 = SE3(zeros(3), wEx1)
# wRx2 = deepcopy(wTx2); wRx2.t = zeros(3);
# Odometries
x1Tx2 = (wTx1\wTx2)
wRlx1Tx2 = wRlx1 * x1Tx2
vEx1_2 = veeEuler(wRlx1Tx2)
XYH1_2 = [vEx1_2[1], vEx1_2[2], vEx1_2[6]]
prz2 = veeEuler(wTx1)[[3;4;5]]

# test with Pose3Pose3
testpp3 = Pose3Pose3(MvNormal(veeEuler(x1Tx2),0.001*eye(6)))
res = zeros(6)
testpp3(res, 1, (veeEuler(x1Tx2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@show res
@test norm(res[1:3]) < 1e-10
@test norm(res[4:6]) < 1e-10

# test with PartialXYH
testppxyh = PartialPose3XYYaw(  MvNormal(XYH1_2,0.001*eye(3))  )
res = zeros(3)
testppxyh(res, 1, (vectoarr2(XYH1_2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@test norm(res[1:2]) < 1e-10
@test abs(res[3]) < 1e-10

println("test pitch and translate case 1")
# different orientation, pitch
qq = convert(CoordinateTransformations.Quat, CoordinateTransformations.AngleAxis(pi/4,0,1.0,0))
wTx[1] = Translation(0.0,0,0) ∘ LinearMap(qq)
wTx[2] = Translation(10.0, 0, 0) ∘ LinearMap(qq)

## Recalculate XYH
# change toolbox
wTx1 = convert(SE3, wTx[1])
wTx2 = convert(SE3, wTx[2])
# get rotation to world frame for local level
wEx1 = convert(Euler, wTx1.R)
wEx1.Y = 0.0
wRlx1 = SE3(zeros(3), wEx1)
# wRx2 = deepcopy(wTx2); wRx2.t = zeros(3);
# Odometries
x1Tx2 = (wTx1\wTx2)
wRlx1Tx2 = wRlx1 * x1Tx2
vEx1_2 = veeEuler(wRlx1Tx2)
XYH1_2 = [vEx1_2[1], vEx1_2[2], vEx1_2[6]]
prz2 = veeEuler(wTx1)[[3;4;5]]

# test with Pose3Pose3
testpp3 = Pose3Pose3(MvNormal(veeEuler(x1Tx2),0.001*eye(6)))
res = zeros(6)
testpp3(res, 1, (veeEuler(x1Tx2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@show res
@test norm(res[1:3]) < 1e-10
@test norm(res[4:6]) < 1e-10

# test with PartialXYH
testppxyh = PartialPose3XYYaw(  MvNormal(XYH1_2,0.001*eye(3))  )
res = zeros(3)
testppxyh(res, 1, (vectoarr2(XYH1_2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@show res
@test norm(res[1:2]) < 1e-10
@test abs(res[3]) < 1e-10

println("test pitch and translate case 2")
# different orientation, pitch
qq = convert(CoordinateTransformations.Quat, CoordinateTransformations.AngleAxis(pi/4,0,1.0,0))
wTx[1] = Translation(0.0,0,0) ∘ LinearMap(qq)
wTx[2] = Translation(0, 10.0, 0) ∘ LinearMap(qq)

## Recalculate XYH
# change toolbox
wTx1 = convert(SE3, wTx[1])
wTx2 = convert(SE3, wTx[2])
# get rotation to world frame for local level
wEx1 = convert(Euler, wTx1.R)
wEx1.Y = 0.0
wRlx1 = SE3(zeros(3), wEx1)
# wRx2 = deepcopy(wTx2); wRx2.t = zeros(3);
# Odometries
x1Tx2 = (wTx1\wTx2)
wRlx1Tx2 = wRlx1 * x1Tx2
vEx1_2 = veeEuler(wRlx1Tx2)
XYH1_2 = [vEx1_2[1], vEx1_2[2], vEx1_2[6]]
prz2 = veeEuler(wTx1)[[3;4;5]]

# test with Pose3Pose3
testpp3 = Pose3Pose3(MvNormal(veeEuler(x1Tx2),0.001*eye(6)))
res = zeros(6)
testpp3(res, 1, (veeEuler(x1Tx2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@show res
@test norm(res[1:3]) < 1e-10
@test norm(res[4:6]) < 1e-10

# test with PartialXYH
testppxyh = PartialPose3XYYaw(  MvNormal(XYH1_2,0.001*eye(3))  )
res = zeros(3)
testppxyh(res, 1, (vectoarr2(XYH1_2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@show res
@test norm(res[1:2]) < 1e-10
@test abs(res[3]) < 1e-10

println("test pitch and translate case 3")
# different orientation, pitch
qq = convert(CoordinateTransformations.Quat, CoordinateTransformations.AngleAxis(pi/4,0,1.0,0))
wTx[1] = Translation(0.0,0,0) ∘ LinearMap(qq)
wTx[2] = Translation(0, 0, 10.0) ∘ LinearMap(qq)

## Recalculate XYH
# change toolbox
wTx1 = convert(SE3, wTx[1])
wTx2 = convert(SE3, wTx[2])
# get rotation to world frame for local level
wEx1 = convert(Euler, wTx1.R)
wEx1.Y = 0.0
wRlx1 = SE3(zeros(3), wEx1)
# wRx2 = deepcopy(wTx2); wRx2.t = zeros(3);
# Odometries
x1Tx2 = (wTx1\wTx2)
wRlx1Tx2 = wRlx1 * x1Tx2
vEx1_2 = veeEuler(wRlx1Tx2)
XYH1_2 = [vEx1_2[1], vEx1_2[2], vEx1_2[6]]
prz2 = veeEuler(wTx1)[[3;4;5]]

# test with Pose3Pose3
testpp3 = Pose3Pose3(MvNormal(veeEuler(x1Tx2),0.001*eye(6)))
res = zeros(6)
testpp3(res, 1, (veeEuler(x1Tx2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@show res
@test norm(res[1:3]) < 1e-10
@test norm(res[4:6]) < 1e-10

# test with PartialXYH
testppxyh = PartialPose3XYYaw(  MvNormal(XYH1_2,0.001*eye(3))  )
res = zeros(3)
testppxyh(res, 1, (vectoarr2(XYH1_2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@show res
@test norm(res[1:2]) < 1e-10
@test abs(res[3]) < 1e-10

println("test pitch and translate case 4")
# different orientation, pitch
qq = convert(CoordinateTransformations.Quat, CoordinateTransformations.AngleAxis(pi/4,0,1.0,0))
wTx[1] = Translation(0.0,0,0) ∘ LinearMap(qq)
wTx[2] = Translation(0, 15.0, 10.0) ∘ LinearMap(qq)

## Recalculate XYH
# change toolbox
wTx1 = convert(SE3, wTx[1])
wTx2 = convert(SE3, wTx[2])
# get rotation to world frame for local level
wEx1 = convert(Euler, wTx1.R)
wEx1.Y = 0.0
wRlx1 = SE3(zeros(3), wEx1)
# wRx2 = deepcopy(wTx2); wRx2.t = zeros(3);
# Odometries
x1Tx2 = (wTx1\wTx2)
wRlx1Tx2 = wRlx1 * x1Tx2
vEx1_2 = veeEuler(wRlx1Tx2)
XYH1_2 = [vEx1_2[1], vEx1_2[2], vEx1_2[6]]
prz2 = veeEuler(wTx1)[[3;4;5]]

# test with Pose3Pose3
testpp3 = Pose3Pose3(MvNormal(veeEuler(x1Tx2),0.001*eye(6)))
res = zeros(6)
testpp3(res, 1, (veeEuler(x1Tx2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@show res
@test norm(res[1:3]) < 1e-10
@test norm(res[4:6]) < 1e-10

# test with PartialXYH
testppxyh = PartialPose3XYYaw(  MvNormal(XYH1_2,0.001*eye(3))  )
res = zeros(3)
testppxyh(res, 1, (vectoarr2(XYH1_2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@show res
@test norm(res[1:2]) < 1e-10
@test abs(res[3]) < 1e-10

println("test yaw and translate case 1")
# different orientation, yaw
qq = convert(CoordinateTransformations.Quat, CoordinateTransformations.AngleAxis(pi/2,0,0,1.0))
wTx[1] = Translation(0.0,0,0) ∘ LinearMap(qq)
wTx[2] = Translation(10.0, 0, 0) ∘ LinearMap(qq)

## Recalculate XYH
# change toolbox
wTx1 = convert(SE3, wTx[1])
wTx2 = convert(SE3, wTx[2])
# get rotation to world frame for local level, free orientation
wEx1 = convert(Euler, wTx1.R)
wEx1.Y = 0.0
wRlx1 = SE3(zeros(3), wEx1)
# wRx2 = deepcopy(wTx2); wRx2.t = zeros(3);
# Odometries
x1Tx2 = (wTx1\wTx2)
wRlx1Tx2 = wRlx1 * x1Tx2
vEx1_2 = veeEuler(wRlx1Tx2)
XYH1_2 = [vEx1_2[1], vEx1_2[2], vEx1_2[6]]
prz2 = veeEuler(wTx1)[[3;4;5]]

# test with Pose3Pose3
testpp3 = Pose3Pose3(MvNormal(veeEuler(x1Tx2),0.001*eye(6)))
res = zeros(6)
testpp3(res, 1, (veeEuler(x1Tx2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@show res
@test norm(res[1:3]) < 1e-10
@test norm(res[4:6]) < 1e-10

# test with PartialXYH
testppxyh = PartialPose3XYYaw(  MvNormal(XYH1_2,0.001*eye(3))  )
res = zeros(3)
testppxyh(res, 1, (vectoarr2(XYH1_2),), vectoarr2(veeEuler(wTx1)), vectoarr2(veeEuler(wTx2)))
@show res
@test norm(res[1:2]) < 1e-10
@test abs(res[3]) < 1e-10

#


text:
# Refinement

@doc raw"""
Return a coefficient matrix A which satisfy
```math
B_{(i,p,k)} = \sum_{j}A_{i,j}B_{(j,p',k')}
```

Assumption:
* ``P ⊆ P′``
"""
function _changebasis_R(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)::Matrix{Float64}
    p = degree(P)
    k = knots(P)
    p′ = degree(P′)
    k′ = knots(P′)
    p₊ = p′ - p

    if p == 0
        n = length(k) - 1
        n′ = length(k′) - p₊ - 1
        A⁰ = Float64[bsplinesupport(typeof(P′)(p₊, k′), j) ⊆ bsplinesupport(typeof(P)(0, k), i) for i in 1:n, j in 1:n′]
        A⁰[:, findall(iszeros(P′))] .= NaN
        return A⁰
    end

    Aᵖ⁻¹ = _changebasis_R(typeof(P)(p-1, k), typeof(P′)(p′-1, k′)) # (n+1) × (n′+1) matrix
    n = dim(P)
    n′ = dim(P′)
    Z = iszeros(typeof(P′)(p′-1, k′))
    W = findall(Z)
    K′ = [k′[i+p′] - k′[i] for i in 1:n′+1]
    K = [ifelse(k[i+p] ≠ k[i], 1 / (k[i+p] - k[i]), 0.0) for i in 1:n+1]
    Δ = (p / p′) * [K′[j] * (K[i] * Aᵖ⁻¹[i, j] - K[i+1] * Aᵖ⁻¹[i+1, j]) for i in 1:n, j in 1:n′+1]
    Aᵖ = zeros(n, n′)
    Aᵖ[:, 1] = Δ[:, 1]
    Aᵖ[:, n′] = -Δ[:, n′+1]

    # split Aᵖ for sub-block
    if length(W) == 0
        Q = [1:n′]
    else
        Q = [1:W[1]-1, [W[i]:W[i+1]-1 for i in 1:length(W)-1]..., W[end]:n′]
    end
    λ = length(Q)
    Λ = length.(Q)
    Ãᵖ = [Aᵖ[:, q] for q in Q]

    for ȷ in 2:λ-1
        if Λ[ȷ] == 1
            # if B(i,p′,k′) = 0
            Ãᵖ[ȷ] .= NaN
        end
    end
    for ȷ in 1:λ-1
        if Λ[ȷ] ≥ 2
            t = k′[W[ȷ]]
            for i in 1:n
                Ãᵖ[ȷ][i, end] = bsplinebasis₋₀(P, i, t)
            end
        end
    end
    for ȷ in 2:λ
        if Λ[ȷ] ≥ 2
            t = k′[W[ȷ-1]+p]
            for i in 1:n
                Ãᵖ[ȷ][i, 1] = bsplinebasis₊₀(P, i, t)
            end
        end
    end
    for ȷ in 1:λ
        if Λ[ȷ] ≥ 3
            r = Q[ȷ]
            A₊ = copy(Ãᵖ[ȷ])
            A₋ = copy(Ãᵖ[ȷ])
            for j in 1:Λ[ȷ]-2
                A₊[:, j+1] = A₊[:, j] + Δ[:, j+r[1]]
                A₋[:, Λ[ȷ]-j] = A₋[:, Λ[ȷ]-j+1] - Δ[:, Λ[ȷ]-j+r[1]]
            end
            Ãᵖ[ȷ] = (A₊ + A₋) / 2
        end
    end
    Aᵖ = hcat(Ãᵖ...) # n × n′ matrix
    return Aᵖ .* Float64[bsplinesupport(P′,j) ⊆ bsplinesupport(P,i) for i in 1:n, j in 1:n′]
end

@doc raw"""
Return a coefficient matrix A which satisfy
```math
B_{(i,p_1,k_1)} = \sum_{j}A_{i,j}B_{(j,p_2,k_2)}
```

Assumption:
* ``P1 ⊑ P2``
* ``P2 ⊑ P1``
"""
function _changebasis_sim(P1::AbstractBSplineSpace, P2::AbstractBSplineSpace)
    # if P1 ⋢ P2
    #     error("P1 ⋢ P2")
    # end
    # if P2 ⋢ P1
    #     error("P2 ⋢ P1")
    # end
    n = dim(P1)
    p = degree(P1)
    v = (knots(P1).vector)[p+1:end-p]

    if length(v) > p
        A = Matrix{Float64}(I, n, n)
        # TODO: Fix below
        vvv = [v[1] * (p - i + 1) / (p + 1) + v[i+1] * (i) / (p + 1) for i in 1:p]
        A1 = [bsplinebasis₊₀(P1, i, t) for i in 1:p, t in vvv]
        A2 = [bsplinebasis₊₀(P2, i, t) for i in 1:p, t in vvv]
        A[1:p, 1:p] = A1 * inv(A2)
        vvv = [v[end-p+i-1] * (p - i + 1) / (p + 1) + v[end] * (i) / (p + 1) for i in 1:p]
        A1 = [bsplinebasis₋₀(P1, i, t) for i in n-p+1:n, t in vvv]
        A2 = [bsplinebasis₋₀(P2, i, t) for i in n-p+1:n, t in vvv]
        A[n-p+1:n, n-p+1:n] = A1 * inv(A2)
        # TODO: Fix above
    else
        # TODO: Fix below
        vvv = [v[1] * (n - i + 1) / (n + 1) + v[end] * (i) / (n + 1) for i in 1:n]
        A1 = [bsplinebasis₋₀(P1, i, t) for i in 1:n, t in vvv]
        A2 = [bsplinebasis₋₀(P2, i, t) for i in 1:n, t in vvv]
        A = A1 * inv(A2)
        # TODO: Fix above
    end
    return A
end

@doc raw"""
Return a coefficient matrix A which satisfy
```math
B_{(i,p,k)} = \sum_{j}A_{i,j}B_{(j,p',k')}
```

Assumption:
* ``P ⊑ P′``
"""
function _changebasis_I(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)::Matrix{Float64}
    I = bsplineunity(P)
    p = degree(P)
    k = knots(P)
    p′ = degree(P′)
    k′ = knots(P′)
    p₊ = p′ - p

    _P = typeof(P)(p, k[1+p:end-p] + p * Knots(k[1+p], k[end-p]))
    # if dim(_P) ≠ dim(P)
    #     error("dim(_P)≠dim(P)")
    # end
    _P′ = typeof(P′)(p′, k′[1+p′:end-p′] + p′ * Knots(k′[1+p′], k′[end-p′]))
    # if dim(_P′) ≠ dim(P′)
    #     error("dim(_P′)≠dim(P′)")
    # end
    _A = _changebasis_R(_P, _P′)
    Asim = _changebasis_sim(P, _P)
    Asim′ = _changebasis_sim(_P′, P′)
    A = Asim * _A * Asim′

    return A
end

function changebasis(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)
    if P ⊆ P′
        return _changebasis_R(P, P′)
    elseif P ⊑ P′
        return _changebasis_I(P, P′)
    else
        throw(DomainError((P, P′),"𝒫[p,k] ⊆ 𝒫[p′,k′] or 𝒫[p,k] ⊑ 𝒫[p′,k′] must hold."))
    end
end

@doc raw"""
Refinement of B-spline manifold with given B-spline spaces.
"""
function refinement(M::AbstractBSplineManifold, Ps′::Vector{<:AbstractBSplineSpace})
    Ps = collect(bsplinespaces(M))
    a = controlpoints(M)
    d = length(Ps)
    n = dim.(Ps)
    n′ = dim.(Ps′)

    A = changebasis.(Ps, Ps′)
    # TODO: general dimension
    if d == 1
        a′ = [sum(A[1][I₁, J₁] * a[I₁] for I₁ in 1:n[1]) for J₁ in 1:n′[1]]
    elseif d == 2
        a′ = [sum(A[1][I₁, J₁] * A[2][I₂, J₂] * a[I₁, I₂] for I₁ in 1:n[1], I₂ in 1:n[2]) for J₁ in 1:n′[1], J₂ in 1:n′[2]]
    elseif d == 3
        a′ = [
            sum(A[1][I₁, J₁] * A[2][I₂, J₂] * A[3][I₃, J₃] * a[I₁, I₂, I₃] for I₁ in 1:n[1], I₂ in 1:n[2], I₃ in 1:n[3])
            for J₁ in 1:n′[1], J₂ in 1:n′[2], J₃ in 1:n′[3]
        ]
    end
    return (typeof(M).name.wrapper)(Ps′, a′)
end

@doc raw"""
Refinement of B-spline manifold with additional degree and knots.
"""
function refinement(M::AbstractBSplineManifold; p₊::Union{Nothing,AbstractVector{<:Integer}} = nothing, k₊::Union{Nothing,Vector{Knots}} = nothing)
    Ps = collect(bsplinespaces(M))
    𝒂 = controlpoints(M)
    d = length(Ps)
    n = dim.(Ps)
    if isnothing(p₊)
        p₊ = zeros(Int, d)
    elseif length(Ps) ≠ length(p₊)
        throw(DimensionMismatch())
    end
    if isnothing(k₊)
        k₊ = zeros(Knots, d)
    elseif length(Ps) ≠ length(k₊)
        throw(DimensionMismatch())
    end

    Ps′ = [(P = Ps[i];
    p = degree(P);
    k = knots(P);
    k_unique = unique(k[1+p:end-p]);
    typeof(P)(p + p₊[i], k + p₊[i] * k_unique + k₊[i])) for i in eachindex(Ps)]

    return refinement(M, Ps′)
end


text:
function tensorvsmatrix(size)
    arraysize = Base.summarysize(fill(0, (size)))
    dict = Dict()
    if Base.summarysize(dict) > arraysize
        return -1
    end
    for i in 1:size
        dict[i] = i*1.0
        if Base.summarysize(dict) > arraysize
            return i
        end
    end
    return size + 1
end

println("Start")

x = 10
while x < 100001
    global x
    println((x, tensorvsmatrix(x)))
    x *= 10
end

#=
10, -1)
(100, 11)
(1000, 158)
(10000, 640)
(100000, 10923)
=#

max_stars_repo_name: Sohva/oolong

max_issues_repo_name: Sohva/oolong

max_forks_repo_name: oolongtensor/oolong


text:
struct CorrelatedOUUnknown{K} <: UnknownAcvTimeSeriesModel{2,Float64}
    σ::Float64
    θ::Float64
    ρ::Float64
    σ²::Float64
    θ²::Float64
    function CorrelatedOUUnknown{K}(σ,θ,ρ) where {K}
        σ > 0 || throw(ArgumentError("OU process requires 0 < σ."))
        θ > 0 || throw(ArgumentError("OU process requires 0 < θ."))
        1>ρ>0 || throw(ArgumentError("OU process requires 0 < ρ < 1."))
        new{K}(σ,θ,ρ,σ^2,θ^2)
    end
function CorrelatedOUUnknown{K}(x::Vector{Float64}) where {K}
        @boundscheck checkparameterlength(x,CorrelatedOUUnknown{K})
        @inbounds CorrelatedOUUnknown{K}(x[1], x[2], x[3])
    end
end
nalias(::CorrelatedOUUnknown{K}) where {K} = K
npars(::Type{CorrelatedOUUnknown{K}}) where {K} = 3
lowerbounds(::Type{CorrelatedOUUnknown{K}}) where {K} = lowerbounds(CorrelatedOU)
upperbounds(::Type{CorrelatedOUUnknown{K}}) where {K} = upperbounds(CorrelatedOU)

@propagate_inbounds function add_sdf!(out, m::CorrelatedOUUnknown, ω)
    @boundscheck checkbounds(out,1:3)
    s = m.σ²/(π*(m.θ²+ω^2))
    @inbounds begin
        out[1] += s
        out[2] += s*m.ρ
        out[3] += s
    end
end

@propagate_inbounds function grad_add_sdf!(out, m::CorrelatedOUUnknown, ω)
    @boundscheck checkbounds(out,1:3,1:3)
    @inbounds begin
        σ, θ, ρ = m.σ, m.θ, m.ρ
        σ_part = 2σ / (π*(m.θ²+ω^2))
        # ∂σ
        out[1,1] += σ_part
        out[2,1] += σ_part * ρ
        out[3,1] += σ_part
        # ∂θ
        θ_part = -π*θ/2*(σ_part)^2
        out[1,2] += θ_part
        out[2,2] += θ_part * ρ
        out[3,2] += θ_part
        # ∂ρ
        # out[1,3] += 0
        out[2,3] += σ_part * σ / 2
        # out[3,3] += 0
    end
    nothing
end

@propagate_inbounds function hess_add_sdf!(out, m::CorrelatedOUUnknown, ω)
    @boundscheck checkbounds(out,1:3,1:5)
    @inbounds begin
        σ, θ, ρ = m.σ, m.θ, m.ρ
        σ_part = 2 / (π*(m.θ²+ω^2))
        θσ_part = -π*θ*σ*(σ_part)^2
        θ_part = -π*σ_part*σ*(σ_part*σ/2 + θ*θσ_part)

        # ∂σ²
        out[1,1] += σ_part
        out[2,1] += σ_part * ρ
        out[3,1] += σ_part

        # ∂σθ
        out[1,2] += θσ_part
        out[2,2] += θσ_part * ρ
        out[3,2] += θσ_part

        #∂σρ
        # out[1,3] += 0
        out[2,3] += σ_part * σ
        # out[3,3] += 0

        # ∂θ²
        out[1,4] += θ_part
        out[2,4] += θ_part * ρ
        out[3,4] += θ_part

        #∂θρ
        # out[1,5] += 0
        out[2,5] += σ/2 * θσ_part
        # out[3,5] += 0

        #∂ρ² all zero
    end

    nothing
end


text:
function testRenderPass(device, glfw, passRef)
    pass = device.createRenderPass(lava.createSimpleForward(lava.RenderPassCreateInfo, glfw.format()))
    pass.setClearColor(vk.VkClearColorValue((0, 0, 0, 0)))
    pass.setClearColors([vk.VkClearColorValue((0, 0, 0, 0))])
    pass.setClearDepthStencil(vk.VkClearDepthStencilValue(0, 0))
    passRef[] = pass
    return true
end


text:
function generateSyntheticMesocosmData(;p = [1e-3, 1e3, 1.6, 0.0, 3.344], hist = 0.0, C₀ = [4.0], t₀ = 0.0, ⏰ = 1e3, reactionKinetics = "firstOrder", shape = "power law", resolution, standardDeviation)
      # Define differential equation
      @sam_model(reactionKinetics, shape)
      # Make a tuple containing the start and end times of the simulation
      tspan = (t₀, ⏰)
      # Specify condition before simulation start
      h(p,t) = hist
      # Build the DDE problem
      prob = DDEProblem(sam_model!, C₀, h, tspan, p)
      # Solve problem
      sol = solve(prob, MethodOfSteps(RK4()), saveat = resolution)
      # Add noise
      synthData = sol + standardDeviation*randn(size(Array(sol)))
      return(synthData)
end

max_stars_repo_name: elizabethjmohr/hyporheicBGC

max_issues_repo_name: elizabethjmohr/hyporheicBGC

max_forks_repo_name: elizabethjmohr/hyporheicBGC


text:
import SKTB
using JuLIP, JuLIP.ASE

TB=SKTB
β, eF, fixed_eF = 1.0, 0.0, true

tests = Any[]

# precompile
at = bulk("Si", pbc = false, cubic=true)
tbm = TB.NRLTB.NRLTBModel(:Si, TB.FermiDiracSmearing(β, eF, fixed_eF),
                           bzquad = TB.GammaPoint() )
energy(tbm, at);
forces(tbm, at);
tbm = TB.NRLTB.NRLTBModel(:Si, TB.FermiDiracSmearing(β, eF, fixed_eF),
                           bzquad = TB.MPGrid(at, (2, 2, 2)) )
rattle!(at, 0.01)
energy(tbm, at);
forces(tbm, at);

test1 = "Large Cluster"
at1 = (3,3,3) * bulk("Si", pbc=false, cubic=true)
tbm1 = TB.NRLTB.NRLTBModel(:Si, TB.FermiDiracSmearing(β, eF, fixed_eF),
                           bzquad = TB.GammaPoint() )
push!(tests, (test1, at1, tbm1))

test2 = "2D Cluster + k in z direction"
at2 = (4,4,1) * bulk("Si", pbc=false, cubic=true)
tbm2 = TB.NRLTB.NRLTBModel(:Si, TB.FermiDiracSmearing(β, eF, fixed_eF),
                           bzquad=TB.MPGrid(at2, (0, 0, 8)) )
push!(tests, (test2, at2, tbm2))

# using ProfileView
# @profile energy(tbm2, at2)
# ProfileView.view()

for (test, at, tbm) in tests
   println("====================================================")
   println(test)
   println("Energy: " )
   for n = 1:4
      n <= 2 && rattle!(at, 0.001)
      @time energy(tbm, at)
   end
   println("Forces: " )
   for n = 1:4
      n <= 2 && rattle!(at, 0.001)
      @time forces(tbm, at)
   end
end
println("====================================================")

####### CONCLUSION OF THIS BENCHMARK:
#
#  Basically, the implementation is cleaner and needs less memory
#  but all the real time is spent in eig - so very little is
#  gained for large systems. The main advantage of the new framework
#  is that we can now readily construct block-space matrices,
#  sparse matrices, full matrices, all without extra cost and less
#  storage.
#
#
# # skh = TB.update!(at2, tbm2.H)
# # @time TB.update!(at2, tbm2.H)
# # @time TB.update!(at2, tbm2.H)
# # @time TB.update!(at2, tbm2.H)
#
# println("Old Hamiltonian > sparse > full")
# @time map(full, hamiltonian(tbm2, at2))
# @time map(full, hamiltonian(tbm2, at2))
# @time map(full, hamiltonian(tbm2, at2))
#
# println("SKH")
# @time skh = TB.SparseSKH(tbm2.H, at2)
# @time skh = TB.SparseSKH(tbm2.H, at2)
# @time skh = TB.SparseSKH(tbm2.H, at2)
#
# println("alloc")
# @time wrk = TB._alloc_full(skh)
# @time wrk = TB._alloc_full(skh)
# @time wrk = TB._alloc_full(skh)
# Hout = wrk[1]
# Mout = wrk[2]
# k = JVecF(zeros(3))
#
# # println("_full!")
# # @time TB._full!(Hout, Mout, skh, k, tbm.H)
# # @time TB._full!(Hout, Mout, skh, k, tbm.H)
# # @time TB._full!(Hout, Mout, skh, k, tbm.H)
#
# println("full!")
# @time TB.full!(wrk, skh)
# @time TB.full!(wrk, skh)
# @time TB.full!(wrk, skh)
#
# println("eig")
# @time eig(Hout, Mout)
# @time eig(Hout, Mout)
# @time eig(Hout, Mout)

# @code_warntype TB._full!(Hout, Mout, skh, k, tbm.H)


text:
function a()
    x = 2
    y = 3
    z = + x y
    print(z)
end

max_stars_repo_name: tylerholmes/Interpreter

max_issues_repo_name: tylerholmes/Interpreter

max_forks_repo_name: tylerholmes/Interpreter


text:
makepositive(x) = exp(x)

invmakepositive(y) = log(y)

sigmoid(x) = 1.0 / (exp(-x) + 1.0)

invsigmoid(y) = - log(1/y -1.0)

transformbetween(x, lower, upper) = sigmoid(x) * (upper - lower) + lower

invtransformbetween(y, lower, upper) = invsigmoid((y - lower)/(upper-lower))


text:
using Pkg
Pkg.activate(joinpath(@__DIR__,".."))
using LinearAlgebra
using BenchmarkTools, Test
using Cthulhu
using Plots ; theme(:dark)
using SmoothedSpectralAbscissa ; const SSA = SmoothedSpectralAbscissa

# using OrdinaryDiffEq
using StochasticDiffEq

##

function easydyn(mat::Matrix{R},x0::Vector{R},tmax::Real,dt::Real=0.01) where R
  ts = range(0,tmax;step=dt)
  ret = Matrix{R}(undef,length(x0),length(ts))
  for (k,t) in enumerate(ts)
    ret[:,k]=exp(mat.*t)*x0
  end
  retnrm = mapslices(norm,ret;dims=1)[:]
  return ts,ret,retnrm
end

function lessnormal(mat::Matrix{R},d::R) where R
  n=size(mat,1)
  sh = schur(mat)
  upd = diagm(0=>fill(1.0,n),1=>fill(d,n-1))
  return sh.vectors*upd*sh.Schur*inv(upd)*sh.vectors'
end

function dyn_noise(mat::Matrix{R},
    x0::Vector{R},noiselevel::Real,t_end::Real;
    stepsize::Real=0.05,verbose=false) where R<:Real
  f(du,u,p,t) = let _du = mat*u ; copy!(du,_du) ; end
  σ_f(du,u,p,t) = fill!(du,noiselevel)
  prob = SDEProblem(f,σ_f,x0,(0.,t_end))
  solv =  solve(prob,EM();verbose=verbose,dt=stepsize)
  ret_u = hcat(solv.u...)
  retnrm = mapslices(norm,ret_u;dims=1)[:]
  return solv.t,ret_u,retnrm
end
##

n = 15
idm =diagm(0=>fill(1.0,n))
matrand=randn(n,n) ./ sqrt(n) - 1.1I
matrandln = lessnormal(matrand,1.00001)
alloc=SSA.SSAAlloc(n)

##
SSA.ssa!(matrand,nothing,alloc)
SSA.ssa!(matrandln,nothing,alloc)

##

@btime  SSA.ssa!($matrand,nothing,$alloc)
@btime  SSA.ssa!($matrandln,nothing,$alloc)
println("\n\n")
@btime  SSA.ssa!($matrand,nothing,$alloc;optim_method=SSA.OptimNewton)
@btime  SSA.ssa!($matrandln,nothing,$alloc;optim_method=SSA.OptimNewton )

##

@descend_code_warntype SSA.ssa!(matrand,nothing,alloc)


text:
test_gage(GeometricCenter,
    Array[
        [1.111, 1., 2.],
        [3.333, 3., 4.],
    ],
    [2.222, 2., 3.],
    "test 2: for array of arrays")


text:
baremodule ThreadLocalCounters

export @tlc

macro tlc end
function list end
function clear end

module Internal

using ..ThreadLocalCounters: ThreadLocalCounters
import ..ThreadLocalCounters: @tlc

import StaticStorages
import RecordArrays

include("internal.jl")

# Use README as the docstring of the module:
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end ThreadLocalCounters

end  # module Internal

end  # baremodule ThreadLocalCounters


text:
@inbounds function _normalize(point::SVector{n}, min_bound::SVector{n}, max_bound::SVector{n}, scale::Real) where {n}
    return (point .- min_bound) ./ scale
end
@inbounds function _normalize(spl::AbstractNormalSpline{n}, point::SVector{n}) where {n}
    return _normalize(point, _get_min_bound(spl), _get_max_bound(spl), _get_scale(spl))
end

@inbounds function _unnormalize(point::SVector{n}, min_bound::SVector{n}, max_bound::SVector{n}, scale::Real) where {n}
    return min_bound .+ scale .* point
end
@inbounds function _unnormalize(spl::AbstractNormalSpline{n}, point::SVector{n}) where {n}
    return _unnormalize(point, _get_min_bound(spl), _get_max_bound(spl), _get_scale(spl))
end

function _normalization_scaling(nodes::AbstractVecOfSVecs)
    min_bound = reduce((x, y) -> min.(x, y), nodes; init = fill(+Inf, eltype(nodes)))
    max_bound = reduce((x, y) -> max.(x, y), nodes; init = fill(-Inf, eltype(nodes)))
    scale = maximum(max_bound .- min_bound)
    return min_bound, max_bound, scale
end

function _normalization_scaling(nodes::AbstractVecOfSVecs, d_nodes::AbstractVecOfSVecs)
    min_bound = min.(reduce((x, y) -> min.(x, y), nodes; init = fill(+Inf, eltype(nodes))), reduce((x, y) -> min.(x, y), d_nodes; init = fill(+Inf, eltype(d_nodes))))
    max_bound = max.(reduce((x, y) -> max.(x, y), nodes; init = fill(-Inf, eltype(nodes))), reduce((x, y) -> max.(x, y), d_nodes; init = fill(-Inf, eltype(d_nodes))))
    scale = maximum(max_bound .- min_bound)
    return min_bound, max_bound, scale
end

function _estimate_accuracy(spl::AbstractNormalSpline{n,T,RK}) where {n,T,RK <: ReproducingKernel_0}
    vmax = maximum(abs, _get_values(spl))
    rmae = zero(T)
    @inbounds for i in 1:length(_get_nodes(spl))
        point = _unnormalize(spl, _get_nodes(spl)[i])
        σ     = _evaluate(spl, point)
        rmae  = max(rmae, abs(_get_values(spl)[i] - σ))
    end
    if vmax > 0
        rmae /= vmax
    end
    rmae = max(rmae, eps(T))
    res  = -floor(log10(rmae)) - 1
    res  = max(res, 0)
    return trunc(Int, res)
end

function _pairwise_sum_norms(nodes::AbstractVecOfSVecs{n,T}) where {n,T}
    ℓ = zero(T)
    @inbounds for i in 1:length(nodes), j in i:length(nodes)
        ℓ += norm(nodes[i] .- nodes[j])
    end
    return ℓ
end

function _pairwise_sum_norms_weighted(nodes::AbstractVecOfSVecs{n,T}, d_nodes::AbstractVecOfSVecs{n,T}, w_d_nodes::T) where {n,T}
    ℓ = zero(T)
    @inbounds for i in 1:length(nodes), j in 1:length(d_nodes)
        ℓ += norm(nodes[i] .- w_d_nodes .* d_nodes[j])
    end
    return ℓ
end

@inline _ε_factor(::RK_H0, ε::T) where {T} = one(T)
@inline _ε_factor(::RK_H1, ε::T) where {T} = T(3)/2
@inline _ε_factor(::RK_H2, ε::T) where {T} = T(2)

@inline _ε_factor_d(::RK_H0, ε::T) where {T} = one(T)
@inline _ε_factor_d(::RK_H1, ε::T) where {T} = T(2)
@inline _ε_factor_d(::RK_H2, ε::T) where {T} = T(5)/2

function _estimate_ε(k::ReproducingKernel_0, nodes)
    ε  = _estimate_ε(nodes)
    ε *= _ε_factor(k, ε)
    k  = typeof(k)(ε)
end

function _estimate_ε(k::ReproducingKernel_0, nodes, d_nodes)
    ε  = _estimate_ε(nodes, d_nodes)
    ε *= _ε_factor_d(k, ε)
    k  = typeof(k)(ε)
end

function _estimate_ε(nodes::AbstractVecOfSVecs{n,T}) where {n,T}
    n₁ = length(nodes)
    ε  = _pairwise_sum_norms(nodes)
    return ε > 0 ? ε * T(n)^T(inv(n)) / T(n₁)^(T(5) / 3) : one(T)
end

function _estimate_ε(nodes::AbstractVecOfSVecs{n,T}, d_nodes::AbstractVecOfSVecs{n,T}, w_d_nodes::T=T(0.1)) where {n,T}
    n₁ = length(nodes)
    n₂ = length(d_nodes)
    ε  = _pairwise_sum_norms(nodes) + _pairwise_sum_norms_weighted(nodes, d_nodes, w_d_nodes) + w_d_nodes * _pairwise_sum_norms(d_nodes)
    return ε > 0 ? ε * T(n)^T(inv(n)) / T(n₁ + n₂)^(T(5) / 3) : one(T)
end

function _estimate_epsilon(nodes::AbstractVecOfSVecs, kernel::ReproducingKernel_0)
    min_bound, max_bound, scale = _normalization_scaling(nodes)
    nodes = _normalize.(nodes, (min_bound,), (max_bound,), scale)
    ε     = _estimate_ε(nodes)
    ε    *= _ε_factor(kernel, ε)
    return ε
end

function _estimate_epsilon(nodes::AbstractVecOfSVecs, d_nodes::AbstractVecOfSVecs, kernel::ReproducingKernel_1)
    min_bound, max_bound, scale = _normalization_scaling(nodes, d_nodes)
    nodes   = _normalize.(nodes, (min_bound,), (max_bound,), scale)
    d_nodes = _normalize.(d_nodes, (min_bound,), (max_bound,), scale)
    ε       = _estimate_ε(nodes, d_nodes)
    ε      *= _ε_factor_d(kernel, ε)
    return ε
end

function _get_gram(nodes::AbstractVecOfSVecs, kernel::ReproducingKernel_0)
    min_bound, max_bound, scale = _normalization_scaling(nodes)
    nodes = _normalize.(nodes, (min_bound,), (max_bound,), scale)
    if kernel.ε == 0
        kernel = _estimate_ε(kernel, nodes)
    end
    return _gram(nodes, kernel)
end

function _get_gram(nodes::AbstractVecOfSVecs, d_nodes::AbstractVecOfSVecs, d_dirs::AbstractVecOfSVecs, kernel::ReproducingKernel_1)
    min_bound, max_bound, scale = _normalization_scaling(nodes, d_nodes)
    nodes   = _normalize.(nodes, (min_bound,), (max_bound,), scale)
    d_nodes = _normalize.(d_nodes, (min_bound,), (max_bound,), scale)
    d_dirs  = d_dirs ./ norm.(d_dirs)
    if kernel.ε == 0
        kernel = _estimate_ε(kernel, nodes, d_nodes)
    end
    return _gram(nodes, d_nodes, d_dirs, kernel)
end

function _get_cond(nodes::AbstractVecOfSVecs, kernel::ReproducingKernel_0)
    T = promote_type(eltype(kernel), eltype(eltype(nodes)))
    gram = _get_gram(nodes, kernel)
    cond = zero(T)
    try
        evs = svdvals!(gram)
        maxevs = maximum(evs)
        minevs = minimum(evs)
        if minevs > 0
            cond = maxevs / minevs
            cond = T(10)^floor(log10(cond))
        end
    catch
    end
    return cond
end

function _get_cond(nodes::AbstractVecOfSVecs, d_nodes::AbstractVecOfSVecs, d_dirs::AbstractVecOfSVecs, kernel::ReproducingKernel_1)
    _get_cond(nodes, kernel)
end

"""
Get estimation of the Gram matrix condition number
Brás, C.P., Hager, W.W. & Júdice, J.J. An investigation of feasible descent algorithms for estimating the condition number of a matrix. TOP 20, 791–809 (2012).
https://link.springer.com/article/10.1007/s11750-010-0161-9
"""
function _estimate_cond(A::AbstractMatrix{T}, F::Factorization, nit = 3) where {T}
    Anorm = norm(A, 1)
    n = size(A, 1)
    x = fill(inv(T(n)), n)
    z = zeros(T, n)
    z′ = zeros(T, n)
    γ = zero(T)
    @inbounds for _ in 1:nit
        ldiv!(z′, F, x)
        γ = zero(T)
        for i in 1:n
            γ += abs(z′[i])
            z′[i] = sign(z′[i])
        end
        ldiv!(z, F, z′)
        zdotx = z ⋅ x
        zmax, imax = T(-Inf), 1
        for i in 1:n
            zᵢ = z[i] = abs(z[i])
            if zᵢ > zmax
                zmax, imax = zᵢ, i
            end
        end
        (zmax <= zdotx) && break
        x .= 0
        x[imax] = 1
    end
    cond = prevpow(T(10), Anorm * γ)
    return cond
end


text:
using BenchmarkTools
using LoopVectorization
using LinearAlgebra
using Plots
using UnPack

using EntropyStableEuler.Fluxes2D
import EntropyStableEuler: γ

using FluxDiffUtils
using StartUpDG
using StartUpDG.ExplicitTimestepUtils

include("HybridizedSBPUtils.jl")
using .HybridizedSBPUtils

N = 2
K1D = 8
T = .25
CFL = .5

# init ref element and mesh
rd = init_reference_tri(N)
VX,VY,EToV = uniform_tri_mesh(K1D)
md = init_DG_mesh(VX,VY,EToV,rd)

# Make domain periodic
LX,LY = (x->maximum(x)-minimum(x)).((VX,VY)) # find lengths of domain
build_periodic_boundary_maps!(md,rd,LX,LY)

#####
##### Define initial coefficients and time-stepping
#####

# "Define the initial conditions by interpolation"
@unpack x,y = md
rho = @. 2 + .5*exp(-100*(x^2+y^2))
u = @. 0*x
v = @. 0*x
p = @. rho^γ

Q = primitive_to_conservative(rho,u,v,p)
Qrhskew,Qshskew,Vh,Ph,VhP = build_hSBP_ops(rd)
SBP_ops = (Matrix(Qrhskew'),Matrix(Qshskew'),Vh,Ph,VhP)

# interpolate geofacs to both vol/surf nodes, pack back into md
@unpack Vq,Vf = rd
@unpack rxJ, sxJ, ryJ, syJ = md
rxJ, sxJ, ryJ, syJ = (x->[Vq;Vf]*x).((rxJ, sxJ, ryJ, syJ)) # interp to hybridized points
@pack! md = rxJ, sxJ, ryJ, syJ

Ax,Ay,Bx,By,B = assemble_global_SBP_matrices_2D(rd, md, Qrhskew, Qshskew)

function euler_fluxes_conservative_vars(UL,UR)
    QL = conservative_to_primitive_beta(UL[1],UL[2],UL[3],UL[4])
    QR = conservative_to_primitive_beta(UR[1],UR[2],UR[3],UR[4])
    return SVector{4}(euler_fluxes(QL[1],QL[2],QL[3],QL[4],
                                   QR[1],QR[2],QR[3],QR[4]))
end

#
# jac = hadamard_jacobian((Ax,Ay),:sym,dF,Q)

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
using DataFrames, MCMCChains, CSV

"""

# read_summary

Read summary output file created by stansummary.

### Method
```julia
read_summary(m)
```

### Required arguments
```julia
* `m`    : A Stan model object, e.g. SampleModel
```

"""
function read_summary(m::T) where {T <: CmdStanModels}

  fname = "$(m.output_base)_summary.csv"
  !isdir(fname) && stan_summary(m)

  df = CSV.read(fname, delim=",", comment="#")

  cnames = lowercase.(convert.(String, String.(names(df))))
  cnames[1] = "parameters"
  cnames[4] = "std"
  cnames[8] = "ess"
  names!(df, Symbol.(cnames))
  df[!, :parameters] = Symbol.(df[!, :parameters])

  ChainDataFrame("CmdStan Summary", df)

end   # end of read_samples


text:
using ApproxFun, Plots
import ApproxFun: factor, component

L = 1.0
T = 0.1
d = Interval(0.,T) × Interval(0.,L)
Dx = Derivative(d,[0,1])
Dt = Derivative(d,[1,0])

α = 1.0

B = [ldirichlet() ⊗ I;
     I ⊗ ldirichlet();
     I ⊗ rdirichlet();
     (I ⊗ rneumann()) - α*(I ⊗ lneumann());]

q0 = Fun(x->x^2*(L-x)^2 , factor(d,2))
f0 = zeros(component(∂(d),1))
g0 = zeros(component(∂(d),3))

u = \([B;Dt + Dx^3] , [q0;f0;g0;0;0] ; tolerance=1E-5)
surface(u)

max_stars_repo_path: PDEs/Linear KdV.jl

max_issues_repo_path: PDEs/Linear KdV.jl

max_forks_repo_path: PDEs/Linear KdV.jl


text:
function timetag()
    n = string(now())
    reduce(replace, ["-" => "", ":" => "", "T" => "", "." => ""]; init=n)
end

function save_model(res, modelpath, val_labels)
    id = timetag()
    f = joinpath(modelpath, "$(id)_layers-$(length(res[1].e)).bson")
    save(f,
        Dict(:model => cpu(res[1]), :hist => res[2],
        :rdata => res[3], :zs => res[4], :val_labels => val_labels))
    f
end

function load_model(id, modelpath)
    f = filter(x->occursin("$id", x), readdir(modelpath))[1]
    modeldata = load(joinpath(modelpath,f))
    return gpu(modeldata[:model]), modeldata[:hist], modeldata[:rdata], modeldata[:zs], modeldata[:val_labels]
end
function load_model_gpu(id, modelpath)
    res = load_model(id, modelpath)
    return gpu(res[1]), res[2], res[3], res[4], res[5]
end

sample_tensor(x::AbstractArray{T,4}, n::Int; kwargs...) where T = x[:,:,:,sample(1:size(x,4), n; kwargs...)]
sample_tensor(x::AbstractArray{T,2}, n::Int; kwargs...) where T = x[:,sample(1:size(x,2), n; kwargs...)]

function batched_loss(lf, x, batchsize)
    @suppress begin # suppress the batchsize warning
       return Flux.mean(map(y->cpu(lf(y)), Flux.Data.DataLoader(x, batchsize=batchsize)))
    end
end

"""
    restart_weights!(model, init = Flux.glorot_uniform)

Restarts the weights of the model using given initialization method.
"""
function restart_weights!(model, init = Flux.glorot_uniform)
    ps = params(model)
    for p in ps
        p .= init(size(p)...)
    end
    model
end

"""
    normalization_coefficients(x::Matrix)

Row-wise normalization coefficients.
"""
normalization_coefficients(x::AbstractArray{T,2}) where T = maximum(x, dims=2), minimum(x, dims=2)

"""
    normalize(x::Matrix[, maxs, mins])

Row-wise normalization of x. If maxima/minima are not supported, they are returned as well.
"""
normalize(x::AbstractArray{T,2}, maxs, mins) where T = (x .- mins) ./ (maxs - mins)
function normalize(x::AbstractArray{T,2}) where T
    maxs, mins = normalization_coeffs(x)
    return normalize(x, maxs, mins), maxs, mins
end
sqnorm(x) = sum(abs2, x)

max_stars_repo_name: vitskvara/HierarchicalAD

max_issues_repo_name: vitskvara/HierarchicalAD

max_forks_repo_name: vitskvara/HierarchicalAD